基于无厚度Goodman单元的粘弹性人工边界模拟

罗 佩， 刘国明， 王文君, 2, 徐辰奎

(1. 福州大学土木工程学院, 福建 福州 350116； 2. 漳州职业技术学院建筑工程系, 福建 漳州 363000)



基于无厚度Goodman单元的粘弹性人工边界模拟

罗 佩1， 刘国明1， 王文君1, 2, 徐辰奎1

(1. 福州大学土木工程学院, 福建 福州 350116； 2. 漳州职业技术学院建筑工程系, 福建 漳州 363000)

基于粘弹性人工边界, 推导了三维一致粘弹性人工边界单元的刚度矩阵和阻尼矩阵. 利用边界单元厚度对计算结果影响不大， 引入了一种新的一致粘弹性单元的模拟方法， 即无厚度Goodman单元. 基于无厚度Goodman单元对一致粘弹性边界单元进行模拟， 并通过半无限空间和均匀半空间的算例进行验证. 结果表明， 使用无厚度Goodman单元实现的一致粘弹性边界单元简单、 实用性强， 并且可以得到足够的精度. 相比于有厚度的边界单元， 无厚度Goodman单元符合边界单元无厚度的实际情况， 理论完善， 且对于可能存在的边界不规则情况适用性强.

粘弹性人工边界； 无厚度Goodman单元； 人工边界单元； 动力计算

0　引言

在坝体-地基系统动力响应分析中， 地基的辐射阻尼的影响是一个备受关注的问题. 用有限化方法对无限地基进行处理的时候， 通过在划定的有限范围地基介质中设置虚拟的边界， 并确定这些虚拟边界所需满足的相应条件， 该边界称之为人工边界. 人工边界条件(artificial boundary condition)的概念是由Alterman[1]最早在1968年提出的. 目前对人工边界的研究主要分为全局人工边界条件和局部人工边界条件两大方面. 局部人工边界条件适应性强， 是时域离散形式的边界条件， 具备时空解耦的特性， 形式简单、 计算量小、 计算耗时少， 便于通过编程实现， 因此， 被广泛应用于有限元方法中.

粘弹性人工边界是一种连续的局部人工边界. 1994年， Deeks等[2]基于粘性边界的基础上提出了粘弹性人工边界. 粘弹性人工边界等效于将连续分布的弹簧元件和阻尼器元件设置于人工截断边界上， 与粘性边界相比， 模拟出了散射波辐射和地基弹性恢复性能. 其概念清晰， 便于有限元实现， 克服了粘性边界的低频失稳问题， 能够模拟地基的弹性恢复能力， 具有较高的计算效率和准确度， 以及具备良好的鲁棒性等优点， 已广泛应用于土-结构动力相互作用问题的分析中. 引入无厚度的Goodman单元模拟二维和三维一致粘弹性人工边界单元， 取得了较好的效果. 无厚度Goodman单元与有厚度的边界单元相比满足了边界单元无厚度的实际情况， 理论完善， 且对于可能存在的边界不规则情况适用性强， 施加简单.

1　粘弹性人工边界

根据球坐标系内的球面波波动方程的推导， 人工边界可以等效为连续分布的并联弹簧-阻尼器系统. 因此， 粘弹性边界的首要问题是选取适当的弹簧刚度和阻尼系数. 切向和法向的弹簧刚度和阻尼系数分别按照式(1)、 (2)来取值.

(1)

(2)

其中： KBT、 KBN分别是切向和法向弹簧刚度； CBT、 CBN分别是切向和法向的阻尼系数； R是人工边界点至波源的距离； cs、 cP分别是S波和P波的波速； G是介质的剪切模量； ρ是介质的质量密度； αT、 αN分别是粘弹性人工边界切向和法向的修正系数， 推荐取0.67、 1.33[3-4].

2　一致粘弹性边界单元

在进行有限元离散时， 可用有限元形函数将连续分布的物理元件转化为耦联的边界单元， 该方法称为一致粘弹性人工边界； 也可以通过对边界单元进行集中处理， 形成解耦的阻尼器和弹簧元件， 该方法称为集中粘弹性人工边界. 刘晶波等[5]将一致粘弹性人工边界与集中粘弹性人工边界进行对比， 得出了采用一致粘弹性人工边界时的计算精度较高的结论.

