摘要

       ABAQUS有限元软件在水利工程中应用广泛，人工粘弹性边界和等效节点荷载法是水工构筑物地震动力响应分析的常用方法。本文基于ABAQUS计算内核，采用Python计算机语言编写了二次开发程序，实现了人工边界处各节点的弹簧、阻尼器及等效节点荷载自动设置的功能，并具体阐述了程序编译流程和关键要素。建立了典型SV波底面入射模型对程序的适用性进行了验证，结果表明地震波的传递及反射规律与理论解答一致，程序合理可行。人工粘弹性边界通过在人工边界处设置弹簧和阻尼器吸收地震波，进而取消地震波的反射现象，分析结果与实际更为相符。地震动的输入方式会对地震分析的结果产生较大影响，刘晶波等将地震波动问题转换为波源问题，将地震波转变为等效节点荷载，再以集中力的方式施加在人工边界的各节点上[1]，这种地震动输入方式因其较好的准确性广泛运用于水利工程中。

       采用有限元模型分析地震问题时，需进行两部分的计算，第一是人工边界处各节点的弹簧刚度和阻尼，第二是人工边界处各节点的等效节点荷载。实际工程中的有限元模型通常网格节点较多，如采用人工手动添加弹簧、阻尼器及节点的等效荷载，在效率性上是不可取的，采用计算机编码实现批量化的前处理功能可根本上解决上述问题。

       ABAQUS是目前较为常用的大型通用有限元软件，提供了二次开发的接口，软件计算内核可根据用户编写的二次开发程序进行定制化的计算。目前ABAQUS软件实现粘弹性边界设置及等效节点荷载输入的程序通常采用Fortran、Python及Matalab编写 [2~4]，但程序自动化程度仍有待提高，具体存在以下问题：

（1）在后处理模块人工提取节点反力并输出归整，操作繁琐复杂，要求使用者有较高的ABAQUS前后处理经验。

（2）在外部编译器内编写边界节点的弹簧刚度、阻尼器阻尼及等效节点荷载的计算程序并进行相应计算，再以INP文件方式输入到ABAQUS计算内核。使用者需对ABAQUS数据格式有较深入的了解，并且涉及多数据源之间的数据传递，处理复杂模型时容错率低。

（3）当采用Fortran和Matalab进行编译时，ABAQUS软件没有提供可交互的软件调试功能，代码的编写效率较低。

       本文基于前人研究成果，以ABAQUS计算结果Odb数据库为中转存储单元，采用Python计算机语言编写了粘弹性边界和等效节点荷载的自动化设置程序，实现了单一数据源在ABAQUS模型内的相互传递功能，全过程无人工操作，极大地提高了前处理的效率和准确性。

图1 程序总体框架

本文程序涉及较多ABAQUS内置Python库的函数命令[5]，为便于读者理解，后文以缩写形式叙述：

（1）模型复制函数Copy_Model：mdb.Model（name='', objectToCopy=mdb.models['']）；

（2）均布力设置函数Pressure_Set：mdb.models[''].Pressure（）；

（3）边界条件设置函数Boundry_Set：mdb.models[''].DisplacementBC（）；

（4）作业等待函数Wait_Job：mdb.jobs[''].waitForCompletion（）；

（5）节点设置函数Node_Set：odb.rootAssembly.instances[''].nodeSets[]

（6）弹簧阻尼器施加函数Spring_Set：mdb.models[''].rootAssembly.engineeringFeatures.SpringDashpotToGround（）；

（7）集中力施加函数Lord_Set：mdb.models[''].ConcentratedForce（）

（8）数据库打开函数Odb_Open：openOdb（Pathname）

1.3弹簧与阻尼器的设置

弹簧刚度与阻尼器阻尼的取值决定了粘弹性边界的合理性，计算公式为[6]：

（1）

      采用式（1）计算得到边界所有节点的法向与切向弹簧系数和阻尼系数后，乘以该节点处的控制面积即可得到地震分析模型内设置的弹簧刚度系数和阻尼系数。既有研究表明，将模型边界处自由度固定、施加法线方向大小为1Pa的均布荷载时，边界各节点的支反力即为该节点的控制面积[7]。

      现有二次开发程序仍需人工提取边界处的节点反力值，输出为txt、csv等格式的文件，再通过外部编译器读取，计算后以INP文件方式传入ABAQUS计算内核。本文开发程序以初始计算支反力的Odb数据库文件为数据存储中转单元，采用编写的Python程序，以循环遍历的方式直接读取数据库内人工边界处各节点信息，按照式（1）计算后导入后续地震分析模型，实现了单一数据源在ABAQUS软件内部的流通与调用，前处理的自动化程度得到了显著提高，具体步骤如下（流程图见图2）：

（1）首先设置各人工边界处的节点集合、表面集合和几何集合，再采用ABAQUS提供的模型复制函数Copy_Model将现有模型复制。编码时需要注意ABAQUS程序内对集合类型的定义，不同的函数调用的集合类型会存在差异。

（2）采用Boundry_Set函数读取几何集合并设置边界条件，采用Pressure_Set函数读取表面集合并设置均布力。

（3）批量建立作业后，运用作业等待函数Wait_Job依次进行计算。

