姚 震

(甘肅省水利水電勘測(cè)設(shè)計(jì)研究院有限責(zé)任公司，蘭州 730030)

1 概 述

土石壩的滲流計(jì)算主要是為了確定壩體浸潤(rùn)線的位置，為壩體的穩(wěn)定分析和布置觀測(cè)設(shè)備提供依據(jù)；同時(shí)確定壩體與壩基的滲透流量，以便估算水庫(kù)的滲漏損失，而且還要確定壩體和壩體滲流區(qū)的滲透坡降，檢查產(chǎn)生滲透變形的可能性，以便采取適當(dāng)?shù)目刂拼胧Ｔ趯?shí)用上，均質(zhì)土壩滲流常采用分段法進(jìn)行計(jì)算，而分段法需要求解非線性方程。Newton迭代法是解非線性方程最常用的方法，適用范圍廣泛，步驟簡(jiǎn)單。但如同所有的迭代法一樣，它需要大量計(jì)算。MATLAB是一種強(qiáng)大的計(jì)算工具，利用MATLAB來(lái)實(shí)施Newton迭代法的計(jì)算可以大大節(jié)省時(shí)間[1]。本文以水平不透水層上均質(zhì)土壩進(jìn)行滲流計(jì)算。

2 滲流分析的水力學(xué)方法原理

2.1 水力學(xué)方法采用的基本假定

1) 壩體的材料為均質(zhì)，且各點(diǎn)在各方向的滲透系數(shù)k相同。

2) 滲流為層流，認(rèn)為滲透流速符合達(dá)西定律v=kJ。

3) 滲透水流是漸變流，認(rèn)為任意過(guò)水?dāng)嗝嫔狭魉賤相同。

2.2 計(jì)算簡(jiǎn)圖

見圖1。

圖1 水平不透水層上均質(zhì)土壩剖面圖

圖1中，H為壩高，m；b為壩頂寬，m；H1為壩前水深，m；H2為下游水深，m；hk為逸出點(diǎn)水深，m；a0為壩體下游逸出點(diǎn)距下游水面高度，m；ΔL為等效矩形體寬度，m；m1為上游壩面邊坡系數(shù)；m2為下游壩面邊坡系數(shù)；k為壩身的滲透系數(shù)，m/s。

2.3 水力學(xué)計(jì)算原理及過(guò)程[2]

在實(shí)用上，土石壩滲流計(jì)算常用分段法，并且又分為三段法和兩段法兩種。三段法是由巴甫洛夫斯基提出的，是將壩體內(nèi)滲流區(qū)劃分為3段，第一段為上游楔形段，第二段為中間段，第三段為下游楔形段。對(duì)每一段應(yīng)用漸變流基本公式建立流量表達(dá)式，然后通過(guò)3段的聯(lián)合求解，即可確定土石壩滲流量及逸出點(diǎn)水深，并可繪出浸潤(rùn)線。兩段法是在三段法的基礎(chǔ)上簡(jiǎn)化而來(lái)的，將上游楔形段和中間段合并，把土石壩滲流區(qū)劃分成上游段和下游段兩端。在下邊的計(jì)算中，將用兩段法來(lái)分析土石壩滲流。

在兩段法中，把上游楔形段用假想的等效矩形體代替，即認(rèn)為水流從垂直面滲入壩體，而矩形體寬度的確定，應(yīng)使在相同的上游水深和單寬流量的作用下，分別通過(guò)矩形體和楔形體的水頭損失相等?，F(xiàn)用兩斷法來(lái)進(jìn)行分析，根據(jù)實(shí)驗(yàn)研究，等效矩形體的寬度可由下式確定：

(1)

式中：m1為壩上游面的邊坡系數(shù)。

2.3.1 上游段的計(jì)算

滲流從過(guò)水?dāng)嗝鍭′B′到CG的水頭差ΔH=H1-hk，兩過(guò)水?dāng)嗝嬷g平均滲透路程Δs=L+ΔL-m2hk。其中L=(H-H1)m1+b+Hm2，m2為壩下游坡面的邊坡系數(shù)，故上游段的平均水力坡度為：

(2)

根據(jù)杜比公式，上游段的平均滲透流速為：

(3)

(4)

很顯然，要用式(4)計(jì)算滲流量還不可能，因?yàn)橐莩鳇c(diǎn)水深hk是未知數(shù)，所以還必須要對(duì)下游建立計(jì)算公式，以便和式(4)進(jìn)行聯(lián)解。

2.3.2 下游段的計(jì)算

