《水环境数学模型》(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）可化成下面形式：

)

(2

1211

11

2

1

1

1111

j

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(

)2

2

1(

11112

1

12

1

12

K 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 ?-

=γ

2

2112

K 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

11

1

1 （4）

边界条件： i=1时 '11

1

11

1

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

11'

1x

E ?+

=δδ

i=n 时 1

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

n j n

n j n n C C δβα=+++-1

'

1

1'

n n n γαα-='

n n n γββ2'

+=

矩阵形式：

?

?????

????????

?---'

'

1112

2

211

n 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 时刻各断面浓度

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)/Dx^2;

b(i)=1/Dt+2*E(i+1)/Dx^2+K/2;

c(i)=-E(i+1)/Dx^2;

d(i)=Cj(i+1)*(1/Dt-U(i+1)/Dx)+Cj(i)*(U(i+1)/Dx-K/2); end

A=zeros(n,n);

for i=1:n

A(i,i)=b(i);

if(i>1)

A(i,i-1)=a(i);

end

if(i

A(i,i+1)=c(i);

end

end

A(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);

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);

end

x=zeros(1,n);

x(n)=b(n)/A(n,n);

for i=n-1:-1:1

x(i)=(b(i)-A(i,i+1)*x(i+1))/A(i,i);

end