The complex behavior of thermal fluids in nuclear reactors require the usage of computational fluid dynamics (CFD) codes for design and analysis. In order to use CFD codes, they require regular benchmark problems to ensure the predictions are reasonable representations of reality. The twin jet water facility (TJWF) designed and built at the University of Tennessee, Knoxville was created for this purpose. The facility features twin planar-like turbulent free shear jets injecting fluid into a transparent tank to study a variety of flow behavior. The experimental work using this facility by Texas A&M University was used for the benchmarking activities. This work was conducted using a steady Reynolds-averaged Navier–Stokes formulation to simulate the flow behavior. It was determined that the standard k–ε and elliptic blending Reynolds stress model (EBRSM) turbulence models can be used to simulate the twin jet behavior with reasonable success for design and analysis activities.

## Introduction

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) [47]. 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.

Fig. 1
Fig. 1
Close modal

### Twin Turbulent Free Shear Jets.

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 [1214]. 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 [1517] 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 = UCL = 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 = UCL,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—Dh) (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 (u0) of the two jets.

### Previous Numerical Studies.

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.

## Experimental Test Facility

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 [47] 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].

Fig. 2
Fig. 2
Close modal

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.

### Recent Simulations Efforts of the Twin Jet Water Facility Experiments.

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].

## Numerical Approach

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.

### Computational Domain and Mesh.

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.

