# A mixed-dimensional formulation for the simulation of slender structures immersed in an incompressible flow

Fabien Lespagnol, Céline Grandmont, Paolo Zunino, Miguel Angel Fernández

► To cite this version:

Fabien Lespagnol, Céline Grandmont, Paolo Zunino, Miguel Angel Fernández. A mixed-dimensional formulation for the simulation of slender structures immersed in an incompressible flow. Computer Methods in Applied Mechanics and Engineering, 2024, 432, pp.117316. 10.1016/j.cma.2024.117316. hal-04318526v2

HAL Id: hal-04318526

https://inria.hal.science/hal-04318526v2

Submitted on 1 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A mixed-dimensional formulation for the simulation of slender structures immersed in an incompressible flow

FABIEN LESPAGNOL

MOX, Department of Mathematics, Politecnico di Milano, Milan, Italy

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France

CÉLINE GRANDMONT

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France

Département de Mathématique, Université libre de Bruxelles, Brussels, Belgium

PAOLO ZUNINO

MOX, Department of Mathematics, Politecnico di Milano, Milan, Italy

MIGUEL A. FERNÁNDEZ

Sorbonne Université, Inria, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Paris, France

## Abstract

We consider the simulation of slender structures immersed in a three-dimensional (3D) flow. By exploiting the special geometric configuration of the slender structures, this particular problem can be modeled by mixed-dimensional coupled equations. Taking advantage of the slenderness of the structure and thus considering 3D/1D coupled problems raise several challenges and difficulties. From a mathematical point of view, these include defining well-posed trace operators of co-dimension two. On the computational standpoint, the non-standard mathematical formulation makes it difficult to ensure the accuracy of the solutions obtained with the mixed-dimensional discrete formulation as compared to a fully resolved one. Here we proposed to circumvent these issues by imposing the fluid-structure coupling conditions on the 2D fluid-structure interface but in a reduced way still taking advantage of the 1D dynamic of the slender structure. We consider the Navier-Stokes equations for the fluid and a Timoshenko beam model for the structure. We complement these models with a mixed-dimensional version of the fluid-structure interface conditions, based on the projection of kinematic coupling conditions on a finite-dimensional Fourier space on each beam cross section. Furthermore, we develop a discrete fictitious domain formulation within the framework of the finite element method, establish the energy stability of the scheme, provide extensive numerical evidence of the accuracy of the discrete formulation, notably with respect to a fully resolved (ALE based) model and a standard reduced modeling approach.

**Keywords:** Fluid-structure interaction, slender structures, mixed-dimensional model, fictitious domain method.

## 1. Introduction

The modeling and simulation of slender structures immersed in three-dimensional (3D) flows is of paramount importance in many engineering fields. With the aim of understanding the behavior of such complex coupled systems, ranging from parachutes and bridge cables [1, 2, 3, 4] to submerged structures or vegetation [5, 6, 7, 8, 9, 10, 11], numerical simulations provide insights which are unattainable by experimentation alone. Its interest also extends to biomedical applications, such as the design of vascular stents [12, 13, 14] and the simulation of cilia, flagella [15] and micro-swimmers [16, 17, 18]. Owing to the geometrical properties of slender structures, many models have exploited the mixed-dimensional nature of these coupled systems with the purpose of developing more efficient computational approaches, in which the solid balance equations are written in a one-dimensional (1D) domain, using for instance geometrically exact beam theory (see, e.g., [19, 20, 21, 22, 23]). This yields a particular case of fluid-structure interaction problems known as fluid-beam interaction. Coupled systems in which 2D flows mechanically interact with 1D beams can be efficiently simulated by identifying the beam centerline with the fluid-structure interface and by neglecting the beam thickness effects in the interface coupling (see, e.g., [24, 25]). This is a widely used modeling simplification when coupling general media featuring a
dimensional gap of one, viz., 2D-1D or 3D-2D coupling (see, e.g., [26, 27]). However, for problems in which 3D flows meet 1D beams, the situation is much more delicate. The difference in dimensionality complicates the formulation of the load and motion transfer conditions and also calls for innovative numerical techniques to bridge the gap between local fluid dynamics near the beam and global behavior of the system. In fact, two distinct modeling approaches can be found in the literature. A first family of methods is based on the formulation of the interface coupling conditions on the physical (2D) surface of the beam (see, e.g., [28, 29, 30]), while the second one reduces the complexity of the problem by considering the interface as the (1D) centerline of the beam (see, e.g., [31, 32, 8, 33, 34]).

In order to reconstruct the coupling conditions on the 2D interface generated from the centerline of the fluid-beam, the fluid mesh then needs to correctly resolve the beam interface, with either matching or non-matching grids. The choice of the method is typically influenced by the degree of displacement exhibited by solid (see, e.g., [25]). Systems involving low- to moderate-interface deflections can be effectively simulated using fitted mesh techniques, which often rely on an arbitrary Lagrangian-Eulerian (ALE) description of the fluid. However, for large interface displacements, the ALE formalism can become cumbersome. In this case, a preferred approach consists in combining a Eulerian description of the fluid with an unfitted mesh discretization. In this approach, the fluid-structure interface deforms independently of a background fluid mesh. Notable among these methodologies is the Immersed Boundary method (IB) (see, e.g., [35, 36, 37]), in which the interaction between the two meshes is accomplished by using approximations of Dirac functions. In the context of 3D-1D fluid-structure interaction (FSI) problems, this method has found widespread application (see, e.g., [31, 32, 38, 8, 39]). Another category comprises methods that use Lagrange multipliers in lieu of Dirac functions, referred to as Fictitious Domain Methods (FD) (see, e.g., [24, 40, 41, 42]).

Since the radius of the slender structures is small with respect to their longitudinal length, another solution is to ignore this 2D interface. A solution mostly adopted in the context of bio-engineering applications involving viscous Stokes flow is the so-called slender-body theory (see, e.g. [43, 44, 45, 46, 47, 48, 49]). Slender body theory is a one-way approach that accounts for the influence of the fluid on the fibers, but not vice versa. Another approach is to directly impose the coupling conditions on the beam centerline. From a numerical perspective, given that the approximation space consists mainly of continuous functions, various techniques can be used such as penalty methods [33] or mortar-type methods [34]. One advantage of this method is that coarse fluid meshes can be used in the vicinity of the beam. However, at the continuous level, two questions naturally emerge. First, is the continuous problem well-posed? Indeed, imposing the kinematic and dynamic coupling conditions on the beam centerline requires the introduction of a trace operator of codimension two, thereby demanding regularity for the solution within the fluid domain. However, these conditions simultaneously introduce a one-dimensional Dirac operator on the right-hand side of the fluid equations, often preventing the solution from reaching the desired regularity. The second question concerns the fidelity of the model: is it a good approximation of the fully resolved problem when the thickness of the slender body tends to 0? This is the case in the context of co-dimension one (see [26]), as indicated above.

In this paper, we introduce a new modeling approach and numerical strategy for fluid-structure interaction systems in which a 3D fluid mechanically interacts with an immersed slender solid, modeled as a 1D solid. The basic idea consists in combining a fictitious domain methodology with a projection of the kinematic constraint onto a finite-dimensional space lying on the beam centerline (see [50] for a similar idea in the context of reinforced materials). More precisely, on each cross section of the beam, the original kinematic condition is replaced by a constraint (enforced via Lagrange multipliers) which acts only on the first *N* modes of the Fourier expansion of the velocity. This approach offers several advantages with respect to the methodologies mentioned above. First, it circumvents the trace issues on manifolds of co-dimension 2 arising when imposing the coupling conditions on the beam centerline (see, e.g., [33, 34]), thus yielding a well-posed mathematical problem and, for *N* ≥ 1, enabling the incorporation of rotational effects in the FSI coupling. Second, compared to full 3D immersed approaches, the 2*N* + 1 Lagrange multiplier unknowns lie in the dimensionally reduced 1D domain. Consequently, at the discrete level, the space of Lagrange multipliers does not depend anymore on the interface mesh, which is introduced here only for integration purposes. This is particularly advantageous for small beam cross-sections where enforcing the kinematic coupling conditions on the two-dimensional interface would require a very fine fluid mesh and a large number of Lagrange multipliers (potentially leading to locking phenomena if the fluid mesh is not sufficiently refined). To some extent, this approach could be seen as a bridge between the aforementioned 2D and 1D representations of the interface, with the purpose of reaching a compromise between accuracy, computational cost and robustness.
A similar approach has been considered in [51] with $N = 0$ for elliptic problems. In particular, the work demonstrated the well-posedness of both the continuous and discrete models within the framework of the finite element method. Moreover, in the 2D case, following a framework similar to [51] but for an arbitrary number of Fourier modes, the authors of [52] study the properties and limitations of the method and show that the difference between the solution of the reduced mixed-dimensional model and the full solution tends to zero when the radius of the slender structure tends to zero. They also obtain the convergence of the Lagrange multiplier provided that $N \ge 1$. In the context of fluid-structure interaction, this point is crucial to correctly capture the fluid-structure dynamics. In order to illustrate the accuracy of the proposed strategy, we compare it with two other approaches, a fully resolved model based on the ALE approach and another reduced mixed-dimensional model with coupling conditions on the beam centerline. In particular, we will investigate the influence of the Reynolds number, of the fluid mesh refinement, or of the number of Fourier modes taken for the Lagrange multiplier space.

In summary, the main contributions of this paper are threefold: (i) the introduction of a 3D-1D FSI model for slender bodies immersed in 3D flow with mathematically sound reduced interface conditions; (ii) the derivation of energy estimates in the continuous and discrete framework; and (iii) the presentation of quantitative and qualitative comparisons with alternative approaches from the literature.

The rest of the paper is organized as follows. In Section 2, we present the geometrical setting and the full 3D coupled fluid-beam models. Section 3 introduces the proposed reduced modeling approach. The discrete formulations of the different approaches are presented in Section 4, including a discussion on the computer implementation of the proposed numerical method. Section 5 is devoted to the numerical experiments. Here, we investigate in a series of numerical examples the accuracy and robustness of the proposed approach with respect to the other two methodologies. Finally, a summary of the results of this work together with some research perspectives are given in Section 6.

## 2. Problem setting and full coupled problem

This section briefly presents the considered beam and fluid models, as well as the full coupled model of the resulting fluid-structure interaction system. We choose a linear Timoshenko beam formulation for the solid while the fluid is described by the incompressible Navier-Stokes equations.

### 2.1. The solid model

Several beam models have been developed in the literature for the simulation of slender structures, on the basis of different assumptions. Beam theories can be introduced through different perspectives, one of which focuses on the inclusion of non-linear effects related to large deformations and strains. One of the simplest and most practical theories, assuming small deformation and strains, introduced by Euler and Bernoulli around 1750 is the linear planar Euler-Bernoulli beam theory. The foundational assumption of this theory is that the beam's cross-sections are rigid and during deformation, the cross-section stay normal to the beam centerline. A possible enrichment of the Euler-Bernoulli model is the linear planar Timoshenko model [53, 54], which still assumes undeformable sections but can eventually not be perpendicular to the beam centerline. Without loss of generality and in order to ease the presentation, in this work we consider an extension of the linear planar Timoshenko beam theory to three dimensions (see [55]) which accounts for bending, shear, torsion, and membrane effects within a small displacements and rotations framework. It should be noted that since this work focuses on the interface coupling aspects of fluid-beam system presented in Section 2.3 below, more complex beam models with non-linear constitutive laws can be considered (see, e.g., [19, 56, 20, 57, 58, 59, 60, 61]). In particular, a non-linear co-rotating beam formulation (see [62]) is considered in the numerical results of Section 5.3.

Let $(e_x, e_y, e_z)$ be a given right-handed Cartesian frame. We consider a reference straight beam of length $L$ and elliptical cross-sections orthogonal to its centerline, denoted by $\hat{\omega}_\varepsilon \stackrel{\text{def}}{=} \varepsilon\mathcal{E} \times (0, L)$, with $\varepsilon > 0$. Here, $\mathcal{E}$ stands for an elliptical region in the $(xOy)$ plane (major axis of length $a$ and a minor axis of length $1$) and $\varepsilon$ is a parameter related to the thickness of the beam. For $s \in (0, L)$, we denote by $\hat{\omega}_\varepsilon(s) \stackrel{\text{def}}{=} \varepsilon\mathcal{E} \times \{s\}$ the cross-section centered in $se_z$, such that $\hat{\omega}_\varepsilon = \bigcup_{s \in (0, L)} \hat{\omega}_\varepsilon(s)$ (see Figure 1). Note that, for the sake of simplicity, a straight beam with tangent vector $e_z$ is chosen as reference configuration of the beam. Nevertheless, a more general setting can be considered with the centerline given as a general curve with arc-length parametrization $r_0: [0, L] \to \mathbb{R}^3$ and unit tangent vector $t \stackrel{\text{def}}{=} r'_0$. We assume
<figure><img src="image_2.png" /><figcaption>(a) Reference configuration of the beam.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Elliptical cross-section of the beam in the reference configuration.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: The deformation map.</figcaption></figure>

that $\hat{\omega}_{\varepsilon}$ corresponds to the position of a beam at the initial time. The motion of the beam is then given in terms of the map $\phi: \hat{\omega}_{\varepsilon} \times \mathbb{R}^{+} \to \mathbb{R}^{3}$, so that the current configuration of the beam $\omega_{\varepsilon}(t)$ is given by $\omega_{\varepsilon}(t) = \phi(\hat{\omega}_{\varepsilon}, t)$ (see Figure 2). As discussed above, in the Timoshenko beam model, the cross-sections are supposed to remain undeformable but may not necessarily remain perpendicular to the beam centerline. The motion of the beam can hence be parametrized in the following way:

$$ \phi(\hat{x}, t) = se_z + \mathbf{r}(s, t) + \mathbf{\Lambda}(s, t)(\hat{x} - se_z), \quad s = \hat{x} \cdot e_z, \quad (\hat{x}, t) \in \hat{\omega}_{\varepsilon} \times \mathbb{R}^{+}, \qquad (2.1) $$

where $\mathbf{r}: (0, L) \times \mathbb{R}^{+} \to \mathbb{R}^{3}$ stands for the displacement of the centerline and $\mathbf{\Lambda}: (0, L) \times \mathbb{R}^{+} \to \mathbf{SO}(3)$ describes the rotation of the cross-sections. Here, $\mathbf{SO}(3)$ stands for the Special Orthogonal group given by

$$ \mathrm{SO}(3) \stackrel{\mathrm{def}}{=} \{ \mathbf{\Lambda} \in \mathbb{R}^{3 \times 3} \mid \mathbf{\Lambda}^{\mathrm{T}} \mathbf{\Lambda} = \mathbf{I}_3, \det \mathbf{\Lambda} = 1 \}, \qquad (2.2) $$

where $\mathbf{I}_3$ denotes the identity matrix in $\mathbb{R}^{3 \times 3}$. Note, in particular, that the centerline $\mathcal{B}(t)$ and the cross-sections $(\omega_{\varepsilon}(s, t))_{s \in (0, L)}$ of the current configuration of the beam are given respectively by

$$ \mathcal{B}(t) \stackrel{\mathrm{def}}{=} \phi(\{0\} \times \{0\} \times (0, L), t), \quad \omega_{\varepsilon}(s, t) \stackrel{\mathrm{def}}{=} \phi(\hat{\omega}_{\varepsilon}(s), t). $$

Due to the orthogonality of the rotation matrix $\mathbf{\Lambda}$, it can be given only in terms of three independent parameters. One option is to consider the rotational vector, defined by

$$ \theta \stackrel{\mathrm{def}}{=} \theta n $$
where **n** is a unit vector defining the axis of rotation and $\theta = \sqrt{\theta \cdot \theta}$ is the angle of rotation. We set $\Theta$ the skew symetric matrix associated to the vector $\theta$ such that for all $v \in \mathbb{R}^3$,

$$ \Theta v = \theta \wedge v, $$

where the symbol $\wedge$ denotes the cross product in three dimensions. Using Rodrigues' rotation formula (see, e.g., [63]), the relation between the rotation matrix and the rotational vector is then given by

$$ \mathbf{\Lambda} = \mathbf{I}_3 + \frac{\sin \theta}{\theta} \mathbf{\Theta} + \frac{1 - \cos \theta}{\theta^2} \mathbf{\Theta}^2. \qquad (2.3) $$

The dynamics of the beam are described, in terms of the centerline displacement **r** and the rotation vector **θ**, by the following system of equations:

$$ \left\{ \begin{array}{ll} \rho^b \mathbf{A} \hat{\partial}_t^2 \mathbf{r} - \mathbf{G} \hat{\partial}_s (\hat{\partial}_s \mathbf{r} - \mathbf{\theta} \wedge \mathbf{e}_z) = \mathbf{f} & \text{in } (0, L), \\ \rho^b \mathbf{I} \hat{\partial}_t^2 \mathbf{\theta} - \mathbf{E} \hat{\partial}_s^2 \mathbf{\theta} - \mathbf{e}_z \wedge \mathbf{G} (\hat{\partial}_s \mathbf{r} - \mathbf{\theta} \wedge \mathbf{e}_z) = \mathbf{m} & \text{in } (0, L), \end{array} \right. \qquad (2.4) $$

where the matrices **A**, **I**, **E** and **G** are given by

$$ \mathbf{A} \stackrel{\text{def}}{=} \begin{bmatrix} A & 0 & 0 \\ 0 & A & 0 \\ 0 & 0 & A \end{bmatrix}, \quad \mathbf{I} \stackrel{\text{def}}{=} \begin{bmatrix} I_x & 0 & 0 \\ 0 & I_y & 0 \\ 0 & 0 & J \end{bmatrix}, \quad \mathbf{E} \stackrel{\text{def}}{=} \begin{bmatrix} E I_x & 0 & 0 \\ 0 & E I_y & 0 \\ 0 & 0 & GJ \end{bmatrix}, \quad \mathbf{G} \stackrel{\text{def}}{=} \begin{bmatrix} G A \kappa & 0 & 0 \\ 0 & G A \kappa & 0 \\ 0 & 0 & EA \end{bmatrix}, $$

while the symbols $\rho^b$, $E$, $G$, $I_x$, $I_y$, $J$, $A$ and $\kappa$, respectively, denote the linear density, Young's modulus, the shear modulus, the principal moments of inertia in $x$ and $y$, the second moment of inertia, the cross-sectional area and the shear correction factor. Furthermore, the symbols **f** and **m** stand for the linear densities of load and torque applied to the beam, respectively. In this paper we shall make extensive use of an extension operator from the beam centerline $\{0\} \times \{0\} \times (0, L)$ to the beam reference domain $\hat{\omega}_\epsilon$, defined as follows: for any vector field $v : (0, L) \times \mathbb{R}^+ \to \mathbb{R}^3$ we define the vector field $\bar{v} : \hat{\omega}_\epsilon \times \mathbb{R}^+ \to \mathbb{R}^3$ by the relation

$$ \bar{v}(\hat{\mathbf{x}}, t) = v(s, t), \quad s = \hat{\mathbf{x}} \cdot \mathbf{e}_z, \quad \forall (\hat{\mathbf{x}}, t) \in \hat{\omega}_\epsilon \times \mathbb{R}^+. $$

In other words, functions defined on the centerline of the beam are lifted to the whole beam reference domain by using a projection operator which links each point of the beam to the cross-section to which it belongs. In what follows, we assume that **r** and **θ** are regular enough and that for all $t > 0$, $\phi_t \stackrel{\text{def}}{=} \phi(\cdot, t) : \hat{\omega}_\epsilon \to \omega_\epsilon(t)$ is one-to-one.

## 2.2. The fluid model

In this paper, the fluid is modeled by the incompressible Navier-Stokes equations. The beam is supposed to move within a fixed and bounded domain $\Omega \subset \mathbb{R}^3$ and we assume that it is surrounded by a fluid which occupies at each time $t$ the time-dependent domain $\Omega_\epsilon^f(t) \stackrel{\text{def}}{=} \Omega \setminus \bar{\omega}_\epsilon(t) \subset \mathbb{R}^3$. We denote by $\Sigma_\epsilon(t) \stackrel{\text{def}}{=} \partial\Omega_\epsilon^f(t) \cap \partial\omega_\epsilon(t)$ the current configuration of the fluid-solid interface. The fluid domain boundary is then partitioned as $\partial\Omega_\epsilon^f(t) = \Sigma_\epsilon(t) \cup \Gamma$ with $\Gamma \subset \partial\Omega$. Furthermore, we can introduce the trajectory of the fluid domain as the non-cylindrical space-time domain: $\mathcal{O} \stackrel{\text{def}}{=} \bigcup_{t \in \mathbb{R}^+} \Omega_\epsilon^f(t) \times \{t\}$. The fluid equations in the Eulerian formalism read: find the fluid velocity $\mathbf{u} : \mathcal{O} \to \mathbb{R}^3$ and fluid pressure $p : \mathcal{O} \to \mathbb{R}$ such that

$$ \left\{ \begin{array}{ll} \rho^f \partial_t \mathbf{u} + \rho^f \mathbf{u} \cdot \nabla \mathbf{u} - \nabla \cdot \boldsymbol{\sigma}(\mathbf{u}, p) = \mathbf{0} & \text{in } \Omega_\epsilon^f(t), \\ \nabla \cdot \mathbf{u} = 0 & \text{in } \Omega_\epsilon^f(t), \\ \mathbf{u} = \mathbf{0} & \text{on } \Gamma, \\ \mathbf{u} = \mathbf{u}^b & \text{on } \Sigma_\epsilon(t). \end{array} \right. \qquad (2.5) $$

Here, $\rho^f$ stands for the fluid density, $\mu^f$ for the dynamic viscosity of the fluid, and $\boldsymbol{\sigma}(\mathbf{u}, p)$ denotes the Cauchy stress tensor given by

$$ \boldsymbol{\sigma}(\mathbf{u}, p) \stackrel{\text{def}}{=} 2\mu^f \boldsymbol{\varepsilon}(\mathbf{u}) - p\mathbf{I}, \quad \boldsymbol{\varepsilon}(\mathbf{u}) \stackrel{\text{def}}{=} \frac{1}{2}(\nabla\mathbf{u} + (\nabla\mathbf{u})^\mathrm{T}). $$

Finally, $\mathbf{u}^b$ denotes the Eulerian velocity of the beam, which will be linked to the beam unknowns in the next section.
## 2.3. The full coupled problem

In this section, we discuss the interface coupling conditions between the beam model (2.4) and the fluid model (2.5). We then present the resulting full coupled system and an alternative model based on a fictitious domain formulation in the fluid.

**Interface coupling conditions.** The first coupling condition enforces the geometrical compatibility between the fluid and solid domains, namely,

$$ \phi = (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z + \bar{\mathbf{r}} + \bar{\Lambda}(\mathbf{I}_3 - (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z) \quad \text{in } \hat{\omega}_\varepsilon, \quad \omega_\varepsilon(t) = \phi(\hat{\omega}_\varepsilon, t), \quad \Omega_\varepsilon^f(t) \stackrel{\text{def}}{=} \Omega \setminus \bar{\omega}_\varepsilon(t). \quad (2.6) $$

The second condition corresponds to the kinematic coupling between the fluid and the beam on the interface $\Sigma_\varepsilon(t)$, it writes

$$ \mathbf{u}^b = \partial_t \phi \circ \phi_t^{-1} \quad \text{on } \Sigma_\varepsilon(t). \quad (2.7) $$

Finally, the third condition expresses the dynamic balance between the fluid and the beam on $\Sigma_\varepsilon(t)$, that is,

$$ \begin{align} (\mathbf{f}, \delta\mathbf{r})_{(0,L)} &= -(\boldsymbol{\sigma}(\mathbf{u}, p)\mathbf{n}, \delta\bar{\mathbf{r}} \circ \boldsymbol{\phi}_t^{-1})_{\Sigma_\varepsilon(t)}, \tag{2.8} \\ (\mathbf{m}, \delta\boldsymbol{\theta})_{(0,L)} &= -(\boldsymbol{\sigma}(\mathbf{u}, p)\mathbf{n}, \delta\bar{\boldsymbol{\theta}} \circ \boldsymbol{\phi}_t^{-1} \wedge (\mathbf{I}_3 - \bar{\mathbf{r}} \circ \boldsymbol{\phi}_t^{-1} - (\boldsymbol{\phi}_t^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\Sigma_\varepsilon(t)} \nonumber \end{align} $$

for all admissibles displacements and rotations $\delta\mathbf{r} : (0, L) \to \mathbb{R}^3$ and $\delta\boldsymbol{\theta} : (0, L) \to \mathbb{R}^3$.

**Full coupled problem in strong form.** By gathering the relations (2.4)-(2.8) we get the following full coupled problem: find the fluid velocity $\mathbf{u} : \mathcal{O} \to \mathbb{R}^3$, the fluid pressure $p : \mathcal{O} \to \mathbb{R}$, the beam displacement $\mathbf{r} : (0, L) \times \mathbb{R}^+ \to \mathbb{R}^3$ and the beam rotation $\boldsymbol{\theta} : (0, L) \times \mathbb{R}^+ \to \mathbb{R}^3$ such that

$$ \left\{ \begin{array}{ll} \rho^f \partial_t \mathbf{u} + \rho^f \mathbf{u} \cdot \nabla \mathbf{u} - \nabla \cdot \boldsymbol{\sigma}(\mathbf{u}, p) = \mathbf{0} & \text{in } \Omega_\varepsilon^f(t), \\ \nabla \cdot \mathbf{u} = \mathbf{0} & \text{in } \Omega_\varepsilon^f(t), \\ \mathbf{u} = \mathbf{0} & \text{on } \Gamma, \end{array} \right. \quad (2.9) $$

$$ \left\{ \begin{array}{ll} \rho^b \mathbf{A} \partial_t^2 \mathbf{r} - \mathbf{G} \partial_s (\partial_s \mathbf{r} - \boldsymbol{\theta} \wedge \mathbf{e}_z) = \mathbf{f} & \text{in } (0, L), \\ \rho^b \mathbf{I} \partial_t^2 \boldsymbol{\theta} - \mathbf{e} \partial_s^2 \boldsymbol{\theta} - \mathbf{e}_z \wedge \mathbf{G} (\partial_s \mathbf{r} - \boldsymbol{\theta} \wedge \mathbf{e}_z) = \mathbf{m} & \text{in } (0, L), \end{array} \right. \quad (2.10) $$

$$ \left\{ \begin{array}{ll} \phi = (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z + \bar{\mathbf{r}} + \bar{\Lambda}(\mathbf{I}_3 - (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z) & \text{in } \hat{\omega}_\varepsilon, \quad \omega_\varepsilon(t) = \phi(\hat{\omega}_\varepsilon, t), \quad \Omega_\varepsilon^f(t) \stackrel{\text{def}}{=} \Omega \setminus \bar{\omega}_\varepsilon(t), \\ \mathbf{u} = \partial_t \phi \circ \phi_t^{-1} & \text{on } \Sigma_\varepsilon(t), \\ (\mathbf{f}, \delta\mathbf{r})_{(0,L)} = -(\boldsymbol{\sigma}(\mathbf{u}, p)\mathbf{n}, \delta\bar{\mathbf{r}} \circ \phi_t^{-1})_{\Sigma_\varepsilon(t)} & \forall \delta\mathbf{r}, \\ (\mathbf{m}, \delta\boldsymbol{\theta})_{(0,L)} = -(\boldsymbol{\sigma}(\mathbf{u}, p)\mathbf{n}, \delta\bar{\boldsymbol{\theta}} \circ \phi_t^{-1} \wedge (\mathbf{I}_3 - \bar{\mathbf{r}} \circ \phi_t^{-1} - (\phi_t^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\Sigma_\varepsilon(t)} & \forall \delta\boldsymbol{\theta}. \end{array} \right. \quad (2.11) $$

The coupled model (2.9)-(2.11) will be discretized with an ALE formalism in the fluid sub-problem (2.9) as presented in Algorithm 2, providing a reference solution to assess the precision of our reduced approach, which will be presented as Algorithm 1. In the case of moderate interface deflections, the ALE formalism in the fluid sub-problem is useful as it facilitates time and spatial approximation using body fitted meshes. However, whenever the beam undergoes arbitrarily large displacements, the ALE formalism becomes cumbersome. In this case, a favored approach zconsists in considering a fictitious domain formulation.

**Full coupled problem with fictitious domain formulation in the fluid.** In the fictitious domain formulation, the fluid equations (2.9) are solved in the entire domain $\Omega$ and the kinematic constraint (2.11)$_2$ on $\Sigma_\varepsilon(t)$ is enforced in a weak fashion through Lagrange multipliers, viz.,

$$ \langle \boldsymbol{\mu}, \mathbf{u} \rangle_{\Sigma_\varepsilon(t)} = \langle \boldsymbol{\mu}, \partial_t \phi \circ \phi_t^{-1} \rangle_{\Sigma_\varepsilon(t)}, \quad \forall \boldsymbol{\mu} \in \mathbf{H}^{-\frac{1}{2}}(\Sigma_\varepsilon(t)), \quad (2.12) $$

where $\langle \cdot, \cdot \rangle_{\Sigma_\varepsilon(t)}$ denotes the duality pairing between $\mathbf{H}^{-\frac{1}{2}}(\Sigma_\varepsilon(t))$ and $\mathbf{H}^{\frac{1}{2}}(\Sigma_\varepsilon(t))$.

**Remark 2.1.** Note that, for the sake of simplicity, the weak equality of the velocities (2.12) assumed that the beam is fully embedded in the fluid and thus $\Sigma_\varepsilon(t) \cap \partial\Omega = \emptyset$. In the case of a beam attached to the fluid cavity then the duality pairing should involve the space $\mathbf{H}_{00}^{\frac{1}{2}}(\Sigma_\varepsilon(t))$ (the trace space of vector valued functions whose extensions by zero on $\Gamma$ belong to $\mathbf{H}^{\frac{1}{2}}(\partial\Omega_\varepsilon(t))$) and its dual.
Let $V \stackrel{\text{def}}{=} H_{\partial\Omega}^1(\Omega)$ and $Q \stackrel{\text{def}}{=} L_0^2(\Omega)$ be the functional spaces for the fluid velocity and pressure, respectively, and $Y$ the space of admissible beam displacements and rotations. The resulting coupled problem with a fictitious domain formulation in the fluid then reads: find $(\mathbf{u}, p, \boldsymbol{\lambda}, \mathbf{r}, \boldsymbol{\theta})$ such that

$$ (\mathbf{u}(t), p(t), \boldsymbol{\lambda}(t), \mathbf{r}(t), \boldsymbol{\theta}(t)) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{H}^{-\frac{1}{2}}(\Sigma_{\varepsilon}(t)) \times \mathbf{Y} \times \mathbf{Y} $$

a.e. $t \in \mathbb{R}^+$ and

$$ \left\{ \begin{array}{l} \phi = (I_3 \cdot e_z)e_z + \bar{r} + \bar{\Lambda}(I_3 - (I_3 e_z)e_z) \quad \text{in } \hat{\omega}_{\varepsilon}, \quad \Sigma_{\varepsilon}(t) = \phi(\hat{\Sigma}_{\varepsilon}, t), \\ \rho^f(\partial_t \mathbf{u}, \mathbf{v})_{\Omega} + a_{\Omega}^f(\mathbf{u}; (\mathbf{u}, p), (\mathbf{v}, q)) \\ \qquad + \rho^b(\mathbf{A}\partial_t^2 \mathbf{r}, \delta\mathbf{r})_{(0,L)} + \rho^b(\mathbf{I}\partial_t^2 \boldsymbol{\theta}, \delta\boldsymbol{\theta})_{(0,L)} + a^b((\mathbf{r}, \boldsymbol{\theta}), (\delta\mathbf{r}, \delta\boldsymbol{\theta})) \\ \qquad - \langle \boldsymbol{\lambda}, \mathbf{v} - \bar{\delta}\mathbf{r} \circ \phi_t^{-1} - \bar{\delta}\boldsymbol{\theta} \circ \phi_t^{-1} \wedge (I_3 - \bar{r} \circ \phi_t^{-1} - (\phi_t^{-1} \cdot e_z)e_z) \rangle_{\Sigma_{\varepsilon}(t)} \\ \qquad + \langle \boldsymbol{\mu}, \mathbf{u} - \partial_t \boldsymbol{\phi} \circ \phi_t^{-1} \rangle_{\Sigma_{\varepsilon}(t)} = 0 \end{array} \right. \quad (2.13) $$

for all $(\mathbf{v}, q, \boldsymbol{\mu}, \delta\mathbf{r}, \delta\boldsymbol{\theta}) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{H}^{-\frac{1}{2}}(\Sigma_{\varepsilon}(t)) \times \mathbf{Y} \times \mathbf{Y}$, with the notations

$$ \begin{aligned} a_{\Omega}^{f}(\mathbf{z}; (\mathbf{u}, p), (\mathbf{v}, q)) &\stackrel{\text{def}}{=} c_{\Omega}(\mathbf{z}; (\mathbf{u}, \mathbf{v})) + a_{\Omega}((\mathbf{u}, p), (\mathbf{v}, q)), & c_{\Omega}(\mathbf{z}; (\mathbf{u}, \mathbf{v})) &\stackrel{\text{def}}{=} \rho^{f}(\mathbf{z} \cdot \nabla \mathbf{u}, \mathbf{v})_{\Omega}, \\ a_{\Omega}((\mathbf{u}, p), (\mathbf{v}, q)) &\stackrel{\text{def}}{=} 2\mu^{f}(\boldsymbol{\varepsilon}(\mathbf{u}), \boldsymbol{\varepsilon}(\mathbf{v}))_{\Omega} - (p, \nabla \cdot \mathbf{v})_{\Omega} + (q, \nabla \cdot \mathbf{u})_{\Omega}, \\ a^{b}((\mathbf{r}, \boldsymbol{\theta}), (\delta\mathbf{r}, \delta\boldsymbol{\theta})) &\stackrel{\text{def}}{=} (\mathbf{E}\partial_{s}\boldsymbol{\theta}, \partial_{s}\delta\boldsymbol{\theta})_{(0,L)} + (\mathbf{G}(\partial_{s}\mathbf{r} - \boldsymbol{\theta} \wedge e_z), (\partial_{s}\delta\mathbf{r} - \delta\boldsymbol{\theta} \wedge e_z))_{(0,L)}. \end{aligned} $$

We conclude this section with a couple of remarks. First, it should be noted that (2.13) induces a modeling error with respect to the coupled problem (2.9)-(2.11). Indeed, since the fluid equations are solved in the entire computational domain $\Omega$, the Lagrange multiplier $\boldsymbol{\lambda}$ corresponds to the stress jump across the interface $\Sigma_{\varepsilon}(t)$. The numerical experiments of Section 5 indicate however that this modeling error is negligible for small values of the thickness parameter $\varepsilon$. Finally, though mathematically sound, the numerical approximation of the coupled model (2.13) poses an important challenge from the computational point of view. Indeed, since the interface conditions are enforced on the physical boundary of the beam $\Sigma_{\varepsilon}(t)$, the fluid mesh needs to be sufficiently fine to correctly resolve this interface, which becomes computationally demanding whenever $\varepsilon$ is small.

### 3. Coupled problem with reduced interface conditions

In this section we present two reduced modeling approaches which mitigate the computational complexity of the full coupled model (2.13) by leveraging the slender nature of the beam. To this purpose, we first briefly review a modeling approach from the literature, which is known of coming with some fundamental drawbacks (Section 3.1). Then we introduce a novel modeling approach which circumvents these issues (Section 3.2).

#### 3.1. Fictitious domain approach with kinematic constraint on the beam centerline

In [33, 34] (see also [31, 32, 8]), the complexity of (2.13) is reduced by assuming that the thickness of the beam can be neglected in the interface coupling. In other words, the interface conditions are enforced directly on the beam centerline $\mathcal{B}(t)$, rather than on the physical interface $\Sigma_{\varepsilon}(t)$. The kinematic constraint, in particular, becomes

$$ (\boldsymbol{\mu}, \mathbf{u})_{\mathcal{B}(t)} = (\boldsymbol{\mu}, \partial_t \mathbf{r} \circ \phi_t^{-1})_{\mathcal{B}(t)} \quad \forall \boldsymbol{\mu} \in \mathbf{L}^2(\mathcal{B}(t)). \quad (3.1) $$

The resulting reduced coupled problem then writes: find $(\mathbf{u}, p, \boldsymbol{\lambda}, \mathbf{r}, \boldsymbol{\theta})$ such that

$$ (\mathbf{u}(t), p(t), \boldsymbol{\lambda}(t), \mathbf{r}(t), \boldsymbol{\theta}(t)) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{L}^2(\mathcal{B}(t)) \times \mathbf{Y} \times \mathbf{Y}, $$

a.e. $t \in \mathbb{R}^+$ and

$$ \left\{ \begin{array}{l} \phi = (I_3 \cdot e_z)e_z + \bar{r} \quad \text{on } \{0\} \times \{0\} \times (0, L), \quad \mathcal{B}(t) = \phi(\{0\} \times \{0\} \times (0, L), t), \\ \rho^f(\partial_t \mathbf{u}, \mathbf{v})_{\Omega} + a_{\Omega}^f(\mathbf{u}; (\mathbf{u}, p), (\mathbf{v}, q)) \\ \qquad + \rho^b(\mathbf{A}\partial_t^2 \mathbf{r}, \delta\mathbf{r})_{(0,L)} + \rho^b(\mathbf{I}\partial_t^2 \boldsymbol{\theta}, \delta\boldsymbol{\theta})_{(0,L)} + a^b((\mathbf{r}, \boldsymbol{\theta}), (\delta\mathbf{r}, \delta\boldsymbol{\theta})) \\ \qquad - (\boldsymbol{\lambda}, \mathbf{v} - \delta\mathbf{r} \circ \phi_t^{-1})_{\mathcal{B}(t)} + (\boldsymbol{\mu}, \mathbf{u} - \partial_t \mathbf{r} \circ \phi_t^{-1})_{\mathcal{B}(t)} = 0 \end{array} \right. \quad (3.2) $$
for all $(v, q, \mu, \delta r, \delta \theta) \in V \times Q \times L^2(B(t)) \times Y \times Y$. The coupled model (3.2) will be discretized as described in Algorithm 3, providing an element for comparison in terms of accuracy and computational efficiency for our reduced approach.

Clearly, the kinematic condition (3.1) reduces the complexity of the coupled problem (3.2) with respect to (2.12), by passing from a Lagrange multiplier defined on a variety of co-dimension one, $\Sigma_{\epsilon}(t)$, to a Lagrange multiplier defined on a variety of co-dimension two, $B(t)$. This reduced computational complexity comes however at a price. First, the rotational effects are neglected in the fluid-beam coupling: no rotational velocity is defined on the beam centerline, so that the fluid torque acting on the beam is null. And second, mathematically, the relation (3.1) is only well-defined for $u \in H^{1+\eta}(\Omega)$, with $\eta > 0$. In the context of problem (3.2), this regularity for $u$ might not be attained, notably due to the singularity induced by the term $(\lambda(t), v)_{B(t)}$. This prevents from establishing both the well-posedness of problem (3.2) and the accuracy of the associated numerical methods, which is a major drawback of this reduced formulation.

## 3.2. A new fictitious domain approach with reduced order kinematic conditions

The basic idea of the new reduced modeling approach introduced in this paper is twofold. First, the interface coupling conditions are enforced on the physical interface $\Sigma_{\epsilon}(t)$. This circumvents the two shortcomings discussed in the previous section. Second, the interface conditions are enforced in a reduced fashion, by replacing the infinite dimensional Lagrange multiplier space $H^{-\frac{1}{2}}(\Sigma_{\epsilon}(t))$ by a low dimensional space $F_N(t)$, which hence reduces the computational complexity of the resulting coupled problem.

In order to define the space $F_N(t)$, we first introduce a local coordinate $\nu$ on each cross-section of the reference configuration of the fluid-solid interface. Any point $\hat{x}$ on $\hat{\Sigma}_{\epsilon} \stackrel{\text{def}}{=} \partial\hat{\Omega}_{\epsilon}^{f} \cap \partial\hat{\omega}_{\epsilon}$ can be parametrized as

$$ \hat{x} = se_z + \epsilon(\cos(\nu)e_x + a\sin(\nu)e_y), \quad \nu \in (0, 2\pi). \qquad (3.3) $$

In particular, any function $v$ defined on $\hat{\Sigma}_{\epsilon}$ can be expressed as follows:

$$ v(s, \nu) \stackrel{\text{def}}{=} v(se_z + \epsilon(\cos(\nu)e_x + a\sin(\nu)e_y)) \quad \forall (s, \nu) \in (0, L) \times (0, 2\pi). $$

We then set for all $N \in \mathbb{N}$ and $t \in \mathbb{R}^{+}$, the space $F_N(t)$ defined by

$$ F_N(t) \stackrel{\text{def}}{=} \left\{ v \circ \phi_t^{-1} \mid v(s, \nu) = a_0(s) + \sum_{k=1}^{N} (a_k(s) \cos(k\nu) + b_k(s) \sin(k\nu)), \ a_k, b_k \in L^2(0, L) \right\}. \quad (3.4) $$

The idea of (3.4) is to take advantage of the tensorization of $\hat{\Sigma}_{\epsilon}$ by $(0, L) \times (0, 2\pi)$ to discretize the space $L^2(0, 2\pi)$ using the finite dimensional Fourier space of degree $N$ (see also [52, 50] for a similar approach in the context of elliptic problems). This reduced representation is then transported to the current configuration of the fluid-solid interface $\Sigma_{\epsilon}(t)$ via composition with $\phi_t^{-1}$. The expression reduced comes from the fact that any function $v \in F_N(t)$ is now uniquely determined by the set of functions $\{a_k\}_{0 \le k \le N}, \{b_k\}_{1 \le k \le N} \subset L^2(0, L)$, thus passing from a functional space $H^{-\frac{1}{2}}(\Sigma_{\epsilon}(t))$, defined in a two-dimensional domain, to a set of $2N + 1$ vector valued functions defined on a 1D domain. In particular, the reduced kinematic coupling conditions in $\Sigma_{\epsilon}(t)$ is given by

$$ (\mu_N, u(t))_{\epsilon,t} = (\mu_N, \partial_t \phi \circ \phi_t^{-1})_{\epsilon,t} \quad \forall \mu_N \in F_N(t), \quad t \in \mathbb{R}^+, \qquad (3.5) $$

where $(\cdot, \cdot)_{\epsilon,t}$ denotes the re-scaled $L^2$ inner-product on $\Sigma_{\epsilon}(t)$:

$$ (u, v)_{\epsilon,t} \stackrel{\text{def}}{=} \int_{\Sigma_{\epsilon}(t)} (u \cdot v)h \circ \phi_t^{-1}, \quad h(s, \nu) \stackrel{\text{def}}{=} \left( \epsilon \sqrt{a^2 \cos^2(\nu) + \sin^2(\nu)} \right)^{-1} $$

for all $(s, \nu) \in (0, L) \times (0, 2\pi)$.

**Remark 3.1.** The choice of the scalar product $(\cdot, \cdot)_{\epsilon,t}$ is justified by the error approximation between the solution of the reduced problem and the full problem, where we consider $H^{-\frac{1}{2}}(\Sigma_{\epsilon}(t))$ as the Lagrange multiplier space. Specifically, in [64] numerical evidence is given on that, for a 2D simplified problem and for $F_N(t)$ chosen as in (3.4), the order of convergence with respect to $\epsilon$ is improved compared to the standard scalar product $(\cdot, \cdot)_{\Sigma_{\epsilon}(t)}$.
The associated reduced coupled model reads as follows: find $(\mathbf{u}, p, \lambda_N, \mathbf{r}, \theta)$ such that

$$ (\mathbf{u}(t), p(t), \lambda_N(t), \mathbf{r}(t), \theta(t)) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{F}_N(t) \times \mathbf{Y} \times \mathbf{Y} $$

a.e. $t \in \mathbb{R}^+$ and

$$ \left\{ \begin{aligned} \phi &= (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z + \bar{\mathbf{r}} + \bar{\mathbf{\Lambda}}(\mathbf{I}_3 - (\mathbf{I}_3 \mathbf{e}_z)\mathbf{e}_z) \quad \text{in } \hat{\omega}_\varepsilon, \quad \Sigma_\varepsilon(t) = \phi(\hat{\Sigma}_\varepsilon, t), \\ \rho^f(\partial_t \mathbf{u}, \mathbf{v})_\Omega + a_\Omega^f(\mathbf{u}; (\mathbf{u}, p), (\mathbf{v}, q)) \\ &\quad + \rho^b(\mathbf{A}\partial_t^2 \mathbf{r}, \delta\mathbf{r})_{(0,L)} + \rho^b(\mathbf{I}\partial_t^2 \theta, \delta\theta)_{(0,L)} + a^b((\mathbf{r}, \theta), (\delta\mathbf{r}, \delta\theta)) \\ &\quad - (\lambda_N, \mathbf{v} - \bar{\delta}\mathbf{r} \circ \phi_t^{-1} - \bar{\delta}\theta \circ \phi_t^{-1} \wedge (\mathbf{I}_3 - \bar{\mathbf{r}} \circ \phi_t^{-1} - (\phi_t^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\varepsilon,t} \\ &\quad + (\mu_N, \mathbf{u} - \partial_t \phi \circ \phi_t^{-1})_{\varepsilon,t} = 0 \end{aligned} \right. \quad (3.6) $$

for all $(\mathbf{v}, q, \mu_N, \delta\mathbf{r}, \delta\theta) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{F}_N(t) \times \mathbf{Y} \times \mathbf{Y}$.

*Energy estimate (small rotational velocity).* We conclude this section by deriving an energy estimate for the reduced coupled model (3.6) valid in the case in which the beam undergoes small rotational velocities. This assumption is needed to cope with the mismatch between the treatment of the geometrical non-linearities in the solid (linear) and in the fluid (non-linear). We hence consider a variant of the kinematic coupling condition (3.5) in which the Eulerian velocity of the beam interface, $\partial_t \phi \circ \phi_t^{-1}$, is taken with a linearized angular velocity, while keeping the non linear reconstruction of the geometry. More precisely, denoting by $\mathbf{w}$ is the instantaneous vector velocity associated to the rotation vector $\theta$, we assume

$$ \mathbf{w} \approx \partial_t \theta, \quad (3.7) $$

which yields

$$ \partial_t \phi \circ \phi_t^{-1} \approx \bar{\partial}_t \mathbf{r} \circ \phi_t^{-1} + \bar{\partial}_t \theta \circ \phi_t^{-1} \wedge (\mathbf{I}_3 - \bar{\mathbf{r}} \circ \phi_t^{-1} - (\phi_t^{-1} \cdot \mathbf{e}_z) \mathbf{e}_z). \quad (3.8) $$

It should be noted that the relations (3.7) and (3.8) are only valid in a small rotational velocity framework. The resulting coupled problem in weak form read as follows: we look for $(\mathbf{u}, p, \lambda_N, \mathbf{r}, \theta)$ such that

$$ (\mathbf{u}(t), p(t), \lambda_N(t), \mathbf{r}(t), \theta(t)) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{F}_N(t) \times \mathbf{Y} \times \mathbf{Y}, $$

a.e. $t \in \mathbb{R}^+$ and

$$ \left\{ \begin{aligned} \phi &= (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z + \bar{\mathbf{r}} + \bar{\mathbf{\Lambda}}(\mathbf{I}_3 - (\mathbf{I}_3 \mathbf{e}_z)\mathbf{e}_z) \quad \text{in } \hat{\omega}_\varepsilon, \quad \Sigma_\varepsilon(t) = \phi(\hat{\Sigma}_\varepsilon, t), \\ \rho^f(\partial_t \mathbf{u}, \mathbf{v})_\Omega + a_\Omega^f(\mathbf{u}; (\mathbf{u}, p), (\mathbf{v}, q)) \\ &\quad + \rho^b(\mathbf{A}\partial_t^2 \mathbf{r}, \delta\mathbf{r})_{(0,L)} + \rho^b(\mathbf{I}\partial_t^2 \theta, \delta\theta)_{(0,L)} + a^b((\mathbf{r}, \theta), (\delta\mathbf{r}, \delta\theta)) \\ &\quad - (\lambda_N, \mathbf{v} - \bar{\delta}\mathbf{r} \circ \phi_t^{-1} - \bar{\delta}\theta \circ \phi_t^{-1} \wedge (\mathbf{I}_3 - \bar{\mathbf{r}} \circ \phi_t^{-1} - (\phi_t^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\varepsilon,t} \\ &\quad + (\mu_N, \mathbf{u} - \partial_t \bar{\mathbf{r}} \circ \phi_t^{-1} - \partial_t \bar{\theta} \circ \phi_t^{-1} \wedge (\mathbf{I}_3 - \bar{\mathbf{r}} \circ \phi_t^{-1} - (\phi_t^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\varepsilon,t} = 0 \end{aligned} \right. \quad (3.9) $$

for all $(\mathbf{v}, q, \mu_N, \delta\mathbf{r}, \delta\theta) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{F}_N(t) \times \mathbf{Y} \times \mathbf{Y}$. The energy stability of (3.9) is stated in the next result.

**Theorem 3.2.** Let $(\mathbf{u}, p, \lambda_N, \mathbf{r}, \theta)$ be a regular enough solution of the reduced coupled problem (3.9). There holds

$$ \mathfrak{E}(t) \le \mathfrak{E}(0) \quad \forall t \in \mathbb{R}^+, \quad (3.10) $$

where the mechanical energy of the system $\mathfrak{E}(t)$ is defined by

$$ \mathfrak{E}(t) \stackrel{\text{def}}{=} \frac{\rho^f}{2} \|\mathbf{u}(t)\|_{0,\Omega}^2 + \frac{\rho^b}{2} \|\partial_t \mathbf{r}(t)\|_{\mathbf{A},(0,L)}^2 + \frac{\rho^b}{2} \|\partial_t \theta(t)\|_{\mathbf{I},(0,L)}^2 + \frac{1}{2} \|(\mathbf{r}(t), \theta(t))\|_{\mathbf{b}}^2, $$

with

$$ \|\cdot\|_{0,\Omega} \stackrel{\text{def}}{=} \sqrt{(\cdot, \cdot)_\Omega}, \quad \|\cdot\|_{\mathbf{A},(0,L)} \stackrel{\text{def}}{=} \sqrt{(\mathbf{A}, \cdot)_\Omega}, \quad \|\cdot\|_{\mathbf{I},(0,L)} \stackrel{\text{def}}{=} \sqrt{(\mathbf{I}, \cdot)_\Omega}, \quad \|\cdot\|_{\mathbf{b}} \stackrel{\text{def}}{=} \sqrt{a^b((\cdot, \cdot), (\cdot, \cdot))}. $$
Proof. By testing (3.6) with $(v, q, \delta r, \delta \theta, \mu_N) = (u(t), p(t), \partial_t r(t), \partial_t \theta(t), \lambda_N(t))$ we get

$$ \rho^f(\partial_t u, u_N)_\Omega + \rho^b(\mathbf{A}\partial_t^2 r, \partial_t r)_{(0,L)} + \rho^b(\mathbf{I}\partial_t^2 \theta, \partial_t \theta)_{(0,L)} \\ + a_\Omega^f(u; (u, p), (u, p)) + a^b((r, \theta), (\partial_t r, \partial_t \theta)) = 0. \quad (3.11) $$

Using integration by parts and the boundary conditions of $u$, we have

$$ a_{\Omega}^{f}(u; (u, p), (u, p)) = 2\mu \|\varepsilon(u)\|_{0,\Omega}^{2}. $$

Since the remaining terms of (3.11) are inner-products, we finally get

$$ \frac{d}{dt} \left( \frac{\rho^f}{2} \|\mathbf{u}\|_{\Omega}^{2} + \frac{\rho^b}{2} \|\partial_t \mathbf{r}\|_{\mathbf{A},(0,L)}^{2} + \frac{\rho^b}{2} \|\partial_t \boldsymbol{\theta}\|_{\mathbf{I},(0,L)}^{2} + \frac{1}{2} \|(\mathbf{r}, \boldsymbol{\theta})\|_{\mathbf{b}}^{2} \right) \le 0. $$

The estimate (3.10) then follows by integrating this bound over $(0, t)$, which completes the proof. $\square$

**Remark 3.3.** The linearization given by (3.7)-(3.8) for the velocity, as assumed in Theorem 3.2, is necessary due to the linear nature of equations (2.4). To get rid of this assumption, we may use a non-linear beam model considering the exact rotation of the beam cross-sections, as for instance the Simo-Reisner geometrically exact beam model [19] or the corotational formulation described in [62, Chapter 5].

**Remark 3.4.** For the full coupled problem (2.9)-(2.11), the linearization of the instantaneous vector velocity $w$ prevents from establishing the energy estimate (3.10). In order to guarantee energy stability for this model, one simple option is to update the fluid domain from a linearized version of (2.1) with $\mathbf{\Lambda} \approx \mathbf{I}_3 + \mathbf{\Theta}$. Another option would be, as discussed in the previous remark, to directly consider a non-linear model of the beam.

## 4. Numerical methods

This section is devoted to the numerical approximation of the reduced coupled problem (3.6). The proposed numerical method is detailed in Section 4.1. Some alternative discrete formulations are briefly discussed in Section 4.2 for comparison purposes. In what follows, we denote by $\tau > 0$ the time-step length and, for $n \in \mathbb{N}$, we set $t^n \stackrel{\text{def}}{=} n\tau$. For any time-dependent function $v$, $v^n$ stands for an approximation of $v(t^n)$. We also introduce the standard notations

$$ \partial_{\tau} \mathbf{v}^{n} \stackrel{\text{def}}{=} \frac{1}{\tau}(\mathbf{v}^{n} - \mathbf{v}^{n-1}), \quad \mathbf{v}^{n-\frac{1}{2}} \stackrel{\text{def}}{=} \frac{1}{2}(\mathbf{v}^{n-1} + \mathbf{v}^{n}), $$

for the first-order backward difference and the mid-point value, respectively. Let $\mathcal{D}$ be an open polygonal convex domain of $\mathbb{R}^3$ and $\mathcal{T}_h(\mathcal{D})$ a triangulation of $\mathcal{D}$ with characteristic size $h$. We consider the following standard continuous Lagrange finite element spaces of degree $k \ge 1$:

$$ \begin{aligned} \mathbf{X}_h^k(\mathcal{D}) &\stackrel{\text{def}}{=} \{v \in \mathbf{H}^1(\mathcal{D}) \mid v|_{\tau} \in [\mathbb{P}_k(K)]^3, \forall K \in \mathcal{T}_h(\mathcal{D})\}, \\ \mathcal{Q}_h^k(\mathcal{D}) &\stackrel{\text{def}}{=} \{v \in \mathbf{H}^1(\mathcal{D}) \mid v|_{\tau} \in \mathbb{P}_k(K), \forall K \in \mathcal{T}_h(\mathcal{D})\}. \end{aligned} $$

We then introduce the finite element spaces $\mathbf{V}_h \stackrel{\text{def}}{=} \mathbf{X}_h^1(\Omega) \cap \mathbf{V}$ and $\mathcal{Q}_h \stackrel{\text{def}}{=} \mathcal{Q}_h^1(\Omega) \cap \mathcal{Q}$ for the approximation of the fluid velocity and pressure, respectively. The approximations of the solid displacement and rotation will be looked for into $\mathbf{Y}_\mathcal{H} \stackrel{\text{def}}{=} \mathbf{X}_\mathcal{H}^1(0, L) \cap \mathbf{Y}$. Here, $h$ and $\mathcal{H}$ respectively denote the characteristic sizes of the fluid and solid (centerline) meshes.

### 4.1. Discrete formulation of the coupled model with reduced interface conditions

For the time discretization of (3.6) we consider a semi-implicit scheme in which the geometrical coupling is treated in an explicit manner, as follows:

$$ \phi_{\mathcal{H}}^{n} = (\mathbf{I}_3 \cdot e_z)e_z + \overline{\mathbf{r}_{\mathcal{H}}^{n-1}} + \overline{\mathbf{\Lambda}_{\mathcal{H}}^{n-1}}(\mathbf{I}_3 - (\mathbf{I}_3 \cdot e_z)e_z) \quad \text{in } \hat{\omega}_{\varepsilon}, \quad \Sigma_{\varepsilon}^{n} = \phi_{\mathcal{H}}^{n}(\hat{\Sigma}_{\varepsilon}). $$
We can hence introduce the fully discrete reduced Lagrange multiplier space $F_{N,\mathcal{H}}^n$, discrete counterpart of $F_N(t_n)$, as

$$ F_{N,\mathcal{H}}^n \stackrel{\text{def}}{=} \left\{ \mathbf{v}_{\mathcal{H}} \circ (\phi_{\mathcal{H}}^n)^{-1} \middle| \mathbf{v}_{\mathcal{H}}(s, \nu) = \mathbf{a}_{0,\mathcal{H}}(s) + \sum_{k=1}^{N} (\mathbf{a}_{k,\mathcal{H}}(s) \cos(k\nu) + \mathbf{b}_{k,\mathcal{H}}(s) \sin(k\nu)), \mathbf{a}_{k,\mathcal{H}}, \mathbf{b}_{k,\mathcal{H}} \in \mathbf{Y}_{\mathcal{H}} \right\}, \quad (4.1) $$

which simply amounts to replace $L^2(0, L)$ by $\mathbf{Y}_{\mathcal{H}}$ in (3.4). Note that here a unique approximation space $\mathbf{Y}_{\mathcal{H}}$ is involved in the approximation of the beam unknowns and of the Lagrange multiplier, but this is not mandatory. We consider a backward-Euler semi-implicit time-discretization for the fluid and a mid-point scheme for the solid. The remaining kinematic and dynamic conditions in (3.6) are discretized with an implicit scheme. The resulting solution procedure is of strongly coupled nature (see, e.g., [65]). By gathering all the above mentioned ingredients, the proposed numerical approximation of (3.6) is detailed in Algorithm 1.

**Algorithm 1** Discrete formulation of the reduced coupled model (3.6).

For $n \ge 1$,

**Step 1:** Update interface position:

$$ \phi_{\mathcal{H}}^n = (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z + \bar{\mathbf{r}}_{\mathcal{H}}^{n-1} + \bar{\Lambda}_{\mathcal{H}}^{n-1} (\mathbf{I}_3 - (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z) \quad \text{in } \hat{\omega}_{\varepsilon}, \quad \Sigma_{\varepsilon}^n = \phi_{\mathcal{H}}^n(\hat{\Sigma}_{\varepsilon}), \quad (4.2) $$

where $\Lambda_{\mathcal{H}}^{n-1}$ is obtained from $\theta_{\mathcal{H}}^{n-1}$ using non linear formula 2.3.

**Step 2:** Find $(\mathbf{u}_h^n, p_h^n, \boldsymbol{\lambda}_{N,\mathcal{H}}, \mathbf{r}_h^n, \boldsymbol{\theta}_h^n) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_{N,\mathcal{H}}^n \times \mathbf{Y}_{\mathcal{H}} \times \mathbf{Y}_{\mathcal{H}}$ with $\dot{\mathbf{r}}_{\mathcal{H}}^{n-1/2} = \partial_{\tau} \mathbf{r}_{\mathcal{H}}^n$ and $\dot{\boldsymbol{\theta}}_{\mathcal{H}}^{n-1/2} = \partial_{\tau} \boldsymbol{\theta}_{\mathcal{H}}^n$ such that

$$ \left\{ \begin{aligned} & \rho^f(\partial_{\tau} \mathbf{u}_h^n, \mathbf{v}_h)_{\Omega} + a_{\Omega,h}^f(\mathbf{u}_h^{n-1}; (\mathbf{u}_h^n, p_h^n, (\mathbf{v}_h, q_h)) \\ & + \rho^b(\mathbf{A}\partial_{\tau}\dot{\mathbf{r}}_{\mathcal{H}}^n, \delta\mathbf{r}_{\mathcal{H}})_{(0,L)} + \rho^b(\mathbf{I}\partial_{\tau}\dot{\boldsymbol{\theta}}_{\mathcal{H}}^n, \delta\boldsymbol{\theta}_{\mathcal{H}})_{(0,L)} + a_{\mathcal{H}}^b((\mathbf{r}_{\mathcal{H}}^{n-1/2}, \boldsymbol{\theta}_{\mathcal{H}}^{n-1/2}), (\delta\mathbf{r}_{\mathcal{H}}, \delta\boldsymbol{\theta}_{\mathcal{H}})) \\ & - (\boldsymbol{\lambda}_{N,\mathcal{H}}, \mathbf{v}_h - \bar{\delta}\mathbf{r}_{\mathcal{H}} \circ (\boldsymbol{\phi}_{\mathcal{H}}^n)^{-1} - \bar{\delta}\boldsymbol{\theta}_{\mathcal{H}} \circ (\boldsymbol{\phi}_{\mathcal{H}}^n)^{-1} \wedge (\mathbf{I}_3 - \bar{\mathbf{r}}_{\mathcal{H}}^{n-1} \circ (\boldsymbol{\phi}_{\mathcal{H}}^n)^{-1} - ((\boldsymbol{\phi}_{\mathcal{H}}^n)^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\varepsilon,n} \\ & + (\boldsymbol{\mu}_{N,\mathcal{H}}, \mathbf{u}_h^n - \partial_{\tau}\boldsymbol{\phi}_{\mathcal{H}}^{n+1} \circ (\boldsymbol{\phi}_{\mathcal{H}}^n)^{-1})_{\varepsilon,n} = 0 \end{aligned} \right. \quad (4.3) $$

for all $(\mathbf{v}_h, q_h, \boldsymbol{\mu}_{N,\mathcal{H}}, \delta\mathbf{r}_{\mathcal{H}}, \boldsymbol{\theta}_{\mathcal{H}}) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_{N,\mathcal{H}}^n \times \mathbf{Y}_{\mathcal{H}} \times \mathbf{Y}_{\mathcal{H}}$.

In Step (4.3) of Algorithm 1, the discrete fluid tri-linear form is given by

$$ \begin{align*} a_{\Omega,h}^f(\mathbf{z}_h; (\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) &\stackrel{\text{def}}{=} c_h^n(\mathbf{z}_h; \mathbf{u}_h, \mathbf{v}_h) + a_{\Omega}((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) + s_{\Omega,h}(\mathbf{z}_h; (\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)), \\ \text{with} \quad c_h^n(\mathbf{z}_h; \mathbf{u}_h, \mathbf{v}_h) &\stackrel{\text{def}}{=} c_{\Omega}(\mathbf{z}_h; \mathbf{u}_h, \mathbf{v}_h) + \frac{\rho^f}{2}((\nabla \cdot \mathbf{z}_h) \mathbf{u}_h, \mathbf{v}_h)_{\Omega}, \\ s_{\Omega,h}(\mathbf{z}_h; (\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) &\stackrel{\text{def}}{=} \sum_{K \in \mathcal{T}_h(\Omega)} \int_K \delta_h (\rho^f(\mathbf{z}_h \cdot \nabla)\mathbf{u}_h + \nabla p_h) \cdot (\rho^f(\mathbf{z}_h \cdot \nabla)\mathbf{v}_h + \nabla q_h), \\ \delta_h &\stackrel{\text{def}}{=} \gamma_S \left( \rho^f \sqrt{\frac{4}{\tau^2} + \frac{16\mu^2}{h^4(\rho^f)^2} + \frac{4|\mathbf{z}_h|^2}{h^2}} \right)^{-1}, \end{align*} $$

where $s_{\Omega,h}$ corresponds to the SUPG/PSPG stabilization (see, e.g., [66]) and $\gamma_S > 0$ is a user-defined parameter. For the discrete beam bi-linear form, we consider the locking-free formulation introduced in [67], viz.,

$$ a_{\mathcal{H}}^b((\mathbf{r}_{\mathcal{H}}, \boldsymbol{\theta}), (\mathbf{r}_{\mathcal{H}}, \delta\mathbf{r}_{\mathcal{H}})) \stackrel{\text{def}}{=} (\mathbf{E}\partial_s\boldsymbol{\theta}_{\mathcal{H}}, \partial_s\delta\boldsymbol{\theta}_{\mathcal{H}})_{(0,L)} + (\mathbf{G}\Pi_{\mathcal{H}}(\partial_s\mathbf{r}_{\mathcal{H}}^n - \boldsymbol{\theta}_{\mathcal{H}}^n \wedge \mathbf{e}_z), \partial_s\delta\mathbf{r}_{\mathcal{H}} - \delta\boldsymbol{\theta}_{\mathcal{H}} \wedge \mathbf{e}_z)_{(0,L)}, $$

where $\Pi_{\mathcal{H}}$ denotes the $L^2$ projection onto $\mathbf{X}_{\mathcal{H}}^0(0, L)$.

*Energy estimate (small rotational velocity).* In the spirit of Section 3.2, we derive an energy estimate for Algorithm 1 under a small rotational velocity framework of (3.6). In this context, Algorithm 1 reduces
to the following discrete coupled problem: find $(\mathbf{u}_h^n, p_h^n, \lambda_{N,\mathcal{H}}^n, \mathbf{r}_\mathcal{H}^n, \boldsymbol{\theta}_\mathcal{H}^n) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_{N,\mathcal{H}}^n \times \mathbf{Y}_\mathcal{H} \times \mathbf{Y}_\mathcal{H}$ with $\dot{\mathbf{r}}_\mathcal{H}^{n-\frac{1}{2}} = \partial_\tau \mathbf{r}_\mathcal{H}^n$ and $\dot{\boldsymbol{\theta}}_\mathcal{H}^{n-\frac{1}{2}} = \partial_\tau \boldsymbol{\theta}_\mathcal{H}^n$ such that

$$
\left\{
\begin{array}{l}
\phi_{\mathcal{H}}^{n} = (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z + \overline{\mathbf{r}_{\mathcal{H}}^{n-1}} + \overline{\Lambda_{\mathcal{H}}^{n-1}}(\mathbf{I}_3 - (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z) \quad \text{in } \hat{\omega}_{\varepsilon}, \quad \Sigma_{\varepsilon}^{n} = \phi_{\mathcal{H}}^{n}(\hat{\Sigma}_{\varepsilon}), \\
\rho^{f}(\partial_{\tau} \mathbf{u}_{h}^{n}, \mathbf{v}_{h})_{\Omega} + a_{\Omega,h}^{f}(\mathbf{u}_{h}^{n-1}; (\mathbf{u}_{h}^{n}, p_{h}^{n}, (\mathbf{v}_{h}, q_{h}))) \\
\qquad + \rho^{b}(\mathbf{A} \partial_{\tau} \dot{\mathbf{r}}_{\mathcal{H}}^{n}, \delta \mathbf{r}_{\mathcal{H}})_{(0,L)} + \rho^{b}(\mathbf{I} \partial_{\tau} \dot{\boldsymbol{\theta}}_{\mathcal{H}}^{n}, \delta \boldsymbol{\theta}_{\mathcal{H}})_{(0,L)} + a_{\mathcal{H}}^{b}((\mathbf{r}_{\mathcal{H}}^{n-\frac{1}{2}}, \boldsymbol{\theta}_{\mathcal{H}}^{n-\frac{1}{2}}), (\delta \mathbf{r}_{\mathcal{H}}, \delta \boldsymbol{\theta}_{\mathcal{H}})) \\
\qquad - (\boldsymbol{\lambda}_{N,\mathcal{H}}^{n}, \mathbf{v}_{h} - \overline{\mathbf{r}_{\mathcal{H}}}) \circ (\phi_{\mathcal{H}}^{n})^{-1} - \overline{\delta \boldsymbol{\theta}_{\mathcal{H}}} \circ (\phi_{\mathcal{H}}^{n})^{-1} \wedge (\mathbf{I}_3 - \overline{\mathbf{r}_{\mathcal{H}}^{n-1}} \circ (\phi_{\mathcal{H}}^{n})^{-1} - ((\phi_{\mathcal{H}}^{n})^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\varepsilon,n} \\
\qquad + (\boldsymbol{\mu}_{N,\mathcal{H}}, \mathbf{u}_{h}^{n} - \overline{\dot{\mathbf{r}}_{\mathcal{H}}^{n-\frac{1}{2}}}) \circ (\phi_{\mathcal{H}}^{n})^{-1} - \overline{\dot{\boldsymbol{\theta}}_{\mathcal{H}}^{n-\frac{1}{2}}} \circ (\phi_{\mathcal{H}}^{n})^{-1} \wedge (\mathbf{I}_3 - \overline{\mathbf{r}_{\mathcal{H}}^{n-1}} \circ (\phi_{\mathcal{H}}^{n})^{-1} - ((\phi_{\mathcal{H}}^{n})^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\varepsilon,n} = 0.
\end{array}
\right.
\quad (4.4)
$$

for all $(\mathbf{v}_h, q_h, \boldsymbol{\mu}_{N,\mathcal{H}}, \delta \mathbf{r}_\mathcal{H}, \boldsymbol{\theta}_\mathcal{H}) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_{N,\mathcal{H}}^n \times \mathbf{Y}_\mathcal{H} \times \mathbf{Y}_\mathcal{H}$.

The energy stability of the numerical scheme (4.4) is stated in the following result.

**Theorem 4.1.** Let $\{(\mathbf{u}_h^n, p_h^n, \lambda_{N,\mathcal{H}}^n, \mathbf{r}_\mathcal{H}^n, \boldsymbol{\theta}_\mathcal{H}^n)\}_{n \ge 1} \subset \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_{N,\mathcal{H}}^n \times \mathbf{Y}_\mathcal{H} \times \mathbf{Y}_\mathcal{H}$ be given by problem (4.4). There holds

$$
\mathfrak{E}^n \le \mathfrak{E}^0, \quad \forall n \ge 0 \qquad (4.5)
$$

where the discrete mechanical energy $\mathfrak{E}^n$ of the system is defined by

$$
\mathfrak{E}^n \stackrel{\text{def}}{=} \frac{\rho^f}{2} \|\mathbf{u}_h^n\|_{0,\Omega}^2 + \frac{\rho^b}{2} \|\dot{\mathbf{r}}_\mathcal{H}^n\|_{\mathbf{A},(0,L)}^2 + \frac{\rho^b}{2} \|\dot{\boldsymbol{\theta}}_\mathcal{H}^n\|_{\mathbf{I},(0,L)}^2 + \frac{1}{2} \|(\mathbf{r}_\mathcal{H}^n, \boldsymbol{\theta}_\mathcal{H}^n)\|_{\mathbf{b},\mathcal{H}}^2,
$$

with

$$
\| (\cdot, \cdot) \|_{\mathbf{b},\mathcal{H}} = \sqrt{a_{\mathcal{H}}^{b}((\cdot, \cdot), (\cdot, \cdot))}.
$$

Proof. By testing (4.4) with $(\mathbf{v}_h, q_h, \boldsymbol{\mu}_{N,\mathcal{H}}, \delta \mathbf{r}_\mathcal{H}, \delta \boldsymbol{\theta}_\mathcal{H}) = (\mathbf{u}_h^n, p_h, \lambda_{N,\mathcal{H}}^n, \overline{\dot{\mathbf{r}}_\mathcal{H}^{n-\frac{1}{2}}}, \overline{\dot{\boldsymbol{\theta}}_\mathcal{H}^{n-\frac{1}{2}}})$, we get

$$
\begin{equation}
\begin{split}
& \rho^f(\partial_\tau \mathbf{u}_h^n, \mathbf{u}_h^n)_\Omega + \rho^b(\mathbf{A} \partial_\tau \dot{\mathbf{r}}_\mathcal{H}^n, \dot{\mathbf{r}}_\mathcal{H}^{n-\frac{1}{2}})_{(0,L)} + \rho^b(\mathbf{I} \partial_\tau \dot{\boldsymbol{\theta}}_\mathcal{H}^n, \dot{\boldsymbol{\theta}}_\mathcal{H}^{n-\frac{1}{2}})_{(0,L)} \\
& \qquad + a_{\Omega,h}^f(\mathbf{u}_h^{n-1}; (\mathbf{u}_h^n, p_h^n), (\mathbf{u}_h^n, p_h^n)) + a_{\mathcal{H}}^b((\mathbf{r}_\mathcal{H}^{n-\frac{1}{2}}, \boldsymbol{\theta}_\mathcal{H}^{n-\frac{1}{2}}), (\dot{\mathbf{r}}_\mathcal{H}^{n-\frac{1}{2}}, \dot{\boldsymbol{\theta}}_\mathcal{H}^{n-\frac{1}{2}})) = 0,
\end{split}
\tag{4.6}
\end{equation}
$$

Using integration by parts and the boundary conditions of $\mathbf{u}_h^{n-1}$, $\mathbf{u}_h^n$, we have

$$
a_{\Omega,h}^{f}(\mathbf{u}_{h}^{n-1}; (\mathbf{u}_{h}^{n}, p_{h}), (\mathbf{u}_{h}^{n}, p_{h})) \geq 2\mu^{f}\|\varepsilon(\mathbf{u}_{h}^{n})\|_{0,\Omega}^{2}.
$$

Moreover for any quantities $\mathbf{v}^n, \mathbf{v}^{n-1}$ and inner product $(\cdot, \cdot)$, it holds

$$
(\mathbf{v}^n - \mathbf{v}^{n-1}, \mathbf{v}^n) = \frac{1}{2} (\|\mathbf{v}^n\|^2 - \|\mathbf{v}^{n-1}\|^2 + \|\mathbf{v}^n - \mathbf{v}^{n-1}\|^2) \ge \frac{1}{2} \partial_\tau \|\mathbf{v}^n\|^2, \quad (\mathbf{v}^n - \mathbf{v}^{n-1}, \mathbf{v}^{n-\frac{1}{2}}) = \frac{1}{2} \partial_\tau \|\mathbf{v}^n\|^2,
$$

where $\| \cdot \|_{\text{def}} = \sqrt{(\cdot, \cdot)}$. We deduce that

$$
\partial_{\tau} \left( \frac{\rho^{f}}{2} \| \mathbf{u}_{h}^{n} \|_{\Omega}^{2} + \frac{\rho^{b}}{2} \| \dot{\mathbf{r}}_{\mathcal{H}}^{n} \|_{\mathbf{A}, (0, L)}^{2} + \frac{\rho^{b}}{2} \| \dot{\boldsymbol{\theta}}_{\mathcal{H}}^{n} \|_{\mathbf{I}, (0, L)}^{2} + \frac{1}{2} \| \mathbf{r}_{\mathcal{H}}^{n}, \boldsymbol{\theta}_{\mathcal{H}}^{n} \|_{\mathbf{b}, \mathcal{H}}^{2} \right) \le 0.
$$

The estimate (4.5) then follows by summing over $\{0...n\}$, which completes the proof. $\square$

**Inf-sup stability.** A crucial point when dealing with Lagrange multipliers is to ensure that the inf-sup condition associated to the bi-linear form $(\lambda_{N,\mathcal{H}}^n, \mathbf{v}_h)_{\varepsilon,n}$ is satisfied on $\mathbf{F}_{N,\mathcal{H}}^n \times \mathbf{V}_h$. In [51], the authors establish the inf-sup condition for the case $N = 0$ under some conditions on the triangulations of $\Omega$, $(0, L)$ and $\Sigma_\varepsilon^n$. To deal with general meshes, they propose a discretization of $(0, L)$ depending on the intersection of the fluid elements and the beam centerline. With the addition of a stabilization term, they established the inf-sup condition on $\mathbf{F}_{0,\mathcal{H}}^n \times \mathbf{V}_h$ where the coefficients $\mathbf{a}_{k,\mathcal{H}}$, $\mathbf{b}_{k,\mathcal{H}}$ introduced in (4.1) and defining the space $\mathbf{F}_{0,\mathcal{H}}^n$ belong to $\mathbf{X}_h^0(0, L)$ instead of $\mathbf{X}_h^1(0, L)$. If $N > 0$, the authors of [52] considered a toy two-dimensional problem and proved a stability result on $\mathbf{F}_{N,\mathcal{H}}^n \times (\mathbf{X}_h^k(\Omega) \cap \mathbf{V})$,
provided $N \le k$. In the same framework, they also introduced a general stabilization term, which in 3D would take the form (see also [68]):

$$ \mathfrak{s}_{\Sigma_{\varepsilon}^{n}, h}(\boldsymbol{\mu}_{N, \mathcal{H}}, \boldsymbol{\lambda}_{N, \mathcal{H}}^{n}) = \gamma \varepsilon^{-1} \left(\frac{h}{\varepsilon}\right) (\boldsymbol{\mu}_{N, \mathcal{H}}, \boldsymbol{\lambda}_{N, \mathcal{H}}^{n})_{\varepsilon, n}, $$

and proved unconditional stability for any pair $(\lambda_{N,\mathcal{H}}^n, v_h) \in F_{N,\mathcal{H}}^n \times (X_h^k(\Omega) \cap V)$ (i.e., with any restriction on $N$). However, no stability proof has been provided in the specific case addressed here. Nevertheless, as we encountered no issues during the numerical simulations, we chose not to include any stabilization term in our numerical study.

*Evaluation of the coupling terms.* In this section, we provide details on the computation of the coupling terms involved in step (4.3) of Algorithm 1. For the sake of conciseness, we only give the expressions of the algebraic counterparts of the bi-linear form $(\lambda_{N,\mathcal{H}}^n, v_h)_{\varepsilon,n}$. The evaluation of the remaining coupling terms follow similarly. Let $\{\psi_i^f\}_{1\le i\le N_h}$ and $\{\psi_j^b\}_{1\le j\le N_h}$ be the basis functions of $V_h$ and $Y_{\mathcal{H}}$, respectively, and let $\{x_i\}_{1\le i\le N_h}$ and $\{s_j\}_{1\le j\le N_h}$ be the corresponding nodes. We hence have

$$ (\lambda_{N,\mathcal{H}}^n, v_h)_{\varepsilon,n} = \sum_{i=1}^{N_h} v_h(x_i) (\lambda_{N,\mathcal{H}}^n, \psi_i^f)_{\varepsilon,n}. $$

Using the definition of $F_{N,\mathcal{H}}^n$ given in (4.1) with tensorized notations, we have

$$ (\lambda_{N,\mathcal{H}}^n, v_h)_{\varepsilon,n} = \sum_{i=1}^{N_h} v_h(x_i) \left( \left( \overline{a_{0,\mathcal{H}}} + \sum_{k=0}^{N} (a_{k,\mathcal{H}} \otimes c_k + b_{k,\mathcal{H}} \otimes s_k) \right) \circ (\phi_{\mathcal{H}}^n)^{-1}, \psi_i^f \right)_{\varepsilon,n}, $$

with $c_k(\nu) \stackrel{\text{def}}{=} \cos(k\nu)$ and $s_k(\nu) \stackrel{\text{def}}{=} \sin(k\nu)$ for $0 \le k \le N$. Now, by expression $a_{k,\mathcal{H}}, b_{k,\mathcal{H}}$ in terms of their basis functions $\psi_j^b$, we get

$$ (\lambda_{N,\mathcal{H}}^n, v_h)_{\varepsilon,n} = \sum_{i=1}^{N_h} v_h(x_i) \sum_{k=0}^{N} \sum_{j=1}^{N_h} \left[ \begin{aligned} & a_{0,\mathcal{H}}(s_j) (\overline{\psi_j^b} \circ (\phi_{\mathcal{H}}^n)^{-1}, \psi_i^f)_{\varepsilon,n} \\ & + a_{k,\mathcal{H}}(s_j) ((\psi_j^b \otimes c_k) \circ (\phi_{\mathcal{H}}^n)^{-1}, \psi_i^f)_{\varepsilon,n} + b_{k,\mathcal{H}}(s_j) ((\psi_j^b \otimes s_k) \circ (\phi_{\mathcal{H}}^n)^{-1}, \psi_i^f)_{\varepsilon,n} \end{aligned} \right]. $$

Let $(\mathbf{B}^N)^T \in \mathbb{R}^{N_{\mathcal{H}}(2N+1) \times N_h}$ denote the algebraic counterpart of the bi-linear term $(\lambda_{N,\mathcal{H}}^n, v)_{\varepsilon,n}$, we get

$$ (\mathbf{B}^N)_{il} = \begin{cases} (\overline{\psi_j^b} \circ (\phi_{\mathcal{H}}^n)^{-1}, \psi_i^f)_{\varepsilon,n}, & l = j \cdot (2N + 1), \\ ((\psi_j^b \otimes c_k) \circ (\phi_{\mathcal{H}}^n)^{-1}, \psi_i^f)_{\varepsilon,n}, & l = j \cdot (2N + 1) + 2k + 1, \quad \forall k \in \{0, \dots, N\}. \\ ((\psi_j^b \otimes s_k) \circ (\phi_{\mathcal{H}}^n)^{-1}, \psi_i^f)_{\varepsilon,n}, & l = j \cdot (2N + 1) + 2k + 2. \end{cases} \quad (4.7) $$

In order to evaluate the coefficients (4.7), we need to evaluate the tensorized basis functions $(\psi_j^b \otimes c_k) \circ (\phi_{\mathcal{H}}^n)^{-1}, (\psi_j^b \otimes s_k) \circ (\phi_{\mathcal{H}}^n)^{-1}$ on $\Sigma_{\varepsilon}^n$. To do so, we first mesh $\hat{\Sigma}_{\varepsilon}$ and evaluate them at the mesh points. We then transport the values to the current configuration via the deformation map $\phi_{\mathcal{H}}^n$. Integration is achieved through $\mathbb{P}_1$-interpolation in the deformed mesh by localizing the quadrature points in the fluid mesh. Note that we choose a mesh of $\hat{\Sigma}_{\varepsilon}$ fine enough so that the error due to the $\mathbb{P}_1$-interpolation is negligible compared to the numerical approximation error of the PDE.

*Matrix formulation.* We consider a partitioned solution procedure for the numerical resolution of problem (4.3) based on [69]. To this purpose, we introduce $U^n, P^n, L^n, R^n$ and $T^n$ as the arrays of degrees of freedom associated with $u_h^n, p_h^n, \lambda_{N,\mathcal{H}}^n, r_h^n$ and $\theta_{\mathcal{H}}^n$ respectively. Taking $(\delta r_{\mathcal{H}}, \delta \theta_{\mathcal{H}}) = (0,0)$ in (4.3), we recover the Navier-Stokes equations with Dirichlet boundary condition on $\Sigma_{\varepsilon}^n$. This problem can be expressed in matrix form as follows:

$$ \begin{bmatrix} A^f & C^T & -(B^N)^T \\ -C & S & 0 \\ B^N & 0 & 0 \end{bmatrix} \begin{bmatrix} U^n \\ P^n \\ L^n \end{bmatrix} = \begin{bmatrix} b^{n-1} \\ 0 \\ B^N B^b (R^n, T^n) \end{bmatrix}, \quad (4.8) $$
with

$$A^f \stackrel{\text{def}}{=} \frac{\rho^f}{\tau} M^f + K^f, \quad b^{n-1} \stackrel{\text{def}}{=} \frac{\rho^f}{\tau} M^f U^{n-1}, \quad U^b(R, T) \stackrel{\text{def}}{=} \frac{1}{\tau} (\phi(R, T) - \phi(R^{n-1}, T^{n-1})).$$

Here, the matrices $M^f$, $[\begin{matrix} K^f & C^T \\ -C & S \end{matrix}]$ and $B^N$ stand for the algebraic counterpart of

$$(u_h^n, v_h)_\Omega, \quad a_{\Omega,h}^f(u_h^{n-1}; (u_h^n, p_h^n), (v_h, q_h)), \quad (\lambda_{N,\mathcal{H}}^n, v)_{\varepsilon,n}.$$

With a slight abuse of notation, the operator $\phi$ allows us to compute the beam deformation given the displacement $R^n$ of the centerline of the beam and its rotation vector $T^n$. To the matrix system (4.8), we can associate the fluid solver operator $\mathcal{F}^n : \mathbb{R}^{N^b} \to \mathbb{R}^{N^t}$ given by

$$\mathcal{F}^n \left( \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} \right) = \mathbf{L}^n.$$

In a similar manner, testing (4.3) with $(v_h, q_h, \mu_{N,\mathcal{H}}) = (0, 0, 0)$ yields the linear Timoshenko sup-problem with fluid source terms. The equivalent matrix system writes as follows:

$$A^b \begin{bmatrix} \mathbf{R}^{n-\frac{1}{2}} \\ \mathbf{T}^{n-\frac{1}{2}} \end{bmatrix} = \mathbf{J}^{n-1} + \begin{bmatrix} \mathbf{F}^n \\ \mathbf{M}^n \end{bmatrix} \qquad (4.9)$$

with $\mathbf{R}^n = 2\mathbf{R}^{n-\frac{1}{2}} - \mathbf{R}^{n-1}$, $\mathbf{T}^n = 2\mathbf{T}^{n-\frac{1}{2}} - \mathbf{T}^{n-1}$ and where

$$\begin{align*}
\mathbf{A}^b \stackrel{\text{def}}{=} \frac{4\rho^b}{\tau^2} \mathbf{M}^b + \mathbf{K}^b, \quad \mathbf{J}^{n-1} \stackrel{\text{def}}{=} \frac{\rho^b}{\tau^2} \mathbf{M}^b \left( 4 \begin{bmatrix} \mathbf{R}^{n-1} \\ \mathbf{T}^{n-1} \end{bmatrix} + 2\tau \begin{bmatrix} \dot{\mathbf{R}}^{n-1} \\ \dot{\mathbf{T}}^{n-1} \end{bmatrix} \right), \\
\mathbf{F}^n = (\mathbf{E}_f^N)^\mathrm{T} \mathbf{L}^n, \quad \mathbf{M}^n = (\mathbf{E}_m^N)^\mathrm{T} \mathbf{L}^n, \\
\dot{\mathbf{R}}^{n-\frac{1}{2}} = \frac{1}{\tau}(\mathbf{R}^n - \mathbf{R}^{n-1}), \quad \dot{\mathbf{T}}^{n-\frac{1}{2}} = \frac{1}{\tau}(\mathbf{T}^n - \mathbf{T}^{n-1}).
\end{align*}$$

Here, the matrices $\mathbf{M}^b, \mathbf{K}^b, \mathbf{E}_f^N, \mathbf{E}_m^N$ stand for the algebraic counterpart of the bi-linear forms

$$\begin{align*}
(\mathbf{A}\mathbf{r}_H^n, \delta\mathbf{r}_H)_{(0,L)} + (\mathbf{I}\boldsymbol{\theta}_H^n, \delta\boldsymbol{\theta}_H)_{(0,L)}, \quad & a_H^b((\mathbf{r}_H^n, \boldsymbol{\theta}_H^n), (\delta\mathbf{r}_H, \delta\boldsymbol{\theta}_H)), \quad (\lambda_{N,\mathcal{H}}^n, \overline{\delta\mathbf{r}} \circ (\boldsymbol{\phi}_H^n)^{-1})_{\varepsilon,n}, \\
(\lambda_{N,\mathcal{H}}^n, \overline{\delta\boldsymbol{\theta}} \circ (\boldsymbol{\phi}_H^n)^{-1} \wedge (\mathbf{I}_3 - \overline{\mathbf{r}_H^{n-1}} \circ (\boldsymbol{\phi}_H^n)^{-1} - ((\boldsymbol{\phi}_H^n)^{-1} \cdot \mathbf{e}_z)\mathbf{e}_z))_{\varepsilon,n}.
\end{align*}$$

The corresponding solid solution operator $\mathcal{S}^n : \mathbb{R}^{N^t} \to \mathbb{R}^{N^b}$ is defined as

$$\mathcal{S}^n \left( \begin{bmatrix} \mathbf{U}^n \\ \mathbf{P}^n \\ \mathbf{L}^n \end{bmatrix} \right) = \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix}.$$

By composition, solving the coupled problem (4.3) is equivalent to compute the roots of the following system

$$\mathcal{R}^n \left( \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} \right) \stackrel{\text{def}}{=} \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} - \mathcal{S}^n \circ \mathcal{F}^n \left( \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} \right) = \mathbf{0}, \qquad (4.10)$$

which can be iteratively approximated using Newton's method below, where $\mathcal{J}^n$ denotes the Jacobian of the residual operator $\mathcal{R}^n$. As in [69], the tangent system (4.11) can be solved in a matrix-free fashion via GMRES iterations, which only invoke the action of the operator $\mathcal{J}^n$ along given directions. Note that the resulting solution procedure is fully partitioned, in the sense that only the fluid and solid operators, and their tangent counterparts, are involved in the numerical resolution of (4.3).
Newton algorithm applied to (4.10)

$$
\begin{array}{ll}
\text{Choose } \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} \in \mathbb{R}^{N^b} & \\
\textbf{while } \left\| \mathcal{R}^n \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} \right\| \ge \text{tol do} & \\
\quad \text{Evaluate the fluid operator } \begin{bmatrix} \mathbf{U}^n \\ \mathbf{P}^n \end{bmatrix} = \mathcal{F}^n \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} & \\
\quad \text{Evaluate the solid operator } \begin{bmatrix} \hat{\mathbf{R}}^n \\ \hat{\mathbf{T}}^n \end{bmatrix} = \mathcal{S}^n \begin{bmatrix} \mathbf{U}^n \\ \mathbf{P}^n \end{bmatrix} & \\
\quad \text{Evaluate the residual } \mathcal{R}^n \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} = \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} - \begin{bmatrix} \hat{\mathbf{R}}^n \\ \hat{\mathbf{T}}^n \end{bmatrix} & \\
\quad \text{Solve tangent problem} & \\
\quad \mathcal{J}^n \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} \begin{bmatrix} \delta \mathbf{R} \\ \delta \mathbf{T} \end{bmatrix} = -\mathcal{R}^n \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} & (4.11) \\
\text{Update } \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} = \begin{bmatrix} \mathbf{R}^n \\ \mathbf{T}^n \end{bmatrix} + \begin{bmatrix} \delta \mathbf{R} \\ \delta \mathbf{T} \end{bmatrix} & \\
\textbf{end while} &
\end{array}
$$

## 4.2. Comparison with other discrete formulations

For the sake of completeness, in this section we briefly present the discrete formulations of the coupled problems (2.9)-(2.11) and (3.2) discussed above. These numerical methods will be used, for comparison purposes, in the numerical experiments of Section 4.

Full coupled problem with ALE formalism in the fluid. When the interface $\Sigma_{\epsilon}(t)$ undergoes moderate displacements, it is convenient (from the numerical point of view) to parametrize the fluid domain $\Omega_{\epsilon}^{f}(t)$ by a one-to-one mapping $\mathcal{A}: \hat{\Omega}_{\epsilon}^{f} \times \mathbb{R}^{+} \rightarrow \mathbb{R}^{3}$, defined as a suitable lifting of the solid motion:

$$ \mathcal{A}_t = \mathcal{L}(\phi_t|_{\hat{\Omega}_\epsilon}), \quad \Omega_\epsilon^f(t) = \mathcal{A}_t(\hat{\Omega}_\epsilon^f), \qquad (4.12) $$

with $\mathcal{A}_t(\cdot) \stackrel{\text{def}}{=} \mathcal{A}(\cdot, t)$ and where $\mathcal{L}$ denotes a given lifting operator from $\hat{\Sigma}_\epsilon$ to $\hat{\Omega}_\epsilon^f$ which vanishes on $\partial\Omega$. Note the extension provided by $\mathcal{L}$ is arbitrary inside $\hat{\Omega}_\epsilon^f$ and the relation (4.12) is simply a reformulation of the geometric compatibility condition, i. e. the first equation of coupling equations (2.11). The fluid equations in ALE formalism are simply obtained from $(2.9)_1$, by replacing the Eulerian time derivative by the standard relation $\partial_t u = \partial_t |_{\mathcal{A}} u - (\hat{w}^f \circ \mathcal{A}_t^{-1}) \cdot \nabla u$, where $\hat{w}^f \stackrel{\text{def}}{=} \partial_t \mathcal{A}$ denotes the fluid domain velocity and $\partial_t |_{\mathcal{A}}$ stands for the so-called ALE time-derivative. This change of variables facilitates the time discretization of quantities defined in moving meshes.

We consider a triangulation of the reference fluid domain $\hat{\Omega}_\epsilon^f$ which is fitted to the triangulation of the interface $\hat{\Sigma}_\epsilon$. We set $X_h^1(\hat{\Omega}_\epsilon)$ as the trace space of $X_h^1(\hat{\Omega}_\epsilon^f)$ and $\Pi_h^b$ the corresponding Lagrange interpolation operator onto $X_h^1(\hat{\Omega}_\epsilon)$. The geometric coupling condition (4.12) is treated explicitly. For a given displacement $r_H^{n-1} \in Y_H$ and rotation vector $\theta_H^{n-1} \in Y_H$ at time step $n-1$, we define the ALE map $\mathcal{A}_h^n$ as

$$ \phi_H^n = (I_3 \cdot e_z) e_z + \overline{r_H^{n-1}} + \overline{\Lambda_H^{n-1}} (I_3 - (I_3 \cdot e_z) e_z) \text{ in } \hat{\omega}_\epsilon, \quad \mathcal{A}_h^n = \mathcal{L}_h (\Pi_h^b \phi_H^n|_{\hat{\Sigma}_\epsilon}), $$

with $\Lambda_H^{n-1}$ given from $\theta_H^{n-1}$, using (2.3), and where $\mathcal{L}_h$ denotes a given discrete lifting operator (e.g., an harmonic extension operator). The one-to-one mapping $\mathcal{A}_h^n$ allows us to consider as fluid unknowns

$$ \hat{\mathbf{u}}_h^n \stackrel{\text{def}}{=} \mathbf{u}_h^n \circ \mathcal{A}_h^n, \quad \hat{p}_h^n \stackrel{\text{def}}{=} p_h^n \circ \mathcal{A}_h^n $$

instead of $\mathbf{u}_h^n = \hat{\mathbf{u}}_h^n \circ (\phi_H^n)^{-1}$ and $p_h^n = \hat{p}_h^n \circ (\phi_H^n)^{-1}$. We then introduce the following finite element space $V_h^f \stackrel{\text{def}}{=} X_h^1(\hat{\Omega}_\epsilon^f) \cap H_{\partial\Omega}^1(\hat{\Omega}_\epsilon^f)$ and $Q_h^f \stackrel{\text{def}}{=} Q_h^1(\hat{\Omega}_\epsilon^f) \cap L^2(\hat{\Omega}_\epsilon^f)$ for the approximation of the fluid unknowns. We write the fluid and structure equations in weak form and apply a finite element discretization with an implicit treatment of the coupling conditions. The resulting discrete formulation is given in Algorithm 2.
**Algorithm 2** Discrete formulation of the coupled problem (2.9)-(2.11) (ALE formalism in the fluid)

For $n \ge 1$,

**Step 1:** Update of the beam deformation and of the fluid domain

$$ \left\{ \begin{array}{l} \phi_{\mathcal{H}}^{n} = (\mathbf{I}_3 \cdot \mathbf{e}_z) \mathbf{e}_z + \overline{\mathbf{r}_{\mathcal{H}}^{n-1}} + \overline{\mathbf{\Lambda}_{\mathcal{H}}^{n-1}} (\mathbf{I}_3 - (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z) \quad \text{in } \hat{\omega}_{\varepsilon}, \\ \mathbf{A}_{h}^{n} = \mathcal{L}_{h}(\Pi_{h}^{\text{b}}\phi_{\mathcal{H}}^{n}|_{\hat{\Omega}_{\varepsilon}}), \quad \Omega_{\varepsilon}^{\text{f},n} = \mathbf{A}_{h}^{n}(\hat{\Omega}_{\varepsilon}^{\text{f}}), \quad \hat{w}_{h}^{\text{f},n} = \partial_{\tau} \mathbf{A}_{h}^{n}, \end{array} \right. $$

with $\mathbf{\Lambda}_{\mathcal{H}}^{n-1}$ given from $\theta_{\mathcal{H}}^{n-1}$ using (2.3).

**Step 2:** Find $(\hat{\mathbf{u}}_h^n, \hat{p}_h^n, \mathbf{r}_\mathcal{H}^n, \theta_\mathcal{H}^n) \in \mathbf{V}_h^f \times \mathbf{Q}_h^f \times \mathbf{Y}_\mathcal{H} \times \mathbf{Y}_\mathcal{H}$ with $\dot{\mathbf{r}}_\mathcal{H}^{n-\frac{1}{2}} = \partial_\tau \mathbf{r}_\mathcal{H}^n, \dot{\theta}_\mathcal{H}^{n-\frac{1}{2}} = \partial_\tau \theta_\mathcal{H}^n$ and such that

$$ \left\{ \begin{array}{l} \hat{\mathbf{u}}_h^n = \Pi_h^b \partial_\tau \phi_\mathcal{H}^n \quad \text{on } \hat{\Sigma}_\varepsilon, \\ \frac{\rho^f}{\tau} \left[ (\mathbf{u}_h^n, \mathbf{v}_h)_{\Omega_\varepsilon^{f,n}} - (\mathbf{u}_h^{n-1}, \mathbf{v}_h)_{\Omega_\varepsilon^{f,n-1}} \right] - \rho^f \left( (\nabla \cdot \mathbf{w}_h^{f,n}) \mathbf{u}_h^n, \mathbf{v}_h \right)_{\Omega_\varepsilon^{f,n}} \\ \qquad + a_{\Omega_\varepsilon^{f,n},h} (\mathbf{u}_h^{n-1} - \mathbf{w}_h^{f,n}; (\mathbf{u}_h^n, p_h^n), (\mathbf{v}_h, q_h)) + \rho^b (\mathbf{A} \partial_\tau \dot{\mathbf{r}}_\mathcal{H}^n, \delta \mathbf{r}_\mathcal{H})_{(0,L)} \\ \qquad + \rho^b (\mathbf{I} \partial_\tau \dot{\theta}_\mathcal{H}^n, \delta \theta_\mathcal{H})_{(0,L)} + a_\mathcal{H}^b ((\mathbf{r}_\mathcal{H}^{n-\frac{1}{2}}, \theta_\mathcal{H}^{n-\frac{1}{2}}), (\delta \mathbf{r}_\mathcal{H}, \delta \theta_\mathcal{H})) = 0 \end{array} \right. \quad (4.13) $$

for all $(\hat{\mathbf{v}}_h, \hat{q}_h, \delta \mathbf{r}_\mathcal{H}, \delta \theta_\mathcal{H}) \in \mathbf{V}_h^f \times \mathbf{Q}_h^f \times \mathbf{Y}_\mathcal{H} \times \mathbf{Y}_\mathcal{H}$ with

$$ \hat{v}_h = \Pi_h^b (\bar{\delta \mathbf{r}}_\mathcal{H} + \bar{\delta \theta}_\mathcal{H} \wedge (\phi_\mathcal{H}^n - \overline{\mathbf{r}_\mathcal{H}^{n-1}} - (\mathbf{I}_3 \cdot \mathbf{e}_z)\mathbf{e}_z)) \quad \text{on } \hat{\Sigma}_\varepsilon. $$

One limitation of Algorithm 2 lies in the capability of the discrete lifting operator $\mathcal{L}_h$ to guarantee mesh quality when dealing with large interface deflections. This problem is even more pronounced in the context of multiple slender structures immersed in a fluid and that can get into contact. Algorithm 1 overcomes this issue by working with a fixed mesh, but it requires the localization of the interfacial mesh within the fluid mesh and the evaluation of the corresponding coupling terms.

*Coupled problem with coupling conditions on the centerline.* We present here the discrete formulation for the reduced coupled problem (3.2). As in the previous discrete formulations, the geometrical coupling is treated in an explicit fashion as follows:

$$ \phi_{\mathcal{H}}^{n} = \mathbf{r}_{\mathcal{H}}^{n-1}, \quad \mathcal{B}^{n} = \phi_{\mathcal{H}}^{n}(\{0\} \times \{0\} \times (0, L)). $$

Then the discrete problem is similar to the one given in Algorithm 1 but with the Lagrange multipliers belonging to $\mathbf{Y}_\mathcal{H}^n = \mathbf{X}_\mathcal{H}^1(\mathcal{B}^n)$ instead of $\mathbf{F}_\mathcal{N},_\mathcal{H}^n$ and the integrals over $\Sigma_\varepsilon^n$ replaced with integrals over $\mathcal{B}^n$. The details of this formulation can be found in Algorithm 3.

**Algorithm 3** Discrete formulation of the reduced coupled problem (3.2).

For $n \ge 1$,

**Step 1:** Update of the beam centerline location

$$ \phi_{\mathcal{H}}^{n} = \mathbf{I}_{3} + \mathbf{r}_{\mathcal{H}}^{n-1} \quad \text{in} \quad \{0\} \times \{0\} \times (0, L), \quad \mathcal{B}^{n} = \phi_{\mathcal{H}}^{n}(\{0\} \times \{0\} \times (0, L)). $$

**Step 2:** Find $(\mathbf{u}_h^n, p_h^n, \lambda_\mathcal{H}^n, \mathbf{r}_\mathcal{H}^n, \theta_\mathcal{H}^n) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{Y}_\mathcal{H}^n \times \mathbf{Y}_\mathcal{H} \times \mathbf{Y}_\mathcal{H}$, with $\dot{\mathbf{r}}_\mathcal{H}^{n-\frac{1}{2}} = \partial_\tau \mathbf{r}_\mathcal{H}^n$ and $\dot{\theta}_\mathcal{H}^{n-\frac{1}{2}} = \partial_\tau \theta_\mathcal{H}^n$, such that

$$ \begin{aligned} & \rho^f (\partial_\tau \mathbf{u}_h^n, \mathbf{v}_h)_{\Omega} + a_{\Omega,h}^f (\mathbf{u}_h^{n-1}; (\mathbf{u}_h^n, p_h^n), (\mathbf{v}_h, q_h)) \\ & + \rho^b (\mathbf{A} \partial_\tau \dot{\mathbf{r}}_\mathcal{H}^n, \delta \mathbf{r}_\mathcal{H})_{(0,L)} + \rho^b (\mathbf{I} \partial_\tau \dot{\theta}_\mathcal{H}^n, \delta \theta_\mathcal{H})_{(0,L)} + a_\mathcal{H}^b ((\mathbf{r}_\mathcal{H}^{n-\frac{1}{2}}, \theta_\mathcal{H}^{n-\frac{1}{2}}), (\delta \mathbf{r}_\mathcal{H}, \delta \theta_\mathcal{H})) \\ & + (\lambda_\mathcal{H}^n, \mathbf{v}_h - \delta \mathbf{r}_\mathcal{H} \circ (\phi_\mathcal{H}^n)^{-1})_{\mathcal{B}^n} + (\mu_\mathcal{H}, \mathbf{u}_h^n - \partial_\tau \mathbf{r}_\mathcal{H}^n \circ (\phi_\mathcal{H}^n)^{-1})_{\mathcal{B}^n} = \mathbf{0} \end{aligned} \quad (4.14) $$

for all $(\mathbf{v}_h, q_h, \mu_\mathcal{H}, \delta \mathbf{r}_\mathcal{H}, \delta \theta_\mathcal{H}) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{Y}_\mathcal{H}^n \times \mathbf{Y}_\mathcal{H} \times \mathbf{Y}_\mathcal{H}$.
Compared to Algorithm 1, Algorithm 3 involves the localization of the beam centerline mesh within the fluid mesh, instead of the interfacial mesh, making the assembly of the matrix B^N, introduced in (4.7), faster. Furthermore, the matrix sparsity is slightly better for the Algorithm 3 as the number of fluid elements intersected by the 1D beam centerline is inherently smaller than when dealing with the 2D interface. For N = 0, passing from one reduced model to the other basically consists in computing the average value of the fluid velocity on each cross section instead of taking it at their center, then the number of degrees of freedom for the Lagrange multipliers space remains unchanged, and consequently, the size of the system matrix is the same.

## 5. Numerical experiments

In this section, we illustrate the accuracy of Algorithm 1 by comparing its numerical solution with those provided by Algorithms 2 and 3. The results provided by Algorithm 2 are taken as the reference solution. We first consider the case of a single beam and investigate the influence of several parameters, such as the mesh size h, the beam thickness ε, the number of Fourier modes N and the Reynolds number Re. In our context, the Reynolds number is given by $Re \stackrel{\text{def}}{=} 2u_{\text{ref}}\rho_f\varepsilon/\mu^l$ where $u_{\text{ref}}$ is a reference velocity whose value will be specified for each test case. To illustrate the robustness of the proposed modeling approach, we also consider the case of large displacements using a non-linear beam model, with multiple beams that can be in contact. All the numerical simulations have been performed with the FELiScE finite element library [70].

### 5.1. Single beam immersed in a Stokesian flow

In the first numerical example, we consider the mechanical interaction between an incompressible viscous fluid and a single beam in a rectangular domain of dimensions $l_x \times l_y \times l_z$ centered at (0.5, 0.5, 0). We assume that the fluid is described by Stokes equations, so that the inertia terms in (2.5) are neglected. Homogeneous Neumann boundary conditions are enforced at the outlet $\Gamma_{\text{out}}$, homogeneous Dirichlet condition on the lower face ($z = 0$) of the domain, and the perfect slip conditions on the other three lateral faces, as illustrated in Figure 3.

<figure><img src="image_5.png" /><figcaption>Figure 3: Geometric setting for the single beam immersed in a rectangular fluid domain.</figcaption></figure>

#### 5.1.1. Convergence study with respect ε and h

In this section, we consider a static problem where the beam is taken as a straight tube of length L=0.06 and the fluid domain parameters are given by $l_x = l_y = l_z = 1$. The following boundary condition is enforced at the inlet boundary $\Gamma_{\text{in}}$:

$$ \mathbf{u}_{\text{in}} = u_{\text{ref}}(1 - (z - 1)(z + 1))e_x, $$

with $u_{\text{ref}}=20$.
<figure><img src="image_6.png" /><figcaption>Figure 4: Cut plane of the fluid mesh for different resolutions of the beam interface and $\epsilon=0.06$.</figcaption></figure>

*Convergence in $\epsilon$.* We first investigate the modeling error between the solution of the reduced and full models, computed with Algorithms 1 and 2 respectively, for a number of Fourier modes $N=0$, and a decreasing beam thickness $\epsilon \in \{0.1, 0.08, 0.06, 0.04, 0.02, 0.01\}$. To solve both problems, we use a fitted mesh with respect to the beam interface where the mesh size $h$ is sufficiently small such that the numerical error of the discrete solutions can be neglected with respect to the model error (see Figures 5a and 5b for the case $\epsilon=0.06$). The $L^2$ and $H^1$ norms of the difference in the fluid solutions are reported in Figure 6. We observe that the $L^2$ error behaves as $\epsilon$ and the $H^1$ error as $\epsilon^{1/2}$.

*Convergence in $h$.* Similarly to the previous paragraph, we consider a static problem where the geometric parameters for the beam are given by $L=0.06$ and $\epsilon=0.1$. We analyze the error between the solution obtained with Algorithm 1, using a uniform non-fitted mesh of size varying in $h \in \{0.1, 0.008, 0.06, 0.04, 0.02\}$, see Figure 5c for the case $h = 0.06$. Figure 5d reports the results obtained with Algorithm 2, using a fitted mesh refined around the beam interface. The $L^2$ and $H^1$ numerical errors are reported in Figure 6. We observe a convergence rate around $1/2$ which corresponds to the classical convergence rate for the fictitious domain method using finite element approximation, see e.g. [71].

### 5.1.2. Impact of fluid mesh refinement

We now consider the full unsteady problem with a fluid domain of size $3 \times 1 \times 1$, where the following time-dependent parabolic velocity profile is enforced at the inlet boundary $\Gamma_{in}$:

$$ \mathbf{u}_{\text{in}}(\mathbf{x}, t) \stackrel{\text{def}}{=} u_{\text{ref}} \left( 1 + \cos \left( \frac{2\pi t}{T} \right) \right) (1 - (z-1)(z+1)) e_x, \quad (5.1) $$

with final time $T=0.06$. The following beam parameters are considered: Young's modulus $E=10^7$, radius $\epsilon=0.06$ and density $\rho^b = 1$. For the fluid, we take a density $\rho^f=1$. We run simulations for different values of $h/\epsilon$, representing the resolution of the beam interface by the fluid mesh. Specifically, we consider two types of fluid meshes, see Figure 4. The first one is refined in a region near the beam interface, with $h/\epsilon=1/4$, and we have $h=0.05$ on the faces of the rectangular domain. The second type correspond to
<figure><img src="image_7.png" /><figcaption>(a) Algorithm 1 with $\varepsilon=0.06$ and $h=0.006$ near the beam interface</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Algorithm 2 with $\varepsilon=0.06$ and $h=0.006$ near the beam interface</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Algorithm 1 with $\varepsilon=0.1$ and $h=0.06$ applied uniformly</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Algorithm 2 with $\varepsilon=0.1$ and $h=0.01$ near the beam interface</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 6: Convergence of the reduced fluid model for $\varepsilon \to 0$ and $h \to 0$.</figcaption></figure>

a uniform mesh in the whole fluid domain, with $h=0.05$. Note that for Algorithm 2, which serves as reference in the comparisons, we only consider the refined fluid mesh fitted to the interface. For the beam and the Lagrange multiplier space, we discretize the interval $(0, L)$ with a uniform mesh such that $H=0.05$. The behavior of the coupled system is studied over the time interval $[0, 0.06]$ with a time-step length of $\tau=5 \cdot 10^{-3}$.

Case $h/\varepsilon \ll 1$. We consider here the refined fluid mesh. In Figure 7, we report the time history of the displacement of the last cross-section of the beam obtained with Algorithms 1, 2 and 3 for different values of viscosity $\mu^f$ and of the reference velocity $u_{\text{ref}}$, deliberately calibrated to keep comparable displacement amplitudes across the different test cases. In Table 1, we also provide the relative error in the $L^\infty$ norm with respect to the solution provided by Algorithm 2. Lastly in Figure 8, we present snapshots of the
fluid velocity magnitude and pressure fields at the final time $t = 0.06$.

<figure><img src="image_12.png" /><figcaption>Figure 7: Time history of the displacement of the last cross-section of the beam along the x-axis for several fluid viscosities $\mu^f$ and reference inlet velocities $u_{ref}$, with $h/\varepsilon \ll 1$.</figcaption></figure>

<table><caption>Table 1: Relative error $\| (\mathbf{r}(L, \cdot) - \mathbf{r}^{\text{ref}}(L, \cdot)) \cdot e_x \|_{\infty} / \| \mathbf{r}^{\text{ref}}(L, \cdot) \cdot e_x \|_{\infty}$ where $\mathbf{r}^{\text{ref}}(L, \cdot)$ is the displacement of the last cross-section of the beam computed with Algorithm 2, with $h/\varepsilon \ll 1$.</caption><thead><tr><th></th><th>μ<sup>f</sup>=1, u<sub>ref</sub>=10.</th><th>μ<sup>f</sup>=0.1, u<sub>ref</sub>=50.</th><th>μ<sup>f</sup>=0.05, u<sub>ref</sub>=100.</th></tr></thead><tbody><tr><td>Algorithm 1, N=0</td><td>2.84 · 10<sup>-2</sup></td><td>2.88 · 10<sup>-2</sup></td><td>3.01 · 10<sup>-2</sup></td></tr><tr><td>Algorithm 3</td><td>4.94 · 10<sup>-1</sup></td><td>5.46 · 10<sup>-1</sup></td><td>5.78 · 10<sup>-1</sup></td></tr></tbody></table>

We can observe that the displacement, velocity and pressure fields obtained with Algorithm 1 provide a very good approximation of the reference given by Algorithm 2. In this context, taking $N=0$ seems to be sufficient for accurately capturing both the beam and fluid dynamics. This can be attributed to the fact that the numerical solution of Algorithm 2 presents a symmetric flow pattern around the beam, which results in a small torque acting onto beam and a subsequent low angular velocity. Additionally, the force transmitted to the beam and the subsequent beam displacement velocity for Algorithm 2 belongs to the functional space $\mathcal{F}_0$.

In contrast, Algorithm 3 does not provide a satisfactory numerical solution. The area of influence of the beam on the fluid appears considerably smaller compared to the solution computed with Algorithm 2. The lower magnitude of the velocity and pressure is noticeable in Figure 8, which is also consistent with differences in displacement observed in Figure 7 and Table 1. As discussed in [33], the extent of the influence of the beam on the surrounding fluid for Algorithm 3 is determined by the size of the fluid elements. When $h/\varepsilon \ll 1$, this difference becomes important and then leads to a considerable deviation from the reference solution.

Case $h/\varepsilon \approx 1$. We now consider the uniform fluid mesh. In this configuration, the thickness of the beam ($\varepsilon=0.06$) is roughly equal to the size of the fluid elements. As indicated above, we keep a refined mesh for Algorithm 2, which serves as reference solution.

<table><caption>Table 2: Relative error $\| (\mathbf{r}(L, \cdot) - \mathbf{r}^{\text{ref}}(L, \cdot)) \cdot e_x \|_{\infty} / \| \mathbf{r}^{\text{ref}}(L, \cdot) \cdot e_x \|_{\infty}$ where $\mathbf{r}^{\text{ref}}(L, \cdot)$ is the displacement of the last cross-section of the beam computed with Algorithm 2, with $h/\varepsilon \approx 1$.</caption><thead><tr><th></th><th>μ<sup>f</sup>=1, u<sub>ref</sub>=10.</th><th>μ<sup>f</sup>=0.1, u<sub>ref</sub>=50.</th><th>μ = 0.05, u<sub>ref</sub> = 100.</th></tr></thead><tbody><tr><td>Algorithm 1, N=0</td><td>6.88 · 10<sup>-2</sup></td><td>8.14 · 10<sup>-2</sup></td><td>9.39 · 10<sup>-2</sup></td></tr><tr><td>Algorithm 3</td><td>2.63 · 10<sup>-1</sup></td><td>2.71 · 10<sup>-1</sup></td><td>2.85 · 10<sup>-1</sup></td></tr></tbody></table>

The results for the displacement of the last beam cross-section are given in Figure 9, relative error in $L^\infty$ norm with respect to the reference solution are presented in Table 2 and snapshots of the velocity magnitude and pressure at time $t=0.06$ are provided in Figure 10. We can observe that Algorithm 1 is able to capture the dynamics of the beam displacement and the main fluid velocity/pressure patterns. However, the comparison with Algorithm 2 degrades a bit with respect to the previous case ($h/\varepsilon \ll 1$). This is visible in the slight deviations observed in the displacement history and in the pressure pattern. As regards Algorithm 3, provides a better approximation of the beam displacement than in the scenario $h/\varepsilon \ll 1$. However, the solution is still less accurate than the one obtained with Algorithm 1. Furthermore, although an improvement can be seen in the solution for beam displacement, Figure 9 shows that the
<figure><img src="image_13.png" /><figcaption>Figure 8: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time $t=0.06$ obtained with Algorithms 1–3, with $h/\varepsilon \ll 1$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9: Time history of the displacement of the last cross-section of the beam along the x-axis for various fluid viscosity $\mu^f$ and reference inlet velocity $u_{ref}$, with $h/\varepsilon \approx 1$.</figcaption></figure>

solution for fluid velocity suffers from the same limitations as in the case $h/\varepsilon \ll 1$. In general, with a coarser mesh, we observe, however, that the pressure jump for both Algorithm 1 and 3 is smaller than the one obtained with Algorithm 2. This can be attributed to the singularity in the pressure, whose approximation using continuous piece-wise affine functions requires a fine mesh to be captured accurately.

Case $h/\varepsilon \gg 1$. Here, we still consider the uniform mesh. However, we decrease the value of $\varepsilon$, in particular, we consider $\varepsilon \in \{0.06, 0.04, 0.02\}$. In Algorithm 2, as before, the mesh is refined around the beam to ensure that $h/\varepsilon=1/4$, which resolves the interface regardless of the beam thickness $\varepsilon$.
<figure><img src="image_15.png" /><figcaption>Figure 10: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time $t=0.06$ obtained with Algorithms 1-3, with $h/\epsilon \approx 1$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 11: Time history of the displacement of the last cross-section of the beam along the x-axis for various fluid viscosity $\mu^f$ and reference inlet velocity $u_{\text{ref}} - h/\epsilon \gg 1$.</figcaption></figure>

<table><caption>Table 3: Relative error $\| (\mathbf{r}(L, \cdot) - \mathbf{r}^{\text{ref}}(L, \cdot)) \cdot e_x \|_{\infty} / \| \mathbf{r}^{\text{ref}}(L, \cdot) \cdot e_x \|_{\infty}$ where $\mathbf{r}^{\text{ref}}(L, \cdot)$ is the displacement of the last cross-section of the beam computed with Algorithm 2, with $h/\epsilon \gg 1$.</caption><thead><tr><th></th><th>&epsilon;=0.06, E=10<sup>7</sup>.</th><th>&epsilon;=0.04, E=4 &middot; 10<sup>7</sup>.</th><th>&epsilon;=0.02, E=4 &middot; 10<sup>8</sup>.</th></tr></thead><tbody><tr><td>Algorithm 1, N=0</td><td>9.39 &middot; 10<sup>-2</sup></td><td>1.13 &middot; 10<sup>-1</sup></td><td>2.20 &middot; 10<sup>-1</sup></td></tr><tr><td>Algorithm 3</td><td>2.85 &middot; 10<sup>-1</sup></td><td>1.45 &middot; 10<sup>-1</sup></td><td>2.03 &middot; 10<sup>-1</sup></td></tr></tbody></table>

To keep comparable amplitudes for the beam displacement independently of the beam thickness, the
<figure><img src="image_17.png" /><figcaption>Figure 12: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time $t=0.06$ obtained with Algorithms 1–3, with $h/\varepsilon \gg 1$.</figcaption></figure>

Young's modulus of the beam is modified accordingly with values $E \in \{10^7, 4 \cdot 10^7, 4 \cdot 10^8\}$. We fix $\mu^f=0.05$ and $u_{\text{ref}}=100$. The results obtained are reported in Figure 11 for the beam displacement, in Table 3 for the relative error in $L^\infty$ norm and in Figure 12 for the velocity and pressure fields. We can clearly observe that by reducing $\varepsilon$, the numerical approximations of the beam displacement, the fluid velocity and the fluid pressure provided by Algorithms 1 and 3 get closer. The price to pay is a significant degradation of the accuracy of Algorithm 1 with respect to the previous cases.

In conclusion, as the ratio $h/\varepsilon$ increases, the resolution of the mesh becomes insufficient to accurately capture the fluid-solid interface. In this context, the numerical comparisons indicate that imposing the coupling conditions on the centerline delivers similar accuracy to enforcing them on interface. However, given the relatively large mesh size compared to the beam thickness, delivering an accurate solution for the fluid pressure and velocity on the interface becomes challenging. The primary advantage of the Algorithm 1 lies in its ability to capture the dynamics of the problem effectively with a reasonably refined fluid mesh near the interface, while maintaining reduced interface conditions in the spirit of Algorithm 2.

### 5.1.3. Influence of the number of Fourier modes

In this section, we consider a variant of the previous numerical example in which the beam is immersed in a shear flow, with the purpose of inducing torsional effects on the beam and illustrate the importance, in this new context, of considering enough modes ($N \ge 1$) in Algorithm 1. We consider a fluid domain of size $3 \times 1 \times 1$ where we enforce homogeneous Neumann boundary conditions on $\Gamma_{\text{in}}$ and $\Gamma_{\text{out}}$, homogeneous Dirichlet boundary conditions on the lower face of the domain, and $u=u_{\text{ref}}(y-0.5)(1-(x-1)(x+1))e_x$
<figure><img src="image_18.png" /><figcaption>Figure 13: Rotation of the last cross-section along the z-axis of the beam immersed in a shear flow.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 14: Snapshots of the beam deformation and the velocity magnitude induced by a shear flow, at time $t \in \{0, 0.02, 0.04, 0.06\}$ (from top to bottom).</figcaption></figure>

on the three other lateral faces, with $u_{\text{ref}}=40$. In order to further induce torsional effects in the beam, we consider elliptical cross-sections with an aspect ratio $a=2$ and a small axis $\epsilon=0.04$. The physical parameters are $\mu^f=1$, $E=4 \cdot 10^5$ and $\rho^b=1$.
In Algorithms 1 and 2, we consider the same type of fluid meshes, refined around the beam ($h/\epsilon=1/4$) and uniform on the faces of the rectangular domain ($h=0.05$). We compare the approximations provided by Algorithm 1 and Algorithm 2 with $N=0, 1$. In Figure 13, we report the time history of the rotation of the last cross-section around the $z$ axis over the time interval $[0, 0.06]$. In Figure 14, we provide some snapshots of the beam deformed configuration at time $t \in \{0, 0.02, 0.04, 0.06\}$. We can clearly observe that Algorithm 1 with $N=0$ cannot capture the torsion induced by the shear fluid flow on the beam. Indeed, since the Lagrange multiplier is constant on each cross-section, no torque is transmitted to the beam. As a result, since the beam is initially at rest, it remains steady throughout the simulation. In contrast, Algorithm 1 with $N=1$ is able to reproduce the rotational motion around the $z$-axis obtained with Algorithm 2. Similarly, the velocity of the beam interface, mainly determined by its torsion velocity around the $z$-axis, is well captured, as shown in Figure 14. In terms of displacement, this favorable outcome comes from the fact that the moment transmitted to the structure, as well as the rotation velocity of the beam interface, given in equations (2.7)-(2.8), belong to the space $F_1$.

## 5.2. Single beam immersed in a Navier-Stokes flow

In this section, we keep the geometric configuration considered in Section 5.1, where the fluid domain parameters are given by $l_x=3$, $l_y=1$ and $l_z=1$, but now with a beam immersed in a Navier-Stokes flow, so that the Reynolds number is not zero anymore. The fluid velocity is initialized with the solution of the Stokes equations, in which the velocity on the beam is set to zero. We consider a fluid density $\rho^f=1$ and a beam density $\rho^b=1$. The thickness of the beam cross-section is set to $\epsilon=0.06$ and the Reynolds number $Re \in \{2.4, 120, 480\}$, with $\mu \in \{1, 0.1, 0.05\}$ and $u_{\text{ref}} \in \{10, 50, 100\}$, respectively. To keep similar displacement amplitudes when varying the Reynolds number, the Young modulus of the beam is chosen as $E \in \{10^7, 1.6 \cdot 10^8, 3 \cdot 10^8\}$. We use a mesh that accurately resolves the beam interface such that $h/\epsilon=1/4$ in a region near the beam and $h = 0.05$ on the faces of the rectangular domain.

<figure><img src="image_20.png" /><figcaption>Figure 15: Time history of the displacement of the last cross-section of the beam along the $x$-axis for various Reynolds number $Re$ and Young's modulus $E$.</figcaption></figure>

<table><caption>Table 4: Relative error $\| (\mathbf{r}(L, \cdot) - \mathbf{r}^{\text{ref}}(L, \cdot)) \cdot e_x \|_{\infty} / \| \mathbf{r}^{\text{ref}}(L, \cdot) \cdot e_x \|_{\infty}$ where $\mathbf{r}^{\text{ref}}(L, \cdot)$ is the displacement of the last cross-section of the beam computed with Algorithm 2 or Algorithm 2-Ω.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Re=2.4, E=4 · 10<sup>7</sup>.</th><th colspan="2">Re=120, E=1.6 · 10<sup>8</sup>.</th><th colspan="2">Re=480, E=3 · 10<sup>8</sup>.</th></tr><tr><th>Algorithm 2</th><th>Algorithm 2-Ω</th><th>Algorithm 2</th><th>Algorithm 2-Ω</th><th>Algorithm 2</th><th>Algorithm 2-Ω</th></tr></thead><tbody><tr><th>Algorithm 1, N=0</th><td>7.05 · 10<sup>-2</sup></td><td>6.49 · 10<sup>-2</sup></td><td>5.17 · 10<sup>-1</sup></td><td>3.41 · 10<sup>-1</sup></td><td>4.31 · 10<sup>-1</sup></td><td>2.90 · 10<sup>-1</sup></td></tr><tr><th>Algorithm 1, N=1</th><td>6.95 · 10<sup>-2</sup></td><td>6.38 · 10<sup>-2</sup></td><td>2.67 · 10<sup>-1</sup></td><td>1.27 · 10<sup>-1</sup></td><td>2.06 · 10<sup>-1</sup></td><td>8.75 · 10<sup>-2</sup></td></tr><tr><th>Algorithm 1, N=2</th><td>7.89 · 10<sup>-2</sup></td><td>7.30 · 10<sup>-2</sup></td><td>2.02 · 10<sup>-1</sup></td><td>7.66 · 10<sup>-2</sup></td><td>2.00 · 10<sup>-1</sup></td><td>8.43 · 10<sup>-2</sup></td></tr></tbody></table>

**Impact of the Reynolds number.** We first investigate the influence of the Reynolds number on the accuracy of Algorithm 1. In Figure 15, we report the time history of the displacement of the last cross-section of the beam with $Re \in \{2.4, 120, 480\}$ and $E \in \{4 \cdot 10^7, 4 \cdot 10^7, 4 \cdot 10^8\}$, and in Table 4 the $L^\infty$ relative error with respect to Algorithm 2. Note that in Figure 15, we introduced an additional algorithm, referred to as Algorithm 2-Ω in which the fluid is solved in the whole domain $\Omega$. The motivation of this algorithm will be discussed in the next paragraph. We observe that the higher the Reynolds number the less accurately the displacement of the beam is captured by Algorithm 1. On the contrary, the accuracy of the results is improved by increasing the number of Fourier modes $N$. In particular, the numerical solution appears to converge when $N$ increases. Indeed, the accuracy gap between $N=0$ and $N=1$ is quite significant, but
decreases significantly between $N=1$ and $N=2$. However, the numerical results indicate that the limit differs from the numerical solution provided by Algorithm 2, notably as the Reynolds number increases. One possible explanation of this mismatch is that, in the fictitious domain approach, the fluid equations are solved in the entire domain $\Omega$, which introduces a non-physical flow in the part occupied by the beam with the subsequent additional artificial stress on the interface (see Section 2.3). It is worth noting that, in the case of a Stokes flow, the numerical results of Section 5.1 indicate that the impact of such spurious stress on the accuracy of Algorithm 1 is rather limited.

<figure><img src="image_21.png" /><figcaption>Figure 16: Time history of the displacement of the last cross-section of the beam along the x-axis for various beam thickness $\epsilon$.</figcaption></figure>

<table><caption>Table 5: Relative error $\| (\mathbf{r}(L, \cdot) - \mathbf{r}^{\text{ref}}(L, \cdot)) \cdot e_x \|_{\infty} / \| \mathbf{r}^{\text{ref}}(L, \cdot) \cdot e_x \|_{\infty}$ where $\mathbf{r}^{\text{ref}}(L, \cdot)$ is the displacement of the last cross-section of the beam computed with Algorithm 2.</caption><thead><tr><th></th><th>Re=800, ε=0.1, E=10<sup>8</sup>.</th><th>Re=480, ε=0.06, E=1.2 · 10<sup>9</sup>.</th><th>Re=160, ε=0.02, E=1.2 · 10<sup>10</sup>.</th></tr></thead><tbody><tr><td>Algorithm 1, N=0</td><td>6.22 · 10<sup>-1</sup></td><td>4.31 · 10<sup>-1</sup></td><td>3.23 · 10<sup>-1</sup></td></tr><tr><td>Algorithm 1, N=1</td><td>2.61 · 10<sup>-1</sup></td><td>2.06 · 10<sup>-1</sup></td><td>1.96 · 10<sup>-1</sup></td></tr><tr><td>Algorithm 1, N=2</td><td>2.54 · 10<sup>-1</sup></td><td>2.00 · 10<sup>-1</sup></td><td>1.63 · 10<sup>-1</sup></td></tr></tbody></table>

*Impact of the fictitious domain method.* In order to investigate the impact of the spurious internal fluid on the accuracy of Algorithm 1, we consider a variant of Algorithm 2 in which the fluid is solved in the whole domain $\Omega$, so that the portion of the domain occupied by the beam is also filled with the fluid. This variant will be termed Algorithm 2-$\Omega$. Figure 15 and Table 4 indicate that the displacement of the last cross-section of the beam obtained with Algorithm 1, with $N \in \{1, 2\}$, is close to the one provided by Algorithm 2-$\Omega$. Since the volume of the region occupied by the beam is proportional to $\epsilon^2$, one would expect that the perturbation induced by the artificial fluid decreases as $\epsilon$ tends to zero. We propose to exhibit this phenomenon by gradually decreasing the thickness of the beam, by taking $\epsilon \in \{0.1, 0.04, 0.02\}$ while keeping the same fluid parameters as in the case with $Re=480$ and $\epsilon=0.06$. To keep similar amplitudes of the beam deflection, we take the Young's modulus $E \in \{10^8, 1.2 \cdot 10^9, 1.2 \cdot 10^{10}\}$. In Figure 16 we report the time history of the displacement of the last cross-section of the beam and in Table 5 the $L^\infty$ relative error with respect to Algorithm 2. We observe that the displacement of the last beam cross-section approaches the reference solution as the thickness decreases.

Another difficulty of standard fictitious domain methods is the approximation of the pressure jump across the beam interface via continuous functions. In Figure 17, we provide snapshots of the magnitude of the velocity and the pressure fields obtained with $Re=480$ and a constant velocity inlet $u_{\text{in}} = 2u_{\text{ref}}e_x$. We can clearly observe that the pressure jump provided by Algorithm 2-$\Omega$ is smaller than the one obtained with Algorithm 2. We can also notice that the continuous pressure approximation does not ensure mass conservation across the interface, as a spurious flow is generated inside the region occupied by the beam. As a result, the vortices that form behind the beam are also much weaker. Although increasing the number of modes in Algorithm 1 significantly reduces the amount of spurious fluid, notably for $N=2$, the obtained numerical solutions is much closer to the one provided by Algorithm 2-$\Omega$ than by Algorithm 2.

Some solutions exist to circumvent these issues. In order to correctly evaluate the fluid stress at the interface, one can force the fluid velocity inside the fictitious region to be equal to the solid velocity using, for instance, distributed Lagrange multipliers [72, 40]. However, this approach often requires the introduction of 3D corrections, which jeopardize the reduced nature of the present methodology. Thus, we have chosen to not include them into the present study.
In order to mitigate the impact of the spurious stress across the interface, we consider the approach which consists in introducing a scalar Lagrange multiplier to guarantee mass conservation within the domain $\Omega_{\epsilon}^f$ (see [73, 74]), viz.,

$$ \int_{\partial \Omega_{\epsilon}^f} u \cdot n = 0. $$

We denote by Algorithm 1 - MC the solution procedure obtained by combining Algorithm 1 with this approach. Figure 17 reports the snapshots of the fluid velocity magnitude and pressure obtained with Algorithms 1 and 2, Algorithm 2-Ω and Algorithm 1-MC at time $t=0.06$. We can observe that the introduction of the scalar mass conservation constraint noticeably reduces the fluid flow both inside and across the beam. Furthermore, the pressure and field around the beam are rendered with greater fidelity, particularly near the beam interface, where we observe an increased upward velocity along the interface, together with vortices forming much closer to the beam.

<figure><img src="image_22.png" /></figure>
<figure><img src="image_23.png" /><figcaption>Figure 17: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time $t=0.06$ with Algorithm 2, Algorithm 2-$\Omega$, Algorithm 1 and Algorithm 1-MC, with $Re=480$.</figcaption></figure>

### 5.3. Multiple-beams immersed in incompressible viscous flow

In order to further assess the robustness of the proposed approach, we consider a more complex scenario involving a swarm of beams made of 75 circular beams, arranged with 5 along the y-axis and 15 along the x-axis, in a fluid domain of size $3 \times 1 \times 1$. The beams have a radius $\epsilon=0.05$ and are spaced 0.05 apart from each other. The boundary conditions are similar to the one described in the previous section with a time dependent inlet velocity (5.1). The beam is modeled by a non-linear co-rotational formulation that allows for large deformations (see [62, Chapter 5]). In this formulation, the total motion of the beam is divided between the motion of the element-based local co-rotating frame of reference and the deformation of the element within this local frame. For the element's deformation within this local frame, the generalized Timoshenko linear beam theory is used as described in Section 2.1. Given the proximity of the beams, they are likely to come in contact. To address this issue, we have incorporated a penalty contact algorithm using a raytracing projection on the 2D beam interface described in Appendix A. All the beams have the same Young's modulus and density, $E=10^8$ and $\rho^b=1$, respectively.

We first consider a Stokes flow with viscosity $\mu^f=1$ and reference inlet velocity $u_{ref}=10$. The size of the mesh is uniform throughout the fluid domain, with $h=0.05$. In Figure 18, we provide snapshots of the fluid velocity magnitude computed with Algorithm 1 at times $t \in \{0.05, 0.1, 0.2, 0.3, 0.45, 0.6\}$. We observe a collective behavior of the beams oscillating from left to right. In particular, they do not get into contact. We then consider a Navier-Stokes flow, with Reynolds number $Re=480$, computed with $\mu^f=0.05$, $u_{ref}=100$ and $\rho^f=1$. The fluid mesh is kept uniform with a smaller mesh size $h=0.035$. Snapshots of the fluid velocity magnitude calculated with Algorithm 1 with $N=2$ are provided in Figure 19. We observe a significant reduction in fluid velocity due to the presence of beams. Additionally, a distinct wave flow pattern appears in the upper region. We also observe orthogonal waves to the inlet flow pushing the beams outward, thereby increasing the fluid velocity within these regions of the domain.

## 6. Conclusions

We have proposed a new approach to modeling and simulation of slender structures immersed in a 3D flow. A salient feature of the proposed methodology is that the 3D fluid model and the reduced 1D solid model are coupled via well-posed trace operators of co-dimension two, based on a suitable Fourier projection of the kinematic constraint on the 2D interface. The method has been extensively compared
<figure><img src="image_24.png" /><figcaption>Figure 18: Snapshots of the fluid velocity magnitude for multiple beams immersed in a Stokes flow with Algorithm 1 and $N=0$ at time $t \in \{0.005, 0.01, 0.02, 0.03, 0.045, 0.06\}$ (note that at time $t=0.03$, the inlet velocity is null).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 19: Snapshots of the fluid magnitude fluid for multiple beams immersed in a Navier-Stokes flow with Algorithm 1 and $N=2$ at time $t \in \{0.005, 0.01, 0.02, 0.03, 0.045, 0.06\}$ (note that at time $t=0.03$, the inlet velocity is null).</figcaption></figure>

with a fully resolved (ALE based) approach and an alternative reduced method. Numerical evidence indicates that at low Reynolds numbers, the method is accurate regardless of the ratio between the size of the background mesh and the radius of the beam. For moderate Reynolds numbers, accurately capturing the fluid flow around the two-dimensional beam interface requires both mesh refinement and a higher number of Fourier modes.

The numerical comparisons also indicated that for high Reynolds numbers, most of the numerical errors arise from the fictitious domain method. As a result, we have proposed strategies to address these issues without resorting to 3D techniques, to guarantee consistency with the underlying reduced nature of the method. Furthermore, numerical evidence has been provided on the fact that the choice of the approximated Fourier space can have a significant influence on the dynamic of the beams, including whether it is possible to capture the rotation of the cross-sections. Finally, to assert the robustness of the method and its potential integration with more general computational frameworks for fluids and solids, we have tested it considering a significant number of beams that can get into contact.
Current investigations are devoted to the design and analysis of stable and accurate loosely coupled schemes. By implementing such schemes, we could strike a balance between computational efficiency and accuracy.

Acknowledgement

Fabien Lespagnol is supported by the project A new computational approach for the fluid-structure interaction of slender bodies immersed in three-dimensional flows granted by the Università Italo-Francese, in the framework Vinci 2019. Fabien Lespagnol and Paolo Zunino acknowledge the partial support of the grant MUR PRIN 2022 No. 2022WKWZA8 Immersed methods for multiscale and multiphysics problems (IMMEDIATE). The present research is part of the activities of project Dipartimento di Eccellenza 2023-2027, funded by MUR. Paolo Zunino is a member of the Gruppo Nazionale per il Calcolo Scientifico (GNCS), Istituto Nazionale di Alta Matematica (INdAM).

Appendix A. Contact algorithm used in Section 5.3

<figure><img src="image_26.png" /><figcaption>Figure A.20: Ray-tracing strategy for contact detection.</figcaption></figure>

In order to model contact between solids, we include an additional term into Algorithm 1. This contact term is based on a ray-tracing strategy with penalization on the reconstructed interface of the beams (see, e.g., [75, Chapter 11]). Since the balance equations are defined on the beam centerline, the penalization term is then integrated on each cross section.

More precisely, let consider $\hat{\omega}_1^\epsilon$ and $\hat{\omega}_2^\epsilon$ denote the reference configuration of two distinct beams, such that $\hat{\Sigma}_1^\epsilon \stackrel{\text{def}}{=} \partial\hat{\omega}_1^\epsilon \setminus \partial\Omega$ is the master surface and $\hat{\Sigma}_2^\epsilon \stackrel{\text{def}}{=} \partial\hat{\omega}_2^\epsilon \setminus \partial\Omega$ is the slave surface. For all $\boldsymbol{x} \in \Sigma_2^{\epsilon,n} \stackrel{\text{def}}{=} \phi^n(\hat{\Sigma}_2^\epsilon)$, denoting by $\boldsymbol{n}_\boldsymbol{x}$ the outer normal of the surface $\Sigma_2^{\epsilon,n} \stackrel{\text{def}}{=} \phi^n(\hat{\Sigma}_2^\epsilon)$ at point $\boldsymbol{x}$, we denote by $\Pi_R(\boldsymbol{x})$ the closest intersection of the master surface $\Sigma_1^{\epsilon,n}$ with the line passing through $\boldsymbol{x}$ and having direction $\boldsymbol{n}_\boldsymbol{x}$. This strategy can be referred to as ray-tracing (see Figure A.20). We then introduce a gap function denoted by $g$ which satisfies

$$g(\boldsymbol{x}) = (\Pi_R(\boldsymbol{x}) - \boldsymbol{x}) \cdot \boldsymbol{n}_\boldsymbol{x} = (\Pi_R(\phi_H^n(\hat{\boldsymbol{x}})) - \phi_H^n(\hat{\boldsymbol{x}})) \cdot \boldsymbol{n}_\boldsymbol{x}$$

where $\hat{\boldsymbol{x}} \stackrel{\text{def}}{=} (\phi_H^n)^{-1}(\boldsymbol{x})$. The non-penetration conditions writes $g \ge 0$ and is enforced through penalization. The contact term for the beam $\hat{\omega}_1^\epsilon$ is then given by

$$\frac{\gamma}{2h^2} \int_{\hat{\Sigma}_2^\epsilon} [(\Pi_R(\phi_H^n(\hat{\boldsymbol{x}})) - \phi_H^n(\hat{\boldsymbol{x}})) \cdot \boldsymbol{n}_\boldsymbol{x}]_+ (\delta\overline{\boldsymbol{r}}_H(\hat{\boldsymbol{x}}) - \delta\overline{\boldsymbol{r}}_H(\hat{\boldsymbol{y}})) \cdot \boldsymbol{n}_\boldsymbol{x} \quad \forall \delta\overline{\boldsymbol{r}}_H \in \mathbf{Y}_H,$$

where $\gamma > 0$ is user-defined parameter and $[x]_+ \stackrel{\text{def}}{=} \max\{0, x\}$.
References

[1] Y. Fan, J. Xia, Simulation of 3d parachute fluid-structure interaction based on nonlinear finite element method and preconditioning finite volume method, chinese Journal of Aeronautics 27 (6) (2014) 1373-1383.

[2] Z. Huang, P. Avery, C. Farhat, J. Rabinovitch, A. Derkevorkian, L. D. Peterson, Simulation of parachute inflation dynamics using an eulerian computational framework for fluid-structure interfaces evolving in high-speed turbulent flows, in: 2018 AIAA Aerospace Sciences Meeting, 2018, p. 1540.

[3] Y. Kim, C. S. Peskin, 3-d parachute simulation by the immersed boundary method, Computers & Fluids 38 (6) (2009) 1080-1090.

[4] S. Sathe, R. Benney, R. Charles, E. Doucette, J. Miletti, M. Senga, K. Stein, T. Tezduyar, Fluid-structure interaction modeling of complex parachute designs with the space-time finite element techniques, Computers & Fluids 36 (1) (2007) 127-135.

[5] D. W. Peaceman, Interpretation of well-block pressures in numerical reservoir simulation (includes associated paper 6988), Society of Petroleum Engineers Journal 18 (03) (1978) 183-194.

[6] D. W. Peaceman, Interpretation of well-block pressures in numerical reservoir simulation with non-square grid blocks and anisotropic permeability, Society of Petroleum Engineers Journal 23 (03) (1983) 531-543.

[7] D. Cerroni, F. Laurino, P. Zunino, Mathematical analysis, finite element approximation and numerical solvers for the interaction of 3d reservoirs with 1d wells, GEM-International Journal on Geomathematics 10 (2019) 1-27.

[8] M. Wang, E. J. Avital, X. Bai, C. Ji, D. Xu, J. J. Williams, A. Munjiza, Fluid-structure interaction of flexible submerged vegetation stems and kinetic turbine blades, Computational Particle Mechanics 7 (2020) 839-848.

[9] M. Luhar, J. Rominger, H. Nepf, Interaction between flow, transport and vegetation spatial structure, Environmental Fluid Mechanics 8 (2008) 423-439.

[10] S. A. Mattis, C. N. Dawson, C. E. Kees, M. W. Farthing, An immersed structure approach for fluid-vegetation interaction, Advances in Water Resources 80 (2015) 1-16.

[11] H. Chen, Q.-P. Zou, Eulerian-lagrangian flow-vegetation interaction model using immersed boundary method and openfoam, Advances in Water Resources 126 (2019) 176-192.

[12] J. Tambača, M. Kosor, S. Čanić, D. P. MD, Mathematical modeling of vascular stents, SIAM journal on applied mathematics 70 (6) (2010) 1922-1952.

[13] J. Tambača, S. Čanić, D. Paniagua, A novel approach to modeling coronary stents using a slender curved rod model: a comparison between fractured xience-like and palmaz-like stents, Applied and Numerical Partial Differential Equations: Scientific Computing in Simulation, Optimization and Control in a Multidisciplinary Context (2010) 41-58.

[14] P. Zunino, J. Tambača, E. Cutrì, S. Čanić, L. Formaggia, F. Migliavacca, Integrated stent models based on dimension reduction: review and future perspectives, Annals of biomedical engineering 44 (2016) 604-617.

[15] A. Decoene, S. Martin, F. Vergnet, A continuum active structure model for the interaction of cilia with a viscous fluid, Zeitschrift für Angewandte Mathematik und Mechanik (2023) e202100534doi: 10.1002/zamm.202100534.

[16] F. Alouges, A. DeSimone, L. Giraldi, M. Zoppello, Self-propulsion of slender micro-swimmers by curvature control: N-link swimmers, International Journal of Non-Linear Mechanics 56 (2013) 132-141, soft Matter: a nonlinear continuum mechanics perspective. doi:https://doi.org/10.1016/j.ijnonlinmec.2013.04.012.
[17] L. Giraldi, P. Martinon, M. Zoppello, Controllability and optimal strokes for n-link microswimmer, in: 52nd IEEE Conference on Decision and Control, IEEE, 2013, pp. 3870–3875.

[18] E. Lauga, T. R. Powers, The hydrodynamics of swimming microorganisms, Reports on progress in physics 72 (9) (2009) 096601.

[19] J. C. Simo, A finite strain beam formulation. the three-dimensional dynamic problem. part i, Computer methods in applied mechanics and engineering 49 (1) (1985) 55–70.

[20] E. Reissner, On one-dimensional finite-strain beam theory: the plane problem, Zeitschrift für angewandte Mathematik und Physik ZAMP 23 (5) (1972) 795–804.

[21] C. Meier, A. Popp, W. A. Wall, An objective 3d large deformation finite element formulation for geometrically exact curved kirchhoff rods, Computer Methods in Applied Mechanics and Engineering 278 (2014) 445–478.

[22] C. Meier, A. Popp, W. A. Wall, A locking-free finite element formulation and reduced models for geometrically exact kirchhoff rods, Computer Methods in Applied Mechanics and Engineering 290 (2015) 314–341.

[23] C. Meier, M. J. Grill, W. A. Wall, A. Popp, Geometrically exact beam elements and smooth contact schemes for the modeling of fiber-based materials and structures, International Journal of Solids and Structures 154 (2018) 124–146.

[24] F. P. Baaijens, A fictitious domain/mortar element method for fluid-structure interaction, International Journal for Numerical Methods in Fluids 35 (7) (2001) 743–761.

[25] L. Boilevin-Kayl, M. Fernández, J.-F. Gerbeau, Numerical methods for immersed FSI with thin-walled structures, Comput. & Fluids 179 (2019) 744–763. doi:10.1016/j.compfluid.2018.05.024.

[26] D. Chapelle, A. Ferent, Modeling of the inclusion of a reinforcing sheet within a 3D medium, Math. Models Methods Appl. Sci. 13 (4) (2003) 573–595. doi:10.1142/S0218202503002635.

[27] M. Landajuela, M. Vidrascu, D. Chapelle, M. Fernández, Coupling schemes for the FSI forward prediction challenge: comparative study and validation, Int. J. Numer. Methods Biomed. Eng. 33 (4) (2017) e02813, 23.

[28] R. Maniyeri, S. Kang, Numerical study on the rotation of an elastic rod in a viscous fluid using an immersed boundary method, Journal of Mechanical Science and Technology 26 (5) (2012) 1515–1522.

[29] D. Huang, P. Avery, C. Farhat, An embedded boundary approach for resolving the contribution of cable subsystems to fully coupled fluid-structure interaction, Internat. J. Numer. Methods Engrg. 122 (19) (2021) 5409–5429.

[30] R. F. Ausas, C. G. Gebhardt, G. C. Buscaglia, A finite element method for simulating soft active non-shearable rods immersed in generalized newtonian fluids, Communications in Nonlinear Science and Numerical Simulation 108 (2022) 106213.

[31] B. E. Griffith, S. Lim, Simulating an elastic ring with bend and twist by an adaptive generalized immersed boundary method, Communications in Computational Physics 12 (2) (2012) 433–461.

[32] J. K. Wiens, J. M. Stockie, Simulating flexible fiber suspensions using a scalable immersed boundary algorithm, Computer Methods in Applied Mechanics and Engineering 290 (2015) 1–18.

[33] N. Hagmeyer, M. Mayr, I. Steinbrecher, A. Popp, One-way coupled fluid-beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa, Adv. Model. Simul. Eng. Sci. 9 (1) (2022) 9. doi:10.1186/s40323-022-00222-y.

[34] N. Hagmeyer, M. Mayr, A. Popp, Fully coupled mortar-type embedding of one-dimensional fibers into three-dimensional fluid flow, CoRR abs/2306.02456 (2023). arXiv:2306.02456, doi:10.48550/arXiv.2306.02456.

[35] C. S. Peskin, The immersed boundary method, Acta numerica 11 (2002) 479–517.
[36] Y. Kim, C. S. Peskin, Penalty immersed boundary method for an elastic boundary with mass, Physics of Fluids 19 (5) (2007).

[37] D. Boffi, N. Cavallini, L. Gastaldi, Finite element approach to immersed boundary method with different fluid and solid densities, Mathematical Models and Methods in Applied Sciences 21 (12) (2011) 2523-2550.

[38] A. P. S. Bhalla, R. Bale, B. E. Griffith, N. A. Patankar, A unified mathematical framework and an adaptive numerical method for fluid-structure interaction with rigid, deforming, and elastic bodies, Journal of Computational Physics 250 (2013) 446-476.

[39] S. Tschisgale, J. Fröhlich, An immersed boundary method for the fluid-structure interaction of slender flexible structures in viscous fluid, Journal of Computational Physics 423 (2020) 109801.

[40] D. Boffi, N. Cavallini, L. Gastaldi, The finite element immersed boundary method with distributed lagrange multiplier, SIAM Journal on Numerical Analysis 53 (6) (2015) 2584-2604.

[41] D. Boffi, L. Gastaldi, A fictitious domain approach with lagrange multiplier for fluid-structure interactions, Numerische Mathematik 135 (2017) 711-732.

[42] R. Glowinski, T.-W. Pan, T. I. Hesla, D. D. Joseph, A distributed lagrange multiplier/fictitious domain method for particulate flows, International Journal of Multiphase Flow 25 (5) (1999) 755-794.

[43] E. L. Bouzarth, M. L. Minion, Modeling slender bodies with the method of regularized stokeslets, Journal of Computational Physics 230 (10) (2011) 3929-3947.

[44] T. T. Bringley, C. S. Peskin, Validation of a simple method for representing spheres and slender bodies in an immersed boundary method for stokes flow on an unbounded domain, Journal of Computational Physics 227 (11) (2008) 5397-5425.

[45] A.-K. Tornberg, M. J. Shelley, Simulating the dynamics and interactions of flexible fibers in stokes flows, Journal of Computational Physics 196 (1) (2004) 8-40.

[46] A.-K. Tornberg, K. Gustavsson, A numerical method for simulations of rigid fiber suspensions, Journal of Computational Physics 215 (1) (2006) 172-196.

[47] R. Cortez, L. Fauci, A. Medovikov, The method of regularized stokeslets in three dimensions: analysis, validation, and application to helical swimming, Physics of Fluids 17 (3) (2005).

[48] Y. Mori, L. Ohm, D. Spirn, Theoretical justification and error analysis for slender body theory, Communications on Pure and Applied Mathematics 73 (6) (2020) 1245-1314. doi:10.1002/cpa.21872.

[49] Y. Mori, L. Ohm, D. Spirn, Theoretical justification and error analysis for slender body theory with free ends, Archive for Rational Mechanics and Analysis 235 (3) (2020) 1905-1978. doi:10.1007/s00205-019-01458-6.

[50] U. Khristenko, S. Schuss, M. Kruger, F. Schmidt, B. Wohlmuth, C. Hesch, Multidimensional coupling: A variationally consistent approach to fiber-reinforced materials, Computer Methods in Applied Mechanics and Engineering 382 (2021) 113869.

[51] M. Kuchta, F. Laurino, K.-A. Mardal, P. Zunino, Analysis and approximation of mixed-dimensional pdes on 3d-1d domains coupled with lagrange multipliers, SIAM Journal on Numerical Analysis 59 (1) (2021) 558-582.

[52] M. Boulakia, C. Grandmont, F. Lespagnol, P. Zunino, Mathematical and numerical analysis of reduced order interface conditions and augmented finite elements for mixed dimensional problems, Available at SSRN 4846685 (2024).

[53] S. P. Timoshenko, On the correction for shear of the differential equation for transverse vibrations of prismatic bars., Phil Mag Ser 41 (1921) 744-764.
[54] S. P. Timoshenko, X. on the transverse vibrations of bars of uniform cross-section, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 43 (253) (1922) 125-131.

[55] K. Arunakirinathar, B. Reddy, Mixed finite element methods for elastic rods of arbitrary geometry, Numerische Mathematik 64 (1993) 13–43.

[56] J. Simo, L. Vu-Quoc, A three-dimensional finite-strain rod model. part ii: Computational aspects, Computer Methods in Applied Mechanics and Engineering 58 (1) (1986) 79–116. doi:https://doi.org/10.1016/0045-7825(86)90079-4.

[57] C. Meier, A. Popp, W. Wall, Geometrically exact finite element formulations for slender beams: Kirchhoff–love theory versus simo–reissner theory, Archives of Computational Methods in Engineering 26 (1) (2019) 163–243. doi:10.1007/s11831-017-9232-5.
URL https://www.scopus.com/inward/record.uri?eid=2-s2.0-85021786336&doi=10.1007%2fs11831-017-9232-5&partnerID=40&md5=ba8dee007e1539004f0b0b97c7eb6b4d

[58] E. Zupan, M. Saje, D. Zupan, On a virtual work consistent three-dimensional reissner-simo beam formulation using the quaternion algebra, Acta Mechanica 224 (2013) 1709–1729.

[59] V. Sonneville, A. Cardona, O. Brüls, Geometrically exact beam finite element formulated on the special euclidean group se (3), Computer Methods in Applied Mechanics and Engineering 268 (2014) 451–474.

[60] S. Eugster, C. Hesch, P. Betsch, C. Glocker, Director-based beam finite elements relying on the geometrically exact beam theory formulated in skew coordinates, International Journal for Numerical Methods in Engineering 97 (2) (2014) 111–129.

[61] F. Demoures, F. Gay-Balmaz, S. Leyendecker, S. Ober-Blöbaum, T. S. Ratiu, Y. Weinand, Discrete variational lie group formulation of geometrically exact beam dynamics, Numerische Mathematik 130 (2015) 73–123.

[62] S. Krenk, Non-linear modeling and analysis of solids and structures, Cambridge University Press, 2009.

[63] J. Dai, Euler-rodrigues formula variations, quaternion conjugation and intrinsic connections, Mechanics and Machine Theory 92 (2015) 144–152. doi:https://doi.org/10.1016/j.mechmachtheory.2015.03.004.

[64] M. Boulakia, C. Grandmont, F. Lespagnol, P. Zunino, Numerical approximation of the Poisson problem with small holes, using augmented finite elements and defective boundary conditions., working paper or preprint (Jan. 2023).
URL https://inria.hal.science/hal-03501521

[65] M. Fernández, J.-F. Gerbeau, Algorithms for fluid-structure interaction problems, in: Cardiovascular mathematics, Vol. 1 of MS&A. Model. Simul. Appl., Springer Italia, Milan, 2009, pp. 307–346.

[66] T. E. Tezduyar, Stabilized finite element formulations for incompressible flow computations, Advances in applied mechanics 28 (1991) 1–44.

[67] D. Chapelle, A locking-free approximation of curved rods by straight beam elements, Numerische Mathematik 77 (1997) 299–322.

[68] G. Barrenechea, C. González, A stabilized finite element method for a fictitious domain problem allowing small inclusions, Numer. Methods Partial Differential Equations 34 (1) (2018) 167–183.

[69] M. Fernández, M. Moubachir, A Newton method using exact Jacobians for solving fluid-structure coupling, Comp. & Struct. 83 (2005) 127–142.

[70] https://gitlab.inria.fr/felisce/felisce (2023).

[71] I. Babuška, The finite element method for elliptic equations with discontinuous coefficients, Computing 5 (3) (1970) 207–213.
[72] V. Girault, R. Glowinski, T. Pan, A fictitious domain method with distributed multiplier for the stokes problem, applied nonlinear analysis, Applied nonlinear analysis (1999) 159–174.

[73] Corti, D.C., Delay, G., Fernández, M.A., Vergnet, F., Vidrascu, M., Low-order fictitious domain method with enhanced mass conservation for an interface stokes problem, ESAIM: M2AN 58 (1) (2024) 303–333.

[74] K. Ohmori, N. Saito, Flux-free finite element method with Lagrange multipliers for two-fluid flows, J. Sci. Comput. 32 (2) (2007) 147–173.

[75] F. Chouly, P. Hild, Y. Renard, Finite Element Approximation of Contact and Friction in Elasticity, Vol. 48, Springer Nature, 2023.