Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.019828

ARTICLE

A Fast Element-Free Galerkin Method for 3D Elasticity Problems

Zhijuan Meng1, Yanan Fang1 and Yumin Cheng2,*

1School of Applied Science, Taiyuan University of Science and Technology, Taiyuan, 030024, China
2Shanghai Key Laboratory of Mechanics in Energy Engineering, Shanghai Institute of Applied Mathematics and Mechanics, School of Mechanics and Engineering Science, Shanghai University, Shanghai, 200072, China
*Corresponding Author: Yumin Cheng. Email: ymcheng@shu.edu.cn
Received: 18 October 2021; Accepted: 15 December 2021

1  Introduction

The 3D elasticity problem is one of the typical mechanical problems in engineering. Because of the complexity of this type of problems, it is difficult to get the analytical solution when the problems are defined in complicated domains. Therefore, the study of the numerical method for 3D elasticity has important theoretical and practical significance. At present, lots of methods have been applied to simulate elasticity problems, such as the finite element method (FEM) [1], finite difference method [2], boundary element method (BEM) [3] and meshless method [4–7].

The meshless method has been an important numerical method for scientific and engineering problems [8–13]. When using the meshless method to solve problems, only discrete nodes need to be distributed in the problem domain and on the boundary. Furthermore, different problems can flexibly distribute the quantity of nodes according to the characteristics. Therefore, the meshless method has good adaptability and high calculation accuracy. It is more suitable for solving those complex problems, for example, crack growth, super large transfiguration and high-speed collision.

The IEFG method is one of the important meshless methods. The shape function is constructed using the improved moving least-squares (IMLS) approximation [14,15]. The IMLS approximation is based on the ordinary least-squares method which has high accuracy and has been applied in various fields [16–18]. Furthermore, the basis functions of the IMLS method are weighted orthogonally. Compared it with the classical element-free Galerkin method, the IEFG method has fewer coefficients, and the computational speed is faster. Zhang et al. used the IEFG method to simulate 3D heat conduction problems [19], 3D potential problems [20] and 3D wave propagation [21]. Ma et al. used the IEFG method to simulate groundwater pollution prevention and control and its application in fluid flow [22]. Peng et al. used it to simulate 3D viscoelasticity problems [23]. Zheng et al. used it to simulate diffusional drug release problems [24]. Yu et al. used it to simulate three-dimensional elastoplasticity problems [25]. Liu et al. used it to simulate elastic large deformation problems [26] and inhomogeneous swelling of polymer gels [27]. Cai et al. used it to simulate elastoplasticity large deformation problems [28]. Wu et al. used it to simulate the elasticity [29]. Zou et al. used it to simulate fracture problems of airport pavement [30]. Cheng et al. used it to simulate the unsteady Schrödinger equation [31]. Cheng et al. used it to simulate nonlinear large deformation [32]. However, for 3D problems, the CPU running time of IEFG method is still long. To improve the computational efficiency, the DSM is introduced to complete the 3D problems.

Using the DSM, the 3D problem can be changed into a series of interrelated 2D problems. Li et al. first chose the DSM to solve the 3D linearly elastic shell [33] and incompressible Navier-Stokes equations [34]. Hou et al. solved a 3D elliptic equation appling the dimension splitting algorithm [35]. Ter Maten EJW applied the splitting method to solve fourth order partial differential equations [36]. Bragin et al. applied dimension splitting to analyze the conservation law [37]. Meng et al. presented the dimension splitting element-free Galerkin method to simulate 3D potential problems [38], 3D wave equations [39], 3D advection-diffusion problems [40] and 3D heat conduction problems [41]. Cheng et al. used the dimension splitting and improved complex variable element-free Galerkin method to simulate 3D problems [42–45]. Wu et al. applied the interpolating dimension splitting element-free Galerkin method to simulate 3D heat conduction and 3D potential problems [46,47]. Peng et al. combined the dimension splitting method and reproduce kernel particle method [48–50] to simulate some classical 3D problems [51–53]. The results of references [38–53] show that the computational efficiency can be greatly improved by the DSM.

