《流函数涡量法的二维方腔流数值模拟(matlab编程).doc》由会员分享，可在线阅读，更多相关《流函数涡量法的二维方腔流数值模拟(matlab编程).doc（13页珍藏版）》请在三一办公上搜索。
1、流函数- 涡量法的二维方腔流数值模拟基本方程： 在直角坐标系下，不可压非定常流体所满足的流函数涡量形式的N-S方程为其中为雷诺数差分格式：采用FTCS格式有：对于本问题，将方腔四边同时分为等分，则有故 在直角坐标系下，不可压定常流体所满足的流函数涡量形式的N-S方程为其中为雷诺数差分格式：采用FTCS格式有：对于本问题，将方腔四边同时分为等分，则有，则有即边界条件：在腔体的两侧和顶边，（第二式由泰勒级数展开得到）在底边（第二式由泰勒级数展开得到）其中代表边界，代表与边界相邻的节点。而即Matlab程序为： 不可压非定常流体clear;%参数设置Re=10; %雷诺数取10,100，500，10
2、00L=1; %空穴几何尺寸n=100;dh=L/n;%delta hdt=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)/dh2; xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1)/dh2; end for j=2:n xi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh2; xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j)/dh2; en
3、d %控制方程 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)*dh2)/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(
4、i-1,j)+xi(i,j+1)+xi(i,j-1)-4*xi(i,j)/Re)/dh2; xi(i,j)=xi(i,j)+rho*err2; temp=max(abs(err1),abs(err2); if errtemp err=temp; end end end if (mod(k,1000)=0) %每千步显示结果 k err contour(psi,100);%contour求迹线 pause(0.5) end if err1e-6 break;endend kerrrhodtcontour(psi,100);时，k=9216，err=9.9957e-07，rho=1，dt=1.00
5、00e-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，1000L=1; %空穴几何尺寸n=100;dh=L/n;%delta hpsi=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*(
6、psi(i,2)-psi(i,1)/dh2; xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1)/dh2; end for j=2:n xi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh2; xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j)/dh2; 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,
7、j-1)+xi(i,j)*dh2)/4-psi(i,j); psi(i,j)=psi(i,j)+rho*err1; err2=(xi(i+1,j)+xi(i-1,j)+xi(i,j+1)+xi(i,j-1)/4 . -Re*dh*(u(i,j)*(xi(i+1,j)-xi(i-1,j)+v(i,j)*(xi(i,j+1)-xi(i,j-1)/8-xi(i,j); xi(i,j)=xi(i,j)+rho*err2; temp=max(abs(err1),abs(err2); if errtemp err=temp; end end end if err1e-6 break;endend kerr
8、rho%psicontour(psi,100);时，k=6445，err=9.9978e-07，rho=1.0;时，k =7533，err=9.9953e-07，rho=1.0;时,k =10707，err =9.9973e-07，rho=1.0；时,在不调节松弛因子时，其发散了，通过减小其松弛因子得到k =27600，err=9.9970e-07，rho=0.6000；最后使用时间相关法再次求解该问题，此时在直角坐标系下，不可压非定常流体所满足的流函数涡量形式的N-S方程为其中为雷诺数差分格式：采用FTCS格式有：对于本问题，将方腔四边同时分为等分，则有故Matlab 程序为：clear;%
9、参数设置Re=1000; %雷诺数取100，500，1000L=1; %空穴几何尺寸n=100;dh=L/n;%delta hdt=5e-5; %时间步长psi=zeros(n+1,n+1);xi=zeros(n+1,n+1);rho=1.0; for k=1:1000000 err=0; for i=2:n xi(i,1)=-2*(psi(i,2)-psi(i,1)/dh2; xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1)/dh2; end for j=2:n xi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh2; xi(n+1,j)=-2*(psi
10、(n,j)-psi(n+1,j)/dh2; 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=dt*(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)-4*psi(i,j)/dh2+xi(i,j); psi(i,j)=psi(i,j)+rho*err1; err2=dt/dh2*(-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); xi(i,j)=xi(i,j)+rho*err2; temp=max(abs(err1),abs(err2); if errtemp err=temp; end end end if (mod(k,1000)=0) k err contour(psi,100); pause(0.5) end if err1e-6 break;endend kerrcontour(psi,100);对于该方法，在此处难以收敛，可能需继续调试，这里就不再讨论了。
