Convergence Analysis of an Unfitted Mesh Semi-implicit Coupling Scheme for Incompressible Fluid-Structure Interaction

The complete numerical analysis of time splitting schemes which avoid strong coupling has rarely been addressed in the literature of unfitted mesh methods for incompressible fluid-structure interaction. In this paper, an error analysis of the semi-implicit scheme recently reported in Fernández and Gerosa (Int. J. Numer. Methods Eng. 122, 5384–5408, 2021) is performed for a linear fluid-structure interaction system. The analysis shows that, under a hyperbolic-CFL condition, the leading term in the energy error scales as O(hr− 1/2), where r = 1,2 stands for the extrapolation order of the solid velocity in the viscous fluid substep. The theoretical findings are illustrated via numerical experiments which show, in particular, that the considered method avoids the spatial non-uniformity issues of standard loosely coupled schemes and that it delivers practically the same accuracy as the strongly coupled scheme.


Introduction
This paper is devoted to the error analysis of a numerical method for a linear fluid-structure coupled system involving the transient Stokes equations (in a fixed domain) and a thin-walled solid elastodynamics model.This system is often used as model problem for the analysis of time-splitting schemes for incompressible fluid-structure interaction (see, e.g., [4,26,19,28,13,9]).Indeed, it retains the fundamental numerical difficulties that have to be faced also in general incompressible fluid-structure systems.A large amount of added-mass in the system is known to severely compromise stability and accuracy in standard explicit coupling schemes (i.e., those which invoke the fluid and solid solvers only once per time-step, see, e.g., [44,22,32,53]).The simplest approach to overcome these issues is to resort to a strongly coupled scheme (i.e., one in which the interface coupling is implicitly discretized in time), but at the expense of a higher computational complexity.
The development and the analysis of time splitting schemes which avoid strong coupling without compromising stability and accuracy has been a very active field of research during the last fifteen years.The vast majority of the studies have been devoted to the case of spatial approximations based on fluid meshes which are fitted to the interface (see, e.g., [27,49,5,4,26,12,13,36,46,6,34,14]).For many applications, such a mesh compatibility can however be cumbersome to maintain in practice (see, e.g., [48,33,51,8,19,42,43]).
The earliest explicit coupling schemes with unfitted meshes have been reported in [8,43], using the immersed boundary method, and in [19,42], using unfitted Nitsche approximations with overlapping meshes.Nevertheless, these methods suffer from major stability/accuracy issues which either require severe time-step restrictions (see [8,19]) or are limited by the amount of added-mass in the system (see [42,43]).A new class of semi-implicit schemes with unfitted meshes has been recently reported in [30,1] for the case of the coupling with thin-walled solids.These methods robustly avoid strong coupling but at the expense of introducing additional unknowns in the fluid sub-problem (intermediate solid velocity).Fully explicit variants of these approaches have been derived in [30] and in [10].Nevertheless, the formulation of the former in the case of immersed solids remains open and the accuracy of the latter relies on a grad-div penalty stabilization (for enhanced mass conservation) which spoils the conditioning of the fluid subsystem.
Additional difficulties, which are not considered in this paper, are the cases of curved and dynamic interfaces.Both subjects have not been particularly analyzed in the literature.We refer to [45] for stability and error analysis accounting also geometric errors.When dynamic interfaces are involved, a further complication derives from the variation of the approximation spaces in time.For the study of parabolic problems with moving interfaces we refer to [54,37].
This paper is devoted to the numerical analysis of the unfitted mesh semi-implicit coupling scheme recently introduced in [29].The method combines a Nitsche based unfitted mesh spatial approximation with a fractional-step time-marching in the fluid.The viscous part of the coupling is treated in an explicit fashion (which avoids strong coupling), while the remaining fluid pressure and solid contributions are treated implicitly (which guarantees added-mass free stability).Robust a priori error estimates are derived for two extrapolated variants (r = 1, 2 stands for the extrapolation of the solid velocity).The analysis highlights the fundamental role played by the time discretization of the Nitsche's penalty term in the stability and accuracy of the splitting.In particular, an O(τ r /h 1 2 ) splitting error is obtained instead of the standard O(τ /h) for the stabilized explicit coupling scheme of [19].The superior accuracy of the method is also supported by numerical experiments in an academic benchmark.
The rest of the paper is organized as follows.Section 2 presents the continuous setting.Its numerical approximation is discussed in Section 3. The numerical analysis of the semi-implicit scheme is reported in Section 4. Numerical evidence on the convergence properties of the methods is given in Section 5. Finally, Section 6 summarizes the main conclusions of this work.

