Effects of uncertainties in positioning of PIV plane on validation of CFD results of a high-head Francis turbine model

The use of Computational Fluid Dynamics (CFD) for the design of modern hydraulic turbines has increased and matured signi ﬁ cantly in the last decades. More recently, CFD is also used to understand how to safely widen the hydraulic turbine operating ranges, and avoid hazardous conditions during transient operation. The accuracy of such CFD results relies on validation with experimental data which contains numerous sources of uncertainties. The present work is focusing on the effects of the un- certainties in the positioning of the experimental Particle Image Velocimetry (PIV) plane on the validation of CFD results of the high-head Francis-99 turbine model. A transient shutdown sequence is considered, where the available experimental and numerical data are considered accurate according to a conventional thorough validation procedure. A part of that validation procedure is the comparison of spatially and temporally varying velocity pro ﬁ les along three lines of the experimental PIV plane. The positioning of this PIV plane is here considered uncertain, using three translational and three rotational stochastic parameters with uniform probability distribution functions. The validated CFD results are used to extract the data that depends on these uncertainties, while this is not possible for the experimental data. The polynomial chaos expansion method is employed for the Uncertainty Quanti ﬁ cation (UQ) study while Sobol ’ indices are utilized for the Sensitivity Analysis (SA). The UQ can be used to show how the considered uncertainties impact the extracted components of the velocity ﬁ eld, and the sensitivity analysis reveals the relative contribution of each uncertain parameter to the quantity of interest. For this particular case it is shown that the so-called horizontal velocity component is most sensitive to the plane-normal positioning of the PIV plane. This is also the velocity component where all the numerical results found in the literature differ most from the experimental results. It is also shown that the probability distribution function of the numerical horizontal velocity is covered by the experimental standard deviation bounds, which means that it is quanti ﬁ ed that the numerical and experimental results are similar within the range of the uncertainties. © 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Computational Fluid Dynamics (CFD) has been employed for decades to study the flow in hydraulic turbines [1]. Advancement of numerical methods as well as significant progress in available computational resources led engineers and researchers to employ CFD methods for design improvement and optimization of hydraulic turbines [2e11]. In recent years, the introduction of intermittent energy sources relies on the flexibility and fast responsiveness of the hydraulic turbines [12,13]. This requires extensive numerical studies of off-design and transient operation of the hydraulic turbines [14]. The performance and reliability of modern hydraulic machines thus depend to a great extent on the accuracy of the CFD results.
Accurate prediction of the flow in draft tubes is one of the most challenging tasks in CFD simulations of hydropower systems [1]. The draft tube flow at off-design is very complex, with a strong swirl in an adverse pressure gradient that leads to vortex breakdown [15,16]. The initiation of the vortex breakdown is highly dependent on the location, velocity distribution, and intensity of the vortex entering the draft tube. The vortex breakdown produces an unstable helical vortical structure, known as Rotating Vortex Rope (RVR). The RVR causes large pressure and force pulsations that could seriously affect the lifetime of the turbine [17]. The importance of reliable computations of the complex draft tube flow field led researchers to establish workshops and benchmark cases, such as the Flow Investigation in Draft Tubes (FLINDT) project [18], the Turbine-99 workshops [19], and the Francis-99 workshops [20]. The outcomes of the workshops showed significant variations in the pressure recovery results produced by the participants [1]. Such large variations can be explained by the numerous uncertainties in the CFD simulations. Modeling a real-world phenomenon using CFD often requires many assumptions that impose uncertainties in the calculations. Verification and validation (V&V) procedures should be used to reduce such uncertainties and to investigate the sensitivity of the predicted or observed flow [21,22].
In most cases, the numerical and experimental results differ to some extent, and engineering experience is used to assess the similarity of the results. However, to further improve the validation process and to better understand the results it is important to investigate the sources of the observed differences. Numerical predictions and experimental observations of hydraulic machinery flows contain uncertainties that arise from different sources. Those can be grouped into aleatory uncertainties (such as randomness in physical properties, geometry, and operational conditions) and epistemic uncertainties (such as inaccurate numerical set-up, or incorrect formulation of turbulence models). Accurate and reliable CFD computations of hydraulic machines require an assessment of the effects of different uncertainties on the system output. In recent years, several authors have studied the effects of different uncertainties on wind turbines [23e27] and gas turbine components [28e31]. Recently, Salehi et al. [32] considered the effects of geometrical and operational uncertainties in a hydraulic machine. The impeller blade shape, volumetric flow rate, and rotational speed of a centrifugal pump were considered uncertain. It was shown that the assumed uncertainties have considerable effects on both the flow field and the performance of the pump. The head variation was observed to be significant while the efficiency of the pump had a robust behavior with respect to the uncertainties. Although the flow field in hydraulic turbines is completely different than in hydraulic pumps, the experiences gained for the UQ analysis are used here.
A high-head Francis turbine test case was presented in the Francis-99 workshop series [20] with experimental data measured at steady and transient operating conditions for validation of numerical simulations. Trivedi [33] performed a V&V study on the Francis-99 test case at steady design and off-design conditions and provided approaches that could be used in V&V of hydraulic turbine simulations. In the test case, two components of the velocity field, i.e., axial and horizontal velocities at certain locations of the draft tube are available from experiments. Therefore, this case could potentially be used for verification of the numerical prediction of swirling draft tube flow. A number of researchers have utilized this test case and tried to predict the velocity field inside the draft tube. However, most of them either did not present and compare the horizontal component with the experiments [34e38], or reported large differences between the numerical and experimental results [39,40].
Recently, Salehi et al. [41,42] presented a detailed numerical analysis of the transient flow field of the Francis-99 turbine during the shutdown and stratup sequences. Different physical aspects of the flow field, such as pressure fluctuations, force analysis, velocity field variation, draft tube vortical structures, etc., were thoroughly analyzed. It was also argued that the extraction of the velocity field near a vortex center could be very sensitive to the PIV plane positioning.
The present work suggests an explanation of the frequently observed deviation between the numerical and experimental velocity field in the draft tube of the Francis-99 turbine. The issue is addressed through an uncertainty quantification and sensitivity analysis at both steady and transient (shutdown) operations. Here, the probe locations of the experimental data are considered stochastic with reasonable uncertainties, and the numerical simulations are sampled through these stochastic probes. The paper describes how sensitive the experimentally and numerically observed properties of the draft tube vortex are to the probe locations used for the data extraction.
The remainder of the current paper is organized as follows. The investigated turbine model, as well as the assumed uncertainties are introduced in Section 2. The governing equations and computational aspects of the CFD simulations are explained in Sections 3 and 4, respectively. A mesh sensitivity study is performed in Section 5. Section 6 describes the utilized probabilistic framework for uncertainty quantification and sensitivity analysis. Section 7 presents and discusses the statistical results. Lastly, a conclusion of the paper is provided in Section 8.

