Abstract

In this paper, a methodology is suggested to evaluate the goodness-of-fit between observed and numerically predicted variables in curved open channel flows. The mathematical model uses the continuity and the Reynolds averaged Navier–Stokes (RANS) equations with the RNG κ-ɛ turbulence model, numerically solved with the commercial software Ansys Inc. The concept of design of experiments (DOE) together with some criteria of goodness-of-fit was used in a new coherent suggested methodology for the calibration and validation of the numerical model. The goodness-of-fit criterion was applied to secondary flow variables like helicity, intensity of the secondary current and the vortex index instead of the primary variables like the velocity field and the free surface. This methodology guarantees that the numerical model reproduces values of the response variables very near to the observed values and it is valid for the mesh independence analysis and the setup of other numerical parameters, including those related to the boundary conditions.

INTRODUCTION

Flow in curved channels is characterized by flow separation, by the presence of secondary currents along the curve and perpendicular to the channel axis generated by the deflexion of the flow stream lines and by free surface variations (Van Balen 2010). These flow structures generate complex tridimensional flow patterns that have been studied by numerical modelling and laboratory experiments, the numerical approach being the one that provides better methodological flexibility to pursue a detailed study despite the need to pursue a model calibration and validation process (Han et al. 2011; Song et al. 2012). Most of the numerical models used to simulate these types of flows are based on the Reynolds averaged mass and momentum conservation laws coupled with a turbulence closure model and in some cases, the ‘rigid lid’ approach is used to simplify the numerical solution (Bai et al. 2014).

No standard methodology has been suggested in the literature to calibrate and validate a numerical model. However, most cases reported in the literature use primary variables like the velocity field at different cross sections of the curve to compare the numerical results to laboratory experiments. Additionally, there is no general consensus on the most suitable variables to represent either the flow patterns or the statistical indicators to be used to evaluate the quality of the model calibration and validation processes. Some authors (Berger & Field 1984; Moffatt & Tsinober 1992; Stranden 2007; Ghobadian & Mohammadi 2011; Han et al. 2011; Ansys Inc. 2016) have successfully used secondary variables like the helicity of the flow, the intensity of the secondary currents and the swirl number for the calibration–validation process. To compare the model results to the laboratory experiments, two statistical methods have been traditionally used: (a) the one factor at a time (OFAT) method, which varies one model parameter at a time to analyse its influence on the results (Lye 2002; Wahid & Nadir 2013) and (b) the design of experiments (DOE) method, which is a systematic method to determine the relationship between factors (variables) of a process affecting the response variables, which has been used mainly for laboratory experiments or observed data Montgomery (2012).

The prediction capability of a model can be evaluated by using different statistical parameters which quantify the goodness-of-fit between the model results and some reference values (laboratory measurements, for example): the determination coefficient (R2); the mean absolute error (MAE); the mean absolute error in percentage (MAPE); the index of agreement (d); the efficiency coefficient of Nash and Sutcliffe (NSE); the root mean square error (RMSE) and the variation coefficient, among others. A detailed discussion of different goodness-of-fit indexes can be found in the following studies (Hsu & Kuo 1999; Legates & McCabe Jr 1999; Toledo & Koutsopoulos 2004; Alexandris et al. 2008; Jara-Rojas et al. 2009; Willmott et al. 2012; Ritter & Muñoz-Carpena 2013).

In this paper, a methodology to calibrate and validate a numerical model is presented and applied to the case of a velocity field flow in a curved channel. We use the DOE method, as a new alterative to the existent methods, to set up the numerical experiments and evaluate the goodness-of-fit between the model results and observed values by using several statistical indexes for the secondary flow variables.

METHODOLOGY

Channel geometry and velocity field experimental data

Coordinate system

The local and global coordinate systems for the velocity field are shown in Figure 1(a), where the global system is defined by the x, y, z coordinates and the local system, whose origin is located at the centre of the circle of the channel curve, is defined by t, r, z coordinates representing the tangential, radial and vertical directions, respectively, whose corresponding velocity components are Vt, Vr and Vz. The relationship between the local and global coordinate systems is obtained from Figure 1(a). The meridional velocity (Vm) is computed from the Vr and Vz components as Vm = (Vr2 + Vz2)0.5 and the direction in the r-z plane is obtained as θP = tan−1(Vz/Vr).

Figure 1

Experimental channel configuration. (a) Channel geometry, three-dimensional view and local and global coordinate systems; (b) velocity measurement grid. Adapted from Bai et al. (2014).

Figure 1

Experimental channel configuration. (a) Channel geometry, three-dimensional view and local and global coordinate systems; (b) velocity measurement grid. Adapted from Bai et al. (2014).

