Thermal transients in a U-Bend

We study numerically the propagation of a hot thermal transient through a U-bend via an ensemble of wall-resolved large eddy simulations. Conjugate heat transfer between fluid and solid domains is accounted for. The flow is in a fully turbulent mixed convection regime, with a bulk Reynolds number of $10,000$, a Richardson number of $2.23$, and water as the working fluid (Prandtl number = $6$). These conditions lead to strong thermal stratification, with buoyancy-induced secondary flows, and the generation of a large and persistent recirculation region. The evolution of Dean vortices as the thermal transient passes is studied. It is found that baroclinic vorticity generation dominates over a large period of the transient, due to the thermal inertia of the wall. Gravitational buoyancy leads to a reversal of the counter-rotating vortex pair. The impact of this reversal on the swirl-switching and secondary-current losses is assessed. It is found that low frequency modes are suppressed in the reversed-vortex state.


Introduction
As a fluid flows around a bend, a local imbalance between the centripetal force and the opposing pressure gradient acts to transport low-inertia near-wall fluid towards the centre of curvature. This induces a counter-rotating vortex pair, known as Dean vortices, as fluid is subsequently transported back along the symmetry plane. The impact these Dean vortices have upon heat and mass transfer within the pipe is generally significant, and can be characterised by the dimensionless Dean number: Dn ≡ Re D(2R c ) −1 (where Re, D and R c are the bulk Reynolds number, the pipe inner-diameter and radius of curvature, respectively).
A number of studies have systematically quantified the behaviour of flow around bends with circular cross section. Improvements in numerical modelling of developing laminar flows in U-shaped pipe were suggested by Humphrey at el. [1], enabling better understanding of the secondary fluid motions. In a turbulent setting, secondary flow development was investigated by Azzozal et al. [2] both experimentally (using a laser Doppler anemometer) and numerically. These studies were continued by Baughn et al. [3] focusing on local heat transfer measurements for similar turbulent configurations. Although the problem of flows in curved pipes has been extensively studied over the past several decades, there is still a relative lack in understanding in the flow physics, and particularly the instabilities involved.
At moderate to high Dean numbers, instabilities in the flow manifest as a cyclical change in the strength and size of the left and right Dean vortices relative to one another [4]. This phenomenon is known as swirl-switching. A number of computational [4,5] and experimental [6,7,8] studies have investigated the physics of swirl-switching, and a review has been conducted by Vester et al. [9].
Transverse forces induced by swirl-switching can have significant impact upon the fatigue-life of plant components [10]. It is generally accepted the switching occurs as a gradual low-frequency shift between states, rather than an abrupt change from one bi-stable state to another [11].
A range of frequencies associated with swirl-switching have been reported.
Although reported values vary (often significantly) between different studies, it is generally accepted the swirl-switching has both low and high frequency modes.
For example, Hellström et al. [7] report two characteristic Strouhal numbers, St, of 0.16 and 0.33 (St ≡ f DU −1 , where U is the bulk velocity and f is the frequency). Recent work by Wang et al. [12] has highlighted that the recorded frequency of swirl-switching is highly sensitive to local flow properties, which may explain the discrepancies in recorded frequencies between different authors. Furthermore, their study highlights the role of large upstream structures on the low-frequency mode, while the higher-frequency mode originates from structures generated within the pipe-bend.
In the present study, we are concerned with the dynamics of flow in a Ubend as a hot thermal transient propagates throughout the domain. This case is relevant to the loop seal of a pressurised water nuclear fission reactor, amongst others. Both hot and cold thermal transients in a U-bend have been studied experimentally by Viollet [13] by applying a linear ramp to the inlet temperature over a short duration. Viollet observed that under the conditions of sufficiently low Reynolds number and Froude number, thermal stratification tends to occur. Stratification can significantly alter the flow, leading to the formation of a large buoyancy-induced recirculation region and steep temperature gradients.
Furthermore, cyclical changes in inlet temperature, alternating between hot and cold states, can lead to thermal fatigue.
The impact of buoyancy upon Dean vortices has been investigated in a number of studies. Lingrani et al. [14] studied experimentally the effect Dean vortices have upon surface heat transfer in a curved channel. They were primarily interested in forced convection flows, but faced challenges in removing buoyancy effects from their experiments. Mixed convection Nusselt numbers were therefore measured, and a correlation was given to convert mixed convection Nusselt numbers into forced convections ones.
Ciofalo et al. [15] investigated computationally the influence of both gravitational and centrifugal buoyancy for laminar flow within a coiled pipe. They performed a parameter study over a range of Richardson numbers with a linear increase in the pipe-wall temperature in the axial direction.
Kurnia et al. [16] performed calculations of laminar flow in straight and helical pipes of various cross section. They fixed the wall temperature and considered three different temperature differentials between wall and inlet. The dilatable working fluid (air) lead to conditions in which both gravitational and centrifugal buoyancy effects were significant. The study highlighted the formation of secondary vortices due to buoyancy in straight pipes, as well as the interaction between buoyancy-driven secondary flows and Dean vortices in the helical pipe cases.
In the aforementioned studies, [14,15,16], the temperature differential between the near-wall fluid and bulk-fluid was generated by heating (or cooling) the wall. To the best of our knowledge, the impact a thermal transient has upon Dean vortex dynamics and swirl-switching has not been previously studied. It is anticipated that a thermal transient, in conjunction with high wall thermal inertia, would lead to conditions in which large radial temperature (and density) gradients are present. In mixed convection flows, the gravitational and centrifugal buoyancy due to these density gradients has the potential to generate secondary flow, similar to those observed by Kurnia et al. [16]. The aim of this work is to test this hypothesis, and assess the impact a hot thermal transient has upon Dean vortex dynamics and swirl-switching within a U-bend configuration under turbulent mixed convection conditions. This paper is structured as follows. In Section 2 we outline our methodology.
In Section 3, key findings are presented. Finally, in Section 4 conclusions are drawn, and recommendations for future work are made.

