In　this　paper,　we　propose　a　subgrid　finite　element　method　for　the　two-dimensional　(2D)　stationary　incompressible　Naver-Stokes　equation　(NSE)　based　on　high　order　finite　element　polynomial　interpolations.　This　method　yields　a　subgrid　eddy　viscosity　which　does　not　act　on　the　large　scale　flow　structures.　The　proposed　eddy　viscous　term　consists　of　the　fluid　flow　fluctuation　stress.　The　fluctuation　stress　can　be　calculated　by　means　of　simple　reduced-order　polynomial　projections.　Assuming　some　regular　results　of　NSE,　we　give　a　complete　error　analysis.　Finally,　in　the　part　of　numerical　tests,　the　numerical　computations　show　that　the　numerical　results　agree　with　some　benchmark　solutions　and　theoretical　analysis　very　well.