- Email: [email protected]
A finite element computation of moderate Reynolds fluid flow using a modified marquardt method
COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 70 (1988) 139-149 NORTH-HOLLAND
A FINITE ELEMENT COMPUTATION OF MODERATE REYNOLDS FLUID FLOW USING A MODIFIED MARQUARDT METHOD D. FAUCHON and P.A. TANGUY DEpartment de genie chimique, UniversitE Laval, Quebec G1K 71'4, Canada
and R.E. HAYES Department of Chemical Engineering, University of Alberta, Edmonton, Alberta T6G 2G6, Canada Received 1 June 1987 Revised manuscript received 28 January 1988
A modification of the Marquardt method is combined with a Galerkin finite element method to produce an efficient algorithm for the solution of the steady-state Navier-Stokes equation at moderate Reynolds numbers. The algorithm is shown to be 3 to 4 times faster than more conventional algorithms based on the modified Newton method with loading on the nonlinear term. Numerical results are presented for a 4:1 contraction at Reynolds numbers in the range 400-1000 and for problems varying in size from 1107 to 3576 equations.
1. Introduction The development of efficient algorithms for the numerical simulation of incompressible fluid flow problems in a wide range of Reynolds number is currently an area of active research. Most of the problems of practical interest are nonlinear, which can lead to difficulties in obtaining the solution of the governing equations at a reasonable cost. Indeed, the classical algorithms may fail to converge when a poor initial estimate is used to start the iterative procedure, or when the problem is very stiff, that is, when steep gradients are present. At the numerical level, the governing differential equations are discretized using (usually) a finite difference or a finite element method in order to convert the partial differential equations into a series of nonlinear algebraic equations. This set of equations must be solved iteratively using a suitable algorithm and an initial guess. Several methods with varying degrees of complexity (first-order and second-order as well as hybrids) have been used, namely: Picard method [1], Newton method [2], modified Newton and quasi-Newton methods [3], with or without Broyden update [4]. Since we are facing an initial value problem, the initial guess is of paramount importance as for the existence and the convergence rate of the numerical solution. These methods can be viewed as various forms of preconditioning. From a practical point of view, the solution of a stiff nonlinear prob!em is generally 0045-7825188/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
140
D. Fauchon et al., A finite element computation of moderate Reynolds
achieved by using some form of incremental loading strategy, either on the boundary conditions or on the nonlinearity itself. This technique is not easy to optimize, hence it can be very costly to implement. In a recent publication [5], a self-adaptive algorithm based on the Marquardt method was introduced for the solution of nonlinear problems. We recall here that the Mazquardt method is a standard technique for parameter estimation [6] and is utilized to seek adequate initial values for the parameters in the iterative process. This new methodology has proven much superior to a loading strategy for obtaining the finite difference solution of a classical stiff chemical engineering problem (coupled heat and mass transfer with exothermic chemical reaction). In this work, we extended this idea to develop an algorithm for the finite element solution of the Navier-Stokes equation at moderate Reynolds numbers. It consists of a combination of a modified Marquardt method and a modified Newton method which exhibits a self-optimizing feature and enables cost-effective solution of very stiff problems. In the following sections, we will outline the nature of the problem to be solved, describe the algorithm developed, and then present some numerical results obtained with it.
2. Problem statement Let us consider the numerical solution of the two-dimensional 4:1 sudden contraction (Fig. 1). This problem, which is a modellers, provides a good test for a solution algorithm because pattern, which is strongly dependent on the Reynolds number. written, at steady state, as
Navier-Stokes equations in a ~avorite one for numerical of the complexity of the flow The governing equations are
p(u .Vu)= -Vp + V2u,
(1)
No slip conditions
J
,rl
S
0 I.J 0. c.) 0 cO
4..I a) r..l' 1=4
Centrellne,
symmetry conditions
Fig. 1. Solution domain and relevant boundary conditions.
Li
D. Fauchon et al., A finite element computation of moderate Reynolds
V.u=O,
141
(2)
where u is the velocity vector, p the pressure, p the fluid density, and tz its viscosity. The boundary conditions are: zero slip condition at the walls, a parabolic velocity profiJe at the inlet, a free outlet, and symmetry conditions at the centerline (Fig. 1). Equations (1) and (2) are put into a weak variational form using a standard Galerkin finite element method [7]. The weak variational form is discretized using the triangular Crouzeix-Raviart element [8], yielding a matrix problem that we solve using a Uzawa algorithm [9]. At the kth iteration, the problem is written as K,(uk)u
~ = F + Btp k-l ,
K , ( U k) = K ( U k) + ( l l r ) B t B
(3) ,
pk = pt-t + (I/r)BU k '
(4) (S)
where K is the Navier-Stokes matrix representing the diffusion and advection contributions, K, the penalty Navier-Stokes matrix, B(B t) the discrete (transpose) divergence matrix, U the global velocity vector, P the global pressure vector, F the boundary condition vector, and r the penalty parameter of order 10 -8. At Reynolds numbers higher than unity, the problem becomes advection driven and an appropriate method must be selected to get the solution. Let us consider for instance the standard Newton method for the solution of this problem. At the kth iteration, (3) is transformed into j k A U ~ = R k = F + B t p k-I - K , ( U k - I ) U k - t ,
(6)
where J is the Jacobian matrix, d K , ( U i ) / d U j, A U is the increment on the unknown vector U, and R is the residual. Provided that a good initial estimate of U is available, the problem will converge rapidly such that A U becomes vanishingly small. If the initial guess is outside of the radius of convergence, however, the problem will diverge rapidly. This behaviour occurs in finite element as well as in finite difference and is related to the eigenvalues of the matrix.
3. Modified Marqmrdt method We have tackled a similar problem with the finite difference method [5] by adapting the Marquardt method such that the initial estimates are within the radius of convergence of a modified Jacobian matrix. Briefly, the method used is as follows. The problem, formulated in terms of a standard Newton-type Jacobian matrix, H, and a standard residual vector, S, is no
= s
.
(7)
A modified problem is then defined as
n,o=s,
(8)
D. Fauchon et al., A finite element computation of moderate Reynolds
142
where HA = H + AI. A is a scaling factor and I the identity matrix. This modified problem is solved iteratively. The values are then updated and the sum of the squares of the components of the residual vector S, O, is computed. If • is reduced, then the new O's are accepted, the value of A is reduced by a factor v, and the modified problem is solved again. If not, then A is increased by the factor v and then the modified problem is resolved. This A-dependent scheme is continued until a suitable convergence criterion is satisfied. When convergence of the solution is attained, the value of • is essentially equal to zero and the value of A should become vanishingly small. Complete details are given in [5]. In the present work, the same methodology is applied in a finite element context. The major point that arises is the determination of an appropriate variational equivalent of the identity matrix. Let us consider the time-dependent Navier-Stokes equation p( OU/ Ot "~" U" VU) • --Vp 4" ~tV2U ,
V.u=O,
(9) (10)
and its variational equivalent
{[(1/At)M + K,(Uk)]Uk}'+~'= {F + B t p k - 1 ) t+At + (1/At)MU' .
(11)
It is apparent that one could use the mass matrix M, which represents the time derivatiive term in the unsteady-state problem. In this case A would be the inverse of the time step. Problem (9) would then correspond to an unsteady-state problem with a variable time step which is automatically scaled so that it maximizes the rate of convergence, yielding a subtle timestepping procedure. Instead of using a time-stepping technique, which is in general expensive since small time steps must be used to achieve convergence (simulation of a start-up), we want to consider directly the steady-state problem (7). To do this, we replace the mass matrix on the left-hand side of (11) with the identity matrix and suppress the contribution of the mass-matrix terms on the right-hand side. Implementation of the Marquardt method in this fashion can now be viewed more as a matrix preconditioning technique whose numerical efficiency will be directly related to the updating strategy on A.
4. Methodology The two algorithms compared in this work are: (a) a straight modified Newton method with an optimized loading strategy, and (b) a modified Newton method incorporating a Marquardtlike method, subsequently referred to as the Newton-Marquardt method. Among the classical algorithms, the modified Newton method [3] was found to be the most efficient solution technique for this type of problem
4.1. Modified Newton algorithm Starting from a zero velocity field, and heuristically setting the viscosity to an optimal value (as small as possible but still yielding convergence):
D. Fauchon et al., A finite element computation of moderate Reynolds
143
Algorithm (a) (modified Newton) Step Step Step Step Step
1. 2. 3. 4. 5.
Compute the Jacobian matrix J (equation (6)). Compute the residual R (equation (6)). Solve the matrix problem for A U. Test for convergence and return to Step 2 if necessary. Update the viscosity and go to Step 1.
Steps 1 to 5 are repeated until a solution is obtained at the required Reynolds number.
4.?. Newton-Marquardt algorithm Starting from either a zero or a nonzero velocity field, and the Reynolds number at which the final solution is desired, A and v being given:
Algorithm (b) (Newton-Marquardt) Compute the sum of the squares of the residuals, call this value • 0. Compute a modified Jacobian (equation (8)). Compute the residual R (equation (6)). Solve the matrix problem for A U. Update the velocity. Compute a sum of the squares of the new residuals, call this value O~. Test for O1 < Oo: if no, then A = A. p and goto Step 2; if yes, goto Step 8. Step 8. Test for (Oo - O1)/Ol < tolerance: if no, ~0 = ~ , A = A/v, and goto Step 2; if yes, A - 0 , goto Algorithm (a) (modified Newton) to finish problem.
Step 1. Step 2. Step 3. Step 4. Step 5. Step 6. Step 7.
The first point which must be considered is the determination of the parameter values A and v, the selection of which is the crux of the method. We have found that the constant A should be selected according to the magnitude of the diagonal elements, with an initial guess of the order of the smallest of these elements being reasonable. It should be recalled that the magnitude of these diagonal elements are a function of the mesh size and the value of the penalty parameter. Such a choice for A arises from stability considerations. Indeed, a value of A based for instance on larger elements would strongly perturbate the diagonal of the elementary nlatrices of the smallest elements, which would ~eld erroneous results. The optimal value of v was found heuristically to lie in the range from 2-5, which is consistent with [6]. Smaller values increased the execution time while larger values caused divergence. This choice for A and ~, is also consistent with the findings for the finite difference method [5]. In the original work by Ma~quardt [6], the stopping criterion was based on the value of O, the sum of the squares of the residual components. In this work, however, we observed that the value of • behaved in an erratic manner, first increasing, then reaching a plateau before decreasing. Oscillations of the value of • were also observed to occur. This behaviour is shown in Fig. 2 and will be explained in more detail in Section 5. We therefore decided to use the change in the value of • between iterations as the stopping criterion, that is, when the rate
144
D. Fauchon et al., A finite element computation of moderate Reynolds 0.31
0.21
"9-
J
0.11
0.01
0
. . . .
I
2
3
4
5
"
6
Number of iterations
Fig. 2. Value of the sum of the square of the residuals as a function of iteration number. Reynolds number 1000, 2303 equations.
of change of • had fallen below a heuristically determined tolerance, the value of A was set to zero. For a A of zero, this algorithm is then mathematically equivalent to using the modified Newton method, hence we used Algorithm (a) for the final iterations, i.e. until the problem converges to the correct solution. This is done for two reasons. The first is to save computation time, for when the rate of change in • has fallen below a certain tolerance, the current estimate of U is sufficiently accurate to achieve convergence of the standard steadystate problem using a modified Newton method. Since this is true, it is pointless as well as time consuming to compute the modified Jacobian and the value of • at each iteration. The second problem which can arise is the previously mentioned oscillatory behaviour about • which was observed to occur when the rate of change of • was small. This may be caused by local instability and can lead to "thrashing" of the algorithm with non convergence. This behaviour which was not studied in depth will be addressed in the near future.
$. Numerical results The preceeding problem was solved on a Control Data Cyber 180-830 computer for Reynolds numbers between 400 and 1000. The number of equations varied from 1107 to 3576.
D. Fauchon a al., A finite element computation of moderate Reynolds
145
Initial estimates of a zero-velocity field and the velocity field resulting from solution of the Stokes problem at Re -0.001 were both tested. The comparison of the two algorithms tested is based on the execution time required to reach the solution. The change of the residual • with number of iterations is illustrated in Fig. 2. It can be seen that, at the beginning of the problem, • is essentially zero since the only contribution to comes from the boundary conditions (U and P were initially taken equal to zero). As the iteration process continues, the algorithm adjusts automatically the value of ,~ so as to search for a good estimate of U and P (Step 7) of Algorithm (b)). This stage corresponds to an increase of O. Once a satisfactorj g,Jess is found, O decreases but is subjected to oscillations while the problem is converging. The large vana~.~.."~':~-of O as well as the observed potential for oscillatory behaviour led to the conclusion that using O as a stopping criterion was not an optimal way of proceeding. Instead we used the relative value of • between iterations, i.e. ~!~ The convergence behaviour of the Newton-Marquardt method is observed to be strongly influenced by the number of equations to be solved, the Reynolds number at which the solution is desired, and the initial guess. Figures 3-5 show the change in value of the residual versus iterations for various Reynolds numbers and the number of equations. In each case the nonzero-velocity initial guess corresponds to a velocity field at Reynolds number of 0.001, which is essentially a Stokes problem. The value of • ~ / O is seen to be smaller in the case of a nonzero-velocity guess at each Reynolds number. The implication of this is t ~ t the maximum e
1.80 1.60 1.40 1.20 %~" 1.00
~' 0.80 0.60
° ~O/oo0"
0.40 0.20 0
I
2
3
4
5
Number of iterotions Fig. 3. Relative residual as a function of the iteration number. Reynolds number 400, 3576 equations.
146
D. Fauchon et al., A finite element computation of moderate Reynolds
s
I
4 b
2
!
I
2
5
4
Humber of iterations Fig. 4. Relative residual as a function of the iteration number. Reynolds number 1000, 1303 equations.
Reynolds number at which the Newton-Marquardt method is effective, is increased for any given number of equations. This is shown in Fig. 6, which is a plot of the maximum Reynolds number at which a solution could be obtained with the Newton-Marquardt method versus the number of equations. We observe that the nonzero-velocity guess leads to superior performance. The price to pay for this improved robustness is the cost of solving the Stokes problem to calculate the initial velocity field. The decrease in the maximum Reynolds number with increasing number of equations can be related to the influence of the A1 matrix on the Jacobian matrix. As the number of equations is increased, the matrix bandwidth is enlarged, therefore the influence of the main diagonal is diminished, thus reducing the effectiveness of this method. The performance of the two algorithms are compared directly in Figs. 7 and 8. In Fig. 7 curves of the execution time versus the number of equations are given for a Reynolds number of 400. The Newton-Marquardt method is consistently better, decreasing the execution time by approximately 75%. The rather unusual result of a decrease in execution time for an increase in the number of equations observed at the beginning of the modified Newton curve results from the tolerance chosen to satisfy convergence. At 1107 equations one additional iteration was required, because AU/U was slightly higher than the tolerance and the final solution has, in fact, a much lower error. Figure 8 shows the execution time versus the Reynolds number for 1303 equations. We observe from this graph that at low Reynolds number (<2) the modified Newton method is superior to the Newton-Marquardt, but as the Reynolds number increases the difference between the two algorithms increases until, at a
1.80
e60
"' '
1.40
I
1.20
~) 1.00 ~'~" 0.80 0.60 0.400.20
0
I
'"
2
3
4
5
Numberof iterations
6
Fig. 5. Relative residual as a function of the iteration number. Reynolds number 1000, 2303 equations.
1200
,
m
....
I guess ,.
•- 800 E w
c o
~. 600
4OO
i
20~)1000
1500
2(XX)
2500 3000 Numberof equations
3500
Fig. 6. Maximum Reynolds number possible as a function of the number of equations.
147
4()00
148
D. Fauchon et al., A finite element computation of moderate Reynolds
- Newton ,o
f ! This work
i 103 U
ioZl
I
IOOO
1500
2O00
25OO
I
3000
3500
4000
Number of equations
Fig. 7. Execution time as a function of the number of equations. Reynolds number 400. 1600
1300
I000
[ u
700
4OO
I00
to-I
j
~ ~ ~,.
~-~'f
I00
......
l
i01 Reynolds number
i02
Fig. 8. Execution time as a function of the Reynolds number, 1303 equations.
I0 3
D. Fauchon
et
al., A finite element computation of moderate Reynolds
149
Reynolds number of 1000, the Newton-Marquardt method gets approximately 3 times faster. The trend is, moreover, clearly in favour of the Newton-Marquardt method, as the slope of the curve is not increasing at so fast a rate.
6. Conclusions and recommendations The Newton-Marquardt method is observed to provide a robust algorithm for the solution of the Navier-Stokes equations by the finite element method up to a Reynolds number of 1000. This algorithm was shown to be more efficient than the classical modified Newton method and can result in considerable cost savings. Furthermore, this algorithm exhibits a degree of self-opt'mization which makes it ideal for implementation. Further work may be desirable in order to expand the range of Reynolds numbers for which the method is effective.
Acknowledgment The financial assistance of Natural Science and Engineering Research Council of Canada and Control Data Corporation is gratefully acknowledged.
References [1] [2] [3] [4]
S.H. Crandail, Engineering Analysis (McGraw-Hill, New York, 1956). R.L. Burden and J.D. Faires, Numerical Analysis (Prindle, Weber and Smidt, Boston, 3rd. ed., 1985). E.J. Dennis and J. Mor~, Quasi-Newton methods, motivation and theory, SIAM Rev. 19 (19/7) 46-89. H. Matthies and G. Strang, The solution of nonlinear finite element equations, lnternat. J. Numer. Meths. Engrg. 14 (1979) 1613-1626. [5] R.E. Hayes and EA. Tanguy, Solution of convection-diffusion differential equations using Marquardt's method, Canad. J. Chem. Engrg. 65 (1987) 522-525. [6] D.W. Marquardt, Estimation of non-linear parameters, J. Soc. Ind. Appl. Math 11 (1963) 431-441. [7] K. Bathe, Finite Element Procedures in Engineering Analysis (Prentice-Hall, Engiewood Cliffs, NJ, 1982). [8] M. Crouzeix and P.A. Raviart, Conforming and non-conforming finite element methods for solving the stationary Stokes equations, RAIRO Anal. Numer. 7 (1973) 33. [9] EA. Tanguy, M. Fortin and L. Choplin, Finite element simulation of dip coating, II: Non-Newtonian fluids, Intemat. J. Numer. Meths. Fluids 4 (1984) 459--475.
A FINITE ELEMENT COMPUTATION OF MODERATE REYNOLDS FLUID FLOW USING A MODIFIED MARQUARDT METHOD D. FAUCHON and P.A. TANGUY DEpartment de genie chimique, UniversitE Laval, Quebec G1K 71'4, Canada
and R.E. HAYES Department of Chemical Engineering, University of Alberta, Edmonton, Alberta T6G 2G6, Canada Received 1 June 1987 Revised manuscript received 28 January 1988
A modification of the Marquardt method is combined with a Galerkin finite element method to produce an efficient algorithm for the solution of the steady-state Navier-Stokes equation at moderate Reynolds numbers. The algorithm is shown to be 3 to 4 times faster than more conventional algorithms based on the modified Newton method with loading on the nonlinear term. Numerical results are presented for a 4:1 contraction at Reynolds numbers in the range 400-1000 and for problems varying in size from 1107 to 3576 equations.
1. Introduction The development of efficient algorithms for the numerical simulation of incompressible fluid flow problems in a wide range of Reynolds number is currently an area of active research. Most of the problems of practical interest are nonlinear, which can lead to difficulties in obtaining the solution of the governing equations at a reasonable cost. Indeed, the classical algorithms may fail to converge when a poor initial estimate is used to start the iterative procedure, or when the problem is very stiff, that is, when steep gradients are present. At the numerical level, the governing differential equations are discretized using (usually) a finite difference or a finite element method in order to convert the partial differential equations into a series of nonlinear algebraic equations. This set of equations must be solved iteratively using a suitable algorithm and an initial guess. Several methods with varying degrees of complexity (first-order and second-order as well as hybrids) have been used, namely: Picard method [1], Newton method [2], modified Newton and quasi-Newton methods [3], with or without Broyden update [4]. Since we are facing an initial value problem, the initial guess is of paramount importance as for the existence and the convergence rate of the numerical solution. These methods can be viewed as various forms of preconditioning. From a practical point of view, the solution of a stiff nonlinear prob!em is generally 0045-7825188/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
140
D. Fauchon et al., A finite element computation of moderate Reynolds
achieved by using some form of incremental loading strategy, either on the boundary conditions or on the nonlinearity itself. This technique is not easy to optimize, hence it can be very costly to implement. In a recent publication [5], a self-adaptive algorithm based on the Marquardt method was introduced for the solution of nonlinear problems. We recall here that the Mazquardt method is a standard technique for parameter estimation [6] and is utilized to seek adequate initial values for the parameters in the iterative process. This new methodology has proven much superior to a loading strategy for obtaining the finite difference solution of a classical stiff chemical engineering problem (coupled heat and mass transfer with exothermic chemical reaction). In this work, we extended this idea to develop an algorithm for the finite element solution of the Navier-Stokes equation at moderate Reynolds numbers. It consists of a combination of a modified Marquardt method and a modified Newton method which exhibits a self-optimizing feature and enables cost-effective solution of very stiff problems. In the following sections, we will outline the nature of the problem to be solved, describe the algorithm developed, and then present some numerical results obtained with it.
2. Problem statement Let us consider the numerical solution of the two-dimensional 4:1 sudden contraction (Fig. 1). This problem, which is a modellers, provides a good test for a solution algorithm because pattern, which is strongly dependent on the Reynolds number. written, at steady state, as
Navier-Stokes equations in a ~avorite one for numerical of the complexity of the flow The governing equations are
p(u .Vu)= -Vp + V2u,
(1)
No slip conditions
J
,rl
S
0 I.J 0. c.) 0 cO
4..I a) r..l' 1=4
Centrellne,
symmetry conditions
Fig. 1. Solution domain and relevant boundary conditions.
Li
D. Fauchon et al., A finite element computation of moderate Reynolds
V.u=O,
141
(2)
where u is the velocity vector, p the pressure, p the fluid density, and tz its viscosity. The boundary conditions are: zero slip condition at the walls, a parabolic velocity profiJe at the inlet, a free outlet, and symmetry conditions at the centerline (Fig. 1). Equations (1) and (2) are put into a weak variational form using a standard Galerkin finite element method [7]. The weak variational form is discretized using the triangular Crouzeix-Raviart element [8], yielding a matrix problem that we solve using a Uzawa algorithm [9]. At the kth iteration, the problem is written as K,(uk)u
~ = F + Btp k-l ,
K , ( U k) = K ( U k) + ( l l r ) B t B
(3) ,
pk = pt-t + (I/r)BU k '
(4) (S)
where K is the Navier-Stokes matrix representing the diffusion and advection contributions, K, the penalty Navier-Stokes matrix, B(B t) the discrete (transpose) divergence matrix, U the global velocity vector, P the global pressure vector, F the boundary condition vector, and r the penalty parameter of order 10 -8. At Reynolds numbers higher than unity, the problem becomes advection driven and an appropriate method must be selected to get the solution. Let us consider for instance the standard Newton method for the solution of this problem. At the kth iteration, (3) is transformed into j k A U ~ = R k = F + B t p k-I - K , ( U k - I ) U k - t ,
(6)
where J is the Jacobian matrix, d K , ( U i ) / d U j, A U is the increment on the unknown vector U, and R is the residual. Provided that a good initial estimate of U is available, the problem will converge rapidly such that A U becomes vanishingly small. If the initial guess is outside of the radius of convergence, however, the problem will diverge rapidly. This behaviour occurs in finite element as well as in finite difference and is related to the eigenvalues of the matrix.
3. Modified Marqmrdt method We have tackled a similar problem with the finite difference method [5] by adapting the Marquardt method such that the initial estimates are within the radius of convergence of a modified Jacobian matrix. Briefly, the method used is as follows. The problem, formulated in terms of a standard Newton-type Jacobian matrix, H, and a standard residual vector, S, is no
= s
.
(7)
A modified problem is then defined as
n,o=s,
(8)
D. Fauchon et al., A finite element computation of moderate Reynolds
142
where HA = H + AI. A is a scaling factor and I the identity matrix. This modified problem is solved iteratively. The values are then updated and the sum of the squares of the components of the residual vector S, O, is computed. If • is reduced, then the new O's are accepted, the value of A is reduced by a factor v, and the modified problem is solved again. If not, then A is increased by the factor v and then the modified problem is resolved. This A-dependent scheme is continued until a suitable convergence criterion is satisfied. When convergence of the solution is attained, the value of • is essentially equal to zero and the value of A should become vanishingly small. Complete details are given in [5]. In the present work, the same methodology is applied in a finite element context. The major point that arises is the determination of an appropriate variational equivalent of the identity matrix. Let us consider the time-dependent Navier-Stokes equation p( OU/ Ot "~" U" VU) • --Vp 4" ~tV2U ,
V.u=O,
(9) (10)
and its variational equivalent
{[(1/At)M + K,(Uk)]Uk}'+~'= {F + B t p k - 1 ) t+At + (1/At)MU' .
(11)
It is apparent that one could use the mass matrix M, which represents the time derivatiive term in the unsteady-state problem. In this case A would be the inverse of the time step. Problem (9) would then correspond to an unsteady-state problem with a variable time step which is automatically scaled so that it maximizes the rate of convergence, yielding a subtle timestepping procedure. Instead of using a time-stepping technique, which is in general expensive since small time steps must be used to achieve convergence (simulation of a start-up), we want to consider directly the steady-state problem (7). To do this, we replace the mass matrix on the left-hand side of (11) with the identity matrix and suppress the contribution of the mass-matrix terms on the right-hand side. Implementation of the Marquardt method in this fashion can now be viewed more as a matrix preconditioning technique whose numerical efficiency will be directly related to the updating strategy on A.
4. Methodology The two algorithms compared in this work are: (a) a straight modified Newton method with an optimized loading strategy, and (b) a modified Newton method incorporating a Marquardtlike method, subsequently referred to as the Newton-Marquardt method. Among the classical algorithms, the modified Newton method [3] was found to be the most efficient solution technique for this type of problem
4.1. Modified Newton algorithm Starting from a zero velocity field, and heuristically setting the viscosity to an optimal value (as small as possible but still yielding convergence):
D. Fauchon et al., A finite element computation of moderate Reynolds
143
Algorithm (a) (modified Newton) Step Step Step Step Step
1. 2. 3. 4. 5.
Compute the Jacobian matrix J (equation (6)). Compute the residual R (equation (6)). Solve the matrix problem for A U. Test for convergence and return to Step 2 if necessary. Update the viscosity and go to Step 1.
Steps 1 to 5 are repeated until a solution is obtained at the required Reynolds number.
4.?. Newton-Marquardt algorithm Starting from either a zero or a nonzero velocity field, and the Reynolds number at which the final solution is desired, A and v being given:
Algorithm (b) (Newton-Marquardt) Compute the sum of the squares of the residuals, call this value • 0. Compute a modified Jacobian (equation (8)). Compute the residual R (equation (6)). Solve the matrix problem for A U. Update the velocity. Compute a sum of the squares of the new residuals, call this value O~. Test for O1 < Oo: if no, then A = A. p and goto Step 2; if yes, goto Step 8. Step 8. Test for (Oo - O1)/Ol < tolerance: if no, ~0 = ~ , A = A/v, and goto Step 2; if yes, A - 0 , goto Algorithm (a) (modified Newton) to finish problem.
Step 1. Step 2. Step 3. Step 4. Step 5. Step 6. Step 7.
The first point which must be considered is the determination of the parameter values A and v, the selection of which is the crux of the method. We have found that the constant A should be selected according to the magnitude of the diagonal elements, with an initial guess of the order of the smallest of these elements being reasonable. It should be recalled that the magnitude of these diagonal elements are a function of the mesh size and the value of the penalty parameter. Such a choice for A arises from stability considerations. Indeed, a value of A based for instance on larger elements would strongly perturbate the diagonal of the elementary nlatrices of the smallest elements, which would ~eld erroneous results. The optimal value of v was found heuristically to lie in the range from 2-5, which is consistent with [6]. Smaller values increased the execution time while larger values caused divergence. This choice for A and ~, is also consistent with the findings for the finite difference method [5]. In the original work by Ma~quardt [6], the stopping criterion was based on the value of O, the sum of the squares of the residual components. In this work, however, we observed that the value of • behaved in an erratic manner, first increasing, then reaching a plateau before decreasing. Oscillations of the value of • were also observed to occur. This behaviour is shown in Fig. 2 and will be explained in more detail in Section 5. We therefore decided to use the change in the value of • between iterations as the stopping criterion, that is, when the rate
144
D. Fauchon et al., A finite element computation of moderate Reynolds 0.31
0.21
"9-
J
0.11
0.01
0
. . . .
I
2
3
4
5
"
6
Number of iterations
Fig. 2. Value of the sum of the square of the residuals as a function of iteration number. Reynolds number 1000, 2303 equations.
of change of • had fallen below a heuristically determined tolerance, the value of A was set to zero. For a A of zero, this algorithm is then mathematically equivalent to using the modified Newton method, hence we used Algorithm (a) for the final iterations, i.e. until the problem converges to the correct solution. This is done for two reasons. The first is to save computation time, for when the rate of change in • has fallen below a certain tolerance, the current estimate of U is sufficiently accurate to achieve convergence of the standard steadystate problem using a modified Newton method. Since this is true, it is pointless as well as time consuming to compute the modified Jacobian and the value of • at each iteration. The second problem which can arise is the previously mentioned oscillatory behaviour about • which was observed to occur when the rate of change of • was small. This may be caused by local instability and can lead to "thrashing" of the algorithm with non convergence. This behaviour which was not studied in depth will be addressed in the near future.
$. Numerical results The preceeding problem was solved on a Control Data Cyber 180-830 computer for Reynolds numbers between 400 and 1000. The number of equations varied from 1107 to 3576.
D. Fauchon a al., A finite element computation of moderate Reynolds
145
Initial estimates of a zero-velocity field and the velocity field resulting from solution of the Stokes problem at Re -0.001 were both tested. The comparison of the two algorithms tested is based on the execution time required to reach the solution. The change of the residual • with number of iterations is illustrated in Fig. 2. It can be seen that, at the beginning of the problem, • is essentially zero since the only contribution to comes from the boundary conditions (U and P were initially taken equal to zero). As the iteration process continues, the algorithm adjusts automatically the value of ,~ so as to search for a good estimate of U and P (Step 7) of Algorithm (b)). This stage corresponds to an increase of O. Once a satisfactorj g,Jess is found, O decreases but is subjected to oscillations while the problem is converging. The large vana~.~.."~':~-of O as well as the observed potential for oscillatory behaviour led to the conclusion that using O as a stopping criterion was not an optimal way of proceeding. Instead we used the relative value of • between iterations, i.e. ~!~ The convergence behaviour of the Newton-Marquardt method is observed to be strongly influenced by the number of equations to be solved, the Reynolds number at which the solution is desired, and the initial guess. Figures 3-5 show the change in value of the residual versus iterations for various Reynolds numbers and the number of equations. In each case the nonzero-velocity initial guess corresponds to a velocity field at Reynolds number of 0.001, which is essentially a Stokes problem. The value of • ~ / O is seen to be smaller in the case of a nonzero-velocity guess at each Reynolds number. The implication of this is t ~ t the maximum e
1.80 1.60 1.40 1.20 %~" 1.00
~' 0.80 0.60
° ~O/oo0"
0.40 0.20 0
I
2
3
4
5
Number of iterotions Fig. 3. Relative residual as a function of the iteration number. Reynolds number 400, 3576 equations.
146
D. Fauchon et al., A finite element computation of moderate Reynolds
s
I
4 b
2
!
I
2
5
4
Humber of iterations Fig. 4. Relative residual as a function of the iteration number. Reynolds number 1000, 1303 equations.
Reynolds number at which the Newton-Marquardt method is effective, is increased for any given number of equations. This is shown in Fig. 6, which is a plot of the maximum Reynolds number at which a solution could be obtained with the Newton-Marquardt method versus the number of equations. We observe that the nonzero-velocity guess leads to superior performance. The price to pay for this improved robustness is the cost of solving the Stokes problem to calculate the initial velocity field. The decrease in the maximum Reynolds number with increasing number of equations can be related to the influence of the A1 matrix on the Jacobian matrix. As the number of equations is increased, the matrix bandwidth is enlarged, therefore the influence of the main diagonal is diminished, thus reducing the effectiveness of this method. The performance of the two algorithms are compared directly in Figs. 7 and 8. In Fig. 7 curves of the execution time versus the number of equations are given for a Reynolds number of 400. The Newton-Marquardt method is consistently better, decreasing the execution time by approximately 75%. The rather unusual result of a decrease in execution time for an increase in the number of equations observed at the beginning of the modified Newton curve results from the tolerance chosen to satisfy convergence. At 1107 equations one additional iteration was required, because AU/U was slightly higher than the tolerance and the final solution has, in fact, a much lower error. Figure 8 shows the execution time versus the Reynolds number for 1303 equations. We observe from this graph that at low Reynolds number (<2) the modified Newton method is superior to the Newton-Marquardt, but as the Reynolds number increases the difference between the two algorithms increases until, at a
1.80
e60
"' '
1.40
I
1.20
~) 1.00 ~'~" 0.80 0.60 0.400.20
0
I
'"
2
3
4
5
Numberof iterations
6
Fig. 5. Relative residual as a function of the iteration number. Reynolds number 1000, 2303 equations.
1200
,
m
....
I guess ,.
•- 800 E w
c o
~. 600
4OO
i
20~)1000
1500
2(XX)
2500 3000 Numberof equations
3500
Fig. 6. Maximum Reynolds number possible as a function of the number of equations.
147
4()00
148
D. Fauchon et al., A finite element computation of moderate Reynolds
- Newton ,o
f ! This work
i 103 U
ioZl
I
IOOO
1500
2O00
25OO
I
3000
3500
4000
Number of equations
Fig. 7. Execution time as a function of the number of equations. Reynolds number 400. 1600
1300
I000
[ u
700
4OO
I00
to-I
j
~ ~ ~,.
~-~'f
I00
......
l
i01 Reynolds number
i02
Fig. 8. Execution time as a function of the Reynolds number, 1303 equations.
I0 3
D. Fauchon
et
al., A finite element computation of moderate Reynolds
149
Reynolds number of 1000, the Newton-Marquardt method gets approximately 3 times faster. The trend is, moreover, clearly in favour of the Newton-Marquardt method, as the slope of the curve is not increasing at so fast a rate.
6. Conclusions and recommendations The Newton-Marquardt method is observed to provide a robust algorithm for the solution of the Navier-Stokes equations by the finite element method up to a Reynolds number of 1000. This algorithm was shown to be more efficient than the classical modified Newton method and can result in considerable cost savings. Furthermore, this algorithm exhibits a degree of self-opt'mization which makes it ideal for implementation. Further work may be desirable in order to expand the range of Reynolds numbers for which the method is effective.
Acknowledgment The financial assistance of Natural Science and Engineering Research Council of Canada and Control Data Corporation is gratefully acknowledged.
References [1] [2] [3] [4]
S.H. Crandail, Engineering Analysis (McGraw-Hill, New York, 1956). R.L. Burden and J.D. Faires, Numerical Analysis (Prindle, Weber and Smidt, Boston, 3rd. ed., 1985). E.J. Dennis and J. Mor~, Quasi-Newton methods, motivation and theory, SIAM Rev. 19 (19/7) 46-89. H. Matthies and G. Strang, The solution of nonlinear finite element equations, lnternat. J. Numer. Meths. Engrg. 14 (1979) 1613-1626. [5] R.E. Hayes and EA. Tanguy, Solution of convection-diffusion differential equations using Marquardt's method, Canad. J. Chem. Engrg. 65 (1987) 522-525. [6] D.W. Marquardt, Estimation of non-linear parameters, J. Soc. Ind. Appl. Math 11 (1963) 431-441. [7] K. Bathe, Finite Element Procedures in Engineering Analysis (Prentice-Hall, Engiewood Cliffs, NJ, 1982). [8] M. Crouzeix and P.A. Raviart, Conforming and non-conforming finite element methods for solving the stationary Stokes equations, RAIRO Anal. Numer. 7 (1973) 33. [9] EA. Tanguy, M. Fortin and L. Choplin, Finite element simulation of dip coating, II: Non-Newtonian fluids, Intemat. J. Numer. Meths. Fluids 4 (1984) 459--475.