On error estimates of the penalty method for the viscoelastic flow problem I: Time discretization

Applied Mathematical Modelling 34 (2010) 4089–4105

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

On error estimates of the penalty method for the viscoelastic flow problem I: Time discretization q Kun Wang a, Yinnian He b,a,*, Xinlong Feng b a b

Faculty of Science, Xi’an Jiaotong University, Xi’an 710049, PR China College of Mathematics and Systems Science, Xinjiang University, Urumqi 830046, PR China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 30 November 2008 Received in revised form 8 April 2010 Accepted 22 April 2010 Available online 29 April 2010

In this paper, we investigate the two-dimensional viscoelastic fluid motion problem arising in the Oldroyd model. By applying the dual method, the Helmholtz decomposition and some other techniques, we deduce the long-time optimal error estimates for its penalty method. Furthermore, the global optimal error bounds in the L2 and H1-norms for the time semi-discrete of the penalty system are derived under the backward Euler implicit scheme, which improves the best estimate available. Ó 2010 Elsevier Inc. All rights reserved.

Keywords: Viscoelastic flow problem Oldroyd model Penalty method Time discretization Global error estimates

1. Introduction The viscoelastic Oldroyd fluid motion problem is governed by the rheological relation (see [1–3])

k0 r þ k1

@r @f ¼ g0 f þ g1 ; @t @t

k1 rðx; 0Þ ¼ g1 fðx; 0Þ:

Here r is the stress tensor, k0, k1, g0, g1 are positive constants and f is the strain tensor with components fij ¼ 12 i; j ¼ 1; . . . ; n, where u = u(x, t) = (u1(x, t), . . . , un(x, t)) is the velocity of the fluid motion and n = 2 or 3. The Cauchy form of the above rheological relation is the following initial-boundary value problem (see [4]):

@u  mDu  @t

Z



@ui @xj

 @u þ @xj ; i

t

bðt  sÞDuds þ ðu  rÞu þ rp ¼ f ;

div u ¼ 0ðx; tÞ 2 X  Rþ ;

ð1Þ

0

uðx; tÞj@ X ¼ 0 t 2 Rþ ;

uðx; 0Þ ¼ u0 ðxÞ x 2 X;

ð2Þ

where

bðtÞ ¼ qedt ;

m¼

g1 k1

;

q¼

g0 k1  k0 g1 2

k1

;

d¼

k0 ; k1

q Supported by the National Natural Science Foundation of China(Nos. 10971166, 10901131), the National Basic Research Program (No. 2005CB321703), the China Scholarship Council (No. 2009628086) and the Natural Science Foundation of Xinjiang Province (No. 2010211B04). * Corresponding author. Address: Faculty of Science, Xi’an Jiaotong University, Xi’an 710049, PR China. Tel./fax: +86 29 82665242. E-mail addresses: [email protected] (K. Wang), [email protected] (Y. He), [email protected] (X. Feng).

0307-904X/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2010.04.008

4090

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

q P 0 is the viscoelastic coefficient, 1/d is the relaxation time, X is an open bounded domain in Rn with smooth boundary oX, p = p(x, t) is the pressure of the fluid, f = f(x, t) is the prescribed external force and u0 = u0(x) is the initial velocity. If

g0k1 = k0g1, the model reduces to Newton’s model of incompressible viscous fluid which is widely known as the Navier– Stokes equations. The system (1) and (2) is based on a structural model for polymeric fluids, suspensions or biological fluids. It was derived under the assumptions that the material can be regarded as a single stationary macroscopic element with small stress and strain rates. For further details of the physical background and its mathematical modeling, the reader is referred to [1–3,5,6]. Eqs. (1) and (2) have been widely investigated. In [7–10], the authors dealt with the existence, uniqueness and continuous dependence of the solution upon the data. Its asymptotic behavior was also analyzed in [11,12]. Moreover, there are lots of works devoted to the numerical approximations of the problem (1) and (2). For the spatial discretization, Akhmatov and Oskolkov [13] considered the difference schemes, Cannon et al. [10] analyzed a modified nonlinear Galerkin method, He et al. [14] and Pani and Yuan [15] investigated the conforming finite element method. For the time discretization, by using the semigroup theory, Pani et al. considered a linearized backward Euler scheme in [16]. Recently, Wang et al. [17] extended the analysis to a fully discrete finite element method. Besides, the linearized viscoelastic Oldroyd flow problem is also considered in [18]. In this article, we consider the problem (1) and (2) in R2 . In the above equations, the velocity u and the pressure p are coupled together by the incompressibility condition ‘‘divu = 0”, which makes the system difficult to solve by using the numerical methods. In order to overcome this difficulty, the penalty method is popular used (see [19–24], and the references therein). The penalty method for the viscoelastic fluid motion problem (1) and (2) is as follows:

@ue  mDue  @t

Z 0

t

e m

divue þ pe ¼ 0ðx; tÞ 2 X  Rþ ;

bðt  sÞDue ds þ ðue  rÞue þ rpe ¼ f ; ue ðx; tÞj@ X ¼ 0t 2 Rþ :

ue ðx; 0Þ ¼ u0 ðxÞx 2 X;

ð3Þ ð4Þ

The penalty method firstly was introduced by Courant (see [25]). In the 70’s of last century, Temam extended it to the Navier–Stokes equations in [24]. Then, lots of works appeared on this subject (see [20,23], and the references therein). Shen derived the optimal error estimates for its penalty system in [23] as follows:

s1=2 ðtn Þkuðtn Þ  ue ðtn ÞkL2 þ sðtn Þkuðtn Þ  ue ðtn ÞkH1 6 ce;

ð5Þ

for tn 2 [0, T] with T being a finite time, where s(tn) = min{1, tn}, c is a general positive constant and u(tn), ue(tn) are the solutions of the Navier–Stokes equations and its penalty system, respectively. Recently, this idea is applied to the finite element method by He [20]. For the viscoelastic flow problem, the literatures are relatively limited. In [21,22], Kotsiolis and Oskolkov proved the existence, uniqueness of the solution of the problem (3) and (4) and lime?0(ue(t), pe(t)) = (u(t), p(t)). They also gave the following error estimate in [22] for all tn P 0:

kuðt n Þ  ue ðt n ÞkL2 þ

Z

tn

0

kuðtÞ  ue ðtÞk2H1 dt

1=2

6c

pffiffiffi

e:

ð6Þ

Compared with the estimate (5) for the unsteady Navier–Stokes equations, (6) is not optimal. Furthermore, when considering the discrete problem for (3) and (4), the estimate (6) is misleading. For example, if the backward Euler scheme (see (107) in Section 5) is applied to (3) and (4), the estimate (6) would lead to n

kuðt n Þ  ue kL2 þ Dt

n X

!1=2 kuðt m Þ 

2 um e kH 1

6 cðDt þ

pffiffiffi

eÞ

ð7Þ

m¼0

for all tn P 0, where 0 < Dt < 1 is the time step size, tn = nDt and une is a penalty approximation of u at the time tn. (7) suggests the choice e = Dt2, which would result in a very ill-conditioned system when we make further spatial discretization (see [23]). The main focus of this paper is to apply the techniques in [23] to the viscoelastic fluid motion equations and derive the optimal error estimates for the penalty system. Moreover, we extend the analysis to the long time behavior and prove that the error estimates are valid uniform in time. If  = m + q/d and f1(x) = limt?1f(x, t) satisfy the following condition:

2Nm1 1 kf1 k1 < 1; where N ¼

sup u;v ;w2X;u;v ;w–0

jbðu;v ;wÞj ; kf1 k1 kukkv kkwk

ð8Þ ¼ sup

v 2X;v –0

jðf1 ;v Þj ;X kv k

¼

ðH10 ð

2

XÞÞ with the norm kuk ¼ krukL2 , and the trilinear form

b(u, v, w) is defined in Section 2, we have the following error estimates for all tn P 0:

s1=2 ðtn Þkuðtn Þ  ue ðtn ÞkL2 þ sðtn Þkuðtn Þ  ue ðtn ÞkH1 6 ce; s1=2 ðtn Þkuðtn Þ  une kL2 þ sðtn Þkuðtn Þ  une kH1 6 cðDt þ eÞ

ð9Þ ð10Þ

for sufficiently small e and Dt, which substantially improves the previous results (6) and (7). In particular, if e = 0, our error estimate based on the backward Euler scheme in H1-norm is sðtn Þkuðt n Þ  un kH1 6 cDt. However, the best error estimate available, to our knowledge, is (see [16])

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

  1 : kuðt n Þ  un kH1 6 cDt t n1=2 þ log Dt

4091

ð11Þ

The remainder of the paper is organized as follows. In the next section, we introduce some notations and preliminary results for the viscoelastic flow problem (1) and (2). We provide error estimates for the linearized penalty system in Section 3 and for the nonlinear penalty system in Section 4. In Section 5, we analyze the backward Euler time discretization scheme for the nonlinear penalty system. Finally, conclusions are given in Section 6. 2. Preliminaries Firstly, we introduce the mathematical setting for the problem (1) and (2). Let X be an open bounded domain in R2 with boundary oX 2 C2 and satisfy the additional condition stated in (A1) below. We introduce the following Hilbert spaces:

X ¼ ðH10 ðXÞÞ2 ;

Y ¼ ðL2 ðXÞÞ2 ;

M ¼ L20 ðXÞ ¼

  Z q 2 L2 ðXÞ; qdx ¼ 0 X

and the Laplace operator