Problem Setting
We consider a linear fluid-structure interaction problem in which the fluid is described by the Stokes equations in a fixed polyhedral bounded domain Ω f ⊂ R d , with d = 2, 3 and the structure by a linear thin-walled solid model with mid-surface given by Σ, also assumed polyhedral.Let the boundary of Ω f be partitioned as ∂Ω f = Σ ∪ Γ and denote the outward unit normal to ∂Ω f by n.In this framework, the considered coupled problem reads as follow: find the fluid velocity and pressure u : with the initial conditions u(0) = u0, d(0) = d0 and .d(0) = .d0.Here, the symbols ρ f and ρ s stand, respectively, for the fluid and solid densities.The thickness of the solid is denoted by ε and the fluid Cauchy stress tensor is given by where µ denotes the fluid dynamic viscosity and I is the identity matrix.The relations in (3) enforce, respectively, the kinematic and dynamic interface coupling conditions.The abstract differential operator L in (2) describes the elastic behaviour of the solid.
In the following, we will make use of the usual Sobolev's spaces H m (Ω)(m 0), with norm • m,Ω and seminorm | • |m,Ω, along with the closed spaces H 1 Γ (Ω), of functions in H 1 (Ω) with zero trace on Γ, and L 2 0 (Ω), of functions in L 2 (Ω) with zero mean in Ω.The scalar product in L 2 (Ω) is denoted by (•, •)Ω.For the weak formulation of the fluid problem we consider V = H 1 Γ (Ω f ) d and Q = L 2 (Ω f ) as the fluid velocity and pressure functional spaces, respectively.The standard Stokes bilinear forms are given by a(u, v) For the solid weak problem we suppose that W ⊂ H for all d ∈ D and w ∈ W .We define the following elastic energy norm on and we assume that the following continuity estimate holds for all w ∈ W , with β s > 0.

Numerical methods
In this section, we discuss three numerical methods for the approximations of the coupled problem ( 1)-( 3).These methods involve an unfitted mesh spatial discretization and different levels of fluid-solid time splitting.

Unfitted mesh spatial approximation
In a standard conforming finite element approximation, typically based on fitted meshes (see, e.g., [23,26]), the kinematic coupling condition (3)1 is strongly enforced.This is no longer feasible in the unfitted mesh setting.We consider the robust and optimal unfitted mesh method with overlapping meshes proposed in [19] (see also the related works [47,2]).Therein, the interfacial fluid-solid coupling is treated in a fully weak fashion via a Nitsche's type mortaring.Let be {T s h } 0<h<1 a family of triangulations of Σ, such that Σ = K∈T s h .We then consider the standard space of continuous piecewise affine functions associated to T s h , namely, For the approximation of the solid discrete space for the displacement and velocity we consider the following space We denote with {T f h } 0<h<1 a family of triangulations that cover the fluid domain Ω f such that: 1.Every T f h is fitted to Γ but, in general, not to Σ; 2. For every simplex K ∈ T f h , we have

In what follows, Ω f
h stands for the domain covered by T f h (i.e., the fluid computational domain).We denote by G h the set of elements of T f h intersected by Σ, by F h the set of the internal edges or faces of T f h , and by FG h the set of edges or faces of the elements of G h that do not lie on ∂Ω f h , namely, We consider the following space of continuous piecewise affine functions defined over For the approximation of the fluid velocity and pressure, we consider the following spaces Since the discrete pair V h /Q h lacks inf-sup stability, we consider the classical Brezzi-Pitkäranta symmetric pressure stabilization (see, e.g., [11]): Note that the pressure stabilization is defined over the computational fluid domain Ω f h .In order to guarantee robustness of the method with respect to the way the interface is cutting the fluid mesh, we consider the following ghost-penalty operator (see [15]): We can hence introduce the following total stabilization operator S h and associated semi-norm: so that the fluid discrete bi-linear form is given by Finally, the considered space semi-discrete unfitted mesh approximation of ( 1)-( 3) reads as follows: for t > 0, find u h (t), p h (t), .
Remark 3.1.It should be noted that in (10) the Nitsche mortaring is taken only from the fluid side (as in [39]).The fundamental reason for this is that there is no dissipative mechanism in the solid model (2)1,2 that can control the perturbation induced by the solid stresses on the interface.
Remark 3.2.In the numerical experiments of Section 5, the second assumption on T f h , that is, all the elements of Ω f h intersect the physical domain Ω f is relaxed.This is achieved by extending the ghost-penalty operator (8) to F h (all internal edges or faces of T f h ), i.e., This ensures the invertibility of the fluid stiffness matrix.It should be noted that the results of the numerical analysis reported in Section 4 below also hold in this case.

