POD-Galerkin modeling of a heated pool

Progress in Nuclear Energy 113 (2019) 196–205

Contents lists available at ScienceDirect

Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene

POD-Galerkin modeling of a heated pool

T

∗

Jorge Yáñez Escanciano , Andreas G. Class Institute for Energy and Nuclear Technique, Karlsruhe Institute of Technology, 76131, Karlsruhe, Germany

A R T I C LE I N FO

A B S T R A C T

Keywords: POD Reduced order model Sobolev Galerkin projection Incompressible fluid mechanic

Spent nuclear fuel elements contain a significant amount of fissile materials that gradually decompose generating heat and radiation. This decomposition occurs inside of deep water pools, where spent elements are cooled through natural convection. CFD calculation of all necessary cases is not feasible. Nevertheless, a light, fast and accurate model of this convection problem can be obtained utilizing Proper Orthogonal Decomposition (POD) and Galerkin Projection methods. Thus, we carry out our modeling as follows: i) Firstly, the high fidelity solver Star-CCM+ is utilized to resolve the incompressible Boussinesq formulation of the pool. The results obtained constitute a set of solutions available at discrete times for each variable. ii) Subsequently, these solutions are utilized to build a special basis, which constitutes the POD. This is built considering an optimal linear combination of the solutions at different times. Notably, each reduced set of components of the basis allows for the reproduction of a maximum of the dynamics of the variables. iii) Thirdly, the system of equations is projected in this basis, Galerkin Projection. A Reduced Order Model (ROM) can be created using a small amount of the components, sufficient for the required accuracy. This simplified model is thus mathematically sound, and derived from first principles. Finally, the results of the model are compared with high fidelity solutions. The assessment includes the capabilities of the model to reproduce transients and to approach the final steady state. Additionally, we evaluate the performance of the MPI-parallelized software generated.

1. Introduction In the nuclear industry, passive heat removal systems, as well as compact pool design reactors of IV generation (Abderrahim et al., 2001, 2010; Alemberti et al., 2011), require the accurate description of pool thermohydraulics (Grishchenko et al., 2015). This especially requires the correct and accurate calculation of the complex transient phenomena caused by natural convection. In the nuclear reactor application, because of the potential disastrous safety consequences of the results, the truthlikeness of the computations is always challenged. The verisimilitude of a study is significantly enhanced by addressing the Uncertainty Analysis, (Smith, 2013). This means considering the possible bias and the unavoidable random variation that occurs connected with the initial conditions as well as the aleatory development of the accident scenario. Typically, this is achieved by a parametric study that may involve a manifold of large dimensions. Such a task can be carried out by means of different methodologies. High Fidelity computational fluid mechanics calculations are nowadays capable of delivering a very high accuracy. Nevertheless, they are numerically intensive and subsequently inadequate to resolve the whole ∗

spectra of conditions of a parametric study. Therefore, the use of a faster methodology is required. Fast calculations are largely done using Low Fidelity simulation, ad hoc modeling and experiments. Lumped Parameter Codes (Housiadas, 2002) or one Dimensional Thermal Hydraulics codes are typical examples of Low Fidelity techniques. In the opinion of the authors, these methods are vulnerable to criticism with respect to the conceptual methodology. This raises concerns about the validity and accuracy of the whole process, especially for those cases involving the resolution of complex flow configurations. Obviously, the use of a methodology based on first principles providing reasonable accuracy with a small computational cost is necessary. Following Smith (2013) several techniques are already available. We may cite the Eigenfunctions or Modal Expansions, Regression Based Models, High Dimensional Model Representation and Surrogate-Based Bayesian Models. Among them, Reduced Order Models based on Proper Orthogonal Decomposition (POD) have developed significantly in the last twenty years. They are mathematically sound and based on first principles. Also, they allow for the calculation of complex geometries drastically reducing the computation time. To summarize, the Proper Orthogonal Decomposition (POD) methodology consists of computing some snapshots of a problem with a CFD

Corresponding author. E-mail address: [email protected] (J.Y. Escanciano).

https://doi.org/10.1016/j.pnucene.2019.01.017 Received 9 October 2018; Received in revised form 19 December 2018; Accepted 17 January 2019 0149-1970/ © 2019 Elsevier Ltd. All rights reserved.

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

finding the most energetic modes, and projecting the governing equations into the manifold generated. The potential of POD techniques as a sound, fast and general methodology for modeling has raised the attention of the applied scientific community. The excellent monographic books of Hesthaven et al. (2016) and Quarteroni et al. (2015) provide a good proof of this interest. These works provide reference and detailed description of the theory of POD and its applications to several branches of the technique. Note that POD has been already successfully implemented in the field of fluid dynamics (Rempfer, 2000). Among others, we may mention the studies of Sirovich (1987) and Aubry et al. (1988) who have investigated turbulence and the turbulent boundary layer with this technique. A more general overview on the applications of the method to fluid dynamics can be obtained in Berkooz et al. (1993) or Rathinam and Petzold (2003). In this article, we devoted ourselves to the study of natural convection in a three dimensional pool heated at its bottom, the so called Rayleigh-Bénard convection, or instability. This problem has received a lot of attention since the germinal work of Rayleigh (1916), see e.g. Chandrasekhar (2013) and Drazin and Reid (2004). Even recently, this problem was the object of some analysis utilizing POD techniques, see Herrero et al. (2013) and Pla et al. (2015). We consider this case an interesting workbench whose successful completion constitutes a first and necessary step for the usage of model order reduction in more realistic problems in the nuclear environment. We would also like to remark that it is a significant case from the point of view of the modeling of complex technical systems involving nonlinear stability analysis. The simplicity of the results obtained from POD technique can be judged observing the final model, eqs. (31) and (32),

Fig. 1. Sketch of the thermal convection system.

x i = x i* / L0 , ui = ui*/ v0 , p = p* /(ρv02) , T = (T * − T0)/(T1 − T0) . We take t0 = L02 / κ and v0 = L0 / t0 where κ is the thermal diffusivity. The nondimensional continuity, momentum, and temperature conservation equations are, following Ferziger and Peric (2012),

∇⋅→ u = 0,

(1)

⎯⎯⎯⎯⎯→ 1 → ⎯⎯⎯→ St → Δ u − ∇P + RT⎯→ ut + → u ⋅ ∇→ u = e⎯1, Re

(2)

⎯⎯⎯⎯→ StTt + → u ⋅ ∇T =

1 ΔT . RePr

(3)

The Strouhal, Reynolds and Rayleigh numbers are St = L0 /(v0 t0) , Re = ρv0 L0 / μ , Ra = ρ2 gβ (T1 − T0) L03 /(μ2 Pr) . In the last equations we have also considered for simplicity the variable R = Ra/(Re2Pr) , where Pr is the Prandtl number, μ the dynamic viscosity and g the acceleration of gravity. u ⋅→ n =0 The problem is subjected to the boundary conditions: i) → → where n is the normal vector in all boundaries; ii) T = Tup and T = Tdown in the upper and lower walls; iii) Lateral X: ∂u y / ∂x⋅n x = 0 , ∂uz / ∂x⋅n x = 0 ; iv) Lateral Z: ∂u x / ∂z⋅n z = 0 , ∂u y / ∂z⋅n z = 0 ; v) Upper and lower walls: ∂u x / ∂y⋅n y = 0 , ∂uz / ∂y⋅n y = 0 ; vi) Thermal slip boundary conditions in laterals; vii) X: ∂T / ∂x⋅n x = 0 ; viii) Z: ∂T / ∂z⋅n z = 0 .

