基于ANSYS的某土石壩滲流數值模擬研究

虎 珀，漆文邦，譚明卓，俞增鑫

（四川大學水利水電學院,四川 成都 610065）

1 引言

在土石壩與水庫失事事故的統計中，約有1/4是由于滲流問題引起的，這表明深入研究滲流問題和設計有效的控制滲流措施是十分重要的。對于土石壩滲流分析問題，有水力學法、流網法和有限元法等方法[1]。有限元方法可以更好地適應復雜的邊界條件和壩體、非均質壩基、各向異性等不同的情況，所以在工程設計中逐漸得到廣泛應用。ANSYS軟件是美國ANSYS公司研制的大型通用有限元分析軟件，其中的熱模塊分析與滲流場分析在基本理論、微分方程、初始條件和邊界條件上基本一致，因此可利用ANSYS軟件對滲流場進行數值模擬，APDL是ANSYS Parmetric Design Language的縮寫，是一種類似FORTRAN的解釋性語言。利用這種語言可以實現一般程序語言的功能，如參數、宏、循環等。本文利用參數化設計語言APDL來實現滲流計算中的循環步驟，大大簡化計算者的操作步驟，同時得到較為準確的結果。

2 基于ANSYS滲流數值模擬的理論基礎

由于本文采用ANSYS軟件中熱模塊進行滲流場的數值模擬，下面需要論證溫度場與滲流場在基礎理論等方面的相似性。

2.1 基礎理論的相似性

對滲流場而言，根據達西定律：

式中：Q為流量；A為斷面面積；H為測壓管水頭；k為滲透系數。

對溫度場而言，根據熱傳導定律：

式中：Q為熱量；A為斷面面積；T為溫度；k為熱傳導傳熱系數。

2.2 微分方程的相似性

滲流場計算的微分方程形式如下：

對不可壓縮各向異性非均質無源穩定滲流，微分方程為：

對可壓縮各向異性非均質非穩定瞬態滲流，微分方程為：

式中：ksx、ksy、ksz為x、y、z方向的滲透系數；h為測壓管水頭；Ss為單位貯藏量。

對無熱源的各向異性非均質穩定熱傳導，微分方程為：

對無熱源的各向異性非均質瞬態熱傳導，微分方程為：

式中：ktx、kty、ktz為x、y、z方向的熱傳導傳熱系數；T為溫度；C為比熱容。

2.3 初始條件及邊界條件的相似性

滲流場的初始條件為：

溫度場的初始條件為：

第一類邊界條件：

第二類邊界條件：

綜上，滲流場與溫度場的計算在基礎理論、微分方程、初始條件和邊界條件等方面具有基本一致的形式，因此以熱模塊來進行滲流場的計算模擬是可行的。

3 基于APDL語言的ANSYS計算土石壩滲流的方法

進行滲流分析，首先要確定壩體內的浸潤線，進而計算后續的滲流比降以及滲流量[2]。利用ANSYS進行滲流計算時，由于無法自動對土體的飽和或非飽和進行判斷，即無法自動賦予單元應有的屬性，因此在計算時需要進行試算，利用ANSYS自帶的生死單元技術來輔助計算。同時需要利用APDL語言來實現循環計算，簡化操作步驟。

3.1 生死單元技術

根據用戶的需要，可以在加載過程中加入或移除材料，使特定的某些單元“不存在”或“存在”。這種技術即為生死單元技術。當需要激活單元的“死”時，ANSYS并非將“殺死”的單元從模型中刪除，而是將其剛度矩陣乘以一個很小的因子，使得被殺死單元的載荷、質量、阻尼和應力等被設置為零。

3.2 APDL語言的運用

計算過程中需要用到的許多數據都存儲于ANSYS軟件的數據庫中，為方便使用可以利用*GET函數將其提取出來，進行定義與賦值[3]。當需要使用自己定義的一組數據時，可以通過*DIM函數來進行賦值。

對于需要進行的循環計算，可使用*IF命令設定執行所滿足的條件，滿足條件即執行不滿足則停止。一般通過比較兩個數的數值來確定當前所滿足的條件值，這就需要把判斷條件均轉化為數值對象。

3.3 浸潤線的計算步驟

1）首先定義上游面水頭值，定義下游面水頭值時建立函數令定義點水頭值等于Y坐標值。

部分APDL語言表示如下：

LSEL,S,,,P51X

NSLL,S,1

D,ALL,TEMP,,!定義上游面水頭

NSEL,R,LOC,Y,,

DO,I,1,NMAX

D,NXY(I),TEMP,NY(NXY(I))!定義下游面水頭等于Y坐標

ENDDO

2）以下游水頭值等于上游進行一次計算，計算之后開啟迭代計算模式，對每個單元的水頭進行計算，以單元水頭值小于Y坐標值為判斷基準殺死這部分單元。再令剩余的活單元的水頭值等于其位置水頭值，以剩余的活單元再次進行滲流計算。再根據結果重復判斷并執行殺死單元步驟，循環往復直至兩步計算的水頭結果之差小于最大允許誤差5×10-5時即可停止。

