黨康寧，蘇晨輝，肖 瑜，張靜宜

(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)證了提出方法的精確性。

1 基于波動(dòng)理論的無(wú)限元?jiǎng)恿吔?/h2>

工程結(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)行研究。

1.1 無(wú)限元?jiǎng)恿吔缭?/h3>

無(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為壓縮波波速。

1.2 無(wú)限元邊界上等效節(jié)點(diǎn)力求解

當(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)去掉即可。

1.3 耦合模型的建立

目前，無(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ì)算模型示意

2 Python腳本實(shí)現(xiàn)無(wú)限元邊界荷載生成及施加

前節(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ǔ)，能有效減少前處理工作量。

2.1 ABAQUS中Python腳本概述

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)格操作。

2.2 無(wú)限元邊界荷載生成和施加實(shí)現(xiàn)的主要流程

模型通過(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 荷載生成和施加程序流程示意

3 算例驗(yàn)證

為驗(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)槟Ｐ拓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)題。

4 結(jié)語(yǔ)

本文進(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)。