Accuracy of boundary layer treatments at different Reynolds scales

Abstract Resistive forces associated to boundary layers (‘friction’) are usually out of scale in physical models of hydraulic structures, especially in the case of hydraulically smooth walls, generating distortions in the model results known as scale effects, that can be problematic in some relevant engineering problems. These scale effects can be quantified and corrected using suitable numerical models. In this paper the accuracy of using numerical simulation through the Reynolds Averaged Navier-Stokes (RANS) approximation in order to represent the head losses introduced by friction in hydraulically smooth walls is evaluated for a wide range of Reynolds scales. This is performed by comparing the numerical results for fully developed flow on circular pipes and between parallel plates against experimental results, using the most popular wall treatments. The associated numerical errors, mesh requirements and ranges of application are established for each treatment. It is shown that, when properly applied, RANS models are able to simulate the head losses produced by smooth wall friction accurately enough as to quantify the scale effects present in physical models. A methodology for upscaling physical model results to prototype scale, free of scale effects, is proposed.


Introduction
The design, optimization and verification of hydraulic structures have been traditionally founded on laboratoryscale physical model studies. Despite the advances in numerical modelling during the last few decades, in particular in CFD simulations, physical models are still one of the main tools supporting the engineering of complex structures [1].
A physical model would be completely similar to its real-world prototype counterpart if the following three similarity criteria were met: geometric, kinematic and dynamic similarity [2]. Geometric similarity requires that the shape of the problem is maintained, with a uniform length scale affecting all dimensions. For dynamic similarity, each force acting at any point over the continuum must scale consistently, with a single uniform and constant ratio, independently of the type of force. Kinematic similarity means that the trajectories of the flow must be preserved, with a uniform scale of time affecting all motions; this condition is naturally achieved when the previous two similarity restrictions are fulfilled [3]. Unfortunately, although geometric similarity can be easily attained, dynamic similarity is virtually impossible to achieve in the real world, as its enforcement would require performing the experiment in a miniature universe [4]. In particular, both the physical model and the prototype share the same (Earth) acceleration of gravity. As a consequence, in problems for which the free surface plays a significant role, then requiring the correct representation of the relationship between inertia and gravity forces (Froude scaling), all the remaining forces acting on the flow are by default out of scale, including viscosity, globally characterized by the Reynolds number. If these out of scale forces have any significance on the flow, then the physical model results no longer correctly represent the flow in the prototype, i.e., scale effects exist [3].
The physical model strategy to overcome this restriction stands on the fact that in most hydraulic engineering flows the Reynolds number is very high in the prototype, so viscosity effects are very small. Hence, to achieve a similar condition on the physical model, a low enough (approach-ing 1) length scale e L (defined as the ratio between length in the prototype and in the model) must be chosen. However, there are two difficulties with this approach. In the first place, physical model costs are approximately proportional to e L −3 [5], so economic limitations may lead to the selection of a smaller length scale than required to attain enough precision. The second problem arises from the assumption that a high Reynolds number (a global parameter) ensures that viscosity effects are negligible throughout the flow. However, this is not the case, for instance, near solid walls, where viscosity forces are dominant. In fact, when the wall roughness height is low enough, i.e., it is completely immersed in the viscous sublayer, the characteristics of this viscosity-dominated zone strongly affect the turbulent boundary layer, and hence eventually the main flow. In other words, hydraulically smooth walls are by default out of scale on physical models with Froude scaling. Scale effects have been broadly studied experimenting on physical models [6][7][8] usually by building a series of similar models at different scales. In this type of experiments, the largest scale model takes the place of the prototype, and scale effects are determined by comparing against it the performance of the remaining models. However, it is not possible to know if the largest model is itself affected by significant scale effects [4]. Studies comparing the behaviour between prototype scale structures and their physical model predictions have also been performed [9].
Many researchers have focused on determining the minimum scale needed to obtain a reasonably scale effects free representation for typical engineering flows, including, among others, flow on broad crest weirs [10], depths and aeration in hydraulic jumps [11,12], surface vortices on intakes [13], and trajectories of ski jumps [14]. However, even when using those scales, significant secondary phenomena might be triggered by out of scale forces [15].
The present work is a contribution to evaluate the current capabilities of Computational Fluid Dynamics (CFD) models to quantify scale effects in physical models, and to upscale the results to the prototype scale free from those effects. It stands on the fact that CFD models can represent, at least in principle, the combined effects of all forces acting both at the prototype and at the physical model scales. Hence, by simulating both the prototype and physical model flows, they should be able to elucidate the supposedly marginal effects of out of scale forces in the physical model. Various papers have reported using numerical models for this purpose [16][17][18][19][20][21][22], but little work has been done to evaluate the accuracy of this procedure. Even with limitations, CFD models should provide a more accurate extrapolation of results from the physical model to the prototype, as reported for the design of the Third Set of Locks of Panama Canal [19]. Now, CFD models are still subject to significant uncertainties, especially those associated to the modelling of turbulence. Various approaches exist to account for the effect of turbulence [23], including Direct Numerical Simulation (DNS), Large Scale Eddy Simulation (LES), and Reynolds-Averaged Navier-Stokes modelling (RANS). In recent years, a promising new line of research based on fractal calculus has been made available [24,25] that may provide a sound conceptual base for studying the properties of turbulence and their interaction with walls of different roughness; it might also be potentially useful to analyse hydraulic scale effects as a whole.
The RANS approach has been used to study a myriad of hydraulic and industrial problems [25,[27][28][29][30][31][32][33][34]. It is still the most commonly used modelling approach for engineering studies, as it combines reasonable computational costs with acceptable accuracy. In RANS models only the ensemble-averaged flow variables are solved explicitly, while turbulence effects are parameterized in terms of those flow variables. Variations in the parameterization of turbulence effects leads to different RANS models, each one claiming some advantages for particular problems. The models have been tested for a variety of flow conditions and problem scales over the years. However, there are little systematic records of RANS models validation for the same problem at different scales, so their accuracy on predicting scale effects is relatively uncertain.
In this paper, the accuracy of the most common RANS models when simulating head losses due to skin friction in hydraulically smooth walls, is evaluated for a wide range of scales, using the fully developed flow in circular pipes and between parallel plates as benchmark problems. The most widely used turbulence models and two alternative treatments to implement the boundary conditions at smooth walls, are considered, establishing their ranges of applicability and margins of error. Their accuracy is quantified by comparing the predicted head losses against experimental results. Some of the models employed in this paper have been validated for boundary layer problems by comparing velocity profiles or friction factors for some Reynolds numbers [35][36][37]. In this paper, the validation is systematically extended for the friction factor over a wide range of conditions. Only skin friction is considered in this paper; scale effects on drag due to flow separation are analysed in a forthcoming paper. Additionally, the methodology to predict prototype scale performance free of scale effects, through the use of numerical modelling, is described.
The plan of the paper is as follows. In the first place the benchmark problems are presented -this includes the formulas and experimental results -, and the errors that arise when simulating them on physical models with Froude scaling are quantified. In the second place, a brief yet comprehensive description of RANS models is presented, including the different wall treatments. In the third place, the errors in the numerical model are evaluated -this includes the formulation of best practices in order to minimize them -and compared to those incurred in physical models at different scales. Finally, the procedure to upscale the results to the prototype scale, is explained.

Description
Two problems are studied in this paper: flow along a pipe with a circular cross section, and flow between parallel flat plates (duct), both under fully developed flow conditions and for hydraulically smooth walls. The geometry for the first case is determined by the conduit diameter D, and for the second case by the distance between the plates H. The effect of the fluid viscosity is characterized by the Reynold number, which is defined as follows: where the subindex c and p indicate circular and flat plate geometries, respectively, U the mean velocity, and ν the kinematic viscosity of the fluid. This somewhat inconsistent definitions of the Reynolds number have been adopted due to their extensive usage in the literature. The definitions of the friction factor for pipes fc, known as Darcy's, and ducts fp, known as Fanning's, are as follows: The friction factors can also be expressed as a function of the pressure difference ∆p between both extremes of a pipe or duct of length L (taking into account the balance between friction and pressure forces): Experimental values for the friction factor published in the literature were employed as a basis for comparison to assess the accuracy of the model results.
For the pipe problem Blasius equation and Colebrook & White curve fit [38,39] where used: The experimental data from the Princeton Superpipe experiments [40], and the transitional data by Patel & Head [41] were also employed as a reference for this case.
In the case of the duct problem, a data fit plotted by Rothfus et al. [42] was taken for comparison. It is pertinent to mention that this plot was built by assuming that the nondimensional velocity profile is the same for the pipe and duct cases. Based on the dispersion of the data, an error of around 7% was estimated.

Errors in Froude scaled physical models
To put into context the precision of the numerical models, the magnitude of expected scale effects in a physical model with Froude scaling is first evaluated. Specifically, the pipe problem is analysed, which could be interpreted as representing a stretch of a hydraulically smooth pipe within a hydraulic system with a free surface, such as the one of a shiplock.
For Froude scaling, the Reynolds number of the physical model relates to the one in the prototype as: R model = R prototype e L −1.5 (5) Assuming that both the prototype and scale model walls are hydraulically smooth, and that the pipes have a circular cross section, the friction factor can be computed with the Colebrook & White equation. These friction factors are presented in Figure 1 as a function of the prototype Reynolds number and the length scale. It is observed that a substantial amplification of the friction factor occurs in the physical model as the scale decreases. The relation between the friction factor of the model and prototype is plotted in Figure 2. Note that in the turbulent flow regime the friction factor for the physical model can be from about 50% to 350% larger than the one for the prototype.
The magnitude of this error highlights the need for methodologies that allow its computation and correction. For this simple problem the experimental fit of Colebrook & White could be employed, but for more complex cases numerical models are the only alternative. This is the motivation to determine how accurately typical RANS models can solve head losses due to friction in smooth walls.

Basic conservation equations
Assuming incompressible monophasic flow, the mass conservation equation is: where U is the Reynolds-averaged velocity. Adopting the Boussinesq approximation to represent Reynolds stresses isotropically using an eddy viscosity ν t , the momentum conservation equation can be written as [43]: where ν is the kinematic viscosity of the fluid, and̃︀ p the Reynolds-averaged value of the modified kinematic pressure, which incorporates the contribution of k, the turbu-lent kinetic energy (TKE):︀ with p being the Reynolds-averaged value of pressure, and ρ the fluid density (assumed as a constant due to incompressibility).

Turbulence models
Several turbulence models belonging to the k-ε family were evaluated in this paper. k-ε models represent turbulence by keeping track of the evolution of the turbulent kinetic energy and its viscous dissipation rate, ε, through their conservation equations. The standard k-ε model [36] reads: where C 1 , C 2 , σ k , σε are constants, and P k is the production rate of TKE, which is expressed as: The two turbulence variables k and ε determine a length and a time scales for the dominant eddies, from which the expression for the eddy viscosity arises: where Cµ is a constant. The values of the constants for the standard k-ε model are: This k-ε model is the most widely validated and used RANS turbulence model, due to its efficiency and relative simplicity [44]. Its performance is considered as particularly satisfactory for confined flows in which Reynolds stresses are dominant, a frequent condition in engineering problems. Its performance is less satisfactory for non-confined flows, boundary layers with large curvature, and rotating flows [45].
The RNG k-ε model [46] was also tested. This is a variant of the standard model derived using Renormalization Group Theory, which states that all scales of motion contribute to the momentum diffusion instead of a single representative scale. The resulting equations are completely similar to those of the standard model, but with different values of the constants, and C 2 expressed in terms of the deformation tensor.
The realizable k-ε model [47] is a variant which introduces a variable Cµ coefficient, depending on both the deformation and rotation tensors. Now, in the boundary layer the velocity decreases to match that of the wall. Hence inertia gradually loses significance when approaching the wall, and viscous forces eventually become dominant. This means that turbulence models able to deal with low Reynolds number conditions are required to solve the boundary layer region explicitly. They are collectively known as Low Reynolds Number models, and are becoming increasingly popular, especially for problems in which the boundary layer is neither in equilibrium nor fully developed, as in fast hydraulic transients [48][49][50][51] and separated flows. The most popular k-ε Low Reynolds Number models were compared by Patel & Rodi [37], who concluded that the one with the best performance was Launder-Sharma's [35,52], followed by that of Chien [53] and Lam & Bremhorst [54]. The Launder-Sharma k-ε model equations are the following: These equations include the following weight functions: fµ = e −3.4/(1+Rt /50) in which R t is the local Reynolds number of the turbulence, defined as:

Numerical code
The conservation equations were solved using OpenFOAM [55,56], a set of applications and libraries to deal with continuum mechanics problems. The capabilities of OpenFOAM are similar to those of Flow-3D [18,57]. The differences lies in accessibility and numerics; while the latter one is a commercial code that uses the Finite Difference Method (FDM) over one or more structured meshes of hexahedra, OpenFOAM is an open-source free software based on the Finite Volume Method (FVM) with a collocated treatment for variables over non-structured meshes of polyhedral elements. In OpenFOAM the temporal discretization is implicit, with a pressure-velocity coupling analogous to that of Rhie & Chow [58].
As the numerical simulations presented in this paper are time independent, the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) [59] iterative procedure was used.

