高毅超，劉　陽(yáng)，梅　真

（1.華僑大學(xué)土木工程學(xué)院，福建廈門(mén) 361021；2.福建省結(jié)構(gòu)工程與防災(zāi)重點(diǎn)實(shí)驗(yàn)室，福建廈門(mén) 361021）

基于OpenSees的二維黏彈性邊界單元開(kāi)發(fā)及應(yīng)用

高毅超1，2，劉陽(yáng)1，2，梅真1，2

（1.華僑大學(xué)土木工程學(xué)院，福建廈門(mén) 361021；2.福建省結(jié)構(gòu)工程與防災(zāi)重點(diǎn)實(shí)驗(yàn)室，福建廈門(mén) 361021）

基于有限元開(kāi)源軟件OpenSees開(kāi)發(fā)了二維黏彈性邊界單元VS2D2Bar，采用數(shù)值算例驗(yàn)證了該二維黏彈性邊界單元在程序?qū)崿F(xiàn)上的正確性，并將其應(yīng)用于重力壩地震響應(yīng)分析。重力壩算例分析結(jié)果表明：與固定邊界無(wú)質(zhì)量地基模型相比，采用黏彈性邊界單元分析得到的壩體動(dòng)力響應(yīng)峰值大幅度減小，在重力壩的地震響應(yīng)分析中必須考慮無(wú)限地基輻射阻尼效應(yīng)的影響；OpenSees適用于重力壩地震響應(yīng)分析，編程開(kāi)發(fā)和調(diào)試工作量較小，易于在復(fù)雜水工結(jié)構(gòu)的靜動(dòng)力分析研究中推廣應(yīng)用。

黏彈性邊界；地震響應(yīng)；重力壩；OpenSees

在高烈度地震作用下，重力壩與無(wú)限地基之間將會(huì)發(fā)生顯著的動(dòng)力相互作用，由此給重力壩的地震響應(yīng)分析結(jié)果帶來(lái)巨大影響。有限元法是重力壩地震響應(yīng)分析中常用的數(shù)值方法，但是有限元法只能模擬有限域，需要在計(jì)算域的截?cái)噙吔缟鲜┘尤斯み吔?，以模擬波動(dòng)在無(wú)限地基中傳播引起的輻射阻尼效應(yīng)。在工程分析中，常用的透射邊界主要有黏性邊界［1］、黏彈性邊界［2-3］和多次透射邊界［4］等。黏彈性邊界具有足夠的計(jì)算精度，并且具有簡(jiǎn)單、易于同有限元法耦合的優(yōu)點(diǎn)，因此在大壩地震響應(yīng)分析中得到廣泛應(yīng)用［5-7］。

OpenSees（open system for earthquake engineering simulation）是由美國(guó)加州大學(xué)伯克利分校研發(fā)的開(kāi)源有限元軟件，該軟件擁有豐富的材料模型庫(kù)和單元庫(kù)，并且集成了多種高效的求解算法，逐漸得到了廣大研究人員的關(guān)注。OpenSees基于面向?qū)ο蟮木幊趟枷?，具有良好的可擴(kuò)展性，非常便于研究者持續(xù)集成最新的研究成果，目前已被應(yīng)用于土木工程中建筑結(jié)構(gòu)的模擬研究［8-9］和混合試驗(yàn)?zāi)M研究［10-11］，但在水工結(jié)構(gòu)中的模擬研究則鮮有報(bào)道。

本文借助OpenSees開(kāi)源代碼和面向?qū)ο缶幊趟枷氲膬?yōu)勢(shì)，在OpenSees中開(kāi)發(fā)了二維黏彈性邊界單元，數(shù)值驗(yàn)證該單元的正確性后將其應(yīng)用于重力壩地震響應(yīng)分析。

1　二維黏彈性邊界單元的開(kāi)發(fā)

黏彈性邊界通過(guò)在截?cái)噙吔缟显O(shè)置靜力彈簧和黏滯阻尼器吸收外行波的能量。在法向邊界，二維黏彈性邊界與物理系統(tǒng)等效的彈簧系數(shù)KBN和阻尼系數(shù)CBN分別為

