# A multi-material CCALE-MOF approach in cylindrical geometry

Marie Billaud-Friess, J. Breil, Stephane Galera, P.-H. Maire, M. Shashkov

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Marie Billaud-Friess, J. Breil, Stephane Galera, P.-H. Maire, M. Shashkov. A multi-material CCALE-MOF approach in cylindrical geometry. Communications in Computational Physics, 2014, 15 (2), pp.330-364. 10.4208/cicp.190912.080513a. hal-01004851

HAL Id: hal-01004851

https://hal.science/hal-01004851

Submitted on 6 Jun 2024

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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A multi-material CCALE-MOF approach in cylindrical geometry

Marie Billaud-Friess, J. Breil, Stephane Galera, P.-H. Maire, M. Shashkov

► To cite this version:

Marie Billaud-Friess, J. Breil, Stephane Galera, P.-H. Maire, M. Shashkov. A multi-material CCALE-MOF approach in cylindrical geometry. Communications in Computational Physics, 2014, 15 (2), pp.330-364. 10.4208/cicp.190912.080513a. hal-01004851

HAL Id: hal-01004851

https://hal.science/hal-01004851v1

Submitted on 6 Jun 2024

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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A multi-material CCALE-MOF approach in cylindrical
geometry

Marie Billaud Friess<sup>a</sup>, Jérôme Breil<sup>b,*</sup>, Stéphane Galera<sup>b</sup>, Pierre-Henri
Maire<sup>c</sup>, Mikhail Shashkov<sup>d</sup>

<sup>a</sup>GEM, UMR CNRS 6183, Ecole Centrale Nantes, 1 rue de la No, BP 92101, 44321
NANTES Cedex 3, France

<sup>b</sup>Univ. Bordeaux, CEA, CNRS, CELIA, UMR5107, F-33400 Talence, France

<sup>c</sup>CEA CESTA, BP 2, 33114 Le Barp Cedex, France

<sup>d</sup>Los Alamos National Laboratory, XCP-4, Los Alamos, NM 87545, USA

Abstract

In this paper we present recent developments concerning a Cell-Centered
Arbitrary Lagrangian Eulerian (CCALE) strategy using the Moment Of
Fluid (MOF) interface reconstruction for the numerical simulation of multi-
material compressible fluid flows on general unstructured grids in cylindri-
cal geometries. Especially, our attention is focused here on the following
points. First, we propose a new formulation of the scheme used during the
Lagrangian phase in the particular case of axisymmetric geometries. Then,
the MOF method is considered for multi-interface reconstruction in cylin-
drical geometry. Subsequently, a method devoted to the rezoning of polar
meshes is detailed. Finally, a generalization of the hybrid remapping to cylin-
drical geometries is presented. These explorations are validated by mean of
several test cases that clearly illustrate the robustness and accuracy of the

<footnote>*Corresponding author
Email addresses: **marie.billaud-friess@ec-nantes.fr** (Marie Billaud Friess),
**breil@celia.u-bordeaux1.fr** (Jérôme Breil)</footnote>
# A multi-material CCALE-MOF approach in cylindrical geometry

Marie Billaud Friess<sup>a</sup>, Jérôme Breil<sup>b,*</sup>, Stéphane Galera<sup>b</sup>, Pierre-Henri Maire<sup>c</sup>, Mikhail Shashkov<sup>d</sup>

<sup>a</sup>GeM, UMR CNRS 6183, Ecole Centrale Nantes, 1 rue de la No, BP 92101, 44321 NANTES Cedex 3, France

<sup>b</sup>Univ. Bordeaux, CEA, CNRS, CELIA, UMR5107, F-33400 Talence, France

<sup>c</sup>CEA CESTA, BP 2, 33114 Le Barp Cedex, France

<sup>d</sup>Los Alamos National Laboratory, XCP-4, Los Alamos, NM 87545, USA

## Abstract

In this paper we present recent developments concerning a Cell-Centered Arbitrary Lagrangian Eulerian (CCALE) strategy using the Moment Of Fluid (MOF) interface reconstruction for the numerical simulation of multi-material compressible fluid flows on general unstructured grids in cylindrical geometries. Especially, our attention is focused here on the following points. First, we propose a new formulation of the scheme used during the Lagrangian phase in the particular case of axisymmetric geometries. Then, the MOF method is considered for multi-interface reconstruction in cylindrical geometry. Subsequently, a method devoted to the rezoning of polar meshes is detailed. Finally, a generalization of the hybrid remapping to cylindrical geometries is presented. These explorations are validated by mean of several test cases that clearly illustrate the robustness and accuracy of the

<footnote>*Corresponding author
Email addresses: **marie.billaud-friess@ec-nantes.fr** (Marie Billaud Friess), **breil@celia.u-bordeaux1.fr** (Jérôme Breil)</footnote>
new method.

*Keywords:* Cell-centered scheme, Lagrangian hydrodynamics, ALE, MOF interface reconstruction, Rezoning algorithm, polar meshes, axisymmetric geometries.

# 1. Introduction

<figure><img src="image_3.png" /><figcaption>Figure 1: Multi-material CCALE algorithm flowchart.</figcaption></figure>

In this work, we consider the simulation of multi-material compressible flows on unstructured meshes in cylindrical geometry. For this, we adopt an ALE description [14] that has the great advantage to combine the best features of both Eulerian and Lagrangian approaches. Indeed, this choice is not only well adapted to naturally track free surfaces and interfaces between different fluids as purely Lagrangian methods, but also to handle flow distortion as Eulerian methods. Here, a CCALE [10, 11] approach is particularly
new method.

*Keywords:* Cell-centered scheme, Lagrangian hydrodynamics, ALE, MOF interface reconstruction, Rezoning algorithm, polar meshes, axisymmetric geometries.

# 1. Introduction

<figure><img src="image_2.png" /><figcaption>Figure 1: Multi-material CCALE algorithm flowchart.</figcaption></figure>

In this work, we consider the simulation of multi-material compressible flows on unstructured meshes in cylindrical geometry. For this, we adopt an ALE description [14] that has the great advantage to combine the best features of both Eulerian and Lagrangian approaches. Indeed, this choice is not only well adapted to naturally track free surfaces and interfaces between different fluids as purely Lagrangian methods, but also to handle flow distortion as Eulerian methods. Here, a CCALE [10, 11] approach is particularly
considered whose the main elements are as follow.

As depicted on figure FIG.1, the first step of the algorithm relies on an explicit Lagrangian phase in which the physical variables and grid are updated thanks to a slightly modified version of the Explicite Unstructured Cell-Centered Lagrangian HYDrodynamics (EUCCLHYD) scheme [19, 20, 21] in cylindrical coordinates. Recently, new investigations have been made about cell-centered Lagrangian schemes [3, 7]. The scheme presented in this paper is a modified version of the area weighted finite volume scheme of [19]. Then, multi-material flows treatment is done thanks to specific interface capturing method. This choice allows to track the volume fraction of each material used for the thermodynamical closure relying on the equal strain rates assumption. This approach is quite simple to implement and to use and remains sufficient in almost cases [11, 24]. This, leads to constant evolution of the volume fraction during the Lagrangian phase. Such an approach allows to reconstruct with accuracy the interface between each material. In this context, many development have been done for 2D Cartesian geometries. First, a previous version of the CCALE algorithm solving two-material compressible flows using a Volume Of Fluid (VOF) have been proposed in [6, 11]. Then an extension to Moment Of Fluid (MOF) approach has been considered to enhance multi-material (more than two components) flows in [8, 10]. Subsequently, a rezoning phase is realized. It consists in moving the Lagrangian nodes to improve the geometric quality of the grid [15]. Finally, the physical variables are conservatively interpolated from the Lagrangian grid onto the new rezoned one during the remapping phase. Here an extension of the hybrib remapping [4] to cylindrical geometries is introduced. We want to
considered whose the main elements are as follow.

As depicted on figure FIG.1, the first step of the algorithm relies on an explicit Lagrangian phase in which the physical variables and grid are updated thanks to a slightly modified version of the Explicite Unstructured Cell-Centered Lagrangian HYDrodynamics (EUCCLHYD) scheme [19, 20, 21] in cylindrical coordinates. Recently, new investigations have been made about cell-centered Lagrangian schemes [3, 7]. The scheme presented in this paper is a modified version of the area weighted finite volume scheme of [19]. Then, multi-material flows treatment is done thanks to specific interface capturing method. This choice allows to track the volume fraction of each material used for the thermodynamical closure relying on the equal strain rates assumption. This approach is quite simple to implement and to use and remains sufficient in almost cases [11, 24]. This, leads to constant evolution of the volume fraction during the Lagrangian phase. Such an approach allows to reconstruct with accuracy the interface between each material. In this context, many development have been done for 2D Cartesian geometries. First, a previous version of the CCALE algorithm solving two-material compressible flows using a Volume Of Fluid (VOF) have been proposed in [6, 11]. Then an extension to Moment Of Fluid (MOF) approach has been considered to enhance multi-material (more than two components) flows in [8, 10]. Subsequently, a rezoning phase is realized. It consists in moving the Lagrangian nodes to improve the geometric quality of the grid [15]. Finally, the physical variables are conservatively interpolated from the Lagrangian grid onto the new rezoned one during the remapping phase. Here an extension of the hybrib remapping [4] to cylindrical geometries is introduced. We want to
notice that in ALE framework using cell-centered formulation, this phase is
straightforward. In the lines of these works, the main goal of this paper is
to extend the CCALE-MOF algorithm to treat both Cartesian and cylindri-
cal geometry. To this end, several modifications are given to the algorithm
previously presented. In a first part, we propose a new formulation of the
numerical scheme introduced in [19] for treating axisymmetric geometries
during the Lagrangian phase. To build this scheme, an area-weighted formu-
lation of the Lagrangian system of equations is proposed. Then, this system
of equations is discretized using a cell-centered finite volume (FV) scheme.
Contrary to [19] in which fluxes are directly deduced from the Geometric
Conservation Law (GCL) constraint, here a simpler formulation that gives
similar results is retained. These two main choices lead to a robust first-
order scheme conservative for the total energy that has the great advantage
to preserve spherical symmetry for one-dimensional flow on uniform angular
polar grids. The high order extension has been performed using the Gener-
alized Riemann Problem (GRP) described in [19]. To treat interface flows, a
MOF interface reconstruction method is retained in the sequel. Once again,
the difficulty here is to propose a natural and consistent adaptation of this
approach able to treat axisymmetric interface flows. To this end, formula-
tions of the moments needed to track interface are revisited for cylindrical
coordinates as in [2]. This leads to an accurate and second order interface
reconstruction method that allows to treat multi-material (more than two)
interfaces in the lines of [8]. The third part of this study is dedicated to
recent enhancement of the rezoning algorithm to improve the mesh quality
during computation especially on polar meshes. As it is done in [10, 11],
notice that in ALE framework using cell-centered formulation, this phase is
straightforward. In the lines of these works, the main goal of this paper is
to extend the CCALE-MOF algorithm to treat both Cartesian and cylindri-
cal geometry. To this end, several modifications are given to the algorithm
previously presented. In a first part, we propose a new formulation of the
numerical scheme introduced in [19] for treating axisymmetric geometries
during the Lagrangian phase. To build this scheme, an area-weighted formu-
lation of the Lagrangian system of equations is proposed. Then, this system
of equations is discretized using a cell-centered finite volume (FV) scheme.
Contrary to [19] in which fluxes are directly deduced from the Geometric
Conservation Law (GCL) constraint, here a simpler formulation that gives
similar results is retained. These two main choices lead to a robust first-
order scheme conservative for the total energy that has the great advantage
to preserve spherical symmetry for one-dimensional flow on uniform angular
polar grids. The high order extension has been performed using the Gener-
alized Riemann Problem (GRP) described in [19]. To treat interface flows, a
MOF interface reconstruction method is retained in the sequel. Once again,
the difficulty here is to propose a natural and consistent adaptation of this
approach able to treat axisymmetric interface flows. To this end, formula-
tions of the moments needed to track interface are revisited for cylindrical
coordinates as in [2]. This leads to an accurate and second order interface
reconstruction method that allows to treat multi-material (more than two)
interfaces in the lines of [8]. The third part of this study is dedicated to
recent enhancement of the rezoning algorithm to improve the mesh quality
during computation especially on polar meshes. As it is done in [10, 11],
mesh rezoning is based on the Condition Number Smoothing (CNS) [15] algorithm on unstructured meshes. Moreover, when used for polar meshes, it is well known that CNS algorithm pushes the nodes toward the origin deteriorating the mesh quality. To avoid this drawback, the main idea developed in this paper is to adapt CNS algorithm to polar grids. Then, extension to unstructured grids (Cartesian-polar) is also explored. Finally, a generalization of the remapping procedure to cylindrical geometries is proposed. Here, an efficient method adapted to multi-material flows is presented. The main idea is to use a hybrid remapping that combine the main advantages of the swept-face and multi-material cell-intersection remapping as in [4, 10]. Finally, a specific attention is done to polynomial integration that preserves the method efficiency.

The paper is structured as follows. We detail in the second section a new formulation of the first-order area weighted Lagrangian scheme used for axisymmetric geometries. Further extensions to high-order are notably detailed in [19]. Afterwards, the extension of the MOF axisymmetric interface reconstruction method is presented for treating multimaterial flows. Then, we describe the General Condition Number Smoothing (GCNS) algorithm for unstructured meshes. Finally, the description of the new hybrid remapping procedure for cylindrical geometry is done. For a complete description of the CCALE-MOF method see [10, 11], except new advances presented in this paper. Then presentation of numerical experiments is made in Section 4. They demonstrate not only the robustness and the accuracy of the present methodology but also its ability to handle successfully complex two-
mesh rezoning is based on the Condition Number Smoothing (CNS) [15] algorithm on unstructured meshes. Moreover, when used for polar meshes, it is well known that CNS algorithm pushes the nodes toward the origin deteriorating the mesh quality. To avoid this drawback, the main idea developed in this paper is to adapt CNS algorithm to polar grids. Then, extension to unstructured grids (Cartesian-polar) is also explored. Finally, a generalization of the remapping procedure to cylindrical geometries is proposed. Here, an efficient method adapted to multi-material flows is presented. The main idea is to use a hybrid remapping that combine the main advantages of the swept-face and multi-material cell-intersection remapping as in [4, 10]. Finally, a specific attention is done to polynomial integration that preserves the method efficiency.

The paper is structured as follows. We detail in the second section a new formulation of the first-order area weighted Lagrangian scheme used for axisymmetric geometries. Further extensions to high-order are notably detailed in [19]. Afterwards, the extension of the MOF axisymmetric interface reconstruction method is presented for treating multimaterial flows. Then, we describe the General Condition Number Smoothing (GCNS) algorithm for unstructured meshes. Finally, the description of the new hybrid remapping procedure for cylindrical geometry is done. For a complete description of the CCALE-MOF method see [10, 11], except new advances presented in this paper. Then presentation of numerical experiments is made in Section 4. They demonstrate not only the robustness and the accuracy of the present methodology but also its ability to handle successfully complex two-
dimensional multi-material fluid flows notably computed for axisymmetric geometries. Finally concluding remarks and perspectives about future works are given in the last section.

## 2. Lagrangian phase in axisymmetric geometry

In this part, an extension of the cell-centered Lagrangian scheme [20, 21] is presented for the numerical simulation of compressible flows in pseudo-Cartesian geometries for unstructured meshes as in [19]. This choice has the great advantage to treat both axisymmetric and Cartesian geometries. In this paper, a new and simple formulation of the scheme introduced in [19] for first-order approximation is proposed. To this end, an area weighted formulation of classical Lagrangian equations is first introduced. Then these equations are discretized with a node-centered approximate Riemann solver.

### 2.1. Governing equations

During the Lagrangian phase, the rates of change of volume, mass, momentum and total energy are computed assuming that discretized volumes move following the flow. Thus, each arbitrary volume $V(t)$ depending on the time $t > 0$ moves satisfying the following system of equations

$$ \frac{d}{dt} \int_V \rho dV = 0, \qquad (1) $$

$$ \frac{d}{dt} \int_V dV - \int_V \nabla \cdot \mathbf{U} dV = 0, \qquad (2) $$

$$ \frac{d}{dt} \int_V \rho \mathbf{U} dV + \int_V \nabla P dV = 0, \qquad (3) $$

$$ \frac{d}{dt} \int_V \rho E dV + \int_V \nabla \cdot (P \mathbf{U}) dV = 0, \qquad (4) $$
dimensional multi-material fluid flows notably computed for axisymmetric geometries. Finally concluding remarks and perspectives about future works are given in the last section.

## 2. Lagrangian phase in axisymmetric geometry

In this part, an extension of the cell-centered Lagrangian scheme [20, 21] is presented for the numerical simulation of compressible flows in pseudo-Cartesian geometries for unstructured meshes as in [19]. This choice has the great advantage to treat both axisymmetric and Cartesian geometries. In this paper, a new and simple formulation of the scheme introduced in [19] for first-order approximation is proposed. To this end, an area weighted formulation of classical Lagrangian equations is first introduced. Then these equations are discretized with a node-centered approximate Riemann solver.

### 2.1. Governing equations

During the Lagrangian phase, the rates of change of volume, mass, momentum and total energy are computed assuming that discretized volumes move following the flow. Thus, each arbitrary volume $V(t)$ depending on the time $t > 0$ moves satisfying the following system of equations

$$ \frac{d}{dt} \int_V \rho dV = 0, \qquad (1) $$

$$ \frac{d}{dt} \int_V dV - \int_V \nabla \cdot \mathbf{U} dV = 0, \qquad (2) $$

$$ \frac{d}{dt} \int_V \rho \mathbf{U} dV + \int_V \nabla P dV = 0, \qquad (3) $$

$$ \frac{d}{dt} \int_V \rho E dV + \int_V \nabla \cdot (P \mathbf{U}) dV = 0, \qquad (4) $$
where $\frac{d}{dt}$ is the Lagrangian derivative and $\rho, U, P, E$ are respectively the density, velocity, pressure and total energy. In addition, this system is closed thanks to an equation of state (EOS) as

