Large Eddy Simulation of Thermo-Hydraulic Mixing in a T-Junction

Unsteady heat transfer problems that are associated with non-isothermal flow mixing in pipe flows have long been the topic of concern in the nuclear engineering community because of the relation to thermal fatigue of nuclear power plant pipe systems. When turbulent flow streams of different velocity and density rapidly mix at the right angle, a contact interface between the two mixing streams oscillates and breaks down because of hydrodynamic instabilities, and large-scale unsteady flow structures emerge (Figure 1). These structures lead to low-frequency oscillations at the scale of the pipe diameter, D, with a period scaling as O(D/U), where U is the characteristic flow velocity. If the mixing flow streams are of different temperatures, the hydrodynamic oscillations are accompanied by thermal fluctuations (thermal striping) on the pipe wall. The latter accelerate thermal-mechanical fatigue, damage the pipe structure and, ultimately, cause its failure.


Introduction
Unsteady heat transfer problems that are associated with non-isothermal flow mixing in pipe flows have long been the topic of concern in the nuclear engineering community because of the relation to thermal fatigue of nuclear power plant pipe systems.When turbulent flow streams of different velocity and density rapidly mix at the right angle, a contact interface between the two mixing streams oscillates and breaks down because of hydrodynamic instabilities, and large-scale unsteady flow structures emerge (Figure 1).These structures lead to low-frequency oscillations at the scale of the pipe diameter, D, with a period scaling as O(D/U), where U is the characteristic flow velocity.If the mixing flow streams are of different temperatures, the hydrodynamic oscillations are accompanied by thermal fluctuations (thermal striping) on the pipe wall.The latter accelerate thermal-mechanical fatigue, damage the pipe structure and, ultimately, cause its failure.
The importance of this phenomenon prompted the nuclear energy modeling and simulation community to establish a common benchmark to test the ability of computational fluid dynamics (CFD) tools to predict the effect of thermal striping.The benchmark is based on the thermal and velocity data measurements obtained in a recent Vattenfall experiment designed specifically for this purpose [1,2].Because thermal striping is associated with large-scale anisotropic mixing, standard engineering modeling tools based on time-averaging approaches, such as Reynolds-averaged Navier-Stokes (RANS) models, are badly suited.Consequently, one must consider unsteady modeling methods such as unsteady RANS (URANS), parabolised stability equations (PSE), or large eddy simulation (LES).Among these choices, the LES approach is most generic and modeling-assumption free, unlike the PSE and URANS methods, which, for instance, assume the existence of a spectral gap between the large and the small scales.On the other hand, owing to recent advances in computing, LES methods are becoming an increasingly affordable tool.In this chapter, we compare the results of three LES codes produced for the Vattenfall benchmark problem: Nek5000, developed at Argonne National Laboratory in the United States, and CABARET and Conv3D, developed at the Moscow Institute for Nuclear Energy Safety (IBRAE) in Russia.Nek5000 is an open source code based on the spectral-element method (SEM), a high-order weighted residual technique that combines the geometric flexibility of the finite-element method with the tensor-product efficiencies of spectral methods [3,4].CABARET, which stands for Compact Accurately Boundary Adjusting high REsolution Technique, is based on the scheme of [5].CABARET was developed as a significant upgrade of the second-order upwind leapfrog scheme for linear advection equation [6,7] to nonlinear conservation laws in one and two dimensions [8][9][10].The CABARET scheme is second-order accurate on nonuniform grids in space and time, is nondissipative, and is low-dispersive.For nonlinear flow calculations, it is equipped with a conservative nonlinear flux correction that is directly based on the maximum principle.For 3D calculations, a new unstructured CABARET solver was developed (e.g.[11,12]) that operates mainly with hexagonal meshes.Conv3D is a well-established solver of IBRAE that has been validated on various experimental and benchmark data [13,14].For ease of the grid generation in the case of complex geometries, Conv3D utilizes the immersed boundary method (IBM) on Cartesian meshes with a cut-cell approach.The underlying scheme of Conv3D [15] is a low-dissipative, second-order approximation in space and a first-order approximation in time.
This chapter is organized as follows.In Section 2 we briefly describe the experimental setup and data collection procedure.In Sections 3-5 we outline the computational models.The results are discussed in Section 6, including a comparison of data from the simulations and experiment.Concluding remarks and a brief discussion of future work are provided in Section 7.

