A quasi-Lagrangian finite element method for the Navier-Stokes equations in a time-dependent domainThis work has been supported by the Russian Science Foundation (RSF) grant 14-31-00024.

# A quasi-Lagrangian finite element method for the Navier-Stokes equations in a time-dependent domain††thanks: This work has been supported by the Russian Science Foundation (RSF) grant 14-31-00024.

Alexander Lozovskiy Institute of Numerical Mathematics RAS; saiya-jin@yandex.ru    Maxim A. Olshanskii Department of Mathematics, University of Houston; molshan@math.uh.edu    Yuri V. Vassilevski Institute of Numerical Mathematics RAS, Moscow Institute of Physics and Technology and First Moscow State Medical University; yuri.vassilevski@gmail.com
###### Abstract

The paper develops a finite element method for the Navier-Stokes equations of incompressible viscous fluid in a time-dependent domain. The method builds on a quasi-Lagrangian formulation of the problem. The paper provides stability and convergence analysis of the fully discrete (finite-difference in time and finite-element in space) method. The analysis does not assume any CFL time-step restriction, it rather needs mild conditions of the form , where depends only on problem data, and , is polynomial degree of velocity finite element space. Both conditions result from a numerical treatment of practically important non-homogeneous boundary conditions. The theoretically predicted convergence rate is confirmed by a set of numerical experiments. Further we apply the method to simulate a flow in a simplified model of the left ventricle of a human heart, where the ventricle wall dynamics is reconstructed from a sequence of contrast enhanced Computed Tomography images.

## 1 Introduction

Fluid flows in time-dependent domains are ubiquitous in nature and engineering. In many cases, finding the domain evolution is part of the problem and the mathematical model couples fluid and structure dynamics. Examples include fluid–structure interaction problems for blood flow in compliant vessels, flows around turbine blades or fish locomotion. In other situations, one may assume that the motion of the domain is given and one has to recover the induced fluid flow. One example of a problem, which often assumes a priori information about the flow domains evolution, is the blood flow simulation in a human heart when the (patient-specific) motion of the heart walls is recovered from a sequence of medical images [1, 2, 3, 4, 5, 6, 7, 8]. Nowadays numerical simulations are commonly used to understand fluid dynamics and predict statistics of practical interest in this and other applications. In the present paper, we develop a finite element (FE) method for a quasi-Lagrangian formulation of the incompressible Navier-Stokes equations in a moving domain. We consider an implicit–explicit method, i.e an implicit method with advection field in the inertia term lagged in time. For the spatial discretization we employ inf-sup stable pressure–velocity elements.

Several techniques have been introduced in the literature to overcome numerical difficulties due to the evolution of the domain. This includes space–time finite element formulations, immersed boundary methods, level-set method, fictitious domain method, unfitted finite elements, and arbitrary Lagrangian–Eulerian (ALE) formulation, see, e.g., [9, 10, 11, 12, 13, 14, 15, 16, 17]. In this paper we analyze a finite element method based on a quasi-Lagrangian formulation of the equations in the reference domain. Related analysis of finite element methods for parabolic or fluid equations in moving domains can be found in several places in the literature. We note that well-posedness of space-time weak saddle-point formulations of the (Navier–)Stokes equations is a subtle question, see the recent treatment in [18] for the case of a steady domain. A rigorous stability and convergence analysis of space–time (FE) methods for fluid problems seems to be largely lacking. Scalar problems have been understood much better; for example, a space–time discontinuous FE method for advection–diffusion problems on time-dependent domains was analyzed in [19]. ALE and Lagrangian finite element methods are more amenable to analysis. The stability of ALE finite element methods for parabolic evolution problems was treated in [16]. The authors of [20] analyzed the convergence of a finite element ALE method for the Stokes equations in a time-dependent domain when the motion of the domain is given. The analysis [20] imposes time step restriction, assumes zero velocity boundary condition and certain smoothness assumptions for the finite element displacement field. A closely related method to the one studied here was considered in [21]. However, that paper introduced an assumption that a divergence free extension of a boundary condition function to the computational domain is given. This assumption is not always practical and the present paper avoids it. Moreover, this paper develops error analysis, while the thrust of [21] was the stability analysis and the numerical recovery procedure of the domain motion from medical images.

