Estimation of Discretization Errors in Viscous Flow Simulations Based on the RANS Equations: Wall-Resolved Versus Wall-Functions

The use of computational simulations in engineering has become almost unavoidable. The complexity of the mathematical models, as , for example, the Navier–Stokes equation in fluid flows, requires the use of numerical solution techniques to determine the quantities of interest. As a consequence, numerical errors have to be estimated (solution verification [1–4]) to guarantee solutions with an acceptable numerical accuracy. Numerical errors are a consequence of the finite precision of computers (round-off error), iterative solution procedures (iterative error), transformation of the continuum equations into a system of algebraic equations (discretization error), effects of initial conditions for unsteady problems and sampling for stochastic problems (statistical errors) [1,2,5].

In this paper, we focus on discretization error estimation for the simulation of wall-bounded viscous flows at high Reynolds numbers, i.e., Reynolds numbers in the range of $107$ to $109$⁠, which are typically found in naval hydrodynamics applications. In such conditions, the Reynolds-averaged Navier–Stokes (RANS) equations are still the most common mathematical model used in computational fluid dynamics (CFD). Furthermore, we will restrict ourselves to statistically steady flows (time-averaging) of single phase, incompressible, Newtonian fluids. Nonetheless, the topics discussed in this study we will also apply to more complex near-wall flows.

One of the main quantities of interest of viscous flow simulations is the shear-stress at the wall $τw$ that can be determined with two different options: directly from its definition, or using wall-functions. The first choice is naturally preferable [6], but it leads to near-wall cells with extremely high aspect ratio [7] that can be numerically challenging at $ReL=109$⁠. Wall-functions are based on the assumption that total shear-stress, mean shear-stress plus Reynolds stress, is approximately constant in the near-wall region of a turbulent flow over a flat plate, which leads to the so-called log-law [8] that relates the mean velocity parallel to the wall with $τw$⁠. Therefore, it allows the use of near-wall cells with a height significantly higher than the previous option leading to a decrease of the numerical difficulties.

Several of the techniques proposed in the open literature to estimate the contribution of the discretization error to the numerical uncertainty [9–14] rely on grid refinement studies and power series expansions that use the typical cell size $hi$ as the independent variable. The need to define $hi$ leads to the requirement of using geometrically similar grids, which are not always easy to generate. However, the generation of these grid families is dependent on the strategy selected to determine $τw$⁠.

When $τw$ is determined from its definition, the near-wall cells should be systematically refined using the same refinement ratio of all the remaining cells. On the other hand, when $τw$ is determined from wall-functions, the height of the near-wall cells should remain fixed for all the grids. In this study, we illustrate the consequences of generating grids sets that do not satisfy these requirements. Grid families with fixed height of the near-wall cells are used with $τw$ determined from its definition, while grid families with a variable height of the near-wall cells are tested with $τw$ determined using wall-functions. Furthermore, two strategies are tested for the generation of grid families for wall-functions simulations:

1. Wall-functions grids generated from grids appropriate for wall-resolved simulations merging near-wall cells in the direction perpendicular to the wall;

2. Grids generated with near-wall expansion ratios that depend on the number of cells in the vertical direction that are tuned to obtain a fixed near-wall cell height.

Naturally, as illustrated in this study, the second option will not lead to geometrically similar grids.

The selected test case is the statistically steady, two-dimensional flow of a single-phase, incompressible, Newtonian fluid over a flat plate at Reynolds numbers (⁠$ReL=V∞L/ν$⁠) based on the plate length L, incoming undisturbed velocity $V∞$ and fluid kinematic viscosity $ν$ of $ReL=107$ and $ReL=109$⁠. Two RANS turbulence models are tested for wall-resolved simulations: the two-equation Shear-Stress Transport $k−ω$ model (SST) [15] and the two-equation $k−kL$ model (KSKL) [16]. The choice of these two models is justified by the wall boundary conditions of the turbulence quantities included in both models [7]. The SST model has $ω$ tending to infinity at a solid wall, whereas both dependent variables of the KSKL model are equal to zero at the wall. Only the SST model has been used with $τw$ determined from wall-functions, which are implemented as “automatic wall functions” [17].

There are several advantages in using this simple test case:

• It is possible to generate sets of exactly geometrical similar grids with a wide range of refinement levels that satisfy the requirements of the two alternative techniques to determine $τw$⁠;

