Large eddy simulation of thermal mixing with conjugate heat transfer at BWR operating conditions

Nuclear Engineering and Design xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Large eddy simulation of thermal mixing with conjugate heat transfer at BWR operating conditions Mattia Bergagioa, , Wenyuan Fana, Roman Thielea, Henryk Anglarta,b ⁎

a b

AlbaNova University Center, Nuclear Engineering Division, Department of Physics, Royal Institute of Technology, 106 91 Stockholm, Sweden Institute of Heat Engineering, Warsaw University of Technology, 21/25 Nowowiejska Street, 00-665 Warsaw, Poland

ARTICLE INFO

ABSTRACT

Keywords: Turbulent mixing Conjugate heat transfer LES WALE SGS model High cycle thermal fatigue

Thermal fatigue occurs in most metals under cyclic heat loads and can threaten the structural integrity of metal parts. Detailed knowledge of these loads is of utter importance to prevent such issues. In this study, a large eddy simulation (LES) with wall-adapting local eddy viscosity (WALE) subgrid model is performed to better understand turbulent thermal mixing in an annulus with a pair of opposing cold inlets at a low axial level (z = 0.15 m ) and with a pair of opposing hot inlets at a higher axial level (z = 0.80 m ). Each inlet pair is 90° from each other in the azimuthal direction. Conjugate heat transfer between fluid and structure is accounted for. The geometry simplifies a control-rod guide tube (CRGT) in a boiling water reactor (BWR). LES results are compared with measurement data. This is one of the first times BWR conditions are met in both experiments and LES: pressure equals 7.2 MPa, while the temperature difference between hot and cold inlets reaches 216 K. LES temperatures at the fluid-structure interface are fairly correlated with their experimental equivalents, with regard to mean values, local variances, and dangerous oscillation modes in fatigue-prone areas (z = 0.65 0.67 m ). An elastic analysis of the structure is performed to evaluate stress intensities there. From them, cumulative fatigue usage factors (CUFs) are estimated and used as screening criteria in the subsequent frequency analysis of temperature time series at the fluid-structure interface. The likelihood of initiating a fatigue crack is linked to the maximum CUF, which is 3.2 × 10−5 for a simulation time of ~10 s .

1. Introduction Thermal fatigue is a common issue in metals under cyclic thermal loads, for example at walls of domains where turbulent mixing of hot and cold fluids occurs. These loads can initiate cracks or expand those already existing. Detailed knowledge of the thermal fluctuations in the material can help to accurately predict the risk of structural failure induced by thermal fatigue and to mitigate thermal fatigue over the lifetime of the structure. This knowledge can either be gained through a long and often expensive series of experiments or through the use of computational fluid dynamics (CFD). 1.1. Previous CFD work Computations and experiments similar to those in this work have recently been carried out by different research groups. Two experimental setups of a control-rod guide tube (CRGT) were used by Angele et al. (2011) to investigate thermal fluctuations at the control-rod stem for two distinct temperature differences between hot

⁎

and cold inlets. Here, this difference is termed T . The test section geometry resembles that of CRGTs in the Swedish boiling water reactors (BWRs) Forsmark 3 and Oskarshamn 3, where thermal fatigue due to turbulent mixing was discovered during routine maintenance (Angele et al., 2011). The first test section was made of plexiglass, which restricted T to about 30 K at ambient pressure. The second test section was made of steel and T was raised to 80 K. Additionally, the authors performed a CFD analysis of the test section using both unsteady Reynolds-averaged Navier-Stokes (URANS) equations with the k- SST turbulence model and scale-adaptive simulations (SAS). Computational and experimental results were shown to agree fairly well; however, spectral peaks in SAS temperatures are larger and less even than in experiments, so caution should be exercised when using these results to predict thermal fatigue at reactor conditions. Angele et al. (2011) suggested that the next step in the investigation be to perform detached eddy simulations (DES) and large eddy simulations (LES). In a DES the core flow is modeled using LES, while URANS is employed close to the wall. Similarly to Angele et al. (2011), Lillberg (2013) performed an LES

Corresponding author. E-mail addresses: [email protected] (M. Bergagio), [email protected] (W. Fan), [email protected] (H. Anglart).

https://doi.org/10.1016/j.nucengdes.2019.110361 Received 18 October 2018; Received in revised form 26 July 2019; Accepted 22 September 2019 0029-5493/ © 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Mattia Bergagio, et al., Nuclear Engineering and Design, https://doi.org/10.1016/j.nucengdes.2019.110361

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Nomenclature

cp Co D Dp g g G h Lp p Pe Pr q (r , , z ) R1 R2 Rio Roi Re s s0 s00 S t T Tf , d u V (x , y , z ) y+

Acronyms CHT conjugate heat transfer CRGT control-rod guide tube CUF cumulative usage factor DES detached eddy simulation EMD empirical mode decomposition HHT Hilbert-Huang transform HSA Hilbert spectral analysis IMF intrinsic mode function LES large eddy simulation MAE mean absolute error PISO pressure-implicit split-operator RANS Reynolds-averaged Navier-Stokes RMSE root mean squared error SGS subgrid scale SIMPLE semi-implicit method for pressure-linked equations TB location in the inner tube, at = 0.035 mm from the inner surface URANS unsteady RANS WALE wall-adapting local eddy viscosity WB location in the water region, at = 0.035 mm from the inner surface Greek symbols

t T

µ

thermal diffusivity times density distance from the wall LES filter width time step absolute temperature difference between cold and hot inlets dimensionless axial coordinate thermal conductivity dynamic viscosity density Cauchy stress tensor shear-rate tensor instantaneous frequency

Subscripts C eff H S

cold inlets effective hot inlets solid

Operators

· ·

Roman symbols c

