# Subcycling Strategy for Finite-Volume Updated-Lagrangian Methods Applied to Fluid-Structure Interaction

Teddy Chantrait, Nicolas Chevaugeon, Stéphane del Pino, Alexandre Gangloff, Emmanuel Labourasse

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

## ► To cite this version:

Teddy Chantrait, Nicolas Chevaugeon, Stéphane del Pino, Alexandre Gangloff, Emmanuel Labourasse. *Subcycling Strategy for Finite-Volume Updated-Lagrangian Methods Applied to Fluid-Structure Interaction*. International Journal for Numerical Methods in Engineering, 2025, 126 (11), pp.e70051. 10.1002/nme.70051. *hal-05237182*

HAL Id: hal-05237182

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

Submitted on 9 Sep 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
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

OPEN ACCESS

# Subcycling Strategy for Finite-Volume Updated-Lagrangian Methods Applied to Fluid-Structure Interaction

Teddy Chantrait<fnref n="1" />,² ID | Nicolas Chevaugeon³ | Stéphane Del Pino<fnref n="1" />,² | Alexandre Gangloff<fnref n="1" />,³ | Emmanuel Labourasse<fnref n="1" />,²

**Correspondence:** Alexandre Gangloff (alexandre.gangloff@cea.fr)

**Received:** 11 March 2025 | **Revised:** 11 March 2025 | **Accepted:** 8 May 2025

**Keywords:** computational efficiency | finite-volume method | fluid-structure interaction (FSI) | partitioned method | shock dynamics | subcycling | updated-Lagrangian method

## ABSTRACT

In this article, we propose and investigate an explicit partitioned method for solving shock dynamics in fluid-structure interaction (FSI) problems. The method is fully conservative, ensuring the local conservation of mass, momentum, and energy, which is crucial for accurately capturing strong shock interactions. Using an updated-Lagrangian finite-volume approach, the method integrates a subcycling strategy to decouple time steps between the fluid and structure, significantly enhancing computational efficiency. Numerical experiments confirm the accuracy and stability of the method, demonstrating that it retains the key properties of monolithic solvers while reducing computational costs. Extensive validation across 1D and 3D FSI problems shows the method's capability for large-scale, fast transient simulations, making it a promising solution for high-performance applications.

## 1 | Introduction

Fluid-structure interaction (FSI) is a complex phenomenon of significant importance in numerous fields of engineering and science, where the interplay between a fluid and a solid structure is crucial for understanding system dynamics. Modeling and simulating such interactions are challenging due to the inherently multi-physical nature of the problem, where both fluid and solid domains must be addressed simultaneously, each governed by distinct behaviors and equations. Applications range from aerodynamics, where airflow interacts with aircraft structures, to biomechanics, where blood flow affects arterial walls. These scenarios necessitate sophisticated and computationally efficient methods to ensure accurate simulations (refer, for instance, to [1] for a review).

Previous research has tackled the complexities of FSI using various modeling approaches, including both monolithic

and partitioned frameworks. In the monolithic approach, the fluid and structure are treated as a unified system and solved concurrently. This strategy, discussed extensively in seminal works [2–5], ensures strong coupling and accurate interaction modeling, particularly in cases involving substantial structural deformation or when the fluid and solid share similar time scales. However, monolithic methods often become computationally prohibitive, especially for large-scale simulations or systems involving multiple interacting structures with vastly differing physical properties.

Partitioned methods, as described in pioneering works such as [3, 4, 6–9], provide more flexibility by solving the fluid and solid equations independently, with coupling conditions applied at the interface. These methods are advantageous because they allow the use of specialized solvers for the fluid and solid components, enhancing computational efficiency. However, maintaining stability and convergence in partitioned methods

This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.

© 2025 The Author(s). *International Journal for Numerical Methods in Engineering* published by John Wiley & Sons Ltd.

<footnote n="1">¹CEA, DAM, DIF, Arpajon, France | ²Université Paris-Saclay, CEA DAM DIF, Laboratoire en Informatique Haute Performance Pour le Calcul et la Simulation, Arpajon, France | ³Université Nantes I, GeM-Institut de Recherche en Génie Civil et Mécanique, Nantes, France</footnote>
frequently requires special techniques to prevent numerical instabilities.

One of the first widely explored strategies permitted by partitioned methods is the use of adaptive time-subcycling, particularly for problems where the material properties or physical behaviors vary significantly between different regions of the computational domain. Subcycling methods are especially effective for managing computational efficiency in multiscale problems, where different regions exhibit distinct time scales. Early works by Belytschko et al. [10, 11] introduced subcycling in the context of explicit time-stepping to manage computational costs associated with refined meshes. Implicit-explicit coupling [12] was also studied very early. More recent works on this subject include [13–17].

In this study, we aim to adapt this kind of strategy (a partitioned framework with time-subcycling) to our specific context involving an explosion interacting with thin structures. Consequently, special attention is given to ensuring the conservation of mass, momentum, and total energy, as well as thermodynamic consistency (e.g., entropy production), to effectively capture strong shock phenomena. This motivated our choice of models: the Euler equations for the fluid and a hyperelastic model for the solid. To solve these models, we leverage recent numerical methods developed within the updated-Lagrangian framework (refer to [18–25]). Other recent works that address FSI with conservative and entropy-stable solvers include [26]. Earlier efforts in this direction include [27–31], which use Godunov schemes to capture FSI in high strain-rate problems. All these studies rely on an Eulerian or direct Arbitrary-Lagrangian-Eulerian (ALE) framework.

In contrast, our work uses an updated-Lagrangian or indirect ALE framework. The updated-Lagrangian approach prevents material fluxes at the discrete level, eliminating the need for mixed-cell closure models. The nodal Riemann solvers we use possess favorable mathematical properties, such as conservation, stability, and consistency. Moreover, the hyperelastic solvers we use [18, 19, 23, 24] can be viewed as generalizations of their fluid counterparts [20, 21, 25], which facilitates the development of a monolithic FSI solver that is locally conservative, even at the interface.

In this article, we propose a subcycling strategy tailored to fast transient problems for these types of schemes. We found very few works in the literature that specifically address this aspect. Brun et al. [32] propose an implicit-explicit time-step management, observing computation time reduction of 60% (see also [33–38] for seminal works in this direction). In the article [39], Faucher explains how to extend this type of methods to contact. All these methods deal with finite-element, and conservation is not strongly enforced.

Our proposed method achieves conservation while using different time steps for the fluid and the solid, with the only requirement being a semi-conformal discretization in time. In other words, each discrete time value $t_f^n$ used for the fluid part must coincide with a discrete time value $t_s^m$ in the solid ($\forall n > 0, \exists m$, such that $t_s^m = t_f^n$). Furthermore, consistency and stability are preserved by our method. In this article, we propose

and investigate an explicit partitioned method to solve shock dynamics in FSI scenarios. The remaining of the article articulates as follows. In Section 2 we describe briefly the model. Then in Section 3 we recall the updated-Lagrangian finite volume nodal solver construction and properties. This is useful to study the proposed partitioned method in Section 4. In Section 5 we present the time-subcycling strategy. This is the main novelty of this article. In short, the main ingredient of the method, which is both simple and effective, consists of fixing the fluxes during the time-subcycling. This approach ensures that all the key properties of the underlying method—conservation and stability (entropy growth)—are preserved. The expected behavior of the method is illustrated through numerical tests in Section 6. We achieve a computational time speedup factor of 8 on the final test problem. Finally, we conclude and draw some perspectives.

## 2 | Model

This section follows the books of Gurtin et al. [40] and Godunov and Romenskii [41].

* Vectors are displayed in bold font x.

* Second-order tensors are displayed in blackboard-bold font X.

* The euclidean norm of x is denoted by $|x|$.

* The scalar product of x and y is denoted by $x \cdot y$.

* The tensor product of x and y is denoted by $x \otimes y$.

* The inner product of second-order tensors X and Y is denoted by $X : Y$.

* $X^T$ is the transposed tensor of X.

* $X^* = X - \frac{1}{d} \text{tr}(X) \mathbf{1}_d$ is the deviatoric part of the tensor X, $\mathbf{1}_d$ the identity tensor and d the spatial dimension of the problem (1, 2 or 3 in our case).

Let $t \in \mathbb{R}^+$ and let $\Omega(t) \subset \mathbb{R}^d$ a Lagrangian open set. We consider the following system of conservation laws

$$ \forall \mathbf{x} \in \Omega, \quad \left\{ \begin{array}{l} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \\ \frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) = \nabla \cdot \mathcal{S}, \\ \frac{\partial \rho E}{\partial t} + \nabla \cdot (\rho E \mathbf{v}) = \nabla \cdot (\mathcal{S} \mathbf{v}), \end{array} \right. \quad (1) $$

with unknowns $\rho(\mathbf{x}, t)$ the density, $v(\mathbf{x}, t)$ the velocity and $E(\mathbf{x}, t)$ the specific total energy. This system is supplemented by valid boundary conditions on $\partial\Omega$, the boundary of $\Omega$. To close the system, one needs to define $\mathcal{S}$, the Cauchy-Stress tensor. For the hyperelastic model, it is shown that $\mathcal{S}$ is a function of $\eta$ and $\mathcal{F}$ (refer for instance to [40, 42]), $\eta$ being the specific entropy and $\mathcal{F}$ the deformation gradient. The deformation gradient $\mathcal{F} = \nabla_{\xi}\mathbf{x}$ is the Jacobian matrix of the transformation mapping the initial configuration $\xi \in \Omega(0)$ to the current configuration $\mathbf{x}(\xi, t) \in \Omega(t)$. To ease the notation, we will use $\nabla \equiv \nabla_{\mathbf{x}}$, meaning that the spatial derivatives are with respect to the current configuration. In the same regard, we omit the time dependency for the
computational domain, writing simply $\Omega$. Hyperelastic models rely to the following modeling relation

$$ \mathcal{S}(\eta, \mathbb{F}) = \rho \nabla_{\mathbb{F}} e^{\mathcal{F}^T} \quad (2) $$

where $e(\eta, \mathbb{F}) = E - \frac{|\mathbb{F}|^2}{2}$ is the internal energy. The entropy $\eta$ satisfies the following Gibbs law

$$ T d\eta = d\epsilon - \nabla_{\mathbb{F}} e : d\mathbb{F} \quad (3) $$

where by definition, $T = \frac{\partial e(\eta, \mathbb{F})}{\partial \eta}$ is the temperature. It can be shown that $-\eta$ is indeed a mathematical entropy for the system (1), that is

$$ \frac{\partial \rho \eta}{\partial t} + \nabla \cdot (\rho \mathbf{v} \eta) \geq 0 \quad (4) $$

with an equality for smooth solutions.

To close the system, one needs a supplementary equation on $\mathbb{F}$

$$ \frac{d}{dt} \mathbb{F} = \nabla_{\mathbb{F}} \mathbb{F} \quad (5) $$

In this article, we will focus on isotropic materials. In this case, $e(\eta, \mathbb{F})$ can be defined from the invariants of the left Cauchy-Green tensor $\mathbb{B} = \mathbb{F}\mathbb{F}^T$. The closure law (2) can be rewritten

$$ \mathcal{S}(\eta, \mathbb{B}) = 2\rho \nabla_{\mathbb{B}} e \mathbb{B} \quad (6) $$

and considering now that $e(\eta, \mathbb{B})$, Gibbs relation reads

$$ T d\eta = d\epsilon - \nabla_{\mathbb{B}} e : d\mathbb{B} \quad (7) $$

Plugging (6) into (7) gives eventually

$$ T d\eta = d\epsilon - \frac{1}{2\rho} \mathcal{S}\mathbb{B}^{-1} : d\mathbb{B} \quad (8) $$

Moreover, we assume that the internal energy can be put in the following separate form

$$ e(\eta, \mathbb{F}) = e_v(\eta, \mathbf{J}) + e_s(I_1, I_2) \quad (9) $$

where $\mathbf{J} = \det(\mathbb{F})$, and where $I_1 = \text{tr}(\overline{\mathbb{B}})$ and $I_2 = \frac{1}{2}(\text{tr}(\overline{\mathbb{B}})^2 - \text{tr}(\overline{\mathbb{B}}^2))$ are two invariants of $\overline{\mathbb{B}} = J^{-2/3}\mathbb{B}$. In this formulation, one separates the hydrodynamic internal energy due to the volume change $e_v(\eta, \mathbf{J})$ from the isochoric deformation energy $e_s(I_1, I_2)$. Putting this decomposition into (6) yields

$$ \mathcal{S} = \rho \mathbf{J} \frac{\partial e_v}{\partial \mathbb{F}}|_{\mathbb{F}} + 2\rho \left( \frac{\partial e_s}{\partial I_1} \overline{\mathbb{B}}^* - \frac{\partial e_s}{\partial I_2} (\overline{\mathbb{B}}^{-1})^* \right) \quad (10) $$

One may consult [40] for details. The pressure $\rho$ is identified with the spherical part of $\mathcal{S}$, that is $\rho := -\rho J \frac{\partial e_v}{\partial J}$. So one has $\mathcal{S}^* = 2\rho \left( \frac{\partial e_s}{\partial I_1} \overline{\mathbb{B}}^* - \frac{\partial e_s}{\partial I_2} (\overline{\mathbb{B}}^{-1})^* \right)$. The hyperbolicity of system (1-5) with the closure (10) has been studied in [24, 43] for instance. It ends up with reasonable assumptions, including the usual ones for fluid, $T > 0$, $e_v^2 = \frac{\partial \rho(\eta, \rho)}{\partial \rho} > 0$.

The fluid regime is obtained by taking $e_s = 0$. In this case, $\mathcal{S} = \rho J \frac{\partial e_v}{\partial J}|_{\mathbb{F}}$, and only the balance of $\mathbf{J}$ is required. It can be deduced from (5)

$$ \frac{d}{dt} \mathbf{J} = \mathbf{J} \nabla \cdot \mathbf{v} \quad (11) $$

**Remark 2.1** (Example of expressions for $e_v$ and $e_s$). The simplest expression for $e_v$ is derived from the well-known perfect gas law. In this case, the specific internal energy is given by:

$$ \epsilon(\rho, \eta) = \rho^{\gamma-1} \exp\left(\frac{\eta}{C_v}\right), $$

where $C_v$ is the specific heat capacity at constant volume, $\rho$ is the density, $\eta$ is the specific entropy, and $\gamma$ is the adiabatic index.

For a simple example of $e_s$, we refer to Appendix A, which discusses Neo-Hookean materials.

