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

河流水质模拟计算思路

河流水质模型采用如下方程：

S dt dC

x C E x C u t C L ++∂∂=∂∂+∂∂2 （1） 中：C —水质指标浓度；

L E —纵向扩散系数；

u —平均流速；

S —外部的源和漏(如支流的影响等)；

t —时间 ；

x —所考察的距离。

以《水环境数学模型》（P39-47）隐式差分法求解如下：

由于dt dC 为水质指标浓度C 的函数，不妨设KC dt dC =，忽略支流影响。则方程（1）可化成下面形式：

)(2121111

2

111111j i j i j

i j i i j i j i j i i j i j i C C K x C C u x C C C E t C C -+-+-+++++-∆--∆+-=∆-

（2） 对方程（2）整理得：

)2()1()221(11112

112112K x u C x u t C C x

E C K x E t C x E i j i i j i j i i j i i j i i -∆+∆-∆=∆-+∆+∆+∆-

-+-+++ （3）

令： 2

x

E i

i ∆-

=γ 221

12K x E t i i +∆+∆=

β 2

x E i

i ∆-

=α )2

()1(

11K x u C x u t C i j

i i j i i -∆+∆-∆=-δ 方程（3）可简化为：

i j i i j i i j i i C C C δγβα=++++++-1

1111 （4）

边界条件： i=1时 '

111

111

1δγβ=+++j j C C

21

1'1

x

E ∆+=δδ

i=n 时 1

11112+-+++-=j n j n j n C C C

n j n

n j n n C C δβα=+++-1'11' n n n γαα-=' n n n γββ2'+=

矩阵形式：

⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣

⎡---''11

12

2

2

11n n

n n n βαγβαγβαγβ

⎥⎥⎥⎥⎥⎥⎦⎤

⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥

⎥⎥

⎥⎥⎦

⎤⎢⎢⎢

⎢⎢⎢⎣⎡-+-n n j n n C C C C δδδδ12'11

121

河流水质模拟计算流程符号说明： Dx —空间步长 Dt —时间步长

Cj —j 时刻各断面的污染物浓度 E —j 时刻各断面的纵向扩散系数 U —j 时刻各断面平均流速 K —污染物的衰减系数 NT —模拟的次数

其中：A,a,b,c,d,x 为计算辅助矩阵。 计算流程图见图一

图一河流水质隐式差分法体系计算框图

算例采用《水环境数学模型》（P39—47）例题：

已知：整个河流看成顺直均匀流，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

表一书中计算结果与自己编程结果对比表

注：（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%左右。 书中的例子是以BOD 5来计算的，其对COD MN 和氨氮也同样成立。因为COD MN 和氨氮的变化均为一级反应，因此将两者可以化为统一的数学表达式：

P KC dt

dC

+-= 其中：K —COD MN 或氨氮的综合衰减系数

P —常数，其为其他物质对COD MN 或氨氮浓度的影响。 附录：MATLAB 程序 主函数：（shuizhi ）

function shuizhi E(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=HLSZAD(E,U,Dt,Dx,K,Cj); j=j+1; if 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 时刻各断面浓度