An enriched meshless finite volume method for the modeling of material discontinuity problems in 2D elasticity

A 2D formulation for incorporating material discontinuities into the meshless finite volume method is proposed. In the proposed formulation, the moving least squares approximation space is enriched by local continuous functions that contain discontinuity in the first derivative at the location of the material interfaces. The formulation utilizes space-filling Voronoishaped finite volumes in order to more intelligently model irregular geometries. Numerical experiments for elastostatic problems in heterogeneous media are presented. The results are compared with the corresponding solutions obtained using the standard meshless finite volume method and element free Galerkin method in order to highlight the improvements achieved by the proposed formulation. It is demonstrated that the enriched meshless finite volume method could alleviate the expecting oscillations in derivative fields around the material discontinuities. The results have revealed the potential of the proposed method in studying the mechanics of heterogeneous media with complex micro-structures.


INTRODUCTION
The finite-volume method FVM is a numerical technique that provides approximate solutions to boundary value problems with partial differential equations. In this method, the computational domain is divided into smaller sub-domains, called as finite volumes, each containing one node point on the discretized geometry and the partial differential equation is satisfied in the integral sense over each finite volume. A distinctive feature of the method is the use of boundary integral instead of the domain integral, for satisfying partial differential equation.
The mesh-based formulation of FVM is a well-established technique in a wide range of problems concerned with the mechanics of solids and structures, see Taylor et al. 2003, Fallah 2004, Cardiff et al. 2014, Escarpini Filho and Marques 2016, Cavalcante and Pindera 2016 . However, for certain classes of problems, the heavy and rigid reliance of the mesh-based FVM on a mesh leads to some difficulties. Some of these difficulties are as follows Liu 2009 : • High computational cost and human-labor in creation of a quality mesh for domain with complex geometry, • Additional computation as well as a degradation of accuracy in remeshing approach to deal with evolution problems, and • Discontinuous nature in calculating stress fields, due to the element-wise continuous nature of the displacement field assumed in the meshbased formulation.
These difficulties could be decreased by using a mesh independent formulation. For this purpose, meshless formulation of FVM, which is also named as MLPG5 by Atluri and Shen 2002 , has been introduced. This method has been accepted by many researchers and widely used for solving problems in solid and structural mechanics, see Batra et al. 2004, Han et al. 2005, Sladek et al. 2008, Hosseini et al. 2011, Soares et al. 2012, Ebrahimnejad et al. 2014, 2017 Interpolation of field variables by moving least squares approximation and integrating weak form of the equilibrium equation over local finite volumes in meshless FVM eliminates the need for mesh generation. Instead, it relies on background local finite volumes where each finite volume is independent of the others, i.e. they can be of any geometric shape and size and even using overlapped finite volumes are allowed. However, they should locate entirely inside the boundaries and their union should cover the whole domain. In addition, the high order of continuity inherits from moving least squares approximation, provides solutions with smooth derivatives throughout the domain. As a result, negligible pre-processing and post-processing costs of meshless FVM are the main attractive features over the prevalent mesh-based methods Shen 2002 . Atluri andShen 2002 have demonstrated that meshless FVM can be a simple and efficient alternative to the all-purpose FEM, due to its speed, accuracy and response stability.
In this paper the meshless FVM is applied to the material discontinuity problem which is an important issue in numerical modeling of heterogeneous media. It should be noted that heterogeneity is a ubiquitous property of solids and structures associated with features such as holes or inclusions that have an important role in the failure Gdoutos 2005 .
The solution of problems in heterogeneous media typically involves discontinuity in gradient fields. In order to capture non-smooth property of the solutions, there are two essentially different techniques. The first technique applies smooth polynomial approximation spaces and trusts in domain discretization that conforms to material discontinuities. The polynomial approximation spaces can be constructed by using either a mesh-based shape function or meshless approximation function. The alternative technique is based on 'enrichment' of smooth approximation space such that the domain discretization being independent of the material discontinuities. The enrichment will be accomplished by adding special weak discontinuous functions to the approximation space with unknown parameters that control the strength of the discontinuity.
The first above mentioned technique has been used by Cavalcante et al. 2007 within the mesh-based FVM framework for functionally graded materials and also by Li et al. 2003 in the 3D formulation of the MLPG method to problems with singularities and material discontinuities. In the latter, the heterogeneous medium is separated into isolated, homogeneous parts and then continuity constraints are enforced at the interface to 'reconnect' the pieces. To accomplish the 'separation', meshless approximation function applied separately within each homogeneous part and therefore domains of influence are truncated at the material interfaces. Also, local subdomains, which partially cut by the interface, are redefined to locate entirely inside the boundaries. In spite of the simple concept of this technique, the presence of truncated domains of influence could introduce some mesh dependency into the solution Cordes andMoran 1996, Li et al. 2000 ; Furthermore, re-definition of local subdomains and domains of influence could be difficult and computationally expensive de Borst 2008 . The socalled enrichment technique has been applied both in the area of mesh-based and meshless methods. Melenk and Babuška 1996 enriched the standard finite element approximation, Belytschko and Black 1999 set up the frame of the extended FEM XFEM and Belytschko et al. 1996 implemented it within the EFG method for introducing discontinuous derivatives into the solution. Batra et al. 2004 utilized this technique in the 1D meshless formulation of FVM to analyze heat conduction in a bimetallic circular disk. Recently, Yoon andSong 2014 andHu et al. 2015 applied this technique in order to handle discontinuities in heterogeneous media. Using the enrichment technique is interesting since a fix and regular domain discretization can be applied and domain discretization cost is reduced to a minimum An et al. 2011 . In this paper, based on the enrichment technique used already by Krongauz and Belytschko 1998 , a 2D formulation within the meshless FVM framework is proposed to handle material discontinuity.
The structure of the present paper is as follows. In the next section, the basic concepts including standard meshless FVM, moving least squares approximation and enrichment of approximation space are briefly described. A 2D formulation in modeling heterogeneous material within the meshless FVM framework is presented in section 3. Numerical examples and results to verify the capability of the formulation are presented in section 4 and concluding remarks are given in the last section.

