A FULLY ADAPTIVE, CONSERVATIVE FRONT TRACKING METHOD FOR THE SIMULATION OF INCOMPRESSIBLE MULTIPHASE FLOWS

This paper presents a fully adaptive formulation of the Front Tracking method for the simulation of incompressible, multiphase, bubbly flows, based on the Tryggvason formulation. The Navier-Stokes equations are discretized using a finite difference scheme, and domain discretization is carried out with Berger & Collela’s structured adaptive mesh refinement (SAMR) algorithm. Time discretization is based on SBDF scheme, with adaptive time stepping. The lagrangian interface is represented using the GTS library, which provides a volumeand shapepreserving remeshing algorithm, therefore minimizing the volume change due to non-conservative interpolation of the eulerian velocity field. Nevertheless, a simple volume recovery algorithm is also provided, along with a subgrid undulation removal algorithm based on the TSUR-3D algorithm [6]. Rising bubble flows were simulated under several regimes, showing small errors when comparing to experimental results.


INTRODUCTION
Front-tracking methods are based on the one-fluid formulation for the simulation of multiphase flows and make use of markers, usually a set of connected points forming an unstructured mesh, to track the interface. Since the interface is explicitly represented by its coordinates, the set of markers is said to use a Lagrangian framework, and the mesh which characterizes the interface is termed Lagrangian mesh. While this interface is based on a Lagrangian framework, an Eulerian framework is used to solve the transport equations. This method is accurate but also complex to be implemented because dynamic remeshing of the Lagrangian interface mesh is required and mapping of the Lagrangian data onto the Eulerian mesh has to be carried out. Difficulties may arise when multiple interfaces interact with each other as in coalescence and breakup, both of which requiring a proper sub-grid model. This feature may be computationally expensive, but can be used to prevent artificial merge or breakup when the number of bubbles in the flow is large.
Although storing the marker points is enough to fully determine the interface position, storing the interface as a mesh is useful for computing geometrical properties such normal, area and the interfacial tension force, not to mention checking against collisions or merge/breakup phenomena.
There are two possible strategic forms of communication. Glim et al. [8] discretized the Navier-Stokes equations using a finite difference method which modiffies the stencil in the vicinity of the mesh overlap, so that the vertices of the two meshes match each other. However, for multiphase flows the most widespread methodology is based on the Immersed Boundary Method (IBM) of Peskin [13,14], which represents the interface by imposing a force field that is spread on the Eulerian mesh. The Navier-Stokes equations are solved based on one-fluid models. Among various methodologies developed after Peskin, one of the most wide spread in the context of multiphase flows is the work of Uverdi and Tryggvason [20] and Trygvason et al. [19]. Similarly to IBM, the interface between the phases is represented by the interfacial tension force field.
In the present work the Front-Tracking method is combined with an AMR (Adaptative Mesh Refinement) framework. The concept of adaptive mesh refinement is based on the works of Berger et all [2] and Jameson [9] and consist in creating a hierarchy of Cartesian meshes with different levels of refinement which cover the entire domain, concentrating the meshes on the regions which require special attention. All meshes in a given level of refinement must have the same characteristic length.

MATHEMATICAL MODEL
The mathematical model used in the present work is composed of the Navier-Stokes and mass equations combined with a front tracking model with which problems of two phase flows are modeled. The model is presented in three parts, as follows.

Mathematical model for the fluid
This model, for a Newtonian fluid and for incompressible flows, is composed by the classical mass and Navier-Stokes Equations.
where u(x, t) is the uid velocity, ρ is the uid density, and µ is the viscosity, p is the pressure and g is the gravity acceleration. The source term models the interfacial force that is different from zero only over an interface and equals zero everywhere. The fluid density and viscosity are calculated with the aid of an indicator function: Usually the indicator function is found by solving a Poisson equation for the surface normal vector [20]. In this work, however, a geometric approach, named CPT (Closest Point Transform) is adopted [12]. Instead of solving an additional elliptic PDE, CPT finds the shortest distance between a given point in the Eulerian domain and the Lagrangian interface. Unlike the Tryggvason approach, in this case the function is calculated only in the neigbourhood of the interface [5].

Mathematical model for the interface
The interface is modeled by an equation for the interface force field and the equations for its Lagrangian transport: where σ is the interface force coefficient, κ is the curvature, n is the interface normal, X(t) is the spatial Lagrangian coordinate, dΩ is an area element and dΓ is the boundary of dΩ.

Mathematical model for the fluid-interface communication
Information exchange between fluid and interface is given by interpolation and distribution equations where X and x are, respectively, positions in the Lagrangian (interface) and the Eulerian (fluid) domains.

Numerical method
The time discretization of the differential model is performed by the general model below: ∇ · u n+1 = 0 where λ = L ∞ (µ) and h depends on the diffusive-, advective-and forcing terms, according to the following equation: The value of α, β and θ are given by: where ω = ∆t n+1 ∆t n is the rate between two consecutive time steps In the present work the above SBDF temporal time discretization, with variable size of time step, was chosen in all cases simulated. To handle the pressure-velocity coupling between Eq. 1 and Eq. 2, the following fractional step method is proposed: The pressure correction q is obtained using the following equation:

