流函数- 涡量法的二维方腔流数值模拟

基本方程：

在直角坐标系下，不可压非定常流体所满足的流函数涡量形式的N-S方程为

其中

为雷诺数

差分格式：

采用FTCS格式有：

对于本问题，将方腔四边同时分为等分，则有

故

在直角坐标系下，不可压定常流体所满足的流函数涡量形式的N-S方程为

其中

为雷诺数

差分格式：

采用FTCS格式有：

对于本问题，将方腔四边同时分为等分，则有，则有即

边界条件：

在腔体的两侧和顶边，

(第二式由泰勒级数展开得到)

在底边

(第二式由泰勒级数展开得到)

其中代表边界，代表与边界相邻的节点。

而

即

Matlab程序为：

不可压非定常流体

clear;

%参数设置

Re=10; %雷诺数取10,100，500，1000

L=1; %空穴几何尺寸

n=100;

dh=L/n;%delta h

dt=1e-4; %时间步长

psi=zeros(n+1,n+1);

xi=zeros(n+1,n+1);

rho=1;

for k=1:1000000

err=0;

%边界条件

for i=2:n

xi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;

xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2;

end

for j=2:n

xi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh^2;

xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j))/dh^2;

end

%控制方程

for i=2:n

for j=2:n

u(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dh);

v(i,j)=-((psi(i+1,j)-psi(i-1,j))/(2*dh));

err1=(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+xi(i,j)*dh^2)/4-psi(i,j);

psi(i,j)=psi(i,j)+rho*err1;

err2=dt*(-dh/2*(u(i,j)*(xi(i+1,j)-xi(i-1,j)) ...

+v(i,j)*(xi(i,j+1)-xi(i,j-1))) ...

+(xi(i+1,j)+xi(i-1,j)+xi(i,j+1)+xi(i,j-1)-4*xi(i,j))/Re)/dh^2;

xi(i,j)=xi(i,j)+rho*err2;

temp=max(abs(err1),abs(err2));

if err

err=temp;

end

end

end

if (mod(k,1000)==0) %每千步显示结果

k

err

contour(psi,100);%contour求迹线

pause(0.5)

end

if err<1e-6

break;

end

end

k

err

rho

dt

contour(psi,100);

时，k=9216，err=9.9957e-07，rho=1，dt=1.0000e-04；

时，k=10043，err=9.9973e-07，rho=1，dt=1.0000e-03；

时，k=11275，err=9.9948e-07，rho=1，dt=0.0100;

时，k=16458，err=9.9983e-07，rho=1，dt=0.0100；

不可压定常流体

clear;

%参数设置

Re=10; %雷诺数取100，500，1000

L=1; %空穴几何尺寸

n=100;

dh=L/n;%delta h

psi=zeros(n+1,n+1);

xi=zeros(n+1,n+1);

rho=1.0;

for k=1:100000

err=0;

for i=2:n

xi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;

xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2;

en