Solving Navier Stokes equation over an obstacle

Hi,

I am trying to solve Navier Stokes equation over an obstacle: http://shrani.si/f/3l/13P/2Tihb3iM/projekt2.png

So I have no-slip wall up and on obstacle, and symmetry down.

This is my current code:

nx=64;ny=32;

dt=0.0001;

nstep=10;

mu=0.05;

maxit=1000;

beta=1.2;

h=1/ny;

u1=zeros(nx+1,ny+2);

v1=zeros(nx+2,ny+1);

p1=zeros(nx+2,ny+2);

ut1=zeros(nx+1,ny+2);

vt1=zeros(nx+2,ny+1);

c1=zeros(nx+2,ny+2)+0.25;

uu1=zeros(nx+1,ny+1);

vv1=zeros(nx+1,ny+1);

w1=zeros(nx+1,ny+1);

c1(nx+1,ny/2+1:ny)=1/3;

c1(1:nx+2,2)=1/3;

c1(1:nx,ny+1)=1/3;

c1(nx+1,ny+1)=1/2;

c1(2, 3:ny) = 1/3;

c1(2, 2) = 1/2;

c1(2, ny+1) = 1/2;

u2=zeros(nx+1,ny+2);

v2=zeros(nx+2,ny+1);

p2=zeros(nx+2,ny+2);

ut2=zeros(nx+1,ny+2);

vt2=zeros(nx+2,ny+1);

c2=zeros(nx+2,ny+2)+0.25;

uu2=zeros(nx+1,ny+1);

vv2=zeros(nx+1,ny+1);

w2=zeros(nx+1,ny+1);

c2(1:nx+2,2)=1/3;

c2(1:nx+2,ny/2+1)=1/3;

%c2(nx+1, ny/2+1) = 1/2;

%c2(nx+1, ny+1) = 1/2;

%c2(2,ny+1)=1/2;

%c2(2,ny/2+1)=1/2;

nxx = 128;

h3 = 1/ny;

u3=zeros(nxx+1,ny+2);

v3=zeros(nxx+2,ny+1);

p3=zeros(nxx+2,ny+2);

ut3=zeros(nxx+1,ny+2);

vt3=zeros(nxx+2,ny+1);

c3=zeros(nxx+2,ny+2)+0.25;

uu3=zeros(nxx+1,ny+1);

vv3=zeros(nxx+1,ny+1);

w3=zeros(nxx+1,ny+1);

c3(1,ny/2+1:ny)=1/3;

c3(1:nxx+2,2)=1/3;

c3(2:nxx+2,ny+1)=1/3;

c3(1,ny+1)=1/2;

c3(nxx+1, 2:ny+1) = 1/3;

c3(nxx+1, ny+1) = 1/2;

c3(nxx+1, 2) = 1/2;

time=0.0;

for is=1:nstep

u1(2:nx+1,1)=-u1(2:nx+1,2); %No slip wall

v1(1:nx+2,1)=0; %Normal velocity

u1(1, 2:ny) = 1; %Inlet

u1(2:nx, ny+2) = u1(2:nx, ny+1); %Symmetry

v1(2:nx+1, ny+1) = 0; %Symmetry

v1(nx+2, ny/2+2:ny+1) = -v1(nx+1, ny/2+2:ny+1); %No slip wall

u1(nx+1, ny/2+2:ny+2) = 0; %Normal velocity

%p1(nx+1, ny/2+2:ny+2)=0;

v1(nx+2, 2:ny/2+1) = v2(1, 2:ny/2+1);

v2(1, ny/2+2:ny+1) = v1(nx+2, ny/2+2:ny+1);

u2(1:nx+1,1)=-u2(1:nx+1,2); %No slip wall

v2(1:nx+2,1)= 0; %Normal velocity

u2(1:nx+1, ny/2+3:ny+2) = 0; % No slip wall

v2(1:nx+2, ny/2+2:ny+2) = 0; %No slip wall

p2(1:nx+2, ny/2+2:ny+2) = 0; %No slip wall

u2(1:nx+1, ny/2+2) = -u2(1:nx+1, ny/2+1); %No slip wall