Wall treatments
Turbulent Boundary Layer theory identifies, adjacent to the wall, the Inner Boundary Layer [60], where the flow (contrary to the Outer Layer) is practically independent of the characteristics of the global flow, the shear stress is rather constant, and the Law of the Wall holds: where y is the distance to the wall, τw the wall shear stress (representative of the shear stress throughout the inner boundary layer), u * = √︀ τw/ρ the shear velocity, u + the dimensionless velocity, and y + the dimensionless wall distance. Three sub-regions can be distinguished, starting from the wall [45]: (i) the viscous sublayer, extending up to approximately y + ≈ 5 [60], where viscous forces are dominant, and the velocity profile is linear; (ii) the transition sublayer, which extends within the interval 5 < y + < 30, where both viscous and Reynolds stresses are significant; and (iii) the inertial sublayer, extending outwards y + > 30 and up to about y + ≈ 500 or less, where viscous stresses are already negligible, and the velocity profile is logarithmic, as observed experimentally and confirmed theoretically through dimensional analysis [60].
Based on this knowledge, different treatments have been developed to account for the boundary layer on CFD models. In this paper, the two most popular treatments are analyzed and compared, namely, Low Reynolds Number (LRN) and Logarithmic Wall Function (LWF).
Regardless of the type of treatment, the goal is to enforce the relation between the velocity parallel to the wall on the first computational node, U, and the wall shear stress. This can be accomplished by either introducing a momentum sink term on that node or locally manipulating the eddy viscosity, . OpenFOAM uses the latter approach. RANS equations near the wall assume that: in which Uw is the wall velocity. Considering static walls (Uw = 0), the following relationship is derived: which provides the value of the eddy viscosity to be assigned a the cells adjacent to the wall (i.e., its boundary condition) once a relationship between u + and y + is known.