In this paper, the FEFG method for 3D elasticity is proposed. The key problem is transforming a 3D elasticity into a series of interrelated 2D elasticities using dimension splitting. Then the 2D problems are simulated with the IEFG method and FDM is chosen in the splitting direction. The final discretized system equations of the 3D elasticity problem are obtained. Some specific examples are provided to illustrate the high accuracy and efficiency of the FEFG method for 3D elasticity by comparing the results of the FEFG method with ones of the IEFG method and analytical solutions.

2  The IMLS Approximation

The trial function uh(x) of the moving Least-squares (MLS) method approximation is

uh(x)=∑i=1mpi(x)ai(x)=pT(x)a(x),(x∈Ω),(1)

where pi(x)(i=1,2,…m) are basis functions, m is the number of basic functions. ai(x)(i=1,2,…m) are the coefficients of pi(x).

The form of pi(x) is often given according to the characteristics of the problem.

To minimize the local approximation error, the coefficients ai(x) are calculated by the weighted least square method. Let the functional

J=∑I=1nw(x−xI)[uh(x,xI)−u(xI)]2=∑I=1nw(x−xI)[∑i=1mpi(xI)ai(x)−u(xI)]2=(Pa−u)TW(x)(Pa−u)(2)

where w(x−xI) means the weight function, xI are points whose influence domain include x,

a(x)=(a1(x),a2(x),…,am(x))T,(3)

u=(u1,u2,…,un)T=(u(x1),u(x2),…,u(xn))T,(4)

P=[p1(x1)p2(x1)⋯pm(x1)p1(x2)p2(x2)⋯pm(x2)⋮⋮⋱⋮p1(xn)p2(xn)⋯pm(xn)],(5)

W(x)=[w(x−x1)0⋯00w(x−x2)⋯0⋮⋮⋱⋮00⋯w(x−xn)].(6)

In order to obtain a(x), from

∂J∂a=A(x)a(x)−B(x)u=0,(7)

we have

a(x)=A−1(x)B(x)u,(8)

where

A(x)=PTW(x)P,(9)

B(x)=PTW(x).(10)

Finally, we can obtain the trial function uh(x) as

uh(x)=Φ(x)u=∑I=1nΦI(x)uI,(11)

where Φ(x) represents the shape function. We have

Φ(x)=(Φ1(x),Φ2(x),…,Φn(x))=pT(x)A−1(x)B(x),(12)

ΦI,i(x)=∑j=1m[pj,i(A−1B)jI+pj((A−1),iB+A−1B,i)jI].(13)

For the IMLS method, a set of orthogonal functions is chosen as the basis function. The ill conditioned or singular equations can be avoided. The inverse of the matrix can be obtained directly, which improves the calculation efficiency.

The orthogonal basis functions can be structured by the Schmidt method. For example, basis functions

q=(qi)=(1,x1,x2,x12,x1x2,x22,…),(14)

are orthogonalized to the following orthogonal basis functions.

pi=qi−∑k=1i−1(qi,pk)(pk,pk)pk,(i=1,2,3,…).(15)

Eq. (7) can be expressed as

[(p1,p1)(p1,p2)⋯(p1,pm)(p2,p1)(p2,p2)⋯(p2,pm)⋮⋮⋱⋮(pm,p1)(pm,p2)⋯(pm,pm)][a1(x)a2(x)⋮am(x)]=[(p1,uI)(p2,uI)⋮(pm,uI)].(16)

If {pi(x)}(i=1,2,…,m) is a set of orthogonal functions, that is

(pi,pj)=0,i≠j,(17)

Eq. (16) becomes to

[(p1,p1)0⋯00(p2,p2)⋯0⋮⋮⋱⋮00⋯(pm,pm)][a1(x)a2(x)⋮am(x)]=[(p1,uI)(p2,uI)⋮(pm,uI)].(18)

We can simply have

ai(x)=(pi,uI)(pi,pi),(i=1,2,…,m).(19)

Then

a(x)=A∗(x)B(x)u,(20)

where

