二维黏弹性人工边界单元及地震波输入在ANSYS中的实现

尹训强,罗 勇,王桂萱

( 大连大学 土木工程技术研究与开发中心，辽宁 大连 116622 )

对无限域或半无限域进行地震响应分析时，往往采用有限元法。该方法通常采取的处理措施是在无限域或半无限域空间中人为地截取一定范围的有限域，进而可以将无限区域转化为有限区域，然后在有限域的基础上，对其进行离散法处理，进而无限域或半无限域最终被转化为有限个单元。同时为了达到近似处理的效果，通常在截取的有限域的边界上人为的施加约束，但是这一处理的缺点是在处理过程中，不仅没有考虑远场无限地基的辐射阻尼效应，同时半空间无限地基弹性恢复性能也被忽略。由于截取的有限域的范围存在不确定性，随着有限域截取范围的增加，自由度也随之增加。基于此离散范围也更加宽广，进而网格划分的单元数及结点数成倍增加，因此，在地震动的动力计算中工作量也进一步被增加，对计算机的要求也更高，进而，在动力学中非线性问题的求解难度进一步增大。目前，针对该关键性问题处理的诸多措施中，较好的处理方法是在截取的有限域边界处设置局部人工边界。目前在工程中较常用的局部边界主要有黏性边界[1]、透射边界[2]、黏弹性边界[3]等。黏性边界的优点是物理概念清晰，在实际动力计算程序中实现起来相对比较容易，对于1阶动力响应精确度较高。基于此，黏性边界单元被嵌入到了很多大型商用软件中以便于工程的实际应用，如LSDYNA、 FIAC 、ABAQUS[4]等。它的缺点是在施加人工边界过程中只着重考虑了人工边界辐射阻尼效应，忽略了人工边界处半无限空间介质的弹性恢复性能。因此在实际工程应用中容易发生整体漂移的现象，进而使计算结果精度不高。透射边界的优点是对于2阶动力响应计算结果存在较高的精度，其缺点是动力计算分析程序的代码编写较复杂、在实际工程应用中可能存在高频失稳的现象。黏弹性边界不仅能有效的模拟半无限空间介质的辐射阻尼效应，同时还能较好的模拟半无限地基的弹性恢复性能。此外还具有良好的频率稳定性和较高的精度。因此在实际工程中得到了广泛的应用[5-10]。

作为大型通用有限元分析软件之一的ANSYS软件，虽然其自身已经存在着丰富的单元库，但其单元库中仍然存在着动力学求解中所需的单元，例如：黏弹元边界单元、非线性单元、损伤单元，为了解决这一关键性问题，本文基于ANSYS平台，采用FORTRAN语言开发了二维黏弹性边界单元用户子程序VSB_UELf，并结合ANSYS二次开发工具UPFs的功能特性，将其成功地嵌入到ANSYS中。最后，通过典型算例验证了用户自定义二维黏弹性边界单元的正确性与可靠性，其计算结果的精度及稳定性均令人满意。

1 二维黏弹性边界单元

为了提高二维黏弹性边界在土-结构相互作用求解过程中的计算精度，对于刚度系数KB和阻尼系数CB的选取显得尤为重要，为了确定合理的KB和CB，其通常需要与人工边界处的材料参数建立起相应的表达式[11-12]，其表达式如下：

(1)

切向阻尼系数：CBt=ρcs

(2)

(3)

(4)

为了便于对表达式中相关参数的理解，现对上式中的相关参数做详细解释，黏弹性人工边界中的切向参数由αt表示，黏弹性人工边界中的法向参数由αn表示。为了保持黏弹性边界在地震动求解过程中的稳定性，依据αn和αt的推荐取值范围，文中αt、αn的取值分别为0.5,1.0。其中二维黏弹性中的αn、αt的推荐取值详见参考文献[13]。

二维黏弹性边界单元的刚度矩阵为：

(5)

