Large-Eddy Simulation of turbulent thermal flow mixing in a vertical T-Junction configuration

International Journal of Thermal Sciences 150 (2020) 106231

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: http://www.elsevier.com/locate/ijts

Large-Eddy Simulation of turbulent thermal flow mixing in a vertical T-Junction configuration Cenk Evrim *, Eckart Laurien Institute of Nuclear Technology and Energy Systems, University of Stuttgart, Pfaffenwaldring 31, 70569, Stuttgart, Germany

A R T I C L E I N F O

A B S T R A C T

Keywords: Large-eddy simulation Turbulent mixing Thermal stratification T-junction Temperature fluctuations Spectral peak

Turbulent mixing of warm and cold water streams within a vertical T-junction configuration of a nuclear power plant piping system causes near-wall temperature fluctuations which may lead to High-Cycle Thermal Fatigue (HCTF) failure of the wall material. To predict frequency and amplitude of such temperature fluctuations accurately the method of Large-Eddy Simulation (LES) in the numerical simulation code OpenFOAM is used. As an experimental test case, the turbulent mixing experiment of Kamide et al. [1], where a vertical branch pipe with a cold water stream is connected from below to a horizontal main pipe with a warm water stream to form a vertical T-junction configuration. From the various observed flow patterns, the temperature and velocity data of the ‘wall jet’ are chosen as a reference because it is known to have the highest potential for thermal fatigue. Incomplete mixing and an instability associated with large turbulent structures near the junction region lead to significant temperature fluctuations close to the wall before a vertical thermal stratification further downstream stabilizes the flow. The highest near-wall fluctuation amplitudes are 26% of the temperature difference between the streams. A spectral peak occurs at a Strouhal number of 0.2. The results of the simulations demonstrate, that the mean velocity and temperature profiles of the experiments are well captured, whereas the root-mean-square (RMS) temperature fluctuations deviate from the measurements in some cases, in particular when coarse grids are used. In order to find a lower limit for the required spatial resolution three different numerical meshes with a variable number of cells up to 28 � 106 are investigated. The results demonstrate that with a sufficient mesh resolution the velocity and temperature distributions as well as the spectral peak can be simulated with good agreement to the experimental data.

1. Introduction 1.1. Motivation In nuclear power plants the phenomenon of safety-relevant thermal fatigue of the piping material is observed. In their coolant circuits pipeline configurations exist, where turbulent water streams with different temperatures are mixed. Such mixing may be incomplete and associated with temperature fluctuations, which may be created near a pipe wall to propagate into the solid structure. Such fluctuations can lead to thermal-fatigue damage of the wall material after a high number of cycles, e.g. to a through crack of a pipe in the nuclear power plant Civaux, see Chapuliot et al. [2]. The fluctuations are considered relevant for thermal fatigue if they are in the frequency range between 0.1 and 10 Hz [2]. The investigation of flow mechanisms leading to such High Cycle

Thermal Fatigue (HCTF) has attracted the attention of researchers, both experimentally and numerically. Turbulent mixing of hot and cold fluids at a 90� T-junction configurations is of primary interest, because this component is particularly prone to HCTF conditions as found by Courtin [3]. T-junction configurations, where a stream from a cold ‘branch’ pipe is injected vertically from above or below into the warmer stream through a larger diameter ‘main’ pipe, are the most common mixing configurations in Light Water Reactors, see Dahlberg et al. [4]. For the analysis and prediction of the frequency and amplitude of the associated temperature fluctuations of turbulence within a wide range of the ma­ terial domain near the T-junction reliable and accurate methods must be developed. Since the fluctuations are generated by turbulence processes, the three-dimensional flow and turbulence simulation and modeling by Computational Fluid Dynamics (CFD) methods is of primary interest.

* Corresponding author. E-mail address: [email protected] (C. Evrim). https://doi.org/10.1016/j.ijthermalsci.2019.106231 Received 21 October 2019; Received in revised form 11 December 2019; Accepted 11 December 2019 Available online 25 December 2019 1290-0729/© 2019 Elsevier Masson SAS. All rights reserved.

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

Nomenclature Cw D f g h PSD p Q Sij Sdij T t u uτ y1

α θ

μ ν ρ

Co MR Re

Pr Prandtl number turbulent Prandtl number Prt St Strouhal number b branch eff effective m main t turbulent RMS root-mean-square þ non-dimensional form * non-dimensional form CFD computational fluid dynamics DILU diagonal incomplete-lu FSI fluid-structure-interaction GAMG geometric-algebraic multi-grid HCTF high-cycle thermal fatigue LDV laser doppler velocity LES large eddy simulation LUST linear-upwind stabilized transport NIST National Institute of Standards and Technology OpenFOAM open source field operation and manipulation PBiGG preconditioned bi-conjugate gradient PIV paricle image velocimetry RANS Reynolds-Averaged Navier-Stokes WATLON water experiment on fluid mixing in t-pipe with long cycle fluctuation

WALE constant [ ] internal diameter of pipe [m] frequency [1/s] acceleration due to gravity [m/s2] enthalpy [J/kg] power spectral density pressure [bar] Q -criterion [1/s2] resolved shear strain rate tensor

traceless symmetric part of the square of the velocity gradient tensor temperature [� C] time [s] velocity [m/s] friction velocity [m/s] nearest node from the wall [m] thermal diffusivity [m2/s] degree [� ] dynamic viscosity [Pa s] kinematic viscosity [m2/s] density [kg/m3] Courant number momentum ratio Reynolds number

1.2. Literature overview

