Numerical instability in the element-free Galerkin method for bi-material modeling

This paper presents numerical results from modeling of bi-material problems using the element-free Galerkin method. The conventional EFG shape functions do not pass through the nodal data, imposition of nodal constraints brings difficulties, thus some special techniques must be employed. In bi-material problems, connecting the two subdomains of different properties may still be possible via Lagrange multipliers, penalty function, Nitsche’s formulation,or direct imposition when using a regularized weighting function to transform the EFG shape functions. All these mentioned techniques are explored and compared in this paper. Some remarks about numerical instabilities are reported. c © 2018 University of West Bohemia. All rights reserved.


Introduction
The meshless methods have been a research focus for about two decades.Among all the meshless methods, the element-free Galerkin (EFG) method has been the most popular and implemented for use in some commercial codes.Despite many advantages, one of the main shortcomings of the EFG is the lack of Kronecker delta property, leading to difficulties when the constraints are to be imposed, either in term of external constraints (boundary conditions) or the internal constraints (tyings).
In the early versions of EFG, special techniques were needed in order to enforce constraint equations.The original EFG version [3] suggested the use of Lagrange multipliers because of its simplicity.Satisfactory results were obtained, but at the same time, it introduced an extra set of unknowns leading to a larger size and bandwidth of the global system of equations.As alternatives to the use of Lagrange multipliers, a modified variational principle [17] was proposed.Banded and positive-definite discrete equations were created, while the solution accuracy was compromised.For two dimensional problems, a modified definition of the discrete L 2 norm proposed in [20] was proven another good choice.Instead of using nodal displacements in the formulation of the discrete L 2 norm, the approximate ones were employed.A modified collocation method [29] was proven to be consistent with the variational basis of the EFG method, provided results with better accuracy than those from [3].In the same work [29], a penalty formulation gave considerably accurate solutions with simpler implementation as compared to the Lagrange multipliers method.The formulation relied on the use of penalty parameter that must be large enough, and the convergence of the solutions could be reached either by the positive or negative parameters [1].Similar to both the Lagrange multipliers and the penalty formulation, the Nitsche's formulation [22] was revisited in [7,15].The method was found to resemble some characteristics of the Lagrange multipliers formulations with better convergence property, and not to suffer from the ill-conditioning problem often found when employing the penalty formulation.
There were also some attempts to directly impose the boundary conditions.In [13], a layer of finite elements was placed along the displacement boundaries, allowing direct imposition of the constraints.Similarly, the use of singular weight functions [10] allowed the approximants to pass through the corresponding nodes, facilitating the direct imposition of essential boundary conditions.In [8], a nodal integration approach was proposed for elliptic equations by discretizing the boundary using a covering point cloud.Another direct imposition approach was proposed in [28] by using the local Kronecker delta property in the EFG shape function.In [9], enforcement of the essential boundary conditions was possible by distributing forces at the essential boundary nodes, which was suitable for explicit time integration problems.The direct imposition was also possible with a transformation matrix in [4].In [21], a regularized singular weight function was used.The enforcement of the boundary conditions was based on the interpolating moving least squares, also satisfying the Kronecker delta property.In [18,19], special weight functions were employed in formulating the EFG shape functions.Interestingly, the non-singular weight functions proposed in [19] leaded to the EFG shape functions that contained the Kronecker delta property.The corresponding results were also much more independent of the nodal configuration and the size of the domain of influence.The technique was later investigated in [27] as the interpolating moving least squares method (IIMLS), reporting higher computational precision than the original EFG method with ease of constraint enforcement.
In modeling of bi-material problems using the EFG method, the Lagrange multipliers were applied to enforce a set of tied constraints along the interface of bi-material problems in [5].Nodes along the interface were doubled and superimposed, belonging solely to each material.The inhomogeneous body was separated into two homogeneous parts.In [12], jump shape functions were added in the formulation.A trial function with a pre-imposed discontinuity in the gradient of the function was added at the location of the material discontinuity for 2D elastic problems.By this technique, solution oscillations were avoided at the surface of discontinuity.However, this method required interpolation in the curve linear coordinates, which became very hard in three dimensional applications.In [4], a method for direct imposition of essential boundary conditions and treatment of material discontinuity in EFG method was presented.By using the actual displacements at the nodes on the displacement boundaries and the material interface in each material domain, the stiffness matrix and load vector at an integral point were rewritten and transformed.Similar errors found when using the Penalty method and Lagrange multiplier method were still observed.In [6], a penalty method was used in modeling composite problems.The continuity of tractions at the interface was weakly satisfied at the variational level.The actual displacements at the interface were determined, and after that, conditions of continuity of displacements at the interface were also directly enforced as in the FE method.The method showed good agreement with the FEM solution, but with a smaller number of degrees of freedom compared to the FEM.
Although many techniques were proposed for the EFG method to deal with constraint enforcement difficulties, there were still some problems to be aware of.It was reported in [23] that some unstable solutions were observed in the displacement field of the Prandtl's punch test, when the constraints were imposed by the use of Lagrange multipliers and penalty functions.The displacement oscillations appeared where there was an abrupt change in boundary conditions, and they became more obvious when the internal (tying) constraint equations were imposed.It could then be a problem in the bi-material modeling where there existed tying constraints.In the finite element context, some alternative variational methods for enforcing a tying constraint on an enriched interface in the context of two-dimensional elasticity was examined in [26].Lagrange multiplier method, penalty method, Nitsche's method and the bubble-stabilization method were used in imposing tying constraints via enriched finite element methods (X-FEM).It was shown in this work some disadvantages of the penalty method including the tendency toward ill-conditioning and the sensitivity of the method to the penalty parameter.It was mentioned in the work that Nitsche's method would be a stable and convergent method to impose tying constraints.All of their numerical examples showed optimal rates of convergence in the interfacial traction by using Nitsche's method.Another approach suggested the material discontinuity at the interface to be modeled using a jump function [25], created by enriching the approximation with a special shape function containing discontinuities.The mixed mode stress intensity factors for bi-material interface cracks were numerically evaluated using the modified domain form of interaction integral approach.In dealing with discontinuous interface, three different approaches namely the domain partitioning, the Lagrange multiplier and the jump function were compared [24].No solution instability was reported.
In this study, we investigate techniques used for modeling bi-material composite problems.To get an accurate solution across interface of the different materials, the constraint enforcement plays an important role.Following [26], we investigate the use of Lagrange multipliers, the penalty formulation, and the Nitsche's formulation in the constraint enforcement.Additionally, it is interesting to add a method of direct constraint imposition via the modified shape functions, based on the regularized weighting functions [18], as a comparison.The last method is different than the other three that the constraints were to be imposed directly in the system of equations, not pre-included in the weak formulation.Our study explores the enforcement techniques for two different kinds of constraints, i.e. the external and the internal constraints along the shared interface of the bi-material problems, using the EFG method.Comparisons in term of numerical instabilities are made and discussed.

