Simulations of Twin Turbulent Planar-Like Jets Injected Into a Large Volume Using RANS

The systems and component designs of advanced reactors can be analyzed through the use of computational fluid dynamics (CFD) codes. In order to justify the usage of these CFD codes and related analysis methods for nuclear reactor design, benchmarking against applicable relevant datasets is needed [1]. There have been instances of previous benchmarking efforts in nuclear thermal hydraulics such as the pandas and mixing tee experiments [2]. The work discussed in this study hopes to build on previous work using the semiclassical problem of the injection of twin planar-like turbulent jets into a large volume (Fig. 1). This includes comparisons to historical insights found in literature and recent experimental efforts at the University of Tennessee, Knoxville [3], and Texas A&M University (TAMU) [4–7]. Utilizing experimental data collected in the twin jet water facility (TJWF), the Star-CCM+ code was benchmarked using the steady Reynolds-averaged Navier–Stokes (RANS) framework. This has not been previously done within all three regions of combining jet behavior with the same amount of resolution at the time of writing.

The first investigation into twin turbulent free shear jets flow physics were by Miller and Comings [8] using hot-wire anemometry (HWA) and static pressure disk probes. Miller and Comings observed the twin jets combine into a single jet due to a subatmospheric pressure region between the two jets with a separating ratio (S/a) of six. In their studies using HWA and static pressure disk probes, Tanaka [9,10] found the subatmospheric pressure region between the jets (with S/a ∼ 8.5 to 25) develops due to the surrounding fluid being entrained. Tanaka observed the combined twin jets compared well to Goerter's analytical velocity profiles derived from the planar jet boundary layer equations [11]. It was noted the turbulent intensity profiles did not have good agreement between a single jet and a combined jet. Tanaka also found that the behavior of the twin jets is independent of Reynolds number (assuming fully developed turbulence before injection). This observation is also supported by a number of later studies by different authors with each using different separation ratios [12–14]. Elbanna [13] investigated twin jets (S/a ∼ 12.5) at a higher Reynolds number and found the three components of turbulent intensity exhibited behavior varying from single jet observations of several authors [15–17] similar to Tanaka. Lin and Sheu found that for a large separating ratio (S/a ∼ 30), the turbulent intensities and Reynolds stress profiles are self-preserving in the combined region [18].

The behavior of twin turbulent free shear jets (Fig. 1) based on literature is summarized for clarity for discussions throughout the study and follows the Anderson vernacular [19]. The flow is characterized by the three regions referred to as the converging, merging, and combining regions. The converging region spans the injection point of the twin jets to the merge point (free stagnation point) of the jets. The region is characterized by the subatmospheric pressure between the two jets that forces the jets together and the time averaged centerline streamwise velocity is primarily negative. The two jets are distinct throughout the merging region. The merging region begins at the merge point which is determined by observing the position where centerline streamwise velocity transitions from negative to positive (MP = U_CL = 0). In merging region, the two jets begin exchanging momentum and the distinct features start to disappear till reaching the combined region. The combined region starts at the combined point where the max centerline streamwise velocity (CP = U_CL,max) and only one jet is observed. In this region, the combined jet will eventually exhibit self-similar behavior to that of a single turbulent planar jet for the velocity profiles. In this region, it is found that there is some disagreement whether the turbulent intensities and Reynolds stresses will be completely similar to a single jet's behavior. This behavior is investigated based on the jet width, a (or hydraulic diameter—D_h) (a is used for this study instead of d used for jet width in the Fig. 1), spacing between the centerline of the two jets (s), and the discharge (initial) velocities (u₀) of the two jets.

