岩石力学的数值模拟(讲义)

精品文档 1欢迎下载 第 10 章 岩石力学的数值模拟 随着计算机软硬件技术的迅速发展 使岩石力学有了长足的进步 特别在岩石力学的 数值计算和模拟方面发展尤为迅速 使得许多岩石力学解析方法难于解决的问题得以重新 认识 正如钱学森在给中国力学学会 力学 迎接 21 世纪新的挑战 的一封信中对力学 发展趋势总结的那样 今日力学是一门用计算机计算去回答一切宏观的实际科学技术问题 计算方法非常重要 岩石力学和其他力学学科一样 需要数值计算方法并推动岩石力学的 发展 岩石介质不同于金属材料 在数值计算方面具有其独特的特点 205 1 岩石介质是赋存于地壳中的各向异性天然介质 2 岩石介质被众多的节理 裂缝等弱面所切割而呈现高度的非均质性 而其物理 化学及力学性质具有随机性特点 3 岩石介质赋存时以受压为主 而且抗压强度远大于抗拉强度 4 岩石力学与工程问题在时空分布上较广 从本质上讲都是三维问题 5 岩石工程一般无法进行原型试验 而实验室测得的数据不能直接应用于工程设计 和计算 6 岩石力学与工程具有数据有限问题 数值计算方法经过几十年的发展 目前已形成许多种岩石力学计算方法 主要有有限 元法 边界元法 有限差分法 离散元法 流形元法 拉格朗日元法 不连续变形法及无 单元法等 它们各有优缺点 有限元的理论基础和应用比较成熟 在金属材料和构件的计 算中应用十分成功 但它是以连续介质为基础 似乎与岩体的非连续性有一定差距 流形 元等数值方法虽然考虑了岩体中节理效应 但其理论基础还不完全成熟 相信在不久的将 来 肯定会出现完全适合于岩体材料和工程的数值计算方法 206 208 10 1 岩石力学的有限元分析 209 213 有限元法 finite element method FEM 是岩石力学数值计算方法中最为广泛应用 的一种 自 20 世纪 50 年代发展至今 有限元已成功地求解了许多复杂的岩石力学与工程 问题 被广大岩石力学研究与工程技术人员喻为解决岩石工程问题的有效工具 有限元法 是根据变分原理求解数学物理方程的一种数值方法 有限元法把连续体离散成有限个单元 每个单元的场函数只包含有限个节点参量的简单场函数 这些有限个单元的场函数集合构 成整个结构连续体场函数 根据能量方程和加权函数方程可建立有限个求解参数的方程组 求解这些离散方程组 就是有限元法的精髓所在 虽然求解时把连续函数转化为求解有限 个离散点处的函数值 但只要单元划分得充分小时 足可以满足计算要求 有限元法求解问题时一般遵循以下步骤 1 有限元计算模型的建立 包括模型单元的划分 确定边界条件 2 对单元体进行力学分析 包括求解节点位移 单元应变和单元应力 3 对计算模型进行分析 4 进行计算分析 10 1 1 线弹性有限元法的基本方程 线弹性有限元是弹塑性有限元 损伤有限元 流变有限元等非线性有限元的基础 线 弹性有限元假定岩石介质连续 均质 小变形和完全弹性 有限元法求解弹性力学问题时通常以位移作为基本未知量 单元位移是以单元节点位 移为基本未知量 选择合理的位移插值函数 将单元位移表达为节点坐标的连续函数 插 值函数也可称为形函数 不同形状的单元具有不同的形函数 图 10 1 为三种最常见单元形式 即三角形 四边形及四面体单元 它们的形函数分别 精品文档 2欢迎下载 为 图 10 1 有限元的三种基本单元形式 a 三角形单元 b 四边形单元 c 四面体单元 三角形的形函数 2 1 mjiiycxba S N iiii 式中 S 为三角形面积 mii yxa kji yyb jki xxc 四边形的形函数 4 3 2 1 1 1 4 1 iN iii 式中 位移量为 iiu Nu iiv Nv iix Nx iiy Ny 四面体的形函数 6 1 mkjiizdycxba V N iiiii kkk mmm jjj i zyx zyx zyx a kk mm jj i zy zy zy b 1 1 1 kk mm jj i zx zx zx c 1 1 1 kk mm jj i yx yx yx d 1 1 1 式中 V 为四面体的体积 单元在直角坐标轴中位移分量分别为u v 因此单元的位移矩阵为w 10 1 eee Nwvuu 式中 为单元位移矩阵 N 为形函数矩阵 为单元节点位移列阵 e u e 根据几何方程 对位移矩阵求偏导数 可以得到应变矩阵 10 2 ee B 式中 B 为连续单元节点位移和单元应变的矩阵 也称为应变矩阵 对于三角形单元 B 为 常数矩阵 元素值取决于单元节点坐标差 根据本构方程 可以得到单元节点位移与单元应力矩阵之间的关系 10 3 eee BDD 式中 D 为弹性矩阵 应用虚功原理和最小势能原理可以推导出单元刚度矩阵的表达式 精品文档 3欢迎下载 10 4 v T e dvBDBK 各单元的体积力和面力按照静力的等效原则移置到各单元的节点上 其等效节点力为 10 5 v T e v T e dsQNQ dvPNP 式中 Pe 为作用于单元体积力 P 的等效节点荷载 Qe 为作用于单元面积力 Q 的等效 节点荷载 设环绕某节点i共有k个单元 则i节点上的外加荷载 Ri 为 10 6 i k i e si k i e vii PQRR 1 1 式中 Pi为作用于i节点上集中力 将各单元节点力与节点位移之间的关系叠加 形成以节点位移列阵为基本未知量的线 性代数方程组 10 7 n e e KK 1 10 8 RK 求解上式有限个线性代数方程组 得到节点位移矩阵 根据相应的节点位移利用式 10 2 和 10 3 计算单元的应力应变值 10 1 2 非线性问题的处理方法 从本质上讲 岩石力学问题都是非线性问题 这是因为一方面岩石材料的应力应变本 构关系绝大多数呈现非线性特征 另一方面岩体的变形又大多是大变形 对于求解岩石力 学的非线性问题 解析方法显得无能为力 而有限元可以求解岩石力学绝大部分的非线性 问题 岩石力学的非线性问题可以分为三大类 材料非线性 即岩石材料的本构关系为 非线性 而变形的几何关系仍是线性的 几何非线性 即岩石几何变形为非线性 而本 构关系仍为线性 两类非线性问题的组合 即岩石材料既是材料本构非线性 又是几何 非线性 这三类非线性问题总体平衡方程组的共同特征都是非线性方程组 可表示为 10 9 K 式中 为刚度矩阵 它是位移的函数 求解这类非线性方程组一般有三种方 K 法 迭代法 增量法和混合法 迭代法又称为牛顿 拉斐逊法 对于一个变量的函数 如图 10 2 它的迭代 fF 过程如下 设函数值F由F0增加至Fa通过切线法做第i次近似值可由下式确定 10 10 11 iaii FFK 式中 Ki 1为第i 1 次迭代时的曲线切线斜率 那么最终的解为 10 11 iii 1 牛顿 拉斐逊法的主要缺点是每次迭代过程中都要重新计算一下切线值 也就是刚度矩 阵及其逆矩阵 因此花费机时较长 为了避免每次计算Ki值 每次迭代时都采用同一个初 精品文档 4欢迎下载 始的K0 如图 10 3 此方法称为修正的牛顿 拉斐逊法 具体的计算公式如下 10 12 10 iai FFK 图 10 2 牛顿 拉斐逊迭代法 图 10 3 修正的牛顿 拉斐逊法 两种迭代法比较而言 修正的牛顿 拉斐逊法的迭代次数多于牛顿 拉斐逊法 但它省 去了每次迭代时重新计算赐度矩阵Ki 计算时间的多少 不仅取决于迭代次数或收敛速度 还取决于每次迭代所花费的时间 在一般情况下常刚度法在总体计算时间上比较合理 对 于非线性很强的材料 常刚度法并不适用 增量法也是把非线性问题线性化的一种处理方法 如图 10 4 把总荷载 F 划分成若干个 增量 逐级施加荷载增量进行求解 有限元计算公式为F 10 13 1iii FK 那么总位移和总荷载为 n i i 1 0 10 14 n i i FFF 1 0 增量法的误差较大 最终误差为各级增量时误差之和 为了减小误差 有时对一级增 量后 加上一个修正正项加以误差矫正 计算公式为 10 15 11 iiii FFK 图 10 4 增量法 由于增量法和迭代法各有其优缺点 目前有限元常常采用增量法和迭代法的结合 即混合法 一般将非线性关系分成若干个荷载增量段 在每一个增量段内用迭代法求解逼 近非线性解 最终累加求得各级荷载作用下的最终应力与应变 10 1 3 岩石弹塑性有限元分析 12 精品文档 5欢迎下载 岩石弹塑性本构关系是岩石主要非线性问题之一 岩石的弹塑性是指岩石材料的应力 应变关系在屈服之前呈线性关系 当应力达到屈服应力时 应力应变关系就变为非线性 由于弹塑性模型中应变不仅依赖于受载的应力状态 而且与加载路径有关 因此一般弹塑 性本构关系不能用应力应变全量关系准确描述 只能用能反映加载路径有关的应力应变增 量关系描述 在岩石非线性弹塑性本构关系有限元分析中一般采用初应力法和初应变法求解非线性 平衡方程组 初应力法是将荷载以微小增量形式逐级加在模型上 每加一级荷载增量 就会产生相应的位移增量 变形增量和应力增量 对于具有初 i dF i d i d i d 应力的弹塑性应力应变本构关系可以写成 10 16 0 ddDd p 式中 为初应力 为塑性矩阵 它与加载前的应力水平有关 而 dDd p p D 与应力增量无关 初应力法是通过对的处理将应力修正到正确的水平上 初应力不仅与加 0 d 0 d 载增量前应力水平有关 还与本级所加荷载增量引起的变形增量有关 如图 10 5 d 增量形式的平衡方程为 10 17 0 FddFdK 图 10 5 初应力法 式中 K0 为线弹性计算中的总刚度矩阵 为矫正荷载项 它由下式决定 F d 10 18 dVdDBFd p T 由于随位移变化而变化 因此计算时必须进行迭代求解 初应力法求解按照下述计 F d 算步骤实现 1 把全部荷载划分成若干个增量 在每一级增量段内 按照增量弹塑性平衡方程进 行求解 2 计算各单元的应力增量及当前的应用 精品文档 6欢迎下载 10 19 jijiji jiji jiji d dDd dBd 1 式中 下标i表示第i级荷载增量 j表示第j次迭代 3 根据岩石的屈服准则 由各单元应力判断单元是否屈服 对于塑性单元 计算应 力修正项并修正应力 10 20 jpjiji jipjp d dDd 4 塑性单元通过修正项计算等效节点力 所有塑性单元的等效节点力叠加 jp d 构成总的修正荷载矢量 10 21 dVdBFd jp T ji 5 在修正荷载作用下进行下次迭代运算 此时基本方程为 10 22 jiji dFdK 重复进行 2 5 步计算直至所有的塑性单元达到收敛精度要求 再进行下一步 的荷载增量计算 6 重新施加下一级荷载增量 重复计算 1 5 步 直至计算完毕 1 i dF 通过累加各级荷载作用的计算结果 求得总位移和总应力 u 一般初应力法的收敛速度比较缓慢 因此通常采用常刚度和变刚度法相结合的方法加 速收敛 初应变法按照弹塑性增量理论 总的应变量可表示为 10 23 pe ddd 式中 为弹性应变增量 为塑性应变增量 类似初应力法可以把塑性应 e d p d 变增量看做初应变 因为在计算过程中的计算格式和弹性系统中的初应变 p d p d 一致 0 d 精品文档 7欢迎下载 图 10 6 初应变法 从图 10 6 中可以看出 荷载增量dF所对应的位移为 按线弹性计算为 两 B du A du 者的位移增量之差称为初始位移 它是由材料塑性引起的附加位移 与模型系统初始位移 对应的单元位移为 那么单元的初应变为 0 d 10 24 00 dBd 10 1 4 岩石节理单元的有限元方法 岩体中常存在大量节理 而节理的变形性质和岩体力学性质有十分明显差别 从某种 程度上讲 节理的存在决定着岩体的力学性质 因此分析节理的变形性质 对于研究岩体 工程的变形情况至关重要 在进行有限元分析时 这种节理一般采用专门的节理单元来处 理 节理单元首先由 Goodman 提出 19 214 Goodman 认为节理不是一个实体 它只在若干 点处接触 因此节理单元也应满足这一特点 图 10 7 为无厚度节理单元 节点 1 与 4 2 与 3 具有相同的坐标 沿节理面的应力分别为法向应力和剪切应力 把节理单元两 n s 侧对应的位移定义为应变 相对的法向位移称为法向应变 相对的切向位移称为切v n 向应变 那么节理单元的本构关系为 s 10 25 D 图 10 7 无厚度节理单元 式中 为单元的刚度矩阵 对于节理单元长度上任何一点s处的应变可定义为界面两 D 侧相应位移之差 10 26 1 2314 vv L s vv L s n 10 27 1 2314 uu L s uu L s s 上式可用矩阵形式表示为 B 精品文档 8欢迎下载 根据一般化公式导出节理单元对应于局部坐标的刚度矩阵 10 28 nnnn ssss nnnn ssss nnnn ssss nnnn ssss e KKKK KKKK KKKK KKKK KKKK KKKK KKKK KKKK K 200020 020002 020200 002020 020200 002020 200020 020002 式中 分别为法向和切向刚度 对于整体坐标的刚度矩阵通过坐标变换矩阵得到 n K s K 具体如下 10 29 TKTK e T e 式中 T 为变换矩阵 它具人分块形式的对角阵 10 30 0000 TTTTT 式中 cossin sincos 0 T 这种无厚度节理单元适用于模拟直接接触的界面 而对于有一定厚度的弱面或夹层不 适宜 10 2 岩石力学的边界元分析 边界元法 boundary element method BEM 是和有限元法同步发展的一种数值计算 方法 与有限元相比有以下特点 边界元法把一个均质区域看做一个大单元 只把它的 边界离散化 区域内不划分单元 场变量处处连续 对于无限区域 场变量自动满足无 穷边界条件及自然表面状态 对于半无限区域 场变量也自动满足无穷远边界条件及自然 表面状态 有限元法是全区域离散化 而边界元是把基本方程转化为边界积分方程 只对 边界离散化建立相应的方程组进行求解 这样边界元使三维问题降为二维问题求解 使二 维问题转化为一维问题求解 当物体的表面积和体积之比较小时 边界元的划分单元数要 比有限元少数倍和十几倍 这样也使待解的方程数目 处理和存储的数据量降低同样的倍 数 大大节省了机时和算题费用 当仅需求解物体内部几个点的应力时 有限元仍不得不 划分整个区域 才能确定这几个点的应力值 而边界元当知道边界的应力解时 就可以根 据需要去求物体内部个别点的解 在应力梯度较高处 有限元法的剖分密度常常受到限制 而边界元可以方便地确定应力梯度的分布 但随着计算机硬件的飞速发展 边界元的优势 逐渐显得不很明显 边界元法矩阵方程中系数阵的元素结构比有限元法刚度阵中的元素复杂 有限元刚度 阵属带状稀疏阵 而边界元法的系数阵为满阵 因此对于面积和体积之比较大的薄壁结构 而言 边界元的优越性就不明显 边界元比较适合求解无限区域和半无限区域问题 如深埋巷道是一个典型的例子 但 边界元在计算非均质介质问题 非线性问题以及模拟工程开挖过程等方面不如有限元方便 有效 有限元与边界元划分单元的比较如图 10 8 精品文档 9欢迎下载 图 10 8 有限元法与边界元法比较 a 力学模型和边界条件 b 有限元单元划分 c 边界元单元划分 求解边界方程组主要有 Massonet 和 Rizze 分别提出的直接和间接两种数值解法 直接 法是直接建立关于边界的积分方程 通过离散化求解边界未知参数 并进而求解计算区域 内场函数值 间接法是设立一个在求解区域内含有若干系数的满足基本方程组的解 代入 边界条件求出系数 进而求得边界及区域内各点的场函数值 10 2 1 边界元分析原理 215 217 在无限的弹性体内作用一单位力引起周围的应力和位移的解析解是边界元求解弹性体 问题的基础 如图 10 9 在平面无限体内i点作用一x方向的单位力 其基本的弹性解在 1848 年由 Kelvin 解出 l k k lkl lklk n x r n x r v x r x r v n r rv 21 21 1 4 1 10 31 10 32 kl lklk x r x r r v vG u 1 ln 43 1 8 1 式中 为i点l方向单位力引起j点k方向的应力和位移 为 Kronecker 符 lk lk u lk 号 当l k时 当时 r为i j两点之间的距离 矢径 1 lkkl 0 lk 为j点作用面外法向对于k和l的方向余弦 为矢径r对外法向n的方向倒数 k n l n n r r x x r l l r x x r k k 当不考虑体力时 根据功的互等定理和以上的 Kelvin 基本解 可以建立弹性体边界上 i j两点之间的应力和位移关系 如图 10 10 10 33 uuuc ii 图 10 9 开尔文基本解条件 图 10 10 边界积分方程的建立 精品文档 10欢迎下载 式中 的大小取决于边界情况 当边界为光滑曲线时 等于 0 5 当边界曲线不光滑 i c i c 时 的值根据刚体位移原则求解 当边界作用分布力时 j点绕行一周 按叠加原理 i c 积分得 10 34 duduuc ii 式 10 34 就是边界元中的边界积分方程 10 2 2 边界单元与边界的离散 边界元是通过离散边界来求解边界积分方程 对于一个确定区域的边界进行离散 划分成有限个线段 称为边界单元 根据单元的节点数目和节点模式 边界单元可分为常 量单元 线性单元 二次单元等 常量单元及线性单元均为直线段 常量单元以单元的中 点为节点 每个单元有一个节点 场变量在单元内是一个常数 而线性单元的场变量在单 元内呈线性变化 单元的模式与有限元相似 可以表示成插值函数形式 即 10 35 m i iiu Nu 1 m i ii N 1 式中 m为单元的节点数目 插值函数由拉格朗日插值公式给出 为节点i处的值 i N i u 如图 10 11 几个常见的插值函数如下 常量单元 1 m1 i N 线性单元 2 m 1 2 1 1 N 1 2 1 2 N 二次单元 3 m 1 2 1 1 N 1 1 2 1 2 N 1 2 1 3 N 图 10 11 常见的几种边界单元 a 常量单元 b 线性单元 c 二次单元 对于常量的边界元 边界积分方程可简化为 10 36 duduuc nn ii 1 1 精品文档 11欢迎下载 对于平面问题 上式有 2n个方程 可写成矩阵形式 10 37 GuH 式中 分别为常量边界单元中点的位移列阵和应力列阵 G H 为系数矩阵 u 边界是元法的边界条件一般都为混合边界条件 即在一部分边界上位移和应力已知 而另一部分未知 为了方程求解 把所有的已知量移到等式的右边 未知量移到等式右边 经过变换后式 10 37 可简化为 10 38 FXA 式中 X 为包含所有位移和应力未知量列阵 F 为已知量列阵 A 为系数矩阵 求解此 方程可求得全部边界节点上未知量 由此进而求得计算域内任意一点的位移和应力值 如图 10 12 根据开尔文基本解和 Betti 的互等定理 可以得到计算域内任意一点的 位移和边界点外力功的关系式 10 39 duudui 图 10 12 计算域内i点与边界j点作用力与位移的关系图 边界离散后 对于常量单元 上式可改写为 10 40 duduu nn i 1 1 欲求在i点l方向的位移分量 可建立边界积分方程 10 41 duduu klkklk i 式中 分别为i点l方向作用单位力于j点k方向引起的应力位移开尔文基本解 lk lk u 欲求计算域内i点的应力可由几何方程和 Lame 方程求得 10 42 i j j i i i ji ij x u x u G x u v Gv 21 2 将边界积分方程代入上式可得 10 43 duSdD kkijkkijij 精品文档 12欢迎下载 对于二维问题 上式中系数 D S 分别为 10 44 2 21 1 4 1 kjikijikjjkikij rrrrrrv rv D 10 45 41 2 21 2 4 21 2 1 4 2 2 ijk lkilkjjikkijkji kjijijkijkij nv nnrrnvrrnrrnv rrrrvrv n r rv G S 式中 i i x r r j j x r r k k x r r 10 2 3 边界元与有限元耦合 如前所述 边界元和有限元在计算时各有优缺点 为了发挥各自的优点 提高求解精 度和解题效率 A Wexler 于 1972 年提出了边界元和有限元耦合求解的思路和方法 以后 J R Osias 和 M Margulies 等对边界元和有限元的耦合求解进行了深入的研究 为了讨论方便 把有限元的基本方程改写为如下形式 10 46 DFUK 式中 F 为边界力的等效节点力 D 为体力的等效节点力 由于边界力可以由面力 P 与面积 之积获得 那么有限元方程可改写为 10 47 DPMUK 一般边界元方程在考虑体力的情况下也可写为 10 48 BPGUH 从式 10 47 和 10 48 可以看出 有限元方程和边界元方程具有类似的形式 因 此从解题方式上提供了建立两者耦合的基本条件 把一个计算域人为地分成两个子域 对两个不同的子域分别采用边界元和有限元进行 求解 如图 10 13 有限元的边界为 S2 边界元的边界为 S1 两者的共同边界为 S3 对两个 区域分别建立有限元方程和边界元方程 再利用两者边界 S3上的平衡和协调条件建立耦合 方程组 耦合方程组有两种方式建立 把边界元方程转化为等价的有限元方程 把有 限元方程转化为等价的边界元方程 一般在耦合求解时前者采用得较多 图 10 13 边界元和有限元耦合的区域划分 边界元方程转化为有限元方程时 首先对边界元方程进行改造 对式 10 48 中 G 进 精品文档 13欢迎下载 行求逆可得 10 49 1 PBuHG 因为边界元法中 P 为分布力 要转化为有限元中的等效节点力 必须在上式两侧左乘一个 M 矩阵 那么可得到与有限元形式一致的方程 10 50 DFUK 式中 10 51 1 KHGM 1 DBGM FPM 需要注意的是有限元法的刚度矩阵为对称矩阵 而由边界元方程转化来的为非对 K 称 为了完全耦合 必须通过最小二乘法使其对称化 对称化后的刚度系数为 10 52 2 1 jiijij kkk 一个建筑结构在地基上沉降问题比较适合耦合求解 如图 10 14 把建筑结构部分划 分为有限元 地基部分为边界元 边界元的离散化有以下两种形式 由地面延伸到到无限 远处 边界元的划分须近似截取一个有限区 或利用半无限平面问题的基本解引入到无限 边界元 本算例采用耦合成有限元方程进行求解 计算结果和有限元结果比较如表 10 1 可以看出耦合法和有限元计算结果十分一致 但计算时间大为减少 表 10 1 结构顶部垂直位移计算结果比较 mm 有限元解 339 97135361600 耦合法解 350 105135370617 图 10 14 有限元和边界元耦合求解的算例 12 10 3 岩石力学有限差分方法 FLAC 有限差分方法 finite difference method FDM 是从一般的物理现象出发建立相应 的微分方程 经离散后得到差分方程 再进行求解的方法 这种方法可能是最古老的求解 方程组的数值方法之一 差分方程在计算机出现以前用一般的手摇计算器也可以求解 随 着计算机的不断发展和其他计算方法的兴起 有限差分法曾一度受到冷遇 但 20 世纪 80 年代末由美国 ITASCA 公司开发的 FLAC fast Lagrangian analysis of continua 程序 广泛采用差分方法进行求解 在岩土工程数值计算中得到了广泛的应用 10 3 1 FLAC 程序的简介 218 FLAC 是为岩土工程应用而开发的连续介质显式有限差分计算机软件 主要适用于模拟 计算岩土类工程地质材料的力学行为 特别是岩土材料达到屈服极限后产生的塑性流动 精品文档 14欢迎下载 材料通过单元和区域表示 根据研究对象的形状构成相应的网络结构 每个单元在外载和 边界约束条件作用下 按照约定的线性和非线性应力 应变关系产生力学响应 FLAC 软件 建立在拉格朗日算法基础上 特别适用于模拟材料的大变形和扭曲转动 FLAC 程序设有多 种本构模型 可模拟地质材料的高度非线性 包括应变软化和硬化 不可逆剪切破坏和压 密 黏弹性 蠕变 孔隙介质的流固耦合 热力耦合以及动力学行为等 另外 程序还设 有边界单元 可以模拟断层 节理和摩擦边界的滑动 张开和闭合行为 支护结构 如衬 砌 锚杆 可缩性支架或板壳等与围岩的相互作用也可以在 FLAC 程序中进行模拟 同时用 户可根据需要在 FLAC 中创建自己的本构模型 进行各种特殊的修正和补充 FLAC 程序具有强大的后处理功能 用户可以直接在屏幕上绘制或以文件形式创建和输 出多种形式的图形 使用者还可以根据需要 将若干个变量合并在同一幅图形中进行研究 分析 由于 FLAC 程序主要为地质工程应用而开发的岩石力学数值计算程序 包括反映地质材 料力学效应的特殊要求 FLAC 设计有 7 种材料本构模型 1 各向同性弹性材料模型 2 横观各向同性弹性材料模型 3 Coulomb Mohr 弹塑性材料模型 4 应变软化 硬化塑性材料模型 5 双屈服塑性材料模型 6 节理材料模型 7 空单元模型 可用来模拟地下开挖和煤层开采 FLAC 程序采用的是快速拉格朗日方法 基于显式差分来获得模型的全部运动方程 包 括内变量 的时间步长解 程序将计算模型划分为若干个不同形状的三维单元 单元之间 用节点相互连接 对某一个节点施加荷载之后 该节点的运动方程可以写成时间步长的有 限差分形式 在某一个微小的时间内 作用于该点的荷载只对周围的若干节点 相邻节点 有影响 根据单元节点的速度变化和时间 程序可求出单元之间的相对位移 进而可以求 出单元应变 根据单元材料的本构方程可求出单元应力 随着时间的推移 这一过程将扩 展到整个计算范围 直到边界 这样程序可以追踪模型从渐进破坏直至整体破坏的全过程 FLAC 程序将计算单元之间的不平衡力 并将此不平衡力重新加到各节点上 再进行下 一步的迭代运算 直到不平衡力足够小或者各节点的位移趋于平衡为止 图 10 15 为 FLAC 程序的计算循环示意图 FLAC 和有限元相比具有以下一些特点 1 由 Maxti 和 Cundall 于 1982 年提出的混合离散法适用于塑性破坏荷载和塑性流 动的模拟 这个方法比有限元中普遍采用的归约积分法更为合理 图 10 15 FLAC 程序中的计算循环示意图 精品文档 15欢迎下载 2 FLAC 虽然具有动力方程 但模拟系统实质上是静止的 这使得 FLAC 不需要数值 上卸载就可以遵循物理的失稳过程 3 FLAC 采用显式表达方式 对于任意非线性应力应变关系计算所用的时间和线性 关系基本一致 而且它并不需要存储任何矩阵 这就意味着计算机一般的内存就可以计算 大量的单元 而且大变形计算所花费的计算和小变形基本一样 因为它没有刚度矩阵 当然 FLAC 也有不足之处 用 FLAC 计算线性问题比同样的有限元要慢 FLAC 对非线 性 大变形问题及岩土物理失稳的计算更有效 FLAC 的计算时间由模拟系统最长固有周 期和最短固有周期之比确定 对于一些特定的模型 如弹性模量和单元尺寸变异较大的模 型 计算效率非常低 10 3 2 有限差分拉格朗日法的基本原理和方程 在有限差分方法中 控制方程组的每一个变量都可以直接用离散点的场变量代数形式 直接表达 其特点是直接求解基本方程和相应的定解条件近似解 具体步骤是首先将求解 域划分成网络 然后在网络节点上用差分方程近似表示微分方程 当采用较多节点时 近 似解的精度可以得到很大的提高 而有限元是将连续求解区域离散成一组有限个相互连接 的单元组合体 利用每一个单元内假设的近似函数来分片求解区域上待求的未知场函数 单元内的近似函数通常由未知场函数及其导数在单元各节点的数值和插值函数来表达 这 样有限元分析中的未知函数及其导数在各个节点上的数值就成为新的未知量 一经求出这 些未知量 就可以通过插值函数计算出各单元内场函数的近似值 进而得到整个求解域上 的近似解 有限差分和有限元两种方法都有一组代数方程需要求解 虽然那些方程是用不同的方 法推导而得 但其最终的形式相同 一般有限元程序把单元矩阵组成一个整体的总刚度矩 阵 而有限差分方法并不需要 因为在每一计算步都产生一个有限差分方程 FLAC 采用拉格朗日分析方法 由于它不需要形成整体刚度矩阵 大变形计算时在每一 个计算步都很容易修正坐标 位移增量施加到坐标上以致网格随着材料的移动和变形 这 就是所谓的拉格朗日法 与此对应的欧拉法 其材料移动和变形是相对于固定的网格 10 3 2 1 连续介质的场方程 对于一个二维的连续介质而言 其运动方程可表示如下 10 53 i ij ij i g xt u 式中 为介质密度 为坐标 为重力加速度 i x i g 在 FLAC 程序中 单元应变率是计算各主要参数的纽带 由于非线性应力应变本构关系 不具有惟一性 一般用增量形式表示 各向同性最简单的弹性本构关系张量差分形式如下 10 54 tGGK ijkkijnijtnij 2 3 2 式中 G K 分别为剪切和体积模量 为时间步 t 10 3 2 2 差分方程 三角形单元的差分方程一般由高斯离散定律推导而得 10 55 A i s i dA x f fdsn 式中 为沿封闭表面边界的积分 f为标量 矢量和张量 xi为位置矢量 为面积 s A 精品文档 16欢迎下载 A上的积分 ni为表面s上的单位法向量 把面积A上f的平均值定义为 10 56 Ad x f Ax f A ii 1 把上式代入式 10 55 就可以得到 10 57 s i i fdsn Ax f1 对于一个三角形单元 其有限差分形式可变为 10 58 Snf Ax f i si 1 式中 为三角形的边长 为每条边上的平均值 由引可以把应变率用节点速度来表S 示 10 59 i j j i ij x u x u 2 1 式中 10 60 s j b i a i j i Snuu Ax u 2 1 式中 标记 a 和 b 为三角形边的两个相邻节点 把应变率代入本构关系式 10 54 就可以获得应力张量 一旦应力计算出来 相应施加到各节点上的力也就确定 如图 10 16 每个三角形节 点由两个相邻边上应力合成得到 因此 10 61 2 1 2 2 1 1 SnSnF jjiji 对于包含两个三角形的四边形单元 由一个三角形形成的节点 其节点力是三角形两 边应力合成之和 有两个三角形形成的节点 其节点力是两者的平均值 划分网络中节点力是其周围所有单元对该节点的矢量和 它包含单元体重力和作 i F 用在节点上的外加荷载 如果单元体处于平衡或稳定状态流动 塑性流动 那么该节点力 将为零 否则该节点有一个加速度 10 62 m t Fuu t i tt i tt i 2 2 精品文档 17欢迎下载 图 10 16 三角形单元节点力的计算模型 a 具有速度矢量的三角形单元 b 节点力矢量 式中 上标表示相应变量计算的时间 对于大变形问题 上式再积分一次以确定网格节点 的新坐标 10 63 tuxx tt i t i tt i 2 有关 FLAC 计算实例请见文献 219 220 10 4 岩石力学的离散元分析 有限元法 边界元法和有限差分法都是基于连续介质力学的数值计算方法 它们都要 求计算对象满足变形的连续性条件 但工程岩体往往被节理和结构面所切割形成明显的节 理岩体 具有明显的不连续性 用连续介质力学的数值计算方法难于处理 针对不连续岩 体的变形和运动的求解 Cundall 于 1971 年提出了一种新的数值计算方法 离散单元法 distinct element method DEM 并用汇编语言编制了计算程序 1978 年 Cundall 又 将最初的汇编语言全部翻译成 FORTRAN 文本 形成离散元的基本程序 到 1985 年 他们完 成了目前广泛采用的离散元数值分析程序 UDEC universal distinct element code 离 散元法由王泳嘉教授等人于 20 世纪 80 年代中期介绍到我国 发展十分迅速 目前在矿业 铁道及水利等行业已被广泛应用 221 223 10 4 1 离散元法基本原理 有限元法的理论基础是基于最小势能原理 边界元法的理论基础是 Betti 互等定理 而离散元法的理论基础则为最简单 最基本的牛顿第二运动定律 它与其他数值方法不同 的是离散元法首先将计算区域划分成有限个独立的多边形块体单元 单元与单元之间具有 一定的初始接触状态 单元在外力或自重的作用下转动和移动 最终达到平衡状态或匀速 运动 离散元法划分的单元 依据力学性质可以分为刚性体和变形体 依据形状可分为多 边形和圆形 由于圆形变形体计算比较复杂 本节只介绍二维刚性体多边形计算模型 离散单元并不像有限单元必须符合三个基本方程 平衡方程 变形协调方程和本构方 程 只要符合平衡方程即可 对于变形体还须符合本构方程 离散单元之间的接触一般有 三种方式 角 角接触 边 角接触或边 边接触 如图 10 17 对于图中块体 B 其相邻 块体分别对块体 B 通过接触点作用一组力 1 5 加上块体 B 的重力 它们 xi F yi Fi 对块体 B 的重心产生合力 F 和合力矩 M 即 10 64 n i cixiciyi n i yiy n i xix yyFxxFM FF FF 1 1 1 精品文档 18欢迎下载 图 10 17 块体的集合及作用于个别块体上的力 式中 yi为块体受力作用点坐标 xc yc为块体重心坐标 如果上式中的合力和合力 i x 矩不等于零 那么块体 B 将按照牛顿第二定律的规律运动 令块体 B 的质量为m 转动惯 量为 I 块体绕其重心转动角度为 那么块体的运动方程为 10 65 I M g m F u m F u y y x x 式中 u为块体位移 g为重力加速度 对上式可采用差分方法进行求解 对x方向的 运动方程采用向前分格式进行数值积分 这样可以得到岩块质心沿x方向的速度和位 移 10 66 tututtu tututtu xxx xxx 00 00 式中 t0为初始时间 为计算时步 对于y方向的运动和转动 都有类似的算式 t 虽然离散单元可被视为不变形的刚体 但在单元接触的边界仍存在变形 如图 10 18 设块体之间的法向力 Fn 两块体之间的 叠合 位移为un 两者之间成正比关系 10 67 nnn uKF 式中 Kn为接触面的法向刚度系数 图 10 18 离散单元间的作用力 a 块体之间压缩 b 块体之间剪切 精品文档 19欢迎下载 因为块体所受的剪切力与块体运动和加载历史无关 所以对于剪切力和剪切位移关系 宜用增量形式表示 10 68 sss uKF 式中 Ks为接触面的剪切刚度系数 为两块体之间的相对位移 s u 上面两式所表示的力和位移是弹性关系 它的极限值符合 Coulomb Mohr 准则 即 10 69 tan ns FCF 式中 C 和分别是接触面的黏结力和内摩擦角 当块体接触面上的剪切力大于极限值时 块体之间就会产生滑动 而当块体受到张力相互分离时 块体产生张拉破坏 作用在块体 表面上的法向力和剪切力随之消失 块体的运动是一个不可逆过程 为了避免岩块在运动过程中发生振动 在离散元法计 算中采用加阻尼的方法来耗散系统在震动过程中的动能 在静力平衡问题中 加阻尼来吸 收系统的动能 以使系统达到稳定状态 阻尼模型一般可采用具有集中质量的 Voigt 模型 如图 10 19 模型阻尼系统的自由振动微分方程为 10 70 0 kuucum 图 10 19 具有集中质量的弹性阻尼系统 式中 c为阻尼系数 根据振动理论 临界阻尼系数为 在实际计算时 阻尼系数mk2 取值应取略小于临界值 以使岩块运动没有回弹 引入阻尼并考虑岩块位移后 离散单元 法的基本运动方程为 10 71 tftkutuctum 式中 m为单元质量 c是黏性阻尼系数 k是为弹性刚度系数 f是单元所受的外荷载 u为位移 t是时间 上式可用动态松弛法进行求解 图 10 20 表示动态松弛法求解力和位 移的计算循环 计算循环是以时步进行差分计算 在每一时步内进行一次迭代 根t t 据前一次迭代所得到的块体位置 求出接触力作为下一次迭代的出发点 如此进行反复迭 代 时步应选得比较小 以使每个单元在一个时步内 只能以很小的位移与其相邻单元作 用 而与较远的单元无关 这一计算循环的物理意义为 相互镶嵌排列的初始块体系统 在空间中有自己的固定位置 处于平衡状态 一旦外界外力或边界位移约束条件发生变化 时 某些块体在重力和外力作用下 按照运动方程产生一定的加速度并产生位移 因而使 块体在空间的位置发生变化 产生位移后块体与相邻块体在空间位置上发生重叠 根据块 精品文档 20欢迎下载 体力和位移关系 使更多的块体由于块体间的重叠而受力 于是又产生新的运动和位移 依次迭代下去 遍及整个块体集合 计算模拟各个块体位移和转动的整个过程 直到最后 收敛达到平衡状态或所需结果为止 图 10 20 时步的计算循环 10 5 岩石力学的流形元分析 实际的工程岩体常被一些节理等结构面所切割 形成一种处于连续与非连续之间的拟 连续介质 因此连续介质力学数值方法 如有限元 边界元 有限差分方法等 处理岩体 节理时采用节理单元方式 而非连续介质力学数值方法 离散元 非连续变形分析 discontinuous deformation analysis DDA 处理事先划分好的块体模型比较适宜 而新近兴起的流形元对于处理连续与非连续介质耦合问题比较适宜 对于岩石材料尤为适 合 流形元方法 manifold element method MEM 是我国著名留美学者石根华 林德璋 博士于 20 世纪 90 年代初提出并开发的一种全新数值计算方法 流形元法以拓扑流形学为 基础 应用有限覆盖技术 通过在计算区域内各物理覆盖上建立通用的覆盖函数和以加权 求和形成总体位移函数 从而把连续和非连续变形力学问题统一起来考虑 它以数值流形 为核心 在 DDA 基础上 联合有限元和解析法的连续分析方法 从而形成在一切空间 包 括有限元 DDA 和解析法 统一的计算形式 224 227 有限覆盖技术是在数学流形分析中经常采用的一种方法 石根华采用这一方法 统一 解决了连续和非连续变形的力学问题 有限覆盖技术构造的覆盖函数具有以下两个基本性 质 局部非零性 覆盖函数只在局部区域内不为零 即在局部区域内有解 在求解区 域内 覆盖函数之和为 1 即局部区域内所有覆盖函数组成的总体位移函数 它与数学流 形的区别在于数学流形在整个求解区域内总体位移函数处处连续并光滑可导 它可完全定 义而与覆盖无关 而数值流形的总体位移函数是在覆盖基础上定义的 且可分片微分 在 覆盖的接触面上是完全非连续的 总之 数值流形的基本特征是覆盖函数在局部区域内连 续 各分片区域之间覆盖函数不连续 通过采用连续和非连续覆盖函数的方法 把连续和 非连续变形的计算统一到数值流形中 10 5 1 数值流形方法的有限覆盖 流形元有两套有限覆盖技术 数学覆盖和物理覆盖 其实覆盖的概念和有限元的网格 基本一致 也就是说流形元有两层独立的网格 数学网格和物理网格 数学网格一般定义 数值解的精度 它由用户根据需要确定 而物理网格由材料的形状 裂隙 边界以及分区 截面等材料的物理和集合性质决定 用户不能随意改动和选择 因此物理覆盖Ui就是计算 边界 材料分区和裂隙等所组成的网格 数学覆盖Vi是数学网格节点所构成的单元 数学 覆盖可以是任何形状的 而物理覆盖完全由材料的物理和几何性质确定 流形元方法就是 把物理和数学两层覆盖重叠在一起 并用物理覆盖把数学覆盖重新划分 把数学覆盖重叠 入物理覆盖 形成新的连续和非连续耦合的有限覆盖系统 在这个系统中 物理覆盖代替 单元的节点 覆盖的并集交线代替单元的边界 覆盖重叠的交集代替单元 力边界条件 力 位移关系 运动方程 a F m 位移边界条件 位移 u力 F 精品文档 21欢迎下载 图 10 21 有限覆盖和流形单元的描述 219 a 两个块体的有限覆盖 b 块体内有一条裂隙的有限覆盖 图 10 21 分别采用四个有限单元网格覆盖一个四边形的材料区域 四个有限单元的节 点辚 1 2 3 2 4 5 2 5 3 3 5 6 图 10 21 a 表示在有 4 个原始 单元的数学覆盖 细线表示 上覆盖有两个块体的物理覆盖 粗线表示 的有限覆盖系统 图 10 21 b 为同一数学覆盖上有一条不连通裂隙的块体情况 图中大数字表示节点覆盖 码 小数字下标表示被不连续的物理覆盖划分的分区码 图中 11 12 21 22等分别表示 节点 1 2 的物理覆盖 10 5 2 流形无网格的覆盖函数和权函数 位移函数与材料边界无关 如果材料只占有单元的一部分 位移函数仍然是相同的 对三角形单元 相应三个节点有三个覆盖包含这个单元 因此每个单元是它三个节点的三 个覆盖公共区域 每个单元需要计算权函数 xi yi 表示为节点i 1 2 3 的坐标 则相应节点 e 1 e 2 e 3 的坐标和位移可表示如下 坐标 位移 1 e 1 1 1 1 eeee vuYX 2 e 2 2 2 2 eeee vuYX 3 e 3 3 3 3 eeee vuYX 流形元中单元任意一点I x y 的权函数为 10 72 y x fff fff fff yxW yxW yxW e e e 1 333231 232221 131211 3 2 1 式中 1 2 2 1 1 2 2 1 3 1 1 3 3 1 1 3 2 3 3 2 2 3 3 2 333231 232221 131211 1 eeeeeeee eeeeeeee eeeeeeee XXYYYXYX XXYYYXYX XXYYYXYX fff fff fff 10 73 10 74 式中 为节点i覆盖的权函数 其意义为加权的平均 它取节点位移的百分数 1 yxWe 精品文档 22欢迎下载 对三节点单元 三个节点的三个权函数之和为 1 权函数相当于有限元中的形函数 物理覆盖Ui上的覆盖函数为 可以是线性的 也可以是非线性的 yxvyxu ii 每一个节点覆盖生成的位移为覆盖函数与权函数的乘积 单元各点的最终位移则为整个物 理覆盖系统上的总体函数 如覆盖函数用标准二阶级数 那么总体函数为 10 75 3 21 3 21 i m j ieijjie i m j ieijjie yxWyxSv yxWyxSu yxv yxu 建立了每个流形单元上的位移函数以后 就可以按照有限单元方法形成单元的单刚及 总刚度矩阵进行求解 10 5 3 数值流形方法的平衡方程式 对