A UNIFORM DISCRETE INF-SUP INEQUALITY FOR FINITE ELEMENT HYDRO-ELASTIC MODELS

. A seminal result concerning ﬁnite element (FEM) approximations of the Stokes equation was the discrete inf-sup inequality that is uniform with respect to the mesh size parameter. This inequality leads to optimal error estimates for the FEM scheme. The original version pertains to the Stokes system with non-slip boundary condition on the entire boundary. On the other hand, in ﬂuid-structure interaction problems, the interface dynamics between the ﬂuid and the solid satisﬁes velocity and stress matching constraints. As a result, the pressure variable is no longer determined up to a constant and becomes subject to non-homogeneous Dirichlet conditions on the common interface. In this context, we establish a uniform discrete inf-sup estimate for a ﬂuid-structure FEM implementation based on Taylor-Hood elements, and use this inequality to verify some stability and error estimates of this numerical scheme. An added beneﬁt of this framework is that it does not require the Poisson-equation approach to solve for the pressure variable.


1.
Introduction. The finite element method (FEM), formulated in the mid-1950's (though with many ideas dating back to early 1940's), was originally introduced for structural analysis problems. During the active development of this technique in the subsequent years, its success with approximations of potential flows (e.g., [54,27]) rendered the FEM as a prospective valuable tool for computational fluid dynamics. Rather than attempting yet another historical overview, for a list of seminal contributions and further developments, we refer the reader to the introduction and the references in [35] and the book [34]. It goes without saying that the subject has vastly expanded in the course of the following decades.
In step with the applications of the FEM to fluid flows, interest in hydro-elastic models was also gaining momentum. From the PDE systems proposed in [23, pp. 328-329] and [44, p. 121], originally motivated by biological applications, fluidelasticity models emerged into a research field of their own; see, for instance, [46,28] and the many references therein. Besides numerous papers accounting for the nonlinear aspects of the dynamics (to name a few: [36,51,24,16,9,42,13,40,38]; for plates and shells in fluid flows see, e.g., [21,47,20,48] and the citations therein), research efforts continue to pursue linearized systems as well; among the relatively more recent contributions are [1,3,4,19,14,2,12].
This work investigates the convergence properties of FEM approximations to fluid-elasticity models. In this context, let us mention a few, albeit by no means exhausting the background literature, foundational articles. A widely popular mixed P 2 /P 1 FEM scheme for the Stokes flow was among the pioneering contributions to the subject developed by Hood and Taylor in [52] (other element types have also been considered starting with the contributions by Fortin [31] and by Crouzeix & Raviart [25]; for subsequent developments and a comprehensive discussion see the tractate Boffi, Brezzi & Fortin [15]). An interesting variation-fundamentally different in several aspects-of the Taylor-Hood approach and the corresponding error estimates were developed by Glowinski and Pironneau in [35], which in part served as a basis for rigorous error analysis of the original Taylor-Hood scheme: see Bercovier and Pironneau [10] and the subsequent simplification by Verfürth [53] (also Fortin [32]).
The quantification of the convergence rates of the FEM for the Stokes system relies on the uniform, with respect to the mesh parameter, discrete "inf-sup" inequality for the associated bilinear form. The inf-sup framework itself could be chronologically attributed to Banach in the form of the Closed Range Theorem, Nečas [49,Thm. 3.1], Ladyzhenskaya in the context of the Stokes problem (though the inequality is not explicitly stated, see for instance the essential technical result [41, Sec 2.1, p. 24]), Babushka [6] and Brezzi [17] with applications to finite element problems.
The same inf-sup estimate is needed in the more involved case of fluid-elasticity interactions. However, this time the functional framework is different. The non-slip condition cannot be enforced on the interface between the fluid and the solid; this feature enlarges the functional fluid-pressure framework from The pressure variable is no longer determined up to a constant and becomes subject to non-homogeneous Dirichlet conditions on the common interface as dictated by the matching of the stresses. The discrete infsup condition for the corresponding finite element scheme was explored by Du et al. [29,Thm. 3.4,p. 12] based on piecewise-linear approximations of the velocity variable. Another numerical algorithm was later developed by Liu [45,Sec. 4,p. 7256] using the associated Poisson equation to determine the pressure. In this article, we establish a discrete inf-sup inequality for a fluid-structure model based on P 2 /P 1 Taylor-Hood elements and without employing the pressure-Poisson-equation approach.
1.1. Goals and outline. The purpose of this note is to establish the uniform infsup inequality as presented below in Theorem 3.1 and apply it to error estimates for the so called "semi-discrete" FEM fluid-structure model-Corollary 1. The argument borrows its strategy from the classical results of [10] and [53] and some of the requisite technical lemmas can actually be reused, such as the "weak" Brezzi inequality in Section 4. However, several adjustments are necessary to complete the proof in the new functional setting as shown in Section 5.
In addition, a fully discrete implementation must also approximate the elliptic bilinear form itself in the variational formulation. This complication allows for stability estimates, however, does not yield convergence in the original norms-see Figure 1. A sample 2D representation of the fluid-structure geometry. Section 3.2. Interestingly, the few particular numerical examples considered in the existing studies [3,11] do reproduce the rates of Corollary 1.