Channel geometry

The channel was built with transparent acrylic plates, in a horizontal bottom with a total deflection angle ϕ of 180°, a channel width W = 0.40 metre (m), a radius of curvature Rm = 1.2 m and two upstream and downstream straight channels, with Le=Ls = 2.0 m. The velocity measurements were obtained for a flowrate Q = 0.01008 m3/s and a water depth at the channel inlet he = 0.15 m (see Figure 1). For these values, a Froude number of 0.148 and a Reynolds number of 68,967 were obtained, where DH is the hydraulic diameter (DH = 4RH), V is the mean velocity at the channel inlet, RH is the hydraulic radius, g is the acceleration due to gravity and ν is the water kinematic viscosity at 298.15 °K. Velocity profiles were measured at seven transverse sections located at 0°, 30°, 60°, 90°, 120°, 150° and 180° from the curve entrance. At every cross section, the velocities were measured, using a micro acoustic Doppler velocimeter (ADV) profiler and reported by Bai et al. (2014), at grid locations shown in Figure 1(b) by Bai et al. (2014).

Experimental data

ADV measurements like those of Bai et al. (2014) are commonly reported in the literature. Although they have some uncertainties, they are reliable as they are consistent with other techniques of measurement such as particle image velocimetry (PIV) (Craig et al. 2011). Among the possible uncertainties in the measurements with ADV are those of random or epistemic origin. The most common uncertainties detected in measurements with ADV are: the Doppler noise, digital filtrations made by the instrument, the presence of anomalous values and incorrect positioning of the instrument. The operator can also introduce errors in the measurements when defining the value of the processing software's input parameters such as the sampling rate, ping algorithm, power level, velocity range and collection time. There are several methods to diminish these errors (Craig et al. 2011). The data collected closer than 5 mm from the bottom were disregarded as the ADV has difficulties measuring close to a solid surface, so 539 data out of the 588 reported by Bai et al. (2014) were used here.

Calibration methodology

A given physical phenomenon, in a domain Ω, can be represented by a group of input parameters (Π) and a group of output (answer) parameters (Y) in a relationship of the form Π → ΩY. In the DOE procedure, those independent parameters are known as factors, the values they can take in the experiment are known as levels and the dependent parameters are known as answer variables (Gómez-Zambrano 2017). The purpose of a DOE is to find a proper combination of the Π parameters, the levels they can take in the experiment and the relevance that any of these Π parameters can take over the Y parameters (Montgomery 2012). For a numerical model calibration, a similar approach can be followed where the equivalent Π parameters may be related to the characteristics of the computational grid (grid size, elements form), the implementation of the boundary conditions, and the size of the computational domain compatible with the boundary conditions, among others. The Y parameters, in this case, are those results of the optimal combination of the Π parameters that make a minimum difference between the numerical results and a set of reference answer variables transformed into statistical indices data (Montgomery 2012).

Any methodology used to compare observed and simulated values must include at least: (i) a graphic procedure that illustrates the comparison between the calculated values and the observed; (ii) a dimensionless index (or relative error indicator) quantifier of the goodness-of-fit like the NSE; and (iii) an error indicator like the RMSE (Legates & McCabe Jr 1999; Ritter & Muñoz-Carpena 2013).

A flow chart of the suggested methodology presented in this work is shown in Figure 2. A critical point of this methodology is whether the goodness-of-fit indexes between the numerical model results and the reference answer variables are accepted or rejected. To resolve this critical point, we chose to follow the recommendations given by Legates & McCabe Jr (1999), and to include in the methodology the goodness-of-fit quantifier index (NSE) and the (RMSE) error indicator.

Figure 2

Flow chart of the suggested calibration methodology using DOE.

Figure 2

Flow chart of the suggested calibration methodology using DOE.

Within the DOE, the central statistical technique is the analysis of variance (ANOVA). Therefore, in this methodology, the ANOVA is used for selecting the significant factors (these are the independent variables) together with the standardized Pareto plot (this chart shows the absolute values of the standardized effects from the largest effect to the smallest effect; it also plots a reference line to indicate which effects are statistically significant to 95%) and half normal probability plot (this plot is a graphical tool that uses these ordered estimated effects to help assess which factors are important and which are unimportant). The ANOVA technique consists of separating the total variation in the variability due to the treatments (a treatment is a factor at a specified level; in our case, one treatment is equal to one simulation of the mathematical model) and the variability due to the error. When the former dominates the latter, it is concluded that the treatments have an effect, or in other words, the population means of the treatments are different. When the treatments do not dominate they contribute equal to or less than the error; it is concluded that the means are equal. In this work, the ANOVA applied to the design of factorial experiments allows determination of whether the simulations show significant differences with respect to the observed data. It also provides a safe criterion for the hydraulic engineer to make decisions about which of the factors has the greater effect on the response variable; in our case, the response variable is the RMSE. Additionally, it allows finding the optimal combination of the input data for the configuration of the numerical model so that the results produce a smaller difference with respect to the observed data. The meanings of these statistical concepts are defined with major profundity in the classical literature of the DOE, which is recommended to be studied for greater understanding of this work (Montgomery 2012).

