MOX 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

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.

Introduction

Modeling and simulating 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 [START_REF] Fan | Simulation of 3d parachute fluid-structure interaction based on nonlinear finite element method and preconditioning finite volume method[END_REF][START_REF] Huang | Simulation of parachute inflation dynamics using an eulerian computational framework for fluid-structure interfaces evolving in high-speed turbulent flows[END_REF][START_REF] Kim | 3-d parachute simulation by the immersed boundary method[END_REF][START_REF] Sathe | Fluidstructure interaction modeling of complex parachute designs with the space-time finite element techniques[END_REF] to submerged structures or vegetation [START_REF] Peaceman | Interpretation of well-block pressures in numerical reservoir simulation (includes associated paper 6988)[END_REF][START_REF] Peaceman | Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability[END_REF][START_REF] Cerroni | Mathematical analysis, finite element approximation and numerical solvers for the interaction of 3d reservoirs with 1d wells[END_REF][START_REF] Wang | Fluid-structure interaction of flexible submerged vegetation stems and kinetic turbine blades[END_REF][START_REF] Luhar | Interaction between flow, transport and vegetation spatial structure[END_REF][START_REF] Mattis | An immersed structure approach for fluid-vegetation interaction[END_REF][START_REF] Chen | Eulerian-lagrangian flow-vegetation interaction model using immersed boundary method and openfoam[END_REF], numerical simulations provide insights which are unattainable by experimentation alone. Its interest also extends to biomedical applications, such as the design of vascular stents [START_REF] Tambača | Mathematical modeling of vascular stents[END_REF][START_REF] Tambača | A novel approach to modeling coronary stents using a slender curved rod model: a comparison between fractured xience-like and palmaz-like stents[END_REF][START_REF] Zunino | Integrated stent models based on dimension reduction: review and future perspectives[END_REF] and the simulation of cilia, flagella [START_REF] Decoene | A continuum active structure model for the interaction of cilia with a viscous fluid[END_REF] and micro-swimmers [START_REF] Alouges | Self-propulsion of slender micro-swimmers by curvature control: N-link swimmers[END_REF][START_REF] Giraldi | Controllability and optimal strokes for n-link microswimmer[END_REF][START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF]. 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., [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF][START_REF] Reissner | On one-dimensional finite-strain beam theory: the plane problem[END_REF][START_REF] Meier | An objective 3d large deformation finite element formulation for geometrically exact curved kirchhoff rods[END_REF][START_REF] Meier | A locking-free finite element formulation and reduced models for geometrically exact kirchhoff rods[END_REF][START_REF] Meier | Geometrically exact beam elements and smooth contact schemes for the modeling of fiber-based materials and structures[END_REF]). 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., [START_REF] Baaijens | A fictitious domain/mortar element method for fluid-structure interaction[END_REF][START_REF] Boilevin-Kayl | Numerical methods for immersed FSI with thinwalled structures[END_REF]). 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., [START_REF] Chapelle | Modeling of the inclusion of a reinforcing sheet within a 3D medium[END_REF][START_REF] Landajuela | Coupling schemes for the FSI forward predication challenge: comparative study and validation[END_REF]). 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., [START_REF] Maniyeri | Numerical study on the rotation of an elastic rod in a viscous fluid using an immersed boundary method[END_REF][START_REF] Huang | An embedded boundary approach for resolving the contribution of cable subsystems to fully coupled fluid-structure interaction[END_REF][START_REF] Ausas | A finite element method for simulating soft active non-shearable rods immersed in generalized newtonian fluids[END_REF]), while the second one reduces the complexity of the problem by considering the interface as the (1D) centerline of the beam (see, e.g., [START_REF] Griffith | Simulating an elastic ring with bend and twist by an adaptive generalized immersed boundary method[END_REF][START_REF] Wiens | Simulating flexible fiber suspensions using a scalable immersed boundary algorithm[END_REF][START_REF] Wang | Fluid-structure interaction of flexible submerged vegetation stems and kinetic turbine blades[END_REF][START_REF] Tschisgale | An immersed boundary method for the fluid-structure interaction of slender flexible structures in viscous fluid[END_REF][START_REF] Hagmeyer | One-way coupled fluid-beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa[END_REF][START_REF] Hagmeyer | Fully coupled mortar-type embedding of one-dimensional fibers into three-dimensional fluid flow[END_REF]).

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 the solid (see, e.g., [START_REF] Boilevin-Kayl | Numerical methods for immersed FSI with thinwalled structures[END_REF]). 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. Noteworthy among these methodologies is the Immersed Boundary method (IB) (see, e.g., [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Kim | Penalty immersed boundary method for an elastic boundary with mass[END_REF][START_REF] Boffi | Finite element approach to immersed boundary method with different fluid and solid densities[END_REF]), 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, as evidenced by works such as [START_REF] Griffith | Simulating an elastic ring with bend and twist by an adaptive generalized immersed boundary method[END_REF][START_REF] Wiens | Simulating flexible fiber suspensions using a scalable immersed boundary algorithm[END_REF][START_REF] Bhalla | A unified mathematical framework and an adaptive numerical method for fluid-structure interaction with rigid, deforming, and elastic bodies[END_REF][START_REF] Wang | Fluid-structure interaction of flexible submerged vegetation stems and kinetic turbine blades[END_REF][START_REF] Tschisgale | An immersed boundary method for the fluid-structure interaction of slender flexible structures in viscous fluid[END_REF]. Another category comprises methods that use Lagrange multipliers in lieu of Dirac functions, referred to as Fictitious Domain Methods (FD) (see, e.g., [START_REF] Baaijens | A fictitious domain/mortar element method for fluid-structure interaction[END_REF][START_REF] Boffi | The finite element immersed boundary method with distributed lagrange multiplier[END_REF][START_REF] Boffi | A fictitious domain approach with lagrange multiplier for fluid-structure interactions[END_REF][START_REF] Glowinski | A distributed lagrange multiplier/fictitious domain method for particulate flows[END_REF]).

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. [START_REF] Bouzarth | Modeling slender bodies with the method of regularized stokeslets[END_REF][START_REF] Bringley | Validation of a simple method for representing spheres and slender bodies in an immersed boundary method for stokes flow on an unbounded domain[END_REF][START_REF] Tornberg | Simulating the dynamics and interactions of flexible fibers in stokes flows[END_REF][START_REF] Tornberg | A numerical method for simulations of rigid fiber suspensions[END_REF][START_REF] Cortez | The method of regularized stokeslets in three dimensions: analysis, validation, and application to helical swimming[END_REF][START_REF] Mori | Theoretical justification and error analysis for slender body theory[END_REF][START_REF] Mori | Theoretical justification and error analysis for slender body theory with free ends[END_REF]). 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 coupling conditions on the beam centerline. From a numerical perspective, given that the subspace used for the approximation of the continuous solution space consists mainly of continuous functions, various techniques can be used to compute solutions for the discrete problem, such as penalty methods [START_REF] Hagmeyer | One-way coupled fluid-beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa[END_REF] or mortar-type methods [START_REF] Hagmeyer | Fully coupled mortar-type embedding of one-dimensional fibers into three-dimensional fluid flow[END_REF]. One advantage of this method is that it allows coarse fluid mesh to be used in the vicinity of the beam. However, on a continuous level, two questions naturally emerge. The first one is if the continuous problem associated with this solution is well-posed. Indeed, imposing 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 is the accuracy of the model, whether imposing coupling conditions on the centerline is a good approximation or not 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 [START_REF] Chapelle | Modeling of the inclusion of a reinforcing sheet within a 3D medium[END_REF] 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 a immersed slender solid, modeled as a 1D structure. 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 [START_REF] Khristenko | Multidimensional coupling: A variationally consistent approach to fiber-reinforced materials[END_REF] for a similar idea in the context of reinforced materials). More precisely, on each cross section of the beam, the original kinematic constraint is replaced by a constraint which acts only on the first N modes of the Fourier expansion of the velocity at the interface. To some extent, this approach serves could be seen as a bridge, rectifying certain disparities among the aforementioned sets of fluid-beam interaction methodologies, whether employing a 2D or a 1D representation of the interface and aims to be a compromise between accuracy, computational cost and robustness.

This method has been used in [START_REF] Kuchta | Analysis and approximation of mixed-dimensional pdes on 3d-1d domains coupled with lagrange multipliers[END_REF] for N " 0 for stationary elliptic problems. Notably, 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 similar framework to [START_REF] Kuchta | Analysis and approximation of mixed-dimensional pdes on 3d-1d domains coupled with lagrange multipliers[END_REF] but for an arbitrary number of Fourier modes, the authors of [START_REF] Boulakia | Numerical approximation of the poisson problem with small holes, using augmented finite elements and defective boundary conditions[END_REF] 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 ě 1. In the context of fluid-structure interaction this point is crucial to correctly capture the fluid-structure dynamic.

In order to illustrate the accuracy of the proposed strategy, we compare it with two others 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.

The rest of the paper is organized as follows. In Section 2, we describe the fluid and solid models, as well as the coupling conditions for the fully resolved ALE coupled model and two reduced mixeddimensional coupled models. The discrete formulation of the proposed model is provided in Section 3. A few computer implementation details of the considered solution procedure are also given. Section 4 is devoted to the numerical experiments, where we investigate how different model parameters impact the accuracy and robustness of the proposed approach. Finally, a summary of the results of this work together with some research perspectives are drawn in Section 5.

Continuous setting and mathematical formulation

In this section, we briefly present the beam and fluid models considered in this paper. We choose a linear Timoshenko beam formulation for the solid and the fluid is assumed to satisfied the incompressible Navier-Stokes equations. Our approach will be presented in Section 2.2.5 for the fluid and in Section 2.3.4 for the coupling conditions.