Investigated case
The high-head Francis turbine model of the Francis-99 workshop series [20] is chosen as a test case and described here.

Turbine model
Francis-99 is a 1:5.1 scaled-down model of a high-head turbine prototype. Fig. 1 displays the model at two sections, namely ynormal and z-normal. The model assembly is comprised of four different regions, i.e., spiral casing, guide vanes, runner, and draft tube (not shown). The spiral casing contains 14 stay vanes and the guide vane region has 28 blades. The runner has 15 blades and 15 splitters.

Available experimental data
During the Francis-99 workshops, a number of measurement locations were determined to compare data between numerical and experimental studies. The velocity measurements were performed using the Particle Image Velocimetry (PIV) flow visualization technique. Therefore, a PIV plane was established in the conical part of the draft tube in which two in-plane velocity components were measured. The PIV plane is shown as a red line in the z-normal section (Fig. 1a) and a gray area in the y-normal section (Fig. 1b). The velocity components were reported on three lines along the PIV plane, which are illustrated in Fig. 1b. Lines 1 and 2 are horizontal, while Line 3 is a vertical line at the draft tube center.
The available measured velocity components are horizontal (U) and axial (W) velocities. The horizontal coordinate is parallel to Lines 1 and 2, similar to the tangential shift shown in Fig. 1a towards the À y direction, whereas the axial direction is the same as the z coordinate (positive axial is upward direction). The normal direction is also normal to the PIV plane, similar to the normal shift shown in Fig. 1a towards the þ x direction.

Uncertainties in the PIV measurements
We have seen in our investigations of the current test case that the experimental validation of the numerical velocity field can be very sensitive to the positioning of the PIV plane. It is therefore expected that some differences between the experimental and numerical results can have their origin in the extraction of the data rather than in the results. Accordingly, in the present paper, the PIV plane positioning is assumed uncertain. A finite plane in 3D space is positioned with six degrees of freedom (three translational and three rotational). The utilized PIV measurement plane has a thickness of 3 mm. This is here used as a reasonable uncertainty in the translational positioning of the PIV plane in the tangential, normal, and axial directions, as shown in Fig. 1. Additionally, a small uncertainty with a standard deviation of 0.5 is considered for the plane Euler angles 4, q, and j, which define the orientation of the plane in 3D space. The angles are shown in Fig. 1. The assumed uncertainties in the PIV plane position are summarized in Table 1. A uniform Probability Distribution Function (PDF) with a specific mean (m) and standard deviation (s) is considered for all uncertain parameters. Hence, based on the integral definition of the standard deviation, the uncertain variable change in the range of ½m À s The data from the numerical results are extracted with these uncertainties, and the sensitivity to the uncertainties is analyzed.

Operating condition
The Francis-99 model has been experimentally studied under different steady operating conditions and during transient sequences. The steady operating conditions are at Best Efficiency Point (BEP), Part Load (PL), and High Load (HL). The model net head at the best efficiency point is H z 12 m. The transient sequences are during load acceptance (BEP to HL), load rejection (BEP to PL), shutdown (from BEP), and startup (to BEP). The guide vane angle, a, which is measured with respect to the completely closed position, may change between 0.80 and 12.43 , corresponding to minimum load and high load, respectively.
In the present study, the transient shutdown sequence is considered. The shutdown sequence is one of the most challenging and complex sequences in terms of both numerical simulations and fluid flow physics. In the current test case, the shutdown sequence starts with 1 s of the steady-state BEP operating condition, continues with 7 s of load reduction, and ends with 2 s of the minimum load condition. During load reduction, the guide vanes rotate (close) around their axes and the flow rate consequently decreases.
The guide vane angles change from a ¼ 9.84 to a ¼ 0.8 , while the flow rate changes from Q ¼ 0.19959 m 3 /s to Q ¼ 0.022 m 3 /s. The runner is rotating at a constant rotational speed of U ¼ 333 rpm during the whole shutdown sequence.

