SUPG Approximation for the Oseen Viscoelastic Fluid Flow with Stabilized Lowest-Equal Order Mixed Finite Element Method

Abstract

In this article, a stabilized mixed finite element (FE) method for the Oseen viscoelastic fluid flow (OVFF) obeying an Oldroyd-B type constitutive law is proposed and investigated by using the Streamline Upwind Petrov–Galerkin (SUPG) method. To find the approximate solution of velocity, pressure and stress tensor, we choose lowest-equal order FE triples $P1$-$P1$-$P1$, respectively. However, it is well known that these elements do not fulfill the $inf$-$sup$ condition. Due to the violation of the main stability condition for mixed FE method, the system becomes unstable. To overcome this difficulty, a standard stabilization term is added in finite element variational formulation. The technique is applied herein possesses attractive features, such as parameter-free, flexible in computation and does not require any higher-order derivatives. The stability analysis and optimal error estimates are obtained. Three benchmark numerical tests are carried out to assess the stability and accuracy of the stabilized lowest-equal order feature of the OVFF.

1. Introduction

The research interest in viscoelastic fluids has increased, due to the connections with industrial applications. Generally, the structure of the viscoelastic fluids are formed by natural complex high molecular weight molecules which may have many internal degrees of freedom [1]. This is indeed a big motivation to the numerical and mathematical analysis of the governing equations [2]. Over the last century, there has been a significant challenge to formulate suitable constitutive model equations to describe the large deformation of the viscoelastic fluid. These were successfully introduced by James G. Oldroyd [3] in 1950 to critically study the behavior of the dilute solution of a polymeric molecule. Moreover, over the past few years, there have been many constitutive equations which have been proposed to describe viscoelastic fluids, e.g., the Phan–Thien–Tanner model, the Johnson–Segalman model, the Maxwell model and so on [4,5,6].

In the FE literature, Baranger and Sandri [7] carried out pioneering work, where they studied an Oldroyd-B fluid in discontinuous approximation and proved existence, uniqueness and error analysis. Later, Sandri extended the idea to the continuous approximation of the stress field for the stationary case [8], and the time-dependent case of the same continuous interpolation was analysed by Baranger in [9].

In Newtonian fluid flow, the Oseen equations are abridged to linearize the system. This is because the Oseen fluid flow model is the reduced linearized form of the Newtonian fluid which is described by the Navier–Stokes equation [10,11]. Moreover, the viscoelastic fluid flow model equations are non-linear equations in many terms. Hence, by taking the assumption of creeping fluid flow, the inertial part $(\mathbf{u}·\nabla \mathbf{u})$ of the momentum equation can be ignored. In this sense, the non-linearity arise only in the constitutive equation [12,13] which may be introduce in a linear form by fixing $\mathbf{u}(x)$ with a known velocity field $\mathbf{b}(x)$. The resulting system of equations can explicitly described the parameter space for $\alpha$, $\lambda$ and ${\parallel \nabla \mathbf{b}\parallel }_{\infty }$, which ensure the well-posedness of the continuous problem and its numerical approximation [14].

In the FE framework, the main difficulty arises in the viscoelastic fluid flow due to its hyperbolic constitutive equation. It needs a stabilization (upwinding ) technique to approximate the FE solution. There are two main approaches that are mainly used to solve the Oldroyd-B model by using the OVFF technique; the discontinuous Galerkin (DG) approximation and the (SUPG) estimate to deal with the spurious oscillations of the hyperbolic constitutive equation; see [15,16,17,18,19,20]. Herein, we consider SUPG method to solve the OVFF as an upwinding technique [21]. To the best of our knowledge, V.J. Ervin et al. introduced the SUPG method in [22] to approximate the model equations with the standard FE or Hood–Taylor FE pair ($P2$ for velocity, $P1$ for pressure and $P1$ for continuous stress), where the existence and uniqueness of a solution to the problem was shown. Later, the same authors studied a defect correction method in [23], where they used the same standard FE ($P2$-$P1$) to approximate the velocity and pressure but for discontinuous stress. They have used the conforming mixed elements, which satisfies the $inf$-$sup$ condition.

Generally, in practice of employing mixed FE methods to solving the model equations, the $inf$-$sup$ condition are considered as mile stones to insure the stability and accuracy of the scheme. Thus, the FE spaces for pressure and velocity are consider stable if they satisfy the $inf$-$sup$ condition, since the $inf$-$sup$ condition is an essential condition which imposes a complete correlation between two FE spaces so that they both have the required properties when utilized for the model equations. On the other hand, some mixed FE pairs which do not satisfy the $inf$-$sup$ condition are also popular in FE literature, i.e., lowest-equal order FE pair.

It is well-known that because of the simple logic and less computational cost, the lowest-equal order FE pairs are the most attractive choice. Therefore, this method is achieving more interest in computational fluid dynamics. Recent studies have focused on stabilization of the lowest-equal order FE pair to solve the Stokes equation [24,25], Navier–Stokes [26], Stokes–Darcy equation [27], and the Oseen equations [28]. Numerical tests show that the violation of the $inf$-$sup$ or Ladyzhenskaya-Babuska-Brezzi $(LBB)$ condition bring about physical pressure oscillations. In order to avoid the instability of model equations, the stabilized FE methods are applied to the incompressible flow. However, different stabilization methods have been developed and analyzed to circumvent the difficulties related with the stability of mixed FE methods which are studied in [29] and also several other ways have been used in the literature to stabilize the lowest-equal order FE pairs, e.g., the penalty methods, the regular methods, the multiscale enrichment methods [30], the local Gauss integration methods [31,32], projection methods [33] and many others [34,35].

This work is motivated by the lack of theoretical analysis and numerical computation of the OVFF with the approximation of lowest-equal order FE triples. Herein, we proved and executed the posednes for the OVFF with the lowest equal order triples. Moreover, the considered model is obtained by imposing a given velocity field $\mathbf{b}(x)$ in the objective derivative instead of the unknown velocity $\mathbf{u}(x)$ in the Oldroyd-B constitutive law, since, the stabilization techniques are often used to overcome the difficulty associated with the stability of the mixed FE method. Especially, the lowest-equal order FE are easily implementable in a scientific computational sense as compared to the high order FE. In this paper, a standard pressure term (stabilized) is added with the incompressibility condition in order to cure the classical $inf$-$sup$ condition between the velocity and the pressure spaces (as in the Stokes equations), and the SUPG scheme is employed in order to treat the hyperbolic nature of viscoelastic constitutive equations. In this context, there is no specific analysis available yet to approximation the OVFF by using the lowest-equal order triples.

Above all, the new approach differs from the existing stabilization techniques in different aspects; most notably, it is free of nonstandard data structures, approximation of higher order derivatives and specification of the mesh dependent parameters; it is optimally accurate and it always lead to symmetric system. Consequently, this new stabilized method can be easily applied by classical FE code with little additional coding effort.

The rest part of the paper is organized as follows. In Section 2, we introduce the viscoelastic fluid flow governing equations and the variational formulation. In Section 3, the new stabilized FE approximations are presented. In Section 4, the stability and error estimates for the stabilized FE solutions are achieved, respectively. In Section 5, we presented three different numerical experiments to verify the stability of the theocratical analysis. Finally, we summarize our work with a short conclusion in Section 6.

2. Model Equations

The steady-state (Johnson–Segalman) model equations can be seen in different viscoelastic articles; for the extensive theoretical modeling and numerical analysis for the considered model, we refer to [9,10,19,23] and the references therein. We considered the following equations:

$$\mathbf{\tau }+\lambda ((\mathbf{u}·\nabla \mathbf{\tau })+g_a(\mathbf{\tau },\nabla \mathbf{u}))-2\alpha D(\mathbf{u})=0.$$

(1)

where $\mathit{\tau }$ denotes the polymeric stress tensor, $\mathbf{u}$ the velocity vector and $\lambda$ denotes the Weissenberg number, which can be defined as the product of the relaxation time and a characteristic strain rate. While $0<\alpha <1$ is represented as the fraction of viscoelastic viscosity [36]. The rate of the strain tensor can be written as

$$D(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+{(\nabla \mathbf{u})}^T),$$

and the $g_a(\mathit{\tau },\nabla \mathbf{u})$ is given by

