# An Immersed Boundary Method for pressure-based compressible solvers with applications to free-convection flows, acoustic wave propagation and thermal plasma

Sergiu Coseru, Sébastien Tanguy, Pierre Freton, Jean Jacques Gonzalez, Annafederica Urbano, Marie Bibal, Gauthier Bourdon

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

## To cite this version:

Sergiu Coseru, Sébastien Tanguy, Pierre Freton, Jean Jacques Gonzalez, Annafederica Urbano, et al. An Immersed Boundary Method for pressure-based compressible solvers with applications to free-convection flows, acoustic wave propagation and thermal plasma. Journal of Computational Physics, 2025, 524, pp.113714. 10.1016/j.jcp.2024.113714. hal-04890185

HAL Id: hal-04890185

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

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# An Immersed Boundary Method for pressure-based compressible solvers with applications to free-convection flows, acoustic wave propagation and thermal plasma

Sergiu Coseru<sup>a,b</sup>, Sébastien Tanguy<sup>a,*</sup>, Pierre Freton<sup>b</sup>, Jean-Jacques Gonzalez<sup>b</sup>, Annafederica Urbano<sup>c</sup>, Marie Bibal<sup>a</sup>, Gauthier Bourdon<sup>a</sup>

## ABSTRACT

Immersed Boundary Methods (IBM) are a practical class of methods that enable fluid computations in complex geometry while keeping a structured mesh. Most of the existing IBM have been developed in the framework of incompressible solvers, despite their significant interest to perform simulations in more complex configurations requiring a compressible solver. In the last years, pressure-based solvers met a growing interest to perform numerical simulations of compressible flows, due to their attractive features, as removing the stability condition on the acoustic time step, and being asymptotically preserving of the incompressible regime when the Mach number tends to zero. As this class of compressible solvers share many common features with classical projection methods for incompressible flows, our objective in this paper is to present an adaptation of an efficient and accurate IBM developed for an incompressible solver by Ng et al. in [1] to a pressure-based compressible solver recently published by Urbano et al. in [2]. The proposed algorithm benefits of the attractive properties of the original IBM proposed in [1] while being able to undertake simulations in much more complex configurations. In particular, we will present validations and illustrations of the proposed solver in various configurations as free-convection flows, acoustic waves propagating in a variable section pipe or interacting with a solid obstacle, as well as the description of a thermal plasma during an electric arc discharge in a gas.

## 1. Introduction

Since the pioneering works of [3], Immersed Boundary Methods became a popular class of methods for fluid computations in various configurations implying a complex geometry, as shown in [4–10], whether an external flow or an internal flow is considered. The assets of these methods are twofold, since they allow computations in complex geometries while keeping structured grids for which highly accurate numerical schemes are available, and they drastically alleviate the preliminary and difficult task of mesh designing. On the other hand, compared to computations on unstructured grids, IBM suffer from a lack of flexibility with regards to the grid adaptation if they are employed with a classical structured mesh. Whereas the latter point can be critical to undertake numerical simulations with high Reynolds number and thin boundary layers around solid obstacles, significant advances in grid adaptability with quad-tree and octree structures coupled to IBM, as presented in [11–14] open new opportunities in regard to both accuracy and convenience and have demonstrated the strong interest of such a coupling. In anterior works, other Adaptive Mesh Refinement strategies, based on block-structured approaches, have also been developed and coupled to embedded solid frontiers in [15] for a fully compressible solver. Next, improvements of this powerful methodology have been presented in [16]. Other couplings between

<footnote><sup>a</sup> IMFT, UMR 5502 CNRS-INPT-UPS, Université de Toulouse 2 All. du Professeur Camille Soula, 31400 Toulouse, France</footnote>

<footnote><sup>b</sup> Laplace, UMR CNRS 5213, 118 Rte de Narbonne, 31077 Toulouse, France</footnote>

<footnote><sup>c</sup> Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, 10 avenue Edouard Belin, Toulouse, BP 54032 - 31055, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: tanguy@imft.fr (S. Tanguy).</footnote>

<footnote>https://doi.org/10.1016/j.jcp.2024.113714</footnote>

<footnote>Received 26 April 2024; Received in revised form 16 November 2024; Accepted 29 December 2024</footnote>
IBM and a compressible solver have been proposed in [17,18], and further developed more recently in the following references [19–26]. Applications of such a numerical strategy to aerodynamic applications with average turbulent models are presented in [27]. However, most of the previously cited works are dedicated to fully explicit compressible solvers which are efficient for supersonic flows but not optimal to study compressible flows in the low Mach number regime, as shown in [28]. Compressible solvers can be of interest in many applications involving a low Mach number, as free-convection flows, compression or expansion of a gas in a closed tank, interaction of acoustic waves with solid frontiers or obstacles and also more specifically the description of thermal plasma during an electric arc discharge. This motivates the development of IBM in the framework of compressible solvers, which are well suited to the low Mach number regime. Such an attempt has been recently published in [29], however the proposed solver is restricted to configurations without acoustic waves. On the other hand, pressure-based compressible solvers [30,31] benefit from a growing interest since the last decade due to their attractive numerical properties. Indeed, such solvers allow stable computations without imposing the time step constraint related to acoustic wave propagation while being able to describe acoustical phenomena. Moreover, they are asymptotically preserving of the incompressible regime when the Mach number tends to zero. Indeed, for pressure-based compressible solvers, the pressure equation which results from the velocity splitting tends strictly to the incompressible Poisson equation for the pressure when the Mach number tends to zero. This property is of utmost importance when considering compressible two-phase flows, as in the following references [32, 33, 28, 34–37, 2, 38], for which low Mach number values are usual in the liquid phase due to the high value of the sound speed. Moreover, as pressure-based compressible solvers present many common features with projection methods for incompressible flows, they can widely benefit from IBM or two-phase flows solvers previously developed in the framework of incompressible flows. One of our purpose in this paper is to adapt a well-known IBM proposed in [1] for incompressible flows to a pressure-based solver for compressible flows proposed in [2]. The proposed coupling maintains the attractive features of the IBM proposed in [1], i.e. a sharp description of the solid-fluid frontier, second order accurate velocity and pressure boundary conditions on the solid-fluid frontier and a resulting scheme which preserves the symmetric definite positive matrix stemming from the spatial discretization of the Helmholtz equation for the pressure. Our long-term aim is to develop an efficient thermal plasma solver able to describe the growth of a plasma bubble during an electric arc discharge in a liquid bath. Indeed, in several industrial applications, such as pulsed metal forming by Electrohydraulic Discharges [39,40] or fracturing rock [41] an electric arc discharge in a liquid pool is carried out to generate intense pressure waves. As such applications involve a very complex coupling between several physical phenomena, efficient computational methods are required in order to improve our knowledge of the overall phenomena and to optimize industrial processes. Therefore, we will present also in this paper additional developments to couple our pressure-based compressible solver with computational modules for the description of a thermal plasma formation during an electric arc discharge, as in [42–44]. In this framework, the proposed IBM for compressible flows will be extended to a plasma model in order to describe the electrodes geometry on which specific boundary conditions have to be imposed. Finally, several original benchmarks are proposed to demonstrate the interest and the relevance of the proposed method to undertake numerical simulations involving physical phenomena typical of compressible or variable density flows, as free convection and acoustic waves propagation. Preliminary simulations involving thermal plasma induced by electric arc discharges between electrodes will also be presented. In particular, these simulations will highlight the ability of the proposed solver to capture pressure waves generation during a transient electric arc discharge while accounting for electrodes geometries by means of the proposed IBM. It is worth mentioning that the present solver is based on a primitive formulation of the set of equations which are not solved in a weak formulation. The overall scheme combines Finite Volume approaches for the discretization of the diffusion operator and Finite Difference approaches for the spatial discretization of convection terms in all equations. As a result, the proposed approach is not strictly conservative and is restricted to the computation of subsonic flows. Extensions of such pressure based compressible solvers to supersonic flows is not straightforward, and it will be considered in future works.

## 2. Fundamental equations for compressible flows and thermal plasma physical model

Fundamental equations, on which are based the pressure-based compressible solver, are presented in this section. The physical model for thermal plasma that will be coupled to our compressible solver is also described. Finally, in a last part, the boundary conditions to impose at a solid frontier are detailed for both the compressible flow and the plasma model, since it is a critical point for the design of an accurate IBM.

### 2.1. Fundamental equations for a pressure-based compressible flow

Fundamental equations are based on primary principles for mass, momentum and energy conservation, to which must be added some constitutive laws and Equations of State (EoS) to close the overall system. The system of equations is similar to the one presented in [2] for compressible two-phase flows. It is worth mentioning that it has been extended to configurations implying liquid-vapor phase change in [45]. However, the present study is only concerned with single phase flows (liquid or gas). Mass conservation is expressed as,

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{V}) = 0, \qquad (1) $$

where $\rho$ and $\vec{V}$ are respectively the density and velocity fields. Momentum conservation is ensured by solving the Navier-Stokes equations for a Newtonian fluid,
$$ \rho \frac{D\vec{V}}{Dt} = -\nabla p + \nabla \cdot \tau + \rho\vec{g}, \qquad (2) $$

with p the pressure, $\vec{g}$ the gravity acceleration and $\tau$ the tensor of viscous constraints defined as,

$$ \tau = \mu (\nabla \vec{V} + \nabla \vec{V}^T) - \frac{2}{3} \mu \nabla \cdot \vec{V} \mathbf{I}, \qquad (3) $$

with $\mu$ being the dynamic viscosity coefficient that can depend on temperature and I the identity matrix. This equation is almost identical to the incompressible formalism, except the additional term in the viscous tensor due to compression or expansion motion. The energy equation is expressed through the pressure variable as follows,

$$ \frac{Dp}{Dt} = -\rho c^2 \nabla \cdot \vec{V} + \frac{c^2 \alpha}{c_p} \left( -\nabla \cdot \vec{q} + \tau \otimes \nabla \vec{V} + \sum_{i=1}^{n} S_i \right), \qquad (4) $$

where c is the sound speed, $c_p$ the heat capacity at constant pressure, $\vec{q}$ the local heat flux defined as $\vec{q} = -k\nabla T$ from the Fourier law with k the thermal conductivity. $\sum_{i=1}^{n} S_i$ is the sum of volumetric source or sink heating terms that can account for heat release or loss due to different phenomena such as Joule effect, combustion, absorption or emission of thermal radiation, and $\alpha$ is the thermal expansion isobaric coefficient, defined as,

$$ \alpha = - \frac{1}{\rho} \frac{\partial \rho}{\partial T} \Big|_{\rho}. \qquad (5) $$

This form of the energy equation is typical of pressure-based solvers for compressible flows. This kind of solver has two benefits since it enables an implicit temporal discretization of acoustic waves, and it is Asymptotically Preserving (AP) of the incompressible regime when the Mach number tends to zero, as pointed out in several anterior works. More details on the derivation of this equation can be found in the following recent reference [2]. Although such a formulation is often associated to semi-implicit temporal resolution, a fully explicit resolution of this set of equations is also possible, as in [46].

## 2.2. Equations of states and thermodynamic relations

As our previous system of partial differential equations contains four unknowns variables $(\rho, \vec{V}, P, T)$ and only three conservation equations, an Equation of State (EoS) is required to close the system. This enables to deduce the temperature field from the pressure and density fields that are computed respectively from energy and mass conservation equations. The simplest EoS is the ideal gas equation,

$$ p = \rho r T \qquad (6) $$

with $r = R/M$, $R$ and $M$ being respectively the ideal gas constant and the molar mass. However, as this equation of state is only valid for gas at moderate pressure and relatively high temperature, more complex EoS, as cubic one, can be used to describe realistically both gases in higher pressure or lower temperature conditions and also liquids. Indeed, cubic EoS, as van der Waals or Peng-Robinson EoS, are able to predict a suitable thermodynamic behavior for the liquid state of many chemical species. In particular, when solving the cubic equation for the density at given pressure and temperature in subcritical conditions, three roots can be found: one corresponds to a dilute state (vapor state), another corresponds to a dense state (liquid state) and a third value, intermediate between gas and liquid, which has no physical meaning. Consequently, such equations are of significant interest to solve compressible flows in liquids for which realistic and thermodynamically consistent equations of states are required. As the overall IBM algorithm will require some extensions if a cubic EoS is used instead of an ideal gas EoS, we briefly remind here the general formalism associated to cubic equations of states. They are commonly expressed under the following generalized form,

$$ p = \frac{\rho RT}{1 - \rho b} - \frac{a \rho^2}{(1 + b_1 \rho)(1 + b_2 \rho)}, \qquad (7) $$

where $a, b, b_1$ and $b_2$ are properties of chemical species which are defined in Appendix A. They can also be expressed, as a cubic equation for the density $\rho$,

$$ \rho^3 + a_1 \rho^2 + a_2 \rho + a_3 = 0, \qquad (8) $$

where the $a_i$ coefficients depend on $T$, $p$, $a$, $b$, $b_1$ and $b_2$. More details on the different EoS can be found in Appendix A or in the following reference [45]. From the EoS, sound speed has to be defined explicitly, since it appears in the right-hand side of the pressure-based energy conservation equation. For an ideal gas EoS, it can be simply written,

$$ c = \sqrt{\gamma r T} \qquad (9) $$

where $\gamma = \frac{c_p}{c_v}$ is the ratio of specific heat capacities. For other equations of state, the following more general definition,

$$ c^2 = \left(\frac{\partial p}{\partial \rho}\right)_S = \left(\frac{\partial p}{\partial \rho}\right)_T + \frac{T}{\rho^2 c_v} \left(\frac{\partial p}{\partial T}\right)_\rho^2 \qquad (10) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. L.T.E. composition of an air thermal plasma at atmospheric pressure. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Mass density of an air thermal plasma at atmospheric pressure.</figcaption></figure>

is used to compute sound speed.

## 2.3. Thermal plasma model

Thermal plasma can be observed in various configurations, as the discharge of an electric arc in a high-voltage circuit breaker or plasma torch, for instance. Thermal plasma modeling is based on chemical equilibrium and the Local Thermal Equilibrium (LTE) assumption between charged ions and free electrons, that is, a single temperature can be defined for all particles: electrons and ions. From this assumption, it is then possible to calculate the composition of the plasma depending on temperature and pressure. This composition is calculated by using equilibrium laws like Saha equations, Dalton equation and neutrality. More details are given in [47]. An example of plasma composition is presented for air at atmospheric pressure in Fig. 1. At low temperature $O_2$ and $N_2$ are majority. Then when the temperature increases they dissociate to create atoms and at temperature over 10kK ionization occurs creating electrons and once ionized species.

From this composition, it is possible to determine thermodynamic properties where equilibrium chemistry is directly taken into account. For example, the mass density obtained from this composition is represented in Fig. 2. Two curves have been plotted. One corresponding to the mass density obtained from a “cold” perfect gas law, and the other corresponding to the one obtained from the composition presented in Fig. 1. We can see that in the case of thermal plasma, the mass density can vary by three order of magnitude on the temperature scale.

The last part concerns the transport coefficients of the plasma. They can be calculated by Chapman-Enskog formulation [47], and an example for the thermal conductivity is given in Fig. 3. In the plasma, the thermal conductivity can be calculated through three component (translation, chemical, internal). The peaks observed are due to the chemical part and can affect strongly the resolution of the system as they lead to a strong variation of the thermal conductivity. All thermodynamic and transport coefficient curves are given in Appendix B (Fig. B.30).

Consequently, the compressible solver presented in [2,45] have to be adapted to simulate the development of thermal plasma and related flows. However, the overall model must be upgraded to account for additional physical electric phenomena. In this framework, the IBM will be particularly convenient to describe the shape of the electrodes that generate the electric discharge. For this purpose, we have to impose the conservation of the electric current,

$$ \nabla \cdot \vec{J} = 0 \qquad (11) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Thermal conductivity of an air thermal plasma at atmospheric pressure.</figcaption></figure>

where $\vec{J}$ is the current density vector, which is expressed using the electric field,

$$ \vec{J} = \sigma \vec{E} \qquad (12) $$

with $\vec{E}$ the electric field and $\sigma$ the electric conductivity calculated from the plasma composition. It is worth mentioning that only electrically neutral plasma will be considered in this study, which results in a stationary equation for the electrical charge conservation without source or sink terms. The electric field can be defined as the sum of the gradient of an electrical potential $U$ and of the temporal derivative of a magnetic potential vector $\vec{A}$,

$$ \vec{E} = -\nabla U - \frac{\partial \vec{A}}{\partial t}. \qquad (13) $$

However, as magnetic effects are neglected in this paper, the following Laplace equation can be easily found for the electrical potential $U$,

$$ \nabla \cdot (\sigma \nabla U) = 0. \qquad (14) $$

