# Time integration of cut-cell surfaces for geometric conservation, applied to scalar transport with moving interfaces

Dorian Dupuy, Adrien Toutant, Antoine Bittard Du Cluzeau, Guillaume Bois

► To cite this version:

Dorian Dupuy, Adrien Toutant, Antoine Bittard Du Cluzeau, Guillaume Bois. Time integration of cut-cell surfaces for geometric conservation, applied to scalar transport with moving interfaces. Computers and Fluids, 2025, 295, pp.106627. 10.1016/j.compfluid.2025.106627. cea-05335832

HAL Id: cea-05335832

https://cea.hal.science/cea-05335832v1

Submitted on 28 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

# Time integration of cut-cell surfaces for geometric conservation, applied to scalar transport with moving interfaces

D. Dupuy<sup>a,b, ID,*</sup>, A. Toutant<sup>b, ID</sup>, A. du Cluzeau<sup>a</sup>, G. Bois<sup>a, ID</sup>

## ARTICLE INFO

**Keywords:**

- Two-phase flow
- Interface-resolved
- Cut-cell method
- Heat transfer

## ABSTRACT

In cut-cell methods, interface-resolved simulations of two-phase flows are performed by cutting a fixed nonconforming mesh at the interface boundary. The cells which are cut to conform to the interface use modified discretisation schemes that account for the modified cell volume and face areas of cut cells, which evolve dynamically with the motion of the interface. This article investigates the effect of the method used for time integration of cut-face areas in a cut-cell method, for the convection-diffusion of a passive scalar in a two-phase flow with moving interfaces. The cut-cell method, based on a finite-volume approach and a three-dimensional staggered Cartesian grid, naturally enforces strict conservation laws and ensures numerical stability in small cells using a flux-redistribution strategy. The simulation of heat diffusion in and around a spherical interface under a uniform velocity field is addressed. A semi-implicit time-integration method taking into account initial and final cut-face areas provides significant improvements at a negligible cost compared to an explicit time-integration method.

## 1. Introduction

This paper presents a method to numerically simulate the transfer of heat in a two-phase flow, such as the flow of two immiscible fluids or around solid particles. The simulation of two-phase flows may be addressed in various modelling frames. This includes e.g. the Eulerian–Eulerian framework of the interpenetrating continua and Eulerian–Lagrangian approaches [1–4]. The present study is concerned with simulations of two-phase flows where the two phases are modelled as a continuous media bound together by an interface. The numerical methods that resolves the interface between the phases are often classified based on whether the computational mesh conforms the interface boundary. The use of interface-conforming meshes is best exemplified by arbitrary Lagrangian–Eulerian methods [5–10], in which mesh nodes may be allowed to move to coincide with the interface. This enables the numerical scheme to prescribe the coupling between the phases as an interface-boundary condition (jump conditions at the interface). However, the intricacies and costs associated with moving-mesh methods [11] have motivated the development of methods using a fixed nonconforming mesh, typically a Cartesian grid. These methods have been collectively referred to as Cartesian-grid

methods, embedded-interface methods or immersed-boundary methods in the context of rigid bodies [12,13].

Since the mesh does not conform to the interface boundary, a different strategy must then be employed to represent the interface. Front-tracking methods track the interface explicitly using a set of Lagrangian markers [14–16]. Front-capturing methods track the interface implicitly using an intermediate variable. The most common front-capturing methods are the volume-of-fluid [17–19], moment-of-fluid [20,21] level-set [22–25], constrained-interpolation-profile [26] and phase-field [27–30] methods. In addition, nonconforming-mesh methods must apply a particular technique to enforce the interface-boundary condition on nonconforming mesh cells. This is most commonly performed within a one-field formulation, that is by modelling the two phases using a single flow field with effective properties varying rapidly at the interface [31]. The continuity of interfacial fluxes is implicitly encoded in the governing equations whereas discontinuities are handled by specific source terms. The process effectively regularises the discontinuity in fluid property between the phases, which introduces first-order errors in the solution [32]. Its relevance has been demonstrated by various authors [24,33–36]. Ghost cells may also be employed to enforce the interfacial condition on a nonconforming mesh

<footnote><sup>a</sup> Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides (STMF), CEA Centre de Saclay, Gif-sur-Yvette, 91191, France</footnote>

<footnote><sup>b</sup> PROMES-CNRS, Université de Perpignan Via Domitia, Rambla de la thermodynamique, Tecnosud, 66100 Perpignan, France</footnote>

<footnote>* Corresponding author at: Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides (STMF), CEA Centre de Saclay, Gif-sur-Yvette, 91191, France.</footnote>

<footnote>E-mail addresses: dorian.dupuy@cea.fr (D. Dupuy), adrien.toutant@promes.cnrs.fr (A. Toutant), Antoine.DUCLUZEAU@cea.fr (A. du Cluzeau), Guillaume.BOIS@cea.fr (G. Bois).</footnote>
while preserving a strict two-field formulation [37–42]. This entails
prescribing a value for phase in cells beyond the interface boundary,
typically using an interpolation scheme that takes into account the
interfacial condition.

Finally, the fluxes at the interface may be imposed using a cut-
cell method, wherein mesh cells are cut to conform to the interface
boundary [43–45]. Although typically classified as a nonconforming-
mesh method, the technique may also be viewed as a particular type of
interface-conforming arbitrary Lagrangian–Eulerian method, that only
modifies the computational cells intersecting the interface [46]. Indeed,
the cut-cell procedure essentially transforms a fixed nonconforming
mesh into a mesh that dynamically conforms to the interface while
preserving the underlying grid structure of the mesh in uncut cells. As a
result, the method can prescribe an interfacial condition on the cut-cell
faces as an interface-boundary condition (jump conditions at the inter-
face) and naturally enforce strict conservation laws. A distinctive fea-
ture of cut-cell methods is the creation of cut cells with arbitrarily small
volumes, which require special handling to avoid restraining the time
step when using an explicit time-integration scheme. The cell-merging
method combines the smaller cells with a neighbour for time ad-
vancement, separating them subsequently to obtain the time-advanced
field on the original cut mesh [47–52]. The cell-linking method is
conceptually similar, but preserves the original cells, allowing simpler
data structures [53–55]. Also related, the state-redistribution method
performs the time advance with the original grid, merging smaller
cells as a post-processing step to stabilise the solution [56–58]. The h-
box method uses embedded boxes extending away from the interface
boundary to extend the domain of dependence of the fluxes [59–62].
The LS-STAG family of methods use an energy-preserving discretisation
scheme to improve the stability of the convection operator, while
diffusion is often treated implicitly [63,64]. The domain-of-dependence
method introduces a penalty term to restore the proper domain of
dependence of the outflow neighbours of cut cells [65,66]. The flux-
redistribution method uses a stable algorithm to advance the solution in
small cells, redistributing any variation as a flux to ensure conservation
properties [67–73]. Finally, it is also possible to use an implicit time-
integration scheme to ensure the stability of small cells, as proposed
by May and co-authors [74,75] as well as Xie and co-authors [76,77].
Cut-cell methods have mostly been applied to simulations of fluid flows
with fixed [78–80] or moving [81–85] rigid boundaries. Simulations of
fluid-fluid flows with a cut-cell method are most commonly performed
with a level-set method to represent the interface. For example, Bai and
co-authors [86–88] combined the level-set method with a real-ghost
mixing method, based on the mixing procedure of Hu et al. [89], to
handle small cells. Shen et al. [90] and Pang et al. [91,92], building
upon Lin et al. [93], combined the level-set method instead with a cell-
merging method. A similar strategy has also been used in Deng and Li
[94] and Chang et al. [95].

With either a rigid interface or a fluid-fluid interface, the area of cut
Cartesian faces varies rapidly with the displacement of the interface.
In fact, the area of a cut face can vary from zero to its maximal value
within a single convective time step of the simulation. Since cut-surface
areas vary more quickly than physical quantities such as velocity
or temperature, an accurate temporal integration of cut surfaces is
required to estimate accurately the fluxes passing through them over
time. Moreover, the temporal integration of cut surfaces plays a role in
the consistency between the volumetric fluxes passing through the cut
surfaces and the variation of the volume of cut cell, which is known
as the geometric conservation law in arbitrary Lagrangian–Eulerian
methods [96–98]. Various approaches have been used to enforce the
geometric conservation law in arbitrary Lagrangian–Eulerian meth-
ods [99,100]: by revising the volume calculation [96], introducing a
source term in the equations, or by constraining the velocity, normals or
areas of the faces [98,101–103]. In all cases, the geometric conservation
law ensures that the numerical method preserves a uniform flow and
the non-linear stability of its fixed-grid counterpart [104]. The issue of

geometric conservation arises when the computational mesh of the sim-
ulation changes during a simulation time step. Cut-cell methods treat
mesh updates differently from arbitrary Lagrangian–Eulerian methods,
because the mesh motion must exactly adhere to the material motion
on the interface, while remaining fixed elsewhere in the domain. As
a result, cut-cell meshes do not deform smoothly: they may present
large distortion and frequently change topology. In addition, the face
normals are fixed in cut-cell methods and the temporal evolution
of the face areas is not continuously differentiable. Finally, an issue
closely related to the temporal integration of cut surfaces is found
in geometric volume-of-fluid methods [105], which use geometrical
operations to approximate the solution of the volume-fraction equation.
A crucial task in geometric volume-of-fluid methods is estimating the
time integral of phase-specific volumetric fluxes, which may involve
determining the evolution of the phase-specific surface area over time.
For instance, the isoAdvector method of Roenby et al. [106] uses a sub-
stepping approach to obtain the temporal evolution within a time step
of the “submerged” face areas. Various other studies have been con-
cerned with the calculation of the geometrical flux volume associated
with the surfaces, in particular in dimensionally unsplit approaches
[107–110]. In the context of cut-cell method, there is to the authors’
knowledge no specific discussion comparing various time-integration
methods of cut-face areas in a given numerical method.

This article examines the effect of the time-integration method
used to determine cut-cell surfaces in a finite-volume cut-cell method
combined with a front-tracking representation of the interface. The
analysis is restricted to the numerical simulation of a scalar trans-
port equation, using heat transport as an example. The two phases
are separated by a moving interface, without phase change nor wall
effects. The numerical simulations are performed using a staggered
grid representation of velocity and temperature. The interface is rep-
resented using a triangular surface mesh [111]. Geometrical properties
of cut cells, including cut surfaces, are calculated based on the in-
tersection of the discretised interface and the Cartesian mesh. The
calculations are generic and could be used for a spherical particle or
a droplet/bubble with complex geometry. The transported variables
are solved on both sides of the interface in a two-field approach. A
physically motivated flux-redistribution method is used in small cells to
maintain numerical stability with an explicit time-integration method.
Three integration methods are compared for the cut Cartesian faces:
an explicit time-integration method, an algebraic semi-implicit method,
and a semi-implicit method that ensures strict geometric conservation.
The relevance of the methods is assessed for the heat diffusion of a
hot spherical particle under a time-dependent spatially uniform veloc-
ity field, using static pure-diffusion for reference, and demonstrated
for different fluid properties, mesh resolutions and time steps of the
simulation.

The paper is organised as follows. The numerical method is pre-
sented in Section 2. The results are discussed in Section 3.

## 2. Numerical method

This article addresses the numerical resolution of scalar transport in
the frame of a finite-volume cut-cell method, with a staggered Cartesian
grid approach. We use heat transport as a practical example to explain
the underlying principles. Section 2.1 provides the exact equations
that underlie the numerical method, whereas Section 2.2 provides the
algorithms and schemes used to numerically solve heat transfers.

### 2.1. Notations and modelling paradigm

Let us consider a system composed of two phases, denoted $\psi_1$ and
$\psi_2$, and separated by an interface $\Gamma(t)$ which moves according to the
velocity field $u(x, t)$. Within a given phase $\psi \in \{\psi_1, \psi_2\}$, the heat con-
ductivity $\lambda(x, t)$ and the isobaric heat capacity per unit volume $\rho c_p(x, t)$
<figure><img src="image_4.png" /><figcaption>Fig. 1. Notations for the volumes and surfaces found in a two-dimensional cut Cartesian cell $\Omega_{i,j}$.</figcaption></figure>

