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

（四川大學(xué)水利水電學(xué)院,四川 成都 610065）

1 引言

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

2 基于ANSYS滲流數(shù)值模擬的理論基礎(chǔ)

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

2.1 基礎(chǔ)理論的相似性

對滲流場而言，根據(jù)達(dá)西定律：

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

對溫度場而言，根據(jù)熱傳導(dǎo)定律：

式中：Q為熱量；A為斷面面積；T為溫度；k為熱傳導(dǎo)傳熱系數(shù)。

2.2 微分方程的相似性

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

對不可壓縮各向異性非均質(zhì)無源穩(wěn)定滲流，微分方程為：

對可壓縮各向異性非均質(zhì)非穩(wěn)定瞬態(tài)滲流，微分方程為：

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

對無熱源的各向異性非均質(zhì)穩(wěn)定熱傳導(dǎo)，微分方程為：

對無熱源的各向異性非均質(zhì)瞬態(tài)熱傳導(dǎo)，微分方程為：

式中：ktx、kty、ktz為x、y、z方向的熱傳導(dǎo)傳熱系數(shù)；T為溫度；C為比熱容。

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

滲流場的初始條件為：

溫度場的初始條件為：

第一類邊界條件：

第二類邊界條件：

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

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

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

3.1 生死單元技術(shù)

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

3.2 APDL語言的運(yùn)用

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

對于需要進(jìn)行的循環(huán)計算，可使用*IF命令設(shè)定執(zhí)行所滿足的條件，滿足條件即執(zhí)行不滿足則停止。一般通過比較兩個數(shù)的數(shù)值來確定當(dāng)前所滿足的條件值，這就需要把判斷條件均轉(zhuǎn)化為數(shù)值對象。

3.3 浸潤線的計算步驟

1）首先定義上游面水頭值，定義下游面水頭值時建立函數(shù)令定義點(diǎn)水頭值等于Y坐標(biāo)值。

部分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坐標(biāo)

ENDDO

2）以下游水頭值等于上游進(jìn)行一次計算，計算之后開啟迭代計算模式，對每個單元的水頭進(jìn)行計算，以單元水頭值小于Y坐標(biāo)值為判斷基準(zhǔn)殺死這部分單元。再令剩余的活單元的水頭值等于其位置水頭值，以剩余的活單元再次進(jìn)行滲流計算。再根據(jù)結(jié)果重復(fù)判斷并執(zhí)行殺死單元步驟，循環(huán)往復(fù)直至兩步計算的水頭結(jié)果之差小于最大允許誤差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值關(guān)系

IF,PT,GT,0,THEN

EALIVE,I

ELSEIF,PT,LE,0,THEN

EKILL,I

ENDIF

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

部分APDL語言表示如下：

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

NSEL,R,LOC,Y,EXITY!選擇最高節(jié)點(diǎn)

IF,CKMAXY,GT,0,THEN

DDELE,ALL,TEMP!刪除出口最高節(jié)點(diǎn)邊界條件

ENDIF

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

4 工程實例

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

4.1 工程概況

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

4.2 計算條件

4.2.1 計算模型

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

圖1 滲流計算模型圖

4.2.2 網(wǎng)格劃分

計算模型采用四邊形單元PLANE55,進(jìn)行網(wǎng)格劃分后，形成3579 個單元，3778 個節(jié)點(diǎn)，見圖2。

圖2 計算模型網(wǎng)格劃分圖

4.2.3 計算參數(shù)

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

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

表1 滲流計算各材料滲透系數(shù)表

4.3 計算結(jié)果

計算結(jié)果見圖3、圖4。

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

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

4.4 與其他軟件對比

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

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

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

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

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

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

5 結(jié)論

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

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

綜上，基于ANSYS熱模塊并利用APDL語言進(jìn)行土石壩滲流數(shù)值模擬是可行的，其計算結(jié)果與專業(yè)巖土計算軟件Seep/W計算結(jié)果很接近。同時ANSYS相對于其他計算滲流的有限元軟件有著可開發(fā)性強(qiáng)的優(yōu)點(diǎn)，深入研究后可以進(jìn)一步提高滲流計算的準(zhǔn)確性及效率。