SIMPLE算法-MATLAB

%1控制体边界x=1;y=1;%密度Dense=1.2;%粘性系数miu=1.85e-2;%2网格（这里可先使用均匀网格）%网格划分要求%主网格划分份数x方向为Nx,y方向为Ny,保证dx=dy,Nx=15;Ny=15;dx=x/(Nx);dy=y/(Ny);%主网格节点分别为X方向和Y方向NX=Nx+1;NY=Ny+1;Npx=NX;Npy=NY;%速度U的网格数Nux=NX+1;Nuy=NY;%速度V的分量Nvx=Nux-1;Nvy=Nuy+1;%3边界值%很重要%计算未知网格归结为计算内部网格，但为写程序方便，编号依然按照常规规定，具体调用ADI时，转化为标准即可%编程原则：边界上的点必须要出现，必须出现在网格标记中。%入口p(1,1:Npy)=1;v(1,1:Nvy)=0;%出口p(Npx,1:Npy)=0;%上边界u(1:Nux,Nuy)=0;fori=2:Npxp(i,Npy)=1-(i-1)*dx;endv(2:Nvx,Nvy)=0;%下边界fori=2:Npxp(i,1)=1-(i-1)*dx;endu(1:Nux,1)=0;v(1:Nvx,1)=0;%4初值,不仅给要计算的点，这主要是为编程方便p(2:Npx-1,2:Npy-1)=0;u(2:Nux,2:Nuy-1)=0;v(2:Nvx,2:Nvy-1)=0.5;p(10,10)=4;pc(1:Npx,1:Npy)=0;sor=0.7;%5核心部分，simple算法的实施fork=1:350u0=u;v0=v;%水平速度离散fori=2:Nux-1forj=2:Nuy-1Feu(i,j)=((Dense*u0(i+1,j))+(Dense*u0(i,j)))/2*dy;Fwu(i,j)=((Dense*u0(i-1,j))+(Dense*u0(i,j)))/2*dy;Fnu(i,j)=((Dense*v0(i-1,j+1))+(Dense*v0(i,j+1)))/2*dx;Fsu(i,j)=((Dense*v0(i-1,j))+(Dense*v0(i,j)))/2*dx;%源项bu(i,j)=(p(i-1,j)-p(i,j))*dy;endendD=miu;fori=2:Nux-1forj=2:Nuy-1Aeu(i,j)=D+max(-Feu(i,j),0);Awu(i,j)=D+max(Fwu(i,j),0);Anu(i,j)=D+max(-Fnu(i,j),0);Asu(i,j)=D+max(Fsu(i,j),0);endendAeu(Nux-1,2:Nuy-1)=0;Awu(2,2:Nuy-1)=0;fori=2:Nux-1forj=2:Nuy-1Apu(i,j)=Aeu(i,j)+Awu(i,j)+Anu(i,j)+Asu(i,j);bu(i,j)=bu(i,j)+(1-sor)*u0(i,j)*Apu(i,j)/sor;endend%6写成ADI可解形%Yakoubibiaozhi=0;fork=1:10000fori=2:Nux-1forj=2:Nuy-1u1(i,j)=(Aeu(i,j)*u(i+1,j)+Awu(i,j)*u(i-1,j)+Anu(i,j)*u(i,j+1)+Asu(i,j)*u(i,j-1)+bu(i,j))/Apu(i,j)*sor;u(i,j)=u1(i,j);endendfori=2:Nx-1forj=2:Ny-1biaozhi=biaozhi+abs(Aeu(i,j)*u(i+1,j)+Awu(i,j)*u(i-1,j)+Anu(i,j)*u(i,j+1)+Asu(i,j)*u(i,j-1)+bu(i,j)-Apu(i,j)/sor*u(i,j));endendbaozhi1=biaozhiifbiaozhi1e-13breakendbiaozhi=0;end%7垂直速度的离散fori=2:Nvx-1forj=2:Nvy-1Fev(i,j)=((Dense*u0(i+1,j-1))+(Dense*u0(i+1,j)))/2*dy;Fwv(i,j)=((Dense*u0(i,j-1))+(Dense*u0(i,j)))/2*dy;Fnv(i,j)=((Dense*v0(i,j))+(Dense*v0(i,j+1)))/2*dx;Fsv(i,j)=((Dense*v0(i,j))+(Dense*v0(i,j-1)))/2*dx;%源项bv(i,j)=(p(i,j-1)-p(i,j))*dx;endendfori=2:Nvx-1forj=2:Nvy-1Aev(i,j)=D+max(-Fev(i,j),0);Awv(i,j)=D+max(Fwv(i,j),0);Anv(i,j)=D+max(-Fnv(i,j),0);Asv(i,j)=D+max(Fsv(i,j),0);endend%局部单项化Aev(Nvx-1,2:Nvy-1)=0;fori=2:Nvx-1forj=2:Nvy-1Apv(i,j)=Aev(i,j)+Awv(i,j)+Anv(i,j)+Asv(i,j);bv(i,j)=(1-sor)*Apv(i,j)/sor*v0(i,j)+bv(i,j);endend%8写成ADI形式biaozhi=0;fork=1:10000fori=2:Nvx-1forj=2:Nvy-1v1(i,j)=(Aev(i,j)*v(i+1,j)+Awv(i,j)*v(i-1,j)+Anv(i,j)*v(i,j+1)+Asv(i,j)*v(i,j-1)+bv(i,j))/Apv(i,j)*sor;v(i,j)=v1(i,j);endendfori=2:Nvx-1forj=2:Nvy-1biaozhi=biaozhi+abs(Aev(i,j)*v(i+1,j)+Awv(i,j)*v(i-1,j)+Anv(i,j)*v(i,j+1)+Asv(i,j)*v(i,j-1)+bv(i,j)-Apv(i,j)/sor*v(i,j));endendbiaozhi2=biaozhiifbiaozhi1e-14breakendbiaozhi=0;end%9压力修正方程%改变速度变量的值%局部单项化forj=2:Nuy-1u(Nux,j)=u(Nux-1,j);u(1,j)=u(2,j);endforj=2:Nvy-1v(Nvx,j)=v(Nvx-1,j);end%10压力修正系数fori=2:Npx-1forj=2:Npy-1Aep(i,j)=Dense*dy*dy/Apu(i+1,j)*sor;Awp(i,j)=Dense*dy*dy/Apu(i,j)*sor;Anp(i,j)=Dense*dx*dx/Apv(i,j+1)*sor;Asp(i,j)=Dense*dx*dx/Apv(i,j)*sor;bp(i,j)=((Dense*u(i,j))-(Dense*u(i+1,j)))*dy+(Dense*v(i,j)-Dense*v(i,j+1))*dx;endend%局部单项化,压力为什么要局部单项化fori=2:Npx-1forj=2:Npy-1App(i,j)=Aep(i,j)+Awp(i,j)+Anp(i,j)+Asp(i,j);endendbiaozhi=0;fork=1:100fori=2:Npx-1forj=2:Npy-1pc1(i,j)=(Aep(i,j)*pc(i+1,j)+Awp(i,j)*pc(i-1,j)+Anp(i,j)*pc(i,j+1)+Asp(i,j)*pc(i,j-1)+bp(i,j))/App(i,j);pc(i,j)=pc1(i,j);endendfori=2:Npx-1forj=2:Npy-1biaozhi=biaozhi+abs(Aep(i,j)*pc(i+1,j)+Awp(i,j)*pc(i-1,j)+Anp(i,j)*pc(i,j+1)+Asp(i,j)*pc(i,j-1)+bp(i,j)-App(i,j)*pc(i,j));endendbiaozhiifbiaozhi1e-14breakendbiaozhi=0;endfori=2:Npx-1forj=2:Npy-1p(i,j)=p(i,j)+0.4*pc(i,j);endend%11修正速度fori=2:Nux-1forj=2:Nuy-1u(i,j)=u(i,j)+sor*dy/Apu(i,j)*(pc(i-1,j)-pc(i,j));endendfori=2:Nvx-1forj=2:Nvy-1v(i,j)=v(i,j)+sor*dx/Apv(i,j)*(pc(i,j-1)-pc(i,j));endendfori=2:Npx-1forj=2:Npy-1bp(i,j)=((Dense*u(i,j))-(Dense*u(i+1,j)))*dy+(Dense*v(i,j)-Dense*v(i,j+1))*dx;endend%局部单项化forj=2:Nuy-1u(Nux,j)=u(Nux-1,j);u(1,j)=u(2,j)endforj=2:Nvy-1v(Nvx,j)=v(Nvx-1,j);endend