黃超 歐陽君（湖南省水利水電勘測設(shè)計研究總院 長沙市 410007）

基于迭代原理求解穩(wěn)定滲流自由面位置的方法探討

黃超 歐陽君
（湖南省水利水電勘測設(shè)計研究總院 長沙市 410007）

為解決水工建筑物穩(wěn)定滲流分析中的自由面求解這一難題，文章提出了基于迭代原理求解穩(wěn)定滲流自由面位置的方法。該方法與傳統(tǒng)方法不同，在迭代求解過程中僅需迭代求解下游逸出點，而不需首先計算自由面的位置，也不用改變網(wǎng)格，即可得到自由面分布。結(jié)果表明，該方法可以便捷、客觀地獲得自由面分布，大大簡化分析工作。

穩(wěn)定滲流 自由面 迭代法

有自由面的滲流是水工建筑物設(shè)計中經(jīng)常遇到的問題，如壩體滲流、邊坡滲流、各種閘壩的繞壩滲流等。由于自由面是待定的、且需同時滿足水頭函數(shù)及壓力條件，它的互補和非線性性質(zhì)給問題的求解帶了極大的困難。在目前一般采用數(shù)值方法對這類具有自由面的滲流問題進(jìn)行求解，有限元法因能處理復(fù)雜邊界、三維滲流和非線性等問題，故在各類數(shù)值方法中脫穎而出[1～3]。目前求解此類問題的有限元解法可以分為移動網(wǎng)格法和固定網(wǎng)格法兩種[1，4～7]。移動網(wǎng)格法是首先假定一個浸潤線，然后對浸潤線以下區(qū)域劃分網(wǎng)格，給定邊界條件計算此假定滲流域，接著根據(jù)計算結(jié)果調(diào)整浸潤線，重新劃分網(wǎng)格，逐步迭代求解，直到自由面位置穩(wěn)定為止。移動網(wǎng)格法在一定意義上可同時滿足自由面上的勢函數(shù)和壓力條件，但不能用同一網(wǎng)格計算滲流場和應(yīng)力場，而且無法計算自由面以上區(qū)域的其它物理量，因此不能進(jìn)行耦合分析。固定網(wǎng)格法的核心就是在計算中保持網(wǎng)格不變，采用了擴大的滲流區(qū)域和固定邊界（通常是全部區(qū)域和邊界）來求解各種各樣的滲流問題，一般可分為剩余流量法、單元傳導(dǎo)矩陣調(diào)整法、子單元法和初流量法等。固定網(wǎng)格法不能同時滿足水頭函數(shù)和壓力條件，在物理概念上有欠妥之處，而且計算較為復(fù)雜，依賴性強。迄今，不管是移動網(wǎng)格法還是固定網(wǎng)格法多采用迭代法來確定自由面，本文提出通過迭代定出逸出點高程來獲得問題的解，為滲流自由面的求解提供了新的研究思路。

1 計算方法

取x，z軸與滲流主方向一致，由水流的連續(xù)條件和Darcy定律可以得到二維滲流控制方程[1]：

式中 kx、kz——分別為水平向與豎向的滲透系數(shù)（cm/s）；

h——總水頭（m）；

SS——單位儲蓄水量（m3/m）。

要確定限定空間內(nèi)發(fā)生的流場分布僅靠滲流基本微分方程還是不夠的，邊值條件的取值也很關(guān)鍵。滲流問題的定解條件中最重要的就是邊界條件，一般存在三類邊界條件[3]：

第一類邊界條件為邊界上給定位勢函數(shù)或水頭分布，或稱水頭邊界條件。即：

第二類邊界條件為邊界上給出位勢函數(shù)或水頭的法向?qū)?shù)，或稱流量邊界條件，在這類邊界條件上流量是已知的。即：

第三類邊界條件為混合邊界條件，是指含水層邊界的內(nèi)外水頭差和交換的流量之間保持一定的線性關(guān)系，即：

式中，α、β為正常數(shù)，都是此類邊界各點的已知數(shù)，無量綱。

對于無壓滲流，內(nèi)部無體積源時，整個滲流區(qū)域內(nèi)水頭函數(shù)為：

式中 φ——飽和流水頭函數(shù)；

P——水壓力；

γ——水容重；

Z——自某基準(zhǔn)算起的高度。