are constant and uniform. Furthermore, the temperature field $T(\mathbf{x}, t)$ is continuous and obeys the heat equation

$$ \frac{\partial \rho c_p T}{\partial t} = -\nabla \cdot \mathcal{C} + \nabla \cdot \mathcal{D}, \quad (1) $$

where $\mathcal{C}(\mathbf{x}, t) = \rho c_p(\mathbf{x}, t)T(\mathbf{x}, t)u(\mathbf{x}, t)$ is the convective flux and $\mathcal{D}(\mathbf{x}, t) = \lambda(\mathbf{x}, t)\nabla T(\mathbf{x}, t)$ the diffusive flux, assumed to be continuous between the phases. We suppose that the velocity field $u(\mathbf{x}, t)$ is divergence-free and continuous, as no phase change occurs.

The numerical resolution of Eq. (1) requires a three-dimensional Cartesian mesh. Let us denote $\Omega_{i,j,k}$ a given cell of this mesh, which may or may not intersect the interface $\Gamma(t)$. The cell $\Omega_{i,j,k}$ is delimited: in the $x$ direction, by the positions $\xi_i = x_i - \Delta\xi_i/2$ and $\xi_{i+1} = x_i + \Delta\xi_i/2$; in the $y$ direction, by the positions $\eta_j = y_j - \Delta\eta_j/2$ and $\eta_{j+1} = y_j + \Delta\eta_j/2$; and in the $z$ direction, by the positions $\zeta_k = z_k - \Delta\zeta_k/2$. The interface $\Gamma(t_n)$ between the two phases is discretised using a triangular surface mesh, using a front-tracking technique [111]. Without loss of generality, we may define the cut cell $\Omega_{i,j,k}^{\psi}(t)$ as the subregion of phase $\psi$ within $\Omega_{i,j,k}$. Indeed, if the Cartesian cell $\Omega_{i,j,k}$ does not intersect the interface $\Gamma(t)$, the subregion $\Omega_{i,j,k}^{\psi}(t)$ would then either vanish or identify with $\Omega_{i,j,k}$. In that context, the numerical method is specified by expressing the time-advance of temperature in the cut cell $\Omega_{i,j,k}^{\psi}(t)$.

Applying the Reynolds transport theorem to the cut cell $\Omega_{i,j,k}^{\psi}(t)$ leads, using Eq. (1), to

$$ \frac{d}{dt} \int_{\Omega_{i,j,k}^{\psi}(t)} \rho c_p(\mathbf{x}, t) T(\mathbf{x}, t) dV = - \int_{\partial \Omega_{i,j,k}^{\psi}(t) \setminus \Gamma(t)} \mathcal{C}^{\psi} \cdot ndS + \int_{\partial \Omega_{i,j,k}^{\psi}(t)} \mathcal{D}^{\psi} \cdot ndS, \quad (2) $$

where $\partial\Omega_{i,j,k}^{\psi}(t)$ is the outer boundary of $\Omega_{i,j,k}^{\psi}(t)$, $n$ is the outward normal to the cell. There is no convective flux on the interface $\Gamma(t)$, as denoted by the set difference $\partial\Omega_{i,j,k}^{\psi}(t) \setminus \Gamma(t)$, since the velocity of the cell boundary and of the flow are equal on the interface. To simplify the discretisation of this equation, let us introduce the following notations for the relevant portion of each face on the outer boundary (Fig. 1): the cut $x$-normal face, $\Sigma_{i,j,k}^{\psi,x}(t) = (\partial\Omega_{i,j,k}^{\psi}(t) \setminus \Gamma(t)) \cap \{\mathbf{x}|x = \xi_i\}$; the cut $y$-normal face, $\Sigma_{i,j,k}^{\psi,y}(t) = (\partial\Omega_{i,j,k}^{\psi}(t) \setminus \Gamma(t)) \cap \{\mathbf{x}|y = \eta_j\}$; as well as the cut $z$-normal face, $\Sigma_{i,j,k}^{\psi,z}(t) = (\partial\Omega_{i,j,k}^{\psi}(t) \setminus \Gamma(t)) \cap \{\mathbf{x}|z = \zeta_k\}$; and finally the intersection between the interface and the control volume, $\Sigma_{i,j,k}^{\Gamma}(t) = \Omega_{i,j,k} \cap \Gamma(t)$. Any of these surfaces might be empty. In addition, we introduce notations for the average convective or diffusive flux over any of these surfaces:

$$ \mathcal{C}_{i,j,k}^{\psi,d}(t) = \frac{1}{S_{i,j,k}^{\psi,d}(t)} \int_{\Sigma_{i,j,k}^{\psi,d}(t)} \mathcal{C}^{\psi} dS, \quad (3) $$

$$ \mathcal{D}_{i,j,k}^{\psi,d}(t) = \frac{1}{S_{i,j,k}^{\psi,d}(t)} \int_{\Sigma_{i,j,k}^{\psi,d}(t)} \mathcal{D}^{\psi} dS, \quad (4) $$

$$ \mathcal{D}_{i,j,k}^{\psi,\Gamma}(t) = \frac{1}{S_{i,j,k}^{\Gamma}(t)} \int_{\Sigma_{i,j,k}^{\Gamma}(t)} \mathcal{D}^{\psi} \cdot ndS, \quad (5) $$

where $d \in \{x, y, z\}$ is an arbitrary coordinate direction, and where $S_{i,j,k}^{\psi,d}(t)$ and $S_{i,j,k}^{\Gamma}(t)$ are respectively the area of $\Sigma_{i,j,k}^{\psi,d}(t)$ and $\Sigma_{i,j,k}^{\Gamma}(t)$. With these notations, Eq. (2) may be expressed as

$$ \frac{d}{dt} \int_{\Omega_{i,j,k}^{\psi}(t)} \rho c_p(\mathbf{x}, t) T(\mathbf{x}, t) dV = \left[ S_{i,j,k}^{\psi,x}(t) (\mathcal{D}_{i,j,k}^{\psi,x}(t) - \mathcal{C}_{i,j,k}^{\psi,x}(t)) - S_{i,j,k}^{\psi,y}(t) (\mathcal{D}_{i,j,k}^{\psi,y}(t) - \mathcal{C}_{i,j,k}^{\psi,y}(t)) \right] \cdot e_x + \left[ S_{i,j,k}^{\psi,y}(t) (\mathcal{D}_{i,j,k}^{\psi,y}(t) - \mathcal{C}_{i,j,k}^{\psi,y}(t)) - S_{i,j,k}^{\psi,z}(t) (\mathcal{D}_{i,j,k}^{\psi,z}(t) - \mathcal{C}_{i,j,k}^{\psi,z}(t)) \right] \cdot e_y + \left[ S_{i,j,k}^{\psi,z}(t) (\mathcal{D}_{i,j,k}^{\psi,z}(t) - \mathcal{C}_{i,j,k}^{\psi,z}(t)) - S_{i,j,k}^{\Gamma}(t) (\mathcal{D}_{i,j,k}^{\Gamma}(t) - \mathcal{C}_{i,j,k}^{\Gamma}(t)) \right] \cdot e_z + S_{i,j,k}^{\Gamma}(t) \mathcal{D}_{i,j,k}^{\psi,\Gamma}(t) \quad (6) $$

by expressing the fluxes on the outer boundary $\partial\Omega_{i,j,k}^{\psi}(t)$ in terms of fluxes on the interface $\Gamma(t)$ and some portions of the outer boundary $\partial\Omega_{i,j,k}^{\Gamma}(t)$. For clarity, this may be formally recast as

$$ \frac{d}{dt} \int_{\Omega_{i,j,k}^{\psi}(t)} \rho c_p(\mathbf{x}, t) T(\mathbf{x}, t) dV = S_{i,j,k}^{\Gamma}(t) \mathcal{D}_{i,j,k}^{\psi,\Gamma}(t) + \sum_{F \in \mathcal{F}_{i,j,k}} S_F(t) (\mathcal{D}_F(t) - \mathcal{C}_F(t)) \cdot n_F, \quad (7) $$

where $\mathcal{F}_{i,j,k}^{\psi}$ is the set of the faces of $\Omega_{i,j,k}^{\psi}$, i.e. excluding the interface.

Let us consider a time step $\bar{\mathcal{T}} = [t_n, t_{n+1}]$, with $t_{n+1} - t_n = \Delta t$. During the time step, the interface is displaced from the location $\Gamma(t_n)$ to $\Gamma(t_{n+1})$, and the geometrical properties of the cut cell $\Omega_{i,j,k}^{\psi}(t)$ evolve accordingly. An equation for the temperature advancement is obtained by integrating Eq. (7) over the time step:

$$ H_{i,j,k}^{\psi}(t_{n+1}) = H_{i,j,k}^{\psi}(t_n) + \int_{t_n}^{t_{n+1}} S_{i,j,k}^{\Gamma}(t) \mathcal{D}_{i,j,k}^{\psi,\Gamma}(t) dt + \int_{t_n}^{t_{n+1}} \sum_{F \in \mathcal{F}_{i,j,k}} S_F(t) (\mathcal{D}_F(t) - \mathcal{C}_F(t)) \cdot n_F dt, \quad (8) $$

where the thermal energy is defined as $H_{i,j,k}^{\psi}(t) = \int_{\Omega_{i,j,k}^{\psi}(t)} \rho c_p T dV$. To simplify this equation, we introduce the notation $(\cdot)$ for time-averaged variables. Namely,

$$ \bar{\mathcal{C}}_{i,j,k}^{\psi,d}(t_n) = \frac{1}{\bar{S}_{i,j,k}^{\psi,d}(t_n) \Delta t} \int_{t_n}^{t_{n+1}} \mathcal{C}_{i,j,k}^{\psi,d}(t) S_{i,j,k}^{\psi,d}(t) dt, \quad (9) $$

$$ \bar{\mathcal{D}}_{i,j,k}^{\psi,d}(t_n) = \frac{1}{\bar{S}_{i,j,k}^{\psi,d}(t_n) \Delta t} \int_{t_n}^{t_{n+1}} \mathcal{D}_{i,j,k}^{\psi,d}(t) S_{i,j,k}^{\psi,d}(t) dt, \quad (10) $$

$$ \bar{\mathcal{D}}_{i,j,k}^{\psi,\Gamma}(t_n) = \frac{1}{\bar{S}_{i,j,k}^{\Gamma}(t_n) \Delta t} \int_{t_n}^{t_{n+1}} \mathcal{D}_{i,j,k}^{\psi,\Gamma}(t) S_{i,j,k}^{\Gamma}(t) dt, \quad (11) $$

where $\bar{S}_{i,j,k}^{\psi,d} = (1/\Delta t) \int_{t_n}^{t_{n+1}} S_{i,j,k}^{\psi,d}(t) dt$ is the effective area of $\Sigma_{i,j,k}^{\psi,d}$ over the interval of time $\bar{\mathcal{T}}$ and $\bar{S}_{i,j,k}^{\Gamma} = (1/\Delta t) \int_{t_n}^{t_{n+1}} S_{i,j,k}^{\Gamma}(t) dt$ is the effective area of $\Sigma_{i,j,k}^{\Gamma}$. With these notations, Eq. (8) may be expressed as

$$ H_{i,j,k}^{\psi}(t_{n+1}) = H_{i,j,k}^{\psi}(t_n) + \bar{S}_{i,j,k}^{\Gamma}(t_n) \bar{\mathcal{D}}_{i,j,k}^{\psi,\Gamma}(t_n) \Delta t + \sum_{F \in \mathcal{F}_{i,j,k}} \bar{S}_F(t_n) (\bar{\mathcal{D}}_F(t_n) - \bar{\mathcal{C}}_F(t_n)) \cdot n_F \Delta t. \quad (12) $$

