Research Article | Open Access

# A Defect-Correction Method for Time-Dependent Viscoelastic Fluid Flow Based on SUPG Formulation

**Academic Editor:**M. De la Sen

#### Abstract

A defect-correction mixed finite element method for solving the time-dependent Johnson-Segalman viscoelastic equations in two dimensions is given. In the defect step, the constitutive equation is computed with the artificially reduced Weissenberg parameter for stability, and the resulting residual is corrected in the correction step on the same grid. A streamline upwind Petrov-Galerkin (SUPG) approximation is used to stabilize the hyperbolic character of the constitutive equation for the stress. We establish a priori error estimates for the defect step and the first correction step of the defect correction method. The derived theoretical results are supported by numerical tests.

#### 1. Introduction

Time-dependent calculations of viscoelastic flows are important to the understanding of many problems in non-Newtonian fluid mechanics, particularity those related to flow instabilities. Error analysis of finite element approximations to time-dependent viscoelastic flow was first analyzed by Baranger and Wardi in [1], using an implicit Euler temporal discretization and a discontinuous Galerkin (DG) approximation for the hyperbolic constitutive equation. In [2, 3], Ervin and Miles have analyzed the problem using an implicit Euler time discretization and an SUPG discretization for the constitutive equation. The key of SUPG method, first introduced by Hughes and Brooks [4], is to stabilize the numerical solution by adding dissipative effects limited to the direction in which their influence is needed. Error analysis of a modified Euler-SUPG approximation to time-dependent viscoelastic flow problem was presented by Bensaada and Esselaoui in [5]. In [6], Ervin and Heuer have proposed a Crank-Nicolson time discretization scheme for this problem. A fractional step -method for the time-dependent problem is described and analyzed in [7]. The reader may find more relevant work, for example, [8, 9] and the reference therein.

In addition, when the Weissenberg number increases, boundary layers for the stress develop. This will add to the difficulty of computing accurate numerical approximations. There are also many papers in developing stable numerical algorithms for high Weissenberg number flows (see [10–15] and the references therein) over the years. In [10–12], a defect-correction method has successfully applied in steady viscoelastic and Oseen-viscoelastic fluid flow model, aiming at high Weissenberg number. The basic idea of defect correction method for viscoelastic fluid flow is as follows. In the defect step, the Weissenberg number is reduced artificially by using a mesh-dependent parameter to obtain better convergence of an iteration scheme. Then, in a correction step, the initial approximation is improved using the residual correction. The defect correction method is an iterative improvement technique for increasing the accuracy of a numerical solution without applying a grid refinement. Due to its good efficiency, there are many works devoted to this method, for example, [16–22].

Inspired by [2, 3, 10–12, 21, 22], in this paper, we will use the defect correction method for solving the time-dependent Johnson-Segalman viscoelastic equations in two dimensions. Approximation in space is made by finite element (FE) method. Approximate stress, velocity, and pressure are, respectively, continuous, continuous, and continuous. As the constitutive equation has hyperbolic character for the stress, we use streamline upwind Petrov-Galerkin (SUPG) formulation to stabilize the constitutive equation. Approximate in time is made by Euler scheme. Finally, we establish a priori error estimates for the defect step and the first correction step of the defect correction method. While in [22], the authors present this problem by DG stabilized style, we know that these two styles have many differences. Choice of the Johnson and Segalman model is not essential, and the results obtained can be extended to more realistic differential models likes those of Phan-Thien and Tanner, Giesekus and Oldroyd, provided a Newtonian viscosity is present.

Detailed descriptions, analysis, and numerical examples of the defect correction method for the viscoelastic fluid flow are presented in the rest of the paper as follows. In the next section, we present the time-dependent viscoelastic fluid flow of the Johnson-Segalman model and its variational formulation. In Section 3, we introduce finite element spaces and present the fully discrete approximation scheme (Euler in time, FE in space). Defect correction algorithm is given in Section 4. Then a priori error estimates for the defect step and the first correction step of the defect correction method are derived in Section 5. Numerical results are presented in Section 6. Finally, a summary and discussion of continuing work are presented in Section 7.

#### 2. The Johnson-Segalman Model and Its Variational Formulation

In this section, we consider a fluid flowing in a connected, bounded polygonal domain with the Lipschitz continuous boundary .

The following notation will be used. is equipped with cartesian coordinates . For a function , is written and is written . Einstein's convention of summation is used. For a scalar function , gradient of is a vector with . For a vector , gradient of is a tensor with . For a vector function , divergence of is a scalar . with . For a tensor function and a vector function , divergence of is a vector with and denotes the operator .

