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

TL;DR

None

Contribution

None

Abstract

In this article, a finite element Galerkin method is applied to the Kelvin-Voigt viscoelastic fluid model, when its forcing function is in $L^{\infty}(\bL^2)$. Some new {\it a priori} bounds for the velocity as well as for the pressure are derived which are independent of inverse powers of the retardation time $\kappa$. Optimal error estimates for the velocity in $L^{\infty} (\bL^2)$ as well as in $L^{\infty}(\bH^1_0)$-norms and for the pressure in $L^{\infty}(L^2)$-norm of the semidiscrete method are discussed which hold uniformly with respect to $\kappa$ as $\kappa\rightarrow 0$ with the initial condition only in $\bH^2\cap\bH_0^1$. Further, under uniqueness condition, these estimates are shown to be uniformly in time as $t \mapsto \infty$. 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 $\kappa\rightarrow 0$.

Tables

Table 1. Table 1: Convergence rates for backward Euler method with k = 𝕆 ​ ( h 2 ) 𝑘 𝕆 superscript ℎ 2 k=\mathbb{O}(h^{2}) and κ → 0 → 𝜅 0 \kappa\rightarrow 0 .

$h$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$
	$\kappa =1$	$\kappa ={10}^{-3}$	$\kappa ={10}^{-6}$	$\kappa ={10}^{-9}$
1/4	1.356182	1.575487	1.575615	1.575615
1/8	1.721053	1.784026	1.784008	1.784008
1/16	1.879115	1.895495	1.895489	1.895489
1/32	1.946837	1.950956	1.950955	1.950955

Table 2. Table 2: Convergence rates for backward Euler method with k = 𝕆 ​ ( h ) 𝑘 𝕆 ℎ k=\mathbb{O}(h) and κ → 0 → 𝜅 0 \kappa\rightarrow 0 .

$h$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐇}^1}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐇}^1}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐇}^1}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐇}^1}$
	$\kappa =1$	$\kappa ={10}^{-3}$	$\kappa ={10}^{-6}$	$\kappa ={10}^{-9}$
1/4	0.666127	0.810563	0.810686	0.810686
1/8	0.856391	0.903887	0.903874	0.903874
1/16	0.939755	0.952619	0.952616	0.952616
1/32	0.973984	0.977304	0.977303	0.977303

Table 3. Table 3: Convergence rates for backward Euler method with k = 𝕆 ​ ( h 2 ) 𝑘 𝕆 superscript ℎ 2 k=\mathbb{O}(h^{2}) and κ → 0 → 𝜅 0 \kappa\rightarrow 0 .

$h$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$
	$\kappa =1$	$\kappa ={10}^{-3}$	$\kappa ={10}^{-6}$	$\kappa ={10}^{-9}$
1/4	0.935384	0.934915	0.934891	0.934891
1/8	0.962732	0.960305	0.960305	0.960305
1/16	0.982119	0.980954	0.980955	0.980955
1/32	0.991049	0.990312	0.990312	0.990312

Table 4. Table 4: Regularization effect on velocity errors for backward Euler method with ν = 0.01 𝜈 0.01 \nu=0.01 , k = 𝕆 ​ ( h 2 ) 𝑘 𝕆 superscript ℎ 2 k=\mathbb{O}(h^{2}) .

$h$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert 𝐮\left(t_n\right)-{𝐔}^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$
	$\kappa =1$	$\kappa ={10}^{-2}$	$\kappa ={10}^{-4}$	$\kappa =0$
1/2	0.156344	5.977491	29.31149	7.060819
1/4	0.083167	3.223425	26.16600	35.799670
1/8	0.027082	1.425267	1.810881	1.798025
1/16	0.007577	0.432046	0.500769	0.494285
1/32	0.001991	0.116965	0.130588	0.128835

Table 5. Table 5: Regularization effect on pressure errors for backward Euler method with ν = 0.01 𝜈 0.01 \nu=0.01 , k = 𝕆 ​ ( h 2 ) 𝑘 𝕆 superscript ℎ 2 k=\mathbb{O}(h^{2}) .

$h$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$	${\Vert p\left(t_n\right)-P^n\Vert }_{{𝐋}^{\mathrm{𝟐}}}$
	$\kappa =1$	$\kappa ={10}^{-2}$	$\kappa ={10}^{-4}$	$\kappa =0$
1/2	3.057302	7.657372	160.344805	885.464339
1/4	1.290707	3.121213	204.368122	419.021951
1/8	0.589169	1.386778	1.1096204	1.089742
1/16	0.280954	0.327177	0.2832502	0.283915
1/32	0.137575	0.139234	0.1419708	0.142181

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}-\kappa \Delta {\mathbf{u}}_t-\nu \Delta \mathbf{u}+\nabla p=\mathbf{f}(x,t),x\in \Omega ,t>0,$$

$$\nabla \cdot \mathbf{u}=0,x\in \Omega ,t>0,$$

$$\mathbf{u}(x,0)={\mathbf{u}}_0\text{mbox}in\Omega ,\mathbf{u}=0,\text{mbox}on\partial \Omega ,t\ge 0.$$

$$\mathbf{J}$$

$${\mathbf{J}}_1$$

$$\Vert \mathbf{v}{\Vert }_2\le C\Vert \tilde{\Delta }\mathbf{v}\Vert \forall \mathbf{v}\in {\mathbf{J}}_1\cap {\mathbf{H}}^2.$$

$$\Vert \mathbf{v}{\Vert }^2\le {\lambda }_1^{-1}\Vert \nabla \mathbf{v}{\Vert }^2\forall \mathbf{v}\in {\mathbf{H}}_0^1(\Omega ),$$

$$\Vert \nabla \mathbf{v}{\Vert }^2\le {\lambda }_1^{-1}\Vert \tilde{\Delta }\mathbf{v}{\Vert }^2\forall \mathbf{v}\in {\mathbf{J}}_1\cap {\mathbf{H}}^2.$$

$${\mathbf{u}}_0\in {\mathbf{H}}^2\cap {\mathbf{J}}_1,\mathbf{f},{\mathbf{f}}_t\in L^{\infty }(0,T;{\mathbf{L}}^2)\mathrm{with}\Vert {\mathbf{u}}_0{\Vert }_2\le \mathrm{M},{\mathrm{esssup}}_{0<\mathrm{t}\le \mathrm{T}}\{\Vert \mathbf{f}(\cdot ,\mathrm{t})\Vert ,\Vert {\mathbf{f}}_{\mathrm{t}}(\cdot ,\mathrm{t})\Vert \}\le M.$$

$$\displaystyle\left.\begin{array}[]{rcl}&&({\bf u}_{t},\mbox{\boldmath$\phi$})+\kappa\,(\nabla{\bf u}_{t},\nabla\mbox{\boldmath$\phi$})+\nu\,(\nabla{\bf u},\nabla\mbox{\boldmath$\phi$})+({\bf u}\cdot\nabla{\bf u},\mbox{\boldmath$\phi$})-(p,\nabla\cdot\mbox{\boldmath$\phi$})=({\bf f},\mbox{\boldmath$\phi$})\;\;\;\;\;\forall\mbox{\boldmath$\phi$}\in{\bf H}_{0}^{1},\\ &&(\nabla\cdot{\bf u},\chi)=0\;\;\;\forall\chi\in L^{2}.\end{array}\right\}$$

$$\displaystyle({\bf u}_{t},\mbox{\boldmath$\phi$})+\kappa\,(\nabla{\bf u}_{t},\nabla\mbox{\boldmath$\phi$})+\nu\,(\nabla{\bf u},\nabla\mbox{\boldmath$\phi$})+({\bf u}\cdot\nabla{\bf u},\mbox{\boldmath$\phi$})=({\bf f},\mbox{\boldmath$\phi$})\;\;\;\forall\mbox{\boldmath$\phi$}\in{\bf J}_{1}.$$

