A mathematical model for turbulent gas-liquid flows is presented.
It is shown that bubble-induced buoyancy resembles natural convection
and can be readily incorporated into an incompressible flow solver
by using  an analog of the Boussinesq approximation. Extra transport
equations  are introduced to describe the evolution of the gas holdup
and compute the turbulent eddy viscosity in the framework of a generalized
k-epsilon model. A robust solution strategy based on nested iterations
is proposed for the numerical treatment of the resulting PDE system.
The incompressible Navier-Stokes equations are solved by a discrete
projection scheme from the family of Pressure Schur Complement methods.
High-resolution finite element schemes designed on the basis of
algebraic flux correction are employed for the discretization of
convective terms. Numerical results for prototypical gas-liquid
reactor configurations illustrate the performance of the developed
simulation tools.