二维黏弹性边界单元的阻尼矩阵为：

(6)

2 用户自定义二维黏弹性边界单元的开发流程

在ANSYS软件中的二次开发工具主要有4个：APDL、UPFs、UIDL、TclTk，其中UPFs是ANSYS二次开发的核心工具，UPFs是基于FORTRAN程序的二次开发，其主要功能是实现对ANSYS功能的扩充，即对ANSYS单元库的扩充。同时ANSYS软件也为用户提供了诸多子程序，其中，UserElem.f子程序不仅可以帮助用户传递自定义单元所需的的必要数据，同时还可以在动力计算过程中自动更新标准ANSYS中的数据库和相关文件。需要说明的是UserElem.f子程序是基于ANSYS单元水平上的接口。基于UserElem.f子程序，用户可以几乎创建各种单元类型而无需直接访问ANSYS数据库和文件。因此,基于二维黏弹性人工边界有关理论,借助于ANSYS二次开发工具中的核心工具UPFs以及接口子程序UserElem的功能特性，用户自定义二维黏弹性人工边界的开发具体流程图见图1。

为了实现将用户自定义的二维黏弹性人工边界单元成功的嵌入到ANSYS中，方便于推广与应用，其中的关键环节之一是需要将接口子程序UserElem与标准程序进行成功连接，基于此用户自定义的开发单元即可完成，其具体连接过程详见参考文献[14]。考虑到用户自定义单元无法参与有限元单元的网格划分，故应先通过ANSYS单元库中的MESH200单元对所建立的二维有限元分析模型进行有限元网格剖分，最后再借助于ANSYS主程序中的内置命令将MESH200单元替换为对应的用户自定义二维黏弹性单元即可。

图1 二维黏弹性边界单元开发流程图

3 地震波的输入方法

在对结构-地基动力相互作用问题进行分析时，为了在黏弹性人工边界处实现地震波的输入，借助于刘晶波等[15]介绍的一种地震波的直接输入方法，即在边界节点处采用力学中脱离体的概念，其主要思想是将输入地震动转化为作用于人工边界节点上的等效荷载的方法来实现地震波动输入。人工边界节点处的等效荷载可表达为：

(7)

4 算例验证

4.1 均质半空间二维自由场算例

如图2所示均质二维弹性半空间，输入荷载方程见式(8)，荷载按照第3节介绍的方法输入，从均质二维弹性半空间的底部垂直输入，进而获得二维均质弹性半空间计算模型中典型位置的位移响应。

(8)

借助于大型通用软件ANSYS平台，在ANSYS中建立二维均质弹性半空间模型，该二维均质弹性有限域模型的计算尺寸：二维平面在长度方向的尺寸为40 m，宽度方向的尺寸为20 m。计算模型的网格剖分采用ANSYS单元库中的四边形单元Mesh200进行有限元网格划分，剖分后的网格尺寸为每份1.0 m，剖分后的网格单元数及节点数分别为880 864。计算模型的相关材料参数如表1所示,网格剖分后的有限元计算模型如图2所示，荷载输入总持时为2 s，每一步的荷载时间步长取为0.005 s。

表1 二维均质弹性半空间计算模型材料参数

图2 均质半空间二维计算模型网格剖分示意图

为了考察在ANSYS中用户自定义单元-二维黏弹性边界单元的准确性与可靠性，选取计算模型中典型位置作为观测点，五个观测点的坐标分别为点A(-20，0)，B(0，0)，C(20，0)，D(20，-20)和E(0，-20)，图3为各个观测点在动荷载作用下的位移时程响应。从观测点的位移时程曲线可以看出各个观测点的位移时程曲线变化趋势基本一致，最大偏差为1.5%；将A、B、C三点与D、E两点的位移时程响应曲线对比来看，可以看出A、B、C三个观测点的位移响应峰值接近D、E两观测点位移峰值的2倍，说明位移时程响应是由两部分作用共同引起，即入射波与反射波共同作用。这些结论均与理论解基本相吻合，因此可以说明在ANSYS中用户自定义单元——二维黏弹性边界单元的可靠性与准确性。