used: (i) bridge the gap between the first cell and the wall by a wall function, or (ii) resolve the unsteady flow near the wall. The method (ii) is often referred to as a ‘wall-resolved LES’. Various simulation for the FSI facility have already been performed by Selvam et al. [14] and Gauder et al. [15]. Jayaraju et al. [16] performed a detailed investigation and showed that with wall functions the velocity and temperature fluctuations tend to be underestimated and have an adverse effect on the heat transfer across the walls. Therefore, it is advisable that the near wall region must be resolved by using the wall-resolved LES. The nature of the flow is then captured explicitly and the near-wall large eddies are directly resolved. The small eddies and the fluctuations associated with them are too fast to be transferred to the structure can therefore be modelled. Increasing computational resources now enable the use of such fine mesh resolu­ tions. However, the wall-resolved LES is still a challenge with increasing Reynolds number Re, because the near-wall layer thickness decreases almost linearly with Re. It is therefore necessary to demonstrate the feasibility of a wall-resolved LES in comparison with experiments, in particular at high Re. Experiments with a dominant peak in the power spectrum are most suitable for this. Miyoshi et al. [17] performed a mixing-tee experiment with a warm side stream injected into a horizontal cold flow in the main branch, denoted as the ‘T-Cubic experiment’. The T-junction installed in the experimental loop was made of steel. A total of 148 thermocouples were installed in the mixing zone within the piping material. The experiments were performed with a temperature in the inlet flow streams of 59.8 � C and 29.7 � C. They investigated also in a fluid-solid coupled simulation its predictive performance for the flow and temperature field. In this work, a dominant peak in the power spectrum at a position in the fluid and solid domain in the HCTF area could be observed, see the numerical work by Utanohara et al. [18]. By Howard and Serre [19] the mixing and thermal fluctuations inducted in another mixing T-junction was analyzed. With the introduction and definition of high-order turbulent statistics this work provides a better understanding of turbulence fea­ tures. It was found that the turbulence transport term is important for

In order to create a database for the development and validation of CFD methods to predict turbulent flow relevant for HCTF various ex­ periments have been performed. By Brückner [5] mixing experiments were carried out in a horizontal T-junction configuration. The mean and instantaneous flow was measured by Particle Image Velocimetry (PIV). The author detected and described vortex structures including counter-rotating vortices directly downstream of the T-junction. Smith et al. [6] carried out thermal mixing experiments in a vertical T-junction configuration. The temperature of the inlet flow streams are set to 36 � C and 19 � C. Well-defined inlet boundary conditions were created by measurements of the velocities in the inlet branches and in the mixing region by PIV and Laser Doppler Velocity (LDV). Addition­ ally, the temperature in the near wall region was measured. With these experimental data as a basis a CFD benchmarking exercise was orga­ nized, see Several numerical works have been performed with different €hne [7] show, that numerical codes [7–9]. The numerical results of Ho simulations on the basis of the Reynolds-Averaged Navier-Stokes Equations (RANS) are unable to predict a realistic mixing behaviour of the two fluids. However, with Large-Eddy Simulation (LES) the velocity profiles and temperature results can be predicted with a good accuracy. Indeed, these studies are performed with rather coarse grids relative to the high Reynold number (Re ¼ 80000 in the main branch and 100000 in the side branch) and the complexity of the flow. To investigate thermal mixing experiments with conditions close to reality in nuclear power plants, the Fluid-Structure-Interaction (FSI) facility was built at the University of Stuttgart, see Kuschewski et al. [10]. The working fluid in the main and branch pipe is water and ex­ periments can be performed with a maximum pressure and temperature about 75 bar and 280 � C. Both streams come together in a T-junction made of austenitic stainless steel [10–13]. In order to capture the near wall flow and heat transfer including the laminar sublayer the heat conduction of the solid material must be taken into account with in coupled model. For the LES two methods can be 2

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

quantifying the turbulence. A T-junction water experiment ‘WATLON’ on fluid mixing in T-pipe with long cycle fluctuations with a dominant peak were performed by Kamide et al. [1]. Here, the Reynolds number in the main branch is 380000 and in the side branch 66000. Georgiou and Papalexandris [20] reported a Direct Numerical Simulation, i.e. a simulation without any turbulence model, of turbulent heat transfer in a T-junction. Here, a dominant peak in the spectral analysis of the temperature signals could also be observed. This inves­ tigation was performed with a T-junction with rectangular cross sections of the main and side branches. However, it should be noted direct nu­ merical simulation is restricted by relatively low Reynolds numbers. The Reynolds number of the cross flow is only Re ¼ 3000. In Ref. [1] the branch pipe of the investigated T-junction is vertically arranged. The thermal mixing experiments were performed with a temperature difference between the branches of ΔT ¼ 15� C. A ther­ mocouple tree was built and used to measure the temperature distri­ bution in the near wall region and in the direct mixing region. The thermocouples are arranged in the radial direction by 5 mm distances from the pipe center and two more thermocouples are installed 1 mm and 3 mm from the wall structure. The velocity was measured by the two-dimensional PIV method. Three different types of flow patterns (wall jet, deflecting jet and impinging jet) were observed in the mixing region. The flow pattern type depends on the momentum ratio between the main and branch pipe. Due to mixing of turbulent flows of different density and velocity at a 90� angle, the mixing interface oscillates and due to hydrodynamic instabilities it breaks down and generates large-scale transient flow structures. The oscillation of these structures leads to low-frequency oscillations of the temperature, refer to Obabko et al. [21]. The solid tends to be insensitive for small-scale turbulent structures which are too rapid high frequency fluctuations. However, large-scale turbulent fluctuations can lead to damaging fluctuations in the solid. In most of the mixing processes discussed above, the effect of gravity or buoyancy on the generation and development of large turbulence structure can be neglected, due to the small temperature difference below 30 K between the branches, e.g. Refs. [1,6,19]. In Ref. [10], however, gravity has a significant effect, leading to a ‘wavy’ density stratification downstream of the mixing region as well as to stabilisation and damping of turbulence. The effects of Ref. [10] investigated sepa­ rately, e.g. in Refs. [11,12], and are not subject of the present investi­ gation. It is therefore unimportant, whether the side branch comes from above or below and which branch is warmer than the other. Due to the considered low temperature differences the thermal coupling between the flow and the solid structure is not considered. The direct mixing region of a T-junction can be separated into fundamental processes as discussed by Nakamura et al. [22]. The ther­ mal mixing of two fluids with a vortex formation generates temperature fluctuations in the mixing region (I). Temperature fluctuations reach and appear in the boundary layer of the flow (II). Thereafter, the fluc­ tuations in the boundary layer are transferred to the inner wall structural surface (III). By heat conduction, the fluctuations on the surface prop­ agate through the wall and generate stress fluctuations inside the structure (IV). This retaking thermal stress fluctuations cause High-Cycle Fatigue Cracking. Thus, thermal fatigue is caused by the following phenomena: turbulent mixing, thermal striping and thermal stratification. Beside investigations where the working medium is liquid water are also several studies on gaseous flows in T-junctions. Hirota et al. [23] performed experiments using PIV on a rectangular T-junction and studied the behaviour of the interface between two flows in the flow-merging region. Additionally, the turbulent heat flux could be obtained by the help of simultaneous measurement of velocity and temperature. They reported that the vertical oscillation motion of the interface is caused by the fluctuation of the streamwise velocity component. Georgiou and Papalexandris [24] performed wall-resolved LES of

