含节理单元的三维P型自适应有限元解法(1)

1、    含节理单元的三维P型自适应有限元解法(1)    本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节理单元升阶谱有限元法的解题步骤，通过具体算例，验证了p型升阶谱有限元法在求解含三维节理单元的有限元问题时的可行性及优越性。 关键词：节理单元 升阶谱有限元 p型有限元    有限元法是求解微分方程数值解的一种重方法，对于一个给定的问题，为改善其有限元解的精度，可以采用以下3种方法。(1)h型有限元法1，这种方法通过减小单元

2、尺寸来提高有限元解的精度。(2)p型有限元法2，这种方法通过增加基底函数的阶次来提高有限元解的精度。(3)hp型有限元法3，这种方法是以上两种方法的综合，它既减小单元尺寸，又增加基底函数的阶次。作者所在的研究小组从1995年开始研究水工结构的h型弹粘塑性有限单元方法，目前已建立了实用的二维分析软件体系4,5，并在三维分析方面取得了进展6。从1999年开始，在水工结构的p型自适应分析方面也有所突破，1999年，程昭7等人针对水工结构分析问题提出了三维升阶谱有限元分析方法。2001年，陈胜宏8等人进一步提出了二维问题的p型自适应分析策略，并将自适应有限元方法归类为全域升阶方法、单元升阶方法和自由度

3、升阶方法等三类。之后，费文平9等人将p型自适应有限元分析方法推广到三维弹粘塑性领域。但是，以上有关p型有限元的研究成果中均未涉及到断层、节理这一类特殊单元。    大坝坝基、坝肩和岩石高边坡等部位总是存在断层、节理和软弱夹层等大规模的不连续面，且对结构的变形和稳定影响巨大，故在有限元分析中应给予高度的重视。古德曼(Goodman)最初运用有限元技术模拟岩体工程中的非线性不连续面问题，并提出了无厚度的节理元的概念10。随后，朱伯芳于1979年提出了等厚度节理元模型，并将其与无厚度节理元模型形成统一的计算公式11，在此基础上，王鸿儒等人提出了变厚度的节理单元

4、的弹塑性模型并将其应用到工程实践中12。目前，国内外在有关p型自适应有限元分析的研究中，尚未涉及到这类特殊单元的处理问题，从而使研究成果的工程应用受到一定程度的限制。    对于常规块体单元的三维p型有限元模型，作者在文献9中已有详尽的论述，本文主给出三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型有限元模型，并给了升阶谱的计算格式。实例分析结果表明，用p型有限元法来求解含三维节理单元的有限元问题具有收敛速度快、计算精度高的优点。1 三维p型无厚度节理单元和等厚度节理单元模型    对厚度很小和厚度变化不大

5、的节理，可以分别采用无厚度的节理单元和等厚度的节理单元进行模拟，可取如图1所示的六面体节理单元，进行单元的网格划分。1.1 形函数及位移函数 对节理单元上下面的相应点、棱和面可取相同的点基函数、棱基函数和面基函数，对无厚度节理单元或等厚度节理单元，不存在连续上、下两面的棱基函数和面基函数，也不存在体基函数。基底函数的具体形式如下13：    点基函数(p1)         (1)式中：，。    棱基函数(p2)  

6、;   ，        (2)    面基函数(p4)     （i j=p, i,j2）    (3)式中：而为Legender多项式。令位移函数为         (4)        (5)  

7、;      (6)    同理可以写出V下，V上，W下，W上及v，w的具体表达式。    将基函数Ni,pEi,pF统一记为i，位移差(uN,i 4-uNi)，(uE,i 4-uEi),(uF2-uF1)统一记为广义结点位移差ui，设单元基底函数个数为fe(p)，上式简化为，同理有,1.2 三维等厚度节理单元或无厚度节理单元升阶过程 当p=1时：N1-1I-2I-3I-4I1I2I3I4I,I为3×3阶的单位阵。当p=2时，N2可在N1的基础上

