Nuclear Science and Techniques

NUCLEAR ENERGY SCIENCE AND ENGINEERING

A Meshless Local Petrov-Galerkin method for solving the neutron diffusion equation

Shima Tayefi，

Ali Pazirandeh，

Mohsen Kheradmand Saadi

Nuclear Science and Techniques

Vol.29, No.11

Article number 169

Published in print 01 Nov 2018

Available online 16 Oct 2018

DOI：10.1007/s41365-018-0506-x

8100

The goal of this study is to solve the neutron diffusion equation by using a meshless method and evaluate its performance compared to traditional methods. This paper proposes a novel method based on coupling the meshless local Petrov-Galerkin approach and the moving least squares approximation. This computational procedure consists of two main steps. The first involved applying the moving least squares approximation to construct the shape function based on the problem domain. Then, the obtained shape function was used in the meshless local Petrov-Galerkin method to solve the neutron diffusion equation. Because the meshless method is based on eliminating the mesh-based topologies, the problem domain was represented by a set of arbitrarily distributed nodes. There is no need to use meshes or elements for field variable interpolation. The process of node generation is simply and fully automated, which can save time. As this method is a local weak form, it does not require any background integration cells and all integrations are performed locally over small quadrature domains. To evaluate the proposed method, several problems were considered. The results were compared with those obtained from the analytical solution and a Galerkin finite element method (GFEM). In addition, the proposed method was used to solve neutronic calculations in the small modular reactor. The results were compared to those of the citation code and reference values. The accuracy and precision of the proposed method were acceptable. Additionally, adding the number of nodes and selecting an appropriate weight function improved the performance of the Meshless local Petrov-Galerkin method. Therefore, the proposed method represents an accurate and alternative method for calculating core neutronic parameters.

Neutron diffusion equationMeshless local Petrov-Galerkin (MLPG)Moving least squares approximation (MLSA)Meshless methods

1 Introduction

The numerical tools of Partial Differential Equations (PDE) play a key role in the accuracy and acceleration of the solution process in engineering problems. In the last few years, significant progress has been achieved in this field due to the development of computer technology. Despite this progress, numerical methods are often generalized to solve differential equations in a heuristic manner [1]. Numerical methods consist of two types based on the discretization method. The first employs a mesh (mesh-based method) and the second uses a local approximation by polynomials (meshless method).

Mesh-based numerical methods are divided into three classical groups for solving PDEs: Finite Difference Methods (FDM), Finite Volume Methods (FVM), and Finite Element Methods (FEM) [2]. These methods have been used to solve engineering problems for more than 60 years and cover a wide scope of problems such as material science, aerospace engineering, computational physics, safety, reliability analysis, solid mechanic, fluid mechanical interactions, molecular interaction, etc. Mesh-based methods are powerful numerical tools for solving problems and evaluating, predicting, and simulating the behavior of systems. However, mesh-based methods have some limitations in that they cannot provide fully acceptable results. Moreover, despite the great progress in mesh generation, creating an accurate mesh remains difficult. In addition, mesh-based methods are typically time consuming [1].

Therefore, there is an interest in developing alternative methods that remove or reduce the need for meshing. A review of these methods shows that meshless methods have achieved considerable success in recent decades and have attracted substantial interest from researchers because of their applications for solving computational problems [3–9]. The most outstanding of these methods are introduced below.

The advent of the meshless concept dates back to the model of astrophysical phenomena without boundaries (Smooth Particle Hydrodynamics (SPH)); for example, exploding stars and dust clouds. This idea was proposed as a theory for non-spherical stars by Gingold and Monaghan [10]. The SPH scheme has been successfully applied to many problems such as heat conduction, explosion phenomena, and free surface [11, 12]. Then, Nayroles et al. (1992) proposed the diffuse element method that has been employed by many researchers. Belytschko et al. (2004) introduced an extended version of Nayroles's method, which included a series of improvements over the diffuse element method formulation (the element-less Galerkin methods) [14]. Liu et al. (1995, 1993) introduced a new method for correcting the lack of consistency in the SPH scheme, known as the reproducing kernel particle method.

Oñate et al. (1996) studied the finite point manner to model elasticity and plate bending and fluid flow problems using the collocation point technique, least squares approximation, and weighted least squares approximation. Since then, research has tended to use the radial basis function approximation technique [18–21]. Several studies on radial basis functions have been applied to solve physical problems such as transport phenomena [22], heat conduction [23], neutron diffusion [24, 25], and analysis of Kirchhoff plates [26]. De and Bathe (2000) proposed the finite spheres approach for constructing the approximation function. Moradipour and Yousefi (2016) used a meshless kernel-based method to solve the Black–Scholes variational inequality of American options. In addition to the analysis of two-dimensional elasticity problems, Ebrahimnejad et al. (2015) presented three schemes of the 2D meshless finite volume method. One of the main methods of meshless schemes is the Meshless Local Petrov-Galerkin (MLPG) method proposed by [30]. The MLPG method has frequently been used to successfully solve many engineering problems [31–35].