The conventional EFG shape functions
The basic concept of the EFG method requires use of the moving least squares (MLS) approximation [14] to form the shape functions.In a two-dimensional problem, an arbitrary function u is approximated at point x = (x, y) by where p(x) contains polynomial basis of order m and thus n = 0.5(m + 1)(m + 2).The vector of coefficients a(x), containing coefficients of the polynomial, can be determined by minimizing a weighted least squares sum J MLS given by The weighting function w i (x) of node i covers the domain of influence centering at the location of node i.In the classical element-free Galerkin method, one of the most common weighting function is the Gaussian (exponential) weighting function [3] defined as where the parameter z is taken as 0.4 as in [11].The variable s is the normalized distance between the interpolation point and the considered supporting point, i.e.
Here, d is the radius of the circular-shape domain of influence, or as known as the influence radius.
Minimizing J MLS with respect to a(x), the approximate function u h (x) can be written in the form of the shape function N(x) as where û contains the nodal parameters u i , and A(x) and B(x) can be determined by where, for m supporting points and n polynomial coefficients, P and W(x) are defined as and

Imposition of constraints
A solid body is defined by the domain Ω, with the outer boundary and the inner boundary that divides the domain into two parts.The outer boundary consists of the displacement boundary Γ d and the stress boundary Γ n , and the inner boundary Γ c is where the degrees of freedom are tied.The governing equilibrium equation is where σ is the stress tensor, b represents the body force vector and ∇ is the gradient operator.
The boundary conditions are given as where t and u are the prescribed tractions on Γ n and the prescribed displacements on Γ d , respectively.The external constraint equations, for enforcement of the prescribed degrees of freedom u, can be written in the form Equation ( 13) represents a set of external constraints.On the other hand, we can also enforce a set of internal constraints, which implies that two sets of degree of freedom displace at the same amount, in a similar format.The internal constraint equations show ties between some slave degrees of freedom u s to some master degrees of freedom u m on a boundary Γ c so that they displace in the same amount.The corresponding internal constraint equations, for enforcement of tied degrees of freedom, can be written as In this study, we consider four techniques to impose constraint equations, including method of Lagrange multipliers, penalty formulation, Nitsche's method, and direct imposition.Formulations and implementation details of the four methods are given as follows.