• It is possible to reduce iterative errors to machine accuracy with all selected turbulence models at both Reynolds numbers, with or without wall functions;

• All post-processing required to determine the quantities of interest can be performed with numerical techniques with an order of accuracy larger or equal than the discretization techniques used in the flow solver;

• It is the only test case where the application of wall-functions does not require any extra assumptions, as neglecting pressure gradient effects and cross-flow mean velocity components in the near-wall region.

The advantages of using sets of geometrically similar grids are also discussed in Ref. [18].

The selected quantities of interest include integral, surface and local flow quantities and numerical uncertainties are estimated using an updated version [19] of the procedure presented in Ref. [12]. Nonetheless, any technique based on power series expansions would lead to equivalent results to those reported in this study.

The remainder of this paper is organized in the following way: for the sake of completeness, Sec. 2 presents the mathematical model and the details of the determination of the shear-stress at the wall; the test case and the grid sets are presented in Sec. 3; the flow solver and the numerical details are presented in Sec. 4; Sec. 5 presents and discusses the results and the main conclusions of this paper are summarized in Sec. 6.

where $ψ˜j$ are the instantaneous dependent variables, $xi$ are the coordinates of a Cartesian coordinate system, and t is the time.

$ρ$ is the fluid density. $Vi$ and P are the mean values of the Cartesian velocity components and relative pressure¹ , respectively. $vi$ are the fluctuating part (turbulence) of the Cartesian velocity components and the overbar designates time-averaging. The Reynolds-stress tensor $−ρvivj¯$ generated by the averaging procedure requires a turbulence model to close the problem.

where k is the turbulence kinetic energy and $δij$ is the Kronecker symbol. The contribution of k to the normal stresses is absorbed in the mean pressure gradient.

The SST model solves transport equations for k and $ω$ and determines $νt$ from the ratio between k and $ω$ and a limiter derived from the Bradshaw hypothesis [20]. The KSKL model solves transport equations for k and $Φ=kL$ and determines the eddy-viscosity from $Φ$ including also a similar limiter to that used in the SST model.

where $μ$ is the dynamic viscosity of the fluid, $Vx$ is the mean velocity component parallel to the wall and y is the direction perpendicular to the wall. The determination of $τw$ from Eq. (4) requires a value of $Vx$ sufficiently close to the wall to obtain an accurate solution [7]. In a grid refinement study, the height of the near-wall cell must decrease with the increase of grid refinement. A fixed height of the near-wall cell will lead to a discretization error that would remain for an “infinitely fine grid.”

where $U+=Vx/uτ$⁠, $y+=uτy/ν$⁠, $κ=0.41$ and $C≃5.2$⁠. $uτ$ is the friction velocity defined by $uτ=τw/ρ$⁠. The upper bound of the log-law layer is supposed to be approximately equal to 15% of the boundary-layer thickness, but in terms of $y+$ it is strongly dependent on the local Reynolds number $Rex=V∞x/ν$ where $x$ is the distance to the leading edge of the plate.

Equation (6) allows the determination of $τw$ (⁠$uτ$⁠) from $Vx$ at any grid location below the top bound of the log-layer region. In order to avoid mixing modeling and numerical errors, the location where Eq. (6) is applied to determine $τw$ should remain fixed in a grid refinement study.

The two-dimensional computational domain is a rectangle with the incoming flow $V∞$ and the plate of length L aligned with the horizontal direction $x$⁠. The leading edge of the plate is located at the origin of the $(x,y)$ Cartesian coordinate system. The length of the domain is 1.5 L with the inlet located at $−0.25L$ and the outlet at 1.25 L. The height of the domain is 0.25 L. The computational domain illustrated in Fig. 1 has been tested in several previous studies, as, for example in Refs. [21,22], and [23].

Two Reynolds numbers $ReL$ are tested: $ReL=107$ and $ReL=109$⁠.

Uniform velocity profiles with the horizontal component $Vx=V∞$ and the vertical component $Vy=0$ are specified at the inlet. Turbulence quantities are also imposed at the inlet and its values are derived from a turbulence intensity $I=0.01V∞$ and an eddy-viscosity $νt=10−8ReLν$⁠. Pressure is extrapolated from the interior at the inlet.

Streamwise derivatives of all dependent variables are set equal to zero at the outlet boundary. Symmetry conditions are applied on the bottom boundary upstream (⁠$−0.25L<x<0$⁠) and downstream (⁠$L<x<1.25L$⁠) of the plate. Pressure is imposed at the top boundary and zero normal derivatives of all remaining dependent variables are set equal to zero.