在切向邊界，彈簧系數(shù)KBT和阻尼系數(shù)CBT分別為

式中：G、ρ分別為介質(zhì)剪切模量和密度；R為波源至人工邊界點(diǎn)的距離；cp、cs分別為P波和S波波速；αN、αT分別為法向和切向參數(shù)，其推薦取值范圍分別為［0.8，1.2］和［0.35，0.65］，本文取αN=1.0、αT=0.5。

采用兩節(jié)點(diǎn)線單元離散截?cái)噙吔?，?yīng)用虛功原理，參照文獻(xiàn)［12］的推導(dǎo)過(guò)程，可以得到二維黏彈性邊界單元的剛度矩陣K和阻尼矩陣C分別為

式中l(wèi)為單元長(zhǎng)度。

筆者基于OpenSees的軟件框架，采用C++語(yǔ)言編寫(xiě)了二維兩節(jié)點(diǎn)黏彈性邊界單元VS2D2Bar，實(shí)現(xiàn)了黏彈性邊界與有限元法的耦合。在OpenSees中新增單元十分容易，開(kāi)發(fā)者只需要從OpenSees提供的Element類(lèi)派生一個(gè)子類(lèi)（本文為VS2D2Bar），在子類(lèi)中需要編程實(shí)現(xiàn)的核心函數(shù)主要包括：獲取單元?jiǎng)偠染仃嘖的函數(shù)getTangentStiff（）、獲取單元質(zhì)量矩陣M的函數(shù)getMass（）、獲取單元阻尼矩陣C的函數(shù)getDamp（）以及動(dòng)力分析中獲取右端荷載向量的函數(shù)getResistingForceIncInertia（）。在黏彈性邊界單元VS2D2Bar中，單元質(zhì)量矩陣為零矩陣，單元輸入?yún)?shù)包括G、ρ、R、αN和αT。

OpenSees基于面向?qū)ο蟮木幊趟枷?，程序模塊化，編程開(kāi)發(fā)和調(diào)試工作量小。以獲取單元?jiǎng)偠染仃嚭瘮?shù)getTangentStiff（）為例，該程序的流程包含以下幾步：①將剛度矩陣K置零；②計(jì)算黏彈性邊界的彈簧系數(shù)KBN和KBT、單元長(zhǎng)度l，并根據(jù)這3個(gè)參數(shù)給剛度矩陣K賦值；③將剛度矩陣K返回主程序。

2　數(shù)值驗(yàn)證

考慮圖1所示的二維半無(wú)限彈性空間，截?cái)噙吔鐚⑵鋭澐殖山鼒?chǎng)和遠(yuǎn)場(chǎng)兩個(gè)部分。為消除量綱的影響，更加客觀地驗(yàn)證黏彈性邊界的精度，本節(jié)算例中所有的參數(shù)值均為無(wú)量綱化處理后的值。近場(chǎng)域的范圍為4×2（寬×高），采用OpenSees自帶的quad單元離散，網(wǎng)格尺寸為0.1×0.1，網(wǎng)格數(shù)為800個(gè)。

圖1　均布荷載作用下的半無(wú)限彈性空間

為了同近場(chǎng)域有限元網(wǎng)格協(xié)調(diào)，截?cái)噙吔绮捎脙晒?jié)點(diǎn)線單元離散，單元類(lèi)型為本文開(kāi)發(fā)的黏彈性邊界單元VS2D2Bar，單元數(shù)量為80個(gè)。材料參數(shù)取值如下：剪切模量為1，密度為1，泊松比為0.25。黏彈性邊界中法向和切向參數(shù)取值分別為αN=1.0和αT=0.5，邊界點(diǎn)至散射波源的距離取平均值。

