Dynamic mode decomposition for the stability analysis of the Molten Salt Fast Reactor core

Nuclear Engineering and Design 362 (2020) 110529

Contents lists available at ScienceDirect

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

Dynamic mode decomposition for the stability analysis of the Molten Salt Fast Reactor core

T

Andrea Di Roncoa, Carolina Introinia, Eric Cervia, Stefano Lorenzia, Yeong Shin Jeongb, ⁎ Seok Bin Seob, In Cheol Bangb, Francesca Giacobboa, Antonio Cammia, a b

Politecnico di Milano, Department of Energy, Via La Masa 34, 20156 Milano, Italy Ulsan National Institute of Science and Technology, Department of Nuclear Energy, Ulsan, South Korea

A R T I C LE I N FO

A B S T R A C T

Keywords: Dynamic mode decomposition Reduced order models Molten salt fast reactor Stability analysis

The study of innovative nuclear reactors involves the use of increasingly complex numerical models. While such models provide a high-fidelity description of many non-linear coupled phenomena, they are not suited for manyquery tasks such as design optimisation, uncertainty quantification, stability analysis or parameter identification due to the required computational effort. For this reason, a variety of techniques have been employed to reduce the complexity and surrogate the response of large nuclear systems. One example is the dynamic mode decomposition (DMD), a data-driven method which builds a low-dimensional eigenvalue-eigenvector representation of the underlying model from numerical data, and allows for non-intrusive analyses of the dynamical properties of the system without knowledge of the model itself. In this work, DMD is applied to the study of a free-dynamics fast transient of the Molten Salt Fast Reactor (MSFR), following a variation of the heat transfer coefficient. The numerical data is provided by a multiphysics model developed using the open-source CFD toolkit OpenFOAM. The aim of this work is to demonstrate the applicability of DMD to the study of large next-generation nuclear systems such as the MSFR. The results show the capabilities of DMD to extract and surrogate the dynamics of the MSFR following perturbation, including the initial non-linear dynamics and the final steadystate. Different values of parameters relevant to the construction of DMD models are tested, to provide some insights on the sensitivity of the method to the selection of the numerical data set and to the size of the reduced model.

1. Introduction Over the last few years, reactor modelling has shifted towards largescale, full-core time-dependent analysis because of the ever-increasing safety requirements. However, the accurate study of the performance of systems with such complex geometries and features poses a great challenge both to numerical simulations, computational algorithms and physical experiments. Despite advancements in computational resources, this will not likely be possible for all those applications in which a large amount of model runs must be performed. It is the case, for example, of design optimisation, parameter identification or uncertainty quantification (i.e. the so-called many-query analyses). Linked to the above, ever-growing importance has been given in the last years to the dynamic monitoring of nuclear reactors for near realtime analysis. This analysis should be robust and quick to identify potential critical situations ahead of time. Such monitoring requires fast and reliable methods to identify and track the evolving dynamics of

⁎

critical system modes (Hauer et al., 2007). For these so-called manyquery analysis and for the dynamic monitoring, relatively inexpensive and yet accurate enough surrogate models are useful tools. Starting from a high-fidelity model and corresponding solutions, or from experimental data, data-driven techniques can generate low-fidelity surrogates without the need of additional simplifications and approximations of the physics of the system, and to extract dominant patterns exhibited by transient processes. The rationale behind the development of model order reduction (MOR) techniques is the existence of low-dimensional coherent structures in physical phenomena (Antoulas, 2005). This allows, in principle, the reduction of the model complexity by reducing its degrees of freedom (DOF). Large models such as computational fluid-dynamics (CFD) and finite element models (FEM) can benefit from reduced order approaches. Reduced order and metamodelling techniques have been extensively employed to analyse and surrogate large-scale models in several fluid dynamics applications, including some in the field of

Corresponding author. E-mail address: [email protected] (A. Cammi).

https://doi.org/10.1016/j.nucengdes.2020.110529 Received 4 July 2019; Received in revised form 7 January 2020; Accepted 21 January 2020 0029-5493/ © 2020 Elsevier B.V. All rights reserved.

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

Nomenclature

A Ã A b bik ci c1̂ k erel , max m m̃ N n Q q‴ q ‴̂ r T T̂

t tk U W Xm 0 x z(t ) zk z jk z (̂ t ) z k̂ z jk̂ Δt Δts λi ϕi ̃ Ψ ψi ̃

Linear operator for the given system ∈ n × n Reduced order matrix of the system ∈ r × r Approximation of the linear operator using the SVD ∈ n × n Vector of the amplitude coefficients bi Amplitude coefficient of the i-th mode at time tk Concentration of neutron precursors of the i-th group Reconstructed concentration of neutron precursors of the 1st group Maximum relative reconstruction error at time tk Number of available snapshots Index of the last snapshot used to build the DMD model Size of the full order computational mesh Size of the full order system state Total reactor thermal power Volumetric power density Reconstructed volumetric power density Truncation value for the SVD decomposition Temperature Reconstructed temperature

