基于ANSYS-UPFs的巖土等效線性模型二次開發

金海 黃明鎮 高宇航 曾勇文

摘要：巖土在動荷載作用下的應力-應變關系非常復雜，在進行土石壩有限元動力分析時，需要選擇合適的巖土本構關系模型。基于黏彈性理論的等效線性本構關系模型因其計算效率高、計算結果合理而得到廣泛使用。使用ANSYS平臺的User Programmable Features（UPFs）開發了基于黏彈性理論的等效線性本構關系模型。使用 SHAKE91說明文件中的經典案例對比分析和土石壩的地震反應分析進行有效性驗證，所得結果符合一般規律，驗證了該模型的可靠性。結果表明：以ANSYS平臺為基礎開發的等效線性模型可用于復雜的三維巖土結構的動力分析。研究成果為大型巖土結構工程的動力分析問題提供了一種新的選擇。

關鍵詞：巖土等效線性模型;二次開發;ANSYS;巖土動力分析

中圖法分類號：TU435文獻標志碼：ADOI：10.15974/j.cnki.slsdkb.2021.08.014

文章編號：1006 - 0081（2021）08 - 0073 - 04

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

本文基于ANSYS-UPFs，結合經沈珠江改進的等效黏彈塑性模型進行二次開發，通過SHAKE91說明文件中的經典案例對比分析和土石壩的地震反應分析進行有效性驗證。

1 基于黏彈性理論的等效線性本構模型

在循環荷載作用下，土體會表現出非線性、滯后性和變形累積等應力-應變關系特點。對于某些簡單問題，可將三者分別加以考慮，在特定范圍內能得到較好的結果。對于復雜問題，則需要將3種特性結合考慮。

1.1 等效線性模型的表達方式

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

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

1.2 子程序編制思路

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

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

1.3 巖土材料液化判別

由于地震激勵的周期和幅值是不固定的，結構模型頻率的取值和邊界約束的施加會對計算結果產生較大影響。在實際計算中，可能會使局部單元的剪應變過大。實際上，當土體等效動剪應力大于抗液化動強度時，土體發生液化。在一定振次下達到某一應變標準時的動應力幅值稱為土的動強度。

要對土體材料進行液化的判別，需要獲取各單元的等效動剪應力和材料動強度。而動強度需要根據材料的動強度與破壞周次關系曲線和地震等效振動次數確定。一般情況下，在后處理中進行液化的判別。本文考慮土體液化對計算結果的影響，在程序中事先設定一個液化應變的判別標準（如取[8%]）、土體液化后剪切模量比（[1%～5%]）和阻尼比（[Dmax～2Dmax]），當剪應變大于設定的液化應變時使用液化后的等效參數。

2 算例驗證

2.1 算例1

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

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

將輸入的地震波加速度峰值調整為0.000 1 m/s2，重復上述計算，得到加速度放大系數與輸入的峰值加速度關系曲線，如圖5所示。

2.2 算例2

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

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

3 結 語

基于ANSYS-UPFs編譯了用于巖土動力分析的等效線性本構模型用戶子程序，并通過SHAKE91說明文件中的案例建立了三維土層進行對比分析，計算結果符合實際情況，證明了本文基于ANSYS-UPFs開發的等效線性模型的可靠性和適用性。通過土石壩算例的計算以及對計算結果的分析證明了本文方法可利用ANSYS的前后處理優勢對大型巖土結構進行應力-應變動力分析。

參考文獻：

[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] 沈珠江，徐剛. 堆石料的動力變形特性[J]. 水利水運科學研究，1996（2）：143-150.

[6] 費康，劉漢龍.ABAQUS的二次開發及在土石壩靜、動力分析中的應用[J].巖土力學，2010，31（3）：881-890.

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

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

（編輯：李 慧）