$$p = p(\rho, \varepsilon), \qquad (5)$$

with the internal energy $\varepsilon$ defined as $\varepsilon = E - |U|^2/2$. At last, we have local kinematic equation

$$\frac{d\mathbf{X}}{dt} = \mathbf{U}, \quad \mathbf{X}(0) = \mathbf{X}^0, \qquad (6)$$

with $\mathbf{X}$ the location of a point of the control volume surface $S(t)$, at time $t > 0$ and $\mathbf{X}^0$ its initial value. This equation is equivalent to (2) also known as geometric conservation of law (GCL).

## 2.2. Area-weighted formulation

<figure><img src="image_4.png" /><figcaption>Figure 2: Notations related to the pseudo-cartesian grid.</figcaption></figure>

For defining the differential operators used in the system of Lagrangian equation (1)-(4) a pseudo-Cartesian reference frame $\{0, X, Y\}$ for the orthonormal basis $(e_X, e_Y)$ is used (see FIG.2). Thus each point is localized
where $\frac{d}{dt}$ is the Lagrangian derivative and $\rho, U, P, E$ are respectively the density, velocity, pressure and total energy. In addition, this system is closed thanks to an equation of state (EOS) as

$$p = p(\rho, \varepsilon), \qquad (5)$$

with the internal energy $\varepsilon$ defined as $\varepsilon = E - |U|^2/2$. At last, we have local kinematic equation

$$\frac{d\mathbf{X}}{dt} = \mathbf{U}, \quad \mathbf{X}(0) = \mathbf{X}^0, \qquad (6)$$

with $\mathbf{X}$ the location of a point of the control volume surface $S(t)$, at time $t > 0$ and $\mathbf{X}^0$ its initial value. This equation is equivalent to (2) also known as geometric conservation of law (GCL).

## 2.2. Area-weighted formulation

<figure><img src="image_3.png" /><figcaption>Figure 2: Notations related to the pseudo-cartesian grid.</figcaption></figure>

For defining the differential operators used in the system of Lagrangian equation (1)-(4) a pseudo-Cartesian reference frame $\{0, X, Y\}$ for the orthonormal basis $(e_X, e_Y)$ is used (see FIG.2). Thus each point is localized
by means of its positions $X$ and $\mathcal{R}(Y) = 1 - \alpha + \alpha Y$ the pseudo-radius. When $\alpha = 0$, the Lagrangian equations for Cartesian geometry are recovered, otherwise for $\alpha = 1$ this corresponds to axisymmetric equations. In this way, axisymmetric geometry is obtained from Cartesian one through a rotational symmetry about the $X$-axis. This implies that the volume $V(t)$ is generated by the rotation of the area $A(t)$ about the $X$-axis. In consequence, the element volume $dV$ writes as $dV = \mathcal{R}dA$ with $dA = dXdY$ the element area in the pseudo-Cartesian frame. In the same manner, the control surface $S(t)$ delimiting $V(t)$ is obtained through the rotation of $L(t)$ the boundary of $A(t)$ and the surface element is given by $dS = \mathcal{R}dL$. Note that we have omitted the $2\pi$ factor in the evaluation of the element volume.

In a such framework, the velocity divergence and the pressure gradient read as follows

$$ \nabla \cdot \mathbf{U} = \frac{1}{\mathcal{R}} \left[ \frac{\partial(\mathcal{R}u)}{\partial X} + \frac{\partial(\mathcal{R}v)}{\partial Y} \right], \quad \text{where } \mathbf{U}^t = (u, v) \qquad (7) $$

and

$$ \nabla P = \left( \frac{\partial P}{\partial X} \mathbf{e}_X + \frac{\partial P}{\partial Y} \mathbf{e}_Y \right). \qquad (8) $$

Using the previous definitions and after some calculations using the Green's formula, it is possible to rewrite (1)-(4) at least in two different ways. The first one, obtained without any approximation is the control volume formulation. When discretized this formulation leads to a conservative scheme for both equations of energy and momentum, and satisfies the local semi-discrete entropy inequality. However, as shown in [19] it does not preserve symmetries. Consequently, an area-weighted formulation is adopted here leading
by means of its positions $X$ and $\mathcal{R}(Y) = 1 - \alpha + \alpha Y$ the pseudo-radius. When $\alpha = 0$, the Lagrangian equations for Cartesian geometry are recovered, otherwise for $\alpha = 1$ this corresponds to axisymmetric equations. In this way, axisymmetric geometry is obtained from Cartesian one through a rotational symmetry about the $X$-axis. This implies that the volume $V(t)$ is generated by the rotation of the area $A(t)$ about the $X$-axis. In consequence, the element volume $dV$ writes as $dV = \mathcal{R}dA$ with $dA = dXdY$ the element area in the pseudo-Cartesian frame. In the same manner, the control surface $S(t)$ delimiting $V(t)$ is obtained through the rotation of $L(t)$ the boundary of $A(t)$ and the surface element is given by $dS = \mathcal{R}dL$. Note that we have omitted the $2\pi$ factor in the evaluation of the element volume.

In a such framework, the velocity divergence and the pressure gradient read as follows

$$ \nabla \cdot \mathbf{U} = \frac{1}{\mathcal{R}} \left[ \frac{\partial(\mathcal{R}u)}{\partial X} + \frac{\partial(\mathcal{R}v)}{\partial Y} \right], \quad \text{where } \mathbf{U}^t = (u, v) \qquad (7) $$

and

$$ \nabla P = \left( \frac{\partial P}{\partial X} \mathbf{e}_X + \frac{\partial P}{\partial Y} \mathbf{e}_Y \right). \qquad (8) $$

Using the previous definitions and after some calculations using the Green's formula, it is possible to rewrite (1)-(4) at least in two different ways. The first one, obtained without any approximation is the control volume formulation. When discretized this formulation leads to a conservative scheme for both equations of energy and momentum, and satisfies the local semi-discrete entropy inequality. However, as shown in [19] it does not preserve symmetries. Consequently, an area-weighted formulation is adopted here leading
to a conservative scheme for energy equation that respect spherical geometries. This formulation is deduced from the control volume one assuming that momentum equation (3) is written in Cartesian geometry. Like this, the area-weighted formulation for the Lagrangian equations reads

$$m \frac{d}{dt} \left\langle \frac{1}{\rho} \right\rangle - \int_L \mathbf{N} \cdot \mathcal{R} \mathbf{U} dL = 0, \qquad (9)$$

$$m \frac{d}{dt} \langle \mathbf{U} \rangle + \bar{\mathcal{R}} \int_L P N dL = 0, \qquad (10)$$

$$m \frac{d}{dt} \langle E \rangle + \int_L P N \cdot \mathcal{R} U dL = 0, \qquad (11)$$

where $m = \int_V \rho dV$ represents the mass of the volume $V$. Each physical variable per unit of mass ($E, U$) is noted as $\phi$, and has its mass density mean value defined by $\langle \phi \rangle = \frac{1}{m} \int_V \rho \phi dV$. The average $\bar{\mathcal{R}}$ corresponds to ratio $\bar{\mathcal{R}} = \frac{V}{A}$. In such case, as $m = \rho V$, the momentum equation is solved in Cartesian geometry. For Cartesian case $V = A$, we recover $\bar{\mathcal{R}} = 1$. Further details on the derivation of this system are available in [19].

## 2.3. Numerical scheme

Thereafter, we recall briefly the first order cell-centered Lagrangian scheme introduced in [19]. To this goal, similar notations as [10, 19, 22] are employed in the sequel. Let us consider a set $\{\Omega_c\}_{c \in \mathbb{N}}$ of non-overlapping polygonal cells that approximates $A(t)$. Each cell noted $\Omega_c$ is assigned a single index $c$. Each vertex of the cell $c$ is labeled with the index $p$ and is localized thanks to its coordinates $\mathbf{X}_p = (X_p, Y_p)^t$ in the pseudo-Cartesian frame. In addition, we introduce $\mathcal{P}(c)$ the list of the vertices belonging to the cell $\Omega_c$ and $\mathcal{C}(p)$ the list of the cells sharing the vertex $p$. These two sets are counterclockwise ordered.
to a conservative scheme for energy equation that respect spherical geometries. This formulation is deduced from the control volume one assuming that momentum equation (3) is written in Cartesian geometry. Like this, the area-weighted formulation for the Lagrangian equations reads

$$m \frac{d}{dt} \left\langle \frac{1}{\rho} \right\rangle - \int_L \mathbf{N} \cdot \mathcal{R} \mathbf{U} dL = 0, \qquad (9)$$

$$m \frac{d}{dt} \langle \mathbf{U} \rangle + \bar{\mathcal{R}} \int_L P N dL = 0, \qquad (10)$$

$$m \frac{d}{dt} \langle E \rangle + \int_L P N \cdot \mathcal{R} U dL = 0, \qquad (11)$$

where $m = \int_V \rho dV$ represents the mass of the volume $V$. Each physical variable per unit of mass ($E, U$) is noted as $\phi$, and has its mass density mean value defined by $\langle \phi \rangle = \frac{1}{m} \int_V \rho \phi dV$. The average $\bar{\mathcal{R}}$ corresponds to ratio $\bar{\mathcal{R}} = \frac{V}{A}$. In such case, as $m = \rho V$, the momentum equation is solved in Cartesian geometry. For Cartesian case $V = A$, we recover $\bar{\mathcal{R}} = 1$. Further details on the derivation of this system are available in [19].

## 2.3. Numerical scheme

Thereafter, we recall briefly the first order cell-centered Lagrangian scheme introduced in [19]. To this goal, similar notations as [10, 19, 22] are employed in the sequel. Let us consider a set $\{\Omega_c\}_{c \in \mathbb{N}}$ of non-overlapping polygonal cells that approximates $A(t)$. Each cell noted $\Omega_c$ is assigned a single index $c$. Each vertex of the cell $c$ is labeled with the index $p$ and is localized thanks to its coordinates $\mathbf{X}_p = (X_p, Y_p)^t$ in the pseudo-Cartesian frame. In addition, we introduce $\mathcal{P}(c)$ the list of the vertices belonging to the cell $\Omega_c$ and $\mathcal{C}(p)$ the list of the cells sharing the vertex $p$. These two sets are counterclockwise ordered.
Let us introduce $p^-$ and $p^+$ the previous and the next nodes with respect to p in $\mathcal{P}(c)$. We denote by $L_{pc}^-$, $L_{pc}^+$ the half length of the edges $[pp^-]$, $[pp^+]$. Similar notations are used for the normals outward $N_{pc}^+$ and $N_{pc}^-$. Finally, the corner normal $L_{pc}N_{pc}$ is given by $L_{pc}N_{pc} = L_{pc}^+N_{pc}^+ + L_{pc}^-N_{pc}^-$. All these notations have been displayed in Fig.3.

<figure><img src="image_5.png" /><figcaption>Figure 3: Notations for the cell-centered scheme.</figcaption></figure>

The first order spatial approximation of (9)-(11) is obtained considering local integrals on each cell $\Omega_c$ rotated about the X-axis. The mass $m_c$ of the cell $\Omega_c$ is $m_c = \int_{\Omega_c} \rho dV$ and each flow variable $\phi$ (as total energy, velocity) is averaged over each cell through the formula

$$ \phi_c = \frac{1}{m_c} \int_{\Omega_c} \rho \phi dV, $$
Let us introduce $p^-$ and $p^+$ the previous and the next nodes with respect to p in $\mathcal{P}(c)$. We denote by $L_{pc}^-$, $L_{pc}^+$ the half length of the edges $[pp^-]$, $[pp^+]$. Similar notations are used for the normals outward $N_{pc}^+$ and $N_{pc}^-$. Finally, the corner normal $L_{pc}N_{pc}$ is given by $L_{pc}N_{pc} = L_{pc}^+N_{pc}^+ + L_{pc}^-N_{pc}^-$. All these notations have been displayed in Fig.3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Notations for the cell-centered scheme.</figcaption></figure>

The first order spatial approximation of (9)-(11) is obtained considering local integrals on each cell $\Omega_c$ rotated about the X-axis. The mass $m_c$ of the cell $\Omega_c$ is $m_c = \int_{\Omega_c} \rho dV$ and each flow variable $\phi$ (as total energy, velocity) is averaged over each cell through the formula

$$ \phi_c = \frac{1}{m_c} \int_{\Omega_c} \rho \phi dV, $$
named cell-centered value. Then, we have

$$m_c \frac{d}{dt} \mathbf{U}_c + \bar{\mathcal{R}}_c \sum_{p \in \mathcal{P}(c)} \mathbf{F}_{pc} = \mathbf{0}, \qquad (12)$$

$$m_c \frac{d}{dt} E_c + \sum_{p \in \mathcal{P}(c)} \mathbf{F}_{pc} \cdot \mathcal{R}_p \mathbf{U}_p = 0. \qquad (13)$$

In addition, the mesh is moved through the local kinematic equation given at each node by

$$\frac{d\mathbf{X}_p}{dt} = \mathbf{U}_p \text{ for } t > 0 \quad \text{and} \quad \mathbf{X}_p(0) = \mathbf{X}_p^0, \qquad (14)$$

with $\mathbf{U}_p$ and $\mathbf{X}_p^0$ respectively the velocity and the position of a node $p$ at initial time. In the previous equations, $\mathbf{F}_{pc}$ is the numerical flux at each node $p$ of each cell $c$ defined by

$$\mathbf{F}_{pc} = L_{pc} P_c \mathbf{N}_{pc} - \mathbf{M}_{pc} (\mathbf{U}_p - \mathbf{U}_c), \qquad (15)$$

with $\mathbf{U}_p$ the velocity at the point $p$ and $P_c$ the mean value of the pressure in the cell $c$. The 2 × 2 matrices $\mathbf{M}_{pc}$ and $\mathbf{M}_p$ are defined as

$$\mathbf{M}_{pc} = \mathbf{Z}_c (L_{pc}^{-} \mathbf{N}_{pc}^{-} \otimes \mathbf{N}_{pc}^{-} + L_{pc}^{+} \mathbf{N}_{pc}^{+} \otimes \mathbf{N}_{pc}^{+}), \quad \text{and} \quad \mathbf{M}_p = \sum_{c \in \mathcal{C}(p)} \mathbf{M}_{pc}. \qquad (16)$$

Where, we introduce the “swept mass flux” [9] associated to the isentropic sound speed $a_c$ that is

$$\mathbf{Z}_c = \rho_c a_c. \qquad (17)$$

This is nothing but the acoustic impedance. As it has been demonstrated in [22] the total energy and momentum conservation is equivalent to

$$\sum_{c \in \mathcal{C}(p)} \mathbf{F}_{pc} = \mathbf{0}. \qquad (18)$$
named cell-centered value. Then, we have

$$m_c \frac{d}{dt} \mathbf{U}_c + \bar{\mathcal{R}}_c \sum_{p \in \mathcal{P}(c)} \mathbf{F}_{pc} = \mathbf{0}, \qquad (12)$$

$$m_c \frac{d}{dt} E_c + \sum_{p \in \mathcal{P}(c)} \mathbf{F}_{pc} \cdot \mathcal{R}_p \mathbf{U}_p = 0. \qquad (13)$$

In addition, the mesh is moved through the local kinematic equation given at each node by

$$\frac{d\mathbf{X}_p}{dt} = \mathbf{U}_p \text{ for } t > 0 \quad \text{and} \quad \mathbf{X}_p(0) = \mathbf{X}_p^0, \qquad (14)$$

with $\mathbf{U}_p$ and $\mathbf{X}_p^0$ respectively the velocity and the position of a node $p$ at initial time. In the previous equations, $\mathbf{F}_{pc}$ is the numerical flux at each node $p$ of each cell $c$ defined by

$$\mathbf{F}_{pc} = L_{pc} P_c \mathbf{N}_{pc} - \mathbf{M}_{pc} (\mathbf{U}_p - \mathbf{U}_c), \qquad (15)$$

with $\mathbf{U}_p$ the velocity at the point $p$ and $P_c$ the mean value of the pressure in the cell $c$. The 2 × 2 matrices $\mathbf{M}_{pc}$ and $\mathbf{M}_p$ are defined as

$$\mathbf{M}_{pc} = \mathbf{Z}_c (L_{pc}^{-} \mathbf{N}_{pc}^{-} \otimes \mathbf{N}_{pc}^{-} + L_{pc}^{+} \mathbf{N}_{pc}^{+} \otimes \mathbf{N}_{pc}^{+}), \quad \text{and} \quad \mathbf{M}_p = \sum_{c \in \mathcal{C}(p)} \mathbf{M}_{pc}. \qquad (16)$$

Where, we introduce the “swept mass flux” [9] associated to the isentropic sound speed $a_c$ that is

$$\mathbf{Z}_c = \rho_c a_c. \qquad (17)$$

This is nothing but the acoustic impedance. As it has been demonstrated in [22] the total energy and momentum conservation is equivalent to

$$\sum_{c \in \mathcal{C}(p)} \mathbf{F}_{pc} = \mathbf{0}. \qquad (18)$$
Finally using (15), the nodal velocity $U_p$ is deduced from (18) by solving the linear system

$$ M_p U_p = \sum_{c \in \mathcal{C}(p)} (L_{pc} P_c N_{pc} + M_{pc} U_c). \quad (19) $$

In [19], the numerical fluxes used in the discretization of (9) and (11) are chosen for satisfying the local GCL constraint (14). Here, we rather adopt a more simple approach that give similar results. Since (14) is explicitly solved for moving the mesh in time, there is no need to solve (9). Thus, each cell volume $V_c$ is directly deduced from (14). Thereby, it is possible to choose for the numerical flux in (11) a simple form as in (13) with $\mathcal{R}_p = 1 - \alpha + \alpha Y_p$. Concerning, the momentum equation the mean value $\mathcal{R}_c$ is equal to the discrete ratio $\bar{\mathcal{R}}_c = \frac{V_c}{A_c}$.