The notations leading to this equation are summarised in Table 1.

## 2.2. Resolution algorithm and numerical schemes

During a time step of the simulation, the Lagrangian markers of the interface are first displaced in order to calculate the geometrical properties of cut cells before and after the displacement. Then, geometrical properties of cut cells (volumes, surface areas, barycentres) are
<table><caption>**Table 1**
Notations used for the cut-cell discretisation.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Cartesian cell Ω<sub>i,j,k</sub></th><th colspan="2">Cut cell Ω<sup>ψ</sup><sub>i,j,k</sub>(t)</th><th rowspan="2">Interface-Γ</th></tr><tr><th>Cell</th><th>Face-d</th><th>Cell</th><th>Face-d</th></tr></thead><tbody><tr><td>Geometrical object ... instantaneous</td><td>Ω<sub>i,j,k</sub></td><td>Σ<sup>d</sup><sub>i,j,k</sub></td><td>Ω<sup>ψ</sup><sub>i,j,k</sub></td><td>Σ<sup>ψ,d</sup><sub>i,j,k</sub></td><td>Σ<sup>Γ</sup><sub>i,j,k</sub></td></tr><tr><td>Size (volume/area) ... instantaneous</td><td>V<sub>i,j,k</sub></td><td>S<sup>d</sup><sub>i,j,k</sub></td><td>V<sup>ψ</sup><sub>i,j,k</sub></td><td>S<sup>ψ,d</sup><sub>i,j,k</sub></td><td>S<sup>Γ</sup><sub>i,j,k</sub></td></tr><tr><td>...... effective</td><td></td><td></td><td>V<sup>ψ</sup><sub>i,j,k</sub></td><td>S<sup>ψ,d</sup><sub>i,j,k</sub></td><td>S<sup>Γ</sup><sub>i,j,k</sub></td></tr><tr><td>Convective flux ... instantaneous</td><td></td><td></td><td></td><td>Σ<sup>ψ,d</sup><sub>i,j,k</sub></td><td>= 0</td></tr><tr><td>...... effective</td><td></td><td></td><td></td><td>Σ<sup>ψ,d</sup><sub>i,j,k</sub></td><td>= 0</td></tr><tr><td>Diffusive flux ... instantaneous</td><td></td><td></td><td></td><td>Σ<sup>ψ,d</sup><sub>i,j,k</sub></td><td>Σ<sup>ψ,Γ</sup><sub>i,j,k</sub></td></tr><tr><td>...... effective</td><td></td><td></td><td></td><td>Σ<sup>ψ,d</sup><sub>i,j,k</sub></td><td>Σ<sup>Γ</sup><sub>i,j,k</sub></td></tr><tr><td>Thermal energy ... instantaneous</td><td></td><td></td><td>H<sup>ψ</sup><sub>i,j,k</sub></td><td></td><td></td></tr><tr><td>Temperature ... instantaneous</td><td></td><td></td><td>T<sup>ψ</sup><sub>i,j,k</sub></td><td></td><td>T<sup>ψ,d</sup><sub>i,j,k</sub></td></tr></tbody></table>

calculated. Finally, the thermal energy in each mesh cell is advanced using Eq. (12). To close the exact Eq. (12), a model is provided for the effective surfaces and the effective convective fluxes. This article considers several time-integration methods for the effective surfaces (Section 2.2.4) in the context of an explicit approximation of the diffusive and convective fluxes, that is $\mathcal{C}^{\psi,x}_{i,j,k}(t_n) = \mathcal{C}^{\psi,x}_{i,j,k}(t_n)$. With such an explicit approximation of fluxes, numerical stability requires an unreasonable restriction of the time step $\Delta t$, based on the size of the smallest non-empty cut cell. This is avoided by correcting the predicted temperature in small cells using a stable method, which can be classified as a flux-redistribution method [67–69].

### 2.2.1. Overview of the algorithm

The resolution algorithm is subdivided in four main steps, which correspond respectively to the main and corrective updates of convection and diffusion (Fig. 2), where the main update is related to the explicit approximation of fluxes in Eq. (12) and the corrective update is related to the special handling of small cells for stability. Namely,

$$H_{i,j,k}^{\psi}(t^{(1)}) = H_{i,j,k}^{\psi}(t_n) - \sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \bar{S}_F(t_n) \mathcal{C}_F(t_n) \cdot n_F \Delta t, \quad (13)$$

$$H_{i,j,k}^{\psi}(t^{(2)}) = H_{i,j,k}^{\psi}(t^{(1)}) + \sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \bar{S}_F(t_n) \mathcal{C}_F^{\psi}(t^{(1)}) \cdot n_F \Delta t, \quad (14)$$

$$H_{i,j,k}^{\psi}(t^{(3)}) = H_{i,j,k}^{\psi}(t^{(2)}) + \bar{S}_{i,j,k}^{\Gamma}(t_n) \bar{\mathcal{D}}_{i,j,k}^{\psi,\Gamma}(t^{(2)}) \Delta t + \sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \bar{S}_F(t_n) \bar{\mathcal{D}}_F(t^{(2)}) \cdot n_F \Delta t, \quad (15)$$

$$H_{i,j,k}^{\psi}(t_{n+1}) = H_{i,j,k}^{\psi}(t^{(3)}) + \sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \bar{S}_F(t_n) \bar{\mathcal{D}}_F^{\psi}(t^{(3)}) \cdot n_F \Delta t, \quad (16)$$

where $(t^{(1)})$, $(t^{(2)})$ and $(t^{(3)})$ formally denote the intermediate times after the main convective update, the corrective convective update and the main diffusive update, respectively. The superscript c denotes corrective fluxes. Corrective fluxes are only applied, for each time step, on the cells that are either

* dying: $V_{i,j,k}^{\psi}(t_{n+1}) = 0$ whereas $V_{i,j,k}^{\psi}(t_n) \neq 0$;

* small: $0 < V_{i,j,k}^{\psi}(t_{n+1}) < \beta_{th}V_{i,j,k}$, while not nascent; or

* nascent: $V_{i,j,k}^{\psi}(t_{n+1}) \neq 0$ whereas $V_{i,j,k}^{\psi}(t_n) = 0$.

The threshold value for small cells, $\beta_{th}$ should be sufficiently large to ensure a stable simulation, but not excessively large as this generally decreases the accuracy of the simulation. It is set to $\beta_{th} = 2.5\%$ of the volume of a Cartesian cell in the present study. Diffusion is considered after convection in the four-step algorithm (13)–(16) to facilitate taking into account the diffusive fluxes in both dying cells, that vanish during the time step, and nascent cells, that appear during the time step.

### 2.2.2. Main schemes

The main schemes (steps (1) and (3) in Fig. 2) advance the temperature in time by providing a discretised expression for the fluxes

$\mathcal{C}^{\psi,x}_{i,j,k}(t_n)$, $\bar{\mathcal{D}}_{i,j,k}^{\psi,x}(t_n)$ and $\bar{\mathcal{D}}_{i,j,k}^{\psi,\Gamma}(t_n)$ in Eq. (12). First, the convective flux on a given face $F \in \bar{\mathcal{F}}_{i,j,k}$ of the cut cell $\Omega_{i,j,k}^{\psi}(t)$ may be expressed as

$$\mathcal{C}_F(t_n) = \hat{T}_F(t_n) u_F(t_n), \quad (17)$$

where $\hat{T}_F(t_n)$ represents the temperature passing through the surface and $u_F(t_n)$ is the velocity over the surface $\Sigma_F$. As stated above, explicit values, at time $t_n$, are used to approximate the effective time-averaged quantities. Given the staggered Cartesian grid approach, we suppose that the velocity $u_F(t_n)$ corresponds to the discretised velocity field associated with the face location. The temperature $\hat{T}_F(t_n)$ is estimated using a QUICK scheme [112], provided that the stencil is available, and an upwind scheme otherwise. For that purpose, a cell is only considered available to the QUICK stencil if its barycentre does not deviate by more than one-tenth grid spacing from the Cartesian cell centre (Fig. 3).

Second, the diffusive flux is calculated using a normal-probe technique on the interface and a finite-difference scheme on cut Cartesian faces. The finite-difference scheme used on Cartesian faces to approximate the temperature gradient is given by,

$$\mathcal{D}_{i,j,k}^{\psi,x}(t_n) = \lambda_{\psi} \frac{T_{i,j,k}^{\psi}(t_n) - T_{i-1,j,k}^{\psi}(t_n)}{x_{i,j,k}^{\psi}(t_n) - x_{i-1,j,k}^{\psi}(t_n)} e_x, \quad (18)$$

$$\mathcal{D}_{i,j,k}^{\psi,y}(t_n) = \lambda_{\psi} \frac{T_{i,j,k}^{\psi}(t_n) - T_{i,j-1,k}^{\psi}(t_n)}{y_{i,j,k}^{\psi}(t_n) - y_{i,j-1,k}^{\psi}(t_n)} e_y, \quad (19)$$

$$\mathcal{D}_{i,j,k}^{\psi,z}(t_n) = \lambda_{\psi} \frac{T_{i,j,k}^{\psi}(t_n) - T_{i,j,k-1}^{\psi}(t_n)}{z_{i,j,k}^{\psi}(t_n) - z_{i,j,k-1}^{\psi}(t_n)} e_z. \quad (20)$$

Note that these equations neglect the fact that cell barycentres are not generally aligned with the Cartesian axes in cut cells. This non-orthogonality issue has been discussed in [113]. Presently, Eqs. (18)–(20) can nevertheless be used for simplicity because the inaccuracy does not concern the flux on the interface and diffusion is well resolved in our test case (Section 3.1). The interfacial flux $\mathcal{D}_{i,j,k}^{\psi,\Gamma}(t_n)$, is expressed from the contributions of the triangular facets of the interface mesh that intersects the cell. Let use denote $x_{\text{fct}}^{(\ell)}(t_n)$ the position vector of the centre of the $\ell$th facet, $T_{\text{fct}}^{(\ell)}(t_n)$ its temperature and $\mathcal{D}_{\text{fct}}^{(\ell)}(t_n)$ its diffusive flux. The cell diffusive flux is

$$\mathcal{D}_{i,j,k}^{\psi,\Gamma}(t_n) = \frac{\sum_{\ell=1}^{N_{\text{fct}}} S_{i,j,k}^{\text{int-}\ell}(t_n) \mathcal{D}_{\text{fct}}^{(\ell)}(t_n)}{\sum_{\ell=1}^{N_{\text{fct}}} S_{i,j,k}^{\text{int-}\ell}(t_n)}, \quad (21)$$

where $S_{i,j,k}^{\text{int-}\ell}(t_n)$ is the surface of the intersection between the $\ell$th facet and the cell $\Omega_{i,j,k}$. The calculation of the facet flux $\mathcal{D}_{\text{fct}}^{(\ell)}(t_n)$ is based on the interpolation $T_1(t_n)$ and $T_2(t_n)$ of temperature on two points at a normal distance $d_1$ and $d_2$ of the facet, towards phase $\psi_1$ and $\psi_2$ respectively (Fig. 4). A tetrahedral interpolation is chosen for simplicity, as the tetrahedron is the simplest three-dimensional shape that allows for three-dimensional interpolation. Since the flux computed from the phases $\psi_1$ and $\psi_2$ are equal, we obtain with a first-order approximation

$$\mathcal{D}_{\text{fct}}^{(\ell)}(t_n) = \lambda_1 \frac{T_1(t_n) - T_{\text{fct}}^{(\ell)}(t_n)}{d_1} = -\lambda_2 \frac{T_2(t_n) - T_{\text{fct}}^{(\ell)}(t_n)}{d_2}, \quad (22)$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Four main steps of the resolution algorithm. The solid, dashed and dotted red lines represent respectively the current, next and previous positions of the interface. Double-headed arrows represent fluxes.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Stencil of the main convective scheme on a cut cell, depending on the interface location, assuming a velocity vector pointing towards the right direction. The bold vertical line represents the target flux location. Upper two rows: an upwind scheme is used. Lower two rows: a QUICK scheme is used [112]. The second and third rows differ only in the extent of the discrepancy between the cut-cell barycentre and the Cartesian cell barycentre.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. To calculate the conductive flux on a facet of the interface, the temperature field is interpolated linearly on two points at a given normal distance of the facet, chosen here to be to the Cartesian grid spacing, by finding tetrahedra (green dashed-dotted lines) enclosing the points.</figcaption></figure>

