A PARTITIONED NUMERICAL SCHEME FOR FLUID–STRUCTURE INTERACTION WITH SLIP

Abstract. We present a loosely coupled, partitioned scheme for solving fluid–structure interaction (FSI) problems with the Navier slip boundary condition. The fluid flow is modeled by the Navier– Stokes equations for an incompressible, viscous fluid, interacting with a thin elastic structure modeled by the membrane or Koiter shell type equations. The fluid and structure are coupled via two sets of coupling conditions: a dynamic coupling condition describing balance of forces, and a kinematic coupling condition describing fluid slipping tangentially to the moving fluid–structure interface, with no penetration in the normal direction. Problems of this type arise in, e.g., FSI with hydrophobic structures or surfaces treated with a no-stick coating, and in biologic FSI involving rough surfaces of elastic tissues or tissue scaffolds. We propose a novel, efficient partitioned scheme where the fluid subproblem is solved separately from the structure sub-problem, and there is no need for sub-iterations at every time step to achieve stability, convergence, and its first-order accuracy. We derive energy estimates, which prove that the proposed scheme is unconditionally stable for the corresponding linear problem. Moreover, we present convergence analysis and show that under a time-step condition, the method is first-order accurate in time and optimally convergent in space for a Finite Element Methodbased spatial discretization. The theoretical rates of convergence in time are confirmed numerically on an example with an explicit solution using the method of manufactured solutions, and on a benchmark problem describing propagation of a pressure pulse in a two-dimensional channel. The effects of the slip rate and fluid viscosity on the FSI solution are numerically investigated in two additional examples: a 2D cylindrical FSI example for which an exact Navier slip Poiseuille-type solution is found and used for comparison, and a squeezed ketchup bottle example with gravity enhanced flow. We show that the Navier-slip boundary condition increases the outflow mass flow rate by 21% for a bottle angled at 45 degrees pointing downward, in the direction of gravity.

used in FSI problems with no slip, the tangential components are split differently, without respecting the so called "added mass effect". Using analysis based on energy estimates, we show that the proposed scheme is still unconditionally stable. This is interesting because it shows that added mass effect is dominant in the normal direction, and that the tangential forces have secondary importance in the added mass effect since they only "reposition" the points on the fluid-structure interface, without significantly changing the displaced fluid mass. This is new. It is because of friction due to slip that our scheme does not require implicit coupling between the fluid and structure inertia in both the normal and tangential direction, and still provides an unconditionally stable partitioned scheme without the need for sub-iterations at every time step.
In this manuscript we also present convergence analysis based on a priori error estimates, which shows that under a time-step condition, the proposed numerical method is first-order accurate in time and optimally convergent in space. Numerical examples are presented to verify the convergence properties and to investigate the effects of the slip rate and fluid viscosity on the solution of the coupled problem. The theoretical rates of convergence in time are confirmed numerically on an example with an explicit solution using the method of manufactured solutions, and on a benchmark problem describing the propagation of a pressure pulse in a two-dimensional channel. The effects of the slip rate and fluid viscosity on the FSI solution are numerically investigated in two additional examples: a 2D cylindrical FSI example for which an exact Navier slip Poiseuilletype solution is found and used for comparison, and a squeezed ketchup bottle example with gravity enhanced flow. We show that the Navier-slip boundary condition increases the outflow mass flow rate by 21% for a bottle angled at 45 • pointing downward, in the direction of gravity. This paper is organized as follows. The FSI model with slip is introduced in Section 2. The proposed numerical scheme is presented in Section 3. Stability analysis is preformed in Section 4, and convergence analysis is presented in Section 5. Numerical examples are presented in Section 6. Finally, conclusions are drawn in Section 7.