This method is based on the idea of the local weak form, which eliminates the requirement of the background cell. In addition, this method provides numerical integration in a meshless sense to simplify the integrand of the weak form [4, 6, 30, 36–38]. In this global weak form method, a rational basis is provided to construct meshless methods with a larger degree of flexibility [36, 39–44]. The MLPG method is an effective method for solving many problems, which is based on a local weak form and the Moving Least Squares Approximation (MLSA) to approximate the shape functions [45, 46]. As this method does not need a "finite element mesh", it is a truly meshless method [2]. All integrals in the formulation can be easily evaluated over regularly shaped domains and their boundaries [30]. The shape functions must be effectively constructed to achieve a desirable order of continuity.

Accordingly, the capability and flexibility of meshless methods provide the motivation for using this approach to solve the neutron diffusion equation in this study.

2 Methodology

In this research, the MLPG method was used to solve the neutron diffusion equation in two-dimensional geometry. First, the solution was approximated by using a local approximation in the problem domain. This approximation scheme was used to formulate the residual of the governing equations on the problem and find the shape function. The obtained shape function was employed in the MLPG algorithm to solve the problem. Then, the problem was modeled by the analytical solutions and the Galerkin finite element method (GFEM) method to evaluate the proposed method. In the following subsections, all steps required to derive the shape function and the weak form of the problem equation are explained. A schematic representation of the local support domains used to create the shape function is demonstrated in Fig. 1. The domain representation of the assumed problem of the mesh-based and meshless approaches is shown in Fig. 2.

Fig. 1

(Color online) A schematic representation of the local support domains in the meshless approach

Fig. 2

(Color online) Comparison of a problem domain in a meshless approach and a mesh-based method

2.1 Determining the shape function by MLSA

Selection of an approximation function has an important role for reaching a sufficient number of equations to determine the unknowns [47, 48]. Some of the most important functions consist of Shepard functions [40, 49], least squares [14, 50], local radial point interpolation [51], and local point interpolation [52]. One of the best schemes for interpolating data with an appropriate accuracy is the moving least squares method, which has more compatibility with MLPG method [53]. The moving least squares approximation was created by mathematic researchers for data fitting and surface construction [54]. The steps for determination of the shape function are as follows:

1) First, to demonstrate the trial function, the meshless approach should use local interpolation. Therefore, a function of the field unknown variables, U(X), is defined in a subdomain, Ω_s, which indicates the support domain of the MLS approximation for U(X). It should be mentioned that the support domain is located in the problem domain Ω_p. X=(x, y), which indicates a two-dimensional problem. To find an approximation distribution of function U in the support domain, the local moving squares approximation of function U(X) for all X in the support domain should be defined as follows [42]:

u^{h} (X) = \sum_{j = 1}^{m} p_{j} {(X)a}_{j} {(X)=P}^{T} (X)a(X) .

(1)

P^T(X) is a complete monomial basis of order m. m is the number of monomial terms, here the linear basis m=3.

P^{T} {(X)=P}^{T} (x, y) = {1, x, y}

(2)

a(X) is a coefficient vector as follow:

a^{T} {(X)=a}^{T} (x, y) = {a_{0} (X), a_{1} (X), \dots, a_{m} (X)} .

(3)

The approximated values of the field function are:

u^{h} {(X,X}_{I}) = {P(X}_{I}) a(X) I = 1, .., n .

(4)

2) In the second step, the weighted discrete L₂ norm is defined as follows:

J(X) = \sum_{(i = 1)}^{n} {W(X-X}_{i}) {[U^{h} {(X,X}_{i}) - U_{i} (X_{i})]}^{2} = {[P^{T} {(X}_{i}) a (X) - U_{i}]}^{2} .

(5)

U_i is the nodal parameter of the field variable at node i. n indicates the number of neighborhood nodes of node i. The weight function $W (X - X_{i})$ is the associated weight function of node i. It is noted that the weight function for all internal and external nodes in the support domain is non-zero and zero, respectively.

3) In the third step, for determining the coefficients matrix, a(X), J should be minimized with respect to a(X) [42]:

\frac{\partial J}{\partial a} = 0.

(6)

This equation leads to a linear relationship (7) between a(X) and $U$ :

{A(X)a(X)=B(X)U}_{s} .

(7)

Solving Eq.(7) for a(X) gives [54]:

{a(X)=A(X)}^{-1} {B(X)U}_{s},

(8)

where U_s is a vector for choosing the field function nodal parameters for interior nodes of the support domain as follows:

U_{s} = {U_{1} {,U}_{2}, \dots ， U_{n}}^{T}

(9)

The matrices A(X) and B(X) are given as [54, 55]:

A (X) = P^{T} W P = \sum_{i = 1}^{n} W_{i} (X - X_{i}) P (X_{i}) P^{T} (X_{i})

(10)

B (X)= P^{T} W =[W_{1} (X) P (X_{1}), \dots, W_{n} (X) P (X_{n})]

(11)

where $W i (X) = W (X - X i)$ .

Hence, for this two-dimensional problem, A(X) is a symmetric matrix according to the linear basis (m=3) that can be explicitly written as:

(12)

少公式(12)

4) The fourth step is the derivation of the shape function. If the above equations are substituted into Eq. (4), the form of an interpolation function that is MLSA of u(X) at X is written as:

u^{h} (X) = \sum_{j = 1}^{m} ϕ_{i} (X) U_{i} (X) = Φ^{T} (X) U_{s} .