Governing Equations
We employ an operator splitting strategy to decompose the problem into fluid and solid domains. The fluid flow is governed by the incompressible fil-tered Navier-Stokes equations (i.e. Large Eddy Simulations, (LES)) with the Boussinesq approximation to account for buoyancy (the impact of this approximation is to be discussed in Section 2.2), and a reduced filtered-energy equation to account for heat transfer: where an overbar denotes a filtered variable, and δ ij is the Kronecker delta.
The filtering operation is performed implicitly by the mesh. The field variables u, p and T denote the velocity, pressure and temperature, respectively, while t denotes the time. The constants are ρ 0 -the reference density, T 0 -the reference temperature, g -the gravitational acceleration, β -the thermal expansion coefficient, ν -the kinematic viscosity and α -the thermal diffusivity. Finally, τ ij and q j are the residual stress tensor and sub-grid heat flux, respectively.
The dynamic Smagorinsky eddy viscosity model provides turbulence closure [17,18,19]: where S ij is the resolved strain-rate tensor, S ≡ 2S ij S ij , ∆ is the local filter width, P r t is the turbulent Prandtl number, and c s is a model constant which is allowed to vary spatially and temporally; i.e. c s = c s (x, t). We dynamically set c s according to the Germano-Lilly procedure [18,19]. This dynamic sub-grid model was chosen since it is valid for relaminarised regions of the flow, as may be encountered in stably stratified turbulence.
In the solid domain, the transport of heat is accounted for via a reduced form of Equation 3 in which the convective and sub-grid terms are zero. The coupling between domains is achieved by enforcing consistency in both the temperature and heat flux at the interface: where subscripts (·) f and (·) s denote the fluid and solid domains, respectively, κ is the thermal conductivity, and n is the outward-pointing interface-normal.
The governing equations are discretised by the finite-volume method with ∆ taken as V 1/3 (where V is the local cell volume), and are solved with the CFD package Code Saturne (Version 5) [20].