$$g_a(\mathbf{\tau },\nabla \mathbf{u})=\frac{1-a}{2}(\mathbf{\tau }(\nabla \mathbf{u})+{(\nabla \mathbf{u})}^T\mathbf{\tau })-\frac{1+a}{2}((\nabla \mathbf{u})\mathbf{\tau }+\mathbf{\tau }{(\nabla \mathbf{u})}^T).$$

(2)

Note: for the assumption of $a\in [-1,1]$ is defined to the material parameter, specifically for the choice of $a=1$ the Oldroyd-B constitutive model is always reducing from the Johnson-Segalman model.

We can write the momentum equation as

$$(\mathbf{u}·\nabla )\mathbf{u}-\nabla ·{\mathbf{\tau }}_{tot}=\mathbf{f},$$

where $\mathbf{f}$ denotes the body force. Here, the total stress tensor is given by ${\mathit{\tau }}_{tot}=-pI+{\mathit{\tau }}_N+\mathit{\tau }$ and it is the combination of the Newtonian part and the viscoelastic part. The Newtonian part can be further defined as ${\mathit{\tau }}_N=2(1-\alpha )D(\mathbf{u})$.

Since the viscoelastic fluid flows are very important to solve many practical (engineering) problems in non-Newtonian fluid mechanics, especially those being part to flow instabilities [9], we substituted the value of the total stress tensor in momentum equation it yields;

$$(\mathbf{u}·\nabla )\mathbf{u}-2(1-\alpha )\nabla ·D(\mathbf{u})-\nabla ·\mathbf{\tau }+\nabla p=\mathbf{f}.$$

(3)

In [5], Guillopé and Saut proved the following model equations with some assumptions for the “slow flow” i.e., in their contribution the $(\mathbf{u}·\nabla )\mathbf{u}$ term in (3) is ignored due to creeping flow. We also motivated with the same assumption and recast the stationary Oldroyd-B viscoelastic fluid flow for a domain $\mathsf{\Omega }$ in ${\mathbb{R}}^2$ with the Lipschitz continuous boundary $\mathsf{\Gamma }$.

$$\begin{array}{c}\hfill \mathbf{\tau }+\lambda \left((\mathbf{u}·\nabla )\mathbf{\tau }+g_a(\mathbf{\tau },\nabla \mathbf{u})\right)-2\alpha D(\mathbf{u})=0\mathrm{in}\mathsf{\Omega },\end{array}$$

(4)

$$\begin{array}{c}\hfill \nabla p-2(1-\alpha )\nabla ·D(\mathbf{u})-\nabla ·\mathbf{\tau }=\mathbf{f}\mathrm{in}\mathsf{\Omega },\end{array}$$

(5)

$$\begin{array}{c}\hfill \nabla ·\mathbf{u}=0\mathrm{in}\mathsf{\Omega },\end{array}$$

(6)

$$\begin{array}{c}\hfill \mathbf{u}=0\mathrm{on}\mathsf{\Gamma }.\end{array}$$

(7)

The dimensionless steady-state model equations under the open, bounded and connected domain $\mathsf{\Omega }$ are considered, with the homogenous Dirichlet boundary condition for the velocity $\mathbf{u}$. Since, in this case there is no inflow boundary, and thus, no boundary condition is required for stress [36]. To this end, we can re-write the OVFF model as

problem $(O)$ Find ($\mathit{\tau }$, $\mathbf{u},p)$ such that:

$$\begin{array}{cc}\hfill \mathbf{\tau }+\lambda (\mathbf{b}·\nabla )\mathbf{\tau }+\lambda g_a(\mathbf{\tau },\nabla \mathbf{b})-2\alpha D(\mathbf{u})& =0\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \nabla p-2(1-\alpha )\nabla ·D(\mathbf{u})-\nabla ·\mathbf{\tau }& =\mathbf{f}\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \nabla ·\mathbf{u}& =0\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \mathbf{u}& =\mathbf{0}\mathrm{on}\mathsf{\Gamma }.\hfill \end{array}$$

(11)

To make the system linear, we modify it with the known velocity $\mathbf{b}(x)$ in nonlinear terms of the constitutive equation instead of the unknown velocity $\mathbf{u}(x)$.

The Variational Formulation

In this section, we introduce some functional basic spaces for the mathematical analysis. The $L^2$ inner product is denoted by $(·,·)$ and while the special case for the $L^2(\mathsf{\Omega })$ and $L^{\infty }(\mathsf{\Omega })$ norms are assumed as $\parallel ·\parallel$ and ${\parallel ·\parallel }_{\infty }$, respectively. Moreover, the sobolev function space $W^{m,p}(\mathsf{\Omega })$ is represented as ${\parallel ·\parallel }_{W^{m,p}}$. The special notation for the sobolev case $W^{m,2}(\mathsf{\Omega })$ is being written as $H^m(\mathsf{\Omega })$. However, the norms and semi norms are represented in classical notations (see [37] ), i.e., ${\parallel ·\parallel }_m$ and ${|·|}_m$. We define the following functional spaces:

$$\begin{array}{ccc}\hfill \mathrm{Velocity}\mathrm{Space}:X:& =& H_0^1{(\mathsf{\Omega })}^2:=\{\mathbf{v}\in H^1{(\mathsf{\Omega })}^2:\mathbf{v}=\mathbf{0}\mathrm{on}\mathsf{\Gamma }\},\hfill \\ \hfill \mathrm{Pressure}\mathrm{Space}:Q:& =& L_0^2(\mathsf{\Omega })=\{q\in L^2(\mathsf{\Omega });{\int }_{\mathsf{\Omega }}qdx=0\},\hfill \\ \hfill \mathrm{Stress}\mathrm{Space}:S:& =& \{\mathbf{\tau }=({\mathbf{\tau }}_{ij});{\mathbf{\tau }}_{ij}={\mathbf{\tau }}_{ji}:{\mathbf{\tau }}_{ij}\in L^2(\mathsf{\Omega });i,j=1,2\}\hfill \\ & & \cap \{\mathbf{\tau }=({\mathbf{\tau }}_{ij});(\mathbf{b}·\nabla )\mathbf{\tau }\in L^2{(\mathsf{\Omega })}^{2\times 2}\}.\hfill \end{array}$$

It is of course well-known that the standard FEs X and Q (velocity and pressure spaces) satisfies the $inf$-$sup$ condition for the well-posedness of mixed formulations. However, there are many FE spaces satisfying $inf$-$sup$ condition, among the others the MINI (P1b,P1) elements, the Hood-Taylor (P2,P1) elements, Raviart–Thomas elements, Brezzi–Douglas–Marini finite element are common for the velocity $\mathbf{u}$ and pressure p, and discontinuous element $(P{1}_{dc})$ for stress tensor $\mathit{\tau }$ [14].

We make the following assumption for b(x), which is consistent with the existence result which has been set up for the viscoelasticity [38] for some $M>0$ as follows:

$$\mathbf{b}\in H_0^1(\mathsf{\Omega })\hspace{0.17em}\hspace{0.17em}\nabla ·\mathbf{b}=0\hspace{0.17em}\hspace{0.17em}{\parallel \mathbf{b}\parallel }_{\infty }\le M\hspace{0.17em}\hspace{0.17em}\parallel \nabla \mathbf{b}{\parallel }_{\infty }\le M<\infty .$$

The corresponding weak formulation of (8)–(11) is then given by

$$\begin{array}{cccc}\hfill (\mathbf{\tau },\mathbf{\sigma }+\lambda {\delta }_1(\mathbf{b}·\nabla \mathbf{\sigma }))+& (\lambda \mathbf{b}·\nabla \mathbf{\tau },\mathbf{\sigma }+\lambda {\delta }_1(\mathbf{b}·\nabla \mathbf{\sigma }))+(\lambda g_a(\mathbf{\tau },\nabla \mathbf{b}),\mathbf{\sigma }\hfill & & \\ & +\lambda {\delta }_1(\mathbf{b}·\nabla \mathbf{\sigma }))-(2\alpha D(\mathbf{u}),\mathbf{\sigma }+\lambda {\delta }_1(\mathbf{b}·\nabla \mathbf{\sigma }))\hfill & =& 0\forall \mathbf{\sigma }\in S,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill -(p,\nabla ·\mathbf{v})+2(1-\alpha )(D(\mathbf{u}),D(\mathbf{v}))+(\mathbf{\tau },D(\mathbf{v}))& =& (\mathbf{f},\mathbf{v})\forall \mathbf{v}\in X,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill (q,\nabla ·\mathbf{u})& =& 0\forall q\in Q.\hfill \end{array}$$