Standard meshless FVM
Most of the meshless methods, e.g. element free Galerkin EFG and smoothed particle hydrodynamics SPH , are based on global weak forms and required the use of background mesh in order to integrate the weak form of the equilibrium equation. On the other hand, the meshless local Petrov-Galerkin MLPG method satisfies the weak form over the local sub-domains. As a special case, MLPG5 Atluri and Shen 2002 chose the Heaviside step function as the test function. Consequently, domain integrals are vanished and integration is performed along the boundary of local sub-domains. Satisfying local weak form denotes the momentum balance law over the local subdomains that resemble the concept of FVM. Therefore, the MLPG5 method could be recognized as the meshless FVM Atluri et al. 2004  where σ is the stress tensor and b is the body force measured per unit volume. The essential and natural boundary conditions are as follows where u and t are defined as the prescribed displacements and tractions, respectively. u  is the global essential boundary and t  is the global natural boundary. Utilizing weighted residual method, weak form of the Equations 1 -2 over I-th local finite volume s  enclosed by s  , can be written as where I v is the test function for I-th local finite volume and 1  >> is a penalty parameter which is applied to impose the essential boundary conditions. In meshless formulations, a local approximation scheme will be applied to approximate the trial function; therefore, the essential boundary conditions cannot be imposed directly, a priori, but superimposed in a weak form. In this paper, the penalty method is applied to superimpose the essential boundary conditions, as the second term in Equation 3 where n is the unit outward normal vector to the boundary s  . By using the Heaviside step function as the test function, imposing natural boundary conditions and doing some algebraic operations, one can obtain the following expression

Moving least squares approximation
In Meshless methods, a local partition of unity approximation is applied to approximate unknown field functions with the nodal parameters of unknown variable at some randomly scattered nodes. Among available local partition of unity approximation schemes, moving least squares MLS approximation has been utilized in meshless FVM due to its reasonable accuracy and simplicity of extension to 3D problems Atluri and Zhu 1998 . Let the neighborhood of a point x, where   In the MLS approximation, the coefficient vector   a x can be determined by minimizing the weighted re- Where n is the number of nodes in the domain of definition for which the weight function   The unknown   a x coefficients can be found by minimizing where the matrices A and V are defined as Solving for   a x from Equation 13 and substituting in Equation 6 , the MLS approximation can be defined as where the shape function of the MLS approximation corresponding to I-th node defined by The It should be noted that according to Equation 17 , continuity of the MLS approximation depends on the continuity of the basis function and also the smoothness of the weight function. Thus, the weight function has a significant role in the efficiency of the MLS approximation. The exponential function, cubic and quartic spline functions are common useful examples of weight functions Liu and Gu 2005 . In this paper, the cubic spline weight function is used where defined as where I I c r r   I x x  and  I x x denotes the distance from I-th node to point x and I c r is the radius of the circular support for the weight function for I-th node. The cubic spline weight function has 2 nd order continuity over the entire domain; therefore, MLS approximant   h u x is C 2 continuous over the entire domain.

