Dordt Digital Collections Dordt Digital Collections

Abstract We present qualitative and numerical results on a partial differential equation (PDE) system which models a certain fluid-structure dynamics. Wellposedness is established by constructing for it a nonstandard semigroup generator representation; this representation is accomplished by an appropriate elimination of the pressure. This coupled PDE model involves the Stokes system which evolves on a three dimensional domain O coupled to a fourth order plate equation, possibly with rotational inertia parameter ρ >0. This plate PDE evolves on a flat portion Ω of the boundary of O. The coupling on Ω is implemented via the Dirichlet trace of the Stokes system fluid variable - and so the no-slip condition is necessarily not in play - and via the Dirichlet boundary trace of the pressure, which essentially acts as a forcing term on Ω . We note that as the Stokes fluid velocity does not vanish on Ω , the pressure variable cannot be eliminated by the classic Leray projector; instead, it is identified as the solution of an elliptic boundary value problem. Eventually, wellposedness of the system is attained through a nonstandard variational (``inf-sup") formulation. Subsequently we show how our constructive proof of wellposedness naturally gives rise to a mixed finite element method for numerically approximating solutions of this fluid-structure dynamics.

1 The PDE and Setting for Wellposedness One of our main objectives in this work is to provide a proof for semigroup wellposedness with respect to the fluid-structure partial differential equation (PDE) model considered in [6] -see also [5] and [7].The proof here will be wholly different than that originally given in [6], and has the virtue of giving insight into a mixed finite element method (FEM) formulation so as to numerically approximate the solution of the fluid and structure variables.A numerical analysis involving this fluid-structure FEM will constitute the second part of this work.Throughout, we will consider situations in which either the "Euler-Bernoulli" or "Kirchhoff" plate PDE is in place to describe the structural component of the fluid-structure model (only the Euler-Bernoulli is considered in [6]).The geometrical situation will be identical to that in [6].We state it here verbatim: O ⊂ R 3 will be a bounded domain with sufficiently smooth boundary.Moreover, ∂O = Ω ∪ S, with Ω ∩ S = ∅, and specifically Ω ⊂ {x = (x 1, x 2 , 0)} , and surface S ⊂ {x = (x 1, x 2 , x 3 ) : x 3 ≤ 0} .
(So when ρ = 0, Euler-Bernoulli plate dynamics are in play; when ρ > 0 we have instead the Kirchhoff plate.)Here, the space of initial data H ρ is defined as follows: Let and Therewith, we then set Moreover, let A D : L 2 (Ω) → L 2 (Ω) be given by If we subsequently make the denotation for all ρ ≥ 0, then the mechanical PDE component ( 2)-( 3) can be written as P ρ w tt + ∆ 2 w = p| Ω on (0, T ).
Using the characterization from [11] that then we can endow the Hilbert space H ρ with norm-inducing inner product ω 1 , ω 2 , f , ω1 , ω2 , f Hρ = (∆ω 1 , ∆ω 1 ) Ω + (P ρ ω 2 , P where (•, •) Ω and (•, •) O are the L 2 -inner products on their respective geometries.We note here, as there was in [6], the necessity for imposing that wave initial displacement and velocity each have zero mean average.To see this: Invoking the boundary condition (6) and the fact that normal vector ν = [0, 0, 1] on Ω, we have then by Green's formula, that for all t ≥ 0, And so we have necessarily, This accounts for the choice of the structural finite energy space components for H ρ , in (10).
As we said, our proof of wellposedness hinges upon demonstrating the existence of a modeling C 0 -semigroup e Aρt t≥0 ⊂ L(H ρ ), for appropriate generator A ρ : H ρ → H ρ .Subsequently, by means of this family, the solution to (2)- (7), for initial data [w 1 , w 2 , u 0 ] ∈ H ρ , will then of course be given via the relation Our particular choice here of generator A ρ : H ρ → H ρ is dictated by the following consideration: If p(t) is a viable pressure variable for ( 2)- (7), then pointwise in time p(t) necessarily satisfies the following boundary value problem: To show the validity of ( 15)-( 17): taking the divergence of both sides of (4) and using the divergence free condition in (5) yields equation (15).Moreover, dotting both sides of (4) with the unit normal vector ν, and then subsequently taking the resulting trace on S will yield the boundary condition (17).(Implicitly, we are also using the fact that u = 0 on S.) Finally, we consider the particular geometry which is in play (where ν = [0, 0, 1] on Ω to establish (16)).Using the equation (2) and boundary condition (6), we have on Ω which gives (16).
The BVP ( 15)-( 17) can be solved through the agency of the following "Robin" maps R ρ and Rρ : We define Therewith, we have that for all real s, (See e.g.[14].We are also using implicity the fact that P −1 ρ is positive definite, self-adjoint on Ω, and moreover manifests elliptic regularity.) Therewith, the pressure variable p(t), as necessarily the solution of ( 15)- (17), can be written pointwise in time as where These relations suggest the following choice for the generator A ρ : H ρ → H ρ .We set (c.f. the generator and domain described in an earlier version of [6] in arXiv:1109.4324.)(Note that as ∆u ∈ L 2 (O) and div(∆u) = 0, then by Theorem 1.2, p. 9 in [16], we have the trace regularity and so the pressure term In short the domain of In what follows, we will have need of solution maps for certain inhomogeneous Stokes flows.To wit: For given φ ∈ H We note that the classic compatibility condition for solvability is satisfied, and that pressure variable π is uniquely defined up to a constant; see e.g., Theorem 2.4, p. 31 of [16].Then by Agmon-Douglis-Nirenberg, we have f , π ∈ L(H (see Proposition 2.2, p. 33 of [16]).
In a similar way, we define for fluid data Again by Agmon-Douglis-Nirenberg we have These two fluid maps will be invoked in the proof of Theorem 1 below, which will yield solvability of the following resolvent equation for λ > 0: Namely, for given After setting the pressure variable p = G ρ,1 (w 1 ) + G ρ,2 (u) in (24), the resolvent equation ( 33) is equivalent to the following system: In particular, it will be seen in Theorem 1 that the solution variable w 1 in (33) can be recovered through finding the unique solution pair [w 1 , c] ∈ H 2 0 (Ω) × R which solves where: Theorem 1 (i) The operator A ρ : H ρ → H ρ is maximal dissipative.Therefore by the Lumer-Phillips Theorem it generates a C 0 -semigroup of contractions e Aρt t≥0 on H ρ .
(ii) Let λ > 0 and [w * 1 , w * 2 , u * ] ∈ H ρ be given.(By part (i), there exists [w 1 , w 2 , u] ∈ D(A ρ ) which solves (33).)Then the structural solution component w 1 and constant component c of the associated pressure term p can be characterized as the solution pair of the variational system (41).Subsequently the remaining unknown terms are given by (after utilizing the maps f and π in (29) and μ and q in (31)).
2 Proof of Theorem 1

Proof of Maximality
In what follows we will make use of the Babuška-Brezzi Theorem.We state it here directly from p. 116 of [13].
Theorem 2 (Babuška-Brezzi) Let Σ, M be Hilbert spaces and Assume that a(•, •) is Z-elliptic, i.e. there exists a constant α > 0 such that Assume further that there exists a constant β > 0 such that Then if κ ∈ Σ and ∈ M , there exists a unique pair b(σ, q) = ( , q), for every q ∈ M.
For λ > 0, we will show that Range(λI Using the structural component (34) and (35), we then have the boundary value problem as well as the fluid system Since [w 1 , w 2 , u] ∈ H ρ , we also have that With this compatibilty condition in mind, by way of "decoupling" the systems ( 52) and (53), we proceed as follows.We apply P ρ to the mechanical equation in (52) and then multiply by test function φ ∈ H 2 0 (Ω).Subsequently integrating over Ω then yields (where (•, •) Ω might also indicate the duality pairing between H 2 0 (Ω) and H −2 (Ω)).Upon integration by parts and using φ ∂Ω = ∂φ ∂ν ∂Ω = 0, we then have (Note that in obtaining this expression we have again used div(u) = 0 on Ω and normal vector Now, from ( 53) and ( 54), we can use the fact that in terms of the maps in ( 29) and (31) above, we can write fluid variables u and p of ( 52)-( 53) as for some (to be determined) constant c, justifying (43).Applying this representation to (55), subsequently integrating by parts, and using the mapping in (29), we then have for every φ ∈ H 2 0 (Ω), Using now the representations in (56), we have then Now using ( 29) and (31) as well as the fact that π, q ∈ L 2 (O) R , we rewrite this expression as This variational relation and the constraint (54) establish now the characterization of the range condition (51) with the mixed variational problem (41).This characterization, along with (34) and (56), establishes Theorem 1(ii).
Note in particular that from the second equation of (41).In turn, we recover w 2 , u, and p via From ( 29) and (31), the variables u and p solve the Stokes system (53). (64) Moreover, from (57) and (29) we have An integration by parts and passing of the adjoint of P 1/2 ρ yields, As variables u and p solve the Stokes system (53), we thus attain the relation (We have also implicitly used (29) and the remark after (55).)In particular, this holds true for φ ∈ D(Ω).Thus we have the distributional relation , and so we infer that w 1 satisfies (52).
Subsequently, we infer by elliptic theory that, as required by the definition of the fluid-structure operator where the (displacement) space is as given in (25).Finally, because u and u * ∈ H fluid , we have a fortiori from (53), Moreover, from (63) and (52), where in the last equality, (53) was again invoked.Thus, from (67) and (68), we have that the pressure variable p we have obtained by Theorem 2 solves As such, where G ρ,i are as given in ( 22) and ( 23).(Note that we are implicitly using the critical regularity from (27).)From ( 62), ( 63), ( 64), ( 65), (66), and (69), we have that the constructed variables [w 1 , w 2 , u] belong to D(A) and solve the resolvent equation (51).This concludes the proof of Theorem 1.

A Numerical Analysis of the Fluid-Structure Dynamics
The objective of this section is to demonstrate how the maximality argument which was given in Section 2.2 can be utilized to approximate solutions to the fluid-structure interactive PDE under present consideration.In particular, the numerical method outlined here solves the static problem resulting from the resolvent equations ( 34) -( 40), but this approach in principle can be modified to solve the time dependent problem in the same way as in [1] (see p. 276).We will outline here a certain numerical implementation of the finite element method (FEM) and provide convergence results for the approximation with respect to "mesh parameter" h.Finally, we will provide an explicit model problem as a numerical example.

Finite Element Formulation
In what follows, the three dimensional body O will be taken to be a polyhedron.Given a positive (and small) parameter h of discretization, we let {e } N h =1 be an FEM "triangulation" of O, where each element e is a tetrahedron (and so, among other properties, [2] and Figure 5 below).
(A) Relative to the "triangulation" of O, V h will denote the classic H 1 -conforming FEM finite dimensional subspace such that (See [2].)Subsequently to handle the inhomogeneity we specify the set (B) In addition Π h will denote the L 2 -FEM finite dimensional subspace for the pressure variable defined by (see [8] and [2]).Moreover, we let {ẽ } Ñh =1 be a FEM triangulation of the two dimensional polygonal region Ω, where each element ẽ is a triangle.(C) Similarly, X h will denote a conforming FEM subspace such that (see e.g.[8] and [15] for details of the explicit construction of these piecewise polynomials.As such, the basis functions which generate X h are "conforming", relative to fourth-order boundary value problems.) For the spaces V h , Π h and X h described above we will have need of the following discrete estimates relative to mesh parameter h: (A ) In regard to the H 1 (O)-conforming FEM space V h in (71) we have the following estimate: (See Theorem 5.6, p. 224, of [2].) (B ) Similarly, in regard to the finite dimensional space Π h , we have the discrete estimate: For (See e.g., Corollary 1.128, p. 70, of [10].)(C ) Finally, with regard to the FEM space X h in (74), we have the following discrete estimates: (See estimate (5.82), p. 225, of [2].) The goal here is to find a finite dimensional approximation [w 1h , w 2h , u h ] ∈ X h × X h × Ṽh to the solution [w 1 , w 2 , u] ∈ D(A ρ ) of (51), as well as an approximation p h of the associated fluid pressure p.We shall see that these particular FEM subspaces are chosen with a view of satisfying the (discrete) Babuška-Brezzi condition relative to a mixed variational formulation, a formulation which is wholly analogous to that in (41) for the static fluid-structure PDE system (34)-(40).We further note that, by way of satisfying said inf-sup condition, it is indispensable that the structural component space X h be H 2 -conforming (see (74)).In addition, this mixed variational formulation for the coupled problem (34)-(40), like the mixed method for uncoupled Stokes or Navier-Stokes flow, allows for the implementation of approximating fluid basis functions (in V h ) which are not divergence free (see [4]).
In line with the maximality argument of Section 2.2, the initial task in the present finite dimensional setting is to numerically resolve the structural solution component of the PDE system (34)-(40).Namely, with reference to the bilinear and linear functionals a λ (•, •), b(•, •) and F(•) of (42), the present discrete problem is to find [w 1h , ch ] ∈ X h × R which solve: Assuming this variational problem can be solved uniquely, w 2h is immediately resolved via the relation (cf. ( 34)).Subsequently we can recover fluid and pressure approximations u h and p h from the discrete solution pair [w 1h , ch ] ∈ X h × R of (79).Indeed, to this end we will invoke the classic mixed variational formulation for Stokes flow, so as to approximate the fluid maps [ f (•), π(•)] and [μ(•), q(•)] of ( 29) and (31), respectively.(See [4].)Let bilinear forms ãλ (•, •) : Moreover, we define the standard Sobolev trace map γ 0 : Since γ 0 (•) is continuous and surjective, then for any φ ∈ H k−1/2 (Ω) we have the existence and uniqueness of an element in H k (O), denoted here as γ + 0 (φ), which satisfies Therewith, the classic mixed FEM for ( 29) is given as follows: With subspaces V h and Π h as given in ( 71) and ( 73) respectively, and given φ ∈ H 1/2 (Ω), find the unique pair Likewise, the classic mixed FEM for (31) is given as follows: For given By the Babuška-Brezzi Theorem, the two discrete variational formulations (84)-( 85) and ( 86)-(87) are well-posed; see [4].(In particular, with the so-called Taylor-Hood formulation in place -i.e, fluid approximation space V h consists of piecewise quadratic functions, and pressure approximation space Π h consists of piecewise linear functions -then the aforesaid inf-sup condition is satisfied uniformly in parameter h.) With the approximating solution maps (84)-(87) in place and assuming the structural component approximation [w 1h , ch ] is known, we then set (c.f.(43).)Now in regard to the variational problem in (79), one will in fact have unique solvability of this discrete problem, via the Babuška-Brezzi Theorem, provided that the following inf-sup condition is satisfied: But what is more, in order to ensure stability and ultimately convergence of the numerical solutions obtained by our particular FEM, it is indispensable that the "discrete" inf-sup condition (90) be uniform of parameter h > 0 (at least for h small enough).
In fact we have the following result: as defined in (42).Then for parameter h > 0 small enough one has the "inf-sup" estimate where C = ξ H 2 0 (Ω) − , and ξ is the solution of the boundary value problem (59).Here, > 0 can be taken arbitrarily small.
Proof of Lemma 3. We resurrect the elliptic variable ξ ∈ H 4 (Ω) ∩ H 2 0 (Ω) from the earlier maximality argument.Namely, ξ solves Then by Green's First Identity we have that ξ solves the following variational problem for all ψ ∈ H 2 0 (Ω): Let now ξ h ∈ X h denote the "energy projection" of ξ on X h .That is, ξ h satisfies the following discrete variational problem: The existence and uniqueness of the discrete solution ξ h ∈ X h follows from the Lax-Milgram Theorem, see e.g., [2], [8].Applying the discrete estimate (77) to the respective variational problems (92) and (93), we then have With these ingredients, ξ and ξ h , we then have sup after using (93) above.Continuing, we then have sup after using the estimate (94).Taking step size parameter now completes the proof.

Error Estimates for the Finite Element Problem
In what follows we will have need of the following result in [10], which will not be stated here in its full generality (see [10], Lemma 2.44, p. 104).
Lemma 4 With reference to the quantities in Theorem 2 above, let Σ h be a subspace of Σ, and let M h be a subspace of M .Suppose further that bilinear form a Also, assume that the following "discrete inf-sup" condition is satisfied: where β h > 0 may depend upon subspaces Σ h and M h .Let moreover (σ h , p h ) ∈ Σ h × M h solve the following (approximating) variational problem: (Note that the existence and uniqueness of the solution pair (σ h , p h ) follows from Theorem 2, in view of (98) and (99).)Then one has the following error estimates: Concerning the efficacy of our FEM for numerically approximating the fluid-structure system (34)-(40), we have the following: Theorem 5 Let h > 0 be the parameter of discretization which gives rise to the FEM subspaces V h , Π h , and X h of (71), (73), and (74), respectively.With respect to the solution variables 34)-(40) and their FEM approximations [w 1h , w 2h , u h , p h ], as given by ( 79) -( 80) and ( 88)-(89), we have the following rates of convergence: Proof of Theorem 5. We first establish parts (i) and (ii) together.Here we will combine Lemma 4 with the bilinear forms in (42).We take in Lemma 4 In addition, by Lemma 3, we can take Subsequently, with reference to the variational system (41) and the approximating system (79) we have from ( 101) (note that that second term from the right hand side of (101) is zero in this case because M h = M = R in this case).Appealing now to estimates (77), (78), and Theorem 1(i) we have the following error estimates: (a) If ρ = 0, This establishes Theorem 5(i).In view of (80), Theorem 5(ii) follows directly.
To achieve the estimates for the fluid variables, we first note that the error in the constant component of the pressure, given by the variational system (41), satisfies: (a ) If ρ = 0, we have upon combining (102) and ( 110) Consequently we can take α ≡ 1 and a = ãλ In addition, since the respective fluid and pressure spaces V h and Π h are piecewise quadratic and piecewise linear -i.e. the so-called Taylor-Hood formulation -then for h > 0 small enough we can take β h = β * , independent of small h > 0. (See e.g., Lemma 4.23,p. 193 of [10].)Thus with u and u h as given in (43) and (88) respectively, we then have associated to points in the mesh and the system is solved in this finite dimensional setting via a matrix/vector equation, see e.g.[2].First consider the plate system (41).This weak formulation takes place over H 2 0 (Ω) and thus the most natural choice for the discretization is a set of H 2 -conforming elements, see [15]; for a MATLAB implementation see [9].We use the quintic Argyris basis functions because they are the lowest order H 2 0 -conforming elements available and they ensure wellposedness of the discrete formulation of (41) because the inf-sup condition on the bilinear form b can still be satisfied.The Argyris basis functions have 21 degrees of freedom(DOF) for each triangle in the mesh, namely Lagrange DOF for function values at each vertex, Hermite DOF for ∂/∂x and ∂/∂y at each vertex, Argyris DOF for ∂ 2 /∂x 2 , ∂ 2 /∂x∂y, and ∂ 2 /∂y 2 at each vertex and one DOF at each edge midpoint for the normal derivative.
On the fluid domain a Stokes system is solved via a mixed formulation.Again one must be careful as to the choice of elements used to guarantee that the discrete problem remains well-posed.In this case we use the popular Taylor-Hood (P 2 /P 1 ) elements, see e.g.[10].
To derive the discretized version of (41) let X h denote the space of Argyris basis functions of dimension n with basis {φ i } n i=1 .Let w h = n i=1 α i φ i .Then for each φ j ∈ X h and r ∈ R we have with a λ , b and F defined in (41) -(42).By writing this for each φ j ∈ X h we build a linear system of the form A B B T 0 Note importantly, that in practice f and μ are not known exactly and so in reality (123) is actually created using a subroutine that numerically approximates the instances of f and μ that occur in a λ using a discrete mixed variational formulation on the fluid domain.The discretized versions of (29) and (31) take the same mixed form as (123), but with ã(µ, ϕ) and b(µ, q) defined as in (81) -(82).
Then for any λ > 0 the functions 40) for data defined by  Since the mesh is refined by a factor of 2 at each step, we compute log Error i Error i+1 / log(2).In the limit this ratio should approach the exponent of convergence, i.e.O(h k ).Now for smooth data (as we have here) the best possible convergence rate one could attain is k = 4 for the H 2 norm of w 1 which the numerical scheme does appear to attain (see Table 2).However, the H 1 and L 2 errors do not appear to improve to k = 5 and k = 6 respectively; this is possibly due to the (unavoidable) approximation of f and μ described above.In Figure 2 we see that already at the 3rd level of mesh the FEM approximation is indistinguishable from the true solution.Similarly, for the fluid approximation we have in Table 3 the errors in the fluid variables for each mesh refinement.In Figure 3 we see a slice of the third component of the 3-D fluid velocity u 3 displaying the test problem's non-trivial boundary interaction with the plate.Figure 4 shows the pressure is converging to the solution as well.

Table 2 :
Computed index k in O(h k ) for structure FEM approximations.

Table 3 :
Errors of fluid FEM approximations.The log error ratios approach what is expected for a P 2 /P 1 implementation, namely k = 3, 2, and 2 respectively as shown in Table4.L 2 (fluid) H 1 (fluid) L 2 (pressure)

Table 4 :
Computed index k in O(h k ) for fluid FEM approximations.