SIMPLE method with matlab

houkensjtu · July 18, 2011, 10:55

Hi,I wrote a Simple method by matlab.
However it's very POOR program and not smart.
I hope it can help someone who need a matlab program.
And I want to know,Why,When I solved out the pressure correction and correct the velocity,the mass is still highly unbalanced(still,the final result seems not bad!)

clear all;
%initialization
nx=100;
ny=30;
%velcity field
u=zeros(ny+2,nx+1,3);
v=zeros(ny+1,nx+2,3);
%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 pressure correction equation
pce=zeros(ny,nx);
pcw=zeros(ny,nx);
pcs=zeros(ny,nx);
pcn=zeros(ny,nx);
aep=zeros(ny,nx);
awp=zeros(ny,nx);
asp=zeros(ny,nx);
anp=zeros(ny,nx);
app=zeros(ny,nx);
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);
errp=zeros(ny,nx);
raw=1000;
dx=0.001/ny;
dy=0.001/ny;
d=0.001;
u(2:ny+1,:,

=0.01;
for q=1:5000
%%%%%%%%%%%%%%%%%%%%%%%%%momentum equation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ll=1:10
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)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i-1,j+1,1)+0.25*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)));
if j==nx
aeu(i,j)=0;
end
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)/apu(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)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i,j+1,1)+0.25*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
if j==nx
aeu(i,j)=0;
end
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)/apu(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)));
if j==nx
aeu(i,j)=0;
end
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)/apu(i,j);
end
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 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)=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)));
if j==nx+1
aev(i,j)=0;
end
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)=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)));
if j==nx+1
aev(i,j)=0;
end
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)));
if j==nx+1
aev(i,j)=0;
end
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

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)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i-1,j+1,1)+0.25*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)));
if j==nx
aeu(i,j)=0;
end
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(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)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i,j+1,1)+0.25*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
if j==nx
aeu(i,j)=0;
end
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(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)));
if j==nx
aeu(i,j)=0;
end
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
end
end
end
end
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)=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)=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
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
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

