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

A Meshless Local Petrov-Galerkin method for solving the neutron diffusion equation 1 first-author 1 corresp Corresponding author, Paziran@ut.ac.ir Corresponding author, Paziran@ut.ac.ir Paziran@ut.ac.ir 1 Department of Nuclear Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran Department of Nuclear Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran 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 equation Meshless local Petrov-Galerkin (MLPG) Moving least squares approximation (MLSA) Meshless methods 1 Introduction 1 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 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 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 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 3 – 9 9 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 10 ]. The SPH scheme has been successfully applied to many problems such as heat conduction, explosion phenomena, and free surface [ 11 11 , 12 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 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 18 – 21 21 18–21 ]. Several studies on radial basis functions have been applied to solve physical problems such as transport phenomena [ 22 22 ], heat conduction [ 23 23 ], neutron diffusion [ 24 24 , 25 25 ], and analysis of Kirchhoff plates [ 26 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 30 ]. The MLPG method has frequently been used to successfully solve many engineering problems [ 31 31 – 35 35 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 4 , 6 6 , 30 30 , 36 36 – 38 38 36–38 ]. In this global weak form method, a rational basis is provided to construct meshless methods with a larger degree of flexibility [ 36 36 , 39 39 – 44 44 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 45 , 46 46 ]. As this method does not need a "finite element mesh", it is a truly meshless method [ 2 2 ]. All integrals in the formulation can be easily evaluated over regularly shaped domains and their boundaries [ 30 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 1 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 Fig. 1 . The domain representation of the assumed problem of the mesh-based and meshless approaches is shown in Fig. 2 Fig. 2 . 2.1 Determining the shape function by MLSA 2 Selection of an approximation function has an important role for reaching a sufficient number of equations to determine the unknowns [ 47 47 , 48 48 ]. Some of the most important functions consist of Shepard functions [ 40 40 , 49 49 ], least squares [ 14 14 , 50 50 ], local radial point interpolation [ 51 51 ], and local point interpolation [ 52 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 53 ]. The moving least squares approximation was created by mathematic researchers for data fitting and surface construction [ 54 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 42 ]: <math display="block" id="M1"><mrow><msup><mtext>u</mtext><mtext>h</mtext></msup><mo stretchy="false">(</mo><mtext>X</mtext><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><mrow><msub><mtext>p</mtext><mi>j</mi></msub><msub><mrow><mtext>(X)a</mtext></mrow><mi>j</mi></msub><msup><mrow><mtext>(X)=P</mtext></mrow><mtext>T</mtext></msup><mtext>(X)a(X)</mtext><mtext>.</mtext></mrow></mstyle></mrow></math> P T ( X ) is a complete monomial basis of order m . m is the number of monomial terms, here the linear basis m =3. <math display="block" id="M2"><mrow><msup><mtext>P</mtext><mtext>T</mtext></msup><msup><mrow><mtext>(X)=P</mtext></mrow><mtext>T</mtext></msup><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mn>1</mn><mo>,</mo><mi>x</mi><mo>,</mo><mi>y</mi></mrow><mo>}</mo></mrow></mrow></math> a ( X ) is a coefficient vector as follow: <math display="block" id="M3"><mrow><msup><mtext>a</mtext><mtext>T</mtext></msup><msup><mrow><mtext>(X)=a</mtext></mrow><mtext>T</mtext></msup><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><msub><mi>a</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>X</mi><mo stretchy="false">)</mo><mo>,</mo><msub><mi>a</mi><mn>1</mn></msub><mtext>(X)</mtext><mo>,</mo><mo>…</mo><mo>,</mo><msub><mi>a</mi><mi>m</mi></msub><mtext>(X)</mtext></mrow><mo>}</mo></mrow><mo>.</mo></mrow></math> The approximated values of the field function are: <math display="block" id="M4"><mrow><msup><mtext>u</mtext><mtext>h</mtext></msup><msub><mrow><mtext>(X,X</mtext></mrow><mi>I</mi></msub><mo stretchy="false">)</mo><mo>=</mo><msub><mrow><mtext>P(X</mtext></mrow><mi>I</mi></msub><mo stretchy="false">)</mo><mtext>a(X)</mtext><mi>I</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>..</mn><mo>,</mo><mi>n</mi><mo>.</mo></mrow></math> 2) In the second step, the weighted discrete L 2 norm is defined as follows: <math display="block" id="M5"><mrow><mtext>J(X)</mtext><mo>=</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mo stretchy="false">(</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mi>n</mi></msubsup><mrow><msub><mrow><mtext>W(X-X</mtext></mrow><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mstyle><msup><mrow><mo stretchy="false">[</mo><msup><mtext>U</mtext><mtext>h</mtext></msup><msub><mrow><mtext>(X,X</mtext></mrow><mi>i</mi></msub><mo stretchy="false">)</mo><mo>−</mo><msub><mtext>U</mtext><mi>i</mi></msub><mo stretchy="false">(</mo><msub><mtext>X</mtext><mi>i</mi></msub><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow><mn>2</mn></msup><mo>=</mo><msup><mrow><mo stretchy="false">[</mo><msup><mtext>P</mtext><mtext>T</mtext></msup><msub><mrow><mtext>(X</mtext></mrow><mi>i</mi></msub><mo stretchy="false">)</mo><mtext>a</mtext><mo stretchy="false">(</mo><mtext>X</mtext><mo stretchy="false">)</mo><mo>−</mo><msub><mtext>U</mtext><mi>i</mi></msub><mo stretchy="false">]</mo></mrow><mn>2</mn></msup><mo>.</mo></mrow></math> 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 <math id="M6"><mi>W</mi><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mtext>X</mtext></mstyle><mo>−</mo><msub><mrow><mstyle mathvariant="bold" mathsize="normal"><mtext>X</mtext></mstyle></mrow><mi>i</mi></msub><mo>)</mo></math> 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 42 ]: <math display="block" id="M7"><mrow><mo> </mo><mfrac><mrow><mo>∂</mo><mi>J</mi></mrow><mrow><mo>∂</mo><mstyle mathvariant="bold" mathsize="normal"><mi>a</mi></mstyle></mrow></mfrac><mo>=</mo><mn>0.</mn></mrow></math> This equation leads to a linear relationship (7) between a ( X ) and <math id="M8"><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle></math> : <math display="block" id="M9"><mrow><mo> </mo><msub><mrow><mtext>A(X)a(X)=B(X)U</mtext></mrow><mi>s</mi></msub><mtext>.</mtext></mrow></math> Solving Eq.(7) for a ( X ) gives [ 54 54 ]: <math display="block" id="M10"><mrow><msup><mrow><mtext>a(X)=A(X)</mtext></mrow><mrow><mtext>-1</mtext></mrow></msup><msub><mrow><mtext>B(X)U</mtext></mrow><mi>s</mi></msub><mtext>,</mtext></mrow></math> where U s is a vector for choosing the field function nodal parameters for interior nodes of the support domain as follows: <math display="block" id="M11"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mtext>s</mtext></msub><mo>=</mo><msup><mrow><mrow><mo>{</mo><mrow><msub><mtext>U</mtext><mn>1</mn></msub><msub><mrow><mtext>,U</mtext></mrow><mn>2</mn></msub><mo>,</mo><mo>…</mo><mi>，</mi><msub><mtext>U</mtext><mi>n</mi></msub></mrow><mo>}</mo></mrow></mrow><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle></msup></mrow></math> The matrices A ( X ) and B ( X ) are given as [ 54 54 , 55 55 ]: <math display="block" id="M12"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle></msup><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi><mi>P</mi></mstyle><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mi>i</mi></msub><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>−</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>i</mi></msub></mrow><mo>)</mo></mrow><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mrow><mo>(</mo><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>i</mi></msub></mrow><mo>)</mo></mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mstyle mathvariant="bold" mathsize="normal"><mi>T</mi></mstyle></msup><mrow><mo>(</mo><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow></math> <math display="block" id="M13"><mrow><mi>B</mi><mtext>(X)=</mtext><msup><mi>P</mi><mtext>T</mtext></msup><mi>W</mi><mtext>=[</mtext><msub><mi>W</mi><mtext>1</mtext></msub><mtext>(</mtext><mi>X</mi><mtext>)</mtext><mi>P</mi><mtext>(</mtext><msub><mi>X</mi><mtext>1</mtext></msub><mtext>),</mtext><mo>…</mo><mtext>,</mtext><msub><mi>W</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>X</mi><mo stretchy="false">)</mo><mi>P</mi><mo stretchy="false">(</mo><msub><mi>X</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></math> where <math id="M14"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mi>i</mi><mtext> </mtext><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>=</mo><mtext> </mtext><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mtext> </mtext><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>i</mi></mrow><mo>)</mo></mrow></mrow></math> . 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: http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10332438&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10332435&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10332438&type= 少公式(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: <math display="block" id="M15"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mtext>h</mtext></msup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><msub><mi>ϕ</mi><mi>i</mi></msub><mrow><mo>(</mo><mtext>X</mtext><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>Φ</mi></mstyle><mtext>T</mtext></msup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mtext>s</mtext></msub><mo>.</mo></mrow></math> <math id="M16"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>Φ</mi></mstyle><mtext>T</mtext></msup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow></mrow></math> can be rewritten as below [ 55 55 ]: <math display="block" id="M17"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>Φ</mi></mstyle><mtext>T</mtext></msup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>=</mo><msub><mrow><mrow><mo>{</mo><mrow><msub><mi>ϕ</mi><mn>1</mn></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><msub><mi>ϕ</mi><mn>2</mn></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo> </mo><mo> </mo><mo>…</mo><mo>…</mo><mo>.</mo><msub><mi>ϕ</mi><mi>n</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow></mrow><mo>}</mo></mrow></mrow><mrow><mn>1</mn><mo>*</mo><mi>n</mi></mrow></msub><mo>=</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>p</mi></mstyle><mtext>T</mtext></msup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mrow><mo>−</mo><mn>1</mn></mrow></msup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M18"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>Φ</mi></mstyle><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow></mrow></math> indicates the MLS shape function vector for n nodes in the support domain of point X . Therefore, the shape function for the i th node is introduced in Eq. (15)[ 55 55 ]: <math display="block" id="M19"><mrow><msub><mi>ϕ</mi><mi>i</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>p</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><msub><mrow><mrow><mo>(</mo><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mrow><mo>−</mo><mn>1</mn></mrow></msup><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow><mrow><mi>j</mi><mi>i</mi></mrow></msub><mo>.</mo></mrow></math> A flowchart of the algorithm for the above steps is shown in Fig. 3 Fig. 3 . 2.2 Solving the neutron diffusion equation by MLPG 2 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 <math id="M20"><mrow><msub><mi>ϕ</mi><mi>j</mi></msub></mrow></math> is a function of X only [ 56 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: <math display="block" id="M21"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>h</mi></msup><mrow><mo>(</mo><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>ϕ</mi><mi>j</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>,</mo></mrow></math> 2.2.1 The temporal derivation of <math id="M22"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>h</mi></msup></mrow></math> is: <math display="block" id="M23"><mrow><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>t</mi><mi>h</mi></msubsup><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>ϕ</mi><mi>j</mi></msub><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>.</mo></mrow></math> 2.2.2 The spatial derivation of <math id="M24"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>h</mi></msup></mrow></math> is: <math display="block" id="M25"><mrow><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>x</mi><mi>h</mi></msubsup><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mover><mrow><msub><mi>ϕ</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>.</mo></mrow></math> 2.2.3 In addition, the second-order spatial domain is: <math display="block" id="M26"><mrow><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mi>x</mi><mi>x</mi></mrow><mi>h</mi></msubsup><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mover><mi>ϕ</mi><mo>¨</mo></mover><mi>j</mi><mrow><mo>(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo>)</mo></mrow><mo>.</mo></mrow></math> The derivatives are denoted by subscripts. In addition, a local quadrature domain, <math id="M27"><mrow><msub><mtext>Ω</mtext><mrow><msub><mi>Q</mi><mi>i</mi></msub></mrow></msub></mrow></math> is defined for each interior node, x i . R indicates the residual function. The integral of the weighted residual is zero. [ 57 57 ]. <math display="block" id="M28"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><msub><mi>Ω</mi><mrow><msub><mi>Q</mi><mi>i</mi></msub></mrow></msub></mrow><mrow/></munderover><mi>R</mi><msub><mi>W</mi><mi>i</mi></msub><mtext>d</mtext><msub><mi>Ω</mi><mrow><msub><mi>Q</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><mn>0</mn></mrow></math> As mentioned earlier, <math id="M29"><mrow><msub><mi>W</mi><mi>i</mi></msub></mrow></math> is a weight function that is non-zero only on <math id="M30"><mrow><msub><mi>Ω</mi><mrow><mi>Q</mi><mi>i</mi></mrow></msub></mrow></math> and similar to the weight function used in the MLSA process. 2.3 Determining the residual of the neutron diffusion equation 2 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): <math display="block" id="M31"><mrow><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mi>v</mi><mrow/></munderover><mfrac><mn>1</mn><mi>v</mi></mfrac><mi>φ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mtext>d</mtext><mi>v</mi></mrow></math> The production rate is considered as a neutron source in Eq. (22): <math display="block" id="M32"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mi>v</mi><mrow/></munderover><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mtext>d</mtext><mi>v</mi></mrow></math> The neutron absorption term is presented in Eq. (23): <math display="block" id="M33"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mi>v</mi><mrow/></munderover><msub><mtext>Σ</mtext><mi>a</mi></msub><mi>φ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mtext>d</mtext><mi>v</mi><mo>.</mo></mrow></math> The leakage term is shown in Eq. (24): <math display="block" id="M34"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mi>v</mi><mrow/></munderover><mo>∇</mo><mo>.</mo><mi>J</mi><mtext>d</mtext><mi>v</mi><mo>=</mo><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mi>v</mi><mrow/></munderover><mo>∇</mo><mo>.</mo><mrow><mo>[</mo><mrow><mi>D</mi><mo>∇</mo><mi>φ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mtext>d</mtext><mi>v</mi><mo>.</mo></mrow></math> The mathematical formulation of the neutron diffusion theory is obtained as follows: <math display="block" id="M35"><mrow><mfrac><mrow><mo>∂</mo><mi>φ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>−</mo><mi>D</mi><msup><mo>∇</mo><mn>2</mn></msup><mi>φ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mtext>Σ</mtext><mi>a</mi></msub><mo> </mo><mi>φ</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mtext> </mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> Equation (25) can be rewritten as equation (26): <math display="block" id="M36"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>t</mi></msub><mo>−</mo><mi>D</mi><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mi>x</mi><mi>x</mi></mrow></msub><mo>+</mo><msub><mtext>Σ</mtext><mi>a</mi></msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mo>=</mo><mn>0.</mn></mrow></math> Where u is the solution, D is the diffusion coefficient, s is a neutron source, and <math id="M37"><mrow><mo> </mo><msub><mtext>Σ</mtext><mi>a</mi></msub></mrow></math> is the absorbtion cross section. The residual for this equation is defined according to an approximation solution as follows: <math display="block" id="M38"><mrow><mtext>Residual</mtext><mo>=</mo><mi>R</mi><mo>=</mo><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>t</mi><mi>h</mi></msubsup><mo>−</mo><mi>D</mi><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mi>x</mi><mi>x</mi></mrow><mi>h</mi></msubsup><mo>+</mo><msub><mtext>Σ</mtext><mi>a</mi></msub><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>h</mi></msup><mo>−</mo><mi>s</mi><mo>=</mo><mn>0.</mn></mrow></math> Equation (27) is substituted into equation (20) to give a system of n linear equations for n nodes: <math display="block" id="M39"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><mtext>Ω</mtext><msub><mi>q</mi><mi>i</mi></msub></mrow><mrow/></munderover><mrow><mo>[</mo><mrow><msubsup><mi>u</mi><mi>t</mi><mi>h</mi></msubsup><mo>−</mo><mi>D</mi><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mi>x</mi><mi>x</mi></mrow><mi>h</mi></msubsup><mo>+</mo><msub><mtext>Σ</mtext><mi>a</mi></msub><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>h</mi></msup><mo>−</mo><mi>s</mi></mrow><mo>]</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mtext>d</mtext><msub><mi>Ω</mi><mrow><mi>Q</mi><msub><mi>q</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><mn>0</mn></mrow></math> 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): <math display="block" id="M40"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><mtext>Ω</mtext><msub><mi>q</mi><mi>i</mi></msub></mrow><mrow/></munderover><mrow><mo>[</mo><mrow><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>t</mi><mi>h</mi></msubsup></mrow><mo>]</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mtext>d</mtext><msub><mi>Ω</mi><mrow><mi>Q</mi><msub><mi>q</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><mo>∫</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>.</mo></mrow></math> b) The second term is obtained by substituting Eq. (19) into Eq. (28): <math display="block" id="M41"><mtable columnalign="left"><mtr><mtd><mo>−</mo><mi>D</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><mtext>Ω</mtext><msub><mi>q</mi><mi>i</mi></msub></mrow><mrow/></munderover><mrow><mo>[</mo><mrow><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mi>x</mi><mi>x</mi></mrow><mi>h</mi></msubsup></mrow><mo>]</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mtext>d</mtext><msub><mi>Ω</mi><mrow><msub><mi>q</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><mtext>Ω</mtext><msub><mi>q</mi><mi>i</mi></msub></mrow><mrow/></munderover><mi>D</mi><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>x</mi><mi>h</mi></msubsup><msub><mi>w</mi><mrow><mi>i</mi><mo>,</mo><mi>X</mi></mrow></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><msub><mi>Ω</mi><mrow><msub><mi>q</mi><mi>i</mi></msub></mrow></msub><mo>−</mo><mi>D</mi><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>x</mi><mi>h</mi></msubsup><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>|</mo><mrow><mtable><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mi>n</mi></mrow></mtd></mtr><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mn>1</mn></mrow></mtd></mtr></mtable></mrow></mrow></mtd></mtr><mtr><mtd><mo>=</mo><mi>D</mi><mo>∫</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mover><mi>φ</mi><mo>´</mo></mover><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mover><mi>w</mi><mo>´</mo></mover><mi>i</mi></msub><mrow><mo>(</mo><mi>X</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>−</mo><mi>D</mi><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> c) The third term is obtained by substituting Eq. (16) into Eq. (28): <math display="block" id="M42"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><mtext>Ω</mtext><msub><mi>q</mi><mi>i</mi></msub></mrow><mrow/></munderover><msub><mtext>Σ</mtext><mi>a</mi></msub><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>h</mi></msup><msub><mi>w</mi><mi>i</mi></msub><mtext>d</mtext><msub><mi>Ω</mi><mrow><msub><mi>q</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><mo>∫</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mtext>Σ</mtext><mi>a</mi></msub><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi><mi>h</mi></msubsup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>X</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>.</mo></mrow></math> d) The fourth term is obtained by substituting the fission effective factors into Eq. (22): <math display="block" id="M43"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><mtext>Ω</mtext><msub><mi>q</mi><mi>i</mi></msub></mrow><mrow/></munderover><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><msub><mi>w</mi><mi>i</mi></msub><mtext>d</mtext><msub><mi>Ω</mi><mrow><msub><mi>q</mi><mi>i</mi></msub></mrow></msub><mo>=</mo><mo>−</mo><mo>∫</mo><mfrac><mi>χ</mi><mi>λ</mi></mfrac><mi>υ</mi><msub><mi>Σ</mi><mi>f</mi></msub><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>.</mo></mrow></math> Therefore, Eq. (28) can be rewritten as: <math display="block" id="M44"><mrow><mo>∫</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>+</mo><mi>D</mi><mo>∫</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mrow><mover><mi>φ</mi><mo>´</mo></mover></mrow><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mrow><mover><mi>w</mi><mo>´</mo></mover></mrow><mi>i</mi></msub><mrow><mo>(</mo><mi>X</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>−</mo><mtext>D</mtext><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>|</mo><mrow><mtable><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mi>n</mi></mrow></mtd></mtr><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mn>1</mn></mrow></mtd></mtr></mtable></mrow></mrow><mo>+</mo><mo>∫</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mtext>Σ</mtext><mi>a</mi></msub><msubsup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi><mi>h</mi></msubsup><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>X</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>dx</mtext><mo>−</mo><mo>∫</mo><mfrac><mi>χ</mi><mi>λ</mi></mfrac><msub><mrow><mtext>υΣ</mtext></mrow><mi>f</mi></msub><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>dx</mtext><mo>=</mo><mn>0.</mn></mrow></math> For clarity, the sums are moved outside the integrals as shown in equation (34): <math display="block" id="M45"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mo>∫</mo><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>+</mo><mi>D</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mo>∫</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mrow><mover><mi>φ</mi><mo>´</mo></mover></mrow><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mrow><mover><mi>w</mi><mo>´</mo></mover></mrow><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>−</mo><mi>D</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>|</mo><mrow><mtable><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mi>n</mi></mrow></mtd></mtr><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mn>1</mn></mrow></mtd></mtr></mtable></mrow></mrow><mo>+</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mo>∫</mo><msub><mtext>Σ</mtext><mi>a</mi></msub><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>−</mo><mo>∫</mo><mfrac><mi>χ</mi><mi>λ</mi></mfrac><msub><mrow><mtext>υΣ</mtext></mrow><mi>f</mi></msub><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>dx</mtext><mo>=</mo><mn>0.</mn></mrow></math> This equation can be rearranged into its final form as shown in Eq.(35): <math display="block" id="M46"><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>∫</mo><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>+</mo><mi>D</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>∫</mo><msub><mrow><mover><mi>φ</mi><mo>´</mo></mover></mrow><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mrow><mover><mi>w</mi><mo>´</mo></mover></mrow><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>−</mo><mi>D</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>|</mo><mrow><mtable><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mi>n</mi></mrow></mtd></mtr><mtr><mtd><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mi>q</mi><mi>i</mi><mn>1</mn></mrow></mtd></mtr></mtable></mrow></mrow><mo>+</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>∫</mo><msub><mtext>Σ</mtext><mi>a</mi></msub><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>d</mtext><mi>x</mi><mo>−</mo><mfrac><mi>χ</mi><mi>λ</mi></mfrac><msub><mrow><mtext>υΣ</mtext></mrow><mi>f</mi></msub><mo>∫</mo><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mtext>dx</mtext><mo>=</mo><mn>0</mn></mrow></math> 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 33 , 45 45 , 58 58 ]. 2.4 Solution of the time-dependent neutron diffusion equation 2 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 <math id="M47"><mrow><msub><mrow><mover><mi>u</mi><mo>´</mo></mover></mrow><mi>j</mi></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi><mo>/</mo><mn>2</mn></mrow><mo>)</mo></mrow></mrow></math> with the central difference approximation [ 59 59 ]: <math display="block" id="M48"><mrow><msub><mrow><mover><mi>u</mi><mo>´</mo></mover></mrow><mi>j</mi></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi><mo>/</mo><mn>2</mn></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mi>u</mi><mi>j</mi></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>−</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mo>.</mo></mrow></math> In this discretization method, <math id="M49"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi><mo>/</mo><mn>2</mn></mrow><mo>)</mo></mrow></mrow></math> is approximated as the average of <math id="M50"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></math> and <math id="M51"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> . 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 <math id="M52"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><msup><mrow/><mrow><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></msup></mrow></math> = ( <math id="M53"><mrow><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mn>1</mn></msub><msup><mrow/><mrow><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></msup><mo>…</mo><mn>..</mn><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>n</mi></msub><msup><mrow/><mrow><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></msup></mrow></math> ), the matrix form of Eq. (35) is derived as below: <math display="block" id="M54"><mrow><mi>M</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><mi>N</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> <math display="block" id="M55"><mrow><mi>M</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mfrac><mrow><mtext>Δ</mtext><mi>t</mi></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mi>N</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mi>j</mi></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mfrac><mrow><mtext>Δ</mtext><mi>t</mi></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> <math display="block" id="M56"><mrow><mi>M</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mrow><mtext>t</mtext><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mi>N</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mrow><mtext>t</mtext><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>+</mo><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>S</mi></mstyle><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> <math display="block" id="M57"><mrow><mn>2</mn><mi>M</mi><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>−</mo><mn>2</mn><mi>M</mi><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>+</mo><mi>N</mi><mtext>Δ</mtext><mi>t</mi><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>+</mo><mi>N</mi><mtext>Δ</mtext><mi>t</mi><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>−</mo><mn>2</mn><mi>t</mi><mtext>Δ</mtext><mi>t</mi><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> <math display="block" id="M58"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mrow><mn>2</mn><mi>M</mi><mo>+</mo><mi>N</mi><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mn>2</mn><mi>t</mi><mtext>Δ</mtext><mi>t</mi><mo>+</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></msup><mrow><mo>(</mo><mrow><mn>2</mn><mi>M</mi><mo>−</mo><mi>N</mi><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where M and N are coefficients of <math id="M59"><mrow><mover><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow><mo>´</mo></mover></mrow></math> and <math id="M60"><mrow><msub><mi>u</mi><mi>j</mi></msub></mrow></math> , respectively. The proposed final relationship, which was used in the computational algorithm as the main formula at each time-step, is: <math display="block" id="M61"><mrow><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mtext>Δ</mtext><mi>t</mi></mrow><mo>)</mo></mrow></mrow></msup><mstyle mathvariant="bold" mathsize="normal"><mi>L</mi></mstyle><mo>=</mo><msup><mi>s</mi><mo>*</mo></msup><mo>+</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>u</mi></mstyle><mrow><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></msup><mstyle mathvariant="bold" mathsize="normal"><mi>K</mi></mstyle><mo>.</mo></mrow></math> 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 <math id="M62"><mrow><msub><mi>ϕ</mi><mi>j</mi></msub><mrow><mo>(</mo><mrow><msub><mi>x</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow></math> . If x j is in the support domain of any point in the quadrature domain of the internal node x i , L ij is <math id="M63"><mrow><mn>2</mn><msub><mi>M</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>+</mo><msub><mi>N</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mtext>Δ</mtext><mi>t</mi></mrow></math> . 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. <math id="M64"><mo> </mo></math> If x j is in the support domain of any point in the quadrature domain of the internal node x i , L ij is <math id="M65"><mrow><mtext> </mtext><mn>2</mn><msub><mi>M</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>−</mo><msub><mi>N</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mtext>Δ</mtext><mi>t</mi></mrow></math> . Otherwise, the value of L ij is zero and <math id="M66"><mrow><msup><mi>S</mi><mtext>*</mtext></msup></mrow></math> 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 48 ]. These conditions consist of having an arbitrary nodal distribution, being numerically stable, satisfying a certain order of consistency, and being compactly supported [ 48 48 ]. The accuracy of interpolation for the arbitrary points depends on the number of nodes in the support domain [ 53 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 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 <math display="block" id="M67"><mtable columnalign="left"><mtr><mtd><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><mo>−</mo><mn>4</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>2</mn></msup><mo>+</mo><mn>4</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>3</mn></msup><mo>,</mo><mtext>               for  </mtext><mn>0</mn><mo>≤</mo><msub><mi>r</mi><mi>i</mi></msub><mo>≤</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><mo>−</mo><mn>4</mn><msub><mi>r</mi><mi>i</mi></msub><mo>+</mo><mn>4</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>2</mn></msup><mo>−</mo><mfrac><mn>4</mn><mn>3</mn></mfrac><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>3</mn></msup><mo>,</mo><mtext>                 for</mtext><mo> </mo><mo> </mo><mo> </mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mo>≤</mo><msub><mi>r</mi><mi>i</mi></msub><mo>≤</mo><mn>1</mn></mtd></mtr><mtr><mtd><mo>=</mo><mn>0</mn><mo>,</mo><mtext>                                                     for </mtext><msub><mi>r</mi><mi>i</mi></msub><mo>></mo><mn>1</mn></mtd></mtr></mtable></math> 2. Quadratic spline weight function <math display="block" id="M68"><mtable columnalign="left"><mtr><mtd><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mn>1</mn><mo>−</mo><mn>6</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>2</mn></msup><mo>+</mo><mn>8</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>3</mn></msup><mo>−</mo><mn>3</mn><msubsup><mi>r</mi><mi>i</mi><mn>4</mn></msubsup><mtext>     for   </mtext><mn>0</mn><mo>≤</mo><msub><mi>r</mi><mi>i</mi></msub><mo>≤</mo><mn>1</mn></mtd></mtr><mtr><mtd><mo>=</mo><mn>0</mn><mo>,</mo><mtext>                                        for   </mtext><msub><mi>r</mi><mi>i</mi></msub><mo>></mo><mn>1</mn></mtd></mtr></mtable></math> 3. Compact support radial basis function (CSRBF) <math display="block" id="M69"><mtable columnalign="left"><mtr><mtd><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><msup><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>r</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mn>6</mn></msup><mo>+</mo><mrow><mo>[</mo><mrow><mn>6</mn><mo>+</mo><mn>36</mn><msub><mi>r</mi><mi>i</mi></msub><mo>+</mo><mn>82</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>2</mn></msup><mo>+</mo><mn>72</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>3</mn></msup><mo>+</mo><mn>30</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>4</mn></msup><mo>+</mo><mn>5</mn><msub><mi>r</mi><mi>i</mi></msub><msup><mrow/><mn>5</mn></msup></mrow><mo>]</mo></mrow><mo>,</mo><mtext>    for     </mtext><mn>0</mn><mo>≤</mo><msub><mi>r</mi><mi>i</mi></msub><mo>≤</mo><mn>1</mn></mtd></mtr><mtr><mtd><mo>=</mo><mn>0</mn><mo>,</mo><mtext>                                                  for    </mtext><msub><mi>r</mi><mi>i</mi></msub><mo>></mo><mn>1</mn></mtd></mtr></mtable></math> 4. Gaussian weight function <math display="block" id="M70"><mrow><msub><mi>w</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mrow><mi>exp</mi><mrow><mo>[</mo><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>d</mi><mi>i</mi></msub></mrow><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow><mo>−</mo><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>r</mi><mi>i</mi></msub></mrow><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow></mrow><mrow><mn>1</mn><mo>−</mo><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>r</mi><mi>i</mi></msub></mrow><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow></mrow></mfrac><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mi>f</mi><mi>o</mi><mi>r</mi><mo> </mo><mo> </mo><mn>0</mn><mo>≤</mo><msub><mi>d</mi><mi>i</mi></msub><mo>≤</mo><msub><mi>r</mi><mi>i</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mn>0</mn><mtext>                             </mtext><mi>f</mi><mi>o</mi><mi>r</mi><mtext>   </mtext><msub><mi>d</mi><mi>i</mi></msub><mo>≥</mo><msub><mi>r</mi><mi>i</mi></msub></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> Where <math id="M71"><mrow><msub><mi>r</mi><mi>i</mi></msub></mrow></math> 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, <math id="M72"><mrow><msub><mi>c</mi><mi>i</mi></msub><mtext> </mtext></mrow></math> 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, <math id="M73"><mrow><msub><mi>c</mi><mi>i</mi></msub><mtext> </mtext></mrow></math> 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, <math id="M74"><mrow><msub><mi>Ω</mi><mrow><msub><mi>Q</mi><mi>i</mi></msub></mrow></msub></mrow></math> , is defined [ 60 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 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 . 3 Results and Discussion 1 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 2 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 Fig. 5 . The logarithm of the error versus the number of nodes is presented in Fig. 6 Fig. 6 . It is clear that increasing the number of nodes decreases the error. Fig. 7 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 Fig. 7b , which demonstrates that the behavior of the proposed method is stable for less than 500 nodes. 3.2 Constant source problem 2 In this test problem, a distributed source problem was considered. This problem consists of two models, a single region ( Fig. 8 Fig. 8 ) and two regions ( Fig. 11 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: CASE-I: The geometry of the single region is presented in Fig. 8 Fig. 8 . The parameters required for this calculation are shown in Table 1 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 Figs. 9 and 10 10 , respectively. 1 1 Applied models 1 1 Energy group 1 1 Absorption cross section 1 1 Nu-Fission cross section 1 1 Diffusion coefficient 1 1 Slab I 1 1 1 1 1 0.03051 1 1 0.000256 1 1 1.440130 1 1   1 1 2 1 1 0.10582 1 1 0.154726 1 1 0.015723 CASE-II: The geometry of the two-region problem is presented in Fig. 11 Fig. 11 and the required parameters are shown in Table 2 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 Figs. 12 and 13 13 , respectively. As shown, the results of the proposed method exhibit good agreement with those of the analytical solution. 1 1 Applied models 1 1 Energy group 1 1 Absorption cross section 1 1 Nu-Fission cross section 1 1 Diffusion coefficient 1 1 Region 1 1 Slab-II 1 1 1 1 1 0.03051 1 1 0.000256 1 1 1.440130 1 1 Fuel 1 1   1 1 2 1 1 0.10582 1 1 0.154726 1 1 0.015723 1 1   1 1 1 1 1 0.035128 1 1 0.000000 1 1 1.852524 1 1 Reflector 1 1   1 1 2 1 1 0.032540 1 1 0.000000 1 1 0.294517 1 1   The average error of these results obtained by applying various weight functions is presented in Table 3 Table 3 , which indicates that the accuracy of the solution is improved by using Gaussian weight functions. In Table 4 Table 4 , the average error and CPU time of this problem are demonstrated for the Gaussian and CSRBF weight functions. Fig. 14 Fig. 14 presents a comparison of the fast neutron flux in CASE-I and CASE-II obtained by inserting 50×50 nodes. 1 1 Error 4 1 MLPG Applied Weight Functions 1 1   1 1 Gaussian 1 1 CSRBF 1 1 Cubic Spline 1 1 Quadratic Spline 1 1 CASE-I 1 1 0.0035 1 1 0.073 1 1 0.094 1 1 0.085 1 1 CASE-II 1 1 0.0038 1 1 0.070 1 1 0.091 1 1 0.083 1 1 Problem No. 4 1 MLPG Applied Weight Functions 1 1   2 1 Gaussian 2 1 CSRBF 1 1   1 1 Error 1 1 CPU Time (s) 1 1 Error 1 1 CPU Time(s) 1 1 CASE-I 1 1 0.0035 1 1 2 1 1 0.073 1 1 2.5 1 1 CASE-II 1 1 0.0038 1 1 3.5 1 1 0.070 1 1 6 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 ). Fig. 15 Fig. 15 indicates the performance of the MLPG method was improved by increasing the number of nodes. Fig. 16 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. 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 Fig. 17 . A comparison of the contours of the thermal neutron flux of this problem is shown in Fig. 18 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. 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 4 . Table 5 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 Fig. 19 . 1 1 Distribution Pattern 4 1 MLPG Applied Weight Functions 1 1   1 1 Gaussian 1 1 CSRBF 1 1 Cubic Spline 1 1 Quadratic Spline 1 1 Normal Distribution 1 1 0.0035 1 1 0.073 1 1 0.094 1 1 0.085 1 1 Scatter Distribution 1 1 0.0030 1 1 0.059 1 1 0.082 1 1 0.065 3.3 Void test problem 2 The void test problem was first studied by Abuzaid and Gashut [ 61 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 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 Fig. 21 , which indicates a good agreement. 3.4 Neutronic calculations for a SMR reactor 2 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 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 62 ]. 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 3 . 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 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 Table 6 demonstrates the effective multiplication factor as a main parameter of the core calculation. Tables 7 Tables 7 and 8 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 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. 1 1 Parameter 1 1 Citation 1 1 Referenced value 4 1 MLPG Applied weight functions 1 1 Nodes Insertion 1 1   1 1   1 1   1 1 Gaussian 1 1 CSRBF 1 1 Cubic Spline 1 1 Quadratic Spline 1 1   1 1 Effective Multiplication Factor 1 1 1.12414 1 1 1.12511 1 1 1.12332 1 1 1.11032 1 1 1.11054 1 1 1.11902 1 1 37×37 1 1   1 1   1 1   1 1 1.12415 1 1 1.11940 1 1 1.11891 1 1 1.11974 1 1 91×91 1 1   1 1   1 1   1 1 1.12496 1 1 1.12439 1 1 1.12475 1 1 1.11996 1 1 289×289 1 1 Parameter 1 1 Citation 1 1 Referenced value 4 1 MLPG Applied weight functions 1 1 Nodes Insertion 1 1   1 1   1 1   1 1 Gaussian 1 1 CSRBF 1 1 Cubic Spline 1 1 Quadratic Spline 1 1   1 1 Maximum to Average Flux Ratio 1 1 2.48520 1 1 2.4980 1 1 2.48213 1 1 2.43191 1 1 2.44521 1 1 2.39120 1 1 37×37 1 1   1 1   1 1   1 1 2.48500 1 1 2.45225 1 1 2.46378 1 1 2.44032 1 1 91×91 1 1   1 1   1 1   1 1 2.49739 1 1 2.48821 1 1 2.49707 1 1 2.45318 1 1 289×289 1 1 Parameter 1 1 Citation 1 1 Referenced value 4 1 MLPG Applied weight functions 1 1 Nodes Insertion 1 1   1 1   1 1   1 1 Gaussian 1 1 CSRBF 1 1 Cubic spline 1 1 Quartic spline 1 1   1 1 Control Rod Worth 1 1 0.15613 1 1 0.1571 1 1 0.15243 1 1 0.14254 1 1 0.1550 1 1 0.14382 1 1 37×37 1 1   1 1   1 1   1 1 0.15615 1 1 0.14367 1 1 0.1561 1 1 0.14630 1 1 91×91 1 1   1 1   1 1   1 1 0.15694 1 1 0.14975 1 1 0.15679 1 1 0.14587 1 1 289×289 1 1 Error 4 1 MLPG Applied weight functions 1 1   1 1 Gaussian 1 1 CSRBF 1 1 Cubic spline 1 1 Quartic spline 1 1 Effective Multiplication Factor 1 1 0.001 1 1 0.0064 1 1 0.032 1 1 0.045 1 1 Maximum to Average Flux Ratio 1 1 0.002 1 1 0.041 1 1 0.037 1 1 0.179 1 1 Control Rod Worth 1 1 0.001 1 1 0.046 1 1 0.020 1 1 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 Table 10 . The errors decreases with increasing node number. In addition, the scatter distribution results are more accurate than those of the normal distribution. 1 1 Parameter 6 1 MLPG Gaussian weight function 1 1   3 1 Nodes Normal Distribution 3 1 Nodes Scatter Distribution 1 1   1 1 Error 1 1 CPU Time(s) 1 1 Nodes 1 1 Error 1 1 CPU Time(s) 1 1 Nodes 1 1 Effective Multiplication Factor 1 1 0.012 1 1 4.0 1 1 37×37 1 1 0.010 1 1 3.5 1 1 37×37 1 1   1 1 0.087 1 1 5.5 1 1 91×91 1 1 0.063 1 1 4.0 1 1 91×91 1 1   1 1 0.001 1 1 10 1 1 289×289 1 1 0.001 1 1 7.5 1 1 289×289 The fast flux and thermal flux of the reactor core is demonstrated in Fig. 23 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. 24 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 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. 4 Conclusion 1 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 ) Computational multi scale modeling of fluids and solids G. R. Liu , Mesh free methods: moving beyond the finite element method . ( Taylor & Francis , Boca Raton , 2009 ) Mesh free methods: moving beyond the finite element method 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 Mesh free methods S. Li , W. K. Liu , Mesh free and particle methods and their applications . Appl. Mech. Rev . 55 , 1 - 34 ( 2002 ). doi: 10.1115/1.1431547 http://doi.org/10.1115/1.1431547 Mesh free and particle methods and their applications Y. T. Gu , Mesh free methods and their comparisons . Int. J. Comput. Methods 2 , 477 - 515 ( 2005 ). doi: 10.1142/S0219876205000673 http://doi.org/10.1142/S0219876205000673 Mesh free methods and their comparisons A. Tayebi , Y. Shekari , M. H. Heydari , A meshless method for solving two-dimensional variable-order time fractional advection-diffusion equation . J. Comput. Phys . 340 , 655 - 669 ( 2017 ). doi: 10.1016/j.jcp.2017.03.061 http://doi.org/10.1016/j.jcp.2017.03.061 A meshless method for solving two-dimensional variable-order time fractional advection-diffusion equation Q. H. Li , S. S. Chen , G. X. Kou , Transient heat conduction analysis using the MLPG method and modified precise time step integration method . J. Comput. Phys . 230 , 2736 - 2750 ( 2011 ). doi: 10.1016/j.jcp.2011.01.019 http://doi.org/10.1016/j.jcp.2011.01.019 Transient heat conduction analysis using the MLPG method and modified precise time step integration method X. H. Wu , W. Q. Tao , S. P. Shen et al ., A stabilized MLPG method for steady state incompressible fluid flow simulation . J. Comput. Phys . 229 , 8564 - 8577 ( 2010 ). doi: 10.1016/j.jcp.2010.08.001 http://doi.org/10.1016/j.jcp.2010.08.001 A stabilized MLPG method for steady state incompressible fluid flow simulation R. Salehi , M. Dehghan , A moving least square reproducing polynomial meshless method . Appl. Numer. Math . 69 , 34 - 58 ( 2013 ). doi: 10.1016/j.apnum.2013.03.001 http://doi.org/10.1016/j.apnum.2013.03.001 A moving least square reproducing polynomial meshless method R. A. Gingold , J. J. Monaghan , Smoothed particle hydrodynamics: theory and application to non-spherical stars . Mon. Not. R. Astron. Soc . 181 , 375 - 389 ( 1977 ). doi: 10.1093/mnras/181.3.375 http://doi.org/10.1093/mnras/181.3.375 Smoothed particle hydrodynamics: theory and application to non-spherical stars P. W. Cleary , Modeling confined multi-material heat and mass flows using SPH . Appl. Math. Model . 22 , 981 - 993 ( 1998 ). doi: 10.1016/S0307-904X(98)10031-8 http://doi.org/10.1016/S0307-904X(98)10031-8 Modeling confined multi-material heat and mass flows using SPH G. R. Liu , K. Y. Dai , K. M. Lim , A point interpolation mesh free method for static and frequency analysis of two-dimensional piezoelectric structures . Comput. Mech . 29 , 510 - 519 ( 2002 ). doi: 10.1007/s00466-002-0360-9 http://doi.org/10.1007/s00466-002-0360-9 A point interpolation mesh free method for static and frequency analysis of two-dimensional piezoelectric structures B. Nayroles , G. Touzot , P. Villon , Generalizing the finite element method: Diffuse approximation and diffuse elements . Comput. Mech . 10 , 307 - 318 ( 1992 ). doi: 10.1007/BF00364252 http://doi.org/10.1007/BF00364252 Generalizing the finite element method: Diffuse approximation and diffuse elements W. K. Liu , S. Li , T. Belytschko , Moving least-square reproducing kernel methods: part I: methodology and convergence . Comput. Method, Appl. Mech . 143 , 113 - 154 ( 1977 ). doi: 10.1016/S0045-7825(96)01132-2 http://doi.org/10.1016/S0045-7825(96)01132-2 Moving least-square reproducing kernel methods: part I: methodology and convergence W. K. Liu , J. Adee , S. Jun et al ., Reproducing kernel particle methods for elastic and plastic problems . in Advanced Computational Methods for Material Modeling , edited by D. A. Siginer, W. E. VanArsdale, C. M. Altan, and A. N. Alexandrou (ASME, New Orleans, 1993), pp. 175-189 Reproducing kernel particle methods for elastic and plastic problems W. K. Liu , S. Jun , Y. F. Zhang , Reproducing kernel particle methods . Numer. Method Fluids 20 , 1081 - 1106 ( 1995 ). doi: 10.1002/fld.1650200824 http://doi.org/10.1002/fld.1650200824 Reproducing kernel particle methods E. Oñate , S. Idelsohn , O. C. Zienkiewicz et al ., A finite point method in computational mechanics - Applications to convective transport and fluid flow . Int. J. Numer. Methods Eng . 39 , 3839 - 3866 ( 1996 ). doi: 10.1002/(SICI)1097-0207(19961130)39:22<3839::AID-NME27>3.0.CO;2-R http://doi.org/10.1002/(SICI)1097-0207(19961130)39:22<3839::AID-NME27>3.0.CO;2-R A finite point method in computational mechanics - Applications to convective transport and fluid flow Y. T. Gu , G. R. Liu , Using radial function basis in a boundrary-type meshless method, boundrary point method . in Int. Conf. Sci. Enrgy Comput . ( Beijing , 2001 ) Using radial function basis in a boundrary-type meshless method, boundrary point method M. D. Buhmann , Radial basis functions: theory and implementions . ( cambridge university , Cambridge , 2003 ) Radial basis functions: theory and implementions H. Mehrabi , B. Voosoghi , On estimating the curvature attributes and strain invariants of deformed surface through radial basis functions . Comput. Appl. Math . 1 - 18 ( 2016 ). doi: 10.1007/s40314-016-0380-2 http://doi.org/10.1007/s40314-016-0380-2 On estimating the curvature attributes and strain invariants of deformed surface through radial basis functions Q. Li , J. Soric , T. Jarak et al ., A locking-free meshless local Petrov-Galerkin formulation for thick and thin plates . J. Comput. Phys . 208 , 116 - 133 ( 2005 ). doi: 10.1016/j.jcp.2005.02.008 http://doi.org/10.1016/j.jcp.2005.02.008 A locking-free meshless local Petrov-Galerkin formulation for thick and thin plates P. Mycek , G. Pinon , G. Germain , Formulation and analysis of a diffusion-velocity particle model for transport-dispersion equations . 35 , 447 - 473 ( 2014 ). doi: 10.1007/s40314-014-0200-5 http://doi.org/10.1007/s40314-014-0200-5 Formulation and analysis of a diffusion-velocity particle model for transport-dispersion equations S. Abbasbandy , A. Shirzadi , MLPG method for two-dimensional diffusion equation with Neumann’s and non-classical boundary conditions . Appl. Numer. Math . 61 , 170 - 180 ( 2011 ). doi: 10.1016/j.apnum.2010.09.002 http://doi.org/10.1016/j.apnum.2010.09.002 MLPG method for two-dimensional diffusion equation with Neumann’s and non-classical boundary conditions B. Rokrok , H. Minuchehr , A. Zolfaghari , Application of radial point interpolation method to neutron diffusion field . Trend Appl. Sci. Res . 7 , 18 - 31 ( 2012 ). doi: 10.3923/tasr.2012.18.31 http://doi.org/10.3923/tasr.2012.18.31 Application of radial point interpolation method to neutron diffusion field T. Tanbay , B. Ozgener , A comparison of the meshless RBF collocation method with finite element and boundary element methods in neutron diffusion calculations . Eng. Anal. Bound. Elem . 46 , 30 - 40 ( 2014 ). doi: 10.1016/j.enganabound.2014.05.005 http://doi.org/10.1016/j.enganabound.2014.05.005 A comparison of the meshless RBF collocation method with finite element and boundary element methods in neutron diffusion calculations J. R. Xiao , R. C. Batra , D. F. Gilhooley et al ., Analysis of thick plates by using a higher-order shear and normal deformable plate theory and MLPG method with radial basis functions . Comput. Methods Appl. Mech. Eng . 196 , 979 - 987 ( 2007 ). doi: 10.1016/j.cma.2006.08.002 http://doi.org/10.1016/j.cma.2006.08.002 Analysis of thick plates by using a higher-order shear and normal deformable plate theory and MLPG method with radial basis functions S. De , K. J. Bathe , The method of finite spheres . Comput. Mech . 25 , 329 - 345 ( 2000 ). doi: 10.1007/s004660050481 http://doi.org/10.1007/s004660050481 The method of finite spheres M. Moradipour , S. Yousefi , Using a meshless kernel-based method to solve the Black - Scholes variational inequality of American options . Comput. Appl. Math . 1 - 13 ( 2016 ). doi: 10.1007/s40314-016-0351-7 http://doi.org/10.1007/s40314-016-0351-7 Using a meshless kernel-based method to solve the Black - Scholes variational inequality of American options M. Ebrahimnejad , N. Fallah , A. K. R , Three types of meshless finite volume method for the analysis of two-dimensional elasticity problems . Comput. Appl. Math . 36 , 971 ( 2015 ). doi: 10.1007/s40314-015-0273-9 http://doi.org/10.1007/s40314-015-0273-9 Three types of meshless finite volume method for the analysis of two-dimensional elasticity problems S. Atluri , T. Zhu , New concepts in meshless methods . Computational Mechanics. Int. J. Numer. Methods Eng . 47 , 117 - 127 ( 2000 ). doi: 10.1002/(SICI)1097-0207(20000110/30)47:1/3<537::AID-NME783>3.0.CO;2-E http://doi.org/10.1002/(SICI)1097-0207(20000110/30)47:1/3<537::AID-NME783>3.0.CO;2-E New concepts in meshless methods M. Dehghan , M. Abbaszadeh , Numerical investigation based on direct meshless local Petrov Galerkin (direct MLPG) method for solving generalized Zakharov system in one and two dimensions and generalized Gross-Pitaevskii equation . Eng. Comput . 33 , 983 - 996 ( 2017 ). doi: 10.1007/s00366-017-0510-5 http://doi.org/10.1007/s00366-017-0510-5 Numerical investigation based on direct meshless local Petrov Galerkin (direct MLPG) method for solving generalized Zakharov system in one and two dimensions and generalized Gross-Pitaevskii equation A. Taleei , M. Dehghan , Direct meshless local Petrov-Galerkin method for elliptic interface problems with applications in electrostatic and elastostatic . Comput. Methods Appl. Mech. Eng . 278 , 479 - 498 ( 2014 ). doi: 10.1016/j.cma.2014.05.016 http://doi.org/10.1016/j.cma.2014.05.016 Direct meshless local Petrov-Galerkin method for elliptic interface problems with applications in electrostatic and elastostatic M. Dehghan , R. Salehi , A meshless local Petrov-Galerkin method for the time-dependent Maxwell equations . J. Comput. Appl. Math . 268 , 93 - 110 ( 2014 ). doi: 10.1016/j.cam.2014.02.013 http://doi.org/10.1016/j.cam.2014.02.013 A meshless local Petrov-Galerkin method for the time-dependent Maxwell equations M. Kamranian , M. Dehghan , M. Tatari , An adaptive meshless local Petrov-Galerkin method based on a posteriori error estimation for the boundary layer problems . Appl. Numer. Math . 111 , 181 - 196 ( 2017 ). doi: 10.1016/j.apnum.2016.09.007 http://doi.org/10.1016/j.apnum.2016.09.007 An adaptive meshless local Petrov-Galerkin method based on a posteriori error estimation for the boundary layer problems M. Dehghan , M. Abbaszadeh , A. Mohebbi , Meshless local Petrov-Galerkin and RBFs collocation methods for solving 2D fractional Klein-Kramers dynamics equation on irregular domains . C. - Comput. Model. Eng. Sci . 107 , 481 - 516 ( 2015 ). Meshless local Petrov-Galerkin and RBFs collocation methods for solving 2D fractional Klein-Kramers dynamics equation on irregular domains W. H. Chen , X. M. Guo , Element Free Galerkin Method for threedimensional structural analysis . Comput. Model. Eng. Sci . 4 , 497 - 508 ( 2001 ). doi: 10.3970/cmes.2001.002.497 http://doi.org/10.3970/cmes.2001.002.497 Element Free Galerkin Method for threedimensional structural analysis Q. W. Ma , Meshless local Petrov-Galerkin method for two-dimensional nonlinear water wave problems . J. Comput. Phys . 205 , 611 - 625 ( 2005 ). doi: 10.1016/j.jcp.2004.11.010 http://doi.org/10.1016/j.jcp.2004.11.010 Meshless local Petrov-Galerkin method for two-dimensional nonlinear water wave problems J. Amani , K. Parand , Numerical pricing of American options under two stochastic factor models with jumps using a meshless local Petrov - Galerkin method . Appl. Numer. Math . 115 , 252 - 274 ( 2017 ). doi: 10.1016/j.apnum.2017.01.015 http://doi.org/10.1016/j.apnum.2017.01.015 Numerical pricing of American options under two stochastic factor models with jumps using a meshless local Petrov - Galerkin method H. Lin , S. N. Atluri , The Meshless Local Petrov-Galerkin (MLPG) method for solving incompressible Navier-Stokes equations . Comput. Model. Eng. Sci . 2 , 117 - 142 ( 2001 ). doi: 10.3970/cmes.2001.002.117 http://doi.org/10.3970/cmes.2001.002.117 The Meshless Local Petrov-Galerkin (MLPG) method for solving incompressible Navier-Stokes equations S. Y. Long , S. N. Atluri , A meshless local Petrov-Galerkin (MLPG) method for solving the bending problems of a thin plate . CMES 3 , 53 - 64 ( 2002 ). doi: 10.3970/cmes.2002.003.053 http://doi.org/10.3970/cmes.2002.003.053 A meshless local Petrov-Galerkin (MLPG) method for solving the bending problems of a thin plate T. Belytschko , Y. Y. Lu , L. Gu , Element-free galerkin method . Int. J. Numer. Methods Eng . 37 , 229 - 256 ( 1994 ). doi: 10.1002/nme.1620370205 http://doi.org/10.1002/nme.1620370205 Element-free galerkin method G. R. Liu , Y. T. Gu , Meshless local Petrov-Galerkin (MLPG) method in combination with finite element and boundary element approaches . Comput. Mech . 26 , 536 - 546 ( 2000 ). doi: 10.1007/s004660000203 http://doi.org/10.1007/s004660000203 Meshless local Petrov-Galerkin (MLPG) method in combination with finite element and boundary element approaches M. Li , F. F. Dou , T. Korakianitis et al ., Boundary node Petrov-Galerkin method in solid structures . Comput. Appl. Math . 1 - 25 ( 2016 ). doi: 10.1007/s40314-016-0335-7 http://doi.org/10.1007/s40314-016-0335-7 Boundary node Petrov-Galerkin method in solid structures I. Debbabi , H. BelhadjSalah , Analysis of thermo-elastic problems using the improved element-free Galerkin method . Comput. Appl. Math . 1 - 16 ( 2016 ). doi: 10.1007/s40314-016-0401-1 http://doi.org/10.1007/s40314-016-0401-1 Analysis of thermo-elastic problems using the improved element-free Galerkin method D. Mirzaei , R. Schaback , Direct Meshless Local Petrov - Galerkin (DMLPG) method : A generalized MLS approximation . Appl. Numer. Math . 68 , 73 - 82 ( 2013 ). doi: 10.1016/j.apnum.2013.01.002 http://doi.org/10.1016/j.apnum.2013.01.002 Direct Meshless Local Petrov - Galerkin (DMLPG) method : A generalized MLS approximation M. Dehghan , D. Mirzaei , Meshless Local Petrov - Galerkin (MLPG) method for the unsteady magnetohydrodynamic (MHD) flow through pipe with arbitrary wall conductivity . Appl. Numer. Math . 59 , 1043 - 1058 ( 2009 ). doi: 10.1016/j.apnum.2008.05.001 http://doi.org/10.1016/j.apnum.2008.05.001 Meshless Local Petrov - Galerkin (MLPG) method for the unsteady magnetohydrodynamic (MHD) flow through pipe with arbitrary wall conductivity S. Atluri , S. Shen , The basis of meshless domain discretization: the meshless local Petrov-Galerkin (MLPG) method . Adv. Comput. Math . 23 , 73 - 93 ( 2005 ). doi: 10.1007/s10444-004-1813-9 http://doi.org/10.1007/s10444-004-1813-9 The basis of meshless domain discretization: the meshless local Petrov-Galerkin (MLPG) method G. R. Liu , Y. T. Gu , An introduction to mesh free methods and their programming . 1 st ed. ( Springer , Dordrecht , 2005 ) An introduction to mesh free methods and their programming Y. T. Gu , G. R. Liu , A meshless local Petrov-Galerkin (MLPG) method for free and forced vibration analyses for solids . Comput. Mech . 27 , 188 - 198 ( 2001 ). doi: 10.1007/s004660100237 http://doi.org/10.1007/s004660100237 A meshless local Petrov-Galerkin (MLPG) method for free and forced vibration analyses for solids P. Lancaster , K. Salkauskas , Surfaces generation by moving least squares methods . Math. Comput . 37 , 141 - 158 ( 1981 ). doi: 10.1090/S0025-5718-1981-0616367-1 http://doi.org/10.1090/S0025-5718-1981-0616367-1 Surfaces generation by moving least squares methods V. R. Hosseini , E. Shivanian , W. Chen , Local radial point interpolation (MLRPI) method for solving time fractional diffusion-wave equation with damping . J. Comput. Phys . 312 , 307 - 332 ( 2016 ). doi: 10.1016/j.jcp.2016.02.030 http://doi.org/10.1016/j.jcp.2016.02.030 Local radial point interpolation (MLRPI) method for solving time fractional diffusion-wave equation with damping Y. T. Gu , G. R. Liu , A local point interpolation method for static and dynamic analysis of thin beams . Comput. Methods Appl. Mech. Eng . 190 , 5515 - 5528 ( 2001 ). doi: 10.1016/S0045-7825(01)00180-3 http://doi.org/10.1016/S0045-7825(01)00180-3 A local point interpolation method for static and dynamic analysis of thin beams G. Liu , Mesh free methods: moving beyond the finite element method . ( CRC press , Boca Raton , 2003 ) Mesh free methods: moving beyond the finite element method M. Dehghan , D. Mirzaei , Numerical solution to the unsteady two-dimensional Schrödinger equation using meshless local boundary integral equation method . Int. J. Numer. Methods Eng . 76 , 501 - 520 ( 2008 ). doi: 10.1002/nme.2338 http://doi.org/10.1002/nme.2338 Numerical solution to the unsteady two-dimensional Schrödinger equation using meshless local boundary integral equation method J. Zheng , S. Long , G. Li , Topology optimization of free vibrating continuum structures based on the element free Galerkin method . Struct. Multidiscip. Optim . 45 , 119 - 127 ( 2012 ). doi: 10.1007/s00158-011-0667-2 http://doi.org/10.1007/s00158-011-0667-2 Topology optimization of free vibrating continuum structures based on the element free Galerkin method M. Sterk , B. Robic , R. Trobec , Mesh free method applied to the diffusion equation . Parallel Numer . 5 , 57 - 66 ( 2005 ). doi: 10.1007/s00158-011-0667-2 http://doi.org/10.1007/s00158-011-0667-2 Mesh free method applied to the diffusion equation T. Tanbay , B. Ozgener , Numerical solution of the multigroup neutron diffusion equation by the meshless RBF collocation method . Math. Comput. Appl . 18 , 399 - 407 ( 2013 ). doi: 10.3390/mca18030399 http://doi.org/10.3390/mca18030399 Numerical solution of the multigroup neutron diffusion equation by the meshless RBF collocation method D. Mirzaei , M. Dehghan , Meshless local Petrov-Galerkin (MLPG) approximation to the two dimensional sine-Gordon equation . J. Comput. Appl. Math . 233 , 2737 - 2754 ( 2010 ). doi: 10.1016/j.cam.2009.11.022 http://doi.org/10.1016/j.cam.2009.11.022 Meshless local Petrov-Galerkin (MLPG) approximation to the two dimensional sine-Gordon equation D. Jeong , J. Kim , A Crank-Nicolson scheme for the Landau-Lifshitz equation without damping . J. Comput. Appl. Math . 234 , 613 - 623 ( 2010 ). doi: 10.1016/j.cam.2010.01.002 http://doi.org/10.1016/j.cam.2010.01.002 A Crank-Nicolson scheme for the Landau-Lifshitz equation without damping O. Abbasi , A. Rostami , G. Karimian , Identification of exonic regions in DNA sequences using cross-correlation and noise suppression by discrete wavelet transform . BMC Bioinformatics 12 , 430 ( 2011 ). doi: 10.1186/1471-2105-12-430 http://doi.org/10.1186/1471-2105-12-430 Identification of exonic regions in DNA sequences using cross-correlation and noise suppression by discrete wavelet transform O. A. Abuzaid , A. H. M. Gashut , Discontinuous Finite Element Methods for Reactor Calculations . in Third Arab Conf. Peac. Uses At. Energy ( Damascus , 1996 ), pp. 39 - 46 Discontinuous Finite Element Methods for Reactor Calculations IAEA , Status of small and medium sized reactors . ( 2011 ) Status of small and medium sized reactors 10.1007/s41365-018-0506-x 2018-05-05 2018-06-08 2018-07-11 2018-10-16 2018-11-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2018 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-018-0506-x http://dx.doi.org/10.1007/s41365-018-0506-x . 2018 Nuclear Science and Techniques 11 29