Analysis of a coupled fluid-structure model with applications to hemodynamics

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. SIAM J. NUMER.ANAL.c 2016 Society for Industrial and Applied Mathematics Vol. 54, No. 2, pp. 994–1019 ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL WITH APPLICATIONS TO HEMODYNAMICS∗ T. CHAC ´ ON REBOLLO†,V.GIRAULT ‡,F.MURAT §,AND O. PIRONNEAU§ Abstract. We propose and analyze a simplified fluid-structure coupled model for flows with compliant walls. As in [F. Nobile and C. Vergara, SIAM J. Sci. Comput., 30 (2008), pp. 731– 763], the wall reaction to the fluid is modeled by a small displacement viscoelastic shell where the tangential stress components and displacements are neglected. We show that within this small displacement approximation a transpiration condition can be used which does not require an update of the geometry at each time step, for pipe flow at least. Such simplifications lead to a model which is well posed and for which a semi-implicit time discretization can be shown to converge. We present some numerical results and a comparison with a standard test case taken from hemodynamics. The model is more stable and less computer demanding than full models with moving mesh. We apply the model to a three-dimensional arterial flow with a stent. Key words. fluid-structure interaction, compliant walls, finite element method, convergence analysis, Navier–Stokes equations, blood flow AMS subject classifications. 35Q30, 74K25, 65M60, 65N12 DOI. 10.1137/140991509 1. Introduction. Fluid-structure interaction (FSI) is computationally challenging because it involves moving geometries and the coupling of Lagrangian and Eulerian models [31, 41]; most popular applications are for biofluid dynamics, hemodynamics, and aerospace. This paper is a contribution to FSI algorithms, not to hemodynamics as such; but since we need to compare solutions we chose this field because it is well documented. Other applications like aircraft design and tires for instance have additional intrinsic difficulties which complicate the comparison with a simplified model. Computational hemodynamics has important applications (see [62, 27] or [61, 46], and the references therein). Modeling flow in a large blood vessel can be done with incompressible Navier–Stokes equations. The blood vessel is more difficult to model as it is a complex material for which the rheology is unclear because it is different in vitro from in vivo [62]. No doubt future computers will be able to handle this complexity and one will use large displacement nonlinear models for the structure [63]. However in the meantime there is a need for fast, well understood, and appropriate though less accurate models. To handle the complexity of moving walls, the method of immersed boundaries has been used—if not invented—by Peskin, the pioneer of computational hemodynamics [49, 48, 50, 64]. The mathematical analysis of this method is difficult [5, 42] and it is also hard to incorporate an elaborate viscoelastic model for the vessels, how- ∗Received by the editors October 14, 2014; accepted for publication (in revised form) January 13, 2016; published electronically April 5, 2016. http://www.siam.org/journals/sinum/54-2/99150.html †IMUS & Departamento de Ecuaciones Diferenciales y An´alisis Num´erico, Universidad de Sevilla, 41012 Sevilla, Spain, and Laboratoire Jacques-Louis Lions (LJLL), UPMC, 75005 Paris, France (c[email protected]). The research of this author was partially funded by the grants MTM2012-36124- C02-01 and MTM2015-64577-C2-1-R of the Spanish Ministerio de Econom´ıa y Competitividad and the FEDER EU Fund. ‡Laboratoire Jacques-Louis Lions (LJLL), Sorbonne Universit´es, UPMC Univ Paris 06, UMR 7598, 75252 Paris Cedex 05, France ([email protected]). §Laboratoire Jacques-Louis Lions UPMC, 75005 Paris, France (francois.[email protected]r, olivier.pir[email protected]r). 994 Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 995 ever, it is in the process of being solved [4]. Most authors prefer to follow the moving boundaries and use separate models for the fluid and the structure and couple them at each time step. On this topic the literature is large; let us just mention some references [13, 10, 16, 65, 27, 2, 34, 24, 43]. Linear elasticity with small displacement for the vessels can be applied either in three dimensions (see [26, 7], for instance) on a thick wall or in two dimensions via a shell model as in [45] and [10, 11] (see also [27]). Still the fact that elasticity is written on a fixed domain while the fluid domain is moving creates a computational difficulty and causes instabilities which also generated a large literature (see, for instance, [45, 26, 23, 24, 20, 21]). Indeed, if time discretization leads to a coupled problem at each time step, the problem is difficult to solve, while a loosely coupled system is slow, cursed by the added mass effect [12]. In this paper our aim is to build a fast and stable algorithm by pursuing two leading principles: (1) avoid remeshing because it is somewhat incompatible with unconditional stability; (2) consistent and energy preserving schemes are likely to be robust and stable. We propose to do so by investigating a system which is derived from Nobile and Vergara’s variational fluid/shell-structure model, simplified by using a boundary condition projected from the moving mesh onto a fixed mesh—an idea that is also developed in [17, 32], but which we validate with the so-called transpiration approximation. The model bears similarities to the one studied in [2] where a monolithic variational formulation with a vector displacement and linear elasticity is proposed on a fixed domain. In this reference existence and regularity of the solution are assumed and space-time error analysis is studied at zero Reynolds number together with the stability of a fractional step algorithm. In contrast, in our study, the displacement of the fluid-structure interface is normal to itself and the transpiration approximation justifies neglecting the time dependence of the fluid domain; existence and regularity are shown with appropriate conditions on the data, unconditional stability, and optimal error in time. The convergence of the spatial approximation by the finite element method is not shown because the condition u×nat the interface leads to difficulties which require a lengthy analysis. Following [47], Nobile and Vergara [45] proposed a second approximation, beside small displacements, namely, that lateral displacements can be neglected in the shell model. Then they showed that Koiter’s model reduces to a scalar equation for the normal displacement ηon the mean position Σ of the vessel’s wall. With a viscoelastic prestress model [20, 10, 25], the normal displacement is governed by (1.1) ρsh∂ttη−∇c·(T∇cη)−∇c·(C∇c∂tη)+a∂tη+bη =σs nn,η,∂ tηgiven at t=0. Here ∇cis the covariant gradient, hdenotes the thickness of the vessel, and ρsits volumetric mass; Tis the stress (possibly prestressed) tensor, Cand aare visco-elastic damping terms and bis also an elastic parameter; the right-hand side is the external normal force on the shell, in our case the normal component of the normal stress at the surface of the solid σnn. Note that in this context and due to the assumption of normal displacements the other components of the normal stress tensor do not appear and hence cannot be matched with those of the fluid. A particular case is the so-called surface pressure model, when all but the last term are neglected in the left-hand side of (1.1): bη =σs nn. In a cylinder, b≈Ehπ A(1−ξ2), where Ais the vessel’s cross section, Ethe Young modulus, ξthe Poisson coefficient. For more complex shapes, bdepends on the coefficients of the map of the cross section to a reference circle, but (up to higher-order terms in η)bis not a function of η.A Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 996 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU general formula is given in [45, (2.2)], assuming that the vessel is shaped like a pipe with smooth and slowly varying cross sections. As in Koiter’s model, this rules out bifurcating pipes. So for general surfaces, bneeds to be measured. The surface pressure model is an interesting prototype with which to understand the complexity of FSIs. The structure model is coupled with the Navier–Stokes equations for the fluid by matching the velocity of the fluid to that of the structure and writing continuity of the normal component of the normal stress from the fluid side and the structure side. When the shell is replaced by a plate (a fourth order operator added to (1.1) above), existence has been established [15, 33, 35, 38]. The fourth order term seems essential to the existence of solutions. The difficulty appears to be connected with the fact that the fluid domain is moving because this paper provides an existence result when the fluid domain motion is neglected. In contrast, the transpiration approximation provides a framework to work on a fixed domain for the fluid. Transpiration is an old idea in CFD (see (4.17) in [52], for instance) and it has been used in the 1990s to analyze wing flutter [22] and for conditioning the fluidstructure coupling algorithm [20]. The idea is simple. If a boundary condition like v=ghas to be imposed on a part of the boundary, Σt={x+ηn:x∈Σ0},where η(x,t) denotes the motion of Σtmeasured in the normal direction n=n(x)ofits position at rest Σ0, then one may write, for all x∈Σ0, (1.2) v(y)=g(y)∀y:= x+ηn∈Σt⇔v(x)+η∇vT(x)n=g(x)+η∇gT(x)n+o(η). Such an approximation is in line with the small displacement hypothesis made in using linear elasticity. Typically for a large vessel like the aorta, a section of length of 5 to 10 cm is considered; the heart pulse is about 1 Hz and the pressure amplitude change is roughly 6 KPa; these numbers induce indeed small displacements [27], except perhaps near the heart where the displacement is due to the heart itself. The resulting model in variational form is analogous to the Navier–Stokes equations with zero tangential velocities on the boundaries plus a surface integral term coupling the pressure and the normal velocities. The model is driven by the initial velocities and possibly the pressure on inflow/outflow surfaces. Compared to [45], the mathematical analysis is greatly simplified by the fact that there is only one variational equation; compared to the Navier–Stokes equations, the difficulties are twofold. First, one has to deal with Sobolev spaces of vector-valued functions with prescribed tangential velocities, i.e., spaces larger than H1 0. Second, the problem is nonlinear and as for most existence proofs, infinite sequences must have converging subsequences in the right spaces for the nonlinear term to have a suitable limit. Despite these difficulties, we prove in this work that the model has a solution and that a time discretization converges to it. As η=t 0u·nthe result gives also the regularity of η. Some hypotheses may be questionable for applications: (1) The boundary is assumed piecewise-C1,1with no obtuse angle at the edges. A triangulated pipe satisfies these, but a bifurcating pipe does not; a rugged wall would not either. (2) The parameter bshould be large compared to the other viscoelastic parameters of (1.1) otherwise the transpiration approximation is poor. The paper is organized as follows. In section 2 we elaborate on (1.1)–(1.2) and introduce the simplified coupled model which is the object of this study. Section 3 is devoted to introducing a semidiscretization in time and it is proved to be well posed in section 4. In section 5 we prove the stability and convergence of the semidiscretization. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 997 Finally section 6 presents some relevant numerical tests and comparison to other models. 2. Modeling simplifications. Assume the flow to be Newtonian in the time varying domain Ωt⊂R3and let Σtbe the part of the boundary at the interface between the solid and the fluid at time t. WedenotebyΣ(resp.,Ω)areference position of Σt(resp., Ωt); Σ could be its initial position Σ0.Letnbe the normal to Σ pointing outside Ω. If η(x,t) denotes the normal displacement of Σtwith respect to Σ at position xand time tthen Σt={x+η(x,t)n(x); x∈Σ}.The Navier–Stokes equations link the fluid velocity uand the pressure pby ρf(∂tu+u·∇u)−∇·σf=0,∇·u= 0 in Ωt,(2.1) where ρfis the volumetric mass density of the fluid, σf=−pI+μ(∇u+∇uT)isthe stress tensor, and μthe viscosity. Equality of normal components of normal stress tensors at the vessel Σtand equality of velocities (the so-called no slip condition) give (2.2) −σs nn|x=σf nn|x+ηn=[p−μnT(∇u+∇uT)n]|x+ηn,∂ tηn|x=u|x+ηn. Recall that σs nn is related to ηby (1.1), or—if simplified—by bη =σnn. Notice that the second condition implies u×n=0on Σt. Recall also that pbeing defined up to a constant, this constant has to be adapted too, at time t= 0, for instance. A variational formulation for (1.1), (2.1), (2.2) has been given for a similar model in [65] and studied in [45] for this model. See also in [17] an approximation which neglects the motion of the moving domain and to which we will compare our numerical results. However it seems that our use of the following is new. Proposition 2.1. Assume that Σis smooth near ˜ x, i.e., there is a neighborhood Vof ˜ x∈Σin which Σis C2. Assume that uis C1,∇·u=0in Ωt∩V,andu×n=0 on Σ∩V,thenat˜ x∈Σ∩V, 1 2nT(∇u+∇uT)n=∂u ∂n ·n=−2u·n R+O(|∇R−1|)+O(|∇r−1|),(2.3) where R, r are the principal radii of curvature and R=(R−1+r−1)−1is the mean radius of curvature of Σat ˜ x. Proof. Let us work with simple toroidal coordinates (r, θ, φ)→(x=˜ Rcos φ, y =˜ Rsin φ, z =rsin θ), where ˜ R=R0+rcos θ,sothat˜ xis at φ=θ=0andΣ tis tangent to the torus (R, r), i.e., rand R:= R0+rare its two principal radii of curvature at ˜ x. Recall that (see, for example, [36, Appendix A]) ∇·u=hrhθhφ∂r ur hθhφ +∂θ uθ hφhr +∂φ uφ hrhθ (2.4) with hr=1,h θ=1 r,h φ=1 ˜ Rbecause, by definition, 1 h2 k=(∂kx)2+(∂ky)2+(∂kz)2, k=r, θ, φ.So∇·u=0,uθ=uφ=0,andu×n=0imply (2.5) ∇·u=∂rur+ur R0+2rcos θ r(R0+rcos θ)=0 ⇒∂rur=−ur1 r+cos θ ˜ R=−2ur Rat ˜ x. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 998 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU Similarly ∇u=ieiThi⊗∂ikekTuk,i, k ∈(r, θ, φ), with er=cosθ(cos φ, sin φ, tan θ), eθ=sinθ(−cos φ, −sin φ, 1 tan θ), eφ=(−sin φ, cos φ, 0). As uθ=uφ=0,wehave nT(∇u)n=er(er⊗erhr∂rur)er=∂rurwhich is also ∂u ∂n·nand also 1 2nT(∇u+∇uT)n, bythefactthatitisascalar. Remark 2.1. Note that the argument is local and requires that locally the surface Σ be tangent to a torus, which is possible when the surface is locally C2;itdoesnot require Σ to be quasi-toroidal globally. On the other hand the surface pressure model for brequires a pipelike shape. If it is not so, then bmay have to be measured or computed by an inverse problem. 2.1. Transpiration approximation. Following (1.2) and using (2.3), (2.2) becomes −σs nn =p+4μ Ru·n+η∂p ∂n +4μ R ∂u ∂n ·n,∂ tη=u+η∂u ∂n·non Σ.(2.6) Thanks to (2.3), setting α=4μ R(1 −2 R), in the case of the surface pressure model, bη =p+αu·n+η∂p ∂n,∂ tη=1−η2 Ru·non Σ.(2.7) The last term is second order and may be dropped. In laminar flows at Reynolds number of a few thousands or less, 1 ρf ∂p ∂n is not large because pdoes not vary much in the boundary layer near Σ. This is not proved mathematically, but is observed experimentally; so when the problem is adimensionalized (a division by ρfis required) and the pressure constant is adapted, then both p ρfand 1 ρf ∂p ∂n are O(1) rather than O(μ−k)forsomek>0. Hence, it is reasonable to assume the following: (2.8) b∂p ∂n|Σ(Hypothesis 1). If this is true, then the quantity (b−∂p ∂n)ηthat appears in the shell model can be approximated by bη. The matching conditions on Σ then reduce to (2.9) u×n=0,U n:= t 0 u·nds, p +αu·n=bUn+cwith c:= p0+αu0·n. For the full shell model (1.1), on Σ, u×n=0,U n:= t 0u·nds,and (2.10) ρsh∂ttUn−∇c·T∇cUn−∇c·C∇c(u·n)+(a−α)u·n+bUn=p−c. 2.2. Other boundary and initial conditions. For a numerical study, the vessel is usually taken out of its environment and two artificial cross sections Γ are introduced. We shall assume that the dynamic pressure or the normal stress are given on Γ; more generally, •on Γp⊂Γ the flow is normal, i.e., u×n=0, and the dynamic pressure is given by p+ρf 2|u|2=pΓ; •on Γf=Γ\Γpthe normal stress is given by −pn+μ(∇u+∇uT)n=g. The difference between the two conditions has been discussed in [3]. Initial conditions are needed for (2.1), for instance, u(0) = u0in Ω,p(0) = p0on Σ.(2.11) Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 999 2.3. Energy considerations. Now there is a problem! The system of equations (2.1)+(2.10)+(2.11) no longer preserves energy because the integral over Ω of u·∇u gives a nonzero surface integral; therefore it is necessary to replace the nonlinear term u·∇uin (2.1) by −u×(∇×u) according to the identity u·∇u=−u×(∇×u)+∇|u|2 2.(2.12) Indeed, this will replace the pressure by the dynamic pressure p+1 2|u|2.OnΣthis change is compatible with the small displacement hypothesis because the change is |u|2 |Σwhich is second order; it is not so inside Ω. Finally recall the identities −∇·(∇u+∇uT)=−Δu=∇×(∇×u)−∇(∇·u).(2.13) The modified Navier–Stokes system for FSIs written in a fixed domain Ω is now, with u×n|Σ=0,U n:= t 0u·nds, (2.14) ∂tu−u×(∇×u)+ν∇×(∇×u)+∇˜p=0,∇·u=0, γ∂ttUn−∇c·˜ T∇cUn−∇c·˜ C∇c(u·n)+(˜a−˜α)u·n+˜ bUn=˜p−c, where cis a compatibility constant, γ:= hρs/ρf,ν:= μ/ρf,˜p:= p/ρf+|u|2/2, and for other variables the ˜.indicates a division by ρf;inparticular,˜ b:= b ρf. In the sequel we drop the tilde over pand band rename a←a−α. Remark 2.2. Notice that on a fixed domain with zero velocity at the walls, the Navier–Stokes equations are strictly equivalent to (2.14). For compliant walls this equivalence holds too but within the approximation of small wall displacements. 2.4. The mathematical problem. Also as is usually the case for the mathematical analysis of Navier–Stokes equations and to some extent without loss of generality [29], we assume that Γ = ∅, i.e., Σ = ∂Ω and we replace gand pΓby a volumetric force f. A remark will be given later concerning the full variational formulation without this hypothesis. Therefore we shall consider the system (2.15) ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂tu−u×(∇×u)+ν∇×(∇×u)+∇p=f,∇·u=0 in QT, u(0) = u0in Ω,u×n|Σ=0,U n:= t 0 u·ndson ST, γ∂ttUn−∇c·T∇cUn−∇c·C∇c(u·n)+au·n+bUn+c=pon ST, where QT=(0,T)×ΩandST=(0,T)×Σ. This complete shell model on a fixed domain is energy conservative for smooth solutions. This may be proved in a standard way from the weak formulation (Definition 5.1). There have been objections to the numerical use of −u×(∇×u) in [40] as not fit for boundary layers. But here again by (2.12) the term is identical to u·∇u−1 2∇|u|2; the two expressions are different only when discretized in space. So this problem will be addressed later. 2.5. Weak form. We use the Sobolev space Wk,p(Ω) and denote its norm by ·k,p,Ω.WedenoteWk,p(Ω) = [Wk,p(Ω)]3,Hk(Ω) = Wk,2(Ω), Hk(Ω) = [Hk(Ω)]3, Lk(Ω) = [Lk(Ω)]3. In this work, we assume that Ω is a bounded, connected, Lipschitz Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 1000 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU domain, but we shall also consider domains which in addition to this property, are curvilinear polyhedra of class C1,1(see [19]). Loosely speaking, let us say that the boundary ∂Ω of Ω consists of a finite number of faces which are open subsets of surfaces of class Ck,1(here k= 1), Ω is locally on one side of its boundary, and any two adjacent faces are connected by smooth edges not forming a cusp. Thus all dihedral angles of ∂Ω are bounded away from zero and 2π. In this paragraph, C denotes a generic constant that is independent of the functions involved. We study problem (2.15) when Σ = ∂Ω to avoid unnecessary technical complexities. To derive a variational formulation of (2.15), let us consider the space  W={w∈L2(Ω) |∇×w∈L2(Ω),∇·w∈L2(Ω),w×n|Σ=0}.(2.16) Then the following lemma holds (cf. Amrouche et al. [1, Proposition 3.7]). Lemma 2.2. The space  Wis well defined and is a Hilbert space endowed with the norm w W=(w2 0,2,Ω+∇×w2 0,2,Ω+∇·w2 0,2,Ω)1/2.Moreover when Ωis either convex or C1,1,then Wis continuously imbedded into H1(Ω) and there exists a constant C>0such that w1,2,Ω≤Cw W∀w∈ W. A straightforward extension of Proposition 3.7 in [1] (see [19]) shows that the continuous imbedding of  Winto H1(Ω) holds when Ω is a bounded, connected, Lipschitz, curvilinear polyhedron of class C1,1without reentrant corner, i.e., with angles bounded away from zero and π. In [29] it is proved that, in a Lipschitz domain, the normal trace w·nis a bounded linear mapping from H(div,Ω) into H−1/2(∂Ω). Let us introduce the velocity space W={w∈ W|w·n∈H1(∂Ω) },(2.17) which is a Hilbert space endowed with the norm w2 W=w2  W+w·n2 1,2,∂Ω.By definition, the functions wof Whave w=gnon ∂Ω,with gin H1(∂Ω). Thus, if w∈Wbelongs to H1(Ω) when ∂Ω has corners, then necessarily gvanishes at the edges of ∂Ω. In other words, if Sdenotes a face of ∂Ω, then g∈H1 0(S). More precisely, we have the following characterization. Lemma 2.3. Let Ωbe a bounded, connected, Lipschitz curvilinear polyhedron of class C1,1without reentrant corner. Let Si,1≤i≤K,bethefacesof∂Ωand let ni be the unit exterior normal to Si. Then each function wof Whas the decomposition (2.18) w=w0+ K  i=1 GiNi, where w0∈H1 0(Ω),andforeachi,1≤i≤K,Gi∈H3/2(Ω),Gi|Si∈H1 0(Si), Gi|Sj=0for j=i,Ni∈W1,∞(Ω),andNi|Si=ni. Conversely, all functions of the form (2.18) belong to W. Proof. Let wbelong to W. The assumption on ∂Ω implies that wis in H1(Ω) and satisfies w=gn∈H1 0(Si) on all faces Siof ∂Ω. Let gi∈H1(∂Ω) denote the restriction of gto Si, extended by zero to all other faces. It has a continuous extension Gi∈H3 2(Ω). Next, since nibelongs to W1,∞(Si), it has a continuous extension in W1,∞(∂Ω), and a continuous extension Ni∈W1,∞(Ω). Then the product GiNi belongs to H1(Ω) and we set w0=w−K i=1 GiNi.It is easy to check that this is the desired decomposition. The converse is obvious. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 1001 The argument of Lemma 2.3 yields the following density result. Theorem 2.4. Under the assumptions of Lemma 2.3, the space W∩W1,∞(Ω) is dense in W. Proof. Let wbelong to W. On one hand, w0can be approximated in H1 0(Ω) by a sequence wm 0in D(Ω)3. On the other hand, gican be approximated in H1 0(Si)by a sequence gm iin W1,∞ 0(Si). As in Lemma 2.3, gm ihas an extension Gm iin W1,∞(Ω) such that Gm i|Sj=0forj=iand Gm i|Si=gm i.Thenwm=wm 0+K i=1 Gm iNi belongs to W∩W1,∞(Ω) and converges to win W. Note that the normal trace is bounded and linear from H1(Ω) into L4(Σ) such that w·n=(γ0w)·nfor all w∈H1(Ω) (γ0denotes the trace; see [14, Lemma A.1]). 3. Semidiscretization in time. In this section we propose a variational formulation for the Navier–Stokes boundary value problem (2.15). Our analysis is inspired by the early works on the Stokes and Navier–Stokes equations with boundary conditions on the pressure (cf. [3, 18]). We shall assume the following. Hypothesis 1. The functions γ, a, b, T,C,csatisfy •γ, a, b ∈L∞(Σ),T,C∈L∞((0,T)×Σ),c∈L2((0,T)×Σ); •γ≥0,a≥0,b≥0a.e. on Σ,andthereexistνT>0,νC>0such that (3.1) αT(w,w)≥νT∇c(w·n)2 0,2,Σ, αC(w,w)≥νC∇c(w·n)2 0,2,Σ∀w∈W, where, for any vand win W,theformsαTand αCare defined by αT(v,w)=Σ∇T c(v·n)T∇c(w·n)dΣ, αC(v,w)=Σ∇T c(v·n)C∇c(w·n)dΣ, and ∇cis the covariant gradient on Σ(intuitively the tangential gradient). The condition (3.1) formulates the ellipticity of forms αTand αC. In (2.10) cis a constant; there are cases where it could be a function of time. As it does not complicate the analysis, here it is assumed to be a function of x, t.We consider the following implicit semidiscretization in time of (2.15) by linearization of the convection term. For the pressure, let M=L2(Ω).Assume f∈L2(W), u0∈W. Let N≥1 be an integer and set δt =T/N.Setu0=u0. Discretization. For n=0,1,...,N −1, find (un+1,p n+1)∈W×Msuch that for any (w,q)∈W×M, un+1 −un δt ,wΩ−(un+1 ×(∇×un),w)Ω+ν(∇×un+1,∇×w)Ω −(pn+1,∇·w)Ω+γun+1 −un δt ·n,w·nΣ +αT(Un+1,w) +αC(un+1,w)+(aun+1 ·n,w·n)Σ+(bUn+1 ·n,w·n)Σ =fn+1,w−(cn+1,w·n)Σ,(3.2) (∇·un+1,q)Ω=0,(3.3) where Un+1 =δt n+1 k=1 uk,fn+1 =1 δt tn+1 tnf(s)ds, and cn+1 =1 δt tn+1 tnc(s)ds. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 1002 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU 4. Analysis of the semidiscrete problem. In this section, for fixed n,we prove existence and uniqueness of the solution (un+1,p n+1) of (3.2)–(3.3). Estimates for un+1 and for a primitive in time of the pressure (instead of the pressure itself) will be obtained in Theorem 5.4 below (see (5.9)). Problem (3.2)–(3.3) is an Oseen-like problem, however, it is nonstandard due to the structure of the convection term, and the presence of the boundary terms issued from the discretization of the vessel shell model. Its well posedness is based upon an inf-sup condition. Lemma 4.1. Assume that the domain Ωis a bounded, connected, Lipschitz, curvilinear polyhedron of class C1,1. Then for some β>0, (4.1) βq0,2,Ω≤sup w∈W (q,∇·w) wW for all q∈M. Proof. Let q∈M. The idea is to use the argument of Boland–Nicolaides (cf. [6]). Thus, we decompose q=q+q0with q=1 |Ω|Ωq,q0=q−q.ThenΩq0=0andthe standard theory [29] ensures that there exists w0∈H1 0(Ω) such that (4.2) ∇·w0=q0,w01,2,Ω≤1 β0q00,2,Ω for some β0>0 independent of q0.Totreatq,wechooseafaceSof ∂Ω and consider a function ρ∈C2,1(Ω) such that ∂Ωρ=1,supp(ρ)∩∂Ω⊂S.Let κ= supp(ρ)∩S. We assume that κis a compact subset of Sthat does not intersect the edges of ∂Ω. Let nbe the normal to S.LetNin W1,∞(Ω) be an extension of nas in Lemma 2.3, and consider the function w1=ρN.Thenw1∈Wand moreover Ω∇·w1=∂Ωw1·n=∂Ωρ=1.Then, if q>0, q0,2,Ω=|Ω|1/2|q|=|Ω|1/2w1W (q,∇·w1)Ω w1W≤β1sup w∈W (q,∇·w)Ω wW , where β1=|Ω|1/2w1W.Ifq<0, it suffices to change the sign of w1in the second equality. To prove the global inf-sup condition, following Boland–Nicolaides, we set z=w0+α q w1with α∈R+to be determined. Then (q,∇·z)Ω=q02 0,2,Ω+α q (q0,∇·w1)+α|q|2 ≥q02 0,2,Ω+α |Ω|q2 0,2,Ω−αC|q|q00,2,Ω ≥q02 0,2,Ω+α |Ω|q2 0,2,Ω−αC 2γq2 0,2,Ω |Ω|+γ−1q02 0,2,Ω ≥1 2q02 0,2,Ω+1 C2q2 0,2,Ω |Ω|≥1 2min 1,1 C2|Ω|q2 0,2,Ω,(4.3) where C=∇·w10,2,Ω,γ=1 C,andα=1 C2. Finally, (4.4) zW≤w01,2,Ω+α|q|w1W≤Cq0,2,Ω, and (4.1) follows easily from (4.3) and (4.4). In what follows, unless otherwise stated, we shall assume the following. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 1009 This lemma implies in particular for all zin L2(L6)andanyr∈(2,6), (5.20) zL2(Lr(Ω)) ≤zθ L2(L2(Ω))z1−θ L2(L6(Ω)), hence the next lemma. Lemma 5.6. Let Ωbe bounded, and let {vδ}δ>0and vin L2(L6)satisfy (5.21) vδL2(L6)≤Cand lim δ→0vδ−vL2(L2)=0 with Cindependent of δ. Then for all real r,1≤r<6,limδ→0vδ−vL2(Lr)=0. Recall also the following result. Lemma 5.7. Let the domain Ωbe bounded. Assume that the sequence {vδ}δ>0⊂ L3(QT)strongly converges to vin L3(QT),letϕ∈C ∞([0,T]),andw∈W.Then vδ(x,t)⊗w(x)ϕ(t)strongly converges to v(x,t)⊗w(x)ϕ(t)in L2(QT)3×3. We also need the following compactness result for space-time functions (cf. [58]). Lemma 5.8. Let X,E,Ybe Banach spaces such that X→E→Ywhere the embedding X→Eis compact. Then the embedding below is compact: Lp(0,T;X)∩Ns,p(0,T;Y)→Lp(0,T;E)with 0<s<1,1≤p≤+∞. We are now in a position to state the convergence result. Theorem 5.9. Assume that the domain Ωsatisfies Hypothesis 2and Σ=∂Ωis connected. Assume that u0∈W,f∈L2(W). Then the sequence (uδ,p δ)δ>0contains a subsequence (uδ,p δ)δ>0that is weakly convergent in L2(W)×H−1(L2)toaweak solution (u,p)of the boundary value problem (2.15) in the sense of Definition 5.1. Moreover (uδ)δ>0is weakly-* convergent in L∞(L2)to u,stronglyinL2(Lr)for 1≤r<6, and the primitives in time (Pδ)δ>0of the pressures (pδ)δ>0are weakly-* convergent in L∞(L2)to a primitive in time Pof the pressure p. If the solution of (5.1) is unique, then the whole sequence converges to the solution. Proof. We proceed in three steps. Step 1. Convergent subsequences: By estimates (5.9) and (5.10), uδis uniformly bounded in L2(H1), in L∞(L2), and in N1/4,2(L2). The embedding H1(Ω) →Lr(Ω) is compact for 1 ≤r<6(cf.Br´ezis [9, Chap. 9]), hence the compactness of the embedding W→Lr(Ω). Applying Lemma 5.8 with X=H1(Ω), E=Lr(Ω), Y= L2(Ω), p=2,ands=1 4, it follows that the sequence (uδ)δ>0is compact in L2(Lr) for 1 ≤r<6. By estimate (5.10), the sequence (Pδ)δ>0is uniformly bounded in L∞(L2). Then the sequence (uδ,P δ)δ>0contains a subsequence (still denoted with δ) such that (uδ)δ>0is strongly convergent in L2(Lr)tosomevfor any 1 ≤r<6, weakly in L2(H1), and weakly-* in L∞(L2), and (Pδ)δ>0is weakly-* convergent in L∞(L2)to some P. We prove in what follows that the pair (v,∂ tP) is a weak solution of Navier– Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 1010 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU Stokes equations (5.1) in the sense of Definition 5.1. To simplify, in the remainder of this proof, all convergences below are up to subsequences. Regarding ˜ uδ, it is easy to prove that the sequence (˜ uδ)δ>0is uniformly bounded in L∞(L2)andL2(H1). Moreover, by summing (5.12) over n,uδ−˜ uδ2 L2(L2)= N−1 n=0 tn+1 tn(tn+1−t δt )2un+1 −un2 0,2,Ωdt=δt 3N−1 n=0 un+1 −un2 0,2,Ω≤Cδt.Since we have, in particular, that limδt→0uδ=vin L2(L2), then this estimate implies that limδt→0˜ uδ=vin L2(L2),and Lemma 5.6 implies for any r,1≤r<6, limδt→0˜ uδ= vin L2(Lr).Clearly, the same limit can be derived for ˜ u− δ,asuδ−˜ u− δ2 L2(L2)= N−1 n=0 tn+1 tn(t−tn δt )2un+1 −un2 0,2,Ωdt≤Cδt. Step 2. Limit of (3.3):Let q∈L2(Ω) and ϕ∈D(0,T) be arbitrary. As ∇·uδweakly converges to ∇·vin L2(L2), we have T 0(∇·v(t),q)Ωϕ(t)dt= limδt→0T 0(∇·uδ(t),q)Ωϕ(t)dt. Consequently, T 0(∇·v(t),q)Ωϕ(t)dt=0∀q∈ L2(Ω),ϕ∈D(0,T).As D(Ω) ⊗D(0,T) is sequentially dense in D(QT), we deduce that ∇·v= 0 a.e. in Ω ×(0,T). Step 3. Limit of (3.2):We reformulate (5.13) as (5.22) −T 0 (uδ(t),w)Ωϕ(t)dt −(u0,w)Ωϕ(0) + T 0 e(˜ u− δ(t); ˜ uδ(t),w)ϕ(t)dt +T 0 d(˜ uδ(t),w)ϕ(t)dt+T 0 (Pδ(t),∇·w)Ωϕ(t)dt −T 0 (γuδ(t)·n,w·n)Σϕ(t)dt−(γu0·n,w·n)Σϕ(0) +T 0 αT(˜ Uδ(t),w)Σϕ(t)dt+T 0 αC(˜ uδ(t),w)Σϕ(t)dt +T 0 (a˜ uδ(t)·n,w·n)Σϕ(t)dt+T 0 (b˜ Uδ(t)·n,w·n)Σϕ(t)dt =T 0˜ fδ(t),wϕ(t)dt−T 0 (˜cδ,w·n)Σϕ(t)dtfor all w∈W, and ϕ∈C∞([0,T]) such that ϕ(T) = 0. The weak convergence of uδin L2(L2) yields limδt→0T 0(uδ(t),w)Ωϕ(t)dt=T 0(v(t),w)Ωϕ(t)dt. Next, we pass to the limit in the convection term. Since wϕbelongs to L∞(L6)and˜ uδconverges to vstrongly in L2(L3), then the product ˜ uδ×(wϕ) converges strongly to v×(wϕ)inL2(L2). In addition, ∇ט u− δtends weakly to ∇×vin L2(L2). Hence lim δt→0T 0 (˜ uδ(t)×(∇ט u− δ(t)),w)Ωϕ(t)dt=T 0 (v(t)×(∇×v(t)),w)Ωϕ(t)dt. As ˜ uδ(t) weakly converges to vin L2(H1), limδt→0T 0d(˜ uδ(t),w)ϕ(t)dt= T 0d(v(t),w)ϕ(t)dt. To treat the boundary terms, the weak convergences of uδand ˜ uδin L2(W) and the continuity of the normal trace for the weak topology in H−1 2(Σ) Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 1011 imply that lim δt→0T 0 (uδ(t)·n,w·n)Σϕ(t)dt=T 0 (v(t)·n,w·n)Σϕ(t)dt, lim δt→0T 0 αC(˜ uδ(t),w)Σϕ(t)dt=T 0 αC(v(t),w)Σϕ(t)dt, lim δt→0T 0 (a˜ uδ(t)·n,w·n)Σϕ(t)dt =T 0 (av(t)·n,w·n)Σϕ(t)dt. Now, we turn to ˜ Uδ. According to (5.9), ˜ Uδ·nconverges weakly-* to some function lin L∞(H1(Σ)). Let w∈W,ϕ∈L1(0,T). Then w(x)·n(x)ϕ(t)∈L1(L2(Σ)), thus limδt→0T 0(˜ Uδ(t)·n,w·n)Σϕ(t)dt=T 0(l(t),w·n)Σϕ(t)dt. To identify the limit lwe first use Green’s formula. As ∇·˜ uδ= 0, we can write that for all σ∈ H1(Ω), (˜ uδ(t)·n,σ)Σ=( ˜ uδ(t),∇σ)Ωa.e. in (0,T). Next, we observe that ∀t∈ (tn,t n+1),˜ Uδ(t)=tn+1 0˜ uδ(s)ds=t 0˜ uδ(s)ds+tn+1 t˜ uδ(s)ds. Thus, for all ϕin L2(0,T)andσin H1(Ω), with R=N−1 n=0 tn+1 tn(tn+1 t(˜ uδ(s),∇σ)Ωds)ϕ(t)dt, T 0 (˜ Uδ(t)·n,σ)Σϕ(t)dt=T 0t 0 ˜ uδ(s)ϕ(t)dsdt, ∇σΩ +R. By interchanging the above integrals, a straightforward computation yields |R|≤ √δt(T 2)1/2∇σL2(Ω)˜ uδ0,2,QT.Thus, taking the limit δt →0, applying Green’s formula and ∇·v=0,weinferT 0(l(t),σ)Σϕ(t)dt=( T 0t 0v(s)ϕ(t)dsdt, ∇σ)Ω= T 0(U(v)(t)·n,σ)Σϕ(t)dt. Then (l(t),σ)Σ=(U(v)(t)·n,σ)Σfor all σ∈H1(Ω) a.e. in (0,T), and hence l=U(v)·n.Sincel∈L∞(H1(Σ)), we also have U(v)·n∈ L∞(H1(Σ)). Hence limδt→0T 0αT(˜ Uδ(t),w)Σϕ(t)dt=T 0αT(U(v)(t),w)Σϕ(t)dt, limδt→0T 0(b˜ Uδ(t)·n,w·n)Σϕ(t)dt=T 0(bU(v)(t)·n,w·n)Σϕ(t)dt. To pass to the limit in the pressure term, observe that owing to (5.10), (Pδ)δ>0is weakly-* convergent in L∞(L2) to some function P, limδt→0T 0(Pδ(t),∇·w)Ωϕ(t)dt= T 0(P(t),∇·w)Ωϕ(t)dt. Finally, as ˜ fδstrongly converges to fin L2(W)and˜cδ strongly converges to cin L2(H1(Σ)), limδt→0T 0˜ fδ(t),wϕ(t)dt=T 0f(t),wϕ(t)dt, limδt→0T 0(˜cδ(t),w·n)Σϕ(t)dt=T 0(c(t),w·n)Σϕ(t)dt. This concludes the proof. 6. Full discretization and numerical tests. 6.1. A finite element method. We introduce in this section a feasible discretization of problem (2.15) in polyhedral domains. Let Thbe a triangulation made of Ktetraedra {Tk}K 1with the usual conformity hypotheses. Consider the P2–P1 element, (or the P1+bubble−P1, see, e.g., [54] or [29]), built from Vh={v∈C0(Ω)3:v|Tk∈(P2)3,∀k∈T h},Q h={q∈C0(Ω) : q|Tk∈P1,∀k∈T h}. Note that the functions of Vhdo not satisfy v×n=0on the boundary. Indeed, as shown in [30] it would be vain to require v×n=0in strong form unless the Nedelec elements of degree at least 2 be used. Therefore the constraint u×n|Σ=0 will be implemented below by penalty. Recall our notation: the compliant boundary is Σ; Γ denotes the inflow/outflow boundaries and we assume that either the dynamic Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 1012 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU pressure is set to pΓand the flow is normal on Γpor the flux is set to gon Γf;we must have Γp∪Γf= Γ. Thus with the surface pressure model a feasible discretization of (2.15) is to find un+1 ∈Vh,p n+1 ∈Qhsuch that for all w∈Vh,q∈Qh, (6.1) Ωw·un+1 −un δt −un+1 ×(∇×un)−pn+1∇·w−q∇·un+1 +νΩ (∇×un+1)·(∇×w)+1 Σ∪Γp (un+1 ×n)·(w×n) +Σ bw·(un+1δt +Un) +Σγun+1 −un δt ·w+T∇c(Un·n)+(C+Tδt)∇c(un+1 ·n) ·∇c(w·n)+aun+1 ·w =Ω fn+1 ·w−Γp pΓw·n+Γf g·w,Un+1 =Un+un+1δt. For a convenient implementation Uis defined everywhere, not just on Σ. Letting w=un+1,q=−pn+1 gives the following energy estimate: (6.2) 1 2δt(un+12 0,2,Ω−un2 0,2,Ω)+ 1 2δtun+1 −un2 0,2,Ω+ν∇×un+12 0,2,Ω +1 un+1 ×n2 0,2,Σ+bδt 2un+1 ·n2 0,2,Σ +1 2δt(Un+1 ·n2 0,2,Σ−Un·n2 0,2,Σ) +1 2δt Σ[∇c((Un+1 −Un)·n)]TT∇c((Un+1 −Un)·n) +Σ [∇c(un+1 ·n)]TC+Tδt 2∇c(un+1 ·n) +γ 2δtun+1 −un2 0,2,Σ+aun+12 0,2,Σ+γ 2δt(un+12 0,2,Σ−un2 0,2,Σ) =Ω fn+1 ·un+1 −Γ pΓun+1 ·n+Γf g·un+1. An analysis similar to that developed in the preceding section establishes the scheme’s stability. Moreover, we deduce un+1 ×n0,2,Σ≤C(fn+1W+pΓ0,2,Γp+g0,2,Γf)√. Curved boundaries are approximated by polyhedral surfaces, thus generating an error of order √hin the approximation of un+1 ×n=0. Then the optimal choice is =h. Remark 6.1. By replacing all exponents n+1 by n+1 2in (6.1) except in the time derivatives and the pressure term, one builds a stable second order scheme in time. But it is nonlinear at each time step. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 1013 6.2. Other discretizations of the nonlinear term. If u·∇u−1 2∇|u|2is preferred, the following approximation is consistent and preserves energy: (6.3) Ωu·∇u−1 2∇|u|2·w≈Ω (un·∇un+1)·w−1 2Σ (un·n)(un+1 ·w). Then we can use un+1 2×(∇×un+1 2)=un+1 2·∇un+1 2−1 2∇|un+1 2|2.Onemayalsouse characteristic-Galerkin schemes [54, 60] applied to (6.3), which we recall here briefly. It is based on the following formula, ∂tu+a·∇u|x,(m+1)δt =1 δt um+1(x)−um(x−am(x)δt)+O(δt).(6.4) While the above is easy to understand, it is better and not harder to use um+1 −um◦ Xm amwith Xm am(x)=Xam(mδt) the solution of dX dτ (τ)=am(X(τ)),X((m+1)δt)=x. Then we need to find un+1 ∈Vh,p n+1 ∈Qhsuch that for all w∈Vh,q∈Qh, (6.5) Ωw·un+1 −un◦Xm um δt −pn+1∇·w−q∇·un+1+νΩ (∇×un+1)·(∇×w) +1 Σ∪Γp (un+1 ×n)·(w×n)+Σ bw·(un+1δt +Un)−1 2Σ (un·n)(un·w) +Σγun+1 −un δt ·w+(T∇c(Un·n)+(C+Tδt)∇c(un+1 ·n)) ·∇ c(w·n)+aun+1 ·w  =Ω fn+1 ·w−Γp pΓw·n+Γf g·w,Un+1 =Un+un+1δt. Each time step requires the solution of a symmetric time independent linear system, so that the matrices need be constructed only once. The stability of the characteristic-Galerkin scheme is addressed in [54, 51, 53]. In practice quadrature rules are needed for computing the integrals of w·(un◦Xm um). The trapezoidal rule is stable [51], but too crude here and, based on our numerical experience, we use the Gauss rule exactly for polynomials of degree 5. Finally note that second order in time extensions have been proposed (see [8, 60]). To summarize, consistent energy preserving first and second order characteristic- Galerkin schemes can be used, but the analysis done for u·n=0needs to be extended. 6.3. Simulation and comparison with a two-dimensional documented case. In Bukaˇc et al. [10] the following test case is proposed: •Flow between two parallel compliant planes at initial distance 2R: the geometry is a rectangle (0,L)×(0,R), L=6,R=0.5; a symmetry condition is imposed on the horizontal axis and the two-dimensional Navier–Stokes equations are used coupled with the full shell model (1.1). •Zero pressure is imposed at {L}×(0,R) on the right; and on the left p=pmax 21−cos 2πt tmif t<t m, else 0; here pmax =210 4,t m=510 −3. •A uniform grid 60 ×10 and a time step δt =10 −4are used. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 1014 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU t=6ms t=8ms t=10ms t=12ms Fig. 1.Zoom near the compliant wall at different times; the motion of the wall is greatly amplified for clarity. The computational domain used in (6.5) is static and below the light blue line. For graphics only the mesh is moved by δtd,bysolving(6.6). The compliant wall is now in green while the one of [10] is displayed with a solid dark blue line. The doted lines are other author’s results (Quaini [56] (dashed red line) and kinematic coupled scheme [10] (dash-dotted green)), documented in [10]). The black lines are the edges of our triangulation moved by δtd. •The constants of the model are as follows: ν=0.035,˜ b=410 5,˜ T=2.510 4, h =0.1,ρs ρf=1.1. •Other constants are zero. The fluid is at rest at time zero. Our test differs slightly from [10] in that our inflow section height on the left is fixed while it changes with the mesh motion in [10]. Note also that the comparison needs some adjustment because our mesh is fixed, while the mesh in [10] moves at each time step. The full model requires that at each time step Σtbe moved along its normal by a quantity δtun+1 ·n. For comparison and for graphic enhancement of our results we use an auxiliary triangulation which is moved at each step by δtdas in [20], where −Δdn+1 =0in Ω,dn+1|Σ=un·n,dn+1|Γ=0.(6.6) More precisely every vertex qjof the triangulation moves by qj→ qj+δt κ dn+1. In theory κ= 1 but can be adjusted for graphic enhancement; however, (6.6) is expensive and slows down the algorithm; it is also a source of instability when κis too large and triangles overlap. Algorithm (6.5) has been used as it is more precise than (6.1) for this test case because of the singularity at the top left corner due to the incompatibility of the condition u×n=0on Σ ∩Γ (see subsection below). The method is programmed using freefem++ [37]; the script is given in Appendix A. The results are shown in Figure 1; they compare rather well with [10] (and those to which they are compared). There is a difference at the inflow section on the left, which is due to the fixed geometry. To draw these graphics we have scanned the picture in [10] and adapted the aspect ratio to theirs. 6.4. Performance. Comparison. To compare (6.1) and (6.5) we make a small change to the previous test case: the mesh is refined near the compliant wall, but the total number of vertices is unchanged. Figure 2 shows the dynamic pressures and the velocity vectors at t=tm computed with the two methods (plots for (6.1) are above those of (6.5) on Figure 2). Notice however the singularity at the top left corner which is due to the boundary Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 1015 Fig. 2.For the two sets of plots, the results using (6.1), i.e., using −u×∇×uare above the same using (6.5), i.e., characteristic-Galerkin based on u·∇u. The level lines of the dynamic pressures are shown (top 2figures) together with the velocity vectors at tm(bottom two figures). Notice the oscillations near the top left corner with (6.1), however, these do not develop into an instability as the convection/diffusion damps them. Fig. 3.The top plot displays velocity vectors with color code proportional to their norm computed by (6.5). The time step is ten times larger than previously used so that tmis reached with 5 iterations (top two plots). The lower plot corresponds to step 50, i.e., t=10tm. conditions. Indeed u×n=0on both sides of the corner implies u=0but the mesh is not fine enough to implement it, consequently uis rapidly varying and ∇×uis infinite. Yet developments of instability are not observed. Stability. The time step is multiplied by 10, so that 5 steps only are needed to reach tm, but we still do 50 steps (hence reaching 10tm) and display, on Figure 3, the results after 5 and 50 time steps by the characteristic-Galerkin method using (6.5). Method (6.1), based on −u×(∇×u), is equally stable and gives similar results but with a small singular region near the top left and right corners. Robustness. When ρs=ρf, the computations with the full schemes are difficult. This is not the case here: both algorithms are stable for very large time steps. Computing speed. The computing time of (6.1) is much larger than that of (6.5) on these meshes. For 50 iterations, it takes 57.4 s on a MacBook pro mid-2012 with (6.1) and only 7.7 s with (6.5), because the linear systems for the latter do not change so we can reuse their factorization. This is done by the library MUMPS in freefem++. 6.5. Conclusion for the two schemes. Because it is much faster and because it is not sensitive to the corner singularity, we conclude that the characteristic-Galerkin algorithm (6.5) is better than (6.1) which uses −u×∇×u. However one must keep in mind that the mathematical analysis of (6.5) is not complete. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 1016 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU Fig. 4.Computation in a three-dimensional geometry with the surface pressure model and a nonconstant bdestined to imitate a stent at t=0.5after 10 time steps. Left: pressure isolines on the surface. Center: color map (mostly between 0.006 (blue) and 0.28 (green)) of |∇ × u|shown on a geometry moved by the physical dblown up by a factor κ=6,solutionof(6.6), and updated from the previous position at each time step. Note that the computational mesh (shown on the left) is not moved while the mesh on the center is deformed only for a graphic display corresponding to the physical reality within the small displacement approximation. Right: isolines of the norm of the (normal) velocity vector on the surface of the geometry moved proportionally to this vector; the effect of the stent is seen. In all three figures the blue color is for zero and the red is for the maximum value. 6.6. A three-dimensional test. From now on we use only (6.5). The aim of the test is to demonstrate that even a fairly complex computation can be done on a standard machine, here an Apple MacBook pro 15 mid-2012 with an intel core i7 at 2.3 MHz and 16 Meg of RAM. The following simulation takes about 10 minutes using freefem++ [37] and medit [28] for the display. The geometry is a pipe coated by a stentlike texture. In practice it means that ˜ b= 100 outside the stent and ˜ b= 1000 on the stent; other structure coefficients in the shell model are all zero (hence it is a surface pressure model). Hecht built the mesh with freefem++ in a quasi-torus with parameters: r =1, R =7, length =10, nonstent region of length 2 at both ends, number of vertices 51448, number of elements 286110, number of boundary elements 26538, number of degrees of freedom of the linear system 1 064122. The other parameters are δt =0.05, ν=0.001, =0.01, no change was observed with =0.001. The pressure difference between the two cross sections as a function of time is p1−p0=6(cos(πt))2. The results are displayed at t=0.5 (after 10 iterations) in Figure 4. On the left the pressure isolines are shown at the surface of the vessel on the fixed mesh used for the computations. In the center the isolines of |∇×u|are shown but the geometry has been deformed by δtκd at each time step by solving (6.6) and with κ= 6. It corresponds to the physical deformation of the geometry computed in the setting of the simplified model but amplified by a factor 6. On the right the isolines of the norm of the velocity vector on the surface are shown on a geometry which is deformed only at t=0.5bymovingall vertices proportionally to u. It is an instantaneous deformation (while in the center it is the physical deformation, exaggerated); at t=0.5 the pressure difference has decreased to zero so the front region shrinks. It is seen also that in the region of the stent, the vessel is slightly dilated and in the stent mesh there is a small deformation too. 7. Summary. By a few minor modifications to the shell model of [45] for FSIs within the small displacement hypothesis, and provided the normal derivative of the pressure is not too large at the compliant boundary (see (2.8)), we have obtained a model which gives numerical results that are similar to the test case of [10] and which can be fully analyzed mathematically in the continuous case and after discretization in time, provided that u×(∇×u) is used for the nonlinear terms. Thus the pres- Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ANALYSIS OF A COUPLED FLUID-STRUCTURE MODEL 1017 sure pis replaced by the dynamic pressure p+ρf 2|u|2. We have also discussed and compared two discretizations by finite elements and concluded that the method which uses characteristic-Galerkin upwinding, though not completely analyzed mathematically (quadrature errors, etc.), is more robust. The method is also computationally not more demanding than a standard Navier–Stokes solver, thus opening the way to computationally viable inverse problems [55, 39, 44]. Acknowledgment. We thank F. Hecht for his very valuable collaboration and for providing the mesh of the pipe with a stent. REFERENCES [1] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault,Vector potentials in threedimensional non-smooth domains, Math. Methods Appl. Sci., 21 (1998), pp. 823–864. [2] M. Astorino and C. Grandmont,Convergence analysis of a projection semi-implicit coupling scheme for fluid-structure interaction problems, Numer. Math., 116 (2010), pp. 721–767. [3] C. B` egue, C. Conca, F. Murat, and O. Pironneau,Les ´equations de Stokes et Navier-Stokes avec des conditions aux limites sur la pression, in Nonlinear Partial Differential Equations and Applications, Coll`egedeFranceSeminar,H.Br´ezis, J.-L. Lions, and D. Cioranescu, eds., vol. IX, 1988, pp. 179–264. [4] D. Boffi, N. Cavallini, and L. Gastaldi,The Finite Element Immersed Boundary Method with Distributed Lagrange Multiplier, SIAM J. Numer. Anal., 53 (2015), pp. 2584–2604. [5] D. Boffi and L. Gastaldi,A finite element approach to the immersed boundary method, Comput. Struct., 81 (2003), pp. 491–501. [6] J. M. Boland and R. Nicolaides,On the stability of bilinear constant velocity-pressure finite elements, Numer. Math., 44 (1984), pp. 219–222. [7] M. Boris, S. ˇ Cani´ c, and M. Bukaˇ c,Fluid-structure interaction in hemodynamics: Modelling, analysis, and numerical simulation, in Fluid-Structure Interaction and Biomedical Applications, T. Bodn´ar, G. Galdi, and ˇ S. Neˇcasov´a, eds., Advances in Mathematical Fluid Mechanics, Birkh¨auser, 2014. [8] K. Boukir, Y. Maday, B. M´ etivet, and E. Razafindrakoto,A high-order characteristics/finite element method for the incompressible Navier-Stokes equations,Internat.J. Numer. Methods Fluids, 25 (1997), pp. 1421–1454. [9] H. Brezis,Analyse fonctionnelle: Th´eorie et Applications, Masson, Paris, 1983. [10] M. Bukaˇ c, S. ˇ Cani´ c, R. Glowinski, J. Tambacac, and A. Quainia,Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement, J. Comput. Phys., 235 (2013), pp. 515–541. [11] M. Bukaˇ c, S. ˇ Cani´ c, and B. Muha,Stability of the Kinematically Coupled β-Scheme for Fluid-Structure Interaction Problems in Hemodynamics, Int. J. Numer. Anal. Model., 12 (2015), pp. 54–80. [12] P. Causin, J. F. Gerbeau, and F. Nobile,Added-mass effect in the design of partitioned algorithms for fluid-structure problems, Comput. Methods Appl. Mech. Engrg., 194 (2005), pp. 42–44. [13] J. Cebral, F. Mut, D. Sforza, R. Lohner, E. Scrivano, P. Lylyk, and C. Putman,Clinical application of image-based CFD for cerebral aneurysms, Internat. J. Numer. Methods Biomed. Eng., 27 (2011), pp. 977–992. [14] T. Chac´ on-Rebollo and R. Lewandowski,Mathematical and Numerical Foundations of Turbulence Models and Applications,Birkh¨auser, Basel, 2014. [15] A. Chambolle, B. Desjardin, M. J. Esteban, and C. Grandmont,Existence of weak solutions for unsteady fluid-plate interaction problem, J. Math. Fluid Mech., 7 (2005), pp. 368– 404. [16] C. Chnafa, S. Mendez, and F. Nicoud,Image-based large-eddy simulation in a realistic left heart, Comput. & Fluids, 94 (2014), pp. 173–187. [17] C. M. Colciago, S. Deparis, and A. Quarteroni,Comparison between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics, J. Comput. Appl. Math., 265 (2014), pp. 120–138. [18] C. Conca, F. Murat, and O. Pironneau,The Stokes and Navier-Stokes equations with boundary conditions on the pressure, Jpn. J. Math., 20 (1994), pp. 297–318. [19] M. Dauge,Neumann and mixed problems on curvilinear polyhedra, Integral Equations Operator Theory, 15 (1992), pp. 227–261. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 1018 CHAC´ ON REBOLLO, GIRAULT, MURAT, AND PIRONNEAU [20] S. Deparis, M. A. Fernandez, and L. Formaggia,Acceleration of a fixed point algorithm for fluid-structure interaction using transpiration conditions, ESAIM Math. Model. Numer. Anal., 37 (2003), pp. 601–616. [21] A. Dervieux,Fluid-Structure Interaction, Innovaive Technology, Kogan Page, London, 2003. [22] T. Fanion, M. Fernandez, and P. Le Tallec,Deriving adequate formulations for fluid structure interaction problems: from ALE to transpiration, Rev. Eur. Elements, 9 (2000), pp. 681–708. [23] C. Farhat, P. Geuzaine, and C. Grandmont,The discrete geometric conservation law and the non-linear stability of ALE schemes for solution of flow problems on moving grids,J. Comput. Phys., 174 (2001), pp. 664–694. [24] M. A. Fernandez, J. Mullaert, and M. Vidrascu,Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures, Comput. Methods Appl. Mech. Engrg., 267 (2013), pp. 566–593. [25] L. Formaggia, J. F. Gerbeau, F. Nobile, and A. Quarteroni,On the coupling of 3Dand1D Navier-Stokes equations for flow problems in compliant vessels, Comput. Methods Appl. Mech. Engrg., 191 (2001), pp. 561–582. [26] L. Formaggia, A. Moura, and F. Nobile,On the stability of the coupling of 3Dand1D fluid-structure interaction models for blood flow simulations, ESAIM Math. Model. Numer. Anal., 41 (2007), pp. 743–769. [27] L. Formaggia, A. Quarteroni, and A. Veneziani,Cardiovascular Mathematics. Modeling and simulation of the circulatory system, Springer MS&A Series, Springer, Milan, Italy, 2009. [28] P. Frey,medit 3.0, OpenGL-Based Scientific Visualization Software,manuscript. [29] V. Girault and P. A. Raviart,Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer Ser. Comput. Math. 5, Springer-Verlag, Berlin, 1986. [30] V. Girault,Incompressible finite element methods for Navier-Stokes equations with nonstandard boundary conditions in R3, Math. Comp., 51 (1988), pp. 55–74. [31] O. Gonzalez and J. C. Simo,On the stability of symplectic and energy-momentum algorithms for non-linear Hamiltonian systems with symmetry, Comput. Methods Appl. Mech. Engrg., 134 (1996), pp. 197–222. [32] K. Gostaff and O. Pironneau,Pressure boundary conditions for blood flows, Chin. Ann. Math. Ser. B, 36 (2015), pp. 829–842. [33] C. Grandmont,Existence of weak solutions for the unsteady interaction of a viscous fluid with an elastic plate, SIAM J. Math. Anal., 40 (2008), pp. 716–737. [34] G. Guidoboni, R. Glowinski, N. Cavallini, S. Lapin, and S. Canic,A kinematically coupled time-splitting scheme for fluid-structure interaction in blood flow, Appl. Math. Lett., 22 (2009), pp. 684–688. [35] G. Guidoboni, M. Guidorzi, and M. Padula,Continuous dependence on initial data in fluid-structure motions, J. Math. Fluid Mech., 14 (2012), pp. 1–32. [36] J. Happel and H. Brenner,Low Reynolds Number Hydrodynamics: With special applications to particulate media, Vol. 1, Springer Science & Business Media, 1983. [37] F. Hecht,New development in freefem++, Numer. Math., 20 (2012), pp. 251–265. [38] A. Hundertmark, M. Luk´ aˇ cov´ a, and S. Neˇ casov´ a,On the weak solution of the fluidstructure interaction problem for shear-dependent fluids, in Recent Developments of Mathematical Fluid Mechanics, H. Amann, Y. Giga, H. Kozono, H. Okamoto, and M. Yamazaki, eds., Adv. Math. Fluid Mech., Birkh¨auser, Basel, 2014, pp. 291–319. [39] L. Iapichino, A. Quarteroni, and G. Rozza,A reduced basis hybrid method for the coupling of parametrized domains represented by fluidic networks, Comput. Methods Appl. Mech. Engrg., 221–222 (2012), pp. 63–82. [40] H. K. Lee, M. A. Olshanskii, and L. G. Rebholz,On error analysis for the 3DNavier– Stokes equations in velocity-vorticity-helicity form, SIAM J. Numer. Anal., 49 (2011), pp. 711–732. [41] P. LeTallec and J. Mouro,Fluid structure interaction with large structural displacements, Comput. Methods Appl. Mech. Engrg., 190 (2001), pp. 3039–3067. [42] Y. Liu and Y. Mori,Lpconvergence of the immersed boundary method for stationary Stokes problems, SIAM J. Numer. Anal., 52 (2014), pp. 496–514. [43] M. Lukacova-Medvidov, G. Rusnakov, and A. Hundertmark-Zauskov,Kinematic splitting algorithm for fluid-structure interaction in hemodynamics, Comput. Methods Appl. Mech. Engrg., 265 (2013), pp. 83–106. [44] A. Manzoni, A. Quarteroni, and G. Rozza,Shape optimization for viscous flows by reduced basis methods and free-form deformation, Internat. J. Numer. Methods Fluids, 70 (2012), pp. 646–712. Downloaded 10/20/16 to 150.214.182.15. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php