{'fraction_non_alphanumeric': 0.053596988814080475, 'fraction_numerical': 0.016593330859354292, 'mean_word_length': 4.009028146574615, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A time-dependent FEM-BEM coupling method for fluid-structure interaction in 3d

We consider the well-posedness and a priori error estimates of a 3d FEM-BEM coupling method for fluid-structure interaction in the time domain. For an elastic body immersed in a fluid, the exterior linear wave equation for the fluid is reduced to an integral equation on the boundary involving the Poincar\'{e}-Steklov operator. The resulting problem is solved using a Galerkin boundary element method in the time domain, coupled to a finite element method for the Lam\'{e} equation inside the elastic body. Based on ideas from the time-independent coupling formulation, we obtain an a priori error estimate and discuss the implementation of the proposed method. Numerical experiments illustrate the performance of our scheme for model problems.


Introduction
Coupled finite and boundary element procedures provide an efficient and extensively investigated tool for the numerical solution of elliptic interface and contact problems, particularly in unbounded domains [21,32]. On the other hand, much of the current interest in boundary element procedures focuses on hyperbolic problems in the time domain, both on Galerkin methods and convolution quadrature [2,9,14,27,31].
To compute the scattering of time-dependent waves by a bounded, penetrable obstacle, the coupling of time domain finite elements (FEM) and boundary elements (BEM) becomes relevant. The recent mathematical analysis of FEM-BEM coupling in the time domain was initiated in [1], coupling discontinuous finite elements to Galerkin boundary elements for the 3d wave equation. A general analysis of the coupling between different discretizations for acoustic wave equations was provided in [6], with a focus on convolution quadrature. Since then, FEM-BEM coupling for convolution quadrature methods has been applied in a variety of applications in 2 dimensions, such as fluid-structure and fluid-thermoelastic problems, as well as nonlinear elastic problems involving piezoelectric scatterers [24,25,26,30]. For time-dependent Galerkin methods and their application to 3d problems, on the other hand, much less is known. In addition to [1], previous related work includes the energy-based formulations investigated by Aimi and collaborators for wave-wave coupling in 3d multidomains and layered media [3,4].
In this article we study a simple space-time Galerkin FEM-BEM coupling method for fluid-structure interaction, describing the transient scattering of waves in an inviscid homogeneous fluid by an elastic obstacle. Based on ideas from time-independent coupling formulations [8,10] and the analysis in the frequency domain [24], we present a basic a priori error estimate for a space-time Galerkin approximation in anisotropic Sobolev spaces [5]. We discuss in detail the numerical implemention of our proposed coupling method in 3d. Numerical experiments for model problems illustrate the performance of the scheme.
To describe the results of this article in more detail, recall the equations for an elastic body submersed in a fluid. Let Ω ⊂ R 3 be a bounded Lipschitz domain, and Ω c = R 3 Ω. The elastic deformation u in Ω is described by the Lamé operator ∆ * u = µ∆u + (λ + µ)∇(divu), with Lamé constants µ ≥ 0 and λ, such that 3λ + 2µ ≥ 0. The deformation is coupled to the wave equation in Ω c , leading to the coupled interface problem v − c 2 ∆v = 0, (x, t) ∈ Ω c × (0, ∞), u(x, 0) =u(x, 0) = 0 in Ω , for a given incident wave v inc in Ω c . Here, the stress is given in terms of the deformation asσ(u) = (λdivu)I + 2µε(u), ε(u) = 1 2 ((∇u) + (∇u) T ), with I the identity matrix. Time derivatives are denoted by a dot, and n is the outwardpointing unit normal vector to ∂Ω. The Neumann trace on Γ from the exterior domain Ω c is denoted by ∂ + n , while ∂ − n is the corresponding Neumann trace from the interior domain Ω. We choose units in which c = ρ 1 = ρ 2 = 1.
To solve this interface problem numerically, we use the Poincaré-Steklov operator for the exterior wave equation to reformulate it as a coupled domain / boundary integral equation in Ω and Γ. The Poincaré-Steklov operator is expressed in terms of layer potentials for the wave equation, as known for timeindependent symmetric FEM-BEM coupling methods. The resulting spacetime weak formulation is approximated using finite elements in Ω and Galerkin boundary elements on Γ, based on tensor products of piecewise polynomial functions on a quasi-uniform mesh in space and a uniform mesh in time. Our a priori estimates assure convergence. We discuss a numerical implementation in detail and study the numerical performance of the method.
The article is organized as follows: Section 2 reformulates the coupled problem (1) as a domain / boundary integral equation in Ω and Γ and discusses its discretization and well-posedness. It provides the basis for the derivation of an a priori error estimate in Section 3. Numerical examples and related algorithmic considerations are the content of Section 4. Two appendices discuss boundary integral operators for the wave equation and the appropriate spacetime anisotropic Sobolev spaces, as well as detailed algorithmic aspects of the proposed scheme. Notation: To simplify notation, we will write f ≲ g, if there exists a constant C > 0 independent of the arguments of the functions f and g such that f ≤ Cg.
We will write f ≲ σ g, if C may depend on a parameter σ. For a function v on R n ∖ Γ, we denote by v ± Γ the trace of v on Γ from the exterior domain Ω c , resp. from the interior domain Ω.

Weak formulation and FEM-BEM coupling
Recall that the fluid-structure interaction problem (1) is well-posed [11]: , which depends continuously on the data. See Appendix A for the definitions of the relevant space-time anisotropic Sobolev spaces, depending on a weight parameter σ > 0.
In this section we reformulate problem (1) as a coupled domain / boundary integral equation and propose a finite element / boundary element coupling method for its numerical solution.
The key ingredient to reduce equation (1a) in Ω c to Γ is the retarded Poincaré-Steklov operator S on Γ, defined as Sv + Γ = ∂ + n v for a solution v of (1a).
Using S, equation (1d) becomes: To compute the Poincaré-Steklov operator S, we use the following formulation in terms of boundary integral operators [21], see Appendix A: In terms of φ = v + Γ and an auxiliary variable Problem (1) is therefore equivalent to the following system on Ω and Γ: The solution v in Ω c is recovered from the representation formula v = Dφ − Sλ, using the single and double layer potentials from Appendix A.
Adding (4) and the weak formulations of (3c) and (3d), one obtains a weak formulation of (3) inX ∶= H 1 σ (R + , H 1 (Ω)) 3 ×H 1 σ (R + , H 1 2 (Γ))×H 1 σ (R + , H −1 2 (Γ)): ) be the conforming discretization spaces from Appendix A in Ω, resp. Γ, based on tensor products of piecewise polynomial functions on a quasi-uniform mesh in space and a uniform mesh in time. Let (△t) . Note that the discretization order is higher in time than in space, in order to be conforming. This corresponds to the loss of one time derivative for the boundary integral operators in Theorem 6, a well-known sub-optimal aspect of the standard functional analytic framework of space-time anisotropic Sobolev spaces.
Practical computations use σ = 0. See [27] for a detailed analysis of the role of the weight σ.
In order to prove the well-posedness of the discrete formulation, we show the equivalence to a coercive formulation. Proposition 2. Let (u, φ, λ) ∈X h,(∆t) be a solution to the weak formulation (6).
Finally, from (6) with Since v satisfies the wave equation (8b), we get (7) from the representation formula: Combining (8d) with the jump relations for ∂ + n v and v Γ , as well as setting φ in the equation (8a), we obtain (6).
where (△t) . Now for all w ∈Z h,(△t) using the second relation in (8c) and ⟨⟦∂ n v⟧, w − Γ ⟩ Γ×R + = 0, Green's formula and (8b) lead to Conversely, assume (11) holds. Using (11) for a test function w ∈Z h,(△t) , with compact support in R 3 Γ we obtain the equation (12): Using integration by parts: Since v satisfies the wave equation on the support of w in R 3 Γ.
An analogous assertion to Propositions 2 and 3 holds for the continuous problem, instead of the finite element discretization.
We now aim to prove coercivity of the problem (11) for (u, v) in a suitable norm, defined as:

Note that
Using integration by parts in time, the zero initial condition and σ > 0, we obtain the coercivity estimate: ,Ω c . This implies, in particular, uniqueness of the solution (11) and therefore also the solutions to (5) and (6).

A priori error estimate
We state an a priori error estimate: Proof. Let (u, φ, λ) ∈X satisfy (5) and (u h , φ h , λ h ) ∈X h,∆t satisfy (6). Then for all (w,φ,λ) ∈X h,∆t By definition of the single and double layer potentials S and D, the functions v, v h ,r all satisfy the wave equation in R 3 Γ. Further, from Propo- Using the definition of A and Green's theorem, we find We estimate the individual terms. Using Young's inequality we have for ǫ > 0 Using the inverse estimate as in (3.182) in [20] φ 1,1 2,Γ ≲ 1 ∆t φ 0,1 2,Γ , we further estimate with the mapping properties of the integral operators For the fourth term, we get: For the last term, analogously the trace theorem and the inverse estimate show: Moving the terms with positive powers of ǫ to the left hand side and choosing a fixed, sufficiently small ǫ > 0 depending on σ, we conclude: and therefore This proves the assertion.

Numerical results
This section presents numerical results for the fluid structure-interaction problem given by (5), in 3d. While finite element discretizations of fluidstructure interaction have attracted significant recent interest, coupled finite and boundary element procedures in the time domain are only beginning to be explored. In 2d, numerical results have been presented in [24], based on time discretization by convolution quadrature of the boundary integral operators. A similar appooach has been demonstrated for the interaction of waves with a thermoelastic solid in 2d [26]. The authors are not aware of any related numerical results in the mathematical literature based on time-domain Galerkin boundary element methods, as presented in this work. On the other hand, such methods are now actively being studied for wave-wave interaction in 2d and 3d, as in [1] and [3].
In the numerical experiments for Problem (1), the variational formulation (5) is solved, for σ = 0, by choosing as ansatz function in the interior domain Ω where {η i h } denotes the basis of piecewise linear hat functions for V 1 h (Ω) and , with H the Heaviside function. The test functions are given byẇ for l = 1, . . . , N o , n = 1, . . . , N t and µ = 1, 2, 3. The ansatz functions on Γ are taken as and where {ξ i h } denotes the basis of piecewise linear hat functions for V 1 h . The corresponding test functions are chosen aṡ for 1 ≤ n ≤ N t and 1 ≤ j ≤ N s ′ . The resulting discretization of the Poincaré-Steklov operator has been tested in [15,16], and corresponding results are obtained for more natural discretizations with piecewise constant λ h,△t . Piecewise linear and higher order test functions are considered in [19]. As shown in Appendix B, this discretization of (5) leads to a time-stepping scheme, which solves a system of the following structure in each time step ≥ 3: The system in the first two time steps is similar, see Appendix B. We solve this system repeatedly until our desired time step N t is reached.
Example. Let Ω = [−1, 1] 3 . Using the discretization described above, we compute the solutions to the discrete system (5) up to time T = 4 for data v inc corresponding to the exact solution We use uniform discretizations by tetrahedra as depicted in Figure 1 and a time step ∆t such that ∆t h ≃ 0.1414. We denote the number of grid points on an edge of the cube by n + 1. The finest mesh is then given by n = 24 and consists of 69120 tetrahedra, corresponding to ∆t = 0.01667. The convergence of the numerical solution to the exact solution is studied as the mesh is refined, and we measure the error in terms of the L 2 -norm in space, resp. space-time. . Ω c , is illustrated in Figure 3, which plots the L 2 -norms of u in Ω, resp. of v Γ = ϕ on Γ, as a function of time for n = 2, 8, 24. These norms are approximated from the solution vectors u n and ϕ n of the discrete system using a trapezoidal rule for the integrals. The L 2 error as a function of time is shown in Figure 4, corresponding to the numerical solutions depicted in Figure 3. All plots show excellent approximation of the simple behavior of the solution for short times and a monotonous convergence on the whole time interval. Note that error in Figure 4 does not seem to grow with time, as expected for a variational method. Figure 5 considers the convergence of the numerical solutions up to n = 24 in terms of the mesh size h. It depicts the L 2 (Ω × [0, T ])-norm of the error in u, as well as the L 2 (Γ × [0, T ])-norm of the error in ϕ. Similar convergence rates are obtained for u and ϕ in these L 2 -norms: 0.76 for u, 0.80 for ϕ. Note that based on Theorem 4 and the trace theorem for Sobolev spaces, one might naively expect slower convergence in ϕ than in u (by a difference of the rates 0.5) in the norms used here. However, as shown in [28] for time-independent FEM-BEM coupling, under mild regularity assumptions the BEM solution ϕ converges at a rate 0.5 faster than predicted from a joint estimate for (u, ϕ) as in Theorem 4. This exactly cancels the above difference of rates and leads to identical convergences rates for u, ϕ in the spacetime L 2 -norms on Ω, resp. Γ. The identical observed convergence rates for u, ϕ are therefore expected and in line with those known for FEM-BEM coupling in time-independent problems [28].