The electric conductivity in a plasma is a function of the local temperature and pressure. It can be extracted from data banks. As already mentioned, the exact calculation of all the transport coefficients: electrical conductivity, thermal conductivity, viscosity inside a plasma can be computed with the Chapman-Enskog theory, as presented in [48,47] and their values depending on the temperature and pressure are then tabulated.

After computing the electric potential field $U$ and deducing the electric current $\vec{J}$, the Joule effect heating is imposed in the energy conservation Eq. (4) through the volumetric power $S_1$ such as,

$$ S_1 = \frac{\vec{J} \cdot \vec{J}}{\sigma}. \qquad (15) $$

This term is responsible for the high temperature during the electric arc discharge which leads to molecules dissociation, ionization, and so to the formation of a thermal plasma. Due to the very high temperature in the plasma medium, an additional heat sink term has to be imposed in Eq. (4) to take into account radiative losses. This can be expressed in a general way as,

$$ S_2 = -4\pi\epsilon, \qquad (16) $$

which corresponds to the divergence of the radiant flux $\vec{d}_{rad}$. The net emission coefficient $\epsilon$ can be computed by considering the part of the total radiation energy leaving an infinite cylinder or a sphere at uniform temperature and Local Thermodynamic Equilibrium, as detailed in [49]. It adds the radiation emitted through the outer shell of the cylinder or the sphere. Inside the plasma, the thermodynamic and transfer coefficients are variable because of changes in temperature and the dissociation of the particles. So the specific heats isobaric $c_p$ and isochoric $c_v$ are variable. The thermal conductivity $\lambda$, the dynamic viscosity coefficient $\mu$, the electric conductivity coefficient $\sigma$, the net emission coefficient $\epsilon$ are all dependent on the temperature and pressure with large relative variations between the extrema, for which quantitative data can be found in Appendix B. The temperature in the plasma state depends on the local density and pressure and is determined from tabulated values. Indeed, in Eq. (6), $r$, defined as $r = \frac{R}{M}$, is a constant for an ideal and single component gas with $R = 8.314 \text{ J.mol}^{-1}.\text{K}^{-1}$. However, in the case where the gas is a mixture of several chemical species, a temporally and spatially varying $r^n$ can be defined as $r^n = \frac{P^n}{\rho^n T^n}$, depending on the local variations of the chemical composition. For instance, that kind of model can be used to describe the density of plasmas, assumed to be a mixture of ideal gases whose composition varies with the temperature.
## 2.4. Immersed boundary condition on the solid frontier

IBM have been initially developed in the framework of incompressible flows in [3] and further developed in [6–10,50]. Although most pioneering methods were based on a smoothed forcing term at the solid/fluid interface, many alternative methods have been proposed thereafter to avoid this interface smoothing, by imposing sharp boundary conditions at the solid-fluid frontier. The work presented in [51] has been a precursor of such approaches, sometimes referred to as Immersed Interface Method (IIM), instead of IBM. Following this guideline, in [1], the authors present an efficient and accurate approach that maintains second order accuracy both on velocity and pressure, while preserving a sharp description of the solid/fluid frontier and a symmetric definite positive matrix resulting from the discretization of the pressure Poisson equation. This approach is based on a classical projection method for incompressible flows, in which a velocity splitting is carried out by introducing a non-solenoidal intermediate velocity $\vec{V}^*$ which does not account for pressure effect. Next, a Poisson equation has to be solved to update the pressure and a final velocity $\vec{V}^{n+1}$, which respects the divergence-free condition, can be computed with the following relation,

$$ \vec{V}^{n+1} = \vec{V}^* - \Delta t \frac{\nabla p}{\rho}. \qquad (17) $$

At a solid/fluid interface $\Gamma$, this will lead to the following scalar relations by projecting along the normal,

$$ V_{n,\Gamma}^{n+1} = V_{n,\Gamma}^* - \frac{\Delta t}{\rho} \left. \frac{\partial p}{\partial n} \right|_{\Gamma}, \qquad (18) $$

and tangential directions,

$$ V_{t,\Gamma}^{n+1} = V_{t,\Gamma}^* - \frac{\Delta t}{\rho} \left. \frac{\partial p}{\partial t} \right|_{\Gamma}, \qquad (19) $$

where $V_{n,\Gamma}$ and $V_{t,\Gamma}$ are respectively the normal and tangential velocity components at the interface $\Gamma$ and $\left. \frac{\partial p}{\partial n} \right|_{\Gamma} = \vec{n} \cdot \nabla p|_{\Gamma}$ and $\left. \frac{\partial p}{\partial t} \right|_{\Gamma} = \vec{t} \cdot \nabla p|_{\Gamma}$ are the pressure derivatives along the normal and tangent directions. Note that in a 3D configuration, two tangent velocity components at the interface have to be defined. At a static solid/fluid interface, the following boundary conditions,

$$ V_{n,\Gamma}^{n+1} = 0, \qquad (20) $$

$$ V_{t,\Gamma}^{n+1} = 0, \qquad (21) $$

have to be imposed to prescribe respectively the impermeability of the solid frontier and a no-slip condition which ensures the fluid adhesion at the solid frontier due to viscous friction. As $\vec{V}^*$ is an intermediate velocity which has no physical meaning, the boundary conditions $V_{n,\Gamma}^*$ and $V_{t,\Gamma}^*$ are not known. However, arbitrary values can be chosen provided Eq. (18) and Eq. (19) are correctly imposed. This means that the boundary condition $\left. \frac{\partial p}{\partial n} \right|_{\Gamma}$ that must be imposed depends on the value $V_{n,\Gamma}^*$ that has been chosen. For example, if one chooses

$$ V_{n,\Gamma}^* = 0, \qquad (22) $$

the corresponding boundary condition that respects Eq. (18) is,

$$ \left. \frac{\partial p}{\partial n} \right|_{\Gamma} = 0 \qquad (23) $$

In a same way, the pressure boundary condition related to,

$$ V_{t,\Gamma}^* = 0, \qquad (24) $$

is,

$$ \left. \frac{\partial p}{\partial t} \right|_{\Gamma} = 0 \qquad (25) $$

Such classical boundary conditions can be imposed with appropriate numerical schemes, as proposed in [52] for a Dirichlet boundary condition and in [1] for a Neumann boundary condition when discretizing a Laplacian at an immersed interface with a classical second order scheme. It is noteworthy that following the same principles, a Robin boundary condition can be also imposed with the method presented in [53] which has been used in [54,55] in the context of droplet evaporation. Our objective in this paper is to propose an extension of the IBM developed in [1] for incompressible flows to a pressure-based solver for compressible flows. Therefore, we take interest, now, to additional boundary conditions that must be imposed to properly account for a solid/fluid interface in cells crossed by a solid frontier in the framework of compressible flows. In particular, as density and temperature fields are computed with the compressible solver, some specific care must be devoted to the computation of density and temperature gradients in the vicinity of the Immersed Interface. For the density field, the following boundary condition has to be imposed,

$$ \left. \frac{\partial p}{\partial n} \right|_{\Gamma} = 0. \qquad (26) $$
For the temperature field, different kinds of boundary conditions can be considered by the present IBM, as an immersed Dirichlet boundary condition on the temperature that simply writes,

$$T_{\Gamma} = T_{wall} \qquad (27)$$

or an Immersed Neumann boundary condition,

$$k \left. \frac{\partial T}{\partial n} \right|_{\Gamma} = q_{wall}, \qquad (28)$$

where variables $T_{wall}$ and $q_{wall}$ can be homogeneous or can vary spatially and/or temporally. It is noteworthy, that the latter condition enables to impose a prescribed local heat flux and becomes an adiabatic condition if $q_{wall} = 0$. We will restrict our analysis to these two kinds of boundary conditions for the sake of simplicity. However, a more general formulation could be considered to take into account conjugate heat transfer between the fluid and the solid domains. Such a formulation would require to solve the heat equation in the solid domain and to couple it with the energy equation in the fluid domain through the following heat flux continuity equation,

$$k_f \left. \frac{\partial T}{\partial n} \right|_{\Gamma,f} = k_s \left. \frac{\partial T}{\partial n} \right|_{\Gamma,s}, \qquad (29)$$

where the subscript f and s refer respectively to the fluid and solid domains. Examples of such a coupling are presented in [56,57] in the framework of nucleate boiling simulations on a plane solid frontier with an incompressible solver. The extension of the present solver to configurations where conjugate heat transfer is taken into account will be the object of future works.

As some extensions related to thermal plasma modeling are also presented in this paper, the boundary conditions related to our thermal plasma model are now presented. When solving Eq. (14), boundary conditions have to be imposed on electrodes to generate an electric potential difference to induce an electric current, and eventually the formation of a plasma. This can be carried out as a Dirichlet condition representing a fixed potential,

$$V|_{\Gamma} = V_{electrode}, \qquad (30)$$

on an electrode. A zero potential will be imposed on the ground. On the other hand, a Neumann boundary condition can be also imposed, if one aims to fix a prescribed current density $J_{n,0}$ in the normal direction at the interface,

$$\sigma \left. \frac{\partial V}{\partial n} \right|_{\Gamma} = J_{n,0}. \qquad (31)$$

Note that $J_{n,0} = 0$ corresponds to an electrically insulated frontier.

## 3. Numerical methods

We present in this section the numerical methods implemented for the coupling between the IBM proposed in [1] and the pressure-based compressible solver recently published in [2]. The overall pressure-based solver is presented in a first subsection. In a second subsection, the spatial schemes that takes into account the solid frontier are presented for cells containing the interface, both for an Immersed Neumann Boundary Condition (INBC) and an Immersed Dirichlet Boundary Condition (IDBC). In the framework of this study, these fundamental schemes are useful for discretizing several terms of the overall system, both for fluids and plasma simulations. Classical extensions processes which are required to populate ghost cells for some fields are then briefly presented. The Finite Volume discretization with INBC and/or IDBC of the mass, momentum and energy conservation equations are detailed in the following sections. Next, description of numerical methods for thermal plasma computations is also addressed. Finally, this section is concluded by some remarks on generic numerical methods and restrictions on the time step.

### 3.1. Overall algorithm of the pressure-based compressible solver

The algorithm resolution of the pressure-based solver for compressible flows proposed in [2] is detailed in this subsection. As, in the framework of IBM, implicit temporal discretization may be required for some terms, since interpolation schemes based on a subcell resolution are employed to impose a given value of a variable at the intersection between the solid interface and a segment of the mesh. Such subcell interpolations may lead to unstable schemes if associated to an explicit temporal scheme. This is especially true for viscous and heat conduction terms, and this is why implicit temporal schemes may be required to impose some Immersed Boundary Conditions.

#### 3.1.1. Explicit discretization of the viscous and heat conduction terms

We remind in this section the algorithm on which is based the pressure-based compressible solver presented in [2,45]. For the sake of simplicity, the overall algorithm is presented with a first order temporal scheme, but can be generalized to higher order as Runge-Kutta schemes. The density field is computed from the mass conservation equation following an explicit discretization,

$$\rho^{n+1} = \rho^n - \Delta t \nabla \cdot (\rho \vec{V})^n. \qquad (32)$$

The temporal discretization of the pressure equation is carried out with a semi-implicit temporal scheme, in which the divergence velocity, $\nabla \cdot \vec{V}$, in the term related to the propagation of the acoustic waves, is taken at time $t^{n+1}$,
$$ \frac{p^{n+1} - p^n}{\Delta t} + \vec{V}^n \cdot \nabla p^n + (\rho c^2)^n \nabla \cdot \vec{V}^{n+1} = \left( \frac{c^2 \alpha}{c_\rho} \right)^n (\tau^n \otimes \nabla \vec{V}^n - \nabla \cdot \vec{q}^n + S), \quad (33) $$

in order to remove the acoustic time step constraint. From the following velocity splitting,

$$ \vec{V}^{n+1} = \vec{V}^* - \Delta t \frac{\nabla p^{n+1}}{\rho^{n+1}}, \quad (34) $$

typical of projection methods for incompressible solvers, we can define the intermediate velocity $V^*$,

$$ \vec{V}^* = \vec{V}^n - \Delta t \left( \vec{V}^n \cdot \nabla \vec{V}^n - \frac{\nabla \cdot \tau^n}{\rho^{n+1}} - \vec{g} \right). \quad (35) $$

An intermediate pressure $p^*$ is also defined to account for the convection of the pressure field,

$$ p^* = p^n - \Delta t \vec{V}^n \cdot \nabla p^n. \quad (36) $$

From the previous two steps, the following pressure equation is obtained,

$$ \frac{p^{n+1}}{(\rho c^2)^n} - \Delta t^2 \nabla \cdot \left( \frac{\nabla p^{n+1}}{\rho^{n+1}} \right) = \frac{p^*}{(\rho c^2)^n} + \Delta t \nabla \cdot \vec{V}^* + \Delta t \left( \frac{\alpha}{\rho c_\rho} \right)^n \left( \tau^n \otimes \nabla \vec{V}^n - \nabla \cdot (k \nabla T^n) + S \right). \quad (37) $$

After computing, $p^{n+1}$, the velocity $V^{n+1}$ can be easily deduced from Eq. (34). Finally, the last step of the algorithm is to compute the temperature from $p^{n+1}$ and $\rho^{n+1}$ by means of the EOS,

$$ T^{n+1} = \frac{p^{n+1}}{\rho^{n+1} r}, \quad (38) $$

for an ideal gas law. In the case of a cubic EoS, as van der Waals or others, the temperature can be simply obtained with the following relation,

$$ T^{n+1} = \frac{1}{\rho^{n+1} R} (M - \rho b) \left( P + \frac{a \alpha \rho^2}{M^2 + 2bM\rho - b^2\rho^2} \right), \quad (39) $$

which will be useful to model a liquid phase or a gas in high pressure and moderate temperature conditions.

### 3.1.2. Implicit discretization of viscous and heat conduction terms

It is worth mentioning that the temporal discretizations of the viscous term in Eq. (35) and of the heat conduction term in Eq. (37) are explicit, and thus implies a time step constraint related respectively to viscosity and heat conduction, that is $\Delta t_{visc} \sim \frac{\rho \Delta \chi^2}{\mu}$ and $\Delta t_{heat} \sim \frac{\rho C_p \Delta \chi^2}{k}$. However, in the framework of the IBM proposed in [1], an implicit discretization of the viscous and heat conduction terms can be required if Dirichlet boundary condition has to be imposed on the solid frontier for the velocity or the temperature. Indeed, as shown in [52], such an immersed boundary condition is imposed through a subcell resolution in which appears a submesh $\theta h$, with $\theta$ ranging between [0, 1] and $h$ the length of a grid cell. This subcell resolution induces more drastic temporal restrictions, i.e. $\Delta t_{visc} \sim \frac{\theta \rho \Delta \chi^2}{\mu}$ and $\Delta t_{heat} \sim \frac{\theta \rho C_p \Delta \chi^2}{k}$, for viscosity and heat conduction respectively. In order to alleviate this time step constraint, an implicit temporal discretization can be mandatory to maintain a reasonably quick time advancement of the computed solution. With respect to the time step constraint $\Delta t_{visc}$, the following implicit temporal scheme can be applied,

$$ \vec{V}^* - \Delta t \frac{1}{\rho^{n+1}} \nabla \cdot \tau^* = \vec{V}^n - \Delta t (\vec{V}^n \cdot \nabla \vec{V}^n - \vec{g}), \quad (40) $$

as proposed in [58] for instance. This scheme involves a coupled linear system between all the velocity components. As the resulting linear system is strongly diagonally dominant in most cases (depending on the constraint on the time step), a simple Gauss-Seidel method is generally efficient to solve this coupled linear system. By applying this implicit scheme, the temporal constraint related to viscosity is removed. If we take interest now to the temporal constraint related to the heat conduction term, an implicit scheme has been proposed in the following references [2,45] to improve the stability computations with regards to $\Delta t_{heat}$. This can be carried out by taking $T$ at time $n+1$ in the pressure equation, such as,

$$ \frac{p^{n+1}}{(\rho c^2)^n} - \Delta t^2 \nabla \cdot \left( \frac{\nabla p^{n+1}}{\rho^{n+1}} \right) = \frac{p^*}{(\rho c^2)^n} + \Delta t \nabla \cdot \vec{V}^* + \Delta t \left( \frac{\alpha}{\rho c_p} \right)^n \left( \tau^n \otimes \nabla \vec{V}^n - \nabla \cdot (k \nabla T^{n+1}) + S \right) \quad (41) $$

One can remark that for any cubic EoS, the following linear general expression between temperature and pressure can be found,

$$ T^{n+1} = C^{n+1} \rho^{n+1} + D^{n+1}, \quad (42) $$

with,

$$ C^{n+1} = \frac{1 - b \rho^{n+1}}{\rho^{n+1} r} \quad \text{and} \quad D^{n+1} = \frac{a \rho^{n+1} (1 - b \rho^{n+1})}{r(1 + b_1 \rho^{n+1})(1 + b_2 \rho^{n+1})}, \quad (43) $$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Schematic of the Immersed Boundary Method on a structured grid.</figcaption></figure>

