# A penalization-based strong partitioned coupling with application to cavitation-induced damage

Lucas Ménez, Marianne Béringhier, Philippe Parnaudeau, Éric Goncalves da Silva

► To cite this version:

Lucas Ménez, Marianne Béringhier, Philippe Parnaudeau, Éric Goncalves da Silva. A penalization-based strong partitioned coupling with application to cavitation-induced damage. Journal of Fluids and Structures, 2025, 137, pp.104364. 10.1016/j.jfluidstructs.2025.104364. hal-05131580

HAL Id: hal-05131580

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

Submitted on 27 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A penalization-based strong partitioned coupling with application to cavitation-induced damage

L. Ménez*<sup>a</sup>, M. Beringhier<sup>a,b</sup>, P. Parnaudeau<sup>a</sup>, E. Goncalves<sup>a</sup>

<sup>a</sup>Pprime, UPR 3346 CNRS, ISAE-ENSMA, Université de Poitiers, 86961, Chasseneuil cedex, France

<sup>b</sup>Arts et Metiers Institute of Technology, LAMPA, HESAM Université, 49035, Angers, France

## Abstract

A novel strong partitioned coupling strategy is developed in order to address Fluid-Structure Interaction (FSI) problems. The Brinkman penalization method is adopted to model the deformable fluid-solid interface on a fixed Cartesian grid. Originally designed for single-phase flows and rigid bodies, the penalization method is extended to compressible multiphase flows and deformable walls. This numerical model is applied to the analysis of cavitation-induced damage at the microscopic scale, focusing on the shock-induced collapse of a single bubble near an elastoplastic material. We examine the effect of initial bubble-wall distance on wall pressure, material damage and permanent wall deformation (i.e., cavitation pit). This parametric study is conducted for different material yield strengths. Both pit depth and area increase rapidly as the bubble-wall distance and yield strength decrease. Whereas closer bubbles generate a deep, circular pit, more distant bubbles can produce a shallower, annular pit. The effect of FSI coupling is thoroughly analyzed across all parametric configurations. Wall deformation results in the damping of wall pressure, leading to differences in material damage between weakly and strongly coupled simulations. At the moment of impact, the damping of wall pressure is initially governed by the ratio of the acoustic impedances of the fluid and solid media. It is then further amplified locally by plasticity or, more generally, in regions of higher deformation. A small reduction in wall pressure leads to a much more significant damping in both pit depth and pit area. While the decrease in wall pressure is locally affected by material deformation, the change in pit size remains approximately constant for all configurations.

**Keywords:** Cavitation-induced damage, fluid-structure interaction, strong coupling, bubble collapse, penalization

## 1. Introduction

Cavitation-induced damage is a major issue for hydraulic and marine applications, as it reduces the lifespan of components such as propellers, turbines and pumps [1]. Beyond its negative effects, cavitation-induced damage has interesting applications in several fields. In medicine, it is used in the destruction

<footnote>*Corresponding author, menez.lucas98@gmail.com</footnote>
of kidney stones through Shock Wave Lithotripsy (SWL) [2] or tumors through High-Intensity Focused Ultrasound (HIFU) [3], and optimization of drug delivery through sonoporation [4]. Drilling and surface cleaning capabilities are improved by self-resonant cavitating jets compared to conventional jets [5]. A comprehensive understanding of the underlying physical mechanisms remains of primary interest to mitigate negative effects of cavitation-induced damage or take advantage of its effects in a controlled manner.

More specifically, cavitation erosion results from the accumulation of microscopic craters on the material's surface, commonly referred to as cavitation pits. These pits are formed due to shock waves and liquid jet impacts on the wall, which are generated by the implosion (or collapse) of bubbles near the solid surface [6, 7]. The small spatio-temporal scales and the intensity of the impacts make experimental measurements challenging. As a result, numerical simulation becomes a crucial tool for understanding the cavitation erosion phenomenon. It is particularly valuable in predicting material damage during the early stage of erosion (i.e., the incubation stage), as it enables a detailed analysis of the fundamental mechanisms of cavitation-induced damage at the bubble scale. However, this complex multiphysics problem remains a major challenge for numerical simulations, as it involves fast dynamics, shock waves, liquid-gas interfaces, small space and time scales, fluid-structure interaction and plasticity.

The collapse of a single bubble near a rigid wall has been extensively studied numerically in the literature [8, 9, 10, 11, 12, 13]. Motivated by experimental studies showing that bubble dynamics can be strongly influenced by wall dynamics [14, 15, 16], various numerical simulations of the collapse of a single bubble near an elastic wall have been proposed [17, 18, 19, 20]. However, modeling the formation of a cavitation pit requires considering an elastoplastic material. Few strong coupling strategies have recently been proposed to address this problem, starting from the past decade. As opposed to weak coupling, strong coupling takes into account the mutual influence between the two physics. To our knowledge, Hsiao et al. [21] were the first to introduce a strong FSI coupling methodology to investigate bubble collapse near an elastoplastic material. Since then, other strongly coupled simulations have been proposed, with a key challenge remaining the numerical modeling of the deformable fluid-solid interface. Coupling methodologies can be categorized into partitioned and monolithic approaches. In contrast to the monolithic approach, the resolution of each subproblem is carried out in separate solvers in a partitioned approach. The solutions of each subproblem are then coupled in both time and space through an exchange of information between the two solvers. Turangan et al. [22] solved the coupled problem monolithically using the Free-Lagrange Method (FLM). However, the scope of that study was limited to elastic-perfectly plastic solid materials. Joshi et al. [23] and Zhang et al. [24] introduced a fluid-solid coupling model based on Smoothed Particle Hydrodynamics (SPH), extended to elastoplastic materials. Among partitioned coupling approaches, Cao et al. [25] combined Finite Volume (FV) and Finite Element (FE) solvers, employing an Embedded Boundary Method (EBM) to track the fluid-solid interface in a non-body-fitted fluid mesh. Paquette et al. [26], Sarkar et al. [27], and Firly et al. [28] introduced an Arbitrary Lagrangian-Eulerian (ALE) formulation of the fluid subproblem in FV or FE
solvers.

The present methodology employs strong coupling between fluid and structure, offering enhanced fidelity compared to weak coupling. For greater flexibility and to take advantage of high-performance solvers specific to each physics, a partitioned approach was chosen over a monolithic one. The fluid subproblem is solved in SCB [29], an in-house compressible multiphase flow solver, based on an FV discretization. The structure subproblem is addressed with the finite element method in the open-source solver FEniCS [30]. The open-source coupling library preCICE [31] acts as a bridge between the two solvers. It manages data communication and ensures the spatio-temporal coupling of each subproblem's solution. The novelty of our approach lies in the use of a penalization method [32, 33] to model the deformable fluid-solid interface. Originally designed for single-phase flows and rigid bodies, the Brinkman penalization method is extended to compressible multiphase flows and deformable walls. This extension can be used in other FSI applications since it does not require prior knowledge of the solid material's constitutive law. The penalization method enables FV discretization on a fixed Cartesian grid that does not conform to the solid geometry. The computational cost is expected to be reduced compared to traditional body-fitted methods, due to faster mesh generation, elimination of remeshing as the interface moves, and the use of efficient Cartesian solvers. The penalization-based partitioned coupling is validated through the simulation of the impact of an incident normal shock wave on a deformable vertical panel.

It has been demonstrated experimentally and numerically that the initial bubble-wall distance, also known as the stand-off distance, plays a key role in the amplitude and distribution of wall pressure [34, 10]. Several experimental studies further analyzed the effect of this parameter on material damage and the resulting pit characteristics [35, 6, 7, 36, 37]. However, only a few numerical simulations have focused on the influence of the stand-off distance on pit characteristics [21, 22]. To address this gap, we present a comprehensive analysis of the influence of stand-off distance on wall pressure, material damage and pit characteristics. To extend the scope of this parametric study, the simulations are conducted for different material yield strengths. Although significant differences between weak and strong coupling have been observed [26, 27], their underlying causes and influence on the solution of the coupled problem remain to be clarified. To shed light on the role of FSI coupling, we thoroughly analyzed differences in wall pressure, material damage, and pit characteristics observed between weakly and strongly coupled simulations.

The present paper is organized as follows. Section 2 presents the governing equations for the fluid-structure problem, their numerical resolution in separate solvers, and the penalization-based strong partitioned coupling. This methodology is then applied to the two-dimensional simulation of cavitation pit formation in Section 4. Weak coupling is first considered in Section 4.2, to focus on shock-induced bubble collapse near a rigid wall and the underlying physical mechanisms. The validity of the penalization method is also verified. Strong coupling is then addressed in Section 4.3, where the effect of FSI strong coupling is introduced, followed by an investigation of the influence of stand-off distance and material yield strength
on wall pressure, material damage, and pit characteristics. Based on this cross-parametric study, insights into the effect of FSI coupling are provided. Finally, concluding remarks and perspectives are presented in Section 5.

## 2. Governing equations and numerical methods

### 2.1. Fluid solver

The two-phase flow dynamics are described by a compressible diffuse-interface formulation based on a four-equation model [38]. It is assumed that the phases evolve at the same speed and are in mechanical and thermal equilibrium. These assumptions lead to a model with one velocity, one pressure, and one temperature. Since the collapse of a single bubble near a wall induces low condensation, mass transfer is neglected. Surface tension effects are dominated by inertial effects and are therefore neglected. Viscous terms are taken into account in the equations. This formulation has already shown excellent results in the study of bubble collapse near a rigid wall [39, 29, 13]. The four-equation model consists of three conservation equations for the conservative variables of the mixture (mass, momentum, total energy) and a transport equation for the void ratio:

$$ \frac{\partial \rho_f}{\partial t} + \nabla \cdot (\rho_f \mathbf{u}) = 0, \qquad (1) $$

$$ \frac{\partial \rho_f \mathbf{u}}{\partial t} + \nabla \cdot (\rho_f \mathbf{u} \otimes \mathbf{u}) + \nabla p + \nabla \cdot \boldsymbol{\tau} = \mathbf{0}, \qquad (2) $$