Definition of the input parameters and the answer variables

The input parameters selected for this application are related to the element size and distribution of the computational grid. The computational domain, formed by a curved channel with two straight branches (upstream and downstream of the curve), is constructed by hexadric elements, in longitudinal and transverse directions. For this grid, we defined nine input factors related to the grid configuration and one input factor related to the velocity at the inlet of the channel as boundary condition, for a total of ten input parameters. This set of parameters is represented by Equation (1): 
formula
(1)

Equation (1) includes two groups of factors, namely, ND factors that are related to the number of grid divisions in every edge (Type and Number of Divisions in meshing of Ansys Inc.) and Se factors that are related to parameters for the grid cells' distribution (Bias Type and Bias Factor in meshing of Ansys Inc.). Each group of factors has a sub-index representing the channel width (W), the water depth in the channel curve (H), the straight channel upstream of the curve (U), the straight channel downstream of the curve (D) and the factors of the curve with sub-index (C). The skewness or distribution factor (Se) represents the ratio Lmax/Lmin, where Lmax is the maximum element length and Lmin is the minimum element length in a given direction. A unit value for SeH and SeW indicates a uniform grid of identical cell size in height and width, respectively, as shown in Figure 3(a), and a value of SeH larger than one indicates that the cells near the solid walls are finer than those away from it, as is observed in Figure 3(b). The NPI factor represents the number of vertical coordinates used to prescribe the logarithmic velocity profile at the channel entrance (Schlichting & Gersten 2000). The answer variable Y is any of the statistical indexes computed for the three velocity components (Vt, Vr and Vz) and for angle

Figure 3

Grid structure for a transverse section of the channel: (a) uniform grid for SeW = 1, SeH = 1; and (b) non-uniform grid SeH = 3 and uniform grid for SeW = 1.

Figure 3

Grid structure for a transverse section of the channel: (a) uniform grid for SeW = 1, SeH = 1; and (b) non-uniform grid SeH = 3 and uniform grid for SeW = 1.

Statistical indexes used to evaluate good fitting for the velocity field

The good-fitting index for the velocity direction (angle) perpendicular to the main flow is computed according to Equation (2): 
formula
(2)
where the subscript P refers to the predicted values and the subscript O refers to the observations, θO_mín and θO_máx are the maximum and minimum observed values of the angle and N is the number of observations.
According to Jorba-Casellas (2005), it is possible to evaluate good fitting for the 2D velocity field on a given plane by using a modified RMSE index given by Equation (3), where the velocities are on the r-z plane: 
formula
(3)
where Vmáx and Vmín are the maximum and minimum values of the velocity, respectively, chosen from the max Vz and the max Vr. For each velocity component, the proper RMSE indexes are given by Equations (4)–(6): 
formula
(4)
 
formula
(5)
 
formula
(6)
One of the disadvantages of the ANOVA analysis is that the methodology is not applicable for simultaneous comparison of multiple answer variables, therefore, to overcome this difficulty, a unique parameter (Xrms) is obtained by grouping the individual indexes as shown by Equation (7): 
formula
(7)
where xi is each one of the statistic parameters that are grouped (Equations (3)–(6)) and Nv is the number of these parameters.
According to Legates & McCabe Jr (1999) and Ritter & Muñoz-Carpena (2013), the NSE parameter (Equation (8)) can be used to complement the good-fitting analysis: 
formula
(8)
where SD is the standard deviation of the observations. The NSE is a non-dimensional parameter for the good-fitting analysis widely used in hydrological modelling and flexible enough to be used in hydrodynamic simulation, as reported by Ritter & Muñoz-Carpena (2013). NSE equal to one indicates perfect fitting.

Definition of secondary flow variables

In order to calibrate and validate the numerical model in a more general context, in addition to the basic variables (velocity field and water depths), two secondary variables were computed from the basic variables, looking for the model to reproduce the general structure of the flow. Some of the secondary variables used in this work are presented below.

Flow helicity

