| Copyright © 2009. National Academy of Sciences. All rights reserved. Terms of Use and Privacy Statement |
Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 568
24th Symposium on Naval Hydrodynamics
Fukuoka, JAPAN, 8-13 July 2002
An Evaluation of Verification Procedures for CFD Applications
L. Epa (Instituto Superior Tecnico, Portugal)
M. Hoekstra (Maritime Research Institute, NetherIands)
Abstract
The study presented in this paper is a contribution
to the establishment of a reliable procedure for ver-
ification of numerical simulations of turbulent flow
around ship hulls. In order to study the behaviour of
the verification process under various circumstances,
we have selected six test cases ranging from a 2D flow
governed by a Laplace equation to the flow around
the KVLCC2 tanker at model and full scale Reynolds
number. For each case, solutions on much more
grids are available than strictly required by the rec-
ommended procedures. Two types of power series
expansions are applied to estimate the error, using un-
known and fixed exponents. A procedure to estimate
the uncertainty when the number of grids available is
greater than the minimum required for Richardson ex-
trapolation is also tested. The results indicate that the
main sources of the scatter observed in the data for
ship flows are the insufficient geometrical similarity
of the grids and the use of numerical interpolation to
obtain flow quantities at locations that do not coin-
cide with grid nodes. In three-grid verification of ship
flow computations conservative uncertainty estimates
are therefore in place.
~ Introduction
The need to prove the credibility of numerical pre-
dictions has led to a broad discussion on Verification
and Validation of CFD-results in several forums, like
the AIAA, [1], the ERCOFTAC, [2] and the ITTC
Resistance Committee, [3]. As a result, a signifi-
cant progress has been achieved in several respects,
but some difficulties remain which frustrate the firm
establishment of standard procedures for Verification
and Validation of CFD for complex turbulent flows.
It is now commonly accepted that the verification
of a numerical simulation is concerned with the nu-
merical error of the prediction, as opposed to vali-
dation which must assess the modeling error. How-
ever, once in the verification stage the error has been
quantified as E, say, the embarrassing question may
be asked: "How do you know?" [41. In other words,
a Verification Procedure should also give a recipe to
estimate the uncertainty in that error, so that the true
solution can be found within the error band with suf-
ficient confidence.
Verification may be based on various concepts, [5],
but in the majority of practical applications grid re-
finement studies, using Richardson Extrapolation, are
chosen. We have also tried single-grid a-posterior)
error estimates, like the ones discussed in [61. The
latter techniques have been developed in the context
of grid adaptivity and we have found their usefulness
in error quantification doubtful, t7]. Therefore, in this
paper we apply Richardson Extrapolation (henceforth
abbreviated as RE) as a basis for error quantification.
The quantification of the numerical error by RE as-
sumes a power series expansion of the error as a func-
tion of the typical cell size. The extrapolation to grid
cell size zero requires numerical solutions on several
grids, which should be geometrically similar (in or-
der for the idea of a typical cell size to make sense).
In practice, one has to accept that a finite number of
terms, usually only one, of the power series expansion
is sufficient to determine the error, i.e. the numerical
solutions must be in the so-called asymptotic range.
It is our experience, [9], that the data of a grid re-
finement study of a complex turbulent flow calcula-
tion do not exhibit exactly the assumed behaviour of
the error. Some scatter in the data may be produced
by several sources:
.
It may be hard to ensure that all the grids are
geometrically similar.
OCR for page 569
· Most how solvers use limiters for convection or
in the turbulence model implementations.
· Even the finest grids are outside the asymptotic
range.
.
. The flow quantities of interest may require post-
processing of the calculated data using numeri-
cal integration or interpolation techniques.
If only solutions are available on the minimum
number of grids, required for RE in any of its vari-
ants, it is extremely difficult to estimate the uncer-
tainty. For example, Roache suggests the use of a
safety factor, [5], the value of which has to be cho-
sen by experience. A more careful verification pro-
cedure would require at least one extra grid to check
the result obtained from the Richardson extrapolation.
The uncertainty in the error estimation could then be
based on the variability of the error estimation from
different sets of grids.
Continuing to add more grids, one can introduce
least squares approximation to achieve good estimates
of the error and its uncertainty, L101. The mean stan-
dard deviation of the fit may then be used to quantify
the uncertainty of the error estimate in addition to the
variability of the error found from several grid com-
binations.
In this paper we focus on Verification Procedures
for steady-flow CFD on structured grids. Our ulti-
mate goal is to be able to make with the least effort
a reliable estimate of the error and its uncertainty in
numerical simulations of the flow around ships. This
paper is a contribution to reaching that goal and it re-
ports work that we have done to find out where and
why current verification procedures fail.
Our basic RE-based procedures are:
.
.
A single-term power series expansion with an
unknown exponent, the observed order of accu-
racy, suggested by Roache, t51.
A multi-term power series expansion with preset
integer exponents, suggested by Oberkampf, [8].
In order to better understand the difficulties that
we have encountered in the application of verifica-
tion procedures to complex turbulent flows, we have
selected six test cases, ranging from a 2-D Laplace
problem in a square domain to the calculation of the
flow around a ship hull. The test cases and their pur-
pose are:
1. 2-D incompressible flow in a square domain. A
simple equation in a simple domain.
2. The Convection-Diffusion transport of a Passive
Scalar in a square domain. Application of differ-
ent order of discretization of the convective and
diffusion terms.
3. 2-D Lid-driven cavity incompressible flow. Non-
linear problem in a simple domain.
4. Laminar flow around a Wigley hull. 3-D test
case with a geometrically similar set of grids.
5. Turbulent flow around a Wigley hull. Introduc-
tion of turbulence models.
6. Turbulent flow around the KVLCC2 tanker.
Practical test case with a complex flow and a dif-
ficult geometry.
The main questions we would like to answer are:
.
What is the origin of the scatter observed in the
data.
· How does one deal with the uncertainty estima-
tion for a set of data which does not follow ex-
actly the assumed behaviour of the error.
In general, the flow quantities of interest in a com-
plex turbulent flow are not directly available from the
flow field. Numerical techniques for integration and
interpolation are usually required to compute resis-
tance coefficients or flow quantities at physical loca-
tions which do not coincide with grid nodes. In this
paper, we have also examined the influence of these
techniques by comparing the results obtained with nu-
merical techniques with different orders of approxi-
mation. These techniques are explored only for the
3-D test cases.
This paper is organized as follows: section 2
presents the tested Verification procedures with their
error and uncertainty estimators. The results of the six
test cases are presented and discussed in section 3 and
section 4 summarizes the conclusions of this study. A
more elaborate account of this work can be found in
[14].
2 Verification Procedures
2.1 Error estimation
All RE-based verification procedures assume the
error of a numerical prediction to be estimated by a
Taylor series expansion:
n
em = hi—ho = ~ Disc ' ( 1 )
j=1
OCR for page 570
where Hi is the numerical solution of any local or inte-
gral scalar quantity on a given grid (designated by the
subscript i), TO is the exact solution, aj are constants,
hi is a parameter which identifies the representative
grid cell size, and pj are exponents related to the or-
der of accuracy of the method, while n is the number
of terms in the expansion. Equation (1) has 1 + 2n
unknowns, and so an equal number of numerical so-
lutions is required for their determination.
The representation of the discretization error by a
Taylor series expansion has some conditions attached
to it, which may or may not be satisfied, like for ex-
ample the geometrical similarity of the different grids,
needed to allow a typical grid cell size to be used.
In the classic approach, [5l, equation (1) is applied
in the asymptotic range, i.e. only the leading term is
retained in equation (1~. In this case, one obtains:
Pi—To = ahP .
(2)
There are three unknowns in equation (21: ¢0, a
and p, where p is the observed order of accuracy.
Therefore, three grids in the asymptotic range are re-
quired to find TO and p. In a grid convergence study
with three suitable grids, the solution is convergent if
1. (02 - ¢~) X (~3 - ¢2) ~ O-
2. p>O.
The first condition states that the changes in ~ are
monotonous while the second implies that a finite
value is obtained for the grid of cell size zero.
An alternative approach, suggested by Oberkampf,
[8], is to assume that all the exponents in equation (1)
are integers. This changes the number of unknowns
to 1 + n, viz. ¢0 and aj. If one retains more than
one term of equation (1), it is possible to take into
account the effect of discretization schemes with dif-
ferent orders of accuracy. The definition of an asymp-
totic range is still required, because the power series
is truncated. Theoretically one would expect that in-
creasing the number of terms considered in equation
(1 ) alleviates the grid refinement requirements and en-
larges the asymptotic range. However, increasing n
has the obvious penalty that also the number of grids
to be handled increases.
In the present paper, we test the following altema-
tive Verification procedures:
· Method p. One-term Taylor series expansion
with p unknown. Requires at least 3 grids.
.
Method tl. One-term Taylor series expansion
with p known and set equal to the lowest order of
accuracy of the discretization schemes applied,
which is 2 in all the test cases. Requires at least
2 grids.
· Method t2. Two-term Taylor series expansion
with fixed exponents 2 and 3. Requires at least 3
grids.
.
Method t3. Three-term Taylor series expansion
with fixed exponents 2, 3 and 4. Requires at least
4 grids.
Aiming for the least amount of work in a Verifica-
tion Procedure one would prefer to use the minimum
number of grids required by the procedure. In our ex-
perience, applying a Verification Procedure with the
minimum number of grids required may easily lead
to erroneous results. If more grids are available differ-
ent sets of grids may be selected to check the extrap-
olated values. Also, with more grids it is natural to
consider the application of the verification procedure
in the least squares sense. We have tried this before,
[9l, [101, [11] and t12l, and will use it again in this
paper, if only to check the quality of the simpler pro-
cedures. The least squares root approach is based on
minimizing the function
ng nj \ 2
S(`po, (/j, pj) = /\ it pi—(¢~0 + ~ ajhiPi)) ~ (3)
\
where ng is the number of grids available. The mini-
mum of the function S is found by setting the deriva-
tives of (3) with respect to oO, pj and at; equal to zero.
This procedure has a straightforward implementation
in the four methods introduced above. It is obvious
that 5(40, Al, pi) = 0 when ng is equal to the num-
ber of grids required by the number of terms retained
in the Taylor series expansion. Therefore, the least
squares root approach demands at least one more grid
than the minimum required to obtain the unknowns
00, pj and ocj.
2.2 Numerical Determination of e(~)
2.2.1 One term expansion with p unknown
The determination of (O. p and a from (2) is per-
formed with a system of three non-linear equations,
where TO and or may be easily eliminated to obtain:
¢3—¢2 h2 P(~) - 1
_ |/ ~ =n
O2 - Of
(he ) 1
a. (4)
Equation (4) is a non-linear equation that defines
p and is solved using a false position method, with
an initial search of possible solutions in the interval
OCR for page 571
O < p < 8. If there are no possible solutions in this in-
terval, the solution is assumed to diverge for the given
triplet.
With p determined from (4), it is easy to obtain
en = ARE = ¢1 ¢° (~;)~ ~
where ARE is no more than the error estimate based
on Richardson extrapolation using the solutions ¢1,
2 and ¢3.
When the least squares method is applied, equation
(3) takes the form
S(¢O, a,p) = ;~ (¢i - (do + ahP)) (6)
Setting the derivatives of 5 with respect to TO, a
and p equal to zero leads to:
ng ng
of, ~ - a I, he
l I
i=1 i=1
Po = -
ng
and
ng ng
, (7)
ng ng \ ng
ng ~ Airy—~ (i | I, hP
i=1 i=1 / i=1
= Q
ng ng ~ ng v
ng I, hi2P—~ hP ~ ~ hP
i=l i=1 J i=l
Hi, ~ihPlog(hi) - TO ~ hPlog(hi) - a Ad, h, Plog(hi) = 0
i=1 i=1 i=1
(9)
Equation (9) is a non-linear equation that defines p
and is solved with a false position method. As in the
triplet approach, we have only considered solutions
with O < p < 8. The number of grids, ng, must be
greater than 3.
2.2.2 Expansions with fi Red exponents
In the present test cases, the lowest order term has
an exponent of 2. In any of the three expansions
tested, with one, two or three terms, the values of TO
and oci are obtained from a system of linear equations,
which may be written for the three terms expansion as
(10)
For the two terms expansion, (10) can be restricted
to a 3 x 3 system, while for the one term expansion a
further reduction to 2 x 2 is appropriate.
The application of the least squares root approach
to the power series expansion with fixed exponents
leads to a system of two, three or four linear equa-
tions, depending on the number of terms retained in
the series.
For the three terms expansion, pO, a2, a3 and a4
are obtained from:
ng ng ng
ng I, 02 I, 03 I, 04
i=1 i=1 i=1
ng ng ng ng
~hi2 hid EhS his
i=1 i=1 i=1 i=1
ng ng ng ng
At; 03 ~~ his >, 06 ~ h7
i=1 i=1 i=1 i=1
ng ng ng ng
04 ~~ 06 ~ h7 >, 08
i=1 i=1 i=1 i=1
ho
a2
a3
a4
_
=
~ ~ (i
i=1
ihi2
, i=1
, ~ ~ihi3
i=1
ng
~ ~ih'
i=1
(1 1)
The solution for the one and two terms expansions
is obtained from an appropriately reduced version of
the system (111. The minimum number of grids, ng,
should be 3, 4 or 5 according to the number of terms
retained in the expansion.
2.3 Uncertainty estimation
An uncertainty has to be defined for the estimation
of the error using Richardson extrapolation, because
some assumptions are implied in its application, like
the requirement that all the grids are in the so-called
asymptotic range.
The objective of the uncertainty definition is to
guarantee that the true error is banded by the error
estimation plus/minus an uncertainty U:
loo—0exact1 = |~1)i Exact ARE| < u (12)
In [5], Roache suggests by his Grid Convergence
Index (GCI) the use of a safety factor, Fs > 1, which
implies the uncertainty estimate1
U = (Fs—1) BERET (13)
A practical problem of this approach is that Fs has to
be chosen by experience. In test cases with known
exact solutions, one may re-write equations (12) and
(13) to obtain the safety factor required to bound the
error:
(FS)min 1 + 1 5RE, (14)
lIn the GCI method, [5], Roache defi nes GCI = U = R|6RE|'
which is completely equivalent to equation 13. The only difference
is that U is deli ned with ~ as the reference solution instead of ¢0.
OCR for page 572
Where possible, (Fs)min' will be evaluated and
compared with a number of alternatives for equation
(13) that we introduce here. We try to represent three
sources of uncertainty:
1. The uncertainty due to the assumption that all the
grids are in the asymptotic range, UO.
value for UO is likely to come out as much bigger than
for the other methods. This does not seem unreason-
able because the tl procedure uses only two grids.
To check the viability of this uncertainty definition
in cases where the analytical solution is known, we
have defined the parameter Fsu given by
2. The variability of ARE with the selection of dif- IXR I (20)
ferent sets of grids, Up.
3. The possible existence of scatter in the data, Us.
It is hard to make a proper guess for UO. In t13l,
Stern et al. have proposed a correction factor based
on the extrapolations with the observed order of accu-
racy and with the theoretical or asymptotic orders of
accuracy. However, the choice of the asymptotic or-
der of the method, may not be as easy as it seems, as
we will demonstrate in particular with the second test
case. Moreover, use of the correction factor is recom-
mended only if the asymptotic and observed orders
differ not too much, while the uncertainty vanishes
when the asymptotic order equals the observed or-
der, so that often recourse must be taken to Roache's
safety factor.
In the present study, we have defined UO as
where TO* is obtained as follows:
· In method p (one-term expansion with p un-
known)
Hi = To +'~hi + /32hi+i , (16)
where I = max(1, int(p)).
· Method If (one-term expansion with fixed expo-
nent)
(i = To +~3~hi3 . (17) \4
· Method t2 (two-term expansion with fixed expo-
nent)
Hi = To +~hi2+~32hi4 . (18)
· Method t3 (three-term expansion with fixed ex-
ponent)
pi = ho + Bohr + p2hi3 + )2hiS . (19)
UO may be computed for grids doublets, triplets or
quadruplets or with the least squares root approach.
For the one term expansion with fixed exponent, the
In the present study, we have in all the six exam-
ples a number of grids which is significantly larger
than the minimum number required to apply any of
the procedures described above. Therefore, the proce-
dures may be applied to different sets of grids leading
to different values for the error estimation. This al-
lows the computation of the range in [iiRE from all the
sets of grids selected. The outcome is clearly related
to the uncertainty in the error estimation. Therefore,
we define
Up = ~ ( ARE )max—( ORE )min ~ (21 )
The value of Up may be computed for any of the pro-
cedures described above.
Associated with UO and Up, we have computed the
parameter FsUp, which is a kind of safety factor de-
fined by
(15) Fsup = 1 + °~` j '' . (22)
The main purpose of the least squares root ap-
proach is to deal with data which show a certain
amount of scatter. This scatter is also a source of un-
certainty in the error estimation. Therefore, we have
computed the standard mean deviation of the fit to try
to quantify the uncertainty due to scatter. Accord-
ingly, we have defined a contribution to the uncer-
tainty, Us' given by
ng ~ n; \ 2
Mini—(¢o+~jhP~)|
Us=\ i=l \ j=l J , (23)
ng—nu
where nu stands for the number of unknowns to be
determined by the fit.
When the least squares root approach is applied, Us
is added to UO and Up to obtain the uncertainty. This
approach may be checked in test cases with known
analytical solutions using the parameter FsUps given
by
F 1 + UO + Up + Us (24)
If the uncertainty computed from UO, Up and Us
works well, the parameters Fsu Fsup and FsUps should
be larger than (Fs)min. Furthermore, these parameters
should increase for coarser grids and tend to (Fs)min
with the grid refinement.
OCR for page 573
3 Verification tests
In the present study we have tested the four verifi-
cation procedures introduced in the previous section
as methods p, tl, t2 and t3. These four possibilities
were explored with the use of the least square root ap-
proach and with different sets of grid doublets, for tl,
triplets for p and t2 and quadruplets for t3.
In all the test cases, we have performed grid de-
pendency studies with a number of grids larger than
the minimum required by any of these techniques. In
these conditions, the number of grid sets to which the
error estimation may be applied is rather large. There-
fore we have limited this number by imposing the fol-
lowing cnteria:
In every case, the most refined grid, which is iden-
tified by hi, is always included in the error estimation.
In method p the three grids required, designated by
hi, hi and 02, have to satisfy the conditions:
(¢i2-¢i) (I-) >0 , ~ <2 , ~ <4
0.5<~<2.
(25)
In method tl, only two grids, designated by h, and
hi, are required. The selected grid doublets have to
satisfy the conditions:
1.5< i<2.
hi
(26)
In method t2, the required grid triplets, identified
by hi, hi and 02, have to satisfy the conditions:
i<2 i2~4. , o.5< i2/ i<2
hi - hi hi/h1 (27)
In method method t3, the required grid quadru-
plets, designated by hi, hi, 02 and 03, have to satisfy
the conditions:
<<2, t<2,
h /h
0.5<~<2 0~<2 (28)
We have applied the least squares approach to sets
of grids ranging from minimum number+ 1 to ng. The
uncertainty, Up, is estimated including the arid set
with the finest grids.
3.1 2-D Incompressible Potential Flow
A well-known result of fluid dynamics is that the
solution of a 2-D incompressible potential flow may
be obtained by solving the Laplace equation under
suitable boundary conditions with the velocity poten-
tial or the stream-function, fir, as the dependent vari-
able.
A representative flow was created by superposing a
uniform flow and eight line-vortices2 located outside
the computational domain, which was chosen as the
unit square O < x < 1 A O < y < 1. Dirichlet boundary
conditions were applied. The same test case has been
used in previous verification studies, [12], L11] and
[71.
Three sets of grids were employed: equally-spaced
rectangular gnds, non-uniform rectangular grids and
non-orthogonal grids. In each case, we generated 24
grids with an equal number of nodes in each direc-
tion with a coarsest grid of 11 x 11 grid nodes and a
finest grid of 241 x 241 grid nodes. These 24 grids
have 81 common grid nodes in the interior of the do-
main. Therefore, there is no need to use any inter-
polation technique to apply the different Verification
procedures for local variables.
All the calculations were performed with 15 dig-
its precision and the solution of the linear system
of equations was carried out to computer precision.
Therefore, iterative errors are negligible compared to
discretization errors.
16.92
16.915~ ) 1 2 3 4
hi/hi
r38lG~0°°O
a
0 Numeric
o Exact
a
a
to
to
to
Figure 1: to at (x= 0.3,y = 0.9) in the 2-D incom-
pressible potential flow.
O O _
Extensive results have been reported in [14], com-
paring the use of the four different power series ex-
2In 2-D9 a line-vortex is represented by a single point.
OCR for page 574
OCR for page 576
OCR for page 577
OCR for page 579
pensions and the use of grid doublets, triplets or
quadruplets versus the least squares root approach. As
an example, figure 1 presents the numerical solution
for the equally-spaced grids at (x = 0.3,y = 0.9~.
Here we summarize the main conclusions of this
exercise as follows:
.
.
The observed order of accuracy, p, deviates from
the theoretical value of 2 by less than 0.04 for all
grid triplets not including the grids coarser than
101 x 101.
Neither the stretching, nor the non-orthogonality
of the grid has a notable influence on the ob-
served order of accuracy.
· UO is not enough to bound the error in many of
the grid doublets, triplets or quadruplets tested.
However, the error is always bounded by UO + Up
in all the cases tested. The data showing minimal
scatter, Us gives the smallest contribution to the
uncertainty and is negligible compared to UO and
Up in any of the four methods tested.
· The least squares root method seems to be a
more robust procedure to estimate the error and
its uncertainty than the use of grid doublets,
triplets or quadruplets. The maximum ratio be-
tween FsUp and the minimum allowable safety
factor, (Fs~mi,~, is closest to one in the least
squares approach.
The safety factors computed tend to 1 with the
grid refinement. In this test case, FsUp is always
larger than (Fs~min, which means that UO + Up is
a safe estimate of the uncertainty. As an exam-
ple, figure 2 presents the values of (Fs~min, FsU,
FsUp and FsUps as a function of the typical grid
cell size of the coarsest grid selected for the p
approach.
3.2 Convection-Diffusion Transport of a
Passive Scalar
with
The second test case is related to the convection-
diffusion transport of a scalar in a square domain of
side L and was taken from E151. Using dimensionless
variables with a reference length of L and a reference
velocity U. the problem is described by
V ~ u ¢)—V ~ ~ V¢) = ~ ¢, (29)
pUL
Rn =
1.15
1.1
.
v
o
(Fs)min
F
su
F
sup
F
sups
Q
v
v
v
v
1
hi/h1
v
Figure 2: Estimation of the safety factors as a function
of the typical cell size of the coarsest grid selected
for or at (x = 0.3,y = 0.9) in the 2-D incompressible
potential flow. Least squares root approach with one
term and an unknown exponent.
The analytical solution is given by:
¢ (x,y) = sin ( 72c ~ ) sin ( 2 r ) (30
with O
equation (29) balances diffusion. This means
that when the Reynolds number tends to zero, the
convection-diffusion problem tends to a Poisson
problem. On the other hand, when the Reynolds
numbers tends to infinity, equation (29) tends to
convection equals zero.
The same set of grids as for the 2-D incompress-
ible potential flow is adopted in this case, and also
here Dirichlet boundary conditions are used at the
four boundaries. An extra layer of virtual grid nodes
is added to the domain to keep the third-order dis-
cretization scheme in all the interior nodes. In these
extra grid nodes, we obtain ~ from the analytical so-
lution.
In order to observe the influence of the Reynolds
number on the observed order of accuracy of the
method, calculations were made at Rn = 1, 103, 104
and 106.
As for the previous example, [141 presents a de-
tailed description of the data for the three sets of grids
and the four Reynolds numbers tested. In the present
paper we restrict ourselves to the main trends found
in the results.
.
In the equally-spaced Cartesian grids, the ob-
served order of accuracy is clearly dependent on
the Reynolds number. At the lowest Reynolds
number, Rn= 1, p is 2 and at the highest
Reynolds number, Rn= 106, the asymptotic or-
der of accuracy is 3. At Rn= 103, p is even
larger than 3, but in this case the variation in p
with the selection of the set of grids is the largest
of the four Reynolds numbers. Table 1 presents
the minimum and maximum observed order of
accuracy estimated from different grid triplets at
the locations that exhibit the maximum error in
the finest grid.
Rn 1 103 104 1 o6
P1 P2 P2 P2
Pmin 1.99 3.18 3.01 2.99
Oman 2.00 3.48 3.03 3.00
Table 1: Minimum and maximum orders of ac-
curacy estimated from different grid triplets in the
convection-diffusion transport of a passive scalar.
Equally-spaced Cartesian grids. PI = (x= 0.6,y =
0.6) and P2 = (x = 0.9,y = 0.4~.
· The observed order of accuracy is independent
of the Reynolds number and equal to 2 for the
stretched and non-orthogonal grids. This does
not mean that the stretching introduces a loss of
accuracy of the method. This result is simply ex-
plained by the numerical approximation of the
metric terms. In the equally-spaced grids the co-
ordinate derivatives are exact, but in a stretched
or non-orthogonal grid they are only second or-
der accurate. Therefore, at the highest Reynolds
numbers the approximation of the metric quan-
tities determines the order of the approximation
of the convective terms.
· In the equally-spaced Cartesian grids, the error
estimation with the one term expansion and a
fixed exponent, tl, does not perform well for the
three largest Reynolds numbers. In these three
cases, the observed order of accuracy is greater
or equal than 3 and the assumed exponent is 2.
The value of (Fs~mi,' does not tend to one with
the grid refinement. However, (Fs~min tends to
one in the other three methods with the grid re-
finement.
· The t2 and t3 series expansion exhibit the ex-
pected behaviour for all the Reynolds numbers.
As in the previous test case, the inclusion of ex-
tra terms in the series expansion improves the fit
to the data.
3.3 Lid-Driven Cavity Flow
The lid-driven cavity flow is a standard test case of
incompressible flow solvers. There are several numer-
ical solutions published for this problem in the open
literature, which are based on finite-difference, finite-
volume, finite elements or spectral methods. Unfor-
tunately, there is no analytical solution available for
such flow.
The flow domain has a very simple geometry, a
simple square, but the flow, even if laminar, has a
complex structure, that may include several vortices
depending on the Reynolds number. It is a completely
bounded flow, for which the no-slip condition defines
all the required boundary conditions. Three of the
four boundaries are fixed walls, while the remaining
boundary (the lid) drives the flow by moving with a
prescribed velocity.
The momentum equations at the boundaries reduce
to
1 ~ 02ui
Rn \< dL~2 J (32)
where ~ = y and i = 2 for the boundaries y = 0 and
y = 1, while ~ = x and i = 1 for the boundaries x = 0
and x = 1. In the present calculation, we have adopted
instead the approximation
UP
b; =0
(33)
to determine the pressure on the boundary. In the nu-
merical solution a second order three-point forward or
backward scheme is applied to equation (33~.
The boundary conditions refer to the velocity com-
ponents only. The impermeability condition requires
that the velocity component normal to the wall is zero.
The no slip condition defines the tangential velocity
component at the boundaries. So the pressure is de-
termined up to a constant, because the equations con-
tain only gradients of the pressure. In order to obtain
a unique solution, the pressure must be fixed at a grid
node. Theoretically the location where the pressure is
fixed is arbitrary. However, in the discretized finite-
difference system the numerical solution may be af-
fected by this choice. In the present calculations, the
pressure is fixed at the nearest node to the top left cor-
ner of the domain, (x = 0,y = 1), on the left vertical
boundary, (x= 0~. However, for the presentation of
the results we have made P= 0 at the centre of the
cavity, (x = 0.5,y = 0.5~.
Details of the numerical solution are given in L12].
In this paper, we only include the list of the discretiza-
tion schemes applied to show the different orders of
accuracy of the schemes adopted.
· Continuity Equation
- Third order four-point scheme with a fixed
bias for both directions.
· Momentum Equations
- Convective Terms
* Upwind third order four-point scheme
for both directions.
-
Pressure Gradient
* Third order four-point scheme with a
fixed bias for both directions.
- Diffusion Terms
* Second order central scheme.
If one assumes that the theoretical order of ac-
curacy is determined by the smallest order of accu-
racy of the discretization schemes applied, the present
method has P'h = 2.
In order to investigate the influence of the Reynolds
number, we have selected two Reynolds numbers to
perform Verification studies: Rn= 100 and Rn=
1000. The Verification studies are performed in
equally-spaced Cartesian grids. For each study, 16
grids with an equal number of nodes in each direction
were generated. The coarsest grid has 11 grid nodes
per direction and the finest grid has 161 x 161 grid
nodes. As in the previous test cases, there are 81 grid
nodes which are common to all 16 grids, which avoid
the use of interpolation techniques in the Verification
procedure. The variables selected for the Verification
studies are the two Cartesian components of the ve-
locity, Ui and U2, the pressure, P and the vorticity,
which in this flow has only one component given by:
SU~ _ ark (34)
by fix
The derivatives in equation (34) are approximated by
second order central-difference schemes.
The calculations were performed with 15 digits
precision and the iterative solution was carried out
until the residual of all the equations has dropped 13
orders of magnitude, from an initial solution of zero
values for all the unknowns. Therefore, iterative er-
rors are negligible compared to discretization errors.
The three power series expansions with fixed expo-
nents do not include the first order term. However,
in this test case, it is possible to obtain an observed
order of accuracy below the theoretical value of 2.
Therefore, when p is below 2, we have also tested
the same fixed exponents power series expansions in-
cluding the first order term. We will designate these
expansions by tl 1, t21, and t31 .
The details of all the comparisons performed for
the two Reynolds numbers tested are given in [141. In
this paper, we summarize the main results obtained
for the lid-driven cavity flow.
.
There is no scatter in the data.
· The asymptotic order of accuracy is more dif-
ficult to establish than in the previous two test
cases. Nevertheless, for the present set of grids,
it is possible to identify the following trends:
· At a given location, the value of p depends
on the flow variable considered and on the
Reynolds number
- For a given flow variable, the observed or-
der of accuracy depends on the location se-
lected.
- The value of p may be larger or smaller
than the theoretical order of 2.
· The comparison of the least squares root ap-
proach with the use of grid doublets, triplets or
quadruplets shows that for the same range of
grids the values of Up and the maximum changes
cp 1 U1 1 UZ 1 ~ 1
Rn = 100 |
2.05 2.01 2.02 0.95
2.31 2.39 2.41 0.99
0.15 2.73 2.16 1.03
1 42 2.96 2.62 1.06
Rn = 1000
P1 Pmin
P1 Pmax
P2 Pmin
P2 Pmauc _
= ~
P2 | Pmin
P2 | pmax
1.77 2.2
2.23 2.61 1 2~42 1 7~66 1
Table 2: Minimum and maximum orders of accuracy
estimated from different grid triplets for the lid-driven
cavity flow. P1 = (x= 0.2,y = 0.2) and P2 = (x=
0.8,y = 0.8).
in p are similar. Nevertheless, the effect of the
use of coarsest grids is easier to identify in the
least square root approach, using the standard
mean deviation of the fit.
.
The performance of the different alternatives of
power series expansion is clearly dependent on
the observed order of accuracy of the flow quan-
tity selected. The results obtained for the present
flow quantities show the following trends:
- For values of p close to the theoretical
value of 2, all the four alternatives are
in agreement, exhibiting the expected be-
haviour of the uncertainty increasing with
the use of coarser grids. In these cases the
p and t2 expansions lead to very similar
results, whereas the tl alternative exhibits
the largest values of the uncertainty for the
finest grids. The t3 approach shows the
weakest dependency of the uncertainty on
the coarsening of the grid.
- For values of p larger than 2, the tl ex-
pansion is not in agreement with the other
three alternatives, which all produce sim-
ilar results. Apparently, the extrapolation
performed with one term and the theoreti-
cal order of accuracy is not correct. As for
p close to 2, the results of the p and t2 ap-
proaches are very similar.
· For values of p smaller than 2, the expan-
sions with fixed exponents are clearly dif-
ferent from the one with an unknown expo-
nent. The results obtained by including the
first order term show a much better agree-
ment with the p approach. In this case,
the tl ~ expansion that includes two terms
gives similar results to the expansion with
unknown exponent.
-0.0209
-0.021. 2
, ;
o
o
a
() 5
C
o
A
O
A
I , , , , I , , , · I
3 4
hill
Figure 3: Cp at (x = 0.8,y = 0.8) as a function of
the typical cell size for the lid-driven cavity flow at
Rn= 100.
As an example of the results obtained in this test
case, figures 3 and 4 present the values of Cp obtained
at (x = 0.8,y = 0.8 ) for Rn = 100 and the extrapolated
cell size zero values with the band of uncertainty as a
function of the typical cell size of the coarsest grid
used in the least squares root approach for the p, tl
andtll approaches.
3.4 Flow around the Wigley Hull
We proceed with a numerical verification for a 3-D
flow, viz. the flow around the Wigley hull in laminar
and turbulent flow. The Reynolds number based on
the undisturbed velocity UOO and the ship length, L, is
7.4 x 103 for the laminar flow and 7.4 x 106 for the
turbulent flow. This test case has the advantage of
referring to a ship-like form defined by an analytical
expression:
Y 2 ILL ( L)1 ~ (H) 1 ' ( )
J L
with x and z in the range
and
O
-0.0208
-0.0209
v
.
~~V
-0.021
-0.0211
-0.0208
-0.0209
-0.021
a
v v
.
v
v
a
o
to~Uo
to+UO
to~(Uo+up)
(o+(Uo+up)
(o-(Ho+up+us)
o to+(up+Up+Us), . . .
2 hill 3 4
1
tl expansion
.
v
v
gs ~ ~ v
(O-UO
(O+UO
0 0
tO-(UO+Up)
(O+(Uo+up)
~o-(uo+up+us)
O tO+(Up I P L S)
-0.0211: 2 3 4
hint
a
l
-0.0208
-0.0209
-0.021
-0.0211
tl1 expansion
.
v
v
~ Q
~~ ~ ~ V
83~@
o
to~Uo
bo+Uo
bo-(uo+up)
to+(Uo+Up)
to-(Uo+Up+Us)
to+(Uo+Up+Us)
a
a
v
.
v
V
.
v
I I I ~ 1
hilt 3
4
Figure 4: Error extrapolation and uncertainty bands as
a function of the typical cell size of the coarsest grid
for Cp at (x = 0.8,y = 0.8) for the lid-driven cavity
flow at Rn = 100. Least squares root approach.
domain extending in lengthwise direction from x =
—0.15L to x = l.l5L for the laminar flow case and
x = - 0.25L to x = 0.25L for the turbulent flow case.
The domain is bounded externally by the relevant part
of the circular cylinder given by:
y2 + Z2 = (0.0937sLy2 -
In the two sets of grids of the laminar case, the
number of nodes varies from 241 x 91 x 61 in the
finest grid to 65 x 25 x 17 in the coarsest grid. For
the turbulent flow, the finest grid has the same number
of grid nodes in the streamwise and girthwise direc-
tions, but it has 121 grid nodes in the normal direc-
tion, which implies 33 in the coarsest grid to preserve
geometrical similarity.
The four sets of grids were generated with the same
technique. The two 161 x 61 x 41 (laminar case) and
161 x 81 x 41 (turbulent case) grids were generated
with a proprietary elliptic PDE grid generator, based
on the GRAPE approach with the required stretch-
ing in the normal direction applied algebraicly, [161.
The remaining grids were all generated with 3-D cu-
bic spline interpolations based on the 161 x 61 x 41
and 161 x 81 x 41 grids. Therefore, all the grids in
each set are geometrically similar. In general, this
procedure may be difficult to apply, but in the simple
geometry of the Wigley hull it is straightforward.
The two sets of grids generated for each Reynolds
number differ only in the streamwise grid line spac-
ing. In set A a weak stretching is applied at the "lead-
ing" and "trailing" edges of the ship, whereas set B
has a stronger stretching at these edges, hence more
variation in grid spacing. The 81 x 31 x 21 grids for
the laminar case at the boundaries are illustrated in
figure 5.
At the ship surface, the no-slip condition is applied.
The boundary conditions imposed on the external and
inlet boundaries were derived from a potential flow
solution. Symmetry boundary conditions are applied
at the free surfaces and at the ship symmetry plane. At
the outlet boundary, the streamwise pressure gradient
is set equal to zero.
All the calculations were performed with PAR-
NASSOS, [17], and were carried out to machine ac-
curacy so that iterative errors may be neglected when
compared to discretization errors. A system of curvi-
linear coordinates ( , ~7, A) is adopted with ~ roughly
aligned with the streamwise direction, rl in the normal
direction to the ship surface and ~ in the girthwise di-
rection.
To illustrate the results of the Verification exercise
we have selected integral and local quantities. The
integral quantities are:
3 Double model flav conditions were assumed.
xTy
Set A
, .
Set B
Figure 5: Illustration of the 81 x 31 x 21 grids for the
calculation of the laminar flow around the Wigley hull
at the boundaries.
· The friction resistance, CF.
Jo= (fir) ~a: xa; dude
C —
F— 2PU2L2
-, (36)
where a:, and a; are the co-variant base vectors.
· The pressure resistance, Cp.
CP = 2 i Cp (a: x a';) · nxd~d`, (37)
· The total resistance, C' = CF + CP.
The integrals are evaluated with Gaussian rules us-
ing two options for the integrand approximation:
1. A first order approximation, which assumes a bi-
linear variation in each grid cell. This option will
be designated be Linear.
2.
. A bi-linear interpolation of a bi-quadratic vari-
ation, which guarantees a third-order interpola-
tion scheme with continuity of the function and
of its first order derivatives. We will refer to this
option as Cubic.
We have selected flow quantities at locations which
always coincide with grid nodes and at selected coor-
dinates that require interpolation from the grid nodes.
Two first order and two third order interpolation
schemes were tested:
Firsts scheme that performs a linear interpolation
in the x direction followed by a bi-linear interpo-
lation in the y—z plane.
The Firstg scheme that is based on a tri-linear
interpolation in the 3-D space.
Thirds scheme which also splits the interpolation
in the x direction and in the y—z plane. The two
interpolations are performed with the blending
of first and second order basis functions.
· The Thirds scheme that is based on 3-D cubic
splines.
As in the previous examples, a complete discussion
of all the data is presented in [141. The Verification
studies for laminar and turbulent flows lead to results
which are qualitatively equivalent. Therefore, we will
restrict ourselves to the discussion of the main results
obtained in the turbulent flow case.
In locations coinciding with a grid node which is
common to all grids there is no scatter in the data.
Nevertheless, for some flow quantities the observed
order of accuracy may be difficult to establish because
the data do not always show a monotonous change.
The selection of the grid node distribution has a
strong influence on attaining the asymptotic range. As
an example, figure 6 presents the pressure coefficient
at the 'trailing edge' of the ship on the free surface,
(x = L,y = 0,z = 0), as function of the typical cell
size for the two sets of grids. The data of set A do
not exhibit a monotonous change, as opposed to the
data of set B which reveal a gradual decrease of Cp
with grid refinement. The local grid line spacing of
the finest grids of set A is larger than the grid line
spacing of most of the grids of set B. Nevertheless, it
is interesting to note that if the three finest grids of set
A are disregarded, one might be tempted to perform a
misleading extrapolation of CP to cell size zero.
The three resistance coefficients, CF, CP and Car,
are plotted as a function of the typical cell size for
sets A and B in figure 7. The variation of the inte-
gral parameters with the typical cell size does not ex-
hibit scatter. The results obtained with the Linear and
Cubic approaches are graphically coincident with the
exception of Cp for set A. The friction resistance, CF,
is clearly influenced by the choice of the streamwise
grid line spacing.
The estimation of the error and its uncertainty is
clearly more difficult than in the previous test cases.
Table 3 summarizes the minimum and maximum val-
ues of p obtained from grid triplets in the two sets
of grids. For Cp and Car of set B. there are no grid
triplets that satisfy the convergence criteria defined
for a given grid triplet.