where $β=0.075$⁠.

with the normal derivative of the mean horizontal velocity component obtained from Eq. (6).

where $Nplate$ is the number of faces on the plate surface. The values of $ri$ and $Nplate$ used in all the grids sets of this study are presented in Table 1. For the wall-resolved simulations, the grids with $ri=1$ include $1.75×106$ cells for $ReL=107$ and $3.98×106$ cells for $ReL=109$⁠.

For each Reynolds number, five sets (RD1, RD2, RD3, RD4, and RD5) of 13 geometrically similar grids were generated with exactly the same number of cells $Ncells$ and $Nfaces$ and a range of refinement ratio of $1≤ri≤8$⁠. For each doubling of $ri$⁠, there are five grids available. The coarsest grid of set RD1 for $ReL=107$ is illustrated in Fig. 2.

The difference between the five grid sets available for each $ReL$ is the size of the near-wall cells that changes according to the selected stretching parameter that is largest for set RD1 (largest height of near-wall cells) and smallest for set RD5 (smallest height of near-wall cells). Although this looks like the wrong choice (RD1 has the largest near-wall spacing), it guarantees that the near-wall cell size of the grid 1 of set RD1 is equal to those used in grid 2 of set RD2, grid 3 of set RD3, grid 4 of set RD4 and grid 5 of set RD5. Therefore, it is possible to combine the five grid sets to obtain five grids covering a grid refinement ratio of 2 that have the same height of the near-wall cells. The same selection can be made for grids 5 to 9 and for grids 9 to 13 to obtain three groups of five grids with a fixed near-wall cell height. We will designate this “combined set” by RDC. Naturally, the grids of this set are not geometrically similar as illustrated in Fig. 3.

An extra grid set (RD0) with significantly smaller near-wall cells than all the other sets was generated to obtain a reference solution for the SST model, due to its dependence on the near-wall cell size [7]. Two extra refinement levels were added to one of the RD sets to obtain the reference solutions with the KSKL model.

When $τw$ is determined from Eq. (6), the size of the near-wall cells is supposed to remain fixed for all the grids of a given set. An efficient way to achieve this goal and preserve geometrical similarity is to merge the near-wall cells of the RD grids to obtain the desired cell height. Two sets (WF1 and WF2) have been generated with this technique. The number of cells merged to generate the WF grids is presented in Table 2. The location where $τw$ is determined is in the buffer-layer for sets WF1 and in the log-law region for sets WF2. The left plot of Fig. 4 illustrates the near-wall region of two grids of the WF2 set.

In complex geometries, it is not easy to merge cells in the near-wall region to obtain a fixed height. Therefore, two alternative techniques (WFA and WFB) were tested: in the WFA grids the stretching parameter was selected in each grid to obtain the desired near-wall cell height (equal to the WF grids); the WFB grids are geometrically similar and the stretching parameter was tuned to obtain the cell height of the WF grids for $ri=2$⁠. Both sets are illustrated in the middle (WFA) and right (WFB) plots of Fig. 4. The WF grids satisfy the two desired properties, geometrical similarity and fixed near-wall cell height, whereas the WFA are not geometrically similar and the WFB grids do not have a fixed near-wall cell height. Furthermore, in the WFA and WFB grids the cell height increases always in the vertical direction, whereas the WF grids have much thinner cells than the near-wall cell above it.

All simulations were performed with ReFRESCO which is a flow solver based on a finite volume discretization of the continuity and momentum equations written in strong conservation form. The solver uses a fully-collocated arrangement and a face-based approach that enables the use of cells with an arbitrary number of faces. Picard linearization and mass conservation is ensured using a simple-like algorithm, [28], and a pressure-weighted interpolation technique to avoid spurious oscillations, [29]. In the present study, a segregated approach was adopted for all simulations, which means that momentum, pressure correction (mass conservation) and turbulence equations are solved sequentially for each nonlinear iteration. Nonetheless, Anderson acceleration [30] is used in most of the simulations performed in this study. Thorough code verification is performed for all releases of ReFRESCO, [31].