(14)

where ${\delta }_1$ is a positive constant. Multiplying Equations (13) and (14) by $2\alpha$ and adding the result to (12). Then for more simplicity, (12)–(14) can be reformulated to find $(\mathit{\tau },\mathbf{u},p)\in S\times X\times Q$ such that

$$\begin{array}{c}\hfill \mathcal{M}((\mathbf{\tau },\mathbf{u},p),(\mathbf{\sigma },\mathbf{v},q))=2\alpha (\mathbf{f},\mathbf{v}).\end{array}$$

(15)

Hence

$$\begin{array}{cc}\hfill \mathcal{M}((\mathbf{\tau },\mathbf{u},p),(\mathbf{\sigma },\mathbf{v},q))& =(\mathbf{\tau },\mathbf{\sigma }+\lambda {\delta }_1(\mathbf{b}·\nabla \mathbf{\sigma }))+(\lambda \mathbf{b}·\nabla \mathbf{\tau },\mathbf{\sigma }+\lambda {\delta }_1(\mathbf{b}·\nabla \mathbf{\sigma }))\hfill \\ & +(\lambda g_a(\mathbf{\tau },\nabla \mathbf{b}),\mathbf{\sigma }+\lambda {\delta }_1(\mathbf{b}·\nabla \mathbf{\sigma }))\hfill \\ & -(2\alpha D(\mathbf{u}),\mathbf{\sigma }+\lambda {\delta }_1(\mathbf{b}·\nabla )\mathbf{\sigma })+2\alpha (\mathbf{\tau },D(\mathbf{v}))\hfill \\ & +4\alpha (1-\alpha )(D(\mathbf{u}),D(\mathbf{v}))-2\alpha (p,\nabla ·\mathbf{v})+2\alpha (q,\nabla ·\mathbf{u})\hfill \end{array}$$

(16)

3. Approximation of the Problem with Finite Element Technique

In this section, we describe the FE framework for the proposed scheme and approximation properties. To this end, we introduce $\mathsf{\Omega }\subset {\mathbb{R}}^d$ (where dimension d = 2,3) as a domain and let $T^h$ is a triangulation of $\mathsf{\Omega }$, such as triangles K: $\overline{\mathsf{\Omega }}=\{\cup K:K\in T^h\}$. We consider $P_1(K)$ denote the space polynomial of degree on $K\in T^h$. Assume that some positive constant ${\upsilon }_1$,${\upsilon }_2$;

$${\upsilon }_{1}h\le h_K\le {\upsilon }_2{R}_K$$

where $h_K$ is the diameter of triangle K, $R_K$ is the diameter of the greatest ball included in K, and $h=ma{x}_{K\in T^h}{h}_K.$ We introduce the corresponding FE spaces to approximate $(\mathit{\tau },\mathbf{u},p)$ in a discrete way.

$$\begin{array}{ccc}\hfill X^h& :=& \{\mathbf{v}\in X\cap C^0{(\overline{\mathsf{\Omega }})}^d;{\mathbf{v}}_{|K}\in P_1{(K)}^d,\forall K\in T^h\},\hfill \\ \hfill Q^h& :=& \{q\in Q\cap C^0(\overline{\mathsf{\Omega }});q_{|K}\in P_1(K);\forall K\in T^h\},\hfill \\ \hfill S^h& :=& \{\mathbf{\sigma }\in S;{\mathbf{\sigma }}_{|K}\in P_1{(K)}^{d\times d};\forall K\in T^h\}.\hfill \end{array}$$

Note: It is well-known that the standard FEs $X^h$ and $Q^h$ does not fulfill the $inf$-$sup$ condition because of the lowest-equal order triples. To ensure stabilization, we introduce the local pressure projection $\Pi :L^2(\mathsf{\Omega })\to {\mathbb{R}}_0$. Where ${\mathbb{R}}_0$ is a piecewise constant space associated with the partition $T^h$. This bilinear, symmetric stabilization term is analysed, proposed and applied in [25,26]. Where, the stabilization term is stated as

$$G(p^h,q^h)=(p^h-\Pi p^h,q^h-\Pi q^h).$$

(17)

We mention the following lemmas which was given in article [24].

Lemma 1.

Let $Q^h$ and $X^h$ be the spaces defined above. There exist positive constants $C_1$ and $C_2$ satisfies:

$$\underset{v^h\in X^h}{sup}\frac{{\int }_{\mathsf{\Omega }}{p}^h\nabla ·v^{h}d\mathsf{\Omega }}{\parallel v^h{\parallel }_1}\ge C_1\parallel p^h{\parallel }_0-C_2\parallel \nabla p^h{\parallel }_0\forall p^h\in Q^h.$$

(18)

Lemma 2.

Let $Q^h$ and $X^h$ be the spaces then there exist a positive constant C satisfies:

$$Ch\parallel \nabla p^h{\parallel }_0\le \parallel p^h-\Pi p^h{\parallel }_0.$$

(19)

Stabilization Scheme

In order to filter the unstable factors, we can add the local stabilization term, which is based on pressure projection. Moreover, the stabilized methods in [24,25,26] are identical from a numerical point of view for the lowest-equal order approximations.

To seek $({\mathbf{u}}^h,{\mathit{\tau }}^h,p^h)$ $\in (X^h\times S^h\times Q^h)$ satisfies $\forall ({\mathbf{v}}^h,{\mathit{\sigma }}^h,q^h)$ $\in (X^h\times S^h\times Q^h)$ such that;

$$\begin{array}{cc}\hfill \tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))& =2\alpha (\mathbf{f},{\mathbf{v}}^h).\end{array}$$

(20)

where

$$\begin{array}{cc}\hfill \tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))& =({\mathbf{\tau }}^h,{\mathbf{\sigma }}^h+\lambda {\delta }_1(\mathbf{b}·\nabla ){\mathbf{\sigma }}^h)+(\lambda \mathbf{b}·\nabla {\mathbf{\tau }}^h,{\mathbf{\sigma }}^h+\lambda {\delta }_1(\mathbf{b}·\nabla ){\mathbf{\sigma }}^h)\hfill \\ & +(\lambda g_a({\mathbf{\tau }}^h,\nabla \mathbf{b}),{\mathbf{\sigma }}^h+\lambda {\delta }_1(\mathbf{b}·\nabla {\mathbf{\sigma }}^h))+2\alpha (\mathbf{\tau },D({\mathbf{v}}^h))\hfill \\ & -2\alpha (D({\mathbf{u}}^h),{\mathbf{\sigma }}^h+\lambda {\delta }_1(\mathbf{b}·\nabla ){\mathbf{\sigma }}^h)+4\alpha (1-\alpha )(D({\mathbf{u}}^h),D({\mathbf{v}}^h))\hfill \\ & -2\alpha (p^h,\nabla ·{\mathbf{v}}^h)+2\alpha (q^h,\nabla ·{\mathbf{u}}^h)+2\alpha G(p^h,q^h).\hfill \end{array}$$

(21)

For our convenience, we use some approximation properties [23], which we will use in our subsequent analysis. If there exists the interpolations (${\tilde{\mathit{\tau }}}^h,{\tilde{\mathbf{u}}}^h,{\tilde{p}}^h$) ∈$(S^h\times X^h\times Q^h)$ such that

$$\begin{array}{ccc}\hfill \parallel \mathbf{u}-{\tilde{\mathbf{u}}}^h{\parallel }_0+h\parallel \nabla (\mathbf{u}-{\tilde{\mathbf{u}}}^h){\parallel }_0& \le & C{h}^2\parallel \mathbf{u}{\parallel }_2,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \parallel \mathbf{\tau }-{\tilde{\mathbf{\tau }}}^h{\parallel }_0+h\mid \mathbf{\tau }-{\tilde{\mathbf{\tau }}}^h{\mid }_1& \le & C{h}^2\parallel \mathbf{\tau }{\parallel }_2,\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \parallel p-{\tilde{p}}^h{\parallel }_0& \le & Ch\parallel p{\parallel }_1.\hfill \end{array}$$