Σ σi ωi

Time instant Discrete time instant of the k-th snapshot −1 Matrix of the economy-size SVD decomposition of X m 0 m−1 Matrix of the economy-size SVD decomposition of X 0 Matrix for the sequence of snapshots ∈ n × (m + 1) Space coordinate vector System state at time t k-th snapshot of the system state (t = tk ) j-th component of the k-th snapshot of the system state DMD approximation at time t DMD approximation at time tk j-th component of the DMD solution at time tk Time step of the snapshots Sampling time step of the snapshots i-th eigenvalue of the reduced matrix A ̃ i-th eigenvector of the reduced matrix A ̃ Matrix of the dynamical modes i-th modal structure extracted from A ̃ Matrix of the economy-size SVD decomposition of X1m − 1 i-th singular value of the SVD decomposition of X1m − 1 i-th “continuous-time” eigenvalue of the reduced matrix A ̃

the presented DMD method to the MSFR case study is presented. Several surrogate models are built based on different reduced model sizes and sampling strategies. Results from the different surrogate models are then compared and discussed on the basis of the approximation error introduced (with respect to the original data) and the time-dependent approximation of the total (i.e. volume-integrated) reactor thermal power during the power transient. Conclusions are finally given in Section 5, together with some considerations regarding possible future developments.

nuclear reactor analysis. Among popular techniques, proper orthogonal decomposition (POD) was used for example to study the coolant pool of a Lead Fast Reactor (LFR) from a control-oriented viewpoint (Lorenzi et al., 2017) and for the parametric multiphysics analysis of the Molten Salt Fast Reactor (MSFR) (German et al., 2019). A study of control rod movement by means of the reduced basis method can be found in Sartori et al. (2016). Other examples of relevant applications of metamodelling techniques include generalised polynomial chaos expansion (GPCE) (Wu and Kozlowski, 2017; Perkó et al., 2014) and high-dimensional model representations (HDMR) (Hu et al., 2014) for sensitivity and uncertainty analysis, machine learning (ML) (Bao et al., 2019; Chang and Dinh, 2019), sparse grid stochastic collocation (SGSC) (Bennett et al., 2019; Wu et al., 2017) and Gaussian processes (GP) (Wu et al., 2018; Marrel et al., 2015; Ikonomopoulos et al., 2013). Dynamic mode decomposition (DMD) is another model order reduction technique which found recent application in works related to nuclear reactor analysis, e.g. dealing with isotope composition evolution (Abdo et al., 2019), decay-heat generation (Alfonsi et al., 2018), pulsed neutron experiments in subcritical systems (Hardy et al., 2019) and time-dependent neutron transport (McClarren, 2019). Interest in DMD has grown over the past years, mainly thanks to its ability to surrogate complex temporal dynamics based only on observed data. The algorithm was first proposed for the analysis of experimental fluid dynamics data (Schmid, 2010), but it was quickly adapted in a variety of fields. A comprehensive description of the method can be found in Tu et al. (2014). DMD is based on singular value decomposition (SVD), and it can extract information from both simulated and experimental data. The method is data-driven in the sense that it does not require knowledge of the underlying governing equations, and only time snapshots of the system state are needed. In this work, the construction of surrogate models based on the DMD method is presented, to approximate time- and space- dependent quantities of large models with a reduced number of DOF. The physical system of interest is the MSFR. What makes the analysis and simulation of such a system challenging is the peculiar coupling between neutronics and thermal–hydraulics, with dynamics acting on very different time scales. The paper is organised as follows. In Section 2 the DMD method is presented, with some brief references to its theoretical foundations and a detailed description of the most common algorithmic implementation. The MSFR system and the adopted full order modelling approach is described in Section 3. In Section 4 the application of

2. The Dynamic Mode Decomposition method Modal decomposition techniques describe the system state as a linear combination of empirical state vectors, also called modes. Depending on the problem at hand, the number of basis vectors required to capture the dominant dynamics of the system can be several orders of magnitudes smaller than the original size. In particular, DMD techniques seek the linear dynamic operator that best approximates the underling dynamics of the system. Because of its theoretical foundation, the DMD algorithm can also characterise non-linear dynamics (Tu et al., 2014). Based on the use of either empirical or simulated ’snapshots’, this approach shares some similarities with other decomposition methods, such as the Proper-Orthogonal Decomposition (POD) (Holmes et al., 2012). Regarding the latter, DMD seeks an eigenvalue-eigenvector decomposition in which the spatial modes (the eigenvectors) are built based on their dynamics properties and thus are associated with a dynamic spectrum (the eigenvalues). The advantages of this approach are several. For instance, DMD can extract dynamically relevant structures regardless of their energy content. DMD also allows an equation-free modelling framework, since it does not require any projection of the modes to get time-dependent coefficients. 2.1. Methodology Given a dynamical system of state size n, the available data is organised in the form of a temporal sequence of snapshots in a matrix X m 0: 0 1 m Xm 0 = {z , z , …, z }

