基于OpenSees的二维黏弹性边界单元开发及应用

高毅超，刘　阳，梅　真

（1.华侨大学土木工程学院，福建厦门 361021；2.福建省结构工程与防灾重点实验室，福建厦门 361021）

基于OpenSees的二维黏弹性边界单元开发及应用

高毅超1，2，刘阳1，2，梅真1，2

（1.华侨大学土木工程学院，福建厦门 361021；2.福建省结构工程与防灾重点实验室，福建厦门 361021）

基于有限元开源软件OpenSees开发了二维黏弹性边界单元VS2D2Bar，采用数值算例验证了该二维黏弹性边界单元在程序实现上的正确性，并将其应用于重力坝地震响应分析。重力坝算例分析结果表明：与固定边界无质量地基模型相比，采用黏弹性边界单元分析得到的坝体动力响应峰值大幅度减小，在重力坝的地震响应分析中必须考虑无限地基辐射阻尼效应的影响；OpenSees适用于重力坝地震响应分析，编程开发和调试工作量较小，易于在复杂水工结构的静动力分析研究中推广应用。

黏弹性边界；地震响应；重力坝；OpenSees

在高烈度地震作用下，重力坝与无限地基之间将会发生显著的动力相互作用，由此给重力坝的地震响应分析结果带来巨大影响。有限元法是重力坝地震响应分析中常用的数值方法，但是有限元法只能模拟有限域，需要在计算域的截断边界上施加人工边界，以模拟波动在无限地基中传播引起的辐射阻尼效应。在工程分析中，常用的透射边界主要有黏性边界［1］、黏弹性边界［2-3］和多次透射边界［4］等。黏弹性边界具有足够的计算精度，并且具有简单、易于同有限元法耦合的优点，因此在大坝地震响应分析中得到广泛应用［5-7］。

OpenSees（open system for earthquake engineering simulation）是由美国加州大学伯克利分校研发的开源有限元软件，该软件拥有丰富的材料模型库和单元库，并且集成了多种高效的求解算法，逐渐得到了广大研究人员的关注。OpenSees基于面向对象的编程思想，具有良好的可扩展性，非常便于研究者持续集成最新的研究成果，目前已被应用于土木工程中建筑结构的模拟研究［8-9］和混合试验模拟研究［10-11］，但在水工结构中的模拟研究则鲜有报道。

本文借助OpenSees开源代码和面向对象编程思想的优势，在OpenSees中开发了二维黏弹性边界单元，数值验证该单元的正确性后将其应用于重力坝地震响应分析。

1　二维黏弹性边界单元的开发

黏弹性边界通过在截断边界上设置静力弹簧和黏滞阻尼器吸收外行波的能量。在法向边界，二维黏弹性边界与物理系统等效的弹簧系数KBN和阻尼系数CBN分别为

在切向边界，弹簧系数KBT和阻尼系数CBT分别为

式中：G、ρ分别为介质剪切模量和密度；R为波源至人工边界点的距离；cp、cs分别为P波和S波波速；αN、αT分别为法向和切向参数，其推荐取值范围分别为［0.8，1.2］和［0.35，0.65］，本文取αN=1.0、αT=0.5。

采用两节点线单元离散截断边界，应用虚功原理，参照文献［12］的推导过程，可以得到二维黏弹性边界单元的刚度矩阵K和阻尼矩阵C分别为

式中l为单元长度。

笔者基于OpenSees的软件框架，采用C++语言编写了二维两节点黏弹性边界单元VS2D2Bar，实现了黏弹性边界与有限元法的耦合。在OpenSees中新增单元十分容易，开发者只需要从OpenSees提供的Element类派生一个子类（本文为VS2D2Bar），在子类中需要编程实现的核心函数主要包括：获取单元刚度矩阵K的函数getTangentStiff（）、获取单元质量矩阵M的函数getMass（）、获取单元阻尼矩阵C的函数getDamp（）以及动力分析中获取右端荷载向量的函数getResistingForceIncInertia（）。在黏弹性边界单元VS2D2Bar中，单元质量矩阵为零矩阵，单元输入参数包括G、ρ、R、αN和αT。

