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

1、20132014学年第二学期Matlab编程技术作业 专业班级 石工博13-02 研究方向 油气田开发 姓 名 王壮壮 学 号 B13020075 评 分 项 目一、选题的创新性或与专业相结合的程度(总分20分)二、选题的意义及难度(总分20分)三、程序的可读性、逻辑性(总分20分)四、程序运行情况（成果、效率等）(总分30分)五、其它（格式等）（10分）综合得分(总分100分)结合自己研究方向，运用Matlab编写科学计算及可视化或其它相关程序。要求：1） 将要解决的问题交代清楚（数学模型、目标等）；2） 编写的程序的关键语句要有注释说明；3） 程序能顺利运行，运行结果和编写的m文件一并提交

2、；4） 独立完成。非线性对流扩散方程不同解法稳定性比较流体力学基本方程组本身就是非线性的对流扩散方程，非线性Burgers方程就是N-S方程很好的模型方程，它的一维形式如下： (1)边界条件为 (2)初始条件是任意可以给出的。我们知道，遇到对流项，我们用迎风格式是绝对没有问题，无论是一阶迎风还是二阶迎风格式都是能够解决非线性对流方程的，如果网格Peclet数允许的话，中心差分也是可以考虑的。不过，对于非线性对流，我们来看看另外两个G-S格式，一个是G-S型迎风半隐格式，另一个是G-S型Samarskii半隐格式，对于任何类型的对流扩散方程，可以收敛到定常解，并且是绝对稳定的，特别适合于解决定常

3、问题。对于式(1)这两个格式分别为 (3) (4)其中式(3)就是G-S型迎风半隐格式，它具有一阶精度，是从一阶迎风格式发展而来的；式(4)是G-S型Samarskii半隐格式，具有二阶精度，它是从Samarskii格式发展而来的。上面说过，它们只适用于求解定常解，因此上标的时层n可以看作是迭代步，可以说它们没有时间精度，如果想用这两个格式求解非定常解，那可是徒劳的。对于上两式，我们可以采用迭代法求解，把它们写成迭代式，分别为 (5) (6)这样我们可以看出来，其实时间步长就充当了一个松弛因子的角色。首先式(1)是有定常精确解的，解为 (7)式中 (8) (9)程序3.3 非线性Burgers

4、方程计算程序%-程序开始->> 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(

5、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*h2*u(i)/(2*h2+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=

6、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*h2*v(i)/(2*h2+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%-

7、牛顿迭代法计算 - 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=5Re=20Re=100式(8)定义的参数就相当于流体力学中的雷诺数，如果雷诺数较大的话，表明对流占优。上三个图就表示了随增大的各个解的情况。我们看到在比较小的时候，G-S迎风格式和G-S型Samarskii格式都与精确解十分接近，Sanaskii格式的精度