Discrete Dynamics in Nature and SocietyVolume 2011 (2011), Article ID 689804, 25 pagesdoi:10.1155/2011/689804

Research Article

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

Yunzhang Zhang,1,2 Yanren Hou,¹ and Ganshan Yang3,4

¹School of Science, Xi'an Jiaotong University, Xi'an 710049, China²School of Mathematics and Statistics, Henan University of Science and Technology, Luoyang 471003, China³Institute of Mathematics Science, Yunnan Normal University, Kunming 650092, China⁴Department of Mathematics, Yunnan Nationalities University, Kunming 650031, China

Received 22 November 2010; Revised 13 March 2011; Accepted 29 April 2011

Academic Editor: M. De la Sen

Copyright © 2011 Yunzhang Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

Table 1: Defect step approximation for .

1/h

r

r

r

r

4
0.06643
—
0.003406
—
0.02316
—
0.01945
—

8
0.01957
1.76
0.001178
1.53
0.006350
1.87
0.004860
2.00

16
0.005087
1.94

1.89
0.001598
1.99
0.001174
2.05

32
0.001469
1.79

1.97

1.90

2.00

Table 2: Correction step approximation for .

1/h

r

r

r

r

4
0.06643
—
0.003410
—
0.02316
—
0.01945
—

8
0.01957
1.76
0.001180
1.53
0.006347
1.87
0.004861
2.00

16
0.005087
1.94

1.89
0.001590
2.00
0.001173
2.05

32
0.001284
1.99

1.97

1.99

2.01

Table 3: Defect step approximation for .

1/h

r

r

r

r

4
0.06710
—
0.003523
—
0.02352
—
0.01959
—

8
0.01962
1.77
0.001210
1.54
0.006419
1.87
0.004885
2.00

16
0.005090
1.95

1.89
0.001613
1.99
0.001178
2.05

32
0.001493
1.77

1.97

1.91

1.99

Table 4: Correction step approximation for .

1/h

r

r

r

r

4
0.06710
—
0.003526
—
0.02352
—
0.01959
—

8
0.01962
1.77
0.001213
1.54
0.006416
1.87
0.004886
2.00

16
0.005090
1.95

1.89
0.001605
2.00
0.001178
2.05

32
0.001284
1.99

1.97

1.99

2.01

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.

Figure 1: Plot of 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: 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 2: Evolution of (a) and (b) in time with time step .

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)

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 using time step .

(a)

(b)

Figure 4: Streamlines and magnitude of velocity contours for . (a): steady problem, (b): time-dependent problem at time using .

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

References

J. Baranger and S. Wardi, “Numerical analysis of a FEM for a transient viscoelastic flow,” Computer Methods in Applied Mechanics and Engineering, vol. 125, no. 1–4, pp. 171–185, 1995.

V. J. Ervin and W. W. Miles, “Approximation of time-dependent viscoelastic fluid flow: SUPG approximation,” SIAM Journal on Numerical Analysis, vol. 41, no. 2, pp. 457–486, 2003.

V. J. Ervin and W. W. Miles, “Approximation of time-dependent, multi-component, viscoelastic fluid flow,” Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 18-20, pp. 2229–2255, 2005.

T. J. Hughes and A. N. Brooks, A multidimensional Upwind Scheme with No Crosswind Diffusion, in Finite Element Methods for Convection Dominated Flows, New York, NY, USA, 1979.

M. Bensaada and D. Esselaoui, “Error estimates for a stabilized finite element method for the Oldroyd B model,” Journal of Mathematical Analysis and Applications, vol. 325, no. 2, pp. 1042–1059, 2007.

V. J. Ervin and N. Heuer, “Approximation of time-dependent, viscoelastic fluid flow: Crank-Nicolson, finite element approximation,” Numerical Methods for Partial Differential Equations, vol. 20, no. 2, pp. 248–283, 2004.

J. C. Chrispell, V. J. Ervin, and E. W. Jenkins, “A fractional step method approximation of time-dependent viscoelastic fluid flow,” Journal of Computational and Applied Mathematics, vol. 232, no. 2, pp. 159–175, 2009.

A. Bonito, P. Clément, and M. Picasso, “Mathematical and numerical analysis of a simplified time-dependent viscoelastic flow,” Numerische Mathematik, vol. 107, no. 2, pp. 213–255, 2007.

P. Nithiarasu, “A fully explicit characteristic based split (CBS) scheme for viscoelastic flow calculations,” International Journal for Numerical Methods in Engineering, vol. 60, no. 5, pp. 949–978, 2004.