Governing equations
The Reynolds-averaged equations of transient incompressible turbulent flow are given by where Àru ī u j is the unknown Reynolds stress tensor. In the current study the SST-based Scale-Adaptive Simulation URANS model (SST-SAS) [43,44] is employed for the calculation of the Reynolds stress tensor. The formulation of this model differs from the k À u SST RANS model [45] mostly by introduction of an additional source term in the turbulence eddy frequency (u) transport equation. The inclusion of this source term enables the model to resolve the turbulent spectrum and eddies of the flows containing transient turbulent instabilities and provide LES-like flow solutions. The model has been utilized and verified in numerous studies on hydraulic turbine flows [33,35,46e50].

Numerical aspects
A brief overview of the numerical aspects of the conducted CFD computations is provided in this section. For more details on the employed methodologies and schemes, please see Salehi et al. [41].
The CFD computations of the present investigation are performed using the OpenFOAM-v1912 open-source CFD solver [51,52], which employs the finite volume method on a collocated grid. The numerical aspects, boundary conditions, and mesh motion approach are here discussed.

Discretization schemes and pressure-velocity coupling
The temporal derivatives are discretized using the implicit second-order backward scheme. The time step is chosen as Dt ¼ 1.25 Â 10 À4 s, which corresponds to 0.25 of runner and 1.625 Â 10 À4 of guide vane rotations. The time step is sufficiently small to resolve the most important flow structures and frequencies (e.g. the propagation of the runner blade wakes).
The convective terms in the momentum equation are discretized using the Linear-Upwind Stabilised Transport (LUST) scheme [53], which blends the central and second-order upwind schemes. A blending factor of 0.75 is employed which means that the face values in the momentum equation are calculated blending 75% of the second-order central and 25% of the second-order upwind schemes. This blended scheme provides better numerical stability while maintaining second-order accuracy. Other convective terms (i.e. for k and u) are approximated using the secondorder upwind discretization scheme.
The pressure field is solved through the PIMPLE algorithm, which combines the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) [54] and the Pressure Implicit with Split Operator (PISO) [55] algorithms. In fact, PIMPLE employs the SIM-PLE algorithm as an outer loop correction while the PISO correction is used as an inner. Ten outer and two inner correction loops are considered in each time step.

Boundary conditions
The guide vane movement for the transient case is applied through a newly developed mesh motion boundary condition that needs the time evolution of the rotational speed and the axis of rotation for each guide vane as input. The guide vane rotational speed is given by the derivative of the guide vane angle (Fig. 2a), shown for the present transient case in Fig. 2b. A smooth transition is applied at the start and stop of the guide vane motion to dampen the pressure surge effects that arise due to the fully incompressible formulation. As seen in these figures, from t ¼ 0 s to t ¼ 1 s, the stationary BEP condition is maintained, where the guide vanes are kept still and the flow rate is fixed at Q BEP ¼ 0.19959 m 3 /s. Also, in order to establish a statistically stationary condition, a 2 s flow time simulation is performed beforehand.
The volume flow rate of the turbine is assumed to vary linearly along with the change of the guide vanes angle, with a similar smooth transition. A time-varying uniform velocity, according to flow rate, is imposed at the inlet of the spiral casing. The inlet values of the turbulent kinetic energy k and the specific rate of dissipation u are also varying with time assuming a fixed turbulence intensity (I ¼ 10%) and viscosity ratio (n t /n ¼ 100). The pressure has a zerogradient assumption at the inlet. At the outlet boundary, the zero-gradient condition is imposed for U, k, and u, while the pressure is assumed to be fixed. A special treatment for the outlet boundary is employed, which prevents backflow. The simulation convergence was monitored throughout the whole transient sequence using inlet and outlet flow rates. For instance, the relative error between the inlet (imposed) and the outlet (calculated) flow rates at BEP was around 0.01%.
A no-slip boundary condition is applied for the velocity at all walls, where the turbulent kinetic energy and specific dissipation rate are calculated using wall functions. The wall pressure is subjected to the zero-gradient assumption. The Cyclic Arbitrary Mesh Interface [56,57] boundary condition is employed for interfaces between the spiral casing, guide vane, runner, and draft tube regions.

Dynamic mesh approach
In the transient operation of Francis turbines, the guide vanes rotate around their own individual axes while the runner is rotating around the turbine axis. Hence, the numerical simulation involves the simultaneous mesh deformation due to guide vane movement at the same time as the runner part of the computational domain is rotating as a solid body. A solid body rotation is applied to the runner region, while a Laplacian displacement mesh morphing solver is used for the mesh deformation of the guide vane region.
In OpenFOAM, mesh morphing is handled through the Laplacian smoothing equation where G is the motion diffusivity and d cell is the displacement vector of the cell centers. Having obtained the solution for d cell , the point displacement d point is calculated through interpolation and consequently the location of new points are computed with More information on the dynamic mesh properties, as well as the open-source case and codes are provided by Salehi and Nilsson [58].

Computational domain and mesh study
The full Francis-99 model turbine is considered in this study, i.e., from the inlet of the spiral casing to the outlet of the draft tube. The computational domain is divided into four regions, namely the spiral casing, guide vanes, runner, and draft tube (see Fig. 1). The corresponding mesh of each region is generated separately. A block-structured mesh is generated for the spiral casing and the draft tube in ICEM-CFD. The guide vane and runner regions are first meshed with block-structured O-grids for a periodic section (containing one blade for the guide vane region and one main blade and one splitter for the runner region), using Turbo-Grid and ICEM-CFD, respectively. The periodic sections are then copied and merged to produce the full guide vane and runner region using ICEM-CFD. The four regions are then merged together in OpenFOAM, constructing the full computational domain. The coupling between the different regions is realized using the cyclicAMI interface.
Three computational meshes with different resolutions, named Coarse, Medium, and Fine, are created to assess the mesh independence at the BEP sequence of the simulation, see Table 2. The impact of mesh resolution on the time-averaged torque of the shaft is also presented in the same table. It is clear that increasing the mesh density improves the accuracy of the shaft torque, compared to the experimental value. Although all meshes have predicted the torque accurately, the minimum error (with respect to the experimental value) corresponds to the finest mesh. The overall view of the turbine assembly is shown in Fig. 4a, while Fig. 4b and c demonstrate close views of the guide vane and runner meshes, respectively. The guide vanes close down during shutdown and their meshes deform. The complete mesh assembly contains almost 16 million cells. The details of the mesh parameters are presented in Table 3, including the maximum and average values of the y þ value at the start and end of the transient simulation.

