# Flow simulation in an elbow diffuser: verification and validation

Sebastiano Mauri, Jean-Louis Kueny, François Avellan

► To cite this version:

Sebastiano Mauri, Jean-Louis Kueny, François Avellan. Flow simulation in an elbow diffuser: verification and validation. 21st IAHR Symposium on hydraulic machinery and systems, Sep 2002, Lausanne, Switzerland. pp.213-240. hal-00262263

HAL Id: hal-00262263

https://hal.science/hal-00262263v1

Submitted on 19 Mar 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# FLOW SIMULATION IN AN ELBOW DIFFUSER:
VERIFICATION AND VALIDATION

Sebastiano MAURI,
Jean-Louis KUENY,
François AVELLAN,

LMH/Swiss Federal Institute of
Technology,
Lausanne, Switzerland

## ABSTRACT

Numerical simulation of the unsteady turbulent flow in a three-dimensional draft tube geometry is performed. The investigation is carried out with a commercial finite volume solver implementing the Reynolds averaged Navier-Stokes equations. The modeling of most practically relevant turbulent flows continues to be based on this equation system. For this reason it is important to evaluate the limitations of this approach. Verification and validation are presented; detailed measurements are compared with computations over a wide range of operating conditions.

## RÉSUMÉ

Une simulation numérique de l'écoulement turbulent instationnaire dans un diffuseur coudé tridimensionnel est effectuée. L'étude est menée avec un code commercial résolvant les équations de Navier Stokes moyennées en formulation volume fini. La modélisation de la plupart des écoulements turbulents rencontrés en pratique continue à être basée sur ce système d'équations. Pour cette raison, il est important d'évaluer les limitations de cette approche, en particulier pour la prédiction des écoulements instationnaires. Des vérifications et validations sont présentées; des mesures détaillées sont comparées aux calculs.

## NOMENCLATURE

<table><thead><tr><th>Term</th><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td>Reference section area</td><td>A<sub>ref</sub></td><td></td></tr><tr><td>Velocity field</td><td>C</td><td></td></tr><tr><td>Axial velocity component</td><td>C<sub>a</sub></td><td></td></tr><tr><td>Normal velocity component</td><td>C<sub>n</sub></td><td></td></tr><tr><td>Mean local normal velocity</td><td>C<sub>o</sub></td><td></td></tr><tr><td>Radial velocity component</td><td>C<sub>r</sub></td><td></td></tr><tr><td>Runner outlet diameter</td><td>D</td><td>0.4 [m]</td></tr><tr><td>Turbulent dissipation rate</td><td>ε</td><td></td></tr><tr><td>Turbulent eddy length scale</td><td>L<sub>ε</sub></td><td></td></tr><tr><td>Turbulent kinetic energy</td><td>k</td><td></td></tr><tr><td>Flow rate</td><td>Q</td><td></td></tr></tbody></table>

<table><thead><tr><th>Term</th><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td>Radius</td><td>r</td><td></td></tr><tr><td>Inlet radius</td><td>R</td><td></td></tr><tr><td>Reynolds number</td><td>Re</td><td>C<sub>o</sub>D/v</td></tr><tr><td>Non-dimensional distance from the wall</td><td>y<sup>+</sup></td><td></td></tr><tr><td>Pressure recovery factor</td><td>χ</td><td>ΔP<sub>dt</sub> / ∫S(Q/A<sub>ref</sub>)<sup>2</sup> dv</td></tr><tr><td>Mean wall pressure difference between draft tube inlet and outlet</td><td>ΔP<sub>dt</sub></td><td></td></tr><tr><td>Water kinematic viscosity</td><td>ν</td><td></td></tr><tr><td>Flow rate coefficient</td><td>φ</td><td></td></tr><tr><td>Flow rate coefficient divided by the φ of the best efficiency point</td><td>φ*</td><td></td></tr><tr><td>Water density</td><td>ρ</td><td></td></tr><tr><td>Half cone opening angle</td><td>θ</td><td></td></tr></tbody></table>
INTRODUCTION

