Numerical Methods in Fluids - 1999 - Zienkiewicz - The Characteristic Based Split Procedure An Efficient and Accurate

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS
Int. J. Numer. Meth. Fluids 31: 359 – 392 (1999)
THE CHARACTERISTIC-BASED-SPLIT PROCEDURE: AN

EFFICIENT AND ACCURATE ALGORITHM FOR FLUID
PROBLEMS
O.C. ZIENKIEWICZa,*, P. NITHIARASUa, R. CODINAb, M. VÁZQUEZb AND

P. ORTIZc
a
Institute for Numerical Methods in Engineering, Department of Ci6il Engineering, Uni6ersity of Wales Swansea,
Swansea SA2 8PP, UK
b
International Center for Numerical Methods in Engineering, Uni6ersidad Politécnica de Catalunya, Barcelona, Spain
c
Centro de Estudios de Técnicas Aplicadas, Cedex Alfonso XII, 3, E-28014, Madrid, Spain
SUMMARY
In 1995 the two senior authors of the present paper introduced a new algorithm designed to replace the
Taylor–Galerkin (or Lax–Wendroff) methods, used by them so far in the solution of compressible flow
problems. The new algorithm was applicable to a wide variety of situations, including fully incompress-
ible flows and shallow water equations, as well as supersonic and hypersonic situations, and has proved
to be always at least as accurate as other algorithms currently used. The algorithm is based on the
solution of conservation equations of fluid mechanics to avoid any possibility of spurious solutions that
may otherwise result. The main aspect of the procedure is to split the equations into two parts, (1) a part
that is a set of simple scalar equations of convective – diffusion type for which it is well known that the
characteristic Galerkin procedure yields an optimal solution; and (2) the part where the equations are
self-adjoint and therefore discretized optimally by the Galerkin procedure. It is possible to solve both the
first and second parts of the system explicitly, retaining there the time step limitations of the Taylor –
Galerkin procedure. But it is also possible to use semi-implicit processes where in the first part we use a
much bigger time step generally governed by the Peclet number of the system while the second part is
solved implicitly and is unconditionally stable. It turns out that the characteristic-based-split (CBS)
process allows equal interpolation to be used for all system variables without difficulties when the
incompressible or nearly incompressible stage is reached. It is hoped that the paper will help to make the
algorithm more widely available and understood by the profession and that its advantages can be widely
realised. Copyright © 1999 John Wiley & Sons, Ltd.
KEY WORDS: CBS algorithm; fluid problems; convective– diffusion-type
1. INTRODUCTION
In recent years, the current authors have published several papers concerning the basis and
applications of a universal algorithm for practically all fluid mechanics problems, but no
specific name has been given to it [1 – 7]. From now on, it shall be referred to as the
* Correspondence to: Institute for Numerical Methods in Engineering, Department of Civil Engineering, University of
Wales Swansea, Swansea SA2 8PP, UK.
CCC 0271–2091/99/170359–35$17.50
Copyright © 1999 John Wiley & Sons, Ltd.
360 O.C. ZIENKIEWICZ ET AL.
characteristic-based-split (CBS) algorithm, which explains the main motivation for its origin.
The present paper summarises the theory of the algorithm as well as demonstrates its practical
applications.
The algorithm can be used in explicit, semi-implicit, nearly implicit and, indeed, in fully
implicit forms. Its advantages can be fully illustrated, however, using explicit and semi-implicit
forms and this paper will concentrate on these. The algorithm was initially based on a series
of preliminary studies conducted between 1990 and 1995 [8–12], but the split was introduced
correctly only later, in Reference [1].
Before the present algorithm was available, most successes of finite difference and finite
element methods in fluid mechanics were based on some variant of the Lax–Wendroff [13]
scheme, which by approximating better the time derivative also introduced a stabilisation of
the convective terms. The method that applied to finite elements become known as the
Taylor–Galerkin method [14,15]. Later it was discovered that for scalar variables, a direct
algorithm utilising the optimal approximation along the characteristics, the characteristic
Galerkin method, could be shown to be identical to the Taylor–Galerkin method [16,17].
The only similarity between the two procedures exists in the scalar case described by the
well-known convection – diffusion equation that is often used as a model for fluid mechanics
problems. However, for problems involving several variables, typical of fluid mechanics, the
application of the characteristic Galerkin method is not possible as only a single characteristic
speed must be involved. For this reason, the Taylor–Galerkin procedure has been used widely
even though giving a suboptimal approximation. Motivation for the development of the
present algorithm came from the fact that the characteristic Galerkin procedure for a scalar
variable is optimal in the sense of approximation and that, with suitable splitting, it can be
applied for the first stage of the solution of the fluid mechanics equation. The reminder of the
split being self-adjoint can then be solved optimally using the Galerkin procedure.
This split follows the general process initially introduced by Chorin [18,19] for incompress-
ible flow problems. Further motivation, which originated in the new algorithm, is based on an
additional benefit coming from the fact that, as already observed, for incompressible situa-
tions, the algorithm permitted equal (and indeed arbitrary) interpolations for all variables to
be used. This sidesteps the Babuska – Brezzi requirement in finite elements, and similarly all
difficulties in standard finite difference schemes. In the present form, the solution turns out to
be fully accurate with arbitrary interpolation for velocity and pressure for full
incompressibility.
This indeed widens the spectrum of solutions using the new algorithm and indeed it is shown
later how, in dealing with problems of solid dynamics, the present algorithm proves to be
capable of generalising such explicit approaches as those given by DYNA [20,21].
Although the original Chorin split [18,19] could never be used in a fully explicit code without
introducing artificial or real compressibility, the split provides a fully explicit algorithm that is
valid at least for steady state problems, even in a fully incompressible range. When compress-
ibility exists, computational advantages of the explicit form compare well with other currently
used schemes and the additional cost due to splitting the operator is insignificant. Generally,
for an identical cost, the results are considerably improved even in standard range of
aerodynamical problems. However, a further advantage is that both subsonic and supersonic
problems can be solved by the same code.
The algorithm involves time integration and, in general, the time step size will be limited by
the nature of the time stepping procedure adopted in each part of the split. If a fully explicit
procedure is adopted in the first part of the split then the time step is governed by a Courant
number defined in terms of the flow velocity u and viscosity n. If an explicit scheme is also
Copyright © 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 31: 359 – 392 (1999)
CHARACTERISTIC-BASED-SPLIT ALGORITHM 361
adopted for the second part of the split, then a Courant number depending on the compress-
ible wave celerity c is invoked and may control the time step size. We shall thus have several
possible categories of problems.
(a) Transonic and supersonic flows, fully explicit process

