一种基于有限单元法的固体连续介质变形的数值模拟方法

摘要

本发明公开了基于有限单元法的固体连续介质变形的数值模拟方法，包括步骤(1)建立固体连续介质的有限元离散模型；步骤(2)在对固体连续介质施加预设边界条件的情况下，计算获得有限元离散模型的每个结点的不平衡力和结点局部刚度系数；步骤(3)筛选出有限元离散模型的最大结点不平衡力，并与预设收敛标准比较，判断固体连续介质是否达到平衡状态或处于不平衡流动状态，以决定进入下一步继续计算或停止计算；步骤(4)计算固体连续介质有限元离散模型中每个结点的位移；根据每个结点的位移，计算有限单元网格每个单元的应力；步骤(5)重复步骤(2)～步骤(4)，直到固体连续介质达到平衡状态或处于不平衡流动状态。

说明书

技术领域

本发明涉及固体连续介质变形分析技术领域，尤其涉及一种基于有限单元法的固体连续 介质变形的数值模拟方法。

背景技术

有限单元法(FEM)产生于上个世纪50年代，当结构力学和弹性力学中的其他近似方法 不适合求解连续结构时被提出。在随后的数十年，有限单元法获得广泛认可，受到密集研究， 并取得丰硕成果，已发展成为高度成功的一般连续介质问题的近似求解技术。有限单元法将 连续介质划分为有限单元组成的离散体系，基于变分原理或加权余量法，将连续域上的偏微 分支配方程转换为线性代数方程组，求解线性代数方程组得到域变量的近似数值解。其中， 有限单元(finiteelement)是指连续介质离散后具有一定尺寸的组成元件，有三角形、四边形、 四面体、八面体单元等。固体连续介质(solidcontinuum)是指在空间中连续无间断无空隙的 两维或三维固体，指自然界的土体和岩体，如土木工程结构的地基和地下洞室围岩等；也包 括人造的固体，如人工填筑的土石体，混凝土体，金属体和塑料体等。

在多种工程问题中，比如土木、水利工程设计中，建筑物岩土地基变形和稳定性分析， 地下空间岩土介质变形和稳定性分析，土石坝、混凝土坝变形、强度和稳定性分析，建筑物 地基和地下洞室岩土介质流变分析；材料工程中材料大变形和材料蠕变分析等，现有的有限 单元方法只给出平衡状态下的结果，不能给出变形和应力的演化过程，不能够追踪介质的准 静态力学运动过程。

发明内容

为了解决现有技术的缺点，本发明提供一种基于有限单元法的固体连续介质变形的数值 模拟方法。该方法将固体连续介质离散成有限单元网格，在准静态荷载作用下，固体连续介 质按运动步运动，从而追踪固体介质运动、变形和应力发展的过程。

本发明采用以下技术方案：

一种基于有限单元法的固体连续介质变形的数值模拟方法，包括：

步骤(1)：将固体连续介质离散成有限单元网格，建立固体连续介质的有限元离散模型；

步骤(2)：在对固体连续介质施加预设边界条件的情况下，根据介质应力和外部荷载计 算获得有限元离散模型的每个结点的不平衡力；根据有限元模型几何特性和介质弹性常数计 算获得有限元离散模型的每个结点的局部刚度系数；

步骤(3)：筛选出有限元离散模型的最大结点不平衡力，将有限元离散模型最大结点不 平衡力与预设收敛标准比较，判断固体连续介质是否达到平衡状态或处于不平衡流动状态， 以决定进入下一步继续计算或停止计算；

步骤(4)：根据有限元离散模型中每个结点的不平衡力与结点局部刚度系数的比值，计 算固体连续介质有限元离散模型中每个结点的位移；根据每个结点的位移，计算有限单元网 格每个单元的应力；

步骤(5)：重复步骤(2)～步骤(4)，直到有限元离散模型最大结点不平衡力小于等于 预设收敛标准，介质达到平衡状态，结束计算，并获得介质的变形和应力状态；或有限元离 散模型的最大结点不平衡力不再减小且保持大于预设收敛标准的某一非零数值，则固体连续 介质处于不平衡流动状态，结束计算或继续计算以追踪介质的流动过程。

所述步骤(2)中，固体连续介质有限元离散模型中结点不平衡力为有限元离散模型的外 荷载等效结点力和介质应力等效结点力的合力。

