孫偉建，侯興民，李遠東，鄭珊珊

(煙臺大學土木工程學院，山東煙臺264005)

基于虛單元法求解滲流自由面的曲線擬合法

孫偉建，侯興民，李遠東，鄭珊珊

(煙臺大學土木工程學院，山東煙臺264005)

滲流在壩坡穩定及壩坡變形中起關鍵作用，滲流自由面的確定是滲流計算的主要內容之一。目前，求解滲流自由面經常采用虛單元法，但是該法需不斷調整通過自由面的單元，迭代效率慢且易造成單元畸形。對有限元法求解的自由面節點水頭做曲線擬合，無需調整單元，即可得到高精度的滲流自由面曲線。應用所提方法求解工程實例，計算結果表明該方法精度高、穩定性好、計算效率高，計算結果與實際自由面形狀更吻合。

滲流自由面；固定網格法；有限元；虛單元法

0 引 言

從20世紀開始，滲流對工程的影響已受到廣泛重視。無壓自由面的滲流研究是水利工程、巖土工程等問題的一個重要研究課題。只有確定出自由面的位置，才能確定滲流區的范圍，才能正確的認識水在壩坡穩定及變形分析中的作用，以減輕水的滲流對結構的危害。由于自由面(二維滲流分析時為浸潤線，三維滲流分析時為浸潤面)和溢出面的位置屬于非線性邊界問題，如何高效簡潔的處理并確定出滲流自由面的位置，以及保證迭代過程收斂的穩定性是該問題的關鍵。

目前，自由面問題的求解方法主要是有限元法，有限元法分為移動網格法和固定網格法兩大類。移動網格法[1]是確定滲流自由面的傳統有限元法，該方法先根據經驗假定一個自由面，并將其作為可移動邊界，在迭代過程中不斷調整網格形狀，修改自由面的位置，直到自由面位置穩定為止。移動網格法雖迭代收斂性好，但存在缺陷，例如事先確定的自由面位置必須和實際自由面位置較為接近，每次需重新形成滲透矩陣，計算量大，同時，在不斷調整網格時易造成單元網格的畸形或重疊，影響計算精度。

為了克服變網格法的局限，學者們開始嘗試采用擴大分析系統的辦法，自Neumann于1973年提出用不變網格分析有自由面滲流的Galerkin法以來，先后出現了多種固定網格法，將滲流區域與非滲流區域結合起來形成總分析系統，并引入邊界條件作分析，以達到求解過程中無需變更單元網格的目的。方法大致分為兩類：一類是調整流量法，如王媛、潘樹來等[2-3]的改進初流量法等；另一類是調整單元滲透矩陣法，如Bathe的單元滲透矩陣調整法、朱軍的改進單元滲透矩陣調整法[4]、王賢能等的高斯點法[5]、陳昌祿的變單元滲透系數法[6]。此外，為了解決或部分解決單元、網格對數值分析的嚴重限制，應用無單元法來求解滲流自由面[7-11]。葛錦宏等將無單元法引入到有自由面的滲流計算中，提出了有自由面滲流分析的無單元法；沈振中在采用罰函數處理滲流邊界條件的同時，給出了罰因子的具體表達式，將無單元法用于求解堤壩滲流問題；M.Darbani提出了自由面淺水方程的無網格法；Yu-xin Jie提出了基于Voronoi形狀函數的自然單元法；Tang Jing采用罰函數無單元法求解復雜土石壩滲流場。以上工作均在自由面求解的精度以及效率方面取得較大的改進。另外，許多學者在此基礎上加以創新，采用不變網格基礎上的變網格法，以提高計算的精度和效率，如吳夢喜等提出的虛單元法[12]，黃蔚的丟單元法[13]，蘇枋的丟棄單元法在三維滲流分析中的應用[14]。

本文在虛單元法的基礎上，提出了一種滲流自由面的曲線擬合法，選取適應能力強的三角形單元進行有限元分析，無需重新調整跨自由面單元，即可得到精度很高的曲線自由面。

1 滲流的基本理論和邊界條件

根據達西定律和地下水運動的連續性條件，不考慮土體和水體的壓縮性，均質各向異性土體的穩定滲流滿足偏微分方程

(1)

式中，h是水頭函數；kx、kz分別為以坐標x、z為滲透主方向的滲透系數。

滲流基本微分方程的邊界條件有兩類。第一類邊界條件又稱水頭邊界條件，為邊界上給定位勢函數或水頭分布。在穩定滲流場中，淹沒水中的滲流邊界為等勢面或等水頭面，即常數，該常數通常為上、下游水位。設Γ1為具有給定水頭的邊界，則邊界條件可寫為