The concept of helicity was first suggested by Berger & Field (1984) for magnetic fields and some years later this concept was used in hydrodynamics to study swirling flows by Moffatt & Tsinober (1992). For secondary currents in a curved channel, the concept of helicity was introduced by Bai et al. (2014), who analysed the secondary currents in a 180° curved channel. The helicity is computed for cross sections (perpendicular to the channel longitudinal axis) of the curved channel as: 
formula
(9)
where He (L4/t2/L) is the absolute helicity density, U is the velocity vector along the cross section and ΔA is a differential area of the cross section.

Intensity of secondary currents

Several authors (Stranden 2007; Ghobadian & Mohammadi 2011) have suggested that the recirculating current can be computed in terms of the intensity of the secondary flow (I), which is computed as the ratio of the kinetic energy in radial direction to the kinetic energy in longitudinal direction, in the curved open channel, and is obtained from Equation (10): 
formula
(10)

Hydrodynamic model

The free surface flow in a curved channel is governed by the mass conservation and the Reynolds averaged Navier–Stokes (RANS) equations (Schlichting & Gersten 2000; Wells et al. 2010; Gholami et al. 2014; Ansys Inc. 2016), which for a Cartesian coordinate system for a steady state flow are shown by Equations (11) and (12): 
formula
(11)
 
formula
(12)
where ui are the time averaged velocity components in the xi Cartesian direction (for i = 1, 2 ,3), μ is the fluid molecular viscosity, ρ is the fluid density, ui′ are the turbulent velocity fluctuations, is the Reynolds stress tensor, p is the pressure and δij is the Kronecker delta tensor. The RNG κ-ɛ model was used for the turbulence closure, as it has been successfully used for simulating the flow in curved channel free surface flows (Matthews et al. 1998; Bai et al. 2014; Ansys Inc. 2016). This RNG κ-ɛ model better captures the side walls' effects (like flow separation) than the standard κ-ɛ model (Kleinstreuer 2003). This set of equations is implemented in the Ansys Inc. software used for the numerical experiments. The RNG κ-ɛ model has a similar form to the standard κ-ɛ model as in Equations (13) and (14): 
formula
(13)
 
formula
(14)
In these equations, κ is the turbulent kinetic energy per unit mass associated with eddies in turbulent flow, ɛ is the rate of dissipation of turbulent kinetic energy, is the effective viscosity, Gκ represents the generation of turbulence kinetic energy due to the mean velocity gradients, Gb is the generation of turbulence kinetic energy due to buoyancy, YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, and Rɛ is an additional term as the main difference between the RNG and standard κ-ɛ models. The quantities σκ and σɛ are the effective Prandtl numbers for κ and ɛ, respectively. Sκ and Sɛ are user-defined source terms. The model constants Clɛ and C2ɛ in Equation (14) have values derived analytically by the RNG theory; values used by a default option in ANSYS Fluent, as Clɛ = 1.42 and C2ɛ = 1.68. For the quantities Rɛ, Gκ, Gb and C3ɛ ANSYS Fluent provides equations for their estimation (Ansys Inc. 2016).

Hydrodynamic simulation

The numerical solutions of Equations (11) and (12) are obtained with the Ansys Inc. software using four types of boundaries: (a) velocity inlet as inlet boundary at entrance, (b) pressure-outlet as open boundary at exit, (c) free surface on the top as symmetry and (d) close boundaries (solid walls). For the open boundary at the entrance, the logarithmic velocity profile (Schlichting & Gersten 2000) was prescribed given the 2 m straight branch of the channel downstream of the curve. For the open boundary at the exit, the zero pressure (atmospheric pressure) condition was used. For the solid walls (channel bottom and lateral walls) the non-slip condition with standard wall functions was prescribed and the free surface was modelled by the ‘rigid lid’ approximation, which is valid for free surface changes smaller than 10% of the still water depth upstream of the curve entrance, according to Bai et al. (2014). The numerical schemes used for the approximate solution of the governing equations, include: (a) the PISO algorithm for coupling velocity and pressure, (b) second order UPWIND scheme for the convective terms in the momentum equations and (c) the PRESTO algorithm for the pressure interpolation.

RESULTS AND DISCUSSION

Factors' detailed selection and optimal model configuration

