An ALE framework for multiphase flows

Abstract

An Arbitrary Lagrangian Eulerian (ALE) framework which combines the advantages of both Lagrangian and Eulerian methods is developed to solve incompressible multiphase flow problems. The div-stable side centered unstructured finite volume formulation is used for the discretization of the incompressible isothermal Navier-Stokes equations along with the isothermal constitutive equations for Oldroyd-B and FENE-CR fluids. In this approach, the velocity vector components are defined at the mid-point of each cell face, while the pressure term and extra stress tensor are defined at element centroids. The present arrangement of the primitive variables leads to exact total mass conservation at machine precision due to the present stable numerical discretization with no ad-hoc modifications. In addition, a special attention is given to satisfy global discrete geometric conservation law (DGCL) at discrete level for the application of the interface kinematic boundary condition in order to conserve the total mass for each species for multiphase flow problems. Furthermore, the pressure field and extra stress field are treated to be discontinuous across the interface with the discontinuous treatment of density and viscosity and jump conditions are satisfied. Surface tension force is treated as a tangent force and discretized in a semi-implicit form. Two different approaches for the computation of unit normal vector have been implemented: the least squares biquadratic surface fitting (LSBSF) and the mean weighted by sine and edge length reciprocals (MWSELR). The combination of MWSELR method and discontinuous treatment of density and viscosity reduced the parasitic currents to the machine precision. The resulting large system of algebraic equations is solved in a fully coupled manner in order to improve the time step restrictions. As a preconditioner, an approximate matrix factorization similar to that of the projection method is employed and the parallel algebraic multigrid solver BoomerAMG provided by the HYPRE library, which is accessed through the PETSc library, has been utilized for the scaled discrete Laplacian of pressure and the diagonal blocks of mesh deformation equations. The present calculations verify that the mass of the bubble can be conserved at machine precision independent of spatial and temporal resolutions. The accuracy of the proposed method is initially validated on the static bubble problem, since the surface tension force is highly sensitive to the accurate evaluation of the unit normal vector and the inaccuracies significantly contribute to unphysical velocities, called parasitic currents. The calculations indicate that the parasitic currents can be reduced to machine precision for the MWSELR method. The MWSELR approach, as far as our knowledge goes, has not been used for the evaluation of normal vectors in multiphase flows. In the second benchmark case, the proposed approach is applied to the single bubble rising in a viscous quiescent liquid for both low and high density ratios. The calculations produce accurate predictions of the bubble shape, center of mass, rise velocity, etc. Furthermore, the mass of each species is conserved at machine precision and discontinuous pressure field is obtained in order to avoid errors due to the incompressibility restriction in the vicinity of liquid-liquid interfaces at large density and viscosity ratios. The third benchmark case is rising of a Taylor bubble in 2D and in 3D. Taylor bubbles are large bullet shaped bubbles whose cross-section almost fill the cross-sectional area of the channel. Therefore, this benchmark case is numerically harder than the previous cases. It is seen, 3D bubble rises faster due to the smaller blockage effect (i.e. cross section of the bubble/cross section of the tube) of the bubble in three dimension with respect to the 2D bubble. In addition, drag force of the bubble decreases due to the three-dimensional relieving effect. The results are compared with the results available in the literature and it is shown that the obtained bubble shape and velocity field in the vicinity of the Taylor bubble are similar to that of the literature. In the fourth test case, rise of a single bubble in a quiescent, viscoelastic fluid due to buoyancy is simulated in 2D and the viscoelastic fluid is modeled as Oldroyd-B. By changing the size of the bubble, domain, placing the bubble to the different locations and changing the fluid properties, many simulations are performed and the change in bubble shape, rise velocity, circularity and sphericity are inspected. It is seen that the existence of the wall highly effects the outcome. In addition, the cusp at the trailing edge of the bubble and negative wake behind the bubble are observed in some cases. Therefore, it is shown that a viscoelastic fluid model that exhibits shear thinning is not essential for negative wake to arise. This result contradicts with the some published papers in the literature but is also supported by the others. The final benchmark case is similar to the previous one but this time viscoelastic fluid is modeled as FENE-CR and the problem is in 3D. Besides, subsequent simulations are performed for Newtonian bubble and Newtonian continuum, Newtonian bubble and viscoelastic continuum, and viscoelastic bubble and Newtonian continuum. It is observed that the bubble has a slight cusp at the trailing edge for Newtonian bubble and viscoelastic continuum. On the other hand, the bubble has a dimple at the trailing edge for the viscoelastic bubble and Newtonian continuum. In addition, it is shown that the results are in a good agreement with the result available in the literature. Finally, the methods used to develop and test the present multiphase solver for both Newtonian and viscoelastic fluids are summarized. Advantages and the drawbacks of the present solver are addressed with possible future applications.