(13)

$Φ^{T} (X)$ can be rewritten as below [55]:

Φ^{T} (X) = {ϕ_{1} (X) ϕ_{2} (X) \dots \dots . ϕ_{n} (X)}_{1 * n} = p^{T} (X) A^{- 1} (X) B (X),

(14)

where $Φ (X)$ indicates the MLS shape function vector for n nodes in the support domain of point X. Therefore, the shape function for the ith node is introduced in Eq. (15)[55]:

ϕ_{i} (X) = \sum_{j = 1}^{m} p_{j} (X) {(A^{- 1} (X) B (X))}_{j i} .

(15)

A flowchart of the algorithm for the above steps is shown in Fig. 3.

Fig. 3

Flowchart of the applied algorithm for deriving the shape function (MLSA algorithm)

2.2 Solving the neutron diffusion equation by MLPG

To solve the neutron diffusion equation, the nodal parameters u_j and its derivatives should be defined. u_j is a function of time, while the shape function $ϕ_{j}$ is a function of X only [56]. The spatial and temporal derivatives of the above shape function with respect to X or Y are extracted from the above equations. Equation (13) can then be rewritten as:

u^{h} (X, t) = \sum_{j = 1}^{n} u_{j} (t) ϕ_{j} (X),

(16)

2.2.1 The temporal derivation of $u^{h}$ is:

u_{t}^{h} = \sum_{j = 1}^{n} \overset{´}{u_{j}} (t) ϕ_{j} (X) .

(17)

2.2.2 The spatial derivation of $u^{h}$ is:

u_{x}^{h} = \sum_{j = 1}^{n} u_{j} (t) \overset{´}{ϕ_{j}} (X) .

(18)

2.2.3 In addition, the second-order spatial domain is:

u_{x x}^{h} = \sum_{j = 1}^{n} u_{j} (t) \ddot{ϕ} j (X) .

(19)

The derivatives are denoted by subscripts. In addition, a local quadrature domain, $Ω_{Q_{i}}$ is defined for each interior node, x_i. R indicates the residual function. The integral of the weighted residual is zero. [57].

\int_{Ω_{Q_{i}}}^{} R W_{i} d Ω_{Q_{i}} = 0

(20)

As mentioned earlier, $W_{i}$ is a weight function that is non-zero only on $Ω_{Q i}$ and similar to the weight function used in the MLSA process.

2.3 Determining the residual of the neutron diffusion equation

In this step, two-dimensional neutron diffusion equations are first introduced. For arbitrary volume, v, the balance neutron equation is considered. This well-known equation is shown below.

[Time rate of change of the number of neutron in a volume]

=[Production rate in the volume]

- [Absorptions in the volume]

- [Net leakage from the surface of volume]

The first term is mathematically expressed in Eq. (21):

\frac{d}{d t} \int_{v}^{} \frac{1}{v} φ (r, t) d v

(21)

The production rate is considered as a neutron source in Eq. (22):

\int_{v}^{} S (r, t) d v

(22)

The neutron absorption term is presented in Eq. (23):

\int_{v}^{} Σ_{a} φ (r, t) d v .

(23)

The leakage term is shown in Eq. (24):

\int_{v}^{} \nabla . J d v = \frac{d}{d t} \int_{v}^{} \nabla . [D \nabla φ (r, t)] d v .

(24)

The mathematical formulation of the neutron diffusion theory is obtained as follows:

\frac{\partial φ}{\partial t} - D \nabla^{2} φ (r, t) + Σ_{a} φ (r, t) = S (r, t) .

(25)

Equation (25) can be rewritten as equation (26):

u_{t} - D u_{x x} + Σ_{a} u - S = 0.

(26)

Where u is the solution, D is the diffusion coefficient, s is a neutron source, and $Σ_{a}$ is the absorbtion cross section. The residual for this equation is defined according to an approximation solution as follows:

Residual = R = u_{t}^{h} - D u_{x x}^{h} + Σ_{a} u^{h} - s = 0.

(27)

Equation (27) is substituted into equation (20) to give a system of n linear equations for n nodes:

\int_{Ω q_{i}}^{} [u_{t}^{h} - D u_{x x}^{h} + Σ_{a} u^{h} - s] w_{i} d Ω_{Q q_{i}} = 0

(28)

Then, each of the neutron diffusion equation terms should be substituted by Eq.(16) into Eq.(19):

a) The first term is obtained by substituting Eq.(17) into Eq. (28):

\int_{Ω q_{i}}^{} [u_{t}^{h}] w_{i} d Ω_{Q q_{i}} = \int \sum_{j = 1}^{n} \overset{´}{u_{j}} φ_{j} (x) w_{i} (x) d x .

(29)

b) The second term is obtained by substituting Eq. (19) into Eq. (28):

\begin{array}{l} - D \int_{Ω q_{i}}^{} [u_{x x}^{h}] w_{i} d Ω_{q_{i}} = \int_{Ω q_{i}}^{} D u_{x}^{h} w_{i, X} (x) d Ω_{q_{i}} - D u_{x}^{h} w_{i} | \begin{matrix} X q i n \\ X q i 1 \end{matrix} \\ = D \int \sum_{j = 1}^{n} u_{j} (t) {\overset{´}{φ}}_{j} (x) {\overset{´}{w}}_{i} (X) d x - D u_{j} (t) φ_{j} (x) w_{i} (x) . \end{array}