with the coefficients $a, b, b_1, b_2$ defined in the Appendix A. The following simplified relations, $C^{n+1} = \frac{1}{\rho^{n+1}r}$ and $D^{n+1} = 0$, is obtained in the simpler case of an ideal gas. Since our energy conservation equation is pressure-based, Eq. (42) can be injected to express $T^{n+1}$ with $p^{n+1}$ in Eq. (41), which leads to the following temporal discretization of the pressure equation,

$$ \frac{p^{n+1}}{(\rho c^2)^n} - \Delta t^2 \nabla \cdot \left( \frac{\nabla p^{n+1}}{\rho^{n+1}} \right) - \Delta t \left( \frac{\alpha}{\rho^n c_p} \right)^n \nabla \cdot \left( k \nabla \left( C^{n+1} p^{n+1} \right) \right) = \frac{p^*}{(\rho c^2)^n} + \Delta t \nabla \cdot V^* + \Delta t \left( \frac{\alpha}{\rho c_p} \right)^n \left( \tau^n \otimes \nabla \vec{V}^n + \nabla \cdot (k \nabla (D^{n+1})) + S \right), \quad (44) $$

in which an additional term appears in the left-hand side, and thus in the pressure matrix. Whereas this additional term presents many similarities with the acoustic term and does not change the number of diagonals in the resulting matrix, this matrix does not remain symmetric due to the inhomogeneous pre-factor which multiplies this term. More details on the matrix coefficients originating from this implicit discretization of the heat conduction term can be found in [45]. However, as pointed out in the latter reference, this implicit scheme does not fully remove the time step constraint $\Delta t_{heat}$ since Eq. (41) is not linear with regards to the heat conduction.

Indeed, the pre-factor $\left( \frac{\alpha}{\rho c_p} \right)^n$ that is multiplying the temperature Laplacian depends also on temperature, which makes non-linear the heat conduction term in a pressure-based formulation. Therefore, a fully implicit resolution of this term would require a solver for non-linear systems involving sub-cycles iterations. Although the algorithm presented in this section is not fully implicit, it permits, however, to alleviate the time step constraint $\Delta t_{heat}$ by multiplying it by a factor 3 with regards to the explicit stability constraints, as reported in [45]. This can be of significant interest in configurations for which an immersed boundary condition on the temperature has to be imposed.

## 3.2. Immersed Boundary Condition following Gibou and co-workers methods

The development of numerical schemes to impose Immersed Boundary Condition inside a structured mesh when solving elliptic equations was subject to many past works, among which pioneering works by Johansen and Colella [59], and further improvements in [60,61]. In a series of papers, Gibou and co-workers have developed second order schemes for an Immersed Dirichlet Boundary Condition (IDBC) in [52], an Immersed Neumann Boundary Condition (INBC) in [1] and an Immersed Robin Boundary Condition (IRBC) in [53]. The latter two schemes are based on a finite volume approach accounting for the interface position inside the cells crossed by the solid frontier, which is briefly reminded in the subsection hereafter. In the coming subsections, all discretization schemes are presented in a 2D configuration for the sake of simplification, but generalizations to 3D configurations are straightforward.

### 3.2.1. Immersed Neumann Boundary Condition scheme (INBC scheme)

For the sake of illustration, a schematic of the IBM on a structured grid and a computational element crossed by a solid frontier are represented in Fig. 4. The part of the cell corresponding to the field which is solved is the white region, which represents the control volume $\Omega_{ij}$ over which integration is performed. The surface delimiting this control volume is denoted $\Sigma_{ij}$, and it is decomposed in five segments $L_{i+1/2j}$, $L_{i-1/2}$, $L_{ij+1/2}$, $L_{ij-1/2}$ and $L_\Gamma$ represented in Fig. 4. The lengths (for a 2D computation) or surfaces (for a 3D computation) delimiting the computational element can be computed from geometric considerations by defining the intersection points of the solid interface with the boundaries of the computational cell. After defining a static Level Set function $\phi$ delimiting the solid frontier (with $\phi > 0$ in the solid and $\phi < 0$ in the fluid), the following procedure can be applied to compute $L_{i-1/2j}$ for example,
$$L_{i-\frac{1}{2},j} = \begin{cases} \Delta y \frac{\phi_{i-\frac{1}{2},j-\frac{1}{2}} - \phi_{i-\frac{1}{2},j+\frac{1}{2}}}{\phi_{i-\frac{1}{2},j-\frac{1}{2}} - \phi_{i-\frac{1}{2},j+\frac{1}{2}}} & \text{if } \phi_{i-\frac{1}{2},j-\frac{1}{2}} < 0 \text{ and } \phi_{i-\frac{1}{2},j+\frac{1}{2}} > 0 \\ \Delta y \frac{\phi_{i-\frac{1}{2},j+\frac{1}{2}} - \phi_{i-\frac{1}{2},j-\frac{1}{2}}}{\phi_{i-\frac{1}{2},j+\frac{1}{2}} - \phi_{i-\frac{1}{2},j-\frac{1}{2}}} & \text{if } \phi_{i-\frac{1}{2},j-\frac{1}{2}} > 0 \text{ and } \phi_{i-\frac{1}{2},j+\frac{1}{2}} < 0 \\ \Delta y & \text{if } \phi_{i-\frac{1}{2},j-\frac{1}{2}} < 0 \text{ and } \phi_{i-\frac{1}{2},j+\frac{1}{2}} < 0 \\ 0 & \text{if } \phi_{i-\frac{1}{2},j-\frac{1}{2}} > 0 \text{ and } \phi_{i-\frac{1}{2},j+\frac{1}{2}} > 0 \end{cases} \quad (45)$$

In a 3D configuration, surfaces must be computed instead of lengths. This more difficult task requires a suitable interface reconstruction method to compute accurately $L_Γ$, as proposed in [62] and [63]. Several examples and validations of 3D computation based on this approach for incompressible flows can be found in [1,58,64,65]. After defining computational cells modified by the immersed solid, divergence and gradient operators in the fundamental equations are discretized following the Finite Volume Method (FVM). For example, the volume integral of the divergence of a given vector $\vec{\phi}$ over a volume cell $\Omega_{ij}$, delimited by the surface $\Sigma_{ij}$, will be transformed as a surface integral following the Green-Ostrogradski theorem, such as,

$$\iiint_{\Omega_{ij}} \nabla \cdot \vec{\phi} \, d\Omega_{ij} = \iiint_{\Sigma_{ij}} \vec{\phi} \cdot \vec{n}_{\Sigma} \, d\Sigma \quad (46)$$

where $d\Sigma$ is a surface element of $\Sigma_{ij}$, and $\vec{n}_{\Sigma}$ the normal vector pointing outside the surface element. If applied to a cut cell as in Fig. 4, the following discrete scheme is obtained,

$$\iiint_{\Sigma_{ij}} \vec{\phi} \cdot \vec{n}_{\Sigma} \, d\Sigma = L_{i+1/2j} (\vec{\phi} \cdot \vec{n}) |_{i+1/2j} + L_{i-1/2j} (\vec{\phi} \cdot \vec{n}) |_{i-1/2j} + L_{ij+1/2} (\vec{\phi} \cdot \vec{n}) |_{ij+1/2} \\ + L_{ij-1/2} (\vec{\phi} \cdot \vec{n}) |_{ij-1/2} + L_{\Gamma} (\vec{\phi} \cdot \vec{n}) |_{\Gamma} \quad (47)$$

which finally gives,

$$\iiint_{\Sigma_{ij}} \vec{\phi} \cdot \vec{n}_{\Sigma} \, d\Sigma = L_{i+1/2j} \varphi_{x,i+1/2j} - L_{i-1/2j} \varphi_{x,i-1/2j} + L_{ij+1/2} \varphi_{y,ij+1/2} - L_{ij-1/2} \varphi_{y,ij-1/2} + L_{\Gamma} (\vec{\phi} \cdot \vec{n}) |_{\Gamma}, \quad (48)$$

where $\varphi_x$ and $\varphi_y$ are respectively the x and y components of $\vec{\phi}$. The Laplacian operator $\nabla \cdot (D \nabla u)$ of a given function $u$ with $D$ a diffusion coefficient is now considered. The following INBC is imposed at the solid interface,

$$D_{\Gamma} \left. \frac{\partial u}{\partial n} \right|_{\Gamma} = b_{\Gamma}, \quad (49)$$

where $\partial \left. \frac{\partial}{\partial n} \right|_{\Gamma}$ is the partial derivative along the normal direction at the interface. After applying the Green-Ostrogradski theorem,

$$\iiint_{\Omega_{ij}} \nabla \cdot (D \nabla u) \, d\Omega_{ij} = \iiint_{\Sigma_{ij}} D \nabla u \cdot \vec{n}_{\Sigma} \, d\Sigma, \quad (50)$$

the following FV scheme can be developed,

$$\iiint_{\Sigma_{ij}} D \nabla u \cdot \vec{n}_{\Sigma} \, d\Sigma = L_{i+1/2j} D_{i+1/2j} \left. \frac{\partial u}{\partial x} \right|_{i+1/2j} - L_{i-1/2j} D_{i-1/2j} \left. \frac{\partial u}{\partial x} \right|_{i-1/2j} \\ + L_{ij+1/2} D_{ij+1/2} \left. \frac{\partial u}{\partial y} \right|_{ij+1/2} - L_{ij-1/2} D_{ij-1/2} \left. \frac{\partial u}{\partial y} \right|_{ij-1/2} + L_{\Gamma} b_{\Gamma}. \quad (51)$$

This leads finally to the following discretization,

$$\iiint_{\Sigma_{ij}} D \nabla u \cdot \vec{n}_{\Sigma} \, d\Sigma = L_{i+1/2j} D_{i+1/2j} \left( \frac{u_{i+1,j} - u_{i,j}}{\Delta x} \right) + L_{i-1/2j} D_{i-1/2j} \left( \frac{u_{i,j} - u_{i-1,j}}{\Delta x} \right) \\ + L_{ij+1/2} D_{ij+1/2} \left( \frac{u_{i,j+1} - u_{i,j}}{\Delta y} \right) + L_{ij-1/2} D_{ij-1/2} \left( \frac{u_{i,j} - u_{i,j-1}}{\Delta y} \right) + L_{\Gamma} b_{\Gamma}. \quad (52)$$

This numerical scheme has attractive features, as maintaining a sharp representation of the solid interface, a second order spatial accuracy and the linear system resulting from this discrete scheme is symmetric definite positive. It can be employed for the discretization of the pressure term or the heat conduction term in Eq. (37) and Eq. (41) or also to compute the electric potential in Eq. (14), in cases involving a Neumann boundary condition at the solid frontier, following Eq. (31). Moreover, it can be directly generalized to IRBC as shown in [53].
### 3.2.2. Immersed Dirichlet Boundary Condition scheme (IDBC scheme)

Earlier, in [52], Gibou and Fedkiw proposed a simple but efficient numerical scheme to impose an Immersed Dirichlet Boundary Condition (IDBC). For the sake of consistency with the previous section, this numerical method is presented in the framework of the Finite Volume method. A Laplacian operator $\nabla \cdot (D\nabla u)$, to which is associated an IDBC at the interface $u_\Gamma = a_\Gamma$, is considered. Applying the general guidelines of the Finite Volume Method (FVM), we remind the following general relation,

$$ \iiint_{\Omega_{ij}} \nabla \cdot (D\nabla u) d\Omega_{ij} = \iint_{\Sigma_{ij}} D \nabla u \cdot \vec{n}_{\Sigma} d\Sigma \quad (53) $$

It is, now, supposed that the segment $[i, j : i + 1, j]$ is crossed by the interface $\Gamma$. By interpolating the first derivative in this cell $\Gamma$ as $D_x u|_{\Gamma} = \left( \frac{u_\Gamma - u_{ij}}{\theta \Delta x} \right)$, we obtain the following discretization scheme for the Laplacian operator at the point $(i, j)$,

$$ \iint_{\Sigma_{ij}} D \nabla u \cdot \vec{n}_{\Sigma} d\Sigma = D_{i+1/2j} \left( \frac{a_{\Gamma} - u_{ij}}{\theta \Delta x} \right) \Delta y - D_{i-1/2j} \left( \frac{u_{ij} - u_{i-1j}}{\Delta x} \right) \Delta y + D_{ij+1/2} \left( \frac{u_{ij+1} - u_{ij}}{\Delta y} \right) \Delta x \\ - D_{ij-1/2} \left( \frac{u_{ij} - u_{ij-1}}{\Delta y} \right) \Delta x, \quad (54) $$

with $\theta \Delta x$ an approximation of the distance between the intersection point of the interface with the segment $[i, j : i + 1, j]$. $\theta$ is defined from the level-set function $\phi$ as follows,

$$ \theta = \frac{|\phi_{ij}|}{|\phi_{ij}| + |\phi_{i+1,j}|}. \quad (55) $$

Configurations involving an interface crossing segments $[i - 1, j : i, j]$, $[i, j - 1 : i, j]$ or $[i, j : i, j + 1]$ or any combinations of such configurations can be approximated similarly. As for the INBC scheme presented in subsection 3.2.1, this scheme maintains a sharp discretization of the interface, allows a second order accuracy in space and the linear system matrix resulting from this scheme conserves symmetric and definite positive properties. It is noteworthy, that interface segments crossed by the interface involves a mesh cell $\sim \theta \Delta x$ with $\theta$ taking any values between 0 and 1. This involves a drastic constraint on the time step if associated to an explicit temporal discretization. This is why this numerical scheme is generally associated to an implicit temporal discretization for which the time step constraint is removed, [52,58]. In some cases, where a fully implicit temporal discretization can be more difficult to implement, see [45] for example, a limiter on $\theta$, such as $\theta \ge \theta_{min}$ with $\theta_{min} = 0.1$ can be employed without significant loss of accuracy, to avoid a too drastic time step constraint.

### 3.3. Finite Volume discretization of the mass, momentum and energy conservation equations

In this subsection, details on the discretization of the conservation equations are presented and numerical schemes related to the IBM are given. In particular, we specify the terms that require the INBC or IDBC numerical schemes. It is noteworthy that advection terms in equations are not formulated in a conservative form, and the resulting scheme is not a Finite Volume scheme for these terms, but a Finite Difference scheme instead. The reason of this choice lies in the energy equation formulated as a pressure equation which is not a conservative formulation.

#### 3.3.1. Finite Volume discretization of the mass conservation equation

From the temporal discretization of the mass conservation equation, Eq. (32), the following formulation is obtained,

$$ \iiint_{\Omega_{ij}} \rho_{ij}^{n+1} d\Omega_{ij} = \iiint_{\Omega_{ij}} \rho_{ij}^{n} d\Omega_{ij} - \Delta t \left( \iiint_{\Omega_{ij}} \vec{V} \cdot \nabla \rho \right)_{ij}^{n} d\Omega_{ij} + \iiint_{\Omega_{ij}} \rho_{ij}^{n} \nabla \cdot \vec{V}_{ij}^{n} d\Omega_{ij} \quad (56) $$

Assuming the density $\rho_{ij}^n$ is uniform inside a computational cell $\Omega_{ij}$ and applying Green-Ostrogradsky theorem, it becomes,

$$ \iiint_{\Omega_{ij}} \rho_{ij}^{n+1} d\Omega_{ij} = \iiint_{\Omega_{ij}} \rho_{ij}^{n} d\Omega_{ij} - \Delta t \left( \iiint_{\Omega_{ij}} \vec{V} \cdot \nabla \rho \right)_{ij}^{n} d\Omega_{ij} + \rho_{ij}^{n} \iint_{\Sigma_{ij}} \vec{V}_{ij}^{n} \cdot \vec{n}_{\Sigma} d\Sigma \quad (57) $$

Finally, by applying the INBC scheme given by Eq. (48), to the surface integral in the right-hand side of Eq. (57), the following equation is obtained,

$$ \rho_{ij}^{n+1} = \rho_{ij}^{n} - \Delta t \left( \vec{V} \cdot \nabla \rho \right)_{ij}^{n} - \frac{\rho_{ij}^{n} \left( L_{i+1/2j} u_{i+1/2j} - L_{i-1/2j} u_{i-1/2j} + L_{ij+1/2} v_{ij+1/2} - L_{ij-1/2} v_{ij-1/2} + L_{\Gamma} (\vec{V} \cdot \vec{n}) |_{\Gamma} \right)}{\Omega_{ij}} \quad (58) $$
in which the last term $(\vec{V} \cdot \vec{n})|_{\Gamma} = 0$ imposes the impermeability condition on the solid interface, as specified by Eq. (20).

### 3.3.2. Discretization of the momentum conservation equation

