- 1 -
中国科技论文在线
无单元伽辽金法在水工结构裂纹仿真中的
应用#
程井1，孙晓勇2，张媛1*
基金项目：水资源与水电工程科学国家重点实验室开放基金(No. 2012B103)；教育部博士点新教师基金(No.
20100094120005)；国家自然科学基金(No. 51009056)
作者简介：程井（1982-），男，讲师，博士后，主要从事高坝结构数值仿真方面的研究
（1. 河海大学水利水电学院，南京 210098； 5
2. 中南电力设计院，湖北 武汉 430071）
摘要：采用节点来离散域的无单元法由于不需要网格，特别适合模拟裂纹扩展问题。本文介
绍了无单元伽辽金法、衍射法模拟不连续面、围线积分法计算裂纹应力强度因子等的基本原
理和方法，编制了裂纹及其扩展计算的 VC＋＋程序 CPTEFG。应用典型单边纹和斜裂纹模10
型对程序进行了验证。以典型重力坝为例介绍了裂纹扩展的模拟过程及结果，最后对无单元
法计算裂纹扩展的一些问题和难题进行了总结。
关键词：大坝裂纹；无单元伽辽金法；修正衍射法；围线积分；裂纹扩展；应力强度因子
中图分类号：TV315
15
Crack propagation simulation via EFGM and its application
in hydraulic structures
CHENG Jing
1
, Sun Xiaoyong
2
, ZHANG Yuan
1
(1. College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing
210098; 20
2. Central SouthernChinaElectricPower DesignInstitute, Wuhan 430071)
Abstract: Due to its no-remeshing property, the meshless method is very suitable for cracking
problems. In this paper the main theory of meshless method is presented as well as the diffraction
method simulating the discontinuity and the contour integral method for SIFs (Stress Intensity
Factors) calculation, and then a VC++ program is developed. Both straight side plane crack and 25
slope side plane crack are simulated to validate the program. Then taking a typical gravity as an
example, the crack propagation is simulated by the proposed approach. At last some key problems
for crack propagation simulation are presented.
Key words: dam crack; Element-free Galerkin Method; modified diffraction method; contour
integral method; crack propagation; stress intensity factor (SIF) 30
0 引言
无单元伽辽金法是在有限差分法、有限单元法、边界元法之后发展起来的一种新型的数
值模拟方法。由于无单元法不需要积分网格，它在流体动力学、断裂力学以及刚塑性力学等
领域有广泛的研究前景。目前，无网格法主要有光滑质点流体动力学方法 SPH、漫射元法、35
无单元伽辽金法 EFGM，重构核粒子法 RKPM，基于配点离散的有限点法 FPM，无网格配
点法 PCM，单位分解有限元法 PUFEM，hp 云图法，局部边界积分方程法 LBIE 和无网格局
部彼得柯夫－伽辽金法 MLPG，点插值法 PIM 等等。这十几种方法中已用于开裂计算的有
EFGM
[1-3]、PUFEM、MLPG[4]、RKPM[5]等。其中尤以美国西北大学 Belytschko 教授提出并
发展的 EFGM 法应用最为广泛。 40
水工混凝土结构裂纹的计算和扩展模拟是一大研究热点，特别是温度裂纹的扩展问题。
- 2 -
中国科技论文在线
本文作者在文献[6]中介绍了无单元法计算温度场的基本原理和方法，在文献[7]中推导了考虑
温度作用的坝体应力计算方法。本文重点介绍二维开裂问题的仿真计算及裂纹扩展的判断和
计算。主要思路为：先采用节点离散整个计算域，对于开裂面即不连续面采用衍射法进行模
拟。然后采用无单元伽辽金法基本原理求解偏微分方程组并通过计算得到整个计算域的位移45
场和应力场。对于型裂纹及复合型裂纹，均可以考虑采用围线积分法计算应力强度因子 K
及 K，并通过判断确定裂纹是否扩展。如果满足扩展条件，则选择合适的裂纹扩展步进行
扩展模拟计算。
1 无单元伽辽金法基本理论
最小二乘法 50
无单元伽辽金法是一种采用伽辽金法来离散控制偏微分方程的数值方法，同有限元法一
样也需要借助背景网格来实现离散域的积分。但无单元伽辽金法的近似方案不是依靠单元内
部固定节点的插值得到，而是通过另外一种称之为“移动最小二乘近似”的方法得到。由于
控制偏微分方程的离散与积分方案基本同有限元法，因此这里不做赘述。本节对移动最小二
乘近似的基本原理进行了简要介绍。 55
静力学问题的基本待求场函数为位移场 ( )u x 。假定求解域内节点 Ix 处的函数值
( )I Iu u x 已知。对于某计算点 0x 局部邻域，构造其近似场函数为 0( , )
hu x x ：
( ) ( , ) ( ) ( , ) ( ) ( )
m
h T
o j j o o
j
u u p a x x x x x x p x a x (1)
其中 ( )p x 为空间坐标的完备基函数，m为基函数的个数。
对于二维问题： 60
2 2
( ) [1, , ], 3
( ) [1, , , , , ], 6
T
T
x y m
x y x xy y m
p x
p x
一次基：
二次基：
(2)
设计算点 0x 的邻域内包含 N 个节点，近似函数通过在这些节点 Ix x 处的函数值误差
的最小加权平均值来求得：
2 2
1 1
( )[ ( , )] ( )[ ( ) ( )]
N N
h T
I o I I o I o I I o
I I
S u u u
x x x x p x a x (3)
式中， ( )I o x 为计算点相对于节点的权函数，它具有紧支性。本文采用如下的权函数： 65
2 3 41 6 8 3 1
( )
0 1
r r r r
r
r
(4)
其中
00
/ Ir d xx x ， 0dx 为计算点 0x 的影响域半径。
令 S取最小值：
0 0
1 10
2 ( )[ ( ) ( ) ] ( ) 0 1,2,...
( )
N m
I i I i I j I
I ij
S
p a u p j m
a
x x x x
x
(5)
由此得： 0 0 0( ) ( ) ( )A a Bx x x u (6) 70
式中 0 0
1
( ) ( ) ( ) ( )
N
T
I I I
I
A
x x p x p x ,
- 3 -
中国科技论文在线
0 1 0 1 2 0 2 0( ) [ ( ) ( ) ( ) ( ) ... ( ) ( )]N NB x x p x x p x x p x , u 为节点位移序列。
由式（6）可以确定待定系数
1
0 0 0( ) ( ) ( )a A B
x x x u (7)
将上式代入(1)式，即 75
1
0 0 0( , ) ( ) ( ) ( ) ( , )
h T
ou A B N
x x p x x x u x x u (8)
其中形函数为：
1
0 0 0( , ) ( ) ( ) ( )
TN A Bx x p x x x (9)
对于求解域中的所有点都可以在其邻域内建立待求函数 ( )u x 的局部最佳近似。这些局
部近似函数 ( , )h ou x x 在点 0x x 处的值的集合就构成了待求函数 ( )u x 在求解域内的全局80
近似函数 h ( )u x ：
0
0( ) ( ) ( , ) ( )
h hu u u
x x
x x x x N x u (10)
其中形函数为：
0
1
0( ) ( , ) ( ) ( ) ( )
T A B
x x
N x N x x p x x x (11)
不连续面的模拟 85
本文采用了文献[8]提出的修正衍射准则考虑不连续面。当计算点 x与节点 Ix 的连线经过
裂纹面时，二者的修正距离计算公式如下：
1 2
1 0
0
( ) ( )
( ) ( )
( )
ps x s x ks
r x s x
s x
(12)
式中， 1 1( ), ( )s x s x 分别表示计算点 x、节点 Ix 到 P3 的距离；
0 ( )s x 为计算点 x与节点 Ix 的空间真实距离； 90
ps 表示在 P2－P3 段之前的裂纹累计长度，在此例中， 3 4ps PP ；
, k 均表示计算常数，取值范围为 ～。本文取 = ； 。
图 1 衍射法模拟不连续面
95
图 2 围线积分示意图
- 4 -
中国科技论文在线
2 围线积分法计算裂纹应力强度因子
围线积
本文采用围线积分法（见图 2）计算应力强度因子：先用无单元法计算离散域的位移场
和应力场，然后对计算式（13）左端的围线积分，式中符号的具体含义参考文献[9]： 100
2
1 2( )i i i i I IIu t u t ds mC m C
(13)
依 Betti 功互等定理可得：
1 22 2
1 1
I II
Gm Gm
K K
， (14)
其中
3
,
1,
2 1
3 4 ,
E
G
平面应力， 为泊松比
=
平面应变
裂纹扩展的判据 105
对于复合型裂纹扩展判据，采用最大拉应力理论[10]。该理论假定裂纹沿最大周向拉应
力的方向起裂。其中，等效应力强度因子为：
0 0sin (3cos 1) 0I IIK K （15）
2 2
0
8
2arctan( )
4
I I II
II
K K K
K
（16）
由此，等效应力强度因子为： 110
20 0
0
3
cos( ) cos ( ) sin
2 2 2
Ieq I IIK K K
（17）
裂纹在最大周向拉应力
max
到达临界值 c 时开始失稳扩展，即
,maxIeq IcK K （18）
,maxIeqK 是与最大周应力 max 对应的等效应力强度因子， IcK 是材料的断裂韧度。
3 基本算例分析 115
依据无单元伽辽金法基本理论和相关断裂力学理论，编制了无单元法计算结构开裂的
VC++程序 CPTEFG。关于算例的验证，参考文献[8]。本文给出了重力坝典型坝踵裂纹的扩
展算例。
典型重力坝剖面，坝高 ，底宽 ，初始计算节点总数 1943，积分单元数 1832。
假定坝体与坝基材料特性相同，均为弹模 E＝20GPa，泊松比为 ，线胀系数 ×10－120
5/℃。坝体上游面距坝基面 处有长度为 的初始裂纹。计算工况：坝体上游水位为
36m，坝体整体温降 10℃。不考虑裂缝水压力。裂纹初始状态下位移及应力计算结果见图 3。
应力强度因子及等效应力强度因子、开裂角如下：
o
. , . , , I II IeqK MPam K MPam K ＝
假定该坝体材料断裂韧度为 ，设置扩展步长为 ，扩展 6 步时的计算结125
- 5 -
中国科技论文在线
果见图 4。
（A）顺流向位移（单位：mm）
（B）铅直向位移（单位：mm）130
（C）顺流向正应力（单位：MPa）
（D）铅直向正应力（单位：MPa） 135
图 3 典型重力坝位移及应力计算结果等值线图
图 4 典型重力坝裂纹扩展路径及扩展后顺流向应力等值线图（应力单位：MPa）
通过以上算例的计算分析发现：
（1）节点及积分单元的尺寸对于计算结果有较大的影响，特别是缝端积分单元及节点140
- 6 -
中国科技论文在线
间距，必须达到裂纹长度的 1/10~1/20 才能得到较好的精度。
（2）裂纹附近节点沿垂直于裂纹扩展方向的尺寸与沿扩展方向的尺寸比例不宜过大，
宜控制在 ～ 之间。
（3）积分围线的半径也要控制在一个合理的范围内。尽量将积分围线布置在较为准确
的应力场处，从而避开缝端奇异及大梯度应力。 145
4 结论
本文简要介绍了无单元法计算裂纹应力强度因子的基本思路和方法，并编制了相应的计
算程序。典型裂纹实例的计算结果表明该方法具有较好的精度，适于水工结构裂纹模拟及其
扩展的仿真计算。计算结果具有较高的精度。但对于局部不规则凹形物理域，计算结果受支
撑域半径影响较大，在仿真模拟时需要特别注意考虑。另外，本文采用线弹性裂纹模型模拟150
坝踵裂纹，实际混凝土裂纹属于弹塑性裂纹，如何将混凝土弹塑性断裂力学应用于坝踵裂纹
仿真及评价，仍有许多工作要做。
致谢
衷心感谢水资源与水电工程科学国家重点实验室开放研究基金 (No. 2012B103)，教育部
博士点新教师基金 (No. 20100094120005) 以及国家自然科学基金 (No. 51009056) 的资金155
的支持。
[参考文献] (References)
[1] Belytschko T, Lu Y Y, Gu L. Crack propagation by element free Galerkin methods[C]. New Orleans, LA, USA:
Publ by ASME, New York, NY, USA, 1993.
[2] Belytschko T, Lu Y Y, Tabbara M. Dynamic fracture by the element-free Galerkin method[C]. Chicago, IL, 160
USA: ASME, New York, NY, USA, 1994.
[3] Belytschko T, Lu Y Y, Gu L. Crack propagation by element-free Galerkin methods[J]. Engineering Fracture
Mechanics. 1995, 51(2): 295-315.
[4] Gao L, Liu K, Liu Y. Applications of MLPG method in dynamic fracture problems[J]. CMES - Computer
Modeling in Engineering and Sciences. 2006, 12(3): 181-195. 165
[5] Tang J, Wang H P. A probe of the meshfree methods in crack propagation life prediction[J]. American Society
of Mechanical Engineers, Pressure Vessels and Piping Division (Publication) PVP. 2000, 400: 99-106.
[6] 程井，常晓林，周伟，等. 应用无单元伽辽金法计算大坝稳定温度场[J]. 武汉大学学报（工学版），2008：
41(3).
[7] 常晓林，程井，周伟. 无单元法在水工混凝土结构温度应力计算中的应用研究[J]. 水利学报，2008，39(7). 170
[8] 程井，常晓林，周伟，卢建华. 基于无单元伽辽金法的水工结构温度应力及温度裂纹扩展计算[J]. 四川
大学学报（工程科学版），2009，41(4)：26-30，75.
[9] 杨晓翔，范家齐，匡震邦. 求解混合型裂纹应力强度因子的围线积分法[J]. 计算力学学报，1996(01).
[10] F E, S G C. On the crack extension in plates under plane loading and transverse shear[J]. J Basic Engng.
1963, 85: 519～527. 175
