黨康寧,蘇晨輝,肖 瑜,張靜宜
(1.陜西省引漢濟(jì)渭工程建設(shè)有限公司,陜西 西安 710010;2. 成都市龍泉驛區(qū)水務(wù)局,四川 成都 610100)
工程結(jié)構(gòu)體系動(dòng)力響應(yīng)分析中,為考慮無(wú)限地基輻射阻尼效應(yīng),各類人工邊界被廣泛研究和應(yīng)用。其中,粘彈性人工邊界得到了學(xué)者們的廣泛采用和認(rèn)可,但在有限元模型建立時(shí),需要在粘彈性人工邊界上建立大量三向彈簧阻尼器,工作量很大。其參數(shù)的設(shè)置和集中等效荷載時(shí)程的添加難度很大。但是,近年來(lái)學(xué)者們采用無(wú)限元技術(shù)與有限元結(jié)合,開(kāi)啟了相關(guān)研究工作。無(wú)限元由Ungless[1]提出,Lysmer[2]、Zienkiewicz[3]、Bettess[4]、Beer等[5]進(jìn)行了發(fā)展和改進(jìn),國(guó)內(nèi)葛修潤(rùn)[6],張楚漢等[7]也較早進(jìn)行了這方面研究。
近年來(lái),一些學(xué)者[8-12]利用ABAQUS有限元軟件中的無(wú)限單元,通過(guò)有限元—無(wú)限元(FE-IE)方式研究了水電站廠房、重力壩、拱橋、地下隧道等結(jié)構(gòu)反應(yīng),取得較好結(jié)果。實(shí)現(xiàn)過(guò)程或采用編制INP輸入文件,或采用FORTRAN單獨(dú)編寫(xiě)程序?qū)崿F(xiàn)地震波動(dòng)荷載的生成,此過(guò)程十分繁瑣并容易出錯(cuò),增加了無(wú)限元?jiǎng)恿θ斯み吔绲膶?shí)現(xiàn)難度。
Python作為目前最流行的腳本語(yǔ)言之一,具有簡(jiǎn)潔、跨平臺(tái)等優(yōu)點(diǎn),是ABAQUS程序前后處理層次的接口語(yǔ)言,可獲取ABAQUS模型中信息,并通過(guò)循環(huán)語(yǔ)句、內(nèi)置函數(shù)等能實(shí)現(xiàn)模型、荷載的高效操縱,從而便捷解決無(wú)限元?jiǎng)恿θ斯み吔绲姆爆崋?wèn)題。
本文基于波動(dòng)理論,利用ABAQUS軟件的無(wú)限單元,建立無(wú)限元—有限元模型,通過(guò)Python編程二次開(kāi)發(fā),實(shí)現(xiàn)各側(cè)人工邊界上節(jié)點(diǎn)的剛度參數(shù)施加、波動(dòng)荷載生成以及荷載的自動(dòng)施加,極大減小了前處理工作量,最后,通過(guò)算例驗(yàn)證了提出方法的精確性。
工程結(jié)構(gòu)地震動(dòng)力響應(yīng)分析模擬時(shí),地震作用作為一種外源輸入,在截取的地基范圍內(nèi)存在地震波入射和反射,基于彈性介質(zhì)的波動(dòng)理論,結(jié)合邊界節(jié)點(diǎn)等效荷載的可疊加原理,可解決外源波的入射問(wèn)題。在有限域人工邊界通過(guò)引入一種特殊的有限元—無(wú)限元,該單元與有限元無(wú)縫銜接,并通過(guò)幾何映射,在局部坐標(biāo)中構(gòu)造插值形狀函數(shù),實(shí)現(xiàn)計(jì)算范圍趨于無(wú)限遠(yuǎn)。無(wú)限元邊界法作為有限元方法的補(bǔ)充,具有良好的“協(xié)調(diào)性”,其對(duì)復(fù)雜散射波動(dòng)的控制能力更優(yōu)[11]。
目前,無(wú)限元邊界實(shí)現(xiàn)較多是無(wú)限元—有限元相結(jié)合的方式,即在工程結(jié)構(gòu)的地基有限元區(qū)域邊界外通過(guò)一層無(wú)限元連接,并通過(guò)一定的衰減函數(shù)實(shí)現(xiàn)能量的吸收。無(wú)限元靜動(dòng)力分析理論分別基于Zienkiewicz[13]和Lysmer[2]等進(jìn)行研究。
無(wú)限單元可以充當(dāng)吸收邊界,單元設(shè)置了阻尼矩陣,根據(jù)荷載情況自動(dòng)計(jì)算值大小。其動(dòng)力分析吸收散射波原理如下。
地震壓縮波(P波)入射時(shí),其在均質(zhì)無(wú)限彈性體中運(yùn)動(dòng)方程為:
(1)
壓縮波沿x軸負(fù)向傳播,從有限區(qū)域進(jìn)入無(wú)限元邊界,則邊界位移應(yīng)為:
(2)
同時(shí),在入射位移波經(jīng)過(guò)邊界節(jié)點(diǎn)后反射的位移波如下:
(3)
入射和反射位移總和:
ux=f1(x-cpt)+f2(x+cpt)
(4)
入射和反射速度總和:
(5)
根據(jù)彈性力學(xué)可知:
(6)
(7)
此時(shí),邊界節(jié)點(diǎn)上阻尼應(yīng)力為:
(8)
式中CBN為阻尼器參數(shù)。
經(jīng)過(guò)人工邊界后,散射波產(chǎn)生的應(yīng)力應(yīng)與阻尼應(yīng)力相等,從而消除散射波影響,即σx=σdamp,所以:
(9)
由上式計(jì)算可得:
(10)
(10)
同理可得剪切波的無(wú)限單元內(nèi)嵌阻尼器系數(shù)CBT:
CBT=ρcs
(11)
式中cs為壓縮波波速。
當(dāng)外部地震荷載傳播進(jìn)入有限元區(qū)域后,地基的運(yùn)動(dòng)由入射波和反射波疊加組成。散射回有限元邊界的波由無(wú)限元吸收,同時(shí)體現(xiàn)地基的彈性作用。假定邊界區(qū)域彈性小變形,可將荷載轉(zhuǎn)化為邊界上節(jié)點(diǎn)等效應(yīng)力,從而解決外源波入射問(wèn)題。
對(duì)于粘彈性人工邊界節(jié)點(diǎn)力公式為:
(12)
無(wú)限元邊界中已自動(dòng)嵌入剛度項(xiàng),因此,取公式中彈簧剛度為0,得到邊界地震動(dòng)輸入的等效節(jié)點(diǎn)力表達(dá)式。對(duì)于從模型底邊界垂直入射的波,依據(jù)一維波動(dòng)理論可分別求得邊界各節(jié)點(diǎn)的等效波動(dòng)荷載,將粘彈性人工邊界等效節(jié)點(diǎn)力公式中剛度項(xiàng)去掉即可。
目前,無(wú)限元人工邊界一般配合有限元模型完成,在有限元地基外部設(shè)置1層無(wú)限單元。由于阻尼項(xiàng)已被無(wú)限單元考慮,因此,只需添加剛度參數(shù)。對(duì)于二維和三維問(wèn)題,無(wú)限單元和有限單元的節(jié)點(diǎn)連接形式如圖1~2所示。
a 二維單元
圖2 二維有限元-無(wú)限元半空間自由場(chǎng)波動(dòng)計(jì)算模型示意
前節(jié)推導(dǎo)了等效荷載的求解過(guò)程,接下來(lái)需要計(jì)算各節(jié)點(diǎn)荷載時(shí)程,并在模型上施加。
作為ABAQUS的內(nèi)核語(yǔ)言,采用Python腳本語(yǔ)言進(jìn)行二次開(kāi)發(fā)具有天然優(yōu)勢(shì),通過(guò)增加特有對(duì)象模型,Python能夠直接與ABAQUS模型交互數(shù)據(jù),讀取模型信息,并更改模型設(shè)置,為實(shí)現(xiàn)有限元、無(wú)限元交界處等效荷載的生成和施加奠定了基礎(chǔ),能有效減少前處理工作量。
Python腳本語(yǔ)言是面向?qū)ο笳Z(yǔ)言,具有對(duì)象(object)、成員(member)、方法(method)、構(gòu)造函數(shù)(constructor)、類(class)、模塊(module)和字典(dictionary)等基本特征。在ABAQUS中Python還有數(shù)據(jù)庫(kù)(database)、容器(Repository)、聲明使用(Access)及路徑(Path)等特有性質(zhì)[14]。
其中數(shù)據(jù)庫(kù)負(fù)責(zé)存儲(chǔ)模型的各種信息,是具有ABAQUS特征的一類特殊的對(duì)象,例如,本文主要對(duì)模型數(shù)據(jù)庫(kù)進(jìn)行操作,就是mdb。mdb對(duì)象是存放ABAQUS有限元模型的根對(duì)象(見(jiàn)圖3),其中Models是倉(cāng)庫(kù)類型,包含有parts、rootAssembly、loads和steps等成員對(duì)象,各成員對(duì)象由含有許多下級(jí)成員對(duì)象。
圖3 ABAQUS中的mdb對(duì)象層次示意
在實(shí)現(xiàn)本文方法時(shí),需要在Python中聲明導(dǎo)入一些基本模塊,以便使用ABAQUS中各對(duì)象。
from abaqus import * #導(dǎo)入ABAQUS模塊所有公共對(duì)象abaqus。
from abaqusConstants import * #導(dǎo)入所有符號(hào)常量abaqusConstants。
from caeModules import * #導(dǎo)入caeModules窗口,實(shí)現(xiàn)ABAQUS窗口中所有對(duì)象模塊的導(dǎo)入。
此外,本文還需要載入io模塊,用于文件操作;載入mesh模塊,用于對(duì)模型網(wǎng)格操作。
模型通過(guò)GUI建立,并準(zhǔn)備好部分前處理文件,然后通過(guò)Python編程實(shí)現(xiàn)無(wú)限元邊界荷載生成和施加,主要流程為:
1)采用前述方法進(jìn)行無(wú)限元-有限元模型的建立,并劃分好網(wǎng)格。通過(guò)計(jì)算節(jié)點(diǎn)反力得到模型四側(cè)及底面上各節(jié)點(diǎn)的影響面積,并分別寫(xiě)入文件。準(zhǔn)備好要輸入的三向荷載的速度時(shí)程。
2) 編寫(xiě)Python腳本,在上述初始化后,通過(guò)變量定義模型名稱、裝配件名稱、分析步名稱、時(shí)間間隔、材料參數(shù)、模型長(zhǎng)寬高等。
3) 利用文件操作函數(shù),并使用循環(huán)將各方向速度時(shí)程、節(jié)點(diǎn)影響面積等文件讀入字典當(dāng)中備用。
4)獲得人工邊界上各節(jié)點(diǎn)的坐標(biāo)值,并進(jìn)行分組處理,以此來(lái)判斷節(jié)點(diǎn)所在模型哪個(gè)側(cè)面。
5)對(duì)每組中節(jié)點(diǎn)通過(guò)前述公式得到模型節(jié)點(diǎn)力及波動(dòng)時(shí)程荷載。定義模型荷載函數(shù),并按方向施加節(jié)點(diǎn)荷載。
程序?qū)崿F(xiàn)流程如圖4所示。
圖4 荷載生成和施加程序流程示意
為驗(yàn)證本文所討論的無(wú)限元邊界法準(zhǔn)確性及所編制程序正確性,采用文獻(xiàn)[15]的算例模型進(jìn)行驗(yàn)證。材料彈性模量為24 MPa,剪切模量為100 MPa,泊松比為0.2,密度為1 000 kg/m3。因此,該材料剪切波速為100 m/s。模型底端作用的荷載為速度脈沖,其表達(dá)式為:
(13)
其中f=4.0,0≤t≤0.25。
建立有限元-無(wú)限元耦合模型如圖5所示,XY平面為水平向,Z坐標(biāo)軸指向?yàn)槟P拓Q向。有限元區(qū)域見(jiàn)圖5a,模型尺寸為6 m×6 m×50 m(長(zhǎng)、寬、高),采用8節(jié)點(diǎn)三維實(shí)體單元離散,網(wǎng)格尺寸為1 m。在有限元模型的底部及四周包裹1層無(wú)限單元(見(jiàn)圖5b),從而完成有限元-無(wú)限元耦合模型的建立(見(jiàn)圖5c)。取有限元模型沿Z軸中軸線上底部、中部和頂部3個(gè)點(diǎn)為位移監(jiān)測(cè)點(diǎn),節(jié)點(diǎn)號(hào)分別為25,1 250,2 475。
a 有限元區(qū)域
圖6~7分別給出了模型沿高度方向底部、中部和頂部的豎向、水平位移。
圖6 模型監(jiān)測(cè)點(diǎn)豎向位移示意
圖7 模型監(jiān)測(cè)點(diǎn)水平向位移示意
由圖6~7可知,無(wú)限元邊界結(jié)果與理論值十分接近,好于粘彈性邊界的計(jì)算值,由此驗(yàn)證了本文提出的建模、荷載生成和施加方法的正確性和精確性,說(shuō)明基于波動(dòng)理論的無(wú)限元邊界能夠很好地解決外源輸入時(shí)的地基輻射阻尼問(wèn)題。
本文進(jìn)行了基于ABAQUS-Python無(wú)限元的動(dòng)力人工邊界研究?;诓ㄔ趶椥跃鶆蚪橘|(zhì)中傳播理論,推導(dǎo)了在無(wú)限元邊界上各節(jié)點(diǎn)荷載,將加速度、位移荷載轉(zhuǎn)化為等效應(yīng)力,并進(jìn)一步得到有限元模型中集中荷載時(shí)程,采用ABAQUS內(nèi)嵌的原生腳本語(yǔ)言進(jìn)行二次編程開(kāi)發(fā),實(shí)現(xiàn)了計(jì)算模型人工邊界節(jié)點(diǎn)上荷載的快速生成和準(zhǔn)確施加。通過(guò)小算例比較了粘彈性人工邊界、本文無(wú)限元邊界和理論值,結(jié)果表明本文無(wú)限元-有限元模型的無(wú)限元人工邊界實(shí)現(xiàn)方法具有很高的精度,且二次開(kāi)發(fā)的程序具有代碼量少、便于遷移應(yīng)用等特點(diǎn)。