Computational results and comparisons with experimental data are shown in order to estimate modeling and discretization errors. A full discussion can be found in Mauri (Ref. 6) and the analysis of the flow is also presented in Mauri et al. (Ref. 5). While the results depend on the interaction between the different parameters, they are here discussed as independent variables for simplicity. Investigated parameters are: mesh density, inlet radial velocity component, near wall velocity profile, and turbulence modeling. Investigations are carried out by varying a parameter one at a time. The other parameters correspond to the best choices resulting from each investigation. The influence of the inlet turbulence and the outlet boundary conditions were already discussed in Mauri et al. (Ref. 1) and more recently in Mauri (Ref. 6). The computations discussed here are performed on the geometry with the simple addition of a downstream channel (Fig. 1) and the inlet eddy length scale is set to $L_ε=0.001$ (0.002D) for all operating points.

<figure><img src="image_2.png" /><figcaption>Fig. 1 a) Investigated geometry. Cross area evolution. b) Outlet boundary conditions c) Velocity and kinetic energy inlet profiles (filled points correspond to the external radius). Machine and draft tube efficiency with the main operating points. $r/R=0$: center, $|r/R|=1$: wall.</figcaption></figure>
## CASE

Experimental studies are carried out on a vertical axis reduced scale model (1:10) of an existing water turbine at the test rig facilities of the EPFL, Laboratory for Hydraulic Machines. A high specific speed ($\nu$=0.56) Francis runner supplies the symmetrical draft tube with a single pier (Fig. 1). The numerical flow analysis is carried out at a constant head for 14 flow rates ranging from 90% to 110% of the best efficiency discharge. The Reynolds number based on the inlet mean velocity and diameter is $Re \approx 1.510^6$.

The inlet conditions are experimentally investigated on the symmetry axis diameter at six operating points by means of the LDA technique. The three components of the velocity and the Reynolds' stress tensor (Fig. 1) are obtained through four positions of a 2D laser probe, details can be found in Ciocan et al. (Ref. 2). The measurements uncertainties are estimated to be less than 3%. The other conditions are linearly interpolated from the measurements.

## MODELING AND VERIFICATION

3D steady Reynolds-averaged Navier-Stokes (RANS) flow computations with several turbulence models and logarithmic wall functions are performed. The commercial code CFX-TASCflow 2.10 is used. Details on the code are given in Raw et al. (Ref. 8) and in Ref. 1. It is assumed that in the operating range considered in this study the main flow characteristics can be described assuming flow steadiness. However, measurements show regions with strong fluctuations in particular for the operating at higher flow rate. Time dependent computations are discussed in Mauri (Ref. 6). The computation is considered converged to the steady solution, when the value of the maximal normalized equation residual is less than $10^{-4}$. No appreciable difference is observed. All operating points have been computed by increasing or decreasing the flow rate using the solution of the previous point as initial state, showing the uniqueness of the solution.

## MESH AND GRID CONVERGENCE

The geometry of the draft tube model is discretized with a structured multiblock mesh. A butterfly topology with a C-shaped grid around the pier is used. The minimal skew angle in the cells is 42° and the maximum aspect ratio of the cells is 14. The $y^+$ values of the first grid points off the wall remain within the range of 20 and 300 for all operating points, where the majority of the points lies between $30 < y^+ < 100$. An a posteriori numerical error estimation based on the generalized Richardson extrapolation, not requiring any restriction to integer refinement and applicable to solution functionals, is carried out. Following Roache (Ref. 8) the more conservative Grid Convergence Index (GCI) is also reported. The recovery factor obtained with four meshes are compared at $\phi^*=0.994$ in Table 1. The number of nodes corresponds to the effective number of points in the draft tube geometry (overlapping nodes at the block interfaces are counted only once and the downstream channel is not considered). The grid refinement is reported simply in terms of the total number of grid points $N$ used in the two meshes as $n_j=(N_i/N_j)^{1/3}$. The meshes have the same topology but the grid refinement is not uniform in the space. Error estimation using unrelated grids poses a challenge. It is
however mostly the case when working with complex three-dimensional meshes. Depending on the mesh resolution at the inlet the resulting flow rate will change and must be corrected to retrieve the measured value by multiplying the velocity field by a factor. Slight differences in the inlet boundary conditions introduce an additional uncertainty. The first three meshes seem to be in the asymptotic range: actual asymptotic rate of convergence $p=2.2$ to be compared with the theoretical order $p=2$.

