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

Abstract

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 Voronoi-shaped 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.

Keywords
finite volume method; meshless method; material discontinuity; enrichment technique; Voronoi tessellation

1 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 mesh-based 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 , 2015 , 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 ( Atluri and Shen 2002 ). Atluri and Shen (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 sub-domains, 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 and Moran 1996 , Li et al. 2000 ); Furthermore, re- definition of local sub-domains and domains of influence could be difficult and computationally expensive ( de Borst 2008 ). The so-called 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 and Song (2014) and Hu 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.

2 BASIC CONCEPTS

2.1 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 sub-domains that resemble the concept of FVM. Therefore, the MLPG5 method could be recognized as the meshless FVM ( Atluri et al. 2004 ).

Let $\Omega$ be a domain bounded by $\partial \Omega$ . Momentum balance principles for static problems are given by

$${\sigma }_{ij,j}+b_i=0forallx\in \Omega ,$$ (1)

where $\sigma$ is the stress tensor and $b$ is the body force measured per unit volume. The essential and natural boundary conditions are as follows

$u_i={\overline{u}}_{i}on{\Gamma }_u$ , (2a)

$t_i={\overline{t}}_{i}on{\Gamma }_t$ , (2b)

where $\overline{u}$ and $\overline{t}$ are defined as the prescribed displacements and tractions, respectively. ${\Gamma }_u$ is the global essential boundary and ${\Gamma }_t$ is the global natural boundary. Utilizing weighted residual method, weak form of the Equations (1) - (2) over I-th local finite volume ${\Omega }_s$ enclosed by $\partial {\Omega }_s$ , can be written as

$${\displaystyle {\int }_{{\Omega }_s}^{}\left({\sigma }_{ij,j}+b_i\right)v^{I}d\Omega }-\alpha {\displaystyle {\int }_{{\Gamma }_{su}}^{}\left(u_i-{\overline{u}}_i\right)v^{I}d\Gamma }=\mathrm{0,}$$ (3)

where $v^I$ is the test function for I-th local finite volume and $\alpha >>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) . Let ${\Gamma }_s$ denotes a part of the $\partial {\Omega }_s$ located on the global boundary, then ${\Gamma }_{su}={\Gamma }_s\cap {\Gamma }_u$ is a part of the ${\Gamma }_s$ that coincides with the global essential boundary. By applying the divergence theorem, Equation (3) could be rewritten as

$${\displaystyle {\int }_{\partial {\Omega }_s}^{}{\sigma }_{ij}{n}_j{v}^{I}d\Gamma }-{\displaystyle {\int }_{{\Omega }_s}^{}\left({\sigma }_{ij}{v}_{,j}^I-b_i{v}^I\right)d\Omega }-\alpha {\displaystyle {\int }_{{\Gamma }_{su}}^{}\left(u_i-{\overline{u}}_i\right)v^{I}d\Gamma }=0$$ , (4)

where $n$ is the unit outward normal vector to the boundary $\partial {\Omega }_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

$-{\displaystyle {\int }_{{\Gamma }_{s0}}^{}{t}_{i}d\Gamma }-{\displaystyle {\int }_{{\Gamma }_{su}}^{}{t}_{i}d\Gamma }-\alpha {\displaystyle {\int }_{{\Gamma }_{su}}^{}{u}_{i}d\Gamma }={\displaystyle {\int }_{{\Gamma }_{st}}^{}{\overline{t}}_{i}d\Gamma }+{\displaystyle {\int }_{{\Omega }_s}^{}{b}_{i}d\Omega }-\alpha {\displaystyle {\int }_{{\Gamma }_{su}}^{}{\overline{u}}_{i}d\Gamma }$ , (5)

where ${\Gamma }_{s0}$ is a part of the $\partial {\Omega }_s$ which is totally inside the global domain and ${\Gamma }_{st}={\Gamma }_s\cap {\Gamma }_t$ is a part of the ${\Gamma }_s$ that coincides with the global natural boundary. Equation (5) represents a weak form of the balance law in the local finite volume ${\Omega }_s$ , with the boundary conditions being enforced. It can be seen that the left hand side of Equation (5) , which leads to the stiffness matrix, does not involve domain integration. In addition, in the absence of body force, domain integration is totally eliminated. It should be mentioned that Equations (1) - (2) would be satisfied, a posteriori, in the global domain $\Omega$ and on its boundary $\partial \Omega$ , respectively, if the union of all local finite volumes covers the global domain ( Atluri and Zhu 1998 ).

2.2 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 $x={\left[x_1,x_2\right]}^T$ contains space coordinates, be a sub-domain ${\Omega }_x$ and called the domain of definition of MLS approximation of the field function $u\left(x\right)$ . The MLS approximation $u^h\left(x\right)$ of the function $u\left(x\right)$ can be defined over a given set of nodes $\left\{x^1,x^2,\dots ,x^n\right\}$ as

