岩石、混凝土类材料断裂破坏有限元数值模拟中的有效方法

Ｘ
（ｂ） 图５
（ｃ）
被分解后的几何体
Ｆｉｇ．６
图６
被分解后的几何体的扫描网格图
ｏｆ
Ｆｉｇ．５
Ａｆｔｅｒ ｄｅｅｏｍｐｏｕｎｄｅｄ ｏｆ ｇｅｏｍｅｔｒｙ
Ａｆｔｅｒ ｄｅｃｏｍｐｏｕｎｄｅｄ
ｇｅｏｍｅｔｒｙ
ｓｃａｎ
ｍｅｓｈ
万 方数据
第２期
周尚志，等：岩石、混凝土类材料断裂破坏有限元数值模拟中的有效方法 ｚ＝≈＋（ｚｉ一≈）“＋（￡Ｉ［孙＋（ｚｌ一‰）“一
第２５卷第２期 ２００８年４月
计算力学学报
Ｃｈｉｎｅｓｅ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｍｅｃｈａｎｉｃｓ
Ｖ０１．２５，Ｎｏ．２
Ａｐｒｉｌ ２００８
文章编号：１００７—４７０８（２００８）０２—０２４８—０６
岩石、混凝土类材料断裂破坏有限元 数值模拟中的有效方法
周尚志“， 党发宁２， 陈厚群２’３，
（３）
式中Ｇ为用户定义的网格畸变的判据。上述判据 同样适用于四面体单元。 决定网格重划的另一个因素是裂缝间的单元 干涉。随着变形的发展，接触区域的单元形状和介 质单元间的接触刚度也发生变化，导致界面间单元 干涉现象的发生。这种现象严重影响计算精度，有 时甚至使计算发散。 对于二维四边形等参元，设某一单元的边进入 了界面中，或对于三维六面体等参元，设某单元的
龙巧玲１
（１．长沙理工大学桥梁与结构学院，长沙４１００７６；２．西安理工大学水利水电学院，西安７１００４８， ３．中国水利水电科学研究院工程抗震中心，北京１０００４４）
摘
要：岩石、混凝土类材料断裂破坏有限元数值模拟中的网格重划，依据单元畸变和裂缝介质问的单元干涉作
为网格重划判据，采用几何体重构技术把几何实体分解成能在ＡＮＳＹＳ上实现六面体网格划分的几个部分，利 用体积判断法确定新结点在旧单元的单元编号，在场量传递上采用基于解析性质的等参有限元逆变换。把旧网 格场量信息传递到新网格中。本文对ＡＮＳＹＳ进行二次开发，实现了三维网格重划，网格重划采用单元畸变和 界面干涉两个判据。在网榕再划分前进行几何体重构，提取变形后的点线面信息重新生成实体．充分利用ＡＮ— ＳＹＳ的函数和体积判断法找到新结点在旧网格中的位置。在新旧网格问的场量传递中采用基于解析逆等参单元 法。在平台上实现了三维有限元网格重划技术，最后利用方料的单轴压缩断裂模拟计算检验了传递前后等效塑 性应变分布用载荷信息的交化。证明了所开发系统的正确性。
整体坐标系下的单元棱边与局部坐标系下的单元 棱边一一对应，且棱边的等比率点亦一一对应，对 于直纹面，整体坐标下棱边的等比率点的连线为直 线。整体坐标中点Ｐ（ｘ，，Ｙ，，都）是３个“正交”的直 纹面的交点，如图７所示，等比率点ｚ；，Ｙ；，≈（ｉ＝ ９，ｌｏ，…，２０）可由结点为，为，≈（．『＝１，２，…，８）表 示；ｑ，ＹＪ，乃（歹＝１，２，…，８）已知。以ｚ坐标为例：
（Ｅ－ｍａｉｌｌｚｈｏｕｓｈａｎｇｚｈｉ一１２３＠１６３．ｃｏｒｎ），
党发宁（１９６２一），男。教授．博士生导师； 陈厚群（１９３２一）．男，中国工程院院士．
万 方数据
第２期
周尚志，等：岩石、混凝土类材料断裂破坏有限元数值模拟中的有效方法
２４９
图１
Ｆｉｇ．１
网格重划单元畸形的判断标准
（ａ）四边形等参元
（ａ）几何体
图３ Ｆｉｇ．３
匍国，圊
（ｂ）网格图 （ａ）几何体 （”网格图 受载荷前的几何体和网格图
ａｎｄ
Hale Waihona Puke 图４ Ｆｉｇ．４变形后的几何体和网格图
ｇｅｏｍｅｔｒｙ
Ｂｅｆｏｒｅ ｌｏａｄｉｎｇ ｏｆ ｇｅｏｍｅｔｒｙ
ｍｅｓｈ ｓｋｅｔｃｈ
Ａｆｔｅｒ ｄｉｓｔｏｒｔｉｏｎ ｏｆ
ａｎｄ ｍｅｓｈ ｓｋｅｔｃｈ
图２ Ｆｉｇ．２
（ｂ）六面体等参兀 边界干涉判据示意图
ｃｒｉｔｅｒｉａ
Ｒｅｍｅｓｈｉｎｇ ｅｌｅｍｅｎｔ ａｂｅｒｒａｔｉｏｎ ｏｆ ｃｒｉｔｅｒｉａ
Ｆｏｕｌ ｌｉｎｅ ｉｎｔｅｒｆｅｒｅｎｃｅ ｏｆ
ｓｋｅｔｃｈ
单元畸变的判据，主要体现在单元内角的变化 和单元边长问的相互关系。这主要是由单元的雅 可比变抽象矩阵的值来决定，因为等参单元的坐标 变换、积分变换以及单元体积的计算等都与雅可比 值接有关，要求其值大于零。 在二维问题分析中，对于四结点四边形单元， 当结点的内角小于０。或超过１８０。时，雅可比变换 矩阵的行列式值为负，使计算无法进行下去；当结 点的内角接近０。或１８０。时，使计算精度降低，甚至 引起不收敛。设Ｐ。为所考虑的结点，如图１（ａ）所 示，分别为从Ｐ。到Ｐ：和Ｐ。到Ｐ。的单位矢量，若 由单位矢量Ａ：和．＝【。构成的平行四边形的面积为负 则认为网格畸变。因此对于二维单元网格畸变的判 据如下： ．：１２×Ａ３≤Ｇ（１） 式中Ｇ为用户定义的网格畸变的判据。 对于三维问题，设Ｐ。为结点，如图１（ｂ）所示， Ａ：和Ａ。分别为从尸ｌ到Ｐ２，Ｐ。到Ｐ。和Ｐ，到Ｐ。单 位矢量，若由单位矢量Ａ：，Ａ。和．；【。构成的平行四边 形的体积为负则认为网格畸变。因此对于三维单元 网格畸变的判据如下：
万 方数据
计
算
力
学
学
报
第２５卷
体可以把它分解为几个可以生成六面体网格的规 则几何体，这样整个物体都可以划分成六面体。 由于ＡＮＳＹＳ是采用Ｌａｇｒａｎｇｅ法描述的有限 元，结点和物体的质点是一致的，结点的位置随着 变形的进行而发生改变，故由结点表示的特征点和 特征线基本能描述物体变形后的几何形状，因此可 以通过表面结点来表示变形后的物体。对于一个原 来可用六面体表示的物体，在经过变形之后，可能 变成不规则，这时需要进行几何重构，首先根据旧 网格模型建立一个参考构形，将其外表面结点作为 描述变形后空间实体的特征点，将其表示的变形体 表面定义为特征面，特征面的交线定义为特征线。 对于六面体网格再划分的几何重构过程可以 分以下步骤［３］： （１）将需要网格重划的旧网格的表面结点全 部提取出来，作为控制新网格结点所在空间位置的 控制点。 （２）将参考构形中特征面的交线定义为特征 线，并将其上的单元结点提取出来，作为建造重构 实体边线所需要的型值点和插值点。 （３）将旧网格中与模具接触的单元结点提取 出来，作为变形体几何重构的控制点和型值点。 （４）根据特征线上的型值点、控制点插值点建
性应变量高出最大值时，首先通过减少加载步长的 方法来减小最大塑性应变增量。如果当加载步长减 小到一定程度后，最大塑性应变增量仍高于给定 值，即认为网格单元需要重划。 ３
网格再划分几何体重构
在ＡＮＳＹＳ中构造六面体的方法有映射法和
扫描法，两者都要有较规则的几何体。如果该物体 不规则，则不能生成六面体网格。对于不规则的物
（Ａ２×Ａ３）Ａ。≤≤Ｃ＿
（２）
面进入了界面，如图２所示，当干涉的大小达到一 定程度时，将导致不精确解。设Ｐ点为干涉边或面 的中点，Ｑ为该单元干涉边的对边中点或干涉面的 对面中点。Ｒ点是ＱＰ连线与界面表面的交点，是 Ｐ点与Ｑ点之间的距离。干涉判断准则为 ｈ；／ｈ。≥Ｑ 式中Ｑ为干涉判据常数。 模拟精度也是考虑网格重划的一个因素［３］，其 解决方案是在每一个施加的增量载荷步内，采用逐 点比较的方法，即对每一个高斯点上所产生的弹塑 性应变增量和应力进行误差的检测，检测的标准是 真实弹塑性应变增量和应力进行误差的检测。一般 而言，真实的应变是很难确定的，因此常常利用光 滑后的应力和应变增量来代替真实的应力和应变 增量。由于这种算法在计算过程中需要逐点判断， 虽然有利于控制有限元的计算精度，但是计算量非 常大，效率较低。对于单元数量较多的问题，一般无 法保证模拟效率。 文献［３３在计算精度的控制方面选用的方法： 在采用力的收敛准则的同时，增加了一个最大弹塑 性应变增量的判据，即在每一个加载步时，要求在 每一个高斯点上得到的最大塑性应变增量的应在 给定的范围之内，一般可取０．０５左右；当最大弹塑
关键词：ＡＮＳＹＳ；非均质材料；断裂；数值模拟；网格重划 中图分类号：ＴＵ３７５；０２４２．２１ 文献标识码：Ａ
１
引
言
止计算，重新划分适合于计算的网格，然后将旧网 格的场变量信息传到新网格上，接着进行进一步的 有限元分析计算。为了更好地完成新旧网格间的 信息传递，人们利用距离插值、面积插值、体各插值 用等参逆变换等方法进行场变量的传递。对于三 维弹塑性有限元，采用等参逆变换法是比较合适 的，这主要是由于六面体单元场变量复杂性和保证 高精度场变量传递所决定的。 本文依据单元畸变和界面间的单元干涉作为 网格重划判据，采用几何体重构技术进行六面体网 格再划分，利用体积判断法确定机关报结点在旧单 元的编号，在场量传递上采用基于解析性质的等参 有限元逆变换，把旧网格场量信息传递到新网格 中。在ＡＮＳＹＳ平台上实现了网格重划技术，最后 利用方料的单轴压缩断裂模拟计算检验了传递前 后等效塑性应变分布用载荷信息的变化，证明了所 开发系统的正确性。
岩石、混凝土类材料内部存在着大量的细、微 观缺陷，如骨料界面、孔洞和微裂隙等，是一种非均 匀的材料。随着试验技术的发展和现代数值分析 技术及高速计算机技术的发展，已经可以对岩石、 混凝土内部结构进行非线性有限元全过程分析，来 预测在加载直至破坏整个过程的应力、应变及最后 的破坏模式。但是，人们发现，计算分析的结果往 往不容易与试验数据符合得很好。其原因之一，是 在对岩石、混凝土试件进行加载作用发生破坏后， 裂缝界面之间不同介质相互干涉，初始过程的有限 元网格将会变形，在变形过程中，网格逐渐畸变，若 把已经畸变的网格形状作为增量计算的参考状态， 将导致计算精度降低，甚至引起不收敛，而不能继 续进行计算。因此对于涉及复杂变形的过程，难以 用一成不变的单元网格把变形过程模拟到底。为 克服上述问题，当网格变形到一定程度后，必须停
２５１
４
新旧网格间的信息传递
对于实现新旧网格间的信息传递首先要确定
≈一（ｚｉ＋麓）ｕ３
（４）