[infikamal5]

Hi Everyone,

I have written a 2D cavity flow problem using staggered grid finite volume method with pressure correction and iterative pressure solver


I tried to change the boundary condition but always ended up with error

Could some please help me to change boundary condition so that the code work as 2D channel Flow with iterative pressure solver.

==========================================
%%%%%%%%%%%%%%%%%%%%

%%% Finite volume method %%%%%%%%%%
%%%% SOR technique for presure %%%%


%%%%%%%%%%% Creating the mesh %%%%%%%%%


clear all;
close all;
clc;

nx = 10 ; %% number of elements
ny = 10 ; %% number of elements

L = 1;
H = 1;


hx = 1/nx;

mu = 10^-6 ;

dt = 0.005 ;

P = zeros(nx+2,ny+2);
U = zeros(ny+2,nx+1);
V = zeros(ny+1,nx+2);

ut = zeros(ny+2,nx+1);
vt = zeros(ny+1,nx+2);

uu = zeros(ny+1,nx+1);
vv = zeros(ny+1,nx+1);
w = zeros(ny+1,nx+1);

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


c(3:ny,2) = 1/3; % Boundary points
c(2,3:nx) = 1/3;
c(3:ny,nx+1) = 1/3;
c(ny+1,3:ny) = 1/3;


c(2,2) =1/2; %corner points
c(2,ny+1) =1/2;
c(nx+1,2) =1/2;
c(nx+1,ny+1) =1/2;

time = 0;

un = 1;
us = 0;
ve = 0;
vw = 0;

beta = 1.2;
mu = 0.1 ;

for is = 1:100



V(1:ny+1,1)=2*vw-V(1:ny+1,2);
V(1:ny+1,nx+2)=2*ve-V(1:ny+1,nx+1);
U(1,1:nx+1)=2*us-U(2,1:nx+1);
U(ny+2,1:nx+1)=2*un-U(ny+1,1:nx+1);



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

ut(j,i) = U(j,i) + dt * ( -(0.25/hx) * ( U(j,i+1) + U(j,i) ).^2 ...
- ( U(j,i-1) + U(j,i) ).^2 ...
+ ( U(j+1,i) + U(j,i) ) * ( V(j,i) + V(j,i+1) ) ...
- ( U(j-1,i) + U(j,i) ) * ( V(j-1,i) + V(j-1,i+1)) ...
+ (mu) * (1/hx^2) * ( U(j,i+1) + U(j,i-1) + U(j+1,i) + U(j-1,i) - 4*U(j,i)));

end
end



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

vt(j,i) = V(j,i) + dt * ( -(0.25/hx) * ( V(j+1,i) + V(j,i) ).^2 ...
- ( V(j-1,i) + V(j,i) ).^2 ...
+ ( V(j,i+1) + V(j,i) ) * ( U(j+1,i) + U(j,i) ) ...
- ( V(j,i-1) + V(j,i) ) * ( U(j,i-1) + U(j+1,i-1) ) ...
+ (mu) * (1/hx^2) * ( V(j,i+1) + V(j,i-1) + V(j+1,i) + V(j-1,i) - 4*V(j,i)));

end
end



for it=1:100

for i=2:nx+1,
for j=2:ny+1 % solve for pressure
P(j,i)= beta*c(j,i)...
* (P(j,i+1)+P(j,i-1)+P(j+1,i)+P(j-1,i)...
-(hx/dt)*(ut(j,i)-ut(j,i-1)+vt(j,i)-vt(j-1,i))) +(1-beta)*P(j,i) ;
end
end
end
%
P(1,2:nx+1) = P(2,2:nx+1); % Bottom (dp/dy=0)
P(ny+2,2:nx+1) = P(ny+1,2:nx+1); % Top (dp/dy=0)
P(2:ny+1,1) = P(2:ny+1,2); % Left (dp/dx=0)
P(2:ny+1,nx+2) = -P(2:ny+1,nx+1);


% correct the velocity
U(2:ny+1,2:nx)=...
ut(2:ny+1,2:nx)-(dt/hx)*(P(2:ny+1,3:nx+1)-P(2:ny+1,2:nx));
V(2:ny,2:nx+1)=...
vt(2:ny,2:nx+1)-(dt/hx)*(P(3:ny+1,2:nx+1)-P(2:ny,2:nx+1));


uu(1:ny+1,1:nx+1)=0.5*(U(2:ny+2,1:nx+1)+U(1:ny+1,1 :nx+1));
vv(1:ny+1,1:nx+1)=0.5*(V(1:ny+1,2:nx+2)+V(1:ny+1,1 :nx+1));
w(1:ny+1,1:nx+1)=(U(2:ny+2,1:nx+1)-U(1:ny+1,1:nx+1)-...
V(1:ny+1,2:nx+2)+V(1:ny+1,1:nx+1))/(2*hx);


time=time+dt ;

figure(2);
hold off;contourf(w,60),axis equal,pause(0.01)
end