An understanding of the hydrodynamics of multiphase processes is essential for their design and operation. Multiphase computational fluid dynamics (CFD) simulations enable researchers to gain insight which is inaccessible experimentally. The model frequently used to simulate these processes is the two-fluid (Euler-Euler) model where fluids are treated as inter-penetrating continua. It is formulated for the multiphase flow regime where one phase is dispersed within another and enables simulation on experimentally relevant scales. Phase fractions are used to describe the composition of the mixture and are bounded quantities. Consequently, numerical solution methods used in simulations must preserve boundedness for accuracy and physical fidelity. In this work, a numerical method for the two-fluid model is developed in which phase fraction constraints are imposed through the use of an nonlinear variational inequality solver which implicitly imposes inequality constraints. The numerical method is verified and compared to an established explicit numerical method.