Enrichment of approximation space
According to the 'reproduction' property of the MLS approximation, it can reproduce any functions that are included in the basis. Clearly Equation 16 is just able to reproduce continuous fields; therefore, only if the basis enriched by a special approximation function, including discontinuity in the derivative, MLS approximation will reflect the discontinuity in the gradient fields Liu and Gu 2005 . The enrichment can be accomplished by adding a special function to the approximation space, extrinsically. In addition, most non-smooth phenomena in solid mechanics have a local character; therefore, it may be useful to employ the enrichment in local regions Fries and Belytschko 2010 . There are two conditions for the special approximation function which is customized for material interface. It should be continuous in the approximation field with a discontinuity in the first derivative at the location of the material interfaces Belytschko et al. 1996 . In addition, it should have compact support to maintain discrete equations banded and sparse Krongauz and Belytschko 1998 . There are several ways to construct such an enrichment approximation function. An et al. 2011 introduced a piecewise linear function with a cusp at the location of material interface. Krongauz and Belytschko 1998 employed a higher-order polynomial function with more constraints and also a localized ramp function. More different enrichment approximation functions could be defined by using the level set function which used within the XFEM framework, see Sukumar et al. 2001, Moës et al. 2003 . In this paper, the spline function introduced already by Krongauz and Belytschko 1998 is applied to enrich MLS approximation space. In a two dimensional domain with d n lines of material discontinuity, the enrichment approximation function for J-th line of material discontinuity is defined as:  In 2D problems, q is expressed in term of the arc length of the material discontinuity and then is discretized over the discontinuity line as follows where K  are one dimensional approximation function, 1 n is the number of enrichment nodes in the support domain of 1D approximation and ˆK q are fictitious nodal values that are considered as unknowns in the discretized equations. In contrast to field nodes, which are employed to define fictitious nodal values of displacement fields, enrichment nodes are distributed along material discontinuities in order to discretize the amplitude parameters; they have no displacement degree of freedom. The definition of r and s in a typical two-dimensional domain is illustrated in Figure 2. It should be noted that, at the points far away from discontinuities the second term on the right-hand side of Equation 20 vanishes due to the compact support property of enrichment approximation function and this equation degenerates to Equation 16 . As a result, using enrichment approximation function, Equation 20 led to a local enriched formulation.

Domain discretization
A promised feature of enriched formulations is the domain discretization independent of the geometry and location of the heterogeneities Sukumar et al. 2001 ; therefore a regular discretization can be employed in heterogeneous media, but these regular sub-domains may be intersected by irregular-shaped external boundaries. In the proposed formulation, the Voronoi tessellation is applied for constructing local background finite volumes for a set of considered field nodes. The resulted finite volumes can be regular inside the domain and irregular for nodes next to the external boundaries.
It should be mentioned that in mesh-based methods usually mesh refinement is applied where there are details. In work presented by Cavalcante et al. 2007 , the need for extensive mesh refinement at the arbitrarily shaped geometries has been eliminated using a parametric mapping. Li et al. 2003 , which applied pure spherical local finite volumes in their meshless FVM formulation, reduced the size of the finite volumes that are close to the weak or strong boundaries so as not to cross over the boundaries; therefore, a dense cloud of nodes is required to completely represent the irregularities.
As mentioned above, the inherent flexibility of the meshless FVM enables the use of the Voronoi tessellation for the domain discretization. In mathematics, for a given set of nodes, partitioning of a domain into sub-domains based on distance to the considered nodes called Voronoi tessellation Okabe et al. 2009 . These sub-domains, so called as finite volumes are space-filling polygons in 2D space and can be constructed corresponding to each node by forming the perpendicular bisector lines between the nodes and connecting the lines around each node. A Voronoi tessellation and the local background finite volume associated with an arbitrary node are illustrated in Fig  In heterogeneous media that contains material interfaces, enriched approximation space can be applied in order to make domain discretization independent of the material discontinuities; to this end, the enriched approximation Equation 20 is utilized as displacement trial function. In order to obtain the discrete equations, approximated displacement function and its derivatives are needed. The enriched approximation is given as By substituting Equations 22 -27 into Equation 5 , finally one obtains the linear algebraic discretized equations for a local finite volume as follows where Û is the vector of unknowns, K and F are the equivalent local stiffness matrix and load vector, respectively, and can be written as follows where t , b and u are the prescribed tractions, the body force measured per unit volume and the prescribed displacements, respectively. In Equation 29 , D is the elasticity matrix and  n contains the entries of the unit outward normal vector and defined as Accordingly, the vector of unknowns, Û , is expressed as It should be noted that both N and B are composed of the standard part and the enriched part, where the standard part is computed at every node but the enriched part is only computed at enriched nodes. Due to the non-polynomial form of the both MLS approximation and enrichment approximation function, exact integration of the weak form of the balance law is difficult. So, efficient numerical integration of the weak form may have an important role in the accuracy and computational efficiency of the method Nguyen et al. 2008;Ventura et al. 2009 . Hence, weak form of the balance law can be numerically integrated by Gaussian quadrature rule over the finite volumes boundaries. Ventura et al. 2009 have demonstrated that the boundary integration scheme is much more accurate than the domain integration for the same computational cost and it is especially efficient when the integrand has discontinuous and/or non-polynomial form.

