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

A defect-correction mixed finite element method for solving the time-dependent JohnsonSegalman 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 PetrovGalerkin 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.


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 of u is a scalar ∇u u i,i .with ∇u ij u i,j .For a tensor function τ and a vector function u, divergence of τ is a vector ∇ • τ with ∇ • τ i τ ij,j and u • ∇ denotes the operator u i ∂/∂x i .
Let T > 0 be a given final time, we consider the time-dependent Johnson-Segalman viscoelastic equations as follows: 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 u 1 , u 2 , 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 ∇u T is the rate of the strain tensor, for all a ∈ −1, 1 , g a σ, ∇u is defined by g a σ, ∇u 1 − a 2 σ∇u ∇u T σ − 1 a 2 ∇u σ σ ∇u T .

2.2
Note that an Oldroyd B constitutive model is obtained when a 1 in g a σ, ∇u .The L 2 Ω inner product is denoted by •, • , and the L p Ω norm by • L p , with the special cases of L 2 Ω and L ∞ Ω norms being written as • and • ∞ .For k ∈ N, we denote the norm associated with the Sobolev space W m,p Ω by • W m,p , with the special case W m,2 Ω being written as H m Ω 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
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 q, ∇ • u 0.

2.4
Using the weak divergence free space V , the weak formulation 2.4 can be written as follows: find σ, u ∈ S, V such that for all τ, v ∈ S, V

2.5
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 C 3 domain, and global existence in time of a unique solution for u and σ under a small data assumption on f, f t , u 0 , σ 0 .For a more complete discussion of existence and uniqueness issues, see 25 .

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 T h is a uniformly regular triangulation of Ω such that Ω {∪K : K ∈ T h } and assume that there exist positive constants C 1 , C 2 such that C 1 h h K C 2 ρ K , where h K is the diameter of K, ρ K is the diameter of the greatest ball included in K, and h max K∈T h h K .Throughout the paper, the constants C, C 1 , C 2 , . . .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 : P 2 -continuous in velocity, P 1 -continuous in pressure; we consider P 1 -continuous approximation of the stress σ.The corresponding FE spaces are

3.1
where P w K denotes the space of polynomials of degree w on K ∈ T h .It is well known 26, 27 that the Taylor-Hood pair X h , Q h 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 N be an integer, we divide the interval 0, T into N intervals of equal length k and denote t n nk.For notational convenience, we denote v n : v •, t n and The following norms are used in the analysis: and when ψ x, t is defined on the entire time interval 0, T , we denote Then, based on SUPG formulation, the fully discrete approximating systems of 2.5 is the following; given u h,0 u 0 , . . ., u h,n ; σ h,0 σ 0 , . . ., σ h,n , n 0, 1, 2, . . ., N − 1, find where τ τ ντ n,u , τ n,u u h,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 .

Defect Correction Method
In this section, our defect correction method used in computing the solution to 3.5 is described as follows.
Step 1. Solve the defected problem: find u where E is chosen such that λ 1 λ − Eh > 0.

4.4
The initial value approximations are taken to be u h,0 1 u h,0 i 1 u 0 , and σ h,0 1 σ h,0 i 1 σ 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 u h,n 1 j , σ h,n 1 j , j 1, 2 at each time step n 1.We use the following induction hypothesis, which simply states that the computed iterates u h,n 1 j , σ h,n 1 j , j 1, 2 are bounded independent of h and 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 , multiplying 4.1 by 2α, and adding to 4.2 , we get where the bilinear form A u, σ; v, τ is defined by

4.7
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 u h,n 1 1 , σ h,n 1 1 .

We now estimate each term of
where we have used the induction hypothesis IH1 and the Korn's lemma, C K is Korn constant.Substituting all above equations into A •, •; •, • , we have 4.9

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 u 1 , p 1 , σ 1 and one correction step approximation u 2 , p 2 , σ 2 .The main results of this section are presented in the following theorem.

5.3
Remark 5.2.As mentioned in Algorithm 1, λ − λ 1 Eh, thus, for i 1, 2, 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 u n u t n , σ n σ t n represent the solution of 2.5 at time and σ n a Clement interpolant of σ n 26, 28 .The following approximating properties are right.

5.5
Let u h,n 1 , σ h,n 1 and u h,n 2 , σ h,n 2 denote the solutions of 4.1 -4.2 and 4.3 -4.4 , respectively.We denote Φ n j , Y n j , e n j , Ψ n j , U n j , n j , j 1, 2 as follows:

5.6
In order to establish Theorem 5.1, we need another induction hypothesis, that is, The induction hypothesis IH2 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 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 Theorem 5.1.
Step 1. Prove the error estimate 5.1 is right.Substep 1.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 ∈ S h × V h , in particular, at time t t n 1 noting that we denote

5.8
Equation 5.8 can also be written as 5.9 Multiplying 4.1 by 2α and adding to 4.2 , for all τ, v ∈ S h × V h , we get

5.10
Subtracting 5.9 from 5.10 , for all τ, v ∈ S h × V h , we obtain the following equation for e n 1 1 and n 1 1 :

5.11
Substituting where

5.13
Using the identity a − b, a

5.14
First multiply 5.14 by k and sum it form n 0 to s − 1, we obtain

5.15
Then, we get, remarking that Y 0 1 0 and

5.16
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 .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, g a terms of F, other terms of F. Firstly, we estimate c •, •, • terms of F;

5.18
For the b

5.19
For the g a terms of

Discrete Dynamics in Nature and Society
According to the approximating properties and the definition 3.5 , we can get

5.30
with ν chosen such that

5.31
Applying the induction hypothesis IH2 and the discrete Gronwall's lemma 29 to 5.30 , we have where

5.33
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
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
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.Substep 1.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
As 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 Step 2. 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 e n 1 2 , n 1 2 .This gives

5.39
Comparing 5.39 with 5.12 , we find that they are only different on the following terms: .

5.40
So we will only deal with these terms.We have

5.41
To conclude, repeat the proof of the first statement of Theorem 5.

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{ 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 , σ 2αD u .

6.1
Let r 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.
To reduce the influence of the time discretization error, the time step is taken to be very small: k O h 2 .
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 P 2 , P 1 , P 1 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    −αλ a − 1 −y/16 2 / a 2 − 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.
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: max where n 1, n denote t n 1 , t n , respectively.In Figure 2, we plot the evolution of the kinetic energy u n 1 h 2 0 /2 and σ n 1 h 2 0 /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.
Figure 3 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.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.

Conclusions
In this paper, we present a defect correction mixed finite element method for solving the timedependent 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.

Figure 2 :of u n 1 h 2 0 /2 a and σ n 1 h 2 0
Figure 2: Evolution of u n 1 h Figure 4 present the streamlines of the fluid with λ 1.3.

Figure 3 :Figure 4 :
Figure 3: 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 time t 24.5 using time step k 0.1.