考慮黏彈性邊界及縫面水壓的重力壩地震開裂模擬

徐 浩 江守燕 杜成斌

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

大壩作為擋水建筑物,在調(diào)節(jié)水資源、航運、發(fā)電等領(lǐng)域都發(fā)揮著舉足輕重的作用,因此如何準(zhǔn)確評估大壩抵御地震荷載的能力是水利水電工程抗震問題的研究重點[1].傳統(tǒng)有限元法模擬裂縫擴(kuò)展等不連續(xù)問題時需不停地進(jìn)行網(wǎng)格重劃分,Belytschko等[2]提出的擴(kuò)展有限元法(extended finite element methods,XFEM)可以使裂紋擴(kuò)展的幾何路徑獨立于網(wǎng)格,避免了網(wǎng)格重劃分,極大地提高了運算效率.杜效鵠等[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)問題,同時可以模擬截斷邊界上彈性恢復(fù)力和無限遠(yuǎn)場散射輻射,具備了較高的精度和穩(wěn)定性,黏彈性邊界最早由Deek和Randolph提出,隨后劉晶波和呂彥東[12]在此基礎(chǔ)上研究并提出了三維時域黏彈性邊界[13].地震荷載輸入方法在ABAQUS中有兩種不同的實現(xiàn)思路,苑舉衛(wèi)等[14]通過ABAQUS中的DLOAD和UTRACLOAD兩個子程序?qū)崿F(xiàn),而何建濤等[15]是計算每個節(jié)點在每一時刻對邊界的作用力,并將其作為振幅數(shù)據(jù)導(dǎo)入ABAQUS中,應(yīng)用到每個相應(yīng)的結(jié)點上.

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

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

1.1 位移模式

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

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

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

1.2 控制方程

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

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

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

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

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

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

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

建立以裂尖為坐標(biāo)原點的(r,θ)極坐標(biāo)系,則裂紋尖端附加函數(shù)Fl(x)由下列4個基函數(shù)組成:結(jié)果的穩(wěn)定性與時間步長無關(guān),時步長的選取依賴于求解精度,文中取δ=-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 黏彈性人工邊界模型及實現(xiàn)

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

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

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

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

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

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

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

根據(jù)波場分離,截斷邊界上的波可分解為散射場與自由場,在黏彈性邊界上設(shè)置自由場時自動滿足散射場,將外源地震波輸入到了有限域地基內(nèi).人工邊界上的等效結(jié)點荷載通過下式計算得到[23]

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

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

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

2.3 程序驗證

假設(shè)一個二維彈性矩形體,從底部邊界垂直入射地震波,則該入射波位移函數(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個四邊形等參元,底邊分布76個黏彈性單元,左、右側(cè)邊各分布38個黏彈性單元.二維彈性體網(wǎng)格圖及計算波形分別如圖2、3所示,計算時間增量Δt為0.005 s,模擬過程總時長t為2 s.

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

圖3 入射波位移時程曲線

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

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

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

3 重力壩地震開裂模擬

以Koyna重力壩為研究對象,圖5(a)為該重力壩的幾何模型及尺寸,大壩壩高103 m,壩頂寬度14.8 m,壩底寬70 m,地震時庫水位為91.75 m,上游面垂直,下游面折坡點以上部分坡比為1∶0.12,以下部分坡比為1∶0.76.壩體混凝土材料彈性模量為31.027 GPa,泊松比為0.15,質(zhì)量密度為2 463 kg/m3,抗拉強度為2.9 MPa,抗壓強度為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é)點平面等參元,共3 794個單元,3 986個結(jié)點,地基范圍向上游面、下游面以及深度方向延伸3倍壩高.

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

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

圖6 輸入的Koyna地震波

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

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

圖7 觀測點A的水平位移時程曲線

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

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

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

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

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

圖10 裂紋長度隨時間變化曲線

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

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

4 結(jié) 論

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

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

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