1.2.
Notation. Vector-value functions with components in a generic Banach space X will be denoted using bold-face symbols, such as L 2 (X) = [L 2 (X)] d , for d = 2, 3, etc. Given the fluid domain Ω f , the form (·, ·) will stand for the standard inner product on L 2 (Ω f ) or L 2 (Ω f ). In turn, (·, ·) s will indicate the same for the solid domain Ω s . The space H 1 Γ f (Ω f ) will denote the subspace of the W 1,2 (Ω f ) whose functions have zero trace on Γ f ⊂ ∂Ω f ; analogously for vector-valued functions Prefix D will stand for the differential, so Du..Dv, for example, is the Frobenius product of the Jacobi matrices Du and Dv of vector-valued functions u, v. The product (Du, Dv) will correspond to Ω f Du..DvdΩ f . Similarly, (Dφ, Dψ) s := Ωs Dφ..DψdΩ s .
The norm · will be the one in L 2 (Ω f ) or L 2 (Ω f ); otherwise, an appropriate subscript will be used to indicate which space the norm pertains to. We will also invoke the semi-norm for the space H 1 (Ω f ) given by Symbol S X will denote the unit sphere in the indicated Banach space X. : We consider the model linearized around rest, whence the domains Ω s and Ω f are time-independent (for a linearization about a solution to a free moving boundary problem see [14]). Denote their boundaries by Γ s := ∂Ω s and Γ f : The vector field n f will denote the outward unit normal field to the fluid domain on ∂Ω f , whereas n s will denote the outward unit normal to the solid domain on Γ s as depicted in Figure 1.

GEORGE AVALOS AND DANIEL TOUNDYKOV
The system of interest is the coupling of the Stokes and the linearized isotropic elasticity equations ( [44, p. 121], [28], [1]) : The stress operator σ is either vectorial Laplacian ∆ or the isotropic elastic stress tensor div Cε(w), ε(w) with the Lamé parameters λ, µ. The difference between the stress operators ultimately amounts to just using a suitable equivalent inner product on H 1 (Ω s ) and the associated normal strain in the boundary conditions, but the analysis remains conceptually identical. Henceforth, for ease of presentation we will focus on the "diagonal" version: The dynamics equations are accompanied by the non-slip conditions on the outer boundary Γ f and the appropriate stress and velocity matching conditions on the interface Γ s : The initial data comes from the state space with We stress again that in the definition of H f l the zero normal component condition is enforced only on the exterior non-slip boundary Γ f , but not on the entire boundary of the fluid domain ∂Ω f . The semigroup well-posedness of (1)-(3) was discussed at length in [1,3,4] (see also [12] for the semigroup result on a related, but more complex model linearized around a steady regime).