The temporal discretization of the predictor step of the momentum equation is given by Eq. (35) or Eq. (40) whether an explicit or an implicit temporal scheme is considered for the viscous term. However, as IDBC, given by Eq. (22) and Eq. (24), have to be imposed at the solid frontier for all velocity components, the implicit temporal scheme given by Eq. (40) must be employed for stability reasons. It can be formulated with Finite Volume as,

$$ \iiint_{\Omega} (\vec{V}^* - \frac{\Delta t}{\rho^{n+1}} \nabla \cdot \boldsymbol{\tau}^*) d\Omega = \iiint_{\Omega} (\vec{V}^n - \Delta t(\vec{V}^n \cdot \nabla \vec{V}^n - \vec{g})) d\Omega. \quad (59) $$

This leads to a coupled linear system between the velocity components, $u_{i+1/2j}$ and $v_{ij+1/2}$ for which all details can be found in [58] for two-phase incompressible flows with an IBM. It is worth mentioning that a staggered grid is used, which implies that velocity components are not expressed on the same computational grids. That is why generic space index can not be specified in Eq. (59).

### 3.3.3. Finite Volume discretization of the energy conservation equation with an explicit heat flux conduction

We take interest now in the pressure-based energy equation, which is a crucial part of the proposed algorithm. A formulation of the semi-discrete equation with an explicit temporal discretization of the heat conduction term, Eq. (37), based on FVM, gives,

$$ \iiint_{\Omega_{ij}} \frac{p^{n+1}}{(\rho c^2)^n} \Big|_{ij} d\Omega_{ij} - \Delta t^2 \iiint_{\Omega_{ij}} \nabla \cdot \left( \frac{\nabla p^{n+1}}{\rho^{n+1}} \right) \Big|_{ij} d\Omega_{ij} = \iiint_{\Omega_{ij}} \frac{p^*}{(\rho c^2)^n} \Big|_{ij} d\Omega_{ij} + \Delta t \iiint_{\Omega_{ij}} \nabla \cdot \vec{V}^* \Big|_{ij} d\Omega_{ij} + \Delta t \iiint_{\Omega_{ij}} \left( \left( \frac{\alpha}{\rho c_p} \right)^n \left( \tau^n \otimes \nabla \vec{V}^n - \nabla \cdot (k \nabla T^n) + S \right) \right) \Big|_{ij} d\Omega_{ij} \quad (60) $$

Remarking that Green-Ostrogradsky theorem can be applied to several terms of this equation, the following linear system on the pressure is obtained,

$$ \begin{aligned} \frac{p^{n+1}}{(\rho c^2)^n} \Big|_{ij} \Omega_{ij} - \Delta t^2 \iint_{\Sigma_{ij}} \left( \frac{\nabla p^{n+1}}{\rho^{n+1}} \right) \cdot \vec{n}_{\Sigma} d\Sigma = & \frac{p^*}{(\rho c^2)^n} \Big|_{ij} \Omega_{ij} + \Delta t \iint_{\Sigma_{ij}} \vec{V}^* \cdot \vec{n}_{\Sigma} d\Sigma \\ & - \Delta t \left( \frac{\alpha}{\rho c_p} \right) \Big|_{ij}^n \iint_{\Sigma_{ij}} k \nabla T^n \cdot \vec{n}_{\Sigma} d\Sigma \\ & + \Delta t \left( \left( \frac{\alpha}{\rho c_p} \right)^n \left( \tau^n \otimes \nabla \vec{V}^n + S \right) \right) \Big|_{ij} \Omega_{ij}. \end{aligned} \quad (61) $$

As the INBC, given by Eq. (23), has to be imposed on the pressure to satisfy the impermeability condition, the surface integral of the pressure gradient in the latter equation is expressed as,

$$ \begin{aligned} \iint_{\Sigma_{ij}} \left( \frac{\nabla p^{n+1}}{\rho^n} \right) \cdot \vec{n}_{\Sigma} d\Sigma = & \frac{L_{i+1/2j}}{\rho_{i+1/2j}} \left( \frac{p_{i+1j}^{n+1} - p_{ij}^{n+1}}{\Delta x} \right) - \frac{L_{i-1/2j}}{\rho_{i-1/2j}} \left( \frac{p_{ij}^{n+1} - p_{i-1j}^{n+1}}{\Delta x} \right) + \frac{L_{ij+1/2}}{\rho_{ij+1/2}} \left( \frac{p_{ij+1}^{n+1} - p_{ij}^{n+1}}{\Delta y} \right) \\ & - \frac{L_{ij-1/2}}{\rho_{ij-1/2}} \left( \frac{p_{ij}^{n+1} - p_{ij-1}^{n+1}}{\Delta y} \right) + \left( \frac{L_{\Gamma}}{\rho} \frac{\partial p}{\partial n} \right) \Big|_{\Gamma}^{n+1}, \end{aligned} \quad (62) $$

in which last term $\left( \frac{L_{\Gamma}}{\rho} \frac{\partial p}{\partial n} \right)_{\Gamma}^{n+1} = 0$, following Eq. (23). Similarly to Eq. (48), the following term is discretized as,

$$ \iint_{\Sigma_{ij}} \vec{V}^* \cdot \vec{n}_{\Sigma} d\Sigma = L_{i+1/2j} u_{i+1/2j}^* - L_{i-1/2j} u_{i-1/2j}^* + L_{ij+1/2} v_{ij+1/2}^* - L_{ij-1/2} v_{ij-1/2}^* + L_{\Gamma} (\vec{V}^* \cdot \vec{n})|_{\Gamma}, \quad (63) $$

with $L_{\Gamma}(\vec{V}^* \cdot \vec{n})|_{\Gamma} = 0$ according to Eq. (22). Finally, the heat conduction term in Eq. (61), can also be formulated as follows,

$$ \begin{aligned} \iint_{\partial \Omega_{ij}} k \nabla T^n \cdot \vec{n}_{\Sigma} d\Sigma = & L_{i+1/2j} k_{i+1/2j} \left( \frac{T_{i+1j}^n - T_{ij}^n}{\Delta x} \right) - L_{i-1/2j} k_{i-1/2j} \left( \frac{T_{ij}^n - T_{i-1j}^n}{\Delta x} \right) + L_{ij+1/2} k_{ij+1/2} \left( \frac{T_{ij+1}^n - T_{ij}^n}{\Delta y} \right) \\ & - L_{ij-1/2} k_{ij-1/2} \left( \frac{T_{ij}^n - T_{ij-1}^n}{\Delta y} \right) + L_{\Gamma} \left( k_{\Gamma} \frac{\partial T}{\partial n} \right) \Big|_{\Gamma}^n, \end{aligned} \quad (64) $$
if an INBC, as Eq. (28), has to be imposed at the solid frontier. On the other hand, if one wants to impose a prescribed temperature on the immersed boundary, Eq. (27), the numerical scheme for IDBC, given by Eq. (54) will be employed. It should be reminded here that such a numerical scheme employed with an explicit temporal scheme will require to define a limiter on the variable $\theta$ in Eq. (54), in order to avoid a too drastic time step constraint related to the heat conduction. Another noteworthy point is that such a limiter is not required with the INBC scheme, for which the classical time step constraint is sufficient. Indeed, unlike the IDBC scheme, the INBC scheme does not define a locally truncated cell size $\theta \Delta x$, and its overall temporal stability is not altered by the Immersed Boundary.

### 3.3.4. Finite Volume discretization of the energy conservation equation with an implicit heat conduction

An implicit temporal discretization of the heat conduction is now considered following Eq. (41). We consider for a first time an INBC, Eq. (28), which provides a relation similar to Eq. (64) but with $T^{n+1}$ instead of $T^n$. Therefore, reminding the following linear relation $T^{n+1} = C^{n+1}p^{n+1} + D^{n+1}$ with $C^{n+1}$ and $D^{n+1}$ given by Eq. (43), the heat conduction term can be expressed with the pressure variable $p^{n+1}$, and splits in two parts. A first part depending on $p^{n+1}$ that will be considered in the matrix of the linear system resulting from Eq. (41) and a second part that depends on $D^{n+1}$ that will be accounted for in the right-hand side of the linear system,

$$ \begin{aligned} \iint_{\partial \Omega_{ij}} k \nabla T^{n+1} \cdot \vec{n} \, d\Sigma = & L_{i+1/2j} k_{i+1/2j} \left( \frac{(C^{n+1} p^{n+1})|_{i+1j} - (C^{n+1} p^{n+1})|_{ij}}{\Delta x} + \frac{D_{i+1j}^{n+1} - D_{ij}^{n+1}}{\Delta x} \right) \\ & - L_{i-1/2j} k_{i-1/2j} \left( \frac{(C^{n+1} p^{n+1})|_{ij} - (C^{n+1} p^{n+1})|_{i-1j}}{\Delta x} + \frac{D_{ij}^{n+1} - D_{i-1j}^{n+1}}{\Delta x} \right) \\ & + L_{ij+1/2} k_{ij+1/2} \left( \frac{(C^{n+1} p^{n+1})|_{ij+1} - (C^{n+1} p^{n+1})|_{ij}}{\Delta y} + \frac{D_{ij+1}^{n+1} - D_{ij}^{n+1}}{\Delta y} \right) \\ & - L_{ij-1/2} k_{ij-1/2} \left( \frac{(C^{n+1} p^{n+1})|_{ij} - (C^{n+1} p^{n+1})|_{ij-1}}{\Delta y} + \frac{D_{ij}^{n+1} - D_{ij-1}^{n+1}}{\Delta y} \right) \\ & + L_{\Gamma} \left( k_{\Gamma} \frac{\partial T}{\partial n} \right) \Bigg|_{\Gamma}^{n+1} \end{aligned} \quad (65) $$

We provide now some general guidelines to consider an IDBC, as Eq. (27). If the solid interface crosses the segment $[i, j : i+1, j]$, the following scheme is obtained,

$$ \begin{aligned} \iint_{\partial \Omega_{ij}} k \nabla T^{n+1} \cdot \vec{n}_{\Sigma} \, d\Sigma = & k_{i+1/2j} \left( \frac{T_{wall} - T_{ij}^{n+1}}{\theta \Delta x} \right) \Delta y - k_{i-1/2j} \left( \frac{T_{ij}^{n+1} - T_{i-1j}^{n+1}}{\Delta x} \right) \Delta y + k_{ij+1/2} \left( \frac{T_{ij+1}^{n+1} - T_{ij}^{n+1}}{\Delta y} \right) \Delta x \\ & - k_{ij-1/2} \left( \frac{T_{ij}^{n+1} - T_{ij-1}^{n+1}}{\Delta y} \right) \Delta x, \end{aligned} \quad (66) $$

which can be solved with a similar approach as INBC by injecting the relation $T^{n+1} = C^{n+1}P^{n+1} + D^{n+1}$. More details can be found on this specific point in [45] in which such an implementation has been carried out to impose the saturation temperature at a liquid/vapor interface.

### 3.4. Laplace equation for the electrical potential

As previously stated, one of the objectives of this paper is to couple our fluid solver with some thermal plasma modules, as a Laplace equation to compute for the electric potential, Eq. (14). From the developments previously presented, applying IDBC or INBC to this Laplace equation is straightforward. Such conditions will be imposed on the electrodes during the discharge of an electric arc to induce an electric current.

### 3.5. Numerical methods and time step

All convective derivatives are computed with WENO-Z (Weighted Essentially Non-Oscillatory-Z) proposed in [66], which combines good stability properties and high accuracy. Other terms are discretized with classical second order finite volume schemes, for which adaptations have been presented previously for cells crossed by the solid frontier. Temporal integration is performed with a second order Runge-Kutta scheme. The pressure linear system resulting from the spatial discretization of Eq. (37) or Eq. (41) is solved with the Black-Box Multi-Grid solver (BBMG) proposed in [67] and benchmarked in [68] against several other solvers for linear systems. Finally, constraints have to be imposed on the time step to ensure numerical stability of the computations. Time step constraints, related to convection terms, viscous terms and heat conduction terms must be defined,
$$ \Delta t_{conv} = \frac{h}{\max||\vec{V}||} \qquad \Delta t_{visc} = \frac{1}{2} \frac{h^2}{\max(v)} \qquad \Delta t_{heat} = \frac{1}{2} \frac{h^2}{\max(D_{th})} \qquad (67) $$

with h the length of a computational cell, v = μ/ρ the kinematic viscosity and Dth = k/ρCρ the thermal diffusivity. Finally, the overall time step can be deduced from the following expression,

$$ \Delta t = \left( \frac{1}{\Delta t_{conv}} + \frac{1}{\Delta t_{visc}} + \frac{1}{\Delta t_{heat}} \right)^{-1} \qquad (68) $$

It is worth mentioning that Δtvisc can be removed with an implicit temporal discretization, as in Eq. (40). On the other hand, regarding the time step constraint related to heat conduction Δtheat, the temporal discretization proposed in Eq. (41), is not fully implicit due to the pre-factor multiplying the temperature Laplacian which makes non-linear this term. Therefore, it does not fully remove this time step constraint, but enables to alleviate it. Regarding the propagation of acoustic waves, no additional time step constraints have to be imposed to ensure the stability of the computations, since the temporal integration of acoustic terms is implicit. However, describing accurately acoustic phenomena requires solving temporal characteristic scales, and therefore a suitable time step.

### 3.5.1. Field extensions

As several variables, as the pressure p or the density ρ are not defined in the solid domain, computations of their convective derivatives $\vec{V} \cdot \nabla p$ or $\vec{V} \cdot \nabla \rho$ in fluid cells close to the solid frontier require a specific care. For this purpose, extrapolation techniques from [69] are commonly used to build continuous ghost fields of the fluid variables on the other side of the solid interface. This can be carried out by solving the following partial differential equation,

$$ \frac{\partial p}{\partial \tau} + \vec{n} \cdot \nabla p = 0 \qquad (69) $$

for the pressure field (or a similar equation for the density) over a fictitious time τ, with $\vec{n}$ the normal vector at the solid surface. An iterative resolution of these equations is carried out until a steady state is reached in a sufficiently broad region around the interface. This allows to maintain the required property of the extrapolated field across the interface, that is, $\vec{n} \cdot \nabla \rho = 0$ and $\vec{n} \cdot \nabla p = 0$. In practice, this is achieved with approximately 10–20 temporal iterations to obtain a relevant extended field inside the solid domain. More details on extrapolation techniques can be found in [70] in the framework of Stefan problems and in [71] for boiling flows. Extensions with an improved accuracy based on the subcell resolution at the interface have also been described in some more recent works [72,73].

## 4. Numerical results

Numerical results and validations of the overall solver are presented in this section. The interest of the coupling between the compressible solver and an IBM is highlighted on different kind of flows for which density variations and compressibility effects are significant. In a first time, external and internal free-convection flows are considered. In order to demonstrate the suitable behavior of the solver to perform numerical simulations of a liquid in the low Mach number regime, a benchmark involving a liquid flow in a channel with a cross-section reduction has been carried out. Numerical simulations of the flow around a cylinder in a weakly compressible regime are also presented and compared to anterior works. Next, computations involving acoustic waves propagation are presented in both configurations, that is, the propagation of a plane wave in a tube with a section reduction and the diffusion of a plane wave by a sphere. Finally, in order to demonstrate the interest of the proposed solver for thermal plasma simulations, preliminary results of an electric arc discharge are also shown. Simulations presented in this paper are carried out in 2D with variables (x, y) or in axisymmetric coordinates with variables (r, z). The extension of the overall formalism to 3D configurations is straightforward, provided a 3D interface reconstruction is available to compute the term $L_Γ$, as detailed in [62,63]. Many 3D examples are presented in the following references [1,58,64,65] for incompressible flows. In configurations involving propagation of acoustic waves, the time step accounts for the following constraint $\Delta t_{acoustic} = \frac{1}{4} \frac{h}{\max(c)}$ which is not required for stability reasons, but enables an accurate description of acoustic waves.

The thermophysical properties of air, at 300 K and ambient pressure, are given here, since it is used in benchmarks described in subsections 4.1, 4.2 and 4.6. The dynamic viscosity is $\mu = 18.1 \times 10^{-6}$ Pa s, the thermal conductivity is $k = 0.026$ W m$^{-1}$ K$^{-1}$, the specific heat at constant pressure is $c_p = 1000$ J kg$^{-1}$ K$^{-1}$, the gravitational acceleration g is 9.81 m s$^{-2}$. From ideal gas EOS, the density of air ($M_{air} \sim 0.029$ kg mol$^{-1}$) can be computed as 1.177 kg m$^{-3}$ at atmospheric pressure and 300 K.

### 4.1. Free convection around a horizontal cylinder in a rectangular closed cavity

In this first benchmark, the free convection flow around a heated horizontal cylinder is considered. This will enable to validate the proposed approach to carry out a correct coupling between heat conduction, buoyancy effects induced by density variations and dynamical effects. The computational configuration is based on the experimental work presented in [74], which consists in a heated
<figure><img src="image_7.png" /><figcaption>Fig. 5. Configuration (a) and temperature field snapshot (b).</figcaption></figure>