<table><caption>Table 1 $\varphi*=0.994$. $N$: number of nodes, $r$: grid refinement ratio, $\chi$: recovery factor, $E$: estimated fractional error, $GCI$: grid convergence index, $p$: scheme order. Values for $E$, $GCI$ are reported in [%] using $p=1-2$.</caption><tr><td>$N_1=1'855'152$</td><td>$N_2=1'107'237$</td><td>$N_3=633'720$</td><td>$N_4=328'360$</td></tr><tr><td>$E_1=0.17-0.08$</td><td>$E_2=0.20-0.11$</td><td>$E_3=0.28-0.16$</td><td></td></tr><tr><td>$GCI_1=0.51-0.23$</td><td>$GCI_2=0.61-0.33$</td><td>$GCI_3=0.83-0.49$</td><td></td></tr><tr><td>$\chi_1=0.7739$</td><td>$\chi_2=0.7737$</td><td>$\chi_3=0.7733$</td><td>$\chi_4=0.7826$</td></tr><tr><td>$n_2=1.187717$</td><td>$n_2=1.204429$</td><td>$n_3=1.430520$</td><td>$n_4=1.781051$</td></tr></table>

This results are expected to be only partially representative for other operating conditions due to the important flow differences. The results for the meshes $N_2$ and $N_4$ are compared also at the points $\varphi*=0.919-1.108$ in Table 2.

<table><caption>Table 2 a) $\varphi*=0.919$, b) $\varphi*=1.108$. See caption in Tab. 1.</caption><tr><td>$N_2=1'107'237$</td><td>$N_4=328'360$</td><td>$N_2=1'107'237$</td><td>$N_4=328'360$</td></tr><tr><td>$r_{24}=1.499559$</td><td></td><td>$r_{24}=1.499559$</td><td></td></tr><tr><td>$\chi_2=0.5783$</td><td>$\chi_4=0.5584$</td><td>$\chi_2=0.5305$</td><td>$\chi_4=0.5385$</td></tr><tr><td>$E_2=6.89-2.75$</td><td>$E_4=10.33-6.20$</td><td>$E_2=3.02-1.21$</td><td>$E_4=4.53-2.72$</td></tr></table>

As observed for $\varphi*=0.994$ we can expect that the solution obtained with $N_4$ is not in the asymptotic range and consequently the estimated fractional error is only indicative. As expected the mesh influence for the operating points lying outside the optimal range increases. The skin friction lines show however only slight differences, visible in particular in the backflow zone for $\varphi*=0.919$. For operating points lying in the optimal range the mesh with $N_3=633'720$ effective nodes insures a $GCI$ smaller than 1% for the recovery factor. While the mesh with $N_4=328'360$ does not lie in the asymptotic range in the prediction of the recovery factor, the results obtained with this mesh show the same flow topology as the finer meshes and compare even slightly better with the measured velocity and pressure profiles.

For these reasons and the limited computational resources the investigations are carried out with the mesh consisting of $N_4=328'360$ nodes.
## INLET RADIAL VELOCITY COMPONENT

Some discrepancies on the radial component of the velocity are observed between the LDA and steady probe measurements. The probe indicates a velocity vector that is parallel to the surface near the walls, with a nearly linear evolution in the interior. LDA results show a slope indicating a more pronounced effect due to the bend. For this reason two different distributions are investigated. While the axial and tangential velocity components are interpolated by cubic splines, this type of interpolation leads to oscillation for the radial component. This is due to the inferior number of measurement points and the poor smoothness of the profile. The profile is therefore interpolated by the best linear fit. The profiles for $\phi*=0.919$, $\phi*=0.994$ and $\phi*=1.108$ are illustrated in Fig. 2. The dotted line describes the evolution of the radial component with the equation $C_r=C_n\tan(\theta r/R)$ (hereafter "geometrical" distribution), $\theta$ being the cone half opening angle. The inclination of the velocity vector in the radial direction is therefore determined by the geometry of the cone (similarly to the steady probe measurements).