Geometry of the domain.
For the well-posedness theory it is sufficient to assume that Ω f and Ω s are of class C 2 , but polygonal boundaries should also be considered when dealing with the FEM framework.
The variational formulation of the problem discussed here only requires the Stokes component of the solution to be compatible with the velocity u in H 1 (Ω f ) and the pressure p in L 2 (Ω f ). For the Stokes operator with non-slip conditions a polyhedral domain is sufficient [26, p. 76]. Note, however, that by the nature of the problem, as shown in Figure 1, it is no longer expected that the fluid domain Ω f is convex, thus, for instance, the results [26, p. 76] do not guarantee that under a smooth forcing one can expect the fluid velocity with regularity u ∈ H 3 (Ω f ) and the pressure p ∈ H 2 (Ω f ). This detail bears possible limitations on the applicability of the convergence results of Corollary 1.
Remark 1 (Elliptic regularity in strong formulation). The definition of the semigroup generator for the full hydro-elastic system involves "pressure-determining" operators [3, p. 264] given by harmonic extension maps for the Dirichlet and Neumann boundary data respectively. By duality, following the transposition argument of Lions and Magenes [43, Ch 2. Sec. 6], or the weak formulation of Babuska et al. [7,8], these operators exist if the Poisson problems with homogeneous Neumann or Dirichlet conditions define topological isomorphisms of L 2 (Ω f ) onto the corresponding subspace of H 2 (Ω f ). However, even in 2D polygonal setting, such elliptic smoothing in general holds only for convex domains [37,Coro. 4.4.3.8 and Lem. 4.4.35,. So the strong semigroup interpretation of the problem holds on any domain which admits such harmonic extension maps, but in general it may, possibly, be lacking regularity in polyhedral approximations of domains like the one depicted in Figure 1. In particular, any polygonal approximation of the boundary interface Γ s will have re-entrant corners with respect to the fluid domain Ω f . (1)-(3) can be recast [3] as a saddle-point problem, where in process one must invoke a different from classical pressureelimination procedure since the Leray projector is not compatible with the space H f l , due to the lack of the boundary conditions on Γ s .

Variational formulation. System
First, express (1)-(3) as an evolution system on state space H, as in (3), for a dissipative semigroup generator A, with appropriate domain D(A) [1,3,4]. The resolvent problem of finding Y ∈ D(A) such that for λ > 0 and given Y * = (w * 1 , w * 2 , u * ) ∈ H, reduces to finding the velocity The symmetric bilinear form a λ is given by [3, (4.24), p. 272]: on where the restriction to Γ s indicates the trace map, D λ , for λ > 0, is the solution map to the Dirichlet harmonic extension problem The second bilinear form is And the functional F is defined by [3, (4.25), p. 273]: wherein Thus ∆ −1 λ is a solution operator for a Poisson-type problem. An important distinction from the non-slip Stokes flow model is that If this form vanishes for all v ∈ H 1 Γ f (Ω f ), then q must be not merely constant, but, in fact, zero in Ω f . Hence solutions to the hydro-elastic system (5) are determined uniquely, whereas the Stokes model determines the pressure only up to a constant.
System (5) has a solution if (e.g., see [39, Thm. 3.1.5, p. 116]) a λ is elliptic (in fact, the ellipticity is only needed on an appropriate subspace with respect to b) and for some c > 0 The latter can be equivalently restated as The ellipticity of a λ is apparent from the definition (6) for λ > 0, and the verification of the continuous inf-sup inequality for (5) can be found in [3].
In the FEM version of (5), one has to find a single c satisfying (10) for every mesh in the appropriate family on domain Ω f .

Main results.
Assume that Ω f is a polyhedron. By the structure of the domain, that implies that Ω s is a polyhedron as well. For a shorthand, denote on Ω = Ω f ∪ Ω s such that every element K ∈ T h resides either in Ω f or in Ω s . Specifically we consider approximations of X , Y and M based on triangular or tetrahedral Taylor-Hood elements P 2 /P 1 . For K ∈ T h let P m (K) denote the set of polynomials of degree ≤ m on K. Then define the approximating spaces by (The slight abuse of notation K ∈ T h ∩ Ω f or s means K ∩ Ω f or s has positive measure.) Note that the space Y h does not explicitly enter the variational formulation (5), though it will come up within the "fully discrete" model discussed below in Section 3.2.
3.1. "Semi-discrete" scheme. To start with, assume that the bilinear form a λ (6) and the functional F (8) can be evaluated exactly. The discretization, thus, shows up only in the function spaces. Accordingly, consider the following approximation of system (5): find u h ∈ X h and p h ∈ M h such that For exposition purposes we present the result in the two-dimensional formulation, but there are no technical obstacles to its extensions to three dimensions. Thus, henceforth, Ω ⊂ R 2 and each K ∈ T h is a triangle.
Theorem 3.1 (Discrete uniform inf-sup inequality). Assume that for each h > 0, every element in T h that is supported in Ω f has at least one vertex not in Γ f . Then there is a constant C * > 0 independent of h ∈ (0, h 0 ), some h 0 > 0, such that Before addressing the proof we can present some convergence estimates for the hydro-elastic model: Corollary 1. If functions u ∈ X and p ∈ M solve (5) and u h ∈ X h , p h ∈ M h solve (11), then under the assumptions of Theorem 3.1, there exists c > 0 independent of h > 0 small, such that (12) Thus, if a priori u ∈ H 3 (Ω f ) and p ∈ H 2 (Ω f ), we recover the classical error estimates, as in [10]: Proof. The first inequality follows directly from [18, Thm. 2.1, p. 60], noting that in this case the pressure is in fact determined in the space L 2 (Ω f ) rather than L 2 (Ω f )/R. The ensuing convergence estimate (13) The bound on ∇(p − q h ) for the piecewise-linear interpolant q h ∈ M h of p ∈ H 2 (Ω f ) is of the order o(h). Hence it suffices to establish a suitable bound on ∇(q h − p h ) . The difference of the equations (11) and (5) with the test function v h ∈ X h gives in combination with the Poincaré inequality, it leads to By the "weak" Brezzi inequality of Lemma 4.2 take the supremum over v h ∈ X h to conclude By the convergence result (13) and according to the interpolation estimate for the interpolant q h ∈ M h of p ∈ H 2 (Ω f ), the right-hand side in (15) is of the order h −1 o(h 2 ) = o(h), which completes the proof.

3.2.
Fully discrete scheme. Another crucial difference from the Stokes problem is that the bilinear form a λ in (6) and functional F in (8) incorporate the Dirichlet harmonic extension operator D λ into the solid domain Ω s , as well as the solution to the Poisson-type problem for ∆ λ . In practice these cannot be exactly evaluated.
Hence the fully discrete problem, rather than having the form (11), actually reads The critical part in estimating the suprema on the right of (17) boils down to considering the difference where D λ is the harmonic extension (7), and D h λ is its finite element approximation. Specifically, g h := D h λ (z h ) for z h ∈ X h is the solution in the space Y h (the finite element subspace of H 1 (Ω s )) to the problem [30,Pro. 3.26,p. 125], but the convergence as h 0 requires a more careful consideration. If we define γ h s (z) to be the trace approximation on Γ s obtained from a suitable interpolant, e.g., [50], of the function z ∈ H 1 (Ω f ), then we have the classical convergence estimate, for instance [30,Coro. 3.29,p. 126], However, with respect to the H 1 (Ω f ) norm only, as the normalization on the righthand side of (17) necessitates, we merely have an upper bound: Thus we can infer that the supremum terms on the right hand side of (17) are bounded uniformly in h, thus verifying a stability estimate for the numerical scheme. But unless one considers convergence in lower-order norms, there is no expectation of decay of these bounds as h 0 in general. Despite this discrepancy in the convergence analysis between the semi-discrete model (5) and the fully discrete model (11), the few particular examples tested in [3] and [11] do exhibit the decay rates of the semi-discrete scenario, as furnished by Corollary 1.
The remainder of the paper is devoted to the proof of Theorem 3.1.

4.
Revisiting the "weak" Brezzi inequality. We will henceforth only deal with the elements in the fluid subdomain Ω f . For such an element K ∈ T h denote by p 1 , p 2 , p 3 its locally numbered vertices. We assume the numbering goes counterclockwise around the element. For the edge from p i to p (i mod 3)+1 , define m i to be its midpoint, and let τ i be a CCW-oriented unit tangent vector along the edge: See a sketch of a sample element in Figure 2.
We start with an elementary result: There exists κ > 0 independent of h > 0 such that for each q h ∈ M h and each element K ∈ T h , the function ∇q h (which is, in fact, constant on K) satisfies In particular, we can choose κ = 2/(sin θ min ) 2 where θ min is a lower bound on the angles between any adjacent element edges for the elements of the mesh family Proof. For q ∈ M h and K ∈ T h the gradient ∇q h is a constant vector field on K so denote it by ∇q h K = (c i ). Pick unit tangential vectors on any two boundaries of the element. For a more compact notation label them as τ andτ . Of course, ∇q h can be reconstructed in terms of these two linearly independent vectors. We only need to remark why the relationship is uniform in h and across all elements. Solving for the first component c 1 of ∇q h amounts to finding α, β such that Thus, from (19) we infer The determinant τ 1τ2 − τ 2τ1 is the area of the parallelogram spanned by the unit vectors τ 1 ,τ 2 , so it is the sine of the acute angle between them. So if θ 0 is the minimal angle in K then we have max{|α|, |β|} ≤ | sin θ 0 | −1 .
Similar analysis holds for the second component c 2 of ∇q h . Thus, when the mesh family is shape-regular with θ min being a lower bound on the angle, the conclusion of the Lemma follows with κ = 2(sin θ min ) −2 , independently of h.
Remark 2. Because the finite element pressure subspace is that of piecewise linear functions, and since in our case the supremum is taken over a larger space X h (now allowing for functions that are not necessarily zero on the boundary portion Γ s ), this Lemma actually follows from [10, Prop 1., p. 214]. We present a sligthly modified proof for completeness and also to point out that the conditions on the mesh can be relaxed a little admitting elements with two edges on Γ s , albeit then at the obvious expense of convexity of Ω s , since with triangular elements it would acquire a reentrant corner. Proof.
Step 1. Test function v h . We employ the fluid test functions, with a slight modification, that were used in the proof of [30,Lemma 4.23,p. 193]. Fix Piecewise quadratic v h is uniquely determined by its values at the vertices and the edge midpoints, hence the definition is consistent. It is helpful to note the following Since the considered mesh family is, in particular, shape-regular, then the isomorphisms between the elements and the reference element are uniformly (operatornorm) bounded in either direction independently of h and K. Hence the Lagrange local basis functions are pointwise bounded independently of h or K. Therefore, we have for C 3 = 3 · C 2 independent of h. We restate this as Step 2. The estimate. Integrate by parts: where recall n f is the exterior normal field to the fluid domain. Label by E j the edges {E j = [d j , f j ]} that comprise Γ s . Then the boundary integral in (22) gives: We could have defined v h to be zero on all the nodes on Γ s which would eliminate this boundary integral. That would impose a minor restriction of each element to have one vertex outside Γ s . Alternatively, note that the given definition works just as well. Specifically, the polynomial q h n f · v h restricted to edge segment (d j , f j ) reduces to a cubic in one variable. For polynomial φ of degree 3, we can exactly obtain the value from the Simpson's rule which gives Thus (22) simplifies to For K ∈ T h the product q h div v h K is a quadratic polynomial. So we can use the exact quadrature the sum being taken over all the edge midpoints of K. Apply this identity to equation (23) and use the definition (20) of v h at the midpoints to obtain By the assumption that each K contains at least one vertex not in Γ f (though it could have its vertices in Γ s ), we know that each sum in the preceding identity involves at least two edge midpoints on the boundary K. According to Lemma 4.1 for each K the sum on the right dominates the constant field ∇q h K . Hence Invoke (21) to arrive at From which the conclusion of the Lemma follows with C = 1/(3κ √ C 3 ) .
where c is independent of h > 0 small. So Lemma 4.2 immediately yields This is not quite the estimate we are after but we will take advantage of this inequality later.
(on a smooth domain we would have simply used n f itself). Now consider the system for an unknown function v Note that the compatibility condition Using the assumption q h L 2 (Ω f ) = 1 we can just write Now we consider a suitable approximation of v in the space X h . Using [22, Thm. 1] there exists a projection operator R h : X → X h such that for some c > 0 and any v ∈ X |v − R h v| k,Ω f ≤ ch 1−k |v| 1,Ω f , k = 0, 1 .
In particular, Let α := Γs µ · n f dΓ s , as the coefficient in (26). Using the property (25) of the vector field µ, the first term on the right-hand side gives Invoking (27), we also obtain Next, the second term on the right-hand side of (29) after integration by parts gives With the help of projection error estimate (28) the term I in (31) can be estimated directly as: For the term II in (31), we first bound the L 2 (Γ s )-norms by H 1/2 (Ω f ) norms and interpolate the latter between L 2 (Ω f ) and H 1 (Ω f ). For the H 1 (Ω f )-norm of (R h v − v) we can, up to a Poincaré constant, use the seminorm | · | 1,Ω f since these functions vanish on Γ f . For estimates of q h use the fact that q h L 2 (Ω f ) = 1. Then invoke the projection error estimate (28): for any ε > 0. Thus, (31) gives us Now plug this estimate, along with (30), into (29): Fix ε < 2C 1 , let h 0 > 0 be small enough so that Further define C * 3 := C 2 + C 3 /(2ε). Then for all h ∈ (0, h 0 ) we have Next, recall from (24) that for small h Directly from the two boxed identities we get, respectively, C * 1 S q h ≥ C * 1 C * 2 − C * 1 C * 3 h ∇q h and C * 3 S q h ≥ C * 3 C * 1 h ∇q h Adding them together yields: S q h ≥ C * := C * 1 C * 2 (C * 1 + C * 3 ) independently of h ∈ (0, h 0 ) and q h ∈ M h , for a suitably small fixed h 0 > 0. This step completes the proof of Theorem 3.1.