Consistent integration schemes for meshfree analysis of strain gradient elasticity

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 357 (2019) 112601 www.elsevier.com/locate/cma

Consistent integration schemes for meshfree analysis of strain gradient elasticity BingBing Wanga , Chunsheng Lub , CuiYing Fana ,∗, MingHao Zhaoa,c ,∗ a School of Mechanics & Engineering Science, Zhengzhou University, Zhengzhou 450001, China School of Civil and Mechanical Engineering, Curtin University, Western Australia 6845, Australia c Henan Key Engineering Laboratory for Anti-fatigue Manufacturing Technology and School of Mechanical Engineering, Zhengzhou University, Zhengzhou, Henan 450001, China b

Received 3 June 2019; received in revised form 17 August 2019; accepted 18 August 2019 Available online xxxx

Highlights • Arbitrary-order integration constraint conditions are derived for meshfree analysis of strain gradient elasticity. • Consistent integration schemes designed to meet constraint conditions are able to pass patch tests. • The consistent integration is superior to the standard Gaussian one based on convergence, accuracy and efficiency.

Abstract Integration schemes with nodal smoothed derivatives, which meet integration constraint conditions, are robust and efficient for use in meshfree Galerkin methods, however, most of them are focussed on the classical elasticity determined by a secondorder partial differential equation. In this paper, arbitrary-order integration constraint conditions are derived for strain gradient elasticity in a fourth-order partial differential equation. These integration constraint conditions provide the discrete forms of nodal shape functions and their first- and second-order derivatives. Furthermore, to meet the integration constraint conditions, consistent integration schemes are designed with nodal smoothed (but not standard) derivatives at evaluating points. It is shown that such nodal smoothed derivatives are able to satisfy the differentiation of approximation consistency. Finally, several case studies are given and the results demonstrate that, based on convergence, accuracy and efficiency, the numerical performance of consistent integration in meshfree analysis of strain gradient elasticity is superior to the standard Gaussian one. c 2019 Published by Elsevier B.V. ⃝ Keywords: Meshfree; Numerical integration; Integration constraint; Smoothed derivatives; Strain gradient; Consistency

1. Introduction Classical continuum theories have been widely used in structure analysis and design. Owing to the lack of an absolute size of microstructures, however, classical continuum models cannot be applied to tackle the scale effect and material instability. In contrast, the strain gradient elasticity, firstly introduced as a GradEla model by Aifantis [1], eliminates strain singularities from dislocations and crack tips, which can describe the size effect of ∗

Corresponding authors. E-mail addresses: [email protected] (C. Fan), [email protected] (M. Zhao).

https://doi.org/10.1016/j.cma.2019.112601 c 2019 Published by Elsevier B.V. 0045-7825/⃝

2

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