Let us note that this new formulation of the area-weighted discretization relies on a node-centered solver which is exactly the same as the one developed in [20] for two-dimensional Cartesian geometry. However, the present spatial discretization does not satisfy rigorously the GCL compatibility requirement. In what follows, we will assess the discrepancy of our discretization to the GCL by analyzing the corresponding discrete divergence operator. The discrete divergence operator that corresponds to the present scheme writes as

$$ (\nabla \cdot \mathbf{U})_c = \frac{1}{V_c} \sum_{p \in \mathcal{P}(c)} \mathcal{R}_p (L_{pc}^{-} N_{pc}^{-} + L_{pc}^{+} N_{pc}^{+}) \cdot \mathbf{U}_p, \quad (20) $$

where $\mathcal{R}_p$ denotes the pseudo-radius of vertex $p$. It is shown in [22, 19] that the discrete divergence operator deduced from the GCL reads

$$ (\nabla \cdot \mathbf{U})_{c}^{\text{GCL}} = \frac{1}{V_c} \sum_{p \in \mathcal{P}(c)} \frac{1}{3} [(2\mathcal{R}_p + \mathcal{R}_p^{-}) L_{pc}^{-} N_{pc}^{-} + (2\mathcal{R}_p + \mathcal{R}_p^{+}) L_{pc}^{+} N_{pc}^{+}] \cdot \mathbf{U}_p. \quad (21) $$
Finally using (15), the nodal velocity $U_p$ is deduced from (18) by solving the linear system

$$ M_p U_p = \sum_{c \in \mathcal{C}(p)} (L_{pc} P_c N_{pc} + M_{pc} U_c). \quad (19) $$

In [19], the numerical fluxes used in the discretization of (9) and (11) are chosen for satisfying the local GCL constraint (14). Here, we rather adopt a more simple approach that give similar results. Since (14) is explicitly solved for moving the mesh in time, there is no need to solve (9). Thus, each cell volume $V_c$ is directly deduced from (14). Thereby, it is possible to choose for the numerical flux in (11) a simple form as in (13) with $\mathcal{R}_p = 1 - \alpha + \alpha Y_p$. Concerning, the momentum equation the mean value $\mathcal{R}_c$ is equal to the discrete ratio $\bar{\mathcal{R}}_c = \frac{V_c}{A_c}$.

Let us note that this new formulation of the area-weighted discretization relies on a node-centered solver which is exactly the same as the one developed in [20] for two-dimensional Cartesian geometry. However, the present spatial discretization does not satisfy rigorously the GCL compatibility requirement. In what follows, we will assess the discrepancy of our discretization to the GCL by analyzing the corresponding discrete divergence operator. The discrete divergence operator that corresponds to the present scheme writes as

$$ (\nabla \cdot \mathbf{U})_c = \frac{1}{V_c} \sum_{p \in \mathcal{P}(c)} \mathcal{R}_p (L_{pc}^{-} N_{pc}^{-} + L_{pc}^{+} N_{pc}^{+}) \cdot \mathbf{U}_p, \quad (20) $$

where $\mathcal{R}_p$ denotes the pseudo-radius of vertex $p$. It is shown in [22, 19] that the discrete divergence operator deduced from the GCL reads

$$ (\nabla \cdot \mathbf{U})_{c}^{\text{GCL}} = \frac{1}{V_c} \sum_{p \in \mathcal{P}(c)} \frac{1}{3} [(2\mathcal{R}_p + \mathcal{R}_p^{-}) L_{pc}^{-} N_{pc}^{-} + (2\mathcal{R}_p + \mathcal{R}_p^{+}) L_{pc}^{+} N_{pc}^{+}] \cdot \mathbf{U}_p. \quad (21) $$
If the time evolution of the position vector, $X_p$, of vertex $p$ is governed by the trajectory equation [14], then one can prove that the time rate of change of the cell volume, $V_c$, satisfies

$$ \frac{1}{V_c} \frac{dV_c}{dt} = (\nabla \cdot U)_c^{GCL}. $$

Subtracting (20) and (21) leads to

$$ (\nabla \cdot U)_c - (\nabla \cdot U)_c^{GCL} = \frac{1}{3V_c} \sum_{p \in \mathcal{P}(c)} [(R_p - R_p^-)L_{pc}^- N_{pc}^- + (R_p - R_p^+) L_{pc}^+ N_{pc}^+] \cdot U_p. \quad (22) $$

Knowing that the summation in the previous equation is cyclic, shifting the index in the second term of the right hand-side yields

$$ (\nabla \cdot U)_c - (\nabla \cdot U)_c^{GCL} = \frac{1}{3V_c} \sum_{p \in \mathcal{P}(c)} [(R_p^+ - R_p)L_{pc}^+ N_{pc}^+] \cdot (U_{p^+} - U_p). \quad (23) $$

In case of a one-dimensional spherical flow on an equi-angular polar grid, the right-hand side of the previous equation is equal to zero. To prove this result, let us consider a quadrangular cell of an equi-angular polar grid. The proof proceeds in the following two steps:

* Either $p$ and $p^+$ are located on the same angular sector and thus the nodal velocity $U_p$ and $U_{p^+}$ are colinear to the direction of the angular sector which is orthogonal to the unit outward normal $N_{pc}^+$. Hence, $(U_{p^+} - U_p) \cdot N_{pc}^+ = 0$.

* Or $p$ and $p^+$ are located on the same cercle of radius $R$, then the Cartesian components of their nodal velocities reads as

$$ U_p = U(R) \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix}, \quad U_{p^+} = U(R) \begin{pmatrix} \cos(\theta + \Delta\theta) \\ \sin(\theta + \Delta\theta) \end{pmatrix}. $$
If the time evolution of the position vector, $X_p$, of vertex $p$ is governed by the trajectory equation [14], then one can prove that the time rate of change of the cell volume, $V_c$, satisfies

$$ \frac{1}{V_c} \frac{dV_c}{dt} = (\nabla \cdot U)_c^{GCL}. $$

Subtracting (20) and (21) leads to

$$ (\nabla \cdot U)_c - (\nabla \cdot U)_c^{GCL} = \frac{1}{3V_c} \sum_{p \in \mathcal{P}(c)} [(R_p - R_p^-) L_{pc}^- N_{pc}^- + (R_p - R_p^+) L_{pc}^+ N_{pc}^+] \cdot U_p. \quad (22) $$

Knowing that the summation in the previous equation is cyclic, shifting the index in the second term of the right hand-side yields

$$ (\nabla \cdot U)_c - (\nabla \cdot U)_c^{GCL} = \frac{1}{3V_c} \sum_{p \in \mathcal{P}(c)} [(R_p^+ - R_p) L_{pc}^+ N_{pc}^+] \cdot (U_{p^+} - U_p). \quad (23) $$

In case of a one-dimensional spherical flow on an equi-angular polar grid, the right-hand side of the previous equation is equal to zero. To prove this result, let us consider a quadrangular cell of an equi-angular polar grid. The proof proceeds in the following two steps:

* Either $p$ and $p^+$ are located on the same angular sector and thus the nodal velocity $U_p$ and $U_{p^+}$ are colinear to the direction of the angular sector which is orthogonal to the unit outward normal $N_{pc}^+$. Hence, $(U_{p^+} - U_p) \cdot N_{pc}^+ = 0$.

* Or $p$ and $p^+$ are located on the same cercle of radius $R$, then the Cartesian components of their nodal velocities reads as

$$ U_p = U(R) \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix}, \quad U_{p^+} = U(R) \begin{pmatrix} \cos(\theta + \Delta\theta) \\ \sin(\theta + \Delta\theta) \end{pmatrix}. $$
Here, $\theta$ denotes the angle of the angular sector, $U(r)$ is the module of the one-dimensional velocity field, and $\Delta\theta$ is size of the angular sector. A straightforward computation shows that

$$ \mathbf{U}_{p^+} - \mathbf{U}_p = 2U(R) \sin\left(\frac{\Delta\theta}{2}\right) \begin{pmatrix} -\sin\left(\theta + \frac{\Delta\theta}{2}\right) \\ \cos\left(\theta + \frac{\Delta\theta}{2}\right) \end{pmatrix}. $$

Knowing that the unit outward normal is given by

$$ \mathbf{N}_{pc}^{+} = \begin{pmatrix} \cos\left(\theta + \frac{\Delta\theta}{2}\right) \\ \sin\left(\theta + \frac{\Delta\theta}{2}\right) \end{pmatrix}, $$

we obtain that $(\mathbf{U}_{p^+} - \mathbf{U}_p) \cdot \mathbf{N}_{pc}^{+} = 0$.

This ends the proof. This result shows that our new area-weighted discretization satisfies rigorously the GCL compatibility requirement for one-dimensional spherical flows on equi-angular polar grids.

### 3. MOF multi-material interface reconstruction phase in axisymmetric geometry

The method used in this work to reconstruct interfaces, is the MOF approach well adapted for treating multi-materials interface problems [1, 8]. Indeed, such a method enables to capture more accurately interfaces than the classical VOF strategy and allows the treatment of general multi-material flows (more than two materials) [10, 17]. This method has been recently extended in cylindrical geometries, for a single interface problem [2]. Here, extension to multi-material interface reconstruction phase to cylindrical coordinates is considered.
Here, $\theta$ denotes the angle of the angular sector, $U(r)$ is the module of the one-dimensional velocity field, and $\Delta\theta$ is size of the angular sector. A straightforward computation shows that

$$ \mathbf{U}_{p^+} - \mathbf{U}_p = 2U(R) \sin\left(\frac{\Delta\theta}{2}\right) \begin{pmatrix} -\sin\left(\theta + \frac{\Delta\theta}{2}\right) \\ \cos\left(\theta + \frac{\Delta\theta}{2}\right) \end{pmatrix}. $$

Knowing that the unit outward normal is given by

$$ \mathbf{N}_{pc}^{+} = \begin{pmatrix} \cos\left(\theta + \frac{\Delta\theta}{2}\right) \\ \sin\left(\theta + \frac{\Delta\theta}{2}\right) \end{pmatrix}, $$

we obtain that $(\mathbf{U}_{p^+} - \mathbf{U}_p) \cdot \mathbf{N}_{pc}^{+} = 0$.

This ends the proof. This result shows that our new area-weighted discretization satisfies rigorously the GCL compatibility requirement for one-dimensional spherical flows on equi-angular polar grids.

### 3. MOF multi-material interface reconstruction phase in axisymmetric geometry

The method used in this work to reconstruct interfaces, is the MOF approach well adapted for treating multi-materials interface problems [1, 8]. Indeed, such a method enables to capture more accurately interfaces than the classical VOF strategy and allows the treatment of general multi-material flows (more than two materials) [10, 17]. This method has been recently extended in cylindrical geometries, for a single interface problem [2]. Here, extension to multi-material interface reconstruction phase to cylindrical coordinates is considered.
### 3.1. Moment of fluid method

The main idea of MOF is to track each fluid in a cell using the zeroth and first moments [8]. Given these two moments, interface is linearly reconstructed insuring volume conservation. To this end, interface update is done minimizing the discrepancy between the given moments and the reconstructed moments of the polygon behind the interface. One should note that no information from neighboring cells is required. This method is exact for linear interfaces and is second order accurate for smoothly curved ones. In the context of multi-material configurations, one has to face to material ordering when reconstructing interface. The method presented here, allows to automatically determine the order of materials by constructing all the possible combination and choosing the sequence that leads to the configuration where the reconstructed moments are the closest to the given ones. The main difference between cylindrical and planar geometry relies in the definition of the different moments. Since the interface reconstruction is done under volume conservative assumption, the zeroth moment $M_{k,c}^0$ of the k-th fluid in each cell c is obviously given by

$$M_{k,c}^0 = \int_{\Omega_{k,c}} \mathcal{R}dA, \qquad (24)$$

from this moment we can deduce the volume fraction

$$\alpha_{k,c} = \frac{M_{k,c}^0}{V_c}, \qquad (25)$$

with the cell volume $V_c = \int_{\Omega_c} \mathcal{R}dA$.

Contrary to the zeroth moment, the first moment can be defined without any specific requirement. Thus, it is possible to compute them in the two
### 3.1. Moment of fluid method

The main idea of MOF is to track each fluid in a cell using the zeroth and first moments [8]. Given these two moments, interface is linearly reconstructed insuring volume conservation. To this end, interface update is done minimizing the discrepancy between the given moments and the reconstructed moments of the polygon behind the interface. One should note that no information from neighboring cells is required. This method is exact for linear interfaces and is second order accurate for smoothly curved ones. In the context of multi-material configurations, one has to face to material ordering when reconstructing interface. The method presented here, allows to automatically determine the order of materials by constructing all the possible combination and choosing the sequence that leads to the configuration where the reconstructed moments are the closest to the given ones. The main difference between cylindrical and planar geometry relies in the definition of the different moments. Since the interface reconstruction is done under volume conservative assumption, the zeroth moment $M_{k,c}^0$ of the k-th fluid in each cell c is obviously given by

$$M_{k,c}^0 = \int_{\Omega_{k,c}} \mathcal{R}dA, \qquad (24)$$

from this moment we can deduce the volume fraction

$$\alpha_{k,c} = \frac{M_{k,c}^0}{V_c}, \qquad (25)$$

with the cell volume $V_c = \int_{\Omega_c} \mathcal{R}dA$.

Contrary to the zeroth moment, the first moment can be defined without any specific requirement. Thus, it is possible to compute them in the two
following different manners. In the one hand we can use the natural extension to axisymmetric geometries

$$M_{k,c}^{1} = \int_{\Omega_{k,c}} \mathcal{R}XdA, \qquad (26)$$

and from this moment we deduce the pseudo-centroid

$$\mathbf{X}_{k,c} = \frac{M_{k,c}^{1}}{V_{k,c}}, \quad \text{with} \quad V_{k,c} = V_c \alpha_{k,c}. \qquad (27)$$

This pseudo-centroid for a matter of simplicity will be called here the axisymmetric centroid.

On the other hand it can also be done with a planar definition as follows

$$M_{k,c}^{1,pl} = \int_{\Omega_{k,c}} \mathbf{X}dA, \qquad (28)$$

and thus planar centroid will be obtained from

$$\mathbf{X}_{k,c}^{pl} = \frac{M_{k,c}^{1,pl}}{A_{k,c}}, \qquad (29)$$

where $A_{k,c}$ is the area of the k-th fluid in the cell c.

Since this interface reconstruction method is coupled to our Lagrangian hydrodynamics scheme it requires to update the volume fractions and material centroids. Using the equal strain assumption, the volume fractions do not evolve during the Lagrangian step (see [11] for more details). However, the centroid locations are given from the Lagrangian step using a barycentric combination of the new positions of the mesh nodes as done in [10].

## 3.2. Numerical validation

The main goal of this section is to compare the results given by both axisymmetric and planar formulations of the centroids on several static test
following different manners. In the one hand we can use the natural extension to axisymmetric geometries

$$M_{k,c}^{1} = \int_{\Omega_{k,c}} \mathcal{R}\mathbf{X}dA, \qquad (26)$$

and from this moment we deduce the pseudo-centroid

$$\mathbf{X}_{k,c} = \frac{M_{k,c}^{1}}{V_{k,c}}, \quad \text{with } V_{k,c} = V_c \alpha_{k,c}. \qquad (27)$$

This pseudo-centroid for a matter of simplicity will be called here the axisymmetric centroid.

On the other hand it can also be done with a planar definition as follows

$$M_{k,c}^{1,pl} = \int_{\Omega_{k,c}} \mathbf{X}dA, \qquad (28)$$

and thus planar centroid will be obtained from

$$\mathbf{X}_{k,c}^{pl} = \frac{M_{k,c}^{1,pl}}{A_{k,c}}, \qquad (29)$$

where $A_{k,c}$ is the area of the k-th fluid in the cell c.

Since this interface reconstruction method is coupled to our Lagrangian hydrodynamics scheme it requires to update the volume fractions and material centroids. Using the equal strain assumption, the volume fractions do not evolve during the Lagrangian step (see [11] for more details). However, the centroid locations are given from the Lagrangian step using a barycentric combination of the new positions of the mesh nodes as done in [10].

## 3.2. Numerical validation

The main goal of this section is to compare the results given by both axisymmetric and planar formulations of the centroids on several static test
cases in one cell. As in [8], we consider three different mixed-cell layouts
that are filament (without junction), T-junction and Y-junction. The first
two configurations correspond to $C^2$-serial partitions whereas the third is not.
In the considered test cases, the parameter $\chi$ corresponds to the radius of
the circles defining the interfaces. Two values are considered with $\chi = 1$ and
$\chi = 64$. In addition, the computation domain is reduced to the cell $[0; 1] \times [0; 1]$
(see figures Fig.4 and Fig.5).

In the first case, with $\chi = 1$, we notice small differences for the filament
case, no notable difference on the T-Junction but the Y-junction results
for axisymmetric and planar formulations present distinct interface positions
due to a different ordering of the materials. For a large radius $\chi = 64$, the
curves are reduced to piecewise linear interfaces. Then, the result using both
formulations are very close to each other. For the two first cases filament
and T-junction, the results are exact. Regarding the Y-junction, it remains
a good approximation. These results illustrate the capability of both planar
and axisymmetric centroid formulation for MOF to treat accurately multi-
material problem. Nevertheless, for consistency with the global cylindrical
coordinate formulation, the axisymmetric formulation for the centroids is
retained in the sequel.
cases in one cell. As in [8], we consider three different mixed-cell layouts
that are filament (without junction), T-junction and Y-junction. The first
two configurations correspond to $C^2$-serial partitions whereas the third is not.
In the considered test cases, the parameter $\chi$ corresponds to the radius of
the circles defining the interfaces. Two values are considered with $\chi = 1$ and
$\chi = 64$. In addition, the computation domain is reduced to the cell $[0; 1] \times [0; 1]$
(see figures Fig.4 and Fig.5).