Following the suggested methodology, we began by selecting the factors that show significant effect on the answer variable. The numerical experiment is designed for a fractional factorial of 210-5 with a IV resolution, allowing evaluation of the most important effects on the answer variable without confused effects or double interactions between factors. The DOE obtained give 32 model runs. After computing the Xrms as the answer variable for each numerical experiment, an ANOVA analysis was performed showing that, for a 95% confidence level, SeH and SeW, in their low level, were the factors showing significant effect on the answer variable (Xrms). This result shows that small variations at their level have a large impact on the answer variable, and produce a minimum value for a uniform size element distribution. This result indicates that a refinement of the elements near the channel walls is not needed as it does not improve the solution. Next in order of significance are the NDH, NDE and NDW factors at their low levels. The standardized Pareto plot (Figure 4(a)), for a 95% significance level (a standardized effect of nearly 2 in Figure 4(a), intermittent vertical line) shows that only the SeH and SeW factors show a significance level higher than 95%. As confirmation of this result, the half-normal probability plot, Figure 4(b) shows the SeH and SeW factors as the more relevant factors by being further away from the origin of the plot.

Figure 4

Statistical analysis results: (a) standardized Pareto plot and (b) half normal probability plot. The symbols (+) and (−) in (a), point out the factors' effect, positive and negative, on the response variable, Xrms, respectively.

Figure 4

Statistical analysis results: (a) standardized Pareto plot and (b) half normal probability plot. The symbols (+) and (−) in (a), point out the factors' effect, positive and negative, on the response variable, Xrms, respectively.

Those factors that have no significant effect on the answer variable can be set to a convenient level in order to minimize the computational effort. In particular, the NPI factor needs to be specified to a small number of points along the vertical to adequately describe the logarithmic velocity profile at the channel entrance.

The model optimal configuration obtained through the ANOVA analysis provides the factor levels that minimized the answer variable (Xrms). Then, those factors related to the computational grid were taken to a non-dimensional factor by scaling them by a characteristic length of the channel, which in our case was the hydraulics radius (RH) at the entrance section in the curve of the channel. Defining Te as the size of the mesh element, the results are: TeW = 0.1166 RH, TeH = 0.0486 RH, TeU = TeD = 0.2333 RH, and NDC = 2° (the subscripts are the same as described above). For those factors related to the element distribution in the domain, it was found that a skewness larger than 1 has no effect on Xrms, a result that indicates that elements are uniformly distributed in the whole domain. When the elements' aspect ratio (in the mesh metric of Ansys Inc.) is checked, a 4.8 value was obtained, which can be used later as a reference value for other cases (computational grids). For this optimal configuration, the model results will be compared to the reference data provided by the laboratory measurements.

Good fitting for the calibration stage

Qualitative comparison

Figure 5 shows the observed (Bai et al. 2014) and modelled structure of the swirling flow in the curved channel. The helicity (Figure 5(a)) at several cross sections along the curve is well reproduced by our numerical model (S1), following closely the experimental data and the simulations reported by Bai et al. (2014), although our results slightly over-predict the observations. Notice how the helicity obtained by the non-optimal grid clearly underestimates the observations, including some differences in the shape of the plot. Regarding the secondary current intensity (Figure 5(b)), our model and the simulations by Bai et al. (2014) are very close. Both models over-predict the experimental data although the shape of the curve is maintained. For both variables, the simulations obtained by our model with the optimal grid correctly predict the location of the peaks (cross section 60° downstream of the entrance), an important feature of this type of flow reported by the literature (Stranden 2007; Bai et al. 2014).

Figure 5

Spatial variation of (a) helicity and (b) intensity of the secondary flow: (O): observed data by Bai et al. (2014); (S1): simulation with the optimal grid; (S2): simulation by Bai et al. (2014); (S3): simulation with a non-optimal grid.

Figure 5

Spatial variation of (a) helicity and (b) intensity of the secondary flow: (O): observed data by Bai et al. (2014); (S1): simulation with the optimal grid; (S2): simulation by Bai et al. (2014); (S3): simulation with a non-optimal grid.

An additional qualitative analysis was performed by comparing Figure 6(a), the streamlines at the cross section 120° from the entrance to the curve, simulation with the optimal grid, simulation by Bai et al. (2014) and simulation with a non-optimal grid. These figures show good agreement among the observed data and S1, S2 simulation, showing a strong and large counter-clockwise vortex with its eye located in the middle of the cross section and a smaller anti-clockwise vortex near the free surface at the inner wall. This last vortex is clearly observable for the simulations (S1 and S2 cases) and not so clearly for the observed data. In contrast, one simulation (S3 case) shows only one large and strong anti-clockwise vortex with its eye located much closer to the inner wall.

Figure 6

Streamlines of the secondary current at the cross section located 120° from the entrance: (a) observed data by Bai et al. (2014); (b) simulation with the optimal grid; (c) simulation by Bai et al. (2014); (d) simulation by non-optimal grid. The velocity vectors have a geometric scale of 1:1.

Figure 6

