This paper is an extension of our previous work on simulation of complex phase front evolution in the diffusion-dominated situation. The Navier-Stokes equations are solved using a finite-volume method based on a second-order accurate central-difference scheme in conjunction with a two-step fractional-step procedure. The key aspects that need to be considered in developing such a solver are imposition of boundary conditions on the immersed boundaries and accurate discretization of the governing equation in cells that are cut by these boundaries. A new interpolation procedure is presented which allows systematic development of a spatial discretization scheme that preserves the second-order spatial accuracy of the underlying solver. The presence of immersed boundaries alters the conditioning of the linear operators and this can slow down the iterative solution of these equations. The convergence is accelerated by using a preconditioned conjugate gradient method where the preconditioner takes advantage of the structured nature of the underlying mesh. The accuracy and fidelity of the solver is validated and the ability of the solver to simulate flows with very complicated immersed boundaries is demonstrated. The method will be useful in studying the effects of fluid flow on the evolution of complex solid-liquid phase boundaries.