Second-order schemes were used for the discretization of the RANS and turbulence models transport equations, with the exception of the convective terms of the turbulence quantities transport equations that were discretized with first-order upwind. It should be mentioned that for the flow over a flat plate first and second-order upwind in the convective terms of the turbulence quantities lead to very similar numerical accuracies for the mean flow quantities as illustrated in Ref. [32]. The same trend is observed in flows around smooth streamlined geometries as, for example, submarines.

All simulations were performed with 14 digits precision and iterative convergence criteria required the maximum value of the residuals of the L_∞ norm of the momentum and continuity equations to be smaller than at least $10−9$⁠, which is almost machine precision for $ReL=109$⁠. Therefore, for all selected quantities of interest, numerical errors are dominated by the discretization error.

The updated version [19] of the procedure proposed in Ref. [12] is used to estimate the discretization error $Ed(ϕ)$ of a quantity of interest $ϕ$ and its contribution to the numerical uncertainty $U(ϕ)$⁠. As for the Grid Convergence Index [9], $U(ϕ)$ is obtained from the multiplication of $|Ed(ϕ)|$ by a safety factor $Fs$⁠, $U(ϕ)=Fs|Ed(ϕ)|$⁠. In Ref. [19], the value of $Fs$ is determined from the standard deviation of the least-squares fits performed to the data of the $ng$ grids that in this study is always equal to $ng=5$⁠.

The fit that leads to the smallest standard deviation selects the equation to determine $ϕo$⁠. However, Eq. (14) (two terms expansion) is only selected when the standard deviation of the fit is smaller than 5% of the mean data change (⁠$(ϕmax−ϕmin)/(ng−1)$⁠). This is the main difference between the Refs. [12] and [19] versions of the procedure.

The numerical uncertainty of the simulations performed with the two turbulence models (SST and KSKL) at the two Reynolds numbers, $ReL=107$ and $ReL=109$⁠, was estimated for an integral quantity, surface quantities and local flow quantities.

obtained from integration of $τw$ on the plate surface. The integral is determined with a midpoint second-order rule. The surface quantities are the skin friction coefficients $Cf$ at 31 equally-spaced locations (⁠$Δx=0.025L$⁠) between $0.1L<x<0.9L$⁠.

where the exponent $q=1.25$ (⁠$ReL=107$⁠) and $q=2$ (⁠$ReL=109$⁠) for the wall resolved simulations. For the simulations using wall functions $q=1$⁠. Naturally, for the wall-functions simulations, some of these locations are inside the near-wall cell and so they will not be used. The values of $Vx$ at the selected locations are determined with one-dimensional cubic interpolations performed sequentially in the horizontal (x) and vertical (y) directions, respectively.

For all quantities of interest and grid sets, three uncertainty estimates are performed: grid with $ri=1$ using five grids included in $1≤ri≤2$⁠; grid with $ri=2$ using five grids included in $2≤ri≤4$ and grid with $ri=4$ using five grids included in $4≤ri≤8$⁠.

The maximum value of the height of the near-wall cells in wall coordinates (⁠$ymax+$⁠) is presented in Fig. 5. The plots contain the results obtained with the SST model. Similar results were obtained with the KSKL model, but the RD0 set was only used with the SST model. Levels of $ymax+$ tested for $ReL=107$ were deliberately chosen larger than what is required by this turbulence model [7], but this choice allows to illustrate also the consequences of too large values of $ymax+$ as indicated by the gray line (⁠$ymax+>1$⁠).

Figure 6 presents $ymax+$ for the grid sets used with wall-functions. Two levels of $ymax+$ are tested for each Reynolds number: for sets WF1 $ymax+$ is in the buffer-layer, whereas sets WF2 have $ymax+$ in the log-law region. We recall that grids WFB are geometrically similar and have a near-wall cell size that changes with grid refinement. As mentioned above, the stretching parameter of this set was tuned to match the near-wall cell height of the other two sets (WF and WFA) for $ri=2$⁠. The most challenging cases are the sets WFB1 that include grids in the buffer-layer and in the log-law region.

In order to compare numerical error estimations for the wall-resolved simulations, it is useful to determine a reference solution. For the SST simulations, the reference solution is obtained from a sixth set RD0 using significantly smaller heights of the near-wall cells, as illustrated in Fig. 5. For both turbulence model and Reynolds numbers, three extra levels of grid refinement were included in the RD sets used to estimate the reference solution. Some of these extrapolations to obtain the reference solution are presented in the figures of the following sections.