In the present paper, we analyze a quasi-Lagrangian FE formulation that is closely related to an ALE formulation, although they are not equivalent. In the present approach we discretize equations in a reference time-independent domain. The geometry evolution is accounted in time-dependent coefficients. Inertia terms are further linearized so that only a system of linear algebraic equations is solved on each time step. We consider practically relevant boundary conditions, which result in non-homogeneous velocity on the boundary. For this method we prove numerical stability and optimal order error estimates in the energy norm without a CFL condition on the time step. Divergence-free condition enforced in the reference domain leads to time dependent functional spaces; this and handling non-homogeneous boundary conditions are two main difficulties that we overcome in the analysis. For the numerical stability bound we shall need the condition on the space mesh size and time step of the form , where is polynomial degree of velocity finite element space and is a constant. We note that if one assumes zero boundary conditions for velocity, which is a standard assumption for FE stability bounds in steady domains, then the results of the paper hold without the above conditions on and . In our opinion, homogeneous boundary conditions is not a suitable assumption for the FE analysis in evolving domains, see discussion in section 3.

Thus the paper advances the known analysis by including inertia effect, removing CFL time-step restriction, handling physically meaningful boundary conditions, and making no further auxiliary assumptions except the following one: the domain evolution is given a priori by a smooth mapping from a reference domain to a physical domain and exact quadrature rules are applied in the reference domain, i.e. we do not analyse possible errors due to inexact numerical integration. The mapping is not necessarily Lagrangian in the internal points, but it has to be Lagrangian for those parts of the boundary, where correct tangential velocity boundary values are important. Theoretical results are illustrated numerically for an example with a synthetic known solution. We further illustrate the performance of the numerical method by applying it to blood flow simulation in a simplified model of the human left ventricle. The domain motion in this example is reconstructed from a sequence of ceCT images of a real patient heart over one cardiac cycle. The reconstruction procedure is described in detail in our preceding paper [21].

The remainder of this paper is organized as follows. In section 2 we review the mathematical model, including governing equations and boundary conditions, and some useful results for this model found in the literature. We recall the energy balance satisfied by smooth solutions. A suitable weak formulation is introduced. Based on the weak formulation, in section 3 we introduce the finite element method. Non-homogeneous boundary conditions are interpolated numerically. Energy stability estimate for the finite element method is shown in section 4. Optimal order error bound for the method is demonstrated in section 5. Section 6 collects results of numerical experiments. Some closing remarks can be found in the summary and outlook section 7.

## 2 Mathematical model

Consider a time-dependent domain , , occupied by fluid. To formulate a flow problem, we introduce the reference domain and a mapping from the space–time cylinder to the physical domain,

 \boldmathξ\unboldmath : Q→Qphys:=⋃t∈[0,T]Ω(t)×{t}.

The mapping is assumed to be level-preserving, i.e. for all . We assume also that the evolution of is sufficiently smooth such that . Denote the spatial gradient matrix of by , and . Furthermore, we assume that there exist such positive reals that

 (1)

The dynamics of incompressible Newtonian fluid can be described in terms of the velocity vector field and the pressure function defined in for . This paper studies a finite element method for fluid equations formulated in the reference domain. For , defined in , the fluid dynamics is given by the following set of equations:

 (2)

with body forces and the initial condition . We assume the fluid to be Newtonian, with the kinematic viscosity parameter . The constitutive relation in the reference domain reads

 ˆ\boldmathσ\unboldmath∘\boldmathξ% \unboldmath=−pI+ν(∇uF−1+F−T(∇u)T)   in Q. (3)

### 2.1 Boundary conditions

We distinguish between the no-slip , Dirichlet and outflow parts of the boundary, and . On we impose no-penetration no-slip boundary condition, i.e. the fluid velocity on is equal to the material velocity of the boundary (see the discussion below),

 ^u=\boldmathξ\unboldmatht∘\boldmathξ% \unboldmath−1on ∂Ωns(t), (4)

while on and we prescribe Dirichlet and Neumann conditions,

 ^u=^uDon ∂ΩD(t),^\boldmathσ\unboldmath^n=^gon ∂ΩN(t). (5)