turbulent mixing and heat transfer in T-junctions and investigated the role of the temperature as an active scalar. Feng et al. [25] studied thermal mixing characteristics accompanied by complicated direct contact condensation. They reported that the condensation rate is proved closely related to the coolant thermody­ namic ratio between the main and branch pipe. In spite of many numerical and experimental investigations, the prediction of flow phenomena with dominant temperature fluctuations associated with local risk areas for HCTF is still a challenge. To analyze of turbulent statistics a large total physical simulation time is necessary, see Tanaka et al. [26]. So far, in almost all numerical investigations fixed value conditions were used as inflow boundary conditions for the ve­ locity, although the flow is turbulent at the inlet. 1.3. Aim of this study A wall-resolved LES with mapped boundary conditions in order to take account of the turbulence at the inlet is performed. The aim of this work is to demonstrate that a reasonable agreement can be achieved if the mesh is fine enough and the relevant scales in the flow can be resolved. To investigate the influence of the mesh, a mesh-refinement study with different resolutions is performed. Also, the influence of the simulation time on the representation of the dominant peak of the temperature fluctuations in the near wall region is investigated. This peak, which can be an HCTF indicator, depends on the flow conditions, mainly on the velocities in the main and branch pipe and on the velocity ratio. It is best expressed by the non-dimensional Strouhal number. This work includes also a variation of the Reynolds number in order to show its effect on the Strouhal number. An interpretation of the un­ steady flow mechanisms in view of vortex formation, penetration of the temperature fluctuations through the near-wall layers and their transfer to the wall is made. 2. Computational domain and numerical grid The geometrical model and the flow conditions are based on the WATLON experiment [1]. Fig. 1 shows a schematic of the computational domain of the T-Junction. The main pipe diameter is Dm ¼ 0.15 m with an inflow length of 12Dm . The branch pipe diameter is Db ¼ 0.05 m and has an inflow length of 12Db . The downstream of the main pipe is 20Dm . In this work, cartesian coordinates are used. The origin of the coordinate system is located at the center of the mixing tee. The gravitational force

g is set along the negative z-direction. The geometry is meshed using the software ANSYS ICEM-CFD and contains hexahedral elements. For the grid sensitivity study, three mesh resolutions (coarse (C), intermediate (I) and fine (F)) are used, respectively. The meshes are consist of approximately 1.5, 6 and 28 million elements. Fig. 2(a–f) shows mesh refinement and resolution in the symmetry plane and in the yz-plane. In Ref. [1] a LES was performed with a total mesh number of around 0.3 millions and a statistical analysis time of 3s. A spectral peak could not be detected. In this study, the simulations are performed with a higher mesh resolution and an increased total simulation time. The non-dimensional parameter yþ 1 defines the distance from the wall measured in wall units. It is defined by Eq. (1) [27]:

⇀

yþ 1 ¼

y1 uτ

ν

(1)

in which y1 defines the distance from the nearest node from the wall, uτ the friction velocity and ν the kinematic viscosity of the fluid. The first cell height within the boundary layer of all three meshes is 2⋅10 5 m. þ The average yþ 1 value is approximately 0.52 and the largest y1 value is around 1.5. Occasionally, yþ exceeds only slightly the recommended 1 value of yþ 1 ¼ 1, so that it can be still considered as a wall-resolved LES. A comparison of the current LES with similar studies in the literature is 3

International Journal of Thermal Sciences 150 (2020) 106231

C. Evrim and E. Laurien

Fig. 1. Computational domain of the T-junction.

Fig. 2. Numerical grid resolutions for coarse (a, b), intermediate (c, d) and fine mesh (e, f).

shown in Table 1.

the same as in the WATLON experiment [1]. The properties of water are temperature dependent and defined as third order polynomials based on the NIST database. Table 2 illustrates the investigated cases with the different Reynolds numbers and the corresponding velocities in the main and branch pipes. Here, the velocity in the main and branch pipe varies but the velocity ratio ðum =ub ¼ 1:46Þ is constant in all three cases. Furthermore, the momentum ratio ðMR ¼ ð4 *Dm *Db *ρm *u2m Þ =ðπ *D2b *ρb *u2b ​ Þ ¼ 8:1Þ is also constant. The density ratio in this study is 0.994 and the water pressure is 1 bar. The laminar Prandtl number for the working fluids is 3.7 (Tm ¼ 48 � C) and 5 (Tb ¼ 33 � C). For the grid sensitivity study and the validation, experimental results from the velocity and temperature measurement of case A are available.

2.1. Boundary conditions To provide a fully developed turbulent flow in the main and branch pipe entrances a recycling procedure is used. The velocity data from an internal plane at 4Dm and 4Db are mapped to the entrance of the main and branch pipe at every time step [31]. To accelerate the initial tur­ bulence development the flow was initialized with a perturbation method similar to Pandey et al. [32,33]. Here, a laminar velocity profile with fluctuations in the near wall region is initialized. The temperature initialization is performed with the constant high fluid temperature in the main pipe and the constant low fluid temperature in the branch pipe. The temperature difference between the two pipes is 15 � C. Warm water flows in the main pipe (Tm ¼ 48 � C) and cold water flows in the branch pipe (Tb ¼ 33 � C). The boundary conditions of the temperature for the main and branch pipe inlet flows are constant and spatially uniform in time. Boundary conditions on the wall are set as adiabatic for the tem­ perature and no-slip conditions for the velocity. At the outlet, the tem­ perature is set to homogenous Neumann condition and an inletOutlet boundary condition for velocity is used. The investigated conditions are