$$u^h\left(x\right)={\displaystyle {\sum }_{J=1}^m{p}^J\left(x\right)a^J\left(x\right)}=p^T\left(x\right)a\left(x\right),\forall x\in {\Omega }_x$$ , (6)

where $p\left(x\right)$ is the complete polynomial basis function of order $m$ and $a\left(x\right)$ is a vector containing associated unknown coefficients. In two dimensions, $p\left(x\right)$ can be expressed as

$p^T(x)=[\mathrm{1,}{x}_1,x_2]linearbasis,m=3$ , (7a)

$p^T\left(x\right)=\left[\mathrm{1,}{x}_1,x_2,{\left(x_1\right)}^2,x_1{x}_2,{\left(x_2\right)}^2\right]quadraticbasis,m=6$ . (7b)

In the MLS approximation, the coefficient vector $a\left(x\right)$ can be determined by minimizing the weighted residual $J\left(x\right)$

$$J\left(x\right)={\displaystyle {\sum }_{I=1}^{n}w\left(x-x^I\right){\left[p^T\left(x^I\right)a\left(x\right)-{\widehat{u}}^I\right]}^2}$$ , (8)

Where $n$ is the number of nodes in the domain of definition for which the weight function $w\left(x-x^I\right)$ associated with the I-th node is positive and ${\widehat{u}}^I$ refers to the fictitious nodal values and not the nodal values of MLS approximant $u^h\left(x\right)$ at the I-th node. Equation (8) can be rewritten in the form

$$J\left(x\right)={\left[P.a\left(x\right)-\widehat{U}\right]}^T.W.\left[P.a\left(x\right)-\widehat{U}\right]$$ , (9)

where

$P={\left[p^T\left(x^1\right),p^T\left(x^2\right),\dots ,p^T\left(x^n\right)\right]}^T$ , (10)

$W=\left[\begin{array}{ccc}w\left(x-x^1\right)& \dots & 0\\ \dots & \dots & \dots \\ 0& \dots & w\left(x-x^n\right)\end{array}\right]$ , (11)

$\widehat{U}={\left[{\widehat{u}}^1,{\widehat{u}}^2,\dots ,{\widehat{u}}^n\right]}^T$ . (12)

The unknown $a\left(x\right)$ coefficients can be found by minimizing of $J\left(x\right)$ with respect to $a\left(x\right)$

$$\frac{\partial J\left(x\right)}{\partial a\left(x\right)}=A\left(x\right)a\left(x\right)-V\left(x\right)\widehat{U}=0$$ , (13)

where the matrices A and V are defined as

$A\left(x\right)={\displaystyle {\sum }_{I=1}^{n}w\left(x-x^I\right)p\left(x^I\right)p^T\left(x^I\right)}=P^{T}WP$ , (14)

$V\left(x\right)=\left[w\left(x-x^1\right)p\left(x^1\right),w\left(x-x^2\right)p\left(x^2\right),\dots ,w\left(x-x^n\right)p\left(x^n\right)\right]=P^{T}W$ . (15)

Solving for $a\left(x\right)$ from Equation (13) and substituting in Equation (6) , the MLS approximation can be defined as

$u^h\left(x\right)={\displaystyle {\sum }_{I=1}^n{\varphi }^I\left(x\right){\widehat{u}}^I}$ , (16)

where the shape function of the MLS approximation corresponding to I-th node defined by

${\varphi }^I\left(x\right)={\displaystyle {\sum }_{J=1}^m{p}^J\left(x\right){\left[A^{-1}\left(x\right)V\left(x\right)\right]}^{JI}}$ . (17)

The MLS approximation would be well-defined if matrix $A$ in Equation (14) be invertible, therefore, domains of definition should be large enough to cover at least $m$ nodes (i.e. $n\ge m$ ) for each point of interest to prevent the singularity of the matrix $A$ . These nodes are also not allowed to lie on a straight line. On the other hand, domains of definition should be small enough in order to preserve the local character of the MLS approximation ( Atluri and Zhu 1998 ).

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