Treatment of additional unknowns
Although in enriched Galerkin-based methods assembling all local equations to the global system leads to as many discrete equations as the number of unknowns, but; the proposed formulation faces to a different situation.
Consider a discretized 2D problem with 1 n enrichment nodes and 2 n field nodes. There are two unknown fictitious nodal values associated with every field node and every enrichment node. The unknowns associated with the field nodes are fictitious displacement fields and with the enrichment nodes are fictitious nodal discontinuities amplitude parameters. Therefore, the vector of unknowns has   2 n discrete equations. This is due to the stack form assembly procedure in finite volume methods.
In order to have the required   1 2 n discrete equations, we need 1 n additional local finite volumes. Thanks to the inherent flexibilities of the meshless FVM, arbitrarily overlapping local finite volumes can be defined, which resolves this difficulty. To this end, we define 1 n new local finite volumes corresponding to the selected field nodes. By doing so, Equation 28 provides two additional discrete equations for each newly-defined local finite volume. These equations are independent from the previous set one.

NUMERICAL EXPERIMENTS AND RESULTS
In order to assess the performance and accuracy of the proposed formulation, a set of numerical experiments are performed. The considered test problems are linear elastostatic within the range of infinitesimal deformations. A one-dimensional bi-material bar and also a two-dimensional infinite plate with a central inclusion both subjected to the prescribed displacement fields are studied. The results are compared with the corresponding solutions obtained by the standard meshless FVM in order to demonstrate the improvements achieved by the proposed formulation. In these experiments, both numerical methods utilized similar domain discretization, but different local approximations. The meshless FVM takes Equation 16 , whereas the proposed formulation takes Equation 20 as displacement trial function. The computations are performed using a domain discretization which is independent of the shape and location of the interface with the appropriate material property at each quadrature point. In order to further highlight the capability of the proposed formulation, some comparisons are also made with the results from the EFG method which has been introduced by Belytschko et al. 1994 .