At the time of writing, only a few published studies have focused on simulating the behavior of twin turbulent free shear jets using CFD. The earliest effort by Behrouzi and McGuirk was done using RANS (time dependence not stated) of twin jets impinging in very weak to strong crossflow [20]. Using the k–ε turbulence model (unspecified variant) in no crossflow, the authors observed that the simulations predicted the general features of the twin jets (using their own HWA data). The authors noted that turbulent quantities are needed to be compared for a stronger conclusion. A more recent effort by Anderson and Spall involved both experimental and CFD efforts for several separating ratios at a Reynolds number of 6000 using an individual jet width of 2.54 cm [19]. They collected centerline streamwise velocity and centerline Reynolds diagonal (or normal) stress in streamwise and one span wise direction using HWA. They utilized a two-dimensional steady RANS methodology with the k–ε and Reynolds stress model (RSM) turbulence models implemented in fluent v 4.4. The authors found agreement of both turbulence models for the centerline streamwise velocity except for the RSM over-prediction of peak velocities in each case. Both turbulence models predicted faster merging of the jets than in the experiments but found the merge and combine points were in agreement with previous and their experimental data. With regards to the combined point, Anderson and Spall found significant scatter in the combined points from literature and their study. They attributed the scatter to varying upstream conditions among the different experimental test setups. The authors suggested that this is the cause of the slight scatter found in the merge point results as well. Durve et al. simulated the twin and triple jets using the experimental setup of Anderson and Spall (twin jets) and others for the triple jets [21]. In their work, they conclude that the merge point is influenced by the spacing ratio in addition to the jet exit conditions which is similar to Anderson and Spall's conclusion. They also propose the combined point is influenced by the merge point location.

The twin jet water facility was designed to investigate and test hydrodynamic and heat transfer instrumentation (ultrasonic velocimetry, thermocouple rakes, and thermal imaging cameras) for use in opaque fluids for advanced reactors such as sodium fast reactors [22]. It was built by the University of Tennessee, Knoxville [3], for initial testing of the previously mentioned instrumentation and usage in CFD validation efforts. The TJWF shown in Fig. 2 was made to be almost entirely transparent for additional experiments using noninvasive measurement techniques, such as laser Doppler velocimetry (LDV), particle image velocity (PIV), and laser-induced fluorescence. The TJWF was moved to Texas A&M University, where LDV and PIV data have been collected for CFD validation efforts [4–7] using the same boundary conditions as Crosskey and Ruggles. The LDV results of Wang et al. includes streamwise velocity and vorticity profiles, Reynolds stress profiles, and merge/combined points. The LDV experimental data were utilized for the benchmark comparisons found in the Results section. The spatial resolution of the LDV experimental data varied depending on the area of interest. The highest resolution of 2 mm along the horizontal (X-axis) was acquired for X/a between −2.414 and 2.241 and the resolution increases linearly to 8 mm outside of this region in both +X and –X-directions. The streamwise or vertical direction (Z axis) had a spatial resolution of 10 mm for Z/a from 1.72 to 17.24. The resolution decreases to 20 mm for Z/a from 20.09 to 48.28. The reported accuracy of the location of each measurement is 0.01 mm or X/a, Y/a, or Z/a ± 0.00172. The experimental uncertainty was calculated (see Appendix) using the data provided by Wang et al. for the ASME V&V Symposium Benchmark Problem in 2016 [23].

The TJWF as-built dimensions are 762 mm by 1016 mm by 1206.5 mm defined as the outside length, width, and height including the 25.4 mm thick walls. The jet head extends 384.175 mm above the bottom of the TJWF and the top of the tank is 822.325 mm above the jet head. The jet head major dimensions are 171.45 mm by 47.752 mm by 279.4 mm in length, width, and height shown in Fig. 2. Each jet is 87.63 mm by 5.8 mm in length and width and is separated by a 12 mm sheet of acrylic. This yields a centerline to centerline separation distance of 17.8 mm.