Study Parameters
The fluid flow can be characterised by the Reynolds number (Re), gravitational Richardson number (Ri) and Prandtl number (P r): where U is the bulk velocity. In the present study, we set Re = 10, 000, Ri = 2.23, and Pr = 6. These parameters lead to mixed convection flow conditions that are based on the study of [13], and relevant to the nuclear industry (amongst others).
Note that by employing the Boussinesq approximation to account for buoyancy (see Section 2.1), the centrifugal buoyancy is neglected. The impact of this approximation can be quantified through consideration of the centrifugal Richardson number: The ratio of Ri c to Ri is then given by From the definition of Re, this can be rewritten as: which for the present study is less than 9 × 10 −4 if D > 0.2m. Similarly, the temperature-induced density difference is less than 0.5% if D > 0.2m, and hence the Bousinessq approximation is reasonable for sufficiently large D.
The final dimensionless groups dictating the rate of conductive heat flow between fluid and solid domains are the ratio of thermal diffusivities, α s /α f , and the ratio of thermal conductivities, κ s /κ f . In the present study, we employ a ratio of α s /α f = 144.8, and κ s /κ f = 123.5, which is representative of water flowing within a steel pipe. This is different to the experiments of [13], where altu-glass pipework was employed, and hence the thermal boundary condition was closer to adiabatic.

Geometry and Mesh
A schematic of the geometry is given in Figure 1. The vertical inlet and outlet sections are 10D in length, while the near-horizontal section is 6D in length with a 1% downward slope. This slope is a feature of the Viollet case [13], and is an approximation of the pipework downstream of the steam-generator outlet of the Superphénix reactor.
The radius of curvature, R c , for both bends is 1.5D, while the wall-thickness is 0.05D. All walls are smooth. This geometry is based on that of [13], but has a longer vertical inlet section. The reason for this change is to allow development of the inflow synthetic turbulence prior to the region of interest, as will be discussed in Section 2.4.
The data presented herein has been computed on a block-structured mesh comprising approximately 47M hexahedral cells (43M and 4M for the fluid and solid domains, respectively). The near-wall grid spacing was such that y + < 0.2 was maintained throughout the domain (with a corresponding T + < 1.2) and hence no additional near-wall modelling or damping terms were required (note the use of a dynamic sub-grid model, which precludes the need for near-wall damping).
A mesh sensitivity study has been conducted with a mesh comprising approximately 24M cells. We observed no appreciable difference in low-order statistics (mean velocity and temperature fields). The ensemble size (see Section 2.5 for details) for the coarse mesh (ten runs) was insufficient for converged higher-order statistics, hence the sensitivity of higher-order terms to the mesh cannot be ruled out.

Initial and Boundary Conditions
At the interface between fluid and solid domains, the no-slip condition provides a boundary condition for the velocity, while a zero-gradient Neumann condition is applied for the pressure. Unless otherwise stated, the interface temperature is set via the conjugate transfer of heat between the domains, as described in Section 2.1. At the external wall of the pipe, a zero gradient Neumann condition was applied for the temperature. Additional runs with an adiabatic wall thermal boundary condition have also been performed for comparison purposes.
Correlated inflow data, with prescribed first-and second-order statistics was generated via the Synthetic Eddy Method (SEM) [21,22]. The statistical input to the SEM was generated via a precursor Reynolds Averaged Navier-Stokes (RANS) simulation of fully-developed pipe-flow, also at Re = 10, 000, in which the Elliptic Blending Reynolds Stress Model (EBRSM) [23] provided turbulence closure. Note that we deliberately do not account for the change in the mean velocity profile and Reynolds stress tensor at the inlet due to upstream thermal stratification. Our aim is to assist reproducibility, and hence a simple station-ary condition for the velocity and Reynolds stresses was favoured. Moreover, tests showed that the synthetic turbulence reached a mature state within ∼ 2D (assessed via the development of wall skin-friction coefficient) -well before the first bend. This development length is in line with previous work at a similar Reynolds number [22]. At the outlet, zero-gradient Neumann conditions were applied for all variables, with the exception of the pressure which had a Dirichlet condition applied.

Statistical Data
Due to the transient nature of the flow, averaging is conducted via an ensemble. We generate a total of forty realisations of the flow, each by seeding the random number generator of the synthetic inflow method differently. This was done at the start of the isothermal initialisation computation. This ensemble size is sufficient for well-converged low-order statistics, and reasonable convergence of higher-order statistics (enough to be useful in the assessment of unsteady-RANS models, for example). The convergence of the ensemble is assessed via Table 1.