$$ \frac{\partial \rho_f \mathcal{E}}{\partial t} + \nabla \cdot ((\rho_f \mathcal{E} + p)\mathbf{u}) + \nabla \cdot (\boldsymbol{\tau} \cdot \mathbf{u}) + \nabla \cdot \mathbf{Q} = 0, \qquad (3) $$

$$ \frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = K \nabla \cdot \mathbf{u}, \qquad K = \left( \frac{\rho_l c_l^2 - \rho_g c_g^2}{\frac{\rho_l c_l^2}{1-\alpha} + \frac{\rho_g c_g^2}{\alpha}} \right), \qquad (4) $$

with $\rho_f$ the mixture density of the fluid, $p$ the fluid pressure, $\mathbf{u}$ the fluid velocity vector, $\boldsymbol{\tau}$ the viscous stress tensor, $\mathcal{E} = e + \mathbf{u}^2/2$ the total energy of the fluid and $e$ its internal energy, $\mathbf{Q}$ the heat flux, $\alpha \in [0, 1]$ the void ratio that captures the topology of the phases, $\rho_l$ and $\rho_g$ the liquid and gas densities, $c_l$ and $c_g$ the speeds of sound in the liquid and gas, respectively. The term $K \nabla \cdot \mathbf{u}$ models the decrease or increase of the void ratio in the mixture during the passage of a compression or rarefaction wave. The system is closed by the convex stiffened gas equation of state (EOS) for pure phases [40]. The pressure and temperature of the mixture are then obtained from the mechanical and thermal equilibrium assumptions. The speed of sound in the mixture is given by the formulation of Wallis [41], because mass transfer is neglected. Equations (1), (2), (3) and (4) are solved using the Finite Volume Method (FVM). The convective flux is computed with a HLLC scheme and the second order in space is reached using a MUSCL-Hancock reconstruction. The
viscous flux is computed with a second-order centered scheme. Time integration of the solution is performed using explicit time-marching. More details about the present in-house solver are given in [29].

In order to solve FSI problems with strong coupling while eliminating the need for remeshing, a penalization method is adopted to model the boundary conditions at the fluid-solid interface [42]. The solid body is modeled as a porous medium with very low permeability [32, 33]. In contrast to the classical approach which requires a body-fitted mesh, the penalization method enables the use of a Cartesian grid which does not necessarily fit with the wall geometry. The boundary conditions are not directly imposed at the interface, but rather modelled by adding a local source term to the right hand side of the four-equation system. The solid body is located by a mask function $\chi$ that is equal to one in the solid and zero elsewhere. In order to model a no-slip boundary condition at the fluid-solid interface, the penalization term reads as:

$$ \mathbf{f} = \begin{pmatrix} f_{\rho} \\ \frac{\chi}{\eta}(\rho_{f}\mathbf{u}_{w} - \rho_{f}\mathbf{u}) \\ \frac{\chi}{\eta}\mathbf{u} \cdot (\rho_{f}\mathbf{u}_{w} - \rho_{f}\mathbf{u}) \\ 0 \end{pmatrix}, \quad \chi(\mathbf{x}, t) = \begin{cases} 1, & \text{if } \mathbf{x} \in \Omega_{s}(t), \\ 0, & \text{otherwise,} \end{cases} \qquad (5) $$

with $\chi$ the mask function, $\eta$ the permeability parameter, $\mathbf{u}_w$ the solid velocity and $\Omega_s$ the solid domain. The expression for the density penalization term, denoted as $f_\rho$, is detailed in the next paragraph. One can note that the source term is only activated within the penalized solid region where $\chi = 1$. The governing equations are therefore not penalized in the fluid area where $\chi = 0$. The numerical error introduced by the penalization method can be divided into two contributions: the forcing error and the discretization error. On the one hand, the forcing error is associated with the penalization term and vanishes as $\eta$ tends to zero [43]. A rigorous error estimate predicts a convergence rate of $O(\eta^{3/4})$ in the incompressible limit [33], while an error of order $O(\eta)$ is observed numerically for compressible flows, in agreement with previous studies [44]. On the other hand, the discretization error of the penalization method is related to the mask function, which approximates the actual solid geometry. In this work, the permeability parameter is set to $\eta = 10^{-10}$, such that the forcing error becomes negligible. The overall penalization error is thus dominated by the discretization error of order $O(\Delta x)$, with $\Delta x$ the mesh size, as demonstrated by Ménez et al. [42]. However, the very low value of $\eta$ makes the source term extremely stiff. For an explicit time integration of the penalization term, the time step must satisfy $\Delta t \le \eta$ to ensure numerical stability. Hence, the penalization term is integrated using an implicit Euler scheme, such that the numerical stability is preserved independently of the value of $\eta$. It results in a semi-implicit overall time integration of the system. Since the penalization term contains no spatial derivatives, the implicit time integration is easily achieved without solving a linear system. Furthermore, the void ratio equation is not penalized. Using a low-dissipation numerical scheme for the convective flux discretization (HLLC-Hancock) results in a good accuracy of the
no-slip condition modeling at the fluid-solid interface. Consequently, void ratio is not transported across the fluid-solid interface.

Two variants of the penalization method are considered, depending on the formulation of the density penalization term. The first one, referred to as the standard penalization method in this paper, does not include any specific treatment for density, i.e., $f_ρ = 0$. The second one, denoted as the impedance penalization method in the following, aims to control the acoustic impedance of the porous medium in order to improve the numerical modeling of wave reflection and transmission at the fluid-solid interface. Liu and Vasilyev [45] first introduced a porosity parameter to penalize the convective term of the mass conservation equation and, consequently, to control the acoustic impedance of the penalized solid region by imposing its porosity. However, it modifies the convective flux and may alter the numerical stability of the solution, leading to a strong restriction on the time step and a significant increase in the computation cost. These aspects were confirmed through a stability analysis and numerical experiments in a previous work [42]. In order to overcome this issue, the present method focuses in controlling the acoustic impedance by imposing a constant density to the penalized solid region. Thus, the density penalization term is easily integrated implicitly and avoids the time step restriction of the standard porosity model. In a similar way to the no-slip condition, the Dirichlet condition for density is modelled with the following penalization term:

$$f_ρ = \frac{\chi}{\eta}(\rho_w - \rho_f), \qquad (6)$$

with $ρ_w$ the imposed density. In the context of FSI, $ρ_w$ is determined such that the acoustic impedance of the penalized solid region matches that of the actual solid material. Assuming plane strain in the material and no penetration of gas into the penalized region, the imposed density is given by:

$$\rho_w = \frac{\rho_s E}{(1 - \nu_s^2) \rho_l c_l^2}. \qquad (7)$$

with $ρ_s$ the actual solid density, $E$ the Young's modulus and $ν_s$ the Poisson's ratio of the solid material, $ρ_l$ the liquid density and $c_l$ the speed of sound in the liquid. More details about the context and motivations behind the development of the impedance penalization method are available in [42].

## 2.2. Structure solver

The structural evolution is investigated by solving an elastoplastic dynamics problem with the Finite Element Method (FEM), implemented within the open-source FEniCS framework [30]. The materials under consideration are exclusively metals. The equation of momentum balance of the structure reads as:

$$\rho_s \ddot{\xi} - \nabla \cdot \sigma_s = F_{\text{ext}}, \qquad (8)$$

with $ρ_s$ the solid density, $ξ$ the displacement vector, $\ddot{\xi} = \partial^2 ξ / \partial t^2$ the acceleration vector, $σ_s$ the stress tensor and $F_{\text{ext}}$ the vector of external forces. In the case of isotropic linear elasticity, stress and strain are
related by the following constitutive law, also known as the Hooke's law:

$$ \sigma_s = \lambda_s \mathrm{tr}(\varepsilon_s) \mathbf{I}_d + 2\mu_s \varepsilon_s, \qquad (9) $$

with $\lambda_s = E\nu_s / ((1 + \nu_s)(1 - 2\nu_s))$ and $\mu_s = E/(2(1 + \nu_s))$ the Lamé coefficients of the material, $E$ the Young's modulus and $\nu_s$ the Poisson's ratio of the material, $\varepsilon_s$ the strain tensor and $\mathbf{I}_d$ the identity tensor. The subscript $s$ refers to the solid. Under the small deformation assumption, the strain tensor and displacement vector are related by the following expression:

$$ \varepsilon_s = \frac{1}{2} (\nabla\xi + \nabla\xi^T), \qquad (10) $$

The combination of the momentum equation (8), the constitutive law (9) and the strain-displacement relation (10) forms the linear elastodynamics equation. When the equivalent stress exceeds the yield stress $\sigma_0$, plasticity occurs. Hence, the previous constitutive law is no longer valid since the material behavior becomes local and nonlinear. In order to be solved numerically, equation (8) is linearized and solved in an iterative Newton-Raphson procedure until global equilibrium for each time step. The linearized system to solve for each sub-iterate $k$ is expressed as:

$$ \left( M \frac{\partial \ddot{\xi}^{(k)}}{\partial \xi} + \sigma_{\text{tang}}^{(k)} \right) \delta\xi^{(k)} = \mathbf{R}^{(k)}. \qquad (11) $$

with $M$ the mass matrix, $\mathbf{R}$ the residual of the nonlinear system, $\sigma_{\text{tang}}$ the tangent matrix and $\delta\xi$ the displacement increment. Acceleration and its derivatives are computed using the Newmark method [46]. In every sub-iteration loop, the elastoplastic constitutive law is integrated using a radial return algorithm [47]. During the global step, an elastic trial stress is first calculated by solving elastodynamics equation with the total strain increment. This is followed by the local step. If the von Mises stress does not exceed the yield stress of the material, no plasticity occurs and the final stress is set to the trial stress. Otherwise, irreversible deformation occurs and the increment of plastic strain is calculated. As plastic strain is assumed to be low in the present application, a linear isotropic hardening law is adopted. The expression of the plastic strain increment can thus be derived analytically from the yield condition. The final stress state is then projected onto the admissible stress according to the incremental plastic strain and the material constitutive law. The present two-dimensional simulations are conducted in the context of plane strain assumption. The structure subproblem is solved using quadratic Lagrange elements (i.e., P2 elements).

