# Remaillage des fronts lagrangiens dans les écoulements bi-fluides

Marcela Cruchaga, Diego Celentano, Piotr Breitkopf, Alain Rassineux, Pierre Villon

► To cite this version:

Marcela Cruchaga, Diego Celentano, Piotr Breitkopf, Alain Rassineux, Pierre Villon. Remaillage des fronts lagrangiens dans les écoulements bi-fluides. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812971

HAL Id: hal-01812971

https://hal.science/hal-01812971

Submitted on 12 Jun 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Remaillage des fronts lagrangiens dans les
écoulements bi-fluides

Marcela Cruchaga* — Diego Celentano* — Piotr Breitkopf** —
Alain Rassineux** — Pierre Villon**

* Departamento de Ingeniería Mecánica

Universidad de Santiago Chile, Av. Bdo. O'Higgins 3363, Santiago, Chile

e-mail: mcruchag@lauca.usach.clr

** Laboratorie Roberval UTC-CNRS, FRE 2833

Université de Téchnologie de Compiègne, Compiègne, France

e-mail: {prénom.nom}@utc.fr

RÉSUMÉ. La technique MLIT (Moving Lagrangian Interface Technique) utilise un maillage fixe pour la résolution du problème de dynamique des fluides, l'interface entre les deux fluides étant modélisée par un ensemble de marqueurs mobiles. Toutefois, comme un schéma lagrangien est employé pour avancer la position de l'interface, la distribution des marqueurs varie au cours du calcul et peut nécessiter des remaillages. Nous proposons ici une technique de placement des nouveaux marqueurs, basée sur l'approximation diffuse de la courbure. La technique proposée intègre la contrainte de préservation du volume et tient compte des restrictions géométriques. Cette approche est illustrée par un exemple de mouvement propre d'un liquide dans un réservoir.

ABSTRACT. In the framework of fixed-mesh finite element approaches for the computation of the fluid dynamics of the problem, a moving lagrangian interface technique (MLIT) was proposed to describe the interface. Nevertheless, as a lagrangian scheme is used to move the interface, its updated position could be represented by a highly distorted distribution of markers. To avoid this problem, we propose in this work a new remeshing technique applied to the interface. The remeshing is performed according to a curvature-based criterion by using a diffuse approximation technique. Moreover, the new distribution of points needs to be volume preserving. This condition is included as a constraint in the remeshing algorithm and it is the novel contribution of this work. Geometric restrictions, size distribution quality and physical aspects are all considered in the proposed formulation. Particular features of their capabilities are evaluated in simple test: a sloshing problem.

MOTS-CLÉS: Interfaces fluide-fluide, remaillage, Approximation Diffuse.

KEYWORDS:Fluid-fluid interfaces, remeshing, Diffuse Approximation.
**1. Introduction**

Many engineering applications require the analysis of time-dependent incompressible flows with moving two-liquid interfaces or free surfaces. Filling processes or open channels flows give an idea of the wide range of problems involving moving boundaries and interfaces. Many investigations have been devoted to numerical developments of formulations able to deal with the treatment of this kind of problems. Either moving or fixed domain discretizations are adopted to describe problems involving interfaces. All of them provide alternative solutions to update the material front and to overcome severe difficulties present in the numerical simulations of moving interfaces. In addition to the inconveniences related to the numerical schemes used to solve the fluid dynamics equations usually coming from the incompressibility constraint and relevant convection terms, there are many others due to the presence of two fluids and to the algorithm used to move the interface. In particular, global mass preserving and discontinuity in material properties are two of the most important aspects to be properly described in order to avoid numerical oscillations and to produce accurate results.

**2. The mass preserving algorithm for the remeshed interface**

A mass conservation algorithm applied to the remeshed interface is presented in this section. A remeshed interface position given by $(x_{\text{rem}}, y_{\text{rem}})$ is not, in general, a mass preserving configuration. Then, a further correction in the front coordinates, denoted as $(^{\text{P}}x_{\text{rem}}, ^{\text{P}}y_{\text{rem}})$, needs to be performed in order to attain a known given volume $\Omega_1$ (Figure 1.a).

<figure><img src="image_1.png" /><figcaption>**Figure 1.** *Mass preserving algorithm: a) interface description, b) move directions.*</figcaption></figure>

The coordinates of each node j are $(x_{\text{rem}}(j), y_{\text{rem}}(j))$ for $j=1,...,n_{\text{rem}}$, where $n_{\text{rem}}$ is the number of nodes. In these conditions, an 1D isoparametric representation of the interface within the interval $[j,j+1]$ is given by:
$$y_{\text{rem}}(\theta) = \sum_{k=1,2} \bar{N}_k(\theta) y_{\text{rem}}^k, \quad x_{\text{rem}}(\theta) = \sum_{k=1,2} \bar{N}_k(\theta) x_{\text{rem}}^k \qquad (1)$$

