鄒 燁

(湖南省交通科學(xué)研究院有限公司,湖南 長(zhǎng)沙 410011)

巖土結(jié)構(gòu)如隧道、高邊坡、壩基等的動(dòng)力穩(wěn)定性一直是工程領(lǐng)域所關(guān)心的問(wèn)題[1],而利用數(shù)值方法進(jìn)行結(jié)構(gòu)動(dòng)力問(wèn)題分析的關(guān)鍵之一就是邊界條件的處理。巖土工程數(shù)值模擬過(guò)程中常使用人工邊界條件,以實(shí)現(xiàn)近場(chǎng)區(qū)域與無(wú)限外部區(qū)域之間波動(dòng)能量的交換[2],讓波動(dòng)自然傳至無(wú)限域內(nèi)。

人工邊界可分為兩種:全局人工邊界和局部人工邊界[3]。全局人工邊界如邊界元法,這類(lèi)邊界滿足了無(wú)限域內(nèi)的所有場(chǎng)方程和物理邊界條件,具有較高的精度;而由于其時(shí)空耦聯(lián)的特性(即某一邊界節(jié)點(diǎn)的運(yùn)動(dòng)與其他所有邊界節(jié)點(diǎn)之前各個(gè)時(shí)刻的運(yùn)動(dòng)有關(guān)),計(jì)算量與內(nèi)存需求較大,難以滿足大型工程問(wèn)題的分析需求;局部人工邊界是在一定的近似條件下導(dǎo)出的人工邊界條件,目前常用的粘性邊界、粘彈性邊界、透射邊界、無(wú)限元邊界均屬于局部人工邊界[4]。局部人工邊界具有時(shí)空解耦的特性(即邊界上某一點(diǎn)的運(yùn)動(dòng)只與相鄰邊界節(jié)點(diǎn)與該節(jié)點(diǎn)上一時(shí)刻的運(yùn)動(dòng)情況有關(guān)),在犧牲一部分精度的條件下具有節(jié)省機(jī)時(shí)與內(nèi)存的作用,在許多大型工程中得到了應(yīng)用。

無(wú)限元是一種幾何上趨于無(wú)窮的單元[5],ABAQUS中的無(wú)限元單元已有許多相關(guān)的研究與應(yīng)用[6-12]。針對(duì)三維情況,僅有的八節(jié)點(diǎn)、十二節(jié)點(diǎn)與十六節(jié)點(diǎn)的六面體三維無(wú)限元單元要求有限域與無(wú)限域連接面為四邊形,這一要求對(duì)于實(shí)際工程中復(fù)雜的地質(zhì)模型比較嚴(yán)苛,因?yàn)閺?fù)雜巖土工程模型為減小計(jì)算成本,獲取更好的收斂性,往往采用四面體單元、六面體無(wú)限元單元的使用由此受到了限制。從這一問(wèn)題出發(fā),在ABAQUS中三維無(wú)限元的基礎(chǔ)上進(jìn)行改進(jìn),提出一種由六面體無(wú)限元單元演變而來(lái)的退化無(wú)限元單元以適應(yīng)有限域與無(wú)限域交界面為三角形的情況,推導(dǎo)退化無(wú)限元單元的等效節(jié)點(diǎn)應(yīng)力計(jì)算公式,并研究其地震動(dòng)輸入方式。通過(guò)建立有限元-無(wú)限元(FEM-INF)耦合分析模型進(jìn)行數(shù)值試驗(yàn),與理論結(jié)果進(jìn)行對(duì)比驗(yàn)證退化無(wú)限元單元的合理性。

1 動(dòng)力無(wú)限元人工邊界

1.1 邊界阻尼

無(wú)限元由Bettess和Zienkiewicz于1977年提出[13],用于解決有限元方法求解無(wú)限域問(wèn)題時(shí)的局限性。ABAQUS中的無(wú)限元理論參考了Lysmer與Kuhlemeyer等的工作,其動(dòng)力無(wú)限元可視為一種黏性吸收邊界[14],通過(guò)引入邊界阻尼d來(lái)消除反射。

w=f1(z-ct)+f2(z+ct)

(1)

式中:z為軸向坐標(biāo);c為波速;t為時(shí)間;w為波動(dòng)所引起的邊界節(jié)點(diǎn)位移,對(duì)于任意瞬時(shí)t,w為z的函數(shù);f1(z-ct)、f2(z+ct)分別為入射波與反射波位移方程。阻尼力的大小可以表示為:

(2)

(3)

式中:G為傳播介質(zhì)的剪切模量;λ為拉梅常數(shù);θ為體應(yīng)變。

由邊界節(jié)點(diǎn)的平衡條件有σ=σd,即:

(4)

(5)

(6)

同樣按照‘邊界上不產(chǎn)生反射’的要求,可得切向阻尼系數(shù):

(7)