Appendix A. Integral operators and finite-boundary elements
We recall basic definitions and properties of boundary integral operators for the wave equation from [31], as well as from [9,22].
Let Γ be the boundary of a polyhedral domain Ω in R 3 , consisting of curved, polygonal boundary faces. In R 3 Γ, a solution v to the homogeneous wave equation may be represented in terms of the jump of the Dirichlet and Neumann data across Γ: v = Dφ − Sλ. Here for x ∈ R 3 Γ and t ≥ 0 are the single, resp. double layer potential for the wave equation defined from the fundamental solution G(τ, x, y) = δ(τ − x−y ) 4π x−y . The coupling method presented in this article relies on the resulting boundary integral operators on Γ × (0, ∞).
They are studied in space-time anisotropic Sobolev spaces H r σ (R + , H s (Γ)) [22]. To define an explicit scale of Sobolev norms, fix a partition of unity α i subordinate to a covering of Γ by open sets B i and diffeomorphisms φ i mapping each B i into the unit cube ⊂ R 2 . They induce a family of norms from R 2 : .
When s ≤ 1 one can show that the spaces are independent of the choice of α i and φ i .
In a bounded Lipschitz domain Ω ⊂ R d , we define space-time anisotropic Sobolev spaces H r σ (R + , H s (Ω)) analogously to above, starting from the standard Sobolev spaces H s (Ω) with norm u s, We state the mapping properties of the boundary integral operators, see e.g. [9,22]: Theorem 6. The following operators are continuous for r ∈ R, σ > 0: By a fundamental observation of Bamberger and Ha-Duong [5], V ∂ t satisfies a coercivity estimate in the norm of H 0 σ (R + , H − 1 2 (Γ)), provided σ > 0: From the mapping properties of Theorem 6 one also has the continuity of the bilinear form associated to V ∂ t in a bigger norm: ⟨V ψ,ψ⟩ ≲ ψ 2 1,− 1 2 ,Γ . Similar estimates hold for W ∂ t : φ 2 0, 1 2 ,Γ ≲ σ ⟨W φ,φ⟩ ≲ φ 2 1, 1 2 ,Γ . Proofs and further information may be found in [22]. For sufficiently regular φ and λ, the following jump relations hold [22]: We consider space-time discretizations based on tensor products of piecewise polynomials: For simplicity, we assume that Ω is a polygonal domain, with a quasi-uniform triangulation T Ω = {T 1 , ⋯, T No }by N o tetrahedra. The induced quasi-uniform triangulation of the boundary Γ, T S = {∆ 1 , ⋯, ∆ N s ′ }, should consist of closed triangular faces ∆ i , such that each ∆ i is a face of one T j and at most one face of T j is contained in Γ. We The space-time approximation spaces are given by tensor products of the approximation spaces in space and time, V q h and V p ∆t , associated to the spacetime meshes T Ω,T = T Ω × T T , respectively T S,T = T S × T T . We write