materials [2–5]. In the simplest form of a strain gradient theory proposed by Altan and Aifantis [2], there is only one additional material parameter in the constitutive law. As discussed in [6], this can be considered as a special case of a generalized continuum [7–9]. The strain gradient is governed by a fourth-order partial differential equation and thus, the C1 continuity of an approximation function is required in a numerical method. This makes it difficult for using a traditional finite element method, such as the complicated C1 finite or mixed elements [10–13] to meet the Ladyzhenskaya–Babuska–Brezzi stability condition [14]. To solve the problem, one strategy is to split a fourth-order partial differential equation into two sets of second-order equations [15], which is further extended to dynamic analysis [16,17]. Meanwhile, meshfree/meshless methods [18–23] developed over the last two decades have an advantage in constructing highorder and smoothed approximations. It is therefore convenient to analyze strain gradient elasticity [24,25] and gradient plasticity [26–28] by using meshfree/meshless methods. In a meshfree method, the widely used moving least-square (MLS) [18] and reproducing kernel (RK) [19] approximations are nonpolynomial rational functions with overlapping nodal support domains [29]. It is thus difficult to integrate a Galerkin weak form. To obtain a stable numerical result, a number of quadrature points are usually employed in a background integration cell but with a lower computational efficiency. The development of stable and efficient domain integration schemes for meshfree Galerkin methods is therefore an active topic [21,30,31], and various approaches such as nodal, support and stress point integration methods, have been proposed [32–37]. In addition to consistency of nodal MLS/RK shape functions, the integration constraint condition should also be followed in linear patch tests [30]. With this motivation, Chen et al. [30] introduced a strain smoothing technique, leading to the establishment of a stabilized conforming nodal integration (SCNI) method. Here it is worth noting that SCNI is able to pass linear patch tests by using minimum quadrature points (nodes). Such a technique has a higher computational efficiency than the standard Gaussian integration (GI) [38–41], but it is only adequately suitable for a linear meshfree method [42]. Duan et al. [42–44] extended SCNI to high-order meshfree methods with relevant consistent integration schemes. Alternatively, an arbitrary-order variationally consistent integration was introduced in a meshfree Galerkin method [45], which belongs to a category of Petrov–Galerkin forms. In addition, Wang and Wu [46] evaluated integration errors of a strain smoothing technique in a second-order meshfree method with two-level triangular subdomains and then, they developed an integration scheme with quadratic exactness. However, further extension to a three-dimensional formulation seems not so straightforward, and an inherently consistent reproducing kernel gradient smoothing framework has thus been proposed [47]. Although integration schemes [30,41–47] with smoothed nodal derivatives in high-order meshfree Galerkin methods are superior to standard GI, they are only suitable to classical elasticity analysis with smoothed first-order nodal derivatives. In the case of strain gradient elasticity, however, stress on a material point depends on strain and its gradient, leading to a fourth-order partial differential equation, and there are thus both first- and secondorder nodal derivatives in stiffness matrices. Therefore, these integration methods [30,39,41–47] cannot be directly applied for strain gradient elasticity. It is thus necessary to derive high-order integration constraint conditions for meshfree strain gradient analysis and to design efficient and stable integration schemes consistent with high-order approximations. The remainder of this paper is organized as follows. In Section 2, the strain gradient formulation is briefly reviewed and on the basis of MLS/RK approximations, discrete equations in a strain gradient theory are proposed in Section 3. Section 4 is then dedicated to discuss high-order integration constraint conditions for strain gradient analysis and relevant consistent integration schemes, with numerical results shown in Section 5. Finally, some concluding remarks are given in Section 6. 2. Governing equations and the Galerkin weak form In a displacement field u i , the linear strain can be expressed as ) 1( u i, j + u j,i , (1) εi j = 2 where the subscript comma means partial differentiation. The strain gradient in type II form [7] is defined as derivatives of a linear strain, or ) 1( κi jk = εi j,k = u i, jk + u j,ik . (2) 2

3

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Stress σi j and double stress τi jk , which are respectively conjugated with strain and strain gradient, are defined as ∂Ψ , ∂εi j ∂Ψ = , ∂κi jk

σi j =

(3)

τi jk

(4)

where Ψ is the strain energy density. Here, the simplest format of a linear strain gradient elasticity is ( ) 1 1 Ψ = λεii ε j j + µεi j εi j + l 2 λκiik κ j jk + µκi jk κi jk , (5) 2 2 where λ and µ are Lame parameters and l is a parameter related to the internal length scale of a material. The corresponding surface and double tractions are given by ( ) ( ) ( ) tk = n i σik − τ jik, j + n i n j τi jk D p n p − D j n i τi jk , (6) rk = n i n j τi jk ,

(7)

where n i is the ith component of a unit surface normal vector and ( ) D j (·) = δ jk − n j n k (·),k

(8)

is a surface gradient operator. The relevant surface normal gradient operator is defined as D (·) = n k (·),k .

(9)

The equilibrium equation of a body with domain Ω and smooth boundary Γ is σik,i − τ jik, ji + bk = 0,

(10)

where bk is the kth direction body force. The essential and conjugated boundary conditions are u i = u i or ti = t i , u i, j n j = wi or ri = r i .

(11) (12)

The standard Galerkin weak form of Eqs. (10)–(12) is ∫ ∫ ∫ ∫ ( ) δu i t i dΓ + δu i bi dΩ + δεi j σi j + δκi jk τi jk dΩ = Ω

Ω

Γ1

( ) δ u i, j n j r i dΓ ,

(13)

Γ2

where Γ1 and Γ2 denote the natural boundaries on which surface and double tractions are given, respectively. 3. Meshfree approximation and discretization 3.1. MLS/RK approximation P Let us consider a solution domain Ω discretized by a set of nodes, {X I } NI =1 . The standard MLS/RK shape function of node I can be written as

N I (x) = pT (X I )w I (x)α(x),

(14)

where p(x) and w I (x) are basis and weight functions given at node I . The unknown vector α(x) is determined by a consistency or reproducibility condition, or p (x) =

NP ∑

p (X I ) N I (x) .

(15)

I =1

Substituting Eq. (14) into Eq. (15) leads to A (x) α (x) = p (x) ,

(16)

where A (x) =

NP ∑ I =1

p (X I ) pT (X I ) w I (x) .

(17)

4

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Here, the vector α (x) can be evaluated by solving Eq. (16), and then the nodal MLS shape function N I (x) is determined. By taking derivatives of Eq. (14), nodal first-order MLS derivatives N I,i (x) are obtained as [ ] (18) N I,i (x) = pT (X I ) w I,i (x) α (x) + w I (x) α,i (x) , where α, i (x) is solved by differentiating Eq. (15), to give A (x) α, i (x) = p, i (x) − A, i (x) α (x) ,

(19)

with A,i (x) =

NP ∑

p (X I ) pT (X I ) w I,i (x) .

(20)

I

By further taking derivatives of Eq. (18), nodal second-order MLS derivatives are obtained as [ N I,i j (x) = pT (X I ) w I,i j (x) α (x) + w I,i (x) ] α, j (x) +w I, j (x) α,i (x) + w I (x) α,i j (x) ,

(21)

where α, i j (x) is solved by differentiating Eq. (19), leading to A (x) α, i j (x) = p, i j (x) − A, i (x) α, j (x) − A, j (x) α,i (x) − A, i j (x) α (x) ,

(22)

with A,i j (x) =

NP ∑

p (X I ) pT (X I ) w I,i j (x) .

(23)

I

3.2. Spatial discretization In terms of MLS shape functions, the displacement field can be approximated as uh (x) = N (x) U =

NP ∑

N I (x) U I ,

(24)

I

[ where U is a nodal displacement parameter vector, and in the case of a two-dimensional problem, U I = u I N (x) is a matrix of nodal MLS shape functions given by [ ] N (x) = N1 (x) N2 (x) · · · N N P (x) ,

vI

]T

.

(25)

where N I (x) = N I (x) I2 , with I2 a two-by-two identity matrix. Substituting Eq. (24) into Eq. (13), a system of discrete equations can be written as (K1 + K2 ) U = f, where stiffness matrices and the external force are given by ∫ ∫ ( T2 ) T K1 = B DBdΩ , K2 = Bx l DBx + BTy l 2 DB y dΩ , Ω ∫ ∫ Ω ∫ ∂NT T T rdΓ , f= N bdΩ + N tdΓ + Γ2 ∂n Ω Γ1

(26)

(27) (28)

with 2

∂ (·) ∑ ∂ (·) = ni , ∂n ∂ xi i=1

(29)

where D is a constitutive matrix based on Lame parameters, and strain and strain gradient matrices are defined as ⎡ ⎤ N I,x 0 [ ] N I,y ⎦ , B = B1 B2 · · · B N P and BI = ⎣ 0 (30) N I,y N I,x

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

⎡

[ Bx = Bx1

Bx2

···

Bx N P

]

and

[ B y = B y1

B y2

···

By N P

]

and

N I,x x Bx I = ⎣ 0 N I,x y ⎡ N I,x y By I = ⎣ 0 N I,yy

5

⎤

0 N I,x y ⎦ , N I,x x ⎤ 0 N I,yy ⎦ . N I,x y

(31)

(32)

Here it is worth noting that, according to Eq. (27), the domain integration is required to evaluate stiffness matrices K1 and K2 . Multiple integration points per a background cell are usually employed to achieve stable numerical results because nodal MLS derivatives are non-polynomial rational functions [42]. In addition, owing to product terms of nodal second-order MLS derivatives in K2 , the domain integration is more difficult than that in a classical elasticity problem. 4. Integration constraints and consistent integration schemes 4.1. Integration constraint for strain gradient elasticity 4.1.1. One-dimensional formulation In the case of strain gradient elasticity, the one-dimensional governing equation can be written as ( 2 4 ) d u 2d u E − l + b = 0 on Ω , dx2 dx4

(33)

u = u on Γ1 , du = w on Γ2 , dx ) ( d 3u du E − l2 3 nx = t dx dx

(34) on Γ3 ,

d 2u n x = r on Γ4 . dx2 Here, E is the elastic modulus. The Galerkin weak form of Eqs. (33) and (34) is ∫ ∫ ∫ ∫ ∫ dδv du dδv d 2 δv 2 d 2 u δvtdΓ δvbdΩ + + r dΓ . El dΩ = E dΩ + 2 2 dx dx Ω dx Γ3 Γ4 d x Ω Ω dx El 2

(35)

As shown in [45], to derive the integration constraint, it is usually supposed that the solution to Eqs. (33) and (34) is a complete monomial, or u = aT p (x) ,

(36)

where p (x) is an arbitrary-order polynomial base and a is a constant vector. Substituting Eq. (36) into Eqs. (33) and (34) leads to ( ) b = −EaT p,x x − l 2 p,x x x x , (37) ( ) T 2 t = Ea p,x − l p,x x x n x on Γ3 , (38) r = El 2 aT p,x x n x on Γ4 . Then, substituting Eqs. (37) and (38) into Eq. (35), we have ∫ ∫ dδv T d 2 δv 2 T Ea p,x dΩ + El a p,x x dΩ 2 Ω ∫d x Ω dx ∫ ( ) ( ) =− δv EaT p,x x − l 2 p,x x x x dΩ + δv EaT p,x − l 2 p,x x x n x dΓ Γ3 ∫ Ω dδv 2 T + El a p,x x n x dΓ . Γ4 d x

(39)

6

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

∑ Here, all of E, l a and δu I are arbitrary and by using δv = NI P N I (x), we thus obtain ∫ ∫ ∫ N I,x (x) p,x dΩ = N I (x) p,x n x dΓ − N I (x) p,x x dΩ , Ω

∫

Γ3

N I,x x (x) p,x x dΩ =

∫

Ω

(40)

Ω

N I,x (x) p,x x n x dΓ −

Γ∫ 4

+

∫

N I (x) p,x x x n x dΓ

Γ3

(41)

N I (x) p,x x x x dΩ .

Ω

Given that the essential boundary conditions are not satisfied, Eqs. (40) and (41) can be rewritten as ∫ ∫ ∫ N I,x (x) p,x dΩ = N I (x) p,x n x dΓ − N I (x) p,x x dΩ , Ω

∫

Γ

N I,x x (x) p,x x dΩ =

Ω

(42)

Ω

∫

N I,x (x) p,x x n x dΓ − Γ∫ + N I (x) p,x x x x dΩ ,

∫

N I (x) p,x x x n x dΓ

Γ

(43)

Ω

which are the arbitrary-order integration constraint conditions of a one-dimensional strain gradient problem. In a quadratic MLS/RK approximation, we have [ ]T 1 2 , (44) p (x) = 1 x x 2 and Eqs. (42) and (43) thus degenerate into { } { } { } ∫ ∫ ∫ 0 1 1 N I (x) dΩ , (45) N I (x) n x dΓ − N I,x (x) dΩ = 1 x x Ω Γ Ω ∫

N I,x x (x) dΩ =

∫

N I,x (x) n x dΓ ,

(46)

Γ

Ω

which are the quadratic integration constraint conditions. In a cubic MLS/RK approximation, we have [ ]T 1 2 1 3 p (x) = 1 x . x x 2 6 Similarly, Eqs. (42) and (43) reduce to ⎧ ⎫ ⎧ ⎫ ⎧ ⎫ ∫ ∫ ∫ ⎨1 ⎬ ⎨1 ⎬ ⎨0 ⎬ N I (x) x n x dΓ − N I (x) 1 dΩ , N I,x (x) x dΩ = ⎩ 2⎭ ⎩ 2⎭ ⎩ ⎭ Ω Γ Ω 2x x x ∫ Ω

(47)

(48)

{ } { } { } ∫ ∫ 1 1 0 N I,x x (x) dΩ = N I,x (x) n x dΓ − N I (x) n dΓ , x x 1 x Γ Γ

(49)

which are the cubic integration constraint conditions in strain gradient analysis. 4.1.2. Two-dimensional formulation If the essential boundaries are not met, Eq. (13) can be further rewritten as ∫ ∫ ∫ ( ) ( ) ( ) σi j δεi j + τi jk δκi jk dΩ = ∂i ∂ j τi jk δu k dΩ − ∂i σi j δu j dΩ Ω Ω∫ ∫ Ω ∫ ( ) + n i σi j δu j dΓ + n i τi jk ∂ j δu k dΓ − n j ∂i τi jk δu k dΓ . Γ

Γ

(50)

Γ

Here, let us suppose that its exact solution is a polynomial of approximate order, or u i = aiT p (x) ,

(51)

7

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

where ai denotes the ith direction polynomial coefficient. Strains are then given by ] 1[ T ai p, j (x) + aTj p,i (x) = biTj q (x) , (52) 2 ] 1[ T κi jk = ai p, jk (x) + aTj p,ik (x) = biTj q,k (x) = biTj q (x) , (53) 2 where q (x) is the basis function that is one order lower than p (x) and while bi j is the corresponding coefficient vector. Similarly, owing to arbitrariness of material parameters and δu I , substituting Eqs. (52) and (53) into Eq. (50) leads to ∫ ∫ ∫ N I,i (x) q (x) dΩ = N I (x) q (x) n i dΓ − N I (x) q,i (x) dΩ (54) εi j =

Ω

Γ

Ω

and ∫ Ω

N I,i j q (x) dΩ

∫

=

] 1[ N I,i q (x) n j + N I, j q (x) n i dΓ Γ∫2 ] 1 [ − N I q, j (x) n i + q,i (x) n j dΓ ∫Γ 2 N I q,i j (x) dΩ . −

(55)

Ω

These are the derived integration constraints for two-dimensional strain gradient elasticity, which express the relationships between the nodal shape functions and their derivatives according to the divergence theory. Therefore, Eqs. (54) and (55) can be automatically satisfied if a numerical integration rule possesses relevant order exactness. 4.2. Consistent integration schemes 4.2.1. One-dimensional case [ Through the numerical integration, discrete forms of Eqs. (42) and (43) in a background integration cell xs are obtained as ng ∑

]

ng

∑ ( ) ( ) ( ) ( ) [ ]⏐x Wg N I x g p,x x x g , Wg N I,x x g p,x x g = N I (x) p,x n x ⏐xts −

(56)

g=1

g=1 ng

∑

xt

( ) ( ) Wg N I,x x x g p,x x x g

=

[

[ ]⏐x ]⏐x N I,x (x) p,x x n x ⏐xts − N I (x) p,x x x n x ⏐xts

g=1

+

ng ∑

(57) ( ) ( ) Wg N I x g p,x x x x x g ,

g=1

where Wg is the domain integration weight, x g is the coordinate integration point, and ng is the number of integration points. For a quadratic basis, Eq. (56) becomes the following two equations, that is ng ∑

( ) Wg N I,x x g = [N I n x ]|xxts ,

g=1 ng ∑

ng

(58)

∑ ( ) Wg N I,x x g x g = [N I xn x ]|xxts − Wg N I ,

g=1

g=1

while Eq. (57) remains one equation, or ng ∑

[ ]⏐x Wg N I,x x = N I,x n x ⏐x1 . 0

g=1

(59)

8

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Fig. 1. Illustration of integration schemes consistent with (a) the quadratic basis for K1 , (b) the quadratic basis for K2 , (c) the cubic basis for K1 , and (d) the cubic basis for K2 , including vertices (squares) of a background cell and domain integration points (circles).

Using nodal smoothed derivatives instead of standard ones on the left sides of Eqs. (58) and (59), we have ng ∑

( ) Wg N˜ I,x x g = [N I n x ]|xxts ,

g=1 ng ∑

ng

(60)

∑ ( ) Wg N˜ I,x x g x g = [N I xn x ]|xxts − Wg N I ,

g=1

g=1

and ng ∑

[ ]⏐x Wg N˜ I,x x = N I,x n x ⏐xts .

(61)

g=1

In solving Eq. (60), a two-point integration scheme (see Fig. 1(a)) is used to integrate K1 . Here, nodal smoothed derivatives are done according to [ ] [ ][ ] W1 W2 N˜ I,x (x1 ) N I (xt ) + N I (xs ) = . (62) W1 x1 W2 x2 N˜ I,x (x2 ) xt N I (xt ) + xs N I (xs ) − W1 N I (x1 ) − W2 N I (x2 ) To integrate K2 , a one-point integration scheme (see Fig. 1(b)) is designed and nodal second-order smoothed derivatives are given by N I,x (xt ) + N I,x (xs ) N I,x (xt ) + N I,x (xs ) N˜ I,x x (x1 ) = = . W1 xt − xs

(63)

For a cubic basis, Eqs. (60) and (61) become ng ∑

( ) Wg N˜ I,x x g = [N I n x ]|xxts ,

g=1 ng ∑

ng

∑ ( ) Wg N I , Wg N˜ I,x x g x g = [N I xn x ]|xxts −

g=1 ng

∑

(64)

g=1 ng

( ) [ ]⏐xt ∑ Wg N˜ I,x x g x g2 = N I x 2 n x ⏐xs − 2x g Wg N I , g=1

g=1

and ng ∑

[ ]⏐x Wg N˜ I,x x = N I,x n x ⏐xts ,

g=1 ng ∑ g=1

(65) Wg x g N˜ I,x x

]⏐x = N I,x xn x ⏐xts − [N I n x ]|xxts . [

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

9

We therefore use the integration schemes illustrated in Fig. 1(c) and (d) to evaluate smoothed derivatives and obtain ⎡ ⎤⎡ ⎤ N˜ I,x (x1 ) W1 W2 W3 ⎣ W1 x1 W2 x2 W3 x3 ⎦ ⎣ N˜ I,x (x2 )⎦ W x 2 W2 x22 W3 x32 N˜ I,x (x3 ) ⎡1 1 ⎤ (66) N I (xt ) + N I (xs ) ⎦, xt N I (xt ) + xs N I (xs ) − W1 N I (x1 ) − W2 N I (x2 ) − W3 N I (x3 ) =⎣ xt2 N I (xt ) + xs2 N I (xs ) − 2x1 W1 N I (x1 ) − 2x2 W2 N I (x2 ) − 2x3 W3 N I (x3 ) and [

W1 W1 x 1

W1 W2 x 2

][

] [ ] N˜ I,x x (x1 ) N I,x (xt ) + N I,x (xs ) . = xt N I,x (xt ) + xs N I,x (xs ) − N I (xt ) − N I (xs ) N˜ I,x x (x2 )

(67)

4.2.2. Two-dimensional case The triangular background integration cells are employed in two-dimensional analysis. With nodal smoothed derivatives instead of standard ones, the discrete integration forms of Eqs. (54) and (55) in a background integration cell can be written as ng ∑

ng

ns 3 ∑ ∑ ( ) ( ) ∑ ( ) ( ) wb N I (xb ) q (xb ) n ie − Wg N˜ I,i xg q xg = Wg N I xg q,i xg ,

g=1

e=1 b=1

(68)

g=1

and ng ∑

3 ∑ ns ∑ ] ( ) ( ) 1 [ wb N I,i (xb ) q (xb ) n ej + N I, j (xb ) q (xb ) n ie Wg N˜ I,i j xg q xg = 2 e=1 b=1 g=1

−

ns 3 ∑ ∑ 1 e=1 b=1 ng

−

∑

2

[ ] wb N I (xb ) q,i (xb ) n ej + N I (xb ) q, j (xb ) n ie

(69)

( ) ( ) Wg N I xg q,i j xg ,

g=1

where xb denotes the position of a boundary integration point with weight wb and ns is the number of integration points on each boundary. Taking the x-direction smoothed derivatives as an example, Eq. (68) becomes the following three equations, that is ng ∑ g=1 ng ∑

3 ∑ ns ( ) ∑ Wg N˜ I,x xg = wb N I (xb ) n ie , e=1 b=1 ng

3 ∑ ns ∑ ∑ ( ) ( ) Wg N˜ I,x xg x g = wb N I (xb ) x g n ie − Wg N I xg ,

g=1 ng

∑ g=1

( ) Wg N˜ I,x xg yg =

e=1 b=1 3 ∑ ns ∑

(70)

g=1

wb N I (xb ) yg n ie .

e=1 b=1

Eq. (69) remains one equation, or ng ∑

3 ∑ ns ( ) ∑ Wg N˜ I,x x xg = wb N I,x (xb ) n ex .

g=1

e=1 b=1

(71)

10

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Fig. 2. Illustration of integration schemes consistent with the quadratic basis based on a triangular background integration cell for (a) K1 and (b) K2 , where squares denote vertices of background cells, circles are domain integration points, and triangles are evaluating points of boundary integration.

Here, a three-point integration scheme is suggested (see Fig. 2(a)) for calculating K1 and a one-point integration scheme (see Fig. 2(b)) for K2 . Therefore, the x-direction smoothed derivatives are obtained by solving ⎡ ⎤ 3 ∑ 2 ∑ e ⎢ ⎥ wb N I (xb ) n i ⎢ ⎥ ⎥ e=1 b=1 ⎡ ⎤⎡ ⎤ ⎢ ⎥ ⎢∑ 3 ∑ 2 3 N˜ I,x (x1 ) W1 W2 W3 ∑ ( ) ⎢ ⎥ e ⎣W1 x1 W2 x1 W3 x1 ⎦ ⎣ N˜ I,x (x2 )⎦ = ⎢ w N W N x x n − (x ) (72) b I b g I g ⎥ g i ⎢ ⎥, ⎢ ⎥ ˜ e=1 b=1 g=1 W1 y1 W2 y1 W3 y1 N I,x (x3 ) ⎢ ⎥ 3 ∑ 2 ⎢ ⎥ ∑ ⎣ ⎦ w N (x ) y n e b

I

b

g i

e=1 b=1

and the xx-direction smoothed derivatives are determined according to 3 2 1 ∑∑ wb N I,x (xb ) n ex . N˜ I,x x (x1 ) = W1 e=1 b=1

(73)

The nodal smoothed derivatives at other directions can be evaluated in a similar way. Those derivatives can then be introduced instead of standard ones to compute stiffness matrices. Here it is worth noting that, in addition to the integration constraint conditions, nodal smoothed derivatives also satisfy differential consistency conditions, or NP ∑

pk (X I ) N˜ I,i (x) = pk,i (x)

(74)

pk (X I ) N˜ I,i j (x) = pk,i j (x) .

(75)

I =1

and NP ∑ I =1

Proofs of Eqs. (74) and (75) are given in Appendix. 5. Numerical results In the following case studies, we denote a standard GI with n points in a cell as GI-n. The essential boundary conditions of fixed displacements and the normal gradients are respectively applied by coupling with a weight function [48] and a penalty method. For simplicity, the units of quantities are omitted in discussion.

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

11

Table 1 Results of quadratic and cubic patch tests. Error

L2 H1 H2

Quadratic patch test

Cubic patch test

GI-2

GI-10

CI

GI-3

GI-10

CI

3.79 × 10−2 5.29 × 10−2 2.69 × 10−1

1.21 × 10−2 4.41 × 10−2 3.22 × 10−2

2.94 × 10−15 1.98 × 10−15 1.30 × 10−13

2.21 × 10−3 9.70 × 10−3 8.38 × 10−2

1.36 × 10−3 5.56 × 10−3 5.93 × 10−2

8.78 × 10−15 1.13 × 10−13 1.66 × 10−12

Fig. 3. Schematic representation of a bar under tension.

5.1. One-dimensional patch tests The patch tests are investigated on a domain Ω = (0, L) ⊂ R1 with 10 irregular nodes. In quadratic and cubic patch tests, the displacement solutions are u = 1 + x + x 2,

(76)

u = x 3.

(77)

Here, all the boundaries are fixed by displacements and double tractions with E = 1 × 106 , L = 10, and l = 0.1. Quadratic and cubic bases are respectively employed in quadratic and cubic patch tests. As listed in Table 1, the consistent integration (CI) passes patch tests whereas GI fails. It is thus confirmed that CI possesses relevant order exactness in strain gradient analysis. 5.2. Bar under tension Let us consider a bar with L = 10, as illustrated in Fig. 3. The boundary conditions at the left end are u (0) = 0, r (0) = 0,

(78)

while those on the right end are u ,x (L) = 0, τ (L) = P. The exact solution [49] can be obtained as [ ] P sinh (x/l) x −l u= . EA cosh (L/l)

(79)

(80)

Here, material parameters are chosen as E = 1 × 106 , L = 10, and l = 0.1L. Four irregular grids with 10, 20, 40 and 100 nodes are applied to calculate the convergence rates by using a least-squares method. Fig. 4 shows that, in the case of a quadratic basis function, the accuracy and convergence of GI-10 are better than those of GI-2. This implies that use of more evaluating points improves the numerical performance in a Gaussian type integration method. The L2 and H1 accuracies of CI are almost one order higher than those of GI-2 with the finest grid. For a quadratic approximation, the theoretical convergence rates of L 2 , H1 and H2 errors in strain gradient analysis are 2, 2 and 1, respectively [49]. The convergence rates of L 2 , H1 and H2 errors of CI are 2.35, 1.93 and 1.02, which reach and even exceed the theoretical convergence rates. In the case of a cubic approximation, the theoretical convergence rates of L 2 , H1 and H2 errors are 4, 3 and 2, respectively [49]. The convergence rates of L 2 , H1 and H2 errors for CI are 3.15, 2.71 and 1.9 (see Fig. 5), which leads to sub-optimal convergence. This might be due to a cubic meshfree approximation that is not as good as that of a quadratic meshfree approximation. The worse convergence of a cubic approximation was reported for

12

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Fig. 4. Convergence results of a bar under tension by a quadratic meshfree method: (a) L2 error, (b) H1 error, and (c) H2 error.

conventional GI in classical elasticity analysis [43,47]. Nevertheless, as shown in Figs. 4 and 5, the accuracy of CI in a cubic meshfree method is almost one order higher than that of CI in a quadratic meshfree method. Additionally, CI performs much better than GI in terms of accuracy and convergence.

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

13

Fig. 5. Convergence results of a bar under tension by a cubic meshfree method: (a) L2 error, (b) H1 error, and (c) H2 error.

Fig. 6 shows that the distributions of u ,x x and u ,x x x along the bar obtained by CI are in agreement with the theoretical solution, while GI-10 has slight oscillation in the field of u ,x x and serious oscillation in the field of u ,x x x even with more integration points.

14

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Fig. 6. Field distributions of (a) u ,x x and (b) u ,x x x along the bar.

5.3. Thick hollow cylinder As illustrated in Fig. 7, let us consider a thick hollow cylinder in a plane strain state. Owing to twofold symmetries, only a quarter of the model is analyzed with an irregular nodal distribution of 3649 nodes in discretization. The elastic modulus and Poisson’s ratio are E = 107 and υ = 0.3. The internal and external radii are R1 = 1 and R2 = 2, respectively. Here, a quadratic MLS/RK basis is employed. As shown in Fig. 8, the numerical results with different internal material lengths (l = 0, 0.2, 0.4, 0.6, 0.8 and 1.0) are consistent with the theoretical solutions [10], which validates the robustness of CI in a two-dimensional problem. In contrast with numerical oscillations of the κ111 fields obtained by GI-3 and GI-16, CI has a much smoother contour (see Fig. 9). 5.4. Manufactured problem Similarly, we consider a square plate on a domain Ω = (0, 2)×(0, 2) ⊂ R2 with material parameters of E = 107 , υ = 0.3, and l = 0.2. Further, let us assume that the plate is manufactured under a body force [ ] E 12x 2 − 24x + 8 − 24g 2 . (81) b=− 1 − υ 2 12y 2 − 24y + 8 − 24g 2 The displacement solution can then be obtained as u = (x 2 − 2x)2 , v = (y 2 − 2y)2 .

(82)

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

15

Fig. 7. Illustration of (a) the solution domain in a thick hollow cylinder and (b) the nodal discretization.

Fig. 8. (a) Radial displacement u r and (b) strain εrr for different g values, where symbols denote the numerical results obtained by using CI and solid lines are the theoretical solutions.

Here, all the boundaries are fixed, that is u i = 0, wi = 0.

(83)

The five grids with 9 × 9, 17 × 17, 21 × 21, 33 × 33, and 81 × 81 irregular nodal configurations are employed for the convergence study with quadratic MLS/RK basis functions. As shown in Fig. 10, CI has much better

16

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Fig. 9. Fields κ111 in a thick hollow cylinder obtained by (a) GI-3, (b) GI-16, and (c) CI.

Fig. 10. Convergence results of (a) L2 error and (b) H1 error in a manufactured problem.

numerical performance than does the standard GI. Based on convergence and accuracy, GI-16 performs a little better than GI-3 but requires much more computational time (see Fig. 11). CI is more efficient than GI. It is also obvious to see in Fig. 12 that the κ111 fields obtained by GI-3 and GI-16 result in serious numerical oscillation. In contrast, CI leads to a smoothed contour even under a non-constant body force. 6. Conclusions In summary, the arbitrary-order integration constraint conditions for strain gradient elasticity have been investigated by using a meshfree Galerkin method. The stable and efficient integration schemes consistent with

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

17

Fig. 11. Computational efficiency results of (a) L2 error and (b) H1 error in a manufactured problem.

Fig. 12. Fields κ111 in a manufactured problem obtained by (a) GI-3, (b) GI-16, and (c) CI.

approximations have been proposed to meet these integration constraint conditions. The main conclusions of this study can be summarized as follows. (1) The integration constraint conditions for strain gradient elasticity are derived by using the divergence theory to establish a discrete formulation of nodal shape functions and their first- and second-order derivatives. (2) The number increase of evaluating points in GI can improve the accuracy and convergence rate but with higher computational cost. (3) In contrast to GI, CI is able to pass relevant patch tests with much better numerical performance based on accuracy, convergence and efficiency.

18

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

Acknowledgments This work has been supported by the National Natural Science Foundation of China (Nos. 11572289, 11702251 and 11702252). Appendix. Proofs of differential consistency conditions (Eqs. (74) and (75)) Multiplying both sides of Eq. (68) by pk (X I ) and summing over I lead to ng ∑

NP ( )∑ ( ) Wg q xg N˜ I,i xg pk (X I )

g=1 3 ∑ ns ∑

=

=

e=1 b=1 3 ∑ ns ∑

I =1

wb q (xb ) n ie

NP ∑

N I (xb ) pk (X I ) −

I =1

ng ∑

NP ( )∑ ( ) Wg q,i xg N I xg pk (X I )

g=1

wb q (xb ) n ie pk (xb ) −

ng ∑

e=1 b=1

(A.1)

I =1

( ) ( ) Wg q,i xg pk xg .

g=1

According to the consistency of the MLS approximation or Eq. (15), we have ng ∑

NP ( )∑ ( ) Wg q xg N˜ I,i xg pk (X I )

g=1 3 ∑ ns ∑

=

I =1

wb q (xb ) n ie pk

(xb ) −

ng ∑

e=1 b=1

(A.2) ( ) ( ) Wg q,i xg pk xg ,

g=1

where the right side is a discrete numerical integration form of obtain 3 ∑ ns ∑

wb q (xb ) n ie pk (xb ) −

e=1 ∫ b=1

ng ∑

∫ Γ

q (x) pk (x) n i dΓ − q,i (x) pk (x) n i dΩ . We thus

( ) ( ) Wg q,i xg pk xg (A.3)

g=1

q (x) pk (x) n i dΓ − q,i (x) pk (x) n i dΩ .

= Γ

Applying the divergence theorem to the right side leads to ∫ ng NP ∑ ( )∑ ( ) pk,i (x)q (x) dΩ . Wg q xg N˜ I,i xg pk (X I ) = g=1

I =1

(A.4)

Ω

The right side of Eq. (A.2) can be exactly integrated and obtained as ng ∑

ng

NP ∑ ( )∑ ( ) ( ) ( ) Wg q xg N˜ I,i xg pk (X I ) = Wg q xg pk,i xg .

g=1

I =1

Then, Eq. (A.2) can be rewritten as [NP ] ng ∑ ( ) ∑ ( ) ( ) ˜ Wg q xg N I,i xg pk (X I ) − pk,i xg = 0, g=1

(A.5)

g=1

(A.6)

I =1

and its matrix–vector form is Qdk = 0,

(A.7)

where [ Q = W1 q (x1 )

W2 q (x2 ) · · ·

( )] Wng q xng ,

(A.8)

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

⎤ ˜ N I,i (x1 ) pk (X I ) − pk,i (x1 ) ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ I =1 ⎥ ⎢ NP ⎥ ⎢ ∑ ⎢ N˜ I,i (x2 ) pk (X I ) − pk,i (x2 ) ⎥ ⎥. dk = ⎢ ⎥ ⎢ I =1 ⎥ ⎢ ⎥ ⎢ ··· ⎥ ⎢ NP ⎢∑ ( ) ( )⎥ ⎦ ⎣ N˜ I,i xng pk (X I ) − pk,i xng ⎡

19

NP ∑

(A.9)

I =1

In the integration schemes, Q is a nonsingular square matrix with only one solution, which yields the differential consistency condition or Eq. (74). Similarly, to prove Eq. (75), we multiply both sides of Eq. (69) by pk (X I ) and sum over I , and obtain ng ∑

NP ( )∑ ( ) Wg q xg N˜ I,i j xg pk (X I ) I =1

g=1

=

ns 3 ∑ ∑

1 wb 2

e=1 b=1 ns 3 ∑ ∑

−

e=1 b=1

−

ng ∑

[NP ∑

1 wb 2

N I, j (xb )

I =1 NP ∑

pk (X I ) N I (xb ) q,i (xb ) n ej +

I =1

] pk (X I ) N I (xb ) q, j (xb ) n ie

I =1

] 1 [ wb pk,i (xb ) q (xb ) n ej + pk, j (xb ) q (xb ) n ie 2

e=1 b=1 3 ∑ ns ∑

−

e=1 b=1

−

] pk (X I ) q (xb ) n ie