采用有限元擴(kuò)展網(wǎng)格解作為參考解，擴(kuò)展網(wǎng)格的離散范圍為16×8（寬×高），P波在該范圍內(nèi)往返傳播的時(shí)間約為10，即計(jì)算時(shí)間前10范圍內(nèi)的計(jì)算結(jié)果有效。網(wǎng)格尺寸為0.1×0.1，網(wǎng)格數(shù)為12800個(gè)，截?cái)噙吔鐬樽杂蛇吔纭?/p>

考慮均布荷載p（t）為圖2所示的三角脈沖荷載，計(jì)算步長(zhǎng)取0.01。三角脈沖荷載作用下，觀測(cè)點(diǎn)A、B、C和D的豎向位移uy的時(shí)程如圖3所示。從圖3可以看出，除了位于截?cái)噙吔绲挠^測(cè)點(diǎn)C外，黏彈性邊界解與擴(kuò)展網(wǎng)格解吻合良好，這與文獻(xiàn)［12］的結(jié)論相吻合，可見(jiàn)本文開(kāi)發(fā)的二維黏彈性邊界單元在程序?qū)崿F(xiàn)上是正確的。

圖2　三角脈沖荷載

3　重力壩算例

某重力壩壩高144 m，地基離散范圍為深度方向約1.4倍壩高，向上、下游各延伸約1.4倍壩高，重力壩-地基系統(tǒng)網(wǎng)格如圖4所示，模型的節(jié)點(diǎn)總數(shù)為3369個(gè)。壩體和地基均采用OpenSees自帶的quad單元模擬，其中壩體采用四節(jié)點(diǎn)平面應(yīng)力單元模擬，單元數(shù)為565個(gè)；地基采用四節(jié)點(diǎn)平面應(yīng)變單元模擬，單元數(shù)為2 665個(gè)。截?cái)噙吔绮捎脙晒?jié)點(diǎn)線單元離散，單元類(lèi)型為本文開(kāi)發(fā)的黏彈性邊界單元VS2D2Bar，單元數(shù)為131個(gè)。

壩體混凝土彈性模量為22 GPa，泊松比為0.167，密度為2400kg/m3；基巖彈性模量為10GPa，泊松比為0.235，密度為2 600 kg/m3。進(jìn)行動(dòng)力時(shí)程分析時(shí)，動(dòng)彈性模量在靜彈性模量的基礎(chǔ)上提高30%。黏彈性邊界單元參數(shù)中，波源至邊界點(diǎn)的距離統(tǒng)一取平均距離300 m。采用El-Centro地震波作為輸入波，算例中僅考慮水平向（順河向）地震作用，地震波的加速度時(shí)程曲線如圖5所示，其速度時(shí)程和位移時(shí)程可以通過(guò)加速度時(shí)程積分得到。

圖3　豎向位移時(shí)程

圖4　重力壩-地基系統(tǒng)網(wǎng)格

圖5　El-Centro地震波加速度時(shí)程

計(jì)算中考慮了兩種地基模型：①傳統(tǒng)的無(wú)質(zhì)量地基模型，邊界為固定邊界，地震荷載以慣性力的方式施加到壩體單元上；②黏彈性邊界模型，考慮地基質(zhì)量的影響，在截?cái)噙吔缟鲜┘羽椥赃吔?，地震荷載按文獻(xiàn)［13］的方式輸入，即在底部和側(cè)邊界換算成等效節(jié)點(diǎn)荷載，實(shí)現(xiàn)波動(dòng)輸入。