$$a({\bf v},\mbox{\boldmath$\phi$})=(\nabla{\bf v},\nabla\mbox{\boldmath$\phi$})$$

$$b({\bf v},{\bf w},\mbox{\boldmath$\phi$})=\frac{1}{2}({\bf v}\cdot\nabla{\bf w},\mbox{\boldmath$\phi$})-\frac{1}{2}({\bf v}\cdot\nabla\mbox{\boldmath$\phi$},{\bf w}).$$

$$b(\mathbf{v},\mathbf{w},\mathbf{w})=0\forall \mathbf{v},\mathbf{w}\in {\mathbf{J}}_1.$$

$$\underset{0<t<\infty }{\sup }\{\Vert \mathbf{u}(t){\Vert }_2+\Vert {\mathbf{u}}_t(t)\Vert +\kappa \Vert \tilde{\Delta }{\mathbf{u}}_t(t)\Vert +\Vert p(t){\Vert }_{{\mathbf{H}}^1/\mathbb{R}}\}\le C,$$

$$e^{-2\alpha t}{\int }_0^t{e}^{2\alpha s}\Vert {\mathbf{u}}_s(s){\Vert }_1^{2}ds\le C,$$

$${\sigma }^{-1}(t){\int }_0^t{e}^{2\alpha s}(\Vert \mathbf{u}(s){\Vert }_2^2+\Vert p(s){\Vert }_{{\mathbf{H}}^1/\mathbb{R}}^2+\kappa \Vert \tilde{\Delta }{\mathbf{u}}_s(s){\Vert }^2)ds\le C,$$

$$\frac{1}{\sigma }{\int }_0^t{e}^{2\alpha s}\Vert \tilde{\Delta }\mathbf{u}(s){\Vert }^{2}ds$$

$$\displaystyle\leq\frac{C}{te^{2\alpha t}}\bigg{(}\frac{e^{2\alpha t}-1}{2\alpha}\bigg{)}$$

$$\displaystyle=\frac{C}{2\alpha te^{2\alpha t}}\bigg{(}\frac{2\alpha t}{1!}+\frac{(2\alpha t)^{2}}{2!}+\frac{(2\alpha t)^{3}}{3!}+\cdots\bigg{)}$$

$$\displaystyle\leq\frac{C}{e^{2\alpha t}}\,\bigg{(}1+\frac{2\alpha t}{1!}+\frac{(2\alpha t)^{2}}{2!}+\cdots\bigg{)}\leq\frac{C}{e^{2\alpha t}}\,e^{2\alpha t}\leq C.$$

$${\sigma }^{-1}(t){\int }_0^t{e}^{2\alpha s}\Vert \mathbf{u}(s){\Vert }_2^{2}ds\le C.$$

$${\sigma }^{-1}(t){\int }_0^t{e}^{2\alpha s}(\Vert p(s){\Vert }_{{\mathbf{H}}^1/\mathbb{R}}^2+\kappa \Vert \tilde{\Delta }{\mathbf{u}}_s(s){\Vert }^2)ds\le C$$

$$\displaystyle\tau(t)\big{(}\|\nabla{\bf u}_{t}\|^{2}+\kappa\,\|\tilde{\Delta}{\bf u}_{t}\|^{2}\big{)}+\nu\,e^{-2\alpha t}\,\displaystyle{\int_{0}^{t}}{\bf\sigma}(s)\,(\|\tilde{\Delta}{\bf u}_{s}(s)\|^{2}+\|\nabla p_{s}(s)\|^{2})ds\leq C,$$

$$e^{-2\alpha t}{\int }_0^t\sigma (s)(\Vert {\mathbf{u}}_{ss}(s){\Vert }^2+\kappa \Vert \nabla {\mathbf{u}}_{ss}(s){\Vert }^2)ds\le C,$$

$$e^{-2\alpha t}{\int }_0^t\sigma (s)(\kappa \Vert \nabla {\mathbf{u}}_{ss}(s){\Vert }^2+{\kappa }^2\Vert \tilde{\Delta }{\mathbf{u}}_{ss}(s){\Vert }^2)ds\le C.$$

$$\displaystyle({\bf u}_{tt},\mbox{\boldmath$\phi$})-\kappa\,(\tilde{\Delta}{\bf u}_{tt},\mbox{\boldmath$\phi$})-\nu\,(\tilde{\Delta}{\bf u}_{t},\mbox{\boldmath$\phi$})+({\bf u}_{t}\cdot\nabla{\bf u}+{\bf u}\cdot\nabla{\bf u}_{t},\mbox{\boldmath$\phi$})=({\bf{f}}_{t},\mbox{\boldmath$\phi$}),\forall\mbox{\boldmath$\phi$}\in{\bf J}_{1}.$$

$$\frac{1}{2}\frac{d}{dt}\sigma \left(\Vert \nabla {\mathbf{u}}_t{\Vert }^2+\kappa \Vert \tilde{\Delta }{\mathbf{u}}_t{\Vert }^2\right)$$

$$+\sigma ({\mathbf{u}}_t\cdot \nabla \mathbf{u},\tilde{\Delta }{\mathbf{u}}_t)+\sigma (\mathbf{u}\cdot \nabla {\mathbf{u}}_t,\tilde{\Delta }{\mathbf{u}}_t)-\sigma ({\mathbf{f}}_t,\tilde{\Delta }{\mathbf{u}}_t)$$

$$\displaystyle=C(\alpha,\lambda_{1})e^{2\alpha t}\big{(}\|\nabla{\bf u}_{t}\|^{2}+\kappa\,\|\tilde{\Delta}{\bf u}_{t}\|^{2}\big{)}+I_{1}+I_{2}+I_{3},(\mbox{say}).$$

$$|I_1|+|I_2|$$

$$\displaystyle\leq C(\epsilon){\bf\sigma}\,\big{(}\|{\bf u}_{t}\|\|\nabla{\bf u}_{t}\|\|\nabla{\bf u}\|\|\tilde{\Delta}{\bf u}\|+\|\nabla{\bf u}\|\|\tilde{\Delta}{\bf u}\|\|\nabla{\bf u}_{t}\|^{2}\big{)}+\epsilon\,{\bf\sigma}\,\|\tilde{\Delta}{\bf u}_{t}\|^{2}$$

$$\le C(ϵ,{\lambda }_1)\sigma \Vert \nabla {\mathbf{u}}_t{\Vert }^2\Vert \tilde{\Delta }\mathbf{u}{\Vert }^2+ϵ\sigma \Vert \tilde{\Delta }{\mathbf{u}}_t{\Vert }^2.$$

$$|I_3|\le C\sigma \Vert {\mathbf{f}}_t\Vert \Vert \tilde{\Delta }{\mathbf{u}}_t\Vert \le C(ϵ)\sigma \Vert {\mathbf{f}}_t{\Vert }^2+ϵ\sigma \Vert \tilde{\Delta }{\mathbf{u}}_t{\Vert }^2.$$

Full text

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 :

[TABLE]

and incompressibility condition

[TABLE]

with initial and boundary conditions

[TABLE]