Au ¼ Du 8u 2 DðAÞ ¼ ðH2 ðXÞÞ2 \ X: We denote by kki the usual norm of the Sobolev space Hi(X) or (Hi(X))2 for i = 1, 2, and by (, ) and jj the inner product and norm on L2(X) or (L2(X))2 (see [26] for more details). The spaces H10 ðXÞ and X are equipped with their usual scalar product and norm

ððu; v ÞÞ ¼ ðru; rv Þ;

kuk ¼ ððu; uÞÞ1=2 :

The frequently used Hilbert spaces V and H are defined by

V ¼ fv 2 X; div

v ¼ 0g;

H ¼ fv 2 Y; div

v ¼ 0; v  nj@X ¼ 0g:

Let Ae u ¼ Du  1e rdivu and the Stokes operator be determined by A ¼ PD with DðAÞ ¼ ðH2 ðXÞÞ2 \ V, where P denotes the L2-orthogonal projection of Y onto H. The power of Aa ; Aae and Aa ða 2 RÞ are well defined. It follows that e 1=2 u; A e 1=2 v Þ 8u 2 DðAÞ; v 2 X, where A e ¼ A or Ae, and ðAu; v Þ ¼ ðA1=2 u; A1=2 v Þ 8u 2 DðAÞ; v 2 V. In particular, there e v Þ ¼ ðA ð Au; holds

    1 1=2 1=2 1=2 A1=2 e u; Ae v ¼ A u; A v þ ðdiv u; div v Þ 8u;

e

v 2 X:

ð12Þ

As mentioned above, we need the additional assumption on X: (A1) Assume that X is regular in the sense that a unique solution (v, q) 2 (X, M) of the Stokes problem

mDv þ rq ¼ g; div

v ¼ 0 in X; v j@X ¼ 0

for any prescribed g 2 Y exists and satisfies

kv k2 þ kqk1 6 cjgj; where c > 0 is a positive constant depending on X and m. In the following, c is a general positive constant independent of e and Dt, which may stand for different values at different occurrences. It is valid that (see [26,27])

jAv j2 6 kv k22 6 cjAv j2 k1 jv j2 6 kv k2

8v 2 X;

8v 2 ðH2 ðXÞÞ2 \ V; k1 kv k2 6 kv k22 ;

ð13Þ kv k22 6 cjAv j 8v 2 DðAÞ;

where k1 is the minimal eigenvalue of the Laplace operator D. Furthermore, we make the following assumption about the prescribed data for the problem (1) and (2): (A2) The initial velocity u0 2 V and the body force f(x, t) satisfy f(x, t), ft(x, t) 2 L1(R+;Y) with

ku0 k2 þ supðjf ðx; tÞj2 þ jft ðx; tÞj2 Þ 6 c: tP0

We define the continuous bilinear forms a(, ) on X  X and d(, ) on X  M, respectively, by

aðu; v Þ ¼ mððu; v ÞÞ 8u;

v 2 X;

dðv ; qÞ ¼ ðq; div v Þ 8v 2 X; q 2 M:

ð14Þ

4092

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

With

1 Bðu; v Þ ¼ ðu  rÞv þ ðdivuÞv ; 2 the trilinear form b(, , ) on X  X  X is determined by

1 1 1 bðu; v ; wÞ ¼ ðBðu; v Þ; wÞ ¼ ððu  rÞv þ ðdivuÞv ; wÞ ¼ ððu  rÞv ; wÞ  ððu  rÞw; v Þ 8u; v ; w 2 X: 2 2 2 It follows that (see [23,28–30])

bðu; v ; wÞ ¼ bðu; w; v Þ 8 u; v ; w 2 X;

ð15Þ

jbðu; v ; wÞj 6 cðkukkv k1=2 jAv j1=2 jwj þ juj1=2 jAuj1=2 kv kjwjÞ 8 u; v ; w 2 X;

ð16Þ

jbðu; v ; wÞj 6 Nkukkv kkwk 8 u; v ; w 2 X;

ð17Þ

jbðu; v ; wÞj 6 cjuj1=2 kuk1=2 ðkv kjwj1=2 kwk1=2 þ kwkjv j1=2 kv k1=2 Þ 8 u; v ; w 2 X:

ð18Þ

With above notations, the variational formulation of the problem (1) and (2) is defined as follows: Find (u, p) 2 (X, M) for all t P 0, such that for all (v, q) 2 (X, M):

ðut ; v Þ þ aðu; v Þ þ Jðt; u; v Þ þ bðu; u; v Þ  dðv ; pÞ þ dðu; qÞ ¼ ðf ; v Þ

ð19Þ

and the variational formulation of the penalty system (3) and (4) reads as: find (ue, pe) 2 (X, M) for all t P 0 such that for all (v, q) 2 (X, M):

e m

ð20Þ

ðuet ; v Þ þ aðue ; v Þ þ Jðt; ue ; v Þ þ bðue ; ue ; v Þ  ðrdivue ; v Þ ¼ ðf ; v Þ

ð21Þ

ðuet ; v Þ þ aðue ; v Þ þ Jðt; ue ; v Þ þ bðue ; ue ; v Þ  dðv ; pe Þ þ dðue ; qÞ þ ðpe ; qÞ ¼ ðf ; v Þ; or find ue 2 X for all t P 0, such that for all

v 2 X: m e

with u(0) = ue(0) = u0, where

    Z t Z t : Jðt; u; v Þ ¼ q edt eds AuðsÞds; v ¼ q edt eds uðsÞds; v 0

0

Before continuing the numerical analysis below, we need to recall the following results (see [14]). Lemma 2.1. Assume that s > 0 and u 2 L1(0, s; X). It is valid that

Z

s

Jðt; u; e

2d0 t

0

Z 2 Z s 2 Z s t 1 2a0 s ds 2a0 t ds uðtÞÞ dt ¼ qe e uðsÞds þ qa0 e e uðsÞds dt: 2 0 0 0

ð22Þ

Moreover, if u, v 2 L1(0, s; D(A)), it is valid that

Z

s

Jðt; u; e

0

2d0 t

Z

2

Z

2 Z s

t



1 2a0 s

s ds

2a0 t

ds

AuðtÞÞ dt ¼ qe e AuðsÞds þ qa0 e e AuðsÞds dt;





2 0 0 0

ð23Þ

where 0 < d0 < 12 minfd; mk1 g; a0 ¼ d  d0 . Lemma 2.2. Assume that s > 0 and u 2 L2(0, s; X). It is valid that

Z

a1 0

s

0

2 Z s Z t 2 Z s ds 2 eds uðsÞds e2a0 t e2d0 t kuðtÞk2 dt: þ e uðsÞds dt 6 a0 0

0

ð24Þ

0

Moreover, if u 2 L2(0, s; D(A)), it is valid that

Z



a1 0

s 0

Z

2

2 Z s Z s

t





eds AuðsÞds

þ e2a0 t

eds AuðsÞds dt 6 a2 e2d0 t jAuðtÞj2 dt: 0

0

0 0

ð25Þ

Lemma 2.3. Assume that s > 0, k P 0, sk/2(t)u 2 L1(0, s; X) and ut 2 L2(0, s; Y). It is valid that

Z

2

0

s





Z

m q 2q sk ðtÞJðt; u; e2d0 t ut ðtÞÞ dt

6 sk ðsÞked0 s uðsÞk2 þ 2 a þdþkþ a0 0 8 m

0

s

e2d0 t kuk2 dt:

ð26Þ

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

4093

Moreover, if sk/2(t)u 2 L1(0, s; D(A)), and ut 2 L2(0, s; Y), it is valid that

Z

2

0

s





Z

m q 2q sk ðtÞJðt; u; e2d0 t Aut ðtÞÞ dt

6 sk ðsÞjed0 s AuðsÞj2 þ 2 a þdþkþ a0 0 8 m

s

e2d0 t jAuj2 dt:

ð27Þ

0

Theorem 2.1. Suppose that assumptions (A1)–(A2) are valid. Then the solution (u, p) of the problem (1) and (2) satisfies the following regularities, for all s P 0,

kuðsÞk2 þ e2d0 s

Z

s

e2d0 t ðjAuj2 þ kpk21 þ jut j2 Þ dt 6 c; Z s sðsÞðjAuðsÞj2 þ jut ðsÞj2 þ kpðsÞk21 Þ þ e2d0 s e2d0 t sðtÞkut k2 dt 6 c; 0 Z s 2 2 2 2d0 s 2d0 t 2 s ðsÞkut ðsÞk þ e e s ðtÞðjAut j þ kpt k21 Þ dt 6 c:

ð28Þ

0

ð29Þ ð30Þ

0

Besides, the Gronwall lemmas will be frequently used too. Lemma 2.4. Gronwall lemma [10,19]Let g, h, y be three locally integrable nonnegative functions on the time interval [t0, 1) that for all t P t0, satisfy

yðtÞ þ GðtÞ 6 C þ

Z

t

hðsÞds þ

t0

Z

t

gðsÞyðsÞds;

t0

where G (t) is a nonnegative function on [0, 1), C P 0 is a constant. Then,

  Z t  Z t yðtÞ þ GðtÞ 6 C þ hðsÞds exp gðsÞds : t0

ð31Þ

t0

Lemma 2.5. discrete Gronwall lemma [30]Let k, C and am, bm, cm, rm, for integers m P 0, be nonnegative numbers such that

an þ k

n X

bm 6 k

m¼0

n X

r m am þ k

m¼0

n X

