MATLAB程序和Excel在大断面水力要素计算中的联合运用
2021-03-25温斯钧张浩安航永
温斯钧,张浩,安航永
(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) 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”工作表,整个计算过程均由计算机完成且可以同时完成多断面多水位的水力要素计算,即提高了计算速度又能确保计算精度。 以市桥水道大断面水力要素计算为例,验证该方法的方便、准确性。市桥水道北自广州市番禺区钟村镇的石壁始,上游为屏山河,南流经西海咀、韦涌,下南山峡,西自龙湾河入口,至南山峡与屏山河相汇称市桥水道,东流经市桥城区和钟村、沙湾、石碁3镇,在观音沙尾汇入沙湾水道后出狮子洋,全长38.60 km。 计算断面取市桥水道大刀沙村委处横断面(见图4),计算间隔dt定为1 m。 图4 市桥水道大刀沙村委处横断面图 根据市桥水道大刀沙村委处横断面实测数据,在录入基础数据后运行MATLAB程序完成计算。程序输出结果包括水位(m)、断面过水面积(m2)、湿周(m)、水力半径(m)、水面宽(m)和水深(m),计算结果见表1。 表1 水力要素计算结果表 以市桥水道大断面水力要素计算为例,联合运用MATLAB程序和Excel,方便、快速地完成了市桥水道大刀沙村委处大断面水力要素计算工作。该方法不仅操作方便、计算迅速,而且避免了传统手工算法的误差,计算精度也得以提高。可见,联合运用MATLAB程序和Excel解决大断面水力要素计算问题,为广大水利工程技术人员提供了一个方便、快捷地解决大断面水力要素计算工作的方法。4 应用实例
5 结语