所述步骤(2)中，固体连续介质有限元离散模型中结点的局部刚度系数为围绕该结点的 所有单元的结点刚度系数之和。

所述步骤(4)中，固体连续介质的有限元离散模型中每个结点的位移，与有限元离散模 型中每个结点的不平衡力与结点局部刚度系数的比值成正比例关系。

所述步骤(4)中，固体连续介质有限元离散模型中结点位移，等于结点的不平衡力与结 点局部刚度系数的比值乘以折减系数。

所述折减系数的取值范围为大于等于0小于等于1.0。

本发明的有益效果为：

(1)本发明作为一种基于介质运动的求解方法，能够追踪介质的准静态力学运动过程， 而现有的有限单元法只给出平衡状态下的结果，不能给出变形和应力演化的过程，而在科学 和工程领域中，有时需要知道介质运动或变形的过程，在这些问题求解中本方法具有价值。

(2)本发明的方法求解高度材料非线性和几何非线性问题的能力较强，对于这类问题的 求解，现有的有限单元法常常遇到计算收敛困难或更耗费计算时间。

(3)本发明的方法不需要建立线性代数方程组和整体刚度矩阵，需要的计算机内存相对 较小，适合求解大规模问题。

(4)本发明的方法在计算模拟固体连续介质运动和变形方面具有重要价值。该方法可应 用于工程领域中的多种问题，如土木、水利工程中，建筑物岩土地基的变形和稳定性分析， 地下空间岩土介质的变形和稳定性分析，土石坝混凝土坝变形、强度和稳定性分析，建筑物 地基和地下空间岩土介质流变分析，材料工程中材料大变形和蠕变分析等。

附图说明

图1是本发明的基于有限单元法的固体连续介质变形的数值模拟方法流程图；

图2是半无限大介质表面受矩形均布荷载作用及计算范围示意图；

图3是岩体中开挖的地下洞室尺寸和介质计算范围示意图；

图4是金属圆柱体压缩示意图；

图5是金属圆柱体有限单元网格示意图；

图6a)是金属圆柱体在第5000运动步，压缩量为26.2772cm的一径向剖面变形示意图；

图6b)是金属圆柱体在第10000运动步，压缩量为38.5205cm的该径向剖面变形示意图；

图6c)是金属圆柱体在第20000运动步，压缩量为50.2398cm的该径向剖面变形示意图；

图6d)是金属圆柱体在第30000运动步，压缩量为54.8277cm的该径向剖面变形示意图；

图6e)是金属圆柱体在第140343运动步，压缩量为58.4678cm的该径向剖面变形示意 图；

图7是地表空间开挖示意图；

图8是地表空间开挖后竖向剖面MN结点位移矢量示意图。

具体实施方式

下面结合附图与实施例对本发明做进一步说明：

如图1所示，本发明的基于有限单元法的固体连续介质变形的数值模拟方法，包括：

步骤(1)：将固体连续介质离散成有限单元网格，建立固体连续介质的有限元离散模型；

步骤(2)：在对固体连续介质施加预设边界条件的情况下，根据介质应力和外部荷载计 算获得有限元离散模型的每个结点的不平衡力；根据有限元模型几何特性和介质弹性常数计 算获得有限元离散模型的每个结点的局部刚度系数；

步骤(3)：筛选出有限元离散模型的最大结点不平衡力，将有限元离散模型最大结点不 平衡力与预设收敛标准比较，判断固体连续介质是否达到平衡状态或处于不平衡流动状态， 以决定进入下一步继续计算或停止计算；

步骤(4)：根据有限元离散模型中每个结点的不平衡力与结点局部刚度系数的比值，计 算固体连续介质有限元离散模型中每个结点的位移；根据每个结点的位移，计算有限元网格 每个单元的应力；

步骤(5)：重复步骤(2)～步骤(4)，直到有限元离散模型最大结点不平衡力小于等于 预设收敛标准，介质达到平衡状态，结束计算，并获得介质的变形和应力状态；或有限元离 散模型的最大结点不平衡力不再减小且保持大于预设收敛标准的某一非零数值，则固体连续 介质处于不平衡流动状态，结束计算或继续计算以追踪介质的流动过程。