→ → αu˙ + β→ u − νθ + ζklm ul um = 0, → → ηθ˙ + ι→ u + ϕθ + μklm ul θm = 0. which we will discuss in detail in section 5.6. Here, we mention only → that → u and θ are the amplitudes of the projected manifold in the basis of POD. They have a dimension of order 10. Greek letters are matrices or second order tensors. Their dimensions are of order 10 × 10 or of order 10 × 10 × 10 . A notable reduction of order when compared to equivalent CFD formulations, of order 105 − 107 degrees of freedom. To achieve our goals, this paper is organized as follows. Firstly, in Section 2 we address the formulation of the problem, its basic equations and boundary conditions. In Section 3 we describe POD and the Method of Snapshots. The High Fidelity computations necessary for the creation of the surrogate model are discussed in Section 4. The Galerkin projection of the system of equations is addressed in Section 5. Finally, the results obtained with our reduced order model are considered in Section 6.

2.2. Quiescent conductive solution The previous problem is known to have a simple stationary conductive solution (Drazin and Reid, 2004), c → u = 0,

Tc =

(4)

Tup − Tdown ymax

P c = P0 + R ⎜⎛ ⎝

2. Formulation of the problem

y + Tdown,

Tup − Tdown y 2 + Tdown y ⎟⎞. ymax 2 ⎠

(5)

(6)

Rewriting the problem subtracting this solution in terms of the variables θ = T − T c and Π= P − P c , we obtain

2.1. Basic formulation We consider a three dimensional parallelepiped fluid domain, Ω, of sizes L0 × L0 × 2 2 L0 confined between the two solid plates marked red in Fig. 1. The walls’ temperatures, Tup and Tdown , are chosen so that Tup − Tdown < 0 . The lateral boundaries, marked in grey, are both vertical, impenetrable and exhibit zero-stress and zero heat-flux. The problem is simplified using the Boussinesq approximation, in which density ρ is a constant, except for the effect of a vertical buoyancy force in which the density is assumed to depend linearly on temperature. ρ = ρ0 [1 − β (T − T0)], where ρ0 is the density at temperature T0 and β is the thermal expansion coefficient. In a Cartesian coordinate system, the dimensional variables (marked with an asterisk) are made non-dimensional considering t = t */ t0 ,

∇⋅→ u = 0,

(7)

⎯⎯⎯⎯⎯→ 1 → ⎯⎯⎯⎯→ ⎯e , St → Δ u − ∇Π + Rθ ⎯→ ut + → u ⋅ ∇→ u = 1 Re

(8)

⎯⎯⎯→ Tup − Tdown → ⎯→ 1 u ⋅ ⎯e1 = St θt + → u ⋅ ∇θ + Δθ , ymax RePr

(9)

with homogeneous boundary conditions canceling. This constitutes a more convenient formulation for our ultimate manipulations. 2.3. Variational formulation At this point, and because it is going to be utilized extensively, we 197

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

(→ u Ju , θ JT , Π JP )(t , x )= → J u = (∑l ul (t ) ϕl (x ), ∑l JT θl (t ) τl (x ), ∑l JP pl (t ) πl (x )),

introduce the variational formulation of the problem (7) to (9). Let us consider the space L2 (Ω) of Lebesgue square integrable functions over Ω. Let us also consider the space H∂1Ω (Ω) of Sobolev functions that verify our Dirichlet boundary conditions. We define our pseudo-velocity space V as [H∂1Ω (Ω)] × [H∂1Ω (Ω)] × [H∂1Ω (Ω)] and our pseudo-temperature space S as L2 (Ω) . The variational formulation of our problem reads, find u ∈ V and θ ∈ S such that

⎯⎯⎯⎯⎯→ 1 St∫Ω →→ v ut dΩ + ∫Ω → v ⎛⎜→ u ⋅ ∇→ u ⎟⎞ dΩ = Re ∫Ω → v Δ→ u dΩ− ⎝ ⎠ ⎯⎯⎯⎯→ ⎯e dΩ, v ∇Π dΩ + R θ→ v ⋅ ⎯→ − →

∫Ω

∫Ω

⎯⎯⎯→ St∫Ω ξθt dΩ + ∫Ω ξ→ u ⋅ ∇θ dΩ+ Tup − Tdown → ⎯→ 1 + ξ u ⋅ ⎯e dΩ =

∫Ω

1

1

→ where ϕl , τl and πl are the reduced basis vectors, ul , θl and pl their respective coefficients and ju , jT , jP the number of components to be considered for each variable.

4. High fidelity calculation. Snapshots obtained In order to obtain the database of snapshots, the commercial code STAR-CCM + version 12.02.011 has been utilized. A hexahedral regular mesh of size 2.5 mm with a total of 179,200 cells has been used for this purpose. Initial condition where that of a completely stagnant flow at 300 K. Shortly summarized, the qualitative distribution of the flow patter is as follows. In the first phase, the flow behaves quiescent. The temperature distribution is growing form the lower plate (similar to the patters described by (4)–(6)). A second phase starts after the onset of the instabilityHerrero et al. (2013). The flow starts to move. Vortexes appear and interacts which each others, reordering. The reordering develops until a steady state is reached. For more details readers are referred to a monograph work, i.e. Koschmieder (1993). In our simulation we have utilized the values L0 = 0.1m λ = 0.0269W m−1K−1, rho = 1.14kg/m−3 , cp = 1000Jkg−1K−1, ν = 2.35⋅10−4m s−2 , β = 1⋅10−5K−1, g = 74m s−2 , Ra = 39990 and R = 398710 . A total of 1000 s were simulated in order to fully achieve a steady state. In Section 3 it has been shown that the reduced order variables are located inside of the span of the snapshots, namely → u Ju ∈ span(u (t1), …, u (tn)) , Π JT ∈ θ JT ∈ span(θ (t1), …, θ (tn)) , span(Π(t1), …, Π(tn)) . The sampled snapshots determine the span of the solution of the surrogate model. In this work we have sampled snapshots at a constant interval of time. Thus, we face the problem of the selection of the sampling frequency. Note that the sampling process constitutes a kind of low-pass filter, excluding any event occurring at a time scale faster than the sampling rate. For its selection, in this paper we have employed the following strategy. The reduced order variables are linear combinations of the base vectors that in turn are linear combinations of the snapshots. Therefore, a sampling frequency which provides consecutive snapshots that are linear or collinear–up to a threshold–cannot be enriched collecting more snapshots. Successive refinement of the sampling frequency have allowed us to find a frequency at which this condition is fulfilled. The variance inflation factor, VIF, see Backhaus et al. (2015), provides a criterion to quantify the degree of multicollinearity employing an ordinary least squares regression analysis. In other words, whether a snapshot can be expressed, up to an arbitrary threshold, is determined as a linear combination of the other snapshots. For our analysis, we have selected the coefficient of determination 0.9999, equivalent to a VIF = 10000 , i.e. an extremely high collinearity. This allows us to assure that the snapshot database is satisfactory. The corresponding sampling rate of 10 Hz was found to satisfy our criteria of multicollinearity for speed, pressure and temperature.

(10)

∫ ξ ΔθdΩ,