<figure><img src="image_3.png" /><figcaption>Fig. 2 a) Inlet radial velocity distribution: 1) $\phi*=0.919$, 2) $\phi*=0.994$, 3) $\phi*=1.108$. Points: LDA measurements, dotted line: best linear fit to the measurements, solid line: "geometrical" distribution. b) Skin friction lines, view from above, $\phi*=1.108$: 1) "geometrical" distribution, 2) best linear fit distribution.</figcaption></figure>

The influence of the inlet radial velocity component is investigated for $\phi*=0.919$, $\phi*=0.994$, $\phi*=1.108$. The recovery factor is affected by this choice as illustrated in Table 3.

<table><caption>Table 3 Influence of the inlet radial velocity component on the recovery factor.</caption><thead><tr><th></th><th>φ*=0.919</th><th>φ*=0.994</th><th>φ*=1.108</th></tr></thead><tbody><tr><td>best linear fit</td><td>χ=0.5642</td><td>χ=0.7955</td><td>χ=0.6224</td></tr><tr><td>"geometrical" distr.</td><td>χ=0.5584</td><td>χ=0.7826</td><td>χ=0.5385</td></tr><tr><td>Measured</td><td>χ=0.4937</td><td>χ=0.7584</td><td>χ=0.5192</td></tr></tbody></table>

The radial velocity component can affect considerably the flow in the draft tube in spite of the small magnitude of this component. While for the operating points at $\phi*=0.919$ and $\phi*=0.994$ the differences using the two distributions are relatively small, at $\phi*=1.108$ the flow differs considerably. The "geometrical" distribution leads to better results in comparison with the measurements for all operating points.

## NEAR WALL VELOCITY PROFILE

The measured velocity profiles are interpolated by cubic splines. The nearest measurement point is at 0.076 inlet diameters from the wall. In order to control the extrapolation at the wall, the velocity is here imposed to be a factor $f$ of the nearby interior measurement point value. This factor plays a role on the wall velocity gradient and consequently on the wall
friction. The swirl is proportional to the squared radius and therefore the velocity near the wall can have an important weight. The influence of this parameter is investigated for the extreme points at $\phi*=0.919$ and $\phi*=1.108$. The recovery factor is clearly affected by this choice, varying almost linearly with f as illustrated in the table of Fig. 3. The importance of the momentum thickness on the pressure recovery is well known. The stronger influence occurs for the point with the highest swirl, where the boundary flow energy is higher. The influence on the skin friction lines at $\phi*=1.108$ is illustrated in Fig. 3. The flow angles in the cone are strongly affected especially at the highest flow rate. The near-wall steep gradients observed experimentally are reproduced in the computations only in the case f=0.900.

a

<table><thead><tr><th></th><th>φ*=0.919</th><th>φ*=1.108</th></tr></thead><tbody><tr><td>f=0.005</td><td>χ=0.6099</td><td>χ=0.5860</td></tr><tr><td>f=0.215</td><td>χ=0.6121</td><td>χ=0.5772</td></tr><tr><td>f=0.900</td><td>χ=0.5584</td><td>χ=0.5385</td></tr><tr><td>meas.</td><td>χ=0.4937</td><td>χ=0.5192</td></tr></tbody></table>

b

<figure><img src="image_4.png" /><figcaption>Fig. 3 a) Influence of the near wall velocity profile on the recovery factor. b) Skin friction lines, view from behind, $\phi*=1.108$. 1) f=0.005, 2) f=0.215 and 3) f=0.900.</figcaption></figure>

## TURBULENCE MODELING

The use of the k-ω and BSL<fnref n="1" /> models does not improve the agreement with the experimental data in comparison with the standard k-ε, even for the operating point at $\phi*=1.108$ where a large separated region is found to play an important role on the recovery factor, as shown in Table 4. Reynolds stress models are characterized by a higher degree of universality and should better describe in particular the effects of streamline curvature and secondary flows as well as the characteristics of swirling flows. The recovery factors obtained using the RSM-LLR model are reproduced in Table 4. It should be mentioned however that the code does not allow specifying the measured stress profiles at the inlet. These are computed from the turbulent kinetic energy. Even if the recovery factor matches better with the use of the RSM model for the operating point at $\phi*=0.994$, the comparisons with the measured velocity and pressure profiles indicate better performances for the standard k-ε for both operating points.