Problem definition
Let Ω ⊂ R d , d = 2, 3 denote the reference fluid domain and assume ∂Ω = Γ ∪ Γ in ∪ Γ out , where Γ represents the elastic part of the boundary while Γ in and Γ out represent artificial inflow and outflow sections, see Figure 1. We assume that Ω is filled with an incompressible, viscous fluid, whose dynamics is modeled using the Navier-Stokes equations for an incompressible, viscous fluid defined on a time-dependent domain Ω(t), which is not known a priori : ρ F (∂ t u + u · ∇u) = ∇ · σ(u, p), ∇ · u = 0, in Ω(t), t ∈ (0, T ), (2.1) where ρ F denotes fluid density, u is the fluid velocity, σ = −pI + 2µ F D(u) is the fluid Cauchy stress tensor, p is the fluid pressure, µ F is the kinematic viscosity coefficient, and D(u) = 1 2 (∇u + ∇ τ u) is the fluid strain rate tensor.
The structure elastodynamics is described by a linearly elastic model, given by a reduced membrane or shell model in general form: ρ S ω∂ tt η + L e η = f , on Γ, t ∈ (0, T ), (2.2) where η is a vector function describing the structure displacement from its reference configuration Γ, L e is a coercive and continuous differential operator on the function space specified later, associated with the elastic energy of the thin structure. For example, for a linearly elastic membrane, the functions space is defined by the Sobolev space H 1 , while for a Koiter shell with radial displacement, the functions space is defined by H 2 . Different models will be used in this manuscript. The notation f denotes force density (load), and ρ S and ω are the structure density and thickness, respectively. The fluid and structure are coupled via two sets of coupling conditions: the dynamic and kinematic coupling conditions. To state the coupling conditions we introduce ξ to denote the structure velocity, ξ = ∂ t η. The coupling conditions are given by: • Dynamic coupling condition: (2.3) stating that the thin structure elastodynamics is driven by the fluid normal stress at the interface. The term J is the Jacobian of the transformation between the Eulerian and Lagrangian formulations of the fluid and structure problems, respectively. • Kinematic coupling condition: ξ · n = u · n| Γ(t) on Γ × (0, T ) (the non-penetration condition) (ξ − u| Γ(t) ) · τ i = αJσn · τ i | Γ(t) , i = 1, . . . , d − 1 on Γ × (0, T ) (slip between fluid and structure) Here, n and τ denote the outward normal and tangential unit vectors to the fluid domain, and α is the slip rate.
At the inlet we prescribe Dirichlet or Neumann boundary conditions, depending on the problem at hand, and at the outlet we prescribe Neumann boundary condition: u = u in (t), or σn = −p in (t)n on Γ in × (0, T ), (2.4) σn = −p out (t)n on Γ out × (0, T ), (2.5) and assume that the structure is clamped at the edges. For membranes, this means η = 0 on ∂Γ.

