A priori error estimates of fully discrete finite element Galerkin method for Kelvin-Voigt viscoelastic fluid flow model

# A priori error estimates of fully discrete finite element Galerkin method for Kelvin-Voigt viscoelastic fluid flow model

Saumya Bajpai
School of Mathematics and Computer Science,
Indian Institute of Technology Goa, Ponda-403401, India
and
Ambit K. Pany
Department of Mathematics, Gandhi Institute for Technological Advancement,
Bhubaneswar-752054, India
###### Abstract

In this article, a finite element Galerkin method is applied to the Kelvin-Voigt viscoelastic fluid model, when its forcing function is in . Some new a priori bounds for the velocity as well as for the pressure are derived which are independent of inverse powers of the retardation time . Optimal error estimates for the velocity in as well as in -norms and for the pressure in -norm of the semidiscrete method are discussed which hold uniformly with respect to as with the initial condition only in . Further, under uniqueness condition, these estimates are shown to be uniformly in time as . For the complete discretization of the semidiscrete system, a first-order accurate backward Euler method is applied and fully discrete optimal error estimates are established. Finally, numerical experiments are conducted to verify the theoretical results. The results derived in this article are sharper than those derived earlier for finite element analysis of the Kelvin-Voigt fluid model in the sense that the error estimates in this article hold true uniformly even as .

Keywords: Kelvin-Voigt viscoelastic model, a priori estimates, semidiscrete finite element Galerkin method, fully discrete optimal error estimates, uniqueness condition.

## 1 Introduction

The equations of motion arising from the Kelvin-Voigt model give rise to the following system of partial differential equations :

 (1.1) ∂u∂t+u⋅∇u−κΔut−νΔu+∇p=f(x,t),x∈Ω,t>0,

and incompressibility condition

 (1.2) ∇⋅u=0,x∈Ω,t>0,

with initial and boundary conditions

 (1.3) u(x,0)=u0inΩ,u=0,on∂Ω,t≥0.

Here, is a bounded convex polygonal or polyhedral domain in with boundary , represents the velocity vector, is the pressure of the fluid, is the external force, denotes the kinematic coefficient of viscosity and is the retardation time. For a more physical description and applications of the model, one may refer [8]-[10], [18] and literature therein. Based on the the proof techniques of Ladyzenskaya [17], Oskolkov and his collaborators [18], [19], [21], [22] have discussed the existence of a unique global “almost “ classical solution for the initial and boundary value problem (1.1)-(1.3) for various assumptions on the right-hand side function and for all time .
There is a considerable amount of literature devoted to the numerical approximations of Kelvin-Voigt fluid flow model, see [2]-[5], [16], [20], [25]-[28]. In [20], Oskolkov has applied the spectral Galerkin approximation to the problem (1.1)-(1.3) and has proved the convergence for with the assumption that the solution is asymptotically stable as . Further, the author established optimal error estimates in -norm, which are local in time, since the constants appearing in error bounds involve exponential in time terms. Later on, as an improvement to the Oskolkov work, Pani et al. [25] have established and -norms optimal error estimates for the spectral Galerkin method applied to (1.1)-(1.3), which are valid uniformly in time under uniqueness assumption. They further applied modified nonlinear Galerkin method to (1.1)-(1.3), and have established optimal uniform in time a priori error estimates with the assumption of uniqueness condition. They have also observed the superconvergence phenomenon in -norm for both spectral Galerkin method and modified nonlinear spectral Galerkin method. Note that, the constants appearing in error estimates derived in [2]-[5], [16], [20], [25] depend on , for which may blow as .
In [26], the authors have applied semidiscrete finite element Galerkin method to the problem (1.1)-(1.3) and have established some new uniform in time a priori bounds for the weak solution. It can be observed that the constants appearing in a priori bounds for the weak solution are independent of inverse powers of which is an improvement over the results derived in earlier articles related to the regularity estimates for the weak solution of this model. Further, using these a priori estimates, they have established optimal error estimates for the velocity in as well as in - norms and for the pressure in -norm, when the forcing function . Here, it can be noted that they have achieved an improvement in the error estimates in powers of as the constants in error bounds depend only on .
As an extension to the work in [26], Pany et al. [27], [28] have employed a linearized first order backward Euler method and a second order backward difference scheme for the time discretization of the problem (1.1)-(1.3) with and have derived a priori bounds for the discrete solution in the Dirichlet norm using a combination of discrete Gronwall’s lemma and Stolz-Cesaro’s classical result for sequences. Then, making use of these a priori estimates for the solution, they have established fully discrete optimal error estimates for the velocity and pressure, which hold true uniformly in time under uniqueness assumption. In [28], the author has also mentioned that assuming the solution is smooth enough, that is, with on , the optimal error estimates independent of can be achieved following the similar analysis as in [26]-[28]. For the articles related to the finite element analysis of the problem (1.1)-(1.3) with the right-hand side forcing function , one may refer to [2]-[5]. For the papers containing the similar results for the Navier-Stokes and Oldroyd models, see [1], [11]-[14], [23], [24], [30], [31] and literature, referred therein.
Since the Kelvin-Voigt fluid is characterized by the fact that after instantaneous removal of the stresses, the velocity of the fluid does not vanish instantaneously but dies out like [19], it is worthwhile to discuss the behavior of the solution as and as . Moreover, this model can be thought of as a regularization of the Navier-Stokes model ([15], [17]). Based on these observations, in this article, we mainly aim at recovering optimal error estimates which are valid uniformly in time as well as in retardation time under realistically assumed minimum regularity assumption on the exact solution with and , .

