《水环境数学模型》(P39-47)隐式差分法求解.doc

河流水质模拟计算思路河流水质模型采用如下方程： （1）中：水质指标浓度；纵向扩散系数；平均流速； 外部的源和漏(如支流的影响等)；时间 ；所考察的距离。以水环境数学模型（P39-47）隐式差分法求解如下：由于为水质指标浓度的函数，不妨设，忽略支流影响。则方程（1）可化成下面形式： （2） 对方程（2）整理得： （3）令： 方程（3）可简化为： （4）边界条件： i=1时 i=n时 矩阵形式： 河流水质模拟计算流程符号说明：Dx空间步长Dt时间步长Cjj时刻各断面的污染物浓度Ej时刻各断面的纵向扩散系数Uj时刻各断面平均流速K污染物的衰减系数NT模拟的次数其中：A,a,b,c,d,x为计算辅助矩阵。计算流程图见图一j=0输入Dx,Dt,U,E,KCj输入j+1时刻C(1)i=1a(i)=-E(i+1)/(Dx2)b(i)=1/Dx+2*E(i+1)/(Dx2)+K/2C(i)=a(i)d(i)=Cj(i+1)*(1/Dt-U(i+1)/Dx)+Cj(i)*(U(i+1)/Dx-K/2)i=length(U)-1i=i+1否是i=2A（i,i）=b(i)A(i-1,i)=a(i)A(i,i+1)=c(i)i=length(U)-2否是i=i+1i=1d(1)=b(1)-a(1)*Cj(1)i=n= length(U)-1A(n,n-1)=a(i)-c(i)A(n,n)=b(i)+2*c(i)i=1m=A(i+1,i)/A(i,i)A(i+1,i:i+1)=A(i+1,i:i+1)-m*A(i,i:i+1);b(i+1)=b(i+1)-m*b(i);i=n=rank(A)-1否i=i+1x(n)=b(n)/A(n,n)i=n-1 x(i)=(b(i)-A(i,i+1)*x(i+ 1)/A(i,i); i=1否是是Cj（2:n）=xjNTj=j+1否结束是 图一河流水质隐式差分法体系计算框图算例采用水环境数学模型（P3947）例题：已知：整个河流看成顺直均匀流，E=2KM2/h，U=5km/h，K=0.151/h，Dt=0.1h,Dx=0.5km在x=0处有一个排放时间为1小时的点源。此程序的设计有一个主程序（shuizhi），两个子程序（HLSZAD和fzhuigan）组成，其中主函数（shuizhi）用于控制模拟次数，初始量的输入以及输出j+1时刻个断面的浓度。其调用子函数（HLSZAD）。子函数（HLSZAD）用于整理三对角方程，为求解做准备，其调用子函数（fzhuigan）。子函数（fzhuigan）是为了求解三对角方程，其返回值为j+1时刻第二个计算断面及其后各断面的浓度。利用上述方法计算结果见表1表一 书中计算结果与自己编程结果对比表ji0246810121416010000000002（理论）106.8781.7060.30750.00730.0010.00010.00.0 实际106.87831.70640.30750.04890.00730.00100.00010.04（理论）109.086.2642.6660.79890.19380.04110.00790.0014实际109.07946.26392.6640.79880.19380.04110.00790.00136（理论）109.6428.4556.0153.1571.2510.4010.110.0269实际109.64198.45446.01503.15721.25110.40100.10970.02458（理论）109.8359.3018.0215.8663.4421.6160.6270.208实际109.83469.30128.02105.86573.44151.61550.62460.193710（理论）09.9119.6529.07.7035.763.6251.9030.8451实际09.91029.65158.9957.70295.75933.62391.89660.8047注：（1）“理论”为书中计算结果，“实际”自己编程结果。 （2）程序中书中边界断面为0断面，而编程中边界断面为第1端面总体上是可以的，只是0.2h末时第8,10,12,14,16断面浓度相差较大，一小时末第16断面浓度相差较大。估计0.2h末时第8,10,12,14,16断面浓度为书中记录错误。一小时末第16断面偏差大现在还不清楚。开始编程时将1小时作为一个计算时间段，即j=10时的瞬间，第0断面其浓度突然由10mg/l突变成0mg/l，这次的突变对后面各断面是没有影响的，而在初次编程时考虑了突变的影响（不合理），导致j=10时前几个断面与书中所给浓度偏差较大，在10%左右。书中的例子是以BOD5来计算的，其对CODMN和氨氮也同样成立。因为CODMN和氨氮的变化均为一级反应，因此将两者可以化为统一的数学表达式：其中：KCODMN或氨氮的综合衰减系数 P常数，其为其他物质对CODMN或氨氮浓度的影响。附录：MATLAB程序主函数：（shuizhi）function shuizhiE(18)=0;E(:)=2;U(18)=0;U(:)=5;K=0.0151;Cj(18)=0;Cj(1)=10;Dt=0.1;Dx=0.5;j=0;while (j=10 Cj(1)=0; end n=length(Cj); fprintf(n第%d次各断面浓度n,j) for i=1:n fprintf(t%f,Cj(i); end end子函数（HLSZAD）function Cj=HLSZAD(E,U,Dt,Dx,K,Cj)%E- j时刻各断面的纵向扩散系数%U-j时刻各断面的平均流速%Dt-时间步长%Dx-空间步长%K-综合衰减系数%Cj-j时刻各断面浓度n=length(U)-1;a(n)=0;b(n)=0;c(n)=0;d(n)=0;for i=1:n a(i)=-E(i+1)/Dx2; b(i)=1/Dt+2*E(i+1)/Dx2+K/2; c(i)=-E(i+1)/Dx2; d(i)=Cj(i+1)*(1/Dt-U(i+1)/Dx)+Cj(i)*(U(i+1)/Dx-K/2);endA=zeros(n,n);for i=1:n A(i,i)=b(i); if(i1) A(i,i-1)=a(i); end if(in) A(i,i+1)=c(i); endendA(n,n-1)=a(n)-c(n);A(n,n)=b(n)+2*c(n);d(1)=d(1)-a(1)*Cj(1);CK=fzhuigan(A,d);%解三对角方程Cj(2:n+1)=CK;子函数（fzhuigan）function x=fzhuigan(A,b)n=rank(A);for i=1:n-1 m=A(i+1,i)/A(i,i)