In recent years, the twin jet problem has become popular to simulate and benchmark CFD codes and methodologies against. The LDV data in the TJWF are being used for the aforementioned ASME V&V Symposium challenge problem. At the time of writing, only a single CFD analysis using RANS, unsteady Reynolds-averaged Navier–Stokes, and Partial-averaged Navier–Stokes was done of the TJWF experiments [24]. The computational domain only included the tank from the point of injection and utilized the realizable k–ε and k–ω shear stress transport turbulence models. Their results agreed well with the provided LDV and PIV data for profiles in the converging and merging regions of the flow for the velocity and Reynolds stress profiles. The current work discussed builds on this effort by extensively investigating steady RANS in the converging, merging, and combined region of flow and by comparison to analytical jet profiles (Tollmien and Goertler) found in literature [11].

The CFD simulations were done using the steady RANS framework within Star-CCM+ v 10.06.010 and v 11.04.012 [25]. The Numerical Approach section outlines the computational domain and mesh, modeling framework and governing equations, turbulence modeling, and boundary/initial conditions.

The TJWF seen previously in Fig. 2 was simplified to reduce computational effort while maintaining critical geometric features. The simplified geometry is provided in Fig. 3 with a separate view of the jet head for clarity. The major simplifications involve the upstream domain from the inlet slender ducts and the free surface of the water overflowing. The upstream geometry of the inlet slender ducts feeding the inlet jets is assumed insignificant on the prediction of flow behavior. The free surface of the water and its air interface is not directly modeled and the domain is extend approximately 50.8 mm for overflow. The exact boundary treatment is discussed in more detail in the Boundary and Initial Conditions section. These assumptions are justified by the benchmark comparisons in the Results section.

The geometric features (or parameters) of interest based on Fig. 1 are defined as follows. The twin jets have an individual jet width (a) of 5.8 mm and a length (w) of 87.63 mm with a height of 279.4 mm. The aspect ratio (w/a) of the jet is 15.1. The aspect ratio for a duct provides a basic means of determining if secondary flows (spanwise velocity components) are still contributing to the “centerline” behavior of jet. The aspect ratio is not high enough to suggest that the TJWF experiments can be simulated without considering the effect of secondary flows [15,26,27]. This requires three-dimensional simulations to be used to appropriately account for the twin jet mixing behaviors. Based on this information, the Reynolds number of 9100 is defined based off the hydraulic diameter (4 × area/perimeter) of each individual jet (D_h = 0.01088 m). This is consistent with the Reynolds number definition used by the experimental works [3,4,6]. The twin jets are separated (based on the centerline of each jet) with a spacing of 17.8 mm which results in a spacing ratio (S/a) of ∼3.07. This is one of the smallest spacing ratios among the twin jet literature.

The domain is meshed using the standard meshing tool in Star-CCM+ [2] using three built-in meshing models. The mean flow was created using the “Trimmer” model, which is based on a hexahedral mesh that gets “cut off” near geometric features and also known as cut-cell mesh. The near wall regions were meshed using the “prism layer” model that creates prism cells at all near wall regions to resolve the viscous gradients. The “extruder” mesher which extrudes mesh cells from a surface mesh were used to mesh the inlet slender ducts and extend the weir overflows. The domain had specific mesh clustering regions to resolve the highest gradients of the jets free stream behavior. The specific mesh clustering regions are broken up into the jet inlet, the core, and expansion regions of isotropic (same for all directions) refinement.

The mesh refinement was done by refining the base size in the entire domain and the isotropic size in the mesh cluster regions. The base refinement is a bulk refinement commonly seen in CFD works. The mesh refinement in the cluster regions is done due to the largest impact on the flow predictions. This impact is due to those areas having the highest gradients, whereas the near-wall region mesh is not critical for the majority of the domain excluding the inlet jets. The lower importance is attributed to the very low flow velocities from being far from the jets. A summary of the meshing for the majority of the computational domain is present in Table 1.

