FLAC3D在土石壩地震反應(yīng)分析中的應(yīng)用

楊 星 ， 余 挺 ， 楊 貴

（1.河海大學(xué)巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2.河海大學(xué)巖土工程科學(xué)研究所，江蘇 南京 210098；3.中國(guó)水電顧問(wèn)集團(tuán)成都勘測(cè)設(shè)計(jì)研究院，四川 成都 610072）

0 引言

目前，土石壩地震反應(yīng)分析普遍采用的擬靜力法無(wú)法反映壩體的動(dòng)力特性、變形特性、輸入地震波的頻率及持續(xù)時(shí)間等特性，且得到的最小安全系數(shù)不能定量地反映大壩的抗震安全度?；诘葍r(jià)粘彈性模型的等效線性分析的有限單元法在參數(shù)確定和工程應(yīng)用方面積累了較豐富的經(jīng)驗(yàn)，得到了廣泛應(yīng)用。然而，該方法的缺點(diǎn)是加荷和卸荷的模量相同，不能直接計(jì)算土體在周期荷載作用下的永久變形；所用割線模量在小應(yīng)變時(shí)與非線性切線模量相近，但在大應(yīng)變時(shí)兩者相差很大，偏于不安全；塑性屈服模擬不合理以及不能計(jì)算大變形。FLAC3D(Fast Lagrangian Analysis of Continua)有限差分程序采用完全非線性動(dòng)力分析方法，可以遵循任何指定的非線性本構(gòu)模型；可以直接計(jì)算永久變形；可以分析巖土體在動(dòng)力荷載作用下發(fā)生的大變形，具有強(qiáng)大的后處理能力等[1]。

本文利用FLAC3D動(dòng)力彈塑性模型對(duì)土石壩進(jìn)行了地震反應(yīng)分析，著重討論了用FLAC3D進(jìn)行土石壩地震反應(yīng)分析時(shí)地震波濾波和基線校正、動(dòng)力邊界條件設(shè)置、瑞利阻尼參數(shù)確定等關(guān)鍵問(wèn)題，并用FISH語(yǔ)言實(shí)現(xiàn)了壩體初始剪切模量和體積模量隨平均主應(yīng)力的非線性變化。在此基礎(chǔ)上，對(duì)土石壩進(jìn)行了地震反應(yīng)分析，為利用FLAC3D進(jìn)行土石壩地震反應(yīng)分析提供了范例。

1 計(jì)算模型的建立

1.1 幾何模型和網(wǎng)格

某心墻堆石壩壩高100 m，壩頂寬10 m，上、下游壩坡坡比1∶2，心墻頂寬6 m，心墻上、下游坡比為1∶0.2。針對(duì)目前 FLAC3D前處理能力較弱，而ABAQUS有限元軟件具有較強(qiáng)建模能力的特點(diǎn)，土石壩模型在ABAQUS有限元軟件中建立。主體模型采用C3D8六面體單元，壩坡邊緣處采用三棱柱單元。運(yùn)用FORTRAN語(yǔ)言編寫了模型轉(zhuǎn)換程序，直接提取ABAQUS生成的數(shù)據(jù)文件（*.inp）中的節(jié)點(diǎn)（Node）、單元 （Element）和集合 （Set）信息， 轉(zhuǎn)化成為FALC3D相應(yīng)的節(jié)點(diǎn)（Gridpoint）、 單元（Zone）和組（Group）信息，并保存為網(wǎng)格數(shù)據(jù)文件（*.flac3d）。壩體幾何模型和轉(zhuǎn)換后的FLAC3D網(wǎng)格見圖1。計(jì)算中沒(méi)有考慮上游水荷載和壩體內(nèi)水的滲流作用。

圖1 壩體幾何模型和有限差分網(wǎng)格

1.2 本構(gòu)模型和材料參數(shù)

壩體堆石料采用摩爾-庫(kù)倫（Mohr-Coulomb）彈塑性本構(gòu)模型，當(dāng)剪應(yīng)力水平在彈性范圍內(nèi)時(shí)，通過(guò)粘滯阻尼達(dá)到能量耗散；當(dāng)剪應(yīng)力達(dá)到屈服強(qiáng)度時(shí)，通過(guò)塑性流動(dòng)實(shí)現(xiàn)能量耗散[2]。該模型能模擬土石壩地震時(shí)的塑性剪切變形，并且由于采用了彈塑性本構(gòu)模型，可以直接計(jì)算壩體的永久變形。壩體堆石料動(dòng)力計(jì)算的初始剪切模量隨平均有效主應(yīng)力的變化關(guān)系為