which corresponds to the following interfacial temperature on facet $\ell$:

$$T_{\text{fct}}^{(\ell)}(t_n) = \frac{\lambda_1 T_1(t_n)/d_1 + \lambda_2 T_2/d_2}{\lambda_1/d_1 + \lambda_2/d_2}. \qquad (23)$$

Diffusive fluxes on cut faces are calculated using an explicit approximation, similarly to convective fluxes. The interfacial diffusive flux is calculated using both the initial temperature field $T(t_n)$ and the convected temperature field $T(t^{(2)})$. Specifically, we use a robust algebraic

approximation of the time average:

$$\bar{\mathcal{D}}_{i,j,k}^{\psi, \Gamma}(t_n) = \begin{cases} (1/2)(\mathcal{D}_{i,j,k}^{\psi, \Gamma}(t_n) + \mathcal{D}_{i,j,k}^{\psi, \Gamma}(t^{(2)})) & \text{if } S_{i,j,k}^{\Gamma}(t_n) \neq 0 \text{ and } S_{i,j,k}^{\Gamma}(t_{n+1}) \neq 0, \\ \mathcal{D}_{i,j,k}^{\psi, \Gamma}(t^{(2)}) & \text{if } S_{i,j,k}^{\Gamma}(t_n) = 0 \text{ and } S_{i,j,k}^{\Gamma}(t_{n+1}) \neq 0, \\ \mathcal{D}_{i,j,k}^{\psi, \Gamma}(t_n) & \text{if } S_{i,j,k}^{\Gamma}(t_n) \neq 0 \text{ and } S_{i,j,k}^{\Gamma}(t_{n+1}) = 0, \\ 0 & \text{if } S_{i,j,k}^{\Gamma}(t_n) = 0 \text{ and } S_{i,j,k}^{\Gamma}(t_{n+1}) = 0, \end{cases} \qquad (24)$$

This allows the method to take into account the interfacial diffusive flux in both dying and nascent cells.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Schematic representation of the convective (left) and diffusive (right) corrections.</figcaption></figure>

### 2.2.3. Corrective schemes

The corrective schemes presuppose a time advancement of temperature, then ensure conservation by providing fluxes consistent with Eqs. (14) and (16). That is to say, assuming that a method to determine the final temperature of the step $T_{i,j,k}^{\psi}(t^{(2)})$ or $T_{i,j,k}^{\psi}(t_{n+1})$ is available, as discussed further below, the fluxes are calculated based on the desired temperature variation:

$$ \mathcal{C}_F^{\mathcal{C}}(t^{(1)}) = \pi_F^{\mathcal{C}}(t^{(1)}) \left( \rho c_{p,\psi} T_{i,j,k}^{\psi}(t^{(2)}) V_{i,j,k}^{\psi}(t_{n+1}) - H_{i,j,k}^{\psi}(t^{(1)}) \right), \quad (25) $$

$$ \mathcal{D}_F^{\mathcal{D}}(t^{(3)}) = \pi_F^{\mathcal{D}}(t^{(3)}) \left( \rho c_{p,\psi} T_{i,j,k}^{\psi}(t_{n+1}) V_{i,j,k}^{\psi}(t_{n+1}) - H_{i,j,k}^{\psi}(t^{(3)}) \right). \quad (26) $$

The vectorial weights $\pi_F^{\mathcal{C}}(t^{(1)})$ and $\pi_F^{\mathcal{D}}(t^{(3)})$ introduced in these equations are coefficients per unit surface and time. Indeed, to ensure that Eqs. (14) and (16) are satisfied, the vectorial weights must obey the following summation properties:

$$ \sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \bar{S}_F(t_n) \pi_F^{\mathcal{C}}(t^{(1)}) \cdot n_F \Delta t = 1, \quad (27) $$

$$ \sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \bar{S}_F(t_n) \pi_F^{\mathcal{D}}(t^{(3)}) \cdot n_F \Delta t = 1. \quad (28) $$

In that frame, a corrective scheme is determined by providing an expression for the final temperature of the step, and the vectorial weights $\pi_F^{\mathcal{C}}(t^{(1)})$ and $\pi_F^{\mathcal{D}}(t^{(3)})$, which indicate the relative magnitude of the various fluxes in a given cell.

In the case of convection (Fig. 5(a)), the corrective fluxes are assumed to be proportional to the discretised convective fluxes (17), that is

$$ \pi_F^{\mathcal{C}}(t^{(1)}) = \frac{\epsilon_F(t_n) \hat{T}_F(t_n) \|\mathbf{u}_F(t_n) \cdot \mathbf{n}_F \| \mathbf{n}_F}{\sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \epsilon_F(t_n) \bar{S}_F(t_n) \hat{T}_F(t_n) \|\mathbf{u}_F(t_n) \cdot \mathbf{n}_F \|}, \quad (29) $$

where $\hat{T}_F(t_n)$ is evaluated using the upwind temperature on the faces of cells, and the neighbouring temperature in small or nascent cells, avoiding the use of the non-initialised temperature of the nascent phase. The boolean variable $\epsilon_F(t_n)$ selects the fluxes that should be used in the computation. To avoid coupling the equations of different dying, small or nascent cells, it is set to zero between a dying cell and a nascent cell, as well as between a small cell and another smaller cell. This essentially means that the corrective variation of temperature is only redistributed in fluxes towards larger cells. In this paper, the corrected temperature $T_{i,j,k}^{\psi}(t^{(2)})$, in small or nascent cells, corresponds to the average of the temperature of the neighbouring cells, weighted by the discretised convective fluxes. Note that it is also possible to determine $T_{i,j,k}^{\psi}(t^{(2)})$ using a semi-Lagrangian approach [114–118] which leads to overall similar results in our numerical tests.

In the case of diffusion (Fig. 5(b)), it is sometimes possible to perform a calculation without any particular correction of diffusive schemes in small and nascent cells, provided that the convective stability time step is limiting. In practical applications, the convective stability time step is for example limiting in direct numerical simulations of bubble swarm [35,119]. If the diffusive stability time step is of the same order or larger than the convective stability time step, a

correction is however necessary. In this paper, the diffusive correction sets the temperature increment to zero, i.e. $H_{i,j,k}^{\psi}(t^{(2)}) = H_{i,j,k}^{\psi}(t_n)$. This means that there is no increase or decrease of temperature within small or nascent cells due to diffusion. This classical assumption is similar to the hypothesis of Grosso et al. [41] in a ghost-fluid method. The goal of the diffusive correction is thus to transfer the interfacial diffusive flux to neighbouring cells in a conservative manner. Since the interfacial flux is normal to the interface, it is natural to set, on a given face $F \in \mathcal{F}_{i,j,k}^{\psi}$ of a cut cell $\Omega_{i,j,k}^{\psi}(t_n)$,

$$ \pi_F^{\mathcal{D}}(t^{(3)}) = \frac{\epsilon_F(t_n) \bar{n}_{i,j,k}^{\psi,\Gamma}(t_n)}{\sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \epsilon_F(t_n) \bar{S}_F(t_n) \bar{n}_{i,j,k}^{\psi,\Gamma}(t_n) \cdot \mathbf{n}_F}, \quad (30) $$

where

$$ \bar{n}_{i,j,k}^{\psi,\Gamma}(t_n) = \begin{cases} (1/2)(n_{i,j,k}^{\psi,\Gamma}(t_n) + n_{i,j,k}^{\psi,\Gamma}(t_{n+1})) & \text{if } V_{i,j,k}^{\psi}(t_n) \neq 0 \text{ and } V_{i,j,k}^{\psi}(t_{n+1}) \neq 0, \\ n_{i,j,k}^{\psi,\Gamma}(t_{n+1}) & \text{if } V_{i,j,k}^{\psi}(t_n) = 0 \text{ and } V_{i,j,k}^{\psi}(t_{n+1}) \neq 0, \\ n_{i,j,k}^{\psi,\Gamma}(t_n) & \text{if } V_{i,j,k}^{\psi}(t_n) \neq 0 \text{ and } V_{i,j,k}^{\psi}(t_{n+1}) = 0, \\ 0 & \text{if } V_{i,j,k}^{\psi}(t_n) = 0 \text{ and } V_{i,j,k}^{\psi}(t_{n+1}) = 0, \end{cases} \quad (31) $$

is the effective normal to the interface. The diffusive correction is for consistency used in all simulations in the present paper.

### 2.2.4. Time integration of cut surfaces

A model is required to evaluate the effective surfaces involved in the calculation of the fluxes, since the exact temporal evolution of the cut surfaces is unknown. An important property for the set of effective surfaces used in the numerical method is the agreement between the variation of cell volume and the volumetric fluxes passing through the cut Cartesian faces, hereafter referred to as the dynamic conservation of volume. The equation for the dynamic conservation of volume may be obtained by following the derivations of Section 2.1 without temperature field ($T(x, t) = 1$), using an explicit approximation of the face velocity, $\bar{u}_F(t_n) = u_F(t_n)$. This leads to

$$ V_{i,j,k}^{\psi}(t_{n+1}) = V_{i,j,k}^{\psi}(t_n) + \sum_{F \in \mathcal{F}_{i,j,k}^{\psi}} \bar{S}_F(t_n) u_F(t_n) \cdot n_F \Delta t. \quad (32) $$

This equation is also referred to as the geometric conservation law in the literature, notably in the context of arbitrary Lagrangian–Eulerian methods [96–98]. A discrepancy on the dynamic conservation of volume (32) results in an erroneous variation of temperature within a cell, as convective fluxes do not match the energy variation due to the change in cell volume. This implies in particular that a uniform temperature field is not preserved during a convective update. To avoid any confusion, let us note however that the total energy of each phase is conserved by the convective update regardless of the time-integration method or whether the dynamic volume conservation (32) is satisfied.

Three models for the effective surface of cut Cartesian faces are evaluated in the present paper. With the explicit time-integration method, cut-face areas are evaluated using the location of the interface at time
$t_n$. This provides an explicit approximation of cut-face areas,

$$ \bar{S}_F^{\text{expl}}(t_n) = \begin{cases} S_F(t_n) & \text{if } S_F(t_n) \neq 0, \\ S_F(t_{n+1}) & \text{if } S_F(t_n) = 0 \text{ and } S_F(t_{n+1}) \neq 0, \\ 0 & \text{if } S_F(t_n) = 0 \text{ and } S_F(t_{n+1}) = 0, \end{cases} \quad (33) $$

An exception is given for nascent surfaces. In that case, the surface at time $t_{n+1}$ is used to avoid a zero cut-face area.

With the algebraic semi-implicit time-integration method, the calculation of cut-face areas takes into account both the initial and final displaced location of the interface. Specifically, the semi-implicit approximation corresponds to the arithmetic average of the initial and final areas of the current time step,

$$ \bar{S}_F^{\text{algb}}(t_n) = \begin{cases} (1/2)(S_F(t_n) + S_F(t_{n+1})) & \text{if } 0 < S_F(t_n) < 1 \text{ and} \\ & \quad 0 < S_F(t_{n+1}) < 1, \\ (1/4)(3S_F(t_n) + S_F(t_{n+1})) & \text{if } S_F(t_n) = 0 \lor S_F(t_n) = 1 \text{ and} \\ & \quad 0 < S_F(t_{n+1}) < 1, \\ (1/4)(S_F(t_n) + 3S_F(t_{n+1})) & \text{if } 0 < S_F(t_n) < 1 \text{ and} \\ & \quad S_F(t_{n+1}) = 0 \lor S_F(t_{n+1}) = 1, \\ S_F(t_n) & \text{otherwise} \end{cases} \quad (34) $$

