Turbulence model performance for ventilation components pressure losses

This study looks to find a suitable turbulence model for calculating pressure losses of ventilation components. In building ventilation, the most relevant Reynolds number range is between 3×104 and 6×105, depending on the duct dimensions and airflow rates. Pressure loss coefficients can increase considerably for some components at Reynolds numbers below 2×105. An initial survey of popular turbulence models was conducted for a selected test case of a bend with such a strong Reynolds number dependence. Most of the turbulence models failed in reproducing this dependence and predicted curve progressions that were too flat and only applicable for higher Reynolds numbers. Viscous effects near walls played an important role in the present simulations. In turbulence modelling, near-wall damping functions are used to account for this influence. A model that implements near-wall modelling is the lag elliptic blending k-ε model. This model gave reasonable predictions for pressure loss coefficients at lower Reynolds numbers. Another example is the low Reynolds number k-ε turbulence model of Wilcox (LRN). The modification uses damping functions and was initially developed for simulating profiles such as aircraft wings. It has not been widely used for internal flows such as air duct flows. Based on selected reference cases, the three closure coefficients of the LRN model were adapted in this work to simulate ventilation components. Improved predictions were obtained with new coefficients (LRNM model). This underlined that low Reynolds number effects are relevant in ventilation ductworks and give first insights for suitable turbulence models for this application. Both the lag elliptic blending model and the modified LRNM model predicted the pressure losses relatively well for the test case where the other tested models failed.


Introduction
Technical pressure loss data of ductwork components is needed for the design, the optimisation and the energy assessment of ventilation systems. Measurements of pressure loss coefficients can be time-consuming and costly, especially for large duct dimensions. Ventilation ductworks often are custom designs for which data is lacking (Shao and Riffat 1995). In consequence, there were many attempts to substitute pressure loss measurements with computational fluid dynamics (CFD) simulations (Sleiti 2013). For this purpose, suitable turbulence models (see Section 1.3) are needed, which are the focus of this study.

