On the Solution of the Problem of a Drop Falling Against a Plane by Using a Level Set – Moving Mesh – Immersed Boundary Method

The interaction of bubbles and drops with solids constitutes a broad research topic with numerous applications. Its study is important in fields as sprays, mineral flotation or cooling of nuclear reactors. For instance, in the study of sprays, the behavior of the atomized drops when approaching surfaces could condition the applicability of a specific technique [1]. In biomedicine, the aforementioned problem has a crucial influence in echocardiography, where bubbles could be used as the contrast medium during the test [2]. These and many other topics make the problem of the approaching of a drop against a wall an essential problem to be deeply understood.

Although the problem of a drop/bubble moving in a non-constricted domain has been thoroughly studied (see for instance [3]), the movement of a bubble/drop through complex geometries is still a relatively unexplored field. Some progresses have been done in this latter case, by experimental, theoretical and numerical methodologies. The work of O’Reilly et al. [4] is among the outstanding experimental works, in which the authors studied sliding elliptical bubbles through inclined planes. There are also valuable theoretical approaches to those problems, for instance, Podvin et al. [5] developed a lubrication model and appropriately tested it against experiments. Numerical methods are also a common approach to study bubbles/drops in complex geometries. In this sense, there are also some remarkable works; for instance, Han and Tryggvason [6] studied the breakup of drop falling against a solid wall by using a Front-Tracking technique.

There are three techniques which play an important role in the method proposed in the present paper: multiphase flows, dynamic meshes and immersed boundaries. First, multiphase flows could be treated numerically by using several approaches. In this paper, we use a Conservative Level Set (CLS) method developed by Balcázar et al. [7], which minimizes the problems of mass conservation of the standard level-set methods. This method was thoroughly verified in Refs [8, 9]. On the other hand, the use of dynamic meshes in the context of multiphase flows has received limited attention. Recently, Gutiérrez et al. [10] successfully used a dynamic mesh method to follow the ascent of bubbles in a CLS framework, thereby studying the challenging problem of the three-dimensional Taylor bubble. The use of this kind of meshes in bubbles/drops problems allows a significant reduction in the computational domain, thus saving computational resources and potentially allowing an improvement in the mesh resolution. Finally, although Immersed boundary (IB) methods have been broadly used in CFD [11], their applicability to multiphase flows remains almost unexplored [12, 13]. Those methods enable the treatment of flows with embedded boundaries on grids that do not conform to the shape of these boundaries. To the best of the author’s knowledge, the present paper is the first work that integrates an IB method in a CLS framework.

By assuming incompressible flow, two Newtonian fluids, an Arbitrary Lagrangian-Eulerian framework, no mass transfer at the interface between fluids, constant surface tension coefficient s and an embedded solid, the Navier–Stokes equations controlling the movement of a multiphase flow are given by the conservation laws for mass and momentum:

where $\mathbf{v}$ is the velocity field, $t$ is the time, $\rho$ is the fluid density, $\mu$ is its viscosity, $\mathbf{v}_{\text {mesh }}$ is the mesh velocity, $p$ is the pressure field, $\mathbf{g}$ is the gravity acceleration, $\mathbf{n}$ is the unit normal vector to the interface, $\kappa$ is its curvature and $\delta$ is the Dirac delta function located at the interface. $\Psi_{\rho_0}=-\rho_0 \mathbf{g}$ represents an extra source term needed to compensate the weight of the fluids within the [8], where $\rho_0=\frac{1}{V_{.}} \int \varphi \rho \mathrm{d} V, V_{\Omega}$ is the domain's volume and $\varphi$ is the smoothed Heavisade function of the Immersed Boundary method (see Section 2.3). Finally, $\Psi_{\mathrm{IB}}$ is another extra source term introduced by the Immersed Boundary method (see Section 2.3).

A CLS method is enforced to tackle the fluids interface [7]. In this method, the interface is represented by the indicator function $\phi$ :

where e is a parameter that defines the interface thickness. The interface between fluids can be located by getting the f = 0.5 isosurface. Thus, density and viscosity are obtained as:

where subscripts 1 and 2 refer to the suspending fluid and to the bubble/drop fluid, respectively. The level set function is advected by the velocity field:

After advection, a reinitialization equation is needed in order to maintain a constant interface thickness [14]:

where $\tau$ is the pseudo-time.

The calculus of the curvature $\kappa$ and the application of the pressure jump to the flow domain is carried out by means of the CSF model [15], which allows the conversion of the singular term $\sigma \kappa \mathbf{n} \delta_{\Gamma}$ into a volume force:

where $\kappa(\phi)$ and $\mathbf{n}$ are given by:

Here, $\nabla \phi$ is computed by means of a least-square method [7].

Governing equations have been discretized in a collocated unstructured grid arrangement by using a finite-volume method. We used the well-known Fractional Step method [7, 16] to solve the pressure-velocity coupling. Thus, momentum equation (Eqn. 2) is calculated by computing the following two steps:

where v* is the predictor velocity, and superscript n denotes that the corresponding variable is evaluated at the node n under consideration. By taking into account the continuity equation (eqn. 1), pressure can be solved by means of the following Poisson equation:

The discretization of this equation leads to a linear system, which is solved by means of a preconditioned conjugate gradient method. The present method has been implemented in the context of a parallel c++/MPI code called Termofluids. See Ref. [7] for further details on the numerical implementation and on the finite volume discretization of the governing equations.

The mesh is linearly moved at the vertical velocity of the drop. Therefore, the drop apparently is not vertically moved, although is deformed. The needed boundary conditions for this arrangement are explained in Section 3. We use the so-called Space Conservation Law (SCL) [17] in order to preserve the total computational volume, which is a needed condition when an Arbitrary Lagrangian-Eulerian framework is used. This law results in a modification of the mass fluxes through the cell faces by taking into account the corresponding swept volume of the mesh movement. As this mesh is linearly moved, the volume correction is straightforwardly computed.

The immersed body is represented by means of a triangular surface mesh in stereo-lithography format (STL). This allows intricate geometries to be handled. For the sake of simplicity, we assume a single solid within the domain, but the changes needed to extend the formulation to the case of multiple solids are straightforward. The solid can be moving at a given velocity V_s. A signed minimum distance field is defined in order to classify the different nodes of simulation domain in interior, exterior and forcing point, as seen in Fig. 1. Additionally, $\varphi(\mathbf{x}, t)$ defined as a smoothed Heavisade function that takes the value of 0 inside the solid, and 1 far in the fluid domain [18]. As the relative position of the mesh and the solid can change at each time-step (either because the solid moves or because the mesh does), this function should be computed at each iteration.

The needed extra source term $\Psi_{\mathrm{IB}}$ of the momentum equation is locally computed by using a direct forcing approach [19]:

where $\mathbf{V}$ is the desired value of the velocity field. It is directly computed at interior nodes, since the velocity of the solid is known. An approximation for the case of forcing points is used. In this case, $\mathbf{V}$ is computed by means of a second-order interpolation among the local velocity of the solid $\mathbf{V}_s$ and the predictor velocity $\mathbf{v}_{\Psi_{\mathrm{IB}}=0}^*$ of neighbor nodes calculated when $\Psi_{\mathrm{IB}}=0$. The reader is referred to Ref. [19] for further details about the calculation of $\Psi_{\mathrm{IB}}$.

A CFL condition is applied to compute the admissible time step at each iteration, in order to get stable simulations. By comparing the different terms of the momentum equation (eqn. 2), the following condition is obtained:

where $h$ is the characteristic cell size calculated as the cubic root of its volume, and $C_{\text {CFL }} \approx 0.1$ is a safety constant. Note that the Immersed Boundary source term $\Psi_{\text {IB }}$ does not introduce an additional restriction to the time step. This is because at the moment of computing the time step, the effect of the embedded body has already been taken into account in the calculus of the velocity field, so the convective restriction (first term of the right-hand part of eqn. (15)) already includes it.

The calculation procedure needed to advance from the current time instant tn to the next one t^{n + 1} is:

1. Calculate the mesh velocity, as explained in Section 2.2.

2. Choose a suitable time step, as explained in Section 2.4.

3. Advect the level set function f by solving eqn. 6.

4. Compress the interfaces between both fluids by solving eqn. 7.

5. Compute j , as explained in Section 2.3.

6. Update the density, viscosity, curvature and normal fields.

7. Solve eqns. 1 and 2 by using the fractional step method:

- Calculate the predictor velocity $\mathbf{v}_{\Psi_{\mathrm{IB}}=0}^*$ without considering the solid.

- Evaluate $\Psi_{\mathrm{IB}}$ as explained in Section 2.3.

- Calculate the final predictor velocity v*.

- Solve the Poisson equation (eqn. 13) to get the pressure field.

- Calculate the cell-face velocity [7]

8. Move the mesh and the solid (if needed).

9. Update mass fluxes by imposing the SCL.

10. Repeat the previous steps to reach the final time.

The method formulated above established a general technique to simulate bubbles/drops interacting with an arbitrary surface. Some of the major advantages are summarized hereunder. First, it constitutes a full three-dimensional method, so the axisymmetric hypothesis is not needed, being potentially able to handle any complex geometry. Second, the fact of including a domain optimization method (i.e. the moving mesh) allows substantial savings in the needed computational resources, both by minimizing the fluid domain and by allowing improved meshes. Finally, by using an STL-represented Immersed Boundary method, intricate surfaces can be easily reproduced. In that sense, the use of the IB method allows to actually compute the geometrical domain, forcing the simulation domain to face the real boundaries has it moves.

However, an important drawback should be beard in mind when more challenging cases are addressed by using this method. The fact of using open boundaries could difficult a proper applicability of this technique. First, the placement of the drop within the domain should not be arbitrary, and it should assure that the open boundaries do not alter the behavior of the drop. Furthermore, the treatment of the solution in the open boundaries could be difficult, for instance, in cases where the pressure profile is not well-defined over these boundaries.