(11) → for all test functions v ∈ V and ξ ∈ S verifying boundary conditions. ymax

RePr Ω

3. Proper orthogonal decomposition. The method of snapshots For the calculation of the reduced base, we utilize the Method of Snapshots. This was introduced by Berkooz et al. (1993) and Sirovich (1987) to study turbulence in incompressible flows. Let us consider the discrete solutions of the system (7)–(9) in a discrete set of times {t1, t2, …, tn} for example obtained from a CFD. These solutions are called snapshots. We consider the matrix formed by the discrete solutions at the set of times, considering one or several variables coupled as Y = [y (t1, x ), y (t2, x ), …, y (tn, x )] where y is any of the variables considered. We utilize as (Quarteroni et al., 2015, chapter 6) the so called Singular Value Decomposition.

V T YG

⎛ σ1 0 = Σ= ⎜ ⎜⎜ ⋮ ⎝0

0 σ2 ⋮ 0

… 0⎞ … 0⎟ , ⋮ ⎟⎟ … σn ⎠

(12)

to obtain a reduced order basis. Here, V, G are orthogonal matrices and Σ contains on its diagonal the singular values σi . The projection in one dimension less produces an error of order ‖σn ‖. This provides the possibility of discarding an arbitrary number of dimensions, j − n , producing a total error ‖σj ‖2 + ‖σj + 1 ‖2 +…+ σn 2 . The reduced system carries a relative information index j n I = ∑i = 1 σi 2 / ∑i = 1 σi 2 . Obtaining the basis vectors is simple if we consider the so called matrix of covariances Y T Y . We obtain Y T Y = G Σ2GT , an eigenvalue problem,

Y T YG = GΣ2 ,

(13)

σi2

where the eigenvalues are and the eigenvectors are the vectors of G. Coming back to Y = V ΣGT , we can resolve using Moore-Penrose inverse V = YGΣ−1 to get the vectors of the basis seek

vi = Ygi / λi .

(15)

(14)

Notably, the vectors constructed in this way (Quarteroni et al., 2015) represent the closest manifold to the data of the original problem with a certain dimension. Fixing a threshold for the relative information index, the number of necessary reduced basis vectors for each of the variables is obtained. At this stage, the adequateness of the system of equations (7)–(9) is evident. As we seek to provide the solution as a linear combination of snapshots, homogeneous boundary conditions allow for a simplification of the boundary treatment. Any linear combination of the snapshots satisfies the boundary conditions. Assuming the variables are affine in t and → x , we expand the variables of the system in terms of the reduced basis, obtaining the reduced order variables

5. Galekin projection 5.1. Basic formulation We continue the narrative of Section 3. We reformulate the variational formulation (10)–(11) in the reduced order basis. That is find → u Ju ∈ span(ϕ1, …, ϕ Ju ) , θ JT ∈ span(τ1, …, τJT ) such that 198

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

⎯⎯⎯⎯⎯⎯→ St∫Ω → v Ju → u Ju t dΩ + ∫Ω → v Ju ⎜⎛→ u Ju⋅ ∇→ u Ju ⎟⎞ dΩ= ⎝ ⎠ 1 = Re ∫Ω → v Ju Δ→ u Ju dΩ− ⎯⎯⎯⎯→ ⎯e dΩ, − ∫Ω → v Ju ∇Π dΩ + R∫Ω θ JT → v Ju⋅ ⎯→ 1 ⎯⎯⎯⎯⎯⎯→ St∫Ω ξ JT θ JT t dΩ + ∫Ω ξ JT → u Ju⋅ ∇θ JT dΩ+ Tup − Tdown → ⎯→ 1 + ∫Ω ξ JT y u Ju⋅ ⎯e1 dΩ = RePr ∫Ω ξ JT Δθ JT dΩ, max

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

(16)

(17)

v Ju ∈ span(ϕ1, …, ϕ Ju ) and all ξ JT ∈ span(τ1, …, τJT ) . for all → → In the Galerkin formulation, v Ju and ξ JT take the form of each of the vectors of the reduced basis.

Foias et al. (1988), found that any system of ordinary differential equations derived from the Galerking projection of a dissipative partial differential equation system appears to be unstable for long term dynamics. This is the case of the reduced model of equations (16) and (17), which exhibit such spurious dynamics (Rapún and Vega, 2010). Although the ultimate cause is still a subject of controversy (Terragni et al., 2011), it appears that lack of sufficient dissipation is the most probable origin. Several strategies have been developed to overcome these difficulties, which restore an adequate level of dissipation. Sirisup and Karniadakis (2004) suggested the inclusion of additional terms to create additional viscosity. Couplet et al. (2005) utilized POD-Galerkin system with modified coefficients. Siegel et al. (2008) utilized a two level proper orthogonal decomposition, in which each mode of the first level can be expressed in terms of the second one. Rapún and Vega (2010), Terragni et al. (2011) and Terragni and Vega (2014) mixed toghether CDF and POD calculation steps to update the ROM basis along the time integration. In this work we follow Kirby (1992), and Iollo et al. (2000) and define the POD in the Sobolev Space H1 (Adams and Fournier, 2003) rather than in L2 . In H1 we define the inner product as

y 1

x

⎯u (t ) … ε ∂y ⎯→ y 1 ⎯→ ⎯ ε ∂z u y (t1) … ⎯→ ⎯u (t ) … z 1

⎯u (t ) … ε ∂x ⎯→ z 1 ⎯→ ⎯ ε ∂y uz (t1) … ⎯u (t ) … ε ∂ ⎯→ z 1

z

(21)

5.4. Variational formulation in H1 In H1, the basic variational formulation of the problem contains numerous additional terms. The system of equations reads

⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯ ⎯ dΩ + Stε ⎯∇ ⎯→dΩ v u v : ∇ ⎯→ u St∫Ω →⎯→ ∫Ω → t t ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ ⎯ ⎯⎯⎯⎯ → ⎯ ⎯⎯⎯ → ⎯⎯⎯⎯⎯→ + ∫Ω → v ⋅⎛⎜→ u ⋅ ∇→ u ⎞⎟ dΩ + ε∫Ω ∇→ v : ∇ ⎛⎜→ u ⋅ ∇→ u ⎞⎟ dΩ= ⎝ ⎠ ⎝ ⎠ ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯→ → 1 1 → v ⋅Δ→ u dΩ + ε ∇→ v : ∇Δ u dΩ

∫

Re Ω

Re

∫Ω

⎯⎯⎯⎯→ ⎯⎯⎯⎯→ − ∫Ω → v ∇Π dΩ − ε∫Ω ∇→ v : ⎯⎯⎯⎯→ → →⎯→ ⎯ + R∫Ω θ v e1 dΩ + Rε∫Ω ∇ v :

⎯⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯→ ∇ ∇Π dΩ ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯e → ∇ (θ ⎯→ 1 ) dΩ,

Ω

(18)

+ ∫Ω

, where ε ≥ 0 is a parameter to be determined.

Ω

Tup − Tdown

+ ε∫Ω 5.3. Matrix of snapshots under Sobolev formulation

ymax

Tup − Tdown

1 RePr Ω

ymax

⎝

⎠

⎯e dΩ+ ξ→ u ⋅ ⎯→ 1 ⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯e → u ⋅ ⎯→ ∇ξ ⋅ ∇ (→ 1 ) dΩ=

∫ ξ Δθd Ω+

1 ε RePr Ω

⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯→