For the sake of consistency with the closure of $\mathcal{S}$, we will consider an evolution equation of $\overline{\mathbb{B}}$ instead of $\mathbb{F}$:

$$ \frac{d}{dt} \overline{\mathbb{B}} = (\nabla \mathbf{v})^* \overline{\mathbb{B}} + \overline{\mathbb{B}}(\nabla \mathbf{v})^{*,T} \quad (12) $$

This evolution equation is established in A. Let us emphasize that $\det(\overline{\mathbb{B}}) = 1$. To contain the same amount of information than in Equation (5), this equation should be supplemented by an evolution equation for $\mathbf{J}$. However, this information is already contained into the first equation of System (1).

The model (1-12) with the closure (10) will be used for both fluid and solid parts of the domain in the following.

### 3 | Finite Volume Updated-Lagrangian Discretization

In this section, one details the discretization of the model in an updated-Lagrangian context. This part takes inspiration in [18, 19, 22–24, 44] in which discrete schemes of this type are proposed. These are extensions to hyperelasticity of the Eulychyd [25, 45] and Glace [20, 21] schemes for the gas dynamics. To induce a finite volume structure, we integrate equations (1-5) on any open subset $\omega(t) \subset \Omega$ moving at the material velocity $\mathbf{v}$. We make use of the Reynold's transport theorem and Green's formula (see for instance [46]). Owing to suitable regularity conditions on the unknowns, we obtain the following integrated form of the equations

$$ \forall \omega(t) \subset \Omega, \quad \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} \int_{\omega(t)} \rho = 0, \\ \displaystyle \frac{d}{dt} \int_{\omega(t)} \rho \mathbf{v} = \int_{\partial\omega(t)} \mathcal{S}\mathbf{n}, \\ \displaystyle \frac{d}{dt} \int_{\omega(t)} \rho E = \int_{\partial\omega(t)} \mathcal{S}\mathbf{v} \cdot \mathbf{n}, \end{array} \right. \quad (13) $$

where $\partial\omega(t)$ is the boundary of $\omega(t)$ and $\mathbf{n}$ is the outward unit normal to $\partial\omega(t)$ with respect to $\omega(t)$.

We now formally apply this form of the equations to the cells of a mesh. Let us consider a regular mesh $\mathcal{M}$ of $\Omega$. We consider here
that $\mathcal{M}$ is defined by a finite collection of cells $j$ that partition $\Omega$.
Specifically, we have

$$1. \forall j \in \mathcal{M}, \quad j \subset \Omega,$$

$$2. \forall \mathbf{x} \in \Omega, \quad \exists j \in \mathcal{M} \text{ s. t. } \mathbf{x} \in \bar{j},$$

$$3. \forall j, k \in \mathcal{M}, \quad \bar{j} \cap \bar{k} \neq \emptyset \iff \left\{ \begin{array}{l} j = k, \\ \text{or} \\ \partial j \cap \partial k \neq \emptyset \text{ and } j \cap k = \emptyset. \end{array} \right.$$

The last relationship indicates that if two distinct cells intersect, it can only be through a piece of their boundary (possibly a vertex, an edge or a face). In the remaining of the article, we use the following notations.

* $\mathcal{J}$ is the set of all cells of the mesh.

* $\mathcal{R}$ is the set of all nodes of the mesh.

* $r \in \mathcal{R}$ denotes a node.

* $V_j$ is the volume of cell $j$.

* For any function $\varphi \in L^1(\Omega)$, we define $\varphi_j := \frac{1}{V_j} \int_j \varphi$ the average of $\varphi$ over the cell $j$.

* By convention, a calligraphic capital letter with a subscript denotes the subset of the set corresponding to the calligraphic capital letter, associated with the subscript. For example $\mathcal{R}_j$ is the subset of elements (nodes) of $\mathcal{R}$ which belong to the cell $j$.

In the updated-Lagrangian framework, we consider the cells as open subsets of $\Omega$ moving at the material velocity $\mathbf{v}$, and apply (13)

$$\forall j \in \mathcal{J}, \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} \int_j \rho = 0, \\[2ex] \displaystyle \frac{d}{dt} \int_j \rho \mathbf{v} = \int_{\partial j} \mathbb{S} \mathbf{n}, \\[2ex] \displaystyle \frac{d}{dt} \int_j \rho E = \int_{\partial j} \mathbb{S} \mathbf{v} \cdot \mathbf{n}, \end{array} \right. \quad (14)$$

where $\partial j$ is the boundary of the cell $j$ and $\mathbf{n}$ is the outward normal to $\partial j$ with respect to $j$. The integrated deformation gradient balance should be obtained either in integrating (5) or in computing the deformation of the cell $j$ at time $t$ with respect to its initial configuration. However, following [18, 19, 23], we rather discretize the balance of $\mathbb{B}$:

$$\frac{d}{dt} \int_j \rho \mathbb{B} = \int_j \rho (\nabla \mathbf{v} \mathbb{B} + \mathbb{B}(\nabla \mathbf{v})^T) \quad (15)$$

We note moreover that, multiplying (11) by $\rho$ and integrating over the cell $j$ gives

$$\frac{d}{dt} V_j = \frac{d}{dt} \int_j 1 = \int_j \nabla \cdot \mathbf{v} \quad (16)$$

Strictly speaking, these integral forms require assumptions of regularity, which are not discussed here (refer for instance to [47] for a general discussion about this subject, or [44] for a specific discussion about the hyperelastic model).

## 3.1 | Discrete Gradient and Divergence Operators

We are now going to define consistent approximations of differential operators based on the geometry of the cells. We assume that a cell $j$ is uniquely defined by the position of vertices $\{\mathbf{x}_r\}_{r \in \mathcal{R}_j}$. It follows that the volume $V_j$ of cell $j$ is also a function of the vertices:

$$V_j = V_j(\{\mathbf{x}_r\}_{r \in \mathcal{R}_j}).$$

Using the chain rule, the variation of $V_j$ is related to the variation of the positions of the vertices by

$$\frac{dV_j}{dt} = \sum_{r \in \mathcal{R}_j} \nabla_{\mathbf{x}_r} V_j \cdot \frac{d\mathbf{x}_r}{dt}.$$

In the following, we define $\mathbf{C}_{jr} := \nabla_{\mathbf{x}_r} V_j$. Obviously, the volume is a $d$-order homogeneous function of the positions of the vertices, that is, $V_j(\{\lambda \mathbf{x}_r\}_{r \in \mathcal{R}_j}) = \lambda^d V_j(\{\mathbf{x}_r\}_{r \in \mathcal{R}_j})$. It follows by the Euler homogeneous function theorem that

$$V_j = \frac{1}{d} \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{x}_r \quad (17)$$

Moreover, $V_j$ is a first-order homogeneous function of each coordinate (the other coordinates being fixed). This yields the more general formula

$$V_j \mathbb{1}_d = \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \otimes \mathbf{x}_r \quad (18)$$

Notice that, because of the translation invariance of the volume

$$\sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} = \mathbf{0} \quad (19)$$

For a given internal node $r$, remarking that the volume $\sum_{j \in \mathcal{J}} V_j$ (with fixed boundary) does not depend on $\mathbf{x}_r$, one has

$$\sum_{j \in \mathcal{J}_r} \mathbf{C}_{jr} = \mathbf{0} \quad (20)$$

Note that since $V_j = \int_j 1 = \frac{1}{d} \int_j \nabla \cdot \mathbf{x}$, we have

$$\int_j \nabla \cdot \mathbf{a} = \sum_{r \in \mathcal{R}_j} \mathbf{a}(\mathbf{x}_r) \cdot \mathbf{C}_{jr} \quad (21)$$

for any affine vector-valued function $\mathbf{a}(\mathbf{x})$. By duality arguments, we also have

$$\int_j \nabla \mathbf{a} = \sum_{r \in \mathcal{R}_j} \mathbf{a}(\mathbf{x}_r) \mathbf{C}_{jr} \quad (22)$$

for any affine scalar-valued function $\mathbf{a}(\mathbf{x})$. We will then use the following approximated operators:

$$\int_j \nabla \mathbf{a} \approx \sum_{r \in \mathcal{R}_j} \mathbf{a}(\mathbf{x}_r) \mathbf{C}_{jr}, \\ \int_j \nabla \cdot \mathbf{a} \approx \sum_{r \in \mathcal{R}_j} \mathbf{a}(\mathbf{x}_r) \cdot \mathbf{C}_{jr},$$
$$
\begin{align}
\int_j \nabla \mathbf{a} &\approx \sum_{r \in \mathcal{R}_j} \mathbf{a}(\mathbf{x}_r) \otimes \mathbf{C}_{jr}, \nonumber \\
\int_j \nabla \cdot \mathcal{A} &\approx \sum_{r \in \mathcal{R}_j} \mathcal{A}(\mathbf{x}_r) \mathbf{C}_{jr}. \tag{23}
\end{align}
$$

It is shown in [48], that these formulas lead to a consistent discretization of the operators for any function in BV($\Omega$).

## 3.2 | Spatial Discretization

Using approximations (23) leads to the following semi-discrete scheme,

$$
\forall j \in \mathcal{J}, \left\{
\begin{array}{l}
\displaystyle \frac{d}{dt} m_j(t) = 0, \\
\displaystyle \frac{d}{dt} V_j(t) = \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{v}_r, \\
m_j \frac{d}{dt} \mathbf{v}_j(t) = - \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}, \\
m_j \frac{d}{dt} E_j(t) = - \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr} \cdot \mathbf{v}_r, \\
m_j \frac{d}{dt} \mathbb{B}_j(t) = \rho_j \left( \sum_{r \in \mathcal{R}_j} \mathbf{v}_r \otimes \mathbf{C}_{jr} \mathbb{B}_j + \mathbb{B}_j \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \otimes \mathbf{v}_r \right),
\end{array}
\right.
\tag{24}
$$

where we defined $m_j := \int_j \rho$.

**Remark 3.1 (Redundancy of the evolution of J).** We can deduce from the second equation a balance equation for $J_j(t) = \frac{\rho_j(0)}{\rho_j(t)}$:

$$m_j \frac{d}{dt} J_j(t) = \rho_j J_j \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{v}_r \quad (25)$$

The same equation is obtained in deducing $\frac{d}{dt} J_j(t)$ from the equation on $\mathbb{B}_j$. Moreover, using the last equation of system (24) and (25) gives us the following semi-discrete balance for $\overline{\mathbb{B}}_j$

$$m_j \frac{d}{dt} \overline{\mathbb{B}}_j(t) = \rho_j \left( \sum_{r \in \mathcal{R}_j} \mathbf{v}_r \otimes \mathbf{C}_{jr} \overline{\mathbb{B}}_j + \overline{\mathbb{B}}_j \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \otimes \mathbf{v}_r - \frac{2}{d} \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{v}_r \overline{\mathbb{B}}_j \right) \quad (26)$$

which is a natural discretization of the continuous equation for $\overline{\mathbb{B}}_j$ (see Appendix A). It shows that at this semi-discrete level, using (26) and the volume equation in (24) is completely equivalent to deducing the evolution of $J_j$ from the equation on $\mathbb{B}_j$.

The fluxes $\mathbf{F}_{jr}$ and $\mathbf{v}_r$ are determined by the linearized Godunov solver

$$\forall r \in \mathcal{R}, \left\{ \begin{array}{l} \displaystyle \sum_{j \in \mathcal{J}_r} \mathbf{F}_{jr} := \mathcal{A}_{jr}(\mathbf{v}_j - \mathbf{v}_r) - \mathcal{S}_j \mathbf{C}_{jr}, \\ \displaystyle \sum_{j \in \mathcal{J}_r} \mathbf{F}_{jr} = \mathbf{0}. \end{array} \right. \quad (27)$$

The matrices $\mathcal{A}_{jr}$ are non-negative symmetric matrices such that $\mathcal{A}_r = \sum_{j \in \mathcal{J}_r} \mathcal{A}_{jr}$ is definite, for the solver to be well-defined, and

for the scheme to be compliant with a discrete Clausius-Duhem inequality. We observe that, summing the first equation of (27) over $j \in \mathcal{J}_r$ cancels the $\mathbf{F}_{jr}$ terms and leads to a problem on $\mathbf{v}_r$ solely:

$$\mathcal{A}_r \mathbf{v}_r = \sum_{j \in \mathcal{J}_r} \mathcal{A}_{jr} \mathbf{v}_j - \sum_{j \in \mathcal{J}_r} \mathcal{S}_j \mathbf{C}_{jr} \quad (28)$$

Since we have assumed that $\mathcal{A}_r$ is invertible, $\mathbf{v}_r$ is uniquely defined and $\mathbf{F}_{jr}$ are then deduced using the first equation of (27). Multiple choices are possible for $\mathcal{A}_{jr}$. However, following [19, 23, 24] we use

$$\mathcal{A}_{jr} := \rho_j \left[ a_j^L \frac{\mathbf{C}_{jr} \otimes \mathbf{C}_{jr}}{\|\mathbf{C}_{jr}\|} + a_j^T \frac{\|\mathbf{C}_{jr}\|^2 \mathbb{I}_d - \mathbf{C}_{jr} \otimes \mathbf{C}_{jr}}{\|\mathbf{C}_{jr}\|} \right] \quad (29)$$

where $a_j^L$ and $a_j^T$ are the speed of the longitudinal and transverse waves (eigenvalues of the acoustic tensor [40]) given by

$$a_j^L = \sqrt{c_j^2 + \frac{4\mu_j}{3\rho_j}}, \quad a_j^T = \sqrt{\frac{\mu_j}{\rho_j}} \quad (30)$$

in which $c_j^2 = \frac{\partial \rho(\rho_j, \eta_j)}{\partial \rho}$ ($\rho_j, \eta_j$) and $\mu_j = \frac{E_j}{2(1+\nu_j)}$ is the mean shear modulus in the cell j with $E_j$ the Young's Modulus and $\nu_j$ the Poisson coefficient (see A). This numerical method requires to be able to compute $\mathcal{S}_j(\mathbb{B}_j, e_j)$, where $e_j := E_j - \frac{|\mathbf{v}_j|^2}{2}$.

**Remark 3.2 (Eucclhyd scheme).** The scheme described above is known as Glace [20, 24]. In the test problems we also use Eucclhyd [23, 25]. Eucclhyd design is very close to Glace. In 2D, the main difference is that $\mathbf{C}_{jr}$ is split into two vectors, $\mathbf{C}_{jr} = \mathbf{N}_{jr-} + \mathbf{N}_{jr+}$ corresponding to the two faces of the cell j which share the node r. Using Eucclhyd instead of Glace changes the definition of $\mathcal{A}_{jr}$ but the global design of the scheme remains unchanged. For the sake of conciseness, we do not describe the Eucclhyd scheme in details, and we refer the interested reader to [18, 23, 25].

**Remark 3.3 (Degrees of freedom localization).** It is possible to design a fully node-centered numerical method based on the same principles (see for instance [49] in the gas dynamics context). However, in such an approach, material interfaces are carried by the nodes, which introduces additional complexity. Specifically, using a node-centered method requires sharing the degrees of freedom at the interface between materials. As a result, the coupling algorithm would need to act directly on the primary variables rather than on the fluxes. This introduces challenges in ensuring key properties of the method, particularly conservation. For this reason, we opted for a cell-centered formulation, which naturally maintains local conservation (see Section 4).

## 3.3 | Semi-Discrete Entropy Growth

We now show that the system (24) enforces a semi-discrete entropy production.

**Proposition 1 (Entropy growth).** Defining the discrete specific internal energy $e_j := E_j - \frac{1}{2}|\mathbf{v}_j|^2$, and $\rho_j := \frac{m_j}{V_j}$, the discrete entropy $\eta_j$ which satisfies a semi-discrete version of (8)
$$T_j \frac{d}{dt} \eta_j = \frac{d}{dt} e_j - \frac{1}{2\rho_j} \mathcal{S}_j \mathbb{B}_j^{-1} : \frac{d}{dt} \mathbb{B}_j \quad (31)$$

is a non-decreasing function of t.

**Proof.** Using the definition of $e_j$, (31) yields

$$m_j T_j \frac{d}{dt} \eta_j = m_j \frac{d}{dt} E_j - \mathbf{v}_j \cdot m_j \frac{d}{dt} \mathbf{v}_j - \frac{1}{2} V_j \mathcal{S}_j \mathbb{B}_j^{-1} : \frac{d}{dt} \mathbb{B}_j.$$

From the $\mathbb{B}_j$ equation of System (24), we infer

$$V_j \mathcal{S}_j \mathbb{B}_j^{-1} : \frac{d}{dt} \mathbb{B}_j = \sum_{r \in \mathcal{R}_j} \mathcal{S}_j \mathbb{B}_j^{-1} : [\mathbf{v}_r \otimes \mathbf{C}_{jr} \mathbb{B}_j + \mathbb{B}_j \mathbf{C}_{jr} \otimes \mathbf{v}_r],$$

which simplifies into

$$V_j \mathcal{S}_j \mathbb{B}_j^{-1} : \frac{d}{dt} \mathbb{B}_j = 2 \sum_{r \in \mathcal{R}_j} \mathcal{S}_j \mathbf{v}_r \cdot \mathbf{C}_{jr}.$$

From the last equation and the system (24), the time derivatives can be substituted, which yields

$$m_j T_j \frac{d}{dt} \eta_j = \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr} \cdot (\mathbf{v}_j - \mathbf{v}_r) - \mathcal{S}_j \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{v}_r.$$