Timoshenko beam model

Let pe x , e y , e z q 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 p ω ε def " εE ˆp0, Lq. Here, E stands for an elliptical region in the pxOyq plane (major axis of length a and a minor axis of length 1) and ε is a parameter related to the thickness of the beam. For s P p0, Lq, we denote by p ω ε psq def " εE ˆtsu the cross-section centered in se z , such that p ω ε " Ť sPp0,Lq p ω ε psq (see Figure 1). We assume that p corresponds to the position of a beam at the initial time. The motion of the beam is given in terms of the map ϕ : p ω ε ˆR`Ñ R 3 , so that the current configuration of the beam ω ε ptq is given by ω ε ptq " ϕpp ω ε , tq (see Figure 2). In what follows, we consider a Timoshenko beam model, so that 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: ϕpp x, tq " se z `rps, tq `Λps, tqpp x ´se z q, s " p x ¨ez , pp x, tq P p

ω ε 0 L se z p ω ε psq Bp ω ε e y e z e x p ω ε " Ť sPp0,Lq
ω ε ˆR`, (2.1)
where r : p0, Lq ˆR`Ñ R 3 stands for the displacement of the centerline and Λ : p0, Lq ˆR`Ñ SOp3q describes the rotation of the cross-sections. Here, SOp3q stands for the Special Orthogonal group given by SOp3q

def

" tΛ P R 3ˆ3 | Λ T Λ " I 3 , detΛ " 1u, (2.2)
where I 3 denotes the identity matrix in R 3ˆ3 . Note, in particular, that the centerline Bptq and the cross-sections pω ε ps, tqq sPp0,Lq of the current configuration of the beam are given respectively by Due to the orthogonality of the rotation matrix Λ, it can be given only in terms of three independent parameters. One option is to consider the rotational vector, defined by

Bptq def " ϕ `t0u ˆt0u ˆp0, Lq , t ˘, ω ε ps, tq def " ϕpp ω ε psq, tq.
θ def " θn,
where n is a unit vector defining the axis of rotation and θ " ? θ ¨θ is the angle of rotation. We set Θ the skew symetric matrix associated to the vector θ such that for all v P R 3 , Θv " θ ^v, where the symbol ^denotes the cross product in three dimensions. Using Rodrigues' rotation formula (see, e.g., [START_REF] Dai | Euler-rodrigues formula variations, quaternion conjugation and intrinsic connections[END_REF]), the relation between the rotation matrix and the rotational vector is then given by

Λ " I 3 `sin θ θ Θ `1 ´cos θ θ 2 Θ 2 . (2.
3)

The dynamics of the beam are hence given in terms of r and θ, which will here be described by a variant of the linear elastic model introduced in [START_REF] Arunakirinathar | Mixed finite element methods for elastic rods of arbitrary geometry[END_REF] for the static case. The model accounts for bending, shear, torsion, and membrane effects within a small displacements and rotations framework. More precisley, the dynamics of the beam are thus described by the following balance equations:

# ρ b AB 2 t r ´GB s pB s r ´θ ^ez q " f in p0, Lq ˆR`, ρ b IB 2 t θ ´EB 2 s θ ´ez ^G pB s r ´θ ^ez q " m in p0, Lq ˆR`, (2.4)
where the matrices A, I, E and G are given by

A def " » - A 0 0 0 A 0 0 0 A fi fl , I def " » - I x 0 0 0 I y 0 0 0 J fi fl , E def " » - EI x 0 0 0 EI y 0 0 0 GJ fi fl , G def " » - GAκ 0 0 0 GAκ 0 0 0 EA fi fl
and the symbols ρ b , E, G, I x , I y , J, A and κ, 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.

Navier-Stokes fluid model

In this paper, the fluid is modeled by the incompressible Navier-Stokes equations, written with full or reduced formulations that will lead to different numerical schemes. For the sake of simplicity we present the different models and strategies in the case where homogeneous Dirichlet boundary conditions are prescribed on the fluid boundary that is not the fluid-structure interface. We assume here that r and Λ are regular enough and that for all t, ϕ t def " ϕp¨, tq : p ω ε Ñ ω ε ptq is one to one. The beam is supposed to move within a fixed and bounded domain Ω Ă R 3 and we assume that it is surrounded by a fluid which occupies at each time t the time-dependent domain

Ω f ε ptq def " ΩzĎ ω ε ptq Ă R 3 .
In what follows, Σ ε ptq def " BΩ f ε ptq X Bω ε ptq denotes the current configuration of the fluid-solid interface. The fluid domain boundary is partitionned as BΩ f ε ptq " Σ ε ptq Y Γ, with Γ Ă BΩ.

Eulerian formalism

The fluid equations in the Eulerian formalism read: find the fluid velocity u : O Ñ R 3 and fluid pressure p

: O Ñ R such that $ ' ' ' ' & ' ' ' ' % ρ f B t u `ρf u ¨∇u ´∇ ¨σpu, pq " 0 in Ω f ε ptq, ∇ ¨u " 0 in Ω f ε ptq, u " 0 on Γ, u " u b on Σ ε ptq. (2.5)
Here, ρ f stands for the fluid density, µ for the dynamic viscosity of the fluid, and σpu, pq denotes the Cauchy stress tensor given by σpu, pq

def " 2µεpuq ´pI, εpuq def " 1 2 `∇u `p∇uq T ˘.
Lastly, u b denotes the Eulerian velocity of the beam, which will be linked to the beam unknowns in Section 2.3 below.

ALE formalism

In the ALE formalism, the fluid domain Ω f ε ptq is parametrized by a one-to-one mapping A : p " Ωzp ω ε denotes a given reference fluid domain configuration. We also set A t p¨q " Ap¨, tq. We can hence defined the ALE time derivative

Ω f ε ˆR`Ñ R 3 such that Ω f ε ptq " Ap p Ω f ε ,
B t | A upx, tq " d dt `upApp x, tq, tq ˘, p x " A ´1 t pxq.
By a simple use of the chain rule, we obtain the following relation

B t u " B t | A u ´p p w f ˝A´1 t q ¨∇u, (2.6)
where p w f def " B t A " B t p d f is the fluid domain velocity. In what follows, for any function p v defined in p Ω f ε , we will denoted by v " p v ˝pA t q ´1 its Eulerian counterpart. The fluid equations (2.5) admit the following equivalent formulation: find the fluid velocity u : O Ñ R 3 and fluid pressure p

: O Ñ R such that $ ' ' ' ' & ' ' ' ' % ρ f B t | A u `ρf pu ´wf q ¨∇u ´∇ ¨σpu, pq " 0 in Ω f ε ptq, ∇ ¨u " 0 in Ω f ε ptq, u " 0 on BΩ, u " u b on Σ ε ptq. (2.7)
The main advantage of (2.7) with respect of (2.6) is that it facilitates the discretization in time.

Fictitious domain approach

In what follows, for any Lipschitz domain D in R 3 , we will make use of the standard Sobolev spaces H s pDq and L m pDq ps P R, m P N X r1, `8rq. We also denote by H 1 S pDq the closed subspace of the functions in H 1 pDq with zero trace on S, and by L 2 0 pDq the subspace of the functions in L 2 pDq of functions with zero mean in D. The scalar product in L 2 pDq is denoted by p¨, ¨qD . The bold notation will be used for vector Sobolev spaces. In the fictitious domain approach, the fluid equations (2.5) are solved in the whole domain Ω while enforcing the kinematic constraint (2.5) 4 on Σ ε ptq through Lagrange multipliers: for t P R `, xµ, uptqy Σεptq " xµ, u b ptqy Σεptq , @ µ P H ´1 2 pΣ ε ptqq, (

where x¨, ¨yD denotes the duality pairing between H ´1 2 pDq and H 1 2 pDq. One of the advantages of the fictitious domain methods is that it allows the use of arbitrary fluid meshes constant in time; for further details on this approach, refer to [START_REF] Girault | Error analysis of a fictitious domain method applied to a dirichlet problem[END_REF][START_REF] Girault | A fictitious domain method with distributed multiplier for the stokes problem, applied nonlinear analysis[END_REF]. Let V def " H 1 BΩ pΩq and Q def " L 2 0 pΩq be the functional spaces for the fluid velocity and pressure, respectively, the considered fictitious domain formulation of the fluid problem (2.5) reads: we look for pu, p, λq such that puptq, pptq, λptqq P V ˆQ ˆH´1

Fictitious domain approach with kinematic constraint on the beam centerline

We present here a first reduced model based on a fictitious domain approach. This formulation is inspired by the reduced models reported in [START_REF] Hagmeyer | One-way coupled fluid-beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa[END_REF][START_REF] Hagmeyer | Fully coupled mortar-type embedding of one-dimensional fibers into three-dimensional fluid flow[END_REF] (see also [START_REF] Griffith | Simulating an elastic ring with bend and twist by an adaptive generalized immersed boundary method[END_REF][START_REF] Wiens | Simulating flexible fiber suspensions using a scalable immersed boundary algorithm[END_REF][START_REF] Wang | Fluid-structure interaction of flexible submerged vegetation stems and kinetic turbine blades[END_REF]) and assumes that, owing to the small thickness of the beam ε, the kinematic constraint can be simply enforced on the beam centerline. The Dirichlet boundary conditions u b are then enforced with Lagrange multipliers defined on the beam centerline as follows: for t P R `, pµ, uptqq Bptq " `µ, u b ptq ˘Bptq @ µ P L 2 pBptqq.

(2.10)