The ALE framework
To deal with the difficulty associated with the fact that the fluid domain changes in time, we adopt the arbitrary Lagranian Eulerian (ALE) framework [33,50]. The ALE approach is based on introducing a family of arbitrary, smooth, homeomorphic mappings A t defined on the reference domain Ω such that, for each t ∈ (t 0 , T ), A t maps the reference domain Ω into the current domain Ω(t): A t : Ω → Ω(t) ⊂ R n , n = 2, 3, x = A t (x) ∈ Ω(t), forx ∈ Ω.
In particular, we consider the ALE mappings A t defined via the harmonic extensions of the boundary data onto the entire fluid domain, i.e., the mappings defined by the following boundary value problems on the reference domain Ω: in Ω, (2.6) Written in the ALE framework, system (2.1) reads as follows: Find u and p, withû( denotes the domain velocity. The fluid domain can be computed as Ω(t) = (I + A t )Ω. Note that ∂f ∂t x denotes the time derivative of f evaluated with respect to the reference domain.

Energy estimate
To state the energy estimate, we introduce the following notation for the projection operator onto the tangent space of the reference configuration Γ [26]: (2.12) Lemma 2.1. Consider the FSI problem modeled by the Navier-Stokes equations in ALE form (2.9)-(2.10) and the structure equation (2.2), with the flow driven by the dynamic pressure inlet and outlet boundary data: Then the following formal energy inequality holds for the coupled problem (2.1)-(2.5): Here, C depends on the initial and boundary data, · Le(Γ) is the norm associated with the coercive, linear continuous operator L e , and the constant c is associated with the coercivity of the structure operator L e .
Proof. We multiply (2.9) by u and (2.10) by p, and integrate over Ω(t). Furthermore, we multiply (2.2) by ∂ t η and integrate over Γ. Using the boundary conditions and adding the equations together, we obtain (2.16) We recast the first term to the reference configuration, use the Euler expansion formula, and recast back, yielding Integrating one half of the convective term by parts (see [63] for details), we obtain Decomposing the interface integral in (2.16) in the normal and tangential components, taking the coupling conditions into account and estimating the forcing term using usual inequalities (see e.g., [63]), we obtain the desired energy estimate.
The energy estimate in Lemma 2.1 shows that the proposed model is reasonable in the sense that the total kinetic and elastic energy of the coupled problem, plus dissipation due to fluid viscosity and slip friction, are all bounded by a constant depending only on the initial and boundary data.

The numerical scheme
The splitting in continuous form. We propose a partitioned scheme which separates the coupled problem into a fluid sub-problem and a structure sub-problem, which are solved only once at every time step, avoiding expensive sub-iterations typically associated with Dirichlet-Neumann partitioned FSI schemes [3-5, 43, 54, 77]. The partitioning proposed here is different from the classical partitioned schemes since we deal with the Navier slip boundary condition, and the coupling between the fluid and structure in the tangential direction to the interface needs to be carefully split in order to keep the accuracy of the scheme first order, and provide a stable loosely coupled scheme.
We start by semi-discretizing the problem in time. Let t n = n∆t for n = 1, . . . , N where T = N ∆t is the final time. We will semi-discretize the coupled evolution problem in time using the Lie operator splitting scheme [44] applied to our coupled problem written as a first-order system in time: Here, A 1 and A 2 correspond to a fluid and structure operators/sub-problems. Following the Lie splitting strategy, this problem is semi-discretized in time in such a way that at every sub-interval (t n , t n+1 ) the sub-problem determined by the operator A 1 is solved (dU/dt = A 1 U ) with the initial data given by the solution from the previous time step, i.e., the solution at t = t n , and then the sub-problem determined by the operator A 2 is solved (dU/dt = A 2 U ) over (t n , t n+1 ) with the initial data given by the just calculated solution of the sub-problem A 1 . This approach has been used to split FSI problems with the no-slip condition in [12,17,59]. Here, we generalize this approach to deal with the Navier slip condition. We propose to split the normal component of the dynamic coupling condition (2.3) into two parts by using the Lie splitting strategy: one part will serve as a boundary condition for the fluid sub-problem, and the other will be a part of the structure sub-problem, such that both parts include the time-derivative term (as defined by the Lie splitting): Additionally, in the fluid sub-problem the normal component of the kinematic coupling condition (nopenetration) will be used to replace the structure velocity by the normal trace of the fluid velocity on Γ(t). Thus, so far we have used the normal component of the dynamic coupling condition and the normal component of the kinematic coupling condition. The remaining tangential components of the dynamic and kinematic coupling conditions are split between the fluid and structure sub-problems as follows. First, the tangential component of the dynamic coupling condition is rewritten by using the tangential component of the kinematic coupling condition (slip) so that: The portion: is then used in the structure sub-problem, and the portion: in the fluid sub-problem. Thus, the structure sub-problem is defined by: ρ S ω∂ t ξ · n + L e η · n = −Jσn · n, defined on Γ, and the coupling conditions that enter the fluid sub-problem as boundary conditions for the Navier-Stokes equations are (from (3.1) and (3.2)): ρ S ω∂ t u · n + Jσn · n = Jσn · n, The splitting scheme in semi-discrete form. To write the splitting in semi-discrete form, we introduce the following notation for the approximate time derivative: By using the Backward Euler scheme, the time-discrete numerical scheme is given by the following.
Step 1 (Structure). Given u n and p n from the previous time step, calculate ξ n+1 and η n+1 = ∆tξ n+1 + η n such that Step 2 (Fluid). Given the thin structure velocity ξ n+1 , calculate u n+1 and p n+1 such that: Step 3 (Fluid domain update). Given the displacement η n+1 of the boundary Γ, we update the fluid domain Ω(t n+1 ) in a classical way by using the harmonic extension Ext(η n+1 ) of the boundary data η n+1 onto the entire domain, and compute the ALE velocity w n+1 which is defined as the time derivative of the ALE mapping . More precisely, we calculate the ALE mapping as: and update ∈ Ω(t n+1 ) and x n = A −1 t n (x) ∈ Ω(t n ), forx ∈ Ω.
Remark 1. We remark that in this scheme the kinematic coupling condition describing continuity of the normal components of the velocity between the fluid and thin structure is satisfied asynchronously, and not identically. It is, in general, not true in this scheme that u n · n F = ξ n · n F . Only in the limit as ∆t → 0, this will be satisfied. In fact, it can be shown that this condition is satisfied to the second-order accuracy in ∆t.

Remark 2.
We further remark that the splitting proposed here is slightly different from the splitting discussed in the existence proof in [63], and it is significantly different from the splitting schemes proposed by the authors in [13,60] to solve FSI problems with the no-slip kinematic coupling condition. One important difference is the form of the Robin boundary condition for the fluid sub-problem. In contrast with the earlier works [10, 13-15, 60, 63], the Robin boundary condition (3.8) ties the fluid and structure inertia implicitly only in the normal component of the inertia. The lack of implicit coupling between the fluid and structure inertia in loosely coupled schemes for problems with no-slip condition typically leads to instabilities due to the added mass effect. We show below that this is not the case here. Our energy estimate for the semi-discretized problem, presented in Theorem 4.1, shows that the energy of the semi-discretized problem is bounded, uniformly in ∆t, indicating that this scheme is also unconditionally stable. In the tangential direction dissipation due to slip friction is sufficient to provide stability of our scheme. It is because of friction due to slip that our scheme does not require implicit coupling between the fluid and structure inertia in both the normal and tangential direction, and still provides an unconditionally stable partitioned scheme without the need for sub-iterations at every time step.

The fully discretized scheme in weak form
We discretize problem (3.5)-(3.9) in space using a finite element method approach. For this purpose we introduce the following spaces: for all t ∈ [0, T ]. These are the spaces associated with the fluid, pressure, and structure problem, respectively. We recall here that the thin structure operator L E , obtained from the elastic energy of the thin structure, is coercive and continuous on the space V S , which defines a bilinear form on V s : and the norm (3.14) The finite element spaces are then defined as the subspaces Step 2 (Fluid): Given the thin structure velocity ξ n+1

Stability analysis
To simplify stability and convergence analysis and focus on the issues related to the splitting strategy, we will assume that the fluid domain is fixed (i.e., geometric nonlinearities are neglected) and neglect the nonlinear fluid advection. These two issues are related since energy estimates for the problem with advection must include moving domain in order to control the cubic terms that arise in the corresponding energy [60,63]. These simplifying assumptions are commonly used in stability analysis of FSI partitioned schemes [11,37,38]. The unconditional stability of the proposed scheme is given by the following theorem: be a solution of (3.15)-(3.16) with negligible nonlinear fluid advection and fixed fluid domain. Then, for each ∆t > 0, there exist the constants C P , C K , C T > 0, independent of ∆t, such that the following a priori energy estimate holds: where E N is the sum of the kinetic energy of the fluid, the kinetic energy of the structure and the elastic energy of the structure: D N denotes dissipation due to fluid viscosity and due to slip between the fluid and structure: and H N and N N denote the terms due to numerical dissipation, defined by: where · E is defined in (3.14) by the elastic operator L E .
Proof. We replace the test functions χ h in (3.15), and (ϕ h , ψ h ) in (3.16), by the following: (3.16). After multiplying the two equations by ∆t, adding them together, and using the identity 2a(a − b) = a 2 − b 2 + (a − b) 2 we get: To express the right hand-side in terms of the L 2 -norms of the fluid and structure quantities we notice that (3.8) implies: Now, the integral I on the right hand-side of (4.2) becomes: After canceling equal terms with opposite signs, we obtain the following energy equality: Using the Cauchy-Schwarz and Young's inequalities with > 0, as well as the trace, Poincare and Korn inequalities, we obtain: Taking into account (4.5) and summing from 0, . . . , N − 1, we obtain the desired estimate.

Convergence analysis
To study convergence of our scheme, we compare the solutions of the fully discretized problem (3.15)-(3.16) obtained using our partitioned scheme, with the weak solution of the monolithic, continuous problem. Since the test functions of the partitioned scheme do not satisfy the no-penetration condition, we will formulate the continuous, monolithic weak form using those same test functions, and then compare the two solutions.
For spatial discretization we use the Lagrangian finite elements of polynomial degree k for all the variables, except for the fluid pressure, for which we use elements of degree s < k. We assume that our finite element spaces satisfy the usual approximation properties, and that the fluid velocity-pressure spaces satisfy the discrete inf-sup condition.
Assume that the continuous solution, in addition to belonging to the space satisfies the following additional regularity assumptions: We introduce the following time discrete norms: Note that they are equivalent to the continuous norms since we use piecewise constant approximations in time. Furthermore, the following inequality holds Here, and in the rest of this text, we use a ( )b to denote the inequality a ≤ (≥)Cb, where constant C is independent of the time step size ∆t, mesh size h and slip rate α.
Let I h be the Lagrangian interpolation operator onto V S h . Similarly as in [17,38], we introduce a Stokes-like projection operator (S h , P h ) : Projection operators S h and I h satisfy the following approximation properties (see [38], Thm. B.5 and [27]): Then, the finite element theory for Ritz projections [27] gives We will use these results to compare the solution of our fully discretized partitioned scheme with the continuous solution. For this purpose we now present the monolithic, continuous, weak formulation, using the test function in V S h × V F h × Q h , which do not necessarily satisfy the no-penetration condition. Monolithic, weak formulation: The error equation. Subtracting (3.15)-(3.16) from (5.17), we obtain the following error equation We will use the error equation (5.18) to obtain the L 2 -error estimates, namely the estimates of the quantities: We split the error of the method as a sum of the approximation error θ n+1 r and the truncation error δ n+1 r , for r ∈ {f, p, η, ξ} as follows: The main result of this section is stated in the following theorem.
Theorem 5.1. Let ∆t > 0 and h > 0. Consider the solution (ξ n h , η n h , u n h , p n h ) of (3.15)-(3.16) at t n = n∆t, with discrete initial data , and let (ξ n , η n , u n , p n ) be the exact solution of the continuous problem (5.17) at t n = n∆t, with the same initial data. Assume that the exact solution satisfies assumptions (5.1)-(5.5) and that Then, for each n ≤ N , the following estimate holds: Proof. We start by rearranging the error equation (5.18) and taking the property (5.15) of the Ritz projection operator into account, to get Thanks to (5.11), the pressure terms simplify as follows: We rewrite the term ∆t Γ L e (δ n+1 η ) · δ n+1 ξ as follows: . Hence, using property (5.15) of the Ritz projection operator, the Cauchy-Schwartz and Young's inequalities, we have To estimate the first term on the right hand side of (5.25), similarly as in [17], we note that (δ n+1 h ) · n. Furthermore, after adding and subtracting the continuous velocity and pressure on the right hand side of (3.8), we see that the following relation holds on Γ: After employing the identity (5.27), we get: = ∆t 2 ρ S ω Γ (σ(e n f , e n p )n · n) σ(e n f , e n p )n · n − σ(e n+1 f , e n+1 p )n · n dS T1 + ∆t 2 ρ S ω Γ (σ(e n f , e n p )n · n) σ(u n+1 − u n , p n+1 − p n )n · n dS. T2. (5.28) The term T 1 can be estimated as follows: f , e n+1 p )n · n 2 L 2 (Γ) + ∆t 2 2ρ S ω σ(e n f , e n p )n · n 2 f , e n+1 p )n · n − σ(e n f , e n p )n · n 2 L 2 (Γ) .

(5.29)
To estimate the last term in (5.29), we again use identity (5.27), (3.8) and Young's inequality as follows: f , e n+1 p )n · n − σ(e n f , e n p )n · n 2 Finally, we estimate T 2 using the Cauchy-Schwartz inequality and Young's inequality to obtain: We bound the remaining terms in (5.25) as follows. First, by using the Cauchy-Schwartz, Young's, Poincaré-Friedrichs, and Korn's inequalities, we obtain: .
(5.32) By rewriting T 3 and using the Cauchy-Schwartz, Poincaré-Friedrichs, Korn's and Young's inequalities, we get Similarly, we have Combining the estimates above with equation (5.25), summing from n = 0, . . . , N − 1 and taking into account the assumption on the initial data, we have σ(e n f , e n p )n · n 2 By using Lemmas A.1, A.2 and A.3, we get the following estimate of the approximation and consistency errors: σ(e n f , e n p )n · n 2 L 2 (Γ) The term ρ S ω∆t 4 in the above inequality can be estimated by adding and subtracting δ n F and using trace-inverse inequality [67], where C T I depends on the angles in the finite element mesh, as follows: We combine the estimates above with (5.35) and recall that the error between the exact and the discrete solution is the sum of the approximation error and the truncation error. Thus, using the triangle inequality, approximation properties (5.12)-(5.16) and the Gronwall lemma we obtain the estimate from Theorem 5.1.
We have shown that the loosely coupled, partitioned scheme (3.15), (3.16) is unconditionally, stable, firstorder accurate in time, and optimally accurate in space. No sub-iterations are needed between the fluid and structure sub-problems to achieve this accuracy and stability. In the next section we present examples which confirm the theoretical findings, and study the influence of the Navier slip boundary condition on the solution of FSI problem as compared with that with no-slip.

Numerical examples
To investigate the computational model and the numerical properties of the proposed partitioned scheme, we present four examples. In the first example, the method of manufactured solutions is used to verify the convergence properties of the scheme on a simplified modeled describing the interaction between a compressible fluid and a thin structure assuming that the domain is fixed. In the second example, we investigate the rates of convergence on a benchmark problem describing a moving domain FSI problem, where the fluid is described using the incompressible Navier-Stokes equations. In the third and fourth example, different flow conditions, fluid domains, and structure models are investigated, including a ketchup bottle model showing the flow difference in the bottle with and without a no-stick coating. All the computations have been performed with FreeFem++ [47].

Example 1
The model. The first numerical example is focused on the computational study of the convergence rates in time using the method of manufactured solutions. We consider a two-dimensional example where the fluid flow is described by the time-dependent Stokes equations: in Ω, (6.1) where g is volumetric force, and both s and g are determined from the manufactured solution.
We assume that the fluid domain is fixed, corresponding to the assumptions made in Section 4. The fluid domain is given by Ω = (0, 1) × (0, 0.5), with the fluid-structure interface corresponding to Γ = (0, 1) × {0.5}. We model the structure using the elastic Koiter shell model [13] defined on Γ, accounting for both tangential (horizontal) and transverse (vertical) displacements η x and η y , given by: The dynamic coupling condition implies: The parameters used in this example are ρ F = µ F = C 0 = C 1 = C 2 = C 3 = α = 1, ρ S = 0.5 and h = 0.1.
The exact solution. The exact solution is given by

4)
p = e t cos(πx) ρ S ω + C 0 π + C 1 π + C 2 π(µ F − C 2 ) ρ S ω + C 3 π 2 sin(πy) + 2µ F cos(πy) , (6.5) We note that while the exact solution satisfies the coupling conditions described in Section 2, the fluid defined by this solution is not incompressible. Nevertheless, this is a good example to test convergence rates in time, which are predominantly affected by the splitting of the coupling conditions at the fluid-structure interface.
Boundary conditions. Dirichlet boundary conditions are imposed on the rigid boundary ∂Ω \ Γ, while the kinematic and dynamic coupling conditions specified in Section 2 are imposed on Γ = (0, 1) × {0.5}.
Computation. Numerical simulations were run using the following set of discretization parameters: corresponding to four levels of mesh refinement. We used P 1 -bubble/P 1 elements to discretize the fluid problem and P 1 elements to discretize the solid problem. Simulations were run until the final time of T = 0.1s. Results. Convergence rates in time were studied by comparing the exact solution with the numerical results for different levels of mesh refinement. Figure 2 shows the relative errors for the fluid velocity (top left), structure velocity (top right), and structure displacement (bottom). We observe that while the fluid and structure velocities converge with the rate greater than 1, the rates of convergence of the solid displacement agree with the theoretical predictions.