(1)

where each zk = z(tk ) ∈ n represents the sampled state vector at time instant tk = k Δt , with k = 0, …, m . The sampling time step Δt needs to be constant for the correct implementation of the DMD method. In 2

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

dimension, from m to r. The subscript res shows the truncated and thus neglected m − r singular values, whereas the superscript ∗ denotes the complex conjugate transpose. Since Σ ̃ is a square matrix, the MoorePenrose pseudo-inverse can be evaluated. A possible way to define r is to study the weight of each of the eigenvalues found by the SVD de−1 composition of the matrix X m . 0 An approximation for the linear operator A can now be calculated as:

addition, its selection is a delicate step in the algorithm definition because it directly influences the information and behaviour that the sampled data can capture. More details about the issues related to the choice of the sampling time are presented SubSection 2.2. Assuming that there is a linear mapping A that connects the state zk taken at time tk to a subsequent one zk + 1 at time tk + 1: (2)

zk + 1 = Azk

If the process is non linear, the mapping A refers to a linear tangent approximation. In both cases, this mapping is taken as constant over the entire data time range [0, mΔt ]. The purpose behind DMD is to extract the dynamical characteristics of the process described by A from the sequence in X m 0 . As the number of snapshots increases, it is reasonable to assume that, after a certain number of snapshots, the sequence of data will capture all the dominant characteristics of the physical process. The addition of further information in terms of additional snapshots will not improve the vector space generated by X m 0 , and therefore the vectors given by Eq. 2 become linearly dependent. The linear mapping A can also be expressed alternatively as: −1 † A = X1m (X m ) 0

̃ ̃ Ũ A ≈ A = X1mWΣ −1

(5)

A represents an approximation of the linear mapping of Eq. 2. An eigenvalue analysis of this matrix produces the dynamic modes and eigenvalues of the system. However, if m ≫ 1, this analysis is computationally very expensive. A much more compact and efficient analysis can be obtained by taking r ≪ m , which amounts to a projection of zk on a linear subspace of r using the general basis transformation Pzk = z k̃ , where P is the projection matrix. A convenient transformation has already been calculated through the SVD for the sequence of data, which results in P = U .̃ The projected reduced order model can now be derived:

(3)

∗ ∗ ̃ −̃ 1 z k̃ = Az̃ k̃ z k̃ + 1 = U ̃ AUz̃ k̃ = U ̃ X1mWΣ