8n P 0:

cm þ C

m¼0

Suppose that krm < 1 for all m, and set rm = (1  krm)1, then,

an þ k

n X

n X

bm 6

m¼0

!

cm þ C exp k

m¼0

n X

!

rm r m

8n P 0:

ð32Þ

m¼0

Finally, for the operator Aeu associated with the penalty method, we recall the following lemma given in [19,23]. Lemma 2.6. There exists a positive constant c0 > 0 depending only on X and such that if ec0 6 1, we have

jA1=2 uj 6 c0 jA1=2 e uj 8u 2 X;

ð33Þ

jAuj 6 c0 jAe uj 8u 2 DðAÞ:

ð34Þ

3. Error estimates for the linearized problem As an intermediate step, in this section, we will investigate the linearized viscoelastic flow problem:

@u  mDu  @t

Z

t

bðt  sÞDuds þ rp ¼ f ;

div u ¼ 0ðx; tÞ 2 X  Rþ ;

ð35Þ

0

uðx; 0Þ ¼ u0 ðxÞ x 2 X;

uðx; tÞj@ X ¼ 0 t 2 Rþ :

ð36Þ

The penalty method of (35) and (36) is as follows:

@ue  mDue  @t

Z

t 0

bðt  sÞDue ds þ rpe ¼ f ;

ue ðx; 0Þ ¼ u0 ðxÞ x 2 X;

e m

div ue þ pe ¼ 0ðx; tÞ 2 X  Rþ ;

ue ðx; tÞj@X ¼ 0 t 2 Rþ :

ð37Þ ð38Þ

Similarly, the variational formulation of the problem (35) and (36) is defined as follows: find (u, p) 2 (X, M) for all t P 0, such that for all (v, q) 2 (X, M):

4094

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

ðut ; v Þ þ aðu; v Þ þ Jðt; u; v Þ  dðv ; pÞ þ dðu; qÞ ¼ ðf ; v Þ

ð39Þ

and the variational formulation of the penalty system (37) and (38) reads as: find (ue, pe) 2 (X, M) for all t P 0 such that for all (v, q) 2 (X, M):

e m

ðuet ; v Þ þ aðue ; v Þ þ Jðt; ue ; v Þ  dðv ; pe Þ þ dðue ; qÞ þ ðpe ; qÞ ¼ ðf ; v Þ;

v 2 X: m ðuet ; v Þ þ aðue ; v Þ þ Jðt; ue ; v Þ  ðrdivue ; v Þ ¼ ðf ; v Þ e

ð40Þ

or find ue 2 X for all t P 0, such that for all

ð41Þ

with u(0) = ue(0) = u0. Setting e = u  ue, l = p  pe, it follows that e(0) = 0. Next, we estimate e and l, which will be used as intermediate results for analyzing the nonlinear problem in Section 4. To get the optimal error estimate for ju  uej, we begin with the parabolic duality argument (see [23]). For any s > 0 and g ¼ e2d0 t e 2 L2 ð0; s; L2 ðXÞ2 Þ, we consider the problem: Find (U(t), W(t)) 2 (X, M) such that for 0 < t < s,

ðv ; Ut Þ  aðv ; UÞ  dðv ; WÞ þ dðU; qÞ  Jðt; v ; UÞ ¼ ðv ; e2d0 t eÞ

ð42Þ

for all (v, q) 2 (X, M) with U(s) = 0. Theorem 2.1 implies that u is sufficiently smooth so that (U, W) is correctly defined for all 0 < t 6 s. Thus (42) is a wellposed problem and has a unique solution (U, W), with

U 2 Cð0; s; VÞ \ L2 ðr; s; H2 ðXÞ2 \ VÞ \ H1 ðr; s; L2 ðXÞ2 Þ; W 2 L2 ðr; s; H1 ðXÞ \ MÞ 8r 2 ð0; sÞ: Lemma 3.1. Let (U, W) be the solution of the problem (42) with the body force function g ¼ e2d0 t e ¼ e2d0 t ðu  ue Þ, then, for any s P 0, it is valid that

sup e2d0 t kUðtÞk21 þ

Z

06t6s

s

0

e2d0 t ðkUk22 þ jUt j2 þ kWk21 Þ dt 6 c

Proof. Applying P on (42), setting

m

e2d0 t jA1=2 Uj2 P

2

jðe; UÞj 6

m 4

mk1 2

Z

s

e2d0 t jej2 dt:

ð43Þ

0

v ¼ e2d t U and considering that 0

e2d0 t jUj2 P d0 e2d0 t jUj2 ;

e2d0 t jA1=2 Uj2 þ ce2d0 t jej2 ;

we have



d 2d0 t 2 m 2d0 t 1=2 2 e jUj þ e jA Uj þ 2Jðt; e2d0 t U; UÞ 6 ce2d0 t jej2 : dt 2

Integrating the above inequality with respect to the time t from r P 0 to s, using Lemma 2.1, we obtain

sup e2d0 r jUðrÞj2 þ 06r6s

Taking

Z

s

e2d0 t jA1=2 Uj2 dt 6 c

Z

0

s

e2d0 t jej2 dt:

ð44Þ

0

v ¼ AU in (42), we get

d  jA1=2 Uj2 þ mjAUj2 þ 2jJðt; AU; UÞj 6 ce4d0 t jej2 : dt

ð45Þ

Multiplying (45) by e2d0 t , integrating with respect to the time from r to s and using Lemma 2.1 and (44), we deduce that

sup e2d0 r jA1=2 UðrÞj2 þ

Z

06r6s

s

e2d0 t jAUj2 dt 6 c

0

Applying P on (42), taking

Z

s

e2d0 t jej2 dt:

0

v ¼ e2d t Ut 0

and making the following estimate

Z t

2



1 jJðt; Ut ; UÞj 6 jUt j2 þ ce2dt

eds AUds

; 4 0

it follows:

e2d0 t jUt j2 

Z t

2



d  2d0 t 1=2 2  e mjA Uj 6 ce2d0 t jej2 þ 2d0 me2d0 t jA1=2 Uj2 þ ce2ðdþd0 Þt

eds AUds

: dt 0

ð46Þ

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

4095

Integrating the above inequality from r to s and using Lemma 2.2, (44) and (46), there holds

sup e2d0 r kUðrÞk21 þ

Z

06r6s

s

Z

e2d0 t jUt j2 dt 6 c

0

s

e2d0 t jej2 dt:

ð47Þ

0

On the other hand, (42) can be rewritten as

Ut þ mDU þ edt

Z

t

eds DUds þ rW ¼ e2d0 t e:

ð48Þ

0

Therefore, combining (48) with (46) and (47), and using Lemma 2.2, we have

Z

s

e2d0 t jrWj2 dt 6 c

Z

0

s

e2d0 t jej2 dt:

ð49Þ

0

Combining (49) with (46) and (47), we complete the proof.

h

Lemma 3.2. Under the assumptions of (A1) and (A2), it is valid, for any s P 0, that

jeðsÞj2 þ e2d0 s

Z

s

e2d0 t kek2 dt þ ee2d0 s

0

e2d0 s

Z

s

Z

s

e2d0 t jlj2 dt 6 ce;

ð50Þ

0

e2d0 t jej2 dt 6 ce2 :

ð51Þ

0

Proof. Subtracting (40) from (39), we obtain

e m

e m

ðet ; v Þ þ aðe; v Þ þ Jðt; e; v Þ  dðv ; lÞ þ dðe; qÞ þ ðl; qÞ ¼ ðp; qÞ:

ð52Þ

Taking ðv ; qÞ ¼ e2d0 t ðe; lÞ in (52), using Cauchy inequality and noting the fact that

1 mk mkek2 P 1 jej2 P d0 jej2 ; 2 2 we deduce that

1 d 2d0 t 2 m e 2e ðe jej Þ þ e2d0 t kek2 þ Jðt; e; e2d0 t eÞ þ e2d0 t jlj2 6 e2d0 t jpj2 : 2 dt 2 2m m

ð53Þ

Integrating (53) with respect to the time from 0 to s and using Lemma 2.1, we get

e2d0 s jeðsÞj2 þ m

Z

s

e2d0 t kek2 dt þ e

0

Z

s

e2d0 t jlj2 dt 6 ce

0

Z

s

e2d0 t jpj2 dt:

0

Multiplying the above inequality by e2d0 s and applying (28), we arrive at

jeðsÞj2 þ e2d0 s

Z

s

e2d0 t kek2 dt þ ee2d0 s

Z

s

e2d0 t jlj2 dt 6 ce;

ð54Þ

e2d0 t jej2 ¼ ðUt ; eÞ  aðe; UÞ  dðe; WÞ þ dðU; lÞ  Jðt; e; UÞ:

ð55Þ

0

0

which yields (50). Setting (v, q) = (e, l) in (42), it follows that:

Taking (v, q) = (U, W) in (52), then adding it to (55), we obtain

e2d0 t jej2 ¼

d e ðU; eÞ  ðW; pe Þ: dt m

ð56Þ

Integrating (56) from 0 to s and using (14), we have

Z 0

s

Z s 1=2 Z s 1=2 e2d0 t jej2 dt 6 ce e2d0 t kWk21 dt e2d0 t jpe j2 dt : 0

0

Applying Lemma 3.1 to the above inequality, then multiplying it by ed0 s , using the triangle inequality, (28) and (50), we deduce that

e2d0 s

Z 0

s

e2d0 t jej2 dt 6 ce2 :

4096

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