In the first case, with $\chi = 1$, we notice small differences for the filament
case, no notable difference on the T-Junction but the Y-junction results
for axisymmetric and planar formulations present distinct interface positions
due to a different ordering of the materials. For a large radius $\chi = 64$, the
curves are reduced to piecewise linear interfaces. Then, the result using both
formulations are very close to each other. For the two first cases filament
and T-junction, the results are exact. Regarding the Y-junction, it remains
a good approximation. These results illustrate the capability of both planar
and axisymmetric centroid formulation for MOF to treat accurately multi-
material problem. Nevertheless, for consistency with the global cylindrical
coordinate formulation, the axisymmetric formulation for the centroids is
retained in the sequel.
<figure><img src="image_6.png" /><figcaption>Figure 4: MOF interface reconstruction test for three materials. From the top to the bottom: the true partitions for $\chi = 1$ and their MOF reconstructions obtained with planar and axisymmetric centroids. From the left to the right: filament, T-junction and Y-junction configurations.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: MOF interface reconstruction test for three materials. From the top to the bottom: the true partitions for $\chi = 1$ and their MOF reconstructions obtained with planar and axisymmetric centroids. From the left to the right: filament, T-junction and Y-junction configurations.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: MOF interface reconstruction test for three materials. From the top to the bottom: the true partitions for $\chi = 64$ and their MOF reconstructions obtained with planar and axisymmetric centroids. From the left to the right: filament, T-junction and Y-junction configurations.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: MOF interface reconstruction test for three materials. From the top to the bottom: the true partitions for $\chi = 64$ and their MOF reconstructions obtained with planar and axisymmetric centroids. From the left to the right: filament, T-junction and Y-junction configurations.</figcaption></figure>
4. Rezoning phase improvement for polar meshes

The rezoning phase introduced in [10, 11] consists in moving the La-
grangian grid to improve its geometric quality. The objective of this part is
to extend this approach to polar meshes. To this end, the proposed proce-
dure relies on two main steps. The first phase is dedicated to compute the
smoothed grid from the Lagrangian one through CNS method. Then the
final mesh is deduced from the smoothed one by a relaxation procedure to
keep the rezoned grid as close as possible to the Lagrangian grid in order
to insure computation accuracy and avoid unphysical mesh rezoning. In the
sequel one should note that rezoning is formulated only for planar geometry
in the frame {$0, X, Y$}.

For the sake of readability, in the rest of the paper the quantities without
any accent $a$ are associated to Lagrangian mesh. After the rezoning step
we use $a^{rez}$, and finally after relaxation the quantities related to the rezoned
mesh are noted with the tilde accent $\tilde{a}$.

4.1. *General condition number smoothing (GCNS)*

As it is pointed out in the introduction, CNS approach is well adapted
to rezone Cartesian meshes but it still suffers from drawbacks for polar ones.
Indeed, in this case the mesh seems to collapse (like an implosion) to the ori-
gin. To circumvent this difficulty, it has been proposed to modify the CNS
algorithm using specific weight associated to the mesh geometry [25] that
controls mesh rezoning with regards to the radius for example. Nevertheless,
this approach is not completely satisfactory. First, it strongly depends on
4. Rezoning phase improvement for polar meshes

The rezoning phase introduced in [10, 11] consists in moving the La-
grangian grid to improve its geometric quality. The objective of this part is
to extend this approach to polar meshes. To this end, the proposed proce-
dure relies on two main steps. The first phase is dedicated to compute the
smoothed grid from the Lagrangian one through CNS method. Then the
final mesh is deduced from the smoothed one by a relaxation procedure to
keep the rezoned grid as close as possible to the Lagrangian grid in order
to insure computation accuracy and avoid unphysical mesh rezoning. In the
sequel one should note that rezoning is formulated only for planar geometry
in the frame {$0, X, Y$}.

For the sake of readability, in the rest of the paper the quantities without
any accent $a$ are associated to Lagrangian mesh. After the rezoning step
we use $a^{rez}$, and finally after relaxation the quantities related to the rezoned
mesh are noted with the tilde accent $\tilde{a}$.

4.1. *General condition number smoothing (GCNS)*

As it is pointed out in the introduction, CNS approach is well adapted
to rezone Cartesian meshes but it still suffers from drawbacks for polar ones.
Indeed, in this case the mesh seems to collapse (like an implosion) to the ori-
gin. To circumvent this difficulty, it has been proposed to modify the CNS
algorithm using specific weight associated to the mesh geometry [25] that
controls mesh rezoning with regards to the radius for example. Nevertheless,
this approach is not completely satisfactory. First, it strongly depends on
the choice of the weight, that may affect the quality of the mesh which can be shifted in the opposite direction to the origin for example. Furthermore, there is still a residual compression near the origin due to singularity at this point. In conclusion, it does not preserve a uniform polar mesh. For this reason, a different strategy is presented here. The main idea developed here is to apply the CNS rezoning algorithm in (r, θ)-coordinate system. In fact, a polar mesh initially expressed using a Cartesian coordinates (X, Y) leads to a structured Cartesian mesh in (r, θ)-coordinates. Here, a general presentation of the algorithm is made for unstructured meshes.

Assuming that the resulting mesh from the Lagrangian phase is unfolded (otherwise untangling procedure is used to correct invalid cells [26]). Thus, the proposed algorithm consists for polar meshes in three different steps as depicted on FIG.17. For the sake of simplicity, we consider in the sequel only the case of Cartesian and polar structured meshes.

<figure><img src="image_8.png" /><figcaption>Figure 6: Rule representation for GCNS algorithm for a polar mesh.</figcaption></figure>

The first step, is dedicated to the mapping between Cartesian and polar coordinates. To this end, consider $c$ a given cell of the Lagrangian grid for $(X, Y)$-coordinates, $p \in \mathcal{P}(c)$ a node of this cell. Notation used in the sequel are depicted on FIG.7. The mapping between a point $p \in c$ of Cartesian
the choice of the weight, that may affect the quality of the mesh which can
be shifted in the opposite direction to the origin for example. Furthermore,
there is still a residual compression near the origin due to singularity at this
point. In conclusion, it does not preserve a uniform polar mesh. For this
reason, a different strategy is presented here. The main idea developed here
is to apply the CNS rezoning algorithm in (r, θ)-coordinate system. In fact, a
polar mesh initially expressed using a Cartesian coordinates (X, Y) leads to a
structured Cartesian mesh in (r, θ)-coordinates. Here, a general presentation
of the algorithm is made for unstructured meshes.

Assuming that the resulting mesh from the Lagrangian phase is unfolded
(otherwise untangling procedure is used to correct invalid cells [26]). Thus,
the proposed algorithm consists for polar meshes in three different steps as
depicted on FIG.17. For the sake of simplicity, we consider in the sequel only
the case of Cartesian and polar structured meshes.

<figure><img src="image_7.png" /><figcaption>Figure 6: Rule representation for GCNS algorithm for a polar mesh.</figcaption></figure>

The first step, is dedicated to the mapping between Cartesian and polar
coordinates. To this end, consider *c* a given cell of the Lagrangian grid for
(X,Y)-coordinates, *p* ∈ *P*(*c*) a node of this cell. Notation used in the sequel
are depicted on FIG.7. The mapping between a point *p* ∈ *c* of Cartesian
<figure><img src="image_9.png" /><figcaption>Figure 7: Notations and mapping between cartesian and polar coordinates.</figcaption></figure>

coordinates $\mathbf{X}_p = (X, Y)^t$ to $\hat{\mathbf{X}}_p = (r_p, \theta_p)^t$ in polar ones is done through the following linear transformation

$$ \hat{\mathbf{X}}_p = T\mathbf{X}_p, \qquad (30) $$

where $T(\mathbf{X}_p) \stackrel{def}{=} I_d - \beta I_d + \beta R(\mathbf{X}_p)$ with $\beta \in \{0, 1\}$. For $\beta = 1$ then $T(\mathbf{X}_p) = R(\mathbf{X}_p)$ with the rotation matrix

$$ R(\mathbf{X}_p) = \begin{pmatrix} \cos(\theta_p) & -\sin(\theta_p) \\ \sin(\theta_p) & \cos(\theta_p) \end{pmatrix}, $$

using the definition $\theta_p = \arctan\left(\frac{Y_p}{X_p}\right)$ and $r_p = \sqrt{X_p^2 + Y_p^2}$. In the case $\beta = 0$, this formula leads to Cartesian rezoning with the identical transformation $\hat{\mathbf{X}}_p = \mathbf{X}_p$. When mapping $(X, Y)$ to $(r, \theta)$, the origin node has to be specifically treated. Indeed the transformation (30) is not defined for this point. Then as it is needed in the rezoning algorithm in the $(r, \theta)$ frame, the origin node is defined by a mapping of the first row on $r = 0$ axis (see FIG.7). Note that these nodes are not used for the final backward mapping.

The second step is the GCNS algorithm. It is based on a minimization problem of a local functional that controls the quality of the mesh. As done
<figure><img src="image_8.png" /><figcaption>Figure 7: Notations and mapping between cartesian and polar coordinates.</figcaption></figure>

coordinates $\mathbf{X}_p = (X, Y)^t$ to $\hat{\mathbf{X}}_p = (r_p, \theta_p)^t$ in polar ones is done through the following linear transformation

$$ \hat{\mathbf{X}}_p = \mathbf{T} \mathbf{X}_p, \qquad (30) $$

where $\mathbf{T}(\mathbf{X}_p) \stackrel{def}{=} I_d - \beta I_d + \beta \mathbf{R}(\mathbf{X}_p)$ with $\beta \in \{0, 1\}$. For $\beta = 1$ then $\mathbf{T}(\mathbf{X}_p) = \mathbf{R}(\mathbf{X}_p)$ with the rotation matrix

$$ \mathbf{R}(\mathbf{X}_p) = \begin{pmatrix} \cos(\theta_p) & -\sin(\theta_p) \\ \sin(\theta_p) & \cos(\theta_p) \end{pmatrix}, $$

using the definition $\theta_p = \arctan\left(\frac{Y_p}{X_p}\right)$ and $r_p = \sqrt{X_p^2 + Y_p^2}$. In the case $\beta = 0$, this formula leads to Cartesian rezoning with the identical transformation $\hat{\mathbf{X}}_p = \mathbf{X}_p$. When mapping $(X, Y)$ to $(r, \theta)$, the origin node has to be specifically treated. Indeed the transformation (30) is not defined for this point. Then as it is needed in the rezoning algorithm in the $(r, \theta)$ frame, the origin node is defined by a mapping of the first row on $r = 0$ axis (see FIG.7). Note that these nodes are not used for the final backward mapping.

The second step is the GCNS algorithm. It is based on a minimization problem of a local functional that controls the quality of the mesh. As done
in [10, 11], one has to distinguish boundary nodes and internal node for which the smoothing procedure is different.

For internal nodes, let us introduce as in [15] the condition number for $(r, \theta)$-coordinates that writes

$$ \kappa(\hat{J}_{cp}) = \frac{\|\hat{\mathbf{X}}_{pp^{+}}\|^2 + \|\hat{\mathbf{X}}_{pp^{-}}\|^2}{\hat{A}_{cp}}, \qquad (31) $$

where $\hat{\mathbf{X}}_{pp^{\pm}} = \hat{\mathbf{X}}_p - \hat{\mathbf{X}}_{p^{\pm}}$, and $\hat{A}_{cp} = \det(\hat{J}_{cp})$ is the area of the triangle delimited by $\{p, p^{+}, p^{-}\}$ in the rezoned grid and $\hat{J}_{cp} = [\hat{\mathbf{X}}_{pp^{+}}| - \hat{\mathbf{X}}_{pp^{-}}]$ the 2 × 2 Jacobian matrix associated to each corner at vertex $p$ of cell $c$. Thanks to this condition number we define the local function associated to the node $p$

$$ F_p(\hat{\mathbf{X}}_p) = \sum_{c \in C(p)} \kappa(\hat{J}_{cp}), \qquad (32) $$

Finally, the new position $\hat{\mathbf{X}}_p^{rez}$ is obtained by the minimization of the local function $F_p$ using the first step of a Newton algorithm. This leads to the formula

$$ \hat{\mathbf{X}}_p^{rez} = \hat{\mathbf{X}}_p - \mathbf{H}_{cp}^{-1}(\hat{\mathbf{X}}_p)\nabla F_p(\hat{\mathbf{X}}_p), \qquad (33) $$

where $\mathbf{H}_{cp}^{-1}$ and $\nabla F_p$ are respectively the Cartesian 2 × 2 Hessian matrix and gradient related to the local functional $F_p$.

For boundary nodes, the rezoned position $\hat{\mathbf{X}}_p^{rez}$ of $p$ is computed in consistent way with the GCNS algorithm. To this end, $\hat{\mathbf{X}}_p^{rez}$ is given thanks to a second-order interpolation Bézier curve [11] leading to

$$ \hat{\mathbf{X}}_p^{rez} = \hat{\mathcal{X}}_p(s^{rez}) = (1 - (s^{rez})^2)\hat{\mathbf{X}}_{p^{-}} + 2(1 - s^{rez})\hat{\mathbf{X}}_i + (s^{rez})^2\hat{\mathbf{X}}_{p^{+}}, \qquad (34) $$

where $\hat{\mathbf{X}}_i$ such that $\hat{\mathcal{X}}_p(1/2) = \hat{\mathbf{X}}_p$. Furthermore, the parameter $s^{rez}$ is computed to minimize $F_p(\hat{\mathcal{X}}_p(s))$ (for more details on this procedure see [11]).
in [10, 11], one has to distinguish boundary nodes and internal node for which the smoothing procedure is different.

For internal nodes, let us introduce as in [15] the condition number for $(r, \theta)$-coordinates that writes

$$ \kappa(\hat{J}_{cp}) = \frac{\|\hat{\mathbf{X}}_{pp^{+}}\|^2 + \|\hat{\mathbf{X}}_{pp^{-}}\|^2}{\hat{A}_{cp}}, \qquad (31) $$

where $\hat{\mathbf{X}}_{pp^{\pm}} = \hat{\mathbf{X}}_p - \hat{\mathbf{X}}_{p^{\pm}}$, and $\hat{A}_{cp} = \det(\hat{J}_{cp})$ is the area of the triangle delimited by $\{p, p^{+}, p^{-}\}$ in the rezoned grid and $\hat{J}_{cp} = [\hat{\mathbf{X}}_{pp^{+}}| - \hat{\mathbf{X}}_{pp^{-}}]$ the 2 × 2 Jacobian matrix associated to each corner at vertex $p$ of cell $c$. Thanks to this condition number we define the local function associated to the node $p$

$$ F_p(\hat{\mathbf{X}}_p) = \sum_{c \in \mathcal{C}(p)} \kappa(\hat{J}_{cp}), \qquad (32) $$

Finally, the new position $\hat{\mathbf{X}}_p^{rez}$ is obtained by the minimization of the local function $F_p$ using the first step of a Newton algorithm. This leads to the formula

$$ \hat{\mathbf{X}}_p^{rez} = \hat{\mathbf{X}}_p - \mathbf{H}_{cp}^{-1}(\hat{\mathbf{X}}_p)\nabla F_p(\hat{\mathbf{X}}_p), \qquad (33) $$

where $\mathbf{H}_{cp}^{-1}$ and $\nabla F_p$ are respectively the Cartesian 2 × 2 Hessian matrix and gradient related to the local functional $F_p$.

For boundary nodes, the rezoned position $\hat{\mathbf{X}}_p^{rez}$ of $p$ is computed in consistent way with the GCNS algorithm. To this end, $\hat{\mathbf{X}}_p^{rez}$ is given thanks to a second-order interpolation Bézier curve [11] leading to

$$ \hat{\mathbf{X}}_p^{rez} = \hat{\mathcal{X}}_p(s^{rez}) = (1 - (s^{rez})^2)\hat{\mathbf{X}}_{p^{-}} + 2(1 - s^{rez})\hat{\mathbf{X}}_i + (s^{rez})^2\hat{\mathbf{X}}_{p^{+}}, \qquad (34) $$

where $\hat{\mathbf{X}}_i$ such that $\hat{\mathcal{X}}_p(1/2) = \hat{\mathbf{X}}_p$. Furthermore, the parameter $s^{rez}$ is computed to minimize $F_p(\hat{\mathcal{X}}_p(s))$ (for more details on this procedure see [11]).
Finally, the third step consists in backward mapping between $\hat{\mathbf{X}}_{p}^{rez}$ and $\mathbf{X}_{p}^{rez}$ using (30), where the inverse of the transformation matrix is taken equal to $\mathbf{T}^{-1}(\mathbf{X}_{p}) \stackrel{def}{=} I_{d} - \beta I_{d} + \beta [\mathbf{R}(\mathbf{X}_{p})]^{-1}$ with $\beta \in \{0, 1\}$ and $[\mathbf{R}(\mathbf{X}_{p})]^{-1}$ the inverse rotation matrix.

## 4.2. Relaxation algorithm

The relaxation algorithm consists in making a convex combination between rezoned grid obtained from GCNS step and its location after Lagrangian step. This reads for each mesh node $p$ by:

$$ \tilde{\mathbf{X}}_{p} = \mathbf{X}_{p} + \omega_{p}(\mathbf{X}_{p}^{rez} - \mathbf{X}_{p}), \text{ with } \omega_{p} \in [0, 1], $$

where $\tilde{\mathbf{X}}_{p}$ is the new mesh node position after the complete rezoning phase. The coefficient $\omega_{p}$ is computed as a function of the right Cauchy-Green tensor associated to the Lagrange grid deformation over a time step (for details see [11, 18]).

## 4.3. Numerical validation

In this section, we compare results obtained by the GCNS algorithm to those obtained for classical CNS for the rezoning of uniform polar and unstructured meshes.