This approach reduces the size of the problem by passing for the boundary conditions from a twodimensional constraint on the beam interface to a one-dimensional constraint on the beam centerline.

The weak form of the fluid equations then writes: we look for pu, p, λq such that puptq, pptq, λptqq P V ˆQ ˆL2 pBptqq a.e. t P R `and ρ f pB t uptq, vq Ω `af Ω `uptq; puptq, pptqq, pv, qq ˘´pλptq, vq Bptq ``µ, uptq ´ub ptq ˘Bptq " 0 (2.11) for all pv, q, µq P V ˆQ ˆL2 pBptqq.

It should be noted that the relation (2.10) is only well-defined for u P H 1`η pΩq, with η ą 0. Nevertheless, in the context of problem (2.11), this minimum regularity for u might not be attained (due to the singularity induced by the term pλptq, vq Bptq ). This prevents us from establishing both the well-posedness of the continuous problem (2.11) and the accuracy of the associated numerical methods, which is a major drawback.

Fictitious domain approach with reduced order kinematic condition

In this section, we introduce our method, which is also based on a fictitious domain approach with Lagrange multipliers, following a similar approach to (2.9) and (2.11). Our objective, similar to (2.11), is to take advantage of the slenderness of the beam to reduce the size of the problem. To do so, we replace the Lagrange multiplier space H ´1 2 pΣ ε ptqqq by a reduced-order space denoted as F N ptq. In order to define the space F N ptq, we first introduce a local coordinate ν on each cross-section of the reference configuration so that any point p

x on Bp ω ε can be parametrized as p

x " se z `ε`c ospνqe x `a sinpνqe y ˘.

(2.12)

In particular, any function v defined on Bp ω ε can be expressed as follows:

vps, νq def " vpse z `ε pcospνqe x `a sinpνqe y qq @ps, νq P p0, Lq ˆp0, 2πq .

Let F N be a finite dimensional subspace of L 2 p0, 2πq, where the dimension varies with N , we set for all t P R `, the space F N ptq defined by

F N ptq def " # v ˝ϕ´1 t ˇˇˇv ps, νq " N ÿ k"0 α k psqg k pνq, α k P L 2 p0, Lq, g k P F N + . (2.13)
The idea of (2.13) is to take advantage of the tensorization of Bp ω ε by p0, Lq ˆp0, 2πq to discretize the space L 2 p0, 2πq using the finite dimensional space F N . This reduced representation is then transported to the current configuration of the beam interface Bω ε ptq via composition with ϕ ´1 t . The expression reduced order comes from the fact that any function v P F N ptq is now uniquely determined by the set of functions tα k u 0ďkďN Ă L 2 p0, Lq, thus passing from a functional space H ´1 2 pBω ε ptqq originally defined in a two-dimensional domain to a set of N `1 functions defined on a 1D domain. In particular, the reduced-order interface conditions on Bω ε ptq read as follows, for t P R for all pv, q, µ N q P V ˆQ ˆFN ptq.

Remark 2.1. The choice of the scalar product p¨, ¨qε,t is justified by the error approximation between the solution of the reduced problem and the full problem, where we consider H ´1 2 pBω ε ptqq as the Lagrange multipliers space. Specifically, in a paper currently in preparation, we give numerical evidence that, for a 2D toy problem and F N chosen as in (2.22), the order of convergence with respect to ε is improved compared to the scalar product p¨, ¨qBωεptq .

Coupled problems: interface coupling conditions

In this section, we provide the interface coupling conditions between the beam model (2.4) and the different fluid modeling options given by (2.7), (2.11) or (2.15). The numerical approximation of the resulting coupled problems is addressed in Section 3.

Evaluation of the beam velocity

In this section, we provide an explicit expression velocity of the beam u b in its current configuration ω ε ptq. First, we consider the Lagrangian velocity p u b def " B t ϕ, so that from (2.1) we get p u b pp x, tq " B t rps, tq `Bt Λps, tqpp x ´se z q, s " p x ¨ez (2.16) for all p x P p ω ε . On the other hand, from (2.1) and the orthogonality of Λ, we have p x ´se z " Λ T ps, tqpx ´rps, tq ´se z q, for all x P ω ε ptq, with p x " ϕ ´1 t pxq and s " p x ¨ez . By inserting this expression into (2.16) we obtain that the Eulerian velocity u b of the beam is given by u b px, tq def " p u b pϕ ´1 t pxq, tq " B t rps, tq `Bt Λps, tqΛ T ps, tq `x ´rps, tq ´se z ˘.

(2.17)

From the orthogonality of Λ, we infert that the matrix B t ΛΛ T is skew symmetric, so that there exists a vector field w in p0, Lq ˆR`, termed angular velocity vector, such that

B t ΛΛ T v " w ^v @v P R 3 .
Hence the relation (2.17) can be equivalently re-written as u b px, tq " B t rps, tq `wps, tq ^px ´rps, tq ´se z q, p

x " ϕ ´1 t pxq, s " p x ¨ez (2.18) for all x P ω ε ptq. The velocity of the beam corresponds to a rigid-body velocity per cross-section.

In this paper we shall make extensive use of an extension operator from the beam centerline, t0u t0u ˆp0, Lq, to the bean reference domain p ω ε , defined as follows: for any vector field v : p0, Lq ˆR`Ñ R 3 we define the vector field s v : p ω ε ˆR`Ñ R 3 by the relation s vpp x, tq " vps, tq, s " p x ¨ez , @ pp x, tq P p ω ε ˆ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. Using (2.18), we can then express the Eulerian velocity of the beam u b as:

u b " Ď B t r ˝ϕ´1 t `s w ˝ϕ´1 t ^`I 3 ´s r ˝ϕ´1 t ´pϕ ´1 t ¨ez qe z ˘. (2.19)
Remark 2.2. Note that the focus of this work is on the formulation of the coupling equations between the beam and the fluid. For simplicity, we have introduced a linear constitutive law defined within the reference configuration (Lagrangian form). However, the choice of beam model does not have an impact on the formulation of the coupling equations. Indeed, a similar formulation of the coupling equations can be obtained for more complex beam models with non-linear constitutive laws and valid for large deformations. Such models can be derived, for example, from the exact geometry beam theory introduced by Simo [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF] and Reissner [START_REF] Reissner | On one-dimensional finite-strain beam theory: the plane problem[END_REF]. In the numerical results presented in Section 4.3, a corotating beam formulation with a linear constitutive law, as described in [START_REF] Krenk | Non-linear modeling and analysis of solids and structures[END_REF], has been considered to address large beam deformations.

ALE formalism for the fluid

In this section, we introduce the coupling conditions associated to the ALE formulation introduced in Section 2.2.2. In what follows, the symbol L denotes a given lifting operator from Bp ω ε to p Ω f ε which vanishes on BΩ. Note that the extension provided by L is arbitrary inside p Ω f ε . The interface coupling conditions between the solid sub-problem (2.4) and the fluid sub-problem (2.7) are given by

$ ' ' ' ' ' & ' ' ' ' ' % ϕ " pI 3 ¨ez qe z `s r `s Λ `I3 ´pI 3 ¨ez qe z ˘in p ω ε , A " Lpϕq, p w f " B t A, Ω f ε ptq " Ap p Ω f ε , tq, u " u b " B t ϕ ˝A´1 t on Σ ε ptq, `f , δr ˘p0,Lq " ´`σpu, pqn, s δr ˝A´1 t ˘Σεptq , `m, δθ ˘p0,Lq " ´`σpu, pqn, Ď δθ ˝A´1 t ^pI 3 ´s r ˝A´1 t q ˘Σεptq , (2.20
) for all δr, δθ P Y, with Y the space of admissible displacements and rotations of the beam (typically H 1 p0, Lq or a subspace of it), and n stands for the exterior unit normal to BΩ f ε ptq. The coupling conditions of (2.20) respectively express the geometric compatibility, kinematic coupling and dynamic balance between the fluid and the beam at the interface.

Let F def " ∇A t and J def " det F denote, respectively, the gradient and Jacobian of A t . Using the Piola transform, the last two equalities of (2.20) reduce to

f psq " ´żBp ωεpsq Jσpu, pq ˝A F ´T p n, mpsq " ´żBp ωεpsq JpA ´s rq ^σpu, pq ˝A F ´T p n,
for any s P p0, Lq and where p n denotes the exterior unit normal to B p Ω f ε . In other words, for a given s P p0, Lq, the forcing terms f psq and mpsq, acting on the beam system (2.4), are given as the crosssection resultant of the fluid force and torque in the reference configuration.

Fictitious domain approach with coupling conditions on the beam centerline

In this approach, the fluid-structure interface is identified with the beam centerline and thickness effects are neglected in the interface coupling (see, e.g., [START_REF] Hagmeyer | One-way coupled fluid-beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa[END_REF][START_REF] Hagmeyer | Fully coupled mortar-type embedding of one-dimensional fibers into three-dimensional fluid flow[END_REF]). In particular, no rotational velocity is defined on the beam's centerline, so that the fluid torque acting on the beam is null. As a result, the interface coupling conditions between the solid sub-problem (2.4) and the fluid sub-problem (2.11) are given by:

$ ' ' ' ' & ' ' ' ' % ϕ " pI 3 ¨ez qe z `s r on t0u ˆt0u ˆp0, Lq , Bptq " ϕ `t0u ˆt0u ˆp0, Lq , t ˘, u b " B t r ˝ϕ´1 t on Bptq,
pf , δrq p0,Lq " ´`λ, δr ˝ϕ´1 t ˘Bptq @ δr P Y, m " 0.