相对于采用弹簧-阻尼器元件， 等效实体边界单元对粘弹性人工边界的模拟更容易实现， 而且对人为划分地基范围时可能存在的边界不规则情况， 这种边界单元的适应性更强. 通过比较不同等效人工边界单元厚度h对计算结果的影响， 发现边界单元厚度对计算结果的影响不大[3, 5]， 厚度h可以在一个较大的范围内灵活取值. 当然， 由于人工边界是没有厚度的， 采用有厚度的实体边界单元模拟毕竟不够完善， 因此在进行有限元动力计算程序的编制时， 引入无厚度的Goodman单元模拟一致粘弹性人工边界单元， 取得了较好的效果. 下面对无厚度Goodman单元对一致粘弹性边界单元的模拟进行介绍.

2.1三维一致粘弹性边界模型单元劲度与阻尼矩阵

采用无厚度Goodman单元进行模拟二维和三维一致粘弹性边界单元. 三维模型单元采用两种无厚度Goodman单元， 即无厚度六面体8结点单元和无厚度五面体6结点单元. 采用等参单元形式表示这两种单元， 如图1所示：

无厚度六面体8结点单元形函数为：

(3)

无厚度五面体6结点单元形函数为：

(4)

位移模式：

(5)

坐标变换式：

(6)

(7)

其中：

上式， 8结点单元时，d=4； 6结点单元时，d=3.

(8)

其中：

三维的劲度矩阵和阻尼矩阵为二维积分， 采用高斯积分确定.

整体坐标系下的单元劲度矩阵k和阻尼矩阵c为:

(9)

其中， 转换矩阵：

结点转换矩阵：

高斯点雅克比行列式为：

局部坐标与整体坐标的关系为：

(10)

(11)

对于二维的模型情况易于采用退化的三维模型得到， 不再进行公式推导.

3　数值算例

采用自主开发的适用于二、 三维混凝土坝的静动力分析的动力线性有限元程序DYCDAM.f. 三维实体单元采用六面体8结点、 五面体6结点和四面体4结点等参单元， 采用无厚度的六面体8结点单元和五面体6结点单元模拟人工边界. 静荷载包括自重、 水荷载、 法向面荷载、 集中力的作用， 静力计算可以考虑分期施工与分期蓄水过程. 采用Wilson-θ法进行动力分析， 动荷载包括地震波自由场输入、 给定结点的动荷载时程以及人工边界的等效地震动输入. 与动水附加质量矩阵结合， 可考虑坝-水动力相互作用. 通过输入2的维数， 自动退化为二维计算程序. 以下通过两个算例， 分别建立二维和三维模型进行动力计算分析， 并与解析解进行对比， 以说明粘弹性人工边界对波动问题处理是有效的. 同时, 通过算例验证本程序的可行性.

3.1二维半无限空间自由边界受冲击荷载作用问题

建立二维半无限空间模型， 模型尺寸L=400 m，H=200 m. 根据对称性， 取一半即L=200 m，H=200 m进行计算， 计算区域如图2(a)所示. 材料的弹性模量为E=0.1 GPa， 泊松比ν=0.22， 密度ρ=2 000 kg·m-3， 介质的水平向波速取为143 m·s-1. 区域顶部为自由表面， 在区域的两侧边和底边处分别施加固定边界条件和等效粘弹性边界条件. 网格尺寸为10 m×10 m， 固定边界条件模型剖分为441个结点， 400个平面四结点四边形等参单元； 等效粘弹性边界模型剖分为482个结点， 440个平面四结点四边形等参单元. 为获得拟精确解， 将截断边界设置在离计算区域足够远的地方， 即将模型扩大为L=1 600 m，H=800 m(计算时取一半大小即L=800 m，H=800 m)， 以保证在计算时间内计算关心的区域不会有冲击荷载传递到截断边界后反射回来的波动影响. 网格尺寸为10 m×10 m， 共剖分6 561个结点， 6 400个平面四结点四边形等参单元.

在A点处施加如图2(b)所示的冲击荷载， 计算分析时间步长为0.01 s， 共计算600步， 历时6 s.A点(0， 0)和B点(190， 0)处的位移时程曲线图如图3所示.

从图3中可以看出， 将截断边界设置在离计算区域足够远时， 在计算时间内计算关心的区域不会有冲击荷载传递到截断边界后反射回来的波动影响， 因此观测点不会受到边界反射波的影响. 采用固定边界时， 由于计算区域的范围较小， 在计算时间内冲击荷载传递到固定边界时会产生反射， 进而对各观测点的位移时程曲线产生较大的影响， 其位移产生了震荡， 与拟精确解相比误差较大. 而在截断边界上采用等效粘弹性人工边界， 在计算区域与固定边界模型相同的情况下， 各个观测点的位移不会受到边界反射波的影响. 通过对比计算得出的位移时程曲线与拟精确解， 可以看出二者非常接近.

算例证明了本程序的合理性， 同时也证明了采用粘弹性边界在处理二维半无限空间的内源波动问题时是有效的.

3.2经典Lamb问题