3. Numerical setup 3.1. Large-Eddy Simulation In an LES, small eddies in the turbulent flow are modelled by using a Table 2 Boundary conditions for the investigated cases, with um =ub ¼ 1:46 and.MR ¼ 8:1

Table 1 Comparison of wall resolution of the current LES with previous studies in the literature. yþ 1;max

Kuhn et al. [28]

Frank et al. [29]

Kl€ oren et al. [30]

Present

5

4.5

1.5

1.5

4

case

Rem ½ �

Reb ½ �

um ½m=s�

ub ½m=s�

A B C

380000 190000 95000

66000 33000 16500

1.46 0.73 0.365

1 0.5 0.25

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

subgrid model and large eddies are directly resolved. Scale separation is performed by a low-pass filter [34]. Z ~ ¼ ρφGðx; ξÞdξ; ρφ ¼ ρφ (2)

standard inputs. The filter parameters are set by cubeRootVolDelta. This model uses the cell volume as filter parameter. The unsteady equations solved with a variable time step and a constant Courant number of Co ¼ 0.5. The total physical simulated time is 30 s and the last 20 s are used for the turbulent statistics.

Ω

in which φ defines any flow variable and Gðx; ξÞ is the filter kernel with the filter width Δ. In case of variable density flows the filtered variables are Favre-filtered (denoted by a tilde), whereas the density and pressure are Reynolds-filtered (denoted by an overbar). The Favre-filtered con­ servation equations of mass, momentum and energy are defined by Eqs. (3)–(5):

∂ρ ∂ρ~ui ¼0 þ ∂t ∂xi

4. Results 4.1. Velocity field The velocity field contains important information about the mixing process and provides valuable insight for the understanding of the mixing of different temperatures in a mixing tee. The mesh study and validation is performed with the experimental measurements of case A (see Table 2). Fig. 3 shows the comparison of the averaged axial velocity ux of case A along the z-coordinate in the main pipe before the direct mixing region at the position x ¼ 0:36Dm . The black line represents the numerical results of the coarse mesh, the red dashed line symbolizes the intermediate mesh and the dotted line represents the simulation results of the fine mesh. The black squares show the experimental results. It can be observed that the simulation result of the averaged axial velocity ux from the different meshes does not differ much from the experimental results and shows a good agreement. The velocity profile is developed and fully turbulent. To sum up, it can be shown that with the mapped boundary conditions the profile of the axial velocity ux in the entrance of the main pipe can be well reproduced. Fig. 4 contains the axial velocity and the corresponding velocity fluctuations in the direct mixing region at the positions x ¼ 0:5Dm and x ¼ 1Dm . The fluid from the branch pipe is bent down by the high fluid velocity in the main pipe. Therefore, it flows along the pipe wall. As a consequence of the direct mixing region of the two pipes, a recirculation and wake area forms, where the axial velocity decreases rapidly and becomes insignificant. The strong distraction of the branch pipe jet leads to a negative axial velocity in the recirculation area at the position x ¼ 0:5Dm . It follows that the axial velocity increases in the upper part of the mixing region. The numerical results of all three meshes show a good agreement with the experimental data at the position x ¼ 0:5Dm . In addition, it is also noticeable, that in the recirculation area the axial velocity fluctuations increase strongly and show also a peak in both investigated positions. At a position x ¼ 0:5Dm the experimental peak reaches the value ux;RMS ¼ 0:52. The results of the middle and fine mesh show a higher peak value, which is in the range of 0:58 0:6. The coarse mesh can reproduce the

(3) �

∂ρ~ui ∂ρ~ui u~j ¼ þ ∂t ∂xj

�

∂p ∂ ∂~ui ∂~uj þ μ þ ∂xj ∂xj eff ∂xj ∂xi �

∂ρ~h ∂ρ~uj ~h ∂ ∂~h ¼ α þ ∂xj eff ∂xj ∂t ∂xj

�

� u 2 ∂~ μeff k δij þ ρgi 3 ∂xk

(4)

� (5)

~ represent the filtered density, the filtered pressure, ~i and h where ρ p, u the Favre-filtered velocity and the Favre-filtered enthalpy, respectively. The acceleration due to gravity is defined as gi , δij represents the Kronecker-Delta and i, j and k are the tensor indices. The turbulent mixing contributes an additive term in the momentum equation and is summarized in the effective dynamic viscosity μeff as the sum of the molecular and turbulent viscosities, μeff ¼ μ þ μt . The mo­ lecular and turbulent diffusivities are summarized in the effective μt μ diffusivity, αeff ¼ Pr þ Pr . Prt is the turbulent Prandtl number with a t

value of 0.85. The turbulent eddy viscosity is modelled using the Wall Adapting Local Eddy (WALE) subgrid-model [35]. The turbulent eddy viscosity is defined as � �3=2 Sdij Sdij 2 μt ¼ ρνt ¼ ρðCw ΔÞ (6) �5=2 � d d �5=4 Sij Sij þ Sij Sij � in which Sij ¼ 12

∂u~i ∂xj

∂u~

þ∂xji

� denotes the resolved shear strain rate tensor.

The traceless symmetric part of the square of the velocity gradient tensor is Sdij ¼ 12 ð~g2ij þ~g2ji Þ

1 ~2 3 ðδij g kk Þ

with the velocity gradient tensor ~gij ¼ ∂∂xu~ij and

Cw ¼ 0:325 represents the WALE constant.

3.2. Numerical discretization In this study, the simulations are performed in OpenFOAM (Open Source Field Operation and Manipulation). Here, the governing equa­ tions are solved by the finite volume method [36]. The linear solver geometric-algebraic multi-grid (GAMG) is used for the pressure. The PIMPLE algorithm was used for the coupling equations for momentum and mass conservation [37,38]. For the other unsteady equations, a preconditioned bi-conjugate gradient (PBiCG) solver for asymmetric lduMatrices with a diagonal incomplete-lu (DILU) pre­ conditioner is used [39]. For the time derivate terms the second order backward implicit time discretization scheme is used. For the convective energy term a second order limitedLinear scheme was selected. The convective momentum terms are discretized using a second order unbounded fixed blended linear-upwind stabilized transport (LUST) scheme with 0.25 linearUp­ wind and 0.75 linear weights [40]. For the enthalpy gradient a cellM­ DLimited Gauss linear scheme has been used. The Laplacian terms are discretized with an unbounded second order scheme. The coefficients for the spatial filter function will be left as the