Here, $\Omega$ is a bounded convex polygonal or polyhedral domain in ${{\rm I{\!}\rm R}}^{d},d=2,3$ with boundary $\partial\Omega$, ${\bf u}=(u_{1},u_{2})$ $(\mbox{or}~{}{\bf u}=(u_{1},u_{2},u_{3}))$ represents the velocity vector, $p$ is the pressure of the fluid, ${\bf{f}}$ is the external force, $\nu>0$ denotes the kinematic coefficient of viscosity and $\kappa$ 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 ${\bf{f}}$ and for all time $t>0$.

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 $t\geq 0$ with the assumption that the solution is asymptotically stable as $t\rightarrow\infty$. Further, the author established optimal error estimates in $L^{\infty}({{\bf{H}}_{0}^{1}})$-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 $L^{\infty}(\bf{L}^{2})$ and $L^{\infty}({\bf{H}}_{0}^{1})$-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 $L^{\infty}({\bf{H}}_{0}^{1})$-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 $\kappa^{-r}$, for $r\geq 2$ which may blow as $\kappa\rightarrow 0$.

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 $\kappa$ 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 $L^{\infty}(\bf{L}^{2})$ as well as in $L^{\infty}({\bf{H}}_{0}^{1})$- norms and for the pressure in $L^{\infty}(L^{2})$-norm, when the forcing function ${\bf{f}}\in L^{\infty}(\bf{L}^{2})$. Here, it can be noted that they have achieved an improvement in the error estimates in powers of $\kappa$ as the constants in error bounds depend only on $\kappa^{-1/2}$.

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 ${\bf{f}}\in L^{\infty}(\bf{L}^{2})$ 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, ${\bf u}_{0}\in{\bf{H}}^{3}\cap{{\bf{H}}_{0}^{1}}$ with $\Delta{\bf u}_{0}=0$ on $\partial\Omega$, the optimal error estimates independent of $\kappa$ 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 ${\bf{f}}=0$, 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 $\exp(-\kappa^{-1}t)$ [19], it is worthwhile to discuss the behavior of the solution as $\kappa\rightarrow 0$ and as $t\rightarrow\infty$. Moreover, this model can be thought of as a $\kappa$ 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 $\kappa$ under realistically assumed minimum regularity assumption on the exact solution with ${\bf u}_{0}\in{\bf{H}}^{2}\cap{\bf{H}}_{0}^{1}$ and ${\bf{f}}$, ${\bf{f}}_{t}\in L^{\infty}(\bf{L}^{2})$.

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 $\kappa$ as $\kappa\rightarrow 0$ under minimum regularity assumptions ${\bf u}_{0}\in{\bf{H}}^{2}\cap{\bf{H}}_{0}^{1}$ and ${\bf{f}}$, ${\bf{f}}_{t}\in L^{\infty}(\bf{L}^{2})$. Here, it can be noted that the introduction of weight function ${\bf\sigma}(t)=\min\{1,\,t\}e^{2\alpha t}$ plays a key role in handling the regularity issues at $t=0$.

(ii) Using the Sobolev-Stokes projection defined earlier in [2], fully discrete optimal error estimates in $L^{\infty}(\bf{L}^{2})$ and $L^{\infty}({\bf{H}}_{0}^{1})$-norms for the finite element velocity approximation and in $L^{\infty}(L^{2})$-norm for the finite element pressure approximation are established. It is further proved that these error estimates hold uniformly as $\kappa\rightarrow 0$. 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 ${\bf\sigma}(t)$ which plays an important role in achieving uniform estimates in terms of $\kappa$.

(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 $\kappa\rightarrow 0$ confirming the results in $(ii).$

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 $\kappa$. As a consequence, the error estimates do not blow up as $\kappa\rightarrow 0$. The main difficulty in making error estimates independent of $\kappa$ arises due to the lack of regularity of solution at $t=0$. In order to overcome this difficulty, we introduce various powers of weight function ${\bf\sigma}(t)$ which takes care of regularity issues of the solution at $t=0$.

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 $\mathbb{R}^{d}$ $(d=2,3)$-valued function spaces using bold face letters, that is, ${\bf H}_{0}^{1}=(H_{0}^{1}(\Omega))^{d}$, ${\bf L}^{2}=(L^{2}(\Omega))^{d}$ and ${\bf H}^{m}=(H^{m}(\Omega))^{d},$ where ${\bf{L}}^{2}(\Omega)$ is the space of square integrable functions defined in $\Omega$ with inner product (\mbox{\boldmath\phi},{\bf\psi})=\displaystyle{\int_{0}^{t}}\mbox{\boldmath\phi}({\bf x}){\bf\psi}({\bf x})dx and norm \|\mbox{\boldmath\phi}\|=\left(\displaystyle{\int_{0}^{t}}|\mbox{\boldmath\phi}({\bf x})|^{2}dx\right)^{1/2}. Further, ${\bf{H}}^{m}(\Omega)$ denotes the standard Hilbert Sobolev space of order $m\in N^{+}$ with norm \|\mbox{\boldmath\phi}\|_{m}=\displaystyle{\sum_{|\alpha|\leq m}}\left(\displaystyle{\int_{0}^{t}}|D^{\alpha}\mbox{\boldmath\phi}|^{2}dx\right)^{1/2}. The space ${\bf{H}}_{0}^{1}$ is equipped with a norm $\|\nabla{\bf v}\|=\left({\sum_{i,j=1}^{d}}(\partial_{j}v_{i},\partial_{j}v_{i})\right)^{1/2}=\left({\sum_{i=1}^{d}}(\nabla v_{i},\nabla v_{i})\right)^{1/2}$. Given a Banach space $X$ endowed with norm ${\parallel\cdot\parallel}_{X}$, let $L^{p}(0,T;X)$ be the space of all strongly measurable functions $\phi:(0,T)\rightarrow X$ satisfying $\displaystyle{\int_{0}^{T}}{\parallel\phi(s)\parallel}^{p}_{X}\,ds<\infty$ if $1\leq p<\infty$ and for $p=\infty$, \rm{\displaystyle{\mathop{ess\sup}_{t\in(0,T)}}{\parallel\mbox{\boldmath\phi}(t)\parallel}_{X}<\infty}. Also, we define the divergence free spaces

[TABLE]

where ${\bf{n}}$ is the unit outward normal to the boundary $\partial\Omega$ and \mbox{\boldmath\phi}\cdot{\bf{n}}|_{\partial\Omega}=0 should be understood in the sense of trace in ${\bf{H}}^{-1/2}(\partial\Omega)$, see [29]. Let $H^{m}/{\rm I\!R}$ be the quotient space with norm $\|\phi\|_{H^{m}/{\rm I\!R}}=\inf_{c\in{\rm I\!R}}\|\phi+c\|_{m}$. For $m=0$, it is denoted by $L^{2}/{\rm I\!R}$. Now, define $P:\bf{L}^{2}\rightarrow{\bf J}$ as the $\bf{L}^{2}$-orthogonal projection.

Throughout this article, we make the following assumptions:

(A1). Setting $-\tilde{\Delta}=-P\Delta:{\bf J}_{1}\cap{\bf H}^{2}\subset{\bf J}\rightarrow{\bf J}$ as the Stokes operator, assume that the following regularity result holds:

[TABLE]

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

[TABLE]

where $\lambda_{1}^{-1}$, is the best possible positive constant depending on the domain $\Omega.$ Further, observe (see, [13]) that

[TABLE]

(A2). There exists a positive constant $M$ such that the initial velocity ${\bf u}_{0}$ and the external force ${\bf{f}}$ satisfy for $t\in(0,\,T]$ with $0<T<\infty$

[TABLE]

The weak formulation of (1.1)-(1.3) is to find $({\bf u}(t),p(t))\in{\bf{H}}_{0}^{1}\times L^{2}/{\rm I\!R}$, such that ${\bf u}(0)={\bf u}_{0}$ and for $t>0$

[TABLE]

Equivalently, find ${\bf u}(t)\in{\bf J}_{1}$ such that for ${\bf u}(0)={\bf u}_{0}$, $t>0$,

[TABLE]

For {\bf v},{\bf w},\mbox{\boldmath\phi}\in{\bf{H}}_{0}^{1}, define the bilinear form $a(\cdot,\cdot)$ as

[TABLE]

and the trilinear form $b(\cdot,\cdot,\cdot)$ as

[TABLE]

Note that, for ${\bf v}\in{\bf J}_{1}$, ${\bf w}$, \mbox{\boldmath\phi}\in{\bf{H}}_{0}^{1}, b({\bf v},{\bf w},\mbox{\boldmath\phi})=({\bf v}\cdot\nabla{\bf w},\mbox{\boldmath\phi}). Because of antisymmetric property of the trilinear form, it is easy to verify that

[TABLE]

We present below in Lemma 2.1, some a priori bounds for the weak solution pair $({\bf u},\,p)$ which will be used in our subsequent error analysis. Since the estimates in (2.8) and (2.9) are already derived in [26], we only provide proof of (2.10).

Lemma 2.1**.**

[[26], pp 241, 244] Let the assumptions (A1)-(A2)* hold. Then, there exists a positive constant $C=C(\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\lambda_{1}\kappa\big{)}}} the following estimates hold true:*