Time splitting schemes
This section is devoted to the time-discretization of (10).We first discuss the strongly coupled and the stabilized explicit coupling schemes reported in [19].Particular attention is paid to the well-known accuracy issues of the latter.We then discuss the semi-implicit projection based coupling scheme reported in [29], whose purpose was precisely to circumvent such difficulties without resorting to strong coupling.In the following, the parameter τ > 0 denotes the time-step length, ∂τ x n stands for the first-order backward difference formula and x ,r represents the r -th order explicit extrapolations to x n , viz.,

Strongly coupled scheme
Traditionally, the natural way to achieve numerical stability has been to consider a strongly coupled scheme, that is, a fully implicit time-discretization of (10).An example of such an approach is reported in Algorithm 1.The method is also known to deliver an optimal O(τ ) + O(h) accuracy in the energy norm (see [30]).The price to pay for this robustness is the resolution (at each time-step) of the hybrid coupled system (13), which can be computationally demanding in practice, particularly, due to its hybrid nature.Indeed, this monolithic system often yields ill-conditioned matrices (see, e.g., [2,52]) which require dedicated solvers.

Stabilized explicit coupling scheme
The stabilized explicit coupling scheme reported in Algorithm 2 enables a fully sequential decoupled time-marching of (10).Energy stability is achieved under a mild CFL-like condition (see [19]), thanks to the specific explicit treatment of the Nitsche penalty interface term and to the addition of an interface pressure stabilization in time (weakly consistent interfacial compressibility).The stability of the method is independent of the added-mass effect and of the considered local fluid and solid time-marching schemes.These features come however at a price: the sub-optimality of the error with respect to Algorithm 1.The interface time splitting (the explicit treatment of the penalty term in the right-hand side of ( 14)) introduces a truncation error whose leading contribution scales as O(τ /h), which prevents convergence when τ = O(h) (see Remark 4.3 and [18]).The pressure stabilization is introduced to control the pressure fluctuations at the interface, produced by the splitting in time of pressure and it has a contribution of the order O(τ ). Correction iterations are thus needed to enhance accuracy, under restrictive constraints on the discretization parameters.
• Fluid sub-step: find Roughly speaking, the lack of spatial uniformity in the splitting error of Algorithm 2 can be explained as follows.After spatial refinement, i.e., whenever h → 0, the solid sub-problem (14)  .In summary, the spatial discretization forces u n h − u n−1 h 0,Σ to be small as h → 0, by amplifying the time-splitting error.This is an essential ingredient of the scheme that guarantees numerical stability but it degrades accuracy.

Projection based semi-implicit coupling scheme
Algorithm 3 reports the projection based semi-implicit scheme of [29].The fundamental idea of this method, borrowed from [27] in the case of fitted mesh approximations, consists in combining a fractionalstep time-marching in the fluid (1) (see, e.g., [35]) with a semi-implicit treatment of the interface coupling (3).In Algorithm 3, the fluid is discretized in time with an incremental pressure-correction method and a backward-Euler method is considered for the solid.We introduce the fluid discrete viscous bi-linear form Note that the fluid viscous-step ( 16) is explicitly coupled with the solid, hence avoiding strong coupling (i.e., reducing computational complexity), whereas the coupled pressure-displacement system (17) guarantees added-mass free stability through the implicit treatment of the fluid incompressibility and solid inertial effects.For r = 2, the Algorithm 3 can be initialized with one step of the scheme with r = 1.
For n ≥ r: for all q h , w h ∈ Q h × W h .