Probabilistic framework
The employed probabilistic framework for statistical analysis is described here. Section 6.1 explains the details of the uncertainty quantification procedure using the polynomial chaos expansion method, while the details of the sensitivity analysis using Sobol' indices are given in Section 6.2. The numerical details of the current UQ and SA computations are provided in Section 6.3.

Uncertainty quantification using polynomial chaos expansion
Suppose y ¼ UðxÞ is a physical or mathematical model, where x is uncertain, with joint probability distribution function f(x), the model response y is stochastic. Using homogeneous chaos theory [59], it has been shown by Soize and Ghanem [60] that if the stochastic input variables x are independent, the model response can be represented as a series of orthogonal polynomial basis functions. The polynomial chaos representation of a random field of order p for d random variables x≡fx i g d i¼1 can be written as where the a 0 s are n-tuples in N d 0 dfða 1 ; …; a d Þ : a j 2N ∪f0gg. The number of terms in the summation is The u a 's are unknown coefficients and the j a (x)'s are multivariate polynomials orthogonal with respect to the joint probability distribution function f(x), i.e., Cj a ðxÞ; j b ðxÞD ¼ Cj a ðxÞ; j a ðxÞDd ab : where C.,.D represents the inner product in the Hilbert space and d is the Kronecker delta. The input random variables x≡fx i g d i¼1 are assumed to be independent, so that the joint probability distribution function can be expressed as where f i (x i ) is the PDF of x i . Each basis function j a (x) is constructed using a tensor product of univariate polynomials. In order to achieve exponential convergence, the basis of the Polynomial Chaos Expansion (PCE) should be a set of polynomials that are orthogonal with respect to the PDF of the input uncertain parameters. These distributions and corresponding polynomials are presented by Xiu and Karniadakis [61]. In the present study, all input uncertain parameters are assumed to have uniform distributions. Hence, Legendre polynomials are employed to construct the orthogonal basis.