The † superscript represents the Moore-Penrose pseudo-inverse, which is found through the following singular value decomposition on the sequence of data: ∗ Σ ̃ 0 ⎤ ⎡ W̃ ⎤ −1 ̃ ̃ ̃∗ ̃ ]⎡ = UΣW∗ = [ U ̃ Ures Xm ∗ ⎥ ≈ UΣW 0 ⎢ ⎢ 0 Σres ⎥ ̃ res ̃ W ⎣ ⎦⎣ ⎦

∗

(6)

By The reduced order model is given by the matrix A ̃ = solving the eigenvalue problem for the matrix A ,̃ the modal structures can be extracted as follows: ∗ ̃ −̃ 1. U ̃ X1mWΣ

̃ i ̃ = ϕi λ̃ i Aϕ

(4)

ψi ̃ =

The matrix U contains the proper orthogonal modes of the data sequence, i.e., the right singular vectors of the snapshot sequence. The above operation thus amounts to a projection of the linear mapping A onto a POD basis, and it allows for a more robust calculation of the lowdimensional representation of the linear operator. In addition, rankdeficiency in the data sequence is taken into account by restricting the projection basis U . The non-zero singular values of the matrix Σ , or above a prescribed threshold, determine this restriction. By defining an appropriate truncation value r of the singular values, it is possible to reduce in an efficient and simple way the data

1 m ̃ −̃ 1 ̃ X1 WΣ ϕi λi

(7) (8)

Since the matrix A ̃ has size r × r , with r ≪ m , this decomposition is expected to be efficient. The discrete system eigenvalues λi can now be mapped on to the corresponding continuous-time ones as follows (Zhu et al., 1985):

ωi =

lnλi Δt

(9)

This step is required to reconstruct for each time instant tk > 0 the

Fig. 1. Molten Salt Fast Reactor concept layout. 3

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

Table 1 Main MSFR design parameters.

Table 2 DMD parameters adopted in the work.

Parameter

Value

Parameter

Values

Nominal thermal power Fuel inlet temperature Fuel outlet temperature Total salt volume Fuel composition (% mol)

3000 MW 650 °C 750 °C 18 m3 LiF (77.5) – ThF4 (20.0) – 233UF4 (2.5)

Reduced model size, r Sampling time step, Δts Last snapshot index, m ̃

10, 25, 50 Δt , 4Δt , 10Δt 500, 750, 1000

general solution: r

z^ (t ) =

∑ i=1

bi0 ψi (x) exp (ωi t ) = Ψ·diag( exp (ωt ))·b0

(10)

bi0

is the initial amplitude of each mode, Ψ is the matrix of dywhere namical modes, whose columns contains the eigenvectors ψi , diag(exp(ωt )) is a diagonal matrix whose elements are the exponential of the eigenvalues, and b is the vector of the amplitude coefficients bi . To evaluate these coefficients, taking z0 as the initial snapshot, the initial amplitude coefficients are evaluated as:

b0 = Ψ †z0

(11)

where Ψ † denotes the Moore-Penrose pseudo-inverse. The algorithm described above requires taking measurements or collecting consecutive snapshots at a constant time interval. This restriction can be easily lifted using a slightly different formulation of the DMD method, which uses consecutive generic data pairs {(x 0 , y0 ), …, (x m , ym )} . Such data is gathered in matrices whose columns have the same index for each pair of data. This method and some alternative implementations of the DMD algorithm are discussed in more detail in Tu et al. (2014). This procedure allows the extraction of the coherent structures Ψ from a sequence of data fields, without knowing the linear mapping A This detachment from the underlying model represents the great advantage of DMD. Since the projection is based only on the snapshots, restricted sub-domains can be processed with this technique by forming a data sequence X m 0 with data only from a smaller part of the full computational domain. For example, the analysis can be restricted on smaller regions where dynamically interesting phenomena are expected or observed. Another application of DMD is the analysis of spatial, rather than temporal, numerical stability. The snapshots are now organised in a

Fig. 2. Structure of the OpenFOAM solver (adapted from Cervi et al., 2019).

Fig. 3. Left: simplified MSFR geometry (heat exchanger and pump region is highlighted in red); right: computational mesh. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 4

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

10 0

10 -2

10 -4

10 -6 0

20

40

60

80

100

−1 Fig. 4. First 100 (normalised) singular values of the SVD decomposition of X m (all available snapshots are considered). The index i refers to the decreasing order, 0 with σ1 being the largest singular value.

spatial sequence, equivalent to taking the transpose matrix of X m 0 . The generated eigenvalues in this case contain information about the spatial dynamics of the process.

numerical solver used in the present work. Originally developed for the transient analysis of the MSFR (Aufiero et al., 2014), it was recently extended to allow for the study of compressibility effects during superprompt-critical transients (Cervi et al., 2019) and the analysis of helium bubbling systems (Cervi et al., 2019; Cervi et al., 2019). At each time step, thermal–hydraulics and neutronics equations are solved in two different cycles, as shown in Fig. 2. The first is based on the standard solver buoyantBoussinesqPimpleFoam, whereas the latter implements the multi-group diffusion equation for the neutron flux and the transport equations of the delayed neutron precursors. Picard iterations are performed until the solution of the thermal–hydraulic part of the problem reaches convergence. Then, the neutronics cycle begins, and once the flux and the decay heat are known, the volumetric power source field is updated and the energy equation is solved again. The newly evaluated fuel temperature and density fields are used to update the cross sections, and the cycle is repeated until convergence is reached. External iterations between the two sub-solvers are also performed. As described in Section 4, the quantities of interest are the power density field q‴ (x, t ) , the fuel temperature T (x, t ) and the concentration of the first neutron precursors group c1 (x, t ) , namely the one with shortest decay time. The multiphysics solver is used to simulate a fast transient in the MSFR, which represents the full order model from which snapshots of the variables of interest are extracted. In particular, the capability of the DMD method to capture both the initial fast dynamics of the transient and the final steady-state are assessed by varying the algorithm parameters, namely the truncation value r, the number of snapshots m and the sampling time Δts . A simplified axisymmetric 2D geometry (a wedge of amplitude 5°) is adopted (Fig. 3) (Fiorina et al., 2014). The mesh counts 22671 hexahedral elements. At t = 0 s the system is in steady-state conditions at a nominal power of about 3000 MWth. The heat transfer coefficient between the reactor and the heat sink is then reduced by 25%, whilst keeping the mass flow rate constant. The response of the system to this variation is simulated. The neutron flux is approximated using six energy groups. The simulation is carried out for 25 s, with a time step of 0.001 s. Snapshots are saved each 0.025 s. To quantify the computational burden of the model, a full order run requires roughly 100 CPU-hours.

2.2. Criterion for the sampling time A fundamental parameter in the DMD algorithm is the sampling time of the snapshots from the data. To extract relevant dynamic information from the system processes, these must be sampled with a sufficiently high frequency. The Nyquist criterion identifies the lower bound for this sampling frequency as at least twice the inherent frequency of the processes of interest (Diniz et al., 2002). A common practice, which gives satisfactory results, is to adopt a sampling frequency of about three times the Nyquist cutoff. By tuning the sampling frequency, fast or slow processes can be discriminated. In fact, slow sampling frequencies samples fast processes at near-random amplitudes, interpreting them as noise. For high sampling rates, slow processes will be almost steady and will not show neither growth/decay nor any oscillatory component. Instead, they are reflected in the mean mode, or a very slow drift mode. 3. Full order modelling approach This section briefly describes the selected case study. The Molten Salt Fast Reactor (MSFR), whose concept layout is shown in Fig. 1 is the system under considered for analysis. Table 1 lists some of its main design parameters (Gerardin et al., 2017). The salt mixture, which is simultaneously fuel and coolant, circulates through sixteen external circuits, where the power produced by fission is removed. The circulating fuel is the main feature of the MSFR, and its employment leads to new design and technological challenges. In particular, the delayed neutron precursors do not decay where they are produced, but they are transported by the fuel mixture flow through the primary circuit. Delayed neutrons can be produced in the peripheral regions of the core or even outside the reactor core, with consequent reactor time constants which depend on the size and geometry of the primary circuit. Therefore, the coupling between system neutronics and thermal fluid-dynamics is even stronger than in traditional reactors because the velocity field directly influences the precursors distribution. The multiphysics open-source toolkit OpenFOAM (Weller et al., 1998; Jasak et al., 2007), based on standard finite-volume methods for the discretization of model equations, was used to develop the 5

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

3000 2900 2800 2700 2600 2500 2400 2300 0

5

10

15

20

25

15

20

25

15

20

25

(a)

3000 2900 2800 2700 2600 2500 2400 2300 0

5

10

(b)

3000 2900 2800 2700 2600 2500 2400 2300 0

Fig. 5. Solution fields of q‴ (top), c1 (middle) and T (bottom) at t = 25 s (end of the simulated transient). Thanks to the symmetries of the problem, the equations are solved only on a 5° wedge.

5

10

(c) Fig. 6. Total reactor thermal power: original transient and DMD reconstruction with different values of r and for Δts = 0.025 s (a), Δts = 0.1 s (b) and Δts = 0.25 s (c).

4. Application 4.1. Setup of the Dynamic Mode Decomposition method

more fields can be considered; however, as the different dynamics are strictly coupled, most of the dynamic information is in most cases redundant. The full order state size n is therefore 68013. It should be remarked that the specific order in which single state variables are included in the state vector is irrelevant, as long as it is kept the same for all the snapshots. Furthermore, even though different fields are processed together as a single state vector, after model reduction and state reconstruction by means of the DMD the single fields can be separated back. In this way, the reconstructed fields q ‴̂ , c1̂ and T ̂ are obtained. Since different variable fields have values in very different ranges,

In this section the application of the DMD technique to the case study presented in Section 3 is described. With Δt equal to 0.025 s, all the snapshots (except for the initial conditions) zk for tk ∈ [0.025, 25] are collected; m = 1000 snapshots are therefore available. The state vector z ∈ n is defined as the concatenation of one or more of the variable fields included in the snapshots, where n is N times the number of scalar fields considered. If vector fields are also considered, they can be included by concatenating the single scalar components as separate scalar fields. In this case, to capture the dominant dynamics of the system, 3 scalar fields are included: the volumetric power q‴, the first group of neutron precursors c1 and the temperature T. In principle, 6

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

Fig. 8. Total thermal power approximation (a) and maximum relative reconstruction error (b) of the DMD reduced order models with r = 50, Δts = Δt and different values of m ̃ (Table 2). Vertical coloured lines in (a) represent the last snapshot used to build the DMD model in each case.

Fig. 7. Maximum relative reconstruction error of the original (a) q‴, (b) c1 and (c) T fields for the r and Δts values listed in Table 2.

each quantity should be suitably normalised before manipulation to provide more homogeneous numeric values. This is done to avoid numerical conditioning issues and also to deal with non-dimensional quantities. In this case, each separate field is divided by the respective maximum value (across all snapshots and all mesh points). Having selected the state variables, a DMD reduced model is built by choosing the size of the reduced model r and the snapshots to be used to build the model. As DMD requires equally time–spaced snapshots, this is specified by means of a sampling time step Δts and the index of the last snapshot used, m .̃ To test the capabilities of the DMD approach in retrieving the dominant dynamics of the MSFR, different values are employed. A summary of the adopted values is listed in Table 2.

Fig. 9. DMD discrete-time eigenvalues λ i with r = 25 modes and m ̃ = 1000. Different colours identify the different sampling frequencies adopted in the work. The unit circle represents the stability region of a discrete-time system.

7

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

Fig. 10. Real parts of the first 4 pairs of spatial modes ψi ,̃ relative to the q‴ field. They are presented in pairs since they the correspond to pairs of complex conjugate eigenvalues. The modes for the q‴ field are simply extracted from ψi ̃ following the definition of the state vector.

The r and Δts values from Table 2 are combined to test the reduced model dynamics with different numbers of retained modes and sampling frequencies. In these cases, snapshots are sampled along the whole transient (m ̃ = 1000). The total thermal power Q, defined as

4.2. Results and discussion The original q‴, c1 and T fields at t = 25 s are shown in Fig. 5. The core region, i.e. the zone of the primary circuit where most of heat generation occurs, is clearly visible. Fig. 5 also shows the effects of recirculation due to the non-optimal design of the reactor shape in the computational model used in this work. Fluid stagnation leads to increased temperatures in the recirculation zone, even though the heat generation is less relevant than in the central region. A qualitative measure of how much reducible is the full order system is given by the −1 singular values σi of the snapshot matrix X m . Rapidly decreasing 0 singular values mean that, from an energy content viewpoint, most information is retained in the first modes. Fig. 4 shows that the singular values fall by several orders of magnitude after few modes.

Q (t ) =

∫V q‴ ⎛x, t ⎞ dV , ⎜

⎟

⎝

⎠

(12)

is used to provide an overall picture of the dynamics of the reduced DMD models compared to the original data. The maximum relative error k erel , max = max j∈J

z jk̂ − z jk z jk

(13)

is given to provide a measure of how accurately the DMD 8

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

Fig. 11. Real parts of the first 4 pairs of spatial modes ψi ,̃ relative to the c1 field. They are presented in pairs since they the correspond to pairs of complex conjugate eigenvalues. The modes for the c1 field are simply extracted from ψi ̃ following the definition of the state vector.

incorrect approximations of the initial dynamics, meaning that a sampling time of 0.25 s is not enough to resolve the fastest time scales. k The reconstruction error erel , max curves (Fig. 7) show that the approximation is in general very good. The maximum error values are of the order of few percents even in worst cases. General trends in the curves can be observed, with accuracy which worsens as r is decreased or Δts is increased. Larger errors are found during the initial part of the transients, as expected. The worst results are observed in the reconstruction of q‴ (Fig. 7a), where values around 10–20% are initially found with Δts = 10Δt . Besides the sampling frequency, it is interesting to assess how many snapshots are needed to capture the final steady-state contribution, i.e. whether or not the full order simulation can be stopped at some point and the previously computed time steps used to propagate the state of  the system with reduced computational effort. In Fig. 8 the usual Q, Q

reconstruction can match the original snapshot data across all the computational domain. To give more consistent results, the index subset J is used to select the different q‴, c1 and T fields separately. Fig. 6 shows the behaviour of Q in the different cases. Even though the DMD method fairly reproduces the overall transient (e.g. the initial undershoot or the final steady-state solution) in all cases, the cases with r = 10 introduce significant disturbances in the power dynamics. Large oscillations are observed during the whole transient. While with r = 25 some discrepancies are still observed, 50 modes are enough to closely match the original data. Regarding the sampling time step, the original snapshot time step Δt is sufficient to capture all dynamic time scales, even during the initial the part of the transient where faster and nonlinear dynamics occurs (Fig. 6a). With Δts = 4Δt (Fig. 6b) the approximation is still good, but some discrepancies start to show at the beginning of the transient. The last cases (Fig. 6c) with Δts = 10Δt show 9

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

Fig. 12. Real parts of the first 4 pairs of spatial modes ψi ,̃ relative to the T field. They are presented in pairs since they the correspond to pairs of complex conjugate eigenvalues. The modes for the T field are simply extracted from ψi ̃ following the definition of the state vector. k (Fig. 8a) and erel , max (Fig. 8b) curves are reported for some selected cases in which r = 50, Δts = Δt and m ̃ assumes the values listed in Table 2. It is evident that the first 750 snapshots are barely enough to propagate the system state without introducing considerable error. Error curves show an increasing trend in the last part of the transient (right after the last used snapshot), meaning that some unstable eigenvalues appear among the dynamic spectrum identified by the λi . This is more evident for m ̃ = 500, as the approximation errors start soon to diverge. As a further proof of how the sampling strategy affects the dynamic properties of the DMD model, Fig. 9 shows the eigenvalues λi of the system for r = 25 and m ̃ = 1000. While no unstable eigenvalues are present (all eigenvalues lie within the unit circle, i.e. the stability region), it is evident that larger sampling time steps lead to oscillating modes. Given a pair of complex eigenvalues, indeed, the more is the angular distance from the real positive semi-axis, the more oscillating is

the associated dynamic mode. When Δts increases, the eigenvalues λi tend to move towards larger angles. For completeness, the real parts of the first 4 pairs of spatial modes ψi ̃ are also reported. They are shown in pairs since complex conjugate eigenvalue pairs correspond to complex conjugate pairs of modes. To allow visual representation, the modes are split back in separate fields corresponding to q‴ (Fig. 10), c1 (Fig. 11) and T (Fig. 12). It is stressed, however, that this holds only for visualisation purposes, as different reduced model fields can not evolve separately due to the coupling between state variables enforced by the definition of the state vector. 5. Conclusions In this paper, the Dynamic Mode Decomposition method was applied to the study of a fast transient of the Molten Salt Fast Reactor 10

Nuclear Engineering and Design 362 (2020) 110529

A. Di Ronco, et al.

following a variation of the heat transfer coefficient. A multiphysics solver developed with the OpenFOAM open-source CFD toolkit was used to solve numerically the model equations on a 2D axisymmetric computational domain and hence produce the full order high-fidelity data set. DMD builds reduced-size surrogate models of the original full order system without the need of projecting its equations on a reduced subspace, but starting only from available data (either numerical snapshots, as in this case, or possibly experimental data). In addition, relevant dynamics structures and stability information can be extracted from a lower-rank approximation of the linear operator of the system. The reduced order models were built using equally spaced time snapshots from a 25 s free-dynamics transient. A total of 1000 snapshots of the full state of the system were available. The volumetric power, the concentration of the first group of precursors and the temperature were selected to build the state vector, with a total number of degrees of freedom of about 7·104. Different sampling strategies and model sizes were adopted to test the capabilities of the method in capturing the dominant dynamics of the MSFR. It was observed that DMD is able to fully characterise the freeevolution dynamics of the reactor, provided that a suitable sampling frequency is adopted. The overall dynamic patterns, such as the final steady-state and the initial fast dynamics are fully caught with a few tenths of modes (i.e. degrees of freedom). The quality of the approximation is degraded only when the too few modes are retained, or when the sampling time step is relaxed too much. This is particularly true for the initial part of the transient, where the non-linearities are stronger and the time-scales are shorter. Results also showed that to correctly retrieve all stable modes, the sampling window must be chosen carefully. In the cases considered, snapshots needed to be sampled along the whole transient to eliminate spurious unstable modes. This aspect needs further investigation, for example, to assess how DMD surrogate models built at runtime can be used to propagate the state of the system with a significant reduction of the computational burden. The selection of the field quantities of interest was made on the basis of a qualitative judgement. Further analysis is needed to assess how to efficiently capture the coupling between different fields and to eliminate possible redundant information from the snapshots. This is important when very large models are considered due to memory considerations. Furthermore, the analysis of forced-dynamics transients is a challenging task. Since more recent developments of the DMD technique (Proctor et al., 2016) allow the analysis of controlled systems, possible application to large nuclear systems could be foreseen.

analysis. Chem. Eng. Sci. 111, 78–90. Bao, H., Dinh, N.T., Lane, J.W., Youngblood, R.W., 2019. A data-driven framework for error estimation and mesh-model optimization in system-level thermal-hydraulic simulation. Nucl. Eng. Des. 349, 27–45. Bennett, A., Martin, N., Avramova, M., 2019. A surrogate model based on sparse grid interpolation for boiling water reactor subchannel void distribution. Ann. Nucl. Energy 131, 51–62. Cervi, E., Lorenzi, S., Cammi, A., Luzzi, L., 2019. Development of a multiphysics model for the study of fuel compressibility effects in the Molten Salt Fast Reactor. Chem. Eng. Sci. 193, 379–393. Cervi, E., Lorenzi, S., Cammi, A., Luzzi, L., 2019. Development of an SP3 neutron transport solver for the analysis of the Molten Salt Fast Reactor. Nucl. Eng. Des. 346, 209–219. Cervi, E., Lorenzi, S., Luzzi, L., Cammi, A., 2019. Multiphysics analysis of the MSFR helium bubbling system: a comparison between neutron diffusion, SP3 neutron transport and Monte Carlo approaches. Ann. Nucl. Energy 132, 227–235. Chang, C.-W., Dinh, N.T., 2019. Classification of machine learning frameworks for datadriven thermal fluid models. Int. J. Therm. Sci. 135, 559–579. Diniz, P.S.R., Silva, E.A.B.D., Netto, S.L., 2002. Digital Signal Processing: System Analysis and Design. Cambridge University Press. Fiorina, C., Lathouwers, D., Aufiero, M., Cammi, A., Guerrieri, C., Kloosterman, J.L., Luzzi, L., Ricotti, M., 2014. Modelling and analysis of the MSFR transient behaviour. Ann. Nucl. Energy 64, 485–498. Gerardin, D., Allibert, M., Heuer, D., Laureau, A., Merle-Lucotte, E., Seuvre, C., 2017. Design evolutions of Molten Salt Fast Reactor. In: International Conference on Fast Reactors and Related Fuel Cycles: Next Generation Nuclear Systems for Sustainable Development (FR17), Yekaterinburg, Russia, June 26–29, 2017. German, P., Ragusa, J., Fiorina, C., Tano Retamales, M., 2019. Reduced-order modeling of parametrized multi-physics computations for the molten salt fast reactor. In: Proceedings of the International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, M&C 2019, pp. 1808–1817. Hardy, Z.K., Morel, J.E., Ahrens, C., 2019. Dynamic mode decomposition for subcritical metal systems. Nucl. Sci. Eng. 193 (11), 1173–1185. Hauer, J.F., Trudnowski, D.J., DeSteese, J.G., 2007. A perspective on WAMS analysis tools for tracking of oscillatory dynamics. In: Proceedings of the IEEE/EPS General Meeting, Tampa, FL, USA. Holmes, P., Lumley, J.L., Berkooz, G., Rowley, C., 2012. Turbulence, Coherent Structures, Dynamical Systems and Symmetry. University Press. Hu, Z., Smith, R.C., Willert, J., Kelley, C.T., 2014. High-dimensional model representations for the neutron transport equation. Nucl. Sci. Eng. 177 (3), 350–360. Ikonomopoulos, A., Alamaniotis, M., Chatzidakis, S., Tsoukalas, L.H., 2013. Gaussian processes for state identification in pressurized water reactors. Nucl. Technol. 182 (1), 1–12. Jasak, H., Jemcov, A., Tukovic, Z., 2007. OpenFOAM: a C++ library for complex physics simulations. In: International Workshop on Coupled Methods in Numerical Dynamics. IUC, Dubrovnik, Croatia. Lorenzi, S., Cammi, A., Luzzi, L., Rozza, G., 2017. A reduced order model for investigating the dynamics of the Gen-IV LFR coolant pool. Appl. Math. Model. 46, 263–284. Marrel, A., Marie, N., Lozzo, M.D., 2015. Advanced surrogate model and sensitivity analysis methods for sodium fast reactor accident assessment. Reliab. Eng. Syst. Saf. 138, 232–241. McClarren, R.G., 2019. Calculating time eigenvalues of the neutron transport equation with dynamic mode decomposition. Nucl. Sci. Eng. 193 (8), 854–867. Perkó, Z., Lathouwers, D., Kloosterman, J.L., van der Hagen, T., 2014. Large scale applicability of a Fully Adaptive Non-Intrusive Spectral Projection technique: sensitivity and uncertainty analysis of a transient. Ann. Nucl. Energy 71, 272–292. Proctor, J.L., Brunton, S.L., Kutz, J.N., 2016. Dynamic mode decomposition with control. SIAM J. Appl. Dyn. Syst. 15, 142–161. Sartori, A., Cammi, A., Luzzi, L., Rozza, G., 2016. A reduced basis approach for modeling the movement of nuclear reactor control rods. J. Nucl. Eng. Radiat. Sci. 2 (2). Schmid, P.J., 2010. Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–28. Tu, J.H., Rowley, C.W., Luchtenburg, D.M., Brunton, S.L., Kutz, J.N., 2014. On dynamic mode decomposition: theory and applications. J. Comput. Dyn. 1, 391–421. Weller, H.G., Tabor, G., Jasak, H., Fureby, C., 1998. A tensorial approach to computational continuum mechanics using object-oriented techniques. Computat. Phys. 12, 620–631. Wu, X., Kozlowski, T., 2017. Inverse uncertainty quantification of reactor simulations under the Bayesian framework using surrogate models constructed by polynomial chaos expansion. Nucl. Eng. Des. 313, 29–52. Wu, X., Mui, T., Hu, G., Meidani, H., Kozlowski, T., 2017. Inverse uncertainty quantification of TRACE physical model parameters using sparse gird stochastic collocation surrogate model. Nucl. Eng. Des. 319, 185–200. Wu, X., Kozlowski, T., Meidani, H., Shirvan, K., 2018. Inverse uncertainty quantification using the modular Bayesian approach based on Gaussian process, Part 1: Theory. Nucl. Eng. Des. 335, 339–355. Zhu, J., Shieh, L.S., Yates, R.E., 1985. Fitting continuous-time and discrete-time models using discrete-time data and their applications. Appl. Math. Model. 9, 93–98.

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. References Abdo, M., Elzohery, R., Roberts, J.A., 2019. Modeling isotopic evolution with surrogates based on dynamic mode decomposition. Ann. Nucl. Energy 129, 280–288. Alfonsi, A., Hummel, A., Sen, S., Strydom, G., Gougar, H., 2018. Decay Heat curve generation for high temperature reactors using exponentials, support vector machines and dynamic mode decomposition within the RAVEN framework. Trans. Am. Nucl. Soc. 118. Antoulas, A.C., 2005. An overview of approximation methods for large-scale dynamical system. Annu. Rev. Control 29, 181–190. Aufiero, M., Cammi, A., Geoffroy, O., Losa, M., Luzzi, L., Ricotti, M.E., Rouch, H., 2014. Development of an OpenFOAM model for the Molten Salt Fast Reactor transient

11