# Poiseulle flow with Pressure boundary conditions

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.

Hello,

Thanks,
Narender.

Hi All,

I am getting velocity in y direction to the order of 1.0e-3 . Also the y velocity in next corner node (Eg: (2,2) for corner node (1,1)) is very high 0.32e-3. These values are after steady state. It is clear that these velocities are generated from the corner nodes. Please help on this.

Thanks,
Narender

Hi All,

This code is working after I incorporate the corner node boundary conditions. The solution is not machine accurate as I am using the standard compressible LB model.

Thanks,
Narender