[TABLE]

where $\tau(t):=\min\{t,\,1\}$ and ${\bf\sigma}(t):=\tau(t)e^{2\alpha t}$.

Proof. We know that ${\bf\sigma}(t)=\tau(t)\,e^{2\alpha t}$, where $\tau(t)=\min\{t,1\}=t\,\text{or}\,1$. Hence, ${\bf\sigma}(t)=e^{2\alpha t}$ or ${\bf\sigma}(t)=t\,e^{2\alpha t}$.

Now, consider the following two cases:

Case 1: ${\bf\sigma}(t)=e^{2\alpha t}$. Then,

[TABLE]

A use of (2.8) in (2.11) leads to

[TABLE]

Case 2: ${\bf\sigma}(t)=t\,e^{2\alpha t}$. Again use (2.8) and well known facts of series to obtain

[TABLE]

Therefore, considering the above two cases, we arrive at

[TABLE]

Following the similar sets of arguments as above, we obtain

[TABLE]

and this completes the remaining part of the proof. $\Box$

In the next lemma, we derive a priori bounds for the highest order time derivatives of weak solution for the problem (2.6).

Lemma 2.2**.**

Let the assumptions (A1)-(A2)* hold. Then, there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\lambda_{1}\kappa\big{)}}} the following estimates hold true:*

[TABLE]

Note that, here and everywhere else in the consecutive analysis the constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ is independent of inverse powers of $\kappa$.

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

[TABLE]

Choose $\phi=-{\bf\sigma}(t)\tilde{\Delta}{\bf u}_{t}$ in (2.17) to obtain

[TABLE]

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

[TABLE]

Once again, apply Cauchy-Schwarz’s inequality and Young’s inequality to bound $|I_{3}|$ as

[TABLE]

After using (2)-(2.20) in (2) with a proper choice of $\epsilon$, integrate the resulting equation with respect to time from [math] to $t$ to arrive at

[TABLE]

Apply Lemma 2.1 and assumption (A2) in (2). Then, multiply the resulting equation by $e^{-2\alpha t}$ to arrive at the desired a priori estimates of ${\bf u}$ in (2.14).

Next, differentiate (2.7) and substitute \mbox{\boldmath\phi}={\bf\sigma}{\bf u}_{tt} in the resulting equation to observe that

[TABLE]

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

[TABLE]

An integration of (2.23) with respect to time from [math] to $t$, a multiplication by $e^{-2\alpha t}$ and a use of (2.14), Lemma 2.1, assumption (A2) complete the proof of (2.15).

Now to derive (2.16), substitute \mbox{\boldmath\phi}=-\tilde{\Delta}{\bf u}_{tt} in (2.7) and use Cauchy-Schwarz’s inequality, Young’s inequality to yield

[TABLE]

Multiply (2.24) by $\kappa\,{\bf\sigma}$ and integrate the resulting equation with respect to time from [math] to $t$ to obtain

[TABLE]

Multiply (2.25) by $e^{-2\alpha t}$ and use (2.14), Lemma 2.1 to arrive at the desired result in (2.16).

Now to prove pressure estimate in (2.14), rewrite (2.6). Then, differentiate the resulting equation with respect to time and obtain

[TABLE]

Choose \mbox{\boldmath\phi}=\nabla p_{t} in (2.26). Then, apply Cauchy-Schwarz’s inequality and generalized Hölder’s inequality to find that

[TABLE]

After squaring both sides of (2.27), multiply it by ${\bf\sigma}$ and integrate with respect to time from [math] to $t$ to arrive at

[TABLE]

A use of estimates of ${\bf u}$ from (2.14), (2.15), Lemma 2.1, assumption (A2) and a multiplication by $e^{-2\alpha t}$ lead to

[TABLE]

This completes the proof of Lemma 2.2. $\Box$

To derive uniform estimates in time, we assume the following uniqueness condition:

[TABLE]

3 Semidiscrete Approximation

Let ${\bf H}_{h}$ and $L_{h}$ be the finite-dimensional subspaces of ${\bf H}_{0}^{1}$ and $L^{2}$, respectively, such that, there exist operators $i_{h}$ and $j_{h}$ satisfying the following approximation properties:

(B1). For each ${\bf w}\in{\bf{J}}_{1}\cap{\bf{H}}^{2}$ and $q\in H^{1}/{\rm I\!R}$, there exist approximations $i_{h}{\bf w}\in{\bf{J}}_{h}$ and $j_{h}q\in L_{h}$ such that

[TABLE]

Note that, $h>0$ be a discretization parameter with $0<h<1$.

Here, it can be noted that the operator $b(\cdot,\cdot,\cdot)$ preserves the antisymmetric properties of the original nonlinear term, i.e.,

[TABLE]

The discrete analogue of the weak formulation (2.6) is as follows:

Find ${\bf u}_{h}(t)\in{\bf H}_{h}$ and $p_{h}(t)\in L_{h}$ such that ${\bf u}_{h}(0)={\bf u}_{0h}$ and for $t>0$,

[TABLE]

where ${\bf u}_{0h}\in{\bf H}_{h}$ is a suitable approximation of ${\bf u}_{0}\in{\bf J}_{1}$.

For subsequent analysis, we define a suitable approximation of ${\bf J}_{1}$ by introducing the discrete incompressibility condition in ${\bf{H}}_{h}$ and call the resulting subspace as ${\bf J}_{h}$. Thus, ${\bf J}_{h}$ is defined as

[TABLE]

Note that, the space ${\bf J}_{h}$ is not a subspace of ${\bf J}_{1}$. Now, an equivalent form of (3) is defined as:

Find ${\bf u}_{h}(t)\in{\bf J}_{h}$ such that ${\bf u}_{h}(0)={\bf u}_{0h}$ and for $t>0$,

[TABLE]

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 $({\bf H}_{h},{L_{h}/N_{h}})$ satisfies a uniform inf-sup condition:

(B2). For every $q_{h}\in L_{h}$, there exist a non-trivial function \mbox{\boldmath\phi}_{h}\in\bf H_{h} and a positive constant $K_{1}$, independent of $h$, such that,

[TABLE]

