基于總勢(shì)能極小原理的滲流溢出點(diǎn)位置計(jì)算方法

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

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

?

基于總勢(shì)能極小原理的滲流溢出點(diǎn)位置計(jì)算方法

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

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

在考慮滲流影響的工程中，自由面的計(jì)算是滲流計(jì)算的主要內(nèi)容,精確確定溢出點(diǎn)位置是準(zhǔn)確計(jì)算自由面的關(guān)鍵?；诳倓?shì)能極小原理的滲流溢出點(diǎn)位置計(jì)算方法，是在溢出點(diǎn)位置調(diào)整的過程中，通過尋找全域總勢(shì)能的極小值確定溢出點(diǎn)位置的一種方法。通過算例分析，得出如下結(jié)論：基于總勢(shì)能極小原理的滲流溢出點(diǎn)位置計(jì)算方法彌補(bǔ)了先前溢出點(diǎn)計(jì)算方法精度低、需迭代計(jì)算、不易收斂的缺點(diǎn)，精確、高效地計(jì)算出溢出點(diǎn)位置。

溢出點(diǎn)；滲流自由面；總勢(shì)能極小值

滲流計(jì)算在土石壩穩(wěn)定性[1]、基坑穩(wěn)定[2]計(jì)算中占有重要地位，其中自由面的確定是滲流計(jì)算的主要內(nèi)容之一。以土石壩穩(wěn)定性計(jì)算為例，通常以壩體的自由面為分割線，自由面以下為飽和區(qū)，取浮重度計(jì)算；自由面以上為非飽和區(qū)，取天然重度計(jì)算。

確定溢出點(diǎn)位置是計(jì)算自由面的關(guān)鍵。在滲流場(chǎng)中，上游最高水位點(diǎn)為自由面的起始點(diǎn)，溢出點(diǎn)為自由面的終止點(diǎn)。上游水位已知，若能準(zhǔn)確確定溢出點(diǎn)位置，無論在應(yīng)用無網(wǎng)格法[3-5]、虛單元法[6-8]、初流量法[9-11]、單元滲透矩陣調(diào)整法[12-14]、或其他方法[15-19]求解滲流自由面時(shí)，均可更準(zhǔn)確地計(jì)算自由面位置。

目前，已經(jīng)提出的確定溢出點(diǎn)的方法主要有沿坡面滑動(dòng)法和二次曲線相交法。沿坡面滑動(dòng)法是將溢出點(diǎn)和滲流自由面上的節(jié)點(diǎn)一同進(jìn)行迭代計(jì)算，沿坡面滑動(dòng)調(diào)整其位置[17]。然而此方法溢出點(diǎn)需要與自由面一同迭代，計(jì)算量大且會(huì)存在不收斂的情況；滲流自由面形狀一般可用二次曲線來描述，二次曲線相交法通過確定自由面上三點(diǎn)做二次曲線，與坡面的交點(diǎn)即為溢出點(diǎn)[17]。但此方法中對(duì)三點(diǎn)的精度要求太高，誤差往往比較大。

針對(duì)上述方法的不足，學(xué)者們已經(jīng)做了一些改進(jìn)與研究。吳夢(mèng)喜等[6]提出了在可能的溢出面范圍內(nèi)，先假設(shè)一排較低的點(diǎn)為溢出點(diǎn)。溢出點(diǎn)以下結(jié)點(diǎn)為已知節(jié)點(diǎn)，以上結(jié)點(diǎn)為未知節(jié)點(diǎn)，計(jì)算區(qū)域各節(jié)點(diǎn)的勢(shì)，反復(fù)迭代直到收斂為止。比較溢出點(diǎn)上部一點(diǎn)的勢(shì)與其位勢(shì)大小，若勢(shì)小于位勢(shì)，則假定的溢出點(diǎn)位置合理，否則應(yīng)把其上一點(diǎn)作為溢出點(diǎn)，重新計(jì)算直至滿足要求為止的方法。此方法為局部變網(wǎng)格法，需將自由面與溢出點(diǎn)共同迭代，計(jì)算量大。朱軍等[13]指出，可將溢出點(diǎn)邊界作為第二類邊界條件，在求出流量后，通過調(diào)整自由項(xiàng)，將其轉(zhuǎn)化為第一類邊界條件迭代計(jì)算。此方法受起始條件的影響較大，且需反復(fù)迭代，計(jì)算量大。黃蔚等[15]提出在處理溢出段時(shí)，根據(jù)水頭與高程的關(guān)系將溢出面上的點(diǎn)定義為固定點(diǎn)與活動(dòng)點(diǎn)，通過迭代與線性插值最終找到溢出點(diǎn)位置。此方法也是將自由面與溢出點(diǎn)一同迭代，需重復(fù)構(gòu)造總滲透矩陣與自由項(xiàng)。謝國海等[18]在等效滲透系數(shù)法中，根據(jù)水頭與高程的差值，將溢出點(diǎn)與自由面一起迭代，直至滿足水頭與高程的差的絕對(duì)值小于收斂值。此方法將溢出點(diǎn)與自由面一起迭代，收斂判別條件復(fù)雜，有時(shí)可能不收斂。