Combining the above inequality with (54), we complete the proof. h Theorem 3.1. Suppose (A1) and (A2) are valid, e ? 0 such that ec0 6 1, where c0 is a positive constant. Then, for any s P 0, we have

sðsÞjeðsÞj2 þ e2d0 s

Z

s

e2d0 t sðtÞkek2 dt þ ee2d0 s

0

s2 ðsÞkeðsÞk2 þ e2d0 s

Z

s

e2d0 t sðtÞjlj2 dt 6 ce2 ;

ð57Þ

0

Z

s

e2d0 t s2 ðtÞjlj2 dt 6 ce2 :

ð58Þ

0

Proof. Let us consider the Helmholtz decomposition (see [29])

H10 ¼ V  V ? ; where V\ = {(D)1rq:q 2 L2(X)} and u(t) 2 V\, such that divu(t) = p(t) with

v = (D)1rq

iff Du = rq and vjoX = 0. Thus, for p(t) 2 M, there exists a unique

kuðtÞk 6 cjpðtÞj 8t P 0:

ð59Þ

Moreover, if pt(t) 2 M, there holds divut(t) = pt(t) with

kut ðtÞk 6 cjpt ðtÞj 8t P 0: Taking ðv ; qÞ ¼ e

2d0 t

ð60Þ

sðtÞðe; lÞ in (52), we have

1 d 2d0 t e e ðe sðtÞjej2 Þ þ me2d0 t sðtÞkek2 þ e2d0 t sðtÞjlj2 þ Jðt; e; e2d0 t sðtÞeÞ 6 ce2d0 t jej2 þ e2d0 t sðtÞðp; lÞ 2 dt m m

e m

e m

¼ ce2d0 t jej2 þ e2d0 t sðtÞðdiv u; lÞ ¼ ce2d0 t jej2  e2d0 t sðtÞðrl; uÞ:

ð61Þ

Since

rl ¼ mDe þ q

Z

t

edðtsÞ Deds  et ;

ð62Þ

0

we arrive at

1 d 2d0 t e ðe sðtÞjej2 Þ þ me2d0 t sðtÞkek2 þ e2d0 t sðtÞjlj2 þ Jðt; e; e2d0 t sðtÞeÞ 2 dt m  Z t  edðtsÞ eds; Au 6 ce2d0 t jej2 þ ee2d0 t sðtÞðet ; uÞ þ mee2d0 t ððe; uÞÞ þ ee2d0 t sðtÞ q 0

d  2d0 t ¼ ce jej þ e e sðtÞðe; uÞ  ee2d0 t sðtÞðe; ut Þ  ee2d0 t ðe; uÞ  2d0 esðtÞðe; uÞ þ mee2d0 t sðtÞððe; uÞÞ dt  Z t  2d0 t edðtsÞ eds; Au : þ ee sðtÞ q 2d0 t

2

ð63Þ

0

Applying (14), (59) and (60) and Theorem 2.1, we obtain

1 2 1 jej þ ce2 juj2 6 jej2 þ ce2 jpj2 ; 4 4 1 2 1 2 jesðtÞðe; uÞj 6 jej þ ce sðtÞjuj2 6 jej2 þ ce2 sðtÞjpj2 ; 4 4 1 2 1 2 2 2 jesðtÞðe; ut Þj 6 jej þ ce s ðtÞjut j 6 jej2 þ ce2 s2 ðtÞjpt j2 ; 4 4 1 1 2 2 mesðtÞjððe; uÞÞj 6 msðtÞkek þ cme sðtÞkuk2 6 msðtÞkek2 þ ce2 sðtÞjpj2 ; 2 2

Z t

2

 Z t 





dt d s 2dt d s ee sðtÞ

q e eds; Au

6 ce

e eds

þ e2 sðtÞkA1=2 uk2 0 0

Z t

2



6 ce2dt

eds eds

þ ce2 sðtÞkpk21 : jceðe; uÞj 6

0

Integrating (63) from 0 to s, using the above estimates, Lemmas 2.1, 2.2, 3.2, (28), (30) and multiplying by e2d0 s , we get

sðsÞjeðsÞj2 þ e2d0 s

Z 0

s

e m

e2d0 t sðtÞðkek2 þ jlj2 Þ dt 6 ce2 þ ce2d0 s

Z 0

s

e2d0 t jej2 dt þ ce2 e2d0 s

Z 0

t

e2d0 t s2 ðtÞjpt j2 ds 6 ce2 :

ð64Þ

4097

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

Differentiating the terms dðe; qÞ þ me ðl; qÞ  me ðp; qÞ with respect to the time t in (52) and taking ðv ; qÞ ¼ e2d0 t s2 ðtÞðet ; lÞ, we have

e2d0 t s2 ðtÞjet j2 þ

 1 d  2d0 t 2 e e s ðtÞðmkek2 þ jlj2 Þ þ Jðt; e; e2d0 t s2 ðtÞet Þ 2 dt m

e m

e m

6 ce2d0 t sðtÞkek2 þ ce2d0 t sðtÞ jlj2 þ e2d0 t s2 ðtÞðpt ; lÞ:

ð65Þ

Considering (14) and (62), it is valid that

 Z

e 2 e e e s ðtÞðpt ; lÞ ¼ s2 ðtÞðdivut ; lÞ ¼  s2 ðtÞðut ; rlÞ ¼ s2 ðtÞ ðet ; ut Þ þ ððe; ut ÞÞ þ q m m m m

0

t

edðtsÞ eds; Aut



2

Z t



1 6 s2 ðtÞjet j2 þ s2 ðtÞkek2 þ ce2 s2 ðtÞkut k2 þ ce2 s2 ðtÞkA1=2 ut k2 þ ce2dt

eds eds

2 0

Z t

2



1 6 s2 ðtÞjet j2 þ s2 ðtÞkek2 þ ce2 s2 ðtÞkpt k21 þ ce2dt

eds eds

: 2 0 Integrating (65) from 0 to s and using the above inequality and Lemmas 2.2, 2.3, 3.2, (30) and (64), we get

e s2 ðsÞkeðsÞk2 þ s2 ðsÞjlðsÞj2 þ e2d0 s m

Z

s

e2d0 t s2 ðtÞjet j2 dt 6 ce2 ;

ð66Þ

0

by a final multiplication of e2d0 s . Finally, due to (62), there holds

jlj2 6 ckrlk21 6 c



mkDek21 þ ket k21 þ kq

Z 0

t

edðtsÞ Dedsk21



  Z t 6 c kek2 þ jet j2 þ qe2dt k eds edsk2 :

ð67Þ

0

Therefore, using Lemmas 2.2, 3.2, (64) and (66), we have

e2d0 s

Z

s

e2d0 t s2 ðtÞjlj2 dt 6 ce2 :

0

Combining (64), (66) with the above inequality, we complete the proof.

h

4. Error estimates for the nonlinear problem The focus of this section is to deduce the global optimal error estimates for the penalty system of the nonlinear viscoelastic flow problem. Preceding to give the error estimate, we firstly derive some regularity results. Lemma 4.1. Under the assumptions of Theorem 2.1, if there holds

1 kue ðsÞk2 þ jdiv ue ðsÞj2 þ e2d0 s

e

Z

s

e2d0 t ð

0

sðsÞjuet ðsÞj2 þ e2d0 s

Z

s

e2d0 t sðtÞ

0

 1

e

1

2Nkf1k1 < 1, then for any s P 0, we have

kdiv ue k2 þ jAue j2 Þ dt 6 c;

e2

 jdiv uet j2 þ kuet k2 dt 6 c;

ð69Þ

lim sup kue ðsÞk 6 1 kf1 k1 ;

ð70Þ

s!1

1

s2 ðsÞðkuet ðsÞk2 þ jdiv uet ðsÞj2 Þ þ e2d0 s e

ð68Þ

Z

s

e2d0 t s2 ðtÞðjAuet j2 þ

0

sðsÞjAue ðsÞj2 þ s2 ðsÞjA1=2 uett ðsÞj2 þ e2d0 s

Z

s

1

e2

kdivuet k2 Þ dt 6 c;

e2d0 t ðsðtÞjA1=2 uett j2 þ s2 ðtÞjuett j2 Þ dt 6 c:

ð71Þ

ð72Þ

0

Proof. (68) and (69) have been given in [21]. Next, we derive the remained bounds. Taking

v ¼ e2d t ue in (21) and using (12), we get 0

1 d 2d0 t ðe jue j2 Þ þ me2d0 t kue k2 þ Jðt; ue ; e2d0 t ue Þ 6 e2d0 t ðf ; ue Þ þ d0 e2d0 t jue j2 : 2 dt

ð73Þ

4098

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

Integrating (73) for the time from 0 to s and multiplying by e2d0 s , we derive

Z s Z s jue ðsÞj2 þ 2me2d0 s e2d0 t kue k2 dt þ 2e2d0 s Jðt; ue ; e2d0 t ue Þ dt 0 0 Z s Z s e2d0 t jue j2 dt þ 2e2d0 s e2d0 t ðf ; ue Þ dt þ e2d0 s ju0 j2 : 6 2d0 e2d0 s 0

ð74Þ

0

Letting s ? 1 in (74), using the L’Hospital rule and considering

lim sup 2e2d0 s

Z

s!1

s

Jðt; ue ; e2d0 t ue Þ dt ¼

0

q

lim sup kue ðsÞk2 ;

dd0

s!1

it follows that:



mþ

q d

