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

王 飛， 宋志強， 劉昱杰， 王 建

(西安理工大學 水利水電學院， 陜西 西安 710048)

1 研究背景

國內外關于各種動力人工邊界做了大量的研究工作，主要集中在局部人工邊界和全局人工邊界兩大類上。常用的局部人工邊界有透射邊界[1]、黏性邊界[2]、黏彈性邊界[3-4]。全局人工邊界的典型代表之一是無限元。1973年Ungless[5]首先提出了無限元思想，用以解決無限域模擬問題。Bettess等[6]根據整體坐標和局部坐標間的映射，提出了映射無限單元的概念，被稱之為Bettess單元。1992年，Bettess[7]將已有的無限元研究成果進行歸納，出版了世界上首部無限元專著《Infinite Elements》。近年一些學者開始應用無限元研究結構-地基動力相互作用問題[8-10]。Kim[11]和Yun[12]創新了動力無限元公式，在頻域和時域內研究了二維、三維層狀土-結構的相互作用。許多學者研究經驗表明，無限元與有限元協調性好，比邊界元等其他數值方法求解無限域問題更具優勢和實用性。戚玉亮等[13]對ABAQUS動力無限元人工邊界進行了改進，但未討論無限元對于靜力人工邊界以及靜-動力統一人工邊界的適用性。

對于一般符合小變形、線彈性的結構，可分別獨立進行靜力荷載和地震動力荷載作用效應分析，將二者疊加即為結構總作用效應。對于具有材料或接觸非線性、強震作用下會發生大變形的結構疊加原理原則上是不適用的，必須進行靜、動力統一分析，由此引出靜、動力邊界的轉換或者統一問題。

文獻[14]正確地給出了一種靜-動力人工邊界轉換的方法，但此方法較為繁瑣，且應用過程中容易由于未施加靜力分析產生的節點約束反力，使得結構在動力零時刻不能保持平衡，導致計算結果失真。為了避免邊界轉換的繁瑣，有些研究者開始進行靜-動力邊界的統一問題研究。

劉晶波等[15]提出了黏彈性靜-動力統一人工邊界，其物理概念清晰、計算公式明確，易于有限元實現。對于半無限空間自由表面在法向集中力和表面點波源振動作用下的靜-動力問題分析結果表明，該靜-動力統一邊界在保證動力問題分析精度和穩定性的同時，兼顧了靜力問題分析的精度，但該方法對于半無限自由表面切向集中力的靜力數值解誤差較大。

高峰等[14]用黏彈性靜-動力統一人工邊界分析半無限空間體應力，同樣認為土體內部各點的橫向應力與理論解之間存在較大偏差。在土木工程結構中，靜力效應主要是結構自重引起的，因此這一不足不會影響黏彈性靜-動力統一人工邊界在該領域中的應用。而在水工結構中，除了結構自重，靜水壓力、淤沙壓力等主要水平向荷載對結構的靜力效應是不容忽視的。

為了既方便又相對準確地計算結構的靜-動力綜合效應，本文建立了基于ABAQUS無限元的靜-動力統一人工邊界?；跓o限元在理論上滿足無窮遠處位移為零、波傳播至無限遠域衰減為零的條件[9]。推導了外源波動輸入下無限元靜-動力統一人工邊界等效節點力計算公式，并基于Python語言編寫了等效節點力計算與施加程序。該方法無需進行靜、動力邊界的轉換和應力場、邊界約束反力的導入，相比黏彈性靜-動力統一邊界條件，靜力計算的精度更高，尤其適用于靜力荷載為切向的情況，此外，基于無限元靜-動力統一人工邊界的有限元模型，地基有限域的模擬范圍可以大幅縮小，而對靜、動力計算的精度及綜合效應并無較大影響，顯著提高了計算效率。

2 基于ABAQUS無限元的外源波動輸入研究

