姚 震

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

1 概 述

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

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

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

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

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

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

2.2 計算簡圖

見圖1。

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

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

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

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

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

(1)

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

2.3.1 上游段的計算

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

(2)

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

(3)

(4)

很顯然，要用式(4)計算滲流量還不可能，因為逸出點水深hk是未知數(shù)，所以還必須要對下游建立計算公式，以便和式(4)進(jìn)行聯(lián)解。

2.3.2 下游段的計算

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

下游段單寬總滲流量為：

(5)

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

2.3.3 浸潤線計算原理及過程

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

(6)

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

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

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

設(shè)非線性方程f(x)=0，且f連續(xù)可微，x*是方程的實根，x(0)是它的一個近似值。在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迭代法有明顯的幾何意義。因為方程f(x)=0的實根是曲線y=f(x)和X軸焦點的橫坐標(biāo)，而x(k+1)是曲線y=f(x)在點Pk(x(k)，f(x(k)))處的切線與X軸焦點的橫坐標(biāo)，因此迭代過程是逐次以切線代替曲線，用切線與X軸交點的橫坐標(biāo)x(k)逼近曲線與X軸交點的橫坐標(biāo)x*的過程。見圖2。所以Newton法又稱為切線法。

圖2 Newton迭代法求根圖示

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

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)實驗測得為0.001 cm/s。

4.2 用MATLAB實現(xiàn)滲流計算

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

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

H=17;

b=6;

H1=15;

H2=2;

m1=3;

m2=2;

k=0.00001;

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

%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、實施迭代

%等價方程：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、繪制浸潤線

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)

運行程序，輸出結(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 與理正巖土滲流分析軟件計算成果的比較分析及建議

5.1 理正巖土滲流分析軟件計算成果

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

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

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

5.2 計算成果比較分析及建議

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

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

6 結(jié) 語

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