本文從全域總勢(shì)能的概念入手，提出了一種確定溢出點(diǎn)位置的新方法，不進(jìn)行迭代，也不重新構(gòu)造總滲透矩陣，依據(jù)全域總勢(shì)能的極小值即可精確確定溢出點(diǎn)位置，并通過與有解析解的算例對(duì)比驗(yàn)證了該方法的精度。

1　滲流的基本理論

1.1支配方程和邊界條件

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

(1)

其中kx，ky為x，y方向的滲透系數(shù)，h為水頭函數(shù)。

第一類邊界條件又稱水頭邊界條件，為邊界上給定的位勢(shì)函數(shù)或水頭分布函數(shù)，可寫為：

(2)

在穩(wěn)定滲流場(chǎng)中，上、下游的已知水頭值，就屬于此類邊界條件。

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

(3)

1.2應(yīng)用有限單元法求解滲流場(chǎng)

在上述邊界條件下，依據(jù)式(1)求解水頭函數(shù)是十分復(fù)雜的，對(duì)邊界的要求十分苛刻。故常采用對(duì)泛函求極值的方法求解微分方程，構(gòu)造如下泛函：

(4)

將全域單元?jiǎng)澐趾螅诿總€(gè)單元中對(duì)節(jié)點(diǎn)水頭求導(dǎo)，得到單元滲透矩陣，組裝為總滲透矩陣，得到求解水頭的線性方程組：

[K]{h}={f}

(5)

[K]為總滲透矩陣，{f}為自由項(xiàng)，當(dāng)沒有流量流入與流出時(shí)，{f}中的元素為0。求解式(5)即可得到滲流場(chǎng)內(nèi)各節(jié)點(diǎn)的水頭值。

2　基于全域總勢(shì)能確定溢出點(diǎn)

2.1全域總勢(shì)能

在二維滲流場(chǎng)中，以上游坡面底部為坐標(biāo)原點(diǎn)，水平方向向左為x軸正方向，鉛垂向上方向?yàn)閥軸正方向，建立直角坐標(biāo)系，對(duì)全域劃分有限單元，定義每延米上的單元總勢(shì)能為：

(6)

其中：Ee為每延米上的單元總勢(shì)能，J/m。ρ為水的密度；g為重力加速度；Ni為單元形函數(shù)；n為單元節(jié)點(diǎn)數(shù)。

每延米上全域總勢(shì)能為：

E=∑Ee

(7)

在式(6)中，ρ、g均為常數(shù)，故單元總勢(shì)能值的變化受積分式的影響，而積分式在數(shù)值上等于單元中所有點(diǎn)水頭值的和。故單元總勢(shì)能受單元中水頭總值的影響，全域總勢(shì)能受全域水頭總值的影響。

2.2基于總勢(shì)能極小原理的求解溢出點(diǎn)方法

圖1給定真實(shí)溢出點(diǎn)情況下的邊界分布情況