Fig. 3. Comparison of time averaged axial velocity at the position x ¼ (case A). 5

0:36Dm

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

Fig. 4. Comparison of time averaged axial velocity (a) x ¼ 0:5Dm and (c) x ¼ 1Dm and time averaged axial velocity fluctuations (b) x ¼ 0:5Dm and (d) x ¼ 1Dm (case A).

peak of the axial velocity fluctuations in both investigated positions. Furthermore, the intermediate and fine mesh shows in the results of the axial velocity fluctuations differences in the position x ¼ 1Dm . Addi­ tionally, the value of the peak varies in this position and the numerical peak value of the middle and fine mesh is smaller than the experimental peak value. The results of the averaged axial velocity of the fine mesh show a good agreement with the experimental curve, mainly in the recirculation and wake area. It is noticeable, that the velocity increases in this region slowly and is becoming positive. The coarse mesh

underrates and the intermediate mesh overrates the axial velocity in the recirculation area at the position x ¼ 1Dm . Although the coarse mesh underestimates the axial velocity in some areas, it predicts the axial velocity fluctuations with a higher accuracy and is closer to the exper­ imental results. 4.2. Temperature field The temperature fluctuations in the near wall region are important

Fig. 5. Comparison of circumferential distribution of the normalized temperature fluctuations T *RMS (a) x ¼ 0:5Dm and (b) x ¼ 1Dm (case A). 6

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

because they result in the thermal stresses in the structure. Fig. 5 shows the circumferential distribution of the normalized temperature fluctu­ ations T*RMS at positions x ¼ 0:5Dm and x ¼ 1Dm in the near wall region (1 mm from the inner wall into the fluid). The temperature fluctuations are normalized with the temperature difference between the main and

area between warm and cold and in the boundary region around the cold fluid entrance. Here is the mixing layer and the impinging jet oscillates with high amplitude. Due to hydrodynamic instabilities they break down and this generates large-scale transient flow structures. The tem­ perature fluctuations intensity decreases rapidly along the axial direc­ tion in the main pipe. The mean temperature shows a stable thermal stratification in the mixing flow. The cold branch pipe fluid remains in the stratification in the lower part in the stratified flow. In the stable thermal stratification area, there are no temperature fluctuations in the mixing region. The stability of stratified mixing flow reduces the po­ tential for thermal fatigue. Fig. 7 presents the comparison of experimental and numerical results

branch pipe. The normalized mean temperature T * and the normalized large-scale RMS temperature fluctuations T*RMS are defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N N u1 X X T T 1 b T* ¼ T* ¼ T * T *RMS ¼ t ðT * T * Þ2 N i¼1 N i¼1 Tm Tb

where N represents the number of sampling points. Here one can see that the temperature fluctuations in the near wall region can also be reproduced with a good accuracy. It can be seen, the peak of the temperature fluctuations at the position x ¼ 0:5Dm is for the experimental results approximately at θ ¼ 35� and for the numerical results, for all three meshes, it is around θ ¼ 30� . It can be determined, that with a high mesh resolution the value of the peak is increasing, the peak from the experiment has the value T*RMS ¼ 0:26. In the numerical results of the coarse and intermediate mesh results a peak value which is smaller than the value from the experiment, with T*RMS ¼ 0:19 for the

of the normalized mean temperature T * and the normalized temperature fluctuations T *RMS in the mixing region. The temperature profile at the position x ¼ 0:5Dm shows that the cold fluid is bent to the pipe wall and then flows near the wall. This is the consequence of a wall jet. The results of the coarse mesh show the biggest difference to the experimental results in the area z � 0:5. Here the numerical results are higher and overpredict the experimental re­ sults. For the temperature fluctuations the numerical results of the coarse mesh agree well with the experiment. However, the value of the peak of the coarse mesh does not agree with the value of the experiment. The experimental peak is approximately at T *RMS ¼ 0:3. In case of the

coarse and T *RMS ¼ 0:23 for the intermediate mesh. However, the peak value of the fine mesh reaches the experimental peak value. For smaller angles (θ � 20� ) it is noticeable, that the normalized temperature fluc­ tuations from the intermediate and fine mesh are above the experi­ mental data and the numerical results from the coarse mesh fit well with the experiment. For higher angles (θ � 55� ) the numerical results of the fine mesh show a good agreement with the experiment. With an increasing distance from the direct mixing region the peak value of the temperature fluctuations intensity decreases. At the position x ¼ 1Dm the experimental peak of the temperature fluctuations is around at θ ¼ 40� . Here is the peak value from the experiment and also from the fine mesh at T *RMS ¼ 0:23. In case of the coarse and intermediate mesh, the intensity reaches a smaller value. This area with high temperature fluctuations is a potential threat for material damage.

coarse mesh, the peak value is at T*RMS ¼ 0:24. The numerical results of the normalized temperature for the intermediate and fine mesh agree better with the experimental results. Mainly, the fine mesh can repro­ duce the trend of the experimental results at best. The temperature fluctuations can also be reproduced well. At the position x ¼ 1Dm the results of the fine mesh shows a clear improvement from the results in z � 0:5 to the coarse and intermediate mesh. The normalized tem­ perature fluctuations show that the intermediate and fine mesh have a higher peak value as the experiment and the coarse mesh. Therefore, the coarse mesh has a better agreement with the highest temperature fluc­ tuations at this position. As already mentioned, the fine mesh has in both investigated positions a better agreement in the temperature profile. The highest fluctuations are in the mixing layer zone. It is also the area, where eddies are generated through the impinging jet from the branch pipe. The horse-shoe vortices are the mechanism which is responsible for the transport of the fluctuations.

In Fig. 6 are the numerical results of the normalized temperature T * and normalized temperature fluctuations T *RMS on circumferential dis­ tribution in the near wall region (1 mm from the inner wall into the fluid) in the direct mixing region. Here it can be seen that the jet, which is penetrating from the branch pipe, is strongly bent to the pipe wall. The high main pipe velocity and the jet type lead in the near wall region to a high temperature difference, mainly in the entrance area of the branch pipe jet. The mixing process results along the downstream and the temperature difference in the wall region becomes less and progressed. In the normalized temperature fluctuations distribution one can see that the high fluctuations in the position x ¼ 0 1Dm are in the transition