Let be a given final time, we consider the time-dependent Johnson-Segalman viscoelastic equations as follows: where is the Weissenberg number, is the Reynolds number, may be considered as the fraction of viscoelastic viscosity ( for Maxwell's model and this case is excluded here, see reference [23] for a description of various models), and the body forces. The unknowns are the fluid velocity vector , the pressure , and the stress , which is the viscoelastic part of the total stress tensor . In (2.1), is the rate of the strain tensor, for all , is defined by Note that an Oldroyd B constitutive model is obtained when in .

The inner product is denoted by , and the norm by , with the special cases of and norms being written as and . For , we denote the norm associated with the Sobolev space by , with the special case being written as with the norm and seminorm . In order to introduce a variational formulation, we set the spaces , , , as follows:

The corresponding weak formulation of problem (2.1) is then obtained: find such that for all

Using the weak divergence free space , the weak formulation (2.4) can be written as follows: find such that for all

Existence results for problem (2.1) are proved in [24] for the "slow flow" model of (2.1) (i.e., the term in momentum equation is ignored). They are of two kinds: local existence in time of strong solutions in space variable in a domain, and global existence (in time) of a unique solution for and under a small data assumption on **, ****, **, . For a more complete discussion of existence and uniqueness issues, see [25].

#### 3. Finite Element Approximation

In this section, we present a fully discrete approximation to (2.5). We begin by describing the finite element approximation framework.

Suppose is a uniformly regular triangulation of such that and assume that there exist positive constants such that , where is the diameter of , is the diameter of the greatest ball included in , and . Throughout the paper, the constants denote different constants which are independent of mesh size and time step . We use the classical Taylor-Hood FE for the approximation in space of in velocity, in pressure; we consider approximation of the stress . The corresponding FE spaces are where denotes the space of polynomials of degree on . It is well known [26, 27] that the Taylor-Hood pair satisfies the discrete inf-sup (or LBB) condition.

To obtain the fully discrete approximation, the time derivatives are replaced by backward Euler differences, and the nonlinear terms are lagged. In order to stabilize the hyperbolic constitutive equation, an SUPG formulation is used to avoid spurious oscillations in the numerical approximation. Let be an integer, we divide the interval into intervals of equal length and denote . For notational convenience, we denote and The following norms are used in the analysis: and when is defined on the entire time interval , we denote

Then, based on SUPG formulation, the fully discrete approximating systems of (2.5) is the following;

given ; , , find such that for all where , is a small positive constant.

The goal of the parameter is used to suppress the production of spurious oscillations in the approximation. The discretization of the constitutive equation is the usual Galerkin finite element method if .

The existence of a unique solution and a priori error estimate to (3.5) can be found in [2, 3, 5].

#### 4. Defect Correction Method

In this section, our defect correction method used in computing the solution to (3.5) is described as follows.

*Algorithm 1 (Defect-correction method). **Step 1. *Solve the defected problem: find such that for all
where is chosen such that .*Step 2. *For , solve the correction problem: find such that for all
The initial value approximations are taken to be , and . In order to ensure computability of the algorithm, we begin by showing that (4.1)-(4.2) and (4.3)-(4.4) are uniquely solvable for at each time step . We use the following induction hypothesis, which simply states that the computed iterates , are bounded independent of and :

In next section, we will prove the induction hypothesis is right.

Lemma 4.1. *Suppose is true. For sufficiently small step size , then Step 1 of Algorithm 1 admits a unique solution .*

*Proof. *Letting , multiplying (4.1) by , and adding to (4.2), we get
where the bilinear form is defined by

As (4.1)-(4.2) represent a finite system of liner equations, the positivity of is a sufficient condition for the existence and uniqueness of .

We now estimate each term of . We have
where we have used the induction hypothesis and the Korn's lemma, is Korn constant.

Substituting all above equations into , we have
Choosing , and , , the bilinear form is positive definite. Hence, for chosen sufficiently small, Step 1 of Algorithm 1 admits a unique solution .

Similarity as Lemma 4.1, we can derive that Step 2 of Algorithm 1 admits a unique solution .

#### 5. A Priori Error Estimates

In this section, we explore the error estimates in approximating the true solution of (2.1) by the defect step approximation and one correction step approximation . The main results of this section are presented in the following theorem.

Theorem 5.1. *Suppose that (2.1) has a solution . In addition, assume that , , and , , , for all , and satisfy (4.1)-(4.2) and (4.3)-(4.4), respectively. Then, there exists a constant , such that
**
where
*

*Remark 5.2. *As mentioned in Algorithm 1, , thus, for , the essence of the estimates (5.1) and (5.2) is

Before deriving the estimates (5.2), (5.9), we first introduce some notation and some approximating properties of finite element spaces. Let represent the solution of (2.5) at time . Assume . Let denote the Stokes projection of into , and a Clement interpolant of [26, 28]. The following approximating properties are right.

Let and denote the solutions of (4.1)-(4.2) and (4.3)-(4.4), respectively. We denote , as follows:

In order to establish Theorem 5.1, we need another induction hypothesis, that is,

The induction hypothesis will be proved later.

*Proof. *We first give the profile of the proof. The proof of Theorem 5.1 is established in two steps.*Step 1. *Prove the error estimate (5.1) is right. We divide Step 1 into two substep.*Substep 1.1. *Under the induction hypothesis and for , we prove the error estimate (5.1) is true.*Substep 1.2. *We prove that the two induction hypotheses and for are right.*Step 2. *Show that the error estimate (5.2) is true.

Now, we start to prove Theorem 5.1.

Step 1. Prove the error estimate (5.1) is right.

Substep 1.1. Under the induction hypotheses and , we prove the error estimate (5.1) is true.

As being the exact solution of (2.1), it satisfies the following consistency equation: for all , in particular, at time (noting that we denote ),
Equation (5.8) can also be written as
Multiplying (4.1) by and adding to (4.2), for all , we get
Subtracting (5.9) from (5.10), for all , we obtain the following equation for and :
Substituting , and into (5.11), we get
where
Using the identity , we obtain
First multiply (5.14) by and sum it form to , we obtain
Then, we get, remarking that and ,

For controlling each term on the right-hand side of (5.16), the assumption and are needed here. We will prove the two induction hypotheses in the next subsection. Let us estimate each term of the of (5.16). For details, please see [2, 3, 6]. We start from the first four terms of the of (5.16);
Now, we estimate each term of . As has many terms and in order to make the proof more clear, we divide the terms of into four parts, that is, terms of , terms of , terms of , other terms of . Firstly, we estimate terms of ;
For the terms of , we have
For the terms of , we can obtain
Now, we estimate other terms of ;
Combining the estimates in (5.18)–(5.21), we obtain the estimate for :
Choosing , and substituting (5.17) and (5.22) into (5.16) yields