Here is a given velocity and is the exterior unit normal vector on . If for some we assume .

In the reference domain, define , , . We assume that , , are independent of .

###### Remark 1.

The normal velocity of the boundary is . However, the material tangential velocity of the boundary is defined by the tangential part of only if is the Lagrangian mapping, i.e. , , defines the material trajectory for (or at least for ). In some applications such Lagrangian mapping is not available, and in this case (4) may produce spurious tangential velocities on the boundary. For example, this may happen if is reconstructed from medical images. Thus, in practice one may or may not amend (4) based on any additional information about the tangential motions for a better model.

### 2.2 An extension result

The solvability of the problem (2)–(3) and the existence of its weak solutions is treated, for example, in [22]. Moreover, it is shown in  [22] that for smoothly evolving the mapping can be chosen in such a way that depends only on . From numerical viewpoint, such a mapping may not be practically available, and so we allow to vary in time and space. However, the following corollary of this result is important for us, see Theorem 4.4. in  [22]: Assume , then there exists such that on and in for . The condition is satisfied in the case of for all . Indeed, the Reynolds transport theorem and the incompressibility assumption for fluid imply

 ddt|Ω(t)|=ddt∫Ω(t)dx=∫∂Ω(t)vΓds=∫∂Ω(t)^n⋅^uds=∫Ω(t)div^udx=0.

For the finite element analysis in this paper we assume that can be taken -smooth. We define smooth function that satisfies

 v1∈C3(Q)d,div(JF−1v1)=0 in Ω0,v1=\boldmathξ% \unboldmatht on ∂Ω0. (6)

We stress that we need the result about existence of for the finite element analysis, but one never needs to know or compute for the implementation of the FE method.

### 2.3 Energy equality

In this section, we assume no-penetration no-slip boundary condition (4) imposed on the whole boundary, i.e. . By we denote the scalar product, and denotes the norm. For vector fields and tensor fields , we use the same notation to denote and , and obviously , . We shall also make use of the identity for all , :

 (w⋅∇u,v)+12((divw)u,v)=12((w⋅∇u,v)−(w⋅∇v,u))+12∫∂Ω0(n⋅w)uvds. (7)

We multiply the first equality in (2) by , integrate it over the reference domain, and employ (7) for integration by parts. We get

here is the exterior unit normal vector on . The mass balance yields the equality

 Jt+div(JF−1(u−\boldmathξ% \unboldmatht))=0in Q. (8)

This identity leads to some cancellations and we get

 12ddt∥J12u∥2+(J(^% \boldmathσ\unboldmath∘\boldmathξ\unboldmath)F−T,∇u)−∫∂Ω0(J(^% \boldmathσ\unboldmath∘\boldmathξ\unboldmath)F−Tn)⋅\boldmathξ\unboldmathtds=(Jf,u).

The Piola identity, , implies the following equality

 div(JF−1u)=J(∇u):F−Tin Q, (9)

where . Using the notation for the rate of deformation tensor in the reference coordinates, we get with the help of (9) and the second equation in (2)

 (J(^\boldmathσ\unboldmath∘\boldmathξ\unboldmath)F−T,∇u)=(J(−pI+ν(∇uF−1+F−T(∇u)T))F−T,∇u)=2ν(JDξ(u)F−T,∇u)=2ν(JDξ(u),∇uF−1)=2ν(JDξ(u),Dξ(u)).

In the last equality we used that for any symmetric tensor and any tensor , it holds . Therefore, the energy balance equality in reference coordinates takes the form

 12ddt∥J12u∥2+2ν∥J12Dξ(u)∥2−∫∂Ω0(J(^\boldmathσ\unboldmath∘\boldmathξ\unboldmath)F−Tn)⋅\boldmathξ\unboldmathtds=(Jf,u). (10)

The mechanical interpretation of (10) is the following one: the work of external forces (right-hand side) is balanced by the change of kinetic energy (the first term), viscous dissipation of energy (the second term), and flow intensification due to the boundary condition (the third term).

### 2.4 Weak formulation