Fig. 3
Fig. 3
Close modal

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 (Dh = 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.

Table 1

Meshing parameters for computational domain

Initial base size (mm)Initial maximum cell size (mm)Default growthBoundary growth
2550MediumMedium
Prism layersPrism stretching factorPrism layer total thickness (mm)
41.12
Initial base size (mm)Initial maximum cell size (mm)Default growthBoundary growth
2550MediumMedium
Prism layersPrism stretching factorPrism layer total thickness (mm)
41.12

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.

Fig. 4
Fig. 4
Close modal
Table 2

Mesh parameters for base size and isotropic region refinement

Mesh identifierBase size (mm)Maximum cell size (mm)Jet inlet size (mm)Core size (mm)Expansion size (mm)Cell count (DOF)
M12550510204.9 × 106
M212.5252.551014.9 × 106
M36.2512.51.252.5519.2 × 106
M43.1256.250.6251.252.553.2 × 106
M51.56253.130.31250.6251.25322 × 106
Mesh identifierBase size (mm)Maximum cell size (mm)Jet inlet size (mm)Core size (mm)Expansion size (mm)Cell count (DOF)
M12550510204.9 × 106
M212.5252.551014.9 × 106
M36.2512.51.252.5519.2 × 106
M43.1256.250.6251.252.553.2 × 106
M51.56253.130.31250.6251.25322 × 106

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.

### Modeling Framework.

Star-CCM+ is a finite volume-based solver using the SIMPLE algorithm to couple the velocity and pressure fields. The steady RANS framework for isothermal flow (i.e., no energy equation solved) was selected for this work. The steady single-phase RANS equations for continuity and momentum used are as follows:
$∂ui¯∂xi=0$
(1)
and
$ρuj¯∂ui¯∂xj=−∂p¯∂xi+∂∂xj2μSij−ρui′uj′¯$
(2)

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/m3 for density and 8.8871 × 10−4 Pa·s for dynamic viscosity of water was used.

#### Turbulence Modeling.

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 standard kε turbulence model of Launder and Spalding [28] with the Rodi's [29] approach of a two-layer wall treatment is used. It is selected for its widespread use in industrial and research applications [3032]. The standard kε turbulence model was used for all meshes in the mesh convergence study. The Reynolds stress tensor is denoted as $ρui′uj′¯$, where $ui′$ is the fluctuating velocity derived from the Reynolds decomposition of the velocity signal. The Reynolds stress tensor is related to the mean strain-rate tensor, $Sij$, though the use of the Boussinesq eddy viscosity assumption in the standard kε turbulence model. This relationship is represented as follows:
$−ui′uj′¯=2ντSij−23kδij$
(3)
where k is the turbulent kinetic energy, $δij$ is the Kronecker delta, and the eddy viscosity, $νt$ is calculated using
$ντ=cμk2ε$
(4)

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

The EBRSM turbulence model is based on the variant of Lardeau and Manceau [33] with the two-layer-wall treatment. This EBRSM model was selected due to previous usage of Reynolds stress model in twin jet flows [19] based on the EBRSM model's ability to resolve anisotropy of the flow. The EBRSM turbulence model directly simulates the Reynolds stress terms through the Reynolds stress transport equations
$−∂ui′uj′¯∂t−Cij+Dijv+DijT+ϕij*+Pij−εij=0$
(5)

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.

#### Numerical Convergence Criteria.

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.

### Boundary and Initial Conditions.

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.

## Results

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.

### Mesh Sensitivity.

The mesh sensitivity was determined based on the centerline velocity profile between the two jets and the velocity and vorticity profiles in the area of highest gradients of the flow. The centerline velocity profile provides a global observation of the meshing impact and sensitivity. Whereas, the velocity and vorticity profiles provides a local observation of the meshing impact with direct dependence on the mesh itself. The vorticity field is calculated by taking the curl of velocity which is shown using
$ω=∇×u=∂uz∂y−∂uy∂z,∂ux∂z−∂uz∂x,∂uy∂x−∂ux∂y$
(6)

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

These quantities are chosen due to the relevance to the physics of interest and the significant impact of meshing resolution. In order to suggest sufficient refinement between meshes, a metric is required. The refinement factor based on the total number of mesh cells of the finer to coarser meshes is used
$r21=N2N113$
(7)
where N1 is the coarser mesh and N2 is the finer mesh. A refinement factor of at least 1.3 between meshes is recommend by the ASME Journal of Fluid Engineering [34]. The refinement factor between the successive meshes used in this analysis are shown in Table 3. Three out of the four refinement factors are above or significantly above 1.3 with only one refinement level below at 1.09. In order to quantify the differences between each mesh and the finest mesh, the absolute error was used
$EM#−M5=QM#−QM5$
(8)
where QM# is the quantity of interest for M# compared to QM5 which is the quantity of interest for M5. For the sake of comparison, what would be considered the “approximate” value of reality is the M5 quantity of interest. Additionally, the mean absolute error (MAE) is used
$MAEM#−M5=1N∑i=1NQM#,i−QM5,i$
(9)

where N is the total number of points compared.

Table 3

Refinement factors for meshes used

Mesh identifierCell count (DOF)Refinement factor
M14.9 × 106n/a
M214.9 × 1061.45
M319.2 × 1061.09
M453.2 × 1061.40
M5322 × 1061.82
Mesh identifierCell count (DOF)Refinement factor
M14.9 × 106n/a
M214.9 × 1061.45
M319.2 × 1061.09
M453.2 × 1061.40
M5322 × 1061.82

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.

Fig. 5
Fig. 5
Close modal
Table 4

Mean absolute error for mesh sensitivity comparisons of centerline velocity, velocity profiles, and vorticity profiles

UCLUω
MAEM1–M50.01870.05500.1932
MAEM2–M50.00730.01970.0847
MAEM3–M50.00250.00670.0298
MAEM4–M50.00130.00530.0130
UCLUω
MAEM1–M50.01870.05500.1932
MAEM2–M50.00730.01970.0847
MAEM3–M50.00250.00670.0298
MAEM4–M50.00130.00530.0130

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.

Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal

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.

### Analytical Self-Similar Turbulent Free Shear Planar Jet Comparisons.

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.

Fig. 8
Fig. 8
Close modal

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.

### Benchmarking of Reynolds-Averaged Navier–Stokes Simulations

#### Merge and Combined Points.

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.

Table 5

Merge and combined points for RANS comparison to experimental benchmark

MPCP
LDV (TAMU)1.72–3.4516.84
PIV (TAMU)2.66–3.515.52
Standard kε3.3021.35
EBRSM3.6922.62
MPCP
LDV (TAMU)1.72–3.4516.84
PIV (TAMU)2.66–3.515.52
Standard kε3.3021.35
EBRSM3.6922.62

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].

#### Velocity Profiles and Characteristics Heights.

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.

