Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
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.
· 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 = hiho = ~ Disc ' ( 1 ) j=1
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: PiTo = 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
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: loo0exact1 = |~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 = (Fs1) 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.
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) ngnu 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.
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.
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<X<1 AO<~<1 L. - L. and IT = sin (:z) cos (~) IJ = cos ~ 2 I; ) sin (2 {) (31) where us and u2 are the Cartesian velocity compo- nents in the x and y direction. Equation (29) is discretized in a general non- orthogonal curvilinear coordinate system with second-order central differencing schemes for the metric terms and for diffusion, which are equivalent to the ones applied in the Laplace problem, and with third-order upwind schemes in the convective term. Compared with the Laplace problem, this test case presents some interesting new features: The equations are discretized with schemes of different order of accuracy: second order for dif- fusion and third order for convection. Theoret- ically, the order of accuracy of the method is 2. However, for very large Reynolds numbers con- vection dominates, and so it is not guaranteed that the asymptotic order of accuracy will be in- dependent of the Reynolds number. · The analytical solution is independent of the Reynolds number. The convective term is zero and the source term on the right-hand side of
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<x<L H<z<O' B = 0.1L H = 0.0625L. For this slender and sharp-keeled body four sets of 16 H-O type grids were generated for a computation
-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 yz 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 yz 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.
0.074` 0.07 To. A O 0.066 11 Fix 0.062 c ~ acoooo 0 00a 0 0 _ 0 o o o o a o o ho o o ° ° e a Set A o Set B 4 0.058 I ~ 1 1 1 , 1 1 ~ 1 , 1 1 ~ 1 4 hi/h1 Figure 6: Cp at the "trailing edge" on the keel line as a function of the typical cell size. Turbulent flow around the Wigley hull. Rn = 7.4 x 106. Set A Set B CF CP ct CF | CP | Ct Pmin 0.15 0.11 0.15 0.16 Pmax 0.79 0.49 1.38 2.41 _ Table 3: Minimum and maximum orders of accuracy estimated from different grid triplets for the resistance coefficients of the turbulent flow around the Wigley hull. Cubic integration. Rn = 7.4 x 106. Nevertheless, in these difficult situations, the least squares root approach seems to be a robust alternative to obtain the error and uncertainty estimations. Figure 8 presents the estimation of the cell size zero solution with its band of uncertainty for the friction resistance obtained from the two sets of grids. The Cubic inte- gration data are adopted and the figure includes the results of the p and tl ~ methods. If more terms are included in the fixed exponents approach the uncer- tainty increases significantly due to the change in ARE, i.e. Up. The extrapolated values from the two sets of grids are consistent and as expected the set A leads to a larger uncertainty than set B. The tl ~ method may also be applied to Cp or Car leading to similar results. However, the p method does not work for Cp and Car, because as for the grid triplets none of the fits has O < p < 8. The location with x = 0.98L,y = 0.002L,z = 0.05L was selected to perform Verification studies for flow quantities that require interpolation from the grid nodes. The flow variables selected are the Carte- sian velocity components, Us, u2 and U3, the pres- 0.234 ., of x 0.233 0.232 0.231 0.7 F 0.6 0 5 0.4~o° 0 Linear, SetA >~,5° o Cubic, Set A 0.3j~ D Linear, SetB 0 Cubic, Set B 02 , , , 1 2 3 4 hill o <a'0 D Ooo ,, .,, ,O, an, ~ , c, . . 1 2 3 4 hi/h1 Linear, Set A o Cubic, Set A Linear, Set B Cubic, Set B a a a o g? hoot ,, 0 0 0 0.238 0.237 or x_ 0.236 0.235 ~~ 0 a ,- O 00OOOOO -Ro°° 5~ 5 a o o Linear, Set A 0 Cubic, Set A D Linear, Set B 0 Cubic, Set B 0.234 2 3 4 hi/h1 Figure 7: Resistance coefficients as a function of the typical cell size. Turbulent flow around the Wigley hull. Rn = 7.4 x 106. sure coefficient, Cp, the axial vorticity, sol, and the eddy-viscosity, vt. The results obtained with the four different types
0.24 0.235 of x 0.23 ~ ; 0.225 - o . ~ ~ . oii ~ 8 8 o ° · o to-(uo+up+ ~o+(uo+up+ o 4)~~ (U~~+Un+Uc) ~ ~ ~ V 0.22 0 tO+(uo+up+us) I , , , I 0.24 ni/n p expansion · > Set A o ~ Set B 0.235 To x 0.23 0.225 0.68 . 0.675 =0.67 0.665 9~ 0 Firsts ~ Firsts O Thirds ~ Thirds _ 0.66 2 3 4 hill 0.055 · 0.0545 Ite ~~i ~ ~ ~ is ~ ico O O O · ° 0 0 0 0 0 0 ° ~ .~. · ~ ~ · ~ 0.054 ~0 to-(uo+up+us) to+(uo+up+us) To ~o-(Uo+Up+Us) 0 22 - - (O+(uo+up+us) 1 , 1 2 3 hill t]1 expansion . 4 Figure 8: Error extrapolation and uncertainty bands as a function of the typical cell size of the coarsest grid for the friction resistance, CF. Least squares root ap- proach. Turbulent flow around the Wigley hull. Cubic integration. Rn = 7.4 x 106. Of interpolation tested showed that splitting the in- terpolation between x and the yz plane leads to results graphically identical to the volume interpola- tions. However, for some flow variables the order of the interpolation has a significant effect on the data. Furthermore, the interpolation introduces scatter in the data, which may be dependent on the variable selected and on the order of the interpolation tech- nique applied. In general, the first order interpolations produce more scatter in the data than the third order schemes. As an example of the problems introduced by the data interpolation, figure 9 presents U ~ and Cp as a function of the grid cell size for set A. , 0.0535 ~ 4~:, ~ 8900 8 ° ~ 0 Firsts ~ Firsts O Thirds O Thirds O ~ 0 Q.UJO . 2 3 4 hi/h1 Figure 9: Ul and Cp at x = 0.98L,y = 0.002L,z = 0.05L as a function of the typical cell size. Set A of grids. Turbulent flow around the Wigley hull. Rn = 7.4x 106. 3.5 Flow around the KVLCC2 Tanker The last test case is the turbulent flow around the KVLCC2 Tanker at model and full scale Reynolds numbers. It is a practical test case, which has been selected for the Workshop on CFD in Ship Hydro- dynamics Gothenburg 2000, [191. The calculations are also performed with the computer code PARNAS- SOS, [17], using the one-equation eddy-viscosity tur- bulence model proposed by Menter, [18], to model the Reynolds stresses. The flow is calculated without taking into account free-surface effects, double-body flow assumptions, in a computational domain bounded by the ship sur- face, the symmetry plane of the ship, the free-surface, which represents a symmetry plane, a vertical plane, x= constant, located 0.25L upstream of the ship,
x= - 0.25L, a vertical plane located 0.25L down- stream of the ship, x = 1.25L, and an elliptical cilinder that defines the outer boundary, given by: ~ Y ~ 2 ~ Z ~ 2 L2 At the outer boundary and at the inlet plane, x= 0.25L, the boundary conditions are derived from a potential flow calculation. 16 H-O type grids were generated for each Reynolds number. For the model scale calculations, the finest grid has 385 x 97 x 49 nodes and the coars- est grid has 145 x 37 x 19 nodes. Essentially the same grids are used at full scale Reynolds numbers, but the number of grid nodes in the direction normal to the hull surface is increased, varying now from 55 (coars- est grid) to 145 (finest grid) nodes. All the volume grids were created with a pro- prietary elliptic PDE grid generator, based on the GRAPE approach [161, using the same control pa- rameters in all the grids. However, unfortunately, this does not guarantee that the grids are strictly geomet- rically similar. In particular, we should remark that we have adopted an H-O topology. We do not intend to discuss the merits and drawbacks of the chosen topology. However, we must recognize that with this choice it is not possible to match the bow and stern contours to grid lines. An illustration of the grid at the boundaries is given in figure 10 for the 225 x 57 x 29 grid, where only every second grid line is plotted. For both Reynolds numbers, the flow quantities that coincide with grid nodes exhibit scatter in its variation with the typical cell size. This result confirms that one of the sources of scatter is the imperfect geometrical similarity of the grids, which is almost unavoidable in turbulent flows around complex geometries. Rn=5.6x 106 Rn=2.03x 109 CF CP C! CF CP Pmin 0.19 0.10 0.13 0.57 0.09 Pmax 6.44 6.02 7.31 6.83 5.22 no 10 79 74 9 68 ct 0.08 6.71 55 Table 4: Minimum and maximum orders of accuracy estimated from different grid triplets for the resistance coefficients of the turbulent flow around the KVLCC2 tanker. Cubic integration. The resistance coefficients, CF and Cp, for both Reynolds numbers are plotted as a function of the typical cell size in figure 11. The results obtained at both Reynolds numbers are similar. Cp shows a monotonous behaviour without any significant influ- ence on the type of integration performed. In contrast, CF exhibits an irregular variation which is dependent on the type of integration performed. In this case, there is scatter in the data obtained for both integral quantities. Table 4 presents the maximum and mini- mum values of p for CF, CP and C' at both Reynolds numbers, obtained from different grid triplets using the data of the Cubic integration. The total number of triplets tested is 92 and no stands for the number of grid triplets that satisfy the convergence criteria. If the scenario was not ideal for the Wigley hull, for the KVLCC2 it is even worse. n' is clearly smaller for CF than for Cp and Car, but the range of possible values of p extracted from the data is too large. The application of the least squares root approach in the p method is illustrated for the CF values at full scale obtained with the Cubic integration in figure 12. The results show that a significant part of the esti- mated uncertainty is due to the mean standard devia- Figure 10: Illustration of the boundary grid for the lion of the fit, Us. calculation of the flow around the KVLCC2 tanker. We have also examined the flow quantities, U i, U2, Us, Cp, At and v, at five different locations in the Apart from the complexity of the flow itself, the main difference of this case compared with the previ- ous one is the loss of the geometrical similarity of all the grids. As for the Wigley test case, we have examined flow quantities at locations which always coincide with grid nodes, integral parameters and variables at se- lected locations where the solution has to be interpo- lated from the grid nodes. . propeller plane: x = 0.9825L,y = 0.0075L,z = - 0.037L, x = 0.9825L,y = 0.0075L,z = - 0.042L, x = 0.9825L,y = 0.0075L,z = - 0.047L, x = 0.9825L,y = 0.0025L,z = - 0.042L
0.225 0.224 2.03 X 109 O 0 0 Linen 0.23 Cubic 0 0 0 0 0 0 0 c 0.223 0 x 0.222 0.221 0.22: 1'5 2 hi/h 0.478 0.477 0.476 To x 0.475 0.474 0.473 o 0 472 , ,11111 , 1 1.5 2 2.5 3 hi/h 1.4 1.2: To x cam 0.8 0.6 0.228 JO 0.226 0 0 x ~ 0.224 a O o o ° 0.222 1 ~ I 2.5 3 5.6x 106 o o o o o . o o a° ° o o o o to a a o o o Linen Cubic o 0 0 0 0 o~O OO 0 ~ D E> 0 Line~,5.6xlO6 Cubic, 5.6x106 Linear,2.03xlO9 O Cubic,2.03xlO9 0.4 . 1.5 - v 0 2 o o o o tO-UO fO+UO -(uO+u; +(uO+u; (uo+up+us) +(uo+up+ o o o ° ,, l\ a 0 22 , .. . . 1 1.5 2 2.5 hi/n1 o 3 Figure 12: Error extrapolation and uncertainty bands as a function of the typical cell size of the coarsest grid for the friction resistance, CF. Least squares root approach for the p method. Cubic integration. KVLCC2 tanker at Rn = 2.03 x 109. As for the Wigley hull test case, the results exhibit scatter, with an intensity that depends on the loca- tion, flow quantity and type of interpolation scheme applied. The comparison between the four different schemes tested is equivalent to the one obtained in the flow around the Wigley hull. Table 5 presents the estimated orders of accuracy at model scale Reynolds number obtained from different grid triplets at these five locations using the data of the Thirds interpolation scheme. 0.4 0.39 0.38 0.37 2 2.5 3 0.36 hi/h1 0.35 Figure 11: Resistance coefficients as a function of the typical cell size. Turbulent flow around the KVLCC2 tanker. and 0 Firsts ~ Firsts o Thirds O Thirdg n o O~ .................... 1 1.5 2 2.5 3 hi/h x = 0.9825L,y = 0.0125L,z = - 0.042L. Figure 13: Us at x = 0.9825L,y = 0.0075L,z = 0.037L as a function of the typical cell size. Turbu- lent flow around the KVLCC2 tanker. Rn = 5.6 x 106.
Pmin Oman nt Pmin Pmax nt Pmin Oman nt Pmin PmaJc nt Pmin Pmax nt 1 61)1 1 At x = 0.; ~ ~ 00 en, = ~ O]L 0.59 0.12 0.38 2.77 0.26 0.19 5.62 6.95 7.09 7.57 5.08 7.82 80 75 20 30 75 49 x = 0.9825L,y = 0.0075L,z = - 0.042L 0.15 0.15 0.87 3.63 0.11 0.76 3.49 7.90 7.82 7.85 5.61 6.75 80 57 12 25 82 15 x = 0.9825L,y = 0.0075L,z = - 0.047L 0.14 0.13 1.18 0.63 0.08 1.14 6.69 0.84 7.22 4.66 7.79 1.14 70 2 43 84 54 1 x = 0.9825L,y = 0.0025L,z = - 0.042L 0.13 0.32 0.13 1.27 0.42 0.10 4.89 7.61 5.68 7.41 6.82 6.66 71 12 80 29 48 54 x = 0.9825L,y = 0.0125L,z = - 0.042L 0.82 0.43 0.37 2.84 0.13 0.35 7.72 7.54 7.46 7.61 6.48 7.81 48 89 70 58 78 35 Table 5: Minimum and maximum orders of accuracy estimated from different grid triplets for Us, U2, U3, Cp, mi and v', at selected locations of the turbulent flow around the KVLCC2 tanker. Thirds interpola- tion. Rn = 5.6 x 106. These results clearly indicate that it is extremely difficult to base the error and uncertainty estimation on grid triplets when there is scatter in the data. As an illustration of the performance of the least squares root approach, the fits to the data of the Thirds inter- polation, obtained with the p, tl and t2 methods for U1 at x = 0.9825L,y = 0.0075L,z = -0.037L in the model scale flow are plotted in figure 14. The data are plotted in figure 13 as a function of the typical cell size. In these cases, which include scatter, the solution with the minimum number of grids required by the power series expansion adopted was disregarded. The estimate of Up for the fixed exponents expansion with more than one term leads to much larger values than the p method. On the other hand, the tl method which uses a single term performs better because Up is dominant in the uncertainty estimation. 4 Conclusions This paper has presented an evaluation of Verifica- tion Procedures for Computational Fluid Dynamics. 0.5 0.45 :~ 0.4 0.35 0.3. 0.5 0.45 0.35 - tl expansion 0 3 1 , , , I , · 1 1.5 high 2-5 3 · ho V (o~Uo V (o+Uo ~ (o-(uo+u) · (o+~o+Up) o (o-(uo+up+us) o (o+~o+up+us) ~~ ~ ~ ~ , e p expansion . . . . . . . . . . . . . . . . 1.5 hill 2.5 3 o Polo +UO o-(uo+u) (o+~o+Up) (o-(uo+up+Us) ¢~+(u^+up+us) e' ~ 8 8 o 5 t2 expansion 0.45 on o 0.4~ .~' 0.35 0.3 e ~ e e I , , , , I , , 1 1 1 ~ ~ I I I 1.5 ho, 2.5 3 Figure 14: Error extrapolation and uncertainty bands as a function of the typical cell size of the coars- est grid for U1 at x = 0.9825L,y = 0.0075L,z = 0.037L. Turbulent flow around the KVLCC2 tanker. Rn=5.6x 106. Six test cases were selected to investigate the diffi- culties of the error and uncertainty estimation in the
calculation of the flow around a ship hull. Two types of power series expansions were applied to estimate the error, using unknown and fixed ex- ponents, and a procedure to estimate the uncertainty when the number of grids available is larger than the minimum required is tested. The uncertainty estima- tion is composed of three quantities related to the vari- ability of the error estimation when different sets of grids are selected, the mean standard deviation of a least squares root fit and to the contribution of the high-order terms. The selected test cases cover a wide range of flows that enabled us to study the effects of the follow- ing aspects: geometrical similarity of the grids, non- orthogonality and stretching of the grids, different or- ders of accuracy in the discretization of the flow equa- tions, introduction of turbulence models and integra- tion and interpolation of the computed data. From the results obtained in this study the follow- ing conclusions have been drawn: · In our practice, it is not yet feasible to attain the asymptotic range with grids of affordable density and so it is difficult to establish the asymptotic order of accuracy. · In the analytical 2-D test cases, the use of stretched or non-orthogonal grids did not show any deterioration of the observed order of accu- racy. · The use of a turbulence model did not show any significant influence on the grid convergence properties of the calculation of the flow around the Wigley hull. · For data without scatter, the following trends were found: · The one term expansion with unknown ex- ponent and the two terms expansion with fixed exponents seems to be two most ef- fective procedures to estimate the error. In the fixed exponents approach, the observed order of accuracy is needed to select the proper exponents. · The proposed estimation of the uncertainty with three contributions performs well for the test cases with a known analytical solu- tion. - The present procedure gives a poor esti- mate for the uncertainty of the single term expansion with a fixed exponent. The esti- mation of the effect of the high-order terms is too crude in this case. · The results of all the test cases show the expected behaviour of the uncertainty with an increase of the uncertainty band with the coarsening of the grids. - The use of grid triplets, doublets or quadru- plets or the least squares root approach leads to similar results. Nevertheless, the latter approach seems to handle better the effect of the coarsest grids. · Two sources of scatter were observed in the present test cases: - The need to use numerical interpolation techniques when a verification procedure is applied to a variable at a location which is not a grid node common to all grids. - The lack of geometrical similarity of the grids. · In complex turbulent flows around ships these two origins of scatter are almost unavoidable. Therefore, in Verification procedures for turbu- lent flows around ship hulls, one has to choose conservative values for the uncertainty or to gen- erate solutions on more grids to be better able to handle the effects of the scatter. . In the calculations performed for the Wigley hull and the KVLCC2 tanker, the following trends were observed in the flow quantities that exhibit scatter: - The intensity of the scatter depends on the location selected, on the flow quantity se- lected and on the order of the interpolation scheme used. The Reynolds number has no significant ef- fect on the intensity of the scatter. The scatter produces an excessive variation in the observed order of accuracy obtained from different grid triplets. Understand- ably, the least squares root approach can cope better with the scatter. The use of power series expansions with fixed exponents and terms of order larger than two leads to a large variation in the er- ror estimation from different sets of grids. In these cases, the use of first (when nec- essary) and second order terms seems to be preferable because the uncertainty is dom- inated by the changes in the error estima- tion.
Acknowledgement We are indebted to Hoyte Raven (MARIN) for many fruitful discussions on the results of our veri- fication studies. References t1] Guide for the Verification and Validation of Computational Fluid Dynamics Simulations, AIAA-G077-1998. [2] Best Practice Guidelines, Version 1.0, ER- COFTAC Special Interest Group on "Quality and Trust in Industrial CFD", January 2000. t3] I11C Quality Manual t4] Easterling R.G. - Measuring the Predictive Ca- pability of Computational Methods: Prin- ciples and Methods, Issues and Illustrations - Sandia Report, SAND2001-0243, February 2001. [5] Roache P.J. - Verification and Validation in Computational Science and Engineering - Hermosa Publishers, 1998. [6] Jasak H. - Error Analysis and Estimation in the Finite Volume Method with Applications to Fluid Flow - PhD Thesis, Imperial College, University of London, 1996. [7] E;a L., Hoekstra M. - On the Application of the Moment Error Estimate to Error Quantifi- cation in CFD - 3st Marnet-CFD Workshop, Creete 2001. [8] Oberkampf W. - Private communication [91 Hoekstra M., E;a L. - An Example of Error Quantification of Ship-Related CFD Results - 7th Numerical Ship Hydrodynamics Confer- ence, Nantes, July 1999. [10] E;a L., Hoekstra M. - On the Numerical Ver- ification of Ship Stern Flow Calculations - At Marnet-CFD Workshop, Barcelona 1999. t11] E;a L., Hoekstra M. - On the Application of Verification Procedures in Computational Fluid Dynamics - 2st Marnet-CFD Workshop, Copenhagen 2000. [12] E;a L., Hoekstra M. - An Evaluation of Ver- ification Procedures for Computational Fluid Dynamics - IST Report D72-7, June 2000. [13] Stern F., Wilson R., Coleman H.W., Paterson E.G. - Comprehemsive Approach to Verifica- tion and Validation of CFD Simulations - Part 1: Methodologies and Procedures - ASME Journal of Fluids Engineering, Vol. 123, pp. 803-810. December 2001. [14] E;a L., Hoekstra M. - Verification Procedures for Computational Fluid Dynamics - IST Re- port D72-14, April 2002. [15] Kobayashi M.H., Pereira J.M.C, Pereira J.C.F. - A Conservative Finite-Volume Second-Order Accurate Projection Method on Hybrid Un- structured Grids. - Journal of Computational Physics 150, pp. 40-75, 1999. [161 E;a L., Hoekstra M., Windt J. - Practical Grid Generation Tools with Applications to Ship Hydrodynamics - 8th International Con- ference in Grid Generation in Computational Field Simulations, Hawaii, June 2002. t17] Hoekstra M., E;a L. - PARNASSOS: An Effi- cient Method for Ship Stern Flow Calculation - Third Osaka Colloquium on Advanced CFD Applications to Ship Flow and Hull Form De- sign, Osaka, Japan, 1998. [18] Menter F.R. - Eddy Uscosity Transport Equa- tions and Their Relation to the k£ Model - Journal of Fluids Engineering, Vol. 119, De- cember 1997, pp. 876-884. t19] Larsson L., Stern F., Bertram V. Eds. - Gothenburg 2000 A Workshop on Numerical Ship Hydrodynamics - September 2000.
DISCUSSION C. Yang George Mason University, USA I believe that authors have done excellent research in the shape optimization. We have done the similar work using CFD tools together with gradient-based optimization techniques to minimize wave drag. We have found that the hull form optimized for a single design speed may yield less desirable hull form for other speeds. Have authors tried multi-speed optimization or checked hydrodynamic characteristics of the optimal hull form obtained at other speeds? I would also like to know what is the geometry constraint and how many design variables are used in the current optimization process? AUTHORS' REPLY For this optimization exercise, only 6 control points moved in the cross direction are considered as design variables. Since the shape modification is quite local, the only constraints included concern the range of displacement of the control points. No constraint involving the volume for instance is taken into account. The comparison of the results obtained for different Reynolds numbers shows for these exercises similar tendencies for the shapes, although some differences are noticed which may have a significant influence on the final performance. More global problems involving wave drag minimization were not tested. Multi- objective optimizations using genetic algorithms may provide an interesting answer for these multi-points problems.