LRN treatment
Under the LRN treatment the flow is explicitly solved throughout the entire boundary layer. Due to the relatively high velocity gradients within the boundary layer, a dense discretization is needed to attain accuracy. The first calculation node, at the center of the cells adjacent to the wall, has to lie inside the viscous sublayer (y + < 5). As viscous effects are significant both in the transitional and viscous sublayers, a Low Reynolds Number turbulence model is required. Now, for the viscous layer the velocity profile is given by [60]: Replacing (14) in (18), the boundary condition for the eddy viscosity arises: Experiments show that the TKE has a peak value in the transition sublayer, associated to the high velocity gradients. However, turbulent eddies rapidly decay inside the viscous sublayer, as shown by Klebanoff [61], so the appropriate boundary condition for k is also a zero value. Due to implementation details in OpenFOAM, null TKE values are not allowed, so a very small number must be assigned (for example, 10 −12 ).
As for the dissipation rate, ε, slightly different boundary conditions are used depending on the Low Reynolds Number model, consisting either of a zero normal gradient or a zero value [37]. In the case of the Launder-Sharma model, a zero value is recommended. As for the TKE, Open-FOAM does not allow a zero value for ε, so a very small one must be used instead.
This Low Reynolds Number treatment assumes from the start that the wall is hydraulically smooth.