The convergence of the friction resistance/drag coefficient $CF$ with grid refinement is presented in Fig. 7. As expected [22], the results are clearly dependent on the selected turbulence model. For the $k−ω$ SST there is a strong dependence of the error level and of the observed order of grid convergence p on the size of the near-wall cells. For the coarsest grids (⁠$ri≥4$⁠) of $ReL=107$⁠, $ymax>1$ and so it is impossible to make a reliable error estimation for the sets RD1 to RD4. For the RDC sets (fixed near-wall cell size) of both Reynolds numbers, the exact solution is estimated with Eq. (11) for the three levels of $ri$ tested, which suggests a reliable error estimate. However, the numerical error is significantly under-estimated for the six cases tested. These results illustrate the risks of using the wrong grid generation strategy in the estimation of discretization errors with grid refinement studies.

On the other hand, for the KSKL model the influence of the near-wall cell size on $CF$ is much weaker than that observed for the SST model. Therefore, the error estimates performed for the RDC sets data are similar to those obtained from the data of the correctly designed grid sets. Nonetheless, for $ri=4$ and $ReL=107$ the inconsistency of the error estimation performed with the data of the RDC is also visible.

The convergence of $CF$ with grid refinement of the simulations performed with the $k−ω$ SST using wall-functions is illustrated in Fig. 8. The results obtained in sets WFB1 show that there is a significant influence on $CF$ of the near-wall cell height in the buffer-layer (⁠$y+<50$⁠). Furthermore, there is also a slight change of $Cf$ in the log-law region for $ReL=107$⁠. Therefore, error estimations performed with grid sets using a variable near-wall cell height can be unreliable.

The results obtained for sets WF and WFA exhibit much less grid dependency than that observed for $CF$ in the wall-resolved simulations. However, the values obtained in the two sets tested are not identical and the WFA grids do not converge to the same result of the WF sets, which are the sets that satisfy all the requirements for the use of wall functions. None of the wall-functions simulations matches the estimated reference solution of the SST model in the wall-resolved simulations. Nonetheless, the differences between the $CF$ values obtained in the finest grids of the wall-resolved and wall-functions simulations are significantly smaller at $ReL=109$ than at $ReL=107$⁠.

For the 31 selected locations, we have determined $E(Cf)ri=(Cf)ri−(Cf)o$ using the reference solution $(Cf)o$ to calculate the ratio between the estimated numerical uncertainty $U(Cf)ri$ and $E(Cf)ri$ that we will designate by $RUE(Cf)ri$⁠. Figure 9 presents the percentage of locations where the uncertainty estimates are unconservative, , i.e., $U(Cf)ri<E(Cf)ri$⁠. Furthermore, four levels of $RUE(Cf)ri$ are tested: $RUE(Cf)ri<1$⁠, $RUE(Cf)ri<0.5$⁠, $RUE(Cf)ri<0.1$⁠, and $RUE(Cf)ri<0.05$⁠. The percentage of locations where each of the error estimation equations (11)–(14) is used is presented in Fig. 10.

It is clear that for the $k−ω$ SST model the RDC grid sets (fixed near-wall cell size) lead to unconservative numerical uncertainty estimates. However, most of these wrong estimates are based on an acceptable order of grid convergence p, as illustrated by the results of Fig. 10. On the other hand, for the RDi sets the error estimation deteriorates when the near-wall cell size corresponds to $ymax+>1$⁠, which is especially visible for the simulations performed for $ReL=107$⁠.

Although the results of the KSKL turbulence model do not show the dependence on the near-wall cell size of the $k−ω$ SST model, there is a significant percentage of locations where the RDC grid sets lead to unconservative numerical uncertainties, especially for the simulations performed at $ReL=107$⁠. In the finest grids of the RDi sets, it is troublesome to determine the observed order of grid convergence p, because the convergence tends to be nonmonotonic. However, for $1≤ri≤4$⁠, the $Cf$ values obtained in the 5 available grids exhibit changes only in the fourth significant digit.

As an example of the grid convergence properties obtained for $Cf$⁠, Fig. 11 presents $Cf$ as a function of $ri$ for the wall-resolved simulations. As for $CF$⁠, the plots also present the least-squares fits and the observed order of grid convergence.