當(dāng)下游水深H2≠0時(shí)，應(yīng)將該段分為Ⅰ、Ⅱ兩部分。第Ⅰ部分位于下游水面以上，為無(wú)壓滲流；第Ⅱ部分位于下游水面以下，為有壓滲流，近似認(rèn)為下游段內(nèi)流線為水平線。

下游段單寬總滲流量為：

(5)

式(5)中a0=hk-H2，聯(lián)解方程式(4)、式(5)，可求得土石壩單寬滲流量q及逸出點(diǎn)高度a0。

2.3.3 浸潤(rùn)線計(jì)算原理及過(guò)程

因用等效的矩形體取代上游三角楔形體，A′B′即為上游的入滲起始斷面，今取以G為坐標(biāo)原點(diǎn)的一組直角坐標(biāo)系來(lái)研究浸潤(rùn)曲線的計(jì)算，X軸以向左為正。浸潤(rùn)曲線的方程式為：

(6)

假定一系列x值，可由式(6)算得一系列相應(yīng)的y值，點(diǎn)繪成浸潤(rùn)線。由式(3)所繪出的浸潤(rùn)曲線，其上游端是從A′點(diǎn)開始的，而實(shí)際上入滲點(diǎn)應(yīng)在A，故曲線的前端A′F應(yīng)加以修正。在實(shí)用上，常采用近似方法來(lái)修正它，即把A點(diǎn)作為曲線的上游端起點(diǎn)，再選擇(用曲線板)適當(dāng)能與后半段曲線光滑連接的曲線AF去代替A′F即可。見圖1。

3 牛頓(Newton)迭代法原理[3]

Newton法是求解非線性方程最常用的方法，它是根據(jù)對(duì)非線性方程逐步線性化建立起來(lái)的一種迭代格式。

設(shè)非線性方程f(x)=0，且f連續(xù)可微，x*是方程的實(shí)根，x(0)是它的一個(gè)近似值。在x(0)附近將f線性化，即用f(x)在x(0)的一階Taylor公式代替f(x)，得

f(x)≈f(x(0))+f′(x(0))(x-x(0))

(7)

則方程近似地表示為線性方程：

f(x(0))+f′(x(0))(x-x(0))=0

(8)

只要f′(x(0))≠0，可得新的近似根，記為x(1)，則：

(9)

一般地，設(shè)x(k)已經(jīng)求出，在x(k)處將方程線性化：

f(x(k))+f′(x(k))(x-x(k))=0

(10)

若f′(x(k))≠0，式(10)有唯一解，記為x(k+1)，則得迭代格式：

(11)

式(11)稱為Newton迭代格式。

Newton迭代法有明顯的幾何意義。因?yàn)榉匠蘤(x)=0的實(shí)根是曲線y=f(x)和X軸焦點(diǎn)的橫坐標(biāo)，而x(k+1)是曲線y=f(x)在點(diǎn)Pk(x(k)，f(x(k)))處的切線與X軸焦點(diǎn)的橫坐標(biāo)，因此迭代過(guò)程是逐次以切線代替曲線，用切線與X軸交點(diǎn)的橫坐標(biāo)x(k)逼近曲線與X軸交點(diǎn)的橫坐標(biāo)x*的過(guò)程。見圖2。所以Newton法又稱為切線法。

圖2 Newton迭代法求根圖示

4 用MATLAB實(shí)現(xiàn)某均質(zhì)土壩的滲流計(jì)算

4.1 工程概況

某均質(zhì)土壩建于不透水地基上，已知壩高H為17 m，上游水深H1為15 m，下游水深H2為2 m，上游邊坡系數(shù)m1為3，下游邊坡系數(shù)m2為2，壩頂寬b為6 m，壩身土的滲透系數(shù)k經(jīng)實(shí)驗(yàn)測(cè)得為0.001 cm/s。

4.2 用MATLAB實(shí)現(xiàn)滲流計(jì)算

%水平不透水層上均質(zhì)土壩滲流計(jì)算

%1、輸入已知參數(shù)

H=17;

b=6;

H1=15;

H2=2;

m1=3;

m2=2;

k=0.00001;

%2、水力學(xué)法滲流計(jì)算過(guò)程

%2.1、上游段

%等效矩形體寬度

ΔL=(m1./(1+2.*m1)).*H1；

L=(H-H1).*m1+b+H.*m2；

%(1)q=k.*((H1.^2-(H2+a0).^2)./(2.*(L+ΔL-m2.*(H2+a0))))；

%2.2、下游段

%(2)q=(k.*a0./m2).*(1+2.3.*log10((H2+a0)./a0))；

%聯(lián)立(1)、(2)可得下式

