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

摘要：本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的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)

式中：，。
　　棱基函数(p≥2)

　　面基函数(p≥4)

式中：而为Legender多项式。
令位移函数为

　　同理可以写出V_下，V_上，W_下，W_上及Δv，Δw的具体表达式。
　　将基函数Ni，^pE_i，^pF统一记为φ_i，位移差(u_N，_i+4-u_Ni)，(u_E，i+4-u_Ei)，(u_F2-u_F1)统一记为广义结点位移差Δu_i，设单元基底函数个数为f_e(p)，上式简化为，同理有，


1.2 三维等厚度节理单元或无厚度节理单元升阶过程 当p=1时：[N]₁＝[-φ₁I-φ₂I-φ₃I-φ₄Iφ₁Iφ₂Iφ₃Iφ₄I]，I为3×3阶的单位阵。当p=2时，[N]₂可在[N]₁的基础上进行扩充，扩充矩阵[ΔN]_p=2为[ΔN]_p=2=[-φ₅I-φ₆I-φ₇I-φ₈Iφ₅Iφ₆Iφ₇Iφ₈I]。依此类推，可得最终的形函数矩阵为

1.3 坐标插值及坐标变换 对于节理上、下面坐标的插值仍采用各面上的四个节点进行插值，即

式中：(x_i，y_i，z_i)_i=1～8为节理间面体单元8个顶点的整体坐标。

　　定义整体坐标系的x轴朝北，y轴朝西，z轴朝上，定义等厚度的节理单元或无厚度的节理的局部坐标系的z′为中面的法线朝上方向，y′指向节理面的倾向，x′轴由右手法则确定，并设等厚度节理的倾角为α，倾向为β。
　　三维等厚节理或无厚节理单元局部坐标与整体坐标的转换矩阵为

1.4 三维等厚度节理的单元刚度矩阵

式中：弹性矩阵；单元应变矩阵[B']＝[L][B]＝1/e[L][N]_p
根据虚功原理，单元刚度矩阵为

1.5 三维无厚度节理的单元刚度矩阵 当等厚度节理单元的厚度e→0时，即形成无厚度的节理单元。此时，可假定单元内应力分量与位移差成正比，同理可得单元刚度矩阵为

式中：[λ′]为单元劲度矩阵。

式中：ξ_i=(-1)ⁱ，η_i=(-1)^[i/2+0.5]，ζ_i=(-1)^[i/4+0.75]。
　　棱基函数(p≥2)

式中：ξ_2i=1-[i/12+0.65]，η_2i＝12(1-(-1)^[i-1/4])；ζ_2i=i-18，将ξ_1i，η_1i，ζ_1i用向量的形式表达为{ξ₁｝=｛0000-1-111-11-11}^T，{η₁}＝{-11-110000-1-11}_T，{ζ₁}＝{-1-111-11-110000}^T。
　　
　　面基函数(p≥4)

式中：i，j≥2，i+j=p。
　　体基函数(p≥6)：

令位移函数

　　将基底函数N_i，^pE_i，^pF_i，^pB统一记为φ_i，位移u_Ni，u_Ei，u_Fi，u_B统一记为广义结点位移u_i，设单元基底函数的个数为f_e(p)，上式简化为，同理。

2.2 三维变厚度节理单元升阶过程 三维变厚度节理单元可按表1的方式进行升阶。

2.3 坐标插值与坐标变换 对于单元任一点的坐标，仍采用单元的8个实结点进行插值，即，此处φ_i=N_i(i=1，2 ，…8)，(x_i，y_i，z_i)为立方体单元的8个顶点坐标。参见图2，取局部坐标系的x′oy′平面与变厚度六面体单元的中面重合，z′轴垂直于中面。根据变厚度节理单元8个实结点的整体坐标，可按如下方法建立坐标转换矩阵及局部坐标系。

　　先求出变厚度六面体单元中面法向向量n＝{x，y，z}，其中x＝[(y₂+y₆-y₁-y₅)(z₃+z₇-z₁-z₅)-(z₂+z₆-z₁-z₅）（y₃+y₇-y₁-y₅)]，y＝[(z₂+z₆-z₁-z₅)(x₃+x₇-x₁-x₅)-(x₂+x₆-x₁-x₅)(z₃+z₇-z₁-z₅)]，z＝[(x₂+x₆-x₁-x₅)(y₃+y₇-y₁-y₅)-（y₂+y₆-y₁-y₅)(x₃+x₇-x₁-x₅)]。

