Numerical Simulation of Non-Isothermal Mixing Flow Characteristics with ELES Method

: Thermal fatigue caused by turbulent thermal mixing in tee pipes is always one of the failure factors of industrial pipes. At present, computational ﬂuid dynamics (CFD) is still the main research method to study the thermal fatigue mechanism. Due to the limitations of the large eddy simulation (LES) model and the classical Reynolds Averaged Navier–Stokes (RANS) model in simulating thermal mixing, an advanced Embedded LES (ELES) model was developed. By comparing the model with data in the open literature, the validity of the ELES model to iso-thermal mixing was evaluated and proven. After that, the ﬂow characteristics of the backﬂow upstream with different momentum ratios ( M R ) were studied using the ELES method, as well as the temperature characteristics near the wall where the backﬂow appears. It was found that the characteristics of the backﬂow and the temperature distribution upstream in the tee were different with different M R values and some regions under speciﬁed M R values are found to be more prone to thermal fatigue at the intersection of the tee upstream. Moreover, the frequency analysis at speciﬁed points near the wall under three different M R values was estimated to evaluate thermal fatigue and the results showed that long-period ﬂuctuations of lower frequencies than 6 Hz upstream were observed. This work helps form a comprehensive understanding of the backﬂow in thermal mixing and the relationship between fatigue and backﬂow in the tee.


Introduction
Since several leakages caused by thermal striping in industrial tee pipes have occurred, notably the accidents that happened in the light water reactor of French PWR Civaux 1, 1998 [1], and in the sodium-cooled fast reactor of PHENIX, 1991 [2], they have aroused wide attention of scholars to focus on the thermal fatigue caused by turbulent thermal mixing. In order to study the relationship between temperature and crack, various experiments have been conducted for temperature fluctuation in the mixing tees. For example, Simoneau et al. [3] and Wakamatsu et al. [4] have conducted a scaled model experiment related to Civaux 1. Furthermore, water experiments (WATLON) have been conducted by Kamide et al. [5] in the Japan Atomic Energy Agency (JAEA) to study the mixing behavior in tee pipes using particle image velocimetry (PIV) and thermocouples. However, due to the limitations of the experimental conditions and the means of observation, computational fluid dynamics (CFD) has become a cheaper and more efficient method to study the flow characteristics of thermal mixing.
So far, many numerical methods have been used to study the turbulence characteristics of mixing flows. Undoubtedly, the Reynolds-Averaged Navier-Stokes (RANS) turbulence models are the most classical and basic. Frank et al. [6], for example, employed RANS turbulence models to simulate the non-isothermal mixing flows in tee pipes using ANSYS CFX software. It was found that the RANS model can accurately simulate the turbulent mixing of isothermal fluids in tee pipes but failed in non-isothermal mixing flows. Additionally, throughout the past few decades, the large eddy simulation (LES) model has been the most used model to simulate free shear flow. For example, Hu and Kazimi [7] carried out a numerical study with LES to investigate the impact of the temperature difference between the mixing flows on temperatures downstream of the tee pipes. Odemark et al. [8] performed an LES analysis of turbulent mixing in a tee pipe with different flow ratios. Tanaka et al. [9] and Utanohara et al. [10] ran numerical simulations with LES approaches for WATLON of a mixing tee in JAEA to investigate the mechanism of thermal fatigue based on the work of Kamide et al. [5]. Other studies using LES were performed by Lee et al. [11], Jayaraju et al. [12], and Kuhn et al. [13] on different aspects of mixing flows. In recent years, several kinds of hybrid RANS-LES models such as Delayed Detached Eddy Simulation (DDES) and Scale-Adaptive Simulation (SAS) were developed for resolving turbulent structures. Zeng et al. [14], Krumbein et al. [15], and Chang et al. [16] employed different hybrid RANS-LES models to simulate the mixing flows in tee pipes. It was found that the results of hybrid RANS-LES models are very close to the experimental data and that of LES.
However, the RANS model, the LES model, and the two hybrid RANS-LES models (DDES and SAS) have their own limitations. Compared with the LES model, there are two main limitations of the RANS model. The first one is the lack of additional information obtained from the RANS simulation, especially for unsteady mixing flows with different temperatures. The second one is related to accuracy. According to the previous study, it is evident that the performance of LES models for free shear flows is significantly more consistent than that of RANS [17,18]. As for the LES model, the high computational cost required for LES at high or even at moderate Reynold numbers is unacceptable [19] because the turbulence scales in LES formulations are of the order of the shear layer thickness, which makes the turbulence length scale extremely small in relation to the thickness of the boundary layer close to the wall [20]. With regard to the hybrid RANS-LES models (DDES and SAS), they may not always be suitable, particularly for flows without strong enough instability to produce turbulent structures on their own [21]. In some conditions, the instability of the separating shear layer cannot be captured by the models, even though resolved turbulence is described at the entrance because the models typically switch back to RANS mode after some boundary layer thicknesses.
Nowadays, a new hybrid RANS-LES model is introduced, which is called the Embedded Large Eddy Simulation (ELES) approach. The ELES model has solved the above problem through the predefinition of interfaces in the geometric model. By this method, the ELES model can not only reduce the computational cost but can also obtain good predictive accuracy of the solution by distinguishing RANS and LES regions and creating synthetic turbulence at RANS-LES interfaces that have been predefined [22]. According to the literature [19], DDES, SAS, and ELES were compared, and it was revealed that better temperature predictions and stronger resolved turbulence activity were found by using the ELES approach than by using the two other hybrid models.
In order to fully understand the flow characteristics and fatigue mechanism in tee pipes, the non-isothermal mixing flow characteristics were investigated using the ELES method with commercial computational Fluid Dynamics (CFD) software ANSYS FLUENT 17.0 in this paper. Firstly, three typical flow patterns of thermal mixing in tee pipes were simulated using the ELES method, and the results were compared with the experimental data in reference [5] to prove the validity of the ELES model. Secondly, the flow characteristics of the backflow upstream with different momentum ratios were studied by carrying out 12 groups of simulations using the ELES method, and some general rules about that were revealed. Thirdly, the temperature characteristics near the wall where the backflow appears were studied. It revealed that some regions were more prone to thermal fatigue at the intersection of the tee upstream under specified momentum ratios due to the mixing induced by the backflow. Moreover, in order to evaluate the thermal fatigue, a frequency analysis of the fluid temperature fluctuation at specified points near the wall of three typical flow patterns was conducted, and the results showed some long-period fluctuations of low frequencies upstream, which may be related to the characteristics of the backflow. The results of this paper can be helpful for the understanding of the thermal fatigue mechanism caused by the backflow, as well as the optimization and improvement of upstream thermal fatigue in tee pipes.