Method of Lagrange multipliers
The Galerkin weak form can be written in the form The displacements and the Lagrange multipliers can be discretized as u = Nû and λ = H λ, where N and H are the shape functions, and û and λ are the nodal degrees of freedom and Lagrange multipliers.Strain tensor is the symmetric tensor of the displacement gradient and can be expressed in the discrete form as ε = ∇u = ∇Nû = Bû.The stress tensor is given by the Hooke's law, i.e. σ = Dε, where D is the tensor of elastic constants.The resulting set of equations becomes where

Penalty formulation
The Galerkin weak form can be written in the form where α is the penalty parameter, generally taken in the range of 10 3 to 10 7 [29].The discretization results in where K and f are previously defined as in ( 17) and (20), respectively.In addition,

Nitsche's method
The weak form for elasticity problem can be written as leading to the discrete form The additional matrices in (28) are defined as Note that the parameter β in ( 28) is a stabilization parameter for Nitsche's method.The number should be sufficiently large to guarantee positive definiteness of the system equation.
The formulations given in ( 16)-(32) are applicable for the bi-material problem, where some degrees of freedom on Material 1 are tied to some degrees of freedom on Material 2. The subscripts s and m refer to the slave and the master degrees of freedom, respectively.

Direct imposition
The constraint equations in the EFG method can be imposed directly via use of the regularized weighting function [18], instead of the conventional Gaussian weighting function, in formulation of the EFG shape functions (cf.Fig. 1).The resulting shape functions satisfy Kronecker delta property, therefore they facilitate imposition of constraints at nodes.The weighting function can be computed from where for compact support.In [18], it was recommended to use the parameters γ = 2 and μ = 10 −5 .

Numerical examples
In this paper, three numerical problems are chosen to study performance of different constraint enforcement techniques.Details of each problem and its results are given in the following subsections.The analytical displacements and strains for this problem [13,16] are

One-dimensional bar
where For this problems, the four constraint enforcement techniques (as mentioned in Section 3) are used.The obtained results are comparatively shown in Fig. 3, also with various sizes of The obviously unstable solutions occur when the domain of influence reaches the size of 3.0 times the nodal distance, when imposing constraints with method of Lagrange multipliers, penalty formulation and Nitsche's method.However, when using direct imposition approach via the regularized weighting function, the obtained solutions (cf.Fig. 3) show less errors, are not sensitive to the size of the domain of influence, and do not exhibit any oscillation in the solution field.

Beam problem
Following [26], a beam of 4-unit depth and 16-unit length is analyzed under plane stress condition.At the left end of the beam, the center node is restrained both in x and in y directions, whereas the other nodes are restrained only in the x-direction (cf.Fig. 4).The normal and shear tractions are prescribed as where t x is considered as the bending stress, t y is the shear stress, and c is the distance from the neutral axis to the top edge of the beam.Here, our input variables include the shear force P = −1 000, the distance c = 2, and the moment of inertia I = (2c) 3 /12.The Young's modulus E = 10 5 is used for both Material 1 (the left half part) and Material 2 (the right half part), and the Poisson's ratio is 0.25 for both subdomains.It should be noted that this problem is not a truly bi-material problem, since the material properties of the two subdomains are the same.For this example, our purpose is to examine effect of constraint enforcement techniques, regardless of material differences in the two domains, on the EFG solution.The problem domain is divided into the left half part and the right half part, each of which contains a set of uniformly distributed 9 × 9 nodes.Totally 18 interfacial nodes at x = 8, 9 on the left and 9 on the right, are tied.The 16 × 8 integration cells are used, with 4 × 4 Gaussian points per cell.The nodes and the integration points are indicated by circular and cross marks, respectively, in Fig. 5.The penalty and the Nitsche parameters are chosen as 10 6 and 10 4 , respectively.The analytical displacements for this problem, as mentioned in [2], are The analytical stresses are In this test, four sizes of the domain of influence are applied as 1.5, 2.5, 3.5 and 4.5 times the shortest distance between nodes.Plots of the shear stress (σ xy ) in the interfacial line (x = L/2) are shown in Fig. 6.Similar to the bar problem, some instabilities of the solutions can be noticed, especially when the domain of influence is large.Comparing the results from all the approaches, the size that gives the most accurate solution of all the computations is the size of 2.5, despite appearance of some yet small oscillations in the traction fields in most approaches.The most sensitive approach is the method of Lagrange multipliers, of which its corresponding solutions show oscillations when the size is larger than 1.5, cf.Fig. 6(a).In Fig. 6(b), the seriously unstable solutions are obtained when the size of the domain of influence is of 3.5 when using the penalty formulation to impose the internal constraints.The relatively serious phenomenon is also noticed when using Nitsche's method.However, the serious oscillation when applying Nitsche's method appears with a larger size of the domain of influence (d = 4.5), cf.Fig. 6(c).Unlike the other techniques and even the finite element solution [26], the results obtained from the EFG method with direct constraint imposition, cf.Fig. 6(d), do not show any oscillations in the traction field for all the chosen sizes of the domain of influence.

