CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

my new matlab code

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By houkensjtu

Reply
 
LinkBack Thread Tools Display Modes
Old   August 17, 2011, 00:16
Default my new matlab code
  #1
Member
 
HouKen
Join Date: Jul 2011
Posts: 65
Rep Power: 5
houkensjtu is on a distinguished road
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 likes this.
houkensjtu is offline   Reply With Quote

Old   August 26, 2011, 00:19
Default please explain the problem u have solved
  #2
Member
 
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 8
prasanthnitt is on a distinguished road
Please specify the problem that you have solved and the algorithm that you have used for the same.
prasanthnitt is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Matlab code for lid-driven cavity M Main CFD Forum 2 April 11, 2010 04:32
matlab code galerkin method? Mink Main CFD Forum 3 May 23, 2007 06:03
Design Integration with CFD? John C. Chien Main CFD Forum 19 May 17, 2001 15:56
What is the Better Way to Do CFD? John C. Chien Main CFD Forum 54 April 23, 2001 08:10
public CFD Code development Heinz Wilkening Main CFD Forum 38 March 5, 1999 12:44


All times are GMT -4. The time now is 13:06.