Governing Equations of LES/WMLES
The ELES method splits the domain into RANS and LES portions. The velocity field U (x, t) is decomposed into a filtered component U (x, t) and a subgrid-scale component u (x, t) in the LES portion as below: Unlike the mean (time-averaged) component U (x, t) in the Reynolds decomposition, here, U (x, t) represents an instantaneous field. The filtered momentum equation can then be derived once the velocity decomposition is brought into the Navier-Stokes equation [23]: where p (x, t) is the filtered pressure field, which contains the body force influence, and the residual stress tensor arising from the subgrid-scale motion is defined as Equation (3): Then, ∂U j ∂t Closure needs to be achieved by additional modeling of τ ij . In addition, the filtered continuity equation is defined as Equation (5): The energy governing equation is defined as Equation (6): where h represents the filtered enthalpy, k is the thermal conductivity, c p is the constant pressure-specific heat of the fluid, Pr t is a turbulent Prandtl number, and µ t is the eddy viscosity.
Back to the closure modeling of the residual stress tensor τ ij , it can be computed based on the Boussinesq hypothesis, and the expression is as Equation (7) [24]: where µ t and S ij are the eddy viscosity and the filtered rate-of-strain tensor, respectively. S ij is defined as Equation (8): The Wall-Modeled LES (WMLES) method is adopted for the entire LES domain and reduces the stringent and Re number-dependent grid resolution requirements of classical wall-resolved LES [19]. The original Algebraic WMLES formulation was put forward by Shur et.al. [25], which incorporates several models, such as a modified Smagorinsky model [26], the wall-damping function of Piomelli [27], and a mixing length model. The eddy viscosity µ t in the Shur et al. model is defined as Equation (9).
where d w is the wall distance, κ = 0.41 is the Von Karman constant, and C Smag = 0.2. S is the strain rate, and y + is the normal wall inner scaling. ∆ is defined as Equation (10).
Here, h max is the maximum edge length for a rectilinear hexahedral cell, h wn is the wallnormal grid spacing, and C w is a constant 0.15.

Governing Equations of SST k-ω Model
The k-omega SST model was used in the RANS portion to obtain the inlet velocity and turbulence quantities. The equations of the turbulence kinetic energy (k) and the specific dissipation rate (ω) for the SST k-ω model are as follows [28,29]: The turbulent viscosity µ t and terms such as G k , Y k , G ω , Y ω , D ω are defined as:

Physical Model
The physical problem was that two streams of water at 48 • C with Pr = 3.7 and 33 • C with Pr = 5.1 flowed into the tee pipe from different inlets and mixed. The water experiments were conducted by Kamide et.al. In the experiment, hot and cold water entered the test section from a horizontal main pipe with a diameter of 150 mm (D m ) and a vertical branch pipe with a diameter of 50 mm (D b ), respectively. The main pipe of the test section had 18 D m of the upstream length to the tee and it was 10 D b in the branch pipe. In this study, the computational domain was established with entrance lengths of 650 mm (4.3 D m ) and 175 mm (3.5 D b ) from the tee with the same diameters as the test section. The downstream length was set as 450 mm (3 D m ) in the computational domain to cover the test results. A schematic diagram of the computational domain is shown in Figure 1. The cross-section at 1 D m upstream of the tee of the main pipe was set as a RANS-LES interface, while the cross-sections at 0.5 D m and 1 D m downstream of the tee identified in Figure 1 were consistent with the test as monitor positions. In the main pipe, the upstream region divided by the RANS-LES interface was the RANS computational domain, while another one was the LES computational domain. Similarly, the cross-section at −2.3 D b upstream of the tee of branch pipe was set as a RANS-LES interface. Besides, the coordinate origin was located at the intersection of pipe axes as shown in Figure 1.
downstream length was set as 450 mm (3 Dm) in the computational domain to cover the test results. A schematic diagram of the computational domain is shown in Figure 1. The cross-section at 1 Dm upstream of the tee of the main pipe was set as a RANS-LES interface, while the cross-sections at 0.5 Dm and 1 Dm downstream of the tee identified in Figure 1 were consistent with the test as monitor positions. In the main pipe, the upstream region divided by the RANS-LES interface was the RANS computational domain, while another one was the LES computational domain. Similarly, the cross-section at −2.3 Db upstream of the tee of branch pipe was set as a RANS-LES interface. Besides, the coordinate origin was located at the intersection of pipe axes as shown in Figure 1.

