MATLAB程序和Excel在大斷面水力要素計算中的聯合運用

溫斯鈞，張浩，安航永

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

0 引言

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

1 基本原理

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

圖1 實測大斷面示意圖

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

圖2 水力要素計算示意圖

2 計算程序及使用方法

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

2.1 數據錄入

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

圖3 斷面數據錄入格式圖

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函數自動調用data.xls中的斷面數據，然后完成水力要素計算，并通過xlswrite函數將演算結果輸出至“output”工作表，整個計算過程均由計算機完成且可以同時完成多斷面多水位的水力要素計算，即提高了計算速度又能確保計算精度。

4 應用實例

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

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

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

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

表1 水力要素計算結果表

5 結語

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