A∗(x)=[1(p1,p1)0⋯001(p2,p2)⋯0⋮⋮⋱⋮00⋯1(pm,pm)].(21)

Substituting Eq. (19) into Eq. (11), we obtain

uh(x)=∑I=1nΦI∗(x)uI=Φ∗(x)u,(22)

where Φ∗(x) mean the shape function, and

Φ∗(x)=(Φ1∗(x),Φ2∗(x),…,Φn∗(x))=PT(x)A∗(x)B(x).(23)

3  The FEFG Method for 3D Elasticity

The equilibrium equation of 3D elasticity is

σ11,1(x)+σ12,2(x)+σ13,3(x)+b1(x)=0,(24)

σ21,1(x)+σ22,2(x)+σ23,3(x)+b2(x)=0,(25)

σ31,1(x)+σ32,2(x)+σ33,3(x)+b3(x)=0,(26)

with the boundary conditions

ui=u¯i,(x∈Γu),(27)

σijnj=t¯i,(i,j=1,2,3x∈Γq),(28)

where σij mean the stress, bi mean the body force, ui mean the displacement, Γ=Γu∪Γq is the boundary of problem domain Ω, and Γu∩Γq=∅. u¯i mean the known displacement on Γu, t¯i mean the known traction on Γq, and nj is the unit external normal vector to Γq.

Then, the elasticity problem composed of Eqs. (24) and (25) and the boundary conditions is solved by the FEFG method.

The 3D problem domain Ω is divided into L+1 2D subdomains along direction x3, which is chosen as the splitting direction, that is, L−1 2D planes are inserted into Ω. And the distance between two adjacent 2D sub-domains is Δx3. That is

Ω=⋃k=0L−1Ω(k)×[x3(k),x3(k+1))∪Ω(L),(29)

where

a=x3(0)<x3(1)<⋯<x3(L)=c,x3∈[a,c],(30)

Δx3=x3(k+1)−x3(k)=(c−a)/(c−a)LL.(31)

For a fixed x3(k), Eqs. (24) and (25), the essential boundary conditions Eq. (27) and the natural boundary conditions Eq. (28) can be transformed as

σ11,1(k)+σ12,2(k)+b1(k)=−σ13,3(k),((x1,x2)∈Ω(k),x3=x3(k)),(32)

σ21,1(k)+σ22,2(k)+b2(k)=−σ23,3(k),((x1,x2)∈Ω(k),x3=x3(k)),(33)

with the boundary conditions

ui(k)(x1,x2)=u¯i(x1,x2,x3(k))=u¯i(k)(x1,x2),((x1,x2)∈Γu(k),i=1,2),(34)

σij(k)(x1,x2)nj(k)=t¯i(x1,x2,x3(k))=t¯i(k)(x1,x2),((x1,x2)∈Γq(k),i,j=1,2),(35)

where ui(k)(x1,x2) mean the displacement in the 2D domain Ω(k), u¯i(k)(x1,x2) mean the known displacement on Γu(k), t¯i(k)(x1,x2) mean the known traction on Γq(k), and Γ(k)=Γu(k)∪Γq(k), Γu(k)∩Γq(k)=∅, nj(k) mean the unit outward normal to boundary Γq(k).

Applying the IEFG method to solve the 2D problem of Eqs. (32)–(35), the trial function uh(x(k),x3(k)) at x(k)=(x1,x2) is related to the value at nodes xI(k) whose domain of influence covers x(k). It can be expressed as

uh(x(k),x3(k))=[u1(k)u2(k)]=[∑I=1nΦI∗(x(k))u1(xI(k),x3(k))∑I=1nΦI∗(x(k))u2(xI(k),x3(k))]=NU(k),(36)

where

N=[Φ1∗0Φ2∗0⋯Φn∗00Φ1∗0Φ2∗⋯0Φn∗],(37)

U(k)=(u1(x1(k)),u2(x1(k)),u1(x2(k)),u2(x2(k)),…,u1(xn(k)),u2(xn(k)))T.(38)

The strain at any point (x(k),x3(k)) in subdomain Ω(k) is