For we introduce the following time-dependent trilinear and bilinear forms:

 c(\boldmathξ\unboldmath;w,u,% \boldmathψ\unboldmath) a(\boldmathξ\unboldmath;u,\boldmathψ\unboldmath) =∫Ω02νJDu:D% \boldmathψ\unboldmathdx,u,% \boldmathψ\unboldmath∈H1(Ω0)d, b(\boldmathξ\unboldmath;p,\boldmathψ% \unboldmath) =∫Ω0pJF−T:∇\boldmathψ% \unboldmathdx,p∈L2(Ω0), % \boldmathψ\unboldmath∈H1(Ω0)d.

The weak formulation of (2)–(3) reads: Find satisfying on , on and

 (Jut,\boldmathψ\unboldmath)+c(% \boldmathξ\unboldmath;u−\boldmathξ\unboldmatht,v,\boldmathψ\unboldmath)+a(\boldmathξ% \unboldmath;u,\boldmathψ\unboldmath)−b(\boldmath% ξ\unboldmath;p,\boldmathψ\unboldmath)+b(\boldmathξ\unboldmath;q,u)=(Jf,\boldmathψ\unboldmath)+∫∂ΩN0Jg⋅\boldmathψ% \unboldmathds (11)

for all on , for all .

## 3 Discretization method

In this section we introduce both time and space discretizations of the formulation (2) in the reference domain. Treating the flow problem in reference coordinates allows us to avoid triangulations and finite element function spaces dependent on time. In this paper, we assume that the mapping is given explicitly and used in the finite element formulation without any further numerical approximation apart from the boundary condition.

Let a collection of simplices (triangles for and tetrahedra for ) form a consistent regular triangulation of the reference domain . We let . Consider conforming FE spaces and ; is a subspace of of functions vanishing on . We assume that and form the LBB-stable finite element pair: There exists a mesh-independent constant , such that

 infqh∈Qhsupvh∈V0h(qh,divvh)∥∇vh∥∥qh∥≥c0>0. (12)

As an example of admissible discretization, we consider the generalized Taylor-Hood finite element spaces,

 Vh ={uh∈C(Ω0)d:uh|T∈[Pm+1(T)]d,∀ T∈Th}, (13) Qh ={qh∈C(Ω0):qh|T∈Pm(T),∀ T∈Th},

where integer is polynomial degree.

Assuming a constant time step , we use the notation , and similar for and . To emphasize the dependence on , denote , , .

For given spatial functions , , denotes the backward finite difference at . For a sufficiently smooth vector function , denote by its nodal Lagrange interpolant.

Let . The finite element discretization of (11) reads: For , find satisfying on , on and the following equations

 (Jk−1[uh]kt,\boldmathψ\unboldmathh)+(12[J]ktukh,\boldmathψ\unboldmathh)+12(div(JkF−1kwkh)ukh,\boldmathψ\unboldmathh)+c(\boldmathξ\unboldmathk;wkh,ukh,\boldmathψ\unboldmathh)+a(\boldmathξ\unboldmathk;ukh,\boldmathψ\unboldmathh)−b(\boldmathξ\unboldmathk;pkh,\boldmathψ% \unboldmathh)+b(\boldmathξ\unboldmathk;qh,ukh)=(Jkfk,\boldmathψ\unboldmathh)+∫∂ΩN0Jkgk⋅\boldmathψ% \unboldmathds (14)

for all with advection velocity .

The second and the third terms in (14) are consistent due to the identity (8) and are added in the FE formulation to enforce the conservation property of the discretization. While our computations show that in practice this term can be skipped, we need these terms for the stability bound in the next section. In the numerical analysis of incompressible Navier-Stokes equations in the Eulerian description, including these terms corresponds to the Temam’s skew-symmetric form of the convective terms [23].

Note that the inertia terms are linearized so that a linear algebraic system should be solved on each time step. In the next section we show that the finite element method is energy stable.

Note that if , then for the boundary condition should be compatible with divergence constraint. Therefore, in this case we let and let on and on .

###### Remark 2.

Condition on those parts of the domain, where no-slip and no-penetration takes place, of an evolving fluid domain is not physically reasonable. Moreover, it leads to strong simplifications of finite element analysis. Indeed, we shall see that constructing suitable FE extensions of boundary conditions to computational domain is not straightforward.

