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

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

1 引言

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

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

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

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

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

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

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

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

2.2 微分方程的相似性

滲流場(chǎng)計(jì)算的微分方程形式如下：

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

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

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

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

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

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

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

滲流場(chǎng)的初始條件為：

溫度場(chǎng)的初始條件為：

第一類邊界條件：

第二類邊界條件：

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

3 基于APDL語(yǔ)言的ANSYS計(jì)算土石壩滲流的方法

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

3.1 生死單元技術(shù)

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

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

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

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

3.3 浸潤(rùn)線的計(jì)算步驟

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

部分APDL語(yǔ)言表示如下：

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

部分APDL語(yǔ)言表示如下：

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

部分APDL語(yǔ)言表示如下：

IF,EXITY2,NE,EXITY,THEN!判斷前后兩次計(jì)算的浸潤(rùn)線出口位置是否相同

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

IF,CKMAXY,GT,0,THEN

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

ENDIF

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

4 工程實(shí)例

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

4.1 工程概況

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

4.2 計(jì)算條件

4.2.1 計(jì)算模型

計(jì)算模型選取該土石壩標(biāo)準(zhǔn)斷面，以上游坡腳為原點(diǎn)，x軸正方向?yàn)轫樅酉蛑赶蛳掠危瑈軸正方向?yàn)殂U直向上方向。計(jì)算范圍選取水平方向距上游壩坡13 m至距下游壩坡19 m處。同時(shí)對(duì)模型進(jìn)行一定簡(jiǎn)化，得到計(jì)算模型見(jiàn)圖1。

圖1 滲流計(jì)算模型圖

4.2.2 網(wǎng)格劃分

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

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

4.2.3 計(jì)算參數(shù)

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

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

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

4.3 計(jì)算結(jié)果

計(jì)算結(jié)果見(jiàn)圖3、圖4。

圖3 ANSYS計(jì)算總水頭成果云圖

圖4 ANSYS計(jì)算壓力水頭成果云圖

4.4 與其他軟件對(duì)比

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

圖5 理正計(jì)算總水頭成果圖

圖6 理正計(jì)算壓力水頭成果圖

圖7 Seep/W計(jì)算總水頭成果云圖

圖8 Seep/W計(jì)算壓力水頭成果云圖

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

5 結(jié)論

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

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

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