The main contributions of the present article are as follows:
(i) Some new regularity results for the higher order time derivatives of the weak solution are derived which are valid uniformly in time. Further, these estimates are shown to be uniformly in as under minimum regularity assumptions and , . Here, it can be noted that the introduction of weight function plays a key role in handling the regularity issues at .
(ii) Using the Sobolev-Stokes projection defined earlier in [2], fully discrete optimal error estimates in and -norms for the finite element velocity approximation and in -norm for the finite element pressure approximation are established. It is further proved that these error estimates hold uniformly as . Here, we would like to highlight an important point that we have resorted to a simple observation (mentioned in Remarks 4.1, 4.2) in order to derive the estimates involving weight function which plays an important role in achieving uniform estimates in terms of .
(iii) Since the error bounds derived in (ii) involve exponential in time terms, it is further established that under the assumption of uniqueness condition, the error estimates are uniformly in time.
(iv) Numerical results are presented to validate our theoretical findings. Moreover, it is depicted that the order of convergence does not degenerate as confirming the results in

Note that, the results in this article are substantial improvements over the results available in literature related to the finite element error analysis of the Kelvin-Voigt model in the sense that we are able to establish error bounds which do not involve inverse powers of . As a consequence, the error estimates do not blow up as . The main difficulty in making error estimates independent of arises due to the lack of regularity of solution at . In order to overcome this difficulty, we introduce various powers of weight function which takes care of regularity issues of the solution at .
The remaining part of the article consists of the following sections. In Section 2, some preliminaries to be used in the subsequent sections are introduced and some new regularity results for the weak solution are derived. In Section 3, assumptions on finite element spaces to determine the discrete solution are presented and semidiscrete finite element approximations are defined. The main results of the article are also stated. Section 4 deals with the optimal error estimates for velocity and pressure. In Section 5, full discretization is achieved by using the backward Euler method. Section 6 presents some numerical results which confirm our theoretical findings. Finally, Section 7 concludes the article by briefly summarizing the results.

## 2 Preliminaries and Weak formulation