8、进行扩充，扩充矩阵Np=2为Np=2=-5I-6I-7I-8I5            摘本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节         本篇论文是由3COME文档频道的网友为您在网络上收集整理饼投稿至本站的，论文版权属原作者，请不用于商业用途或者抄袭，仅供参考学习之用，否者后果自负，如果此文侵犯您的合法权益，请联系我们。I6I7

9、I8I。依此类推，可得最终的形函数矩阵为         (7)1.3 坐标插值及坐标变换 对于节理上、下面坐标的插值仍采用各面上的四个节点进行插值，即         (8)        (9)式中：(xi,yi,zi)i=18为节理间面体单元8个顶点的整体坐标。    定义整体坐标系的x轴朝北，y轴朝西，z轴朝上

10、，定义等厚度的节理单元或无厚度的节理的局部坐标系的z为中面的法线朝上方向，y指向节理面的倾向，x轴由右手法则确定，并设等厚度节理的倾角为，倾向为。    三维等厚节理或无厚节理单元局部坐标与整体坐标的转换矩阵为         (10)1.4 三维等厚度节理的单元刚度矩阵         (11)        (12

11、)式中：弹性矩阵；单元应变矩阵B'LB1/eLNp根据虚功原理，单元刚度矩阵为         (13)1.5 三维无厚度节理的单元刚度矩阵 当等厚度节理单元的厚度e0时，即形成无厚度的节理单元。此时，可假定单元内应力分量与位移差成正比，同理可得单元刚度矩阵为         (14)式中：为单元劲度矩阵。    2 三维p型变厚度的节理单元模型    

12、60;当节理单元的厚度变化较大时，应将等厚度节理单元推广得到变厚度的节理单元，如图2所示建立局部坐标系。2.1 形函数及位移函数 六面体变厚度的节理单元的基底函数，由点基函数、棱基函数、面基函数和体基函数组成，各基底函数的具体形式如下：点基函数(p1)     Ni(,)=1/8(1i)(1 i)(1 i) (i= 1，2，8)    (15)式中：i=(-1)i,i=(-1)i/2 0.5，i=(-1)i/4 0.75。    棱基函数(p2)  

13、60;  pEi=1/4(1 1i)(1 1i)(1 1i)(2ip() 2ip() 2ip()，(i,1,2,12)    (16)式中：2i=1-i/12 0.65，2i12(1-(-1)i-1/4)；2i=i-18，将1i，1i,1i用向量的形式表达为1=0000-1-111-11-11T，1-11-110000-1-11T，1-1-111-11-110000T。        面基函数(p4)      &#

14、160;  (17)        (18)        (19)式中：i,j2，i j=p。    体基函数(p6)：     (i,j,k2,i j k=p)    (20)令位移函数         (21)&#

15、160;           摘本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节         本篇论文是由3COME文档频道的网友为您在网络上收集整理饼投稿至本站的，论文版权属原作者，请不用于商业用途或者抄袭，仅供参考学习之用，否者后果自负，如果此文侵犯您的合法权益，请联系我们。      &#

16、160; 将基底函数Ni,pEi,pFi,pB统一记为i，位移uNi,uEi,uFi,uB统一记为广义结点位移ui，设单元基底函数的个数为fe(p)，上式简化为，同理。2.2 三维变厚度节理单元升阶过程 三维变厚度节理单元可按表1的方式进行升阶。    表1 基底函数与阶次的关系    阶次    1    2    3    4 

17、60;  5    6        基底函数    1-8    1-20    1-32    1-50    1-74    1-105    2.3 坐标插值与坐标变换 对

18、于单元任一点的坐标，仍采用单元的8个实结点进行插值，即，此处i=Ni(i=1,2 ,8),(xi,yi,zi)为立方体单元的8个顶点坐标。参见图2，取局部坐标系的xoy平面与变厚度六面体单元的中面重合，z轴垂直于中面。根据变厚度节理单元8个实结点的整体坐标，可按如下方法建立坐标转换矩阵及局部坐标系。    先求出变厚度六面体单元中面法向向量nx,y,z，其中x(y2 y6-y1-y5)(z3 z7-z1-z5)-(z2 z6-z1-z5）（y3 y7-y1-y5),y(z2 z6-z1-z5)(x3 x7-x1-x5)-(x2 x6-x1-x5)(z3

19、z7-z1-z5)，z(x2 x6-x1-x5)(y3 y7-y1-y5)-（y2 y6-y1-y5)(x3 x7-x1-x5)。    转换矩阵L的形成，令     l1=(x2 x6-x1-x5)/l,m1=(y2 y6-y1-y5)/l,n1=(z2 z6-z1-z5)/l    (22)     ,    (23)    向量l2 m2 n2T可根据