(30)

c) The third term is obtained by substituting Eq. (16) into Eq. (28):

\int_{Ω q_{i}}^{} Σ_{a} u^{h} w_{i} d Ω_{q_{i}} = \int \sum_{j = 1}^{n} Σ_{a} u_{j}^{h} (t) φ_{j} (X) w_{i} (x) d x .

(31)

d) The fourth term is obtained by substituting the fission effective factors into Eq. (22):

\int_{Ω q_{i}}^{} - S w_{i} d Ω_{q_{i}} = - \int \frac{χ}{λ} υ Σ_{f} w_{i} (x) d x .

(32)

Therefore, Eq. (28) can be rewritten as:

\int \sum_{j = 1}^{n} \overset{´}{u_{j}} φ_{j} (x) w_{i} (x) d x + D \int \sum_{j = 1}^{n} u_{j} (t) {\overset{´}{φ}}_{j} (x) {\overset{´}{w}}_{i} (X) d x - D \sum_{j = 1}^{n} u_{j} (t) φ_{j} (x) w_{i} (x) | \begin{matrix} X q i n \\ X q i 1 \end{matrix} + \int \sum_{j = 1}^{n} Σ_{a} u_{j}^{h} (t) φ_{j} (X) w_{i} (x) dx - \int \frac{χ}{λ} {υΣ}_{f} w_{i} (x) dx = 0.

(33)

For clarity, the sums are moved outside the integrals as shown in equation (34):

\sum_{j = 1}^{n} \int \overset{´}{u_{j}} (t) φ_{j} (x) w_{i} (x) d x + D \sum_{j = 1}^{n} \int u_{j} (t) {\overset{´}{φ}}_{j} (x) {\overset{´}{w}}_{i} (x) d x - D \sum_{j = 1}^{n} u_{j} (t) φ_{j} (x) w_{i} (x) | \begin{matrix} X q i n \\ X q i 1 \end{matrix} + \sum_{j = 1}^{n} \int Σ_{a} u_{j} (t) φ_{j} (x) w_{i} (x) d x - \int \frac{χ}{λ} {υΣ}_{f} w_{i} (x) dx = 0.

(34)

This equation can be rearranged into its final form as shown in Eq.(35):

\sum_{j = 1}^{n} \overset{´}{u_{j}} (t) \int φ_{j} (x) w_{i} (x) d x + D \sum_{j = 1}^{n} u_{j} (t) \int {\overset{´}{φ}}_{j} (x) {\overset{´}{w}}_{i} (x) d x - D \sum_{j = 1}^{n} u_{j} (t) φ_{j} (x) w_{i} (x) | \begin{matrix} X q i n \\ X q i 1 \end{matrix} + \sum_{j = 1}^{n} u_{j} (t) \int Σ_{a} φ_{j} (x) w_{i} (x) d x - \frac{χ}{λ} {υΣ}_{f} \int w_{i} (x) dx = 0

(35)

It should be noted that the integrals obtained in discretized equations are regular integrals. To calculate the integrals, the numerical integration method based on the 20 point Gauss–Legendre quadrature was applied [33, 45, 58].

2.4 Solution of the time-dependent neutron diffusion equation

The first stage in the time-dependent solution is the discretization of space and time variables. The space and time variables are essentially discretized to gain a series of (nonlinear) algebraic systems. Therefore, it is necessary to decouple the space and time variables to solve this problem.

In this study, three different approaches to time discretization were applied, the forward Euler, Crank-Nicolson, and backward Euler scheme. It is notable that the results of the second scheme were more acceptable than other methods because the method is unconditionally stable. In addition, the time required for this process was shorter than other methods. In the method, time is replaced by the time derivative at the half-step ${\overset{´}{u}}_{j} (t + Δ t / 2)$ with the central difference approximation [59]:

{\overset{´}{u}}_{j} (t + Δ t / 2) = \frac{u_{j} (t + Δ t) - u_{j} (t)}{Δ t} .

(36)

In this discretization method, $u_{j} (t + Δ t / 2)$ is approximated as the average of $u_{j} (t)$ and $u_{j} (t + Δ t)$ . It is noted that the second-order accuracy in time and space are provided by the Crank-Nicolson method. In addition, this method is unconditionally stable.

By collecting all the nodal parameters into a vector $u_{j}^{(t)}$ = ( $u_{1}^{(t)} \dots .. u_{n}^{(t)}$ ), the matrix form of Eq. (35) is derived as below:

M \sum_{j = 1}^{n} \overset{´}{u_{j}} (t) + N \sum_{j = 1}^{n} u_{j} (t) - S = 0,

(37)

M \sum_{j = 1}^{n} \overset{´}{u_{j}} (t + \frac{Δ t}{2}) + N \sum_{j = 1}^{n} u_{j} (t + \frac{Δ t}{2}) - S = 0,

(38)

M (\frac{u (t + Δ t) - u (t)}{Δ t}) + N (\frac{u (t + Δ t) + u (t)}{2}) - S = 0,

(39)

