无网格法计算中一种新的积分点搜索算法

岩石力学与工程学报 Vol.28 Supp.2

2009年9月 Chinese Journal of Rock Mechanics and Engineering Sept.，2009

无网格法计算中一种新的积分点搜索算法

刘红生1，邢忠文1，杨玉英2

(1. 哈尔滨工业大学 机电工程学院，黑龙江 哈尔滨 150001，2. 哈尔滨工业大学 材料科学与工程学院，黑龙江 哈尔滨 150001)

摘要：无网格法计算中常采用背景网格积分法进行数值积分，需判断积分点与分析域的关系，当工程问题域边界较复杂时，虽可采用常用的射线法判断积分点与分析域的关系，但当分析域边界线段处于同一直线时，算法的复杂性大大增加，使其效率不高。为此，提出采用矩阵法判断积分点与分析域的关系，在此基础上，为加快搜索速度，提出并建立一种新的数据结构——稀疏矩阵，该矩阵中的元素记录各节点的信息，包括节点与其相邻节点的距离、节点编号和节点影响域尺寸等，通过稀疏矩阵可快速计算矩阵的行列式，从而给出一种积分点搜索速度更快的算法。算例表明，该算法大大加快积分点的搜索速度，提高无网格法的计算效率。 关键词：岩土工程；无网格法；背景网格；射线法；矩阵法；稀疏矩阵

中图分类号：TU 443 文献标识码：A 文章编号：1000–6915(2009)增2–3659–07

A NEW QUADRATURE POINT DETECTION ALGORITHM IN MESHLESS

CALCULATION

LIU Hongsheng1，XING Zhongwen1，YANG Yuying2

(1. School of Mechatronics Engineering，Harbin Institute of Technology，Harbin，Heilongjiang 150001，China；2. School of

Materials Science and Engineering，Harbin Institute of Technology，Harbin，Heilongjiang 150001，China)

Abstract：Background cell integration method is widely used to implement numerical integration in meshless computation，in this procedure，the relationship between quadrature points and the problem domain must be assured. Radial method is commonly used when the boundary is complicated. However，the complexity of the algorithm is increased while several line sections are parallel，which results in low efficiency. In order to fasten the quadrature point detection，the matrix method is proposed to implement the detection through calculating the determinant of matrix，and furthermore，a new data structure—sparse matrix is presented；each element in this sparse matrix stores the following information associated with the nodes such as the distance between the node and its neighbor，the number of the node and radius of the node influence domain and the computational cost of determination of matrix can be decreased dramatically by sparse matrix. As a result，an advanced algorithm，which can fasten the quadrature point detection more quickly，is generated. Numerical results show that this algorithm can fasten the quadrature point detection sharply and therefore improve the efficiency of meshless numerical computation. Key words：geomechanical engineering；meshless method；background cells；radial method；matrix method；sparse matrix

1 引 言

有限元法在工程领域得到了广泛地应用，岩土

收稿日期：2008–10–19；修回日期：2008–12–30

工程中常采用有限元法对岩土及地下工程进行数值分析。目前，开发了一些专门用于地下工程的有限元软件，如同济大学地下工程系开发的曙光软件，日本的2D-σ 和3D-σ。但由于岩土及地下工程问题

基金项目：哈尔滨工业大学优秀青年老师培养计划资助项目(HITQNJS–2009–014)

作者简介：刘红生(1976–)，男，博士，2001年毕业于哈尔滨工业大学材料学院材料成型及控制专业，现任讲师，主要从事计算力学方面的教学与研究工作。Email：hs_liu_hit@163.com

• 3660 • 岩石力学与工程学报 2009年

的复杂性和特殊性不同于其他结构，不仅涉及复杂材料、锚杆支护、梁支撑以及衬砌等的计算，而且还需分步模拟施工建造的全过程，要完整地建立一个岩土及地下工程的有限元分析模型，需大量的网格划分工作，并且要求技术人员有较高的力学素质和丰富的网格划分经验。

无网格法的出现为岩土及地下工程的数值计算提供了一种较为理想的数值方法。该方法是一种新的数值方法，相对传统有限元法，无网格法不具备网格，分析工程问题不必生成网格，且其形函数连续性高，理论上具有更高的模拟精度。无网格法在岩土力学中的应用得到了初步的研究[1

～6]

的数据结构——稀疏矩阵加快积分点的搜索，通过稀疏矩阵提高无网格法的计算效率。