where $-1 \le \theta \le 1$ the isoparametric coordinate associated to the shape function matrix $\bar{N}$ ($\bar{N}_1 = 0.5(1-\theta)$ and $\bar{N}_2 = 0.5(1+\theta)$) such that $x_{\text{rem}}^1 = x_{\text{rem}}(j)$, $x_{\text{rem}}^2 = x_{\text{rem}}(j+1)$ and $y_{\text{rem}}^1 = y_{\text{rem}}(j)$, $y_{\text{rem}}^2 = y_{\text{rem}}(j+1)$.

The computation of the grey area in Figure 1.a is performed through the Green-Riemann formula

$$\text{Area} = \sum_{e=1}^{n_{\text{rem}}-1} \int_{-1}^{1} y_{\text{rem}}(\theta) x'_{\text{rem}}(\theta) d\theta = X_{\text{rem}}^{\text{T}} Q Y_{\text{rem}} \qquad (2)$$

where $X_{\text{rem}}^{\text{T}} = (x_{\text{rem}}(1), \dots, x_{\text{rem}}(n_{\text{rem}}))$, $Y_{\text{rem}}^{\text{T}} = (y_{\text{rem}}(1), \dots, y_{\text{rem}}(n_{\text{rem}}))$ and the matrix Q is computed as the assembly of the $(n_{\text{rem}}-1)$ segment contributions individually given by:

$$Q_{ij}^{e} = \int_{-1}^{1} \bar{N}_{i} \bar{N}_{j}^{\prime} d\theta \qquad (3)$$

resulting, for the adopted linear shape functions, in the following expression:

$$Q = \begin{bmatrix} -1/2 & 1/2 & 0 & \dots & 0 & 0 & 0 \\ -1/2 & 0 & 1/2 & \dots & 0 & 0 & 0 \\ 0 & \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & \dots & \dots & \dots & -1/2 & 0 & 1/2 \\ 0 & \dots & \dots & \dots & 0 & -1/2 & 1/2 \end{bmatrix} \qquad (4)$$

of a tridiagonal matrix whose components are independent of the distance between markers.

The volume to preserve can be written, according to equation (2), as:

$$^{\text{vp}}X_{\text{rem}}^{\text{T}} Q \quad ^{\text{vp}}Y_{\text{rem}} = \Omega_1 \qquad (5)$$

considering that the volume preserving new coordinates are given by:

$$^{\text{vp}}X_{\text{rem}} = X_{\text{rem}} + \beta \delta X_{\text{rem}}, \quad ^{\text{vp}}Y_{\text{rem}} = Y_{\text{rem}} + \beta \delta Y_{\text{rem}} \qquad (6)$$

$\delta X_{\text{new}}$ and $\delta Y_{\text{new}}$ being the admissible variations of the coordinates and $\beta$ is a scalar parameter. These variations are assumed to be described by:

$$\delta x_{\text{rem}}(j) = -\text{sign}[\mathbf{V}(j) \cdot \mathbf{n}(j)]u(j)\delta t, \quad \delta y_{\text{rem}}(j) = -\text{sign}[\mathbf{V}(j) \cdot \mathbf{n}(j)]v(j)\delta t \qquad (7)$$

where $\mathbf{V}=(u,v)$. The expressions (7) are proposed on the basis that the admissible motion of the markers is in the direction of their velocities and, besides, they must
compensate the loss (or gain) of mass flow through the interface. It is important to remark that with this procedure a physically-based correction is attained. Then, by introducing equations (6) and (7) into (5), we have the following quadratic equation on the $\beta$ parameter:

$$a\beta^2 + b\beta + c = 0 \qquad (8)$$

where the values of a, b and c are:

$$ 
\begin{aligned}
a &= \left\{ -\delta x_{\text{rem}}^{1} \left[ \delta y_{\text{rem}}^{1} + \delta y_{\text{rem}}^{2} \right] + \delta x_{\text{rem}}^{n_{\text{rem}}} \left[ \delta y_{\text{rem}}^{n_{\text{rem}}-1} + \delta y_{\text{rem}}^{n_{\text{rem}}} \right] + \sum_{j=2}^{n_{\text{rem}}-1} \delta x_{\text{rem}}^{j} \left[ \delta y_{\text{rem}}^{j-1} - \delta y_{\text{rem}}^{j+1} \right] \right\} / 2 \\
b &= \left\{ x_{\text{rem}}^{n_{\text{rem}}} \left[ \delta y_{\text{rem}}^{n_{\text{rem}}-1} + \delta y_{\text{rem}}^{n_{\text{rem}}} \right] + \delta x_{\text{rem}}^{n_{\text{rem}}} \left[ y_{\text{rem}}^{n_{\text{rem}}-1} + y_{\text{rem}}^{n_{\text{rem}}} \right] - x_{\text{rem}}^{1} \left[ \delta y_{\text{rem}}^{1} + \delta y_{\text{rem}}^{2} \right] \right. \\
& \quad \left. - \delta x_{\text{rem}}^{1} \left[ y_{\text{rem}}^{1} + y_{\text{rem}}^{2} \right] + \sum_{j=2}^{n_{\text{rem}}-1} \left( x_{\text{rem}}^{j} \left[ \delta y_{\text{rem}}^{j-1} - \delta y_{\text{rem}}^{j+1} \right] + \delta x_{\text{rem}}^{j} \left[ y_{\text{rem}}^{j-1} - y_{\text{rem}}^{j+1} \right] \right) \right\} / 2 \\
c &= \left\{ -x_{\text{rem}}^{1} \left[ \delta y_{\text{rem}}^{1} + \delta y_{\text{rem}}^{2} \right] + x_{\text{rem}}^{n_{\text{rem}}} \left[ \delta y_{\text{rem}}^{n_{\text{rem}}-1} + \delta y_{\text{rem}}^{n_{\text{rem}}} \right] \right. \\
& \quad \left. + \sum_{j=2}^{n_{\text{rem}}-1} \left( x_{\text{rem}}^{j} \left[ \delta y_{\text{rem}}^{j-1} - \delta y_{\text{rem}}^{j+1} \right] + \delta x_{\text{rem}}^{j} \left[ y_{\text{rem}}^{j-1} - y_{\text{rem}}^{j+1} \right] \right) \right\} / 2 - \Omega_{1}
\end{aligned}
\qquad (9)
$$

