徐 浩 江守燕 杜成斌

(河海大學(xué) 工程力學(xué)系,南京 210098)

大壩作為擋水建筑物,在調(diào)節(jié)水資源、航運(yùn)、發(fā)電等領(lǐng)域都發(fā)揮著舉足輕重的作用,因此如何準(zhǔn)確評估大壩抵御地震荷載的能力是水利水電工程抗震問題的研究重點(diǎn)[1].傳統(tǒng)有限元法模擬裂縫擴(kuò)展等不連續(xù)問題時(shí)需不停地進(jìn)行網(wǎng)格重劃分,Belytschko等[2]提出的擴(kuò)展有限元法(extended finite element methods,XFEM)可以使裂紋擴(kuò)展的幾何路徑獨(dú)立于網(wǎng)格,避免了網(wǎng)格重劃分,極大地提高了運(yùn)算效率.杜效鵠等[3]結(jié)合擴(kuò)展有限元法對重力壩動態(tài)響應(yīng)進(jìn)行了數(shù)值模擬;董玉文等[4]推導(dǎo)建立了滿足大壩水力壓裂問題的擴(kuò)展有限元法;高景泉等[5]研究了靜力作用下重力壩踵預(yù)設(shè)裂縫的水力壓裂問題;龐林等[6]基于比例邊界有限元法研究了縫水壓力和裂縫面接觸條件對重力壩裂縫擴(kuò)展的影響;李宗利等[7]研究了混凝土重力壩裂縫面上施加不同水壓分布形式的穩(wěn)定性情況,這些研究均為地震作用下壩體水力劈裂研究奠定了重要基礎(chǔ).

重力壩地基邊界模擬和地震荷載輸入方式是大壩地震響應(yīng)分析的關(guān)鍵,對數(shù)值模擬結(jié)果的準(zhǔn)確性有著直接的影響.目前,動力人工邊界主要有透射邊界[8]、無限元邊界[9]、黏性邊界[10]、黏彈性邊界[11]等.透射邊界應(yīng)用較復(fù)雜,且通常只限于二階精度以內(nèi),可能存在高頻失穩(wěn)的問題;黏性人工邊界忽略了彈性恢復(fù)力,只考慮阻尼吸收作用,容易出現(xiàn)結(jié)構(gòu)整體漂移,不宜應(yīng)用在較復(fù)雜的大型結(jié)構(gòu)中.黏彈性人工邊界較于黏性邊界避免了低頻失穩(wěn)問題,同時(shí)可以模擬截?cái)噙吔缟蠌椥曰謴?fù)力和無限遠(yuǎn)場散射輻射,具備了較高的精度和穩(wěn)定性,黏彈性邊界最早由Deek和Randolph提出,隨后劉晶波和呂彥東[12]在此基礎(chǔ)上研究并提出了三維時(shí)域黏彈性邊界[13].地震荷載輸入方法在ABAQUS中有兩種不同的實(shí)現(xiàn)思路,苑舉衛(wèi)等[14]通過ABAQUS中的DLOAD和UTRACLOAD兩個(gè)子程序?qū)崿F(xiàn),而何建濤等[15]是計(jì)算每個(gè)節(jié)點(diǎn)在每一時(shí)刻對邊界的作用力,并將其作為振幅數(shù)據(jù)導(dǎo)入ABAQUS中,應(yīng)用到每個(gè)相應(yīng)的結(jié)點(diǎn)上.

本文采用擴(kuò)展有限元法仿真模擬混凝土重力壩的地震開裂過程,根據(jù)波場分離理論和與黏彈性人工邊界相適應(yīng)的地震動輸入公式,結(jié)合ABAQUS二次開發(fā)功能,利用Matlab編寫相應(yīng)的邊界條件自動施加和地震動輸入程序,分析探討不同的抗震分析模型對大壩地震開裂過程的影響.

1 動力擴(kuò)展有限元法

1.1 位移模式

對于裂紋問題,XFEM位移模式可以表示為:

式中:Ni(x)為常規(guī)有限元插值形函數(shù);ui為常規(guī)結(jié)點(diǎn)自由度;ai為與跳躍函數(shù)有關(guān)的附加自由度;bli為裂尖附加自由度;Ω為所有離散節(jié)點(diǎn)的集合;ΩΓ為裂紋貫穿單元的結(jié)點(diǎn)集;ΩΛ為裂尖單元的結(jié)點(diǎn)集.H(x)為Heaviside函數(shù),是反映裂紋面位移不連續(xù)的一個(gè)附加函數(shù),分別取裂紋面上、下側(cè)的+1和-1來反映位移不連續(xù)性,定義為

式中:r、θ分別為裂尖局部極坐標(biāo)半徑、夾角.

1.2 控制方程

因縫水壓力為結(jié)構(gòu)所受外荷載,故作為面力施加到裂縫面,且作用于增強(qiáng)自由度上.因此參照文獻(xiàn)[16],可以假定縫面水壓力P的分布形式為:

式中:r為開裂縫面上節(jié)點(diǎn)到裂紋尖端的距離;a為裂縫長度;P0為裂紋開口處的靜水壓力;論文采用均勻分布水壓,故n取0.

考慮阻尼的影響,由彈性動力學(xué)的基本方程[17]及XFEM位移模式,運(yùn)用虛功原理[18],可得動力XFEM的控制方程為:

采用Newmark隱式時(shí)間積分算法進(jìn)行時(shí)間離散[19],假設(shè)

式中:α和β為積分常數(shù),與積分精度和數(shù)值穩(wěn)定性有關(guān),α=0.25(1-δ)2,β=0.5-δ,且δ∈[-1/3,0].

Newmark隱式時(shí)間積分算法是無條件穩(wěn)定的,

式中:x為觀測點(diǎn);x*為觀測點(diǎn)在裂紋面上的投影點(diǎn);n為x*處裂紋面的單位外法向.

建立以裂尖為坐標(biāo)原點(diǎn)的(r,θ)極坐標(biāo)系,則裂紋尖端附加函數(shù)Fl(x)由下列4個(gè)基函數(shù)組成:結(jié)果的穩(wěn)定性與時(shí)間步長無關(guān),時(shí)步長的選取依賴于求解精度,文中取δ=-0.05.

1.3 裂紋開裂準(zhǔn)則

ABAQUS中基于拉伸分離損傷法的初始損傷準(zhǔn)則有6種,但只有最大主應(yīng)力準(zhǔn)則和最大主應(yīng)變準(zhǔn)則可用于模擬裂紋擴(kuò)展等不連續(xù)問題.選擇最大主應(yīng)力準(zhǔn)則[20]作為裂紋開裂的判據(jù),可表示為

2 黏彈性人工邊界及地震動輸入模型

2.1 黏彈性人工邊界模型及實(shí)現(xiàn)

通過在邊界上設(shè)置連續(xù)分布的并聯(lián)彈簧阻尼單元來實(shí)現(xiàn)黏彈性人工邊界(如圖1所示).

圖1 黏彈性邊界彈簧-阻尼器單元

其中彈簧和阻尼器的力學(xué)參數(shù)由圍巖材料參數(shù)所決定,可通過下式計(jì)算[21]

式中:KBT、KBN分別為彈簧在切向和法向的剛度系數(shù);CBT、CBN分別為阻尼器在切向和法向的阻尼系數(shù);αT、αN為修正系數(shù),一般情況下,αT、αN的取值參考劉晶波等[22]的研究,取αT=0.5,αN=1.0;R為邊界點(diǎn)與波源的距離;ρ、G分別為材料的密度和剪切模量;Cs和Cp分別為材料的剪切波速和縱波波速,可表示為

式中:E、λ、μ分別為介質(zhì)的彈性模量、Lame常數(shù)和泊松比.

利用有限元軟件ABAQUS中提供的彈簧單元和阻尼器單元,可以實(shí)現(xiàn)截?cái)噙吔缟细鹘Y(jié)點(diǎn)彈簧單元和阻尼單元的施加.文中基于MATLAB語言編寫了批量自動施加黏彈性邊界的程序,首先讀取邊界結(jié)點(diǎn)坐標(biāo),然后根據(jù)式(8)計(jì)算出彈簧剛度系數(shù)和阻尼系數(shù),最后通過修改ABAQUS中的輸入數(shù)據(jù)文件完成所有邊界結(jié)點(diǎn)處彈簧和阻尼單元的施加.

