“另类单元生死”在大坝等结构施工仿真中的应用

王 旭，赖国伟，王志鹏

（武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉 430072）

0 引言

采用隐式有限元方法求解复杂工程问题已经十分普及，但只适用于小位移、小转动、小变形的界限内。隐式有限元方法会经常因局部材料破坏导致迭代不收敛，所以一旦需要模拟施工过程并且还要进行高度非线性分析时，隐式算法会产生严重的收敛困难问题。显式有限元是求解高度非线性问题的一个非常有效的工具，在工程界和学术界已经取得高度的认可。Hooputra， H.等[1]利用显式准静态方法分析了薄壁铝材逐步失效破坏，并与实验数据进行了对比验证；Linde P等[2]模拟纤维金属成型得到了理想的结果；罗州[3]采用显式有限元分析程序ABAQUS/Explicit提供的动力显式算法，并加入与实际工程更相似的约束条件，从而得到了更为准确的结果。在混凝土工程方面，王素裹等[4]为了考察ABAQUS/Explicit显式分析方法在钢筋混凝土结构上应用的适用性和准确性，用显式与隐式两种分析方法对迭代算法和求解时间进行了对比，计算结果表明，显式分析方法能够在解决材料失效和破坏导致的收敛问题的基础上，较为准确地模拟钢筋混凝土结构响应。

ABAQUS/Explicit显式有限元方法最显著的特点是不需要在隐式方法中用到的整体刚度矩阵，显式地前推模型的状态，不需要迭代和收敛准则；而且在相同的模拟条件下，显式方法比隐式方法需要更少的磁盘空间和内存。但是ABAQUS/Explicit显式有限元没有类似于隐式的 “单元生死”，所以无法模拟施工过程。如果不模拟施工过程，仿真模拟得到的结果会与现实的结果差别很大。

为了使用ABAQUS/Explicit显式有限元强大的求解器来解决工程项目的诸多高度非线性问题并模拟施工过程，本文提出通过数据传递的方式来模拟单元生死，并将此方式定义为 “另类单元生死”法。基于该方法，模拟了小湾拱坝的施工过程，并将计算结果与Ansys隐式计算结果进行对比。

1 “另类单元生死”和显式算法

1.1 “另类单元生死”

本文自定义的 “另类单元生死”原理是[5]，通过ABAQUS/Explicit模块中有限元模型和应力变形场等数据传递的方法，亦即将结构前一荷载步的有限元模型和计算应力变形场等数据作为下一荷载步计算的初始输入数据，来解决ABAQUS/Explicit显式有限元计算的单元生死问题，达到模拟结构施工过程的目的。

“另类单元生死”法与隐式有限元的 “单元生死”法不同。隐式有限元的 “单元生死”是根据施工步的顺序而改变相应单元刚度并且设定与单元相联系的质量、阻尼、比热及其他类似的效果的值。而 “另类单元生死”则是保留当前施工步仅有的有限元模型，删除后面所有施工步出现的单元和结点，然后对当前施工步进行求解；之后，将当前施工步的有限元模型和计算结果传递到下一个施工步，作为下一个施工步定义初始状态场，下一个施工步即在当前施工步的有限元模型和计算结果基础上，增加相应的单元和结点，接着进行求解。

1.2 “另类单元生死”的实现方法

基于ABAQUS/Explicit显式有限元，用 “另类单元生死”模拟施工步，把当前施工步的有限元模型和计算结果传递到下一个施工步的ABAQUS命令流文件，亦即inp文件格式如下所示：

*IMPORT，STEP=N，UPDATE=no：把当前施工步的计算结果传递到下一个施工步，其中，N表示当前施工步为第几个施工步；

*IMPORT NSET：把当前施工步的结点组件传递到下一个施工步；

*IMPORT ELSET：把当前施工步的单元组件传递到下一个施工步。

采用 “另类单元生死”法需注意的是，ABAQUS计算所产生的后缀名为abq、stt、pac、prt和odb的文件必须保存在同一工作目录下，才能保证计算顺利进行。

1.3 显式算法

ABAQUS/Explicit应用中心差分方法对运动方程进行显式的时间积分[6]，由一个增量步的动力学条件计算下一个增量步的动力学条件。计算公式如下：

1)哪一种图形密码最安全；2)图形密码与传统的六位数字密码，哪一个更安全；3)改变构造图案的规则，密码的安全程度会有哪些改变。

该算法不需要迭代和收敛准则，为了得到精确的结果，时间增量步必须相当小，这样在增量步中加速度几乎为常数。由于时间增量步必须很小，一个典型的分析需要成千上万个增量步，但在求解过程中，不必同时求解联立方程组，故每一个增量步的计算成本并不高。

2 工程实例

2.1 小湾拱坝及计算模型

小湾大坝为混凝土双曲拱坝，坝顶高程1 245 m，坝基面高程953 m，最大坝高292 m，坝顶宽12 m，最大坝底宽72.912 m，坝顶弧长901.771 m，大坝混凝土量842万m3，弧高比3.056，弦高比2.74，厚高比0.25。拱坝承受的总水推力为1.66×108kN。有限元网格划分见图1。其中，坝体及岩体单元以六面体单元为主，四面体单元为辅。整个计算模型有48 609个单元，其中坝体有25 136个单元。模型整体坐标系X正向指向右岸，Y正向指向顺水流向，Z正向为铅直向上。

图1 小湾拱坝计算模型

为检验 “另类单元生死”法的正确性，ABAQUS/Explicit显式有限元和Ansys隐式有限元采用同一有限元模型。

2.2 基本资料