<table><caption>Table 4 Influence of the turbulence model on the recovery factor.</caption><thead><tr><th></th><th>φ*=0.919</th><th>φ*=0.994</th><th>φ*=1.108</th></tr></thead><tbody><tr><td>k-ω</td><td>-</td><td>χ=0.8012</td><td>χ=0.6318</td></tr><tr><td>BSL</td><td>-</td><td>χ=0.8021</td><td>-</td></tr><tr><td>k-ε</td><td>χ=0.5584</td><td>χ=0.7826</td><td>χ=0.5385</td></tr><tr><td>RSM-LLR</td><td>χ=0.5675</td><td>χ=0.7652</td><td>-</td></tr><tr><td>meas.</td><td>χ=0.4937</td><td>χ=0.7584</td><td>χ=0.5192</td></tr></tbody></table>

## VALIDATION

In this section the flow field obtained with the modeling choices previously discussed and summarized in the table of Fig. 4., are compared with the experimental data.

<footnote n="1">¹blends between the k-ω model, applied near the surface and the k-ε model applied outside the boundary layer</footnote>
## RECOVERY FACTOR

<table><thead><tr><th>domain</th><th>box</th></tr></thead><tbody><tr><td>mesh density</td><td>328'360</td></tr><tr><td>radial velocity</td><td>geometrical</td></tr><tr><td>near wall vel. extrapolation</td><td>f=0.900</td></tr><tr><td>inlet turbulent dissipation</td><td>L<sub>ε</sub>=0.001</td></tr><tr><td>advection scheme</td><td>MLPS</td></tr><tr><td>turbulence model</td><td>k-ε</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig. 4 a) Computational parameters. b) Static pressure recovery. Comparison measurement computation. The GCI is reported for $\phi*=0.994$.</figcaption></figure>

The recovery factor is compared with the measurements in Fig. 4. The overall agreement is fairly good. The pressure recovery drop takes place at the same flow rate which is observed experimentally. However, the computations overestimate the recovery factor over the whole range, with a maximal difference of 14% of the measured value. Near the best efficiency conditions the maximal difference is 6%.

## VELOCITY AND PRESSURE PROFILES, SECTIONS 1.75 AND 15.5

The velocity and pressure profiles acquired with a five-sensor steady probe on 16 measurement axes on sections 1.75 and 15.5 are compared with the computations for an operating condition in Fig. 5. On the whole the flow is fairly well predicted. Locally important differences occur for all velocity components. Despite the short distance from the inlet, at the section 1.75 for the operating point at higher swirl corresponding to $\phi*=0.919$, the maximum velocity difference reaches $0.5C_o$. The averaged difference value is $0.1C_o$. This is explained by the known difficulties of the k-ε model to correctly simulate swirling flows. Just upstream the pier on section 15.5, the averaged difference values range from .1 to $0.3C_o$ and the maximum difference is $0.5C_o$. The reference pressure corresponds to the value at the outlet ring manifold. A similar agreement is observed also at the others operating points.

## WALL PRESSURE, SECTIONS 1.3, 1.75, 6.5, 9.5, 12.5

The wall static pressure is compared at six sections in Fig. 5 for three operating conditions. The global agreement is fairly good. Locally the differences reach 80% of the measured value. The computations clearly overestimate the bend influence at the section 1.75, but the differences decrease in the following sections.
<figure><img src="image_6.png" /><figcaption>Fig. 5 a) Wall static pressure, sections 1.3, 1.75, 6.5, 9.5, 12.5. $\phi^*=0.919$, $\phi^*=0.994$, $\phi^*=1.108$ (from the left to the right). The pressure values are divided by the mean value on section 1.3 and the reference is set at the same section. Solid line: computation, dotted lines: experiments. b) Sections 1.75 (left) and 15.5 (right), $\phi^*=0.919$: comparisons measurements (points) - computations (line). Section 1.75: from the wall to the center, section 15.5: from the upper wall to the bottom. p: static pressure, $c_n$: normal velocity component, $c_r$: radial velocity component, $c_{\theta}$, $c_{t1}$, $c_{t2}$: tangential velocity component, $E=p+1/2\rho c^2$: total pressure. The velocity is divided by the mean normal velocity $C_o$ of the investigated section and the pressure by $1/2\rho C_o^2$. The reference pressure corresponds to the value at the outlet ring manifold.</figcaption></figure>