## 2.3. Fluid-structure coupling

When a fluid and a solid structure are in contact, the force exerted by the fluid induce stress and deformation within the structure. The movement of the structure, in turn, modifies the fluid flow. The two physical systems are coupled at the fluid-structure interface through the following kinematic and dynamic conditions:

$$ u = \dot{\xi}, \qquad (12) $$
$$ \sigma_f \cdot n = \sigma_s \cdot n. \qquad (13) $$

with $u$ the fluid velocity vector, $\dot{\xi}$ the structure velocity vector, $\sigma$ the stress tensor, and $n$ the normal vector to the fluid-solid interface. The subscripts $f$ and $s$ refer to fluid and solid, respectively. For simplicity, the normal fluid stress at the wall, $\sigma_f \cdot n$, is just denoted as $\sigma_f$ in the following, as there is no ambiguity with the fluid stress tensor.

<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic diagram of the penalization-based strong partitioned coupling SCB-FEniCS-preCICE.</figcaption></figure>

In the present strategy, the kinematic condition (12) and the dynamic condition (13) are met by a strong coupling between fluid and structure. A Dirichlet-Neumann partitioned coupling is considered: the fluid solver communicates wall stress to the solid solver, which returns displacement of the coupled interface to the fluid solver in order to update the interface position and interface velocity, as well as the mask function of the penalization method. In other words, fluid loading is applied to the solid as a Neumann boundary condition and interface velocity is applied to the fluid as a Dirichlet boundary condition. A schematic diagram of the present coupling strategy is presented in Figure 1. Data communication, space interpolation at the coupled interface and time coupling of the solution between the two solvers are performed by the coupling library preCICE [31] and the FEniCS-preCICE adapter [48]. preCICE enables the use of serial/parallel or explicit/implicit coupling schemes. Because of the very small time steps required for the bubble collapse simulation, resulting from the fast dynamics of the phenomenon, a serial explicit coupling scheme is adopted in the present work. Furthermore, fluid and structure dynamics do not operate at the same temporal scales. In order to avoid unnecessary iterations and to save some computation time, fluid subcycling is performed: within one single time window $\Delta T$, the fluid solution is advanced by several time steps, $\Delta t_f < \Delta T$, before exchanging information with the solid solver. This means that equations (1)–(4) are solved multiple times within one single time window, whereas equations (8)–(11) are solved only once, since the time step of the solid solver is set to $\Delta t_s = \Delta T$ in a subcycling configuration. Fluid and structure dynamics also evolve on
different spatial scales. Very high resolution is required to properly capture the fluid flow discontinuities such as liquid-gas interfaces and shocks. This results in a much finer grid for the fluid solver than for the solid solver and, therefore, in non-matching mesh nodes at the coupled interface. Even if the mesh size is the same between fluid and solid, the difference in the discretization methods used by both solvers leads to non-matching meshes. The coupled variables are thus exchanged at the fluid-solid interface using a linear interpolation.

The penalization method is extended to a deformable interface, whose deformation is determined by the resolution of the structure subproblem using the FEM. Since the fluid equations are discretized on a fixed, non-body-fitted Cartesian grid, the extended configuration of the penalization method includes a Lagrangian description of the coupled interface (see Figure 2a). The set of $N_p$ Lagrangian points that describes the coupled interface in the fluid solver is denoted as $\mathbf{X}_f$. It is also referred to as the Immersed Boundary (IB). An additional set of Lagrangian points, composed of the sampling points $\mathbf{X}$, is placed in the fluid at a distance of $\delta = \sqrt{\Delta x^2 + \Delta y^2}$ from the coupled interface, along the normal direction (see Figure 2a). The sampling points are used in the wall stress reconstruction.

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) The penalization method: schematic diagram of the extended configuration used near a deformable interface $\Gamma_f$ separating a fluid domain $\Omega_f$ and a solid domain $\Omega_s$. Symbols: □ penalized cells; ● nodes of the coupled interface $\mathbf{X}_f$; • Eulerian points $\mathbf{x}$ (cell centers); ● sampling points $\mathbf{X}$; (b) Schematic representation of the ray casting algorithm.</figcaption></figure>

The different steps of the coupling strategy are summarized in Figure 3, considering a single time window $\Delta T = [t^n, t^{n+1}]$. These are presented in the simplified configuration where the fluid and solid time steps are equal, i.e., $\Delta T = \Delta t_f = \Delta t_s$.

At the beginning of the time window, the interface displacement $\boldsymbol{\xi}^n$ is linearly interpolated from the solid mesh to the fluid mesh using preCICE (see step ① in Figure 3). The interpolated displacement from the solid mesh $\mathbf{X}_s$ to the fluid mesh $\mathbf{X}_f$, denoted as $\tilde{\boldsymbol{\xi}}^n$, is then employed to update the position of the IB
as follows (see step ② in Figure 3):

$$X_f^n = X_f^0 + \tilde{\xi}^n, \qquad (14)$$

with $X_f^0$ the initial position of the IB.

Once the new position of the interface is known, the mask function is updated using a ray casting algorithm [49] (see step ③ in Figure 3). This technique involves casting rays horizontally from each Eulerian point, defined as the centers of the Cartesian grid cells. If the ray intersects the fluid-solid interface an even number of times, the corresponding cell is considered to be outside the solid, and the mask function is set to zero. Conversely, if the ray intersects the interface an odd number of times, the Eulerian point lies inside the solid, and the mask function is set to one (see Figure 2b).

After updating the mask function, the interface velocity is extrapolated into each cell of the penalized solid region using inverse distance weighting (see step ④ in Figure 3):

$$u_w^n(\mathbf{x}) = \left( \sum_{k=1}^{N_p} ||\mathbf{x} - \mathbf{X}_{f,k}^n||_2^{-1} \dot{\mathbf{X}}_{f,k}^n \right) / \left( \sum_{k=1}^{N_p} ||\mathbf{x} - \mathbf{X}_{f,k}^n||_2^{-1} \right), \qquad (15)$$

with $||.||_2$ the Euclidean distance and $\dot{\mathbf{X}}_f^n$ the interface velocity at time $t = t^n$, approximated using a first-order backward Euler method. The extrapolated solid velocity, $u_w$, is subsequently utilized to compute the penalization term in equation (5).

The fluid solution $U = (\rho_f, \rho_f u, \rho_f \mathcal{E}, \alpha)^T$ is then updated by solving the penalized four-equation system with the FVM (see section 2.1 and step ⑤ in Figure 3).

Furthermore, when the interface deforms, some penalized cells can suddenly become fluid cells, and vice versa. The value of the mask function thus switches abruptly in these cells. The characteristic velocity of this transition is of the order of $\Delta x/\Delta t_f$, which can be one or two orders of magnitude higher than the interface velocity. In order to ensure that the fluid state is correctly convected near the boundary, a correction step is applied to the solution at the vicinity of the interface where the mask function has been modified (see step ⑥ in Figure 3). The correction is computed as a first-order accurate average of the neighboring cells, weighted by $\chi$ or $1 - \chi$ depending on the transition of the mask function. Indeed, second-order numerical schemes are used in the discretization of the convective and viscous fluxes. However, the penalization method reduces the accuracy of the solution to first order near the interface [42]. Therefore, it is unnecessary to adopt a more sophisticated correction method here.

Finally, the fluid stress is reconstructed at the IB (see step ⑦ in Figure 3). In the present application, wall shear stress is dominated by the very high impact pressure. Hence, the wall stress is reduced to the total wall pressure. The total pressure is first interpolated at the sampling point $\mathbf{X} = (X, Y)$ from the values at the four neighboring Eulerian points, using a bilinear interpolation. The resulting pressure is given by the following equation:

$$p_{tot}(\mathbf{X}) = a_1 X + a_2 Y + a_3 XY + a_4, \qquad (16)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: The different steps involved in the penalization-based strong partitioned coupling.</figcaption></figure>

with $p_{tot} = p + \rho_f |\mathbf{u}|^2/2$ the total pressure and $p$ the static pressure of the fluid, and $a_1, a_2, a_3, a_4$ the coefficients of the bilinear interpolation. These coefficients depend on the distance between $\mathbf{X}$ and the neighboring points $\mathbf{x}_l = (x_l, y_l)$ with $l \in [1; 4]$, as well as on the pressure values at these points. They are obtained analytically by solving the following system of equations:

$$p_{tot}(\mathbf{x}_l) = a_1 x_l + a_2 y_l + a_3 x_l y_l + a_4, \quad l \in [1; 4], \qquad (17)$$

The choice $\delta = \sqrt{\Delta x^2 + \Delta y^2}$ is optimal, as it ensures that $\mathbf{X}$ is close to $\mathbf{X}_f$, while preventing the bilinear interpolation from using Eulerian points that belong to $\Omega_s$. The wall pressure is then reconstructed
using a first-order Lagrange interpolation along the normal direction, i.e., $p_{tot}(\mathbf{X}_f) = p_{tot}(\mathbf{X})$. Numerical experiments have shown that a second-order Lagrange interpolation may lead to negative wall pressure when applied across a discontinuity induced by a shock wave propagating along the wall. This scenario is often encountered during bubble collapse near a solid wall. Once the pressure is known at the IB, the wall normal stress is obtained as:

$$ \sigma_f \cdot \mathbf{n} = -p_{tot} \mathbf{I}_d \cdot \mathbf{n}, \qquad (18) $$

with $\mathbf{I}_d$ the identity matrix and $\mathbf{n}$ the wall normal vector. The normal stress is then interpolated from the fluid mesh to the structural mesh using preCICE (see step ⑧ in Figure 3). It is applied as a Neumann boundary condition for the elastoplastic problem, solved using the FEM (see section 2.2 and step ⑨ in Figure 3).

