In the present study a more complete numerical solution approach using parallel computing technology is provided for the three-dimensional modeling of mold insert polymer injection molding process by considering the effects of phase-change and compressibility for non-Newtonian fluid flow conditions. A volume of fluid (VOF) method coupled with a finite volume approach is used to simulate the mold-filling stage of the injection molding process. The variations in viscosity and density in the polymer melt flow are successfully resolved in the present VOF method to more accurately represent the rheological behavior of the polymer melt flow during the mold filling. A comprehensive high-resolution differencing scheme (compressive interface capturing scheme for arbitrary meshes or CICSAM) is successfully utilized to capture moving interfaces and the pressure-implicit with splitting operators pressure-velocity coupling algorithm is employed to enable a higher degree of approximate relation between corrections for pressure and velocity. The capabilities of the proposed numerical methodology in modeling real molding flow conditions are verified through quantitative and qualitative comparisons with other simulation programs and the data obtained from the experimental study conducted. The present numerical results are also compared with each other for a polypropylene female threaded adaptor pipe fitting model with a metallic insert for varying governing process conditions/parameters to assess the modeling constraints and enhancements of the present numerical procedure and the effects of these conditions to optimize the polymer melt flow for mold insert polymer injection molding process. The numerical results suggest that the present numerical solution approach can be used with a confidence for further studies of optimization of design of mold insert polymer injection molding processes.