We denote -valued function spaces using bold face letters, that is, , and where is the space of square integrable functions defined in with inner product and norm . Further, denotes the standard Hilbert Sobolev space of order with norm . The space is equipped with a norm . Given a Banach space endowed with norm , let be the space of all strongly measurable functions satisfying if and for , . Also, we define the divergence free spaces

 J = {\boldmathϕ∈L2:∇⋅% \boldmathϕ=0\rm inΩ,\boldmathϕ⋅n|∂Ω=0\rm holdsweakly}, J1 = {\boldmathϕ∈H10:∇⋅% \boldmathϕ=0},

where is the unit outward normal to the boundary and should be understood in the sense of trace in , see [29]. Let be the quotient space with norm . For , it is denoted by . Now, define as the -orthogonal projection.
(A1). Setting as the Stokes operator, assume that the following regularity result holds:

 (2.1) ∥v∥2≤C∥~Δv∥∀v∈J1∩H2.

The above assumption is valid as the domain is a convex polygon or convex polyhedron. It can be noted that the following Poincaré inequality [13] holds true:

 (2.2) ∥v∥2≤λ−11∥∇v∥2∀v∈H10(Ω),

where , is the best possible positive constant depending on the domain Further, observe (see, [13]) that

 (2.3) ∥∇v∥2≤λ−11∥~Δv∥2∀v∈J1∩H2.

(A2). There exists a positive constant such that the initial velocity and the external force satisfy for with

 u0∈H2∩J1,f,ft∈L∞(0,T;L2)with∥u0∥2≤M,esssup0

The weak formulation of (1.1)-(1.3) is to find , such that and for

 (2.4) (ut,\boldmathϕ)+κ(∇ut,∇\boldmathϕ)+ν(∇u,∇\boldmathϕ)+(u⋅∇u,\boldmathϕ% )−(p,∇⋅\boldmathϕ)=(f,\boldmathϕ)∀\boldmathϕ∈H10,(∇⋅u,χ)=0∀χ∈L2.}

Equivalently, find such that for , ,

 (2.5) (ut,\boldmathϕ)+κ(∇ut,∇\boldmathϕ)+ν(∇u,∇\boldmathϕ)+(u⋅∇u,\boldmathϕ)=(f,\boldmathϕ)∀\boldmathϕ∈J1.

For , define the bilinear form as

 a(v,\boldmathϕ)=(∇v,∇\boldmathϕ)

and the trilinear form as

 b(v,w,\boldmathϕ)=12(v⋅∇w,\boldmathϕ)−12(v⋅∇\boldmathϕ,w).

Note that, for , , , . Because of antisymmetric property of the trilinear form, it is easy to verify that

 b(v,w,w)=0∀v,w∈J1.

We present below in Lemma 2.1, some a priori bounds for the weak solution pair which will be used in our subsequent error analysis. Since the estimates in (2.6) and (2.7) are already derived in [26], we only provide proof of (2.8).

###### Lemma 2.1.

[[26], pp 241, 244] Let the assumptions (A1)-(A2) hold. Then, there exists a positive constant such that for the following estimates hold true:

 (2.6) sup0

where and .

Proof. We know that , where . Hence, or .
Now, consider the following two cases:
Case 1: . Then,

 (2.9) 1σ∫t0e2αs∥~Δu(s)∥2ds ≤1e2αtsup0

A use of (2.6) in (2.9) leads to

 (2.10) 1σ∫t0e2αs∥~Δu(s)∥2ds

Case 2: . Again use (2.6) and well known facts of series to obtain

 1σ∫t0e2αs∥~Δu(s)∥2ds ≤1te2αtsup0

Therefore, considering the above two cases, we arrive at

 σ−1(t)∫t0e2αs∥u(s)∥22ds≤C.

Following the similar sets of arguments as above, we obtain

 σ−1(t)∫t0e2αs(∥p(s)∥2H1/R+κ∥~Δus(s)∥2)ds≤C

and this completes the remaining part of the proof.
In the next lemma, we derive a priori bounds for the highest order time derivatives of weak solution for the problem (2.4).

###### Lemma 2.2.