對于有壓滲流，由于不存在自由面，邊界位置是確定的，求解方程（1）即可得到問題的解。對于有自由面的滲流問題，因為滲流邊界位置不能預(yù)先確定，因此滲流邊界是解題的關(guān)鍵。一般在土石壩平面穩(wěn)定滲流問題求解時，上游水位以下的上游壩坡、下游水位以下的下游壩坡及下游逸出點到下游水位之間的下游壩坡為給定水頭邊界，其余壩體輪廓取為不透水邊界，同時在滲流自由面上和下游逸出段各點的孔隙水壓都為0。即在自由面邊界上應(yīng)該同時滿足式（1）和式（2），且在式（2）中取H=z。但是滲流自由面的位置事先并不知道，需要在求解過程中通過不斷迭代來確定。目前采用迭代法進(jìn)行求解時，不管是固定網(wǎng)格法還是移動網(wǎng)格法都存在著各種各樣的問題[1]，在求解自由面的時候事先假定滲流自由面初始位置與實際位置相差過大，并且迭代計算比較復(fù)雜、易導(dǎo)致嚴(yán)重的誤差。

為了彌補上述兩種方法中求解滲流自由面中的不足，可以利用另一部分完全已知的邊界條件來求解。對土石壩平面穩(wěn)定滲流問題來說，整個滲流區(qū)域中下游壩坡存在著幾種邊值條件：壩頂?shù)较掠我莩鳇c的壩坡為不透水邊界，下游逸出點到下游水位之間的壩坡為零孔壓邊界，下游水位以下的壩坡為給定水頭邊界，這幾種邊值條件無法完全確定是因為下游逸出點高程不確定。該方法僅對下游逸出點高程進(jìn)行迭代求解，不需要事先假定自由面位置和求解整個自由面上的一系列節(jié)點的位置，自由面位置通過節(jié)點孔壓內(nèi)插求得；在迭代過程中不需要對滲流區(qū)域內(nèi)部網(wǎng)格進(jìn)行調(diào)整，僅要求在劃分初始網(wǎng)格時對下游可能逸出點附近網(wǎng)格局部加密，并且網(wǎng)格劃分也不需要考慮滲流自由面的可能位置，可以和一般應(yīng)力應(yīng)變分析及邊坡穩(wěn)定性分析共用同一套網(wǎng)格進(jìn)行耦合分析。

在用該方法求解滲流問題時，首先需假定下游逸出點的初始高程，再求解方程（1～5），得到滲流區(qū)域的滲流場，然后驗算零孔壓等值線與下游壩坡的交點是否與假定逸出點高程是否在允許誤差范圍內(nèi)，若在，計算結(jié)束；否則，修正逸出點的位置，并重復(fù)上述迭代過程，直到滿足精度要求為止。

綜上所述，該方法是在滲流邊界條件上開展，出于解決問題的方便，把自由面上下的滲流采用同一個模型、同一套網(wǎng)格進(jìn)行描述，并把整個壩體作為滲流區(qū)域進(jìn)行分析；核心是把滲流自由面看作一條幾乎不透水邊界而不是嚴(yán)格的完全不透水邊界，即在自由面以上，滲透流量幾乎為0，而不完全為0。

2 計算實例及驗證

為驗證本文方法的可靠性，下面就上述方法做出一個簡單算例。假設(shè)某均質(zhì)壩壩高50 m，壩頂寬10 m，上游水深40 m，下游水深5 m，上下游坡率為2.0和1.5，未設(shè)任何排水設(shè)施。

圖1 算例的有限元網(wǎng)格圖

該算例的有限元網(wǎng)格如圖1所示，它完全不受自由面位置的限制。同時，在假定逸出點的初始位置時，也不必?fù)?dān)心它與實際位置的差別太大而引起計算失真，一般都可假設(shè)下游水位為逸出點高程。用該方法對算例從下游水位開始迭代求解，采用二分法步步逼近，經(jīng)過4次迭代即可得到問題的解。算例逸出點計算的迭代過程如表1所示。從下游水位開始，總能使迭代出來的零孔壓點高程大于假定逸出點高程，迭代結(jié)果慢慢逼近實際逸出點。在第四次迭代時，逸出點的假定高程和計算高程一致，故可以認(rèn)為18.76 m為真實逸出點高程，計算結(jié)束。壩體孔隙水壓力的等值線分布如圖2所示，根據(jù)其物理意義，孔隙壓力為0的等值線即為浸潤線。這些結(jié)果與文獻(xiàn)[4]的結(jié)果是相吻合的；從計算結(jié)果中還直接地反映出，自由面以上的土體內(nèi)具有負(fù)孔壓，且自由面以上的滲流很小，這些都符合傳統(tǒng)滲流理論的定性解釋。