在对固体连续介质施加预设边界条件的情况下，连续介质有限元离散模型的外荷载等效 结点力设为^fF_i^l，其中，下标i＝1,2,3，代表三个空间坐标轴的方向；上标l表示结点编号； 上标f表示外部荷载。

有限单元法中，基于虚功原理或加权余量法，给出外荷载产生的单元等效结点力计算公 式。对于三维等参单元，计算公式为：

${\left\{F\right\}}^e=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\left[N\right]}_{\zeta =\pm 1}^T\left\{p\right\}Ad\xi d\eta ---\left(1\right)$

公式(1)和(2)分别给出单元边界面力{p}和体积力{f_b}产生的单元等效结点力。式中， ξ,η,和ζ为单元局部坐标；[N]为单元形函数矩阵；|J|是雅可比行列式；A是一个系数， 定义为：

$A=\sqrt{{(\frac{\partial x}{\partial \xi }\frac{\partial y}{\partial \eta }-\frac{\partial x}{\partial \eta }\frac{\partial y}{\partial \xi })}^2+{(\frac{\partial y}{\partial \xi }\frac{\partial z}{\partial \eta }-\frac{\partial y}{\partial \eta }\frac{\partial z}{\partial \xi })}^2+{(\frac{\partial z}{\partial \xi }\frac{\partial x}{\partial \eta }-\frac{\partial z}{\partial \eta }\frac{\partial x}{\partial \xi })}^2},$

式中，x,y和z为空间整体坐标。

连续介质有限元离散模型的外荷载等效结点力^fF_i^l，是单元外荷载等效结点力的叠加。

在步骤(2)中，连续介质有限元离散模型的介质应力等效结点力设为^sF_i^l，其中，下标 i和上标l的意义同前面；上标s表示介质应力。

有限单元法中，基于虚功原理或加权余量法，给出单元应力对应的单元等效结点力计算 公式。对于三维等参单元，计算公式为：

公式中，[B]是单元应变矩阵；其余各量的意义同前面。

连续介质有限元模型的介质应力等效结点力^sF_i^l，是单元应力对应的单元等效结点力的 叠加。

介质有限元模型的结点不平衡力F_i^l，是外荷载等效结点力^fF_i^l和介质应力等效结点力 ^sF_i^l的合力：

F_i^l＝^fF_i^l+^sF_i^l(4)

在步骤(2)中，介质有限元模型中结点的局部刚度系数，为围绕该结点的所有单元的结 点刚度系数之和。

有限单元体系结点局部刚度系数定义为：

$K_i^l=\underset{N}{\Sigma }{}^{e}k_i^l---\left(5\right)$

式中，为结点l在单元e中在i方向的结点单元刚度系数，其为单元刚度矩阵主对角线 上的元素；N为围绕结点l的单元数。对于几何非线性问题，有限元模型结点局部刚度系数需要在每运动步后更新。

在步骤(4)中，设有限元模型结点位移为

有限元模型结点位移计算公式为：

$u_i^l=r\times \frac{F_i^l}{K_i^l}---\left(6\right)$

公式中，r为一个折减系数。计算稳定性要求r的取值小于等于1.0，一般在0.1到1.0 之间。过大的r值，意味着一个运动步中过大的结点位移，将导致介质振荡并逐渐反射放大， 最大结点不平衡力将快速增加，最终导致计算机数据溢出和计算中断。

弹塑性力学中的“试探-修正”算法可用于弹塑性介质求解。在一个运动步中，先计算介质 应力的弹性试探解，然后根据塑性屈服准则判断。如未发生塑性屈服，则介质应力的弹性试 探解即为正确的介质应力解；如果发生塑性屈服，则修正弹性试探解得到正确的介质应力。

基于本发明的该方法，下面给出工程算例：

(1)小变形、弹塑性问题：半无限大介质表面受矩形均布荷载作用。

如图2所示，竖向均布压力为6.5×10⁵Pa，矩形作用区域为3m×6m；取介质计算范围 87m×90m×42m。建立的有限单元网格包含50880个八结点六面体单元和54825个结点。介 质力学指标取值见表1。采用Drucker-Prager理想塑性屈服准则。

表1半无限大介质力学指标

在矩形荷载作用面下方的塑性屈服区，和离荷载作用面相对较远的弹性区，选取若干结 点和单元，这些结点位移和单元应力计算结果，列于表2和3。有限单元法的计算结果也列 于表中。本发明的方法和有限单元法的结果很接近。