式中，C為動(dòng)剪切模量系數(shù)；n為動(dòng)剪切模量指數(shù)；兩者值由試驗(yàn)確定。σm為平均有效主應(yīng)力；Pa為大氣壓。

在FLAC3D中，根據(jù)式（1）通過(guò)FISH語(yǔ)言編程實(shí)現(xiàn)初始剪切模量隨平均有效主應(yīng)力的非線性變化。

堆石料的動(dòng)泊松比υd一般在0.3～0.4之間。動(dòng)泊松比υd取為0.35，則體積模量K=3G。計(jì)算中采用的材料參數(shù)見表1。壩體初始剪切模量計(jì)算結(jié)果見圖2。

1.3 地震波濾波和基線校正

模型動(dòng)力計(jì)算的地震波采用1989年美國(guó)加利福尼亞洛馬普列塔（Loma Prieta）地震期間列克星頓（Lexington）大壩的場(chǎng)地實(shí)測(cè)波，地震波峰值加速度（PGA）為0.17 g，持續(xù)時(shí)間約為40 s。輸入的地震波加速度時(shí)程曲線見圖3。本計(jì)算只考慮水平向地震波輸入。

表1 堆石料的計(jì)算參數(shù)

圖2 壩體初始剪切模量分布（單位：Pa）

圖3 輸入的地震波加速度時(shí)程曲線

模型的網(wǎng)格尺寸會(huì)影響地震波傳播的準(zhǔn)確性。庫(kù)勒邁耶（Kuhlemeyer）和賴斯默（Lysmer）的研究表明[3]，要想精確描述模型中地震波的傳播，模型網(wǎng)格尺寸 必須小于輸入波最高頻率對(duì)應(yīng)的波長(zhǎng)的1/10～1/8。對(duì)加速度做譜分析，其功率譜見圖4。由圖4分析可知，地震波95%以上的能量集中在5 Hz范圍內(nèi)。因此，地震波濾波的截?cái)囝l率確定為5 Hz。根據(jù)截?cái)囝l率和初始剪切模量，劃分網(wǎng)格時(shí)控制模型主體網(wǎng)格最大尺寸為5 m。

在FLAC3D地震計(jì)算中，若對(duì)時(shí)程加速度積分得到的最終速度和位移不為0，則在動(dòng)力計(jì)算結(jié)束后模型會(huì)出現(xiàn)繼續(xù)的殘余位移，此時(shí)需要對(duì)加速度時(shí)程曲線作基線校正。使用SeismoSignal地震波處理軟件對(duì)上述地震波進(jìn)行濾波和基線校正，修正后的速度和位移時(shí)程曲線分別見圖5、6。

1.4 動(dòng)力邊界條件設(shè)置

圖4 加速度時(shí)程的功率譜

圖5 修正后的速度時(shí)程曲線

圖6 未修正和修正后的位移時(shí)程曲線

在動(dòng)力計(jì)算中，模型邊界條件的選取是一項(xiàng)重要內(nèi)容，因?yàn)檫吔缟喜ǖ姆瓷鋾?huì)對(duì)動(dòng)力計(jì)算結(jié)果產(chǎn)生影響。FLAC3D提供了靜態(tài)（粘滯）邊界和自由場(chǎng)邊界兩種邊界條件。自由場(chǎng)邊界條件要求把模型設(shè)置得足夠遠(yuǎn)，模型的范圍取得足夠大，才能使模型邊界上波的反射盡可能的小，并且自由場(chǎng)邊界要求模型4個(gè)側(cè)面必須垂直。對(duì)于本文單獨(dú)以壩體為研究對(duì)象的模型，只能采用靜態(tài)（粘滯）邊界條件。

靜態(tài)（粘滯）邊界是Lysmer和Kuhlemeyer提出來(lái)的[4]，通過(guò)在模型邊界的法向和切向分別設(shè)置自由的阻尼器來(lái)實(shí)現(xiàn)吸收反射波。若在靜態(tài)（粘滯）邊界上施加動(dòng)荷載，則只能施加時(shí)程應(yīng)力。時(shí)程速度轉(zhuǎn)換為時(shí)程應(yīng)力的公式如下