Pressure loss coefficients
Pressure losses ΔP for straight ducts and components are usually expressed as non-dimensional pressure loss coefficients λ and ζ according to Eqs. (1) and (2) (Idelchick 2008 (Idelchick 2008): Figure 1 schematically shows the dependence of pressure loss coefficients on the Reynolds number. For low Reynolds numbers (laminar flow), pressure losses increase linearly with the velocity (ΔP~u) and ζ(Re) asymptotically approaches the form ζ → C 0 /Re (Wagner 2012). For high Reynolds numbers (fully turbulent flow), pressure losses increase quadratically with the velocity (ΔP~u) 2 and a constant value ζ → ζ ∞ is approached that depends on the wall roughness (Wagner 2012). In the intermediate transition region, the curve progression depends on the component and the inlet conditions (Idelchick 2008

Bends
In a rough estimate, 90% of the pressure losses in ventilation systems are typically due to components and only about 10% due to straight ducts. Within the category of components, bends are particularly important with a relative share of about 20% (Kriegel et al. 2018). They are also an often-found fluid dynamic disturbance in front of heat exchangers, fans, flow meters, and volume flow controllers. Sharp bends or deflections are also often found in air handling units, e.g., in the vicinity of cross flow heat exchangers or internal bypasses. In current practice, mainly constant pressure loss coefficients are used for bends, independent of the Reynolds number (VDI 2006;CIBSE 2007; ASHRAE duct fitting database). However, it is known that a Reynolds number dependency exists for bends especially in the range Re < 2×10 5 (Koch 2006;Idelchick 2008).
The transition region for curved channels and bends extends to higher Reynolds numbers than for straight ducts (Kalpakli Vester et al. 2016). In practice, it is often less well known that for bends, Reynolds number effects can be relevant up to relatively high values such as Re = 10 5 . Sprenger (1969) measured pressure loss coefficients of 17 variants of air duct bends in the Reynolds number range from 3×10 4 to 6×10 5 . The ducts had a rectangular cross-section of 200 mm × 100 mm. The pressure loss coefficients and their dependence on the Reynolds number were presented for different inner and outer curvature radii of the bends.
Technical tables, e.g., VDI (2006), often only list a single value ζ(Re = 2×10 5 ) and lack information on the dependence on the Reynolds number (Koch 2006). The work of Sprenger is a notable exception in this respect. Figure 2 shows measured  Sprenger (1969). Light markers indicate the full data set of all components values of ζ(Re) for four selected bends from the dataset. The full dataset is indicated by light markers. The largest variation of ζ(Re) occurred for the bend 25/05 (labelled by bend outer radius/inner radius in cm, see Figure 3) and the smallest variation occurred for the bend 20/00. The curve progressions depended most notably on the inner curvature radius of the bend.

Turbulence modelling
A formidable task for most CFD simulations is that of modelling turbulence. Examples of previous CFD studies on ventilation components are listed in Table 1. Depending on the research focus, the considered physical phenomena and the validation approach, the authors of above-mentioned studies selected different turbulence models. Criteria taken into consideration were, for example, the prediction of velocity distributions, turbulent quantities, pressure losses, heat transfer, secondary flows and rotation, flow separation and reattachment, model availability, prior performance in other applications and computational effort. Differences in priorities between these criteria led to diverging recommendations and criticism of different models, depending on the selected focus. As there exist a multitude of turbulence models and variants, some more explanations are useful on the selection of a turbulence model for the present task.
An important factor for air duct flow simulations is wall boundary layers. They determine the wall friction and surface heat transfer and thus pressure and heat losses. Presently, the computational effort of time-resolved simulations for wall-bounded flows at high Reynolds numbers is prohibitive and turbulence models are widely used for wallbounded flows (Durbin 2018). Two-equation turbulence models based on the turbulent kinetic energy k and the specific dissipation ω were specifically developed for this purpose (Wilcox 2006). In eddy viscosity models, the Boussinesq approach is used (Eq. (4)). A linear relation between the Reynolds stresses ij τ and the mean strain rate tensor (composed from gradients of the mean velocity u i ) is assumed, with the turbulent viscosity μ t as the parameter, viz.:  In the popular SST model (Menter 2009), μ t is modelled as μ t ~ k/ω. Transport equations are solved for k and ω.
Flows in bends are complex due to secondary flows and using the Boussinesq approach for the simulation of bends can be criticised (Kalpakli Vester et al. 2016). Nevertheless, eddy viscosity models are widely used in industrial applications, despite their known shortcomings (Argyropoulos and Markatos 2015). In non-circular ducts, Reynolds stresses cause secondary motions that cannot be simulated using conventional isotropic two-equation eddy viscosity models. Detailed DNS simulations of square duct flows are now available (Pirozzoli et al. 2018) for Re ≤ 4×10 4 that clarify the importance of this influence. The intensity of the secondary flows was only in the order of 1%-2% of the bulk velocity and their influence on the pressure loss was minor. Consequently, neglecting the secondary flows may be justified for straight rectangular ducts if the intention is the simulation of pressure losses. Another common criticism of two-equation models is their lack of sensitivity to streamline curvature and rotation. Modifications exist that attempt to improve the predictive capabilities for this influence, one of which (Arolla and Durbin 2013) was also tested in this work. As an alternative, Reynolds Stress Models (RSM) such as the EBRS (Lardeau and Manceau 2014) are available that solve for the six independent components of ij τ . In comparison to two-equation models, RSMs increase the computational effort and are known to be more challenging to solve numerically (Argyropoulos and Markatos 2015).

Low Reynolds number turbulence modelling
Conventional turbulence models tend to be most accurate when viscous stresses are much smaller than turbulent stresses. Viscous stresses can become relevant where turbulent fluctuations are relatively weak such as near walls. Besides walls, mean flow acceleration, deceleration and rotation can weaken or increase turbulent stresses (Durbin 2018) and thus cause a change of the flow regime. Near-wall modelling is an important factor when simulating surface heat transfer (Kriegel 2005;Mangeon et al. 2020) and flow separation from curved surfaces .
The Reynolds number can be interpreted as a ratio of inertial to viscous forces. It is also used as a parameter for the relative strength of turbulent stresses because they are driven by inertial forces of the main flow and tend to be damped by viscous stresses. Flows where viscous stresses are relatively important are therefore termed low Reynolds number flows. The definition of a Reynolds number is, however, case-dependent on the geometry and flow problem and the definition of a low Reynolds number flow can also be only case-dependent.
The SST model of Menter (2009) was used successfully in many engineering applications because it combines unique predictive capabilities of other models in a single formulation (Argyropoulos and Markatos 2015). It also ranks high in the prediction of pressure losses and heat transfer in internal flows. The predicted dependency of pressure loss coefficients on the Reynolds number was in some cases found unsatisfactory for lower Reynolds numbers in simulations of bends (Kriegel et al. 2018;Tawackolian et al. 2016). This motivated an investigation on how the SST model can be adapted or calibrated to better reproduce ζ(Re).
The first idea was to sensitise the turbulence model closure coefficients directly to the Reynolds number Re, specifically the production term in the transport equation of k. Although first promising results were obtained in this way, the limits of such an approach soon were apparent. The Reynolds number, as defined in Eq. (3) is a non-local parameter that is defined differently depending on the geometry. Such a definition becomes unsuitable for simulating complex geometries where a single global Reynolds number can often not be defined. In contrast to this, the differently defined turbulent Reynolds number Ret is a local parameter: This parameter can be calculated locally without additional information about the geometry such as a reference length D. Since Re t is a field quantity that varies in the computational domain, the formulation of closure coefficients in dependence of Re t becomes more intricate. Wilcox (1992Wilcox ( , 1994 proposed a modification of his k-ω turbulence model for low Reynolds number flows (LRN model) with damping functions dependent on Re t ,which is available in popular CFD software (ANSYS FLUENT, OpenFOAM, Siemens CCM+, etc.). The LRN model was mainly used for simulating aerofoils in aircraft and turbomachinery applications (Genç et al. 2012) and not for internal flows.
The formulation of the LRN model is described in the Appendix, which is in the Electronic Supplementary Material (ESM) of the online version of this article. It introduces three closure coefficients that were selected based on comparison with channel flow at a relatively low Reynolds number Re = 3300 and a transitional flat-plate boundary layer (Wilcox 2006). In the present work, a first test with the LRN model resulted in only a minor change of the ζ(Re) prediction for a selected test case (see Section 3). To the authors' knowledge, there is yet no published research where the closure coefficients of the original LRN model were adapted for specific applications different from the initial calibration, apart from cases where altogether different functional terms were proposed (Peng et al. 1997;Bredberg et al. 2002).
Later, a different modification on the production term in the transport equation of k (Eq. (A5), Appendix in the ESM) was proposed by Menter et al. (2015), again targeting transition on aerofoils. Instead of an algebraic damping function, a transport equation was introduced to obtain the damping coefficient. The term was introduced with different physical reasoning as intermittency γ where γ = 0 corresponds to laminar and γ = 1 to fully turbulent flow. The behaviour of a transport equation is more complex than an algebraic damping function and the calibration of this approach for new applications can therefore be expected to be more intricate. A transport equation intrinsically includes history effects, i.e., influences from upstream flow conditions. That is why it can be expected to be superior in predicting extended regions of transition from laminar to turbulent flow that can occur on aerofoils and flat plate boundary layers. Abraham et al. (2019) suggested that the transition model of Menter can be used for internal flows if different closure coefficients are used.
As summarised by Patel et al. (1985) and Hrenya et al. (1995), many variants of k-ε models were proposed that incorporate damping functions based on Re t or the wall distance. Durbin (1991) introduced an elliptic blending parameter f for near-wall modelling, obviating algebraic damping functions. An additional transport equation is solved to introduce a parameter for the anisotropy of the Reynolds stress tensor. The approach was modified by Davidson (2003) and evolved into the elliptic-blending model by Billard and Laurence (2012). Later, in lag-elliptic blending models (Lardeau and Billard 2016), a lag parameter φ was used. A transport equation is solved for φ and an elliptic function for a parameter α that is used for blending in dependence on the wall distance. Revell et al. (2011) derived the lag parameter, based on the comparison of the Boussinesq model with exact terms from the Reynolds stress transport equations. The lag-elliptic blending approach was recently also applied to the k-ω model ) and the SST model (Shang and Agarwal 2020).

Geometry
The dimensions of the components investigated here are based on an ongoing measurement campaign at our institute. The measurements are based on the Reynolds-similarity principle. The dimensions were scaled down and the velocity was correspondingly increased to obtain the Reynolds number range of interest. In this way, pressure losses are more easily measured in the low Reynolds number range Re < 2×10 5 .
A rectangular duct with internal dimensions of 150 mm × 75 mm (hydraulic diameter D = 100 mm) was used. Four bends were selected from the study by Sprenger (1969): 25/1.4, 25/05, 25/10 and 30/00 (Figure 3). They have an outer radius of 25 cm or 30 cm and inner curvature radii of 1.4 cm, 5 cm, 10 cm or 0 cm respectively. Case 25/05 is a simple bend with a curvature radius R k /D of 0.75 and was used as a reference for the calibration of the LRN model. Since the duct dimensions of the reference were 200 mm × 100 mm all geometric dimensions were scaled by a factor of 0.75 for the corresponding simulations.

Mesh
Meshes were created using the trimmer mesher in Siemens CCM+ with ten prism layers at the wall. The wall distance of near-wall cell centres was 5 μm corresponding to a nondimensional wall distance y + of 1 for the highest Reynolds number 5×10 5 . The inlet and outlet surfaces of the bends were extruded 1 m (10D) upstream and 8 m (80D) downstream of the component. Based on an initial mesh study of case 25/05 (see Appendix in the ESM and Figure 4) the mesh discretisation was selected. The resulting cell count for each case is listed in Table 2. Coarse meshes (1.1 m cells) were used for the optimisation studies of the closure coefficients for the LRN model (Section 2.7 and Section 3.3) because of the high computational effort. Simulations of one test case for all Reynolds numbers on the fine mesh required approximately 4 days on a server node with 16 CPUs.

Fluid physics and boundary conditions
Computational fluid dynamics (CFD) simulations were carried out using Siemens STAR-CCM+ version 15.04 (Siemens 2020). The main simulation parameters are listed in Table 3. They were selected in accordance with the planned experiments.
The air was modelled as a constant density and constant temperature Newtonian fluid at a temperature T = 26 °C. Heat transfer was neglected. Fully developed flow conditions were described at the inlet. They were generated with precursory simulations of a straight duct segment and periodic boundary conditions for each turbulence model, respectively. Fixed pressure boundary conditions were used at the outlet. All walls (upper, bottom, lateral) were defined as no slip boundaries and were considered as hydraulically smooth. Steady RANS simulations were carried out for 5000 iterations and the convergence was checked by monitoring the pressure drop between the inlet and the outlet of the domain.

Turbulence model survey
Several popular turbulence models were considered for a survey to identify their performance in simulating ventilation components. The models were all available through Siemens CCM+ and are listed in Table 4. The list includes two Reynolds Stress Models (EBRS and LPRS) and 13 eddy viscosity models. Case 25/05 was selected for the survey because of its distinct dependence on ζ(Re).

Straight duct simulation
In the first step, straight duct simulations were carried out to generate inlet boundary conditions for the main simulations. Periodic boundary conditions were used on a duct segment of 0.5 m (5D) to simulate fully developed flows. A structured mesh of 75 × 100 cells in the two cross-stream and 25 cells in the streamwise direction was used. The near-wall spacing was 5 μm, corresponding to a non-dimensional wall distance y + of 1 for the highest investigated Reynolds number 5 × 10 5 .

Calculation of pressure loss coefficients
In analogy to the measurement procedure according to ASHRAE (2017)

Optimisation study for the LRN model
The design manager of Siemens CCM+ was used to carry out an optimisation study with the three closure coefficients of the LRN model as input parameters. Case 25/05 was again used as a reference because of the distinct variation of ζ(Re) for this case in the experimental dataset. For each realisation, simulations at the four Reynolds numbers 4×10 4 , 1×10 5 , 3×10 5 and 5×10 5 were carried out. These Reynolds numbers were selected based on the slope of the reference curve ( Figure 2). The pressure loss coefficients were calculated from the simulation results and the deviations from the experimental reference values were determined for each Reynolds number. The target function for the optimiser was to minimise the weighted average of the deviation at the four Reynolds numbers by varying the closure coefficients.
Since there was only one reference value at the lowest Reynolds number 4×10 4 , the deviation at this point was weighted by a factor of two and the other deviations were weighted equally by a factor of 1. The SHERPA algorithm was used for the optimisation (Red Cedar Tech. 2008) which automatically selected between multiple optimisation algorithms during runtime to vary the input parameters. In total, 200 simulations were carried out during the optimisation run. Figure 5 depicts the calculated friction coefficients λ in dependence on the Reynolds number for selected models. For comparison, the Kármán-Prandtl friction law for pipe flows and the friction law given by Schultz and Flack (2013) for turbulent channel flow is included. The Kármán-Prandtl friction law is valid for smooth pipes but was also found to apply to turbulent duct flows if the hydraulic diameter is used (Leutheusser 1984). All tested models had a similar prediction of the friction coefficient, mostly within ±10% of the Kármán-Prandtl friction law. The wall damping formulation of the LRNM model slightly influenced the prediction of λ in comparison to the SST model. The friction coefficient prediction for straight ducts is of minor practical importance for air ductworks since it is often considered to be known and the influence due to the low Reynolds modification was relatively small here.  Sprenger (1969) and small markers indicate corresponding error bounds of ±10%. Most conventional models did not correctly predict the curve slope in the transition region ( Figure 6) that takes up nearly the full shown Reynolds number range for this case. Instead, they predicted a flatly descending curve progression that was only valid for a fully turbulent flow at higher Reynolds numbers, outside the depicted range. Generally, due to the large differences in the predicted curve progressions, a match of ζ(Re) between simulations and experiments was only occasional. The difference between LRN and SST was minor. The case LRNM depicts the LRN model with different closure coefficients that were obtained from the optimisation study. With the new coefficients, a similar prediction to EBL was obtained. The curvature correction (Arolla and Durbin 2013) only had a marginal influence on the results of the SST model (SSTCC). The curve progressions of RKE and RKEN differed markedly from the measurements. SKELR agreed well at low Reynolds numbers. Further tests with SKELR which are not shown here showed that it was sensitive to the near wall mesh spacing. LPRS had initial convergence issues and it was initialised with a solution of EBRS.

Initial turbulence model survey
A quantitative comparison of RMS deviations is shown in Figure 10. Based on the availability of experimental data, RMS deviations of ζ(Re) between the experiments and simulations were calculated for 3×10 4 < Re < 5×10 5 . Based on this metric, LPRS and LRNM ranked highest and RKE(N) ranked lowest in the prediction performance of ζ(Re). A limitation of this metric is that it only includes absolute deviations and not the slope of ζ(Re). KOM therefore ranked high at the third position, although the curve slope differed from the measurements.

Modification of the LRN model
The three closure coefficients which proved to be the best variants by the optimisation study are listed in Table 5.

Validation simulations
To test the modified LRN model it was applied to the other representative cases from the dataset from Sprenger (1969). For case 30/00 the sharpness of the bend inner corner is an uncertainty factor and a sharp edge was assumed in the simulation. The results are shown in Figure 11. The blue curves depict the simulation results obtained with the new modified closure coefficients. An improvement of the predicted curve slope in comparison to the SST model is noticeable not only for calibration case 25/05 but also for the other test cases. A quantitative comparison is given in Figure 12 that shows RMS deviations of ζ(Re) predictions for all validation test cases.

Discussion
The damping functions of the original LRN model only resulted in a minor influence on the predicted pressure loss coefficients. This may explain why this modification was not widely used for internal flows. The prediction ζ(Re) was better with the LRNM and the EBL model than with SST and RKE(N). For the LRNM model, it is encouraging that an improved prediction was not limited to the calibration test case but also noticeable for the other validation test cases. This is an indication of the function and importance of the near-wall model for these test cases. The closure coefficients of the LRNM model differ by a large factor (3.9 to 7.3) from the original LRN model. This may raise doubts about their validity for general applications and motivates a closer inspection.

Further work
The parameter optimisation study for the new closure coefficients of the LRNM model was only conducted for the single test case 25/05. It should be extended to include a more comprehensive set of test cases to cover the relevant range for ventilation components. The present results indicated that the prediction of ζ(Re) was improved with turbulence models which include near-wall modelling. Different functional forms of damping functions were proposed in the past and it is currently unclear which of these is most suitable for the present application. A study with the goal of investigating the best functional formulation of the blending function seems interesting. Considering different functional forms would make a parameter study more time-consuming. In the models EBL and GTR, the role of the damping function is replaced by differential equations. Closure coefficients also occur in the differential equations and it should be assessed if these can also be adapted. Since elliptical lag model variants were proposed for the k-ω ) and the SST (Shang and Agarwal 2020) model it will be interesting to also evaluate their performance for the present test cases.
In two-equation models, a ratio of quantities is used for the calculation of the turbulent viscosity (e.g., t~/ μ k ω). Therefore, the assessment of near-wall models can be under-determined if solely a statistic mean quantity, such as here the pressure loss, is considered as a criterion. Similar solutions are obtained if k and ω change by the same factor and the individual values can be ambiguous. It may therefore be useful to include additional statistic quantities in the optimisation target function that correlate specifically with k or ω. It is important that the reference data used for calibration is reliable and that the uncertainty is acceptably small. In the past, canonical DNS simulations were the primary choice for such reference data. We wish to include our own high precision measurements specifically for calibration. Experimental data for the selected test case was available for Re > 3×10 4 . It will be interesting to have measurements also for slightly lower Re due to the relatively large slope of ζ(Re) at low Re for the validation of the turbulence model. The present choice of case 25/05 for the calibration was due to its distinct variation of ζ(Re). Flow separation at the inner wall of the bend was important for this case. The simulation of separated flows is a separate challenging topic for RANS models. It may therefore be useful to find a calibration reference case that has strong viscous effects but is not affected by flow separation to facilitate the calibration. Once a validated model is available, it can be applied to more ductwork components where current data on ζ(Re) is lacking.

Conclusions
Viscous effects at low Reynolds numbers influence the prediction of pressure loss coefficients ζ(Re) of duct components. Caution should be taken as not all turbulence models are suitable in this regime and this can lead to incorrect curve progressions of ζ(Re) which are not applicable for lower Reynolds numbers. Both the EBL model and the LRNM model performed well for the simulation of several test cases. They are presently considered as the best candidates for the calculation of pressure losses in ventilation systems.