溫斯鈞，張浩，安航永

(1.華北水利水電大學(xué)，河南 鄭州 450045；2. 濮陽黃河河務(wù)局范縣黃河河務(wù)局，河南 濮陽 457506；3.廣州禺山水務(wù)勘測設(shè)計股份有限公司，廣東 廣州 511400)

0 引言

大斷面水力要素計算是水文與水力學(xué)工作中經(jīng)常遇到的問題，其目的是根據(jù)已知的大斷面實測數(shù)據(jù)，計算特定水位所對應(yīng)的過水面積、濕周、水力半徑、水面寬等水力要素。傳統(tǒng)手算法求解斷面水力要素費時、費力，且當(dāng)需要計算多個斷面多個水位對應(yīng)的水力要素時大量的重復(fù)性工作極易出錯，計算精度難以保證。MATLAB強大的插值計算和數(shù)據(jù)處理能力將傳統(tǒng)水力要素計算方法程序化，使其能夠更加高效、準(zhǔn)確地解決大斷面水力要素計算問題。

1 基本原理

實測大斷面示意圖如圖1所示，計算水位為h。

圖1 實測大斷面示意圖

水力要素計算步驟如下：根據(jù)計算水位h，在大斷面上插值得a、b點；去掉高于a、b點的實測點，見圖2；過實測點在水位h和地形線間做輔助線，分別求出各單元對應(yīng)的河寬Bi、過水面積Ai、濕周Si；累加Bi得水位h對應(yīng)的水面寬度，累加Ai得水位h對應(yīng)的過水面積，累加Si得水位h對應(yīng)的濕周；水力半徑R=過水面積A/濕周S。

圖2 水力要素計算示意圖

2 計算程序及使用方法

由Excel完成實測斷面數(shù)據(jù)的錄入工作，然后通過MATLAB程序調(diào)用斷面數(shù)據(jù)完成大斷面水力要素計算并將計算結(jié)果輸出至Excel。

2.1 數(shù)據(jù)錄入

首先，根據(jù)程序要求創(chuàng)建名稱為data的Excel文件。在data.xls內(nèi)建立2個工作表，分別用來存儲斷面數(shù)據(jù)(input)和計算結(jié)果(output)，工作表數(shù)據(jù)錄入格式如圖3。

圖3 斷面數(shù)據(jù)錄入格式圖

2.2 計算程序

大斷面水力要素計算MATLAB程序如下：

[ipd,T]=xlsread(′data.xls′,′input′);

[hang,lie]=size(ipd);lie1=lie/2;opd=[];

for sc=1:lie1;

sc0=sc*2;dtA=ipd(:,sc0-1:sc0);

data=dtA(～any(isnan(dtA),2),:);

qdj=data(:,1);gc=data(:,2);

migc=min(gc);masw=floor(max(gc));

misw=ceil(min(gc));dt=0.1;

fdt=-1*dt;misw2=misw-dt;

if misw2>migc

sw0=[migc misw2:fdt:migc misw:dt:masw];

sw=sort(sw0);

elseif misw2<=migc&&misw～=migc

sw=[migc misw:dt:masw];

elseif misw2<=migc&&misw==migc

sw=(misw:dt:masw);

end

dto=[sc 0 0 0 0 sc];n1=length(sw);

for i=1:n1

n2=length(data)-1;h=sw(i);dt0=data;

for j=1:n2

if gc(j)>h&&gc(j+1)

a1=[qdj(j) qdj(j+1)];b1=[gc(j) gc(j+1)];

Jl=interp1(b1,a1,sw(i),′linear′);

dti=[Jl h];dt0=cat(1,dt0,dti);

elseif gc(j)h

a1=[qdj(j) qdj(j+1)];b1=[gc(j) gc(j+1)];

Jl=interp1(b1,a1,h,′linear′);

dti=[Jl h];dt0=cat(1,dt0,dti);

end

end

dt2=sortrows(dt0,1);qdj2=dt2(:,1);

gc2=dt2(:,2);n3=length(dt2)-2;

for L=1:n3

if qdj2(L)==qdj2(L+1)&&qdj2(L)==qdj2(L+2)&&gc2(L)

gc3=[gc2(L) gc2(L+1) gc2(L+2)];gc2(L)=min(gc3);

gc2(L+1)=h;gc2(L+2)=max(gc3);

end

if qdj2(L)==qdj2(L+1)&&qdj2(L)==qdj2(L+2)&&gc2(L)>gc2(L+2)

gc3=[gc2(L) gc2(L+1) gc2(L+2)];gc2(L)=max(gc3);

gc2(L+1)=h;gc2(L+2)=min(gc3);

end

end

a=qdj2;b=gc2;n4=length(a)-1;A=0;S=0;B=0;H=h-migc;

for k=1:n4

A1=(a(k+1)-a(k))*(h-migc);

A2=(a(k+1)-a(k))*(b(k)+b(k+1)-2*migc)/2;

A0=A1-A2;

S0=sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2);

B0=a(k+1)-a(k);c0=[b(k) b(k+1)];

c1=max(c0);c2=min(c0);

if c1>h

A0=0;S0=0;B0=0;

elseif c1==c2&&c1==h

A0=0;S0=0;B0=0;

end

A=A+A0;S=S+S0;B=B+B0;

end

R=A/S;dt4=[h A S R B H];

dto=cat(1,dto,dt4);

end

opd=[opd;dto];

end

xlswrite(′data.xls′,opd,′output′,′A2′);

程序通過MATLAB中xlsread函數(shù)自動調(diào)用data.xls中的斷面數(shù)據(jù)，然后完成水力要素計算，并通過xlswrite函數(shù)將演算結(jié)果輸出至“output”工作表，整個計算過程均由計算機完成且可以同時完成多斷面多水位的水力要素計算，即提高了計算速度又能確保計算精度。

4 應(yīng)用實例

以市橋水道大斷面水力要素計算為例，驗證該方法的方便、準(zhǔn)確性。市橋水道北自廣州市番禺區(qū)鐘村鎮(zhèn)的石壁始，上游為屏山河，南流經(jīng)西海咀、韋涌，下南山峽，西自龍灣河入口，至南山峽與屏山河相匯稱市橋水道，東流經(jīng)市橋城區(qū)和鐘村、沙灣、石碁3鎮(zhèn)，在觀音沙尾匯入沙灣水道后出獅子洋，全長38.60 km。

計算斷面取市橋水道大刀沙村委處橫斷面(見圖4)，計算間隔dt定為1 m。

圖4 市橋水道大刀沙村委處橫斷面圖

根據(jù)市橋水道大刀沙村委處橫斷面實測數(shù)據(jù)，在錄入基礎(chǔ)數(shù)據(jù)后運行MATLAB程序完成計算。程序輸出結(jié)果包括水位(m)、斷面過水面積(m2)、濕周(m)、水力半徑(m)、水面寬(m)和水深(m)，計算結(jié)果見表1。

表1 水力要素計算結(jié)果表

5 結(jié)語

以市橋水道大斷面水力要素計算為例，聯(lián)合運用MATLAB程序和Excel，方便、快速地完成了市橋水道大刀沙村委處大斷面水力要素計算工作。該方法不僅操作方便、計算迅速，而且避免了傳統(tǒng)手工算法的誤差，計算精度也得以提高?？梢?，聯(lián)合運用MATLAB程序和Excel解決大斷面水力要素計算問題，為廣大水利工程技術(shù)人員提供了一個方便、快捷地解決大斷面水力要素計算工作的方法。