**Uniform mesh.** First, we consider an uniform polar mesh made of 20 × 10 elements see Fig.8-(a). Results obtained after 100 iterations for the classical and general smoothing are presented on Fig.8. For each method the relaxation coefficient $\omega_{p}$ is taken equal to 1. As already mentioned, the classical smoothing does not converge on polar mesh and implies the collapse of cell layers to the origins (see Fig.8-(b)). However, for the GCNS, the
Finally, the third step consists in backward mapping between $\hat{\mathbf{X}}_p^{rez}$ and $\mathbf{X}_p^{rez}$ using (30), where the inverse of the transformation matrix is taken equal to $\mathbf{T}^{-1}(\mathbf{X}_p) \stackrel{def}{=} I_d - \beta I_d + \beta [\mathbf{R}(\mathbf{X}_p)]^{-1}$ with $\beta \in \{0, 1\}$ and $[\mathbf{R}(\mathbf{X}_p)]^{-1}$ the inverse rotation matrix.

## 4.2. Relaxation algorithm

The relaxation algorithm consists in making a convex combination between rezoned grid obtained from GCNS step and its location after Lagrangian step. This reads for each mesh node $p$ by:

$$ \tilde{\mathbf{X}}_p = \mathbf{X}_p + \omega_p(\mathbf{X}_p^{rez} - \mathbf{X}_p), \text{ with } \omega_p \in [0, 1], $$

where $\tilde{\mathbf{X}}_p$ is the new mesh node position after the complete rezoning phase. The coefficient $\omega_p$ is computed as a function of the right Cauchy-Green tensor associated to the Lagrange grid deformation over a time step (for details see [11, 18]).

## 4.3. Numerical validation

In this section, we compare results obtained by the GCNS algorithm to those obtained for classical CNS for the rezoning of uniform polar and unstructured meshes.

**Uniform mesh.** First, we consider an uniform polar mesh made of 20 × 10 elements see Fig.8-(a). Results obtained after 100 iterations for the classical and general smoothing are presented on Fig.8. For each method the relaxation coefficient $\omega_p$ is taken equal to 1. As already mentioned, the classical smoothing does not converge on polar mesh and implies the collapse of cell layers to the origins (see Fig.8-(b)). However, for the GCNS, the
result obtained (see FIG.8-(c)) is converged. The mesh initially uniform, is not modified at the end of the computation. This clearly illustrates the good behavior of our smoothing algorithm.

<figure><img src="image_10.png" /><figcaption>Figure 8: Smoothing of a static polar grid 16 × 10: (a) initial grid; Smoothed grids after 100 iterations: (b) CNS, (c) GCNS.</figcaption></figure>

**Unstructured mesh.** Now, rezoning for an unstructured mesh is studied. Let us consider a mesh made of 175 quadrangular cells as depicted on FIG.9-(a). When applying the full Cartesian rezoning to the mesh, similar observations as previously can be made. It suffers from an implosion of central cells to the origin and does not converge (see FIG.9-(b)). For the full GCNS algorithm, one can see after convergence, the formation of mesh distortion on the square region and a polar mesh far from the center (see FIG.9-(c)). Nevertheless, it is possible to improve this rezoning. Thus, the main idea developed in the sequel is to apply the GCNS rezoning algorithm differently for a node belonging initially to a Cartesian or polar region of the mesh. To this end, the transformation $T$ between $(X, Y)$ and $(r, \theta)$ coordinates is
result obtained (see FIG.8-(c)) is converged. The mesh initially uniform, is not modified at the end of the computation. This clearly illustrates the good behavior of our smoothing algorithm.

<figure><img src="image_9.png" /><figcaption>Figure 8: Smoothing of a static polar grid 16 × 10: (a) initial grid; Smoothed grids after 100 iterations: (b) CNS, (c) GCNS.</figcaption></figure>

**Unstructured mesh.** Now, rezoning for an unstructured mesh is studied. Let us consider a mesh made of 175 quadrangular cells as depicted on FIG.9-(a). When applying the full Cartesian rezoning to the mesh, similar observations as previously can be made. It suffers from an implosion of central cells to the origin and does not converge (see FIG.9-(b)). For the full GCNS algorithm, one can see after convergence, the formation of mesh distortion on the square region and a polar mesh far from the center (see FIG.9-(c)). Nevertheless, it is possible to improve this rezoning. Thus, the main idea developed in the sequel is to apply the GCNS rezoning algorithm differently for a node belonging initially to a Cartesian or polar region of the mesh. To this end, the transformation $T$ between $(X, Y)$ and $(r, \theta)$ coordinates is
<figure><img src="image_11.png" /><figcaption>Figure 9: Smoothing of a static unstructured grid: (a) initial grid; Smoothed grids after 100 iterations (b) CNS, (c) GCNS after 100 iterations.</figcaption></figure>

modified in the following way

$$T(\mathbf{X}_p) \stackrel{\text{def}}{=} I_d - \beta(\mathbf{X}_p)I_d + \beta(\mathbf{X}_p)R(\mathbf{X}_p), \qquad (35)$$

with

$$\beta(\mathbf{X}_p) = \begin{cases} 1 & \text{if } \mathbf{X}_p^0 \in \mathcal{P}^{\text{pol}}, \\ 0 & \text{if } \mathbf{X}_p^0 \in \mathcal{P}^{\text{car}}, \end{cases}$$

where $\mathcal{P}^{\text{car}}$ and $\mathcal{P}^{\text{pol}}$ are the sets of nodes that belong to the Cartesian and respectively polar region of the mesh at the initial time. These regions are represented thanks to red and blue color (see Fig.10-(a)) for the considered mesh. Nodes localized at the frontier between the polar and Cartesian meshes (black nodes on Fig.10-(a)) can be considered either polar, or Cartesian. As represented on Fig.10-(b,c), both possibilities are tested. The obtained results illustrate that the Cartesian choice remains better contrary to the polar one that introduce mesh distortion.
<figure><img src="image_10.png" /><figcaption>Figure 9: Smoothing of a static unstructured grid: (a) initial grid; Smoothed grids after 100 iterations (b) CNS, (c) GCNS after 100 iterations.</figcaption></figure>

modified in the following way

$$T(\mathbf{X}_p) \stackrel{\text{def}}{=} I_d - \beta(\mathbf{X}_p)I_d + \beta(\mathbf{X}_p)R(\mathbf{X}_p), \quad (35)$$

with

$$\beta(\mathbf{X}_p) = \begin{cases} 1 & \text{if } \mathbf{X}_p^0 \in \mathcal{P}^{\text{pol}}, \\ 0 & \text{if } \mathbf{X}_p^0 \in \mathcal{P}^{\text{car}}, \end{cases}$$

where $\mathcal{P}^{\text{car}}$ and $\mathcal{P}^{\text{pol}}$ are the sets of nodes that belong to the Cartesian and respectively polar region of the mesh at the initial time. These regions are represented thanks to red and blue color (see Fig.10-(a)) for the considered mesh. Nodes localized at the frontier between the polar and Cartesian meshes (black nodes on Fig.10-(a)) can be considered either polar, or Cartesian. As represented on Fig.10-(b,c), both possibilities are tested. The obtained results illustrate that the Cartesian choice remains better contrary to the polar one that introduce mesh distortion.
<figure><img src="image_12.png" /><figcaption>Figure 10: Smoothing of a static unstructured grid: (a) initial grid with Cartesian (blue) and polar (red) rezoning regions; Smoothed grids after 100 iterations (b) GCNS with interfacial polar rezoning, (c) GCNS with interfacial Cartesian rezoning.</figcaption></figure>

5. Hybrid remapping in axisymmetric geometry

During the remapping phase, the physical unknowns (density, velocity, total energy) computed thanks to the Lagrangian step are conservatively remapped from the Lagrangian mesh to the rezoned one. To this end, an extension of the *Hybrid Remapping Algorithm* for multi-material flows [10, 16, 4] to cylindrical geometry is proposed here. This strategy consists in the following two steps. First a *swept-faced remapping* is used to treat cells and nodes localized far from the interface. Then, a *cell-intersection-based method* [11] is applied to the cells and nodes in the neighborhood of the interface. In this way, this approach combines the ability of the cell-intersection method to remap the interface and the efficiency of the swept flux approach for the other cells that significantly reduce the global computing cost of the method. As done previously, in the perspective of general use of the method, a global formulation including both Cartesian and axisymmetric framework is presented.
<figure><img src="image_11.png" /><figcaption>Figure 10: Smoothing of a static unstructured grid: (a) initial grid with Cartesian (blue) and polar (red) rezoning regions; Smoothed grids after 100 iterations (b) GCNS with interfacial polar rezoning, (c) GCNS with interfacial Cartesian rezoning.</figcaption></figure>

5. Hybrid remapping in axisymmetric geometry

During the remapping phase, the physical unknowns (density, velocity, total energy) computed thanks to the Lagrangian step are conservatively remapped from the Lagrangian mesh to the rezoned one. To this end, an extension of the *Hybrid Remapping Algorithm* for multi-material flows [10, 16, 4] to cylindrical geometry is proposed here. This strategy consists in the following two steps. First a *swept-faced remapping* is used to treat cells and nodes localized far from the interface. Then, a *cell-intersection-based method* [11] is applied to the cells and nodes in the neighborhood of the interface. In this way, this approach combines the ability of the cell-intersection method to remap the interface and the efficiency of the swept flux approach for the other cells that significantly reduce the global computing cost of the method. As done previously, in the perspective of general use of the method, a global formulation including both Cartesian and axisymmetric framework is presented.
We assume in the sequel, that there is no topology change of the mesh, the cells of the Lagrangian and rezoned grids are respectively designed by $\Omega_c$ and $\tilde{\Omega}_c$.

### 5.1. Multi-material cell-intersection-based (MCIB) remapping

<figure><img src="image_13.png" /><figcaption>Figure 11: Notations for MCIB method.</figcaption></figure>

The main goal of remapping is as follows. Given the piecewise constant representation of the physical variables per unit of volume $(\rho, \rho U, \rho E)$ noted $\psi_c = \rho_c \phi_c$ in each cell of the Lagrangian grid, we want to compute its equivalent $\tilde{\psi}_c$ in each cell of the rezoned grid given as

$$ \tilde{\psi}_c = \frac{1}{\tilde{V}_c} \int_{\tilde{\Omega}_c} \rho \tilde{\phi} R dA, \qquad (36) $$

with $\tilde{V}_c$ the volume of the cell $\tilde{\Omega}_c$. Contrary to single fluid approach, here the rezoned values $\tilde{\psi}_c$ can not be computed directly in each cell c. In fact, one has to take into account multi-material aspects.
We assume in the sequel, that there is no topology change of the mesh, the cells of the Lagrangian and rezoned grids are respectively designed by $\Omega_c$ and $\tilde{\Omega}_c$.

### 5.1. Multi-material cell-intersection-based (MCIB) remapping

<figure><img src="image_12.png" /><figcaption>Figure 11: Notations for MCIB method.</figcaption></figure>

The main goal of remapping is as follows. Given the piecewise constant representation of the physical variables per unit of volume $(\rho, \rho U, \rho E)$ noted $\psi_c = \rho_c \phi_c$ in each cell of the Lagrangian grid, we want to compute its equivalent $\tilde{\psi}_c$ in each cell of the rezoned grid given as

$$ \tilde{\psi}_c = \frac{1}{\tilde{V}_c} \int_{\tilde{\Omega}_c} \rho \tilde{\phi} R dA, \qquad (36) $$

with $\tilde{V}_c$ the volume of the cell $\tilde{\Omega}_c$. Contrary to single fluid approach, here the rezoned values $\tilde{\psi}_c$ can not be computed directly in each cell c. In fact, one has to take into account multi-material aspects.
First of all, let us introduce some notations. Each material of the flow noted $k$ occupies the polygon $\Omega_{k,c} \subset \Omega_c$, within the MOF framework, such that $\Omega_c = \bigcup_k \Omega_{k,c}$ and is characterized by its partial mass, density, pressure, internal energy and variables per unit of volume (total energy, momentum) whose averaged values in each sub-cell are respectively $m_{k,c}$, $\rho_{k,c}$, $P_{k,c}$, $\varepsilon_{k,c}$ and $\psi_{k,c} = \rho_{k,c}\phi_{k,c}$ with $\phi_{k,c}$ the partial velocity or energy per unit of mass.

Thus, for multi-material flow, the main idea of remapping is not to directly compute the global rezoned quantities $\tilde{\psi}_c$ but the partial rezoned ones noted $\tilde{\psi}_{k,c}$. This is particularly true for the MCIB method that is dedicated to treat cell in the interface neighborhood. To this end, we first propose a second order reconstruction $\Psi_{k,c}(\mathbf{X})$ of $\psi_{k,c}$ over each Lagrangian cell $c$ through the piecewise linear function

$$ \Psi_{k,c}(\mathbf{X}) = \psi_{k,c} + (\nabla \Psi_k)_c(\mathbf{X} - \mathbf{X}_{k,c}), \quad (37) $$

where $(\nabla \Psi_k)_c$ denotes the constant gradient of $\Psi_{k,c}$ within cell $c$ computed thanks to a least-squares approach. Finally $\mathbf{X}_{k,c}$ is the centroid related to the $k$-th fluid in the cell $c$ given by

$$ \mathbf{X}_{k,c} = \frac{1}{V_{k,c}} \int_{\Omega_{k,c}} \mathcal{R} \mathbf{X} dA. \quad (38) $$

Thanks to these notations, the remapped value for MCIB is given by