The calculation is intentionally biased towards the value of the uncut surface in the case of a dying or nascent surface. Indeed, it is clear that, in that case, the face area does not evolve linearly during the entire time step. The present approximation is obtained by assuming that the death (or equivalently birth) of the surface occurs at the time $(t_n+t_{n+1})/2$. and that the surface evolves linearly during a half time step.

With the volume-conserving semi-implicit time-integration method, a correction to the algebraic semi-implicit surface (Eq. (34)) is sought, in order to ensure the dynamic conservation of volume. This may be achieved e.g. by seeking the minimum-norm solution for the increment of surface $\delta \bar{S}_F^{\text{cons}}(t_n)$ in the linear system (35),

$$ V_{i,j,k}^{\psi}(t_{n+1}) = V_{i,j,k}^{\psi}(t_n) + \sum_{F \in \bar{S}_F^{\psi}} (\bar{S}_F^{\text{algb}}(t_n) + \delta \bar{S}_F^{\text{cons}}(t_n)) u_F(t_n) \cdot n_F \Delta t. \quad (35) $$

Assuming that the cut surface $S_F(t)$ evolves smoothly with time, the explicit time integration provides a first-order accurate approximation of the mean surface $\bar{S}_F(t_n)$, while the corresponding volume variation $\Delta V = V_{i,j,k}^{\psi}(t_{n+1}) - V_{i,j,k}^{\psi}(t_n)$ is approximated at the second order in time. On the other hand, the algebraic semi-implicit time integration provides a second-order accurate approximation of the mean surface $\bar{S}_F(t_n)$. In practice, the function $S_F(t)$ is however not continuously differentiable. Indeed, the slope of $S_F(t)$ changes abruptly when the interface-face intersection crosses a vertex of the Cartesian grid or of the interface, since this modifies the shape of the intersection. This includes all cases where crossing a vertex vanishes (or brings forth) a surface, i.e. dying (or nascent) surfaces. A continuously differentiable evolution of $S_F(t)$ is thus only possible without dying or nascent surfaces and with a smooth (non-discretised) interface.

As opposed to the time-integration of cut Cartesian faces, the effect of the time-integration method of the interfacial effective area is generally negligible. In the present paper, an algebraic semi-implicit time integration (Eq. (34)) is used to compute the effective area $\bar{S}_{i,j,k}^R(t_n)$ of the interface in all simulations. This allows the analysis to focus on the effect of time integration on cut Cartesian faces, which is directly related to the dynamic conservation of volume and hence the accurate resolution of scalar transport.

# 3. Numerical experiment

## 3.1. Case definition

This section investigates the heat diffusion from an initially hot spherical particle into the surrounding fluid, as seen from a moving observer. This test case is specifically designed to assess the relevance of the numerical method. Indeed, the apparent motion of the particle modifies the calculation of convective and diffusive heat fluxes, numerically, without changing the underlying physics. Specifically, the apparent motion of the particle as seen by the observer follows a cyclic ellipsoidal motion. The velocity vector $u(t) = (u, v, w)(t)$ is given by

$$ \begin{align} u(t) &= u_0 \cos(2\pi t/\tau), \notag \\ v(t) &= 0, \tag{36} \\ w(t) &= (u_0/2) \sin(2\pi t/\tau). \notag \end{align} $$

where $u_0$ is a constant velocity, $t$ is time and $\tau = 0.5$ s is the period of oscillation. The resulting trajectory is an ellipse with semimajor axis length $a_0 = u_0\tau/(2\pi)$ and eccentricity $e_0 = \sqrt{3}/2$.

Two sets of phase properties (P1 and P2) are considered, as reported in Table 2. The two cases are selected to provide a distinctive physical behaviour. In case P1, the particle has an artificially high density and thermal conductivity; the temperature within the particle is quickly homogenised. In case P2, the particle has a low density and a moderate thermal conductivity; which is representative of a bubble of pressurised air in a coolant liquid. In both cases, the diameter of the spherical particle is $d_p = 2$ mm and its radius denoted $R = d_p/2$. At the initial time $t = 0$, the temperature field is set either, in case P2, to a Gaussian decay $T_0 + \Delta T_0 \exp(-(\|x\|/d_p)^2/2)$, with $T_0 = 293$ K the background temperature and $\Delta T_0 = 80$ K the amplitude of the step or, in case P1, to a Heaviside step $T(x, 0) = T_0 + \Delta T_0 [ \|x\| < d_p/2 ]$, where $[ \cdot ]$ are Iverson brackets. In case P1, the difference in initial condition has little effect on the resulting profiles, as the temperature would regardless quickly become uniform inside the particle. For clarity, all results are presented in terms of the nondimensionalised temperature $T^* = (T - T_0)/\Delta T_0$.

## 3.2. Numerical setting

The computational domain is a cube with side length $L = 5d_p$. A periodic boundary condition is used on the domain boundaries, representing an infinite three-dimensional array of particles with spacing $L$. The particle is initially placed at the centre of the domain, conventionally the location $(0, 0, 0)$. The domain is discretised using a Cartesian grid of uniform spacing $\Delta = d_p/m$ in the three directions, where $m$ is the number of cells per particle diameter. Several mesh resolutions are considered, ranging from $m = 5.2$ to $m = 41.6$ (Table 3). The spherical interface between the two phases is represented by a triangular mesh with typical edge length $\Delta_F/d_p = 0.196$ in case P1 and $\Delta_F/d_p = 0.035$ in case P2, regardless of the Cartesian mesh resolution. The geometrical properties of cut cells, such as volumes, surface areas or barycentres, are calculated based on the intersection of the discretised interface and the Cartesian mesh. The calculations are generic, and do not make use of the particular geometry of this test case.

The numerical simulation of the case is addressed using the thermohydraulics solver TRUST/ TrioCFD [120], which has been used to simulate various fluid phenomena [41,119,121,122]. Simulations are performed using the cut-cell formulation presented in Section 2 and a baseline one-field formulation. The baseline one-field method solves the heat Eq. (1) in temperature form using a staggered grid system, a third-order QUICK scheme for convection [112] and a second-order centred scheme for diffusion. It uses a harmonic averaging of heat conductivity and heat capacity in mixed cells, as presented in [123] (TOF formulation). This regularises the discontinuity between the phases and introduces first-order errors in the solution [32]. The cut-cell method
<caption>**Table 2**</caption>

Operating parameters of the four considered cases. Dotted rows (...) in a sub-variant P2-slow or P2-fast indicates that the parameter is not modified compared to P2 case.

<table><caption>**Table 3**</caption><thead><tr><th>Case</th><th>P1</th><th>P2</th><th>P2-slow</th><th>P2-fast</th></tr></thead><tbody><tr><td>Properties in the outer domain:</td><td></td><td></td><td></td><td></td></tr><tr><td>Heat conductivity, λ₁ [W m<sup>-1</sup>K<sup>-1</sup>]</td><td>5.500 × 10<sup>-2</sup></td><td>6.898 × 10<sup>-1</sup></td><td>...</td><td>...</td></tr><tr><td>Heat capacity, ρc<sub>p,1</sub> [J K<sup>-1</sup> m<sup>-3</sup>]</td><td>7.028 × 10<sup>+4</sup></td><td>6.009 × 10<sup>+6</sup></td><td>...</td><td>...</td></tr><tr><td>Peclet number, Pe<sub>1</sub> = ρ<sub>1</sub>c<sub>p,1</sub>u<sub>0</sub>d<sub>p</sub>/λ<sub>1</sub></td><td>102</td><td>697</td><td>139</td><td>3485</td></tr><tr><td>Properties within the particle:</td><td></td><td></td><td></td><td></td></tr><tr><td>Heat conductivity, λ₂ [W m<sup>-1</sup>K<sup>-1</sup>]</td><td>1.550 × 10<sup>+1</sup></td><td>2.210 × 10<sup>-2</sup></td><td>...</td><td>...</td></tr><tr><td>Heat capacity, ρc<sub>p,2</sub> [J K<sup>-1</sup> m<sup>-3</sup>]</td><td>7.028 × 10<sup>+5</sup></td><td>2.526 × 10<sup>+5</sup></td><td>...</td><td>...</td></tr><tr><td>Peclet number, Pe<sub>2</sub> = ρ<sub>2</sub>c<sub>p,2</sub>u<sub>0</sub>d<sub>p</sub>/λ<sub>2</sub></td><td>3.6</td><td>914</td><td>183</td><td>4572</td></tr><tr><td>Number of elliptic cycles, N<sub>0</sub></td><td>5</td><td>5</td><td>1</td><td>25</td></tr><tr><td>Travelled distance, ϵ<sub>0</sub> ≈ 4.84N<sub>0</sub>u<sub>0</sub>τ/(2πd<sub>p</sub>)</td><td>39</td><td>39</td><td>8</td><td>193</td></tr><tr><td>Initial temperature field, T*(x, 0)</td><td>"Heaviside": [||x|| &lt; d<sub>p</sub>/2]</td><td>"Gaussian": exp(-(||x||/d<sub>p</sub>)<sup>2</sup>/2)</td><td>...</td><td>...</td></tr></tbody></table>

List of mesh resolutions considered, with corresponding number of cells per particle diameter *m* and numerical Peclet number Pe<sub>Δ</sub> = (Δ/d<sub>p</sub>) min(Pe<sub>1</sub>, Pe<sub>2</sub>).

<table><thead><tr><th rowspan="2">Total number of cells</th><th rowspan="2">m = d<sub>p</sub>/Δ</th><th colspan="5">Numerical Peclet number, Pe<sub>Δ</sub></th></tr><tr><th>Case P1</th><th>Case P2</th><th>Case P2-slow</th><th>Case P2-fast</th></tr></thead><tbody><tr><td>26 × 26 × 26</td><td>5.2</td><td>0.70</td><td>134</td><td>26.8</td><td>670</td></tr><tr><td>36 × 36 × 36</td><td>7.2</td><td>0.50</td><td>97</td><td>19.4</td><td>484</td></tr><tr><td>52 × 52 × 52</td><td>10.4</td><td>0.35</td><td>67</td><td>13.4</td><td>335</td></tr><tr><td>72 × 72 × 72</td><td>14.4</td><td>0.25</td><td>48</td><td>9.7</td><td>241</td></tr><tr><td>104 × 104 × 104</td><td>20.8</td><td>0.17</td><td>34</td><td>6.7</td><td>168</td></tr><tr><td>144 × 144 × 144</td><td>28.8</td><td>0.13</td><td>24</td><td>4.8</td><td>121</td></tr><tr><td>208 × 208 × 208</td><td>41.6</td><td>0.09</td><td>17</td><td>3.3</td><td>84</td></tr></tbody></table>

uses the schemes detailed in Section 2.2. The necessary geometric information are obtained based on the discretised interface mesh. On each time step, the Lagrangian markers of the interface are displaced according to the velocity (36). This is achieved without any interpolation, since the velocity field is uniform in space. The volume of each cut cell and the area of each cut surface are deduced for each time. This allows then the computation of the effective surfaces using one of the different algorithms presented in Section 2.2.4. All simulations are performed with a forward Euler time scheme. The time step is set according to the stability criterion Δt = F<sub>t</sub> min(Δt<sup>Δ</sup>, Δt<sup>Δ</sup>), where F<sub>t</sub> is a coefficient of time stepping, Δt<sup>Δ</sup> = (1/2)Δ/ max(u, v, w) is the convective stability time step and Δt<sup>Δ</sup> = (1/6)Δ<sup>2</sup>/ max(λ<sub>ψ</sub>/(ρ<sub>ψ</sub>c<sub>p,ψ</sub>)) the diffusive stability time step. From empirical observations, a time-stepping coefficient F<sub>t</sub> = 0.1 has been used in direct numerical simulations of bubble swarm [35, 119], while values ranging from F<sub>t</sub> = 0.02 to 0.8 are compared in this study. The numerical Peclet number Pe<sub>Δ</sub> = (Δ/d<sub>p</sub>) min(Pe<sub>1</sub>, Pe<sub>2</sub>) compares the convective and diffusive stability time steps. With the mesh resolutions considered (Table 3), diffusion is limiting in case P1 (Pe<sub>Δ</sub> < 1) whereas convection is limiting in case P2.

