The evolution of numerical methods and computational facilities allow re- searchers to explore complex physical phenomenons such as multiphase flows. The specific regime of incompressible, turbulent, bubbly two-phase flow (where a car- rier fluid is infused with bubbles or particles) is also receiving increased attention due to it’s appearance in major industrial processes. The main challenges arise from coupling individual aspects of the physics into a unified model and to provide a robust numerical framework. The presented work aimed at to achieve the second part by employing the most frequently used dispersed two-phase flow model and another incompressible, turbulent single phase solver as a base flow provider for coupled Lagrangian or surface tracking tools. Among the numerical techniques, the finite element method is a powerful can- didate when the need arises for multiphysics simulations (for example coupling with an electrochemical module) where the counterpart has a node based ap- proach. Stabilization schemes such as PSPG/SUPG/BULK provide remedies for the pressure decoupling and the inherent instability of the central discretization when applied for convective flow problems. As an alternative to unsteady solvers based upon an explicit or a fully im- plicit nonlinear treatment of the convective terms, a semi-implicit scheme results in a method of second order accurate in both space and time, has absolute linear stability and requires only a single or two linear system solution per time step. The application of the skew symmetric approach to the convective term further stabilizes the solution procedure and in some cases it even prevents divergence. The Eulerian-Eulerian two-phase flow model poses various issues to be over- come. The major difficulty is the density ratio between the phases; for an ordinary engineering problem it is in the order of thousands or more. The seemingly minus- cule differences in the formulation of the stabilizations can cause very different end results and require careful analysis. Volume fraction boundedness is of concern as well, but it is treatable by solving for its logarithm. Since the equations allow jumps (even separation of the phases) in the volume fraction field, discontinuity capturing techniques are also needed. Besides the standard ’spatial’ stabilization temporal smoothing is also necessary, otherwise the limitation in time step size becomes too stringent. Designing a flow solver is one side of the adventure, but verification is equally important. Comparison against analytical solution (such as the single and two- phase Taylor-Green testcase) provides insight and confirmation about the mathe- matical and physical properties. Meanwhile comparing with real life experiments prove the industrialization and usability of a code, dealing with low quality meshes and effective utilization of computer clusters. / Doctorat en Sciences de l'ingénieur et technologie / info:eu-repo/semantics/nonPublished