滲流自由面數(shù)值模擬方法比較

蔣勝銀，李連俠，廖華勝，楊 華，鄒 俊

（1.四川大學(xué)水力學(xué)及山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗室，成都 610065；2.珠江水利科學(xué)研究院，廣州 510611）

滲流自由面數(shù)值模擬方法比較

蔣勝銀1，李連俠1，廖華勝1，楊 華1，鄒 俊2

（1.四川大學(xué)水力學(xué)及山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗室，成都 610065；2.珠江水利科學(xué)研究院，廣州 510611）

結(jié)合最新提出的自由面適應(yīng)網(wǎng)格法，與現(xiàn)行的2種固定網(wǎng)格法即變單元滲透系數(shù)法和初流量法進(jìn)行數(shù)值模擬比較，以探討各種方法在自由面滲流計算中的優(yōu)缺點(diǎn)及其應(yīng)用情況。試探性地選取3種固定網(wǎng)格方法進(jìn)行系統(tǒng)的比較，即在精度、迭代步數(shù)、網(wǎng)格敏感度、初設(shè)自由面敏感性上進(jìn)行多方位比較，分析各方法的靈活性、魯棒性、數(shù)值穩(wěn)定性、計算效率和可靠性。結(jié)果表明，在同一網(wǎng)格系統(tǒng)下，自由面適應(yīng)網(wǎng)格法較其它2種方法具有計算速度快、精度較高等優(yōu)點(diǎn)，其迭代步數(shù)和結(jié)果精度對初設(shè)自由面和網(wǎng)格密度敏感度均較其它2種方法低。

滲流自由面；數(shù)值模擬；固定網(wǎng)格方法；自由面適應(yīng)網(wǎng)格方法

1 概 述

由于與水相聯(lián)系，很多水工建筑物都有滲流問題。滲流會影響建筑物的受力情況以及結(jié)構(gòu)物的材料性質(zhì)等，所以，滲流關(guān)乎水工建筑物的安全性和可靠性。尤其在各類擋水建筑物中，例如：攔河壩、深基坑圍堰、閘壩等。如果不能有效應(yīng)對滲流，則可能發(fā)生建筑物滑動、土體滲透破壞、滲漏水量過大等問題。

國內(nèi)外以往滲流研究的內(nèi)容一般是過流體中浸潤線的幾何位置，過流體的水頭場、水力坡降場和流速分布，滲流量等。以這些信息為基礎(chǔ)，進(jìn)行滲流的整體或局部滲流穩(wěn)定性分析、質(zhì)量的評估等。

滲流計算的數(shù)值模擬方法按網(wǎng)格性質(zhì)的不同還分為變網(wǎng)格法、固定網(wǎng)格法和無網(wǎng)格法。本文探討的主要是固定網(wǎng)格法。固定網(wǎng)格法即在計算中保持網(wǎng)格不變，采用擴(kuò)大的滲流區(qū)域和固定邊界（通常是全部區(qū)域和邊界）來求解各種各樣的滲流問題。自Neuman于1973年提出用不變網(wǎng)格分析帶有自由面滲流問題的Galerkin方法［1］以來，出現(xiàn)了多種固定網(wǎng)格法，比較有影響的有Desai（1976）［2］穩(wěn)定滲流計算的剩余流量法，并于1983年推廣到非穩(wěn)定滲流計算；Bathe（1979）［3］的單元滲透矩陣調(diào)整法；張有天［4］（1988）的初流量法。另外還有變分不等式法［5］、結(jié)點(diǎn)虛流量法［6］、虛單元法［7］、穩(wěn)定滲流計算的截止負(fù)壓法［8］、高斯點(diǎn)法［9］、變單元滲透系數(shù)法［10］等。

本文結(jié)合最新提出的自由面適應(yīng)網(wǎng)格法［11］，與現(xiàn)行的幾種固定網(wǎng)格法如變單元滲透系數(shù)法和初流量法進(jìn)行數(shù)值模擬比較，以探討各種方法在自由面滲流計算中的優(yōu)缺點(diǎn)及其應(yīng)用情況。

2 常見滲流自由面模擬方法

