# my new matlab code

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 August 17, 2011, 00:16 my new matlab code #1 Member   HouKen Join Date: Jul 2011 Posts: 66 Rep Power: 6 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 arshiya and Yoking Wang like this.

 August 26, 2011, 00:19 please explain the problem u have solved #2 Member   Prasanth P Join Date: May 2009 Posts: 30 Rep Power: 8 Please specify the problem that you have solved and the algorithm that you have used for the same.

 September 7, 2015, 07:32 Bugs #3 New Member   mgedwards Join Date: Aug 2015 Posts: 2 Rep Power: 0 Best explain the problem You have left some bugs

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post M Main CFD Forum 2 April 11, 2010 04:32 Mink Main CFD Forum 3 May 23, 2007 06:03 John C. Chien Main CFD Forum 19 May 17, 2001 15:56 John C. Chien Main CFD Forum 54 April 23, 2001 08:10 Heinz Wilkening Main CFD Forum 38 March 5, 1999 12:44

All times are GMT -4. The time now is 08:14.

 Contact Us - CFD Online - Top