金海 黃明鎮(zhèn) 高宇航 曾勇文

摘要：巖土在動(dòng)荷載作用下的應(yīng)力-應(yīng)變關(guān)系非常復(fù)雜，在進(jìn)行土石壩有限元?jiǎng)恿Ψ治鰰r(shí)，需要選擇合適的巖土本構(gòu)關(guān)系模型?；陴椥岳碚摰牡刃Ь€性本構(gòu)關(guān)系模型因其計(jì)算效率高、計(jì)算結(jié)果合理而得到廣泛使用。使用ANSYS平臺(tái)的User Programmable Features（UPFs）開(kāi)發(fā)了基于黏彈性理論的等效線性本構(gòu)關(guān)系模型。使用 SHAKE91說(shuō)明文件中的經(jīng)典案例對(duì)比分析和土石壩的地震反應(yīng)分析進(jìn)行有效性驗(yàn)證，所得結(jié)果符合一般規(guī)律，驗(yàn)證了該模型的可靠性。結(jié)果表明：以ANSYS平臺(tái)為基礎(chǔ)開(kāi)發(fā)的等效線性模型可用于復(fù)雜的三維巖土結(jié)構(gòu)的動(dòng)力分析。研究成果為大型巖土結(jié)構(gòu)工程的動(dòng)力分析問(wèn)題提供了一種新的選擇。

關(guān)鍵詞：巖土等效線性模型;二次開(kāi)發(fā);ANSYS;巖土動(dòng)力分析

中圖法分類號(hào)：TU435文獻(xiàn)標(biāo)志碼：ADOI：10.15974/j.cnki.slsdkb.2021.08.014

文章編號(hào)：1006 - 0081（2021）08 - 0073 - 04

土石壩等巖土結(jié)構(gòu)動(dòng)力分析研究的主要難點(diǎn)在于構(gòu)造有效的土體應(yīng)力-應(yīng)變模型來(lái)模擬巖土結(jié)構(gòu)在動(dòng)力作用下的響應(yīng)。1966年Clough[1]將有限單元法引入巖土力學(xué)領(lǐng)域后，Dibaj[2]、Idriss[3]等相繼提出了巖土結(jié)構(gòu)的非線性分析方法，直到1972年，Hardin和Drnevich[4]把土體作為黏彈性材料，提出了可以反映土體材料非線性和黏滯性的Hardin-Drnevich等效線性模型。1996年，沈珠江[5]結(jié)合等效線性黏彈性模型以及殘余應(yīng)變經(jīng)驗(yàn)公式，提出了沈珠江等效黏彈塑性模型，可直接計(jì)算土石壩在地震作用下的殘余變形，又根據(jù)大型振動(dòng)三軸試驗(yàn)結(jié)果，對(duì)堆石料動(dòng)力本構(gòu)模型作了系統(tǒng)研究，改進(jìn)了堆石料的等效黏彈性模型，在后續(xù)土石壩的動(dòng)力研究中被廣泛應(yīng)用。

本文基于ANSYS-UPFs，結(jié)合經(jīng)沈珠江改進(jìn)的等效黏彈塑性模型進(jìn)行二次開(kāi)發(fā)，通過(guò)SHAKE91說(shuō)明文件中的經(jīng)典案例對(duì)比分析和土石壩的地震反應(yīng)分析進(jìn)行有效性驗(yàn)證。

1 基于黏彈性理論的等效線性本構(gòu)模型

在循環(huán)荷載作用下，土體會(huì)表現(xiàn)出非線性、滯后性和變形累積等應(yīng)力-應(yīng)變關(guān)系特點(diǎn)。對(duì)于某些簡(jiǎn)單問(wèn)題，可將三者分別加以考慮，在特定范圍內(nèi)能得到較好的結(jié)果。對(duì)于復(fù)雜問(wèn)題，則需要將3種特性結(jié)合考慮。

1.1 等效線性模型的表達(dá)方式

在巖土動(dòng)力分析中需要考慮動(dòng)剪切模量和阻尼比這兩個(gè)參數(shù)，在循環(huán)荷載作用下，動(dòng)剪切模量和阻尼比會(huì)隨著土體剪應(yīng)變的變化而改變。因此，在等效線性模型中引入等效剪切模量[G]和等效阻尼比[D]的概念，將滯回特性用等效阻尼比[D]與剪應(yīng)變幅值[γa]的關(guān)系反映;而土體的骨架曲線則用等效剪切模量[G]與剪應(yīng)變幅值[γa]的關(guān)系反映。

