基于ABAQUS-Python无限元的动力人工边界研究

党康宁，苏晨辉，肖 瑜，张静宜

(1.陕西省引汉济渭工程建设有限公司，陕西 西安 710010；2. 成都市龙泉驿区水务局，四川 成都 610100)

工程结构体系动力响应分析中，为考虑无限地基辐射阻尼效应，各类人工边界被广泛研究和应用。其中，粘弹性人工边界得到了学者们的广泛采用和认可，但在有限元模型建立时，需要在粘弹性人工边界上建立大量三向弹簧阻尼器，工作量很大。其参数的设置和集中等效荷载时程的添加难度很大。但是，近年来学者们采用无限元技术与有限元结合，开启了相关研究工作。无限元由Ungless[1]提出，Lysmer[2]、Zienkiewicz[3]、Bettess[4]、Beer等[5]进行了发展和改进，国内葛修润[6]，张楚汉等[7]也较早进行了这方面研究。

近年来，一些学者[8-12]利用ABAQUS有限元软件中的无限单元，通过有限元—无限元(FE-IE)方式研究了水电站厂房、重力坝、拱桥、地下隧道等结构反应，取得较好结果。实现过程或采用编制INP输入文件，或采用FORTRAN单独编写程序实现地震波动荷载的生成，此过程十分繁琐并容易出错，增加了无限元动力人工边界的实现难度。

Python作为目前最流行的脚本语言之一，具有简洁、跨平台等优点，是ABAQUS程序前后处理层次的接口语言，可获取ABAQUS模型中信息，并通过循环语句、内置函数等能实现模型、荷载的高效操纵，从而便捷解决无限元动力人工边界的繁琐问题。

本文基于波动理论，利用ABAQUS软件的无限单元，建立无限元—有限元模型，通过Python编程二次开发，实现各侧人工边界上节点的刚度参数施加、波动荷载生成以及荷载的自动施加，极大减小了前处理工作量，最后，通过算例验证了提出方法的精确性。

1 基于波动理论的无限元动力边界

工程结构地震动力响应分析模拟时，地震作用作为一种外源输入，在截取的地基范围内存在地震波入射和反射，基于弹性介质的波动理论，结合边界节点等效荷载的可叠加原理，可解决外源波的入射问题。在有限域人工边界通过引入一种特殊的有限元—无限元，该单元与有限元无缝衔接，并通过几何映射，在局部坐标中构造插值形状函数，实现计算范围趋于无限远。无限元边界法作为有限元方法的补充，具有良好的“协调性”，其对复杂散射波动的控制能力更优[11]。

目前，无限元边界实现较多是无限元—有限元相结合的方式，即在工程结构的地基有限元区域边界外通过一层无限元连接，并通过一定的衰减函数实现能量的吸收。无限元静动力分析理论分别基于Zienkiewicz[13]和Lysmer[2]等进行研究。

1.1 无限元动力边界原理

无限单元可以充当吸收边界，单元设置了阻尼矩阵，根据荷载情况自动计算值大小。其动力分析吸收散射波原理如下。

地震压缩波(P波)入射时，其在均质无限弹性体中运动方程为：

(1)

压缩波沿x轴负向传播，从有限区域进入无限元边界，则边界位移应为：

(2)

同时，在入射位移波经过边界节点后反射的位移波如下：

(3)

入射和反射位移总和：

ux=f1(x-cpt)+f2(x+cpt)

(4)

入射和反射速度总和：

(5)

根据弹性力学可知：

(6)

(7)

此时，边界节点上阻尼应力为：

(8)

式中CBN为阻尼器参数。

经过人工边界后，散射波产生的应力应与阻尼应力相等，从而消除散射波影响，即σx=σdamp，所以：

(9)

由上式计算可得：

(10)

(10)

同理可得剪切波的无限单元内嵌阻尼器系数CBT：

CBT=ρcs

(11)

式中cs为压缩波波速。

1.2 无限元边界上等效节点力求解

当外部地震荷载传播进入有限元区域后，地基的运动由入射波和反射波叠加组成。散射回有限元边界的波由无限元吸收，同时体现地基的弹性作用。假定边界区域弹性小变形，可将荷载转化为边界上节点等效应力，从而解决外源波入射问题。

对于粘弹性人工边界节点力公式为：

(12)

无限元边界中已自动嵌入刚度项，因此，取公式中弹簧刚度为0，得到边界地震动输入的等效节点力表达式。对于从模型底边界垂直入射的波，依据一维波动理论可分别求得边界各节点的等效波动荷载，将粘弹性人工边界等效节点力公式中刚度项去掉即可。

1.3 耦合模型的建立

目前，无限元人工边界一般配合有限元模型完成，在有限元地基外部设置1层无限单元。由于阻尼项已被无限单元考虑，因此，只需添加刚度参数。对于二维和三维问题，无限单元和有限单元的节点连接形式如图1～2所示。

