Hi,
I am doing the Poiseulle flow with pressure boundary conditions. My code is as follows. It is blowing up. I could not figure out where I am going wrong. Anybody please check the code.
TAU=0.75;
OMEGA=1/TAU;
XMAX=17; % Mesh size in x-direction
YMAX=5; % Mesh size in y-direction
pin = 1.1;
pout = 1;
cx=[1 0 -1 0 1 -1 -1 1 0]; % components of lattice velocities in x-direction
cy=[0 1 0 -1 1 1 -1 -1 0]; % components of lattice velocities in y-direction
jx=zeros(XMAX,YMAX); % creating an array for storing momentum in x-direction
jy=zeros(XMAX,YMAX); % creating an array for storing momentum in y-direction
rho=ones(XMAX,YMAX); % creating an array for storing densities
f=zeros(XMAX,YMAX); % creating an array for storing distributions
fprop=zeros(XMAX,YMAX); % creating an array for storing propagating distributions
image(1:XMAX,1:YMAX)=0; % Fluid
image(1:XMAX,1)= 1; % Wall
image(1:XMAX,YMAX)= 1; % Wall
image(1,2:YMAX-1)= 9; % Inlet
image(XMAX,2:YMAX-1)= 10; % Outlet
MAXIT=input('Enter the number of iterations:');
%Initializataion of distributions
u=jx./rho;
v=jy./rho;
f = function_feq(rho, u, v);
%Initialization completed
%Assign first fprop to equilibrium distributions.
fprop=f;
feq=zeros(XMAX,YMAX,9);
for iter=1:MAXIT
feq = function_feq(rho, u, v); % Gives Equilibrium Distribution function..
fprop = (1.0-OMEGA).* f + OMEGA.* feq;
for x=1:XMAX
for y=1:YMAX
for i=1:9
xnew = x + cx(i);
ynew = y + cy(i);
if (xnew>=1 && xnew <= XMAX) && (ynew>=1 && ynew <= YMAX)
f(xnew, ynew, i) = fprop(x, y, i);
end
end
end
end
% Boundary Conditions
for x=1:XMAX
for y=1:YMAX
if image(x,y)==9 % Inlet
u(x,y) = 1 - (f(x,y,9)+f(x,y,2)+f(x,y,4)+2*(f(x,y,3)+f(x,y,6)+f(x,y,7)))/3*pin;
f(x,y,1) = f(x,y,3) + 2/3 * (3*pin*u(x,y));
f(x,y,5) = f(x,y,7) - 1/2 * (f(x,y,2)-f(x,y,4)) + 1/6 * (3*pin*u(x,y));
f(x,y,8) = f(x,y,6) + 1/2 * (f(x,y,2)-f(x,y,4)) + 1/6 * (3*pin*u(x,y));
elseif image(x,y)==10 %Outlet
u(x,y) = (f(x,y,9)+f(x,y,2)+f(x,y,4)+2*(f(x,y,5)+f(x,y,1)+f(x,y,8)))/3*pout - 1;
f(x,y,3) = f(x,y,1) - 2/3 * (3*pout*u(x,y));
f(x,y,6) = f(x,y,8) - 1/2 * (f(x,y,2)-f(x,y,4)) - 1/6 * (3*pout*u(x,y));
f(x,y,7) = f(x,y,5) + 1/2 * (f(x,y,2)-f(x,y,4)) - 1/6 * (3*pout*u(x,y));
elseif image(x,y) == 1 % Wall
temp = f(x,y,2);
f(x,y,2) = f(x,y,4);
f(x,y,4) = temp;
temp = f(x,y,6);
f(x,y,6) = f(x,y,8);
f(x,y,8) = temp;
temp = f(x,y,5);
f(x,y,7) = f(x,y,5);
f(x,y,5) = temp;
end
end
end
rho = f(:,:,1)+f(:,:,2)+f(:,:,3)+f(:,:,4)+f(:,:,5)+f(:,:,6)+f(:,:,7)+f(:,:,8)+f(:,:,9);
jx = f(:,:,1)-f(:,:,3)+f(:,:,5)-f(:,:,6)-f(:,:,7)+f(:,:,8);%Distributions multiplied by lattice velocities in x-directions
jy = f(:,:,2)-f(:,:,4)+f(:,:,5)+f(:,:,6)-f(:,:,7)-f(:,:,8);%Distributions multiplied by lattice velocities in y-directions
u = jx./rho;
v = jy./rho;
end
Thanks,
Narender.