式中，σn、σs分別為施加在靜態(tài)（粘滯）邊界上的法向應(yīng)力和切向應(yīng)力；ρ為密度；Cp和Cs為通過(guò)介質(zhì)的p波和s波的波速；υn和υs為邊界上的法向和切向速度分量。

1.5 阻尼的選取

FLAC3D動(dòng)力計(jì)算提供了瑞利阻尼、局部阻尼和滯后阻尼等3種阻尼形式。瑞利阻尼其理論與常規(guī)動(dòng)力分析方法類似，計(jì)算得到的加速度響應(yīng)規(guī)律比較符合實(shí)際[1]。因此，本文動(dòng)力計(jì)算采用瑞利阻尼。

在FALC3D中，使用瑞利阻尼必須確定最小中心頻率fmin和最小臨界阻尼比ξmin。瑞利阻尼是與頻率相關(guān)的，對(duì)于土石壩這種材料組成較為復(fù)雜的結(jié)構(gòu)，中心頻率既不是模型的自振頻率，也不是輸入地震波的主頻，而與兩者都有關(guān)系。本文采用對(duì)模型做無(wú)阻尼動(dòng)力計(jì)算的方法確定中心頻率。對(duì)土石壩進(jìn)行無(wú)阻尼地震動(dòng)力計(jì)算，得到壩頂點(diǎn)x方向的時(shí)程速度（見圖7）。對(duì)時(shí)程速度作傅里葉變換，得到該點(diǎn)時(shí)程速度的功率譜（見圖8）。確定該土石壩的中心頻率fmin約為1.05 Hz。對(duì)于巖土材料，臨界阻尼比的范圍一般是2%～5%，在FLAC3D動(dòng)力計(jì)算中一般可直接選取為5%[5]。

圖7 速度時(shí)程曲線

圖8 速度時(shí)程功率譜

2 計(jì)算結(jié)果分析

壩頂中心點(diǎn)的豎向永久變形見圖9。從圖9可知，壩體豎向永久變形達(dá)到22.3 cm。

輸入加速度和壩頂中心點(diǎn)地震反應(yīng)加速度的時(shí)程見圖10。從圖10可知，壩頂加速度最大值為5.03 m/s2，與輸入地震波的PGA相比，放大倍數(shù)達(dá)到2.9。在實(shí)際工程中，壩頂部位應(yīng)采取抗震加固措施，并預(yù)留足夠的安全超高，以保證大壩的安全。

圖9 壩頂豎向永久位移時(shí)程曲線

圖10 加速度時(shí)程曲線

3 結(jié)語(yǔ)

利用FLAC3D動(dòng)力彈塑性模型對(duì)土石壩進(jìn)行了地震反應(yīng)分析，結(jié)果表明，在地震荷載作用下，土石壩豎向永久變形達(dá)到22.3 cm，壩頂加速度最大值為5.03 m/s2，與輸入地震波的PGA相比，放大了2.9倍，計(jì)算結(jié)果符合已有工程研究規(guī)律。在實(shí)際工程中，壩頂部位應(yīng)采取抗震加固措施，并預(yù)留足夠的安全超高，以保證大壩的安全。

［1］陳育民，徐鼎平.FLAC/FLAC3D基礎(chǔ)與工程實(shí)例［M］.北京：中國(guó)水利水電出版社，2009.

［2］Noorzad R,Omidvar M.Seismic displacement analysis of embankment dams with reinforced cohesive shell［J］.Soil Dynamics and Earthquake Engineering,2010,30（11）:1149-1157.

［3］Kuhlemeyer R L,Lysmer J.Finite element method accuracy for wave propagation problems［J］.Journal of Soil Mechanics and Foundations,1973,99（SM4）:421-427.

［4］Lysmer J,Kuhlemyer R L.Finite dynamic model for infinite media［J］.Journal of Engineering Mechanics,1969,95（EM4）:859-877.

［5］Chakraborty D,Choudhury D.Investigation of the behavior of tailings earth dam under seismic conditions［J］.American Journal of Engineering and Applied Sciences,2009,2（3）:559-564.