Using the approximate Riemann solver (27), we infer

$$m_j T_j \frac{d}{dt} \eta_j = \sum_{r \in \mathcal{R}_j} A_{jr}(\mathbf{v}_j - \mathbf{v}_r) \cdot (\mathbf{v}_j - \mathbf{v}_r) - \sum_{r \in \mathcal{R}_j} \mathcal{S}_j \mathbf{C}_{jr} \cdot (\mathbf{v}_j - \mathbf{v}_r) \\ - \mathcal{S}_j \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{v}_r.$$

Using (20), we remark that $\sum_{r \in \mathcal{R}_j} \mathcal{S}_j \mathbf{C}_{jr} \cdot \mathbf{v}_j = 0$, and get

$$m_j T_j \frac{d}{dt} \eta_j = \sum_{r \in \mathcal{R}_j} (\mathbf{v}_j - \mathbf{v}_r)^T A_{jr} (\mathbf{v}_j - \mathbf{v}_r).$$

Since $A_{jr}$ are non-negative matrices, $m_j T_j \frac{d}{dt} \eta_j \ge 0$. □

## 3.4 | Time Discretization

For the discrete-in-time counterpart of (24), we use a forward-Euler integration:

$$\forall j \in \mathcal{J}, \quad \left\{ \begin{array}{l} V_j^{n+1} = V_j^n + \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n, \\ m_j \mathbf{v}_j^{n+1} = m_j \mathbf{v}_j^n - \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^n, \\ m_j E_j^{n+1} = m_j E_j^n - \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^n \cdot \mathbf{v}_r^n, \\ m_j \mathbb{B}_j^{n+1} = m_j \mathbb{B}_j^n + \Delta t \rho_j^n \left( \sum_{r \in \mathcal{R}_j} \mathbf{v}_r^n \otimes \mathbf{C}_{jr}^n \mathbb{B}_j^n + \mathbb{B}_j^n \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^n \otimes \mathbf{v}_r^n \right), \\ \mathbf{x}_r^{n+1} = \mathbf{x}_r^n + \Delta t \mathbf{v}_r^n. \end{array} \right. \quad (32)$$

**Remark 3.4** (Coherence between $V^{n+1}$ and $\sqrt{\det \mathbb{B}_j^{n+1}}$). Note that, at this fully discrete level, solving the equation on $\mathbb{B}_j$

does not allow to recover exactly the equation on $V_j$. It can be shown ([18, 24]) that $\sqrt{\det \mathbb{B}_j^{n+1}} = \frac{V_j^{n+1}}{V_j^0} + O(\Delta t)$. In other words, updating $\overline{\mathbb{B}}_j^{n+1}$ with

$$m_j \overline{\mathbb{B}}_j^{n+1} = m_j \overline{\mathbb{B}}_j^n + \Delta t \rho_j^n \left( \sum_{r \in \mathcal{R}_j} \mathbf{v}_r^n \otimes \mathbf{C}_{jr}^n \overline{\mathbb{B}}_j^n + \overline{\mathbb{B}}_j^n \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^n \otimes \mathbf{v}_r^n - \frac{2}{d} \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n \overline{\mathbb{B}}_j^n \right),$$

leads to $\det(\overline{\mathbb{B}}_j^{n+1}) = \det(\overline{\mathbb{B}}_j^n) + O(\Delta t^2)$, preventing $\det(\overline{\mathbb{B}}_j^n)$ to remain equal to 1. A possible fix for this issue is proposed in [50] or in [51]. For the applications we are considering, we regard this error as acceptable.

The explicit time integration is conditionally stable. We use an estimate of the CFL condition (see for instance [19, 24]), which leads to the following formula for $\Delta t$

$$\Delta t = \gamma \min_j \left( \frac{L_j}{a_j^L + a_j^T} \right) \quad (33)$$

where $L_j = \frac{V_j}{\sum_{r \in \mathcal{R}_j} \|C_{jr}\|}$, and $\gamma < 1$ is a safety coefficient. The time-step is then proportional to the inverse of the speed of the waves, which can be very large in dense and rigid bodies.

Let us now establish a few properties of the scheme.

**Proposition 2** (Conservation). The scheme (32) is conservative in volume, mass, momentum and total energy.

**Proof.** The proof is obvious and follows for instance [20]. □

**Proposition 3** (Fluid solver). If $e_s = 0$ then the scheme (32) is the first-order Glace scheme [20].

**Proof.** In this case, the strain tensor reduces to

$$\mathcal{S}_j = -p_j \mathbb{d},$$

where we assume $p_j(e_j, \rho_j)$. Since the knowledge of $\overline{\mathbb{B}}_j$ is no more required, the system (32) reduces to

$$\forall j \in \mathcal{J} \left\{ \begin{array}{l} V_j^{n+1} = V_j^n + \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n, \\ m_j \mathbf{v}_j^{n+1} = m_j \mathbf{v}_j^n - \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^n, \\ m_j E_j^{n+1} = m_j E_j^n - \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^n \cdot \mathbf{v}_r^n, \\ \mathbf{x}_r^{n+1} = \mathbf{x}_r^n + \Delta t \mathbf{v}_r^n. \end{array} \right. \quad (34)$$

The approximate solver (27) writes

$$\forall r \in \mathcal{R}_j, \quad \left\{ \begin{array}{l} \mathbf{F}_{jr} := A_{jr}(\mathbf{v}_j - \mathbf{v}_r) + p_j \mathbf{C}_{jr}, \\ \sum_{j \in \mathcal{J}_r} \mathbf{F}_{jr} = \mathbf{0}. \end{array} \right. \quad (35)$$
Moreover, the speed of acoustic waves is $c = \sqrt{\frac{\partial p(\rho, \eta)}{\partial \rho}}$ and $\mu = 0$, which yields

$$A_{jr} := \rho_j c_j \frac{\mathbf{C}_{jr} \otimes \mathbf{C}_{jr}}{\|\mathbf{C}_{jr}\|} \qquad (36)$$

As for the hyperelastic solver, the hydrodynamic solver leads to the following equation on $\mathbf{v}_r$:

$$A_r \mathbf{v}_r = \sum_{j \in J_r} A_{jr} \mathbf{v}_j + \sum_{j \in J_r} p_j \mathbf{C}_{jr} \qquad (37)$$

The time-step estimate (33) writes now

$$\Delta t = \gamma \min_{j \in J} \left( \frac{L_j}{c_j} \right) \qquad (38)$$

## 3.5 | Monolithic Solver

To conclude this Section, a direct consequence of Proposition 3 is that we have described a Lagrangian finite-volume scheme that naturally treats the fluid-structure interaction by setting locally (in each cell) the strain tensor $\mathcal{S}$. In the remaining of the article, this solver will be referred as the monolithic solver. It is consistent, entropy stable and conservative. This approach is not new and in the case of Lagrangian finite-volume methods, it is a direct consequence of the seminal works [19, 23, 24, 44, 45].

As stated in the introduction the main issue of this explicit scheme is that the time-step is defined globally, for all materials. This may results into expensive calculations: the time-step may be over constrained in some regions.

To address this issue, in the following Section 4, we propose a partitioned approach that consists in rewriting the continuous in time scheme in a more flexible way. In Section 5, this more flexible framework is used to implement a time-subcycling method which preserves all the properties of the scheme at the fully discrete level. A special attention being given to the local conservation. In order to fix ideas, the time step will be different in the fluid and in the solid. This is the main objective of this article.

## 4 | Locally Conservative Partitioned Method

As we just wrote, in view of dealing with fluid-structure interaction using local time stepping, we recast the previous scheme into a partitioned framework. Even if in this article we use the same scheme to treat all components, our long term objective is to be able to use different methods or even different codes to solve each piece in a locally conservative fashion.

Thus, without any loss of generality we decompose the fluid-structure interaction problem into two parts: a solver for the solid and a solver for the fluid. Let us denote by $M^s$ the mesh of the solid domain $\Omega^s$ and by $M^f$ the mesh of the fluid domain $\Omega^f$. We consider here a global conforming mesh $M = M^s \cup M^f$. As illustrated in Figure 1, the blue part accounts for the solid mesh $M^s$. The green part accounts for the fluid mesh

<figure><img src="image_4.png" /><figcaption>FIGURE 1 | Partitioned computation domain.</figcaption></figure>

$M^f$. The coupling interface is the red line and it is defined as the intersection $\Gamma = M^s \cap M^f$ which consists in the common boundaries of these meshes.

### 4.1 | Abstract Coupling Procedure

The three next subsections make intensive use of optimization techniques, and the interested reader is invited to refer to [52], Chap. 9 and 10 for more details. Even if we will, in the end, consider the discrete algorithms presented before, we now propose an abstract framework for the coupling that can be applied to a wider scope of discretizations. Since the balance of $\mathcal{B}$ is restricted to the solid domain, it will not be considered in the coupling procedure. We denote $\mathcal{R}^s$ (resp. $\mathcal{J}^s$) the set of nodes (resp. cells) belonging to $M^s$ and $\mathcal{R}^f$ (resp. $\mathcal{J}^f$) the set of nodes (resp. cells) belonging to $M^f$. In the following, we call the number of elements of the set $A$, $n_A = \text{card}(A)$. Let us write a general form of the system (24)

$$\frac{d}{dt} \mathbf{U}^i + \mathcal{M}^i \mathbf{G}^i = \mathbf{0} \quad i \in \{s, f\} \qquad (39)$$

where we denote $\forall i \in \{s, f\}$, $\mathbf{U}^i \in \mathbb{R}^{(d+2)n_{ji}}$ the state vector defined on mesh $M^i$, $\mathbf{G}^i \in \mathbb{R}^{(d+2)n_{gi}}$ the vector of the fluxes, and $\mathcal{M}^i \in \mathbb{R}^{(d+2)n_{ji} \times (d+2)n_{gi}}$. The set $\mathcal{E}$ corresponds to the locations where the fluxes are evaluated, for instance the set of the faces or of the nodes of the mesh. We use conservative schemes for both domains, which means that $\forall k$, $\sum_j \mathcal{M}_{jk}^i = 0$, except for boundary cells. We assume moreover that the fluxes can be computed as the unique solution of

$$\mathbf{G}^i = \underset{\mathbf{H}^i}{\operatorname{argmin}} J^i(\mathbf{H}^i) \quad i \in \{s, f\} \qquad (40)$$

where $J^i: \mathbb{R}^{(d+2)n_{gi}} \to \mathbb{R}$ is an $\alpha$-convex differentiable function¹. The fact that $J^i$ is $\alpha$-convex ensures that each of the minimization problems has a unique solution (refer for instance to [52]). The abstract algorithm consists in solving at each time iteration

$$\left\{ \begin{array}{l} \mathbf{G} = (\mathbf{G}^s, \mathbf{G}^f)^T = \underset{\mathbf{H}=(\mathbf{H}^s, \mathbf{H}^f)^T}{\operatorname{argmin}} (J^s(\mathbf{H}^s) + J^f(\mathbf{H}^f)), \\ \qquad \text{s. t. } \mathbf{G}_e^s = \mathbf{G}_e^f \text{ if } e \in \Gamma. \end{array} \right. \qquad (41)$$
We denote $\mathcal{K} = \{ \mathbf{G} \in \mathbb{R}^{(d+2)(n_{\mathrm{e}}+n_{\mathrm{f}})} , \mathbf{G}_{\mathrm{e}}^{s} = \mathbf{G}_{\mathrm{e}}^{f} \text{ on } \Gamma \}$, the set of vectors satisfying the constraint.

One easily checks that the set $\mathcal{K}$ has the following properties.

**Proposition 4** (Properties of $\mathcal{K}$).

1. $\mathcal{K}$ is nonempty ($\mathbf{0}$ is in the set) and closed.

2. $\mathcal{K}$ is convex.

3. If $\mathbf{G} \in \mathcal{K}$, $\mathbf{G} + \mathbf{G}_{\mathrm{a}} \in \mathcal{K}$, where $\mathbf{G}_{\mathrm{a}} = (\mathbf{a}, \dots, \mathbf{a})^T$ is the constant vector of $\mathbb{R}^{(d+2)(n_{\mathrm{e}}+n_{\mathrm{f}})}$ of component $\mathbf{a} \in \mathbb{R}^{d+2}$.

4. $\mathcal{K}$ is a cone, that is, if $\mathbf{G} \in \mathcal{K}$, $\mu \mathbf{G} \in \mathcal{K}$, $\forall \mu \ge 0$.

The first two properties, associated with the $\alpha$-convexity of $J^s$ and $J^f$, ensure the existence and uniqueness of a minimum for $J^s + J^f$ on $\mathcal{K}$ (see [52]). The two other properties are of interest, in the following, to demonstrate the momentum and total energy conservation of the global algorithm.

*Remark 4.1* (Different sizes of the systems). Strictly speaking, systems (32) and (34) have not the same size, because of the equation on $\overline{\mathbb{B}}$ in the system (32), and then the above framework cannot be directly applied. We will see in the next section how we circumvent this difficulty, in considering the evolution of $\overline{\mathbb{B}}$ as a post-treatment.

## 4.2 | Applications to the Schemes (27–32) and (34, 35)

We have seen previously that the solvers (27) and (35) are equivalent to the resolution of local problems (28) and (37) on the nodal velocity $\mathbf{v}_r$, the fluxes $\mathbf{F}_{jr}$ being deduced from the nodal velocities using (27) and (35). We denote $\mathbf{V}^s \in \mathbb{R}^{dn_{\mathrm{R}}^s}$, $\mathbf{A}^s \in \mathbb{R}^{dn_{\mathrm{R}}^s} \times^{dn_{\mathrm{R}}^s}$ and $\mathbf{B}^s \in \mathbb{R}^{dn_{\mathrm{R}}^s}$ (resp. $\mathbf{V}^f$, $\mathbf{A}^f$ and $\mathbf{B}^f$) the vector of the nodal velocities $\mathbf{v}_r$, $r \in \mathbb{R}^s$, the bloc diagonal matrix of bloc $\mathbf{A}_{\mathrm{r}}$, $r \in \mathbb{R}^s$ and the vector of the right-hand-side of (27) for $r \in \mathbb{R}^s$ (resp. $\mathbf{v}_r$, $r \in \mathbb{R}^f$ and $\mathbf{A}_{\mathrm{r}}$, $r \in \mathbb{R}^f$ and the right-hand-side of (35) for $r \in \mathbb{R}^f$). For instance for $r \in \mathbb{R}^s$,

$$ \mathcal{A}^s = \left( \begin{array}{cccc} \mathcal{A}_1 & & & \\ & \ddots & & \\ & & \mathcal{A}_r & \\ & & & \mathcal{A}_{R^s} \end{array} \right) \quad \text{and} \quad \mathbf{B}^s = \left( \begin{array}{c} \sum_{j \in J_1} \mathcal{A}_{j1} \mathbf{v}_j - \sum_{j \in J_1} \mathbb{S}_j \mathbf{C}_{j1} \\ \vdots \\ \sum_{j \in J_r} \mathcal{A}_{jr} \mathbf{v}_j - \sum_{j \in J_r} \mathbb{S}_j \mathbf{C}_{jr} \\ \vdots \\ \sum_{j \in J_{R^s}} \mathcal{A}_{jR^s} \mathbf{v}_j - \sum_{j \in J_{R^s}} \mathbb{S}_j \mathbf{C}_{jR^s} \end{array} \right) \quad (42) $$

**Proposition 5** (Solid minimization problem). Solving (27) is equivalent to minimizing the function

$$ \left\{ \begin{array}{l} J^s : \mathbb{R}^{dn_{\mathrm{R}}^s} \rightarrow \mathbb{R}, \\ \mathbf{V} \rightarrow J^s(\mathbf{V}) = \frac{1}{2} \mathcal{A}^s \mathbf{V} \cdot \mathbf{V} - \mathbf{B}^s \cdot \mathbf{V}. \end{array} \right. \quad (43) $$

*Proof.* As a bloc diagonal matrix of positive symmetric definite matrices $\mathcal{A}_{\mathrm{r}}$, $\mathcal{A}^s$ is positive symmetric and definite. Thus, $J^s$ is

$\alpha$-convex (see for instance [52]). And so, the problem (43) has a unique solution given by $\nabla J^s = \mathbf{0}$, that is

$$ \mathcal{A}^s \mathbf{V} = \mathbf{B}^s \quad (44) $$

Since $\mathbf{V}^s$ is solution of this system, it is the unique solution. $\square$

For the same reason, the following proposition holds.

**Proposition 6** (Fluid minimization problem). Solving (35) is equivalent to minimizing the function

$$ \left\{ \begin{array}{l} J^f : \mathbb{R}^{dn_{\mathrm{R}}^f} \rightarrow \mathbb{R}, \\ \mathbf{V} \rightarrow J^f(\mathbf{V}) = \frac{1}{2} \mathcal{A}^f \mathbf{V} \cdot \mathbf{V} - \mathbf{B}^f \cdot \mathbf{V}. \end{array} \right. \quad (45) $$

Let us introduce $\mathbf{V}_{\Gamma}^s \in \mathbb{R}^{dn_{\mathrm{R}}^f}$ (resp. $\mathbf{V}_{\Gamma}^f \in \mathbb{R}^{dn_{\mathrm{R}}^f}$) the vector of the nodal velocities $\mathbf{v}_r$, $r \in \mathbb{R}^s \cap \Gamma$ (resp. $\mathbf{v}_r$, $r \in \mathbb{R}^f \cap \Gamma$). Our coupling problem can be formulated as, solve

$$ \mathbf{V} = (\mathbf{V}^s, \mathbf{V}^f)^T = \underset{\mathbf{W}=(\mathbf{W}^r, \mathbf{W}^f)^T \in \mathcal{K}}{\operatorname{argmin}} (J^s(\mathbf{W}^s) + J^f(\mathbf{W}^f)) \quad (46) $$

where $\mathcal{K} = \{ \mathbf{V} \in \mathbb{R}^{d(n_{\mathrm{R}}^s + n_{\mathrm{f}})} , \mathbf{V}_{\Gamma}^s = \mathbf{V}_{\Gamma}^f \}$. One easily checks that $\mathcal{K}$ satisfies the conditions of Proposition 4. An important remark is that solving problem (46) can be reformulated in minimizing

$$ \left\{ \begin{array}{l} J : \mathbb{R}^{d(n_{\mathrm{R}}^s + n_{\mathrm{f}})} \rightarrow \mathbb{R}, \\ \mathbf{V} \rightarrow J(\mathbf{V}) = \frac{1}{2} \mathcal{A} \mathbf{V} \cdot \mathbf{V} - \mathbf{B} \cdot \mathbf{V}, \end{array} \right. \quad (47) $$

with

$$ \mathcal{A} = \begin{pmatrix} \mathcal{A}^s & \mathbf{O} \\ \mathbf{O} & \mathcal{A}^f \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} \mathbf{B}^s \\ \mathbf{B}^f \end{pmatrix} \quad (48) $$

on $\mathcal{K}$. We now prove the conservation of momentum and total energy using the fluxes solution of the problem (47).

**Proposition 7** (Volume conservation). The coupled problem which uses the discretization (32) in $\mathcal{M}^s$, (34) in $\mathcal{M}^f$ with the fluxes obtained in solving the minimization problem (46), is conservative in volume, that is:

$$ \sum_{j \in J^s} V_j^{n+1} + \sum_{j \in J^f} V_j^{n+1} = \sum_{j \in J^s} V_j^n + \sum_{j \in J^f} V_j^n \quad (49) $$

*Proof.* This is a direct consequence of the definition of $\mathcal{K}$ and property (20). $\square$

**Proposition 8** (Momentum conservation). The coupled problem which uses the discretization (32) in $\mathcal{M}^s$, (34) in $\mathcal{M}^f$ with the fluxes obtained in solving the minimization problem (46), is conservative in momentum, that is:

$$ \sum_{j \in J^s} m_j \mathbf{V}_j^{n+1} + \sum_{j \in J^f} m_j \mathbf{V}_j^{n+1} = \sum_{j \in J^s} m_j \mathbf{V}_j^n + \sum_{j \in J^f} m_j \mathbf{V}_j^n \quad (50) $$
Proof. Let $\mathbf{V} \in \mathcal{K}$, be solution of problem (46). We recall the Euler inequality:

$$ \nabla J(\mathbf{V}) \cdot (\mathbf{V} - \mathbf{U}) \geq 0, \quad \forall \mathbf{U} \in \mathcal{K} \qquad (51) $$

Since $\mathcal{K}$ is translational invariant, $(\mathbf{V} + \mathbf{V}_a) \in \mathcal{K}$, $\forall \mathbf{a} \in \mathbb{R}^d$, where $\mathbb{R}^{d(n_R, n_{R_f})} \ni \mathbf{V}_a = (\mathbf{a}, \dots, \mathbf{a})^T$. Using $\mathbf{U} = \mathbf{V} + \mathbf{V}_a$ and $\mathbf{U} = \mathbf{V} + \mathbf{V}_{-a}$ in (51) yields

$$ \nabla J(\mathbf{V}) \cdot \mathbf{V}_a = 0, \quad \forall \mathbf{a} \in \mathbb{R}^d \qquad (52) $$

Expanding this expression and using first lines of the Riemann solvers (27) and (35) yields

$$ \forall \mathbf{a} \in \mathbb{R}^d, \quad \left( \sum_{r \in R^s} \sum_{j \in J_r} \mathbf{A}_{jr}(\mathbf{v}_r - \mathbf{v}_j) - \mathbf{S}_j \mathbf{C}_{jr} + \sum_{r \in R^f} \sum_{j \in J_r} \mathbf{A}_{jr}(\mathbf{v}_r - \mathbf{v}_j) - p_j \mathbf{C}_{jr} \right) \cdot \mathbf{a} \\ = - \left( \sum_{r \in R^s} \sum_{j \in J_r} \mathbf{F}_{jr} - \sum_{r \in R^f} \sum_{j \in J_r} \mathbf{F}_{jr} \right) \cdot \mathbf{a} = 0 \qquad (53) $$

which obviously implies, that $\sum_{r \in R^s} \sum_{j \in J_r} \mathbf{F}_{jr} + \sum_{r \in R^f} \sum_{j \in J_r} \mathbf{F}_{jr} = \mathbf{0}$. Inverting the summations, and using this with (32) and (34) completes the proof. $\square$

**Proposition 9** (Total energy conservation). The coupled problem which uses the scheme (32) in $\mathcal{M}^s$, (34) in $\mathcal{M}^f$ with the flux obtained in solving the minimization problem (46), is conservative in total energy.

Proof. Let $\mathbf{V} \in \mathcal{K}$, be solution of problem (46). Since, $\mathcal{K}$ is a cone (meaning that $\mathbf{U} \in \mathcal{K} \Rightarrow \mu \mathbf{U} \in \mathcal{K}$, $\forall \mu \in \mathbb{R}^+$), we can use Euler inequality (51) with $\mathbf{U} = \mathbf{0}$ and $\mathbf{U} = 2\mathbf{V}$, which yields

$$ \nabla J(\mathbf{V}) \cdot \mathbf{V} = 0 \qquad (54) $$

In expanding this expression and using first lines of the Riemann solvers (27) and (35) we infer

$$ \sum_{r \in R^s} \sum_{j \in J_r} (\mathbf{A}_{jr}(\mathbf{v}_r - \mathbf{v}_j) - \mathbf{S}_j \mathbf{C}_{jr}) \cdot \mathbf{v}_r + \sum_{r \in R^f} \sum_{j \in J_r} (\mathbf{A}_{jr}(\mathbf{v}_r - \mathbf{v}_j) - p_j \mathbf{C}_{jr}) \cdot \mathbf{v}_r \\ = - \sum_{r \in R^s} \sum_{j \in J_r} \mathbf{F}_{jr} \cdot \mathbf{v}_r - \sum_{r \in R^f} \sum_{j \in J_r} \mathbf{F}_{jr} \cdot \mathbf{v}_r = 0. \qquad (55) $$

Inverting the summations, and using this with (32) and (34) completes the proof. $\square$

## 4.3 | Saddle Point Problem Formulation

Let us denote $rs \in \Gamma^s$ and $rf \in \Gamma^f$, two coinciding vertices of $\mathcal{M}^s$ and $\mathcal{M}^f$. The constraints $\mathbf{K}_r(\mathbf{V}) = \mathbf{v}_{rs} - \mathbf{v}_{rf} = \mathbf{0}$ are linear, thus derivable and convex. Moreover, the matrices $\frac{d}{d\mathbf{V}}\mathbf{K}_r(\mathbf{V})$ are linearly independent. Consequently, we can benefit of the framework of Lagrange multipliers to solve the minimization problem (46). Let us define $\lambda \in \mathbb{R}^{n_{R_f}}$, and the Lagrangian $\mathcal{L}(\mathbf{V}, \lambda) = J(\mathbf{V}) + \lambda \cdot \mathbf{K}(\mathbf{V})$, where $\mathbf{K}(\mathbf{V})$ is the vector of component $\mathbf{K}_r(\mathbf{V})$. Then solving (46) is equivalent to finding $(\mathbf{V}, \lambda)$, s.t.

$$ \nabla \mathcal{L} = \mathbf{0} \Leftrightarrow \left\{ \begin{array}{l} \mathbf{A}\mathbf{V} - \mathbf{B} + \mathbf{K}^T \lambda = \mathbf{0}, \\ \mathbf{K}\mathbf{V} = \mathbf{0}, \end{array} \right. \qquad (56) $$

where $\mathbf{K} \in \mathbb{R}^{n_{R_f} \times d n_{R_f}}$ has the line-vector block structure

$$ \mathbf{K} = \begin{pmatrix} \mathbf{K}_{1,1}^T & \dots & \mathbf{K}_{1,r}^T & \dots & \mathbf{K}_{1,n_R}^T \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \mathbf{K}_{k,1}^T & \dots & \mathbf{K}_{k,r}^T & \dots & \mathbf{K}_{k,n_R}^T \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \mathbf{K}_{n_{R_f},1}^T & \dots & \mathbf{K}_{n_{R_f},r}^T & \dots & \mathbf{K}_{n_{R_f},n_R}^T \end{pmatrix} \qquad (57) $$

and $\mathbf{K}_{k,r} \in \mathbb{R}^d$ is defined (when $d=3$) by

$$ \mathbf{K}_{k,r}^T = \begin{cases} (1, 1, 1) & \text{if } r \text{ is the } k\text{th vertex of } \Gamma^s, \\ (-1, -1, -1) & \text{if } r \text{ is the } k\text{th vertex of } \Gamma^f, \\ (0, 0, 0) & \text{else.} \end{cases} \qquad (58) $$

Since $\mathbf{A}$ is invertible, we can express $\mathbf{V}$ as a function of $\lambda$,

$$ \mathbf{V} = \mathbf{A}^{-1}(\mathbf{B} - \mathbf{K}^T \lambda) \qquad (59) $$

and inject it in the second equation of (56) to obtain a problem on $\lambda$ only,

$$ \mathbf{K}\mathbf{A}^{-1}\mathbf{K}^T\lambda = \mathbf{K}\mathbf{A}^{-1}\mathbf{B} \qquad (60) $$

**Remark 4.2** (Physical meaning of $\lambda$). Remembering the block diagonal structure of $\mathbf{A}$, for all vertex $r$ of $\Gamma$ the right-hand-side of (60) can be rewritten into

$$ (\mathbf{K}\mathbf{A}^{-1}\mathbf{B})_k = (\mathbf{A}_{rs})^{-1}\mathbf{B}_{rs} - (\mathbf{A}_{rf})^{-1}\mathbf{B}_{rf} \qquad (61) $$

where $k$ is the line number corresponding to the common vertex $rs \in \mathcal{M}^s$ and $rf \in \mathcal{M}^f$ in $\Gamma$. Note that $\mathbf{v}_{r,free}^s = (\mathbf{A}_{rs})^{-1}\mathbf{B}_{rs}$ and $\mathbf{v}_{r,free} = (\mathbf{A}_{rf})^{-1}\mathbf{B}_{rf}$ are the constraints-free velocities on $\Gamma$ for the solid and the fluid parts. It means that $\lambda$ is the cohesion-force enforcing the equality of the velocities.

**Algorithm 1** (Constrained problem solver). A valid algorithm to solve this problem consists in

* solving separately the two problems $\mathbf{V}_{free}^s = (\mathbf{A}_{\Gamma}^s)^{-1}\mathbf{B}^s$ and $\mathbf{V}_{free}^f = (\mathbf{A}_{\Gamma}^f)^{-1}\mathbf{B}^f$,

* for all $r \notin \Gamma$, set $\mathbf{v}_r = \mathbf{v}_{r,free}$,

* for all couple of common vertices $rs, rf \in \Gamma$, compute

$$ \lambda_k = ((\mathbf{A}_{rs})^{-1} + (\mathbf{A}_{rf})^{-1})^{-1}(\mathbf{v}_{rs,free} - \mathbf{v}_{rf,free}) \qquad (62) $$

* correct only the interface velocities using (59).

Note that $(\mathbf{A}_{rs})^{-1} + (\mathbf{A}_{rf})^{-1}$ is symmetric definite positive, thus invertible.

**Proposition 10** (Equivalence with monolithic solver). Solving (56) is equivalent to solving the monolithic problem $\mathbf{A}_m\mathbf{V}_m = \mathbf{b}_m$, where $\mathbf{A}_m \in \mathbb{R}^{d(n_R - n_{R_f}) \times d(n_R - n_{R_f})}$ is a bloc diagonal matrix of bloc $\mathbf{A}_{r,m} \in \mathbb{R}^d \times \mathbb{R}^d$ and $\mathbf{b}_m \in \mathbb{R}^{d(n_R - n_{R_f})}$ is a vector of vector-bloc component $\mathbf{b}_{r,m} \in \mathbb{R}^d$, defined by

$$ \left\{ \begin{array}{l} \mathbf{A}_{r,m} = \mathbf{A}_{rs} + \mathbf{A}_{rf}, \quad \mathbf{b}_{r,m} = \mathbf{b}_{rs} + \mathbf{b}_{rf}, \quad \text{if } r \in \Gamma, \\ \mathbf{A}_{r,m} = \mathbf{A}_r, \quad \mathbf{b}_{r,m} = \mathbf{b}_r \quad \text{else.} \end{array} \right. \qquad (63) $$
**Proof.** Remembering that $\mathbf{v}_{rs,free} = (\mathbf{A}_{rs})^{-1}\mathbf{b}_{rs}$ and $\mathbf{v}_{rf,free} = (\mathbf{A}_{rf})^{-1}\mathbf{b}_{rf}$ and using identity (62) into (59), we infer for a vertex of $\Gamma_s$

$$ \mathbf{v}_{rs} = (\mathbf{A}_{rs})^{-1} (\mathbf{b}_{rs} - ((\mathbf{A}_{rs})^{-1} + (\mathbf{A}_{rf})^{-1})^{-1} ((\mathbf{A}_{rs})^{-1}\mathbf{b}_{rs} - (\mathbf{A}_{rf})^{-1}\mathbf{b}_{rf})). $$

Multiplying all terms by $\mathbf{A}_{rf}((\mathbf{A}_{rs})^{-1} + (\mathbf{A}_{rf})^{-1})\mathbf{A}_{rs}$ yields

$$ (\mathbf{A}_{rs} + \mathbf{A}_{rf})\mathbf{v}_{rs} = \mathbf{b}_{rs} + \mathbf{b}_{rf}, $$

which is the local monolithic system for a vertex on $\Gamma$. The fact that the local systems are unchanged outside of $\Gamma$ completes the proof. □

**Remark 4.3 (No-slip constraint).** In this work, we impose a no-slip condition between the solid and the fluid domains. It is possible to weaken this condition in enforcing only the sliding (refer to [53–55] in the context of two fluids). In this case, the constraints become $\mathbf{K}_r(\mathbf{V}) = (\mathbf{v}_{rs} - \mathbf{v}_{rf}) \cdot \mathbf{n}_\Gamma = 0$, where $\mathbf{n}_\Gamma$ is the normal vector to $\Gamma$. While the constraints remain linear, keeping most of the analysis unchanged, the definition of the interface $\Gamma$ is then a difficult problem, which is not the subject of this work.

## 5 | Time-Subcycling

As mentioned in the introduction, our goal is to compute the interaction between a fluid defined in the domain $\Omega^f$ and a solid defined in $\Omega^s$. Since the time-step can be several orders of magnitude larger in the fluid domain, we want to be able to handle different time-steps in the two calculation domains.

However, we want to keep the main properties that we established in the previous Section when using a common time-step: the conservation of mass, momentum and total energy, and the semi-discrete entropy growth. In the following, we denote by $\Delta t^n$ the (large) time-step in the fluid domain, and by $\delta t^{n,k}$ the (small) time-step in the solid domain as depicted in Figure 2. The superscript $k \ge 1$ accounts for a subiteration index. Actually, we enforce a semi-conformal grid in time, so we define

$$ \Delta t^n := \gamma_f \min_{j \in \Omega^f} \left( \frac{L_j^n}{c_j^n} \right), $$

$$ \delta t^{n,k} := \begin{cases} \gamma_s \min_{j \in J^s} \left( \frac{L_j^{n,k}}{a_j^{l,n,k} + a_j^{T,n,k}} \right), & \text{if } \sum_{l=0}^{k-1} \delta t^{l,n} \le \Delta t^n, \\ \Delta t^n - \sum_{l=0}^{k-1} \delta t^{l,n}, & \text{else.} \end{cases} \quad (64) $$

We denote by $q$ the smallest natural integer, such that $\sum_{k=0}^{q-1} \delta t^{n,k} = \Delta t^n$. It is actually the number of subtime steps, that is, if $q = 1$ then $\delta t^{n,0} = \delta t^n = \Delta t^n$.

In the following, to simplify the notations, we omit the superscript $n$ for the time-steps ($\Delta t \equiv \Delta t^n$ and $\delta t^k \equiv \delta t^{n,k}$). Also, in the sake of simplicity, we use the same time index for both fluid and solid quantities. It means that a quantity defined at time $t^n$ in the fluid will be updated using $\Delta t$ at time $t^{n+q}$.

According to this local time-step, the discrete schemes associated to these two domains write

$$ \begin{array}{l} V_j^{n+k+1} = V_j^{n+k} + \delta t^k \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^{n+k} \cdot \mathbf{v}_r^{n+k}, \\ m_j \mathbf{v}_j^{n+k+1} = m_j \mathbf{v}_j^{n+k} - \delta t^k \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^{n+k}, \\ m_j E_j^{n+k+1} = m_j E_j^{n+k} - \delta t^k \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^{n+k} \cdot \mathbf{v}_r^{n+k}, \\ m_j \mathbf{B}_j^{n+k+1} = m_j \mathbf{B}_j^{n+k} + \delta t^k \rho_j^{n+k} \left( \sum_{r \in \mathcal{R}_j} \mathbf{v}_r^{n+k} \otimes \mathbf{C}_{jr}^{n+k} \mathbf{B}_j^{n+k} + \mathbf{B}_j^{n+k} \sum_{r \in \mathcal{R}_j} \mathbf{v}_r^{n+k} \right), \end{array} \quad (65) $$

and

$$ \forall j \in \mathcal{J}^f, \quad \left\{ \begin{array}{l} V_j^{n+q} = V_j^n + \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n, \\ m_j \mathbf{v}_j^{n+q} = m_j \mathbf{v}_j^n - \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^n, \\ m_j E_j^{n+q} = m_j E_j^n - \Delta t \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}^n \cdot \mathbf{v}_r^n. \end{array} \right. \quad (66) $$

The solid mesh $\mathcal{M}^s$ is displaced after each sub time-step by

$$ \forall k \in \{0, \dots, q-1\}, \forall r \in \mathcal{R}^s, \quad \mathbf{x}_r^{n+k+1} = \mathbf{x}_r^{n+k} + \delta t^k \mathbf{v}_r^{n+k}, $$

and the fluid mesh $\mathcal{M}^f$ is updated at once using

$$ \forall r \in \mathcal{R}^f, \quad \mathbf{x}_r^{n+q} = \mathbf{x}_r^n + \Delta t \mathbf{v}_r^n. $$

**Proposition 11 (Conservation constraint).** The global system (65) and (66) is conservative in volume, momentum and total energy over the time-step $\Delta t$, if and only if

$$ \left\{ \begin{array}{l} \Delta t \sum_{r \in \Gamma^f} \sum_{j \in J_r} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n + \sum_{k=0}^{q-1} \delta t^k \sum_{r \in \Gamma^f} \sum_{j \in J_r} \mathbf{C}_{jr}^{n+k} \cdot \mathbf{v}_r^{n+k} = 0, \\ \Delta t \sum_{r \in \Gamma^f} \sum_{j \in J_r} \mathbf{F}_{jr}^n + \sum_{k=0}^{q-1} \delta t^k \sum_{r \in \Gamma^f} \sum_{j \in J_r} \mathbf{F}_{jr}^{n+k} = \mathbf{0}, \\ \Delta t \sum_{r \in \Gamma^f} \sum_{j \in J_r} \mathbf{F}_{jr}^n \cdot \mathbf{v}_r^n + \sum_{k=0}^{q-1} \delta t^k \sum_{r \in \Gamma^f} \sum_{j \in J_r} \mathbf{F}_{jr}^{n+k} \cdot \mathbf{v}_r^{n+k} = 0. \end{array} \right. \quad (67) $$

<figure><img src="image_5.png" /><figcaption>**FIGURE 2** | Illustration of time steps.</figcaption></figure>
*Proof.* It is sufficient to sum the system (65) over all cells $j \in \mathcal{J}^s$, and all steps $k \in \{0, \dots, q-1\}$ and the system (66) over all cells $j \in \mathcal{J}^f$. Then, inverting the summation symbols and noticing that the sums over the cells cancel away from $\Gamma$ gives the result. □

Since our method is designed for conformal meshes, it is pertinent to enforce the condition at each vertex of $\Gamma$. Thus, a local version of (67) is a sufficient conservation condition

$$ \forall r \in \Gamma, \quad \left\{ \begin{array}{l} \Delta t \sum_{j \in \mathcal{J}_r^f} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n + \sum_{k=0}^{q-1} \delta t^k \sum_{j \in \mathcal{J}_r^f} \mathbf{C}_{jr}^{n+k} \cdot \mathbf{v}_r^{n+k} = 0, \\ \Delta t \sum_{j \in \mathcal{J}_r^f} \mathbf{F}_{jr}^n + \sum_{k=0}^{q-1} \delta t^k \sum_{j \in \mathcal{J}_r^f} \mathbf{F}_{jr}^{n+k} = \mathbf{0}, \quad (68) \\ \Delta t \sum_{j \in \mathcal{J}_r^f} \mathbf{F}_{jr}^n \cdot \mathbf{v}_r^n + \sum_{k=0}^{q-1} \delta t^k \sum_{j \in \mathcal{J}_r^f} \mathbf{F}_{jr}^{n+k} \cdot \mathbf{v}_r^{n+k} = 0. \end{array} \right. $$

Finally, since our time discretization is first-order accurate, we notice that a sufficient condition of conservation consists in fixing the fluxes in the solid domain at the beginning of the large time-step $\Delta t$, leading to the conditions

$$ \forall r \in \Gamma, \quad \left\{ \begin{array}{l} \displaystyle \sum_{j \in \mathcal{J}_r^f} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n + \sum_{j \in \mathcal{J}_r^s} \mathbf{C}_{jr}^n \cdot \mathbf{v}_r^n = 0, \\ \displaystyle \sum_{j \in \mathcal{J}_r^f} \mathbf{F}_{jr}^n + \sum_{j \in \mathcal{J}_r^f} \mathbf{F}_{jr}^n = \mathbf{0}, \\ \displaystyle \sum_{j \in \mathcal{J}_r^f} \mathbf{F}_{jr}^n \cdot \mathbf{v}_r^n + \sum_{j \in \mathcal{J}_r^s} \mathbf{F}_{jr}^n \cdot \mathbf{v}_r^n = 0. \end{array} \right. \quad (69) $$

This fits exactly the framework of (41), and we saw in Section 4.3 that it can be solved in our case in considering the saddle point problem (56). Thus, we have defined a fully conservative coupling method.

**Proposition 12** (Stability of the time-subcycling). *The time-subcycling strategy has no impact on the stability conditions.*

*Proof.* The first remark is that the semi-discrete schemes are unchanged. Thus, the semi-discrete entropy growth Property 1 still holds. Moreover, the time-step in the fluid domain is still computed using (38). In the solid domain, the time-step is lesser or equal to the one computed by (33). □

# 6 | Numerical Experiments

The goal of this section is to perform numerical experiments on the method explained in the previous sections. We start with a

verification test problem which is 1D and fluid only to clearly show the impact of the time-subcycling strategy on a well-known solution. We then run a 1D fluid-structure interaction Riemann problem, including a transverse velocity. Eventually, we show the potential of the time-subcycling strategy on a 3D fluid-structure interaction problem.

## 6.1 | Elementary 1D Study on the Sod Shock Tube

Let $\Omega =]0, 1[$, we consider the following initial state of the perfect gas ($\gamma = 1.4$):

$$ \begin{align*} \rho_0 &= \begin{cases} 1, & \text{if } x < 0.5, \\ 0.125, & \text{elsewhere,} \end{cases} \\ u_0 &= 0, \\ \rho_0 &= \begin{cases} 1, & \text{if } x < 0.5, \\ 0.1, & \text{elsewhere.} \end{cases} \end{align*} $$

An exact solution can be computed for this problem [56]. Our goal for this test is to assess the method and investigate potential spurious features arising from the abrupt change in time-step in a well-known Riemann problem for the Euler system. To verify our method in this test, we impose a fictitious interface within the computational domain. We conducted several experiments by placing the interface at different locations within the domain. However, for the sake of conciseness, we present only a few results with the interface positioned on the contact discontinuity (initially at $x = 0.5$). Left part of the domain ($x < 0.5$) is arbitrarily considered as $\Omega^f$ and right part as $\Omega^s$. Other configuration choices (positions of the interface or ordering of the domains) lead to the same conclusions. Because the wave speeds are very close in the two parts of the domain, we enforce manually the time-step ratio $q$, meaning that $\Delta t$ is calculated using (64), applied to the whole computation domain $\Omega^f \cup \Omega^s$, and $\delta t = \frac{\Delta t}{q}$. First, we performed a convergence study.

The order of convergence observed is the same as for a single time-step calculation. The results are reported on Table 1 and plotted in Figure 3. The asymptotic convergence order is $\approx 0.85$.

Then, for a 1000 cells mesh, we increase gradually $q$ from 1 to $2 \times 10^4$, and report the error with respect to the exact solution in Figure 4.

We observe that, the error slightly increases from $q = 1$ to $q = 10$ and then stagnates, showing that the method has almost no impact on the quality of the numerical solution.

<table><caption>**TABLE 1** | Sod shock tube problem: convergence study for $q = 2$ at $t = 0.2$. $\rho_e$ accounts for the exact value of the density.</caption><tr><td>Number of cells</td><td>100</td><td>200</td><td>400</td><td>800</td><td>1600</td><td>3200</td><td>6400</td><td>12800</td><td>25600</td></tr><tr><td>||ρ - ρ<sub>e</sub>||<sub>L<sub>1</sub></sub></td><td>0.120</td><td>0.0075</td><td>0.0046</td><td>0.0027</td><td>0.0027</td><td>0.0016</td><td>0.0009</td><td>0.0005</td><td>0.0001</td></tr><tr><td>Convergence rate</td><td>—</td><td>0.6744</td><td>0.7151</td><td>0.7503</td><td>0.7802</td><td>0.8057</td><td>0.8266</td><td>0.8419</td><td>0.8581</td></tr></table>
<figure><img src="image_6.png" /><figcaption>FIGURE 3 | Sod shock tube problem at $t = 0.2$. Convergence study for $q = 2$. $\rho_e$ accounts for the exact value of the density.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 4 | Sod shock tube problem at $t = 0.2$. Error versus the time-step ratio $q$. $\rho_e$ accounts for the exact value of the density.</figcaption></figure>

## 6.2 | 1D Gas-Copper Interaction Problem

To assess the method on a fluid structure coupled problem, we consider a 1D hyperelastic problem given in [57] (test 4.1). It is 1D in the sense that the unknowns depend only on $x$, but the transversal velocity $v_y \neq 0$. The initial configuration is depicted in Figure 5 using the SI base units. The central part of the problem ($x(m) \in [0.2 \ 1.8]$) is exactly the elastic problem proposed in [57].

A stiffened gas equation of state is used in the fluid part

$$p = \rho(\gamma^f - 1)e - \gamma^f p_{\infty}^f,$$

with $\gamma^f = 4.22$, and $p_{\infty}^f = 6.8 \ 10^8 \text{ kg.m}^{-1}\text{s}^{-2}$. The internal energy for the solid constitutive law is given by

$$e = \frac{p + \gamma^s p_{\infty}^s}{(\gamma^s - 1)\rho} + \frac{\mu}{\rho J}(\text{tr}(\overline{\mathbb{B}}) - 2),$$

which is in the separated form (9) with $e_v = \frac{p+\gamma p_{\infty}}{(\gamma-1)\rho}$ and $e_s = \frac{\mu}{\rho J}(\text{tr}(\overline{\mathbb{B}}) - 2)$. Parameters are $\gamma^s = 4.22$, $p_{\infty}^s = 3.42 \ 10^{10} \text{ kg.m}^{-1}\text{s}^{-2}$ and $\mu = 5 \ 10^{10} \text{ kg.m}^{-1}\text{s}^{-2}$. Equation (10) yields (refer to Appendix A)

$$\mathbb{S} = -p^s \mathbb{I}_d + 2 \frac{\mu}{J} \overline{\mathbb{B}}^*.$$

We compute the reference solution on a 204800 cells mesh and compare it to the reference solution of [57]. The solutions are identical. We then performed two calculations with 2000 cells and 4000 cells until $t = 10^{-4}$s. Plots of density and pressure are displayed in Figure 6, the components of $\mathbb{S}$ in Figure 7, and of the velocity in Figure 8. Our profiles are oscillations free, and sharper than those of [57], despite our scheme being formally only first-order accurate. This improvement is a key advantage of using Lagrangian schemes, as they allow the interface to remain Lagrangian. Consequently, there is no numerical mixing between the solid and the fluid, and no subgrid model is required to
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** | Gas-copper problem: initial state. The SI units are used, densities ρ are expressed in (kg.m⁻³), pressures p are expressed in (kg.m⁻¹.s⁻²) and velocities v and its components are expressed in (m.s⁻¹).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 6** | Gas-copper interaction problem at t = 10⁻⁴s. Orange line: 2000 cells, green line: 4000 cells, and dashed line: reference.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 7** | Gas-copper interaction problem at t = 10⁻⁴s. Orange line: 2000 cells, green line: 4000 cells, and dashed line: reference.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**FIGURE 8** | Gas-copper interaction problem at t = 10⁻⁴s. Orange line: 2000 cells, green line: 4000 cells, and dashed line: reference.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**FIGURE 9** | Gas-copper interaction problem at $t = 1.5 \times 10^{-4}$s. Orange line: 2000 cells, solid domain; green line: 4000 cells, solid domain; purple line: 2000 cells, fluid domain; blue line: 4000 cells, fluid domain; and dashed line: reference.</figcaption></figure>

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

<figure><img src="image_15.png" /><figcaption>**FIGURE 10** | Gas-copper interaction problem at $t = 1.5 \times 10^{-4}$s. Orange line: 2000 cells, solid domain; green line: 4000 cells, solid domain; purple line: 2000 cells, fluid domain; blue line: 4000 cells, fluid domain; and dashed line: reference.</figcaption></figure>

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

<figure><img src="image_17.png" /><figcaption>**FIGURE 11** | Gas-copper interaction problem at $t = 1.5 \times 10^{-4}$s. Orange line: 2000 cells, solid domain; green line: 4000 cells, solid domain; purple line: 2000 cells, fluid domain; blue line: 4000 cells, fluid domain; and dashed line: reference.</figcaption></figure>

account for this mixing, as such models inevitably degrade accuracy at the interface.

Then the problem is run until $t = 1.5 \times 10^{-4}$s, so that the shock waves generated by the Riemann problem in the copper interact with the fluid. Corresponding plots are displayed in Figures 9–11. The time-step is computed independently in the fluid and the solid domains using (64). We observe a mean ratio, $\frac{\Delta t}{\delta t} \approx 2.81$, during the whole calculation.

We performed a convergence study on this problem. Because we do not have the exact solution, the convergence rate is estimated with:

$$ \text{Convergence rate} = \frac{\ln(\|\rho_{4\Delta x} - \rho_{2\Delta x}\|_{L_1}) - \ln(\|\rho_{2\Delta x} - \rho_{\Delta x}\|_{L_1})}{\ln 2} \quad (70) $$

where $\rho_{\Delta x}$ accounts for the solution on a $\frac{2}{\Delta x}$ cells mesh. The results are reported in Table 2 and plotted in Figure 12. We
<table><caption>TABLE 2 | Gas-copper interaction problem: convergence.</caption><tr><td>Number of cells</td><td>1600</td><td>3200</td><td>6400</td><td>12800</td><td>25600</td><td>51200</td><td>102400</td><td>204800</td></tr><tr><td>||ρ<sub>Δx</sub> - ρ<sub>2Δx</sub>||<sub>L<sub>1</sub></sub></td><td>—</td><td>2.5945</td><td>1.3144</td><td>0.6667</td><td>0.3430</td><td>0.1819</td><td>0.0986</td><td>0.0543</td></tr><tr><td>Convergence rate</td><td>—</td><td>—</td><td>0.9811</td><td>0.9794</td><td>0.9587</td><td>0.9152</td><td>0.8830</td><td>0.8618</td></tr></table>

<figure><img src="image_18.png" /><figcaption>FIGURE 12 | Gas-copper interaction problem: convergence plot.</figcaption></figure>

observe an asymptotic convergence rate of ≈ 0.85, which is consistent with what is observed for the Sod shock tube problem (see Table 1).

## 6.3 | 3D Experiment

We now perform a 3D experiment to evaluate the efficiency of the method in this context, and report on the computation time benefit. Our benchmark consists in a square elastic plate clamped on all edges and impacted by a shock wave. From a numerical point of view, we aim at reporting on the benefits in terms of computational time that we can expect from the time-subcycling algorithm. The initial configuration is depicted in Figure 13.

To solve this problem, we used an Arbitrary-Lagrangian-Eulerian (ALE) strategy. This approach involves performing the Lagrangian step as described in this article. Next, the fluid domain grid is smoothed to prevent excessively deformed cells, and the conserved quantities $ρ$, $ρv$, $ρE$ are remapped on the new mesh. To this end, we use a fluxing algorithm, described in Appendix B. The solid domain is kept Lagrangian. The fluid is endowed with a perfect gas equation of state $p = (γ - 1)ρe$, with $γ = 1.4$. The initial state is $(ρ^{fs}, ρ^{fs}, v^{fs}) = (10^9, 1000, 0)$ in the shocked domain and $(ρ^{fr}, ρ^{fr}, v^{fr}) = (10^5, 1000, 0)$, in the domain at rest. The solid is endowed with the same constitutive law as for the previous test 6.2. The calculation is run until $t = 0.005$. The shocked fluid mesh (respectively the first part of the mesh of the fluid at rest, the solid part mesh, and the second part of the mesh of the fluid at rest) is split on $30 \times 50 \times 50$ cells (resp. $10 \times 50 \times 50$, $7 \times 50 \times 50$ and $30 \times 50 \times 50$ cells). The total number of cells is 217500, with only 17500 in the solid part, which is favorable to our method.

<figure><img src="image_19.png" /><figcaption>FIGURE 13 | 3D test problem: Initial configuration. The blue domain is a shocked gas at $p^{fs} = 10^9$ kg.m$^{-1}$.s$^{-2}$. Light-blue domains are filled with gas at rest ($p^{fr} = 10^5$ kg.m$^{-1}$.s$^{-2}$). Green domain is the solid initially at rest $S^s = 0$.</figcaption></figure>

We run three configurations using the Euclyhd scheme. The first-one, called MONO, uses the monolithic solver. The second-one, called PART-1, uses the partitioned solver (described in 1), but with a global time-stepping strategy. The third-one, called PART-2, uses the partitioned solver, and the local time-subcycling strategy described in Section 5. We expect to recover exactly the same results for MONO and PART-1, which is a sanity check.

We display in Figure 14 the copper domain at final time for the three calculations. The deformations are very similar. It is
MONO

PART-1

PART-2

<figure><img src="image_20.png" /><figcaption>**FIGURE 14** | 3D test problem: copper domain at *t* = 0.005s for the monolithic (MONO), partitioned (PART-1) and partitioned with time-subcycling (PART-2) calculations. Color map is *ρ* (kg.m⁻³).</figcaption></figure>

more obvious in Figure 15. The maximum displacement of the right boundary of the copper domain is recorded versus time. The result of the monolithic calculation and the partitioned calculation match perfectly, while there is a small difference with the time-subcycling strategy calculation. This is expected, since the calculation is not converged. In Figure 16, we display the time-step in both domains, and the number of subcycles for the time-subcycling strategy. We observe that the time-step ratio decreases when the shock interacts with the copper domain, due to smaller cells in the fluid domain. Then it grows back after the rebound of the shock. The number of subcycles is related to this effect. The calculation is performed on four processors using a shared-memory parallelism algorithm. The global computational time of the time-subcycling calculation is around 1/8th of those of the monolithic and partitioned calculations. No over-cost is observed for the partitioned calculation comparing to the monolithic one.

**Remark 6.1 (Initial conditions).** The initial conditions chosen for this problem are somewhat artificial. As enlighten by a reviewer, an initial density of $\rho^{fr} = 1$ would be more realistic for this test problem. Using $\rho^{fr} = 1$ lead to a similar conclusion regarding computational time savings. However, in that case, the fluid at rest in the left part of the domain is heavily compressed, significantly increasing the computational effort required for remapping to maintain a satisfactory mesh quality. As a result, the remapping step dominates the total computational cost, shifting the focus away from what we intended to demonstrate. With the initial conditions chosen here, the time spent in the remapping stage is negligible compared to the time spent in

<figure><img src="image_21.png" /><figcaption>**FIGURE 15** | 3D test problem: maximum displacement of the left boundary of the copper domain versus time for the monolithic (black line), partitioned (red cross), and partitioned with time-subcycling (green line) calculations.</figcaption></figure>

the Lagrangian step, making it easier to isolate and analyze the impact of the subcycling strategy.

7 | Conclusion

In this article, we studied a new partitioned method for solving compressible fluid-structure interaction problems. Given the
<figure><img src="image_22.png" /><figcaption>FIGURE 16 | 3D test problem: Left: time-step computed from the CFL in the fluid domain (red curve) and in the solid domain (black curve) versus time. Right: number of subcycles for the PART-2 calculation versus time.</figcaption></figure>

presence of shocks, we focused on the conservation properties of the coupling method. As a result, we obtained a locally conservative method that allows decoupled time steps for each part of the calculation (i.e., in the fluid and the structure).

As an initial step, we considered only first-order approximations, which are sufficient to demonstrate the capabilities of the approach in one, two, or three dimensions. We showed that the proposed method is promising: the numerical properties of the monolithic method are preserved in the partitioned version. Furthermore, all the numerical experiments we conducted showed no deterioration in the quality of the approximation, with no spurious artifacts observed. Additionally, the method's efficiency aligns with expectations, particularly in cases where the structure dictates the time step and the majority of the mesh cells belong to the fluid. In such scenarios, we observed a speedup of eight when comparing the time-subcycling algorithm with the reference approach.

There are numerous possible extensions to this work, and we outline a few of them here. A natural extension is to develop the method to achieve second-order accuracy in both space and time. Standard second-order extensions in space, such as the MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) reconstruction [58], and in time, such as Strong-Stability-Preserving Runge-Kutta (SSP-RK) schemes [59], could be used, ensuring improved accuracy. However, two key challenges must be addressed. The first is stability, particularly in designing efficient and Galilean-invariant limiters for the stress tensor, which is far from trivial. The second is adapting our subcycling strategy to a high-order time integration, requiring a careful reformulation to ensure conservation is maintained. Another interesting extension could be related to the ALE framework of the method. Handling very large motions (as rotation) or deformation of the solid is not feasible with the ALE approach used in this study. The use of remeshing techniques, such as those proposed in [60, 61], can considerably increase the scope of our method, allowing much wider movements. Another direction is to use different numerical schemes for treating the fluid and the structure, as proposed in Section 4. Along similar lines, one could consider Lagrange-Euler coupling using different

codes. Another natural extension is the implementation of a dynamic coupling interface, where regions with small and large time steps change dynamically during the calculation.

## Acknowledgments

The authors would like to thank François Hermeline for his careful proof-reading and insightful comments.

## Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Endnotes

¹An α-convex function $J$ satisfies $J(\frac{u+v}{2}) \le \frac{J(u)+J(v)}{2} - \alpha\|u-v\|^2$, for some $\alpha > 0$, see [52] definition 9.2.4.

## References

1. Y. Bazilevs, K. Takizawa, and T. Tezduyar, *Computational Fluid-Structure Interaction: Methods and Applications* (John Wiley & Sons, 2013).

2. J. Donea, S. Giuliani, and J.-P. Halleux, “An Arbitrary Lagrangian-Eulerian Finite Element Method for Transient Dynamic Fluid-Structure Interactions,” *Computer Methods in Applied Mechanics and Engineering* 33, no. 1–3 (1982): 689–723.

3. T. Hughes, W.-K. Liu, and T. Zimmermann, “Lagrangian-Eulerian Finite Element Formulation for Incompressible Viscous Flows,” *Computer Methods in Applied Mechanics and Engineering* 29, no. 3 (1981): 329–349.

4. T. Tezduyar, “Stabilized Finite Element Formulations for Incompressible Flow Computations,” *Advances in Applied Mechanics* 28 (1992): 1–44.

5. W. Wall and E. Ramm, “Fluid-Structure Interaction Based on a Stabilized (ALE) Finite Element Method,” in *Proceedings of the European Conference on Computational Mechanics* (New Trends and Applications, 1998).

6. T. Belytschko, J. Kennedy, and D. Schoeberle, “Finite Element and Finite Difference Methods for Transient Analysis of Fluid-Structure Interaction,” *Nuclear Engineering and Design* 42, no. 2 (1977): 345–354.
7. C. Farhat, M. Lesoinne, and P. Le Tallec, "Load and Motion Transfer Algorithms for Fluid/Structure Interaction Problems With Non-Matching Discrete Interfaces: Momentum and Energy Conservation, Optimal Discretization, and Application to Aeroelasticity," *Computer Methods in Applied Mechanics and Engineering* 157, no. 1–2 (1998): 95–114.

8. C. Felippa, K. Park, and C. Farhat, "Partitioned Analysis of Coupled Mechanical Systems," *Computer Methods in Applied Mechanics and Engineering* 190, no. 24–25 (2001): 3247–3270.

9. H. Matthies and J. Steindorf, "Partitioned but Strongly Coupled Iterative Schemes for Nonlinear Fluid-Structure Interaction," *Computers & Structures* 80, no. 27–30 (2002): 1991–1999.

10. T. Belytschko and J. Kennedy, "Computer Methods for Transient Analysis of Fluid-Structure Interaction," *Nuclear Engineering and Design* 47, no. 1 (1978): 183–195.

11. T. Belytschko, H.-J. Yen, and R. Mullen, "Mixed Methods for Time Integration," *Computer Methods in Applied Mechanics and Engineering* 17-18 (1979): 259–275.

12. T. Hughes and W.-K. Liu, "Implicit-Explicit Finite Elements in Transient Analysis: Stability Theory," *Journal of Applied Mechanics* 45 (1978): 375–378.

13. S. Badia, A. Martín, and J. Principe, "Partitioned Fluid-Structure Interaction Scheme With Dynamic Time Stepping and Mesh Adaptation," *Computational Mechanics* 67 (2021): 537–554.

14. L. Heltai and D. Boffi, "Numerical Analysis of a Partitioned Method for Fluid-Structure Interaction Based on a Nitsche-Type Coupling," *Journal of Computational Physics* 426 (2021): 109903.

15. M. López and C. Kees, "A Stabilized, Fully-Implicit, Strongly-Coupled Partitioned Method for Fsi Problems With Strong Added-Mass Effect," *International Journal for Numerical Methods in Engineering* 113, no. 13 (2018): 1930–1954.

16. W. Wall, M. Gee, and E. van Brummelen, "Challenges in Computational Fluid-Structure Interaction," *Acta Numerica* 27 (2018): 215–344.

17. K. Yeo and R. Li, "Implicit-Explicit Temporal Schemes for Partitioned Fluid-Structure Interaction With Sub-Cycling," *International Journal for Numerical Methods in Engineering* 123, no. 14 (2022): 6299–6320.

18. W. Boscheri, R. Loubère, and P.-H. Maire, "A 3D Cell-Centered ADER MOOD Finite Volume Method for Solving Updated Lagrangian Hyperelasticity on Unstructured Grids," *Journal of Computational Physics* 449 (2022): 110779.

19. J. Breil, G. Georges, and P.-H. Maire, "3D Cell-Centered Lagrangian Second Order Scheme for the Numerical Modeling of Hyperelasticity System," *Computers & Fluids* 207 (2020): 104523.

20. G. Carré, S. Del Pino, B. Després, and E. Labourasse, "A Cell-Centered Lagrangian Hydrodynamics Scheme in Arbitrary Dimension," *Journal of Computational Physics* 228, no. 14 (2009): 5160–5183.

21. B. Després and C. Mazeran, "Lagrangian Gas Dynamics in Two Dimensions and Lagrangian Systems," *Archive for Rational Mechanics and Analysis* 178 (2005): 327–372.

22. G. Georges, "Développement D'un Schéma Aux Volumes Finis Centré Lagrangien Pour la Résolution 3D Des Équations de l'Hydrodynamique et de l'Hyperélasticité," *Theses, Université de Bordeaux* (2016).

23. G. Georges, J. Breil, and P.-H. Maire, "A 3D Finite Volume Scheme for Solving the Updated Lagrangian Form of Hyperelasticity," *International Journal for Numerical Methods in Fluids* 84, no. 1 (2017): 41–54.

24. G. Kluth and B. Després, "Discretization of Hyperelasticity on Unstructured Mesh With a Cell-Centered Lagrangian Scheme," *Journal of Computational Physics* 229, no. 24 (2010): 9092–9118.

25. P.-H. Maire, R. Abgrall, J. Breil, and J. Ovadia, "A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems," *SIAM Journal on Scientific Computing* 29, no. 4 (2007): 1781–1824.

26. T. Di Giusto, C.-H. Lee, A. Gil, J. Bonet, C. Wood, and M. Giacomini, "A First-Order Hyperbolic Arbitrary Lagrangian Eulerian Conservation Formulation for Nonlinear Solid Dynamics in Irreversible Processes," *Journal of Computational Physics* 518 (2024): 113322.

27. P. Barton, R. Deiterding, D. Meiron, and D. Pullin, "Eulerian Adaptive Finite-Difference Method for High-Velocity Impact and Penetration Problems," *Journal of Computational Physics* 240 (2013): 76–99.

28. A. de Brauer, A. Iollo, and T. Milcent, "A Cartesian Scheme for Compressible Multimaterial Models in 3d," *Journal of Computational Physics* 313 (2016): 121–143.

29. N. Favrie, S. L. Gavrilyuk, and R. Saurel, "Solid-Fluid Diffuse Interface Model in Cases of Extreme Deformations," *Journal of Computational Physics* 228, no. 16 (2009): 6037–6077.

30. G. Miller and P. Colella, "A Conservative Three-Dimensional Eulerian Method for Coupled Solid-Fluid Shock Capturing," *Journal of Computational Physics* 183, no. 1 (2002): 26–82.

31. S. Schoch, K. Nordin-Bates, and N. Nikiforakis, "An Eulerian Algorithm for Coupled Simulations of Elastoplastic-Solids and Condensed-Phase Explosives," *Journal of Computational Physics* 252 (2013): 163–194.

32. M. Brun, A. Batti, A. Limam, and A. Gravouil, "Explicit/Implicit Multi-Time Step Co-Computations for Blast Analyses on a Reinforced Concrete Frame Structure," *Finite Elements in Analysis and Design* 52 (2012): 41–59.

33. T. Chantrait, J. Rannou, and A. Gravouil, "Low Intrusive Coupling of Implicit and Explicit Time Integration Schemes for Structural Dynamics: Application to Low Energy Impacts on Composite Structures," *Finite Elements in Analysis and Design* 86 (2014): 23–33.

34. V. Faucher and A. Combescure, "A Time and Space Mortar Method for Coupling Linear Modal Subdomains and Non-Linear Subdomains in Explicit Structural Dynamics," *Computer Methods in Applied Mechanics and Engineering* 192, no. 5 (2003): 509–533.

35. F. E. Fekak, M. Brun, A. Gravouil, and B. Depale, "A New Heterogeneous Asynchronous Explicit-Implicit Time Integrator for Nonsmooth Dynamics," *Computational Mechanics* 60 (2017): 1–21.

36. A. Gravouil and A. Combescure, "Multi-Time-Step Explicit-Implicit Method for Non-Linear Structural Dynamics," *International Journal for Numerical Methods in Engineering* 50, no. 1 (2001): 199–225.

37. N. Mahjoubi, A. Gravouil, A. Combescure, and N. Greffet, "A Monolithic Energy Conserving Method to Couple Heterogeneous Time Integrators With Incompatible Time Steps in Structural Dynamics," *Computer Methods in Applied Mechanics and Engineering* 200, no. 9–12 (2011): 1069–1086.

38. A. Prakash and K. D. Hjelmstad, "A Feti-Based Multi-Time-Step Coupling Method for Newmark Schemes in Structural Dynamics," *International Journal for Numerical Methods in Engineering* 61, no. 13 (2004): 2183–2204.

39. V. Faucher, "Advanced Parallel Strategy for Strongly Coupled Fast Transient Fluid-Structure Dynamics With Dual Management of Kinematic Constraints," *Advances in Engineering Software* 67 (2014): 70–89.

40. M. Gurtin, E. Fried, and L. Anand, *The Mechanics and Thermodynamics of Continua* (Cambridge University Press, 2010).

41. S. Godunov and E. Romenskii, *Elements of Continuum Mechanics and Conservation Laws* (Springer Science & Business Media, 2003).

42. J. Bonet, A. Gil, and R. Wood, *Nonlinear Solid Mechanics for Finite Element Analysis: Statics* (Cambridge University Press, 2016).
43. S. Ndanou, N. Favrie, and S. Gavrilyuk, "Criterion of Hyperbolicity in Hyperelasticity in the Case of the Stored Energy in Separable Form," *Journal of Elasticity* 115, no. 1 (2014): 1–25.

44. G. Kluth, "Analyse Mathématique et Numérique de Systèmes Hyperélastiques et Introduction à la Plasticité," PhD Thesis, Université Pierre et Marie Curie (2008).

45. G. Georges, J. Breil, and P.-H. Maire, "A 3D GCL Compatible Cell-Centered Lagrangian Scheme for Solving Gas Dynamics Equations," *Journal of Computational Physics* 305 (2016): 921–941.

46. G. Duvaut, *Mécanique Des Milieux Continus* (Dunod, 2001).

47. E. Godlewski and P.-A. Raviart, *Numerical Approximation of Hyperbolic Systems of Conservation Laws*, vol. 118 (Springer Science & Business Media, 2013).

48. B. Després, "Weak Consistency of the Cell-Centered Lagrangian GLACE Scheme on General Meshes in Any Dimension," *Computer Methods in Applied Mechanics and Engineering* 199, no. 41 (2010): 2669–2679.

49. N. R. Morgan, K. N. Lipnikov, D. E. Burton, and M. A. Kenamond, "A Lagrangian Staggered Grid Godunov-Like Approach for Hydrodynamics," *Journal of Computational Physics* 259 (2014): 568–597.

50. W. Boscheri, R. Loubère, J.-P. Braeunig, and P.-H. Maire, "A Geometrically and Thermodynamically Compatible Finite Volume Scheme for Continuum Mechanics on Unstructured Polygonal Meshes," *Journal of Computational Physics* 507 (2024): 112957.

51. G. Miller and P. Colella, "A High-Order Eulerian Godunov Method for Elastic–Plastic Flow in Solids," *Journal of Computational Physics* 167, no. 1 (2001): 131–176.

52. G. Allaire, *Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and Numerical Simulation* (OUP Oxford, 2007).

53. S. Bertoluzza, S. Del Pino, and E. Labourasse, "A Conservative Slide Method for Cell-Centered Semi-Lagrangian and ALE Schemes in 2D," *ESAIM. Mathematical Modelling and Numerical Analysis* 50, no. 1 (2016): 187–214.

54. G. Clair, B. Després, and E. Labourasse, "A New Method to Introduce Constraints in Cell-Centered Lagrangian Schemes," *Computer Methods in Applied Mechanics and Engineering* 261–262 (2013): 56–65.

55. G. Clair, B. Després, and E. Labourasse, "A One-Mesh Method for the Cell-Centered Discretization of Sliding," *Computer Methods in Applied Mechanics and Engineering* 269 (2014): 315–333.

56. G. Sod, "A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws," *Journal of Computational Physics* 27, no. 1 (1978): 1–31.

57. A. Thomann, A. Iollo, and G. Puppo, "Implicit Relaxed All Mach Number Schemes for Gases and Compressible Materials," *SIAM Journal on Scientific Computing* 45, no. 5 (2023): A2632–A2656.

58. B. van Leer, "Towards the Ultimate Conservative Difference Scheme. v. A Second-Order Sequel to Godunov's Method," *Journal of Computational Physics* 32, no. 1 (1979): 101–136.

59. S. Gottlieb, C.-W. Shu, and E. Tadmor, "Strong Stability-Preserving High-Order Time Discretization Methods," *SIAM Review* 43, no. 1 (2001): 89–112.

60. S. Del Pino, "Metric-Based Mesh Adaptation for 2d Lagrangian Compressible Flows," *Journal of Computational Physics* 230, no. 5 (2011): 1793–1821.

61. S. Del Pino and I. Marmajou, "Triangular Metric-Based Mesh Adaptation for Compressible Multi-Material Flows in Semi-Lagrangian Coordinates," *Journal of Computational Physics* 478 (2023): 111975.

62. A. M. Winslow, "Numerical Solution of the Quasilinear Poisson Equation in a Nonuniform Triangle Mesh," *Journal of Computational Physics* 1, no. 2 (1966): 149–172.

63. J. M. Escobar, E. Rodríguez, R. Montenegro, G. Montero, and J. M. González-Yuste, "Simultaneous Untangling and Smoothing of Tetrahedral Meshes," *Computer Methods in Applied Mechanics and Engineering* 192, no. 25 (2003): 2775–2787.

64. P. M. Knupp, "Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities. Part I—A Framework for Surface Mesh Optimization," *International Journal for Numerical Methods in Engineering* 48, no. 3 (2000): 401–420.

65. D. Benson, "Computational Methods in Lagrangian and Eulerian Hydrocodes," *Computer Methods in Applied Mechanics and Engineering* 99, no. 2–3 (1992): 235–394.

66. P. Hoch and E. Labourasse, "A Frame Invariant and Maximum Principle Enforcing Second-Order Extension for Cell-Centered ALE Schemes Based on Local Convex Hull Preservation," *International Journal for Numerical Methods in Fluids* 76 (2014): 1043–1063.

## Appendix A

### Tensor Calculus

For ease of reading, we provide in this appendix some tensor formulas that can be found in any good continuum mechanics textbook, such as [40].

### Some Tensor Relationships

Let $\mathbb{A}$ be a non-singular matrix, then

$$d\varphi(\mathbb{A}) = \frac{\partial}{\partial \mathbb{A}} \varphi(\mathbb{A}) : d\mathbb{A} \qquad (A1)$$

It yields,

$$\frac{\partial |\mathbb{A}|}{\partial \mathbb{A}} = \frac{\mathbb{A}}{|\mathbb{A}|} \qquad (A2)$$

$$\frac{\partial \det \mathbb{A}}{\partial \mathbb{A}} = \det \mathbb{A} \mathbb{A}^{-T} \qquad (A3)$$

$$\frac{d}{dt}(\det \mathbb{A}) = \det \mathbb{A} \operatorname{tr}\left(\frac{d}{dt}\mathbb{A}\mathbb{A}^{-1}\right) \qquad (A4)$$

Let $\mathbb{S}$ be a symmetric definite matrix, and $I_1(\mathbb{S}) = \operatorname{tr}(\mathbb{S})$, $I_2(\mathbb{S}) = \frac{1}{2}(\operatorname{tr}(\mathbb{S})^2 - \operatorname{tr}(\mathbb{S}^2))$ and $I_3(\mathbb{S}) = \det \mathbb{S}$ be the principal invariants of $\mathbb{S}$ then

$$\frac{\partial I_1(\mathbb{S})}{\partial \mathbb{S}} = \mathbb{I}_d \qquad (A5)$$

$$\frac{\partial I_2(\mathbb{S})}{\partial \mathbb{S}} = I_1(\mathbb{S}) \mathbb{I}_d - \mathbb{S} \qquad (A6)$$

$$\frac{\partial I_3(\mathbb{S})}{\partial \mathbb{S}} = I_3(\mathbb{S}) \mathbb{S}^{-T} \qquad (A7)$$

### Neo-Hookean Constitutive Law

Let us consider a Neo-Hookean material with the strain internal energy defined by:

$$e_s = \frac{\mu}{2\rho_0}(\bar{I}_1 - 3) \qquad (A8)$$

where $\mu$ is the shear modulus, $\bar{I}_1 = \operatorname{tr}(\bar{\mathbb{B}})$, and $\bar{\mathbb{B}} = J^{-\frac{2}{3}}\mathbb{B}$ is the isochoric counterpart of $\mathbb{B}$. The $-3$ is just a cosmetic term, which ensures that if the material is at rest, $e_s = 0$. Now, let us compute $\mathbb{S} = 2\rho \frac{\det \mathbb{B}}{\partial \mathbb{B}} \mathbb{B}$. We have $\bar{I}_1 = J^{-\frac{2}{3}} I_1$. Thus using the chain rule, and the tensor identities

$$\frac{\partial e_s}{\partial \mathbb{B}} = \frac{\mu}{2\rho_0} \left( -\frac{2}{3} J^{-\frac{5}{3}} \left( \frac{1}{2} J \mathbb{B}^{-1} \right) I_1 + J^{-\frac{2}{3}} \mathbb{I}_d \right) \qquad (A9)$$
It yields

$$
\begin{align}
\frac{\partial e_s}{\partial \mathbb{B}} &= \frac{\mu}{2\rho_0} \left( -\frac{1}{3} J^{-\frac{2}{3}} I_1 \mathbb{1}_d + J^{-\frac{2}{3}} \mathbb{B} \right), \nonumber \\
&= \frac{\mu}{2\rho_0} \left( -\frac{1}{3} \bar{I}_1 \mathbb{1}_d + \bar{\mathbb{B}} \right), \tag{A10} \\
&= \frac{\mu}{2\rho_0} \bar{\mathbb{B}}^*. \nonumber
\end{align}
$$

Using $\rho J = \rho_0$, we deduce

$$
\mathbb{S} = \frac{\mu}{J} \bar{\mathbb{B}}^* \quad (A11)
$$

Derivative of $\overline{\mathbb{B}}$

We start with the balance equation of F in total derivative form (see [40])

$$
\frac{d}{dt} \mathbb{F} = \nabla \mathbf{v} \mathbb{F} \quad (A12)
$$

The chain rule gives

$$
\frac{d}{dt} \mathbb{B} = \nabla \mathbf{v} \mathbb{B} + \mathbb{B} \nabla \mathbf{v}^T \quad (A13)
$$

We deduce

$$
\frac{d}{dt} \overline{\mathbb{B}} = \frac{d}{dt} (J^{-\frac{2}{3}} \mathbb{B}) = J^{-\frac{2}{3}} \frac{d}{dt} \mathbb{B} + \mathbb{B} \frac{d}{dt} (J^{-\frac{2}{3}}) \quad (A14)
$$

Expanding $\frac{d}{dt}(J^{-\frac{2}{3}})$ yields

$$
\begin{align}
\frac{d}{dt}(J^{-\frac{2}{3}}) &= \frac{d}{dt}\left((J^2)^{-\frac{1}{3}}\right), \tag{A15} \\
&= -\frac{1}{3}J^{-\frac{2}{3}}\operatorname{tr}\left(\frac{d}{dt}\mathbb{B}\mathbb{B}^{-1}\right), \quad (\det(\mathbb{B}) = J^2). \nonumber
\end{align}
$$

Using (A13) one gets

$$
\begin{align}
\operatorname{tr}\left(\frac{d}{dt}\mathbb{B}\mathbb{B}^{-1}\right) &= \left(\nabla\mathbf{v}\mathbb{B} + \mathbb{B}\nabla\mathbf{v}^T\right) : \mathbb{B}^{-1}, \tag{A16} \\
&= \operatorname{tr}\left(\nabla\mathbf{v} + \nabla\mathbf{v}^T\right), \quad \text{using the symmetry of } \mathbb{B}. \nonumber
\end{align}
$$

Moreover,

$$
J^{-\frac{2}{3}} \frac{d}{dt} \mathbb{B} = \nabla \mathbf{v} \overline{\mathbb{B}} + \overline{\mathbb{B}} \nabla \mathbf{v}^T \quad (A17)
$$

Plugging (A16) and (A17) into (A14) yields

$$
\frac{d}{dt} \overline{\mathbb{B}} = (\nabla \mathbf{v})^* \overline{\mathbb{B}} + \overline{\mathbb{B}} (\nabla \mathbf{v})^{*,T} \quad (A18)
$$

Appendix B

Remap Step

In this appendix, we provide a short description of the remap step of
the method. Consider the space integration of System (1) over a domain
$\omega(t) \subset \Omega$ moving with an arbitrary velocity $\mathbf{w}$. Using the Reynold's trans-
port theorem and the Green's formula, we infer:

$$
\forall \omega(t) \subset \Omega, \left\{
\begin{array}{l}
\displaystyle \frac{d}{dt} \int_{\omega(t)} \rho = \int_{\partial\omega(t)} \rho(\mathbf{w} - \mathbf{v}) \cdot \mathbf{n}, \\[2ex]
\displaystyle \frac{d}{dt} \int_{\omega(t)} \rho \mathbf{v} = \int_{\partial\omega(t)} (\mathbb{S} + \rho(\mathbf{w} - \mathbf{v}) \otimes \mathbf{v}) \mathbf{n}, \\[2ex]
\displaystyle \frac{d}{dt} \int_{\omega(t)} \rho E = \int_{\partial\omega(t)} (\mathbb{S}\mathbf{v} + \rho E(\mathbf{w} - \mathbf{v})) \cdot \mathbf{n}.
\end{array}
\right.
\quad (B1)
$$

Obviously, Sytem (13) is recovered in taking **w** = **v** in Sytem (B1). To solve
this system, we consider the cells as open subsets of Ω moving at the

velocity **w**, and we split **w** = **v** + (**w** − **v**) into two parts. A first tempo-
ral integration is performed on cells moving at the fluid velocity **v**, and
is called Lagrangian step. For this step one solves (14) with the scheme
described in Section 3. At the end of this step, all physical quantities
are updated. The second step is called remap step and consist in solving
∂_t φ = 0, φ = {ρ, ρ**v**, ρE} on cells moving at the velocity **u** = **w** − **v**, which
gives

$$
\forall j \subset \mathcal{M}, \quad \frac{d}{dt} \int_j \varphi = \int_{\partial j} \varphi \mathbf{u} \cdot \mathbf{n}, \quad (B2)
$$

Since this step purely involves mesh adaptation and does not relies to the physics of the problem, it is referred to as the remap step. The velocity **u** is determined by the grid-smoother algorithm. In this work, we used a very basic smoother called Laplacian smoother [62]. More sophisticated smoothers (as described in [63, 64]) could be implemented to increase the robustness. These algorithms move the nodes of the mesh to increase the mesh quality. The principles of the method used to solve these equations, known as swept, are described in [65]. By integrating the previous system over time, and using Green’s theorem, we obtain

$$
\int_{t}^{t+\Delta t} \frac{d}{dt} \int_{j} \varphi - \int_{t}^{t+\Delta t} \int_{\partial j} \varphi \mathbf{u} \cdot \mathbf{n} = 0 \quad (B3)
$$

where $\varphi$ is 1, $\rho$, $\rho\mathbf{u}$ or $\rho E$. In this equation, the velocity $\mathbf{u}$ and the time step $\Delta t$ can be chosen arbitrarily as long as they respect for each vertex of the mesh the constraint

$$
\int_{t}^{t+\Delta t} \mathbf{u}_r = \mathbf{x}_r^R - \mathbf{x}_r^L \quad (B4)
$$

where $\mathbf{x}_r^R$ and $\mathbf{x}_r^L$ are respectively the rezoned position and the position at the end of the Lagrangian phase of the vertex $r$. Since this step transforms a polyhedral (or polygonal) mesh into another polyhedral (or polygonal) mesh, $\mathbf{u}$ can conveniently be chosen as affine on each edge of the elements and constant in time. The previous system can then be rewritten exactly

$$
V_j^R \varphi_j^R - V_j^L \varphi_j^L = - \sum_{f \in \mathcal{F}_j} \int_{t}^{t+\Delta t} \int_{f} \varphi \mathbf{u}_f \cdot \mathbf{n} \quad (B5)
$$

where $f$ refers to the faces and $f \in \mathcal{F}_j$ refers to all the edges of the cell $j$. We can see that if we apply this formula to the function $\varphi = 1$, the right-hand-side corresponds to the algebraic volume $\delta V_f$ shown in 2D in Figure B1. An exact formula for the volume conservation equation is deduced from it

$$
V_j^R - V_j^L = \sum_{f \in \mathcal{F}_j} \delta V_f \quad (B6)
$$

<figure><img src="image_23.png" /><figcaption>FIGURE B1 | Swept method.</figcaption></figure>
Similarly, the fluxes of $\varphi = \rho, \rho \mathbf{u}$ or $\rho E$

$$ \varphi_j^R V_j^R - V_j^L V_j^L = \sum_{f \in F_j} \varphi_f \delta V_f \quad (B7) $$

At order 1, we choose $\varphi_f$ according to the sign of $\delta V_f$, so as to ensure stability. The upwind scheme is written

$$ \varphi_f = \begin{cases} \varphi_j^L & \text{if } \delta V_f > 0, \\ \varphi_j^L & \text{else,} \end{cases} \quad (B8) $$

with the notations of the Figure B1.

The first-order version of this algorithm respects the maximum principle on $\rho, \mathbf{u}$ (in the sense of the convex hull [66]) and $E$. In this work, only the fluid part of the mesh is smoothed. Extending this to solids requires addressing the nontrivial problem of remapping $\mathcal{B}$.