Some Remarks on the Method of Fundamental Solutions for Two-Dimensional Elasticity

: In this paper, some remarks for more efficient analysis of two-dimensional elastostatic problems using the method of fundamental solutions are made. First, the effects of the distance between pseudo and main boundaries on the solution are investigated and by a numerical study a lower bound for the distance of each source point to the main boundary is suggested. In some cases, the resulting system of equations becomes ill-conditioned for which, the truncated singular value decomposition with a criterion based on the accuracy of the imposition of boundary conditions is used. Moreover, a procedure for normalizing the shear modulus is presented that significantly reduces the condition number of the system of equations. By solving two example problems with stress concentration, the effectiveness of the proposed methods is demonstrated.


Introduction
Nowadays, the method of fundamental solutions (MFS) is widely used by engineers and scientists in a variety of applications. The MFS has a great capability to provide very accurate solutions [Alves (2009) ;Tsai and Young (2013)] without considering any unknown variable within the domain of the problem. A few different boundary-type meshfree techniques have also been proposed and developed which are different but related to the MFS, see for example [Zhang, Gu, Hua et al. (2018) ;Liu, Wang and Qu (2018); Gu, Fan and Xu (2019)]. Many important meshfree methods based on weak form formulation have been proposed, in which efficient evaluation of domain integrals requires special techniques [Chen, Yoon and Wu (2002); Hematiyan, Khosravifard and Liu (2014) ;Jamshidi, Hematiyan and Mahzoon (2019)]; however, the MFS is an integration-free strong-form meshfree method. Besides these advantages, the MFS involves some difficulties, which have restricted its fast development, especially in engineering fields. The proper determination of the location of source points and solving the resulting ill-conditioned system of equations are important issues, which should be dealt with in the MFS. Until now, several researches have been conducted for proper determination of the location of source points in the MFS. Tsai et al. [Tsai, Lin, Young et al. (2006)] investigated on the MFS for Laplace, Helmholtz, modified Helmholtz, and biharmonic equations and found out that more accurate solutions can be obtained by considering source points farther from the boundary. They also mentioned that the condition number of the system of equations increases by increasing the distance between main and pseudo boundaries and concluded that this distance should be limited with an attention to the capability of the system solver. Gorzelańczyk et al. [Gorzelańczyk and Kołodziej (2008)] investigated on the analysis of the torsion problem with the Poisson's governing equation using the MFS, where they employed the fundamental solution of the two-dimensional Laplace equation. They examined circular pseudo boundaries and source contours geometrically similar to the main boundary and concluded that the latter leads to much more accurate results. Karageorghis [Karageorghis (2009)] proposed an optimization method with one variable for determining the suitable distance between pseudo and main boundaries in the MFS. He examined harmonic and biharmonic boundary value problems with different boundary conditions. Wong et al. [Wong and Ling (2011)] proposed an optimization method for determining the suitable configuration of the pseudo boundary in the MFS for Laplace equation. They considered two cases with one optimization variable. In the first case, the distance between main and pseudo boundaries was the unknown of the optimization problem, while in the second case, this distance was fixed and the number of source points was the unknown of the optimization problem. Liu [Liu (2012)] introduced an optimization method for location of source points in the MFS for the Laplace equation. He considered one unknown distance factor for each source point and could find the unknowns by solving an uncoupled system of nonlinear equations. Li et al. [Li, Chen and Karageorghis (2013)] investigated on different configurations of collocation and source points in the MFS for 2D Laplace equation. They observed that placing the source points on a circle leads to poor results and recommended pseudo boundaries similar to the main boundary of the problem. Chen et al. [Chen, Karageorghis and Li (2016)] investigated on the configuration of source points in the MFS for 2D and 3D boundary value problems governed by the Laplace and biharmonic equations. They used an optimization algorithm based on the accuracy of the imposition of boundary conditions for proper determination of the location of source points. Hematiyan et al. [Hematiyan, Haghighi and Khosravifard (2018)] proposed a method for determining the location of collocation and source points in the MFS for the 2D Laplace equation without any optimization process. They introduced a new index, namely "the location parameter" for source points that expresses the location of the source point relative to its neighboring source and collocation points. They showed that values greater than 0.8 are suitable for the location parameter of source points in the MFS for the 2D Laplace equation. This method was successfully used by Mohammadi et al. [Mohammadi, Hematiyan and Khosravifard (2019)] for selecting the position of source points in the MFS for steady-state heat conduction analysis. Grabski et al. [Grabski and Karageorghis (2019)] applied the MFS for potential boundary value problems, in which both intensity and position of source points were unknown. This approach yields a system of nonlinear equations. Recently, Grabski [Grabski (2019)] examined four different configurations of source points for analysis of transient heat conduction problems using the MFS. He observed that the results were more accurate in the cases that source points were located further from the main boundary of the problem. As it is seen, most of the above-mentioned researches on the location of source points in the MFS are devoted to boundary value problems with the Laplace, Helmholtz, and biharmonic equations. The investigations on the configuration of source points in the MFS for elastostatic problems have been very rare; however, the MFS formulation for 2D and 3D elastostatic problems have been considered by several researches that are briefly reviewed here. Berger et al. [Berger and Karageorghis (2001)] used the MFS for analysis of elastostatic problems in 2D isotropic and anisotropic domains. They analyzed domains with one or two materials. Poullikkas et al. [Poullikkas, Karageorghis and Georgiou (2002)] applied the MFS to 3D isotropic elastostatic problems. Similar to the work by Berger et al. [Berger and Karageorghis (2001)], components of pseudo forces and their locations were considered as unknowns of the problem in this work. De Medeiros et al. [De Medeiros, Partridge and Brandão (2004)] formulated the MFS with the dual reciprocity method for analysis of 2D elastostatic problems involving body forces. They considered circular pseudo boundaries and used singular value decomposition for solving the MFS system of equations. Marin et al. [Marin and Lesnic (2004)] and Marin [Marin (2005)] used the MFS for the Cauchy problem in 2D and 3D elasticity, respectively. They observed that the error of solutions decreases by increasing the distance between pseudo and main boundaries. Fam et al. [Fam and Rashed (2005)] used the MFS with a new set of particular solutions to solve 3D elastostatic problems involving body forces. They derived particular solutions in explicit forms for displacements and stresses using a decoupling technique. They located the source points on a pseudo boundary similar to the main boundary of the problem. Tsai [Tsai (2009)] developed the MFS with the dual reciprocity for analysis of 3D thermoelastic problems involving body forces. He observed that using truncated singular value decomposition (TSVD) could improve the accuracy of the obtained results. Marin et al. [Marin and Karageorghis (2013)] formulated the MFS with the method of particular solutions for analysis of thermoelastic problems. They derived new particular solutions for the nonhomogeneous equilibrium equations. Karageorghis (2013a, 2013b); Marin, Karageorghis and Lesnic (2015)] used the MFS for inverse analysis of inverse thermoelastic problems too. Sun et al. [Sun and Marin (2017)] proposed an invariant method of fundamental solutions for analysis of 2D elastostatic problems. They used the Tikhonov regularization method and the discrepancy principle and obtained very accurate solutions. Askour et al. [Askour, Tri, Braikat et al. (2018)], by combination of the MFS, the asymptotic numerical method, and the analog equation method, developed a technique for analysis of nonlinear elastic problems. They also used the TSVD and Tikhonov regularization methods for solving the resulting ill-conditioned system of equations. As it can be seen from the above review, several researches have been conducted on the configuration of source points in the MFS analysis of boundary value problems with the Laplace, Helmholtz, and biharmonic equations; however, to the authors' best knowledge, no research devoted to the configuration of source points for the MFS analysis of elastostatic problems has been reported yet. In the present work, a procedure for determining a suitable configuration of source points in the MFS for 2D elastostatic problems is presented. A criterion based on the accuracy of imposition of boundary conditions is employed for solving the resulting ill-conditioned system of equations using the TSVD. Moreover, a method for reducing the condition number of system of equations by normalizing the shear modulus is suggested. By analyzing two example problems with stress concentration, it is shown that the proposed methods can be efficiently used for analysis of complicated 2D elastostatic problems.