## 4 Stability of FEM solution

From now on we assume for all . To show the stability, we need several preparatory steps which help us to handle non-homogeneous boundary conditions and time-dependent bilinear forms. First we note the Korn’s-type inequality in the reference domain,

 ∥∇u∥≤CK∥J12Dξ(u)∥∀ u∈H10(Ω0)d,t∈[0,T], (15)

with uniformly bounded with respect to . The estimate (15) easily follows from the standard Korn inequality and assumptions in (1), see [21]. Thanks to (15) the bilinear form is coercive on uniformly in time. Also due to (1) the bilinear forms and are continuous uniformly in time on and , respectively.

Unlike the continuous case, the finite element solution does not satisfy a strong formulation and so the arguments from section 2.3 do not apply directly. To show the proper energy balance for the finite element solution, we split it into a part vanishing on the boundary and another a priori defined (and so stable) part, which satisfy correct boundary conditions. Therefore, we consider decomposition , such that

 vkh,1=ukh  on ∂Ω0×[0,T],b(\boldmathξ\unboldmathk;qh,vh,1)=0  ∀ qh∈Qh,  k=1,2,…,N, (16)

and

 ∥vkh,1∥W1,∞≤C,(Jk−1(vkh,1−vk−1h,1),vkh)≤CΔt(∥vkh∥+hm+2∥[vh]kt∥), (17)

with some real depending only on data and independent of , . The existence of such decomposition will be explicitly demonstrated in the next section.

With the help of this decomposition, the finite element method can be re-formulated as follows: find satisfying for all

 (Jk−1[vh]kt,\boldmathψ\unboldmathh)+(12[J]ktvkh,\boldmathψ\unboldmathh)+12(div(JkF−1kwkh)vkh,\boldmathψ\unboldmathh)+c(\boldmathξ\unboldmathk;wkh,vkh,\boldmathψ\unboldmathh)+c(\boldmathξ\unboldmathk;vkh,vkh,1,\boldmathψ\unboldmathh)+a(\boldmathξ\unboldmathk;vkh,\boldmathψ\unboldmathh)−b(\boldmathξ\unboldmathk;pkh,% \boldmathψ\unboldmathh)+b(\boldmathξ\unboldmathk;qh,vkh)=⟨˜fkh,\boldmathψ\unboldmathh⟩ (18)

with . We stress that the formulation (18) appears here only for the purpose of analysis. Although equivalent to (14), the formulation (18) is not practical, since it requires an explicit knowledge of . For the analysis, it is sufficient to know that exists.

We test (18) with , . We handle each resulting term separately and start with the first term in (18):

 (Jk−1[vh]kt,vkh)=12Δt(∥J12kvkh∥2−∥J12k−1vk−1h∥2)−12([J]ktvkh,vkh)+Δt2∥J12k−1[vh]kt∥2. (19)

The term in (19) cancels with the second term in (18). Applying (7) to the fourth (inertia) term in (18) and using boundary conditions give

 (Jk(∇vkhF−1kwkh),vkh)=−12(div(JkF−1kwkh)vkh,vkh). (20)

This cancels with the third term in (18). We keep the fifth term as it is. The sixth term in (18) gives

 a(\boldmathξ\unboldmathk;vkh,vkh)=2ν(JkDk(vkh),Dk(vkh))=2ν∥∥∥J12kDk(vkh)∥∥∥2.

The -terms cancel out for . Substituting all equalities back into (18), we obtain the following energy balance for the -part of finite element solution :

 12Δt(∥J12kvkh∥2−∥J12k−1vk−1h∥2)+2ν∥∥∥J12kDk(vkh)∥∥∥2+Δt2∥∥∥J12k−1[vh]kt∥∥∥2+(Jk(∇vkh,1F−1k)vkh,vkh)=⟨˜fkh,vkh⟩. (21)

We deduce an energy stability estimate for the finite element method from the balance in (21) and a priori estimates in (17). For the sake of notation, we introduce , which defines a -dependent norm uniformly equivalent to the -norm. Using estimates (17) and the definition of one shows that the forcing term is bounded,

