孫偉建，侯興民，李遠(yuǎn)東，鄭珊珊

(煙臺(tái)大學(xué)土木工程學(xué)院，山東煙臺(tái)264005)

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

孫偉建，侯興民，李遠(yuǎn)東，鄭珊珊

(煙臺(tái)大學(xué)土木工程學(xué)院，山東煙臺(tái)264005)

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

滲流自由面；固定網(wǎng)格法；有限元；虛單元法

0 引 言

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

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

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

本文在虛單元法的基礎(chǔ)上，提出了一種滲流自由面的曲線擬合法，選取適應(yīng)能力強(qiáng)的三角形單元進(jìn)行有限元分析，無(wú)需重新調(diào)整跨自由面單元，即可得到精度很高的曲線自由面。

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

根據(jù)達(dá)西定律和地下水運(yùn)動(dòng)的連續(xù)性條件，不考慮土體和水體的壓縮性，均質(zhì)各向異性土體的穩(wěn)定滲流滿足偏微分方程

(1)

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

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

(2)

第二類邊界條件又稱流量邊界條件，為在邊界上給出位勢(shì)函數(shù)或水頭的法向?qū)?shù)，表示為

(3)

2 曲線擬合法原理

虛單元法以上一次有限元計(jì)算求得的節(jié)點(diǎn)勢(shì)為基礎(chǔ)，用插值法與逼近法求出自由面與單元邊線的交點(diǎn)，將自由面穿越的母單元作為變形單元，自由面將單元分成上下兩個(gè)區(qū)域，自由面以上的區(qū)域在下一次計(jì)算時(shí)不參與形成滲透矩陣，自由面以下為滲流計(jì)算區(qū)域，該域逐步逼近滲流實(shí)際區(qū)域，所得流場(chǎng)逐步逼近真實(shí)解。

(1)對(duì)土石壩進(jìn)行全域規(guī)則劃分。首先對(duì)壩體進(jìn)行豎向和橫向剖分，得到規(guī)則的四邊形單元，再對(duì)各四邊形單元進(jìn)行同方向?qū)蔷€連接，剖分為三角形單元，節(jié)點(diǎn)編號(hào)按照由左向右的順序，由上而下依次編號(hào)。

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

(4)

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

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

(5)

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

進(jìn)行曲線擬合時(shí)，縱向水頭差值不宜過(guò)大，以確保擬合曲線為順次遞減曲線，要求在劃分網(wǎng)格時(shí)，單元適當(dāng)?shù)膭澐侄嘁恍?，即可確保水頭求解的精確性，又可確保擬合的曲線更加符合自由面形狀。為確保擬合曲線的準(zhǔn)確性，本文采用基于全域總勢(shì)能確定溢出點(diǎn)，此方法可以較為準(zhǔn)確地確定溢出點(diǎn)。

表1 解析解與數(shù)值解自由面位置數(shù)據(jù)對(duì)比 m

圖1 均質(zhì)矩形壩網(wǎng)格剖分示意

圖2 自由面線形示意

3 算例

3.1 有解析解矩形均質(zhì)土壩

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

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

3.2 有解析解的均質(zhì)梯形壩

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

圖3 均質(zhì)梯形壩大網(wǎng)格剖分

圖4 均質(zhì)梯形壩大網(wǎng)格剖分下三種自由面的比較

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

4 結(jié) 論

滲流在壩坡穩(wěn)定及壩坡變形中起關(guān)鍵性作用，滲流自由面的確定是滲流計(jì)算的關(guān)鍵內(nèi)容之一。本文提出了一種曲線擬合法，并編制有限元程序，應(yīng)用于計(jì)算模型及工程實(shí)例中得到如下結(jié)論：

表2 兩種數(shù)值方法計(jì)算不同網(wǎng)格劃分自由面位置數(shù)據(jù)對(duì)比 m

(1)本文在全域求解的基礎(chǔ)上，應(yīng)用最小二乘法理論對(duì)全域求解后的自由面節(jié)點(diǎn)進(jìn)行逐次遞減的高次曲線擬合，無(wú)需多次迭代，簡(jiǎn)化了求解過(guò)程。

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

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

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

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

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

[4]朱軍， 劉光廷. 改進(jìn)的單元滲透矩陣調(diào)整法求解無(wú)壓滲流場(chǎng)[J]. 水利學(xué)報(bào)， 2001(8)： 49-52．

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

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

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

[8]沈振中， 陳小虎， 吳越建. 求解堤壩滲流場(chǎng)的罰函數(shù)無(wú)單元法[J].河海大學(xué)學(xué)報(bào)：自然科學(xué)版，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]吳夢(mèng)喜， 張雪勤. 有自由面滲流分析的虛單元法[J]. 水利學(xué)報(bào)， 1994(8)： 64-71．

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

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

(責(zé)任編輯 王 琪)

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

山東省自然科學(xué)基金資助項(xiàng)目(ZR2012EEL31)

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

TV139.14

A

0559-9342(2016)11-0050-04