王 飛， 宋志強(qiáng)， 劉昱杰， 王 建

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

1 研究背景

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

對于一般符合小變形、線彈性的結(jié)構(gòu)，可分別獨立進(jìn)行靜力荷載和地震動力荷載作用效應(yīng)分析，將二者疊加即為結(jié)構(gòu)總作用效應(yīng)。對于具有材料或接觸非線性、強(qiáng)震作用下會發(fā)生大變形的結(jié)構(gòu)疊加原理原則上是不適用的，必須進(jìn)行靜、動力統(tǒng)一分析，由此引出靜、動力邊界的轉(zhuǎn)換或者統(tǒng)一問題。

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

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

高峰等[14]用黏彈性靜-動力統(tǒng)一人工邊界分析半無限空間體應(yīng)力，同樣認(rèn)為土體內(nèi)部各點的橫向應(yīng)力與理論解之間存在較大偏差。在土木工程結(jié)構(gòu)中，靜力效應(yīng)主要是結(jié)構(gòu)自重引起的，因此這一不足不會影響?zhàn)椥造o-動力統(tǒng)一人工邊界在該領(lǐng)域中的應(yīng)用。而在水工結(jié)構(gòu)中，除了結(jié)構(gòu)自重，靜水壓力、淤沙壓力等主要水平向荷載對結(jié)構(gòu)的靜力效應(yīng)是不容忽視的。

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

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

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

圖1 黏彈性邊界示意圖

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

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

(1)

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

(2)

與y軸平行時：

(3)

與z軸平行時：

(4)

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

對于三維模型具體表達(dá)式如下：

(1)對于底面

(5)

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

(2)對于x負(fù)向邊界面

(6)

(3)對于x正向邊界面

(7)

(4)對于y負(fù)向邊界面

(8)

(5)對于y正向邊界面

(9)

3 算例分析

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

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

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

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

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

(10)

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

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

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

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

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

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

4 分析某混凝土重力壩的靜-動力綜合效應(yīng)

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

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

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

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

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

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

圖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(yuǎn)=6.0 m處x向相對位移時程曲線

5 結(jié) 論

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

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

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