4.3. Instantaneous and mean flow field Fig. 8 (a) shows the xy-plane at 0.065 m below the coordinate center of the main pipe. Here can be clearly seen that two decisive eddies (twineddies) exist in the recirculation area. Fig. 8 (b) represents the time averaged streamlines in xz-plane. The temperature transport and tem­ perature fluctuations in the recirculation and wake area are generated

Fig. 6. Circumferential distribution of the normalized temperature (top) and normalized temperature fluctuations (bottom) along the downstream (1 mm from the main pipe wall) for case A. 7

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

Fig. 7. Comparison of time averaged normalized temperature (a) x ¼ 0:5Dm and (c) x ¼ 1Dm and time averaged normalized temperature fluctuations (b) x ¼ 0:5Dm and (d) x ¼ 1Dm (case A).

Fig. 8. Streamlines of the mean velocity field in the (a) xy-plane at the position z ¼

by large eddy structures in the flow. The counter-rotating vortex pair is similar to the observations made during the T-junction mixing in­ vestigations by Ref. [1]. The Q-criterion (second invariant of the velocity gradient tensor) where Q ¼ 0:5ð∂ui =∂xj Þð∂uj =∂xi Þ is used to investigate

0:065 m and (b) xz-plane (case A).

the relation and mechanism between temperature diffusion and large-scale transient vortex structures [41]. Fig. 9 shows time series of instantaneous isosurface at Q ¼ 4500 s12 visualized with the temperature distribution in the mixing region along

Fig. 9. Instantaneous flow characteristic visualized by the Q-criterion colored with instantaneous temperature (case A). 8

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

the downstream direction. This kind of vortex structures exist at 0:25Dm � x � 2Dm . The development of these structures can be separated into funda­ mental processes. Due to the impinging jet from the branch pipe eddies are generated (I). Thereafter, a detachment and a downward movement of the vortex structures can be seen (II). The top of the horseshoe vortices corresponds to the cold locations in the temperature distribution. After 2Dm downstream, the horseshoe vortex breaks up into small vortices and thus disappears (III). Therefore, a part of the cold fluid is segregated from the branch pipe fluid jet due to the movement of the horseshoe vortex structure and transports cold fluid to the downstream. The structures of the flow are well predicted with the present LES. In addi­ tion, other authors such as [13,26] have found similar vortices and mechanisms of coherent turbulent structures. Fig. 10 shows the prediction of the recirculation area on the down wall and the mixing layer. The acceleration area in Fig. 10 (d) starts in the connection edge of the two branches. It starts before the recirculation area and encloses it (see Fig. 10 (d)). The angle of the mixing layer is around 40� for this Tjunction with a velocity ratio of 1.46. The LES simulation performed in a T-junction by Ref. [19] has also shown that an acceleration area can be detected. In close proximity to the mixing layer the fluid temperature is equal to the average temperature of the fluid from the main and branch pipe. The horseshoe vortices form and cover the edges of the accelera­ tion. This, leads to a formation of the constrained recirculation area. The horseshoe vortex structures also form the wave-shaped temperature layer in the center of the main pipe (see Fig. 10 (a)). This leads to high temperature fluctuations in the center of the main pipe. By the transient velocity vector field the horseshoe vortex structure provides the cold fluid to the recirculation and wake area with its rotational motion as it passes this zone. Fig. 11(a–d) shows the temperature and velocity distribution at different streamwise positions downstream of the T-junction for the investigated case. The mixing process downstream of the T-junction shows that the penetrating cold fluid is heated in the downward part of the pipe and at the same time it cools this part while the upper part of the pipe remains hot. As a result, the majority of thermal mixing is in the downward part as compared to the upper pipe part.

The formation of the impinging jet leads to interaction between the two fluid streams across the entire main pipe cross-section. The mixing process results along the downstream (see Fig. 11). Fig. 12 shows the mean normalized temperature distribution in the xz-plane for the different cases (see Table 2.) The jet formation looks similar for all three cases. Hence, we can be conclude that with a constant velocity ratio and a constant momentum ratio the flow pattern type in the direct mixing region does not change. 4.4. Spectral analysis To analyze the grid sensitivity, the energy spectrum of the velocity has to be considered (see Fig. 13). It can be seen that with an increasing mesh resolution, the inertial subrange with a slope of 5/3 is increasing. The fine mesh resolves an adequate amount of the inertial subrange. In comparison with the intermediate mesh, it shows an improvement. The coarse mesh can hardly reproduce the slope in the inertial subrange. Important turbulent scales are resolved. The results of the spectral analysis at the position x ¼ 1Dm and θ ¼ 30� can be seen in Fig. 14. One significant factor to estimate High-Cycle Thermal Fatigue is a frequency analysis of temperature fluctuations. The conversion of tem­ perature fluctuations in the fluid to stress in the solid depends on the fluctuations frequency [42]. The solid wall tends to be inured to quick high frequency fluctuations affiliated with small-scale turbulent struc­ tures. However, large-scale turbulent structures and the corresponding fluctuations can lead to fluctuations in the structure wall. In Ref. [2], the critical frequencies of temperature fluctuations were probably respon­ sible for possible thermal fatigue damage in the solid wall. These critical frequencies are in the range of 0.1 and 10 Hz. The Vattenfall T-junction investigation looked in detail for the existence of a spectral peak (dominant frequency) of temperature fluctuations in the previously defined frequency range [6]. A frequency analysis for the temperature fluctuations in the near wall region (1 mm from the wall into the fluid) was performed for the different mesh resolutions. The dominant frequency of temperature fluctuations can be observed in all three mesh resolutions at 6 Hz. In the frequency analysis of the intermediate and fine mesh the spectral peak is clearly visible and the energy content is also high.

Fig. 10. Instantaneous (a, b) and mean (c, d) normalized temperature (a, c) and velocity (b, d) for case A. 9

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

