  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.