All these steps are repeated for each time window.

# 3. Validation of the computational framework

## 3.1. Fluid solver

The fluid solver (SCB) has been extensively validated on shock-bubble interaction against both experimental and numerical data in previous work [39, 29, 13].

## 3.2. Structure solver

The structure solver (FEiCS) is validated in section 4.2.2 against numerical data obtained with the commercial solver Abaqus, for the problem of cavitation pit formation.

## 3.3. Fluid-structure coupling

The penalization-based strong partitioned coupling is validated through the simulation of an elastic panel deformation induced by the impact of an incident normal shock wave. The goal is to validate the computational framework for predicting wall deformation under shock wave impact, as well as shock wave transmission and reflection at the fluid-solid interface, which are key aspects of cavitation-induced damage. A vertical panel of height $h = 50$ mm and width $w = 1$ mm, made of steel, is placed in a shock tube on a rigid base (see Figure 4). Displacement and pressure sensors are positioned at the tip of the panel and on the top boundary, respectively. They are represented by the symbols (●) and (◇) in Figure 4. The dimensions, flow properties and material characteristics are detailed in Giordano et al. [50], who designed this FSI test case. The fluid mesh is refined near the panel, where $\Delta x_{min,f} = w/16$ and $\Delta y_{min,f} = w/8$. The fluid time step is set to $\Delta t_f = 5.75 \times 10^{-8}$ s. The coupling time window is set to $\Delta T = 1.15 \times 10^{-5}$ s, corresponding to approximately $T_s/250$, where $T_s$ is the oscillation period of the structure. In order to emphasize the effect of FSI coupling, both weakly (1-way) and strongly (2-way) coupled simulations are conducted.

Horizontal displacement of the panel tip at sensor (●) is illustrated in Figure 5a. The present results show an excellent agreement with the reference solution of Giordano et al. [50], which is also obtained
<figure><img src="image_5.png" /><figcaption>Figure 4: Shock-panel interaction: initial configuration of the problem. Pressure and displacement sensors are illustrated by symbols (◇) and (●), respectively.</figcaption></figure>

using a 2-way coupling. The pressure signal obtained at sensor (◇) with a 2-way coupling is presented in Figure 5b. It shows a very good agreement with the results of Giordano *et al.* [50] and Pasquariello *et al.* [51], indicating that the reflection and transmission of shock waves at the fluid-solid interface are correctly reproduced by the penalization-based strong partitioned coupling.

<figure><img src="image_6.png" /><figcaption>Figure 5: Numerical validation of the penalization-based strong partitioned coupling: (a) horizontal displacement of the panel at sensor (●); (b) pressure signal at sensor (◇).</figcaption></figure>

# 4. Cavitation pit formation

## 4.1. Problem setup

The present coupling strategy is applied to the numerical simulation of shock-induced bubble collapse near an elastoplastic material, aiming to reproduce cavitation damage at the microscopic scale. The initial conditions are similar to those presented by Paquette *et al.* [26]. A gas bubble, with an initial radius of $R_0 = 0.05$ mm, is immersed in liquid water at rest with initial conditions $p_0 = 10^5$ Pa and $T_0 = 298$ K.
The bubble collapse is triggered by an incident normal shock wave traveling at Mach number $M_{sh} = 1.04$, corresponding to an amplitude of $p_{sh} = 120$ MPa. We consider an elastoplastic material (Al7075) with the following parameters: Young's modulus $E = 70$ GPa, density $\rho_s = 2700$ kg/m$^3$, Poisson's ratio $\nu_s = 0.36$ and yield strength $\sigma_0 = 500$ MPa. Roy et al. [52] conducted compression tests to determine material properties of Al7075 in the plastic regime, characterized by a Johnson-Cook plasticity model [53]. Assuming low plastic strain in the present application ($\epsilon_{eq}^p < 0.1$, with $\epsilon_{eq}^p$ the equivalent plastic strain), the stress-strain relation can be approximated by a linear slope in the plastic regime. A linear isotropic hardening law is therefore adopted, where stress and strain are linearly related by the tangent modulus $E_t$. The latter is set to $E_t = E/100$ in this study, such that the stress-strain slope matches relatively well with experimental results of Roy et al. [52] and data from Paquette et al. [26]. This value of $E_t$ for Al7075 is also used by Hsiao et al. [21].

<figure><img src="image_7.png" /><figcaption>Figure 6: Schematic view of the initial configuration of the fluid-structure problem.</figcaption></figure>

The initial configuration of the fluid and solid states is shown in Figure 6. The non-dimensional stand-off distance between the bubble center and the wall is defined as $\gamma_0 = H_0/R_0$, where $H_0$ is the dimensional stand-off distance. This parameter is first set to $\gamma_0 = 2$. Only the upper half of the physical domain is computed, and a symmetry boundary condition is applied at $y/R_0 = 0$, taking advantage of the flow symmetry assumption. A non-reflecting condition is applied to the left boundary of the fluid domain, while a slip condition is imposed on its top boundary. Zero displacement in the normal direction is enforced on the top and right boundaries of the solid domain. Therefore, the solid domain must be large enough to prevent
wave reflections from influencing the solution near the coupled interface. The fluid and solid computational domain sizes are $[L_{x,f} \times L_{y,f}] = [6R_0 \times 20R_0]$ and $[L_{x,s} \times L_{y,s}] = [40R_0 \times 20R_0]$, respectively. Pressure and displacement probes are placed at the solid wall at points P0 (◇), P1 (◇), and P2 (◇), corresponding to vertical positions $y/R_0 = \{0, 2.5, 5\}$, respectively. The fluid domain is discretized using a stretched Cartesian grid with a uniform grid size $\Delta x_{min,f} = \Delta y_{min,f} = R_0/100$ in the area of interest. It results in a total number of 700 × 1053 cells. The CFL number is set to 0.2 by fixing the fluid time step to $\Delta t_f = 7.5 \times 10^{-11}$ s. The rectangular solid domain is discretized into a structured mesh of triangular elements with quadratic interpolation. Similarly to the fluid mesh, the node distribution is refined within the zone of interest, where $\Delta x_{min,s} = \Delta y_{min,s} = R_0/20$. The total number of degrees of freedom is about 350,000. In the present subcycling configuration, the structure time step and the coupling time window are set to $\Delta T = \Delta t_s = 5\Delta t_f$. The fluid mesh size, the structure mesh size and the coupling time window are determined based on numerical convergence analyses, which are presented in Appendix A.

## 4.2. Weak coupling

Weak coupling, or 1-way coupling, is first considered. This means that the fluid dynamics are calculated with the wall treated as rigid.

### 4.2.1. Shock-induced bubble collapse near a wall

The major events of the bubble collapse are illustrated in Figure 7 through instantaneous Schlieren visualizations (top) and pressure fields (bottom) at times $t = \{2.16, 2.46, 2.88\} \times 10^{-7}$ s. Due to the pressure difference induced by the passage of the incident normal shock wave, the bubble contracts asymmetrically. A high-speed liquid jet forms and penetrates the bubble horizontally towards the wall. When this jet impacts the opposite side of the bubble, an intense shock wave is emitted, also known as the water hammer shock. Meanwhile, the remaining bubble ring continues to collapse, leading to the emission of a secondary shock wave. The two shocks propagate towards the wall, impact it, and are reflected. The impact pressure exceeds 900 MPa, which is significantly higher than the yield strength of the material. Consequently, the material may undergo plastic deformation, explaining the formation of a cavitation pit at its surface.

The temporal evolution of the wall pressure at probes P0, P1, and P2 is shown in Figure 8. The incident shock impacts the wall at $t \sim 1.25 \times 10^{-7}$ s. This event is marked by an initial pressure jump of about 250 MPa, which is weaker at P0 as the bubble acts as a protective screen. The amplitude of the incident shock is much smaller than that of the shocks emitted during the bubble collapse. The water hammer shock impacts the wall at $t \sim 2.80 \times 10^{-7}$ s at point P0, generating a pressure peak of about 915 MPa. It is followed by the secondary shock, with slightly lower amplitude. They then continue to propagate along the wall. The water hammer shock and the secondary shock merge to form a single shock, which explains the presence of a single pressure peak at points P1 and P2. The shock reaches point P1 at $t \sim 3.25 \times 10^{-7}$ s, then point P2 at $t \sim 3.95 \times 10^{-7}$ s. The intensity and speed of the shock decrease as it propagates along the wall.
<figure><img src="image_8.png" /><figcaption>Figure 7: Instantaneous Schlieren visualizations (top) and pressure fields (bottom) at times $t = \{2.16, 2.46, 2.88\} \times 10^{-7}$ s.</figcaption></figure>

These aspects are marked by a decrease in the amplitude of the pressure peaks and by an increasing time gap between the peaks at points P0, P1, and P2. The pressure peaks exceed the yield strength of Al7075, $\sigma_0 = 500$ MPa, shown by a horizontal gray line in Figure 8. This indicates a potential material damage.

<figure><img src="image_9.png" /><figcaption>Figure 8: Temporal evolution of the wall pressure at probes P0, P1, and P2.</figcaption></figure>

The present results are compared to those proposed by Paquette *et al.* [54]. The amplitude of the pressure peaks is higher in the present study. The spatial resolution is similar between the two simulations, but the modeling of the multiphase flow is different. We use a four-equation formulation, whereas Paquette *et al.* [54] adopt the multicomponent model proposed by Johnsen *et al.* [10]. In their comparative study
of various one-fluid diffuse-interface formulations, Goncalves et al. [39] observe that the pressure peaks
are underestimated with the multicomponent model, compared to the present four-equation model and the
five-equation model of Kapila et al. [55].

4.2.2. Validation of the structure solver

