The present numerical study deals with the stabilized finite element solution of the steady natural convection flow in terms of primitive variables. The domain is discretized into a set of regular linear triangular elements for all unknowns containing a single subgrid node, which has a crucial effect in the stabilization with a cheap computational cost. The proposed stabilized finite element method enables obtainment of a stable solution and avoids oscillations, especially in the pressure. The algorithm is presented for different benchmark problems considered in different geometries and the solutions are obtained for high values of Rayleigh number (Ra) (103Ra106). Furthermore, the solutions show the well-known characteristics of the natural convection flow and agree well with the results in the literature.