(24)

Throughout the paper we apply C as generic positive constant. In that case, we suppose this positive constant depends only on domain and always independent of the mesh size h, whose value may change from place to place. To show that (20) is a stable variational problem, we assume

$$\begin{array}{c}\hfill \parallel \Pi p{\parallel }_0\le C\parallel p{\parallel }_0,\parallel p-\Pi p{\parallel }_0\le Ch\parallel p{\parallel }_1.\end{array}$$

(25)

To end this section, we give some facts, which are easy to understand and also applicable in the analysis, i.e., let us consider the incompressibility condition $\nabla ·\mathbf{u}=0$ and $\mathbf{u}=0$ on $\mathsf{\Gamma }$, then it is not difficult to see that $2(D(\mathbf{u}),D(\mathbf{v}))=(\nabla \mathbf{u},\nabla \mathbf{v})$, and also we have $\parallel D(\mathbf{u})\parallel \le \parallel \nabla \mathbf{u}\parallel$.

4. Existence, Uniqueness and Error Bounds of the Problem

In this section, we prove the existence and uniqueness of the new proposed FE scheme for the approximation of the OVFF which satisfies the solution for all the positive $\lambda ,M,\alpha ,d$ and ${\delta }_1$ [10,21].

Theorem 1.

(Continuity) Given $\mathit{f}\in H^{-1}(\mathsf{\Omega })$, if $1-2\lambda Md>0,{\delta }_1>0$, there exist a unique solution $({\mathit{\tau }}^h,{\mathit{u}}^h,p^h)\in (S^h\times X^h\times Q^h)$ satisfying (20).

Proof. 

We know

$$\begin{array}{ccc}\hfill \tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,q^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))& =& F({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\forall ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\in (S^h,X^h,Q^h).\hfill \end{array}$$

(26)

The right hand side of (26 can be stated as $F(·):S^h\times X^h\times Q^h⟶\mathbb{R}$ is a functional defined by

$$\begin{array}{ccc}\hfill F({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)& =& 2\alpha (\mathbf{f},{\mathbf{v}}^h).\hfill \end{array}$$

Hence,

$$\begin{array}{ccc}\hfill \mid F({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid & \le & 2\alpha \parallel \mathbf{f}{\parallel }_{-1}\parallel {\mathbf{v}}^h{\parallel }_1\hfill \\ \hfill & \le & 2\alpha \parallel \mathbf{f}{\parallel }_{-1}\mid \parallel ({\mathbf{\tau }}^h,{\mathbf{v}}^h,q^h)\mid {\parallel }_{(S^h\times X^h\times Q^h)}.\hfill \end{array}$$

(27)

where

$$\begin{array}{ccc}\hfill \mid \parallel ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid {\parallel }_{(S^h\times X^h\times Q^h)}& =& {(\parallel {\mathbf{\sigma }}^h{\parallel }_0^2+\parallel \mathbf{b}·\nabla {\mathbf{\sigma }}^h{\parallel }_0^2+\parallel {\mathbf{v}}^h{\parallel }_1^2+\parallel q^h{\parallel }_0^2)}^{\frac{1}{2}}.\hfill \end{array}$$

(28)

We will show the bilinear form $\tilde{\mathcal{M}}(·,·,·)$ is continuous and coercive in $(S^h\times X^h\times Q^h)$ × $(S^h\times X^h\times Q^h)$ such that

$$\begin{array}{cc}\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,\hfill & {\mathbf{u}}^h,p^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))\hfill \\ & \le \parallel {\mathbf{\tau }}^h+\lambda \mathbf{b}·\nabla {\mathbf{\tau }}^h+\lambda g_a({\mathbf{\tau }}^h,\nabla \mathbf{b})-2\alpha D({\mathbf{u}}^h){\parallel }_0\parallel {\mathbf{\sigma }}^h+\lambda {\delta }_1\mathbf{b}·\nabla {\mathbf{\sigma }}^h{\parallel }_0\hfill \\ & +\parallel 4\alpha (1-\alpha )D({\mathbf{u}}^h)+2\alpha {\mathbf{\tau }}^h{\parallel }_0\parallel D({\mathbf{v}}^h){\parallel }_0+2\alpha d\parallel q^h{\parallel }_0\parallel \nabla {\mathbf{u}}^h{\parallel }_0\hfill \\ & +2\alpha d\parallel p^h{\parallel }_0\parallel \nabla {\mathbf{v}}^h{\parallel }_0+2\alpha \parallel p^h{\parallel }_0\parallel q^h{\parallel }_0,\hfill \\ & \le \left(\parallel {\mathbf{\tau }}^h{\parallel }_0+\lambda \parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0+2\lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0+2\alpha \parallel D({\mathbf{u}}^h){\parallel }_0\right)\parallel {\mathbf{\sigma }}^h+\lambda {\delta }_1\mathbf{b}·\nabla {\mathbf{\sigma }}^h{\parallel }_0\hfill \\ & +\left(4\alpha (1-\alpha )\parallel D({\mathbf{u}}^h){\parallel }_0+2\alpha \parallel {\mathbf{\tau }}^h{\parallel }_0\right)\parallel D({\mathbf{v}}^h){\parallel }_0+2\alpha d\parallel q^h{\parallel }_0\parallel \nabla {\mathbf{u}}^h{\parallel }_0\hfill \\ & +2\alpha d\parallel p^h{\parallel }_0\parallel \nabla {\mathbf{v}}^h{\parallel }_0+2\alpha \parallel p^h{\parallel }_0\parallel q^h{\parallel }_0,\hfill \\ & \le C_1\left((1+2\alpha )\parallel {\mathbf{\tau }}^h{\parallel }_0+\lambda \parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0+2\lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0+6\alpha \parallel D({\mathbf{u}}^h){\parallel }_0\right)\parallel ({\mathbf{v}}^h,{\mathbf{\sigma }}^h){\parallel }_{(S^h\times X^h)}\hfill \\ & +2\alpha d\parallel q^h{\parallel }_0\parallel \nabla {\mathbf{u}}^h{\parallel }_0+2\alpha d\parallel p^h{\parallel }_0\parallel \nabla {\mathbf{v}}^h{\parallel }_0+2\alpha \parallel p^h{\parallel }_0\parallel q^h{\parallel }_0,\hfill \\ & \le C\mid \parallel ({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h)\mid {\parallel }_{(S^h\times X^h\times Q^h)}\mid \parallel ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid {\parallel }_{(S^h\times X^h\times Q^h)}.\hfill \end{array}$$

(29)

Thus, it suffices to show the continuous property of the stabilized method. □

Theorem 2.

(coercivity) Let us consider a constant $1-2\lambda Md,\mathrm{Y}>0$ exist in such a way that the following inequality holds for all $({\mathit{\tau }}^h,{\mathit{u}}^h,p^h)\in (S^h\times X^h\times Q^h)$

$$\begin{array}{c}\hfill \underset{({\mathit{\sigma }}^h,{\mathit{v}}^h,q^h)\in S^h\times X^h\times Q^h}{sup}\frac{\tilde{\mathcal{M}}(({\mathit{\tau }}^h,{\mathit{u}}^h,p^h),({\mathit{\sigma }}^h,{\mathit{v}}^h,q^h))}{\mid \parallel ({\mathit{\sigma }}^h,{\mathit{v}}^h,q^h)\mid \parallel }\ge \mathrm{Y}\mid \parallel ({\mathit{\tau }}^h,{\mathit{u}}^h,p^h)\mid \parallel .\end{array}$$

(30)

Proof. 

To find weak coercivity of $\tilde{\mathcal{M}}$, we suppose a positive constant $\mathrm{Y}$ which is independent of h, for all $p^h\in Q^h\subset Q$, and interpolation ${\mathbf{w}}^h\in X^h$ of $\mathbf{w}$, there exists a positive constant ${\mathrm{Y}}_0$ and $\mathbf{w}\in X$ such that

$$\begin{array}{ccc}\hfill (div\mathbf{w},p^h)& =& \parallel p^h{\parallel }_0^2,\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill \parallel {\mathbf{w}}^h{\parallel }_1& \le & {\mathrm{Y}}_0\parallel p^h{\parallel }_0.\hfill \end{array}$$

(32)

Thanks to [24]

$$\begin{array}{ccc}\hfill {\int }_{\mathsf{\Omega }}{p}^h\nabla ·{\mathbf{w}}^{h}d\mathsf{\Omega }& \ge & C_1\parallel p^h{\parallel }_0^2-C_2\parallel (I-\Pi )p^h{\parallel }_0\parallel p^h{\parallel }_0.\hfill \end{array}$$

(33)

Using the fact about the velocity field $\mathbf{b}=0$ on boundary and $\nabla ·\mathbf{b}=0$, and integration by parts gives the following results [21],

$$({\mathit{\tau }}^h,\mathbf{b}·\nabla {\mathit{\tau }}^h)=-({\mathit{\tau }}^h,\mathbf{b}·\nabla {\mathit{\tau }}^h)\Rightarrow ({\mathit{\tau }}^h,\mathbf{b}·\nabla {\mathit{\tau }}^h)=0.$$

Now by setting the values $({\mathit{\sigma }}^h={\mathit{\tau }}^h,{\mathbf{v}}^h={\mathbf{u}}^h-\xi {\mathbf{w}}^h,q^h=p^h,)$ in (20), where $\xi$ is a real and positive parameter, we have

$$\begin{array}{cc}\hfill \tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\tau }}^h,& {\mathbf{u}}^h-\xi {\mathbf{w}}^h,p^h))\hfill \\ & =\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h))-\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),(0,\xi {\mathbf{w}}^h,0)).\hfill \end{array}$$

