Low-Reynolds number mixing ventilation flows: Impact of physical and numerical diffusion on flow and dispersion

Quality assurance in computational fluid dynamics (CFD) is essential for an accurate and reliable assessment of complex indoor airflow. Two important aspects are the limitation of numerical diffusion and the appropriate choice of inlet conditions to ensure the correct amount of physical diffusion. This paper presents an assessment of the impact of both numerical and physical diffusion on the predicted flow patterns and contaminant distribution in steady Reynolds-averaged Navier–Stokes (RANS) CFD simulations of mixing ventilation at a low slot Reynolds number (Re≈2,500). The simulations are performed on five different grids and with three different spatial discretization schemes; i.e. first-order upwind (FOU), second-order upwind (SOU) and QUICK. The impact of physical diffusion is assessed by varying the inlet turbulence intensity (TI) that is often less known in practice. The analysis shows that: (1) excessive numerical and physical diffusion leads to erroneous results in terms of delayed detachment of the wall jet and locally decreased velocity gradients; (2) excessive numerical diffusion by FOU schemes leads to deviations (up to 100%) in mean velocity and concentration, even on very high-resolution grids; (3) difference between SOU and FOU on the coarsest grid is larger than difference between SOU on coarsest grid and SOU on 22 times finer grid; (4) imposing TI values from 1% to 100% at the inlet results in very different flow patterns (enhanced or delayed detachment of wall jet) and different contaminant concentrations (deviations up to 40%); (5) impact of physical diffusion on contaminant transport can markedly differ from that of numerical diffusion.


