Analysis of a velocity–stress–pressure formulation for a fluid–structure interaction problem

Applied Numerical Mathematics 123 (2018) 275–299

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Analysis of a velocity–stress–pressure formulation for a fluid–structure interaction problem María González a,∗ , Virginia Selgas b a b

Departamento de Matemáticas, Universidade da Coruña, Campus de Elviña s/n, 15071, A Coruña, Spain Departamento de Matemáticas, Universidad de Oviedo, Escuela Politécnica de Ingeniería, 33203, Gijón, Spain

a r t i c l e

i n f o

Article history: Received 29 November 2016 Received in revised form 14 September 2017 Accepted 23 September 2017 Available online 28 September 2017 Keywords: Fluid–structure interaction Stokes equations Incompressible fluid Linear elastodynamics Finite element Semidiscrete scheme

a b s t r a c t We consider a fluid–structure interaction problem consisting of the time-dependent Stokes equations in the fluid domain coupled with the equations of linear elastodynamics in the solid domain. For simplicity, all changes of geometry are neglected. We propose a new method in terms of the fluid velocity, the fluid pressure, the structural velocity and the Cauchy stress tensor. We show that the new weak formulation is well-posed. Then, we propose a new semidiscrete problem where the velocities and the fluid pressure are approximated using a stable pair for the Stokes problem in the fluid domain and compatible finite elements in the solid domain. We obtain a priori estimates for the solution of the semidiscrete problem, prove the convergence of these solutions to the solution of the weak formulation and obtain error estimates. A time discretization based on the backward Euler method leads to a fully discrete scheme in which the computation of the approximated Cauchy stress tensor can be decoupled from that of the remaining unknowns at each time step. The displacements in the structure (if needed) can be recovered by quadrature. Finally, some numerical experiments showing the performance of the method are provided. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Fluid–structure interaction problems appear in many branches of science and engineering; see, for instance, the books [4,24,5]. There are different mathematical models to study fluid–structure interactions; for a classification depending on the fluid model, see [11]. In this paper, we consider the interaction of a viscous incompressible fluid and a linear elastic structure in a two- or three—dimensional bounded domain. The motion of the fluid is given by the time-dependent Stokes equations while the motion of the structure is described by the equations of linear elastodynamics. The coupling is expressed through the equilibrium of surface forces and the continuity of velocities through a fixed interface. This kind of model can be used when the structure undergoes only infinitesimal elastic displacements but the structural velocity is large enough so that the fluid and the structure remain fully coupled. This setting arises, for instance, when the elastic structure is subject to high frequency, small displacement oscillations; we refer to [11] for more details on the physical validity of the model. This model has been studied also in [23,10–12,3,25,9,15]; other similar models can be found in [22,14]. In [23] the eigenmodes of the coupled system considered here were analyzed, while its homogenization was developed in [10]. In [11], the authors define a divergence-free weak formulation that does not involve the fluid pressure and prove the existence

*

Corresponding author. E-mail addresses: [email protected], [email protected] (M. González), [email protected] (V. Selgas).

https://doi.org/10.1016/j.apnum.2017.09.011 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.

276

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

F

F



-

S

S

n

S

 S

F n

S w

F

 S

F nS o

Fig. 1. Domain of the problem. Cases  F = ∅ =  S (above),  F = ∅ (below, left) and  S = ∅ (below, right).

and uniqueness of the fluid velocity and solid displacement. Under some additional smoothness assumptions on the data, they proved a regularity result for the weak solution and the existence of a square-integrable pressure field. In their later work [12], the authors define semidiscrete finite element approximations, prove the convergence of finite element solutions and derive error estimates for the finite element approximations. In [3] an error analysis of a full discretization of the linear coupled problem is performed. A finite element approximation is considered for the space discretization while a semi-implicit time stepping strategy is used for the time discretization. On the other hand, the motion of a Stokes fluid through an elastic half cylinder with thickness, where the fluid is driven by a small time-dependent pressure drop between the outflow and the inflow sections is analyzed in [25]. There, existence and uniqueness of a weak solution is proved and, under appropriate regularity assumptions on the data, strong energy estimates and the existence of a pressure are shown. The model considered in [25] can be viewed as a first step in modeling dynamic blood flow in an arterial segment. In [9] an extension of the Nitsche method to fluid–structure interaction problems on unfitted meshes is considered. Optimal convergence in energy norm is obtained for smooth solutions for the semi-discretization in space. Besides, [15] deals with the convergence analysis of generalized Robin–Neumann explicit coupling schemes. In [11,12,25] a global velocity together with the corresponding initial condition are defined in order to prove existence and uniqueness. In this work, we follow [20] and propose a new formulation of the problem in terms of the global velocity, together with the fluid pressure and the structure stress tensor. Displacements in the structure (if needed) can be recovered by means of a quadrature formula. Let us emphasize that our new weak formulation has a non-standard structure. We show that it has a unique solution. Then, making use of a stable pair for the Stokes problem in the fluid domain and compatible finite elements in the solid domain, we propose a semidiscrete problem for which we obtain a priori estimates of its solution, prove the convergence to the solution of the weak formulation and establish a priori error estimates. Finally, using the backward Euler method for the time discretization, we arrive at a fully discrete scheme for which the computation of the structure stress tensor can be decoupled in each iteration of that of the remaining unknowns. In order to compute approximations of the global velocity and the pressure field, a saddle-point problem has to be solved. We remark that this new approach is easy to implement as compared to other methods available in the literature. Moreover, the numerical experiments we carried out confirm the robustness and good convergence properties of the proposed scheme. The rest of the paper is organized as follows. In Section 2, we describe the model problem. In Section 3, we derive a new variational formulation in terms of the fluid velocity, the fluid pressure, the structural velocity and the structure Cauchy stress tensor, and show that it is well-posed. Then, in Section 4 we propose and analyze a semidiscrete scheme using a stable pair for the Stokes problem in the fluid domain and compatible finite elements in the solid domain; we obtain a priori estimates for the solution of the semidiscrete problem, prove its convergence to the solution of the weak formulation and derive a priori error estimates. Finally, in Section 5 we deduce a fully discrete scheme, where the computation of an approximation of the Cauchy stress tensor can be decoupled from that of the remaining unknowns. We also provide some examples of finite element subspaces verifying the hypotheses made in the analysis. Finally, we illustrate the good behavior of our method by providing some numerical results. Throughout the paper, we use the standard notations for Sobolev and Bochner spaces and the corresponding norms; see, for instance, [1,21]. 2. The model problem Let d ∈ {2, 3} be the spatial dimension,  F ⊂ Rd be the fluid region and  S ⊂ Rd be the solid region. We assume that  F and  S are bounded Lipschitz domains such that  F ∩  S = ∅. For simplicity, we assume that  := ∂  F ∩ ∂  S , the interface between the fluid and the structure, is fixed and with meas() > 0. We also assume that meas( F ) > 0 and/or meas( S ) > 0, where  F := ∂  F \  and  S := ∂  S \  are the remainder parts of the fluid and solid boundaries, respectively. Finally, we denote by  :=  F ∪  ∪  S the entire domain.

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

277

We remark that the method proposed in this paper can be applied in two and three spatial dimensions with any of the following configurations: the fluid and the solid can be adjacent (meas( F ) > 0 and meas( S ) > 0), the fluid can be contained in the solid ( F = ∅) or the solid can be surrounded by the fluid ( S = ∅); see Fig. 1. In the sequel, we denote by I the identity matrix and, given a sufficiently smooth vector field w, we denote by ε (w) := 12 ∇ w + (∇ w)t the corresponding strain tensor of small deformations. In the fluid domain we consider the time-dependent Stokes equations:

⎧ ρ F ∂t v F − div(σ F ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ σF ⎪ ⎨ div(v F ) ⎪ ⎪ ⎪ ⎪ vF ⎪ ⎪ ⎪ ⎩ v F (0)

= fF

in (0, T ) ×  F

= − p I + 2 ν ε (v F ) in (0, T ) ×  F = 0

in (0, T ) ×  F

= 0

on (0, T ) ×  F

=

in  F

v0F

(1)

where ρ F is the fluid density and ν is the kinematic viscosity of the fluid (both assumed to be constant), and T > 0 denotes the final time. Besides, f F is a given body force and v0F the initial fluid velocity. Further, the unknowns v F , σ F and p represent the fluid velocity, the fluid stress tensor and the fluid pressure, respectively. On the other hand, in the solid domain we consider the equations of linear elastodynamics:

⎧ ρ S ∂tt u S − div(σ S ) = f S in (0, T ) ×  S ⎪ ⎪ ⎪ ⎪ ⎨ σ S = λ div(u S ) I + 2 μ ε (u S ) in (0, T ) ×  S ⎪ ⎪ ⎪ ⎪ ⎩

uS = 0 u S (0) = u0S ,

on (0, T ) ×  S

∂t u S (0) = v0S

(2)

in  S

where ρ S is the solid density (assumed constant), and λ and μ are the Lamé parameters. Besides, f S is the given loading force, and u0S , v0S are the initial displacement and initial velocity of the solid, respectively. The unknowns u S and σ S stand for the displacement and the Cauchy stress tensor, given by Hooke’s law. Notice that we suppose, without loss of generality, that the solid is fixed on  S . Finally, we assume that the velocity and the normal stresses are continuous across :

v F = ∂t u S ,

σ F n = σ S n

on (0, T ) ×  ,

(3)

where n is the unit normal vector to  pointing into  S (see Fig. 1). For compatibility reasons, we also suppose that the initial velocities satisfy

v0F = v0S

on  .

(4)

Typically, the fluid stress tensor and the Cauchy stress tensor are substituted in the corresponding equilibrium equations, and problem (1)–(3) is formulated in terms of the fluid velocity v F , the fluid pressure p and the solid displacement u S (see [11,12,25]). However, doing so, the continuity of the velocity through the interface  (see (3)) produces some difficulties in the analysis and the numerical solution of the corresponding discrete scheme. In this work, we follow [20] and formulate the problem in the solid domain in terms of the structural velocity v S := ∂t u S (instead of the displacement u S ) and the Cauchy stress tensor σ S . To this end, we take the time derivative in the constitutive law and the Dirichlet boundary condition, and provide initial data for the new unknowns. Then, the equations of linear elastodynamics can be written as a first order system in time:

⎧ ρ S ∂t v S − div(σ S ) = f S in (0, T ) ×  S ⎪ ⎪ ⎪ ⎪ ⎨ ∂t σ S = λ div(v S ) I + 2 μ ε (v S ) in (0, T ) ×  S ⎪ ⎪ ⎪ ⎪ ⎩

vS = 0 v S (0) =

v0S

,

where we define σ := can be recovered using that

σ S (0) = σ

λ div(u0S ) I + 2

0 S

on (0, T ) ×  S

με

0 S

(u0S ),

in  S assuming u0S is smooth enough. Notice that the solid displacement (if needed)

t u S (t ) = u0S +

v S ( z) dz .

(5)

0

Substituting the fluid stress tensor in the momentum equation in  F and rewriting the transmission conditions on , we obtain an equivalent formulation of the model problem in terms of the fluid velocity v F , the fluid pressure p, the structural velocity v S and the Cauchy stress tensor σ S as follows:

278

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

⎧ ρ F ∂t v F + ∇ p − 2 ν div(ε (v F )) ⎪ ⎪ ⎪ ⎪ ⎪ div(v F ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ρ S ∂t v S − div(σ S ) ⎪ ⎪ ⎪ ⎪ ⎪ ∂t σ S − λ div(v S ) I − 2 μ ε (v S ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ vF ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

= fF

in (0, T ) ×  F

= 0

in (0, T ) ×  F

= fS

in (0, T ) ×  S

= 0

in (0, T ) ×  S

= vS

on (0, T ) × 

(− p I + 2 ν ε (v F )) n = σ S n vF = 0

on (0, T ) ×  F

vS = 0

on (0, T ) ×  S

v F (0) = v0F v S (0) =

v0S

in  F

σ S (0) = σ

,

(6)

on (0, T ) × 

0 S

in  S

Notice that problems (1)–(3) and (6) are equivalent in the sense of the following proposition. Proposition 1. If (v F , p , u S ) is a solution to problem (1)–(3), then (v F , p , v S , σ S ), with v S := ∂t u S and σ S := λ div(u S ) I + 2 μ ε (u S ), is a solution to problem (6). Conversely, if (v F , p , v S , σ S ) is a solution to problem (6), then (v F , p , u S ), with u S given by (5), is a solution to problem (1)–(3). 3. Weak formulation Let us work formally with the system of equations (6) for a.e. t ∈ (0, T ). Multiplying the first equation in (6) by a test function w :  F → Rd such that w = 0 on  F , and integrating by parts, we obtain



d

ρF





v F (t ) · w −

dt F



−

p (t ) div(w) + 2ν F

ε (v F (t )): ε (w) 

F

(− p (t ) I + 2 ν ε (v F (t )))n · w = 

f F (t ) · w

∀ w ∈ H 1 F ( F ) ,

(7)

F

with H  F ( F ) := {w ∈ [ H ( F )] ; w = 0 on  F }. On the other hand, the weak formulation of the incompressibility condition reads 1

1

d



∀ q ∈ L 2 ( F ) .

q div(v F (t )) = 0 F

Multiplying the third equation in (6) by a test function w :  S → Rd such that w = 0 on  S and integrating by parts, we obtain

ρS



d





v S (t ) · w +

dt S

σ S (t ) : ε (w) +



σ S (t )n · w = 

S

f S (t ) · w

∀ w ∈ H 1 S ( S ) ,

S

where H  S ( S ) := {w ∈ [ H ( S )] ; w = 0 on  S }. Then, multiplying the fourth equation in (6) by a test function Rd×d and integrating in  S , we get 1

d

1



dt

d





σ S (t ) : τ = λ S

div(v S (t ))tr τ + 2 μ S

d×d

(8)

ε (v S (t )) : τ

τ : S →

×d ∀ τ ∈[ L 2 ( S )]dsym ,

S d×d

with [ L ( S )]sym := {τ ∈ [ L ( S )] ; = τ in  S }. Taking into account the continuity of the velocities across  (fifth equation in (6)), we can look for a global velocity v(t ) ∈ [ H 01 ()]d defined by 2

2

v := v F in (0, T ) ×  F ,

τt

v := v S in (0, T ) ×  S ,

and consider the initial global velocity v0 defined by

v0 := v0F in  F ,

v := v0S in  S ,

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

which is continuous across  thanks to (4). Similarly, we use body force, independently of the domain:

ρ and f to denote, respectively, the density and the given

ρ := ρ F in  F , ρ := ρ S in  S , f := f F in (0, T ) ×  F ,

279

f := f S in (0, T ) ×  S .

(9)

Then, adding equations (7) and (8), and using the continuity of the normal stresses across  (sixth equation in (6)), we ×d × L 2 ( )) arrive at the following variational formulation of problem (6): find (v, σ S , p ) ∈ L 2 (0, T ; [ H 01 ()]d × [ L 2 ( S )]dsym F such that

⎧ 0 0 ⎪ ⎪ v(0) = v in  , σ S (0) = σ S in  S ⎪ ⎪     ⎪ ⎪ d  ⎪ ⎪ ⎪ ρ v(t )· w + 2ν ε (v(t )): ε (w) + σ S (t ): ε (w) − p (t ) div(w) = f(t )· w ⎪ ⎪ dt ⎪ ⎪ ⎨  S  F    F d

⎪ σ S (t ): τ = λ div(v(t ))tr(τ ) + 2μ ε (v(t )) : τ ⎪ ⎪ dt ⎪ ⎪ ⎪ S S ⎪  S ⎪ ⎪ ⎪ ⎪ q div(v(t )) = 0 ⎪ ⎪ ⎩

(10)

F

×d × L 2 ( ) and for a.e. t ∈ (0, T ). for any (w, τ , q) ∈ [ H 01 ()]d × [ L 2 ( S )]dsym F Now, in order to prove that problem (10) is well-posed, we assume that

f F ∈ H 1 (0, T ; [ L 2 ( F )]d ) ,

f S ∈ H 1 (0, T ; [ L 2 ( S )]d ) ,

1 v0F ∈ [ H  ( F ) ∩ H 2 ( F )]d F

with div(v0F ) = 0 in  F ,

(11)

1 u0S , v0S ∈ [ H  ( S ) ∩ H 2 ( S )]d , S

v0F = v0S

on  ,

and that there exists p 0 ∈ H 1 ( F ) such that



 − p 0 I + 2 ν ε (v0F ) n = σ 0S n .

(12)

1 Theorem 1. Let (v F , p , u S ) ∈ L 2 (0, T ; [ H  ( F )]d × L 2 ( F ) × [ H 1 S ( S )]d ) be the unique solution to problem F

⎧ v F (0) = v0F in  F , u S (0) = u0S in  S , ∂t u S (0) = v0S in  S , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ t ⎪ ⎪ ⎪ ⎪ ⎪ v F ( z) dz = u S (t ) − u0S on  ⎪ ⎪ ⎪ ⎪ 0 ⎪ ⎪    ⎪ ⎪ d d2 ⎪ ⎪ ⎪ ρF v F (t )· w + 2 ν ε (v F (t )): ε (w) + ρ S 2 u S (t )· w ⎪ ⎪ dt dt ⎨ F F S   ⎪ ⎪ +2 μ ε (u S (t )): ε (w) + λ div(u S (t )) div(w) ⎪ ⎪ ⎪ ⎪ ⎪ S S  ⎪   ⎪ ⎪ ⎪ ⎪ ⎪ − p ( t ) div ( w ) = f ( t )· w + f S (t )· w F ⎪ ⎪ ⎪ ⎪ ⎪ F F S ⎪  ⎪ ⎪ ⎪ ⎪ q div(v F (t )) = 0 ⎪ ⎪ ⎩

(13)

F

for any (w, q) ∈ [ H 01 ()]d × L 2 ( F ) and a.e. t ∈ (0, T ). Then, (v, σ S , p ), with

v := v F in (0, T ) ×  F ,

v := ∂t u S in (0, T ) ×  S ,

σ S = λ div(u S ) I + 2 μ ε (u S ) in (0, T ) ×  S , is a solution to problem (10). Conversely, if (v, σ S , p ) is a solution to problem (10), then (v F , p , u S ), with

(14)

280

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

t v F := v|(0, T )× F ,

v( z)| S dz + u0S ,

u S (t ) :=

(15)

0

is the unique solution to problem (13). Moreover,

σ S (t ) = λ div(u S (t )) I + 2 μ ε (u S (t ))

a.e. t ∈ (0, T ) .

In particular, problem (10) has a unique solution. Proof. First of all, we recall from Theorem 3.4 in [11] that, under hypotheses (11) and (12), problem (13) has a unique 1 solution (v F , p , u S ) ∈ L 2 (0, T ; [ H  ( F )]d × L 2 ( F ) × [ H 1 S ( S )]d ), with the regularity F 1 v F , ∂t v F ∈ L ∞ (0, T ; [ L 2 ( F )]d ) ∩ L 2 (0, T ; [ H  ( F )]d ) , F 1 u S , ∂t u S ∈ L ∞ (0, T ; [ H  ( S )]d ) , S

∂tt u S ∈ L ∞ (0, T ; [ L 2 ( S )]d ) ,

(16)

p ∈ L 2 (0, T ; L 2 ( F )) . 1 Let v and σ S be defined by (14). Notice that, by (16), v F ∈ L 2 (0, T ; [ H  ( F )]d ) and ∂t u S ∈ L ∞ (0, T ; [ H 1 S ( S )]d ). MoreF over, taking the time derivative of the identity on  in (13) we have ∂t u S (t ) = v F (t ) on  and a.e. t ∈ (0, T ), so that

(v(t )| F )| = v F (t )| = ∂t u S (t )| = (v(t )| S )|

a.e. t ∈ (0, T ) .

1 Therefore, v ∈ L 2 (0, T ; [ H 01 ()]d ). On the other hand, u S ∈ L ∞ (0, T ; [ H  ( S )]d ) guarantees that S

Now, let w ∈ nitions of v and that

d

[ H 01 ()]d and work for a.e. t ∈ (0, T ). Decomposing S provided in (14) and the global notation for

σ

the integral over  into  F and  S , using the defi-

ρ (see (9)), and that ρ F and ρ S are constant, we have







ρ v(t )· w + 2 ν

dt

×d ). σ S ∈ L 2 (0, T ; [ L 2 ( S )]dsym



F

= ρF

S



d

v F (t )· w + ρ S

dt



F

+λ  =



ε (v(t )): ε (w) + σ S (t ): ε (w) −



u S (t )· w + 2 ν

dt 2



S

div(u S (t )) div(w) + 2 μ S

f F (t )· w +

F



d2

S

p (t ) div(w) =

F

 f S (t )· w =

ε (v F (t )): ε (w) 

F

ε (u S (t )): ε (w) − S

p (t ) div(w)

F

f(t )· w , 

where in the penultimate identity we used that (v F , p , u S ) is the solution to problem (13) and in the last step we used the definition of f (see (9)). Thus, the first equation of problem (10) is satisfied. ×d . Then, using the definitions of σ and v in  , we have Next, let τ ∈ [ L 2 ( S )]dsym S S

d



σ S (t ): τ = λ

dt

d

div(u S (t ))tr(τ ) + 2 μ

dt



S



=λ

S



div(v(t ))tr(τ ) + 2 μ S

d



ε (u S (t )) : τ

dt S

ε (v(t )) : τ . S

Concerning the last equation of problem (10), notice that it is exactly the same as the last one of problem (13). Finally, it is also straightforward to verify that v(0) = v0 in  and σ S (0) = σ 0S in  S . Therefore, (v, σ S , p ) is a solution to problem (10) and, in particular, problem (10) has at least one solution. Conversely, let (v, σ S , p ) be a solution to problem (10) and define v F and u S (t ) by (15). Since v ∈ L 2 (0, T ; [ H 01 ()]d ), 1 we have that v F ∈ L 2 (0, T ; [ H  ( F )]d ) and v S := v|(0, T )× S ∈ L 2 (0, T ; [ H 1 S ( S )]d ). In particular, taking into account that F 1 u0S ∈ [ H  ( S )]d , it follows from (15) that u S ∈ H 1 (0, T ; [ H 1 S ( S )]d ). S From the penultimate equation of problem (10), we deduce that

∂t σ S (t ) = λ div(v(t )| S ) I + 2 μ ε (v(t )| S ) Then, we integrate in [0, t ] for t ∈ (0, T ] and use that

a.e. t ∈ (0, T ) .

σ S (0) = σ 0S in  S to obtain that

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

t

σ S (t ) = σ 0S + λ div(

t v( z)| S dz) I + 2 μ ε (

0

v( z)| S dz) , 0

which is equivalent to the following, thanks to the definitions of

σ S (t ) = λ div(u S (t )) I + 2 μ ε (u S (t ))

σ 0S and u S (see (15)):

a.e. t ∈ (0, T ) .

(17)

Now, let w ∈ [ H 01 ()]d and t ∈ (0, T ). Then, using the definition of

ρF

d



 v F (t )· w + 2 ν

dt



F

ε (v F (t )): ε (w) + ρ S 

F

+2 μ

ε (u S (t )): ε (w) + λ S

=

d



 =







u S (t )· w

dt 2



S



S

f(t )· w =

d2

ρ and (17), we have



div(u S (t )): div(w) −

ρ v(t )· w + 2 ν

dt

281

p (t ) div(w)

 F



ε (v(t )): ε (w) + σ S (t ): ε (w) − 

F

f F (t )· w +

F

S

p (t ) div(w) =

F

a.e. t ∈ (0, T ) ,

f S (t )· w

S

where in the penultimate step we used that (v, σ S , p ) is a solution to problem (10) and in the last step we just employed the definition of f. Thus, the first variational equation of problem (13) is satisfied. Here again, we notice that the last equation of problem of (13) is exactly the same as that of problem (10) thanks to (15). Finally, from (15) and (10), it is straightforward that

v F (0) = v(0)| F = v0 | F = v0F u S (0) =

u0S

,

in  F

∂t u S (0) = v(0)| S = v | S = v0S 0

in  S

and for each t ∈ (0, T ) we know by the transmission condition that

t

t v F ( z)| dz =

0

t (v( z)| F )| dz =

0

t ∂t u S ( z)| dz = u S (t )| − u0S | .

(v( z)| S )| dz = 0

0

Therefore, (v F , p , u S ) is the solution to problem (13).

2

Remark. Assumption (12) is used in the previous Theorem to prove the existence of a unique solution to the variational formulation (10). It will also be used to obtain a priori estimates on the solution to the semidiscrete problem; more precisely, it is used in Lemma 2 to obtain a priori estimates on the derivative of the velocity at time t = 0. 4. Semidiscretization in space In what follows, we assume (11) and (12). From now on, we also suppose that  F and  S are two-dimensional polygons or three-dimensional polyhedra. We let h be a discretization parameter associated with a mesh Th of . We assume that the mesh Th consists of triangles in two-dimensions or tetrahedra in three dimensions, and that it is compatible with the interface  (i.e. the elements of Th do not cross the interface ). We also assume that there exists a mesh Th0 such that ×d and Q ⊂ L 2 ( ) be finite ¯ ) ∩ H 1 ()]d , S h ⊂ [ L 2 ( S )]dsym for all h < h0 , Th is a refinement of Th0 . For each h, let V h ⊂ [C ( F h 0 dimensional subspaces of piecewise polynomials over the mesh Th such that

ε ( V h | S ) ⊆ S h .

(18)

[ H 01 ( F )]d ,

Further, we suppose that V h contains piecewise linear functions and that the pair ( V h | F ∩ stable for the Stokes equations, that is to say, there exists some β F > 0, independent of h, such that

Qh ∩

L 20 ( F ))

is



qh div(wh ) sup wh ∈ V h | F ∩[ H 01 ( F )]d , wh =0

F

wh 1, F

≥ β F qh 0, F

∀ qh ∈ Q h ∩ L 20 ( F ) .

(19)

We remark that, under the previous hypotheses, Theorem 3.5 in [12] states that the stability assumption (19) guarantees that the pair ( V h , Q h ) satisfies the following inf–sup condition:

282

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

 qh div(wh ) F

sup

≥ β˜ F qh 0, F

wh 1,

wh ∈ V h \{0}

∀qh ∈ Q h ,

(20)

for some β˜ F > 0, independent of h. We also assume that the spaces V h , S h and Q h are suitable to approximate the corresponding continuous ones. In other words, that they satisfy the following approximation properties for any r ∈ [0, k], k ∈ N:

∀w ∈ [ H 01 () ∩ H 1+r ()]d

∃{wh }h with wh ∈ V h s.t. wh − w1, ≤ C hr w1+r , ,

×d ∩ [ H r ( )]d×d ∀τ ∈ [ L 2 ( S )]dsym S

∀q ∈ L ( F ) ∩ H ( F ) 2

r

∃{τ h }h with τ h ∈ S h s.t. τ h − τ 0, S ≤ C hr τ r , S ,

(21)

∃{qh }h with qh ∈ Q h s.t. qh − q0, F ≤ C h qr , F . r

Finally, we introduce a suitable approximation of the initial conditions, assuming that the initial data have the regularity v0 ∈ [ H 1+r ()]d , σ 0S ∈ [ H r ( S )]d×d and p 0 ∈ H r ( F ). To this end, we first consider (vh0, F , p h0 ) ∈ V h | F × Q h the finite element solution of the following Stokes system:

   ⎧  0 0 0 ⎪ 2 ν ε ( v ) : ε ( w ) − p div ( w ) = 2 ν ε ( v ) : ε ( w ) − p 0 div(wh ) ⎪ h h h h, F h ⎪ ⎪ ⎪ ⎪ F F F ⎨  F 0 q div ( v ) = 0 h ⎪ h, F ⎪ ⎪ ⎪ ⎪ F ⎪ ⎩ 0 (vh, F , wh )0, = (v0 , wh )0,

(22)

for all (wh , qh ) ∈ V h | F × Q h . Notice that the discrete inf–sup condition (20) guarantees that



qh div(wh ) sup

F

wh ∈ V h | F \{0}

wh 1, F

≥ β˜ F qh 0, F

∀qh ∈ Q h ,

which is the discrete counterpart of the following result proved in Theorem 3.3 in [11]:



q div(w) F

sup w∈ H 01 ()| F

\{0}

w1, F

≥ α˜ F q0, F

∀q ∈ L 2 ( F ) ,

˜ F > 0. Consequently, the discrete problem (22) is well-posed and, under our regularity assumptions on the initial for some α data, it holds

vh0, F − v0 1, F +  ph0 − p 0 0, F ≤ C hr (v0 1+r , F +  p 0 r , F ) , (see Theorem 1.1 of Chapter II in [16]). Besides, we define vh0, S ∈ V h | S as the solution of the following problem:

⎧ ⎨ (∇ vh0, S , ∇ wh )0, S = (∇ v0 , ∇ wh )0, S ⎩

(vh0, S , wh )0, = (v0 , wh )0,

(23)

for all wh ∈ V h | S . Then, under our regularity assumption on the initial datum v0 , we have that

vh0, S − v0 1, S ≤ C hr v0 1+r , S . Further, notice that vh0 defined piecewise by

vh0 = vh0, F in  F , belongs to



V˜ h := wh ∈ V h ;

vh0 = vh0, S in  S ,

 qh div(wh ) = 0 F

and also satisfies

∀qh ∈ Q h ,

(24)

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

vh0 − v0 1, ≤ C hr (v0 1+r , F + v0 1+r , S ) . On the other hand, let





σ h0 : τ h = S

283

(25)

σ ∈ S h be the orthogonal projection of σ with respect to the [ L 2 ( S )]d×d -inner product, that is, 0 h

0 S

σ 0S : τ h ∀ τ h ∈ S h .

(26)

S

Then, under our regularity assumption on the initial datum

σ 0S ,

σ h0 − σ 0S 0, S ≤ C hr σ 0S r , S .

(27)

In particular, we have that

vh0 1, +  ph0 0, F ≤ C 0 (v0 1+r , F + v0 1+r , S +  p 0 r , F ) ,

σ h0 0, S ≤ C 0 σ 0S r , S .

(28)

Making use of the above approximation of the initial conditions, we propose the following semidiscrete problem: find

((vh , σ h ), ph ) ∈ C 1 ([0, T ]; V h × S h ) × C 0 ([0, T ]; Q h ) such that vh (0) = vh0 in , σ h (0) = σ h0 in  S and

⎧    d ⎪ ⎪ ρ v ( t )· w + 2 ν ε ( v ( t )): ε ( w ) + σ h (t ): ε (wh ) ⎪ h h h h ⎪ dt ⎪ ⎪ ⎪   F  S ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ − ph (t ) div(wh ) = f(t )· wh ⎪ ⎪ ⎨   F   d ⎪ ⎪ σ ( t ) : τ = λ div ( v ( t )) tr ( τ ) + 2 μ ε (vh (t )) : τ h ⎪ h h h h ⎪ ⎪ dt ⎪ ⎪ S S ⎪  S ⎪ ⎪ ⎪ ⎪ ⎪ qh div(vh (t )) = 0 ⎪ ⎪ ⎩

(29)

F

for any (wh , τ h , qh ) ∈ V h × S h × Q h and t ∈ (0, T ]. In the next Theorem we show that problem (29) is well-posed. Theorem 2. Problem (29) has a unique solution ((vh , σ h ), p h ) ∈ C 1 ([0, T ]; V h × S h ) × C 0 ([0, T ]; Q h ). Proof. Let us consider the following auxiliary problem: find (vh , σ h ) ∈ C 1 ([0, T ]; V˜ h × S h ) such that for all t ∈ (0, T ],

   ⎧ d  ⎪ ρ v ( t ) · w + 2 ν ε ( v ( t )) : ε ( w ) + σ ( t ) : ε ( w ) = f(t ) · wh ∀wh ∈ V˜ h ⎪ h h h h h h ⎪ ⎪ dt ⎪ ⎪ ⎨     S  F d σ h (t ) : τ h = λ div(vh (t ))tr τ h + 2 μ ε (vh (t )) : τ h ∀ τ h ∈ S h ⎪ ⎪ dt ⎪ ⎪ ⎪ S S ⎪ S ⎩ vh (0) = vh0 in  , σ h (0) = σ h0 in  S ,

(30)

where V˜ h is defined in (24). Let {wi }iN=1 and {τ k }kM=1 be bases of V˜ h and S h , respectively. Then, since vh0 ∈ V˜ h , we can write

vh (t ) =

N  j =1

v j (t )w j ,

σ h (t ) =

M  l =1

σl (t )τ l , vh0 =

N 

v 0j w j ,

j =1

σ h0 =

M 

σl0 τ l ,

l =1

and problem (30) is equivalent to the following initial value problem:

⎧     N N M    ⎪ ⎪  ⎪ v j (t ) ρ w j · wi = f(t ) · wi − v j (t ) 2ν ε (w j ) : ε (wi ) − σl (t ) τ l : ε (wi ) ⎪ ⎪ ⎪ ⎪ ⎪ j =1 l =1 ⎨ j =1   F S   M N 

  ⎪ ⎪ σl (t ) τ l : τ k = v j (t ) λ div(w j )tr τ k + 2μ ε (w j ) : τ k ⎪ ⎪ ⎪ ⎪ j = 1 l = 1 ⎪ S S S ⎪ ⎩ v i (0) = v 0i , σk (0) = σk0

for all i = 1, . . . , N and k = 1, . . . , M. By virtue of Picard–Lindelöf Theorem, the previous system has a unique solution in C 1 ([0, T ]), and thus, problem (30) has a unique solution (vh , σ h ) ∈ C 1 ([0, T ]; V˜ h × S h ).

284

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

Now, we proceed as in [12] and consider the following problem: find p h ∈ C ([0, T ]; Q h ) such that





d

p h (t ) div(wh ) =



ρ vh (t )· wh + 2 ν

dt

F







ε(vh (t )): ε (wh ) + σ h (t ): ε (wh ) − F

S

f(t )· wh 

for all wh ∈ V h and t ∈ [0, T ], where (vh , σ h ) ∈ C 1 ([0, T ]; V˜ h × S h ) is the unique solution to problem (30). Using Lemma 4.1 of Chapter I in [16] and the inf–sup condition (20), we can ensure the existence and uniqueness of a p h ∈ C ([0, T ]; Q h ) satisfying the first equation of (29). On the other hand, we recall that since vh ∈ C 1 ([0, T ]; V˜ h ), the third equation in (29) is automatically satisfied. Therefore, problem (29) has a unique solution. 2 We dedicate the rest of this section to the analysis of the semidiscrete problem (29). More precisely, in the next subsection we obtain some a priori estimates on the semidiscrete solution. After that, we take advantage of such estimates to prove convergence of the sequence of semidiscrete solutions to the unique solution of the weak formulation (10) as h → 0+ . Finally, we provide error estimates on the semidiscrete solutions. 4.1. A priori estimates In this subsection we derive some a priori estimates on the semidiscrete solution and its derivatives. A key ingredient is the following version of Gronwall’s Lemma (see Appendix B in [13]): Let g ∈ L 1 (0, T ) such that

t 0 ≤ g (t ) ≤ a1

g (s) ds + a2

a.e. t ∈ [0, T ]

0

for some a1 , a2 ≥ 0. Then

g (t ) ≤ a2 (1 + a1 t ea1 t )

a.e. t ∈ [0, T ] .

(31)

Lemma 1. Let us denote ρmin := min(ρ F , ρ S ), ρmax := max(ρ F , ρ S ),

C ∗ :=

T

ρmin

+1+(

T

ρmin

− 1) e T /ρmin ,

C 1 := f2L 2 (0, T ;[ L 2 ()]d ) +



C 02 4μ

σ 0S r2, S + C 02 ρmax (v0 1+r , F + v0 1+r , S +  p 0 r , F )2 ,

C 2 := C 1 1 + C ∗ . There holds

vh (t )20, ≤

C1

ρmin

1+

t

ρmin

et /ρmin

in (0, T ) ;

(32)

in particular,

vh 2L ∞ (0,T ;[ L 2 ()]d ) ≤

C1

ρmin

1+

T

ρmin

e T /ρmin ,

vh 2L 2 (0,T ;[ L 2 ()]d ) ≤ C 1 C ∗ .

(33)

Besides,

ε (vh )2L 2 (0,T ;[ L 2 (

d×d ) F )]

≤

C2 4ν

t C2 div( vh ( z) dz)20, S ≤ , λ 0

(34)

, t C2 ε ( vh ( z) dz)20, S ≤

μ

in (0, T ) .

(35)

0

Proof. We proceed in the usual way and, for t ∈ (0, T ), choose wh = vh (t ) ∈ V h and qh = p h (t ) ∈ Q h in the first and third equations of (29), respectively. Then, using that ∂t vh (t ) · vh (t ) = 12 ∂t |vh (t )|2 , we have

1 d





ρ |vh (t )|2 + 2 ν ε (vh (t ))20, F +

2 dt 



σ h (t ) : ε (vh (t )) = S

f(t ) · vh (t ) . 

(36)

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

285

To deal with the third term on the left hand side of (36), we integrate with respect to t the second equation of the semidiscrete problem (29) and then take τ h = ε (vh (t )) ∈ S h :



 0 h:

σ h (t ): ε (vh (t )) = S

σ ε (vh (t )) +



λd 2 dt

S

t  t d 2 |div( vh ( z) dz)| + μ |ε ( vh ( z) dz)|2 . dt

S

0

S

0

Substituting the above identity in (36), integrating in [0, t ], and then using the initial condition vh (0) = vh0 as well as the Cauchy–Schwarz inequality and Young’s inequality, we obtain

t

t ε (vh (s))20, F ds + λ div( vh ( z) dz)20, S

ρmin vh (t )20, + 4ν 0

0

(37)

t t +μ ε ( vh ( z) dz)20, S ≤ C 1 + vh ( z)20, dz . 0

0

Notice that, in particular,

t 2 min vh (t )0,

ρ

vh ( z)20, dz ;

≤ C1 + 0

then, applying Gronwall’s Lemma (31), we deduce (32), from which (33) follows straightforwardly. Finally, using (32) in inequality (37), we obtain (34) and (35). 2 In the following two lemmas we use similar techniques to deduce bounds for the time derivative of vh , first at t = 0 and later in (0, T ]. Lemma 2. There holds



ρ |∂t vh (0)|2 ≤

4

C 0 ,

ρmin



where

C 0 := 4 ν 2 div ε (v0 )20, F + div σ 0S 20, S +  p 0 21, F + f(0)20, . Proof. Taking the first equation in (29) at t = 0 and choosing wh = ∂t vh (0) ∈ V h , we have







ρ |∂t vh (0)|2 + 2 ν 



−

ε (vh0 ): ε (∂t vh (0)) + σ h0 : ε (∂t vh (0)) 

F

p h (0) div(∂t vh (0)) =

F

S

f(0)·∂t vh (0) . 

Taking the time derivative of the third equation in (29), we have



div(∂t vh (t )) qh = 0 ∀qh ∈ Q h .

(38)

F

In particular,



 p h0 div(∂t vh (0)) = 0 .

p h (0) div(∂t vh (0)) = F

F

Therefore, using the definition of the approximated initial conditions (22) and (26), and integrating by parts,





ρ |∂t vh (0)|2 = 2 ν 



+



div ε (v0 )·∂t vh (0) + F



∇ p 0 ·∂t vh (0)

F

(−2 ν ε (v0 )n + σ 0S n + p 0 n )·∂t vh (0) .

f(0)·∂t vh (0) + 

S



div σ 0S ·∂t vh (0) −



286

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

We remark that the last term vanishes thanks to the transmission condition assumed to be satisfied by the initial data (see (12)). Consequently, applying the Cauchy–Schwarz inequality, we conclude that

√





ρmin ( ρ |∂t vh (0)|2 )1/2 ≤ 2 ν div ε (v0 )0, F + div σ 0S 0, S +  p 0 1, F + f(0)0, ≤ 2 C 0 . 2 

  4 + (λ + 2 μ) C 0 (v0 1+r , F + v0 1+r , S +  p 0 r , F ) + ρmin C 0 and C 2 := C 1 1 + C ∗ .

Lemma 3. Let C 1 := ∂t f22

L (0, T ;[ L 2 ()]d )

There holds

∂t vh (t )20, ≤

C 1

ρmin

t

1+

ρmin

et /ρmin

in [0, T ] .

(39)

In particular,

∂t vh 2L 2 (0,T ;[ L 2 ()]d ) ≤ C 1 C ∗ ,

(40)

and

 ε ∂t vh )2L 2 (0,T ;[ L 2 ( div(vh (t ))20, S ≤

≤

d×d ) F )]

C 2

C 2 4ν

C 2

ε (vh (t ))20, S ≤

,

λ

(41)

, 2μ

(42)

.

Proof. Taking the time derivative of the first equation in (29), choosing wh = ∂t vh (t ) ∈ V h and making use of (38) for qh = ∂t p h (t ) ∈ Q h , we have



1





ρ ∂t (|∂t vh (t )|2 ) + 2 ν ε (∂t vh (t ))20, F + ∂t σ h (t ): ε (∂t vh (t )) =

2 

S 1 2

∂t f(t )·∂t vh (t ) ,

(43)



where we used that ∂t (∂t vh (t )) · ∂t vh (t ) = ∂t (|∂t vh (t )| ). To deal with the third term in (43), we take τ h = ε (∂t vh (t )) ∈ S h in the second equation of (29):



λ d

∂t σ h (t ): ε (∂t vh (t )) =

2 dt

2

div(vh (t ))20, S + μ

d dt

ε (vh (t ))20, S .

S

Substituting the previous identity in (43), we obtain



1 d

ρ |∂t vh (t )|2 + 2 ν ε (∂t vh (t ))20, F +

2 dt 

+μ

d dt

 ε (vh (t ))20, S =

λ d 2 dt

div(vh (t ))20, S (44)

∂t f(t )·∂t vh (t ) . 

Integrating (44) in [0, t ], using the Cauchy–Schwarz inequality, Young’s inequality, Lemma 2 and (28), we have

1



t 2

ε (∂t vh ( z))20, F dz +

ρ |∂t vh (t )| + 2 ν

2 

0

t  =

∂t f( z)·∂t vh ( z) dz + 0 

≤

C 1 2

+

1

t

1

λ 2

div(vh (t ))20, S + μ ε (vh (t ))20, S



ρ |∂t vh (0)|2 +

2

λ 2

div(vh0 )20, S + μ ε (vh0 )20, S

(45)



∂t vh ( z)20, dz .

2 0

In particular, 2 min ∂t vh (t )0,

ρ



t ∂t vh ( z)20, dz .

≤ C1 + 0

Then, applying Gronwall’s Lemma (31), we deduce (39); inequality (40) can be obtained by integrating in [0, T ]. Inequalities (41) and (42) follow then straightforwardly from (45). 2

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

287

Next we obtain an a priori bound for the semidiscrete approximation of the Cauchy stress tensor and its time derivative. Lemma 4. There holds

 √ √ σ h (t )0, S ≤ C 0 σ 0S r , S + C 2 ( λ d + 2 μ)  √  ∂t σ h (t )0, S ≤ C 2 ( λ d + 2 μ) in (0, T ) .

In particular,

σ h 2L 2 (0,T ;[ L 2 (

in (0, T ) ,

(46) (47)

√   √ ≤ 2 T C 02 σ 0S r2, S + C 2 ( λ d + 2 μ)2 , √  ≤ T C 2 ( λ d + 2 μ)2 . )]d×d )

(48)

d×d ) S )]