I =1

3 ∑ ns ∑

ng ∑

+

NP ∑

NP ( )∑ ( ) N I xg pk (X I ) Wg q,i j xg

g=1

=

N I,i (xb )

I =1 [NP ∑

pk (X I ) q (xb ) n ej

] 1 [ wb pk (xb ) q,i (xb ) n ej + pk (xb ) q, j (xb ) n ie 2

( ) ( ) Wg q,i j xg pk xg

g=1

∫ =

] 1[ pk,i (x) q (x) n j + pk, j (x) q (x) n i dΓ Γ∫2 [ ] 1 pk (x) q, j (x) n i + q,i (x) n j dΓ − 2 ∫Γ − Ω

(A.10)

pk (x) q,i j (x) dΩ .

Applying the divergence theorem to the right side twice, we have ∫ ng NP ∑ ( )∑ ( ) ˜ Wg q xg N I,i j xg pk (X I ) = pk,i j (x)q (x) dΩ . g=1

I =1

Then, the right side can be exactly integrated and further written as [NP ] ng ∑ ( ) ∑ ( ) ( ) Wg q xg N˜ I,i j xg pk (X I ) − pk,i j xg = 0, g=1

(A.11)

Ω

(A.12)

I =1

and its matrix–vector form is Qdk = 0,

(A.13)

20

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

where [ Q = W1 q (x1 ) W2 q (x2 ) · · · ⎡ NP ∑ N˜ I,i j (x1 ) pk (X I ) − ⎢ ⎢ ⎢ I =1 ⎢ NP ⎢ ∑ ⎢ N˜ I,i j (x2 ) pk (X I ) − dk = ⎢ ⎢ I =1 ⎢ ⎢ ··· ⎢ NP ⎢∑ ( ) ⎣ N˜ I,i j xng pk (X I ) −

( )] Wng q xng , ⎤ pk,i j (x1 ) ⎥ ⎥ ⎥ ⎥ ⎥ pk,i j (x2 ) ⎥ ⎥. ⎥ ⎥ ⎥ ⎥ ( )⎥ ⎦ pk,i j xng

(A.14)

(A.15)

I =1

The nonsingular property of square matrix Q leads to only one solution; i.e., the differential consistency condition or Eq. (75). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

E.C. Aifantis, On the role of gradients in the localization of deformation and fracture, Internat. J. Engrg. Sci. 30 (1992) 1279–1299. B. Altan, E.C. Aifantis, On some aspects in the special theory of gradients elasticity, J. Mech. Behav. Mater. 8 (1997) 231–282. E.C. Aifantis, Strain gradient interpretation of size effects, Int. J. Fract. 95 (1999) 299–314. H.T. Zhu, H.M. Zbib, E.C. Aifantis, Strain gradients and continuum modeling of size effect in metal matrix composites, Acta Mech. 121 (1997) 165–176. E.C. Aifantis, Higher order gradients and size effects, in: A. Carpinteri (Ed.), Size-Scale Effects in the Failure Mechanisms of Materials and Structures, Chapman & Hall, London, 1996. E.C. Aifantis, On the gradient approach - Relation to Eringen’s nonlocal theory, Internat. J. Engrg. Sci. 49 (2011) 1367–1377. R.D. Mindlin, Micro-structure in linear elasticity, Arch. Ration. Mech. Anal. 16 (1964) 51–78. A.C. Eringen, E.S. Suhubi, Nonlinear theory of simple micro-elastic solids-I, Internat. J. Engrg. Sci. 2 (1964) 189–203. E.C. Aifantis, Internal length gradient (ILG) material mechanics across scales and disciplines, Adv. Appl. Mech. 49 (2016) 1–110. A. Zervos, S.A. Papanicolopulos, I. Vardoulakis, Two finite element discretizations for gradient elasticity, J. Eng. Mech. 13 (2009) 203–213. S.A. Papanicolopulos, A. Zervos, I. Vardoulakis, A three dimensional C1 finite element for gradient elasticity, Internat. J. Numer. Methods Engrg. 77 (2009) 1396–1415. E. Amanatidou, N. Aravas, Mixed finite element formulations of strain-gradient elasticity problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 1723–1751. V. Phunpeng, P.M. Baiz, Mixed finite element formulations for strain-gradient elasticity problems using the FEniCS environment, Finite Elem. Ana. Des. 96 (2015) 23–40. M. Malagu, E. Benvenuti, C.A. Duarte, A. Simone, One-dimensional nonlocal and gradient elasticity: Assessment of high order approximation schemes, Comput. Methods Appl. Mech. Engrg. 275 (2014) 138–158. C.Q. Ru, E.C. Aifantis, A simple approach to solve boundary-value problems in gradient elasticity, Acta Mech. 101 (1993) 59–68. H. Askes, E.C. Aifantis, Gradient elasticity theories in statics and dynamics - A unification of approaches, Int. J. Fract. 139 (2006) 297–304. H. Askes, E.C. Aifantis, Gradient elasticity in statics and dynamics: An overview of formulations, length scale identification procedures, finite element implementations and new results, Int. J. Solids Struct. 48 (2011) 1962–1990. T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Internat. J. Numer. Methods Engrg. 37 (1994) 229–256. W.K. Liu, S. Jun, Y.F. Zhang, Reproducing kernel particle methods, Int. J. Numer. Methods Fluids 20 (1994) 1081–1106. S.N. Atluri, T. Zhu, A new Meshless Local Petrov–Galerkin (MLPG) approach in computational mechanics, Comput. Mech. 22 (1998) 117–127. J.S. Chen, M. Hillman, S.W. Chi, Meshfree methods: Progress made after 20 years, J. Eng. Mech. 143 (2017) 04017001. N. Sukumar, R.W. Wright, Overview and construction of meshfree basis functions: from moving least squares to entropy approximants, Internat. J. Numer. Methods Engrg. 70 (2007) 181–205. D.D. Wang, J.R. Wang, J.C. Wu, Superconvergent gradient smoothing meshfree collocation method, Comput. Methods Appl. Mech. Engrg. 340 (2018) 728–766. H. Askes, E.C. Aifantis, Numerical modeling of size effects with gradient elasticity-Formulation, meshless discretization and examples, Int. J. Frcture 117 (2002) 347–358. C. Sansour, S. Skatulla, A strain gradient generalized continuum approach for modelling elastic scale effects, Comput. Methods Appl. Mech. Engrg. 198 (2009) 1401–1412. I. Tsagrakis, E.C. Aifantis, Element-free Galerkin implementation of gradient plasticity, Part I: Formulation and application to 1D strain localization, J. Mech. Behav. Mat. 14 (2003) 199–231. I. Tsagrakis, E.C. Aifantis, Element-free Galerkin implementation of gradient plasticity, Part II: Applications to 2D strain localization and size effects, J. Mech. Behav. Mat. 14 (2003) 233–253.

B. Wang, C. Lu, C. Fan et al. / Computer Methods in Applied Mechanics and Engineering 357 (2019) 112601

21

[28] I. Tsagrakis, E.C. Aifantis, Recent developments in gradient plasticity. Part I: Formulation and size effects, J. Eng. Mat. Tech. 124 (2002) 352–357. [29] J. Dolbow, T. Belytschko, Numerical integration of the Galerkin weak form in meshfree methods, Comput. Mech. 23 (3) (1999) 219–230. [30] J.S. Chen, C.T. Wu, S. Yoon, Y. You, A stabilized conforming nodal integration for Galerkin mesh-free methods, Internat. J. Numer. Methods Engrg. 50 (2001) 435–466. [31] I. Babuska, U. Banerjee, J.E. Osborn, Q.L. Li, Quadrature for meshless methods, Internat. J. Numer. Methods Engrg. 76 (2008) 1434–1470. [32] S. Beissel, T. Belytschko, Nodal integration of the element-free Galerkin method, Comput. Methods Appl. Mech. Engrg. 139 (1996) 49–74. [33] G.R. Liu, G.Y. Zhang, Y.Y. Wang, Z.H. Zhong, G.Y. Li, X. Han, A nodal integration technique for meshfree radial point interpolation method (NI-RPIM), Int. J. Solids Struct. 44 (2007) 3840–3860. [34] K.C. Kwon, S.H. Park, S.K. Youn, The support integration scheme in the least-square meshfree method, Finite Elem. Anal. Des. 43 (2006) 127–144. [35] Y. Liu, T. Belytschko, A new support integration scheme for the weakform in mesh-free methods, Internat. J. Numer. Methods Engrg. 82 (2010) 699–725. [36] C.T. Dyka, P.W. Randles, R.P. Ingel, Stress points for tension instability in SPH, Internat. J. Numer. Methods Engrg. 40 (1997) 2325–2341. [37] T.P. Fries, T. Belytschko, Convergence and stabilization of stress-point integration in mesh-free and particle methods, Internat. J. Numer. Methods Engrg. 74 (2008) 1067–1087. [38] D.D. Wang, J.S. Chen, Locking-free stabilized conforming nodal integration for meshfree Mindlin-Reissner plate formulation, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1065–1083. [39] D.D. Wang, J.S. Chen, A Hermite reproducing kernel approximation for thin-plate analysis with sub-domain stabilized conforming integration, Internat. J. Numer. Methods Engrg. 74 (2008) 368–390. [40] J.S. Chen, D.D. Wang, A constrained reproducing kernel particle formulation for shear deformable shell in Cartesian coordinate, Internat. J. Numer. Methods Engrg. 68 (2006) 151–172. [41] J.S. Chen, S. Yoon, C.T. Wu, Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods, Internat. J. Numer. Methods Engrg. 53 (2002) 2587–2615. [42] Q.L. Duan, X.K. Li, H.W. Zhang, T. Belytschko, Second-order accurate derivatives and integration schemes for meshfree methods, Internat. J. Numer. Methods Engrg. 92 (2012) 399–424. [43] Q.L. Duan, X. Gao, B.B. Wang, X.K. Li, H.W. Zhang, T. Belytschko, Y.L. Shao, Consistent element-free Galerkin method, Internat. J. Numer. Methods Engrg. 99 (2014) 79–101. [44] Q.L. Duan, X. Gao, B.B. Wang X.K. Li, H.W. Zhang, A four-point integration scheme with quadratic exactness for three-dimensional element-free Galerkin method based on variationally consistent formulation, Comput. Methods Appl. Mech. Engrg. 280 (2014) 84–116. [45] J.S. Chen, M. Hillman, M. Ruter, An arbitrary order variationally consistent integration for Galerkin meshfree methods, Internat. J. Numer. Methods Engrg. 95 (2013) 387–418. [46] D.D. Wang, J. Wu, An efficient nesting sub-domain gradient smoothing integration algorithm with quadratic exactness for Galerkin meshfree methods, Comput. Methods Appl. Mech. Engrg. 298 (2016) 485–519. [47] D.D. Wang, J. Wu, An inherently consistent reproducing kernel gradient smoothing framework toward efficient Galerkin meshfree formulation with explicit quadrature, Comput. Methods Appl. Mech. Engrg. 349 (2019) 628–672. [48] J.S. Chen, W.M. Han, Y. You, X.P. Meng, A reproducing kernel method with interpolation property, Internat. J. Numer. Methods Engrg. 56 (2003) 935–960. [49] J. Niiranen, S. Khakalo, V. Balobanov, A.H. Niemi, Variational formulation and isogeometric analysis for fourth-order boundary value problems of gradient-elastic bar and plane strain/stress problems, Comput. Methods Appl. Mech. Engrg. 308 (2016) 182–211.