Introduction
Computational fluid dynamics (CFD) can be a powerful tool to analyze complex indoor airflows. The use of CFD for ventilation problems has increased tremendously over the last decades (e.g. Chen 2009, Blocken 2014. CFD is-among others-used for natural ventilation studies (e.g. Kato et al. 1992;Jiang and Chen 2002;Kurabuchi et al. 2004;Hu et al. 2005;van Hooff and Blocken 2010a,b;Novoselac 2011, 2013;Ramponi and Blocken 2012a,b;Ai and Mak 2014a,b;Perén et al. 2015;Tong et al. 2016a,b), for indoor airflow studies (e.g. Gan and Awbi 1994;Chen 1995;Nielsen 1998;Zhang et al. 2007; Wang and Chen 2009;Cao and Meyers 2013;Liu et al. 2013Liu et al. , 2016aYou et al. 2016) and for indoor pollutant concentration studies (e.g. Chung and Hsu 2001;Rouaud and Havet 2005;Chen et al. 2014Chen et al. , 2015. Quality assurance is vital for accurate and reliable CFD simulations and verification and validation studies are therefore imperative. The Special Interest Group on Quality and Trust in Industrial CFD of the European Research Community on Flow, Turbulence and Combustion (ERCOFTAC) published general guidelines on the use of CFD by industry (Casey and Wintergerste 2000). For the modeling of indoor airflow (ventilation flows) several similar guidelines as well as journal publications on specific topics have been published (e.g. Jones and Whittle 1992;Chen and Srebric 2002;Sørensen and Nielsen 2003;Nielsen 2004;Nielsen et al. 2007).
In CFD simulations diffusion plays an important role. A distinction can be made between physical diffusion, consisting of (1) laminar (molecular) diffusion and (2) turbulent diffusion, and numerical diffusion-also called artificial, or false diffusion, or artificial viscosity, -which is caused by truncation errors due to the discretization of the governing flow equations. Numerical diffusion is related to the spatial resolution of the computational grid and the employed spatial discretization schemes. The effect of physical and numerical diffusion is similar, since both can result in reduced gradients of velocities, turbulence parameters, and scalars such as temperatures and mass concentrations, also reported as "smeared out" predictions (e.g. Anderson 1995;Ferziger and Peric 1996;Nielsen 2004;Nielsen et al. 2007;Versteeg and Malalasekera 2007). Therefore, CFD users should apply caution in selecting applying appropriate computational grids, spatial discretization schemes and turbulence boundary conditions.

Grid resolution
The grid resolution determines to a large extent the amount of numerical diffusion, with coarser grids leading to a higher level of numerical diffusion (e.g. Skovgaard and Nielsen 1990;Jones and Whittle 1992;Ramponi and Blocken 2012a). In 1992, Jones and Whittle (1992) stated in their overview paper that sufficient grid resolution should be applied in jet regions, boundary layers and in thermal plumes, and that: "A failure to provide enough mesh in these areas will result in the supply jet or the boundary layer flow being insufficiently resolved (where the momentum is rapidly and artificially diffused) resulting in an unrealistic local (and possibly global) solution.". In addition, they explicitly mentioned the increased presence of numerical diffusion when the flow direction is not aligned to the wall surface in case of a structured grid. Li and Nielsen (2011) stated that in theory, numerical errors can be limited to an acceptable level when the number of grid points is sufficient and the grid resolution is sufficiently high. In addition, several previously published guidelines stated that a grid-sensitivity analysis should be conducted to make sure the results obtained are (nearly) grid independent and do not suffer from excessive amounts of numerical diffusion (e.g. Jones and Whittle 1992;Chen and Srebric 2002;Chen and Zhai 2004;Roache 1997;Li and Nielsen 2011;Blocken 2015). Sørensen and Nielsen (2003) stated that-at that time-obtaining a grid-independent solution for 3D cases was very difficult due to the lack of computational resources. Therefore, they recommended to perform tests to make sure at least grid convergence is present; i.e. "further refinement of the computational grid will change the solution, but not significantly" (Sørensen and Nielsen 2003).
In general, a systematic grid refinement with a constant factor is recommended to find the grid that provides gridindependent results, or at least to prove that grid convergence is present. The possible errors when a grid-sensitivity study has not been performed with sufficient scrutiny, or when only grid convergence has been obtained, can be important however and are studied and presented in this paper.

Spatial discretization schemes
Best practice guidelines (e.g. Casey and Wintergerste 2000;Chen and Srebric 2002;Nielsen 2004;Nielsen et al. 2007) recommend at least the use of second-order spatial discretization schemes to limit the amount of numerical diffusion. Note that hereafter "spatial" will be omitted when referring to discretization schemes. Information on the diffusive character of first-order upwind (FOU) schemes and the impact on the numerical solution was already published in the 1970s and 1980s (e.g. Leonard and Mokhtari 1990). Jones and Whittle (1992) indicated the influence of the discretization scheme on the results of CFD simulations of ventilation flows and mentioned the potential of higherorder schemes to limit numerical diffusion. Casey and Wintergerste (2000) recommended to avoid FOU schemes and to use at least second-order schemes for all transported quantities. Chen and Srebric (2002), Awbi (2003), Sørensen and Nielsen (2003), Nielsen (2004Nielsen ( , 2009), Ramponi and Blocken (2012a,b), and Goethals and Janssens (2013) also demonstrated the potentially detrimental impact of FOU schemes. For example, Sørensen and Nielsen (2003) presented results of numerical simulations of the Smith and Hutton problem, which resembles a recirculating flow as often encountered in indoor airflows. They mentioned that FOU schemes can accurately predict the flow as long as it is aligned with the grid lines, which however is hardly the case in ventilation flows in which recirculation zones, detachment of wall jets and thermal plumes are often important flow features. As a result, they counter advised the use of first-order schemes. Also the REHVA guide book by Nielsen et al. (2007) counter advises on the use of FOU discretization schemes based on a study of the impact of FOU schemes on the prediction of the flow from a wall-mounted opening. Ramponi and Blocken (2012a,b) demonstrated the detrimental impact of FOU schemes for the prediction of cross-ventilation flows, where the numerical diffusion led to flow patterns strongly deviating from the wind-tunnel experiments by Karava et al. (2011), while second-order upwind (SOU) schemes provided a very close agreement with those experiments. The potentially negative impact of FOU schemes on the simulation accuracy is also recognized by several scientific journals (e.g. Journal of Fluids Engineering, AIAA journal (for "spatially smooth solutions"; AIAA 2014)), who do not consider articles presenting numerical fluid flow studies using FOU discretization schemes for publication (e.g. Freitas 2002).
However, despite all guidelines and recommendations issued in the past, still new publications keep appearing in the scientific literature in which FOU discretization schemes are used, which is sometimes justified by using a highresolution grid or by an observed lack of convergence using SOU or other higher-order schemes. The errors made with FOU depend on the grid resolution and the flow being studied. To the best knowledge of the authors, the impact of the chosen discretization schemes in combination with the grid resolution employed on the calculated velocities and concentrations for low-Reynolds number ventilation flows has not yet been investigated in detail.

Physical diffusion
Past publications (e.g. Sørensen and Nielsen 2003;Ramponi and Blocken 2012a) reported that the effect of numerical (artificial) diffusion can be compared to the physical diffusion caused by the turbulence in CFD simulations. Both physical and numerical diffusion smear out the results and a high level of numerical diffusion will negatively impact the numerical accuracy. Ramponi and Blocken (2012a) showed that the proper choice of the inlet TI was equally important as the choice of a proper grid resolution and higher-order discretization scheme. In the past, some numerical studies focused on the effect of inlet TI on the flow pattern. Saïd et al. (1993) performed CFD simulations of a mixing ventilation flow for the geometry of the International Energy Agency (IEA) Annex 20 case (Nielsen 1990) and found that a variation in inlet TI between 4% and 37.4% did not influence the resulting velocity pattern inside the enclosure. Joubert et al. (1996) performed CFD simulations for the same case and for the same range of TI (between 4% and 37.4%) and also found no effect on the results. However, Abdilghanie et al. (2009) assessed the influence of the inlet TI for a generic enclosure and found local differences of more than 100% for the mean velocities and the velocity fluctuations when comparing the results for inlet turbulence intensities of 5% and 13%. Cao and Meyers (2013) studied mixing ventilation flow in a cubical enclosure and found a clear effect of inlet turbulence conditions on the flow pattern calculated; detachment of the incoming wall jet was delayed with increasing TI (range between 5% and 30%). The effect of the turbulence intensity imposed at the inlet appears to be case/flow dependent and requires additional research.

Objective and scope of the present paper
Given the lack of knowledge on the impact of numerical and physical diffusion on mixing ventilation flow as outlined above, this paper presents 3D steady RANS CFD simulations of indoor airflow in which the impact of several grid resolutions, discretization schemes, and inlet turbulence intensities on the results is analyzed in detail. To the best knowledge of the authors this is the first paper that systematically addresses the impact of both numerical diffusion (grid resolution and spatial discretization scheme) and physical diffusion (inlet turbulence intensity) for low-Reynolds number ventilation flows. In addition, the analysis includes both velocities and contaminant concentrations. CFD simulations of transitional ventilation flow (Re ≈ 2,500) are performed which are compared with experimental data from particle image velocimetry (PIV) measurements in a reduced-scale water-filled model (van Hooff et al. 2012a). The computational models are presented in Section 2. The results of the verification study and the validation study using PIV measurements are shown in Section 3, after which the analyses of the impact of numerical diffusion (Section 4), and physical diffusion (Section 5) are presented. A discussion (Section 6) and conclusions (Section 7) conclude this paper.

Computational geometry and grid
The computational model is a replica of a part of the experimental setup used in previous publications by the authors (van Hooff et al. 2012a(van Hooff et al. , 2014 (Fig. 1a). The cubical enclosure (part of the model) has dimensions 0.3 m × 0.3 m × 0.3 m (L × W × H) and an inlet height of hinlet/L = 0.1. To simplify the computational geometry the large conditioning section upstream of the contraction is not included; only the contraction itself is modeled (Fig. 1a). The outlet is extended in the x-direction to enhance convergence of the simulations (no or limited recirculating flow at outlet). The surface grid-extrusion technique by van Hooff and Blocken (2010a) is employed to create a fully structured computational grid. The resolution is increased in regions with high velocity gradients, i.e. the boundary layers near the wall and the shear layer of the wall jet. Five different grids are constructed using a linear refinement factor of about 2 in each direction. The total number of cells for each of the grids is: (1) coarsest grid; 84,669 cells; (2) coarse grid; 303,924 cells; (3) basic grid; 664,240 cells; (4) fine grid; 1,899,184 cells (Fig. 1b); (5) finest grid; 5,314,560 cells (Fig. 2). The dimensionless wall distance (y*) for all grids is well below the recommended value of y* < 5 to enable low-Reynolds number modeling (LRNM), which implies solving the flow all the way down to the viscous sublayer. For example, the value of y* at the top surface for the coarsest grid is between 0.025 < y* < 0.6. Note that the number of cells is different from those in previous numerical studies (van Hooff et al. , 2014 due to the inclusion of two local contaminant sources inside the computational domain.

Boundary conditions
The boundary conditions are chosen to replicate those of the experiments as closely as possible. The surfaces are modeled as smooth no-slip walls. A uniform velocity is imposed at the CFD inlet (see Fig. 1a), which was based on the Reynolds number at the actual ventilation inlet during the experiments and the ratio between the height of the CFD inlet (h/L = 0.3) and the actual ventilation inlet (h/L = 0.1); Uinlet;CFD = 0.025 m/s for the Reynolds number used in this study: Re ≈ 2,500. The turbulence parameters are specified based on the hydraulic diameter and the turbulence intensity (TI = 21%). The measured turbulence intensity (u RMS /U M ) in the wall jet region at x/L = 0.2 was around 3%-4% for each value of Re (van Hooff et al. 2012a). Due to the contraction the resulting turbulence intensity at the entrance of the cubic test section corresponded with the measured values when, after an iterative process, TI = 18% was imposed by van Hooff et al. ( , 2014. Note that in the current paper TI = 21% is imposed, resulting in a slightly better agreement for mean velocities in the wall jet and the remainder of the enclosure than with TI = 18%. Zero static pressure is imposed at the outlet. To study the impact of grid resolution, chosen discretization scheme, and TI imposed at the inlet, on contaminant transport inside the ventilated enclosure, two passive contaminant sources are introduced in the vertical center plane (z/L = 0.5) of the enclosure (see Fig. 1a). The middle points of the cubical sources with edge length xp/L = 0.016 are located at: Source 1: x/L = 0.33; y/L = 0.33; Source 2: x/L = 0.93; y/L = 0.33. The sources emit a contaminant with a density equal to water (ρ c = 998.3 kg/m 3 ) at a rate S c = 0.1 kg/(m 3 ·s). The volume of the sources is equal to 1.25 × 10 −7 m 3 . As in the experiments, the setup is a reduced-scale model filled with water.

Solver settings
The 3D steady RANS equations are solved using ANSYS Fluent 12 (2009). The low-Reynolds number k-ε turbulence model by Chang et al. (1995) is used for this study due to its good performance in previous studies (van Hooff et al. , 2014. Pollutant concentrations are obtained using an advection-diffusion equation (Eulerian approach) and turbulent mass transport is calculated using the standard- by Sc t =  t /D t , and is taken equal to Sc t = 0.7.
Pressure-velocity coupling is taken care of by the SIMPLEC algorithm, pressure interpolation is second order and SOU discretization schemes are used for both the convection terms and the viscous terms of the governing equations. In addition, simulations are conducted using FOU and QUICK (Leonard 1979) discretization schemes for the convection terms. Note that ANSYS Fluent by default uses second-order accurate schemes for the viscous terms of the governing equations (ANSYS 2009). Convergence is monitored carefully and is declared when the scaled residuals do not show any further reduction with increasing number of iterations. For a limited number of simulations oscillatory convergence is present. As discussed earlier in several publications addressing oscillatory convergence (e.g. Blocken 2015), an "average" solution is obtained by averaging the results over a sufficient amount of iterations (difference when averaging longer becomes less than 1%), over an interval of 1,000-10,000 iterations, depending on grid size, discretization scheme and inlet TI.

Grid-sensitivity analysis
A grid-sensitivity analysis is conducted to obtain a reference case with results that are (nearly) grid independent. In total five grids are used for the grid-sensitivity analysis (Fig. 2). The grids are compared based on the dimensionless mean streamwise (U/UM) and vertical (V/U M ) velocity along three lines in the vertical center plane, located at x/L = 0.2, x/L = 0.5 and x/L = 0.8 (see Fig. 2f). U M is the local maximum streamwise velocity along the vertical lines and thus varies for each of the three lines. Figure 3 shows that the differences for U/U M between the five grids are quite small, with the main differences present at x/L = 0.8 (Fig. 3c). The profiles of V/U M indicate that only the fine grid provides nearly grid-independent results (Figs. 3d-f). Note that the horizontal axis range in Figs. 3d-f varies for each location to increase the readability of the graphs. Figure 4 shows the results of the grid-sensitivity analysis for the dimensionless concentrations K for two contaminant source locations. K is calculated from with S c the contaminant source rate (= 0.1 kg/(m 3 ·s)), V the volume of the contaminant source (=1.25 × 10 −7 m 3 ), h inlet  Figure 4 shows that there are large differences in K between the five grids, especially for source location 1 (Figs. 4a-c). For source location 2 (Figs. 4d-f) the results from the finest grid and the fine grid are nearly identical.
To quantify the grid-sensitivity results, the gridconvergence index (GCI) by Roache (1997) is used to determine the error when the fine grid is used: with F s the safety factor, taken to be 1.25, which is the recommended value when three or more grids are considered in the grid-sensitivity analysis (Roache 1997), r the linear grid refinement factor (r = 2 ), p the formal order of accuracy which is assigned a value of 2 due to the SOU discretization schemes used for the reference case simulations, and f 2 and f 1 the solutions obtained on the fine and finest grid, respectively, which correspond to U/U M , V/U M or K.
Note that the formal order of accuracy is used to calculate the GCI as suggested by Roache (1997). Celik et al. (2008) suggested using the observed (apparent) order of accuracy. However, a very slight non-monotonic convergence due to oscillations around the expected grid-independent solution is present, i.e. ε 32 /ε 21 < 0, with ε 32 = f 3 − f 2 and ε 21 = f 2 − f 1 , which inhibits the calculation of the observed order of accuracy. Figure 5 illustrates the results on the fine grid with indication of the GCI. Table 1 lists the average GCI value and the maximum GCI value along the three lines. The average GCI values for U/U M and V/U M are very low; below 1.2% at x/L = 0.2 and x/L = 0.5, and below 2.85% at x/L = 0.8. The maximum values are below 15.71% and 9.64% for U/U M and V/U M , respectively, with the largest values present at x/L = 0.8. Also the average values for K for source location 2 are very low; below 0.33% at all three vertical lines, while larger values are present for K for source location 1; up to 12.22%. The maximum GCI value for source location 1 is 26.16% (at x/L = 0.5) and for source location 2 it is 2.06% (at x/L = 0.5).   maximum GCI values in Tables 1 and 2, the fine grid is retained in the validation study in the next section.

Validation
To establish a benchmark for the comparison of the different discretization schemes and inlet TIs, the results of the reference case (fine grid, SOU, TI = 21%) are compared with the PIV measurements by van Hooff et al. (2012a).

Experimental setup
The reduced-scale model consisted of a water column, a conditioning section and a cubic test section with edge length 0.3 m. A 2D PIV system was used to conduct the measurements. It consisted of a Nd:Yag (532 nm) doublecavity laser (2 × 200 mJ, repetition rate < 10 Hz) used to illuminate the field of view, and one CCD (Charge Coupled Device) camera (1376 pixel × 1040 pixel resolution, 10 frames/s) for image acquisition. The laser was mounted on a translation stage and was positioned above the cubic test section to create a laser sheet in the vertical center plane of the cube (z/L = 0.5); the camera was positioned perpendicular to the target area. Seeding of the water was provided by hollow glass micro spheres (3M; type K1) with diameters in the range of 30 -115 μm. More information on the experimental setup can be found in van Hooff et al. (2012b). Two sets of PIV measurements were performed. The first set focused on the entire vertical center plane of the cube, i.e. a target area of 0.3 m × 0.3 m (= ROI1) (Fig. 6). The second set focused on a smaller target area of 0.18 m × 0.12 m (W × H) in the proximity of the inlet, enabling a higher measurement resolution (= ROI2) (Fig. 6). The required measuring frequency was estimated from the integral length scale (= inlet height) and the characteristic velocity (i.e. inlet velocity) and was set at 2 Hz. Each measurement set consisted of 360 uncorrelated samples (i.e. averaging time = 180 seconds). For the validation study in this paper the timeaveraged velocities along three vertical lines inside the enclosure are used (Fig. 2f). Note that the results for y/L < 0.05 are not used because they are less accurate due to reflections of the laser sheet on the glass bottom of the cube. More information can be found in van Hooff et al. (2012a,b).

Results
with P i the predicted (CFD) value, O i the observed (measured) value and n the number of measurement locations. An optimal score for FAC2 and FAC1.5 is equal to 1. Along  Figure 9 shows the dimensionless concentration coefficient K (see Eq. (1)). Clear deviations are present between the results at x/L = 0.2 and x/L = 0.5 for source location 1 (Figs. 9a-c); the difference at mid-height between K obtained with FOU and SOU schemes is 14% and 10% at x/L = 0.2 and x/L = 0.5, respectively. At x/L = 0.8, the largest deviations are present in the region close to the top of the enclosure; at y/L = 0.8 the deviation is around 100%, although it must be noted that the absolute values of the concentration are low in this area. The results obtained with SOU schemes are generally within 2% of those obtained with QUICK schemes. Figures 9d-f show that differences between the results by FOU versus SOU and QUICK are around 100% over a large part of the height of the enclosure along all three lines (between 0.2 < y/L < 0.85 at x/L = 0.2 and x/L = 0.5, and between 0.2 < y/L < 0.65 at x/L = 0.8). The fact that differences for source location 2 are much larger than for source location 1 can be attributed to the fact that source 2, in contrast to source 1, lies in an area with high velocity gradients (edge of downward directed jet flow), and the velocity gradients are strongly affected by the discretization scheme employed. Excessive numerical diffusion by FOU leads to lower velocity gradients, which lead to a lower production of turbulent kinetic energy, a decrease in turbulent viscosity (t = ρC μ f μ k 2 /ε; see (Chang et al. 1995)) and thus also to a decrease in turbulent mass fluxes due to a decrease in turbulent mass diffusivity (D t =  t /Sc t ), compared to SOU or QUICK. Apparently, this decrease in turbulent kinetic energy and thus turbulent viscosity is stronger than the effect of the increase in numerical (artificial) viscosity. The limited amount of turbulent mass transport for FOU in the area of the contaminant source results in less contaminant transport to the recirculation zone and thus lower contaminant concentrations. Figure 3 (Section 3.1) indicated a close agreement between the results for U/UM on the coarsest, coarse and basic grid. If the grid-sensitivity would include only these three grids and if a decision would be made based on the streamwise velocities only, one might conclude that sufficient grid convergence can already be achieved using the coarsest grid. To illustrate the importance of a proper grid-sensitivity analysis, and to show the impact of the discretization schemes on the results obtained on coarser grids, this section presents the results of the flow field and the contaminant distribution obtained with the three different discretization schemes on the coarsest grid (84,669 cells). Figure 10 displays the dimensionless velocities along the three vertical lines obtained on the coarsest grid. As expected, the impact of the chosen discretization scheme becomes larger for coarser grids. Along each of the three vertical lines the FOU results clearly differ from the SOU results, whereas the SOU results are nearly identical to those by QUICK. Figure 10c shows that FOU on the coarsest grid predicts a delayed detachment of the jet; i.e. y/L = 0.90 for the maximum value of U/U M instead of y/L = 0.83. Figure 10c shows that for this particular case the impact of FOU versus SOU discretization scheme is larger than that of the grid resolution (compare "SOU" with "Fine | SOU" (= reference case) and with "FOU"), even with a grid that has 22 times less cells than the grid in the reference case (fine grid). A comparison of Figs. 10d,e with Figs. 8d,e shows that larger differences are present at x/L = 0.2 and x/L = 0.5 for V/U M on the coarsest grid. Figure 10f illustrates the large diffusive effect of FOU on the coarsest grid. For FOU, from 0.3 < y/L < 0.9, V/U M has a fairly constant value (≈ −0.55) over the height. For SOU, between 0.65 < y/L < 0.8, V/U M on the fine grid ranges from −0.37 to −0.8 with a strong velocity gradient. The velocity gradient at this location predicted on the coarsest grid with SOU discretization schemes is also much larger than with FOU schemes. Figure 11 depicts profiles of K obtained on the coarsest grid with different discretization schemes. For comparison, the results obtained by SOU on the fine grid are shown as well (= reference case: dotted line). Figure 11a shows the largest differences between FOU and SOU below the source location; i.e. about 25% for 0.1 < y/L < 0.2. Figure 11b shows a difference of more than 100% between FOU and SOU results in the middle of the enclosure (x/L = 0.5, y/L = 0.5) for source location 1. Also at x/L = 0.8 pronounced differences in K are observed, with the largest differences located in the upper part of the enclosure; e.g. up to 100% at y/L = 0.8; however, it must be noted that the concentrations in this region are low. For source location 2, differences up to 100% occur between FOU and SOU along all three lines. The SOU results on the coarsest grid are again closer to the SOU results on the fine grid (= reference case) than to the FOU results on the coarse grid, indicating that the choice of discretization schemes can have a larger effect than the grid resolution. Figure 12 shows the contours of dimensionless mean velocity magnitude (|V|/U inlet;CFD ) (Figs. 12a,d,g,j) and K (Figs. 12b,c,e,f,h,i,k,l) in the vertical center plane obtained with four different simulations: (a-c) FOU on coarsest grid; (d-f) SOU on coarsest grid; (g-i) FOU on fine grid; (j-l) SOU on fine grid. The velocity contours show that the detachment of the wall jet is affected (delayed) by numerical diffusion due to a coarse grid resolution (Fig. 12d vs.  Fig. 12j). Numerical diffusion on the coarsest grid is further enhanced with FOU yielding further delay of jet detachment (Fig. 12a vs. Fig. 12d). In addition, the velocities in the jet after detachment from the top surface clearly depend on the chosen grid resolution and discretization schemes. Figure 12b clearly shows a smeared out distribution of K compared to Fig. 12e and Fig. 12k. Due to numerical diffusion a much larger area with high concentrations (K = 2) is visible. For source location 2 (Figs. 12c,f,i,l), FOU on the coarsest grid leads to lowest concentrations near the inlet wall (left wall), closely followed by FOU on fine grid, due to the lower level of concentration transported along the floor from the source location. As mentioned in Section 4.2, this can be attributed to reduced turbulent mass transport when velocity gradients are smeared out due to numerical diffusion. The higher concentrations along the floor and the inlet wall with SOU (Figs. 12f,l) affect the concentrations in the recirculation cell as well; the highest concentrations are present for SOU on the fine grid (Fig. 12l) and the lowest for FOU on the coarsest grid (Fig. 12c). Comparison of Figs. 12f,i,l again illustrates that the choice for FOU instead of SOU leads to larger deviations than SOU on a much coarser (22 times) grid.

Velocity field
The impact of physical diffusion is examined by changing the TI imposed at the inlet. In the reference case simulation, TI at the inlet (beginning of contraction) was set to 21%. In this section, the results of CFD simulations with TI equal to (1) TI = 1%; (2) TI = 10%; (3) TI = 50%; and (4) TI = 100% are presented together with the results from the reference case (TI = 21%). Figure 13 shows U/U M and V/U M for different TI along the three lines in the vertical center plane. TI has a strong impact on the results in general and on the detachment of the wall jet from the top surface in particular (Fig. 13c). The higher turbulence intensities lead to a delayed detachment of the wall jet, due to the better resilience of the boundary layer to withstand adverse pressure gradients when turbulence levels are higher. The location of maximum velocity ranges from y/L ≈ 0.9 for TI = 100% to y/L ≈ 0.72 for TI = 1% and TI = 10%. The differences in profiles of V/U M are also pronounced; at x/L = 0.5 the maximum value is around V/U M = −0.33 for TI = 1% and 10%, while V/U M is only −0.04 for TI = 100%. The profiles of V/U M at x/L = 0.8 also clearly indicate that an increased TI alters the results and is responsible for a smeared out prediction (decrease of velocity gradients). The effects found (delayed detachment and smaller velocity gradients) are very similar to those found in Figs. 3c,f (lower grid resolutions), in Figs. 8c,f (FOU schemes), and in Figs. 10c,f (lower grid resolution in combination with FOU schemes). Figure 14 shows K along the three vertical lines. Higher TI results in lower K in the majority of the locations studied for source location 1 (Figs. 14a-c). The difference between K for TI = 1% and TI = 21% at x/L = 0.5 can be as large as Note that Section 4.2 showed that numerical diffusion (FOU) actually decreased the concentrations for source location 2; there is a different effect of increased physical diffusion compared to increased numerical diffusion on contaminant transport for this source location. Increasing TI at the inlet not only affects the mean velocity field (smeared out velocity gradients and delayed detachment of the wall jet (see Fig. 13)) and thus convective mass transport, but also indirectly reduces turbulent mass transport due to reduced production of turbulent kinetic energy (smaller velocity gradients). However, on the other hand, the increase of inlet TI, and thus of inlet turbulent kinetic energy, also directly increases the turbulent viscosity ( t = ρC μ f μ k 2 /ε) and thus the turbulent mass fluxes by the increase in turbulent mass diffusivity (D t =  t /Sc t ). A net increase of turbulent mass transport with increasing TI for source location 2 is observed, while Section 4.2 showed that numerical diffusion results in decreased turbulent mass transport. Figure 15 shows the contours of |V|/U inlet;CFD and K in the vertical center plane obtained for (a-c) TI = 1%; (d-f) TI = 21% (= reference case); (g-i) TI = 100%. Inlet TI has a strong impact on the detachment of the wall jet, which occurs further downstream with increasing TI. Inlet TI also impacts the contaminant distribution, with larger concentrations for lower TI for source location 1, while only relatively small differences in concentrations are present for source location 2 (Figs. 15c,f,i).

Discussion
This study focused on the occurrence of numerical and physical diffusion in CFD simulations of mixing ventilation driven by a transitional wall jet in a cubical enclosure. The simulations are validated using reduced-scale experiments by the authors, after which the impact of grid resolution and discretization schemes (numerical diffusion) and inlet turbulence intensity (physical diffusion) is systematically assessed. The study has a few limitations that can give rise to future research on the following topics:  Application for different boundary conditions beyond the transitional flow in the present study with low turbulence intensities. The conclusions might be different when fully turbulent indoor airflow is assessed. For example, the impact of inlet turbulence intensity appeared to be quite significant in the presented study, which is in contradiction with the conclusions from some of the earlier studies (e.g. Saïd et al. 1993;Joubert et al. 1996), who found no strong influence of the inlet turbulence intensity on the results.  Application for other grid types, such as unstructured and non-conformal grids, and different cell types, such as tetrahedral cells. Tetrahedral cells result in an increased level of numerical diffusion compared to structured grids with hexahedral cells.  Assessment of numerical diffusion due to temporal discretization.  Assessment of numerical diffusion for transient simulations such as unsteady RANS (URANS) and Large Eddy Simulation (LES).

Summary and conclusions
Quality assurance in CFD is essential for an accurate and reliable assessment of indoor airflows. Guidelines are available with-among others-recommendations on (I) grid resolution; (II) spatial discretization schemes; (III) inlet turbulence intensity. However, additional research on the impact of the three aforementioned aspects on CFD simulations of low-Reynolds number mixing ventilation flow was needed for the following reasons: I. In general, a systematic grid refinement with a constant factor is recommended to find the grid that provides grid-independent results, or at least to prove that grid convergence is present. The possible errors when a grid-sensitivity study has not been performed with sufficient scrutiny, or when only grid convergence has been obtained, can be important however and need to be assessed. II. To the best knowledge of the authors, the impact of the chosen discretization schemes in combination with the grid resolution employed on the calculated velocities and concentrations for low-Reynolds number ventilation flows has not yet been investigated in detail. III. The effect of the turbulence intensity imposed at the inlet appears to be case/flow dependent and requires additional research, since no consensus is found in literature. Therefore, this paper presented 3D steady RANS CFD simulations with the aim to assess the numerical and physical diffusion in the mean velocity and concentration patterns of transitional (Re ≈ 2,500) mixing ventilation flow in a cubical enclosure. The CFD simulations were validated using reduced-scale experiments by the authors, after which the impact of grid resolution and discretization schemes (numerical diffusion) and inlet turbulence intensity (physical diffusion) was systematically analyzed. The grid-sensitivity analysis showed that the fine grid with 1,899,184 cells provided nearly grid-independent results. In addition, the results show that the vertical velocities were more sensitive to the applied grid resolution than the streamwise velocity and that caution is needed to prevent a possibly erroneous decision when it is solely based on one single parameter, such as the streamwise velocities. The following main conclusions can be drawn: (1) Excessive numerical and physical diffusion leads to erroneous results in terms of delayed detachment of the wall jet and locally decreased velocity gradients.
(2) Excessive numerical diffusion by FOU schemes leads to deviations (up to 100%) in mean velocity and concentration, even on very high-resolution grids. Using a too coarse grid (coarsest grid; 84,669 cells) leads to larger differences between the results obtained with different discretization schemes. However, the velocity gradients in the wall jet predicted on the coarsest grid with SOU are much larger and show a much better resemblance with the measurement data than with FOU. (3) The differences between SOU and FOU on the coarsest grid are larger than the differences between SOU on the coarsest grid and SOU on a 22 times finer grid. Therefore, it is concluded that the choice for a FOU discretization scheme has a larger effect on the results than the chosen grid resolution. (4) Imposing TI values from 1% to 100% at the inlet results in very different velocity fields (enhanced or delayed detachment of wall jet) and different contaminant concentrations (deviations up to 40%). The effects on the velocity field found when increased values of inlet TI (physical diffusion) are imposed (delayed detachment and smaller velocity gradients) are very similar to those found for increased levels of numerical diffusion. (5) The impact of physical diffusion on contaminant transport can markedly differ from that of numerical diffusion. Increased numerical diffusion results in reduced turbulent mass transport due to smeared out velocity gradients and thus net lower turbulent kinetic energy levels due to less production, while increasing the physical diffusion results in a net increase of turbulent kinetic energy and thus an increase of turbulent mass transport.

Recommendations
It is shown that diffusion can have a large impact on indoor airflow patterns and contaminant distribution, therefore, an appropriate choice needs to be made for the grid resolution and discretization schemes (numerical diffusion) and inlet turbulence intensity (physical diffusion) to obtain accurate and trustworthy results. Based on this study the following notes and recommendations are provided:  A grid-sensitivity analysis should be based on the parameter of interest, and possibly for multiple parameters to ensure a proper grid resolution.  FOU discretization schemes lead to too much numerical diffusion, even on very fine grids, and should thus be avoided at all cost.  The difference between FOU and SOU schemes is larger than the difference between two simulations with SOU and different grid resolutions; therefore, enhancing grid quality is better than to resort to FOU when convergence cannot be achieved.  The correct choice for TI at the inlet is important; imposing an erroneous value can lead to very different results with respect to velocities and contaminant concentrations.  Physical and numerical diffusion can cancel each other out to a certain extent for velocities, but the results are not necessarily realistic or trustworthy.