This section focuses on the validation of the resolution of the structure subproblem, i.e., the dynamic elastoplasticity problem (see section 2.2). The present results, obtained with FEniCS, are compared with those from the commercial solver Abaqus. The solution from each solver is computed using the same linear isotropic hardening law, the same spatial resolution ($\Delta x_{min,s} = R_0/20$) and the same external loading. The contours of the cumulative plastic strain at the final time, obtained with FEniCS and Abaqus, are shown in Figures 9a and 9b. The topology of the damaged zone and the amplitude of the cumulative plastic strain are in excellent agreement. Quantitatively, a vertical cross-section of the cumulative plastic strain passing through the maximum is presented in Figure 10 for each solver. Excellent agreement is observed between the two solutions.

<figure><img src="image_10.png" /><figcaption>Figure 9: Numerical validation of the structure solver: contours of the cumulative plastic strain at the final time, obtained with (a) FEniCS and (b) Abaqus.</figcaption></figure>

4.2.3. Standard/Impedance penalization methods

The standard penalization method and its variant, the impedance penalization method (see section 2.1), are compared to the body-fitted method. The temporal evolution of the wall pressure is shown in Figure 11 for the three methods. The body-fitted method and the impedance penalization method produce very similar results. This indicates that the acoustic impedance, implicitly modeled by the body-fitted method, closely matches that of the actual material. Compared to the body-fitted method, the pressure peak obtained with
<figure><img src="image_11.png" /><figcaption>Figure 10: Numerical validation of the structure solver: vertical cross-section of the cumulative plastic strain passing through the maximum.</figcaption></figure>

the impedance penalization is slightly more intense at the impact point P0. However, this trend reverses at points P1 and P2. With the standard penalization method, the pressure peak is attenuated at P0 and P1 and delayed at P2. This behavior arises because the acoustic impedance of the penalized solid region is not controlled, as no penalization term is applied to the mass conservation equation. As a result, part of the shock wave is transmitted into the penalized region rather than being reflected.

<figure><img src="image_12.png" /><figcaption>Figure 11: Standard/Impedance penalization methods: temporal evolution of the wall pressure at probes P0, P1 and P2.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Standard/Impedance penalization methods: (a) contours of the cumulative plastic strain at the final time, obtained with the standard penalization method; (b) vertical cross-section of the cumulative plastic strain passing through the maximum.</figcaption></figure>

The contours of the cumulative plastic strain at the final time, obtained with the standard penalization method, are illustrated in Figure 12a. The damping of the impact observed with the standard penalization results in a significant reduction in plastic strain compared to the body-fitted method. This is shown in Figure 12b with the vertical cross-section of the cumulative plastic strain passing through the maximum. Very small variations are observed between the body-fitted method and the impedance penalization method. In contrast, the differences are much larger with the standard penalization method. The maximum plastic strain is reduced by nearly 30% compared to the body-fitted method. The discrepancy between the standard penalization and the other methods is more pronounced at the core of the damaged region. Indeed, Figure 11 shows that the attenuation of the pressure peak is also maximal near the damaged zone, at point P1. The impedance penalization method is retained for the remainder of this study.

### 4.3. Strong coupling

A strong coupling, or 2-way coupling, is now considered between the fluid and the structure. The influence of the wall deformation on the fluid flow is taken into account by updating the position of the IB and the mask function within the fluid solver (see Section 2.3).

#### 4.3.1. Preliminary results on the coupling effect

The time evolution of the wall pressure at points P0, P1, and P2 is shown in Figure 13. The results obtained with strong coupling are compared to those from weak coupling. Damping of the wall pressure
is observed at each point of the interface when 2-way coupling is taken into account. A reduction of 17%, 24%, and 10% in the peak pressure is observed at points P0, P1, and P2, respectively.

<figure><img src="image_14.png" /><figcaption>Figure 13: Preliminary results on the coupling effect: temporal evolution of the wall pressure at probes P0, P1 and P2.</figcaption></figure>

The contours of the cumulative plastic strain at the final time, obtained with the 1-way and 2-way approaches, are compared in Figure 14. A significant decrease in the plastic strain is observed when strong coupling is considered. Unlike points P0 and P2, the vicinity of point P1 is affected by plasticity. This explains the greater reduction in the pressure peak at point P1: the 1-way/2-way damping is enhanced by plastic deformation. Although the damping of the pressure peaks is between 10% and 24%, those of the maximum plastic strain and damaged surface area reach 41% and 54%, respectively. For the same damage criterion ($\epsilon_{eq}^p \ge 10^{-3}$), Paquette et al. [26] also observe a halving of the damaged area when strong coupling is adopted, compared to weak coupling. These significant differences highlight the importance of FSI coupling when characterizing material damage induced by bubble collapse.

### 4.3.2. Effect of stand-off distance

The initial distance between the bubble and the wall, also known as the stand-off distance, plays a key role in the intensity of the shock impact [34, 6, 10]. In the previous sections, the center of the bubble was initially located at a distance of $2R_0$ from the solid wall, i.e., $\gamma_0 = 2$ (see section 4.1). In this section, a parametric study is carried out with respect to the stand-off distance. Ten different initial configurations are considered such that $\gamma_0 = \{1.1, 1.2, \dots, 1.9, 2\}$. The collapse of the bubble is illustrated for $\gamma_0 = 1.1$ in Figure 15 through instantaneous Schlieren visualizations (top) and pressure field visualizations (bottom). The view is magnified by a factor of two around the bubble compared to Figure 7 for $\gamma_0 = 2$. In comparison with the bubble collapse presented in Figure 7, the water hammer occurs closer to the wall. The resulting
<figure><img src="image_15.png" /><figcaption>Figure 14: Preliminary results on the coupling effect: contours of the cumulative plastic strain at the final time.</figcaption></figure>

water hammer shock is much more intense and immediately impacts the fluid-solid interface. The complex interaction between the reflected water hammer shock and the collapsing bubble leads to a second impact on the wall.

<figure><img src="image_16.png" /><figcaption>Figure 15: Effect of stand-off distance, $\gamma_0 = 1.1$: instantaneous Schlieren visualizations (top) and pressure fields (bottom) at times $t = \{2.28, 2.52, 2.76\} \times 10^{-7}$ s.</figcaption></figure>

The time evolution of the wall pressure at points P0, P1, and P2 is presented in Figure 16 for $\gamma_0 = \{1.1, 1.4, 1.7, 2\}$. The impact occurs earlier for a closer bubble, as the shock wave has a shorter distance to travel before reaching the wall. The amplitude of the pressure peak at point P0 increases as the bubble gets closer to the wall. It is nearly doubled when $\gamma_0$ decreases from 2 to 1.1. A very important aspect concerns
the loading duration. Due to shock reflections between the wall and the bubble, the temporal extent of the pressure peak is greater for smaller stand-off distances. In other words, the characteristic loading time increases as $γ_0$ decreases. This is a fundamental factor in the material's response to the loading. The pressure peaks at P1 and P2 are much less affected by the stand-off distance.

<figure><img src="image_17.png" /><figcaption>Figure 16: Effect of stand-off distance: temporal evolution of the wall pressure at probes P0, P1 and P2.</figcaption></figure>

The contours of the cumulative plastic strain at the final time are illustrated in Figure 17 for $γ_0 = \{1.1, 1.4, 1.7\}$. The color scale is different for $γ_0 = 1.1$ because the amplitude of the plastic strain is much higher than for the other values of $γ_0$. The stand-off distance strongly influences both the maximum amplitude of the plastic strain and the topology of the damaged surface. The latter is localized near the impact position for a closer bubble ($γ_0 = 1.1$), then is progressively shifted away from the symmetry axis as the stand-off distance increases. Consequently, the damaged surface no longer coincides with the impact position, where the loading is nevertheless maximal (see Figure 16). Therefore, it is expected that a closer bubble forms a deeper, circular pit, while a more distant bubble creates a shallower, annular pit. This physical mechanism is also observed by previous authors [54, 23]. A physical explanation is provided by Joshi *et al.* [23], who attribute this mechanism to inertial effects within the solid. This interpretation aligns with the observations of Roy *et al.* [56], who studied the influence of loading duration on pit depth. When the impact duration is short, the loading frequency is very high, and the deformation is governed by the material's inertia rather than its stiffness. Most of the energy is therefore converted into kinetic
<figure><img src="image_18.png" /><figcaption>Figure 17: Effect of stand-off distance: contours of the cumulative plastic strain at the final time.</figcaption></figure>

energy, leading to minimal plastic deformation. For a distant bubble, the characteristic loading time is much shorter than the material's characteristic time during impact. As a result, the material does not immediately respond to the impact. Roy *et al.* [56] define the material's characteristic time *a posteriori* as the ratio of the damaged zone depth to the speed of sound in the solid. In practice, the material's characteristic time is difficult to estimate because the depth of the damaged zone depends on a damage criterion. There is no consensus on the choice of this criterion in the literature.

<figure><img src="image_19.png" /><figcaption>Figure 18: Effect of stand-off distance: maximum cumulative plastic strain as a function of the stand-off distance.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Effect of yield strength: stress-strain curve for each material considered.</figcaption></figure>

The maximum plastic strain is illustrated as a function of the stand-off distance in Figure 18. Two main regimes are identified. In the case of a distant bubble ($\gamma_0 \ge 1.4$), the maximum plastic strain is quite low.
According to Figure 16, maximum wall pressure exhibits a rapid decrease as the shock front propagates along the wall. Since the material's entry into plasticity is delayed for distant bubbles, plasticity occurs only after the loading has been significantly damped, resulting in a notably reduced plastic strain. Maximum plastic strain also shows very small variation for $\gamma_0 \ge 1.4$. Figure 16 highlights that the decay in peak pressure becomes more pronounced as the stand-off distance decreases, due to the shock reflection occurring at a smaller angle. As a result, the difference in wall pressure for varying stand-off distances decreases as we approach point P1. Combined with the delayed onset of plasticity as $\gamma_0$ increases, this leads to minimal variation in the maximum plastic strain for $\gamma_0 \ge 1.4$. For a closer bubble ($\gamma_0 < 1.4$), the maximum plastic strain increases exponentially as the stand-off distance decreases. Since the loading duration increases as the stand-off distance decreases, plasticity occurs earlier for a closer bubble, when the shock amplitude has not yet been damped. Additionally, the loading amplitude increases as the stand-off distance decreases (see Figure 16). This results in a much higher loading at the onset of material yielding, and consequently, a much greater maximum plastic strain.