∫ ∇ξ ⋅ ∇Δθ dΩ,

Under this formulation, we may form three matrices of snapshots, for each of the variables

where we have dropped Ju and JT sub-indexes.

→ ⎛ θ (t1) ⎜ → ⎜ ε ∂x θ (t1) → ⎜ ε ∂y θ (t1) ⎜ ⎜ ε∂ → z θ (t1) ⎝

5.5. Modified Sobolev formulation

→ ⎛ Π (t1) ⎜ → ε ∂x Π (t1) ⎜ → ⎜ ε ∂y Π (t1) ⎜⎜ → ⎝ ε ∂z Π (t1)

… … … …

… … … …

→ θ (tn ) ⎞ ⎟ → ε ∂x θ (tn ) ⎟ , → ε ∂y θ (tn ) ⎟ ⎟ → ε ∂z θ (tn ) ⎟⎠

→ Π (tn ) ⎞ → ⎟ ε ∂x Π (tn ) ⎟, → ε ∂y Π (tn ) ⎟ ⎟ → ε ∂z Π (tn ) ⎟⎠

(22)

⎯⎯⎯→ ⎯⎯⎯⎯→ St∫Ω ξθt dΩ + Stε∫Ω ∇ξ ⋅ ∇θt dΩ ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ ⎯⎯⎯→ ⎯⎯⎯→ ⎯⎯⎯→ + ∫ ξ→ u ⋅ ∇θ dΩ + ε∫ ∇ξ ⋅ ∇ ⎜⎛→ u ⋅ ∇θ ⎟⎞ dΩ

⎯⎯⎯→ ⎯⎯⎯→

∫Ω fgdΩ + ε∫Ω ∇f ⋅ ∇g dΩ

x 1

y

⎯u (t ) … ε ∂z ⎯→ x 1 ⎯→ ⎯u (t ) … y 1 ⎯u (t ) … ε ∂ ⎯→

⎯→ ⎯u (t ) ⎞ x n ⎯u (t ) ⎟ ε ∂x ⎯→ x n ⎯u (t ) ⎟ ε ∂y ⎯→ x n ⎟ ⎯u (t ) ⎟ ε ∂z ⎯→ x n ⎯→ ⎯u (t ) ⎟ y n ⎟ ⎯u (t ) ε ∂x ⎯→ x n ⎟ ⎯u (t ) ⎟ ε ∂y ⎯→ x n ⎟ ⎯→ ⎯ ε ∂z u y (tn ) ⎟ ⎯→ ⎯u (t ) ⎟ z n ⎯u (t ) ⎟ ε ∂x ⎯→ x n ⎟ ⎯u (t ) ⎟ ε ∂y ⎯→ x n ⎯u (t ) ⎟ ε ∂z ⎯→ z n ⎠

and subsequently execute the rest of the content of Section 3. We have opted for a formulation in which the three components of velocity are lumped together while the rest of the variables are considered independently.

5.2. Sobolev Space formulation

〈f , g〉ε =

⎯→ ⎯u (t ) … x 1 ⎯→ ⎯ ε ∂x u x (t1) … ⎯u (t ) … ε ∂ ⎯→

(23)

Substitution of the simplifications of Appendix A into equations (22) and (23) yield, after some algebra,

⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯ ⎯ dΩ + Stε ⎯∇ ⎯→dΩ St∫Ω →⎯→ v u v ⋅ ∇ ⎯→ u ∫Ω → t t ⎯⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯⎯⎯⎯⎯ → ⎯ ⎯⎯⎯ → ⎯ ⎯⎯⎯⎯ → ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯→ 1 1 + Re ∫Ω ∇→ v : ∇→ u dΩ + Re ε∫Ω ∇ ∇→ v : ∇ ∇→ u dΩ ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯e → − R∫Ω θ→⎯→ v ⎯e1 dΩ − Rε∫Ω ∇→ v : ∇ (θ ⎯→ 1 ) dΩ=

(19)

⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯→ − ∫Ω → v ⋅⎛⎜→ u ⋅ ∇→ u ⎟⎞ dΩ − ε∫Ω ∇→ v ⋅ ∇→ u ⋅ ∇→ u dΩ ⎝ ⎠ → ⎯⎯⎯⎯→ → ⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯⎯⎯⎯⎯→ − ε∫Ω ∇→ v ⋅ u ⋅ ∇ ∇→ u dΩ,

(20)

199

(24)

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

⎯⎯⎯→ ⎯⎯⎯⎯→ St∫Ω τθt d Ω+ Stε∫Ω ∇τ ⋅ ∇θt dΩ + ∫Ω

Tup − Tdown

+ ε∫Ω +

ymax

1 RePr Ω

(

⎯e dΩ τ→ u ⋅ ⎯→ 1

Tup − Tdown ymax

→ (St A + St εB ) u˙ + → 1 1 + Re B + Re εS → u − (RC + RεG ) θ =

⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯e → ∇τ ⋅ ∇ (→ u ⋅ ⎯→ 1 ) dΩ

⎯⎯⎯→ ⎯⎯⎯→

∫ ∇τ ⋅ ∇θ d Ω +

1 ε RePr Ω

)

= − (Ωklm + εψklm + εχklm ) ul um , → (St H + St εL) θ˙ +

→ ⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯⎯⎯ ⎯⎯⎯→ ⎯⎯⎯→

∫ ∇ ∇τ : ∇ ∇θ dΩ=

⎯⎯⎯→ ⎯⎯⎯→ ⎯⎯⎯⎯⎯→ ⎯⎯⎯→ − ∫Ω τ→ u ⋅ ∇θ dΩ − ε∫Ω ∇τ ⋅ ∇→ u ⋅ ∇θ d Ω ⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯→ ⎯⎯⎯→ − ε∫Ω ∇τ ⋅→ u ⋅ ∇ ∇θ dΩ.

(29)

+

Tup − Tdown ymax

(J + εM ) → u +

1 (L RePr

→ + εN ) θ =

= − (Kklm + εγklm + εOklm) ul θm . (25)

(30)

Collecting terms, we obtain the very compact formulation

→ → αu˙ + β→ u − νθ + ζklm ul um = 0,

(31)

→ → ηθ˙ + ι→ u + ϕθ + μklm ul θm = 0.

(32)

5.6. Matrix analogous formulation The L2 inner products of two arbitrary variables ℵ1, ℵ2 is approximated

5.7. Numerical details

n

∫Ω ℵ1ℵ2 dΩ ≈ ∑ ℵ1 (xi) ℵ2 (xi).

(ℵ1, ℵ2) =

(26)

i

The calculation of the terms of Table 1 involves the first and second derivatives of the reduced base vectors. Taking into account the simple geometry of Fig. 1, an appropriate way of calculating the derivatives could be through the use of the Fourier Transform, as carried out by Kirby (1992). Nevertheless, as we intend to apply our development to complex geometries with an arbitrary shape, we have calculated the first and second derivatives with a second order finite differences scheme applicable for arbitrarily spaced grids, see Fornberg (1988). At the boundaries, we lower the order of approximation to first order. The system of equations (31) and (32) is integrated employing a Crank-Nicolson method (Endre and Mayers, 2003). Expressing ⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯→ → → ⎯→ ⎯ ⎯→ ⎯ ut + 1 = δu + ut and θt + 1 = δθ + θt for a discrete time increment δt we obtain,