v2(1:nx+2, ny/2+1) = 0; %Normal velocity

v2(nx+1, ny/2+2:ny+1) = -v3(1, ny/2+2:ny+1);

v2(nx+2, 2:ny/2+2) = v3(1, 2:ny/2+2);

u3(1:nxx+1,1)=-u3(1:nxx+1,2); %no slip

v3(1:nxx+2,1)= 0;%

u3(1:nxx+1, ny+2) = u3(1:nxx+1, ny+1); %

v3(1:nxx+2, ny+1) = 0; % wall

u3(1, ny/2+2: ny+1) = 0;% wall

u3(nxx+1, 2:ny+1) = 1; % outlet

for i=2:nx,for j=2:ny+1 %

ut1(i,j)=u1(i,j)+dt*(-(0.25/h)*((u1(i+1,j)+u1(i,j))^2-(u1(i,j)+...

u1(i-1,j))^2+(u1(i,j+1)+u1(i,j))*(v1(i+1,j)+...

v1(i,j))-(u1(i,j)+u1(i,j-1))*(v1(i+1,j-1)+v1(i,j-1)))+...

(mu/h^2)*(u1(i+1,j)+u1(i-1,j)+u1(i,j+1)+u1(i,j-1)-4*u1(i,j)));

end,end

for j=2:ny/2+1 %

ut1(nx+1,j)=u1(nx+1, j)+dt*(-(0.25/h)*((u2(1,j)+u1(nx+1,j))^2-(u1(nx+1,j)+...

u1(nx,j))^2+(u1(nx+1,j+1)+u1(nx+1,j))*(v2(1,j)+...

v1(nx+1,j))-(u1(nx+1,j)+u1(nx+1,j-1))*(v2(1,j-1)+v1(nx+1,j-1)))+...

(mu/h^2)*(u2(1,j)+u1(nx,j)+u1(nx+1,j+1)+u1(nx+1,j-1)-4*u1(nx+1,j)));

end

ut1(nx+1, ny/2+2:ny+1) = 0;

for i=2:nx,for j=2:ny %

vt1(i,j)=v1(i,j)+dt*(-(0.25/h)*((u1(i,j+1)+u1(i,j))*(v1(i+1,j)+...

v1(i,j))-(u1(i-1,j+1)+u1(i-1,j))*(v1(i,j)+v1(i-1,j))+...

(v1(i,j+1)+v1(i,j))^2-(v1(i,j)+v1(i,j-1))^2)+...

(mu/h^2)*(v1(i+1,j)+v1(i-1,j)+v1(i,j+1)+v1(i,j-1)-4*v1(i,j)));

end,end

for j = 2:ny/2+1

vt1(nx+1,j)=v1(nx+1,j)+dt*(-(0.25/h)*((u1(nx+1,j+1)+u1(nx+1,j))*(v2(1,j)+...

v1(nx+1,j))-(u1(nx,j+1)+u1(nx,j))*(v1(nx+1,j)+v1(nx,j))+...

(v1(nx+1,j+1)+v1(nx+1,j))^2-(v1(nx+1,j)+v1(nx+1,j-1))^2)+...

(mu/h^2)*(v2(1,j)+v1(nx,j)+v1(nx+1,j+1)+v1(nx+1,j-1)-4*v1(nx+1,j)));

end

i = nx+1;

for j = ny/2+2:ny

vt1(i,j)=v1(i,j)+dt*(-(0.25/h)*((u1(i,j+1)+u1(i,j))*(v1(i+1,j)+...

v1(i,j))-(u1(i-1,j+1)+u1(i-1,j))*(v1(i,j)+v1(i-1,j))+...

(v1(i,j+1)+v1(i,j))^2-(v1(i,j)+v1(i,j-1))^2)+...

(mu/h^2)*(v1(i+1,j)+v1(i-1,j)+v1(i,j+1)+v1(i,j-1)-4*v1(i,j)));

end

for j=2:ny/2+1