（1）水位。上游正常蓄水位1 240.0 m；下游水位1 004.0 m。

（2）泥沙。淤沙高程1 097 m，淤沙容重9.5 kN/m3， 内摩擦角 24°。

（3）依据混凝土设计参数为：容重γ=24 kN/m3，泊松比 γc=0．2， 线膨胀系数 α=8．26×10-6/℃， 抗剪强度 f′=1．4， c′=1．6 MPa。 不同混凝土分区 （A 区、 A0区， B 区， C 区） 180 d弹性模量分别为 23．1、 22．4、21．7 GPa，已按工程经验考虑了混凝土的徐变影响。

（4）坝基采用均质模型，弹性模量21 GPa，泊松比 0．26， 基础岩体抗剪强度 f′=1．4， c′=1．8 MPa。

2.3 材料本构关系

材料本构关系采用理想弹塑性模型[7，8]，屈服准则采用常用的Drucker-Prager准则，即

式中，I1为应力张量的第一不变量；J2为应力偏张量的第二不变量；α、k均为材料常数，与材料粘结力c、内摩擦角φ的关系取决于Drucker-Prager圆锥面与摩尔库仑六棱锥面之间的相互关系。本文中α、k采用外接圆锥公式计算

2.4 荷载及工况

计算中考虑的荷载有拱坝自重、上下游坝面水压力、上游泥沙压力、运行期温度荷载、水库库盆水压力。自重和上游坝面水压力按预备工况分期施加 （见表1）；下游坝面水压力、上游泥沙压力、运行期温度荷载、水库库盆水压力均一次性施加，考虑由面力代替模拟渗流场的体积力，乘以折减系数，折减系数取 0．8。

表1 施工过程模拟

3 计算结果及分析

考虑到Ansys隐式有限元计算[9]通过隐式的 “单元生死”法能够模拟施工过程，为验证 “另类单元生死”法的正确性，本文将分别采用ABAQUS/Explicit显式有限元与Ansys隐式有限元进行对比计算分析，以拉应力为正。顺河向位移等值线见图2，第一、第三主应力等值线见图3。坝体最大位移和最大应力及其位置见表2。

图2 顺河向位移等值线 （单位：m）

表2 坝体最大位移和最大应力及其位置

从图2、3和表2可以看出：①顺河向位移分布大致相同，而且都呈左右对称分布，用两种软件计算的顺河向最大位移最大值差别不大，并且出现的位置基本一致；②第一主应力在坝顶左右有少许不同，其中用ABAQUS/Explicit计算的数值稍大，但是整体的第一主应力分布基本相似，且都呈左右对称分布，用两种软件计算的第一主应力最大值差别不大，并且出现的位置基本一致；③第三主应力在坝体1/4坝高处有少许差别，其中用ABAQUS/Explicit计算的数值稍偏大。虽然两者局部有所不同，但是从整体上来看，第三主应力分布大致类似，且都呈左右对称，除此之外，用两种软件计算的第三主应力最大值差别不大，并且出现的位置基本一致。

综上所述，对比ABAQUS/Explicit显式有限元与Ansys隐式有限元的位移、应力及其分布情况可以看出，基于 “另类单元生死”法的显式有限元计算结果与传统隐式有限元区别很小，表明ABAQUS/Explicit显式有限元采用 “另类单元生死”法，通过数据传递命令能够传递结构施工荷载步的位移、应力等信息，计算结果正确。

4 结论

图3 第一、第三主应力等值线 （单位：MPa）

（1）“另类单元生死” 法的原理是， 通过ABAQUS/Explicit模块中有限元模型和应力变形场等数据传递的方法，亦即将结构前一荷载步的有限元模型和计算应力变形场等数据作为下一荷载步计算的初始输入数据，来解决ABAQUS/Explicit显式有限元计算的单元生死问题，达到模拟结构施工过程的目的。

（2）在塑性破坏的情况下，ABAQUS/Explicit显式有限元计算结果与Ansys隐式计算结果差别很小，表明ABAQUS/Explicit显式有限元通过数据传递方式能够模拟施工过程， “另类单元生死”法可以达到隐式的 “单元生死”要求。

（3）借助 “另类单元生死”法，广大科研人员可以用ABAQUS/Explicit显式有限元强大的非线性求解器对大坝等结构进行高度非线性仿真分析。

.

［1］ Hooputra H,Gese H,Dell H,etal.A Comprehensive Failure Model for Crashworthiness Simulation of Aluminium Extrusions［J］.International Journal of Crashworthiness,2004（9）:449-463.

［2］ Linde P,Pleitner J.Modelling and Simulation of Fiber Metal Laminates［A］.ABAQUS Users'Conference,2004.

［3］ 罗州.金属环件轧制过程刚粘塑性动力显式算法有限元模拟［D］.武汉：武汉理工大学，2004.

［4］ 王素裹，韩小雷，季静.ABAQUS显式分析方法在钢筋混凝土结构中的应用［J］.科学技术与工程， 2009， 9（16）:4688-4692.

［5］ 曹金凤，石亦平.ABAQUS有限元分析常见问题解答［M］.北京：机械工业出版社，2009.

［6］ 庄茁，由小川，等.基于ABAQUS的有限元分析和应用［M］.北京：清华大学出版社，2009.

［7］ 王金昌，陈叶开.ABAQUS在土木工程中的应用［M］.杭州：浙江大学出版社，2006.

［8］ 陈明祥.弹塑性力学［M］.北京：科学出版社，2007.

［9］ 王国强.实用工程数值模拟技术及其在ANSYS上的实践［M］.西安：西北工业大学出版社，1999.