在數(shù)學(xué)范疇，自由面滲流和無自由面滲流的數(shù)學(xué)模型基本是一樣的，但是無自由面滲流沒有自由面邊界和逸出邊界，在數(shù)值模擬時，是一個線性問題，求解方便。反之，自由面滲流則因為自由面位置也為求解變量之一，對其求解則成為非線性問題，使得滲流計算變得復(fù)雜和困難。通常用迭代逼近的方法來求其近似解。具有自由面三維滲流問題的數(shù)學(xué)模型可用下式表示：

式中，K為滲透系數(shù)張量；H＝H（x，y，z，t）為水頭；Ω為滲流計算區(qū)域；ΩH為源（匯）項；SS為儲水率；t0為計算初始時刻，H0（x，y，z）為t0時刻Ω域的水頭分布；Γ1為第1類（給定水頭）邊界，H1（x，y，z，t）為Γ1上的水頭分布；Γ2為第2類（給定流量）邊界；n2為Γ2的法向方向；q（x，y，z，t）為Γ2上的流量分布；Γ3為第3類邊界，即滲流自由面邊界；ε為Γ3上的補(bǔ)給強(qiáng)度；μ為給水度。

針對這一問題，為求解浸潤面提出了諸多方法。本文主要針對滲流數(shù)值模擬方法中的固定網(wǎng)格法進(jìn)行討論。對于固定網(wǎng)格法而言，一般分為3類：第1類是弱化滲流自由面以上區(qū)域的影響，比如變單元滲透系數(shù)法、丟單元法；第2類是補(bǔ)償流量的方法，比如剩余流量法、初流量法；第3類是自由面適應(yīng)網(wǎng)格方法。為對比這3類方法的優(yōu)劣，本文選取了3種方法：變單元滲透系數(shù)法、初流量法、自由面適應(yīng)網(wǎng)格法進(jìn)行比較。

2.1 變單元滲透系數(shù)法

為求解滲流自由面邊界的變單元滲透系數(shù)法，具體做法如下：

（1）在上游水位與下游水位之間先估計出初始自由面的區(qū)域。

（2）對全區(qū)域進(jìn)行網(wǎng)格劃分，在預(yù)估的初始自由面域內(nèi)將網(wǎng)格劃分得很密。

（3）自由面上結(jié)點(diǎn)的水頭等于其位置勢，將自由面之上的所有結(jié)點(diǎn)進(jìn)行標(biāo)識并判斷處于自由面上的所有單元。

（4）把位于自由面上單元的滲透系數(shù)改為一個很小的數(shù)。由于自由面之上的單元在第2次以后的迭代計算過程中，不需要參與計算，為了讓自由面之上的單元對以后的計算不產(chǎn)生影響，方便編程序，所以把自由面以上的單元在迭代過程中都乘以一個很小的數(shù)，比如為10×10－3，足以在計算過程中把自由面以上的單元滲流結(jié)果忽略，不再考慮自由面以上單元對滲流的影響。

（5）進(jìn)行下一次求解，將本次求出的結(jié)點(diǎn)勢與上一次求出的結(jié)點(diǎn)勢比較，當(dāng)滲流域內(nèi)結(jié)點(diǎn)均滿足收斂條件時，滲流場中所有滿足｜H－Z｜≤ε2的結(jié)點(diǎn)的連線即為自由面。

2.2 初流量法

初流量法是在達(dá)西定律中增加一初流量項q0i，通過對初流量值的調(diào)整，將一系列非線性分析化為線性分析。若以整個區(qū)域為考察對象（包括非飽和區(qū)），可將達(dá)西定律改寫為

該式采用了張量記法，kijhj為介質(zhì)滲透張量；q0i表示初流量值。

在計算中，通過第一步迭代的結(jié)果找出自由面以上單元和包含自由面的單元，通過單元高斯點(diǎn)的水頭與高程的關(guān)系，判斷是否計算初流量。形成了新的Q項（等效結(jié)點(diǎn)流量列陣）。解矩陣方程，通過這次的結(jié)果計算高斯點(diǎn)的流量，又形成新的Q項，迭代計算，直到滿足收斂條件。