2 无网格法的基本原理

仅用节点离散分析域Ω，如图1所示。每个节

点k具有一个影响域Ωk和权函数wk(x)，二维情况下，影响域为圆形，如Ωk，或者为矩形，如Ωj。若x∈Ωk，则wk(x)≥0，否则wk(x)= 0。若采用

Shepard函数建立无网格形函数Φ(x)，则节点k的形函数Φk可表示为

，但不具备网

格的特性给无网格法中的数值积分带来了困难。众多学者对无网格法中的数值积分做了大量工作。T. Belytschko等[7]提出了背景网格积分法，背景网格有别于有限元法中的单元，背景网格只是用来进行数值积分的一种工具，不具近似功能，与有限元法中的单元有本质的区别。为减少计算量和摆脱网格的约束，提出了节点积分法[8]和无网格法的局部Petrov- Galerkin(MLPG)法[9]。节点积分法不必借助网格，但节点积分易使数值解不稳定[10]，而MLPG法会破坏数值模型的对称性

[11]

Φk=wk(x)/∑wi(x) (1)

i=1

NP

式中：NP为分析域内的节点总数。

ΩΩk节点kΩj 节点j 。

图1 无网格法计算模型

Fig.1 The meshless computational model

背景网格积分法是无网格法中最常用的积分法，该方法具有数据结构简单和对称的数值模型等优点。背景网格有两种：Otree型和单元型。单元型背景网格类似于有限单元，即把分析域分割成为网格，进行数值积分时不需判断积分点与分析域的关系，但单元型背景网格的生成很费时。Otree型网格生成简单，只需用一张带有网格的方形块(二维情况下)完全覆盖分析域即可，但采用Otree型背景网格进行数值积分需判断积分点与分析域的关系。因此，如何快速判断积分点与分析域的关系是Otree型背景网格积分法中的关键，直接关系到无网格法的计算效率。判断点与多边形域关系的方法有多种，如叉积判断法、角度法、射线法[12]等，射线法是最常用的方法，但射线法的主要问题是多边形的多条边共线时很难判断射线与多边形的交点个数。为此，提出了改进的射线法[13]，但该方法过于复杂，当分析域形状复杂时，需用一多边形近似分析域的边界，当多边形的边数过大，射线法的效率不高。针对无网格法计算中积分点与分析域关系的判断，为提高无网格法的计算效率，本文提出采用矩阵法判断点与分析域的关系，在此基础上，提出并建立一种新

[7]

因此，x处的场量近似可以表示为

NP

u=∑Φiui* (2)

h

i=1

式中：ui*为节点i的节点虚值。

权函数采用三次样条函数[14]，可表示为 ⎧21⎞⎛23

4s4s s−+≤⎪⎜⎟ 32⎝⎠⎪⎪w(s)=⎨443⎛1⎞ (3) 2

⎪−4s+4s−s⎜2＜s≤1⎟33⎝⎠⎪⎪⎩0 (s＞1)

其中，

s=di/rinf，di=||x−xi||

式中：rinf为节点影响域半径(影响域为圆形)。

3 积分点搜索

3.1 积分点与分析域关系判断

进行点与多边形域的关系判断前，需定义分析

第28卷 增2 刘红生，等. 无网格法计算中一种新的积分点搜索算法 • 3661 •

域。由于无网格法不具备网格，S. F. Li等[15]在用无网格法分析接触碰撞问题时采用矩阵法定义分析域。该矩阵定义为

NP

R0便可确定P0∈Ω。