Here, in general, fully explicit computations are preferred as the time step limitation for both
parts of the split are similar.
(b) Low Mach number with low 6iscosity, semi-implicit process
Here, the flow is nearly incompressible and the time step used will be governed by first part
of the split, with the Courant number not affected much by viscosity.
(c) Low Mach number with high 6iscosity, nearly implicit process
In this range of flow, the limiting time step again is governed by the first part of the split, but
here viscosity may pose very severe limits. In such a regime, a ‘nearly’ implicit procedure is
recommended in which the viscous terms are treated implicitly.
(d) Shallow water equations, semi-implicit or nearly implicit
Here, fully explicit computations have been efficiently used by Peraire et al. [22]. However,
semi-implicit treatment can increase the time step in tidal computations by one (or even two)
orders of magnitude [23 – 25]. In such problems, if horizontal viscosity is large, a ‘nearly’
implicit process may be considered.
Typical examples for each of the categories will be shown later. Indeed, the procedure used
in (e) fully explicit dynamic computation in solids will be demonstrated. However, many further
applications of the CBS algorithm are possible and will be reported below.
2. THE SCALAR CONVECTION – DIFFUSION PROBLEM AND THE

CHARACTERISTIC GALERKIN EXPLICIT APPROXIMATION
Before proceeding with the description of the full algorithm, we shall recall the application of
the characteristic Galerkin method in the explicit form to a typical convection–diffusion
process with a scalar dependent variable f.
The governing equations can be written in a conservation form as
(V (Fi (Gi
+ + +Q = 0, (1)
(t (xi (xi
where xi is the ith co-ordinate (i =1, 2, 3),
Fi =ujf (2)
is the convected flux,
(f
Gi = − k (3)
(xi
is the diffusion flux,
Q = Q(x, t) (4)
is the source term, and
u= u(x, t) (5)
is the velocity field that is assumed to be known.
The full equation can thus be alternatively written as
(f
= −uj
(f
+
(
k
(f (u
−Q − f j =R(f), (6)
(t (xj (xi (xi (xj
in which only the first term on the right-hand-side is not self-adjoint. As that term corresponds
precisely to an advection wave moving with a velocity u, a change of co-ordinates to the
characteristic ones, given by
dx%i = dxi −ui dt, (7)
makes the offending term vanish leaving a fully self-adjoint system.
For such a self-adjoint system it is known that the standard Galerkin approximation in
space is optimal, but the inconvenience of a moving co-ordinate system is introduced.
However, this can be overcome with suitable remeshing and the procedure has been used in
many early solutions of the above equation. Here, Adye–Brrebia (1974) [26] appears to be first
followed by much later work [27 – 31]. While the exact co-ordinate transformation introduced
no error, simplified procedures using a Taylor approximation within the time step eliminated
the costly process of remeshing (or mesh interpolation), introducing, however, a time step
limitation [32]. This process is fully described in References [1,2] and the explicit form can be
written in fixed co-ordinates as
Dt 2 (
Df =f n + 1 − f n =DtR(f)n + u3 − ui (R(f))n + O(Dt 3), (8)
2 (xi
where 05 u3 5 1 and R(f) is defined by Equation (6).
This, as already mentioned, is of an identical form to that resulting from the Taylor–
Galerkin procedure and the second term adds to the stabilizing diffusion in the streamline
direction. Indeed, a similar form can be obtained here in a variety of ways, and recently Oñate
[33] introduced an interesting form of stability control by using so-called ‘finite increment’
equations. However, only the characteristic form ensures the optimal approximation, as
Equation (8) is derived from a self-adjoint form in which spatial discretization by the Galerkin
method can be optimally used. We can thus write the approximation
f = Nf( (9)
and use the weighting N T in the integrated residual expression. Thus, we obtain
M(f( n + 1 −f( n) =Dt[(Cf( n +Kf( n +f n) − Dt(Kuf( n + f ns )], (10)
if we assume u3 =0 as is generally done in an explicit form and omit higher derivatives and
source terms. In the above equation,
M=
& N TN dV, C=
& NT
(
(u N) dV,
(xi i
& &
V V
(N T (N
K= k dV, f= N TQ dV + b.t. (11)
V (xi (xi V
and Ku and f ns come from the new term introduced by the discretization along the characteris-
tics. After integration by parts, the expression of Ku and fs is
Ku = −
1 & (
(uiN T)
(
(u N) dV,
V (xi (xi i
&
2
1 (
fs = − (uiN T)Q dV +b.t., (12)
2 V (x i
where b.t denotes integrals along region boundaries [32].

The approximation is valid for any scalar convected quantity even if it is the velocity ui
component itself, as is the case with the momentum conservation equations. For this reason we
have elaborated the full details of the spatial approximation above as the matrices will be
repeatedly used.
It is of interest that the explicit form of Equation (10) is only conditionally stable. For
'
one-dimensional problems, the stability condition is given as (neglecting the effect of sources)
h 1 1 1
Dt 5Dtcrit = + − (13)
u Pe 2 3 Pe
for linear elements. In this equation, the Peclet number Pe is defined as
uh
Pe= . (14)
2k
In two-dimensional problems, the critical time step may be computed as [34]
Dts Dtn
Dtcrit = , (15)
Dts +Dtn
where Dts is given by Equation (13) and Dtn = h 2/2k is the diffusive limit for the critical
one-dimensional time step.
Further, with Dt =Dtcrit, the steady state solution results in a diffusion change almost
identical to that obtained using the optimal streamline upwinding procedures [32]. Thus, if
steady state solutions are the main objective of the computation, such a value of Dt should be
used in connection with the Ku term.
In Equation (8) we have introduced a u3 parameter which, when \ 12, stabilizes the explicit
solution. A fully implicit form is, however, expensive, involving non-linear, unsymmetric
matrices. However, it is often convenient to apply u3 ] 12 to the diffusive term only. We call this
a nearly implicit form and if employed we return to the stability condition
h
Dtcrit = , (16)
u
which can present an appreciable benefit.
3. THE GENERAL FRACTIONAL STEP ALGORITHM FOR THE NAVIER–STOKES

EQUATIONS
3.1. The equations of flow— Na6ier – Stokes problem

The full conservation form of the Navier–Stokes equations for compressible flow is
traditionally written as
(V (Fi (Gi
+ + +Q = 0, (17)
(t (xi (xi
with
V T = (r, ru1, ru2, rE) (18)
being the independent variable vector,
F Ti =(rui, rui u1, di1p, rui u2 +di2p, ui (rE + p)) (19)
defining the convective flux vector, and
G Ti = (0, −ti1, − ti2, qi −tij uj ) (20)
defining the diffusion fluxes. Finally,
Q T =(0, rg1, rg2, r(gi ui +q§)) (21)
gives the source terms.
In the above, the stress components tij are related to velocity gradients by
tij =m
(ui (uj 2 (uk
+ − d ,
(22)
(xj (xi 3 (xk ij
where ui are the velocity components; r is the density; E is the specific energy; q is the heat
flux; q§ is the heat generation per unit mass; p is the pressure; and gi is the acceleration due
to gravity.
The equations are completed by the universal gas law
p = rRT, (23)
where R is the gas constant.
The sound velocity is defined assuming constant entropy as
(p gp
c2= = , (24)
(r r
where g is the ratio of specific heats.
Further, we can write conveniently
(r (r (p 1 (p
= = , (25)
(t (p (t c 2 (t
though this expression assumes again constant entropy and is therefore only an approxima-
tion. In the following, we shall use Equation (25) but elsewhere we discuss the possibility of
correcting any errors involved by amendment of the algorithm [7].
While in the gas flow all the equations are fully coupled, for incompressible flows in which
c =, the energy equations can be solved independently after the velocity field has been
established. Nevertheless, a single algorithm for the solution of both problems is possible as we
shall now show.
Although the form of Equation (17) is identical to that of the convection–diffusion problem
of Equation (1), three wave speeds exist and the characteristic Galerkin procedure can not be
directly applied. In the next section we show how this can be employed with a fractional step.
In Appendix A, we show the depth averaged, shallow water equations. The format of these
is precisely the same as that of the Navier –Stokes equations and all of the procedures
discussed above can be again used for that class of problems.
3.2. Characteristic-based-split algorithm

For convenience we shall rewrite Equation (17) in a more direct form, omitting initially the
energy equation. These equations can be solved completely in a time increment Dt as the only
coupling that exists is through the speed of sound c for which we shall simply use the value at
time tn due to the explicit nature of the time stepping algorithm.
We thus write the first form of Equation (17), i.e. the mass flow continuity, as
(r 1 (p (Ui
= 2 =− , (26)
(t c (t (xi
in which we use Equation (25). Similarily, for each of the momentum conservation equations,
we write
(Ui ( (tij (p
=− (uj Ui ) + − −gi. (27)
(t (xj (xj (xi
In the above we define the mass flow fluxes as
Ui =rui (28)
and note that we can discretize Equation (27) in time using the characteristic Galerkin process
as, except for the pressure term, this equation is similar to the convection–diffusion 6iz.
n
Equation (8). Therefore, with the variable being Ui, we have
( (t Dt ( ( n
U ni + 1 − U ni = Dt (uj Ui ) + ij + uk (u U )+ gi
(xj (xj 2 (xk (xj j i
−(1 −u2)Dt
(p Dt
− uk
( (p n n
− u2Dt
n
(p n+1
. (29)
(xi 2 (xk (xi (xi
It is observed that the pressure term has now been evaluated at t n + u2Dt and is not treated
explicitly. Before proceeding further it is convenient to introduce an auxiliary variable U *i such
n
that
( (t Dt ( ( n
DU*i = U*i n + 1 −U ni =Dt − (u U ) + ij − gi uk (uj Ui )+ gi (30)
(xj j i (xj 2 (xk (xj
represents the first explicit part of the split and the ‘correction’ given below is available once
n
the pressure increment is evaluated.
DUi =U ni + 1 − U ni =DU*i −Dt(1 − u2)

(p Dt
− uk
( (p n n
− u2Dt
(p n+1
(xi 2 (xk (xi (xi
=DU*− Dt

(p n
+u2
(Dp n
+(1 −u2)
Dt 2
uk
( (p n
. (31)
(xi (xi (xk (xi
i
2
In the above, Dp =p n + 1 −p n and the last term, as before, represents the ‘source’ correction.
n
From Equation (26) we have, omitting third-order terms
1 n (U ni + u1 (U ni (DU*i ( 2p n ( 2Dp
Dr = Dp = −Dt = − Dt + u − Dtu + u .
(xi (xi (xi (xi (xi (xi (xi
1 1 2
c2
(32)
It is clear that the equations can be solved after spatial discretization in the following order:
Equation (30) to obtain DU*i ,
Equation (32) to obtain Dp,

Equation (31) to obtain DUi, thus establishing the values at t n + l.
It is important to remark that this sequence allows one to solve in an efficient manner and
with the adequate numerical damping of the continuous equations (26) and (27), i.e. the
conservation form of the Navier – Stokes equations. Therefore, this algorithm is well suited for
dealing with supersonic and hypersonic problems in which the conservation form ensures that
shocks will be placed at right position.
In all of the equations given below the standard Galerkin procedure is used with a
discretization
Ui =NU( i, DUi =NDU( i, DU*i =NDU( *i (33)
and
p= Npp̄. (34)
This gives from Equation (30) the solution for U( *i as
Step 1
DU( *i = − M − 1Dt[(CU( +KU( −f) −Dt(KuU( + f)]n, (35)
where all the discretization matrices are the same as those defined by Equations (11) and (12).
Similarily, the discretization of Equation (32) gives
Step 2
(M0 + Dt 2u1u2H)Dp̄ =Dt[Q(U( n +u1DU( *) −Dtu1Hp̄−fp ]n, (36)
which can be solved for Dp̄.
The new matrices arising here are
H=
& (N Tp (Np
dV,
V (xi (xi
M0 =
& N Tp
1 n
Np dV, (37)
c2
&
V
(N Tp
Q= N dV.
V (xi
The question of establishing the boundary conditions for the pressure is discussed in detail in
Appendix B.
The final stage of the computation of the mass flow vector U ni + 1 is completed by the
discretization of Equation (31) and we have now simply
Step 3

DU( = DU( * − M − 1Dt Q T(p̄ n +u2Dp̄) +
Dt n
Pp̄ n , (38)
2
where
P= (1− u2)
& (
(uiN Tp )
(N Tp
dV. (39)
V (xi (xi
At the completion of this stage, the values of U( n + 1 and p̄ n + 1 are fully determined but the
computation of energy rE n + 1 is needed so that new values of c n + 1, the speed of sound, can
be determined.
The last of Equation (17), i.e. the energy conservation equation, can be written as
((rE)
=−
(
(ujrE) +
(
k
(T
−
( (
(uj p)+ (tijuj ). (40)
(t (xj (xi (xi (xj (xi
Once again the equation is identical in form to that of the scalar problem of convection–
diffusion (Equation (6)) if we observe that p, Ui, etc., have been determined. Now the last term
can be evaluated at time (n + u3) for improved accuracy but in what follows we shall take
u3 = 0 for simplicity.
Using the characteristic Galerkin approximation of Equation (8) and discretizing as
rE= NEE( , (41)
we have
Step 4
DE( = −Dt[CE( n +KT n +f ne −Dt(KuE( n +f nes)], (42)
where E( contains the nodal values of rE and again the matrices are identical to those
previously obtained (assuming that rE and T can be suitably scaled in the conduction term).
Again, the forcing vectors can be appropriately defined as
fe =
& N TE
(
(u pt u ) dV +b.t.,
(xj j ij i
&
V
1 ( (
fes = − (uiN TE) (ui ptijui ) dV. (43)
2 V (x i (x j
It is of interest to observe that the process of Step 4 can be extended to include, in an

identical way, the equations describing the transport of such quantities as turbulence parame-
ters, chemical concentrations, etc., once the first essential Steps 1–3 have been completed.
At this point it is clear why we call our scheme ‘characteristic based split’. Observe that
Equation (29) has been discretized along the characteristic of the total derivative applied to the
momentum components Ui, which is not a characteristic of the whole problem due to the
presence of the pressure. However, this pressure does not appear in Equation (30) for the
intermediate momentum. Therefore, the characteristic of the total time derivative applied to
the momentum can in fact be thought of as the characteristic of the problem for the
intermediate momentum. It is in this sense that we call our scheme ‘characteristic based’.
4. SEMI-IMPLICIT AND EXPLICIT FORMS OF THE ALGORITHM
The algorithm described can be used in a semi-implicit form and indeed only in this form can
incompressible problems in which c = and M0 = 0 be solved. Taking
1
5 u1 51,
2
1
5 u2 51, (44)
2
the algorithm is conditionally stable. The permissible time step is governed by the critical step
of the explicit relation solved in Step 1 of the algorithm. This is the standard convection–
diffusion problem discussed in Section 2 and the same stability limits apply, reaching a value
close to
h
Dtcrit = (45)
u
for an inviscid fluid.

For slightly compressible or incompressible problems in which M0 is small or zero, the
semi-implicit form is efficient and it should be noted that the matrix H of Equation (36) does
not vary during the computation process and can be partially inverted.
In other semi-implicit forms, when compressibility exists, the question of the correctness of
the approximation of Equation (25) can be questionable, although an iterative correction can
always be applied.
However, it is possible to revert to a fully explicit form by putting u2 = 0 but still keeping
u1 in the interval presented in Equation (44). Now of course the critical step will be reduced
to the order of
h
(46)
c+ u
and this is indeed the same limit as that encountered in other explicit forms of Euler or
Navier–Stokes computational schemes currently used.
Further, Equations (35), (36), (38) and (42) can be solved simultaneously if the correction
step on the right-hand-side of Equation (36) is omitted. This is an additional approximation
and is not necessary but is introduced here to mimic artificial diffusions previously used with
the standard Galerkin form.
The use of the approximation of Equation (25) is not necessary in any fully explicit scheme
expected as the density increment is directly obtained if we note that
M0 Dp̄ = MDr̄, (47)
With the above simplifications we can return to the original equation (17) and using the
Galerkin approximation on this we can write directly
DV( = −M − 1Dt
& NT
(Fi (Gi
+
1
dV − Dt
& N TD dV
n n
(48)
V (xi (xi 2 V
omitting the source terms for clarity. The added diffusion terms D are defined below and have
to be integrated by parts in the usual manner.
Á (2 Â
2 p
Ã (xi (xi

Ãui ( ( (ujru1) + (p n Ã
Ã

Ã (x( i (x( i (x1
(p n Ã
Íui (xi (xi (ujru2) + (x2 Ì
Ã ( (
u
(u ru ) +
(p n Ã
Ã i (xi (xi j 3 (x3 Ã

Ãui ( ( (ujrE +uj p) n Ã.
Ä (xi (xi Å (49)
The ‘diffusions’ added are simple and are streamline-oriented (for divergence-free velocities),
thus not masking the true effects of viscosity as happens in some schemes.
If only steady state results only are sought it would appear that Dt in the definition of the
matrix D should be set at its optimal value of Dtcrit : h/u and we generally recommend this
value for the full scheme.
However, the oversimplified scheme of Equation (49) loses some accuracy and, even when
steady state is reached, will give slightly different results to those obtained using the full
sequential updating. The small additional cost involved in computing the sequence DU( *
Dp DU( DE( will have to be balanced against the accuracy increase. In general, we have
found that the two-step version is preferable.
It is of interest to note here in passing that the full sequential scheme introduces a so-called
‘forth-order’ diffusion proportional to Dt(Q TM − 1Q− H)p̄ in addition to the second-order
diffusion proportional to DtHp̄ into the computation. We shall indicate how this arises and
how keeping it is useful in the next section.
5. THE UNEXPECTED BONUS OF ‘CIRCUMVENTING’ THE BB RESTRICTIONS
We examine here the structure of equations reached in steady conditions. For simplicity we
shall consider only the Stokes form of governing equations in which the convective terms
disappear. Further, we shall take the fluid as incompressible and thus uncoupled from the
energy equations. Now, the three steps of Equations (35), (36) and (38) are written as
DU( * = − DtM − 1[KU n −f],
1
Dp̄ = H − 1[Q(U( n +u1DU*) − Dtu1Hp̄ n − fp ], (50)
Dtu1u2
DU( = DU( *− DtM − 1Q T(p̄n +u2Dp̄).
In the steady state we have Dp̄ =DU( =0 and eliminating DU( * we can write (dropping the
superscript n)
KU( +Q Tp̄= f (51)
from the first and third equations of (50) and
QU( + u1DtQM − 1Q Tp − Dtu1Hp −fp =0 (52)

from the second and third.
n! " ! "
We finally have a system that can be written in the form of
K QT U( f
= 1 , (53)
−Q Dtu 1[H −QM − 1Q T] p̄ f2
where f1 and f2 arise from the forcing terms.
This system has a non-zero diagonal that is proportional to Dt and which, as already
mentioned, is very similar to the forms suggested by other reasoning [32] to avoid matrix
singularity and ensure stability. Further, the system is always positive definite.
Also, it should be noticed that if the additional simplification introduced into Equations (48)
and (49) is made to avoid the sequential operations, the diagonal term becomes Dtu1H and we
still avoid BB conditions.
However, it can be easily verified that if the pressure gradient term is retained in Equation
(30) (which would seem to give a better approximation), the diagonal term of Equation (53) is
identically zero and the BB conditions in the full scheme can not be avoided.
6. SOME SOLUTION OF TYPICAL EXAMPLES
In this section we illustrate the applications and show the advantages gained by the use of the
CBS algorithm in various classes of problems.
We shall follow the classification of problems into various categories stated in the introduc-
tion, but some general remarks will be made as these arise. In all of the problems discussed in
this paper, the same computer coding was used and only the simplest element, i.e. the linear
triangular, is used in all examples.
The application of higher-order elements and elements of different nature is of course
applicable under certain circumstances but we shall not make detail comments on those here.
6.1. Fully explicit procedure, subsonic, transonic and supersonic flows

The use of a fully explicit procedure is preferred in aeronautical computations and the CBS
algorithm is well suited for this.
We have used in all of the present examples the simplex, linear, triangular elements in
two-dimensional form or simple tetrahedral elements in three-dimensional form. However, the
use of quadratic (or higher-order) elements is possible, especially when steady state conditions
are only sought and when the ‘lumped’ approximation of the mass matrix causes no significant
errors. Examples of the use of such elements are, however, not numerous, with some in [3].
Example 1
NACA 0012 aerofoil with zero angle of attack, M= 0.5, 1.2.
In this example we assess the performance of the CBS in a fully explicit form in subsonic
and supersonic situations.
In Figure 1(a) and (b), a typical graded mesh for the study of the NACA 0012 aerofoil is
shown. In Figure 1(c) – (e), the results of two Mach numbers (0.5 and 1.2) with the fully explicit
form of the multistep and single step forms of the algorithm are presented.
It is of interest to note that for the lower (subsonic) Mach number, almost no difference in
performance is observed, while in the supersonic range, the single step algorithm gives some
Figure 1. Inviscid flow past NACA 0012 aerofoil, M = 0.5, a =0°. (a) Unstructured mesh of 1824 elements and 969
nodes; (b) details of mesh near leading-edge; (c) convergence for M=0.5 with multi- and single step, fully explicit
form; (d) convergence for M = 1.2 for multi step scheme; (e) convergence for M= 1.2 for single step scheme.
overshoots and oscillations. The conclusion, confirmed by other solutions, is that it is always
safe and more accurate to perform the full algorithm procedure.
Figure 2 compares the performance of the new algorithm with the performance of the
previously used Taylor – Galerkin scheme.
It is noticed that the CBS algorithm improves the results dramatically near the stagnation
solution point without the use of any additional artificial diffusion (which is essential to get
any reasonable result using the Taylor – Galerkin scheme).
Indeed, it is rather surprising that even the single step procedure gives nearly the correct
pressures/densities at the stagnation without the addition of further stabilizations.
Figure 2. Subsonic inviscid flow past NACA 0012 aerofoil with a =0° and M = 0.5. (a) Density contours with T – G
scheme with Cs =0 (no additional viscosity); (b) density contours with T – G scheme with Cs =0.5 (with additional
viscosity); (c) density contours with CBS scheme with Cs = 0 (no additional viscosity); (d) comparison of density along
the stagnation line.
Example 2
In6iscid transonic flow o6er NACA 0012 aerofoil with 1.2 ° angle of attack, M= 0.85.
In this example a total lift is generated and in addition some shock developments occur.
With triangular elements and on an unstructured mesh some adaptivity is used to remesh the
problem, as in Figure 3(a) and (b), where the adapted mesh as well as the pressure contours
are shown.
In Figure 3(c), the pressure distribution resulting in a lifting force is given and the results are
compared with those obtained in an AGARD workshop [35].
Examples 3 and 4
Supersonic and hypersonic in6iscid flow past cylinder, M= 3, 10.
These two examples are designed to plot the performance of the CBS algorithm in
supersonic (or hypersonic) computation for which often very specialized schemes are used. Of
course, with the use of the CBS algorithm, an addition of a shock capturing viscosity is needed
and many suitable schemes are discussed in [36].
Figure 4 shows the mesh used and the Mach number contours obtained for supersonic
inviscid flow past a full cylinder. Here the inlet Mach number is 3. The Mach contour in
Figure 4(b) is the result obtained from the MUSCL scheme [5]. Figure 4(c) and (d) are
obtained using the CBS and different shock capturing viscosities. Using the appropriate shock
capturing viscosity along with the CBS algorithm gives an excellent comparison with the
special schemes [5]. Figure 5 shows the quantitative comparison of the present predictions with
those of the special schemes [5]. It is seen that the present scheme agrees well with the special
schemes except behind the cylinder, where all the schemes differ from each other slightly.
Figure 6 shows the mesh generated and the results obtained for hypersonic inviscid flow past
a half cylinder. The inlet Mach number here is 10. The sensitive zone with recirculation behind
the cylinder is well captured by this algorithm.
Examples 5 and 6
Viscous compressible flow past a plate (M = 3.0, Re= 1000) and an NACA 0012 aerofoil
(M= 0.95, Re= 5000).
In both examples fully explicit schemes are used to solve supersonic viscous flow problems.
In the first example (the Carter problem), the Mach number at the inflow is 3.0, the inlet
Reynolds number based on the length of the plate is 1000. The temperature of the plate is

assumed to be constant and equal to the stagnation temperature given as
g−1 2
Ts =T 1+ M . (54)
2

The temperature dependence of viscosity is accounted through Sutherland’s law
m Tr +S0 T 1.5
= , (55)
mr T+S0 Tr
where S0 is Sutherland’s constant and is equal to 198.6° Rankine.
A uniform mesh with 13045 nodes and 6680 elements is used in this computation. The
results obtained are shown in Figures 7 and 8. It is seen that the contours are smooth and
shocks are in the proper location. A little underprediction in the pressure distribution (Figure
8(a)) can be attributed to the unstructured uniform mesh used in the study. However, the
present results appear to be free from any oscillations, unlike Carter [37]. The outlet velocity
profile agrees excellently (Figure 8(b)) with the Carter [37] results except near the shock.
Figure 3. Inviscid flow past an NACA 0012 aerofoil with an angle of attack, a =1.2°, M= 0.85: (a) adapted mesh,
5512 elements and 2821 nodes; (b) pressure contours; (c) coefficient of pressure distribution along the upper and lower
surfaces.
Figure 4. Supersonic inviscid flow past a full cylinder, a = 0°, M =3.0. (a) Adapted mesh, 24979 elements and 12651
nodes; (b) Mach contours with MUSCL scheme; (c) Mach contours with CBS scheme (M and ML shock capturing);
(d) Mach contours with CBS scheme (anisotropic shock capturing).
Figure 5. Comparison of coefficient of pressure and Mach number along the middle in horizontal direction of the
domain.
In the second problem, an adaptively graded mesh is used to resolve the wake region of the
viscous, transonic flow past an NACA 0012 aerofoil. The total domain diameter is 25 times the
chord of the aerofoil considered. The inlet Reynolds number is equal to 5000. The adaptive
procedures are performed based on References [38,39]. Figure 9 and Plate 1 shows respectively
Figure 6. Hypersonic inviscid flow past a half cylinder, a =0°, M =10.0. Adapted mesh and solution obtained.
Figure 7. Supersonic viscous flow past a plate a = 0°, M = 3.0, Re= 1000: (a) Mach contours; (b) pressure contours;
(c) density contours.
Figure 8. Comparison of pressure distribution and outlet velocity profile with Carter [37]. (a) Pressure distribution
along the surface, (b) outlet velocity profile.
the mesh used and the density contours. It is seen that the instability predicted in this example
well matches the one predicted using more complicated adaptive procedures [40].
6.2. Semi-implicit procedure

Below, three examples of a semi-implicit nature are presented.
Example 1
High Reynolds number flow in a lid-dri6en ca6ity, M= 0.
The first example is a fully viscous incompressible flow (M= 0) in a lid-driven cavity. This
is a well-known test case used by many authors. A high-Reynolds number of 5000 is used in
this study. An adapted mesh taken from [41] is used to study the problem. The number of
nodes and elements used are much lower than the one used in the normal computations [4,6].
Figure 9. Adapted mesh for transonic, viscous flow past an NACA 0012 aerofoil a =0°, M =0.95, Re= 5000
(elements = 32522, nodes = 16388).
Figure 10. Viscous incompressible flow in a lid-driven cavity, Re =5000. (a) Adapted mesh, elements =1962,
nodes = 1034; (b) streamlines; (c) pressure contours.
Figure 11. Comparison of present velocity distribution with the benchmark solution of Ghia et al. [42].
The mesh used and the streamline and pressure contours are shown in Figure 10. The
recirculation near the corners are well predicted by the scheme and compares well with the
benchmark solution of Giha et al. [42]. The velocity distribution across the cavity is compared
with the benchmark solution in Figure 11. The agreement is excellent.
Example 2
Incompressible flow past a backward-facing step, M= 0.
This is another benchmark problem for which experimental results are available. Here also,
the adapted mesh is used to solve the flow. The inlet velocity profile is assumed to be parabolic
and the inlet Reynolds number is equal to 229. The total length of the domain is taken equal
to 40 times the step height so that traction-free conditions are valid at the outlet.
Figure 12 shows the mesh used and the solution obtained. The pressure distribution is
smooth and the comparison of velocity profiles with the experiments [43] is good. The small
deviation from the experimental prediction can be attributed to the inlet parabolic profile used
in the present study. The referred experimental results are generated from an inlet profile
slightly different from the parabolic profile.
Example 3
In6iscid transonic compressible flow past an NACA 0012 aerofoil, M= 0.5.
In this example, the inlet Mach number is assumed to be equal to 0.5. The velocity is
prescribed at the inlet and the density is prescribed at the exit. The density is calculated
analytically using the following relation
Figure 12. Incompressible viscous flow past a backward-facing step, Re = 229, finite element mesh and solution
obtained. (a) Mesh hmin = 0.05, hmax = 0.5, C = 0.15, number of nodes=1820, number of elements= 3441; (b)
comparison of velocity profiles; (c) streamline contours; (d) pressure contours; (e) velocity vector plot.
Figure 13. Subsonic inviscid flow past an NACA 0012 aerofoil, semi-implicit solution, a =0°, M =0.5. (a) Pressure
contours; (b) pressure distribution near stagnation point.
r0 =r 1+
g −1 2
M
1/(g − 1)
, (56)
2
which gives a value of 1.1297 and the present computation gives 1.1320 and the difference is
negligibly small.
Figure 13 shows the pressure contours predicted by the present scheme. The oscillation-free
pressure distribution near the stagnation point shows the excellent performance of the CBS
algorithm in its semi-implicit form.
6.3. Implicit procedure

An example of an implicit nature is presented.
Figure 14. Incompressible inviscid flow past an NACA 0012 aerofoil, implicit solution, a = 5°, M= 0. (a) Mesh near
stagnation point; (b) pressure contours.
Figure 15. The Severn Estuary. (a) Analysis domain; (b) finite element mesh; (c) bathymetry.
Example 1
In6iscid incompressible flow past an NACA 0012 aerofoil with an angle of attack, M= 0.
In this example a fully implicit form of the solution procedure is used to solve the inviscid
flow past an NACA 0012 aerofoil with an angle of attack equal to 5°. Here, although an
additional system of equations needs to be solved per time step, the time step may be taken
much larger, making the total CPU time needed much smaller. The mesh used near the
stagnation point and the pressure contours obtained are shown in Figure 14. The comparison
between the analytical stagnation and the present prediction is very close and the deviation is
only about 3%. Further details of the implicit procedure can be found elsewhere [7].
Figure 16. (a) Comparison of computed and measured surface elevation; (b) comparison of computed results for
explicit and semi-implicit schemes; (c) comparison of computed and measured elevations at two points in Severn Bore.
6.4. Shallow water problems

We have mentioned that for this class of problem, the first efficient use of the fully explicit
form was made by Peraire et al. [22]. For tidal problems, however, the semi-implicit form is
highly desirable and much larger time steps can be used.
Table I. Comparison of present results with experiments
Observations (HRS) Model (d.t. =50 s)
Station a (cm) g (d) a (cm) g (d) % (a)
Boundary (mid) 235 160 238 160 (0) 1.2

Tenby 262 170 254 170 (0) 3
Swansea 315 173 307 176 (+3) 2.5
Ilfracombe 308 162 300 165 (+5) 2.6
Barry 382 182 387 186 (+4) 1.3
Porthcawl 317 173 327 178 (+5) 3
Avonmouth 422 202 432 191 (−11) 2.3
In Figures 15 and 16 we show an analysis of tidal flows in the Severn Estuary and both
results of the explicit and semi-implicit procedures are given. Table I shows the comparison of
the predicted and experimental values.
Plate 2 shows a steady state supercritical flow in an open channel in colour. Here, just like
supersonic flows, a fully explicit scheme is desirable.
6.5. Solid dynamic problems

We conclude this section of examples by a somewhat unusual elesto-plastic solid and the
CBS scheme is used entirely to make use of its capabilities of overcoming the BB-type
restrictions. Such restrictions allowed so far only to use rectangles and hexahedra and a single
point integration to be used as triangles and tetrahedra were not effective.
The present scheme allows both triangles and tetrahedra to be used in such codes as DYNA
3D, as Figures 17 and 18 show the use of such elements in a dynamic problem.
7. CONCLUDING REMARKS
We hope that in this paper we have introduced the logical background and the excellent
performance of the CBS algorithm.
The development of the algorithm has taken some time and the present authors would like
to acknowledge the earlier works which pointed the way to the final algorithm. In our original
paper [1], we discuss the SEARCH process and mentioned the early papers that were only
partly successful. However, we borrow from the title of the first of these early ventures a title
that well describes the present procedure as
‘AN ALGORITHM FOR ALL SEASONS’
ACKNOWLEDGMENTS
This research has been partially supported by NASA grant NAGW/2127, AMES Control
Number 90-144.
Figure 17. Impact of plastic bar—axisymmetric solution. (a) Initial shape with the finite element discretization; (b)
linear triangles — CBS algorithm; (c) linear triangles—displacement algorithm; (d) bilinear Quadrilaterals— DYNA
algorithm.
Figure 18. Impact of a plastic bar—three-dimensional solution. (a) CBS algorithm; (b) standard displacement
algorithm.
APPENDIX A. SHALLOW WATER DEPTH-AVERAGED EQUATIONS
These equations reduce the problem to a two-dimensional one considering only averaged
velocities in the vertical direction.
If x1 and x2 are the horizontal co-ordinates, Figure A1 defines the basic variables.
The equations are identical to those of a two-dimensional compressible flow without the
energy equation. Thus, we can write
Figure A1. Shallow water formulation.
1 (p (Ui
+ =0,
c 2 (t (xi
(Ui ( (p
+ (uj Uk ) + +Qi =0 (i, j= 1, 2), (57)
(t (xj (xi
where
Ui =hui (58)
denotes the mass flow,
c 2 = gh (59)
is the wave speed, and
1
p= g(h 2 − H 2) (60)
2
corresponds to the pressure, with h, H and g defined in Figure A1.
The source term is defined below
(H u u h (pc
Qi = − g(h− H) +g i 2 i + ri −ti + . (61)
(xi C h r (xi
Here, the bed friction force is defined by the Chezy formula (with C being the Chezy constant,
r is the Coriolis force), with the vector r is defined as
r=
n
− fu 2
, (62)
fu 1
with f being the Coriolis friction. Finally, ti is the wind surface traction.
Once again sub and supersonic phenomena can occur as shown in the examples.
APPENDIX B. BOUNDARY CONDITIONS
In most of the problems we generally have ‘open’ boundaries limiting the extent of the analysis
domain. On such boundaries, various approximations are conventionally used. It is not our
purpose to discuss this type of boundary but the reader should refer to procedures used
conventionally used in the downstream region [8,44].
Our interest, however, is the boundary conditions on solid walls and we discuss these below.
B.1. Boundary conditions on solid walls

(a) On solid boundaries
un =0 normal velocity zero (63)
and either
ts =0 tangential traction zero
or
us =0 tangential velocity zero (64)
(b) When computing DU( * we do not impose any part at DU( * = 0, but leave values to be
determined. However, traction conditions on deviatoric components must now be applied to
this and appropriate values of ts and tn must be calculated on the whole boundary. This
requires, in viscous problems, the evaluation of tij.
(c) When computing Dr, we integrate by parts obtaining
& N TDr dV = −Dt

& NT
(
U( ni +u1DU( i − u1Dt
(p n + u2
dV
(xi (xi
& &
V V
( n (N T (p n + u2
= − Dt NT U( i +Dtu1 DU( *i − Dt dV
(xi V (xi (xi
&
V
(p n + u2
− Dtu1 N Tni DU( *i −Dt . (65)
G (xi
Here, ni is the outward drawn normal. The last term in the above equation is equal to zero,
as is the condition of Equation (63), Stokes flow

un =ni DU( i = ni DU( *i −Dt
(p =0 (66)
(xi
to the order of Dt 2. This point seems to have baffled some investigators who simply assume
(p
=0 (67)
(xi
on solid boundaries. Note that this is not exactly true.
REFERENCES
1. O.C. Zienkiewicz and R. Codina, ‘Search for a general fluid mechanics algorithm’, in D.A. Caughey and M.M.
Hafez (eds.), Frontiers of Computational Fluid Dynamics, Wiley, New York, 1995, pp. 101 – 113.
2. O.C. Zienkiewicz and R. Codina, ‘A general algorithm for compressible and incompressible flow, part I. The split
characteristic based scheme’, Int. J. Numer. Methods Fluids, 20, 869 – 885 (1995).
3. O.C. Zienkiewicz, B.V.K. Satya Sai, K. Morgan, R. Codina and M. Vázquez, ‘A general algorithm for
compressible and incompressible flow part II. Tests on the explicit form’, Int. J. Numer. Methods Fluids, 20,
887 – 913 (1995).
4. R. Codina, M. Vázquez and O.C. Zienkiewicz, ‘General algorithm for compressible and incompressible flows, part
III. A semi-implicit form’, Int. J. Numer. Methods Fluids, 27, 13 – 32 (1998).
5. B.V.K. Satya Sai, O.C. Zienkiewicz, M.T. Manzari, P.R.M. Lyra, and K. Morgan, ‘General purpose vs. special
algorithms for high speed flows with shocks’, Int. J. Numer. Methods Fluids, 27, 57 – 80 (1998).
6. O.C. Zienkiewicz, B.V.K. Satya Sai, K. Morgan and R. Codina, ‘Split characteristic based semi-implicit algorithm
for laminar/turbulent incompressible flows’, Int. J. Numer. Methods Fluids, 23, 1 – 23 (1996).
7. R. Codina, M. Vázquez and O.C. Zienkiewicz, ‘A fractional step method for the solution of compressible
Navier – Stokes equations’, in M. Hafez and K. Oshima (eds.), Computational Fluid Dynamics Re6iew 1998, Vol.
1, World Scientific, Singapore, 1998, pp. 331–347.
8. O.C. Zienkiewicz, J. Szmelter and J. Peraire, ‘Compressible and incompressible flow: an algorithm for all seasons’,
Comput. Method Appl. Mech. Eng., 78, 105–121 (1990).
9. O.C. Zienkiewicz, ‘Explicit or semi-explicit general algorithm for compressible and incompressible flows with
equal finite element interpolation’, Report 90 /5, Chalmers University of Technology, 1990.
10. O.C. Zienkiewicz, ‘Finite elements and computational fluid mechanics’, Metodos Numericos en Ingenieria, SEMNI
Congress, Gran Canaria, 1991, pp. 56–61.
11. O.C. Zienkiewicz and J. Wu, ‘Incompressibility without tears! How to avoid restrictions of mixed formulation’,
Int. J. Numer. Methods Eng., 32, 1184–1203 (1991).
12. O.C. Zienkiewicz and J. Wu, ‘A general explicit of semi-explicit algorithm for compressible and incompressible
flows’, Int. J. Numer. Methods Eng., 35, 457–479 (1992).
13. P.D. Lax and B. Wendroff, ‘Systems of conservation laws’, Comm. Pure. Appl. Math., 13, 217 – 237 (1960).
14. J. Donea, S. Giuliani, H. Laval and L. Quartapelle, ‘Time-accurate solution of advection – diffusion problems by
finite elements’, Comput. Methods Appl. Mech. Eng., 45, 123 – 145 (1984).
15. J. Donea, ‘A Taylor–Galerkin method for convective transport problems’, Int. J. Numer. Methods Eng., 20,
101 – 119 (1984).
16. R. Löhner, K. Morgan and O.C. Zienkiewicz, ‘The solution of non-linear hyperbolic equation systems by the
finite element method’, Int. J. Numer. Methods Fluids, 4, 1043 – 1063 (1984).
17. O.C. Zienkiewicz, R. Löhner, K. Morgan and J. Peraire, ‘High speed compressible flow and other advection
dominated problems of fluid mechanics’, in R.H. Gallagher, G.F. Carey, J.T. Oden and O.C. Zienkiewicz (eds.),
Finite Elements in Fluids, vol. 6, Chap. 2, Wiley, New York, 1985, pp. 41 – 88.
18. A.J. Chorin, ‘Numerical solution of Navier–Stokes equations’, Math. Comput., 22, 745 – 762 (1968).
19. A.J. Chorin, ‘On the convergence of discrete approximation to the Navier – Stokes equations’, Math. Comput., 23,
341 – 353 (1969).
20. G.L. Goudreau and J.O. Hallquist, ‘Recent developments in large scale Lagrangian hydrocodes’, Comput.
Methods Appl. Mech. Eng., 33, 725–757 (1982).
21. O.C. Zienkiewicz, J. Rojek, R.L. Taylor and M. Pastor, ‘Triangles and tetrahedra in explicit dynamic codes for
solids’, Int. J. Numer. Methods Eng., 43, 565–583 (1998).
22. J. Peraire, O.C. Zienkiewicz and K. Morgan, ‘Shallow water problems: a general explicit formulation’, Int. J.
Numer. Methods Eng., 22, 547–574 (1986).
23. O.C. Zienkiewicz and P. Ortiz, ‘A split-characteristic based finite element model for shallow water equations’, Int.
J. Numer. Methods Fluids, 20, 1061–1080 (1995).
24. O.C. Zienkiewicz and P. Ortiz, ‘An improved finite element model for shallow water problems’, in G.F. Carey
(ed.), Finite Element Modelling of En6ironmental Problems, Wiley, New York, 1996, pp. 61 – 84.
25. P. Ortiz and O.C. Zienkiewicz, ‘Tide and bore propagation in the Severn Estuary by a new fluid algorithm’, Proc.
Int. Conf. on Finite Elements Fluids — New Trends and Applications, Venezia, 1995, pp. 1543 – 1552.
26. R.A. Adey and C.A. Brebbia, ‘Finite element solution of effluent dispersion’, in C.A. Brebbia and J.J. Connor
(eds.), Numerical Methods in Fluid Mechanics, Pentech Press, 1974, pp. 325 – 354.
27. R.E. Ewing and T.F. Russell, ‘Multistep Galerkin methods along characteristics of convection – diffusion prob-
lems’, in R. Vichnevetsky and R.S. Stepleman (eds.), Ad6anced Computer Methods in PDEs IV, IMACS, Rutgers
University, 1981, pp. 28–36.
28. O. Pironneau, ‘On the transport diffusion algorithm and its application to the Navier – Stokes equation’, Numer.
Math., 38, 309 – 332 (1982).
29. M. Bercovier, O. Pironneau, Y. Harbani and E. Livne, ‘Characteristics and finite element methods applied to
equations of fluids’, in J.R. Whiteman (ed.), The Mathematics of Finite Elements and Applications, Vol. V,
Academic Press, New York, 1982, pp. 471–478.
30. M. Bercovier, O. Pironneau and V. Sastri, ‘Finite elements and characteristics for some parabolic – hyperbolic
problems’, Appl. Math. Model., 7, 89–96 (1983).
31. J.P. Benque, J.P. Gregoire, A. Hauguel and M. Maxant, ‘Application des methodes du decomposition aux calculs
numeriques an hydraulique industrielle’, INRIA 4th Coll. Methodes de Calcul Sci. et Techn., Versailles, 1983.
32. O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, Vol. 2, 4th edition, McGraw-Hill, New York,
1991.
33. P. Oñate, ‘Derivation of stabilized equations for numerical solution of advective – diffusive transport and fluid flow
problems’, Comput. Methods Appl. Mech. Eng., 151, 233 – 265 (1998).
34. R. Codina, ‘Stability analysis of forward Euler scheme for the convection – diffusion equation using the SUPG
formulation in space’, Int. J. Numer. Methods Eng., 36, 1445 – 1464 (1993).
35. T.H. Pulliam and J.T. Barton, ‘Euler computations of AGARD working group 07 aerofoil test cases’, AIAA 23rd
Aerospace Sciences Meeting, Reno, NV, 14–17 January 1985.
36. P. Nithiarasu, O.C. Zienkiewicz, B.V.K. Satya Sai, K. Morgan, R. Codina and M. Vázquez, ‘Shock capturing
viscosities for the general fluid mechanics algorithm’, Int. J. Numer. Methods Fluids, 28, 1325 – 1353 (1998).
37. J.E. Carter, ‘Numerical solutions of the Navier–Stokes equations for the supersonic laminar flow over a
two-dimensional compression corner’, NASA TR-R-385, 1972.
38. J. Peraire, M. Vahdati, K. Morgan and O.C. Zienkiewicz, ‘Adaptive remeshing for compressible flow computa-
tions’, J. Comput. Phys., 72, 449–466 (1987).
39. O.C. Zienkiewicz and J. Wu, ‘Automatic directional refinement in adaptive analysis of compressible flows’, Int. J.
Numer. Methods Eng., 37, 2189–2210 (1994).
40. H. Borouchaki, P.L. George, F. Hecht, P. Laug and E. Saltel, ‘Delaunay mesh generation governed by metric
specifications. Part I. Algorithms’, Finite Elem. Anal. Des., 25, 61 – 83 (1997).
41. P. Nithiarasu and O.C. Zienkiewicz, ‘Towards better choices of adaptive mesh generation for fluid mechanics
problems’, 6th ACME-UK Conference, 6–7 April 1998, pp. 129 – 132.
42. U. Ghia, K.N. Ghia and C.T. Shin, ‘High-Re solution for incompressible flow using the Navier – Stokes equations
and multigrid method’, J. Comput. Phys., 48, 387–411 (1982).
43. M.K. Denham and M.A. Patrick, ‘Laminar flow over a downstream facing step in a two-dimensional channel’,
Trans. Inst. Chem. Eng., 52, 361–367 (1974).
44. T.C. Papanastasiou, N. Malamataris and K. Ellwood, ‘A new outflow boundary condition’, Int. J. Numer.
Methods Fluids, 14, 587–608 (1992).
Plate 1. Density contours.
Plate 2. Symmetric channel of variable width, contours of h, Froude number = 2.5.
Copyright © 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 31(1) (1999)

Numerical Methods in Fluids - 1999 - Zienkiewicz - The Characteristic Based Split Procedure An Efficient and Accurate

Uploaded by

Copyright:

Available Formats

Numerical Methods in Fluids - 1999 - Zienkiewicz - The Characteristic Based Split Procedure An Efficient and Accurate

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numerical Methods in Fluids - 1999 - Zienkiewicz - The Characteristic Based Split Procedure An Efficient and Accurate

Uploaded by

Copyright:

Available Formats

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS

Int. J. Numer. Meth. Fluids 31: 359 – 392 (1999)

THE CHARACTERISTIC-BASED-SPLIT PROCEDURE: AN

O.C. ZIENKIEWICZa,*, P. NITHIARASUa, R. CODINAb, M. VÁZQUEZb AND

KEY WORDS: CBS algorithm; fluid problems; convective– diffusion-type

(a) Transonic and supersonic flows, fully explicit process

2. THE SCALAR CONVECTION – DIFFUSION PROBLEM AND THE

where b.t denotes integrals along region boundaries [32].

3. THE GENERAL FRACTIONAL STEP ALGORITHM FOR THE NAVIER–STOKES

3.1. The equations of flow— Na6ier – Stokes problem

3.2. Characteristic-based-split algorithm

DUi =U ni + 1 − U ni =DU*i −Dt(1 − u2)

(xi 2 (xk (xi (xi

Equation (32) to obtain Dp,

It is of interest to observe that the process of Step 4 can be extended to include, in an

4. SEMI-IMPLICIT AND EXPLICIT FORMS OF THE ALGORITHM

for an inviscid fluid.

M0 Dp̄ = MDr̄, (47)

5. THE UNEXPECTED BONUS OF ‘CIRCUMVENTING’ THE BB RESTRICTIONS

DU( * = − DtM − 1[KU n −f],

DU( = DU( *− DtM − 1Q T(p̄n +u2Dp̄).

KU( +Q Tp̄= f (51)

from the first and third equations of (50) and

QU( + u1DtQM − 1Q Tp − Dtu1Hp −fp =0 (52)

6. SOME SOLUTION OF TYPICAL EXAMPLES

6.1. Fully explicit procedure, subsonic, transonic and supersonic flows

6.2. Semi-implicit procedure

6.3. Implicit procedure

6.4. Shallow water problems

Table I. Comparison of present results with experiments

Observations (HRS) Model (d.t. =50 s)

Station a (cm) g (d) a (cm) g (d) % (a)

Boundary (mid) 235 160 238 160 (0) 1.2

6.5. Solid dynamic problems

APPENDIX A. SHALLOW WATER DEPTH-AVERAGED EQUATIONS

Figure A1. Shallow water formulation.

APPENDIX B. BOUNDARY CONDITIONS

B.1. Boundary conditions on solid walls

& N TDr dV = −Dt

on solid boundaries. Note that this is not exactly true.

Plate 2. Symmetric channel of variable width, contours of h, Froude number = 2.5.

You might also like