[Jan995]

Hello there,

I am writing my own MATLAB CFD code for a certain problem.



So up I have no-slip wall and down I have a full-slip wall or symmetry. I have divided problem into three parts, left, step, and right.
I dont know if I am making the boundaries between these regions right, because my results are not good and are actually not converging if I apply everything.

What I am actually doing is this:


u2(1, 1:ny/2) = u1(nx, 1:ny/2);
v2(1, 1:ny/2) = v1(nx, 1:ny/2);
u1(nx+1, 1:ny/2) = u2(2, 1:ny/2);
v1(nx+1, 1:ny/2) = v2(2, 1:ny/2);
p2(1, 1:ny/2+1) = p1(nx+1, 1:ny/2+1);
p1(nx+2, 1:ny/2+1) = p2(2, 1:ny/2+1);

And I am not sure if this is the right way. If anyone has any recommendations, I am really going to be grateful for any help, because I am almost out of ideas.
Here is my code:

nx=32;ny=32;
dt=0.00001;
nstep=100;
mu=0.05;
maxit=100;
beta=1.2;
h=1/nx;
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(2:nx+1,2)=1/3;
c1(2:nx,ny+1)=1/3;
c1(nx+1, ny/2+1) = 1/2;
c1(nx+1,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(nx+1,ny/2+1:ny)=1/3;
c2(2,ny/2+1:ny)=1/3;
c2(2:nx+1,2)=1/3;
c2(2:nx,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 = 64;
h3 = 1/nxx;
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(2,ny/2+1:ny)=1/3;
c3(2:nxx+1,2)=1/3;
c3(2:nxx,ny+1)=1/3;
c3(2, ny/2+1) = 1/2;
c3(2,ny+1)=1/2;

un=0;
us=0;
time=0.0;
for is=1:nstep
%No-slip wall up
u1(1:nx+1,1)=2*us-u1(1:nx+1,2);
u2(1:nx+1,1)=2*un-u2(1:nx+1,2);
u3(1:nxx+1,1)=2*un-u3(1:nxx+1,2);
%Left slip wall left, full-slip wall down and no-slip wall right by u1.
u1(1, 1:ny+1) = 1;
u1(1:nx+1, ny+1) = u1(1:nx+1, ny);
v1(1:nx+1, ny+1) = 0;
v1(nx+1, ny/2:ny+1) = -v1(nx, ny/2:ny+1);
%At u2 there is no-slip wall.
u2(1:nx+1, ny/2+2:ny+1) = 0;
v2(1:nx+1, ny/2+2:ny+1) = 0;
u2(1:nx+1, ny/2+1) = -u2(1:nx+1, ny/2);
u2(1, 1:ny/2) = u1(nx, 1:ny/2);
v2(1, 1:ny/2) = v1(nx, 1:ny/2);
u1(nx+1, 1:ny/2) = u2(2, 1:ny/2);
v1(nx+1, 1:ny/2) = v2(2, 1:ny/2);
%p2(1, 1:ny/2+1) = p1(nx+1, 1:ny/2+1);
%p1(nx+2, 1:ny/2+1) = p2(2, 1:ny/2+1);
%Še u3 in v3
%u3(1, 1:ny/2) = u2(nx, 1:ny/2);
%v3(1, 1:ny/2) = v2(nx, 1:ny/2);
%p3(1, 1:ny/2) = p2(nx+1, 1:ny/2);
%p2(nx+1, 1:ny/2) = p3(2, 1:ny/2);
%u2(nx+1, 1:ny/2) = u3(2, 1:ny/2);
%v2(nx+1, 1:ny/2) = v3(2, 1:ny/2);
v3(1, ny/2+1: ny+1) = -v3(2, ny/2+1:ny+1);
u3(nxx+1, 1:ny+1) = 0.5;

for i=2:nx,for j=2:ny+1 % temporary u-velocity
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 i=2:nx+1,for j=2:ny % temporary v-velocity
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 i=2:nx,for j=2:ny/2+1 % temporary u-velocity
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 i=2:nx+1,for j=2:ny/2 % temporary v-velocity
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 i=2:nxx,for j=2:ny+1 % temporary u-velocity
ut3(i,j)=u3(i,j)+dt*(-(0.25/h3)*((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/h3^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 % temporary v-velocity
vt3(i,j)=v3(i,j)+dt*(-(0.25/h3)*((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/h3^2)*(v3(i+1,j)+v3(i-1,j)+v3(i,j+1)+v3(i,j-1)-4*v3(i,j)));
end,end

for it=1:maxit % solve for pressure
for i=2:nx+1,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
end
for it=1:maxit % solve for pressure
for i=2:nx+1,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
end
for it=1:maxit % solve for pressure
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)-...
(h3/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));

u2(2:nx,2:ny+1)=...
ut2(2:nx,2:ny+1)-(dt/h)*(p2(3:nx+1,2:ny+1)-p2(2:nx,2:ny+1));
v2(2:nx+1,2:ny)=...
vt2(2:nx+1,2:ny)-(dt/h)*(p2(2:nx+1,3:ny+1)-p2(2:nx+1,2:ny));
u2(1:nx+1, ny/2+2:ny+1) = 0;
v2(1:nx+1, ny/2+2:ny+1) = 0;

u3(2:nxx,2:ny+1)=...
ut3(2:nxx,2:ny+1)-(dt/h3)*(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/h3)*(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

[FMDenaro]

Symmetry condition is prescribed on the same vector, why are you using index 1 and 2?

[Jan995]

Originally Posted by FMDenaro

Symmetry condition is prescribed on the same vector, why are you using index 1 and 2?

Hi,



I am applying symmetry on the same vector.



I am using index 1 and 2 at the boundary of the regions 1/2 and 2/3 to get proper updated versions of vectors in all three regions. Region 2 needs last vector from region 1 and uses it like an inflow and region 1 needs second vector from region 2 to properly calculate last vector or outflow. Is this way of thinking wrong? How should I apply boundary conditions here then?

[Jan995]

Hello again,

I have modified my code and I think the boundaries are now done correctly. I still have 3 regions, where the middle region is the region with an obstacle.

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.

[FMDenaro]

I don't think that anyone will look into all your code...


I can say what I always done in my code. You define the matrix for the variable that considers a unique region, with the obstacle. Then you have just a mask that marks the node of the solid obstacle.

The cycles are along the indexes of the matrix,


For example a very generic cycle



J=Ny.. Nytop_ob,-1
I = 1, Nx


...


End

End


J=Nytop_ob,1,-1
I = 1, Nxfront_ob
...
End do
Endo do


You do not need to define different arrays for each region. Bcs are prescribed in a straightforward way

[Jan995]

Thank you for your reply, it was really helpful.


I think I have now fairly good code, I just have one more problem left. The inlet and outlet boundary conditions. I just prescribe velocity at the inlet as u(inlet) = 1 and the same for outlet u(outlet) = 1. I do not prescribe any pressure. And if I prescribe it like that, the results I get are not okay. u(inlet+1) is usually fairly low number and if I iterate for long time I even get backflow which then diverges. Although, if I prescribe some p(inlet) and/or p(outlet) I get some fairly good results, but I dont believe this is the way to do it.
Any more guidance in that direction?


Thank you very much!



My code can now be found here:
https://github.com/JankoTheBuck/Comp...ter/NSObstacle

[FMDenaro]

Originally Posted by Jan995

Thank you for your reply, it was really helpful.


I think I have now fairly good code, I just have one more problem left. The inlet and outlet boundary conditions. I just prescribe velocity at the inlet as u(inlet) = 1 and the same for outlet u(outlet) = 1. I do not prescribe any pressure. And if I prescribe it like that, the results I get are not okay. u(inlet+1) is usually fairly low number and if I iterate for long time I even get backflow which then diverges. Although, if I prescribe some p(inlet) and/or p(outlet) I get some fairly good results, but I dont believe this is the way to do it.
Any more guidance in that direction?


Thank you very much!



My code can now be found here:
https://github.com/JankoTheBuck/Comp...ter/NSObstacle

Your prescription of the velocity at the outlet is not physical. The velocity is influenced by the body, cannot be exactly as same as at the inlet.
You should set some convection-type BC