图3 A点—E点x向位移时程

4.2 层状均质半空间二维自由场算例

为了考察在均质层状土条件下的准确性与可靠性，借助于大型通用软件ANSYS平台，在ANSYS中建立了二维层状均质弹性半空间模型(见图4)，二维层状均质弹性半空间计算模型在长边方向上与宽度方向上的尺寸与与4.1节算例中模型的计算尺寸相同，并将计算模型分为两层来进行计算，两土层在ANSYS中用户自定义单元——二维黏弹性单元边界处的材料参数如表2所示,网格剖分后的有限元计算模型如图4所示。

表2 二维层状均质弹性半空间模型材料参数

(9)

式中：G(τ)=τ3H(τ)，τ=t/T，H(τ)为Heaviside函数分别按照黏弹性边界及黏性边界对模型进行计算，并与解析解[16]对比，观测点的竖向位移响应对比结果如图5所示。

从两观测点的竖向位移响应曲线(见图5)可以看出，观测点A与观测点B的竖向位移响应曲线变化趋势基本一致，仅位移峰值有一定的幅度差，A点的位移峰值较B点位移峰值所有减小，减小幅度为58%，其原因在于较波源的距离而言，B点比A点距离波源更近。此外，将三种不同解对比分析来看，黏弹性边界解的变化趋势及峰值与解析解的变化趋势及峰值基本一致，说明黏弹性边界计算结果具有高度的准确性与可靠性，而黏性边界解较解析解以及黏弹性边界解而言，在0.5 s前其变化趋势与峰值与其它两种解基本一致，而在0.5 s后出现了明显的整体漂移现象，进而导致计算结果失真。如上所述，黏弹性边界在考虑了无限半空间介质的弹性恢复性能后，很好地模拟了地震波在人工边界处能量的传播过程，同时也证明了在ANSYS中用户自定义单元——二维黏弹性边界单元具有一定的可靠性与准确性。

图4 层状半空间二维计算模型网格剖分示意图

图5 观测点处竖向位移响应时程曲线

4.3 地震动输入验证

为了提高地震波以等效荷载的输入形式计算结果可对比性，将考虑行波效应后的入射波位移时程作为该计算模型地表位移时程的理论解(见图6)。从观测点B的水平位移时程响应(见图7、图8)可以看出在两种不同波即P波、S波的入射下，观测点B的水平位移时程响应曲线的变化趋势及位移峰值基本一致。同时得到的观测点B水平位移时程响应与理论解的位移时程响应基本完全吻合，进一步可以说明，在ANSYS中用户自定义单元——二维黏弹性边界单元实现了在单元级别上的地震动等效荷载的施加，其施加方法及计算结果还具有较高的可靠性与准确性。

图6 入射波位移时程

图7 S波入射下B点水平位移响应

图8 P波入射下B点水平位移响应

5 结 论

(1) 通过4.1节、4.2节的计算结果可以看出，借助于大型通用有限元有限元分析软件-ANSYS-平台，基于Fortran语言开发的用户单元子程序VSB_UELf，在均质地基条件下及层状土地基条件下均表现出计算结果的稳定性，同时计算精度也能满足要求，说明了本文所开发的二维黏弹性边界单元具有较好的准确性与可靠性。

(2) 通过4.3节的计算结果可以看出，在P波及S波的输入下，计算结果均有着较好的稳定性与准确性，进而说明了本文所开发的二维黏弹性边界单元中波动输入方法是准确与可靠的。

(3) 另外，本文所开发的二维黏弹性边界单元很容易扩展至三维，对于实际工程中的土-结构相互作用问题的求解有一定的使用价值，能更好的方便于实际工程的应用。