The following properties of the $L^{2}$ projection $P_{h}:\bf{L}^{2}\rightarrow{\bf J}_{h}$ can be derived using conditions (B1)-(B2) ( for a proof, see ([11], [13]):

[TABLE]

and

[TABLE]

We may define the discrete operator $\Delta_{h}:{\bf{H}}_{h}\rightarrow{\bf{H}}_{h}$ through the bilinear form $a(\cdot,\cdot)$ as

[TABLE]

Set the discrete analogue of the Stokes operator $\tilde{\Delta}=P\Delta$ as $\tilde{\Delta}_{h}=P_{h}\Delta_{h}$. Examples of subspaces $\bf H_{h}$ and $L_{h}$ satisfying assumptions $(\rm{\bf{B1}})$ and $(\rm{\bf{B2}})$ 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 ${\bf u}_{h}$ 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 $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\lambda_{1}\kappa\big{)}}} the following estimates hold true:*

[TABLE]

$\Box$

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 ${\bf 5}$.

Theorem 3.1**.**

Let the assumptions (A1)-(A2)* and (B1)-(B2) be satisfied. Let ${\bf u}_{0h}=P_{h}{\bf u}_{0}$, then, there exists a positive constant $C$ depending on $\kappa$, $\lambda_{1}$, $\nu$, $\alpha$ and $M$, such that, for fixed $T>0$ with $t\in(0,\,T)$ and for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\lambda_{1}\kappa\big{)}}}, the following estimates hold true:*

[TABLE]

$K(t)=Ce^{Ct}$. Under the uniqueness condition (2.30), $K(t)=C$, that is, the estimates are uniform in time.

Theorem 3.2**.**

Under the hypotheses of Theorem 3.1, there exists a positive constant $C$ depending on $\kappa$, $\lambda_{1}$, $\nu$, $\alpha$ and $M$, such that, for all $t>0$, the following holds true:

[TABLE]

Here again, under the uniqueness condition (2.30), $K(t)=C$, 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 ${\bf J}_{h}$ is not a subspace of ${\bf J}_{1}$, the weak solution ${\bf u}$ satisfies

[TABLE]

Set ${\bf e}={\bf u}-{\bf u}_{h}$. Then, subtract (4.1) from (3.3) to arrive at

[TABLE]

where {\bf\Lambda}({\mbox{\boldmath\phi}}_{h})=-b({\bf u},{\bf u},{\mbox{\boldmath\phi}}_{h})+b({\bf u}_{h},{\bf u}_{h},{\mbox{\boldmath\phi}}_{h}). Below, we derive the optimal error estimates of $||{\bf e}(t)||$ and $||\nabla{\bf e}(t)||$, for $t>0$.

In order to deal with the nonlinearity, an intermediate solution ${\bf v}_{h}$ is introduced which is a finite element Galerkin approximation to a linearized Kelvin-Voigt equation. The solution ${\bf v}_{h}$ satisfies

[TABLE]

with ${\bf v}_{h}(0)=P_{h}{\bf u}_{0}.$

Now, we split ${\bf{e}}$ as

[TABLE]

Here, $\xi$ is the error due to the approximation using a linearized Kelvin-Voigt equation (4.3), whereas $\eta$ denotes the error due to the non-linearity in the equation. A subtraction of (4.3) from (4.1) leads to the equation in $\xi$ as

[TABLE]

In order to derive optimal error estimates of $\xi$ in $L^{\infty}(\bf{L}^{2})$ and $L^{\infty}({\bf{H}}^{1})$-norms, we introduce the following auxiliary projection $V_{h}$ such that $V_{h}{\bf u}:[0,\infty)\rightarrow{\bf J}_{h}$ satisfying

[TABLE]

where $V_{h}{\bf u}(0)=P_{h}{\bf u}_{0}.$

With $V_{h}{\bf u}$ defined as above, we now split $\xi$ as

[TABLE]

To obtain estimates for ${{\bf e}}$, first of all, we establish a few estimates of $\zeta$ in Lemmas 4.1-4.7. Then with the help of $\zeta$ estimates, we derive various estimates of $\rho$ and \nabla\mbox{\boldmath\rho} in Lemmas 4.8 and 4.10. Finally, in Lemma 4.11, we derive estimates for $\eta$ and complete the proof of Theorem 3.1.

Lemma 4.1**.**

Assume that assumptions (A1)-(A2)* and (B1)-(B2) are satisfied. Then, there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimate holds true:*

[TABLE]

Proof. Multiply (4.5) by $e^{\alpha t}$ with \mbox{\boldmath\zeta}={\bf u}-V_{h}{\bf u}, use e^{\alpha t}\mbox{\boldmath\zeta}_{t}=\hat{}\mbox{\boldmath\zeta}_{t}-\alpha\hat{}\mbox{\boldmath\zeta} and substitute {\mbox{\boldmath\phi}}_{h}=P_{h}\hat{\mbox{\boldmath\zeta}}=\hat{}\mbox{\boldmath\zeta}+(P_{h}\hat{\bf u}-\hat{\bf u}) to arrive at

[TABLE]

Integrate (4) with respect to time from [math] to $t$ and apply (3.4) along with Young’s inequality. A simplification of resulting equation with a use of $\|{\bf u}_{0}-P_{h}{\bf u}_{0}\|=\|{\bf u}(0)-P_{h}{\bf u}(0)\|$ yields

[TABLE]

After applying Cauchy-Schwarz’s inequality in the first term of right-hand side, use Young’s inequality with $p=1/2$, $q=1/2$ to obtain

[TABLE]

A use of (4.8) in (4 ) leads to

[TABLE]

The first term on both sides will cancel out. To deal with the second term on right-hand side, rewrite it as

[TABLE]

Apply (4) in (4) along with Young’s inequality, (3.5) and (B1) to arrive at

[TABLE]

A use of a priori bounds for ${\bf u}$ and $p$ stated in Lemma 2.1 completes the proof. $\Box$

Remark 4.1**.**

Note that using (4), we rewrite the second term on the right-hand side of (4) and thereby write the entire right-hand side of (4) as an integration. This plays an important role in achieving weight function ${\bf\sigma}$ in the desired estimates. The presence of ${\bf\sigma}$ in the estimates is crucial in order to deal with the regularity issues at $t=0$ while making error estimates independent of $\kappa$.

Next, we prove the estimates for the time derivative of $\zeta$.

Lemma 4.2**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimates hold true:*

[TABLE]

Proof. Recall (4.5) now with {\mbox{\boldmath\phi}}_{h}=P_{h}{\mbox{\boldmath\zeta}}_{t}=\mbox{\boldmath\zeta}_{t}+(P_{h}{\bf u}_{t}-{\bf u}_{t}) to find that

[TABLE]

Rewrite the first term on the right-hand side of (4.13) as

[TABLE]

and substitute in (4.13) to obtain

[TABLE]

After using the Cauchy-Schwarz inequality and discrete incompressibility condition in (4), multiply the resulting equation by ${\bf\sigma}(t)$ to arrive at

[TABLE]

Integrate (4) with respect to time from [math] to $t$ and apply Young’s inequality, (3.5), $\rm({\bf B1})$ to obtain

[TABLE]

Now, the desired results in (4.11) follows by using Lemmas 2.1 and 4.1 in (4).

Next to prove (4.12), substitute \mbox{\boldmath\phi}_{h}=P_{h}\mbox{\boldmath\zeta}_{t} with \mbox{\boldmath\zeta}={\bf u}-V_{h}{\bf u} in (4.5) and arrive at

[TABLE]

A use of Cauchy-Schwarz’s inequality, (3.5) and discrete incompressibility condition in (4.18) yield

[TABLE]

After taking a square of (4.19) on both sides, multiply the resulting equation by $e^{2\alpha t}$. Then, integrate with respect to time from [math] to $t$ and use bounds from Lemmas 2.1, 4.1 to arrive at the desired result. This completes the rest of the proof. $\Box$

Lemma 4.3**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimate holds true:*

[TABLE]