a 二维单元

图2 二维有限元-无限元半空间自由场波动计算模型示意

2 Python脚本实现无限元边界荷载生成及施加

前节推导了等效荷载的求解过程，接下来需要计算各节点荷载时程，并在模型上施加。

作为ABAQUS的内核语言，采用Python脚本语言进行二次开发具有天然优势，通过增加特有对象模型，Python能够直接与ABAQUS模型交互数据，读取模型信息，并更改模型设置，为实现有限元、无限元交界处等效荷载的生成和施加奠定了基础，能有效减少前处理工作量。

2.1 ABAQUS中Python脚本概述

Python脚本语言是面向对象语言，具有对象(object)、成员(member)、方法(method)、构造函数(constructor)、类(class)、模块(module)和字典(dictionary)等基本特征。在ABAQUS中Python还有数据库(database)、容器(Repository)、声明使用(Access)及路径(Path)等特有性质[14]。

其中数据库负责存储模型的各种信息，是具有ABAQUS特征的一类特殊的对象，例如，本文主要对模型数据库进行操作，就是mdb。mdb对象是存放ABAQUS有限元模型的根对象(见图3)，其中Models是仓库类型，包含有parts、rootAssembly、loads和steps等成员对象，各成员对象由含有许多下级成员对象。

图3 ABAQUS中的mdb对象层次示意

在实现本文方法时，需要在Python中声明导入一些基本模块，以便使用ABAQUS中各对象。

from abaqus import * #导入ABAQUS模块所有公共对象abaqus。

from abaqusConstants import * #导入所有符号常量abaqusConstants。

from caeModules import * #导入caeModules窗口，实现ABAQUS窗口中所有对象模块的导入。

此外，本文还需要载入io模块，用于文件操作；载入mesh模块，用于对模型网格操作。

2.2 无限元边界荷载生成和施加实现的主要流程

模型通过GUI建立，并准备好部分前处理文件，然后通过Python编程实现无限元边界荷载生成和施加，主要流程为：

1)采用前述方法进行无限元-有限元模型的建立，并划分好网格。通过计算节点反力得到模型四侧及底面上各节点的影响面积，并分别写入文件。准备好要输入的三向荷载的速度时程。

2) 编写Python脚本，在上述初始化后，通过变量定义模型名称、装配件名称、分析步名称、时间间隔、材料参数、模型长宽高等。

3) 利用文件操作函数，并使用循环将各方向速度时程、节点影响面积等文件读入字典当中备用。

4)获得人工边界上各节点的坐标值，并进行分组处理，以此来判断节点所在模型哪个侧面。

5)对每组中节点通过前述公式得到模型节点力及波动时程荷载。定义模型荷载函数，并按方向施加节点荷载。

程序实现流程如图4所示。

图4 荷载生成和施加程序流程示意

3 算例验证

为验证本文所讨论的无限元边界法准确性及所编制程序正确性，采用文献[15]的算例模型进行验证。材料弹性模量为24 MPa，剪切模量为100 MPa，泊松比为0.2，密度为1 000 kg/m3。因此，该材料剪切波速为100 m/s。模型底端作用的荷载为速度脉冲，其表达式为：

(13)

其中f=4.0，0≤t≤0.25。

建立有限元-无限元耦合模型如图5所示，XY平面为水平向，Z坐标轴指向为模型竖向。有限元区域见图5a，模型尺寸为6 m×6 m×50 m(长、宽、高)，采用8节点三维实体单元离散，网格尺寸为1 m。在有限元模型的底部及四周包裹1层无限单元(见图5b)，从而完成有限元-无限元耦合模型的建立(见图5c)。取有限元模型沿Z轴中轴线上底部、中部和顶部3个点为位移监测点，节点号分别为25，1 250，2 475。

a 有限元区域

图6～7分别给出了模型沿高度方向底部、中部和顶部的竖向、水平位移。

图6 模型监测点竖向位移示意

图7 模型监测点水平向位移示意

由图6～7可知，无限元边界结果与理论值十分接近，好于粘弹性边界的计算值，由此验证了本文提出的建模、荷载生成和施加方法的正确性和精确性，说明基于波动理论的无限元边界能够很好地解决外源输入时的地基辐射阻尼问题。

4 结语

本文进行了基于ABAQUS-Python无限元的动力人工边界研究。基于波在弹性均匀介质中传播理论，推导了在无限元边界上各节点荷载，将加速度、位移荷载转化为等效应力，并进一步得到有限元模型中集中荷载时程，采用ABAQUS内嵌的原生脚本语言进行二次编程开发，实现了计算模型人工边界节点上荷载的快速生成和准确施加。通过小算例比较了粘弹性人工边界、本文无限元边界和理论值，结果表明本文无限元-有限元模型的无限元人工边界实现方法具有很高的精度，且二次开发的程序具有代码量少、便于迁移应用等特点。