滲流自由面數值模擬方法比較

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

（1.四川大學水力學及山區河流開發保護國家重點實驗室，成都 610065；2.珠江水利科學研究院，廣州 510611）

滲流自由面數值模擬方法比較

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

（1.四川大學水力學及山區河流開發保護國家重點實驗室，成都 610065；2.珠江水利科學研究院，廣州 510611）

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

滲流自由面；數值模擬；固定網格方法；自由面適應網格方法

1 概 述

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

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

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

本文結合最新提出的自由面適應網格法［11］，與現行的幾種固定網格法如變單元滲透系數法和初流量法進行數值模擬比較，以探討各種方法在自由面滲流計算中的優缺點及其應用情況。

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

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

式中，K為滲透系數張量；H＝H（x，y，z，t）為水頭；Ω為滲流計算區域；Ω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上的補給強度；μ為給水度。

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

2.1 變單元滲透系數法

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

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

（2）對全區域進行網格劃分，在預估的初始自由面域內將網格劃分得很密。

（3）自由面上結點的水頭等于其位置勢，將自由面之上的所有結點進行標識并判斷處于自由面上的所有單元。

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

（5）進行下一次求解，將本次求出的結點勢與上一次求出的結點勢比較，當滲流域內結點均滿足收斂條件時，滲流場中所有滿足｜H－Z｜≤ε2的結點的連線即為自由面。

2.2 初流量法

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

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

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

初流量法有3種收斂標準：①結點初流量的絕對值不超過某一允許量；②計算出的自由面穩定，若干次迭代后h＜z的高斯點數不變；③飽和區結點水頭增量不超過某一允許量。本文選用的是后2個收斂標準。

2.3 自由面適應網格方法［11］

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

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

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

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

圖1 自由面適應網格方法實現原理示意圖Fig.1 Principle of adaptive grid for free surface

3 對比算例及比較參數擬定

3.1 對比算例及計算工況

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

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

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

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

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

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

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

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

假定一系列的x值，由（3）式可得到相應的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種數值模擬方法進行相應計算，并與上述理論解進行計算精度、初設自由面敏感度、網格敏感度和CPU時間的比較。計算工況分為2種網格和3種初設自由面，其計算工況見表2。

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

圖4 初設自由面示意Fig.4 Initial free surface

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

3.2 比較參數擬定

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

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

4 模擬結果分析及比較

4.1 結果精度

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

可見，在同一工況下自由面適應網格法的結果總體較為穩定；初流量法的求解結果偏大，變滲透系數法結果一般。

4.2 初設自由面敏感度

對初設自由面敏感度的對比見圖6。每種方法用3種初設自由面下的結果，計算各點數值的標準差。每種方法得到3種初設自由面下的3組標準差后，與其他的方法進行比較。

由圖可見，這個算例的初設自由面敏感度三者大致相當。

4.3 網格敏感度

對網格敏感度的對比見圖7，采用計算結果的標準差進行量化比較。

由圖可見，自由面適應網格法的標準差最小（個別點較大），變滲透系數法次之，初流量法最大。這說明，自由面適應網格法和變滲透系數法的網格敏感度較小，具有較好的適應性和魯棒性，并且自由面適應網格法對于尺寸較粗的網格也能得到較精確的結果，初流量法表現一般。

4.4 CPU時間

同種工況下CPU時間對比見表3。同種初設自由面下，當網格數目較少時，自由面適應網格法和變滲透系數法用時明顯少于初流量法；且隨著網格數目的增加，自由面適應網格法的用時增加速率遠小于其它2種方法，CPU時間遠小于其它2種方法，且網格數越多這種優勢越明顯。可見自由面適應網格法計算效率更高，更適用于大區域和復雜滲流問題的計算。

圖6 不同初設自由面下3種方法計算結果的標準差Fig.6 Standard deviation of the results of three methods w ith different initial free surfaces

圖7 不同網格下3種方法計算結果的標準差Fig.7 Standard deviation of the results of threemethods in different grids

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

綜上所述，3種方法的綜合對比見表4，可看出自由面適應網格法綜合能力優于其他2種方法：在大型計算時，變滲透系數的時間花費太大，小型計算時其效果較好；初流量法的計算指標居中。

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

5 結 論

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

（1）自由面適應網格方法和變滲透系數法一樣，具有良好的可靠性和計算精度。

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

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

（4）自由面適應網格法進行自由面調整時，方法簡單，除了映射計算，沒有涉及過多的額外計算，而且該方法收斂較快，其CPU時間沒有因為網格數目的增大而有數量級的增大，與其他方法相比優勢明顯，說明自由面適應網格方法具有良好的計算效率，更適用于大區域和復雜滲流計算中。

［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］.水利學報，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］ 速寶玉，朱岳明.不變網格確定滲流自由面的節點虛流量法［J］.河海大學學報，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］ 吳夢喜，張學勤.有自由面滲流分析的虛單元法［J］.水利學報，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］ 速寶玉，沈振中，趙 堅.用變分不等式理論求解滲流問題的截止負壓法［J］.水利學報，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］ 王賢能，黃潤秋.有自由面滲流分析的高斯點法［J］.水文地質工程地質，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］黨發寧，王曉章，鄭忠安，等.有自由面滲流分析的變單元滲透系數法［J］.西北水力發電，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］李連俠，廖華勝，劉 達，等.滲流計算中求解浸潤線的自由面適應網格方法［J］.四川大學學報（工程科學版），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］吳持恭.水力學（下冊）［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

教育部博士點新教師基金（20090181120013），四川大學青年基金（0030614132002）

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

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