Application of multiphysics model order reduction to doppler/neutronic feedback

Chia sẻ: Huỳnh Lê Ngọc Thy | Ngày: | Loại File: PDF | Số trang:9

Thêm vào BST

Báo xấu

5
lượt xem 0
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

In this paper, a proper orthogonal decomposition based reduced-order model is presented for parametrized multiphysics computations. Our application physics is Doppler feedback in a simplified model of the molten salt fast reactor concept.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Application of multiphysics model order reduction to doppler/neutronic feedback

EPJ Nuclear Sci. Technol. 5, 17 (2019) Nuclear Sciences c P. German et al., published by EDP Sciences, 2019 & Technologies https://doi.org/10.1051/epjn/2019034 Available online at: https://www.epj-n.org REGULAR ARTICLE Application of multiphysics model order reduction to doppler/neutronic feedback Peter German 1 , Jean C. Ragusa 1, * , and Carlo Fiorina 2 1 Texas A&M University, Department of Nuclear Engineering, College Station, TX 77840, USA 2 Ecole Polytechnique F´ed´erale de Lausanne, Laboratory of Reactor Physics and Systems Behaviour, PH D3 465 (Batiment PH), Station 3, 1015 Lausanne, Switzerland Received: 19 May 2019 / Received in final form: 26 September 2019 / Accepted: 27 September 2019 Abstract. In this paper, a proper orthogonal decomposition based reduced-order model is presented for parametrized multiphysics computations. Our application physics is Doppler feedback in a simplified model of the molten salt fast reactor concept. The reduced model is created using the method of snapshots where the offline training set is obtained by exercising a full-order model created with the OpenFOAM based multiphysics solver, GeN-Foam. The steady state models solve the multi-group diffusion k-eigenvalue equations with moving precursors together with the energy equation. A fixed velocity field is assumed throughout the computations, hence the momentum and continuity equations are not solved. The discrete empirical interpolation method is used for the efficient coupling of the ROM solvers, while the input parameter space is surveyed using the improved distributed latin hypercube sampling algorithm. 1 Introduction reactivity feedbacks; they operate at atmospheric pres- sure; they allow for an online removal of gaseous fission Molten salt reactor (MSR) designs were originally devel- products; and they give the possibility to drain the fuel oped in the mid-1950s at Oak Ridge National Laboratory salt in passively cooled and critically-safe tanks in case (ORNL, USA) [1–3]. In MSRs, the nuclear fuel is in liq- of emergency. Most fast-spectrum MSRs currently under uid form, dissolved in a salt. Salt compositions vary, but development are based on pumped-loop designs, where are typically based on fluorides or chlorides and include the fuel salt is pumped outside of the primary vessel one or more of the following compounds: LiF, NaF, BeF2 , and transfers heat to a secondary coolant in separate ZrF4 , KF, NaCl, MgCl2 . Diverse variations on that reac- heat exchangers. Examples of fast-spectrum molten salt tor concept were investigated in the 1960s and 1970s, reactor designs include: the MOlten Salt Actinide Recy- including graphite-moderated thermal-spectrum reactors cler and Transforming (MOSART) project [9], the Molten at ORNL [4,5] as well as fast-spectrum burner reac- Salt Fast Reactor (MSFR) concept based on fluoride tors at Argonne National Laboratory [6,7]. In the early salt, developed in the EVOL (Evaluation and Viability 1970s, MSR research was in competition for U.S. fed- of Liquid Fuel Fast Reactors) [10–12] and then SAMO- eral funding with sodium-cooled fast reactor systems and FAR (Safety Assessment of the Molten Salt Fast Reactor) the MSR research in the USA dwindled down to low- [13] programs under the auspices of EURATOM, and the priority, low-funding efforts over the following decades, Molten Chloride Fast Reactor (MCFR), currently devel- while the thermal-spectrum light water-cooled pressur- oped by Terrapower [14]. Loop-type fast-spectrum molten ized and boiling water reactors and the sodium-cooled fast salt reactors present new modeling challenges: reactors became the reference baseline worldwide for ther- mal and fast spectrum systems, respectively. Nonetheless, – the fuel is in liquid form, yielding a more complex electricity-production priorities and safety requirements in level of multiphysics coupling than traditional light the nuclear sector have evolved over the last 50 years and water reactors (e.g., velocity fields needed to assess MSRs are now one of the six concepts selected for further space/time location of fuel and delayed neutron investigation in the frame of the Generation-4 Interna- precursors); tional Forum [8]. MSRs have highly promising features in – turbulent fuel-salt flow, leading to a large impact terms of sustainability and safety features. Indeed liquid- of turbulence modeling (thermal flow mixing in the fueled MSRs can be designed to have strong negative-only core, effects of nozzle inlets,... [15]); – presence of gas bubbles in the salt, leading to * e-mail: jean.ragusa@tamu.edu compressibility and reactivity effects; This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
2 P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) – high solidification temperature of the salt, requir- connected to reactor physics has already been developed ing the modeling of solidification/melting zones (wall in [29–33]. solidification, frozen drain functionality, see [15]); However, the utilization of POD-based ROM for multi- – high operating temperature in the salt, necessitating physics applications can be challenging because the oper- the study of thermal radiation heat transfer. ators in the full-order models often do not have an affine decomposition, meaning that the ROMs involve costly, full-order operations, possibly yielding negligible savings The modeling challenges in fast-spectrum MSRs largely in computation time. A method for the reduced-order prohibit the use of simulation packages developed and modeling of multiphysics problems in nuclear engineer- tailored specifically for Light Water Reactors (LWRs). ing has also been derived in [30]. The approach described For instance, the Virtual Environment for Reactor Appli- in the present work serves as an alternative. Instead of cations (VERA), developed as part of the Consortium assuming an overall linear temperature dependence for for Advanced Simulation of Light water reactors (CASL) the cross sections, here we opt for a hyper-reduction is being leveraged for MSR modeling but is only at an technique, namely, the Discrete Empirical Interpolation early stage of development for MSR and mostly focused Method (DEIM) [34]. This allows the handling of an on thermal-spectrum reactors, where salt flows mostly arbitrary, non-linear temperature dependence of the cross uni-directionally in graphite channels [16,17]. Hence, high- sections (in our case, in the logarithm of temperature). fidelity models, based on first-principle physics, become Our approach is exemplified using a liquid nuclear fuel the only available tool to gain understanding of the sig- system, more specifically, for a simplified 2D model of nificance of foreseeable phenomena unique to molten salt the Molten Salt Fast Reactor. The methodology has been and circulating fuel systems. Driven by the need for integrated into GeN-Foam. For different applications of high-fidelity computational fluid dynamics (CFD), many DEIM in other fields of engineering, we refer the reader research teams developing fast-spectrum MSRs rely on the to [35–37]. open-source OpenFOAM CFD platform, either directly The rest of the paper is organized as follows. In [18,19] or embedded in reactor physics packages such as Section 2, projection-based model order reduction is GeN-Foam (Generalized Nuclear Foam) [20,21]. reviewed, first for linear systems, then for nonlinear However, the simulation of such complex systems systems, such as the ones encountered in multiphysics requires the solution of coupled partial differential equa- applications. In Section 3, a full-order simulation model tions, which can be computationally expensive to obtain. is provided for a fast-spectrum molten salt reactor and Model-order reduction comprises a set of empirical and its associated reduced-order model is derived. Section 4 mathematical techniques that can be used for lowering describes the chosen MSFR computational model and its the computational complexity of the Full-Order Models parametric input space. Results comparing the FOM and (FOM) by creating Reduced-Order Models (ROM) that ROM models are provided in Section 5. We conclude and sacrifice a modest amount of accuracy for large gains in propose future work in Section 6. compute time. An effective mathematical technique for model-order reduction is the Reduced Basis (RB) method, which is based on the assumption that the solution of 2 Model order reduction: background a complex model lives in a relatively small subspace. When parametric studies are to be performed, one should In this section, we review some of the basics of projection- expand that subspace to account for solution variability. based model-order reduction. For brevity of the exposi- One manner by which such a subspace is “discovered” is tion, the reduced basis will always be sought through a through multiple full-order solutions, with adequate sam- POD decomposition and the reduced system will always pling of the parameter (input) space. This is known as be obtained using Galerkin projection. We refer the reader the method of snapshots. The information contained in to [38,39] for other reduced basis approaches and to [40] these snapshots is then extracted by means of Proper for Petrov-Galerkin projection-based ROM. Orthogonal Decomposition (POD) [22,23] (e.g., via cor- relation matrix or singular value decompositions). The 2.1 Model order reduction for linear systems FOM is then Galerkin-projected on the obtained basis functions to yield the ROM. In this work, a POD-based First, we consider a linear parameterized steady-state ROM is created for parametrized multiphysics computa- FOM. After discretization, the FOM can be written as tions on the Doppler feedback effect in the Molten Salt Fast Reactor (MSFR) [24,25]. In this technique, the FOM A(µ)x(µ) = b(µ), (1) is projected onto an appropriate subspace obtained via a POD of solutions determined by exercising the FOM where A(µ) ∈ RN ×N is the discretized linear operator itself. In other words, snapshots are taken at different (a possibly large system of dimension N ), x(µ) ∈ RN is states of the FOM and a POD is used to build a suit- the solution vector and the parametric dependence of the able basis for projection-based reduction. In the nuclear model is denoted by the d-dimensional input parameter engineering community, applications of POD-based ROMs vector µ = [µ1 , . . . , µd ]T . In order to discover the lower can be noted in reactor kinetics problems [26,27], fixed dimension manifold in which the parametric solutions can source, steady state neutral particle transport applications be adequately represented, the FOM is exercised for a cer- [28]. An approach for POD-ROM of eigenvalue problems tain number of realizations of the input parameter vector
P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) 3 µi , with 1 ≤ i ≤ NS and NS denoting the number of snap- reduced order model, compared to the original full-order shots. Letting S = [x(µ1 ), . . . , x(µNS )] ∈ RN ×NS denote model. To remedy this, the Discrete Empirical Interpola- the snapshot matrix, we perform a POD on S, which tion Method (DEIM) [34] is used in order to approximate constructs an orthonormal basis V = [v 1 , . . . , v r ] ∈ RN ×r F in a low-dimensional space by sampling it at only m where only the r most dominant modes are retained. A N components. Hence, during the offline stage, a second computationally effective manner to do this consists in matrix of snapshots is collected for the nonlinear function generating the correlation matrix Q of the snapshots as values, SF = [F (x(µ1 )), . . . , F (x(µNS ))] ∈ RN ×NS . Sim- ilarly to the method described in equations (2)–(4), the Q = S T S. (2) POD of SF is computed and m of the basis functions for the nonlinear vector-valued function are retained to subse- Then, the eigenvalue decomposition of the correlation quently interpolate F : [u1 , . . . , um ] = U ∈ RN ×m . DEIM matrix is obtained also selects m distinct interpolation points p1 , . . . , pm ∈ Q = WΛW T , (3) [1, N ] in order to assemble the DEIM interpolation point matrix P = [ep1 , . . . , epm ] ∈ RN ×m , where ei is the ith where W = [w1 , . . . , wNs ] is the matrix of eigenvectors canonical unit vector. Finally, the DEIM interpolant and Λ is a diagonal matrix containing the Λi eigenval- of F is ues. The orthonormal basis vectors in V can then be constructed using the snapshots by U (P T U )−1 P T F (x(µ)), Ns 1 X and the resulting nonlinear reduced-order model system vj = √ x(µi )wj,i . (4) Λi i=1 is Once this offline stage has been completed, the reduced Ar xr (µ) + V T U (P T U )−1 P T F (V xr (µ)) = br . (9) order model is obtained by projection of the full order system: Several remarks are in order: Ar (µ)xr (µ) = br (µ), (5) – V T U (P T U )−1 ∈ Rr×m is a small matrix that can with the reduced linear operator Ar (µ) = V T A(µ)V ∈ be pre-computed once for all; – P T F extracts only a few (actually m, with m Rr×r , reduced right-hand side br (µ) = V T b(µ) ∈ Rr , and N ) components of F that need to be evaluated. the reduced state vector xr (µ) ∈ Rn . This linear system This is a large savings afforded by the use of the of size r N , is solved for xr , which are the coefficients of DEIM. Because discretization of partial differen- the full-order solution in the basis V . Hence, the full-order tial equations usually results in local connectivity solution is reconstructed as between an unknown and its neighbors, evalu- r X ating P T F (V x(µ)) is therefore computationally x≈ v i xr,i = V xr . (6) inexpensive when m N . i=1 For additional information on DEIM, we refer the reader For additional information, we refer the reader to [41,42]. to [34,43]. 2.2 Model order reduction for nonlinear systems Next, we consider a nonlinear full-order model (FOM) 3 Governing laws which, after discretization, can be written as In this section, we select mathematical models (partial dif- A(µ)x(µ) + F (x(µ), µ) = b, (7) ferential equations) that will be used as FOM and derive their reduced-order counterparts. These models are said where a nonlinear vector-valued function has been added, to be parametric in the sense that certain parameters are F (x(µ), µ) ∈ RN . Henceforth, for the sake of easier read- uncertain in the governing equations and will be sam- ability, the dependence of the operators on µ is only pled from appropriate probability density functions. In assumed and not shown explicitly. If one carries out this work, the following parameter are assumed uncer- the procedures described in Section 2.1, the resulting tain: the total reactor power (Pth ), the heat exchange nonlinear reduced-order model is coefficient in the heat sink, i.e., in the heat exchanger Ar xr (µ) + V T F (V xr (µ)) = br . (8) (αext ), the ultimate heat sink temperature (Text ), and the volumetric area of the heat sink (AV ). The vector con- Even though this equation is expressed in terms of the vec- taining all the uncertain parameters will be denoted by tor of reduced unknowns, xr , solving it requires evaluating µ = (αext , AV , Text , Pth )T . the full-order nonlinear vector-valued function, of size The mathematical models are selected to be representa- N r. This can be computational very expensive, often tive of a simplified molten salt reactor design. The molten resulting in limited CPU time savings when solving the salt reactor application is discussed in a later section.
4 P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) 3.1 Full-Order Model (FOM) of the uncertain parameters (the RANS continuity and momentum equations were solved using k− model for In order to represent the physics of a MSFR, we define a turbulence [21]). Thus, ~u and αeff are known fields dur- multiphysics model comprising of ing the generation of the reduced operators. We leave – neutronics: six-group k-eigenvalue diffusion theory, model-order reduction of turbulent hydrodynamics for with delayed neutron precursor balance equation subsequent work. The multigroup neutron fluxes are nor- with a drift (advection) term, and malized using the group-wise power cross section Σp,g to – thermal-hydraulics: here, we only consider an energy ensure a certain total reactor power, Pth : balance equation, with a given (but spatial vary- G ing) flow field computed using nominal values of the Z X 3 uncertain parameters. For the types of parametric d r Σp,g (r)Φg (r) = Pth . (12) domain g=1 modifications employed here, the effect of flow per- turbations should be small. For the reduced-order modeling for turbulent flows, we refer the reader to The cross sections are assumed to be temperature- [44,45]. dependent. For fast spectrum reactors, as shown in [47], a logarithmic interpolation between pre-computed data- The cross sections in the neutronics models are bases (for example at 900 K and 1200 K) of the group temperature-dependent. The multi-group diffusion k- constants yields good results: eigenvalue problem [46] can be described as a coupled partial differential equation system expressing the balance Σ1200 − Σ900 T of neutrons in each energy bin as Σ(T ) = Σ900 + log log 1200 900 900 G ~ g ) + Σr,g Φg = (1 − β)χp,g = Σconst + Σvar log (T ). (13) X ~ · (Dg ∇Φ −∇ νΣf,g0 Φg0 k 0 g =1 I G where T denotes the temperature field. X X Finally, to be able to account for the temperature +χd,g λz Cz + Σg0 →g Φg0 , (10) feedback, the following energy balance equation is solved: i=z g 0 6=g G where G denotes the number of energy groups, I the X ~ · (~uρcp T ) = ∇ ∇ ~ · (kT ∇T ~ ) − αAV (T − Text ) + Σp,g Φg , number of delayed neutron precursor groups, Φg is the neutron scalar flux, Dg is the diffusion coefficient, Σr,g g=1 is the macroscopic removal cross-section and νΣf,g is (14) the total fission neutron yield times the macroscopic fis- where ρ is the density, cp is the heat capacity, kT the sion cross-section in energy group g. Furthermore, Σg0 →g effective thermal conductivity of the fluid, while α is the is the macroscopic scattering cross-section from group heat transfer coefficient and AV is the volumetric area g 0 to g, χp,g and χd,g describe the fraction of prompt of the heat sink. All of the parameters in the equation and delayed neutrons released in group g, while β is above are assumed to be constant. Furthermore, the heat the effective delayed neutron yield. Moreover, λz denotes transfer coefficient α is computed as the harmonic mean the decay constant corresponding to the delayed neutron of two coefficients, one that characterizes the heat transfer precursor concentration Cz of group precursor group z. between the salt and the structure of the heat exchanger The problem is supplemented with boundary conditions: (αsalt ) and a second describing the heat flow between the reflective boundary conditions (∇Φ ~ g · ~n = 0, g = 1, . . . , G) heat exchanger and the external heat sink (αext ). For sim- are applied on the symmetry planes and zero-incoming plicity, zero gradient boundary conditions (∇T~ · ~n = 0) are current boundaries are applied on the other boundaries used for every surface in the model. The iteration scheme (−Dg ∇Φ~ g · ~n = 1 Φg , g = 1, . . . , G). These equations are used for the solution of the coupled problem is discussed in 2 coupled with the steady-state balance equations for the details in [21]. αext , AV and Text are assumed uncertain; delayed neutron precursors, commonly expressed as hence this model is parametric in those input parameters. G ~ z ) = βz X ~ · (~uCz ) − ∇ ∇ ~ · (αeff ∇C νΣf,g Φg − λz Cz , (11) 3.2 Reduced-Order Model (ROM) k g=1 The ROMs are constructed by physics-wise (equation- where ~u is a precomputed (stationary) velocity field, βz wise) projection of the FOM equations onto suitable is the delayed neutron yield for precursor group z and reduced bases obtained applying POD to the solution αeff is an effective diffusion coefficient that takes into snapshots. During the offline phase, the solution fields account the effects of turbulent mixing as well. For sim- are collected into corresponding snapshot matrices and a plicity, zero gradient boundary conditions (∇C ~ z · ~n = 0 Proper Orthogonal Decomposition is carried out for each for z = 1, . . . , I) are used for each of the precursor equa- snapshot matrix separately. This segregated approach has tions on every wall. We stress that ~u is not uniform, but been proven to be effective for k-eigenvalue multigroup taken from a CFD simulation using the nominal values problems [33]. The approximate solutions in the reduced
P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) 5 spaces can be written as as follows: rg X rz X l = (P T Lr )−1 log(P T T r t). Φg ≈ ψg,i fg,i = Ψr,g fg , Cz ≈ Cz,i cz,i = Cr,d cd , i=1 i=1 It should be noted that, in this case, it is enough the carry (15) out the DEIM in terms of the logarithm of the temper- rT X rl X ature due to the fact that the nonlinearity in the model T ≈ τi ti = Tr t, log (T ) ≈ Li li = Lr l is caused by the temperature-dependence of the cross sec- i=1 i=1 tions, and the non-linear function involves the logarithm (16) of the temperature. As noted in the previous Section, we where ψg,i is the ith basis vector of the subspace selected recall that the reduced operator can be precomputed in a for the neutron flux in group g, Cz,i is the ith basis tensor form and every time the operator has to be recon- vector for the precursor group z, τi denotes the basis vec- structed due to the changing temperature field, it only tor of temperature and Li is the ith basis vector of the requires the summation of reduced matrices making the logarithmic temperature. Moreover, fg , cd , t and l vec- coupling of the reduced models extremely efficient. The tors contain the coordinates of the approximated flux in same treatment is applied to the other reduced operators group g, approximated precursor concentration in group d, containing cross-sections, even though in the following def- approximated temperature and logarithmic temperature initions it is not shown explicitly. The additional reduced within their corresponding reduced spaces. Using these, operators in equation (17) are computed as the ROMs for the multi-group diffusion equations can be described as (S r,g )i,j = hψg,i , Σr,g ψg,j i, (S f,g0 )i,j = hψg,i , (1 − β)χp,g νΣf,g0 ψg0 ,j i, (20) G I 1 X X (Lz )i,j = hψg,i , χd,g λz Cz,j i, Rg f g + S r,g f g = S f,g0 f g0 + Lz cz kr 0 z=1 (S s,g0 )i,j = hψg,i , Σg0 →g ψg0 ,j i. (21) g =1 G X One can notice that these reduced matrices may be rect- + S s,g0 f g0 , (17) angular depending on the number of POD modes used for g 0 6=g the different reduced bases. Similarly, the reduced form of the precursor equations can be expressed as where kr is the largest eigenvalue of the reduced sys- tem and the reduced operators are computed using the G 1 X approximation mentioned in equation (15) together with F z cz − M z cz = E z,g f g − L∗z cz , (22) a Galerkin projection onto the corresponding subspace. kr g=1 Thus, the elements of the reduced diffusion operator can be expressed as where the entries of the reduced operators are computed as ~ · (Dg ∇ψ (Rg )i,j = hψg,i , −∇ ~ g,j )i, (18) ~ · (~uCz,j )i, (F z )i,j = hCz,i , ∇ where h · i denotes the volumetric integral of the given (M z )i,j = hCz,i , ∇~ · (αeff ∇C ~ z,i )i, (23) scalar fields that can be carried out numerically. Even though it is not shown explicitly, this reduced operator (E z,g )i,j = hCz,i , βz νΣf,g ψg,j i, takes into account the boundary conditions imposed on (L∗ z )i,j = hCz,i , λz Cz,j i. (24) the scalar flux by incorporating a hψg,i , 12 ψg,j iΓ term (in case of vacuum boundary) for every boundary face Γ of As a last step, the reduction of the energy equation is the computational domain. Again, it must be mentioned, carried out as that Dg depends on the approximate logarithmic temper- G ature (see Eq. (16)). By translating this linear dependence X Ht = Kt − At − a + S p,g f g , (25) into the ROM, the expressions of the reduced operators g=1 become slightly more convoluted and can be written as where the entries of the reduced matrices and sink vector ~ · (Dgconst ∇ψg,j )i (Rg )i,j = hψg,i , −∇ are given by rl X ~ · (Dgvar Lk ∇ψg,j )ilk , (19) ~ · (~uρcp τj )i (K)i,j = hτi , ∇ (H)i,j = hτi , ∇ ~ · (kT ∇τj )i + hψg,i , −∇ k=1 (26) (A)i,j = hτi , −αAV τj i (a)i = hτi , αAV Text i (27) where the values of Dgconst and Dgvar can be determined using equation (13) and the coefficients of the logarith- (S p,g )i,j = hτi , Σp,g ψg,j i. (28) mic temperature (li ) are computed using the coefficients of the temperature (ti ) through the Discrete Empirical Finally, we note that the elements of input parameter Interpolation Method, described in Section 2.2 in detail, vector µ = (αext , AV , Text , Pth )T are simply factors in the
6 P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) Table 1. The energy group structure used for the computations [52]. Energy group Lower bound (MeV) Upper bound (MeV) 1 2.231E−00 2 4.979E−01 2.231E−00 3 2.479E−02 4.979E−01 4 5.531E−03 2.479E−02 5 2.485E−04 5.531E−03 6 2.485E−04 Fig. 1. The geometry and flow field used for the snapshot gener- Table 2. The dimensions of the sub-spaces used for the ation with the FOM. Dimensions are shown in mm. (P – Pump, approximation of the fields and the projection of the HX – Heat Exchanger). equations. Field dim. Field dim. Field dim. Field dim. reduced order equations as well, thus the reconstruction Φ1 2 Φ5 2 C3 1 C7 2 of the reduced matrices is not necessary for repeating sim- Φ2 2 Φ6 2 C4 1 C8 2 ulations with the ROM. For the solution of the coupled Φ3 2 C1 1 C5 1 T 3 problem, the standard GeN-Foam iteration scheme is used Φ4 2 C2 1 C6 1 log(T ) 2 which is discussed in paper [48] in detail. The same itera- tive scheme is used to solve both the FOM and the ROM models. To compare the solutions of the ROM with those Altogether 20 parameter vectors are sampled for of the FOM we use the following two error indicators: snapshot generation using the Improved Distributed Latin Hypercube Sampling (IHS) method [53]. The ∆k = |k − kr |, (29) uncertain parameters in these vectors are varied in a ±20% interval around their mean values µ ¯ = that describes the absolute difference between the largest (105 W m−2 K, 100 m2 m−3 , 900 K, 1440 MWth ). Using eigenvalue of FOM and the ROM, and these snapshots, altogether 16 reduced spaces are created ||ζFOM − ζROM ||L2 using POD. One for the neutron flux in each energy group, eζ = (30) one for the precursor concentration in each group, one for ||ζFOM ||L2 the temperature and one for the logarithmic temperature. The decay in the eigenvalue of the correlation matrices that gives the relative L2 error in the field variables, where (defined in Eq. (2)) are presented in Figure 2; the eigen- ζ can be Φg (g = 1, . . . , G), Cz (z = 1, . . . , I) or T . values are normalized to their largest one. It is visible that for all field variables, the eigenvalues decay rapidly, sug- 4 MSFR computational model gesting that only a few modes are enough to approximate the full-order solution. It can also be observed that, as A 2D axisymmetric model of the MSFR has been created the half-life of the precursor group decreases (from group using the available information in [18,48,49]. The dimen- 1 to 8), the decay in the eigenvalues of the corresponding sions of the geometry are provided in Figure 1 together correlation matrices is slower. with the pre-computed velocity field. The velocity field The dimensions of the reduced spaces for the differ- has been obtained by a standalone steady state solve of ent solution fields are summarized in Table 2. These the incompressible porous Navier-Stokes equations with numbers are acquired using a strategy based on the k- turbulence model [21]. The heat exchanger (HX) is energy-retention limit defined as modeled as a porous medium responsible for flow resis- ri P tance and heat sink, while the pump (P) is modeled as a λi simple volumetric momentum source. For more informa- i Ns < 1 − lim , tion about the semi-empirical expressions used to compute P λi the parameters of the flow resistance and heat sink, the i reader may refer to [21]. The mesh used for the computations has been created where λi are the eigenvalues of the correlation matrix built using SALOME [50] and contains 11,064 cells. For the using the snapshots of a selected field, r is the number neutronics computations, six energy and eight precur- of POD modes retained for this field and (1 − lim ) is sor groups are used, meaning that the full-order model the energy-retention limit. This characterizes the error has 165,960 degrees of freedom. The corresponding group in the reconstruction of the snapshots originating from constants are generated using Serpent 2 Monte Carlo discarding the rest of the POD modes (r + 1, . . . , Ns ). Transport code [51] for two salt temperatures, 900 K and When lim = 0, all of the POD modes are used, hence 1200 K. The bounds of the used energy group structure the snapshots can be reconstructed exactly. In this work, are presented in Table 1. lim = 10−7 was used.
P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) 7 Fig. 2. The decay in the eigenvalues of correlation matrices for the neutron flux (left), precursor concentration (middle) and temperature (right). Fig. 3. The reduced order solution for the scalar flux in energy Fig. 4. The reduced order solution for the precursor concentra- group one together with its absolute deviation from the full order tion in group eight together with its absolute deviation from the solution. full order solution. The relatively low dimensions can be explained by the fact that all the uncertain parameters appear in the thermal balance equation and changing them does not influence the distributions of flux and temperature fields considerably. This also means that the 165,960 degrees of freedom in the FOM can be reduced to 27 unknowns in the ROM. 5 Results After building the reduced-order model, a new real- ization of the input parameter vector, µ∗ = (9.3 × Fig. 5. The reduced order solution for the temperature together 104 W m−2 K, 111.1 m2 m−3 , 833.3 K, 1320 MWth ) is with its absolute deviation from the full order solution. selected and simulations are performed with both the ROM and the FOM in order to test the accuracy and efficacy of the ROM. We stress that µ∗ was not included in the training set. Figure 3 shows the reconstructed scalar factors coming from the ROM and FOM was 8.5 pcm, flux in energy group one from the ROM together with its which is also satisfactory. It is worth noting that loosen- absolute deviation from the full order solution. It is vis- ing the energy-retention limit to 10−6 resulted a 77.6 pcm ible that the maximum error is more than three orders difference in the eigenvalues, while tightening it to 10−8 of magnitude lower than the maximum value of the full gave 22.3 pcm due to the inclusion of a third POD mode order solution. for the logarithmic temperature that decreased the accu- Figure 4 compares the precursor concentration in pre- racy in the temperature coefficients. By further increasing cursor group eight. Again, it is visible that the maximum the energy-retention limit, we observed no further change absolute deviation is approximately three orders of mag- from the 22 pcm difference level. nitude lower than the maximum of the original solution. Furthermore, the relative L2 errors in the field variables, Figure 5 shows the reconstructed temperature profile of defined at equation (30), are summarized in Table 3. It the ROM with its relative deviation from the FOM. Again, can be observed that none of the fields have a relative L2 the maximum relative error is slightly above 1%. Further- error above 1%. The speed-up factor for the solution of more, the deviation between the effective multiplication the problem was 1550 for this specific example.
8 P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) Table 3. The dimensions of the sub-spaces used in the 9. V. Ignatiev, O. Feynberg, I. Gnidoi, A. Merzlyakov, projections. A. Surenkov, V. Uglov, A. Zagnitko, V. Subbotin, I. San- nikov, A. Toropov, et al., Molten salt actinide recycler and ζ eζ (%) ζ eζ (%) ζ eζ (%) ζ eζ (%) transforming system without and with Th-U support: fuel Φ1 0.034 Φ5 0.032 C3 0.006 C7 0.023 cycle flexibility and key material properties, Ann. Nucl. Φ2 0.033 Φ6 0.085 C4 0.003 C8 0.028 Energy 64, 408 (2014) Φ3 0.032 C1 0.007 C5 0.006 T 0.596 10. L. Mathieu, D. Heuer, E. Merle-Lucotte, R. Brissot, Φ4 0.031 C2 0.007 C6 0.013 C. Le Brun, E. Liatard, J.-M. Loiseaux, O. M´eplan, A. Nuttin, D. Lecarpentier, Possible configurations for the thorium molten salt reactor and advantages of the fast nonmoderated version, Nucl. Sci. Eng. 161, 78 (2009) 6 Conclusions 11. D. Heuer, E. Merle-Lucotte, M. Allibert, M. Brovchenko, V. Ghetta, P. Rubiolo, Towards the thorium fuel cycle with An efficient way of coupling parametrized neutronics molten salt fast reactors, Ann. Nucl. Energy 64, 421 (2014) and energy equations has been developed for POD-based 12. EVOL (Project no 249696) Final Reportr, https://cordis. ROMs. It utilizes the group-wise projection of the multi- europa.eu/docs/results/249/249696/final1-final-report-f. group diffusion equations, including the precursor balance pdf, Accessed: 04-10-2019 equations, and the enthalpy equation in moving fuel sys- 13. SAMOFAR (Safety Assessment of the Molten Salt Fast tems. In the current state of the development, the velocity Reactor), http://samofar.eu/, Accessed: 04-10-2019 equations are not solved; a constant velocity field is 14. MCFR Solutions: Nuclear Innovation for New Options in assumed. The non-linear temperature dependence of the American Industry, https://terrapower.com/technologies/ group constants is handled using DEIM. To demonstrate mcfr, Accessed: 04-10-2019 the viability of the method, it has been implemented in 15. M. Tano-Retamales, Mod´elisation multi-physique multi-´echelle de caloporteurs sels fondus et validation OpenFOAM based multiphysics solver, GeN-Foam. More- exp´erimentale, Ph.D. thesis, Universit´e Grenoble-Alpes, over, a ROM has been created for the MSFR and the 2018 corresponding reduced spaces have been generated by 16. B.S. Collins, C.A. Gentry, A.J. Wysocki, R.K. Salko, Molten exercising a simplified 2D axisymmetric FOM. The lim- salt reactor simulations using MPACT-CTF, Tech. rep., ited number of parameters and the underlying physics Oak Ridge National Lab.(ORNL), Oak Ridge, TN (United resulted in reduced spaces with dimensions between 1-3. States), 2017 A test has been performed with four uncertain parameters 17. C.A. Gentry, B.R. Betzler, B.S. Collins, Initial Benchmark- in the energy equation and the results of the ROM have ing of ChemTriton and MPACT MSR Modeling Capabili- been compared to those of the FOM. A good agreement ties, Tech. rep., Oak Ridge National Laboratory (ORNL), is observed both in terms of the solution vectors and the Oak Ridge, TN (United States), 2017 effective multiplication factors. The experienced speed up 18. M. Aufiero, A. Cammi, O. Geoffroy, M. Losa, L. Luzzi, M.E. factor was above 1,500. Ricotti, H. Rouch, Development of an OpenFOAM model for the Molten Salt Fast Reactor transient analysis, Chem. Eng. Sci. 111, 390 (2014) This material is based upon work supported under an NEUP- 19. A. Laureau, D. Heuer, E. Merle-Lucotte, P. Rubiolo, IRP Award of the U.S. Department of Energy, Office of Nuclear M. Allibert, M. Aufiero, Transient coupled calculations of Energy (contract reference DE-NE0008651). the Molten Salt Fast Reactor using the transient fission matrix approach, Nucl. Eng. Des. 316, 112 (2017) 20. J. Bao, Development of the model for the multiphysics anal- References ysis of Molten Salt Reactor Experiment using GeN-Foam code, Master’s thesis, EPFL & ETH Zurich, 2016, https: 1. M.W. Rosenthal, P.R. Kasten, R.B. Briggs, Molten-salt //www.psi.ch/fast/PublicationsEN/FB-DOC-16-014.pdf Reactors – history, status, and potential, Nucl. Appl. 21. C. Fiorina, I. Clifford, M. Aufiero, K. Mikityuk, GeN-Foam: Technol. 8, 107 (1970) a novel OpenFOAM R based multi-physics solver for 2D/3D 2. D. LeBlanc, Molten salt reactors: a new beginning for an transient analysis of nuclear reactors, Nucl. Eng. Des. 294, old idea, Nucl. Eng. Des. 240, 1644 (2010) 24 (2015) 3. H.G. MacPherson, The molten salt reactor adventure, Nucl. 22. K. Pearson, LIII. On lines and planes of closest fit to systems Sci. Eng. 90, 374 (1985) of points in space, Philos. Mag. 2, 559 (1901) 4. R. Robertson, MSRE design and operations report. Part 23. R. Pinnau, Model reduction via proper orthogonal decompo- I. Description of reactor design, Tech. Rep. ORNL-TM- sition, in Model Order Reduction: Theory, Research Aspects 728, 1965, DOI:10.2172/4654707 and Applications (Springer, 2008), p. 95 5. M.W. Rosentahl, The development status of molten-salt 24. E. Merle-Lucotte, D. Heuer, M. Allibert, M. Brovchenko, breeder reactors, Tech. Rep. ORNL-4812, 1972 N. Capellan, V. Ghetta, Launching the thorium fuel cycle 6. L.G. Alexander, Molten-Salt Fast Reactors, Tech. Rep. with the molten salt fast reactor, in Proceedings of ICAPP ANL-6792, 1963 (2011), p. 2 7. D. Holcomb, et al., Fast Spectrum Molten Salt Reactor 25. M. Allibert, M. Aufiero, M. Brovchenko, S. Delpech, Options, Tech. Rep. ORNL/TM-2011/105, 2011 V. Ghetta, D. Heuer, A. Laureau, E. Merle-Lucotte, Molten 8. The Generation IV International Forum (GIF), https:// salt fast reactors, in Handbook of Generation IV Nuclear www.gen-4.org/gif/, Accessed: 04-10-2019 Reactors (Elsevier, 2016), p. 157
P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) 9 26. A. Sartori, D. Baroli, A. Cammi, D. Chiesa, L. Luzzi, 39. Z. Bai, Krylov subspace techniques for reduced-order mod- R. Ponciroli, E. Previtali, M.E. Ricotti, G. Rozza, eling of large-scale dynamical systems, Appl. Numer. Math. M. Sisti, Comparison of a Modal Method and a Proper 43, 9 (2002) Orthogonal Decomposition approach for multi-group time- 40. S. Lorenzi, An adjoint proper orthogonal decomposition dependent reactor spatial kinetics, Ann. Nucl. Energy 71, method for a neutronics reduced order model, Ann. Nucl. 217 (2014) Energy 114, 245 (2018) 27. D.P. Prill, A.G. Class, Semi-automated proper orthogo- 41. M. Rathinam, L.R. Petzold, A new look at proper orthog- nal decomposition reduced order model non-linear anal- onal decomposition, SIAM J. Numer. Anal. 41, 1893 ysis for future BWR stability, Ann. Nucl. Energy 67, (2003) 70 (2014) 42. Y. Liang, H. Lee, S. Lim, W. Lin, K. Lee, C. Wu, 28. A.G. Buchan, A. Calloo, M.G. Goffin, S. Dargaville, F. Fang, Proper orthogonal decomposition and its applications–part C.C. Pain, I.M. Navon, A POD reduced order model for I: theory, J. Sound Vib. 252, 527 (2002) resolving angular direction in neutron/photon transport 43. B. Peherstorfer, D. Butnaru, K. Willcox, H. Bungartz, problems, J. Comput. Phys. 296, 138 (2015) Localized discrete empirical interpolation method, SIAM J. 29. A. Buchan, C. Pain, F. Fang, I. Navon, A POD reduced- Sci. Comput. 36, A168 (2014) order model for eigenvalue problems with application to 44. S. Lorenzi, A. Cammi, L. Luzzi, G. Rozza, POD-Galerkin reactor physics, Int. J. Numer. Methods Eng. 95, 1011 method for finite volume approximation of Navier–Stokes (2013) and RANS equations, Comput. Methods Appl. Mech. Eng. 30. A. Sartori, A. Cammi, L. Luzzi, G. Rozza, A multi-physics 311, 151 (2016) reduced order model for the analysis of Lead Fast Reactor 45. S. Hijazi, G. Stabile, A. Mola, G. Rozza, Data-driven single channel, Ann. Nucl. Energy 87, 198 (2016) POD-Galerkin reduced order model for turbulent flows, 31. C. Wang, H.S. Abdel-Khalik, Construction of accuracy- arXiv:1907.09909 preserving surrogate for the eigenvalue radiation diffusion 46. J.J. Duderstadt, L.J. Hamilton, in Nuclear reactor analysis and/or transport problem, Tech. rep., American Nuclear (Wiley-Interscience, Ann Arbor, Michigan, 1976), p. 650 Society, Inc., 555 N. Kensington Avenue, La Grange Park, 47. M. Aufiero, A. Cammi, C. Fiorina, L. Luzzi, A. Sartori, Illinois 60526, United States, 2012 A multi-physics time-dependent model for the Lead Fast 32. C. Wang, H.S. Abdel-Khalik, U. Mertyurek, Crane: A new Reactor single-channel analysis, Nucl. Eng. Des. 256, 14 scale super-sequence for neutron transport calculations, (2013) in Proceeding of MC 2015, Nashville, TN, April 19|23, 48. C. Fiorina, D. Lathouwers, M. Aufiero, A. Cammi, C. 2015 Guerrieri, J.L. Kloosterman, L. Luzzi, M.E. Ricotti, Mod- 33. P. German, J.C. Ragusa, Reduced-order modeling of param- elling and analysis of the MSFR transient behaviour, Ann. eterized multi-group diffusion k-eigenvalue problems, Ann. Nucl. Energy 64, 485 (2014) Nucl. Energy 134, 144 (2019) 49. E. Merle, Concept of Molten Salt Fast Reactor, in Molten 34. S. Chaturantabut, D.C. Sorensen, Nonlinear model reduc- Salt Reactor Workshow, PSI, 2017 tion via discrete empirical interpolation, SIAM J. Sci. 50. SALOME, http://www.salome-platform.org, Accessed: 04- Comput. 32, 2737 (2010) 10-2019 35. A. Hochman, B.N. Bond, J.K. White, A stabilized discrete 51. J. Lepp¨ anen, M. Pusa, T. Viitanen, V. Valtavirta, T. empirical interpolation method for model reduction of elec- Kaltiaisenaho, The serpent monte carlo code: Status, devel- trical, thermal, and microelectromechanical systems, in 2011 opment and applications in 2013, in SNA+ MC 2013-Joint 48th ACM/EDAC/IEEE Design Automation Conference International Conference on Supercomputing in Nuclear (DAC), IEEE (2011), p. 540 Applications+ Monte Carlo (2014), p. 6 36. A. Radermacher, S. Reese, POD-based model reduction with 52. C. Fiorina, The Molten Salt Fast Reactor as a Fast- empirical interpolation applied to nonlinear elasticity, Int. J. Spectrum Candidate for Thorium Implementation, Ph.D. Numer. Methods Eng. 107, 477 (2016) thesis, Politecnico di Milano, https://www.politesi.polimi. 37. W. Yao, S. Marques, Nonlinear aerodynamic and aeroelas- it/bitstream/10589/74324/1/2013 03 PhD Fiorina.pdf, tic model reduction using a discrete empirical interpolation 2013 method, AIAA J. 55, 624 (2017) 53. B. Beachkofski, R. Grandhi, Improved distributed hyper- 38. S. Gugercin, A.C. Antoulas, A survey of model reduction by cube sampling, in 43rd AIAA/ASME/ASCE/AHS/ASC balanced truncation and some new results, Int. J. Control Structures, Structural Dynamics, and Materials Conference 77, 748 (2004) (2002), p. 1274 Cite this article as: Peter German, Jean C. Ragusa and Carlo Fiorina, Multiphysics model order reduction and application to molten salt fast reactor simulation, EPJ Nuclear Sci. Technol. 5, 17 (2019)