Intermediate velocity step: find
for all v h ∈ V h .
It is worth noting that the discrete interface stresses in the (17)2 involve the same penalty term as in (16).In other words, the viscous stresses in (17)2 correspond to the variationally consistent residual of ( 16).This constitutes a fundamental difference with respect to Algorithm 2 (and also with respect to [3] with fitted meshes).
The next section provides an error estimate for Algorithm 3 which shows superior accuracy with respect to Algorithm 2, namely: O(τ r /h 1 2 ), with r = 1, 2, instead of O(τ /h).Furthermore, the numerical evidence reported in Section 5 suggests Algorithm 3 delivers practically the same accuracy as Algorithm 1, which is uniform with respect to h.The price to pay for this superior accuracy is threefold: • An additional CFL-like condition for stability (see [29,Theorem 1] and Theorem 4.1 below), with respect to Algorithm 1 ; • The solution of the coupled pressure-displacement system (17), not required in Algorithm 2; • A limited flexibility in the choice of the time-stepping for the fluid and solid sub-systems.
An additional appealing characteristic of Algorithm 3, is the reduced computational cost with respect to Algorithm 1.We refer to [29] for a run time study in the context of Nitsche-XFEM for immersed thin-walled structures.
Remark 3.3.In practice, we can avoid solving the third step, by inserting (18) into the viscous step (16), obtaining With the purpose of further investigating the space uniformity of the splitting error, we introduce a variant of the projection-based semi-implicit scheme with second-order backward difference discretization.The resulting method is shown in Algorithm 4 by choosing l = 2 and with the following definition of the l-order backward difference formulas: Observe that, for l = 1 we retrieve Algorithm 3 from [29].
In the following, Algorithm 4 is introduced only for comparison pourposes conducted in Section 5. We do not convey any stability analysis or error estimates with second-order in time discretization.
Algorithm 4 Projection-based semi-implicit scheme with l -th order backward difference.
For n ≥ r:

Numerical analysis
This section is devoted to the numerical analysis of Algorithm 3. We first recall the main ingredients for the energy stability of the method reported in [29] and extend the proof to cover the case of a second-order extrapolation (r = 2).An a priori error estimate is derived in Section 4.2.

Energy stability
We assume (see, e.g., [38,20]) that the following trace inequality holds for all v ∈ H 1 (K), K ∈ T f h and CT depending only on Σ.By combining (19) with a discrete inverse inequality, it follows for all v h ∈ V h .This estimates are fundamental for the energy stability of the method.Let define the discrete total energy E n h by the following expression: The following result states the conditional energy based stability of the approximation provided by Algorithm 3.
n≥1 be given by Algorithm 3 with r = 1, 2. Under the following conditions with α > 0, the discrete energy estimate presented below holds: for all n ≥ 1.As a result, Algorithm 3 is conditionally stable in the energy norm (22).
Proof.The proof for r = 1 was reported in [29].We recall here some of the steps and provide some details that include also the case r = 2, which follows similarly.We proceed by testing ( 16)- (18) with . By proceeding like in [29] (the sole difference lies in the choice of .d ,r h ), we get the following energy estimate for n ≥ r.Term T1 can be bounded by adding and subtracting .d n h , using the Cauchy−Schwarz, Young's and trace inequalities (20), as follows: with α1, α2 > 0. Similarly, for the second term, we have with α3 > 0. By inserting ( 27) and ( 28) into (26) we get for n ≥ r.We now set so that (29) yields for n ≥ r.We now replace the upper index n by m and sum over m = 2 . . .n and multiply by τ .This yields for n ≥ r.By rearranging the terms in the first sums, we get for n ≥ r.
In the case r = 1, the previous bound yields the energy estimate provided in [29].For r = 2, we need to control the contributions coming form the initialization step, namely, which are obtained from one step of Algorithm 3 with r = 1.We hence consider (30) with (n = 1, r = 1), which yields Hence, by adding this expression to (30), we finally get for n ≥ 1.The energy estimate (22) hence follows from (31) under the conditions ( 23)-( 24), which completes the proof.
Remark 4.1.A similar stability analysis can be derived in the case of a thick-walled solid.The solid quantities appearing on the interface Σ, such as .
, are controlled on the whole solid domain using element-wise trace inequalities.This yields to a parabolic CFL-type stability condition, namely, which is more restrictive than in the case of a thin-walled solid.An analogous stability result is reported in [27] for the non-incremental version of Algorithm 3 within the framework of fitted mesh.Stability is guaranteed under the CFL-like condition ρ f h 2 + 2µτ ρ s εh for a thin-walled solid and ρ f h 2 + 2µτ ρ s h 2 in the case of thick-walled solid.