LWF treatment
The LWF treatment consists of extending the flow simulation down to a certain point contained within the inertial sublayer, where a wall function is imposed.
On the inertial sublayer, the velocity profile can be expressed as [60]: in which E is a constant (= 9.8) for large enough Reynolds numbers and smooth wall. The first computational node must lie within the inertial sublayer, i.e., it has to verify the condition . Compared to the LRN treatment, using wall functions greatly reduces the computational cost due to the very significant reduction in the number of computational nodes. This explains its popularity.
Experiments indicate that in the inertial sublayer TKE production and dissipation rate are in equilibrium [45]. Therefore, from equations (9) and (10), and taking into account the logarithmic profile for U, the following equations can be derived: Equation (17) is used to impose the boundary conditions for ε. In the case of k, instead of (17) a zero gradient condition, ∂k/∂y = 0 (based on the equality of production and dissipation) is usually adopted due to numerical stability reasons. Now, equation (21) must be solved iteratively, as τw is not known a priori. An additional complication is that τw may become locally null at points of flow detachment, so this special case has to be dealt with specifically in the code in order to prevent from spurious results. In order to avoid these difficulties, Launder & Spalding [36] proposed an alternative computational procedure, consisting in deriving the shear velocity from the TKE using equation (17):

Implementation of the numerical tests
To test the accuracy of the different wall treatments and turbulence models, RANS models were implemented for each condition on the two flow cases: flow along a pipe with a circular cross section, and flow between parallel flat plates (duct). Each model was run to convergence for different Reynolds numbers, and the resulting friction factors were computed and compared to experimental values. Both the diameter of the conduit D and the distance between the parallel plates H were fixed at 0.10 m. A mean velocity of 1.00 m/s was imposed at the upstream boundaries. Periodic boundary conditions were applied for the velocity distribution and the turbulent properties, in order to achieve fully developed flow conditions independently of the model length. The fluid viscosity was varied according to equation (1) in order to get different values for the Reynolds number the same model geometry and mean velocity.
Due to the radial and lateral symmetries of the pipe and duct problems, respectively, both of them are amenable to a two-dimensional approach. Rectangular elements were selected for the computational grids. The size of the first layer of elements adjacent to the wall was adjusted, for every Reynolds number, in order to achieve a correct y + value at the first node (y + 1 ), as required by the adopted wall treatment. As the LRN model is quite sensitive to the spatial discretization close to the wall, the cell size was kept constant up to y + ≈ 20, and then increased as a geometric progression, with an expansion factor R = 1.1. The same geometric progression was used for the case of the wall function treatment, starting from the cell closest to the wall.
The experimental values for the friction factor published in the literature were used to estimate a priori the required distance of the first node to the wall (hence avoiding an iterative procedure).
The testing methodology is summarized in Figure 3. Further details of the model implementation are provided in Table 1