for ll=1:1000
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
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
for i=1:ny
for j=1:nx
p(i,j,1)=p(i,j,1)+0.01*pc(i+1,j+1,1);
end
end

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
err=norm(mp(:,

)
if err<10^-5
break
end
for i=2:ny-1
for j=2:nx-1
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
end

contourf(u(2:31,:,1),30);axis equal

kid · July 22, 2011, 08:44

hi
Can you explain SIMPLE ?

VincentD · July 23, 2011, 19:51

If you want to improve your code you should make more use of matrices. Matrix calculation is what matlab does best (MATrix LABoratory).

For instance let's say we want to take a second derivative of our velocity field U (2D). In formula form:

(U(i+1,j)-2*U(i,j)+U(i-1,j))/dx^2

We could also write this down in matrix form, namely:

[ 2 -1 0 ; -1 2 -1 ; 0 -1 2 ] or to create a nxn matrix of this type we can write in matlab

A = diag(ones(n-1,1),1)
E = 2*eye(n) - A -A'

Where E is now our desired second difference matrix.

Multiply U*E will give our second derivate at each location (but beware of correctly imposing BCs).

You should check the speed increase, it's huge (Matlabs for loops are very slow). And whats even better it is easier to use and gives more compact coding.

A next step could be to use the spare matrix identity. If you let matlab know if a matrix is spare (like a large nxn difference matrix) it will result in quicker computation and less memory use.

As for your question:
A well posed pressure correction should make your solution divergence free. If this is not the case then something must be off. I would advice to write down (on a piece of paper) the mathematical equations and see for yourself how you want to pressure-correction scheme to work. Using a staggered grid (defining velocity on cell faces, pressure in the center) is the way to go for a proper pressure-velocity coupling.

Good luck!

houkensjtu · July 23, 2011, 21:57

Quote:

Originally Posted by VincentD

If you want to improve your code you should make more use of matrices. Matrix calculation is what matlab does best (MATrix LABoratory).

For instance let's say we want to take a second derivative of our velocity field U (2D). In formula form:

(U(i+1,j)-2*U(i,j)+U(i-1,j))/dx^2

We could also write this down in matrix form, namely:

[ 2 -1 0 ; -1 2 -1 ; 0 -1 2 ] or to create a nxn matrix of this type we can write in matlab

A = diag(ones(n-1,1),1)
E = 2*eye(n) - A -A'

Where E is now our desired second difference matrix.

Multiply U*E will give our second derivate at each location (but beware of correctly imposing BCs).

You should check the speed increase, it's huge (Matlabs for loops are very slow). And whats even better it is easier to use and gives more compact coding.

A next step could be to use the spare matrix identity. If you let matlab know if a matrix is spare (like a large nxn difference matrix) it will result in quicker computation and less memory use.

As for your question:
A well posed pressure correction should make your solution divergence free. If this is not the case then something must be off. I would advice to write down (on a piece of paper) the mathematical equations and see for yourself how you want to pressure-correction scheme to work. Using a staggered grid (defining velocity on cell faces, pressure in the center) is the way to go for a proper pressure-velocity coupling.

Good luck!

Thank u very very much for your reply!It really helped me since I'm beginner to matlab.I will take your advice to improve my code.

As for the pressure correction,yes,I studied it again and I notice BECAUSE in SIMPLE algorithm,it usually drop the temporary mass unbalance when doing the velocity correction,thus,even I get a perfect solution for pressure correction equation,it will not totally make the velocity field to be divergence free.Am I right??

vxv · July 23, 2011, 22:55

The whole idea of solving the pressure correction equation is to make the velocity divergence free. After solving the pressure correction equation and correcting the velocities, theoretically speaking continuity equation should be exaclty satisfied i.e. velocity field should be exactly divergence free (enough though momentum equation is not). But usually, the velocity corrections (after pressure correction equation is solved) use approximate formulas. Hence, there is small residual error in the continuity equation. However, usually we iterate over momentum and pressure-correction equations and eventually the residual errors can be reduced to machine zero.

houkensjtu · July 24, 2011, 02:06

Quote:

Originally Posted by vxv

The whole idea of solving the pressure correction equation is to make the velocity divergence free. After solving the pressure correction equation and correcting the velocities, theoretically speaking continuity equation should be exaclty satisfied i.e. velocity field should be exactly divergence free (enough though momentum equation is not). But usually, the velocity corrections (after pressure correction equation is solved) use approximate formulas. Hence, there is small residual error in the continuity equation. However, usually we iterate over momentum and pressure-correction equations and eventually the residual errors can be reduced to machine zero.

yeah,just as u said,in SIMPLE method,both velocity correction and pressure correction is approximated.Therefore,although finally it will go to a convergence,during the loop it's not necessary velocity field will be divergence free after correction.In fact I think it will not be divergence free(after a velocity correction everytime) during the loop because during the loop,the term dropped in velocity correction equation is certainly NOT zero.

VincentD · July 24, 2011, 06:07

I will explain how steps to take when running the SIMPLE method and also what could be the problem.

1. We guess an initial pressure field p*
2. Use the momentum equations to determine the velocity field u*,v*

This field does not correctly represent the flow field since it's not divergence free. Our pressure field is bound to be the wrong guess.

3. We take the following correction the the velocity and pressure:
u** = u* + u'
v** = v* + v'
p** = p* + p'

where ue' = dy/ae*(Pp-Pe) (Pp is pressure in center cell, Pe is pressure in ajacent cell, ae is the cell coefficient).

4. Substituting these corrected velocity into a continuity equation yields us our pressure correction.

5. This process is repeated untill convergence

For stability a underrelaxation coefficient is introduced.

p** = p* + A*p'

This coefficient depends on reynolds number and the pressure correction should be small in respect to the original pressure.

So a couple of questions:

Are you indeed using a staggered grid? - Otherwise a phenomena know as checkerboarding can occur.
Did you try changing the value of your relaxation coefficient?
Are you positive your solution is converged?

Good luck!

dsmith · September 8, 2011, 04:25

Quote:

Originally Posted by VincentD

If you want to improve your code you should make more use of matrices. Matrix calculation is what matlab does best (MATrix LABoratory).

For instance let's say we want to take a second derivative of our velocity field U (2D). In formula form:

(U(i+1,j)-2*U(i,j)+U(i-1,j))/dx^2

We could also write this down in matrix form, namely:

[ 2 -1 0 ; -1 2 -1 ; 0 -1 2 ] or to create a nxn matrix of this type we can write in matlab

A = diag(ones(n-1,1),1)
E = 2*eye(n) - A -A'

Where E is now our desired second difference matrix.

Multiply U*E will give our second derivate at each location (but beware of correctly imposing BCs).

You should check the speed increase, it's huge (Matlabs for loops are very slow). And whats even better it is easier to use and gives more compact coding.

A next step could be to use the spare matrix identity. If you let matlab know if a matrix is spare (like a large nxn difference matrix) it will result in quicker computation and less memory use.

As for your question:
A well posed pressure correction should make your solution divergence free. If this is not the case then something must be off. I would advice to write down (on a piece of paper) the mathematical equations and see for yourself how you want to pressure-correction scheme to work. Using a staggered grid (defining velocity on cell faces, pressure in the center) is the way to go for a proper pressure-velocity coupling.

Good luck!

Vincent,

Can you please explain how you got from your equation for the 2nd derivative of U to the 3x3 matrix? You have me pretty stumped here!

Cheers,

VincentD · September 8, 2011, 06:24

Quote:

Vincent,

Can you please explain how you got from your equation for the 2nd derivative of U to the 3x3 matrix? You have me pretty stumped here!

Cheers,

Sure not a problem. First of all I would like to apologize for a couple of mistakes I made in that post, especially if that was what had you stumped.

(U(i+1,j)-2*U(i,j)+U(i-1,j))/dx^2 is the 2nd derivative

If we want to write this as a matrix multiplication we need to build a matrix so that

A*U = (U(i+1,j)-2*U(i,j)+U(i-1,j))/dx^2

In order to do so we need -2's on the diagonal and 1's on the offdiagonals (I switched the sign before, sorry!).
Using a function like the one below we actaully build this matrix:

function A = K1(n,h,a11,a22)
% a11,a22: Neumann=-1, Dirichlet=-2, Dirichlet mid=-3;
A = spdiags([1 a11 0;ones(n-2,1)*[1 -2 1];0 a22 1],-1:1,n,n)'/h^2;

Here n is equal to the size of our matrix U and h is the delta x. In order to properly enforce the boundary conditions of our system we need to use the correct values at the points (1,1) and (n,n).

It might be insightfull to type:
A = spdiags([1 -1 0;ones(5-2,1)*[1 -2 1];0 -1 1],-1:1,5,5)'
B = full(A)

Then you have a nice picture of what happens in case of -1 values at the side (1,1 and n,n). The 2nd difference for the first cell will be equal to -U(i,j) +U(i+1,j), which is the same as taking U(i,j)=U(i-1,j) so du/dx=0 (Neumann) boundary condition.

Good luck!

majesty · August 15, 2013, 08:28

Hi,

I could not figure out your grid system. Did you use staggered grid or something else? What are your boundary conditions? I am working on 2D simulation of the flow inside a heat pipe and there is something wrong with my code (most probably in the implementation of the BCs) so I am stucked and have been trying to figure out whats wrong for a long while. Moreover why did you re-calculate the same coefficients in pressure correction equation?

Thanks in advance for the response.

July 18, 2011, 10:55	SIMPLE method with matlab	#1
houkensjtu Member HouKen Join Date: Jul 2011 Posts: 67 Rep Power: 15	Hi,I wrote a Simple method by matlab. However it's very POOR program and not smart. I hope it can help someone who need a matlab program. And I want to know,Why,When I solved out the pressure correction and correct the velocity,the mass is still highly unbalanced(still,the final result seems not bad!) clear all; %initialization nx=100; ny=30; %velcity field u=zeros(ny+2,nx+1,3); v=zeros(ny+1,nx+2,3); %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 pressure correction equation pce=zeros(ny,nx); pcw=zeros(ny,nx); pcs=zeros(ny,nx); pcn=zeros(ny,nx); aep=zeros(ny,nx); awp=zeros(ny,nx); asp=zeros(ny,nx); anp=zeros(ny,nx); app=zeros(ny,nx); 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); errp=zeros(ny,nx); raw=1000; dx=0.001/ny; dy=0.001/ny; d=0.001; u(2:ny+1,:,=0.01; for q=1:5000 %%%%%%%%%%%%%%%%%%%%%%%%%momentum equation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ll=1:10 for i=2:ny+1 for j=2:nx if i==2 aeu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))); awu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))); asu(i,j)=2dmax(0,(1-0.1abs(rawdy(0.25v(i-1,j+1,1)+0.25v(i-1,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))); anu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))); if j==nx aeu(i,j)=0; end 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)/apu(i,j); else if i==ny+1 aeu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))); awu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))); asu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))); anu(i,j)=2dmax(0,(1-0.1abs(rawdy(0.25v(i,j+1,1)+0.25v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))); if j==nx aeu(i,j)=0; end 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)/apu(i,j); else aeu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))); awu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))); asu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))); anu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))); if j==nx aeu(i,j)=0; end 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)/apu(i,j); end 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 i=2:ny for j=2:nx+1 if j==2 aev(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))); awv(i,j)=2dmax(0,(1-0.1abs(rawdy(0.25u(i,j-1,1)+0.25u(i+1,j-1,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))); asv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))); anv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))); if j==nx+1 aev(i,j)=0; end 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)=2dmax(0,(1-0.1abs(rawdy(0.25u(i,j,1)+0.25u(i+1,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))); awv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))); asv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))); anv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))); if j==nx+1 aev(i,j)=0; end 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)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))); awv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))); asv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))); anv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))); if j==nx+1 aev(i,j)=0; end 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 end %%%%%%%%%%%%%%%%%%%pressure correction%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %co. for pressure correction for i=2:ny+1 for j=2:nx if i==2 aeu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))); awu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))); asu(i,j)=2dmax(0,(1-0.1abs(rawdy(0.25v(i-1,j+1,1)+0.25v(i-1,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))); anu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))); if j==nx aeu(i,j)=0; end apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j); else if i==ny+1 aeu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))); awu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))); asu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))); anu(i,j)=2dmax(0,(1-0.1abs(rawdy(0.25v(i,j+1,1)+0.25v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))); if j==nx aeu(i,j)=0; end apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j); else aeu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j+1,1)+0.5u(i,j,1))); awu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i,j,1))); asu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j+1,1)+0.5v(i-1,j,1))); anu(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i,j+1,1)+0.5v(i,j,1))); if j==nx aeu(i,j)=0; end apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j); end end end end for i=2:ny for j=2:nx+1 if j==2 aev(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))); awv(i,j)=2dmax(0,(1-0.1abs(rawdy(0.25u(i,j-1,1)+0.25u(i+1,j-1,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))); asv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))); anv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i+1,j,1)+0.5v(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)=2dmax(0,(1-0.1abs(rawdy(0.25u(i,j,1)+0.25u(i+1,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))); awv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))); asv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))); anv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); else aev(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))/d))^5)+max(0,-rawdy(0.5u(i,j,1)+0.5u(i+1,j,1))); awv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))/d))^5)+max(0,rawdy(0.5u(i,j-1,1)+0.5u(i+1,j-1,1))); asv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))/d))^5)+max(0,rawdy(0.5v(i-1,j,1)+0.5v(i,j,1))); anv(i,j)=dmax(0,(1-0.1abs(rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))/d))^5)+max(0,-rawdy(0.5v(i+1,j,1)+0.5v(i,j,1))); apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j); end end end end for i=1:ny for j=2:nx-1 aepc(i,j)=rawdy^2/apu(i+1,j+1); awpc(i,j)=rawdy^2/apu(i+1,j); end end for i=2:ny-1 for j=1:nx anpc(i,j)=rawdx^2/apv(i+1,j+1); aspc(i,j)=rawdx^2/apv(i,j+1); end end for i=1:ny awpc(i,1)=0; aepc(i,1)=rawdy^2/apu(i+1,2); awpc(i,nx)=rawdy^2/apu(i+1,nx); aepc(i,nx)=0; end for j=1:nx anpc(1,j)=rawdx^2/apv(2,j+1); aspc(1,j)=0; anpc(ny,j)=0; aspc(ny,j)=rawdx^2/apv(ny,j+1); end 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)=rawdx((u(i+1,j+1,1)-u(i+1,j,1))+(v(i+1,j+1,1)-v(i,j+1,1))); end end for ll=1:1000 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 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 for i=1:ny for j=1:nx p(i,j,1)=p(i,j,1)+0.01pc(i+1,j+1,1); end end 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 err=norm(mp(:,) if err<10^-5 break end for i=2:ny-1 for j=2:nx-1 mp(i,j)=rawdx((u(i,j+1,1)-u(i,j,1))+(v(i+1,j,1)-v(i,j,1))); end end end contourf(u(2:31,:,1),30);axis equal albalu1366 and Salttz like this.

July 23, 2011, 22:55		#5
vxv New Member vxv Join Date: Jul 2011 Posts: 6 Rep Power: 15	The whole idea of solving the pressure correction equation is to make the velocity divergence free. After solving the pressure correction equation and correcting the velocities, theoretically speaking continuity equation should be exaclty satisfied i.e. velocity field should be exactly divergence free (enough though momentum equation is not). But usually, the velocity corrections (after pressure correction equation is solved) use approximate formulas. Hence, there is small residual error in the continuity equation. However, usually we iterate over momentum and pressure-correction equations and eventually the residual errors can be reduced to machine zero. __________________ vxv

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
How to calculate drag/lift forces, using SIMPLE method.	abcdef123	Main CFD Forum	4	December 7, 2009 14:16
Apply Multigrid method to SIMPLE Algorithm	yang	Main CFD Forum	1	February 25, 2006 12:28
need help in fotran code for simple method	saritha	Main CFD Forum	3	April 30, 2004 13:55
Output Boundary settingon Simple Method	Yoshi	Main CFD Forum	1	March 22, 2004 11:56
SIMPLE method for compressible flow in pipe	Reza	Main CFD Forum	0	August 27, 2002 23:43

July 22, 2011, 08:44		#2
kid Senior Member cfdkid Join Date: Mar 2009 Posts: 133 Rep Power: 17	hi Can you explain SIMPLE ?

July 23, 2011, 19:51		#3
VincentD New Member Vincent Join Date: Jul 2011 Posts: 29 Rep Power: 15	If you want to improve your code you should make more use of matrices. Matrix calculation is what matlab does best (MATrix LABoratory). For instance let's say we want to take a second derivative of our velocity field U (2D). In formula form: (U(i+1,j)-2U(i,j)+U(i-1,j))/dx^2 We could also write this down in matrix form, namely: [ 2 -1 0 ; -1 2 -1 ; 0 -1 2 ] or to create a nxn matrix of this type we can write in matlab A = diag(ones(n-1,1),1) E = 2eye(n) - A -A' Where E is now our desired second difference matrix. Multiply U*E will give our second derivate at each location (but beware of correctly imposing BCs). You should check the speed increase, it's huge (Matlabs for loops are very slow). And whats even better it is easier to use and gives more compact coding. A next step could be to use the spare matrix identity. If you let matlab know if a matrix is spare (like a large nxn difference matrix) it will result in quicker computation and less memory use. As for your question: A well posed pressure correction should make your solution divergence free. If this is not the case then something must be off. I would advice to write down (on a piece of paper) the mathematical equations and see for yourself how you want to pressure-correction scheme to work. Using a staggered grid (defining velocity on cell faces, pressure in the center) is the way to go for a proper pressure-velocity coupling. Good luck!

July 24, 2011, 06:07		#7
VincentD New Member Vincent Join Date: Jul 2011 Posts: 29 Rep Power: 15	I will explain how steps to take when running the SIMPLE method and also what could be the problem. 1. We guess an initial pressure field p* 2. Use the momentum equations to determine the velocity field u,v This field does not correctly represent the flow field since it's not divergence free. Our pressure field is bound to be the wrong guess. 3. We take the following correction the the velocity and pressure: u** = u* + u' v** = v* + v' p** = p* + p' where ue' = dy/ae(Pp-Pe) (Pp is pressure in center cell, Pe is pressure in ajacent cell, ae is the cell coefficient). 4. Substituting these corrected velocity into a continuity equation yields us our pressure correction. 5. This process is repeated untill convergence For stability a underrelaxation coefficient is introduced. p* = p* + A*p' This coefficient depends on reynolds number and the pressure correction should be small in respect to the original pressure. So a couple of questions: Are you indeed using a staggered grid? - Otherwise a phenomena know as checkerboarding can occur. Did you try changing the value of your relaxation coefficient? Are you positive your solution is converged? Good luck!

August 15, 2013, 08:28		#10
majesty New Member Join Date: Jul 2013 Posts: 2 Rep Power: 0	Hi, I could not figure out your grid system. Did you use staggered grid or something else? What are your boundary conditions? I am working on 2D simulation of the flow inside a heat pipe and there is something wrong with my code (most probably in the implementation of the BCs) so I am stucked and have been trying to figure out whats wrong for a long while. Moreover why did you re-calculate the same coefficients in pressure correction equation? Thanks in advance for the response.