(34)

The right hand side of the Equation (34) can be bounded as

term (34)₁

$$\begin{array}{cc}\hfill \tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h))& =\left({\mathbf{\tau }}^h,{\mathbf{\tau }}^h+\lambda {\delta }_1(\mathbf{b}·\nabla ){\mathbf{\tau }}^h\right)+\lambda \left((\mathbf{b}·\nabla ){\mathbf{\tau }}^h,{\mathbf{\tau }}^h+{\delta }_1(\mathbf{b}·\nabla ){\mathbf{\tau }}^h\right)\hfill \\ & +\lambda \left(g_a({\mathbf{\tau }}^h,\nabla \mathbf{b}),{\mathbf{\tau }}^h+{\delta }_1(\mathbf{b}·\nabla ){\mathbf{\tau }}^h\right)+2\alpha (\mathbf{\tau },D({\mathbf{u}}^h))\hfill \\ & -2\alpha \left(D({\mathbf{u}}^h),{\mathbf{\tau }}^h+\lambda {\delta }_1(\mathbf{b}·\nabla ){\mathbf{\tau }}^h\right)+4\alpha (1-\alpha )\left(D({\mathbf{u}}^h),D({\mathbf{u}}^h)\right)\hfill \\ & -2\alpha \left(p^h,\nabla ·{\mathbf{u}}^h\right)+2\alpha \left(p^h,\nabla ·{\mathbf{u}}^h\right)+2\alpha \parallel (I-\Pi )p^h{\parallel }_0^2,\hfill \\ & =\parallel {\mathbf{\tau }}^h{\parallel }_0^2+{\delta }_1^2{\lambda }^2\parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0^2+\lambda \left(g_a({\mathbf{\tau }}^h,\nabla \mathbf{b}),{\mathbf{\tau }}^h+{\delta }_1(\mathbf{b}·\nabla ){\mathbf{\tau }}^h\right)\hfill \\ & -2\alpha \lambda (D({\mathbf{u}}^h),\mathbf{b}·\nabla {\mathbf{\tau }}^h)+4\alpha (1-\alpha )\parallel D({\mathbf{u}}^h){\parallel }_0^2+2\alpha \parallel (I-\Pi )p^h{\parallel }_0^2.\hfill \end{array}$$

A straight forward computation here is,

$$\begin{array}{ccc}\hfill \mid (\lambda {\mathbf{\tau }}^h\nabla \mathbf{b},{\mathbf{\tau }}^h)\mid & \le & \lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0^2.\hfill \end{array}$$

Similarly

$$\begin{array}{c}\hfill \lambda \mid (\nabla {\mathbf{b}}^T{\mathbf{\tau }}^h,{\mathbf{\tau }}^h)\mid ,\lambda \mid (\nabla \mathbf{b}{\mathbf{\tau }}^h,{\mathbf{\tau }}^h)\mid ,\lambda \mid ({\mathbf{\tau }}^h\nabla {\mathbf{b}}^T,{\mathbf{\tau }}^h)\mid \le \lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0^2.\end{array}$$

Hence

$$\begin{array}{ccc}\hfill \lambda \mid (g_a({\mathbf{\tau }}^h,\nabla \mathbf{b}),{\mathbf{\tau }}^h)\mid & \le & 2\lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0^2,\hfill \\ \hfill \lambda \mid (g_a({\mathbf{\tau }}^h,\nabla \mathbf{b}),{\delta }_1\mathbf{b}·\nabla {\mathbf{\tau }}^h)\mid & \le & 2\lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0{\delta }_1\parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0.\hfill \end{array}$$

(35)

By using Young’s inequality

$$\begin{array}{cc}\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h))\hfill & \ge \parallel {\mathbf{\tau }}^h{\parallel }_0^2+{\delta }_1^2{\lambda }^2\parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0^2+4\alpha (1-\alpha )\parallel D({\mathbf{u}}^h){\parallel }_0^2\hfill \\ & -2\lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0^2-2\lambda Md\parallel {\mathbf{\tau }}^h{\parallel }_0{\delta }_1\parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0\hfill \\ & -2\alpha \parallel D({\mathbf{u}}^h){\parallel }_0{\delta }_1\lambda \parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0+2\alpha \parallel (I-\Pi )p^h{\parallel }_0^2,\hfill \\ & \ge (1-2\lambda Md-{ϵ}_1{\delta }_1\lambda Md)\parallel {\mathbf{\tau }}^h{\parallel }_0^2+(1-\frac{Md}{{ϵ}_1}-{ϵ}_2){\delta }_1^2{\lambda }^2\parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0^2\hfill \\ & +(4\alpha (1-\alpha )-\frac{{\alpha }^2{\delta }_1}{{ϵ}_2})\parallel D({\mathbf{u}}^h){\parallel }_0^2+2\alpha \parallel (I-\Pi )p^h{\parallel }_0^2.\hfill \end{array}$$

(36)

term(34)₂

$$\begin{array}{cc}\hfill -\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),(0,\xi {\mathbf{w}}^h,0))& \ge -4\alpha (1-\alpha )\xi (D({\mathbf{u}}^h),D({\mathbf{w}}^h))-2\alpha \xi ({\mathbf{\tau }}^h,D({\mathbf{w}}^h))\hfill \\ & +2\alpha \xi (p^h,\nabla ·{\mathbf{w}}^h).\hfill \end{array}$$

(37)

By using (31)–(33) and Young’s inequality, the right side of Equation (37), can be estimate as

$$\begin{array}{cc}\hfill -4\alpha (1-\alpha )\xi (D({\mathbf{u}}^h),D({\mathbf{w}}^h))& \ge -4\alpha (1-\alpha )\xi \parallel D({\mathbf{u}}^h){\parallel }_0\parallel D({\mathbf{w}}^h){\parallel }_0\hfill \\ & \ge -4\alpha (1-\alpha )\xi \parallel D({\mathbf{u}}^h)\parallel {\mathrm{Y}}_0\parallel p^h{\parallel }_0\hfill \\ & \ge -4\alpha {(1-\alpha )}^2\xi \parallel D({\mathbf{u}}^h){\parallel }_0^2-\alpha \xi {\mathrm{Y}}_0^2\parallel p^h{\parallel }_0^2,\hfill \end{array}$$

$$\begin{array}{cc}\hfill -2\alpha \xi ({\mathbf{\tau }}^h,D({\mathbf{w}}^h))& \ge -2\alpha \xi \parallel {\mathbf{\tau }}^h{\parallel }_0{\mathrm{Y}}_o\parallel p^h{\parallel }_0\hfill \\ & \ge -2\alpha \xi {\mathrm{Y}}_o\parallel {\mathbf{\tau }}^h{\parallel }_0\parallel p^h{\parallel }_0\hfill \\ & \ge -\alpha \xi \parallel {\mathbf{\tau }}^h{\parallel }_0^2-\alpha \xi {\mathrm{Y}}_0^2\parallel p^h{\parallel }_0^2,\hfill \end{array}$$

