基于ABAQUS-Python無限元的動力人工邊界研究

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

(1.陜西省引漢濟(jì)渭工程建設(shè)有限公司，陜西 西安 710010；2. 成都市龍泉驛區(qū)水務(wù)局，四川 成都 610100)

工程結(jié)構(gòu)體系動力響應(yīng)分析中，為考慮無限地基輻射阻尼效應(yīng)，各類人工邊界被廣泛研究和應(yīng)用。其中，粘彈性人工邊界得到了學(xué)者們的廣泛采用和認(rèn)可，但在有限元模型建立時，需要在粘彈性人工邊界上建立大量三向彈簧阻尼器，工作量很大。其參數(shù)的設(shè)置和集中等效荷載時程的添加難度很大。但是，近年來學(xué)者們采用無限元技術(shù)與有限元結(jié)合，開啟了相關(guān)研究工作。無限元由Ungless[1]提出，Lysmer[2]、Zienkiewicz[3]、Bettess[4]、Beer等[5]進(jìn)行了發(fā)展和改進(jìn)，國內(nèi)葛修潤[6]，張楚漢等[7]也較早進(jìn)行了這方面研究。

近年來，一些學(xué)者[8-12]利用ABAQUS有限元軟件中的無限單元，通過有限元—無限元(FE-IE)方式研究了水電站廠房、重力壩、拱橋、地下隧道等結(jié)構(gòu)反應(yīng)，取得較好結(jié)果。實現(xiàn)過程或采用編制INP輸入文件，或采用FORTRAN單獨編寫程序?qū)崿F(xiàn)地震波動荷載的生成，此過程十分繁瑣并容易出錯，增加了無限元動力人工邊界的實現(xiàn)難度。

Python作為目前最流行的腳本語言之一，具有簡潔、跨平臺等優(yōu)點，是ABAQUS程序前后處理層次的接口語言，可獲取ABAQUS模型中信息，并通過循環(huán)語句、內(nèi)置函數(shù)等能實現(xiàn)模型、荷載的高效操縱，從而便捷解決無限元動力人工邊界的繁瑣問題。

本文基于波動理論，利用ABAQUS軟件的無限單元，建立無限元—有限元模型，通過Python編程二次開發(fā)，實現(xiàn)各側(cè)人工邊界上節(jié)點的剛度參數(shù)施加、波動荷載生成以及荷載的自動施加，極大減小了前處理工作量，最后，通過算例驗證了提出方法的精確性。

1 基于波動理論的無限元動力邊界

工程結(jié)構(gòu)地震動力響應(yīng)分析模擬時，地震作用作為一種外源輸入，在截取的地基范圍內(nèi)存在地震波入射和反射，基于彈性介質(zhì)的波動理論，結(jié)合邊界節(jié)點等效荷載的可疊加原理，可解決外源波的入射問題。在有限域人工邊界通過引入一種特殊的有限元—無限元，該單元與有限元無縫銜接，并通過幾何映射，在局部坐標(biāo)中構(gòu)造插值形狀函數(shù)，實現(xiàn)計算范圍趨于無限遠(yuǎn)。無限元邊界法作為有限元方法的補(bǔ)充，具有良好的“協(xié)調(diào)性”，其對復(fù)雜散射波動的控制能力更優(yōu)[11]。

目前，無限元邊界實現(xiàn)較多是無限元—有限元相結(jié)合的方式，即在工程結(jié)構(gòu)的地基有限元區(qū)域邊界外通過一層無限元連接，并通過一定的衰減函數(shù)實現(xiàn)能量的吸收。無限元靜動力分析理論分別基于Zienkiewicz[13]和Lysmer[2]等進(jìn)行研究。

1.1 無限元動力邊界原理

無限單元可以充當(dāng)吸收邊界，單元設(shè)置了阻尼矩陣，根據(jù)荷載情況自動計算值大小。其動力分析吸收散射波原理如下。

地震壓縮波(P波)入射時，其在均質(zhì)無限彈性體中運動方程為：

(1)

壓縮波沿x軸負(fù)向傳播，從有限區(qū)域進(jìn)入無限元邊界，則邊界位移應(yīng)為：