The sum is extended to all nodes n coming from the high fidelity calculation. We simultaneously consider the set of all test functions. Each corresponds to a vector of the base. We denote the grouped vectors of the velocity and temperature-reduced base 0

0

⎛Q ⎞ Q = ⎜ Q1 ⎟ = ⎜ 2⎟ ⎝Q ⎠

⎛ ϕ0 1 ⎜ ϕ0 ⎜ϕ 2 ⎝ 0

… ϕ J0u

⎞

… ϕ J1u ⎟, … ϕ J2u

⎟ ⎠

(27)

Z = ( τ0 … τJT ).

(28)

The processing of each of the terms of the equations (24) and (25) provide the matrix, or tensor equivalent terms described in Table 1. In this table the symbol ⊗ denotes the outer product. Their substitution in equations (24) and (25) yield equations (29) and (30). Table 1 Equivalent terms. Comma index notation denotes derivative. Integral term

Equivalent terms

v ⋅→ ut dΩ ∫Ω →

Akl = ∑i Qki T Qli

⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯ ⎯→dΩ v : ∇ ⎯→ u ∫Ω ∇→ t → ⎯⎯⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯ → → u dΩ ∫Ω ∇ ∇ v : ∇ ∇→ θ→⎯→ v ⎯e dΩ

Bkl = ∑i ∑j Qki ,Tj Qli, j

∫Ω

1

⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯e → v : ∇ (θ ⎯→ ∫Ω ∇→ 1 ) dΩ ⎯⎯⎯⎯⎯→

Gkl = ∑j Qk1, j Zl, j

∫Ω

ψklm = ∑n ∑i ∑j ⎜⎛ϕki , n , ϕl,jn ⊗ ϕmi , j ⎟⎞ ⎝ ⎠

⎯⎯⎯⎯⎯⎯⎯⎯⎯→

→ ⎯⎯⎯⎯⎯→ → ⎯⎯⎯⎯⎯→ v ⋅ u ⋅ ∇ ∇ u dΩ ∫Ω ∇→

∫Ω τθt dΩ ⎯⎯⎯→ ⎯⎯⎯⎯→

∫Ω ∇τ ⋅ ∇θt dΩ ⎯e dΩ u ⋅ ⎯→ ∫Ω τ→ 1 ⎯⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯→ ⎯⎯⎯→ ∫Ω ∇ ∇τ : ∇ ∇θ dΩ ⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ ⎯e → u ⋅ ⎯→ ∫Ω ∇τ ⋅ ∇ (→ 1 ) dΩ ⎯⎯⎯→ u ⋅ ∇θ dΩ ∫Ω τ→ ⎯⎯⎯→ ⎯⎯⎯⎯⎯→ ⎯⎯⎯→ u ⋅ ∇θ dΩ ∫Ω ∇τ ⋅ ∇→ ⎯⎯⎯⎯⎯⎯⎯⎯→ ⎯⎯⎯→ → ⎯⎯⎯→ ∫Ω ∇τ ⋅ u ⋅ ∇ ∇θ dΩ

→ ⎛ −1 ⎛ →t ⎞ α ⎜−βu + νθt − ζklm ult umt⎞⎟ ⎜ ⎟ δt ⎝ ⎠ − 2⎜ ⎟. →t →t t t ⎞⎟ ⎜ η−1 ⎛⎜−ιu ϕθ μ u θ − − klm l m⎟ ⎟ ⎜ ⎝ ⎠⎠ ⎝

χklm = ∑n ∑i ∑j ⎛⎜ϕki , n , ϕl j ⊗ ϕmi , nj ⎞⎟ ⎝ ⎠ Hkl = ZkT Zl

(34)

5.8. Pressure equation

Lkl = ∑j ZkT, j Zl, j Jkl =

(33)

⎯⎯⎯⎯⎯→ ⎛ −1 ⎛ ⎯⎯⎯⎯⎯→ ⎞ α ⎜−β ut + 1 + ν θt + 1 − ζklm ult + 1 umt+ 1⎞⎟ ⎜ ⎟ δt ⎝ ⎠ → en + 1 = 2 ⎜ ⎟ → ⎯⎯⎯⎯⎯→ t + 1 t + 1⎞ ⎟ t + 1 − ϕ θt + 1 − μ ⎜ η−1 ⎛⎜−ι ⎯⎯⎯⎯⎯ u u θ m ⎟⎟ klm l ⎜ ⎝ ⎠⎠ ⎝

Ckl = Qk1 T Zl

Ωklm = ∑i ∑j ⎛⎜ϕki , ϕl j ⊗ ϕmi , j ⎞⎟ ⎝ ⎠

⎝ ⎠ → ⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯→ v ⋅ ∇ u ⋅ ∇→ u dΩ ∇→

δt

−2ν

For the step n + 1, the integration error is,

Skl = ∑i ∑j ∑m Qki ,Tjm Qli, jm

v ⎜⎛→ u ⋅ ∇→ u ⎟⎞ dΩ ∫Ω →

δt (β + ζklm umt + ζklm ult ) 2 δt (ι + μklm θmt) 2

⎯ ⎞ ⎛ ⎯→ δu ⎞ ⎟ ⎯ ⎟+ ⎜ ⎯→ δt η + 2 (ϕ + μklm ult ) ⎟ ⎝ δθ ⎠ ⎠ → → ⎛−βut + νθt − ζklm ult umt ⎞ δul δum ⎞ δt ζ + 2 ⎛⎜ klm ⎟ = δt ⎜ ⎟. →t →t ⎜ −ιu t t ⎟ ⎝ μklm δul δθm ⎠ ϕθ μ u θ − − l m klm ⎠ ⎝ ⎛α + ⎜⎜ ⎝

Zk Ql1

Pressure depends on both the magnitude of velocity and temperature. Its amplitudes in the reduced base can be calculated explicitly as a function of the amplitudes of the modes of velocity and temperature. We start from the pressure differential equation,

Nkl = ∑i ∑j ZkT, ij Zl, ij

Mkl = ∑j ZkT, j Ql1, j K klm = ∑i (τ k , ϕli ⊗ τm, i )

:

⎯⎯⎯→ ⎯⎯⎯⎯⎯→ ΔΠ = R→ e1⋅ ∇θ − ∇→ u

ϒklm = ∑n ∑i (τ k, n, ϕli, n ⊗ τm, i ) Oklm = ∑n ∑i (τ k, n, ϕli ⊗ τm, ni )

⎯⎯⎯⎯⎯ → ∇→ u.

(35)

L2 (Ω) .

We define our pseudo-pressure space X as The variational formulation of equation (35) corresponds to find Π∈ X so that 200

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

→ →

∫Ω hΔΠdΩ + ε∫Ω ∇ h⋅∇ ΔΠdΩ= → → → ⎯ → ⎯e ⋅∇ R [∫Ω h⎯→ θd Ω+ ε∫Ω ∇ h⋅∇ (⎯→ e y ⋅∇ θ) dΩ]− y

:

:

⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯ → → → → ⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯ ⎡ ⎤ − ⎢∫Ω h ∇→ u ∇→ u dΩ + ε∫Ω ∇ h⋅∇ ⎛⎜ ∇→ u ∇→ u ⎞⎟ dΩ⎥. ⎝ ⎠ ⎦ ⎣

(36)

for all test functions h ∈ X . Taking into account the simplifications of Appendix A, after some algebra we obtain