The skin friction coefficient $Cf$ in the region of the 31 selected locations is illustrated in Fig. 12 for the three levels of grid refinement tested. The plots also present the wall-resolved solution obtained in the finest grid of the RD0 set for comparison purposes. The results obtained in the WFB1 and WFB2 grids sets are not included in the plots of Fig. 12, because they lead to very large numerical uncertainties and significant changes with grid refinement. As discussed above, these results are a consequence of near-wall cells of different sizes including several grids with the near-wall cell size in the buffer-layer (⁠$y+<50$⁠).

All the simulations performed with wall-functions in grid sets WF and WFA exhibit a very weak dependence on the grid refinement level that was not observed for the wall-resolved simulations. At $ReL=107$ there is a larger influence of the near-wall cell size than at $ReL=109$⁠. Furthermore, for $ReL=109$⁠, with the exception of the WFA2 grid, all results at $ReL=109$ are graphically similar and in agreement with the wall-resolved solution. We note that the plots of Fig. 12 do not include the leading edge region, where the wall-resolved and wall-functions are inevitably different. We also recall that in the grids of set WFA the cell height increases with the distance to the wall and so the number of cells in the boundary-layer region decreases significantly between the WF sets and the WFA sets.

Figure 13 presents the grid convergence of $Cf$ at $x=0.5L$ for the wall-functions simulations. The trends observed in the data are similar to those discussed above for $CF$⁠.

For each Reynolds number and grid set, numerical uncertainties $U(Vx)ri$ and error estimates $E(Vx)ri$ were performed at 1550 (⁠$31×50$⁠) locations in the boundary-layer region. The percentage of locations that leads to unconservative numerical uncertainty estimates is presented in Fig. 14. As for $Cf$⁠, four levels of $RUE(Vx)ri$ are presented in the plots. Figure 15 presents the percentage of locations where each of the error estimation equations (11)–(14) is used.

The results confirm the trends discussed above for $Cf$ with the $k−ω$ SST data obtained in the RDC sets leading again to severely under-estimated numerical uncertainties. For the $k−ω$ SST, it is also visible the deterioration of the uncertainty estimates in the RDi grids that exhibit $ymax+>1$⁠. This tendency is confirmed by the percentage of locations that lead to an acceptable order of grid convergence, as illustrated in Fig. 15.

The results obtained with the KSKL turbulence model show that the numerical uncertainty estimation in the RDi sets is more troublesome for the finest grids than for the coarsest grids. However, as for $Cf$⁠, the differences between the $Vx$ results obtained in the finest grids are significantly smaller than those obtained in the coarsest grids. Furthermore, for $ReL=107$ the observed order of grid refinement exhibits unacceptable levels or is impossible to determine in many of the selected locations.

The trends observed in the grid convergence properties of $CF$ and $Cf$ are also observed in Fig. 16 that presents the grid convergence of the mean horizontal velocity component $Vx$ at an y location with $x=0.5L$⁠.

We have selected the velocity profile at $x=0.5L$ to illustrate the results obtained using wall-functions. $yo+$ is calculated from the friction velocity obtained from the reference solution calculated with the wall-resolved simulations, $yo+=(y/L)ReL(Cf)o/2$⁠. Figure 17 depicts the $Vx$ profiles at $x=0.5L$ obtained in the wall-functions simulations. As for $Cf$⁠, the results obtained in the WFB grid sets are not presented. The solution obtained in the finest grid of set RD0 of the wall-resolved simulations is also presented in the plots for comparison purposes.

The estimated numerical uncertainties of $Vx$ at the selected locations are significantly smaller in all the wall-functions simulations (WF and WFA grid sets lead to error bars smaller than the size of the symbols) than in the wall-resolved simulations. The results obtained for $ri=1$⁠, $ri=2$⁠, and $ri=4$ are graphically identical.

With the exception of the WFA1 results in the buffer-layer, there is graphical agreement between all the wall-functions simulations and the wall-resolved simulation. However, a quantitative assessment of the differences between the solutions at the selected locations shows two main trends: the WFA2 solutions are slightly above the data obtained in the WF1, WF2 and WFA1 sets; the agreement between the wall-resolved and wall-functions solutions is better for $ReL=109$ than for $ReL=107$⁠.

Figure 18 illustrates the grid convergence of the mean horizontal velocity component $Vx$ at an y location with $x=0.5L$⁠. As for the wall-resolved simulations, similar trends are observed in the grid convergence properties of integral, surface and local flow quantities.