lim sup kue ðsÞk2 ¼ lim supðf ðsÞ; ue ðsÞÞ 6 kf1 k1 lim sup kue ðsÞk;

s!1

s!1

s!1

which implies (70). Rewrite (21) as

ðuet ; v Þ þ mðAe ue ; v Þ þ Jðt; ue ; v Þ þ bðue ; ue ; v Þ ¼ ðf ; v Þ:

ð75Þ

Differentiating (75) with respect to the time t, we have

ðuett ; v Þ þ mðAe uet ; v Þ þ qðAue ; v Þ þ bðuet ; ue ; v Þ þ bðue ; uet ; v Þ ¼ dJðt; ue ; v Þ þ ðft ; v Þ: Taking

ð76Þ

v = Aeuet in (76), we obtain

1 d 1 ðkuet k2 þ jdivuet j2 Þ þ mjAe uet j2 þ qðAue ; Ae uet Þ ¼ dJðt; ue ; Ae uet Þ  bðuet ; ue ; Ae uet Þ  bðue ; uet ; Ae uet Þ þ ðft ; Ae uet Þ: 2 dt e ð77Þ Using (16) and Lemma 2.6, there hold

jbðue ; uet ; Ae uet Þj 6 kue kkuet k1=2 jAuet j1=2 jAe uet j 6 kue kkuet k1=2 jAe uet j3=2 6

m

m 8

jAe uet j2 þ ckue k4 kuet k2 ;

jAe uet j2 þ ckue k4 kuet k2 ;

Z t

2



m jdJðt; ue ; Ae uet Þj 6 jAe uet j2 þ ce2dt

eds Aue ds

; 8 0

jbðuet ; ue ; Ae uet Þj 6

jðft ; Ae uet Þj 6

m

16

jqðAue ; Ae uet Þj 6

8

jAe uet j2 þ cjft j2 ;

m 16

jAe uet j2 þ cjAue j2 :

Taking these estimates into (77) and multiplying it by e2d0 t s2 ðtÞ, we get

Z t

2





d 2d0 t 2 1 1 e s ðtÞðkuet k2 þ jdiv uet j2 Þ þ me2d0 t s2 ðtÞjAe uet j2 6 ce2d0 t sðtÞðkuet k2 þ jdiv uet j2 Þ þ ce2a0 t

eds Aue ds

dt e e 0 þ ce2d0 t sðtÞjft j2 þ ce2d0 t jAue j2 þ ce2d0 t sðtÞkue k4 kuet k2 :

ð78Þ

Integrating (78) for the time from 0 to s, using Lemma 2.2, (68) and (69) and noting the fact that (see [19,21])

cðjAuet j2 þ

1

e2

kdivuet k2 Þ 6 jAe uet j2 ;

which is valid in a domain X 2 R2 with boundary oX 2 C2, we derive that





1

Z

s2 ðsÞ kuet ðsÞk2 þ jdivuet ðsÞj2 þ e2d0 s e

s

0

  1 e2d0 t s2 ðtÞ jAuet j2 þ 2 kdivuet k2 dt 6 c

ð79Þ

e

by a final multiplication of e2d0 s , which yields (71). Since (75) implies

mAe ue ¼ uet  qedt

Z

t

eds Aue ds  Bðue ; ue Þ þ f ;

0

applying (15) and the H ölder inequality, it follows

sðtÞm2 jAe ue ðtÞj2 6 sðtÞjuet ðtÞj2 þ ce2a0 t

Z 0

t

e2a0 s ds  e2d0 t

Z 0

t

e2d0 s jAue j2 ds þ csðtÞkue ðtÞk2 þ

m2 2

sðtÞjAue ðtÞj2 þ sðtÞjf ðtÞj2 :

4099

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

Combining the above inequality with (68)–(70), we arrive at

sðtÞjAue ðtÞj2 6 c: Taking

v ¼e

ð80Þ

s ðtÞuett in (76), we get

2d0 t 2

1 d 2d0 t 2 2 2d0 t 2 ðe s ðtÞmjA1=2 s ðtÞuett Þ þ bðue ; uet ; e2d0 t s2 ðtÞuett Þ e uet j Þ þ bðuet ; ue ; e 2 dt 2 2d0 t 2 s ðtÞðft ; uett Þ  qe2d0 t ðAue ; s2 ðtÞuett Þ þ de2d0 t Jðt; u; s2 ðtÞuett Þ: 6 ce2d0 t sðtÞjA1=2 e uet j þ e

e2d0 t s2 ðtÞjuett j2 þ

ð81Þ

It is valid that

jbðue ; uet ; s2 ðtÞuett Þj 6 cs2 ðtÞjAue jkuet kjuett j 1 6 s2 ðtÞjuett j2 þ cs2 ðtÞjAue j2 kuet k2 ; 8 1 jbðuet ; ue ; s2 ðtÞuett Þj 6 s2 ðtÞjuett j2 þ cs2 ðtÞjAue j2 kuet k2 ; 8 1 2 2 js ðtÞðft ; uett Þj 6 s ðtÞjuett j2 þ cs2 ðtÞjft j2 ; 16 1 2 jqðAue ; s2 ðtÞuett Þj 6 s ðtÞjuett j2 þ cs2 ðtÞjAue j2 ; 16

Z t

2



1 jJðt; u; s2 ðtÞuett Þj 6 s2 ðtÞjuett j2 þ ce2dt

eds Aue ds

: 8 0 Combining these estimates with (81) and using (71), we obtain

e2d0 t s2 ðtÞjuett j2 þ

Z t

2



1 d 2d0 t 2 2 1=2 2 2 2d0 t 2d0 t 2 2d0 t 2a0 t

ds

u j Þ 6 ce s ðtÞjA u j þ ce s ðtÞjf j þ ce jAu j þ e e Au d s ðe s ðtÞmjA1=2 et et t e e : e e

2 dt 0 ð82Þ 2d0 s

Integrating (82) from 0 to s, multiplying by e

e2d0 s

Z

s

and using (12), Theorem 2.1 and (68) and (69), we derive that

e2d0 t s2 ðtÞjuett j2 dt 6 c:

ð83Þ

0

Moreover, thanks to (76), there holds

uett ¼ mAe uet  qAue  Bðuet ; ue Þ  Bðue ; uet Þ þ qdedt

Z

t

eds Aue ds þ ft :

0

Applying (16), (68) and (69) and Lemma 2.2 to the above inequality, we have

s2 ðsÞjA1=2 uett ðsÞj2 þ e2d0 s

Z

s

e2d0 t sðtÞjA1=2 uett j2 dt 6 c;

0

which and (80), (83) yield (72). We complete the proof. h Next, we consider the following intermediate linear equations:

gt  mDg  q

Z

t

edðtsÞ Dgds þ rw ¼ f  Bðu; uÞ;

ð84Þ

0

e m

div g þ w ¼ 0 gð0Þ ¼ u0 ;

ð85Þ

where u is the solution of the viscoelastic fluid motion problem (1) and (2). Setting n = g  u, / = w  p, and subtracting (1) from (84), we get

nt  mDn  q

e m

Z 0

t

edðtsÞ Dnds þ r/ ¼ 0;

ð86Þ

e m

ð87Þ

div n þ / ¼  p;

nð0Þ ¼ 0:

Lemma 4.2. Under the assumptions of Theorem 2.1, it is valid for any s P 0 that

e2d0 s

Z 0

s

e2d0 t jnj2 dt þ sðsÞjnðsÞj2 þ s2 ðsÞknðsÞk2 þ e2d0 s

Z

s

e2d0 t s2 ðtÞj/j2 dt 6 ce2 :

ð88Þ

0

Proof. From Section 3, we note that the assumption (A2) for a linear problem can be replaced by the weaker condition f(x, t), s(t)ft(x, t) 2 L2(R+, Y) (or see Lemma 4.1 in [23]). Considering Theorem 2.1 and (16), we have

4100

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

e2d0 s

Z

s

e2d0 t jf  Bðu; uÞj2 dt 6 c þ e2d0 s

0

Z

s

e2d0 t jAuj2 kuk2 dt 6 c

0

and

e2d0 s

Z

s

e2d0 t s2 ðtÞj

0

@ ðf  Bðu; uÞÞj2 dt 6 c þ 2e2d0 s @t

Z

s

e2d0 t s2 ðtÞkut k2 jAuj2 dt 6 c:

0

Thus, Lemma 4.2 follows from Lemma 3.2 and Theorem 3.1. We complete the proof.

h

Letting w = ue  g, r = pe  w. Subtracting (84) from (3), we have that

wt  mDw  q

Z

t

edðtsÞ Dwds þ Bðue ; n þ wÞ þ Bðn þ w; uÞ þ rr ¼ 0;

ð89Þ

0

e m

div w þ r ¼ 0;

wð0Þ ¼ 0:

ð90Þ

The variational formulation of (89) reads as: Find (w, r) 2 X  M for all t P 0, such that

e m

ðwt ; v Þ þ aðw; v Þ þ Jðt; w; v Þ þ bðue ; n þ w; v Þ þ bðn þ w; u; v Þ  dðv ; rÞ þ dðw; qÞ þ ðr; qÞ ¼ 0:

ð91Þ

Theorem 4.1. Suppose (A1) and (A2) are valid, e ? 0 such that ec0 6 1, where c0 is a positive constant. If there holds 2Nkf1k1 < 1, then, for any s P 0, we have

sðsÞjuðsÞ  ue ðsÞj2 þ s2 ðsÞkuðsÞ  ue ðsÞk2 þ e2d0 s

