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

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

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

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

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

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

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

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

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

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

1.1 無限元動力邊界原理

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

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

(1)

壓縮波沿x軸負向傳播，從有限區域進入無限元邊界，則邊界位移應為：

(2)

同時，在入射位移波經過邊界節點后反射的位移波如下：

(3)

入射和反射位移總和：

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

(4)

入射和反射速度總和：

(5)

根據彈性力學可知：

(6)

(7)

此時，邊界節點上阻尼應力為：

(8)

式中CBN為阻尼器參數。

經過人工邊界后，散射波產生的應力應與阻尼應力相等，從而消除散射波影響，即σx=σdamp，所以：

(9)

由上式計算可得：

(10)

(10)

同理可得剪切波的無限單元內嵌阻尼器系數CBT：

CBT=ρcs

(11)

式中cs為壓縮波波速。

1.2 無限元邊界上等效節點力求解

當外部地震荷載傳播進入有限元區域后，地基的運動由入射波和反射波疊加組成。散射回有限元邊界的波由無限元吸收，同時體現地基的彈性作用。假定邊界區域彈性小變形，可將荷載轉化為邊界上節點等效應力，從而解決外源波入射問題。

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

(12)

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

1.3 耦合模型的建立

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

a 二維單元

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

2 Python腳本實現無限元邊界荷載生成及施加

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

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

2.1 ABAQUS中Python腳本概述

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

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

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

在實現本文方法時，需要在Python中聲明導入一些基本模塊，以便使用ABAQUS中各對象。

from abaqus import * #導入ABAQUS模塊所有公共對象abaqus。

from abaqusConstants import * #導入所有符號常量abaqusConstants。

from caeModules import * #導入caeModules窗口，實現ABAQUS窗口中所有對象模塊的導入。

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

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

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

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

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

3) 利用文件操作函數，并使用循環將各方向速度時程、節點影響面積等文件讀入字典當中備用。

4)獲得人工邊界上各節點的坐標值，并進行分組處理，以此來判斷節點所在模型哪個側面。

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

程序實現流程如圖4所示。

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

3 算例驗證

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

(13)

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

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

a 有限元區域

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

圖6 模型監測點豎向位移示意

圖7 模型監測點水平向位移示意

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

4 結語

本文進行了基于ABAQUS-Python無限元的動力人工邊界研究。基于波在彈性均勻介質中傳播理論，推導了在無限元邊界上各節點荷載，將加速度、位移荷載轉化為等效應力，并進一步得到有限元模型中集中荷載時程，采用ABAQUS內嵌的原生腳本語言進行二次編程開發，實現了計算模型人工邊界節點上荷載的快速生成和準確施加。通過小算例比較了粘彈性人工邊界、本文無限元邊界和理論值，結果表明本文無限元-有限元模型的無限元人工邊界實現方法具有很高的精度，且二次開發的程序具有代碼量少、便于遷移應用等特點。