Results for the LRN treatment 4.2.1 Model performance
Tests were undertaken using the LRN treatment, described in equations (19) and (20), with the Launder-Sharma k-ε turbulence model for a wide range of Reynolds numbers, from 2 × 10 2 to 1 × 10 7 . Different mesh configurations were tested until mesh convergence was achieved, with a tolerance of 1% on the friction factor. The thickness of the first set of layers for the converged mesh was ∆y + = 0.25 for 0 ≤ y + ≤ 20 (hence, y + 1 = ∆y + /2 = 0.125). The friction factor was computed from the model results using equation (3). Figure 4 shows the comparison with the experimental curves. An overall satisfactory agreement is observed for the two problems, and for both the laminar and the turbulent regimes. However, it is observed that the model is not able to represent the sudden changes of the friction factor within the transition range, but instead it smoothly links the curves associated to laminar and turbulent flow. This is not considered as a practical limitation, as the usual engineering studies do not deal with transitional flows. Figure 5 shows the relation between the computed and experimental friction factors. For the circular pipes the relative error is below 1% in the laminar regime, increases up to 9% in the turbulent regime, then decreasing with the Reynolds number. For the parallel plates the trend is qualitatively similar, but a higher maximum error is observed for the laminar regime, of about 2%, and lower one for turbulent flow, of about 6%. However, it should be kept in mind that the reference friction factor curve for the parallel plates flow is relatively uncertain (7% error), so the deviation trends quantified for pipes should be considered as more indicative of the actual model performance.
When these errors are compared with those associated to a Froude scaled physical model, as presented in Figure 2, it is observed that the latter ones are practically negligible, so the numerical model results can be considered as very accurate from the engineering point of view.

