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

2013—2014学年第二学期 《Matlab编程技术》作业

专业班级 石工博13-02

研究方向 油气田开发

姓 名 王壮壮

学 号 B********

评 分 项 目 一、选题的创新性或二、选题的意与专业相结合的程度 义及难度 (总分20分) (总分20分) 三、程序的可读性、逻辑性 (总分20分) 四、程序运行情况五、其它 综合得分 （成果、效率等） （格式等） (总分100分) (总分30分) （10分）

结合自己研究方向，运用Matlab编写科学计算及可视化或其它相关程序。 要求：

1） 将要解决的问题交代清楚（数学模型、目标等）； 2） 编写的程序的关键语句要有注释说明；

3） 程序能顺利运行，运行结果和编写的m文件一并提交； 4） 独立完成。

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

流体力学基本方程组本身就是非线性的对流扩散方程，非线性Burgers方程就是N-S方程很好的模型方程，它的一维形式如下：

uu2uu20xL (1) txx边界条件为

x0,uu0 (2) xL,u0初始条件是任意可以给出的。

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

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

对于式(1)这两个格式分别为

uin1uin1nn11uin1uinuinnui12ui11u1Ri (3) 22hhniuin1uin其中

1nn111uin1uinuinnui12ui11uRi (4) n22hh1RiniRinuinh2

式(3)就是G-S型迎风半隐格式，它具有一阶精度，是从一阶迎风格式发展而来的；式(4)是G-S型Samarskii半隐格式，具有二阶精度，它是从Samarskii格式发展而来的。上面说过，它们只适用于求解定常解，因此上标的时层n可以看作是迭代步，可以说它们没有时间精度，如果想用这两个格式求解非定常解，那可是徒劳的。

对于上两式，我们可以采用迭代法求解，把它们写成迭代式，分别为

un1i112nhuinuin1uin21Rinuin1uin112hui (5) 2h241Rin

uin111nnn12nhuinuin1uin2Ruu2hui1ii1i11Rni (6) 1n2h24Ri1Rni这样我们可以看出来，其实时间步长就充当了一个松弛因子的角色。 首先式(1)是有定常精确解的，解为

xL1expuReLL (7) uu0uxL1expuReLL式中

ReLu0L (8)

u1expuReL (9) u1

程序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,'-Re=5

Re=20

Re=100

式(8)定义的参数ReL就相当于流体力学中的雷诺数，如果雷诺数较大的话，表明对流占优。上三个图就表示了随ReL增大的各个解的情况。我们看到在ReL比较小的时候，G-S迎风格式和G-S型Samarskii格式都与精确解十分接近，Sanaskii格式的精度更高一些。在ReL较大的情况下，G-S迎风格式脱离精确解的情况比较严重，可以明显看出它伴随着很强的假扩散效应，而Samarskii格式却不存在这种现象，或者说这种现象不明显。

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