考虑经典Lamb问题， 建立三维半无限空间的计算模型， 计算时各单位为无量纲. 模型尺寸为1×1×0.5， 根据对称性， 在水平方向上取其四分之一大小即0.5×0.5×0.5， 如图4(a)所示， 共划分为1 662个结点和1 300个单元. 在x、y、z方向上划分为10×10×10的网格， 网格尺寸为0.05. 模型的剪切模量G=16， 泊松比ν=0.25， 介质密度ρ=1， 时间步长为0.01. 等效一致粘弹性边界单元通过在模型的四周和底面设置厚度为0的三维Goodman单元， 并将人工边界单元的最外层结点固定. 输入的荷载为作用于半空间表面的集中荷载， 如图4(b)所示.

计算时， 人工边界参数αT=0.67，αN=1.33. S波和P波波速均取4 m·s-1. 计算结果如图5所示， 取模型上的观测点分别距加载中心O点距离为r=0.20和r=0.40， 分别计算出各个观测点在解析解、 粘弹性人工边界和固定边界下的z向位移反应.

从图5所示的位移时程曲线可以看出， 对应不同观测点， 当采用固定边界时， 位移产生了较大幅度的波动， 相对于解析解存在着很大的误差； 而使用一致粘弹性人工边界时， 在接近加载点附近的结点和靠近边界的结点上， 其计算结果与解析解非常接近， 计算的精度较高. 同时也可以看出， 在对粘弹性人工边界的模拟上， 采用等效实体粘弹性边界单元代替弹簧-阻尼器， 将使施加更为简单， 而且可以在边界不规则的情况下使用.

4　结论

基于粘弹性人工边界推导了三维一致粘弹性人工边界单元的刚度矩阵和阻尼矩阵. 通过理论推导引入了一种新的一致粘弹性单元的模拟方法， 即采用无厚度Goodman单元实现粘弹性边界单元的模拟. 数值算例表明， 使用无厚度Goodman单元实现的一致粘弹性边界单元施加简单、 实用性强， 并且具有足够的精度. 无厚度Goodman单元与有厚度的边界单元相比满足了边界单元无厚度的实际情况， 理论完善， 且对于可能存在的边界不规则情况适用性强， 施加简便， 体现出一致粘弹性人工边界在进行动力计算时的优越性.

[1] ALTERMAN Z S, KARAL F C. Propagation of elastic waves in layered media by finite-difference methods[J]. Bull Seism Soc Am, 1968, 58(1): 367-398.

[2] DEEKS A J, RANDOLPH M F. Axisymmetric time-domain transmitting boundaries[J]. Journal of Engineering Mechanics, 1994, 120(1): 25-42.

[3] 谷音, 刘晶波, 杜义欣. 三维一致粘弹性人工边界及等效粘弹性边界单元[J]. 工程力学, 2007, 24(12): 31-37.

[4] 赵建峰, 杜修力, 韩强, 等. 外源波动问题数值模拟的一种实现方式[J]. 工程力学, 2007, 24(4): 52-58.

[5] 刘晶波, 谷音, 杜义欣. 一致粘弹性人工边界及粘弹性边界单元[J]. 岩土工程学报,2006, 28(9): 1 070-1 075.

(责任编辑： 洪江星)

Simulation of viscous-spring artificial boundaries on the basis of zero-thickness Goodman element

LUO Pei1, LIU Guoming1, WANG Wenjun1, 2, XU Chenkui1

(1. College of Civil Engineering， Fuzhou University， Fuzhou， Fujian 350116， China；2. Department of Architectural Engineering, Zhangzhou Vocational Technical College, Zhangzhou, Fujian 363000， China)

The stiffness and damping matrixes of 3D consistent artificial boundary are deduced through viscous-spring artificial boundary theory. Using thickness of the boundary element has little effect on the results, the zero-thickness Goodman element is developed as a new simulation method of consistent viscous-spring artificial boundary. The viscous-spring boundary element is simulated by the zero-thickness Goodman element, which is verified through the example of the half infinite space and uniform half space problems. The results show that this method is simple, practical, and can get enough accuracy. Compared with thickness of the boundary element, the zero-thickness Goodman element conforms to the actual situation, which the boundary element is no thick, have perfect theory and well adaptability to the possible irregular boundary conditions.

viscous-spring artificial boundary; zero-thickness Goodman element; artificial boundary element; dynamic calculation

10.7631/issn.1000-2243.2016.01.0104

1000-2243(2016)01-0104-06

2014-07-17

刘国明(1963-)， 博士， 教授， 主要从事水工结构工程专业研究， lgm6379@163.com

福建省自然科学基金资助项目(2011J01309)

TV312

A