(2)

第二類邊界條件又稱流量邊界條件，為在邊界上給出位勢函數或水頭的法向導數，表示為

(3)

2 曲線擬合法原理

虛單元法以上一次有限元計算求得的節點勢為基礎，用插值法與逼近法求出自由面與單元邊線的交點，將自由面穿越的母單元作為變形單元，自由面將單元分成上下兩個區域，自由面以上的區域在下一次計算時不參與形成滲透矩陣，自由面以下為滲流計算區域，該域逐步逼近滲流實際區域，所得流場逐步逼近真實解。

(1)對土石壩進行全域規則劃分。首先對壩體進行豎向和橫向剖分，得到規則的四邊形單元，再對各四邊形單元進行同方向對角線連接，剖分為三角形單元，節點編號按照由左向右的順序，由上而下依次編號。

(2)釋放第一類邊界條件，得到所有節點的水頭值，即求解矩陣線性代數方程

(4)

(3)由自由面滿足的邊界條件h=z，即在跨自由面單元中其邊界條件h=z與h=N1h1+N2h2+N3h3聯立求得自由面的位置，其中Ni為形函數。

(4)根據最小二乘法，由全域劃分求解出的自由面的水頭值擬合一條高次遞減曲線，使曲線滿足目標函數

(5)

式中，ωi為自由面節點的權重；φ*(xi)為待求近似函數；m為過自由面節點數目；φ為函數類型。由于最小二乘法不能保證近似函數通過已知點，可通過選取適當的權函數將此誤差控制在一定的范圍。本文所用的方法為三次曲線擬合，自由面節點的權重取1，所擬合的高次曲線即為精確的自由面。

進行曲線擬合時，縱向水頭差值不宜過大，以確保擬合曲線為順次遞減曲線，要求在劃分網格時，單元適當的劃分多一些，即可確保水頭求解的精確性，又可確保擬合的曲線更加符合自由面形狀。為確保擬合曲線的準確性，本文采用基于全域總勢能確定溢出點，此方法可以較為準確地確定溢出點。

表1 解析解與數值解自由面位置數據對比 m

圖1 均質矩形壩網格剖分示意

圖2 自由面線形示意

3 算例

3.1 有解析解矩形均質土壩

以10 m×10 m均質土壩作分析，土壩上下游水位分別為10 m和2 m，底部為不透水邊界，下游滲出面邊界取為h=z。該問題自由面的解析表達式為y2=100-8x。建立3個單元尺寸不同的計算模型(見圖1)用于有限元計算：①模型1。網格劃分為200個直角邊均為1m的三節點直角三角形平面單元；②模型2。網格劃分為50個直角邊均為2m的三節點直角三角形平面單元；③模型3。網格劃分為18個三結點直角三角形平面單元。利用自編的二維有限元程序計算3個模型，自由面線形分別見圖2，用本文方法求得的自由面位置與虛單元法精度到0.01m求得的自由面位置、解析解結果對比，結果見表1。

由表1中的數據可以看出，虛單元法(精度為0.01m)隨著網格的加密，計算精度也越來越高。在進行有限元計算時網格劃分的越密，計算所需的時間會越長。本文方法不需重新劃分網格，僅需一次全域求解，通過最小二乘法進行高次曲線擬合即可得到精確解。本文方法與解析解最大誤差分別為0.18、0.17、0.26m，精度比虛單元法高。由圖2可看出本文方法擬合的高次曲線與解析解曲線更加接近，圖形更加吻合。

3.2 有解析解的均質梯形壩

梯形壩滲流在工程中比較常見，如用于擋水的土石壩和基坑開挖的圍堰堰體等?，F以不透水地基上均質土石壩為例，壩高17m、長97m，上、下游水位分別為15、4m。本算例中網格劃分為240和714個三角形單元，大網格劃分單元及求解曲線分別見圖3、圖4。利用本文編制的有限元程序求解壩體自由面的位置與虛單元法求解的自由面位置進行比較見表2。

圖3 均質梯形壩大網格剖分

圖4 均質梯形壩大網格剖分下三種自由面的比較

由表2及圖4可以看出，本文方法在求解自由面時與虛單元法求解自由面位置相差不大，虛單元法為精確到0.01m時的自由面位置。虛單元法在求解自由面是需不斷重新調整跨自由面單元，容易造成網格畸形，嚴重影響計算精度，且需不斷調整自由面節點位置，計算效率較低。本文方法僅需劃分一次網格，不需對跨自由面單元進行處理，對全域求解后，用最小二乘法對跨自由面單元中高程等于水頭的點進行高次曲線擬合，形成一條近似真實自由面曲線。由此算例可看出，本文方法在求解自由面時無需重新劃分網格，求解效率高，穩定性好，精度高。