Underprediction for turbulent flow
The results plotted in Figures 4 and 5 show as systematic underprediction of the friction factor over the entire range of turbulent conditions for the circular pipe, and over most of the range for the parallel plates. This is consistent with the results published by Etemad et al. [62], in which various LRN models were shown to underpredict the friction factor for a square duct at R = 5.6 × 10 5 . It is therefore of interest to assess why that bias is observed. Figure 6 presents the results of the LRN model with the thickness of the first set of layers (up to y + = 20) increased from ∆y + = 0.25 to ∆y + = 2. A consistent increase of the friction factor is observed for both the pipe and duct problems. This higher energy loss better matches the experimental curve for the circular pipe, as shown in Figure 7.
These results suggest that the increase in truncation error arising from using a coarser mesh might help in reducing the under-prediction of the friction factor. To investigate the reason behind this effect, u + and k profiles were plotted on Figure 8 for the two mesh sizes corresponding to the pipe problem, for R C = 10 6 . It is observed that the coarse mesh provides an insufficient resolution to accurately represent the TKE profile 1 < y + < 10 range, leading to a higher peak value in the transitional sublayer. This extra TKE requires an extra dissipation, and hence leads to a higher friction factor, which better matches the experimental data. This argument would then suggest that the Launder-Sharma LRN model would need some tuning in order to better represent the correct TKE profile on fine meshes.
The previous test also indicates that the results for the friction factor from the LRN treatment are very sensitive to mesh discretization in the transitional sub-layer,  besides the restriction for the location of the first node (which is perfectly fulfilled by both meshes). It seems that a higher computational efficiency (in terms of computing time) could be achieved if the mesh were densified without the need of moving the first node closer to the wall. In this respect, the collocated arrangement of variables used by OpenFOAM, among many other computer packages based on the FVM, is limited, as the computational nodes can at most be arranged with a fixed spacing equal to 2y 1 . Furthermore, if the cell size is geometrically increased, as per usual practice, the resolution on the transitional sublayer is further decreased.

Simulation of different Reynolds numbers with the same mesh
An interesting property of the LRN treatment is that, theoretically, if a sufficiently fine mesh is used it should be possible to solve a wide range of Reynolds conditions with the same mesh. This is important in cases where flow transients occur, and the discharge changes significantly over time. Figure 9 presents, for the pipe problem, a comparison between the results with the converged mesh for each Reynolds number, and the ones obtained with a single mesh for the whole Reynolds number range, namely, the converged mesh for the highest tested Reynolds number,

Model performance
Tests were carried out for the two geometrical problems using the LWF treatment, described in equations (21) to (23), with three different turbulence models: Standard k-ε, RNG k-ε and realizable k-ε.
The LWF treatment only applies when the flow is in the turbulent regime. Hence, the tests were performed for the range of Reynolds numbers above R = 3000, which is the approximate threshold at least for the circular smooth pipe problem.
For every test the mesh was generated so as the first node close to the wall lied at a nondimensional distance y + 1 = 30. From thereon, the cell size was increased geometrically with an expansion factor R = 1.1. Figure 10 shows the comparison of the numerical simulations with the experimental curves, indicating an overall satisfactory agreement. The two methodologies above described to solve for the shear velocity were tested for the standard k-ε model, i.e., iterating over y + and using the method by Launder & Spalding [36], but the results turned out to be identical, which could be expected due to the absence of flow separation in the test problems. Figure 11 presents the errors relative to the experimental curves. In the case of the circular conduit it is observed that the LWF treatment tends to under-predict the friction factor with the three turbulence models. The standard kε model shows the better overall performance, with differences below 4% for R C > 3×10 4 , and reaching a maximum of 7% around R C ≈ 10 4 . For the realizable and RNG turbulence models the errors are similar to those of the standard model at low Reynolds numbers, but significantly larger at higher R C .
For the case of parallel plates, the three models overpredict the friction factor for R P < 2 × 10 4 , with the realizable turbulence model showing the lowest errors. However, for higher Reynolds numbers the standard model performs much better, with a maximum error of around 3%.
In any case, as for the LRN treatment, these errors are very much smaller than the ones associated to physical models.