long cylinder of diameter D placed inside a rectangular cavity of height H and width W. The temperature of the heated cylinder is $T_h$, and the Rayleigh number can be defined as,

$$Ra = \frac{\rho g (T_h - T_c) D^3}{\nu D_{th}}, \qquad (70)$$

where $T_h$ is the cylinder temperature which is considered uniform, $T_c$ is the external temperature imposed on the lateral sides. The height and the width of the computational domain are taken respectively at $H = 57$ mm and $W = 30$ mm. The fluid in the computational domain is air and is considered as an ideal gas. Boundary conditions on the lateral sides maintain a constant fluid temperature $T_c$. Wall boundary conditions are considered on all sides of the computational domain and on the Immersed Boundary. From the point of view of the thermal solver, the temperature $T_h$ on the Immersed Boundary is imposed with the IDBC scheme as defined previously for the temperature Laplacian, as in Eq. (66). An adiabatic boundary condition, $k \frac{\partial T}{\partial n}|_w = 0$, is imposed at the bottom of the computational domain. On the upper side a flux continuity is imposed $k \frac{\partial T}{\partial n}|_w = h(T_w - T_c)$ where $T_w$ is the local wall temperature and $h = 10 \text{ W m}^{-2} \text{ K}^{-1}$ a heat transfer coefficient related to external convection. As the wall temperature $T_w$ is not known a priori, the previous condition is a Robin boundary condition. This set of boundary conditions is similar to those proposed in the numerical part of the study presented in [74]. A schematic of the computational domain is plotted in Fig. 5a. To investigate on grid sensitivity, simulations have been performed on four successive computational grids, 64 × 128, 128 × 256, 256 × 512, 512 × 1024. Simulations are integrated until reaching a steady state at $t \sim 5s$. The presented configuration corresponds to a Rayleigh number $Ra = 1300$.

The Nusselt number which characterizes the interaction between heat conduction and free convection is extracted from the simulations and compared to experimental data presented in [74]. It is defined as the following integral along the surface $S_Γ$ of the whole cylinder,

$$Nu = \frac{1}{\pi} \iint_{S_Γ} \frac{\nabla T}{T_h - T_c} dS \qquad (71)$$

Snapshots of the temperature field and streamlines around the heated cylinder are plotted in Fig. 5b at steady state. This highlights the good qualitative behavior of the flow around the cylinder with tangent streamlines to the solid frontier, and thermal boundary layer around the solid obstacle. In Table 1, the values of the Nusselt number extracted from the simulations are reported for the four different computational grids. Space convergence is clearly observed with a numerical result approaching around ~5% to the expected experimental value, which is a reasonable agreement between a numerical simulation and experimental results in the framework of free-convection flow. More detailed investigations on the velocity components and temperature profiles along the x-axis and the y-axis (respectively $y = 0$ and $x = 0$) are plotted in Fig. 6, 7 and 8 for the four different grids, except for the temperature profile for which only the results on the coarsest and thinnest grids are plotted, since all curves are superimposed. This clearly shows the space convergence of the numerical solution. Although no analytical solutions exist on this configuration, the order of convergence of the solver has been approximately estimated by considering the solution on the more refined grid 512 × 1024 as a reference solution. For this purpose, the error has been analyzed at a point of coordinates ($x = 12.4$ mm, $y = 0$) where the x-component velocity $u$ reaches a peak, as observed in Fig. 6b. These results are presented in Table 2, and show an approximately second order accuracy, which is relevant with the IDBC and INBC schemes which are both second order.
<table><caption>**Table 1**

Comparisons between the computed Nusselt number and the experimental Nusselt number for two Rayleigh numbers and different meshes.</caption><thead><tr><th>Mesh</th><th>Nu for Ra = 1300</th></tr></thead><tbody><tr><td>64 × 128</td><td>2.102</td></tr><tr><td>128 × 256</td><td>2.234</td></tr><tr><td>256 × 512</td><td>2.298</td></tr><tr><td>512 × 1024</td><td>2.324</td></tr><tr><td>experimental data [74]</td><td>2.46</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>(a) x-velocity component along x = 0</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) x-velocity component along y = 0</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(a) y-velocity component along x = 0</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) y-velocity component along y = 0</figcaption></figure>

<caption>**Fig. 7.** y-velocity component profiles along x and y axis for four different grids.</caption>

## 4.2. Free-convection inside a horizontal hollow cylinder

The previous benchmark was focused on a free-convection flow around a heated solid. In this test-case, an internal free-convection flow inside an infinite horizontal cylinder is considered. A schematic of the computational domain is plotted in Fig. 9 which enables to visualize orientation of the gravity vector with regards to the cylinder orientation. Free-convection inside the horizontal cylinder
<figure><img src="image_12.png" /><figcaption>Fig. 8. Temperature profiles along x and y axis for two different grids (the coarsest and the thinnest).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 9. Schematic of the computational domain (in blue) with angle $\theta$ and the gravity vector $\vec{g}$.</figcaption></figure>

<caption>**Table 2**</caption>

Error analysis on the x-component velocity $u$ at point $x = 12.4$ mm and $y = 0$ from numerical results on four different grids.

<table><thead><tr><th>Mesh</th><th>u (mm/s)</th><th>% relative error</th><th>approximate order</th></tr></thead><tbody><tr><td>64 × 128</td><td>0.2087</td><td>17.18%</td><td>-</td></tr><tr><td>128 × 256</td><td>0.1848</td><td>3.76%</td><td>2.2</td></tr><tr><td>256 × 512</td><td>0.1801</td><td>1.12%</td><td>1.75</td></tr><tr><td>512 × 1024</td><td>0.1781</td><td>-</td><td>reference</td></tr></tbody></table>

can be induced by imposing higher temperatures at the bottom of the cylinder than at the top of the cylinder. For this purpose, the following angular distribution of the wall temperature can be imposed,

$$T_w(\theta) = T_0 + A_0 \cos\theta, \qquad (72)$$

with $T_0 = 300$ K, $A_0 = 30$ K and $\theta$ is defined in Fig. 9. The internal free-convection flow can also be induced with an angular distribution of the local heat flux at the solid frontier,

$$\varphi_w(\theta) = B_0 \cos\theta, \qquad (73)$$

with $B_0$ the amplitude of the imposed heat flux. $D_c$ is the diameter of the cylinder for which three different values have been considered, $D_c = 0.01$ m, 0.015 m and 0.02 m, corresponding respectively to three Rayleigh numbers, $Ra = 6\,800, 23\,000$ and 54 000, based on a temperature variation $\Delta T = 2A_0$ in configurations involving a prescribed temperature boundary condition. If a local heat flux is imposed on the solid frontier, Eq. (73), identical Rayleigh numbers are obtained for the same values of $D_c$ by taking the following expression on the amplitude of the local heat flux,

$$B_0 = \frac{\lambda A_0}{D_c}. \qquad (74)$$
<figure><img src="image_14.png" /><figcaption>Fig. 10. Temperature field and streamlines at steady state for the three different Rayleigh numbers with an imposed temperature (IDBC) for the grid 256 × 256.</figcaption></figure>

Both kind of boundary conditions are respectively imposed with the IDBC or INBC schemes. As in the previous section, the fluid inside the cylinder is air, and it is considered as an ideal gas. The initial temperature is taken at $T_0 = 300$ K and the initial pressure is atmospheric pressure $P_0 = 101325$ Pa.

Numerical simulations have been performed with three different meshes, 64 × 64, 128 × 128 and 256 × 256 in each configuration, as well as a computation with a much more refined grid 1024 × 1024 to serve as a reference solution for the configuration with $Ra = 23 000$ and a Dirichlet boundary condition on the temperature at the solid frontier. Snapshots of the temperature field and streamlines of the free-convection flow are represented in Fig. 10 at steady state for the four Rayleigh numbers. For the lower Rayleigh number, $Ra = 6 800$, we observe a thermal stratification, with the heated gas which remains at the bottom of the computational domain, although streamlines highlight a symmetric recirculating free-convection motion in the cylinder. This illustrates that conduction heat transfer has a prominent effect on the thermal field in this first configuration. By increasing the Rayleigh number to $Ra = 23 000$, the thermal stratification is not maintained, and the impact of the free-convection symmetric recirculating flow on the thermal field is clearly visualized. Finally, considering the higher Rayleigh number, $Ra = 54 000$, one can observe the symmetry breaking of the free-convection flow and the transition to a flow rotating all along the cylinder frontier. In this last case, the strong coupling between the thermal field and the free-convection motion is also clearly observed. In Fig. 11, snapshots of the temperature field and streamlines are shown in the second configuration in which a local heat flux boundary condition is imposed on the solid frontier for $Ra = 23 000$. We take interest now to the conservation of mass and energy in the whole computational domain to assess the correct formulation of the boundary conditions of our set of equations on the solid frontier. The total mass in the whole domain is computed from the following volume integral,

$$m = \int_{\Omega_f} \rho d\Omega_f, \qquad (75)$$

where $\Omega_f$ is the total fluid volume. The total energy of the gas is calculated as the sum of the volume integral of the kinetic, potential and internal energies and the surface and temporal integral of the local heat flux, such as,

$$E = \iiint_{\Omega_f} \frac{\rho \vec{u}^2}{2} dV + \iiint_{\Omega_f} \rho e dV + \iiint_{\Omega_f} \rho g y dV + \int_{t=0}^{t_f} \left( \iiint_{\Sigma_f} k \nabla T \cdot \vec{n}_{\Sigma_f} d\Sigma \right) dt \qquad (76)$$

where $\Sigma_f$ is the surface delimiting the volume $\Omega_f$, $\vec{n}_{\Sigma_f}$ the normal vector pointing outward and $y$ is the vertical coordinate and $e$ the internal energy defined by $de = C_v dT$ for an ideal gas. It is noteworthy that the last term which expresses the heat flux transferred outward of the computational domain should be zero if Eq. (73) is imposed on the boundary condition, whereas it may not be the case if a Dirichlet boundary condition, Eq. (72), is imposed on the temperature.
<figure><img src="image_15.png" /><figcaption>Fig. 11. Temperature field and streamlines at steady state for the Rayleigh numbers Ra = 23000 with an imposed local heat flux (INBC) for the grid 256 × 256.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Temporal evolution of the total energy and mass for different grids 64x64, 128x128, 256x256 with Ra = 23000 for INBC.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Temperature, x-component and y-component velocity profiles along the y = 0 axis for two different grids (the coarsest and the thinnest) in the case Ra = 23000.</figcaption></figure>

Dimensionless mass and energy are defined respectively as $\bar{m}(t) = \frac{m(t)}{m_0}$ and $\bar{E}(t) = \frac{E(t)}{E_0}$, where $m_0$ and $E_0$ are the total mass and energy at the initial time in the computational domain. In Fig. 12 the temporal evolutions of $\bar{m}(t)$ and $\bar{E}(t)$ are plotted for the three different meshes and with INBC, that is, a heat flux boundary condition. This highlights that conservation errors are less than 0.02% for the coarsest grid both on mass and energy and decrease adequately by refining the grid. This lack of conservation can be related to the Finite Difference schemes used to discretize the convection terms in a primitive formulation. Indeed, it is well known that such formulation is not strictly conservative. On the other hand, the discretizations of diffusion terms can be expressed as Finite Volume schemes and do not alter conservation properties. To investigate further on the overall accuracy of the proposed solver, temperature, x-component and y-component velocity profiles have been plotted along the x-axis ($y = 0$) in Fig. 13. As all curves are almost superimposed, only the thinnest and coarsest grids have been plotted. An error analysis on the temperature and y-component velocity at ($x = 0, y = 0$) is presented in Table 3 and show an accuracy order around 1.4 – 1.5 on this test-case.
<table><caption>**Table 3**
Error analysis on the y-component velocity $v$ and temperature at point of coordinates $(x = 0, y = 0)$ from numerical results on four different grids.</caption><thead><tr><th>Mesh</th><th>v (mm/s)</th><th>% relative error v</th><th>order</th><th>T (K)</th><th>% relative error T</th><th>order</th></tr></thead><tbody><tr><td>64 × 64</td><td>29.735</td><td>4.33%</td><td>-</td><td>309.293</td><td>0.083%</td><td>-</td></tr><tr><td>128 × 128</td><td>30.623</td><td>1.48%</td><td>1.55</td><td>309.452</td><td>0.032%</td><td>1.37</td></tr><tr><td>256 × 256</td><td>30.909</td><td>0.557%</td><td>1.41</td><td>309.512</td><td>0.012%</td><td>1.41</td></tr><tr><td>1024 × 1024</td><td>31.082</td><td>-</td><td>reference</td><td>309.55</td><td>-</td><td>reference</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 14.** Pressure field and streamlines in the convergent pipe.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Density field in the convergent pipe.</figcaption></figure>

### 4.3. Liquid flow in a convergent pipe in the incompressible regime

We propose in this subsection to carry out numerical simulations of a liquid flow in a convergent pipe in the incompressible regime, to demonstrate the relevance of the solver to configurations with a more complex EoS than the simple ideal gas one. Indeed, the consistency of the proposed IBM will be assessed for a cubic EoS which is well suited to describe liquid thermodynamic states. It is worth mentioning, that this configuration is similar to the one presented in [45] in which bubble cavitation has been computed, although we just focus here on a single-phase flow. The solid frontier defined by the Immersed Boundary is included in a cylindrical computational domain of height $H = 0.8$ mm and radius $R = 0.1$ mm. The inner space inside the solid frontier is filled with liquid $CH_4$ at $T = 120$ K and $p_0 = 1.87372$ bar. In addition, the following physical properties are considered: $C_v = 1643$ J kg$^{-1}$ K$^{-1}$, $k = 0.172$ W m$^{-1}$ K$^{-1}$ and $\mu = 98.68$ $\mu$ Pa s. A pressure drop of 0.4 bar is imposed between the inlet and outlet of the nozzle, with $p_{inlet} = 1.87372$ bar at the inlet and $p_{outlet} = 1.47372$ bar at the outlet. The radius of the entry section is $R_{inlet} = 0.112$ mm, and at the throat, the section is $R_{outlet} = 0.310^{-4}$ mm, which corresponds to a velocity ratio $\sim 14$ between the outlet velocity and the inlet velocity. At steady state, the pressure and velocity fields are established, the outlet velocity is about $V_{outlet} \sim 14$ m s$^{-1}$, and the inlet velocity is about $V_{inlet} \sim 1$ m s$^{-1}$. In the present range of temperature and pressure, the density of the liquid methane according to the SRK EoS remains nearly constant at approximately 406 kg m$^{-3}$. Considering the liquid viscosity and channel dimensions, the corresponding Reynolds number is low enough in order to ensure that the flow will remain in the laminar regime. The simulations have been performed with three successive grids, containing respectively 256 × 2048, 512 × 4096 and 1024 × 8192 elements in the $r-z$ directions. In the reduced section, the pressure drop is observed in Fig. 14 from $p_{inlet} \sim 1.8 \times 10^5$ Pa to $p_{outlet} \sim 1.4 \times 10^5$ Pa. It is induced by both the velocity increase in the reduced section and the viscous dissipation. In Fig. 15, the density field is plotted and shows density variations related to the pressure field. As expected for a liquid flow in the incompressible regime, these variations are very weak around $\sim 0.1\%$. The temporal variation of the inlet mass flow rate $q_{in}$, which results from this pressure drop, is plotted in Fig. 16a for different grids 256 × 2048, 512 × 4096 and 1024 × 8192. In this figure, one can observe that the inlet mass flow rate reaches a steady value after a preliminary increase in the initial transient state. This figure also clearly demonstrates the convergence of the computed solutions, with very close results for the three grids. Finally, in Fig. 16b, the inlet and outlet mass flow rate $q_{out}$ are both plotted in order to check the correct mass conservation of the liquid between the inlet section and the outlet section. The figure clearly shows the good agreement between $q_{in}$ and $q_{out}$ at each time, which validates the correct behavior of the proposed compressible IBM for an internal liquid flow in the incompressible regime.

### 4.4. Flow around a cylinder at $Re = 100$ and $Ma = 0.25$

For the sake of comparisons of the proposed solver with anterior studies, we present in this subsection some results of a flow around a cylinder in a weakly compressible regime, that is for a Mach number $Ma = 0.25$, and for a moderate Reynolds number $Re = 100$. In
<figure><img src="image_20.png" /><figcaption>Fig. 16. Inlet and outlet mass flow rate, grid sensitivity and conservation.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 17. Velocity magnitude and pressure fields for the flow around a cylinder at Ma = 0.25 and Re = 100.</figcaption></figure>