The solution of equation (8) results in two real roots (the variations (7) were found to always guarantee the condition $b^2 - 4ac \ge 0$). In practice, $\beta$ is chosen as the root with a minimum absolute value. With this parameter, equation (6) gives a mass preserving interface position. This approach provides a remeshed interface position with zero error in the volume. No iterations are needed to achieve convergence. This mass conservation corrector algorithm can be applied between two different configurations with different number of markers. The technique can be also extended to 3D where in this case the expression equivalent to (8) is of third degree on $\beta$.

## 2 Numerical example: Two-liquid interface problem

A closed container (Figure 2.a) with dimensions 0.8 m x 0.6 m is filled with two fluids with the lighter one being on top of the heavier one. The initial, inclined interface is linear with an average height of 0.3 m and slope 1:4. The fluid properties [1] used in the present analysis are $\rho_1 = 2$ kg/m$^3$, $\rho_2 = 1$ kg/m$^3$, $\mu_1 = \mu_2 = 0.001$ kg/ms with g=0.294 m/s$^2$ (the effect of different set of properties has been assessed in [49]). The boundary conditions consist of no-slip conditions at the upper and lower walls, and slip conditions at the side walls. Both fluids are initially at rest. Two different meshes are used (Figures 2.b and 2.c) to test the performance of the improved formulation when, in particular, coarse meshes are used. Note that the meshes in Figure 2 are refined in the zone where the movement of the interface is
expected to take place. The computations are performed with two different time
steps of sizes 0.5 s and 0.1 s.

The time history of the relative wave height at the left side of the container (i.e., (hA – 0.3)/0.3, hA being the instantaneous wave height at point A) obtained using different numerical approaches is presented in Figure 3. The results computed with both methods using the finer mesh and the smaller time step are practically identical. These results also adjust to those computed in [7,8] that are considered as the reference ones. Moreover, when the coarse mesh and the smaller time step are considered, the MLIRT provides slightly better results than those corresponding to ETILT. Nevertheless, erroneous results are obtained when using the ETILT with greater time step for both meshes. On the other hand, the relative wave height at point A predicted by the MLIRT exhibits a good response also for the greater time step. This satisfactory behavior in relatively coarse both time and space discretizations is an attractive aspect of the proposed methodology. Further investigation on this subject needs to be done in more severe problems.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Two-fluid interface problem: a) problem description and meshes used in the analysis: b) fine and c) coarse (respectively identified as “fm” and “cm”).</figcaption></figure>

The interface profiles at different instants of the analysis with MLIRT are
presented in Figure 4. Although high interface curvatures are not developed, their
changes during the analysis are properly captured.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Two-fluid interface problem. Analysis of time and mesh size dependency: a) ETILT solutions, b) MLIRT solutions using the fine (fm) and coarse (cm) meshes and different time discretizations.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** *Two-fluid interface problem. MLIRT solutions for the interface position at different instants of the analysis.*</figcaption></figure>

**10. Conclusions**

In this work, a volume preserving remeshing technique based on Diffuse Approximation of curvatures aimed at redefining the markers’ distribution on the interface has been described. The remeshing scheme is embedded in the two-fluid dynamic computation. The numerical predictions of the interface position during time have been found to be nearly independent of the number of initial markers and their distribution.

**11. Acknowledgements**

The support provided by the Chilean Council for Research and Technology CONICYT (FONDECYT Project N° 1020029), the Cooperation Project CONICYT -CNRS, the Department of Technological and Scientific Research at the Universidad de Santiago de Chile and the Université de Téchnologie de Compiègne (Visitor Professor Program) are all gratefully acknowledged.

**12. References**

[1] Tezduyar T.E., Behr M., and Liu J., “A new strategy for finite element computations involving moving boundaries and interfaces – the deforming-spatial-domain/space-time procedure, part I and II”, *Comp. Meth. in Appl. Mech. and Eng.*, 94(3), 339-371 (1992).

[2] Cruchaga M., Celentano D., and Tezduyar T., “A moving Lagrangian interface technique for flow computations over fixed meshes”, *Comp. Meth. in Appl. Mech. and Eng.*, 191, 525-543 (2001).