初流量法有3種收斂標(biāo)準(zhǔn)：①結(jié)點(diǎn)初流量的絕對值不超過某一允許量；②計算出的自由面穩(wěn)定，若干次迭代后h＜z的高斯點(diǎn)數(shù)不變；③飽和區(qū)結(jié)點(diǎn)水頭增量不超過某一允許量。本文選用的是后2個收斂標(biāo)準(zhǔn)。

2.3 自由面適應(yīng)網(wǎng)格方法［11］

滲流計算中浸潤線計算一般迭代過程可描述如下：

（2）求解滲流場，求出滲流自由面上的水頭分布；

（4）重復(fù)第（2）和第（3）步，直到第n＋1和n次迭代所求出的水頭Hn＋1和Hn的相對誤差值在精度要求范圍內(nèi)；

與其他方法不同，自由面適應(yīng)網(wǎng)格法在浸潤線（自由面）的調(diào)整過程中，不需要將網(wǎng)格變形，而是將浸潤線按其所在位置，近似映射到最近的在既定網(wǎng)格系統(tǒng)中求解，其中的映射既可以是簡單的零階映射，也可是高階映射，它最大限度地利用了浸潤線上H＝z的邊界條件特性。其具體方法如下（見圖1）：圖中a，b，c，d，e，f為其與縱向網(wǎng)格線的交點(diǎn)，abcdef為計算過程中可能的浸潤線中間結(jié)果。對于浸潤線上每一分段，按照距離節(jié)點(diǎn)位置最近原則把該段映射到計算節(jié)點(diǎn)上，如ab段就被映射到節(jié)點(diǎn)A1，A2，B上，依此類推，計算浸潤線abcdef最后就按該原則被映射在固定網(wǎng)格系統(tǒng)中的A1A2BC1C2DE1E2F，該折線被稱為映射浸潤線，實(shí)現(xiàn)了自由面適應(yīng)網(wǎng)格的過程。設(shè)｜A1a｜表示節(jié)點(diǎn)A1到點(diǎn)a的距離，｜Bb｜表示節(jié)點(diǎn)B到點(diǎn)b的距離，依此類推，則如果｜A1a｜，｜A2a｜，｜Bb｜，｜Cc｜，｜Dd｜，｜E1e｜，｜E2e｜，｜Ff｜的最小值滿足精度要求，即計算浸潤線和映射浸潤線足夠接近，就認(rèn)為迭代收斂，否則進(jìn)行下一點(diǎn)的迭代，直至收斂為止。

圖1 自由面適應(yīng)網(wǎng)格方法實(shí)現(xiàn)原理示意圖Fig.1 Principle of adaptive grid for free surface

3 對比算例及比較參數(shù)擬定

3.1 對比算例及計算工況

梯形壩滲流是工程中常見的滲流現(xiàn)象，比如用于擋水的土石壩和基坑開挖的圍堰堰體等建筑物其體型均屬梯形壩。均勻介質(zhì)的梯形壩的自由面滲流是有理論解的，通過這個理論解與3種模擬方法結(jié)果進(jìn)行比較，可以進(jìn)一步說明問題和發(fā)現(xiàn)問題。

選取的梯形壩體型如圖2，上游坡面邊坡系數(shù)m1＝3，下游坡面邊坡系數(shù)為m2＝2。上游水位H1＝15 m，下游水位H2＝4 m。該梯形壩模型屬于水平不透水層上均質(zhì)壩。上游液體將通過邊界AB滲入壩體，在壩內(nèi)形成自由表面（浸潤面）AC，C點(diǎn)稱為逸出點(diǎn)，ABDC區(qū)域為滲流區(qū)。

圖2 梯形壩模型及自由面理論計算示意圖Fig.2 Trapezoidal dam model and theoretical solution of free surface

上述問題浸潤曲線的解析式［12］為

其中逸出點(diǎn)水深hk由（4）、（5）兩式聯(lián)立求解得出。（3）至（6）式中符號參見圖2。

上游滲流段A′B′GC所通過的單寬滲流量：

通過下游出滲段CGD的單寬滲流量為

等效的矩形體的寬度ΔL由下式確定：

假定一系列的x值，由（3）式可得到相應(yīng)的y值，從而描繪出梯形壩滲流浸潤線如圖3。表1給出該自由面解析解幾個典型位置的解。