This paper focuses on the simulation of high Reynolds numbers (⁠$ReL=107$⁠) viscous flows using the RANS equations with two eddy-viscosity turbulence models: $k−ω$ SST and $k−kL$⁠. The main goal of this study is to discuss the grid generation strategy required by the estimation of the contribution of the discretization error to the numerical uncertainty of integral, surface and local flow quantities using grid refinement studies and power series expansions. Two different options are explored:

• wall-resolved simulations that calculate the shear-stress at the wall from its definition;

• wall-functions simulations (only for the $k−ω$ SST model) that determine the shear-stress at the wall from a prescribed mean velocity profile that blends the linear sub-layer with the log-law region.

The selected test case is the statistically steady, two-dimensional flow of a single-phase, Newtonian, incompressible fluid over a flat plate at Reynolds numbers of $ReL=107$ and $ReL=109$⁠. This choice is justified by three main reasons: (1) the ability to reduce round-off and iterative errors to negligible levels when compared to the discretization error; (2) the generation of sets of geometrically similar grids that satisfy the requirements of wall-resolved and wall-functions simulations; (3) it is possible to determine the integral, surface and local quantities of interest using techniques with an accuracy larger or equal than that used in the flow solver.

For wall-resolved simulations, all grid cells should be systematically refined, which means that the height of the near-wall cells decreases with grid refinement. In this study, we have generated five sets of grids that allow combinations of five nonsimilar grids that have a fixed near-wall size, which is a technique often found in the open literature. On the other hand, wall-functions simulations require sets of grids with near-wall cells with a fixed height to avoid mixing of modeling and numerical errors. Three different grid generation strategies were tested: (1) sets with geometrically similar grids and a fixed near-wall cell height obtained from merging cells of the wall-resolved grids in the near-wall region; (2) sets with nonsimilar grids and fixed near-wall cell size obtained by tuning the stretching parameter used in the near-wall region; (3) sets with geometrical similar grids that reduce the near-wall cell height with grid refinement.

The results obtained for the friction resistance/drag coefficient, skin friction coefficient along the plate (31 locations) and mean horizontal velocity component (1550 locations in the boundary-layer) exhibited similar trends that suggest the following conclusions:

• Wall-resolved simulations

  1. The convergence properties of the selected quantities of interest are strongly dependent on the choice of turbulence model. The $ω$ boundary condition at the wall (⁠$ω→∞$⁠) makes the results obtained with the $k−ω$ SST model strongly dependent on the near-wall cell size. As a consequence, numerical uncertainty estimations performed for grid sets with a fixed size of the near-wall cells are not reliable. In such cases, the estimated numerical uncertainty significantly under-estimates the discretization error (more than one order of magnitude in many cases) for error estimates that exhibit an acceptable order of grid convergence. On the other hand, the grid convergence properties in the sets with geometrically similar grids deteriorate significantly for near-wall cell sizes with $y+>1$⁠, which makes the discretization error estimation troublesome (non-monotonic and noisy data).

  2. For the $k−kL$ model, the grid convergence properties of the selected quantities of interest exhibit a much weaker dependence on the near-wall cell size than the $k−ω$ SST model. Therefore, the use of the wrong grid generation strategy is mainly observed for the coarsest grids.

• Wall-functions simulations

  1. The results depend on the height of the near-wall cells, especially if the near-wall cells are in the buffer-layer (⁠$y+<50$⁠). Therefore, performing error estimates with grid sets that do not have a fixed near-wall cell size may lead to awkward results that are mainly a consequence of modeling errors.

  2. Numerical uncertainties estimated for the quantities of interest of the simulations performed in all the sets with a fixed near-wall cell size are significantly smaller than those obtained for the wall-resolved simulations.

  3. Differences between the solutions obtained in the two sets that have the same near-wall cell size grow with the increase of the near-wall cell height. This is not a surprising result because one of the sets reduces the number of cells in the boundary-layer with the increase of the near-wall cell height. Therefore, in grids with the smallest cell height for the near-wall cell it is best to adjust the near-wall cell height to fall within the range of $50<y+<100$⁠.

  4. Excellent agreement is obtained between the wall-functions and wall-resolved simulations for $ReL=109$⁠.

For wall-resolved simulations, most of these conclusions will also apply to RANS simulations of viscous flows around complex geometries in the same range of Reynolds numbers. Naturally, the trends obtained for wall-functions simulations are not general, because the modeling accuracy of wall-functions in three-dimensional geometries may be questionable, especially if there is flow separation.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.