Spatial discretization
Space is discretized using a structured adaptive mesh refinement, which is based on the Immersed Boundary method introduced by Roma et al. [16], and in the hierarchical grid structure proposed by Berger and Colella [3]. In this scheme, regions of the flow being of special interest are covered by block-structured grids, defined as a hierarchical sequence of nested, progressively finer levels (composite grids). Each level is formed by a set of disjoint rectangular grids and the ratio between two successive refinement levels is constant and equal to two. Ghost cells are employed around each grid, for all the levels, and underneath finer grid patches to formally prevent the finite difference operators from being redefined at grid borders and at interior regions which are covered by finer levels.
Values defined in these cells are obtained from interpolation schemes, usually with second or third order accuracy, and not from solving the equations of the problem. A staggered composite grid is used, i.e., pressure and others scalar variables are placed at the centers of the computational cells, while vector variables are located at the respective cell faces. The discretization of the Laplacian, gradient and divergence differential operators are performed by standard, cell-centered second order finite difference stecils. Although a variable can be defined or initialized in any level, in the current work the interface must be completely covered by the finest level, ensuring that the most important physical phenomena are being captured.
The initial mesh for the Lagrangian interface is generated using GMSH [7]. Volume conservation is ensured in the surface remeshing during the simulations by applying the memoryless simplification technique by Lindstrom and Turk [11], as provided by the GTS Library [15]. Nevertheless, the interpolation of the velocity field from the eulerian framework onto the Lagrangian interface usually is non-conservative. Although the volume change in a single time step is nearly negligible, it is a cumulative process that may lead to considerable volume change in long term simulations. To overcome this issue, a simple volume recovery procedure applied at the end of the advection step will be described next.

Volume Recovery Algorithm
The mass conservation issue in Front Tracking arises after interpolating the Eulerian velocity field to the Lagrangian interface, because the interpolation functions are not conservative. This problem can be tackled using two different approaches: (i) implementing a conservative interpolation scheme or (ii) applying a correction step after interpolating the velocity field in order to correct the surface volume. This second approach will be followed here, since it is much simpler than (i).
Average based volume corrections are not new in Front Tracking Methods [17,18] and usually are carried out appart from the remeshing procedure, such that no topological change is introduced in the mesh. If the volume change is small, it is safe to assume that all variation occurs in the normal direction. Since the interface is discretized by triangular elements, the volume change can be seen as the sum of small prisms and formed by the mesh elements and a small constant height η: where A i is the area of the ith-element of the mesh, η is the average length of the normal and N e is the total number of elements in the mesh. Since h is constant, eq. 18 can be written as where A is total surface area. Also, the volume change may also be computed as the difference between the initial and current interface volume, i.e.: 20) and the length of the normal will be just That is, the volume correction can be achieved by just expanding or shrinking the surface in the normal direction by a factor equal to η. The correction is then given by eq. 22, where x old and x new are the original and new coordinates of the mesh vertices.
Since this movement will be applied to the vertices, an accurate pseudo-normal computation must be provided. A good approximation is given by an weighted average of the normals of the elements sharing a given vertex [10] according to Eq. 23: Here, N t is the number of triangles which share the vertex, n i is the normal of element i, e i1 and e i2 are the edges of the ith-element which are employed to compute this element's normal, and α i is the angle between them.

RESULTS
Firstly a verification of the volume-conserving abilities of the front tracking algorithm will be made, based on an analytical shear flow test. Then, results showing the simulation of rising bubbles will be presented.