$$ \tilde{\psi}_{k,c} = \frac{1}{\tilde{V}_{k,c}} \sum_{d \in \mathcal{C}(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R} \Psi_{k,c} dA, \quad (39) $$

where the intersection polygons $\Omega_{k,d} \cap \tilde{\Omega}_c$ are computed thanks to a specific triangulation of the mesh. The procedure is detailed in [11]. The set $\mathcal{C}(c)$ contains the cells including $c$ that share at least one node with the cell $c$. At last,
First of all, let us introduce some notations. Each material of the flow noted $k$ occupies the polygon $\Omega_{k,c} \subset \Omega_c$, within the MOF framework, such that $\Omega_c = \bigcup_k \Omega_{k,c}$ and is characterized by its partial mass, density, pressure, internal energy and variables per unit of volume (total energy, momentum) whose averaged values in each sub-cell are respectively $m_{k,c}$, $\rho_{k,c}$, $P_{k,c}$, $\varepsilon_{k,c}$ and $\psi_{k,c} = \rho_{k,c}\phi_{k,c}$ with $\phi_{k,c}$ the partial velocity or energy per unit of mass.

Thus, for multi-material flow, the main idea of remapping is not to directly compute the global rezoned quantities $\tilde{\psi}_c$ but the partial rezoned ones noted $\tilde{\psi}_{k,c}$. This is particularly true for the MCIB method that is dedicated to treat cell in the interface neighborhood. To this end, we first propose a second order reconstruction $\Psi_{k,c}(\mathbf{X})$ of $\psi_{k,c}$ over each Lagrangian cell $c$ through the piecewise linear function

$$ \Psi_{k,c}(\mathbf{X}) = \psi_{k,c} + (\nabla \Psi_k)_c(\mathbf{X} - \mathbf{X}_{k,c}), \quad (37) $$

where $(\nabla \Psi_k)_c$ denotes the constant gradient of $\Psi_{k,c}$ within cell $c$ computed thanks to a least-squares approach. Finally $\mathbf{X}_{k,c}$ is the centroid related to the $k$-th fluid in the cell $c$ given by

$$ \mathbf{X}_{k,c} = \frac{1}{V_{k,c}} \int_{\Omega_{k,c}} \mathcal{R} \mathbf{X} dA. \quad (38) $$

Thanks to these notations, the remapped value for MCIB is given by

$$ \tilde{\psi}_{k,c} = \frac{1}{\tilde{V}_{k,c}} \sum_{d \in \mathcal{C}(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R} \Psi_{k,c} dA, \quad (39) $$

where the intersection polygons $\Omega_{k,d} \cap \tilde{\Omega}_c$ are computed thanks to a specific triangulation of the mesh. The procedure is detailed in [11]. The set $\mathcal{C}(c)$ contains the cells including $c$ that share at least one node with the cell $c$. At last,
the partial volume defined on the rezoned cell is $\tilde{V}_{k,c} = \sum_{d \in C(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R}dA$.

In the context of MOF reconstruction, one has to define additional quantities as the partial remapped mass corresponding to material k. It is computed as $\tilde{m}_{k,c} = \tilde{\rho}_{k,c} \tilde{\alpha}_{k,c} \tilde{V}_{k,c}$ with the volume fraction

$$ \tilde{\alpha}_{k,c} = \frac{1}{\tilde{V}_c} \sum_{d \in C(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R}dA, \quad (40) $$

thus the partial volume can be also expressed as $\tilde{V}_{k,c} = \tilde{V}_c \tilde{\alpha}_{k,c}$. In addition, each material centroid position is defined thanks to

$$ \tilde{\mathbf{X}}_{k,c} = \frac{1}{\tilde{V}_{k,c}} \sum_{d \in C(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R}\mathbf{X}dA. \quad (41) $$

## 5.2. Pure cell swept-face (PCSF) remapping

<figure><img src="image_14.png" /><figcaption>Figure 12: Notations for swept face-based method.</figcaption></figure>

As explained before, the PCSF remapping is used only to treat single fluid cells. In this context, one should remark that $\Omega_c = \Omega_{k,c}$, thus the mean
the partial volume defined on the rezoned cell is $\tilde{V}_{k,c} = \sum_{d \in \mathcal{C}(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R}dA$.

In the context of MOF reconstruction, one has to define additional quantities as the partial remapped mass corresponding to material k. It is computed as $\tilde{m}_{k,c} = \tilde{\rho}_{k,c} \tilde{\alpha}_{k,c} \tilde{V}_{k,c}$ with the volume fraction

$$ \tilde{\alpha}_{k,c} = \frac{1}{\tilde{V}_c} \sum_{d \in \mathcal{C}(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R}dA, \quad (40) $$

thus the partial volume can be also expressed as $\tilde{V}_{k,c} = \tilde{V}_c \tilde{\alpha}_{k,c}$. In addition, each material centroid position is defined thanks to

$$ \tilde{\mathbf{X}}_{k,c} = \frac{1}{\tilde{V}_{k,c}} \sum_{d \in \mathcal{C}(c)} \int_{\Omega_{k,d} \cap \tilde{\Omega}_c} \mathcal{R}\mathbf{X}dA. \quad (41) $$

## 5.2. Pure cell swept-face (PCSF) remapping

<figure><img src="image_13.png" /><figcaption>Figure 12: Notations for swept face-based method.</figcaption></figure>

As explained before, the PCSF remapping is used only to treat single fluid cells. In this context, one should remark that $\Omega_c = \Omega_{k,c}$, thus the mean
value $\tilde{\psi}_{c,k}$ is given through

$$ \tilde{\psi}_{k,c} = \psi_{k,c} + \sum_{f \in \mathcal{F}(c)} \int_{A_f} \mathcal{R} \Psi_{k,f} dA, \quad (42) $$

with $A_f$ the quadrangular signed area swept by the face $f$ of a cell $c$ between the Lagrangian grid and the rezoned grid delimited by the ordered nodes of coordinates $\{\mathbf{X}_p, \tilde{\mathbf{X}}_p, \tilde{\mathbf{X}}_{p+}, \mathbf{X}_{p+}\}$ (refer to fig FIG.12). We note $\mathcal{F}(c)$ the set of the faces $f$ of a cell $c$. In addition, $\Psi_{k,f}$ is the upwind value given by

$$ \Psi_{k,f} = \begin{cases} \Psi_{k,c^+} & \text{if } A_f > 0 \\ \Psi_{k,c} & \text{otherwise.} \end{cases} \quad (43) $$

with $c^+$ the neighbor cell of $c$ through the face $f$. During this step the volume fractions $\tilde{\alpha}_{k,c} = \alpha_{k,c}$ do not change as we consider single fluid cells and the material centroid can be updated directly from the geometry $\tilde{\mathbf{X}}_{k,c} = \tilde{\mathbf{X}}_c$ where $\tilde{\mathbf{X}}_c$ is the centroid of the cell $\tilde{\Omega}_c$.

### 5.3. Integration strategy

For both PCSF and MCIB remapping, one has to compute several surface integrals, on polygons where the integrand is a polynomial function of $(X, Y)$. This can be done using a triangulation of these areas. Nevertheless, this is expensive. Here, we rather adopt a more efficient method as in [23]. In this context, integrals are simplified using Taylor decomposition of the polynomial integrand and Green's formula leading to compute circular integrals over the edges of the polygons defining the integration areas. For further details on integral computations see [23].
value $\tilde{\psi}_{c,k}$ is given through

$$ \tilde{\psi}_{k,c} = \psi_{k,c} + \sum_{f \in \mathcal{F}(c)} \int_{A_f} \mathcal{R} \Psi_{k,f} dA, \qquad (42) $$

with $A_f$ the quadrangular signed area swept by the face $f$ of a cell $c$ between the Lagrangian grid and the rezoned grid delimited by the ordered nodes of coordinates $\{\mathbf{X}_p, \tilde{\mathbf{X}}_p, \tilde{\mathbf{X}}_{p+}, \mathbf{X}_{p+}\}$ (refer to fig FIG.12). We note $\mathcal{F}(c)$ the set of the faces $f$ of a cell $c$. In addition, $\Psi_{k,f}$ is the upwind value given by

$$ \Psi_{k,f} = \begin{cases} \Psi_{k,c^+} & \text{if } A_f > 0 \\ \Psi_{k,c} & \text{otherwise.} \end{cases} \qquad (43) $$

with $c^+$ the neighbor cell of $c$ through the face $f$. During this step the volume fractions $\tilde{\alpha}_{k,c} = \alpha_{k,c}$ do not change as we consider single fluid cells and the material centroid can be updated directly from the geometry $\tilde{\mathbf{X}}_{k,c} = \tilde{\mathbf{X}}_c$ where $\tilde{\mathbf{X}}_c$ is the centroid of the cell $\tilde{\Omega}_c$.

### 5.3. Integration strategy

For both PCSF and MCIB remapping, one has to compute several surface integrals, on polygons where the integrand is a polynomial function of $(X, Y)$. This can be done using a triangulation of these areas. Nevertheless, this is expensive. Here, we rather adopt a more efficient method as in [23]. In this context, integrals are simplified using Taylor decomposition of the polynomial integrand and Green's formula leading to compute circular integrals over the edges of the polygons defining the integration areas. For further details on integral computations see [23].
<figure><img src="image_15.png" /><figcaption>Figure 13: Hybrid remapping principle in one-dimension case.</figcaption></figure>

5.4. *Hybrid remapping algorithm*

In this part, we detail the hybrid remapping algorithm that is summarized
on FIG.13. To this end, let us introduce $\mathcal{N}^P$ and $\mathcal{N}^M$ the sets of nodes and
in the same manner $\mathcal{C}^P$ and $\mathcal{C}^M$ the sets of cells respectively used for PCSF
and MCIB remapping. Here $\mathcal{N}^M$ collects mixed nodes belonging to cells
that contain the interface or are on this interface (white nodes on FIG.13)
despite $\mathcal{N}^P$ contains the pure ones (black nodes on FIG.13). In addition, $\mathcal{C}^M$
is the set of mixed cell that include cells intersected by the interfaces and
their neighbors by nodes. Finally, $\mathcal{C}^P$ contains the cells that have at least
one node in $\mathcal{N}^P$.

The hybrid remapping procedure consists in performing the following steps.
<figure><img src="image_14.png" /><figcaption>Figure 13: Hybrid remapping principle in one-dimension case.</figcaption></figure>

5.4. *Hybrid remapping algorithm*

In this part, we detail the hybrid remapping algorithm that is summarized
on FIG.13. To this end, let us introduce $\mathcal{N}^P$ and $\mathcal{N}^M$ the sets of nodes and
in the same manner $\mathcal{C}^P$ and $\mathcal{C}^M$ the sets of cells respectively used for PCSF
and MCIB remapping. Here $\mathcal{N}^M$ collects mixed nodes belonging to cells
that contain the interface or are on this interface (white nodes on FIG.13)
despite $\mathcal{N}^P$ contains the pure ones (black nodes on FIG.13). In addition, $\mathcal{C}^M$
is the set of mixed cell that include cells intersected by the interfaces and
their neighbors by nodes. Finally, $\mathcal{C}^P$ contains the cells that have at least
one node in $\mathcal{N}^P$.

The hybrid remapping procedure consists in performing the following steps.
1. **PCSF step.** In this step we first move the pure nodes included in $\mathcal{N}^P$ and we remap the quantities in cell $c$ belonging to $\mathcal{C}^P$. Thus, we have $\tilde{\psi}_{k,c} = (\tilde{\rho}_{k,c}, \tilde{\rho} \tilde{E}_{k,c}, \tilde{\rho} \tilde{U}_{k,c})$ using relation (42) and $\tilde{m}_{k,c}, \tilde{\alpha}_{k,c}, \tilde{\mathbf{X}}_{k,c}$ for each cell $c \in \mathcal{C}^P$.

2. **MCIB step.** Now, the mixed nodes in $\mathcal{N}^M$ are moved and the $\tilde{\psi}_{k,c} = (\tilde{\rho}_{k,c}, \tilde{\rho} \tilde{E}_{k,c}, \tilde{\rho} \tilde{U}_{k,c})$ are remapped thanks to (39) and $\tilde{m}_{k,c}, \tilde{\alpha}_{k,c}, \tilde{\mathbf{X}}_{k,c}$ are computed for cells $c \in \mathcal{C}^M$.

Since $\mathcal{C}^M \cap \mathcal{C}^P \neq \{\emptyset\}$, one should note that cell included in this intersection are remapped at each step of the algorithm.

At the end of remapping, only the partial values of the physical variables per unit of volume are known. A this step, a first point is to compute the physical variables per unit of mass. The remapped partial total energy is given using $\tilde{E}_{k,c} = (\tilde{\rho} \tilde{E})_{k,c} / \tilde{\rho}_{k,c}$. However, this is different for the remapped partial velocity $\tilde{\mathbf{U}}_{k,c}$. Indeed, as explained in the second part of this paper, the Lagrangian computation of the velocity is done in Cartesian geometry. For this reason, the remapped velocity is deduced from the $(\tilde{\rho} \tilde{\mathbf{U}})_{k,c}$ through $\tilde{\mathbf{U}}_{k,c} = (\tilde{\rho} \tilde{\mathbf{U}})_{k,c}^{pl} / \tilde{\rho}_{k,c}^{pl}$ using the planar remapped density and momentum given through (39) and (42) with $\mathcal{R} = 1$. The second point is dedicated to the reconstruction of the global values required for the next Lagrangian step. To this end, a classical procedure is to use specific averages

$$ \tilde{\phi}_c = \frac{1}{\tilde{m}_c} \sum_k \tilde{m}_{k,c} \tilde{\phi}_{k,c}, \qquad (44) $$

with the global mass and density deduced from

$$ \tilde{m}_c = \sum_k \tilde{m}_{k,c} \tilde{\alpha}_{k,c} \quad \text{and} \quad \tilde{\rho}_c = \sum_k \tilde{\rho}_{k,c} \tilde{\alpha}_{k,c}. \qquad (45) $$
1. **PCSF step.** In this step we first move the pure nodes included in $\mathcal{N}^P$ and we remap the quantities in cell $c$ belonging to $\mathcal{C}^P$. Thus, we have $\tilde{\psi}_{k,c} = (\tilde{\rho}_{k,c}, \tilde{\rho} \tilde{E}_{k,c}, \tilde{\rho} \tilde{U}_{k,c})$ using relation (42) and $\tilde{m}_{k,c}, \tilde{\alpha}_{k,c}, \tilde{\mathbf{X}}_{k,c}$ for each cell $c \in \mathcal{C}^P$.

2. **MCIB step.** Now, the mixed nodes in $\mathcal{N}^M$ are moved and the $\tilde{\psi}_{k,c} = (\tilde{\rho}_{k,c}, \tilde{\rho} \tilde{E}_{k,c}, \tilde{\rho} \tilde{U}_{k,c})$ are remapped thanks to (39) and $\tilde{m}_{k,c}, \tilde{\alpha}_{k,c}, \tilde{\mathbf{X}}_{k,c}$ are computed for cells $c \in \mathcal{C}^M$.

Since $\mathcal{C}^M \cap \mathcal{C}^P \neq \{\emptyset\}$, one should note that cell included in this intersection are remapped at each step of the algorithm.

At the end of remapping, only the partial values of the physical variables per unit of volume are known. A this step, a first point is to compute the physical variables per unit of mass. The remapped partial total energy is given using $\tilde{E}_{k,c} = (\tilde{\rho} \tilde{E})_{k,c}/\tilde{\rho}_{k,c}$. However, this is different for the remapped partial velocity $\tilde{\mathbf{U}}_{k,c}$. Indeed, as explained in the second part of this paper, the Lagrangian computation of the velocity is done in Cartesian geometry. For this reason, the remapped velocity is deduced from the $(\tilde{\rho} \tilde{\mathbf{U}})_{k,c}$ through $\tilde{\mathbf{U}}_{k,c} = (\tilde{\rho} \tilde{\mathbf{U}})_{k,c}^{pl}/\tilde{\rho}_{k,c}^{pl}$ using the *planar* remapped density and momentum given through (39) and (42) with $\mathcal{R} = 1$. The second point is dedicated to the reconstruction of the global values required for the next Lagrangian step. To this end, a classical procedure is to use specific averages

$$ \tilde{\phi}_c = \frac{1}{\tilde{m}_c} \sum_k \tilde{m}_{k,c} \tilde{\phi}_{k,c}, \qquad (44) $$

with the global mass and density deduced from

$$ \tilde{m}_c = \sum_k \tilde{m}_{k,c} \tilde{\alpha}_{k,c} \quad \text{and} \quad \tilde{\rho}_c = \sum_k \tilde{\rho}_{k,c} \tilde{\alpha}_{k,c}. \qquad (45) $$
At last, thermodynamical variables as pressure $P$ and internal energy $\varepsilon$ are
obtained thanks to specific thermodynamical closures as done in [11].

**6. Numerical results**

We present in this section several numerical test cases performed using
the CCALE-MOF computing procedure detailed in [10, 11] including the
various development proposed in this paper. In the sequel, all the materials
are governed by perfect gas equation of state $p = \rho e(\gamma - 1)$, where $\gamma$ stands
for the polytropic index of gas.

6.1. Sedov problem

<figure><img src="image_16.png" /><figcaption>Figure 14: Initial grid and material positions for the Sedov problem.</figcaption></figure>

We present in this first section a Sedov problem for a point blast in a uni-
form medium with spherical symmetry. We use this test case to compare our
new formulation with the original EUCCLHYD scheme in pure Lagrangian
and coupled to the CCALE-MOF procedure. The initial conditions are given
At last, thermodynamical variables as pressure *P* and internal energy *ε* are
obtained thanks to specific thermodynamical closures as done in [11].

**6. Numerical results**

We present in this section several numerical test cases performed using
the CCALE-MOF computing procedure detailed in [10, 11] including the
various development proposed in this paper. In the sequel, all the materials
are governed by perfect gas equation of state p = ρe(γ − 1), where γ stands
for the polytropic index of gas.

6.1. Sedov problem

<figure><img src="image_15.png" /><figcaption>Figure 14: Initial grid and material positions for the Sedov problem.</figcaption></figure>

We present in this first section a Sedov problem for a point blast in a uni-
form medium with spherical symmetry. We use this test case to compare our
new formulation with the original EUCCLHYD scheme in pure Lagrangian
and coupled to the CCALE-MOF procedure. The initial conditions are given
by $(\rho^0, P^0, U^0) = (1, 10^{-6}, 0)$ in a spherical domain of radius 1.2 except in the cell at the origin $(0, 0)$ where an initial delta-function energy source is set through the pressure

$$P_{or} = (\gamma - 1)\rho_{or}\frac{\mathcal{E}_0}{V_{or}},$$

with $V_{or}$ the volume of the origin cell and $\mathcal{E}_0 = 0.851072$ is the total amount of released energy. The fluid has its polytropic index $\gamma$ equal to $\frac{7}{5}$. Contrary to the original single material test case, we add here three artificial interfaces, to test our multi-material CCALE-MOF algorithm. These interfaces are initially located for a radius equals to 0.1, 0.2 and 0.3 (see FIG.14).

Here we consider both Lagrangian and ALE computations for an initial unstructured mesh depicted on FIG.14. This grid is obtained after one rezoning step, with $\omega_p = 1$ of an unstructured mesh initially paved with 500 quadrangular cells. Numerical results are depicted on FIG.15 and FIG.16 for a final time of $t_{end} = 1$ and compared to the analytical solution computed using self-similar arguments as done in [11]. It consists of a diverging shock wave whose front is exactly localized at radius $R = 1$. As it is illustrated on FIG.15, the pure Lagrangian solutions are in good agreement with the analytical one for both approaches. We can notice that the new formulation is less dissipative as we reach a higher density level in the shock region. Indeed for the Lagrangian method as for the CCALE-MOF one the shock location is well resolved without any spurious oscillation (FIG.16). In addition, this simple problem underlines the robustness (better mesh quality near the origin) and accuracy (shock location) of the axisymmetric CCALE-MOF approach especially when considering multi-material flows whose interfaces are well captured thanks to the MOF reconstruction (see FIG.16).
by $(\rho^0, P^0, U^0) = (1, 10^{-6}, 0)$ in a spherical domain of radius 1.2 except in the cell at the origin $(0, 0)$ where an initial delta-function energy source is set through the pressure

$$P_{or} = (\gamma - 1)\rho_{or}\frac{\mathcal{E}_0}{V_{or}},$$

with $V_{or}$ the volume of the origin cell and $\mathcal{E}_0 = 0.851072$ is the total amount of released energy. The fluid has its polytropic index $\gamma$ equal to $\frac{7}{5}$. Contrary to the original single material test case, we add here three artificial interfaces, to test our multi-material CCALE-MOF algorithm. These interfaces are initially located for a radius equals to 0.1, 0.2 and 0.3 (see FIG.14).

Here we consider both Lagrangian and ALE computations for an initial unstructured mesh depicted on FIG.14. This grid is obtained after one rezoning step, with $\omega_p = 1$ of an unstructured mesh initially paved with 500 quadrangular cells. Numerical results are depicted on FIG.15 and FIG.16 for a final time of $t_{end} = 1$ and compared to the analytical solution computed using self-similar arguments as done in [11]. It consists of a diverging shock wave whose front is exactly localized at radius $R = 1$. As it is illustrated on FIG.15, the pure Lagrangian solutions are in good agreement with the analytical one for both approaches. We can notice that the new formulation is less dissipative as we reach a higher density level in the shock region. Indeed for the Lagrangian method as for the CCALE-MOF one the shock location is well resolved without any spurious oscillation (FIG.16). In addition, this simple problem underlines the robustness (better mesh quality near the origin) and accuracy (shock location) of the axisymmetric CCALE-MOF approach especially when considering multi-material flows whose interfaces are well captured thanks to the MOF reconstruction (see FIG.16).
<figure><img src="image_17.png" /><figcaption>Figure 15: Sedov problem. From the top to the bottom: Interface positions, density maps, density profiles defined as a function of the cell center radius compared to the analytical solution at final time step for the pure Lagrangian computation using new scheme (left) and the original EUCCLHYD scheme (right).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: Sedov problem. From the top to the bottom: Interface positions, density maps, density profiles defined as a function of the cell center radius compared to the analytical solution at final time step for the pure Lagrangian computation using new scheme (left) and the original EUCCLHYD scheme (right).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 16: Sedov problem. From the top to the bottom: Interface positions, density maps, density profiles defined as a function of the cell center radius compared to the analytical solution at final time step for the new CCALE-MOF procedure (left) and the EUCCLHYD CCALE-MOF procedure (right).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Sedov problem. From the top to the bottom: Interface positions, density maps, density profiles defined as a function of the cell center radius compared to the analytical solution at final time step for the new CCALE-MOF procedure (left) and the EUCCLHYD CCALE-MOF procedure (right).</figcaption></figure>
We point out that during the Lagrangian computation, non-convex cells appeared. This may lead to interface reconstruction failure when considering multi-material flows. As illustrated by the previous numerical results, the proposed CCALE-MOF algorithm remains adapted to treat such configuration without any difficulty demonstrating once again its robustness.

## 6.2. Axi-symmetric triple point problem

<figure><img src="image_19.png" /><figcaption>Figure 17: Axi-symmetric triple point problem: geometry and initial data.</figcaption></figure>

We consider in this part a three-material problem that corresponds to a three-state Riemann problem in an axisymmetric geometry. This problem has been wisely studied in Cartesian geometry and here we propose new results for cylindrical geometry. The computational domain is rectangular and composed of three regions (blue, green, red) whose dimensions are depicted on FIG.17. The top, left and right boundaries are closed thanks to walls. A symmetry condition is applied to the bottom boundary corresponding to the X-axis axi-symmetry. Initially, the blue region contains a fluid with high pressure and density taken equal to $(\rho_1, p_1) = (1, 1)$. The green region contains a low density and pressure fluid whose initial state is $(\rho_2, p_2) = (0.1, 0.125)$. The third fluid in the red region, initially has a low pressure and an high density equal to $(\rho_3, p_3) = (1, 0.1)$. At the beginning of the computation, all
We point out that during the Lagrangian computation, non-convex cells appeared. This may lead to interface reconstruction failure when considering multi-material flows. As illustrated by the previous numerical results, the proposed CCALE-MOF algorithm remains adapted to treat such configuration without any difficulty demonstrating once again its robustness.

## 6.2. Axi-symmetric triple point problem

<figure><img src="image_18.png" /><figcaption>Figure 17: Axi-symmetric triple point problem: geometry and initial data.</figcaption></figure>

We consider in this part a three-material problem that corresponds to a three-state Riemann problem in an axisymmetric geometry. This problem has been wisely studied in Cartesian geometry and here we propose new results for cylindrical geometry. The computational domain is rectangular and composed of three regions (blue, green, red) whose dimensions are depicted on FIG.17. The top, left and right boundaries are closed thanks to walls. A symmetry condition is applied to the bottom boundary corresponding to the X-axis axi-symmetry. Initially, the blue region contains a fluid with high pressure and density taken equal to $(\rho_1, p_1) = (1, 1)$. The green region contains a low density and pressure fluid whose initial state is $(\rho_2, p_2) = (0.1, 0.125)$. The third fluid in the red region, initially has a low pressure and an high density equal to $(\rho_3, p_3) = (1, 0.1)$. At the beginning of the computation, all
fluids are supposed to be at rest then $U_1 = U_2 = U_3 = 0$. The blue and green material have the same polytropic index $\gamma_1 = \gamma_2 = 1.5$, despite the red one has $\gamma_3 = 1.4$.

<figure><img src="image_20.png" /><figcaption>Figure 18: Axi-symmetric triple point problem. Mesh and material positions at $t = 5$ for Eulerian computation.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19: Axi-symmetric triple point problem. Mesh and material positions at $t = 5$ for ALE computation.</figcaption></figure>

The computation using the presented axisymmetric extension of the CCALE-MOF algorithm is made on a grid initially paved with 140 × 60 square cells until a final time $t_f = 5$. For this simulation, comparison with a full Lagrangian computation can not be performed since it suffers from important
fluids are supposed to be at rest then $U_1 = U_2 = U_3 = 0$. The blue and green material have the same polytropic index $\gamma_1 = \gamma_2 = 1.5$, despite the red one has $\gamma_3 = 1.4$.

<figure><img src="image_19.png" /><figcaption>Figure 18: Axi-symmetric triple point problem. Mesh and material positions at $t = 5$ for Eulerian computation.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Axi-symmetric triple point problem. Mesh and material positions at $t = 5$ for ALE computation.</figcaption></figure>

The computation using the presented axisymmetric extension of the CCALE-MOF algorithm is made on a grid initially paved with 140 × 60 square cells until a final time $t_f = 5$. For this simulation, comparison with a full Lagrangian computation can not be performed since it suffers from important
mesh tangling as shown in [18]. However comparison to full Eulerian simulations is done. In this case, nodes are moved to their initial positions during the rezoning step. Numerical results for both ALE and Eulerian methods representing interfaces and meshes are depicted on FIG.18-19. As expected, since there is a shock wave with high speed that propagates from the heavy material (blue) to the light one (red), the interface is sheared at the triple point producing a Kelvin-Helmholtz like instability. Here, comparison to planar 2D computations [10] demonstrates that axisymmetric geometry particularly affects the vortex shape that is 3D. Although the global behavior of the solutions is very similar comparing ALE approach to the Eulerian one.

### 6.3. Spherical Air-Helium shock/bubble interaction test

<figure><img src="image_22.png" /><figcaption>Figure 20: Air-Helium shock/bubble interaction. Initial geometry and data.</figcaption></figure>

We deal in this part with the numerical simulation of the experiment of [12] concerning the impact of a Mach 1.25 shock travelling through the air onto a spherical bubble of Helium. To this goal, let us consider a rectangular domain of dimensions [0, 0.65] x [0, 0.0445] initially full of Air of data ($ρ_1, P_1$) = (0.182, 10⁵) except in a half disc centered in (0, 0.32) of radius
mesh tangling as shown in [18]. However comparison to full Eulerian simulations is done. In this case, nodes are moved to their initial positions during the rezoning step. Numerical results for both ALE and Eulerian methods representing interfaces and meshes are depicted on FIG.18-19. As expected, since there is a shock wave with high speed that propagates from the heavy material (blue) to the light one (red), the interface is sheared at the triple point producing a Kelvin-Helmholtz like instability. Here, comparison to planar 2D computations [10] demonstrates that axisymmetric geometry particularly affects the vortex shape that is 3D. Although the global behavior of the solutions is very similar comparing ALE approach to the Eulerian one.

### 6.3. Spherical Air-Helium shock/bubble interaction test

<figure><img src="image_21.png" /><figcaption>Figure 20: Air-Helium shock/bubble interaction. Initial geometry and data.</figcaption></figure>

We deal in this part with the numerical simulation of the experiment of [12] concerning the impact of a Mach 1.25 shock travelling through the air onto a spherical bubble of Helium. To this goal, let us consider a rectangular domain of dimensions [0, 0.65] x [0, 0.0445] initially full of Air of data ($\rho_1, P_1$) = (0.182, 10⁵) except in a half disc centered in (0, 0.32) of radius
<table><caption>Table 1: Air-Helium shock/bubble interaction: EOS parameters.</caption><thead><tr><th>Fluid</th><th>Polytropic index γ</th><th>Molar mass <i>M</i></th></tr></thead><tbody><tr><td>Air</td><td>1.4</td><td>28.963</td></tr><tr><td>Helium</td><td>1.648</td><td>5.269 × 10<sup>-3</sup></td></tr></tbody></table>

0.0225 that contains Helium characterized by $(\rho_2, P_2) = (1, 10^5)$ as depicted on FIG.20. Here, spherical geometry is obtained thanks to a rotation around the X-axis. EOS parameters for each fluids are stated on TAB.1. Wall boundary and symmetry conditions are respectively chosen for the left, top boundaries. Despite, a piston-like condition is imposed to the right one for an incoming velocity equal to $U^* = (u^*, 0)$. Here, the horizontal velocity $u^*$ is computed thanks to Rankine-Hugoniot conditions and is given by $u^* = -140.312$ corresponding to an incident shock moving at the velocity $D_c = -467.707$.

The domain is initially paved with a structured cartesian grid composed of 520 × 72 cells. Here, the bubble is directly initialized through the volume fraction on this mesh. Computations are done for the multi-material axisymmetric CCALE-MOF for a final time chosen equal to $t_f = t_i + 600 \times 10^{-6}$ where $t_i = 657.463 \times 10^{-6}$ corresponds to the time of the shock/bubble interaction. Here once again, simulations can not be achieved using pure Lagrangian framework due to the apparition of important mesh distortion. Numerical results associated to the Schlieren density profiles [13] and interface positions deduced from the MOF reconstruction are respectively depicted on FIG.22 and FIG.23. Let us note that each pictures are obtained thanks to an axial symmetry with respect to the X-axis. Comparisons between the Schlieren density profiles and the shadow-graphs of the experiment show a good agree-
<table><caption>Table 1: Air-Helium shock/bubble interaction: EOS parameters.</caption><thead><tr><th>Fluid</th><th>Polytropic index γ</th><th>Molar mass <i>M</i></th></tr></thead><tbody><tr><td>Air</td><td>1.4</td><td>28.963</td></tr><tr><td>Helium</td><td>1.648</td><td>5.269 × 10<sup>-3</sup></td></tr></tbody></table>

0.0225 that contains Helium characterized by $(\rho_2, P_2) = (1, 10^5)$ as depicted on Fig. 20. Here, spherical geometry is obtained thanks to a rotation around the X-axis. EOS parameters for each fluids are stated on Tab. 1. Wall boundary and symmetry conditions are respectively chosen for the left, top boundaries. Despite, a piston-like condition is imposed to the right one for an incoming velocity equal to $\mathbf{U}^* = (u^*, 0)$. Here, the horizontal velocity $u^*$ is computed thanks to Rankine-Hugoniot conditions and is given by $u^* = -140.312$ corresponding to an incident shock moving at the velocity $D_c = -467.707$.

The domain is initially paved with a structured cartesian grid composed of 520 × 72 cells. Here, the bubble is directly initialized through the volume fraction on this mesh. Computations are done for the multi-material axisymmetric CCALE-MOF for a final time chosen equal to $t_f = t_i + 600 \times 10^{-6}$ where $t_i = 657.463 \times 10^{-6}$ corresponds to the time of the shock/bubble interaction. Here once again, simulations can not be achieved using pure Lagrangian framework due to the apparition of important mesh distortion. Numerical results associated to the Schlieren density profiles [13] and interface positions deduced from the MOF reconstruction are respectively depicted on Fig. 22 and Fig. 23. Let us note that each pictures are obtained thanks to an axial symmetry with respect to the X-axis. Comparisons between the Schlieren density profiles and the shadow-graphs of the experiment show a good agree-
ment, especially when observing the bubble shape deformations. Moreover,
waves generated by the initial shock are well localized and illustrate multiple
reflections and refractions especially on the bubble and the domain bound-
aries. These main points clearly demonstrate the accuracy and the robustness
of the method and validate the axisymmetric CCALE-MOF approach when
computing spherical test-cases coming from experiment.

6.4. Spherical implosion

<figure><img src="image_23.png" /><figcaption>Figure 21: Multi-mode implosion in spherical geometry. Initial geometry and data.</figcaption></figure>

The last test-case of this paper deals with the numerical computation of a
spherical implosion as initially treated in [27]. The interest of this simulation
is twofold. First, this is a realistic problem quite close to those encountered
in Ignition Confinement Fusion (ICF) simulation. Then, it allows to test
the capability of the multi-material CCALE-MOF algorithm with hybrid
rezoning.
ment, especially when observing the bubble shape deformations. Moreover,
waves generated by the initial shock are well localized and illustrate multiple
reflections and refractions especially on the bubble and the domain bound-
aries. These main points clearly demonstrate the accuracy and the robustness
of the method and validate the axisymmetric CCALE-MOF approach when
computing spherical test-cases coming from experiment.

6.4. Spherical implosion

<figure><img src="image_22.png" /><figcaption>Figure 21: Multi-mode implosion in spherical geometry. Initial geometry and data.</figcaption></figure>

The last test-case of this paper deals with the numerical computation of a
spherical implosion as initially treated in [27]. The interest of this simulation
is twofold. First, this is a realistic problem quite close to those encountered
in Ignition Confinement Fusion (ICF) simulation. Then, it allows to test
the capability of the multi-material CCALE-MOF algorithm with hybrid
rezoning.
<figure><img src="image_24.png" /><figcaption>Figure 22: Spherical Air-Helium shock/bubble interaction. Schlieren diagram of density. Axi-symmetric CCALE-MOF results (on the left) compared to experimental results (on the right) [12].</figcaption></figure>
<figure><img src="image_23.png" /></figure>

$$t = t_i + 20 \times 10^{-6}$$

<figure><img src="image_24.png" /></figure>

$$t = t_i + 145 \times 10^{-6}$$

<figure><img src="image_25.png" /></figure>

$$t = t_i + 223 \times 10^{-6}$$

<figure><img src="image_26.png" /></figure>

$$t = t_i + 350 \times 10^{-6}$$

<figure><img src="image_27.png" /></figure>

$$t = t_i + 600 \times 10^{-6}$$

<caption>Figure 22: Spherical Air-Helium shock/bubble interaction. Schlieren diagram of density. Axi-symmetric CCALE-MOF results (on the left) compared to experimental results (on the right) [12].</caption>
<figure><img src="image_25.png" /></figure>

$$t = t_i + 20 \times 10^{-6}$$

<figure><img src="image_26.png" /></figure>

$$t = t_i + 145 \times 10^{-6}$$

<figure><img src="image_27.png" /></figure>

$$t = t_i + 223 \times 10^{-6}$$

<figure><img src="image_28.png" /></figure>

$$t = t_i + 350 \times 10^{-6}$$

<figure><img src="image_29.png" /></figure>

$$t = t_i + 600 \times 10^{-6}$$

<caption>Figure 23: Spherical Air-Helium shock/bubble interaction. Mesh and material interface evolution after the shock hits the bubble at time $t_i = 657.463 \times 10^{-6}$.</caption>
<figure><img src="image_28.png" /><figcaption>Figure 23: Spherical Air-Helium shock/bubble interaction. Mesh and material interface evolution after the shock hits the bubble at time $t_i = 657.463 \times 10^{-6}$.</figcaption></figure>
Here we focus on the treatment of perturbed interfaces where compressible Rayleigh-Taylor instabilities occur.

Let us consider a spherical ball of light fluid ($r \in [0, 10]$) initially surrounded by a shell of heavy fluid ($R \in [10, 12]$) as depicted on Fig.21. For both fluid the polytropic index is the same $\gamma_l = \gamma_h = \frac{5}{3}$. The initial pressures and densities are $(\rho_l, p_l) = (0.05, 0.1)$ and $(\rho_h, p_h) = (1, 0.1)$. The implosion is driven by imposing the following pressure law on the dense shell boundary

$$p^*(t) = \begin{cases} 10 & \text{if } t \in [0, 0.5], \\ 12 - 4t & \text{if } t \in [0.5, 3]. \end{cases}$$

Finally, the interface between the light and the heavy fluids is initially perturbed according to the law

$$r_p^{per} = r_p(1 + a_0 \mathcal{D}(r_p) P_l(\cos(\theta_p)))$$

with the damping factor

$$\mathcal{D}(r_p) = \begin{cases} 1 - \frac{r_p - r_i}{r_e - r_i} & \text{if } r_p \in [r_i, r_e], \\ 1 - \frac{r_i - r_p}{r_i} & \text{if } r_p \in [0, r_i]. \end{cases}$$

where $r_i^{per}$ denotes the perturbed radius and $a_0$ is the amplitude of the perturbation. Finally, $P_l$ is the $l^{th}$ Legendre polynomial. In the sequel $l = 10$ and several values of $a_0$ are considered from the non-perturbed case $a_0 = 0$, to weakly and strongly perturbed one with respectively $a_0 = 2 \times 10^{-4}$ and $a_0 = 1 \times 10^{-3}$.

Computations are made for two different meshes until the final time $t_f = 3$. The first one is a polar grid displayed on Fig.24-(left) composed of 90 × 40 cells. Size of cells in the radial direction have been chosen respecting a mass
Here we focus on the treatment of perturbed interfaces where compressible Rayleigh-Taylor instabilities occur.

Let us consider a spherical ball of light fluid ($r \in [0, 10]$) initially surrounded by a shell of heavy fluid ($R \in [10, 12]$) as depicted on Fig.21. For both fluid the polytropic index is the same $\gamma_l = \gamma_h = \frac{5}{3}$. The initial pressures and densities are $(\rho_l, p_l) = (0.05, 0.1)$ and $(\rho_h, p_h) = (1, 0.1)$. The implosion is driven by imposing the following pressure law on the dense shell boundary

$$p^*(t) = \begin{cases} 10 & \text{if } t \in [0, 0.5], \\ 12 - 4t & \text{if } t \in [0.5, 3]. \end{cases}$$

Finally, the interface between the light and the heavy fluids is initially perturbed according to the law

$$r_p^{per} = r_p(1 + a_0 \mathcal{D}(r_p) P_l(\cos(\theta_p)))$$

with the damping factor

$$\mathcal{D}(r_p) = \begin{cases} 1 - \frac{r_p - r_i}{r_e - r_i} & \text{if } r_p \in [r_i, r_e], \\ 1 - \frac{r_i - r_p}{r_i} & \text{if } r_p \in [0, r_i]. \end{cases}$$

where $r_i^{per}$ denotes the perturbed radius and $a_0$ is the amplitude of the perturbation. Finally, $P_l$ is the $l^{th}$ Legendre polynomial. In the sequel $l = 10$ and several values of $a_0$ are considered from the non-perturbed case $a_0 = 0$, to weakly and strongly perturbed one with respectively $a_0 = 2 \times 10^{-4}$ and $a_0 = 1 \times 10^{-3}$.

Computations are made for two different meshes until the final time $t_f = 3$. The first one is a polar grid displayed on Fig.24-(left) composed of 90 × 40 cells. Size of cells in the radial direction have been chosen respecting a mass
radial spacing deduced from the equivalent one-dimensional test case. The other grid, is obtained after an hybrid regularization for $\omega_p = 1$ of an un-structured mesh initially paved with 3200 quadrangular cells respecting the mass radial spacing (see FIG.24-(right)).

*Non-perturbed case with $a_0 = 0$. As a first study, we test the behavior of our algorithm in axisymmetric geometries in pure Lagrange computation for both meshes. As shown on FIG.25, numerical results for both meshes are similar. Nevertheless, the method remains faster on the unstructured mesh. Indeed, it has the advantage to not impose a drastic time step for computation due to triangular cells with high aspect ratio in the polar mesh as shown in [11].*

<figure><img src="image_30.png" /><figcaption>Figure 24: Spherical implosion. Initial polar (left) and unstructured (right) grids.</figcaption></figure>

*Weakly perturbed case with $a_0 = 2 \times 10^{-4}$. Now, we investigate the capability of our CCALE-MOF algorithm to treat perturbed interfaces on both non-structured and polar meshes. To this end, comparisons with pure Lagrangian results are first achieved for weakly perturbed interfaces imposing*
radial spacing deduced from the equivalent one-dimensional test case. The other grid, is obtained after an hybrid regularization for $\omega_p = 1$ of an unstructured mesh initially paved with 3200 quadrangular cells respecting the mass radial spacing (see FIG.24-(right)).

*Non-perturbed case with $a_0 = 0$. As a first study, we test the behavior of our algorithm in axisymmetric geometries in pure Lagrange computation for both meshes. As shown on FIG.25, numerical results for both meshes are similar. Nevertheless, the method remains faster on the unstructured mesh. Indeed, it has the advantage to not impose a drastic time step for computation due to triangular cells with high aspect ratio in the polar mesh as shown in [11].*

<figure><img src="image_29.png" /><figcaption>Figure 24: Spherical implosion. Initial polar (left) and unstructured (right) grids.</figcaption></figure>

*Weakly perturbed case with $a_0 = 2 \times 10^{-4}$. Now, we investigate the capability of our CCALE-MOF algorithm to treat perturbed interfaces on both non-structured and polar meshes. To this end, comparisons with pure Lagrangian results are first achieved for weakly perturbed interfaces imposing*
<figure><img src="image_31.png" /><figcaption>Figure 25: Spherical implosion without deformation. Mesh and density for polar (left) and unstructured (right) grids at final time $t_f = 3$.</figcaption></figure>

$a_0 = 2 \times 10^{-4}$. Here for both polar and hybrid meshes, the GCNS is used. As demonstrated on FIG. 26, for the polar mesh as well as for the non-structured mesh, ALE results, especially concerning the interface deformation, are in very good agreement to those obtained thanks to pure Lagrangian computations. Furthermore, one should note that for the ALE computation on polar grid the quality of the mesh is improved near the origin. Indeed, the central cells are not systematically shifted to the origin contrary to computations achieved using CNS rezoning.

Strongly perturbed case with $a_0 = 1 \times 10^{-3}$. Finally, we perform a computation of this implosion for a more perturbed interface choosing $a_0$ five times greater than previously with $a_0 = 1 \times 10^{-3}$. Due to mesh tangling, this is not possible to purchase such a test case using only Lagrangian method whose computation fails for $t > t_{fail} = 2.6$. Here, only results obtained thanks to our axisymmetric multi-material CCALE-MOF are presented. Contrary, to
<figure><img src="image_30.png" /><figcaption>Figure 25: Spherical implosion without deformation. Mesh and density for polar (left) and unstructured (right) grids at final time $t_f = 3$.</figcaption></figure>

$a_0 = 2 \times 10^{-4}$. Here for both polar and hybrid meshes, the GCNS is used. As demonstrated on FIG. 26, for the polar mesh as well as for the non-structured mesh, ALE results, especially concerning the interface deformation, are in very good agreement to those obtained thanks to pure Lagrangian computations. Furthermore, one should note that for the ALE computation on polar grid the quality of the mesh is improved near the origin. Indeed, the central cells are not systematically shifted to the origin contrary to computations achieved using CNS rezoning.

Strongly perturbed case with $a_0 = 1 \times 10^{-3}$. Finally, we perform a computation of this implosion for a more perturbed interface choosing $a_0$ five times greater than previously with $a_0 = 1 \times 10^{-3}$. Due to mesh tangling, this is not possible to purchase such a test case using only Lagrangian method whose computation fails for $t > t_{fail} = 2.6$. Here, only results obtained thanks to our axisymmetric multi-material CCALE-MOF are presented. Contrary, to
<figure><img src="image_32.png" /><figcaption>Figure 26: Spherical implosion with small deformation. Mesh and density for Lagrangian (top) and ALE (bottom) computations at final time $t_f = 3$ for both polar (left) and unstructured (right) grids.</figcaption></figure>

Lagrangian computations, the multi-material ALE simulations run without any difficulties thanks to specific rezoning. For both grids, final results (see Fig. 27) are very close. In particular we note the Rayleigh-Taylor instability has grown in a same way leading to similar interface shape deformation at final time.
<figure><img src="image_31.png" /><figcaption>Figure 26: Spherical implosion with small deformation. Mesh and density for Lagrangian (top) and ALE (bottom) computations at final time $t_f = 3$ for both polar (left) and unstructured (right) grids.</figcaption></figure>

Lagrangian computations, the multi-material ALE simulations run without any difficulties thanks to specific rezoning. For both grids, final results (see FIG.27) are very close. In particular we note the Rayleigh-Taylor instability has grown in a same way leading to similar interface shape deformation at final time.
<figure><img src="image_33.png" /><figcaption>Figure 27: Spherical implosion with important deformation. Mesh and density for ALE computation for both polar (left) and unstructured (right) grids at final time $t_f = 3$.</figcaption></figure>

## 7. Conclusion and future work

In this paper, we have presented several extensions concerning a Cell-Centered Arbitrary Lagrangian-Eulerian (CCALE) strategy using the Moment of Fluid (MOF) interface reconstruction devoted to the numerical simulation of multi-material compressible flows especially in axisymmetric geometry on both polar and Cartesian unstructured meshes. To this end, we have introduced a simple and unified formulation of the Lagrangian scheme relying on an area-weighted formulation, a multi-material MOF interface reconstruction, a new formulation of rezoning for both polar and Cartesian grids and finally a general hybrid remap procedure for both axisymmetric and Cartesian geometry. As demonstrated on several academical as well as ICF-like test cases, the proposed method remains accurate and robust.

As future work, we plan to incorporate the proposed method in the multi-physic code CHIC dedicated to the simulation of ICF experiment. The main
<figure><img src="image_32.png" /><figcaption>Figure 27: Spherical implosion with important deformation. Mesh and density for ALE computation for both polar (left) and unstructured (right) grids at final time $t_f = 3$.</figcaption></figure>

## 7. Conclusion and future work

In this paper, we have presented several extensions concerning a Cell-Centered Arbitrary Lagrangian-Eulerian (CCALE) strategy using the Moment of Fluid (MOF) interface reconstruction devoted to the numerical simulation of multi-material compressible flows especially in axisymmetric geometry on both polar and Cartesian unstructured meshes. To this end, we have introduced a simple and unified formulation of the Lagrangian scheme relying on an area-weighted formulation, a multi-material MOF interface reconstruction, a new formulation of rezoning for both polar and Cartesian grids and finally a general hybrid remap procedure for both axisymmetric and Cartesian geometry. As demonstrated on several academical as well as ICF-like test cases, the proposed method remains accurate and robust.

As future work, we plan to incorporate the proposed method in the multi-physic code CHIC dedicated to the simulation of ICF experiment. The main
goal is to treat eventually more general configurations notably coupling re-
alistic EOS, laser energy deposition, with multi-material hydrodynamics in
the lines of [5].
goal is to treat eventually more general configurations notably coupling re-
alistic EOS, laser energy deposition, with multi-material hydrodynamics in
the lines of [5].
References

[1] H.T. Ahn, M.J. Shashkov: *Multi-material interface reconstruction on generalized polyhedral meshes*, J. Comput. Phys., 226(2):2096-2132, 2007.

[2] H.R. Anbarlooei, K. Mazaheri: *'Moment of fluid' interface reconstruction method in axisymmetric coordinates*, Int. J. Numer. Meth. Biomed. Engng., 27(10):1640-1651, 2011.

[3] A.J. Barlow, P.L. Roe: *A cell centred Lagrangian Godunov scheme for shock hydrodynamics*, Comput. Fluids, 46(1):133-136, 2011.

[4] M. Berndt, J. Breil, S. Galera, M. Kucharik, P.-H. Maire, M. Shashkov: *Two-step hybrid conservative remapping for multi-material arbitrary Lagrangian-Eulerian methods*, J. Comput. Phys., 230(17):6664-6687, 2011.

[5] J. Breil, S. Galera, P.-H. Maire: *Multi-material ALE computation in Inertial Confinement Fusion CHIC*, Comput. Fluids, 46(1):161-167, 2011.

[6] J. Breil, S. Galera, P.-H. Maire: *A two-dimensional VOF interface reconstruction in a multi-material cell-centered ALE scheme*, Int. J. Numer. Meth. Fluids, 65(11-12):1351-1364, 2011.

[7] G. Carr, S. Del Pino, B. Desprs: *A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension*, J. Comput. Phys., 228(14):5160-5183, 2009.
References

[1] H.T. Ahn, M.J. Shashkov: *Multi-material interface reconstruction on generalized polyhedral meshes*, J. Comput. Phys., 226(2):2096-2132, 2007.

[2] H.R. Anbarlooei, K. Mazaheri: *'Moment of fluid' interface reconstruction method in axisymmetric coordinates*, Int. J. Numer. Meth. Biomed. Engng., 27(10):1640-1651, 2011.

[3] A.J. Barlow, P.L. Roe: *A cell centred Lagrangian Godunov scheme for shock hydrodynamics*, Comput. Fluids, 46(1):133-136, 2011.

[4] M. Berndt, J. Breil, S. Galera, M. Kucharik, P.-H. Maire, M. Shashkov: *Two-step hybrid conservative remapping for multi-material arbitrary Lagrangian-Eulerian methods*, J. Comput. Phys., 230(17):6664-6687, 2011.

[5] J. Breil, S. Galera, P.-H. Maire: *Multi-material ALE computation in Inertial Confinement Fusion CHIC*, Comput. Fluids, 46(1):161-167, 2011.

[6] J. Breil, S. Galera, P.-H. Maire: *A two-dimensional VOF interface reconstruction in a multi-material cell-centered ALE scheme*, Int. J. Numer. Meth. Fluids, 65(11-12):1351-1364, 2011.

[7] G. Carr, S. Del Pino, B. Desprs: *A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension*, J. Comput. Phys., 228(14):5160-5183, 2009.
[8] V. Dyadechko, M. Shashkov: *Reconstruction of Multi-material Interfaces from Moment Data*, J. Comput. Phys., 227(11):5361-5384, 2008.

[9] J. K. Dukowicz: *A general, non-iterative Riemann solver for Godunov's method*, J. Comput. Phys., 61(1):119-137, 1985.

[10] S. Galera, J. Breil and P.-H. Maire: *A 2D unstructured multi-material Cell-Centered Arbitrary Lagrangian-Eulerian (CCALE) scheme using MOF interface reconstruction*, Comput. Fluids, 46(1):237-244, 2011.

[11] S. Galera, P.-H. Maire and J. Breil: *A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction*, J. Comput. Phys., 229(16):5755-5787, 2010.

[12] J.-F. Haas, B. Sturtevant: *Interaction of weak shock wave with cylindrical and spherical gas inhomogeneities*, J. Fluid. Mech., 181:41-76, 1987.

[13] A. Hadjadj, A. Kudryavtsev *Computation and flow visualization in high-speed aerodynamics*, Journal of Turbulence, 6(16), 2005.

[14] C.W. Hirt, A. Amsden, and J.L. Cook: *An arbitrary Lagrangian-Eulerian computing method for all flow speeds*, J. Comput. Phys., 14:227-253, 1974.

[15] P. Knupp: *Achieving finite element mesh quality via optimization of the Jacobian matrix norm and associated quantities. Part I– a framework for surface mesh optimization*, Int. J. Numer. Meth. Engng, 48:401-420, 2000.
[8] V. Dyadechko, M. Shashkov: *Reconstruction of Multi-material Interfaces from Moment Data*, J. Comput. Phys., 227(11):5361-5384, 2008.

[9] J. K. Dukowicz: *A general, non-iterative Riemann solver for Godunov's method*, J. Comput. Phys., 61(1):119-137, 1985.

[10] S. Galera, J. Breil and P.-H. Maire: *A 2D unstructured multi-material Cell-Centered Arbitrary Lagrangian-Eulerian (CCALE) scheme using MOF interface reconstruction*, Comput. Fluids, 46(1):237-244, 2011.

[11] S. Galera, P.-H. Maire and J. Breil: *A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction*, J. Comput. Phys., 229(16):5755-5787, 2010.

[12] J.-F. Haas, B. Sturtevant: *Interaction of weak shock wave with cylindrical and spherical gas inhomogeneities*, J. Fluid. Mech., 181:41-76, 1987.

[13] A. Hadjadj, A. Kudryavtsev *Computation and flow visualization in high-speed aerodynamics*, Journal of Turbulence, 6(16), 2005.

[14] C.W. Hirt, A. Amsden, and J.L. Cook: *An arbitrary Lagrangian-Eulerian computing method for all flow speeds*, J. Comput. Phys., 14:227-253, 1974.

[15] P. Knupp: *Achieving finite element mesh quality via optimization of the Jacobian matrix norm and associated quantities. Part I– a framework for surface mesh optimization*, Int. J. Numer. Meth. Engng, 48:401-420, 2000.
[16] M. Kucharik, J. Breil, S. Galera, P.-H. Maire, M. Berndt, M. Shashkov: *Hybrid remap for multi-material ALE*, Comput. Fluids, 46(1):293-297, 2011.

[17] M. Kucharik, R.V. Garimella, S.P. Schofield, M.J. Shashkov: *A comparative study of interface reconstruction methods for multi-material ALE simulations*, J. Comput. Phys., 229(7):2432-2452, 2010.

[18] R. Loubère, P.-H. Maire, M. Shashkov, J. Breil, S. Galera: *ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method*, J. Comput. Phys., 229(12):4724-4761, 2010.

[19] P.-H. Maire: *A high-order cell-centered Lagrangian scheme for compressible fluid flows in two-dimensional cylindrical geometry*, J. Comput. Phys., 228(18):6882-6915, 2009.

[20] P.-H. Maire: *A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes*, J. Comput. Phys., 228(7):2391-2425, 2009.

[21] P.-H. Maire, R. Abgrall, J. Breil, J. Ovadia: *A cell-centered Lagrangian scheme for two-dimensional compressible flow problems*, SIAM Journal of Scientific Computing, 29(4):1781-1824, 2007.

[22] P.-H. Maire: *Contribution to the numerical modeling of Inertial Confinement Fusion*, Habilitation à Diriger des Recherches, Bordeaux University, 2011; Available at: http://tel.archives-ouvertes.fr/docs/00/58/97/58/PDF/hdr_main.pdf.
[16] M. Kucharik, J. Breil, S. Galera, P.-H. Maire, M. Berndt, M. Shashkov: *Hybrid remap for multi-material ALE*, Comput. Fluids, 46(1):293-297, 2011.

[17] M. Kucharik, R.V. Garimella, S.P. Schofield, M.J. Shashkov: *A comparative study of interface reconstruction methods for multi-material ALE simulations*, J. Comput. Phys., 229(7):2432-2452, 2010.

[18] R. Loubère, P.-H. Maire, M. Shashkov, J. Breil, S. Galera: *ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method*, J. Comput. Phys., 229(12):4724-4761, 2010.

[19] P.-H. Maire: *A high-order cell-centered Lagrangian scheme for compressible fluid flows in two-dimensional cylindrical geometry*, J. Comput. Phys., 228(18):6882-6915, 2009.

[20] P.-H. Maire: *A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes*, J. Comput. Phys., 228(7):2391-2425, 2009.

[21] P.-H. Maire, R. Abgrall, J. Breil, J. Ovadia: *A cell-centered Lagrangian scheme for two-dimensional compressible flow problems*, SIAM Journal of Scientific Computing, 29(4):1781-1824, 2007.

[22] P.-H. Maire: *Contribution to the numerical modeling of Inertial Confinement Fusion*, Habilitation à Diriger des Recherches, Bordeaux University, 2011; Available at: http://tel.archives-ouvertes.fr/docs/00/58/97/58/PDF/hdr_main.pdf.
[23] L. G. Margolin, M. Shashkov: *Second-order sign-preserving conservative interpolation (remapping) on general grids*, J. Comput. Phys., 184(1):266-298, 2003.

[24] M.Shashkov: *Closure models for multidimensional cells in arbitrary Lagrangian-Eulerian hydrocodes*, Int. J. Numer. Meth. Fluids 56:1497-1504, 2008.

[25] P. Vachal, P.-H. Maire: *Discretizations for weighted condition number smoothing on general unstructured meshes*, Comput. Fluids, 46(1):479-485, 2011.

[26] P. Vachal, R.V. Garimella, M.J. Shashkov: *Untangling of 2D meshes in ALE simulations*, J. Comput. Phys., 196:627-644, 2004.

[27] D. L. Youngs: *Multi-mode implosion in cylindrical 3D geometry*, 11th International Workshop on the Physics of Compressible Turbulent Mixing (IWPCTM11), Santa Fe, 2008.
[23] L. G. Margolin, M. Shashkov: *Second-order sign-preserving conservative interpolation (remapping) on general grids*, J. Comput. Phys., 184(1):266-298, 2003.

[24] M.Shashkov: *Closure models for multidimensional cells in arbitrary Lagrangian-Eulerian hydrocodes*, Int. J. Numer. Meth. Fluids 56:1497-1504, 2008.

[25] P. Vachal, P.-H. Maire: *Discretizations for weighted condition number smoothing on general unstructured meshes*, Comput. Fluids, 46(1):479-485, 2011.

[26] P. Vachal, R.V. Garimella, M.J. Shashkov: *Untangling of 2D meshes in ALE simulations*, J. Comput. Phys., 196:627-644, 2004.

[27] D. L. Youngs: *Multi-mode implosion in cylindrical 3D geometry*, 11th International Workshop on the Physics of Compressible Turbulent Mixing (IWPCTM11), Santa Fe, 2008.