Precursor Simulations
When the ELES simulations are carried out, unsteady fluctuations at the RANS-LES interface are required to be provided to the LES domain [19]. Therefore, precursor RANS simulations are conducted for the profiles of inlet turbulence quantities and the pre-estimation of mesh size in the central mixing zone, and the K-omega SST turbulence model is employed before carrying out the ELES. In the precursor RANS simulations, turbulent flows are conducted in pipes with a diameter of 1 Dm and 1 Db, respectively. In order to obtain the fully developed inlet profiles, the pipe lengths are set as 18 Dm and 10 Db, respectively, as the requirement stated in reference [19]. Through the precursor simulations, the profiles of inlet turbulence quantities can be obtained.
Precursor simulations of three cases with typical flow patterns are simulated in the pipes with the K-omega SST turbulence model. Details of the three cases are listed in Table  1. According to reference [5], the flow pattern can be predicted by the momentum ratio MR: MR < 0.35: impinging jet, 0.35 < MR < 1.35: deflecting jet, MR > 1.35: wall jet. Furthermore, the momentum ratio MR is defined as Equation (17):

Precursor Simulations
When the ELES simulations are carried out, unsteady fluctuations at the RANS-LES interface are required to be provided to the LES domain [19]. Therefore, precursor RANS simulations are conducted for the profiles of inlet turbulence quantities and the pre-estimation of mesh size in the central mixing zone, and the K-omega SST turbulence model is employed before carrying out the ELES. In the precursor RANS simulations, turbulent flows are conducted in pipes with a diameter of 1 D m and 1 D b , respectively. In order to obtain the fully developed inlet profiles, the pipe lengths are set as 18 D m and 10 D b , respectively, as the requirement stated in reference [19]. Through the precursor simulations, the profiles of inlet turbulence quantities can be obtained.
Precursor simulations of three cases with typical flow patterns are simulated in the pipes with the K-omega SST turbulence model. Details of the three cases are listed in Table 1. According to reference [5], the flow pattern can be predicted by the momentum ratio M R : M R < 0.35: impinging jet, 0.35 < M R < 1.35: deflecting jet, M R > 1.35: wall jet. Furthermore, the momentum ratio M R is defined as Equation (17): Quite a number of suggestions on the meshing requirements of LES have been given in the related literature. NJ Georgiadis et al. [30] gave the usual values of the grid spacing in three directions: ∆z + = 50-150, ∆x + = 15-40, and N y = 60-80 beginning with ∆y + ≤ 1, where x, y, and z are the spanwise direction, the wall normal, and the streamwise, ∆i + is the non-dimensional grid spacing defined as ∆i + = u τ ∆i/v (i refers to x, y, and z), and N y is the number of cells to cover the boundary layer in the wall-normal direction. This is a typical requirement for classic wall-resolved LES, and the values of the actual grid spacing (∆i) would be Re-dependent in practice.
Meanwhile, in applications with a wall-modeled LES formulation, the requirement can be relaxed. In the wall-normal direction, although a coarser ∆y + value such as ∆y + ≈ 20 [31] can be tolerable, a near-wall resolution of ∆y + = 1 was taken in the current work and the number of cells across the boundary layer was set as N y = 60 as recommended by ANSYS. A typical wall-normal grid spacing ∆y of the main pipe was set as 0.02 mm to ensure the near-wall resolution in the current work had a Re number of approximately 3.9 × 10 5 . In the other two directions, the grid spacing ∆x ≈ ∆z ≈ (0.05-0.1)δ was adopted in the current work as recommended by ANSYS, where δ (mm) is the per boundary layer thickness. The ∆x/∆z of the main pipe was fixed at 2.25 mm/3.75 mm in the current work, and (∆x + , ∆z + ) ≈ (113, 188) in wall units with a Re number of approximately 3.9 × 10 5 .

Meshing Requirement in the Mixing Zone Downstream in the LES Portion
To meet the mesh resolution of the separating shear layer in the mixing zone downstream, Addad et al. [32] recommended a minimal LES mesh resolution with which the mesh size ∆ was taken as max (λ R , L t R /10), where the Taylor microscale λ R and the turbulent energy length L t R were both obtained from RANS estimations, as seen in Equations (18) and (19).
Kuczaj et al. [33] concluded that the mesh size of ∆ ≈ λ/3, where λ was the Taylor microscale obtained from a posteriori LES analysis, can give sufficiently accurate results, and the value of λ/3 was found to be equivalent to the RANS Taylor microscale λ R . X Fang et al. [31] recommended a mesh size of ∆/η R < (50-100) for LES mixing layers, where the Kolmogorov scale η R can be calculated as Equation (20): Here, when non-equilateral cells are used, the mesh scale ∆ refers to max (∆x, ∆y, ∆z); in addition, k, ε, and ω are the turbulent kinetic energy, the turbulent dissipation energy, and the specific dissipation rate, respectively. Eventually, the mesh scale of the free shear portion was set as Equation (21): As precursor RANS simulations were performed, a typical mesh size was obtained by the above estimator, that is ∆ = 2.5 mm in the free shear portion. Then, hexahedral structure grid elements were constructed in the current work.