Example 2
The model. In this example, we investigate the rates of convergence on a benchmark problem describing the interaction between an incompressible fluid and a thin, elastic structure. As opposed to the first example, we consider an incompressible fluid described by the Navier-Stokes equations in a moving domain, as specified in Section 2. Here, we do not have an explicit solution as in Example 1, so we compare the results for different discretizations with the solution obtained using the finest discretization.
The structure is described using the same model as in Example 1, given by (6.2)-(6.3), with where E denotes the Young's modulus of elasticity, and ν denotes the Poisson's ratio of the elastic structure material. The parameters used in this example are given in Table 1.
The inlet and outlet data. The inlet and outlet forcing terms, p in and p out , as defined in (2.4)-(2.5), are given by where t max = 0.03 s with maximum pressure p max = 1.333 × 10 4 dyne/cm 2 . The final time is T = 10 ms. This type of inlet and outlet data is common in testing numerical solvers for FSI problems.
Initial data. Initially, zero fluid and structure velocity, and zero structure displacement were prescribed.
Computation. We used P 1 -bubble/P 1 elements to discretize the fluid problem and P 1 elements to discretize the solid problem. Numerical simulations were run using the following set of discretization parameters:  corresponding to six levels of mesh refinement. The solution obtained using the finest set of parameters was taken as a reference solution.
Results. The pressure in the domain at t = 2, 4, 6 and 8 ms, as well as the magnitude of the structure displacement are shown in Figure 3. Figure 4 shows the relative errors for the fluid velocity (top left), structure velocity (top right), and structure displacement (bottom) obtained at T = 10 ms. We note that the convergence rates for all variables are close to one, agreeing with theoretical predictions.