∂t σ h 2L 2 (0,T ;[ L 2 (

S

Proof. We integrate the second equation in (29) in [0, t ] and then take

 σ h (t )20, S =

t



(49)

τ h = σ h (t ) ∈ S h :



t

σ h0 : σ h (t ) + λ div( vh (z) dz) tr σ h (t ) + 2μ ε ( vh (z) dz): σ h (t ) . S

0

S

S

0

Then, inequality (46) follows using the Cauchy–Schwarz inequality, that tr τ 0, S ≤ On the other hand, the second equation in (29) for τ h = ∂t σ h (t ) ∈ S h states that

  ∂t σ h (t )20, S = λ div(vh (t )) tr ∂t σ h (t ) + 2μ ε (vh (t )):∂t σ h (t ) S

√

dτ 0, S , (28) and (35).

in (0, T ) .

S

√

As above, we deduce (47) just using the Cauchy–Schwarz inequality, that tr τ 0, S ≤ dτ 0, S and (42). Inequalities (48) and (49) follow straightforwardly integrating in [0, T ] the squares of (46) and (47), respectively. 2 Finally, we obtain an a priori bound for the pressure. Lemma 5. It holds

 ph 2L 2 (0,T ; L 2 (

F

≤ ))

4

β˜ F2

2 C ∗ C 1 + ν C 2 f2L 2 (0,T ;[ L 2 ()]d ) + ρmax

(50)

√ √ + 2 T (C 02 σ 0S r2, S + C 2 ( λ d + 2 μ)2 ) .

Proof. For t ∈ (0, T ), we take qh := p h (t ) ∈ Q h in (20). Then,



p h (t ) div(wh )

 ph (t )0, F ≤

1

β˜ F

F

sup

wh 1,

wh ∈ V h \{0}

Now, from the first equation in (29),





p h (t ) div(wh ) = − F





f(t )· wh + 

. 

ρ ∂t vh (t )· wh + 2ν ε (vh (t )): ε (wh ) + σ h (t ): ε (wh ) , 

F

S

for all wh ∈ V h and t ∈ (0, T ). Thus, using the Cauchy–Schwarz inequality, we have

 ph (t )0, F ≤

1

β˜ F

f(t )0, + ρmax ∂t vh (t )0, + 2ν ε (vh (t ))0, F + σ h (t )0, S

in (0, T ) .

In particular, taking squares and integrating in [0, T ], we obtain

 ph 2L 2 (0,T ; L 2 (

F ))

≤

4

β˜ F2

2 ∂t vh 2L 2 (0,T ;[ L 2 ()]d ) f2L 2 (0,T ;[ L 2 ()]d ) + ρmax

+ 4 ν 2 ε (vh )2L 2 (0,T ;[ L 2 ( so that inequality (50) follows using (40), (34) and (48).

2

F

]d×d ))

+ σ h 2L 2 (0,T ;[ L 2 (

S

)]d×d )

,

288

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

4.2. Convergence In this subsection we prove the following convergence result. Theorem 3. The sequence of solutions {(vh , σ h , p h )}h>0 to problems (29) converges to the unique solution to problem (10), (v, σ S , p ), in H 1 (0, T ; [ L 2 ()]d ) × H 1 (0, T ; [ L 2sym ( S )]d×d ) × L 2 (0, T ; L 2 ( F )) as h → 0+ . Proof. Using the a priori estimates on the semidiscrete solution, convergence of the sequence of solutions of the semidiscrete problems to the exact solution of (10) follows from standard arguments. Indeed, in the previous subsection we have shown ×d ) × L 2 (0, T ; L 2 ( )). Moreover, by that {(vh , σ h , p h )}h is uniformly bounded in H 1 (0, T ; [ L 2 ()]d ) × H 1 (0, T ; [ L 2 ( S )]dsym F 1 (34), {vh |(0, T )× F }h is uniformly bounded in L ∞ (0, T ; [ H  ( F )]d ) and by (42), {vh |(0, T )× S }h is uniformly bounded in F 1 L ∞ (0, T ; [ H  ( S )]d ). Therefore, there exists a subsequence of {(vh , σ h , ph )}h , that we will denote identically, such that S

vh v in H 1 (0, T ; [ L 2 ()]d ) , 1

(51)

d×d

2

σ h σ in H (0, T ; [ L ( S )]sym ) , ph p

(52)

in L 2 (0, T ; L 2 ( F )) ,

vh |(0, T )× F v|(0, T )× F vh |(0, T )× S v|(0, T )× S

(53)

1 in L 2 (0, T ; [ H  ( F )]d ) , F 2

1

(54)

d

in L (0, T ; [ H  S ( S )] ) ,

(55)

1 ×d ) × L 2 (0, T ; L 2 ( )) with v| 2 d for some (v, σ , p ) ∈ H 1 (0, T ; [ L 2 ()]d ) × H 1 (0, T ; [ L 2 ( S )]dsym F (0, T )× F ∈ L (0, T ; [ H  F ( F )] ) and 1 v|(0, T )× S ∈ L 2 (0, T ; [ H  ( S )]d ). S

×d ) × L 2 (0, T ; L 2 ( )) solves problem (10), we In order to show that (v, σ , p ) ∈ H 1 (0, T ; [ L 2 ()]d ) × H 1 (0, T ; [ L 2 ( S )]dsym F

×d × L 2 ( )). start studying the variational equations. To this end, let us consider (w, τ , q) ∈ C 0 ([0, T ]; [ H 01 ()]d × [ L 2 ( S )]dsym F Notice that, thanks to our assumption (21) about the density of the approximation spaces, we can take (wh , τ h , qh ) ∈ C 0 ([0, T ]; V h × S h × Q h ) such that (wh , τ h , qh ) → (w, τ , q) in C 0 ([0, T ]; V h × S h × Q h ). Using (29) and integrating in [0, T ], we have

⎧ T T T ⎪ ⎪ ⎪ ⎪ ⎪ ρ ∂t vh (t )· wh (t ) dt + 2ν ε (vh (t )): ε (wh (t )) dt + σ h (t ): ε (wh (t )) dt ⎪ ⎪ ⎪ ⎪ ⎪ 0  0  0  F S ⎪ ⎪ ⎪ T T ⎪ ⎪ ⎪ ⎪ ⎪ − p h (t ) div(wh (t )) dt = f(t )· wh (t ) dt ⎪ ⎪ ⎪ ⎪ ⎨ 0 0  F

T T T ⎪ ⎪ ⎪ ⎪ ⎪ ∂ σ ( t ): τ ( t ) dt = λ div ( v ( t )) tr ( τ ( t )) dt + 2 μ ε (vh (t )): τ h (t ) dt t h h h h ⎪ ⎪ ⎪ ⎪ ⎪ 0 S 0 S 0 S ⎪ ⎪ ⎪ ⎪ ⎪ T ⎪ ⎪  ⎪ ⎪ ⎪ qh (t ) div(vh (t )) dt = 0 . ⎪ ⎪ ⎩

(56)

0 F

We may pass to the limit as h → 0+ in the first equation of (56) making use of the weak convergences (51), (54), (52) and (53), so that

T

T

ρ ∂t v(t )· w(t ) dt + 2ν 0 

T −

T

ε(v(t )): ε (w(t )) dt + 0 F

0 F

0 S

T

p (t ) div(w(t )) dt =

σ (t ): ε (w(t )) dt

f(t )· w(t ) dt , 0 

for any w ∈ C 0 ([0, T ]; H 01 ()), and by density, for any w ∈ L 2 (0, T ; H 01 ()). Therefore,

(57)

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299









ρ ∂t v(t )· w + 2ν ε (v(t )): ε (w) + σ (t ): ε (w) − 

F

S

289

 p (t ) div(w) =

F

f(t )· w , 

for all w ∈ and a.e. t ∈ (0, T ). Similarly, we may pass to the limit as h → 0+ in the second equation of (56) applying the weak convergences (52) and (55), so that H 01 ()

T

T ∂t σ (t ): τ (t ) dt = λ

0 S

T div(v(t )) tr(τ (t )) dt + 2μ

0 S

ε (v(t )): τ (t ) dt ,

(58)

0 S

×d ). Therefore, τ ∈ L 2 (0, T ; [ L 2 ( S )]dsym    ∂t σ (t ): τ = λ div(v(t )) tr(τ ) + 2μ ε (v(t )): τ ,

for any

S

S

S

d×d

for all τ ∈ [ L ( S )]sym and a.e. t ∈ (0, T ). Finally, we may pass to the limit as h → 0+ in the third equation of (56) taking advantage of (54) to obtain 2

T  ∀q ∈ L 2 (0, T ; L 2 ( F )).

q(t ) div(v(t )) dt = 0 , 0 F

Therefore,



q div(v(t )) dt = 0 , F

for all q ∈ L 2 ( F ) and a.e. t ∈ (0, T ). ×d ) × L 2 (0, T ; L 2 ( )), with Summing up, we have proved that (v, σ , p ) ∈ H 1 (0, T ; [ L 2 ()]d ) × H 1 (0, T ; [ L 2 ( S )]dsym F 1 v|(0, T )× F ∈ L 2 (0, T ; [ H  ( F )]d ) and v|(0, T )× S ∈ L 2 (0, T ; [ H 1 S ( S )]d ), satisfies the equations in (10). F

1 Notice that v ∈ H 1 (0, T ; [ L 2 ()]d ), with v|(0, T )× F ∈ L 2 (0, T ; [ H  ( F )]d ) and v|(0, T )× S ∈ L 2 (0, T ; [ H 1 S ( S )]d ), so that F

it remains to clarify that v ∈ H 1 (0, T ; [ H 01 ()]d ). To this end, we simply notice that, since vh (t ) ∈ V h , (vh |(0, T )× F )|(0, T )× = (vh |(0, T )× S )|(0, T )× for all h > 0. Then, passing to the limit as h → 0+ and using (54) and (55), we deduce that

(vh |(0,T )× F )|(0,T )× (v|(0,T )× F )|(0,T )× in L 2 (0, T ; [ H 1/2 ()]d ) , (vh |(0,T )× S )|(0,T )× (v|(0,T )× S )|(0,T )× in L 2 (0, T ; [ H 1/2 ()]d ) . Therefore, we have the following transmission condition across :

(v|(0,T )× F )|(0,T )× = (v|(0,T )× S )|(0,T )× . Let us prove now that v and

σ satisfy the initial conditions of problem (10). To this end, we first take w ∈

C 1 ([0, T ]; [ H 01 ()]d ) such that w( T ) = 0 in (57) and integrate by parts in the first term on the left hand side:

T −

T



ρ v(t )·∂t w(t ) dt − ρ v(0)· w(0) + 2ν 0  T

+



0 S

0 F

T

σ (t ): ε (w(t )) dt −

ε (v(t )): ε (w(t )) dt

p (t ) div(w(t )) dt = 0 F

(59)

T  f(t )· w(t ) dt . 0 

Besides, if we take wh ∈ C ([0, T ]; V h ) with wh ( T ) = 0 and wh → w in C 0 ([0, T ]; [ H 01 ()]d ) in the first equation of (56), integrating by parts in the first term on the left hand side and using that vh (0) = vh0 in , we have 1

T −

T



ρ vh (t )·∂t wh (t ) dt − ρ 0  T

+



vh0 · wh (0) + 2

T

σ h (t ): ε (wh (t )) dt − 0 S

ε (vh (t )): ε (wh (t )) dt

ν 0 F

T

p h (t ) div(wh (t )) dt = 0 F

f(t )· wh (t ) dt . 0 

290

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

If we pass to the limit as h → 0+ making use of the weak convergences (51), (54), (52) and (53) as well as the approximation properties of the initial conditions (25), we obtain

T

T

 0

−

ρ v(t )·∂t w(t ) dt − ρ v · w(0) + 2ν 0 



T

+

ε (v(t )): ε (w(t )) dt 0 F

T

σ (t ): ε (w(t )) dt − 0 S

p (t ) div(w(t )) dt = 0 F

(60)

T f(t )· w(t ) dt . 0 

From (59) and (60), we deduce that





ρ v(0) · w(0) = ρ v0 · w(0) 

∀w ∈ C 1 ([0, T ]; [ H 01 ()]d ) with w( T ) = 0 ,



which implies that v(0) = v0 in . ×d ) such that Similarly, let us consider τ ∈ C 1 ([0, T ]; [ L 2 ( S )]dsym

τ ( T ) = 0, as well as an approximating sequence τ h ∈ ×d ). We use the second equation in (56) and integrate C 1 ([0, T ]; S h ) such that τ h ( T ) = 0 and τ h → τ in C 0 ([0, T ]; [ L 2 ( S )]dsym by parts its left hand side, so that

T

T

 0 h:

−

σ h (t ):∂t τ h (t ) dt − σ τ h (0) = λ 0 S

T div(vh (t )) tr(τ h (t )) dt + 2 μ

0 S

S

ε (vh (t )): τ h (t ) dt . 0 S

Then, we may pass to the limit as h → 0+ using the weak convergences (52) and (55) as well as the approximation properties of the initial conditions (27), to obtain

T

T

 0 S:

−

σ (t ):∂t τ (t ) dt − σ τ (0) = λ 0 S

T

+2 μ

div(v(t )) tr(τ (t )) dt 0 S

S

(61)

ε(v(t )): τ (t ) dt . 0 S

On the other hand, integrating by parts in the left hand side of (58), we obtain

T −

T



σ (t ):∂t τ (t ) dt − 0 S

S

T

+2 μ

σ (0): τ (0) = λ

div(v(t )) tr(τ (t )) dt 0 S

(62)

ε(v(t )): τ (t ) dt . 0 S

From (61) and (62) we deduce that





σ 0S : τ (0) = S

×d σ (0): τ (0) , ∀ τ ∈ C 1 ([0, T ]; [ L 2 ( S ]dsym ) with τ ( T ) = 0 .

S

Therefore,

σ also satisfies the initial condition of problem (10) in  S . 2

4.3. Error estimates In this subsection we obtain error estimates on the global velocity and the Cauchy stress tensor under the same regularity assumptions on the domains  F and  S as those made in [12]. More precisely, we assume that the problem

⎧ find (¯v, p¯ ) ∈ [ H 01 ( F )]d × L 20 ( F ) such that ⎪ ⎪ ⎪ ⎪ ⎪    ⎪ ⎪ ⎪ ⎨ ρ ∇ v¯ · ∇ z − p¯ div(z) = f¯ F · z ∀ z ∈ [ H 01 ( F )]d F F F ⎪ ⎪ ⎪ ⎪ 2 ⎪ ⎪ q div(v) = 0 ∀ q ∈ L ( F ) ⎪ ⎪ ⎩ F

(63)

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

291

is H 2− 1 -regular for some 1 ∈ (0, 1), that is, for any f¯ F ∈ [ L 2 ( F )]d , the solution to problem (63) shows the regularity (¯v, p¯ ) ∈ [ H 2− 1 ( F )]d × H 1− 1 ( F ), − p¯ n + 2ε (¯v)n ∈ [ H 1/2− 1 ()]d , and it satisfies

¯v2− 1 , F +  p¯ 1− 1 , F +  − p¯ n + 2ε (¯v)n 1/2− 1 , ≤ C f¯ F 0, F . We also assume that the problem

⎧ ¯ ∈ [ H 01 ( S )]d such that ⎪ find u ⎪ ⎪ ⎨  ⎪ ∇ u¯ : ∇ w = f¯ S · w ∀w ∈ [ H 01 ( S )]d ⎪ ⎪ ⎩ S

(64)

S

¯ ∈ [ H 2− 2 ( S )]d , is H 2− 2 -regular for some 2 ∈ (0, 1), that is, for any f¯ S ∈ [ L 2 ( S )]d , the solution to problem (64) satisfies u 1 / 2 −

d 2 ∇ u¯ · n ∈ [ H ()] and

u¯ 2− 2 , S + ∇ u¯ · n 1/2− 2 , ≤ C f¯ S 0, S . The previous hypotheses on the regularity of problems (63) and (64) are equivalent to angle conditions on  F and  S , respectively; see [16] and [18]. Now, let h : [ L 2 ()]d → V˜ h denote the projection operator characterized by





ρ h w · zh = 

ρ w · zh ∀zh ∈ V˜ h , ∀ w ∈ [ L 2 ()]d .



In particular, by the definition of V˜ h ,

 div(h w) qh = 0

∀ qh ∈ Q h ,

∀ w ∈ [ L 2 ()]d .

F

Under the previous hypotheses, there exists a positive constant C , independent of h, such that

  w − h w1, ≤ C hr − wr +1, F + wr +1, S ,

(65)

for all w ∈ V˜ := {w ∈ [ H 01 ()]d ; div(w) = 0 in  F } with w| ∈ [ H r +1 ( F )]d and w| ∈ [ H r +1 ( S )]d , where r ∈ [0, k] and F S

:= max( 1 , 2 ) (see Theorem A.3 in [12]). In the next Theorem we obtain a priori estimates for the error in both the global velocity and the Cauchy stress tensor. Theorem 4. Assume that problem (63) is H 2− 1 -regular for some 1 ∈ (0, 1) and that problem (64) is H 2− 2 -regular for some 2 ∈ (0, 1). Assume further that v(t )| F ∈ [ H r +1 ( F )]d , v(t )| S ∈ [ H r +1 ( S )]d , σ S (t ) ∈ [ H r ( S )]d×d and p (t ) ∈ H r ( F ) for t ∈ [0, T ]. Then, there exists C > 0, independent of h, such that

v(t ) − vh (t )20, + σ S (t ) − σ h (t ))20, S ≤ C h2(r − ) , where := max( 1 , 2 ). Proof. Let us consider the error equations:

⎧    d ⎪ ⎪ ρ (v(t ) − vh (t ))· wh + 2ν (ε (v(t )) − ε(vh (t ))): ε (wh ) + (σ S (t ) − σ h (t )): ε (wh ) ⎪ ⎪ dt ⎪ ⎪ ⎪  F S ⎪  ⎪ ⎪ ⎪ ⎪ ⎪ − ( p (t ) − ph (t )) div(wh ) = 0 ⎪ ⎪ ⎨  F   d ⎪ ⎪ (σ S (t ) − σ h (t )): τ h = λ div(v(t ) − vh (t )) tr(τ h ) + 2μ ε (v(t ) − vh (t )): τ h ⎪ ⎪ ⎪ dt ⎪ ⎪  S S S ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ qh div(v(t ) − vh (t )) = 0 ⎪ ⎪ ⎩ F

for any (wh , τ h , qh ) ∈ V h × S h × Q h and a.e. t ∈ (0, T ].

(66)

292

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

We first derive an estimate for the error in the approximation of the global velocity. From the first equation in (66), with wh = h v(t ) − vh (t ) ∈ V h , we have

1



ρ ∂t |v(t ) − vh (t )|2 + 2ν ε (v(t ) − vh (t ))20, F =

2







ρ ∂t (v(t ) − vh (t )) · (v(t ) − h v(t )) + 2ν

=





+  =

ε (v(t ) − vh (t )): ε (v(t ) − h v(t )) 

F

ρ ∂t (v(t ) − vh (t )) · (h v(t ) − vh (t )) + 2ν ε (v(t ) − vh (t )): ε (h v(t ) − vh (t )) = 



(67)

F

ρ ∂t (v(t ) − vh (t )) · (v(t ) − h v(t )) + 2ν ε (v(t ) − vh (t )): ε (v(t ) − h v(t )) 



+

 F

( p (t ) − ph (t )) div(h v(t ) − vh (t )) − F

(σ S (t ) − σ h (t )): ε (h v(t ) − vh (t )) . S

To deal with the last term on the right hand side of (67), we integrate the second equation of (66) in [0, t ], use the initial conditions and take τ h = ε ((h v(t ) − vh (t ))| S ) ∈ S h . Substituting in (67), we obtain that:

1



ρ ∂t |v(t ) − vh (t )|2 + 2ν ε (v(t ) − vh (t ))20, F =

2





=



ρ ∂t (v(t ) − vh (t )) · (v(t ) − h v(t )) + 2ν ε (v(t ) − vh (t )): ε (v(t ) − h v(t )) 



+

 F ( p (t ) − ph (t )) div(h v(t ) − vh (t )) − (σ 0S − σ h0 ): ε (h v(t ) − vh (t ))

F

S

t − λ div( (v(s) − vh (s)) ds) div(h v(t ) − vh (t )) 

S



0

t − 2μ ε ( (v(s) − vh (s)) ds): ε (h v(t ) − vh (t )) . S

0

Notice that ∂t (vh (t ) − h v(t )) ∈ V˜ h , so that





ρ ∂t (v(t ) − vh (t )) · (v(t ) − h v(t )) = 

ρ ∂t (v(t ) − h v(t )) · (v(t ) − h v(t )) , 

and we have that

1



ρ ∂t |v(t ) − vh (t )|2 + 2ν ε (v(t ) − vh (t ))20, F =

2 

=

1d





ρ |v(t ) − h v(t )|2 + 2ν ε (v(t ) − vh (t )): ε (v(t ) − h v(t ))

2 dt 

F



+

 (σ 0S − σ h0 ): ε (h v(t ) − vh (t ))

( p (t ) − qh ) div(h v(t ) − vh (t )) − F

S



t − λ div( (v(s) − vh (s)) ds) div(h v(t ) − v(t )) 0

S



t − 2μ ε ( (v(s) − vh (s)) ds): ε (h v(t ) − v(t )) S

0

(68)

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

293

t  − λ div( (v(s) − vh (s)) ds) div(v(t ) − vh (t )) 0

S



t − 2μ ε ( (v(s) − vh (s)) ds): ε (v(t ) − vh (t )) . 0

S

We can rewrite the penultimate term using that div(v(t ) − vh (t )) =

d

t div(v(s) − vh (s)) ds, so that

dt 0



t  t 1d div( (v(s) − vh (s)) ds) div(v(t ) − vh (t )) = |div( (v(s) − vh (s)) ds)|2 . 2 dt

0

S

0

S

Rewriting the last term in (68) similarly, integrating (68) in [0, t ] and using the initial conditions, we deduce that

ρmin 2

t v(t ) − vh (t )20,

ε (v(s) − vh (s))20, F ds +

+ 2ν 0

t ρmin 0 ρmax v − vh0 20, + v(t ) − h v(t )20, + με ( (v(s) − vh (s)) ds)20, S ≤ 2

0

−

ρmax 2

v0 − h v0 20, + 2 ν

ε (v(s) − vh (s)) : ε (v(s) − h v(s)) ds 0 F

t +

2

t 

 0 S

( p (s) − qh )div(h v(s) − vh (s)) ds − 0 F t

(σ − σ

0 h)

t : ε ( (h v(s) − vh (s))) ds 0

S

s div( (v( z) − vh ( z)) dz)div(h v(s) − v(s)) ds

−λ

0 S t

− 2μ

0

s ε( (v(z) − vh (z)) dz) : ε (h v(s) − v(s)) ds .

0 S

0

Now, using the Cauchy–Schwarz inequality and Young’s inequality repeatedly, we obtain

ρmin 2

v(t ) − vh (t )20,

+

ν

t ε (v(s) − vh (s))20, F ds

2 0

+

μ 2

t ρmin 0 ρmax ε ( (v(s) − vh (s)) ds)20, S ≤ v − vh0 20, + v(t ) − h v(t )20, 2

0

−

ρmax 2

v

0

− h v0 20,

+

3ν

ε

2 0

+

1

μ

λ 3μ σ 0S − σ h0 20, S + ( d + )

λ + ( d + μ)

2

t

2

2

t (v(s) − h v(s))20, F

+

d

t  p (s) − qh 20, F

ν 0

t

ε (h v(s) − v(s))20, S ds

2

0

s ε ( (v( z) − vh ( z)) dz)20, S ds .

0

0

Thus, assuming that p (t ) ∈ H ( F ) for t ∈ [0, T ], using the approximation properties (65) of the projection operator h , (21) and the approximation properties (25) and (27) of the initial data, we have r

294

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

ρmin 2

v(t ) − vh (t )20, +

ν

t ε (v(s) − vh (s))20, F ds +

2 0

μ 2

t ε ( (v(s) − vh (s)) ds)20, S ≤

t

s λd ≤ C h2(r − ) + ( + μ) ε ( (v( z) − vh ( z)) dz)20, S ds .

0

(69)

2

0

0

In particular,

t t  s  λd 2 2(r − ) ε ( (v(s) − vh (s)) ds)0, S ≤ C h + 2+ ε ( (v( z) − vh ( z)) dz)20, S ds .

μ

0

0

0

Therefore, by Gronwall’s lemma (31), we deduce that

t λd + 2μ λd+μ2μ t ε ( (v(s) − vh (s)) ds)20, S ≤ C h2(r − ) (1 + te ) ≤ C h2(r − ) .

μ

(70)

0

Using this bound in (69), we have

ρmin 2

v(t ) − vh (t )20, +

ν

t

2

ε (v(s) − vh (s))20, F ds ≤ C h2(r − ) .

0

×d →  h : [ L 2 ( S )]dsym Now we obtain a bound for the error in the approximation of the Cauchy stress tensor. To this end, let 

S h be the orthogonal projection operator with respect to the usual [ L 2 ( S )]d×d -inner product. In particular, for any ×d , τ ∈ [ L 2 ( S )]dsym

 h τ 0, S ≤ C hr . τ − 

(71)

Then, we can write

  h σ S (t )):(σ S (t ) −   h σ S (t )) σ S (t ) − σ h (t )20, S = (σ S (t ) −  S    h σ S (t ) − σ h (t )):(σ S (t ) −   h σ S (t )) + (σ S (t ) − σ h (t )):(  h σ S (t ) − σ h (t )) . + ( S

S

Besides, integrating in [0, t ] the second equation of (66) and taking



S

 h σ S (t ) − σ h (t ) ∈ S h , τh = 

 h σ S (t ) − σ h (t )) = (σ S (t ) − σ h (t )):(  0 S

= (σ − σ S

0  h ):(h

t  h σ S (t ) − σ h (t )) σ S (t ) − σ h (t )) + λ div( (v(s) − vh (s))ds) tr( 

S



0

t  h σ S (t ) − σ h (t )) . + 2μ ε ( (v(s) − vh (s))ds):( 0

S

Consequently,

 h σ S (t )20, + σ S (t ) − σ h (t )20, S = σ S (t ) −  S  0 S

+ (σ − σ S

0  h ):(h

 h σ S (t ) − σ h (t )):(σ S (t ) −   h σ S (t )) (

S

t  h σ S (t ) − σ h (t )) σ S (t ) − σ h (t )) + λ div( (v(s) − vh (s))ds) tr( 

S





0

t  h σ S (t ) − σ h (t )) . + 2μ ε ( (v(s) − vh (s))ds):( S

0

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

295

Using the Cauchy–Schwarz inequality, (71), (27), (70), Young’s inequality and the triangle inequality, we conclude

σ S (t ) − σ h (t ))20, S ≤ C h2(r − ) . 2 5. A fully discrete scheme. Numerical experiments In this section we use the backward Euler method to deduce a fully discrete scheme associated to the semidiscrete problem (29). We also give examples of finite element subspaces satisfying the hypotheses of the previous section. Finally, we provide some numerical results. 5.1. Fully discrete scheme Let N be a given positive integer and let {t i }iN=0 be a uniform partition of [0, T ] of mesh size t := T / N. Then, for n = 1, 2, . . . , N, we approximate

∂t vh (tn ) ≈

vh (tn ) − vh (tn−1 )

t

,

∂t σ h (tn ) ≈

σ h (tn ) − σ h (tn−1 ) t

We arrive at the following fully discrete scheme: Given vh0 ∈ V h and S h × Q h such that

.

σ h0 ∈ S h , for n = 1, 2, . . . , N, find (vnh , σ nh , pnh ) ∈ V h ×

  ⎧ n n ⎪ ⎪ ρ v · w + 2 ν  t ε ( v ): ε ( w ) +  t σ nh : ε (wh ) h h ⎪ h h ⎪ ⎪ ⎪ ⎪  F S  ⎪   ⎪ ⎪ ⎪ n ⎪ −t ph div(wh ) = t f(tn ) · wh + ρ vnh−1 · wh ∀ wh ∈ V h ⎪ ⎪ ⎪ ⎪ ⎪ F   ⎪   ⎪ ⎪ ⎨ n n σ h : τ h = λ t div(vh ) tr(τ h ) ⎪ ⎪ S ⎪  S  ⎪ ⎪ ⎪ n ⎪ + 2 μ t ε (vh ): τ h + σ nh−1 : τ h ∀ τ h ∈ S h ⎪ ⎪ ⎪ ⎪ ⎪ S S ⎪  ⎪ ⎪ ⎪ n ⎪ ⎪ qh div(vh ) = 0 ∀ qh ∈ Q h ⎪ ⎪ ⎩ F

where for n = 1, 2, . . . , N, vnh , σ nh and pnh are approximations of vh (tn ), σ h (tn ) and p h (tn ), respectively. Now, we observe that taking τ h = ε (wh | S ), it is possible to decouple at each time step the computation of vnh and pnh from that of σ nh . That is, in each time step, we first compute (vnh , pnh ) by solving a problem with a saddle point structure. Then, we compute σ nh by solving a simple equation: 1. Given vnh−1 ∈ V h and

σ nh−1 ∈ S h , find (vnh , pnh ) ∈ V h × Q h such that   ⎧ ⎪ n n 2 ⎪ ρ v · w + 2 ν  t ε ( v ): ε ( w ) + 2 μ  t ε (vnh ): ε (wh ) ⎪ h h h h ⎪ ⎪ ⎪ ⎪  F S ⎪   ⎪ ⎪ ⎪ 2 n n ⎪ +λ  t div ( v ) div ( w ) −  t p div ( w ⎪ h h) = ⎪ h h ⎪ ⎨ S F    ⎪ n −1 ⎪ = t f(tn ) · wh + ρ vh · wh − t σ nh−1 : ε (wh ) ∀ wh ∈ V h ⎪ ⎪ ⎪ ⎪ ⎪  S ⎪   ⎪ ⎪ ⎪ n ⎪ ⎪ t qh div(vh ) = 0 ∀ qh ∈ Q h ⎪ ⎪ ⎩ F

2. Given

vnh



∈ V h and σ nh−1 ∈ S h , compute σ nh ∈ S h such that for any τ h ∈ S h



σ nh : τ h = t S

S



 λ div(vnh ) I + 2 μ ε (vnh ) : τ h +



σ nh−1 : τ h

S

It is clear that, given vnh ∈ V h and σ nh−1 ∈ S h , the second problem has a unique solution problem, we remark that, thanks to Korn’s inequality, the bilinear form

σ nh ∈ S h . Concerning the first

296

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299







ρ v · w + 2 ν t ε (v): ε (w) + 2 μ t 2

a(v, w) := 

F



ε (v): ε (w) + λ t 2 S

div(v) div(w) S

is elliptic in [ H 01 ()]d and, in particular, on V h . On the other hand, using the global inf–sup condition (20) and the Babuska– Brezzi theory (see [8,16]), we deduce that the first problem has a unique solution (vnh , pnh ) ∈ V h × Q h . 5.2. Examples of finite element subspaces Next we describe some finite element subspaces that can be used to solve the fully discrete scheme in practice. We consider continuous approximations of the pressures. In what follows, we assume that Th is regular. Moreover, given an element T ∈ Th and an integer l ≥ 0, we denote by Pl ( T ) the space of polynomials of total degree at most l defined on T . 5.2.1. A first order method based on the mini element Let d = 2 and denote by λi , i = 1, 2, 3, the i-th barycentric coordinate of a triangle T ∈ Th . First, we consider a finite element method based on the mini-element introduced by Arnold, Brezzi and Fortin in [2]. Indeed, given a triangle T ∈ Th , we define



P˜ 1 ( T ) :=

P1 ( T ) ⊕ λ1 λ2 λ3  if T ⊂  F , if T ⊂  S .

P1 ( T )

Then, we choose the following finite element spaces:

¯ ) ∩ H 1 ()]2 ; vh |T ∈ [P˜ 1 ( T )]2 , V h := {vh ∈ [C 0 ( 0 ¯ F ) ; q h | T ∈ P1 ( T ) , Q h := {qh ∈ C ( 0

×2 ; S h := {τ h ∈ [ L 2 ( S )]2sym

∀ T ∈ Th } ,

∀ T ∈ Th | F } ,

(72) (73)

τ h |T ∈ [P0 ( T )]2×2 } .

(74)

We recall that the degrees of freedom of the mini-element are the values of the velocity at the vertices and center of a triangle T ∈ Th | F and the values of the pressure at the vertices of T ∈ Th | F . Then, the finite element subspace V h is well-defined. Moreover, the pair ( V h , Q h ) satisfies the inf–sup condition (19) (see [2,16]). Finally, we also remark that the space S h satisfies condition (18). 5.2.2. An example based on the generalized Hood–Taylor finite element Now we consider the generalization of the Hood–Taylor finite element (cf. [19]) to higher orders in two and three space dimensions. Let k ≥ 1 be a natural number. We choose the following finite element spaces:

¯ ) ∩ H 01 ()]d ; vh |T ∈ [Pk+1 ( T )]d , V h := {vh ∈ [C 0 ( ¯ F ) ; qh |T ∈ Pk ( T ) , Q h := {qh ∈ C ( 0

d×d

2

S h := {τ h ∈ [ L ( S )]sym ;

∀ T ∈ Th } ,

∀ T ∈ Th | F } ,

τ h |T ∈ [Pk ( T )]d×d } .

Here again, the finite element pair ( V h , Q h ) satisfies the inf–sup condition (19) (see [16,6,7]), and the space S h satisfies condition (18). 5.3. Numerical experiments In this subsection, we provide some numerical results. We implemented the fully discrete scheme described in subsection 5.1 in a MatLab code using the first order method based on the mini element, that is, we used the finite element spaces (72), (73) and (74). We remark that this approach is easy to implement. For the sake of simplicity, here we just provide some results for a simple test and refer to [17] for some more examples and details about the implementation. We let  F := {x ∈ R2 ; (x1 − 0.2)2 + (x2 − 0.2)2 < 0.0025} and  S := ((0, 2.5) × (0, 0.41)) \  F . We remark that in this example  F = ∅. Moreover, we take T = 1, ρ F = ρ S = 1, ν = 12 , λ = 1 and μ = 12 , and choose f and the initial data so that the exact solution is

v 1 (t , x) = 2 x1 x2 et , v 2 (t , x) =

(x21

−

x22 )

in (0, T ) × , t

e,

in (0, T ) × ,

p (t , x) = (|x − (0.2, 0.2)| − 0.05 ) e , 2

2

t

(75) in (0, T ) ×  F .

Concerning the approximation of the initial data, we need to define vh0 ∈ V h and σ h0 ∈ S h . Instead of solving problems (22), (23) and (26), we consider continuous piecewise linear approximations of the initial velocity and of the displacements (in order to initialize the stress tensor).

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

297

Table 1 Errors for successive meshes and the fixed time step t = 10−1 . h

Nh

dof

eht ( v 1 )

eht ( v 2 )

eht ( p )

eht (σ )

eht (v, p , σ )

l l/2 l/4 l/8 l/16

128 471 1805 7065 27953

205 879 3643 14835 59875

0.322069 0.156695 0.078325 0.041833 0.026045

0.344253 0.161394 0.081523 0.046970 0.033670

0.006639 0.001317 0.000478 0.000288 0.000231

0.879222 0.441115 0.222592 0.114478 0.063047

0.997655 0.495163 0.249656 0.130620 0.076072

Table 2 Errors for successive meshes and the fixed time step t = 10−2 . h

Nh

dof

eht ( v 1 )

eht ( v 2 )

eht ( p )

eht (σ )

eht (v, p , σ )

l l/2 l/4 l/8 l/16

128 471 1805 7065 27953

205 879 3643 14835 59875

0.348366 0.165381 0.077098 0.037246 0.018558

0.426616 0.172998 0.077168 0.037030 0.018822

0.006181 0.001249 0.000420 0.000186 0.000075

0.841012 0.421156 0.211057 0.105666 0.052931

1.005336 0.484410 0.237580 0.118000 0.059164

Table 3 Errors for successive meshes and the fixed time step t = 10−3 . h

Nh

dof

eht ( v 1 )

eht ( v 2 )

eht ( p )

eht (σ )

eht (v, p , σ )

l l/2 l/4 l/8 l/16

128 471 1805 7065 27953

205 879 3643 14835 59875

0.501274 0.237086 0.097198 0.042409 0.019519

0.639969 0.251344 0.097492 0.041071 0.018816

0.006380 0.001297 0.000427 0.000186 0.000072

0.837936 0.419378 0.210107 0.105146 0.052591

1.167481 0.543381 0.251192 0.120586 0.059167

Table 4 Errors for successive meshes and the fixed time step t = 10−4 . h

Nh

dof

eht ( v 1 )

eht ( v 2 )

eht ( p )

eht (σ )

eht (v, p , σ )

l l/2 l/4 l/8 l/16

128 471 1805 7065 27953

205 879 3643 14835 59875

0.650843 0.322121 0.139601 0.061127 0.026386

0.789583 0.344105 0.142698 0.059762 0.024756

0.006692 0.001384 0.000445 0.000188 0.000072

0.837875 0.419269 0.210028 0.105101 0.052566

1.322542 0.630840 0.289764 0.135478 0.063815

We build a sequence of meshes by uniformly refining an initial mesh; we denote by h = l, l/2, . . . , l/16 the corresponding mesh-sizes, N h is the number of vertices, and dof stands for the degrees of freedom in the global velocity and fluid pressure unknowns. To determine the order of convergence, we compare the computed solutions (vht , p ht , σ ht ) for a fixed time step t, with the exact solution (v, p , σ ) in the corresponding norms: t eht ( v j ) := || v j − v  j ,h || L 2 ((0, T ), H 1 ())

t

j = 1, 2,

t

eh ( p ) := || p − p h || L 2 ((0, T ), L 2 ( F )) , eht (σ ) := ||σ − σ ht || L 2 ((0, T ), L 2 ( S )2×2 ) ,

1/2  eht (v, p , σ ) := eht ( v 1 )2 + eht ( v 2 )2 + eht ( p )2 + eht (σ )2 . The experimental convergence rate is then computed for each pair of consecutive meshes as the slope in log–log scale of the error versus the number of vertices of the associated mesh:

rh := −2

log eh/2 − log eh log N h/2 − log N h

.

(76)

In Tables 1–4 we show the total number of vertices, the corresponding dof and the errors for each unknown obtained with the time steps t = 10− j ( j = 1, . . . , 4), respectively. The corresponding experimental convergence rates for each unknown are shown in Tables 5–8. We observe in Table 5 that, for t = 10−1 and the coarse meshes, experimental convergence rates are around 1 for the two components of the velocity and the stress tensor; however, experimental convergence rates degrade as the meshes become finer, specially for the pressure variable.

298

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

Table 5 Experimental convergence rates for the fixed time step t = 10−1 . rh ( v 1 ) rh ( v 2 ) rh ( p ) rh (σ ) rh (v, p , σ )

1.105997 1.162899 22.483089 1.058822

1.032317 1.016717 1.508998 1.018218

0.919218 0.808115 0.744440 0.974589

0.689075 0.484083 0.321952 0.867403

1.075385

1.019461

0.949434

0.786120

Table 6 Experimental convergence rates for the fixed time step t = 10−2 . rh ( v 1 ) rh ( v 2 ) rh ( p ) rh (σ )

1.143671 1.385613 2.454835 1.061695

1.136135 1.201811 1.622782 1.028504

1.066281 1.076159 1.194629 1.013989

1.013050 0.984037 1.324761 1.005257

rh (v, p , σ )

1.120863

1.060588

1.025686

1.003911

Table 7 Experimental convergence rates for the fixed time step t = 10−3 . rh ( v 1 ) rh ( v 2 ) rh ( p ) rh (σ )

1.149389 1.434720 2.445925 1.062566

1.327431 1.409867 1.654835 1.028922

1.215576 1.267008 1.220486 1.014612

1.128409 1.135122 1.384567 1.007453

rh (v, p , σ )

1.174052

1.148670

1.075559

1.035347

Table 8 Experimental convergence rates for the fixed time step t = 10−4 . rh ( v 1 ) rh ( v 2 ) rh ( p ) rh (σ )

1.079714 1.275007 2.418891 1.062850

1.244754 1.310380 1.688830 1.029092

1.210386 1.275632 1.262575 1.014693

1.221631 1.281573 1.397115 1.007503

rh (v, p , σ )

1.136387

1.158183

1.114267

1.094722

Nevertheless, this issue can be overcome just by taking a smaller time step t. Indeed, we can observe in Tables 6–8 that for smaller t, optimal convergence rates are attained with respect to the mesh size h for all the unknowns even when we use the finest meshes. In particular, we observe that in this example, the pressure variable converges superlinearly. On the other hand, for small enough time steps t, the error is dominated by the mesh size h, as can be seen comparing the rows of Tables 2–4 which correspond to the absolute errors for each fixed mesh size h = l, . . . , l/16. We observed a similar behavior in other tests. 6. Conclusions In summary, we have introduced a new formulation of a fluid–structure interaction problem in terms of the global velocity, the fluid pressure and the Cauchy stress tensor. Displacements in the solid can be recovered by quadrature. We have shown the existence and uniqueness of a solution to the new weak formulation, a priori estimates for the solution of a semidiscrete problem, convergence and error estimates provided that the finite element pair in the fluid domain is stable, compatible finite elements are used in the solid and the domains satisfy certain regularity conditions. We have proposed a fully discrete scheme based on the backward Euler method. This scheme involves solving at each time step a saddle point problem for the velocity and the pressure; the structure stress tensor is computed by solving a simple equation. We implemented the method using the mini-element in the fluid domain, continuous piecewise affine functions to approximate the structural velocity and piecewise constant symmetric tensors to approximate the Cauchy stress tensor. Numerical experiments show that the method is robust, and optimal convergence rates can be attained with respect to the mesh size if the time discretization is fine enough. The derivation of this new formulation depends on the linear nature of the constitutive relation in the solid domain. In the case of a non-linear constitutive law in the solid, it seems that one cannot eliminate the displacements from the formulation. This type of problems will be addressed in a future work. We also remark that the analysis of the semidiscrete problem presented in this paper is valid only for a class of initial discrete fluid velocities satisfying the relation (22), which involves the initial fluid pressure. However, for the implementation of the fully discrete scheme described in subsection 5.1, we considered continuous piecewise linear approximations of the initial velocity and displacements, avoiding the use of any initial pressure.

M. González, V. Selgas / Applied Numerical Mathematics 123 (2018) 275–299

299

Acknowledgements The authors of the manuscript would like to thank the anonymous referees for their valuable comments. This research has been partially supported by MICINN grant MTM2010-21135-C02-01. Besides, the research of M. González was partially supported by MINECO grant MTM2013-47800-C2-1-P while the research of V. Selgas was partially supported by MINECO grant MTM2013-43671-P. References [1] R.A. Adams, J.J.F. Fournier, Sobolev Spaces, second edition, Pure and Applied Mathematics, Elsevier/Academic Press, Amsterdam, 2003. [2] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337–344. [3] M. Astorino, C. Grandmont, Convergence analysis of a projection semi-implicit scheme for fluid–structure interaction problems, Numer. Math. 116 (2010) 721–767. [4] Y. Bazilevs, K. Takizawa, T.E. Tezduyar, Computational Fluid-Structure Interaction, Methods and Applications, Wiley, 2013. [5] T. Bodnár, G.P. Galdi, S. Neˇcasová (Eds.), Fluid-Structure Interaction and Biomedical Applications, Birkhäuser, 2014. [6] D. Boffi, Stability of higher-order triangular Hood–Taylor methods for the stationary Stokes equations, Math. Models Methods Appl. Sci. 4 (1994) 223–235. [7] D. Boffi, Three-dimensional finite element methods for the Stokes problem, SIAM J. Numer. Anal. 34 (1997) 664–670. [8] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, 1991. [9] E. Burman, M.A. Fernández, An unfitted Nitsche method for incompressible fluid–structure interaction using overlapping meshes, Comput. Methods Appl. Mech. Eng. 279 (2014) 497–514. [10] S. Dasser, A penalization method for the homogenization of a mixed fluid–structure problem, C. R. Acad. Sci. Paris, Ser. I Math. 320 (1995) 759–764. [11] Q. Du, M.D. Gunzburger, L.S. Hou, J. Lee, Analysis of a linear fluid–structure interaction problem, Discrete Contin. Dyn. Syst. 9 (3) (2003) 633–650. [12] Q. Du, M.D. Gunzburger, L.S. Hou, J. Lee, Semidiscrete finite element approximations of a linear fluid–structure interaction problem, SIAM J. Numer. Anal. 42 (2004) 1–29. [13] L.C. Evans, Partial Differential Equations, AMS Graduate Studies in Mathematics, vol. 19, 1994. [14] M.A. Fernández, Incremental displacement correction schemes for incompressible fluid–structure interaction. Stability and convergence analysis, Numer. Math. 123 (2013) 21–65. [15] M.A. Fernández, J. Mullaert, Convergence and error analysis for a class of splitting schemes in incompressible fluid–structure interaction, IMA J. Numer. Anal. 36 (2016) 1748–1782. [16] V. Girault, P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations, Theory and Algorithms, Springer Series in Computational Mathematics, Springer-Verlag, 1986. [17] M. González-Taboada, V. Selgas, A new numerical scheme for a linear fluid–structure interaction problem, in: IV International Conference on Computational Methods for Coupled Problems in Science and Engineering, International Center for Numerical Methods in Engineering, 2011. [18] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, 1985. [19] P. Hood, C. Taylor, A numerical solution of the Navier–Stokes equations using the finite element technique, Comput. Fluids 1 (1973) 73–100. [20] O. Kayser-Herold, H.G. Matthies, A unified least-squares formulation for fluid–structure interaction problems, Comput. Struct. 85 (2007) 998–1011. [21] A. Kufner, O. John, S. Fucík, Function Spaces, Academia, Praga, 1977. [22] P. Le Tallec, S. Mani, Numerical analysis of a linearised fluid–structure interaction problem, Numer. Math. 87 (2000) 317–354. [23] R. Schulkes, Interactions of an elastic solid with a viscous fluid: eigenmode analysis, J. Comput. Phys. 100 (1992) 270–283. [24] M. Souli, D.J. Benson (Eds.), Arbitrary Lagrangian–Eulerian and Fluid–Structure Interaction. Numerical Simulation, Wiley, 2010. [25] C. Surulescu, On a time-dependent fluid–solid coupling in 3D with nonstandard boundary conditions. Steps towards modeling blood flow, Acta Appl. Math. 110 (2010) 1087–1104.