A numerical study on natural convection induced by free surface heat flux and cold left and hot right walls in glass melting furnaces has been performed. A function of heat flux derived from the combustion environments of actual glass melting furnace is applied to thermal boundary condition at free surface. Fundamentally there exist two flow cells in cavity (left counterclockwise one and right clockwise one). The effects of heat flux and Rayleigh number are investigated through two-dimensional steady-state assumption. The convection strength of two flow cell located in left region continuously increases. In the mean time the strength of flow cell in right region increases and then decreases. Critical Rayleigh number in which two flow cells take place above and below show linear dependence on the free surface heat flux. To maintain the traditional flow pattern (left and right flow cells) in glass melting furnace, Rayleigh number is recommended to be below 10$^{5}$ .