Resolution limit
Though the previous results indicate that the LWF treatment is suitable down to relatively low Reynolds numbers, even below 1 × 10 4 for the circular pipe, there exists a numerical resolution problem for that low range. For exam-  ple, for the pipe problem the spatial range of the y coordinate is limited by the restriction 30 ≤ y + ≤ r + , where r + = (︀ D/2 )︀ u * ν −1 is the nondimensional radius of the conduit. Figure 12 shows the lower (y + = 30) and upper (y + = 500) limits of the inertial sublayer in terms of the nondimensional coordinate y ≡ y/D, as a function of R. It is observed that below R = 2 × 10 4 the outer layer disappears, so the model domain extends completely within the inertial sublayer. As R decreases further, the lower limit of the inertial sublayer grows steadily, reducing the effective model domain. Figure 12 also shows the number of grid points in the cross section when using equal size cells with y + 1 = 30. As observed, the discretization becomes increasingly coarser when R decreases. For example, for R = 2 × 10 4 the entire cross section of the conduit ends up being discretized in only 20 cells. This quantity might be considered the minimum amount of grid points required to give a reasonable representation of the velocity profile, and hence the lower limit of R for the application of the LWF treatment.
Though for the particular problems solved in this paper (fully developed flow in conduits) the results were found to be reasonably accurate for Reynolds numbers below that limit, this is attributed to the simplicity of the problems.

Discussion on numerical model errors
As a guidance for 'best management practices', Table 2 presents a synthesis of the conditions applicable to RANS models of the k-ε family to simulate wall-bounded flows using LRN or LWF treatments. They include the restrictions and error bounds for the friction factor found in the present paper. The LWF treatment has been shown to be quite competitive in relation to the LRN treatment for high Reynolds numbers, as it provides similar errors for the friction factor while keeping a much coarser grid (hence requiring a much lower computing time). However, for R < 2 × 10 4 the requirement of locating the first node within the inertial sublayer leads to a substantial loss of resolution of the flow profiles.
Using these best practices, it has been shown that numerical models based on the RANS approach with k-ε turbulence models are able to represent the friction factor for a wide range of Reynolds numbers with errors below 10%. This means that for a typical hydraulic project, these models can accurately represent smooth wall frictional head losses for both the prototype scale and any reduced scale adopted for the physical modelling.

Methodology to avoid scale effects
It has been shown that the friction factor in a physical model amplifies relative to the one in the prototype, with increments from about 50% to 350%, then introducing significant scale effects. On the other hand, it has also been demonstrated that RANS models can be used as a relatively accurate tool to quantify the scale effects introduced by smooth wall friction on physical models, by numerically simulating the problem at both scales and comparing the results. Moving further, this allows the use of the numerical model to pre- dict prototype scale performance free of scale effect with relative accuracy, by applying the following methodology, outlined on Figure 13: 1. Implement a numerical model representing the reduced scale physical model. 2. Validate this reduced scale numerical model by comparing its results to those measured in the physical model. Different turbulence models should be tested, in order to find the one that best matches the problem. Eventually, some turbulence model parameters could be adjusted for better calibration. 3. Implement a numerical model representing the prototype at full scale, using the same turbulence model and parameters as the ones adopted for the physical model.
This procedure was already applied by the authors when dealing with the design of the filling/emptying system for the Third Set of Locks of Panama Canal [15], albeit without evaluating the errors of the numerical model.

Conclusions
From analyses undertaken for fully developed flow in smooth pipes and ducts in the range of Reynolds number from 10 4 to 10 8 , it has been shown that the friction factor in Froude-scaled physical models is significantly over-predicted, leading to the introduction of scale effects. Meanwhile, numerical models based on the RANS approach, using k-ε turbulence models, are able to represent the friction factor for both scales with errors below 10%, which can be considered as practically negligible from the engineering point of view. RANS models can therefore be used as a relatively accurate tool to quantify the scale effects introduced by smooth wall friction on physical models. In order to attain the reported accuracy, the computa-tional mesh must fulfil some conditions (synthetized in Table 2). In the case of the LRN treatment, a minimum cell size and a band of constant cell sizes within the inertial sublayer is recommended. If the LWF treatment is applied, the condition on the location of the first node close to the wall restricts both the model domain and the grid resolution, so a lower limit of 2 × 10 4 is recommended for the Reynolds number. Digging deeper into the characteristics of the numerical error, it has been shown that the consistent underprediction of the friction factor of Launder-Sharma LRN model for the full turbulent regime, could be explained as a misrepresentation of the peak value of TKE occurring in the inertial sublayer, suggesting that some tuning of this model would be necessary in order to better represent this effect. To attain results free of scale effects, a two-step methodology using numerical models is proposed. In the first place, the physical model scale is numerically simulated in order to validate the numerical treatment. In the second place, the prototype scale is numerically simulated with the validated numerical treatment.