LINEARIZED HYDRO-ELASTICITY: A NUMERICAL STUDY

. In view of control and stability theory, a recently obtained linearization around a steady state of a ﬂuid-structure interaction is considered. The linearization was performed with respect to an external forcing term and was derived in an earlier paper via shape optimization techniques. In con- trast to other approaches, like transporting to a ﬁxed reference conﬁguration, or using transpiration techniques, the shape optimization route is most suited to incorporating the geometry of the problem into the analysis. This reﬁned description brings up new terms—missing in the classical coupling of linear Stokes ﬂow and linear elasticity—in the matching of the normal stresses and the velocities on the interface. Later, it was demonstrated that this linear PDE system generates a C 0 semigroup, however, unlike in the standard Stokes-elasticity coupling, the wellposedness result depended on the ﬂuid’s viscosity and the new boundary terms which, among other things, involve the curvature of the interface. Here, we implement a ﬁnite element scheme for approximat- ing solutions of this ﬂuid-elasticity dynamics and numerically investigate the dependence of the discretized model on the “new” terms present therein, in contrast with the classical Stokes-linear elasticity system.

1. Introduction. The interaction of a flexible structure with a gas or fluid flow (either surrounding or contained within the structure) arises in a variety of physical phenomena and, thus, in many applied problems: the stability and response of aircraft wings, vibration of turbine and compressor blades, blood flow in arteries, interocular dynamics, etc. The model belongs in the category of free and/or moving boundary problems. Free boundary problems are typically associated with a setup where the geometry is unknown, but remains stationary and a steady-state regime is considered. A moving boundary scenario instead provides a priori information on the evolution of the geometry. In turn, the free and moving combination requires the unknown boundary to be determined as a function of time and space. In this context, both steady-state and dynamical cases have been the focus of intense them; with the smallness assumption in place the regularity issues still have to be addressed. In [5] the well-posedness of this model is investigated via semigroup theory in the special case when the linearization is performed near a sufficiently smooth low-velocity steady-state regime. The maximality result for the evolution generator was demonstrated using a variational approach and the Babuška-Brezzi theorem following the strategy in [1]. This argument also provides a convenient framework for a numerical study of this system via the finite element method. The Babuška-Brezzi approach eliminates the need for divergence-free fluid basis functions in FEM. Moreover, once the algorithm to solve the maximality system is in place, one can use it to solve the dynamical problem with prescribed Cauchy data by using the classical formula for the flow exponential ( [39,Thm. 8.3,p. 33]) and corresponding implicit schemes.
One of the traits of the linear semigroup flow approach is that the solution lives in the associated "finite-energy" state space associated in this case to elastodynamics and Stokes flow. This regularity is much weaker than employed in the study of fully nonlinear models. For example, [13] takes advantage of five additional orders of differentiability in the elastic displacement initial data (total of H 6 ) on top of the finite energy H 1 of solutions to isotropic elasticity. Of course, the gain partly follows from the assertion that linearization is performed around a smooth solution to the nonlinear problem in the first place, such as linearization near equilibrium. The analysis in [5], however, shows that if the linearization takes into account geometry of the domain and is performed near a nonzero steady regime then existence of finite-energy solutions likely requires a sufficiently large viscosity-this aspect was irrelevant in the simple classical Stokes-elasticity coupling. The analytic study in [5] does not not provide precise quantitative estimates on the relative parameter sizes, and the current paper aims at investigating this aspect further.
1.2. Goals and challenges. The primary goal of this paper is to numerically investigate the linearized model [7] and study the effect of the new, when compared to the classical model, trace terms on the discretization. In particular, we numerically examine: 1. how the ellipticity constants and matrix condition numbers in the associated discretized system are affected by the magnitude of the boundary terms present in the matchings of the normal stress tensors and of the velocities on the common boundary (see (7)-(8) below). 2. whether the size of the boundary terms may affect the convergence rates in the numerical scheme.
The structure of the interface conditions ( (7) and (8) below) is quite involved and the coefficients depend on the exact solution to the associated steady-regime problem around which the linearization is performed. A direct investigation would thus require an exact solution to a nonlinear hydro-elastic system. In addition, elements with curved boundaries would be needed to consider the curvature terms. We can partly circumvent these difficulties by noting that the semigroup-wellposedness analysis [5] critically depends only on the following properties of the linearized system: (i) The norms (in appropriate spaces) of the fluid velocity, the pressure, and of the deformation in the steady regime around which the linearization is performed, and the size of the fluid viscosity relative to these norms.
(ii) The fact that the additional terms (when compared to classical coupling of Stokes and linear elastodynamics) on the interface can be decomposed into zero-order and first-order tangential derivatives only, without normal components. (iii) The regularity and magnitude of the variable coefficients (accounting for the geometry of the common boundary in the interaction) of the terms on the interface.
Based on these observations, we consider a modified system which admits simplified interface conditions and piecewise-flat geometry. At the same time, this model partly recreates the original dynamics by incorporating zero and first-order perturbations on the interface. Thus a large-curvature scenario is simulated by directly adjusting the size of the coefficients of certain trace terms.
The system is numerically implemented using a finite element scheme with Taylor-Hood elements and illustrates optimal convergence rates [4] for the considered analytic test-solutions. The dependence of the numerical system on the size of the parameters specific to this model is then investigated.
The main findings are summarized here and can be found in Sections 6 and 7. The primary observations are that: (i) the system is more sensitive with respect to the lower-order terms on the common boundary Γ c which depend on the curvature of the interface; (ii) large zero-order terms on the boundary and large velocity of the steady regime around which linearization is performed lead to singularities independently of the size of the viscosity (confirming the observations in [5]); (iii) larger viscosity renders the condition numbers of the matrices in the discretization more robust with respect to the parameter change; (iv) within the admissible parameter range, the FEM convergence rates in a test case do not seem to be influenced by the dynamics on the interface and match (slightly better than) those derived in [4]. All these considerations, of course, do not apply to the standard coupling of linear Stokes flow and system of linear elasticity, where the boundary terms arising from the linearization of the common moving boundary are absent.

