含节理单元的三维p型自适应有限元解法_费文平

3003 年 4 月
水利
SHUILI
学报
XUEBAO
文章编号: 0559- 9350 ( 2003) 04-0058- 09
第4期
含节理单元的三维 p 型自适应有限元解法
费文平1, 陈胜宏1
( 1 武汉大学 水利水电学院, 湖北 武汉 430072)
摘要: 本文提出了含三维无厚度节理单元、等厚度 节理单 元和变 厚度节 理单元 的 p 型自 适应有 限元模 型, 给出
外在有关 p 型自适应有限元分析的研究中, 尚未涉及到这类特殊单元的处理问题, 从而使研究成果的 工程应用受到一定程度的限制。对于常规块体单元的三维 p 型有限元模型, 作者在文献[ 9] 中已有详
尽的论述, 本文主要给出三维无厚度节理单元、等厚度节理单元和变厚度节理单元的 p 型有限元模
型, 并给了升阶谱的计算格式。实例分析结果表明, 用 p 型有限元法来求解含三维节理单元的有限元
研究小组从 1995 年开始研究水工结构的 h 型弹粘塑性有限单元方法, 目前已建立了实用的二维分析 软件体系[ 4, 5] , 并在三维分析方面取得了进展[ 6] 。从 1999 年开始, 在水工结构的 p 型自适应分析方面 也有所突破, 1999 年, 程昭[ 7] 等人针对水工结构分析问题提出了三维升阶谱有限元分析方法。2001 年, 陈胜宏[ 8] 等人进一步提出了二维问题的 p 型自适应分析策略, 并将自适应有限元方法归类为全域 升阶方法、单元升阶方法和自由度升阶方法等三类。之后, 费文平[ 9] 等人将 p 型自适应有限元分析方
式中: i = ( - 1) i , i = ( - 1) [ i 2+ 0 5] , i = ( - 1) [ i 4+ 0 75] 。
图 2 变厚度的节理单元
棱基函数 ( p 2)
pEi =
1 4
(
1+
1i ) ( 1+
1i ) ( 1+
1i ) ( 2 i p ( ) +
2i p ( ) +
2i p ( ) ) ,
( i, j) F2
( i, j) F1
( 6)
i= 1
p= 2 i= 1
p = 4 i, j 2
同理可以写出 v 下 , v上 , w 下 , w 上 及 v, w 的具体表达式。
将基函数 N i , p Ei , p F 统一记为 i , 位移差 ( uN, i+ 4 - uNi ) , ( uE, i+ 4 - uEi ) , ( uF2 - uF1 ) 统一记为
,
)=
1 2
(
1
) p( )
面基函数 ( p 4)
式中: p ( x ) = 项式。
pF ( i, j ) ( , ) = i ( ) j ( ) ( i + j = p , i , j 2)
x
2p 2
1
Lp- 1 ( t ) dt , 而 Ln ( t ) =
-1
1 dn 2nn! dt n
( t 2 - 1) n
了三维节理单元升阶谱有限元法的解题步骤, 通 过具体 算例, 验证 了 p 型升 阶谱有 限元法 在求解 含三维 节理单
元的有限元问题时的可行性及优越性。
关键词: 节理单元; 升阶谱有限元; p 型有限元
中图分类号: TU457
文献标识码: A
有限元法是求解微分方程数值解的一种重要方法, 对于一个给定的问题, 为改善其有限元解的精 度, 可以采用以下 3 种方法。( 1) h 型有限元法[ 1] , 这种方法通过减小单元尺寸来提高有限元解的精 度。( 2) p 型有限元法[ 2] , 这种方法通过增加基底函数的阶次来提高有限元解的精度。( 3) hp 型有限 元法[ 3] , 这种方法是以上两种方法的综合, 它既减小单元尺寸, 又增加基底函数的阶次。作者所在的
( i , 1, 2, 12)
( 16)
式中:
2i = 1- [ i 12+ 0 65] , 2 i =
1 2
(
1-
(-
1) [ i- 1 4] ) ;
2i =
i- 1 8
,将
1i ,
1i ,
1i 用向量的形式
表达为 { 1} = { 0 0 0 0 - 1 - 1 1 1- 1 1 - 1 1} T , { 1} = { - 1 1- 1 1 0 0 0 0 - 1 - 1 1 1} T,
{
}= {
xz
yz
z }T =
1 e
{
u
v
w }T = [ B ]{ } e
( 11)
{ } = { xz yz z } T = [ D ] { } = [ D ] [ B ] { } e
( 12)
G0 0 式中: 弹性矩阵 [ D ] = 0 G 0
0 0 + 2G 根据虚功原理, 单元刚度矩阵为
58
拟, 可取如图 1 所示的六面体节理单元, 进行单元的网格划分。 1 1 形函数及位移函数 对节理单元上下面的相应点、棱和面可取 相同的点基函数、棱基函数和面基函数, 对无厚度节理单元或等厚 度节理单元, 不存在连续上、下两面的棱基函数和面基函数, 也不 存在体基函数。基底函数的具体形式如下[ 13] :
法推广到三维弹粘塑性领域。但是, 以上有关 p 型有限元的研究成果中均未涉及到断层、节理这一类 特殊单元。
大坝坝基、坝肩和岩石高边坡等部位总是存在断层、节理和软弱夹层等大规模的不连续面, 且对
结构的变形和稳定影响巨大, 故在有限元分析中应给予高度的重视。古德曼 ( Goodman) 最初运用有 限元技术模拟岩体工程中的非线性不连续面问题, 并提出了无厚度的节理元的概念[ 10] 。随后, 朱伯 芳于 1979 年提出了等厚度节理元模型, 并将其与无厚度节理元模型形成统一的计算公式[ 11] , 在此基 础上, 王鸿儒等人提出了变厚度的节理单元的弹塑性模型并将其应用到工程实践中[ 12] 。目前, 国内
, (p =
2, 3,
( 2)
( 3) ) 为 Legender 多
令位移函数为
4
p4
p i+ j = p
u下=
N iuNi +
pEi pu Ei +
F u p ( i, j ) p ( i, j ) F1
( 4)
i= 1
p= 2 i= 1
p= 4 j 2
4
p4
p i+ j = p
u上=
N iuN , i+ 4 +
f (p)
e
广义结点位移差
ui , 设单元基底函数个数为 f e ( p ) , 上式简化为
u= i= 1
i
ui ,
f ( p) e
f ( p) e
v=
i vi,
i=1
w=
i wi
i= 1
1 2 三维等厚度节理单元或无厚度节理单元升阶过程 当 p = 1 时: [ N ] 1= [ -
同理有 1I - 2I -
{ 1 } = { - 1 - 1 1 1- 1 1- 1 1 0 0 0 0} T。
面基函数 ( p 4)
F ( p ( i, j) 1, 2
,
,
,) =
1 2
(
1
) i( ) j ( )
( 17)
60
F ( p ( i, j) 3, 4
,
,
,) =
1 2
(
1
) i( ) j( )
( 18)
F ( p ( i, j) 5, 6
问题具有收敛速度快、计算精度高的优点。
1 三维 p 型无厚度节理单元和等厚度节理单元模型
对厚度很小和厚度变化不大的节理, 可以分别采用无厚度的节理单元和等厚度的节理单元进行模
收稿日期: 2001-11-30 基金项目: 国家自然科学基金 ( 59979021) 和教育骨干教师基金资助 作者简介: 费文平 ( 1972- ) , 男, 湖北天门人, 武汉大学博士生, 主要从事水工结构数值计算的研究。
部坐标系的 z 为中面的法线朝上方向, y 指向节理面的倾向, x 轴由右手法则确定, 并设等厚度节理
的倾角为 , 倾向为 。
59
三维等厚节理或无厚节理单元局部坐标与整体坐标的转换矩阵为
- sin
- cos
0
[ L] = cos cos - cos sin - sin
( 10)
sin cos - sin cos cos 2 4 三维等厚度节理的单元刚度矩阵
pEi
up E, i+ 4
+
F u p ( i, j ) p ( i, j ) F2
( 5)
i= 1
p= 2 i= 1
p = 4 i, j 2
4
p4
p i+ j= p
u=
N i ( uN , i+ 4 - uNi ) +
pEi p ( uE, i+ 4 - uEi ) +
F ( u - u ) p ( i, j) p
的形函数矩阵为
[ N ] p = [ [ N ] 1[ N ] p= 2[ N] p = 3 [ N] p = p ]
( 7)
1 3 坐标插值及坐标变换 对于节理上、下面坐标的插值仍采用各面上的四个节点进行插值, 即
4
4
4
x下 =
Nix i , y 下 =
N iy i , z 下 =
Nizi ,
( 8)
p i+ j= p 6
p i+ j + k= p