- Email: [email protected]
Surrogate modelling of transient stratified flow in a u-shaped pipe
Nuclear Engineering and Design 368 (2020) 110787
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Surrogate modelling of transient stratified flow in a u-shaped pipe a,⁎
b,c
b
Małgorzata J. Zimoń , Alex Skillen , Wendi Liu , Robert Sawko a b c
T
a
IBM Research Europe, Daresbury Laboratory, UK Scientific Computing Department, STFC Daresbury Laboratory, UK Research and Innovation (UKRI), UK Department of Mechanical, Aerospace and Civil Engineering (MACE), School of Engineering, University of Manchester, UK
A R T I C LE I N FO
A B S T R A C T
Keywords: U-shaped pipe URANS Thermal transients Uncertainty quantification Surrogate modelling
In this work, we compute the mixed convection flow through a u-shaped bend, where transient thermal boundary conditions are applied at the inlet. The flow is governed by the Unsteady Reynolds Averaged NavierStokes equations. The Reynolds number is set to 10 000, while the Prandtl number is 6 (water). Polynomial chaos with parametric space discretisation is employed to assess the influence of Froude number (Fr), related to the magnitude of the thermal transient, uniformly distributed over a range of Fr = 0.402 to 0.938. Skewed and bimodal responses in the strength of the stratification were observed, which are related to the flow being at locally sub-critical or supercritical Fr. The application of polynomial expansion enables wider exploration of parameter space to gain better understanding of flow behaviour. For Fr < 0.67 , we observe a more persistent thermal stratification as the secondary flow motion confines the hotter core fluid to the centres of the two counter-rotating vortices away from the pipe wall.
1. Introduction Heat transfer by turbulent pipe flows occurs in variety of industrial applications, including heating and cooling networks, cryogenic processes, chemical reactors as well as components of internal combustion engines (Vashisth et al., 2008; Noorani et al., 2013; Giannakopoulos et al., 2017). As fluid flows through a bend, the local imbalance between centripetal forces and the opposing pressure gradient can lead to development of counter-rotating Dean vortices. Under certain conditions (sufficiently high Reynolds number, and sufficiently low radius of curvature), a low frequency oscillation in the size and strength of the Dean vortices relative to one another, can occur. This is known as swirlswitching. A number of computational (Rütten et al., 2005; Carlsson et al., 2015) and experimental studies (Sakakibara and Machida, 2012; Hellström et al., 2013; Kalpakli and Örlü, 2013) have identified this oscillatory secondary flow motion as a source of mechanical fatigue. Thermal transients, where the bulk temperature of the flow increases or decreases with time, are another potential source of fatigue in pipes. Thermal transients tend to induce thermal stratification in mixed convection conditions. This is a phenomenon that has been identified to occur frequently in nuclear power plants, e.g. in the pressuriser surge line, and other industrial processes (da Silva et al.,
⁎
2011; Qiao et al., 2014; Kumar et al., 2014). It is one of the significant contributors to material stresses in safety-related lines where there is potential for hot and cold fluids to come in contact (Schuler et al., 2004; Kumar et al., 2014). When a hot (lighter) fluid flows over a cold (heavier) one, large temperature gradients can occur, particularly at high Prandtl numbers. Both swirl-switching and thermal stratification, and hence both mechanical and thermal fatigue, have the potential to be significantly altered by secondary transport of heat and mass due to thermal transients. Skillen et al. (2020) showed that a thermal transient in mixed convection flows can induce secondary flow due to the thermal inertia of the wall. Since the wall takes some time to reach its equilibrium temperature, the near-wall fluid will be colder (and heavier) than that further from the wall (for a hot transient). This radial temperature gradient generates a baroclinic torque, due to gravitational buoyancy. In Skillen et al. (2020), the u-shaped piping at the Steam-Generator outlet of the Superphénix reactor was studied. The generation of a large and persistent recirculation region due to thermal stratification was observed, and it was shown that the wall’s thermal inertia causes a reversal of the counter-rotating Dean vortex pair. However, the study was limited to only one point in parameter space. The aim of this paper is to explore the secondary flow reversal over a range of mixed convection flow conditions, described in sections that follow, of relevance
Corresponding author. E-mail address: [email protected] (M.J. Zimoń).
https://doi.org/10.1016/j.nucengdes.2020.110787 Received 2 December 2019; Received in revised form 11 May 2020; Accepted 29 July 2020 Available online 19 August 2020 0029-5493/ © 2020 Elsevier B.V. All rights reserved.
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Reynolds-Stress Model (EBRSM) to close the equations (Manceau and Hanjalić, 2002). Turbulent thermal fluxes are accounted for via the Generalised Gradient Diffusion Hypothesis (GGDH) (Daly and Harlow, 1970). The buoyancy effects are included via a Boussinesq approximation in both the gravity source term and the turbulence production (see Viollet, 1987). The density difference at the smallest Fr is ~1.5% for the u-shaped pipe of diameter D = 0.2 , hence the Boussinesq approximation is reasonable. The conjugate heat transfer (CHT) simulation is performed with Code_Saturne version 5.0 (EDF, 2019) which couples the solid and fluid domains internally. In this coupling, the temperature and the heat flux are consistent at the interface:
Table 1 Summary of important dimensionless parameters. Parameter Reynolds Number Prandtl Number Reduced Froude Number
Equation
Re ≡
UD ν
Pr ≡ ν/ α Fr ≡
U gβ ΔTD
Value
10000 6 0.67 ± 40%
to the nuclear industry. In order to explore the parameter space efficiently, a minimal representation of the model can be constructed which allows us to perform multiple executions of different scenarios. Hereafter, we will refer to these approximations as surrogates. In the field of uncertainty quantification (UQ), surrogate-based techniques are used to measure the impact of varying input parameters on the numerical predictions. A key dimensionless parameter characterising stably stratified flow with internal gravity waves is the Froude number (see Table 1) (Viollet, 1987; Iyer and Theofanous, 1991; Rezende et al., 2011; Rezende et al., 2012). We treat Fr as random variable over a broad range. The practical utility of UQ relies heavily on the ability to perform a certain number of model runs, constraining significantly the choice of useful methods. Traditional Monte Carlo (MC) is often infeasible for computationally intensive systems due to its slow convergence rate, O (M−1/2) , where M is a number of simulations. Some techniques have been introduced to alleviate the cost of the brute-force approach (McKay et al., 1979; Dick et al., 2013; Giles, 2013). However, their usefulness is still limited to high-dimensional problems and non-smooth integrands, for which the Monte Carlo-based methods are the preferred solutions. An alternative to statistical methods for UQ are spectral expansion techniques. At their core, these are orthogonal decomposition methods, e.g. polynomial chaos expansion (PCE) or Karhunen-Loève, in which a random variable over a probability space, a quantity of interest (QoI), is expanded with respect to appropriate basis functions, e.g. polynomials or wavelets. The article focuses on the use of non-intrusive strategies which compute/estimate spectral coefficients via a set of deterministic model solutions. Such approaches enable us to treat the simulator as a black box. For practical applications, some of the constraints of the PCE, such as moderate number of random inputs, or presence of correlation between input processes, do not pose significant problems. In case of lowdimensional parameter spaces, PCE works very well. Some authors (Augustin et al., 2008) even argue that, if industrial relevance is considered, PCE is the best option, which in combination with the adaptive Gauss quadrature is a state-of-the-art method for solving stochastic Finite Elements. In the simulated flow, the presence of the stratification due to buoyancy force generates complex temperature distributions at different locations in the pipe. Therefore, high polynomial orders are needed to obtain convergence of stochastic moments. We show that, if the QoI is defined as a point solution of the considered model for a given input, then PCE with global basis functions can be inefficient. To improve the quality of surrogate models, we use a multi-element formulation of polynomial expression which adaptively decomposes the random space, enabling localised approximations. This article is organised as follows. The simulation details are provided in Section 2. Section 3 gives the theory behind surrogate modelling. Results and discussion are given in Section 4. Finally, we draw some conclusions and discuss future work.
Ts = Tf ,
κs
∂T ∂n
(1)
= −κf s
∂T ∂n
, f
(2)
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 interface conditions are updated explicitly at each time-step. The geometry and simulation parameters are based on the study performed by Viollet (Viollet, 1987) and are the same as in Skillen et al. (2020) (with the exception of Fr, which is considered a varying parameter or random variable). The pipework depicted in Fig. 1 has wall thickness of 0.05D . The length of both vertical sections is 10D , while the radius of curvature of both elbows is 1.5D . The near-horizontal section is 6D in length, with a slope of 1%. The properties of the system are expressed in terms of non-dimensional numbers summarised in Table 1,
2. Numerical modelling In this study the flow is governed by the Unsteady Reynolds Averaging Navier-Stokes (URANS) equations with the Elliptic Blending
Fig. 1. Dimensional details of the u-shaped pipe. 2
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
where U , g , ρ , ν , α, β , and ΔT are the bulk flow velocity, gravitational acceleration, fluid density, kinematic viscosity, thermal diffusivity, thermal expansion coefficient, and temperature difference between the initial and final states, respectively. We are considering water flowing within a steel pipe, where the ratio of thermal diffusivities and thermal conductivities are αsolid/ αfluid = 144.8 and κsolid/ κfluid = 123.5, respectively. The problem boundary conditions were set to be consistent with the reference Large Eddy Simulation (LES). 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. The interface temperature is set via the conjugate transfer of heat between the domains, as described previously. At the external wall of the pipe, a zero gradient Neumann condition was applied for the temperature. The inlet conditions were generated via a precursor Reynolds Averaged Navier-Stokes (RANS) simulation of fully-developed isothermal pipe-flow, also at Re = 10 000 , in which EBRSM model 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 as the thermal transient commences. This modelling decision was made in order to remain consistent with the reference LES (Skillen et al., 2020), which had constant fully-developed statistics at the inlet for the velocity and Reynolds stresses at all times. The flow was initialised by conducting an isothermal EBRSM computation at reference temperature T0 . The solution was time-marched to a steady state, to provide initial conditions to the main computation. At time t (̃ ≡Ut / D) = 0 , we initiated the thermal transient by linearly increasing the inlet temperature to a value T1. This linear ramp acts until t ̃ = 7.5. For all times t ̃ > 7.5, the inlet temperature remains fixed at T1. The magnitude of the temperature difference, ΔT ≡ (T1 − T0) , is dictated by the Froude number, Fr. At the outlet, zero-gradient Neumann conditions were applied for all variables, with the exception of the pressure which had a Dirichlet condition applied. For the URANS calculations, mesh sensitivity tests have been carried out with 1.2 million and 0.3 million cells for the fluid and solid domains, respectively. Fig. 2 shows the coarse (1.2 million fluid cells) and fine (4 million fluid cells) give broadly similar results. The results of the coarse mesh were also obtained with a time-step size twice that of the fine mesh (Δt = 0.02D / U , compared with 0.01D / U for the fine simulation), hence this test also serves as a sensitivity test of the time-step size. In both cases, a second order time-scheme was employed. In all figures the temperature is normalised as follows:
T − T0 , T1 − T0
(3)
A quantitative comparison is made in Fig. 3. It can be seen that the agreement is reasonable, and the fine solution can be considered mesh and time-step independent (to within a tolerance of a few percent). All subsequent calculations presented herein are performed with the fine mesh and time-step. Reference LES data at Fr = 0.67 from Skillen et al. (2020) is compared against our URANS results in Fig. 4. It can be seen that the temperature is predicted well by the EBRSM model. For the velocity profiles, the agreement between URANS and LES is slightly less satisfactory, although still broadly acceptable. The region of separated flow at the lower wall is well captured. Comparison with higher-order statistics has not been possible, due to the noise in higherorder statistics of the reference data. It is possible that the discrepancy may be remedied through the use of a more advanced thermal flux model. However, we found it was not possible to test this with the current implementation in Code_Saturne version 5.0. One of the primary goals of this study is to be able to accurately predict in-wall temperatures with a surrogate model. This is motivated by the fact thermal stratification can lead to large thermal stresses, and therefore an accurate prediction of the wall-temperature is of utmost importance. We therefore further test the predictive capabilities of our model against the reference LES for the fluid-solid interface temperature. Fig. 5 shows this comparison. It can be seen that the interface temperature is well predicted by the model prior to the unstably stratified region at the downstream vertical pipe (at z ′/ D ≳ 20 ).
3. Theoretical background In this section, we provide a brief description of generalised polynomial chaos expansion and its non-intrusive variant with adaptive random space decomposition. The classical PCE has already been applied to CFD; some examples can be found in Knio and Le Maître (2006); Najm (2009); Couaillier and Savin (2019). The modification of PCE for low-regular problems was first proposed by Wan and Karniadakis (2005) as multi-element generalised polynomial chaos, which yields approximate local statistics for uniformly distributed variables. Non-intrusive multi-element surrogate modelling was later described by Foo et al. (2008). This approach performs the probabilistic stochastic collocation in a decomposed random space. Here, we apply a similar logic as in the multi-element probabilistic collocation method (ME-PCM); in each subspace (element) we perform a pseudo-spectral (discrete) projection (Xiu, 2010) and then patch the local approximations together to model the global response.
Fig. 2. Temperature contours in the symmetry plane for coarse and fine URANS meshes at time Ut / D = 40 with Fr = 0.67 . 3
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 3. Profiles in the symmetry plane. Temperature (top) and streamwise velocity (bottom); Fine mesh (black) is compared against coarse mesh (red) at Ut / D = 30 and Fr = 0.67 . The solid domain is denoted with the blue background.
Fig. 4. Profiles in the symmetry plane. Temperature (top) and streamwise velocity (bottom); EBRSM model (black) is compared against an ensemble of 40 LES runs (red) at Ut / D = 30 and Fr = 0.67 . The solid domain is denoted with the blue background.
polynomial chaos expansion coefficients. The basis functions are orthogonal with respect to measure μ , i.e.
3.1. Generalised polynomial chaos We are considering the Froude number Fr (ξ ) to be a random variable, where ξ : Ω → has a certain probability measure μ . We want to quantify its effect on the temperature distribution along the pipe. The quantity of interest, Y (ω) , can be represented as a sum of weighted polynomials:
〈ψk , ψj 〉 =
∑ k=0
yk ψk (ξ (ω)),
for k , j = 0, …, Np
(5)
where 〈. ,.〉 is the inner product, δkj is the Kronecker delta. Therefore, the choice of the basis is imposed by the probability distribution of uncertain inputs. The normalisation constants γk are unique to the chosen basis and are known. The list of polynomials whose weights are associated with a particular random variable can be found in Xiu and Karniadakis (2002). In practice, the series is truncated at Np + 1 terms which represent
∞
Y (ω) =
∫Ω ψk (ξ ) ψj (ξ ) μ (ξ ) dξ = γk δkj
(4)
where ω is the random event and {yk }∞ k = 0 is referred to as a set of
Fig. 5. Flattened interface normalised temperature contours (cylindrical coordinate system). Dashed lines denote locations of the two bends. The inlet is at z ′/ D = 0 , where z ′ is the cylindrical streamwise direction. EBRSM model (top) is compared against an ensemble of 40 LES runs (bottom) at Ut / D = 30 and Fr = 0.67 . 4
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
threshold ∊t : abs(yn,̃ Np / yn,0 ̃ ) < ∊t .
the main contribution to the output variance. Construction of orthogonal polynomials for arbitrary probability density functions (PDFs) can be performed by exploiting the three-term recurrence relation (Sullivan, 2015). Due to the orthogonality of the basis, and assuming that first polynomial is a unit constant, the stochastic moments, i.e., the expected value, , and the variance, , can be expressed in terms of polynomial coefficients:
4. Surrogate modelling This section discusses the results of applying polynomial chaos approach to URANS model of thermal hydraulics of hot fluid flow in ushaped pipe. In order to accurately describe the characteristic of thermal stratification and reversal of secondary flow, we use a more sophisticated variant of PCE described in Section 3.3 to obtain the surrogate model. The approximation allows us, e.g., to predict temperature distributions for unseen values of Froude number from a given distribution, i.e., input values for which the simulations were not executed. To efficiently perform the study, we developed a wrapper (Zimoń et al., 2018; Zimoń et al., 2019) around the PCE library ChaosPy (Feinberg and Langtangen, 2015).
Np
[Y ] < Ψ0, Y > =
∑ yk
< Ψ0, Ψk> = y0 ,
(6)
k=0
Np
[Y ] = [∣Y − [Y ] ∣2 ] =
∑ yk2
< Ψk , Ψk> ,
(7)
k=0
where Np =
(p + d) ! p!d!
− 1 and Ψ(ξ ) are multivariate polynomials with
maximum degree p, ξ = {ξ1, …, ξd} and Ψ0 = 1.
4.1. Application of polynomial chaos expansion
3.2. Non-intrusive approach
In the present study, the Froude number has been defined as a random variable distributed uniformly on [0.402, 0.938]. The choice of the distribution suggests that we do not favour any particular parametric value within the domain of interest. In addition, we are looking at wide range of initial flow conditions which are more likely to result in complex response. To analyse different flow developments and potential risk related to induced thermal fluctuations at low frequencies, we measure the wall temperature distribution along the rings shown in Fig. 6. At first, the results are sampled at non-dimensionalised t ̃ ≡ Ut / D = 50 . In general, there is no standard procedure for choosing the polynomial order a priori, if the characteristics of the response function are unknown. When using global polynomials, one approach is to gradually increase p until the stochastic convergence is obtained, i.e. ∣ SDp − SDp − 1 ∣ < εSD , where ∣. ∣ is the Euclidean norm, SDp denotes the ∣ SD ∣
The PCE coefficients can be obtained non-intrusively (without reformulation of the model) either through regression, using least-squares solution to a linear system (Owen et al., 2017), or non-intrusive spectral projection (NISP) which projects the simulator output onto each polynomial term:
yk =
< Y , Ψk > 1 = γk < Ψ2k >
∫ Y (ξ )Ψk (ξ ) dμ (ξ ),
k = 0, …, Np.
(8)
To approximate these integrals, we can invoke quadrature rules. The realisations of Y are used to determine a surrogate model ∼ N Y = ∑k =p 0 yk̃ Ψk ≈ Y for which
1 y~k = γk
Q
∑ Y (ξ (j) )Ψk (ξ (j) ) w (j) . (9)
j=1
p−1
3.3. Multi-element polynomial chaos
standard deviation estimated with pth order polynomials and εSD is a predefined threshold. For low stochastic moments, a good estimate can be reached fairly quickly, particularly if the underlying function is smooth. However, if the accuracy of a surrogate in replicating the simulator is of importance, higher polynomial terms might need to be retained. In addition, if the response of the system is irregular, contains sharp transitions, steep gradients, or discontinuities, using global basis may never lead to accurate approximations. To be effective, refinement strategies could be employed, which minimise the number of grid points (in random space) based on local error estimations. To illustrate
Simulations of thermally stratified flows can lead to rapidly varying, irregular response functions. This makes global PCE inefficient when constructing accurate approximations, particularly when QoI is the state of the system at a given time and location. An adaptive framework utilising local basis function can offer greater accuracy in this case. The main idea behind the multi-element variant of PCE is simple: discretise the parametric space D into Ne non-overlapping subsets {Dn }nN=e 1 (i.e., D = ⋃nNe Dn and Dn1 ∩ Dn2 = ∅ if n1 ≠ n2 ), prescribe a set of collocation points (quadrature nodes) and estimate local surrogates. Since the degree of irregularity in each element is reduced, a relatively low polynomial order p can be maintained. In Foo et al. (2008) it was shown that such approach can achieve h − p convergence, where h denotes the size of random elements. The moments of Y (ξ ) on the entire domain are defined by the Bayes’ Rule and the law of total probability as a new random variable in each element is subject to a conditional PDF (Wan and Karniadakis, 2006; Foo et al., 2008). The decomposition of the current element is taking place if scaled value of a relative error of the PCE is higher than some threshold. In this study, for simplicity, we focus solely on the absolute value of the ratio of the last polynomial coefficient to the zeroth term; all of the ME-PCM results presented in Section 4 are obtained adaptively with a predefined
Fig. 6. Graphical representation of sampling regions. The rings are located at x / D ≈ 1.5, 3, 4.5, 6 and 7.5 for rings 1 to 5 respectively, and have radius 0.525D .
ξ (j)
and weights In this expression, the integration points (nodes) w (j) are computed using the theory of orthogonal polynomials. In the present work, we limit our analysis to Gaussian quadrature grids, which are not nested. In polynomial chaos, using Q nodes and weights, this quadrature formula can be made accurate for all polynomials of degree at most 2Q − 1, and no other formula on Q nodes has higher order of accuracy.
5
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 7. Error in standard deviation of temperature at t ̃ ≡ Ut / D = 50 in the 1st (a and b) and 2nd ring (c and d) for increasing polynomial order.
the estimation. For validation, assuming that the interpolation model ∼ generated with 250 runs is true, we plot the approximation error of Yp for the second ring in Fig. 10. The surrogates are constructed using p + 1 quadrature points (simulation outputs) and the results are based on 103 points sampled from the whole distribution of input Froude numbers. We can see that higher polynomial orders are required to represent the underlying function with global polynomials, i.e. εY20 < 0.02 . The resulting low-order models and corresponding PDFs at 3rd, 4th and 5th ring are plotted in Fig. 11. We can see that the functions are more noisy close to the second bend; using global polynomial approach is less successful in reproducing the simulation results. Stochastic moments are also slow-converging as shown in Fig. 12. This is because the hot plume has not yet affected the wall temperature in this location.
this, we initiate the study by constructing surrogates with global polynomials. To address the drawbacks of classical PCE, we then compare the approach with results given by using piece-wise basis. Knowing that the nature of a fluid response can be highly sensitive to the Froude number, for practical reasons, higher order basis functions should be considered. However, for the sake of the study, we initiate the NISP analysis with a small number of simulator executions (using polynomial orders as low as p = 1). Setting the threshold to εSD = 0.01 for convergence in both expected value E and standard deviation SD of temperatures in the first and second ring, the PDFs of temperature at each spatial point have been obtained. The plots in Fig. 7 show the relative difference of SD for different polynomial resolutions. As it can be seen, the standard deviation converges at p = 7 and p = 13 for the first and second ring, respectively, which means that p + 1 = 8 and p + 1 = 14 simulation evaluations were needed to obtain the stochastic response. The resulting PDF for the first ring in Fig. 8 highlights that the extreme values of wall temperature may appear to be more likely due to varying flow conditions. However, for the second ring the mean profile represents well the most probable temperature distribution. It should be noted that the PDFs in Fig. 8 are normalised based on maximum and minimum values at each radial location. The relative approximation errors obtained for a given position are low, i.e.
4.2. Multi-element surrogate modelling Although the system does not produce discontinuous responses, achieving converged statistics and accurate model representations requires large number of simulations even in one-dimension, i.e. d = 1. To address this problem, following the approach presented by Foo and Karniadakis (2010), we discretise the parametric space based on the behaviour of polynomial coefficients and construct local surrogates which are used to estimate the global response. The goal is to obtain an accurate approximation of the model for limited number of simulation results. For simplicity, we follow the procedure described in Foo et al. (2008), where we split the elements by half based on a given criterion. In our study, we investigate a local surrogate and refine the corresponding sub-domain if a ratio of the last coefficient to the zeroth one is higher than a pre-determined threshold. Such approach results in a more systematic study than analysis based purely on stochastic convergence. By keeping the polynomial order low for a given set of
∼ ∣ Ytrue − Yp ∣
εYp = ∣ Y ∣ < 0.05, where Ytrue is a response function approximated true with 250 runs for different samples of input Fr (marked as true model in ∼ Fig. 9) and Yp is a surrogate obtained with NISP using polynomials of order p. Note that the red circles in Fig. 9 indicate the model outputs (at Gaussian quadrature nodes) that are used to construct PCEs. We can see that due to more complex response for the second ring, the global polynomial expansion is not representing the simulation very accurately. Having performed the validation runs, we can continue increasing the polynomial order of the global surrogate for the second ring to assess how many expansion terms would have been required to improve 6
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 8. Resulting temperature distribution at t ̃ = 50 in the sections located near the first bend of the pipe. The PDFs are scaled using the minimum and maximum value of each profile individually. Therefore, the darker colour denotes local maxima of density functions.
threshold values, we generate multi-element polynomial representations. Together local surrogates reconstruct the response functions more accurately than global PCE for the same number of simulation runs. In Fig. 13 we plot the approximation error (computed assuming the interpolation of 250 URANS evaluations as true model) of solutions at t ̃ = 50 , obtained with PCE and multi-element approach, referred to as ME-PCM, for increasing number of simulation outputs used for expansion. The reconstructed functions describe the change of temperature at θ = 180° in selected rings due to varying Froude number. The multielement approach results in more efficient modelling as it adaptively divides the parametric space into elements that can be approximated with low-order basis as shown in Fig. 14. 4.3. Transient surrogate modelling At each time step, the system solutions can be evaluated only at the defined collocation points in the current set of elements. If at a given time one of the sub-domains needs to be divided, the conditions of the new element can be either directly transferred from the original element, if the random variables are equivalent, or the system has to be initialised for new set of nodal points. For transient turbulent flows, the complexity of the response function may vary substantially across the time domain. The more complex the response, the more improvement can be gained by splitting the parametric domain into elements with high regularity. To show this, we fix the threshold for decomposition to ∊t = 0.1 and the order of local bases (all of the local results are obtained with p = 2). We then allow the system to evolve from the initial flow
Fig. 10. Approximation error of PCE surrogates of the temperature in the 2nd ring constructed with increasing polynomial order p and p + 1 simulation runs.
condition and sample the temperature along the rings (see Fig. 6). Due to the transient nature of the flow, different responses are observed and, therefore, the same threshold may lead to higher number of sub-domains (more simulations required) at later instances. The comparison of approximation errors of surrogates for the first, second and third ring built with PCE and ME-PCM is shown in Fig. 15. For the first two rings,
Fig. 9. Surrogate representation of the response function for the first two rings at t ̃ = 50 at position θ = 180°. Note, red circles are the simulator outputs for which the surrogates are generated. The black line is the result of interpolated data from 250 simulations used to validate the reduced models. 7
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 11. Examples of surrogates of normalised temperature and the resulting PDFs for 3rd, 4th and 5th ring in the u-shaped pipe at t ̃ = 50 .
Fig. 12. Error in standard deviation of temperature in the 5th ring for increasing polynomial order. 8
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 13. Approximation errors of temperature response at t ̃ = 50 and position θ = 180° for given Froude number. Two methods are compared: surrogate modelling based on global expansion (PCE) and adaptive domain decomposition (ME-PCM). The improvement for the 5th ring for small number of evaluations is less substantial as the underlying function is noisy. In case of noisy results, it may be beneficial to filter the outputs before performing surrogate modelling.
the initial decomposition takes place at t ̃ = 35 and further refinement is performed at following times t ̃ = {40, 45, 50} if the last polynomial term is decaying too slowly. Although the division of parametric space performed at earlier times might not be optimal at later instances, the resulting error is still lower than in case of global approximation. For the 3rd ring we are showing the decomposition only at t ̃ = {45, 50} as there is not much change in the wall temperature at earlier times. As the study shows, the application of ME-PCM can improve the accuracy of surrogates, particularly when only small number of runs can be performed. For each local projection we use again the Gaussian quadrature rule. In the non-intrusive time-dependent analysis, the simulations are re-started for new nodal points when the domain element is split. 5. Results and discussion
t = 30 , for three In Fig. 16, we show the flow streamlines at time ~ different Froude numbers. It can be seen that the flow stratification that was observed in the LES of Skillen et al. (2020) is qualitatively well captured at Fr = 0.67 . At lower Froude numbers, the stratification is
Fig. 14. Reconstructed temperature response at θ = 180° in ring 3 for uncertain Froude number. The global surrogate is assembled from Ne = 5 approximations with p = 3. 9
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 15. Approximation errors of surrogates. The ME-PCM minimal models are constructed based on fixed threshold for convergence of polynomial terms and order of bases, p = 2 . The number of elements is equal to Ne = number of runs/(p + 1) . If Ne is increased when moving forward in time, the simulations in the new subdomains are evaluated from initial time for the new set of inputs. This is a limitation of using NISP with Gaussian quadratures. To reduce the cost, nested rules should be considered.
stronger, and hence the stratified recirculation region extends further upstream. On the other hand, at the highest Froude number under consideration, the stratified layer has almost completely disappeared by t = 30 . the time ~ Fig. 17 shows the fluid-solid interface temperature for different Froude numbers at time t ̃ = 50 . It can be seen that the temperature at the interface for the downstream vertical section is higher for lower values of Fr. The fluid in this region is unstably stratified, and hence the positive buoyancy production of turbulent kinetic energy will be enhanced at lower Froude number (since at low Fr, gravitational forces are large relative to the inertial forces). This enhanced buoyancy production leads to enhanced turbulent mixing, and hence higher average wall temperatures. In the near-horizontal section, the interface temperature towards the upper portion of the pipe (θ = 0 or 2π ) decreases with increasing Fr. This is due to the increased buoyancy effect at low Fr driving the hot transient upwards towards the upper wall (see Fig. 18). In the upstream vertical pipework, the interface temperature also decreases with increasing Fr. In this region, the thermal inertia of the wall plays a role, with the cold near-wall fluid being accelerated due to gravity in the streamwise direction, relative to the hot core fluid. This acceleration (more prominent at lower Fr since buoyancy forces are greater) draws the hotter fluid from the inlet further into the domain, leading to a thinning of both the thermal and momentum boundary layers and subsequent heating of the wall. We now consider the surrogate response for the second ring obtained with ME-PCM using Ne = 2 and Ne = 4 elements and p = 4 . In Fig. 19 we plot the Reynolds averaged temperature, as a function of Fr at locations θ = {0°, 180°} . We observe a complex response at the bottom of the ring (θ = 180°). As we saw in Fig. 16, the stratified recirculation region extends further upstream at lower Fr. If a given streamwise coordinate at θ = 180° is inside or outside the stratified recirculating flow will therefore depend on both the Froude number, and the time, t .̃ In Fig. 19b, we can see the bottom of the second ring is
t = 30 for three different Froude numbers, Fig. 16. Flow streamlines at time ~ coloured by normalised temperature (see legend in Fig. 17). In all cases, the seeding of the streamlines is at the same location.
10
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 17. Flattened interface temperature contours (cylindrical coordinate system). Dashed lines denote locations of the two bends. The inlet is at z ′/ D = 0 . Time t ̃ = 50 .
strength of the stratified recirculation region is decreasing with increasing Fr, hence the temperature increases. In Fig. 19a, we observe that the higher gravitational forces at lower Fr cause the hot fluid to rise towards the upper wall. This causes a decreasing wall temperature at θ = 0 with increasing Fr. In Fig. 20, we plot the secondary flow and temperature distribution at various Fr and streamwise locations in the near-horizontal pipework. At low Fr, the strength of the buoyancy-induced secondary flow due to the thermal inertia of the wall is relatively high, due to the greater relative strength of the gravitational forces at lower Fr. Strong counterrotating vortex pairs have the potential to confine the hot fluid to the vortex cores, away from the wall. This effect can be observed in Fig. 20 (e.g. at Fr = 0.458 and x / D = 4.5) where entrainment of the cold stratified fluid by the secondary flow leads to a rising column of cold fluid at the core of the pipe. The subsequent conversion of kinetic energy (of the cold secondary flow) to gravitational potential energy causes this
Fig. 18. Profiles of Reynolds averaged normalised temperature at different Fr. Red – Fr = 0.458; Black – Fr = 0.670 ; Green – Fr = 0.885. Solid domain denoted with a blue colour.
inside the recirculation region (at time t ̃ = 50 ) for Fr ≲ 0.55, and outside for higher Fr. This leads to the sharp jump in temperature at this location, since being outside the recirculation region leads to higher heating of the wall. At low Fr (prior to the jump), the temperature is increasing with increasing Fr. The location at the bottom of the second ring is inside the recirculation region at low Fr, but the extent and
Fig. 19. Surrogates of temperature in the second ring constructed with ME-PCM using polynomials of order p = 4 and N2 = 2 (a) and Ne = 4 (b) elements at different locations. Time t ̃ = 50 . 11
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
t = 30, Fr ∈ {0.458, 0.67, 0.885} . Outline of the solid Fig. 20. Secondary flow and temperature distribution at cross-sections through the near-horizontal section, at ~ domain denoted with a red line.
cold fluid to sink, creating the observed tight vortex pairs. As Fr increases further, the strength of the secondary flow decreases (since gravitational forces become less significant). The transport of heat towards the upper wall is therefore reduced, and hence the temperature at θ = 0° also reduces. At sufficiently high Fr, we would expect the secondary flow to revoke back to Dean vortices, which have the opposite direction to the buoyancy-induced vortices shown in Fig. 20. Although the simulations were not run at sufficiently high Fr to observe this effect fully, the onset of Dean vortex formation can be observed in Fig. 20, at Fr = 0.885 and x / D = 4.5, where there are two pairs of vortices (one small pair of Dean vortices, and another pair due to the buoyancy induced secondary flow). The effect this has is a thickening of the thermal boundary layer around the top of the pipe as the Dean vortex secondary flow transports fluid away from the upper part of the wall.
6. Concluding remarks We have used surrogate modelling in order to enable a systematic approach for analysing thermal transients and their impact upon vortex dynamics over a broad range of Fr. The results have shown that such study can be efficiently performed with a limited number of model evaluations. This is particularly beneficial when single simulation is computationally expensive. To improve the applicability of surrogates, more sophisticated adaptive criteria for space decomposition will be considered in future work. We observe secondary flows due to the thermal inertia of the wall; the near-wall fluid heats relatively slowly, which causes a radial temperature (and hence density) gradient. This was observed by Skillen et al. (2020) for a single point in parameter space. This is the first study, to the best of our knowledge, that extends the parameter space and assesses the physics of this secondary flow over a broad range of mixed 12
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
convection conditions that are of relevance to the nuclear industry. We have found that, the secondary flow has a tendency to confine the hotter core fluid to the centres of the two counter-rotating vortices, away from the pipe wall. For larger values of Fr, the opposite has been observed and the weak gravitational forces cause only small levels of flow recirculation and comparatively weak secondary flow. The polynomial chaos approach with parameter decomposition allows us to identify the most probable temperatures in time and space due to a large range of possible initial flow conditions. In some cases, the resulting distribution of temperature is bimodal in nature. The extreme values of the variable are more likely to occur than the average profile. This observation highlights misinterpretation of the system response in finite order statistics such as the mean and standard deviation.
Foo, J., Wan, X., Karniadakis, G.E., 2008. The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications. J. Comput. Phys. 227 (22), 9572–9595. Giannakopoulos, G.K., Frouzakis, C.E., Boulouchos, K., Fischer, P.F., Tomboulides, A.G., 2017. Direct numerical simulation of the flow in the intake pipe of an internal combustion engine. Int. J. Heat Fluid Flow 68, 257–268. Giles, M.B., 2013. Multilevel Monte Carlo methods. In: Monte Carlo and Quasi-Monte Carlo Methods 2012. Springer, pp. 83–103. Hellström, L.H., Zlatinov, M.B., Cao, G., Smits, A.J. Turbulent pipe flow downstream of a 90° bend. J. Fluid Mech. 735. Iyer, K.N., Theofanous, T.G., 1991. Decay of buoyancy-driven stratified layers with applications to pressurized thermal shock: reactor predictions. Nucl. Sci. Eng. 108 (2), 184–197. Kalpakli, A., Örlü, R., 2013. Turbulent pipe flow downstream a 90° pipe bend with and without superimposed swirl. Int. J. Heat Fluid Flow 41, 103–111. Knio, O., Le Maître, O., 2006. Uncertainty propagation in CFD using polynomial chaos decomposition. Fluid Dyn. Res. Kumar, R., Jadhav, P.A., Gupta, S.K., Gaikwad, A.J., 2014. Evaluation of thermal stratification induced stress in pipe and its impact on fatigue design. Proc. Eng. 86, 342–349. Manceau, R., Hanjalić, K., 2002. Elliptic blending model: a new near-wall Reynolds-stress turbulence closure. Phys. Fluids 14 (2), 744–754. McKay, M.D., Beckman, R.J., Conover, W.J., 1979. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245. Najm, H.N., 2009. Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics. Annu. Rev. Fluid Mech. 41, 35–52. Noorani, A., El Khoury, G.K., Schlatter, P., 2013. Evolution of turbulence characteristics from straight to curved pipes. Int. J. Heat Fluid Flow 41, 16–26. Owen, N.E., Challenor, P., Menon, P., Bennani, S., 2017. Comparison of surrogate-based uncertainty quantification methods for computationally expensive simulators. SIAM/ ASA J. Uncertainty Quantif. 5 (1), 403–435. Qiao, S., Gu, H., Wang, H., Luo, Y., Wang, D., Liu, P., Wang, Q., Mao, Q., 2014. Experimental investigation of thermal stratification in a pressurizer surge line. Ann. Nucl. Energy 73, 211–217. Rezende, H.C., Navarro, M.A., Mesquita, A.Z., Campagnole-dos Santos, A.A., Jordão, E., 2011. Experiments on one-phase thermally stratified flows in nuclear reactor pipe lines. Científica 15 (1), 17–24. Rezende, H.C., Santos, A.A.C., Navarro, M.A., Jordão, E., 2012. Verification and Validation of a thermal stratification experiment CFD simulation. Nucl. Eng. Des. 248, 72–81. Rütten, F., Schröder, W., Meinke, M., 2005. Large-eddy simulation of low frequency oscillations of the dean vortices in turbulent pipe bend flows. Phys. Fluids 17 (3) 035107. Sakakibara, J., Machida, N., 2012. Measurement of turbulent flow upstream and downstream of a circular pipe bend. Phys. Fluids 24 (4) 041702. Schuler, X., Herter, K.-H., 2004. Thermal fatigue due to stratification and thermal shock loading of piping. In: 30th MPA–Seminar in Conjunction with the 9th GermanJapanese Seminar, Stuttgart, Germany. pp. 6–1. Skillen, A., Zimoń, M.J., Sawko, R., Tunstall, R., Moulinec, C., Emerson, D.R., 2020. Thermal transients in a U-bend. Int. J. Heat Mass Transf. 148 119039. Sullivan, T.J., 2015. Introduction to Uncertainty Quantification, vol. 63 Springer. Vashisth, S., Kumar, V., Nigam, K.D.P., 2008. A review on the potential applications of curved geometries in process industry. Ind. Eng. Chem. Res. 47 (10), 3291–3337. Viollet, P.L., 1987. Observation and numerical modelling of density currents resulting from thermal transients in a non rectilinear pipe. J. Hydraul. Res. 25 (2), 235–261. Wan, X., Karniadakis, G.E., 2005. An adaptive multi-element generalized polynomial chaos method for stochastic differential equations. J. Comput. Phys. 209 (2), 617–642. Wan, X., Karniadakis, G.E., 2006. Multi-element generalized polynomial chaos for arbitrary probability measures. SIAM J. Scientific Comput. 28 (3), 901–928. Xiu, D., 2010. Numerical Methods for Stochastic Computations: A Apectral Method Approach. Princeton University Press. Xiu, D., Karniadakis, G.E., 2002. The Wiener-Askey polynomial chaos for stochastic differential equations. J. Scientific Comput. 24 (2), 619–644. Zimoń, M.J., Elisseev, V., Sawko, R., Antão, S., Jordan, K., 2018. Uncertainty quantification-as-a-service. In: Proceedings of the 28th Annual International Conference on Computer Science and Software Engineering. IBM Corp. pp. 331–337. Zimoń, M.J., Antão, S., Sawko, R., Skillen, A., Elisseev, V., 2019. Enabling UQ for complex modelling workflows. In: International Conference on Computational Science. Springer. pp. 269–281.
CRediT authorship contribution statement Małgorzata J. Zimoń: Methodology, Investigation, Writing - original draft. Alex Skillen: Formal analysis, Investigation, Writing - original draft. Wendi Liu: Writing - review & editing. Robert Sawko: Conceptualization, Software. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was supported by the STFC Hartree Centre’s Innovation Return on Research programme, funded by the Department for Business, Energy & Industrial Strategy. We thank Dr. Samuel Antão (IBM Research Europe) for leading the development of QUTE software that was used to generate discussed surrogate results. We acknowledge Dr. Charles Moulinec (STFC Daresbury Laboratory) for his help in modelling with Code_Saturne. References Augustin, F., Gilg, A., Paffrath, M., Rentrop, P., Wever, U., 2008. Polynomial chaos for the approximation of uncertainties: chances and limits. Eur. J. Appl. Math. 19 (2), 149–190. Carlsson, C., Alenius, E., Fuchs, L., 2015. Swirl switching in turbulent flow through 90 pipe bends. Phys. Fluids 27 (8) 085112. Couaillier, V., Savin, É., 2019. Generalized polynomial chaos for non-intrusive uncertainty quantification in computational fluid dynamics. In: Uncertainty Management for Robust Industrial Design in Aeronautics. Springer, pp. 123–141. Daly, B.J., Harlow, F.H., 1970. Transport equations in turbulence. Phys. Fluids 13 (11), 2634–2649. da Silva, L.L., Mansur, T.R., Cimini, C.A., 2011. Thermal fatigue damage evaluation of a pwr npp steam generator injection nozzle model subjected to thermal stratification phenomenon. Nucl. Eng. Des. 241 (3), 672–680. Dick, J., Kuo, F.Y., Sloan, I.H., 2013. High-dimensional integration: the quasi-Monte Carlo way. Acta Numerica 22, 133–288. EDF, Code_Saturne, 2019.https://www.code-saturne.org. Feinberg, J., Langtangen, H.P., 2015. Chaospy: an open source tool for designing methods of uncertainty quantification. J. Comput. Sci. 11, 46–57. Foo, J., Karniadakis, G.E., 2010. Multi-element probabilistic collocation method in high dimensions. J. Comput. Phys. 229 (5), 1536–1557.
13
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Surrogate modelling of transient stratified flow in a u-shaped pipe a,⁎
b,c
b
Małgorzata J. Zimoń , Alex Skillen , Wendi Liu , Robert Sawko a b c
T
a
IBM Research Europe, Daresbury Laboratory, UK Scientific Computing Department, STFC Daresbury Laboratory, UK Research and Innovation (UKRI), UK Department of Mechanical, Aerospace and Civil Engineering (MACE), School of Engineering, University of Manchester, UK
A R T I C LE I N FO
A B S T R A C T
Keywords: U-shaped pipe URANS Thermal transients Uncertainty quantification Surrogate modelling
In this work, we compute the mixed convection flow through a u-shaped bend, where transient thermal boundary conditions are applied at the inlet. The flow is governed by the Unsteady Reynolds Averaged NavierStokes equations. The Reynolds number is set to 10 000, while the Prandtl number is 6 (water). Polynomial chaos with parametric space discretisation is employed to assess the influence of Froude number (Fr), related to the magnitude of the thermal transient, uniformly distributed over a range of Fr = 0.402 to 0.938. Skewed and bimodal responses in the strength of the stratification were observed, which are related to the flow being at locally sub-critical or supercritical Fr. The application of polynomial expansion enables wider exploration of parameter space to gain better understanding of flow behaviour. For Fr < 0.67 , we observe a more persistent thermal stratification as the secondary flow motion confines the hotter core fluid to the centres of the two counter-rotating vortices away from the pipe wall.
1. Introduction Heat transfer by turbulent pipe flows occurs in variety of industrial applications, including heating and cooling networks, cryogenic processes, chemical reactors as well as components of internal combustion engines (Vashisth et al., 2008; Noorani et al., 2013; Giannakopoulos et al., 2017). As fluid flows through a bend, the local imbalance between centripetal forces and the opposing pressure gradient can lead to development of counter-rotating Dean vortices. Under certain conditions (sufficiently high Reynolds number, and sufficiently low radius of curvature), a low frequency oscillation in the size and strength of the Dean vortices relative to one another, can occur. This is known as swirlswitching. A number of computational (Rütten et al., 2005; Carlsson et al., 2015) and experimental studies (Sakakibara and Machida, 2012; Hellström et al., 2013; Kalpakli and Örlü, 2013) have identified this oscillatory secondary flow motion as a source of mechanical fatigue. Thermal transients, where the bulk temperature of the flow increases or decreases with time, are another potential source of fatigue in pipes. Thermal transients tend to induce thermal stratification in mixed convection conditions. This is a phenomenon that has been identified to occur frequently in nuclear power plants, e.g. in the pressuriser surge line, and other industrial processes (da Silva et al.,
⁎
2011; Qiao et al., 2014; Kumar et al., 2014). It is one of the significant contributors to material stresses in safety-related lines where there is potential for hot and cold fluids to come in contact (Schuler et al., 2004; Kumar et al., 2014). When a hot (lighter) fluid flows over a cold (heavier) one, large temperature gradients can occur, particularly at high Prandtl numbers. Both swirl-switching and thermal stratification, and hence both mechanical and thermal fatigue, have the potential to be significantly altered by secondary transport of heat and mass due to thermal transients. Skillen et al. (2020) showed that a thermal transient in mixed convection flows can induce secondary flow due to the thermal inertia of the wall. Since the wall takes some time to reach its equilibrium temperature, the near-wall fluid will be colder (and heavier) than that further from the wall (for a hot transient). This radial temperature gradient generates a baroclinic torque, due to gravitational buoyancy. In Skillen et al. (2020), the u-shaped piping at the Steam-Generator outlet of the Superphénix reactor was studied. The generation of a large and persistent recirculation region due to thermal stratification was observed, and it was shown that the wall’s thermal inertia causes a reversal of the counter-rotating Dean vortex pair. However, the study was limited to only one point in parameter space. The aim of this paper is to explore the secondary flow reversal over a range of mixed convection flow conditions, described in sections that follow, of relevance
Corresponding author. E-mail address: [email protected] (M.J. Zimoń).
https://doi.org/10.1016/j.nucengdes.2020.110787 Received 2 December 2019; Received in revised form 11 May 2020; Accepted 29 July 2020 Available online 19 August 2020 0029-5493/ © 2020 Elsevier B.V. All rights reserved.
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Reynolds-Stress Model (EBRSM) to close the equations (Manceau and Hanjalić, 2002). Turbulent thermal fluxes are accounted for via the Generalised Gradient Diffusion Hypothesis (GGDH) (Daly and Harlow, 1970). The buoyancy effects are included via a Boussinesq approximation in both the gravity source term and the turbulence production (see Viollet, 1987). The density difference at the smallest Fr is ~1.5% for the u-shaped pipe of diameter D = 0.2 , hence the Boussinesq approximation is reasonable. The conjugate heat transfer (CHT) simulation is performed with Code_Saturne version 5.0 (EDF, 2019) which couples the solid and fluid domains internally. In this coupling, the temperature and the heat flux are consistent at the interface:
Table 1 Summary of important dimensionless parameters. Parameter Reynolds Number Prandtl Number Reduced Froude Number
Equation
Re ≡
UD ν
Pr ≡ ν/ α Fr ≡
U gβ ΔTD
Value
10000 6 0.67 ± 40%
to the nuclear industry. In order to explore the parameter space efficiently, a minimal representation of the model can be constructed which allows us to perform multiple executions of different scenarios. Hereafter, we will refer to these approximations as surrogates. In the field of uncertainty quantification (UQ), surrogate-based techniques are used to measure the impact of varying input parameters on the numerical predictions. A key dimensionless parameter characterising stably stratified flow with internal gravity waves is the Froude number (see Table 1) (Viollet, 1987; Iyer and Theofanous, 1991; Rezende et al., 2011; Rezende et al., 2012). We treat Fr as random variable over a broad range. The practical utility of UQ relies heavily on the ability to perform a certain number of model runs, constraining significantly the choice of useful methods. Traditional Monte Carlo (MC) is often infeasible for computationally intensive systems due to its slow convergence rate, O (M−1/2) , where M is a number of simulations. Some techniques have been introduced to alleviate the cost of the brute-force approach (McKay et al., 1979; Dick et al., 2013; Giles, 2013). However, their usefulness is still limited to high-dimensional problems and non-smooth integrands, for which the Monte Carlo-based methods are the preferred solutions. An alternative to statistical methods for UQ are spectral expansion techniques. At their core, these are orthogonal decomposition methods, e.g. polynomial chaos expansion (PCE) or Karhunen-Loève, in which a random variable over a probability space, a quantity of interest (QoI), is expanded with respect to appropriate basis functions, e.g. polynomials or wavelets. The article focuses on the use of non-intrusive strategies which compute/estimate spectral coefficients via a set of deterministic model solutions. Such approaches enable us to treat the simulator as a black box. For practical applications, some of the constraints of the PCE, such as moderate number of random inputs, or presence of correlation between input processes, do not pose significant problems. In case of lowdimensional parameter spaces, PCE works very well. Some authors (Augustin et al., 2008) even argue that, if industrial relevance is considered, PCE is the best option, which in combination with the adaptive Gauss quadrature is a state-of-the-art method for solving stochastic Finite Elements. In the simulated flow, the presence of the stratification due to buoyancy force generates complex temperature distributions at different locations in the pipe. Therefore, high polynomial orders are needed to obtain convergence of stochastic moments. We show that, if the QoI is defined as a point solution of the considered model for a given input, then PCE with global basis functions can be inefficient. To improve the quality of surrogate models, we use a multi-element formulation of polynomial expression which adaptively decomposes the random space, enabling localised approximations. This article is organised as follows. The simulation details are provided in Section 2. Section 3 gives the theory behind surrogate modelling. Results and discussion are given in Section 4. Finally, we draw some conclusions and discuss future work.
Ts = Tf ,
κs
∂T ∂n
(1)
= −κf s
∂T ∂n
, f
(2)
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 interface conditions are updated explicitly at each time-step. The geometry and simulation parameters are based on the study performed by Viollet (Viollet, 1987) and are the same as in Skillen et al. (2020) (with the exception of Fr, which is considered a varying parameter or random variable). The pipework depicted in Fig. 1 has wall thickness of 0.05D . The length of both vertical sections is 10D , while the radius of curvature of both elbows is 1.5D . The near-horizontal section is 6D in length, with a slope of 1%. The properties of the system are expressed in terms of non-dimensional numbers summarised in Table 1,
2. Numerical modelling In this study the flow is governed by the Unsteady Reynolds Averaging Navier-Stokes (URANS) equations with the Elliptic Blending
Fig. 1. Dimensional details of the u-shaped pipe. 2
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
where U , g , ρ , ν , α, β , and ΔT are the bulk flow velocity, gravitational acceleration, fluid density, kinematic viscosity, thermal diffusivity, thermal expansion coefficient, and temperature difference between the initial and final states, respectively. We are considering water flowing within a steel pipe, where the ratio of thermal diffusivities and thermal conductivities are αsolid/ αfluid = 144.8 and κsolid/ κfluid = 123.5, respectively. The problem boundary conditions were set to be consistent with the reference Large Eddy Simulation (LES). 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. The interface temperature is set via the conjugate transfer of heat between the domains, as described previously. At the external wall of the pipe, a zero gradient Neumann condition was applied for the temperature. The inlet conditions were generated via a precursor Reynolds Averaged Navier-Stokes (RANS) simulation of fully-developed isothermal pipe-flow, also at Re = 10 000 , in which EBRSM model 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 as the thermal transient commences. This modelling decision was made in order to remain consistent with the reference LES (Skillen et al., 2020), which had constant fully-developed statistics at the inlet for the velocity and Reynolds stresses at all times. The flow was initialised by conducting an isothermal EBRSM computation at reference temperature T0 . The solution was time-marched to a steady state, to provide initial conditions to the main computation. At time t (̃ ≡Ut / D) = 0 , we initiated the thermal transient by linearly increasing the inlet temperature to a value T1. This linear ramp acts until t ̃ = 7.5. For all times t ̃ > 7.5, the inlet temperature remains fixed at T1. The magnitude of the temperature difference, ΔT ≡ (T1 − T0) , is dictated by the Froude number, Fr. At the outlet, zero-gradient Neumann conditions were applied for all variables, with the exception of the pressure which had a Dirichlet condition applied. For the URANS calculations, mesh sensitivity tests have been carried out with 1.2 million and 0.3 million cells for the fluid and solid domains, respectively. Fig. 2 shows the coarse (1.2 million fluid cells) and fine (4 million fluid cells) give broadly similar results. The results of the coarse mesh were also obtained with a time-step size twice that of the fine mesh (Δt = 0.02D / U , compared with 0.01D / U for the fine simulation), hence this test also serves as a sensitivity test of the time-step size. In both cases, a second order time-scheme was employed. In all figures the temperature is normalised as follows:
T − T0 , T1 − T0
(3)
A quantitative comparison is made in Fig. 3. It can be seen that the agreement is reasonable, and the fine solution can be considered mesh and time-step independent (to within a tolerance of a few percent). All subsequent calculations presented herein are performed with the fine mesh and time-step. Reference LES data at Fr = 0.67 from Skillen et al. (2020) is compared against our URANS results in Fig. 4. It can be seen that the temperature is predicted well by the EBRSM model. For the velocity profiles, the agreement between URANS and LES is slightly less satisfactory, although still broadly acceptable. The region of separated flow at the lower wall is well captured. Comparison with higher-order statistics has not been possible, due to the noise in higherorder statistics of the reference data. It is possible that the discrepancy may be remedied through the use of a more advanced thermal flux model. However, we found it was not possible to test this with the current implementation in Code_Saturne version 5.0. One of the primary goals of this study is to be able to accurately predict in-wall temperatures with a surrogate model. This is motivated by the fact thermal stratification can lead to large thermal stresses, and therefore an accurate prediction of the wall-temperature is of utmost importance. We therefore further test the predictive capabilities of our model against the reference LES for the fluid-solid interface temperature. Fig. 5 shows this comparison. It can be seen that the interface temperature is well predicted by the model prior to the unstably stratified region at the downstream vertical pipe (at z ′/ D ≳ 20 ).
3. Theoretical background In this section, we provide a brief description of generalised polynomial chaos expansion and its non-intrusive variant with adaptive random space decomposition. The classical PCE has already been applied to CFD; some examples can be found in Knio and Le Maître (2006); Najm (2009); Couaillier and Savin (2019). The modification of PCE for low-regular problems was first proposed by Wan and Karniadakis (2005) as multi-element generalised polynomial chaos, which yields approximate local statistics for uniformly distributed variables. Non-intrusive multi-element surrogate modelling was later described by Foo et al. (2008). This approach performs the probabilistic stochastic collocation in a decomposed random space. Here, we apply a similar logic as in the multi-element probabilistic collocation method (ME-PCM); in each subspace (element) we perform a pseudo-spectral (discrete) projection (Xiu, 2010) and then patch the local approximations together to model the global response.
Fig. 2. Temperature contours in the symmetry plane for coarse and fine URANS meshes at time Ut / D = 40 with Fr = 0.67 . 3
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 3. Profiles in the symmetry plane. Temperature (top) and streamwise velocity (bottom); Fine mesh (black) is compared against coarse mesh (red) at Ut / D = 30 and Fr = 0.67 . The solid domain is denoted with the blue background.
Fig. 4. Profiles in the symmetry plane. Temperature (top) and streamwise velocity (bottom); EBRSM model (black) is compared against an ensemble of 40 LES runs (red) at Ut / D = 30 and Fr = 0.67 . The solid domain is denoted with the blue background.
polynomial chaos expansion coefficients. The basis functions are orthogonal with respect to measure μ , i.e.
3.1. Generalised polynomial chaos We are considering the Froude number Fr (ξ ) to be a random variable, where ξ : Ω → has a certain probability measure μ . We want to quantify its effect on the temperature distribution along the pipe. The quantity of interest, Y (ω) , can be represented as a sum of weighted polynomials:
〈ψk , ψj 〉 =
∑ k=0
yk ψk (ξ (ω)),
for k , j = 0, …, Np
(5)
where 〈. ,.〉 is the inner product, δkj is the Kronecker delta. Therefore, the choice of the basis is imposed by the probability distribution of uncertain inputs. The normalisation constants γk are unique to the chosen basis and are known. The list of polynomials whose weights are associated with a particular random variable can be found in Xiu and Karniadakis (2002). In practice, the series is truncated at Np + 1 terms which represent
∞
Y (ω) =
∫Ω ψk (ξ ) ψj (ξ ) μ (ξ ) dξ = γk δkj
(4)
where ω is the random event and {yk }∞ k = 0 is referred to as a set of
Fig. 5. Flattened interface normalised temperature contours (cylindrical coordinate system). Dashed lines denote locations of the two bends. The inlet is at z ′/ D = 0 , where z ′ is the cylindrical streamwise direction. EBRSM model (top) is compared against an ensemble of 40 LES runs (bottom) at Ut / D = 30 and Fr = 0.67 . 4
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
threshold ∊t : abs(yn,̃ Np / yn,0 ̃ ) < ∊t .
the main contribution to the output variance. Construction of orthogonal polynomials for arbitrary probability density functions (PDFs) can be performed by exploiting the three-term recurrence relation (Sullivan, 2015). Due to the orthogonality of the basis, and assuming that first polynomial is a unit constant, the stochastic moments, i.e., the expected value, , and the variance, , can be expressed in terms of polynomial coefficients:
4. Surrogate modelling This section discusses the results of applying polynomial chaos approach to URANS model of thermal hydraulics of hot fluid flow in ushaped pipe. In order to accurately describe the characteristic of thermal stratification and reversal of secondary flow, we use a more sophisticated variant of PCE described in Section 3.3 to obtain the surrogate model. The approximation allows us, e.g., to predict temperature distributions for unseen values of Froude number from a given distribution, i.e., input values for which the simulations were not executed. To efficiently perform the study, we developed a wrapper (Zimoń et al., 2018; Zimoń et al., 2019) around the PCE library ChaosPy (Feinberg and Langtangen, 2015).
Np
[Y ] < Ψ0, Y > =
∑ yk
< Ψ0, Ψk> = y0 ,
(6)
k=0
Np
[Y ] = [∣Y − [Y ] ∣2 ] =
∑ yk2
< Ψk , Ψk> ,
(7)
k=0
where Np =
(p + d) ! p!d!
− 1 and Ψ(ξ ) are multivariate polynomials with
maximum degree p, ξ = {ξ1, …, ξd} and Ψ0 = 1.
4.1. Application of polynomial chaos expansion
3.2. Non-intrusive approach
In the present study, the Froude number has been defined as a random variable distributed uniformly on [0.402, 0.938]. The choice of the distribution suggests that we do not favour any particular parametric value within the domain of interest. In addition, we are looking at wide range of initial flow conditions which are more likely to result in complex response. To analyse different flow developments and potential risk related to induced thermal fluctuations at low frequencies, we measure the wall temperature distribution along the rings shown in Fig. 6. At first, the results are sampled at non-dimensionalised t ̃ ≡ Ut / D = 50 . In general, there is no standard procedure for choosing the polynomial order a priori, if the characteristics of the response function are unknown. When using global polynomials, one approach is to gradually increase p until the stochastic convergence is obtained, i.e. ∣ SDp − SDp − 1 ∣ < εSD , where ∣. ∣ is the Euclidean norm, SDp denotes the ∣ SD ∣
The PCE coefficients can be obtained non-intrusively (without reformulation of the model) either through regression, using least-squares solution to a linear system (Owen et al., 2017), or non-intrusive spectral projection (NISP) which projects the simulator output onto each polynomial term:
yk =
< Y , Ψk > 1 = γk < Ψ2k >
∫ Y (ξ )Ψk (ξ ) dμ (ξ ),
k = 0, …, Np.
(8)
To approximate these integrals, we can invoke quadrature rules. The realisations of Y are used to determine a surrogate model ∼ N Y = ∑k =p 0 yk̃ Ψk ≈ Y for which
1 y~k = γk
Q
∑ Y (ξ (j) )Ψk (ξ (j) ) w (j) . (9)
j=1
p−1
3.3. Multi-element polynomial chaos
standard deviation estimated with pth order polynomials and εSD is a predefined threshold. For low stochastic moments, a good estimate can be reached fairly quickly, particularly if the underlying function is smooth. However, if the accuracy of a surrogate in replicating the simulator is of importance, higher polynomial terms might need to be retained. In addition, if the response of the system is irregular, contains sharp transitions, steep gradients, or discontinuities, using global basis may never lead to accurate approximations. To be effective, refinement strategies could be employed, which minimise the number of grid points (in random space) based on local error estimations. To illustrate
Simulations of thermally stratified flows can lead to rapidly varying, irregular response functions. This makes global PCE inefficient when constructing accurate approximations, particularly when QoI is the state of the system at a given time and location. An adaptive framework utilising local basis function can offer greater accuracy in this case. The main idea behind the multi-element variant of PCE is simple: discretise the parametric space D into Ne non-overlapping subsets {Dn }nN=e 1 (i.e., D = ⋃nNe Dn and Dn1 ∩ Dn2 = ∅ if n1 ≠ n2 ), prescribe a set of collocation points (quadrature nodes) and estimate local surrogates. Since the degree of irregularity in each element is reduced, a relatively low polynomial order p can be maintained. In Foo et al. (2008) it was shown that such approach can achieve h − p convergence, where h denotes the size of random elements. The moments of Y (ξ ) on the entire domain are defined by the Bayes’ Rule and the law of total probability as a new random variable in each element is subject to a conditional PDF (Wan and Karniadakis, 2006; Foo et al., 2008). The decomposition of the current element is taking place if scaled value of a relative error of the PCE is higher than some threshold. In this study, for simplicity, we focus solely on the absolute value of the ratio of the last polynomial coefficient to the zeroth term; all of the ME-PCM results presented in Section 4 are obtained adaptively with a predefined
Fig. 6. Graphical representation of sampling regions. The rings are located at x / D ≈ 1.5, 3, 4.5, 6 and 7.5 for rings 1 to 5 respectively, and have radius 0.525D .
ξ (j)
and weights In this expression, the integration points (nodes) w (j) are computed using the theory of orthogonal polynomials. In the present work, we limit our analysis to Gaussian quadrature grids, which are not nested. In polynomial chaos, using Q nodes and weights, this quadrature formula can be made accurate for all polynomials of degree at most 2Q − 1, and no other formula on Q nodes has higher order of accuracy.
5
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 7. Error in standard deviation of temperature at t ̃ ≡ Ut / D = 50 in the 1st (a and b) and 2nd ring (c and d) for increasing polynomial order.
the estimation. For validation, assuming that the interpolation model ∼ generated with 250 runs is true, we plot the approximation error of Yp for the second ring in Fig. 10. The surrogates are constructed using p + 1 quadrature points (simulation outputs) and the results are based on 103 points sampled from the whole distribution of input Froude numbers. We can see that higher polynomial orders are required to represent the underlying function with global polynomials, i.e. εY20 < 0.02 . The resulting low-order models and corresponding PDFs at 3rd, 4th and 5th ring are plotted in Fig. 11. We can see that the functions are more noisy close to the second bend; using global polynomial approach is less successful in reproducing the simulation results. Stochastic moments are also slow-converging as shown in Fig. 12. This is because the hot plume has not yet affected the wall temperature in this location.
this, we initiate the study by constructing surrogates with global polynomials. To address the drawbacks of classical PCE, we then compare the approach with results given by using piece-wise basis. Knowing that the nature of a fluid response can be highly sensitive to the Froude number, for practical reasons, higher order basis functions should be considered. However, for the sake of the study, we initiate the NISP analysis with a small number of simulator executions (using polynomial orders as low as p = 1). Setting the threshold to εSD = 0.01 for convergence in both expected value E and standard deviation SD of temperatures in the first and second ring, the PDFs of temperature at each spatial point have been obtained. The plots in Fig. 7 show the relative difference of SD for different polynomial resolutions. As it can be seen, the standard deviation converges at p = 7 and p = 13 for the first and second ring, respectively, which means that p + 1 = 8 and p + 1 = 14 simulation evaluations were needed to obtain the stochastic response. The resulting PDF for the first ring in Fig. 8 highlights that the extreme values of wall temperature may appear to be more likely due to varying flow conditions. However, for the second ring the mean profile represents well the most probable temperature distribution. It should be noted that the PDFs in Fig. 8 are normalised based on maximum and minimum values at each radial location. The relative approximation errors obtained for a given position are low, i.e.
4.2. Multi-element surrogate modelling Although the system does not produce discontinuous responses, achieving converged statistics and accurate model representations requires large number of simulations even in one-dimension, i.e. d = 1. To address this problem, following the approach presented by Foo and Karniadakis (2010), we discretise the parametric space based on the behaviour of polynomial coefficients and construct local surrogates which are used to estimate the global response. The goal is to obtain an accurate approximation of the model for limited number of simulation results. For simplicity, we follow the procedure described in Foo et al. (2008), where we split the elements by half based on a given criterion. In our study, we investigate a local surrogate and refine the corresponding sub-domain if a ratio of the last coefficient to the zeroth one is higher than a pre-determined threshold. Such approach results in a more systematic study than analysis based purely on stochastic convergence. By keeping the polynomial order low for a given set of
∼ ∣ Ytrue − Yp ∣
εYp = ∣ Y ∣ < 0.05, where Ytrue is a response function approximated true with 250 runs for different samples of input Fr (marked as true model in ∼ Fig. 9) and Yp is a surrogate obtained with NISP using polynomials of order p. Note that the red circles in Fig. 9 indicate the model outputs (at Gaussian quadrature nodes) that are used to construct PCEs. We can see that due to more complex response for the second ring, the global polynomial expansion is not representing the simulation very accurately. Having performed the validation runs, we can continue increasing the polynomial order of the global surrogate for the second ring to assess how many expansion terms would have been required to improve 6
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 8. Resulting temperature distribution at t ̃ = 50 in the sections located near the first bend of the pipe. The PDFs are scaled using the minimum and maximum value of each profile individually. Therefore, the darker colour denotes local maxima of density functions.
threshold values, we generate multi-element polynomial representations. Together local surrogates reconstruct the response functions more accurately than global PCE for the same number of simulation runs. In Fig. 13 we plot the approximation error (computed assuming the interpolation of 250 URANS evaluations as true model) of solutions at t ̃ = 50 , obtained with PCE and multi-element approach, referred to as ME-PCM, for increasing number of simulation outputs used for expansion. The reconstructed functions describe the change of temperature at θ = 180° in selected rings due to varying Froude number. The multielement approach results in more efficient modelling as it adaptively divides the parametric space into elements that can be approximated with low-order basis as shown in Fig. 14. 4.3. Transient surrogate modelling At each time step, the system solutions can be evaluated only at the defined collocation points in the current set of elements. If at a given time one of the sub-domains needs to be divided, the conditions of the new element can be either directly transferred from the original element, if the random variables are equivalent, or the system has to be initialised for new set of nodal points. For transient turbulent flows, the complexity of the response function may vary substantially across the time domain. The more complex the response, the more improvement can be gained by splitting the parametric domain into elements with high regularity. To show this, we fix the threshold for decomposition to ∊t = 0.1 and the order of local bases (all of the local results are obtained with p = 2). We then allow the system to evolve from the initial flow
Fig. 10. Approximation error of PCE surrogates of the temperature in the 2nd ring constructed with increasing polynomial order p and p + 1 simulation runs.
condition and sample the temperature along the rings (see Fig. 6). Due to the transient nature of the flow, different responses are observed and, therefore, the same threshold may lead to higher number of sub-domains (more simulations required) at later instances. The comparison of approximation errors of surrogates for the first, second and third ring built with PCE and ME-PCM is shown in Fig. 15. For the first two rings,
Fig. 9. Surrogate representation of the response function for the first two rings at t ̃ = 50 at position θ = 180°. Note, red circles are the simulator outputs for which the surrogates are generated. The black line is the result of interpolated data from 250 simulations used to validate the reduced models. 7
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 11. Examples of surrogates of normalised temperature and the resulting PDFs for 3rd, 4th and 5th ring in the u-shaped pipe at t ̃ = 50 .
Fig. 12. Error in standard deviation of temperature in the 5th ring for increasing polynomial order. 8
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 13. Approximation errors of temperature response at t ̃ = 50 and position θ = 180° for given Froude number. Two methods are compared: surrogate modelling based on global expansion (PCE) and adaptive domain decomposition (ME-PCM). The improvement for the 5th ring for small number of evaluations is less substantial as the underlying function is noisy. In case of noisy results, it may be beneficial to filter the outputs before performing surrogate modelling.
the initial decomposition takes place at t ̃ = 35 and further refinement is performed at following times t ̃ = {40, 45, 50} if the last polynomial term is decaying too slowly. Although the division of parametric space performed at earlier times might not be optimal at later instances, the resulting error is still lower than in case of global approximation. For the 3rd ring we are showing the decomposition only at t ̃ = {45, 50} as there is not much change in the wall temperature at earlier times. As the study shows, the application of ME-PCM can improve the accuracy of surrogates, particularly when only small number of runs can be performed. For each local projection we use again the Gaussian quadrature rule. In the non-intrusive time-dependent analysis, the simulations are re-started for new nodal points when the domain element is split. 5. Results and discussion
t = 30 , for three In Fig. 16, we show the flow streamlines at time ~ different Froude numbers. It can be seen that the flow stratification that was observed in the LES of Skillen et al. (2020) is qualitatively well captured at Fr = 0.67 . At lower Froude numbers, the stratification is
Fig. 14. Reconstructed temperature response at θ = 180° in ring 3 for uncertain Froude number. The global surrogate is assembled from Ne = 5 approximations with p = 3. 9
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 15. Approximation errors of surrogates. The ME-PCM minimal models are constructed based on fixed threshold for convergence of polynomial terms and order of bases, p = 2 . The number of elements is equal to Ne = number of runs/(p + 1) . If Ne is increased when moving forward in time, the simulations in the new subdomains are evaluated from initial time for the new set of inputs. This is a limitation of using NISP with Gaussian quadratures. To reduce the cost, nested rules should be considered.
stronger, and hence the stratified recirculation region extends further upstream. On the other hand, at the highest Froude number under consideration, the stratified layer has almost completely disappeared by t = 30 . the time ~ Fig. 17 shows the fluid-solid interface temperature for different Froude numbers at time t ̃ = 50 . It can be seen that the temperature at the interface for the downstream vertical section is higher for lower values of Fr. The fluid in this region is unstably stratified, and hence the positive buoyancy production of turbulent kinetic energy will be enhanced at lower Froude number (since at low Fr, gravitational forces are large relative to the inertial forces). This enhanced buoyancy production leads to enhanced turbulent mixing, and hence higher average wall temperatures. In the near-horizontal section, the interface temperature towards the upper portion of the pipe (θ = 0 or 2π ) decreases with increasing Fr. This is due to the increased buoyancy effect at low Fr driving the hot transient upwards towards the upper wall (see Fig. 18). In the upstream vertical pipework, the interface temperature also decreases with increasing Fr. In this region, the thermal inertia of the wall plays a role, with the cold near-wall fluid being accelerated due to gravity in the streamwise direction, relative to the hot core fluid. This acceleration (more prominent at lower Fr since buoyancy forces are greater) draws the hotter fluid from the inlet further into the domain, leading to a thinning of both the thermal and momentum boundary layers and subsequent heating of the wall. We now consider the surrogate response for the second ring obtained with ME-PCM using Ne = 2 and Ne = 4 elements and p = 4 . In Fig. 19 we plot the Reynolds averaged temperature, as a function of Fr at locations θ = {0°, 180°} . We observe a complex response at the bottom of the ring (θ = 180°). As we saw in Fig. 16, the stratified recirculation region extends further upstream at lower Fr. If a given streamwise coordinate at θ = 180° is inside or outside the stratified recirculating flow will therefore depend on both the Froude number, and the time, t .̃ In Fig. 19b, we can see the bottom of the second ring is
t = 30 for three different Froude numbers, Fig. 16. Flow streamlines at time ~ coloured by normalised temperature (see legend in Fig. 17). In all cases, the seeding of the streamlines is at the same location.
10
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
Fig. 17. Flattened interface temperature contours (cylindrical coordinate system). Dashed lines denote locations of the two bends. The inlet is at z ′/ D = 0 . Time t ̃ = 50 .
strength of the stratified recirculation region is decreasing with increasing Fr, hence the temperature increases. In Fig. 19a, we observe that the higher gravitational forces at lower Fr cause the hot fluid to rise towards the upper wall. This causes a decreasing wall temperature at θ = 0 with increasing Fr. In Fig. 20, we plot the secondary flow and temperature distribution at various Fr and streamwise locations in the near-horizontal pipework. At low Fr, the strength of the buoyancy-induced secondary flow due to the thermal inertia of the wall is relatively high, due to the greater relative strength of the gravitational forces at lower Fr. Strong counterrotating vortex pairs have the potential to confine the hot fluid to the vortex cores, away from the wall. This effect can be observed in Fig. 20 (e.g. at Fr = 0.458 and x / D = 4.5) where entrainment of the cold stratified fluid by the secondary flow leads to a rising column of cold fluid at the core of the pipe. The subsequent conversion of kinetic energy (of the cold secondary flow) to gravitational potential energy causes this
Fig. 18. Profiles of Reynolds averaged normalised temperature at different Fr. Red – Fr = 0.458; Black – Fr = 0.670 ; Green – Fr = 0.885. Solid domain denoted with a blue colour.
inside the recirculation region (at time t ̃ = 50 ) for Fr ≲ 0.55, and outside for higher Fr. This leads to the sharp jump in temperature at this location, since being outside the recirculation region leads to higher heating of the wall. At low Fr (prior to the jump), the temperature is increasing with increasing Fr. The location at the bottom of the second ring is inside the recirculation region at low Fr, but the extent and
Fig. 19. Surrogates of temperature in the second ring constructed with ME-PCM using polynomials of order p = 4 and N2 = 2 (a) and Ne = 4 (b) elements at different locations. Time t ̃ = 50 . 11
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
t = 30, Fr ∈ {0.458, 0.67, 0.885} . Outline of the solid Fig. 20. Secondary flow and temperature distribution at cross-sections through the near-horizontal section, at ~ domain denoted with a red line.
cold fluid to sink, creating the observed tight vortex pairs. As Fr increases further, the strength of the secondary flow decreases (since gravitational forces become less significant). The transport of heat towards the upper wall is therefore reduced, and hence the temperature at θ = 0° also reduces. At sufficiently high Fr, we would expect the secondary flow to revoke back to Dean vortices, which have the opposite direction to the buoyancy-induced vortices shown in Fig. 20. Although the simulations were not run at sufficiently high Fr to observe this effect fully, the onset of Dean vortex formation can be observed in Fig. 20, at Fr = 0.885 and x / D = 4.5, where there are two pairs of vortices (one small pair of Dean vortices, and another pair due to the buoyancy induced secondary flow). The effect this has is a thickening of the thermal boundary layer around the top of the pipe as the Dean vortex secondary flow transports fluid away from the upper part of the wall.
6. Concluding remarks We have used surrogate modelling in order to enable a systematic approach for analysing thermal transients and their impact upon vortex dynamics over a broad range of Fr. The results have shown that such study can be efficiently performed with a limited number of model evaluations. This is particularly beneficial when single simulation is computationally expensive. To improve the applicability of surrogates, more sophisticated adaptive criteria for space decomposition will be considered in future work. We observe secondary flows due to the thermal inertia of the wall; the near-wall fluid heats relatively slowly, which causes a radial temperature (and hence density) gradient. This was observed by Skillen et al. (2020) for a single point in parameter space. This is the first study, to the best of our knowledge, that extends the parameter space and assesses the physics of this secondary flow over a broad range of mixed 12
Nuclear Engineering and Design 368 (2020) 110787
M.J. Zimoń, et al.
convection conditions that are of relevance to the nuclear industry. We have found that, the secondary flow has a tendency to confine the hotter core fluid to the centres of the two counter-rotating vortices, away from the pipe wall. For larger values of Fr, the opposite has been observed and the weak gravitational forces cause only small levels of flow recirculation and comparatively weak secondary flow. The polynomial chaos approach with parameter decomposition allows us to identify the most probable temperatures in time and space due to a large range of possible initial flow conditions. In some cases, the resulting distribution of temperature is bimodal in nature. The extreme values of the variable are more likely to occur than the average profile. This observation highlights misinterpretation of the system response in finite order statistics such as the mean and standard deviation.
Foo, J., Wan, X., Karniadakis, G.E., 2008. The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications. J. Comput. Phys. 227 (22), 9572–9595. Giannakopoulos, G.K., Frouzakis, C.E., Boulouchos, K., Fischer, P.F., Tomboulides, A.G., 2017. Direct numerical simulation of the flow in the intake pipe of an internal combustion engine. Int. J. Heat Fluid Flow 68, 257–268. Giles, M.B., 2013. Multilevel Monte Carlo methods. In: Monte Carlo and Quasi-Monte Carlo Methods 2012. Springer, pp. 83–103. Hellström, L.H., Zlatinov, M.B., Cao, G., Smits, A.J. Turbulent pipe flow downstream of a 90° bend. J. Fluid Mech. 735. Iyer, K.N., Theofanous, T.G., 1991. Decay of buoyancy-driven stratified layers with applications to pressurized thermal shock: reactor predictions. Nucl. Sci. Eng. 108 (2), 184–197. Kalpakli, A., Örlü, R., 2013. Turbulent pipe flow downstream a 90° pipe bend with and without superimposed swirl. Int. J. Heat Fluid Flow 41, 103–111. Knio, O., Le Maître, O., 2006. Uncertainty propagation in CFD using polynomial chaos decomposition. Fluid Dyn. Res. Kumar, R., Jadhav, P.A., Gupta, S.K., Gaikwad, A.J., 2014. Evaluation of thermal stratification induced stress in pipe and its impact on fatigue design. Proc. Eng. 86, 342–349. Manceau, R., Hanjalić, K., 2002. Elliptic blending model: a new near-wall Reynolds-stress turbulence closure. Phys. Fluids 14 (2), 744–754. McKay, M.D., Beckman, R.J., Conover, W.J., 1979. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245. Najm, H.N., 2009. Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics. Annu. Rev. Fluid Mech. 41, 35–52. Noorani, A., El Khoury, G.K., Schlatter, P., 2013. Evolution of turbulence characteristics from straight to curved pipes. Int. J. Heat Fluid Flow 41, 16–26. Owen, N.E., Challenor, P., Menon, P., Bennani, S., 2017. Comparison of surrogate-based uncertainty quantification methods for computationally expensive simulators. SIAM/ ASA J. Uncertainty Quantif. 5 (1), 403–435. Qiao, S., Gu, H., Wang, H., Luo, Y., Wang, D., Liu, P., Wang, Q., Mao, Q., 2014. Experimental investigation of thermal stratification in a pressurizer surge line. Ann. Nucl. Energy 73, 211–217. Rezende, H.C., Navarro, M.A., Mesquita, A.Z., Campagnole-dos Santos, A.A., Jordão, E., 2011. Experiments on one-phase thermally stratified flows in nuclear reactor pipe lines. Científica 15 (1), 17–24. Rezende, H.C., Santos, A.A.C., Navarro, M.A., Jordão, E., 2012. Verification and Validation of a thermal stratification experiment CFD simulation. Nucl. Eng. Des. 248, 72–81. Rütten, F., Schröder, W., Meinke, M., 2005. Large-eddy simulation of low frequency oscillations of the dean vortices in turbulent pipe bend flows. Phys. Fluids 17 (3) 035107. Sakakibara, J., Machida, N., 2012. Measurement of turbulent flow upstream and downstream of a circular pipe bend. Phys. Fluids 24 (4) 041702. Schuler, X., Herter, K.-H., 2004. Thermal fatigue due to stratification and thermal shock loading of piping. In: 30th MPA–Seminar in Conjunction with the 9th GermanJapanese Seminar, Stuttgart, Germany. pp. 6–1. Skillen, A., Zimoń, M.J., Sawko, R., Tunstall, R., Moulinec, C., Emerson, D.R., 2020. Thermal transients in a U-bend. Int. J. Heat Mass Transf. 148 119039. Sullivan, T.J., 2015. Introduction to Uncertainty Quantification, vol. 63 Springer. Vashisth, S., Kumar, V., Nigam, K.D.P., 2008. A review on the potential applications of curved geometries in process industry. Ind. Eng. Chem. Res. 47 (10), 3291–3337. Viollet, P.L., 1987. Observation and numerical modelling of density currents resulting from thermal transients in a non rectilinear pipe. J. Hydraul. Res. 25 (2), 235–261. Wan, X., Karniadakis, G.E., 2005. An adaptive multi-element generalized polynomial chaos method for stochastic differential equations. J. Comput. Phys. 209 (2), 617–642. Wan, X., Karniadakis, G.E., 2006. Multi-element generalized polynomial chaos for arbitrary probability measures. SIAM J. Scientific Comput. 28 (3), 901–928. Xiu, D., 2010. Numerical Methods for Stochastic Computations: A Apectral Method Approach. Princeton University Press. Xiu, D., Karniadakis, G.E., 2002. The Wiener-Askey polynomial chaos for stochastic differential equations. J. Scientific Comput. 24 (2), 619–644. Zimoń, M.J., Elisseev, V., Sawko, R., Antão, S., Jordan, K., 2018. Uncertainty quantification-as-a-service. In: Proceedings of the 28th Annual International Conference on Computer Science and Software Engineering. IBM Corp. pp. 331–337. Zimoń, M.J., Antão, S., Sawko, R., Skillen, A., Elisseev, V., 2019. Enabling UQ for complex modelling workflows. In: International Conference on Computational Science. Springer. pp. 269–281.
CRediT authorship contribution statement Małgorzata J. Zimoń: Methodology, Investigation, Writing - original draft. Alex Skillen: Formal analysis, Investigation, Writing - original draft. Wendi Liu: Writing - review & editing. Robert Sawko: Conceptualization, Software. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was supported by the STFC Hartree Centre’s Innovation Return on Research programme, funded by the Department for Business, Energy & Industrial Strategy. We thank Dr. Samuel Antão (IBM Research Europe) for leading the development of QUTE software that was used to generate discussed surrogate results. We acknowledge Dr. Charles Moulinec (STFC Daresbury Laboratory) for his help in modelling with Code_Saturne. References Augustin, F., Gilg, A., Paffrath, M., Rentrop, P., Wever, U., 2008. Polynomial chaos for the approximation of uncertainties: chances and limits. Eur. J. Appl. Math. 19 (2), 149–190. Carlsson, C., Alenius, E., Fuchs, L., 2015. Swirl switching in turbulent flow through 90 pipe bends. Phys. Fluids 27 (8) 085112. Couaillier, V., Savin, É., 2019. Generalized polynomial chaos for non-intrusive uncertainty quantification in computational fluid dynamics. In: Uncertainty Management for Robust Industrial Design in Aeronautics. Springer, pp. 123–141. Daly, B.J., Harlow, F.H., 1970. Transport equations in turbulence. Phys. Fluids 13 (11), 2634–2649. da Silva, L.L., Mansur, T.R., Cimini, C.A., 2011. Thermal fatigue damage evaluation of a pwr npp steam generator injection nozzle model subjected to thermal stratification phenomenon. Nucl. Eng. Des. 241 (3), 672–680. Dick, J., Kuo, F.Y., Sloan, I.H., 2013. High-dimensional integration: the quasi-Monte Carlo way. Acta Numerica 22, 133–288. EDF, Code_Saturne, 2019.https://www.code-saturne.org. Feinberg, J., Langtangen, H.P., 2015. Chaospy: an open source tool for designing methods of uncertainty quantification. J. Comput. Sci. 11, 46–57. Foo, J., Karniadakis, G.E., 2010. Multi-element probabilistic collocation method in high dimensions. J. Comput. Phys. 229 (5), 1536–1557.
13