Z 0

s

e2d0 t s2 ðtÞjp  pe j2 dt 6 ce2 :

ð92Þ

Remark 4.1. Under the assumptions of Theorem 4.1, by applying the similar precess as that used in [22,28], the solution ue can be expanded with respect to the penalty parameter e. Proof. Taking ðv ; qÞ ¼ e2d0 t ðw; rÞ in (91), we have

1 d 2d0 t 2 e e jwj þ me2d0 t kwk2 þ e2d0 t bðue ; n þ w; wÞ þ e2d0 t bðn þ w; u; wÞ þ Jðt; w; e2d0 t wÞ þ jrj2 ¼ d0 e2d0 t jwj2 : 2 dt m

ð93Þ

Thanks to (17), it is valid that

jbðn; u; wÞj þ jbðue ; n; wÞj 6 Nðkuk þ kue kÞknkkwk; jbðw; ue ; wÞj 6 Nkue kkwk2 : Combining these estimates with (93) and using Theorem 2.1, we obtain

d 2d0 t 2 e jwj þ 2ðm  Nkue kÞe2d0 t kwk2 þ 2Jðt; w; e2d0 t wÞ 6 2d0 e2d0 t jwj2 þ ce2d0 t ð1 þ kue kÞknkkwk: dt



ð94Þ

Integrating the above inequality from 0 to s, multiplying by e2d0 s and using Lemma 4.2, we arrive at

Z

jwðsÞj2 þ e2d0 s 6 2d0 e2d0 s

Z

s

2e2d0 t ðm  Nkue kÞkwk2 dt þ 2e2d0 s

0 s

Z

s

Jðt; w; e2d0 t wÞ dt

0

 1=2 Z s e2d0 t jwj2 dt þ ce e2d0 s e2d0 t ð1 þ kue k2 Þkwk2 dt :

0

ð95Þ

0

Letting s ? 1 in (95), using the L’Hospital rule and Theorem 2.1, Lemma 4.1, we get

ð  N 1 kf1 k1 Þ lim sup kwðsÞk2 6 ce lim sup kwðsÞk: s!1

Because of



2

s!1

Nkf1k1 < 1, there holds

lim sup jwðsÞj2 6 c lim sup kwðsÞk2 6 ce2 :

s!1

ð96Þ

s!1

On the other hand, multiplying (93) by s(t), applying (16) and Lemma 2.6, we obtain

m sðtÞjbðue ; n þ w; wÞj ¼ sðtÞjbðue ; w; n þ wÞj 6 sðtÞjAue jkwkje þ wj 6 sðtÞkwk2 þ csðtÞjAue j2 ðjnj2 þ jwj2 Þ; 4

m

sðtÞjbðn þ w; u; wÞj 6 sðtÞkwk2 þ csðtÞjAuj2 ðjnj2 þ jwj2 Þ: 4

4101

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

Integrating (93) for time from 0 to s, using Theorem 2.1, Lemmas 2.1, 2.4, 4.1, 4.2 and the above inequalities, we have

sðsÞe2d0 s jwðsÞj2 þ

Z

s

e m

e2d0 t sðtÞðkwk2 þ jrj2 Þ dt 6 ce2d0 s e2 þ

0

Z

s

e2d0 t jwj2 dt:

ð97Þ

0

When s 2 (0, T] with T being a finite time, applying Lemma 2.4 to (98) and multiplying by e2d0 s , it is valid that

sðsÞjwðsÞj2 þ e2d0 s

Z

s

e m

e2d0 t sðtÞðkwk2 þ jrj2 Þdt 6 ce2 :

0

ð98Þ

When s ? 1, considering (96) and (97), inequality (98) holds too. Differentiating dðw; qÞ þ me ðr; qÞ in (91) and taking ðv ; qÞ ¼ e2d0 t s2 ðtÞðwt ; rÞ, we obtain

e2d0 t s2 ðtÞjwt j2 þ

1 d 2d0 t 2 e e e s ðtÞðmkwk2 þ jrj2 Þ 6 ce2d0 t sðtÞðkwk2 þ jrj2 Þ  e2d0 t s2 ðtÞbðue ; w þ n; wt Þ 2 dt m m  e2d0 t s2 ðtÞbðn þ w; u; wt Þ  Jðt; w; e2d0 t s2 ðtÞwt Þ:

ð99Þ

Due to (16), it is valid

js2 ðtÞbðue ; w þ n; wt Þj 6 s2 ðtÞjAue jkw þ nkjwt j 1 6 s2 ðtÞjwt j2 þ cs2 ðtÞjAue j2 ðkwk2 þ knk2 Þ; 4 1 2 js ðtÞbðw þ n; u; wt Þj 6 s2 ðtÞjwt j2 þ cs2 ðtÞjAuj2 ðkwk2 þ knk2 Þ: 4 Integrating (99) for the time from 0 to s, using Theorem 2.1 and Lemmas 2.3, 4.1, 4.2 and (98), we obtain that

s2 ðsÞkwðsÞk2 þ e2d0 s

Z

s

e2d0 t s2 ðtÞjwt j2 dt 6 ce2 ;

ð100Þ

0

by a final multiplication by e2d0 s . From (89), there holds

rr ¼ wt þ mDw þ q

Z

t

edðtsÞ Dwds  Bðue ; n þ wÞ  Bðw þ n; uÞ:

ð101Þ

0

Considering (17), we have

kBðue ; n þ wÞk1 6 ckue kkw þ nk 6 ckue kðkwk þ knkÞ; kBðw þ n; uÞk1 6 ckukkw þ nk 6 ckukðkwk þ knkÞ: Therefore, using (98), (100) and the above equations, we derive

e2d0 s

Z

s

e2d0 t s2 ðtÞjrj2 dt 6 e2d0 s

Z

0

s

0

e2d0 t s2 ðtÞkrrk21 dt 6 ce2 :

Combining (98), (100), (102) with Lemma 4.2 and using the triangle inequality, we complete the proof.

ð102Þ h

5. Time discretization of the penalty method In this section, we will analyze the backward Euler time discretization scheme for the nonlinear penalty system. Let 0 < Dt < 1 denote the time step size and tn = nDt. For the smooth function / defined on [0, 1), set /n = /(tn). For the integral term, we apply the right rectangle rule as (see [16,31])

Mn ð/Þ ¼ Dt

n X

bnj /j 

Z

tn

bðt n  tÞ/ðtÞdt;

0

j¼1

where bnj = b(tn  tj). Let dn be the quadrature error associated with the quadrature rule. For / 2 C1[0, tn], it is defined by

dn ð/Þ :¼

Z

tn

bðt n  tÞ/ðtÞ dt  M n ð/Þ:

ð103Þ

0

It is valid that

sðtnþ1 Þ 6 sðtn Þ þ Dt; Dt 6 sðtn Þ; t  tn 6 sðtÞ 8t 2 ½tn ; tnþ1  and

/ðt nþ1 Þ 

1 Dt

Z

t nþ1

tn

/ðtÞdt ¼

1 Dt

Z

tnþ1

tn

ðt  t n Þ/t ðtÞdt;

ð104Þ

4102

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

In the case of 0 < tj < 1, s(tj) = tj = jDt. Therefore, using Lemma 4.1, Taylor formula and (104), we get

2



X nþ1 Z tj   2 @  dt 1=2 

nþ1 1=2

2

dt nþ1 e A ue dt

A ue ¼ q

e ðt  t j1 Þ

d

j¼1 tj1 @t



X nþ1 Z t j   2



1=2 1=2 dðt nþ1 tÞ 6 c

e ðt  t j1 Þ dA ue þ A uet dt



j¼1 tj1 

Z Z

2

2  nþ1 nþ1 tj tj X X e2a0 ðtnþ1 tÞ





6 c Dt 2 e2d0 t ðsðtj1 Þ þ DtÞ A1=2 ue þ A1=2 uet dt dt  e2d0 tnþ1 sðtj Þ t j1 tj1 j¼1 j¼1 

nþ1 2a0 t j nþ1 Z t j

2

2  X X e





6 cDt 2 e2a0 tnþ1 e2d0 t sðtÞ A1=2 ue þ A1=2 uet dt e2d0 tnþ1 j t j1 j¼1 j¼1 !  nþ1  nþ1 X X 1 2 2a0 t nþ1 2a0 tj 2 2a0 t nþ1 2a0 t nþ1 2a0 tj 6 c Dt 2 ; 6 c Dt e 6 c Dt e þ 2a0 Dte ðlnðn þ 1Þ þ rÞ þ cDte e j j¼1 j¼1 ð105Þ where r is the Euler constant. Similarly, for jdn + 1(Aue)j, we have nþ1  a0 t j 2 nþ1 Z t j   X X

nþ1

e

d ðAue Þ 2 6 cDt2 e2a0 tnþ1 e2d0 tnþ1 e2d0 t s2 ðtÞ jAue j2 þ jAuet j2 dt j tj1 j¼1 j¼1 ! nþ1 X 1 6 cDt2 e2a0 tnþ1 þ a20 Dt2 e2a0 tj 6 cDt 2 : 2 j j¼1

ð106Þ

While in the case of tj P 1, s(tj) = 1, (105) and (106) hold, too. The backward Euler time discretization scheme of (21) is as follows:

 nþ1      nþ1   m  ue  une nþ1 nþ1 ðAue Þ; v þ b unþ1 rdiv unþ1 ; v Þ; ; v þ a unþ1 e ;v þ M e ; ue ; v  e ; v ¼ ðf Dt e

