弹性动力学的相似边界元法

2
重庆建筑大学学报
第 22 卷
对式（l）~（4）作用 LapIace 或 Fourier 变换，即得相应的变换域中的方程，再应用加权残数法， 得到弹性动力学问题在变换域中的边界积分方程
Ck（i 0 ）U（i 0 ）=
Ueki t i c! -
teki Ui c! +
U
e ki
f
i
c"
!
!
"
（5）
其中
可求得其它单元的相应矩阵，然后通过迭加建立代数方程组系数矩阵。与通常的每个单
元都各自进行积分计算相比，本文方法可大幅度减少计算量。
关 键 词：弹性动力学；边界积分方程；相似单元；相似边界元法
中图分类号：0343
文献标识码：A
在边界元法中，要建立最终可求解的代数方程组，需要在所有单元上进行大量的积分运算。当 单元数目较多时，单元上积分的计算量将大大增加计算时间。
eiement method
弹性动力学的相似边界元法
作者： 作者单位：
刊名： 英文刊名： 年，卷(期)： 被引用次数：
程玉民， 彭妙娟 程玉民(上海大学上海市应用数学和力学研究所, 上海 200072)， 彭妙娟(上海大学 土木工程系, 上海 200072)
CHENC Y U-minl，PENC Miao- Uan2
（l . Shanghai Institute of Appiied Mathematics and Mechanics，Shanghai University，Shanghai 200072，Chian；2 . Department of Civii Engineering，Shanghai University，Shanghai 200072，China）
Abstract：For the boundary eiement method of eiastodynamics，some properties of matrices are discussed in case of simiiar boundary eiements and the simiiar boundary eiement method is presented . In a series of simiiar boundary eiements，when the corresponding matrices of a boundary eiement are obtained，the ones of other boundary eiements in the series can be obtained by proportion . Then the coefficient matrix of the iast system of iinear aigebraic eguations can be obtained by the method of superposition . Compared with the generai boundary eiement method，the computing speed can be raised by the simiiar boundary eiement method given in this paper . Keywords：eiastodynamics；boundary integrai eguation；simiiar boundary eiements；simiiar boundary
2 相似边界元法
由方程（9）可知，只要求得矩阵〔 H 〕和〔 G 〕，求解方程组即可求解。而〔 H 〕和〔 G 〕是分别由
〔 HI〕和〔 GI〕组装而成的。要求得〔 HI〕和〔 GI〕，需要做大量的如下形式的积分
+l +l
-l -l teki N m（#，$）（#，$）c#c$ = 0
（l0）
+l +l
0
为边界点；U
e ki
和
t
eki 为基本解；Ck（i
0
）是自由项，与
0
点处边界的几何特征有关；U i
和 ti
分
别表示变换域中的位移和面力分量。
在无体力的情况下，边界积分方程（5）可离散为
NM
+l +l
Ck（i 0 ）U（i 0 ）=
!
!
t
m i
I
I=l m=l
-l
-
l
U
e ki
N
m（#，$）（#，$）c#c$
在有限元法的研究中，文［1］讨论了单元相似时，相似单元之间单元刚度矩阵的关系。在无限 元法的研究中，文［2］也讨论了单元相似时的情况。在边界元法中，若能建立相似单元之间相应矩 阵的关系，则可以减少大量的积分运算。
为此，本文对弹性动力学问题的边界元法讨论了边界单元相似时单元之间的一些矩阵关系，建 立了弹性动力学的相似边界元法。在一组相似单元中，只要求得一个单元的相应矩阵，通过比例关 系即可求得其它单元的相应矩阵，然后通过迭加建立线性代数方程组的系数矩阵。与通常的每个 单元都各自进行积分计算相比，本文方法可大幅度减少计算量。
下面来推导相似边界元法的公式。
第6期
程玉民等：弹性动力学的相似边界元法
3
将一局部区域剖分为一组相似单元，即有
{x
I i
=
Ox
I i
-1
y
I i
=
Oy
I i
-1
z
I i
=
Oz
I i
-（1 如果需要的话）
其中：x
Ii ，y
I i
和
z
I i
为此局部区域中第
I
个单元的第
i
个节点的坐标，O
为比例常数。那么
I = O2 I -1
应力和位移满足如下边界条件
} t（i x ，t）= Ti I = p（i x ，t） x "TO
U（i x ，t）= g（i x ，t）
x " TU
（2）
和初始条件
} U（i x ，0+ ）= Ui（0 x ）
U'（i x ，0+ ）= 1 i（0 x ） x " O
（3）
本构关系为
T（i x ，t）=
弹性动力学的相似边界元法!
程玉民1， 彭妙娟2
（1 . 上海大学 上海市应用数学和力学研究所，上海 200072；2 . 上海大学 土木工程系，上海 200072）
摘要：讨论了弹性动力学边界元法中边界单元相似时单元之间的一些矩阵关系，建立了相
似边界元法的公式。在一组相似单源自文库中，只要求得一个单元的相应矩阵，通过比例关系即
Jacobi 行列式
（#，$）=
3r 3#
X
3r 3$
（7）
对式（6）进行数值求解即得变换域中边界节点的位移和面力分量，然后由数值反变换求得时间
域中的解。
将边界积分方程（6）写成矩阵形式
N
N
〔 C〕〔 U 〕= !〔 GI〕〔 T I〕- !〔 HI〕〔 U I〕
I=l
I=l
（8）
其中〔 U 〕=｛Ui｝，i = l，2，…，K 为节点位移向量；〔 U I 〕=｛U Ii ｝，i = l，2，…，M 为单元 I 的节点位
P〔（
c
2 1
-
2 c22）Um，m（ x ，t）Um，m（ x ，t）Si
+ c（22 Ui，（ x ，t）+ U ，（i x ，t））〕
x "O
（4）
其中 TO 和TU 分别为已知面力和位移的边界，TO#TU =T；Si 为 Kroneker DeIta。
! 收稿日期：2000 - 05 - 01 基金项目：国家自然科学基金资助项目（59608005） 作者简介：程玉民（1965 - ），男，山西人，教授，博士导师，博士，主要从事计算力学、大跨空间结构研究。
-l -l Ueki N m（#，$）（#，$）c#c$
（ll）
在非奇异单元上采用 Gauss 积分；在奇异单元上需对奇异积分做消除奇异性或降低奇异性阶数处
理，或利用刚体位移特解来求奇异积分。
为减少形如（l0）和（ll）的积分的计算量，本文提出相似单元的概念，建立了相似边界元法。
相似边界元法的基本思想是，将弹性体所在的区域 " 的边界! 按照边界条件或几何形状划分 为若干局部区域，然后将每个局部区域剖分为若干相似单元，建立相似单元上矩阵〔 HI〕及〔 GI 〕的 关系。在一组相似单元中，只要求得某个单元上的〔 HI 〕和〔 GI 〕，其它单元上的〔 HI 〕和〔 GI 〕即可 由比例关系得到而不需再做积分运算。
NM
+l +l
-
!
!
U
m i
I
I=l m=l
-l
-l teki N m（#，$）（#，$）c#c$
（6）
其中
N
为边界单元数，每个单元均为有
M
个节点的等参单元，边界节点总数为
K
；U
m i
I
，tim
I
分别
表示第 I 个单元第 m 个节点的 i 方向的位移和面力；N（#，$）， = l，2，…，M 是形函数；（#，$）为
可得
（12） （13）
〔 HI〕= O〔2 HI -1〕
（14）
〔 GI〕= O〔2 GI -1〕
（15）
这样就建立了相似单元之间矩阵〔 HI〕及〔 GI〕的关系。则在一局部区域中，不论有多少单元，
只要通过积分求得〔 H 1〕及〔 G 1〕，由式（14）和式（15）可求得此局部区域中所有单元上的〔 HI 〕和
! K !（ I）=
2P k+1
2LT〔4U（B I）- V （C I）〕
其中对平面应变问题 k = 3 -
（19）
4U，对 平 面 应 力 问 题 k =（3 U）（/ 1 +U）；U（B I ）和UC（ I ）分别 为裂纹尖端单元上 B 点和 C 点
在 I 时刻 y 方向的位移。 计算 所 得 的 正 则 应 力 强 度
［2］ 应隆安 . 无限元方法〔 M〕. 北京：北京大学出版社，1992 ［3］ 嵇醒，臧跃龙，程玉民 . 边界元法进展及通用程序〔 M〕. 上海：同济大学出版社，1997
4
重庆建筑大学学报
第 22 卷
Similar Boundary Element Method in Elastodynamics
因子（即 K !（ I ）/ K 1，K 1 为 静 态应力强度因子）及其与前人所
做工 作 的 比 较 见 图 2。可 以 看
出，本文方法的计算结果与其它
方法基本吻合。
参考文献：
图 1 中心裂纹板
图 2 应力强度因子
［1］ Leung AYT and Su RKL . Mode I Crack Probiems by Fractai T wo Levei Finite Eiement Method〔s J〕. Engineering Fracture Mechanics，1994，4（8 6）：847 ~ 856
第 22 卷 第 6 期
重庆建筑大学学报
VoI. 22 No . 6
2000 年 12 月
JournaI of Chongging Jianzhu University
Dec. 2000
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
文章编号：1006 - 7329（2000）06 - 0001 - 04
〔 GI〕而不必再做大量的奇异或非奇异积分运算，从而大幅度减少计算量。
对二维问题，有
I = O I-1 〔 HI〕= O〔 HI -1〕 〔 GI〕= O〔 GI -1〕
（16） （17） （18）
3 算例
以下对突加载荷6 作用下的二维矩形中心裂纹板的应力强度因子作了计算。 矩形中心裂纹板的模型如图 1 所示，板长 4 cm，宽 2 cm，裂纹长度为 2 a，a = 0 . 24 cm。材料的 弹性模量 E = 2 . 0 X 105 M N / m2，泊松比U= 0 . 3，质量密度P = 0 . 005 M N·s2 / m4。 对!型对称裂纹，其动态应力强度因子的计算公式为
1 弹性动力学边界元法
对于均匀、各向同性的线弹性材料，在小变形条件下，运动微分方程为
（ c21
-
c22）Ui，（i x ，t）+
c
2 2
U
，（ii
x
，t ）+
f（ x ，t）=
U（ x ，t）
x "O
（1）
其中 c1 和 c2 分别为弹性体内膨胀波和畸变波的传播速度；f 为体力分量；O 为弹性体所在的域， 其边界为 T；x =（ x 1，x 2，x 3）（对三维问题）或 x =（ x 1，x 2）（对二维问题）。
移向量； 〔 T I〕=〔 t Ii 〕，i = l，2，…，M 为单元 I 的节点面力向量；〔 C〕为与〔 U 〕对应的自由项矩阵；
〔 GI〕、 〔 HI〕分别为与〔 T I〕和〔 U I〕对应的系数矩阵。
进一步，方程（8）可写为
（〔 C〕+〔 H 〕）〔 U 〕=〔 G 〕〔 T 〕
（9）
其中〔 T 〕=｛ti｝，i = l，2，…，K 为节点面力向量；〔 H 〕、〔 G 〕是分别由〔 HI 〕、〔 GI 〕组装而成的对应 于〔 U 〕和〔 T 〕的系数矩阵。