Here, $\tau(t):=\min\{t,\,1\}$ and ${\bf\sigma}_{1}(t):=\tau^{2}(t)\,e^{2\alpha t}$.

Proof. Differentiate (4.5) with respect to time and substitute \mbox{\boldmath\phi}_{h}=P_{h}\mbox{\boldmath\zeta}_{t} to observe

[TABLE]

Note that, a use of

[TABLE]

(B1), (3.5), Cauchy-Schwarz’s inequality and Young’s inequality in (4.20) lead to

[TABLE]

Multiplication of (4.22) by ${\bf\sigma}_{1}$ and integration of the resulting equation from [math] to $t$ yield

[TABLE]

A use of estimates from Lemmas 2.2, 4.2 and a multiplication of resulting equation by ${\bf\sigma}^{-1}(t)$ complete the proof of Lemma 4.3. $\Box$

Below, in Lemma 4.4 we discuss the $L^{2}(\bf{L}^{2})$-estimate of $\zeta$. The similar kind of estimate has already been discussed in Lemma 5.3 of [26]. The difference between the estimate of $\zeta$ in Lemma 5.3 of [26] and Lemma 4.4 in this article is the presence of weight function ${\bf\sigma}$ in Lemma 4.4 which will be very helpful in making the error estimates independent of $\kappa$. Therefore, in order to justify the presence of ${\bf\sigma}$, we present a short proof highlighting only the modifications.

Lemma 4.4**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimate holds true for $t>0$:*

[TABLE]

Proof. For obtaining the desired estimates of $\zeta$, we appeal to the Aubin-Nitsche duality argument by assuming $({\bf w},q)$ to be the unique solution of the steady state Stokes system:

[TABLE]

satisfying the following regularity result:

[TABLE]

Form $L^{2}$-inner product between (4.23) and \hat{\mbox{\boldmath\zeta}} and use discrete incompressibility condition. Then, apply (4.5) with {\mbox{\boldmath\phi}}_{h} replaced by $P_{h}{\bf w}$ to obtain

[TABLE]

Once again, form an $L^{2}$-inner product between (4.23) and e^{\alpha t}{\mbox{\boldmath\zeta}}_{t} and use it to replace the last term in (4) as follows

[TABLE]

Apply Cauchy-Schwarz’s inequality, assumption (B1) and regularity estimates (4.26) along with Young’s inequality in (4). Then, integrate the resulting equation with respect to time from [math] to $t$ to obtain

[TABLE]

Using \|\mbox{\boldmath\zeta}(0)\|^{2}=\|{\bf u}_{0}-P_{h}{\bf u}_{0}\|^{2}, we write

[TABLE]

A use of orthogonality property of $P_{h}$ and Cauchy-Schwarz’s inequality yield

[TABLE]

A simplification of (4) leads to

[TABLE]

An application of (4) and (4.32) in (4.29) yield

[TABLE]

By using (3.5), we arrive at

[TABLE]

Since 0\leq\alpha<\displaystyle{\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, $(\nu-\alpha\kappa)>0.$ Then, use estimates from Lemmas 2.1 and 4.1 to complete the rest of the proof. $\Box$

Remark 4.2**.**

Here again, using (4)-(4.32), we tackle the first term on the right-hand side of (4.29) and express the entire right-hand side of (4.29) as an integration. As mentioned earlier in Remark 4.1, this provides ${\bf\sigma}(t)$ in the estimate which is used to handle regularity issues of the solution at $t=0$ in the process of making error bounds independent of $\kappa$.

Now, Lemma 4.5 provides the estimate for the time derivative \mbox{\boldmath\zeta}_{t}. Here again, the estimate differs from the estimate of \mbox{\boldmath\zeta}_{t} in Lemma 5.3 of [26] in terms of involvement of weight function ${\bf\sigma}$ and additional power of $\kappa$. As stated earlier, the presence of ${\bf\sigma}(t)$ and additional power of $\kappa$ in the estimate play a crucial role in making error estimates independent of $\kappa$. The proof proceeds in an exactly similar manner as the proof of Lemma 4.4 with the right-hand side of (4.23) replaced by e^{\alpha t}\mbox{\boldmath\zeta}_{t}. But in order to justify the presence of ${\bf\sigma}$ in the estimate, we present a short proof.

Lemma 4.5**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following holds true:*

[TABLE]

Proof. For obtaining the desired estimate of \mbox{\boldmath\zeta}_{t}, we replace the right-hand side of (4.23) by e^{\alpha t}\mbox{\boldmath\zeta}_{t} and form an $\bf{L}^{2}$-inner product of resulting equation with e^{\alpha t}\mbox{\boldmath\zeta}_{t}. Then, use (4.5) with {\mbox{\boldmath\phi}}_{h}=e^{\alpha t}P_{h}{\bf w} in a similar way as in the $\bf{L}^{2}$-estimate of $\zeta$ to obtain

[TABLE]

Multiply (4) by $\kappa$ and use Cauchy-Schwarz’s inequality with (3.5), approximation property (B1), regularity result (4.26) with right-hand side e^{\alpha t}\mbox{\boldmath\zeta}_{t}. Then, after squaring both sides of the resulting equation, perform an integration with respect to time from [math] to $t$ to obtain

[TABLE]

An application of Lemmas 2.1, 4.1, 4.2, 4.4 would lead to the desired estimates. $\Box$

Lemma 4.6**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimate holds true for $t>0$:*

[TABLE]

Proof. Once again, we apply the Aubin-Nitsche duality argument. Let $({\bf w},q)$ be the unique solution of the following steady state Stokes system:

[TABLE]

Now, using assumption (A1), $({\bf w},q)$ satisfies the following regularity result:

[TABLE]

Taking an $L^{2}$-inner product between (4.37) and \mbox{\boldmath\zeta}_{t} and using the discrete incompressibility condition, we obtain

[TABLE]

Now, by using (4.5) with {\mbox{\boldmath\phi}}_{h} replaced by $P_{h}{\bf w}$ and (4.38), the last term in (4.41) can be rewritten as

[TABLE]

Use (4.37) to rewrite last term in (4.42) as

[TABLE]

Apply (4.42), (4.43) in (4.41) to obtain

[TABLE]

A simplification of (4) yields

[TABLE]

After multiplying (4) by $\kappa\,{\bf\sigma}$, rewrite the resulting equation as

[TABLE]

An integration of (4) with respect to time from [math] to $t$ along with a use of Cauchy-Schwarz’s inequality, Young’s inequality leads to

[TABLE]

Apply the bounds from Lemmas 2.1, 4.1, 4.2, 4.4, 4.5 and the regularity estimates (4.40) to arrive at the desired result.

Next, to derive $\L^{\infty}(\bf{L}^{2})$ estimates of $\zeta$, follow the similar steps as in Lemma 4.4 with \hat{}\mbox{\boldmath\zeta} replaced by $\zeta$ and arrive at (4). Then, use Cauchy-Schwarz’s inequality to obtain

[TABLE]

Apply regularity estimates (4.26) with right-hand side as $\zeta$ along with (4.19), Lemmas 2.1, 4.2, 4.6 to compelte the rest part of the proof. $\Box$

Lemma 4.7**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following holds true:*

[TABLE]

where ${\bf\sigma}_{1}(t):=\tau^{2}(t)e^{2\alpha t}$ with $\tau(t):=\min\{t,\,1\}$.

Proof. A use of Cauchy-Schwaz’s inequality, (B1) and (3.5) in (4) yield

[TABLE]

Multiply (4) by ${\bf\sigma}_{1}$, apply (4.21) and integrate the resulting equation from [math] to $t$ to arrive at

[TABLE]

An application of Young’s inequality along with regularity estimates (4.40) leads to

[TABLE]

A use of Lemmas 2.2, 4.3, 4.6 would lead us to the desired result. $\Box$

Since \mbox{\boldmath\xi}=\mbox{\boldmath\zeta}+\mbox{\boldmath\rho} and the estimates of $\zeta$ are already derived, it suffices to derive the estimates of $\rho$ to obtain estimates for $\xi$. Below, in Lemma 4.8, we state without proof estimates of $\rho$. We skip the proof as it follows the similar lines as in the proofs of Lemma 5.6 ([2]) and Lemma 4.1 in this article. We also present a couple of estimates of $\xi$ which can be easily derived using the estimates of $\zeta$ and $\rho$.

Lemma 4.8**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimates hold true:*

[TABLE]

$\Box$

Lemma 4.9**.**

Under the assumptions (A1)-(A2)* and (B1)-(B2), there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimate holds true:*

[TABLE]

Proof. To estimate $L^{2}$-error, we use the following duality argument: For fixed $t>0$ with $t\in(0,T)$, let ${\bf w}(\tau)\in{\bf J}_{1}$, $q(\tau)\in L^{2}/{\rm I\!R}$ be the unique solution of the backward Stokes problem

[TABLE]

The pair $({\bf w},q)$ satisfy the following regularity estimates

[TABLE]

Form an $L^{2}$ inner product between (4.52) and $\xi$ to arrive at

[TABLE]

A use of (4.4) with \mbox{\boldmath\phi}_{h} replaced by $P_{h}{\bf w}$ in (4) yields

[TABLE]

Note that,

[TABLE]

A simplification of (4), using (4.56) leads to

[TABLE]

An integration of (4) with respect to time from [math] to $t$ along with Cauchy-Schwarz’s inequality yields

[TABLE]

The first term in (4) vanishes due to ${\bf w}(t)=0$ and the second term disappears due to the orthogonality property of $P_{h}$. Now, a use of Young’s inequality along with the regularity estimates (4.53) leads to

[TABLE]

Apply estimates from Lemmas 2.1, 4.8 to arrive at the desired result. $\Box$

Lemma 4.10**.**

Let the assumptions (A1)-(A2)* and (B1)-(B2) be satisfied. Then, there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimates hold true:*

[TABLE]

where ${\bf\sigma}_{1}(t):=\tau^{2}(t)e^{2\alpha t}$ with $\tau(t):=\min\{t,\,1\}$.

Proof. Subtracting (4.5) from (4.4), we find that

[TABLE]

Multiply (4.62) by $\kappa\,e^{2\alpha t}$, substitute \mbox{\boldmath\phi}_{h}=\mbox{\boldmath\rho} and use Cauchy-Schwarz’s inequality, Young’s inequality in the resulting equation. Then, integrate the equation from [math] to $t$ to arrive at

[TABLE]

Note that, \mbox{\boldmath\rho}=\mbox{\boldmath\xi}-\mbox{\boldmath\zeta}. A use of the triangle inequality along with Lemmas 4.4 and 4.9 yields

[TABLE]

An application of the results from Lemma 4.5 and (4.64) in (4.63) and a multiplication of the resulting equation by $e^{-2\alpha t}$ complete the proof of (4.60).

Next to prove (4.61), substitute \mbox{\boldmath\phi}_{h}=\mbox{\boldmath\rho} in (4.62) and multiply the resulting equation by ${\bf\sigma}_{1}$ to arrive at

[TABLE]

After applying Cauchy-Schwarz’s inequality and Young’s inequality, integrate the resulting equation with respect to time from [math] to $t$ to obtain

[TABLE]

Use estimates from (4.60), (4.64) and Lemma 4.7 to arrive at (4.61) and this completes the proof of Lemma 4.10. $\Box$

Now, we derive the proof of the main Theorem 3.1.

Note that {\bf e}={\bf u}-{\bf u}_{h}=({\bf u}-{\bf v_{h}})+({\bf v_{h}}-{\bf u}_{h})=\mbox{\boldmath\xi}+\mbox{\boldmath\eta}. A use of the triangle inequality, the inverse inequality and Lemmas 4.6, 4.10 lead to the following estimates of $\xi$.

[TABLE]

In Lemma 4.11, we present the estimates of $\eta$. For a proof, one may refer to [26] (Theorem 5.1, pp. 249 - 250).

Lemma 4.11**.**

Let the assumptions (A1)-(A2)* and (B1)-(B2) be satisfied. Then, there exists a positive constant $C=C(\kappa,\nu,\alpha,\lambda_{1},M)$ such that for \displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4\big{(}1+\kappa\lambda_{1}\big{)}}}, the following estimates hold true:*

