Newton Linearization of The Incompressible Navier-Stokes Equations
Newton Linearization of The Incompressible Navier-Stokes Equations
Newton Linearization of The Incompressible Navier-Stokes Equations
Correspondence to: Tony W. H. Sheu, Department of Engineering Science and Ocean Engineering, National Taiwan
University, 73 Chow-Shan Road, Taipei, Taiwan.
E-mail: [email protected]
Ph.D. Candidate.
Contract=grant sponsor: National Science Council; contract=grant number: 90-2811-E-002-008
Received 27 May 2002
Copyright
?
2004 John Wiley & Sons, Ltd. Revised 14 April 2003
298 T. W. H. SHEU AND R. K. LIN
has been successfully applied to solve for non-linear equations for uid dynamics and heat
transfer. In Newtons method, equation Ax =b, where x is the state vector and A is a function
of x, is approximated by a rst-order Taylor series to render Jx =R. Here, the component of
the Jacobian matrix (or the Frechlet derivative of R [3]) is dened by J(i; j) = @R(i)=@x(j),
where R is known as the residual vector R=b Ax. The matrix equation is then solved for
x
k
from Jx
k
=R
k
, followed by obtaining x
k+1
=x
k
+ x
k
. This is done until the residual
norm falls below a users specied tolerance, with J and R being computed each time using
the mostly updated solution x. Some practical issues of inverting the huge Jacobian matrix J
were discussed in [4].
Newtons method is potentially attractive in accelerating convergence due to its ability to
oer q-quadratic convergence [2]. The implementation of this method involves, however, time-
consuming manipulation of the Frechlet derivative of R at the current solution vector x
k
. This
disadvantage has been addressed previously by Hunt [5] in his three-dimensional viscous ow
analysis. The NewtonKrylov method was proposed to avoid calculation and, thus, storage
of the Frechlet derivative [3]. The underlying principle of this method is to minimize the
residual in a Krylov space by linearizing the equation using the Newton method and solve
the resulting linear system of algebraic equations with a Krylov method. Both Arnoldi- and
Lanczos-based iterative matrix solution solvers can, thus, be applied [6]. Since the Jacobian
matrix is neither formed nor stored, the NewtonKrylov method can be specically refereed
to as a matrix-free NewtonKrylov method [7]. An assessment study of several matrix-free
NewtonKrylov methods can be seen in the work of McHugh and Knoll [8]. For a further
minimization of residual, a multigrid preconditioner can be used together with the Newton
Krylov linearization procedure [7].
Another major challenge in using Newtons method, assuming the required memory is
available, is the increasing radius of convergence. The variable secant procedure, known
more generally as the NewtonRaphson method, was proposed to overcome this diculty
by replacing the derivative with a secant line through two points. This potentially attractive
linearization method, however, requires factorization of tangent matrix at each iteration [9]. In
fact, the need to solve a large-scale system of linear equations at each iteration is considered as
a major shortcoming of the Newton-family methods. To reduce memory requirement and com-
putational cost when performing a classical NewtonRaphson linearization method, one can
introduce some iterative means to approximately solve the linear system of Newton lineariza-
tion equations. We refer to this class of methods as the inexact Newton methods [1012]. In
the modied NewtonRaphson method, the tangent matrix is factorized only once for a num-
ber of steps. This occasionally updating strategy may lead to poor convergence for a highly
non-linear system. As a means of partly circumventing this problem, the asymptotic Newton
method was proposed [9]. The Newton-relaxation method [13], with Newton being the pri-
mary iteration and relaxation the secondary iteration, is another useful method to avoid direct
solution of a large-scale linear system of three-dimensional equations, which must be solved
at each iteration of Newtons method. In the light of the above literature survey, Newtons
linearization procedure which is suited to be used together with the discretization scheme and
the solution algorithm is proposed.
The rest of this paper is organized as follows. In Section 2 the Newton linearization pro-
cedure, used together with the iterative solution algorithm, is detailed. This is followed by
proposing a discretization scheme mostly suited for solving the resulting linearized momentum
equations. Discretization of incompressible NavierStokes equations on non-staggered grids is
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
NEWTON LINEARIZATION OF INCOMPRESSIBLE NAVIERSTOKES 299
derailed in Section 4. An assessment study on the present method and the simple iterative
update coecient method is given in Section 5. The last section provides the concluding
remarks.
2. LINEARIZATION OF NAVIERSTOKES EQUATIONS
We consider in the paper a two-dimensional steady-state ow equations. Subject to the in-
compressible constraint condition, the transport equations for the viscous uid ow in are
as follows:
u
x
+ v
y
=0 (1)
uu
x
+ vu
y
=p
x
+ (u
xx
+ u
yy
) (2)
uv
x
+ vv
y
=p
y
+ (v
xx
+ v
yy
) (3)
In this paper, the subscript denotes the partial derivative. Here, u =(u; v) and p are velocity
vector and pressure, respectively. Note that Equations (2) and (3) serve as the transport
equations for u and v, respectively. Working equation for p can be obtained as follows by
summing the @=@x (2) and @=@y (3) and employing (1):
p
xx
+ p
yy
=[(u
x
)
2
+ 2u
y
v
x
+ (v
y
)
2
] (4)
Note that the inaccuracy stemming from approximation of terms shown in the right-hand side
of (4) for a given velocity eld may limit the convergence rate, as discussed [14]. The above
pressure Poisson equation needs to be supplemented by the boundary condition given by
p
n
= [(u )u +
2
u + f] n (5)
where n is the outward-directed unit vector normal to the boundary of . In what follows the
dynamic viscosity of the uid ow is considered uniform for simplicity.
Linearization of convective terms on the left-hand sides of momentum equations (2) and
(3) starts from rewriting them as
(u
2
)
x
+ (uv)
y
=p
x
+ (u
xx
+ u
yy
) (6)
(uv)
x
+ (v
2
)
y
=p
y
+ (v
xx
+ v
yy
) (7)
Consider a function st, we can expand it in a Taylor series about the current value and
terminate the series expansion after the rst-derivative terms. The result is as follows:
s
k+1
t
k+1
=s
k
t
k
+
_
@
@s
(st)
k
_
(s
k+1
s
k
) +
_
@
@t
(st)
k
_
(t
k+1
t
k
) + H:O:T
=s
k+1
t
k
+ s
k
t
k+1
s
k
t
k
+ H:O:T (8)
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
300 T. W. H. SHEU AND R. K. LIN
In the derivation that follows, all variables denoted by the superscript k are evaluated using
solutions obtained at the previous iteration counter. As for terms with the superscript k + 1,
they are evaluated at the most updated iteration and are, therefore, referred to as the active
quantities. According to Equation (8), we can linearize (u
2
)
k+1
x
and (uv)
k+1
y
as
(u
2
)
k+1
x
=(u
k+1
u
k
+ u
k
u
k+1
u
k
u
k
)
x
=u
k+1
x
u
k
+ u
k+1
u
k
x
+ u
k
x
u
k+1
+ u
k
u
k+1
x
u
k
x
u
k
u
k
u
k
x
(9a)
(uv)
k+1
y
=(u
k+1
v
k
+ u
k
v
k+1
u
k
v
k
)
y
=u
k+1
y
v
k
+ u
k+1
v
k
y
+ u
k
y
v
k+1
+ u
k
v
k+1
y
u
k
y
v
k
u
k
v
k
y
(9b)
Substituting (9a) and (9b) into (6) led us to derive the linearized x-momentum equation as
follows:
u
k
u
k+1
x
+ v
k
u
k+1
y
(u
k+1
xx
+ u
k+1
yy
) + u
k
x
u
k+1
=p
k+1
x
+ u
k
u
k
x
+ v
k
u
k
y
u
k
y
v
k+1
(10)
Similarly, one can derive the following convectiondiusionreaction (CDR) equation for v
u
k
v
k+1
x
+ v
k
v
k+1
y
(v
k+1
xx
+ v
k+1
yy
) + v
k
y
v
k+1
=p
k+1
y
+ u
k
v
k
x
+ v
k
v
k
y
v
k
x
u
k+1
(11)
Neglect of the underlined terms from the Newton linearization Equations (10) and (11) results
in the conventional lagging coecient linearized equations.
For computational eciency, we can solve for Equation (10), for example, iteratively by
virtue of the following Alternating Direction Implicit (ADI) steps [15]:
u
k
u
k+1
x
u
k+1
xx
+ u
k
x
u
k+1
=p
k+1
x
+ v
k
u
k+1
y
u
k+1
yy
+ f
1
(12a)
v
k
u
k+1
y
u
k+1
yy
+ u
k
x
u
k+1
=p
k+1
x
+ u
k
u
k+1
x
u
k+1
xx
+ f
2
(12b)
where
f
1
=u
k
u
k
x
+ v
k
u
k
y
u
k
y
v
k+1
(13a)
f
2
=u
k
u
k
x
+ v
k
u
k
y
u
k
y
v
k+1
(13b)
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
NEWTON LINEARIZATION OF INCOMPRESSIBLE NAVIERSTOKES 301
Note that when an ADI method is applied together with the pseudo-transient approach, where
pseudo-time derivative terms are added to the steady NavierStokes equations so as to be
able to parabolize the elliptic dierential system for time marching, its strong convergence
property discussed in Reference [16] deteriorates considerably for cases involving complex ge-
ometries and high Reynolds numbers [17]. Under these circumstances, use of an ADI method
is constrained by a small time increment to maintain convergence solutions to steady-state.
3. NUMERICAL MODEL
As Equations (12a) and (12b) show, the prototype equation takes the following form:
u
x
+ v
y
k(
xx
+
yy
) + c=f (14)
For simplicity, the model equation is solved subject to a specied boundary value of . In
the above, k and c denote the diusion coecient and the reaction coecient, respectively.
In what follows u, v, k and c are assumed to have constant values.
By virtue of the operator splitting method of Peaceman and Rachford [15], solutions to
Equation (14) are sought from the predictor and corrector steps, respectively:
u
x
k
xx
+ c
=f
1
(15a)
v
n+1
y
k
n+1
yy
+ c
n+1
=f
2
(15b)
In the above, f
1
=f
v
n
y
k
n
yy
and f
2
=f
n+1
u
x
k
xx
. As Equations (15a) and
(15b) reveal, a key concern in the analysis of the two-dimensional CDR equation (14) is the
discretization of the following one-dimensional equation:
u
x
k
xx
+ c=f (15c)
For illustrative purposes, f is assumed to be a known constant.
Our strategy of approximating (15c) is to employ its general solution
=c
1
e
1
x
+ c
2
e
2
x
+
f
c
(16)
where (
1
;
2
) =(u +
u
2
+ 4ck=2k; u
u
2
+ 4ck=2k). In Equation (16), c
1
and c
2
are con-
stants. Terms other than the diusive term shown in Equation (15c) are approximated by the
centre-like scheme. Therefore, the discrete equation at an interior node i can be expressed as
_
u
2h
m
h
2
+
c
6
_
i1
+ 2
_
m
h
2
+
c
3
_
i
+
_
u
2h
m
h
2
+
c
6
_
i+1
=f (17)
where h is the mesh size. We then substitute the exact solutions
i
=c
1
e
1
x
i
+ c
2
e
2
x
i
+ f=c,
i+1
=c
1
e
1
h
e
1
x
i
+c
2
e
2
h
e
2
x
i
+f=c, and
i1
=c
1
e
1
h
e
1
x
i
+c
2
e
2
h
e
2
x
i
+f=c into Equation (17).
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
302 T. W. H. SHEU AND R. K. LIN
The closed-form of m can be derived as [18]
m=h
2
_
c=3 + c=6 cosh(
1
) cosh(
2
) + u=2h sinh(
1
) cosh(
2
)
cosh(
1
) cosh(
2
) 1
_
(18)
where (
1
;
2
) =(uh=2k;
_
(uh=2k)
2
+ ch
2
=k). Note that the predicted inaccuracy stems solely
from the approximated f.
4. INCOMPRESSIBLE NAVIERSTOKES CALCULATION
ON NON-STAGGERED GRIDS
On physical grounds, p in the equations of motion is discretized by a centred scheme.
However, centre approximation of @p=@x and @p=@y on non-staggered grids engenders spurious
evenodd oscillations [19, 20]. Therefore, one has to suppress these erroneous checkerboard-
ing pressures when simulating the incompressible ow equations on grids of the simplest
form [21].
For overcoming diculty with the evenodd decoupling, we calculate F
j
(h
x
) and
G
j
(h
2
xx
) implicitly from
0
F
j+1
+
0
F
j
+
0
F
j1
=a
0
(
j+2
j+1
) + b
0
(
j+1
j
)
+c
0
(
j
j1
) + d
0
(
j1
j2
) (19)
and
1
G
j+1
+
1
G
j
+
1
G
j1
=a
1
j+2
+ b
1
j+1
+ c
1
j
+ d
1
j1
+ e
1
j2
(20)
Provided that (
0
;
0
;
0
; a
0
; b
0
; c
0
; d
0
) =(
1
5
;
3
5
;
1
5
;
1
60
;
29
60
;
29
60
;
1
60
) and (
1
;
1
;
1
; a
1
; b
1
; c
1
; d
1
; e
1
) =
(1;
11
2
; 1;
3
8
; 6;
51
4
; 6;
3
8
), both
x
and
xx
accommodate sixth-order accuracy.
The implicit equations for F and G at nodes immediately adjacent to the boundary can be
derived by specifying d
0
=e
1
=0 and a
0
=a
1
=0 at nodes next to the left and right boundaries,
respectively. By performing Taylor series expansion, the coecients can be analytically de-
rived as (
0
;
0
;
0
; a
0
; b
0
; c
0
; d
0
) =(
3
10
;
3
5
;
1
10
;
1
30
;
19
30
;
1
3
; 0) and (
1
10
;
3
5
;
3
10
; 0;
1
3
;
19
30
;
1
30
) at nodes next
to the left and right boundaries, respectively. In addition, coecients for evaluating G
j
are
exactly derived as (
1
;
1
;
1
; a
1
; b
1
; c
1
; d
1
; e
1
) =(1; 10; 1; 0; 12; 24; 12; 0).
5. NUMERICAL RESULTS
5.1. Non-linear advectiondiusion scalar equation
To verify the proposed Newton linearization method, the following two-dimensional non-linear
convectiondiusion equation for u is considered in 06x; y61:
uu
x
+ vu
y
k(u
xx
+ u
yy
) =f (21)
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
NEWTON LINEARIZATION OF INCOMPRESSIBLE NAVIERSTOKES 303
k-th nonlinear iteration
10 20 30 40 50 60
10
-12
10
-11
10
-10
10
-9
10
-8
10
-7
10
-6
10
-5
10
-4
Relaxation; = 0.2
= 0.4
= 0.6
= 0.8
Newton
{
R
e
s
i
d
u
a
l
Figure 1. Comparison of the convergence histories, with the initial guess value u =0:5, for solving the
non-linear advectiondiusion equation (21).
Figure 2. Plots of the ADI iteration number against the outer (non-linear) iteration number
for the two investigated linearization methods.
The validation test was performed at k =x
2
, v =y and f=2x
3
(y
4
x). The solution to
equation (21) was exactly derived as u =x
2
y
2
.
We assess our proposed model and then the standard relaxation method based on u
new
=
u
new
+ (1 )u
old
, where 0661. As Figure 1 shows, a considerable amount of non-linear
iterations has been saved in view of the iteration numbers needed for the results obtained
at =0:2, 0.4, 0.6 and 0.8. The tolerance, dened as [1=N
(u
new
u
old
)
2
]
1=2
, set for each
calculation is 10
13
. Here, N denotes the number of solution points. We also compare the
iteration numbers needed to reach the convergent ADI solution at each non-linear iteration.
It is seen from Figure 2 that much fewer ADI iterations are needed when using the present
Newton linearization method.
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
304 T. W. H. SHEU AND R. K. LIN
5.2. Non-linear NavierStokes equations
For the sake of validation, a problem with f =0 is considered. The analytic pressure in the
unit square is
p=
2
(1 + x)
2
+ (1 + y)
2
(22)
provided that the boundary velocities are analytically specied as
u =
2(1 + y)
(1 + x)
2
+ (1 + y)
2
(23a)
v =
2(1 + x)
(1 + x)
2
+ (1 + y)
2
(23b)
In Figure 3, we plot the computed rates of convergence for u, v, and p according to
C=
ln ||err
1
|| ln ||err
2
||
ln |h
1
| ln |h
2
|
(24)
The error measure is cast in the discrete L
2
-norm
E=
_
1
M
M
1=2
i; j=1
(u
ij
U
ij
)
2
_
1=2
(25)
In the above equation, U
ij
denotes the exact solution at an interior nodal point (i; j) and u
ij
is the corresponding nite-dierence solution. As the computed rates of convergence show
in Figure 3, the validity of the method is justied. Like the scalar problem, it is found that
Figure 3. The computed rates of convergence for u, v and p.
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
NEWTON LINEARIZATION OF INCOMPRESSIBLE NAVIERSTOKES 305
Figure 4. Comparison of the convergent histories for solving the non-linear NavierStokes problem,
which has the analytic solutions given in (22)(23b) at Re =1000. The initial guess solutions for u
and v are u =v =0:5: (a) convergence histories for u; and (b) convergence histories for v.
Figure 5. Plots of the ADI iteration number against the outer (non-linear) iteration number for the two
investigated linearization methods: (a) u; and (b) v.
a considerable amount of non-linear and ADI iterations can be saved, as seen in Figures 4
and 5.
The second problem to be investigated is known as the Kovasznay ow problem [22],
which is amenable to the analytic solutions given below
u =1 e
x
cos(2y) (26a)
v =
2
e
x
sin(2y) (26b)
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
306 T. W. H. SHEU AND R. K. LIN
Figure 6. The computed rates of convergence for u, v and p.
Figure 7. Comparison of the convergent histories for solving the non-linear NavierStokes problem,
which has the analytic solutions given in (27) at Re =1000. The initial guess solutions for u and v are
u =v =0:5: (a) convergence histories for u; and (b) convergence histories for v.
p=
1
2
(1 e
2x
) (26c)
where =Re=2 (Re
2
=4 + 4
2
)
1=2
. Numerical calculations have been carried out in a square
which is covered with uniform grids. For the test Reynolds number 1000, both pressure and
velocity elds are well-predicted, as seen from the predicted errors shown in Figure 6. As the
error reduction plot shows in Figure 7, non-linear iteration numbers have been considerably
saved. The inner iteration number for each non-linear iteration is also largely reduced, as seen
in Figure 8. This clearly shows the advantage of using the Newton linearization method.
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
NEWTON LINEARIZATION OF INCOMPRESSIBLE NAVIERSTOKES 307
Figure 8. Plots of the ADI iteration number against the outer (non-linear) iteration number for the two
investigated linearization methods: (a) u; and (b) v.
Figure 9. Comparison of the convergent histories for solving the non-linear NavierStokes problem,
which has the analytic solutions given in (28)(30) at Re =1000. The initial guess solutions for u
and v are u =0:5: (a) convergence histories for u; and (b) convergence histories for v.
The validation and assessment are followed by considering another analytic lid-driven cavity
ow problem [23]. In a square domain, the NavierStokes equations are solved subject to the
following boundary conditions for u and v at x =0, 1 and y =0, 1:
u =8(x
4
2x
3
+ x
2
)(4y
3
2y) (27a)
v =8(4x
3
6x
2
+ 2x)(y
4
y
2
) (27b)
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
308 T. W. H. SHEU AND R. K. LIN
k-th nonlinear iteration
A
D
I
i
t
e
r
a
t
i
o
n
25 50 75
10
0
10
1
10
2
10
3
Relaxation; = 0.2
= 0.4
= 0.6
= 0.8
Newton
{
(a) (b)
Figure 10. Plots of the ADI iteration number against the outer (non-linear) iteration number for the two
investigated linearization methods: (a) u; and (b) v.
h
L
2
-
e
r
r
o
r
n
o
r
m
0.04 0.05 0.06 0.07 0.08 0.09 0.1
0.0001
0.0002
0.0003
0.0004
0.0005
0.0006
u
v
p
ra
te
o
f
c
o
n
v
e
rg
e
n
c
e
=
2
.3
0
4
ra
te
o
f
c
o
n
v
e
rg
e
n
c
e
=
1
.9
9
4
ra
te
o
f
c
o
n
v
e
rg
e
n
c
e
=
2
.3
6
9
Figure 11. The computed rates of convergence for u, v and p.
Moreover, as the body force f =(f
1
; f
2
) is given as
f
1
=0 (28a)
f
2
=
Re
8
[24J
1
(x) + 2I
1
(x)I
2
(y) + I
1
(x)I
2
(y)] + 64[J
3
(x)J
4
(y) I
2
(y)I
2
(y)J
2
(x)]
(28b)
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
NEWTON LINEARIZATION OF INCOMPRESSIBLE NAVIERSTOKES 309
Figure 12. The computed solutions at Re =1000: (a) streamlines; (b) pressure contours; and (c)
mid-sectional velocity proles for u and v.
the pressure solution takes the following analytic form:
p=
8
Re
[J
1
(x)I
2
(y) + I
1
(x)I
2
(y)] + 64J
3
(x)[I
2
(y)I
2
(y) (I
2
(y))
2
] (29)
where
I
1
(x) =x
4
2x
3
+ x
2
I
2
(y) =y
4
y
2
J
1
(x) =0:2x
5
0:5x
4
+
1
3
x
3
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
310 T. W. H. SHEU AND R. K. LIN
Figure 13. A comparison of the computed and Ghias velocity proles
for u(x; 0:5) and v(0:5; y) at Re =5000.
Figure 14. Comparison of the convergent histories for solving the lid-driven cavity ow prob-
lem at Re =5000. The initial guess solutions for u and v are u =v =0:5: (a) convergence
histories for u; and (b) convergence histories for v.
J
2
(x) =4x
6
+ 12x
5
14x
4
+ 8x
3
2x
2
J
3
(x) =0:5(x
4
2x
3
+ x
2
)
2
J
4
(y) =24y
5
+ 8y
3
4y
Considering the case with Re =1000, employment of the Newton linearization method
renders a much faster convergent solution. The evidences are given in Figures 9 and 10,
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
NEWTON LINEARIZATION OF INCOMPRESSIBLE NAVIERSTOKES 311
Figure 15. Plots of the ADI iteration number against the outer (non-linear) iteration number for the two
investigated linearization methods: (a) u; and (b) v.
from which a considerable amount of CPU time is saved. The computed solutions given in
Figure 11 reveal that the proposed model oers good accuracy but not at the cost of de-
teriorated convergence. For completeness, the contour levels for streamline and pressure are
plotted in Figure 12, together with the mid-sectional velocity proles for u and v.
5.3. Lid-driven cavity ow problem
The next NavierStokes problem considers the cavity ow driven by a constant upper lid
velocity u
lid
. The geometrical simplicity and physical complexity have made this problem
an attractive test for benchmarking the incompressible NavierStokes models. With as the
characteristic length, u
lid
the characteristic velocity, the Reynolds number under investigation
is 5000.
We continuously rene the mesh and plot the grid independent mid-plane velocity proles
u(0:5; y) and v(x; 0:5) in Figure 13. For the sake of comparison, the steady-state benchmark
solutions obtained by Ghia [24] are also given in the same gure. The agreement between the
two numerical solutions is extremely good. Most importantly, much improved convergence
histories are seen from Figures 14 and 15 and these conrm the applicability of the proposed
scheme.
6. CONCLUDING REMARKS
The objective of this study is to show the eectiveness of using the Newton linearization
method to solve for the incompressible NavierStokes equations. Revealed from this study
is that the linearized equations can be eciently solved on non-staggered grids using the
computationally very accurate CDR scheme. Numerical study of several problems shows the
eectiveness of Newtons method in oering much faster outer iteration (or non-linear itera-
tion) and inner iteration (ADI iteration) convergence to the convergent solutions is seen for
all test problems.
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312
312 T. W. H. SHEU AND R. K. LIN
ACKNOWLEDGEMENTS
We acknowledge the nancial support from National Science Council under NSC 90-2811-E-002-008.
This work was partially accomplished in the course of the rst authors sabbatical leave in University of
Paris 6. The facilities provided by Professors Oliver Pironneau and Yvon Maday are highly appreciated.
A list of useful references suggested by the anonymous reviewer is gratefully acknowledged.
REFERENCES
1. Galpin PF, Raithby GD. Treatment of non-linearities in the numerical solution of the incompressible Navier
Stokes equations. International Journal for Numerical Methods in Fluids 1986; 6:409426.
2. Dennis Jr JE, Schmabel RB. Numerical Methods for Unconstrained Optimization and Nonlinear Equations.
Prentice-Hall: Englewood Clis, NJ, 1983.
3. Pereira JMC, Kobayahi MH, Pereira JCF. A fourth-order-accurate nite volume compact method for the
incompressible NavierStokes solutions. Journal of Computational Physics 2001; 167:217243.
4. Hunt R. The numerical solution of the ow in a general bifurcating channel at moderately high Reynolds
number using boundary-led co-ordinates, primitive variables and Newton iteration. International Journal for
Numerical Methods in Fluids 1993; 17:711729.
5. Hunt R. Three-dimensional steady ow in a divergent channel using nite and pseudo spectral dierences.
SIAM Journal on Scientic Computing 1996; 17(3):561578.
6. McHugh PR, Knoll DA. Fully coupled nite volume solutions of the incompressible NavierStokes and energy
equations using an inexact Newton method. International Journal for Numerical Methods in Fluids 1994;
19:439455.
7. Pernice M, Tocci MD. A multigrid-preconditioned NewtonKrylov method for the incompressible NavierStokes
equations. SIAM Journal on Scientic Computing 2001; 23(2):398418.
8. McHugh PR, Knoll DA. Comparison of standard and matrix-free implementations of several NewtonKrylov
solvers. AIAA Journal 1994; 32(12):23942400.
9. Hadji S, Dhatt G. Asymptotic-Newton method for solving incompressible ows. International Journal for
Numerical Methods in Fluids 1997; 25:861878.
10. Dembo RS, Eisenstat SC, Steihaug T. Inexact Newton methods. SIAM Journal on Numerical Analysis 1982;
19:400408.
11. Eisenstat SC, Walker HF. Globally convergent inexact Newton methods. SIAM Journal on Optimization 1994;
4:393422.
12. Eisenstat SC, Walker HF. Choosing the Newton method. SIAM Journal on Scientic Computing 1996; 17:
1632.
13. Sheng C, Taylor LK, Whiteld DL. Multigrid algorithm for three-dimensional incompressible high-Reynolds
number turbulent ows. AIAA Journal 1995; 33(11):20732079.
14. Chaviaropoulos P, Giannakoglou K. A vorticity-streamfunction formulation for steady incompressible two-
dimensional ows. International Journal for Numerical Methods in Fluids 1996; 23:431444.
15. Peaceman DW, Rachford HH. The numerical solution of parabolic and elliptic dierential equations. Journal of
the Society for Industrial and Applied Mathematics 1955; 3:2841.
16. Stone HL. Iterative solution of implicit approximation of multi-dimensional partial dierential equations. SIAM
Journal on Numerical Analysis 1968; 5:530558.
17. Jordan SA. An iterative scheme for numerical solution of steady incompressible viscous ows. Computers and
Fluids 1992; 21(4):503517.
18. Sheu TWH, Wang SK, Lin RK. An implicit scheme for solving the convectiondiusionreaction equation in
two dimensions. Journal of Computational Physics 2000; 164:123142.
19. Patankar SV. Numerical Heat Transfer and Fluid Flow. Hemisphere: New York, 1990.
20. Harlow FW, Welch JE. Numerical calculation of time-dependent viscous incompressible ow of uid with free
surfaces. The Physics of Fluids 1965; 8:21822189.
21. Rhie CM, Chow WL. Numerical study of the turbulent ow past an airfoil with trailing edge separation. AIAA
Journal 1983; 21:15251532.
22. Kovasznay LIG. Laminar ow behind a two-dimensional grid. Proceedings of Cambridge Philosophical Society
1948, 44.
23. Shih TM, Tan CH, Hwang BC. Eects of grid staggering on numerical schemes. International Journal for
Numerical Methods in Fluids 1989; 8:193212.
24. Ghia U, Ghia KN, Shin CT. High-Re solutions for incompressible ow using the NavierStokes equations and
a multigrid method. Journal of Computational Physics 1982; 48:387411.
Copyright ? 2004 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2004; 44:297312