Matlab 一维洪水模型（圣维南方程的离散库朗差分解法）程序

lilongduzhi

% 洪水模型（圣维南方程的差分求解方法）；
% All Right Reserved 城市公共安全应急系统中一维洪水演进及其可视化研究 ;
% Modified by lilongduzhi@yahoo.com.cn hehe:)
% RiverLength = 5089.2
clear all ;
zd = [11.018 ,10.915 ,10.810 ,10.710 ,10.610 , ...
      10.510 ,10.410 ,10.300 ,10.200 ,10.100 ,10.000 ];
Zo = [19.518 ,19.515 ,19.510 ,19.510 ,19.510 , ...
      19.510 ,19.510 ,19.500 ,19.500 ,19.500 ,19.500 ];
Qo = [30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ] ;
b = 5;   % 底宽 ；
m = 3;   % 边坡 ；
I = 0.0002;          % 底坡 ；
n = 0.013 ;          % 河床糙率 ；
g = 9.8  ;
% 初始赋值　；
j = 1 ;
z(j ,:) = Zo ;
Q(j ,:) = Qo ;
% 时间步长 ；
t = 0 ;
Times = 30 ;
numOfs = 11 ;
delta_s = 500;
delta_t = 60 ;
% 测试：
Wm =[];
Np = [];
zpp = [];
zll =[];
while j < Times
     t = t + delta_t;
     % 计算各断面在t时刻的水面宽和过水面积(河道断面为梯形) ；
     % 计算断面面积及水面宽 ；
     h = z(j ,:) - zd ;       % 水深　；
     A = (b + m.*h).*h ;      % 过水断面面积 ；
     p = b +(1+m^2.*h).^0.5;  % 湿周 ；
     B = b + 2*m.*h  ;        % 水面宽 ；
     j = j +1 ;               % 时间数 ；
      for i=1:numOfs
        % 计算M点的值（应当是P点的前一时间段的值） ；
        AM = A(i);   
        BM = B(i);
        QM = Q(j-1 ,i) ;
        zM = z(j-1 ,i) ;
        RM = A(i)/p(i) ;
        vM = QM/AM  ;
        M = [vM ,AM ,BM ,QM ,RM ,n ,g];
  wminus = vM - (g*AM/BM)^0.5;
  wplus  = vM + (g*AM/BM)^0.5;
  N = BM*(QM/AM)^2*I-g*n^2*QM^2/(AM*RM^(4/3));
        % 计算P点的左右节点的值 ；
        if i == 1   % 上游 ；
            zP = 19.518 ;  % 上游边界条件，定值 ；
            QE = Q(j-1 ,i+1); zE = z(j-1 ,i+1) ;
            zR = (delta_t/delta_s)*wminus*(zM-zE)+zM;
            QR = (delta_t/delta_s)*wminus*(QM-QE)+QM ;
            QP = QR + BM*wplus*(zP-zR) + N*delta_t ;
        elseif i == numOfs % 下游 ；
            QD = Q(j-1 ,i-1); zD = z(j-1 ,i-1);
            zL = (delta_t/delta_s)*wplus*(zD-zM)+zM ;
            QL = (delta_t/delta_s)*wplus*(QD-QM)+QM ;
            QP = 30+ t/60*(630-30)/20 ;
            QP = min(QP ,630) ;  
            % 水位计算公式：Bw-(M)*(zP - zL)-QP+QL = N*delta_t ;
            zP = (N*delta_t -QL +QP)/(BM*wminus) + zL ;
        else               % 中游节点
        QD = Q(j-1 ,i-1); zD = z(j-1 ,i-1);
        QE = Q(j-1 ,i+1); zE = z(j-1 ,i+1);
        zL = (delta_t/delta_s)*wplus*(zD-zM)+zM ;
        QL = (delta_t/delta_s)*wplus*(QD-QM)+QM ;
        zR = (delta_t/delta_s)*wminus*(zM-zE)+zM;
        QR = (delta_t/delta_s)*wminus*(QM-QE)+QM;
        %  计算P点的值 ；
        zP = (QR - QL + BM*wminus*zL-BM*wplus*zR)...
           / (BM *wminus - BM *wplus);
        QP = QR + BM*wplus*(zP-zR) + N*delta_t;
        end
        % 记录数据：
        Q(j ,i) = QP ;
        z(j ,i) = zP ;
    end
end
i = 7 ;        % 取第七个切面做观测 ，可任选：
subplot(1,2,1)
plot(z(: ,i)); % 上游流量变化 ；
title(['第',num2str(i),'断面水位变化']);
xlabel('时间');
ylabel('水位');
subplot(1,2,2)
plot(Q(: ,i)); % 上游流量变化 ；
title(['第',num2str(i),'断面流量变化']);
xlabel('时间');
ylabel('流量');
% ------------------------------------------- %
PS： 文献所述的有限，当初编这个程序有很多的疑问，后解决了，相信对大家了解圣维南方程的求解会有帮助
研究洪水模型源于参加的研究生数模，因文献的算例有一定的特殊性（给出上下游水位或流量过程线），没有用到提交的论文上
该程序使用文献的数据以及模型进行计算，断面为梯形（已经够特殊了），原来搜索时发现关于洪水模型以及圣维南方程解法方面的资料较少，
现在发给大家一起探讨一下，希望对大家有所帮助，
（有兴趣的同学参看文献《城市公共安全应急系统中一维洪水演进及其可视化研究》同时还有一篇
《应急平台中一维洪水演进模型研究》可以参考)


[ 本帖最后由 lilongduzhi 于 2008-10-13 23:29 编辑 ]

aimeehoney

回复 1 # lilongduzhi 的帖子

非常感谢，留下研究~~

西门牛

回复 2 # aimeehoney 的帖子

谢谢分享

linhao0987

求洪水淹没分析中流向分析程序。