　　转换矩阵[L]的形成，
令

　　向量{l₂ m₂ n₂}^T可根据右手法则，由向量{l₁ m₁ n₁}^T和{l₃ m₃ n₃}^T确定，即

2.4 三维变厚度节理的单元刚度矩阵 变厚度节理单元的单元应变和应力分别为

　　变厚度节理单元的单元刚度矩阵为

　　若用分块矩阵的形式来表示单元刚度矩阵，则

3 含节理单元的三维p型自适应有限元算法

3.1 控制方程与离散 根据实际情况分类进行单元的网格剖分，形成三种特殊单元及常规块体单元的单元刚度矩阵，并组合成总体刚度矩阵。形成荷载列向量，即

　　式中右端三项分别为体力、面力和集中力，并假设面荷载作用在ξ＝1面上(余类推)，而E_η＝x²，_η+y²，_η+z²，_η，E_ζ＝x2，_ζ+y²，_ζ+x²，_ζ，E_ηζ＝x，_ηx，_ζ+y，_ηy，_ζ+z，_ηz，_ζ，可得到控制方程

3.2 升阶方法与误差估计^[9] 常用的升阶分析方法有3种，即全域升阶法、自由度升阶法和单元升阶法。此处选用单元升阶法进行升阶计算，即对精度不满足要求的单元进行升阶。在误差估计分析中，取高阶的应力和应变为应力和应变精确值的“最佳估计”。定义误差能范为

总能范数为

3.3 含节理单元的三维p型自适应有限元分析过程 设相对控制误差能范为e_t，对每一个单元进行误差估计，按下式计算单元的误差参数

　　若ξ_i&>1，则将该单元进行升阶，误差估计结束后，重新进行弹粘塑性的有限元分析，如此反复，直到所有的单元均满足精度要求。

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

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

5 结 论

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

［1］O C Zienkiewicz， D V Phillips. An automatic mesh generation scheme for plane and curved surface by isoparametric coordinates[J]. Int.J.Num. Meth. Eng.，1971，3:519-528.
［2］O C Zieniewicz， J P DE， R. Gago，D. W. Kelly. The hierarchical concept in finite element analysis[J]. Computers && structures，1983， 16(1-4):53-65.
［3］Zienkiewicz O C， Zhu J Z， Gong N G. Effective and practical hp version adaptive analysis procedures for the finite element method [J].Int. J. Numer. Meth. Eng.，1989，28(1-6):879-891.
［4］陈尚法，陈胜宏.弹粘塑性自适应有限元在滑坡稳定分析中的应用［J］.岩石力学与工程学报，2001，20(增2)：1596-1599.
［5］陈胜宏，王劲松，张君禄.水工结构的弹粘塑性自适应有限元分析［J］.水利学报，1996，(2)：68-75.
［6］Xia H X， Chen S H. 3-D adaptive FEM in rock slope stabilityanalysis ［J］. In:Proc. of the 10th Int. Conf. on Compu. Meth. and Advs. in Geomech.，Arizona， U.S.A:A.Bolkema，2001.
［7］程昭，陈胜宏.水工结构的三维阶谱有限元分析［J］.水利学报，1999，(12)：53-58.
［8］陈胜宏，程昭.水工结构分析的p型自适应有限单元法研究［J］.水利学报，2001，(11)：62-69.
［9］费文平，陈胜宏.水工结构的三维p型弹粘塑性有限单元法［J］.水利学报，2003，(2)：.
［10］Goodman R E. Methods of geological engineering in discontinuous rocks ［M］. St， Paul:West Publishing Co.，1976.
［11］朱伯芳.有限单元法原理及应用［M］.北京：水利电力出版社，1979.
［12］王鸿儒，李步娟，丁德平.复杂岩基三维弹塑性分析［J］.岩土工程学报，1988，(2)：19-30.
［13］Babuska I， Griebel M， Pitkaranta J. The problem of selecting the shape functions for a ptyefinite element[J]. Int. J. Numer. Meth. Eng.，1989，28(7-9):1891-1908.