A reference for the simulations is obtained by resolving heat diffusion in the static case (u<sub>0</sub> = 0), without convection. A high-fidelity solution is provided by solving the heat equation in the one-dimensional spherical case, i.e.

$$ \frac{\partial \rho c_p T}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left( \lambda r^2 \frac{\partial T}{\partial r} \right), \qquad (37) $$

where *r* = ||*x*|| is the radial coordinate, on a mesh fitted to the interface boundary *r* = *d*<sub>p</sub>/2. This is performed with a one-dimensional finite difference scheme with a spacing *d*<sub>p</sub>/80.

### 3.3. Numerical results in case P1

The numerical simulations of case P1 with the baseline and cut-cell methods are compared in Fig. 6. Within the range of mesh resolutions studied, the cut-cell method leads to a more accurate temperature profile in the static pure-diffusion case than the baseline one-field

method for a given mesh resolution (Fig. 6(a)). Furthermore, when convective effects are taken into account, the baseline one-field formulation overestimates the effective heat transfers between the two phases, leading to a dramatically inaccurate temperature profile (Fig. 6(b)). This is foremost due to the numerical diffusion between the two phases induced by the convective scheme, as the temperature profile is more satisfactory without convection. With a cut-cell method, numerical diffusion cannot transfer heat from one phase to another. As a result, the cut-cell method is able to predict a reasonable temperature profile with convection as well as without convection.

Fig. 7(a) compares several time-integration methods for cut surfaces in the cut-cell method using F<sub>t</sub> = 0.2. The difference between the different time-integration methods is clearly negligible. This can be explained by the fact that diffusion is limiting in case P1. Indeed, the time step is dominated by the low Peclet number of the particle and the numerical Peclet number of the simulation is firmly below 1 (Table 3). As a result, convective transfers are finely resolved in time. This observation holds true even for large values of the time-stepping coefficient F<sub>t</sub>, up to the stability limit (Fig. 7(b)). This case is useful to show that an explicit time-integration method can be used effectively to compute cut surfaces in two-phase flow simulations, provided that the simulation time step is limited by diffusion.

### 3.4. Numerical results in case P2

Convective transfers are more important in case P2 because the Peclet number is larger, both in the outer domain and within the particle. The temperature profile within the particle is not uniformised, as the centre retains a higher temperature. In the static pure-diffusion case, the reference temperature profile is well reproduced with both the baseline one-field method and the cut-cell method (Fig. 8(a)). For a given mesh resolution, the temperature profile is closer to the reference with the cut-cell method than with the baseline method. In the moving-particle setup (Fig. 8(b)), the mesh convergence of the temperature profile is obtained with the cut-cell method, using F<sub>t</sub> = 0.05. It requires a finer mesh resolution, with *m* = 41.6 cells per particle diameter, to prevent the cumulative errors due to the convective schemes to affect the solution. In the baseline one-field formulation, the solution does not agree with the reference and the temperature-gradient discontinuity at the interface is not reproduced.

The temporal convergence of the profile predicted by the cut-cell method is examined in Fig. 9(a) with the finer mesh resolution *m* = 41.6. The results are converged in time using the selected time-stepping coefficient F<sub>t</sub> = 0.05 with algebraic semi-implicit time integration of the cut surfaces. With an explicit time-integration method, the temperature profile is not yet converged in time and deviates from the reference profile. We may also observe a small discrepancy between the algebraic and volume-conserving semi-implicit time-integration methods very near to the interface, which is overall negligible. With a larger
<figure><img src="image_9.png" /><figcaption>Fig. 6. Temperature profile along the x-axis for increasing mesh resolutions, in case P1, using $F_t = 0.2$, with the baseline one-field method and the cut-cell method with algebraic semi-implicit time integration. Aberrant values of the temperature with the cut-cell method occur near the interface without particle motion because the temperature of small cells is not updated during the simulation, which does not affect the validity of the simulations. The solid line (—) is the reference static one-dimensional spherical solution.</figcaption></figure>

Baseline: —○— $m = 5.2$, ·□· $m = 7.2$, ·△· $m = 10.4$, ·∇· $m = 14.4$, ·◇· $m = 20.8$;

Cut-cell: —●— $m = 5.2$, ·■· $m = 7.2$, ▲· $m = 10.4$, —▼— $m = 14.4$, ·♦· $m = 20.8$.

<figure><img src="image_10.png" /><figcaption>Fig. 7. Effect of the cut-surface time-integration method of the cut cell method on the temperature profile along the x-axis, in case P1, using a mesh resolution $m = 20.8$.</figcaption></figure>

— Reference; — Explicit; — ∇— Algebraic semi-implicit; — ◊— Volume-conserving semi-implicit.

time-stepping coefficient (Figs. 9(b)–(d)), the profile obtained with explicit time integration deviates further from the reference profile. The temperature variations within the particle are smoothened, producing a large temperature discontinuity at the interface and an underestimation of the interfacial diffusive flux. With an algebraic or volume-conserving semi-implicit time integration, the results are closer to the reference profile. Using $F_t = 0.4$, the difference between the algebraic and volume-conserving semi-implicit methods becomes visible. However, the difference between the algebraic and volume-conserving methods is still very small compared to the discrepancy between the temperature profile and the reference profile.

With a lower particle velocity (sub-variant P2-slow), the numerical error associated with convection is reduced and accurate results are obtained with a lower mesh resolution $m = 28.8$. Using $F_t = 0.05$, accurate results are obtained with the three cut-surface time-integration methods investigated (Fig. 10(a)). However, the algebraic and volume-conserving semi-implicit methods allow the simulation of the case with a larger value $F_t = 0.8$ of the time-stepping coefficient (Fig. 10(b)), whereas the temperature profile obtained with the explicit method is in that case disturbed near the interface.

With a larger particle velocity (sub-variant P2-fast), the effect of the cut-surface time-integration method is increased further. In that case,

the simulation is firmly out of temporal convergence as early as $F_t = 0.1$ (Fig. 11(c)). Along the x-axis direction, the poor temporal resolution of convective schemes flattens the temperature profile. With the explicit time-integration method, the effect extends up to and beyond interface. By contrast, with the algebraic and volume-conserving semi-implicit time-integration methods, the temperature profile is flattened near the interface while preserving a region of higher temperature at the centre of the particle. The effects depend on the direction considered, since the resolved temperature field is anisotropic, contrary to the reference solution. The temperature field for higher time-stepping coefficients is deformed along the x direction (Fig. 12(b, d, f)), which is the direction of motion at the time of snapshot. An accurate representation of the temperature profile along the x-axis is thus generally representative of an accurate simulation. With regard to the comparison between the algebraic and volume-conserving semi-implicit methods at $F_t = 0.1$, we may remark that the volume-conserving method provide slightly more accurate results. The difference between the two methods is however very small compared to the error with the reference temperature profile. Finally, as in the other cases, the difference between the three methods becomes smaller as the time-stepping coefficient $F_t$ is reduced, since the simulation approaches temporal convergence. This leads to a relatively accurate temperature profile at $F_t = 0.02$ (Fig. 11(a)) with
<figure><img src="image_11.png" /><figcaption>Fig. 8. Temperature profile along the x-axis for increasing mesh resolutions, in case P2, using $F_t = 0.05$, with the baseline one-field method and the cut-cell method with algebraic semi-implicit time integration. The solid line (—) is the reference static one-dimensional spherical solution.</figcaption></figure>

Baseline: —○— $m = 10.4$, - □- $m = 14.4$, △- $m = 20.8$, ▽- $m = 28.8$, ◄- $m = 41.6$;

Cut-cell: —●— $m = 10.4$, - ■- $m = 14.4$, ▲- $m = 20.8$, - ▽- $m = 28.8$, - ◄- $m = 41.6$.

<figure><img src="image_12.png" /><figcaption>Fig. 9. Effect of the cut-surface time-integration method of the cut-cell method on the temperature profile along the x-axis, in case P2, using a mesh resolution $m = 41.6$. — Reference; —●— Explicit; - ▽- Algebraic semi-implicit; - ◄- Volume-conserving semi-implicit.</figcaption></figure>

an isotropic distribution (Fig. 12(a, c, e)). To summarise these results, Fig. 13 shows the norm of the global error with the different time-integration methods as a function of the time-stepping coefficient $F_t$. Overall, the present cut-cell method has a linear order of convergence in the limit of small time steps for the prediction of the temperature field. This is consistent with the expected rate of convergence with

a forward Euler time scheme. The use of an explicit time-integration method increases the norm of the error, compared to an algebraic semi-implicit time-integration method, but does not modify the rate of convergence of the solution. Based on the linear limiting behaviour, the effect on the global error is comparable in magnitude to a halving of the simulation time step. For a larger time-stepping coefficient ($F_t >$
<figure><img src="image_13.png" /><figcaption>Fig. 10. Effect of the cut-surface time-integration method of the cut cell method on the temperature profile along the x-axis, in case P2-slow, using a mesh resolution m = 28.8. Reference; — Explicit; - - - Algebraic semi-implicit; ◦ - Volume-conserving semi-implicit.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Effect of the cut-surface time-integration method of the cut cell method on the temperature profile along the x-axis, in case P2-fast, using a mesh resolution m = 41.6. Reference; — Explicit; - - - Algebraic semi-implicit; ◦ - Volume-conserving semi-implicit.</figcaption></figure>

0.1), the improvement provided by the algebraic semi-implicit time integration might be larger than suggested by this limiting behaviour.

The improvements obtained with a semi-implicit time-integration method, compared to an explicit time-integration method, coincide with a more accurate satisfaction of the dynamic conservation of volume (35). Fig. 14 compares the error made on dynamic volume conservation with the explicit and algebraic semi-implicit time-integration methods, which do not ensure this property. This is evaluated using

the mean relative error over cut cells of the dynamic conservation of volume (35),

$$ \epsilon_V = \frac{1}{\sum_n N_{\text{cut}}(t_n)} \times \sum_n \sum_{i,j,k} \frac{\|V_{i,j,k}^{\psi}(t_{n+1}) - V_{i,j,k}^{\psi}(t_n) - \sum_F \in \tilde{\mathcal{S}}_{i,j,k}^{\psi} \overline{S}_F(t_n) u_F(t_n) \cdot n_F \Delta t\|}{\|V_{i,j,k}^{\psi}(t_{n+1}) - V_{i,j,k}^{\psi}(t_n)\|}, \quad (38) $$
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Nondimensionalised temperature $T^*$ on the plane $y = 0$ with the cut-cell method for different cut-surface time-integration methods, in case P2-fast, using a mesh resolution *m* = 41.6. The black circle indicates the interface between the two phases. Contour lines are incremented by 0.05.</figcaption></figure>

where $N_{\text{cut}}(t_n)$ is the total number of cut cells at time $t_n$. The results show that volume-conservation error is relatively independent of the mesh resolution for a given time-stepping coefficient $F_t$ and increases with $F_t$ at a mildly superlinear rate. For a given time-stepping coefficient $F_t$, the algebraic semi-implicit time-integration of cut-surfaces