Let the assumptions (A1)-(A2) hold. Then, there exists a positive constant such that for the following estimates hold true:

 (2.12) τ(t)(∥∇ut∥2+κ∥~Δut∥2)+νe−2αt∫t0σ(s)(∥~Δus(s)∥2+∥∇ps(s)∥2)ds≤C, (2.13) e−2αt∫t0σ(s)(∥uss(s)∥2+κ∥∇uss(s)∥2)ds≤C, (2.14) e−2αt∫t0σ(s)(κ∥∇uss(s)∥2+κ2∥~Δuss(s)∥2)ds≤C.

Note that, here and everywhere else in the consecutive analysis the constant is independent of inverse powers of .

Proof. Rewrite (2.5) and differentiate the resulting equation with respect to time to arrive at

 (2.15) (utt,\boldmathϕ)−κ(~Δutt,\boldmathϕ)−ν(~Δut,\boldmath% ϕ)+(ut⋅∇u+u⋅∇ut,% \boldmathϕ)=(ft,\boldmathϕ),∀\boldmath% ϕ∈J1.

Choose in (2.15) to obtain

 12ddtσ(∥∇ut∥2+κ∥~Δut∥2) +νσ∥~Δut∥2=C(α,λ1)e2αt(∥∇ut∥2+κ∥~Δut∥2) +σ(ut⋅∇u,~Δut)+σ(u⋅∇ut,~Δut)−σ(ft,~Δut) (2.16) =C(α,λ1)e2αt(∥∇ut∥2+κ∥~Δut∥2)+I1+I2+I3,(% say).

A use of Cauchy-Schwarz’s inequality and Young’s inequality lead to

 |I1|+|I2| ≤C(ϵ)σ(∥ut∥∥∇ut∥∥∇u∥∥~Δu∥+∥∇u∥∥~Δu∥∥∇ut∥2)+ϵσ∥~Δut∥2 (2.17) ≤C(ϵ,λ1)σ∥∇ut∥2∥~Δu∥2+ϵσ∥~Δut∥2.

Once again, apply Cauchy-Schwarz’s inequality and Young’s inequality to bound as

 (2.18) |I3|≤Cσ∥ft∥∥~Δut∥≤C(ϵ)σ∥ft∥2+ϵσ∥~Δut∥2.

After using (2)-(2.18) in (2) with a proper choice of , integrate the resulting equation with respect to time from to to arrive at

 σ(∥∇ut∥2 +κ∥~Δut∥2)+ν∫t0σ(s)∥~Δus(s)∥2ds≤C(α,λ1,ν)(∫t0e2αs(∥∇us(s)∥2+κ∥~Δus(s)∥2)ds (2.19) +∫t0σ(s)∥∇us(s)∥2∥~Δu(s)∥2ds+∫t0σ(s)∥fs(s)∥2ds).

Apply Lemma 2.1 and assumption (A2) in (2). Then, multiply the resulting equation by to arrive at the desired a priori estimates of in (2.12).
Next, differentiate (2.5) and substitute in the resulting equation to observe that

 (2.20) σ(∥utt∥2+κ∥∇utt∥2)=−νσ(∇ut,∇utt)−σ(ut⋅∇u+u⋅∇ut,utt)+σ(ft,utt).

After rewriting the first term on the right-hand side of (2.20), apply Cauchy-Schwarz’s inequality, Young’s inequality and obtain

 (2.21) σ(∥utt∥2+κ∥∇utt∥2)≤C(ν)σ(∥~Δut∥2+∥∇ut∥2+∥~Δu∥2+∥ft∥2).

An integration of (2.21) with respect to time from to , a multiplication by and a use of (2.12), Lemma 2.1, assumption (A2) complete the proof of (2.13).
Now to derive (2.14), substitute in (2.5) and use Cauchy-Schwarz’s inequality, Young’s inequality to yield

 (2.22) ∥∇utt∥2+κ∥~Δutt∥2≤C(ν)κ(∥~Δut∥2+∥∇ut∥2+∥~Δu∥2+∥ft∥2).