specific heat capacity at constant pressure Courant number cumulative usage factor, or CUF inner diameter of the inlet pipes intrinsic mode function, or IMF gravitational acceleration number of IMFs extracted from a time series specific enthalpy distance between inlet and sampling plane pressure Péclet number Prandtl number heat flux cylindrical coordinate system attached to the inner tube 225°} region defined as R1 = {( , z ): 135° region defined as outer radius of the inner tube inner radius of the outer tube Reynolds number displacement at the current time displacement at the last time displacement at the time previous to the last rate-of-strain tensor time temperature inner-surface temperature instantaneous velocity cell volume Cartesian coordinate system attached to the inner tube dimensionless distance from the wall

low-pass filtered, unless otherwise specified Favre-filtered, unless otherwise specified

specific heat capacity

of turbulent thermal mixing in a CRGT. The simulation was run with the open-source toolbox OpenFOAM using a one-equation eddy viscosity LES subgrid-scale (SGS) model without wall functions for turbulence closure of the Navier-Stokes equations. Wall functions were not employed because their use would lead to inadequate results: most wall functions are developed for steady and fully developed flow, whereas the boundary layer in the area of interest quickly builds up and breaks down. Conjugate heat transfer (CHT) was included: the temperature equations were also solved in the solid domain and coupled to the fluid energy equation at the wall boundary. Acceptable agreement between experiment and simulation was found for the most important variables, such as temperature frequencies and amplitudes, at a location where RMS temperatures are high. In preparation for the experiments in Bergagio et al. (2017), Pegonen et al. (2014) performed an LES of turbulent mixing at reactor conditions in a simplified model of a BWR CRGT using the dynamic Smagorinsky model ( T = 216 K ). This simulation was performed to define the size of the mixing zone and to determine dominant

frequencies. In the most dangerous region, these frequencies were found to be about 0.1 Hz; however, null flux was prescribed at the walls, which were not modeled. Setups which are comparable to, but less complex than, the aforementioned CRGT setups are experiments on thermal mixing in T-junctions. Westin et al. (2008) investigated experiments of this kind at the Vattenfall test facility ( T = 15 K ) using LES with an adapting local eddy-viscosity (WALE) SGS model (Nicoud and Ducros, 1999). Westin et al. (2008) did not only focus on inlet conditions, but also on how different turbulence modeling approaches and mesh resolutions influence LES results. Three inlet velocity boundary conditions were considered: generation of synthetic turbulence through a 2D vortex method; no perturbation; and isotropic inlet velocity from data files. Westin et al. (2008) showed that LES results are not very sensitive to the inlet conditions tested. A DES of the same case was then conducted. Its results were compared with those from an LES featuring the same inlet velocity boundary conditions. It was shown that both modeling strategies deliver acceptable mean temperatures at the left and right 2

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

sides of the pipe wall downstream of the T-junction, but only the LES could deliver accurate RMS temperatures there, because it is less dissipative than the DES. Specifically, RMS values from DES were reported to be 150-200% the experimental RMS values at the same coordinates, between 4 and 10 diameters downstream of the T-junction. Overall, although the coarser mesh could capture large-scale unsteady flow in the mixing region, the finer mesh yielded better results, which is expected because a finer mesh resolves more turbulent scales. The same wall-resolved WALE model was successfully employed by Jayaraju et al. (2010) to calculate flow and heat transfer in the Vattenfall T-junction. Here, the wall-modeled approach with wall functions was found to underestimate gradients of RMS velocities and temperatures close to the walls, as well as RMS heat fluxes along the streamwise direction. These shortcomings could be addressed by adopting a wallresolved approach. Jayaraju et al. (2010) were among the 29 research groups involved in the OECD/NEA blind benchmark exercise (Smith et al., 2013) to assess how well CFD can predict turbulent mixing in T-junctions. CFD results were compared against the Vattenfall experimental data ( T = 17 K ). CHT was neglected. 19 groups chose LES, 3 chose DES, and the rest chose RANS. LES held the first 7 or 9 ranks in terms of agreement with experimental velocities, temperatures, and peaks in smoothed Fourier velocity spectra. It was suggested that the number of cells be high enough to accurately predict velocity; however, a high number of cells does not imply correct temperature predictions, as their distribution is also critical. LES dominant frequencies in both velocity and temperature spectra were found to be subject to user input: fluctuation modes might be susceptible to boundary conditions and nodalization. SAS performed much worse than LES, especially in terms of velocity RMS values and spectra. Timperi (2014) investigated different inlet conditions and the influence of CHT on LES. Here, results from the above simulations were compared with the Vattenfall experimental data. The Smagorinsky SGS model was adopted. Two different inlet conditions were investigated: steady profiles and synthetic turbulent inflow data based on a 2D vortex method. Both inlet conditions led to good agreement between experimental and numerical results; however, analogously to Westin et al. (2008), the case with the turbulent inlet showed slightly worse agreement in some areas. CHT was confirmed to have a strong influence on the temperature fluctuations at the wall surface and to be important for accurate results, because of the thermal inertia of the wall. Karthick Selvam et al. (2017) performed LES of thermal mixing in Tjunctions with low branch velocities using the WALE model. LES results were compared with measurement data from the T-junction test facility at the University of Stuttgart. Here, two cases were considered: (1) T = 65 K ; and (2) T = 143 K . The average mesh size was determined according to Addad et al. (2008); that is, by conducting a precursor RANS on the same geometry to determine the turbulence length scales in the flow. Thus, a highly optimized LES mesh could be obtained. Additionally, CHT was applied at the wall. Overall, experimental and numerical temperature time series were shown to agree reasonably well with each other, in terms of means and RMS values. However, RMS values were found to be over- or underestimated at certain locations. In the areas examined, Karthick Selvam et al. (2017) could not find dominant frequencies between 0.1 and 10 Hz; i.e., in the frequency range which is deemed more relevant to thermal fatigue in T-junctions (Chapuliot et al., 2005). Overall, LES with WALE SGS model appear to predict thermal mixing better than RANS. A wall-resolved approach should be preferred to a wall-modelled one. LES results do not seem to be strongly affected by inlet conditions, but synthetic turbulent inflow data might yield worse results if the mesh resolution is too low. The number of cells does not have to be extremely high: the optimal mesh size should be computed from turbulence lengthscales, which can be estimated based on precursor RANS. CHT has to be included to accurately model the interaction between fluid and surrounding walls. The simulation time

should be longer than 10 s. To the authors’ knowledge, LES of thermal mixing with CHT and T > 200 K (i.e., at BWR conditions) have not been performed before. These simulations need validation against experimental data, especially near the wall. Furthermore, temperature fluctuation modes leading to thermal fatigue have not yet been properly identified. This research tries to address such issues. 1.2. Previous work on CFD-FEA coupling So far high-cycle thermal fatigue has been investigated from various perspectives. Here, our interest lies in coupling CFD and FEA (finite element analysis). On this matter, the Network for Evaluation of Structural Components (NESC) outlined a four-level approach to quantify thermal loads. In Level 3, NESC recommended assessing fatigue by examining the full local load spectra along with fatigue S-N curves (Dahlberg et al., 2007). Other researchers have tried to assess fatigue damage starting from CFD results. In the studies summarized here, unless otherwise specified, coupled CFD-FEA simulations of turbulent, thermal mixing in T-junctions are performed; CHT is accounted for; linear elastic analyses are conducted under CFD thermal loads; and cumulative usage factors (CUFs) are evaluated according to the rainflow-counting procedure together with Miner’s linear damage rule. Chapuliot et al. (2005) performed a very large eddy simulation (VLES) of flow mixing in a T-junction with four upstream bends. Water temperatures in the hot and cold legs differed by T = 160 K . The Smagorinsky SGS model was applied. Adiabaticity was assumed at the water-steel interface. The VLES simulation time reached 20 s. The elastic analysis was conducted using austenitic steel properties at room temperature, thermal loads from the last 10 s of the CFD simulation, and a constant heat transfer coefficient. The CFD and FEA meshes were conformal. Crack propagation was studied based on Paris’ law. Jhung (2013) performed a URANS of flow mixing in a T-junction with T = 117 K . The SST model was applied. CFD and FEA meshes were conformal. Their respective simulation times were 300 s and 200 s. CUFs were deemed very low in the seven locations examined. Kim et al. (2013) performed a DES of flow mixing in a T-junction. The boundary conditions were similar to those in Jhung (2013). In the one-way separate analysis, stresses were only calculated in the downstream part of the T-junction, considering austenitic steel properties at an average temperature. CFD and FEA simulation times equaled 15 s and 5 s, respectively. To predict significant CUFs from the one-way analysis, a mean stress correction was included in one of the S-N curves under scrutiny through the modified Goodman equation. Hurrell et al. (2015) performed a URANS of flow mixing in a Tjunction with internal sleeve. Here, T = 200 K . An LES was performed on a part of the domain using an adiabatic heat transfer model. The WALE model was applied. Temperature fluctuations on the liquid-solid boundary and area-averaged heat transfer coefficients were interpolated from the CFD mesh to the FEA one. The simulation time amounted to 27 s. Similarly to Kim et al. (2013), stress amplitudes were estimated from the modified Goodman equation. Fatigue damage was also assessed through in-phase sinusoidal cycling and modified maximum stress approaches, which were regarded as more conservative than the rainflow approach. Zhang and Lu (2016) performed LES simulations of flow mixing in a T-junction with T = 80 K . The Smagorinsky model was applied. Although different simulation times were investigated (namely, 5, 10, and 20 s), the elastic analysis was conducted using thermal loads from 5 s of the CFD simulation. The FEA mesh was made coarser than the CFD one. Similarly to Kim et al. (2013) and Hurrell et al. (2015), the modified Goodman equation was applied. Von Mises equivalent stresses reduced the multiaxial stress state to a uniaxial one. Although the uniaxial S-N curve was extrapolated below the fatigue limit, CUFs were found to be very low. 3

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Wilson et al. (2016) performed LES simulations of flow mixing in a T-junction for two inlet temperature differences: T = 195 K (Transient 1) and T = 155 K (Transient 2). Although the simulation time reached 20 s, the first 2 s were neglected. CFD and FEA meshes featured the same element size, at least in the regions of interest. In order to improve statistics, the whole stress history was split into four segments of 2.5 s each. A simplified rainflow-counting method for a repeating history was then applied to each segment. Fatigue contour plots were generated after blending the fatigue analysis results at each node with those from several neighbor nodes for better confidence. The final CUF was much higher in Transient 1. Under BWR conditions, Reynolds number Re and molecular Prandtl number Pr differ from those in experiments from Angele et al. (2011). Lower viscosities at higher temperatures can lead to different penetration depths and therefore to dissimilar results in the mixing zone. As evidenced in Karthick Selvam et al. (2017), so far turbulent thermal mixing at T > 100 K has been sparingly investigated at test facilities. We circumvented this issue by conducting one of the first experiments on turbulent thermal mixing under BWR conditions (Bergagio and Anglart, 2017; Bergagio et al., 2017) before starting the current LES. In the present work, an LES is performed to predict temperature fluctuations and fatigue damage at the wall in a simplified model of a CRGT under BWR conditions. CHT between fluid and adjoining walls is included. LES and experimental data are compared in terms of means, variances, and instantaneous frequencies in critical areas. Furthermore, transient stresses from an elastic analysis of the inner tube, which echoes the control-rod stem, are computed. Fatigue damage factors are also estimated.

|u|2

resolved kinetic energy per unit mass; that is, K = 2 . Term p is net glected to avoid divergence when density changes with temperature. Heat flux q eff is written as

q eff =

Mechanical source ·( · u) , which describes the work rate done by viscous forces on the fluid, is not added to the energy equation. In the solid regions, Eq. (5) reduces to

( S hS )

with S = and hS : =hS (TS ) . Concerning CHT, boundary temperature Tf , d is estimated by linearly approximating the gradients in Eqs. (6) and (8):

(T1 ) 1

(Tf , d

T2),

(9)

eff

p +

·

eff

+ g,

SGS

( K) + t

· q eff + u·g,

PrSGS

SGS

(11)

is simply given by the Reynolds analogy

,

(12)

(13)

kSGS ,

kSGS = Cw2

(4)

·( uK ) =

=

µSGS

SGS

where SGS kinetic energy kSGS is obtained as

The Favre-filtered energy equation is given by

·( uh ) +

+

µSGS = Ck

(3)

1 ( u + ( u)T ). 2

(T ) c p (T )

SGS

where the turbulent Prandtl number, PrSGS , is set to 0.85. µSGS in Eqs. (10) and (12) is calculated using the WALE model (Nicoud and Ducros, 1999), which reproduces the turbulence damping caused by the no-slip condition. Specifically, it captures the 3 behavior of µSGS near the wall. It also guarantees zero µSGS in pure shear regions and isotropic contractions/expansions, without damping functions or dynamic procedures (Nicoud et al., 2011). More precisely, it is the simplest model that automatically ensures no µSGS in laminar shear flows, which plays a key role in modeling transitional flows (Menter, 2015). Furthermore, the WALE model is less sensitive to the model constant value than the Smagorinsky one (STAR-CCM+, 2018). In the WALE model, µSGS is given by

(2)

1 ( ·u) I , 3

=

from Eq. (5).

(1)

u) =

(10)

from Eq. (2) and an SGS thermal diffusivity

where the rate-of-strain tensor S is

( h) + t

2

µeff = µ (T ) + µSGS

p + tr( SGS) , p being where modified pressure p is defined as p the low-pass filtered pressure, and g is the gravitational acceleration. The effective stress tensor eff is

S ( u) =

S (T2 )

Tf , d)

Through the filtering process we obtain an SGS eddy viscosity µSGS

1 3

= 2 µeff S