ABAQUS分析軟件提供的無限元動力人工邊界應用了Lysmer和Kuhlemeyer提出的黏性邊界理論，區別在于阻尼器是內嵌均勻分布在ABAQUS提供的無限單元中?；诖?，無限元人工邊界下的自由場波動輸入方式可參考黏性人工邊界理論或黏彈性人工邊界理論，即將地震波在邊界處的自由場運動轉換為作用于人工邊界結點上的等效結點力進行施加如圖1所示。

圖1 黏彈性邊界示意圖

近些年來，黏彈性邊界法的應用與改進較多，在模擬無限遠域地基的輻射阻尼效應時，其精度和穩定性較黏性人工邊界高[16-18]，故本文在參考部分黏彈性邊界理論[3]的基礎上推導出無限元動力邊界法的地震動輸入方式。

黏彈性邊界法人工邊界結點上的等效結點力公式為[16]：

(1)

當公式(1)中的Kb=0時，便得到可應用于ABAQUS軟件提供的動力無限元邊界上地震動輸入的等效結點力。Cb在不同人工邊界面上有其特定的表達形式，當人工邊界面外法線方向與x軸(如圖1)平行時：

(2)

與y軸平行時：

(3)

與z軸平行時：

(4)

式中：CBN=cpρ、CBT=csρ，即內嵌均勻分布于動力無限單元內部的阻尼器的阻尼系數，與黏彈性人工邊界的阻尼系數相同，因其計算與賦值由ABAQUS軟件自動完成，與黏彈性人工邊界法比較省去了大量的前處理工作；cp、cs分別為P波和S波波速，ρ為介質質量密度。

對于三維模型具體表達式如下：

(1)對于底面

(5)

式中：等效節點力上標代表節點所在人工邊界面的外法線方向，規定與坐標軸正方向一致為正，反之為負；下標表示等效節點力分量的方向；H為邊界底面到地表的距離，h為人工邊界節點至底邊界面的距離。

(2)對于x負向邊界面

(6)

(3)對于x正向邊界面

(7)

(4)對于y負向邊界面

(8)

(5)對于y正向邊界面

(9)

3 算例分析

為了比較本文提出的無限元靜-動力統一邊界在處理靜-動力綜合問題的優越性，本文對4種方法進行對比分析。為避免混淆，規定各種方法中人工邊界的名稱，并對各種人工邊界做出簡要說明，如表1。

表1 人工邊界名稱及施加方法

3.1 半無限空間體受表面法向均布靜荷載和單位脈沖荷載聯合作用

在半無限域中截取范圍為-30.0 m≤x,y≤30.0 m、-50.0 m≤z≤0 m的模型模擬半無限空間體，坐標原點位于模型頂部中心點，如圖2。模型介質質量密度ρ=1 000 kg/m3，彈性模量E=2.4×107Pa，泊松比ν=0.2。首先在模型表面中心區域作用矩形均布靜荷載q=1×106Pa，矩形的邊長為6.0 m，在模型的底面和側面分別施加黏彈性動力人工邊界(方法1)、黏彈性靜-動力統一人工邊界(方法2)和無限元靜-動力統一人工邊界(方法3)，如圖3；接著在模型的底面垂直向上輸入沿z向的單位脈沖P波以及沿x、y向的單位脈沖S波，其位移表達式見公式(10)，位移波形圖見圖4。由介質的彈性性質可得波速cp=163.30 m/s、cs=100.00 m/s。

另建立有限元-無限元耦合模型，其中內部離散有限域范圍為-6.0 m≤x,y≤6.0 m，-50.0 m≤z≤0，在有限單元的外部包裹著無限單元，加載方式等同方法3，即方法4。采用方法1、方法2、方法3和方法4對該靜-動綜合力問題進行分析。計算時長為2.0 s，時間步長為0.01 s。

(10)

由彈性力學理論公式[19]和疊加原理，可得到模型頂部中心點O(0，0，0)、中部中心點B(0，0，-25)和底部中心點C(0，0，-50)的z向位移的理論解，圖5～7給出了4種方法下點O、B和C的z向位移數值解。