The MFS for 2D elastostatic problems
In the absence of body forces, the governing equations for elastostatic problems are expressed as follows: where σ and ε are the stress and strain tensors, respectively; u is the displacement vector, ν and G are the Poisson's ratio and shear modulus, respectively, and x is a point within the problem domain Ω or over its boundary Γ . The shear modulus G can be expressed in terms of the elastic modulus E as ) 1 ( 2 / ν + = E G . The unknowns in a 2D elastostatic problem are 1 u , 2 u , 11 ε , 12 ε , 22 ε , 11 σ , 12 σ , and 22 σ .
Boundary conditions at a boundary point of a 2D elasto-static problem can be expressed as follows: where n u and t u are normal and tangential displacement, respectively, and n t and t t are normal and tangential components of the traction vector at the boundary point. n u , t u , n t , and t t are given functions. As it is seen from Eqs. (4) and (5) , while other functions of α f and α g are zero. We describe the MFS formulation for plane strain problems only. A plane stress problem formulated in terms of the elastic constants ν and G can be solved by a solver of plane strain problems by changing ν into ) 1 ( ν ν + , while the shear modulus is not changed. In the MFS with N sources (fictitious concentrated force), the displacement solution of the 2D elastostatic problem is expressed as follows: In Eq. (9), the indices i and j can be 1 or 2. By substituting Eq. (8) into Eq. (2) the components of the strain tensor are found as follows: The stress components are found by substituting Eq. (10) There are 2N unknowns, i.e., 1 k a and 2 ..., , , 1 1 over the boundary Γ . By using Eqs. (8) and (12) and satisfying Eqs. (6) and (7) for the M collocation points, the following system of 2M equations are obtained: The system of equations given in Eqs. (14) and (15) can be expressed as follows: Eq. (20) represents a system of 2M linear equations with 2N unknowns. For the cases with M=N, the system may be solved by a standard system solving method; however, the system should be solved using an ordinary least-squares method for the cases with For elastostatic problems with a body force or elastodynamic problems, a nonhomogeneous term appears in Eq. (1). In the MFS and other methods that use point collocation, the nonhomogeneous term can be treated using particular solutions [Liu, Wang and Qu (2018)].

A suitable configuration for source points in the MFS for 2D elasticity
The configuration of source points in the MFS has a noticeable effect on the accuracy of the solution. If source points are located very close to the boundary, the solution would experience undesired oscillations on the boundary. On the other hand, if source points are located very far from the boundary, we confront an ill-conditioned system of equations. Hematiyan et al. Khosravifard (2015, 2018)] have introduced a parameter for each source point that describes the position of the source point relative to the physical boundary and relative to its neighbouring source points. This parameter has been called the location parameter, which is defined as follows: of the location parameter of source points means that the distance between pseudo and physical boundaries is large in comparison with the distance between two adjacent source points. For a problem with a large value of the location parameter of source points, the local effect of each source (pseudo force) on the boundary reduces and therefore, the condition number of the coefficient matrix becomes large. For simple geometries such as a circular domain, with regular boundary conditions, larger values of the location parameter lead to more accurate solutions if a sufficiently accurate system solver is used [Hematiyan, Haghighi and Khosravifard (2018)]. It has been observed that a location parameter greater than 0.80 is suitable for source points in the MFS for problems with the Laplace fundamental solutions [Hematiyan, Haghighi and Khosravifard (2018)]. In this section, by a numerical study it is observed that a location parameter greater than 0.85 is suitable for source points in the 2D elastostatic problems. Suppose that we want to solve an elastostatic problem in a circular domain as shown in   The obtained solutions for the deformed shape in two different cases with location parameters of 0.2 and 0.4 for source points are shown in Fig. 3. The exact solution for r u is uniform over the boundary; however, as it can be seen in Fig. 3, the solution obtained by the MFS has a local oscillation. Also it is observed that the oscillation in the case with For further study, we consider different cases with different numbers of source points (N=16, 64, 256, 1024) and different values of the location parameter. The angle covered by two adjacent source points is denoted by θ ∆ as shown in Fig. 4. By increasing the number of source points, the value of θ ∆ decreases. For the case with N=16, the value of θ ∆ is 8 / π (Fig. 4), where the boundary between two adjacent base points has a curved shape. For the case with N=1024, the value of θ ∆ is / 512 π , which is small angle and it means the boundary between two adjacent base points in this case is     however, the present study shows that the location parameter for source points should be more than 0.85 for 2D elastostatic problems. Moreover, it should be mentioned that a configuration of source points with 85 . 0 ≥ K is a suitable configuration for 2D elastostatic problems but not necessarily the optimum one.

