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

　　含节理单元的三维P型自适应有限元解法论文 摘要：本文提出了含三维无厚度节理单元、等厚度节理单元和变厚度节理单元的p型自适应有限元模型，给出了三维节理单元升阶谱有限元法的解题步骤，通过具体算例，验证了p型升阶谱有限元法在求解含三维节理单元的有限元问题时的可行性及优越性。 关键词：节理单元 升阶谱有限元 p型有限元 有限元法是求解微分方程数值解的一种重要方法，对于一个给定的问题，为改善其有限元解的精度，可以采用以下3种方法。(1)h型有限元法1，这种方法通过减小单元尺寸来提高有限元解的精度。(2)p型有限元法2.freelan)最初运用有限元技术模拟岩体工程中的非线性不连续面问题，并提出了无厚度的节理元的概念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上，Pa，底部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 ，坝顶宽10m，上游坡面垂直，下游坡面斜率1∶0.75。坝基上游取1.5倍的坝高，下游取2倍的坝高，坝基的深度取2倍的坝高。再取10m的坝段宽构成三维六面体网格，单元网格如图5所示(x-z平面投影)，共160个单元，384个实结点，各部位的材料参数见表4。 表4 材料参数 参数 岩基 坝体 变厚度节理 等厚度节理 无厚度节理 容重/(MN/m3) 弹模/GPa 泊松比 凝聚力/MPa 内摩擦角/(°) 流动系数 抗拉强力/MPa 倾角/(°) 倾向度/(°) 法向刚度/(MN/m3) 切向刚度/(MN/m3) 0.026 25 0.2 1.2 45 1.0 0.2 / / / / 0.024 28 0.167 1.6 40 1.0 0.18 / / / / 0.025 24 0.2 1.0 45 1.0 0.2 / / / / 0.025 24 0.2 1.0 45 0.1 0.2 33.69 180 / / / / / 2.0 45 1.0 0.2 0 180 25 10 以坝踵为坐标原点，x轴指向下游，z轴垂直向上，y轴由右手法则确定，两侧坝面无y向位移，坝基上下游面无x向位移，坝基底面无z向位移。假设坝基地应力场由坝基重力场产生，坝体一次浇筑完成，水库一次蓄满，下游无水，取相对控制误差能范为5.0%。计算结果表明，当p=3时，计算过程收敛，升阶过程的