⎯⎯⎯→ ⎯⎯⎯⎯→ − ∫Ω ∇h ∇Π dΩ − ε∫Ω ΔhΔΠdΩ= → → → ⎯ → ⎯e ⋅∇ R∫Ω h⎯→ θdΩ + Rε∫Ω ∇ h⋅∇ (⎯→ e1 ⋅∇ θ) dΩ 1 ⎯⎯⎯⎯⎯ → ⎯ ⎯⎯⎯⎯ → ⎯ ⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯→ − ∫Ω h ∇→ u ∇→ u dΩ + ε∫Ω Δh ∇→ u ∇→ u dΩ.

:

:

(37)

We use W to denote the matrix formed by the vectors of the reduced base, in order to simultaneously consider all vectors of the reduced base,

W = ( π0 … π JP ).

Fig. 2. Successive values of the ε obtained in the minimization procedure.

(38)

We process each of the terms of equation (37) in order to obtain a discrete formulation. The resulting matrix and tensor terms are contained in Table 2. Substituting the terms of the table into equation (37) we obtain,

procedure itself was carried out utilizing Brent's method, see Atkinson (2008). The convergence of the method has been depicted in Fig. 2 resulting in a value of ε = 1.95⋅10−4 .

→ * (E + εE *)→ p = −R (F + εF *) θ + (Λklm − εΛklm ) ul um .

6.2. Modes of proper orthogonal decomposition. Analysis

(39)

And collecting terms,

→ E ′→ p = −F ' θ + Λ′klmul um .

The theory of Sections 3 and 5.3 was applied to the database of snapshots described in Section 4. It is concluded that for a relative information index, I, of 99.9%, 13 Temperature, 9 Velocity and 3 Pressure modes should be taken into account. We call this our Low POD Model. We have depicted the first ten temperature modes in Fig. 3. In Fig. 4 we have depicted the evolution of the four first temperature modes compared with the values obtained from projecting the results of the high fidelity calculation into the same reduced basis. We observe a high degree of accuracy with some minor deviations. The dynamics of the problem are well reproduced. Also the quasi-steady state that is practically reached at 0.18 dimensionless seconds is well captured. The method provides thus a successful simulation. A second, larger, relative information index threshold has also been considered. This allows, see Terragni et al. (2011) and Terragni and Vega (2012), the estimation of the spatial relative error associated with the Galerkin approximation. This estimate is available apriori in contrast to the method suggested by Quarteroni et al. (2015) and does not require knowledge about the discretization schemes of the high fidelity solver. This knowledge obviously is inaccessible if a commercial CFD solver is utilized. The method is based on the well known fact (Hesthaven et al., 2016) (Quarteroni et al., 2015) that the magnitude of the amplitudes of the PODs exponentially decreases with the increasing number of PODS. This can also be seen in Fig. 5. The magnitude of the amplitudes of the modes contained between two sufficiently separated information thresholds is thus far a fairy good estimator of the error in the case that we were only using the Low POD Model,

(40)

6. Results of the numerical experiment 6.1. Minimization of Sobolev parameter In the theory contained in Sections 5.3 to 5.8 we have described our problem as a function of a free parameter ε. This parameter determines not only the results of the integration of the system of equations (31) and (32), but it also plays its role in the finding of the POD modes. ε can be determined by minimizing the distance between the amplitudes of the modes of the reduced order model and the equivalent amplitudes of the high fidelity calculation (HFC) in the same reduced base JT

n

l

i=1

min ⎧∑ ∑ |θl (ti ) − θlHFC (ti )|⎫.

ε ∈ [0,1]

⎨ ⎩

⎬ ⎭

(41)

After several trials, it was found that most accurate results were obtained when the sole amplitudes of a single variable, temperature or velocity, were taken into account. Other formulations, considering combination of variables proved unpractical. The minimization Table 2 Equivalent terms for pressure equation. Integral term

⎯⎯⎯→ ⎯⎯⎯⎯→

∫Ω ∇h ∇Π dΩ ∫Ω ΔhΔΠdΩ ⎯⎯⎯→

e1⋅ ∇θ dΩ ∫Ω h→ → →

→