Fig. 9
Fig. 9
Close modal

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.

#### ReynoldsStressu′v′¯ProfilesatCharacteristicsHeights.

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.

Fig. 10
Fig. 10
Close modal

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.

## Conclusions

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.

## Acknowledgment

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

## Funding Data

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

## Nomenclature

a =

jet width (mm)

b =

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

Dh =

hydraulic diameter (m)

EM#M5 =

absolute error between M# and M5

M# =

mesh level used for refinement factor

N =

total number of samples or points

QM# =

quantity of M# used for absolute error

r12 =

refinement factor

S =

separation distance between jet centerlines (mm)

u0 =

discharge (initial) velocity of jets

UCL =

centerline streamwise velocity (ms−1)

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

### Appendix

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.

Fig. 11
Fig. 11
Close modal
Each of the points shown were sampled 3000 times randomly in time for the stream-wise and span-wise velocity signals with an example shown in Fig. 12. The mean, standard deviation, and standard error of the mean (SEM) of each velocity signal (both stream-wise and span-wise components) were calculated using the following equations:
$μu¯=∑i=1NuiN$
(A1)
$σu¯=1N−1∑i=1Nui−μu¯2$
(A2)
$αu¯=σu¯N$
(A3)
Fig. 12
Fig. 12
Close modal
where N is the number of samples and u is the signal. The Reynolds stress off-diagonal component used for comparison is calculated by taking the covariance of the two signal components. This is calculated using the following equation:
$u′v′¯=covui,vi=Euivi−EuiEvi=∑i=1NuiviN−μu¯μv¯$
(A4)
The SEM for the $u′v′¯$ was calculated by first calculating the cross-correlation of the stream-wise and span-wise components
$ρuv=covui,viσu¯σv¯$
(A5)
Then using above, the SEM was calculated using the following equation found in Ref. [35]:
$αu′v′¯=σu¯σv¯1+ρuv2N−1$
(A6)

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$.

Table 6

Mean values and fractional SEM at points in the converging, merging, and combined regions

$X/a$$Z/a$$μu¯m/s$$αu¯/μu¯$$μv¯m/s$$αv¯/μv¯$$u′v′¯m/s$2$αu′v′¯/u′v′$
1.553.450.7510.00287−0.03810.04160.003430.0575
−1.5515.520.5090.003630.00062.81−0.00360.0513
0.5234.480.4950.003160.01970.07160.001480.0836
$X/a$$Z/a$$μu¯m/s$$αu¯/μu¯$$μv¯m/s$$αv¯/μv¯$$u′v′¯m/s$2$αu′v′¯/u′v′$
1.553.450.7510.00287−0.03810.04160.003430.0575
−1.5515.520.5090.003630.00062.81−0.00360.0513
0.5234.480.4950.003160.01970.07160.001480.0836

## References