The resulting mesh with the clustering regions is provided in Fig. 4 for the coarsest mesh for clarity. The jet inlet region (dotted line in Fig. 4) is a rectangular refinement volume that extends above the expected merge point and the length of the jets (into the page—Y-direction). The core region (dash-dot line in Fig. 4) of the expected combined jet is refined in a cylindrical volume that spans from the jet exhausts to the top of the domain. The expansion region (red dotted line) surrounds the jet inlet and core regions in a conical volume to resolve the expected jet spreading and entrainment of ambient fluid. The base size of the overall mesh and the isotropic sizes of the refinement regions are provided in Table 2. The resulting mesh cell count or total degrees-of-freedom (DOF) for each mesh are provided in Table 2 as well.

The extruded regions of the inlet slender ducts and the weir overflows are not changed for each successful mesh refinement. The inlet slender ducts that exhaust the inlet jets need to be appropriately modeled (i.e., a duct is modeled for each jet using the experimental dimensions) to provide the inlet conditions similar to the experiment using the boundary conditions discussed in the Boundary and Initial Conditions section. The inlet slender ducts was created using the as-constructed length of the slender duct (0.2794 m) and 50 “layers” of mesh cells with a stretching ratio of 1.3 in the streamwise direction. The surface mesh regions that the jets were extruded from had minimum and maximum surface sizes of 0.15 and 0.25 mm. The outlets are extruded from the surface where the original computer-aided design model weir overflows end. The outlets are extended a total distance of 600 mm using a 1.4 stretching ratio with 25 layers of mesh. This was done to strongly reduce any backflow/reversed flow issues that would cause numerical instabilities.

where $ui¯$ and $p¯$ are the mean velocity and pressure fields. The Reynolds stress tensor, $ρui′uj′¯$⁠, and the mean strain-rate tensor, $Sij$⁠, are discussed in the Turbulence Modeling section.

The fluid properties of 997.561 kg/m³ for density and 8.8871 × 10⁻⁴ Pa·s for dynamic viscosity of water was used.

The turbulence modeling was approached using two turbulence models, the first being a two-equation standard k–ε with a two-layer wall treatment, and the second being an eight equation elliptic blending Reynolds stress model (EBRSM). These two turbulence models were found to be used in both past twin jet literature and generic CFD analysis.

The equations for the turbulent kinetic energy (k) and dissipation (ε) and details for other modeling constants are found in Rodi's original work [29].

where $Cij$ through $εij$ are the convection, viscous diffusion, turbulent diffusion, redistribution, production, and dissipation terms. The details regarding additional equations (such as the dissipation and elliptic blending equations), these terms, and another model constants can be found in the work by Lardeau and Manceau [33].

The M4 mesh was used with the EBRSM model and the usage is explained in the Results section. It is worth noting that the EBRSM model was not able to reach the same convergence as the standard k–ε turbulence model. The results shown for it have to be taken with an amount of skepticism.

The convergence criteria was based on the residuals of the x, y, and z momentum, continuity, and turbulence parameters. The simulations were run until a reduction of three to four orders of magnitude in the residuals was observed. This was checked by observing the change in parameters and the lack of significant change at points of interest for each simulation. The lack of significant change was quantified as less than 0.5% change from one iteration to the next for each quantity of interest. The observed parameters were the maximum velocity in the entire domain and velocity at points in the converging, merging, and combining regions.

The boundary conditions for the simulations are highlighted and labeled in Fig. 3. The inlet boundary conditions were placed at the physical entrance of each slender ducts that eventually exhaust the jets into the enclosure. The inlet velocity condition that was used for each jet was a uniform velocity profile. A uniform velocity profile was used for each jet since the velocity profile at the beginning of each slender duct was not measured. Further, the slender jets are downstream of stagnation boxes that were not modeled. The entrance length (height of the slender duct simulated) corresponded to approximately 25 hydraulic diameters or 48 jet widths. This is considered sufficient to have fully developed flow for a turbulent regime or at least flow behavior based on the physical geometry of the experimental setup. The uniform velocity was set to 0.75 m/s corresponding to an individual jet Reynolds number of 9100. The inlet turbulent specification is done using the turbulence intensity of 0.053 and a length scale of 0.752 mm for both turbulence models used. The inflow conditions for the EBRSM turbulence model such as the six Reynolds stress components were calculated directly by Star-CCM+ using the previously mentioned turbulent intensity and length scale. The free surface of the water at the top of the domain is defined as a symmetric boundary condition. The symmetry boundary acts as a slip wall condition to reduce the computational requirements with modeling a free surface. The free surface is assumed to experience minor movement and insignificant impact on the results discussed in the Results section. The outlet boundaries were placed on the ends of the overflows and were static pressure outlets. These were set to a static pressure of 0 Pa with a reference domain pressure of atmospheric pressure. The domain was initialized with zero velocity in all cells before the simulations were run.