(2.21) Note that the last relation completely neglects the rotational dynamics. In the next section, we propose an alternative reduced coupled model which overcomes this issue, as illustrated in the numerical experiments of Section 3.

Fictitious domain approach with reduced order interface coupling conditions

In this section, we introduce the coupling conditions associated to our method, which corresponds to the fictitious domain formulation introduced in Section 2.2.5. As finite dimensional space F N Ă L 2 p0, 2πq in (2.13), we consider the following vector counterpart of the Fourier space reported for a Poisson problem in [START_REF] Boulakia | Numerical approximation of the poisson problem with small holes, using augmented finite elements and defective boundary conditions[END_REF] (see also [START_REF] Khristenko | Multidimensional coupling: A variationally consistent approach to fiber-reinforced materials[END_REF]):

F N " span ␣ sinpkνqe j , cospkνqe j ( 0ďkďN,jPtx,y,zu . We hence get the following reduced Lagrange multiplier space

F N ptq " # v ˝ϕ´1 t ˇˇˇv ps, νq " a 0 psq `N ÿ k"1 pa k psq cospkνq `bk psq sinpkνqq , a k , b k P L 2 p0, Lq + (2.22)
for t P R `. The interface coupling conditions between the solid sub-problem (2.4) and the fluid subproblem (2.15) read as

$ ' ' ' ' ' & ' ' ' ' ' %
ϕ " pI 3 ¨ez qe z `s r `s Λ `I3 ´pI 3 ¨ez q e z ˘in p ω ε , Σ ε ptq " ϕp p Σ ε , tq,

u b " B t ϕ ˝ϕ´1 t on Σ ε ptq,
pf , δrq p0,Lq " ´`λ N , s δr ˝ϕ´1 t ˘ε,t @ δr P Y, pm, δθq p0,Lq " ´`λ N , Ď δθ ˝ϕ´1 t ^pI 3 ´s r ˝ϕ´1 t q ˘ε,t @ δθ P Y.

(2.23)

Remark 2.3. It should be noted that, given that the fluid computational domain is the whole domain Ω, the relation λ N " σpu N , p N qn does not hold. Instead, λ N corresponds to the jump of the fluid stress across Σ ε ptq, viz., λ N " ´`σpu, pq `n`´σ pu, pq ´n´˘o n Bω ε ptq.

Here, the indices `and ´denote trace values on Bω ε ptq taken from Ω f ε ptq and ω ε ptq, respectively. Nevertheless, numerical evidence provided in Section 3, suggest that the contribution of the internal fluid stress is negligible when the size ε of the cross-section is small. Energy estimate (small rotational velocity). The purpose of the last part of this section is to provide an energy estimate for the coupled reduced model (2.4), (2.15), (2.23) by assuming that the beam undergoes small rotations. 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 (2.23) 2 in which the Eulerian velocity of the beam, u b given by (2.19), is taken with a linearized angular velocity, namely,

w " B t θ, (2.24)
which yields

u b " Ď B t r ˝ϕ´1 t `Ď B t θ ˝ϕ´1 t ^`I 3 ´s r ˝ϕ´1 t ˘. (2.25)
Note that in this case the relation u b " B t ϕ ˝ϕ´1 t does not hold anymore. The relations (2.24) and (2.25) are only valid in a small rotational velocity framework.

The resulting problem in weak form read as follows: we look for pu, p, λ N q such that puptq, pptq, λ N ptq, rptq, θptqq P V ˆQ ˆFN ptq ˆY ˆY a.e. t P R `and

$ ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' % ϕ " s r `s Λ `I3 ´pI 3 ¨ez q e z ˘in p ω ε , Σ ε ptq " ϕp p Σ ε ,
where the mechanical energy of the system Eptq is defined by Proof. By testing (2.26) with pv, q, δr, δθ, µ N q " puptq, pptq, B t rptq, B t θptq, λ N ptqq we get ). In order to guarantee energy stability for this model, one simple option is to to update the fluid domain from a linearized version of (2.1) with Λ « I 3 `Θ. Another option would be to directly consider a non-linear model for the balance of the beam, such as the corotational formulation (see, e.g., [START_REF] Krenk | Non-linear modeling and analysis of solids and structures[END_REF]).

Eptq def " ρ f 2 }uptq} 2 0,Ω `ρb 2 }B t rptq} 2 A,p0,Lq `ρb 2 }B t θptq} 2 I,p0,Lq `1 2 › › prptq, θptqq
ρ f pB t u, u N q Ω `ρb pAB

Numerical method

In this section, we introduce a numerical scheme for the approximation of the coupled problem (2.4), (2.15), (2.23). In what follows, we denote by τ ą 0 the time-step length and, for n P N, we set t n def " nτ . For any time-dependent function v, v n stands for an approximation of vpt n q. We also introduce the standard notations

B τ v n def " 1 τ `vn ´vn´1 ˘, v n´1 2 def " 1 2 `vn´1 `vn ˘,
for the first-order backward difference and the mid-point value, respectively. Let D be an open polygonal convex domain of R 3 and T h pDq a triangulation of D with characteristic size h. We consider the following standard continuous Lagrange finite element spaces of degree k:

X k h pDq def " ␣ v P H 1 pDq | v| τ P rP k pKqs 3 , @ K P T h pDq ( , Q k h pDq def " ␣ v P H 1 pDq | v| τ P P k pKq, @ K P T h pDq ( .
We then introduce the finite element spaces

V h def " X 1 h pΩqXV and Q h def " Q 1 h
pΩqXQ for the approximation of the fluid velocity u and pressure p, respectively, while for the solid displacement r and rotation θ their approximations will be looked for into

Y H def " X 1 H p0, Lq X Y.
Here, h and H respectively denote the characteristic sizes of the fluid and solid (centerline) meshes.

Discrete formulation of the coupled model with reduced order interface conditions

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

ϕ n H " pI 3 ¨ez qe z `Ę r n´1 H `Ğ Λ n´1 H `I3 ´pI 3 ¨ez q e z ˘in p ω ε , Σ n ε " ϕ n H p p Σ ε q.
We can hence introduce the fully discrete Lagrange multiplier space

F n N,H , discrete counterpart of F N pt n q, as F n N,H def " # v h ˝pϕ n H q ´1 ˇˇˇv h ps, νq " a 0,H psq `N ÿ k"1 `ak,H psq cospkνq `bk,H psq sinpkνq ˘, a k,H , b k,H P Y H + , (3.1)
which simply amounts to replace L 2 p0, Lq by Y H in (2.22). Note that here a unique approximation space Y 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 sub-problem (2.15) and a mid-point scheme for the solid sub-problem (2.4). The remaining kinematic and dynamic conditions in (2.23) are discretized with an implicit scheme. The resulting solution procedure is of strongly coupled nature (see, e.g., [START_REF] Fernández | Algorithms for fluid-structure interaction problems[END_REF]). By gathering all the above mentioned ingredients, the proposed numerical approximation of (2.4), (2.15), (2.23) is detailed in Algorithm 1. In (3.3), the discrete fluid tri-linear form is given by

a f Ω,h `zh ; pu h , p h q, pv h , q h q ˘def " c n h pz h ; u h , v h q `aΩ `pu h , p h q, pv h , q h q ˘`s Ω,h `zh ; pu h , p h q, pv h , q h q ˘, c n h pz h ; u h , v h q def " c Ω pz h ; u h , v h q `ρf 2 `p∇ ¨zh q u h , v h ˘Ω, s Ω,h `zh ; pu h , p h q, pv h , q h q ˘def " ÿ KPT h pΩq ż K δ h `ρf pz h ¨∇qu h `∇p h ˘¨`ρ f pz h ¨∇qv h `∇q h ˘, δ h def " γ S ˜ρf d 4 τ 2 `16µ 2 h 4 pρ f q 2 `4|z h | 2 h 2 ¸´1 ,
where s Ω,h corresponds to the SUPG/PSPG stabilization (see, e.g., [START_REF] Tezduyar | Stabilized finite element formulations for incompressible flow computations[END_REF]) and γ S ą 0 is a user-defined parameter. For the discrete beam bi-linear form, we consider the locking-free formulation introduced in [ For n ě 1,

Step 1: Update interface position:

ϕ n H " pI 3 ¨ez qe z `Ę r n´1 H `Ğ Λ n´1 H `I3 ´pI 3 ¨ez q e z ˘in p ω ε , Σ n ε " ϕ n H p p Σ ε q, (3.2)
where Λ n´1 H is obtained from θ n´1 H using non linear formula 2.3.

Step 2:

Find pu n h , p n h , λ n N,H , r n H , θ n H q P V h ˆQh ˆFn N,H ˆYH ˆYH with 9 r n´1 2 H " B τ r n H and 9 θ n´1 2 H " B τ θ n H such that $ ' ' ' ' ' ' & ' ' ' ' ' ' % ρ f pB τ u n h , v h q Ω `af Ω,h `un´1 h ; pu n h , p n h q, pv h , q h q ρb pAB τ 9 r n H , δr H q p0,Lq `ρb `IB τ 9 θ n H , δθ H ˘p0,Lq `ab H `pr n´1 2 H , θ n´1 2 H q, pδr H , δθ H q λn N,H , v h ´Ě δr H ˝pϕ n H q ´1 ´Ě δθ H ˝pϕ n H q ´1 ^´I 3 ´Ę r n´1 H ˝pϕ n H q ´1¯¯ε ,n ``µ N,H , u h ´Bτ ϕ n`1 H ˝pϕ n H q ´1˘ε ,n " 0 (3.3)
for all pv h , q h , µ N,H , δr H , θ H q P V h ˆQh ˆFn N,H ˆYH ˆYH .

Energy estimate (small rotational velocity). In the spirit of Section 2.3.4, we provide an energy estimate for Algorithm 1 under a small rotational velocity framework of (2.26). In this context, Algorithm 1 reduces to the following discrete coupled problem (3.3) now writes: find

pu n h , p n h , λ n N,H , r n H , θ n H q P V h ˆQh ˆFn N,H ˆYH ˆYH with 9 r n´1 2 H " B τ r n H and 9 θ n´1 2 H " B τ θ n H such that $ ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' % ϕ n H " pI 3 ¨ez qe z `Ę r n´1 H `Ğ Λ n´1 H `I3 ´pI 3 ¨ez q e z ˘in p ω ε , Σ n ε " ϕ n H p p Σ ε q, ρ f pB τ u n h , v h q Ω `af Ω,h `un´1 h ; pu n h , p n h q, pv h , q h q ρb `AB τ 9 r n H , δr H ˘p0,Lq `ρb `IB τ 9 θ n H , δθ H ˘p0,Lq `ab H `pr n´1 2 H , θ n´1 2 H q, pδr H , δθ H q λn N,H , v h ´Ě δr H ˝pϕ n H q ´1 ´Ě δθ H ˝pϕ n H q ´1 ^´I 3 ´Ę r n´1 H ˝pϕ n H q ´1¯¯ε ,n `ˆµ N,H , u n h ´9 r n´1 2 H ˝pϕ n H q ´1 ´9 θ n´1 2 H ˝pϕ n H q ´1 ^´I 3 ´Ę r n´1 H ˝pϕ n H q ´1¯˙ε ,n " 0.
(3.4) for all pv h , q h , µ N,H , δr H , θ H q P V h ˆQh ˆFn N,H ˆYH ˆYH . The energy stability of this method is stated in the following result.

Theorem 3.1. Let ␣ pu n h , p n h , λ n N,H , r n H , θ n H q ( ně1 Ă V h ˆQh ˆFn N,H
ˆYH ˆYH be given by problem (3.4). There holds

E n ď E 0 , @n ě 0 (3.5)
where the discrete mechanical energy E n of the system is defined by

E n def " ρ f 2 }u n h } 2 0,Ω `ρb 2 }9 r n H } 2 A,p0,Lq `ρb 2 } 9 θ n H } 2 I,p0,Lq `1 2 › › pr n H , θ n H q › › 2 b,H , with }p¨, ¨q} b,H " b a b H `p¨, ¨q, p¨, ¨q˘.
Proof. By testing (3.4) with pv h , q h , µ N,H , δr H , δθ H q " pu n h , p h , λ n N,H , 9 r q, we get

ρ f pB τ u n h , u n h q Ω `ρb ´AB τ 9 r n H , 9 r n´1 2 H ¯p0,Lq `ρb ˆIB τ 9 θ n H , 9 θ n´1 2 H ˙p0,Lq `af Ω,h p `un´1 h ; pu n h , p n h q, pu n h , p n h q ˘`a b H `pr n´1 2 H , θ n´1 2 H q, p9 r n´1 2 H , 9 θ n´1 2 H q ˘" 0, (3.6)
Using integration by parts and the boundary conditions of u n´1 h , u n h , we have

a f Ω,h pu n´1 h ; pu n h , p h q, pu n h , p h qq ě 2µ}εpu n h q} 2 0,Ω .
Moreover for any quantities v n , v n´1 and inner product p¨, ¨q, it holds

$ ' ' & ' ' % `vn ´vn´1 , v n ˘" 1 2 p}v n } 2 ´}v n´1 } 2 `}v n ´vn´1 } 2 q ě 1 2 B τ }v n } 2 , ˆvn ´vn´1 , v n `vn´1 2
˙" 1 2 B τ }v n } 2
where } ¨} def " a p¨, ¨q. We deduce that