在順河向El-Centro地震波作用下，壩體關(guān)鍵部位A點(diǎn)（壩頂）的動(dòng)位移（水平位移ux、豎向位移uy）時(shí)程如圖6所示。由于固定邊界的無(wú)質(zhì)量地基模型得到的是相對(duì)位移，因此將黏彈性邊界模型的計(jì)算位移扣除邊界位移后，再與無(wú)質(zhì)量地基模型的計(jì)算結(jié)果進(jìn)行對(duì)比。從圖6可以看出，與傳統(tǒng)的固定邊界無(wú)質(zhì)量地基模型相比，考慮了黏彈性邊界后，壩體的位移響應(yīng)峰值大幅降低。B點(diǎn)（壩踵）和C點(diǎn)（壩址）的動(dòng)應(yīng)力響應(yīng)（最大主應(yīng)力σ1、最小主應(yīng)力σ3）時(shí)程分別如圖7和圖8所示，同樣可以看出黏彈性邊界模型的計(jì)算結(jié)果均遠(yuǎn)小于無(wú)質(zhì)量地基模型的計(jì)算結(jié)果，可見(jiàn)傳統(tǒng)的固定邊界無(wú)質(zhì)量地基模型高估了壩體的動(dòng)力響應(yīng)。因此，在重力壩的地震響應(yīng)分析中必須考慮無(wú)限地基輻射阻尼效應(yīng)的影響，這與文獻(xiàn)［14］的結(jié)論是一致的。

圖6　A點(diǎn)位移時(shí)程

圖7　B點(diǎn)最大主應(yīng)力時(shí)程

同時(shí)算例結(jié)果表明，OpenSees適用于重力壩等水工結(jié)構(gòu)的地震響應(yīng)分析。此外，在OpenSees中集成最新研究成果的編程開(kāi)發(fā)、調(diào)試等工作量較小，非常易于將OpenSees推廣應(yīng)用到復(fù)雜水工結(jié)構(gòu)的靜動(dòng)力分析中。

圖8　C點(diǎn)最小主應(yīng)力時(shí)程

4　結(jié) 語(yǔ)

以開(kāi)源有限元軟件OpenSees為平臺(tái)，借助其開(kāi)源代碼和面向?qū)ο缶幊趟枷氲膬?yōu)勢(shì)，開(kāi)發(fā)了二維黏彈性邊界單元VS2D2Bar。數(shù)值驗(yàn)證了二維黏彈性邊界單元VS2D2Bar程序的正確性，并將該單元應(yīng)用于重力壩地震響應(yīng)分析。與固定邊界無(wú)質(zhì)量地基模型相比，黏彈性邊界模型計(jì)算得到的壩體動(dòng)力響應(yīng)峰值大幅降低，因此在重力壩的地震響應(yīng)分析中必須考慮無(wú)限地基輻射阻尼效應(yīng)的影響。OpenSees適用于重力壩地震響應(yīng)分析，并且在該軟件平臺(tái)上集成最新研究成果編程開(kāi)發(fā)和調(diào)試的工作量較小，可在復(fù)雜水工結(jié)構(gòu)的靜動(dòng)力分析研究中推廣應(yīng)用。

［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］陳燈紅，杜成斌，苑舉衛(wèi).基于ABAQUS的黏彈性邊界單元及在重力壩抗震分析中的應(yīng)用［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］邱流潮，劉樺，金峰.混凝土重力壩地震響應(yīng)數(shù)值分析方法［J］.河海大學(xué)學(xué)報(bào)（自然科學(xué)版），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］韓小雷，陳學(xué)偉，戴金華，等.基于OpenSees的剪力墻低周往復(fù)試驗(yàn)的數(shù)值分析［J］.華南理工大學(xué)學(xué)報(bào)（自然科學(xué)版），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］.工程力學(xué)，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混合模擬試驗(yàn)技術(shù)發(fā)展與應(yīng)用［J］.力學(xué)進(jìn)展，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］許國(guó)山，郝偉，陳永盛，等.基于有限元軟件OpenSees和OpenFresco的混合試驗(yàn)［J］.土木工程學(xué)報(bào)，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］.巖土工程學(xué)報(bào)，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］杜修力，趙密.基于黏彈性邊界的拱壩地震反應(yīng)分析方法［J］.水利學(xué)報(bào)，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］陳厚群.壩址地震動(dòng)輸入機(jī)制探討［J］.水利學(xué)報(bào)，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

國(guó)家自然科學(xué)基金（51409107，51378228）；中國(guó)博士后科學(xué)基金（2014M551832）

高毅超（1985—），男，講師，博士，主要從事結(jié)構(gòu)仿真研究。E-mail：gyc@hqu.edu.cn

（20150823 編輯：熊水斌）