|
[Sponsors] |
August 17, 2011, 01:16 |
my new matlab code
|
#1 |
Member
HouKen
Join Date: Jul 2011
Posts: 67
Rep Power: 15 |
hey guys
i have made a new matlab code on laminar flow.It uses a staggered grid and SIP for iteration method.If u r also a matlab user working on cfd,u can download it and just run it on your pc. And also i hope u can make some advice. Thank u %%%%%%%%%%%%%%%%%2D-cartesian coordinate SIMPLE Method%%%%%%%%%%%%%%%%%%%%% clc clear all; format long; %initialization errdata=zeros(1,100000); nx=120; ny=30; %velcity field u=zeros(ny+2,nx+1,3); v=zeros(ny+1,nx+2,3); ua=zeros(ny+2,nx+1); va=zeros(ny+1,nx+2); ua(16:17,40:42)=10^25; va(15:17,41:42)=10^25; %pressure p=zeros(ny,nx,3); %co. for momentum equation aeu=zeros(ny+2,nx+1); awu=zeros(ny+2,nx+1); asu=zeros(ny+2,nx+1); anu=zeros(ny+2,nx+1); apu=zeros(ny+2,nx+1); aev=zeros(ny+1,nx+2); awv=zeros(ny+1,nx+2); asv=zeros(ny+1,nx+2); anv=zeros(ny+1,nx+2); apv=zeros(ny+1,nx+2); %co. for SIP method ls=zeros(ny+2,nx+1); lw=zeros(ny+2,nx+1); un=zeros(ny+2,nx+1); ue=zeros(ny+2,nx+1); lpr=zeros(ny+2,nx+1); res=zeros(ny+2,nx+1); resl=0; %co. for pressure correction equation aepc=zeros(ny,nx); awpc=zeros(ny,nx); anpc=zeros(ny,nx); aspc=zeros(ny,nx); appc=zeros(ny,nx); pc=zeros(ny+2,nx+2,3); %mass unbalance mp=zeros(ny,nx); %dimensionless density,grid size,viscosity raw=1; dx=1/ny; dy=1/ny; d=0.1; alpha=0.0002; %the inlet velocity=RE for i=2:ny+1 u(i,:,=-6*((0.5/ny+(i-2)/ny)-0.5)^2+1.5; end u(16:17,40:42)=0; v(15:17,41:42)=0; for j=1:nx p(:,j,=(1.2*nx/ny)/nx*(1-j); end %at this time mp is all zero for i=1:ny for j=1:nx mp(i,j)=raw*dx*((u(i,j+1,1)-u(i,j,1))+(v(i+1,j,1)-v(i,j,1))); end end %main loop for q=1:100000 %%%%%%%%%%%%%%%%%%%%%%%%%momentum equation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ll=1:1 u(16:17,40:42)=0; v(15:17,41:42)=0; for i=2:ny+1 for j=2:nx if i==2 aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))); awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))); asu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))); anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))); %awu(:,2)=0; %aeu(:,nx)=0; apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%(dy*(1.2*nx/ny)/nx+0.1*(u(i+1,j,1)-u(i,j,1)))/u(i,j,1); %u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy+ua(i,j)*u(i,j,1))/(apu(i,j)+ua(i,j)); else if i==ny+1 aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))); awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))); asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))); anu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))); %awu(:,2)=0; %aeu(:,nx)=0; apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%+(dy*(1.2*nx/ny)/nx+0.1*(u(i-1,j,1)-u(i,j,1)))/u(i,j,1); %u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy+ua(i,j)*u(i,j,1))/(apu(i,j)+ua(i,j)); else aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))); awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))); asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))); anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))); %awu(:,2)=0; %aeu(:,nx)=0; apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j); %u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy+ua(i,j)*u(i,j,1))/(apu(i,j)+ua(i,j)); end end end end %SIP METHOD for u velocity su=zeros(ny+2,nx+1); urf=0.02; %co. for SIP method ls=zeros(ny+2,nx+1); lw=zeros(ny+2,nx+1); un=zeros(ny+2,nx+1); ue=zeros(ny+2,nx+1); lpr=zeros(ny+2,nx+1); res=zeros(ny+2,nx+1); resl=0; for i=2:ny+1 for j=2:nx apu(i,j)=-(apu(i,j)+ua(i,j))/0.8; su(i,j)=0.2*apu(i,j)*u(i,j,1)+dy*(p(i-1,j,1)-p(i-1,j-1,1))-ua(i,j)*u(i,j,1); end end for j=2:nx for i=2:ny+1 lw(i,j)=awu(i,j)/(1+0.9*un(i,j-1)); ls(i,j)=asu(i,j)/(1+0.9*ue(i-1,j)); p1=0.9*lw(i,j)*un(i,j-1); p2=0.9*ls(i,j)*ue(i-1,j); lpr(i,j)=1/(apu(i,j)+p1+p2-lw(i,j)*ue(i,j-1)-ls(i,j)*un(i-1,j)+10^-25); un(i,j)=(anu(i,j)-p1)*lpr(i,j); ue(i,j)=(aeu(i,j)-p2)*lpr(i,j); end end for l=1:3 resl=0; for i=2:ny+1 for j=2:nx res(i,j)=su(i,j)-(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+apu(i,j)*u(i,j,1)); resl=resl+abs(res(i,j)); res(i,j)=(res(i,j)-ls(i,j)*res(i-1,j)-lw(i,j)*res(i,j-1))*lpr(i,j); end end if l==1 resor=resl; end % rsm=resl/(resor+10^-25) resl; for i=(ny+1):-1:2 for j=nx:-1:2 res(i,j)=res(i,j)-un(i,j)*res(i+1,j)-ue(i,j)*res(i,j+1); u(i,j,1)=u(i,j,1)+res(i,j); end end end %u(:,nx,1)=u(:,nx,1)*(sum(u(:,1,1))/sum(u(:,nx,1))); u(:,nx+1,1)=u(:,nx,1);%for a outlet condition for i=2:ny for j=2:nx+1 if j==2 aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))); awv(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j-1,1)+0.25*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))); asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))); anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); else if j==nx+1 aev(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j,1)+0.25*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))); awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))); asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))); anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); else aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))); awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))); asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))); anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); end end v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy+va(i,j)*v(i,j,1))/(apv(i,j)+va(i,j)); end end %SIP METHOD for v velocity sv=zeros(ny+1,nx+2); urf=0.02; %co. for SIP method ls=zeros(ny+1,nx+2); lw=zeros(ny+1,nx+2); un=zeros(ny+1,nx+2); ue=zeros(ny+1,nx+2); lpr=zeros(ny+1,nx+2); res=zeros(ny+1,nx+2); resl=0; for i=2:ny for j=2:nx+1 apv(i,j)=-(apv(i,j)+va(i,j))/0.8; sv(i,j)=0.2*apv(i,j)*v(i,j,1)+dx*(p(i,j-1,1)-p(i-1,j-1,1))-va(i,j)*v(i,j,1); end end for j=2:nx+1 for i=2:ny lw(i,j)=awv(i,j)/(1+0.9*un(i,j-1)); ls(i,j)=asv(i,j)/(1+0.9*ue(i-1,j)); p1=0.9*lw(i,j)*un(i,j-1); p2=0.9*ls(i,j)*ue(i-1,j); lpr(i,j)=1/(apv(i,j)+p1+p2-lw(i,j)*ue(i,j-1)-ls(i,j)*un(i-1,j)+10^-25); un(i,j)=(anv(i,j)-p1)*lpr(i,j); ue(i,j)=(aev(i,j)-p2)*lpr(i,j); end end for l=1:3 resl=0; for i=2:ny for j=2:nx+1 res(i,j)=sv(i,j)-(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+apv(i,j)*v(i,j,1)); resl=resl+abs(res(i,j)); res(i,j)=(res(i,j)-ls(i,j)*res(i-1,j)-lw(i,j)*res(i,j-1))*lpr(i,j); end end if l==1 resor=resl; end % rsm=resl/(resor+10^-25) resl; for i=(ny):-1:2 for j=(nx+1):-1:2 res(i,j)=res(i,j)-un(i,j)*res(i+1,j)-ue(i,j)*res(i,j+1); v(i,j,1)=v(i,j,1)+res(i,j); end end end for i=1:ny for j=1:nx %mp(i,j)=raw*dx*((u(i+1,j+1,1)-u(i+1,j,1))+(v(i+1,j+1,1)-v(i,j+1,1))); end end %norm(mp);%mp become nonzero! end %%%%%%%%%%%%%%%%%%%pressure correction%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %co. for pressure correction for i=2:ny+1 for j=2:nx if i==2 aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))); awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))); asu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))); anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))); %awu(:,2)=0; %aeu(:,nx)=0; apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%(dy*(1.2*nx/ny)/nx+0.1*(u(i+1,j,1)-u(i,j,1)))/u(i,j,1); else if i==ny+1 aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))); awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))); asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))); anu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))); %awu(:,2)=0; %aeu(:,nx)=0; apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%+(dy*(1.2*nx/ny)/nx+0.1*(u(i-1,j,1)-u(i,j,1)))/u(i,j,1); else aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))); awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))); asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))); anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))); %awu(:,2)=0; %aeu(:,nx)=0; apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j); end end %u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy)/apu(i,j); end end %u(:,nx,1)=u(:,nx,1)*(sum(u(:,1,1))/sum(u(:,nx,1))); %u(:,nx+1,1)=u(:,nx,1); for i=2:ny for j=2:nx+1 if j==2 aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))); awv(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j-1,1)+0.25*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))); asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))); anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); %v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j); else if j==nx+1 aev(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j,1)+0.25*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))); awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))); asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))); anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); %v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j); else aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))); awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))); asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))); anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); %v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j); end end end end for i=1:ny for j=2:nx-1 aepc(i,j)=raw*dy^2/apu(i+1,j+1); awpc(i,j)=raw*dy^2/apu(i+1,j); end end for i=2:ny-1 for j=1:nx anpc(i,j)=raw*dx^2/apv(i+1,j+1); aspc(i,j)=raw*dx^2/apv(i,j+1); end end for i=1:ny awpc(i,1)=0; aepc(i,1)=raw*dy^2/apu(i+1,2); awpc(i,nx)=raw*dy^2/apu(i+1,nx); aepc(i,nx)=0; end for j=1:nx anpc(1,j)=raw*dx^2/apv(2,j+1); aspc(1,j)=0; anpc(ny,j)=0; aspc(ny,j)=raw*dx^2/apv(ny,j+1); end aepc(15:16,39)=0; awpc(15:16,42)=0; anpc(14,40:41)=0; aspc(17,40:41)=0; for i=1:ny for j=1:nx appc(i,j)=aepc(i,j)+awpc(i,j)+anpc(i,j)+aspc(i,j); end end for i=1:ny for j=1:nx mp(i,j)=raw*dx*((u(i+1,j+1,1)-u(i+1,j,1))+(v(i+1,j+1,1)-v(i,j+1,1))); end end %SIP METHOD for pressure correction sp=zeros(ny+2,nx+2); urf=0.02; %co. for SIP method ls=zeros(ny+2,nx+2); lw=zeros(ny+2,nx+2); un=zeros(ny+2,nx+2); ue=zeros(ny+2,nx+2); lpr=zeros(ny+2,nx+2); res=zeros(ny+2,nx+2); resl=0; for i=2:ny+1 for j=2:nx+1 appc(i-1,j-1)=-appc(i-1,j-1)/0.2; sp(i,j)=0.8*appc(i-1,j-1)*pc(i,j,1)+mp(i-1,j-1); end end for j=2:nx+1 for i=2:ny+1 lw(i,j)=awpc(i-1,j-1)/(1+0.9*un(i,j-1)); ls(i,j)=aspc(i-1,j-1)/(1+0.9*ue(i-1,j)); p1=0.9*lw(i,j)*un(i,j-1); p2=0.9*ls(i,j)*ue(i-1,j); lpr(i,j)=1/(appc(i-1,j-1)+p1+p2-lw(i,j)*ue(i,j-1)-ls(i,j)*un(i-1,j)+10^-25); un(i,j)=(anpc(i-1,j-1)-p1)*lpr(i,j); ue(i,j)=(aepc(i-1,j-1)-p2)*lpr(i,j); end end for l=1:5 resl=0; for i=2:ny+1 for j=2:nx+1 res(i,j)=sp(i,j)-(aepc(i-1,j-1)*pc(i,j+1,1)+awpc(i-1,j-1)*pc(i,j-1,1)+anpc(i-1,j-1)*pc(i+1,j,1)+aspc(i-1,j-1)*pc(i-1,j,1)+appc(i-1,j-1)*pc(i,j,1)); resl=resl+abs(res(i,j)); res(i,j)=(res(i,j)-ls(i,j)*res(i-1,j)-lw(i,j)*res(i,j-1))*lpr(i,j); end end if l==1 resor=resl; end % rsm=resl/(resor+10^-25) resl; for i=(ny+1):-1:2 for j=(nx+1):-1:2 res(i,j)=res(i,j)-un(i,j)*res(i+1,j)-ue(i,j)*res(i,j+1); pc(i,j,1)=pc(i,j,1)+res(i,j); end end end %keep the corner pressure tobe 0 ref=pc(2,2,1); for i=2:ny+1 for j=2:nx+1 pc(i,j,1)=pc(i,j,1)-ref; end end %correct the pressure field for ll=1:200 for i=2:ny+1 for j=2:nx+1 % pc(i,j,1)=(pc(i,j+1,1)*aepc(i-1,j-1)+pc(i,j-1,1)*awpc(i-1,j-1)+pc(i+1,j,1)*anpc(i-1,j-1)+pc(i-1,j,1)*aspc(i-1,j-1)-mp(i-1,j-1))/appc(i-1,j-1); end end end for i=1:ny for j=1:nx p(i,j,1)=p(i,j,1)+alpha*pc(i+1,j+1,1); end end %correct velocity field for i=2:ny+1 for j=2:nx u(i,j,1)=u(i,j,1)+dx/apu(i,j)*(pc(i,j,1)-pc(i,j+1,1)); end end for i=2:ny for j=2:nx v(i,j,1)=v(i,j,1)+dy/apv(i,j)*(pc(i,j,1)-pc(i+1,j,1)); end end %check for residual err=norm(mp(:,) errdata(q)=err; u(16:17,40:42)=0; v(15:17,41:42)=0; if err<10^-8 break end clf subplot(2,2,1),contourf(u(:,:,1),30);axis equal subplot(2,2,2), plot(errdata(1,1:q)); subplot(2,2,3),contourf(pc(:,:,1),30);axis equal subplot(2,2,4),contourf(p(:,:,1),30);axis equal hold off; drawnow end |
|
August 26, 2011, 01:19 |
please explain the problem u have solved
|
#2 |
Member
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 17 |
Please specify the problem that you have solved and the algorithm that you have used for the same.
|
|
September 7, 2015, 08:32 |
Bugs
|
#3 |
New Member
mgedwards
Join Date: Aug 2015
Posts: 2
Rep Power: 0 |
Best explain the problem
You have left some bugs |
|
July 1, 2016, 17:28 |
|
#4 |
New Member
Kalyan
Join Date: Jun 2016
Posts: 5
Rep Power: 10 |
Hey,
Can you explain how have you set up the grid? I'm having problems implementing SIMPLE in MATLAB for flow between plates. Thanks in advance, Kalyan |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Matlab code for lid-driven cavity | M | Main CFD Forum | 2 | April 11, 2010 05:32 |
matlab code galerkin method? | Mink | Main CFD Forum | 3 | May 23, 2007 07:03 |
Design Integration with CFD? | John C. Chien | Main CFD Forum | 19 | May 17, 2001 16:56 |
What is the Better Way to Do CFD? | John C. Chien | Main CFD Forum | 54 | April 23, 2001 09:10 |
public CFD Code development | Heinz Wilkening | Main CFD Forum | 38 | March 5, 1999 12:44 |