The problem of a drop falling against a plane is addressed in subsequent sections, aiming to validate the method previously posed. As reference case, we chose one of the cases studied by Han and Tryggvason [6]. This reference case is defined by the following dimensionless numbers: $\eta_\rho=1.15, \eta_\mu=1, E o=12, O h_d=0.0466 . \eta_\rho=\rho_2 / \rho_1$ and $\eta_\mu=\mu_2 / \mu_1$ are the density and the viscosity ratios, respectively, $E o=g d^2\left(\rho_2-\rho_1\right) / \sigma$ is the Eötvös number, $O h_d=\mu_2 / \sqrt{\rho_2 d \sigma}$ is the Ohnesorge number, and $d$ is the initial droplet diameter.

A sketch of the initial arrangement is presented in Fig. 2a. The initial shape of the drop is a sphere of diameter d . The initial distance from the drop centre to the solid is set to 12d . We used that value since, based on previous studies available in the literature [20], it is enough in order to assure that the drop achieves its steady state before interacting with the solid. The lateral distance from the drop centroid to the lateral boundaries is fixed to 5d , since it gives rise to enough accurate results, as studied in Ref. [21]. The values of the distances from the drop to the inlet h_i = 2.9d and to outlet h_o = 5.1d are founded on a compromise between domain size and disturbance of the solution due to the proximity of those boundaries to the drop. See Ref. [10] for further notes about the setting of these magnitudes. No-slip boundary condition is imposed at the lateral walls. The used inflow and outflow boundary conditions are described in Ref. [10].

We used an unstructured mesh composed by tetrahedral control volumes. The drop is going to stay vertically steady at its initial position, and based on previous studies, it is known that at the selected regime the lateral movement of the drop is negligible. Therefore, we designed a mesh with a dense core of radius d , and a radial exponential growing in the size of the control volumes (see Fig. 2b). A mesh independence study is included in the results description, where the meshes presented in Table 1 were used.

In the present section, results of the gravity-driven falling drop problem are summarized. First, a comparison of the dimensionless terminal velocity $U_T^*=U_T(d g)^{-1 / 2}$ and the deformation parameter $\Delta=(L-B) /(L+B)$ at the falling state is presented in Table 2, where $L$ is the average length of the drop and $B$ is its average width. Good agreement were found in those results, especially for the finer mesh $M 3$, with errors of less than $5 \%$ in both magnitudes in comparison with reference data [6]. Furthermore, Fig. 3a shows the time evolution of the dimensionless velocity $U_T^*$. As this figure reveals, the proposed method properly capture the time in which the drop achieves to the solid (at a dimensionless time $t^*=\operatorname{tg}^{1 / 2} d^{-1 / 2}$ of around 45).

The time evolution of the vertical distance from the drop centre to the plane surface is plotted in Fig. 3b. The last part of this evolution reveals that the mesh M 3 is the sole one capable to maintain a constant distance from the drop interface to the wall. This shows that the coarser meshes do not properly capture the interaction of the drop with the plane, even though the falling behavior is correctly reproduced.

The mass conservation is evaluated in Table 3 by means of the percentage change in the drop volume. This table shows that, in that sense, the proposed CLS+MM+IB method behaves better than other methods available in the literature. This may be attributed to the smaller numerical errors of the CLS method, to a proper design and placement of the open boundaries, and to a bigger resolution in the region of interests. Note also that the references use axisymmetric solvers.

Streamlines and vorticity field are plotted in Figs. 4 and 5, obtained by using the results from M 3 mesh. First, the Fig. 4a shows the streamlines and vorticity in a plane perpendicular to the surface toward which the drop heads for. A single vortex is observed in the wake of the drop, close to the drop interface. Figure 4b shows that this vortex moves upwards and loses intensity when the drop is close to the solid. Additionally, two low-intensity counter-rotating vortices appear on the periphery of the drop. Figure 5a and b show the streamlines and vortic- ity field in a plane parallel to the floor, through the center of the drop. Both graphs highlight that the problem is intrinsically axisymmetric. Moreover, a clear pattern of pairs of low-intensity counter-rotating vortices can be observed in the vicinities of the interface, for both states.

Finally, the profiles evolution is presented in Fig. 6 in comparison with reference data [6]. Those profiles where obtained by using the mesh M 3. As can be observed, results from the present method qualitatively match the reference results.

In the present paper, a CLS-MM-IB method is formulated. This approach aims to couple the flexibility to represent intricate geometries of the Immersed Boundary method [19], the optimized domain obtained by using a Moving Mesh method [10], and the intrinsic advantages of mass conservation and robustness of the employed CLS method [7]. Therefore, the movement of bubbles/drops through full three-dimensional complex geometries could be addressed, with a very reasonable computational cost.

This method has been validated by solving the gravity-driven falling drop problem. Obtained results have been successfully compared with the reference data [6]. Reported results include velocity evolution, mesh independence study, evolution of the vertical position of the drop, streamlines and vorticity fields, and profiles evolution. These results present a very reasonable level of accuracy in comparison with reference results, and a consistent physical behavior.