圖3 梯形壩理論自由面示意圖Fig.3 Theoreticalwater free surface of trapezoidal dam

表1 梯形壩自由面理論解Table 1 Theoretical solution of water free surface of trapezoidal dam

求得該梯形壩的理論解后，用3種數(shù)值模擬方法進(jìn)行相應(yīng)計算，并與上述理論解進(jìn)行計算精度、初設(shè)自由面敏感度、網(wǎng)格敏感度和CPU時間的比較。計算工況分為2種網(wǎng)格和3種初設(shè)自由面，其計算工況見表2。

表2 計算工況Table 2 Com putational cases m

圖4 初設(shè)自由面示意Fig.4 Initial free surface

初設(shè)自由面均為直線，其分布見圖4。取值如下：初設(shè)自由面1，上游15m，下游15m；初設(shè)自由面2，上游15m，下游4m；初設(shè)自由面3，上游15m，下游8 m。

3.2 比較參數(shù)擬定

為分析相應(yīng)方法的靈活性、魯棒性、數(shù)值穩(wěn)定性、計算效率和可靠性，針對以上選定方法和算例，擬定以下幾個參數(shù)進(jìn)行量化比較：①計算結(jié)果精確度；②初設(shè)自由面敏感度；③網(wǎng)格變化敏感度；④計算的CPU時間。

圖5 不同工況下3種方法的計算結(jié)果精度對比Fig.5 Comparison of calculation accuracy of threemethods in different cases

4 模擬結(jié)果分析及比較

4.1 結(jié)果精度

3種模擬方法的計算結(jié)果精度對比見圖5。一共6個不同的工況。3種方法結(jié)果與理論解的偏差程度進(jìn)行對比，可看出它們求解準(zhǔn)確性的相對差異。

可見，在同一工況下自由面適應(yīng)網(wǎng)格法的結(jié)果總體較為穩(wěn)定；初流量法的求解結(jié)果偏大，變滲透系數(shù)法結(jié)果一般。

4.2 初設(shè)自由面敏感度

對初設(shè)自由面敏感度的對比見圖6。每種方法用3種初設(shè)自由面下的結(jié)果，計算各點(diǎn)數(shù)值的標(biāo)準(zhǔn)差。每種方法得到3種初設(shè)自由面下的3組標(biāo)準(zhǔn)差后，與其他的方法進(jìn)行比較。

由圖可見，這個算例的初設(shè)自由面敏感度三者大致相當(dāng)。

4.3 網(wǎng)格敏感度

對網(wǎng)格敏感度的對比見圖7，采用計算結(jié)果的標(biāo)準(zhǔn)差進(jìn)行量化比較。

由圖可見，自由面適應(yīng)網(wǎng)格法的標(biāo)準(zhǔn)差最?。▊€別點(diǎn)較大），變滲透系數(shù)法次之，初流量法最大。這說明，自由面適應(yīng)網(wǎng)格法和變滲透系數(shù)法的網(wǎng)格敏感度較小，具有較好的適應(yīng)性和魯棒性，并且自由面適應(yīng)網(wǎng)格法對于尺寸較粗的網(wǎng)格也能得到較精確的結(jié)果，初流量法表現(xiàn)一般。

4.4 CPU時間

同種工況下CPU時間對比見表3。同種初設(shè)自由面下，當(dāng)網(wǎng)格數(shù)目較少時，自由面適應(yīng)網(wǎng)格法和變滲透系數(shù)法用時明顯少于初流量法；且隨著網(wǎng)格數(shù)目的增加，自由面適應(yīng)網(wǎng)格法的用時增加速率遠(yuǎn)小于其它2種方法，CPU時間遠(yuǎn)小于其它2種方法，且網(wǎng)格數(shù)越多這種優(yōu)勢越明顯?？梢娮杂擅孢m應(yīng)網(wǎng)格法計算效率更高，更適用于大區(qū)域和復(fù)雜滲流問題的計算。

圖6 不同初設(shè)自由面下3種方法計算結(jié)果的標(biāo)準(zhǔn)差Fig.6 Standard deviation of the results of three methods w ith different initial free surfaces

圖7 不同網(wǎng)格下3種方法計算結(jié)果的標(biāo)準(zhǔn)差Fig.7 Standard deviation of the results of threemethods in different grids

