- Email: [email protected]
A preconditioned generalized minimal residual algorithm for the numerical solution of viscoelastic fluid flows
277
Journal of Non-Newtonian Fluid Mechanics, 36 (1990) 277-288 Elsevier Science Publishers B.V., Amsterdam
A PRECONDITIONED GENERALIZED ALGORITHM FOR THE NUMERICAL OF VISCOELASTIC FLUID FLOWS
ANDRI?
MINIMAL RESIDUAL SOLUTION
FORTIN
Department of Applied Mathematics, .&ole Polytechnique Succursale A, Montrkal H3C 3A7 (Canada) MICHEL Department (Received
de MontrPal C. P. 6079,
FORTIN of Applied Mathematics,
UniversitP Lava1 Quebec GIK 7P4 (Canada)
August 24, 1989; in revised form January
10, 1990)
Abstract The numerical simulation of the flow of a viscoelastic fluid by a generalized minimal residual (GMRES) method is investigated and tested to determine its convergence properties. A method based on the LesaintRaviart technique is used as a preconditioner for the GMRES method to enhance the convergence rate and robustness of the algorithm. Numerical results are presented for the stick-slip problem. Keywords: decoupling; finite elements; GMRES; Oldroyd-B Model; viscoelasticity
1. Introduction Numerical simulations of viscoelastic fluid flows require the solution of very large strongly non-linear systems of equations. This is the main reason Ghy most successful numerical methods are based on the Newton-Raphson algorithm which offers robustness and quadratic convergence in most cases. This means that the velocities, pressure and stresses are obtained as the solution of a global system of non-linear equations that is solved by a Newton method. This requires the assembly and factorization of huge systems of equations that are computationally expensive. This is why we believe that new iterative schemes based on decoupled methods are needed and should receive more attention. In contrast, decoupled methods usually need the resolution of a larger number of smaller systems of equations. The basic idea is to solve separately 0377-0257/90/$03.50
0 1990 - Elsevier Science Publishers
B.V.
278 the constitutive equation for the stresses and the conservation equations for the velocities and pressure. Because of the strong non-linearities of viscoelastic fluid flow problems, decoupled methods have not been very successful. Convergence difficulties are quickly encountered (that is at relatively low Weissenberg numbers), precluding their use in practical situations. It is, however, generally accepted [l] that appropriate decoupled methods must be developed for the solution of industrial bidimensional applications. This statement is even stronger for three-dimensional problems. The difficulty in designing a robust decoupled method is to avoid constructing a large global matrix while keeping the good convergence properties of the Newton method. We present in this paper a generalized minimal residual (GMRES) algorithm using a decoupled method described by Fortin and Fortin [2] as a preconditioner. The resulting algorithm is proved to be efficient but of course not as robust as the full Newton method. Numerical results are presented for the stick-slip problem showing convergence of the iterative scheme up to De = 4. 2. Statement of the problem 2. I. Basic equations We briefly recall in this section the statement of the problem. An Oldroyd B model is used as a constitutive equation relating the stress tensor to the velocity field. This model has been chosen for the purpose of comparison with existing results only. The numerical solution of the flow of a viscoelastic fluid requires the computation of the velocity field u, the pressure p and the stress tensor 7 satisfying the following partial differential equations. Momentum equations: *[2q2i)(u)]
-v
-v
*7+vp=.f.
(1)
Mass conservation: V *u=o. Constitutive
(2) equation
7+x[u.v7-vu.7-T~(vu)t] where
(Oldroyd
B model): =2Tjtj+),
(3)
279
One of the major difficulties in solving the system (l)-(3) is the numerical treatment of the constitutive equation (3). The hyperbolicity of this equation requires special care as clearly demonstrated by Marchal and Crochet [3] and Fortin and Fortin [2]. Indeed, some kind of upwinding seems necessary to stabilize most numerical algorithms. In Ref. 3, a streamline upwind method is presented while the Lesaint-Raviart method is prefered in Ref. 2. This last method will be used again and we refer to Ref. 2 for a detailed description. 2.2. Discretization It is now accepted that the discretization of the system of equations (l)-(3) requires special attention. The difficulties come from the (u, p, T) formulation of the basic equations. In the Newtonian case, it is possible to discard 7 from the equations and we are left with the classical velocity-pressure formulation. Even in this simpler case, it is well known that the discretizations of the velocities and pressure must satisfy the so-called Brezzi condition (see Ref. 4). For viscoelastic fluids, it is not possible to discard r (unless h = 0) and consequently it is necessary to develop consistent (u, p, 7) discretizations. In a recent paper, Fortin and Pierre [5] have studied this particular formulation in the case of the element introduced by Marchal and Crochet [3]. Their analysis allowed them to justify the 4 by 4 subdivision of the elements for the stresses. It is possible to generalize their analysis to develop various consistent discretizations. We shall concentrate on discretizations involving discontinuous approximations for the stresses. In that case, Fortin and Pierre came to the conclusion that if V,, Q,, and IV, denote respectively the discrete spaces for the velocities, pressure and stresses, then V, and W, must satisfy the condition vv,c
w,,
(4)
i.e. the gradient of any function belonging to V, must be included in the space W,. Condition (4) is complementary to the Brezzi condition relating V, and Qh. In this paper, we use a continuous biquadratic (nine nodes) approximation for the velocities and a piecewise linear pressure (no continuity is required at element interfaces for the pressure). For the stresses, we set IV, = VV,, i.e. the space IV, consists of the gradient of all the functions belonging to V,. In other words, IV, is made up of incomplete biquadratic (eight nodes) polynomials. Since the functions in V, are only continuous (not C’) at element interfaces, the gradient of such functions will be discontinuous on
280 those interfaces. This is illustrated in Fig. 3 by setting the stress nodes inside the element. We now give the variational formulation for the system (l)-(3). Find (u, p, T) E V,?x Q,, X W, such that
=
JK
29,?(u):
cp dx dy
for every element
K and V’cpE V$.
(7)
Problem (7) can be solved on an element-by-element basis, thanks to the discontinuous approximation of the stresses and to the Lesaint-Raviart method. aK_ (aK+) denotes the inflow (outflow) part of the boundary of element K and [T] is the jump in 7 across the element boundary (see Ref. 2 for a complete justification). The parameter CYis the upwinding factor, usually taking a value between 0 and 1. However, as we shall see, a larger value can help to stabilize the numerical results. 3. Numerical method 3.1. Decoupled methods The basis of most decoupled methods is summarized by the following steps: 1. Solve the constitutive equation (7) for 7 using the velocity field obtained from a previous iteration. 2. Solve (5) and (6) for u and p using the value of 7 calculated in step 1. 3. Check for convergence and return to step 1 if necessary. Many variants of the above algorithm are possible. Most of them are based on the method of characteristics [6-81. As already mentioned, the major drawback of this method is the loss of convergence of the iterative scheme at relatively low Weissenberg numbers. For example, when using the Lesaint-Raviart method, the above algorithm failed to converge for any value of the Weissenberg number. A steady solution could, however, be obtained as a long time solution of transient simulations. For a complete
281 review of these difficulties, we refer to Keunings [9]. It is worth mentioning that appropriate treatment of the constitutive equation improves the convergence properties of decoupled methods as clearly demonstrated by Luo and Tanner [lo]. 3.2. The generalized
minimal residual method
The GMRES method was first introduced by Saad and Schultz [ll] for the resolution of non-symmetric linear systems. This method was then extended to the solution of large non-linear systems of equations. The GMRES algorithm has recently received a lot of attention. The reason for this is obvious. In its simplest form, this method requires no assembly and factorization of a global matrix. Only the residual has to be computed and stored to construct a Krylov space. The economy in memory space can be very important and depends on the dimension of that Krylov space. Let us recall the main features of the GMRES algorithm which. as we will verify, can be seen as a quasi-Newton method. We refer the reader to Ref. 11 for a more complete discussion. The basic idea is to start with the Newton algorithm for the solution of the non-linear system: R(X)=O.
(8)
The Newton algorithm is then X0 being given For nk0: 1. solve J, AX = -R(X,,) (J, 2. update Xn+l = X, + AX, 3. check for convergence.
is the Jacobian matrix at X,),
The linear system in step 1 is usually solved by LU decomposition. In the case of viscoelastic fluid flow problems, this matrix is extremely large and should be constructed and factorized at each iteration. The main feature of the GMRES algorithm is to avoid the construction of J,. To achieve this goal, two important remarks must be made. Firstly, iterative methods for solving a linear system of the form AX = b such as the conjugate gradient method (for symmetric matrices) and the minimal residual method do not require the matrix A itself but only the product of A with different direction vectors d. Moreover, in the case of a Jacobian matrix, we have lim
Jd=
n
h-0
WC,+hd) -R(X) h
(9)
282
Use of eqn. (9) eliminates the need for the explicit construction of J,. Finally, it is worth noticing that the linear system in step 1 does not need to be solved exactly. This results in a quasi-Newton method since the Jacobian matrix is only approximated by the finite difference (9) and the related linear systems are not solved exactly. We refer to Ref. 11 for a detailed description of the construction of the Krylov space for the resolution of linear systems by the GMRES method. 3.3. Application to viscoelastic jluid flow problems As already mentioned, viscoelastic fluid flow problems are strongly nonlinear and it is not likely that the GMRES method will converge, at least in its simplest form. However, a good preconditioner could possibly enhance the convergence of the algorithm. The idea is to use a standard decoupled method as a preconditioner for GMRES. In matrix form, the variational problem (5)-(7) can be written as Au + stp + CT =f, = 0, 23u + Eg=
-Du
(10)
0.
All the above matrices are independent of the unknown vector (u, p, T) except E, which is related to the convected derivative and is therefore dependent on u. The residual is then defined as
(11)
The GMRES R(u,
algorithm
applied
to the resolution
P, 7) =o
of W)
is not likely to converge because of a poor conditioning. We propose therefore to use a preconditioned residual P( u, p, T) and to solve p(u,
P, 4
= 0,
where for a given vector (u, p, T),
03)
283 with ASu+B’Sp=f-Au-B’p-C7, B 6u
= -Bu,
(15)
E, 8r
=Du-E,7.
06)
It is clear that if (6rr, Sp, 8~) is zero, the corresponding vector (u, p, 7) is the desired solution. Equations (15) and (16) can be solved separately. The GMRES algorithm is then applied to the solution of (13). Note that we now have to compute and factorize a symmetric matrix for the solution of (15) but this is done just once since that matrix does not depend on previous iterations. Moreover, this symmetric system is much smaller since it involves only the velocity degrees of freedom (pressure degrees of freedom can be eliminated by penalization or by using an Uzawa algorithm [12]). Finally, the Lesaint-Raviart method [2] allows the solution of (16) on an elementby-element basis, avoiding the construction of another large matrix. This is done by solving a small linear system on each element and eventually sweeping the domain a few times to get convergence. All of this results in an almost totally decoupled method: velocities, pressure and stress tensor (on each element) are computed separately and the GMRES method provides the recoupling of the equations. Some details concerning the practical implementation of the GMRES algorithm are worth a few comments. The most important difficulty that we have met during implementation was related to the finite difference (9) and its effective computation. Indeed, we used an Uzawa algorithm [12] to solve the linear Stokes problem (15) which consists in replacing the matrix A in (15) by A, defined by A,.=A
+ rB’B,
07)
where Y is a large penalization coefficient (Y = log). The introduction of this supplementary term makes the system (15) very poorly conditioned. The consequence of this bad conditioning was that the evaluation of the finite difference (9) suffered from important round-off errors resulting in a loss of convergence of the algorithm. The way around this is quite simple but requires a few developments. It is easily seen from the definition of the preconditioned residual P that PI and P2 are both linear with respect to u, p and 7 ( P3 is non-linear because of the matrix E,). Indeed, we can write
284 where M=A [B
1
B’ 0 ’
(19)
If we now perform the finite difference (9) in the (d,, d,, d,), we obtain for a fixed (but small) value of 12: P(u+hd,,
p+hd,,
Jd-
T+hd,) h
(T + hd,) - E,-:,,,,D(
-P(u,
P,
direction
n=
T>
u + hd,) - T + E,-‘Du h
In other words, if we denote by Jd, and Jdz the first two components Jd, they can be computed by solving the Stokes problem: A( Jd,) + B’( Jd,) B(
Jd,)
= -Ad,
= -Bd,,
of
- B’d, - Cd,, (21)
instead of using (9). The last component Jd, must, however, be computed using the finite difference (9) but this calculation is not particularly ill-conditioned, although some care must be taken in choosing the value of 12. An initial value of lo-’ was used in the reported calculations. This choice was made on the basis that the computations were performed in double precision and it is well known from elementary numerical analysis that a smaller value would perform poorly. Our last comment concerns the dimension k of the Krylov space that has to be constructed in the GMRES method. In all our tests. the value of k was between 10 and 25. Recall that this requires the construction and storage of at least k vectors of length ndf, the total number of degrees of freedom. This also means that, for one global iteration of the GMRES method, k Stokes-like problems must be solved in addition to li pseudo-solutions of the constitutive equations. It is difficult to make any comparison with the Newton method since this method is not applicable owing to the non-differentiable boundary term in eqn. (7). Finally, let us mention that the Picard method was divergent for any value of h. 4. Numerical results In order to verify the efficiency of our method. we have performed simulations for the well known stick-slip problem described in Fig. 1. The
285 (-20,l) 1
1
(2)
cw)
(3)
i (1) (4)
ww) 3
(5)
: (20,O)
(W)
(1)
(2) symmetry (3)
v = 0
(4)
u = v = 0
(5)
Fig. 1. The stick-slip
axis
v = 0
problem.
difficulty comes from the singularity at the point (0, 0) resulting from the transition from a Dirichlet to a Neumann boundary condition. This type of singularity is severe even for Newtonian flows. In all our simulations, a 600-element mesh was used, yielding roughly 20000 degrees of freedom. The mesh was concentrated in the vicinity of the singularity (see Fig. 2). This large number of degrees of freedom comes from the discontinuous Q2
Fig. 2. Finite 0.058 x 0.028).
element
mesh
near
the
singularity
(the
area
of the smallest
element
is
286
Fig. 3. Discretization
of u. p and T (p and T discontinuous).
approximation of the three components of the stress tensor. The velocity and pressure were interpolated by a nine-noded quadratic and a piecewise linear (discontinuous) approximation respectively (see Fig. 3). We also refer to Ref. 2 for a complete description and justification. We present numerical results up to De = 4.0 (the Deborah number is defined by De = A+,). Figure 4 shows the velocity profiles on the line y = 0 passing through the singularity. Figure 5 gives the 7XXprofile on the same axis. The GMRES algorithm failed to converge at a Deborah number of approximately 4.5. For Deborah numbers greater than 2, small oscillations occurred in the velocity and stress fields near the singularity. These results were obtained with (Y= 1 in eqn. (7). The oscillations were eliminated by increasing (Yup to 10. Increasing the value of (Yis equivalent, in some sense, to a modification of the constitutive equation similar to the one introduced by the streamline upwind method [3]. 0 ur numerical results are at least qualitatively similar to those of Ref. 3 although we have not been able to reach the same high values of the Deborah number. The slope of the velocity profile decreases while the maximum value of 7,., on the singularity increases when De increases. Finally, let us say a few words on the convergence behavior of GMRES. Table 1 shows the evolution of the Euclidean norm of the residual P with respect to the iteration number at De = 1. The starting solution was the
De=0.0 2.0 4.0 6-
Fig. 4. Velocity profiles on the line _Y= 0 at De = 0,2 and 4. Fig. 5. ?,A on the line y = 0 at De = 0, 2 and 4.
287 TABLE 1 Evolution of residual with iteration number Iteration 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
IIPI1 3.39 x lo2 1.76 x lo2 5.2 x10’ 1.09x10* 4.9 x10’ 2.5 x10’ 1.3 x10’ 8.0 xlOO 1.0 xlo” 9.2 x10-l 1.8 x10-l 1.1 x10-l 3.5 x1o-2 1.2 x10-* 4.8 x1O-3 1.8 x1o-3
solution at De = 0.8. 15 iterations were enough to reduce the norm of the residual by 5 orders of magnitude. Of course, convergence is slower when starting with a bad initial guess and for higher values of the Deborah number. 5. Conclusion In this paper, we have introduced a decoupled method for the simulation of viscoelastic fluid flow problems that offers convergence properties similar to a full Newton method. The main advantage of the method is that the global matrix is smaller in size since it involves only the velocity degrees of freedom. The stresses are computed on an element-by-element basis, thanks to the Lesaint-Raviart method. The preconditioned GMRES method provides a much greater reliability of the algorithm if we compare with the usual Picard- type scheme. This kind of approach is certainly of great interest, particularly in view of three-dimensional simulations. References 1 M.J. Crochet, Numerical simulation of highly viscoelastic flows, Proc. 10th Int. Congr. on Rheology, Sydney, 1988, pp. 19-24.
288 2 M. Fortin
3 4 5 6 7 8
9 10 11 12
and A. Fortin, A new approach for the FEM simulation of viscoelastic flows, J. Non-Newtonian Fluid Mech., 32 (1989) 295310. J.M. Marchal and M.J. Crochet, A new element for calculating viscoelastic flows, J. Non-Newtonian Fluid Mech., 26 (1987) 77-114. F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multiplier, RAIRO, Anal. Numer., 8 (R2) (1974) 129-151. M. Fortin and R. Pierre, On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows, Comput. Methods Appl. Mech. Eng., 73 (1989) 341-350. X.L. Luo and RI. Tanner, A streamline element scheme for solving viscoelastic flow problems, J. Non-Newtonian Fluid Mech., 22 (1986) 61-89. M. Fortin and D. Esselaoui, A finite element procedure for viscoelastic flows, Int. J. Num. Meth. Fluids, 7 (10) (1987) 1035-1052. M. El Hadj and P.A. Tanguy, A finite element procedure coupled with the method of characteristics for simulation of the viscoelastic fluid Bow. J. Non-Newtonian Fluid Mech.. 36 (1990) 333--349. R. Keunings, Simulation of viscoelastic flow, in: C.L. Tucker III (Ed.), Fundamentals of Computer Modeling for Polymer Processing, Carl Hanser, 1987. X.L. Luo and R.I. Tanner, A decoupled finite element streamline-upwind scheme for viscoelastic flow problems, J. Non Newtonian Fluid Mech., 31 (1989) 143-162. Y. Saad and M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput.. 7 (3) (1986) 856-869. M. Fortin and A. Fortin, A generalization of Uzawa’s algorithm for the solution of the Navier-Stokes equations, Commun. Appl. Numer. Methods, 1 (1985) 2055208.
Journal of Non-Newtonian Fluid Mechanics, 36 (1990) 277-288 Elsevier Science Publishers B.V., Amsterdam
A PRECONDITIONED GENERALIZED ALGORITHM FOR THE NUMERICAL OF VISCOELASTIC FLUID FLOWS
ANDRI?
MINIMAL RESIDUAL SOLUTION
FORTIN
Department of Applied Mathematics, .&ole Polytechnique Succursale A, Montrkal H3C 3A7 (Canada) MICHEL Department (Received
de MontrPal C. P. 6079,
FORTIN of Applied Mathematics,
UniversitP Lava1 Quebec GIK 7P4 (Canada)
August 24, 1989; in revised form January
10, 1990)
Abstract The numerical simulation of the flow of a viscoelastic fluid by a generalized minimal residual (GMRES) method is investigated and tested to determine its convergence properties. A method based on the LesaintRaviart technique is used as a preconditioner for the GMRES method to enhance the convergence rate and robustness of the algorithm. Numerical results are presented for the stick-slip problem. Keywords: decoupling; finite elements; GMRES; Oldroyd-B Model; viscoelasticity
1. Introduction Numerical simulations of viscoelastic fluid flows require the solution of very large strongly non-linear systems of equations. This is the main reason Ghy most successful numerical methods are based on the Newton-Raphson algorithm which offers robustness and quadratic convergence in most cases. This means that the velocities, pressure and stresses are obtained as the solution of a global system of non-linear equations that is solved by a Newton method. This requires the assembly and factorization of huge systems of equations that are computationally expensive. This is why we believe that new iterative schemes based on decoupled methods are needed and should receive more attention. In contrast, decoupled methods usually need the resolution of a larger number of smaller systems of equations. The basic idea is to solve separately 0377-0257/90/$03.50
0 1990 - Elsevier Science Publishers
B.V.
278 the constitutive equation for the stresses and the conservation equations for the velocities and pressure. Because of the strong non-linearities of viscoelastic fluid flow problems, decoupled methods have not been very successful. Convergence difficulties are quickly encountered (that is at relatively low Weissenberg numbers), precluding their use in practical situations. It is, however, generally accepted [l] that appropriate decoupled methods must be developed for the solution of industrial bidimensional applications. This statement is even stronger for three-dimensional problems. The difficulty in designing a robust decoupled method is to avoid constructing a large global matrix while keeping the good convergence properties of the Newton method. We present in this paper a generalized minimal residual (GMRES) algorithm using a decoupled method described by Fortin and Fortin [2] as a preconditioner. The resulting algorithm is proved to be efficient but of course not as robust as the full Newton method. Numerical results are presented for the stick-slip problem showing convergence of the iterative scheme up to De = 4. 2. Statement of the problem 2. I. Basic equations We briefly recall in this section the statement of the problem. An Oldroyd B model is used as a constitutive equation relating the stress tensor to the velocity field. This model has been chosen for the purpose of comparison with existing results only. The numerical solution of the flow of a viscoelastic fluid requires the computation of the velocity field u, the pressure p and the stress tensor 7 satisfying the following partial differential equations. Momentum equations: *[2q2i)(u)]
-v
-v
*7+vp=.f.
(1)
Mass conservation: V *u=o. Constitutive
(2) equation
7+x[u.v7-vu.7-T~(vu)t] where
(Oldroyd
B model): =2Tjtj+),
(3)
279
One of the major difficulties in solving the system (l)-(3) is the numerical treatment of the constitutive equation (3). The hyperbolicity of this equation requires special care as clearly demonstrated by Marchal and Crochet [3] and Fortin and Fortin [2]. Indeed, some kind of upwinding seems necessary to stabilize most numerical algorithms. In Ref. 3, a streamline upwind method is presented while the Lesaint-Raviart method is prefered in Ref. 2. This last method will be used again and we refer to Ref. 2 for a detailed description. 2.2. Discretization It is now accepted that the discretization of the system of equations (l)-(3) requires special attention. The difficulties come from the (u, p, T) formulation of the basic equations. In the Newtonian case, it is possible to discard 7 from the equations and we are left with the classical velocity-pressure formulation. Even in this simpler case, it is well known that the discretizations of the velocities and pressure must satisfy the so-called Brezzi condition (see Ref. 4). For viscoelastic fluids, it is not possible to discard r (unless h = 0) and consequently it is necessary to develop consistent (u, p, 7) discretizations. In a recent paper, Fortin and Pierre [5] have studied this particular formulation in the case of the element introduced by Marchal and Crochet [3]. Their analysis allowed them to justify the 4 by 4 subdivision of the elements for the stresses. It is possible to generalize their analysis to develop various consistent discretizations. We shall concentrate on discretizations involving discontinuous approximations for the stresses. In that case, Fortin and Pierre came to the conclusion that if V,, Q,, and IV, denote respectively the discrete spaces for the velocities, pressure and stresses, then V, and W, must satisfy the condition vv,c
w,,
(4)
i.e. the gradient of any function belonging to V, must be included in the space W,. Condition (4) is complementary to the Brezzi condition relating V, and Qh. In this paper, we use a continuous biquadratic (nine nodes) approximation for the velocities and a piecewise linear pressure (no continuity is required at element interfaces for the pressure). For the stresses, we set IV, = VV,, i.e. the space IV, consists of the gradient of all the functions belonging to V,. In other words, IV, is made up of incomplete biquadratic (eight nodes) polynomials. Since the functions in V, are only continuous (not C’) at element interfaces, the gradient of such functions will be discontinuous on
280 those interfaces. This is illustrated in Fig. 3 by setting the stress nodes inside the element. We now give the variational formulation for the system (l)-(3). Find (u, p, T) E V,?x Q,, X W, such that
=
JK
29,?(u):
cp dx dy
for every element
K and V’cpE V$.
(7)
Problem (7) can be solved on an element-by-element basis, thanks to the discontinuous approximation of the stresses and to the Lesaint-Raviart method. aK_ (aK+) denotes the inflow (outflow) part of the boundary of element K and [T] is the jump in 7 across the element boundary (see Ref. 2 for a complete justification). The parameter CYis the upwinding factor, usually taking a value between 0 and 1. However, as we shall see, a larger value can help to stabilize the numerical results. 3. Numerical method 3.1. Decoupled methods The basis of most decoupled methods is summarized by the following steps: 1. Solve the constitutive equation (7) for 7 using the velocity field obtained from a previous iteration. 2. Solve (5) and (6) for u and p using the value of 7 calculated in step 1. 3. Check for convergence and return to step 1 if necessary. Many variants of the above algorithm are possible. Most of them are based on the method of characteristics [6-81. As already mentioned, the major drawback of this method is the loss of convergence of the iterative scheme at relatively low Weissenberg numbers. For example, when using the Lesaint-Raviart method, the above algorithm failed to converge for any value of the Weissenberg number. A steady solution could, however, be obtained as a long time solution of transient simulations. For a complete
281 review of these difficulties, we refer to Keunings [9]. It is worth mentioning that appropriate treatment of the constitutive equation improves the convergence properties of decoupled methods as clearly demonstrated by Luo and Tanner [lo]. 3.2. The generalized
minimal residual method
The GMRES method was first introduced by Saad and Schultz [ll] for the resolution of non-symmetric linear systems. This method was then extended to the solution of large non-linear systems of equations. The GMRES algorithm has recently received a lot of attention. The reason for this is obvious. In its simplest form, this method requires no assembly and factorization of a global matrix. Only the residual has to be computed and stored to construct a Krylov space. The economy in memory space can be very important and depends on the dimension of that Krylov space. Let us recall the main features of the GMRES algorithm which. as we will verify, can be seen as a quasi-Newton method. We refer the reader to Ref. 11 for a more complete discussion. The basic idea is to start with the Newton algorithm for the solution of the non-linear system: R(X)=O.
(8)
The Newton algorithm is then X0 being given For nk0: 1. solve J, AX = -R(X,,) (J, 2. update Xn+l = X, + AX, 3. check for convergence.
is the Jacobian matrix at X,),
The linear system in step 1 is usually solved by LU decomposition. In the case of viscoelastic fluid flow problems, this matrix is extremely large and should be constructed and factorized at each iteration. The main feature of the GMRES algorithm is to avoid the construction of J,. To achieve this goal, two important remarks must be made. Firstly, iterative methods for solving a linear system of the form AX = b such as the conjugate gradient method (for symmetric matrices) and the minimal residual method do not require the matrix A itself but only the product of A with different direction vectors d. Moreover, in the case of a Jacobian matrix, we have lim
Jd=
n
h-0
WC,+hd) -R(X) h
(9)
282
Use of eqn. (9) eliminates the need for the explicit construction of J,. Finally, it is worth noticing that the linear system in step 1 does not need to be solved exactly. This results in a quasi-Newton method since the Jacobian matrix is only approximated by the finite difference (9) and the related linear systems are not solved exactly. We refer to Ref. 11 for a detailed description of the construction of the Krylov space for the resolution of linear systems by the GMRES method. 3.3. Application to viscoelastic jluid flow problems As already mentioned, viscoelastic fluid flow problems are strongly nonlinear and it is not likely that the GMRES method will converge, at least in its simplest form. However, a good preconditioner could possibly enhance the convergence of the algorithm. The idea is to use a standard decoupled method as a preconditioner for GMRES. In matrix form, the variational problem (5)-(7) can be written as Au + stp + CT =f, = 0, 23u + Eg=
-Du
(10)
0.
All the above matrices are independent of the unknown vector (u, p, T) except E, which is related to the convected derivative and is therefore dependent on u. The residual is then defined as
(11)
The GMRES R(u,
algorithm
applied
to the resolution
P, 7) =o
of W)
is not likely to converge because of a poor conditioning. We propose therefore to use a preconditioned residual P( u, p, T) and to solve p(u,
P, 4
= 0,
where for a given vector (u, p, T),
03)
283 with ASu+B’Sp=f-Au-B’p-C7, B 6u
= -Bu,
(15)
E, 8r
=Du-E,7.
06)
It is clear that if (6rr, Sp, 8~) is zero, the corresponding vector (u, p, 7) is the desired solution. Equations (15) and (16) can be solved separately. The GMRES algorithm is then applied to the solution of (13). Note that we now have to compute and factorize a symmetric matrix for the solution of (15) but this is done just once since that matrix does not depend on previous iterations. Moreover, this symmetric system is much smaller since it involves only the velocity degrees of freedom (pressure degrees of freedom can be eliminated by penalization or by using an Uzawa algorithm [12]). Finally, the Lesaint-Raviart method [2] allows the solution of (16) on an elementby-element basis, avoiding the construction of another large matrix. This is done by solving a small linear system on each element and eventually sweeping the domain a few times to get convergence. All of this results in an almost totally decoupled method: velocities, pressure and stress tensor (on each element) are computed separately and the GMRES method provides the recoupling of the equations. Some details concerning the practical implementation of the GMRES algorithm are worth a few comments. The most important difficulty that we have met during implementation was related to the finite difference (9) and its effective computation. Indeed, we used an Uzawa algorithm [12] to solve the linear Stokes problem (15) which consists in replacing the matrix A in (15) by A, defined by A,.=A
+ rB’B,
07)
where Y is a large penalization coefficient (Y = log). The introduction of this supplementary term makes the system (15) very poorly conditioned. The consequence of this bad conditioning was that the evaluation of the finite difference (9) suffered from important round-off errors resulting in a loss of convergence of the algorithm. The way around this is quite simple but requires a few developments. It is easily seen from the definition of the preconditioned residual P that PI and P2 are both linear with respect to u, p and 7 ( P3 is non-linear because of the matrix E,). Indeed, we can write
284 where M=A [B
1
B’ 0 ’
(19)
If we now perform the finite difference (9) in the (d,, d,, d,), we obtain for a fixed (but small) value of 12: P(u+hd,,
p+hd,,
Jd-
T+hd,) h
(T + hd,) - E,-:,,,,D(
-P(u,
P,
direction
n=
T>
u + hd,) - T + E,-‘Du h
In other words, if we denote by Jd, and Jdz the first two components Jd, they can be computed by solving the Stokes problem: A( Jd,) + B’( Jd,) B(
Jd,)
= -Ad,
= -Bd,,
of
- B’d, - Cd,, (21)
instead of using (9). The last component Jd, must, however, be computed using the finite difference (9) but this calculation is not particularly ill-conditioned, although some care must be taken in choosing the value of 12. An initial value of lo-’ was used in the reported calculations. This choice was made on the basis that the computations were performed in double precision and it is well known from elementary numerical analysis that a smaller value would perform poorly. Our last comment concerns the dimension k of the Krylov space that has to be constructed in the GMRES method. In all our tests. the value of k was between 10 and 25. Recall that this requires the construction and storage of at least k vectors of length ndf, the total number of degrees of freedom. This also means that, for one global iteration of the GMRES method, k Stokes-like problems must be solved in addition to li pseudo-solutions of the constitutive equations. It is difficult to make any comparison with the Newton method since this method is not applicable owing to the non-differentiable boundary term in eqn. (7). Finally, let us mention that the Picard method was divergent for any value of h. 4. Numerical results In order to verify the efficiency of our method. we have performed simulations for the well known stick-slip problem described in Fig. 1. The
285 (-20,l) 1
1
(2)
cw)
(3)
i (1) (4)
ww) 3
(5)
: (20,O)
(W)
(1)
(2) symmetry (3)
v = 0
(4)
u = v = 0
(5)
Fig. 1. The stick-slip
axis
v = 0
problem.
difficulty comes from the singularity at the point (0, 0) resulting from the transition from a Dirichlet to a Neumann boundary condition. This type of singularity is severe even for Newtonian flows. In all our simulations, a 600-element mesh was used, yielding roughly 20000 degrees of freedom. The mesh was concentrated in the vicinity of the singularity (see Fig. 2). This large number of degrees of freedom comes from the discontinuous Q2
Fig. 2. Finite 0.058 x 0.028).
element
mesh
near
the
singularity
(the
area
of the smallest
element
is
286
Fig. 3. Discretization
of u. p and T (p and T discontinuous).
approximation of the three components of the stress tensor. The velocity and pressure were interpolated by a nine-noded quadratic and a piecewise linear (discontinuous) approximation respectively (see Fig. 3). We also refer to Ref. 2 for a complete description and justification. We present numerical results up to De = 4.0 (the Deborah number is defined by De = A+,). Figure 4 shows the velocity profiles on the line y = 0 passing through the singularity. Figure 5 gives the 7XXprofile on the same axis. The GMRES algorithm failed to converge at a Deborah number of approximately 4.5. For Deborah numbers greater than 2, small oscillations occurred in the velocity and stress fields near the singularity. These results were obtained with (Y= 1 in eqn. (7). The oscillations were eliminated by increasing (Yup to 10. Increasing the value of (Yis equivalent, in some sense, to a modification of the constitutive equation similar to the one introduced by the streamline upwind method [3]. 0 ur numerical results are at least qualitatively similar to those of Ref. 3 although we have not been able to reach the same high values of the Deborah number. The slope of the velocity profile decreases while the maximum value of 7,., on the singularity increases when De increases. Finally, let us say a few words on the convergence behavior of GMRES. Table 1 shows the evolution of the Euclidean norm of the residual P with respect to the iteration number at De = 1. The starting solution was the
De=0.0 2.0 4.0 6-
Fig. 4. Velocity profiles on the line _Y= 0 at De = 0,2 and 4. Fig. 5. ?,A on the line y = 0 at De = 0, 2 and 4.
287 TABLE 1 Evolution of residual with iteration number Iteration 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
IIPI1 3.39 x lo2 1.76 x lo2 5.2 x10’ 1.09x10* 4.9 x10’ 2.5 x10’ 1.3 x10’ 8.0 xlOO 1.0 xlo” 9.2 x10-l 1.8 x10-l 1.1 x10-l 3.5 x1o-2 1.2 x10-* 4.8 x1O-3 1.8 x1o-3
solution at De = 0.8. 15 iterations were enough to reduce the norm of the residual by 5 orders of magnitude. Of course, convergence is slower when starting with a bad initial guess and for higher values of the Deborah number. 5. Conclusion In this paper, we have introduced a decoupled method for the simulation of viscoelastic fluid flow problems that offers convergence properties similar to a full Newton method. The main advantage of the method is that the global matrix is smaller in size since it involves only the velocity degrees of freedom. The stresses are computed on an element-by-element basis, thanks to the Lesaint-Raviart method. The preconditioned GMRES method provides a much greater reliability of the algorithm if we compare with the usual Picard- type scheme. This kind of approach is certainly of great interest, particularly in view of three-dimensional simulations. References 1 M.J. Crochet, Numerical simulation of highly viscoelastic flows, Proc. 10th Int. Congr. on Rheology, Sydney, 1988, pp. 19-24.
288 2 M. Fortin
3 4 5 6 7 8
9 10 11 12
and A. Fortin, A new approach for the FEM simulation of viscoelastic flows, J. Non-Newtonian Fluid Mech., 32 (1989) 295310. J.M. Marchal and M.J. Crochet, A new element for calculating viscoelastic flows, J. Non-Newtonian Fluid Mech., 26 (1987) 77-114. F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multiplier, RAIRO, Anal. Numer., 8 (R2) (1974) 129-151. M. Fortin and R. Pierre, On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows, Comput. Methods Appl. Mech. Eng., 73 (1989) 341-350. X.L. Luo and RI. Tanner, A streamline element scheme for solving viscoelastic flow problems, J. Non-Newtonian Fluid Mech., 22 (1986) 61-89. M. Fortin and D. Esselaoui, A finite element procedure for viscoelastic flows, Int. J. Num. Meth. Fluids, 7 (10) (1987) 1035-1052. M. El Hadj and P.A. Tanguy, A finite element procedure coupled with the method of characteristics for simulation of the viscoelastic fluid Bow. J. Non-Newtonian Fluid Mech.. 36 (1990) 333--349. R. Keunings, Simulation of viscoelastic flow, in: C.L. Tucker III (Ed.), Fundamentals of Computer Modeling for Polymer Processing, Carl Hanser, 1987. X.L. Luo and R.I. Tanner, A decoupled finite element streamline-upwind scheme for viscoelastic flow problems, J. Non Newtonian Fluid Mech., 31 (1989) 143-162. Y. Saad and M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput.. 7 (3) (1986) 856-869. M. Fortin and A. Fortin, A generalization of Uzawa’s algorithm for the solution of the Navier-Stokes equations, Commun. Appl. Numer. Methods, 1 (1985) 2055208.