B τ ˆρf 2 }u n h } 2 Ω `ρb 2 }9 r n H } 2 A,p0,Lq `ρb 2 } 9 θ n H } 2 I,p0,Lq `1 2 }r n H , θ n H } 2 b,H ˙ď0.
The estimate (3.5) then follows by summing over t0 . . . nu, which completes the proof.

Remark 3.2. Once again, to obtain the energy estimate and because we consider a linear Timoshenko beam formulation, we have to assume small rotational velocity which induces a mismatch between the structure velocity at the boundary and the velocity associated to the domain motion. It is nevertheless not an issue to derive the stability of the scheme due to the fictitious domain approach. It is not the case for the discrete counterpart of the ALE formulation where the velocities should match at the interface and where one should either consider a linear reconstruction of the geometry or a non linear beam equation.

Inf-sup stability. A crucial point when dealing with Lagrange multipliers is to ensure that the inf-sup condition associated to the bilinear form pλ n N,H , v h q ε,n is satisfied on F n N,H ˆVh . In [START_REF] Kuchta | Analysis and approximation of mixed-dimensional pdes on 3d-1d domains coupled with lagrange multipliers[END_REF], the authors establish the inf-sup condition for the case N " 0 under some conditions on the triangulations of Ω, p0, Lq and Σ n ε . To deal with general meshes, they propose a discretization of p0, Lq depending on the intersection of the fluid elements and the beam centerline. With the the addition of a stabilization term, they established the inf-sup condition on F n 0,H ˆVh where the coefficients a k,H , b k,H introduced in (3.1) and defining the space F n 0,H belong to X 0 h p0, Lq instead of X 1 h p0, Lq. If N ą 0, the authors of [52] considered a toy two-dimensional problem and proved a stability result on F n N,H ˆ`X k h pΩq X V ˘, provided N ď k. In the same framework, they also introduced a general stabilization term, which in 3D would take the form (see also [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF])

s Bω n ε ,h pλ n N,H , µ N,H q " γh `λn N,H , µ N,H ˘ε,n ,
and proved unconditional stability for any pair pλ n N,H , v h q P F n N,H ˆ`X k h pΩq X 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 (3.3) of Algorithm 1. For the sake of conciseness, we only give the expressions of the algebraic counterparts of the bilinear form `λn N,H , v h ˘ε,n . The evaluation of the remaining coupling terms can follows similarly. Let tψ f i u 1ďiďN h and tψ b j u 1ďjďN H be the basis functions of V h and Y H , respectively, and let tx i u 1ďiďN h and ts j u 1ďjďN H be the corresponding nodes. We hence have

`λn

N,H , v h ˘ε,n "

N h ÿ i"1 v h px i qpλ n N,H , ψ f i q ε,n .
Using the definition of F n N,H given in (3.1) with tensorized notations, we have

`λn N,H , v h ˘ε,n " N h ÿ i"1 v h px i q ˆˆa 0,H `N ÿ k"0 pa k,H c k `bk,H s k q ˙˝pϕ n H q ´1 , ψ f i ˙ε,n ,
`λn N,H , v h ˘ε,n " N h ÿ i"1 v h px i q N ÿ k"0 N H ÿ j"1 " a 0,H ps j q ´ψb j ˝pϕ n H q ´1, ψ f i ¯ε,n `ak,H ps j q ´´ψ b j b c k ¯˝pϕ n H q ´1 , ψ f i ¯ε,n `bk,H ps j q ´´ψ b j b s k ¯˝pϕ n H q ´1 , ψ f i ¯ε,n ȷ .
Let pB N q T P R N H p2N `1qˆN h denote the algebraic counterpart of the bilinar term `λn N,H , v ˘ε,n , we get

pB N q il " $ ' ' ' ' ' & ' ' ' ' ' % ´ψb j ˝pϕ n H q ´1 , ψ f i ¯ε,n , l " j ¨p2N `1q, ´´ψ b j b c k ¯˝pϕ n H q ´1 , ψ f i ¯ε,n , l " j ¨p2N `1q `2k `1, ´´ψ b j b s k ¯˝pϕ n H q ´1 , ψ f i ¯ε,n , l " j ¨p2N `1q `2k `2. @ k P t0 . . . N u, (3.7)
In order to evaluate the coefficients (3.7), we need to determine the value of the tensorized basis functions pψ b j b c k q ˝pϕ n H q ´1, pψ b j b s k q ˝pϕ n H q ´1 on Σ n ε . To do so, we first mesh p Σ ε in the reference configuration and evaluate them at the mesh points. We then transport their values to the current configuration via the deformation map ϕ n H . Integration is achieved through P 1 -interpolation in the deformed mesh by localizing the quadrature points in the fluid mesh. Note that we choose a mesh of p Σ ε fine enough so that the error due to the 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 (3.3) based on [START_REF] Fernández | A Newton method using exact Jacobians for solving fluid-structure coupling[END_REF]. To this purpose, we introduce U n , P n , L n , R n and Θ n as the arrays of degrees of freedom associated with u n h , p n h , λ n N,H , r n H and θ n H respectively. Taking pδr H , δθ H q " p0, 0q in (3.3), we recover the Navier-Stokes equations with Dirichlet boundary condition on Σ n ε . This problem can be expressed in matrix form as follows:

» - A f C T ´pB N q T ´C S 0 B N 0 0 fi fl » - U n P n L n fi fl " » - b n´1 0 B N U b pR n , Θ n q fi fl , (3.8)
with

A f def " ρ f τ M f `Kf , b n´1 def " ρ f τ M f U n´1 , U b pR, Θq def " 1 τ `ϕpR, Θq ´ϕpR n´1 , Θ n´1 q ˘.
Here, the matrices M f ,

" K f C T ´C S
ȷ and B N stand for the algebraic counterpart of

pu n h , v h q Ω , a f Ω,h `un´1 h ; pu n h , p n h q, pv h , q h q ˘, `λn N,H , v ˘ε,n .
With a slight abuse of notation, the operator ϕ allows us to compute the beam deformation given the displacement R n of the centerline of the beam and its rotation vector Θ n . To the matrix system (3.8), we can associate the fluid solver operator F n : R N b Ñ R N f given by

F n ˆ"R n Θ n ȷ˙" L n .
In a similar manner, testing (3.3) with pv h , q h , µ N,H q " p0, 0, 0q yields the linear Timoshenko supproblem with fluid source terms. The equivalent matrix system writes as follows:

A b « R n´1 2 Θ n´1 2 ff " J n´1 `" F n M n ȷ (3.9)
with R n " 2R n´1 2 ´Rn´1 , Θ n " 2Θ n´1 2 ´Θn´1 and where