表3 網(wǎng)格1和網(wǎng)格2情況下3種初設(shè)自由面計算時間Table 3 The calculation time of three differentmethods w ith three initial free surfaces in grid 1 and grid 2 s

綜上所述，3種方法的綜合對比見表4，可看出自由面適應(yīng)網(wǎng)格法綜合能力優(yōu)于其他2種方法：在大型計算時，變滲透系數(shù)的時間花費(fèi)太大，小型計算時其效果較好；初流量法的計算指標(biāo)居中。

表4 3種方法綜合對比Table 4 Com prehensive com parison of threemethods

5 結(jié) 論

本文選取了3種計算滲流浸潤線固定網(wǎng)格方法，并擬定了方法間用于比較的參數(shù)：靈活性、魯棒性、數(shù)值穩(wěn)定性、計算效率及可靠性等。比較結(jié)果表明：

（1）自由面適應(yīng)網(wǎng)格方法和變滲透系數(shù)法一樣，具有良好的可靠性和計算精度。

（2）初設(shè)自由面對自由面適應(yīng)網(wǎng)格方法影響較小，該方法具有較好的適應(yīng)性和魯棒性；初流量法受初設(shè)自由面的影響相對較大，變滲透系數(shù)法居于二者之間。

（3）初流量法和變滲透系數(shù)法迭代步數(shù)受網(wǎng)格疏密影響較大，隨著網(wǎng)格數(shù)的增加而遞增；而自由面適應(yīng)網(wǎng)格法的計算結(jié)果受網(wǎng)格數(shù)目的影響很小，隨著網(wǎng)格數(shù)的遞增，迭代步數(shù)沒有發(fā)生劇增，說明自由面適應(yīng)網(wǎng)格方法具有良好的數(shù)值穩(wěn)定性，在大區(qū)域和復(fù)雜滲流計算中具有良好的應(yīng)用前景。

（4）自由面適應(yīng)網(wǎng)格法進(jìn)行自由面調(diào)整時，方法簡單，除了映射計算，沒有涉及過多的額外計算，而且該方法收斂較快，其CPU時間沒有因為網(wǎng)格數(shù)目的增大而有數(shù)量級的增大，與其他方法相比優(yōu)勢明顯，說明自由面適應(yīng)網(wǎng)格方法具有良好的計算效率，更適用于大區(qū)域和復(fù)雜滲流計算中。

［1］ NEUMAN S P.Saturated-Unsaturated Seepage by Finite Elements［J］.Journal of Hydraulic Division，ASCE，1973，99（12）：2233－2250.

［2］ DESAIC S.Finite Element Residual Shemes for Unconfined Flow［J］.International Journal for Numerical Methods in Engineering，1976，10（6）：1415－1418.

［3］ BATHE K J.Finite Element Free Surface Seepage Analy-sisWithout Mesh Iteration［J］.International Journal for Numerical and Analytical Methods in Geomechanics，1979，3（1）：13－22.

［4］ 張有天，陳 平，王 鐳.有自由面滲流分析的初流量法［J］.水利學(xué)報，1988，（8）：18.（ZHANG You-tian，CHEN Ping，WANG Lei.Initial Flow Method for Seepage Analysiswith Free Surface［J］.Journal of Hydraulic Engineering，1988，（8）：18.（in Chinese））

［5］ BAIOCCHI C，COMINCIOLI V，MAGENES E，et al.Free Boundary Problems in the Theory of Fluid Flow Through Porous Media：Existence and Uniqueness Theoremes［J］.Annali Di Matematica Pura ed Applicata，1973，97（4）：1－82.

［6］ 速寶玉，朱岳明.不變網(wǎng)格確定滲流自由面的節(jié)點(diǎn)虛流量法［J］.河海大學(xué)學(xué)報，1991，19（5）：113－117.（SU Bao-yu，ZHU Yue-ming.Nodal Virtual Discharge Method for Seepage Problems with Free Surface on Fixed Grids［J］.Journal of Hohai University，1991，19（5）：113－117.（in Chinese））