2 M u (t + Δ t) - 2 M u (t + Δ t) + N Δ t u (t) + N Δ t u (t + Δ t) - 2 t Δ t = 0,

(40)

u^{(t + Δ t)} (2 M + N Δ t) = 2 t Δ t + u^{(t)} (2 M - N Δ t),

(41)

where M and N are coefficients of $\overset{´}{u_{j}}$ and $u_{j}$ , respectively. The proposed final relationship, which was used in the computational algorithm as the main formula at each time-step, is:

u^{(t + Δ t)} L = s^{*} + u^{(t)} K .

(42)

All approaches were calculated by a code written entirely in FORTRAN 95. To enforce the boundary conditions in the problem, the following provisions should be considered. Each element of matrix L, L_ij, is defined based on the three conditions. If x_j is in the support domain of boundary node x_i, L_ij is $ϕ_{j} (x_{i})$ . If x_j is in the support domain of any point in the quadrature domain of the internal node x_i, L_ij is $2 M_{i j} + N_{i j} Δ t$ . Otherwise, the value of L_ij is zero.

In addition, each element of matrix K, K_ij, is defined based on the three conditions. If x_j is in the support domain of boundary node x_i, L_ij is zero. If x_j is in the support domain of any point in the quadrature domain of the internal node x_i, L_ij is $2 M_{i j} - N_{i j} Δ t$ . Otherwise, the value of L_ij is zero and $S^{*}$ is based on the problem conditions.

As mentioned previously, the shape function is an important requirement for applying meshless methods. A qualified approach should satisfy several conditions to provide the shape functions of meshless methods [48]. These conditions consist of having an arbitrary nodal distribution, being numerically stable, satisfying a certain order of consistency, and being compactly supported [48]. The accuracy of interpolation for the arbitrary points depends on the number of nodes in the support domain [53]. The number of nodes in each interpolation domain should be used to determine the appropriate shape function and local interpolation schemes to match the sensitivity of any meshless interpolation methods to a variable. Therefore, an appropriate support domain should be selected to apply an accurate approximation [4]. It should be noted that, in order to evaluate the performance of the proposed method, several weight functions such as the cubic spline, the quadratic spline, the compact support radial basis, and Gaussian weight functions were studied. These weight functions are presented below:

1. Cubic spline weight function

\begin{array}{l} w_{i} (x) = \frac{2}{3} - 4 r_{i}^{2} + 4 r_{i}^{3}, for 0 \leq r_{i} \leq \frac{1}{2} \\ = \frac{4}{3} - 4 r_{i} + 4 r_{i}^{2} - \frac{4}{3} r_{i}^{3}, for \frac{1}{2} \leq r_{i} \leq 1 \\ = 0, for r_{i} > 1 \end{array}

(43)

2. Quadratic spline weight function

\begin{array}{l} w_{i} (x) = 1 - 6 r_{i}^{2} + 8 r_{i}^{3} - 3 r_{i}^{4} for 0 \leq r_{i} \leq 1 \\ = 0, for r_{i} > 1 \end{array}

(44)

3. Compact support radial basis function (CSRBF)

\begin{array}{l} w_{i} (x) = {(1 - r_{i})}^{6} + [6 + 36 r_{i} + 82 r_{i}^{2} + 72 r_{i}^{3} + 30 r_{i}^{4} + 5 r_{i}^{5}], for 0 \leq r_{i} \leq 1 \\ = 0, for r_{i} > 1 \end{array}

(45)

4. Gaussian weight function