The results are separated into the mesh sensitivity study, analytical self-similar free shear planar jet comparisons, and benchmarking of the RANS results with the two turbulence models to the available experimental data. The results are based on velocity, vorticity, and Reynolds stress information in characteristic regions or points in the twin jet domain.

where the second component or y-vorticity was used for this study.

where N is the total number of points compared.

The centerline streamwise velocity and associated absolute error is shown in Fig. 5 and Table 4 for each mesh with the standard k–ε turbulence model. It is observed that for even two orders of refinement, the global behavior converges to a consistent behavior and the MAE reduces two order of magnitude for the M2 result. It is worth noting that even the least refined mesh (M1) does capture the behavior but is clearly under resolved up to Z/a ∼ 15. The notable features are the negative velocity region between Z/a ∼ 0 to 3.6 where the converging region exist, the merging region where the maximum velocity of the combined jet is reached, and the combined region where single jet behavior exists.

In order to investigate how well resolved the areas of interest with the highest gradients (close to where the jets exhaust and in the converging region), the streamwise velocity profile and y-vorticity at Z/a ∼ 1.72 are presented in Figs. 6 and 7. The corresponding MAE for each is provided in Table 4. The streamwise velocity provides a clear picture that the M3 mesh is needed as a minimum to resolve the velocity field with a corresponding two order reduction of the MAE. There is minimal visual change in the velocity field at this profile for the meshes above the M3 mesh.

The y-vorticity profiles for increasing levels of refinement show a need for a higher meshing refinement at a minimum of the M4 mesh to properly resolve the gradient term. The MAE does show a two order reduction starting at the M2 mesh but the M4 mesh is clearly needed based on the visual results in Fig. 7. This is not surprising due to the sensitivity of the meshing refinement on the gradient and provides a strong basis for using the M4 mesh for further comparisons. In order to calculate the Reynolds off-diagonal or shear stress profiles for the standard k–ε turbulence model discussed in the Reynolds Stress Profiles at Characteristics Heights section, the gradients of velocity are needed to be properly resolved. The requirement is due to the a-posteri calculation of Reynolds stress using the eddy viscosity term and planar velocity gradients. It is worth noting that for comparing profiles such as velocity, the usage of the M3 mesh would be sufficient.

The single turbulent free shear planar jet velocity profile can be analytically solved using the boundary layer equations which serve as a comparison to the single jet behavior in the TJWF simulations. It provides a verification that the expected behavior is captured and similarly the behavior that should not be captured is not predicted by the CFD simulations. The analytical solutions of Tollmien, derived using mixing length hypothesis, and Goertler, derived using the Prandtl turbulent shear stress relationship, are used for this comparison [11].

The CFD results for the M4 mesh with the standard k–ε turbulence model are compared to the analytical solutions in Fig. 8. The parameter, b, used for the spanwise coordinate in the x-direction is the distance from datum in either direction where the velocity is one-half that of the maximum velocity at the datum. This is done to provide a direct comparison to the analytical solutions of Tollmien and Goertler which used this parameter.