Fig. 11. Instantaneous (a, b) and mean (c, d) normalized temperature (a, c) and velocity (b, d) in the cross-section at 1Dm , 2Dm , 3Dm and 4Dm (case A).

Fig. 12. Mean normalized temperature in xz-plane for (a) case A, (b) case B and (c) case C.

Fig. 13. Velocity spectrum at the centerline for different meshes at x ¼ 4Dm (case A).

Fig. 14. Comparison of power spectrum density for different meshes at the position x ¼ 1Dm and θ ¼ 30� (1 mm from the main pipe wall) for case A.

Fig. 15a) illustrates the frequency analysis for case A, B and C with the different velocities in the main and branch pipe, however, with a constant velocity ratio and momentum ratio MR (see Table 2). In

Fig. 15b) the normalized power spectral and the Strouhal ðPSD* ¼ ðPSD *um Þ =ðΔT2 *Dm Þ ​ Þ ðSt ¼ f *Db =um Þ is used instead of dimensional parameters. 10

density number

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231

Fig. 15. Comparison of power spectrum density (left) and non-dimensional power spectrum density (right) for the cases A, B and C.

The profiles of the non-dimensional power spectral density PSD* are in good agreement with the velocity variated cases. The spectral peak is at a Strouhal number of 0:2. It is demonstrated that the frequency analysis of the temperature fluctuations shows in a non-dimensional diagram (see Fig. 15 (b)) an identical profile for cases with the same momentum ratio and velocity ratio. It can be considered that the PSD* profiles can be used to account the frequency characteristics of the temperature fluctuations under the condition of wide range velocities.

Acknowledgments The presented work is funded by the German Federal Ministry of Economic Affairs and Energy (BMWi, project no. 1501575) on the basis of a decision by the German Bundestag. The authors would like to convey their sincere thanks to the High Performance Computing Center Stuttgart (HLRS) for providing the necessary computational resources to perform LES calculations. References

5. Conclusions

[1] H. Kamide, M. Igarashi, S. Kawashima, N. Kimura, K. Hayashi, Study on mixing behavior in a tee piping and numerical analyses for evaluation of thermal striping, Nucl. Eng. Des. 239 (2009) 58–67. [2] S. Chapuliot, C. Gourdin, T. Payen, J.P. Magnaud, A. Monavon, Hydro-thermalmechanical analysis of thermal fatigue in a mixing tee, Nucl. Eng. Des. 235 (2005) 575–596. [3] S. Courtin, High cycle thermal fatigue damage prediction in mixing zones of nuclear power plants. Engineering issues illustrated on the father case, Proceed. Eng. 66 (2013) 240–249. [4] M. Dahlberg, K.F. Nilsson, N. Taylor, C. Faidy, U. Wilke, S. Chapuliot, D. Kalkhof, I. Bretherton, J.M. Church, J. Solin, J. Catalano, Development of a European Procedure for Assessment of High Cycle Thermal Fatigue in Light Water Reactors. Final Report of the NESC-Thermal Fatigue Project, EUR 22763 EN, ISSN, 2007, pp. 1018–5593. [5] C. Brückner, Study of the three-dimensional flow in a T-junction using a dualscanning method for three-dimensional scanning-particle-image velocimetry, Exp. Therm. Fluid Sci. 14 (1997) 35–44. [6] B.L. Smith, J.H. Mahaffy, K. Angele, J. Westin, Report of the OECD/NEA-Vattenfall T-Junction Benchmark Exercise, 2011. [7] T. H€ ohne, Scale resolved simulations of the OECD/NEA–Vattenfall T-junction benchmark, Nucl. Eng. Des. 269 (2014) 149–154. [8] H. Ayhan, C.N. S€ okmen, CFD modeling of thermal mixing in a T-junction geometry using LES model, Nucl. Eng. Des. 253 (2012) 183–191. [9] B.L. Smith, J.H. Mahaffy, K. Angele, A CFD benchmarking exercise based on flow mixing in a T-junction, Nucl. Eng. Des. 264 (2013) 80–88. [10] M. Kuschewski, R. Kulenovic, E. Laurien, Experimental setup for the investigation of fluid–structure interactions in a T-junction, Nucl. Eng. Des. 264 (2013) 223–230. [11] P.K. Selvam, R. Kulenovic, E. Laurien, J. Kickhofel, H.-M. Prasser, Thermal mixing of flows in horizontal T-junctions with low branch velocities, Nucl. Eng. Des. 322 (2017) 32–54. [12] M. Zhou, R. Kulenovic, E. Laurien, T-junction experiment with high temperature and high pressure to investigate flow rate influence on mixing characteristics, Int. J. Heat Fluid Flow 71 (2018) 451–459. [13] P.K. Selvam, R. Kulenovic, E. Laurien, Experimental and numerical analyses on the effect of increasing inflow temperatures on the flow mixing behavior in a Tjunction, Int. J. Heat Fluid Flow 61 (2016) 323–342. [14] P.K. Selvam, R. Kulenovic, E. Laurien, Large eddy simulation on thermal mixing of fluids in a T-junction with conjugate heat transfer, Nucl. Eng. Des. 284 (2015) 238–246. [15] P. Gauder, P. Karthick Selvam, R. Kulenovic, E. Laurien, Large eddy simulation studies on the influence of turbulent inlet conditions on the flow behavior in a mixing tee, Nucl. Eng. Des. 298 (2016) 51–63. [16] S.T. Jayaraju, E.M.J. Komen, E. Baglietto, Suitability of wall-functions in large eddy simulation for thermal fatigue in a T-junction, Nucl. Eng. Des. 240 (2010) 2544–2554.