improves the relative error on dynamic volume conservation by almost one order of magnitude compared to the explicit time integration. With the algebraic semi-implicit time-integration method, the relative volume-conservation error is above 20% at $F_t = 0.4$, whereas this threshold is reached near $F_t = 0.075$ with the explicit time-integration
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Evolution of the global error of the simulation with the time-stepping coefficient $F_t$, with the cut-cell method, in case P2-fast, using a mesh resolution $m = 41.6$. The global error corresponds to the mean squared difference, over all computational cells, between the temperature field of the simulation and the temperature field of the corresponding simulation without convection (in the static pure-diffusion case) using $F_t = 0.2$. The dash-dotted lines correspond to a linear growth. • Explicit, ▽ Algebraic semi-implicit, ◆ Conservative semi-implicit.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Evolution of the mean relative volume-conservation error $\epsilon_v$ of cut cells with the time-stepping coefficient $F_t$, with the cut-cell method, in case P2-fast, using mesh resolutions ranging from $m = 7.2$ to $m = 41.6$. Blue symbols (○, ▼, ▽, △) are used for the explicit time-integration method, whereas magenta symbols (□, △, ▲, ●, ○) are used for the algebraic semi-implicit time-integration method. The dash-dotted lines correspond to a linear growth. The dashed lines are power-law fit with exponent 1.2.</figcaption></figure>

method. Since the dynamic volume conservation equation corresponds to the transport equation of a uniform scalar field, we might expect that the discrepancy affects to the accuracy of convective transfers. Negligible further improvements are obtained in practice with the volume-conserving semi-implicit time-integration method, which enforces the dynamic volume conservation by construction. This demonstrates that the dynamic volume conservation is satisfied with sufficient accuracy with the algebraic semi-implicit time-integration method, in all the examined cases that are also able to resolve convective and diffusive transfers with some accuracy.

## 4. Conclusion

The time integration of cut surfaces in cut-cell methods should account for the displacement of the interface during a simulation time step. Assuming a given displacement of the interface, the area of cut faces resulting from the intersection of the background grid and the interface boundary may be provided by a semi-implicit method. This approximates the effective area of cut faces during a simulation time

step, and as a result satisfies more accurately the dynamic conservation of volume in cut cells. The practical effect of the cut-surface time-integration method is demonstrated in a cut-cell method for the convection-diffusion of a passive scalar. The cut-cell method combines the explicit or semi-implicit treatment of cut surfaces with an explicit time-integration of fluxes, and an explicit displacement of the interface. The finite-volume approach, on a staggered Cartesian grid system, ensures strict conservation properties and maintains numerical stability using a flux-redistribution strategy. The simulation of heat diffusion from an initially hot spherical particle into the surrounding fluid, with a cyclic ellipsoidal apparent motion, validates the accuracy of the cut-cell method compared to a baseline one-field method. The time-integration method for cut surfaces may have a significant effect on the numerical prediction provided that the numerical Peclet number is larger than one. Compared to an explicit cut-surface time integration, the algebraic and volume-conserving semi-implicit time-integration methods may increase the critical time step of temporal convergence, hence improving the computational efficiency of the simulation. These results provide insight for the development of cut-cell methods in moving rigid bodies and two-phase flows with complex deformable interfaces, such as bubbles.

## CRediT authorship contribution statement

**D. Dupuy:** Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing. **A. Toutant:** Conceptualization, Investigation, Methodology, Writing – review & editing. **A. du Cluzeau:** Conceptualization, Investigation, Methodology, Writing – review & editing. **G. Bois:** Conceptualization, Investigation, Methodology, Writing – review & editing.

## 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

This project was provided with computer and storage resources by GENCI, France at TGCC thanks to the grant A0162B07712 on the supercomputer Joliot Curie's ROME partition.

## Data availability

Data will be made available on request.

## References

[1] Ishii M, Mishima K. Two-fluid model and hydrodynamic constitutive relations. Nucl Eng Des 1984;82(2-3):107-26.

[2] Tsuji Y, Kawaguchi T, Tanaka T. Discrete particle simulation of two-dimensional fluidized bed. Powder Technol 1993;77(1):79-87.

[3] Crowe CT, Troutt TR, Chung JN. Numerical models for two-phase turbulent flows. Annu Rev Fluid Mech 1996;(28):11-43.

[4] Sommerfeld M. Numerical methods for dispersed multiphase flows. Part Flows 2017;327-96.

[5] Hirt CW, Amsden AA, Cook JL. An arbitrary Lagrangian-Eulerian computing method for all flow speeds. J Comput Phys 1974;14(3):227-53.

[6] Hu HH, Patankar NA, Zhu MY. Direct numerical simulations of fluid-solid systems using the arbitrary Lagrangian-Eulerian technique. J Comput Phys 2001;169(2):427-62.

[7] Donea J, Huerta A, Ponthot J-P, Rodríguez-Ferran A. Arbitrary Lagrangian-Eulerian methods. Encycl Comput Mech 2004.

[8] Quan S, Schmidt DP. A moving mesh interface tracking method for 3D incompressible two-phase flows. J Comput Phys 2007;221(2):761-80.

[9] Anderson RW, Dobrev VA, Kolev TV, Rieben RN, Tomov VZ. High-order multi-material ALE hydrodynamics. SIAM J Sci Comput 2018;40(1):B32-58.
[10] Anjos GR, Mangiavacchi N, Thome JR. An ALE-FE method for two-phase flows with dynamic boundaries. Comput Methods Appl Mech Engrg 2020;362:112820.

[11] Haeri S, Shrimpton JS. On the application of immersed boundary, fictitious domain and body-conformal mesh methods to many particle multiphase flows. Int J Multiph Flow 2012;40:38–55.

[12] Mittal R, Iaccarino G. Immersed boundary methods. Annu Rev Fluid Mech 2005;37:239–61.

[13] Dolbow J, Farhat C, Harari I, Lew A. Special issue: Advances in embedded interface methods. Internat J Numer Methods Engrg 2015;104(7):469–71.

[14] Unverdi SO, Tryggvason G. A front-tracking method for viscous, incompressible, multi-fluid flows. J Comput Phys 1992;100(1):25–37.

[15] Tryggvason G, Bunner B, Esmaeli A, Juric D, Al-Rawahi N, Tauber W, et al. A front-tracking method for the computations of multiphase flow. J Comput Phys 2001;169(2):708–59.

[16] Koffi Bi D-A, Tavares M, Chénier E, Vincent S. Accuracy and convergence of the curvature and normal vector discretizations for 3D static and dynamic front-tracking interfaces. J Comput Phys 2022;461:111197.

[17] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39(1):201–25.

[18] Rider WJ, Kothe DB. Reconstructing volume tracking. J Comput Phys 1998;141(2):112–52.

[19] Scapin N, Costa P, Brandt L. A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows. J Comput Phys 2020;407:109251.

[20] Dyadechko V, Shashkov M. Moment-of-fluid interface reconstruction. Los Alamos Report LA-UR-05-7571, 49, 2005.

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

[22] Dervieux A, Thomasset F. A finite element method for the simulation of a Rayleigh-Taylor instability. Approx Methods Navier-Stokes Probl 1980;145–58.

[23] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. J Comput Phys 1988;79(1):12–49.

[24] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114(1):146–59.

[25] Gibou F, Fedkiw R, Osher S. A review of level-set methods and some recent applications. J Comput Phys 2018;353:82–109.

[26] Yabe T, Xiao F, Utsumi T. The constrained interpolation profile method for multiphase analysis. J Comput Phys 2001;169(2):556–93.

[27] Fix GJ. Phase field methods for free boundary problems. In: Free boundary problems: theory and applications. vol. 1, Pitman Boston; 1983, p. 580–9.

[28] Langer J. Models of pattern formation in first-order phase transitions. In: Distributions in condensed matter physics: memorial volume in honor of shang-keung ma. World Scientific; 1986, p. 165–86.

[29] Anderson DM, McFadden GB, Wheeler AA. Diffuse-interface methods in fluid mechanics. Annu Rev Fluid Mech 1998;30(1):139–65.

[30] Li J, Zheng D, Zhang W. Advances of phase-field model in the numerical simulation of multiphase flows: A review. Atmosphere 2023;14(8):1311.

[31] Tryggvason G, Bunner B, Ebrat O, Tauber W. Computations of multiphase flows by a finite difference/front tracking method. I. Multi-fluid flows. In: Lecture series-von Karman institute for fluid dynamics, Von Karman Institute for Fluid Dynamics; 1998, p. 7.

[32] Fuster D, Mimoh Y. Regularization errors introduced by the one-fluid formulation in the solution of two-phase elliptic problems. J Comput Phys 2024;113202.

[33] Prosperetti A, Tryggvason G. Computational methods for multiphase flow. Cambridge University Press; 2009.

[34] Tenneti S, Subramaniam S. Particle-resolved direct numerical simulation for gas-solid flow model development. Annu Rev Fluid Mech 2014;46(1):199–230.

[35] Du Cluzeau A, Bois G, Leoni N, Toutant A. Analysis and modeling of bubble-induced agitation from direct numerical simulation of homogeneous bubbly flows. Phys Rev Fluids 2022;7(4):044604.

[36] Hamidi MS, Toutant A, Mer S, Bataille F. Assessment of a coupled VOF-Front-Tracking/DEM method for simulating fluid-particles flows. Int J Multiph Flow 2023;165:10467.

[37] Fedkiw RP, Aslam T, Merriman B, Osher S. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J Comput Phys 1999;152(2):457–92.

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

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

[40] Bois G, du Cluzeau A. DNS of turbulent bubbly flows in plane channels using the Front-Tracking algorithm of TrioCFD. In: Fluids engineering division summer meeting. Vol. 58066, American Society of Mechanical Engineers; 2017, V01CT16A005.

[41] Grosso M, Bois G, Toutant A. Thermal boundary layer modelling for heat flux prediction of bubbles at saturation: A priori analysis based on fully-resolved simulations. Int J Heat Mass Transfer 2024;222:124980.

[42] Lopez SM, Christou A, Pan S, Stoesser T, Xie Z. A new ghost-cell/level-set method for three-dimensional flows. J Comput Phys 2024;499:112710.

[43] Clarke DK, Salas MD, Hassan HA. Euler calculations for multielement airfoils using Cartesian grids. AIAA J 1986;24(3):353–8.

[44] Berger M, Leveque R. An adaptive Cartesian mesh algorithm for the Euler equations in arbitrary geometries. In: 9th computational fluid dynamics conference. 1989, p. 1930.

[45] Udaykumar HS, Shyy W, Rao MM. Elafint: a mixed Eulerian-Lagrangian method for fluid flows with complex and moving boundaries. Internat J Numer Methods Fluids 1996;22(8):691–712.

[46] Mittal R, Seo JH. Origin and evolution of immersed boundary methods in computational fluid dynamics. Phys Rev Fluids 2023;8(10):100501.

[47] Chiang Y-L, van Leer B, Powell KG. Simulation of unsteady inviscid flow on an adaptively refined Cartesian grid. In: 30th aerospace sciences meeting and exhibit. 1992, p. 443.

[48] Udaykumar HS, Mittal R, Rampunggoon P, Khanna A. A sharp interface Cartesian grid method for simulating flows with complex moving boundaries. J Comput Phys 2001;174(1):345–80.

[49] Lahur PR, Nakamura Y. Flow calculation of moving body on Cartesian grid by using new cell-merging method. Trans Japan Soc Aeronaut Space Sci 2001;44(145):171–8.

[50] Ingram DM, Causon DM, Mingham CG. Developments in Cartesian cut cell methods. Math Comput Simulation 2003;61(3–6):561–72.

[51] Chung M-H. Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape. Comput & Fluids 2006;35(6):607–23.

[52] Ji H, Lien F-S, Yee E. Numerical simulation of detonation using an adaptive Cartesian cut-cell method combined with a cell-merging technique. Comput & Fluids 2010;39(6):1041–57.

[53] Kirkpatrick MP, Armfield SW, Kent JH. A representation of curved boundaries for the solution of the Navier-Stokes equations on a staggered three-dimensional Cartesian grid. J Comput Phys 2003;184(1):1–36.

[54] Hartmann D, Meinke M, Schröder W. An adaptive multilevel multigrid formulation for Cartesian hierarchical grid methods. Comput & Fluids 2008;37(9):1103–25.

[55] Hartmann D, Meinke M, Schröder W. A strictly conservative Cartesian cut-cell method for compressible viscous flows on adaptive grids. Comput Methods Appl Mech Engrg 2011;200(9–12):1038–52.