ε(k)=ε(x(k),x3(k))=[ε11(k)ε22(k)ε12(k)]=[(∑I=1nΦI∗(x(k))u1(xI(k),x3(k))),1(∑I=1nΦI∗(x(k))u2(xI(k),x3(k))),2(∑I=1nΦI∗(x(k))u1(xI(k),x3(k))),2+(∑I=1nΦI∗(x(k))u2(xI(k),x3(k))),1].(39)

It can be written as

ε(x(k),x3(k))=BU(k),(40)

where

B=[Φ1,1∗0Φ2,1∗0⋯Φn,1∗00Φ1,2∗0Φ2,2∗⋯0Φn,2∗Φ1,2∗Φ1,1∗Φ2,2∗Φ2,1∗⋯Φn,2∗Φn,1∗].(41)

The stress at any point (x(k),x3(k)) in subdomain Ω(k) is

σ(k)=σ(x(k),x3(k))=[σ1,1(k)σ2,2(k)σ1,2(k)]=DBU(k)=Dε(k).(42)

For plane stress problems,

D=E1−ν2[1ν0ν10001−ν2];(43)

and for plane strain problems,

D=E(1+ν)(1−2v)[1−vν0ν1−v0001−2ν2],(44)

where E means Young's modulus, ν means Poisson's ratio, and D means the elasticity matrix [54].

The penalty method is applied to apply the essential boundary conditions, and the weak form of the Galerkin integral of Eqs. (32)–(35) is

∫Ω(k)δ ε(k)T⋅σ(k)dΩ(k)−G∫Ω(k)δ u(k)T⋅u″(k)dΩ(k)−∫Ω(k)δ u(k)T⋅b(k)dΩ(k)−∫Γq(k)δ u(k)T⋅t¯(k)dΓ(k)+α∫Γu(k)δ u(k)T⋅S⋅(u(k)−u¯(k))dΓ(k)=0                     ,(45)

where

G=E2(1+ν),(46)

is the shear modulus, and

u′′(k)=(∂2u1(k)∂x32,∂2u2(k)∂x32)T=Φ∗(x(k))U′′(k),(47)

U′′(k)=(∂2u1(x1(k))∂x32,∂2u2(x1(k))∂x32,∂2u1(x2(k))∂x32,∂2u2(x2(k))∂x32,…,∂2u1(xn(k))∂x32,∂2u2(xn(k))∂x32)T,(48)

b(k)=(b1(k),b2(k))T,(49)

t¯(k)=(t¯1(k),t¯2(k))T,(50)

u¯(k)=(u¯1(k),u¯2(k))T,(51)

S=[s100s2].(52)

If ui(k) is the displacement on Γu(k) along direction xi, si=1, otherwise si=0.

Eqs. (36), (40), (42) and (47) are substituted into Eq. (45). Then, we have

∫Ω(k)δ[BU(k)]T⋅[DBU(k)]dΩ(k)−G∫Ω(k)δ[Φ*(x(k))U(k)]T⋅[Φ*(x(k))U″(k)]dΩ(k) −∫Ω(k)δ[Φ*(x(k))U(k)]T⋅b(k)dΩ(k)−∫Γq(k)δ[Φ*(x(k))U(k)]T⋅t¯(k)dΓ(k) +α∫Γu(k)δ[Φ*(x(k))U(k)]T⋅S⋅[Φ*(x(k))U(k)]dΓ(k) −α∫Γu(k)δ[Φ*(x(k))U(k)]T⋅S⋅u¯(k)dΓ(k)=0(53)

To get the solvable algebraic equations, all integrals in Eq. (53) need to be analyzed separately.

The first integral term in Eq. (53) is

∫Ω(k)δ[BU(k)]T⋅[DBU(k)]dΩ(k)=δU(k)T⋅[∫Ω(k)BTDBdΩ(k)]⋅U(k)=δU(k)T⋅K⋅U(k),(54)

where

K=∫Ω(k)BTDBdΩ(k).(55)

The second integral term in Eq. (53) is

