In this article, two-phase slug flow is simulated numerically in a horizontal duct with rectangular cross-section using Volume Of Fluid (VOF) method. Conservation equations of mass, momentum and advection equation are solved in open source OpenFOAM code accompanying k-ω SST turbulence equations. Simulation is conducted based on the experimental results in the duct with rectangular cross-section. The results shows, due to Kelvin-Helmholtz (K-H) instability criteria slug initiation forms in the air-water interface during three dimensional turbulence modeling. Water level was increased slightly at interface in both numerical simulation and experiment. This level increase satisfies the K-H instability to generate a slug at interface. During slug initiation, the pressure behind slug is increased significantly. Big pressure gradient at the beginning of the slug in compare to the end of it causes the slug length to be increased as propagate along the duct. The numerical simulation of present research is capable of predicting the slug length accurately in accordance with experiment; however, the slug position with 22% inaccuracy was obtained. Comparison of the results with the numerical and experimental results of other researchers confirms higher accuracy of flow prediction in the present work.