A b def " 4ρ b τ 2 M b `Kb , J n´1 def " ρ b τ 2 M b ˜4 " R n´1 Θ n´1 ȷ `2τ « 9 R n´1 9 Θ n´1 ff¸, F n " `EN f ˘T L n , M n " `EN m ˘T L n , 9 R n´1 2 " 1 τ pR n ´Rn´1 q, 9 Θ n´1 2 " 1 τ pΘ n ´Θn´1 q.
Here, the matrices M b , K b , E N f , E N m stand for the algebraic counterpart of the bilinear forms

pAr n H , δr H q `pIθ n H , δθ H q , a b H `pr n H , θ n H q, pδr H , δθ H q ˘, ´λn N,H , s δr ˝pϕ n H q ´1¯, ´λn N,H , Ď δθ ˝pϕ n H q ´1 ^`I 3 ´Ę r n´1 H ˝pϕ n H q ´1 ˘¯.
The corresponding solid solution operator S n : R N f Ñ R N b is defiend as

S n ¨» - U n P n L n fi fl '" " R n Θ n ȷ .
By composition, solving the coupled problem (3.3) is equivalent to compute the roots of the following system

R n ˆ"R n Θ n ȷ˙d ef " " R n Θ n ȷ ´Sn ˝Fn ˆ"R n Θ n ȷ˙" 0, (3.10)
which can be iteratively approximated using Newton's method below, where J n denotes the Jacobian of the residual operator R n . As in [START_REF] Fernández | A Newton method using exact Jacobians for solving fluid-structure coupling[END_REF], the tangent system (3.11) can be solved in a matrix-free fashion via GMRES iterations, which only invoke the action of the operator 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 (3.3).

Newton algorithm applied to (3.10)

Choose " R n Θ n ȷ P R N b while › › › › R n ˆ"R n Θ n ȷ˙› › › › ě tol do Evaluate the fluid operator " U n P n ȷ " F n ˆ"R n Θ n ȷĖ valuate the solid operator « Rn Θn ff " S n ˆ"U n P n ȷĖ valuate the residual R n ˆ"R n Θ n ȷ˙" " R n Θ n ȷ ´« Rn Θn ff Solve tangent problem J n ˆ"R n Θ n ȷ˙" δR δΘ ȷ " ´Rn ˆ"R n Θ n ȷ˙( 3.11) Update " R n Θ n ȷ " " R n Θ n ȷ `"δR δΘ ȷ end while

Comparison with other discrete formulations

For the sake of completeness, in this section we provide the discrete formulations of the coupled problems (2.4), (2.7), (2.20) and (2.4), (2.11), (2.21) introduced above. These numerical methods will be used in the numerical experiments of Section 3 for comparison purposes.

Coupled problem with ALE formalism in the fluid. Since for the the coupled problem (2.4), (2.7), (2.20) the kinematic coupling is enforced in a strong fashion, we consider a triangulation of the reference fluid domain p Ω f ε which is fitted to the triangulation of the interface p Σ ε . We set X 1 h p p Σ ε q as the trace space of X 1 h p p Ω f ε q and Π b h the corresponding Lagrange interpolation operator onto X 1 h p p Σ ε q. The geometric coupling condition (2.20) 1 is treated explicitly. For a given displacement r n´1 H P Y h and rotation vector θ n´1 H P Y h at time step n ´1, we define the ALE map A n h as

ϕ n H " pI 3 ¨ez q e z `Ę r n´1 H `Ğ Λ n´1 H `I3 ´pI 3 ¨ez q e z ˘in p ω ε , A n h " L h `Πb h ϕ n H ˘, with Λ n´1 H
given from θ n´1 H , using (2.3), and where L h denotes a given discrete lifting operator (e.g., an harmonic extension operator). The one-to-one mapping A n h allows us to consider as fluid unknowns ûn

h def " u n h ˝An h , pn h def " p n h ˝An h instead of u n h " ûn h ˝pϕ n H q ´1 and p n h " pn h ˝pϕ n H q ´1.
We then introduce the following finite element space

V f h def " X 1 h p p Ω f ε q X H 1 BΩ p p Ω f ε q and Q f h def " Q 1 h p p Ω f ε q X L 2 p p Ω f ε q
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.4), (2.7), (2.20).

For n ě 1,

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

$ & % ϕ n H " pI 3 ¨ez q e z `Ę r n´1 H `Ğ Λ n´1 H `I3 ´pI 3 ¨ez q e z ˘in p ω ε , A n h " L h `Πb h ϕ n H ˘, Ω f,n ε " A n h ´p Ω f ε ¯, p w f,n h " B τ A n h , with Λ n´1 H given from θ n´1 H using (2.

3).

Step 2:

Find pp u n h , p p n h , r n H , θ n H q P V f h ˆQf h ˆYH ˆYH with 9 r n´1 2 H " B τ r n H , 9 θ n´1 2 H " B τ θ n H and such that $ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % p u n h " Π b h B τ ϕ n H on p Σ ε , ρ f τ " pu n h , v h q Ω f,n ε ´`u n´1 h , v h ˘Ωf,n´1 ε ı ´ρf ´p∇ ¨wf,n h qu n h , v h ¯Ωf,n ε `aΩ f,n ε ,h `un´1 h ´wf,n h ; pu n h , p n h q , pv h , q h q ˘`ρ b pAB τ 9 r n H , δr H q p0,Lq `ρb `IB τ 9 θ n H , δθ H ˘p0,Lq `ab H `pr n´1 2 H , θ n´1 2 H q, pδr H , δθ H q ˘" 0 (3.12) for all pp v h , p q h , µ, δr H , θ H q P V f h ˆQf h ˆYH ˆYH with p v h " Π b h `Ě δr H `Ě δθ H ^pϕ n H ´rH n q ˘on p Σ ε .
One limitation of Algorithm 2 lies in the capability of the discrete lifting operator 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 coupled problem (2.4), (2.11), (2.21). As in the previous discrete formulations, the geometrical coupling is treated in an explicit fashion as follows:

ϕ n H " r n´1 H , B n " ϕ n H `t0u ˆt0u ˆp0, Lq ˘.
Then the discrete problem is similar to the one given in Algorithm 1 but with the Lagrange multipliers belonging to Y n H " X 1 h pB n q instead of F n N,H and the integrals over Σ n ε replaced with integrals over B n . The details of this formulation can be found in Algorithm 3.

Algorithm 3 Discrete formulation of the coupled problem (2.4), (2.11), (2.21).

For n ě 1,

Step 1: Update of the beam centerline location

ϕ n H " r n´1 H , B n " ϕ n H `t0u ˆt0u ˆp0, Lq ˘.
Step 2:

Find pu n h , p n h , λ n H , r n H , θ n H q P V h ˆQh ˆYn H ˆYH ˆYH , with 9 r n´1 2 H " B τ r n H and 9 θ n´1 2 H " B τ θ n H , such that ρ f pB τ u n h , v h q Ω `af Ω,h `un´1 h ; pu n h , p n h q , pv h , q h q ρb pAB τ 9 r n H , δr H q p0,Lq `ρb ´IB τ 9 θ n H , δθ H ¯p0,Lq `ab H `pr n´1 2 H , θ n´1 2 H q, pδr H , δθ H q λn H , v h ´δr H ˝pϕ n H q ´1¯B n `´µ H , u n h ´Bτ r n H ˝pϕ n H q ´1¯B n (3.13)
for all pv h , q h , µ H , δr H , δθ H q P V h ˆQh ˆYn H ˆYH ˆYH .