表1 迭代計算結(jié)果

圖2 壩體孔隙水壓力分布圖

算例計算還表明：

（1）下游可能逸出點附近的網(wǎng)格加密對結(jié)果的精度有一定影響。在迭代過程中，需要不斷修正逸出點的位置，因該處的計算結(jié)果是迭代收斂的控制因素，如果網(wǎng)格劃分使單元變得過于奇異，會導(dǎo)致計算無法收斂，對計算結(jié)果有一定的影響。

（2）迭代的精度可適當(dāng)降低。假定的逸出點高程與其實際值相差不太大時，孔壓分布和浸潤線位置相差都不大，只是在其下游段有一定的差別。

為了驗證本文方法的實用性，下面以某一大型心墻堆石壩為例進(jìn)行算例分析。大壩典型斷面、材料分區(qū)見圖3、單元剖分見圖4、材料特性見表2。計算工況為穩(wěn)定滲流期，上游水位為313 m，下游水位為221.57 m。

在計算中，將整個壩體、壩基和排水棱體進(jìn)行網(wǎng)格劃分，但對排水棱體的透水性賦予一個較大的值。按照本文方法，首先假定逸出點高程為下游水位與下游壩坡的交點，而計算得到的零孔壓等值線與排水棱體下游壩面的交點仍在此位置，因此可以認(rèn)為第一次試算的結(jié)果即為滿意的計算結(jié)果。圖5給出了計算的自由面位置和總水頭分布圖，圖中零孔壓線即為浸潤線，零孔壓線與下游壩坡的交點即為逸出點。計算結(jié)果表明該土石壩在心墻處零孔壓下出現(xiàn)了明顯的驟降，計算結(jié)果符合一般計算規(guī)律，計算合理，同時表明該堆石心墻壩的心墻防滲效果很明顯，有效地降低了下游壩坡的孔隙水壓力；從計算過程中還直接反映出這樣一個事實，當(dāng)下游水位與下游壩坡交點在排水棱體上時，在排水棱體的上游坡面不存在逸出段。

圖3 某心墻堆石壩典型斷面示意圖（單位：m）

圖4 某心墻堆石壩典型斷面有限元模型及網(wǎng)格剖分圖（單位：m）

表2 壩體材料滲透系數(shù)指標(biāo)

圖5 某心墻堆石壩滲流期壩體孔隙水壓力分布圖（單位：m）

3 結(jié) 論

（1）在迭代求解滲流自由面時，無需改動網(wǎng)格，而是用求解下游逸出點的方法來求解自由面。原理可靠，工作量小，實現(xiàn)方便，可以推廣到復(fù)雜土工邊值問題的滲流分析。

（2）經(jīng)過對實際工程的算例分析論證，計算結(jié)果符合計算規(guī)律，滿足工程滲流計算精度要求。

（3）該方法最大限度地利用了通用軟件使用方便、功能齊全、運算速度快和計算精度高等優(yōu)點，為實際工程問題的分析提供了新的研究思路。

[1]李廣信.高等土力學(xué)[M].北京:清華大學(xué)出版社，2002.

[2]錢家歡，殷宗澤.土工原理與計算[M].北京:水利電力出版社，1994.

[3]杜延齡，許國安.滲流分析的有限元法和電網(wǎng)絡(luò)法[M].北京:水利電力出版社，1992.

[4]張有天，陳平，王鐳.有自由面滲流分析的初參數(shù)法[J].水利學(xué)報，1988，（8）:18-26.

[5]柴軍瑞，件彥卿.巖體裂隙網(wǎng)絡(luò)滲流自由面的確定方法[J].工程勘察，2000，（01）:20-21.

[6]徐千軍，張建紅.確定穩(wěn)定滲流自由面位置的一種簡便方法[J].水動力學(xué)研究與進(jìn)展，1999，（12）:418-423.

[7]陳洪凱，唐紅梅.滲流自由面求解的基本方法、修正及應(yīng)用[J].重慶交通學(xué)院學(xué)報，1997，16（3）:5-10.

[8]劉杰.土的滲透穩(wěn)定及滲流控制[M].北京：水利電力出版社，1992.

2016-05-03）

黃超（1981-），湖南寧鄉(xiāng)人，大學(xué)本科，工程師，從事水利水電工程建筑設(shè)計工作，手機：13975807897。