with a constant depending only on problem data. We substitute this in (21) and further use the bound . This yields the estimate

 12Δt(∥vkh∥2k−∥vk−1h∥2k−1)+ν∥∥Dk(vkh)∥∥2k+(Jk(∇vkh,1F−1k)vkh,vkh)≤C(1+(Δt)−1h2m+4)), (22)

where the constant depends on the problem data, but not on and .

Thanks to Sobolev’s embedding inequalities as well as (1) and (15), we bound the fourth term in (22) resulting from the boundary motion in two ways,

If the factor is not too large such that it holds

 C1∥∇vkh,1∥≤ν/2, (23)

then the intensification term can be absorbed by the viscous dissipation term. So we obtain from (22)

 12∥vkh∥2k+νΔt∥Dk(vkh)∥2k≤12∥vk−1h∥2k−1+CΔt(1+(Δt)−1h2m+4)), (24)

with some depending only on problem data.

Summing up inequalities (24) over , , gives

 (25)

Otherwise, if (23) does not hold, we estimate

 12(1−2C2Δt)∥vkh∥2k+νΔt∥Dk(vkh)∥2k≤12∥vk−1h∥2k−1+CΔt(1+(Δt)−1h2m+4)). (26)

Now we assume that is sufficiently small such that . Summing over , gives

 12∥vnh∥2n+νn∑k=1Δt∥Dk(vkh)∥2k≤C2n∑k=1Δt∥vkh∥2k+12∥v0∥20+C(1+(Δt)−1h2m+4). (27)

Applying discrete Gronwall’s inequality yields

 (28)

Finally, we apply the triangle inequality and (17) one more time to show the stability bound for the velocity solution to (14),

 (29)

where constant depends on the problem data, i.e. , , , , but not on and .

Note that the restriction has resulted from handling non-homogeneous boundary conditions. Indeed, one gets the extra -dependent term in (17) while estimating the difference . If one defines as a nodal interpolant to the smooth divergence-free function then the estimate as in (17) follows without this extra term. However, the nodal interpolant does not satisfy the “divergence-free” condition in (16) and is defined as a suitable projection of .

## 5 Convergence analysis

In this section, we demonstrate optimal order convergence of to . From now on, we assume that is a convex polyhedron (polygon if ).

### 5.1 Preliminaries

For the Navier-Stokes equations solution in the reference domain we recall the notations . The finite element solution to (14) is . Furthermore,

 {ek,ek}:={uk−ukh,pk−pkh}

denotes the finite element error.

For the fixed time instance , we define the subspace of ,

 Xkh:={\boldmathψ\unboldmathh∈Vh:b(\boldmathξ\unboldmathk;qh,\boldmathψ% \unboldmathh)=0  ∀ qh∈Qh}.

In Lemma 5.1 below we prove an important technical result: if the pair of spaces is inf-sup stable in the sense of (12), and satisfies certain assumptions, then the subspace has full approximation properties. Here and further we assume that for all the domains are such that the Stokes problem in the physical domain satisfies -regularity property: for any the unique solution of

 −Δ^\boldmathϕ\unboldmath+∇^λ=f,div^\boldmathϕ\unboldmath=0  in Ω(t),^\boldmathϕ\unboldmath=0  on ∂Ω(t) (30)

satisfies , and

 ∥^\boldmathϕ\unboldmath∥H2(Ω(t))+∥^λ∥H1(Ω(t))≤CR∥f∥L2(Ω(t)), (31)

with some uniformly bounded for .

We recall that denotes the nodal interpolant of and .

###### Lemma 5.1.

Assume (12) and

 supQ∥I−F∥F≤ε, (32)

with sufficiently small and the identity matrix . Then for satisfying , , we have

 (33)

where is the polynomial degree from (13) and depends on , but not on , or .

###### Proof.

From (32) it follows that

 supQ∥I−F−1∥F≤ε1−ε,supQ|1−J|≤cε (34)

holds with some depending only on space dimension . The inf-sup condition (12) implies that for a given there exists such that and . Thanks to (34) we estimate

 b(\boldmathξ\unboldmathk;qh,\boldmathψ\unboldmathh) =