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
According to the approximating properties and the definition (3.5), we can get
Combining the inequalities (5.24)–(5.29) with (5.23), we yield
with chosen such that
Applying the induction hypothesis and the discrete Gronwall's lemma [29] to (5.30), we have
where
That is to say, under the induction hypotheses and , we establish that the inequality (5.32) is right for all , and consequently for all , the inequality (5.32) is right. In the next subsection, we will prove that the induction hypotheses and 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
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
Thus, we complete the proof of Substep 1.1. That is the inequality (5.1) is right if the induction hypotheses and are right.
Substep 1.2. In this subsection, we will prove that the two induction hypotheses and for are right. We first verify the induction hypothesis .
Assume is right for any , we will prove that is right for . By approximating properties (5.5), inverse inequality [26], and (5.32), we have
As independent of , we can choose , and get .
We can get by the same method. Define , we confirm is right for .
Now we establish that the induction hypothesis is right.
Assume is right for any , we will prove that is right for . Using inverse inequality [26], Korn's inequality, and (5.30), we have
with the choosing and .
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 and for . 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 . This gives
Substituting , and into (5.38), we obtain
Comparing (5.39) with (5.12), we find that they are only different on the following terms:
So we will only deal with these terms. We have
To conclude, repeat the proof of the first statement of Theorem 5.1, replacing , , , , and by , , , , , and , respectively, using Korn inequality, approximation properties (5.5), (5.24)–(5.29), and the bound . Hence, we can obtain
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 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 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], , and chosen functions are added to the right-hand sides of (2.1) such that the true solution to the problem is given by
Let be the experimental global rate of convergence given by , where and denote two consecutive mesh sizes with corresponding global errors and . In this example, we select , , . 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: .
For , and the final time , the calculated convergence rates in Tables 1–4 confirm what is predicted by Theorem 5.1 for continuous 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 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 . The computations of the mesh are also shown in Figure 1 with and . We denote and . On this domain the velocity boundary conditions are and . On , specified boundary conditions for are given as follows: , , . Symmetry conditions are imposed on the bottom of the computational domain. In this example, the parameters , , , and are set to 1, , 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: where denote , respectively.
In Figure 2, we plot the evolution of the kinetic energy and using time step until it reaches its steady state, where we observe they convergence towards a steady state and also the absence of oscillations along the process.
(a)
(b)
Figure 3 presents the horizontal and vertical velocities near the reentrance corner along the vertical line for . We observe that the horizontal velocity is almost continuous, while the vertical velocity has high gradients near and 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 . We plot the streamlines for both the steady problem and the time-dependent problem at final time . It is easy to observe that these two figures are almost alike.
(a)
(b)
(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 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.
Acknowledgments
The authors would like to thank the editor and the anonymous referees for their criticism and valuable comments, which led to the improvement of this paper. The authors want to thank Professor J. S. Howell for meaty the discussion of part program code. This paper was supported by the National Natural Science Foundation of China (no. 10871156, 10861014) and the Fund of Xi'an Jiaotong University (no. 2009xjtujc30).