Analytical Shear Flow
The test consists in subjecting a surface, usually spherical, to a solenoidal, periodic, analytic velocity field and assessing its volume change as the flow evolves. Assuming a velocity field with a period T , the surface is stretched until t = T /2 and then returns to its original position and shape. The velocity field may be calculated either directly on the Lagrangian vertices or interpolated from an auxiliary grid, allowing to assess the influence of the interpolation functions on the volume change during the simulation.
Regarding the interface advection, the two most common time integration schemes are first order Euler and second order trapezoidal rule, respectively Eqs. 24a and 24b. Note that, as long as there is a time history of the Eulerian velocity field, it is possible to implement second order methods for the advection of the Lagrangian interface.
Surface stretching requires remeshing, which also may also interfere on the conservation process. Regardless its volume conservation properties, any remeshing algorithm destroys the Lagrangian velocity history, requiring a new interpolation of the velocity field. Since interpolation schemes are usually non-conservative, a conservative remeshing algorithm could still lead to a non-conserving outcome owing to the need for velocity interpolation.
In order to assess the influence of such variables on the flow outcome, a series of tests was performed, in which three main parameters were analysed: (i) the velocity field interpolation, (ii) the surface remeshing and (iii) the volume recovery algorithm. These cases are summarized in The velocity profile is given by equation 25 [17] and the interpolation of the velocity field was performed using equation 26 [21].
u(x, y, z, t) = cos(πt/T )sin 2 (πx)(sin(2πz) − sin(2πy)) v(x, y, z, t) = cos(πt/T )sin 2 (πy)(sin(2πx) − sin(2πz)) w(x, y, z, t) = cos(πt/T )sin 2 (πz)(sin(2πy) − sin(2πx)) (25) δ(r) = 1 4 (1 + cos( π Figure 1 shows snapshots of the sphere shape during a typical simulation at t = 0, t = T /4, t = T /2, t = 3T /4 and t = T . The flow stretches the surface until T = T /2 and then symmetrically returns it to its initial position, so that the shape and position are the same not only t = 0 and t = T , but also at t = T /4 and t = 3T /4. The period is T = 4s, the same which was used in the simulations presented here. Figure 1. Geometry deformation when subjected to the analytical velocity field given by Eq. 25. Snapshots taken at t = 0, t = T 4 , t = T 2 , t = 3T 4 e t = T Figure 2 shows the time evolution of the relative volume error (ε vol = (V f /V i − 1) × 100) for all cases simulated. Graphs on the left column show the cases in which volume recovery was not applied, while the right column shows the results achieved with volume recovery. Notice that, in this case, the error was kept within the specified threshold, that is ±0.01%. Regarding the non-conservative cases, both interpolation process and interface remeshing have shown to play little influence on the final volume error when the time integration is carried out with Euler scheme. Although velocity interpolation changes the time history of volume error, the final value achieved is similar in both cases. The trapezoidal rule, although being more sensitive to both surface remesh and velocity interpolation, leads to considerably smaller errors.

Rising Bubbles
This section presents the simulation of rising bubbles in various flow regimes. Firstly, a case study is performed in order to assess the influence of some geometry-correction tools on the outcome of the flow. Then a set of low Reynolds cases is simulated, in order to assess the influnce of the Morton number on the prediction of the terminal Re number. All cases in this set have the same Eo number. Finally, a wobbling case is simulated. In all cases simulated in this paper, the initial distance from the center of the bubble to the domain boundaries is four times the bubble diameter. Also, the maximum edge length allowed for the Lagrangian mesh is half the length of the Eulerian grid spacing at its finest level.

Geometry-Correction Tools
The physical parameters defining the rising bubble flow regime in this section are: Eo = 40, M = 0.056, Re = 20.6, ρ C = 1000kg/m 3 , µ C = 0.273556P a · s, λ ρ = 100, λ µ = 100. The bubble diameter is φ = 0.02m. Eo is the Eotvos number, M is the Morton number, ρ C and µ C are the density and dynamic viscosity of the continuous phase and λ ρ and λ µ are the density ratio and viscosity ratio between the continous and dispersed phase, respectively.
Besides the well-known non-conservative behaviour of the velocity interpolation process, another common problem in front tracking methods is the development of small un- dulations or wrinkles that occurs at the rear of rising bubbles. The TSUR-3D algorithm [6] eliminates such undulations by requiring that the normal of a given edge in the mesh be parallel to the resultant of the approximate normals of the vertices of this edge.
Four cases were run in order to assess the role of the geometric correction tools on the flow simulation. In case A1 the standard Front Tracking Method was used, without any further intervention. In case A2 a volume correction step was applied at the end of the interface advection, while in case A3 only TSUR-3D was used. Case A4 applies both, volume recovery and TSUR-3D algorithms. Table 2 shows the terminal Re achieved in each simulation, as well as the relative error.

Case Description
Re Error (%)
The effect of the geometric tools on the Re number and terminal velocity is shown in figure 3. Notice that the volume recovery algorithm has little influence at the beginning of the flow. In fact, only after reaching the terminal regime a small decrease in the Re number can be seen. Surface wrinkle removal, on the other hand, affects the early stages of the flow by smoothing the transient region and decreasing the Re maximum at the overshoot region. The final position of the bubbles, depicted on the right, shows that the algorithm affects the rise velocity, not the bubble volume.  Figure 4 shows details of the bubble shape at the end of the flow for cases A1 and A3. Figure 5 show that the presence of the wrinkles lead to the formation of artificial vortices on the bubble rear, influencing the flow inside and outside the bubble.

Low Reynolds flows
After assessing the influence of the geometry correction tools on the bubble behaviour, a set of low Re flows was simulated using the volume recovery and TSUR-3D algorithms. For all cases simulated, Eo=116 while M ranged between 1.31 and 848. Figure 6 shows the time history of the Reynolds number for all cases on the left and the error on the terminal Re on the right. Table 3 shows the results of the simulations and compares the terminal bubble shape and Re number to the experimental work of Bhaga [4]. Despite the difference in the density ratio (1370 in the experiment vs 100 here), the errors were small even for high Morton numbers, which usually demands more effort on the solver. Notice also that in all cases the predicted shapes are in good agreement with the experiments.

CONCLUSIONS
The numerical method to solve the Navier-Stokes equations and the front-tracking method for the two phase flow interfaces treatment was implemented in an AMR framework. The methodology was verified and validated under comparison with experimental results. The preliminaries conclusions are that the methodology represents well the main physical characteristics of this very complex problem. The next steps for the present project is to improve this methodology in order to capture coalescence and fragmentation of bubbles in turbulent two phase flows as well as to make numerical experiments to obtain statistical parameters as drag, lift, added mass and Basset force.