对于非均匀分布的节点，如图4所示，则矩阵行列式的分布如图5所示。矩阵法具有如下规律：对任意一节点Pi，若Pi∈Ω则det(M(xi)＞R，R= min[det(M(xi)](i = 1，2，\"， NP)。

M=∑PwiP∆Vi (4)

T

i=1

其中，

P={p1(x)，p2(x)，pm(x)}

0NP

式中：m为多项式基函数P(x)的项数；∆Vi=∆xi∆yi，∆xi，∆yi分别为节点i在x，y方向上的节点间距。矩阵M的元素Mkl为

Mkl=∑sk+l−2wi∆Vi (5)

i=1NP

图4 节点不均匀分布(一维)

Fig.4 Non-uniform distribution of nodes (one dimension)

积分点搜索时，可采用矩阵法判断积分点与分析域的关系。若分析域内的节点均匀分布(见图2)，若P0∈Ω则则矩阵法具有如下规律：任意一点P0，det[M(P0)] = R，否则det[M(P0)]→0，R为常数(见

det(M)

节点位置

图3)。图2为采用无网格节点离散的区域Ω，图3为在区域内及区域外det(M)的分布。

(0，b)

(a，b)

图5 矩阵M行列式的分布(一维) Fig.5 Distribution of determination of matrix M(one

dimension)

Ω

3.2 算 法

为获得矩阵M的行列式，需进行大量的计算，导致计算时间过长，严重影响无网格法的计算效率。

(0，a)

(0，0)

图2 区域的无网格离散

为此，本文提出并建立一种新的数据结构——稀疏矩阵减少矩阵行列式的计算量，提高无网格法的计算效率。

3.2.1 稀疏矩阵的建立

Fig.2 A meshless discretization of problem domain

分析域外det(M) = 0

分析域内det(M) = R

分析域内 0＜det(M)＜R

为提高积分点的搜索速度，基于分析域内分布的节点(见图6)建立稀疏矩阵。把节点坐标值xi按小到大的顺序存储在矩阵中的某一行，多个节点的xi值相同时，均存放在该行的同一列，如果节点坐标xi有k个不同值，则矩阵列数为k列。把xi从小

到大的顺序存储成行矩阵x(k)。同样地，对节点坐标yi也进行相同的存储，节点坐标yi有m个不同的

值，则矩阵的行数为m，存放yi值的矩阵为y(m)。最后得

1

5

2

4

6

图3 在区域内及区域外det(M)的分布

Fig.3 Distribution of the det(M) inside and outside the

problem domain

分析域边界处有0＜det(M)＜R，矩阵M行列式在边界处出现了一个缓冲区，见图3。因此，需设定

3

图6 无网格节点分布 Fig.6 Distribution of meshless nodes

一门槛值R0 = (0.7～0.8)R。任意点P0若det[M(P0)]≥

• 3662 • 岩石力学与工程学报 2009年

到k×m的稀疏矩阵N 1)⎡N(1，

⎢N(2， 1)N=⎢

⎢#⎢

⎣N(m，1)

N(1， 2)N(2， 2)#N(m，2)

\"\"%\"

N(1， k)⎤N(2， k)⎥⎥ (6)

⎥#

⎥

N(m，k)⎦

p1

p0

矩阵N便对每个节点建立了一空间向量N(i，j)，每个向量记录了节点的多项信息：节点编号、节点坐标和节点与相邻节点的距离等。需要指出的是，并不是每个向量N(i，j)代表一个节点编号，有些向量为空，即N(i，j) = 0。为存储稀疏矩阵虽需较大的内存空间，但通过稀疏矩阵N可快速建立矩阵和快速计算行列式，使得积分的点搜索大大加快。建立稀疏矩阵的思想是用空间换取时间。用图6所示的节点建立的稀疏矩阵为

图7 点p0和p1的覆盖域

Fig.7 The domains covered by influence area of p0 and p1

根据行列式的值判断该点是否在分析域内。以图7所分布的节点为例，nmax= 16。

具体算法描述如下： (1) 在分析域内布置节点；

(2) 在x和y方向分别建立节点坐标序列x(k)和y(m)；

(3) 设定节点影响域半径rinf； (4) 用分析域内的节点建立矩阵M； (5) 在分析域内的建立稀疏矩阵N； (6) 计算R0的值； (7) 设定nmax；

2，\"，N) (N(8) 循环要判断的积分点pi(i = 1，为要判断积分点数，点pi的坐标(xp，yp))：

① 如果xp＞x(N(1，k))或xp＜x(N(1，1))或yp＞y(N(1，m))或yp＜y(N(1，1))转至

)；

② 以点p为中心建立正方形Ωp，正方形Ωp的左边、右边、上边、下边分别定义为x1=lc，x2=rc，y1=uc，y2=dc(dc为正方形下边的纵坐标值)；

③ 找出lc≥x(N(1，j))≥rc的元素个数nr，设j)与分析域边界t+nr－1≥j≥t(t为横向上节点元素(i，间的节点数)；

④ 找出uc≥y(N(i，1))≥dc的元素个数n1；设q+n1－1≥i≥q(q为纵向上节点元素(i，j)与分析域边界间的节点数)；

⑤ 找出矩阵N(i，j)(t+nr−1≥j≥t，q+n1−1≥ i≥q)中非空元素，设总数为nz；

⑥ 如果nz≥nmax，转至

；

⑦ 如果nz≤nmax，记录在正方形Ωp内的节点号，用这些节点坐标建立矩阵M；

⑧ 计算点pi的行列式det[M(xp，xp)]； ⑨ 如果det[M(xp，xp)]≥R0，转至⑩ 如果det[M(xp，xp)]＜R0，转至 点pi在分析域内；

； ；

1)⎡K(1，

⎢K(2， 1)⎢

K=⎢K(3， 1)

⎢

1)⎢K(4，⎢ 1)⎣K(5，

K(1， 2)

K(2， 2)K(3， 2)K(4， 2)K(5， 2)K(1， 3)

K(2， 3)K(3， 3)K(4， 3)K(5， 3)K(1， 4)

K(2， 4)K(3， 4)K(4， 4)K(5， 4)

K(1， 5)⎤

K(2， 5)⎥⎥K(3， 5)⎥

⎥

K(4， 5)⎥K(5， 5)⎥⎦ (7)

式中：K(4，1) = 1，K(3，2) = 2，K(1，2) = 3，K(2，3) = 4，K(5，4) = 5，K(2，5) = 6，其他元素均为0。 3.2.2 算法描述

各积分点的判断均需计算式4中矩阵M的行列式，但由于计算M的行列式很费时，从图7可知，并不是所有点的搜索均需计算M的行列式。以图7为例，对于分析域内部的点p0，以点p0为中心，以2rinf(rinf为节点影响域的尺寸表征，对于非均匀分布的节点，rinf可根据自适应法[16]确定)为边长的正方形影响域所覆盖的节点数n均满足n≥16。因此，对于某些点的搜索并不需计算M的行列式便可确定该点与分析域的关系。分析域内的节点分布须满足可容性，即须保证所分布的节点能构造形函数。二维情况下，为保证节点分布具可容性，分析域内任意点至少需被3个节点的影响所覆盖[15]，反之，分析域内任意点的影响域至少需覆盖3个节点。为保证点在分析域内的充分性，二维情况下，该点影响域所覆盖的节点数为nmax，须满足nmax>>3，可取nmax=5～8×3。图7中的点p0影响域所覆盖的节点数L，若L≥nmax，不需计算det[M(p0)]便可确定

p0在分析域内，反之，则需计算该点处的行列式，

第28卷 增2 刘红生，等. 无网格法计算中一种新的积分点搜索算法 • 3663 •

点pi不在分析域内； (9) 计算结果。

背景网络

4 算 例

4.1 点与复杂边界区域关系的判断

σ

为验证矩阵法的适用性及效率，对点(积分点)与复杂边界的分析区域(见图8)的关系进行判断，即判断网格内的点与分析域的关系。分析域边界近似为17条线段构成的多边形，背景网格总数为25×29个，每个背景网格内布置4个点(积分点)。分别用射线法、矩阵法和矩阵法+稀疏矩阵的方法判断点(积分点)与分析域的关系。所用时间见表1。

条线段构成的多边形。分别采用射线法、矩阵法和矩阵法+稀疏矩阵的方法判断网格内的积分点与分析域的关系。其无网格法为无网格迦辽金法[7]。比较无网格法得到的x = 0截面上的x方向的正应力与解析解，比较结果如图10所示，从图10可看出两种方法的模拟结果均接近解析解。射线法，矩阵法和矩阵法+稀疏矩阵所用计算时间见表2。

Q

图9 圆孔板内节点布置

Fig.9 The arrangement of nodes in plate with a hole

图8 复杂形状区域

Fig.8 Domain with complex boundary

表1 射线法、矩阵法和矩阵法+稀疏矩阵判断点与分析域

的关系所用时间比较

Table 1 Comparison of time consumed by radial method，

moment method and moment method and sparse matrix for judging the relationship of point and analysis domain

判断方法 射线法 矩阵法 矩阵法+稀疏矩阵

所耗费时间/s

时间节省率/%

精确解

无网格迦辽金法(矩阵法) 无网格迦辽金法(射线法)

σ /Pa

y/m

图10 x = 0处截面的x向应力σx的比较 Fig.10 Comparison of the stress σx at x = 0

疏矩阵所用时间

and moment method and sparse matrix for tension circular plate with holes

判断方法

所耗费时间/s

时间节省率/%

9.32 0.0 表2 受拉带孔圆孔板计算时射线法、矩阵法和矩阵法+稀 6.26 32.8

2.59 72.1 Table 2 Time consumed by radial method，moment method

4.2 受拉的带孔圆孔板

设一个在左侧x方向承受均匀拉伸的平板(见图9)，其拉应力为σ = 1 Pa，板中半径为a = 1 m的小圆弧，材料的弹性模量E = 30 MPa；泊松比ν = 0.3。在平面内布置20×20个背景积分网格，每个网格内布置2×2个高斯积分点。圆孔板平面内布置777个节点。对于边界离散，圆弧用32段线段近似，其他边界用4段线段近似，分析域便可近似成由36

射线法 523.8 0.0 矩阵法 394.4 24.7 矩阵法+稀疏矩阵 270.8 48.3

4.3 隧道开挖算例

开挖的模拟是岩土工程数值方法研究的重要问

• 3664 • 岩石力学与工程学报 2009年

题之一。典型的隧道开挖模型如图11所示。岩体力学参数：重度γ=26 kN/m3，变形模量E = 3.6 GPa，泊松比µ= 0.2。根据对称性，取其隧道开挖模型的一半进行计算。开挖过程分为3步，模型的下边界受双向固定约束，左右边界受水平方向固定约束。按平面应变进行分析。对该开挖模型进行无网格法数值计算，采用背景网格进行数值积分，背景网格和积分点如图11所示，每个网格内分布积分点(无网格节点未显示)。图12为各开挖步中积分域边界的变化，图13为各开挖步中的主应力分布。采用射

线法、矩阵法和矩阵法+稀疏矩阵的方法积分点进行搜索。所用计算时间见表3。

表3 隧道开挖计算时用射线法、矩阵法和矩阵法+稀疏矩

阵所用时间

Table 3 Time consumed by radial method，moment method and moment method andsparse matrix in tunnel excavation

判断方法

所耗费时间/s

时间节省率/%

射线法 306.4 0.0 矩阵法 227.3 25.8 矩阵法+稀疏矩阵 155.3 49.3

背景网格

5 讨 论

采用传统有限元法分析岩石工程问题时，要求技术人员有较高的力学素质和丰富的网格划分经验，且需耗费大量时间划分网格。无网格方法为这

一困难的解决提供了有效的处理方法，它们只需节点而不需单元或网格信息，前处理特别简单，可免去网格划分。但现有各种无网格方法均存在计算时间长、效率低，不同材料的交界区域或不连续面的处理十分复杂，难以准确地施加边界条件等缺点。本文通过提出加快积分点搜索的算法，从而可在一定程度上提高无网格法的计算效率，还可从其他方面，如建立简化的形函数或材料形函数，提高无网格法的计算效率，这些问题是以后的研究方

积分点

图11 开挖隧道模型的背景网格和积分点

Fig.11 Background cells and integration points of calculation

model of tunnel

向。

(a) 第1步 (b) 第2步 (c) 第3步

图12 各开挖步中积分域边界的变化

Fig.12 Variation of boundary of integration domain in steps of

excavation

6 结 论

无网格计算中采用背景网格积分法进行数值积

1.0 MPa分时，提出了用矩阵法判断积分点与分析域的关系，在此基础上，为提高积分点的搜索速度，提出并建立了稀疏矩阵，给出了一种更快捷的积分点搜索算法。数值算例表明，相对于射线法，矩阵法可提高积分点的搜索效率。稀疏矩阵中各元素可记录各节点的多种信息，通过稀疏矩阵可快速确定积分点影

(a) 开挖前 (b) 第1步 (c) 第2步 (d) 第3步

响域内所覆盖的节点及其影响域，从而可快速建立矩阵，减少行列式的计算量，改善了矩阵行列式的计算效率，进一步地提高了积分点的搜索效率，从而提高了无网格法的计算效率。算法具有易于实现和搜索速度快的特点。对复杂形状的分析域，算法

图13 各开挖步中的主应力分布(σ = 2.0 MPa) Fig.13 Distribution of principal stress in steps of excavation

(σ = 2.0 MPa)

第28卷 增2 刘红生，等. 无网格法计算中一种新的积分点搜索算法 • 3665 •

的优势尤为明显，算法不需利用边界的离散线段信息，只需分析域内的节点信息。

参考文献(References)：

[1] COETZEE C J，BASSON A H，VERMEER P A. Discrete and

continuum modelling of excavator bucket filling[J]. Journal of Terramechanics，2007，44(2)：177–186. [2]

ANDREW J D，CHARLES E A. A hybrid meshless local Petrov- Galerkin method for unbounded domains[J]. Computer Methods in Applied Mechanics and Engineering，2007，196(4)：843–852. [3] WANG W D，WANG J G，WANG Z L，et al. An unequal-order radial

interpolation meshless method for Biot's consolidation theory[J]. Computers and Geotechnics，2007，33(1)：61–70.

[4] WANG Z L，WANG J G，SHENC R F. The application of a meshless

method to consolidation analysis of saturated soils with anisotropic damage[J]. Computers and Geosciences，2008，34(5)：849–859. [5]

卢 波，丁秀丽，邬爱清. 无网格法对岩体不连续面的模拟[J]. 岩土力学与工程学报，2008，27：2 108–2 117.(LU Bo，DING Xiuli，WU Aiqing. Modeling of rock discontinuity with meshless method[J]. Chinese Journal of Rock Mechanics and Engineering，2008，27(10)：2 108–2 117.(in Chinese)) [6]

李树沈，周锦添，李术才. 新型边界积分方程的无单元解及其应用[J]. 岩土力学，2007，28(12)：2 549–2 559.(LI Shushen，ZHOU Jintian，LI Shucai. An element-free solution of the new boundary integral equation and its application[J]. Rock and Soil Mechanics，2007，28(12)：2 549–2 559.(in Chinese))

[7] BELYSTCHKO T，KRONGAUZ，ORGAN Y D，et al. Meshless

methods：an overview and recent developments[J]. Computer Methods in Applied Mechanics and Engineering，1996，185(1)：3–47. [8]

CHEN J S，WANG D D.A constrained reproducing kernel particle formulation for shear deformable shell in Cartesian coordinates[J]. International Journal for Numerical Methods in Engineering，2006，68(2)：151–172.

[9] 赵美玲，聂玉峰，袁占斌. 无网格局部Petrov-Galerkin方法中本质边界条件的处理[J]. 应用力学学报，2006，23(2)：493–495.(ZHAO Meiling，NIE Yufeng，YUAN Zhanbin. Imposing essential boundary conditions in meshless local Petrov-Galerkin (MLPG) method[J]. Chinese Journal of Applied Mechanics，2006，23(2)：493–495.(in Chinese))

[10] YOO J W，BRIAN M，CHEN J S. Stabilized conforming nodal

ntegration in the natural element method[J]. International Journal for Numerical Methods in Engineering，2004，60(5)：861–890. [11] LI Q，SORIC J，JARAK T，et al. A locking-free meshless local

petrov-galerkin formulation for thick and thin plates[J]. Journal of Computational Physics，2005，28(1)：116–133.

[12] 童劲松，蔡 青，张博锋. 一种快速的点–多边形/多面体关系分

类算法[J]. 西北工业大学学报，1996，14(1)：130–133.(TONG Jinsong，CAI Qing，ZHANG Bofeng. A fast method for classifying point to polygon polyhedron relationships[J]. Journal of Northwestern Polytechnical University，1996，14(1)：130–133.(in Chinese)) [13] 王泽根. 射线法判定点与多边形包含关系的改进[J]. 解放军测绘

学院学报，1999，16(2)：130–132.(WANG Zegeng. The improvement of method of the ray to judge the relation of polygon and point[J]. Journal of the PLA Institute of Surveying and Mapping，1999，16(2)：130–132.(in Chinese))

[14] ZHOU J X，WANG X M，ZHANG Z Q，et al. Explicit 3D RKPM

shape functions in terms of kernel function moments for accelerated computation[J]. Computer Methods in Applied Mechanics and Engineering，2005，194(8)：1 027–1 035.

[15] LI S F，DONG Q，LIU W K，et al. A meshfree contact-detection

algorithm[J]. Computer Methods in Applied Mechanics and Engineering，2001，190(10)：3 271–3 292.

[16] 何沛祥，孙昌玲. 一种自适应影响域半径无网格Galerkin法[J]. 合

肥工业大学学报，2004，27(1)：88–91.(HE Peixiang，SUN Changlin. Element free Galerkin method in consideration of self adaptive radius of influential domain[J]. Journal of Hefei University of Technology，2004，27(1)：88–91.(in Chinese))