lilongduzhi 发表于 2008-10-13 23:25

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

% 洪水模型(圣维南方程的差分求解方法);
% 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 = ;
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 = ;
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 发表于 2012-4-3 20:04

回复 1 # lilongduzhi 的帖子

非常感谢,留下研究~~

西门牛 发表于 2012-4-26 19:11

回复 2 # aimeehoney 的帖子

谢谢分享{:{23}:}

linhao0987 发表于 2012-5-8 08:41

求洪水淹没分析中流向分析程序。{:{11}:}
页: [1]
查看完整版本: Matlab 一维洪水模型(圣维南方程的离散库朗差分解法)程序