SENSITIVITY ANALYSIS FOR A FREE BOUNDARY FLUID-ELASTICITY INTERACTION

. In this paper a total linearization is derived for the free boundary nonlinear elasticity - incompressible ﬂuid interaction. The equations and the free boundary are linearized together and the new linearization turns out to be diﬀerent from the usual coupling of classical linear models. New extra terms are present on the common interface, some of them involving the boundary curvatures. These terms play an important role in the ﬁnal linearized system and can not be neglected.


1.
Introduction. Interaction between a fluid and a structure via a common interface is a basic coupling mechanism in continuum mechanics. There are basically two different scenarios: the first one involves an elastic solid which is fully immersed in a fluid (like a submarine under water, or adherence/detachment of leukocytes in blood flow), and the second one is represented by a fluid flowing inside a body, like the blood flowing in a human artery. Both models have attracted a lot of attention over the years, and mathematical and numerical subtleties are common to both cases.
Historically, the problem of rigid body motion in both compressible and incompressible flows was first considered. This entailed coupling the Navier-Stokes equations with systems of ordinary differential equations [15,22,25,27,28,37]. More recently, there have been works that treat the motion of an elastic body in an incompressible flow, by coupling the Navier-Stokes equation with hyperbolic elasticity equations on fixed domains: for linear models [1,23], and for nonlinear models [2,3,31,32,35]. Other types of structures have been studied as well, such as plates and beams, in addition to strongly damped elastic bodies [12,21,29].
In most real-life problems of interaction between a fluid and an elastic structure, the appropriate case is that of a "free boundary", i.e. a moving fluid-solid interface [23], which is generally not known, and must be found as part of the solution process. Mathematically, the interaction is described by a partial differential equations (PDEs) system that couples the parabolic and hyperbolic phases, where the key issue is that the traces of the elastic (wave/hyperbolic) component at the energy level are not defined via the standard trace theory. The loss of regularity induced by the hyperbolic component left the basic question of existence open until recently. [16,17,31,32] (for the case of a linear or quasi-linear elastic body flowing within a fluid). The local existence and uniqueness of solution were obtained by analyzing the fluid part in a hyperbolic-type functional framework that, at the price of increased compatibility assessments between the initial and boundary conditions, allows the increased regularization of the elastic behavior to be discarded. To our best knowledge, these existence results have not yet been extended to the configuration of a fluid inside an elastic envelope, which is of great interest in the modeling and analysis of the cardiovascular system. This paper will focus on performing a sensitivity analysis for the dynamical case of fluid moving inside a compliant elastic body. A specific example of this model is represented by blood flow in human arteries. However, the theory and the techniques also apply to the counterpart case of an elastic body floating within a fluid. When linearizing any coupling of fluid and structure, such as those occurring in arteries, submarines, and bridges, a major challenge arises in just how to linearize the free boundary. Linear models for the coupling present in the literature were obtained by linearizing each equation separately, then coupling the two linear problems. However, such an approach misses important coupling terms on the boundary, which would be present if the system were linearized as a whole.
Our approach in this paper is to couple the systems (Navier-Stokes and nonlinear elasticity) first, yielding a complete non-linear model, and then linearize the system consistently. This is necessary due to the large-deformation regime considered. The method is much more challenging mathematically, resulting in a free boundary which moves with the given perturbation parameter. The clear advantage of this idea is that it highlights the geometrical aspects of the problem, like curvature and boundary acceleration, which are unaccounted for in current linear models. This is particularly important in the case where the boundary oscillates, sending the mean curvature to infinity. Modeling of this geometrical aspect is critical for a correct physical interpretation of the fluid-structure interaction [9,10].
We successfully addressed this problem in the case of steady-state fluid-elasticity interactions in [9,10]. We adapted recent results developed in shape optimization and image analysis [11,18,24,38,41] in order to fully characterize the linearization of the steady-state Navier-Stokes/potential fluid-nonlinear elasticity coupling [9], and then obtain a good setting for the linearized problem where a variational approach can be used to prove existence and uniqueness of solution [10]. Even though the system is not dynamical in this case and thus the common boundary does not move with time, the problem still belongs to the category of "free boundary problems", due to the fact that the stationary interface is unknown and has to be determined by the steady-state problem. The linearized model around fluid rest obtained for the case of steady-state interaction and zero applied body forces [9] shows the role of the free common interface in the coupling. While the linearization preserves the linear equations for fluid and elasticity in their domains, it is quite different on the boundary from the usual coupling of linear elasticity and linear fluid.
Our method highlights the presence of the curvature terms on the common interface, captured in the form of the Hessian D 2 b Ω , where b Ω is the oriented distance function [9].
The goal of this paper is to provide the first linear model for fluid-elasticity interactions that takes into account the common interface and its curvatures, which are critical for a correct physical interpretation of the coupled system. Specific areas of application are represented by (i) the study and simulation of blood flow in compliant (large and medium sized) vessels, (ii) stability analysis performed on linear systems, or optimal control problems, such as Bolza control problems [34], (iii) control problems in fluid-structure interactions, and (iv) numerical analysis for fluid-elasticity interactions.
The rest of the paper is organized as follows. Section 2 will introduce the notation used throughout the paper, the problem and the PDE model. Section 3 will present the main results and discuss their significance. Lastly, Section 4 will be devoted to the proofs of the theorems introduced in Section 3.