Example 3
In this example, we investigate the effects of the slip rate α on the fluid flow in a straight channel, for two different values of fluid viscosity. Here we no longer assume that the fluid domain is fixed. Namely, the full nonlinear coupling is implemented.
The model. We consider the full FSI model described by the Navier-Stokes equations, as specified in Section 2, in a straight channel of reference radius R = 0.5 and length L = 5. Since the problem is axially symmetric, we consider only the top part of the fluid domain. We set the reference fluid domain to be Ω = (0, 5) × (0, 0.5), with the top boundary representing a thin, elastic structure, whose reference configuration is Γ = (0, 5) × {0.5}.
As in Example 2, we prescribe symmetry boundary conditions at the bottom boundary Γ b . On Γ, the elastodynamics of the structure is modeled by the Koiter shell model (6.2)-(6.3) described in Example 1, allowing both horizontal and vertical displacement. In this example, the coefficients of the structural model are , The parameters used in this example are given in Table 2.
The inlet and outlet data. At the inlet x = 0, we prescribe Dirichlet data with the parabolic velocity profile: At the outlet x = 5, we prescribe Neumann data corresponding to the normal fluid stress equal to zero.
Initial data. Initially, zero fluid and structure velocity, with zero structure displacement were prescribed.
Computation. The simulations were run using ∆t = 10 −3 , until steady state is achieved for the prescribed inlet and outlet data. The domain is discretized using a uniform mesh with mesh size h = 0.025. As in Example 1, P 1 -bubble/P 1 elements were used to discretize the fluid problem, and P 1 elements to discretize the solid problem.
Results. We considered two different viscosities, µ F = 0.035 poise and µ F = 100 poise, corresponding to the viscosity of water and honey, respectively. In each of these cases, we found a solution corresponding to the slip rates α = 1, 10 −1 and 10 −2 . Figure 5 shows the 2D velocity magnitude for the solutions corresponding to α = 1 and the no-slip case, for the two values of viscosity µ = 0.035 (top) and µ = 100 (bottom). To obtain the no-slip solution, we used the kinematically coupled β scheme [13,17] with β = 1. The two solutions, shown in Figure 5, show the difference in the solution, especially near the boundary, where faster flow is observed near the boundary in the slip case. The difference between the two solutions is more prominent for the more viscous flow, as shown in the bottom panel of Figure 5. When the slip condition is prescribed, significant amount of flow occurs close to the fluid-structure interface. See top panel in Figure 6. Because of the slip, the flow near the fluid-structure interface causes significantly less structure displacement than in the no-slip case, as can be seen in the bottom panel of Figure 6. Thus, in the slip case, less of the fluid kinetic energy is spent on structure displacement, and more on dissipation due to slip friction. Indeed, Figure 6 shows the slip at the boundary, and structure displacement, respectively. Figure 6 shows also show that as α approaches zero, the slip solutions become closer to the solution for the no-slip case. We remark that significant boundary layer can be observed in these figures, especially near the inlet boundary, which is due to the clamped shell boundary condition, used to solve the structure equations. Finally, we compared the velocity profiles at the mid-point (x = 2.5) of the fluid domain obtained with the no-slip and slip conditions. Figure 7 shows the velocity profiles obtained with µ F = 0.035 poise (left) and µ F = 100 poise (right). When µ F = 0.035 poise, the velocity profile is parabolic and it does not change significantly when no-slip or slip conditions are imposed. However, for µ F = 100 poise the velocity profile obtained with no-slip condition is still parabolic, but it gets closer to plug velocity profile as the slip rate increases. We compare our numerical solution to an exact solution for the 2D Poiseuille flow in a rigid tube with the the Navier slip condition, and compare the exact solution to the numerically calculated solutions for different values of µ and α. We calculate the exact 2D Poiseuille solution in a rigid tube, satisfying the Navier slip condition to obtain: where α is the slip rate, p 0 and p L are the inlet and outlet pressure, µ is the fluid viscosity and R and L are the tube's radius and length, respectively. Since we do not have the inlet and outlet pressure data, to calculate the steady-state solution of our problem in which Dirichlet data are prescribed at the inlet, and zero normal stress at the outlet, we use the conservation of mass principle to recover the pressure drop (p 0 − p L ) from the inlet and outlet data, and from the form of the solution we just obtained. Thus, we set the flow rate determined by the inlet Dirichlet data to be equal to the flow rate of the steady-state solution inside the tube: and calculate an expression that determines the pressure gradient p 0 − p L : Once the pressure drop is calculated from u max , R, L, µ and α, we can plug it into the formula for the exact solution, to obtain the following results, plotted in Figure 8.
One can see that the exact solutions are in a very good agreement with the numerically calculated solutions. The biggest difference can be observed for µ = 100 poise in the slip at the wall, which is slightly bigger in the rigid case. This makes sense, since in the case µ = 100 poise, the structure displacement is bigger than in the case µ = 0.035 poise, and so more of the kinetic energy of the fluid is spent on the kinetic energy of the structure, and less on slip friction, which gives a smaller slip at the wall in the numerically calculated solution, as compared to the rigid wall exact solution.  Figure 7 to observe that for µ = 100 poise, the slip at the wall in the exact solution is visibly larger than in the moving wall case when more energy is spent on fluid domain motion, rather than on slip friction.
Based on the comparison with an exact solution in this example, as well as the convergence study in Examples 1 and 2, we note that the operator splitting error does not contribute significantly to the differences between the results obtained using different viscosities shown in this example.