[TABLE]

Moreover, under the assumptions of Theorem 3.1 and the uniqueness condition (2.30), the constant $K(t)=C$. That is, the estimates are valid uniformly with respect to time.

Proof of Theorem 3.1. The proof follows by using the triangle inequality, inverse inequality, (4.66) and Lemma 4.11. $\Box$

Following the similar steps as in [26] (Theorem 6.1) and using $\kappa$ independent estimates derived earlier, we arrive at the desired pressure error estimates in Theorem 3.2 and this completes the proof. $\Box$

5 Fully Discrete Approximation

In this section, we apply a backward Euler method for time discretization of the finite element Galerkin approximation (3) of (1.1)-(1.3). Let $\{t_{n}\}_{n=0}^{N}$ be a uniform partition of $[0,T]$, and $t_{n}=nk$, with time step $k>0$. For smooth function $\phi$ defined on $[0,T]$, set \mbox{\boldmath\phi}^{n}=\mbox{\boldmath\phi}(t_{n}) and \bar{\partial}_{t}\mbox{\boldmath\phi}^{n}=\frac{(\mbox{\boldmath\phi}^{n}-\mbox{\boldmath\phi}^{n-1})}{k}.

The backward Euler method applied to (3) determines a sequence of functions ${\{{\bf U}^{n}\}}_{n\geq 1}\in{\bf{H}}_{h}$ and ${\{P^{n}\}}_{n\geq 1}\in L_{h}$ as solutions of the following recursive nonlinear algebraic equations:

[TABLE]

Equivalently, we seek ${\{{\bf U}^{n}\}}_{n\geq 1}\in{\bf J}_{h}$ such that

[TABLE]

Next, in Lemma 5.1, we state a priori bounds for the discrete solution $\{{\bf U}^{n}\}_{n\geq 1}$. We skip the proof as it will be an imitation of the proof of Lemma 4.1 in [27].

Lemma 5.1**.**

With $\displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4(1+\lambda_{1}\kappa)}}$, choose $k_{0}$ so that for $0<k\leq k_{0}$

[TABLE]

Then the discrete solution ${\bf U}^{N}$, $N\geq 1$ of (5.2) satisfies

[TABLE]

$\Box$**

Next, we proceed to derive fully discrete estimates for the velocity error ${\bf e}^{n}={\bf U}^{n}-{\bf u}_{h}(t_{n})={\bf U}^{n}-{\bf u}_{h}^{n}$ and for the pressure error $\rho^{n}=P^{n}-p_{h}(t_{n})=P^{n}-p_{h}^{n}$. Below, in Lemma 5.2, we present the various estimates of ${\bf e}^{n}$. The proof of (5.4) follows the similar lines as in the proof of Theorem 5.1 of [27]. Therefore, we skip the proof. The estimates of $\|\bar{\partial}_{t}{\bf e}^{n}\|$ and $\kappa\|\bar{\partial}_{t}\nabla{\bf e}^{n}\|$ are also discussed in Lemma 5.1 of [27], but, these estimates involve $\kappa^{-1}$ term. Therefore, here we provide a short proof of (5.5) by only highlighting the steps involved in making estimates independent of the inverse power of $\kappa$.