由圖5～7可知，模型在表面法向均布靜荷載和從底面垂直向上輸入的z向脈沖荷載的聯合作用下，4種方法的z向位移與理論解的偏差都較小。圖8給出了圖5中點O在0.55s之后的位移時程曲線放大圖，結果表明靜-動力人工邊界轉換法和黏彈性統一人工邊界法在外行波穿過邊界時有小幅值的震蕩，相比較而言，無限元統一人工邊界法在邊界處比較穩定。點B和點C的位移時程曲線也存在類似的現象。

3.2 半無限空間體受表面切向集中力和單位脈沖荷載聯合作用

將3.1節中的矩形均布靜力換為在坐標原點施加q=1×107N的切向集中力，其他不變。根據彈性力學的理論公式[20]和疊加原理，求得半無限空間體表面x軸方向上觀測點D(2,0,0)、E(4,0,0)、F(6,0,0)的x向位移理論解。同樣采用上述4種方法計算得到x軸方向上相應各點的x向位移數值解，如圖9～11。

圖9～11表明，模型在表面切向集中力和底面垂直向上輸入的x、y向脈沖荷載聯合作用下，無限元靜-動力統一人工邊界法的精度最高，即使是有限元-無限元耦合模型中的有限域范圍縮減25倍后，數值精度并沒有明顯的降低。雖然在切向集中力附近，無限元靜-動力統一人工邊界法的x向位移數值解與理論解之間有一定誤差，但是相比其他兩種人工邊界法，其x向位移數值解是最精確的，且距切向集中力施加點一定距離后，其x向位移數值解與理論解越接近。

當模型只受脈沖荷載作用時，黏彈性動力人工邊界法、黏彈性靜-動力統一人工邊界法、無限元靜-動力統一人工邊界法的精度非常接近且較高，動力響應的誤差對綜合響應的誤差貢獻較小，說明綜合響應的誤差大小在很大程度上取決于靜力效應的誤差大小。

4 分析某混凝土重力壩的靜-動力綜合效應

取混凝土重力壩的某一標準壩段建立三維有限元模型。該壩段最大壩高H=100 m,壩底寬b1=66.5 m,壩頂寬b2=7.0 m,下游壩面坡度m=0.7,上游壩面為鉛垂面。重力壩混凝土質量密度ρ=2 400 kg/m3,彈性模量E=24.0 GPa，泊松比ν=0.17；基巖質量密度ρ=2 600 kg/m3，彈性模量E=15.0 GPa，泊松比ν=0.25。上游正常蓄水位H0=95.0 m,下游正常蓄水位H1=9.7 m。規定坐標原點位于大壩上游面底部中心，x為順水流向，指向下游為正，y向為垂直水流向，指向左岸為正，z向為豎直方向，指向上為正。

分析大壩在自重、上下游靜水壓力、揚壓力以及地震作用下綜合響應，其中動水壓力在相應位置以附加質量的形式模擬。重力壩的兩側為自由邊界，忽略相臨壩段間的靜、動力相互作用。

對于靜力問題，在距結構一定遠處地基的位移已很小，遠域地基對結構的影響較小，可以忽略。通常情況，地基范圍至少取1.0倍結構高度[21]，因此靜-動力人工邊界轉換法(方法1)和黏彈性靜-動力統一人工邊界法(方法2)的有限元模型地基計算范圍為：深度方向自建基面向下延伸1.0倍壩高，上、下游方向和左、右岸方向同樣延伸1.0倍壩高，如圖12。無限元靜-動力統一人工邊界法(方法3)的模型：在有限單元的(單元類型C3D8R)外圍包裹一層無限單元(單元類型CIN3D8)，形成有限元-無限元耦合模型，以無限元模擬有限域地基外的無限域地基，如圖13。

另建一有限元-無限元耦合模型，其中有限地基范圍在深度、上下游和左右岸方向均取0.5倍壩高，地基范圍縮小了7/8(方法4)。