The presented CFD results are in the combined jet region and show an excellent comparison to both analytical solutions at all heights shown. This suggests the CFD results are exhibiting self-similar or preserving behavior expected of the single jet. The CFD results are sufficiently resolving the twin jet physics to properly predict single jet behavior. These results in conjunction with the mesh sensitivity study suggest the M4 mesh prediction is resolved enough to be used for comparison to the experimental results.

The merge and combined points for the LDV and PIV measurements are compared against the standard k–ε and EBRSM turbulence models in Table 5. The reported range for the MP by Wang et al. [4] using LDV is due to the difficulty of obtaining additional data points between Z/a ∼ 1.72 and 3.45. Separately, the range reported for the PIV is attributed to a using two different magnifications of the images acquired with the higher merge point from the higher magnification. The RANS results for both turbulence models predict the MP that is toward the higher end of the experimental measurements. The standard k–ε prediction of the MP is shown to be well bounded by both LDV and PIV measurements. Whereas the EBRSM over predicts by 6.72% and 5.29% difference for the higher end of the LDV and PIV measurements.

The standard k–ε prediction for the CP is higher than the LDV and PIV measurements by 23.62% and 31.62% difference. Similarity, the EBRSM results have a higher over-prediction of 29.30% and 37.23% difference for the LDV and PIV results. The scatter of the predicted CPs and the experimentally observed CPs is expected based on the discussion by Anderson and Spall [19]. They suggest that this scatter is due to the CP being sensitive to the upstream effects while not having a meaningful impact on the combined jet. The higher predictions of the CP and MP by the EBRSM simulation is similar to the higher predictions by Anderson and Spall using a similar RSM [19].

The velocity profiles at different characteristic heights of the twin jet experiments are used for comparison to the simulations of both turbulence models. The comparisons in Fig. 9 are considered to be a representative snapshot of the behavior in the converging, merging, and combined regions. In general, it is observed that the velocity profiles are captured well with a few exceptions. It is worth noting that the error bars for the experimental data are included, which is discussed in detail in the Appendix. It is noted that the error band is not easily observed as a result of being quite small and the provided error bars might not be a full characterization of the experimental uncertainty.

In the converging region (Z/a ∼ 1.72), jet profiles were properly predicted except for a slight experimentally observed asymmetry in the left jet (X/a ∼ −2.3 to 0) not captured by either turbulence model. The velocity profiles of both turbulence models at the height of Z/a ∼ 3.45 or close to the MP overlays almost on top of the experimental results. It is only at the experimental observed height of Z/a ∼ 15.52 or close to the CP that the most notable deviation is found between X/a ∼ −2 to 2 for both turbulence models. The standard k–ε model results show somewhat distinct jets instead of the single combined profile. EBRSM results shows similar behavior while under-predicting the experimental results for X/a ∼ −2 to −4 and X/a ∼ 2 to 4. At the height of Z/a ∼ 48.28, both turbulence models have the expected combined jet profile with some deficiencies. The jet profile of the standard k–ε model between X/a ∼ −2 to 2 is well captured whereas the tails in both directions are slightly under predicted. The jet profile for the EBRSM result is over-predicted between X/a ∼ −2 to 2 and show a larger under-prediction of the tails. The difference in the simulation and experimental results is likely due to the velocity at the top of the combined jet region being quite low. This results in a much longer data collection time for each point collected. The data collection time for these points may not have fully captured the behavior as compared to the results of the higher velocities measured closer to the lower jet regions. The EBRSM result does capture the general behavior, but clear deficiencies are observed which might be due to the meshing deficiencies in this region.

The Reynolds stress (⁠$u′v′¯$⁠) profiles at different characteristic heights of the twin jets are used as a higher order comparison of the experimental and simulation results. The results shown in Fig. 10 are provided as a representative snapshot of the behavior observed and predicted in the converging, merging, and combined regions of the twin jets. In general, the standard k–ε and EBRSM turbulence models are able to capture the shape of the Reynolds stress profiles, but there is deficiencies observed.