本文使用沈珠江[5]提出的如下形式確定等效剪切模量[G]和等效阻尼比[D]：

1.2 子程序編制思路

由式（2）可知，壩體的初始剪切模量[Gmax]與震前的平均有效應(yīng)力有關(guān)，因此需在動(dòng)力分析前進(jìn)行靜力分析，得到各單元的初始應(yīng)力狀態(tài)并計(jì)算平均有效應(yīng)力，可通過(guò)APDL命令“INISTATE，WRITE，1，，，，，S”將應(yīng)力狀態(tài)輸出為IST文件。在第一步計(jì)算時(shí)，通過(guò)INISTATE，READ‘保存的文件名，‘IST命令讀入單元的初始應(yīng)力，計(jì)算出平均有效應(yīng)力，進(jìn)而計(jì)算出初始剪切模量[Gmax]。由于在后續(xù)計(jì)算中需要用到初始剪切模量，且該變量在地震過(guò)程中保持不變，因此使用子程序中的狀態(tài)變量statev（1）保存，在第一步計(jì)算的起始步更新 ，后續(xù)過(guò)程無(wú)需修改。

在第一次迭代前假定初始動(dòng)剪應(yīng)變?yōu)?，并根據(jù)式（1）～（4） 確定初始剪切模量比[G1Gmax]和阻尼比[D1]。然后，通過(guò)線性分析方法求解得到各單元的剪應(yīng)變并更新。在一次迭代結(jié)束后，通過(guò)更新的應(yīng)變水平按照式（5）～（7）確定各單元的剪切模量和阻尼比，并進(jìn)入下一次迭代。由于剪切模量比和阻尼比在每次計(jì)算中保持不變，在子程序中用statev（2）和statev（3）保存，并在每次計(jì)算起始步更新。最大剪應(yīng)變[γdmax]用statev（4）保存，并在每個(gè)增量步進(jìn)行更新。等效線性模型結(jié)構(gòu)框架見(jiàn)圖1。

1.3 巖土材料液化判別

由于地震激勵(lì)的周期和幅值是不固定的，結(jié)構(gòu)模型頻率的取值和邊界約束的施加會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生較大影響。在實(shí)際計(jì)算中，可能會(huì)使局部單元的剪應(yīng)變過(guò)大。實(shí)際上，當(dāng)土體等效動(dòng)剪應(yīng)力大于抗液化動(dòng)強(qiáng)度時(shí)，土體發(fā)生液化。在一定振次下達(dá)到某一應(yīng)變標(biāo)準(zhǔn)時(shí)的動(dòng)應(yīng)力幅值稱為土的動(dòng)強(qiáng)度。

要對(duì)土體材料進(jìn)行液化的判別，需要獲取各單元的等效動(dòng)剪應(yīng)力和材料動(dòng)強(qiáng)度。而動(dòng)強(qiáng)度需要根據(jù)材料的動(dòng)強(qiáng)度與破壞周次關(guān)系曲線和地震等效振動(dòng)次數(shù)確定。一般情況下，在后處理中進(jìn)行液化的判別。本文考慮土體液化對(duì)計(jì)算結(jié)果的影響，在程序中事先設(shè)定一個(gè)液化應(yīng)變的判別標(biāo)準(zhǔn)（如取[8%]）、土體液化后剪切模量比（[1%～5%]）和阻尼比（[Dmax～2Dmax]），當(dāng)剪應(yīng)變大于設(shè)定的液化應(yīng)變時(shí)使用液化后的等效參數(shù)。

2 算例驗(yàn)證

2.1 算例1

本例使用費(fèi)康[6]、謝倫武[7]等的SHAKE91說(shuō)明文件[8]中的經(jīng)典算例來(lái)驗(yàn)證子程序的可行性和可靠性。算例為一坐落于基巖上的水平地基，由砂土和黏土組成。土體材料參數(shù)及相應(yīng)的等效線性參數(shù)見(jiàn)表1。砂土和黏土的剪切模量比、阻尼比與剪應(yīng)變關(guān)系見(jiàn)圖2。輸入的地震加速度時(shí)程曲線見(jiàn)圖3。

