DECOUPLED TIME STEPPING METHODS FOR ... - Semantic Scholar
DECOUPLED TIME STEPPING METHODS FOR FLUID-FLUID INTERACTION JEFFREY M. CONNORS⇤ , JASON S. HOWELL† , AND WILLIAM J. LAYTON‡ Abstract. A model of two incompressible, Newtonian fluids coupled across a common interface is studied. The nonlinearity of the coupling condition exacerbates the problem of decoupling the fluid calculations in each subdomain, a natural parallelization strategy employed in current climate models. A specialized partitioned time stepping method is studied which decouples the discrete fluid equations without sacrificing stability and maintaining convergence. This is accomplished through explicit updating of the size of the jump in tangential velocities across the fluid-fluid interface by a geometric averaging of this data over the previous two time levels. A full numerical analysis is presented and computational tests are performed demonstrating the robustness of this method. Key words. semi-implicit, fluid-fluid interaction, atmosphere-ocean, implicit-explicit method.
1. Introduction. The dynamic core in atmosphere-ocean interaction is a critical component of climatology models and has attracted the interest of many mathematical researchers due to the richness of the theory and technical complexities of efficiently computing approximations to the coupled system using only (uncoupled) atmosphere and ocean solves, (see e.g. [4, 6, 17, 18, 19]). Motivated by this problem we consider uncoupled or “partitioned” methods for two fluids coupled across their shared interface I by a rigid-lid coupling condition, i.e. no penetration and a slip with friction condition allowing a jump in the tangential velocities across I. This rigid-lid assumption is used in many oceanography models, [20, 21]. Partitioned methods allow existing, highly optimized codes for each subproblem to be used in parallel as black boxes at each time step to solve the coupled problem. The coupling between the atmosphere and the ocean dynamics occurs in the nonlinear interface condition (1.2) below. It is not clear (see [9]) if the nonlinearity in (1.2) can be treated in the most natural way (which is stable if (1.2) is linear). We give in (2.8)-(2.9) below an uncoupled and surprising (to us) discretization for the nonlinear condition (1.2) which has no analog for the linear version of (1.2) and which yields an unconditionally stable partitioned method for the fully coupled problem (1.1) - (1.6) below. The modification is based on using the geometric averaging at two time levels of the slip velocity at the interface to compute the friction coefficient rather than simply lagging it. To reduce the dynamic core of the coupled atmosphere-ocean problem to its simplest form which still retains the essential difficulty of the coupling condition, let the domain consist of two subdomains ⌦1 and ⌦2 of Rd , d = 2, 3 with outward unit normal vectors n ˆ 1 and n ˆ 2 , respectively, coupled across an interface I (example in Figure 1.1 below). The problem is: given ⌫i > 0, fi : [0, T ] ! H 1 (⌦i ), ui (0) 2 H 1 (⌦i ) and ⇤ Lawrence Livermore National Laboratory, Center for Applied Scientific Computing, Livermore, CA, 94550. Partially supported by NSF Grants DMS 0508260 and 0810385. This work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344, LLNL-JRNL-490264. email: [email protected] † Department of Mathematics and Computer Science, Clarkson University, Potsdam, NY 136995815. email: [email protected]. This material is partially based upon work supported by the Center for Nonlinear Analysis (CNA) under the National Science Foundation Grant No. DMS– 0635983. ‡ Department of Mathematics, University of Pittsburgh, Pittsburgh, PA, 15260, USA. email: [email protected]. Partially supported by NSF Grants DMS 0508260 and 0810385 .
1
2 R, find (for i = 1, 2) ui : ⌦i ⇥ [0, T ] ! Rd and pi : ⌦i ⇥ [0, T ] ! R satisfying ui,t + ui · rui
⌫i ui + rpi = fi ,
in ⌦i ,
⌫i n ˆ i · rui · ⌧ = |ui uj |(ui uj ) · ⌧, on I , for i, j = 1, 2 , i = 6 j, ui · n ˆ i = 0 on I, i = 1, 2 r · ui = 0 , in ⌦i
ui (x, 0) = u0i (x), ui = 0,
on
in ⌦i , i
= @⌦i \ I.
(1.1) (1.2) (1.3) (1.4) (1.5) (1.6)
The lateral boundary conditions on ⌦i , (1.6), are not essential for our study. We focus on (1.2). Let Xi : = {vi 2 H 1 (⌦i )d : vi = 0 on i = @⌦i \ I, vi · n ˆ i = 0 on I} ⇢ Z qi d⌦i = 0 . Qi : = qi 2 L2 (⌦i )d :
(1.7) (1.8)
⌦i
For ui 2 Xi we denote u = (u1 , u2 ) and X := {v = (v1 , v2 ) : vi 2 Xi , i = 1, 2}. Similarly, for qi 2 Qi we denote q = (q1 , q2 ) and Q := {q = (q1 , q2 ) : qi 2 Qi , i = 1, 2}. A natural subdomain variational formulation for (1.1)-(1.6), obtained by multiplying (1.1) by vi and (1.4) by qi , integrating and applying the divergence theorem, is to find (for i, j = 1, 2, i 6= j) ui : [0, T ] ! Xi and pi : [0, T ] ! Qi satisfying (ui,t , vi )⌦i + (⌫i rui , rvi )⌦i + (ui · rui , vi )⌦i (pi , r · vi )⌦i Z + |ui uj |(ui uj )vi ds = (fi , vi )⌦i , 8 vi 2 Xi
(1.9)
I
(r · ui , qi )⌦i = 0 , 8qi 2 Qi .
Let [·] denote the jump of the indicated quantity across the interface I , (·, ·) the L2 (⌦1 [ ⌦2 ) inner product and ⌫ = ⌫i and f = fi in ⌦i . The natural coupled or monolithic variational formulation for (1.1)-(1.6) is found by summing (1.9) over i, j = 1, 2 and i 6= j and is to find u : [0, T ] ! X and p : [0, T ] ! Q satisfying (ut , v) + (⌫ru, rv) + (u · ru, v) (p, r · v) Z + |[u]| [u][v]ds = (f , v), 8v 2 X,
(1.10)
I
(r · u, q) = 0 , 8q 2 Q.
Figure 1.1 illustrates the subdomains considered here, sometimes used in tests of fluid-fluid interaction, [6]. Comparing (1.10) and (1.9) we see that the monolithic problem (1.10) has a global energy that is exactly conserved, (in the appropriate sense), (set v = u and q = p in (1.10)). The subdomain sub-problems (1.9) do not posses a subdomain energy which behaves similarly due to energy transfer back and forth across the interface I. It is possible for decoupling strategies to be unstable due to the input of non-physical energy as a numerical artifact. We note that a realistic atmosphere-ocean model would contain many more terms including (typically) some eddy-viscocity term representing the e↵ects of unresolved turbulent fluctuations on resolved scales. Without such a term there is an interesting theoretical question as to whether weak solutions of (1.1) - (1.6) possess enough regularity for their traces to be well defined in L3 (I). We focus on algorithmic issues herein and suppose that the solution to (1.1) - (1.6) approximated is a strong solution. 2
n ˆ2 6
⌦1
? n ˆ1
⌦2
Fig. 1.1. Example adjoining subdomains.
1.1. Partitioned stability. The simplest and most natural partitioned method for (1.1) - (1.6) is to extrapolate from past time values the speed |ui uj | in the friction coefficient, pass the previous, known interface velocity uj across I to the problem in ⌦i and solve the problem in ⌦i . In the case of a linear coupling (without |ui uj | on the right hand side of (1.2)), this is known to be unconditionally stable, [9]; in the nonlinear case it is given as follows. For final time T and time step size t = T /N for some N 2 N, find un+1 2 Xih ⇢ Xi , pn+1 2 Qhi ⇢ Qi for i, j = 1, 2, i 6= j satisfying i i ✓
un+1 i
uni t
(pin+1 , r
, vi
◆
⌦i
+ (⌫i run+1 , rvi )⌦i + (un+1 · run+1 , vi )⌦i i i i
· vi )⌦i +
Z
I
|uni
unj |(un+1 i
unj )vi ds = (fi (tn+1 ), vi )⌦i , 8 vi 2 Xih (r · un+1 , qi )⌦i = 0 , 8qi 2 Qhi . (1.11) i
In the case of linear coupling (without |uni unj |), (1.11) has been proven in [9] to be unconditionally stable. Stability of (1.11) is an open question; while we cannot prove either stability or instability, in Section 4 we present one test where (1.11) exhibits some problems with energy dissipation. The development of a partitioned discretization method which is unconditionally and nonlinearly stable for large time steps for the nonlinear interface condition (1.2) is the main motivation of this report. We show that replacing the term unj |uni unj | in (1.11) by using a geometric averaging of the jump, unj |uni
unj |1/2 |uni
1
unj
1 1/2
|
,
gives this partitioned method. The remainder of this work is organized as follows: in Section 2, the mathematical setting and notation are described, and the partitioned time-stepping algorithm is presented. Theoretical results regarding the stability and convergence of the method are presented in Section 3. Unconditional stability of the proposed numerical method is proven and an optimal error estimate is derived. A remark is included discussing the expected dissipation of the decoupled algorithm. Numerical results that support the theoretical results are presented in section 4. Convergence rates are demonstrated for an academic test. Also, an example is provided where the unmodified IMEX discretization (1.11) exhibits a faster growth in global relative error than the geometric averaging method. 3
1.2. Related work. A standard simplification of the interface condition (1.5) used in some climate models, (e.g. [21]), is to represent the “wind shear stress” or “momentum flux” by a relationship of the form ⇢O µO ruO · n ˆ O = CD ⇢a |ua |ua where ⇢O , ⇢a are densities of the ocean and atmosphere, respectively, µO is the dynamic viscocity of the ocean, uO and ua are the velocities of the ocean and atmosphere, respectively, n ˆ O is the unit outward normal vector of the ocean domain and CD is a (dimensionless) drag coefficient. The velocity jump ua uO does not appear as it would using (1.2) since the ocean velocity is assumed to be negligible compared to the air velocity near the interface. The atmosphere and ocean are updated by an explicit update of the ocean to compute the momentum flux, which is then mapped conservatively to the atmosphere mesh. An atmosphere-ocean system, more physically refined than (1.1) - (1.6) but sharing the same mathematical structure for the dynamic core and coupling, is studied by Lions, Temam and Wang in a series of papers performing a thorough mathematical analysis [17, 18, 19]. The coupling condition (1.2) is a more physically relevant representation of the momentum flux across the interface, taking into account when the ocean surface flow may be influencing the surface air current. Bresch and Koko [6] introduce some decoupling algorithms for a model akin to (1.1)-(1.6). In it, the monolithic, coupled problem at each step is solved by preconditioned iterative methods, uncoupling can occur in the residual calculation and in the preconditioning step. Often systems of ODEs or PDEs may be decomposed into a sum of linear, sti↵ terms and other less sti↵ (possibly nonlinear) terms. Stable, consistent methods for such ODEs were derived by using implicit treatment of sti↵ terms and explicit treatment of the remaining terms, the so called additive Runge-Kutta methods, [10]. Kennedy and Carpenter exhibited the utility of this idea for convection-di↵usionreaction equations, [14], and in [3] implicit-explicit (IMEX) methods of up to fourth order have been developed for some classes of PDEs. These methods seek to reduce the time step restriction on an explicit method for stability by implicit treatment of sti↵ terms while reducing computational cost through explicit updating of other terms, particularly nonlinear terms. One application of the IMEX approach is to use explicit updates of the terms responsible for coupling across I, giving semi-implicit partitioned methods. A straightforward partitioning approach for some coupled problems may result in an unstable or conditionally stable method, as demonstrated for the fluid-structure interaction problem, see [7, 8]. 2. Method Description, Notation and Preliminaries. This section presents the numerical scheme for (1.1)-(1.6), and provides the necessary definitions and lemmas for the stability and convergence analysis. For D ⇢ ⌦, the Sobolev space H k (D) = W k,2 (D) is equipped with the usual norm k·kH k (D) , and semi-norm |·|H k (D) , for 1 k < 1, e.g. Adams [2]. The L2 norm is denoted by k·kD . For functions v(x, t) defined for almost every t 2 (0, T ) on a function space V (D), we define the norms (1 p 1) !1/p Z kvkL1 (0,T ;V ) = ess sup kv(·, t)kV
and
0 0, f 2 H 1 (⌦). Given n 1 u , un 2 Vh , n 2 {1, 2, · · · , N 1}, find un+1 2 V1,h satisfying 1 ✓
un+1 1
un1 t
, v1
◆
⌦1
+ ⌫1 run+1 , rv1 1 +
Z
I
⌦1
+ c1 un+1 ; un+1 , v1 1 1
|[un ]|un+1 · v1 ds 1
= f1 (tn+1 ), v1
⌦1
Z
I
|[un ]|1/2 |[un
1 1/2 n ]| u2
, 8v1 2 V1,h .
(2.8)
and un+1 2 V2,h satisfying 2 ✓
un+1 2
un2 t
, v2
◆
⌦2
+ ⌫2 run+1 , rv2 2 +
Z
I
⌦2
+ c2 un+1 ; un+1 , v2 2 2
|[un ]|un+1 · v2 ds 2
= f2 (tn+1 ), v2
⌦2
· v1 ds
Z
I
|[un ]|1/2 |[un
1 1/2 n ]| u1
· v2 ds
, 8v2 2 V2,h .
(2.9) The main feature of this algorithm is the use of the geometric average |[un ]|1/2 |[un 1 ]|1/2 in the coupling terms. This will be shown to have a stabilizing e↵ect for the algorithm. Hence Algorithm 2.1 will often be referred to as the geometric averaging method, or GA method. 3. Numerical Analysis. Stability of the partitioned scheme (Algorithm 2.1) is established below. The polarization identity is used in the stability and convergence proofs: (wj+1
wj ) · wj+1 =
1 |wj+1 |2 + |wj+1 2 7
wj |2
|wj |2
(3.1)
Lemma 3.1. (Partitioned Stability) Let uj1 2 X1,h satisfy (2.8) for each j 2 {0, 1, 2, · · · , n 1, n}, and uj2 2 X2,h satisfy (2.9) for each j 2 {0, 1, 2, · · · , n 1, n}. Then Algorithm 2.1 satisfies the following energy equality at time step n + 1: un+1
2
+
n X
uj+1
j=1
+
t
+
2
uj
t
+2 t
j=1
n Z X
j=1 I n Z X j=1
Z
n ✓ X
I
2
⌫1 ruj+1 1
⌦1
2
+ ⌫2 ruj+1 2
|[uj ]|1/2 uj+1 1
|[uj
1 1/2 j ]| u2
|[uj ]|1/2 uj+1 2
|[uj
1 1/2 j ]| u1
2
⌦2
◆
ds
2
ds
⇣ ⌘ 2 2 |[un ]| un+1 + un+1 ds 1 2 I Z ⇣ ⌘ 2 2 2 + t [u0 ] u11 + u12 ds
+ t = u1 +2 t
n ⇢⇣ X
I
f1 (tj+1 ), uj+1 1
j=1
⌘
⌦1
⇣ ⌘ + f2 (tj+1 ), uj+1 2
. (3.2)
⌦2
Furthermore, the algorithm is unconditionally stable, satisfying: n+1 2
u
+
n X
j 2
u
+
u
j+1
j=1
ruj+1 1
⌫1
j=1
Z
2 ⌦1
ruj+1 2
+ ⌫2
⇣ ⌘ 2 2 |[un ]| un+1 + un+1 ds 1 2 I Z ⇣ ⌘ 2 2 2 + t [u0 ] u11 + u12 ds
+ t u1
t
n ✓ X
2 ⌦2
◆
I
n ✓ X t + f1 (tj+1 ) ⌫ 1 j=1
2 H
1 (⌦
1)
+
t
⌫2
f2 (t
j+1
)
2 H
1 (⌦
2)
◆
. (3.3)
Remark 3.1 (Discrete energy). The discrete energy equation (3.2) includes the following exact “kinetic energy” of the discretization, denoted by KE n+1 : Z ⇣ ⌘ 2 2 2 KE n+1 = un+1 + t |[un ]| un+1 + un+1 ds. 1 2 I
The corresponding “energy dissipation” of the method, denoted by ✏n+1 is ✏n+1 = 2 t
n ✓ X j=1
⌫1 ruj+1 1
+
t
2 ⌦1
n Z X j=1
I
2
+ ⌫2 ruj+1 2
⌦2
|[uj ]|1/2 uj+1 1 +
t
|[uj
n Z X j=1
8
◆
I
1 1/2 j ]| u2
|[uj ]|1/2 uj+1 2
2
ds |[uj
1 1/2 j ]| u1
2
ds.
Proof. . Let j 2 {1, 2, . . . , N v1 = uj+1 1 , (2.8) yields
uj+1 1
uj1 t
, u1j+1
!
+ ⌫1 ruj+1 1
⌦1
1} and consider (2.8) with n = j. Setting
Z
I
2 ⌦1
+
Z
2 j |uj+1 1 | |[u ]| ds
I
j 1/2 uj2 · uj+1 |[uj 1 |[u ]|
]|
1 1/2
⇣ ⌘ ds = f1 (tj+1 ), uj+1 1
⌦1
. (3.4)
Similarly, choosing v2 = uj+1 in (2.9) yields 2
uj+1 2
uj2 t
, u2j+1
!
ruj+1 2
+ ⌫2
⌦2
Z
I
2 ⌦2
+
Z
2 j |uj+1 2 | |[u ]| ds
I
j 1/2 uj1 · uj+1 |[uj 2 |[u ]|
]|
1 1/2
⇣ ⌘ ds = f2 (tj+1 ), uj+1 2
⌦2
. (3.5)
Applying the identity (3.1) to the first terms on the left hand side of (3.4)-(3.5) and adding, it follows that
1 2 t
✓
+ ⌫1
2
uj+1 1
⌦1
ruj+1 1
Z
I
2 ⌦1
uj1 + ⌫2
uj2 · uj+1 [uj ] 1
2 ⌦1
+
uj+1 1
ruj+1 2 1/2
[uj
2 ⌦2
uj1
+
1 1/2
Z
I
2 ⌦1
uj+1 2
+
uj+1 1
2
2 ⌦2
2
uj2
[uj ] ds +
Z
⌦2
I
Z
+
uj+1 2
uj+1 2 2
1/2
uj1 · uj+1 [uj ] [uj 1 ] 2 ⇣ ⌘ ⇣ ⌘ = f1 (tj+1 ), uj+1 + f2 (tj+1 ), uj+1 1 2 ]
ds
2 ⌦2
[uj ] ds 1/2
ds
I
⌦1
The following two interface integrals may be expressed in a di↵erent way: Z 2 1/2 1/2 j uj+1 [u ] ds uj2 · uj+1 [uj ] [uj 1 ] ds 1 1 I I Z ⇣ ⌘ j 1/2 j 1/2 = uj+1 · uj+1 uj2 |[uj 1 ]|1/2 ds 1 |[u ]| 1 |[u ]| ZI Z 2 2 j = uj+1 [u ] ds uj2 [uj 1 ] ds 1 2 I 2 I Z 2 j 1/2 + uj+1 uj2 |[uj 1 ]|1/2 ds. 1 |[u ]| 2 I
uj2
Z
9
⌦2
. (3.6)
◆
This expression is inserted into (3.6) and the remaining two interface integrals are treated analogously:
◆ 2 2 2 2 j j+1 j uj+1 u + u u 1 1 2 2 2 t ⌦1 ⌦1 ⌦2 ⌦2 ◆ ✓ 2 2 2 1 + uj+1 uj1 + uj+1 uj2 + ⌫1 ruj+1 + ⌫2 ruj+1 1 2 1 2 2 t ⌦1 ⌦2 ⌦1 ⇢ Z ⇢ Z 2 2 2 2 [uj ] ds [uj 1 ] ds + uj+1 + uj+1 uj1 + uj2 1 2 2 I 2 I Z 2 j 1/2 + uj+1 uj2 |[uj 1 ]|1/2 ds 1 |[u ]| 2 I Z 2 j 1/2 + uj+1 uj1 |[uj 1 ]|1/2 ds 2 |[u ]| 2 ⇣ I ⌘ ⇣ ⌘ j+1 j+1 = f1 (tj+1 ), uj+1 + f (t ), u . 2 1 2 1
✓
⌦1
2 ⌦2
(3.7)
⌦2
Multiply through (3.7) by 2 t, and summing over j = 1, 2, . . . , n proves the energy equation (3.2). To prove (3.3), bound the right hand side of (3.7) by:
⇣
f1 (tj+1 ), uj+1 1
⌘
⌦1
⇣ ⌘ + f2 (tj+1 ), uj+1 2
1 f1 (tj+1 ) 2⌫1
2 H
⌦2
1 (⌦
1)
1 2 f2 (tj+1 ) H 1 (⌦ ) 2 2⌫2 2 ⌫1 ⌫2 + ruj+1 + ruj+1 1 2 2 2 ⌦1 +
2 ⌦2
. (3.8)
Insert (3.8) into (3.7), and the desired result follows by multiplying through by 2 t and summing over j = 1, 2, . . . , n. Theorem 3.2 (Convergence of Algorithm 2.1). Let uk1 2 X1,h satisfy (2.8) for each k 2 {2, · · · , n N 1}, and uk2 2 X2,h satisfy (2.9) for each k 2 {2, · · · , n N 1}. Let ⌫˜ = max{⌫1 1 , ⌫2 1 }, ⌫ˆ = max{⌫1 , ⌫2 }, and let ⇣ Dn+1 = ⌫˜3 1 + 4 E n+1 + ru(tn+1 )
4
⌘
,
where E n+1 = maxj=0,1,...,n+1 {max{ku(tj )k4I , kuj k4I }}. Assume t 1/Dn+1 , and that (u, p) is a strong solution of the coupled NSE system (1.1)-(1.6) with ut 2 L2 (0, T ; X) and utt 2 L2 (0, T ; L2 (⌦)). Then the solution un+1 of Algorithm 2.1 10
satisfies: u(tn+1 )
un+1
t + 2
⌫1
+ + +
t
j X
n=1
t ⌫1 +
n X
k=1
C exp +
2
✓
1
r(u1 (t
k+1
Dn+1 t Dn+1
r(u1 (t1 )
)
!(
u11 )
v11 2V1,h
inf
v21 2V2,h
r(u2 (t ) 1
2 ⌦1
2 ⌦2
u(t1 )
+
+
2 v21 ) ⌦ 2
2
+ 1inf
v 2Vh
k+1
2 uk+1 ) ⌦ 2 2
)
u(t1 )
1 2 r(u2 (t0 ) u02 ) ⌦ 2 2 1 + inf r(u2 (t0 ) 2 v20 2V2,h
v2Vh
2
2
u1
r(u2 (t
1 2 r(u1 (t0 ) u01 ) ⌦ 1 2 1 + inf r(u1 (t0 ) 2 v10 2V1,h
t2 kutt kL2 (0,T ;L2 (⌦)) + inf k(u + 2 t2 kut kL2 (0,T ;X) + T
+ ⌫2
n X
k=1
2 v11 ) ⌦ 1
inf r(u1 (t ) ✓ t ⌫2 r(u2 (t1 ) u12 ) 1
2 uk+1 ) ⌦ 1 1
v1
2 v10 ) ⌦ 1
◆
2 v20 ) ⌦ 2
◆
2
v)t kL2 (0,T ;L2 (⌦)) + inf kp
max
k=2,...,n+1
✓
q2Qh
inf
vk 2Vh
r(u(tk )
vk )
2
2
2
qkL2 (0,T ;L2 (⌦)) ◆) ,
(3.9)
where C has the following dependence on , ⌫1 , and ⌫2 : C = O max{ˆ ⌫ , ⌫˜, (1 +
!
t4 )˜ ⌫ 3 , 2 } .
Remark 3.2. The timestep restriction is for the error estimate and not stability. It comes from the discrete Gronwall inequality and is typical for nonlinearly implicit methods. Our stability and error analysis can readily be extended to the linearly implicit method (where un+1 · run+1 is replaced by uni · run+1 ). See Ingram [13] for a i i i detailed analysis of this e↵ect. Proof. To prove convergence we will consider the error equations on each subdomain independently first. The resulting error bounds from the two subdomains will then be summed. Begin by taking the variational formulation for the true solution u1 (tn+1 ) on ⌦1 and subtracting (2.8), 1 (u1 (tn+1 ) un+1 ) (u1 (tn ) un1 ), v1 ⌦ + ⌫1 r(u1 (tn+1 ) un+1 ), rv1 1 1 1 t n+1 n+1 n+1 n+1 n+1 + c1 u1 (t ); u1 (t ), v1 c1 u1 ; u1 , v1 p1 (t ), r · v1 ⌦ 1 Z Z n+1 n+1 n+1 n + u1 (t )|[u(t )]| · v1 ds u1 |[u ]| · v1 ds I ZI Z n n 1/2 n 1 1/2 + u2 |[u ]| |[u ]| · v1 ds u2 (tn+1 )|[u(tn+1 )]| · v1 ds I ✓ ◆I u1 (tn+1 ) u1 (tn ) @u1 (tn+1 ) , v1 , 8v1 2 V1,h . = t @t ⌦1
Errors are decomposed using the equation ui (tj )
uji = (ui (tj )
u ˜ji ) + (˜ uji 11
uji ) = ⌘ij +
j i
⌦1
where ⌘ij 2 Vi and ji 2 Vh,i with u ˜ji 2 Vh,i arbitrarily chosen. We also subtract an j arbitrarily chosen qi 2 Qh,i from the pressure for free. Thus choosing v1 = n+1 the 1 above error equation is rewritten as 1 ⌘ n+1 ⌘1n , n+1 1 t 1 n+1 2 + ⌫1 kr 1 k⌦1
⌦1
+
1 t
n+1 1
p1 (tn+1 )
n 1,
n+1 1 ⌦1
+ ⌫1 r⌘1n+1 , r
q1n+1 , r ·
n+1 1 ⌦1 n+1 n+1 c1 u1 ; u1 ,
n+1 + c1 u1 (tn+1 ); u1 (tn+1 ), n+1 1 1 Z Z + u1 (tn+1 )|[u(tn+1 )]| · n+1 ds un+1 |[un ]| · n+1 ds 1 1 1 I I Z Z + un2 |[un ]|1/2 |[un 1 ]|1/2 · n+1 ds u2 (tn+1 )|[u(tn+1 )]| · 1
=
I R1n+1 ,
I
n+1 1 ⌦1
n+1 1 ⌦1
n+1 1
ds
,
u (tn+1 ) u (tn )
@u (tn+1 )
j j , for j = 1, 2. The goal is to bound n+1 , where Rjn+1 = j 1 t @t hence we move many terms to the right hand side. Apply (3.1) on the left hand side:
1 n+1 2 n 2 k n+1 k2 k n1 k2⌦1 + k n+1 k⌦1 1 k⌦1 + ⌫1 kr 1 1 2 t ⇢1 Z ⌦1 Z + u1 (tn+1 )|[u(tn+1 )]| · n+1 ds un+1 |[un ]| · n+1 ds 1 1 1 I I ⇢ Z Z + un2 |[un ]|1/2 |[un 1 ]|1/2 · n+1 ds u2 (tn+1 )|[u(tn+1 )]| · 1 I
1 ⌘ n+1 ⌘1n , n+1 1 ⌦1 t 1 n+1 n+1 n+1 + p1 (t ) q1 , r · 1 =
c1 u1 (t
n+1
); u1 (t
n+1
),
n+1 1
I
⌫1 r⌘1n+1 , r
n+1 1
ds
n+1 1 ⌦1
⌦1
+ c1 un+1 ; un+1 , 1 1
n+1 1
+ R1n+1 ,
n+1 1 ⌦1
. (3.10) Bounds for the non-interface terms follow as in the standard NSE case, (see e.g. [15]). The interface terms require special treatment. First, note that by adding and subtracting terms in sequence one may derive the following expression: Z Z u1 (tn+1 )|[u(tn+1 )]| · n+1 ds un+1 |[un ]| · n+1 ds 1 1 1 I I Z = u1 (tn+1 ) |[u(tn+1 )]| |[u(tn )]| · n+1 ds 1 ZI + u1 (tn+1 ) (|[u(tn )]| |[˜ un ]|) · n+1 ds 1 I Z + u1 (tn+1 ) (|[˜ un ]| |[un ]|) · n+1 ds 1 I Z Z + ⌘1n+1 |[un ]| · n+1 ds + | n+1 |2 |[un ]| ds. 1 1 I
I
(3.11) The splitting used for the remaining pair of interface integrals is more complicated due to the product |[un ]|1/2 |[un 1 ]|1/2 . The key is to recognize an error of O( t) is committed in replacing the product |[u(tn )]|1/2 |[u(tn 1 )]|1/2 with the average 21 (|[u(tn )]|+ 12
|[u(tn
1
)]|). We proceed as follows.
Z un2 |[un ]|1/2 |[un 1 ]|1/2 · n+1 ds u2 (tn+1 )|[u(tn+1 )]| · n+1 ds 1 1 I I ⇢ Z 1 (|[un ]| + |[un 1 ]|) · n+1 ds = un2 |[un ]|1/2 |[un 1 ]|1/2 1 2 I Z 1 + un2 (|[un ]| + |[un 1 ]|) (|[˜ un ]| + |[˜ un 1 ]|) · n+1 ds 1 2 I Z 1 + un2 (|[˜ un ]| + |[˜ un 1 ]|) (|[u(tn )]| + |[u(tn 1 )]|) · n+1 ds 1 2 I ⇢ Z 1 + un2 (|[u(tn )]| + |[u(tn 1 )]|) |[u(tn+1 )]| · n+1 ds 1 2 I Z Z n+1 n n+1 |[u(t )]| · ds ⌘2n |[u(tn+1 )]| · n+1 ds 2 1 1 I I Z + (u2 (tn ) u2 (tn+1 ))|[u(tn+1 )]| · n+1 ds 1
Z
(3.12)
I
Using (3.11)–(3.12), substitution for the interface integrals in (3.10) is performed. The last term of (3.11) is positive and may be kept on the left hand side of (3.10), but all other interface terms are now moved to the right hand side, and bounded. Applying the reverse triangle inequality and the equality [a] [b] = [a b] for a, b 2 X, we first note |[u(tj )]|
|[˜ uj ]| |[⌘ j ]|
and
|[˜ uj ]|
j
|[uj ]| |[
]|.
Thus (3.10)–(3.12) yield the following error inequality: 1 2 t
k
n+1 2 k⌦ 1 1
k
n 2 1 k⌦ 1
1 ⌘ n+1 ⌘1n , n+1 1 t 1 n+1 n+1 c1 u1 (t ); u1 (t ),
+k
n+1 1
n 2 1 k⌦1
+ ⌫1 kr
⌦1
⌫1 r⌘1n+1 , r
n+1 1
+ c1 un+1 ; un+1 , 1 1
n+1 2 k⌦ 1 1
+
Z
+ p1 (tn+1 )
n+1 1 ⌦1
n+1 1
+ R1n+1 ,
I
|
n+1 2 | |[un ]| ds 1
q1n+1 , r ·
n+1 1 ⌦1
n+1 1 ⌦1
+ I1 + I2 + I3 + I4 ,
(3.13) where the Ij -terms denote the following interface integrals. Z Z n+1 n+1 n+1 n I1 = |u1 (t )||[u(t ) u(t )]|| 1 | ds + |u1 (tn+1 )||[⌘ n ]|| n+1 | ds 1 I Z ZI + |u2 (tn ) u2 (tn+1 )||[u(tn+1 )]|| n+1 | ds + |⌘1n+1 ||[un ]|| n+1 | ds 1 1 I I Z Z 1 + |un2 |(|[⌘ n ]| + |[⌘ n 1 ]|)| n+1 | ds + |⌘2n ||[u(tn+1 )]|| n+1 | ds 1 1 2 I I Z 1 |un2 |(|[ n ]| + |[ n 1 ]|)| 2 I Z + | n2 ||[u(tn+1 )]|| n+1 | ds 1
I2 =
n+1 | ds 1
I
13
+
Z
I
|u1 (tn+1 )||[
n
]||
n+1 | ds 1
Z
I3 =
I
I4 =
Z
I
|un2 |
1 (|[u(tn )]| + |[u(tn 2
|un2 | |[un ]|1/2 |[un
]|
1 1/2
1
)]|)
|[u(tn+1 )]| |
1 (|[un ]| + |[un 2
1
n+1 | ds 1
]|) |
n+1 | ds 1
The application of H¨ older’s inequality, Young’s inequality, and (2.2) to the noninterface terms of (3.13) gives the estimate Z 1 n+1 n+1 2 2 n 2 n 2 k n+1 k k k + k k + ⌫ kr k + | n+1 |2 |[un ]| ds 1 ⌦1 1 ⌦1 1 ⌦1 ⌦1 1 1 1 1 2 t I ⌫1 1 d 2 4 n+1 2 r⌘1n+1 ⌦ + 3 3 ru1 (tn+1 ) ⌦ + kp1 (tn+1 ) q1n+1 k2L2 (⌦1 ) 1 ⌦1 1 1 2"2 4"3 ⌫1 2"2 ⌫1 ✓ ◆ ✓ n+1 2 3"3 C ⌘1 ⌘1n 2 2 + 2"1 + "2 + ⌫1 r n+1 + + R1n+1 V 0 1 ⌦1 1 4 2"1 ⌫1 t 0 V1 ◆ 2 run+1 r⌘1n+1 ⌦1 + kru1 (tn+1 )k2⌦1 kr⌘1n+1 k2⌦1 + un+1 1 1 ⌦1 ⌦1 + I1 + I2 + I3 + I4 ,
(3.14) for constants C and "i , i = 1, 2, 3 independent of ⌫1 . The interface integrals are bounded using Lemma 2.2, with particular attention paid to the constants involving 2 n+1 2 , ⌫1 , ⌫2 and terms with r n+1 or . In the term I1 , each piece is 1 1 ⌦1 ⌦1 bounded using (2.5) with ↵1 = 192. Then the first two terms of I2 are bounded using (2.7) with ↵1 = 16, ↵2 = 32, 1 = 96 and the last term using (2.6) with ↵2 = 8 and 1 = 48. Some basic inequalities can be used to simplify I3 , I4 and derive corresponding error estimates. Applying the reverse triangle inequality for I3 yields 1 (|[u(tn )]| + |[u(tn 1 )]|) |[u(tn+1 )]| 2 1 1 |[u(tn )]| |[u(tn+1 )] + |[u(tn 1 )]| |[u(tn+1 )]| 2 2 1 1 n n+1 [u(t ) u(t )] + [u(tn 1 ) u(tn+1 )] . 2 2 I4 is bounded by first noting |[un ]|1/2 |[un
]|
1 1/2
Add and subtract u(tn ) |[un ]|1/2 |[un
]|
1 1/2
(3.15)
2 1 1 (|[un ]| + |[un 1 ]|) = |[un ]|1/2 |[un 1 ]|1/2 2 2 1 |[un ]|1/2 |[un 1 ]|1/2 |[un ]|1/2 + |[un 1 ]|1/2 2 1 1 = |[un ]| |[un 1 ]| [un un 1 ] . 2 2
u(tn
1
),
1 (|[un ]| + |[un 2
1
]|)
1 1 [(un u(tn 1 )) + (u(tn 1 ) un 1 )] + [u(tn ) u(tn 1 )] 2 2 o 1n n |[ ]| + [ n 1 ] + |[⌘ n ]| + [⌘ n 1 ] + [u(tn ) u(tn 1 )] . 2
14
(3.16)
The results (3.15)–(3.16) are then bounded using the same approach as particular terms in I1 and I2 . For example, I3 can be bounded in the same manner as the first term from I1 , while the first term from I4 is bounded the same way as the first term of I2 . Thus, the interface terms are bounded by I1 + I2 + I3 + I4 ⇣ ⇣ ⌫1 2 4 3 4 r n+1 + C⌫ 1 + 48 kun2 kI + 12 [u(tn+1 )] 1 1 ⌦1 16 ⇣ ⌘ ⌫ ⇣ ⌫1 2 2 2 2 + kr n1 k⌦1 + r n1 1 ⌦ + kr n2 k⌦2 + r n2 1 1 16 16 n ⇣ ⌘ ⇣ n 2 1 k⌦ 1
+ C4 M n+1 ⌫1 3 k
n 1 2 1 ⌦1
+
+ ⌫2 3 k
n 2 2 k⌦ 2
+
4 I
⌘⌘
2 ⌦2
n 2
n+1 2 1 ⌦1
⌘
+ C2 Ln+1 P n+1 ⌘o 1 2 , ⌦ 2
(3.17) where C is a constant independent of , ⌫1 , and ⌫2 . The quantities Ln+1 , M n+1 and P n+1 are defined by: Ln+1 = u1 (tn+1 ) M n+1 = u1 (tn+1 )
2 I 4 I
and P n+1 = u1 (tn+1 )
+ u2 (tn+1 ) + u2 (tn+1 ) 2 I
u1 (tn )
+ u1 (tn+1 )
2
2
4
4
+ kun1 kI + kun2 kI ,
+ kun1 kI + kun2 kI ,
+ u2 (tn+1 )
u1 (tn
2
2 I 4 I
1
2 I
)
2
+ k⌘1n kI + k⌘2n kI + ⌘1n
u2 (tn )
2 I
+ u2 (tn+1 )
1 2 I
1 2 I
+ ⌘2n
u2 (tn
1
+ ⌘1n+1
)
2 I
2 I
+ ⌘2n+1
2 I
.
(The last term in the definition of P n+1 is auxiliary, it is included in the definition for use in the estimate on ⌦2 .) Then, with choices for "1 , "2 , "3 satisfying 2"1 +"2 +3"3 /4 = 1/16, (3.14) can be bounded as follows: 1 2 t
n+1 2 k⌦1 1
k
k
n 2 1 k⌦ 1
+k
n+1 1
7⌫1 kr 8
+
n 2 1 k⌦ 1
n ⇣ 4 C⌫1 3 1 + kru1 (tn+1 )k4⌦1 + 4 kun2 kI + [u(tn+1 )]
+ C⌫1
1
2
⌘1n
+
t
V10
kR1n+1 k2V10
+ kp1 (t
n+1
)
4 I
⌘o
n+1 2 1 ⌦1
q1n+1 k2⌦1
!
⌘ ⌫ ⇣ ⌘ 2 2 2 2 + r n1 1 ⌦1 + kr n2 k⌦2 + r 2n 1 ⌦2 + C2 Ln+1 P n+1 16 ⌘ n ⇣ ⇣ ⌘o 3 n 1 2 3 n 1 2 4 n+1 n 2 n 2 + C M ⌫1 k 1 k⌦ 1 + 1 + ⌫ k k + 2 ⌦2 2 2 ⌦1 ⌦2 ⇣ ⌘ 1 1 n+1 n+1 n+1 2 n+1 2 + C ⌫1 + ⌫1 kru1 (t )k⌦1 + ⌫1 u1 ru1 kr⌘1 k⌦1 . ⌦1 ⌦1 +
⌫1 ⇣ kr 16
⌘1n+1
n+1 2 k⌦1 1
n 2 1 k⌦ 1
(3.18)
Let Sn+1 denote Sn+1 = u1
2
+ t
Z
I
[u0 ]
⇣
u11
2
+ u12
2
⌘
ds
n ✓ X t + f1 (tj+1 ) ⌫ 1 j=1
15
2 V10
+
t ⌫2
f2 (t
j+1
)
2 V20
◆
, (3.19)
i.e., the right-hand side of the stability estimate (3.3). Also, let ⌫ˆ := max{⌫1 , ⌫2 } and ⌫˜ := max{⌫1 1 , ⌫2 1 }. Then (3.18) can be simplified to 1 7⌫1 n 2 k 1n+1 k2⌦1 k n1 k2⌦1 + k n+1 kr n+1 k2⌦1 1 k⌦ 1 + 1 1 2 t 8 n ⇣ ⌘o 4 4 C⌫1 3 1 + kru1 (tn+1 )k4⌦1 + 4 kun2 kI + [u(tn+1 )] I + C⌫1
⌘1n+1
1
2
⌘1n
+
t
⌫1 ⇣ + kr 16
V10
kR1n+1 k2V10
)
q1n+1 k2⌦1
!
(3.20)
⌘ ⌫2 ⇣ n 1 2 n 2 n 2 k + r + kr k + r 1 ⌦1 2 ⌦2 2 ⌦2 ⇣ 16 ⌘ n 2 n 1 2 2 n+1 n+1 4 n+1 3 + C L P + C M ⌫˜ k k + ⇣ ⌘ 1/2 + C ⌫1 + ⌫1 1 kru1 (tn+1 )k2⌦1 + ⌫1 1 Sn+1 run+1 kr⌘1n+1 k2⌦1 . 1 ⌦1 n 1 2 1 ⌦1
⌘
+ kp1 (t
n+1
n+1 2 1 ⌦1
Proceeding further requires estimates for the solution on ⌦2 . This analysis is done in the same manner as above for ⌦1 using (2.9), and a similar estimate is derived: 1 2 t
n+1 2 k⌦ 2 2
k
n 2 2 k⌦2
k
+k
n+1 2
n 2 2 k⌦ 2
+
7⌫2 kr 8
n ⇣ 4 C⌫2 3 1 + kru2 (tn+1 )k4⌦2 + 4 kun1 kI + [u(tn+1 )]
+ C⌫2
⌘2n+1
1
⌫1 ⇣ + kr 16
n+1 2 k⌦2 2 4 I
2
⌘2n
+
t
V20
kR2n+1 k2V20
)
q2n+1 k2⌦2
n+1 2 2 ⌦2
!
(3.21)
⌘ ⌫2 ⇣ 2 2 + r + kr n2 k⌦2 + r n2 1 ⌦ 2 ⇣ 16 ⌘ n 2 n 1 2 2 n+1 n+1 4 n+1 3 + C L P + C M ⌫˜ k k + ⇣ ⌘ 1/2 + C ⌫2 + ⌫2 1 kru2 (tn+1 )k2⌦2 + ⌫2 1 Sn+1 run+1 kr⌘2n+1 k2⌦2 . 2 ⌦2 n 1 2 1 ⌦1
n 2 1 k⌦ 1
⌘
+ kp2 (t
n+1
⌘o
Combining (3.20) and (3.21) gives the global estimate 1 k 2 t
n+1 2
k
k
n 2
k +k
n+1
n 2
k
+ C ⌫˜
t
n+1 2 k⌦ 1 1
n 2 1 k⌦ 1
+
7⌫2 kr 8
n+1 2 k⌦2 2
2
+ V0
kR1n+1 k2V10
+
+ r
n 1 2 1 ⌦1
⌘
kR2n+1 k2V20
+ kp(t
n+1
⌫2 ⇣ 2 + kr n2 k⌦2 + r 2n 8 ⇣ ⌘ 2 2 + C2 Ln+1 P n+1 + C4 M n+1 ⌫˜3 k n k + n 1 ⇣ ⌘ 1/2 + C ⌫ˆ + ⌫˜kru(tn+1 )k2 + ⌫˜Sn+1 run+1 kr⌘ n+1 k2 . +
⌫1 ⇣ kr 8
⌘n
7⌫1 kr 8
n+1 2
C ⌫˜3 1 + kru(tn+1 )k4 + 4 M n+1 ⌘ n+1
+
16
)
1 2 ⌦2
q
n+1 2
⌘
k
!
(3.22)
Rearranging terms gives 1
⌫1 ⌫2 2 k + r n+1 + r 1 ⌦ 1 2 2⌘ ⌘ ⌫ ⇣ 2 1 2 2 kr n1 k⌦1 + kr n1 k⌦1 r 1n 1 ⌦ 1 8 ⌘ ⌫ ⇣ ⌘ 2 n 1 2 n 2 n 2 kr 2 k⌦2 + kr 2 k⌦2 r 2 ⌦2 8
k n+1 k2 k 2 t ⇣ ⌫1 2 + r n+1 1 ⌦1 4 ⇣ ⌫2 2 + r n+1 2 ⌦2 4
n 2
n+1
k +k
n 2
n+1 2 2 ⌦2
n+1 2
C An+1 + C B n+1 + C ⌫˜3 1 + kru(tn+1 )k4 + 4 M n+1 ⇣ ⌘ 2 2 + C2 Ln+1 P n+1 + C4 M n+1 ⌫˜3 k n k + n 1
(3.23)
where An+1 and B n+1 are defined by
⌘ n+1
An+1 = ⌫˜
⌘n t
2
+ kR1n+1 k2V10 + kR2n+1 k2V20 + kp(tn+1 )
V0
qn+1 k2
⇣ ⌘ 1/2 B n+1 = ⌫ˆ + ⌫˜kru(tn+1 )k2 + ⌫˜Sn+1 run+1 kr⌘ n+1 k2 .
!
(3.24) Summing over n = 1, 2, . . . , p 1, multiplying through by 2 t, and rearranging gives a global error estimate at time tp :
k
p 2
k +
+
p 1 X
n 2
n+1
+
⌫1
t ⌫1 4
✓
n=1
2 kr
p 2 1 k⌦ 1
+ r
1 + C t ⌫˜ (M + M ) ✓ t ⌫1 1 2 + r 11 ⌦ + r 1 2 2 +C
t
n=1
t
p 1 X
p 1 X
3 4
2
p 1 1
2 ⌦1
1 2
3
0 1
◆
+
r
n+1 2 1 ⌦1
t ⌫2 4
✓
2 kr
+ r
1 2 2 ⌦2 p 1 X
1 + r 2
⌦2
◆
(3.25)
where Dn+1 is defined by 4
◆
n+1 2
Dn+1
n=1
⇣ Dn+1 = ⌫˜3 1 + 4 E n+1 + ru(tn+1 )
0 2 2 ⌦2
!
2
p 1 2
0 2
2
An+1 + B n+1 + 2 Ln+1 P n+1 + C t
n=1
p 2 2 k⌦ 2
n+1 2 2 ⌦2
r
n=1
+ C t ⌫˜ M ✓ ◆ t ⌫2 2 + r ⌦1 2 3 4
+ ⌫2
p 1 X
⌘
(3.26)
,
with
E
n+1
8 n+3 > + M n+2 + M n+1 <M = M n+2 + M n+1 > : n+1 M 17
if 1 n p if n = p 2 if n = p 1
3 .
(3.27)
To use the discrete Gronwall’s inequality (see [15]), note that for j = 1, 2, . . . , N , write (3.25) as aj+1 +
j X
bn+1 +
t
n=1
j X
cn+1
n=1
J +C t
j X
Dn+1 an+1 + C
t
n=1
j X
An+1 + B n+1 + 2 Ln+1 P n+1
(3.28)
n=1
where aj+1 =
j+1 2
j+1
bj+1 =
,
j 2
cj+1 = ⌫1 r
,
j+1 1
2 ⌦1
+ ⌫2 r
2
j+1 2
⌦2
,
and J = 1 + C t ⌫˜3 4 (M 2 + M 3 ) ✓ t ⌫1 1 2 + r 11 ⌦ + r 1 2 2
1 2 0 1
+ C t ⌫˜3 4 M 2 ✓ ◆ t ⌫2 2 + r 12 ⌦1 2
0 2 2 ⌦2
1 + r 2
0 2 2 ⌦2
◆
.
Applying Gronwall’s inequality, (3.28) gives j+1
a
+
j X
b
n+1
+
n=1
exp C t
t
j X
cn+1
n=1 j X
n=1
n+1
!(
J +C t
j X
A
n+1
+B
n+1
+ L
n=1
2
n+1
P
n+1
)
(3.29)
where n+1
=
1
Dn+1 . t Dn+1
(3.30)
It remains to bound the individual terms on the right-hand side of (3.29). The terms in J are bounded using the triangle inequality. The terms in An+1 are bounded in the usual way (see the proof of Theorem 4.1 of [9], for example). In B n+1 , the summed factor run+1 is bounded using the stability estimate (3.3) and the ru(tn+1 ) term is bounded by a standard a priori estimate on the solution to the continuous problem. The terms in Ln+1 and M n+1 are bounded by a standard a priori estimate on u(tn+1 ) or using S n for the discrete solutions, (the bound (2.4), and the stability estimate (3.3)). The terms in P n+1 are bounded the same way. The estimate (3.9) is derived through the application of the triangle inequality to the left hand side of (3.29). Corollary 3.1 (Convergence rate for MINI-element,[1]). Let (u, p) be a solution of (1.1)-(1.6) satisfying the assumptions of Lemma 3.1 and Theorem 3.2, approximated using the MINI-element on a mesh of maximum element width h by Algorithm 2.1. Assume the velocity data u0 , u1 satisfies ku(t0 )
u0 kX + ku(t1 ) 18
u1 kX C1 h,
for a generic constant C1 independent of t, h. Then there exists C > 0 independent of h, t such that the discrete velocity satisfies the optimal error estimate t
⌫1
n X
k=1 2
r(u1 (tk+1 )
C (h +
uk+1 ) 1
2 ⌦1
+ ⌫2
n X
k=1
t2 ).
r(u2 (tk+1 )
uk+1 ) 2
2 ⌦2
!
Proof. This follows by using standard interpolation error estimates for the MINIelement, e.g. [1, 5, 15], in the right hand side of the error bound (3.9). 4. Numerical Experiments. In this section numerical experiments are presented that support the theory described in this report. An example is provided to verify the asymptotic properties predicted by Corollary 3.1 using a manufactured solution. A subsequent example in Section 4.2 demonstrates a benefit of the unconditional stability property for the proposed algorithm. All computations were performed using the software package FreeFem++, [12]. 4.1. [0, 1] ⇥ [ n ˆ 2 = [0, function
Problem 1: Analytic Solution. Assume ⌦1 = [0, 1] ⇥ [0, 1] and ⌦2 = 1, 0], so I is the portion of the x-axis from 0 to 1. Then n ˆ 1 = [0, 1]T and T 1] . For a, ⌫1 , ⌫2 , and all arbitrary positive constants, the right hand side f1 is chosen to ensure that u1,1 (t, x, y) = ax2 (1
x)2 (1
y)e
u1,2 (t, x, y) = axy( 2 + y + 6x p1 (t, x, y) = e
t
t
,
3xy
4x2 + 2x2 y)e
t
,
cos(⇡x) sin(⇡y) .
One may verify that the above velocity field is divergence-free and that u1,2 is zero on I. The function f2 is chosen to ensure that 1) u2,1 satisfies the interface condition, 2) u2,2 is zero on I, 3) u2 is a divergence-free velocity field, and 4) p2 (t, x, y) = e t cos(⇡x) sin(⇡y). The spatial discretization is accomplished using the conforming P1 +bubble/P1 element pair, (the MINI-element, [1]). A convergence study is performed using the parameter values a = 1, ⌫1 = 0.5, ⌫2 = 0.05, = 100, T = 1. Table 4.1 gives the errors produced Algorithm 2.1. A compact notation is used to denote the various errors. Letting ui (t), pi (t) denote the NSE solution in ⌦i , i = 1, 2, approximated at t = tn by uni , pni respectively, the errors are Err(ui ) =
t
N X
n=0
Err(pi ) =
t
N X
n=0
kui (t )
uni k2H 1 (⌦i )
kpi (t )
pni k2L2 (⌦i )
n
n
!1/2
!1/2
.
The expected first-order convergence rates for the velocity are achieved. The pressure rates may be limited to first-order with the MINI-element; here the pressure convergence rates are slightly above the predicted value but decreasing. 19
t 1/4 1/8 1/16 1/32 1/64
h 1/4 1/8 1/16 1/32 1/64
Err(u1 ) 0.196536 0.071197 0.031561 0.015434 0.007684
rate 1.46 1.17 1.03 1.01
Err(u2 ) 1.673680 0.479848 0.160358 0.069638 0.032969
rate 1.80 1.58 1.20 1.08
Err(p1 ) 0.174935 0.034112 0.008706 0.002815 0.000978
rate 2.36 1.97 1.63 1.53
Err(p2 ) 0.172981 0.033365 0.008108 0.002893 0.001089
rate 2.37 2.04 1.49 1.41
Table 4.1 Errors for the manufactured solution test problem
4.2. Growth of error with decoupled methods. In addition to the GA scheme previously described, the following algorithms will be compared: Algorithm 4.1 (Fully Implicit Scheme). Let t > 0, f 2 H 1 (⌦). For each M 2 N, M Tt , given un 2 Vh , n = 0, 1, 2, · · · , M 1, find un+1 2 Vh satisfying ✓ n+1 ◆ X u un , v + ⌫ run+1 , v + ci (un+1 ; un+1 , vi ) i i t i=1,2 (4.1) Z n n+1 n+1 + |[u ]|[u ][v] = (f (t ), v) . I
Algorithm 4.2 (First-order IMEX Scheme). Let t > 0, f 2 H 1 (⌦). For each M 2 N, M Tt , given un 2 Vh , n = 0, 1, 2, · · · , M 1, find un+1 2 V1,h satisfying 1 ✓
un+1 1
un1 t
, v1
◆
⌦1
+ ⌫1 run+1 , rv1 1 +
Z
I
⌦1
+ c1 (un+1 ; un+1 , v1 ) 1 1
|[un ]|un+1 · v1 ds 1
= f1 (tn+1 ), v1
⌦1
Z
I
|[un ]|un2 · v1 ds
, 8v1 2 V1,h .
and un+1 2 V2,h satisfying 2 ✓ n+1 ◆ u2 un2 , v2 + ⌫2 run+1 , rv2 ⌦ + c2 (un+1 ; un+1 , v2 ) 2 2 2 2 t ⌦2 Z Z + |[un ]|un+1 · v ds |[un ]|un1 · v2 ds 2 2 I
(4.2)
(4.3)
I
= f2 (tn+1 ), v2
⌦2
, 8v2 2 V2,h .
The fully implicit scheme is fully coupled at the interface. A fully coupled calculation of this type is not currently used for a motivating application like climate modelling due to the expense of solving the large system of equations, and is used here only to generate reference data since it does not introduce the decoupling errors associated with the partitioned time stepping methods. The implicit-explicit (IMEX) scheme is, perhaps, a more intuitive method of partitioning the calculations than the GA method. However, Lemma 3.1 predicts energetic stability of the GA method regardless of the choice of problem parameters or discretization parameters. A similar result for the IMEX method is not known, and in fact this method may introduce energy numerically under some conditions. We proceed to investigate this e↵ect in the following example. Let ⌦1 = [0, 10] ⇥ [0, 1] and ⌦2 = [0, 10] ⇥ [ 1, 0], so I is the portion of the x-axis from 0 to 10, n ˆ 1 = [0, 1]T and n ˆ 2 = [0, 1]T . The kinematic viscosities are ⌫1 = 1 20
and ⌫2 = 0.2, with = 100. Consider the following functions: g1,1 (t, x, y) = 1.005 · x(10 g1,2 (t, x, y) = 1.005 · y(1 g2,1 (t, x, y) = x(10
x)(1
2y) exp( 3t) ,
y)(2x
10) exp( 3t) ,
x)(1 + 2y) exp( 2t) ,
g2,2 (t, x, y) = y(1 + y)(2x
10) exp( 2t) .
Dirichlet data is imposed by choosing ui,j = gi,j on @⌦i \I for i, j = 1, 2. A divergencefree initial condition is chosen such that ui · n ˆ i = 0 by setting ui,j (0, x, y) = gi,j (0, x, y) on ⌦i , for i, j = 1, 2. A mesh is constructed as in the previous section, consisting of triangles of horizontal width 10/32 and vertical width 1/32. The problem is run with a time step size of 1/40 on the time interval [0, 1]. The right hand side functions f1 and f2 are set to zero. Then Algorithm 4.1 is run. The implicit solution at the first time step and the initial data are used to start the two partitioned algorithms at the second time step. The partitioned algorithms are compared at time steps starting at t2 = 2 t. The test is constructed to demonstrate conditions under which the IMEX method will artificially introduce energy into the calculations due to the numerical treatment of the coupling terms. The fully coupled, first order implicit method, while not a practical scheme for the application, is useful as a point of comparison in the absence of an exact solution. In the first test we compute the kinetic energy of the velocity at the interface for each time step, as predicted by the 3 methods. That is, we plot 0.5 · kunj k2|L2 (I) , for j = 1, 2; there are two values as interface velocities have di↵erent values on ⌦1 and ⌦2 . These values are plotted in Figure 4.1 below. 1000 900 800
1000 Coupled IMEX method GA method
900 800 700
0.5 ku2 k2L2 (I)
0.5 ku1 k2L2 (I)
700 600 500 400
600 500 400
300
300
200
200
100
100
0 0
Coupled IMEX method GA method
0.5 time
1
0 0
0.5 time
1
Fig. 4.1. Relative solution sizes on I. Left: data on ⌦1 . Right: data on ⌦2 .
The damping terms exp( 2t) and exp( 3t) were included in the boundary conditions so that the true solution would decay to zero rapidly. Thus any growth in an approximate solution is an instability. (If these multipliers are not included we observed a slow, long term growth in the energy of the coupled solution.) The kinetic energy in the implicit methods decayed to zero rapidly. The energy in the GA scheme exhibited oscillations, decaying to zero over time. The energy in the IMEX scheme exhibited a similar behavior, but with a significantly higher kinetic energy. It is natural to investigate the growth of relative error for the partitioned methods, in light of the excess energy they exhibit. To distinguish the fully coupled data from 21
the partitioned methods, the solution to (4.1) on ⌦1 is denoted by un1,imp . Similarly, the solution to (4.1) on ⌦2 is denoted by un2,imp . Define the relative errors for the partitioned algorithms as RelErr(unj ) =
kunj unj,imp kH 1 (⌦j ) , for j = 1, 2. kunj,imp kH 1 (⌦j )
These errors are plotted in Figure 4.2. The relative errors for the IMEX method grow faster in both subdomains than those of the GA method. 8 7
3 IMEX method GA method
2.5
IMEX method GA method
6 2
RelErr(un2 )
RelErr(un1 )
5 4
1.5
3
1
2 0.5
1 0 0
0.5 time
1
0 0
0.5 time
1
Fig. 4.2. Relative solution errors. Left: errors on ⌦1 . Right: errors on ⌦2 .
5. Conclusions. For nonlinear interface conditions, stability of partitioned methods depends critically on their precise treatment. We have shown that, perhaps surprisingly, geometric averaging yields an unconditionally stable partitioned method while the most natural treatment may excessively introduce artificial energy. The kinetic energy and energy dissipation predicted by the model has been proven to deviate from values obtained using a fully implicit discretization by an amount proportional (asymptotically) to the time step size. Convergence as time and spatial discretization parameters go to zero was proven. REFERENCES [1] D. N. Arnold, F. Brezzi and M. Fortin, A stable finite element for the Stokes equations, Calcolo, Vol. 21, No. 4, 1984, pp. 337-344. [2] R.A. Adams and J.J.F. Fournier, Sobolev Spaces, vol. 140 of Pure and Applied Mathematics, Academic Press, 2003. [3] U. M. Ascher, S. J. Ruuth and B. T. R. Wetton, Implicit-explicit methods for time-dependent partial di↵erential equations, SIAM J. Num. Anal. 32(3), 1995. [4] C. Bernardi, T. Chac´ on-Rebollo, R. Lewandowski, F. Murat, A model of two coupled turbulent fluids, Part II: Numerical approximations by spectral discretization, SIAM Jour. Num. Analysis, Vol. 40, No. 6, pp. 2368-2394, 2003. [5] S.C. Brenner and L.R. Scott, The Mathematical Theory of Finite Element Methods, SpringerVerlag, 2002. [6] D. Bresch and J. Koko, Operator-Splitting and Lagrange Multiplier Domain Decomposition Methods for Numerical Simulation of Two Coupled Navier-Stokes Fluids, Int. J. Appl. Math. Comput. Sci., Vol. 16, No. 4, 2006, pp. 419-429. [7] E. Burman, Miguel A. Fern´ andez, Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility, INRIA Research Rept. RR-6455, Feb. 2008. 22
[8] P. Causin, J.-F. Gerbeau and F. Nobile, Added-mass e↵ect in the design of partitioned algorithms for fluid-structure problems, INRIA Rept. 5084, 2004. [9] J. Connors, J. Howell and W. Layton, Partitioned timestepping for a parabolic two domain problem, to appear, SIAM Jour. Num. Analysis, 2009. [10] G.T. Cooper and A. Sayfy, Additive Runge-Kutta methods for sti↵ ordinary di↵erential equations, Mathematics of Computation, Vol. 40, No. 161, pp. 207-218. [11] G. P. Galdi, An Introduction to the Mathematical Theory of the Navier-Stokes Equations, Springer Tracts in Natural Philosophy, Volume I, Springer-Verlag, New York, 1994. [12] F. Hecht, A. LeHyaric and O. Pironneau. Freefem++ version 2.24-1, 2008. http://www.freefem.org/ff++. [13] R. Ingram, Unconditional convergence of high-order extrapolations of the Crank-Nicolson, nite element method for the Navier-Stokes equations, University of Pittsburgh technical report TR-MATH-10-09, 2010. [14] C.A. Kennedy and M. H. Carpenter, Additive Runge-Kutta schemes for convection-di↵usionreaction equations, Applied Numerical Math., Vol. 44, pp. 139-181, 2003. [15] W. Layton, Introduction to the Numerical Analysis of Incompressible Viscous Flows, SIAM Computational Science and Engineering 6, 2008. [16] R. Lewandowski, Analyse Math´ ematique et Oc´ eanographic, collection RMA, Masson, 1997. [17] J. -L. Lions, R. Temam and S. Wang, Models for the coupled atmosphere and ocean (CAO I), Computational Mechanics Advances 1, 1993, pp. 1-54. [18] J. -L. Lions, R. Temam and S. Wang, Numerical analysis of the coupled atmosphere-ocean models (CAO II), Computational Mechanics Advances 1, 1993, pp. 55-119. [19] J. -L. Lions, R. Temam and S. Wang, Mathematical theory for the coupled atmosphere-ocean models (CAO III), J. Math. Pures. Appl. (74), 1995, pp. 105-163. [20] P. M¨ uller, The equations of oceanic motions, Cambridge University Press, Cambridge, UK, 2006. [21] W. M. Washington and C. L. Parkinson, An introduction to three-dimensional climate modeling, second edition, University Science Books, California, 2005.
23
1. Introduction. The dynamic core in atmosphere-ocean interaction is a critical component of climatology models and has attracted the interest of many mathematical researchers due to the richness of the theory and technical complexities of efficiently computing approximations to the coupled system using only (uncoupled) atmosphere and ocean solves, (see e.g. [4, 6, 17, 18, 19]). Motivated by this problem we consider uncoupled or “partitioned” methods for two fluids coupled across their shared interface I by a rigid-lid coupling condition, i.e. no penetration and a slip with friction condition allowing a jump in the tangential velocities across I. This rigid-lid assumption is used in many oceanography models, [20, 21]. Partitioned methods allow existing, highly optimized codes for each subproblem to be used in parallel as black boxes at each time step to solve the coupled problem. The coupling between the atmosphere and the ocean dynamics occurs in the nonlinear interface condition (1.2) below. It is not clear (see [9]) if the nonlinearity in (1.2) can be treated in the most natural way (which is stable if (1.2) is linear). We give in (2.8)-(2.9) below an uncoupled and surprising (to us) discretization for the nonlinear condition (1.2) which has no analog for the linear version of (1.2) and which yields an unconditionally stable partitioned method for the fully coupled problem (1.1) - (1.6) below. The modification is based on using the geometric averaging at two time levels of the slip velocity at the interface to compute the friction coefficient rather than simply lagging it. To reduce the dynamic core of the coupled atmosphere-ocean problem to its simplest form which still retains the essential difficulty of the coupling condition, let the domain consist of two subdomains ⌦1 and ⌦2 of Rd , d = 2, 3 with outward unit normal vectors n ˆ 1 and n ˆ 2 , respectively, coupled across an interface I (example in Figure 1.1 below). The problem is: given ⌫i > 0, fi : [0, T ] ! H 1 (⌦i ), ui (0) 2 H 1 (⌦i ) and ⇤ Lawrence Livermore National Laboratory, Center for Applied Scientific Computing, Livermore, CA, 94550. Partially supported by NSF Grants DMS 0508260 and 0810385. This work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344, LLNL-JRNL-490264. email: [email protected] † Department of Mathematics and Computer Science, Clarkson University, Potsdam, NY 136995815. email: [email protected]. This material is partially based upon work supported by the Center for Nonlinear Analysis (CNA) under the National Science Foundation Grant No. DMS– 0635983. ‡ Department of Mathematics, University of Pittsburgh, Pittsburgh, PA, 15260, USA. email: [email protected]. Partially supported by NSF Grants DMS 0508260 and 0810385 .
1
2 R, find (for i = 1, 2) ui : ⌦i ⇥ [0, T ] ! Rd and pi : ⌦i ⇥ [0, T ] ! R satisfying ui,t + ui · rui
⌫i ui + rpi = fi ,
in ⌦i ,
⌫i n ˆ i · rui · ⌧ = |ui uj |(ui uj ) · ⌧, on I , for i, j = 1, 2 , i = 6 j, ui · n ˆ i = 0 on I, i = 1, 2 r · ui = 0 , in ⌦i
ui (x, 0) = u0i (x), ui = 0,
on
in ⌦i , i
= @⌦i \ I.
(1.1) (1.2) (1.3) (1.4) (1.5) (1.6)
The lateral boundary conditions on ⌦i , (1.6), are not essential for our study. We focus on (1.2). Let Xi : = {vi 2 H 1 (⌦i )d : vi = 0 on i = @⌦i \ I, vi · n ˆ i = 0 on I} ⇢ Z qi d⌦i = 0 . Qi : = qi 2 L2 (⌦i )d :
(1.7) (1.8)
⌦i
For ui 2 Xi we denote u = (u1 , u2 ) and X := {v = (v1 , v2 ) : vi 2 Xi , i = 1, 2}. Similarly, for qi 2 Qi we denote q = (q1 , q2 ) and Q := {q = (q1 , q2 ) : qi 2 Qi , i = 1, 2}. A natural subdomain variational formulation for (1.1)-(1.6), obtained by multiplying (1.1) by vi and (1.4) by qi , integrating and applying the divergence theorem, is to find (for i, j = 1, 2, i 6= j) ui : [0, T ] ! Xi and pi : [0, T ] ! Qi satisfying (ui,t , vi )⌦i + (⌫i rui , rvi )⌦i + (ui · rui , vi )⌦i (pi , r · vi )⌦i Z + |ui uj |(ui uj )vi ds = (fi , vi )⌦i , 8 vi 2 Xi
(1.9)
I
(r · ui , qi )⌦i = 0 , 8qi 2 Qi .
Let [·] denote the jump of the indicated quantity across the interface I , (·, ·) the L2 (⌦1 [ ⌦2 ) inner product and ⌫ = ⌫i and f = fi in ⌦i . The natural coupled or monolithic variational formulation for (1.1)-(1.6) is found by summing (1.9) over i, j = 1, 2 and i 6= j and is to find u : [0, T ] ! X and p : [0, T ] ! Q satisfying (ut , v) + (⌫ru, rv) + (u · ru, v) (p, r · v) Z + |[u]| [u][v]ds = (f , v), 8v 2 X,
(1.10)
I
(r · u, q) = 0 , 8q 2 Q.
Figure 1.1 illustrates the subdomains considered here, sometimes used in tests of fluid-fluid interaction, [6]. Comparing (1.10) and (1.9) we see that the monolithic problem (1.10) has a global energy that is exactly conserved, (in the appropriate sense), (set v = u and q = p in (1.10)). The subdomain sub-problems (1.9) do not posses a subdomain energy which behaves similarly due to energy transfer back and forth across the interface I. It is possible for decoupling strategies to be unstable due to the input of non-physical energy as a numerical artifact. We note that a realistic atmosphere-ocean model would contain many more terms including (typically) some eddy-viscocity term representing the e↵ects of unresolved turbulent fluctuations on resolved scales. Without such a term there is an interesting theoretical question as to whether weak solutions of (1.1) - (1.6) possess enough regularity for their traces to be well defined in L3 (I). We focus on algorithmic issues herein and suppose that the solution to (1.1) - (1.6) approximated is a strong solution. 2
n ˆ2 6
⌦1
? n ˆ1
⌦2
Fig. 1.1. Example adjoining subdomains.
1.1. Partitioned stability. The simplest and most natural partitioned method for (1.1) - (1.6) is to extrapolate from past time values the speed |ui uj | in the friction coefficient, pass the previous, known interface velocity uj across I to the problem in ⌦i and solve the problem in ⌦i . In the case of a linear coupling (without |ui uj | on the right hand side of (1.2)), this is known to be unconditionally stable, [9]; in the nonlinear case it is given as follows. For final time T and time step size t = T /N for some N 2 N, find un+1 2 Xih ⇢ Xi , pn+1 2 Qhi ⇢ Qi for i, j = 1, 2, i 6= j satisfying i i ✓
un+1 i
uni t
(pin+1 , r
, vi
◆
⌦i
+ (⌫i run+1 , rvi )⌦i + (un+1 · run+1 , vi )⌦i i i i
· vi )⌦i +
Z
I
|uni
unj |(un+1 i
unj )vi ds = (fi (tn+1 ), vi )⌦i , 8 vi 2 Xih (r · un+1 , qi )⌦i = 0 , 8qi 2 Qhi . (1.11) i
In the case of linear coupling (without |uni unj |), (1.11) has been proven in [9] to be unconditionally stable. Stability of (1.11) is an open question; while we cannot prove either stability or instability, in Section 4 we present one test where (1.11) exhibits some problems with energy dissipation. The development of a partitioned discretization method which is unconditionally and nonlinearly stable for large time steps for the nonlinear interface condition (1.2) is the main motivation of this report. We show that replacing the term unj |uni unj | in (1.11) by using a geometric averaging of the jump, unj |uni
unj |1/2 |uni
1
unj
1 1/2
|
,
gives this partitioned method. The remainder of this work is organized as follows: in Section 2, the mathematical setting and notation are described, and the partitioned time-stepping algorithm is presented. Theoretical results regarding the stability and convergence of the method are presented in Section 3. Unconditional stability of the proposed numerical method is proven and an optimal error estimate is derived. A remark is included discussing the expected dissipation of the decoupled algorithm. Numerical results that support the theoretical results are presented in section 4. Convergence rates are demonstrated for an academic test. Also, an example is provided where the unmodified IMEX discretization (1.11) exhibits a faster growth in global relative error than the geometric averaging method. 3
1.2. Related work. A standard simplification of the interface condition (1.5) used in some climate models, (e.g. [21]), is to represent the “wind shear stress” or “momentum flux” by a relationship of the form ⇢O µO ruO · n ˆ O = CD ⇢a |ua |ua where ⇢O , ⇢a are densities of the ocean and atmosphere, respectively, µO is the dynamic viscocity of the ocean, uO and ua are the velocities of the ocean and atmosphere, respectively, n ˆ O is the unit outward normal vector of the ocean domain and CD is a (dimensionless) drag coefficient. The velocity jump ua uO does not appear as it would using (1.2) since the ocean velocity is assumed to be negligible compared to the air velocity near the interface. The atmosphere and ocean are updated by an explicit update of the ocean to compute the momentum flux, which is then mapped conservatively to the atmosphere mesh. An atmosphere-ocean system, more physically refined than (1.1) - (1.6) but sharing the same mathematical structure for the dynamic core and coupling, is studied by Lions, Temam and Wang in a series of papers performing a thorough mathematical analysis [17, 18, 19]. The coupling condition (1.2) is a more physically relevant representation of the momentum flux across the interface, taking into account when the ocean surface flow may be influencing the surface air current. Bresch and Koko [6] introduce some decoupling algorithms for a model akin to (1.1)-(1.6). In it, the monolithic, coupled problem at each step is solved by preconditioned iterative methods, uncoupling can occur in the residual calculation and in the preconditioning step. Often systems of ODEs or PDEs may be decomposed into a sum of linear, sti↵ terms and other less sti↵ (possibly nonlinear) terms. Stable, consistent methods for such ODEs were derived by using implicit treatment of sti↵ terms and explicit treatment of the remaining terms, the so called additive Runge-Kutta methods, [10]. Kennedy and Carpenter exhibited the utility of this idea for convection-di↵usionreaction equations, [14], and in [3] implicit-explicit (IMEX) methods of up to fourth order have been developed for some classes of PDEs. These methods seek to reduce the time step restriction on an explicit method for stability by implicit treatment of sti↵ terms while reducing computational cost through explicit updating of other terms, particularly nonlinear terms. One application of the IMEX approach is to use explicit updates of the terms responsible for coupling across I, giving semi-implicit partitioned methods. A straightforward partitioning approach for some coupled problems may result in an unstable or conditionally stable method, as demonstrated for the fluid-structure interaction problem, see [7, 8]. 2. Method Description, Notation and Preliminaries. This section presents the numerical scheme for (1.1)-(1.6), and provides the necessary definitions and lemmas for the stability and convergence analysis. For D ⇢ ⌦, the Sobolev space H k (D) = W k,2 (D) is equipped with the usual norm k·kH k (D) , and semi-norm |·|H k (D) , for 1 k < 1, e.g. Adams [2]. The L2 norm is denoted by k·kD . For functions v(x, t) defined for almost every t 2 (0, T ) on a function space V (D), we define the norms (1 p 1) !1/p Z kvkL1 (0,T ;V ) = ess sup kv(·, t)kV
and
0 0, f 2 H 1 (⌦). Given n 1 u , un 2 Vh , n 2 {1, 2, · · · , N 1}, find un+1 2 V1,h satisfying 1 ✓
un+1 1
un1 t
, v1
◆
⌦1
+ ⌫1 run+1 , rv1 1 +
Z
I
⌦1
+ c1 un+1 ; un+1 , v1 1 1
|[un ]|un+1 · v1 ds 1
= f1 (tn+1 ), v1
⌦1
Z
I
|[un ]|1/2 |[un
1 1/2 n ]| u2
, 8v1 2 V1,h .
(2.8)
and un+1 2 V2,h satisfying 2 ✓
un+1 2
un2 t
, v2
◆
⌦2
+ ⌫2 run+1 , rv2 2 +
Z
I
⌦2
+ c2 un+1 ; un+1 , v2 2 2
|[un ]|un+1 · v2 ds 2
= f2 (tn+1 ), v2
⌦2
· v1 ds
Z
I
|[un ]|1/2 |[un
1 1/2 n ]| u1
· v2 ds
, 8v2 2 V2,h .
(2.9) The main feature of this algorithm is the use of the geometric average |[un ]|1/2 |[un 1 ]|1/2 in the coupling terms. This will be shown to have a stabilizing e↵ect for the algorithm. Hence Algorithm 2.1 will often be referred to as the geometric averaging method, or GA method. 3. Numerical Analysis. Stability of the partitioned scheme (Algorithm 2.1) is established below. The polarization identity is used in the stability and convergence proofs: (wj+1
wj ) · wj+1 =
1 |wj+1 |2 + |wj+1 2 7
wj |2
|wj |2
(3.1)
Lemma 3.1. (Partitioned Stability) Let uj1 2 X1,h satisfy (2.8) for each j 2 {0, 1, 2, · · · , n 1, n}, and uj2 2 X2,h satisfy (2.9) for each j 2 {0, 1, 2, · · · , n 1, n}. Then Algorithm 2.1 satisfies the following energy equality at time step n + 1: un+1
2
+
n X
uj+1
j=1
+
t
+
2
uj
t
+2 t
j=1
n Z X
j=1 I n Z X j=1
Z
n ✓ X
I
2
⌫1 ruj+1 1
⌦1
2
+ ⌫2 ruj+1 2
|[uj ]|1/2 uj+1 1
|[uj
1 1/2 j ]| u2
|[uj ]|1/2 uj+1 2
|[uj
1 1/2 j ]| u1
2
⌦2
◆
ds
2
ds
⇣ ⌘ 2 2 |[un ]| un+1 + un+1 ds 1 2 I Z ⇣ ⌘ 2 2 2 + t [u0 ] u11 + u12 ds
+ t = u1 +2 t
n ⇢⇣ X
I
f1 (tj+1 ), uj+1 1
j=1
⌘
⌦1
⇣ ⌘ + f2 (tj+1 ), uj+1 2
. (3.2)
⌦2
Furthermore, the algorithm is unconditionally stable, satisfying: n+1 2
u
+
n X
j 2
u
+
u
j+1
j=1
ruj+1 1
⌫1
j=1
Z
2 ⌦1
ruj+1 2
+ ⌫2
⇣ ⌘ 2 2 |[un ]| un+1 + un+1 ds 1 2 I Z ⇣ ⌘ 2 2 2 + t [u0 ] u11 + u12 ds
+ t u1
t
n ✓ X
2 ⌦2
◆
I
n ✓ X t + f1 (tj+1 ) ⌫ 1 j=1
2 H
1 (⌦
1)
+
t
⌫2
f2 (t
j+1
)
2 H
1 (⌦
2)
◆
. (3.3)
Remark 3.1 (Discrete energy). The discrete energy equation (3.2) includes the following exact “kinetic energy” of the discretization, denoted by KE n+1 : Z ⇣ ⌘ 2 2 2 KE n+1 = un+1 + t |[un ]| un+1 + un+1 ds. 1 2 I
The corresponding “energy dissipation” of the method, denoted by ✏n+1 is ✏n+1 = 2 t
n ✓ X j=1
⌫1 ruj+1 1
+
t
2 ⌦1
n Z X j=1
I
2
+ ⌫2 ruj+1 2
⌦2
|[uj ]|1/2 uj+1 1 +
t
|[uj
n Z X j=1
8
◆
I
1 1/2 j ]| u2
|[uj ]|1/2 uj+1 2
2
ds |[uj
1 1/2 j ]| u1
2
ds.
Proof. . Let j 2 {1, 2, . . . , N v1 = uj+1 1 , (2.8) yields
uj+1 1
uj1 t
, u1j+1
!
+ ⌫1 ruj+1 1
⌦1
1} and consider (2.8) with n = j. Setting
Z
I
2 ⌦1
+
Z
2 j |uj+1 1 | |[u ]| ds
I
j 1/2 uj2 · uj+1 |[uj 1 |[u ]|
]|
1 1/2
⇣ ⌘ ds = f1 (tj+1 ), uj+1 1
⌦1
. (3.4)
Similarly, choosing v2 = uj+1 in (2.9) yields 2
uj+1 2
uj2 t
, u2j+1
!
ruj+1 2
+ ⌫2
⌦2
Z
I
2 ⌦2
+
Z
2 j |uj+1 2 | |[u ]| ds
I
j 1/2 uj1 · uj+1 |[uj 2 |[u ]|
]|
1 1/2
⇣ ⌘ ds = f2 (tj+1 ), uj+1 2
⌦2
. (3.5)
Applying the identity (3.1) to the first terms on the left hand side of (3.4)-(3.5) and adding, it follows that
1 2 t
✓
+ ⌫1
2
uj+1 1
⌦1
ruj+1 1
Z
I
2 ⌦1
uj1 + ⌫2
uj2 · uj+1 [uj ] 1
2 ⌦1
+
uj+1 1
ruj+1 2 1/2
[uj
2 ⌦2
uj1
+
1 1/2
Z
I
2 ⌦1
uj+1 2
+
uj+1 1
2
2 ⌦2
2
uj2
[uj ] ds +
Z
⌦2
I
Z
+
uj+1 2
uj+1 2 2
1/2
uj1 · uj+1 [uj ] [uj 1 ] 2 ⇣ ⌘ ⇣ ⌘ = f1 (tj+1 ), uj+1 + f2 (tj+1 ), uj+1 1 2 ]
ds
2 ⌦2
[uj ] ds 1/2
ds
I
⌦1
The following two interface integrals may be expressed in a di↵erent way: Z 2 1/2 1/2 j uj+1 [u ] ds uj2 · uj+1 [uj ] [uj 1 ] ds 1 1 I I Z ⇣ ⌘ j 1/2 j 1/2 = uj+1 · uj+1 uj2 |[uj 1 ]|1/2 ds 1 |[u ]| 1 |[u ]| ZI Z 2 2 j = uj+1 [u ] ds uj2 [uj 1 ] ds 1 2 I 2 I Z 2 j 1/2 + uj+1 uj2 |[uj 1 ]|1/2 ds. 1 |[u ]| 2 I
uj2
Z
9
⌦2
. (3.6)
◆
This expression is inserted into (3.6) and the remaining two interface integrals are treated analogously:
◆ 2 2 2 2 j j+1 j uj+1 u + u u 1 1 2 2 2 t ⌦1 ⌦1 ⌦2 ⌦2 ◆ ✓ 2 2 2 1 + uj+1 uj1 + uj+1 uj2 + ⌫1 ruj+1 + ⌫2 ruj+1 1 2 1 2 2 t ⌦1 ⌦2 ⌦1 ⇢ Z ⇢ Z 2 2 2 2 [uj ] ds [uj 1 ] ds + uj+1 + uj+1 uj1 + uj2 1 2 2 I 2 I Z 2 j 1/2 + uj+1 uj2 |[uj 1 ]|1/2 ds 1 |[u ]| 2 I Z 2 j 1/2 + uj+1 uj1 |[uj 1 ]|1/2 ds 2 |[u ]| 2 ⇣ I ⌘ ⇣ ⌘ j+1 j+1 = f1 (tj+1 ), uj+1 + f (t ), u . 2 1 2 1
✓
⌦1
2 ⌦2
(3.7)
⌦2
Multiply through (3.7) by 2 t, and summing over j = 1, 2, . . . , n proves the energy equation (3.2). To prove (3.3), bound the right hand side of (3.7) by:
⇣
f1 (tj+1 ), uj+1 1
⌘
⌦1
⇣ ⌘ + f2 (tj+1 ), uj+1 2
1 f1 (tj+1 ) 2⌫1
2 H
⌦2
1 (⌦
1)
1 2 f2 (tj+1 ) H 1 (⌦ ) 2 2⌫2 2 ⌫1 ⌫2 + ruj+1 + ruj+1 1 2 2 2 ⌦1 +
2 ⌦2
. (3.8)
Insert (3.8) into (3.7), and the desired result follows by multiplying through by 2 t and summing over j = 1, 2, . . . , n. Theorem 3.2 (Convergence of Algorithm 2.1). Let uk1 2 X1,h satisfy (2.8) for each k 2 {2, · · · , n N 1}, and uk2 2 X2,h satisfy (2.9) for each k 2 {2, · · · , n N 1}. Let ⌫˜ = max{⌫1 1 , ⌫2 1 }, ⌫ˆ = max{⌫1 , ⌫2 }, and let ⇣ Dn+1 = ⌫˜3 1 + 4 E n+1 + ru(tn+1 )
4
⌘
,
where E n+1 = maxj=0,1,...,n+1 {max{ku(tj )k4I , kuj k4I }}. Assume t 1/Dn+1 , and that (u, p) is a strong solution of the coupled NSE system (1.1)-(1.6) with ut 2 L2 (0, T ; X) and utt 2 L2 (0, T ; L2 (⌦)). Then the solution un+1 of Algorithm 2.1 10
satisfies: u(tn+1 )
un+1
t + 2
⌫1
+ + +
t
j X
n=1
t ⌫1 +
n X
k=1
C exp +
2
✓
1
r(u1 (t
k+1
Dn+1 t Dn+1
r(u1 (t1 )
)
!(
u11 )
v11 2V1,h
inf
v21 2V2,h
r(u2 (t ) 1
2 ⌦1
2 ⌦2
u(t1 )
+
+
2 v21 ) ⌦ 2
2
+ 1inf
v 2Vh
k+1
2 uk+1 ) ⌦ 2 2
)
u(t1 )
1 2 r(u2 (t0 ) u02 ) ⌦ 2 2 1 + inf r(u2 (t0 ) 2 v20 2V2,h
v2Vh
2
2
u1
r(u2 (t
1 2 r(u1 (t0 ) u01 ) ⌦ 1 2 1 + inf r(u1 (t0 ) 2 v10 2V1,h
t2 kutt kL2 (0,T ;L2 (⌦)) + inf k(u + 2 t2 kut kL2 (0,T ;X) + T
+ ⌫2
n X
k=1
2 v11 ) ⌦ 1
inf r(u1 (t ) ✓ t ⌫2 r(u2 (t1 ) u12 ) 1
2 uk+1 ) ⌦ 1 1
v1
2 v10 ) ⌦ 1
◆
2 v20 ) ⌦ 2
◆
2
v)t kL2 (0,T ;L2 (⌦)) + inf kp
max
k=2,...,n+1
✓
q2Qh
inf
vk 2Vh
r(u(tk )
vk )
2
2
2
qkL2 (0,T ;L2 (⌦)) ◆) ,
(3.9)
where C has the following dependence on , ⌫1 , and ⌫2 : C = O max{ˆ ⌫ , ⌫˜, (1 +
!
t4 )˜ ⌫ 3 , 2 } .
Remark 3.2. The timestep restriction is for the error estimate and not stability. It comes from the discrete Gronwall inequality and is typical for nonlinearly implicit methods. Our stability and error analysis can readily be extended to the linearly implicit method (where un+1 · run+1 is replaced by uni · run+1 ). See Ingram [13] for a i i i detailed analysis of this e↵ect. Proof. To prove convergence we will consider the error equations on each subdomain independently first. The resulting error bounds from the two subdomains will then be summed. Begin by taking the variational formulation for the true solution u1 (tn+1 ) on ⌦1 and subtracting (2.8), 1 (u1 (tn+1 ) un+1 ) (u1 (tn ) un1 ), v1 ⌦ + ⌫1 r(u1 (tn+1 ) un+1 ), rv1 1 1 1 t n+1 n+1 n+1 n+1 n+1 + c1 u1 (t ); u1 (t ), v1 c1 u1 ; u1 , v1 p1 (t ), r · v1 ⌦ 1 Z Z n+1 n+1 n+1 n + u1 (t )|[u(t )]| · v1 ds u1 |[u ]| · v1 ds I ZI Z n n 1/2 n 1 1/2 + u2 |[u ]| |[u ]| · v1 ds u2 (tn+1 )|[u(tn+1 )]| · v1 ds I ✓ ◆I u1 (tn+1 ) u1 (tn ) @u1 (tn+1 ) , v1 , 8v1 2 V1,h . = t @t ⌦1
Errors are decomposed using the equation ui (tj )
uji = (ui (tj )
u ˜ji ) + (˜ uji 11
uji ) = ⌘ij +
j i
⌦1
where ⌘ij 2 Vi and ji 2 Vh,i with u ˜ji 2 Vh,i arbitrarily chosen. We also subtract an j arbitrarily chosen qi 2 Qh,i from the pressure for free. Thus choosing v1 = n+1 the 1 above error equation is rewritten as 1 ⌘ n+1 ⌘1n , n+1 1 t 1 n+1 2 + ⌫1 kr 1 k⌦1
⌦1
+
1 t
n+1 1
p1 (tn+1 )
n 1,
n+1 1 ⌦1
+ ⌫1 r⌘1n+1 , r
q1n+1 , r ·
n+1 1 ⌦1 n+1 n+1 c1 u1 ; u1 ,
n+1 + c1 u1 (tn+1 ); u1 (tn+1 ), n+1 1 1 Z Z + u1 (tn+1 )|[u(tn+1 )]| · n+1 ds un+1 |[un ]| · n+1 ds 1 1 1 I I Z Z + un2 |[un ]|1/2 |[un 1 ]|1/2 · n+1 ds u2 (tn+1 )|[u(tn+1 )]| · 1
=
I R1n+1 ,
I
n+1 1 ⌦1
n+1 1 ⌦1
n+1 1
ds
,
u (tn+1 ) u (tn )
@u (tn+1 )
j j , for j = 1, 2. The goal is to bound n+1 , where Rjn+1 = j 1 t @t hence we move many terms to the right hand side. Apply (3.1) on the left hand side:
1 n+1 2 n 2 k n+1 k2 k n1 k2⌦1 + k n+1 k⌦1 1 k⌦1 + ⌫1 kr 1 1 2 t ⇢1 Z ⌦1 Z + u1 (tn+1 )|[u(tn+1 )]| · n+1 ds un+1 |[un ]| · n+1 ds 1 1 1 I I ⇢ Z Z + un2 |[un ]|1/2 |[un 1 ]|1/2 · n+1 ds u2 (tn+1 )|[u(tn+1 )]| · 1 I
1 ⌘ n+1 ⌘1n , n+1 1 ⌦1 t 1 n+1 n+1 n+1 + p1 (t ) q1 , r · 1 =
c1 u1 (t
n+1
); u1 (t
n+1
),
n+1 1
I
⌫1 r⌘1n+1 , r
n+1 1
ds
n+1 1 ⌦1
⌦1
+ c1 un+1 ; un+1 , 1 1
n+1 1
+ R1n+1 ,
n+1 1 ⌦1
. (3.10) Bounds for the non-interface terms follow as in the standard NSE case, (see e.g. [15]). The interface terms require special treatment. First, note that by adding and subtracting terms in sequence one may derive the following expression: Z Z u1 (tn+1 )|[u(tn+1 )]| · n+1 ds un+1 |[un ]| · n+1 ds 1 1 1 I I Z = u1 (tn+1 ) |[u(tn+1 )]| |[u(tn )]| · n+1 ds 1 ZI + u1 (tn+1 ) (|[u(tn )]| |[˜ un ]|) · n+1 ds 1 I Z + u1 (tn+1 ) (|[˜ un ]| |[un ]|) · n+1 ds 1 I Z Z + ⌘1n+1 |[un ]| · n+1 ds + | n+1 |2 |[un ]| ds. 1 1 I
I
(3.11) The splitting used for the remaining pair of interface integrals is more complicated due to the product |[un ]|1/2 |[un 1 ]|1/2 . The key is to recognize an error of O( t) is committed in replacing the product |[u(tn )]|1/2 |[u(tn 1 )]|1/2 with the average 21 (|[u(tn )]|+ 12
|[u(tn
1
)]|). We proceed as follows.
Z un2 |[un ]|1/2 |[un 1 ]|1/2 · n+1 ds u2 (tn+1 )|[u(tn+1 )]| · n+1 ds 1 1 I I ⇢ Z 1 (|[un ]| + |[un 1 ]|) · n+1 ds = un2 |[un ]|1/2 |[un 1 ]|1/2 1 2 I Z 1 + un2 (|[un ]| + |[un 1 ]|) (|[˜ un ]| + |[˜ un 1 ]|) · n+1 ds 1 2 I Z 1 + un2 (|[˜ un ]| + |[˜ un 1 ]|) (|[u(tn )]| + |[u(tn 1 )]|) · n+1 ds 1 2 I ⇢ Z 1 + un2 (|[u(tn )]| + |[u(tn 1 )]|) |[u(tn+1 )]| · n+1 ds 1 2 I Z Z n+1 n n+1 |[u(t )]| · ds ⌘2n |[u(tn+1 )]| · n+1 ds 2 1 1 I I Z + (u2 (tn ) u2 (tn+1 ))|[u(tn+1 )]| · n+1 ds 1
Z
(3.12)
I
Using (3.11)–(3.12), substitution for the interface integrals in (3.10) is performed. The last term of (3.11) is positive and may be kept on the left hand side of (3.10), but all other interface terms are now moved to the right hand side, and bounded. Applying the reverse triangle inequality and the equality [a] [b] = [a b] for a, b 2 X, we first note |[u(tj )]|
|[˜ uj ]| |[⌘ j ]|
and
|[˜ uj ]|
j
|[uj ]| |[
]|.
Thus (3.10)–(3.12) yield the following error inequality: 1 2 t
k
n+1 2 k⌦ 1 1
k
n 2 1 k⌦ 1
1 ⌘ n+1 ⌘1n , n+1 1 t 1 n+1 n+1 c1 u1 (t ); u1 (t ),
+k
n+1 1
n 2 1 k⌦1
+ ⌫1 kr
⌦1
⌫1 r⌘1n+1 , r
n+1 1
+ c1 un+1 ; un+1 , 1 1
n+1 2 k⌦ 1 1
+
Z
+ p1 (tn+1 )
n+1 1 ⌦1
n+1 1
+ R1n+1 ,
I
|
n+1 2 | |[un ]| ds 1
q1n+1 , r ·
n+1 1 ⌦1
n+1 1 ⌦1
+ I1 + I2 + I3 + I4 ,
(3.13) where the Ij -terms denote the following interface integrals. Z Z n+1 n+1 n+1 n I1 = |u1 (t )||[u(t ) u(t )]|| 1 | ds + |u1 (tn+1 )||[⌘ n ]|| n+1 | ds 1 I Z ZI + |u2 (tn ) u2 (tn+1 )||[u(tn+1 )]|| n+1 | ds + |⌘1n+1 ||[un ]|| n+1 | ds 1 1 I I Z Z 1 + |un2 |(|[⌘ n ]| + |[⌘ n 1 ]|)| n+1 | ds + |⌘2n ||[u(tn+1 )]|| n+1 | ds 1 1 2 I I Z 1 |un2 |(|[ n ]| + |[ n 1 ]|)| 2 I Z + | n2 ||[u(tn+1 )]|| n+1 | ds 1
I2 =
n+1 | ds 1
I
13
+
Z
I
|u1 (tn+1 )||[
n
]||
n+1 | ds 1
Z
I3 =
I
I4 =
Z
I
|un2 |
1 (|[u(tn )]| + |[u(tn 2
|un2 | |[un ]|1/2 |[un
]|
1 1/2
1
)]|)
|[u(tn+1 )]| |
1 (|[un ]| + |[un 2
1
n+1 | ds 1
]|) |
n+1 | ds 1
The application of H¨ older’s inequality, Young’s inequality, and (2.2) to the noninterface terms of (3.13) gives the estimate Z 1 n+1 n+1 2 2 n 2 n 2 k n+1 k k k + k k + ⌫ kr k + | n+1 |2 |[un ]| ds 1 ⌦1 1 ⌦1 1 ⌦1 ⌦1 1 1 1 1 2 t I ⌫1 1 d 2 4 n+1 2 r⌘1n+1 ⌦ + 3 3 ru1 (tn+1 ) ⌦ + kp1 (tn+1 ) q1n+1 k2L2 (⌦1 ) 1 ⌦1 1 1 2"2 4"3 ⌫1 2"2 ⌫1 ✓ ◆ ✓ n+1 2 3"3 C ⌘1 ⌘1n 2 2 + 2"1 + "2 + ⌫1 r n+1 + + R1n+1 V 0 1 ⌦1 1 4 2"1 ⌫1 t 0 V1 ◆ 2 run+1 r⌘1n+1 ⌦1 + kru1 (tn+1 )k2⌦1 kr⌘1n+1 k2⌦1 + un+1 1 1 ⌦1 ⌦1 + I1 + I2 + I3 + I4 ,
(3.14) for constants C and "i , i = 1, 2, 3 independent of ⌫1 . The interface integrals are bounded using Lemma 2.2, with particular attention paid to the constants involving 2 n+1 2 , ⌫1 , ⌫2 and terms with r n+1 or . In the term I1 , each piece is 1 1 ⌦1 ⌦1 bounded using (2.5) with ↵1 = 192. Then the first two terms of I2 are bounded using (2.7) with ↵1 = 16, ↵2 = 32, 1 = 96 and the last term using (2.6) with ↵2 = 8 and 1 = 48. Some basic inequalities can be used to simplify I3 , I4 and derive corresponding error estimates. Applying the reverse triangle inequality for I3 yields 1 (|[u(tn )]| + |[u(tn 1 )]|) |[u(tn+1 )]| 2 1 1 |[u(tn )]| |[u(tn+1 )] + |[u(tn 1 )]| |[u(tn+1 )]| 2 2 1 1 n n+1 [u(t ) u(t )] + [u(tn 1 ) u(tn+1 )] . 2 2 I4 is bounded by first noting |[un ]|1/2 |[un
]|
1 1/2
Add and subtract u(tn ) |[un ]|1/2 |[un
]|
1 1/2
(3.15)
2 1 1 (|[un ]| + |[un 1 ]|) = |[un ]|1/2 |[un 1 ]|1/2 2 2 1 |[un ]|1/2 |[un 1 ]|1/2 |[un ]|1/2 + |[un 1 ]|1/2 2 1 1 = |[un ]| |[un 1 ]| [un un 1 ] . 2 2
u(tn
1
),
1 (|[un ]| + |[un 2
1
]|)
1 1 [(un u(tn 1 )) + (u(tn 1 ) un 1 )] + [u(tn ) u(tn 1 )] 2 2 o 1n n |[ ]| + [ n 1 ] + |[⌘ n ]| + [⌘ n 1 ] + [u(tn ) u(tn 1 )] . 2
14
(3.16)
The results (3.15)–(3.16) are then bounded using the same approach as particular terms in I1 and I2 . For example, I3 can be bounded in the same manner as the first term from I1 , while the first term from I4 is bounded the same way as the first term of I2 . Thus, the interface terms are bounded by I1 + I2 + I3 + I4 ⇣ ⇣ ⌫1 2 4 3 4 r n+1 + C⌫ 1 + 48 kun2 kI + 12 [u(tn+1 )] 1 1 ⌦1 16 ⇣ ⌘ ⌫ ⇣ ⌫1 2 2 2 2 + kr n1 k⌦1 + r n1 1 ⌦ + kr n2 k⌦2 + r n2 1 1 16 16 n ⇣ ⌘ ⇣ n 2 1 k⌦ 1
+ C4 M n+1 ⌫1 3 k
n 1 2 1 ⌦1
+
+ ⌫2 3 k
n 2 2 k⌦ 2
+
4 I
⌘⌘
2 ⌦2
n 2
n+1 2 1 ⌦1
⌘
+ C2 Ln+1 P n+1 ⌘o 1 2 , ⌦ 2
(3.17) where C is a constant independent of , ⌫1 , and ⌫2 . The quantities Ln+1 , M n+1 and P n+1 are defined by: Ln+1 = u1 (tn+1 ) M n+1 = u1 (tn+1 )
2 I 4 I
and P n+1 = u1 (tn+1 )
+ u2 (tn+1 ) + u2 (tn+1 ) 2 I
u1 (tn )
+ u1 (tn+1 )
2
2
4
4
+ kun1 kI + kun2 kI ,
+ kun1 kI + kun2 kI ,
+ u2 (tn+1 )
u1 (tn
2
2 I 4 I
1
2 I
)
2
+ k⌘1n kI + k⌘2n kI + ⌘1n
u2 (tn )
2 I
+ u2 (tn+1 )
1 2 I
1 2 I
+ ⌘2n
u2 (tn
1
+ ⌘1n+1
)
2 I
2 I
+ ⌘2n+1
2 I
.
(The last term in the definition of P n+1 is auxiliary, it is included in the definition for use in the estimate on ⌦2 .) Then, with choices for "1 , "2 , "3 satisfying 2"1 +"2 +3"3 /4 = 1/16, (3.14) can be bounded as follows: 1 2 t
n+1 2 k⌦1 1
k
k
n 2 1 k⌦ 1
+k
n+1 1
7⌫1 kr 8
+
n 2 1 k⌦ 1
n ⇣ 4 C⌫1 3 1 + kru1 (tn+1 )k4⌦1 + 4 kun2 kI + [u(tn+1 )]
+ C⌫1
1
2
⌘1n
+
t
V10
kR1n+1 k2V10
+ kp1 (t
n+1
)
4 I
⌘o
n+1 2 1 ⌦1
q1n+1 k2⌦1
!
⌘ ⌫ ⇣ ⌘ 2 2 2 2 + r n1 1 ⌦1 + kr n2 k⌦2 + r 2n 1 ⌦2 + C2 Ln+1 P n+1 16 ⌘ n ⇣ ⇣ ⌘o 3 n 1 2 3 n 1 2 4 n+1 n 2 n 2 + C M ⌫1 k 1 k⌦ 1 + 1 + ⌫ k k + 2 ⌦2 2 2 ⌦1 ⌦2 ⇣ ⌘ 1 1 n+1 n+1 n+1 2 n+1 2 + C ⌫1 + ⌫1 kru1 (t )k⌦1 + ⌫1 u1 ru1 kr⌘1 k⌦1 . ⌦1 ⌦1 +
⌫1 ⇣ kr 16
⌘1n+1
n+1 2 k⌦1 1
n 2 1 k⌦ 1
(3.18)
Let Sn+1 denote Sn+1 = u1
2
+ t
Z
I
[u0 ]
⇣
u11
2
+ u12
2
⌘
ds
n ✓ X t + f1 (tj+1 ) ⌫ 1 j=1
15
2 V10
+
t ⌫2
f2 (t
j+1
)
2 V20
◆
, (3.19)
i.e., the right-hand side of the stability estimate (3.3). Also, let ⌫ˆ := max{⌫1 , ⌫2 } and ⌫˜ := max{⌫1 1 , ⌫2 1 }. Then (3.18) can be simplified to 1 7⌫1 n 2 k 1n+1 k2⌦1 k n1 k2⌦1 + k n+1 kr n+1 k2⌦1 1 k⌦ 1 + 1 1 2 t 8 n ⇣ ⌘o 4 4 C⌫1 3 1 + kru1 (tn+1 )k4⌦1 + 4 kun2 kI + [u(tn+1 )] I + C⌫1
⌘1n+1
1
2
⌘1n
+
t
⌫1 ⇣ + kr 16
V10
kR1n+1 k2V10
)
q1n+1 k2⌦1
!
(3.20)
⌘ ⌫2 ⇣ n 1 2 n 2 n 2 k + r + kr k + r 1 ⌦1 2 ⌦2 2 ⌦2 ⇣ 16 ⌘ n 2 n 1 2 2 n+1 n+1 4 n+1 3 + C L P + C M ⌫˜ k k + ⇣ ⌘ 1/2 + C ⌫1 + ⌫1 1 kru1 (tn+1 )k2⌦1 + ⌫1 1 Sn+1 run+1 kr⌘1n+1 k2⌦1 . 1 ⌦1 n 1 2 1 ⌦1
⌘
+ kp1 (t
n+1
n+1 2 1 ⌦1
Proceeding further requires estimates for the solution on ⌦2 . This analysis is done in the same manner as above for ⌦1 using (2.9), and a similar estimate is derived: 1 2 t
n+1 2 k⌦ 2 2
k
n 2 2 k⌦2
k
+k
n+1 2
n 2 2 k⌦ 2
+
7⌫2 kr 8
n ⇣ 4 C⌫2 3 1 + kru2 (tn+1 )k4⌦2 + 4 kun1 kI + [u(tn+1 )]
+ C⌫2
⌘2n+1
1
⌫1 ⇣ + kr 16
n+1 2 k⌦2 2 4 I
2
⌘2n
+
t
V20
kR2n+1 k2V20
)
q2n+1 k2⌦2
n+1 2 2 ⌦2
!
(3.21)
⌘ ⌫2 ⇣ 2 2 + r + kr n2 k⌦2 + r n2 1 ⌦ 2 ⇣ 16 ⌘ n 2 n 1 2 2 n+1 n+1 4 n+1 3 + C L P + C M ⌫˜ k k + ⇣ ⌘ 1/2 + C ⌫2 + ⌫2 1 kru2 (tn+1 )k2⌦2 + ⌫2 1 Sn+1 run+1 kr⌘2n+1 k2⌦2 . 2 ⌦2 n 1 2 1 ⌦1
n 2 1 k⌦ 1
⌘
+ kp2 (t
n+1
⌘o
Combining (3.20) and (3.21) gives the global estimate 1 k 2 t
n+1 2
k
k
n 2
k +k
n+1
n 2
k
+ C ⌫˜
t
n+1 2 k⌦ 1 1
n 2 1 k⌦ 1
+
7⌫2 kr 8
n+1 2 k⌦2 2
2
+ V0
kR1n+1 k2V10
+
+ r
n 1 2 1 ⌦1
⌘
kR2n+1 k2V20
+ kp(t
n+1
⌫2 ⇣ 2 + kr n2 k⌦2 + r 2n 8 ⇣ ⌘ 2 2 + C2 Ln+1 P n+1 + C4 M n+1 ⌫˜3 k n k + n 1 ⇣ ⌘ 1/2 + C ⌫ˆ + ⌫˜kru(tn+1 )k2 + ⌫˜Sn+1 run+1 kr⌘ n+1 k2 . +
⌫1 ⇣ kr 8
⌘n
7⌫1 kr 8
n+1 2
C ⌫˜3 1 + kru(tn+1 )k4 + 4 M n+1 ⌘ n+1
+
16
)
1 2 ⌦2
q
n+1 2
⌘
k
!
(3.22)
Rearranging terms gives 1
⌫1 ⌫2 2 k + r n+1 + r 1 ⌦ 1 2 2⌘ ⌘ ⌫ ⇣ 2 1 2 2 kr n1 k⌦1 + kr n1 k⌦1 r 1n 1 ⌦ 1 8 ⌘ ⌫ ⇣ ⌘ 2 n 1 2 n 2 n 2 kr 2 k⌦2 + kr 2 k⌦2 r 2 ⌦2 8
k n+1 k2 k 2 t ⇣ ⌫1 2 + r n+1 1 ⌦1 4 ⇣ ⌫2 2 + r n+1 2 ⌦2 4
n 2
n+1
k +k
n 2
n+1 2 2 ⌦2
n+1 2
C An+1 + C B n+1 + C ⌫˜3 1 + kru(tn+1 )k4 + 4 M n+1 ⇣ ⌘ 2 2 + C2 Ln+1 P n+1 + C4 M n+1 ⌫˜3 k n k + n 1
(3.23)
where An+1 and B n+1 are defined by
⌘ n+1
An+1 = ⌫˜
⌘n t
2
+ kR1n+1 k2V10 + kR2n+1 k2V20 + kp(tn+1 )
V0
qn+1 k2
⇣ ⌘ 1/2 B n+1 = ⌫ˆ + ⌫˜kru(tn+1 )k2 + ⌫˜Sn+1 run+1 kr⌘ n+1 k2 .
!
(3.24) Summing over n = 1, 2, . . . , p 1, multiplying through by 2 t, and rearranging gives a global error estimate at time tp :
k
p 2
k +
+
p 1 X
n 2
n+1
+
⌫1
t ⌫1 4
✓
n=1
2 kr
p 2 1 k⌦ 1
+ r
1 + C t ⌫˜ (M + M ) ✓ t ⌫1 1 2 + r 11 ⌦ + r 1 2 2 +C
t
n=1
t
p 1 X
p 1 X
3 4
2
p 1 1
2 ⌦1
1 2
3
0 1
◆
+
r
n+1 2 1 ⌦1
t ⌫2 4
✓
2 kr
+ r
1 2 2 ⌦2 p 1 X
1 + r 2
⌦2
◆
(3.25)
where Dn+1 is defined by 4
◆
n+1 2
Dn+1
n=1
⇣ Dn+1 = ⌫˜3 1 + 4 E n+1 + ru(tn+1 )
0 2 2 ⌦2
!
2
p 1 2
0 2
2
An+1 + B n+1 + 2 Ln+1 P n+1 + C t
n=1
p 2 2 k⌦ 2
n+1 2 2 ⌦2
r
n=1
+ C t ⌫˜ M ✓ ◆ t ⌫2 2 + r ⌦1 2 3 4
+ ⌫2
p 1 X
⌘
(3.26)
,
with
E
n+1
8 n+3 > + M n+2 + M n+1 <M = M n+2 + M n+1 > : n+1 M 17
if 1 n p if n = p 2 if n = p 1
3 .
(3.27)
To use the discrete Gronwall’s inequality (see [15]), note that for j = 1, 2, . . . , N , write (3.25) as aj+1 +
j X
bn+1 +
t
n=1
j X
cn+1
n=1
J +C t
j X
Dn+1 an+1 + C
t
n=1
j X
An+1 + B n+1 + 2 Ln+1 P n+1
(3.28)
n=1
where aj+1 =
j+1 2
j+1
bj+1 =
,
j 2
cj+1 = ⌫1 r
,
j+1 1
2 ⌦1
+ ⌫2 r
2
j+1 2
⌦2
,
and J = 1 + C t ⌫˜3 4 (M 2 + M 3 ) ✓ t ⌫1 1 2 + r 11 ⌦ + r 1 2 2
1 2 0 1
+ C t ⌫˜3 4 M 2 ✓ ◆ t ⌫2 2 + r 12 ⌦1 2
0 2 2 ⌦2
1 + r 2
0 2 2 ⌦2
◆
.
Applying Gronwall’s inequality, (3.28) gives j+1
a
+
j X
b
n+1
+
n=1
exp C t
t
j X
cn+1
n=1 j X
n=1
n+1
!(
J +C t
j X
A
n+1
+B
n+1
+ L
n=1
2
n+1
P
n+1
)
(3.29)
where n+1
=
1
Dn+1 . t Dn+1
(3.30)
It remains to bound the individual terms on the right-hand side of (3.29). The terms in J are bounded using the triangle inequality. The terms in An+1 are bounded in the usual way (see the proof of Theorem 4.1 of [9], for example). In B n+1 , the summed factor run+1 is bounded using the stability estimate (3.3) and the ru(tn+1 ) term is bounded by a standard a priori estimate on the solution to the continuous problem. The terms in Ln+1 and M n+1 are bounded by a standard a priori estimate on u(tn+1 ) or using S n for the discrete solutions, (the bound (2.4), and the stability estimate (3.3)). The terms in P n+1 are bounded the same way. The estimate (3.9) is derived through the application of the triangle inequality to the left hand side of (3.29). Corollary 3.1 (Convergence rate for MINI-element,[1]). Let (u, p) be a solution of (1.1)-(1.6) satisfying the assumptions of Lemma 3.1 and Theorem 3.2, approximated using the MINI-element on a mesh of maximum element width h by Algorithm 2.1. Assume the velocity data u0 , u1 satisfies ku(t0 )
u0 kX + ku(t1 ) 18
u1 kX C1 h,
for a generic constant C1 independent of t, h. Then there exists C > 0 independent of h, t such that the discrete velocity satisfies the optimal error estimate t
⌫1
n X
k=1 2
r(u1 (tk+1 )
C (h +
uk+1 ) 1
2 ⌦1
+ ⌫2
n X
k=1
t2 ).
r(u2 (tk+1 )
uk+1 ) 2
2 ⌦2
!
Proof. This follows by using standard interpolation error estimates for the MINIelement, e.g. [1, 5, 15], in the right hand side of the error bound (3.9). 4. Numerical Experiments. In this section numerical experiments are presented that support the theory described in this report. An example is provided to verify the asymptotic properties predicted by Corollary 3.1 using a manufactured solution. A subsequent example in Section 4.2 demonstrates a benefit of the unconditional stability property for the proposed algorithm. All computations were performed using the software package FreeFem++, [12]. 4.1. [0, 1] ⇥ [ n ˆ 2 = [0, function
Problem 1: Analytic Solution. Assume ⌦1 = [0, 1] ⇥ [0, 1] and ⌦2 = 1, 0], so I is the portion of the x-axis from 0 to 1. Then n ˆ 1 = [0, 1]T and T 1] . For a, ⌫1 , ⌫2 , and all arbitrary positive constants, the right hand side f1 is chosen to ensure that u1,1 (t, x, y) = ax2 (1
x)2 (1
y)e
u1,2 (t, x, y) = axy( 2 + y + 6x p1 (t, x, y) = e
t
t
,
3xy
4x2 + 2x2 y)e
t
,
cos(⇡x) sin(⇡y) .
One may verify that the above velocity field is divergence-free and that u1,2 is zero on I. The function f2 is chosen to ensure that 1) u2,1 satisfies the interface condition, 2) u2,2 is zero on I, 3) u2 is a divergence-free velocity field, and 4) p2 (t, x, y) = e t cos(⇡x) sin(⇡y). The spatial discretization is accomplished using the conforming P1 +bubble/P1 element pair, (the MINI-element, [1]). A convergence study is performed using the parameter values a = 1, ⌫1 = 0.5, ⌫2 = 0.05, = 100, T = 1. Table 4.1 gives the errors produced Algorithm 2.1. A compact notation is used to denote the various errors. Letting ui (t), pi (t) denote the NSE solution in ⌦i , i = 1, 2, approximated at t = tn by uni , pni respectively, the errors are Err(ui ) =
t
N X
n=0
Err(pi ) =
t
N X
n=0
kui (t )
uni k2H 1 (⌦i )
kpi (t )
pni k2L2 (⌦i )
n
n
!1/2
!1/2
.
The expected first-order convergence rates for the velocity are achieved. The pressure rates may be limited to first-order with the MINI-element; here the pressure convergence rates are slightly above the predicted value but decreasing. 19
t 1/4 1/8 1/16 1/32 1/64
h 1/4 1/8 1/16 1/32 1/64
Err(u1 ) 0.196536 0.071197 0.031561 0.015434 0.007684
rate 1.46 1.17 1.03 1.01
Err(u2 ) 1.673680 0.479848 0.160358 0.069638 0.032969
rate 1.80 1.58 1.20 1.08
Err(p1 ) 0.174935 0.034112 0.008706 0.002815 0.000978
rate 2.36 1.97 1.63 1.53
Err(p2 ) 0.172981 0.033365 0.008108 0.002893 0.001089
rate 2.37 2.04 1.49 1.41
Table 4.1 Errors for the manufactured solution test problem
4.2. Growth of error with decoupled methods. In addition to the GA scheme previously described, the following algorithms will be compared: Algorithm 4.1 (Fully Implicit Scheme). Let t > 0, f 2 H 1 (⌦). For each M 2 N, M Tt , given un 2 Vh , n = 0, 1, 2, · · · , M 1, find un+1 2 Vh satisfying ✓ n+1 ◆ X u un , v + ⌫ run+1 , v + ci (un+1 ; un+1 , vi ) i i t i=1,2 (4.1) Z n n+1 n+1 + |[u ]|[u ][v] = (f (t ), v) . I
Algorithm 4.2 (First-order IMEX Scheme). Let t > 0, f 2 H 1 (⌦). For each M 2 N, M Tt , given un 2 Vh , n = 0, 1, 2, · · · , M 1, find un+1 2 V1,h satisfying 1 ✓
un+1 1
un1 t
, v1
◆
⌦1
+ ⌫1 run+1 , rv1 1 +
Z
I
⌦1
+ c1 (un+1 ; un+1 , v1 ) 1 1
|[un ]|un+1 · v1 ds 1
= f1 (tn+1 ), v1
⌦1
Z
I
|[un ]|un2 · v1 ds
, 8v1 2 V1,h .
and un+1 2 V2,h satisfying 2 ✓ n+1 ◆ u2 un2 , v2 + ⌫2 run+1 , rv2 ⌦ + c2 (un+1 ; un+1 , v2 ) 2 2 2 2 t ⌦2 Z Z + |[un ]|un+1 · v ds |[un ]|un1 · v2 ds 2 2 I
(4.2)
(4.3)
I
= f2 (tn+1 ), v2
⌦2
, 8v2 2 V2,h .
The fully implicit scheme is fully coupled at the interface. A fully coupled calculation of this type is not currently used for a motivating application like climate modelling due to the expense of solving the large system of equations, and is used here only to generate reference data since it does not introduce the decoupling errors associated with the partitioned time stepping methods. The implicit-explicit (IMEX) scheme is, perhaps, a more intuitive method of partitioning the calculations than the GA method. However, Lemma 3.1 predicts energetic stability of the GA method regardless of the choice of problem parameters or discretization parameters. A similar result for the IMEX method is not known, and in fact this method may introduce energy numerically under some conditions. We proceed to investigate this e↵ect in the following example. Let ⌦1 = [0, 10] ⇥ [0, 1] and ⌦2 = [0, 10] ⇥ [ 1, 0], so I is the portion of the x-axis from 0 to 10, n ˆ 1 = [0, 1]T and n ˆ 2 = [0, 1]T . The kinematic viscosities are ⌫1 = 1 20
and ⌫2 = 0.2, with = 100. Consider the following functions: g1,1 (t, x, y) = 1.005 · x(10 g1,2 (t, x, y) = 1.005 · y(1 g2,1 (t, x, y) = x(10
x)(1
2y) exp( 3t) ,
y)(2x
10) exp( 3t) ,
x)(1 + 2y) exp( 2t) ,
g2,2 (t, x, y) = y(1 + y)(2x
10) exp( 2t) .
Dirichlet data is imposed by choosing ui,j = gi,j on @⌦i \I for i, j = 1, 2. A divergencefree initial condition is chosen such that ui · n ˆ i = 0 by setting ui,j (0, x, y) = gi,j (0, x, y) on ⌦i , for i, j = 1, 2. A mesh is constructed as in the previous section, consisting of triangles of horizontal width 10/32 and vertical width 1/32. The problem is run with a time step size of 1/40 on the time interval [0, 1]. The right hand side functions f1 and f2 are set to zero. Then Algorithm 4.1 is run. The implicit solution at the first time step and the initial data are used to start the two partitioned algorithms at the second time step. The partitioned algorithms are compared at time steps starting at t2 = 2 t. The test is constructed to demonstrate conditions under which the IMEX method will artificially introduce energy into the calculations due to the numerical treatment of the coupling terms. The fully coupled, first order implicit method, while not a practical scheme for the application, is useful as a point of comparison in the absence of an exact solution. In the first test we compute the kinetic energy of the velocity at the interface for each time step, as predicted by the 3 methods. That is, we plot 0.5 · kunj k2|L2 (I) , for j = 1, 2; there are two values as interface velocities have di↵erent values on ⌦1 and ⌦2 . These values are plotted in Figure 4.1 below. 1000 900 800
1000 Coupled IMEX method GA method
900 800 700
0.5 ku2 k2L2 (I)
0.5 ku1 k2L2 (I)
700 600 500 400
600 500 400
300
300
200
200
100
100
0 0
Coupled IMEX method GA method
0.5 time
1
0 0
0.5 time
1
Fig. 4.1. Relative solution sizes on I. Left: data on ⌦1 . Right: data on ⌦2 .
The damping terms exp( 2t) and exp( 3t) were included in the boundary conditions so that the true solution would decay to zero rapidly. Thus any growth in an approximate solution is an instability. (If these multipliers are not included we observed a slow, long term growth in the energy of the coupled solution.) The kinetic energy in the implicit methods decayed to zero rapidly. The energy in the GA scheme exhibited oscillations, decaying to zero over time. The energy in the IMEX scheme exhibited a similar behavior, but with a significantly higher kinetic energy. It is natural to investigate the growth of relative error for the partitioned methods, in light of the excess energy they exhibit. To distinguish the fully coupled data from 21
the partitioned methods, the solution to (4.1) on ⌦1 is denoted by un1,imp . Similarly, the solution to (4.1) on ⌦2 is denoted by un2,imp . Define the relative errors for the partitioned algorithms as RelErr(unj ) =
kunj unj,imp kH 1 (⌦j ) , for j = 1, 2. kunj,imp kH 1 (⌦j )
These errors are plotted in Figure 4.2. The relative errors for the IMEX method grow faster in both subdomains than those of the GA method. 8 7
3 IMEX method GA method
2.5
IMEX method GA method
6 2
RelErr(un2 )
RelErr(un1 )
5 4
1.5
3
1
2 0.5
1 0 0
0.5 time
1
0 0
0.5 time
1
Fig. 4.2. Relative solution errors. Left: errors on ⌦1 . Right: errors on ⌦2 .
5. Conclusions. For nonlinear interface conditions, stability of partitioned methods depends critically on their precise treatment. We have shown that, perhaps surprisingly, geometric averaging yields an unconditionally stable partitioned method while the most natural treatment may excessively introduce artificial energy. The kinetic energy and energy dissipation predicted by the model has been proven to deviate from values obtained using a fully implicit discretization by an amount proportional (asymptotically) to the time step size. Convergence as time and spatial discretization parameters go to zero was proven. REFERENCES [1] D. N. Arnold, F. Brezzi and M. Fortin, A stable finite element for the Stokes equations, Calcolo, Vol. 21, No. 4, 1984, pp. 337-344. [2] R.A. Adams and J.J.F. Fournier, Sobolev Spaces, vol. 140 of Pure and Applied Mathematics, Academic Press, 2003. [3] U. M. Ascher, S. J. Ruuth and B. T. R. Wetton, Implicit-explicit methods for time-dependent partial di↵erential equations, SIAM J. Num. Anal. 32(3), 1995. [4] C. Bernardi, T. Chac´ on-Rebollo, R. Lewandowski, F. Murat, A model of two coupled turbulent fluids, Part II: Numerical approximations by spectral discretization, SIAM Jour. Num. Analysis, Vol. 40, No. 6, pp. 2368-2394, 2003. [5] S.C. Brenner and L.R. Scott, The Mathematical Theory of Finite Element Methods, SpringerVerlag, 2002. [6] D. Bresch and J. Koko, Operator-Splitting and Lagrange Multiplier Domain Decomposition Methods for Numerical Simulation of Two Coupled Navier-Stokes Fluids, Int. J. Appl. Math. Comput. Sci., Vol. 16, No. 4, 2006, pp. 419-429. [7] E. Burman, Miguel A. Fern´ andez, Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility, INRIA Research Rept. RR-6455, Feb. 2008. 22
[8] P. Causin, J.-F. Gerbeau and F. Nobile, Added-mass e↵ect in the design of partitioned algorithms for fluid-structure problems, INRIA Rept. 5084, 2004. [9] J. Connors, J. Howell and W. Layton, Partitioned timestepping for a parabolic two domain problem, to appear, SIAM Jour. Num. Analysis, 2009. [10] G.T. Cooper and A. Sayfy, Additive Runge-Kutta methods for sti↵ ordinary di↵erential equations, Mathematics of Computation, Vol. 40, No. 161, pp. 207-218. [11] G. P. Galdi, An Introduction to the Mathematical Theory of the Navier-Stokes Equations, Springer Tracts in Natural Philosophy, Volume I, Springer-Verlag, New York, 1994. [12] F. Hecht, A. LeHyaric and O. Pironneau. Freefem++ version 2.24-1, 2008. http://www.freefem.org/ff++. [13] R. Ingram, Unconditional convergence of high-order extrapolations of the Crank-Nicolson, nite element method for the Navier-Stokes equations, University of Pittsburgh technical report TR-MATH-10-09, 2010. [14] C.A. Kennedy and M. H. Carpenter, Additive Runge-Kutta schemes for convection-di↵usionreaction equations, Applied Numerical Math., Vol. 44, pp. 139-181, 2003. [15] W. Layton, Introduction to the Numerical Analysis of Incompressible Viscous Flows, SIAM Computational Science and Engineering 6, 2008. [16] R. Lewandowski, Analyse Math´ ematique et Oc´ eanographic, collection RMA, Masson, 1997. [17] J. -L. Lions, R. Temam and S. Wang, Models for the coupled atmosphere and ocean (CAO I), Computational Mechanics Advances 1, 1993, pp. 1-54. [18] J. -L. Lions, R. Temam and S. Wang, Numerical analysis of the coupled atmosphere-ocean models (CAO II), Computational Mechanics Advances 1, 1993, pp. 55-119. [19] J. -L. Lions, R. Temam and S. Wang, Mathematical theory for the coupled atmosphere-ocean models (CAO III), J. Math. Pures. Appl. (74), 1995, pp. 105-163. [20] P. M¨ uller, The equations of oceanic motions, Cambridge University Press, Cambridge, UK, 2006. [21] W. M. Washington and C. L. Parkinson, An introduction to three-dimensional climate modeling, second edition, University Science Books, California, 2005.
23