部分APDL語言表示如下：

DO,I,1,EMAX

DO,NN,1,4

GET,NP(KK),ELEM,I,NODE,NN

ENDDO

STY=(NST(NP(1))+NST(NP(2))+NST(NP(3))+NST(NP(4)))/4!得到單元平均水頭值

PT=STY-CENTRY(I)!判斷水頭值與Y值關系

IF,PT,GT,0,THEN

EALIVE,I

ELSEIF,PT,LE,0,THEN

EKILL,I

ENDIF

3）每步對單元死活進行定義后，都應該重新定義出口邊界條件：將計算出的浸潤線出口位置與前一次的出口位置進行比較，若計算出的出口位置低于前一次，則刪除原出口處邊界條件，用計算出的滲流出口位置作為下游水頭邊界。隨著計算循環，直至兩次出口邊界相等，則計算結束。利用這樣的方法可以準確的定義隨著循環進行動態變化的下游邊界條件。

部分APDL語言表示如下：

IF,EXITY2,NE,EXITY,THEN!判斷前后兩次計算的浸潤線出口位置是否相同

NSEL,R,LOC,Y,EXITY!選擇最高節點

IF,CKMAXY,GT,0,THEN

DDELE,ALL,TEMP!刪除出口最高節點邊界條件

ENDIF

計算結束后，得到總水頭、壓力水頭、流速等計算結果。將流速映射在出口邊界路徑上得到斷面的滲流量。

4 工程實例

以某?。?）型水庫土石壩的滲流計算分析為例，利用ANSYS、Seep/W與理正巖土分別進行計算。

4.1 工程概況

某小（2）型水庫是一座以灌溉供水為主，兼顧下游防洪要求的水庫。控制集水面積6.4 km3，總庫容12.1 萬m3。樞紐主要建筑物由大壩、溢洪道及輸水涵洞三部分組成。大壩為黏土心墻壩，壩頂高程1904.4 m，正常蓄水位1900.5 m；設計洪水位1902.8 m；校核洪水位1903.73 m。壩軸線長約104.2 m，壩頂平均寬3.2 m；上游迎水坡平均坡比1∶2.3，下游平均坡比1∶2.0。

4.2 計算條件

4.2.1 計算模型

計算模型選取該土石壩標準斷面，以上游坡腳為原點，x軸正方向為順河向指向下游，y軸正方向為鉛直向上方向。計算范圍選取水平方向距上游壩坡13 m至距下游壩坡19 m處。同時對模型進行一定簡化，得到計算模型見圖1。

圖1 滲流計算模型圖

4.2.2 網格劃分

計算模型采用四邊形單元PLANE55,進行網格劃分后，形成3579 個單元，3778 個節點，見圖2。

圖2 計算模型網格劃分圖

4.2.3 計算參數

計算工況為上游水位為校核洪水位1903.73 m，邊界條件為上游坡面采用水頭邊界，下游坡面水頭邊界通過迭代計算得出。

滲流計算所采用各材料的參數見表1。

表1 滲流計算各材料滲透系數表

4.3 計算結果

計算結果見圖3、圖4。

圖3 ANSYS計算總水頭成果云圖

圖4 ANSYS計算壓力水頭成果云圖

4.4 與其他軟件對比

以相同的模型與計算參數，分別在理正巖土、Seep/W中進行計算，計算出的總水頭和壓力水頭結果見圖5～ 圖8。

圖5 理正計算總水頭成果圖

圖6 理正計算壓力水頭成果圖

圖7 Seep/W計算總水頭成果云圖

圖8 Seep/W計算壓力水頭成果云圖

經計算，ANSYS計算得到的大壩單寬滲流量為1.718×10-5m3/s，理正巖土計算得到的大壩單寬滲流量為1.851×10-5m3/s，Seep/W計算得到的大壩單寬滲流量為1.810×10-5m3/s。根據以上結果分析可知，大壩滲流量較小，粘土心墻防滲效果較好。

5 結論

1）ANSYS、理正巖土、Seep/W三種軟件在計算滲流的結果基本一致，理正巖土計算的浸潤線與其他二者結果稍有不同，可能是因為其需要設定下游水位而缺少了對滲流出逸點及浸潤線的試算而造成的。

2）三種軟件計算出的滲流量結果基本一致，理正巖土計算的最大滲透坡降出現在心墻與覆蓋層連接處，而ANSYS與Seep/W計算的最大滲透坡降出現在心墻內部。

綜上，基于ANSYS熱模塊并利用APDL語言進行土石壩滲流數值模擬是可行的，其計算結果與專業巖土計算軟件Seep/W計算結果很接近。同時ANSYS相對于其他計算滲流的有限元軟件有著可開發性強的優點，深入研究后可以進一步提高滲流計算的準確性及效率。