⎯e ⋅∇ θ) dΩ ∫Ω ∇ h⋅∇ (⎯→ 1

: :

⎯⎯⎯⎯⎯→ u ∫Ω h ∇→

⎯⎯⎯⎯⎯→

u ∫Ω Δh ∇→

n

Equivalent terms

⎯⎯⎯⎯⎯→ → ∇ u dΩ ⎯⎯⎯⎯⎯→ → ∇ u dΩ

Enn1 =

Ekl = WkT, j Wl, j

Ekl* = ∑i ∑j WkT, ii Wl, jj

∑ j 1= n + 1 aj2 n

∑ j 1= 1 aj2

, (42)

where the aj are the amplitudes of an arbitrary variable, and n and n1 are the number of modes associated with the first and second relative information index threshold respectively. Terragni et al. (2011) found that seeking an error of the relative information index, 1 − I , a hundred times smaller than the original one, (1 − I )/100 , provides consistent results for the estimation of the error associated with the Galerkin approximation. This resulted into 42 Temperature, 26 Velocity and 13 Pressure modes. Some of the

Fkl = WkT Zl,1

Fkl* = ∑i WkT, i Zl, i1 Λklm = ∑i ∑j ⎛⎜πk , ϕl,ji ⊗ ϕmi , j ⎞⎟ ⎝ ⎠

Λ*klm = ∑n ∑i ∑j ⎛⎜πk, nn, ϕl,ji ⊗ ϕmi , j ⎞⎟ ⎝ ⎠

201

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

Fig. 4. Evolution of the amplitudes of the first four temperature modes in time compared with the values resulting from the CFD high fidelity calculation.

Fig. 5. Normalized eigenvectors of Temperature and velocity.

frequency fluid structures than higher ones. This suggests that higher modes contain a larger amount of numerical noise. This provably originates from the CFD calculation and the calculation of the derivatives. Therefore, lower modes are physically more significant than higher modes, that may even lack any physical meaning. Thus, high Enn1 values imply high numerical error.

Fig. 3. Isosurfaces of the first ten dimensionless temperature POD modes.

7. Concluding remarks We have described a reduced order model which is mathematically sound, and derived from first principles. Our construction involves the solution of a non-linear model consisting of low dimensional matrices and second order tensors. Integration of the ROM can be accomplished within a fraction of the time needed for the high fidelity simulation, so that several orders of speed increments can be achieved. Our creation is especially intended for uncertainty qualification and is suitable for problems in which massive analysis of multiple initial conditions is required. In this sense we may shortly analyze the performance of our construct.

amplitudes of the corresponding modes are shown in Fig. 5. We have depicted the values of Enn1 in Fig. 6. The aprioristic estimation remains under 5% for most of the calculation. The comparatively large relative error of the velocities at the beginning of the simulation is due to the low velocity values at this stage. More significant are the deviations located between 0.02 and 0.04 dimensionless seconds. Those correspond to the error in the peak values and the ulterior phase. We can also interpret the two thresholds’ approximation in terms of the shape of the POD modes, see Fig. 3. Lower modes contain lower 202

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

This involves the calculation of the reduced basis generating the POD modes and the computing of the terms appearing in Table 1. With the very extensive database of this problem with a sampling rate of 10 Hz and considering 300 s of simulated time this phase involve roughly 14 min in 32 processors. (b) The integration of the ROM system of equations, eqs. (31) and (32), requires 100 s. (c) Optimization phase of ε parameter. This requires the completion of points 7 and 7 several times. Some savings by storage can be carried out leaving each iteration in 9 min. In the naive optimization of Fig. 2, 16 iterations where necessary totaling 2.4 h. By simply choosing a lower initial value, e.g. ε = 0.001 this phase can be strongly reduced. (d) Now the reduced order model is fully operational. Any new calculation, for an eventual parameter study, etc. will require 300 s. We intend to expand the methodology to more complex and applied systems including heat sources, turbulence, coupled equations, generic boundary conditions and complex geometries.

Fig. 6. Relative error.

1. High fidelity CFD calculations involve approximately 3 days calculation. 2. The calculation POD can be divided into tasks: (a) The calculation of the matrices of the model eqs. (29) and (30).

Acknowledgement We kindly give thanks to the Framatome Professional School and the European Union Horizon 2020 SESAME project.

Appendix A. Simplification of the formulation We seek for some simplification of the equations (22) and (23). Firstly, we realize that pressure terms in equations (22) and (23) cancel. This is due to the boundary conditions and to the soleinoidal velocity field. (a)

⎯⎯⎯⎯→

n dΩ − ∫ Π∇⋅→ v dΩ. ∫Ω →v ⋅ ∇Π dΩ = ∫∂Ω Π→v ⋅→ Ω

(A.1)

u ⋅→ n = 0, → v ⋅→ n = 0 in ∂Ω. Therefore, Due to the homogeneous Dirichlet boundary conditions and because → n = 0. ∫∂Ω Π→v ⋅→

(A.2)

→ Also, as → v is a linear combination of snapshots it is soleinoidal. Thus, ∇⋅ v = 0 , so that ⎯⎯⎯⎯→

∫Ω →v ⋅ ∇Π d Ω= 0.

(A.3)

(b)

→ ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯⎯ ⎯⎯⎯⎯→

∫Ω ∇→v : ∇ ∇Π dΩ = ∫Ω ∂x vi ∂2x x ΠdΩ. j

(A.4)

j i

Here → v is soleinoidal. Thus,

∂2xj xi vi = ∂ xj (∂ xi vi) = 0.

(A.5)

So that

∫Ω ∂xj vi ∂2xj xi ΠdΩ = ∫Ω ∂xi (∂xj vi ∂xj Π) dΩ= = ∫∂Ω ni ∂ xj vi ∂ xj ΠdΩ= = ∫∂Ω ∂ xj (ni vi ∂ xj Π) dΩ − ∫∂Ω ni vi ∂2 2 ΠdΩ. x

(A.6)

j

Moreover the boundary conditions ni vi = 0 imply that

→ ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯⎯ ⎯⎯⎯⎯→

∫Ω ∇→v : ∇ ∇Π d Ω= 0.

(A.7)

Secondly, we seek for a decrease of the order of derivatives thus minimizing the needed manipulation of the snapshot data. The original data was generated by a finite volume solver utilizing no more that first order derivatives. Therefore, it is advantageous to reduce derivative order whenever possible. (c)

203

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

⎯⎯⎯⎯⎯→

⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯→

u dΩ = ∫ ⎛→ v ⋅ ∇→ u ⎞⋅→ n dΩ − ∫ ∇→ v : ∇→ u dΩ. ∫Ω →v ⋅Δ→ ∂Ω Ω

∫Ω

⎜

⎟

⎝

⎠

(A.8)

⎯⎯⎯⎯⎯→ → v ⋅ ∇→ u ⎟⎞⋅ n ≡ vi nj ∂ xj ui , and boundary conditions imply that nj ∂ xj ui |∂Ω = 0 . Thus, Since ⎛⎜→ ⎝ ⎠ ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯→ → v ⋅Δ→ u dΩ = − ∇→ v : ∇→ u dΩ.

∫Ω

(A.9)

(d)

⎯⎯⎯→

⎯⎯⎯→ ⎯⎯⎯→

∫Ω τ ΔθdΩ = ∫Ω ∇⋅⎛τ ∇θ ⎞ dΩ − ∫Ω ∇τ ⋅ ∇θ dΩ. ⎜

⎟

⎝

⎠

(A.10)

Here,

⎯⎯⎯→

⎯⎯⎯→

n d Ω= 0, ∫Ω ∇⋅⎛τ ∇θ ⎞ dΩ = ∫∂Ω τ ∇θ ⋅→ ⎜

⎟

⎝

⎠

(A.11)

due to the boundary conditions. Therefore, (e)

⎯⎯⎯→ ⎯⎯⎯→

∫Ω τ ΔθdΩ = −∫Ω ∇τ ⋅ ∇θ dΩ.

(A.12)

⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯→ → v : ∇Δ u dΩ = ∫Ω ∂ xk vi ∂ xk ∂2xj xj ui dΩ= ∫Ω ∇→ = ∫∂Ω nj ∂ xk vi ∂2xk xj ui dΩ − ∫Ω ∂2xk xj vi ∂2xk xj ui dΩ.

(A.13)

The first term,

∫∂Ω nj ∂x vi ∂2x x ui dΩ = ∫∂Ω ∂x (nj ∂x ui) ∂x vi dΩ. k

j

k

k j

(A.14)

k

The slip boundary conditions imply that nj ∂ xj ui = 0 . Thus

∫∂Ω nj ∂x vi ∂2x x ui dΩ = 0. k

(A.15)

k j

Finally,

⎯⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯→

⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯⎯→

u dΩ = −∫ ∇ ∇→ v : ∇ ∇→ u dΩ. ∫Ω ∇→v : ∇Δ→ Ω

(A.16)

(f)

⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯→

∫Ω ∇τ ⋅ ∇Δθ dΩ = ∫Ω ∂xi τ ∂xi ⎜⎛∂2xj xj θ⎟⎞ dΩ=

⎝ ⎠ = ∫δ Ω ∂ xi τ ∂ xi (nj ∂ xj θ) dΩ − ∫Ω ∂2xi xj τ ∂2xi xj θdΩ.

(A.17)

But from the boundary conditions, nj ∂ xj θ . Thus,

→ ⎯⎯⎯⎯⎯⎯⎯ → ⎯⎯⎯⎯⎯⎯⎯ ⎯⎯⎯→ ⎯⎯⎯→

⎯⎯⎯→ ⎯⎯⎯⎯⎯⎯→

∫Ω ∇τ ⋅ ∇Δθ dΩ = −∫Ω ∇ ∇τ : ∇ ∇θ dΩ.

(A.18)

We simplify equation (36) analogously to the procedure shown for velocity and temperature. We look for a decrease in the order of the deri⎯⎯⎯→ ⎯⎯⎯⎯→ vatives. Additionally, to ∫Ω hΔΠdΩ = −∫Ω ∇h ∇Π dΩ we also consider: (i)

→ →

∫Ω ∇ h⋅∇ ΔΠdΩ = ∫Ω ∂xn h∂xn (∂ xi2 Π) dΩ= = ∫∂Ω nn ∂ xn h∂ xi2 ΠdΩ − ∫Ω ∂2x 2 h∂2x 2 ΠdΩ. n

(A.19)

i

In Galerkin Formulation, h is any of the πl . At any time instance this can be expressed as a linear combination of dimensionless homogeneous pressure snapshots, but the boundary conditions of the snapshots verify ni ∂ xi Π(tk )|∂Ω = 0 . Therefore

→ →

∫Ω ∇ h⋅∇ ΔΠdΩ = −∫Ω ∂2x

2

2 h∂ 2 ΠdΩ. xi n

(A.20)

(ii)

:

→ → → ⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯ ∇→ u ⎟⎞ dΩ ≡ ∫Ω ∂ xn h∂ xn (∂ xi uj ∂ xj ui ) dΩ ⎝ ⎠ = ∫Ω nn ∂ xn h∂ xi uj ∂ xj ui dΩ − ∫Ω ∂2 2 h∂ xi uj ∂ xj ui dΩ.

u ∫Ω ∇ h⋅∇ ⎜⎛ ∇→

(A.21)

xn

Following the same reasoning (i), we obtain

204

Progress in Nuclear Energy 113 (2019) 196–205

J.Y. Escanciano, A.G. Class

:

→ → → ⎯⎯⎯⎯⎯→ ⎯⎯⎯⎯⎯ ∇→ u ⎟⎞ dΩ = − ⎜ ⎝ ⎠

u ∫Ω ∇ h⋅∇ ⎛ ∇→

∫Ω ∂2x

2 h∂ xi uj ∂ xj ui dΩ. n

(A.22)

Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.pnucene.2019.01.017.

Housiadas, C., 2002. Lumped parameters analysis of coupled kinetics and thermal-hydraulics for small reactors. Ann. Nucl. Energy 29 (11), 1315–1325. Iollo, A., Lanteri, S., Désidéri, J.-A., 2000. Stability properties of pod–galerkin approximations for the compressible Navier–Stokes equations. Theor. Comput. Fluid Dynam. 13 (6), 377–396. Kirby, M., 1992. Minimal dynamical systems from pdes using sobolev eigenfunctions. Phys. Nonlinear Phenom. 57 (3–4), 466–475. Koschmieder, E.L., 1993. Bénard Cells and Taylor Vortices. Cambridge University Press. Pla, F., Herrero, H., Vega, J.M., 2015. A flexible symmetry-preserving galerkin/pod reduced order model applied to a convective instability problem. Comput. Fluids 119, 162–175. Quarteroni, A., Manzoni, A., Negri, F., 2015. Reduced Basis Methods for Partial Differential Equations: an Introduction, vol. 92 Springer. Rapún, M.-L., Vega, J.M., 2010. Reduced order models based on local pod plus galerkin projection. J. Comput. Phys. 229 (8), 3046–3063. Rathinam, M., Petzold, L.R., 2003. A new look at proper orthogonal decomposition. SIAM J. Numer. Anal. 41 (5), 1893–1925. Rayleigh, L., 1916. Lix. on convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. The London, Edinburgh Dublin Philosoph. Magaz. J. Sci. 32 (192), 529–546. Rempfer, D., 2000. On low-dimensional galerkin models for fluid flow. Theor. Comput. Fluid Dynam. 14 (2), 75–88. Siegel, S.G., Seidel, J., Fagley, C., Luchtenburg, D., Cohen, K., McLaughlin, T., 2008. Lowdimensional modelling of a transient cylinder wake using double proper orthogonal decomposition. J. Fluid Mech. 610, 1–42. Sirisup, S., Karniadakis, G., 2004. A spectral viscosity method for correcting the long-term behavior of pod models. J. Comput. Phys. 194 (1), 92–116. Sirovich, L., 1987. Turbulence and the dynamics of coherent structures. i. coherent structures. Q. Appl. Math. 45 (3), 561–571. Smith, R.C., 2013. Uncertainty Quantification: Theory, Implementation, and Applications, vol. 12 Siam. Terragni, F., Valero, E., Vega, J.M., 2011. Local pod plus galerkin projection in the unsteady lid-driven cavity problem. SIAM J. Sci. Comput. 33 (6), 3538–3561. Terragni, F., Vega, J.M., 2012. On the use of pod-based roms to analyze bifurcations in some dissipative systems. Phys. Nonlinear Phenom. 241 (17), 1393–1405. Terragni, F., Vega, J.M., 2014. Construction of bifurcation diagrams using pod on the fly. SIAM J. Appl. Dyn. Syst. 13 (1), 339–365.

References Abderrahim, H.A., Baeten, P., De Bruyn, D., Heyse, J., Schuurmans, P., Wagemans, J., 2010. Myrrha, a multipurpose hybrid research reactor for high-end applications. Nucl. Phys. News 20 (1), 24–28. Abderrahim, H.A., Kupschus, P., Malambu, E., Benoit, P., Van Tichelen, K., Arien, B., Vermeersch, F., Dhondt, P., Jongen, Y., Ternier, S., et al., 2001. Myrrha: a multipurpose accelerator driven system for research & development. Nucl. Instrum. Methods Phys. Res. Sect. A Accel. Spectrom. Detect. Assoc. Equip. 463 (3), 487–494. Adams, R.A., Fournier, J.J., 2003. Sobolev Spaces, vol. 140 Academic press. Alemberti, A., Carlsson, J., Malambu, E., Orden, A., Cinotti, L., Struwe, D., Agostini, P., Monti, S., 2011. Elsyeuropean lfr activities. J. Nucl. Sci. Technol. 48 (4), 479–482. Atkinson, K.E., 2008. An Introduction to Numerical Analysis. John Wiley & Sons. Aubry, N., Holmes, P., Lumley, J.L., Stone, E., 1988. The dynamics of coherent structures in the wall region of a turbulent boundary layer. J. Fluid Mech. 192, 115–173. Backhaus, K., Erichson, B., Plinke, W., Weiber, R., 2015. Multivariate analysemethoden: eine anwendungsorientierte einführung. Springer-Verlag. Berkooz, G., Holmes, P., Lumley, J.L., 1993. The proper orthogonal decomposition in the analysis of turbulent flows. Annu. Rev. Fluid Mech. 25 (1), 539–575. Chandrasekhar, S., 2013. Hydrodynamic and Hydromagnetic Stability. Courier Corporation. Couplet, M., Basdevant, C., Sagaut, P., 2005. Calibrated reduced-order pod-galerkin system for fluid flow modelling. J. Comput. Phys. 207 (1), 192–220. Drazin, P.G., Reid, W.H., 2004. Hydrodynamic Stability. Cambridge university press. Endre, S., Mayers, D., 2003. An Introduction to Numerical Analysis. (Cambridge, UK). Ferziger, J.H., Peric, M., 2012. Computational Methods for Fluid Dynamics. Springer Science & Business Media. Foias, C., Jolly, M., Kevrekidis, I., Sell, G., Titi, E., 1988. On the computation of inertial manifolds. Phys. Lett. 131 (7–8), 433–436. Fornberg, B., 1988. Generation of finite difference formulas on arbitrarily spaced grids. Math. Comput. 51 (184), 699–706. Grishchenko, D., Jeltsov, M., Kööp, K., Karbojian, A., Villanueva, W., Kudinov, P., 2015. The tall-3d facility design and commissioning tests for validation of coupled sth and cfd codes. Nucl. Eng. Des. 290, 144–153. Herrero, H., Maday, Y., Pla, F., 2013. RB (reduced basis) for RB (Rayleigh–Bénard). Comput. Methods Appl. Mech. Eng. 261, 132–141. Hesthaven, J.S., Rozza, G., Stamm, B., et al., 2016. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. Springer.

205