Multiply (2.22) by and integrate the resulting equation with respect to time from to to obtain

 (2.23) ∫t0σ(s)(κ∥∇uss(s)∥2+κ2∥~Δuss(s)∥2)ds≤C∫t0σ(s)(ν∥~Δus(s)∥2+∥∇us(s)∥2+∥~Δu(s)∥2+∥fs(s)∥2)ds.

Multiply (2.23) by and use (2.12), Lemma 2.1 to arrive at the desired result in (2.14).
Now to prove pressure estimate in (2.12), rewrite (2.4). Then, differentiate the resulting equation with respect to time and obtain

 (2.24) (∇pt,\boldmathϕ)=(utt,% \boldmathϕ)−κ(~Δutt,\boldmathϕ)−ν(~Δut,\boldmathϕ)+(ut⋅∇u,\boldmathϕ)+(u⋅∇ut,\boldmathϕ)−(ft,\boldmathϕ).

Choose in (2.24). Then, apply Cauchy-Schwarz’s inequality and generalized Hölder’s inequality to find that

 (2.25) ∥∇pt∥≤C(∥utt∥+κ∥~Δutt∥+ν∥~Δut∥+∥∇ut∥+∥~Δu∥+∥ft∥).

After squaring both sides of (2.25), multiply it by and integrate with respect to time from to to arrive at

 ∫t0σ(s)∥∇ps(s)∥2ds≤C∫t0 σ(s)(∥uss(s)∥2+κ2∥~Δuss(s)∥2+ν∥~Δus(s)∥2 (2.26) +∥∇us(s)∥2+∥~Δu(s)∥2+∥fs(s)∥2)ds.

A use of estimates of from (2.12), (2.13), Lemma 2.1, assumption (A2) and a multiplication by lead to

 (2.27) e−2αt∫t0σ(s)∥∇ps(s)∥2ds≤C.

This completes the proof of Lemma 2.2.
To derive uniform estimates in time, we assume the following uniqueness condition:

 (2.28) Nν2∥f∥L∞(0,∞;L2)<1andN=supu,v,w∈H10(Ω)b(u,v,w)∥∇u∥∥∇v∥∥∇w∥.

## 3 Semidiscrete Approximation

Let and be the finite-dimensional subspaces of and , respectively, such that, there exist operators and satisfying the following approximation properties:
(B1). For each and , there exist approximations and such that

 ∥w−ihw∥+h∥∇(w−ihw)∥≤K0h2∥w∥2,∥q−jhq∥L2/IR≤K0h∥q∥H1/IR.

Note that, be a discretization parameter with .
Here, it can be noted that the operator preserves the antisymmetric properties of the original nonlinear term, i.e.,

 (3.1) b(vh,wh,wh)=0∀vh,wh∈Hh.

The discrete analogue of the weak formulation (2.4) is as follows:
Find and such that and for ,

 (uht,\boldmathϕh)+κa(uht,\boldmathϕh)+νa(uh,\boldmathϕh) +b(uh,uh,\boldmathϕh)−(ph,∇⋅\boldmathϕh)=0∀\boldmathϕh∈Hh, (3.2) (∇⋅uh,χh) =0∀χh∈Lh,

where is a suitable approximation of .
For subsequent analysis, we define a suitable approximation of by introducing the discrete incompressibility condition in and call the resulting subspace as . Thus, is defined as

 Jh={vh∈Hh:(χh,∇⋅vh)=0∀χh∈Lh}.

Note that, the space is not a subspace of . Now, an equivalent form of (3) is defined as:
Find such that and for ,

 (3.3) (uht,\boldmathϕh)+κa(uht,\boldmathϕh)+νa(uh,\boldmathϕh)+b(uh,uh,\boldmathϕh)=0∀% \boldmathϕh∈Jh.