System solving using the singular value decomposition and a suitable truncation criterion
As it was demonstrated in the previous section, for avoiding undesired local oscillation of the solution on the boundary, it is necessary to consider a sufficient distance between pseudo and physical boundaries. This leads to an ill-conditioned system of equations in some problems [Liu (2008)]. Therefore, a reliable method should be used to solve the system of equations. The TSVD method has been widely used to solve ill-conditioned systems of equations (e.g., [Ramachandran (2002) where the integer number T N represents the number of truncated singular values. The suitable value for T N should be selected based on a suitable criterion. Some criteria based on generalized cross validation [Golub, Heath and Wahba (1979)] and L-curve method [Hansen (1992); Hansen and O'Leary (1993)] have been proposed. In this work, an alternative criterion based on the accuracy of the MFS solution on the boundary is used for this purpose. In the MFS for elastostatic problems, the governing equations of the problem are exactly satisfied. Boundary conditions at collocation points are exactly satisfied too. However, the boundary conditions are approximately satisfied at other boundary points between collocation points. We consider a few test points between two neighbouring collocation points as shown in Fig. 9. A value of the truncation parameter, i.e., T N , which minimizes the summation of the errors of boundary conditions at test points ( Γ E ) is selected. We can start with 0 = T N and then increase T N and compute the corresponding Γ E . For an ill-conditions system of equations, the value of Γ E first decreases and then increases with increase of T N and therefore the optimum value of T N can be simply found.

Reducing the condition number of the system of equations by normalizing the shear modulus
In applied engineering problems with stiff materials such as metals, the magnitude of the shear modulus is very large. In these problems, the magnitude of the stress over a part of the boundary with natural boundary condition can be very large, e.g., Pa 10 8 . On the other hand, in the same problems, the magnitude of the displacement is very small, e.g., m 10 3 − , on the boundaries with essential boundary conditions. Considering Eqs. (6) and (7), it is observed that equations corresponding to essential boundary conditions are expressed in terms of small numbers, while equations corresponding to natural boundary conditions are expressed in terms of very large numbers. It means that the numbers in some lines of the coefficient matrix are very small and in some other lines are very large. This configuration of numbers makes the coefficient matrix highly ill-conditioned. To avoid this issue, before solving the problem using the MFS, we set the shear modulus to the value of 1.0. Moreover, prescribed stresses in the natural boundary conditions are divided by the shear modulus. Then the problem is solved using the MFS. After that, the computed values for stresses should be multiplied by the shear modulus. For better explanation of this approach, we can rewrite the governing equations, i.e., Eqs. (1), (2) and (3), as follows: , which shows that prescribed stresses on the boundary should be divided by G in the equivalent problem. By this approach, the condition number of the coefficient matrix significantly reduces and the problem can be solved more accurately. The effectiveness of normalizing shear modulus is demonstrated in Section 6.

Examples
In this section, two examples with stress concentration are considered. In the first example, the stress concentration is due to the shape of the member, while in the second example the stress concentration is due to the presence of a concentrated load on the boundary. The locations of source points are determined based on the limitation proposed in this work. The system of equations in each example is solved by the TSVD using the criterion proposed in this study. In an elastostatic problem, two boundary conditions are defined for each boundary point. The boundary data computed using the MFS at the test points have a difference with the prescribed one (exact values). There is usually a difference between the boundary data computed using the MFS with the exact ones (prescribed values) at test points. The error of boundary data at test points, i.e., Γ E , is defined as follows: t , respectively. The total number of displacement and traction boundary data at test points are denoted by u N and t N , respectively. Therefore, t u N N + is twice the total number of test points. In this work, only one test point is considered between two adjacent collocation points. Each test point is located at the midpoint of two adjacent collocation points.

A plate with a circular hole
Consider an infinite plate with a circular hole of radius a under the uniform far-field tension T . Assuming the center of the hole is placed at the origin of the coordinates system and T is applied in the x-direction, the solution for stress field can be expressed as follows [Sadd (2009) To solve the problem by the MFS, two cases with 58 and 118 collocation points are considered. The configurations of source and collocation points for these two cases are shown in Figs. 11 and 12, respectively. In the first case with 58 = N , the distance between each source point to the main boundary of the problem is 0.75. In this case, the minimum and average values of the location parameter for source points are 0.85 and 0.95, respectively. In the second case with 118 = N , the distance between each source point to the main boundary of the problem is 0.65. In this case, the minimum and average values of the location parameter for source points are 0.93 and 0.98, respectively. In the two cases the minimum value of the location parameter is greater than 0.85. In the MFS analysis of the problem, the method described in Section 4 was used to solve the system of equations and it was observed that there was no need for truncating any singular value, i.e., The effects of the shear modulus normalizing (described in Section 5) on the results are also investigated in this example. The results for normal stresses at points A and B with and without shear modulus normalizing are given in Tab. 3. From this table, it is observed that the condition number significantly reduces by performing the shear modulus normalizing procedure. In the case without shear modulus normalizing, 3 singular values are truncated; however, the error of the obtained results is more than the case with shear modulus normalizing.   100.14 (6.77%) Exact ------300.00 107.41

A rectangular domain under a pressure with a sharp variation
In the previous example, the complexity in the geometry resulted in a stress concentration at a part of the domain. In the present example, a simple domain is considered; however, a load with a sharp variation causes a stress concentration at a part of the domain. The geometry and boundary conditions of the problem are shown in Fig. 14. The pressure To solve the problem by the MFS, two cases with 56 and 111 collocation points are considered. The configurations of source and collocation points for these two cases are shown in Figs. 15 and 16, respectively. In both cases, near the critical point (0,0.5), a smaller distance between collocation points has been considered. In both cases, the location parameter for all source points around the critical point is 0.95. Therefore, the source points are located closer to the boundary around the critical point.   Fig. 19. The numerical results are listed in Tab. 4. From this table it is observed that a large number of finite elements are required to accurately model the problem. However, it is possible to obtain very accurate solutions using the MFS with a relatively small number of collocation points.
In the previous analysis, we considered a configuration of collocation points with nonuniform spacing. For better clarification, we also examined a case with 116 collocation points and with uniform spacing as shown in Fig. 20. In this case, the distance between the pseudo and main boundaries is uniform and equal to 0.285. This distance has been found by solving an optimization problem for minimizing Γ E (see Eq. (36)). The obtained value for x σ at point (0,0.5) in this case is -191.8 MPa. Comparing this value with those reported in Tab. 4, it is observed that the results obtained by the configurations with non-uniform spacing of collocation points are much better than the case with uniform spacing between collocation points. Figure 19: Discretization of the rectangle using 1200 finite elements

Conclusions
In this work, by a numerical study, the effects of the distance between pseudo and main boundaries in the MFS for 2D elastostatic problems were investigated. It was observed that in the cases where the source points were near to the main boundary, the obtained solutions for displacements and stresses had unacceptable oscillations on the boundary. It was suggested to locate the source points in a way that their location parameters would be greater than 0.85. Greater values of the location parameter results in more smooth solutions; however, it makes the resulting system of equations ill-conditioned. To solve the ill-conditioned system of equations, the TSVD with a criterion based on the accuracy of imposition of boundary conditions was employed. The solutions of elastostatic problems obtained by the MFS exactly satisfy the governing equations of elasticity; however, the boundary conditions are only satisfied at collocation points. By making sure that boundary conditions are also accurately satisfied at some test points between collocation points, we can obtain very accurate solutions.
In real engineering elastostatic problems, the values of the elastic constants are very large and the values of displacements are very small. This is another reason for making the system of equations ill-conditioned. A procedure for normalizing the shear modulus was suggested in this work and it was observed that it is very effective in reducing the condition number of the system of equations and improving the results. It was found that problems involving tractions with a local concentration on the boundary can be suitably modeled and accurately solved by the MFS. In these cases, one can consider condensed collocation points near the critical point on the boundary. The spacing between collocation points should be gradually changed from a small value at the critical point to a larger value at the other parts of the boundary. We can use a specific value of the location parameter (e.g., 95 . 0 = K ) for source points around the critical point. By this approach, the source points around the critical points would be closer to the boundary (see Figs. 15 and 16). The results of the presented numerical examples showed that by considering the remarks made in this work, we are able to accurately solve very complicated 2D elastostatic problems using the MFS with a small number of collocation points. Since the MFS is a boundary-type truly meshfree method, it can be efficiently used for analysis of shape and structural optimization problems. It should be mentioned that the application of the presented method to three-dimensional problems is not straightforward and needs more investigations.