1.2 退化的動(dòng)力無(wú)限元單元

ABAQUS提供了三種適用于連續(xù)實(shí)體介質(zhì)的三維無(wú)限元單元,均為六面體單元,分別為八節(jié)點(diǎn)、十二節(jié)點(diǎn)和十八節(jié)點(diǎn)三種類(lèi)型(見(jiàn)圖1),單元形狀決定了使用時(shí)要求有限元與無(wú)限元的交界面為四邊形。對(duì)一些模型復(fù)雜的實(shí)際工程問(wèn)題,網(wǎng)格劃分時(shí)常采用適用性更強(qiáng)、計(jì)算規(guī)模更低的四面體單元,這時(shí)將不能使用僅提供的三維六面體無(wú)限元單元。

為解決這一問(wèn)題,在現(xiàn)有單元的基礎(chǔ)上提出了一種退化的無(wú)限元單元。具體做法是,讓單元交界面上的任意一條邊上的節(jié)點(diǎn)重合而退化成一個(gè)點(diǎn),使交界面由(曲邊)四邊形退化成(曲邊)三角形。交界面的退化過(guò)程見(jiàn)圖2。

圖2 交界面的退化過(guò)程

2 退化單元的等效節(jié)點(diǎn)力

(8)

2.1 等效節(jié)點(diǎn)力的計(jì)算

以地震波從底部垂直入射的情況為例,見(jiàn)圖3。

圖3 等效節(jié)點(diǎn)力計(jì)算

(1)在P波作用下,假設(shè)輸入波的速度時(shí)程為Vp(t)。

(9)

(10)

由對(duì)稱性,右側(cè)、前側(cè)、后側(cè)邊界的法向與切向等效節(jié)點(diǎn)力分別為式(11)、式(12)、式(13):

(11)

(12)

(13)

(2)在S波作用下,假設(shè)輸入波的速度時(shí)程為Vs(t)。

(14)

(15)

由反對(duì)稱條件,右側(cè)邊界:

(16)

(17)

由對(duì)稱條件,后側(cè)邊界:

(18)

各節(jié)點(diǎn)力的方向見(jiàn)圖3。

上列各式中,Vp-1(t)、Vp-2(t)、Vs-1(t)、Vs-2(t)為相應(yīng)節(jié)點(diǎn)在t時(shí)刻的振動(dòng)速度,需要由入射波時(shí)程與反射波時(shí)程在計(jì)入對(duì)應(yīng)的延時(shí)后疊加求得。其中:

(19)

式中:L為底面至頂面的距離,l為節(jié)點(diǎn)到底面的距離;當(dāng)頂面非平面時(shí),L取節(jié)點(diǎn)i在底面與頂面投影點(diǎn)之間的距離。

2.2 退化單元的等效節(jié)點(diǎn)力

退化單元等效節(jié)點(diǎn)力的計(jì)算原理與2.1節(jié)所示一樣,在計(jì)算因退化而重合的節(jié)點(diǎn)的等效力時(shí),把這些點(diǎn)仍然看成單獨(dú)的點(diǎn)。以退化的八節(jié)點(diǎn)三維無(wú)限元單元為例,如圖4所示。

圖4 單元的轉(zhuǎn)換

重合點(diǎn)的影響面積為:

(20)

非重合點(diǎn)的影響面積為:

(21)

式中:n為受節(jié)點(diǎn)i影響的單元的個(gè)數(shù);Ak為單元的面積。

各邊界上的等效節(jié)點(diǎn)力計(jì)算公式見(jiàn)表1。

表1 退化單元等效節(jié)點(diǎn)力計(jì)算公式

當(dāng)節(jié)點(diǎn)為重合點(diǎn)時(shí),q=2;當(dāng)節(jié)點(diǎn)不是重合點(diǎn)時(shí),q=4。

3 算例驗(yàn)證

為說(shuō)明退化三維動(dòng)力無(wú)限元單元的合理性,下面用數(shù)值試驗(yàn)加以驗(yàn)證。建立兩個(gè)100×100×100的立方塊,將這兩個(gè)立方塊分別劃分為八節(jié)點(diǎn)六面體模型(C3D8)與四節(jié)點(diǎn)四面體模型(C3D4),以分別模擬無(wú)限域采用常規(guī)三維無(wú)限元單元與退化無(wú)限元單元的情況。立方塊體模型除頂面外,其余5個(gè)面均采用無(wú)限元邊界,無(wú)限元邊界分別采用常規(guī)的八節(jié)點(diǎn)無(wú)限元單元與退化的八節(jié)點(diǎn)無(wú)限元單元。圖5為有限元區(qū)域劃分為六面體單元C3D8時(shí)的FEM-INF耦合模型;圖6為有限元區(qū)域劃分四面體單元C3D4時(shí)的FEM-INF耦合模型。邊界上的節(jié)點(diǎn)一般較多,等效節(jié)點(diǎn)力的計(jì)算與輸入要通過(guò)自編程序?qū)崿F(xiàn)。