A priori error estimate
In the following we use the notation v n def = v (tn) (without subscript h) to denote the temporal value v (tn) of a time dependent function v.For conciseness, an abuse of notation will be committed by denoting (∂tv) n with ∂tv n .Furthermore, the symbol indicates inequalities up to a multiplicative constant (independent of the discretization parameter h and of the physical parameters).We consider the following mesh-dependent seminorms for functions defined on the interface Σ: where ΣK denotes the part of the interface intersecting the simplex K, i.e., ΣK def = K ∩ Σ.For the sake of simplicity, in the error analysis we assume that the interface Σ is flat.Furthermore, the elements of the solid mesh are supposed to be grouped in disjoint macropatches Pi, with meas(Pi) = O(h d ).Each (d − 1)-dimensional macro patch Pi is assumed to contain at least one interior node and the union of the Pi is assumed to cover Σ, viz., ∪iFi = Σ.
The discrete interpolation operators are those introduced in [19] (see also [30]).For the solid displacement, we consider the elastic Ritz-projection operator π s h : W → W h defined by the relation for all w h ∈ W h , and for which there holds for all w ∈ H 2 (Σ) d ∩ W .For the solid velocity, we consider the operator I h : W → W h which is defined as a correction of the operator π s h by the relation with αi ∈ R to be fixed with a constraint detailed below.The ϕi are functions with support in the macropatches Pi, such that 0 ϕi 1, ϕi 0,P i h and take the value 1, component-wise, in the interior nodes of the associated patch Pi.The scalars αi are chosen so that the following condition holds: This orthogonality condition is used in the error analysis to control the interface terms coupling the fluid pressure and the solid velocity.We refer to [7] for the detailed construction of such an operator.It can be shown (see [19,Lemma 3.3]) that for all w ∈ H 2 (Σ) Since the physical solution and the discrete one, are defined on different domains, namely Ω f and Ω f h , with Ω f ⊂ Ω f h , we assume that there exist two linear continuous lifting operators E2 : , [25,50]).To interpolate the resulting extended fluid solution we consider the Scott−Zhang operator isz, see [24] for extra details.Then it holds (see [19,Lemma 3.3]): for all v ∈ H 2 (Ω f ) d , q ∈ H 1 (Ω f ) and w ∈ H 2 (Σ) d ∩ W .Moreover, using an inverse inequality, (36) and the stability of the extension operator we have the following stability result for the gradient projection For the pressure and ghost-penalty stabilization operators ( 7)-( 8), the following consistency properties hold (see, e.g., [19,21]): and g h (iszE2v, iszE2v) In the following we will make use of the discrete Grönwall lemma (see, e.g., [41]), which we collect here without a proof.For the a priori error estimate, we assume that the exact solution of problem ( 1)-( 3) has the following regularity, for a given final time T τ : We define the energy norm of the error at time tn as .
We can then state the following a priori error estimate.
for r = 1, 2. Suppose that the exact solution has the regularity (40) and that the stability conditions (23)-( 24) hold.Then, for n ≥ r and nτ < T , we have the following discrete error estimate: where {ci} 3 i=1 denote positive constants independent of h and τ , but which depend on the physical parameters and on the regularity of the exact solution.
Proof.The proof combines some of the arguments reported in [16,19].Note however that analysis of [19] focuses on the spatial semi-discrete problem (10) and the work of [16] is limited to a pure fluid problem.Multiplying (1)-( 2) by (v h , q h ) ∈ V h × Q h and w h ∈ W h , integrating by parts over Ω f and using (1)3 and (3)2 we obtain 1. Fluid sub-problem: 2. Solid sub-problem: Note that only the viscous term has been integrated by parts in the fluid.
On the other hand, owing to the kinematic coupling condition (3)1, we also have 1.Fluid sub-problem: 2. Solid sub-problem: for all w h ∈ W h .
Thereafter, using the lifting operators (component-wise) we introduce the following decomposition of the errors for the fluid: and for the solid: By adding and subtracting ∂τ π s h d n , we can rewrite ξn We also introduce the following notations: In particular, owing to (48), we have Similar, from (49), one straightforwardly gets the following useful relations: The essential part of the proof focuses on deriving an a priori estimate for the discrete errors in terms of the following energy norm .
To this purpose, we first focus on the fluid subsystem.By subtracting ( 16) from the momentum equation of ( 43) at t = tn, with n ≥ r, we get Owing to the error decompositions ( 45) -( 46) and using ( 50)-( 51), the identity ( 52) can be rewritten as for n ≥ r.
For the pressure, subtracting the pressure-projection step of ( 17) from the mass conservation equation ( 43) at t = tn, with n ≥ r, we get the following relation Again, using the definition of error decomposition ( 45)-( 46), the coupling kinematic condition (3)1 and (51), from ( 54) we obtain for n ≥ r.
Finally, adding and subtracting iszE2u n , iszE1p n , iszE1p n−1 in (18) and using (51), we obtain the following relation for the end-of-step velocity error for n ≥ r.
Subtracting the solid problem of ( 17) from ( 44) t = tn, with n ≥ r, and using the relation (50), we obtain Thus, using ( 45) -( 46) and ( 50), we finally get the equation for the solid discrete errors: for n ≥ r.Note that term a s ξ n π , w h vanishes due to the definition of the solid velocity projection operator (32).
Owing to ( 53), ( 55), ( 56) and (57), we have that the the discrete errors (θ n h , θ n h , y n h , ξ n h , ξn h ) satisfy a time-stepping scheme similar to Algorithm 3, but with a modified right-hand side and pressure increment (i.e., we have ).Therefore, we can leverage the stability arguments of Theorem 4.1 to derive an a priori error estimate.We proceed by testing ( 53), (55), ( 56) and ( 57) with By adding the resulting expressions, using the steps of Theorem 4.1 under condition (23), we obtain the following energy inequality for the discrete errors: for n ≥ r, and with the right-hand side G ,r h defined by Considering condition (23), equation ( 58) can be written as: The lack of telescoping sum on the pressure terms ∇y n h is not an issue (see, e.g., [16]).Indeed, using (37) we have so that by inserting this expression into (58), we have for n ≥ r.We now replace the upper index n by m and sum over m = r . . .n, this yields for n ≥ r.
We proceed by estimating G ,r h , by treating each term in (59) separately.The first term can be bound in a standard fashion using a Taylor expansion, (36), the Cauchy-Schwarz and the Poincaré's inequalities with constant CP.This yields with ε1 > 0. Observe that the last term can be absorbed in the left-hand side of (60) with ε1 small enough.For term T9 we proceed in a similar fashion.Using (36), we get with ε9 > 0 and where the last term can can be controlled in (60) using a Grönwall argument (Lemma 4.2).For term T2, using (36), we have The second term can be absorbed in the left-hand side of (60) for ε2 > 0 small enough.
The last term can be, once again, absorbed in the left-hand side of (60), for ε3 > 0 sufficiently small.Terms T4 and T5 involve the Nitsche splitting error, namely θ n π − χ ,r π 1 2 ,h,Σ .Using ( 35), ( 36) and a Taylor expansion we have To estimate T4, we follow the same idea of [19].Using the robust trace inequality (21) combined with (65), we get Once more, the last term can be absorbed in the left-hand side of (60), for ε4 > 0 sufficiently small.Similarly, for T5 we have Note that the last term can be included in the left-hand side of (60) for ε5 > 0 small enough.Term T6 can be handled using (36) as follows: Again, the last term can be absorbed in the left-hand side of (60), for ε6 > 0 small enough.
Integrating by parts in T7, we have Term T7,1 can be easily handled by using (36) as follows: The second term can be absorbed in the left-hand side of (60), for ε7 > 0 sufficiently small.For the second term of (69), we proceed as in [19] (see also [30]).We denote by ȳn i ∈ R the average of y n h over the interface patch Pi.Combining the trace inequality (19) with the orthogonality property (34) of the interpolation operator I h and the standard estimate y n h − ȳn i 0,P i h ∇y n h 0,P i , term T7,2 can be estimated as follow: It should be noted here we have assumed that the solid mesh step has an asymptotic regime similar to the fluid mesh step, namely, h s = O(h f ).As for T7,1, the last term in (71) can be, once again, absorbed in the left-hand side of (60), for ε7 > 0 sufficiently small.For term T8, using the weak consistency of the stabilization operators (38) and (39), we have Again, the last term can be absorbed in the left-hand side of (60), for ε8 > 0 small enough.Term T10 can be bounded using the continuity estimate for the elastic bilinear form ( 5), ( 32), (36) and a triangular inequality.This yields Note that the first term can controlled via Lemma 4.2 in (60).In summary, the term n m=r G m,r h in the right-hand side of (60) can be estimated by collecting the estimates (61)-(73) and by inserting them into (59), for n ≥ r.The desired estimate (41) hence follows from (60) together with the stability condition (24) and Lemma 4.2 with and by noting that, owing to the initial data, we have for r = 1, 2. This completes proof.
Corollary 4.1.Assume that Algorithm 3 with r = 2 is initialized with one step of the method with r = 1.Then, under the assumptions of Theorem 4.3, for n ≥ 1 and nτ < T , the following discrete error estimate holds for the scheme with r = 2: where {ci} 4 i=1 denote positive constants independent of h and τ , but which depend on the physical parameters and on the regularity of the exact solution.
Proof.For r = 2, we have to bound the contributions from the initialization step in the right-hand side of (60), namely, To this purpose, we use the fact that the initialization of Algorithm 3 with r = 2 is provided by the first step of the scheme with r = 1.We can hence use the estimate provided by (60) with r = 1 and n = 1 to control (77).More precisely, using (75), this yields Hence, by inserting this estimate in (60), we get for n ≥ 1. Owing to the initialization procedure, the bounds provided in (66)-(67) for terms T4 and T5 of G 1,1 h yield a O(τ 3 /h) splitting error, by noting that The estimate (41) for r = 2 hence follows from (78) together with the stability condition (24) and Lemma 4.2 with (74).This completes proof.
We conclude this section with a series of remarks.
Remark 4.2.For r = 2, the last term in (41) comes from the bound of the first step of Algorithm 3 with r = 1, that is, the estimate given by (79).This bound is quasi-optimal in time because the Taylor expansions are evaluated in L 2 (0, T ) instead of L 1 (0, T ).Alternatively, one could avoid this term by initializing Algorithm 3 with the first-step of Algorithm 1.
Remark 4.3.Note that Algorithm 2 introduces the following perturbations terms in the discrete error equation with w h = τ ξn h .The first term leads to the following bound: The second term can be controlled via Lemma 4.2 while the first yields the above mentioned O(τ /h) sub-optimal splitting error.
Remark 4.4.As shown in Theorem 4.3, the discrete error estimates of Algorithm 3 contains terms of order O(τ r /h 1/2 ), which are not visible numerically (see Section 5).To fully understand the impact of selecting the same penalty term in the viscous step as in the solid sub-step in Algorithm 3, we consider the coupling of a parabolic equation with and an hyperbolic one.The considered coupled problem reads as follow: find u : with the respective initial conditions.We propose to discretize the problem via a loosely coupled scheme, inspired by the semi-implicit scheme of Algorithm 3. The fully discrete approximation results in the following (explicit) scheme: For n ≥ 1: When considering loosely coupled schemes with Nitsche's coupling, the sub-optimal terms come typically from the fact that we introduce a time-splitting error inside the Nitsche's penalty term, which is scaled with an h −1 .A possible way to overcome this issue, is to remove the time-splitting error from the Nitsche's penalty term, by introducing an error in time within the definition of the projection errors.Thus, considering the following decomposition of the errors for the parabolic-hyperbolic explicit scheme: it can be proven that the scheme delivers optimal space and time accuracy.More in detail, using similar arguments of the proof of Theorem 4.3, we will get the following terms inside the Nitsche's penalty part: which does not contain error in time, in fact the arising terms involving θ n h − ξn−1 h are controlled via the stability result and terms involving θ n π − ξn π have optimal convergence order.The only terms which contain ξn τ are the corresponding T9 and T10 terms of (59) and their optimality can be proved.A similar strategy fails when considered for the semi-implicit scheme of Algorithm 3. In particular we will retrieve terms of order O(τ /h 1/2 ) when controlling the pressure term T7,2 of (59).