particular, this benchmark has already been presented in anterior studies, as for example in [75,76]. The sizes of the computational
domain are $l_x = l_y = 40$ mm. The radius of the cylinder is $D_{cyl} = l_x/20$. Numerical simulations have been performed with three grids
of increasing resolution, $512 \times 512$, $1024 \times 1024$ and $2048 \times 2048$. The physical properties are properties of air at $T_\infty = 300$ K, except
for the viscosity $\mu = 2.10^{-3}$ Pa.s which has been increased in order to impose both a rather high velocity $V_\infty = 85$ m.s$^{-1}$ ($Ma = 0.25$)
and a low Reynolds number $Re = 100$. Snapshots of the velocity magnitude and pressure fields are plotted in Fig. 17a and 17b
respectively. This allows to visualize the complex flow structures that appear in the wake of the cylinder. Appearance of so-called
von-Karman streets is relevant with the Reynolds number $Re = 100$, since such flow instabilities develop in the wake of a circular
cylinder for $Re > 50$. Comparisons of our numerical results with anterior works on the average drag coefficient $C_{D,mean}$, the maximum
lift coefficient $C_{L,max}$ and Strouhal number, defined as $St = \frac{fD}{V_\infty}$ with f the oscillation frequency, are presented in Table 4. These
results show a reasonable agreement with anterior works with deviations between 5% to 10% for $C_{D,mean}$ and $C_{L,max}$. It is worth
mentioning that for a Reynolds number $Re = 100$, the expected value of the drag coefficient of a circular cylinder is $C_{D,mean} = 1.35$
in the incompressible regime. However, it is well known that the drag coefficient increases with the Mach number even for subsonic
flows ($0.1 < Ma < 1$) for which compressibility effects are weak but not negligible. With a $C_{D,mean} = 1.48$, our numerical results are
relevant with this trend. Moreover, the Strouhal number, which is expected at $St = 0.17$ for $Re = 100$, is accurately predicted in our
numerical results with a deviation to the expected value around 0.6%.
<caption>**Table 4**</caption>

Average drag coefficient, maximum lift coefficient and Strouhal number for different grids and comparisons with anterior studies.

<table><thead><tr><th>Mesh</th><th>C<sub>D,mean</sub></th><th>C<sub>L,max</sub></th><th>S<sub>t</sub></th></tr></thead><tbody><tr><td>512 × 512</td><td>1.52</td><td>0.332</td><td>0.173</td></tr><tr><td>1024 × 1024</td><td>1.49</td><td>0.315</td><td>0.171</td></tr><tr><td>2048 × 2048</td><td>1.48</td><td>0.343</td><td>0.172</td></tr><tr><td>Canuto & al [76]</td><td>1.38</td><td>0.325</td><td>0.163</td></tr><tr><td>Karagözis & al [75]</td><td>1.34</td><td>0.319</td><td>0.168</td></tr></tbody></table>

<figure><img src="image_22.png" /><figcaption>**Fig. 18.** Schematic representation of the test case with 2 surfaces $S_1$ and $S_2$ and the pressure wave on the right with a pressure drop $p_0 = p_{atm} - p$.</figcaption></figure>

### 4.5. Propagation of an acoustic wave in a rigid tube with a sudden reduction of the cross section

The objective of this benchmark is now to investigate on the ability of the present solver to perform numerical simulations of the interaction between acoustic waves and an immersed solid frontier. The computational domain consists in a 2D rigid tube with a sudden reduction in the cross-section, as illustrated in Fig. 18.

Its total dimensions are $l_x = 3$ m and $l_y = 0.3$ m, the cross-section in the narrower part of the tube is $S_1 = 0.1$ m and $S_2 = 0.2$ m in the larger part of the tube. An acoustic wave is initialized on the right side of the domain, where the cross-section is larger and propagates towards the left direction. The initial profiles of the pressure, density and x-component velocity are given by the following relations,