$$\begin{array}{cc}\hfill 2\alpha \xi (p^h,\nabla ·{\mathbf{w}}^h)& \ge 2\alpha \xi (C_1\parallel p^h{\parallel }_0^2-C_2\parallel (I-\Pi )p^h{\parallel }_0\parallel p^h{\parallel }_0)\hfill \\ & \ge 2\alpha \xi C_1\parallel p^h{\parallel }_0^2-2\alpha \xi C_2\parallel (I-\Pi )p^h{\parallel }_0\parallel p^h{\parallel }_0\hfill \\ & \ge \alpha \xi (2C_1-C_2)\parallel p^h{\parallel }_0^2-\alpha \xi C_2\parallel (I-\Pi )p^h{\parallel }_0^2\hfill \end{array}$$

now Equation (37) gives,

$$\begin{array}{cc}\hfill -\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),(0,\xi {\mathbf{w}}^h,0))& \ge -\alpha \xi \parallel {\mathbf{\tau }}^h{\parallel }_0^2-4\alpha {(1-\alpha )}^2\xi \parallel D({\mathbf{u}}^h){\parallel }_0^2\hfill \\ & +\alpha \xi (2C_1-C_2-{\mathrm{Y}}_0^2)\parallel p^h{\parallel }_0^2-\alpha \xi C_2\parallel (I-\Pi )p^h{\parallel }_0^2.\hfill \end{array}$$

By summarize all the bounded values, we derive