OpenSees基于面向对象的编程思想，程序模块化，编程开发和调试工作量小。以获取单元刚度矩阵函数getTangentStiff（）为例，该程序的流程包含以下几步：①将刚度矩阵K置零；②计算黏弹性边界的弹簧系数KBN和KBT、单元长度l，并根据这3个参数给刚度矩阵K赋值；③将刚度矩阵K返回主程序。

2　数值验证

考虑图1所示的二维半无限弹性空间，截断边界将其划分成近场和远场两个部分。为消除量纲的影响，更加客观地验证黏弹性边界的精度，本节算例中所有的参数值均为无量纲化处理后的值。近场域的范围为4×2（宽×高），采用OpenSees自带的quad单元离散，网格尺寸为0.1×0.1，网格数为800个。

图1　均布荷载作用下的半无限弹性空间

为了同近场域有限元网格协调，截断边界采用两节点线单元离散，单元类型为本文开发的黏弹性边界单元VS2D2Bar，单元数量为80个。材料参数取值如下：剪切模量为1，密度为1，泊松比为0.25。黏弹性边界中法向和切向参数取值分别为αN=1.0和αT=0.5，边界点至散射波源的距离取平均值。

采用有限元扩展网格解作为参考解，扩展网格的离散范围为16×8（宽×高），P波在该范围内往返传播的时间约为10，即计算时间前10范围内的计算结果有效。网格尺寸为0.1×0.1，网格数为12800个，截断边界为自由边界。

考虑均布荷载p（t）为图2所示的三角脉冲荷载，计算步长取0.01。三角脉冲荷载作用下，观测点A、B、C和D的竖向位移uy的时程如图3所示。从图3可以看出，除了位于截断边界的观测点C外，黏弹性边界解与扩展网格解吻合良好，这与文献［12］的结论相吻合，可见本文开发的二维黏弹性边界单元在程序实现上是正确的。

图2　三角脉冲荷载

3　重力坝算例

某重力坝坝高144 m，地基离散范围为深度方向约1.4倍坝高，向上、下游各延伸约1.4倍坝高，重力坝-地基系统网格如图4所示，模型的节点总数为3369个。坝体和地基均采用OpenSees自带的quad单元模拟，其中坝体采用四节点平面应力单元模拟，单元数为565个；地基采用四节点平面应变单元模拟，单元数为2 665个。截断边界采用两节点线单元离散，单元类型为本文开发的黏弹性边界单元VS2D2Bar，单元数为131个。

坝体混凝土弹性模量为22 GPa，泊松比为0.167，密度为2400kg/m3；基岩弹性模量为10GPa，泊松比为0.235，密度为2 600 kg/m3。进行动力时程分析时，动弹性模量在静弹性模量的基础上提高30%。黏弹性边界单元参数中，波源至边界点的距离统一取平均距离300 m。采用El-Centro地震波作为输入波，算例中仅考虑水平向（顺河向）地震作用，地震波的加速度时程曲线如图5所示，其速度时程和位移时程可以通过加速度时程积分得到。

图3　竖向位移时程

图4　重力坝-地基系统网格

图5　El-Centro地震波加速度时程

计算中考虑了两种地基模型：①传统的无质量地基模型，边界为固定边界，地震荷载以惯性力的方式施加到坝体单元上；②黏弹性边界模型，考虑地基质量的影响，在截断边界上施加黏弹性边界，地震荷载按文献［13］的方式输入，即在底部和侧边界换算成等效节点荷载，实现波动输入。

