Abstract: Mesh free methods are applied in many new researches in the areas of approximation theory and numerical analysis as an alternative to the mesh based methods such as finite element and finite difference. The recent new developed mesh free methods have the advantage of not requiring a mesh of elements while providing comparable results to the other numerical methods. In this study, a mesh free scheme based on the radial point interpolation method was used to solve the one dimensional neutron diffusion equation. The applied method uses the Galerkin weak form of the differential equations. Radial basis functions which are powerful functions in the field of function approximation, were used to construct the shape functions. Gauss quadrature scheme was applied in order to calculate the integrations of the weak form of the equations. The efficiency and accuracy of the method was evaluated through a number of examples. The Reed test problem was also solved through the applied method to demonstrate the capabilities of the developed program. The applied Mesh free method results were in good agreement with the analytical solutions. Comparing with the results of the finite element, one concludes the applied mesh free method improves the computational accuracy but at a cost of requiring more computational time.
INTRODUCTION
Recently, there has been great interest in the development of computational methods able to avoid mesh dependency, provided the results to be comparable with the currently established methods. Common approximation methods such as finite elements require underlying mesh to define basis functions. This is usually a difficult process to implement in space dimensions. A group of techniques which are called mesh free methods, have been developed in the last decades and successfully applied in a wide variety of computational fields (Leitao et al., 2007). In mesh free methods a set of scattered nodes are applied to represent the problem domain and its boundaries. These sets of nodes which are called field nodes, do not form a mesh of elements. On the other hand, for the interpolation or approximation of unknown functions no priori information on the relationship between the nodes is required (Liu and Gu, 2005).
Mesh free methods have a number of interesting properties. They allow an accurate representation of complex problem domain and there is no pre-defined connectivity between the nodes (Avila and Perez, 2008; Ding et al., 2004). Another exciting property of mesh free methods is their potential in adaptive techniques where a user in an interactive mode could simply add a number of nodes into the portion of problem domain where more accuracy is required. Then, the method with a simple implementation process enhances the results without need to construct a new mesh of elements (Belytschko et al., 1994).
Several mesh free methods have been developed up to now, among them the Smoothed Particle Hydrodynamics (Liu and Liu, 2003), the Element Free Galerkin Method (Belytschko et al., 1994), Point Interpolation Method (Liu and Gu, 2001), the hp-clouds method (Duarte and Oden, 1996), the Reproducing Kernel Particle Method (Liu et al., 1995), Generalized Finite Element Method (Babuska et al., 2004) and the Partition of Unity Method (Babuska and Melenk, 1997) are the better known methods.
Mesh free methods are classified based on the form of the governing equations (strong, weak or a combination of both week and strong) and also on the applied function approximation/interpolation scheme. Mesh free methods which are based on a strong form of equations do not require an integration process. These techniques are truly mesh free methods (Liu, 2003).
Mesh free methods based on the weak form of equations need the integration process. These methods require background cells to perform integration. In this case, contrary to the Finite Element Method (FEM), there is no relation between nodes and background cells. Currently, the mesh free weak form is the most widely used due to its excellent stability (Liu et al., 2006). Regarding classification based on the applied function approximation/interpolation scheme, the Moving Least Square (MLS) approximation and Point Interpolation Method (PIM) are the most widely used methods (Liu et al., 2004). The Element Free Galerkin (EFG) and Mesh Less Petrov-Galerkin (MLPG) methods have been developed based on MLS approximation. Contrary to the MLS approximation, PIM uses interpolation to construct shape functions that possess the Kronecker delta function property which is useful for implementing the boundary conditions. (Cui et al., 2010).
Although, the mesh free methods have been applied successfully in vast variety of areas, most are still under development. According to the authors’ search results, to date, the mesh free methods have not been applied in the field of reactor physics calculations. In this study, a mesh free method based on the Radial Point Interpolation Method (RPIM) is implemented to solve the neutron diffusion equation in one dimension. The Galerkin method was applied to discretise the neutron diffusion equation and the Gauss-Legendre scheme is used to carry out the integrations of the weak form of the equations. A computer code was developed in Matlab environment for implementation of the method and several one dimensional test problems were solved through this computer code. The results were compared with the FEM results and the analytical solutions.
MESH FREE RADIAL POINT INTERPOLATION METHOD
Since the mesh free methods do not use mesh of elements, the field variable φ(x) at a point of interest x is interpolated in the problem domain using the function values at neighbor nodes (which are called local support domain) of the point x:
| (1) |
where, n is the number of neighbor nodes considered in the local support domain of the point x, φi is the nodal field variable at field node i and Ni(x) is the shape function corresponding to field node i. A local support domain of a point x determines the number of nodes(n) to be used to approximate the function value at x.
Similar to FEM, the discrete equations of a Mesh free method can be formulated using the shape functions and the strong or weak form of the problem equations.
Radial basis functions: Over the last decades, Radial Basis Functions (RBFs) have gained popularity and found to be widely successful in the field of scatter data interpolation (Sarra, 2006). More recently, methods based on the RBFs have been applied for the numerical solution of partial differential equations (Kansa, 1990). RBFs were also successfully applied to neural networks (Ozyilmaz et al., 2002; Farivar et al., 2009; Qasem and Shamsuddin, 2010).
In radial basis functions, the variable is only the distance between the point of interest x and a node at xi (r = ||x-xi||). There are a variety of RBF types and their characteristics have been widely investigated (Buhmann, 2003). Radial basis functions may or may not contain free parameters called the shape parameters. The most widely used RBFs are listed in Table 1 (Liu and Gu, 2005).
The mentioned RBFs in the above table which are called classic RBFs, have global support. In addition, so-called compactly supported radial basis functions such as the Wendland functions (Wendland, 1995) have also been developed. Some study on mechanics problems failed to find clear advantages of compactly supported RBFs over classic RBFs (Liu and Gu, 2005).
Radial point interpolation shape function: One of the main differences of the Mesh free methods and FEM is at the stage of shape function construction. Interpolation of function φ(x) through the radial basis functions can be written as:
| (2) |
where, Ri (r) is a RBF, n is the number of RBFs, Pk(x) is a monomial in the space coordinate and m is the number of polynomial basis functions, ai and bi are the coefficient to be determined.
The polynomial term in Eq. 2 is not always necessary. A number of advantages such as improvement in accuracy and stability of interpolations have been found for adding the polynomial terms (Liu and Gu, 2005). Due to these advantages, in this study RBFs augmented with linear order polynomial term have been used to construct the shape functions.
In order to obtain the shape functions, Coefficients a and b in Eq. 2 should be determined. For this purpose a support domain including a number of field nodes is formed for the point of interest x. Coefficients ai and bk in Eq. 2 can be determined by enforcing Eq. 2 to be satisfied at the n nodes surrounding the point of interest x (local support domain). This leads to n linear equations, one for each node. The matrix form of Eq. 2 can be expressed as:
| (3) |
where, φs is the function values at the nodes.
| Table 1: | The most widely used radial basis functions |
| dc: Average distance of failed nodes in the local support domin | |
The moment matrix of RBFs and the polynomial moment matrix can be expressed as:
| (4) |
| (5) |
In Eq. 4, ri in Rj (ri) is: |xi-xj|.
However, there are n + m unknowns in Eq. 3 (a1,…,an and b1,…,bm). The additional m equations can be added using the following m constraint conditions (Liu et al., 2007):
| (6) |
Combining Eq. 3 and 6 yields the following set of equations in the matrix form:
| (7) |
| (8) |
From Eq. 2 and 7 we can obtain:
| (9) |
Where:
| (10) |
Considering that:
Eq. 9 can be re-written as:
| (11) |
According to Eq. 11, the shape functions corresponding to the field nodes are as follows:
| (12) |
These shape functions are found through Eq. 10 by omitting the extra elements of the obtained vector. More details and discussion about the properties of the above shape functions which are called RPIM shape functions can be found by Liu and Gu (2005).
In this research, different types of RBFs including MQ, EXP and TPS augmented with the linear order polynomial term have been used to construct the shape functions.
WEAK FORM OF THE NEUTRON DIFFUSION EQUATION
The standard weak form of the neutron diffusion equation is derived by multiplying the governing equation to the weight function v(x) and integration over the domain:
| (13) |
Using the integration by part, Eq. 13 is reduced to the following form:
| (14) |
v (x) should satisfy the essential boundary condition.
IMPLEMENTATION OF THE METHOD
In order to implement the described method, the Galerkin method is used to
discretise the weak form of the diffusion equation. For this purpose v (x),
φ (x) and
| (15) |
| (16) |
| (17) |
where, n is the number of nodes in the local support domain, N(x) is the RPIM shape function (Eq. 10, 12) and φj is the field variable value (neutron flux) at point j.
The integration in the left hand side of Eq. 14, considering the above replacements, can be written as:
| (18) |
where, n is the number of field nodes in the local support domain and:
is called the nodal stiffness matrix.
The integration in the right hand side of Eq. 14, similar to that shown above can be written as:
| (19) |
where fi is the nodal body force vector.
The last two terms in Eq. 14 are specified by the boundary conditions defined in the problem. As the RPIM shape functions possess the Kronecker delta function property, the boundary conditions can be imposed directly and accurately without any additional treatment compared with FEM (Liu, 2003).
In order to calculate the stiffness and body force integrations, the problem domain is discretised into a set of background cells and the global integrations can be expressed as a summation of integrals over these cells:
| (20) |
where, nc is the number of background cells, G represents the integrand and Ωk is the domain of the kth background cell.
Contrary to FEM, there is no relation between the background cells and the field nodes. This property simplifies the node insertion process in mesh free methods compare to that in FEM.
Several methods to perform the numerical integrations in the Mesh free methods have been suggested up to now (Liu et al., 2007; Dolbow and Belytschko, 1999). However, in this research the Gauss quadrature scheme which is commonly used in the FEM and possess good properties regarding the accuracy and stability, was employed to perform the numerical integrations over the cells. When ng Gauss points are used in each background cell, the integration for each background cell in one dimension is:
| (21) |
where, wi is the Gauss weighting factor for the ith Gauss point at xQi.
NUMERICAL RESULTS AND DISCUSSION
In order to demonstrate the ability of applied mesh free method to deal with neutron diffusion equation and consequently to demonstrate the method upon which it is based, a number of test cases in one dimension was solved. The existence of analytical solutions was the main criteria in selection of the problems. A computer code was developed in Matlab software to implement the method. The mesh free results were compared with the FEM and analytical solutions. In the implementation of the mesh free method three different RBF (MQ, EXP and TPS as introduced in Table 1) were used. Furthermore, two different sizes of local support domain were considered to construct shape functions. In the first situation which is called “2n MFM”, 2 nodes were regarded in the local support domain and in the second situation which is called “4n MFM”, 4 nodes considered in the local support domains. In order to keep the similarity of the solutions to compare the results with the FEM, linear (2 nodes in each element; which is called 2nFEM) and cubic (4 nodes in each element; which is called 4nFEM) shape functions are used for the FEM treatment of problems. In addition FEM based on the Galerkin weak form of the equations is used. The following definition of the error was applied in the calculations:
| (22) |
Case one: Consider a 2 region finite slab of non-multiplying with 5 cm length of each region and a constant plane source (S = 1) at the left side (at x = 0), subject to the following boundary conditions:
| • | Left B.C.: |
| • | Right B.C.: φ = 0 at x = 10 cm |
And considering the continuity of the neutron flux and current at the interface (x = 5) and the following constants and node numbers:
Case I:
| • | Region 1: | Σα = 0.5 cm-1, D = 0.33 cm, ne = 6 |
| • | Region 2: | Σα = 0.1 cm-1, D = 0.56 cm, ne = 6 |
Case II:
| • | Region 1: | Σα = 1.0 cm-1, D = 0.22 cm, ne = 11 |
| • | Region 2: | Σα = 0.1 cm-1, D = 0.56 cm, ne = 11 |
Case III:
| • | Region 1: | Σα = 0.01 cm-1, D = 0.65 cm, ne = 9 |
| • | Region 2: | Σα = 1.0 cm-1, D = 0.22 cm, ne = 9 |
Figure 1a shows the obtained results from the Mesh free method, FEM and the analytical solution for case5-1-III. Table 2 shows the error of the Mesh free and FEM methods defined by Eq. 22.
Case two: Consider a 2 region finite slab of non-multiplying with 5 cm length of each region and a distributed source (S = 1) on the entire length of region 1 of slab, subject to the following boundary conditions:
| • | Left B.C.: Perfect Reflector: ( |
| • | Right B.C.: Bare condition or equally Extrapolated length B.C. (δ = 2.13D) |
| Fig. 1: | Calculated flux for cases (a) case one-III and (b) case two-III |
| Table 2: | Calculated error for the cases one and two |
And considering the continuity of the neutron flux and current at the interface (x = 5 cm) and the following constants and node numbers:
Case I:
| • | Region 1: | Σα = 0.5 cm-1, D = 0.33 cm, ne = 6 |
| • | Region 2: | Σα = 0.1 cm-1, D = 0.56 cm, ne = 6 |
Case II:
| • | Region 1: | Σα = 1.0 cm-1, D = 0.22 cm, ne = 6 |
| • | Region 2: | Σα = 0.1 cm-1, D = 0.56 cm, ne = 6 |
Case III:
| • | Region 1: | Σα = 1.5 cm-1, D = 0.17 cm, ne = 9 |
| • | Region 2: | Σα = 1.0 cm-1, D = 0.22 cm, ne = 9 |
Figure 1b shows the obtained results from the Mesh free method, FEM and the analytical solution for case 5-2-III. Table 2 shows the error of the Mesh free and FEM methods defined by Eq. 22:
Calculated error for the Cases one and two: The only difference between ‘2nMFM’ and ‘4nMFM’ methods is the size of local support domains or in other words the number of nodes in the local support domains. In 2nMFM, 2 nodes were used in each local support domain while in 4nMFM, 4 nodes were applied. The abbreviations MQ, EXP and TPS in Table 1 are the RBF types. The 4nMFM indicated results on Fig. 1 obtained using EXP-RBFs. Since the 2nMFM implementation through different RBF types show the same results in all cases, the RBF type extension of the methods was omitted and only one column selected to indicated it. Both 2nMFM and 2nFEM use 2 nodes to construct the shape function and as it is clear in Table 2, the outcome of both methods are very close to each other. This is due to the fact that in the 2 nodes simulation, the linear polynomial which is augmented to the RBFs plays a dominant role in the interpolations. As indicated in Table 2, increasing nodes in the local support domains which is easy to implement process in mesh free methods, enhances the accuracy of the results. The 4nMFM results through different types of RBFs in all cases show higher accuracy in comparison to 4nFEM results. This is due to the power of RBFs in the interpolation field. Among the applied mesh free methods with 4 nodes in the local support domain, the 4nMFM-EXP results are more accurate than the others. The parameter ‘e’ in Table 2 is defined as the ratio of the 4nMFM-EXP error to the 4nFEM error. Illustrated results in Fig. 1 and Table 2 show that the applied mesh free method gives comparable results to 2nMFM and more accurate results are estimated using 4nMFM. Increasing the number of nodes to construct the shape function in mesh free methods is carried out without causing complexity of procedure while this process in FEM is performed through modifying the order of polynomial.
The effect of node numbers in the problem domain: The impact of applying different node numbers in the problem domain on the accuracy of the results were studied through examining a number of the above mentioned case studies. Fig. 2a and b show the solution errors (defined according to Eq. 22) for the case studies case one-II and case two-III, respectively.
| Fig. 2(a-b): | The effect of node numbers on the accuracy of results (a) case one-II and (b) case two-III |
| Fig. 3(a-b): | Calculated flux for problem (a) case one-II and (b) case two-III with irregular nodes |
| Table 3: | Calculated errors for a number of case studies with the irregular node distribution |
Using less node numbers in problem domain, one find mesh free results are more accurate than FEM. By increasing the number of nodes in the problem domain, the FEM results approach to the mesh free results.
Irregular field node distribution: In order to evaluate the robustness of the method to the irregularity of the field nodes distribution, some of the above problems were solved with a different distribution of field nodes. The robustness of the Mesh free method was evaluated as comparable with the FEM results. The obtained results for case studies case one-III and case two-III with irregular distribution of the field nodes are indicated in Fig. 3a and b, respectively and also in Table 3. These case studies were also solved through the FEM method with the same irregularity of nodes and the results are indicated in Fig. 3a and b accordingly and also in Table 3 for comparison with the Mesh free results. In the regular case the nodes were regularly distributed between 0 and 10 with the interval of 0.625 cm and in the irregular case the nodes were at: {0 0.4 1.3 1.7 2.1 3.0 3.7 4.1 5.0 5.3 6.2 6.7 7.6 8.4 8.7 9.2 10}.
Calculated errors for a number of case studies with the irregular node distribution: It is clear from the results, in the irregular cases, the errors of 2nMFM and 2nFEM solutions are the same and the 4nMFM errors are less than 4nFEM. Comparing the error of methods and the parameter ‘e’ in Table 2 and 3 indicate that whether the irregularity of nodes increases the errors in a number of cases, the ratio of errors is decreased in most of the cases. This shows the excellent behavior of the applied mesh free method at the irregularity of node distribution compared with the applied FEM.
Node insertion: Node insertion is one of the exciting characteristics of Mesh free methods where a user in an interactive mode could simply insert a large number of nodes into the portion of problem domain where more accuracy is required and the method would precede to enhance the results.
| Fig. 4: | The effect of node insertion in the problem domain |
| Table 4: | Specification of the 1D reed problem |
Figure 4 shows the obtained results for case study 5-2 with the following constants:
| • | Region 1: | Σα = 1.5 cm-1, D = 1/6 cm |
| • | Region 2: | Σα = 1.0 cm-1, D = 1/4.5 cm |
For this case, adding only two more nodes to the problem domain besides interface of the regions (at x = 4.5 and 5.5 cm), enhance the accuracy of the results as shown on Fig. 4 (the problem was solved through 4nMFM-EXP implementation). The process of node insertion in Mesh free methods is done simply be defining the position of new nodes and without any other change in the other parts of solution.
The reed problem: This one group multi-region problem for an idealized reactor lattice cell that is located at the edge of a bare core provides a severe test for transport theory codes. This problem is usually considered to check the validity of transport codes but the diffusion solution of this problem could be comparable to the P1 solution. It is a problem of four different material regions with a reflecting boundary condition on the left and a vacuum boundary condition on the right. This problem was solved to demonstrate the capabilities of the developed computer program. The specification of the Reed problem including region width, material properties and source strength of each region is given in Table 4 (Abuzaid, 1994). The number of nodes which were used in different regions of the problem domain is indicated in Table 5.
| Fig. 5: | The obtained results of the 1D Reed Problem |
| Table 5: | The number of nodes applied to solve the 1D Reed Problem |
Figure 5 shows the obtained results of solution of the Reed Problem through the applied mesh free and FEM methods. The results are in good agreement with the P1 solution of this problem obtained by Abuzaid (1994).
CONCLUSION
The applied Mesh free method gave the same results as FEM where 2 nodes were used in the local support domains and led to more accurate results than FEM where 4 nodes were used in the local support domains. The error of the applied Mesh free methods with 2 nodes in local support domains using different types of RBFs were the same but among the applied Mesh free methods with 4 nodes in the support domains, the 4nMFM-EXP results were more accurate than the others. Adding more nodes in the problem domain (near to the interfaces and also everywhere which more accuracy is required) enhances the accuracy of results. This is an easy to implement process in mesh free methods compare to mesh based methods, especially in space dimensions. Robustness of the applied method to irregular node distribution was examined and good performance was revealed in this area. The main disadvantage of the Mesh free methods is that the shape function construction is a more complicated and more time consuming process compared to FEM but it should be considered that nowadays, increases in the processing speed of computers and also parallel computation can compensate this drawback of the Mesh free methods.