Calculation of PCE coefficients
The orthogonal basis of the PCE is chosen based on the input PDFs. Therefore, if the PCE coefficients (u a 's) are known, the PCE of the response is also known. Thus, the problem is to find the coefficients. In the current study, the regression method has been used to compute the response expansion coefficients. The regression-based Non-Intrusive Polynomial Chaos (NIPC) method [62] starts with Eq. (5). First, one should generate N realizations to fill the input stochastic vector X ¼ {x (1) , x (2) , …, x (N) } and evaluate the QoI for each realization Y ¼ fy ð1Þ ; y ð2Þ ; …; y ðNÞ g T , where y ðiÞ ¼ Uðx ðiÞ Þ. The problem is computationally ill-conditioned and the number of realizations is commonly more than the number of unknowns for numerical stability. Following Hosder et al. [62], N ¼ 2(P þ 1) vectors may be chosen in the stochastic space and the stochastic function is evaluated at these sampling points using the   deterministic solver. The coefficients can then be obtained from the solution of the over-determined linear system given by using the least-square approach where J is the measurement matrix In the current study, the Sobol' quasi-random sequence [63] is used to generate the sample points. The Sobol' sequence is a base-2 digital sequence that fills space in a highly uniform manner.

Postprocessing of PCE
Due to the orthogonality of the basis, the mean and variance of the QoI read Hence, if the PC basis is orthonormal, the variance will be easily calculated by summing the squares of PCE coefficients except for the first coefficient.

Sensitivity analysis using Sobol' indices
Sensitivity Analysis (SA) methods aim at quantifying the respective effects of input random variables (or combinations thereof) on the variance of the response of a system. Among the abundant literature on sensitivity measures, the Sobol' indices have received much attention because they provide accurate information for most models [64]. Each Sobol' index S i1;…;is is a sensitivity measure that describes which amount of the total variance is due to the uncertainties in the set of input parameters {i 1 , …, i s }. Using Sobol' decomposition of PCE [64,65] into summands of increasing dimension, it can be shown that the PC-based Sobol' indices are defined as S i1;…;is ¼ 1 where I i1;…;is is the set of a tuples such that only the indices (i 1 , …, i s ) are nonzero, i.e. I i1;…;is ¼ a : a k > 0 ck ¼ 1; …; n; k2ði 1 ; …; i s Þ a k ¼ 0 ck ¼ 1; …; n; k;ði 1 ; …; i s Þ : In addition, the total PC-based sensitivity indices read S T j1;…;js ¼ X ði1;…;isÞ2I j 1 ;…;js S i1;…;is : (17)

UQ and SA computations
All the UQ and SA computations were carried out through an inhouse uncertainty analysis program, which has been previously developed and validated [32,66e70]. First, the CFD computation is conducted and then the UQ and SA calculations are carried out as a post-processing step.
As discussed in Section 2.3, the positioning of the PIV measurement plane is assumed uncertain with respect to six stochastic variables (d ¼ 6) that are uniformly distributed over certain ranges (see Table 1). A second-order Legendre polynomial basis (p ¼ 2) is utilized for the construction of the PCE. Therefore, considering over-sampling of n p ¼ 2, 2  (13) and (14). When the PCE coefficients are known, the PDF of the response could also be computed using the calculated PCE. Here, all the presented PDFs are constructed using a large number of samples (10 6 ) on the PCE.

Results and discussion
The current section is devoted to the uncertainty quantification and sensitivity analysis results. Firstly, the results during the steady phase at BEP are presented and discussed in Section 7.1 and then the results during the transient shutdown sequence are described in Section 7.2.

Steady at best efficiency point
It was described in Section 4.2 that the shutdown sequence starts with BEP condition for 1 s flow time (see Fig. 2). The numerical results for the BEP condition are time-averaged throughout this 1 s stationary period.
Time-averaged horizontal (U), axial (W), and normal (V) velocity components from the deterministic simulation and UQ analysis at BEP are presented and compared to measured PIV data in Fig. 5. It should be recalled that the location of the PIV measurement plane and its corresponding uncertainties were previously described in Section 2 (See Fig. 1 and Table 1). The horizontal and normal velocity components are tangential and normal components to the PIV plane, and the axial direction is the same as the z-direction. The contours of normalized two-dimensional PDFs (PDF/max(PDF)) are presented alongside the mean and deterministic (without considering uncertainties on the PIV plane positioning) values. The available experimental measurements are also presented in Fig. 5 for comparison. The normal velocity component to the plane (V) was not measured in the experiments. The experimental data are reported at BEP for 200 different snapshots in time. Therefore, the ensemble average and standard deviation (StD) of the experimental data are calculated and presented as ± s bounds. It is seen that the experimental horizontal velocity (U) has a larger standard and 0.065, respectively. This could be due to the fact that the horizontal component has much lower values than the axial. Therefore, compared to the axial component, the measured horizontal velocity is much more sensitive to turbulent flow fluctuations, as well as different experimental errors (e.g. measurement accuracy and human errors). Additionally, Line 3 is located at the center of the draft tube cone downstream of the runner crown and thus inside the recirculation region. Hence, the standard deviation of the velocity components, especially the horizontal velocity, is larger on Line 3. The numerical mean (m, calculated through Eq. (13)) and deterministic velocity profiles are quite similar, which can be explained by the fact that all stochastic variables are assumed to have symmetrical PDFs where the mean values are identical to the deterministic ones. The PDF bounds of the numerical horizontal velocity are much larger than for the other components. The average CoV of the different velocity components on the three PIV lines is presented in Table 4. The average standard deviations of the horizontal and normal velocities are moderately larger than for the axial velocity. Thus, the U and V components are more affected by the considered uncertainties. However, when the components are normalized by their corresponding mean values, the CoV of U is significantly larger which is mostly due to small mean values, whereas W shows a negligible coefficient of variation. For instance, the average CoV of the U, W, and V velocities on Line 3 are 1.022, 0.0035, and 0.5512, respectively. Accordingly, the horizontal velocity is extremely sensitive to the uncertainty in the PIV plane positioning. The normal velocity has a smaller CoV and is moderately affected by the uncertainties, while the assumed uncertainties  The comparison of numerical and experimental results shows similar trends for both horizontal and axial velocities. While the deviation between the axial components is small, that of the horizontal components is rather large compared to the actual values. It is one of the purposes of the present work to further investigate this difference, which has not yet been sufficiently done in the literature. The decrease of axial velocity close to the rotational axis, caused by the wake behind the flat surface of the runner at the center, is accurately predicted by the computations. On the other hand, for the horizontal velocity, the deviations between the experimental and numerical data are larger. In fact, an accurate prediction of the horizontal (or radial) velocity in the draft tube is more challenging and requires a precise capturing of the location and intensity of the draft tube vortex. As previously discussed, most researchers that studied the current test case either did not compare the radial component with the experiments [34e38] or reported large differences between the numerical and experimental results [39,40]. The uncertainty analysis in the present paper helps to understand and address this large deviation.
A centered vortex would produce zero radial velocity on the midpoint of the horizontal lines. The non-zero horizontal velocity at the center of Lines 1 and 2, and along the entire Line 3, indicates a slightly off-centered vortex in both the numerical and experimental results. The numerical results show a time-averaged vortex that is more centered at the PIV plane. The fact that the vortex in the experiments is more off-centered explains the high standard deviation of the U velocity observed at the center in the measurements. Even though there are some differences between the numerical and experimental values of time-averaged horizontal velocity, interestingly the numerical PDF bounds are always covered by the experimental standard deviation. Consequently, the discrepancies between the numerical and experimental U components can be partly explained by the uncertainties in the experimental measurements (for instance the uncertainties in the PIV plane positioning). One should note that a small standard deviation is considered for the input random variables in the present study (0.5 for the angles and 1.5 mm for the shifts). Nevertheless, the standard deviation of the horizontal velocity is significant. It is seen in Table 4 that the average standard deviations of the U component on Lines 1, 2 and 3 are 0.024, 0.022, and 0.042 m/s, which corresponds to coefficients of variation of 2.78, 4.76, and 1.02, respectively. This comparison highlights the importance of considering uncertainties in the prediction of the horizontal (radial) velocity in turbine draft tubes.
A sensitivity analysis using Sobol' indices is performed to assess the importance of each input stochastic variable. The total Sobol' index of a random variable specifies its total contribution to the standard deviation of the stochastic response. The total Sobol' indices of all considered uncertain parameters on the different velocity components at the three lines are plotted in Fig. 6. Here, the curves are smoothened separately for better illustration. It is obvious that the effect of the 4 angle and the axial shift on the standard deviation of all three velocity components at the different PIV lines are negligible. This means that the gradient of the velocity components with respect to the circumferential and axial directions is small. For each velocity component, similar trends are seen for Lines 1 and 2, while Line 3 shows different behavior, which is expected because it is aligned in another direction.
The j angle (rotation around the normal vector of the plane, see Fig. 1) is mostly responsible for the uncertainty of the horizontal velocity, except at the center where the normal shift plays the most important role. This could be explained by the fact that the horizontal velocity at the center is mainly affected by the location of the secondary flow vortex core. Therefore, the uncertainty in the normal shift influences the distance between the PIV line centers and the vortex core and hence changes the horizontal velocity. This effect is also visible in the S U Sobol' index on Line 3, which is located at the center of the draft tube. Accordingly, the normal shift is also the most important uncertain variable in the standard deviation of the horizontal velocity on Line 3. After the normal shift, the j angle has the highest impact on the standard deviation of the horizontal velocity because it changes the direction of the U unit vector.
Moreover, the q angle effects are only visible at both ends of Line 3 where the PIV plane location is more influenced by this uncertain variable.
The second row of Fig. 6 displays Sobol' indices of the axial velocity (W). The tangential shift is the main stochastic parameter that determines the uncertainty of the W velocity on Lines 1 and 2. This can be explained by the sharp variation of the axial velocity along the line (tangential direction) in Fig. 5. The normal shift parameter only affects W in the middle of the line. Here again, the normal shift varies the distance between the line center and the vortex core. Although it is seen that at the top of Line 3 the tangential shift is the dominant variable and then it decreases along the line, one should note that the standard deviation (or CoV) of the axial velocity along this line is small.
Finally, the Sobol' indices of the normal velocity (V) are plotted in the third row of Fig. 6. Here the q angle is most significant on Lines 1 and 2, mainly because it changes the direction of the normal velocity unit vector. The importance of the tangential shift at the center and also on Line 3 can be explained by the same argument about the distance to the vortex core.

Shutdown sequence
After the BEP steady-state phase, the turbine proceeds with the transient shutdown sequence. The guide vane closure initiates and thus the flow rate decreases until the minimum load condition and then the turbine works at this condition for 2 s (See Fig. 2). The present section provides and discusses the UQ and SA analysis of the flow field during the turbine transient shutdown sequence. The 56 realizations created for the PIV plane are sampled throughout the entire shutdown sequence and statistical moments, PDFs, and Sobol' sensitivity indices are computed employing the UQ and SA methodology described in Section 6. For each probe location the UQ analysis is performed for all time steps of the transient sequence and the time-variation of the statistical moments, namely, mean (m, Eq. (13)), standard deviation (s, Eq. (14)), and Sobol' indices (S, Eqs. The U velocity (Fig. 7) has a small instantaneous value and is mostly fluctuating around zero. The deterministic PIV plane is located at the center of the draft tube cone. As previously described, if the center of the vortex is positioned at the draft tube center, the radial velocity component (here called horizontal) should be zero on both Lines 1 and 2. Therefore, the small oscillating horizontal velocity at the center of Lines 1 and 2 at the BEP condition suggests a presence of a weak and slightly off-centered swirling flow. When the turbine goes through the shutdown sequence, a triangle-shaped fluctuating zone develops in both the numerical and experimental U velocity distributions after 2 s, which indicates an expansion of the off-centered vortex.
The axial velocity component, W, (Fig. 8) indicates that a reversed flow region is created with the start of the transient process at t ¼ 2 s. The region is first formed at the center of the draft tube and then develops in time until it covers the whole PIV lines.
The comparison of the statistical UQ mean results with the measured data suggests that the trends of both the U and W velocities are adequately predicted by the present numerical simulation. The comparison outcomes are compatible with the steady results shown in Fig. 5. As explained in Section 7.1, the U component is quite sensitive to the assumed uncertainties and thus the deviations between the numerical and experimental results are large. The axial velocity component demonstrates a more robust behavior with respect to the uncertainties, and hence the experimental and numerical results of W are more similar.
No experimental measurements are available for the normal velocity component, V, which represents the tangential (or swirl) velocity (with a switch of the sign at the center). The sign of the normal velocity on Lines 1 and 2 at the BEP condition (t ¼ 0e1 s) is negative for 0 < S/S max < 0.5 and positive for 0.5 < S/S max < 1, indicating the presence of a counter-rotating vortex at BEP. Moreover, the magnitude of the non-zero V velocity at BEP is not significant and therefore the vortex is not very strong. However one can see that when the turbine load decreases, after t > 2 s the direction of the swirl changes, and a vortex in the same direction as the runner develops in the draft tube. In fact, after t > 2 s, the sign of V velocity becomes positive and negative for 0 < S/S max < 0.5 and 0.5 < S/S max < 1, receptively. The strength of the vortex amplifies in time. It has been recently shown that the swirl number (swirl intensity) in the draft tube of Francis-99 significantly grows during load rejection operation [71]. The enhancement of the swirl velocity is explained with the velocity triangles of the runner outlet, where the swirl component is increased with a flow rate reduction. The fluctuations seen for the normal velocity between 4 and 6.5 s show a periodic change of sign, implying that the location of the center of the vortex is changing in time. Large oscillations are noticed in all three velocity components between t ¼ 4 s and t ¼ 6 s when the RVR is formed. For more details on the variation of the velocity field and the mechanism of formation and collapse of the RVR during the shutdown sequence, readers are referred to Salehi et al. [41].
Performing the previously described UQ analysis, the standard deviation and consequently the coefficient of variation of the QoI can be obtained. The time-variation of the StD and CoV of the three velocity components on Line 1 are contoured in Fig. 10. Similar results for Lines 2 and 3 are not shown for briefness. Before the transient sequence starts (t ¼ 0e1 s), a fluctuating high standard deviation region is observed around the center of the PIV lines (i.e., the center of the draft tube). In fact, there is a slender counterrotating vortical structure around the center of the draft tube at the steady BEP condition (the vortical structures are later shown using the l 2 criterion in Fig. 15). The presence of such a vortex makes the velocity field sensitive to the assumed uncertainties in the sampling plane position. However, the strength of this vortex at BEP is not significant. As seen in Fig. 9, after t ¼ 2 s, a strong vortex in the same direction as the runner forms at the draft tube center. The strong central vortex considerably increases the standard deviation of all the velocity components near the center. The vortex temporarily enlarges until t ¼ 3 s which results in the development of a V-shaped high standard deviation region. Hereafter, the central vortex becomes unstable and a reversed flow region develops which results in the decrease of velocity standard deviation near the center.
An analysis of the CoV provides additional information as it clarifies the significance of the StD compared to the mean. As seen in Fig. 10, the CoV of the U velocity at BEP is the largest, due to the presence of a nearly axisymmetric flow with a small fluctuating radial component. However, in the part load condition, the unstable vortical structures in the draft tube provide a non-zero horizontal velocity along Line 1 and therefore the CoV decreases. The CoV of the axial velocity component (W) is only significant at the edge of the developing reversed region (also seen in the mean contours in Fig. 8), where the mean value is small and consequently the CoV is large. The CoV of the V velocity is remarkable around the center of the draft tube at the BEP condition. When the turbine goes into the transient shutdown sequence, the swirl direction changes, and the normal velocity momentarily becomes zero. This can be clearly seen as a high values region in the contours of CoV around t ¼ 2 s. After that (for t > 2 s) the CoV is first large around the center due to the formation of the central vortex. When the turbine reaches the part load condition the vortex becomes unstable and moves in a helical pattern in the draft tube.
The impact of the considered uncertainties on the transient velocity field at different time instances is further analyzed in Figs. 11 and 12. The figures illustrate normalized 2D PDFs of different velocity components of Lines 1 and 2 compared to the UQ mean, deterministic curves, and experimental data. Four different time instances of the shutdown sequence are presented, namely, t ¼ 1.5 s, t ¼ 2.0 s, t ¼ 4.0 s, and t ¼ 6.0 s, corresponding to guide vane opening angles of a ¼ 9.35 , a ¼ 8.70 , a ¼ 6.10 , and a ¼ 3.50 and turbine loads of 95.09%, 88.55%, 62.38%, and 36.21%, respectively.
As previously observed, the horizontal velocity component exhibits larger PDF bounds indicating that the U velocity is more sensitive to the PIV plane positioning uncertainties. At t ¼ 2 s the standard deviation of U on both lines is significant around the center, denoting the formation of a strong central vortex, whereas at t ¼ 4 s the vortex has momentarily shifted to the right side of Line 1 and the PDF bounds are more significant in that region. Here again, the mean and deterministic curves are similar for all velocity components at different times. Generally, the deviation between the numerical results and the experimental measurements during the transient shutdown sequence is more noticeable compared to the BEP results presented in Fig. 5, especially for the horizontal velocity. This is explained through the random turbulent fluctuations at low load conditions. The highly swirling turbulent flow produces flow instabilities that result in random local fluctuations. This is in accordance with a previous study by Goyal et al. [72], which performed several experimental repetitions of transient sequences (load rejection and load acceptance) on the same test     case (Francis-99) and showed that during PL the velocity components are not satisfactorily repeatable and vary within a remarkable range.
The sudden increase in the standard deviation seen in Fig. 10 is more clear in Fig. 13a, where the time-variation of the UQ results are presented for the midpoint of Line 1 (the point location is shown in Fig. 1b). At the start of the shutdown sequence, the StD of the velocity components on the Line 1 midpoint drastically rises from 0.05 m/s at t ¼ 1.5 s to more than 0.2 m/s at t ¼ 2 s and then it reduces promptly. The increase is remarkably more pronounced for the U velocity when it is normalized and presented as CoV in Fig. 13b, due to its small mean values. Hereafter, some local peaks are seen in the StD of the velocity components (e.g. at t ¼ 3.5 s and t ¼ 6 s) that indicate a presence of local vortical structures that are temporarily located around the center of the draft tube cone.
In order to further study this effect (sudden rise in StD), the surface streamlines on the z ¼ À0.3386 m ¼ À0.976D plane of the draft tube cone, where Line 1 is positioned, are shown in Fig. 14.  Line 1 is seen in these pictures as a white strip. The surface is also colored by the secondary flow magnitude (magnitude of the inplane velocity vector, ). The change in the swirl direction, which was observed and explained in Figs. 9 and 10, is clearly visible in the streamlines between t ¼ 1.6 s and t ¼ 2.2 s. Line 1 passes through the draft tube center and thus as previously argued, the counter-rotating vortex at BEP is slightly off-centered at BEP (t ¼ 1 s, Fig. 14a). Although a snapshot at BEP is shown in this figure, it was previously shown that even the time-averaged vortex is not perfectly centered. When the turbine load decreases, a vortex originates at the draft tube center and its strength increases over time (See Fig. 14c). Then after t ¼ 2.2 s, the central vortex becomes unstable and moves from the center. The abrupt rise in the StD of the midpoint velocity, observed in Fig. 13, can also be explained in this figure. It is seen that after t ¼ 2 s a strong secondary flow is formed at the center of the draft tube which makes the midpoint velocity components more sensitive to the uncertainties in the PIV plane position. However, after t ¼ 3 s this vortex becomes unstable and starts moving around while the secondary flow magnitude at the center decreases drastically and a central reversed flow region is formed (Fig. 14def).
Furthermore, Fig. 15 displays the draft tube vortical structures and reversed flow region using iso-surfaces of l 2 ¼ 750 s À2 and W ¼ 0, respectively. The l 2 value was selected to get the best visualization of the vortical structures. When the turbine moves into the part load condition, the residual tangential velocity leaving the runner increases. Consequently, a central vortex structure forms and swells (Fig. 15aec). The development of this strong and straight vortex, between t ¼ 1.5 s and t ¼ 2 s, explains the sudden increase in StD of the velocity components in Figs. 10 and 13. The presence of such a strong and straight vortex amplifies the velocity gradient and thus increases the sensitivity of velocities to the probe location. Thereafter, a reversed flow region initiates at the center of the draft tube, demonstrated by the red color in Fig. 15def. Accordingly, similar to previous experimental studies in the literature [73], the central vortex becomes unstable and collapses into different smaller vortical structures that are wrapped around the central reversed flow region and thus the velocity StD reduces. A Sobol' sensitivity analysis is carried out on all time steps throughout the shutdown sequence and the time-variation of the total Sobol' indices of different velocity components on Line 1 are presented in Fig. 16. The total Sobol' indices of the 4 angle and axial shift uncertain variables are omitted here as their effects are negligible throughout the entire transient operation. The j angle is largely accountable for the standard deviation of the U velocity at BEP and the initial part of the transient procedure (t < 3 s). After the creation of the wide reversed flow region and formation of the rotating vortex rope, the Sobol' index of normal shift, which affects the distance between PIV lines and the vortex center, becomes dominant. The tangential shift is the most significant Sobol' index for the W velocity, which is due to the steep variation of the axial velocity along the PIV line. The q angle that changes the direction of the normal velocity unit vector is mostly responsible for changes in the V component at BEP and the initial part of the transient operation. However, in the part load condition, where the draft tube vortex is unstable the tangential shift becomes the dominant uncertain variable due to the fact that it changes the distance to the vortex core.

Conclusion
An in-depth analysis of the impact of the uncertainties in the positioning of an experimental PIV plane on the validation of the numerical results was presented. The high-head Francis-99 turbine model during shutdown sequence was considered as the investigated test case. The flow computations were carried out employing the OpenFOAM-v1912 open-source CFD code. The uncertainty quantification was conducted using the polynomial chaos expansion method employing second-order Legendre polynomials, while the Sobol' indices were utilized for the sensitivity analysis. The positioning of the experimental PIV plane in the draft tube cone was assumed uncertain with three translational and three rotational stochastic parameters with uniform probability distribution functions.
At the best efficiency point condition, the horizontal velocity profiles in the draft tube were highly sensitive to the uncertainties in the PIV plane positioning. The average CoV of the U component on Lines 1, 2, and 3 was 2.78, 4.76, and 1.02, respectively. On the other hand, small PDF bounds of the axial velocity on all three lines with average CoV of 0.0058, 0.0059, and 0.0045 showed a negligible impact of the considered uncertainties on this component. The PDFs of the horizontal velocity were covered by the experimental standard deviation bounds. Therefore, some of the discrepancies between the numerical and experimental results can be explained by the uncertainties in the PIV plane positioning and/or its relation to the position of the draft tube vortex. In a deterministic simulation and validation of similar cases, one might observe imperfect agreement between numerical and experimental horizontal velocity, and the current work showed the importance of considering the effects of such uncertainties on the validation of the results. Besides, through a Sobol' indices sensitivity analysis, it was shown that the uncertainties in the rotation around the normal vector of the plane and the normal shift are mostly responsible for the changes in the horizontal velocity, which was the most sensitive component.
During the transient shutdown sequence, the UQ mean and deterministic velocity components match, which was explained through the symmetrical PDFs of the input random variables with identical mean and deterministic values. Moreover, the UQ mean velocities acceptably agreed with the experimental data. At the beginning of the shutdown sequence, a sudden rise was observed The present work shows that the uncertainty due to the positioning of the PIV plane is largest for a strong and straight vortex that is supposed to be centered at the PIV plane. The time-averaged flow in the draft tube cone is similar to such a vortex, and the uncertainty is thus present for such data. The uncertainty is largest for the component that for the Francis-99 case is referred to as the horizontal velocity. It is the normal shift of the PIV plane that mostly influences the extracted horizontal velocity at the center of the vortex. This can be explained by the fact that the horizontal velocity has a large gradient normal to the PIV plane at the center of the vortex. The misalignment of the PIV plane and the vortex center can be due to either the positioning of the PIV plane or a slight shift in the location of the vortex center. It is important to be aware of this uncertainty when comparing the experimental and numerical results, since large differences in the time-averaged horizontal velocity may not necessarily mean that the results are very different. In fact, from the experimental data of the Francis-99 case, it is clear that the time-averaged vortex is not centered at the PIV plane, since the horizontal velocity is negative throughout all of the PIV lines. The present numerical results show a time-averaged vortex that is more centered at the PIV plane. Thus, the extraction of the numerical horizontal velocity also captures more of the details at the actual center of the vortex.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.