### 4.3.3. Effect of yield strength

The effect of material's yield strength is now investigated, considering $\sigma_0 = \{400, 500, 750\}$ MPa. The other mechanical properties remain unchanged. The stress-strain curve is illustrated in Figure 19 for each material. As described in section 4.1, stress and strain follow a linear isotropic hardening law in the plastic regime. The combined effect of yield strength with the stand-off distance is further analyzed for $\gamma_0 = \{1.1, 1.4, 1.7, 2\}$.

<figure><img src="image_21.png" /><figcaption>Figure 20: Effect of yield strength, $\gamma_0 = 2$: contours of the cumulative plastic strain at the final time.</figcaption></figure>

The temporal evolution of the wall pressure is not affected by the yield strength of the material and is
therefore not represented here. In contrast, the material's damage is strongly influenced by this parameter. The contours of the cumulative plastic strain for $\gamma_0 = 2$ are presented in Figure 20. They show that material damage increases as the yield strength decreases. For $\sigma_0 = 750$ MPa, the damage is almost eliminated. The area of the damaged surface increases dramatically as $\sigma_0$ decreases, growing twentyfold when the yield strength drops from 750 MPa to 400 MPa. Its evolution is almost independent of $\gamma_0$, as the yield strength's effect on the damaged surface area dominates over that of the stand-off distance. However, the stand-off distance alters the impact of $\sigma_0$ on the maximum plastic strain. For $\gamma_0 \ge 1.4$, the maximum plastic strain exhibits relatively small variations with respect to $\sigma_0$. In contrast, for $\gamma_0 = 1.1$, the yield strength has a much stronger influence on the maximum plastic strain, which increases significantly faster as the yield strength decreases. This mechanism arises from the same principle described in the previous section: when plasticity occurs early, the loading amplitude at its onset is higher due to limited damping. This leads to a significant increase in the maximum plastic strain. This phenomenon is particularly observed for low yield strengths and closer bubbles.

The plastic deformation of the structure results in the formation of a pit on the material surface. The deformed interface for each $(\sigma_0, \gamma_0)$ configuration is shown in Figure 21. The depth of the pits is much smaller than their radius. In the most severe case ($\sigma_0 = 400$ MPa, $\gamma_0 = 1.1$), the pit reaches a depth of approximately $R_0/33$, with its radius slightly smaller than $R_0$. Convex deformation of the material surface appears in the region surrounding the pit. This phenomenon typically occurs in metals and is commonly referred to as *pile-up* [57]. During cavitation impact, the material is pushed radially outward from the impact point, causing the region surrounding the pit to slightly rise. This mechanism is observed in both experiments [36, 37] and numerical simulations [21]. Figure 21 shows that the extent of *pile-up* increases with pit depth.

The pit depth increases as the yield strength decreases. It is approximately multiplied by two when the yield strength is reduced from 500 MPa to 400 MPa. No noticeable pit is formed for $\sigma_0 = 750$ MPa. As the stand-off distance increases, the pit depth decreases rapidly, while its radius increases. The pit even takes on an annular shape for $\gamma_0 = 2$. This confirms that closer bubbles create a deep, circular pit, while the most distant bubbles can generate a shallower, annular pit.

The effect of yield strength on pit depth and pit area is presented in more detail in Figures 22a and 22b, respectively. The pit depth, denoted as $d$, is defined as the maximum permanent displacement of the wall. The pit area, denoted as $A$, is defined as the area of the newly formed fluid region where the permanent displacement of the interface is positive. The definition of these geometric characteristics is illustrated in Figure 21a. Figure 22a shows a rapid increase in pit depth as yield strength decreases. This increase becomes more pronounced as the stand-off distance decreases. Similarly, Figure 22b shows a rapid increase in pit area with decreasing yield strength. However, the influence of stand-off distance on the growth rate of pit area is less significant compared to that on pit depth.
<figure><img src="image_22.png" /><figcaption>Figure 21: Geometry of the cavitation pits. The initial horizontal position of the wall is denoted as $x_0$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 22: Effect of yield strength on the pit characteristics for different stand-off distances.</figcaption></figure>

### 4.3.4. Insights into the coupling effect

It is well known that strong coupling attenuates the impact pressure and the resulting material damage compared to weak coupling [26, 27]. However, the extent and underlying mechanisms of this damping remain poorly understood. This section aims to provide a better understanding of the differences between weak
and strong coupling, in terms of impact pressure and pit characteristics. This is achieved by assembling the results obtained for the different stand-off distances and yield strengths. The following notation is used to express the relative difference between a quantity q obtained with strong coupling (2-way) and weak coupling (1-way):

$$ \Delta q = \frac{q_{2w} - q_{1w}}{q_{1w}}, \qquad (19) $$

where $q_{1w}$ and $q_{2w}$ are obtained from a 1-way coupling and a 2-way coupling, respectively.

Let us first investigate the damping of the wall pressure between 1-way and 2-way computations. It is denoted as $\Delta p_w$ and follows the definition of equation (19). In the simplified one-dimensional problem of an incident shock wave of relatively low amplitude impacting a purely elastic material in its normal direction, the damping of the impact pressure is supposed to be determined by the ratio of the acoustic impedances of the liquid and solid. Moreover, in a weakly coupled computation, the wall is assumed to be rigid, which means that its acoustic impedance is infinite. Under the above-mentioned assumptions, the acoustic wave theory of Brekhovskikh and Godin [58] states that the damping of the impact pressure between a rigid wall and a deformable wall with a finite impedance $Z_s$ is given by:

$$ \Delta p_w^{th} = -\frac{1}{1 + Z_s/Z_f}. \qquad (20) $$

with $Z_s = \sqrt{\rho_s E/(1 - \nu_s^2)}$ the acoustic impedance of the solid under plane strain assumption and $Z_f = \rho_l c_l$ the acoustic impedance of the ambient liquid. In the present configuration, equation (20) provides an analytical value of $\Delta p_w^{th} = -9\%$. In the shock-induced bubble collapse, these assumptions hold true at the very beginning of the collapse, when the normal incident shock wave impacts the wall after passing through the bubble. At this moment, the impact pressure is lower than the material's yield strength and no plasticity occurs: the material's behavior is fully elastic. In addition, the amplitude of the incident shock wave is much lower than that of the water hammer shock. The impact pressures obtained numerically with weak and strong coupling are 247 MPa and 227 MPa, respectively (see Figure 13). This leads to a numerical reduction in the impact pressure of 8%, which agrees well with the analytical value of 9%.

However, these assumptions are no longer valid during and after the impact of the water hammer shock on the wall. The latter is not normally incident, the shock amplitude is much higher and plasticity occurs in the material. Figure 23 illustrates the damping of the peak pressure, denoted as $\Delta p_w^{max}$, at points P0 and P1 for all values of $\gamma_0$ and $\sigma_0$. The damping behavior is very different between the two wall positions P0 and P1. The reduction in peak pressure remains constant at point P0 (see Figure 23a), whereas it varies with both parameters at point P1 (see Figure 23b). Since different stand-off distances and yield strengths lead to significant differences in the material's behavior, the constant damping observed in Figure 23a still appears to be governed by the ratio of the acoustic impedances of the liquid and solid. Cao et al. [25] investigated
numerically the effect of acoustic impedance on the shock-induced bubble collapse near solid materials. They observed that the same bubble collapsing near materials with the same acoustic impedance produces a similar impact pressure, even though the materials differed in density, Young's modulus, and Poisson's ratio. Furthermore, they reported a reduction of 14.5% in the maximum impact pressure between a rigid wall and a deformable wall made of Al7075. This corresponds relatively well with the present result, as our damping in the impact pressure at P0 is approximately 17% for each set of parameters. The present damping may be slightly increased by the discrete-time evolution of the mask function related to the penalization method.

<figure><img src="image_24.png" /><figcaption>Figure 23: Effect of coupling on the maximum wall pressure at points P0 and P1.</figcaption></figure>

Further from the point of impact, the reduction in peak pressure is locally influenced by the material's deformation. The damping is more significant in areas with greater deformation and less pronounced in regions with minimal deformation. Indeed, Figure 23b shows that the damping at point P1 increases with decreasing yield strength and increasing stand-off distance. The present parametric study has shown that plasticity is more pronounced at low yield strengths and that the damaged surface shifts radially towards point P1 as the stand-off distance increases. The deformation is thus increased because the slope of the stress-strain curve is much lower in the plastic regime than in the elastic regime, as illustrated in Figure 19.

In summary, wall deformation leads to the damping of the wall pressure in the case of strong coupling, compared to weak coupling, where the wall is treated as rigid. This attenuation is primarily determined by the ratio of acoustic impedances between the fluid and solid media. It is further amplified locally by plasticity or, more generally, in regions of greater deformation.

In turn, material damage is significantly reduced when strong coupling is considered. Figure 24 illustrates the pit geometry for 1-way and 2-way coupling for different stand-off distances, with $\sigma_0 = 500$ MPa. It shows
that the pit depth and the pit area are significantly decreased when strong coupling is considered.

The non-dimensional pit depth is shown in Figure 25a for all values of $\gamma_0$, with $\sigma_0 = 500$ MPa. Its relative difference between weakly and strongly coupled computations, denoted as $\Delta d$, is plotted in blue and associated with the right axis. Although the pit depth increases rapidly with decreasing stand-off distance, its attenuation remains almost constant, ranging from 55% to 59%. Similarly, the non-dimensional pit area and its associated damping, denoted as $\Delta A$, are presented in Figure 25b. A nearly constant damping is also observed for $\gamma_0 = \{1.1, 1.4, 2\}$, ranging from 70% to 73%. For $\gamma_0 = 1.7$, the 1-way/2-way reduction in the pit area is slightly less pronounced, with $\Delta A = -60\%$.