$$\begin{array}{cc}\hfill \tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),& ({\mathbf{\tau }}^h,\mathbf{u}-\xi {\mathbf{w}}^h,p^h))\hfill \\ & \ge (1-2\lambda Md-{ϵ}_1{\delta }_1\lambda Md-\alpha \xi )\parallel {\mathbf{\tau }}^h{\parallel }_0^2\hfill \\ & +(4\alpha (1-\alpha )-\frac{{\alpha }^2{\delta }_1^2}{{ϵ}_2}-4\alpha {(1-\alpha )}^2\xi )\parallel D({\mathbf{u}}^h){\parallel }_0^2\hfill \\ & +(1-\frac{\lambda Md}{{ϵ}_1}-{ϵ}_2){\delta }_1{\lambda }^2\parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0^2+\alpha \xi (2C_1-C_2-{\mathrm{Y}}_0^2)\parallel p^h{\parallel }_0^2\hfill \\ & +\alpha (2-\xi C_2)\parallel (I-\Pi )p^h{\parallel }_0^2\hfill \\ & \ge C_3\parallel {\mathbf{\tau }}^h{\parallel }_0^2+C_4\parallel \mathbf{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0^2+C_5\parallel D({\mathbf{u}}^h){\parallel }_0^2+C_6\parallel p^h{\parallel }_0^2\hfill \\ & \ge \beta \mid \parallel ({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h)\mid {\parallel }_{S^h\times X^h\times Q^h}^2.\hfill \end{array}$$

(38)

By (32) and triangle inequality, we get

$$\begin{array}{cc}\hfill \mid \mid \mid ({\mathbf{\tau }}^h,{\mathbf{u}}^h-\xi {\mathbf{w}}^h,p^h)\mid \mid \mid & \le \parallel {\mathbf{\tau }}^h{\parallel }_0+\parallel \mathbf{u}-\xi {\mathbf{w}}^h{\parallel }_1+\parallel p^h{\parallel }_0\hfill \\ & \le \parallel {\mathbf{\tau }}^h{\parallel }_0+\parallel {\mathbf{u}}^h{\parallel }_1+\xi \parallel {\mathbf{w}}^h{\parallel }_0+\parallel p^h{\parallel }_0\hfill \\ & \le \parallel {\mathbf{\tau }}^h{\parallel }_0+\parallel {\mathbf{u}}^h{\parallel }_1+\xi {\mathrm{Y}}_0\parallel p^h{\parallel }_0+\parallel p^h{\parallel }_0\hfill \\ & \le C\mid \parallel ({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h)\mid \parallel .\hfill \end{array}$$

(39)

By combining (30), (38) and (39), we arrive at

$$\begin{array}{cc}\hfill \underset{({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\in S^h\times X^h\times Q^h}{sup}\frac{\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))}{\mid \parallel ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid \parallel }\ge & \frac{\tilde{\mathcal{M}}(({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h),({\mathbf{\tau }}^h,\mathbf{u}-\xi {\mathbf{w}}^h,p^h))}{\mid \parallel ({\mathbf{\tau }}^h,\mathbf{u}-\xi {\mathbf{w}}^h,p^h)\mid \parallel }\hfill \\ \hfill \ge & \frac{\beta \mid \parallel ({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h)\mid {\parallel }^2}{C\mid \parallel ({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h)\mid \parallel }\ge \frac{\beta }{C}\mid \parallel ({\mathbf{\tau }}^h,{\mathbf{u}}^h,p^h)\mid \parallel .\hfill \end{array}$$

By setting $\mathrm{Y}=\frac{\beta }{C}$, we can complete the proof of (30). □

Together, (29) and (30) imply that (20) is a stable variational problem

Theorem 3.

Let $(\mathit{\tau },\mathit{u},p)$ and $({\mathit{\tau }}^h,{\mathit{u}}^h,p^h)$ be the solution of (15) and (20), then the following error estimate holds:

$$\begin{array}{ccc}\hfill \parallel \mathbf{\tau }-{\mathbf{\tau }}^h{\parallel }_0+\parallel \mathit{b}·\nabla {\mathbf{\tau }}^h{\parallel }_0+\parallel \mathit{u}-{\mathit{u}}^h{\parallel }_0+\parallel p-p^h{\parallel }_0& \le & Ch.\hfill \end{array}$$

(40)

Proof. 

Subtracting (20) from (15) it yields

$$\begin{array}{ccc}\hfill \tilde{\mathcal{M}}((\mathbf{\tau }-{\mathbf{\tau }}^h,\mathbf{u}-{\mathbf{u}}^h,p-p^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))& =& 2\alpha G(p,q^h).\hfill \end{array}$$

(41)

By adding and subtracting $({\tilde{\mathit{\tau }}}^h,{\tilde{\mathbf{u}}}^h,{\tilde{p}}^h)$ and using the orthogonality gives

$$\begin{array}{cc}\hfill \tilde{\mathcal{M}}(({\tilde{\mathbf{\tau }}}^h-{\mathbf{\tau }}^h,{\tilde{\mathbf{u}}}^h-{\mathbf{u}}^h,{\tilde{p}}^h-p^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))=& \tilde{\mathcal{M}}(({\tilde{\mathbf{\tau }}}^h-\mathbf{\tau },{\tilde{\mathbf{u}}}^h-\mathbf{u},{\tilde{p}}^h-p),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))\hfill \\ & +2\alpha G(p,q^h).\hfill \end{array}$$

(42)

Using the weak coercivity bound (30), error orthogonality and (42)

$$\begin{array}{cc}\hfill \mathrm{Y}\mid \parallel & ({\tilde{\mathbf{\tau }}}^h-{\mathbf{\tau }}^h,{\tilde{\mathbf{u}}}^h-{\mathbf{u}}^h,{\tilde{p}}^h-p^h)\mid \parallel \hfill \\ & \le \underset{({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\in S^h\times X^h\times Q^h}{sup}\frac{\tilde{\mathcal{M}}(({\tilde{\mathbf{\tau }}}^h-{\mathbf{\tau }}^h,{\tilde{\mathbf{u}}}^h-{\mathbf{u}}^h,{\tilde{p}}^h-p^h),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))}{\mid \parallel ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid \parallel }\hfill \\ & =\underset{({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\in S^h\times X^h\times Q^h}{sup}\frac{\tilde{\mathcal{M}}(({\tilde{\mathbf{\tau }}}^h-\mathbf{\tau },{\tilde{\mathbf{u}}}^h-\mathbf{u},{\tilde{p}}^h-p),({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h))+2\alpha G(p,q^h)}{\mid \parallel ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid \parallel },\hfill \end{array}$$

from (25), we can bound

$$\begin{array}{c}\hfill 2\alpha G(p,q^h)\le C2\alpha G{(p,p)}^{1/2}\parallel q^h{\parallel }_0.\end{array}$$

(43)

From (30), we have

$$\begin{array}{cc}\hfill \tilde{\mathcal{M}}(({\tilde{\mathbf{\tau }}}^h-& \mathbf{\tau },{\tilde{\mathbf{u}}}^h-\mathbf{u},{\tilde{p}}^h-p),({\tilde{\mathbf{\tau }}}^h-{\mathbf{\tau }}^h,{\tilde{\mathbf{u}}}^h-{\mathbf{u}}^h,{\tilde{p}}^h-p^h))\hfill \\ & \le C(\parallel {\tilde{\mathbf{\tau }}}^h-\mathbf{\tau }{\parallel }_0+\parallel {\tilde{\mathbf{u}}}^h-\mathbf{u}{\parallel }_1+\parallel {\tilde{p}}^h-p{\parallel }_0+2\alpha \parallel (I-\Pi )p{\parallel }_0)\hfill \\ & \mid \mid \mid ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid \mid \mid \hfill \end{array}$$

(44)

As a result

$$\begin{array}{c}\mathrm{Y}\mid \parallel ({\tilde{\mathbf{\tau }}}^h-{\mathbf{\tau }}^h,{\tilde{\mathbf{u}}}^h-{\mathbf{u}}^h,{\tilde{p}}^h-p^h)\mid \parallel \hfill \\ \le \underset{({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\in S^h\times X^h\times Q^h}{sup}\times \hfill \\ \frac{C(\parallel {\tilde{\mathbf{\tau }}}^h-\mathbf{\tau }{\parallel }_0+\parallel {\tilde{\mathbf{u}}}^h-\mathbf{u}{\parallel }_1+\parallel {\tilde{p}}^h-p{\parallel }_0+2\alpha \parallel (I-\Pi )p{\parallel }_0)\mid \mid \mid ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid \mid \mid }{\mid \mid \mid ({\mathbf{\sigma }}^h,{\mathbf{v}}^h,q^h)\mid \mid \mid }\hfill \\ \le \frac{C}{\mathrm{Y}}(\parallel {\tilde{\mathbf{\tau }}}^h-\mathbf{\tau }{\parallel }_0+\parallel {\tilde{\mathbf{u}}}^h-\mathbf{u}{\parallel }_1+\parallel {\tilde{p}}^h-p{\parallel }_0+\parallel (I-\Pi )p{\parallel }_0)\hfill \end{array}$$

By the triangle inequality, we obtain Equation (40). □

5. Numerical Tests

In this section, we implement three different numerical tests for the support of our theoretical analysis. Numerical simulation for an analytic solution is studied in the first test which verifies the convergence order. In the second numerical test, we derive a viscoelastic cavity flow to show the pressure oscillation with the stabilization and without the stabilization term. In the third numerical test, to show the streamlines and contours of the pressure, we perform the well-known “4-to-1 abrupt contraction channel” fluid flow simulation. The numerical tests verify the accuracy and correctness of FE approximation for OVFF with lowest-equal order FE triples $P1$-$P1$-$P1$. In order to show the distinguished features of the new stabilized model, we compare numerical tests among three different FE schemes i.e., the standard $P1b$-$P1$-$P1$ FE (MINI elements), $P1$-$P1$-$P1$ without stabilization FE and $P1$-$P1$-$P1$ with stabilization FE. All numerical tests are performed by a public domain free software Freefem++ [39]. Furthermore, we have also drawn graphs and figures by MATLAB software and Tecplot 360 package.

5.1. Analytical Solution Test

In this example, the theoretical convergence rates are examined by applying fluid flow across a domain $\mathsf{\Omega }=[0,1]\times [0,1]$ with parameters $a=0$, $\lambda =5.0$ and $\alpha =0.5$, respectively. Different authors used this experimental pattern for Stokes and Navier–Stokes equations [23,36,40,41], where the function b(x) was chosen to be the exact solution of velocity $\mathbf{u}$. In this context, a right hand-side function is added to the (8) and $\mathbf{f}$ in (9) is studied with the help of following true solution.

$$\left\{\begin{array}{ccc}\hfill \mathbf{u}& =& \left(\begin{array}{c}-10(x^4-2x^3+x^2)(2y^3-3y^2+y)\\ 10(2x^3-3x^2+x)(y^4-2y^3+y^2)\end{array}\right),\hfill \\ \hfill p& =& -10.0(2x-1)(2y-1),\hfill \\ \hfill \mathbf{\tau }& =& 2\alpha D(\mathbf{u}).\hfill \end{array}\right.$$

For convenience we introduce some denotations

Scheme 1:

$P1b$-$P1$-$P1$ standard FE considered as-$p1bp1p1$ (MINI element spaces)

Scheme 2:

$P1$-$P1$-$P1$ with-out stabilization FE considered as-$p1p1p1unstb$

Scheme 3:

$P1$-$P1$-$P1$ with stabilization FE considered as-$p1p1p1stb$

The numerical results are presented in different tables;

Table 1. The error results with standard FE $P1b$-$P1$-$P1$.

h	$\left|\right|\mathit{u}-{\mathit{u}}^{\mathit{h}}{\left|\right|}_0$	$\left|\right|\mathit{u}-{\mathit{u}}^{\mathit{h}}{\left|\right|}_1$	$\left|\right|\mathit{p}-{\mathit{p}}^{\mathit{h}}{\left|\right|}_0$	$\left|\right|\mathit{\tau }-{\mathit{\tau }}^{\mathit{h}}{\left|\right|}_0$	CPU
$1/4$	0.017056200	0.22666200	0.22581000	0.1092200	0.078
$1/8$	0.005001480	0.10252500	0.06270120	0.0334345	0.312
$1/16$	0.001269810	0.04869610	0.01921590	0.0101251	1.235
$1/32$	0.000315101	0.02386040	0.00617860	0.0030970	5.156
$1/64$	0.000078142	0.01183820	0.00207000	0.0010032	19.454

,

Table 2. The error results before addition of stabilized term with FE $P1$-$P1$-$P1$.

h	$\left|\right|\mathit{u}-{\mathit{u}}^{\mathit{h}}{\left|\right|}_0$	$\left|\right|\mathit{u}-{\mathit{u}}^{\mathit{h}}{\left|\right|}_1$	$\left|\right|\mathit{p}-{\mathit{p}}^{\mathit{h}}{\left|\right|}_0$	$\left|\right|\mathit{\tau }-{\mathit{\tau }}^{\mathit{h}}{\left|\right|}_0$	CPU
$1/4$	0.02691260	0.235550	1.802590	0.2065030	0.062
$1/8$	0.00801367	0.122772	0.970023	0.0660082	0.25
$1/16$	0.00204691	0.060951	0.540790	0.0212788	0.938
$1/32$	0.00051251	0.030274	0.332033	0.0066949	3.813
$1/64$	0.00012819	0.015081	0.247405	0.0022083	17.422

and

Table 3. The error results after addition of a stabilized term with FE $P1$-$P1$-$P1$.

h	$\left|\right|\mathit{u}-{\mathit{u}}^{\mathit{h}}{\left|\right|}_0$	$\left|\right|\mathit{u}-{\mathit{u}}^{\mathit{h}}{\left|\right|}_1$	$\left|\right|\mathit{p}-{\mathit{p}}^{\mathit{h}}{\left|\right|}_0$	$\left|\right|\mathit{\tau }-{\mathit{\tau }}^{\mathit{h}}{\left|\right|}_0$	CPU
$1/4$	0.049930	0.452665	1.277330	0.310693	0.062
$1/8$	0.015867	0.182065	0.410506	0.119536	0.235
$1/16$	0.004361	0.071193	0.128865	0.041918	0.937
$1/32$	0.001129	0.029959	0.040235	0.014405	3.75
$1/64$	0.000286	0.013607	0.012769	0.004951	15.11

illustrates, the numerical values for the standard FE $p1bp1p1$, $p1p1p1unstb$ and $p1p1p1stb$, respectively. The $L^2$-norm for velocity $\mathbf{u}$, $H^1$-norm for velocity $\mathbf{u}$, $L^2$- norm for pressure p and $L^2$-norm for stress $\mathit{\tau }$, respectively are computed for different values of $h=1/8,1/16,1/32,1/64$. We compare the formulated results given in Table 3 with the standard result presented in Table 1. It is well-known that standard $p1bp1p1$ FE have better approximation which is also true for our current problem, since we compare the result of Table 2 and Table 3. The better numerical results of our scheme can be observed as expected in Table 3. Figure 1 and Figure 2 demonstrate the convergence orders for $p1bp1p1$, $p1p1p1unstb$ and $p1p1p1stb$ in $\left|\right|u-u^h{\left|\right|}_0$, $\left|\right|u-u^h{\left|\right|}_1$, $\left|\right|p-p^h{\left|\right|}_0$ and $\left|\right|\mathit{\tau }-{\mathit{\tau }}^h{\left|\right|}_0$, respectively. As to verify our desirable feature of the method, we compare the convergence orders where the experimental result shows that the order of convergence of velocity and stress are optimal order as h decreases for all three methods. However, the pressure for $p1p1p1unstb$ without stabilization can not obtain optimal error order.

5.2. The Viscoelastic Driven Cavity Flow

In this example, we apply the scheme in the given traditional benchmark problem for testing a numerical scheme known as “driven cavity flow”. Because, cavity flows have been used in many test cases for the testing of incompressible fluid dynamics algorithm for Stokes flow [31,42], for this example we set a domain $\mathsf{\Omega }=[0,1]\times [0,1]$ with the following boundary condition, i.e., on the top side $\{(x,1):0<x<1\}$ the velocity equal to $\mathbf{u}=(1,0)$, other boundaries are considered as zero Dirichlet condition. We consider the viscoelastic cavity problem for our numerical data comparison using the parameters $a=0$, $\lambda =5.0$, $\mathbf{f}=0$ and $\alpha =0.5$. The solution for the Oseen viscoelastic flow for the cavity flow is shown by the velocity field and pressure level lines. In Figure 3, we display cavity flow for three different schemes as: $p1bp1p1$ standard Galerkin FE approximation, $p1p1p1unstb$ without stabilized FE ( unstable) and $p1p1p1stb$ (new scheme) stabilized FE approximation. To sum up, we can see that the standard $p1bp1p1$ Galerkin method and $p1p1p1stb$ stabilized method completely resembles the physical rule for pressure. On the other hand, the result $p1p1p1unstb$ displays the pressure oscillations. As a result, the scheme is verified. Typically, through the finite element method, we have the following information as:

total triangles or elements = 8192,

nb boundary edges = 256,

number of nodes P1 = 4225.

5.3. 4-to-1 Contraction Channel Flow

We consider a third example to test the proposed scheme, which is a well-known problem for viscoelastic flow “4-to-1 contraction channel flow problem” which has a huge application in polymeric liquid industries and studied by many authors [20,43]. Moreover, 4-to-1 has been widely used in the literature to show the convergence, stability, and behavior of the streamlines of the contraction channel and the behavior of pressure [36,41]. The geometry of the 4-to-1 contraction commonly occurs in the forming of ‘die’ for the viscoelastic fluid. The domain is constructed in such a way that the channel lengths are sufficiently long for a fully developed Poiseuille flow at both the inflow and outflow boundaries. We presented a domain in Figure 4 for the shape of typical geometry for physical representation. Primarily, due to the sudden reduction in width, in the corner region, a vortex appears. In this example, we also study the behaviour of the fluid flow through the proposed scheme. The OVFF problem is based on the linear form of viscoelastic fluids. To apply the known function given in the theoretical part $\mathbf{b}(x)=(b1,b2)$ in numerical simulation, we perform following steps in the computational code.

• We first execute out put data of the approximate solution from the non-linear velocity for true solution $\mathbf{u}=(u_1,u_2)$.
• We use the executed solution of $(u_1)$ and $(u_2)$ as a known solution ($u_1=b_1$ and $u_2=b_2$). As a result, now the solution for the approximation considers for the linear one and it will be known for the velocity field $\mathbf{b}=(b_1,b_2)$, respectively.

As a matter of fact, the computations of the mesh are depicted in Figure 4 with $▵x_{min}=0.06250$ and $▵y_{min}=0.0156250$, with ${\mathsf{\Gamma }}_{in}=\{(x,y):x=0,0⩽y⩽1.0\}$ and ${\mathsf{\Gamma }}_{out}=\{(x,y):x=8.0,0⩽y⩽0.25\}$. The lower left corner of the domain corresponds to $x=y=0$. In this domain the boundary conditions for velocity are specifically designed as,

$$\begin{array}{c}{u}_1=\frac{1}{32}(1-y)(1+y),u_2=0\mathrm{on}{\mathsf{\Gamma }}_{\mathrm{in}},\hfill \\ u_1=2(\frac{1}{16}-y^2),u_2=0\mathrm{on}{\mathsf{\Gamma }}_{\mathrm{out}}.\hfill \end{array}$$

For stress $\mathit{\tau }$, on ${\mathsf{\Gamma }}_{in}$,

$$\begin{array}{c}{\mathbf{\tau }}_{11}=\frac{-\alpha \lambda (a+1){(-y/16)}^2}{(a+1)(a-1){\lambda }^2{(-y/16)}^2-1},\hfill \\ {\mathbf{\tau }}_{12}={\mathbf{\tau }}_{21}=\frac{-\alpha (-y/16)}{(a+1)(a-1){\lambda }^2{(-y/16)}^2-1},\hfill \\ {\mathbf{\tau }}_{22}=\frac{-\alpha \lambda (a-1){(-y/16)}^2}{(a+1)(a-1){\lambda }^2{(-y/16)}^2-1}.\hfill \end{array}$$

Here, the symmetry conditions are supposed on the bottom of the computational domain. In addition, the physical parameters $\alpha$, $\lambda$, and a are chosen to 1, $8/9$, $0.7$ and 1, respectively.

In Figure 5, we illustrate the streamlines of $p1bp1p1$ for standard FE. The left second figure depicts the streamlines of $p1p1p1unstb$ without stabilization and the third figure demonstrates the streamlines of $p1p1p1stb$ with stabilization for the physical parameter $\alpha$, $\lambda$, and a are chosen consequently as 1, $8/9$, $0.7$ and 1, respectively. It can be clearly seen that without the stabilization term for $p1p1p1unstb$, FE streamline behaviours are different and show some irregular shape than standard FE $p1bp1p1$ and with stabilization of $p1p1p1stb$ FE. In this perspective, the behavior of the streamlines with $p1bp1p1$ and $p1p1p1stb$ appear in a similar traditional $4:1$ abrupt contraction channel which resembles the physical rule completely. The behavior of the streamlines confirms the theoretical results for the stabilization by lowest-equal order FE in approximation of the OVFF.

In Figure 6, we have presented three figures; where pressure with $p1bp1p1$ for standard Galerkin FE is shown in the first figure (from left), $p1p1p1unstb$ FE without stabilization in the second figure, and $p1p1p1stb$ FE with stabilization in the third figure. A difference can be seen between the figure of $p1p1p1stb$ FE with stabilization being almost similar to the figure of $p1bp1p1$ for standard FE. On the other hand, the figure of $p1p1p1unstb$ FE without stabilization is completely different from the rest of the figures. Furthermore, the oscillation of pressure in Figure 6 part b, is clearly seen. These results also affirm that the theoretical analysis for the pressure stabilization under lowest-equal order FE triples for approximation of the OVFF.

6. Conclusions and Future Work

In this contribution, we studied an approximation of the Oseen viscoelastic fluid flow (OVFF) problem with lowest-equal order FE $P1$-$P1$-$P1$ method. In the standard Galerkin FE method, the $inf$-$sup$ (or $LBB$) conditions are essential to finding the well-posedness of the model equations, while the lowest-equal order FE triples violate the necessary condition $inf$-$sup$. This deficiency of the important condition makes the system unstable. To deal with the instability, we added symmetric and parameter-free extra pressure terms in the discrete variational form. We proved the well-posedness of the model equation with the lowest-equal order triples. Moreover, the optimal error estimates are derived, and three different numerical examples are performed which are consistent with the theoretical prediction. Most importantly, this method is simple and computationally efficient compared to other stabilized methods. In future studies, we can extend this method to approximate the solution with $P1$-$P0$-$P1$ triples.