Numerical experiments
In this section, we illustrate via numerical experiments the convergence properties of Algorithm 3 with r = 1, 2 (semi-implicit scheme) in an academic numerical example.The obtained results are compared with those of Algorithms 1 (strongly coupled scheme) and Algorithm 2 (stabilized explicit coupling scheme).The considered test case is the well-known fluid-structure interaction benchmark describing the propagation of a pressure wave within a straight two-dimensional elastic tube (see, e.g., [31,17,19]).In the following, all the units are given in the CGS system.
The fluid domain is defined as As expected, all the considered considered methods deliver a numerical solution with a stable pressurewave propagation.For illustration purposes, Figure 4 provides the snapshots of the fluid pressure and solid deformation at time t = 5 × 10 −3 , 10 −2 and 1.5 × 10 −2 , obtained with τ = 10 −4 and h = 0.05 using respectively Algorithms 1-3.The solid displacement has been amplified by a factor 20. A very good agreement between Algorithm 1 and Algorithm 3 (r = 1, 2) is clearly visible, while a difference on the solid displacement is noticeable with Algorithm 2.
In order to quantify the accuracy properties of each coupling scheme we have evaluated the convergence histories by uniformly refining in space and in time (h, τ ) ∈ 0.1/2 i , 2 × 10 −4 /2 i 4 i=0 . (81) Figure 5 shows the corresponding solid displacement at t = 1.5 × 10 −2 for i = 0, .., 3 and the different coupling schemes.As in Figure 4, a very good fit is observed between Algorithm 1 and Algorithm 3 (r = 1, 2), while a degradation of accuracy is visible for Algorithm 2 under space-time refinement.The  depicted reference solution has been generated using the strongly coupled fitted method with a high space-time grid resolution (h = 3.125 × 10 −3 and τ = 10 −6 ).Aiming to investigate the space uniformity of the splitting error, we compare Algorithms 1−3 with the variant of the projection-based semi-implicit scheme with second-order backward difference discretization, shown in Algorithm 4. Figure 6 reports the convergence history of the solid displacement at time t = 1.5 × 10 −2 , in the relative elastic energy-norm.Note that, from the choice of space and time discretization parameters in (81), we have τ = O(h).The results show that Algorithm 3 with r = 1, 2 and Algorithm 1 retrieve the overall optimal first-order accuracy O(h) of Algorithm 1.As expected, Algorithm 2 shows a non-convergence behaviour.This points out the expected /h) sub-optimality of the splitting error (see Remark 4.3).Algorithm 4 exhibits better accuracy than Algorithm 3. Full second-order accuracy is not achieved due to the first-order extrapolation of the pressure.Finally, it is worth noting that no effect from the O(τ /h

Conclusion
In this paper, an a priori error analysis of the unfitted mesh based semi-implicit coupling scheme introduced in [29] has been performed in the context of linear-fluid structure interaction system involving a thin-walled solid.The considered method combines a Nitsche based unfitted mesh spatial approximation with a fractional-step time-splitting in the fluid.Strong coupling is avoided by treating explicitly the  The stability analysis of [29] has been extended to the case of second order extrapolation of the solid velocity (r = 2) under a similar hyberbolic CFL condition.The proposed error analysis predicts superior accuracy of the semi-implicit scheme of [29] with respect to the stabilized explicit scheme of [19], namely O(τ r /h 1 2 ), r = 1, 2, instead of O(τ /h).Extensive numerical evidence indicated that the semi-implicit algorithm and the strongly coupled (from [19]) deliver practically the same accuracy.In particular, no effect from the non-uniformity in space is visible, for both the first and second-order backward difference variants.

Figure 4 :
Figure 4: Snapshots of the fluid pressure and deformation (magnified) at different time instants.