$w\left(x-x^I\right)=\{\begin{array}{ll}2/3-4{\left({\overline{r}}^I\right)}^2+4{\left({\overline{r}}^I\right)}^3\hfill & {\overline{r}}^I\le 0.5\hfill \\ 4/3-4{\overline{r}}^I+4{\left({\overline{r}}^I\right)}^2-4/3{\left({\overline{r}}^I\right)}^3\hfill & 0.5<{\overline{r}}^I<1\hfill \\ 0\hfill & {\overline{r}}^I\ge 1\hfill \end{array}$ , (18)

where ${\overline{r}}^I=\left|x-x^I\right|∕r_c^I$ and $\left|x-x^I\right|$ denotes the distance from I-th node to point $x$ and $r_c^I$ 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 $u^h\left(x\right)$ is C² continuous over the entire domain.

2.3 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 $n_d$ lines of material discontinuity, the enrichment approximation function for J -th line of material discontinuity is defined as:

$${}_{}{}^J\psi \left({\overline{r}}^J\right)=\{\begin{array}{ll}-\frac{1}{6}{\left({\overline{r}}^J\right)}^3+\frac{1}{2}{\left({\overline{r}}^J\right)}^2-\frac{1}{2}{\overline{r}}^J+\frac{1}{6}\hfill & {\overline{r}}^J\le 1\hfill \\ 0\hfill & {\overline{r}}^J\ge 1\hfill \end{array}$$ , (19)

where ${\overline{r}}^J=r^J∕r_c^J$ and $r^J$ is the distance to the closest point on the J-th line of material discontinuity and $r_c^J$ is the length of the support of ${}_{}{}^J\Psi$ , i.e. the distance over which ${}_{}{}^J\Psi \ne 0$ . The plot of the enrichment approximation function ( Equation 19 ) and its derivative in one dimension are illustrated in Figure 1 , where the kink of the enrichment function and the jump in its derivative across the material discontinuity are clearly seen.

Figure 1
Enrichment function and its derivative in one dimension.

The enriched approximation $u^h\left(x\right)$ for the function $u\left(x\right)$ , is then given by

$$u^h\left(x\right)={\displaystyle {\sum }_{I=1}^{n_2}{\varphi }^I\left(x\right){\widehat{u}}^I}+{\displaystyle {\sum }_{J=1}^{n_d}{}_{}{}^{J}q{}_{}{}^J\psi \left({\overline{r}}^J\right)}$$ , (20)

where the first term on the right-hand side of this equation represents the standard meshless FVM approximation, see Equation (16) , $n_2$ is the number of field nodes in the support domain of 2D MLS approximation, $n_d$ is the number of lines of material discontinuity and ${}_{}{}^{J}q$ are amplitude parameters that govern the strength of the discontinuity across the material interfaces.

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

$q\left(s\right)={\displaystyle {\sum }_{K=1}^{n_1}{\varphi }^K\left(s\right){\widehat{q}}^K}$ , (21)

where ${\varphi }^K$ are one dimensional approximation function, $n_1$ is the number of enrichment nodes in the support domain of 1D approximation and ${\widehat{q}}^K$ 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.

Figure 2
Illustration of r and s values for a material discontinuity. Line of material discontinuity is illustrated by a solid line; solid and open circles are denoted the field nodes and the enrichment nodes, respectively and squares are arbitrary points. The values of r and s for point K are shown by dashed lines and those for point L are illustrated by dash-dot lines.

3 THE PROPOSED FORMULATION BASED ON THE MESHLESS FVM FRAMEWORK

3.1 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 Figure 3 .

Figure 3
Voronoi shaped finite volumes. Voronoi tessellation for a given set of nodes (dots) is illustrated by solid lines and the bold polygon is represented the local background finite volume associated with an arbitrary node.

3.2 Enrichment of approximation space

In order to discretize the governing equations, formulation of the meshless FVM can be used. It satisfies weak form of the balance law over local finite volumes which are formed around the field nodes, see Equation (5) . For a linear elastic solid within the range of infinitesimal deformations, the traction in Equation (5) can be written as follows

$t_i={\sigma }_{ij}{n}_j=C_{ijkl}{e}_{kl}{n}_j$ , (22)

where $C$ is the 4^th order elasticity tensor, $n$ is the unit outward normal vector and $e$ is the strain tensor which is defined as

$e_{kl}=\left(u_{k,l}+u_{l,k}\right)/2$ . (23)

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

$u_i^h\left(x\right)={\displaystyle {\sum }_{I=1}^{n_2}{\varphi }^I\left(x\right){\widehat{u}}_i^I}+{\displaystyle {\sum }_{J=1}^{n_d}{\displaystyle {\sum }_{K=1}^{n_1}{\varphi }^K\left(s\right){}_{}{}^J\widehat{q}{}_i^K{}_{}{}^J\psi \left(r\right)}}$ , (24)

where ${\widehat{u}}_i^I$ and ${}_{}{}^J\widehat{q}{}_i^K$ are the unknown fictitious nodal values for displacement fields and discontinuities amplitude parameters. By differentiating this equation one can obtain

$u_{i,j}^h\left(x\right)={\displaystyle {\sum }_{I=1}^{n_2}{\varphi }_{,j}^I\left(x\right){\widehat{u}}_i^I}+{\displaystyle {\sum }_{J=1}^{n_d}{\displaystyle {\sum }_{K=1}^{n_1}{\varphi }_{,j}^K\left(s\right){}_{}{}^J\psi \left(r\right){}_{}{}^J\widehat{q}{}_i^K}}+{\displaystyle {\sum }_{J=1}^{n_d}{\displaystyle {\sum }_{K=1}^{n_1}{\varphi }^K\left(s\right){}_{}{}^J\psi {}_{,j}\left(r\right){}_{}{}^J\widehat{q}{}_i^K}}$ (25)

where ${\varphi }_{,j}^I$ and ${\varphi }_{,j}^K$ are the derivatives of 2D and 1D MLS approximation functions, respectively, and ${}_{}{}^J\psi {}_{,j}$ are derivative of the enrichment approximation functions. Using the chain rule of differentiation

$${\varphi }_{,j}^K\left(s\right)=\frac{\partial {\varphi }^K\left(s\right)}{\partial x_j}=\frac{\partial {\varphi }^K\left(s\right)}{\partial s}\frac{\partial s}{\partial x_j}$$ , (26)

${}_{}{}^J\psi {}_{,j}\left(r\right)=\frac{\partial {}_{}{}^J\psi {\left(r\right)}^{}}{\partial x_j}=\frac{\partial {}_{}{}^J\psi {\left(r\right)}^{}}{\partial r}\frac{\partial r}{\partial x_j}$ . (27)

By substituting Equations (22) - (27) into Equation (5) , finally one obtains the linear algebraic discretized equations for a local finite volume as follows

$$K\widehat{U}=F$$ , (28)

where $\widehat{U}$ is the vector of unknowns, $K$ and $F$ are the equivalent local stiffness matrix and load vector, respectively, and can be written as follows

$$K=-{\displaystyle {\int }_{{\Gamma }_{s0}}^{}\tilde{n}DBd\Gamma }-{\displaystyle {\int }_{{\Gamma }_{su}}^{}\tilde{n}DBd\Gamma }-\alpha {\displaystyle {\int }_{{\Gamma }_{su}}^{}Nd\Gamma }$$ , (29)

$F={\displaystyle {\int }_{{\Gamma }_{st}}^{}\overline{t}d\Gamma }+{\displaystyle {\int }_{{\Omega }_s}^{}bd\Omega }-\alpha {\displaystyle {\int }_{{\Gamma }_{su}}^{}\overline{u}d\Gamma }$ , (30)

where $\overline{t}$ , $b$ and $\overline{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 $\tilde{n}$ contains the entries of the unit outward normal vector and defined as

$\tilde{n}=\left[\begin{array}{ccc}{n}_1& 0& n_2\\ 0& n_2& n_1\end{array}\right]$ . (31)

Also, the element contribution to $N$ and strain–displacement matrix $B$ are as follows

$N=\left[\begin{array}{cc}{N}_{std}^I& {}_{}{}^{J}N{}_{enr}^K\end{array}\right]$ , (32)

$B=\left[\begin{array}{cc}{B}_{std}^I& {}_{}{}^{J}B{}_{enr}^K\end{array}\right]$ , (33)

with their entries being

$\begin{array}{cc}{N}_{std}^I=\left[\begin{array}{cc}{\varphi }^I\left(x\right)& 0\\ 0& {\varphi }^I\left(x\right)\end{array}\right],& {}_{}{}^{J}N{}_{enr}^K=\left[\begin{array}{cc}{\varphi }^K\left(s\right){}_{}{}^J\psi \left(r\right)& 0\\ 0& {\varphi }^K\left(s\right){}_{}{}^J\psi \left(r\right)\end{array}\right]\end{array}$ , (34)

$$B_{std}^I=\left[\begin{array}{cc}{\varphi }_{\mathrm{,1}}^I\left(x\right)& 0\\ 0& {\varphi }_{\mathrm{,2}}^I\left(x\right)\\ {\varphi }_{\mathrm{,2}}^I\left(x\right)& {\varphi }_{\mathrm{,1}}^I\left(x\right)\end{array}\right],{}_{}{}^{J}B{}_{enr}^K=\left[\begin{array}{cc}{\varphi }_{\mathrm{,1}}^K\left(s\right){}_{}{}^J\psi \left(r\right)+{\varphi }^K\left(s\right){}_{}{}^J\psi {}_{\mathrm{,1}}\left(r\right)& 0\\ 0& {\varphi }_{\mathrm{,2}}^K\left(s\right){}_{}{}^J\psi \left(r\right)+{\varphi }^K\left(s\right){}_{}{}^J\psi {}_{\mathrm{,2}}\left(r\right)\\ {\varphi }_{\mathrm{,2}}^K\left(s\right){}_{}{}^J\psi \left(r\right)+{\varphi }^K\left(s\right){}_{}{}^J\psi {}_{\mathrm{,2}}\left(r\right)& {\varphi }_{\mathrm{,1}}^K\left(s\right){}_{}{}^J\psi \left(r\right)+{\varphi }^K\left(s\right){}_{}{}^J\psi {}_{\mathrm{,1}}\left(r\right)\end{array}\right]$$ _. (35)

Accordingly, the vector of unknowns, $\widehat{U}$ , is expressed as

$\widehat{U}=\left[\begin{array}{cc}\begin{array}{cc}{\widehat{u}}_1^I& {\widehat{u}}_2^I\end{array}& \begin{array}{cc}{}_{}{}^J\widehat{q}{}_1^K& {}_{}{}^J\widehat{q}{}_2^K\end{array}\end{array}\right]$ . (36)

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.

3.3 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 $n_1$ enrichment nodes and $n_2$ 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\left(n_1+n_2\right)$ entries. On the other hand, Equation (28) provides two discrete equations for each individual field node or corresponding finite volume, whereas there are no equations corresponding to enrichment nodes, because they have no displacement degrees of freedom. In order to satisfy the momentum balance law over the entire domain, Equation (28) should be applied over all finite volumes. The resulted equations are assembled to form the global system equations. By doing so, we have only $2\left(n_2\right)$ discrete equations. This is due to the stack form assembly procedure in finite volume methods.

In order to have the required $2\left(n_1\right)$ discrete equations, we need $n_1$ 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 $n_1$ 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.

4 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) .

4.1 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. $\Omega ={\Omega }_1\cup {\Omega }_2$ . Young's moduli in ${\Omega }_1=\left(\mathrm{0,}\Gamma \right)$ and ${\Omega }_2=\left(\Gamma ,L\right)$ are denoted by $E_1$ and $E_2$ , respectively. There is a material discontinuity at $x=\Gamma$ . The bar is fixed at one end and the other end is subjected to a prescribed displacement $\overline{u}$ . The problem is depicted in Figure 4 .

Figure 4
1D bi-material bar. Problem description. The bar is composed of two different materials, i.e. $\Omega ={\Omega }_1\cup {\Omega }_2$ , with a material discontinuity at $x=\Gamma$ . The bar is fixed at one end and the other end is subjected to a prescribed displacement $\overline{u}$ .

In the computations, let L = 10, A = 1, $\Gamma$ = 5, $E_1$ = 1, $E_2$ = 0.5 and $\overline{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.

Figure 5
1D bi-material bar. Domain discretization and corresponding finite volumes. Blue and red dots are denoted the enrichment node and field nodes, respectively. 1D local background finite volumes are illustrated by red lines.

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 )

$u\left(x\right)=\frac{1}{5\left(E_1+E_2\right)}\{\begin{array}{ll}{E}_{2}x\hfill & x\le 5\hfill \\ E_1\left(x-5\right)+5E_2\hfill & x\ge 5\hfill \end{array}$ , (37)

$\sigma \left(x\right)=\frac{E_1{E}_2}{5\left(E_1+E_2\right)}forallx$ . (38)

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 $\sqrt{{\left({\sigma }^{numerical}-{\sigma }^{exact}\right)}^2}/{\sigma }^{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%.

Figure 6
1D bi-material bar. Displacement variation along the bar.

Figure 7
1D bi-material bar. Stress variation along the bar.

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.

Figure 8
1D bi-material bar. Convergence results. (a) Meshless FVM, (b) Proposed formulation.

In the calculation of error, let ${()}^{num}$ and ${()}^{exact}$ represent approximated and exact analytical solutions, respectively. Then, the relative error norms in displacement and stress can be defined as follows

$e_d=\sqrt{\frac{{\displaystyle {\int }_{\Omega }^{}{\left(u^{num}-u^{exact}\right)}^{2}d\Omega }}{{\displaystyle {\int }_{\Omega }^{}{\left(u^{exact}\right)}^{2}d\Omega }}}$ , (39)

$e_{\sigma }=\sqrt{\frac{{\displaystyle {\int }_{\Omega }^{}{\left({\sigma }^{num}-{\sigma }^{exact}\right)}^{2}d\Omega }}{{\displaystyle {\int }_{\Omega }^{}{\left({\sigma }^{exact}\right)}^{2}d\Omega }}}$ , (40)

where $u$ and $\sigma$ are the displacement and stress at each node, respectively, and $\Omega$ 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.

Table 1
1D bi-material bar. The relative error norms.

Error norms	Proposed formulation	Meshless FVM	EFG method
$e_d$	0.37%	0.38%	0.05%
$e_{\sigma }$	3.12%	6.96%	1.03%

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.

4.2 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 $r_1$ is considered, i.e. $\Omega ={\Omega }_1\cup {\Omega }_2$ . Young's modulus and Poisson's ratio of the inclusion ${\Omega }_1$ are denoted by $E_1$ and ${\nu }_1$ and those of the matrix ${\Omega }_2$ are denoted by $E_2$ and ${\nu }_2$ , respectively. There is a material discontinuity across the interface ${\Gamma }_1$ ( $r=r_1$ ) where tractions and displacements are continuous across it. The problem is depicted in Figure 9 .

Figure 9
2D infinite plate with an inclusion. Problem description. The plate is composed of a square domain of length L with a centered circular inclusion of radius $r_1$ , i.e. $\Omega ={\Omega }_1\cup {\Omega }_2$ , with a material discontinuity across the interface ${\text{Γ}}_1$ . The model is subjected to a prescribed displacement field on the external boundary ${\Gamma }_2$ .

In the polar coordinate system, the essential boundary conditions on the external boundary

${\Gamma }_2$ are given by

$u_r=r,u_{\theta }=0$ . (41)

In the computations, L = 2, $r_1$ = 0.4, $E_1$ = 2, ${\nu }_1$ = 0.2, $E_2$ = 10, ${\nu }_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.

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.

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.

In the polar coordinate system, the exact analytical solution for displacements is as follows ( Sukumar et al. 2001 )

$u_r=\{\begin{array}{ll}\left[\left(1-\frac{L^2}{r_1^2}\right)\beta +\frac{L^2}{r_1^2}\right]r\hfill & 0\le r\le r_1\hfill \\ \left(r-\frac{L^2}{r}\right)\beta +\frac{L^2}{r}\hfill & r_1\le r\le L\hfill \end{array}$ , (42)

$u_{\theta }=0$ , (43)

and the exact analytical stress fields are given by the following equations where the shear stress field is zero ( Sukumar et al. 2001 ).

${\sigma }_{rr}=\{\begin{array}{ll}\left[\left(1-\frac{L^2}{r_1^2}\right)\beta +\frac{L^2}{r_1^2}\right]\left(2\mu +2\lambda \right)\hfill & 0\le r\le r_1\hfill \\ \left[\left(1+\frac{L^2}{r^2}\right)\beta -\frac{L^2}{r^2}\right]\left(2\mu +2\lambda \beta \right)\hfill & r_1\le r\le L\hfill \end{array}$ , (44)

${\sigma }_{\theta \theta }=\{\begin{array}{ll}\left[\left(1-\frac{L^2}{r_1^2}\right)\beta +\frac{L^2}{r_1^2}\right]\left(2\mu +2\lambda \right)\hfill & 0\le r\le r_1\hfill \\ \left[\left(1-\frac{L^2}{r^2}\right)\beta +\frac{L^2}{r^2}\right]\left(2\mu +2\lambda \beta \right)\hfill & r_1\le r\le L\hfill \end{array}$ , (45)

where

$\lambda$ and $\mu$ are appropriate Lame constants at each material and $\beta$ is defined as

$\beta =\frac{\left({\lambda }_1+{\mu }_1+{\mu }_2\right)L^2}{\left({\lambda }_2+{\mu }_2\right)r_1^2+\left({\lambda }_1+{\mu }_1\right)\left(L^2-r_1^2\right)+{\mu }_2{L}^2}$ . (46)

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 $u_x$ and the normal stress ${\sigma }_{xx}$ both along the line $y=x$ , respectively. A magnified view around the discontinuities is presented in Figure 12 for better observation.

Figure 11
2D infinite plate with an inclusion. Variation of the $u_x$ along the line $y=x$ .

Figure 12
2D infinite plate with an inclusion. Variation of the ${\sigma }_{xx}$ along the line $y=x$ . A magnified around the discontinuities view.

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 ${\sigma }_{xx}$ is constant within the inclusion. The L2 relative error norm of the predicted ${\sigma }_{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

$e_d=\sqrt{\frac{{\displaystyle {\int }_{\Omega }^{}\left(u^{num}-u^{exact}\right).\left(u^{num}-u^{exact}\right)d\Omega }}{{\displaystyle {\int }_{\Omega }^{}{u}^{exact}.u^{exact}d\Omega }}}$ , (47)

$$e_{\sigma }=\sqrt{\frac{{\displaystyle {\int }_{\Omega }^{}\left({\sigma }^{num}-{\sigma }^{exact}\right).\left({\sigma }^{num}-{\sigma }^{exact}\right)d\Omega }}{{\displaystyle {\int }_{\Omega }^{}{\sigma }^{exact}.{\sigma }^{exact}d\Omega }}}$$ , (48)

where $u$ and $\sigma$ are the vectors of displacement and stress at each point, respectively, and $\Omega$ 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.

Table 2
2D infinite plate with an inclusion. The relative error norms.

Error norms	Proposed formulation	Meshless FVM	EFG method
$e_d$	0.48%	2.56%	1.50%
$e_{\sigma }$	1.96%	20.20%	8.05%

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 .

Figure 13
2D infinite plate with an inclusion. Convergence results.

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).

4.3 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 $E_1$ = 0.02 while $E_2$ 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 $u_x$ and the normal stress ${\sigma }_{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
2D infinite plate with an extremely soft inclusion. Variation of the $u_x$ along the line $y=x$ .

Figure 15
2D infinite plate with an extremely soft inclusion. Variation of the ${\sigma }_{xx}$ along the line $y=x$ . A magnified around the discontinuities view.

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 ${\sigma }_{xx}$ is constant within the inclusion. The L2 relative error norm of the predicted ${\sigma }_{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 .

Table 3
2D infinite plate with an extremely soft inclusion. The relative error norms.

Error norms	Proposed formulation	Meshless FVM	EFG method
$e_d$	0.97%	3.77%	3.23%
$e_{\sigma }$	3.40%	31.09%	7.00%

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.

Figure 16
2D infinite plate with an extremely soft inclusion. Convergence results.

5 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.

Appendix A Transformation of field variables

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 .

Figure A.1
Cartesian and polar coordinate systems.

Displacement transformation

$$u_x=\cos \theta u_r-\sin \theta u_{\theta }$$ (A.1)

$u_y=\sin \theta u_r+\cos \theta u_{\theta }$ (A.2)

Stress transformation

$${\sigma }_x=\frac{1}{2}\left[\left(1+\cos 2\theta \right){\sigma }_r+\left(1-\cos 2\theta \right){\sigma }_{\theta }-\left(2\sin 2\theta \right){\sigma }_{r\theta }\right]$$ (A.3)

${\sigma }_y=\frac{1}{2}\left[\left(1-\cos 2\theta \right){\sigma }_r+\left(1+\cos 2\theta \right){\sigma }_{\theta }+\left(2\sin 2\theta \right){\sigma }_{r\theta }\right]$ (A.4)

${\sigma }_{xy}=\frac{1}{2}\left[\left(\sin 2\theta \right){\sigma }_r-\left(\sin 2\theta \right){\sigma }_{\theta }+\left(2\cos 2\theta \right){\sigma }_{r\theta }\right]$ (A.5)

• http://dx.doi.org/10.1590/1679-78254121

• Erratum

  In the article An enriched meshless finite volume method for the modeling of material discontinuity problems in 2D elasticity, doi number: http://dx.doi.org/10.1590/1679-78254121, published in journal Latin American Journal of Solids and Structures, 15(2):e13, page 1:
  Where it reads:
  Abdullah Davoudi-Kia^a
  N. Fallah^a,b,*
  ^a Department of Civil Engineering, University of Guilan, P.O. Box 3756, Rasht, Iran.
  ^b Faculty of Civil Engineering, Babol Noshirvani University of Technology, Babol, Iran.
  [e-mail]
  *[e-mail]
  It should read:
  Abdullah Davoudi-Kia^a
  N. Fallah^a,b,*
  ^aDepartment of Civil Engineering, University of Guilan, P.O. Box 3756, Rasht, Iran. E-mail: [e-mail]
  ^bFaculty of Civil Engineering, Babol Noshirvani University of Technology, Babol, Iran. E-mail: [e-mail]
  *Corresponding author

References

• An, X., Ma, G., Cai, Y., Zhu, H. (2011). A new way to treat material discontinuities in the numerical manifold method. Computer Methods in Applied Mechanics and Engineering, 200(47): 3296-3308.
• Atluri, S., Han, Z., Rajendran, A. (2004). A new implementation of the meshless finite volume method, through the MLPG “mixed” approach. CMES: Computer Modeling in Engineering & Sciences, 6(6): 491-513.
• Atluri, S.N., Shen, S. (2002). The Meshless Local Petrov-Galerkin (MLPG) Method: A Simple\ & Less-costly Alternative to the Finite Element and Boundary Element Methods. CMES: Computer Modeling in Engineering & Sciences, 3(1): 11-52.
• Atluri, S.N., Zhu, T. (1998). A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics. Computational mechanics, 22(2): 117-127.
• Batra, R., Porfiri, M., Spinello, D. (2004). Treatment of material discontinuity in two meshless local Petrov-Galerkin (MLPG) formulations of axisymmetric transient heat conduction. International Journal for Numerical Methods in Engineering, 61(14): 2461-2479.
• Belytschko, T., Black, T. (1999). Elastic crack growth in finite elements with minimal remeshing. International journal for numerical methods in engineering, 45(5): 601-620.
• Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., Krysl, P. (1996). Meshless methods: an overview and recent developments. Computer methods in applied mechanics and engineering, 139(1): 3-47.
• Belytschko, T., Lu, Y.Y., Gu, L. (1994). Element-Free Galerkin Methods. International journal for numerical methods in engineering, 37(2): 229-256.
• Cardiff, P., Karač, A., Ivanković, A. (2014). A large strain finite volume method for orthotropic bodies with general material orientations. Computer Methods in Applied Mechanics and Engineering, 268 318-335.
• Cavalcante, M.A., Marques, S.P., Pindera, M.-J. (2007). Parametric formulation of the finite-volume theory for functionally graded materials—part I: analysis. Journal of Applied Mechanics, 74(5): 935-945.
• Cavalcante, M.A., Pindera, M.-J. (2016). Generalized FVDAM theory for elastic–plastic periodic materials. International Journal of Plasticity, 77 90-117.
• Cordes, L., Moran, B. (1996). Treatment of material discontinuity in the element-free Galerkin method. Computer Methods in Applied Mechanics and Engineering, 139(1): 75-89.
• de Borst, R. (2008). Challenges in computational materials science: multiple scales, multi-physics and evolving discontinuities. Computational Materials Science, 43(1): 1-15.
• Ebrahimnejad, M., Fallah, N., Khoei, A. (2014). New approximation functions in the meshless finite volume method for 2D elasticity problems. Engineering Analysis with Boundary Elements, 46 10-22.
• Ebrahimnejad, M., Fallah, N., Khoei, A. (2015). Adaptive refinement in the meshless finite volume method for elasticity problems. Computers & Mathematics with Applications, 69(12): 1420-1443.
• Ebrahimnejad, M., Fallah, N., Khoei, A.R. (2017). Three types of meshless finite volume method for the analysis of two-dimensional elasticity problems. Computational and Applied Mathematics, 36(2): 971-990.
• Escarpini Filho, R.d.S., Marques, S.P.C. (2016). A model for homogenization of linear viscoelastic periodic composite materials with imperfect interface. Latin American Journal of Solids and Structures, 13(14): 2706-2735.
• Fallah, N. (2004). A cell vertex and cell centred finite volume method for plate bending analysis. Computer methods in applied Mechanics and Engineering, 193(33): 3457-3470.
• Fries, T.-P., Belytschko, T. (2010). The extended/generalized finite element method: an overview of the method and its applications. International Journal for Numerical Methods in Engineering, 84(3): 253-304.
• Gdoutos, E.E. (2005). Fracture mechanics: an introduction. Springer Science & Business Media, Netherlands.
• Han, Z., Rajendran, A., Atluri, S. (2005). Meshless Local Petrov-Galerkin (MLPG) Approaches for Solving Nonlinear Problems with Large Deformations and Rotations. CMES: Computer Modeling in Engineering & Sciences, 10(1): 1-12.
• Hosseini, S.M., Sladek, J., Sladek, V. (2011). Meshless local Petrov–Galerkin method for coupled thermoelasticity analysis of a functionally graded thick hollow cylinder. Engineering Analysis with Boundary Elements, 35(6): 827-835.
• Hu, M., Wang, Y., Rutqvist, J. (2015). An effective approach for modeling fluid flow in heterogeneous media using numerical manifold method. International Journal for Numerical Methods in Fluids, 77(8): 459-476.
• Krongauz, Y., Belytschko, T. (1998). EFG approximation with discontinuous derivatives. International Journal for Numerical Methods in Engineering, 41(7): 1215-1233.
• Li, Q., Shen, S., Han, Z., Atluri, S. (2003). Application of meshless local Petrov-Galerkin (MLPG) to problems with singularities, and material discontinuities, in 3-D elasticity. Computer Modeling in Engineering and Sciences, 4(5): 571-586.
• Li, S., Hao, W., Liu, W.K. (2000). Mesh-free simulations of shear banding in large deformation. International Journal of solids and structures, 37(48): 7185-7206.
• Liu, G.-R. (2009). Meshfree methods: moving beyond the finite element method. CRC Press, Taylor & Francis, Boca Raton.
• Liu, G.-R., Gu, Y.-T. (2005). An introduction to meshfree methods and their programming. Springer Science & Business Media, Netherlands.
• Melenk, J.M., Babuška, I. (1996). The partition of unity finite element method: basic theory and applications. Computer methods in applied mechanics and engineering, 139(1): 289-314.
• Moës, N., Cloirec, M., Cartraud, P., Remacle, J.-F. (2003). A computational approach to handle complex microstructure geometries. Computer methods in applied mechanics and engineering, 192(28): 3163-3177.
• Nguyen, V.P., Rabczuk, T., Bordas, S., Duflot, M. (2008). Meshless methods: a review and computer implementation aspects. Mathematics and computers in simulation, 79(3): 763-813.
• Okabe, A., Boots, B., Sugihara, K., Chiu, S.N. (2009). Spatial tessellations: concepts and applications of Voronoi diagrams. John Wiley & Sons, England.
• Sladek, J., Sladek, V., Solek, P., Saez, A. (2008). Dynamic 3D axisymmetric problems in continuously non-homogeneous piezoelectric solids. International Journal of Solids and Structures, 45(16): 4523-4542.
• Soares, D., Sladek, V., Sladek, J. (2012). Modified meshless local Petrov–Galerkin formulations for elastodynamics. International Journal for Numerical Methods in Engineering, 90(12): 1508-1528.
• Sukumar, N., Chopp, D.L., Moës, N., Belytschko, T. (2001). Modeling holes and inclusions by level sets in the extended finite-element method. Computer methods in applied mechanics and engineering, 190(46): 6183-6200.
• Taylor, G., Bailey, C., Cross, M. (2003). A vertex‐based finite volume method applied to non‐linear material problems in computational solid mechanics. International journal for numerical methods in engineering, 56(4): 507-529.
• Ventura, G., Gracie, R., Belytschko, T. (2009). Fast integration and weight function blending in the extended finite element method. International Journal for Numerical Methods in Engineering, 77(1): 1-29.
• Yoon, Y.-C., Song, J.-H. (2014). Extended particle difference method for weak and strong discontinuity problems: part I. Derivation of the extended particle derivative approximation for the representation of weak and strong discontinuities. Computational Mechanics, 53(6): 1087-1103.

Publication Dates

• Publication in this collection
  2018

History

• Received
  14 June 2017
• Reviewed
  01 Dec 2017
• Accepted
  01 Dec 2017