(2)

同時，在入射位移波經(jīng)過邊界節(jié)點后反射的位移波如下：

(3)

入射和反射位移總和：

ux=f1(x-cpt)+f2(x+cpt)

(4)

入射和反射速度總和：

(5)

根據(jù)彈性力學(xué)可知：

(6)

(7)

此時，邊界節(jié)點上阻尼應(yīng)力為：

(8)

式中CBN為阻尼器參數(shù)。

經(jīng)過人工邊界后，散射波產(chǎn)生的應(yīng)力應(yīng)與阻尼應(yīng)力相等，從而消除散射波影響，即σx=σdamp，所以：

(9)

由上式計算可得：

(10)

(10)

同理可得剪切波的無限單元內(nèi)嵌阻尼器系數(shù)CBT：

CBT=ρcs

(11)

式中cs為壓縮波波速。

1.2 無限元邊界上等效節(jié)點力求解

當(dāng)外部地震荷載傳播進(jìn)入有限元區(qū)域后，地基的運動由入射波和反射波疊加組成。散射回有限元邊界的波由無限元吸收，同時體現(xiàn)地基的彈性作用。假定邊界區(qū)域彈性小變形，可將荷載轉(zhuǎn)化為邊界上節(jié)點等效應(yīng)力，從而解決外源波入射問題。

對于粘彈性人工邊界節(jié)點力公式為：

(12)

無限元邊界中已自動嵌入剛度項，因此，取公式中彈簧剛度為0，得到邊界地震動輸入的等效節(jié)點力表達(dá)式。對于從模型底邊界垂直入射的波，依據(jù)一維波動理論可分別求得邊界各節(jié)點的等效波動荷載，將粘彈性人工邊界等效節(jié)點力公式中剛度項去掉即可。

1.3 耦合模型的建立

目前，無限元人工邊界一般配合有限元模型完成，在有限元地基外部設(shè)置1層無限單元。由于阻尼項已被無限單元考慮，因此，只需添加剛度參數(shù)。對于二維和三維問題，無限單元和有限單元的節(jié)點連接形式如圖1～2所示。

a 二維單元

圖2 二維有限元-無限元半空間自由場波動計算模型示意

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

前節(jié)推導(dǎo)了等效荷載的求解過程，接下來需要計算各節(jié)點荷載時程，并在模型上施加。

作為ABAQUS的內(nèi)核語言，采用Python腳本語言進(jìn)行二次開發(fā)具有天然優(yōu)勢，通過增加特有對象模型，Python能夠直接與ABAQUS模型交互數(shù)據(jù)，讀取模型信息，并更改模型設(shè)置，為實現(xiàn)有限元、無限元交界處等效荷載的生成和施加奠定了基礎(chǔ)，能有效減少前處理工作量。

2.1 ABAQUS中Python腳本概述

Python腳本語言是面向?qū)ο笳Z言，具有對象(object)、成員(member)、方法(method)、構(gòu)造函數(shù)(constructor)、類(class)、模塊(module)和字典(dictionary)等基本特征。在ABAQUS中Python還有數(shù)據(jù)庫(database)、容器(Repository)、聲明使用(Access)及路徑(Path)等特有性質(zhì)[14]。

其中數(shù)據(jù)庫負(fù)責(zé)存儲模型的各種信息，是具有ABAQUS特征的一類特殊的對象，例如，本文主要對模型數(shù)據(jù)庫進(jìn)行操作，就是mdb。mdb對象是存放ABAQUS有限元模型的根對象(見圖3)，其中Models是倉庫類型，包含有parts、rootAssembly、loads和steps等成員對象，各成員對象由含有許多下級成員對象。

圖3 ABAQUS中的mdb對象層次示意

在實現(xiàn)本文方法時，需要在Python中聲明導(dǎo)入一些基本模塊，以便使用ABAQUS中各對象。

from abaqus import * #導(dǎo)入ABAQUS模塊所有公共對象abaqus。

from abaqusConstants import * #導(dǎo)入所有符號常量abaqusConstants。