w_{i} (x) = {\begin{matrix} \frac{\exp [- {(\frac{d_{i}}{c_{i}})}^{2}] - exp [- {(\frac{r_{i}}{c_{i}})}^{2}]}{1 - exp [- {(\frac{r_{i}}{c_{i}})}^{2}]} f o r 0 \leq d_{i} \leq r_{i} \\ 0 f o r d_{i} \geq r_{i} \end{matrix}

(46)

Where $r_{i}$ is the size of the support for the weight function, w_i, and d_i=│x-x_i│is the distance between node x_i and x. The distance between nodes is equal to half the mean free path of the neutron. The value of the mean free path for a thermal neutron in a reactor core is 0.7. In addition, the radius of the support domain for the weight function is approximately three mean free paths of the neutron. Moreover, $c_{i}$ is a constant controlling the shape of the weight function, w_i., which is an effective parameter for the stability and accuracy of the method that is obtained through experiments. In this problem, $c_{i}$ in the range of 0.3-0.8 leads to more stable and accurate results. Then, for each interior node, x_i, a local quadrature domain, $Ω_{Q_{i}}$ , is defined [60].

The main steps of the proposed algorithm include:

1. Defining the initial condition of the problem;

2. Defining the boundary condition of the problem;

3. Setting up the cloud points;

4. Finding the affect domain of each node;

5. Looping over the time steps;

6. Appling the MLSA algorithm:

6.1 Selecting the neighboring nodes;

6.2 Deriving the shape functions and its derivative for the nodes;

6.3 Evaluating the nodal matrices and vectors;

6.4 Setting the nodal portion to the global matrices and vectors;

6.5 End loop (see Fig. 3);

7. Solving the government equations at each node;

8. Appling the moving least squares shape function for interior nodes of a local domain for recalculation;

9. Determining the coefficients at each node to obtain higher accuracy;

10. Determining and evaluating the obtained results;

11. Evaluating the effects of increasing the amount of nodes in the problem domain;

12. Recording the history of the state variables and their derivatives;

13. Back to 5;

14. End.

The flowchart of the MLPG algorithm is demonstrated in Fig. 4.

Fig. 4

Flowchart of the MLPG algorithm

3 Results and Discussion

Several example problems were studied to validate the proposed method. The results of the MLPG method were then compared to the results of the analytical solution and the Galerkin finite element method.

3.1 Convergence degree test

In order to calculate the convergence degree of the MLPG method, a slab is considered. The boundary conditions are a perfect reflector on the left and a bare surface on the right of the slab. The model geometry is presented in Fig. 5. The logarithm of the error versus the number of nodes is presented in Fig. 6. It is clear that increasing the number of nodes decreases the error. Fig. 7 illustrates the logarithm of the error versus the logarithm of the inverse of the number of nodes. It is well known that the slope of this diagram depicts the convergence rate. For further clarification, the convergence rate versus the number of nodes is presented in Fig. 7b, which demonstrates that the behavior of the proposed method is stable for less than 500 nodes.

Fig. 5

Geometry of the convergence degree test

Fig. 6

Error logarithm versus the number of nodes

Fig. 7

(a) Logarithm of the error versus the logarithm of the inverse of the number of nodes and (b) the convergence rate versus the number of nodes

3.2 Constant source problem

In this test problem, a distributed source problem was considered. This problem consists of two models, a single region (Fig. 8) and two regions (Fig. 11) including a source in region 1 and a reflector in region 2 with dimensions of 150 cm × 150 cm. The boundary conditions include a perfect reflector on the left side and a bare surface on the right side. In this section, three types of problem were considered:

Fig. 8

(Color online) Geometry of the single region

Fig. 11

(Color online) Geometry of the two-region problem

CASE-I: The geometry of the single region is presented in Fig. 8. The parameters required for this calculation are shown in Table 1. The fast and thermal flux of the MLPG method compared to that of GFEM and the analytical solution is presented in Figs. 9 and 10, respectively.

Applied neutronic parameters for CASE-I

Applied models	Energygroup	Absorptioncross section	Nu-Fissioncross section	Diffusion coefficient
Slab I	1	0.03051	0.000256	1.440130
	2	0.10582	0.154726	0.015723

Fig. 9

Thermal neutron flux obtained from the MLPG and GFEM and the analytical solution (CASE-I)

Fig. 10

Fast neutron flux obtained from the MLPG and GFEM methods and the analytical solution (CASE-I)

CASE-II: The geometry of the two-region problem is presented in Fig. 11 and the required parameters are shown in Table 2. A comparison of the fast and the thermal flux of the MLPG method with that of the GFEM and the analytical solution is presented in Figs. 12 and 13, respectively. As shown, the results of the proposed method exhibit good agreement with those of the analytical solution.

Applied neutronic parameters for CASE-II

Applied models	Energy group	Absorption cross section	Nu-Fission cross section	Diffusion coefficient	Region
Slab-II	1	0.03051	0.000256	1.440130	Fuel
	2	0.10582	0.154726	0.015723
	1	0.035128	0.000000	1.852524	Reflector
	2	0.032540	0.000000	0.294517

Fig. 12

Thermal neutron flux obtained from the MLPG and GFEM and the analytical solution (CASE-II)

Fig. 13

Fast neutron flux obtained from the MLPG and GFEM and the analytical solution (CASE-II)

The average error of these results obtained by applying various weight functions is presented in Table 3, which indicates that the accuracy of the solution is improved by using Gaussian weight functions. In Table 4, the average error and CPU time of this problem are demonstrated for the Gaussian and CSRBF weight functions. Fig. 14 presents a comparison of the fast neutron flux in CASE-I and CASE-II obtained by inserting 50×50 nodes.

Comparison between average errors obtained by applying various weight functions

Error	MLPG Applied Weight Functions
	Gaussian	CSRBF	Cubic Spline	Quadratic Spline
CASE-I	0.0035	0.073	0.094	0.085
CASE-II	0.0038	0.070	0.091	0.083

Comparison between average errors and CPU time obtained by applying Gaussian and CSRBF weight functions

Problem No.	MLPG Applied Weight Functions
	Gaussian		CSRBF
	Error	CPU Time (s)	Error	CPU Time(s)
CASE-I	0.0035	2	0.073	2.5
CASE-II	0.0038	3.5	0.070	6

Fig. 14

(Color online) Comparison of the fast neutron flux in CASE-I (left) and CASE-II (right)

In addition, to evaluate the performance of the number of inserted nodes, MLPG methods with 20×20 nodes and 40×40 nodes were compared. As shown, an accurate value of the thermal flux was achieved by the MLPG with 40×40 nodes for CASE-I (Fig. 15). Fig. 15 indicates the performance of the MLPG method was improved by increasing the number of nodes. Fig. 16 demonstrates the contour of the fast and the thermal neutron flux. The best results for CASE-II were obtained with 40×40 nodes.

Fig. 15

(Color online) Comparison of the effect of the number of nodes on the thermal flux

Fig. 16

(Color online) Contour of the fast (left) and the thermal (right) neutron flux

CASE-III: In addition, in order to evaluate the performance of the proposed method for dealing with curve boundaries, a problem with a curved domain was considered. The geometry of this problem is presented in Fig. 17. A comparison of the contours of the thermal neutron flux of this problem is shown in Fig. 18. Two MLPG methods with various numbers of nodes were employed in this problem. The results show good agreement between the proposed method and analytical solutions. In addition, the obtained results illustrate that, to achieve the required accuracy, a large number of nodes can be inserted in the problem domain. Furthermore, increasing the number of nodes in a curved domain by the MLPG method is much simpler than mesh generation for mesh-based methods. However, the integration process for a curved domain is challenging. As mentioned, a parallel algorithm was used in this study to improve the calculation.

Fig. 17

(Color online) Geometry of CASE-III

Fig. 18

(Color online) Thermal neutron flux contour for CASE-III

CASE-IV: Moreover, in order to evaluate the performance of the proposed method, two types of node distribution, normal and scatter models, were applied to CASE-I. The average error of these results obtained by applying various weight functions is presented in Table 4. Table 5 indicates that the accuracy of the solution was improved by using Gaussian weight functions and the scatter distribution. The patterns of the node distributions are presented in Fig. 19.

Comparison between average errors obtained by applying normal and scatter node distributions for CASE-I

Distribution Pattern	MLPG Applied Weight Functions
	Gaussian	CSRBF	Cubic Spline	Quadratic Spline
Normal Distribution	0.0035	0.073	0.094	0.085
Scatter Distribution	0.0030	0.059	0.082	0.065

Fig. 19

(a) Normal and (b) scatter node distributions

3.3 Void test problem

The void test problem was first studied by Abuzaid and Gashut [61]. The problem involves a slab including multiple regions with a perfect reflector on the left side and a bare surface on the right side. Fig. 20 presents the geometry and properties of the problem. A comparison between the results of the proposed method and the analytical solution is provided in Fig. 21, which indicates a good agreement.

Fig. 20

Geometry of the void test problem

Fig. 21

Comparison of the MLPG method and the analytical solution in the void test problem

3.4 Neutronic calculations for a SMR reactor

The proposed method was applied to calculate the neutronic parameters in a Small Modular Reactor (SMR). The core used for the analysis was Westinghouse’s SMR, with an electrical output of 200 MW and an active core height of 2.4 m. The core consists of 89 assemblies contained within a core barrel and reactor vessel. These assemblies include 52 fuel assemblies and 37 control rod assemblies. The SMR core map is presented in Fig. 22. The reactor vessel components were based on an AP1000 design but modified to a reduced diameter and height of 3.5 m and 24.7 m, respectively [62].

Fig. 22

(Color online) SMR core map

The fuel assembly was a square lattice in a standard 17×17 layout with 264 fuel rod locations, 24 guide tube locations, and 1 central location for instrumentation, incorporating the standard Westinghouse design specifications. The fuel of the core was uranium-oxide (UOX) with less than 5% U-235 enrichment (2.5w/o, 3.5w/o, 4.2%w/o U-235 fuel) and a density of 10.36 g/cm³. The fuel rod consisted of a uniform cylindrical pellet stacked together within a Zircaloy clad tube. Between the fuel stack and the clad, clearance was provided to accommodate fuel swelling due to accumulation of fission products, thereby preventing clad rupture. The gap was filled with helium gas to improve heat conduction from fuel to cladding. The guide tubes in the fuel assembly served as a location for the insertion of the rod cluster control assembly (RCCA). The RCCA is a spider assembly consisting of evenly spaced control rods of either silver indium cadmium or boron carbide rods based on the type of fuel used. Detailed specifications for the fuel rod, cladding, structure, control rod, and burnable poison (i.e. discrete and integral) were taken from CASL (Consortium for Advanced Simulation of LWRs) VERA core physics benchmark specifications [62].

This study reveals the ability of the MLPG method to solve the neutronic diffusion equation and evaluates the performance of this method when different selectable parameters were incorporated into the solution. In order to automate the procedure of the proposed method for this problem, a computer program was developed. The first stage of computational procedures consists of creating the cross section database by WIMSD4 code in three-group energy. Then, some of the main neutronic parameters of the core were calculated by the proposed method. The results of the MLPG method were compared with the results of the citation code and the reference values. In addition, as the selection of the weight function and the number of nodes in the problem domain plays an important role in solution accuracy, the performance of the MLPG method was examined with respect to these parameters. Table 6 demonstrates the effective multiplication factor as a main parameter of the core calculation. Tables 7 and 8 present the maximum to average flux ratio and the control rod worth, respectively. The choice of weight function is clearly a significant factor for obtaining an accurate solution. The most important frequently used weight functions are the cubic spline, quadratic spline, CSRBF, and Gaussian. Table 9 compares the average errors of results obtained from the various weight functions and reveals that the Gaussian function produces better accuracy than other functions.

Comparison of the effective multiplication factor between MLPG methods and citation code

Parameter	Citation	Referenced value	MLPG Applied weight functions				Nodes Insertion
			Gaussian	CSRBF	Cubic Spline	Quadratic Spline
Effective Multiplication Factor	1.12414	1.12511	1.12332	1.11032	1.11054	1.11902	37×37
			1.12415	1.11940	1.11891	1.11974	91×91
			1.12496	1.12439	1.12475	1.11996	289×289

Comparison of the maximum to average flux ratio between MLPG methods and citation code

Parameter	Citation	Referenced value	MLPG Applied weight functions				Nodes Insertion
			Gaussian	CSRBF	Cubic Spline	Quadratic Spline
Maximum to Average Flux Ratio	2.48520	2.4980	2.48213	2.43191	2.44521	2.39120	37×37
			2.48500	2.45225	2.46378	2.44032	91×91
			2.49739	2.48821	2.49707	2.45318	289×289

Comparison of the control rod worth between MLPG methods and citation code

Parameter	Citation	Referenced value	MLPG Applied weight functions				Nodes Insertion
			Gaussian	CSRBF	Cubic spline	Quartic spline
Control Rod Worth	0.15613	0.1571	0.15243	0.14254	0.1550	0.14382	37×37
			0.15615	0.14367	0.1561	0.14630	91×91
			0.15694	0.14975	0.15679	0.14587	289×289

Comparison between average errors of results obtained by applying various weight functions (for 289×289 nodes)

Error	MLPG Applied weight functions
	Gaussian	CSRBF	Cubic spline	Quartic spline
Effective Multiplication Factor	0.001	0.0064	0.032	0.045
Maximum to Average Flux Ratio	0.002	0.041	0.037	0.179
Control Rod Worth	0.001	0.046	0.020	0.071

In addition, the performance of the MLPG method with respect to the number of nodes was evaluated for the above parameters. The values obtained for these parameters are promisingly close to the citation code results and reference values for 289×289 nodes. Moreover, a comparison between the average errors and the CPU time by applying Gaussian weight functions and various node distributions to derive an effective multiplication factor are listed in Table 10. The errors decreases with increasing node number. In addition, the scatter distribution results are more accurate than those of the normal distribution.

Comparison between average errors of results and CPU time obtained by applying Gaussian weight function and various node distributions

Parameter	MLPG Gaussian weight function
	Nodes Normal Distribution			Nodes Scatter Distribution
	Error	CPU Time(s)	Nodes	Error	CPU Time(s)	Nodes
Effective Multiplication Factor	0.012	4.0	37×37	0.010	3.5	37×37
	0.087	5.5	91×91	0.063	4.0	91×91
	0.001	10	289×289	0.001	7.5	289×289

The fast flux and thermal flux of the reactor core is demonstrated in Fig. 23 for the MLPG method with various node numbers. In addition, to evaluate the performance of the MLPG method with respect to the number of nodes, 37×37, 91× 91, and 289×289 nodes in the radial direction and 21 layers in the axial direction were employed in the problem domain for any of the fuel assembly types.

Fig. 23

Comparison of a) the fast flux and b) the thermal flux of the applied MLPG method for 37×37, 91× 91, and 289×289 nodes with the citation code