由于真實(shí)溢出點(diǎn)的位置是未知的，計(jì)算時(shí)需要對(duì)溢出點(diǎn)的位置進(jìn)行設(shè)置，即將h=y的邊界條件賦予可能的溢出段。當(dāng)設(shè)置的溢出點(diǎn)低于真實(shí)的溢出點(diǎn)時(shí)，部分真實(shí)溢出段邊界將被定義為不透水邊界(如圖2中FF0所示)，這阻止了部分溢出段邊界處的流量外泄，導(dǎo)致此時(shí)的全域水頭總值比真實(shí)情況的全域水頭總值高，使得此時(shí)全域總勢(shì)能比真實(shí)情況的全域總勢(shì)能高；當(dāng)設(shè)置的溢出點(diǎn)高于真實(shí)的溢出點(diǎn)時(shí)，原本位于BF邊界之上，應(yīng)滿足h

圖2　給定溢出點(diǎn)低于真實(shí)的溢出點(diǎn)情況下的邊界分布情況

圖3給定溢出點(diǎn)高于真實(shí)的溢出點(diǎn)情況下的邊界分布情況

因此，全域總勢(shì)能在賦予溢出點(diǎn)位置升高的過程中，會(huì)呈現(xiàn)先減小后增大的變化過程，當(dāng)全域總勢(shì)能取極小值時(shí)，對(duì)應(yīng)的溢出點(diǎn)位置就是真實(shí)的溢出點(diǎn)位置(如圖4中F點(diǎn)所示)。

圖4全域總勢(shì)能與溢出點(diǎn)位置關(guān)系示意圖

2.3計(jì)算單元的選取

根據(jù)上述全域總勢(shì)能計(jì)算方法及溢出點(diǎn)求解方式，進(jìn)行有限元數(shù)值計(jì)算時(shí)需要注意以下幾點(diǎn)：

(1)在計(jì)算全域總勢(shì)能時(shí)，需在求出單元的水頭函數(shù)后，將水頭函數(shù)對(duì)單元面積積分，得到單元的總勢(shì)能，再對(duì)所有單元的總勢(shì)能求和。為提高計(jì)算效率，推薦采用水頭函數(shù)及單元形函數(shù)較為簡(jiǎn)單的單元類型；

(2) 為溢出段及上、下游水位賦予第一類邊界條件時(shí)，應(yīng)盡可能的準(zhǔn)確。對(duì)于直線邊界，邊界上水頭函數(shù)為線性時(shí)最佳。

綜上所述，本文選取三角形單元作為計(jì)算區(qū)域劃分的基本單元。水頭函數(shù)的表達(dá)式為

h(x,y)=α1+α2x+α3y

(8)

α1，α2，α3為與單元相關(guān)的常數(shù)。

此類單元可以很好地模擬第一類邊界條件，水頭函數(shù)對(duì)面積的積分式也較為簡(jiǎn)單。

2.4計(jì)算步驟

根據(jù)全域總勢(shì)能確定溢出點(diǎn)位置的程序框圖如圖5所示。

圖5溢出點(diǎn)計(jì)算程序框圖

下面以斜坡溢出面邊界為例，具體解釋精確計(jì)算溢出點(diǎn)位置的程序步驟：

(1) 全域劃分有限單元，釋放上、下游水位邊界條件，依據(jù)式(5)，求解各節(jié)點(diǎn)水頭值，再依據(jù)式(6)、式(7)求解在沒有賦予溢出點(diǎn)情況下的全域總勢(shì)能E0。

(2) 將溢出點(diǎn)升高一個(gè)節(jié)點(diǎn)，依據(jù)式(5)、式(6)、式(7)求解全域總勢(shì)能Ep，p為執(zhí)行此步驟的序號(hào)。

(3) 判斷Ep與Ep-1的大小，若Ep>Ep-1，執(zhí)行下一步；若Ep-1>Ep，回到步驟(2)。

(4) 設(shè)此時(shí)的p值為m，當(dāng)p=m-2，m-1，m時(shí)，溢出點(diǎn)的坐標(biāo)分別為(x1，y1)，(x2，y2)，(x3，y3)，令E10=Em-2。

(5) 設(shè)最小誤差為α，令y2=y1+qα，x2=x1+(x3-x1)qα/(y3-y1)，根據(jù)(x2,y2)的新坐標(biāo)，微調(diào)總滲透矩陣[K]，并將(x2，y2)作為新的溢出點(diǎn)，依據(jù)式(5)、式(6)、式(7)求解全域總勢(shì)能E1q，q為執(zhí)行此步驟的序號(hào)。