在顺河向El-Centro地震波作用下，坝体关键部位A点（坝顶）的动位移（水平位移ux、竖向位移uy）时程如图6所示。由于固定边界的无质量地基模型得到的是相对位移，因此将黏弹性边界模型的计算位移扣除边界位移后，再与无质量地基模型的计算结果进行对比。从图6可以看出，与传统的固定边界无质量地基模型相比，考虑了黏弹性边界后，坝体的位移响应峰值大幅降低。B点（坝踵）和C点（坝址）的动应力响应（最大主应力σ1、最小主应力σ3）时程分别如图7和图8所示，同样可以看出黏弹性边界模型的计算结果均远小于无质量地基模型的计算结果，可见传统的固定边界无质量地基模型高估了坝体的动力响应。因此，在重力坝的地震响应分析中必须考虑无限地基辐射阻尼效应的影响，这与文献［14］的结论是一致的。

图6　A点位移时程

图7　B点最大主应力时程

同时算例结果表明，OpenSees适用于重力坝等水工结构的地震响应分析。此外，在OpenSees中集成最新研究成果的编程开发、调试等工作量较小，非常易于将OpenSees推广应用到复杂水工结构的静动力分析中。

图8　C点最小主应力时程

4　结 语

以开源有限元软件OpenSees为平台，借助其开源代码和面向对象编程思想的优势，开发了二维黏弹性边界单元VS2D2Bar。数值验证了二维黏弹性边界单元VS2D2Bar程序的正确性，并将该单元应用于重力坝地震响应分析。与固定边界无质量地基模型相比，黏弹性边界模型计算得到的坝体动力响应峰值大幅降低，因此在重力坝的地震响应分析中必须考虑无限地基辐射阻尼效应的影响。OpenSees适用于重力坝地震响应分析，并且在该软件平台上集成最新研究成果编程开发和调试的工作量较小，可在复杂水工结构的静动力分析研究中推广应用。

［1］LYSMER J，KUHLEMEYER R L.Finite dynamic model for infinitemedia［J］.JournaloftheEngineering Mechanics Division，1969，95（4）：859-877.

［2］DEEKS A J，RANDOLPH M F.Axisymmetric time-domain transmitting boundaries［J］.JournalofEngineering Mechanics，1994，120（1）：25-42.

［3］LIU J，DU Y，DU X，et al.3D viscous-spring artificial boundary in time domain［J］.Earthquake Engineering and Engineering Vibration，2006，5（1）：93-102.

［4］LIAO Z P，WONG H L.A transmitting boundary for the numerical simulation of elastic wave propagation［J］. International Journal of Soil Dynamics and Earthquake Engineering，1984，3（4）：174-183.

［5］ZHANG Chuhan，PAN Jianwen，WANG Jinting.Influence of seismic input mechanisms and radiation damping on arch dam response［J］.Soil Dynamics and Earthquake Engineering，2009，29（9）：1282-1293.

［6］陈灯红，杜成斌，苑举卫.基于ABAQUS的黏弹性边界单元及在重力坝抗震分析中的应用［J］.世界地震工程，2010，26（3）：127-132.（CHEN Denghong，DU Chengbin，YUAN Juwei.Viscous-spring boundary element based on ABAQUS and application to dynamic analysis of a gravity dam［J］.World Earthquake Engineering，2010，26（3）：127-132.（in Chinese））

［7］邱流潮，刘桦，金峰.混凝土重力坝地震响应数值分析方法［J］.河海大学学报（自然科学版），2006，34（1）：65-68.（QIU Liuchao，LIU Hua，JIN Feng.Numerical method for analysis of seismic response of concrete gravity dams［J］.JournalofHohaiUniversity（Natural Sciences），2006，34（1）：65-68.（in Chinese））

［8］韩小雷，陈学伟，戴金华，等.基于OpenSees的剪力墙低周往复试验的数值分析［J］.华南理工大学学报（自然科学版），2008，36（12）：7-12.（HAN Xiaolei，CHEN Xuewei，DAI Jinhua，et al.Numerical analysis for low cyclic loading test of shear walls based on OpenSees［J］. Journal of South China University of Technology（Natural Science Edition），2008，36（12）：7-12.（in Chinese））