Experimental configuration
The Vatenfall experiment [16] is based on water flow in a main pipe of diameter D=140 mm with a side branch of diameter D H =100 mm adjoining the main at a 90 degree angle.The pipes and the T-junction, which is made from a plexiglass block, are transparent so that the velocity can be measured with laser Doppler anemometry (LDA).Velocity data was taken under isothermal conditions with both flows entering at 19 • C. In order to measure thermal striping, time-dependent temperature data was collected from thermocouples downstream of the T-junction with flow at 19  The flow enters the cold (main) branch from a stagnation chamber located 80 D upstream of the T-junction and is assumed to be a fully developed turbulent flow by the time it reaches the T-junction.The hot branch flow enters at 20D H upstream and is not quite fully developed as it enters the T-junction.The inlet flow rates are 9 and 6 liters per second (l/s), respectively in the cold and hot branches, which corresponds to a Reynolds number of Re ≈ 80, 000 and 100, 000, respectively.The key dimensional parameters are summarized in Table 1 and their nondimensional counterparts in Table 2.More detailed description of the T-junction benchmark experiment can be found in [1,2].

Nek5000 simulations
The Nek5000 simulations are based on the spectral element method developed by Patera [3].Nek5000 supports two formulations for spatial and temporal discretization of the Navier-Stokes equations.The first is the lP N − lP N−2 method [17][18][19], in which the velocity/pressure spaces are based on tensor-product polynomials of degree N and N − 2, respectively, in the reference element Ω := [−1, 1] d , for d = 2 or 3.The computational domain consists of the union of E elements Ω e , each of which is parametrically mapped from Ω to yield a body-fitted mesh.The second is the low-Mach number formulation due to Tomboulides and Orszag [20,21], which uses consistent order-N approximation spaces for both the velocity and pressure.The low-Mach number formulation is also valid in the zero-Mach (incompressible) limit [22].Both formulations yield a decoupled set of elliptic problems to be solved at each timestep.In d=3 space dimensions, one has three Helmholtz solves of the form (βI d, and a pressure Poisson solve of the form Ap n = g n at each timestep, t n , n = 1, . . . .Here, A is the symmetric positive-definite Laplace operator, and β is an order-unity coefficient coming from a third-order backward-difference approximation to the temporal derivative.(For the lP N − lP N−2 method, the Laplace operator A is replaced by a spectrally equivalent matrix arising from the unsteady Stokes equations [4,19].)For marginally resolved LES cases, we find that the higher-order pressure approximation of the lP N − lP N formulation tends to yield improved skin-friction estimates, and this is consequently the formulation considered here.Consisting of E = 62, 176 elements, the computational mesh for the Nek5000 simulations was generated by using CUBIT and read through the Nek-MOAB coupling interface [23].Within each element, velocity and pressure are represented as Lagrange interpolating polynomials on tensor products of Nth-order Gauss-Lobatto-Legendre points.Unless otherwise noted, all the simulations were run with polynomial order N=7, corresponding to a total number of mesh points n ≈ EN 3 ≈ 21 million.Figure 2 shows a closeup of the mesh in the vicinity of the origin, which is located at the intersection of the branch centerlines.The inlet for the main branch is at x=-9.2143 and for the hot branch at z=6.4286.These lengths permitted generation of fully developed turbulence upstream of the T-junction, as described below.The outlet at x=22 was chosen to allow downstream tracking of temperature data at locations provided in the experiment.Away from the origin, the axial extent of the spectral elements in the main branch is 0.18 D, corresponding to a maximum axial mesh spacing of δx max = 0.0377.At the wall, the wall-normal element size is 0.01222, corresponding to a minimum spacing of δy n ≈ 0.0008.The submitted simulations were run with Re = 40, 000 in the inlet branches, yielding Re = 60, 000 in the outlet.Downstream of the T-junction, the first grid point away from the wall is thus at y + ≈ 2.5 in wall units.
Inlet flow conditions in the main branch are based on a recycling technique in which the inlet velocity at time t n is given by αu( x, y, z, t n−1 ), where x = −3 and α is chosen to ensure that the mass flow rate at inflow is constant ( u n (−9, y, z)dy dz ≡ π/4).Recycling is also used for the hot branch, save that the inflow condition is 0.8βu(x, y, z, t n−1 ) − 0.2U H , with β chosen so that the average inflow velocity is -U H and z = 2.1.The 0.8 multiplier was added in order to give a flatter profile characteristic of the non-fully-developed flow realized in the experiment, but a systematic study of this parameter choice was not performed (see also Section 6.1.1).
Initial conditions for all branches were taken from fully developed turbulent pipe flow simulations at Re = 40, 000.The timestep size was ∆t = 2.5 × 10 −4 in convective time units, or about 6 × 10 −5 seconds in physical units corresponding to the experiment.The simulation was run to a time of t=28 convective time units prior to acquiring data.Data was then collected over the interval t ∈ [28,58] (in convective time units) for the benchmark submission 1 and over t ∈ [28,104] subsequently (longer average) both with N = 7.In addition to N = 7 results (i.e., with n ≈ EN 3 ≈ 2.1 × 10 7 points), Nek5000 runs were conducted with N = 5 (n ≈ 7.7 × 10 6 points) and with N = 9 (n ≈ 4.5 × 10 7 points).The case with N = 5 started with the N = 7 results and was run and averaged over about 110 convective time units (i.e., about 26 seconds).The timestep size for N = 5 runs was twice as big (i.e,∆t = 5 × 10 −4 ).The benchmark submission results required approximately 3.4 × 10 5 CPU hours on Intrepid (IBM BlueGene/P) while the follow-up study for a longer time average results took additional 5.9 × 10 5 CPU hours.Note that all Nek5000 results reported here were obtained with constant density equal to the nondimensional value of 1.000 (see Table 2).

CABARET simulations
The system of Navier-Stokes equations in a slightly compressible form (i.e., for a constant sound speed) is solved with a new code based on CABARET.In order to lessen the computational grid requirements, the code uses hybrid unstructured hexahedral/tetrahedral grids.CABARET is an extension of the original second-order, upwind leapfrog scheme [6] to nonlinear conservation laws [5,24] and to multiple dimensions [9,25].In summary, CABARET is an explicit, nondissipative, conservative finite-difference scheme of second-order approximation in space and time.In addition to having low numerical dispersion, CABARET has a very compact computational stencil because of the use of separate variables for fluxes and conservation.The stencil is staggered in space and time and for advection includes only one computational cell.For nonlinear flows, CABARET uses a low-dissipative, conservative flux correction method directly based on the maximum principle [8].In the LES framework, the nonlinear flux correction plays the role of implicit turbulence closure following the MILES approach of Grinstein and Fureby [26, ch.5], which was discussed in the ocean modeling context in [10].
A detailed description of the CABARET code on a mixed unstructured grid will be the subject of a future publication; an outline of the method on a structured Cartesian grid is given below.Let us consider the governing equations written in the standard conservation form: where the sources in the right-hand side include viscous terms.By mapping the physical domain to the grid space-time coordinates (x, y, z, t) → (i, j, k, n) and referring control volumes to the cell centers (fractional indices) and fluxes to the cell faces (integer indices), the algorithm proceeds from the known solution at time level (n) to the next timestep (n+1) as follows: • Conservation predictor step: • Upwind extrapolation based on the characteristic decomposition: • For each cell face, decompose the conservation and flux variables into local Riemann fields, U → R q , q = 1, . . ., N, that correspond to the local cell-face-normal coordinate basis, where N is the dimension of the system.
• For each cell face at the new timestep, compute a dual set of preliminary local Riemann variables that correspond to the upwind and downwind extrapolation of the characteristic fields, for example, Rn+1 for the upwind/downwind conservation volumes.
• Correct the characteristic flux fields if they lie outside the monotonicity bounds For reconstructing a single set of flux variables at the cell face, use an approximate Riemann solver.
• Conservation corrector step: where a second-order central approximation is used for the viscous term.
For the T-junction problem with the CABARET method, two hybrid computational grids are used with 0.53 and 4 million cells.Figure 3 shows the layout of (a) the computational domain, (b,d) the hexahedral grid with a uniform Cartesian block in the pipe center, and (c,e) a small collar area of the pipe junction covered by the mixed hexahedral/tetrahedral elements for the smaller and larger grid, respectively.For specifying inlet boundary conditions, a recycling technique is used similar to that for Nek5000.The outflow boundary is prescribed by using characteristic boundary conditions.For the inlet boundaries at the main and the hot branch, laminar inflow conditions are specified based on prescribing the mean flow velocity profiles.The length of the pipe upstream of the junction was sufficiently far from the junction to permit an adequate turbulent flow upstream of the junction.The outflow boundary is imposed at 20 jet diameters downstream of the junction, where characteristic boundary conditions are set.
In order to speed statistical convergence, the LES CABARET solution was started from a precursor RANS k-epsilon calculation.The CABARET simulation was then run for 10 seconds to allow the statistics to settle down.This was followed by the production run during which the required solution fields were stored for a duration of 5 and 10 seconds.The computations on grids with 0.53 and 4 million cells required approximately 1.8 × 10 4 and 2.9 × 10 5 CPU hours on Lomonosov (Intel Xeon X5570 / X5670 processors).

Conv3D simulations
Researchers at IBRAE have been developing a 3D, unified, numerical thermal-hydraulic technique for safety analysis of the nuclear power plants, which includes (1) methods, algorithms, and software for automatically generating computing grids with local refinement near body borders; (2) methods and algorithms for solving heat and mass transfer compressible/incompressible problems for research of 3D thermal-hydraulic phenomena; and (3) approaches for modeling turbulence.
For grid construction at IBRAE an automatic technology using CAD systems for designing the computational domain has been developed.A generation of structured orthogonal/Cartesian grids with a local refinement near boundaries is incorporated into a specially developed program that has a user-friendly interface and can be utilized on parallel computers [27].The computational technique is based on the developed algorithms with small-scheme diffusion, for which discrete approximations are constructed with the use of finite-volume methods and fully staggered grids.For modeling 3D turbulent single-phase flows, the LES approach (commutative filters) and a quasi-direct numerical simulation (QDNS) approach are used.For simulation of 3D turbulent two-phase flows by means of DNS, detailed grids and effective numerical methods developed at IBRAE for solving CFD problems are applied.For observing the interface of two-phase flow, a modified level set (LS) method and multidimensional advection/convection schemas of total variation diminishing (TVD) type with small scheme diffusion involving subgrid simulation (with local resolution) are used.The Conv3D code is fully parallelized and highly effective on high-performance computers.The developed modules were validated on a series of well-known tests with Rayleigh numbers ranging from 10 6 to 10 16 and Reynolds numbers ranging from 10 3 to 10 5 .
In order to simulate thermal-hydraulic phenomena in incompressible media, the time-dependent incompressible Navier-Stokes equations in the primitive variables [13] coupled with the energy equation are used: where p is pressure, normalized by the density.The basic features of the numerical algorithm [14,28] are the following.An operator-splitting scheme for the Navier-Stokes equations is used as the predictor-corrector procedure with correction for the pressure δp: In order to construct the time-integration scheme for the energy equation, its operators are decomposed into two parts associated with the enthalpy and temperature, respectively; the following two-step procedure results: In the momentum equation, the operators are also split into two parts.The first part is associated with the velocity transport by convection/diffusion written in linearized form as The second part deals with the pressure gradient A 2 = grad.We note that the gradient and divergence operators are adjoints of each other (i.e., A * 2 = −div).The additive scheme of splitting looks like the following: where f n is the right-hand side.This numerical scheme is used as the predictor-corrector procedure.That is, introducing the pressure correction in Equations 14-15 leads to the well-known Poisson equation and the equation for velocity correction in the form of Equations 9-10.
In computational mathematics two variants of fictitious domain methods are recognized: continuation of coefficients at lower-order derivatives and continuation of coefficients at the highest-order derivatives.Both approaches are commonly used in computational fluid dynamics involving phase change processes.Here the first variant is employed, which in a physical sense can be considered as inclusion into the momentum equations of the model of a porous medium: Large Eddy Simulation of Thermo-Hydraulic Mixing in a T-Junction 9 Various formulae of c ǫ can be employed for the flow resistance term in these equations.For Equations 16-17, the modified predictor-corrector procedure taking into account the fictitious domain method looks like the following: An iterative method with a Chebyshev set of parameters using a fast Fourier transform (FFT) solver for the Laplace operator as a preconditioner can serve as an alternative to the conjugate gradient method.The application of this approach for solving the elliptical equations with variable coefficients allows one to reach 50 times the acceleration of the commonly used method of conjugate gradients.For solving the convection problem, a regularized nonlinear monotonic operator-splitting scheme was developed [15].A special treatment of approximation of convection terms C(v) results in the discrete convective operator, which is skew-symmetric and does not contribute to the kinetic energy (i.e., is energetically neutral [15]).The numerical scheme is second, and first-order accurate in space and time, correspondingly.The algorithm is stable at a large-enough integration timestep.Details of the validation for the approach on a wide set of both 2D and 3D tests are reported in [29].
In the next section, we present the numerical results computed with Conv3D on a uniform mesh with 40 million nodes.For the sensitivity study (Section 6.2), we provide results from computation on a uniform mesh with 12 million nodes and on a nonuniform mesh with 3 and 40 million nodes with near-wall refinement.The computations on a uniform mesh with 40 million nodes required approximately 4.3 × 10 4 CPU hours on Lomonosov (Intel Xeon X5570 / X5670 processors).

Results
We focus here on a comparison of the experimental data with the numerical results from three simulation codes: Nek5000, CABARET, and Conv3D.In Section 6.2 we study the velocity field sensitivity to the computational mesh and integration time interval.A similar sensitivity study is undertaken in Section 6.3, where we investigate the effects of grid resolution and time integration interval on the reattachment region immediately downstream of the T-junction.A comparison of velocity and temperature spectra for these codes will be reported elsewhere [30].
All results presented here are nondimensionalized with the cold inlet parameters according to Table 1.For reference, the mean and rms quantities for a set of temporal values u n = u(t = t n ), n = 1, . . ., M, are defined as usual: (21)

Comparison with experiment
In this section, we compare the results of the three codes with experimental data.First we look at the inlet profiles in the cold and hot branches of the T-junction and then at the vertical and horizontal profiles of the mean and rms axial velocity downstream of the junction.

Cold and Hot Inlets
Figure 4 shows the profiles of the mean (top) and rms (bottom) streamwise velocity in the cold branch at x = −3.1 (left) and in the hot inlet branch at z = 2.14286 (i.e., z * /D H = 3) (right).The experimental data are plotted with symbols, the blue solid line represents the Nek5000 simulation, and red dashed line are the inlet profiles from the CABARET simulation.Note that the streamwise component of velocity for cold inlet coincides with the x-direction, while the hot flow in the hot inlet is in the direction opposite the z coordinate; hence, −w is plotted for the hot inlet.Also note that experimental data is plotted in the same style as in [2, Figures 5 and 6 or Figures A5 and A8].
We have observed that the normalized experimental data for the cold inlet does not integrate to unity, indicating a discrepancy between the reported flow rate and the LDA measurements of approximately 6%.On the contrary, using the trapezoidal rule, the integrals of Nek5000 profiles for the cold and hot inlets are equal to 1.0027 and 1.3080, respectively, once averaged over and normalized by the corresponding inlet cross-sections.The same integration procedure for CABARET profiles in Figure 4 gives the normalized inlet flow rates equal to 1.014 and 1.330 for the cold and hot inlets, respectively.These values are in a good agreement with the nondimensional values for inlet velocities U = 1.000 and U H = 1.307 given in Table 2.Note that the normalized inlet mass flow rates for Nek5000 were averaged over the ξ = z line (y = 0) and ξ = y line (z = 0) for the cold inlet and over the ξ = x line (y = 0) and ξ = y line (x = 0) for the hot inlet (Figure 4).
Moreover, a comparison of the shape of the inlet simulation profiles with the experimental data reveals a few differences.The shape of the hot inlet profile for the mean and rms velocity for Nek5000 simulation is not as flat as the experimental or CABARET profiles.This difference can be attributed to a particular modeling of the non-fully-developed flow with the recycling technique described in Section 3. Note that the same technique works nearly perfectly for the cold inlet, where the profile of the mean velocity for Nek5000 is in excellent agreement with the experimental data points after multiplication by 0.94 factor to account for the mass conservation uncertainty of 6%.Similarly, the cold inlet rms velocity for Nek5000 data shows good agreement, diverging from the experimental data points only in the near-wall region, which can be explained by the lower Reynolds number of the simulation, Re = 4 × 10 5 (cf.Table 2).These results indicate that further investigation is needed into the effects of the non-fully-developed flow in the hot inlet for Nek5000 simulation.On the contrary, the CABARET simulation models the flat profile at the hot inlet well (Figure 4, right).The cold inlet profiles, however, deviate substantially from the experimental data, especially in the case of rms profiles (Figure 4, left).

Downstream of the T-Junction
We next look at the axial mean u and rms u ′ velocity downstream of the T-junction.
Figure 5 shows a "bird's-eye" view of the mean (u) and rms (u ′ ) velocity profiles downstream of the T-junction at x=0.6, 1.6, 2.6, 3.6, and 4.6.Here the experimental data (black symbols) are contrasted with numerical simulations with CABARET (red), Conv3D (magenta), and Nek5000 for the benchmark submission results (green) and for longer time integration/averaging (blue).Note the T-junction geometry outline and the equal unit scale for the mean and rms axial velocity equal to the velocity scale, namely, the axial mean velocity in the cold inlet branch U = 1.000 (Table 2).For reference, we provide a detailed  horizontal profiles are plotted similarly in Figures 10-13.Detailed comparison with Nek5000 benchmark submission results is reported in Section 6.2 (see Figure 14).
It is encouraging that overall agreement of simulation results for the three codes with the experimential data is good for the both mean and rms axial velocity.The Nek5000 profiles (blue lines) of the rms velocity match experimental data points well (Figures 5 and 6 right).Moreover, the agreement between the simulation and experiment for the mean velocity (left figures) is best at x = 2.6 (Figures 7 and 11) and at x = 1.6 (Figures 6 and 10) apart of two near-wall data points close to at z = 0.5 (Figure 6).This discrepancy is the focus of an investigation described in Section 6.3.The deviation of the mean velocity for Nek5000 results from experimental data further downstream at x = 3.6 and 4.6 in Figures 5, 8 12-13 can be attributed to the lower Reynolds number used in the simulation because of the time constraints and will be the subject of a follow-up study.
On the contrary, CABARET simulation results (red lines) agree with experiment remarkably well further downstream, at x = 3.6 (Figures 8 and 12) and x = 4.6 (Figures 9 and 13), with the notable exception of the best match of rms data with experimental points in the Large Eddy Simulation of Thermo-Hydraulic Mixing in a T-Junction 17 recirculation region at x = 1.6 and z > 0 (Figure 6).However, close to the T-junction, at x=1.6 and 2.6, the CABARET profiles of the mean axial velocity deviate from experimental points at z > 0 (Figures 6-7) and near the centerline y = 0 (Figures [10][11].Similar to CABARET profiles, the Conv3D results (magenta lines) show the most deviation from experimental data for the mean axial velocity at z > 0 (Figures 6-9) and near the centerline y = 0 (Figures 10-13).However, the agreement of Conv3D simulation with experiment is best at x = 4.6 and 0.4 < |y| < 0.5 (Figure 13).

Sensitivity study
To study the effects of increasing resolution and time averaging interval, we performed an additional set of simulations with Nek5000, CABARET, and Conv3D.
The Figures 14-16 highlight the mesh sensitivity of velocity profiles extracted from results of numerical simulations with Nek5000, CABARET and Conv3D, correspondingly.Each figure compares the vertical and horizontal profiles of the mean axial velocity u and its rms u ′ for different meshes with the experimental data.More detailed comparisons can be found in [31].
In addition to the benchmark submission results with N = 7 (i.e., n ≈ EN 3 ≈ 2.1 × 10 7 points), Nek5000 runs were conducted with N = 5 (n ≈ 7.7 × 10 6 points).This case started with the N = 7 results and was run with ∆t = 5 × 10 −4 and averaged over about 110 convective time units (i.e., about 26 seconds).Figure 14 shows the vertical (left) and horizontal (right) profiles of the mean (top) and rms (bottom) axial velocity profiles at x=0.6 (magenta), 1.6 (black), 2.6 (blue), 3.6 (green), and 4.6 (red) with N = 5 (dashed) and N = 7 (solid).The benchmark submission results (i.e., with N=7 and shorter time average) plotted with dash-dotted line at x = 1.6 . . .4.6 are in the excellent agreement with the longer time average run (solid).All profiles agree well with the experimental data points (symbols), especially considering the fact that the Reynolds number of these simulations, Re = 4 × 10 5 , is two times less than that of the experiment (see Table 2).In general, the profiles from the coarse mesh simulation (i.e., N = 5, dashed) are close to the solution for the finer mesh (solid and dash-dotted lines).Note the excellent agreement between the profiles for N = 5 and N = 7 at x = 0.6, where the strength of the reversed flow is close to its peak (Figure 17).The largest deviation in the profiles (up to about 0.2) is observed at x = 1.6 and 0.1 < z < 0.25; thus, a further study is warranted with even finer mesh, say, N = 8 or N = 9.
Similarly, for the CABARET simulations, the vertical and horizontal profiles of u and u ′ are plotted in Figure 15 at x = 1.6 . . .4.6.These figures show results from CABARET calculations on a coarser mesh with 0.5 million points averaged over a half time interval (red dotted) and the full time interval (blue dash-dotted) and on a finer mesh with 4 million points averaged over a half time interval (green dashed) and the full time interval (solid black line).Note the excellent convergence for the mean axial velocity profiles at x = 3.6 and x = 4.6.
For the Conv3D simulations, the vertical and horizontal profiles of u and u ′ are plotted in Figure 16 at x = 1.6 . . .4.6.This figure shows results from Conv3D simulations on a 40 million node uniform (black solid) and nonuniform (blue dash-dotted) mesh, on a 12 million node uniform mesh (green dashed), and on a 3 million node nonuniform mesh (red dotted line).Note the good convergence for the mean axial velocity profiles at x = 1.6 . . .3.6 and z < 0.

Study of reversed flow region
To address the discrepancy between experimental data and simulations near the upper wall at x = 1.6, we have also undertaken a series of Nek5000 simulations to investigate the sensitivity of the reattachment of the recirculation region downstream of the T-junction.
Summarizing the effects of Reynolds number, grid resolution, and time averaging on the reversed flow region, Figure 17 shows time-averaged velocity profiles near the upper wall at z=0.005 and y=0 for Nek5000 simulations at Re = 9 × 10 4 with N = 11 (red), at Re = 6 × 10 with N = 9 (cyan) and at Re = 4 × 10 4 with N = 9 (magenta), N = 5 (black) and N = 7 for benchmark submission results (green) and longer time averaging (blue).Note the dotted lines that correspond to experimental profile measurements and u = 0 value.
Based on this study, we conclude that the reattachment region is insensitive mainly to an increase in Reynolds number and grid resolution.Moreover, results at Re = 4 × 10 4 indicate that the recirculation region is not sensitive to the averaging/integration time interval.Ideally, one should conduct a further investigation at higher Reynolds number  with longer integration/averaging time to resolve the discrepancy between the simulations and experimental data points near the upper wall at x = 1.6.For now, we quote that in the experimental measurements of velocity, ". . . the focus . . . is not the near-wall region" [2, p. 15].This underscores a necessity of conducting concurrent experiments and validation simulations in which discrepancies like the extent of reciculation region or flow rate uncertainty arising here could be easier to detect and additional data acquisitions and simulations could be performed in order to resolve them.

Conclusions and future work
Several fully unsteady computational models in the framework of large eddy simulations (LES) are implemented for a thermohydraulic transport problem relevant to the design of nuclear power plant pipe systems.Specifically, numerical simulations for the recent 2010 OECD/NEA Vattenfall T-junction benchmark problem concerned with thermal stripping in a T-junction have been conducted with three numerical techniques based on finite-difference implicit LES (CABARET code), finite-volume LES (CONV3D code) and spectral-element (Nek5000 code) approaches.The simulation results of all three methods, including the blind test submission results of Nek5000, show encouraging agreement with the experiment despite some differences in the operating conditions between the simulations and the experiment.In particular, the Nek5000 results tend to closely match the experiment data close to the T-junction, and the CABARET solution well captures the experimental profiles farther downstream of the junction.The differences in the operating conditions between the simulation and the experiment include the uncertainty of the mass flow rate reported in the experiment (about 6%), a reduction in the effective Reynolds number used in some of the simulations that was needed to speedup the computations (Nek5000), and some discrepancies in the inflow boundary conditions.Nevertheless, the good level of agreement between the current simulations and the experiment despite these discrepancies indicates that the flow dynamics in the T-junction experiment is driven mainly by large-scale mixing effects that were not very sensitive to the differences in the operating conditions.
For further investigation, the improvement of the quality of turbulent inflow conditions (most notably, at the hot inlet boundary for Nek5000 and at the cold inlet for CABARET) is planned as well as running additional high-resolution simulations to capture the high Reynolds number regimes typical of the experiment.A further study also will be devoted to a detailed analysis of the temporal spectra of velocity and temperature fluctuations obtained from all three LES solutions.

Figure 3 .
Figure 3. Computational grid used with the CABARET method: full domaim (a), pipe inlet (b,d), and mixed-grid elements in the vicinity of the junction (c,e) for the 0.5 and 4 million cell grid, respectively.
Mean profile in the cold inlet (b) Mean profile in the hot inlet (c) Rms profile in the cold inlet (d) Rms profile in the hot inlet

Figure 4 .
Figure 4. Mean and rms velocity profiles in the cold inlet branch and in the hot inlet branch vs a centerline coordinate ξ for the experimental data (symbols) and simulation results with Nek5000 (blue solid line) and with CABARET (red dashed line).

Figure 5 .
Figure5.T-junction experimental data (•), CABARET (red), and Conv3D (magenta) results, and Nek5000 simulation for the benchmark submission results (green) and for the calculation with a longer time integration (blue line).From top to bottom: vertical profiles of axial mean (u) and rms (u ′ ) velocity and their horizontal profiles.

Figure 7 . 6 Figure 8 . 6 Figure 9 .
Figure 7. Axial mean velocity and rms vertical profiles at x=2.6 for experimental data (circles) and simulation with Nek5000 (blue), CABARET (red), and Conv3D (magenta line).comparison of numerical simulations with experimental data points of the mean and rms axial velocity for each cross-section separately (Figures6-13).Each figure shows the results of numerical simulations with Nek5000 (blue), CABARET (red), and Conv3D (magenta) against the experimental data points (symbols).The vertical profiles of the mean u (left) and rms u ′ (right) axial velocity are shown in Figures6-9at x = 1.6 . . .4.6, correspondingly, while the

Table 1 .
• C entering the main branch and flow at 36 • C entering the side branch.Dimensional Parameters for Vatenfall T-junction Experiment

Table 2 .
Nondimensional Parameters for Vatenfall T-junction Experiment 4Mean and rms profiles of axial velocity at x=1.6 (black), 2.6 (blue), 3.6 (green), and 4.6 (red) for the experimental data (symbols) and CABARET simulation on a coarser mesh with 0.5 million ponts averaged over a half (dotted) and full (dash-dotted) interval, and on a finer mesh with 4 million points averaged over a half (dashed) and full (solid line) time interval.