Meshing Requirement in the RANS Portion
After the details of meshing requirements in the LES portion are emphasized, the resolution requirements in the RANS portion should be determined. It is recommended that the wall boundary layer should be covered with (20-30) cells with y +~1 , and the minimum resolution across free shear flows can follow Equation (22).
Therefore, we set ∆ = 7.5 mm in the RANS portion as the coarser mesh can be accepted in this domain.
Finally, a typical grid structure with a total of approximately 4 million elements is shown in Figure 2. The A-A cross-section is in the LES portion and the B-B cross-section is in the RANS portion. It can be seen that the mesh in the RANS portion is much coarser than that in the LES portion, which will reduce the number of grids in the whole computational domain, as well as the computation cost. Therefore, we set Δ = 7.5 mm in the RANS portion as the coarser mesh can be accepted in this domain.
Finally, a typical grid structure with a total of approximately 4 million elements is shown in Figure 2. The A-A cross-section is in the LES portion and the B-B cross-section is in the RANS portion. It can be seen that the mesh in the RANS portion is much coarser than that in the LES portion, which will reduce the number of grids in the whole computational domain, as well as the computation cost. In order to assess the grid independence, four types of grids shown in Figure 3 were made according to the above criteria. They had approximately 0.5, 1.6, 4.0, and 6.4 million cells, and the information is listed in Table 2 (cases are defined as Case-A, Case-B, Case-C, and Case-D in ascending order of the grid's number). The cell thickness of the first layer from the pipe's inner surface is marked as Δy1. Additionally, the grid independence is covered in more detail in Section 3.1. Table 2. Information for grid independence.

Treatment near the Wall Treatment in the Free Shear Layer Criteria
Ny = 30, Δy1 = 0.02, y + ~ 1 Δ < 5 More relaxed requirements put forward by NJ Georgiadis et al. [28] and Equation (21) Ny = 60, Δy1 = 1, y + ~ 50 Δ < 2.5 Strict requirements put forward by NJ Georgiadis et al. [28] and Equation (21) Ny = 60, Δy1 = 0.02, y + ~ 1 Δ < 2.5 More relaxed requirements put forward by NJ Georgiadis et al. [28] and Equation (21) Ny = 80, Δy1 = 0.02, y + ~ 1 Δ < 1.5 Strict requirements put forward by NJ Georgiadis et al. [28] and Equation (21) In order to assess the grid independence, four types of grids shown in Figure 3 were made according to the above criteria. They had approximately 0.5, 1.6, 4.0, and 6.4 million cells, and the information is listed in Table 2 (cases are defined as Case-A, Case-B, Case-C, and Case-D in ascending order of the grid's number). The cell thickness of the first layer from the pipe's inner surface is marked as ∆y 1 . Additionally, the grid independence is covered in more detail in Section 3.1.

Boundary Conditions and Numerical Setting
The vortex method [34] was adopted to generate turbulent fluctuations at RANS-LES interfaces in the simulations. The number of vortices was set to be large enough to affect all spots on the face zone. For the inlet boundaries, the profiles of velocity and turbulence quantities obtained from precursor RANS simulations were specified. For the solid boundaries, the pipe walls were set under non-slip and adiabatic conditions.
As recommended in reference [19], the time steps should be chosen to achieve a Courant number of CFL~1 in the LES portion. Therefore, the time step was set as 0.001 s in our simulations, which leads to CFL~1 in the central mixing zone. The mass conservation was most concerned in the process of convergence, and the mass residuals decreased by more than 10 orders of magnitude each time step with 20 inner iteration loops. Other numerical settings and related information in the current simulations are listed in Table  3.
As for the physical properties of the water, the density, specific heat capacity, and coefficient of thermal conductivity change slightly with the temperature, hence, the three physical properties are treated as constant. However, the dynamic viscosity changes a great deal with the temperature and it seems that the dynamic viscosity is approximately

Case Treatment near the Wall Treatment in the Free Shear Layer Criteria
Case-C N y = 60, ∆y 1 = 0.02, y +~1 ∆ < 2.5 More relaxed requirements put forward by NJ Georgiadis et al. [28] and Equation (21) Case-D N y = 80, ∆y 1 = 0.02, y +~1 ∆ < 1.5 Strict requirements put forward by NJ Georgiadis et al. [28] and Equation (21) 2.6. Boundary Conditions and Numerical Setting The vortex method [34] was adopted to generate turbulent fluctuations at RANS-LES interfaces in the simulations. The number of vortices was set to be large enough to affect all spots on the face zone. For the inlet boundaries, the profiles of velocity and turbulence quantities obtained from precursor RANS simulations were specified. For the solid boundaries, the pipe walls were set under non-slip and adiabatic conditions.
As recommended in reference [19], the time steps should be chosen to achieve a Courant number of CFL~1 in the LES portion. Therefore, the time step was set as 0.001 s in our simulations, which leads to CFL~1 in the central mixing zone. The mass conservation was most concerned in the process of convergence, and the mass residuals decreased by more than 10 orders of magnitude each time step with 20 inner iteration loops. Other numerical settings and related information in the current simulations are listed in Table 3. As for the physical properties of the water, the density, specific heat capacity, and coefficient of thermal conductivity change slightly with the temperature, hence, the three physical properties are treated as constant. However, the dynamic viscosity changes a great deal with the temperature and it seems that the dynamic viscosity is approximately linear with a temperature between 33 • C and 48 • C; therefore, a UDF was constructed, and the detailed code can be seen in Appendix A.

Meshing and Independence Validations for Grid
The different mesh requirements used in simulation may influence the computational accuracy and the computational time. Therefore, a grid independence check was performed before the numerical computation to establish an appropriate grid structure in order to maintain a balance between computational accuracy and time. Four cases using the different mesh requirements listed in Table 2 were selected for comparison under the same computational conditions. Additionally, two statistical variables, which were mean axial velocity V Z and axial velocity fluctuation intensity at position Z = −0.53 D m in the W-ref case, were taken as a comparison. Both the mean velocity and the intensity were normalized by the bulk velocity in the main pipe, V m . Definitions are as below: where N is the total number of sampling data and V z i is an instantaneous value at a certain position along the vertical line.
Considering that the frequency of the particle image velocimetry (PIV) used in reference [5] was 15 Hz (0.066 s) and approximately 17 s for the measurement, just for comparison, the sampling interval was set as 0.066 s (66 time-steps) and a period of approximately 10 s was taken for statistics after calculation convergence in the current work.
The simulation results of four grids are shown in Figure 4. It shows that the results reproduce the tendency of the distribution pretty well. At −0.53 D m upstream of the tee, there is no big difference between the mean axial velocity among the grids. Although the results of the four cases vary to some degree, the calculated mean axial velocity is in good agreement with the measured data. However, the results of the axial velocity fluctuation intensity in Case-C and Case-D show better than that of Case-A and Case-B. Case-A underestimated the experimental data, while Case-B performed worse predictions near the wall than others. Considering the computational cost, Case-C seems the best choice among the four grids since it not only has high enough accuracy but also saves computational resources. Therefore, the grid structure of Case-C was selected in the present research.

Turbulence Generation and Structure
Three cases with typical flow patterns listed in Table 1 were simulated with ELES. Before processing the validation of the ELES model, some work needs to be completed such as checking the turbulence generation and structure.
By vorticity post-processing, Figure 5a shows that the resolved turbulences were well-generated from the upstream interfaces in the three reference cases. Here, the vorticity (omega) has projections in the Cartesian coordinate system as seen in Equation (25): Additionally, the turbulence structures in the current ELES simulations need to be checked. Figure 5b shows iso-surfaces of the Q-criterion colored with V z in three cases. Turbulence structures are depicted to be well-developed and show no damping or disruption downstream. The definition of Q is defined as Equation (26): where S is the strain rate and Ω is the vorticity, and both are absolute values. The isosurfaces of the Q-criterion shown in Figure 5b are equal to 200 [s −2 ]. The above analysis reveals that the current ELES simulations provide high resolution on turbulence structures in LES zones with the vortex method adopted at the upstream interface.

Turbulence Generation and Structure
Three cases with typical flow patterns listed in Table 1 were simulated with ELES Before processing the validation of the ELES model, some work needs to be completed such as checking the turbulence generation and structure.
By vorticity post-processing, Figure 5a shows that the resolved turbulences were well-generated from the upstream interfaces in the three reference cases. Here, the vorti city (omega) has projections in the Cartesian coordinate system as seen in Equation (25): Additionally, the turbulence structures in the current ELES simulations need to be checked. Figure 5b shows iso-surfaces of the Q-criterion colored with Vz in three cases Turbulence structures are depicted to be well-developed and show no damping or dis ruption downstream. The definition of Q is defined as Equation (26): where S is the strain rate and Ω is the vorticity, and both are absolute values. The iso surfaces of the Q-criterion shown in Figure 5b are equal to 200 [s −2 ]. The above analysi reveals that the current ELES simulations provide high resolution on turbulence struc tures in LES zones with the vortex method adopted at the upstream interface.

The Contours of Temperature and Velocity Distributions Downstream
The distributions of temperature and fluctuation intensity on a half-cylindrical surface 1 mm apart from the main pipe wall and two cross-sections downstream (as shown in Figure 1) were obtained by simulation. Figures 6-8 show comparisons of these statistical results and the experimental ones. The left sides of Figures 6-8 are the LES results of Kamde et al. [5]. using the AQUA code, which is an in-house thermal hydraulic simulation code [35]. The right sides are the ELES simulation results using FLUENT 17.0.
The data displayed in Figures 6-8 were normalized with the temperature difference as seen in Equations (21) and (22): Here, Tb is the branch temperature and (Tm − Tb) is the temperature difference fixed at 15  The distributions of temperature and fluctuation intensity on a half-cylindrical surface 1 mm apart from the main pipe wall and two cross-sections downstream (as shown in Figure 1) were obtained by simulation. Figures 6-8 show comparisons of these statistical results and the experimental ones. The left sides of Figures 6-8 are the LES results of Kamde et al. [5]. using the AQUA code, which is an in-house thermal hydraulic simulation code [35]. The right sides are the ELES simulation results using FLUENT 17.0. flow direction downstream in all three cases, the temperature range narrowed and maximum fluctuation intensity decreased, which suggested gradually sufficient mi was obtained with the flow.
The main difference between the ELES results and the LES ones conducted by mide et al. [5] shown in Figures 6-8 is that the temperature distributions near the w ELES are much clearer than that of LES, which indicates that ELES is better able to cap temperature fluctuations near the wall. As can be seen from the distribution results in radial cross-sections, along the main flow direction downstream in all three cases, the temperature range narrowed and the maximum fluctuation intensity decreased, which suggested gradually sufficient mixing was obtained with the flow.
The main difference between the ELES results and the LES ones conducted by Kamide et al. [5] shown in Figures 6-8 is that the temperature distributions near the wall of ELES are much clearer than that of LES, which indicates that ELES is better able to capture temperature fluctuations near the wall.    Figure 1).
The calculated distributions were in good agreement with those in the experim The fluctuation intensities were well-simulated to obtain peak values at positions w the physical quantity had a sudden change along the vertical direction. The distrib trends of temperature were similar to that of the axial velocity component, and both The data displayed in Figures 6-8 were normalized with the temperature difference as seen in Equations (21) and (22): Here, T b is the branch temperature and (T m − T b ) is the temperature difference fixed at 15 • C, T is the time-averaged temperature at a certain position on a surface, and N is the total number of sampling data. Considering that the sampling frequency was 100 Hz and the measurement time period was 480 s in reference [5], just for comparison, the sampling interval was set as 0.01 s (10 time-steps), and a period of approximately 100 s was taken for statistics after calculation convergence in the current work. It can be found that the numerical values and distributions in Figures 6-8 show good agreement. It can be clearly seen in Figures 6-8 that the thermal stripping phenomenon appears and the high fluctuation intensity occurs around the cold fluid in contact with the hot. In the W-ref case, the high-intensity region touched the pipe wall to form a half-circle big-value area at each radial cross-section, while in the deflecting jet case, the high-intensity region was separate from the pipe wall because the cold fluid from the branch pipe flowed into the central part of the main pipe. On the other hand, the cold jet from the branch pipe reached the opposite side pipe wall in the impinging jet case, and the high-intensity region then appeared in the upper half of the main pipe and decreased rapidly in the axial direction.
As can be seen from the distribution results in radial cross-sections, along the main flow direction downstream in all three cases, the temperature range narrowed and the maximum fluctuation intensity decreased, which suggested gradually sufficient mixing was obtained with the flow.
The main difference between the ELES results and the LES ones conducted by Kamide et al. [5] shown in Figures 6-8 is that the temperature distributions near the wall of ELES are much clearer than that of LES, which indicates that ELES is better able to capture temperature fluctuations near the wall. In addition, the peak values of the temperature fluctuation intensity were overestimated in ELES analyses, especially at the downstream position Z = 1 Dm (the same in AQUA analyses performed by Kamide et al. [5]), which indicated that it became thermal mixing in the actual flow likely because of the enhancement of the velocity field disturbance downstream during the experiment. In addition, it obtained more appropriate results near the wall with ELES, especially for the capture of the velocity fluctuation intensity, than that in AQUA analyses performed by Kamide et al. [5]. The calculated distributions were in good agreement with those in the experiment. The fluctuation intensities were well-simulated to obtain peak values at positions where the physical quantity had a sudden change along the vertical direction. The distribution trends of temperature were similar to that of the axial velocity component, and both their obvious features appeared at the position around Y = −0.2 D m . The velocity fluctuation intensities were slightly underestimated within ELES simulations, and it was considered that the fluctuations were actually enhanced in the experiments due to the dense arrangement of quite a number of thermocouples.
In addition, the peak values of the temperature fluctuation intensity were overestimated in ELES analyses, especially at the downstream position Z = 1 D m (the same in AQUA analyses performed by Kamide et al. [5]), which indicated that it became thermal mixing in the actual flow likely because of the enhancement of the velocity field disturbance downstream during the experiment. In addition, it obtained more appropriate results near the wall with ELES, especially for the capture of the velocity fluctuation intensity, than that in AQUA analyses performed by Kamide et al. [5].

The Phenomena of Backflow
As a comparison to the LES results conducted by Kamide et al. [5] shown in Figure 10a, mean axial and vertical velocities were used to depict the time-averaged streamlines in the y-z plane in the W-ref case, and the characteristics of the jet bending and eddies' formation obtained with ELES are shown in Figure 10b, which show good agreement. On the other hand, due to the prediction ability near the wall of ELES, it was found that backflow happened near the wall around the tee region, which had not been mentioned in the LES results conducted by Kamde et al. [5]. It was noteworthy that these phenomena were not special cases at certain moments as they were mean results. It can be seen as a vortex generated by backflow, which drives the mixing of hot and cold flows.

The Phenomena of Backflow
As a comparison to the LES results conducted by Kamide et al. [5] shown in Figure  10a, mean axial and vertical velocities were used to depict the time-averaged streamlines in the y-z plane in the W-ref case, and the characteristics of the jet bending and eddies' formation obtained with ELES are shown in Figure 10b, which show good agreement. On the other hand, due to the prediction ability near the wall of ELES, it was found that backflow happened near the wall around the tee region, which had not been mentioned in the LES results conducted by Kamde et al. [5]. It was noteworthy that these phenomena were not special cases at certain moments as they were mean results. It can be seen as a vortex generated by backflow, which drives the mixing of hot and cold flows. In order to study the influence of temperature fluctuation caused by the backflow near the wall around the tee region, a y-z plane of W-ref case at X = 0 was taken to form the contour of temperature fluctuation intensity, the result of which is shown in Figure  11. Figure 11a shows the overview of temperature fluctuations around the right-angle structure upstream of the tee, and temperature fluctuations in the region approximately 1 mm away from the wall are shown in Figure 11b. It can be found from Figure 11 that the temperature fluctuation intensity is close to 0.3, which is almost the same downstream of the pipe shown in Figures 6-8. Obviously, the right-angle structure upstream where the welding exists is exposed to temperature fluctuations, which makes the thermal stresssensitive area more prone to fatigue failure. Therefore, the phenomena of the backflow in the tee region are vital to the study of thermal fatigue and need more attention.
In the following, more detailed quantitative analyses are further carried out to form a comprehensive understanding of the backflow in thermal mixing and the relationship between fatigue and backflow in the tee. In order to study the influence of temperature fluctuation caused by the backflow near the wall around the tee region, a y-z plane of W-ref case at X = 0 was taken to form the contour of temperature fluctuation intensity, the result of which is shown in Figure 11. Figure 11a shows the overview of temperature fluctuations around the right-angle structure upstream of the tee, and temperature fluctuations in the region approximately 1 mm away from the wall are shown in Figure 11b. It can be found from Figure 11 that the temperature fluctuation intensity is close to 0.3, which is almost the same downstream of the pipe shown in Figures 6-8. Obviously, the right-angle structure upstream where the welding exists is exposed to temperature fluctuations, which makes the thermal stress-sensitive area more prone to fatigue failure. Therefore, the phenomena of the backflow in the tee region are vital to the study of thermal fatigue and need more attention.

Characteristics of Backflow in the Tee
Backflow is widespread in non-isothermal mixing in the tee, and these phenomena are closely related to the momentum ratio of the main to branch flow. Therefore, to fully study these phenomena and understand the mechanism of the backflow in the tee, 12 In the following, more detailed quantitative analyses are further carried out to form a comprehensive understanding of the backflow in thermal mixing and the relationship between fatigue and backflow in the tee.

Characteristics of Backflow in the Tee
Backflow is widespread in non-isothermal mixing in the tee, and these phenomena are closely related to the momentum ratio of the main to branch flow. Therefore, to fully study these phenomena and understand the mechanism of the backflow in the tee, 12 groups of simulations with different momentum ratios were carried out and the flow field characteristics were analyzed. Detailed information on the simulations is listed in Table 4. The streamline diagrams of the three impinging jet cases were compared and the results are shown in Figure 12a-c. The streamline diagrams were drawn in the y-z plane at X = 0. Figure 12a-c shows that when the flow rate of branch flow is high enough, the upper fluid of the main flow is squeezed into a vortex and the phenomena of the backflow appear near the opposite sidewall. The higher the flow rate, the larger the eddy and the more obvious the phenomena of the backflow. When M R ≤ 0.1, the backflow with large eddies occurs near the opposite sidewall upstream. The results are similar to the experiments conducted by several researchers, e.g., Mei-Shiue Chen et al. [36], C.H. Lin et al. [37], G.Y. Chuang and Y.M. Ferng [38,39], and Wu et al. [40] In their research, the phenomena of the backflow occurred when M R ≤ 0.05, M R ≤ 0.062, M R ≤ 0.06, and M R ≤ 0.05, respectively. The difference may be the limit of observation methods in the experiment regarding the inadequate sensitivity of temperature sensors.
As the velocity of the branch flow decreases, so does the vortex, and when M R = 0.2, the phenomena of the backflow that appear near the opposite sidewall almost disappeared and gradually transitioned to the place of the right-angle structure upstream of the main flow, as shown in Figure 12d. As the velocity of the branch flow decreases and that of the main flow increases, the vortex formed by backflow becomes larger at the right-angle structure upstream of the branch flow, as shown in Figure 12e-f. Compared with the size of the eddies near the opposite sidewall upstream when M R ≤ 0.1, the eddies shown in Figure 12d-f are relatively small, which may be the main reason that researchers have rarely observed the phenomena.

Temperature Characteristics near the Wall Where Backflow Appears
It can be found from the above research that different flow patterns with different momentum ratios have different characteristics of the backflow. In order to study the temperature distributions related to the backflow, detailed studies have been carried out under typical conditions in each flow pattern. From Figure 12, three typical cases with the most obvious characteristics of the backflow were chosen as the representative of each flow pattern, which were Case1 with M R = 0.05, Case7 with M R = 0.8, and Case12 with M R = 32. Additionally, the temperature was recorded at every time step with a sampling frequency of 1000 Hz for a duration of 74.8 s.   Axial distributions of the temperature fluctuation intensity of the above cases along corresponding lines are depicted in Figure 13. The corresponding lines of the three cases described in Figure 13 are 1 mm away from the wall. From the distribution of temperature fluctuation along line1 in Figure 13, we can see that the axial range affected by thermal striping is very wide, and it also seems that there exists more than one peak, with a maximum of approximately at Z = −0.5 D m upstream, which indicates that the thermal mixing caused by the backflow is complex. In addition, the temperature fluctuation intensity at Z = −0.5 D m shown in Figure 13 is close to that of 0.5 D m downstream plotted in Figures 6-8. In order to further study the mixing behavior, two contours were plotted at the cross-sections of Z = −0.5 and −0.2 D m upstream as shown in Figure 14a,b. As we can see from Figure 14a  In order to analyze the frequency characteristics (PSD), the FFT method was applied to the data extracted over a period of 74.8 s and presented in Figure 15. The prominent frequency components were observed in all analyses in Figure 15. At Z = −0.5 Dm, θ = 180° From the distribution of temperature fluctuation along line2 in Figure 13, we can see that the temperature fluctuation intensity is relatively low and the axial range affected by thermal striping is narrow, which indicates that the thermal mixing is not fully developed although there exists a relatively larger eddy as shown in Figure 12d. The main reason for this phenomenon may be that the momentums of the main and branch flows are similar as M R = 0.8. When each flow blocks the other one in the tee with M R~1 , there is not enough power for one to penetrate upstream of the other one. A contour plotted at the cross-section of Z = −0.18 D m upstream as shown in Figure 14c reveals the above view. Therefore, the depth of penetration of one to another is dependent on the momentum ratio.
From the distribution of temperature fluctuation along line3 in Figure 13, we can see that although the axial range affected by the thermal striping is narrow, the temperature fluctuation intensity is relatively higher compared with that of Case1 and Case7. With the penetration depth, the temperature fluctuation intensity tends to increase and then decrease, and a peak of approximately 0.32 appears at the location of Y = −1.55 D b . Figure 14d shows the temperature distributions at Y = −1.55 D b . It can be seen that the phenomenon of thermal striping is very significant, and the radial range affected by the thermal striping is wide as it is from −90 • to 90 • . Meanwhile, the fluctuation peak appears at θ = 45 • in Figure 14d.  In order to analyze the frequency characteristics (PSD), the FFT method was app to the data extracted over a period of 74.8 s and presented in Figure 15. The promi frequency components were observed in all analyses in Figure 15. At Z = −0.5 Dm, θ = of Case 1 in Figure 15, frequency peaks appeared at 0.5 Hz, 0. Although the peaks did not exactly match, were similar in the three locations. Besides, it is generally known that the frequency proaching 6 Hz in Figure 15 corresponds to the frequency of the Karman vortex stree addition, the characteristic frequency was not very significant, at approximately 6 downstream shown in references [5,10], but long-period temperature fluctuation In order to analyze the frequency characteristics (PSD), the FFT method was applied to the data extracted over a period of 74.8 s and presented in Figure 15. The prominent frequency components were observed in all analyses in Figure 15. At Z = −0.5 D m , θ = 180 • of Case 1 in Figure 15, frequency peaks appeared at 0.5 Hz, 0. Although the peaks did not exactly match, they were similar in the three locations. Besides, it is generally known that the frequency approaching 6 Hz in Figure 15 corresponds to the frequency of the Karman vortex street. In addition, the characteristic frequency was not very significant, at approximately 6 Hz downstream shown in references [5,10], but long-period temperature fluctuations of frequencies lower than 6 Hz upstream from the tee were found, which had also been observed by Utanohara et al. [10] downstream. As described in Smith et al. [41] and Thomas [42], long-period temperature fluctuations of a low frequency (0.1-10 Hz) would cause the most harmful thermal loads. Hence, the conditions may require more attention.
At Z = −0.18 D m , θ = 0 • of Case 7 in Figure 15, frequency peaks appeared at 0.49 Hz, 0.89 Hz, 1.47 Hz, 2.75 Hz, 6.0 Hz, etc., while at Y = −1.5 D b , θ = 45 • of Case 12, frequency peaks appeared at 0.95 Hz, 1.55 Hz, 3.2 Hz, 5.76 Hz, etc. The characteristic frequency of approximately 6 Hz in the two cases seems more obvious than that in case1, which is because the phenomena that occurred at the points in case7 and case12 are more similar to the phenomena of the Karman vortex street. Moreover, a long-period temperature fluctuation of a frequency lower than 6 Hz upstream (e.g., approximately 0.5 Hz, 1.2 Hz, 3.0 Hz, etc.) from the tee was observed, too. This may be related to the periodic variation of the backflow around the attached wall. The long-period temperature fluctuation of a frequency lower than 6 Hz downstream from the T-junction was also seen in the experimental investigation by Miyoshi et al. [43] and the numerical study by Utanohara et al. [10]. However, it requires further study to validate the long-period fluctuation of a frequency lower than 6 Hz upstream through more experiments.
Hz, etc.) from the tee was observed, too. This may be related to the periodic variation of the backflow around the attached wall. The long-period temperature fluctuation of a frequency lower than 6 Hz downstream from the T-junction was also seen in the experimental investigation by Miyoshi et al. [43] and the numerical study by Utanohara et al. [10]. However, it requires further study to validate the long-period fluctuation of a frequency lower than 6 Hz upstream through more experiments. In sum, by analyzing the temperature distribution of three different flow patterns with different MR values, three more regions are found to be prone to thermal fatigue at the intersection of the tee upstream. The effects of fatigue are different in the three regions. It seems that the impinging jet with a smaller MR and the wall jet with a larger MR have a greater possibility of fatigue crack on the wall around the intersection, while the effect of fatigue caused by the deflecting jet on the wall seems less significant.