Fig. 24 shows the power factor distributions for a 1/8 symmetry of this core. In this figure, the reference values are compared with the values of the proposed method (MLPG, 289×289 nodes). The reference value exhibits good agreement with the calculated value. The differences between these parameters are also demonstrated in Fig. 24. The comparison of these methods with the citation code and reference values indicates that increasing the number of nodes leads to a more accurate performance of the MLPG method.

Fig. 24

(Color online) Power factor distributions in a 1/8 symmetry of the SMR core

4 Conclusion

In this work, the meshless local Petrov-Galerkin method and moving least squares approximation were coupled to solve the neutron diffusion equation. Several problems were examined to demonstrate the degree of accuracy that can be obtained by employing this method as well as the applicability of the proposed method with respect to effective factors. The results of the proposed method were compared with the analytical solution and the GFEM method, indicating good agreement. Additionally, the proposed method was used to solve neutronic calculations in the SMR reactor. The results were compared to those obtained from the citation code and reference values.

The accuracy and precision of the proposed method were both acceptable and further improved by increasing the number of nodes and through selection of an appropriate weight function. It should be noted that the run-time varied from 0.5 to 10 s, depending mainly on the weight function type, the number of nodes per cell domain, the number of nodes in the boundary, the pattern of node distribution, the type of boundary conditions applied, and the particular computer system on which the code is run. As mentioned above, the numerical proposed method removes the mesh generation process by inserting nodes into the problem domain. Therefore, this study eliminated the mesh generation time plus the interactive and CPU time and decreased the time costs.

In summary, using the MLPG method to solve the neutron diffusion equations is promising due to improved computational accuracy and execution time. Accordingly, the outstanding performance of the MLPG method has been confirmed for nuclear engineering calculations, indicating that it can be successfully applied to the desired nuclear computational codes.

References

M. Steinhauser, Computational multi scale modeling of fluids and solids. (Springer, Berlin, 2017)

G. R. Liu, Mesh free methods: moving beyond the finite element method. (Taylor & Francis, Boca Raton, 2009)

T. Belytschko, T. Rabczuk, A. Huerta et al.,

Mesh free methods

. in Encyclopedia of Computational Mechanics, edited by E. Stein, R. de Borst, and J. R. Hughes (John Wiley and Sons, Chichester, 2004), pp. 1-48