Compared to Algorithm 1, Algorithm 3 involves the localization of the beam centerline mesh wihtin the fluid mesh, instead of the interfacial mesh, making the assembly of the matrix B N , introduced in (3.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-order 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.

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 def " 2u ref ρε{µ where u 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 [63].

Single beam immersed in a Stokesian flow

In the first numerical example, we consider the mechanical interaction between an incopmressible viscous fluid fluid and a single beam in a rectangular domain of dimensions 3ˆ1ˆ1 centered at p0.5, 0.5, 0q. 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 Γ 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.

The following parabolic velocity profile is enforced at the inlet boundary Γ in :

u in px, tq def " u ref ˜1 `cos ˆ2πt T ˙¸`1 ´pz ´1qpz `1q ˘ex , (4.1)
with final time T =0.06. The following beam parameters are considered: Young's modulus E=10 7 , radius ε=0.06 and density ρ b " 1. For the fluid, we take a density ρ f =1. We run simulations for different values of h{ε, 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{ε " 1{4, and we have h " 0.05 on the faces of the rectangular domain. The second type correspond to 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 p0, Lq with a uniform mesh such that H=0.05. The behavior of the coupled system is studied over the time interval r0, 0.06s with a time-step length of τ =5 ¨10 ´3.

l x l y l z Γ in u " u in

Impact of fluid mesh refinement

Case h{ε ! 1. We consider here the refined fluid mesh. In Figure 5, 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 µ and of the reference velocity u 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 8 norm with respect to the solution provided by Algorithm 2. Lastly in Figure 6, we present snapshots of the fluid velocity magnitude and pressure fields at the final time t " 0.06. µ " 1, u ref " 10. µ " 0. 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 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 Algorithm 2. The lower magnitude of the velocity and pressure is noticeable in Figure 6, which is also consistent with differences in displacement observed in Figure 5 and Table 1. As discussed in [START_REF] Hagmeyer | One-way coupled fluid-beam interaction: capturing the effect of embedded slender bodies on global fluid flow and vice versa[END_REF], 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{ε ! 1, this difference becomes important and then leads to a considerable deviation from the reference solution.

Case h{ε « 1. We now consider the uniform fluid mesh. In this configuration, the thickness of the beam (ε " 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.

The results for the displacement of the last beam cross-section are given in Figure 7, relative error in L 8 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 8. 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{ε ! 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{ε ! 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 7 shows that the solution for fluid velocity suffers from the same limitations as in the case h{ε ! 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 afine functions requires a fine mesh to be captured accurately.

Case h{ε " 1. Here, we still consider the uniform mesh. However, we decrease the value of ε, in particular, we consider ε P t0.06, 0.04, 0.02u. In Algorithm 2, as before, the mesh is refined around the beam to ensure that h{ε " 1{4, which resolves the interface regardless of the beam thickness ε. ε=0.06, E=10 To keep comparable amplitudes for the beam displacement independently of the beam thickness, the Young's modulus of the beam is modified accordingly with values E P t10 7 , 4 ¨10 7 , 4 ¨10 8 u. We fix µ=0.05 and u ref =100. The results obtained are reported in Figure 9 for the beam displacement, in Table 3 for the relative error in L 8 norm and in Figure 10 for the velocity and pressure fields. We can clearly observe that by reducing ε, 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{ε 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 reasonnably refined fluid mesh near the interface, while maintaining reduced order interface conditions in the spirit of Algorithm 2. 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 ě 1) in Algorithm 1. We consider homogeneous Neumann boundary conditions on Γ in and Γ out , homogeneous Dirichlet boundary conditions on the lower face of the domain, and u " u ref py ´0.5qp1 ´px ´1qpx `1qqe x on the three other lateral faces, with u 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 ε=0.04. The physical parameters are µ " 1, E=4 ¨10 5 and ρ b =1.

In Algorithms 1 and 2, we consider the same type of fluid meshes, refined around the beam (h{ε= 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 11, we report the time history of the rotation of the last cross-section around the z axis over the time interval r0, 0.06s. In Figure 12, we provide some snapshots of the beam deformed configuration at time t P t0, 0.02, 0.04, 0.06u. 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 12. 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.20), belong to the space F 1 .

Single beam immersed in a Navier-Stokes flow

In this section, we keep the geometric configuration considered in Section 4.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 ρ f =1 and a beam density ρ b =1. The thickness of the beam cross-section is set to ε=0.06 and the Reynolds number Re P t2.4, 120, 480u, with µ P t1, 0.1, 0.05u and u ref P t10, 50, 100u, respectively. To keep similar displacement amplitudes when varying the Reynolds number, the Young modulus of the beam is chosen as E P t10 7 , 1.6 ¨10 8 , 3 ¨10 8 u. We use a mesh that accurately resolves the beam interface such that h{ε " 1{4 in a region near the beam and h " 0.05 on the faces of the rectangular domain. Impact of the Reynolds number. We first investigate the influence of the Reynolds number on the accuracy of Algorithm 1. In Figure 13, we report the time history of the displacement of the lass crosssection of the beam with Re P t2.4, 120, 480u and E P t4 ¨10 relative error with respect to Algorithm 2. Note that in Figure 13, we introduced an additional algorithm, referred to as Algorithm 2-Ω in in which the fluid is solved in the whole domain Ω. 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 Ω, which introduces a non-physical flow in the part occupied by the beam with the subsequent additional artificial stress on the interface (see Remark 2.3). It is worth noting that, in the case of a Stokes flow, the numerical results of Section 4.1 indicate that the impact of such spurious stress on the accuracy of Algorithm 1 is rather limited. 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 Ω, so that the portion of the domain occupied by the beam is also filled with the fluid. This variant will be termed Algorithm 2-Ω. Figure 13 and Table 4 indicate that the displacement of the last cross-section of the beam obtained with Algorithm 1, with N P t1, 2u, is close to the one provided by Algorithm 2-Ω. Since the volume of the region occupied by the beam is proportional to ε 2 , one would expect that the perturbation induced by the artificial fluid decreases as ε tends to zero. We propose to exhibit this phenomenon by gradually decreasing the thickness of the beam, by taking ε P t0.1, 0.04, 0.02u while keeping the same fluid parameters as in the case with Re=480 and ε " 0.06. To keep similar amplitudes of the beam deflection, we take the Young's modulus E P t10 8 , 1.2¨10 9 , 1.2¨10 10 u. In Figure 14 we report the time history of the displacement of the last cross-section of the beam and in Table 5 the L 8 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 15, we provide snapshots of the magnitude of the velocity and the pressure fields obtained with Re=480 and a constant velocity inlet u in " 2u ref e x . We can clearly observe that the pressure jump provided by Algorithm 2-Ω 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-Ω 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 [START_REF] Girault | A fictitious domain method with distributed multiplier for the stokes problem, applied nonlinear analysis[END_REF][START_REF] Boffi | The finite element immersed boundary method with distributed lagrange multiplier[END_REF]. 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 Ω f ε (see [START_REF] Corti | Low-order fictitious domain method with enhanced mass conservation for an interface stokes problem[END_REF][START_REF] Ohmori | Flux-free finite element method with Lagrange multipliers for two-fluid flows[END_REF]), viz., ż

BΩ f ε u ¨n " 0.
We denote by Algorithm 1 -MC the solution procedure obtained by combining Algorithm 1 with this approach. Figure 15 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.

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. The beams have a radius ε=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 (4.1). The beam is modeled by a non-linear co-rotational formulation that allows for large deformations (see [START_REF] Krenk | Non-linear modeling and analysis of solids and structures[END_REF]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 ρ b =1, respectively.

We first consider a Stokes flow with viscosity µ=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 16, we provide snapshots of the fluid velocity magnitude computed with Algorithm 1 at times t P t0.05, 0.1, 0.2, 0.3, 0.45, 0.6u. 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 µ=0.05, u ref =100 and ρ 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 17. 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.

Conclusion

We have proposed a new approach for the 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 with a fully resolved (ALE based) approach and an alternative reduced based method. Numerical evidence indicates that, at low Reynolds numbers, the method is accurate regardless of the ratio between the background mesh size 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. To handle contact interactions, we incorporate a contact term into Algorithm 1. We define this term in the context of two distinct beams ω n,1 ε and ω n,2 ε . Let for all x P Bω n,1 ε , y its projection on Bω n,2 ε using ray-tracing strategy, we introduce a gap function denoted as g which satisfies the relationship y " x`gn x , with n x representing the outer normal at point x on Bω n,1 ε . The contact term for the beam ω n, x " pϕ n H q ´1 pxq and p y " pϕ n H q ´1 pyq. Further details on the ray-tracing strategy and comprehensive insights into the contact methodology can be found in [START_REF] Chouly | Finite Element Approximation of Contact and Friction in Elasticity[END_REF]Chapter 11].

p ω ε 1 (

 1 psq (a) Reference configuration of the beam. εa ε p ω ε psq Bp ω ε psq se z p ω ε psq " εE ˆtsu (b) Elliptical cross-section of the beam in the reference configuration.E a c) Reference ellipse.

Figure 1 :

 1 Figure 1: Geometrical configuration of the beam and notations.

Figure 2 :

 2 Figure 2: The deformation map.

  , } ¨}A,p0,Lq def " b pA ¨, ¨qp0,Lq , } ¨}I,p0,Lq def " b pI ¨, ¨qp0,Lq , }¨, ¨}b def " b a b `p¨, ¨q, p¨, ¨q˘.

with c k pνq def "

 def cospkνq and s k pνq def " sinpkνq for 0 ď k ď N . Now, by expression a k,H , b k,H in terms of their basis functions ψ b j , we get

LFigure 3 :

 3 Figure 3: Geometric setting for the single beam immersed in a rectangular fluid domain.

( a )

 a Refined fluid mesh -Algorithms 1 and 3 (b) Uniform fluid mesh -Algorithms 1 and 3 (c) Refined fluid mesh -Algorithm 2

Figure 4 :

 4 Figure 4: Cut plane of the fluid mesh for different resolutions of the beam interface and ε " 0.06.

Figure 5 :

 5 Figure 5: Time history of the displacement of the last cross-section of the beam along the x-axis for several fluid viscosities µ and reference inlet velocities u ref , with h{ε ! 1.

(a) Algorithm 2 -

 2 Reference Solution. (b) Algorithm 1, N =0. (c) Algorithm 3.

Figure 6 :

 6 Figure 6: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time t=0.06 obtained with Algorithms 1-3, with h{ε ! 1.

Figure 7 :

 7 Figure 7: Time history of the displacement of the last cross-section of the beam along the x-axis for various fluid viscosity µ and reference inlet velocity u ref , with h{ε « 1.

(a) Algorithm 2 -

 2 Reference solution. (b) Algorithm 1 with N =0. (c) Algorithm 3.

Figure 8 :

 8 Figure 8: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time t=0.06 obtained with Algorithms 1-3, with h{ε « 1.

Figure 9 :

 9 Figure 9: Time history of the displacement of the last cross-section of the beam along the x-axis for various fluid viscosity µ and reference inlet velocity u ref -h{ε " 1.

(a) Algorithm 2 ,

 2 Reference solution. (b) Algorithm 1 with N =0. (c) Algorithm 3.

Figure 10 :

 10 Figure 10: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time t=0.06 obtained with Algorithms 1-3, with h{ε " 1.

4. 1 . 2 ." 1 Figure 11 :

 12111 Figure 11: Rotation of the last cross-section along the z-axis of the beam immersed in a shear flow.

Figure 12 :

 12 Figure 12: Snapshots of the beam deformation and the velocity magnitude induced by a shear flow, at time t P t0, 0.02, 0.04, 0.06u (from top to bottom).

Figure 13 :

 13 Figure 13: 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.

  Re=160, ε=0.02, E=1.2 ¨10 10 .

Figure 14 :

 14 Figure 14: Time history of the displacement of the last cross-section of the beam along the x-axis for various beam thickness ε.

  Algorithm 2-Ω. (c) Algorithm 1, N =0. (d) Algorithm 1, N =1.

  (e) Algorithm 1, N =2. (f) Algorithm 1-MC, N =2.

Figure 15 :

 15 Figure 15: Snapshots of the fluid velocity magnitude (left) and pressure (right) at time t=0.06 with Algorithm 2, Algorithm 2-Ω, Algorithm 1 and Algorithm 1-MC, with Re=480.

Figure 16 :

 16 Figure 16: Snapshots of the fluid velocity magnitude for multiple beams immersed in a Stokes flow with Algorithm 1 and N =0 at time t P t0.005, 0.01, 0.02, 0.03, 0.045, 0.06u (note that at time t " 0.03, the inlet velocity is null).

Figure 17 :

 17 Figure 17: Snapshots of the fluid magnitude fluid for multiple beams immersed in a Navier-Stokes flow with Algorithm 1and N =2 at time t P t0.005, 0.01, 0.02, 0.03, 0.045, 0.06u (note that at time t=0.03, the inlet velocity is null).

  2 pΣ ε ptqq a.e. t P R εpvqq Ω ´pp, ∇ ¨vq Ω `pq, ∇ ¨uq Ω .

				ànd
	ρ f pB t uptq, vq Ω	`af Ω `uptq; puptq, pptqq, pv, qq ˘´xλptq, vy Σεptq `xµ, uptq ´ub ptqy Σεptq " 0	(2.9)
	for all pv, q, µq P V ˆQ ˆH´1 2 pΣ ε ptqq, and where
	a f Ω `z; pu, pq, pv, qq	˘def " c Ω pz, u, vq `aΩ `pu, pq, pv, qq ˘, c Ω pz, u, vq	def " ρ f pz ¨∇u, vq Ω ,
	a Ω `pu, pq, pv, qq ˘def " 2µpεpuq,

  Lq

						`,
			pµ N , uptqq ε,t " `µN , u b ptq ˘ε,t @ µ N P F N ptq,	(2.14)
	where p¨, ¨qε,t denotes the re-scaled L 2 inner-product on Bω ε ptq:
			ż		ˆεb	˙´1
	pu, vq ε,t	def "	Σεptq	pu ¨vq h ˝ϕ´1 t , hps, νq	def "	a 2 cos 2 pνq `sin 2 pνq
	for all ps, νq P p0,

ˆp0, 2πq . Note that (2.14) reduces (2.8) by enforcing N `1 scalar constraints on each cross-section of the interface when (2.10) reduces (2.8) by neglecting the thickness of the beam and enforcing a single vector constraint on the beam centerline. The associated reduced weak formulation reads: we look for pu, p, λ N q such that puptq, pptq, λ N ptqq P V ˆQ ˆFN ptq a.e. t P R `and ρ f pB t uptq, vq Ω `af Ω `uptq; puptq, pptqq, pv, qq ˘´pλ N ptq, vq ε,t ``µ N , uptq ´ub ptq ˘ε,t " 0 (2.15)

  N , δr, δθq P V ˆQ ˆFN ptq ˆY ˆY. Here, the bilinear form a b represents the weak form of the beam elastic operator given by a b `pr, θq, pδr, δθq ˘def " pEB s θ, B s δθq p0,Lq `pG pB s r ´θ ^ez q , pB s δr ´δθ ^ez qq p0,Lq .

				tq,
	ρ f pB t u, vq Ω	`af Ω `u; pu, pq, pv, qq ρb
	`AB 2 t r, δr ˘p0,Lq	`ρb `IB 2 t θ, δθ ˘p0,Lq	`ab `pr, θq, pδr, δθq	λN (2.26)
	, v ´s δr ˝ϕ´1 t ´Ď δθ ˝ϕ´1 t ^`I 3 ´s r ˝ϕ´1 t ˘˘ε,t
	``µ N , u ´Ď B t r ˝ϕ´1 t ´Ď B t θ ˝ϕ´1 t ^pI 3 ´s r ˝ϕ´1 t q ˘ε,t " 0
	for all pv, q, µ

The energy stability of (2.26) is stated in the next result.

Theorem 2.4. Let pu, p, λ N , r, θq be a regular enough solution of the reduced coupled problem (2.26).

There holds

Eptq ď Ep0q @ t P R `,

  Discrete formulation of the reduced coupled model (2.4), (2.15), (2.23).

60], viz., a b H `pr H , θq , pr H , δr H q ˘def " `EB s θ H , B s δθ H ˘p0,Lq ``GΠ h pB s r n H ´θn H ^ez q, B s δr H ´δθ H ^ez ˘p0,Lq ,

where Π h denotes the L 2 projection onto X 0 h p0, Lq.

Algorithm 1

Table 1 :

 1 Relative error } `rpL, ¨q ´rref pL, ¨q˘¨e x}8{}r ref pL, ¨q ¨ex}8 where r ref pL, ¨q is the displacement of the last crosssection of the beam computed with Algorithm 2, with h{ε ! 1.

	Algorithm 1, N =0	2.84 ¨10 ´2	2.88 ¨10 ´2	3.01 ¨10	´2
	Algorithm 3	4.94 ¨10 ´1	5.46 ¨10 ´1	5.78 ¨10	´1

1, u ref " 50. µ " 0.05, u ref " 100.

  µ=1, u ref =10. µ=0.1, u ref =50. µ " 0.05, u ref " 100.

	Algorithm 1, N =0	6.88 ¨10 ´2	8.14 ¨10 ´2	9.39 ¨10	´2
	Algorithm 3	2.63 ¨10 ´1	2.71 ¨10 ´1	2.85 ¨10

´1

Table 2: Relative error } `rpL, ¨q ´rref pL, ¨q˘¨e x}8{}r ref pL, ¨q ¨ex}8 where r ref pL, ¨q is the displacement of the last crosssection of the beam computed with Algorithm 2, with h{ε « 1.

Table 3 :

 3 7 . ε=0.04, E=4 ¨10 7 . ε=0.02, E=4 ¨10 8 . Relative error } `rpL, ¨q ´rref pL, ¨q˘¨e x}8{}r ref pL, ¨q ¨ex}8 where r ref pL, ¨q is the displacement of the last crosssection of the beam computed with Algorithm 2, with h{ε " 1.

	Algorithm 1, N " 0	9.39 ¨10 ´2	1.13 ¨10 ´1	2.20 ¨10	´1
	Algorithm 3	2.85 ¨10 ´1	1.45 ¨10 ´1	2.03 ¨10

´1

Table 4 :

 4 7 , 4 ¨10 7 , 4 ¨10 8 u, and in Table 4 the L 8 Relative error } `rpL, ¨q ´rref pL, ¨q˘¨e x}8{}r ref pL, ¨q ¨ex}8 where r ref pL, ¨q is the displacement of the last crosssection of the beam computed with Algorithm 2 or Algorithm 2-Ω.

		Re=2.4, E=4 ¨10 7	Re=120, E=1.6 ¨10 8	Re=480, E=3 ¨10 8
		Algorithm 2 Algorithm 2-Ω Algorithm 2 Algorithm 2-Ω Algorithm 2 Algorithm 2-Ω
	Algorithm 1, N =0	7.05 ¨10 ´2	6.49 ¨10 ´2	5.17 ¨10 ´1	3.41 ¨10 ´1	4.31 ¨10 ´1	2.90 ¨10	´1
	Algorithm 1, N =1	6.95 ¨10 ´2	6.38 ¨10 ´2	2.67 ¨10 ´1	1.27 ¨10 ´1	2.06 ¨10 ´1	8.75 ¨10	´2
	Algorithm 1, N =2	7.89 ¨10 ´2	7.30 ¨10 ´2	2.02 ¨10 ´1	7.66 ¨10 ´2	2.00 ¨10 ´1	8.43 ¨10

´2

8

  Re=480, ε=0.06, E=1.2 ¨10 9 Re=160, ε=0.02, E=1.2 ¨10 10

	Algorithm 1, N =0	6.22 ¨10 ´1	4.31 ¨10 ´1	3.23 ¨10 ´1
	Algorithm 1, N =1	2.61 ¨10 ´1	2.06 ¨10 ´1	1.96 ¨10 ´1
	Algorithm 1, N =2	2.54 ¨10 ´1	2.00 ¨10 ´1	1.63 ¨10 ´1

Table 5 :

 5 Relative error } `rpL, ¨q ´rref pL, ¨q˘¨e x}8{}r ref pL, ¨q ¨ex}8 where r ref pL, ¨q is the displacement of the last crosssection of the beam computed with Algorithm 2.