2.
The problem and the model.

2.1.
Notation. For the rest of the paper, we use the repeated index convention for summation whenever the same Latin index appears twice, and the following notation: • Div T (a) = ∂ j T ij e i ∈ R 3 is the divergence of any second-order tensor field T = (T ij ) : X ⊂ R 3 → M 3 at a ∈ X.
Since we ignore the distinction between covariant and contravariant components, we will identify the set of all second-order tensors with the set M 3 of all square matrices of order three.
, ∀x ∈ R n is the oriented distance function from x to Ω, for any Ω ⊂ R n . • H = ∆b Ω = Tr(D 2 b Ω ) is the additive curvature of Γ = ∂Ω.  We consider the following PDE model for the fluid-structure interaction: where n(t) is the unit outer normal vector along Γ(t) with respect to Ω(t).
• The investigation is focused on flow in large and medium sized vessels, suitably modeled in the Newtonian framework. Hence the flow is governed by the Navier-Stokes equations in terms of velocity of the fluid w and its pressure p.
The fluid occupies domain Ω c (t) = D \Ω(t) (called the lumen of the vessel), which is moving with time, and has piecewise smooth boundary Γ(t) ∪ Γ in ∪ Γ out . The viscosity of the fluid is ν > 0, and the fluid strain tensor is given by ε(w s ) = [Dw + (Dw) * ]/2 (where Dw is the gradient matrix of w, and (Dw) * represents the transpose of Dw).
• In a fluid-structure interaction framework the evolution of the fluid domain Ω c (t) is induced by the structural deformation through the common interface Γ(t). Thus we describe Ω c (t) according to a map acting in a fixed reference domain. This approach is usually used for the solid domain Ω(t), by means of the Lagrangian formulation. Given a material reference configuration for the solid O ⊂ D with boundary S ∪ Γ , we take O c = D \Ō as the reference fluid configuration . Then the control volume D is described by a smooth, injective map: This means that ϕ(x, t) = ϕ e (x, t), ∀x ∈Ō, and ϕ(x, t) = Ext(ϕ e | S )(x, t), ∀x ∈Ō c .
Let J(ϕ) define the Jacobian of the deformation ϕ(t). We take ϕ to be orientation-preserving, which implies that We consider the St.Venant -Kirchhoff equations, which model large displacement, small deformation elasticity. The 3D problem is challenging from the mathematical point of view, and interesting from the physical point of view, because it cannot be reduced to its boundary, like in the case of a membrane or a shell.
The elasticity evolution is given by its motion ϕ and the second Piola-Kirchoff stress tensor, related to ϕ through the following constitutive law: where λ > 0 and µ > 0 are the Lamé constants of the material, and σ(ϕ) is the Green-St Venant nonlinear strain tensor, given by Remark: In the case of linear elasticity, the strain tensor reduces to where u = ϕ − I.
In general, the equilibrium equations for elasticity are written on the reference configuration O, in terms of the Lagrange variable x, using the Piola transform of the Cauchy stress tensor In our case, since the coupling takes place on Γ(t), which is the boundary of the deformed configuration, we wrote all the equations in the system on the moving domains Ω(t) and Ω c (t). In order to provide the accurate matching of the two dynamics on the common interface, we need the elasticity equations in terms of the Cauchy stress tensor T . Using the connection between P and T provided in [14], we obtained the following formula for the Cauchy stress tensor T [9]: T :Ω → S 3 , • The interaction takes place on the common boundary Γ(t) and is realized via suitable transmission boundary conditions by requiring continuity of both fluid and boundary velocities, and of the normal stress tensors across the interface Γ(t). A velocity profile is prescribed at the proximal boundary Γ in , i.e. the Dirichlet datum is given by a smooth, compactly supported function The compatibility condition with respect to the divergence free property of the fluid, and the matching of the fluid velocity and the speed of the common interface Γ(t), give the following identity: where n c is the unit outer normal vector along Γ(t) w.r.t Ω c (t).