(6) 判斷E1q與E1q-1的大小，若E1q>E1q-1，執(zhí)行下一步；若E1q-1>E1q，回到步驟(5)。

(7) 取y=y2-α，x=x2-(x3-x1)α/(y3-y1)，(x，y)即為所求的溢出點(diǎn)坐標(biāo)。

3　算　例

3.1算例1

以10m×10m均質(zhì)土壩為例進(jìn)行分析，土壩上下游水位分別為10m和2m，底部為不透水邊界，該問題的自由面解析表達(dá)式為y2=100-8x。其溢出點(diǎn)位置的解析解為4.472m。采取兩種不同的單元剖分進(jìn)行計(jì)算，計(jì)算模型1共200個(gè)單元，以直角邊長(zhǎng)為1m的三角形進(jìn)行劃分，如圖6(a)；計(jì)算模型2共50個(gè)單元，以直角邊長(zhǎng)為2m的三角形進(jìn)行劃分，如圖6(b)。應(yīng)用上文所述方法求解溢出點(diǎn)位置。全域總勢(shì)能的變化情況如表1所示。

圖6　算例1網(wǎng)格剖分圖

由表1知，在計(jì)算模型1中，求出的溢出點(diǎn)位置為4.55m，與解析解相差0.08m，相對(duì)誤差為1.74%；在計(jì)算模型2中，求出的溢出點(diǎn)位置為4.35m，與解析解相差-0.12m，相對(duì)誤差為-2.73%，滿足規(guī)范的誤差要求，在大網(wǎng)格劃分下，仍然能得到接近實(shí)際位置的溢出點(diǎn)。

由表2知，同樣采取計(jì)算模型1時(shí)，本文方法算得的溢出點(diǎn)，精度比節(jié)點(diǎn)虛流量法、初流量法、改進(jìn)的初流量法、改進(jìn)的截至負(fù)壓法、改進(jìn)的丟單元法等方法都高。

表2　溢出點(diǎn)精度對(duì)比

3.2算例2

為了進(jìn)一步驗(yàn)證本文方法，選取帶有坡面的土石壩壩體進(jìn)行溢出點(diǎn)位置的計(jì)算，壩體上、下游水位分別為15m和4m，上、下游坡比分別為1∶3和1∶2，最高處高17m，底部橫向?qū)挾葹?7m，底部為不透水邊界，解析解求得的溢出點(diǎn)的位置為6.23m。

如圖7所示，分別以高為1m的三角形與高為2m的三角形對(duì)此壩體進(jìn)行單元剖分，應(yīng)用本文方法編制有限元程序，對(duì)溢出點(diǎn)位置進(jìn)行計(jì)算，全域總勢(shì)能如表3所示。

圖7　算例2網(wǎng)格剖分圖

由表3知，1m網(wǎng)格劃分下，求出溢出點(diǎn)位置為6.40m，與解析解的差為0.17m，相對(duì)誤差為2.72%；在2m網(wǎng)格劃分下，求出溢出點(diǎn)位置為6.00m，與解析解的差為-0.23m，相對(duì)誤差為-3.69%，滿足工程精度要求。

4　結(jié)　論

本文提出了一種以全域總勢(shì)能極小值精確計(jì)算溢出點(diǎn)的方法，應(yīng)用有限元程序求解了二維計(jì)算模型與土石壩壩體的溢出點(diǎn)位置并與其他方法作了對(duì)比，主要結(jié)論如下：

(1) 應(yīng)用本文基于全域總勢(shì)能計(jì)算溢出點(diǎn)的方法，相比以前的一些方法，可更精確地確定溢出點(diǎn)位置。

(2) 本文方法避免了同類研究由于將溢出點(diǎn)與自由面一同迭代，結(jié)果易不收斂的問題；在求解溢出點(diǎn)時(shí)，不必重新構(gòu)造總滲透矩陣，減少計(jì)算量，提高計(jì)算效率。

(3) 在較大單元網(wǎng)格劃分中，本文方法仍然具有很高的精度。

[1]劉娟奇，王志強(qiáng)，梁收運(yùn).庫水位下降對(duì)新集水庫均質(zhì)土壩滲流及穩(wěn)定性影響分析[J]. 水利與建筑工程學(xué)報(bào),2014,12(6):38-43.

