Boiling flows are central to many industrial and natural processes. The high heat transfer rate and the ability of fluids to store large amount of energy in the form of latent heat make boiling particularly important in large scale energy generation and thermal energy storage. Because of the large volume change and the high temperatures involved, however, consequences of design or operational errors can be catastrophic, and accurate predictions of heat transfer and fluid flow are fundamental for safe operations. Although boiling has been studied extensively, a more basic understanding of the subject has been hindered by the small spatial and temporal scales which prevent accurate experimental measurements. The nucleate pool boiling process is known to be a very efficient mode of heat transfer. It is desirable to operate many engineering applications in this mode since high heat transfer rates and convection coefficients are associated with small values of the excess temperature. Despite its importance, nucleate pool boiling has not been fully numerically simulated until very recently because of the complexity of dealing with the phase change problem in addition to the existence of the interface between two phases. Nucleate pool boiling is typically characterized by cyclic growth and departure of vapor bubbles from a heated wall. In this study, direct numerical simulation of saturated nucleate pool boiling is performed using a front tracking method. Here, the simulated phenomenon occurred in a two dimensional rectangular domain under the constant heated wall temperature and laminar flow. The complete cycles of nucleate pool boiling including nucleation seed growth, rising and departing from the heated wall is captured using Euliar-Lagrangian point of view associated with hybrid front tracking method. Effect of different non-dimensional parameters like, Gr, Pr and Ja numbers as well as liquid/vapor properties on space averaged heat flux, in a given nucleation site, are studied. Effect of initial nucleation sites is studied as well. The obtained results shows that, an increase in Gr number leads to rize in the magnitude of vertical velocity component, and therefore a change in the shape of the detached and rising bubble, and decreasing in bubble detachment frequency. Analysis of heat flux reveals that under the considered condition, by increasing Gr and Pr numbers, the amount of space averaged heat flux is reduced. And there is less bubble departure on heated wall where there is lower bubble departure frequency. Also, results of space averaged heat flux with respect to Ja number shows that by increasing Ja number, the amount of space averaged heatflux is increased. By increasing Ja number, amount of bubble formation and frequency are increased. The obtained results are compared with that predicted by published theoretical correlation showing acceptable agreement in different simulated cases. Also for multiple nucleation sites, simulations were performed and their subsequent results are presented.The amont of space averaged heat flux of these simulations were calculated for nucleate pool boiling and show that by increasing nucleation sites, the amount of space averaged heat flux is increased. A grid refinement test for nucleate pool boiling on a horizontal surface for different grids shows the convergance of the results. Keywords: Heat transfer, Nucleate boiling, Front tracking, Pool boiling, Heat flux, Multi phase, Numerical simulation