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

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