［7］ 吳夢喜，張學(xué)勤.有自由面滲流分析的虛單元法［J］.水利學(xué)報，1994，（8）：64－71.（WU Meng-xi，ZHANG Xue-qin.Imaginary Element Method for Numerical Analysis of Seepagewith Free Surface［J］.Journal of Hydraulic Engineering，1994，（8）：64－67.（in Chinese））

［8］ 速寶玉，沈振中，趙 堅.用變分不等式理論求解滲流問題的截止負(fù)壓法［J］.水利學(xué)報，1996，（3）：22－29，35.（SU Bao-yu，SHEN Zhen-zhong，ZHAO Jian.The Cut-off Negative Pressure Method for Soling Filtration Problems Based on the Theory of Variational Inequalities［J］.Journal of Hydraulic Engineering，1996，（3）：22－29，35.（in Chinese））

［9］ 王賢能，黃潤秋.有自由面滲流分析的高斯點(diǎn)法［J］.水文地質(zhì)工程地質(zhì)，1997，（6）：1－4.（WANG Xianneng，HUANG Run-qiu.Gavss PointMethod for Seepage Analysis with Free Surface［J］.Hydrogeology and Engineering Geology，1997，（6）：1－4.（in Chinese））

［10］黨發(fā)寧，王曉章，鄭忠安，等.有自由面滲流分析的變單元滲透系數(shù)法［J］.西北水力發(fā)電，2004，（1）：1－3.（DANG Fa-ning，WANG Xiao-zhang，ZHENG Zhongan，et al.Variable Element Seepage Coefficient Method for Seepage Numerical Analysis of Seepagewith Free Surface［J］.Journal of Northwest Hydroelectric Power，2004，（1）：1－3.（in Chinese））

［11］李連俠，廖華勝，劉 達(dá)，等.滲流計算中求解浸潤線的自由面適應(yīng)網(wǎng)格方法［J］.四川大學(xué)學(xué)報（工程科學(xué)版），2006，（5）：76－81.（LI Lian-xia，LIAO Hua-sheng LIU Da，etal.A Method ofWater Free Surface Adapting Grid System to Simulate Phreatic Surface in Seepage Flow［J］.Journal of Sichuan University（Engineering Science Edition），2006，（5）：76－81.（in Chinese））

［12］吳持恭.水力學(xué)（下冊）［M］.北京：高等教育出版社，2008：234－238.（WU Chi-gong.Hydraulics（Volume Two）［M］.Beijing：Higher Education Press，2008：234－238.（in Chinese） ）

（編輯：王 慰）

Comparison of Numerical Approaches of Simulating Seepage Flow w ith Free Surface

JIANG Sheng-yin1，LILian-xia1，LIAO Hua-sheng1，YANG Hua1，ZOU Jun2
（1.State Key Laboratory of Hydraulics and Mountain River Engineering，Sichuan University，Sichuan 610065，China；2.Pearl River Hydraulics Research Institute，Guangzhou 510611，China）

The fixed grid methods can be classified into initial flow method，variable seepage coefficientmethod，variational inequalitymethod，node imaginary flow method，imaginary elementsmethod and the water free surface adaptive grid method which was put forward recently.Each method has its advantages and disadvantages.To ultimately guide us to choose propermethods when processing the practical seepage problems in large scale，it is of great value to evaluate their applicability including flexibility，robustness，numerical stability，computational efficiency and reliability by comparing the accuracy，iterative steps，grid sensitivity，and initial free surface sensitivity of three different fixed grid methods.The results show that the water free surface adaptive grid method is of higher computing efficiency and accuracy than the other twomethods in the same grid system.The initiative steps are less，and the accuracy of computing results is less sensitive to the initial free surface and grid density than the other two methods.

free surface of seepage flow；numerical simulation；fixed grid methods；method of water free surface adaptive grid

TV139.1

A

1001－5485（2011）07－0037－06

2010-08-30

教育部博士點(diǎn)新教師基金（20090181120013），四川大學(xué)青年基金（0030614132002）

蔣勝銀（1986-），男，四川成都人，碩士研究生，主要從事水工水力學(xué)方面的研究，（電話）15528330512（電子信箱）scujsy＠163.com。

李連俠（1978-），男，河南南陽人，副教授，主要從事水工水力學(xué)方面的研究，（電話）13308021514（電子信箱）lianxiali＠yahoo.com.cn。