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

my new matlab code

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

Like Tree2Likes
  • 2 Post By houkensjtu

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 17, 2011, 01:16
Default my new matlab code
  #1
Member
 
HouKen
Join Date: Jul 2011
Posts: 67
Rep Power: 15
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 and Yoking Wang like this.
houkensjtu is offline   Reply With Quote

Old   August 26, 2011, 01:19
Default please explain the problem u have solved
  #2
Member
 
Prasanth P
Join Date: May 2009
Posts: 30
Rep Power: 17
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

Old   September 7, 2015, 08:32
Default Bugs
  #3
New Member
 
mgedwards
Join Date: Aug 2015
Posts: 2
Rep Power: 0
mgedwards is on a distinguished road
Best explain the problem
You have left some bugs
mgedwards is offline   Reply With Quote

Old   July 1, 2016, 17:28
Default
  #4
New Member
 
Kalyan
Join Date: Jun 2016
Posts: 5
Rep Power: 10
kalanpi4 is on a distinguished road
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
kalanpi4 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
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 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


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