% k.*((H1.^2-(H2+a0).^2)./(2.*(L+ΔL-m2.*(H2+a0))))…

% =(k.*a0./m2).*(1+2.3.*log10((H2+a0)./a0))

%3、用牛頓迭代法求解滲流量

%3.1、m文件定義函數(shù)

function y=fun1(a0，k，H1，H2，L，ΔL，m2)

y= k.*((H1.^2-(H2+a0).^2)./(2.*(L+ΔL-m2.*(H2+a0))))…

-(k.*a0./m2).*(1+2.3.*log10((H2+a0)./a0))；

end

function Y=fun2(a0，k，H1，H2，L，ΔL，m2)

Y=a0-fun1(a0，k，H1，H2，L，ΔL，m2)./ diff(fun1(a0，k，H1，H2，L，ΔL，m2))；

end

%3.2、選定初始值

a0=0:0.000001:5；

y=fun1(a0，k，H1，H2，L，ΔL，m2)；

plot(a0，y)；

grid on；

%3.3、求出迭代函數(shù)

syms a0

%Y=a0-y./ diff(y)；

Y=fun2(a0，k，H1，H2，L，ΔL，m2)；

%3.4、實(shí)施迭代

%等價(jià)方程：a0=Y

a0=3；

i=1;

while i<=10

a0+1=subs(Y，a0)；

if (abs(a0+1-a0)>1e-5)

a0=a0+1；

else break

end

i=i+1；

end

fprintf(’a0=%f’，a0)

%3.5、求出單寬滲流量

q=(k.*a0./m2).*(1+2.3.*log10((H2+a0)./a0))

%4、繪制浸潤(rùn)線

x=[0：5：L+ΔL-m2.*(a0+H2)，L+ΔL-m2.*(a0+H2)]

y=sqrt((x./(L+ΔL-m2.*(H2+a0))).*(H1.^2-(H2+a0).^2)+(H2+a0).^2)

運(yùn)行程序，輸出結(jié)果：

a0=3.162856

q=2.36e-005

x=0 5.0000 10.0000 15.0000 20.0000 25.0000 30.0000 35.0000 40.0000 42.1029

y=5.1629 7.0859 8.5886 9.8651 10.9943 12.0179 12.9609 13.8398 14.6661 15.0000

5 與理正巖土滲流分析軟件計(jì)算成果的比較分析及建議

5.1 理正巖土滲流分析軟件計(jì)算成果

將上述工程計(jì)算參數(shù)輸入理正巖土滲流分析軟件滲流問(wèn)題公式法模塊，計(jì)算結(jié)果整理如下：

下游出逸點(diǎn)高度：a0=3.993(m)

單位寬度滲流量：q=2.34×10-5(m3/s·m)

5.2 計(jì)算成果比較分析及建議

比較本文水力學(xué)計(jì)算成果和理正巖土滲流分析軟件計(jì)算成果，前者下游出逸點(diǎn)高度為3.163 m，后者為3.993 m；前者單位寬度滲流量為2.36×10-5m3/s·m，后者為2.34×10-5m3/s·m。造成兩種計(jì)算方法下游出逸點(diǎn)高度差值較大的原因是：本文采用水力學(xué)滲流計(jì)算公式，理正巖土滲流分析軟件依據(jù)《堤防工程設(shè)計(jì)規(guī)范》(GB 50286-2013)不透水堤基均質(zhì)土堤滲流計(jì)算公式[7]，兩種公式計(jì)算結(jié)果有差異。

依據(jù)《碾壓式土石壩設(shè)計(jì)規(guī)范》(SL 274-2001)，貼坡排水頂部高程應(yīng)高于壩體浸潤(rùn)線出逸點(diǎn)[8]。理正巖土滲流分析軟件計(jì)算的出逸點(diǎn)高度較本文水力學(xué)滲流計(jì)算公式高，相應(yīng)貼坡排水頂部高程也較高，理正計(jì)算結(jié)果對(duì)于工程來(lái)說(shuō)更加保守。

6 結(jié) 語(yǔ)

本文結(jié)合實(shí)例，闡述了水平不透水層上均質(zhì)土壩的滲流計(jì)算，可以看出利用MATLAB來(lái)實(shí)施Newton迭代法可以大大節(jié)省時(shí)間。MATLAB滲流計(jì)算程序具有通用性，對(duì)其它工程均質(zhì)土壩滲流計(jì)算具有參考意義。通過(guò)分析比較本文水力學(xué)滲流計(jì)算成果和理正巖土滲流分析軟件計(jì)算成果可知，理正計(jì)算結(jié)果對(duì)于工程來(lái)說(shuō)更加保守。