Conclusions
In this paper, the mixing flow of non-isothermal fluids in tee pipes upstream is numerically studied in detail. Firstly, three typical flow patterns of thermal mixing in tee pipes (wall jet with MR = 8.1, deflecting jet with MR = 0.8, and impinging jet with MR = 0.2) were simulated using the ELES method, and the results have been compared with the experimental data in reference [5] to prove the validity of the model. Secondly, the flow In sum, by analyzing the temperature distribution of three different flow patterns with different M R values, three more regions are found to be prone to thermal fatigue at the intersection of the tee upstream. The effects of fatigue are different in the three regions. It seems that the impinging jet with a smaller M R and the wall jet with a larger M R have a greater possibility of fatigue crack on the wall around the intersection, while the effect of fatigue caused by the deflecting jet on the wall seems less significant.

Conclusions
In this paper, the mixing flow of non-isothermal fluids in tee pipes upstream is numerically studied in detail. Firstly, three typical flow patterns of thermal mixing in tee pipes (wall jet with M R = 8.1, deflecting jet with M R = 0.8, and impinging jet with M R = 0.2) were simulated using the ELES method, and the results have been compared with the experimental data in reference [5] to prove the validity of the model. Secondly, the flow characteristics of the backflow upstream with different momentum ratios were studied by carrying out 12 groups of simulations using ELES, as well as the temperature characteristics near the wall where the backflow appears. Moreover, in order to evaluate the thermal fatigue, the frequency analysis at specified points near the wall of three typical flow patterns was estimated. The conclusions are summarized below: (1) ELES is a valid tee pipe numerical turbulent-solving method and aptly predicts the flow patterns of non-isothermal mixing in tee pipes. (2) The characteristics of the backflow in the tee are obviously different with different momentum ratios M R . The relatively general rules are that the eddies caused by the backflow appear near the opposite sidewall of the branch pipe for cases of the impinging jet when M R ≤ 0.1, the eddies appear near the wall at the right-angle structure upstream of the main pipe for cases of the deflecting jet when M R~0 .8, while the eddies appear near the wall at the right-angle structure upstream of the branch pipe for cases of the wall jet when M R ≥ 11.5. In addition, the size of the eddies caused by the backflow is dependent on M R .
(3) The characteristics of temperature distribution and the frequency of temperature fluctuations are different under different backflow patterns with different momentum ratios M R . Besides the pipe wall downstream, three more regions are found to be prone to thermal fatigue at the intersection of the tee upstream. The effects of fatigue are different in the three regions. It seems that an impinging jet with smaller M R and a wall jet with larger M R have a greater possibility of fatigue cracks on the wall around the intersection, while the effect of fatigue caused by the deflecting jet on the wall seems less significant. (4) Long-period temperature fluctuations of frequencies lower than 6 Hz upstream (e.g., approximately 0.5 Hz, 1.2 Hz, 3.0 Hz, etc.) from the tee were observed, and they may be related to the periodic variation of the backflow around the attached wall. Besides, this requires further study to validate the long-period fluctuation of a frequency lower than 3 Hz upstream through further experiments. Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.  Figure A1. The relationship between dynamic viscosity and temperature. Figure A1. The relationship between dynamic viscosity and temperature.