G∫Ω(k)δ[Φ*(x(k))U(k)]T⋅[Φ*(x(k))U″(k)]dΩ(k)=δU(k)T⋅G[∫Ω(k)Φ*T(x(k))Φ*(x(k))dΩ(k)]⋅U″(k)=δU(k)T⋅C⋅U″(k),(56)

where

C=G∫Ω(k)Φ∗T(x(k))Φ∗(x(k))dΩ(k).(57)

The third integral term in Eq. (53) is

∫Ω(k)δ[Φ∗(x(k))U(k)]T⋅b(k)dΩ(k)=δU(k)T⋅∫Ω(k)Φ∗T(x(k))b(k)dΩ(k)=δU(k)T⋅F1,(58)

where

F1=∫Ω(k)Φ∗T(x(k))b(k)dΩ(k).(59)

The fourth integral term in Eq. (53) is

∫Γq(k)δ[Φ∗(x(k))U(k)]T⋅t¯(k)dΓ(k)=δU(k)T⋅∫Γq(k)Φ∗T(x(k))t¯(k)dΓ(k)=δU(k)T⋅F2,(60)

where

F2=∫Γq(k)Φ∗T(x(k))t¯(k)dΓ(k).(61)

The fifth integral term in Eq. (53) is

α∫Γu(k)δ[Φ∗(x(k))U(k)]T⋅S⋅[Φ∗(x(k))U(k)]dΓ(k)

=δU(k)T⋅[α∫Γu(k)Φ∗T(x(k))SΦ∗(x(k))dΓ(k)]⋅U(k)=δU(k)T⋅Kα⋅U(k),(62)

where

Kα=α∫Γu(k)Φ∗T(x(k))SΦ∗(x(k))dΓ(k).(63)

The sixth integral term in Eq. (53) is

α∫Γu(k)δ[Φ*(x(k))U(k)]T⋅S⋅u¯(k)dΓ(k)=δU(k)T⋅α∫Γu(k)Φ*T(x(k))Su¯(k)dΓ(k)=δU(k)T⋅Fα,(64)

where

Fα=α∫Γu(k)Φ∗T(x(k))Su¯(k)dΓ(k).(65)

Substituting Eqs. (54), (56), (58) (60), (62) and (64) into Eq. (53), we obtain

δU(k)T⋅(KU(k)−CU′′(k)−F1−F2+KαU(k)−Fα)=0.(66)

For the δU(k)T is arbitrary, we let

CU′′(k)+K¯U(k)=F,(67)

where

K¯=−K−Kα,(68)

F=−F1−F2−Fα.(69)

In order to obtain the solution of Eq. (67), points x3(1),x3(2),…,x3(L−1) are uniformly inserted along direction x3 in domain [a,c]. Let

U(x3(0))=U(0)=U(0)(a),(70)

U(x3(1))=U(1),(71)

U(x3(2))=U(2),(72)

⋮

U(x3(L−1))=U(L−1),(73)

U(x3(L))=U(L)=U(L)(c).(74)

Using the FDM in the splitting direction x3, we have

U′′(k)≈U(k−1)−2U(k)+U(k+1)(Δx3)2,(k=1,2,…,L−1).(75)

Then, Eq. (67) is written as

C⋅U(0)−2U(1)+U(2)(Δx3)2+K¯U(1)=F(1),(76)

C⋅U(1)−2U(2)+U(3)(Δx3)2+K¯U(2)=F(2),(77)

C⋅U(2)−2U(3)+U(4)(Δx3)2+K¯U(3)=F(3),(78)

⋮

C⋅U(L−2)−2U(L−1)+U(L)(Δx3)2+K¯U(L−1)=F(L−1).(79)

The matrix form of Eqs. (76)–(79) is

1(Δx3)2[HCCHCCHC⋱⋱⋱CHCCH][U(1)U(2)U(3)⋮U(L−2)U(L−1)]=[F(1)−CU(0)(Δx3)2F(2)F(3)⋮F(L−2)F(L−1)−CU(L)(Δx3)2],(80)

where

H=−2C+(Δx3)2K¯.(81)

Let

M=1(Δx3)2[HCCHCCHC⋱⋱⋱CHCCH],(82)