CC BY
11Математические заметки СВФУ Апрель—июнь, 2021. Том 28, № 2
UDC 519.63
MULTISCALE MODEL REDUCTION FOR NEUTRON DIFFUSION EQUATION
A. O. Vasilev, D. A. Spiridonov, and A. V. Avvakumov
Abstract. Modelling of dynamic processes in nuclear reactors is carried out, mainly, on the basis of the multigroup diffusion approximation for the neutron flux. The neutron diffusion approximation is widely used for reactor analysis and applied in most engineering calculation codes.
In this paper, we attempt to employ a model reduction technique based on the multiscale method for neutron diffusion equation. The proposed method is based on the use of a generalized multiscale finite element method. The main idea is to create multiscale basis functions that can be used to effectively solve on a coarse grid. From calculation results, we obtain that multiscale basis functions can properly take into account the small-scale characteristics of the medium and provide accurate solutions. The results calculated with the GMsFEM are compared with the reference fine-grid calculation results.
DOI: 10.25587/SVFU.2021.27.52.008
Keywords: parabolic equation, neutron diffusion, multiscale simulation, generalized multiscale finite element method (GMsFEM).
1. Introduction
Among the basic physical processes occurring in a nuclear reactor [1], the focus is on neutron transport. To predict the neutron distribution, a rather complex integro-differential equation is involved; besides, the neutron flux distribution depends on time, energy, spatial and angular variables [2]. In practical calculations of nuclear reactors, as a rule, different approximative forms of the neutron transport equation are used. For example, the most neutronics codes are based on the multigroup diffusion equations. The computational model is a system of coupled second-order parabolic equations, which is supplemented by a system of ordinary differential equations to take into account delayed neutrons.
The stationary state of neutron flux, which is related to the critical state of the reactor, is characterised by local balancing of neutron absorption and generation. This boundary state is usually described by solution of a spectral problem (Lambda Modes problem [3-5]) provided that the fundamental eigenvalue (maximal eigenvalue) that is called k-effective of the reactor core, is equal to unity. In this case, the
This work was supported by the grant of the Russian Science Foundation (#19—71—00008).
© 2021 A. O. Vasilev, D. A. Spiridonov, and A. V. Avvakumov
stationary neutron field is related with the corresponding eigenfunction. Computations of k-effective of the reactor using the Lambda Modes problem are obligatory for reactor design calculations.
Real dynamic neutronic calculations require the use of large computational grids during large characteristic times. For these reasons, the reactor calculations can be very expensive. When solving these problems, it is common to use numerical averaging techniques or different multiscale methods that can significantly reduce the number of unknowns, through the construction of coarse-scale model. Also note a class of methods for modelling the non-stationary neutron diffusion, which is associated with the multiplicative representation of the space-time factorization methods and the quasistatic method [6, 7]. The general strategy for the approximate solution of the non-stationary problems of neutron transport, which is oriented to fast calculations using the state change modal method was formulated in [8].
Homogenization methods are used for coarse-grid approximation and allow calculating the effective properties of the domain. As a result, the methods of numerical homogenization establish macroscopic laws on the basis of local computations, but such approaches are usually based on a priori formulated assumptions, see e.g. [9-11]. On the other hand, multiscale methods construct multiscale basis functions in each coarse-grid block and couple these basis functions in a global formulation. These methods have a two-way information exchange between micro-and macro- scales. There are many variations of multiscale finite element methods (MsFEM). A detailed description of several variations of the method is presented in [12].
In this paper, we propose a Generalized Multiscale Finite Element Method (GMsFEM [13,14]) to solve the neutron diffusion equation. The main idea of the GMsFEM is the dimensionality reduction of the problem. GMsFEM involves two basic steps: the construction of multiscale basis functions that take into account small scale heterogeneities in local domains and the construction of the coarse scale approximation. For the construction of the multiscale basis functions we solve spectral problems in local domains. Spectral problems allow identifying the most important characteristics of the solution. In contrast to the available techniques, this method allows to avoid the limitations associated with idealization and the applicability of the method. This method also more general technology that takes into account the different scale processes. The construction of basis functions occurs independently for each local domain and doesn't require information exchanges between processors, and the method is effective in a high parallelization. Using constructed multiscale basis functions, we construct a mathematical model on a coarse grid that allows to significant reduce the solution time, the amount of used memory, and can be used to perform calculations for a given configuration of heterogeneous properties.
In this paper, we present numerical results for 2D small PWR test. We compare the GMsFEM results with the fine-grid simulations. In these calculations, we use various number of multiscale basis functions in each local subdomain. As can be
seen, as the number of basis functions increases, the multiscale solution becomes more accurate. In particular, when we use few functions, the result is not accurate, therefore it necessary to select more basis functions. The errors of the method are shown in L2 and H1 norms.
2. Problem statement
Let's consider a one-group diffusion approximation [2] for the neutron transport modelling. We consider a bounded 2D domain Q (x = {xi,x2} G Q) with a convex boundary dQ. The one-group transport equation in the diffusion approximation can be written as follows:
IËÉ. _ y- DV0 + £r0= i^i/E/0 + Ac. (1)
v dt Kef f
Here 0(x, t) is neutron flux at point x and time t, c(x, t) — density of source of delayed neutrons, v — effective neutron velocity, D(x,t) — diffusion coefficient, Xr(x, t) — removal cross-section, Xf (x, t) — fission cross-section, ft — effective fraction of delayed neutrons, A — decay constant of the delayed neutron source. The density of source of delayed neutrons is described by the equation
de B
The albedo-type conditions are set at the boundary dQ:
where n — outer normal to the boundary dQ, 7 — albedo constant. System of equations (1) and (2) is supplemented with boundary conditions (3) and corresponding initial conditions
0(x, 0) = 0°(x), e(x, 0)= e0(x). (4)
Let's propose that at the initial time t = 0, the reactor is in steady-state critical condition.
Discretization. Define a uniform grid
w = {tn = nr, n = 0,1,... ,N, tN = T}
and use next notations <f>n = ^(x,tn), en = e(x,tn). We discretize the time derivatives of equation (1) using finite-difference scheme. For approximation in time, we use a fully implicit scheme with time step t . For delayed neutrons source equation we use numerical-analytical method for the construction of approximations in time. The equation (2) in the equivalent form
dpxte
After integration on time interval [tn,tn+1] one can obtained
tn+1
cn+1 = e-ATcn + 13 J eA(t-t"+1)vE/9$g dt. (5)
tn
For spatial approximation we use the finite element method [15]. To specify the variational formulation, we multiply equations by the test function q and integrate over the domain Q. Using the integration by parts, we obtain the following variational formulation: let's find $n+1 G V such that
/1 $n+1 r
-qdx+ / Dn+1V<j)n+1 ■ Vqdx
n n
n n
/(' 1 $n i' \rn
ry<j)n+1q ds = -^-qdx+ -¡—¡-qdx Vq G V, (6)
dn n n
where V = H 1(Q) is the Sobolev space consisting of scalar functions u such that u2 and |Vu2| have a finite integral in Q.
Further, it's necessary to pass from the continuous variational problem (6) to the discrete problem. We introduce finite-dimensional space of finite elements Vh C V and formulate a discrete variational problem. We use standard linear basis functions as basis functions to solve the problem on the fine grid and
Nf
$ = 53
i= 1
where Nf is the number of vertices on the fine grid.
The problem is solving a system of linear algebraic equations
Af $ = bf, (7)
where the operator Af corresponds to the left side of equation (6), and the vector bf corresponds to the right side of equation (6).
3. Multiscale method
For the discretization on the coarse grid we use GMsFEM [13]. We construct two grids: fine grid (<^h) and coarse grid (JH) (see Fig. 1). For multiscale basis construction, we define local domains wi, where i = 1,..., Nv and Nv is the number of coarse grid nodes. We assume that ¿Ph is a refinement of JH, where h and H represent the fine and coarse grid sizes, respectively. We assume that the fine-scale grid ¿Ph is sufficiently fine to fully resolve the small-scale information of the domain
Fig. 1. Coarse grid and local domain wi with Kj.
while is a coarse grid containing many fine-scale features. A local domain ^ is obtained by the combining all the coarse cells around one vertex of the coarse grid. We construct the multiscale function space
VOff = span{yj }f=1,
where N is the number of coarse basis functions. Each y is supported in local domain w.
Basis functions are designed to capture the multiscale features of the solution. Important multiscale features of the solution are incorporated into localized basis functions which contain information about the scales that are smaller (as well as larger) than the local numerical scale defined by the basis functions.
Multiscale space. The computation of basis functions use local spectral problems to reduce the dimension of the local problem. In order to construct conforming basis functions, we multiply eigenvectors related to dominant eigenvalues to the partition of unity functions. We use following spectral problem in ^
= A S<p*,
(8)
where the elements of the matrices A = {aij} and S = {sij} are defined as follow
+ at - p
I
DV0 •Vqdx + S r 0q dx —
J Tir(j)q dx — J Y
I
Keff (1 + AT )
vSf 4>q dx,
(9)
Duq dx.
Then, we choose eigenvectors corresponding to dominant Mi eigenvalues from (8) and use them to construct the multiscale basis functions.
As partition of unity functions, we use linear functions in each domain Ui. Partitions of unity are calculated in the domain Kj as a linear function from T to the vertex A, and 0 is assigned to the entire segment T, and at point A is assigned the value 1. Thus, we obtain a linear function from 0 to 1 over the entire domain
r r
Fig. 2. Partition of unity functions on the ^ (right) and Kj (left).
Kj. Partitions of unity are shown in Fig. 2. Domain Kj is one element from a coarse grid.
The multiscale space is defined as the span of yi = Xi , where Xi is the usual nodal basis function for the node i (linear partition of unity functions). The number of bases can be different, the accuracy of the solution can be improved when we increase the number of bases (see [12]).
Coarse-scale approximation. Next, we create the following matrix for each
Ui
Ri = [yi,... ,yMi-1 ,VMi].
and define the transition matrix R (transition from a fine grid to a coarse grid) to reduce the dimension of the problem
R = [R1, R2,..., RNv],
where Nv is the number of local domains ui.
Then using the transition matrix R and fine grid system (7), we construct the coarse grid approximation
Ac фс = bc, Ac = RAf RT and hc = Rbf, (10)
and using the coarse-scale solution фc, we can reconstruct the fine grid solution
ф.ms = RT фc.
4. Numerical results
Let's consider the 2D test problem for small PWR reactor [16] (Q is the reactor core area). The geometrical model of the small PWR reactor core is presented in Fig. 3. The diameter of the fuel rod is 0.82 cm, the cell width is 1.26 cm. Diffusion neutronics constants in the common notations are given in Table 1. There are two types of cassettes, with fuel 1% UO2 and 2% UO2. The reflective boundary
Fig. 3. Geometrical model of the small PWR-2D reactor core.
conditions (3) are used (7 = 0). The following delayed neutrons parameters are used: 3 = 6.5 • 10-3, A = 0.08 s-1 and v = 5 • 105 cm/s.
We define the next scenario of the process:
• The A-spectral problem is solved and the solution is taken as the initial condition;
• The non-stationary model is calculated until t = 0.4 sec;
• At t = 0.1 sec and t = 0.3 sec Xr for fuel in the zone 1 changes to +2% and -3%, respectively (simulation of insertion or withdrawal of control rods).
At each time step the integrated power is calculated as
where a is the normalization coefficient, which corresponds to a given value of the integrated power.
The coarse grid contains 49 vertices. The fine grid contains 115891 vertices. The time step for both grids is t = 0.001. As an exact solution, we take the fine-grid solution. The initial value of Keff is 1.183280.
The integral power for the fine grid and relative errors of integral powers are shown in Fig. 4. When using single basis, the error does not exceed 1%, and for using 4 or more bases it does not exceed 0.01%.
In Fig. 5, we present relative L2 and Hi errors of the multiscale solution vs. time for different number of multiscale basis functions. The numerical results show good convergence behaviour, provided that we take sufficient number of the multiscale basis functions.
In Table 2, we present relative L2 and H1 errors at final time for different number of the multiscale basis functions. For example, when we use 8 spectral basis functions, we obtain 0.34% for L2 error and 7.18% for H1 error. Our calculations show that it is necessary to use 4 or more basis functions. In Fig. 6, we present first four (of 16) multiscale basis functions in local domain u. The fine-grid solution and the multiscale solution (16 basis functions on each local domain) ui are shown in Fig. 7. Relative errors are 0.13% for L2 and 4.43% for H1.
Q
Table 1. Diffusion neutronics constants for the small PWR-2D reactor core
Zone 1 2
coolant fuel coolant fuel
D 0.34473872445 0.77002585054 0.31679441560 0.80236505122
Er 5.3858400E-03 8.9337900E-02 6.0670900E-03 6.6279500E-02
E/ 0.0 5.4731800E-02 0.0 3.3377700E-02
V 0.0 2.44844861671 0.0 2.45482762443
Fig. 4. Integral power (fine grid) and relative errors (%) of the multiscale solution power.
Fig. 5. Relative L2 and H1 errors (%) of the multiscale solution.
5. Conclusions
A Generalized Multiscale Finite Element method was developed successfully for modelling neutron transport in one-group diffusion approximation. We presented an implementation of GMsFEM. We considered each step of GMsFEM algorithm. The results showed that GMsFEM performed with a good accuracy in all considered cases.
In the current work, we considered the most popular and simplest model of neutron transport equation. Computational expenses are always an issue even for modern computers. In the future, we will consider more complex models of neutron transport, such as SP3 approximation.
Table 2. Relative L2 and H1 errors (%) of the solution at final time.
Number of bases Number of DOF L2 error H1 error Calc time
1 49 8.09 34.36 0.015
2 98 4.32 26.73 0.018
4 196 0.56 9.24 0.026
8 392 0.34 7.18 0.056
12 588 0.17 5.17 0.102
16 784 0.13 4.43 0.239
fine 115891 - - 6.816
Fig. 6. The first four multiscale basis functions.
flux
909. 915 920 925 930 935 940 945 950 955 960 965 970 975 902.
Fig. 7. Fine grid and multiscale solution with 16 basis functions at final time.
ЛИТЕРАТУРА
1. Duderstadt J. J. Nuclear reactor analysis. New York: Wiley, 1976.
2. Stacey W. M. Nuclear reactor physics. Weinheim : wiley-vch, 2007. V. 2.
3. Avvakumov A. V. et al. Spectral properties of dynamic processes in a nuclear reactor // Ann. Nuclear Energy. 2017. V. 99. P. 68-79.
4. Verd G. et al. 3D A-modes of the neutron-diffusion equation // Ann. Nuclear Energy. 1994. V. 21. N 7. P. 405-421.
5. Munoz-Cobo J. L. et al. 3D calculation of the lambda eigenvalues and eigenmodes of the two-group neutron diffusion equation by coarse-mesh nodal methods // Progr. Nuclear Energy. 2019. V. 110. P. 393-409.
6. Dulla S., Mund E. H., Ravetto P. The quasi-static method revisited // Progr. Nuclear Energy. 2008. V. 50. N 8. P. 908-920.
7. Carre o A. et al. Modal methods for the neutron diffusion equation using different spatial
modes // Progr. Nuclear Energy. 2019. V. 115. P. 181-193.
8. Avvakumov A. V. et al. State change modal method for numerical simulation of dynamic processes in a nuclear reactor // Progr. Nuclear Energy. 2018. V. 106. P. 240-261.
9. Vasilyeva M., Stalnov D. A. Numerical averaging for the heat conduction problem in inhomogeneous and perforated media // Herald of MK Ammosov Northeastern Federal University. 2017.
10. Bakhvalov N. S., Panasenko G. Homogenisation: averaging processes in periodic media: mathematical problems in the mechanics of composite materials. Springer Sci. & Business Media, 2012. V. 36.
11. Vidal-Ferra ndiz A. et al. Pin-wise homogenization for SPN neutron transport approximation using the finite element method // J. Comput. Appl. Math. 2018. V. 330. P. 806-821.
12. Efendiev Y., Hou T. Y. Multiscale finite element methods: theory and applications. Springer Sci. & Business Media, 2009. V. 4.
13. Efendiev Y., Galvis J., Hou T. Y. Generalized multiscale finite element methods (GMsFEM) // J. Comput. Physics. 2013. V. 251. P. 116-135.
14. Spiridonov D., Vasilyeva M., Leung W. T. A Generalized Multiscale Finite Element Method (GMsFEM) for perforated domain flows with Robin boundary conditions //J. Comput. Appl. Math. 2019. V. 357. P. 319-328.
15. Brenner S., Scott R. The mathematical theory of finite element methods. Springer Sci. & Business Media, 2007. V. 15.
16. Forget B., Rahnema F. 3-D heterogeneous coarse-mesh transport method. France: N. p., 2005.
Submitted January 18, 2021 Revised April 22, 2021 Accepted May 26, 2021
Aleksandr O. Vasilev
International Research Laboratory "Multisclae Model Reduction" Ammosov North-Eastern Federal University, 42 Kulakovsky Street, Yakutsk, 677980, Russia haska87@gmail.com Denis A. Spiridonov
International Research Laboratory "Multisclae Model Reduction" Ammosov North-Eastern Federal University, 42 Kulakovsky Street, Yakutsk, 677980, Russia d.stalnov@mail.ru Alexander V. Avvakumov
National Research Center "Kurchatov Institute" Moscow, Russia 1, Akademika Kurchatova pl., Moscow, 123182, Russia avvakumov2009@rambler.ru