Lemma 5.2**.**

Let $\displaystyle{0\leq\alpha<\frac{\nu\lambda_{1}}{4(1+\kappa\lambda_{1})}}$ and $k_{0}>0$ be such that for $0<k\leq k_{0}$, (5.3) is satisfied. For some fixed $h>0$, let ${\bf u}_{h}(t)$ satisfies (3.3). Then, there is a positive constant $C_{T}$ that depends on $T$ such that

[TABLE]

Proof. To prove (5.5), consider (3.3) at $t=t_{i}$ and subtract it from (5.2) to obtain

[TABLE]

where ${\bf\sigma}_{1}^{n}={\bf u}_{ht}^{n}-\bar{\partial}_{t}{\bf u}_{h}^{n}$ and \Lambda_{h}(\mbox{\boldmath\phi}_{h})=b({\bf u}_{h}^{n},{\bf u}_{h}^{n},\mbox{\boldmath\phi}_{h})-b({\bf U}^{n},{\bf U}^{n},\mbox{\boldmath\phi}_{h}).

Note that, applying Taylor’s series expansion in the interval $(t_{i-1},t_{i})$, Cauchy-Schwarz’s inequality, Young’s inequality and estimates from Lemma 3.1, we arrive at

[TABLE]

and

[TABLE]

Rewrite the nonlinear term and apply generalized Hölder’s inequality to observe that

[TABLE]

Now, substitute \mbox{\boldmath\phi}_{h}=\bar{\partial}_{t}{\bf e}^{n} in (5.6), drop the first term from left hand side and use (5.7)-(5) to observe that

[TABLE]

A use of (5.4), Lemmas 3.1, 5.1 yield

[TABLE]

Now, following the steps involved in arriving at the equation (107) from (106) in the proof Lemma 5.1 of [27], we arrive at

[TABLE]

A combination of (5.10) and (5.11) completes the rest of the proof. $\Box$

Remark 5.1**.**

Note that in the proof of Theorem 5.1 of [27], the presence of $\kappa^{-1}$ in the first term of right hand side of equation (93) is a typo, as the first term is a combination of equation (88) and (89), in which estimates are independent of $\kappa^{-1}$.

To prove the pressure error estimates, subtract (5) from (3) and write \mbox{\boldmath\rho}^{n}=P^{n}-p_{h}^{n} to obtain

[TABLE]

A use of Cauchy-Schwarz’s inequality along with (5.7)-(5), Lemmas 3.1, 5.1, 5.2 yields

[TABLE]

A combination of (5.12), Lemma 5.2 and Theorems 3.1, 3.2 lead to the following fully discrete error estimates.

Theorem 5.1**.**

Under the assumptions of Theorem 3.1 and Lemma 5.2, the following hold true:

[TABLE]

6 Numerical Experiments

This section conducts numerical experiments to validate our theoretical results obtained in Theorem 5.1 for finite element Galerkin approximations of (1.1)-(1.3). We apply mixed finite element $P_{2}$-$P_{0}$ for space discretization and backward Euler method for time discretization.

Example 6.1**.**

In this example, we choose right-hand side function ${\bf{f}}$ in such a way that the exact solution $({\bf u},p)=((u_{1},u_{2}),p)$ takes the following form:

[TABLE]

with $(x,y)\in(0,1)\times(0,1)$ along with the Dirichlet boundary condition. Here, the fluid viscosity $\nu$=1, time interval $(0,1]$ with final time $T=1$.

Tables 1, 2 represent the convergence rates for velocity in $L^{\infty}(\bf{L}^{2})$, $L^{\infty}({\bf{H}}^{1})$-norms, respectively, and Table 3 depicts the convergence rates for pressure in $L^{\infty}(L^{2})$-norm for different values of $\kappa$. The numerical convergence rates presented in tables validate the theoretical findings obtained in Theorem 5.1. Moreover, it can be inferred that the numerical results still hold true as $\kappa\rightarrow 0$.

In Tables 4, 5, we present the velocity error in $L^{\infty}(\bf{L}^{2})$-norm and the pressure error in $L^{\infty}(L^{2})$-norm, respectively, for different values of $\kappa$ and fixed $\nu=0.01$. It can be observed from the tables that the velocity and pressure errors are quite high and are not stable for the mesh size $h=1/2,\,1/4$ for the Navier-Stokes system with $\kappa=0$ and $\nu=0.01$. Therefore, more mesh refinement is needed to achieve the desired accuracy. To overcome this issue, we introduce a reasonably small presence of $\kappa$ to the Navier-Stokes system and make the system more regularized. Therefore, in this case by introducing a significantly small value of $\kappa$, $\kappa\in\{10^{-2},\,1\}$ to the Navier-Stokes system, the errors of the desired accuracy are achieved at a coarser mesh $h\in\{1/2,\,1/4\}$ with much less computational efforts. Note that, a significantly small value of $\kappa$ means that here the presence of $\kappa=10^{-4}$ in the Navier-Stokes system does not provide the desired accuracy as the errors are still quite high and are not stable for the mesh size $h=1/2,\,1/4$. The results in tables 4, 5 validate the fact that the Kelvin-Voigt model can be thought of as a $\kappa$ regularization of the Navier-Stokes model.

Example 6.2**.**

In this example, we take right-hand side function ${\bf{f}}=0$, initial condition ${\bf u}_{0}=(10x^{2}(x-1)^{2}y(y-1)(2y-1),-10y^{2}(y-1)^{2}x(x-1)(2x-1),y)$, $\nu=1$ and $\kappa=1$. Figure 1 represents velocity plots of the Kelvin-Voigt model and the Navier-Stokes model for final time $T=4$. We observe that the Kelvin-Voigt fluid velocity tends to zero at a slower rate in comparison to the Navier-Stokes fluid velocity. This confirms the fact that after instantaneous removal of the forces, the velocity of the Kelvin-Voigt fluid does not vanish instantaneously as in the case of the Navier-Stokes fluid.

Example 6.3**.**

This example deals with the benchmark problem lid-driven cavity flow on a unit square with zero body force. Here, no-slip boundary conditions are considered everywhere except non zero velocity $(u_{1},u_{2})=(1,0)$ on the upper part of the boundary, that is, the lid of the square is moving with a velocity $(1,\,0)$. For numerical results, we have considered lines $(0.5,y)$ and $(x,0.5)$, final time $T=40$, $h=1/32$ and $\nu=1$. In figure 2, we have compared the Kelvin-Voigt velocity to the steady-state Navier-Stokes velocity for large time and different values of $\kappa=\{1,0.001,0.00001\}$. The plots depict that the Kelvin-Voigt solution converges to the steady state solution for large time and as $\kappa\rightarrow 0$.

7 Summary

The article discusses some new higher order regularity estimates for the weak solution which are valid for all time $t>0$ and as $\kappa\rightarrow 0$. Semidiscrete optimal error estimates are derived for the velocity in $L^{\infty}(\bf{L}^{2})$, $L^{\infty}({\bf{H}}_{0}^{1})$-norms and for the pressure in $L^{\infty}(\bf{L}^{2})$-norm. Further, under uniqueness condition, these estimates are shown uniformly in time. Note that, the constants appearing in a priori error bounds are made independent of inverse powers of $\kappa$ by introducing weight functions in powers of $t$. In fact, an introduction of these weight functions takes care of regularity issues at time $t=0$. Further, the backward Euler method is applied for the complete discretization of the model and fully discrete optimal error estimates are derived. Finally, the article is concluded by presenting some numerical results which validate our theoretical observations.

Acknowledgement: The authors thank Professor Amiya K. Pani for his valuable comments and suggestions regarding this.