For proof of the global existence of a unique solution of (3.3), one may refer to [2].
In order to deal with the pressure estimates in subsequent analysis, we assume the pair satisfies a uniform inf-sup condition:
(B2). For every , there exist a non-trivial function and a positive constant , independent of , such that,

 |(qh,∇⋅\boldmathϕh)|≥K1∥∇\boldmathϕh∥∥qh∥L2/Nh.

The following properties of the projection can be derived using conditions (B1)-(B2) ( for a proof, see ([11], [13]):

 (3.4) ∥\boldmathϕ−Ph\boldmathϕ∥+h∥∇Ph\boldmathϕ∥≤Ch∥∇\boldmathϕ∥∀\boldmathϕ∈J1,

and

 (3.5) ∥\boldmathϕ−Ph\boldmathϕ∥+h∥∇(\boldmathϕ−Ph\boldmathϕ)∥≤Ch2∥~Δ\boldmathϕ∥∀\boldmathϕ∈J1∩H2.

We may define the discrete operator through the bilinear form as

 (3.6) a(vh,\boldmathϕh)=(−Δhvh,\boldmathϕh)∀vh,\boldmathϕh∈Hh.

Set the discrete analogue of the Stokes operator as . Examples of subspaces and satisfying assumptions and in the context of both conforming and non-conforming analysis can be found in [6], [7] and [13].
We recall below in Lemma 3.1, some a priori bounds of which will be used in the derivation of fully discrete error estimates in the subsequent section. For proof, one may refer to [26] (Lemma 4.2), [28] (Lemma 3.2).

###### Lemma 3.1.

Let the assumptions (A1)-(A2) hold. Then, there exists a positive constant such that for the following estimates hold true:

 ∥uht(t)∥2+∥~Δhuh(t)∥2+e−2αt∫t0e2αs(∥∇uh(s)∥2+∥~Δhuh(s)∥2+∥∇uhs(s)∥2)ds≤C, e−2αt∫t0e2αs(∥uhss(s)∥2−1+κ∥∇uhss(s)∥2)ds≤C.

Now, in Theorem 3.1, the main results of the article are stated in which we present the semidiscrete optimal error estimates of the velocity and pressure. The proofs are established in Sections .

###### Theorem 3.1.

Let the assumptions (A1)-(A2) and (B1)-(B2) be satisfied. Let , then, there exists a positive constant depending on , , , and , such that, for fixed with and for , the following estimates hold true:

 ∥(u−uh)(t)∥+h∥∇(u−uh)(t)∥≤K(t)h2,

. Under the uniqueness condition (2.28), , that is, the estimates are uniform in time.

###### Theorem 3.2.

Under the hypotheses of Theorem 3.1, there exists a positive constant depending on , , , and , such that, for all , the following holds true:

 ∥(p−ph)(t)∥L2/Nh≤K(t)h.

Here again, under the uniqueness condition (2.28), , that is, the estimate holds uniformly with respect to time.

## 4 Semidiscrete Finite Element Error Estimates

This section deals with the optimal error estimates of velocity and pressure. Note that, since is not a subspace of , the weak solution satisfies

 (4.1) (ut,\boldmathϕh)+κa(ut,\boldmathϕh)+νa(u,\boldmathϕh)=−b(u,u,\boldmathϕh)+(p,∇⋅\boldmathϕh)+(f,\boldmathϕh)∀% \boldmathϕh∈Jh.

Set . Then, subtract (4.1) from (3.3) to arrive at

 (4.2) (et,\boldmathϕh)+κa(et,\boldmathϕh)+νa(e,\boldmathϕh)=Λ(\boldmathϕh)+(p,∇⋅\boldmathϕh),

where . Below, we derive the optimal error estimates of and , for .
In order to deal with the nonlinearity, an intermediate solution is introduced which is a finite element Galerkin approximation to a linearized Kelvin-Voigt equation. The solution satisfies

 (4.3) (vht,\boldmathϕh)+κa(vht,\boldmathϕh)+νa(vh,\boldmathϕh)=−b(u,u,\boldmathϕh)∀% \boldmathϕ∈Jh,

with