Streamlines of the secondary current at the cross section located 120° from the entrance: (a) observed data by Bai et al. (2014); (b) simulation with the optimal grid; (c) simulation by Bai et al. (2014); (d) simulation by non-optimal grid. The velocity vectors have a geometric scale of 1:1.

The next challenge is to quantify the good-fitting behaviour of the simulations considering the characteristics of the vortex structure, characteristics like the location of the vortex eye and the vortex intensity.

Measure of the absolute error

Based on the obtained optimal factor configuration and the results from the ANOVA analysis, the model results are compared, in terms of the goodness-of-fit, to those reported in the specialized literature. The comparison is performed with the data reported by Bai et al. (2014) and Patil et al. (2014) and by using Equations (2), (4)–(6), and whose results are presented in Figure 7(a). It is clear that the magnitude of the RMSE error for the three velocity components and for the angle θ for all cases and for all reference data are of the same order of magnitude. Patil et al. (2014) did not report error values for Vt because they considered a 2D simulation. Note that the RMSE values calculated with the poor mesh simulation (S3) and the observed data (Figure 7(a)) are very large compared to the other options; this indicates that the proposed methodology is able to reduce errors between predicted and observed values in an efficient manner. In Figure 7(b), we compare the RMSE2D and Xrms, obtained with our simulation in optimal mesh conditions (S1) and those obtained with the simulation performed by Bai et al. (2014); note in this figure the minimum difference between them, which indicates that the configuration of the numerical model resulting from applying to the proposed methodology points in the right direction.

Figure 7

Statistical indexes obtained with our model against cases reported in the literature. (a) RMSE index and (b) Xrms and RMSE2D indexes.

Figure 7

Statistical indexes obtained with our model against cases reported in the literature. (a) RMSE index and (b) Xrms and RMSE2D indexes.

Measure of the relative error

As a measurement of the relative error, the NSE coefficient for the axial velocity component is computed for the Bai et al. (2014) simulations (Figure 8(a)) and for our simulation (Figure 8(b)). The comparison was calculated using the FITEVAL software developed by Ritter & Muñoz-Carpena (2013). We compared 539 data of axial velocity to the observed data for the seven cross sections of the curved open channel. For both simulations, the NSE coefficient is of the same order of magnitude: 0.721 for Bai et al. (2014) and 0.765 for our simulation. The goodness-of-fit of the observed data to the simulated data by Bai et al. (2014), Figure 8(a), achieves 62.8% probability of obtaining values of the efficiency (NSE) coefficient greater than 65% (acceptable). For our simulated data, we achieve 81% probability of obtaining values of the efficiency (NSE) coefficient greater than 65% (Figure 8(b)) (acceptable). The NSE criterion is stricter than the RMSE index and none of the simulations (including those reported in the literature) reach a value of 95% corresponding to a qualification of ‘very good agreement’, this fact becoming a challenge in hydrodynamic simulations.

Figure 8

Goodness-of-fit for the axial velocity component by the NSE coefficient, generated with FITEVAL: (a) simulations by Bai et al. (2014) and (b) our simulation.

Figure 8

Goodness-of-fit for the axial velocity component by the NSE coefficient, generated with FITEVAL: (a) simulations by Bai et al. (2014) and (b) our simulation.

It is relevant to consider that in order to obtain a ‘very good agreement’ between the measurements and the results of a numerical model, it is necessary to maintain a coordinated procedure between the observed and the numerical simulations given the difficulties of applying exactly the same conditions for both procedures. It is not always possible to specify the same conditions for the numerical model, for example, irregularities in the geometry of the experimental channel presented after the channel construction (walls along the channel curve, junctions between the straight branches of the channel and the channel curve, etc.), irregularities that may induce some flow conditions (flow bursts) that cannot be captured by the numerical model because the irregularities are not introduced in the model geometry.

Meaning and use of the results

The literature does not report calibration of mathematical models to simulate vortical flow by using representative variables of this flow, variables such as helicity and intensity of secondary currents, showing that they have practical applicability despite being a theoretical concept. In open channel flows the meaning of these results is important for the hydraulic engineer because in the methodology presented here, it is possible to incorporate several factors that, together with DOE, allow study of secondary flow in curved channels.

The proposed methodology contributes to saving time in the calibration of a hydrodynamic numerical model since it follows a systematic procedure, which allows a significant reduction in the number of simulations compared, for example, with the classic OFAT methodology.

This methodology has applicability in projects of hydraulic engineering where it is necessary to calibrate hydrodynamic models, whose comparison variables correspond to the velocity field of a flow, to study the spatial distribution of velocities, vortices within the flow, transport of solid particles and sediments, and vortex hydrodynamic separators in particle removal.