2.2 地震動輸入模型及實(shí)現(xiàn)

根據(jù)波場分離,截?cái)噙吔缟系牟煞纸鉃樯⑸鋱雠c自由場,在黏彈性邊界上設(shè)置自由場時(shí)自動滿足散射場,將外源地震波輸入到了有限域地基內(nèi).人工邊界上的等效結(jié)點(diǎn)荷載通過下式計(jì)算得到[23]

式中:KB為剛度系數(shù);CB為阻尼系數(shù)為該人工邊界點(diǎn)的自由場位移為該點(diǎn)的自由場速度為該點(diǎn)應(yīng)力;n為該點(diǎn)外法線方向余弦.

根據(jù)彈性力學(xué)理論,結(jié)合底邊位移邊界條件和波動理論,首先計(jì)算應(yīng)變,然后由應(yīng)力-應(yīng)變關(guān)系求出應(yīng)力,再由邊界結(jié)點(diǎn)平衡關(guān)系得到邊界結(jié)點(diǎn)上的等效結(jié)點(diǎn)力,具體邊界結(jié)點(diǎn)力計(jì)算公式可參見文獻(xiàn)[24].

文中采用的方法是提前計(jì)算人工邊界上各點(diǎn)隨時(shí)間變化的切向與法向荷載數(shù)值,利用ABAQUS軟件中的指令在邊界點(diǎn)上輸入隨時(shí)間變化的等效荷載.具體做法是:首先提取壩體模型的尺寸大小、材料參數(shù)值、單元結(jié)點(diǎn)坐標(biāo)和編號、結(jié)點(diǎn)控制面積大小、地震波的位移和速度時(shí)程曲線;然后利用式(8)中的彈簧剛度系數(shù)和阻尼系數(shù)計(jì)算值,結(jié)合地震波處理軟件SeismoSignal分時(shí)段積分得到的速度與位移時(shí)程,推導(dǎo)出地震在自由場產(chǎn)生的應(yīng)力以及彈簧和阻尼單元產(chǎn)生的附加應(yīng)力;最后與提取得到的結(jié)點(diǎn)控制面積相乘,得到邊界上各結(jié)點(diǎn)的等效荷載時(shí)程情況,通過在原有模型文件中插入新的inp計(jì)算文件,在ABAQUS求解器提交運(yùn)算即可.

2.3 程序驗(yàn)證

假設(shè)一個(gè)二維彈性矩形體,從底部邊界垂直入射地震波,則該入射波位移函數(shù)為:

彈性體長為762 m,寬為381 m,質(zhì)量密度ρ=2.7 kg/m3,彈性模量E=13.23 MPa,泊松比ν=0.25,P波波速cP=2 425 m/s,SV波波速csv=1 400 m/s.假定為平面應(yīng)變問題,彈性體被離散成76×38個(gè)四邊形等參元,底邊分布76個(gè)黏彈性單元,左、右側(cè)邊各分布38個(gè)黏彈性單元.二維彈性體網(wǎng)格圖及計(jì)算波形分別如圖2、3所示,計(jì)算時(shí)間增量Δt為0.005 s,模擬過程總時(shí)長t為2 s.

圖2 二維彈性體網(wǎng)格圖

圖3 入射波位移時(shí)程曲線

數(shù)值計(jì)算結(jié)果如圖4所示,上邊界A、B、C3個(gè)觀測點(diǎn)位移幅值接近于入射波位移幅值的兩倍,這一現(xiàn)象證實(shí)了入射波傳播至上邊界后形成反射波,波形疊加形成了上邊界各點(diǎn)的位移時(shí)程曲線,0.75 s后3個(gè)觀測點(diǎn)都不再有位移變化,這表明反射波抵達(dá)底部邊界后可以被黏彈性邊界充分吸收,之后不再反射.

圖4 理論解與數(shù)值解的比較