[56] Berger M, Giuliani A. A state redistribution algorithm for finite volume schemes on cut cell meshes. J Comput Phys 2021;428:109820.

[57] Giuliani A, Almgren AS, Bell JB, Berger MJ, de Frahan MTH, Rangarajan D. A weighted state redistribution algorithm for embedded boundary grids. J Comput Phys 2022;464:111305.

[58] Berger M, Giuliani A. A new provably stable weighted state redistribution algorithm. SIAM J Sci. Comput. 2024;46(5):A2848–73.

[59] Berger MJ, LeVeque RJ. Stable boundary conditions for Cartesian grid calculations. Comput Syst Eng 1990;1(2–4):305–11.

[60] Berger MJ, LeVeque RJ. A rotated difference scheme for Cartesian grids in complex geometries. AIAA Paper CP-91-1602, 1991, p. 1–9.

[61] Helzel C, Berger MJ, LeVeque RJ. A high-resolution rotated grid method for conservation laws with embedded geometries. SIAM J Sci Comput 2005;26(3):785–809.

[62] Berger M, Helzel C. A simplified h-box method for embedded boundary grids. SIAM J Comput 2012;34(2):A861–88.

[63] Cheny Y, Botella O. The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties. J Comput Phys 2010;229(4):1043–76.

[64] Nikfarjam F, Cheny Y, Botella O. The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries. Comput Phys Comm 2018;226:67–80.

[65] Engwer C, May S, Nüßing A, Streitbürger F. A stabilized DG cut cell method for discretizing the linear transport equation. SIAM J Sci Comput 2020;42(6):A3677–703.

[66] May S, Streitbürger F. DoD stabilization for non-linear hyperbolic conservation laws on cut cell meshes in one dimension. Appl Math Comput 2022;419:126854.

[67] Chern I-L, Colella P. A conservative front tracking method for hyperbolic conservation laws. LLNL Rep. No. UCRL-97200, 51, Lawrence Livermore National Laboratory; 1987, p. 83–110.

[68] Bell J, Welcome M, Colella P. Conservative front-tracking for inviscid compressible flow. In: 10th computational fluid dynamics conference. 1991, p. 1599.

[69] Pember RB, Bell JB, Colella P, Curtchfield WY, Welcome ML. An adaptive Cartesian grid method for unsteady compressible flow in irregular regions. J Comput Phys 1995;120(2):278–304.

[70] Colella P, Graves DT, Keen BJ, Modiano D. A Cartesian grid embedded boundary method for hyperbolic conservation laws. J Comput Phys 2006;211(1):347–66.

[71] Graves D, Colella P, Modiano D, Johnson J, Sjogreen B, Gao X. A cartesian grid embedded boundary method for the compressible Navier-Stokes equations. Commun Appl Math Comput Sci 2013;8(1):99–122.
[72] Schneiders L, Hartmann D, Meinke M, Schröder W. An accurate moving boundary formulation in cut-cell methods. J Comput Phys 2013;235:786-809.

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

[74] May S, Berger M. An explicit implicit scheme for cut cells in embedded boundary meshes. J Sci Comput 2017;71:919-43.

[75] May S, Laakmann F. Accuracy analysis for explicit-implicit finite volume schemes on cut cell meshes. Commun Appl Math Comput 2024;1-26.

[76] Xie Z. An implicit Cartesian cut-cell method for incompressible viscous flows with complex geometries. Comput Methods Appl Mech Engrg 2022;399:115449.

[77] Xie Z, Lin P, Stoesser T. A conservative and consistent implicit Cartesian cut-cell method for moving geometries with reduced spurious pressure oscillations. J Comput Phys 2022;459:111124.

[78] Sverdrup K, Almgren A, Nikiforakis N. An embedded boundary approach for efficient simulations of viscoplastic fluids in three dimensions. Phys Fluids 2019;31(9).

[79] Fröhlich K, Meinke M, Schröder W. Correlations for inclined prolates based on highly resolved simulations. J Fluid Mech 2020;901:A5.

[80] Wang T, Song Z, Zhou G, Jiang C, Shi F. A cell-based smoothed finite element method for three-dimensional incompressible flows using Cartesian cut-cell meshes. Phys Fluids 2024;36(6).

[81] Falcovitz J, Alfandary G, Hanoch G. A two-dimensional conservation laws scheme for compressible flows with moving boundaries. J Comput Phys 1997;138(1):83-102.

[82] Yang G, Causon DM, Ingram DM. Calculation of compressible flows about complex moving geometries using a three-dimensional Cartesian cut cell method. Internat J Numer Methods Fluids 2000;33(8):1121-51.

[83] Schneiders L, Günther C, Meinke M, Schröder W. An efficient conservative cut-cell method for rigid bodies interacting with viscous compressible flows. J Comput Phys 2016;311:62-86.

[84] Bennett WP, Nikiforakis N, Klein R. A moving boundary flux stabilization method for Cartesian cut-cell grids using directional operator splitting. J Comput Phys 2018;368:333-58.

[85] Ghigo AR, Popinet S, Wachs A. A conservative finite volume cut-cell method on an adaptive Cartesian tree grid for moving rigid bodies in incompressible flows. 2021.

[86] Bai X, Deng X. A sharp interface method for compressible multi-phase flows based on the cut cell and ghost fluid methods. Adv Appl Math Mech 2017;9(5):1052-75.

[87] Bai X, Li M. Simulating compressible two-phase flows with sharp-interface discontinuous Galerkin methods based on ghost fluid method and cut cell scheme. J Comput Phys 2022;459:111107.

[88] Bai X, Li M. A conservative sharp-interface numerical method for two-dimensional compressible two-phase flows. J Sci Comput 2023;97(2):30.

[89] Hu XY, Khoo BC, Adams NA, Huang FL. A conservative interface method for compressible flows. J Comput Phys 2006;219(2):553-78.

[90] Shen Y, Ren Y, Ding H. A 3D conservative sharp interface method for simulation of compressible two-phase flows. J Comput Phys 2020;403:109107.

[91] Pang B, Ren Y, Shen Y, Liu H-R, Ding H. A conservative sharp interface method for two-dimensional incompressible two-phase flows with phase change. Phys Fluids 2023;35(12).

[92] Pang B, Ren Y, Shen Y, Liu H-R, Ding H. A conservative sharp interface method for incompressible viscous two-phase flows with topology changes and large density ratio. Comput & Fluids 2024;274:106212.

[93] Lin J-Y, Shen Y, Ding H, Liu N-S, Lu X-Y. Simulation of compressible two-phase flows with topology change of fluid-fluid interface by a robust cut-cell method. J Comput Phys 2017;328:140-59.

[94] Deng X-L, Li M. Simulating compressible two-medium flows with sharp-interface adaptive Runge-Kutta discontinuous Galerkin methods. J Sci Comput 2018;74(3):1347-68.

[95] Chang C-H, Deng X, Theofanous TG. Direct numerical simulation of interfacial instabilities: a consistent, conservative, all-speed, sharp-interface method. J Comput Phys 2013;242:946-90.

[96] Thomas PD, Lombard CK. Geometric conservation law and its application to flow computations on moving grids. AIAA J 1979;17(10):1030-7.

[97] Zhang H, Reggio M, Trepanier JY, Camarero R. Discrete form of the GCL for moving meshes and its implementation in CFD schemes. Comput & Fluids 1993;22(1):9-23.

[98] Guillard H, Farhat C. On the significance of the geometric conservation law for flow computations on moving meshes. Comput Methods Appl Mech Engrg 2000;190(11-12):1467-82.

[99] Chang X, Ma R, Zhang L, He X, Li M. Further study on the geometric conservation law for finite volume method on dynamic unstructured mesh. Comput & Fluids 2015;120:98-110.

[100] Ma R, Chang X, Zhang L, He X, Li M. On the geometric conservation law for unsteady flow simulations on moving mesh. Procedia Eng 2015;126:639-44.

[101] Wang Z, Yang H. Unsteady flow simulation using a zonal multi-grid approach with moving boundaries. In: 32nd aerospace sciences meeting and exhibit. 1994, p. 57.

[102] Mavripis DJ, Yang Z. Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes. J Comput Phys 2006;213(2):557-73.

[103] Étienne S, Garon A, Pelletier D. Perspective on the geometric conservation law and finite element methods for ALE simulations of incompressible flow. J Comput Phys 2009;228(7):2313-33.

[104] Farhat C, Geuzaine P, Grandmont C. The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids. J Comput Phys 2001;174(2):669-94.

[105] Mulbah C, Kang C, Mao N, Zhang W, Shaikh AR, Teng S. A review of VOF methods for simulating bubble dynamics. Prog Nucl Energy 2022;154:104478.

[106] Roenby J, Bredmose H, Jasak H. A computational method for sharp interface advection. R Soc Open Sci 2016;3(11):160405.

[107] Owkes M, Desjardins O. A computational framework for conservative, three-dimensional, unsplit, geometric transport with application to the volume-of-fluid (VOF) method. J Comput Phys 2014;270:587-612.

[108] Ivey CB, Moin P. Conservative and bounded volume-of-fluid advection on unstructured grids. J Comput Phys 2017;350:387-419.

[109] Marić T, Marschall H, Bothe D. An enhanced un-split face-vertex flux-based VoF method. J Comput Phys 2018;371:967-93.

[110] Marić T, Kothe DB, Bothe D. Unstructured un-split geometrical volume-of-fluid methods-a review. J Comput Phys 2020;420:109695.

[111] Mathieu B. Demonstration of a 3D parallel implementation of the Front-tracking method: Simulation of fluid mixing with a moving boundary and a free surface. In: Proceedings of the advances in the modelling methodologies of two-phase flows meeting. Lyon, France; 2004.

[112] Leonard BP. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput Methods Appl Mech Engrg 1979;19(1):59-98.

[113] Portelenelle B, Botella O, Cheny Y. Accurate discretization of diffusion in the LS-STAG cut-cell method using diamond cell techniques. Comput & Fluids 2019;189:34-45.

[114] Giraldo FX, Neta B. A comparison of a family of Eulerian and semi-Lagrangian finite element methods for the advection-diffusion equation. WIT Trans Built Environ 1997;30.

[115] Bonaventura L. An introduction to semi-Lagrangian methods for geophysical scale flows. In: Lecture notes, ERCOFTAC Leonhard Euler Lectures, SAM-ETH Zurich, Vol. 140, Citeseer; 2004.

[116] Lentine M, Grétarsson JT, Fedkiw R. An unconditionally stable fully conservative semi-Lagrangian method. J Comput Phys 2011;230(8):2857-79.

[117] Liu L, Becerra M. An efficient semi-Lagrangian algorithm for simulation of corona discharges: the position-state separation method. IEEE Trans Plasma Sci 2016;44(11):2822-31.

[118] Ingelsten S, Mark A, Kádár R, Edelvik F. A backwards-tracking Lagrangian-Eulerian method for viscoelastic two-fluid flows. Appl Sci 2021;11(1):439.

[119] Du Cluzeau A, Bois G, Toutant A, Martinez J-M. On bubble forces in turbulent channel flows from direct numerical simulations. J Fluid Mech 2020;882:A27.

[120] Calvin C, Cueto O, Emonot P. An object-oriented approach to the design of fluid mechanics software. ESAIM Math Model Numer Anal 2002;36(05):907-21.

[121] Chandesris M, d'Hueppe A, Mathieu B, Jamet D, Goyeau B. Direct numerical simulation of turbulent heat transfer in a fluid-porous domain. Phys Fluids 2013;25(12):125110.

[122] Dupuy D, Toutant A, Bataille F. Turbulence kinetic energy exchanges in flows with highly variable fluid properties. J Fluid Mech 2018;834:5-54.

[123] Sonolet A, Grosso M, Bois G, Toutant A. Conservative formulations for heat transfer simulation with sharp interfaces and discontinuous properties in two-phase flows. Comput & Fluids 2023;254:105798.