ð107Þ

with u0e ¼ u0 and fn+1 = f(tn+1). Lemma 5.1. Under the assumptions of Theorem 2.1, if

2Nm1 1 kf1 k1 < 1;

ð108Þ

then, for any n > 0, it is valid that

sðtn Þjen j2 þ Dte2d0 tn s2 ðtn Þken k2 þ Dte

n X



2

m

2 e2d0 tm sðt m Þ A1=2 e e 6 cDt ;

m¼1 n X 2d0 t n 2d0 t m

e

ð109Þ

s2 ðtm ÞjAe em j2 6 cDt2 ;

ð110Þ

m¼1

where en ¼ ue ðtn Þ  une with e0 = 0. Proof. Subtracting (107) from (21) at t = tn+1, we obtain

 nþ1            e  en nþ1 ; v þ b enþ1 ; ue ; v ¼ Rnþ1 ; v þ m Ae enþ1 ; v þ dnþ1 ðAue Þ; v þ b unþ1 e ;e e ;v ; Dt

ð111Þ

where

Rnþ1 ¼ uet ðtnþ1 Þ  e Taking

1 1 ðue ðtnþ1 Þ  ue ðt n ÞÞ ¼ Dt Dt

Z

tnþ1

ðt  tn Þuett dt:

ð112Þ

tn

v = 2Dts(tn+1)en+1 in (111), we derive that 





2

nþ1

sðtnþ1 Þ jenþ1 j2  jen j2 þ jenþ1  en j2 þ 2Dtsðtnþ1 Þm

A1=2

þ 2Dt sðt nþ1 Þðdnþ1 ðAue Þ; enþ1 Þ e e

    nþ1 þ 2Dt sðt nþ1 Þb enþ1 ; ue ; enþ1 ¼ 2Dt sðt nþ1 Þ Rnþ1 e ;e

ð113Þ

4103

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

or



2



nþ1

sðtnþ1 Þjenþ1 j2  sðtn Þjen j2  Dtjen j2 þ 2Dtsðtnþ1 Þm

A1=2

þ 2Dtsðtnþ1 Þ dnþ1 ðAue Þ; enþ1 e e



  þ 2Dtsðtnþ1 Þb enþ1 ; ue ; enþ1   nþ1 : 6 2Dt sðt nþ1 Þ Rnþ1 e ;e

ð114Þ

Using (18), (105), Lemma 2.6 and Lemma 4.1, it follows that

 

 



2Dt sðt nþ1 Þ dnþ1 ðAue Þ; enþ1 ¼ 2Dtsðtnþ1 Þ dnþ1 ðA1=2 ue Þ; A1=2 enþ1



2   2 m Dt



1=2 nþ1

6 sðtnþ1 Þ

A1=2

þ cDt sðt nþ1 Þ dnþ1 A ue

e e 4

2 m Dt nþ1

6 sðtnþ1 Þ

A1=2

þ c Dt 3 ; e e 4





  



1=2 nþ1 1=2 nþ1

1=2 nþ1 1=2 nþ1 nþ1

Re ; A e Re jA e j 2Dt sðt nþ1 Þ Rnþ1

¼ 2Dt sðt nþ1 Þ A

6 2Dt sðt nþ1 ÞjA e ;e

Z tnþ1

2

mDt csðtnþ1 Þ



nþ1 2 6 sðtnþ1 Þ

A1=2 j þ ðt  tn Þ A1=2 uett

e e 2 Dt tn Z tnþ1 Z

2

2 mDt sðtnþ1 Þ tnþ1 nþ1

6 sðtnþ1 Þ

A1=2 e þ c ðt  t Þdt  sðtn Þ

A1=2 uett

dt

n e 2 sðtn Þ tn tn Z tnþ1





mDt

1=2 nþ1 2

1=2 2 2 sðtnþ1 Þ Ae e þ cDt sðtÞ A uett dt 6 2 tn and

 



2Dt sðt nþ1 Þ bðenþ1 ; ue ; enþ1 Þ 6 cDtsðtnþ1 Þ jenþ1 jkenþ1 kkue k þ jenþ1 j1=2 kenþ1 k3=2 kue k

2 mDt 2 nþ1

6 sðtnþ1 Þ

A1=2

þ cDtjenþ1 j : e e 4

ð115Þ

ð116Þ

For m 6 M with M being a finite number, taking the above estimate into (114), summing it from 0 to m and using Lemma 4.1, we have

sðtmþ1 Þjemþ1 j2 þ Dt

m X n¼0

2



nþ1

sðtnþ1 Þ

A1=2

6 c Dt 2 þ c Dt e e

m X

jenþ1 j2 :

n¼0

Applying Lemma 2.5 to the above inequality, we have

sðtmþ1 Þjemþ1 j2 6 cDt2 8m 2 ð0; M:

ð117Þ

Furthermore, when n ? 1, there holds s(tn+1) = 1. We estimate the trilinear term in (113) as

 nþ1 

b e ; ue ; enþ1 6 Nkue kkenþ1 k2 6 N1 kf1 k kenþ1 k2 ; 1

ð118Þ

by using (70). Taking (118) into (113) and considering (72), it follows that

jenþ1 j2  jen j2 þ jenþ1  en j2 þ



Thanks to (108), it follows c: = m  2N

jenþ1 j2 6



m  2N1 kf1 k1 Dtkenþ1 k2 6 cDt3 :

ð119Þ

1

kf1k1 > 0. Using the Poincare inequality, we deduce

n X 1 c Dt 3 1 1 je0 j2 þ cDt 3 jen j2 þ 6 nþ1 mþ1 1 þ k1 cDt 1 þ k1 cDt ð1 þ k1 cDtÞ ð1 þ k 1 cDtÞ m¼0

6 c Dt 3

1  ð1 þ k1 cDtÞðnþ1Þ 6 c Dt 2 k1 cDt

8n 2 ðM; 1Þ;

ð120Þ

which and (117) imply

sðtmþ1 Þjemþ1 j2 6 cDt2 8m > 0: 2d0 t nþ1

Multiplying (114) by e

ð121Þ 2d0 Dt

and noting the fact that c P e

P 1 þ 2d0 Dt, we get

2

nþ1

e2d0 tnþ1 sðt nþ1 Þjenþ1 j2  e2d0 tn sðtn Þjen j2 þ 2Dte2d0 tnþ1 sðt nþ1 Þm A1=2

þ 2Dte2d0 tnþ1 sðtnþ1 Þðdnþ1 ðAue Þ; enþ1 Þ e e     nþ1 : þ 2Dte2d0 tnþ1 sðt nþ1 Þb enþ1 ; ue ; enþ1 6 cDte2d0 tn jen j2 þ 2Dte2d0 tnþ1 sðt nþ1 Þ Rnþ1 e ;e

ð122Þ

4104

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

Taking the summation of (122) from 0 to m, using (121), Lemma 4.1 and multiplying by e2d0 tmþ1 , we derive that

Dte2d0 tmþ1

m X n¼0



2 m m X X

nþ1

e2d0 tnþ1 sðtnþ1 Þ A1=2 e2d0 tnþ1 jenþ1 j2 þ cDt3 e2d0 tmþ1 e2d0 tnþ1

6 cDte2d0 tmþ1 e e n¼0

þ cDt 2 e2d0 tmþ1

Z

n¼0

t mþ1

0



2



e2d0 t sðtÞ A1=2 uett dt 6 cDt 2 :

ð123Þ

Combining (121) with (123), we obtain (109). Thus, using Lemma 4.1, (109) and the triangle inequality, it is valid that



2

6 sðt nþ1 Þkenþ1 k2 þ kue ðtnþ1 Þk2 6 cDt þ c 6 2c: sðtnþ1 Þ unþ1 e Taking

ð124Þ

v ¼ 2Dte2d t s2 ðtnþ1 ÞAe enþ1 in (111), we obtain 0 nþ1





2 

2



2



1=2 2 1=2 nþ1

e2d0 tnþ1 s2 ðt nþ1 Þ A1=2

 Ae en þ Ae ðenþ1  en Þ þ 2Dte2d0 tnþ1 s2 ðt nþ1 Þm Ae enþ1

e e     nþ1 þ 2Dte2d0 tnþ1 s2 ðtnþ1 Þ dnþ1 ðAue Þ; Ae enþ1 þ 2Dte2d0 tnþ1 s2 ðt nþ1 Þb unþ1 ; Ae enþ1 e ;e     nþ1 þ 2Dte2d0 tnþ1 s2 ðtnþ1 Þb enþ1 ; ue ; Ae enþ1 ¼ 2Dte2d0 tnþ1 s2 ðt nþ1 Þ Rnþ1 : e ; Ae e

ð125Þ

Considering (16), (106), (115), (124), Lemma 2.6 and Lemma 4.1, we have

 





 

2Dts2 ðt nþ1 Þ dnþ1 ðAue Þ; Ae enþ1 6 mDt s2 ðt nþ1 Þ Ae enþ1 2 þ cDt3 ;

2Dte2d0 tnþ1 s2 ðt nþ1 Þ Rnþ1 ; Ae enþ1

e 4 Z tnþ1

2 mDt 2d0 tnþ1 2 6 e s ðtnþ1 Þ Ae enþ1 þ cDt2 e2d0 t s2 ðtÞjuett j2 dt; 2 tn  nþ1 1=2 nþ1 1=2 nþ1

       nþ1 ke k jAe j Ae e