CONCLUSIONS

A new methodology for the calibration of numerical models of flow in curved channels is presented and validated. The implemented numerical model is based on the mass conservation, the RANS equations and the RNG κ-ɛ model for turbulence closure. The numerical experiments were designed by following the factorial fractional DOE procedure to study the effects of the computational grid and other factors on the model results. The quality of the numerical results compared to the observed data was studied using goodness-of-fit criteria.

It was shown that the suggested methodology was successfully used to study the computational grid characteristics and the configuration of several factors related to the model boundary conditions. The results of this study showed that, after the proper factor selection, it is possible to define the optimal model configuration. The study showed also that the grid factors related to the longitudinal characteristics of the grid (NDE, NDS, NDC) have no significance regarding the answer variable. On the contrary, those factors related to the grid element distribution (SeE, SeS, SeW, SeH) have relevance on the model results in such a way that the optimal grid resulted in a uniformly distributed grid. This somehow unexpected result can be explained by the use of the RNG κ-ɛ turbulence closure model because near the solid walls a wall function law is used, avoiding the need for a finer grid in these regions.

The values of the statistical indexes RMSE, Xrms and RMSE2D, computed for the velocity field during the calibration stage, were of the same order of magnitude of those reported in the literature. This result was possible because of the numerical model optimal grid configuration obtained by applying the suggested methodology. Results were shown for the model non-optimal configuration, and its corresponding statistical error indexes were far inferior to those corresponding to the model optimal grid configuration.

Given the complexity of the secondary flow in the channel curve due to its tridimensional behaviour, and given the simulation conditions (steady state) used in this work to calibrate the model, the flow structure (helicity and intensity of the secondary flow) was used instead of the simulated velocity field. In some sense, the comparisons with measurements have an average meaning. Non-instant values were considered.

The methodology presented in this article provides a new tool for the calibration of numerical models, which facilitates the selection of the significant factors with safe statistical criteria to find the optimal combination of these factors, so that the results of the numerical simulation are similar to the measurements.

Despite the reasonable results obtained with this model configuration, more challenges need to be tackled in further research. When the numerical simulation is focused on the flow hydrodynamic characteristics like vortices, their size, their intensity, their rotation frequency, the suggested methodology may be improved with more sophisticated measuring and calibration techniques.

REFERENCES