H. K. Lee, “Analysis of a defect correction method for viscoelastic fluid flow,” Computers and Mathematics with Applications, vol. 48, no. 7-8, pp. 1213–1229, 2004.

V. J. Ervin and H. Lee, “Defect correction method for viscoelastic fluid flows at high Weissenberg number,” Numerical Methods for Partial Differential Equations, vol. 22, no. 1, pp. 145–164, 2006.

V. J. Ervin, J. S. Howell, and H. Lee, “A two-parameter defect-correction method for computation of steady-state viscoelastic fluid flow,” Applied Mathematics and Computation, vol. 196, no. 2, pp. 818–834, 2008.

X. L. Luo, “An incremental difference formulation for viscoelastic flows and high resolution FEM solutions at high Weissenberg numbers,” Journal of Non-Newtonian Fluid Mechanics, vol. 79, no. 1, pp. 57–75, 1998.

J. Petera, “A new finite element scheme using the Lagrangian framework for simulation of viscoelastic fluid flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 103, no. 1, pp. 1–43, 2002.

J. S. Howell, “Computation of viscoelastic fluid flows using continuation methods,” Journal of Computational and Applied Mathematics, vol. 225, no. 1, pp. 187–201, 2009.

K. Bhmer and H. J. Stetter, Defect Correction Methods-Theory and Applications, vol. 5 of Computing Supplementum, Springer, Vienna, Austria, 1984.

R. Minero, M. J. H. Anthonissen, and R. M. M. Mattheij, “A local defect correction technique for time-dependent problems,” Numerical Methods for Partial Differential Equations, vol. 22, no. 1, pp. 128–144, 2006.

V. J. Ervin, W. J. Layton, and J. M. Maubach, “Adaptive defect-correction methods for viscous incompressible flow problems,” SIAM Journal on Numerical Analysis, vol. 37, no. 4, pp. 1165–1185, 2000.

W. Layton, H. Lee, and J. A. Peterson, “A defect-correction method for the incompressible Navier-Stokes equations,” Applied Mathematics and Computation, vol. 129, no. 1, pp. 1–19, 2002.

Q. Liu and Y. R. Hou, “A two-level defect-correction method for Navier-Stokes equations,” Bulletin of the Australian Mathematical Society, vol. 81, no. 3, pp. 442–454, 2010.

A. Labovschii, “A defect correction method for the time-dependent Navier-Stokes equations,” Numerical Methods for Partial Differential Equations, vol. 25, no. 1, pp. 1–25, 2009.

Y. Zhang, Y. Hou, and B. Mu, “Defect correction method for time-dependent viscoelastic fluid flow,” International Journal of Computer Mathematics, vol. 88, no. 7, pp. 1546–1563, 2011.

R. B. Bird, R.C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids, John Wiley and Sons, New York, NY, USA, 1987.

C. Guillopé and J. C. Saut, “Existence results for the flow of viscoelastic fluids with a differential constitutive law,” Nonlinear Analysis. Theory, Methods & Applications, vol. 15, no. 9, pp. 849–869, 1990.

M. Renardy, Mathematical Analysis of Viscoelastic Flows, vol. 73 of CBMS-NSF Regional Conference Series in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pa, USA, 2000.

S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, vol. 15 of Texts in Applied Mathematics, Springer, New York, NY, USA, 1994.

V. Girault and P. A. Raviart, Finite Element methods for Navier-Stokes Equations, vol. 5 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1986.

P. Clément, “Approximation by finite element functions using local regularization,” RAIRO Analyse Numerique, vol. 9, no. R-2, pp. 77–84, 1975.

J. G. Heywood and R. Rannacher, “Finite-element approximation of the nonstationary Navier-Stokes problem. IV. Error analysis for second-order time discretization,” SIAM Journal on Numerical Analysis, vol. 27, no. 2, pp. 353–384, 1990.

F. Hecht, O. Pironneau, A. Le Hyaric, and K. Ohtsuka, “FreeFem++,” 2009, http://www.freefem.org/ff++.

H. Lee, “A multigrid method for viscoelastic fluid flow,” SIAM Journal on Numerical Analysis, vol. 42, no. 1, pp. 109–129, 2004.

F. P.T. Baaijens, “Mixed finite element methods for viscoelastic flow analysis: a review,” Journal of Non-Newtonian Fluid Mechanics, vol. 79, no. 2-3, pp. 361–385, 1998.

R. G. Owens and T. N. Phillips, Computational Rheology, Imperial College Press, London, UK, 2002.