1D Bi-material bar
Consider a one-dimensional bar of length L and constant cross sectional area, A, composed of two different materials, i.e. tively. There is a material discontinuity at x  . The bar is fixed at one end and the other end is subjected to a prescribed displacement u . The problem is depicted in Figure 4. 10/22   , with a material discontinuity at x  . The bar is fixed at one end and the other end is subjected to a prescribed displacement u .
In the computations, let L 10, A 1,  5, 1 E 1, 2 E 0.5 and u 1. The computation is performed using 22 uniformly distributed field nodes. One-dimensional finite volumes are constructed with faces located at the middle of the field nodes. In the proposed formulation one enrichment node at the location of the material discontinuity is required. The domain discretization and corresponding finite volumes used in this test are shown in Figure 5. It should be noted that the computations of the EFG method are also performed by using the similar nodal arrangement applied in the standard meshless FVM. A linear basis and the cubic spline weight function are used in the MLS approximation. The essential boundary conditions are imposed by the penalty method using a value of 1.00E 07 as penalty parameter. The exact analytical solutions for this problem are as follows Krongauz and Belytschko 1998 The nodal displacement and stress values are obtained from the proposed formulation and compared with the exact analytical solutions in Figures 6-7, respectively. In addition, the numerical solutions obtained from the meshless FVM and EFG method are also shown in the same figures for comparison purposes. Figure 6 shows a very good agreement between the nodal displacements obtained from the proposed formulation and the exact analytical solution. In Figure 7, all these numerical methods result in a stress concentration in the vicinity of the material discontinuity. In order to compare local accuracy of these methods, a local relative error norm is defined as the maximum value of   2 / numerical exact exact     which is calculated at each node. According to Figure 7, the EFG method yields a local error of 4.95%, the meshless FVM produces the local error 47.53% where the proposed formulation improved it to 17.24%.  In order to further compare the proposed formulation with the meshless FVM, a convergence study is conducted using 22, 34 and 50 uniformly distributed field nodes. As mentioned above, in these tests, the proposed formulation employed one enrichment node at the location of the material discontinuity. The results of the numerical nodal stresses are obtained and a magnified view around the interface is presented in Figure 8. It is clear that the maximum error in nodal values has remained roughly at the same level while the node number increases; therefore, the accuracy of the proposed formulation can't be achieved by the meshless FVM even with a dense nodal distribution. where u and  are the displacement and stress at each node, respectively, and  is the computational domain of the bar. Table 1 represents the numerical results for the relative error norms obtained with 50 uniform field nodes and one enrichment node. It is clear that both of the proposed formulation and the meshless FVM give almost equal relative error norms in displacement. On the other hand, the proposed formulation yields to a more accurate solution for stress; the relative error norm in stress obtained from the proposed formulation is about half of the meshless FVM result. In addition, both error norms obtained from the EFG method are significantly smaller than those of other methods.
Although in this numerical experiment, the EFG method demonstrates good results, however, it should be mentioned that the observed performance of EFG method is due to the well-adjusted quadrature points which cannot be accommodated in 1D formulation of FV-based methods. It is worth mentioning that in the 1D formulation of FV-based methods quadrature points are restricted to be located at two endpoints of each 1D local finite volume. However, as demonstrated, the effect of this drawback has been reduced by the proposed formulation through the enrichment of approximation space.

2D infinite plate with an inclusion
Consider a two-dimensional infinite plate with a circular inclusion subjected to a prescribed displacement field. In the numerical model, a bi-material square domain of length L with a centered circular inclusion of radius 1 r is considered, i.e. Ω Ω Ω   , with a material discontinuity across the interface 1 Γ .
The model is subjected to a prescribed displacement field on the external boundary 2 Γ .
In the polar coordinate system, the essential boundary conditions on the external boundary 2  are given by , 0 r u r u q = = .

41
In the computations, L 2, 1 r 0.4, 1 E 2, 1  0.2, 2 E 10, 2  0.3 and the plane strain state are assumed. The computation is made using 452 field nodes. The proposed formulation added 36 enrichment nodes along the line of material discontinuity. The Voronoi tessellation is employed to construct finite volumes. The domain discretization and corresponding finite volumes used in this test are shown in Figure 10. It should be noted that similar nodal arrangement used in the standard meshless FVM is also utilized in the EFG method. A four-point Gaussian quadrature rule is applied in order to integrate the weak form of governing equation. A quadratic basis, circular domain of influence and the cubic spline weight function are used in the MLS approximation. The essential boundary conditions are imposed by the penalty method using a value of 1.00E 07 as penalty parameter. a b Figure 10: 2D infinite plate with an inclusion. Domain discretization and corresponding finite volumes. Blue and red dots are denoted the enrichment nodes and field nodes, respectively. 1D and 2D local background finite volumes are illustrated by blue and red lines, respectively. Due to symmetry, only a quarter of the domain is shown. a Meshless FVM, b Proposed formulation.
In the polar coordinate system, the exact analytical solution for displacements is as follows Sukumar et al. 2001 2 2 1 2 2 1 1 2 2 where  and  are appropriate Lame constants at each material and  is defined as The analytical solutions for this problem are given in the polar coordinate system, in which they can be transformed to the Cartesian coordinate system using equations provided in Appendix A.
The solutions obtained from the proposed formulation, meshless FVM, EFG method and also exact analytical one are shown in Figures 11-12 for the horizontal displacement x u and the normal stress xx  both along the line y x  , respectively. A magnified view around the discontinuities is presented in Figure 12 for better observation. It is clear from Figure 11 that the nodal displacement curve obtained from the proposed formulation agrees well with the exact analytical solution, especially in the vicinity of the material discontinuity. In addition, it is seen in Figure 12 that a significantly higher accuracy in the nodal stress curve is achieved by the proposed formulation as compared with the other meshless methods. The meshless FVM and EFG method fail to reproduce the expected jumps in the nodal stress curve. However, oscillations around the material discontinuity are vanished in the solution obtained from the proposed formulation.
Another notable observation is that normal stress xx  is constant within the inclusion. The L2 relative error norm of the predicted xx  inside the inclusion for meshless FVM and EFG method is 0.012 and 0.008, respectively, while the proposed formulation reduces it to 0.002.
In the calculation of error, the numerical solutions are compared to the exact analytical ones. Let   num and   exact represent approximated numerical and exact analytical solutions, respectively, then the relative error norms in displacement and stress could be defined as follows where u and σ are the vectors of displacement and stress at each point, respectively, and  is the computational domain of the problem. Table 2 represents the results for the relative error norms. Clearly, proposed formulation yields to more accurate results than both meshless FVM and EFG method. In order to further evaluate the proposed formulation, a convergence study is conducted using 452, 1068 and 1548 field nodes. In the above mentioned tests, the proposed formulation employs 36, 84 and 140 enrichment nodes along the material interface, respectively. Computations of the meshless FVM are carried out using similar field nodes. The numerical results for the relative error norm in displacement predictions are plotted against the field node number in Figure 13. The proposed formulation yields to solutions with less error than meshless FVM and the error of the proposed formulation falls below 0.2% when the node number increases to 1688 including 1548 field nodes and 140 enrichment nodes .