20、右手法则，由向量l1 m1 n1T和l3 m3 n3T确定，即     l2=m3n1-m1n3,m2=n3l1-n1l3,n2=l3m1-l1m3     (24)2.4 三维变厚度节理的单元刚度矩阵 变厚度节理单元的单元应变和应力分别为         (25)        (26)    变厚度节理单元的单元刚

21、度矩阵为         (27)    若用分块矩阵的形式来表示单元刚度矩阵，则         (28)3 含节理单元的三维p型自适应有限元算法3.1 控制方程与离散 根据实际情况分类进行单元的网格剖分，形成三种特殊单元及常规块体单元的单元刚度矩阵，并组合成总体刚度矩阵。形成荷载列向量，即         

22、;   摘本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节         本篇论文是由3COME文档频道的网友为您在网络上收集整理饼投稿至本站的，论文版权属原作者，请不用于商业用途或者抄袭，仅供参考学习之用，否者后果自负，如果此文侵犯您的合法权益，请联系我们。         (29)    式中右端三项分别为体力、

23、面力和集中力，并假设面荷载作用在1面上(余类推)，而Ex2, y2, z2,Ex2, y2, x2,Ex,x, y,y, z,z,可得到控制方程     KU=Rp    (30)3.2 升阶方法与误差估计9 常用的升阶分析方法有3种，即全域升阶法、自由度升阶法和单元升阶法。此处选用单元升阶法进行升阶计算，即对精度不满足求的单元进行升阶。在误差估计分析中，取高阶的应力和应变为应力和应变精确值的“最佳估计”。定义误差能范为         

24、;(31)总能范数为         (32)3.3 含节理单元的三维p型自适应有限元分析过程 设相对控制误差能范为et，对每一个单元进行误差估计，按下式计算单元的误差参数     (i=1,2,Ne)    (33)    若i>1，则将该单元进行升阶，误差估计结束后，重新进行弹粘塑性的有限元分析，如此反复，直到所有的单元均满足精度求。    4 算

25、 例4.1 算例1 为验证含节理单元的三维p型自适应有限元理论及程序的可行性和优越性，考察具有代表性的三维变厚度节理单元，取如图3所示的计算试件，对程序进行考核，单元及结点编号如图3。其中，单元、为块体单元，单元为变厚度节理单元，试件尺寸1m×1m×2m试件顶部作用一面力荷载-10.MPa，底部4结点固定结束，相对控制误差取3.0%，材料参数见表2。同时，将该试件进行网格细分，得到如图4所示的计算试件细分网格模型，进行常规有限元计算，并将p型有限元和各阶计算结果与细网格的常规有限元计算结点位移值进行比较，见表3。计算结果表明，在p=4时，p型有限单元法达到求的计算精度，随着

26、阶次的逐步升高，p型有限元的计算结果与细网格的常规有限元计算结果越来越接近，在p=4时两种有限单元法的计算结果相差甚微。这说明，p型有限元单元法在处理三维节理单元时是切实可行的。 同时可以看出，与常规有限元法相比，p型有限单元法对网格划分的求很低，用较少的单元数就能同等效果地模拟较复杂的常规有限单元网格，且具有很快的收敛速度。         表2 材料参数    单元    弹模/GPa    

