含节理单元的三维P型自适应有限元解法
时间：2022-10-10 11:00:30
摘要：本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节理单元升阶谱有限元法的解题步骤，通过具体算例，验证了p型升阶谱有限元法在求解含三维节理单元的有限元问题时的可行性及优越性。
关键词：节理单元 升阶谱有限元 p型有限元
有限元法是求解微分方程数值解的一种重要方法，对于一个给定的问题，为改善其有限元解的精度，可以采用以下3种方法。(1)h型有限元法[1]，这种方法通过减小单元尺寸来提高有限元解的精度。(2)p型有限元法[2]，这种方法通过增加基底函数的阶次来提高有限元解的精度。(3)hp型有限元法[3]，这种方法是以上两种方法的综合，它既减小单元尺寸，又增加基底函数的阶次。作者所在的研究小组从1995年开始研究水工结构的h型弹粘塑性有限单元方法，目前已建立了实用的二维分析软件体系[4,5]，并在三维分析方面取得了进展[6]。从1999年开始，在水工结构的p型自适应分析方面也有所突破，1999年，程昭[7]等人针对水工结构分析问题提出了三维升阶谱有限元分析方法。2001年，陈胜宏[8]等人进一步提出了二维问题的p型自适应分析策略，并将自适应有限元方法归类为全域升阶方法、单元升阶方法和自由度升阶方法等三类。之后，费文平[9]等人将p型自适应有限元分析方法推广到三维弹粘塑性领域。但是，以上有关p型有限元的研究成果中均未涉及到断层、节理这一类特殊单元。
大坝坝基、坝肩和岩石高边坡等部位总是存在断层、节理和软弱夹层等大规模的不连续面，且对结构的变形和稳定影响巨大，故在有限元分析中应给予高度的重视。古德曼(Goodman)最初运用有限元技术模拟岩体工程中的非线性不连续面问题，并提出了无厚度的节理元的概念[10]。随后，朱伯芳于1979年提出了等厚度节理元模型，并将其与无厚
度节理元模型形成统一的计算公式[11]，在此基础上，王鸿儒等人提出了变厚度的节理单元的弹塑性模型并将其应用到工程实践中[12]。目前，国内外在有关p型自适应有限元分析的研究中，尚未涉及到这类特殊单元的处理问题，从而使研究成果的工程应用受到一定程度的限制。
对于常规块体单元的三维p型有限元模型，作者在文献[9]中已有详尽的论述，本文主要给出三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型有限元模型，并给了升阶谱的计算格式。实例分析结果表明，用p型有限元法来求解含三维节理单元的有限元问题具有收敛速度快、计算精度高的优点。
1 三维p型无厚度节理单元和等厚度节理单元模型
对厚度很小和厚度变化不大的节理，可以分别采用无厚度的节理单元和等厚度的节理单元进行模拟，可取如图1所示的六面体节理单元，进行单元的网格划分。
1.1 形函数及位移函数 对节理单元上下面的相应点、棱和面可取相同的点基函数、棱基函数和面基函数，对无厚度节理单元或等厚度节理单元，不存在连续上、下两面的棱基函数和面基函数，也不存在体基函数。基底函数的具体形式如下[13]：
点基函数(p≥1)
(1)
式中：，。
棱基函数(p≥2)
，
(2)
面基函数(p≥4)
（i+j=p, i,j≥2）
(3)
式中：而为Legender多项式。
令位移函数为
(4)
(5)
(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时：[N]1＝[-φ1I-φ2I-φ3I-φ4Iφ1Iφ2Iφ3Iφ4I],I为3×3阶的单位阵。当p=2时，[N]2可在[N]1的基础上进行扩充，扩充矩阵[ΔN]p=2为[ΔN]p=2=[-φ5I-φ6I-φ7I-φ8Iφ5Iφ6Iφ7Iφ8I]。依此类推，可得最终的形函数矩阵为
(7)
1.3 坐标插值及坐标变换 对于节理上、下面坐标的插值仍采用各面上的四个节点进行插值，即
(8)
(9)
式中：(xi,yi,zi)i=1～8为节理间面体单元8个顶点的整体坐标。
定义整体坐标系的x轴朝北，y轴朝西，z轴朝上，定义等厚度的节理单元或无厚度的节理的局部坐标系的z′为中面的法线朝上方向，y′指向节理面的倾向，x′轴由右手法则确定，并设等厚度节理的倾角为α，倾向为β。
三维等厚节理或无厚节理单元局部坐标与整体坐标的转换矩阵为
(10)
1.4 三维等厚度节理的单元刚度矩阵
(11)
(12)
式中：弹性矩阵；单元应变矩阵[B']＝[L][B]＝1/e[L][N]p
根据虚功原理，单元刚度矩阵为
(13)
1.5 三维无厚度节理的单元刚度矩阵 当等厚度节理单元的厚度e0时，即形成无厚度的节理单元。此时，可假定单元内应力分量与位移差成正比，同理可得单元刚度矩阵为
(14)
式中：[λ′]为单元劲度矩阵。
2 三维p型变厚度的节理单元模型
当节理单元的厚度变化较大时，应将等厚度节理单元推广得到变厚度的节理单元，如图2所示建立局部坐标系。
2.1 形函数及位移函数 六面体变厚度的节理单元的基底函数，由点基函数、棱基函数、面基函数和体基函数组成，各基底函数的具体形式如下：点基函数(p≥1) Ni(ξ,η,ζ)=1/8(1＋ξiξ)(1+ηiη)(1+ζiζ) (i= 1，2，…，8)
(15)
式中：ξi=(-1)i,ηi=(-1)[i/2+0.5]，ζi=(-1)[i/4+0.75]。
棱基函数(p≥2)
pEi=1/4(1+ξ1iξ)(1+η1iη)(1+ζ1iζ)(ξ2iΦp(ξ)+η2iΦp(η)+ζ2iΦp(ζ))，(i,1,2,…12)
(16)
式中：ξ2i=1-[i/12+0.65]，η2i＝12(1-(-1)[i-1/4])；ζ2i=i-18，将ξ1i，η1i,ζ1i用向量的形式表达为{ξ1｝=｛0000-1-111-11-11}T，{η1}＝{-11-110000-1-11}T，{ζ1}＝{-1-111-11-110000}T。
面基函数(p≥4)
(17)
(18)
(19)
式中：i,j≥2，i+j=p。
体基函数(p≥6)：
(i,j,k≥2,i+j+k=p)
(20)
令位移函数
(21)
将基底函数Ni,pEi,pFi,pB统一记为φi，位移uNi,uEi,uFi,uB统一记为广义结点位移ui，设单元基底函数的个数为fe(p)，上式简化为，同理。
2.2 三维变厚度节理单元升阶过程 三维变厚度节理单元可按表1的方式进行升阶。
表1 基底函数与阶次的关系
阶次
1
2
3
4
5
6
……
基底函数
1-8
1-20
1-32
1-50
1-74
1-105
……
(25)
(26)
变厚度节理单元的单元刚度矩阵为?
(27)
若用分块矩阵的形式来表示单元刚度矩阵，则
(28)
3 含节理单元的三维p型自适应有限元算法
3.1 控制方程与离散 根据实际情况分类进行单元的网格剖分，形成三种特殊单元及常规块体单元的单元刚度矩阵，并组合成总体刚度矩阵。形成荷载列向量，即
(29)
式中右端三项分别为体力、面力和集中力，并假设面荷载作用在ξ＝1面上(余类推)，而Eη＝x2,η+y2,η+z2,η,Eζ＝x2,ζ+y2,ζ+x2,ζ,Eηζ＝x,ηx,ζ+y,ηy,ζ+z,ηz,ζ,可得到控制方程
[K]{U}={Rp}
(30)
3.2 升阶方法与误差估计[9] 常用的升阶分析方法有3种，即全域升阶法、自由度升阶法和单元升阶法。此处选用单元升阶法进行升阶计算，即对精度不满足要求的单元进行升阶。在误差估计分析中，取高阶的应力和应变为应力和应变精确值的“最佳估计”。定义误差能范为
(31)
总能范数为
(32)
3.3 含节理单元的三维p型自适应有限元分析过程 设相对控制误差能范为et，对每一个单元进行误差估计，按下式计算单元的误差参数
(i=1,2,……Ne)
(33)
若ξi>1，则将该单元进行升阶，误差估计结束后，重新进行弹粘塑性的有限元分析，如此反复，直到所有的单元均满足精度要求。
4 算 例
4.1 算例1 为验证含节理单元的三维p型自适应有限元理论及程序的可行性和优越性，考察具有代表性的三维变厚度节理单元，取如图3所示的计算试件，对程序进行考核，单元及结点编号如图3。其中，单元①、②为块体单元，单元③为变厚度节理单元，试件尺寸1m×1m×2m?试件顶部作用一面力荷载＝-10.MPa，底部4结点固定结束，相对控制误差取3.0%，材料参数见表2。同时，将该试件进行网格细分，得到如图4所示的计算试件细分网格模型，进行常规有限元计算，并将p型有限元和各阶计算结果与细网格的常规有限元计算结点位移值进行比较，见表3。
计算结果表明，在p=4时，p型有限单元法达到要求的计算精度，随着阶次的逐步升高，p型有限元的计算结果与细网格的常规有限元计算结果越来越接近，在p=4时两种有限单元法的计算结果相差甚微。这说明，p型有限元单元法在处理三维节理单元时是切实可行的。 同时可以看出，与常规有限元法相比，p型有限单元法对网格划分的要求很低，用较少的单元数就能同等效果地模拟较复杂的常规有限单元网格，且具有很快的收敛速度。
表2 材料参数
单元
弹模/GPa
泊松比
凝聚力/MPa
内摩擦角/(°)
块体
节理
30
15
0.2
0.15
2.8
1.6
40
45
表3 两种有限单元法位移值比较(单位：mm)
粗网格p型有限单元法
细网格常规
有限单元法
p=1
p=2
p=3
p=4
计算误差(%)
/
7.12
3.44
2.02
5
u
v
w
-0.004325
-0.004255
-0.032110
-0.003284
-0.003139
-0.033100
-0.003646
-0.003505
-0.033330
-0.003519
-0.003377
-0.032760
-0.003515
-0.003383
-0.032930
6
u
v
w
0.004437
-0.004339
-0.032230
0.003210
-0.003126
-0.033030
0.003646
-0.003553
-0.033370
0.003480
-0.003374
-0.032720
0.003485
-0.003393
-0.032910
9
u
v
w
-0.004322
-0.003632
-0.038410
-0.003915
-0.003213
-0.039310
-0.004036
-0.003337
-0.039360
-0.004077
-0.003368
-0.039150
-0.004055
-0.003361
-0.039220
10
u
v
w
0.002973
-0.003545
-0.044920
0.002580
-0.003242
-0.045500
0.002714
-0.003341
-0.045610
0.002728
-0.003341
-0.045610
0.002720
-0.003355
-0.045500
4.2 算例2 考察一混凝土重力坝，坝基中含一个变厚度的节理面节理1和一个等厚度的节理面节理2，并将坝体与坝基的胶结面视为一个无厚度的节理面节理3。坝高100m，坝顶宽10m，上游坡面垂直，下游坡面斜率1∶0.75。坝基上游取1.5倍的坝高，下游取2倍的坝高，坝基的深度取2倍的坝高。再取10m的坝段宽构成三维六面体网格，单元网格如图5所示(x-z平面投影)，共160个单元，384个实结点，各部位的材料参数见表4。