Inclusion in an infinite plate
Following the work of [5], the infinite plate problem (cf.Fig. 7) is modeled under plane stress condition.Only a quarter part is used due to symmetry of its domain.The problem domain is divided into two materials; the Young's moduli of the inner and the outer parts are E 1 = 800 and E 2 = 1 000, respectively.The Poisson's ratio is ν = 0.28, for both materials.Only one-quarter is modeled using totally 126 nodes, with 18 nodes tied on the interface at r = 2, as shown in Fig. 8 The infinite plate is modeled by applying the exact displacements on all boundaries of both materials.The penalty and the Nitsche parameters are chosen as 10 6 and 10 4 , respectively.Following [16], the analytical displacement is provided in cylindrical coordinates as The corresponding strains are and where and the Lame ´parameters of the material i are defined as The eigenstrain is assumed as a constant dilational strain ε * 1 = 0.01, following [5].Again, size of the domain of influence is varied to test instability of EFG solutions.The sizes of 2.0, 2.5, 3.0 and 3.5 times the shortest distance between nodes are used.In Fig. 9, the hoop strains on the interface between two materials are plotted.It can be seen that the best solutions from all approaches are those using size of the domain of influence d = 2.5.Using the smaller and the larger domains of influence, some instabilities in the hoop strain fields appear.Among all the EFG methods in this test, the least sensitive approach to the change of the size of the domain is obviously the EFG method based on regularized weighting function, where the constraints are strongly imposed.

Conclusions
This paper investigates techniques of constraint imposition in the EFG method for the bi-material problems.External constraints and internal constraints are enforced in the EFG formulation.Four techniques are chosen, namely use of Lagrange multipliers, penalty formulation, Nitsche's method and direct imposition via use of regularized weighting function.The first three techniques apply weak imposition of the constraint equations in the EFG formulation, whereas strong (pointwise) imposition of the constraints can be employed in the last technique as the modified shape function contains a Kronecker delta property.The results obtained from all the techniques are close but not exactly the same.For some cases, oscillations in the solution field appear, especially when larger domains of influence are used.The appearance of such numerical instability is due to the smooth character of the EFG shape function.With larger domain of influence, the EFG shape functions are smoothed out thus its nodal value goes further from one.The shape functions are, therefore, not local enough to provide an exact constraint imposition, especially when weak imposition techniques are used.This problem becomes apparent in the bi-material problems where constraints are tied at the interface.Among all the approaches, the solutions from the EFG method based on regularized weighting function appears to suffer least from such numerical instabilities; as the constraints can be directly imposed, the solutions do not exhibit large oscillations even when large domain of influence is used.The technique is also the least sensitive to change of the size of the domain of influence.

Fig. 1 .
Fig. 1.The EFG shape functions at the node x = 2 based on Gaussian weighting function (left) and based on the regularized weighting function (right)

Fig. 3 .
Fig. 3. Comparison of displacement functions in the bar problem

Fig. 5 .
Fig. 5. Nodal discretization, background integration cells and integration points used in the beam problem

Fig. 6 .
Fig. 6.Comparison of shear stress at the interface in the beam problem

Fig. 8 .
Fig. 8. Nodal discretization, background integration cells and integration points used in the plate problem

5 (Fig. 9 .
Fig. 9. Comparison of hoop strains at the interface in the plate problem