The presented work is dedicated to the mathematical and numerical modeling of unsteady single- and two-phase flows using finite volume and penalty methods. Two classes of Navier-Stokes solvers are considered. Their accuracy and robustness are compared to identify their respective strengths and weaknesses. Exact (also referred to as monolythic) solvers such as the Augmented Lagrangian and the Fully Coupled methods address the saddle-point structure on the pressure-velocity couple of the discretized system by means of a penalization term or even directly, whereas approximate (segregated) solvers such as the Standard Projection method rely on operator splitting to break the problem down to decoupled systems. The objective is to compare all approaches in the context of two-phase flows at high viscosity and density ratios. To characterize the interface location, a volume-of-fluid (VOF) approach is used based on a Piecewise Linear Interface Construction (PLIC). Various 2D simulations are performed on single- and two-phase flows to characterize the behavior and performances of the various solvers.