ut2(1,j)=u1(1, j)+dt*(-(0.25/h)*((u2(2,j)+u2(1,j))^2-(u2(1,j)+...

u1(nx+1,j))^2+(u2(1,j+1)+u2(1,j))*(v2(2,j)+...

v2(1,j))-(u2(1,j)+u2(1,j-1))*(v2(2,j-1)+v2(1,j-1)))+...

(mu/h^2)*(u2(2,j)+u1(nx+1,j)+u2(1,j+1)+u2(1,j-1)-4*u2(1,j)));

end

for j = 2:ny/2

vt2(1,j)=v2(1,j)+dt*(-(0.25/h)*((u2(1,j+1)+u2(1,j))*(v2(2,j)+...

v2(1,j))-(u1(nx+1,j+1)+u1(nx+1,j))*(v2(1,j)+v1(nx+1,j))+...

(v2(1,j+1)+v2(1,j))^2-(v2(1,j)+v2(1,j-1))^2)+...

(mu/h^2)*(v2(2,j)+v1(nx+1,j)+v2(1,j+1)+v2(1,j-1)-4*v2(1,j)));

end

for i=2:nx,for j=2:ny/2+1 %

ut2(i,j)=u2(i,j)+dt*(-(0.25/h)*((u2(i+1,j)+u2(i,j))^2-(u2(i,j)+...

u2(i-1,j))^2+(u2(i,j+1)+u2(i,j))*(v2(i+1,j)+...

v2(i,j))-(u2(i,j)+u2(i,j-1))*(v2(i+1,j-1)+v2(i,j-1)))+...

(mu/h^2)*(u2(i+1,j)+u2(i-1,j)+u2(i,j+1)+u1(i,j-1)-4*u2(i,j)));

end,end

for j=2:ny/2 %

ut2(nx+1,j)=u2(nx+1, j)+dt*(-(0.25/h)*((u3(1,j)+u2(nx+1,j))^2-(u2(nx+1,j)+...

u2(nx,j))^2+(u2(nx+1,j+1)+u2(nx+1,j))*(v3(1,j)+...

v2(nx+1,j))-(u2(nx+1,j)+u2(nx+1,j-1))*(v3(1,j-1)+v2(nx+1,j-1)))+...

(mu/h^2)*(u3(2,j)+u2(nx,j)+u2(nx+1,j+1)+u2(nx+1,j-1)-4*u2(nx+1,j)));

end

for i=2:nx,for j=2:ny/2 %

vt2(i,j)=v2(i,j)+dt*(-(0.25/h)*((u2(i,j+1)+u2(i,j))*(v2(i+1,j)+...

v2(i,j))-(u2(i-1,j+1)+u2(i-1,j))*(v2(i,j)+v2(i-1,j))+...

(v2(i,j+1)+v2(i,j))^2-(v2(i,j)+v2(i,j-1))^2)+...

(mu/h^2)*(v2(i+1,j)+v2(i-1,j)+v2(i,j+1)+v2(i,j-1)-4*v2(i,j)));

end,end

for j = 2:ny/2 %

vt2(nx+1,j)=v2(nx+1,j)+dt*(-(0.25/h)*((u2(nx+1,j+1)+u2(nx+1,j))*(v3(1,j)+...

v2(nx+1,j))-(u2(nx,j+1)+u2(nx,j))*(v2(nx+1,j)+v2(nx,j))+...

(v2(nx+1,j+1)+v2(nx+1,j))^2-(v2(nx+1,j)+v2(nx+1,j-1))^2)+...

(mu/h^2)*(v3(1,j)+v2(nx,j)+v2(nx+1,j+1)+v2(nx+1,j-1)-4*v2(nx+1,j)));

end

vt1(nx+2, 2:ny/2+1) = vt2(1, 2:ny/2+1);

%Še tretja regija.

for j=2:ny/2+1

ut3(1,j)=u3(1, j)+dt*(-(0.25/h)*((u3(2,j)+u3(1,j))^2-(u3(1,j)+...

u2(nx+1,j))^2+(u3(1,j+1)+u3(1,j))*(v3(2,j)+...

v3(1,j))-(u3(1,j)+u3(1,j-1))*(v3(2,j-1)+v3(1,j-1)))+...

(mu/h^2)*(u3(2,j)+u2(nx+1,j)+u3(1,j+1)+u3(1,j-1)-4*u3(1,j)));

