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

流函数-涡量法的二维方腔流数值模拟基本方程：①在直角坐标系下，不可压非定常流体所满足的流函数涡量形式的N-S方程为221Reuvtxy其中,uvyxRe为雷诺数差分格式：采用FTCS格式有：1,,1,1,,1,11,,1,,1,,1,,2222122Re()()nnnnnnnnnnnnijijijijijijijijijijijijnnijijuvtxyxy1,,1,,1,,1,2222=0()()ijijijijijijijxy,1,11,1,,,,22nnnnijijijijnijijuvyx对于本问题，将方腔四边同时分为n等分，则有xyh故1,1,,1,1,1,,,1,1,,,1,124+()2Rennnnnijijijijijnnnnnnnnijijijijijijijijthuvh21,1,,1,1,1,,,()4nnnnnijijijijijnnnijijijh,1,11,1,,,,22nnnnijijijijnijijuvhh②在直角坐标系下，不可压定常流体所满足的流函数涡量形式的N-S方程为221Re=0uvxy其中,uvyxRe为雷诺数差分格式：采用FTCS格式有：1,1,,1,11,,1,,1,,1,22221+=22Re()()ijijijijijijijijijijnijijuvxyxy1,,1,,1,,1,2222=0()()ijijijijijijijxy,1,11,1,,,,22ijijijijijijuvyx对于本问题，将方腔四边同时分为n等分，则有xyh，则有即1,1,,1,11,,,1,1,,,1,1,Re=+48nnnnijijijijnnnnnnnnnijijijijijijijijijhuv21,1,,1,1,1,,,()4nnnnnijijijijijnnnijijijh,1,11,1,,,,22nnnnijijijijnijijuvhh边界条件：在腔体的两侧和顶边，*202()()ssssh（第二式由泰勒级数展开得到）在底边*202()()sssshh（第二式由泰勒级数展开得到）其中s代表边界，*s代表与边界相邻的节点。而,,xshys当代表左右两边当代表上下两边xyh即2,1,1,22(+)()jjjhh,2,1,122()()iiih,1,1,22()()IjIjIjh,,1,122()()iJiJiJhMatlab程序为：①不可压非定常流体clear;%参数设置Re=10;%雷诺数取10,100，500，1000L=1;%空穴几何尺寸n=100;dh=L/n;%deltahdt=1e-4;%时间步长psi=zeros(n+1,n+1);xi=zeros(n+1,n+1);rho=1;fork=1:1000000err=0;%边界条件fori=2:nxi(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;endforj=2:nxi(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%控制方程fori=2:nforj=2:nu(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));iferrtemperr=temp;endendendif(mod(k,1000)==0)%每千步显示结果kerrcontour(psi,100);%contour求迹线pause(0.5)endiferr1e-6break;endendkerrrhodtcontour(psi,100);Re=10时，k=9216，err=9.9957e-07，rho=1，dt=1.0000e-04；Re=100时，k=10043，err=9.9973e-07，rho=1，dt=1.0000e-03；Re=500时，k=11275，err=9.9948e-07，rho=1，dt=0.0100;Re=1000时，k=16458，err=9.9983e-07，rho=1，dt=0.0100；②不可压定常流体clear;%参数设置Re=10;%雷诺数取100，500，1000L=1;%空穴几何尺寸n=100;dh=L/n;%deltahpsi=zeros(n+1,n+1);xi=zeros(n+1,n+1);rho=1.0;fork=1:100000err=0;fori=2:nxi(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;endforj=2:nxi(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;endfori=2:nforj=2:nu(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=(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));iferrtemperr=temp;endendendiferr1e-6break;endendkerrrho%psicontour(psi,100);Re=10时，k=6445，err=9.9978e-07，rho=1.0;Re=100时，k=7533，err=9.9953e-07，rho=1.0;Re=500时,k=10707，err=9.9973e-07，rho=1.0；Re=1000时,在不调节松弛因子时，其发散了，通过减小其松弛因子得到k=27600，err=9.9970e-07，rho=0.6000；最后使用时间相关法再次求解该问题，此时在直角坐标系下，不可压非定常流体所满足的流函数涡量形式的N-S方程为221Reuvtxyt其中,uvyxRe为雷诺数差分格式：采用FTCS格式有：1,,1,1,,1,11,,1,,1,,1,,2222122Re()()nnnnnnnnnnnnijijijijijijijijijijijijnnijijuvtxyxy1,,1,,1,,1,,1,2222()()nnnnnnnnijijijijijijijijnijtxy,1,11,1,,,,22nnnnijijijijnijijuvyx对于本问题，将方腔四边同时分为n等分，则有xyh故1,1,,1,1,1,,,1,1,,,1,124+()2Rennnnnijijijijijnnnnnnnnijijijijijijijijthuvh1,1,,1,1,1,,,24()nnnnnijijijijijnnnijijijth,1,11,1,,,,22nnnnijijijijnijijuvhhMatlab程序为：clear;%参数设置Re=1000;%雷诺数取100，500，1000L=1;%空穴几何尺寸n=100;dh=L/n;%deltahdt=5e-5;%时间步长psi=zeros(n+1,n+1);xi=zeros(n+1,n+1);rho=1.0;fork=1:1000000err=0;fori=2:nxi(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;endforj=2:nxi(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;endfori=2:nforj=2:nu(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))/dh^2+xi(i,j));psi(i,j)=psi(i,j)+rho*err1;err2=dt/dh^2*(-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));iferrtemperr=temp;endendendif(mod(k,1000)==0)kerrcontour(psi,100);pause(0.5)endiferr1e-6break;e