［9］解琳琳，黄羽立，陆新征，等.基于OpenSees的RC框架核心筒超高层建筑抗震弹塑性分析［J］.工程力学，2014，31（1）：64-71.（XIE Linlin，HUANG Yuli，LU Xinzheng，et al.Elasto-plastic analysis for super tall RC frame-core tube structures based on OpenSees［J］. Engineering Mechanics，2014，31（1）：64-71.（in Chinese））

［10］何政，蒋碧聪.OpenSees混合模拟试验技术发展与应用［J］.力学进展，2012，42（6）：804-820.（HE Zheng，JIANG Bicong.Development and applications of hybrid simulation and experimentation using OpenSees［J］. Advances in Mechanics，2012，42（6）：804-820.（in Chinese））

［11］许国山，郝伟，陈永盛，等.基于有限元软件OpenSees和OpenFresco的混合试验［J］.土木工程学报，2012，45（增刊1）：36-41.（XU Guoshan，HAO Wei，CHEN Yongsheng，et al.Hybrid testing based on finite element software OpenSees and OpenFresco［J］.China Civil Engineering Journal，2012，45（Sup1）：36-41.（in Chinese））

［12］刘晶波，谷音，杜义欣.一致黏弹性人工边界及黏弹性边界单元［J］.岩土工程学报，2006，28（9）：1070-1075.（LIU Jingbo，GU Yin，DU Yixin.Consistent viscousspring artificial boundaries and viscous-spring boundary elements［J］.ChineseJournalofGeotechnical Engineering，2006，28（9）：1070-1075.（in Chinese））

［13］杜修力，赵密.基于黏弹性边界的拱坝地震反应分析方法［J］.水利学报，2006，37（9）：1063-1069.（DU Xiuli，ZHAO Mi.Analysis method for seismic response of arch dams in time domain based on viscous-spring artificial boundary condition［J］.Journal of Hydraulic Engineering，2006，37（9）：1063-1069.（in Chinese））

［14］陈厚群.坝址地震动输入机制探讨［J］.水利学报，2006，37（12）：1417-1423.（CHEN Houqun.Discussion on seismic input mechanism at dam site［J］.Journal of Hydraulic Engineering，2006，37（12）：1417-1423.（in Chinese））

Development of two-dimensional viscoelastic boundary element based on OpenSees and its application

GAO Yichao1，2，LIU Yang1，2，MEI Zhen1，2
（1.College of Civil Engineering，Huaqiao University，Xiamen 361021，China；2.Key Laboratory for Structural Engineering and Disaster Prevention of Fujian Province，Xiamen 361021，China）

A two-dimensional（2D）viscoelastic boundary element，VS2D2Bar，was developed based on the open source finite element package OpenSees.The applicability of the 2D viscoelastic boundary element was verified with a numerical example，and it was applied to the seismic response analysis of a gravity dam.Numerical results demonstrate that the peak dynamic response of the dam body according to the viscoelastic boundary element model greatly decreases in contrast to that obtained from the massless foundation model with fixed boundary conditions.As a result，the radiation damping effect of the infinite foundation should be considered in seismic response analysis of gravity dams.The results also show that OpenSees is applicable to the seismic response analysis of gravity dams.With a high efficiency in programming and debugging，OpenSees can be easily extended to the static and dynamic analysis of complex hydraulic structures.

viscoelastic boundary；seismic response；gravity dam；OpenSees

TV312

A

10067647（2016）05006504

10.3880/j.issn.10067647.2016.05.012

国家自然科学基金（51409107，51378228）；中国博士后科学基金（2014M551832）

高毅超（1985—），男，讲师，博士，主要从事结构仿真研究。E-mail：gyc@hqu.edu.cn

（20150823 编辑：熊水斌）