FLAC3D在土石壩地震反應分析中的應用

楊 星 ， 余 挺 ， 楊 貴

（1.河海大學巖土力學與堤壩工程教育部重點實驗室，江蘇 南京 210098；2.河海大學巖土工程科學研究所，江蘇 南京 210098；3.中國水電顧問集團成都勘測設計研究院，四川 成都 610072）

0 引言

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

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

1 計算模型的建立

1.1 幾何模型和網格

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

圖1 壩體幾何模型和有限差分網格

1.2 本構模型和材料參數

壩體堆石料采用摩爾-庫倫（Mohr-Coulomb）彈塑性本構模型，當剪應力水平在彈性范圍內時，通過粘滯阻尼達到能量耗散；當剪應力達到屈服強度時，通過塑性流動實現能量耗散[2]。該模型能模擬土石壩地震時的塑性剪切變形，并且由于采用了彈塑性本構模型，可以直接計算壩體的永久變形。壩體堆石料動力計算的初始剪切模量隨平均有效主應力的變化關系為

式中，C為動剪切模量系數；n為動剪切模量指數；兩者值由試驗確定。σm為平均有效主應力；Pa為大氣壓。

在FLAC3D中，根據式（1）通過FISH語言編程實現初始剪切模量隨平均有效主應力的非線性變化。

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

1.3 地震波濾波和基線校正

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

表1 堆石料的計算參數

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

圖3 輸入的地震波加速度時程曲線

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

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

1.4 動力邊界條件設置

圖4 加速度時程的功率譜

圖5 修正后的速度時程曲線

圖6 未修正和修正后的位移時程曲線

在動力計算中，模型邊界條件的選取是一項重要內容，因為邊界上波的反射會對動力計算結果產生影響。FLAC3D提供了靜態（粘滯）邊界和自由場邊界兩種邊界條件。自由場邊界條件要求把模型設置得足夠遠，模型的范圍取得足夠大，才能使模型邊界上波的反射盡可能的小，并且自由場邊界要求模型4個側面必須垂直。對于本文單獨以壩體為研究對象的模型，只能采用靜態（粘滯）邊界條件。

靜態（粘滯）邊界是Lysmer和Kuhlemeyer提出來的[4]，通過在模型邊界的法向和切向分別設置自由的阻尼器來實現吸收反射波。若在靜態（粘滯）邊界上施加動荷載，則只能施加時程應力。時程速度轉換為時程應力的公式如下

式中，σn、σs分別為施加在靜態（粘滯）邊界上的法向應力和切向應力；ρ為密度；Cp和Cs為通過介質的p波和s波的波速；υn和υs為邊界上的法向和切向速度分量。

1.5 阻尼的選取

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

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

圖7 速度時程曲線

圖8 速度時程功率譜

2 計算結果分析

壩頂中心點的豎向永久變形見圖9。從圖9可知，壩體豎向永久變形達到22.3 cm。

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

圖9 壩頂豎向永久位移時程曲線

圖10 加速度時程曲線

3 結語

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

［1］陳育民，徐鼎平.FLAC/FLAC3D基礎與工程實例［M］.北京：中國水利水電出版社，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.