Numerical study of the shape parameter dependence of the local radial point interpolation method in linear elasticity

Highlights • The LRPIM is derived from the local weak form of the equilibrium equations for solving a thin elastic plate.• The method LRPIM is used the trial and test functions in the weak form.• Convergence of the LRPIM depends on number of parameters derived from local weak form and sub-domains.• The effect of distributions nodes number by varying nature of material and the RBF-TPS.• The calculated results are compared with the analytical solution of the deflection.


Method detail
Meshless method has attracted more and more attention from researchers in recent years.It is regarded as a potential numerical method in computational mechanics. Several meshless methods, such as smooth particle hydrodynamics (SPH) method [1,2], element free Galerkin (EFG) method [3], meshless local Petrov-Galerkin (MLPG) method [4][5][6][7][8], point interpolation method (PIM) [9,10], local point interpolation method (LPIM) [11] and local radial point interpolation method (LRPIM) has been proposed by Liu et al. [10,12,13]. In LRPIM, the point interpolation using the radial basis function to construct the shape functions which have the delta function property. The radial basis function (RBFs) is thin plate spline (TPS) [14,15]. Local weak forms are developed using weighted residual method locally from the partial differential equation of linear elasticity of 2D solids. The number of numerical examples will be presented to demonstrate the convergence and accuracy, validity and efficiency of the present methods. The local radial point interpolation method LRPIM is a Meshless method with properties of simple implementation of the essential boundary conditions and with the lower cost than moving least squares (MLS) methods. This paper deals with the effect of sizing parameter of subdomains on the convergence and accuracy of the methods. Numerical values will be presented to specify the convergence domain by precising maximum and minimum values as a function of distribution nodes number and by using the radial basis functions TPS (thin plate spline). It also will present a comparison with numerical results for different materials. The analytical solution of the deflection confirms the numerical results. The LRPIM method will be developed to solve the problem of a thin elastic homogenous plate. The local weak form and numerical implementation are presented in section 3, numerical example for 2D problem are given in section 4. Then, the paper will end with the results, discussions and finally the conclusions.

RPIM shape functions in meshless method
u h ðxÞ is composed of two parts: P j ðxÞPolynomial basis functions and R i ðxÞ the radial basis functions RBFs [16][17][18][19][20]: n is the number of field nodes in the local support domain and m is the number of polynomial terms. Radial basis is a function of distance r: The above Eq. (1) can be expressed in the matrix form [17] where U 1 is the vector of function values: U 1 ¼ ½u 1 ; u 2 ; u 3 ; . . . ; u n T R The moment matrix of RBFs, P is the moment matrix of Polynomial basis function and a, b are the values of unknowns coefficients (Radial and Polynomial).
We note that, to obtain the unique solutions of Eq. (2), the constraint conditions should be applied as follows [21]: the combining of Eqs. (3) and (4) yields a set of equations in the matrix form: The unknown vector in the equation above can be obtained by the inversion of the matrix: G ¼ Substitution of the vector obtained by inversion of matrix G into Eq. (1) leads to: Local weak form method LRPIM Let us consider a two-dimensional problem of solid mechanics in domain V bounded by Γ whose strong-form of governing equation and the essential boundary conditions are given by: where in V: s T ¼ ½s xx ; s yy ; t xy is the stress vector, b T ¼ ½b x ; b y the body force vector. n ¼ ðn 1 ; n 2 Þ is the vector of unit outward normal at a point on the natural boundaries. t 0 is the prescribed effort, ½u 1 ; u 2 are the displacement components in the plan and ½u 0 1 ; u 0 2 are the prescribed displacement on the essential boundaries.
In the local Petrov-Galerkin approache [3], one may use a weak form over V Q a local quadrature domain (for node I), which may have an arbitrary shape, and contain the point x Q in question, see Fig. 1. The generalized local weak form of the differential Eq. (7) is obtained by: where V Q is the local domain of quadrature for node I and y I is the weight or test function (y I 2 C K ðVÞ) [4]. Generally, in meshfree methods, the representation of field nodes in the domain will be associated to other repartitions of problem domain: V I is the influence domain for nodes interpolation, V S is the support domain for accuracy. For each node V y is the weight function domain and V Q is the quadrature domain for local integration.
Using the divergence theorem [4] in Eq. (10), we obtain: The internal boundary of the quadrature domain G Qt : The part of the natural boundary that intersects with the quadrature domain G Qu : The part of the essential boundary that intersects with the quadrature domain We can then change the expression of Eq. (11): Using the radial point interpolation method (RPIM) shape functions (see sub-section 2), we can approximate the trial function for the displacement at a point x (8x 2 V S ) as Eq. (6) The stress vector is defined by: where C is the symmetric elasticity tensor of the material where V I ¼ y I;x 0 0 y I;y y I;y y I;x 0 @ 1 A is a matrix that contains the derivatives of the weight functions and V ¼ y I 0 0 y I is the matrix of weight function.
Substituting the differential operator L d ¼ @=@x 0 0 @=@y @=@y @=@x 0 @ 1 A into Eq. (13) we obtain: Substituting Eqs. (15) and (16) into Eq. (14), we obtain the discret systems of linear equations for the node I. (17) The matrix form of Eq. (17) can be written as in matrix form: where K I is a 2 Â 2 matrix called a nodal stiffness matrix, given by and nodal force vector with contributions from body forces applied in the problem domains: where n 0 denote the set of the nodes in the support domain V S of point x Q .
Two independent linear equations can be obtained for each node in the entire problem domain and by assembling all these n Â 2 equations to obtain the final global system equations: To solve the precedent system, the standard Gauss quadrature formula is applied with 16 G points [3,22] for calculating integrals in Eqs. (19) and (20) on both boundary and domain.
The size of quadrature domain is specified by setting a Q ¼ 2 and a regular distribution of nodes on the mid-surface of plate in (x, y) plane is employed.