REFERENCES
Alexandris
S.
,
Stricevic
R.
&
Petkovic
S.
2008
Comparative analysis of reference evapotranspiration from the surface of rainfed grass in central Serbia, calculated by six empirical methods against the Penman-Monteith formula
.
Eur. Water
21
,
17
28
.
Ansys Inc.
2016
ANSYS Help Viewer Fluent User's Guide
.
Berger
M. A.
&
Field
G. B.
1984
The topological properties of magnetic helicity
.
J. Fluid Mech.
147
,
133
148
.
https://doi.org/10.1017/S0022112084002019
.
Craig
R. G. A.
,
Loadman
C.
,
Clement
B.
,
Rusello
P. J.
&
Siegel
E.
2011
Characterization and testing of a new bistatic profiling acoustic Doppler velocimeter: the Vectrino-II
.
Presented at the 2011 IEEE/OES 10th Current, Waves and Turbulence Measurements (CWTM) Conference
.
Monterey, CA
,
USA
, pp.
246
252
. .
Ghobadian
R.
&
Mohammadi
K.
2011
Simulation of subcritical flow pattern in 180° uniform and convergent open-channel bends using SSIIM 3-D model
.
Water Sci. Eng.
4
,
270
283
.
Gholami
A.
,
Akhtari
A. A.
,
Minatour
Y.
,
Bonakdari
H.
&
Javadi
A. A.
2014
Experimental and numerical study on velocity fields and water surface profile in a strongly-curved 90° open channel bend
.
Eng. Appl. Comput. Fluid Mech.
8
,
447
461
.
https://doi.org/10.1080/19942060.2014.11015528
.
Gómez-Zambrano
H. J.
2017
Investigación Experimental Y Numérica del Proceso de Separación Fluido-Sólidos Usando Conductos Curvos (Experimental and Numerical Research of the Fluid-Solids Separation in Curved Open Channels)
.
PhD thesis
,
Universidad Nacional de Colombia
,
Sede Medellín
,
Colombia
.
Han
S. S.
,
Ramamurthy
A. S.
&
Biron
P. M.
2011
Characteristics of flow around open channel 90° bends with vanes
.
J. Irrig. Drain. Eng.
137
,
668
676
.
Jara-Rojas
F.
,
Ortega-Farías
S.
,
Valdés-Gómez
H.
,
Poblete
C.
&
del Pozo
A.
2009
Model validation for estimating the leaf stomatal conductance in cabernet sauvignon grapevines
.
Chil. J. Agric. Res.
69
,
88
96
.
Jorba-Casellas
O
, .
2005
Simulación de los Campos de Viento de la Península Ibérica Y el área Geográfica de Catalunya con Alta Resolución Espacial Para Distintas Situaciones Meteorológicas Típicas (Simulation of the Wind Fields of the Iberian Peninsula and the Geographical Area of Catalonia with High Spatial Resolution for Different Typical Meteorological Situations)
.
Doctoral thesis
,
Universitat Politècnica de Catalunya
,
Barcelona
,
Spain
.
Kleinstreuer
C.
2003
Two-Phase Flow: Theory and Applications
.
Taylor & Francis
,
New York
,
USA
.
Legates
D. R.
&
McCabe
G. J.
Jr
1999
Evaluating the use of ‘goodness-of-fit’ measures in hydrologic and hydroclimatic model validation
.
Water Resour. Res.
35
,
233
241
.
https://doi.org/10.1029/1998WR900018
.
Lye
L. M.
2002
Design of experiments in civil engineering: are we still in the 1920's?
Presented at the Annual Conference of the Canadian Society for Civil Engineering
,
Montréal, Québec, Canada
.
Matthews
B. W.
,
Fletcher
C. A. J.
&
Partridge
A. C.
1998
Computational simulation of fluid and dilute particulate flows on spiral concentrators
.
Appl. Math. Model.
22
,
965
979
.
https://doi.org/10.1016/S0307-904X(98)10030-6
.
Moffatt
H. K.
&
Tsinober
A.
1992
Helicity in laminar and turbulent flow
.
Annu. Rev. Fluid Mech.
24
,
281
312
.
https://doi.org/10.1146/annurev.fl.24.010192.001433
.
Montgomery
D. C.
2012
Design and Analysis of Experiments
,
8th edn
.
Wiley
,
Hoboken, NJ
,
USA
.
Patil
V.
,
Finn
J.
,
He
X.
,
Ziazi
R.
,
Apte
S. V.
,
Liburdy
J. A.
&
Wood
B.
2014
Experimental versus computational methods in the study of flow in porous media
. In:
ASME 2014 4th Joint US-European Fluids Engineering Division Summer Meeting Collocated with the ASME 2014 12th International Conference on Nanochannels, Microchannels, and Minichannels
,
Chicago, IL, USA
.
https://doi.org/10.1115/FEDSM2014-21886
Ritter
A.
&
Muñoz-Carpena
R.
2013
Performance evaluation of hydrological models: statistical significance for reducing subjectivity in goodness-of-fit assessments
.
J. Hydrol.
480
,
33
45
.
https://doi.org/10.1016/j.jhydrol.2012.12.004
.
Schlichting
H.
&
Gersten
K.
2000
Boundary-Layer Theory
,
8th edn
.
Springer Science & Business Media
,
Berlin
,
Germany
.
Song
C. G.
,
Seo
I. W.
&
Kim
Y. D.
2012
Analysis of secondary current effect in the modeling of shallow flow in open channels
.
Adv. Water Resour.
41
,
29
48
.
https://doi.org/10.1016/j.advwatres.2012.02.003
.
Stranden
D. T.
2007
Numerical Studies of Flow in Curved Channels
.
PhD thesis
,
University of Bergen
,
Bergen
,
Norway
.
Toledo
T.
&
Koutsopoulos
H. N.
2004
Statistical validation of traffic simulation models
.
Transp. Res. Rec. J. Transp. Res. Board
1876
,
142
150
.
Van Balen
W
, .
2010
Curved Open-Channel Flows. A Numerical Study
.
PhD thesis
,
Technische Universiteit Delft
,
Delft
,
The Netherlands
.
Wahid
Z.
&
Nadir
N.
2013
Improvement of one factor at a time through design of experiments
.
World Appl. Sci. J.
21
,
56
61
.
https://doi.org/10.5829/idosi.wasj.2013.21.mae.99919
.
Wells
J.
,
Salem Said
A.-H.
&
Saad
R.
2010
Effects of turbulence modeling on RANS simulations of tip vortices
.
Presented at the 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Aerospace Sciences Meetings
,
Orlando, FL, USA
.
https://doi.org/10.2514/6.2010-1104
Willmott
C. J.
,
Robeson
S. M.
&
Matsuura
K.
2012
A refined index of model performance
.
Int. J. Climatol.
32
,
2088
2094
.
https://doi.org/10.1002/joc.2419
.