## VELOCITY FIELD, SECTIONS 20.75

Two velocity components acquired with the LDA system at 950 points on sections 20.75 are compared with the computations in Fig. 6. Better agreement is observed for the extreme points. The central points are less satisfactory, probably due to rapid flow changes occurring in this region.
<figure><img src="image_7.png" /><figcaption>Fig. 6 a) Outlet velocity field, comparisons with measurement $\phi^*=0.919$, $\phi^*=0.973$, $\phi^*=0.994$
b) $\phi^*=1.027$, $\phi^*=1.054$, $\phi^*=1.108$.</figcaption></figure>

## SUMMARY

The flow in the draft tube is complex because of its unsteady three-dimensional and rotating nature. This work attempts to define the prediction capability of RANS computations. Comparisons with detailed experimental data over an extended range of operating conditions reveal the limitations of the approach and the sensitivity on the inlet boundary conditions. After the parameters are calibrated, the mean flow field is however correctly captured with main trends. Quantitatively the static pressure recovery is reasonably well predicted while locally the flow details are partially missed. The numerical optimization of this component without measurements could be therefore hazardous. In spite of the detailed measurements, it is clear that insufficient data is present for the definition of the mean flow boundary conditions. A great effort is underway in the development of advanced turbulence models such as LES methods. These are today applicable to complex geometries. The use of advanced models however, requires even more detailed information at the boundaries, making this approach even more difficult to properly apply.

## ACKNOWLEDGEMENT

This study is established within of the Flindt project (Eureka N° 1625) whose participants are: EPFL, Electricite de France, Alstom Power, General Electric Canada, VA Tech Escher Wyss, Voith Hydro who are gratefully acknowledged. The project is also supported by the Swiss Federal Commission for Technology and Innovation (PSEL) and the German Ministry of Science and Technology (BMBF). The authors wish to thank G. Ciocan and J. Arpe for the
experimental data and the staff of the Laboratory for Hydraulic Machines for providing
support in the experimental and numerical work.

REFERENCES

Ref. 1 AEA Technology Engineering Software Limited, 2000, "CFX-TASCflow Theory Documentation Version 2.1", Waterloo, Ontario, Canada N2L 5Z4.

Ref. 2 Ciocan G. D., Kueny J.-L., Avellan F., 2000, "Optical Measurement Techniques for Experimental Analysis of Hydraulic Turbines Rotor-Stator Interaction", *Proceedings of 2000 ASME Fluids Engineering Conference*, Boston, June 11-15.

Ref. 3 Idelcik I. E., 1969, "Memento des pertes de charge", Eyrolles, Paris.

Ref. 4 Mauri S., Kueny J.-L., Avellan F., 2000, "Numerical Prediction of the Flow in a Turbine Draft Tube Influence of the Boundary Conditions", *ASME Fluids Engineering Conference*, Boston, June 11-15.

Ref. 5 Mauri S., Kueny J.-L., Avellan F., 2002, "Werlé-Legendre Separation in a Hydraulic Machine Draft Tube", *Proceedings of 2002 ASME Fluids Engineering Conference*, Montreal, Quebec, Canada, July 14-18.

Ref. 6 Mauri S., 2002, "Numerical Simulation and Flow Analysis of an Elbow Diffuser", *EPFL thesis N°* 2527, Lausanne, http://lmhwww.epfl.ch/Publications/Thesis/Mauri/MauriPhD.html.

Ref. 7 Raw M. J., Galpin P. F., Hutchinson B. R., 1989, "A Colocated Finite Volume Method for Solving the Navier-Stokes Equations for Incompressible and Compressible Flows in Turbomachinery: Results and Applications", *Thirty-Sixth Annual General Meeting Canadian Aeronautics and Space Institute*, Engineering Conference, Boston, 13.1-13, Ottawa, Ontario.

Ref. 8 Roache P. J., 1997, "Quantification of Uncertainty in Computational Fluid Dynamics", *Ann. Rev. Fluid. Mech.* 29, pp. 123-160.