Appendix B. Discretization and MOT-Algorithm
This appendix discusses the details of the discretization (6), where we set σ = 0. The resulting formulas for the entries of the Galerkin matrices reduce their assembly to numerical quadratures over certain light cones E l below. This structure is crucial for the practical implementation in standard time-domain boundary element codes, see [14,33], as well as the recent Ph.D. thesis [29] of the second author.
We choose the finite element ansatz and test functions as in (13) - (18). For the discretization of the boundary integral operators, we begin with the retarded hypersingular operator.
We choose the ansatz function as in (15) and the test function as in (17).
Here, for l ∈ N 0 we define the light cone The matrix W n−m is therefore a sum of integrals over the three light cones E n−m , E n−m−1 and E n−m−2 .
We now consider the discretization of the single layer potential. For the ansatz function we choose (16) and as test function we choose (18).
After some computations, we obtain We next consider the retarded adjoint double layer potential: The related term for the mass matrix is given by It remains to consider the coupling contributions. For j = 1, . . . , N s ′ ⟨u h,△t Γ ⋅ n,ẇ h,∆t ⟩ Γ×R+ = For the second coupling term: For completeness we mention the right hand side: Set h =v inc n and g = ∂ + n v inc . We approximate the time integral by the trapezoidal rule, so that: − ⟨v inc n,ẇ Γ ⟩ Γ×R+ = −⟨h,ẇ Γ ⟩ Γ×R+ = − (△t) 2 Γ (h n +h n−1 )η j h (x)e µ ds x =∶ H n + H n−1 ⟨∂ + n v inc ,ẇ⟩ Γ = ⟨g,ẇ⟩ Γ = (△t) 2 Γ (g n + g n−1 )ξ j h (x)ds x =∶ G n + G n−1 , where h n = h(x, t n ) and g n = g(x, t n ). The matrices W k , K k , K T k , V k vanish if the index k is negative. Therefore we get for (B.1) in the first time step (n = 1): Note the zero block in [0, −RIn x ], as u h,∆t Γ only depends on the values in the nodes on the boundary Γ. Similarly, one obtains a zero block in [0, n x RI] T , corresponding to the vanishing contribution of nodes in the interior of Ω to the trace of the test function w h,∆t Γ . For the second time step (n = 2) we obtain using u 1 , ϕ 1 and λ 1 from above. For later time steps n ≥ 3 we conclude: This system is solved repeatedly until reaching time step N t ≥ 3.