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

Authors: Yunzhang Zhang; Yanren Hou; Ganshan Yang
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## 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.


---

## Body

## 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, P1 continuous, P2 continuous, and P1 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 Section3, 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Ω⊂ℝd′,d′=2 with the Lipschitz continuous boundary Γ.The following notation will be used.R2 is equipped with cartesian coordinates xi,i=1,2. For a function u, ∂u/∂xi is written u,i and ∂u/∂t is written ut. Einstein's convention of summation is used. For a scalar function p, gradient of p is a vector ∇p with (∇p)i=p,i. For a vector u, gradient of u is a tensor ∇u with (∇u)ij=ui,j. For a vector function u, divergence of u is a scalar ∇u=ui,i. with (∇u)ij=ui,j. For a tensor function τ and a vector function u, divergence of τ is a vector ∇·τ with (∇·τ)i=τij,j and u·∇ denotes the operator ui(∂/∂xi).LetT>0 be a given final time, we consider the time-dependent Johnson-Segalman viscoelastic equations as follows:(2.1)Re(ut+u⋅∇u)-∇⋅σ-2(1-α)∇⋅D(u)+∇p=f,inΩ×(0,T),λσt+σ+λ(u⋅∇)σ+λga(σ,∇u)-2αD(u)=0,inΩ×(0,T),div⁡u=0,inΩ×(0,T),u=0,on⁡Γ×(0,T),u(0,x)=u0(x),σ(0,x)=σ0(x),on⁡Ω×{0},
where λ is the Weissenberg number, Re is the Reynolds number, 0<α<1 may be considered as the fraction of viscoelastic viscosity (α=1 for Maxwell's model and this case is excluded here, see reference [23] for a description of various models), and f the body forces. The unknowns are the fluid velocity vector u=(u1,u2), the pressure p, and the stress σ, which is the viscoelastic part of the total stress tensor σtot=σ+2(1-α)D(u)-pI. In (2.1), D(u)=(1/2)(∇u+∇uT) is the rate of the strain tensor, for all a∈[-1,1], ga(σ,∇u) is defined by(2.2)ga(σ,∇u)=1-a2(σ∇u+(∇u)Tσ)-1+a2((∇u)σ+σ(∇u)T).
Note that an Oldroyd B constitutive model is obtained when a=1 in ga(σ,∇u).TheL2(Ω) inner product is denoted by (·,·), and the Lp(Ω) norm by ∥·∥Lp, with the special cases of L2(Ω) and L∞(Ω) norms being written as ∥·∥ and ∥·∥∞. For k∈N, we denote the norm associated with the Sobolev space Wm,p(Ω) by ∥·∥Wm,p, with the special case Wm,2(Ω) being written as Hm(Ω) with the norm ∥·∥m and seminorm |·|m. In order to introduce a variational formulation, we set the spaces X, Q, S, V as follows:(2.3)X:=H01(Ω)2:={v∈H1(Ω)2:v=0on⁡Γ},Q:=L02(Ω)={q∈L2(Ω);∫Ωqdx=0},S:={τ=(τij);τij=τji;τij∈L2(Ω);i,j=1,2}∩{τ=(τij);v⋅∇τ∈L2(Ω),∀v∈X},V:={v∈X;∫Ωq(∇⋅v)dx=0,∀q∈Q}.The corresponding weak formulation of problem (2.1) is then obtained: find (σ,u,p)∈(S,X,Q) such that for all (τ,v,q)∈(S,X,Q)(2.4)(Reut+Reu⋅∇u,v)+(σ,D(v))+2(1-α)(D(u),D(v))+(p,∇⋅v)=(f,v),λ(σt,τ)+(σ,τ)+λ((u⋅∇)σ,τ)+λ(ga(σ,∇u),τ)-2α(D(u),τ)=0,(q,∇⋅u)=0.Using the weak divergence free spaceV, the weak formulation (2.4) can be written as follows: find (σ,u)∈(S,V) such that for all (τ,v)∈(S,V)(2.5)(Reut+Reu⋅∇u,v)+(σ,D(v))+2(1-α)(D(u),D(v))=(f,v),λ(σt,τ)+(σ,τ)+λ((u⋅∇)σ,τ)+λ(ga(σ,∇u),τ)-2α(D(u),τ)=0.Existence results for problem (2.1) are proved in [24] for the "slow flow" model of (2.1) (i.e., the u·∇u term in momentum equation is ignored). They are of two kinds: local existence in time of strong solutions in space variable in a C3 domain, and global existence (in time) of a unique solution for u and σ under a small data assumption on f,ft,u0, σ0. 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.SupposeTh is a uniformly regular triangulation of Ω such that Ω={∪K:K∈Th} and assume that there exist positive constants C1,C2 such that C1h⩽hK⩽C2ρK, where hK is the diameter of K, ρK is the diameter of the greatest ball included in K, and h=maxK∈ThhK. Throughout the paper, the constants C,C1,C2,… denote different constants which are independent of mesh size h and time step k. We use the classical Taylor-Hood FE for the approximation in space of (u,p):P2-continuous in velocity, P1-continuous in pressure; we consider P1-continuous approximation of the stress σ. The corresponding FE spaces are (3.1)Xh={v∈X∩C0(Ω)2;v|K∈P2(K)2,∀K∈Th},Sh={τ∈S∩C0(Ω)2×2;τ|K∈P1(K)2×2;∀K∈Th},Qh={q∈Q∩C0(Ω);q|K∈P1(K);∀K∈Th},Vh={v∈Xh;(q,∇⋅v)=0,∀q∈Qh},
where Pw(K) denotes the space of polynomials of degree ⩽w on K∈Th. It is well known [26, 27] that the Taylor-Hood pair (Xh,Qh) 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. LetN be an integer, we divide the interval [0,T] into N intervals of equal length k and denote tn=nk. For notational convenience, we denote vn:=v(·,tn) and (3.2)dtgs:=g(ts)-g(ts-1)k,b(u,σ,τ)=(u⋅∇σ,τ),c(u,v,w)=(u⋅∇v,w).
The following norms are used in the analysis: (3.3)|‖g‖|∞,m:=max⁡1⩽n⩽N⁡‖gn‖m,|‖g‖|0,m:=[∑n=1Nk‖gn‖m2]1/2,
and when ψ(x,t) is defined on the entire time interval (0,T), we denote (3.4)‖ψ‖∞,m:=sup⁡0⩽t⩽T⁡‖φ(⋅,t)‖m,‖ψ‖0,m2:=∫0T‖ψ(⋅,t)‖m2dt,‖ψ‖(t)=‖ψ(⋅,t)‖.Then, based on SUPG formulation, the fully discrete approximating systems of (2.5) is the following;givenuh,0=u0,…,uh,n; σh,0=σ0,…,σh,n, n=0,1,2,…,N-1, find (uh,n+1,σh,n+1)∈(Vh,Sh) such that for all (v,τ)∈(Vh,Sh)(3.5)Re(dtuh,n+1,v)+Rec(uh,n,uh,n+1,v)+2(1-α)(D(uh,n+1),D(v))+(σh,n+1,D(v))=(fn+1,v),λ(dtσh,n+1,τ)+(σh,n+1,τ¯)+λb(uh,n,σh,n+1,τ¯)-2α(D(uh,n+1),τ¯)=-λ(ga(σh,n,∇uh,n),τ¯),
where τ¯=τ+ντn,u,τn,u=uh,n·∇τ, ν 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 ν=0.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(u1h,n+1,σ1h,n+1)∈(Vh,Sh) such that for all (v,τ)∈(Vh,Sh)(4.1)Re(dtu1h,n+1,v)+Rec(u1h,n,u1h,n+1,v)+2(1-α)(D(u1h,n+1),D(v))+(σ1h,n+1,D(v))=(fn+1,v),(4.2)λ(dtσ1h,n+1,τ)+(σ1h,n+1,τ¯)+λ1b(u1h,n,σ1h,n+1,τ¯)-2α(D(u1h,n+1),τ¯)=-λ(ga(σ1h,n,∇u1h,n),τ¯),
where E is chosen such that λ1=λ-Eh>0.Step 2.
Fori=1,2,…, solve the correction problem: find (ui+1h,n+1,σi+1h,n+1)∈(Vh,Sh) such that for all (v,τ)∈(Vh,Sh)(4.3)Re(dtui+1h,n+1,v)+Rec(ui+1h,n,ui+1h,n+1,v)+2(1-α)(D(ui+1h,n+1),D(v))+(σi+1h,n+1,D(v))=(fn+1,v),(4.4)λ(dtσi+1h,n+1,τ)+(σi+1h,n+1,τ¯)+λ1b(ui+1h,n,σi+1h,n+1,τ¯)-2α(D(ui+1h,n+1),τ¯)=-λ(ga(σi+1h,n,∇ui+1h,n),τ¯)-(λ-λ1)b(ui+1h,n,σih,n+1,τ¯).
The initial value approximations are taken to be u1h,0=ui+1h,0=u0, and σ1h,0=σi+1h,0=σ0. 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 (ujh,n+1,σjh,n+1),j=1,2 at each time step n+1. We use the following induction hypothesis, which simply states that the computed iterates (ujh,n+1,σjh,n+1), j=1,2 are bounded independent of h and n:
(4.5)(IH1):‖ujh,s‖∞,‖σjh,s‖∞⩽K,j=1,2,0⩽s⩽n.In next section, we will prove the induction hypothesis(IH1) is right.Lemma 4.1.
Suppose(IH1) is true. For sufficiently small step size k, then Step 1 of Algorithm 1 admits a unique solution (u1h,n+1,σ1h,n+1)∈(Vh,Sh).Proof.
Lettingv=u1h,n+1,τ=σ1h,n+1, multiplying (4.1) by 2α, and adding to (4.2), we get
(4.6)A(u1h,n+1,σ1h,n+1;u1h,n+1,σ1h,n+1)=2α(fn+1,u1h,n+1)+2αRek(u1h,n,u1h,n+1)-λ(ga(σ1h,n,∇u1h,n),σ¯1h,n+1)+λk(σ1h,n,σ1h,n+1),
where the bilinear form A(u,σ;v,τ) is defined by
(4.7)A(u,σ;v,τ)=2αRek(u,v)+2αRec(un,u,v)+4α(1-α)(D(u),D(v))+2α(σ,D(v))+(σ,τ¯)+λk(σ,τ)+λ1b(un,σ,τ)+λ1b(un,σ,νun⋅∇τ)-2α(D(u),σ)-2α(D(u),νun⋅∇τ).
As (4.1)-(4.2) represent a finite system of liner equations, the positivity of A(u1h,n+1,σ1h,n+1;u1h,n+1,σ1h,n+1) is a sufficient condition for the existence and uniqueness of (u1h,n+1,σ1h,n+1).
We now estimate each term ofA(u1h,n+1,σ1h,n+1;u1h,n+1,σ1h,n+1). We have
(4.8)|2αRec(u1h,n,u1h,n+1,u1h,n+1)|≤ϵ1‖D(u1h,n+1)‖2+α2d′K2CK2Re2ϵ1‖u1h,n+1‖2,|(σ1h,n+1,νσ1n,u)|≤‖σ1h,n+1‖02+ν24‖σ1n,u‖2,|λ1b(u1h,n,σ1h,n+1,σ1h,n+1)|≤ϵ2‖σ1n,u‖2+λ124ϵ2‖σ1h,n+1‖2,λ1b(u1h,n,σ1h,n+1,νσ1n,u)=λ1ν‖σ1n,u‖2,|2α(D(u1h,n+1),νσ1n,u)|≤ϵ3‖D(u1h,n+1)‖2+α2ν2ϵ3‖σ1n,u‖2,
where we have used the induction hypothesis (IH1)and the Korn's lemma, CK is Korn constant.
Substituting all above equations intoA(·,·;·,·), we have
(4.9)A(u1h,n+1,σ1h,n+1,u1h,n+1,σ1h,n+1)≥(2αRek-α2d′K2CK2Re2ϵ1)‖u1h,n+1‖2+(λk-λ124ϵ2)‖σ1h,n+1‖2+(4α(1-α)-ϵ1-ϵ3)∥D(u1h,n+1)∥2+(λ1ν-ϵ2-ν24-α2ν2ϵ3)∥σ1n,u∥2.
Choosing ν≤(2λ1(1-α))/(1+3α), k≤min{2λν/λ1,2(1-α)/d′K2CK2Re} and ϵ1=ϵ3=α(1-α), ϵ2=λ1ν/2, the bilinear form A(·,·;·,·) is positive definite. Hence, for k chosen sufficiently small, Step 1 of Algorithm 1 admits a unique solution (u1h,n+1,σ1h,n+1)∈(Vh,Sh).Similarity as Lemma4.1, we can derive that Step 2 of Algorithm 1 admits a unique solution (u2h,n+1,σ2h,n+1)∈(Vh,Sh).


## 5. A Priori Error Estimates

In this section, we explore the error estimates in approximating the true solution(u,p,σ) of (2.1) by the defect step approximation (u1,p1,σ1) and one correction step approximation (u2,p2,σ2). The main results of this section are presented in the following theorem.Theorem 5.1.
Suppose that (2.1) has a solution (u,σ,p)∈C2(0,T;H3)×C2(0,T;H2)×C1(0,T;H2). In addition, assume that k, ν⩽chd′/2, and ∥u∥∞, ∥∇u∥∞, ∥σ∥∞, ∥∇σ∥∞⩽M for all t∈[0,T], (u1h,σ1h) and (u2h,σ2h) satisfy (4.1)-(4.2) and (4.3)-(4.4), respectively. Then, there exists a constant C=C(Re,α,Ω,T,u,p,σ,f,λ1,λ), such that
(5.1)|‖u-u1h‖|∞,0+|‖σ-σ1h‖|∞,0+|‖u-u1h‖|0,1+|‖σ-σ1h‖|0,0⩽B(k,h,ν),(5.2)|‖u-u2h‖|∞,0+|‖σ-σ2h‖|∞,0+|‖u-u2h‖|0,1+|‖σ-σ2h‖|0,0⩽B(k,h,ν),
where
(5.3)B(k,h,ν)=C(h2|‖u‖|0,3+h3|‖ut‖|0,3+h|‖σ‖|0,2+h2|‖σt‖|0,2+h2|‖p‖|0,2)+C(h3|‖u‖|∞,3+h2|‖σ‖|∞,2)+Ck(‖ut‖0,1+‖utt‖0,0+‖σt‖0,1+‖σtt‖0,0)+Cν(|‖σt‖|0,1+|‖σt‖|∞,0)+CT(λ-λ1).Remark 5.2.
As mentioned in Algorithm1, λ-λ1=Eh, thus, for i=1,2, the essence of the estimates (5.1) and (5.2) is
(5.4)|‖u-uih‖|∞,0+|‖σ-σih‖|∞,0+|‖u-uih‖|0,1+|‖σ-σih‖|0,0⩽C(h+k+ν).Before deriving the estimates (5.2), (5.9), we first introduce some notation and some approximating properties of finite element spaces. Let un=u(tn),σn=σ(tn) represent the solution of (2.5) at time tn=nk. Assume (u,σ,p)∈C2(0,T;H3)×C2(0,T;H2)×C1(0,T;H2). Let (ũn,p̃n) denote the Stokes projection of (un,pn) into (Vh,Qh), and σ̃n a Clement interpolant of σn [26, 28]. The following approximating properties are right.(5.5)‖un-ũn‖+h‖∇(un-ũn)‖⩽Ch3‖un‖3,‖pn-p̃n‖⩽Ch2‖pn‖2,‖σn-σ̃n‖+h‖∇(σn-σ̃n)‖⩽Ch2‖σn‖2.Let(u1h,n,σ1h,n) and (u2h,n,σ2h,n) denote the solutions of (4.1)-(4.2) and (4.3)-(4.4), respectively. We denote Φjn,Yjn,ejn,Ψjn,Ujn,ϵjn, j=1,2 as follows:(5.6)Φjn=un-ũn,Yjn=ũn-ujh,n,ejn=Φjn+Yjn=un-ujh,n,Ψjn=σn-σ̃n,Ujn=σ̃n-σjh,n,ϵjn=Ψjn+Ujn=σn-σjh,n.In order to establish Theorem5.1, we need another induction hypothesis, that is, (5.7)(IH2):∑n=0s-1k∥∇Yjn∥∞⩽1,j=1,2.The induction hypothesis(IH2) will be proved later.Proof.
We first give the profile of the proof. The proof of Theorem5.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(IH1) and (IH2) for j=1, we prove the error estimate (5.1) is true.Substep 1.2.
We prove that the two induction hypotheses(IH1) and (IH2) for j=1 are right.Step 2.
Show that the error estimate (5.2) is true.
Now, we start to prove Theorem5.1.
Step1. Prove the error estimate (5.1) is right.
Substep1.1. Under the induction hypotheses (IH1) and (IH2), we prove the error estimate (5.1) is true.
As(σ,u,p) being the exact solution of (2.1), it satisfies the following consistency equation: for all (τ,v)∈Sh×Vh, in particular, at time t=tn+1 (noting that we denote vn+1:=v(·,tn+1)),
(5.8)2αRe(utn+1,v)+(λσtn+1,τ¯)+2αc(Reun+1,un+1,v)+2α(σn+1,D(v))+4α(1-α)(D(un+1),D(v))+λb(un+1,σn+1,τ¯)+(σn+1,τ¯)-2α(D(un+1),τ¯)=2α(fn+1,v)+2α(pn+1,∇⋅v)-λ(ga(σn+1,∇un+1),τ¯).
Equation (5.8) can also be written as
(5.9)2αRe(dtun+1,v)+(λdtσn+1,τ)+2αc(Reu1h,n,un+1,v)+4α(1-α)(D(un+1),D(v))+2α(σn+1,D(v))+λ1b(u1h,n,σn+1,τ¯)+(σn+1,τ¯)-2α(D(un+1),τ¯)=2α(fn+1,v)+2α(pn+1,∇⋅v)-λ(ga(σ1h,n,∇u1h,n),τ¯)+2αRe(dtun+1,v)-2αRe(utn+1,v)+2αc(Reu1h,n,un+1,v)-2αc(Reun+1,un+1,v)+λ(ga(σ1h,n,∇u1h,n),τ¯)-λ(ga(σn+1,∇un+1),τ¯)+(λdtσn+1,τ)-(λσtn+1,τ¯)+λ1b(u1h,n,σn+1,τ¯)-λb(un+1,σn+1,τ¯).
Multiplying (4.1) by 2α and adding to (4.2), for all (τ,v)∈Sh×Vh, we get
(5.10)2αRe(dtu1h,n+1,v)+(λdtσ1h,n+1,τ)+2αc(Reu1h,n,u1h,n+1,v)+2α(σ1h,n+1,D(v))+4α(1-α)(D(u1h,n+1),D(v))+λ1b(u1h,n,σ1h,n+1,τ¯)+(σ1h,n+1,τ¯)-2α(D(u1h,n+1),τ¯)=2α(fn+1,v)-λ(ga(σ1h,n,∇u1h,n),τ¯).
Subtracting (5.9) from (5.10), for all (τ,v)∈Sh×Vh, we obtain the following equation for e1n+1 and ϵ1n+1:
(5.11)2αRe(dte1n+1,v)+(λdtϵ1n+1,τ)+2αc(Reu1h,n,e1n+1,v)+4α(1-α)(D(e1n+1),D(v))+2α(ϵ1n+1,D(v))+λ1b(u1h,n,ϵ1n+1,τ¯)+(ϵ1n+1,τ¯)-2α(D(e1n+1),τ¯)=2αRe(dtun+1-utn+1,v)+2αRec(u1h,n-un+1,un+1,v)-(λσtn+1,τ¯)+(λdtσn+1,τ)-λb(un+1,σn+1,τ¯)+λ1b(u1h,n,σn+1,τ¯)+λ(ga(σ1h,n,∇u1h,n),τ¯)-λ(ga(σn+1,∇un+1),τ¯)+2α(pn+1,∇⋅v).
Substituting e1n+1=Y1n+1+Φ1n+1,ϵ1n+1=U1n+1+Ψ1n+1, v=Y1n+1 and τ=U1n+1 into (5.11), we get
(5.12)2αRe(dtY1n+1,Y1n+1)+λ(dtU1n+1,U1n+1)+2α(U1n+1,D(Y1n+1))+2αc(Reu1h,n,Y1n+1,Y1n+1)+4α(1-α)(D(Y1n+1),D(Y1n+1))+(U1n+1,U¯1n+1)+λ1b(u1h,n,U1n+1,U¯1n+1)-2α(D(Y1n+1),U¯1n+1)=F(Y1n+1,U1n+1),
where
(5.13)F(Y1n+1,U1n+1)=2α(pn+1,∇⋅Y1n+1)+2αRe(dtun+1-utn+1,Y1n+1)-(λσtn+1,U¯1n+1)+2αRec(u1h,n-un+1,un+1,Y1n+1)+λ(ga(σ1h,n,∇u1h,n),U¯1n+1)+(λdtσn+1,U1n+1)-λ(ga(σn+1,∇un+1),U¯1n+1)-2αc(Reu1h,n,Φ1n+1,Y1n+1)-4α(1-α)(D(Φ1n+1),D(Y1n+1))-2α(Ψ1n+1,D(Y1n+1))-(Ψ1n+1,U¯1n+1)-λb(un+1,σn+1,U¯1n+1)+λ1b(u1h,n,σn+1,U¯1n+1)-λ1b(u1h,n,Ψ1n+1,U¯1n+1)+2α(D(Φ1n+1),U¯1n+1)-2αRe(dtΦ1n+1,Y1n+1)-λ(dtΨ1n+1,U1n+1).
Using the identity (a-b,a)=(1/2)(∥a∥2-∥b∥2+∥a-b∥2), we obtain
(5.14)αRek(‖Y1n+1‖2-‖Y1n‖2+‖Y1n+1-Y1n‖2)+4α(1-α)‖D(Y1n+1)‖2+λ2k(‖U1n+1‖2-‖U1n‖2+‖U1n+1-U1n‖2)+‖U1n+1‖2+νλ1‖U1n,u‖2⩽2αRe|c(u1h,n,Y1n+1,Y1n+1)|+ν|(U1n+1,U1n,u)|+λ1|b(u1h,n,U1n+1,U1n+1)|+2αν|(D(Y1n+1),U1n,u)|+F(Y1n+1,U1n+1).
First multiply (5.14) by k and sum it form n=0 to s-1, we obtain
(5.15)αRe∑n=0s-1[‖Y1n+1‖2-‖Y1n‖2+‖Y1n+1-Y1n‖2]+4α(1-α)k∑n=0s-1‖D(Y1n+1)‖2+λ2∑n=0s-1[‖U1n+1‖2-‖U1n‖2+‖U1n+1-U1n‖2]+k∑n=0s-1‖U1n+1‖2+νλ1k∑n=0s-1‖U1n,u‖2⩽k∑n=0s-1[2αRe|c(u1h,n,Y1n+1,Y1n+1)|+ν|(U1n+1,U1n,u)|+λ1|b(u1h,n,U1n+1,U1n+1)|+2αν|(D(Y1n+1),U1n,u)|]+k∑n=0s-1F(Y1n+1,U1n+1).
Then, we get, remarking that Y10=0 and U10=0,
(5.16)αRe[‖Y1s‖2+∑n=0s-1‖Y1n+1-Y1n‖2]+λ2[‖U1s‖2+∑n=0s-1‖U1n+1-U1n‖2]+k∑n=0s-1‖U1n+1‖2+4α(1-α)k∑n=0s-1‖D(Y1n+1)‖2+νλ1k∑n=0s-1‖U1n,u‖2⩽k∑n=0s-1[2αRe|c(u1h,n,Y1n+1,Y1n+1)|+ν|(U1n+1,U1n,u)|+λ1|b(u1h,n,U1n+1,U1n+1)|+2αν|(D(Y1n+1),U1n,u)|]+k∑n=0s-1F(Y1n+1,U1n+1).
For controlling each term on the right-hand side(RHS) of (5.16), the assumption (IH1) and (IH2) are needed here. We will prove the two induction hypotheses in the next subsection. Let us estimate each term of the RHS of (5.16). For details, please see [2, 3, 6]. We start from the first four terms of the RHS of (5.16);
(5.17)2αRe|c(u1h,n,Y1n+1,Y1n+1)|⩽4α2Re2CK2ϵ1‖D(Y1n+1)‖2+d′K24ϵ1‖Y1n+1‖2,ν|(U1n+1,U1n,u)|⩽‖U1n+1‖2+ν24‖U1n,u‖2,|λ1b(u1h,n,U1n+1,U1n+1)|⩽λ12(d′∥∇Y1n∥∞∥U1n+1∥2+M∥U1n+1∥2),2αν|(D(Y1n+1),U1n,u)|⩽4α2ϵ2‖D(Y1n+1)‖2+ν24ϵ2‖U1n,u‖2.
Now, we estimate each term of F(Y1n+1,U1n+1). As F has many terms and in order to make the proof more clear, we divide the terms of F into four parts, that is, c(·,·,·) terms of F, b(·,·,·) terms of F, ga terms of F, other terms of F. Firstly, we estimate c(·,·,·) terms of F;
(5.18)2αRe|c(u1h,n,Φ1n+1,Y1n+1)|⩽2α(Re2‖Y1n+1‖2+K2d′4‖∇Φ1n+1‖2).2αRe|c(u1h,n-un+1,un+1,Y1n+1)|⩽2α(Re2d′2M22‖Y1n‖2+32α‖Y1n+1‖2+Re2d′2M22‖Φ1n‖2+Re22d′2M2k∫tntn+1‖ut‖2dt).
For the b(·,·,·) terms of F(Y1n+1,U1n+1), we have
(5.19)λ1|b(u1h,n,Ψ1n+1,U¯1n+1)|⩽λ12‖U1n+1‖2+ν2‖U1n,u‖2+K2d′2‖∇Ψ1n+1‖2,λ1|(b(u1h,n-un+1,σn+1,U¯1n+1))|⩽λ12‖U1n+1‖2+λ12ν2‖U1n,u‖2+1.5d′3M2‖Y1n‖2+1.5d′3M2‖Φ1n‖2+1.5d′3M2k∫tntn+1‖u‖2dt,|(λ-λ1)b(un+1,σn+1,U¯1n+1)|⩽12(λ-λ1)2M4+‖U1n+1‖2+ν2‖U1n,u‖2.
For the ga terms of F(Y1n+1,U1n+1), we can obtain
(5.20)λ|ga(σ1h,n,∇u1h,n)-ga(σn+1,∇un+1),U¯1n+1|⩽8d′2K2λ2ϵ6‖U1n+1‖2+ϵ3‖∇Y1n‖2+8d′2K2λ2ϵ3ν2‖U1n,u‖2+8d′2K2‖∇Φ1n‖2+5λ2‖U1n+1‖2+5λ2ν2‖F1n,u‖2+8d′2K2λ2k∫tntn+1‖∇ut‖2dt+8d′2M2‖U1n‖2+8d′2M2‖Ψ1n‖2+8d′2M2k∫tntn+1‖σt‖2dt.
Now, we estimate other terms of F(Y1n+1,U1n+1);
(5.21)2α|(pn+1,∇⋅Y1n+1)|⩽2α(CK2ϵ4‖D(Y1n+1)‖2+d′4ϵ4‖pn+1-p̃n+1‖2),2αRe|(dtΦ1n+1,Y1n+1)|⩽2α(Re2‖Y1n+1‖2+14‖dtΦ1n+1‖2),4α(1-α)|D(Φ1n+1),D(Y1n+1)|⩽2α(1-α)(ϵ5‖D(Y1n+1)‖2+‖D(Φ1n+1)‖2ϵ5),2α|(Ψ1n+1,D(Y1n+1))|⩽2α(ϵ6‖D(Y1n+1)‖2+14ϵ6‖Ψ1n+1‖2),2αRe(dtun+1-utn+1,Y1n+1)⩽2α(Re2‖Y1n+1‖2+14‖dtun+1-utn+1‖2),λ|(dtΨ1n+1,U1n+1)|⩽λ2‖U1n+1‖2+14‖dtΨ1n+1‖2,2α|(D(Φ1n+1),U¯1n+1)|⩽‖U1n+1‖2+ν2‖U1n,u‖2+2α2‖∇Φ1n+1‖2,|Ψ1n+1,U¯1n+1|⩽‖U1n+1‖2+ν2‖U1n,u‖2+12‖Ψ1n+1‖2,λ|(dtσn+1-σtn+1,U1n+1)|⩽λ2‖U1n+1‖2+14‖dtσn+1-σtn+1‖2,|(λσtn+1,νU1n,u)|⩽λ2(2+d′‖∇Y1n‖∞)‖U1n+1‖2+ν24(d′2M2+d′‖∇Y1n‖∞)‖σtn+1‖2+ν24K2d′‖∇σtn+1‖2.
Combining the estimates in (5.18)–(5.21), we obtain the estimate for F(Y1n+1,U1n+1):
(5.22)F(Y1n+1,U1n+1)⩽2α(CK2ϵ4+(1-α)ϵ5+ϵ6)‖D(Y1n+1)‖2+CK2ϵ3‖D(Y1n+1)‖2+2α(3Re2+1.5)‖Y1n+1‖2+(3d′M22+αRe2d′2M2)‖Y1n‖2+{d′‖∇Y1n‖∞+3+9λ2+2λ12+8d′2K2λ2ϵ3}‖U1n+1‖2+8d′2M2‖U1n‖2+ν2(5λ2+λ12+4+8d′2K2λ2ϵ3)‖U1n,u‖2+2αd′4ϵ4‖pn+1-p̃n+1‖2+(αRe2d′2M2+3d′3M22)‖Φ1n‖2+2α(K2d′4+1-αϵ5+2α2)‖∇Φ1n+1‖2+8d′2K2‖∇Φ1n‖2+α2‖dtΦ1n+1‖2+14‖dtΨ1n+1‖2+(α2ϵ6+12)‖Ψ1n+1‖2+d′K22‖∇Ψ1n+1‖2+8d′2M2‖Ψ1n‖2+α2‖dtun+1-utn+1‖2+14‖dtσn+1-σtn+1‖2+ν24(d′‖∇Y1n‖∞+d′2M2)⋅‖σtn+1‖2+(αRe2d′2M2k+1.5d′3M2k)×∫tntn+1‖ut‖2dt+ν24K2d′‖∇σtn+1‖2+8d′2M2k∫tntn+1‖σt‖2dt+8d′2K2k∫tntn+1‖∇ut‖2dt+12(λ-λ1)2M4.
Choosing ϵ1=(1-α/12CK2Re3α),ϵ2=(1-α/12α),ϵ3=(2α(1-α)/6CK2),ϵ4=(1-α/6CK2),ϵ5=1/6,ϵ6=(1-α)/6, and substituting (5.17) and (5.22) into (5.16) yields
(5.23)αRe‖Y1s‖2+λ2‖U1s‖2+2α(1-α)∑n=0s-1k‖D(Y1n+1)‖2+k∑n=0s-1‖U1n+1‖2+(νλ1-ν2(7α2(1-α)+28d′2CK2K2λ2α(1-α)+174+5λ2+λ12))k∑n=0s-1‖U1n,u‖2⩽C3k∑n=0s-1‖Y1n+1‖2+k∑n=0s-1(C4(1+‖∇Y1n‖∞))‖U1n+1‖2+C5k∑n=0s-1‖Φ1n+1‖2+C6k∑n=0s-1‖∇Φ1n+1‖2+α2k∑n=0s-1‖dtΦ1n+1‖2+14k∑n=0s-1‖dtΨ1n+1‖2+C7k∑n=0s-1‖Ψ1n+1‖2+C8k∑n=0s-1‖∇Ψ1n+1‖2+α2k∑n=0s-1‖dtun+1-utn+1‖2+14k∑n=0s-1‖dtσn+1-σtn+1‖2+3αd′CK21-αk∑n=0s-1‖pn+1-p̃n+1‖2+k∑n=0s-1ν4(d′2M2+d′‖∇Y1n‖∞)‖σtn+1‖2+k2(αRe2d′2M2+1.5d′3M2)‖ut‖0,02+K2d′ν24|‖σt‖|0,12+8d′2k2(M2‖σt‖0,02+K2‖ut‖0,12)+k∑n=0s-112(λ-λ1)2M4.
According to the approximating properties and the definition (3.5), we can get
(5.24)k∑n=0s-1‖∇Φ1n+1‖2+k∑n=0s-1‖∇Ψ1n+1‖2⩽C(k∑n=0s-1(h4‖un+1‖32+h2‖σn+1‖22))⩽C(h4|‖u‖|0,32+h2|‖σ‖|0,22),(5.25)k∑n=0s-1‖Φ1n+1‖2+k∑n=0s-1‖Ψ1n+1‖2+k∑n=0s-1‖pn+1-p̃n+1‖2⩽C(h6k∑n=0s-1‖un+1‖32+h4k∑n=0s-1‖σn+1‖22+h4k∑n=0s-1‖pn+1‖22)⩽C(h6|‖u‖|0,32+h4|‖σ‖|0,22+h4|‖p‖|0,22),(5.26)k∑n=0s-1‖dtΦ1n+1‖2⩽Ch6|‖ut‖|0,32,(5.27)k∑n=0s-1‖dtΨ1n+1‖2⩽Ch4|‖ut‖|0,22,(5.28)k∑n=0s-1‖dtun+1-utn+1‖2⩽13k2∥utt∥0,02,(5.29)k∑n=0s-1‖dtσn+1-σtn+1‖2⩽13k2∥σtt∥0,02.
Combining the inequalities (5.24)–(5.29) with (5.23), we yield
(5.30)αRe‖Y1s‖2+λ2‖U1s‖2+k∑n=0s-1(2α(1-α)‖D(Y1n+1)‖2+ν2λ1‖U1n,u‖2)⩽Ck∑n=0s-1(‖Y1n+1‖2+‖U1n+1‖2+‖∇Y1n‖∞‖U1n+1‖2)+Cν2|‖σt‖|0,12+ν24k∑n=0s-1‖∇Y1n‖∞‖σtn+1‖2+Ck2(‖ut‖0,12+‖σt‖0,02+‖utt‖0,02+‖σtt‖0,02)+Ch6|‖u‖|0,32+Ch4|‖σ‖|0,22+Ch4|‖p‖|0,22+Ch4|‖u‖|0,32+Ch6|‖ut‖|0,32+Ch2|‖σ‖|0,22+Ch4|‖σt‖|0,22+12(λ-λ1)2TM4,
with ν chosen such that
(5.31)ν⩽λ12(7α2(1-α)+28d′2CK2K2λ2α(1-α)+174+5λ2+λ12)-1.
Applying the induction hypothesis (IH2) and the discrete Gronwall's lemma [29] to (5.30), we have
(5.32)‖Y1s‖2+‖U1s‖2⩽H(k,h,ν),
where
(5.33)H(k,h,ν)=C(h4|‖u‖|0,32+h6|‖ut‖|0,32+h2|‖σ‖|0,22+h4|‖σt‖|0,22+h4|‖p‖|0,22)+Ck2(‖ut‖0,12+‖utt‖0,02+‖σt‖0,12+‖σtt‖0,02)+Cν2(|‖σt‖|0,12+|‖σt‖|∞,02)+C(λ-λ1)2T.
That is to say, under the induction hypotheses(IH1) and (IH2), we establish that the inequality (5.32) is right for all 1⩽s⩽n, and consequently for all n:0⩽n⩽T/k, the inequality (5.32) is right. In the next subsection, we will prove that the induction hypotheses (IH1) and (IH2) are right.
Let us continue to prove the inequality (5.1).
Using triangle inequality, the estimate (5.32), approximation properties (5.5), and (5.24)-(5.25), we have
(5.34)|‖u-u1h‖|∞,02+|‖σ-σ1h‖|∞,02⩽|‖Y1‖|∞,02+|‖Φ1‖|∞,02+|‖U1‖|∞,02+|‖Ψ1‖|∞,02⩽H(k,h,ν)+C(h6|‖u‖|∞,32+h4|‖σ‖|∞,22).
Similarly, by means of triangle inequality, the estimate (5.32) and (5.30), approximation properties (5.5) and (5.24)-(5.25), we can obtain
(5.35)|‖u-u1h‖|0,12+|‖σ-σ1h‖|0,02⩽|‖Y1‖|0,12+|‖Φ1‖|0,12+|‖U1‖|0,02+|‖Ψ1‖|0,02⩽C(T)(H(k,h,ν)+h6|‖u‖|0,32+h4|‖σ‖|0,22).
Thus, we complete the proof of Substep 1.1. That is the inequality (5.1) is right if the induction hypotheses (IH1) and (IH2) are right.
Substep1.2. In this subsection, we will prove that the two induction hypotheses (IH1) and (IH2) for j=1 are right. We first verify the induction hypothesis (IH1).
Assume(IH1) is right for any n=0,1,…,s-1, we will prove that (IH1) is right for n=s. By approximating properties (5.5), inverse inequality [26], and (5.32), we have
(5.36)‖u1h,s‖∞⩽‖u1h,s-us‖∞+‖us‖∞⩽‖Y1s‖∞+‖Φ1s‖∞+M⩽Ch-d′/2‖Y1s‖+Ch-d′/2‖Φ1s‖+M⩽C(kh-d′/2+νh-d′/2+h(1-d′)/2+h(2-d′)/2+h(3-d′)/2+h-d′/2(λ-λ1))+M.
As C(kh-d′/2+νh-d′/2+h(1-d′)/2+h(2-d′)/2+h(3-d′)/2+h-d′/2(λ-λ1)) independent of s, we can choose k,ν⩽hd′/2/C,λ-λ1⩽hd′/2/C, and get ∥u1h,s+1∥∞⩽M+6.
We can get∥σ1h,s+1∥∞⩽M+6 by the same method. Define K=M+6, we confirm (IH1) is right for n=s.
Now we establish that the induction hypothesis(IH2) is right.
Assume(IH2) is right for any n=0,1,…,s-1, we will prove that (IH2) is right for n=s. Using inverse inequality [26], Korn's inequality, and (5.30), we have
(5.37)∑n=0sk∥∇Y1n∥∞=∑n=0s-1k‖∇Y1n+1‖∞(asY10=0)⩽Ch-d′/2∑n=0s-1k‖∇Y1n+1‖⩽Ch-d′/2ks(∑n=0s-1k‖∇Y1n+1‖2)1/2⩽CT(kh-d′/2+νh-d′/2+h(1-d′)/2+h(2-d′)/2+h-d′/2(λ-λ1))⩽1,
with the choosing k,ν,λ-λ1⩽hd′/2/6CT and h(2-d′)/2,h(1-d′)/2,h(2-d′)/2⩽1/6CT.
Step2. We will show that the inequality (5.2) is true.
In order to get the inequality (5.2), we also need induction hypotheses (IH1) and (IH2) for j=2. As the procedure of proof is almost same as Step 1, we only give the different places with Step 1.
Now, we combine the correction problem (4.3)-(4.4) with (5.8) and introduce the approximation error e2n+1,ϵ2n+1. This gives
(5.38)2αRe(dte2n+1,v)+(λdtϵ2n+1,τ)+2αc(Reu2h,n,e2n+1,v)+2α(ϵ2n+1,D(v))+4α(1-α)(D(e2n+1),D(v))+λ1b(u2h,n,ϵ2n+1,τ¯)+(ϵ2n+1,τ¯)-2α(D(e2n+1),τ¯)=2α(pn+1,∇⋅v)+2αRe(dtun+1-utn+1,v)+2αRec(u2h,n-un+1,un+1,v)-(λσtn+1,τ¯)+(λdtσn+1,τ)+(λ-λ1)b(u2h,n,σ1h,n+1-σn+1,τ¯)+(λ-λ1)b(un+1,σn+1,τ¯)+λ(ga(σ2h,n,∇u2h,n),τ¯)-λ(ga(σn+1,∇un+1),τ¯).
Substituting e2n+1=Y2n+1+Φ2n+1,ϵ2n+1=U2n+1+Ψ2n+1, v=Y2n+1 and τ=U2n+1 into (5.38), we obtain
(5.39)αRek(‖Y2n+1‖2-‖Y2n‖2)+λ2k(‖U2n+1‖2-‖U2n‖2)+4α(1-α)(D(Y2n+1),D(Y2n+1))+2α(U2n+1,D(Y2n+1))+(U2n+1,U¯2n+1)+λ1b(u2h,n,U2n+1,U¯2n+1)-2α(D(Y2n+1),U¯2n+1)+2αc(Reu2h,n,Y2n+1,Y2n+1)⩽2α(pn+1,∇⋅Y2n+1)+2αRe(dtun+1-utn+1,Y2n+1)+(λdtσn+1,U2n+1)+2αRec(u2h,n-un+1,un+1,Y2n+1)-λ1b(u2h,n,Ψ2n+1,U¯2n+1)-(λσtn+1,U¯2n+1)+λ(ga(σ2h,n,∇u2h,n),U¯2n+1)-λ(ga(σn+1,∇un+1),U¯2n+1)-2αc(Reu2h,n,Φ2n+1,Y2n+1)-4α(1-α)(D(Φ2n+1),D(Y2n+1))-2α(Ψ2n+1,D(Y2n+1))-(λ-λ1)b(u2h,n,Ψ1n+1,U¯2n+1)+(λ-λ1)b(un+1,σn+1,U¯2n+1)+2α(D(Φ2n+1),U¯2n+1)-(Ψ2n+1,U¯2n+1)-2αRe(dtΦ2n+1,Y2n+1)-(λdtΨ2n+1,U2n+1)-(λ-λ1)b(u2h,n,U1n+1,U¯2n+1).
Comparing (5.39) with (5.12), we find that they are only different on the following terms:
(5.40)-(λ-λ1)b(u2h,n,Ψ1n+1,U¯2n+1)-(λ-λ1)b(u2h,n,U1n+1,U¯2n+1).
So we will only deal with these terms. We have
(5.41)|b(u2h,n,Ψ1n+1,U¯2n+1)|⩽d′k22‖∇Ψ1n+1‖2+‖U2n+1‖2+v2‖U2n,u‖2⩽d′k22ch2‖σn+1‖22+‖U2n+1‖2+v2‖U2n,u‖2,|b(u2h,n,U1n+1,U¯2n+1)|⩽d′k22‖∇U1n+1‖2+‖U2n+1‖2+v2‖U2n,u‖2⩽d′k2ch-22+‖U1n+1‖2+‖U2n+1‖2+v2‖U2n,u‖2,
To conclude, repeat the proof of the first statement of Theorem5.1, replacing u1h,n, σ1h,n, Y1n, Φ1n, U1n and Ψ1n by u2h,n, σ2h,n, Y2n, Φ2n, U2n, and Ψ2n, respectively, using Korn inequality, approximation properties (5.5), (5.24)–(5.29), and the bound ∥Y1n+1∥2+∥U1n+1∥2⩽H(k,h,ν). Hence, we can obtain
(5.42)|‖u-u2h‖|∞,0+|‖σ-σ2h‖|∞,0+|‖u-u2h‖|0,1+|‖σ-σ2h‖|0,0⩽B(k,h,ν).


## 6. Numerical Results

In this section, numerical results for the defect correction method applied to viscoelastic fluid flow are presented using two test problems. The first example is a known analytical solution to verify numerical convergence rates for the defect correction method. The second example simulates viscoelastic flow through a four-to-one contraction flow, a prototypical problem for viscoelastic fluid flow. As mentioned above, continuous piecewise quadratic elements were used for modeling the velocity, and continuous piecewise linear elements were used for the pressure and stress. The constitutive equation was stabilized using an SUPG discretization with parameterν. In this paper, we will not investigate the influence of the parameter ν, thus, we set ν=0.6 h.The defect correction algorithms are implemented using public domain finite element software [30]. Linear systems are solved using the UMFPACK solver. We use the stopping criterion defined by max{∥uih-ui-1h∥,∥σih-σi-1h∥}⩽10-8 for the iterative solver in both the defect step and the correction step of the method. We also set the maximum number of iteration equal to 15.Example 6.1 (Analytical solution).
The theoretical convergence rates were verified by considering fluid flow across a unit square with a known solution. As in [7, 31], Ω=[0,1]2, and chosen functions are added to the right-hand sides of (2.1) such that the true solution to the problem is given by
(6.1)u=(10x2(x-1)2y(y-1)(2y-1)e-t-10x(x-1)(2x-1)y2(y-1)2e-t),p=10(2x-1)(2y-1)e-t,σ=2αD(u).
Letr be the experimental global rate of convergence given by r=log(Er/Er′)/log(h/h′), where h and h′ denote two consecutive mesh sizes with corresponding global errors Er and Er′. In this example, we select Re=1, a=0, α=0.5. The numerical results for Example 6.1 are presented in Tables 1, 2, 3, and 4.Table 1
Defect step approximation forλ=5.0.
1/h|∥ϵ1∥|∞,0r|∥e1∥|∞,0r|∥ϵ1∥|0,0r|∥e1∥|0,1r40.06643—0.003406—0.02316—0.01945—80.019571.760.0011781.530.0063501.870.0048602.00160.0050871.943.181e-41.890.0015981.990.0011742.05320.0014691.798.135e-51.974.268e-41.902.937e-42.00Table 2
Correction step approximation forλ=5.0.
1/h|∥ϵ2∥|∞,0r|∥e2∥|∞,0r|∥ϵ2∥|0,0r|∥e2∥|0,1r40.06643—0.003410—0.02316—0.01945—80.019571.760.0011801.530.0063471.870.0048612.00160.0050871.943.191e-41.890.0015902.000.0011732.05320.0012841.998.151e-51.973.995e-41.992.917e-42.01Table 3
Defect step approximation forλ=1000.
1/h|∥ϵ1∥|∞,0r|∥e1∥|∞,0r|∥ϵ1∥|0,0r|∥e1∥|0,1r40.06710—0.003523—0.02352—0.01959—80.019621.770.0012101.540.0064191.870.0048852.00160.0050901.953.255e-41.890.0016131.990.0011782.05320.0014931.778.319e-51.974.306e-41.912.949e-41.99Table 4
Correction step approximation forλ=1000.
1/h|∥ϵ2∥|∞,0r|∥e2∥|∞,0r|∥ϵ2∥|0,0r|∥e2∥|0,1r40.06710—0.003526—0.02352—0.01959—80.019621.770.0012131.540.0064161.870.0048862.00160.0050901.953.266e-41.890.0016052.000.0011782.05320.0012841.998.336e-51.974.030e-41.992.927e-42.01To reduce the influence of the time discretization error, the time step is taken to be very small:k=O(h2).Forλ=5, λ=1000 and the final time T=0.1, the calculated convergence rates in Tables 1–4 confirm what is predicted by Theorem 5.1 for continuous (P2,P1,P1) elements in space. In fact, our numerical convergence rates are better than the theoretical ones. We will find the reason in future work.Example 6.2 (Four-to-one contraction flow).
Numerical simulations of viscoelastic flow through a planar or axisymmetric contraction have been widely studied [32, 33]. Here the case of planar flow through a contraction geometry with a ratio of 4 : 1 with respect to upstream and downstream channel widths is considered. The contraction angle is a fixed 3π/2 and the channel lengths are sufficiently long to impose a fully developed Poiseuille flow in the inflow and outflow channels. The geometry of the computational domain is illustrated in Figure 1. The lower left corner of the domain corresponds to x=y=0. The computations of the mesh are also shown in Figure 1 with Δxmin=0.0625 and Δymin=0.015625. We denote Γin={(x,y):x=0,0⩽y⩽1} and Γout={(x,y):x=8,0⩽y⩽0.25}. On this domain the velocity boundary conditions are u1=(1/32)(1-y2),u2=0,onΓin and u1=2((1/16)-y2),u2=0,onΓout. On Γin, specified boundary conditions for σ are given as follows: σ11=(-αλ(a+1)(-y/16)2)/((a2-1)λ2(-y/16)2-1),  σ12=σ21=-α(-y/16)/((a2-1)λ2(-y/16)2-1),  σ22=-αλ(a-1)(-y/16)2/((a2-1)λ2(-y/16)2-1). Symmetry conditions are imposed on the bottom of the computational domain. In this example, the parameters Re, α, λ, and a are set to 1, 8/9, 1.3, and 0, respectively.Figure 1
Plot of4:1 contraction domain geometry and sample mesh.We performed the following study: starting from rest, we measured the time that the approximation solution reaches a steady state. The criterion to stop this process is the following:(6.2)max⁡⁡{‖uhn+1-uhn‖L2(Ω)‖uhn+1‖L2(Ω),‖σhn+1-σhn‖L2(Ω)‖σhn+1‖L2(Ω)}⩽10-5,
where n+1,n denote tn+1,tn, respectively.In Figure2, we plot the evolution of the kinetic energy ∥uhn+1∥02/2 and ∥σhn+1∥02/2 using time step k=0.1 until it reaches its steady state, where we observe they convergence towards a steady state and also the absence of oscillations along the process.Evolution of∥uhn+1∥02/2 (a) and ∥σhn+1∥02/2 (b) in time with time step k=0.1.
(a)(b)Figure3 presents the horizontal and vertical velocities near the reentrance corner along the vertical line x=4.0625 for λ=1.3. We observe that the horizontal velocity is almost continuous, while the vertical velocity has high gradients near y=0.11 and y=0.23 from Figure 3. However, we find that the solutions of the time-dependent problem can converge to the solutions of the steady problem. Figure 4 present the streamlines of the fluid with λ=1.3. We plot the streamlines for both the steady problem and the time-dependent problem at final time t=24.5. It is easy to observe that these two figures are almost alike.Horizontal (a) and vertical (b) velocity near reentrant corner. The marks “+” indicate results for steady problem and “–” indicate results for the time-dependent problem at timet=24.5 using time step k=0.1.
(a)(b)Streamlines and magnitude of velocity contours foru. (a): steady problem, (b): time-dependent problem at time t=24.5 using k=0.1.
(a)(b)


## 7. Conclusions

In this paper, we present a defect correction mixed finite element method for solving the time-dependent Johnson-Segalman viscoelastic equations. A priori error estimates for the defect step and the first correction step of the defect correction method are derived. Finally, we present two numerical examples. One is a known problem and the other is a benchmark problem.For simplicity, in this paper, we choose the parameterν=0.6 h. More appropriate choice of parameter ν is currently under investigation. Further developments should be extending the method to other non-Newtonian flow or other discrete scheme.



---

*Source: 689804-2011-07-10.xml*