相比于上邊界各點(diǎn)波形疊加只形成了單個(gè)波峰,底部邊界上D、E兩點(diǎn)均有一前一后兩個(gè)波峰,可以看出入射波在形成第一個(gè)波峰之后,經(jīng)上邊界反射之后在下邊界又形成了一個(gè)波峰.另外,同一水平位置觀測點(diǎn)的位移曲線規(guī)律基本相同,說明波的傳播是均勻的.該數(shù)值結(jié)果與波傳播規(guī)律的理論解基本一致,說明文中所應(yīng)用的程序與模擬方法是正確可行的.

3 重力壩地震開裂模擬

以Koyna重力壩為研究對象,圖5(a)為該重力壩的幾何模型及尺寸,大壩壩高103 m,壩頂寬度14.8 m,壩底寬70 m,地震時(shí)庫水位為91.75 m,上游面垂直,下游面折坡點(diǎn)以上部分坡比為1∶0.12,以下部分坡比為1∶0.76.壩體混凝土材料彈性模量為31.027 GPa,泊松比為0.15,質(zhì)量密度為2 463 kg/m3,抗拉強(qiáng)度為2.9 MPa,抗壓強(qiáng)度為24.1 MPa,斷裂能為250 N/m.壩基巖體的彈性模量為24 GPa,泊松比0.25,質(zhì)量密度為2 500 kg/m3.參照文獻(xiàn)[17],為便于研究和對比,假定壩體上游面距離壩踵20 m處有一水平預(yù)設(shè)裂紋,長度設(shè)為8 m.圖5(b)為壩體-地基系統(tǒng)網(wǎng)格,壩基系統(tǒng)離散成四結(jié)點(diǎn)平面等參元,共3 794個(gè)單元,3 986個(gè)結(jié)點(diǎn),地基范圍向上游面、下游面以及深度方向延伸3倍壩高.

圖5 重力壩壩體幾何模型及壩基系統(tǒng)網(wǎng)格

在進(jìn)行數(shù)值計(jì)算時(shí),施加的荷載包括壩體自重、上游面靜水壓力、地震荷載以及動水壓力.地震波采用Koyna地震實(shí)測地震波,水平向峰值加速度為0.49g,豎向峰值加速度為0.34g,如圖6所示.動水壓力采用附加質(zhì)量法近似模擬.計(jì)算總時(shí)間為10 s,計(jì)算時(shí)間增量為0.01 s.根據(jù)加速度時(shí)程曲線分時(shí)段積分得到速度與位移時(shí)程,再根據(jù)底邊與側(cè)邊各點(diǎn)速度和位移計(jì)算出等效荷載實(shí)現(xiàn)地震波輸入.

圖6 輸入的Koyna地震波

文中分析探討了3種模型對大壩地震開裂的影響,模型1為無質(zhì)量地基模型且縫面無水壓;模型2為無質(zhì)量地基模型縫面有水壓;模型3考慮縫面有水壓的同時(shí),地基采用有質(zhì)量彈性地基并設(shè)置黏彈性人工邊界,地震波輸入采用等效荷載形式.采用黏彈性人工邊界時(shí),底邊與側(cè)邊共有黏彈性單元83個(gè).

圖7給出了壩頂觀測點(diǎn)A在3種模型下的水平位移時(shí)程曲線,可以看出縫面有無水壓兩種情況的位移變化規(guī)律大致一樣,有水壓的位移幅值比無水壓的位移幅值高出26.1%,且二者出現(xiàn)位移幅值的時(shí)間與預(yù)設(shè)裂紋長度驟然增加的時(shí)間相吻合,這表明無質(zhì)量地基模型在8～9 s之間地震響應(yīng)最明顯.考慮了黏彈性邊界的位移幅值較不考慮黏彈性的位移幅值顯著降低,說明采用無質(zhì)量剛性地基重力壩抗震模型相較于黏彈性地基模型,放大了地震動響應(yīng).

圖7 觀測點(diǎn)A的水平位移時(shí)程曲線

圖8給出了3種模型下的裂紋擴(kuò)展路徑,圖9對比了3種模型下的裂紋開裂擴(kuò)展路徑.

圖8 地震開裂擴(kuò)展路徑