（4）采用OdbSet类读取预先设置好的边界节点集合，再通过GetSubset函数读取OdbSet集合内各节点的节点编号、坐标以及支反力，并返回到节点信息存储序列Array_Save中。

（5）读取Array_Save序列中节点支反力、节点编号以及坐标，进行弹簧刚度和阻尼器阻尼的计算，再采用Spring_Set函数设置。

图2 弹簧与阻尼器设置流程图

1.4等效节点荷载的设置

       将地震波动问题转换为波源问题时，需依据波场分离理论计算边界处各节点的等效节点荷载，具体推导过程可参考中国地震局马笙杰所撰文章[8]。两侧边界各节点处由入射波和反射波引起的位移与速度响应计算公式为：

（2）

 底部边界处各节点的等效节点荷载计算公式为：

（3）左边界各节点等效节点荷载计算公式为：

（4）右边界各节点等效节点荷载计算公式为：

（5）

        由式（2）~（5）可知，人工边界节点的弹簧刚度和阻尼器阻尼已知的情况下，计算节点等效荷载的关键是模型侧边界节点在地震时程内的位移响应和速度响应（下文统称响应），它们的响应是由底部入射波和顶部自由表面反射波的叠加所引起。而任意节点处对入射波和反射波的响应可表示为底部入射波函数u0（t）在t轴方向大小为△t的相位差u0（t-Δt）。

       当u0（t）函数形式已知时，可直接采用函数计算节点的响应。但入射波的数据通常以离散的数据点方式提供，很难在保证精度的前提下采用一个函数去拟合，本文建议采用拉格朗日插值法[9]去计算节点的响应（程序流程见图3），即在计算边界节点在t=tn时所对应的响应前，采用For循环遍历入射波的时间数据，找寻到ti<tn-△t<ti+1这两个数据点，计算公式如下：

（6）

        式中，tn为时间参数（s）；为入射波或反射波在b节点的延时（s）；ui和vi是与入射波数据中与ti相对应的位移值（m）和速度值（m·s-1）；其余参数同前。

 

图3 插值计算流程图

        读者需注意，当模型边界处节点数量较多或入射波数据量较大时，建议采用二分法定位插值区间，可极大地提高循环的遍历速度。计算好边界各节点入射波和反射波的响应后，即可采用式（3）~（5）计算边界各节点的节点等效荷载，通过节点编号定位后以集中荷载函数Lord_Set施加到地震分析模型（具体程序流程见图4）。

图4 等效节点荷载程序流程图

02 适用性验证

       为验证程序的合理性与正确性，现以文献8工况为例进行二维平面底边界垂直入射SV波的地震响应分析。在二维无限弹性半空间截取 50 m×50 m 的有限区域作为计算区域，材料密度ρ=2000kg·m-3；弹性模量E=200MPa；泊松比v=0.25；剪切波速Cs=200 m/s。采用尺寸为1m的平面应变单元（CPE8）对区域进行离散，取模型顶面中点A和底部入射边界中点B以及相同纵坐标的C、D、E、F为监测点（见图5），入射波位移时程函数采用有限差分近似（见图6）：

式中，T为脉冲持时（s），为Heaviside函数。

 

图5 计算图示

图6 入射波位移时程曲线

        提取入射面各监测点位移与时程关系曲线（见图7），提取顶面自由表面处各监测点位移与时程关系曲线（见图8）；提取不同时刻模型位移云图（见图9）。

        由图7~9可知，地震波在0.25s时传至自由表面并发生反射，自由表面处位移响应是入射波位移的2倍，位移的时程变化曲线形态也与入射波一致。底部边界的位移响应分为入射波和反射波的两次响应，当入射波传入，边界的变形与波形一致，0.2s时位移为零，当0.5s时反射波到达底部边界，底部边界再次发生与波形一致的位移响应，0.7s后位移再次为零，粘弹性边界能较好地吸收反射波，并未再次发生折射。等效节点荷载能够较好地模拟地震动输入，SV波的传递均匀，同一纵坐标下各监测点的位移响应基本一致，

图7 入射面节点位移响应曲线

图8 自由表面节点位移响应曲线

（b）t=0.45s

图9 结构位移响应云图

03 结论

       本文采用Python计算机语言编写了ABAQUS粘弹性边界及地震动输入的二次开发程序，阐述了程序编写的具体流程和关键因素，并对比了二维SV波入射情况下模型特征节点的位移响应与理论解答，得到了以下重要结论：

（1）Python计算机语言与ABAQUS软件有较好的兼容性，编写的程序以Odb结果数据库作为中转存储单元，能够直接访问Odb数据库读取节点关键信息，实现了单一数据源在ABAQUS软件内部的传递，无需人工提取数据并进行数据归整。

（2）等效节点荷载可采用ABAQUS库内的Lord_Set函数直接施加到边界节点，无需在程序外部编写INP文件后导入。

（3）等效节点荷载计算的关键是侧边界各节点对入射波和反射波的位移响应和速度响应，可依据入射边界处的位移和速度时程数据，采用拉格朗日插值函数计算相应节点的响应。

（4）二次开发程序能够准确地实现粘弹性边界设置及等效节点荷载输入的功能，计算结果与理论值基本一致，能够运用于工程实践。

（5）对于考虑结构-土相互作用的地震分析，读者可依据文献[10]中的计算方法，采用本文程序框架编码实现。