Pulsating Heat Pipes (PHPs) are passive two-phase thermal transport devices, offering compact and high-performance solutions for electronics cooling. Yet, numerical modeling of PHPs using Volume of Fluid (VOF)-based Computational Fluid Dynamics (CFD) methods remains challenging due to limited clarity in prior studies regarding key modeling assumptions, particularly the selection of mass relaxation coefficients in the Lee phase-change model and a lack of systematic analysis on numerical stability. This study introduces a structured computational framework for conducting time-efficient and numerically stable VOF simulations of PHPs. Although previous studies have independently employed both 2D and 3D simulations, a direct comparison between these approaches is lacking, and the influence of phase-change coefficients on mass balance within the system has not been reported. This study aims to fill that gap by directly comparing 2D and 3D simulations conducted on the same geometry. While 3D simulations with relatively coarser mesh tend to overpredict thermal resistance, the 2D simulations yield stable results with thermal resistance values in good agreement with experimental data owing to the sufficient mesh resolution near the wall, which enabled accurately capturing the thin liquid film around vapor bubbles. By monitoring mass variation and systematically analyzing the influence of parameters such as time step size and phase-change coefficients, the framework identifies conditions under which stable and physically consistent results are achieved. The proposed framework provides a practical and adaptable approach for modeling PHPs, offering clear guidance for achieving reliable, convergent, and physically meaningful simulations in two-phase flow systems.