Flow and Heat Transfer Overview
In Figure 2, ensemble averaged streamlines showing the temporal evolution of the flow are presented. The formation of a large buoyancy-induced recirculation region is observed, qualitatively matching the observations of [13]. It can be seen that the thermal stratification persists for several dimensionless time-units, indicating limited heat transfer into the recirculation region.

Vorticity Budgets
We study the budgets of the vorticity transport equation in order to elucidate further the mechanism of secondary-flow reversal. The vorticity transport equation can be obtained by applying the curl operator to the momentum transport equation, and is written as where ijk is the Levi-Civita symbol. On the right hand side, II is the vortex stretching term, III is the baroclinic torque, and IV is the viscous diffusion of vorticity. The baroclinic term is active when the density gradient is misaligned with the pressure gradient. In our case, at the exit of the first bend, the pressure gradient balances the centrifugal and hydrostatic forces, and hence has a large vertical component. Meanwhile, the density gradient acts in approximately the wall-normal direction, due to the thermal inertia of the wall. At θ = π/2, we would thus expect the baroclinic term to be large as the transient passes, but before the wall has fully heated. In Figure 9, we plot the vorticity budgets in the wall-normal direction at θ = π/2. It can be seen that the baroclinic term, III, is the dominant source of vorticity production as the thermal transient passes.
This vorticity is redistributed by viscous diffusion (IV) which balances the baroclinic production term. Note that the residual of Equation (7) is not zero, as is apparent from studying Figure 9. This is due to the lack of convergence of the higher-order terms (the convective term in I, and the vortex stretching term, II) due to the sample size available for ensemble averaging. Despite this, it is apparent that the baroclinic term (which is well converged) is the dominant mode of vorticity production for times where the thermal transient is present.

Swirl-Switching
Swirl-switching can affect the fatigue life of pipes where the flow is subjected to a bend. We are interested in assessing how the thermal transient affects this phenomenon. Swirl-switching can be characterised through the swirl number, Sw, which is defined as the ratio of the flux of angular momentum in the axial direction and the flux of axial momentum in the axial direction, normalised by pipe radius: where U θ and U ax are the tangential and axial velocity components, respectively. In order to assess the frequency of the swirl-switching, we perform a continuous wavelet transform (CWT) to the Sw(t) signal, using a Morlet mother wavelet [24]. This allows us to transform the time signal to the frequency domain, and is applicable to a non-stationary signal. To reduce the noise in the spectral power estimation, we ensemble average the CWT coefficients over an ensemble size of three runs (note the smaller ensemble size than for the flow statistics, due to the availability of data). This process is analogous to Welch's method for spectral density estimation for a stationary signal. The resulting scalogram is presented in Figure 11. It can be seen from the figure that prior The low frequency mode that was present fort 10 has been suppressed. Since this low frequency mode is thought to originate from the large-scale structures upstream of the bend [12], one possible explanation for this observation is the turbulence-suppressing effects of stable stratification, which reduce the kinetic energy contained within the large structures. Figures 12 and 13 show contours of the normalised axial velocity 0.05D from the wall, for timest = 0 andt = 30, respectively. Also presented in the plots is a reconstruction using the first five POD modes, in order to highlight the most energetic structures. From Figure 12, long streak-like structures are clearly visible, typical of wall-bounded turbulence. Comparing this with Figure 13, it is apparent that the long structures are indeed significantly less pronounced due to the stable stratification, potentially explaining the lack of the low-frequency mode in the swirl-switching as the thermal transient passes.

Conclusion
In this paper, we studied the flow through a U-bend that is subject to a

Supplementary material
Full 3D volumetric datasets containing ensemble averaged data can be downloaded from [25].
D Pipe inner diameter.

Dn
Dean number.
f Frequency.
g Acceleration due to gravity.
R c Bend radius of curvature.

Ri
Gravitational Richardson number.
Ri c Centrifugal Richardson number.

Re
Reynolds number.
p Pressure.
P r Prandtl number.
P r t Turbulent Prandtl number.

Sw
Swirl number.
T Temperature.