<figure><img src="image_25.png" /><figcaption>Figure 24: Effect of coupling on the pit geometry for different stand-off distances, with $\sigma_0 = 500$ MPa. The initial horizontal position of the wall is denoted as $x_0$.</figcaption></figure>

Two key observations can be made regarding the 1-way/2-way damping of the wall pressure and the resulting pit characteristics:

* The reduction in the pit depth and area is much more important than the damping in the wall pressure (by a factor of three).

* While the reduction in wall pressure varies locally with stand-off distance and yield strength, the decrease in pit size remains relatively unchanged.

## 5. Conclusions

We developed a novel strong partitioned coupling strategy in order to address FSI problems. The penalization method is extended to compressible multiphase flows and deformable walls and is used to model the deformable fluid-solid interface. This paper aims to apply this FSI coupling methodology to the
<figure><img src="image_26.png" /><figcaption>Figure 25: Effect of coupling on the pit characteristics for different stand-off distances, with $\sigma_0$ = 500 MPa.</figcaption></figure>

numerical analysis of cavitation-induced damage. We focused on the formation of a cavitation pit at the microscopic scale, generated by the shock-induced collapse of a single bubble near an elastoplastic material. The present paper offers insights into the fundamental physical mechanisms responsible for cavitation-induced damage, as well as a deeper understanding of the numerical modeling of this complex multiphysics problem. This is achieved through two crossed parametric studies on the initial bubble-wall distance (or stand-off distance) and the material's yield strength. We examined the effect of these parameters on wall pressure, material damage and pit characteristics.

The amplitude and duration of the wall pressure strongly depend on the stand-off distance, resulting in significant variations in the material damage. Closer bubbles produce damage concentrated near the impact point, while increasing the stand-off distance causes the damaged surface to shift radially and reduces the amplitude of the maximum plastic strain. This mechanism, consistent with the findings of Joshi et al. [23], is attributed to inertial effects in the material. As the stand-off distance increases, the impact duration decreases and the deformation is thus governed by the material's inertia rather than its stiffness. Additionally, for closer bubbles ($\gamma_0$ < 1.4), the maximum plastic strain grows exponentially as the stand-off distance decreases. Since the maximum wall pressure decreases rapidly as the shock front propagates along the wall, the material undergoes much higher loading at the onset of material yielding for closer bubbles. Consequently, closer bubbles generate a deeper, circular pit, whereas more distant bubbles can produce a shallower, annular pit. Both pit depth and pit area increase rapidly as the stand-off distance and yield strength decrease.

The standard penalization is not able to properly capture the wall pressure, resulting in an error of 30% in the maximum plastic strain. To overcome this issue, we developed a novel approach, called the impedance
penalization method [42]. It provides an accurate modeling of the shock wave transmission and reflection at the fluid-solid interface by controlling the acoustic impedance of the penalized solid region. Furthermore, this work offers new insights into the role of FSI coupling in the context of cavitation-induced damage. In strong coupling, wall deformation results in the damping of wall pressure. At the moment of impact, this damping is first governed by the ratio of acoustic impedances of fluid and solid media. It is further enhanced locally by plasticity or, more generally, in regions undergoing higher deformation. A slight damping in wall pressure results in a much more important reduction in material damage, and consequently, in pit depth and pit area. Although reduction in wall pressure is locally influenced by solid deformation, the decrease in pit size shows very small variation regarding stand-off distance and yield strength.

Future work will be dedicated to the investigation of three-dimensional (3D) and multi-bubble configurations. The fluid solver is already efficient and accurate in 3D [29, 59], and the extension of the penalization-based FSI coupling is straightforward. The main challenge thus lies in extending the structure solver, with the migration from FEniCS to FEniCSx [60], whose architecture is better suited for high-performance computing.

Beyond this study, the present numerical tool could contribute to the development of an innovative solution to mitigate cavitation erosion.

# 6. Acknowledgements

This work was supported by the France 2030 program (LABEX INTERACTIFS, reference ANR-11-LABX-0017-01) and the Agence Innovation Defense (AID).

# Appendix A. Numerical verification

**Fluid mesh size.** Numerical convergence analysis of the fluid solver for the shock-induced bubble collapse near a rigid wall is available in [29, 13].

**Structure mesh size.** This section examines the influence of the structure mesh size on material damage. The plastic strain of the material induced by the bubble collapse is analyzed by considering four progressively finer meshes. The mesh size is based on the initial bubble radius. The different meshes are presented in Table A.1. The coarse, medium, fine and very fine meshes are composed of $N_{x,s} = \{25, 50, 100, 200\}$ and $N_{y,s} = \{50, 100, 200, 400\}$ elements in each direction within the refined area, respectively.

The contours of the cumulative plastic strain are presented in Figure A.26 at the final time for both the coarse and fine meshes. The amplitude of the cumulative plastic strain at the final time is observed to increase with spatial resolution.

More quantitatively, the maximum cumulative plastic strain and the damaged surface area, $A_p$, are presented in Table A.1 for each mesh. The relative error with respect to the reference value, obtained from the
<figure><img src="image_27.png" /><figcaption>Figure A.26: Numerical spatial convergence of the structure solver: contours of the cumulative plastic strain at the final time.</figcaption></figure>

<table><caption>Table A.1: Numerical spatial convergence of the structure solver: maximum cumulative plastic strain, damaged surface area and computation time.</caption><thead><tr><th>Mesh</th><th>Coarse</th><th>Medium</th><th>Fine</th><th>Very fine</th></tr></thead><tbody><tr><td>Δx<sub>min,s</sub>, Δy<sub>min,s</sub></td><td>R<sub>0</sub>/5</td><td>R<sub>0</sub>/10</td><td>R<sub>0</sub>/20</td><td>R<sub>0</sub>/40</td></tr><tr><td>ε<sup>p,max</sup><sub>eq</sub> [-]</td><td>0.74 × 10<sup>-2</sup></td><td>0.99 × 10<sup>-2</sup></td><td>1.08 × 10<sup>-2</sup></td><td>1.11 × 10<sup>-2</sup></td></tr><tr><td>A<sub>p</sub> [mm<sup>2</sup>]</td><td>1.45 × 10<sup>-2</sup></td><td>1.53 × 10<sup>-2</sup></td><td>1.49 × 10<sup>-2</sup></td><td>1.48 × 10<sup>-2</sup></td></tr><tr><td>Err(ε<sup>p,max</sup><sub>eq</sub>)</td><td>33.3%</td><td>10.8%</td><td>2.7%</td><td>-</td></tr><tr><td>Err(A<sub>p</sub>)</td><td>2.0%</td><td>3.4%</td><td>0.7%</td><td>-</td></tr><tr><td>Comput. time (structure)</td><td>18 min</td><td>3 h 25 min</td><td>12 h 20 min</td><td>36 h 55 min</td></tr></tbody></table>

very fine mesh, is also provided. The maximum cumulative plastic strain increases with spatial resolution. Its relative error decreases as the resolution increases, reaching only 2.7% for the fine mesh. The damaged surface area is computed based on the damage criterion $ε_{eq}^p ≥ 10^{-3}$, consistent with the study of Paquette et al. [54]. The evolution of the damaged surface area with respect to spatial resolution is not monotonic, but it converges toward $A_p = 1.48 × 10^{-2}$ mm$^2$, obtained with the very fine mesh. The relative error in the damaged surface area is reduced to 0.7% for the fine mesh. The computation time of the structure solver is also presented in Table A.1 for each mesh. It becomes extremely high and even prohibitive for the very fine mesh. Based on these observations, the fine mesh is retained in this paper.

**Coupling time window.** The independence of the solution with respect to the time window is verified using the structure fine mesh. The material's response is analyzed for two distinct time windows, $ΔT = 5Δt_f$ and $ΔT = 2Δt_f$. The cumulative plastic strain remains relatively unaffected as $ΔT$ varies, with variations
in $\epsilon_{eq}^{p,max}$ and $A_p$ being only 2.9% and 0.7%, respectively. However, reducing the time window to $2\Delta t_f$
increases the computational cost by approximately a factor of 2.5. Based on these observations, the time
window is set to $\Delta T = 5\Delta t_f$ in the present paper.

References
---

[1] G. L. Chahine, J.-P. Franc, A. Karimi, Cavitation and Cavitation Erosion, Springer Netherlands, Dordrecht, 2014, pp. 3–20.

[2] T. G. Leighton, R. O. Cleveland, Lithotripsy, Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 224 (2) (2010) 317–342.

[3] C. E. Brennen, Cavitation in medicine, Interface Focus 5(5) (2015) 20150022.

[4] C. D. Ohl, M. Arora, R. Ikink, N. de Jong, M. Versluis, M. Delius, D. Lohse, Sonoporation from jetting cavitation bubbles, Biophysical Journal 91 (11) (2006) 4285–4295.

[5] H. Soyama, Cavitating jet: A review, Applied Sciences 10 (20) (2020) 7280.

[6] Y. Tomita, A. Shima, Mechanisms of impulsive pressure generation and damage pit formation by bubble collapse, Journal of Fluid Mechanics 169 (1986) 535–564.

[7] A. Philipp, W. Lauterborn, Cavitation erosion by single laser-produced bubbles, Journal of Fluid Mechanics 361 (1998) 75–116.

[8] M. S. Plesset, R. B. Chapman, Collapse of an initially spherical vapour cavity in the neighbourhood of a solid boundary, Journal of Fluid Mechanics 47 (1971) 283–290.

[9] C. K. Turangan, A. R. Jamaluddin, G. J. Ball, T. G. Leighton, Free-lagrange simulations of the expansion and jetting collapse of air bubbles in water, Journal of Fluid Mechanics 598 (2008) 1–25.

[10] E. Johnsen, T. Colonius, Numerical simulations of non-spherical bubble collapse, Journal of Fluid Mechanics 629 (2009) 231–262.

[11] N. A. Hawker, Y. Ventikos, Interaction of a strong shockwave with a gas bubble in a liquid medium: a numerical study, Journal of Fluid Mechanics 701 (2012) 59–97.

