MATLAB程序和Excel在大斷面水力要素計(jì)算中的聯(lián)合運(yùn)用

溫斯鈞，張浩，安航永

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

0 引言

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

1 基本原理

實(shí)測(cè)大斷面示意圖如圖1所示，計(jì)算水位為h。

圖1 實(shí)測(cè)大斷面示意圖

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

圖2 水力要素計(jì)算示意圖

2 計(jì)算程序及使用方法

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

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

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

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

2.2 計(jì)算程序

大斷面水力要素計(jì)算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′);

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

4 應(yīng)用實(shí)例

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

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

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

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

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

5 結(jié)語(yǔ)

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