(T1

where represents the distance from the wall, while indices 1 and 2 denote the first cell centers off the wall in the fluid and solid regions, respectively. Many methods can be implemented to estimate Tf , d – see, e.g., Tuominen (2015). While for some applications, such as supercritical fluids, Eq. (9) does not provide accurate results, a preliminary investigation shows that, in our case, Tf , d from Eq. (9) is comparable to that from Tuominen (2015).

where is the low-pass filtered density and u is the Favre-filtered velocity. The Favre-filtered momentum equation is given by

eff

(8)

hS ,

S

2.2. SGS model

·( u) = 0,

·( u

(7)

S (TS ) cS (TS )

All equations are discretized and solved using OpenFOAM 4.1. A transient solver for compressible flows is specified, which computes CHT between fluid and solid regions; namely, chtMultiRegionFoam. Previous URANS underestimated temperature changes (GallegoMarcos, 2013; Haces Manzano, 2013), so an LES approach is adopted here to validate computational results against experimental data. Unlike URANS, which apply Reynolds averaging to the Navier-Stokes equations, LES filter them. The Favre-filtered mass continuity equation is given by

( u) + t

·q S ,

where

qS =

2.1. Governing equations

t

=

t

2. Methods

+

(6)

h.

eff

2

Ck

( Sd : Sd)3

( (S : S)

5 2

)

5 2

+ ( Sd : Sd) 4

, (14)

which makes use of both the rate-of-strain tensor S and the deviatoric symmetric part of the square of the velocity gradient tensor Sd

(5)

Sd = B

where h indicates the Favre-filtered enthalpy, while K denotes the 4

1 tr(B) I , 3

(15)

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

where B = 2 (A + AT ), A being equal to A = u · u . In Eq. (14), model constants Cw and Ck equal 0.5 and 0.094, respectively. Nie and Xu (2019) and Renze and Akermann (2019) have recently selected the same Ck . In the framework of LES with implicit 1 filtering, filter width is given by = V 3 , V being the cell volume; thus, it is appropriate for isotropic hexahedral meshes. 1

Table 1 Main mesh quality parameters.

2.3. Numerical methods Pressure-velocity coupling is achieved by means of the PIMPLE (merged PISO-SIMPLE) algorithm (OpenFOAM, 2017). PIMPLE ensures stability even for time steps exceeding the Courant-Friedrichs-Lewy 2 , which implies t = 8 × 10 6 s . Time integration is limit: here, Co performed using a Crank-Nicolson scheme. Divergence terms are discretized using the limitedLinear and limitedLinearV schemes in scalar and vector transport equations, respectively. In both cases, a Sweby limiter is applied to a linear (i.e., central differencing) scheme (Sweby, 1984) in order to fully conform to the total variation diminishing (TVD) criterion: pure central differencing triggers oscillations in the solutions to convective-diffusion equations if cell Pe > 2 (Patankar, 1980); moreover, Timperi (2014) stressed that pure central differencing yields unphysical oscillations of the velocity field at the inlets of a mixing Tjunction with all hexahedral cells. Gradient terms are calculated by linear interpolation, with their components equally clipped for better stability. An explicit non-orthogonal correction improves the discretization of surface-normal gradient terms. A preconditioned conjugate gradient (PCG) linear solver with a diagonal incomplete-Cholesky (DIC) preconditioner is specified for the pressure equation. The final residual at each time step is set to 1 × 10−9. A stabilized preconditioned bi-conjugate gradient (PBiCGStab) solver with a diagonal incomplete-LU (DILU) preconditioner is selected for the other scalar equations and the momentum equation. The final residuals at each time step are set to 1 × 10−11.

Region

Maximum aspect ratio

Maximum nonorthogonality (°)

Average nonorthogonality (°)

Maximum skewness

Water Inner tube Outer tube

30.63 25.63 18.16

57.12 57.73 46.08

6.28 9.03 5.39

0.93 2.43 1.02

3

LRM =

k2

(17)

and (18) RM

=

15

kµ

,

(18)

respectively (Addad et al., 2008), where k is the turbulence kinetic energy, while is its rate of dissipation. This RANS uses a standard k turbulence model for compressible flows, with a rapid distortion theory-based compression term and wall functions. A multi-block structured mesh is generated by ICEM CFD 18.0. The mesh size is set based on Eq. (16). Fig. 3 depicts the mesh near the hot inlets. Water cells smaller than 1 × 10−11 m3 can be mostly found in proximity to the inner tube, inlets, and outlets. The water region contains 17.32 × 106 cells, while the inner and outer tube regions contain 2.61 × 106 and 4.70 × 106 cells, respectively. By introducing O-grids at the inlets and outlets, the fluid region inside these tubes as well as intersections between these tubes and the outer tube are discretized by high-quality hexahedral cells – see Table 1. To smoothly reach convergence and improve stability, it is advisable to keep the maximum internal face skewness and non-orthogonality below 4 and 65°, respectively (OpenFOAM, 2017). The mesh is conformal between the fluid and solid regions. No wall functions are specified. In the mixing zone, the thickness of the first cell off the inner tube is 3.46 × 10−5 m; consequently, y+ 5 on the inner surface (that is, on the water-inner tube interface). More precisely, as shown in Fig. 1, the time-averaged y+ ranges between 0.09 and 9.40 at 0.60 m z 0.73 m .

2.4. Properties Viscosity µ , density , specific heat capacities cp and cS , and conductivities and S change with temperature. Water properties are temperature-dependent polynomial fits of the IAPWS-IF97 data (Wagner et al., 2000), while conductivity and specific heat of 316LN stainless steel are temperature-dependent polynomials taken from the MPDB software (JAHM Software Inc, 2015). Steel density is assumed constant with temperature. 2.5. Space domain and mesh The meshed regions – i.e., water, inner and outer tubes – are shown in Fig. 2. Here, the axis of the hot inlets is located at the same azimuthal angle as that of the outlets, whereas the axis of the cold inlets is positioned at 90° of azimuth from that of the outlets. All four inlet pipes have inner and outer diameters of 7.5 mm and 16 mm, respectively. The two outlet pipes have inner and outer diameters of 14 mm and 22.5 mm, respectively. Inlets and outlets are modeled as 30 mm long, straight, and perpendicular to z. According to common practice, an LES should resolve about 80% of the turbulent energy spectrum of a flow (Pope, 2000), whereas the remaining 20% is modeled using an appropriate SGS model. Therefore, the LES mesh is based on integral length scale LRM and Taylor length scale RM . In accordance with Addad et al. (2008), the cell size is determined using these two length scales:

= max

RM ,

LRM . 10

(16)

These two length scales are calculated a priori using a RANS for the same geometry and experimental boundary conditions; namely, the integral length and the Taylor length scale are approximated by Eq. (17),

Fig. 1. Time-averaged y+ on the inner surface, at 0.60 m time t0 = 19.2 s . 5

z

0.73 m . Start

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 2. Views of the meshed regions. Gray (inner region): inner tube. Red (outer region): outer tube. Green: water. The outer tube is 1.08 m long. It extends from z = 0.036 m to z = 1.116 m . Plane 0 (that is, the midplane of the cold inlets) at z = 0.15 m . Plane 1 (that is, the midplane of the hot inlets) at z = 0.8 m . Plane 2 (that is, the midplane of the outlets) at z = 1 m . Rii = 12.5 mm . Rio = 17.5 mm . R oi = 40 mm . R oo = 50 mm . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

This range of axial levels is inspected because LES and experimental temperatures are compared at 0.63 m z 0.72 m . The cell thickness gradually increases in the water layer adjacent to the wall. A growth rate of 1.15 is applied to the first two cells. In the same zone, the thickness of the first cell off the outer tube is 6.92 × 10−5 m; consequently, y+ 6 on the outer surface (that is, on the water-outer tube interface). The cell thickness gradually increases in the water layer adjacent to the wall. A growth rate of 1.15 is applied to the first two cells. In the mixing zone, the axial resolution reaches 0.8 mm. It is smaller than the axial distance between thermocouple tips on the left thermocouple disc, which equals 4.24 mm (Bergagio and Anglart, 2017).

not be kept constant throughout the experiments, nominal inlet conditions differ from the experimental values. LES inlet conditions match the experimental data. Remarkably, the above nominal values ensure azimuthal symmetry in the LES results, which is, however, broken when the actual boundary conditions are examined. The velocity inlet boundary condition can have a large impact on the accuracy of the LES results; therefore, it should be chosen with caution. According to Tabor and Baba-Ahmadi (2010) and Dhamankar et al. (2015), who investigated several velocity inlet conditions, the generated velocity field should represent the fully developed turbulent flow, be realistic – by, e.g., ensuring that the stochastically-varying velocity components conform to a given energy spectrum – and free from spatial discretization effects; additionally, it should not depend on the grid type or add any unphysical features, such as artificial periodicities. At the same time, the computational overhead for computing the stochastically-varying velocity components should be very small. The geometry in this study allows the generation of inflow data by using an internal-mapping mode, a method exploiting recycling and rescaling: the velocity field is sampled on a plane downstream of the inlet and mapped back to the inlet. Mapped velocities are scaled to average to the target mean inlet velocities in Table 2, which match their experimental analogues over a certain time interval. Specifically, each

2.6. Boundary and initial conditions In the present simulation, we apply the same boundary conditions as in Case 1 from our experimental matrix (Bergagio and Anglart, 2017). Ideally, a water mass flow rate of 0.4 kg s−1 enters each hot inlet at 549 K, while a water mass flow rate of 0.035 kg s−1 enters each cold inlet at 333 K. However, the actual boundary conditions differ from these nominal values, as emphasized in Table 2. Nominal values were selected before the experiments began. As boundary conditions could 6

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

condition is assumed. The outlet velocity gradient is set to zero, but a user-specified velocity value is imposed if reverse flow occurs. The water-steel interfaces need special treatment. Thermal fatigue triggered by fluid temperature oscillations is primarily a CHT problem; consequently, it depends on the specific combination of fluid and wall material (Tinoco, 2013). CHT is implemented at the water-steel interfaces as per Eq. (9). This method was validated in the SAFIR2010 project (Peltola et al., 2011). The innermost and outermost walls as well as top and bottom are assumed to be adiabatic. The outlets are treated likewise, but a user-specified temperature value is imposed if reverse flow occurs. A pressure of 7.2 MPa is prescribed at the outlets. At the inlets, top, bottom, and water-steel interfaces, the pressure gradient is not set to zero; instead, it is adjusted according to the velocity boundary condition. Thus, gravity is also accounted for. simulation (see Section 2.5) is carried out on First, the a priori k a coarse mesh, to find the optimal mesh size. Then, a simplified LES is performed on an intermediate mesh (Thiele, 2015, Bergagio et al., 2017), the initial temperature and velocities being mapped from the k simulation. Finally, the LES at hand is conducted on the optimal mesh. Although the initial temperature and velocities are mostly mapped from the simplified LES, temperatures in the mixing region are initialized by spatially interpolating time-averaged experimental temperatures. 2.7. Linear elastic analysis and fatigue damage assessment A finite-element analysis is performed on the inner tube to estimate both thermal stresses caused by mixing and fatigue damage. To this purpose, the Cauchy momentum equation is examined 2s S

t2

=

· + 2s

S g,

(19)

2 s0 + s00) . ( t )2

where 2 Linear elasticity is considered; accordingly, t stresses are computed as

Fig. 3. Views of the LES mesh. Gray (inner region): inner tube. Red (outer region): outer tube. Green: water. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

=

S(

(s

· s) I + 2µS (s)

(3

S

+ 2 µS )

S (TS

Tref ) I ,

(20)

where S and µS represent the first and second Lamé parameters, S is the mean coefficient of thermal expansion over temperature range (Tref , TS ) , and strain tensor (s) is defined in the same way as S (s) (see Eq. (4)). As in Section 2.4, properties of 316LN stainless steel are polynomials in temperature. Their expressions are taken from the MPDB software (JAHM Software Inc, 2015). The top surface is fixed in the axial direction (i.e., sz = 0 ), while displacement is prevented at the bottom surface (i.e., s = 0 ). No displacement is assumed as initial condition; in other terms, s0 = s00 = 0 at start time t 0 . The weak form of Eq. (19) is derived. Its numerical solution is calculated using the opensource finite-element package FEniCS (Logg et al., 2012; Alnæs et al., 2015), version 2017.2. LES temperatures are entered in Eq. (20), while LES pressure is set as boundary condition on the inner surface. The FEA mesh is coarser than the LES one: given that one of the primary objectives of the present research is to develop a methodology for estimating high-cycle fatigue damage from LES data, in future work the FEA mesh will match the LES one or the accuracy, boundedness, and

inlet pipe is 30 mm long – as already stated in Section 2.5 –, while its sampling plane is located 25.5 mm downstream along the pipe axis; that Lp is, at D = 3.4 , where Lp is the distance between inlet and sampling p

plane, while Dp is the pipe diameter. Internal mapping is recognized to have a very low computational overhead and produce high-quality turbulence (Tabor and Baba-Ahmadi, 2010, Dhamankar et al., 2015). Despite some common traits, mapping methods based on concurrent simulations, whose results are mapped directly into the main simulation, are not of interest here: these methods require additional domains, so their adoption would increase the overhead costs of calculating, storing, and reading velocities back. Westin et al. (2008),Odemark et al. (2009), and Kuczaj and Komen (2010) remarked that turbulence is mainly generated by the mixing of fluid streams; hence, similarly to Gauder (2016) and Karthick Selvam et al. (2017), the fluid flow at the inlets is not perturbed. At water-steel interfaces, top, and bottom, a no-slip boundary

Table 2 Inlet boundary conditions. u and T denote the mean velocity and temperature over each inlet cross section, respectively. The coordinate system is shown in Fig. 2. Nominal values are reported for completeness. Inlet

Hot inlet at x = 0.08 m Hot inlet at x = + 0.08 m Cold inlet at y = 0.08 m Cold inlet at y = + 0.08 m

u (ms−1) (12.21, 0, 0) ( 12.33, 0, 0) (0, 0.83, 0) (0, 0.80, 0)

Nominal u (ms−1)

(11.93, 0, 0) ( 11.93, 0, 0) (0, 0.80, 0) (0, 0.80, 0)

Nominal

T (K) 549.74 549.47 334.16 334.15

7

T (K)

Re (–)

Nominal Re (–)

Pr (–)

Nominal Pr (–)

549 549 333 333

728687 735550 13263 12892

711367 711367 12696 12696

0.85 0.85 2.93 2.93

0.85 0.85 2.98 2.98

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

conservation properties of different interpolation schemes will be tested. In order to determine stress amplitudes and number of cycles, a rainflow-counting algorithm (DTU Wind Energy, 2018) is applied to each nodal stress time history in the mixing region. Concerning the identification of extreme stress conditions (ASME, 2015b), the signed stress intensity from Green and Ferrari (2016) alleviates some shortcomings of the stress intensity based approach from the ASME Code; for example, it allows to determine whether extreme stress states are dominated by tension or compression without introducing spurious extreme stress states. The above method complies with ASME (2015b); thus, it is based on the maximum shear stress criterion. Using this method, changes in stress conditions are plotted on an equivalent uniaxial curve. The number of cycles to failure is determined from a uniaxial S-N curve for austenitic steel in air (ASME, 2015a). The elastic modulus correction from Section NB-3222.4 of ASME (2015b) is also evaluated, to compensate for the difference between the elastic modulus given on the S-N curve and that used in our elastic analysis. Similarly to Wilson et al. (2016), a CUF is computed based on Miner’s rule (linear cumulative damage hypothesis) for each node in the mixing region. Damage at stress amplitudes lower than the fatigue limit is included by extrapolating the S-N curve below the fatigue limit (Schijve, 2009). A time-step independence study is performed (see Section 4). FEA stresses and the reduced stress history are verified against results in Radu et al. (2008) and Green and Ferrari (2016), respectively.

of the experimental matrix (Bergagio and Anglart, 2017). Fig. 6 displays the locations at which LES temperatures are sampled. Although measurement locations in the inner tube, on the inner surface, and in the water region are radially aligned to determine the radial heat flux distribution, this analysis is postponed to future studies. Experimental inner-surface temperatures were measured at the same azimuthal coordinates. Fig. 7 shows streamwise velocity profiles at the centerlines of two cross sections along the inlet pipes. Considering that profiles at Lp Lp = 1.33 tally with those at D = 3.4 , Figs. 7b and d confirm that inD p

p

ternal mapping generates fully developed flow. As information on the experimental velocity profile is unavailable, the empirical turbulent velocity profile

ux ux , max

(

= 1

y

2D

p

)

1 n

is plotted for comparison. Here,

n= , Reynolds number Re being based on the mean velocity (Qian et al., 2015). The streamwise velocity profiles created by the internal mapping technique at the hot inlets are consistent with each Lp Lp other: the profiles at D = 1.33 match those at D = 3.4 . The same holds 3.45 Re 0.07

p

p

for the cold inlets. Moreover, these profiles seem to be compatible with physical velocity profiles. However, Figs. 7c and e highlight that streamwise velocity profiles at the cold inlets resemble laminar profiles. This issue is being investigated. It might be due to relaminarization, which can occur if initial conditions are not appropriate (Pamiès et al., 2009). A self-sustaining turbulent flow can be established by adding enough turbulent fluctuations to the initial conditions (Dhamankar et al., 2015). However, as stressed in Section 2.6, inlet conditions are not expected to impact on the results of LES of thermal mixing. Moreover, the cold inlets are at z = 0.15 m , far from the mixing region; therefore, their influence on mixing is expected to be low. Fig. 8 illustrates LES temperatures in the fluid at z = 0.67 m (see

3. Experimental setup The results of the current LES are compared with experiments performed in the HWAT (High-pressure WAter Test) loop at the Royal Institute of Technology in Stockholm, Sweden (Bergagio and Anglart, 2017). Of the ten runs in the experimental matrix, only one is examined here, given that an LES requires long computational times. The HWAT test section (see Fig. 4) resembles the CRGTs in Swedish BWRs Forsmark 3 and Oskarshamn 3: although it simplifies their geometry, accurate and repeatable measurements are taken in this test section without marring phenomena occurring in the above CRGTs. These experiments are performed under reactor conditions ( T = 216 K , p = 7.2 MPa ). Design and dimensions of the HWAT test section correspond to those in Fig. 2. Inner and mean outer diameters of inlets and outlets match those in Section 2.5. Top and bottom walls, which are omitted from Fig. 2, bring the overall length of the HWAT test section to 1.15 m. Temperatures at the interface between water and inner tube are termed inner-surface temperatures. They are recorded at 1000 Hz over 120 s and low-pass filtered at 4 Hz. Their uncertainty equals 1.58 K , as detailed in Bergagio and Anglart (2017), where the placement of the thermocouples sampling these temperatures is also described. Fig. 5 shows that the above thermocouples are inserted into holes with a diameter of 0.7 mm, so that their tips lie at r = Rio . The thermocouples have a diameter of 0.5 mm. The holes are drilled into so-called thermocouple discs, which are mounted into the inner-tube wall by welding. 40 measurement locations are sampled by translating the inner tube and rotating it around its axis. Boundary conditions, such as inlet mass flow rates, inlet temperatures, and pressure, are recorded at 1 Hz. 4. Results The duration of the LES transient is ~13 s , from start time t 0 = 16.3 s to end time tend = 29.2 s . The simulation required 130 days of computational time on a high-performance computing cluster. When postprocessing the LES results, the first 2.9 s is usually neglected, so start time t 0 = 19.2 s is mostly selected. The following section shows some of the LES results and compares them with the data obtained from Case 1

Fig. 4. The HWAT test section. 8

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 5. Longitudinal section of an inner-tube thermocouple at z = constant . 1 (green): thermocouple disc; 2 (yellow): thermocouple; 3 (gray): casing; 4 (red): solder. The thermocouple is soldered to the disc. L 4.5 mm . From Bergagio and Anglart (2017). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Locations where LES temperatures are sampled. Gray (inner region): inner tube. Red (outer region): outer tube. Green: water. Pink dot (TB): In the inner tube, at = 0.035 mm from the inner surface. Black dot: On the inner surface; that is, at Rio . Orange dot (WB): In the water region, at = 0.035 mm from the inner surface. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2) and different times – a snapshot is taken every ~1.4 s . During this time the temperature field in this cross section changes dramatically, often cycling between 400 K and 500 K throughout the water domain. Large, unstable hot structures penetrate the cross section at different times and locations; for example, at 24.98 s, a large hot structure reaches the cross section at ~135°, only to be replaced by colder water 1.45 s later, when hot patches border the inner tube between 180° and 270°. It can be observed that, on average, inner-surface temperatures at 90° are higher than at 360°, although the hot inlets are located at 180° and 360°. This incongruity is discussed later. Full-size snapshots are available in folder U_and_T (Bergagio, 2019). Fig. 9 shows the LES velocity magnitude at the same axial level and times as in Fig. 8. Full-size snapshots are available in folder U_and_T (Bergagio, 2019). In Fig. 9, the velocity magnitude is higher in hot areas, especially if water is unmixed or poorly mixed. Consequently, it can be inferred that, similarly to Fig. 8, the velocity magnitude at 90° and 270° is higher than at 180° and 360°. The largest in-plane velocity components are often observed in the shear layer between hot and cold areas. Even when mixing is scarce, at this level (i.e., z = 0.67 m ) the highest velocity magnitude of the hot streams is significantly decreased, to approximately one tenth that at the hot inlets (i.e., z = 0.80 m ). Despite that, as the density difference between hot and cold streams is high and the cold stream flowing upwards is slower and heavier than the hot streams, velocity redistribution seems sporadic. Vortices are shown in Figs. 9 and 10. Despite their chaotic nature and different sizes, it can be noted that most vortices appear to initiate or enhance mixing, since Figs. 8 and 9 suggest a general absence of vortices in zones of extreme temperatures. As mentioned before, mean WB and inner-surface temperatures at 90° and 270° are higher than at 180° and 360°. As Fig. 10 and figures in

Fig. 7. Time-averaged streamwise velocity profile downstream of the inlets.

folder Q_and_Uz (Bergagio, 2019) show, this increase in mean temperatures could be due to two factors: (1) vortices traveling downwards at 90° and 270° in the near-wall region are hotter and propagate farther down than those at 180° and 360°; (2) colder, often elongated vortices travel upwards at 90° and 270° from as low as 0.60 m. These elongated vortices determine more dissipation of turbulent kinetic energy into heat than the 2D structures characterizing the flow at 180° and 360° near the bottom of the contour plot. In Fig. 10, vortices are identified as isocontours of Q = 0.5 s 2 , Q being the second invariant of the velocity gradient tensor (OpenFOAM, 2017); equivalently,

Q= 9

1 (( · u)2 + 2

:

S : S ),

(21)

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 9. Velocity magnitude and in-plane velocity vectors at different times. z = 0.67 m .

Fig. 8. Temperature at different times. Water region. z = 0.67 m .

layer, on the inner surface, and in the inner tube span about 80 K, 35 K, and 37 K, respectively; in Fig. 11b, a similar behavior is observed, but ranges are now sensibly reduced, with temperatures in the water layer, on the inner surface, and in the inner tube varying over about 22 K, 7 K, and 6 K, respectively. The shrunken ranges hint that the temperature range strongly depends on the measurement region and axial level. Fig. 11 also proves that the radial heat flux changes sign sporadically, when the inner tube is hotter than water. Altogether, this figure evidences that the time series under study are highly intermittent; therefore, frequency analysis techniques should handle them properly. Broadly speaking, temperature ranges and frequencies change with azimuthal coordinate, as highlighted in Fig. 12. Here, computational and experimental inner-surface temperatures are presented, at different angles, 0.13 m below the midplane of the hot inlets. While LES temperature ranges and means match their experimental analogues at

where the rate-of-strain tensor S (u) is defined in Eq. (4), while the 1 vorticity tensor is expressed as (u) = 2 ( u ( u)T ) . Interpreting Eq. (21) is more complicated than in incompressible flows, where a positive Q indicates that the vorticity magnitude exceeds the rate-of-strain magnitude. For the sake of completeness, the distribution of y+ on the inner surface also implies more turbulence at 90° and 270° (see Fig. 1). Fig. 11 exemplifies temperature time series at two locations in the mixing region. Although inner-surface temperatures were sampled at these very locations in the HWAT test section, the LES provides more information, in that it returns temperatures in the water layer adjacent to the wall, on the inner surface, and in the inner tube, all at the same ( , z ) and time. Fig. 11 confirms that the inner tube presents a low-pass filter behavior, since it damps high-frequency temperature fluctuations thanks to its thermal inertia: in Fig. 11a, temperatures in the water 10

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 10. Time-averaged Q-criterion isocontours colored by time-averaged temperature. Q = 0.5 s 2 . Start time t0 = 19.2 s . 17.5e 0.60 m z 0.73 m .

(90°, 0.67 m) (where mean ~488 K ) and (180°, 0.67 m) (where mean ~427 K ), they do not at (270°, 0.67 m) and (360°, 0.67 m) : LES temperatures average to 472 K at the former position and to 434 K at the latter, whereas, within the same time window, experimental temperatures average to 432 K at the former position and to 487 K at the latter; that is, their means seem to be roughly reversed. It is worth emphasizing that LES inlet velocities and temperatures in Table 2 average

3m

x2 + y2

22.5e

3 m.

experimental values obtained when inner-surface temperatures were sampled at 0.65 m. From this and similar results, such as Figs. mean (dimensionless T)_exp.pdf and mean(dimensionless T) _LES.pdf (Bergagio, 2019), it becomes clear that LES and experimental temperatures can be compared in terms of means at the axial, not azimuthal, coordinates of interest. We also remark that symmetry about the z-axis is more significant in LES results than in the 11

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 11. LES temperature time series on the inner surface, at WB, and at TB.

experimental ones, where frequency spectra of temperatures at (90°, 0.67 m) correlate more with those at (360°, 0.67 m) than with those at (270°, 0.67 m) . This correlation is clear if Kendall’s correlation coefficient is computed. measures the monotonic relationship between datasets. No assumption is made about their distributions. between Hilbert-Huang marginal spectra of temperatures at (90°, 0.67 m) and (270°, 0.67 m) is 0.805, while between Hilbert-Huang marginal spectra of temperatures at (90°, 0.67 m) and (360°, 0.67 m) reaches 0.857, thus suggesting higher correlation between spectra at (90°, 0.67 m) and (360°, 0.67 m) – see Fig. HMS1.pdf (Bergagio, 2019). Hilbert-Huang marginal spectra are described later. The larger experimental asymmetry can arise from many reasons, including (1) test-section deformation due to uneven heating and insulation of the annulus walls; (2) imperfect test-section design; and (3) uneven velocities and temperatures for each inlet pair, even more than in Table 2. Fig. 13 shows mean, highest, and lowest inner-surface temperatures along z, in the mixing region. Here, LES and experimental temperatures are normalized according to Eq. (22)

Tf , d =

Tf , d

TC

TH

TC

,

and lowest LES temperatures are consistent with their experimental equivalents. On the matter of differences between LES and experimental results, at = 3.333 (i.e., z = 0.65 m) , experimental temperatures exhibit a much broader spread around their mean. It is worth reminding that inner-surface temperatures were measured in the water domain using thermocouples whose tips are aligned with the inner surface, so larger deviations from the mean are to be expected in the experimental case. This difference in standard deviation is further explored in the comment on Fig. 14. Turbulent mixing is extremely chaotic and therefore sensitive to minute discrepancies in inlet conditions and flow geometry; accordingly, only a statistical comparison between computational and experimental temperatures is grounded. To this end, a mixing estimator is introduced, whose value is computed as the mean, outlier-free standard deviation of the inner-surface temperature time series at a given location (Bergagio et al., 2017). Computational and experimental mixing inhomogeneity is assessed 225°} and , which enclose in Fig. 14. Regions R1 = {( , z ): 135° the hot inlets, are colored orange and yellow, respectively. Both Figs. 14a and b suggest that mixing is slightly non-uniform at 0.63 m, 0.70 m, and 0.72 m. while its inhomogeneity exceeds 0.75 at (90°, 0.65 m) , (180°, 0.67 m) , (45°, 0.67 m) , and (225°, 0.67 m) in the LES case and at (45°, 0.65 m) , (360°, 0.65 m) , (90°, 0.65 m) , (225°, 0.67 m) , and (270°, 0.67 m) in the experimental case. Thus, mixing inhomogeneity often seems to stretch farther down in the experiment, especially in region R2 . Moreover, experimental peaks of mixing inhomogeneity are more comparable to one another than their LES analogues. Both these deviations could arise from the dissipative character of the SGS model and numerical schemes of choice. In addition, experimental mixing is more uneven in R2 than in R1 , whereas, in the LES case, mixing non-uniformity in R1 is slightly more consistent with that in R2 . Nevertheless, as previously observed, inasmuch as mass flow rates through the hot inlets are dissimilar, symmetry about the z-axis is compromised; to clarify, both LES and experimental mixing estimators 180°. For the sake of at differ from those at + 180° , where 0° <

(22)

where inner-surface temperatures are called Tf , d . Besides being lowpass filtered at 4 Hz (see Section 3), experimental temperatures are detrended using empirical mode decomposition (EMD). Contrastingly, LES temperatures are not filtered or detrended. TC and TH represent cold and hot inlet temperatures TC and TH , weighted over the respective mass flow rates and interpolated over time. Computational TC and TH are constant over time, so interpolation is not required. Dimensionless co(z z) ordinate is calculated as = 2(R H R ) , the midplane of the hot inlets oi io being located at zH = 0.8 m . In conformity with Fig. 10, below = 3.777 , or z = 0.63 m , and above = 1.777 , or z = 0.72 m , the spread around mean temperatures is so small that nonuniform mixing is deemed negligible. For ease of comparison, the same coordinates are examined for both LES and experiments. Overall, the mean, highest,

Fig. 12. Computational and experimental inner-surface temperature time series, at different angles and z = 0.67 m . 12

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 13. Axial distribution of the dimensionless inner-surface temperatures T f , d . Temperatures are sampled at the same azimuthal positions in both cases. The area between lower and upper envelopes is colored yellow. Here, T f , d denotes T f , d averaged over time and azimuthal positions.

thoroughness, a similar issue affects the distribution of y+ on the inner surface (see Fig. 1). This flaw in LES mixing estimators is not imputable to the mesh, in that it is perfectly symmetric. Experimental mixing nonuniformity appears even more distorted, possibly owing, as already mentioned, to asymmetries in the test-section geometry or variable boundary conditions. Overall, despite some disagreement between computational and experimental results, Figs. 13 and 14 confirm acceptable predictive capabilities of the LES model at hand. To this end, two metrics are considered: root mean squared error (RMSE) and mean absolute error (MAE). RMSE is defined as

1 N

RMSE =

sample sizes of axial and azimuthal distributions differ. MAE = 0.15 and RMSE = 0.20 for mixing inhomogeneity in the mixing region (N = 40 ). MAE and RMSE confirm that LES can predict time-averaged T f , d more accurately than mixing inhomogeneity, which is related to the standard deviation of the inner -surface temperatures. Concerning the elastic analysis of the inner tube, Fig. 15 shows timeaveraged stresses, close to the inner surface. Radial stresses of substantial magnitude are generated where mixing is most inhomogeneous. However, they are usually smaller in absolute value than axial and hoop stresses at the same coordinates. Overall, significant hoop stress gradients emerge from 0.63 m to 0.68 m, mainly in the axial direction; there, hoop stresses fluctuate from compressive to tensile (see Fig. 15b). Conversely, areas where the hot inlet jets impinge on the inner surface (z~0.80 m ), cold water is unmixed (z < 0.63 m ), or mixing is uniform, so that the temperature spread around the mean is small (0.70 m z 0.72 m , as shown in Fig. 13a), are associated with compressive hoop stresses of comparatively low magnitude. As shown in Fig. 15c, time-averaged axial stresses are one order of magnitude larger than time-averaged hoop stresses, most likely because the “stripe constraint” (Miyoshi et al., 2014; Kamaya and Miyoshi, 2017), prevails over the “circumferential constraint”. In essence, both constraints limit the deformation due to mixing non-uniformity, but the stripe constraint acts axially, thus increasing axial stresses, whereas the circumferential constraint acts azimuthally, thus enhancing hoop stresses instead. As inferable from Fig. 12a, hot vortices at 90° and 270° reach farther down than those at 180° and 360°, which is echoed in the distribution of the time-averaged axial stresses. As can be seen by comparing Fig. 14a to Fig. 15d, axial stress ranges are significant at locations where the mixing estimator is higher; in simple terms, where inner-surface

N

(MLES, i

Mexp, i ) 2 ,

i

(23)

where MLES , i and Mexp, i denote LES and experimental scalars, respectively. N is the sample size. MAE is defined as

MAE =

1 N

N

|MLES, i i

Mexp, i|.

(24)

RMSE is more sensitive to outliers than MAE. Here, RMSE and MAE are only computed for temperatures in the mixing region (0.63 m z 0.72 m ). MAE = 0.07 and RMSE = 0.10 for time-averaged dimensionless inner-surface temperatures T f , d (N = 40 ). Concerning T f , d – that is, T f , d averaged over time and azimuthal positions –, MAE = 0.03 and RMSE = 0.04 (N = 5). If T f , d is averaged over time and axial positions, MAE = 0.06 and RMSE = 0.07 (N = 8). Thus, the axial distribution of temperature T f , d might be better predicted than the azimuthal one; however, this remark should be taken cautiously, as the

Fig. 14. Mixing estimator as the average standard deviation of the dimensionless inner-surface temperatures T f , d after outlier rejection. All values are normalized. 13

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 15. Stresses in the vicinity of the inner surface, at 0.993 Rio . Values in Pa. Start time t0 = 19.2 s .

temperature time series reach higher variances. This is also true for Tjunctions, as noted by Kamaya and Nakamura (2011) and Costa Garrido et al. (2015). Axial and hoop stresses show similar ranges – see Fig. sttt_range.pdf (Bergagio, 2019). In order to identify which frequencies dominate inner-surface temperature spectra in the most damaged areas, more information on fatigue damage is necessary. Fig. 16 highlights that the CUF is highest between 0.65 m and 0.67 m, particularly from 45° to 135° and from 225° to 315°. The areas with the highest CUF are small, so fine meshing of the inner tube is necessary to capture them. The maximum CUF, which is related to the likelihood of initiating a fatigue crack, equals 3.2 × 10−5 over ~10 s .

Results in Fig. 17 are obtained by comparing CUFs on the inner surface at 0.6 m z 1 m . More specifically, Fig. 17a demonstrates that CUFs calculated using different start times t 0 and the same t typically converge with increasing start time. However, the selection of an appropriate start time cannot be exclusively based on convergence, since, with increasing start time, more load cycles from a given time history are lost. Concerning time-step independence, Fig. 17b stresses that CUFs calculated using the same t 0 and different time steps converge if the smaller t falls below 0.10 s. Further studies ought to evaluate mesh independence. In order to determine which frequencies dominate the spectra of the 14

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

termed gi , with i = 1 , …, G , while the final residual is called rG . The IMFs in subset {g1, …, gc0 1} form the decomposed time series, while the sum of rG and IMFs in subset {gc 0, …, gG } constitutes the underlying trend. The HSA is then applied to each IMF in the decomposed time series to determine its instantaneous frequencies and amplitudes. Both these quantities depend on time. Thus, the HHT circumvents the HeisenbergGabor principle, which limits the time-frequency resolution in Fourierbased representations. Our implementation of the HHT is detailed in Bergagio et al. (2017). The Hilbert-Huang marginal spectrum is derived from the HHT spectrum and quantifies how much each contributes to the total amplitude, over the whole time sequence. Index c0 is determined from the Hilbert-Huang marginal spectra of two successive IMFs, according to Yang et al. (2013). The Hilbert-Huang spectra of IMFs in subset S = {gc 0 4 , …, gc0 1} are also plotted. Fig. 18a depicts the Hilbert-Huang marginal spectrum of the LES temperature time series at the location where the CUF in Fig. 16 is highest. Since the total simulation time is much shorter than the experiment time (that is, 120 s), comparing the computational spectrum to the experimental one in Fig. 18b is challenging: the computational and experimental marginal spectra peak at 0.23 Hz and 0.04 Hz , respectively, so longer timescales dominate the experimental spectrum. Notwithstanding, in both cases, the largest peaks are associated with IMFs corresponding to long timescales, close to the underlying trends: the largest and second largest peaks in the computational spectrum are related to IMF g8 , while the largest and second largest peaks in the experimental spectrum are related to g9 and g8 , respectively. In both cases, c0 = 10 . Comparably to Fig. 18a, Fig. 19a portrays the Hilbert-Huang marginal spectrum of the LES temperature time series at the position where the CUF in Fig. 16 is second highest. In the computational spectrum, the dominant frequencies are 0.33 Hz and 0.51 Hz. They are mainly triggered by IMFs g7 and g6 , respectively, c0 being equal to 8. In the corresponding experimental spectrum (see Fig. 18b), the dominant frequencies are 0.04 Hz and 0.06 Hz. They are mainly caused by g9 , while c0 reaches 11. Even here, long timescales seem to contribute the most to the largest spectral peaks. Accordingly, our spectral analysis could benefit from running the LES longer. In Figs. 18 and 19, the LES spectra exhibit more pronounced peaks, whereas, overall, experimental temperatures show smoother transitions between adjacent peaks. This discrepancy could be due to the short simulation time. The histograms in Fig. 20a display which IMFs are most responsible for the two highest spectral peaks at locations most vulnerable to fatigue damage. Specifically, 20 inner-surface temperature time series are analyzed. Similarly to experiments (Bergagio et al., 2017), the IMFs behind these peaks come from the rightmost elements in subset S; in other words, they neighbor the underlying trend, which is especially true for the highest peaks. Fig. 20b indicates that the highest peaks at locations most prone to fatigue damage lie between 0.2 and 0.6 Hz, while the second highest peaks surface between 0.6 and 1 Hz. However, given the short simulation time, it might be safer to conclude that inner-surface temperatures at the locations most sensitive to fatigue damage exhibit rather low dominant frequencies, which are two-four times the inverse of the simulation time. Likewise, experimental inner-surface temperatures at the locations with strongest mixing inhomogeneity exhibit dominant frequencies between 0.03 and 0.06 Hz; i.e., three-seven times the inverse of the experiment time after low-pass filtering at 4 Hz (i.e., 109.96 s ; see Bergagio et al. (2017) and Bergagio and Anglart (2017)).

Fig. 16. CUF on the inner surface. Start time t0 = 19.2 s .

Fig. 17.

2

-norm of the difference between new and old CUFs.

LES inner-surface temperature time series in damage-susceptible areas, the measurement locations closest to areas where damage is at least 0.7 times the maximum CUF are examined. According to Fig. 16, these locations, sorted by CUF, are (45°, 0.67 m) , (90°, 0.65 m) , (90°, 0.67 m), (135°, 0.67 m) , and (315°, 0.67 m) . Interestingly, the highest and third highest mixing estimators in the LES case are found at (90°, 0.65 m) and (45°, 0.67 m) , respectively, which hints at a correlation between CUF and mixing estimator. As Fourier transforms and its variants, such as short-time Fourier transforms and wavelets, simulate non-stationary, nonlinear signals by introducing spurious harmonic components (see, e.g., Huang et al. (1998)), Hilbert-Huang marginal spectra are preferred. The HilbertHuang transform (HHT), a time-frequency-energy representation of a time series, combines EMD and Hilbert spectral analysis (HSA). More precisely, the EMD sifting procedure (see, e.g., Rilling et al. (2003)) decomposes the time series into a set of intrinsic, locally orthogonal oscillatory modes (IMFs), identified by their characteristic timescales in the series; hence, no a priori basis function is required. Here, each IMF is

5. Conclusions and suggestions for future work In this work, a CFD analysis of turbulent thermal mixing in an annulus has been successfully coupled to an FEA of the inner tube and to a fatigue damage assessment there. CFD simulation (i.e., an LES with the WALE SGS model), FEA, and fatigue damage assessment have been performed using free, open-source code. Temperatures and velocities at 15

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 18. Hilbert-Huang marginal spectra of computational and experimental inner-surface temperatures at (45°, 0.67 m) .

the inlets, together with pressure, are compatible with those from one of the first experiments on turbulent mixing under BWR conditions ( T = 216 K, p = 7.2 MPa ). Our analysis has focused on temperatures at the inner tube-water interface, here termed inner surface. LES innersurface temperatures agree with the experimental ones in terms of mixing inhomogeneity (acceptable agreement, with MAE = 0.15 and RMSE = 0.20 ), time-averaged values (good agreement, with MAE = 0.07 and RMSE = 0.10 ), and critical fluctuation components at the locations most susceptible to fatigue damage (good agreement). The LES dominant frequencies there are one order of magnitude larger than the experimental dominant frequencies at the locations where mixing inhomogeneity is strongest, probably because the simulation time (~10 s ) is shorter than the net experiment time (109.96 s ). This issue does not invalidate the present study: both LES and experimental dominant frequencies at the locations inspected are less than ten times the inverse of the respective time windows. In the current LES, strong mixing inhomogeneity reaches farther down at 90° and 270° than at 180° and 360°, possibly (1) because hotter vortices transport higher momentum fluid there from above, or (2) because colder, elongated vortices moving upwards induce more dissipation of turbulent kinetic energy there than the 2D structures at 180°

and 360° typically do, for the same axial levels. Time-averaged axial stresses are one order of magnitude larger than time-averaged hoop stresses: similarly to T-junctions, (1) the “stripe constraint” (Miyoshi et al., 2014; Kamaya and Miyoshi, 2017) overcomes the “circumferential constraint”, and (2) axial and hoop stresses show relevant ranges where mixing inhomogeneity is high. This analogy supports the validity of our CFD-FEA coupling method. Concerning fatigue damage, the highest CUFs appear in small areas, often where mixing inhomogeneity is highest, which strengthens the relevance of the inhomogeneity estimator. The maximum CUF, which indicates the likelihood of initiating a fatigue crack, is 3.2 × 10−5 over ~10 s . The numerical setup provides more data than the experimental setup. LES results can be used for further investigation, such as studies on buoyancy; evaluation of heat transfer coefficients at the inner surface; and testing of less dissipative SGS models and inverse heat transfer algorithms. y+, albeit acceptable in the mixing region, is deemed high at higher axial levels, so future simulations might benefit from increasing the mesh resolution. DES could be tested in order to reduce the computational time. However, mesh refinement would mainly decrease y+ in 16

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

Fig. 19. Hilbert-Huang marginal spectra of computational and experimental inner-surface temperatures at (90°, 0.65 m) .

areas far from those most prone to fatigue damage; consequently, dominant frequencies in critical areas are expected to be more sensitive to longer simulation times than to mesh refinement. Furthermore, because of Courant number constraints, refining the mesh would lead to smaller time steps and to a longer computational time. Thus, a longer

simulation time should be prioritized. Postprocessing can be enhanced by, e.g., addressing the mode mixing issue (Bergagio et al., 2017) and comparing the HHT with more traditional methods, such as wavelets. The fatigue assessment method can be improved by, e.g.,

Fig. 20. Peaks in the Hilbert-Huang marginal spectra of LES inner-surface temperatures at locations with the highest CUF. Start time t0 = 19.2 s . 17

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al.

introducing plasticity and environmental fatigue correction factors (Fen); building an LES-conformal mesh or testing different interpolation schemes from the CFD mesh to the FEA one; and performing a mesh independence study.

Jayaraju, S.T., Komen, E.M.J., Baglietto, E., 2010. Suitability of wall-functions in Large Eddy Simulation for thermal fatigue in a T-junction. Nucl. Eng. Des. 240 (10), 2544–2554. https://doi.org/10.1016/j.nucengdes.2010.05.026. ISSN 00295493. Jhung, M.J., 2013. Assessment of thermal fatigue in mixing tee by FSI analysis. Nucl. Eng. Technol. 45 (1), 99–106. https://doi.org/10.5516/NET.09.2012.026. Kamaya, M., Miyoshi, K., 2017. Thermal fatigue damage assessment at mixing tees (elastic-plastic deformation effect on stress and strain fluctuations). Nucl. Eng. Des. 318, 202–212. https://doi.org/10.1016/j.nucengdes.2017.04.022. Kamaya, M., Nakamura, A., 2011. Thermal stress analysis for fatigue damage evaluation at a mixing tee. Nucl. Eng. Des. 241 (8), 2674–2687. https://doi.org/10.1016/j. nucengdes.2011.05.029. Kim, S.-H., Huh, N.-S., Kim, M.-K., Cho, D.-G., Choi, Y.-H., Lee, J.-H., Choi, J.-B., 2013. Hydro-thermo-mechanical analysis on high cycle thermal fatigue induced by thermal striping in a T-junction. J. Mech. Sci. Technol. 27 (10), 3087–3095. https://doi.org/ 10.1007/s12206-013-0827-y. Kuczaj, A.K., Komen, E.M.J., 2010. An assessment of large-eddy simulation toward thermal fatigue prediction. Nucl. Technol. 170 (1), 2–15. https://doi.org/10.13182/ NT10-A9441. Lillberg, E., 2013. Predicting thermal mixing and fatigue inside control rod guide tubes. In: 21st International Conference on Nuclear Engineering, volume 3: Nuclear Safety and Security; Codes, Standards, Licensing and Regulatory Issues; Computational Fluid Dynamics and Coupled Codes. ASME, Chengdu, China. https://doi.org/10.1115/ ICONE21-16632. page V003T10A042. Logg, A., Mardal, K.-A., Wells, G., 2012. Automated solution of differential equations by the finite element method: The FEniCS book. Lecture Notes in Computational Science and Engineering, vol. 84 Springer978-3-642-23099-8https://doi.org/10.1007/978-3642-23099-8. Menter, F.R., 2015. Best practice: scale-resolving simulations in ANSYS CFD. ANSYS GmbH Technical report. Miyoshi, K., Nakamura, A., Utanohara, Y., Takenaka, N., 2014. An investigation of wall temperature characteristics to evaluate thermal fatigue at a T-junction pipe. Mech. Eng. J. 1 (5). https://doi.org/10.1299/mej.2014tep0050. TEP0050–TEP0050. Nicoud, F., Ducros, F., 1999. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence Combustion 62 (3), 183–200. https://doi. org/10.1023/A:1009995426001. Nicoud, F., Baya Toda, H., Cabrit, O., Bose, S., Lee, J., 2011. Using singular values to build a subgrid-scale model for large eddy simulations. Phys. Fluids 23 (8). https://doi.org/ 10.1063/1.3623274. 085106. Nie, B., Xu, F., 2019. Scales of natural convection on a convectively heated vertical wall. Phys. Fluids 31 (2), 024107. https://doi.org/10.1063/1.5083671. Odemark, Y., Green, T.M., Angele, K., Westin, J., Alavyoon, F., Lundström, S., 2009. Highcycle thermal fatigue in mixing tees: new large-eddy simulations validated against new data obtained by PIV in the vattenfall experiment. In: Proceedings of the 17th International Conference on Nuclear Engineering ICONE17. ASME, pp. 775–785. https://doi.org/10.1115/ICONE17-75962. OpenFOAM. User Guide version 5.0. The OpenFOAM Foundation; 2017. Pamiès, M., Weiss, P.-É., Garnier, E., Deck, S., Sagaut, P., 2009. Generation of synthetic turbulent inflow data for large eddy simulation of spatially evolving wall-bounded flows. Phys. Fluids 21 (4), 045103. https://doi.org/10.1063/1.3103881. Patankar, S., 1980. Numerical Heat Transfer and Fluid Flow. CRC Press9780891165224. Pegonen, R., Edh, N., Angele, K., Anglart, H., 2014. Investigation of thermal mixing in the control rod top tube using large eddy simulations. J. Power Technol. 94 (1), 67–78. Peltola, J., Pättikangas, T., Brockmann, T., Siikonen, T., Toppila, T., Brandt, T., 2011. Adaptation and validation of OpenFOAM CFD-solvers for nuclear safety related flow simulations. In: SAFIR2010, The Finnish Research Programme on Nuclear Power Plant Safety 2007-2010. Final Report. Espoo, Finland, pp. 283–294. Pope, S.B., 2000. Turbulent Flows. Cambridge Univ. Press, Cambridge [u.a.]. https://doi. org/10.1088/0957-0233/12/11/705. ISBN 9780511840531. Qian, S., Kanamaru, S., Kasahara, N., 2015. High-accuracy CFD prediction methods for fluid and structure temperature fluctuations at T-junction for thermal fatigue evaluation. Nucl. Eng. Des. 288, 98–109. https://doi.org/10.1016/j.nucengdes.2015.04. 006. Radu, V., Taylor, N., Paffumi, E., 2008. Development of new analytical solutions for elastic thermal stress components in a hollow cylinder under sinusoidal transient thermal loading. Int. J. Press. Vessels Pip. 85 (12), 885–893. https://doi.org/10. 1016/j.ijpvp.2008.04.010. Renze, P., Akermann, K., 2019. Simulation of conjugate heat transfer in thermal processes with open source CFD. ChemEngineering 3 (2), 59. https://doi.org/10.3390/ chemengineering3020059. Rilling, G., Flandrin, P., Gonçalvès, P., 2003. On empirical mode decomposition and its algorithms. In: Proceedings of IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing NSIP-03. Schijve, J., 2009. Fatigue of structures and materials, second ed. Springer, Dordrecht9781-4020-6808-9https://doi.org/10.1007/978-1-4020-6808-9. Karthick Selvam, P., Kulenovic, R., Laurien, E., Kickhofel, J., Prasser, H.-M., 2017. Thermal mixing of flows in horizontal T-junctions with low branch velocities. Nucl. Eng. Des. 322, 32–54. https://doi.org/10.1016/j.nucengdes.2017.06.041. Smith, B.L., Mahaffy, J.H., Angele, K., 2013. A CFD benchmarking exercise based on flow mixing in a T-junction. Nucl. Eng. Des. 264, 80–88. https://doi.org/10.1016/j. nucengdes.2013.02.030. JAHM Software Inc, 2015. Material Properties Database (MPDB) v. 8.21. North Reading, USA.https://www.jahm.com/, 2015. STAR-CCM+, 2018. User Manual. V13.02.011. CD-adapco. Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J. Numer. Anal. 21 (5), 995–1011. https://doi.org/10.1137/ 0721062.

Acknowledgments Mattia Bergagio and Roman Thiele acknowledge financial support through the THEMFA project by the Swedish Radiation Authority (SSM) and the Swedish Centre for Nuclear Technology (SKC). Wenyuan Fan is grateful for the support of China Scholarship Council (CSC no. 201600160035). The current LES exploited resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC Centre for High Performance Computing (PDC-HPC). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.nucengdes.2019.110361. References Addad, Y., Gaitonde, U., Laurence, D., Rolfo, S., 2008. Optimal unstructured meshing for large eddy simulations. In: Quality and Reliability of Large-Eddy Simulations. Ercoftac Series, vol. 12. Springer, pp. 93–103. Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N., 2015. The FEniCS project version 1.5. Arch. Numer. Software 3 (100), 9–23. https://doi.org/10.11588/ans.2015.100.20553. Angele, K., Odemark, Y., Cehlin, M., Hemström, B., Högström, C.M., Henriksson, M., Tinoco, H., Lindqvist, H., 2011. Flow mixing inside a control-rod guide tube – experimental tests and CFD simulations. Nucl. Eng. Des. 241 (12), 4803–4812. https:// doi.org/10.1016/j.nucengdes.2011.04.034. ASME, 2015. ASME Boiler & Pressure Vessel Code 2015 – Section III, Appendices.. ASME, 2015. ASME Boiler & Pressure Vessel Code 2015 - Section III, Division 1. Bergagio, M., 2019. Large eddy simulation of thermal mixing under boiling water reactor conditions. Mendeley Data v2. https://doi.org/10.17632/pdmymp8sxj.2. Bergagio, M., Anglart, H., 2017. Experimental investigation of mixing of non-isothermal water streams at BWR operating conditions. Nucl. Eng. Des. 317, 158–176. https:// doi.org/10.1016/j.nucengdes.2017.03.034. Bergagio, M., Thiele, R., Anglart, H., 2017. Analysis of temperature fluctuations caused by mixing of non-isothermal water streams at elevated pressure. Int. J. Heat Mass Transf. 104, 979–992. https://doi.org/10.1016/j.ijheatmasstransfer.2016.08.082. Chapuliot, S., Gourdin, C., Payen, T., Magnaud, J.P., Monavon, A., 2005. Hydro-thermalmechanical analysis of thermal fatigue in a mixing tee. Nucl. Eng. Des. 235 (5), 575–596. https://doi.org/10.1016/j.nucengdes.2004.09.011. Costa Garrido, O., El Shawish, S., Cizelj, L., 2015. Stress assessment in piping under synthetic thermal loads emulating turbulent fluid mixing. Nucl. Eng. Des. 283, 114–130. https://doi.org/10.1016/j.nucengdes.2014.10.016. Dahlberg, M., Nilsson, K.F., Taylor, N., Faidy, C., Wilke, U., Chapuliot, S., Kalkhof, D., Bretherton, I., Church, J.M., Solin, J., et al., 2007. Development of a European procedure for assessment of high cycle thermal fatigue in light water reactors. Final report of the NESC-Thermal fatigue project. European Commission Technical Report EUR 22763 EN. Dhamankar, N.S., Blaisdell, G.A., Lyrintzis, A.S., 2015. An Overview of turbulent inflow boundary conditions for large eddy simulations (Invited). In: 22nd AIAA Computational Fluid Dynamics Conference. American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/6.2015-3213. DTU Wind Energy, 2018. WindEnergyToolbox. Collection of HAWC2, HAWCSTable 2 tools etc from DTU Wind Energy.https://gitlab.windenergy.dtu.dk/toolbox/ WindEnergyToolbox. Online; (accessed September 3, 2018). Gallego-Marcos, I., 2013. Thermal Mixing CHT Simulations with OpenFOAM: URANS and LES. KTH Royal Institute of Technology (Master’s thesis). Gauder, P., Karthick Selvam, P., Kulenovic, R., Laurien, E., 2016. Large eddy simulation studies on the influence of turbulent inlet conditions on the flow behavior in a mixing tee. Nucl. Eng. Des. 298, 51–63. https://doi.org/10.1016/j.nucengdes.2015.12.001. Green, M., Ferrari, C.G., 2016. Desarrollo de un Procedimiento de Cálculo de Fatiga Multiaxial según ASME Sección III Incluyendo la Metodología Fen. Mecánica Computacional Vol XXXIV. pp. 1753–1772 (in Spanish). Haces Manzano, D., 2013. U-RANS and LES conjugate heat transfer analysis of the thermal fluctuations in a high difference temperature mixing zone using OpenFOAM. KTH Royal Institute of Technology (Master’s thesis). Huang, N.E., Shen, Z., Long, S.R., Wu, M.C., Shih, H.H., Zheng, Q., Yen, N.-C., Tung, C.C., Liu, H.H., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. A: Math., Phys. Eng. Sci. 454 (1971), 903–995. https://doi.org/10.1098/rspa.1998.0193. Hurrell, P.R., Philipson, R., Carter, K.F., Davenport, L., Wright, K., 2015. High cycle thermal fatigue analysis of pipe mixing tee with internal sleeve. In: Transactions, SMiRT-23, Manchester, 2015, p. 10.

18

Nuclear Engineering and Design xxx (xxxx) xxxx

M. Bergagio, et al. Tabor, G.R., Baba-Ahmadi, M.H., 2010. Inlet conditions for large eddy simulation: a review. Comput. Fluids 39 (4), 553–567. https://doi.org/10.1016/j.compfluid.2009. 10.007. Thiele, R., 2015. Mechanistic Modeling of Wall-Fluid Thermal Interactions for Innovative Nuclear Systems (Ph.D. thesis). KTH Royal Institute of Technology. Timperi, A., 2014. Conjugate heat transfer LES of thermal mixing in a T-junction. Nucl. Eng. Des. 273, 483–496. https://doi.org/10.1016/j.nucengdes.2014.02.031. Tinoco, H., 2013. CFD as a tool for the analysis of the mechanical integrity of light water nuclear reactors. In: Nuclear Reactor Thermal Hydraulics and Other Applications. InTech. https://doi.org/10.5772/52691. ISBN 978-953-51-0987-7. Tuominen, R., 2015. Coupling Serpent and OpenFOAM for neutronics-CFD Multi-physics Calculations. Master’s thesis. Aalto University. Wagner, W., Cooper, J.R., Dittmann, A., Kijima, J., Kretzschmar, H.-J., Kruse, A., Mareš, R., Oguchi, K., Sato, H., Stöcker, I., Šifner, O., Takaishi, Y., Tanishita, I., Trübenbach, J., Willkommen, Th., 2000. The IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam. J. Eng. Gas Turbines Power 122 (1), 150–184. https://doi.org/10.1115/1.483186. Westin, J., Veber, P., Andersson, L., ‘t Mannetje, C., Andersson, U., Eriksson, J., Henriksson, M.E., Alavyoon, F., Andersson, C., 2008. High-cycle thermal fatigue in

mixing tees: large-eddy simulations compared to a new validation experiment. In: 16th International Conference on Nuclear Engineering, volume 2: Fuel Cycle and High Level Waste Management; Computational Fluid Dynamics, Neutronics Methods and Coupled Codes; Student Paper Competition, 2008, pp. 515–525.https://doi.org/ 10.1115/ICONE16-48731. Wilson, J., Currie, C., Jones, M., Davenport, L., 2016. A case study evaluating the effects of high cycle thermal loading within a pressurised water reactor mixing tee using conjugate CFD/FE methods. In: ASME 2016 Pressure Vessels and Piping Conference. American Society of Mechanical Engineers, Vancouver, British Columbia, Canada. https://doi.org/10.1115/PVP2016-63233. Yang, Z., Ling, B.W.-K., Bingham, C., 2013. Trend extraction based on separations of consecutive empirical mode decomposition components in Hilbert marginal spectrum. Measurement 46 (8), 2481–2491. https://doi.org/10.1016/j.measurement. 2013.04.071. Zhang, Y., Lu, T., 2016. Study of the quantitative assessment method for high-cycle thermal fatigue of a T-pipe under turbulent fluid mixing based on the coupled CFDFEM method and the rainflow counting method. Nucl. Eng. Des. 309, 175–196. https://doi.org/10.1016/j.nucengdes.2016.09.021.

19