[12] F. Wermelinger, U. Rasthofer, P. E. Hadjidoukas, P. Koumoutsakos, Petascale simulations of compressible flows with interfaces, Journal of Computational Science 26 (2018) 217–225.

[13] E. Goncalves, P. Parnaudeau, Numerical study of pressure loads generated by a shock-induced bubble collapse, Physics of Fluids 33 (2021) 113318.

[14] D. C. Gibson, J. R. Blake, The growth and collapse of bubbles near deformable surfaces, Applied Scientific Research 38 (1982) 215–224.

[15] A. Shima, Y. Tomita, D. C. Gibson, J. R. Blake, The growth and collapse of cavitation bubbles near composite surfaces, Journal of Fluid Mechanics 203 (1989) 199–214.

[16] E. A. Brujan, K. Nahen, P. Schmidt, A. Vogel, Dynamics of laser-induced cavitation bubbles near an elastic boundary, Journal of Fluid Mechanics 433 (2001) 251–281.

[17] J. H. Duncan, C. D. Milligan, S. Zhang, On the interaction between a bubble and a submerged compliant structure, Journal of Sound and Vibration 197 (1) (1996) 17–44.

[18] G. N. Sankin, P. Zhong, Interaction between shock wave and single inertial bubbles near an elastic boundary, Physical Review E 74 (2006) 046304.

[19] H. Miao, S. M. Gracewski, Coupled FEM and BEM code for simulating acoustically excited bubbles near deformable structures, Computational Mechanics 42 (2008) 95–106.

[20] T. Zhao, X. Zhao, X. Ma, B. Huang, Numerical investigations of the bubble collapse near an elastic plate, Journal of Fluids and Structures 106 (2021) 103372.

[21] C. T. Hsiao, A. Jayaprakash, A. Kapahi, J. K. Choi, G. L. Chahine, Modelling of material pitting from cavitation bubble collapse, Journal of Fluid Mechanics 755 (2014) 142–175.

[22] C. K. Turangan, G. J. Ball, A. R. Jamaluddin, T. G. Leighton, Numerical studies of cavitation erosion on an elastic-plastic material caused by shock-induced bubble collapse, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473 (2017) 20170315.

[23] S. Joshi, J.-P. Franc, G. Ghigliotti, M. Fivel, SPH modelling of a cavitation bubble collapse near an elasto-visco-plastic material, Journal of the Mechanics and Physics of Solids 125 (2019) 420–439.

[24] Q. Zhang, X. Liu, X. Dong, L. Yin, Z. Cheng, A volume-adaptive mesh-free model for FSI simulation of cavitation erosion with bubble collapse, Computational Particle Mechanics 11 (2024) 2325–2351.

[25] S. Cao, G. Wang, O. Coutier-Delgosha, K. Wang, Shock-induced bubble collapse near solid materials: effect of acoustic impedance, Journal of Fluid Mechanics 907 (2021) A17.

[26] Y. Paquette, M. Fivel, G. Ghigliotti, E. Johnsen, J.-P. Franc, Fluid-structure interaction in cavitation erosion, in: Proceedings of the 10th Int. Symposium on Cavitation, CAV2018, Baltimore, USA, 2018.

[27] P. Sarkar, G. Ghigliotti, J.-P. Franc, M. Fivel, Fluid-structure modelling for material deformation during cavitation bubble collapse, Journal of Fluids and Structures 106 (2021) 103370.

[28] R. Firly, K. Inaba, F. Triawan, K. Kishimoto, H. Nakamoto, Numerical study of impact phenomena due to cavitation bubble collapse on metals and polymers, European Journal of Mechanics - B/Fluids 101 (2023) 257–272.

[29] R. Dubois, E. Goncalves, P. Parnaudeau, High performance computing of stiff bubble collapse on CPU-GPU heterogeneous platform, Computers & Mathematics with Applications 99 (2021) 246–256.
735

736

737

738

739

740

741

742

743

744

745

746

747

748

749

750

751

752

753

754

755

756

757

758

759

760

761

762

763

764

765

766

767

768

769

770

771

772

773

774

775

776

777

778

779

780

781

782

783

784

785

786

787

788

789

790

791

792

793

794

795

[30] A. Logg, G. N. Wells, DOLFIN: Automated finite element computing, ACM Transactions on Mathematical Software 37 (2) (2010) 1–28.

[31] G. Chourdakis, K. Davis, B. Rodenberg, M. Schulte, F. Simonis, B. Uekermann, G. Abrams, H. J. Bungartz, L. Cheung Yau, I. Desai, K. Eder, R. Hertrich, F. Lindner, A. Rusch, D. Sashko, D. Schneider, A. Totounferoush, D. Volland, P. Vollmer, O. Z. Koseomur, preCICE v2: A sustainable and user-friendly coupling library [version 2; peer review: 2 approved], Open Research Europe 2 (51) (2022).

[32] E. Arquis, J. P. Caltagirone, Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieu poreux: Application à la convection naturelle, Comptes Rendus de l'Académie des Sciences de Paris 299 (1984) 1–4.

[33] P. Angot, C. H. Bruneau, P. Fabrie, A penalization method to take into account obstacles in incompressible viscous flows, Numerische Mathematik 81 (1999) 497–520.

[34] A. Shima, K. Takayama, Y. Tomita, N. Ohsawa, Mechanism of impact pressure generation from spark-generated bubble collapse near a wall, AIAA Journal 21 (1) (1983) 55–59.

[35] C. F. Naudé, A. T. Ellis, On the mechanism of cavitation damage by nonhemispherical cavities collapsing in contact with a solid boundary, Journal of Basic Engineering 83 (1961) 648–656.

[36] M. Dular, T. Požar, J. Zevnik, R. Petkovšek, High speed observation of damage created by a collapse of a single cavitation bubble, Wear 418-419 (2019) 13–23.

[37] H. J. Sagar, O. el Moctar, Dynamics of a cavitation bubble near a solid surface and the induced damage, Journal of Fluids and Structures 92 (2020) 102799.

[38] E. Goncalves, B. Charrière, Modelling for isothermal cavitation with a four-equation model, International Journal of Multiphase Flow 59 (2014) 54–72.

[39] E. Goncalves, P. Parnaudeau, Comparison of multiphase models for computing shock-induced bubble collapse, International Journal of Numerical Methods for Heat & Fluid Flow 30 (2020) 3845–3877.

[40] O. Le Metayer, J. Massoni, R. Saurel, Elaborating equations of state of a liquid and its vapor for two-phase flow models, International Journal of Thermal Sciences 43 (2004) 265–276.

[41] G. Wallis, One-dimensional two-phase flow, New York: McGraw-Hill (1967).

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

[43] P. Angot, Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows, Mathematical Methods in The Applied Sciences 22 (1999) 1395–1412.

[44] A. Piquet, O. Roussel, A. Hadjadj, A comparative study of Brinkman penalization and direct-forcing immersed boundary methods for compressible viscous flows, Computers & Fluids 136 (2016) 272–284.

[45] Q. Liu, O. V. Vasilyev, A Brinkman penalization method for compressible flows in complex geometries, Journal of Computational Physics 227 (2007) 946–966.

[46] N. M. Newmark, A method of computation for structural dynamics, Journal of the Engineering Mechanics Division 85 (3) (1959) 67–94.

[47] R. D. Krieg, D. B. Krieg, Accuracies of numerical solution methods for the elastic-perfectly plastic model, Journal of Pressure Vessel Technology 99 (4) (1977) 510–515.

[48] B. Rodenberg, I. Desai, R. Hertrich, A. Jaust, B. Uekermann, FEniCS–preCICE: Coupling FEniCS to other simulation software, SoftwareX 16 (2021) 100807.

[49] S. D. Roth, Ray casting for modeling solids, Computer Graphics and Image Processing 18 (2) (1982) 109–144.

[50] J. Giordano, G. Jourdan, Y. Burtschell, M. Medale, D. E. Zeitoun, L. Houas, Shock wave impacts on deforming panel, an application of fluid-structure interaction, Shock Waves 14 (2005) 103–110.

[51] V. Pasquariello, G. Hammerl, F. Örley, S. Hickel, C. Danowski, A. Popp, W. A. Wall, N. A. Adams, A cut-cell finite volume–finite element coupling approach for fluid–structure interaction in compressible flow, Journal of Computational Physics 307 (2016) 670–695.

[52] S. C. Roy, J.-P. Franc, M. Fivel, Cavitation erosion: Using the target material as a pressure sensor, Journal of Applied Physics 118 (16) (2015) 164905.

[53] G. R. Johnson, W. H. Cook, Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures, Engineering Fracture Mechanics 21 (1) (1985) 31–48.

[54] Y. Paquette, Interaction fluide-structure et Érosion de cavitation, Thèse de doctorat, Université Grenoble Alpes (2017).

[55] A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son, D. S. Stewart, Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations, Physics of Fluids 13 (2001) 3002–3024.

[56] S. C. Roy, J.-P. Franc, N. Ranc, M. Fivel, Determination of cavitation load spectra—part 2: Dynamic finite element approach, Wear 344-345 (2015) 120–129.

[57] J.-P. Franc, G. L. Chahine, A. Karimi, Pitting and Incubation Period, Springer Netherlands, Dordrecht, 2014, pp. 37–69.

[58] L. M. Brekhovskikh, O. A. Godin, Acoustics of Layered Media I, 1st Edition, Springer, Berlin, Heidelberg, 1990.

[59] E. Goncalves da Silva, P. Parnaudeau, Wall pressure amplification of shock-collapsed multi-bubble arrays near a rigid wall, Physics of Fluids 37 (1) (2025) 013319.

[60] I. A. Baratta, J. P. Dean, J. S. Dokken, M. Habera, J. S. Hale, C. N. Richardson, M. E. Rognes, M. W. Scroggs, N. Sime, G. N. Wells, DOLFINx: The next generation FEniCS problem solving environment (2023). doi:10.5281/zenodo.10447666.