4 結 論

滲流在壩坡穩定及壩坡變形中起關鍵性作用，滲流自由面的確定是滲流計算的關鍵內容之一。本文提出了一種曲線擬合法，并編制有限元程序，應用于計算模型及工程實例中得到如下結論：

表2 兩種數值方法計算不同網格劃分自由面位置數據對比 m

(1)本文在全域求解的基礎上，應用最小二乘法理論對全域求解后的自由面節點進行逐次遞減的高次曲線擬合，無需多次迭代，簡化了求解過程。

(2)本文方法與虛單元法相比，無需重新劃分網格，只需進行一次網格劃分即可求解出精確的自由面曲線，避免了全域求解后重新調整跨浸潤面單元造成的網格畸形。

(3)通過實際工程算例表明，即使采用大網格劃分，本文方法也可較高精度地求解出自由面曲線。

[1]ZIENKIEWIEZ O C， TAYLOR R L. The finite element method[M]. New York： Mcgraw-Hill， 1991．

[2]王媛. 求解有自由面滲流問題的初流量法的改進[J]. 水利學報， 1998(3)： 68-73．

[3]潘樹來， 王全鳳， 俞縉. 利用初流量法分析有自由面滲流問題之改進[J]. 巖土工程學報， 2012， 34(2)： 202-209．

[4]朱軍， 劉光廷. 改進的單元滲透矩陣調整法求解無壓滲流場[J]. 水利學報， 2001(8)： 49-52．

[5]王賢能， 黃潤秋. 有自由面滲流分析的高斯點法[J]. 水文地質工程地質， 1997(6)： 1-4．

[6]陳昌祿， 潘文彥， 王曉章. 求解滲流自由面的變單元法[J]. 巖土工程技術， 2005， 19(4)： 166-169．

[7]李廣信， 葛錦宏， 介玉新. 有自由面滲流的無單元法[J]. 清華大學學報： 自然科學版， 2002， 42(11)： 1552-1555．

[8]沈振中， 陳小虎， 吳越建. 求解堤壩滲流場的罰函數無單元法[J].河海大學學報：自然科學版，2008，26(1)：44-48．

[9]DARBANI M， OUAHSINE A， VILLON P， et al. Meshless method for shallow water equations with free surface flow[J]. Applied Mathematics and Computation， 2011(217)： 5113-5124．

[10]JIE Yuxin， LIU Lizhen， XU Wenjie， et al. Application of NEM in seepage analysis with a free surface[J]. Mathematics and Computers in Simulation， 2013(89)： 23-37．

[11]TANG Jing， LIU Yongbiao. Penalty Function Element Free Method to Solve Complex Seepage Field of Earth Fill Dam[C]. IERI Procedia， 2012(1)： 117-123．

[12]吳夢喜， 張雪勤. 有自由面滲流分析的虛單元法[J]. 水利學報， 1994(8)： 64-71．

[13]黃蔚， 劉迎曦， 周承芳. 三維無壓滲流場的有限元算法研究[J]. 水利學報， 2001(6)： 33-36．

[14]蘇枋， 王建祥， 熱依汗. 丟棄單元法在三維有限元滲流分析中的應用[J]. 水力發電， 2009， 35(3)： 26-28.

(責任編輯 王 琪)

A Curve Fitting Method Based on Virtual Element Method for Solving Seepage Free Surface

SUN Weijian, HOU Xingmin, LI Yuandong, ZHENG Shanshan

(School of Civil Engineering, Yantai University, Yantai 264005, Shandong, China)

The seepage plays a key role in the stability and deformation of dam slope, so the free surface determination is one of main contents in seepage calculation. At present, the virtual element method is often used to solve this problem, but the method needs to adjust the element passing through the free surface continuously, which makes the iteration efficiency slow and causes the element deformity. A curve fitting method is proposed herein to obtain a high-accurate free surface without re-adjusting the elements. Practical engineering examples show that the method is efficient, stable and accurate．

seepage free surface; fixed grid method; FEA; virtual element method

2016-01-27

山東省自然科學基金資助項目(ZR2012EEL31)

孫偉建(1991—)，男，山東濰坊人，碩士研究生，主要從事巖土工程研究．

TV139.14

A

0559-9342(2016)11-0050-04