如圖4（a）所示，在ANSYS前處理中建立高45 m（y方向）的土柱進(jìn)行分析。土柱使用三維八結(jié)點(diǎn)SOLID185單元，共剖分為5 632個(gè)單元。土體材料使用自定義的等效線性材料參數(shù)。最大剪切模量按表1給定值輸入。土柱地面施加全約束，頂面無(wú)約束，側(cè)面只允許發(fā)生x向變形，并沿x向輸入地震波，輸入的地震波加速度峰值設(shè)置為0。圖4（b）為最大地震反應(yīng)加速度沿地基深度分布。本文計(jì)算方法與SHAKE91計(jì)算方法相比較，加速度沿深度方向的變化符合變化趨勢(shì)，土層底部加速度的變化很小，在土層中上部加速度變化明顯，且上部的加速度放大較為明顯。本文方法和SHAKE91求得的土層頂面絕對(duì)加速度值分別為0.315 g和0.29 g，誤差在容許范圍內(nèi)，約為8.6%。

將輸入的地震波加速度峰值調(diào)整為0.000 1 m/s2，重復(fù)上述計(jì)算，得到加速度放大系數(shù)與輸入的峰值加速度關(guān)系曲線，如圖5所示。

2.2 算例2

對(duì)于一個(gè)高100 m、壩頂寬10 m、心墻頂高6 m心墻土石壩，上下游對(duì)稱，坡比為1∶2，心墻坡比為1∶0.2，如圖6所示。土石壩使用三維八結(jié)點(diǎn)SOLID185單元，共剖分為14 000個(gè)單元，等效線性材料參數(shù)如表2所示。輸入的河道向加速度按如圖7所示加速度時(shí)程曲線取最大加速度為0，豎向加速度取河道向加速度的2/3。

觀察圖8（a）心墻中心線上最大加速度分布，在70 m以下的最大加速度變化不明顯，而在70 m以上部分最大加速度隨著壩高急速上升，并在壩頂處達(dá)到最大，最大加速度分布存在著鞭梢效應(yīng)。觀察圖8（b）心墻中心線上最大剪應(yīng)變分布，在60 m以下變化較小，在60～80 m急速增大，在80 m以上又迅速減小。最大剪應(yīng)變出現(xiàn)在壩體2/3處，與土石壩受地震作用后，常出現(xiàn)裂縫或產(chǎn)生較大變形的位置相符。

3 結(jié) 語(yǔ)

基于ANSYS-UPFs編譯了用于巖土動(dòng)力分析的等效線性本構(gòu)模型用戶子程序，并通過(guò)SHAKE91說(shuō)明文件中的案例建立了三維土層進(jìn)行對(duì)比分析，計(jì)算結(jié)果符合實(shí)際情況，證明了本文基于ANSYS-UPFs開(kāi)發(fā)的等效線性模型的可靠性和適用性。通過(guò)土石壩算例的計(jì)算以及對(duì)計(jì)算結(jié)果的分析證明了本文方法可利用ANSYS的前后處理優(yōu)勢(shì)對(duì)大型巖土結(jié)構(gòu)進(jìn)行應(yīng)力-應(yīng)變動(dòng)力分析。

參考文獻(xiàn)：

[1] CLOUGH R W，CHOPRA A K. Earthquake Stress Analysis in Earth Dams[C]// Proc.ASCE，No.EM2，1996.

[2] DIBAJ? M， PENZIEN J. Nonlinear Seismic Response of Earth Structures[R]. Report EERC-69-2，1969.

[3] IDRISS I M，SEED H B，SERFF N. Seismic Response by vari-ble Damping Elements[C]// Proc.ASCE，No.GT1，1974.

[4] HARDIN B O，DRNEVICH V P. Shear modulus and damping in soils design equations and curves[J]. Journal of Soils Mechanics and Foundation Division，1972，98（7）：667-692.

[5] 沈珠江，徐剛. 堆石料的動(dòng)力變形特性[J]. 水利水運(yùn)科學(xué)研究，1996（2）：143-150.

[6] 費(fèi)康，劉漢龍.ABAQUS的二次開(kāi)發(fā)及在土石壩靜、動(dòng)力分析中的應(yīng)用[J].巖土力學(xué)，2010，31（3）：881-890.

[7] 謝倫武，熊峰，姚梓渝，等.基于MATLAB和ABAQUS的土體等效線性化方法二次開(kāi)發(fā)[J].地震工程與工程振動(dòng)，2015，35（1）：135-142.

[8] IDRISSIM，SUN J I. Users manual for SHAKEE91[M].Davis： University of California，1992.

（編輯：李 慧）