Example 4
In the last example of this manuscript, we model the flow of a viscous, incompressible fluid through a bottleshaped domain, while the bottle is being squeezed using a time-periodic forcing, and the flow is additionally driven by gravity ρ F g. The case with no-slip and the case with slip prescribed at the lateral wall of the bottle are compared. The slip condition at the wall can be used to model the no-stick surface coating. We observe significantly different solutions, and show that the mass flow rate is, indeed, increased with the application of the no-stick coating. This approach can be used to, e.g., optimize surface coating depending on the type of viscous fluid considered. The bottle-shaped fluid domain is then rotated by −π/4 about the x−axis to obtain a tilted bottle, with gravity acting in the vertical, y direction. See Figure 9.
To describe the deformation of the wall, we use the linearly elastic membrane model proposed in [28,39]. The model, written in weak formulation, is given by the following: (6.12) Here, Π γ denotes the (surface) stress tensor, which, for a linearly elastic, isotropic structure reads: and ∇ γ (·) denotes the surface gradient, which can be computed using [9,28]: The symbol ⊗ denotes the tensor product, and I is the identity operator.
To model the time-periodic squeezing of the bottle at the lower, wider half, we use the following external forcing applied to the external surface of the linearly elastic membrane in the normal direction: n, x ∈ (0, 5). (6.13) Now, the dynamic coupling condition is given by equation (6.12), with the forcing term f replaced by the jump in the normal stress exerted by the fluid on one side, and the external forcing on the other: (6.14) The dynamic condition is split between the fluid and structure sub-problems as described in (3.4)-(3.9), with the integral containing the forcing term f ext used in the structure sub-problem (3.15) as the forcing in (3.4). The other parameters used in this example are described in Table 3.
The inlet and outlet data. At the inlet, zero fluid velocity is prescribed u max = 2.5 cm/s, while at the outlet we prescribe "free outflow" by imposing the zero normal stress data, as in the previous examples.
Computation. The mesh is discretized using 4000 elements; P 1 -bubble/P 1 elements are used to discretize the fluid problem and P 1 elements are used to discretize the solid problem. The simulations are run until the final  time T = 5 s is reached using ∆t = 5 × 10 −4 . We solve the problem by imposing both the no-slip condition, and the Navier slip condition with slip rate α = 10 −1 .
Results. Figure 9 shows the pressure in the bottle obtained using simulations with the no-slip and slip lateral boundary condition at times t = 4, 4.25, 4.5 and 4.75 s. As expected, we observe higher pressure during the squeezing part of the external force cycle, which, based on formula (6.13) corresponds to 4 < t ≤ 4.5. Notice, however, that since the squeezing rate is higher at t = 4.25 than at t = 4.5, the pressure at t = 4.25 is higher than at t = 4.5. Figure 10 shows the corresponding velocity magnitudes. This figure should be compared to Figure 11 where the velocity profiles along the bottle are shown. We observe that during the squeezing part of the cycle, the velocity is positive, while during the relaxation part of the cycle, the velocities are negative. Higher velocity occurs in the throat of the bottle, as expected. Velocity profiles at x = 1, 2.5 and 4 cm, at times t = 4, 4.25, 4.5 and 4.75 s, obtained using slip and no-slip conditions are shown in Figure 11. We observe that larger deviations in velocity close to the wall and in the middle of domain occur with no-slip condition, while slip condition gives rise to profiles more similar to plug velocity profile.
We calculated the corresponding mass flow rates at the outlet of the bottle for both the no-slip and the slip case, and obtained that the mass flow rate for the Navier-slip case is bigger than for the no-slip case. In particular, the positive mass flow rate at the outlet of the bottle for the slip case was 1.55 cm 2 /s and for the no-slip case 1.23 cm 2 /s. Thus the Navier slip boundary condition increases the mass flow rate at the outlet by 21%, for the flow mostly driven by gravity in the bottle angled at 45 degrees, with a slightly squeezed wall near the bottom, causing the change in the radius at the squeezing location of less than 15%.

