非线性对流扩散方程不同解法稳定性比较

2013—2014学年第二学期 《Matlab 编程技术》作业
专业班级 石工博13-02
研究方向 油气田开发
姓 名 王壮壮
学 号
B********
结合自己研究方向，运用Matlab编写科学计算及可视化或其它相关程序。

要求：
1）将要解决的问题交代清楚（数学模型、目标等）；
2）编写的程序的关键语句要有注释说明；
3）程序能顺利运行，运行结果和编写的m文件一并提交；
4）独立完成。

非线性对流扩散方程不同解法稳定性比较
流体力学基本方程组本身就是非线性的对流扩散方程，非线性Burgers 方程就是N-S 方程很好的模型方程，它的一维形式如下：
L x x
u x u u t u ≤≤∂∂=∂∂+∂∂022μ (1) 边界条件为
⎩
⎨
⎧====0,,00
u L x u u x (2) 初始条件是任意可以给出的。

我们知道，遇到对流项，我们用迎风格式是绝对没有问题，无论是一阶迎风还是二阶迎风格式都是能够解决非线性对流方程的，如果网格Peclet 数允许的话，中心差分也是可以考虑的。

不过，对于非线性对流，我们来看看另外两个G-S 格式，一个是G-S 型迎风半隐格式，另一个是G-S 型Samarskii 半隐格式，对于任何类型的对流扩散方程，可以收敛到定常解，并且是绝对稳定的，特别适合于解决定常问题。

对于式(1)这两个格式分别为
()
2
11
111111212h u u u R h u u u u u n i n i n i n i n i n i n
i
n
i n i +-+++-+++-+=-+-μτ (3) 21
1
1111112112h u u u R R h u u u u u n i n i n i n i n i n i n i n
i
n
i n i +-+++-+++-⎪⎪⎭⎫ ⎝⎛++=-+-μτ
(4) 其中
μ
2h u R n i n i =
式(3)就是G-S 型迎风半隐格式，它具有一阶精度，是从一阶迎风格式发展而来的；式(4)是G-S 型Samarskii 半隐格式，具有二阶精度，它是从Samarskii 格式发展而来的。

上面说过，它们只适用于求解定常解，因此上标的时层n 可以看作是迭代步，可以说它们没有时间精度，如果想用这两个格式求解非定常解，那可是徒劳的。

对于上两式，我们可以采用迭代法求解，把它们写成迭代式，分别为
()[](
)(
)
(
)
n
i n
i
n i n i n i n i n i n i n i
R h u h u u R u u u h u
++++++--=
+-++-++142212*********τμτμτ (5)
()[]()
⎪⎪⎭
⎫ ⎝⎛+++++⎪
⎪⎭⎫ ⎝⎛+++--=
+-++-++n i n
i n i n i n i n i n i n i n i n i n i R R h u h u u R R u u u h u 11422112221111111
τμτμτ (6) 这样我们可以看出来，其实时间步长τ就充当了一个松弛因子的角色。

首先式(1)是有定常精确解的，解为
⎥⎥⎥⎥⎦⎤⎢⎢
⎢
⎢⎣
⎡⎪⎭⎫ ⎝⎛-+⎪⎭⎫ ⎝⎛
--=L L x u L L x u u u u L L Re exp 1Re exp 10 (7) 式中
μ
L
u L 0Re =
(8)
()
L u u u Re ex p 1
1
-=+- (9)
程序3.3 非线性Burgers 方程计算程序
%-------------------------------------程序开始---------------------------------
>> clear
>> Re=100; %式(8) L=5; %计算域，可以随意设置 u0=1; %内边界条件，u0可随意 miu=u0*L/Re; %扩散系数 N=41; %网格节点数
t=2; %时间步长(松弛因子) h=L/(N-1); %网格步长 for i=1:N %网格划分 x(i)=(i-1)*h; end
%---------------------G-S 型迎风半隐格式计算----------------------- u=zeros(N,1); %迭代初值 u(1)=u0; %内边界条件 u1=u; tol=1;p=0;
while tol>=1e-5 %收敛精度，0.00001 for i=2:N-1
R(i)=abs(u(i))*h/2/miu;
u1(i)=(-t*h*(u(i)*(u(i+1)-u1(i-1)))+2*t*miu*(1+R(i))*(u(i+1)+u1(i-1))
+2*h^2*u(i))/(2*h^2+4*t*miu*(1+R(i)));
end
tol=max(abs(u-u1)); %计算两个迭代层的误差
u=u1;
p=p+1;
if p>=5000
disp('G-S型迎风半隐格式不收敛！')
break
end
end
%------------------------G-S型Samarskii型半隐格式计算--------------- tol=1;p=0;
v=zeros(N,1); %迭代初值
v(1)=u0; %内边界条件
v1=v;
while tol>=1e-5 %收敛精度，0.00001
for i=2:N-1
R(i)=abs(v(i))*h/2/miu;
v1(i)=(-t*h*(v(i)*(v(i+1)-v1(i-1)))+2*t*miu*(1/(1+R(i))+R(i))*(v(i+1) +v1(i-1))+2*h^2*v(i))/(2*h^2+4*t*miu*(1/(1+R(i))+R(i)));
end
tol=max(abs(v1-v)); %计算量迭代层误差
v=v1;
p=p+1;
if p>=5000
disp('G-S型Samarskii型半隐格式不收敛！')
break
end
end
%-----------------------牛顿迭代法计算 --------------------
tol=1;U=2; %误差；迭代初值
while tol>=1e-10 %收敛精度0.0000000001
f=(U-1)/(U+1)-exp(-U*Re); %式(9)建立的函数
fd=1/(U+1)-(U-1)/(U+1)^2+Re/exp(U*Re); %式(9)建立的函数的导数
U1=U-f/fd;
tol=abs(U-U1); %计算误差
U=U1;
end
%--------------------------计算精确解-----------------------------------
for i=1:N
ue(i)=u0*U*((1-exp(U*Re*(x(i)/L-1)))/(1+exp(U*Re*(x(i)/L-1)))); end
plot(x,ue,'-*',x,u,'-or',x,v,'-<g') %作图
Re=5
Re=20
Re=100
式(8)定义的参数L Re 就相当于流体力学中的雷诺数，如果雷诺数较大的话，表明对流占优。

上三个图就表示了随L Re 增大的各个解的情况。

我们看到在L Re 比较小的时候，G -S 迎风格式和G -S 型Samarskii 格式都与精确解十分接近，Sanaskii 格式的精度更高一些。

在L Re 较大的情况下，G -S 迎风格式脱离精确解的情况比较严重，可以明显看出它伴随着很强的假扩散效应，而Samarskii 格式却不存在这种现象，或者说这种现象不明显。

所以，我们说这两种差分格式特别适用于各种对流扩散方程，计算也十分稳定。