表2半无限大介质不同位置若干结点的位移(mm)

表中，有下划线的数值为有限单元法(FEM)的计算结果。

表3半无限大介质不同位置若干单元的应力(10⁵Pa)

该表单元应力为单元一个高斯积分点处的应力，有下划线的数值为有限单元法(FEM) 的计算结果。

本算例中，当介质的强度指标凝聚力减小为0.01MPa，内摩擦角减小为25⁰(膨胀角也相 应取为25⁰)，本发明的方法给出了更大规模的介质塑性流动。发生在矩形荷载作用面中心点 的最大沉降量为103.9166mm。该种情况有限单元法计算遇到收敛困难，难以给出高精度的 结果。

(2)小变形、弹塑性问题：岩体中地下洞室开挖。

洞室尺寸如图3所示。取介质竖向边界离洞室侧边120m，顶部边界离洞顶点100m，底 部边界离洞底边100m，介质计算范围如图3所示。建立的有限单元网格包含6596个六面体 单元和13438个结点。作为平面变形问题，介质厚度取2m。初始地应力取自重应力γh为竖 向分量，0.7γh为横向和厚度方向分量，这里γ为岩体容重，h为某点在介质表面以下的深度。 岩体力学指标取值见表4。采用理想塑性Mohr-Coulomb屈服准则，并考虑拉伸屈服。

表4岩体力学指标

计算了岩体凝聚力0.8MPa，1.0MPa和1.3MPa三种情况。计算给出了洞室开挖引起的围 岩位移，围岩塑性屈服情况和围岩应力。三种情况围岩位移矢量形式相似。岩体凝聚力0.8MPa 情况的围岩位移矢量的情况，岩体最大位移发生在洞室底边中点，方向向上，量值4.405mm。 岩体凝聚力1.0MPa和1.3MPa情况，最大岩体位移分别为4.395mm和4.394mm，发生的位置 和方向与岩体凝聚力0.8MPa情况相同。

计算表明，岩体的凝聚力在1.3MPa以上时，该规模的开挖理论上不产生围岩塑性屈服破 坏；凝聚力0.8MPa和1.0MPa情况发生不同程度的岩体塑性破坏。

(3)大变形、弹塑性问题：金属圆柱体压缩。

如图4所示，一个高度1.2m直径1.0m的金属圆柱体置于刚性垫块上，顶端受到刚性加载 块的压缩。柱体材料视为理想塑性，采用Von-Mises屈服准则。假定压缩过程中柱体顶面和 底面与加载块和垫块之间没有相对位移。不考虑随变形介质密度的变化。

建立的有限单元网格如图5所示，包含1464个六面体单元和1807个结点。高温下材料 力学指标取值见表5。

表5金属圆柱体力学指标

计算给出了准静态荷载作用下圆柱体的变形和应力发展过程。图6a)-图6e)显示了准 静态荷载作用下金属圆柱体一个竖向径向剖面在不同时刻的变形。图6a)-图6e)的圆柱体 压缩量分别为26.2772cm、38.5205cm、50.2398cm、54.8277cm和58.4678cm。

(4)大变形、弹塑性问题：地表空间开挖。

如图7所示，在地表开挖一空间，深度20m，边坡坡度1:1.5，矩形底面10m×20m。取 介质的四个竖向边界离其相应的空间顶边70m，底部边界离空间底面50m，由此介质的计算 范围为220m×210m×70m。建立的有限单元网格包含58600个六面体单元和63321个结点。 初始地应力取自重应力γh为竖向分量，0.55γh为两个水平坐标轴方向的分量，这里γ为土体 容重，h为某点在介质表面以下的深度。

土体力学指标取值见表6。采用理想塑性Mohr-Coulomb屈服准则，并考虑拉伸屈服。不 考虑随变形介质密度的变化。

表6土体力学指标

计算表明，空间开挖后，伴随着较大变形，土体达到平衡状态。最大结点合位移达到 1.917m。图8显示了地表空间开挖后竖向剖面MN的结点位移矢量，其中最大位移矢量为 1.917m。本算例中，当土体的凝聚力取零值、其他指标保持不变时，计算表明，空间开挖后 土体进入塑性流动状态，并从坡脚处涌入空间内。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限 制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付 出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。