Conclusions
We presented a novel numerical method for the interaction between an incompressible, viscous fluid and a thin structure (described in a lower-dimensional space), where the two are coupled by imposing balance of contact forces, continuity of velocities in the normal direction (no penetration), and the Navier slip condition in the tangential direction. The energy estimates showed that the proposed numerical method is unconditionally stable. We also presented error analysis showing first-order convergence in time and optimal convergence in space, obtained under a time-step condition. The theoretical results were computationally verified in Example 1, where the predicted rates of convergence in time were obtained by comparing the numerical solution with an exact, time-dependent solution of a simplified manufactured problem, and in Example 2, where the rates were computed on a moving domain problem. The third example was focused on a study of a viscous, incompressible fluid flow in a straight 2D cylinder with deformable walls, and the Navier slip condition at the lateral boundary, with different slip rates α for two different fluid viscosities (µ f = 0.035 poise (water), and µ f = 100 poise (honey)). We showed, among other things, that our numerical method produces solutions that approach, as the slip rate α approaches zero, the solution of the FSI problem with no-slip. Furthermore, in this case we derived the exact solution in a 2D rigid tube with the Navier slip condition, and showed that our numerical solution compares very well with the exact solution, with the differences attributed to structure displacement (which is more pronounced in the case µ f = 100 poise), captured by the numerical solution. In the fourth example we studied the flow of a viscous fluid in a ketchup bottle with and without a no-sick coating, i.e., with the Navier slip and with the no-slip boundary condition. The flow was driven mostly by gravity, with a slight push due to the squeezing of the bottle at the wider end. We showed how the solutions differ in pressure and velocity, and showed that at the outlet of the bottle, the mass flow rate with the Navier slip is higher (by 21%) than the mass flow rate in the bottle with the no-slip boundary condition. Proof. See [17] for proof. Proof. The last three inequalities follow directly from approximation properties (5.12)-(5.16). See [17] for more details.