圖9 3種模型下裂紋開裂擴(kuò)展路徑對比

從圖中不難看出裂紋面上有水壓時(shí)的偏轉(zhuǎn)角要小于無水壓時(shí)的結(jié)果,這個(gè)結(jié)果與文獻(xiàn)[24]中得到的結(jié)果相一致,而考慮黏彈性邊界且采用有質(zhì)量彈性地基時(shí)的裂紋擴(kuò)展路徑則比較平緩,偏轉(zhuǎn)角也要小于不考慮黏彈性邊界時(shí)的結(jié)果.分析裂紋路徑成因,可能由于慣性力周期性變化,因此在開裂初期不斷由上游面向下游面擴(kuò)展,在裂紋擴(kuò)展前期,自重引起的壓力作用導(dǎo)致裂紋閉合,故裂紋沿水平方向延伸,保持了壩體的基本穩(wěn)定.而在后期裂紋在擴(kuò)展中開裂路徑并非水平延伸,而是斜向下延伸,推測原因是裂紋尖端區(qū)應(yīng)力需受慣性力和壩體自重疊加影響,而在裂紋擴(kuò)展后期裂紋尖端區(qū)應(yīng)力受慣性力影響減弱,壩體自重影響增大,故裂紋向右下方延伸.

圖10對比了3種模型下的裂紋開裂擴(kuò)展路徑,給出了裂紋長度隨時(shí)間的變化曲線,容易看出縫面有水壓時(shí)裂紋的開裂距離大于無水壓的開裂距離,分析其原因,是由于外法向水壓作用于裂紋表面,相當(dāng)于增加了Ⅰ型裂紋開裂的驅(qū)動力.此外考慮黏彈性邊界時(shí)的裂紋開裂長度小于無質(zhì)量地基模型的結(jié)果.

圖10 裂紋長度隨時(shí)間變化曲線

圖11是3種情況的裂縫寬度變化情況,不難發(fā)現(xiàn)三者的擴(kuò)展規(guī)律相一致,都是裂紋寬度越接近裂尖位置越小,縫面有水壓時(shí)裂紋擴(kuò)展寬度明顯大于無水壓情況,這是因?yàn)榱鸭y面上施加了外法線方向的水壓力,水壓的驅(qū)動力必然引起法線方向的擴(kuò)展寬度變大.考慮黏彈性邊界時(shí)的最大裂紋寬度小于無質(zhì)量地基模型的結(jié)果,減小了約19%.

圖11 3種模型下裂縫寬度變化對比

4 結(jié) 論

采用擴(kuò)展有限元法仿真模擬了混凝土重力壩的地震開裂過程.分析探討了3種模型對大壩地震開裂過程的影響:1)無質(zhì)量地基模型且縫面無水壓;2)無質(zhì)量地基模型且縫面有水壓;3)縫面有水壓且考慮黏彈性人工邊界有質(zhì)量地基及地震動等效荷載輸入的抗震計(jì)算模型.

仿真模擬結(jié)果表明:1)黏彈性人工邊界模型計(jì)算得到的位移幅值較無質(zhì)量地基模型的結(jié)果小,即無質(zhì)量剛性地基抗震分析模型放大了大壩的地震動響應(yīng);2)縫面有水壓時(shí)的裂紋偏轉(zhuǎn)角小于無水壓時(shí)的結(jié)果,而考慮黏彈性邊界且采用有質(zhì)量彈性地基時(shí)的裂紋擴(kuò)展路徑則更平緩;3)縫面有水壓時(shí)裂紋的開裂距離大于無水壓時(shí)的開裂距離,考慮黏彈性人工邊界時(shí)的裂紋擴(kuò)展長度小于無質(zhì)量地基模型的結(jié)果;4)縫面有水壓時(shí)裂紋擴(kuò)展寬度明顯大于無水壓情況,考慮黏彈性邊界時(shí)的最大裂紋寬度小于無質(zhì)量地基模型的結(jié)果,減小了約19%.

此外,庫水的可壓縮性對重力壩地震開裂也有一定影響,在今后的研究中將進(jìn)一步考慮大壩-地基的流固耦合效應(yīng),探討無限庫水的影響.