In particular, the under prediction of the Reynolds stress is seen in the merging region (Z/a ∼ 1.72) spikes at X/a ∼ −2, −1, 1, and 2. The behavior in this region is likely due to the isotropy assumptions for the Reynolds stress tensor for the standard k–ε turbulence model, whereas the deficiencies in the EBRSM results could be due to the sensitivity of the inlet turbulent boundary conditions on the results. The behavior for the Z/a ∼ 3.45 is captured by the standard k–ε turbulence model while the EBRSM under-predicts the spikes. At the combined point (Z/a ∼ 15.52), the two models capture the experimental data excluding the X/a ∼ −1 to 1 region. The EBRSM model under-predicts the extent of the curves between X/a ∼ −4 to −2 and X/a ∼ 2 to 4. In the combined region (Z/a ∼ 48.28), the under-prediction of the curves for X/a < −2 and X/a > 2 likely corresponds to the same issue seen for the experimental velocity profile related to the data collection time at the same height.

The application of steady RANS modeling using both the standard k–ε and EBRSM turbulence models was done as a benchmark against the twin planar-like turbulent free shear jet experiments in the TJWF. It was found that a sufficient meshing strategy involving both bulk and specific refinement regions yields a sufficiently resolved set of simulations for the standard k–ε turbulence model. This is supported by observing the velocity and vorticity fields to suggest a sufficient refinement level for comparison to analytical solutions (as available) and experimental data. The simulations were compared to the analytical planar jet solutions of Goertler and Tollmien and are viewed to compare favorably. For the experimental data sets, both turbulence models resulted in reasonable predictions of the MPs while the predictions of the CPs were found to be significantly over-predicted. The velocity profile comparison of the two models to the LDV experiment show the behavior is simulated with reasonable success with some deficiencies. Finally, the Reynolds stress profiles of the off-diagonal component of the two turbulence models are observed to capture the trends reasonably well, but the magnitudes in certain regions are not.

The predictions of the converging, merging, and the higher up combined regions (further from the combined point) are shown to compare well against the LDV experimental data. In general, the flow behavior around the combined point is a needed area of improvement for RANS modeling of this approach. Based on the conclusions earlier, this methodology can be used in other similar analysis with reasonable caveats.

Part of this research was conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.³ The authors would like to thank Landon Brockmeyer and Jonathan Lai for their assistance in preparing the data for this paper.

• This material is based upon work supported under a Department of Energy Nuclear Energy University Programs Graduate Fellowship.

a =

jet width (mm)

b =

distance from datum where the velocity is one-half that of the maximum velocity at the datum (mm)

D_h =

hydraulic diameter (m)

E_M#–M5 =

absolute error between M# and M5

M# =

mesh level used for refinement factor

N =

total number of samples or points

Q_M# =

quantity of M# used for absolute error

r₁₂ =

refinement factor

S =

separation distance between jet centerlines (mm)

u₀ =

discharge (initial) velocity of jets

U_CL =

centerline streamwise velocity (ms⁻¹)

w =

jet length (mm)

x =

span-wise coordinate (horizontal component) of measurement plane

y =

span-wise coordinate (horizontal component) of measurement plane

z =

streamwise coordinate (vertical component) of measurement plane

The experimental uncertainty for the velocity in the streamwise and spanwise directions, and off-diagonal Reynolds stress component (⁠$u′v′¯$) on the measurement plane were calculated using data measured at three different points in the domain. These data were provided by Wang and Hassan [23] for the ASME V&V Symposium Benchmark Problem in 2016 for the different participants to construct their own error bands. The data points were provided in locations shown in Fig. 11 with each point in either the converging, merging, or combined region of flow. This allows for an uncertainty representative of each region to be applied to the quantities of interest used for comparison.

The fractional uncertainties (fractional SEM) are calculated by normalizing by the mean value of the quantity of interest. The mean and SEM of each quality are summarized in Table 6 and applied as $x±2αx$⁠.