$$p'(x, y, 0) = p_0 \exp\left(-\left(\frac{x - x_0}{\sigma_a}\right)^2\right) \qquad \rho(x, y, 0) = \frac{p'(x, y, 0)}{c^2} \qquad u(x, y, 0) = \frac{p'(x, y, 0)}{\rho c} \qquad (77)$$

which describe a Gaussian acoustic pulse, centered at $x = x_0$ with an amplitude $p_0 = -100$ Pa, and which propagates towards the left direction. $\sigma_a$ is taken at 0.1 m. The solid walls are adiabatic, and the inner gas is steam water ($M_{H_2O} = 0.018 \text{ kg mol}^{-1}$), considered as an ideal gas, with an initial temperature equal to 375.26 K. The steam thermophysical properties are: $\mu = 12 \times 10^{-6}$ Pas, $k = 0.0246 \text{ W m}^{-1} \text{ K}^{-1}$, $c_p = 2000 \text{ J kg}^{-1} \text{ K}^{-1}$ and the density at 375 K is $0.59 \text{ kg m}^{-3}$ which is obtained from the ideal gas EOS. The numerical simulations have been carried out with four different meshes, $64 \times 384$, $128 \times 768$, $256 \times 1536$ and $512 \times 3072$. In such a configuration, it is expected that when the acoustical pulse interacts with the sudden reduction of the cross-section, a part of the acoustic wave is transmitted and the other part is reflected. In Fig. 19, the pressure profile is plotted at three different times for the most refined grid. At $t = 1$ ms, the pressure profile can be visualized before the impact of the acoustic wave with the cross-section reduction. At $t = 2$ ms the acoustical wave is impacting the cross-section reduction, and $t = 3$ ms shows the pressure profile after the impact. These snapshots demonstrate that the numerical simulation reproduces the expected behavior, namely the splitting of the initial wave in a transmitted wave and a reflected wave after the impact with the cross-section reduction, as observed in Fig. 19c. To assess the conservation properties of the solver in a case involving acoustic wave propagation, the temporal evolution of the acoustic energy has been plotted in Fig. 20a for the four different meshes. For this purpose, the total acoustic energy $E_a$ has been computed as the sum of $E_p$, the energy related to pressure effect, and of the kinetic energy $E_k$, with the following expression,

$$E_a = E_p + E_k = \int_{\Omega_f} \left( \frac{p'^2}{\rho_0 c_0^2} + \frac{\rho_0 u^2}{2} \right) d\Omega_f \qquad (78)$$

where $p' = p - p_0$. We define also $\tilde{E}_a = \frac{E_a}{E_{a,0}}$, as the acoustic energy normalized by the initial acoustic energy $E_{a,0}$. As the flux of acoustic waves across solid boundaries is zero and viscous dissipation and heat conduction effects are negligible in this case, the total acoustic energy should be constant in time. Fig. 20a demonstrates clearly that conservation errors on the acoustic energy along the computation decrease when refining the computational grid. Error analysis on $\tilde{E}_a$ at final time $t = 4$ s is presented in Table 5. Accuracy order around ~1 is reported for this benchmark. This low order accuracy can be explained considering that the viscous effects are not fully resolved in this configuration. Indeed, even if refined grids have been used, it is worth mentioning that the grid cells remain too coarse to fully resolve the very thin acoustic boundary layers related to viscous effects at very low scale. However, this is relevant with the inviscid theory to which are compared our simulations.
<figure><img src="image_23.png" /><figcaption>Fig. 19. Pressure profile in the tube at different times for the grid 512 × 3072.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20. Temporal evolution of dimensionless energies $\tilde{E}_a$.</figcaption></figure>

**Table 5**

<table><caption>Error analysis on $\tilde{E}_a$ at final time $t = 4$ s for four different grids.</caption><thead><tr><td>Mesh</td><td>|1 - E&#773;<sub>a</sub>|</td><td>% relative error</td><td>order</td></tr></thead><tbody><tr><td>64 × 384</td><td>0.84</td><td>16%</td><td>-</td></tr><tr><td>128 × 768</td><td>0.918</td><td>8.2%</td><td>0.96</td></tr><tr><td>256 × 1536</td><td>0.978</td><td>4.2%</td><td>0.96</td></tr><tr><td>512 × 3072</td><td>0.982</td><td>1.8%</td><td>1.2</td></tr></tbody></table>

In Fig. 20b, the two components of the total energy have also been plotted. In this figure, one can remark that both contributions to the total energy are almost equal during the propagation of the acoustic wave, except during the impact with the cross-section reduction, where $E_p$ increases and $E_k$ decreases, while $E_a$ remains constant. After the impact, both contributions recover their initial value.
<figure><img src="image_25.png" /><figcaption>Fig. 21. Comparison between numerical simulations and theory on the variations of $\tilde{p}'_R$ and $\tilde{p}'_T$ with $\tilde{S}$ for the grid 512 × 3072. The blue curve is the transmitted pressure and the orange curve is the reflected pressure ratio. Squares and triangles are for horizontal and inclined tubes, respectively.</figcaption></figure>

As an existing theory enables to determine the respective intensities of the transmitted and reflected waves, we have performed a parametric study with other cross-sections $S_1$, such as the ratio $\frac{S_1}{S_2}$ ranges between 0.1 to 0.9, in order to compare the numerical simulations to the theoretically predicted values.

This theory based on linear non-viscous acoustics specifies that the pressure varies as follows in a tube of changing cross-sections,

$$p'_{R} = p'_{I} \frac{Y_{1} - Y_{2}}{Y_{1} + Y_{2}} \qquad (79)$$

$$p'_{T} = p'_{I} \frac{2Y_{1}}{Y_{1} + Y_{2}} \qquad (80)$$

where $Y_1$ and $Y_2$ are the admittance of the tube in the two different regions, $P_I$, $P_R$ and $P_T$ are, respectively, the incident, reflected and transmitted wave pressure. The characteristic admittance of a tube is defined as,

$$Y = \frac{S}{\rho c} \qquad (81)$$

The acoustic admittance ratio is given by the following relation,

$$\frac{Y_1}{Y_2} = \frac{S_1}{S_2}, \qquad (82)$$

if the fluid density and sound speed are uniform in the whole domain. From these relations, the variations of the reflected and the transmitted pressure $\tilde{p}'_R$ and $\tilde{p}'_T$ with the cross-section ratio $\tilde{S} = \frac{S_1}{S_2}$ can be determined.

In Fig. 21, a comparison between the theoretical predictions and the computed values of $\tilde{p}'_R$ and $\tilde{p}'_T$ is plotted for several cross-section ratios. This shows a satisfactory agreement between the theory and the numerical simulations and validates the interaction of an acoustical wave with a cross-section reduction. Finally, to assess the proposed solver in a configuration for which the solid frontiers are not aligned with the mesh axis, additional simulations have been performed with a tube axis inclined of an angle 30° with the mesh axis. These simulations have been carried out for two cross-section ratios, that is $\tilde{S} = 0.5$ and $\tilde{S} = 0.7$. Values of $\tilde{p}'_R$ and $\tilde{p}'_T$ are also reported in Fig. 21 with triangle. For these configurations, the overall domain has been modified to fit with the tube geometry. The lengths of the computational domain are $l_x = 3$ m and $l_y = 2$ m and the number of grid points is 3072 × 2048. This corresponds to a size of the computational cells $\Delta x \approx 1$ mm. Snapshots of the pressure field are plotted in Fig. 22 at four different times to illustrate the suitable behavior of the acoustic wave when impacting the cross-section reduction, with the splitting of the incident wave in a reflected wave and a transmitted wave.

## 4.6. Diffusion of an acoustic wave by a sphere

After considering the propagation of an internal acoustic wave in a closed tube with a cross-section reduction, we take interest in this subsection in the interaction of an acoustic wave with a solid sphere. This will enable to assess the accuracy of the presented solver
<figure><img src="image_26.png" /><figcaption>Fig. 22. Pressure profile in the tube with an inclination of 30° at different times for a grid 3072 × 2048.</figcaption></figure>

for an acoustical wave interacting with a curved solid surface. The numerical simulations are performed in axisymmetric coordinates. At the initial time, a plane wave, defined with a Gaussian pressure fluctuation,

$$p'(r, z, 0) = p'_{0} \exp \left( -\left( \frac{z - z_{0}}{\sigma_{a}} \right)^{2} \right), \qquad (83)$$

with $p'_{0} = 100$ Pa, $\sigma = 0.0001$ m and $z_{0} = -l_{z}/5$, is propagating in the direction of a solid sphere of radius $R_{S} = 0.42$ mm. $\rho(r, z, 0)$ and $v(r, z, 0)$ are deduced from the initial pressure according to Eq. (77). The lengths of the computational field are $l_{r} = 0.0025$ m and $l_{z} = 0.005$ m. The gas is air, still considered as an ideal gas. The initial temperature of the gas is 300 K and the initial pressure is $P_{i} = 101325$ Pa. To demonstrate the space convergence of the simulations, all spatial and temporal scales of the flow have to be accurately resolved. In particular, viscous effects will develop a boundary layer around the solid obstacle, which thickness $\eta$ can be approximated by the following relation [77,78],

$$\eta \sim \sqrt{\frac{\nu}{f}}, \qquad (84)$$

where $\nu$ is the kinematic viscosity, $f$ is the frequency of the acoustic wave. As acoustic boundary layers are very thin, kinematic viscosity has been increased to $\nu = 2.85 \cdot 10^{-4}$ m$^2$s$^{-1}$ in order to be able to capture the acoustic boundary layer. For a Gaussian
<figure><img src="image_27.png" /><figcaption>Fig. 23. Pressure field at different times during the diffusion of an acoustical wave by a sphere for the grid 2048 × 4096.</figcaption></figure>

acoustic pulse, wavelength $\lambda_w$ is approximately taken equal to $\sigma_a$, which gives the following frequency, $f = 3.47$ MHz (corresponding to an ultrasonic acoustic wave), for a sound speed $c = 347$ m s$^{-1}$. From this, we obtain an acoustic boundary layer thickness $\eta \sim 10$ µm. Numerical simulations have been performed with the five following meshes of increasing resolution, 128 × 256, 256 × 512, 512 × 1024, 1024 × 2048 and 2048 × 4096 which correspond respectively to the following sizes of the computational cells, 19.6 µm, 9.8 µm, 4.9 µm, 2.4 µm, 1.2 µm. As at least 4 computational cells are required inside a boundary layer to achieve a satisfactory convergence, this first estimation suggests that only the two more refined grids will provide sufficiently accurate results. In Fig. 23, the pressure field has been plotted at different times, $t_0 = 0$ s, $t_1 = 3.10^{-6}$ s, $t_2 = 5.10^{-6}$ s and $t_3 = 7.10^{-6}$ s from the computation with the more refined grid 2048 × 4096. At the initial time $t_0$, the initial wave is propagating toward the solid sphere. At time $t_1$, the impact of the acoustical wave with the sphere can be visualized, with, in particular, the formation of a spherical wave due to the diffusion of the incident acoustic wave. At the time $t_2$, the incident acoustic wave has almost crossed the sphere, and a larger pressure is observed at the top of the sphere where both parts of the incident wave are impacting. At the last time, $t_3$, the incident wave moves away from the sphere and a secondary circular acoustic wave is formed after the impact of both parts of the incident wave at the top of the sphere, as previously observed in Fig. 23c. To illustrate the grid independence of the numerical solution, the pressure and velocity profiles have been plotted at different times along r-axis for two possible values of z, that is $z_1 = 0$ which corresponds to the center of the cylinder and $z_2 = R_c$ which corresponds to the top of the cylinder. In Fig. 24, the profiles are plotted at $z = z_1$ and $t = t_1$, in Fig. 25 at $z = z_2$ and $t = t_2$ and in Fig. 26 at $z = z_2$ and $t = t_3$. All these curves demonstrate the appropriate convergence of the numerical solution for the more refined grids 1024 × 2048 and 2048 × 4096, which validate the present algorithm for the interaction of an acoustical wave
<figure><img src="image_28.png" /><figcaption>Fig. 24. Pressure and z-velocity component profiles along r-axis at z = z₁ and t = t₁ for five different grids.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 25. Pressure and z-velocity component profiles along r-axis at z = z₂ and t = t₂ for five different grids.</figcaption></figure>

<caption>**Table 6**</caption>

Error analysis on the pressure at the solid frontier $p_r$ at z = z₁ and t = t₁ and at z = z₂ and t = t₂ from computations on five different grids.

<table><thead><tr><th>Mesh</th><th>p<sub>r</sub>(z<sub>1</sub>, t<sub>1</sub>) (Pa)</th><th>δp<sub>r</sub>(z<sub>1</sub>, t<sub>1</sub>) (Pa)</th><th>order</th><th>p<sub>r</sub>(z<sub>2</sub>, t<sub>2</sub>) (Pa)</th><th>δp<sub>r</sub>(z<sub>2</sub>, t<sub>2</sub>) (Pa)</th><th>order</th></tr></thead><tbody><tr><td>128 × 256</td><td>101402.5</td><td>21.7</td><td>-</td><td>101384</td><td>19.0</td><td>-</td></tr><tr><td>256 × 512</td><td>101412.0</td><td>12.2</td><td>0.83</td><td>101392.5</td><td>10.5</td><td>0.86</td></tr><tr><td>512 × 1024</td><td>101420.3</td><td>3.9</td><td>1.64</td><td>101399.5</td><td>3.5</td><td>1.58</td></tr><tr><td>1024 × 2048</td><td>101423.3</td><td>0.9</td><td>2.11</td><td>101402</td><td>1.0</td><td>1.8</td></tr><tr><td>2048 × 4096</td><td>101424.2</td><td>-</td><td>reference</td><td>101403</td><td>-</td><td>reference</td></tr></tbody></table>

with a curved solid surface. In particular, in Fig. 24, the acoustic boundary layer around the cylinder can be clearly observed on the z-velocity component profile along the radial direction with a measured boundary layer thickness $\eta \sim 13$ µm which is consistent with the theoretical estimation given by Eq. (84). As expected, only the two more refined grids enable an accurate computation of this boundary layer. Other snapshots confirm this trend with a satisfactory space convergence, which is clearly observed on all cases with the two more refined grids.

To investigate more thoroughly on the accuracy order, we present an error analysis in Table 6. As previously done, the solution on the most refined grid is used as a reference solution to determine the accuracy order. These analyses are performed on the pressure
<figure><img src="image_30.png" /><figcaption>Fig. 26. Pressure and z-velocity component profiles along r-axis at z = z₂ and t = t₃ for five different grids.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 27. Schematic representation of the axisymmetric domain. The outer walls are isolated thermally and electrically. The imposed electric flux is represented on the zone 1. The constant electric potential which absorb the current is the zone 2. All the other surfaces are isolated electrically.</figcaption></figure>

at the solid frontier $p_{\Gamma}$ at two points $z = z_1$ and $z = z_2$, and at times $t = t_1$ and $t = t_2$ for which the acoustical wave is interacting with the solid frontier at the considered point. For this purpose we define the error, as $\delta p_{\Gamma}(z,t) = |p_{\Gamma}^{ref}(z,t) - p_{\Gamma}(z,t)|$, with $p_{\Gamma}^{ref}(z,t)$ the reference solution obtained on the thinnest grid 2048 × 4096. Table of errors show a convergence rate around 1 for the coarsest grid which next tends to 2 for higher resolution. Such evolution of the accuracy order for increasing resolution can be explained by considering that a computation on a too coarse grid cannot capture the formation of the acoustical boundary layer, and only first order can be obtained in a configuration where the grid spacing is larger than the typical spatial scale of the phenomenon. On the other hand, on a more refined grid, this acoustical boundary layer can be fully resolved, and then the correct order of the numerical solver can be recovered.

### 4.7. Electric arc

To highlight the interest of the proposed solver for thermal plasma modeling, preliminary results on an arc electric discharge between two electrodes are presented in this section. We consider an axisymmetric computational domain represented in Fig. 27, on which one can visualize the cylindrical electrodes and the related boundary conditions described with the immersed boundary method. $R_{elec} = 1.5$ mm is the radius of the electrodes and the space between the two electrodes is $H_{elec} = 3$ mm. The lengths of the computational domain are $l_r = 0.01$ m and $l_z = 0.02$ m, and the computational grid contains 256 × 512 points, which corresponds to a
<figure><img src="image_32.png" /><figcaption>(a) Temperature field and streamlines of the velocity at t = 20 µs (b) Electric potential and streamlines of the electric current at t = 20 µs</figcaption></figure>

<caption>Fig. 28. Closer view of the computed solution between the electrodes.</caption>

cell size ~ 50 µm. As a first simplification, we simply apply Neumann boundary conditions for the temperature field on the electrodes boundaries, $\bar{n} \cdot \nabla T = 0$. Indeed, with this benchmark we do not take interest in the electrodes heating by the electrical arc which is a more complex problem, and will be the topic of future works. As it can be visualized in Fig. 27, the electrical voltage that generates the electric arc is computed by imposing a Neumann boundary condition on the electric potential of the upper electrode such as,

$$ \sigma \left. \frac{\partial V}{\partial n} \right|_{\Gamma} = a \left( 1 - \left( \frac{2r}{R_{elec}} \right)^5 \right) \quad \text{if } z = z_1 \text{ and } 0 < r < R_1 \text{ and } \sigma \left. \frac{\partial V}{\partial n} \right|_{\Gamma} = 0 \text{ elsewhere,} \quad (85) $$

with $z_1 = 0.0015$ m, $R_1 = R_{elec}/2$ and $a = 4 \cdot 10^7$ A m$^{-2}$ is chosen so that the electric current $I$ is equal to 50 A. For the bottom electrode, the following conditions are imposed to receive the emitted current from the upper electrode,

$$ V_{\Gamma} = 0 \quad \text{if } z = z_2 \text{ and } 0 < r < R_2 \text{ and } \sigma \left. \frac{\partial V}{\partial n} \right|_{\Gamma} = 0 \quad \text{elsewhere,} \quad (86) $$

with $z_2 = -0.0015$ m and $R_2 = 0.9 R_{elec}$.

The fluid considered is air, and the physical properties are those of an air plasma in the plasma region. Strong variations of the physical properties with the temperature are reported in Appendix B, where tabulated values have been plotted for $C_p$, $C_v$, $\rho$, $k$, $\mu$, $\epsilon$ and $\sigma_e$. These variations of physical properties are taken into account in our numerical approach. The space between the two electrodes is initialized as a hot channel of 7000 K (at atmospheric pressure) at which electrical conductivity is high enough to allow electric current circulation in the domain, as it can be observed in Fig. B.31d. A uniform temperature $T_0 = 300$ K is imposed elsewhere. The temperature field and streamlines of the flow are plotted in Fig. 28a at $t = 20$ µs. One can observe an expanding motion of the plasma and high temperature values varying in the range 12000 – 25000K between the electrodes. In Fig. 28b, the electric potential and the streamlines of the electric current are plotted. This enables to visualize the formation of the electric current between both electrodes due to the voltage induced by the immersed boundary conditions given by Eq. (85) and Eq. (86). In Fig. 28, the propagation of the pressure wave induced during the first time of the electric discharge can be visualized at different times. A larger view of the temperature field at $t = 100$ µs has also been plotted in this figure in order to visualize the plasma expansion related to the electrodes position at a longer time.

## 5. Conclusion

In this article, a well-known IBM, initially developed for incompressible flows in [1], has been adapted to the pressure-based compressible solver recently proposed in [2]. The coupling uses extensively second order numerical schemes developed in anterior works for Immersed Dirichlet Boundary Condition (IDBC) [52] and Immersed Neumann Boundary Condition (INBC) [1]. This provides a consistent framework to impose accurately suitable boundary conditions in computational cells crossed by the solid interface. Several benchmarks involving density variations and an interaction with a complex geometry are proposed and successfully performed. External flows around the immersed solid are considered, as well as internal flows inside the solid. Validations based on space convergence, mass and/or energy conservation and comparisons with reference solutions are proposed. The proposed benchmarks sweep several kind of applications such as free-convection, acoustical wave propagation and a preliminary benchmark of the formation of a thermal plasma between two electrodes during the discharge of an electric arc. This illustrates the potential of the present solver for future works on thermal plasma modeling and electric arc discharge. The developments proposed in this paper will be easily
<figure><img src="image_33.png" /><figcaption>Fig. 29. Acoustic wave on 29a, 29b and 29c. The Temperature field on 29d.</figcaption></figure>

applicable to other pressure based compressible solvers as for instance, in [35,36,79]. As the proposed developments are based on a primitive formulation of the equations, they are restricted to subsonic flows. Future works will be devoted to a generalization of the overall formalism to a conservative formulation of the equation, to be able to perform numerical simulation of supersonic flows. The present work is a step forward that will enable to perform numerical simulations of a plasma bubble formation during an electric arc discharge in a liquid pool. For this purpose, this solver will be coupled in future works to the pressure based compressible solver for two-phase flows with liquid-vapor phase change recently presented in [45].

CRediT authorship contribution statement

**Sergiu Coseru:** Writing – review & editing, Visualization, Validation, Software, Conceptualization. **Sébastien Tanguy:** Writing – original draft, Validation, Supervision, Software, Project administration, Methodology, Funding acquisition, Conceptualization. **Pierre Freton:** Writing – review & editing, Validation, Software, Project administration, Methodology, Funding acquisition. **Jean-Jacques Gonzalez:** Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization. **Annafederica Urbano:** Writing – review & editing, Validation, Software, Conceptualization. **Marie Bibal:** Writing – review & editing, Validation, Software, Conceptualization. **Gauthier Bourdon:** Writing – review & editing, Validation, Software, Conceptualization.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors gratefully acknowledge French region Occitanie and UFTMIP for funding the PhD thesis of Sergiu Coseru.

The authors gratefully acknowledge AID (Agence Innovation Défense) for funding the PhD thesis of Marie Bibal (No 2020213).

The authors acknowledge for the financial support of the French government and BPI France in the framework of France Relance 2020-2023 for the project French Fab n° DOS0152116 in collaboration with Framatome for funding the PhD thesis of Gauthier Bourdon.

This work was granted access to the HPC resources of Calmip under the allocation P23023.

Appendix A. Cubic equations of state

A cubic EoS can be expressed with a cubic equation for the density $\rho$:

$$ \rho^3 + a_1 \rho^2 + a_2 \rho + a_3 = 0, \qquad (A.1) $$

where the $a_i$ coefficients are a function of $T$, $p$ and depend on the fluid. Eq. (A.1) is commonly expressed under the following generalized form:

$$ p = \frac{\rho RT}{1 - \rho b} - \frac{a \rho^2}{(1 + b_1 \rho)(1 + b_2 \rho)}, \qquad (A.2) $$

where $a$, $b$, $b_1$ and $b_2$ are related with the $a_i$ coefficients of Eq. (A.1) by the following relations:

$$ a_1 = - \frac{RT(b_1 + b_2) - a + (b(b_1 + b_2) - b_1 b_2)p}{b_1 b_2 RT + (a + p b_1 b_2)b}, \qquad (A.3) $$

$$ a_2 = \frac{RT + (b - b_1 - b_2)p}{b_1 b_2 RT + (a + p b_1 b_2)b}, \qquad (A.4) $$

$$ a_3 = \frac{-p}{b_1 b_2 RT + (a + p b_1 b_2)b}, \qquad (A.5) $$

where $R = R/W$ is the gas constant, with $R = 8314.46$ J/kg/K universal gas constant and $W$ molecular weight. Different cubic EoS are characterized by different coefficients $a_i$. Specifically, the coefficients $b_1$ and $b_2$ are reported in Table A.7 for the different EoS as function of the coefficients $a$, $b$ and $c$ that are defined as:

$$ a = \Omega_a \frac{R^2 T_c^2}{p_c} \alpha(T), \qquad (A.6) $$

$$ b = \Omega_b \frac{RT_c}{p_c}, \qquad (A.7) $$

$$ c = \Omega_c \frac{RT_c}{p_c}, \qquad (A.8) $$

where $\alpha(T)$ is the isochoric heat capacity correction factor, $p_c$ and $T_c$ are the critical pressure and temperature. The parameters $\Omega_a$, $\Omega_b$ and $\Omega_c$ are reported in Table A.8. For the PTV EoS, these are based on the compressibility factor $Z_c$ at the critical point:

$$ Z_c = \frac{p_c V_c}{RT_c}, \qquad (A.9) $$

<caption>Table A.7</caption>

$b_1$ and $b_2$ coefficients of Eq. (A.2) for different cubic EoS: Van der Waals (VdW), Soave-Redlich-Kwong (SRK), Peng-Robinson (PR) and Patel-Teja-Valderrama (PTV).

<table><thead><tr><th>EoS</th><th>b<sub>1</sub></th><th>b<sub>2</sub></th></tr></thead><tbody><tr><td>VdW</td><td>0</td><td>0</td></tr><tr><td>SRK</td><td>0</td><td>b</td></tr><tr><td>PR</td><td>b(1 - &radic;2)</td><td>b(1 + &radic;2)</td></tr><tr><td>PTV</td><td>b/2 (1 + c/b + &radic;(1 + 6c/b + c<sup>2</sup>/b))</td><td>b/2 (1 + c/b - &radic;(1 + 6c/b + c<sup>2</sup>/b))</td></tr></tbody></table>
<caption>**Table A.8**</caption>

$\Omega_a$, $\Omega_b$ and $\Omega_c$ coefficients for the computations of the coefficients Eqs. (A.6)-(A.8) for different cubic EoS: Van der Waals (VdW), Soave-Redlich-Kwong (SRK), Peng-Robinson (PR) and Patel-Teja-Valderrama (PTV).

<table><caption>**Table A.9**</caption><thead><tr><th>EoS</th><th>Ω<sub>a</sub></th><th>Ω<sub>b</sub></th><th>Ω<sub>c</sub></th></tr></thead><tbody><tr><td>VdW</td><td>27</td><td>1</td><td>0</td></tr><tr><td>SRK</td><td>64</td><td>8</td><td>0</td></tr><tr><td>PR</td><td>0.42748</td><td>0.08664</td><td>0</td></tr><tr><td>PTV</td><td>0.45724</td><td>0.07780</td><td>0</td></tr><tr><td></td><td>0.66121 – 0.76105Z<sub>c</sub></td><td>0.02207 + 0.20868Z<sub>c</sub></td><td>0.57765 – 1.87080Z<sub>c</sub></td></tr></tbody></table>

$m_0$, $m_1$ and $m_2$ coefficients used in Eqs. (A.11) and (A.12) for different cubic EoS: Van der Waals (VdW), Soave-Redlich-Kwong (SRK), Peng-Robinson (PR) and Patel-Teja-Valderrama (PTV).

<table><thead><tr><th>EoS</th><th>m<sub>0</sub></th><th>m<sub>1</sub></th><th>m<sub>2</sub></th></tr></thead><tbody><tr><td>VdW</td><td>0</td><td>0</td><td>0</td></tr><tr><td>SRK</td><td>0.488</td><td>1.574</td><td>-0.176</td></tr><tr><td>PR</td><td>0.37464</td><td>1.54226</td><td>-0.26992</td></tr><tr><td>PTV</td><td>0.46283</td><td>3.58230</td><td>8.19417</td></tr></tbody></table>

with $V_c$ molar volume at the critical point.

The capacity correction factor $\alpha(T)$ is given by:

$$ \alpha(T) = \left[ 1 + m \left( 1 - \sqrt{\frac{T}{T_c}} \right) \right]^2, \qquad (A.10) $$

where the parameter $m$ depends on the EoS and is given by:

$$ m = m_0 + m_1(\omega Z_c) + m_2(\omega Z_c)^2, \qquad (A.11) $$

for PTV and by:

$$ m = m_0 + m_1 \omega + m_2 \omega^2, \qquad (A.12) $$

for all the other EoS. The $m_i$ coefficients are reported in Table A.9. Finally, $\omega$ is the acentric factor:

$$ \omega = -10 \log \left( \frac{p_{sat}}{p_c} \right), \qquad (A.13) $$

at $T = 0.7T_c$.

In order to evaluate the density $\rho$ for a given $p$ and $T$, both for single phase and two phase regions, the roots of the third order polynomial Eq. (A.1) have to be evaluated. This can be done analytically, for instance using the Cardan method (see [2] for more details). On the other hand, the equation is linear with respect to the temperature which can be directly evaluated from $\rho$ and $p$ from Eq. (A.2):

$$ T = \frac{a - b\rho}{\rho R} p + \frac{a\rho(1 - b\rho)}{R(1 + b_1\rho)(1 + b_2\rho)}. \qquad (A.14) $$

Other thermodynamic variables are required to close the system of Eqs. (1)-(4) as $c_\rho$, $c^2$ and $\alpha$ which can be computed in each phase from the EoS. In particular, $c_\rho$ can be deduced from the specific heat at constant volume $c_v$ with the following relation [2]:

$$ c_{\rho} = c_{v} + \frac{T\alpha^{2}}{\rho\beta}, \qquad (A.15) $$

where $\alpha$ and $\beta$ are the isobaric and isothermal expansion coefficients:

$$ \alpha = -\frac{1}{\rho} \left( \frac{\partial p}{\partial T} \right)_\rho, \qquad \beta = \frac{1}{\rho} \left( \frac{\partial p}{\partial \rho} \right)_T. \qquad (A.16) $$
Note that within the framework of cubic equations of state $c_v$ does not depend on the pressure, and it is a function of temperature only. In the present study for which temperature variations are in the range of a few Kelvin, the $c_v$ is kept constant during the simulation and equal to the ideal gas value at the average temperature. From the generalized EoS form of Eq. (A.2) the pressure thermodynamic derivatives are given by:

$$ \left(\frac{\partial p}{\partial T}\right)_{\rho} = \frac{\rho R}{1 - b\rho} - \frac{a' \rho^2}{(1 + b_1 \rho)(1 + b_2 \rho)}, \qquad (A.17) $$

$$ \left(\frac{\partial p}{\partial \rho}\right)_T = \frac{\frac{RT}{(b\rho-1)^2} - a\rho(\rho(b_1+b_2)+2)}{[(1+b_1\rho)(1+b_2\rho)]^2}, \qquad (A.18) $$

where $a' = \frac{da}{dT}$ is the derivative of $a$ with respect to temperature $T$ (see Eq. (A.10)).

## Appendix B. Thermodynamic and transport coefficients of a plasma

The plasma thermodynamic properties like the specific heats $c_p$, $c_v$, and the density $\rho$ are computed using the partition function. As for the transport properties like the thermal conductivity $k$, electrical conductivity $\sigma_e$, the dynamic viscosity $\mu$, they are computed using the Chapman-Enskog method from the Boltzmann equation. The thermophysical variations of these properties are plotted in this appendix for an air plasma. These tabulated values have been included in the proposed plasma solver.

### Data availability

Data will be made available on request.

<figure><img src="image_34.png" /><figcaption>Fig. B.30. The graphs of Thermodynamic properties of air. The red line for the density graph represents the perfect gas law, and the black line is the computed plasma density.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Fig. B.31. The graphs of transport properties of air.</figcaption></figure>

References

[1] Y.T. Ng, C. Min, F. Gibou, An efficient fluid-solid coupling algorithm for single-phase flows, J. Comput. Phys. 228 (2009) 8807-8829.

[2] A. Urbano, M. Bibal, S. Tanguy, A semi implicit compressible solver for two-phase flows of real fluids, J. Comput. Phys. 456 (2022) 111034.

[3] C.S. Peskin, Flow patterns around heart valves: a numerical method, J. Comput. Phys. 10 (1972) 252-271.

[4] E.Y. Tau, A second-order projection method for the incompressible Navier-Stokes equations in arbitrary domains, J. Comput. Phys. 115 (1994) 147-152.

[5] A. Almgren, J.B. Bell, P. Colella, A Cartesian grid projection method for the incompressible Euler equations in complex geometries, SIAM J. Sci. Comput. 18 (5) (1997) 1289-1309.

[6] H. Udaykumar, R. Mittal, P. Rampumggoon, A. Khanna, A sharp interface Cartesian grid method for simulating flows with complex moving boundaries, J. Comput. Phys. 174 (2001) 345-380.

[7] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479-517.

[8] R. Mittal, I. Gianluca, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239-261.

[9] M. Ulhmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209 (2005) 448-476.

[10] R. Mittal, H. Dong, M. Bozkurttas, F. Najjar, A. Vargas, A. von Loebbecke, A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries, J. Comput. Phys. 227 (2008) 4825-4852.

[11] S. Popinet, Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries, J. Comput. Phys. 190 (2003) 572-600.

[12] A. Guittet, M. Theillard, F. Gibou, A stable projection method for the incompressible Navier-Stokes equations on arbitrary geometries and adaptive quad/octrees, J. Comput. Phys. 292 (2015) 215-238.

[13] C. Selçuk, A.R. Ghigo, S. Popinet, A. Wachs, A fictitious domain method with distributed Lagrange multipliers on adaptive quad/octrees for the direct numerical simulation of particle-laden flows, J. Comput. Phys. 430 (2021) 109954.

[14] R. Egan, A. Guittet, F. Temprano-Coleto, T. Isaac, F.J. Peaudecerf, J.R. Landel, P. Luzzatto-Fegiz, C. Burstedde, F. Gibou, Direct numerical simulation of incompressible flows on parallel octree grids, J. Comput. Phys. 428 (2021) 110084.

[15] R.B. Pember, J.B. Bell, P. Colella, W.Y. Crutchfield, M.L. Welcome, An adaptive Cartesian grid method for unsteady compressible flow in irregular dregions, J. Comput. Phys. 120 (1995) 278-304.

[16] P. Colella, D.P. Graves, B.J. Keen, D. Modiano, A Cartesian grid embedded boundary method for hyperbolic conservation laws, J. Comput. Phys. 211 (2006) 347-366.

[17] R. Ghias, R. Mittal, H. Dong, A sharp interface immersed boundary method for compressible viscous flows, J. Comput. Phys. 225 (2007) 528-553.

[18] M. de Tullio, P. De Palma, G. Iaccarino, G. Pascazio, M. Napolitano, An immersed boundary method for compressible flows using local grid refinement, J. Comput. Phys. 225 (2007) 2098-2117.

[19] F. Örley, V. Pasquariello, S. Hickel, N.A. Adams, Cut-element based immersed boundary method for moving geometries in compressible liquid flows with cavitation, J. Comput. Phys. 283 (2015) 1-22.

[20] B. Muralidharan, S. Menon, A high-order adaptive Cartesian cut-cell method for simulation of compressible viscous flow over immersed bodies, J. Comput. Phys. 321 (2016) 342-368.

[21] Y. Qu, R. Shi, R.C. Batra, An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids, J. Comput. Phys. 354 (2018) 672-691.
[22] N. Jianguo, Z. Hetao, X. Xiangzhao, M. Tianbao, A novel fluid-structure interaction algorithm for compressible flows and deformable structures, J. Comput. Phys. 426 (2021) 109921.

[23] M. Natarajan, R. Grout, W. Zhang, M. Day, A moving embedded boundary approach for the compressible Navier-Stokes equations in a block-structured adaptive refinement framework, J. Comput. Phys. 465 (2022) 111315.

[24] G. Shallcross, J. Capecelatro, An explicit characteristic-based immersed boundary method for compressible flows, J. Comput. Phys. 449 (2022) 110804.

[25] H. Yu, M. Pantano, An immersed boundary method with implicit body force for compressible viscous flow, J. Comput. Phys. 459 (2022) 111125.

[26] L. Ménez, P. Parnaudeau, M. Beringhier, E. Goncalves Da Silva, Assessment of volume penalization and immersed boundary methods for compressible flows with various thermal boundary conditions, J. Comput. Phys. 493 (2023) 112465.

[27] M. Berger, M. Aftosmis, Progress towards a Cartesian cut-cell method for viscous compressible flow, in: 50th AIAA Aerospace Sciences Meeting, vol. 131, 2012.

[28] G. Huber, S. Tanguy, J.-C. Béra, B. Gilles, A time splitting projection scheme for compressible two-phase flows. Application to the interaction of bubbles with ultrasound waves, J. Comput. Phys. 302 (2015) 439–468.

[29] Z. Ou, C. Chi, L. Guo, D. Thevenin, A directional ghost-cell immersed boundary method for low Mach number reacting flows with interphase heat and mass transfer, J. Comput. Phys. 468 (2022) 11447.

[30] N. Kwatra, J. Su, J.T. Grétarsson, R. Fedkiw, A method for avoiding the acoustic time step restriction in compressible flow, J. Comput. Phys. 228 (2009) 4146–4161.

[31] F. Xiao, R. Akoh, S. Ii, Unified formulation for compressible and incompressible flows by using multi-integrated moments II: multi-dimensional version for compressible and incompressible flows, J. Comput. Phys. 213 (2006) 31–56.

[32] J.-P. Caltagirone, S. Vincent, C. Caruyer, A multiphase compressible model for the simulation of multiphase flows, Comput. Fluids 50 (2011) 24–34.

[33] M. Jemison, M. Sussman, M. Arienti, Compressible, multiphase semi-implicit method with moment of fluid interface representation, J. Comput. Phys. 279 (2014) 182–217.

[34] M. Arienti, M. Sussman, An embedded level set method for sharp-interface multiphase simulations of diesel injectors, Int. J. Multiph. Flow 59 (2014) 1–14.

[35] M. Arienti, M. Sussman, A numerical study of the thermal transient in high-pressure diesel injection, Int. J. Multiph. Flow 88 (2017) 205–221.

[36] D. Fuster, S. Popinet, An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension, J. Comput. Phys. 374 (2018) 752–768.

[37] M.B. Kuhn, O. Desjardins, An all-Mach, low-dissipation strategy for simulating multiphase flows, J. Comput. Phys. 445 (2021) 110602.

[38] Y. Saade, D. Lhose, D. Fuster, A multigrid solver for the coupled pressure-temperature equations in an all-Mach solver with vof, J. Comput. Phys. 476 (2023) 11865.

[39] S.F. Golovashchenko, A.J. Gillard, A.V. Mamutov, R. Ibrahim, Pulsed electrohydraulic springback calibration of parts stamped from advanced high strength steel, J. Mater. Process. Technol. 214 (2014) 2796–2810.

[40] M. Zia, A. Fazli, M. Soltanpour, Warm electrohydraulic forming: a novel high speed forming process, Proc. Eng. 207 (2017) 323–328.

[41] Q. Yu, H. Zhang, R. Yang, Z. Cai, K. Liu, Effects of confining pressure and hydrostatic pressure on the fracturing of rock under cyclic electrohydraulic shock waves, Energies 15 (16) (2022) 6032, Multidisciplinary Digital Publishing Institute.

[42] F. Lago, J.-J. Gonzalez, P. Freton, A. Gleizes, A numerical modelling of an electric arcand its interaction with the anode: Part I. The two-dimensional model, J. Phys. D, Appl. Phys. 37 (2004) 883–897.

[43] J.-J. Gonzalez, F. Lago, M. Masquère, X. Francies, Numerical modelling of an electric arc and its interaction with the anode: Part II. The three-dimensional model-influence of external forces on the arc column, J. Phys. D, Appl. Phys. 38 (2005) 306–318.

[44] A. Gleizes, J.-J. Gonzalez, P. Freton, Thermal plasma modelling, J. Phys. D, Appl. Phys. 38 (2005) R153–R183.

[45] M. Bibal, M. Deferrez, S. Tanguy, A. Urbano, A compressible solver for two phase-flows with liquid-vapour phase change for bubble cavitation, J. Comput. Phys. 500 (2024) 112750.

[46] H. Terashima, G. Tryggvason, A front-tracking/ghost-fluid method for fluid interfaces in compressible flows, J. Comput. Phys. 228 (2009) 4012–4037.

[47] H. Solo, M. Benmouffok, P. Freton, J.-J. Gonzalez, Stoichiometry air - ch4 mixture: composition, thermodynamic properties and transport coefficients, Plasmas Phys. Technol. 7 (2020) 21–29.

[48] H. Solo, P. Freton, J.-J. Gonzalez, Compositions chimiques et propriétés thermodynamiques à l'ét l'un mélange air-ch4, J. Int. Technol, Innov. Phys. Energ. Environ. 5 (2) (2019) 21–29.

[49] T. Billoux, Y. Cressault, P. Teulet, A. Gleizes, Calculation of the net emission coefficient of an air thermal plasma at very high pressure, J. Phys. Conf. Ser. 406 (2012) 012010.

[50] S. Vincent, J. Brandle de Motta, A. Sarthou, J.-L. Estivalezes, O. Simonin, E. Climent, A Lagrangian vof tensorial penalty method for the dns of resolved particle-laden flows, J. Comput. Phys. 256 (2014) 582–614.

[51] Y.-H. Tseng, J.H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, J. Comput. Phys. 192 (2003) 593–623.

[52] F. Gibou, R. Fedkiw, L.-T. Cheng, M. Kang, A second-order-accurate symmetric discretization of the Poisson equation on irregular domains, J. Comput. Phys. 176 (2002) 205–227.

[53] J. Papac, F. Gibou, C. Ratsch, Efficient symmetric discretization for the Poisson, heat and Stefan-type problems with Robin boundary conditions, J. Comput. Phys. 229 (3) (2010) 875–889.

[54] L. Rueda Villegas, R. Alis, M. Lepilliez, S. Tanguy, A ghost fluid/level set method for boiling flows and liquid evaporation: application to the Leidenfrost effect, J. Comput. Phys. 316 (2016) 789–813.

[55] G. Mialhe, S. Tanguy, L. Tranier, E.-R. Popescu, D. Legendre, An extended model for the direct numerical simulation of droplet evaporation. Influence of the Marangoni convection on Leidenfrost droplet, J. Comput. Phys. (2023) 112366.

[56] A. Urbano, S. Tanguy, C. Colin, Direct numerical simulation of nucleate boiling in zero gravity conditions, Int. J. Heat Mass Transf. 143 (2019) 118521.

[57] L. Torres, A. Urbano, C. Colin, S. Tanguy, On the coupling between direct numerical simulation of nucleate boiling and a micro-region model at the contact line, J. Comput. Phys. 497 (2024) 112602.

[58] M. Lepilliez, E.-R. Popescu, F. Gibou, S. Tanguy, On two-phase flow solvers in irregular domains with contact line, J. Comput. Phys. 321 (2016) 1217–1251.

[59] E.Y. Tau, A Cartesian grid embedded boundary method for Poisson's equation on irregular domains, J. Comput. Phys. 147 (1998) 60–85.

[60] P. McCorquodale, P. Colella, H. Johansen, A Cartesian grid embedded boundary method for the heat equation on irregular domains, J. Comput. Phys. 173 (2001) 620–635.

[61] P. Schwartz, M. Barad, P. Colella, T. Ligocki, A Cartesian grid embedded boundary method for the heat equation and Poissons equation in three dimensions, J. Comput. Phys. 211 (2006) 531–550.

[62] C. Min, F. Gibou, Geometric integration over irregular domains with application to level-set methods, J. Comput. Phys. 226 (2) (2007) 1432–1443.

[63] C. Min, F. Gibou, Robust second-order accurate discretizations of the multi-dimensional Heaviside and Dirac delta functions, J. Comput. Phys. 227 (22) (2008) 9686–9695.

[64] A. Dalmon, M. Lepilliez, S. Tanguy, A. Pedrono, B. Busset, H. Bavestrello, J. Mignot, Direct numerical simulation of a bubble motion in a spherical tank under external forces and microgravity conditions, J. Fluid Mech. 849 (2018) 467–497.

[65] A. Dalmon, M. Lepilliez, S. Tanguy, R. Alis, E.-R. Popescu, R. Roumigué, T. Miquel, B. Busset, H. Bavestrello, J. Mignot, Comparison between the fluidics experiment and direct numerical simulations of fluid sloshing in spherical tanks under microgravity conditions, Microgravity Sci. Technol. 31 (1) (2019) 123–138.

[66] R. Borges, M. Carmona, B. Costa, W.S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (6) (2008) 3191–3211.
[67] J. Dendy, Black box multigrid, J. Comput. Phys. 48 (1982) 366-386.

[68] S. MacLachlan, J. Tang, C. Vuik, Fast and robust solvers for pressure-correction in bubbly flow problems, J. Comput. Phys. 227 (2008) 9742-9761.

[69] T.D. Aslam, A partial differential equation approach to multidimensional extrapolation, J. Comput. Phys. 193 (1) (2004) 349-355.

[70] F. Gibou, C. Min, R. Fedkiw, High resolution sharp computational methods for elliptic and parabolic problems in complex geometries, J. Sci. Comput. 54 (2) (2013) 369-413.

[71] S. Tanguy, M. Sagan, B. Lalanne, F. Couderc, C. Colin, Benchmarks and numerical methods for the simulation of boiling flows, J. Comput. Phys. 264 (2014) 1-22.

[72] A. Dalmon, K. Kenteswaran, G. Mialhe, B. Lalanne, S. Tanguy, Fluids-membrane interaction with a full Eulerian approach based on the level set method, J. Comput. Phys. 406 (2020) 109171.

[73] R. Egan, F. Gibou, xgfm: recovering convergence of fluxes in the ghost fluid method, J. Comput. Phys. 409 (2020) 109351.

[74] G. Cesini, M. Paroncini, G. Cortella, M. Manzan, Natural convection from a horizontal cylinder in a rectangular cavity, Int. J. Heat Mass Transf. 42 (10) (1999) 1801-1811.

[75] K. Karagözis, R. Kamakoti, C. Pantano, A low numerical dissipation immersed interface method for the compressible Navier-Stokes equations, J. Comput. Phys. 229 (2010) 701-727.

[76] D. Canuto, K. Taira, Two-dimensional compressible viscous flow around a circular cylinder, J. Fluid Mech. 785 (2015) 349-371.

[77] M. Berggren, A. Bernland, D. Noreland, Acoustic boundary layers as boundary conditions, J. Comput. Phys. 371 (2018) 633-650.

[78] G.P. Ward, R.K. Lovelock, A.R.J. Murray, A.P. Hibbins, J.R. Sambles, J.D. Smith, Boundary-layer effects on acoustic transmission through narrow slit cavities, Phys. Rev. Lett. 115 (2015) 044302.

[79] L. Germes Martinez, B. Duret, J. Reveillon, F.-X. Demoulin, A new dns formalism dedicated to turbulent two-phase flows with phase change, Int. J. Multiph. Flow 143 (2021) 103762.