Numerical 2D elastostatic example
This section is about numerical results for a cantilever rectangular plate see (Fig. 2). First, the effects of the size of support and quadrature domains were investigated and the convergence of LRPIM method for several materials was examined numerically; then, comparisons will be made with the analytic solution for several materials [23] We choose: steel, zinc, aluminium and copper with:  The sizes of support domain V s (quadrature domain V Q resp.) are defined by: d s ¼ a s d c (r Q ¼ a Q d cI resp.) where d c (d cI resp) is the nodal spacing near node I (see Fig. 3) and a s (a Q resp) is the size of the support domain V s (local quadrature domains resp) for node I. The sizes of support domain V s (quadrature domains resp) will be respectively determined in x and y directions. For simplicity, we put: a Sx ¼ a Sy ¼ a S (a Qx ¼ a Qy ¼ a Q resp) are used for V s (V Q resp).

Results and discussions
The standard Gaussian quadrature formula with 16 G points and the radial point interpolation method (RPIM) approximation, linear polynomial basis functions are applied. The cubic or quadratic spline functions a used as the test functions in the LRPIM local weak-form.
Throughout this section and for all calculations a Q ¼ 2 was fixed.

Simulation results and comparison with the radial basis function RBF-TPS (of a thin plate spline)
The use of radial basis function RBF-TPS have not been extensively studied on literature. We give in this paragraph numerical results for different materials.    For n t ¼ 175, the convergence domain is: 0:5 < h < 6:5 For n t ¼ 189, the convergence domain is: 1 < h < 6:5 The results illustrated in the Figs. 5-7 to the number of nodes n t = 55, 91, 175 and 189 with the radial basis RBF-TPS give the domain of convergence large enough for the different types of materials studied. This shows better convergence than that given in references [10,24,25]. These authors give a single value (h ¼ 4:001). We show here that the values h can be varied in a broader segment  depending on the number n t .We found that the domain of convergence is more larger than the value given by Liu et al. Fig. 8 shows the variations of displacement as a function of x 1 for x 2 ¼ 0 with the shape parameter h (h = 5 and 6) of the radial basis function RBF-TPS. We used two studied materials: Steel and Zinc. We note that there is a coincidence between the analytical solution and the result obtained by the LRPIM method of the radial basis RBF-TPS. It shows the convergence of LRPIM method (see Fig. 7)

Conclusion
In this paper, the meshless LRPIM method is employed for solving a 2D elastostatic problem. The LRPIM method depends on the sizing parameter a S , which is associated to different parameters coming from the weak form formulation. In the study a Q ¼ 2 is fixed, the nature of convergence domain as a function of a S and the effect of distributions nodes number n t by varying nature of material and the radial basis functions RBF-TPS. We conclude that for small value of n t ð55Þ lead to the upper extremity of convergence domain which is limited to a S ¼ 5. The results obtained for the number of nodes n t = 55, 91, 175, 189 and for the radial basis RBF-TPS, give a domain of convergence which is large enough for the different types of materials studied. This shows a better convergence more than the result that is found in the references [10,24,25]. There is a good agreement with the analytical solution of the deflection.