end

for j = 2:ny %

vt3(1,j)=v3(1,j)+dt*(-(0.25/h)*((u3(1,j+1)+u3(1,j))*(v3(2,j)+...

v3(1,j))-(u2(nx+1,j+1)+u2(nx+1,j))*(v3(1,j)+v2(nx+2,j))+...

(v3(1,j+1)+v3(1,j))^2-(v3(1,j)+v3(1,j-1))^2)+...

(mu/h^2)*(v3(2,j)+v2(nx+2,j)+v3(1,j+1)+v3(1,j-1)-4*v3(1,j)));

end

for i=2:nxx,for j=2:ny %

ut3(i,j)=u3(i,j)+dt*(-(0.25/h)*((u3(i+1,j)+u3(i,j))^2-(u3(i,j)+...

u3(i-1,j))^2+(u3(i,j+1)+u3(i,j))*(v3(i+1,j)+...

v3(i,j))-(u3(i,j)+u3(i,j-1))*(v3(i+1,j-1)+v3(i,j-1)))+...

(mu/h^2)*(u3(i+1,j)+u3(i-1,j)+u3(i,j+1)+u3(i,j-1)-4*u3(i,j)));

end,end

for i=2:nxx+1,for j=2:ny %

vt3(i,j)=v3(i,j)+dt*(-(0.25/h)*((u3(i,j+1)+u3(i,j))*(v3(i+1,j)+...

v3(i,j))-(u3(i-1,j+1)+u3(i-1,j))*(v3(i,j)+v3(i-1,j))+...

(v3(i,j+1)+v3(i,j))^2-(v3(i,j)+v3(i,j-1))^2)+...

(mu/h^2)*(v3(i+1,j)+v3(i-1,j)+v3(i,j+1)+v3(i,j-1)-4*v3(i,j)));

end,end

vt2(nx+2, 2:ny/2+1) = vt3(1, 2:ny/2+1);

for it=1:maxit % solve for pressure

for i=2:nx,for j=2:ny+1

p1(i,j)=beta*c1(i,j)*(p1(i+1,j)+p1(i-1,j)+p1(i,j+1)+p1(i,j-1)-...

(h/dt)*(ut1(i,j)-ut1(i-1,j)+vt1(i,j)-vt1(i,j-1)))+(1-beta)*p1(i,j);

end,end

for j =2:ny+1

p1(nx+1,j)=beta*c1(nx+1,j)*(p2(1,j)+p1(nx,j)+p1(nx+1,j+1)+p1(nx+1,j-1)-...

(h/dt)*(ut1(nx+1,j)-ut1(nx,j)+vt1(nx+1,j)-vt1(nx+1,j-1)))+(1-beta)*p1(nx+1,j);

end

for j =2:ny/2+1

p2(1,j)=beta*c2(1,j)*(p2(2,j)+p1(nx+1,j)+p2(1,j+1)+p2(1,j-1)-...

(h/dt)*(ut2(1,j)-ut1(nx+1,j)+vt2(1,j)-vt2(1,j-1)))+(1-beta)*p2(1,j);

end

for i=2:nx,for j=2:ny/2+1

p2(i,j)=beta*c2(i,j)*(p2(i+1,j)+p2(i-1,j)+p2(i,j+1)+p2(i,j-1)-...

(h/dt)*(ut2(i,j)-ut2(i-1,j)+vt2(i,j)-vt2(i,j-1)))+(1-beta)*p2(i,j);

end,end

p2(1:nx+2, ny/2+2:ny+1) = 0;

for j=2:ny/2+1

p2(nx+1,j)=beta*c2(nx+1,j)*(p3(1,j)+p2(nx,j)+p2(nx+1,j+1)+p2(nx,j-1)-...

(h/dt)*(ut2(nx,j)-ut2(nx,j)+vt2(nx+1,j)-vt2(nx+1,j-1)))+(1-beta)*p2(nx+1,j);

end

for j =2:ny+1