Dts2 ðt nþ1 Þ b unþ1 ; Ae enþ1 þ b enþ1 ; ue ; Ae enþ1 6 cDts2 ðt nþ1 Þ kue k þ unþ1 e ;e e

2

2 mDt 2 nþ1

s ðtnþ1 Þ Ae enþ1 þ cDtsðtnþ1 Þ

A1=2 6

: e e 4 Taking the summation of (125) for n from 0 to m, using the above inequalities and the relation that





2







2



1=2 2

1=2 2

1=2

1=2 2 nþ1

s2 ðtnþ1 Þ

A1=2

 s2 ðtn Þ Ae en  2Dt Ae en 6 s2 ðt nþ1 Þ Ae enþ1  Ae en

e e



and making the similar analysis as in (123), we get, after a final multiplication by e2d0 tmþ1 , that



2

mþ1

s2 ðtmþ1 Þ

A1=2

þ Dte2d0 tmþ1 e e

m X n¼0

m

2 X

2

nþ1

e2d0 tnþ1 s2 ðt nþ1 Þ Ae enþ1 6 cDt 2 þ cDte2d0 tmþ1 e2d0 tnþ1 sðtnþ1 Þ A1=2

: e e

ð126Þ

n¼0

Thanks to (12) and (109), we complete the proof. h Theorem 5.1. Suppose that assumptions (A1)-(A2) are valid, e ? 0 such that ec0 6 1, where c0 is a positive constant. (u(tn),p(tn)) and ðune ; pne Þ are the solutions of the problem (1) and (2) and the system (107), respectively. If there holds 2Nm11kf1k1 < 1, then, for any tn P 0, we have









sðtn Þ uðtn Þ  une 2 þ s2 ðtn Þ uðtn Þ  une 2 6 cðDt2 þ e2 Þ; Dte2d0 tn

n X



2

6 cðDt2 þ e2 Þ: e2d0 tm s2 ðt m Þ pðt m Þ  pm e

ð127Þ ð128Þ

m¼1

Proof. (127) follows from Theorem 4.1 and Lemma 5.1. To estimate (128), taking

v ¼ e2d t s2 ðtnþ1 Þðenþ1  en Þ in (111), we have 0 nþ1

2

nþ1 



2 





e  en

1=2 nþ1 2 1=2 n 2 1=2 nþ1 2d0 t nþ1 2 n

Dte2d0 tnþ1 s2 ðt nþ1 Þ

þ e s ðt Þ m A e  A e þ A ðe  e Þ







nþ1 e e e Dt

    nþ1 ¼ e2d0 tnþ1 s2 ðtnþ1 Þ Rnþ1  en Þ  e2d0 tnþ1 s2 ðtnþ1 Þ dnþ1 ðAue Þ; ðenþ1  en Þ e ; ðe     nþ1 nþ1  e2d0 tnþ1 s2 ðtnþ1 Þb unþ1 ;e  en  e2d0 tnþ1 s2 ðt nþ1 Þb enþ1 ; ue ; enþ1  en : e ;e

ð129Þ

By applying a similar approach as that used in Lemma 5.1, and using (109) and (110), we see that

Dte2d0 tmþ1

nþ1

2

e  en

e2d0 tnþ1 s2 ðt nþ1 Þ

6 c Dt 2 : Dt

n¼0

m X

ð130Þ

K. Wang et al. / Applied Mathematical Modelling 34 (2010) 4089–4105

4105

Then using (111), Lemma 5.1 and Theorem 4.1, we can prove

Dte2d0 tmþ1

m X



2

6 cðDt2 þ e2 Þ: e2d0 tnþ1 s2 ðt nþ1 Þ pðtnþ1 Þ  pnþ1 e



ð131Þ

n¼0

6. Conclusions In this paper, we deduced the optimal error estimates for the penalty system and its time semi-discrete scheme of the two-dimensional viscoelastic fluid motion problem under some assumptions of the data (u0, f) and the parameters e and Dt. It suggests proper choices of e to solve numerically. All of these results are uniform in time. Obviously, we can easily extend the present analysis to a fully discrete scheme by combining it with the finite element approximation results in [14]. Acknowledgments The authors thank the editor and reviewers for their criticism, valuable comments, and suggestions which helped to improve the results of this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

[22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

R. Bird, R. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Fluid Mechanics, vol. 1, John Wiley & Sons, New York, 1987. D. Joseph, Fluid Dynamics of Viscoelastic Liquids, Springer-Verlag, New York, 1990. J. Oldroyd, On the formulation of the rheological equations of state, Proc. Roy. Soc. Lond. Math. Phys. Sci. 200 (1950) 523–541. A. Oskolkov, Initial boundary value problems for the equations of motion of Kelvin–Voigt fluids and Oldroyd fluids, Proc. Steklov Inst. Math. 2 (1989) 137–182. V. Ortov, P. Sobolevskii, On mathematical models of a viscoelastic with a memory, Differ. Inte. Equat. 4 (1991) 103–115. P. Sobolevskii, Stabilization of viscoelastic fluid motion (Oldroyd’s mathematical model), Differ. Inte. Equat. 7 (1994) 1579–1621. G. Araújo, S. Menezes, A. Marinho, Existence of solutions for an Oldroyd model of viscoelastic fluids, J. Differ. Equat. 2009 (2009) 1–16. Y. Agranovich, P. Sobolevskii, Investigation of viscoelastic fluid mathematical model, RAC. Ukranian SSR. Ser. A 10 (1989) 71–74. Y. Agranovich, P. Sobolevskii, Motion of non-linear viscoelastic fluid, Nonlinear Anal. 32 (1998) 755–760. J. Cannon, R. Ewing, Y. He, Y. Lin, A modified nonlinear Galerkin method for the viscoelastic fluid motion equations, Int. J. Eng. Sci. 37 (1999) 1643– 1662. Y.N. He, Y.P. Lin, S.P. Shen, R. Tait, Convergence to steady state of the viscoelastic fluid motion, Adv. Differ. Equat. 7 (2002) 717–742. P. Sobolevskii, Asymptotic of stable viscoelastic fluid motion (Oldroyd’s mathematical model), Math. Nachr. 177 (1996) 281–305. M. Akhmatov, A. Oskolkov, On convergence difference schemes for the equations of motion of an Oldroyd fluid, J. Soviet. Math. 47 (1989) 2926–2933. Y.N. He, Y.P. Lin, S.P. Shen, W.W. Sun, R. Tait, Finite element approximation for the viscoelastic fluid motion problem, J. Comput. Appl. Math. 155 (2003) 201–222. A. Pani, J. Yuan, Semidiscrete finite element Galerkin approximation to the equations of motion arising in the Oldroyd model, IMA J. Numer. Anal. 25 (2005) 750–782. A. Pani, J.Y. Yuan, P. Damazio, On a linearized backward Euler method for the equations of motion of Oldroyd fluids of order one, SIAM J. Numer. Anal. 44 (2006) 804–825. K. Wang, Y.N. He, Y.Q. Shang, Fully discrete finite element method for the viscoelastic fluid motion equations, Discrete Contin. Dyn. Syst. -Ser. B 13 (2010) 665–684. Y.N. He, Y. Li, Asymptotic behavior of linearized viscoelastic flow problem, Discrete Contin. Dyn. Syst. -Ser. B 10 (2008) 843–856. B. Brefort, J. Ghidaglia, R. Temam, Attractors for the penalized Navier–Stokes equations, SIAM J. Math. Anal. 19 (1988) 1–21. Y.N. He, Optimal error estimate of the penalty finite element method for the time-dependent Navier–Stokes equations, Math. Comp. 74 (2005) 1201– 1216. A. Katsiolis, A. Oskolkov, Initial boundary-value problems for equations of slightly compressible Jeffreys–Oldroyd fluids, Zap. Nauchn. Semin. POMI 208 (1993) 200–218; A. Katsiolis, A. Oskolkov, Initial boundary-value problems for equations of slightly compressible Jeffreys–Oldroyd fluids, J. Math. Sci. 81 (1996) 2578– 2588. A. Oskolkov, The penalty method for equations of viscoelastic media, Zap. Nauchn. Semin. POMI 224 (1995) 267–278; A. Oskolkov, The penalty method for equations of viscoelastic media, J. Math. Sci. 88 (1998) 283–291. J. Shen, On error estimates of the penalty method for unsteady Navier–Stokes equations, SIAM J. Numer. Anal. 32 (1995) 386–403. R. Temam, Une méthode d’approximation des solutions deséquations de Navier–Stokes, Bull. Soc. 98 (1968) 115–152. France. R. Courant, Variational methods for the solution of problems of equilibrium and vibrations, Bull. Am. Math. Soc. 49 (1943) 1–23. R. Adams, Sobolev Spaces, Academic Press, New York, 1975. J. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier–Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization, SIAM J. Numer. Anal. 19 (1982) 275–311. R. Temam, Navier–Stokes Equations, Theory and Numerical Analysis, North-Holland, Amsterdam, 1984. V. Giraut, P. Raviart, Finite Element Approximation of the Navier–Stokes Equations, Springer-Verlag, Berlin, 1979. J. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier–Stokes problem part IV: error analysis for second-order time discretization, SIAM J. Numer. Anal. 27 (1990) 353–384. W. Allegretto, Y.P. Lin, A.H. Zhou, Long-time stability of finite element approximations for parabolic equations with memory, Numer. Methods Part. Differ. Equat. 15 (1999) 333–354.