27、泊松比    凝聚力/MPa    内摩擦角/(°)    块体节理    3015            4045     表3 两种有限单元法位移值比较(单位：mm)        粗网

28、格p型有限单元法    细网格常规有限单元法    p=1    p=2    p=3    p=4    计算误差(%)        /           

29、0;                摘本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节         本篇论文是由3COME文档频道的网友为您在网络上收集整理饼投稿至本站的，论文版权属原作者，请不用于商业用途或者抄袭，仅供参考学习之用，否者后果自负，如果此文侵犯您的合法权益，请联系我们。 

30、0;  5    uvw                        6    uvw              &#

31、160;         9    uvw                        10    uvw       

32、;             4.2 算例2 考察一混凝土重力坝，坝基中含一个变厚度的节理面节理1和一个等厚度的节理面节理2，并将坝体与坝基的胶结面视为一个无厚度的节理面节理3。坝高100m，坝顶宽10m，上游坡面垂直，下游坡面斜率10.75。坝基上游取1.5倍的坝高，下游取2倍的坝高，坝基的深度取2倍的坝高。再取10m的坝段宽构成三维六面体网格，单元网格如图5所示(x-z平面投影)，共160个单元，384个实结点，各部位的材料参数见表4。  &#

33、160;      表4 材料参数    参数    岩基    坝体    变厚度节理    等厚度节理    无厚度节理    容重/(MN/m3)弹模/GPa泊松比凝聚力/MPa内摩擦角/(°)流动系数抗拉强力/MPa倾角/(°)倾向度/(&#

34、176;)法向刚度/(MN/m3)切向刚度/(MN/m3)    0.026250.21.2451.00.2/    0.024280.1671.6401.00.18/    0.025240.21.0451.00.2/    0.025240.21.0450.10.233.69180/    /2.0451.00.201802510    以坝踵为坐标原点

35、，x轴指向下游，z轴垂直向上，y轴由右手法则确定，两侧坝面无y向位移，坝基上下游面无x向位移，坝基底面无z向位移。假设坝基地应力场由坝基重力场产生，坝体一次浇筑完成，水库一次蓄满，下游无水，取相对控制误差能范为5.0%。计算结果表明，当p=3时，计算过程收敛，升阶过程的主指标如表5所示。    整个升阶过程结束后的单元阶次分布图、屈服单元分布图、第一主应力、第三主应力、最大剪应力等值线图(单位：MPa)、位移矢量图和等安全系数图分别见图6图12所示。         表5 自

36、适应计算的主指标    升阶次数    升阶单元个数            摘本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节         本篇论文是由3COME文档频道的网友为您在网络上收集整理饼投稿至本站的，论文版权属原作者，请不用于商业用途或者抄袭，仅供参考学

37、习之用，否者后果自负，如果此文侵犯您的合法权益，请联系我们。    整体广义结点数    整体自由度    刚度矩阵稀疏度(%)    整体误差能范(%)    123    /16097    38412781939    70625994058  

38、0;         5 结 论    (1)在求解含三维节理单元的有限元问题时，p型有限单元法是一种切实有效的方法，它具有很高的收敛速度。(2)用p型有限单元法对含三维节理单元的水工结构有限元问题进行自 适应分析时，不必对单元网格重新剖分加密，在整个计算过程中网格可以保持不变。(3)用p型有限单元法求解含三维节理单元的水工结构有限元问题时，用较为简单的网格就模拟常规有限单元法中较为复杂的网格。(4)由于采用了升阶谱的单元模型，使得后一步的计算可以承袭前一步计算的结果，从而大大地减小了计算工作量。    参 考 文 献： 1O C Zienkiewicz, D