1.
Harvego
,
E. A.
,
Schultz
,
R. R.
, and
Crane
,
R. L.
,
2011
, “
Development of a Consensus Standard for Verification and Validation of Nuclear System Thermal-Fluids Software
,”
Nucl. Eng. Des.
,
241
(
12
), pp.
4691
4696
.
2.
Obabko
,
A.
,
Fischer
,
P.
,
Tautges
,
T.
,
Karabasov
,
S.
,
Goloviznin
,
V.
, and
Zaytsev
,
M.
,
2011
, “
CFD Validation in OECD/NEA T—Junction Benchmark
,” Argonne National Laboratory, Argonne, IL, Report No.
ANL/NE-11/25
.
3.
Crosskey
,
M.
, and
Ruggles
,
A.
,
2014
, “
UTK Twin Jet Water Facility Computational Fluid Dynamics Validation Data Set
,”
International Congress on Advances in Nuclear Power Plants (ICAPP '14)
, Charlotte, NC, Apr. 6–9, pp.
1940
1945
.
4.
Wang
,
H.
,
Lee
,
S.
,
Hassan
,
Y. A.
, and
Ruggles
,
A. E.
,
2016
, “
Laser-Doppler Measurements of the Turbulent Mixing of Two Rectangular Water Jets Impinging on a Stationary Pool
,”
Int. J. Heat Mass Transfer
,
92
, pp.
206
227
.
5.
Wang
,
H.
,
Lee
,
S.
,
Hassan
,
Y. A.
, and
Ruggles
,
A. E.
,
2018
, “
Corrigendum to ‘Laser-Doppler Measurements of the Turbulent Mixing of Two Rectangular Water Jets Impinging on a Stationary Pool’ [Int. J. Heat Mass Transfer 92 (2016) 206–227
,”
Int. J. Heat Mass Transfer
,
118
, p. 1373.
6.
Wang
,
H.
,
Lee
,
S.
, and
Hassan
,
Y. A.
,
2015
, “
Particle Image Velocimetry Measurements of the Flow in the Converging Region of Two Parallel Jets
,”
Nucl. Eng. Des.
,
306
, pp.
89
97
.
7.
Wang
,
H.
,
Lee
,
S.
, and
Hassan
,
Y. A.
,
2017
, “
Corrigendum to Particle Image Velocimetry Measurements of the Flow in the Converging Region of Two Parallel Jets [Nuclear Engineering and Design 306, (2016) 89–97]
,”
Nucl. Eng. Des.
,
324
, p.
416
.
8.
Miller
,
D. R.
, and
Comings
,
E. W.
,
1960
, “
Force-Momentum Fields in a Dual-Jet Flow
,”
J. Fluid Mech.
,
7
(
2
), pp.
237
256
.
9.
Tanaka
,
E.
,
1970
, “
The Interference of Two-Dimensional Parallel Jets (1st Report, Experiments on Dual Jets)
,”
Bull. JSME
,
13
(
56
), pp.
272
280
.
10.
Tanaka
,
E.
,
1974
, “
The Interference of Two-Dimensional Parallel Jets (2nd Report, Experiments on the Combined Flow of Dual Jets)
,”
Bull. JSME
,
17
(
109
), pp.
920
927
.
11.
Rajaratnam
,
N.
,
1976
,
Turblent Jets
,
Elsevier
, Amsterdam, The Netherlands.
12.
Marsters
,
G.
,
1977
, “
The Interaction of Two Plane Parallel Jets
,”
AIAA J.
,
15
(
12
), pp.
1756
1762
.
13.
Elbanna
,
H.
,
Gahin
,
S.
, and
Rashed
,
M. I. I.
,
1982
, “
Investigation of Two Plane Parallel Jets
,”
AIAA J.
,
21
(
7
), pp.
986
991
.
14.
Nasr
,
A.
, and
Lai
,
J. C. S.
,
1997
, “
Two Parallel Plane Jets: Mean Flow and Effects of Acoustic Excitation
,”
Exp. Fluids
,
22
(
3
), pp.
251
260
.
15.
Gutmark
,
E.
, and
Wygnanski
,
I.
,
1976
, “
The Planar Turbulent Jet
,”
J. Fluid Mech.
,
73
(
3
), pp.
465
495
.
16.
,
L. J. S.
,
1965
, “
The Structure of a Self-Preserving Turbulent Plane Jet
,”
J. Fluid Mech.
,
23
(
1
), pp.
31
64
.
17.
,
G.
,
1965
, “
Hot-Wire Measurements in a Plane Turbulent Jet
,”
ASME J. Appl. Mech
,
32
(
4
), pp.
721
734
.
18.
Lin
,
Y. F.
, and
Sheu
,
M. J.
,
1990
, “
Investigation of Two Plane Parallel Unventilated Jets
,”
Exp. Fluids
,
10
(
1
), pp.
17
22
.
19.
Anderson
,
E. A.
, and
Spall
,
R. E.
,
2001
, “
Experimental and Numerical Investigation of Two-Dimensional Parallel Jets
,”
ASME J. Fluids Eng.
,
123
(
2
), pp.
401
406
.
20.
Behrouzi
,
P.
, and
McGuirk
,
J. J.
,
1998
, “
Laser Doppler Velocimetry Measurements of Twin-Jet Impingement Flow for Validation of Computational Models
,”
Opt. Lasers Eng.
,
30
(
3–4
), pp.
265
277
.
21.
Durve
,
A.
,
Patwardhan
,
A. W.
,
Banarjee
,
I.
,
,
G.
, and
Vaidyanathan
,
G.
,
2012
, “
Numerical Investigation of Mixing in Parallel Jets
,”
Nucl. Eng. Des.
,
242
, pp.
78
90
.
22.
Tokuhiro
,
A.
, and
Omotowa
,
O.
,
2012
,
Data Collection Methods for Validation of Advanced Multi-Resolution Fast Reactor Simulations
, University of Idaho,
Idaho Falls, ID
.
23.
Wang
,
H.
, and
Hassan
,
Y. A.
,
2015
,
Benchmark Data for Validation of Computational Simulations of Nuclear System Thermal Fluids Behavior
, Texas A&M University,
College Station, TX
.
24.
Li
,
H.
,
Anand
,
N.
, and
Hassan
,
Y. A.
,
2018
, “
Computational Study of Turbulent Flow Interaction Between Twin Rectangular Jets
,”
Int. J. Heat Mass Transfer
,
119
, pp.
752
767
.
25.
2016
, “
Star-CCM+ 10.06.010 User's Guide
26.
Vinuesa
,
R.
,
Noorani
,
A.
,
Lozano-Durán
,
A.
,
Khoury
,
G. E.
,
Schlatter
,
P.
,
Fischer
,
P. F.
, and
Nagib
,
H. M.
,
2014
, “
Aspect Ratio Effects in Turbulent Duct Flows Studied Through Direct Numerical Simulation
,”
J. Turbul.
,
15
(
10
), pp.
677
706
.
27.
Vinuesa
,
R.
,
Schlatter
,
P.
, and
Nagib
,
H. M.
,
2015
, “
Characterization of the Secondary Flow in Turbulent Rectangular Ducts With Varying Aspect Ratio
,”
International Symposium on Turbulence and Shear Flow Phenomena
, Melbourne, Australia, June 30–July 3, pp.
1
6
.
28.
Launder
,
B. E.
, and
Spalding
,
D. B.
,
1972
,
Mathematical Models of Turbulence
,
,
London
.
29.
Rodi
,
W.
,
1991
, “
Experience With Two-Layer Models Combining the K-Epsilon Model With a One-Equation Model Near the Wall
,”
AIAA
Paper No. 91-0216.
30.
Iaccarino
,
G.
,
2001
, “
Predictions of a Turbulent Separated Flow Using Commercial CFD Codes
,”
ASME J. Fluids Eng.
,
123
(
4
), pp.
819
828
.
31.
Fu
,
C.
,
Uddin
,
M.
, and
Curley
,
A.
,
2016
, “
Insights Derived From CFD Studies on the Evolution of Planar Wall Jets
,”
Eng. Appl. Comput. Fluid Mech.
,
10
(
1
), pp.
44
56
.
32.
Kimouche
,
A.
,
Mataoui
,
A.
,
Oztop
,
H. F.
, and
Abu-Hamdeh
,
N.
,
2017
, “
Analysis of Heat Transfer of Different Nanofluids Flow Through an Abrupt Expansion Pipe
,”
Appl. Therm. Eng.
,
112
, pp.
965
974
.
33.
Lardeau
,
S.
, and
Manceau
,
R.
,
2014
, “
Computations of Complex Flow Configurations Using a Modified Elliptic-Blending Reynolds-Stress Model
,”
Tenth ERCOFTAC International Symposium on Engineering Turbulence Modelling and Measurements
, Marbella, Spain, Sept. 17–19, pp. 1–8.
34.
Celik
,
I. B.
,
Ghia
,
U.
,
Roache
,
P. J.
,
Freitas
,
C. J.
,
Coleman
,
H.
, and
,
P. E.
,
2008
, “
Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications
,”
ASME J. Fluids Eng.
,
130
(
7
), p. 078001.
35.
Sciacchitano
,
A.
, and
Wieneke
,
B.
,
2016
, “
PIV Uncertainty Propagation
,”
Meas. Sci. Technol.
,
27
(
8
), p.
84006
.