圖5 有限域?yàn)榱骟w單元時(shí)的FEM-INF耦合模型

圖6 有限域?yàn)樗拿骟w單元時(shí)的FEM-INF耦合模型

以P波的形式從模型底部邊界垂直輸入圖7所示的波形,這是一個(gè)峰值加速度為4 m/s2的單峰波,其加速度為:

圖7 入射波波形

A=4e-30π(t-0.5)2

(22)

模型的介質(zhì)材料參數(shù)設(shè)置為:彈性模量E=50 MPa,泊松比ν=0.3;重度γ=25 kN/m3。模型除受重力與地震力作用外,無(wú)其他形式的力作用。

根據(jù)波速的計(jì)算公式,可得縱波波速為:cp=164.08 m/s;橫波波速為:cs=87.71 m/s。因此,P波從模型底面?zhèn)鞑ブ另斆娴臅r(shí)間延遲約為0.61 s,S波從模型底面?zhèn)鞑ブ另斆娴臅r(shí)間延遲約為1.14 s。

考慮P波從底面垂直入射的情況,算例的計(jì)算結(jié)果見(jiàn)圖8(a)。若不計(jì)立方塊體材料介質(zhì)的阻尼,按照理論計(jì)算的結(jié)果,P波從立方塊體底部入射后,經(jīng)過(guò)一定的時(shí)間延遲后到達(dá)頂面,由于這是一個(gè)與空氣的界面,因此,P波將會(huì)在頂面發(fā)生反射,入射波與反射波疊加會(huì)在頂面出現(xiàn)一個(gè)2倍的放大效應(yīng),即頂面P波反射加速度的理論解為:

圖8 加速度響應(yīng)

Ap-頂=8e-30π(t-1.11)2

(23)

反射波在此傳回地面后,在無(wú)限元邊界的作用下,不再發(fā)生反射,傳向了無(wú)窮遠(yuǎn)處,這也就是人工邊界條件的作用,避免了波動(dòng)多次反射現(xiàn)象的發(fā)生。從圖8(a)可以看出,ABAQUS中的常規(guī)三維無(wú)限元單元在動(dòng)力分析中取得了與理論解幾乎一致的效果,而采用退化的三維無(wú)限元單元,雖在后續(xù)有小幅的波動(dòng)響應(yīng),但其結(jié)果與理論解仍具有很好的吻合度。

S波從底部垂直入射的情況采用與P波入射時(shí)類(lèi)似的方式進(jìn)行。頂面S波反射加速度的理論解為:

As-頂=8e-30π(t-1.64)2

(24)

算例驗(yàn)證結(jié)果見(jiàn)圖8(b)。從圖可以看出,采用常規(guī)無(wú)限元單元與退化無(wú)限元單元對(duì)地震波均具有良好的吸收效果;采用退化無(wú)限元單元雖精度要求略低于常規(guī)無(wú)限元單元,當(dāng)反射波離開(kāi)頂面后,有小幅的波動(dòng)現(xiàn)象發(fā)生,但很快趨于穩(wěn)定,能基本滿足動(dòng)力計(jì)算上的要求。

4 結(jié) 論

在ABAQUS中無(wú)限元單元的基礎(chǔ)上,提出了一種退化無(wú)限元單元,并利用數(shù)值試驗(yàn)進(jìn)行驗(yàn)證,可以得到以下結(jié)論:

(1) 退化無(wú)限元單元能使有限域與無(wú)限域之間的交界面由四邊形退化為三角形,對(duì)有限域?yàn)樗拿骟w單元的情況也能較好適用,提升了三維無(wú)限元單元對(duì)主體網(wǎng)格剖分方式的適應(yīng)能力,且無(wú)需進(jìn)行二次開(kāi)發(fā),方便工程應(yīng)用。

(2) 退化單元數(shù)值模擬結(jié)果與理論解吻合良好,表明退化單元具有較好的精度,研究?jī)?nèi)容為構(gòu)建復(fù)雜的FEM-INF耦合分析模型提供了便利,可應(yīng)用于工程動(dòng)力問(wèn)題的分析。

(3) ABAQUS中采用動(dòng)力無(wú)限元邊界條件時(shí),需要將地震動(dòng)以等效節(jié)點(diǎn)力的形式進(jìn)行輸入。邊界節(jié)點(diǎn)較多時(shí),需要通過(guò)編程實(shí)現(xiàn)等效節(jié)點(diǎn)力的自動(dòng)施加。

值得指出的是,研究?jī)?nèi)容雖然采用ABAQUS進(jìn)行,但其中單元退化過(guò)程的處理方法對(duì)其他類(lèi)似問(wèn)題仍具有較好的參考意義。