[2]位俊俊，張利偉，孔德志.基坑滲流穩(wěn)定分析[J].水利與建筑工程學(xué)報(bào),2012,10(3):79-87.

[3]DarbaniM,OuahsineA,VillonP,etal.Meshlessmethodforshallowwaterequationswithfreesurfaceflow[J].AppliedMathematicsandComputation, 2011,217(11):5113-5124.

[4]JieYX,LiuWJ,LiGX.ApplicationofNEMinseepageanalysiswithafreesurface[J].MathematicsandComputersinSimulation, 2013,89(2):23-37.

[5]TangJing,LiuYongbiao.Penaltyfunctionelementfreemethodtosolvecomplexseepagefieldofearthfilldam[J].IERIProcedia, 2012(1):117-123.

[6]吳夢(mèng)喜,張雪勤.有自由面滲流分析的虛單元法[J].水利學(xué)報(bào),1994(8):67-71.

[7]凌道盛.有自由面滲流分析的虛節(jié)點(diǎn)法[J].浙江大學(xué)學(xué)報(bào)(工學(xué)版),2002,36(3):243-246.

[8]崔皓東,朱岳明.有自由面滲流分析的改進(jìn)節(jié)點(diǎn)虛流量全域迭代法[J].武漢理工大學(xué)學(xué)報(bào)(交通科學(xué)與工程版),2009,33(2):238-241.

[9]張有天,陳平,王鐳.有自由面滲流分析的初流量法[J].水利學(xué)報(bào),1988(8):18-26.

[10]王媛.求解有自由面滲流問題的初流量法的改進(jìn)[J].水利學(xué)報(bào),1998(3):68-73.[11]潘樹來,王全鳳,俞 縉.利用初流量法分析有自由面滲流問題之改進(jìn)[J].巖土工程學(xué)報(bào),2012,34(2):202-209.[12]王賢能.有自由面滲流分析的高斯點(diǎn)法[J].水文地質(zhì)工程地質(zhì),1997(6):1-4.

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

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

[15]梁業(yè)國,熊文林,周創(chuàng)兵.有自由面滲流分析的子單元法[J].水利學(xué)報(bào),1997(8):34-38.

[16]王均星,吳雅峰,白呈富.有自由面滲流分析的流行單元法[J].水電能源科學(xué),2003,21(4):23-26.

[17]毛昶熙,段祥寶,李祖貽，等.滲流數(shù)值計(jì)算與程序應(yīng)用[M].南京:河海大學(xué)出版社,1999.

[18]謝國海,章廣成.二維無壓穩(wěn)定滲流分析的改進(jìn)方法[J].地球與環(huán)境,2005,33(S1):52-56.

[19]張宜虎.采用等效滲透系數(shù)法搜索滲流自由面[J].工程地質(zhì)學(xué)報(bào),2004,12(2):177-192.

Calculation Method of Seepage Overflow Point Based on the Principle of Minimum Total Potential Energy

LI Yuandong, HOU Xingmin, ZHENG Shanshan, SUN Weijian

(SchoolofCivilEngineering,YantaiUniversity,Yantai,Shandong264005,China)

In engineering that needs to consider the effect of seepage, determining the overflow point precisely is the key to the calculation of free surface, which is the main content of seepage calculation. Based on the principle of minimum total potential energy to calculate the overflow point is a method that the location of overflow point is determined through calculating the minimum of global total potential energy by adjusting its position. Example analysis showed that this method can calculate the location of overflow point precisely and efficiently, which overcomes the shortcomings of previous methods such as low accuracy, requiring iterative calculation and difficult to converge.

overflow point; seepage free surface; minimum total potential energy

10.3969/j.issn.1672-1144.2016.04.017

2016-04-17

2016-05-14

國家自然科學(xué)基金項(xiàng)目(51479174)

李遠(yuǎn)東(1990—)，男，山東青島人，碩士研究生，研究方向?yàn)榻ㄖc土木工程。 E-mail：wddd2013@163.com

侯興民(1970—)，男，山東壽光人，博士，教授，主要從事巖土工程方面教學(xué)與科研工作。 E-mail：houxm@ytu.edu.cn

TV139.14

A

1672—1144(2016)04—0084—05