1.3.
Outline. The rest of the paper proceeds as follows: • in Section 2 we present the original nonlinear model for a steady regime; Subsection 2.2 discusses the linearization of a dynamic problem around this steady regime. • Section 3 summarizes the known semigroup well-posedness result for the linearized system and Subsection 3.1 details the formulation of the associated resolvent problem as a constrained minimization in the context of the Babuška-Brezzi theorem. • In Section 4 we state a simplified system that captures some of the key properties of the original linearized model and is subsequently used for numerical simulation. • Section 5 outlines a few specifics of the finite element implementation of the Babuška-Brezzi system. In Subsection 5.2 we examine errors and convergence rates with respect to the mesh parameter of the FEM scheme for some available closed-form solutions. • Section 6 investigates the properties of the discretized problem with respect to the new parameters-which quantify the difference between the presented linearization and its classical Stokes-elasticity analogue. 2. PDE model. Motivated by the investigation of the stability of an equilibrium in hydro-elasticity systems, the steady-state solution was perturbed by a timedependent distributed forcing term, and the "total linearization" around the steady regime was obtained in [7]. Here we outline only the final results, and refer the interested reader to [7] and [5] for additional details.
The investigation focuses on fluid that is Newtonian, viscous, homogeneous and incompressible. Its behavior is described by its velocity field w and pressure p. The steady state PDE model is given by Let us elaborate on the notation used in (1): • The viscosity of the fluid is ν > 0, and the fluid strain tensor is given by ε(w) = [Dw + (Dw) * ]/2, where Dw is the gradient matrix of w, and (Dw) * represents the transpose of Dw. The fluid state satisfies the Navier-Stokes equations with no-slip condition on the outer boundary Γ f . • The elasticity domain Ω e is described according to a map acting in a fixed reference domain (Lagrangian formulation). Then the evolution of the fluid domain Ω f is induced by the structural deformation through the common interface Γ c . Given a material reference configuration for the solid Ω e ⊂ D with boundary Γ c , let Ω f = D \ Ω e be the reference fluid configuration. Then the volume D is described by a smooth, injective, orientation-preserving map: For any x ∈ Ω e , ϕ(x) represents the position of the material point x.
On Ω f , ϕ(x) is defined as an arbitrary extension of the restriction of ϕ to Γ c , which preserves the boundary Γ f , i.e., ϕ = I Γ f on Γ f . We define J(ϕ) = det(Dϕ) > 0 as the Jacobian of the deformation ϕ. where n is the unit outer normal vector along Γ c with respect to Ω e . Equivalently, the equilibrium equations for the elasticity can be written on the reference configuration Ω e in terms of the Piola transform of the Cauchy stress tensor. This translates into a more complicated matching on the boundary of the reference configuration: Pn = J(ϕ)(σ(p, w) • ϕ)(Dϕ) − * n on Γ c . • We assume that the distributed force v acts on the entire domain D, and therefore v 1 = v| Ω f , and v 2 = v| Ωe .
The well-posedness analysis for steady-state interactions of this form appeared in [24] are recently, in [50]. Accordingly, we assume that he nonlinear system (1) is well-posed, and possesses sufficient regularity to take advantage of the existence result in [5], discussed below, for linearized model.
2.1.1. Notation. Before proceeding further we summarize the basic notation that will be employed henceforth: • By M m×n denote the real m × n matrices, M m := M m×m .
• A..B = Tr(A * · B) ∈ R is the Frobenius product in M 3 .
• (Df (a)) ij = ∂ j f i (a) ∈ M 3 is the Jacobi matrix at a ∈ X of any vector field f = (f i ) : X ⊂ R 3 → R 3 . • If T is a 2-tensor, then DT is a 3-tensor whose contraction with a vector gives For a given vector n, domain Ω ⊂ R k . Thus ∆b Ω = Tr(D 2 b Ω ) is the additive curvature of Γ = ∂Ω. Also note the following bilinear forms, norms and duality pairings: • The bilinear form (·, ·) X will indicate the inner product on a given Hilbert space X. • Product spaces [X] k will be abbreviated with bold X. Likewise the dual [X ] k ∼ = ([X] k ) will be denoted X . The value of k will be clear from the context. • · X will denote the norm on X.