from caeModules import * #導(dǎo)入caeModules窗口，實現(xiàn)ABAQUS窗口中所有對象模塊的導(dǎo)入。

此外，本文還需要載入io模塊，用于文件操作；載入mesh模塊，用于對模型網(wǎng)格操作。

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

模型通過GUI建立，并準(zhǔn)備好部分前處理文件，然后通過Python編程實現(xiàn)無限元邊界荷載生成和施加，主要流程為：

1)采用前述方法進(jìn)行無限元-有限元模型的建立，并劃分好網(wǎng)格。通過計算節(jié)點反力得到模型四側(cè)及底面上各節(jié)點的影響面積，并分別寫入文件。準(zhǔn)備好要輸入的三向荷載的速度時程。

2) 編寫Python腳本，在上述初始化后，通過變量定義模型名稱、裝配件名稱、分析步名稱、時間間隔、材料參數(shù)、模型長寬高等。

3) 利用文件操作函數(shù)，并使用循環(huán)將各方向速度時程、節(jié)點影響面積等文件讀入字典當(dāng)中備用。

4)獲得人工邊界上各節(jié)點的坐標(biāo)值，并進(jìn)行分組處理，以此來判斷節(jié)點所在模型哪個側(cè)面。

5)對每組中節(jié)點通過前述公式得到模型節(jié)點力及波動時程荷載。定義模型荷載函數(shù)，并按方向施加節(jié)點荷載。

程序?qū)崿F(xiàn)流程如圖4所示。

圖4 荷載生成和施加程序流程示意

3 算例驗證

為驗證本文所討論的無限元邊界法準(zhǔn)確性及所編制程序正確性，采用文獻(xiàn)[15]的算例模型進(jìn)行驗證。材料彈性模量為24 MPa，剪切模量為100 MPa，泊松比為0.2，密度為1 000 kg/m3。因此，該材料剪切波速為100 m/s。模型底端作用的荷載為速度脈沖，其表達(dá)式為：

(13)

其中f=4.0，0≤t≤0.25。

建立有限元-無限元耦合模型如圖5所示，XY平面為水平向，Z坐標(biāo)軸指向為模型豎向。有限元區(qū)域見圖5a，模型尺寸為6 m×6 m×50 m(長、寬、高)，采用8節(jié)點三維實體單元離散，網(wǎng)格尺寸為1 m。在有限元模型的底部及四周包裹1層無限單元(見圖5b)，從而完成有限元-無限元耦合模型的建立(見圖5c)。取有限元模型沿Z軸中軸線上底部、中部和頂部3個點為位移監(jiān)測點，節(jié)點號分別為25，1 250，2 475。

a 有限元區(qū)域

圖6～7分別給出了模型沿高度方向底部、中部和頂部的豎向、水平位移。

圖6 模型監(jiān)測點豎向位移示意

圖7 模型監(jiān)測點水平向位移示意

由圖6～7可知，無限元邊界結(jié)果與理論值十分接近，好于粘彈性邊界的計算值，由此驗證了本文提出的建模、荷載生成和施加方法的正確性和精確性，說明基于波動理論的無限元邊界能夠很好地解決外源輸入時的地基輻射阻尼問題。

4 結(jié)語

本文進(jìn)行了基于ABAQUS-Python無限元的動力人工邊界研究?；诓ㄔ趶椥跃鶆蚪橘|(zhì)中傳播理論，推導(dǎo)了在無限元邊界上各節(jié)點荷載，將加速度、位移荷載轉(zhuǎn)化為等效應(yīng)力，并進(jìn)一步得到有限元模型中集中荷載時程，采用ABAQUS內(nèi)嵌的原生腳本語言進(jìn)行二次編程開發(fā)，實現(xiàn)了計算模型人工邊界節(jié)點上荷載的快速生成和準(zhǔn)確施加。通過小算例比較了粘彈性人工邊界、本文無限元邊界和理論值，結(jié)果表明本文無限元-有限元模型的無限元人工邊界實現(xiàn)方法具有很高的精度，且二次開發(fā)的程序具有代碼量少、便于遷移應(yīng)用等特點。