This study presents the development of a front tracking method for melting of phase change materials (PCMs) inside horizontal shell and tube heat exchangers. Two numerical models, i.e. pure conduction (PC) model as well as combined conduction and natural convection (CCNC) model, are used to develop the method. Governing equations are numerically solved by ANSYS Fluent v17.2. The PC model benefits from simplicity but its prediction is far from reality, whereas CCNC model's prediction is more realistic but its modeling is complicated. Generally, during the melting process, the upper half of the system is affected by the upward buoyancy-driven melted PCM motion. To consider this phenomenon, the front tracking method assumes that the upper and lower halves of the system have two separate melting fronts. Therefore, it is assumed that the natural convection contributes only to the upper half until the upper half liquid fraction value reaches unity. Meanwhile, the lower half melting front is assumed to be the same as that of the PC model. Once the upper half is totally melted, the method attributes the rest of the natural convection to the lower half of the system. Using three different PCMs and three different geometries, two correlations have been developed for each half based on two dimensionless numbers; i.e. the shell-to-tube radius ratio and PC model liquid fraction. The method is then verified using another PCM, which has not been included during the correlation development stage to guarantee the methods validity. These correlations provide results within +/- 15% discrepancy range.