p3(1,j)=beta*c3(1,j)*(p3(2,j)+p2(nx+1,j)+p3(1,j+1)+p3(1,j-1)-...

(h/dt)*(ut3(1,j)-ut2(nx+1,j)+vt3(1,j)-vt3(1,j-1)))+(1-beta)*p3(1,j);

end

for i=2:nxx+1,for j=2:ny+1

p3(i,j)=beta*c3(i,j)*(p3(i+1,j)+p3(i-1,j)+p3(i,j+1)+p3(i,j-1)-...

(h/dt)*(ut3(i,j)-ut3(i-1,j)+vt3(i,j)-vt3(i,j-1)))+(1-beta)*p3(i,j);

end,end

end

% correct the velocity

u1(2:nx,2:ny+1)=...

ut1(2:nx,2:ny+1)-(dt/h)*(p1(3:nx+1,2:ny+1)-p1(2:nx,2:ny+1));

v1(2:nx+1,2:ny)=...

vt1(2:nx+1,2:ny)-(dt/h)*(p1(2:nx+1,3:ny+1)-p1(2:nx+1,2:ny));

u1(nx+1,2:ny/2+1)=...

ut1(nx+1,2:ny/2+1)-(dt/h)*(p2(1,2:ny/2+1)-p1(nx+1,2:ny/2+1));

u2(1,2:ny/2+1)=...

ut2(1,2:ny/2+1)-(dt/h)*(p2(1,2:ny/2+1)-p1(nx+1,2:ny/2+1));

v2(1,2:ny/2)=...

vt2(1,2:ny/2)-(dt/h)*(p2(1,3:ny/2+1)-p2(1,2:ny/2));

u2(2:nx,2:ny/2+1)=...

ut2(2:nx,2:ny/2+1)-(dt/h)*(p2(3:nx+1,2:ny/2+1)-p2(2:nx,2:ny/2+1));

v2(2:nx+1,2:ny/2)=...

vt2(2:nx+1,2:ny/2)-(dt/h)*(p2(2:nx+1,3:ny/2+1)-p2(2:nx+1,2:ny/2));

u2(1:nx+1, ny/2+3:ny+1) = 0;

v2(1:nx+1, ny/2+2:ny+1) = 0;

u2(nx+1,2:ny/2+1)=...

ut2(nx+1,2:ny/2+1)-(dt/h)*(p2(nx+1,2:ny/2+1)-p2(nx,2:ny/2+1));

u3(1,2:ny/2+1)=...

ut3(1,2:ny/2+1)-(dt/h)*(p3(1,2:ny/2+1)-p2(nx+1,2:ny/2+1));

v3(1,2:ny)=...

vt3(1,2:ny)-(dt/h)*(p3(1,3:ny+1)-p3(1,2:ny));

u3(2:nxx,2:ny+1)=...

ut3(2:nxx,2:ny+1)-(dt/h)*(p3(3:nxx+1,2:ny+1)-p3(2:nxx,2:ny+1));

v3(2:nxx+1,2:ny)=...

vt3(2:nxx+1,2:ny)-(dt/h)*(p3(2:nxx+1,3:ny+1)-p3(2:nxx+1,2:ny));

time=time+dt; % plot the results

uu1(1:nx+1,1:ny+1)=0.5*(u1(1:nx+1,2:ny+2)+u1(1:nx+1,1:ny+1));

vv1(1:nx+1,1:ny+1)=0.5*(v1(2:nx+2,1:ny+1)+v1(1:nx+1,1:ny+1));

w1(1:nx+1,1:ny+1)=(u1(1:nx+1,2:ny+2)-u1(1:nx+1,1:ny+1)-...

v1(2:nx+2,1:ny+1)+v1(1:nx+1,1:ny+1))/(2*h);

hold off,quiver(flipud(rot90(uu1)),flipud(rot90(vv1)),'r');

hold on;contour(flipud(rot90(w1)),20),axis equal,pause(0.01)

end

I believe the problem is at pressure part or maybe at dividing it into three regions, I should have specified boundaries differently. If anyone has any recommendations or maybe has a code that is solving N-S in multiple regions, please let me know.

Thanks for all the help.