2.2.
Linearization of the free boundary interaction. The system is perturbed in the forcing v, by assuming that the perturbation v s (x, t) depends linearly on a small parameter of variation s ∈ [0, This perturbation arises in optimal control problems. Of particular interest is the case when the state of the entire hydro-elastic elastic system is controlled from the fluid domain or its outer layer. Functions (ϕ s , w s , p s ) solve a non-cylindrical evolution problem (that is, the associated space-time domain is not the cylinder given by a Cartesian product of the form (0, T ) × X, but rather ∪ t∈(0,T ) {t} × X(t) since the space domains, in this case Ω f and Ω e , change with t): where n s (t) is the unit outer normal vector along Γ s c (t) with respect to Ω s e (t), J s (t) is the Jacobian of the deformation ϕ s (t), T s is given by We will use " " for the s-derivatives at 0, i.e., .
Define the elastic displacement vector (for linearized elasticity on the deformed configuration) as U := ϕ • ϕ −1 , and let Then, in particular, it follows that The nonlinear strain associated with ϕ is given by and the linearized strain associated with U becomes Implicitly, we have that Tr(σ (ϕ) • ϕ −1 ) = Tr(Θ * (DU )Θ). For a shorthand define and then rewrite the strain as the symmetrized gradient in this form: The associated stress tensor then takes the following form: For a shorthand let's also define The Cauchy stress tensor can be expressed as: The shape derivative of T w.r.t. s at s = 0, calculated in [7], is given by To formulate the boundary conditions it is convenient to use the following modified stress tensor instead of T : With the above notation at hand, we are ready to restate the existence theorem.  (2) is s-differentiable and that the initial condition (ϕ 0 , w 0 , p 0 ) associated with (2) solves the steady regime problem (1). Then (ϕ , w , p ) satisfy the following cylindrical (i.e. on a time-independent space domain) evolution problem: where the boundary operator B is given by (Recall that b Ω is the oriented distance function for domain Ω e ; in particular, D 2 b Ω carries information on the curvature of Γ c .) 3. Wellposedness result. The well-posedness of (7) was studied in [5] using semigroup theory methods. For a convenient functional framework, equip the space L 2 ( Ω e ) with an equivalent inner product (assuming det Θ is positive bounded below on Ω) For the space of elastic displacements H 1 ( Ω e ), we define the following equivalent (under appropriate assumptions) inner product Lastly, we set the fluid velocity space to The above amounts to the definition of the state space for the entire coupled system: where w is an extension of the steady "linearization velocity" w to the entire control volume domain D. Under the assumption that w is small and regular enough to admit such an extension with small Dw L ∞ (D) , this bilinear form defines an equivalent inner product on H .
In the above functional setting, system (7) can be recast as an evolution problem y = Ay on the state space H , with the generator action A : The mapping π is a solution to the elliptic problem which expresses the pressure p via w , U [5, (5.7), p. 1976], and the elasticity operator S stands for (For convenience in [5] a bounded perturbation-the identity-is added to S to ensure the ellipticity of the associated operator on H 1 ( Ω e ); but a bounded perturbation does not affect semigroup wellposedness).
A complete description can be found in [5, Sec. 6.3, pp. 1986-87]. Ultimately the definition ensures that sufficiently regular states satisfy the appropriate trace conditions on the interface Γ.
The solution to the steady regime problem (1) is assumed sufficiently smooth (see Theorem 2.1), as it appears in the coefficients of the linearized system. More precise sufficient conditions are given in [5, Assumption 3(a), p. 1972]. In turn, sufficient smallness assumptions on the steady regime around which the linearization is performed are more carefuly quantified in [5, Assumption 3(b), p. 1973]. In this setting, when the linearization state is small and regular, it is demonstrated in [5] that a bounded perturbation of A yields a maximal dissipative operator, whence: . Under suitable regularity and smallness assumptions, system (7) generates a C 0 semigroup on the state space H . In addition, the unique weak solution (U , U t , w ) possesses additional regularity from which the velocity matching boundary condition on the interface also gives U t ∈ L 2 (0, T ; H 1/2 ( Γ c )).

3.1.
Babuška-Brezzi formulation. The proof of Theorem 3.1 in [5] verifies that A is an ω-dissipative operator on H modulo a bounded perturbation. The wellposedness analysis reduces to the study of a bounded perturbation of A with parameter Λ, for some y = [U , V , w ] * ∈ D(A) and given y 1 = [U 1 , V 1 , w 1 ] * ∈ H . The resolvent problem essentially determines the evolution generator for the original system (e.g. [39, Thm. 8.3, p. 33]), so the time-dependent case can subsequently be studied. One can follow the approach of [1] to reformulate (15) into a Babuška-Brezzi system Define the extension map and let a S be a bilinear form associated with the elastic operator S in (14): Γc , Then the bilinear forms a λ , b and the functional F in (16) are given by 4. A prototype problem. The theoretical analysis presented in [5] confirms the well-posedness of the linearized model under certain assumptions on the size of the viscosity and the magnitude of the new terms that arise in (7), but do not show up in the classical Stokes-elasticity coupling. The theoretical framework so far, however, offers little insight into the exact quantitative relation between the parameters of the system. At the moment there is no exact data, whether theoretical or numerical, available on what exactly "sufficiently small" and "sufficiently large" mean when investigating the well-posedness of this system, other than the fact that such parameter ranges exist. The "small" terms in question ultimately pertain to the curvature of the contact interface as well as the velocity of the steady regime around which the linearization is performed . The "large" parameters are primarily the viscosity of the fluid and the (real part of) the resolvent value λ. These features need to be addressed in further studies, but to lay numerical groundwork for this analysis we attempt to look into the following features: 1. Determine whether a sufficiently high viscosity can counteract the influence of the terms specific to model (7): (a) Zero-order terms on the interface (including to curvature) (b) First-order terms on the boundary (c) Size of the state around which linearization is performed. 2. Check the need for the additional parameter Λ used in (15). It does not affect semigroup generation since it represents merely a bounded perturbation. However, it was instrumental in proving the solvability of the Babushka-Brezzi variational formulation in [5]. Since numerical semigroup approximations will require Λ = 0 and would only vary λ, this setting will be attempted here numerically. 3. Verify if the system is necessarily ω-dissipative. In particular, check, at least numerically, what happens as λ → 0 + relative to the viscosity. As evident in [5], the parabolic smoothing from the viscous term was a major driver of the well-posedness argument, and the magnitude of the viscosity may be correlated with the lowest admissible bound on the resolvent parameter λ.
A direct simulation of (7) would in general involve a domain with non-constant boundary curvature (hence curvilinear elements) and an exact solution to the steady regime problem introduced in (1). However, some features of the system can be investigated without resorting to such a costly implementation. The argument in [5] demonstrated that the precise structure of the coupling in (7) is not critical to the semigroup well-posedness; however, the following aspects of the model are important: • The fact that the components of the boundary operator B(V ) defined in (8) are (variable) linear combinations of V and of the first-order tangential derivatives of V (that is not obvious from (8); see [5, Prop. 6.2, p. 1979]). • The magnitude of the coefficients that perturb the elastic stress tensor in the interior, in particular T in the second summand of (6), and the size of the coefficients of the zero and first-order terms of V in (8). In addition, the magnitude of these coefficients relative to the viscosity ν plays a crucial role. • The norm of the differential Dw, of the velocity in the linearization state, appearing in the velocity matching condition in (7).
The semigroup wellposedness relies exclusively on the magnitudes and smoothness of the above terms. For instance, the nature of w being a steady regime solution played no role in the proof. Accordingly, the prototype problem will be designed to capture only these essential features of the model.
We consider a simplified two-dimensional model on a rectangular domain: where D Γc V (x) = (DV )τ (x) with an oriented tangent field τ on ∂ Ω e . In particular, the following adjustments have been introduced in (21): 1. Since the focus is on the normal stress and velocity matching conditions on the interface, the stress in the interior of the elastic component is modeled by the standard linear isotropic tensor. 2. We define w = ε 3 w to replace the velocity of the steady state solution. Its norms are controlled by parameter ε 3 . Finding an admissible w would require to solve a nonlinear steady-regime free boundary problem. However, since the semigroup investigations turn out to not depend on the precise structure of this term, we will not require w to be such a solution. Instead, w will be chosen an arbitrary function, though we will consider cases when it is divergence-free. 3. The effect of the curvature and of the first-order tangential derivatives in (8) is simulated by redefining the boundary operator B (8) to incorporate such terms explicitly. The (positive) sign of the zero-order term is chosen so that it does not contribute to the dissipativity of the uncoupled elastic problem: for some c(x). We will test it with c(x) = 1 and c(x) = |x| 2 . With these simplifications, the maximality system (15) becomes and it is solved using the Babuška-Brezzi system (16)-(20) where 5. FEM framework. A very successful and popular finite element framework for Stokes flow is based on the Taylor-Hood scheme, developed in [48] and subsequently studied in [4,49]. The model uses mixed {P 2 , P 1 } finite element basis functions for the fluid and pressure spaces respectively to approximate solutions of the system (16). The solvability of the discretized model ultimately hinges on the inf-sup condition for the associated bilinear "constraint" form b in (16). Moreover, the convergence rates of the scheme require the associated inf-sup inequality to hold uniformly with respect to all sufficiently small values of the mesh parameter. The extension of this particular FEM model to fluid-structure interactions was first carried out in [1], albeit without analysis of the discrete inf-sup condition. This latter aspect differs from its Stokes counterpart because the lack of the non-slip boundary conditions on the interface between the fluid and the solid enlarges the underlying function spaces. Specifically, the pressure function is determined uniquely instead of up to a constant and the fluid state space no longer has any constraints on the interface boundary. The rigorous well-posedness analysis for the FEM approximation of Stokeselasticity coupling was investigated in [2]. It has been shown that the FEM fluidstructure model is well-posed, numerically stable, and, moreover, if the elliptic bilinear forms are evaluated exactly on the FEM subspaces, then the scheme converges at the classical optimal rates of [4]. In addition, the analysis relies only on the ellipticity of the associated bilinear form a λ , as in (18). Thus, the presence of tangential boundary terms on the fluid-solid interface, as those in (22), does not affect the well-posedness of the FEM scheme, provided the ellipticity of a λ is preserved.

LORENA BOCIU, STEVEN DEROCHERS AND DANIEL TOUNDYKOV
The study in [2] still leaves some questions unanswered. In general, only the stability of the scheme is guaranteed because the elliptic form (18) is itself approximated using an FEM solver on the solid domain. Likewise, the existing studies can shed no light on the dependence between the size of the boundary perturbations and the viscosity parameter. Here we attempt to pursue this question numerically.
The approximation of function spaces is based on a quasi-uniform family of triangular meshes. The "first-level" mesh considered for the discretization is shown on Figure 2  All subsequent mesh levels are obtained recursively by dividing each current element into four triangles by connecting the edge midpoints with line segments. In the {P 2 , P 1 } Taylor-Hood scheme, the fluid velocity w is interpolated by a piecewisequadratic basis, and pressure p is approximated by a piecewise-linear function. Quadratic basis elements use the vertices as well as the edge midpoints for the nodes, while the linear ones only use the vertices of the nodes.
Let R denote the reference element with vertices (0, 0), (1, 0), and (0, 1) with standard reference shape functions ϕ k ref , k = 1, 2, . . . 12 (two sets of 6: one set for each coordinate of the 2D fluid velocity or the elastic displacement vectors). The affine map L K maps R onto the mesh element K (assume it preserves CCW orientation of the local node numbers). Therefore the local basis functions are given by ϕ k The L 2 and energy products are as usual obtained by Gaussian quadrature is used to evaluate these exactly. The requisite boundary integrals are computed by evaluating similar forms through parameterizations of the element edges in the boundary of the domain mesh.

Details of the implementation.
We give a brief overview of the FEM implementation. Consider global piecewise-quadratic basis {ϕ m } (piecewise defined through the local ϕ k K 's introduced above) and piecewise-linear basis functions η m over the mesh. Then the Babuška-Brezzi system (16) translates into the following linear system: The approximations of the fluid velocity and pressure are given respectively by where N 1 is the total number of nodes, including the edge midpoints, and N 2 is the number of vertex nodes.
5.1.1. The stiffness and mass matrices. First, we need to solve the Dirichlet extension problems D S,λ 2 [F, R] described in (17). The size of the corresponding finite element system is not equal to the size of the entire system since here we are only using the vertex and midpoint nodes in the solid domain. This means that initially, in the algorithm, each node should be identified with both a global node number over the entire domain D as well as a global node number for the elastic sub-problem on Ω e . The stiffness matrix for the extension problem is given by 5.1.2. Assemble the system. Once the approximations to solutions of the extension problems have been obtained, we build the matrices for A i,j = a λ (ϕ i , ϕ j ), B i,j = b(ϕ i , η j ), and "load" vector F i = F (ϕ i ) for the Babuška-Brezzi problem. The last two terms terms arising in the bilinear form (18), namely can now be calculated as u * Sv where u is the coefficient vector The same applies to the analogous terms in (20): . Finally, to solve for the coefficient vectors for w and π , we solve the linear system (25); again, however, we first replace the appropriate rows in A B B * 0 with the rows of the identity matrix and the corresponding entries of the right-hand side vector with 0, ensuring that the no-slip conditions are enforced. Let V 1 = 0 and w 1 = 0. Consider a solution where the pressure is constant, say π = 1, and the fluid is at rest w = 0. Then the elastic deformation U must solve the following system (recall that n is the unit normal field outward to the solid domain): This U can be approximated numerically using FEM on the elastic domain. For V 1 = 0 in (23) if the resolvent parameter λ is set to 1, then U and U 1 must coincide. Hence we set U 1 to U and solve, now the full the system instead of just FEM for elasticity, expecting to recover the intended true solution w = 0 and π = 1.
With viscosity ν = 1 at the very first mesh level we recover "machine epsilon" in the error:

Mesh Lvl
No. Fluid Elems Note that this is a purely numerical check without an analytic solution at hand.

5.2.
Convergence rates. The next stage is to test the implementation on some closed-form solutions. The structure of the elastic stress tensor and the presence of the tangential terms on the boundary make it very difficult to find an exact analytic solution compatible with this simplified geometry. The first option is to consider the case when there is no interaction on the interface at all. The true solution in that case simply decouples the fluid and solid components. But the numerical implementation, however, will pick up some error data on the interface. and consider a deformation whose values and first-order derivatives vanish on the boundary: (The 8-th degree polynomial is pointwise very small, hence the 10e6 factor is introduced). Then directly compute We expect to get approximations close to w = 0 and π = 0. Running the code for λ = 2, ν = 1, and ε 1 = ε 2 = ε 3 = 0.1 we obtain the following error estimates The decay exponent for the error in the table is displayed in Figure 3. In this particular case the convergence rate for H 1 ( Ω f ) norm is slightly better than the predicted rate of 2 for Taylor-Hood elements [4,Thm. 1,p.219]. Note that the scheme also shows improved convergence rates for the L 2 ( Ω f ) norm of the solution despite the fact that the fluid domain is polygonal and non-convex [4,Prop. 3,p. 222].

5.2.2.
Test solution 2. Next, consider a closed-form solution for which the 1st-order elasticity traces do not entirely vanish on the interface. In (23), (24) again let along with the following right hand side data: Then for U = λU 1 and V = 0, we see that w = 0 and constant pressure π = ε 2 satisfy the interface boundary conditions. For λ = 2, ν = 1, ε 1 = 0 and ε 2 = ε 3 = 0.1 we outright recover the "machine epsilon", which is not surprising since the solution consists of affine functions only.

Mesh Lvl
No. Fluid Elem 6. Numerical sensitivity analysis. Now we address some aspects of the sensitivity of the numerical approximation to (23), (24) with respect to parameters ε 1 , ε 2 , and ε 3 . Recall that these appear in the "steady regime" fluid velocity and in the boundary operator or its version with a variable coefficient at the first-order term: 6.1. Dependence on ε i : Constant coefficient case. Figures 4, 5 and 6 illustrate the properties of the system with respect to the aforementioned parameters: • ε 1 : the coefficient of the zero-order elasticity term in (26) on the interface. This term simulates the presence of the zero order terms from (8). Among them is the curvature of the boundary as it appears in (D 2 b Ω )V in (8) • ε 2 : the coefficient of the first-order tangential terms in (26). The tangential terms themselves correspond to the quantity div(V )T · n − T · (DV ) * · n in (8). As shown in [5,Prop. 6.2] this expression depends on the tangential derivatives only. • ε 3 : the coefficient of the first-order interior term in the fluid equation in (21).
It governs the magnitude of the state w which simulates the state w around which the linearization is performed in the original system (7).
(In the classical Stokes-elastodynamics coupling all ε i would be zero). The approximating mesh level is 2, that is, one refinement by subdivision of the mesh presented in Figure 2. As a result there are total of 128 fluid elements and 16 solid elements. The figures occupy the subsequent three pages. Each figure showcases multiple graphs and within each graph it considers several viscosity values. The plotted data displays: • The minimal eigenvalue of A + A * in (25), which quantifies the ellipticity of the discretization of the bilinear form a λ in (18). • The logarithm of the condition number for the full matrix in (25) and the relative change in the condition number with respect to the indicated parameter ε i . • The logarithm of the condition number for the elliptic part A in (25) and the relative change in this condition number with respect to ε i .
Based on these graphs the following observations can be made: 1. Large values of coefficients ε 1 and ε 3 both lead to the loss of ellipticity in the discretized bilinear form a λ of the Babuška-Brezzi system. 2. Tangential derivatives with constant coefficient ε 2 have no quantifiable impact on the solvability of the system, at least in flat geometry of our test problem. The variable coefficient case is considered below. 3. Larger values of viscosity ν yield an increase in the condition numbers of the matrices (the full matrix and sub-matrix A) in (25). However, the system becomes more robust with respect to the interface dynamics. Thus, relative change of the condition numbers with respect to ε i are negligible at higher viscosity values. 4. Viscosity does extend the range of values for ε 1 and ε 3 under which the discretized form a λ remains elliptic. However, it can do so only up to a certain extent. For example, in Figure 4 changing the viscosity from 0.1 to 1 notably extends the range of admissible values for ε 1 before the eigenvalue becomes negative. But there is relatively no gain out of the change of viscosity from 1 to 10.
6.2. Dependence on ε 2 : Variable coefficient case. The interesting aspect of the system showing up on Figure 5 is the independence of the ellipticity of the discretized bilinear form a λ of the value of ε 2 , which corresponds to the 1st-order tangential terms in the boundary condition (26). This behavior may be caused by considering the constant-coefficient tangential derivatives only. In this run we repeat the data set for parameter ε 2 , however, we use the model (27) with a variable coefficient next to the 1st-order terms, which would be more faithful to the original (8).
As evident from Figure 7, the variable coefficient plays a key role and the eigenvalues degenerate into negative as ε 2 grows. As before, increasing the viscosity beyond a certain threshold does not seem to have much effect, since the ellipticity still diminishes about the same value for ε 2 .
6.3. Dependence on viscosity ν and resolvent value λ. Figures 4, 6 and 7 suggest that that for a fixed resolvent value λ, even large viscosity cannot enable us to use arbitrarily large perturbations. That is no matter how big the value of ν is, this particular numerical approximation cannot accommodate too big ε i . This behavior is more apparent on the following Figure 8 where the zero eigenvalues of A + A * trace out an asymptote at ε 1 ≈ 0.83.
Larger resolvent values λ > 0, on the other hand, may extend the admissible range of the parameters ε i , as shown on the next Figure 9. There, for a fixed viscosity, the intersection of the minimal eigenvalue graph of A + A * with the zero plane takes place at higher values of ε 1 , as λ increases. The numerical data for the same figure also shows a "spike" near zero, since when ε 1 = 0 we are back to the classical model which is dissipative ((0, ∞) is in the resolvent set).
On the other hand, smaller values of the viscosity, e.g., ν < 0.05 with λ = 0.01 are not admissible at all because matrix A cannot be formed, since the harmonic extension problem becomes ill-posed. As Figure 9 shows, the value λ = 0.01 would require a much higher viscosity ν, e.g., equal to 1. However, testing separately the case, ν = 1, ε 1 = 0.01 and λ = 0.001 on the 2nd-level mesh likewise yields an illposed elliptic problem for the elastic component. This data suggests a discontinuity near the origin of the graph on Figure 9 which would be in agreement with the asserted strict ω-dissipativity (ω > 0) of the evolution generator: namely, for any ε 1 > 0 it is likely not possible to admit arbitrarily small λ > 0.

7.
Conclusions. This note provides a preliminary numerical study of the solvability and the main features of the solution to a linearized fluid-elasticity interaction around a steady-state regime. The linearization is performed with respect to an external forcing and takes into account the shape and movement of the boundary interface between the solid and the fluid. The investigated test problem uses a simplified domain, but following the theoretical analysis in [5] it includes special terms that simulate the key aspects of the original model.
The numerical experiments underline the following traits: • The system degenerates with respect to zero-order terms of elasticity in (8), and with respect to the norm of the state around which the linearization is performed.
• The dependence on the oblique derivative condition (8) may not be apparent in polygonal geometry unless variable-coefficients are used. Thus, constant coefficient approximations of (8) may turn out highly inaccurate. • Large enough positive value of viscosity ν is needed to ensure solvability of this finite element model. However, when solving the resolvent problem for a particular λ, the admissible size of the perturbations is limited independently of the viscosity. The additional bounded perturbation parameter Λw, for Λ > 0, in the fluid equation may not be necessary to solve the variational problem.
• The system appears to be strictly ω-dissipative independently of the parabolic effect. That is, even for small perturbation values and relatively large viscosity the spectrum of the generator may overlap a neighborhood of the origin. These observations serve as tentative guidelines for further numerical and theoretical studies of the original model (7). For example, the dependence on the coefficient of the zero-order terms on the common boundary, confirms numerically that the geometry of the interface in the original model plays a critical role in the dynamics of the phenomenon. The additional dissipative term "+Λw" added to the fluid equation in order to assure the maximality of the generator in [5], may not be necessary to solve the resolvent problem in this variational formulation. Last, but not least, the numerical evidence points towards strict ω-dissipativity, thus suggesting that the parabolic smoothing from the fluid cannot ensure the non-expansive property of the evolution, even for viscous flows and small perturbation values. Thus, stability studies of the linear model should consider additional stabilizing feedbacks.