根據工程實際情況以及NB 35047-2015《水電工程水工建筑物抗震設計規范》[22]中的標準設計反應譜人工擬合地震動加速度時程。重力壩所在地為I0類場地，場地特征周期Ts=0.2s,標準設計反應譜最大值得代表值βmax=2.0，水平地震動峰值加速度PGA=0.1g,豎向地震動峰值加速度取水平向的2/3，加速度總時長t=20 s，時間步長為0.01 s。地基-大壩體系采用Rayleigh阻尼模型，按結構體系的前兩階自振頻率(f1=1.614、f2=2.521)確定阻尼常數，兩階阻尼比均為0.1，阻尼常數ɑ0、ɑ1分別為1.259和0.0077。利用SeismoSignal或Matlab等軟件通過數值積分求得各向速度、位移時程，即可將地震動輸入轉換為等效節點力加載到有限元模型和有限元-無限元耦合邊界面節點上。

圖14和15 給出了不同邊界法上游壩面壩頂中心點的x向和z向相對靜動力綜合位移時程，位移參考點為坐標原點，可見不同的邊界法對壩頂位置的順河向和豎向靜動綜合位移影響較小。圖16和17給出了上游壩面不同高度處的x向和z向相對位移最大值，經過歸一化可以看出，3種不同邊界法的z向最大相對位移的差別沿壩高基本不變，說明3種邊界法對于壩體的豎向靜動力綜合位移的計算精度是相近的。3種邊界法的x向最大相對位移的差別隨著壩高的增加逐漸減小，即對于水平向靜動綜合位移，3種邊界法的精度在壩高度較低的位置差別較大。根據以上半無限體的分析結論，3種邊界法的動力位移精度是相近的，從圖18同樣可以看出，3種邊界法的靜動力綜合效應的波型相似，只是均值水平上相差一個固定值，即由靜力效應引起的差別。造成這種結果的原因是前2種邊界對于靜力荷載為切向時的誤差較大。而無限元邊界方法則克服了這一缺點?？梢?，對于常常承受水平切向荷載的水工結構來說，水平向位移的精度是至關重要的，因此，對于水工結構的靜動力綜合問題，本文提出的無限元方法是適合的，另外2種方法較為繁瑣，精度也不理想，并且將無限元靜-動力統一人工邊界法的地基計算范圍縮小了7/8后，x向相對位移的誤差仍然小于黏彈性靜-動力統一人工邊界法和靜-動力人工邊界轉換法。

圖2有限元模型圖3有限元-無限元耦合模型

圖4入射位移波波形圖圖5模型頂部O點的z向位移解

圖6模型中部B點的z向位移解圖7模型底部C點z向位移解

圖8 O點0.55s之后z向位移解圖9模型頂部D點x向位移解

圖10模型頂部E點x向位移解圖11模型頂部F點x向位移解

圖12有限元模型圖13有限元-無限元耦合模型

圖14壩頂x方向相對位移時程曲線圖15壩頂z向相對位移時程曲線

圖16不同壩高位置歸一化x向相對位移圖17不同壩高位置歸一化z向相對位移

圖18 壩高h=6.0 m處x向相對位移時程曲線

5 結 論

本文根據黏彈性動力邊界理論，推導出有限元-無限元耦合界面上適用于外源波動輸入的等效節點力計算公式，基于Python語言，編制了等效節點力計算與施加為一體的程序，從而建立了基于ABAQUS無限元的靜-動力統一人工邊界。這種邊界在計算結構的綜合效應時，不需要進行靜、動力邊界轉換、應力場的導入以及節點約束反力的施加。通過對算例和實際工程的分析得出以下結論：

(1)ABAQUS無限元靜-動力統一人工邊界法實現相對簡單，其計算精度較靜、動力人工邊界轉換法和黏彈性靜-動力統一人工邊界法的精度要高，尤其當靜力作用方向為水平向時，其精度優勢能夠得到突顯。應用ABAQUS無限元靜-動力統一人工邊界條件分析重力壩-地基體系受靜-動力聯合作用的地震響應，驗證了該邊界對水工結構靜-動力分析的適用性。

(2)有限元-無限元耦合模型中的有限域地基范圍被大幅度縮減后，ABAQUS無限元靜-動力統一人工邊界法的計算精度仍然較高。表明該方法在提高計算效率的同時仍然能夠保持較高的精度，這對于大型非線性問題分析有良好的應用前景。