Assumption 1.
• Initial Data: We assume that the initial conditions have the following regularity: Remark 2.1. In [9,10] we considered the steady regime associated with F = 0. In the present case, F (x, t) can be seen as a random variable (introduced as a new step in the modeling process) that forces the regime to be time dependent.
3. Perturbed model and linearization technique. The system is perturbed through the boundary flux c in (x, t), by assuming that c in (s, x, t) depends on a small parameter of variation s ∈ [0, s 0 ]. We will assume that the dependance on s is linear, and take c in (s, . As a consequence, all the unknowns become functions of s (w s , p s , ϕ s , T s ), and moreover, the geometrical domains are perturbed through s as well (Ω s , Γ s , and Ω c s ). This produces a free boundary that moves with the given parameter s and thus makes the problem challenging mathematically. Therefore system (1) becomes: where n s (t) is the unit outer normal vector along Γ s (t) with respect to Ω s (t), J s (t) is the Jacobian of the deformation ϕ s (t), i.e. J s (t) = det(Dϕ s ), and T s is given by: We associate the following assumptions with system (3): We assume the following existence of weak derivative w.r.t. s: and we denote by ϕ ,ẇ andṗ those derivatives at s = 0.
3.1. The s derivatives of w s , p s , and ϕ s . We assume that system (3) has unique solution (ϕ s , w s , p s ) for any s in some small interval [0, s 0 ], having the same regularity as the solution to (1) mentioned in Assumption 1.
Our goal is to differentiate system (3) with respect to s at s = 0 and obtain a linearized system for the shape derivatives of w s , p s , and ϕ s (which become the new unknown variables). This will represent the total linearization of (1) around a non-steady regime.
We will use the usual " " notation for the s derivatives. Thus , and ϕ = ∂ ∂s ϕ s s=0 stand for the derivatives of w, p, and ϕ, w.r.t. s, at s = 0. For each time t, the perturbation Γ s (t) of the boundary Γ(t) is built by the flow of the vector field When computing the shape derivatives, the vector field associated with the flow T s will appear only at s = 0, in the form of From [18], we have the following formulas, which will be useful in the next sections: For each value of the parameter s, the moving boundary t → Γ s (t) is built by the flow of the transverse speed vector field Using the flow map T s (V t ), we can transport the unknowns in the system on fixed domains (with respect to s), and then compute the material derivatives as usual: Remark 3.1. Note that the s derivatives introduced above are in fact shape derivatives with respect to V t , which is a vector field that depends on ϕ s and thus it is not given a-priori. We perturb the geometry through the parameter s, hence the geometry moves with the flow of a vector field that depends on ϕ s . This is different from the standard theory on shape derivatives, where the domain is perturbed by an a-priori given vector field V and then the speed method is applied. Nevertheless, the s derivatives are pseudo shape derivatives, and we will use the two names interchangeably from now on.
Shape derivatives are then obtained using the well known relation between material and shape derivatives: 3.2. The shape derivative of T s . At this point, we need to compute the shape derivative for the Cauchy stress tensor T s . The aim is to express the elastic equation in term of the new variable U (t) = ϕ (t) • ϕ −1 (t). Using the above expression (4) for T s and basic rules of differentiation, we obtain where the shape derivative of P is given by with Combining (10) with (6), we obtain the shape derivative for J s (t): Now using (11) back into formula (7), we obtain the formula for the shape derivative of T s • ϕ s : On the other hand, direct differentiation of T s • ϕ s w.r.t. s leads to where (DT ) is a three entries tensor representing the gradient matrix of T . Its contraction with ϕ gives the second order tensor [(DT ) • ϕ] ϕ .
Solving for T in (13) and using (12), we obtain (14) above provides the relationship between T and P . If one wants to see the clear dependance of T on U , this can be obtained as follows.
For sake of exposition, we introduce the following notation. Let Now combining (12) with (14), and using the newly introduced variable Θ, we obtain

LORENA BOCIU AND JEAN-PAUL ZOLÉSIO
In order to simplify the preceeding expression, we will use the following lemma: Proof. The proof is immediate. Starting on the right hand side, we have: As an immediate consequence of Lemma 3.1, we obtain the following relations: For sake of exposition, we let Using the notation and the relations introduced above in (15), we obtain our final formula for T :

4.
Main results. Below we present the main results in this paper. The first theorem describes the total linearization of the free boundary Navier-Stokes fluidelasticity coupled system, around an arbitrary smooth solution (ϕ, w, p) of (1). , and p = ∂ ∂s p s (t) s=0 satisfy the following linear system: T and T are given by (16) and (2), respectively, and (ϕ, w, p) is a smooth solution for (1).
Remark 4.1. Note that the linearization obtained in Theorem (4.1) is more complicated than just the coupling of the associated linear problems in variables (ϕ , w , p ), and it shows the contribution and influence of the common interface Γ (through terms like (T − pI + 2νε(w)) · [(D Γ U ) * n + (D 2 b Ω )U Γ ]) in the linearized coupled system. Remark 4.2. In the linearized system (17) we wrote the the equation of elastodynamics on Ω(t), in terms of the shape derivative of the Cauchy stress tensor T . In this form, one can see exactly the new terms (depending on the variable U ) that are present in the equation.
We would like to point out that the linearized elasticity equation is equivalent to the following one on the reference configuration O: Moreover, if the exact dependence on T is not needed, the linearized elasticity equation on Ω(t) can be simplified to the following expression: 4.1. The specific situation when (ϕ, w, p, Γ) solves the steady regime associated with (1). Let (ϕ, w, p, Γ) be the unique smooth solution of the steady Navier-Stokes and nonlinear elasticity coupled system: where n is the unit outer normal vector along Γ with respect to Ω, F (x) ∈ L 2 (D), and c in (x) ∈ H 2 (Γ in ) ∩ H 1 0 (Γ in ). Assume that (18) . Then (ϕ s , w s , p s ) solve the following non-cylindrical evolution problem: where n s (t) is the unit outer normal vector along Γ s (t) with respect to Ω s (t), J s (t) is the Jacobian of the deformation ϕ s (t), and T s is given by (4).  , and p = ∂ ∂s p s (t) s=0 satisfy the following cylindrical evolution problem: where U = ϕ •ϕ −1 , Θ = Dϕ•ϕ −1 , T and T are given by (16) and (2), respectively, and (ϕ, w, p) is a smooth solution for (18).  (20) is equivalent to the following form: The following Corollary is an immediate consequence of Theorem (4.2), where the system is linearized around fluid rest, i.e. w = 0.
Remark 4.5. Note that (20) and (21) are simplified versions of (17). Nevertheless, the terms involving the curvatures of the boundary and the boundary acceleration are still present in the two linearizations, proving again that the common interface can not be neglected while performing sensitivity analysis on the nonlinear coupled system (1).

5.
Proof of Theorem (4.1). Let I T = (0, T ) be the time interval. Recall the perturbed system: where n s (t) is the unit outer normal vector along Γ s (t) with respect to Ω s (t), J s is the Jacobian of the deformation ϕ s (t), i.e. J s (t) = det(Dϕ s ), and T s is given by:

The linearized sticking condition on Γ(t). The sticking condition [∀x
can be written in a weak form as follows: Using the change of variable T s (V t ), (23) becomes where the density ω s is given as Note that ω 0 = 1 and from [18] we know that where H(t) = Tr(D 2 b Ω(t) ) is the additive curvature of Γ(t), and b Ω(t) is the oriented distance function associated with Ω(t) [18].
In (24) we take derivative w.r.t. s at s = 0 and obtain: Since (26) is true ∀Φ ∈ C 1 (D), then it is equivalent to the following boundary condition on Γ(t): Combining (27) with the sticking condition on Γ(t) present in the coupled system (1), we obtain the new, final boundary condition on Γ(t):
Therefore we obtain the following linearized Navier-Stokes equation on Ω(t) in variables w and p :

5.3.
Elasticity derivative. First, we write the equation of elastodynamics on Ω s (t) in the following equivalent form: This means that ∀R ∈ C 1 (I × D), we have τ 0 Ωs On the first term on the LHS of (30) we use the change of variable ϕ s , and obtain τ 0 Ωs where R s = R • ϕ s . For the second term on the LHS of (30) we use integration by parts in space and obtain Proof. We start on the LHS and we use the change of variable ϕ s = I + u s . Thus we obtain τ 0 Ωs Recall the relation between P s and T s : Also, since R s = R • ϕ s , then Using (35) and (36) Now by the definition of matrix inner product, and the fact that trace is invariant under cyclic permutations, we have that Combining (37) with (38), we obtain the desired equality (33).

Shape derivative on O.
We will deal first with the term on the LHS of (39).
where R represents the shape derivative of R s w.r.t s at s = 0, i.e.
For sake of notation, letR = R • ϕ(t), and thus we rewrite (41) and (40) as Our goal is to find the equation that ϕ (t) satisfies on O. In the second term on the RHS of (42) we integrate by parts in space and obtain that ( P · n S ,R + < P · n S , [(DR)(Dϕ(t)) −1 ]ϕ (t) ) dXdt.
Now we make use of the following identities: Let a, and f be two vector fields in R 3 , and letR ∈ C 1 (I × D). Then the following identity is satisfied: Moreover, if we assume that DR.n S = 0 on S, then Using identities (44) and (45) in (42) and (43), we obtain that Combining (42) with (46) and (47), we obtain the final expression for the shape derivative of the domain integrals: Next, we calculate the shape derivative on the RHS of (39).

5.3.2.
Shape derivative of the boundary integral. We can not directly calculate the shape derivative w.r.t. s of the term τ 0 Γs(t) p s n s − 2νε(w s ) · n s , R dΓ s (t) dt, since Γ s (t) is not a full boundary. If we choose the test function R such that R = 0 on Γ in ∪ Γ out , then the boundary integral turns into the following volume integral: where the speed v is given by v = V t (0) · n.
On the second term on the RHS of (50) we use the tangential Green's formula and obtain: In order to simplify the last term of the RHS of (51), we choose the test function R such that DR · n S = 0 on S. This is equivalent to DR · n(t) = 0 on Γ(t). Using the choice mentioned above, we obtain the following identities: D( pR) · n, n = p DR · n, n + ( ∇pR * ) · n, n = ∂ ∂n p R, n , and (52) Combining (50) with (51), (52), and (53), we obtain our final expression for the shape derivative of the boundary integral: (48) and (54), we obtain that u satisfies the following weak form identity:

Equation for u (t) on O. From
Therefore we recover the following linearized equation for u on O: We can simplify this further to: But since u solves the elastodynamics equation on O, we have that ∂ 2 ϕ(t) ∂t 2 − DivP = 0 on O, and therefore we obtain the desired equation for u on O: Using the relation between P and T that we obtained in Section 3.2, we can also find the equivalent linearized elasticity equation on Ω, in terms of T : . Equation for u on Γ(t). In order to recover the equation for u (t) on Γ(t), one way would consist in transporting all the boundary terms on S present in (55) back on Γ, namely the following expression: This tedious work can be avoided as follows. We use an equivalent weak formulation for the elastodynamics equation: ∀R ∈ C 1 (I × D), and with R s = R • ϕ s , we have We take shape derivative w.r.t. s at s = 0 in (57) and obtain: We already calculated the term on the RHS of (58). Its derivative is given by (54). We also have an expression for the first term on the LHS of (58): The only term on the boundary is τ 0 S ∂ 2 ∂t 2 ϕ(t),R (Dϕ(t)) −1 ϕ (t) · n S dsdt. However this term is on S, and thus we need to use the change of variable ϕ −1 to bring it back on Γ(t): S ∂ 2 ∂t 2 ϕ(t),R (Dϕ(t)) −1 ϕ (t) · n S dsdt = S u , (Dϕ) − * n S ϕ tt ,R ds where the density ω is given by ω = det(Dϕ −1 )|(Dϕ −1 ) − * n S |.