The turbulent flow in a T-junction has been investigated by using Large-Eddy Simulation with the open-source code OpenFOAM. The investigated conditions are the same as in the WATLON experiment. The hot main pipe has a Reynolds number of approximately 380000 and the cold branch pipe has a Reynolds number of around 66000. The turbulent mixing process is produced through the hot and cold streams. To ensure a sufficient resolution of the numerical grid a grid sensi­ tivity study was performed. The simulation results were compared with the experimental data and a good agreement was shown. The axial ve­ locity and its fluctuations can be reproduced well. The RMS values of the temperature fluctuations and the axial ve­ locity show that the temperature fluctuations have their high values in the mixing layer while the velocity fluctuations are very high in the recirculation area. In the near-wall region of the pipe, high fluid tem­ perature fluctuations were observed. Due to the turbulent mixing, horseshoe vortices appear which encourage the transport of the thermal fluctuations. Furthermore, a recirculation area develops in the direct mixing region. Decisive eddies are existing in the mixing region. Large-scale horseshoe vortices are visualized by the isosurface of the second invariant of the velocity gradient tensor. This kind of vortices caused high fluctuations during the mixing process. The frequency analysis of the thermal fluctuations showed a domi­ nant frequency in the HCTF range. The prominent frequency is constant at St ¼ 0:2. The presence of peak is less sensitive to the grid resolution. However, for a clear representation of a spectral peak in the frequency analysis, a high mesh resolution is required. Declaration of competing interest The authors declare that they have no conflict of interest.

11

C. Evrim and E. Laurien

International Journal of Thermal Sciences 150 (2020) 106231 [29] T. Frank, C. Lifante, H.M. Prasser, F. Menter, Simulation of turbulent and thermal mixing in T-junctions using URANS and scale-resolving turbulence models in ANSYS CFX6, Nucl. Eng. Des. 240 (9) (2010) 2313–2328. [30] D. Kl€ oren, E. Laurien, Large-eddy simulations of stratified and non-stratified Tjunction mixing flows, in: W.E. Nagel, D.H. Kr€ oner, M.M. Resch (Eds.), High Performance Computing in Science and Engineering ‘12, Springer Berlin Heidelberg, 2013, pp. 311–324. [31] T.S. Lund, X. Wu, K.D. Squires, Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations, J. Comput. Phys. 140 (1998) 233–258. [32] S. Pandey, X. Chu, E. Laurien, Investigation of in-tube cooling of carbon dioxide at supercritical pressure by means of direct numerical simulation, Int. J. Heat Mass Transf. 114 (2017) 944–957. [33] S. Pandey, X. Chu, E. Laurien, B. Weigand, Buoyancy induced turbulence modulation in pipe flow at supercritical pressure under cooling conditions, Phys. Fluids 30 (6) (2018) 65105. [34] E. Garnier, N. Adams, P. Sagaut, Large Eddy Simulation for Compressible Flows, first ed., Springer Netherlands, 2009. [35] F. Nicoud, F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor, Flow, Turbul. Combust. 62 (1999) 183–200. [36] F. Moukalled, L. Mangani, M. Darwish (Eds.), The Finite Volume Method in Computational Fluid Dynamics, Springer, Heidelberg, 2016. [37] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operatorsplitting, J. Comput. Phys. 62 (1) (1986) 40–65. [38] S.V. Satankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980. [39] M.R. Hestenes, E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand. 29 (1952) 409–436. [40] J. Martínez, F. Piscaglia, A. Montorfano, A. Onorati, S.M. Aithal, Influence of spatial discretization schemes on accuracy of explicit LES. Canonical problems to engine-like geometries, Comput. Fluids 117 (2015) 62–78. [41] J.C.R. Hunt, A.A. Wray, P. Moin, Eddies, streams, and convergence zones in turbulent flows, in: Proceedings 1988 Summer Program. Center for Turbulent Research, Stanford University, 1988, pp. 193–207. [42] N. Kasahara, H. Takasho, A. Yacumpai, Structural response function approach for evaluation of thermal striping phenomena, Nucl. Eng. Des. 212 (2002) 281–292.

[17] K. Miyoshi, A. Nakamura, Y. Utanohara, N. Takenaka, An investigation of wall temperature characteristics to evaluate thermal fatigue at a T-junction pipe, Mech. Eng. J. 1 (2014) 1–12. [18] Y. Utanohara, K. Miyoshi, A. Nakamura, Conjugate numerical simulation of wall temperature fluctuation at a T-junction pipe, Mech. Eng. J. 5 (2018) 1–23. [19] R.J.A. Howard, E. Serre, Large-eddy simulation in a mixing tee junction: high-order turbulent statistics analysis, Int. J. Heat Fluid Flow 51 (2015) 65–77. [20] M. Georgiou, M.V. Papalexandris, Direct numerical simulation of turbulent heat transfer in a T-junction, J. Fluid Mech. 845 (2018) 581–614. [21] V.A. Obabko, P.F. Fischer, J.T. Tautges, M.V. Goloviznin, A.M. Zaytsev, V. V. Chudanov, A.V. Pervichko, E.A. Aksenova, A.S. Karabasov, Large eddy simulation of thermo-hydraulic mixing in a T-junction, Nucl. React. Therm. Hydraul. Other Appl. (2012) 45–69. [22] A. Nakamura, Y. Utanohara, K. Miyoshi, N. Kasahara, A review of evaluation methods developed for numerical simulation of the temperature fluctuation contributing to thermal fatigue of a T-junction pipe, E-J. Adv. Mainten. 6 (2015) 118–130. [23] M. Hirota, E. Mohri, H. Asano, H. Goto, Experimental study on turbulent mixing process in cross-flow type T-junction, Int. J. Heat Fluid Flow 31 (5) (2010) 776–784. [24] M. Georgiou, M.V. Papalexandris, Turbulent mixing in T-junctions. The role of the temperature as an active scalar, Int. J. Heat Mass Transf. 115 (2017) 793–809. [25] T. Feng, M. Wang, P. Song, L. Liu, W. Tian, G.H. Su, S. Qiu, Numerical research on thermal mixing characteristics in a 45-degree T-junction for two-phase stratified flow during the emergency core cooling safety injection, Prog. Nucl. Energy 114 (2019) 91–104. [26] M. Tanaka, O. Hiroyuki, M. Hideaki, Thermal mixing in T-junction piping system related to high-cycle thermal fatigue in structure, J. Nucl. Sci. Technol. 47 (2010) 790–801. [27] S.B. Pope, Turbulent Flows, Cambridge University Press, Cambridge, United Kingdom, 2000. [28] S. Kuhn, O. Braillard, B. Ni�ceno, H.-M. Prasser, Computational study of conjugate heat transfer in T-junctions, Nucl. Eng. Des. 240 (6) (2010) 1548–1557.

12