流函數- 渦量法的二維方腔流數值模擬

基本方程：

在直角坐標系下，不可壓非定常流體所滿足的流函數渦量形式的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

總結

以上是生活随笔為你收集整理的matlab计算涡度的函数_流函数涡量法的二维方腔流数值模拟matlab编程.doc的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。