2D infinite plate with an extremely soft inclusion
In order to demonstrate the capability of the proposed formulation in the extreme cases, above mentioned two-dimensional infinite plate is considered in a different state regarding to inclusion material property. In this case we let 1 E 0.02 while 2 E remains equal to 10. Similar geometry, boundary conditions, nodal arrangement and numerical techniques are used as presented in the previous test.
The solutions for the horizontal displacement x u and the normal stress xx σ both along the line y x  , obtained from the proposed formulation, meshless FVM, EFG method and also exact analytical one are shown in Figures 14-15, respectively.   Figure 14 shows a very good agreement between the solution of the proposed formulation and the exact analytical one. The kinks across the material discontinuities are accurately captured. In Figure 15, a magnified view around the discontinuities is presented. As can be seen in Figure 15, while nodal stress curves obtained from meshless FVM and EFG method show oscillations in the vicinity of the material discontinuity, they are alleviated in the proposed formulation results. As mentioned already, the normal stress xx σ is constant within the inclusion. The L2 relative error norm of the predicted xx σ inside the extremely soft inclusion for the meshless FVM, EFG method and proposed formulation are 0.019, 0.024 and 0.001, respectively.
In addition, the results for the relative error norms calculated using Equations 47 -48 are reported in Table 3. Again, the convergence study is conducted where three different tests with the same nodal arrangement as the problem of 2D infinite plate with an inclusion are considered. The numerical results for the relative error norm in displacement predictions are plotted against the field node number in Figure 16. It is clear that the proposed formulation yields to more accurate results as compared to the meshless FVM.

CONCLUSIONS
The aim of this paper is to propose a 2D formulation based on the meshless FVM framework for solving problems with material discontinuity. In order to capture the discontinuity in the gradient fields, the MLS approximation space has been enriched by a continuous function, having discontinuity in its derivatives. In contrast to the standard meshless methods, in the proposed formulation, the domain discretization task is done independently of the material discontinuities and also displacement and traction continuity constraints at the interface is automatically satisfied without additional treatment. Utilizing enrichment and also space-filling Voronoi-shaped local finite volumes, greatly simplifies the domain discretization task.
Several numerical experiments have been performed. The comparisons with the analytical solutions have shown an excellent agreement in both displacement and stress fields predictions. It has demonstrated that the proposed formulation could alleviate the expecting oscillations in derivative fields in the vicinity of the material discontinuities appeared in the standard meshless methods. The formulation has been also successfully applied to the case of extremely soft inclusion problem. In comparison with the Galerkin-based methods, the proposed formulation is particularly attractive as it employs boundary integration over the finite volumes boundaries instead of domain integration. It should be mentioned that, as stated in Ventura et al. 2009 , the boundary integration scheme is much more accurate than the domain integration for the same computational cost and it is especially efficient when the integrand has non-polynomial form.
The results indicate the potential and promise feathers of the proposed formulation for studies of the mechanics of heterogeneous media. Two-dimensional displacement and stress components in the polar coordinate system can be transformed to the Cartesian coordinate system through the equations provided here. The coordinate systems are illustrated in Figure A.1.