Strange boundary condition problem (poiseuille flow)

jmazo · April 2, 2010, 10:47pm

Hi there
I´m trying to simulate a simple 2d channel using 2dq9 bgk

I took as guide the code available on this site (cylinder.m) and rewrote it in a way I could understand it.

If I simulate the flow using the inlet boundary condition the poiseuille flow it works as expected., but If I use something different like a constant speed
It shows an strange profile (NORTH WEST) and (SOUTH EAST), I´ve being trying to debug the problem but I just cant get it working right.

If you want to see what I´m talking about just change

  vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);   
  %vxin = 0.04;

for
% vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);
vxin = 0.04;

on initialaizing the fluid and when setting the inlet boundary condition

thanks in advance.

% mycavity
clear
iMax = 100000000000;

%grid size

lx = 400;
ly = 100;

ux = 0.04; % horizontal lid velocity
uy = 0; % vertical lid velocity
rhoOut = 1.0; %outlet pressure
obst_r = ly/10+1;
Re = 100; % Reynolds number

uMax = 0.1; % maximum velocity of Poiseuille inflow
nu = uMax * 2.obst_r / Re; % kinematic viscosity
tau = 1. / (3nu+1./2.);

% D2Q9 LATTICE CONSTANTS

w0 = 4/9;
wm = 1/9;
wd = 1/36;

fc1 = 3.0;
fc2 = 9/2;
fc3 = 3/2;

%Arrays
% rho
rho = zeros(ly,lx);
u_x = zeros(ly,lx);
u_y = zeros(ly,lx);
vor = zeros(ly,lx);

% Equilibrium functions arrays
f0 = zeros(ly,lx);
f1 = zeros(ly,lx);
f2 = zeros(ly,lx);
f3 = zeros(ly,lx);
f4 = zeros(ly,lx);
f5 = zeros(ly,lx);
f6 = zeros(ly,lx);
f7 = zeros(ly,lx);
f8 = zeros(ly,lx);

% temporal storage
tmpf0 = zeros(ly,lx);
tmpf1 = zeros(ly,lx);
tmpf2 = zeros(ly,lx);
tmpf3 = zeros(ly,lx);
tmpf4 = zeros(ly,lx);
tmpf5 = zeros(ly,lx);
tmpf6 = zeros(ly,lx);
tmpf7 = zeros(ly,lx);
tmpf8 = zeros(ly,lx);

%---------------Create domain--------------------
% 1 = fluid
% 0 = Wall

D = ones(ly,lx);

for x=1:lx;
D(1,x) = 0; %top wall
D(ly,x)= 0; %bottom wall
end

roout = rhoOut;

% set the fluid to rest (rho=0, u=0) ==> fi(i) = t(i)
for y=1:ly;
L = ly-2;
y_phys = y-1.5;
for x=1:lx;

  vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);  
    
  
  %vxin = 0.04;
  
  f0(y,x)  = w0 * rhoOut * (1                             - (3/2)*vxin*vxin);
  f1(y,x)  = wm * rhoOut * (1 + 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
  f2(y,x)  = wm * rhoOut * (1                             - (3/2)*vxin*vxin);
  f3(y,x)  = wm * rhoOut * (1 - 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
  f4(y,x)  = wm * rhoOut * (1                             - (3/2)*vxin*vxin);
  f5(y,x)  = wd * rhoOut * (1 + 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
  f6(y,x)  = wd * rhoOut * (1 - 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
  f7(y,x)  = wd * rhoOut * (1 - 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
  f8(y,x)  = wd * rhoOut * (1 + 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
     
end

end

for k=1:iMax% MAIN LOOP

%calculate macroscopic variables
%rho, u_x, u_y
for i=1:ly;
 for j=1:lx;
     
     if D(i,j) ~= 0;
         rom =  f0(i,j) + f1(i,j)+ f2(i,j)+ f3(i,j)+ f4(i,j)+ f5(i,j)+ f6(i,j)+ f7(i,j)+ f8(i,j);
        
         uxm =            f1(i,j)         - f3(i,j)         + f5(i,j)- f6(i,j)- f7(i,j)+ f8(i,j);
         uym =                   + f2(i,j)         - f4(i,j)+ f5(i,j)+ f6(i,j)- f7(i,j)- f8(i,j);
         
         rho(i,j) = rom;
         u_x(i,j) = uxm / rom;
         u_y(i,j) = uym / rom;
         
         
     end
 end
end 


%INLET
for i=2:ly-1
    L = ly-2; 
    y_phys = i-1.5;
    vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);
    %vxin = 0.04;


    
    %macro (FLUX)
    
    u_x(i,1) = vxin;
    u_y(i,1) = uy;
 
    %micro
    rho0 = ((f0(i,1) + f2(i,1) + f4(i,1)) + 2*(f3(i,1) + f7(i,1) + f6(i,1))) / (1 - vxin);
    rho(i,1) = rho0;
    
    ru = rho0 * vxin;
    
    f1(i,1) = f3(i,1) + (2/3)*ru;
    f5(i,1) = f7(i,1) + (1/6)*ru + (1/2)*(f2(i,1) - f4(i,1));
    f8(i,1) = f6(i,1) + (1/6)*ru + (1/2)*(f4(i,1) - f2(i,1));
     
end


%OUTLET
for i=2:ly-1
    %macro (PRESSURE)
    rho(i,lx) = rhoOut;
    u_y(i,lx) = uy; 
    %micro
    
    ux0 = -1 + ((f0(i,lx) + f2(i,lx) + f4(i,lx))  + 2*(f1(i,lx) + f5(i,lx) + f8(i,lx)))/ rhoOut;
    
    ru = rhoOut * ux0;
    
    f3(i,lx) = f1(i,lx) - (2/3) * ru;
    f7(i,lx) = f5(i,lx) - (1/6) * ru + (1/2)*(f2(i,lx) - f4(i,lx));
    f6(i,lx) = f8(i,lx) - (1/6) * ru + (1/2)*(f4(i,lx) - f2(i,lx));
     
end




%=====================================================================
%============== COLLISION ============================================
%=====================================================================


 %Compute equilibrium functions
 for i=1:ly;
   for j=1:lx;
         
        %get needed values
        rq = rho(i,j);
        vx = u_x(i,j);
        vy = u_y(i,j);
         
        rt0 = w0 * rq;
        rt1 = wm * rq;
        rt2 = wd * rq;
        
        
        v_sq_term = fc3*(vx*vx + vy*vy);
        
        % Evaluate the local equilibrium f values in all directions
        f0eq = rt0 * (1 - v_sq_term);
        f1eq = rt1 * (1 + fc1*vx +        fc2*vx*vx -                v_sq_term);
        f2eq = rt1 * (1 + fc1*vy +        fc2*vy*vy -                v_sq_term);
        f3eq = rt1 * (1 - fc1*vx +        fc2*vx*vx -                v_sq_term);
        f4eq = rt1 * (1 - fc1*vy +        fc2*vy*vy -                v_sq_term);
        f5eq = rt2 * (1 + fc1*(vx + vy) + fc2*(vx + vy)*(vx + vy) -  v_sq_term);
        f6eq = rt2 * (1 + fc1*(-vx + vy)+ fc2*(-vx + vy)*(-vx + vy)- v_sq_term);
        f7eq = rt2 * (1 + fc1*(-vx - vy)+ fc2*(-vx - vy)*(-vx - vy)- v_sq_term);
        f8eq = rt2 * (1 + fc1*(vx - vy) + fc2*(vx - vy)*(vx - vy) -  v_sq_term);
        


        tmpf0(i,j) = f0(i,j) - (f0(i,j) - f0eq)*tau;
        tmpf1(i,j) = f1(i,j) - (f1(i,j) - f1eq)*tau;
        tmpf2(i,j) = f2(i,j) - (f2(i,j) - f2eq)*tau;
        tmpf3(i,j) = f3(i,j) - (f3(i,j) - f3eq)*tau;
        tmpf4(i,j) = f4(i,j) - (f4(i,j) - f4eq)*tau;
        tmpf5(i,j) = f5(i,j) - (f5(i,j) - f5eq)*tau;
        tmpf6(i,j) = f6(i,j) - (f6(i,j) - f6eq)*tau;
        tmpf7(i,j) = f7(i,j) - (f7(i,j) - f7eq)*tau;
        tmpf8(i,j) = f8(i,j) - (f8(i,j) - f8eq)*tau;

     
    end
 end
     
 
 % Bounce-Back
   for i=1:ly;
     for j=1:lx;
         if D(i,j) == 0;
        f1old = f1(i,j);
        f2old = f2(i,j);
        f3old = f3(i,j);
        f4old = f4(i,j);
        f5old = f5(i,j);
        f6old = f6(i,j);
        f7old = f7(i,j);
        f8old = f8(i,j);

        tmpf1(i,j) = f3old;
        tmpf2(i,j) = f4old;
        tmpf3(i,j) = f1old;
        tmpf4(i,j) = f2old;
        tmpf5(i,j) = f7old;
        tmpf6(i,j) = f8old;
        tmpf7(i,j) = f5old;
        tmpf8(i,j) = f6old;
         end
     end
   end
 
 
 %=====================================================================
 %=====================================================================
 
% Streaming
for i=1:ly
    in = i - 1;
    is = i + 1;
    
    if i == 1
        in = 1;
    end
    
    if i == ly
        is = ly;
    end
        
    
    
    
    for j=1:lx
        je = j + 1;
        jw = j - 1;
        
        
        if j == 1
            jw = 1;
        end

        if j == lx
            je = lx;
        end
        
         f0(i,j) = tmpf0(i,j);
         f1(i,je) = tmpf1(i,j);
         f2(in,j) = tmpf2(i,j);
         f3(i,jw) = tmpf3(i,j);
         f4(is,j) = tmpf4(i,j);
         f5(in,je) = tmpf5(i,j);
         f6(in,jw) = tmpf6(i,j);
         f7(is,jw) = tmpf7(i,j);
         f8(is,je) = tmpf8(i,j);
        
    end

  
end 
 
 

% VISUALIZATION
if (mod(k,2)==0)
    u = sqrt(u_x.^2+u_y.^2);

    imagesc(u);
    %contour(u);

    axis equal off; drawnow
end

end

jmazo · April 3, 2010, 2:26am

Hi there
I have modified the cylinder.m example and set a inlet speed inSpeedx of 0.04;
the same phenomenon is observed!
it looks like a boundary layer, but I dont think it is one, but rather problem with the inlet boundary.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cylinder.m: Channel flow past a cylinderical
% obstacle, using a LB method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample in Matlab
% Copyright © 2006-2008 Jonas Latt
% Address: EPFL, 1015 Lausanne, Switzerland
% E-mail: jonas@lbmethod.org
% Get the most recent version of this file on LBMethod.org:
% http://www.lbmethod.org/_media/numerics:cylinder.m
%
% Original implementaion of Zou/He boundary condition by
% Adriano Sciacovelli (see example “cavity.m”)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% You should have received a copy of the GNU General Public
% License along with this program; if not, write to the Free
% Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
% Boston, MA 02110-1301, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% GENERAL FLOW CONSTANTS
lx = 400; % number of cells in x-direction
ly = 100; % number of cells in y-direction
obst_x = lx/5+1; % position of the cylinder; (exact
obst_y = ly/2+3; % y-symmetry is avoided)
obst_r = ly/20+1; % radius of the cylinder
uMax = 0.1; % maximum velocity of Poiseuille inflow
Re = 100; % Reynolds number
nu = uMax * 2.obst_r / Re; % kinematic viscosity
omega = 1. / (3nu+1./2.); % relaxation parameter
%omega = 0.51;
maxT = 400000; % total number of iterations
tPlot = 5; % cycles

% D2Q9 LATTICE CONSTANTS
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
col = [2:(ly-1)];
in = 1; % position of inlet
out = lx; % position of outlet

%init Speeds
inSpeedx = 0.04;
inLetCol = zeros(1, ly-2);
inLetCol = inLetCol + inSpeedx;

[y,x] = meshgrid(1:ly,1:lx); % get coordinate of matrix indices

obst = … % Location of cylinder
(x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2;
obst(:,[1,ly]) = 1; % Location of top/bottom boundary
bbRegion = find(obst); % Boolean mask for bounce-back cells

% INITIAL CONDITION: Poiseuille profile at equilibrium
L = ly-2; y_phys = y-1.5;
%ux = 4 * uMax / (L*L) * (y_phys.L-y_phys.y_phys);
ux = zeros(lx,ly);
ux = ux + inSpeedx;
uy = zeros(lx,ly);
rho = 1;
for i=1:9
cu = 3(cx(i)ux+cy(i)uy);
fIn(i,:, = rho . t(i) . …
( 1 + cu + 1/2(cu.cu) - 3/2(ux.^2+uy.^2) );
end

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT

% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux  = reshape ( (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
uy  = reshape ( (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
  
% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
  % Inlet: Poiseuille profile
y_phys = col-1.5;
%ux(:,in,col) = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);

ux(:,in,col) = inLetCol;
uy(:,in,col) = 0;
rho(:,in,col) = 1 ./ (1-ux(:,in,col)) .* ( ...
    sum(fIn([1,3,5],in,col)) + 2*sum(fIn([4,7,8],in,col)) );
  % Outlet: Constant pressure
rho(:,out,col) = 1;
ux(:,out,col) = -1 + 1 ./ (rho(:,out,col)) .* ( ...
    sum(fIn([1,3,5],out,col)) + 2*sum(fIn([2,6,9],out,col)) );
uy(:,out,col)  = 0;

% MICROSCOPIC BOUNDARY CONDITIONS: INLET (Zou/He BC)
fIn(2,in,col) = fIn(4,in,col) + 2/3*rho(:,in,col).*ux(:,in,col); 
fIn(6,in,col) = fIn(8,in,col) + 1/2*(fIn(5,in,col)-fIn(3,in,col)) ... 
                                + 1/2*rho(:,in,col).*uy(:,in,col) ...
                                + 1/6*rho(:,in,col).*ux(:,in,col); 
fIn(9,in,col) = fIn(7,in,col) + 1/2*(fIn(3,in,col)-fIn(5,in,col)) ... 
                                - 1/2*rho(:,in,col).*uy(:,in,col) ...
                                + 1/6*rho(:,in,col).*ux(:,in,col); 

% MICROSCOPIC BOUNDARY CONDITIONS: OUTLET (Zou/He BC)
fIn(4,out,col) = fIn(2,out,col) - 2/3*rho(:,out,col).*ux(:,out,col); 
fIn(8,out,col) = fIn(6,out,col) + 1/2*(fIn(3,out,col)-fIn(5,out,col)) ... 
                                  - 1/2*rho(:,out,col).*uy(:,out,col) ...
                                  - 1/6*rho(:,out,col).*ux(:,out,col); 
fIn(7,out,col) = fIn(9,out,col) + 1/2*(fIn(5,out,col)-fIn(3,out,col)) ... 
                                  + 1/2*rho(:,out,col).*uy(:,out,col) ...
                                  - 1/6*rho(:,out,col).*ux(:,out,col); 

% COLLISION STEP
for i=1:9
   cu = 3*(cx(i)*ux+cy(i)*uy);
   fEq(i,:,:)  = rho .* t(i) .* ...
                   ( 1 + cu + 1/2*(cu.*cu)  - 3/2*(ux.^2+uy.^2) );
   fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:));
end

% OBSTACLE (BOUNCE-BACK)
for i=1:9
     fOut(i,bbRegion) = fIn(opp(i),bbRegion);
end

% STREAMING STEP
for i=1:9
   fIn(i,:,:) = circshift(fOut(i,:,:), [0,cx(i),cy(i)]);
end

% VISUALIZATION
if (mod(cycle,tPlot)==1)
    u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
    %p = reshape(rho/3,lx,ly);
    u(bbRegion) = nan;
    imagesc(u');
    %p(bbRegion) = nan;
    %imagesc(p');
    axis equal off; drawnow
end

end

ydt · April 12, 2010, 5:47am

Hi
can inlet boundary conditon be changed from velocity to pressure?, at the same way, can initial condition be changed the pressure ?.

ydt · April 12, 2010, 5:53am

ux = 4 * uMax / (L*L) * (y_phys.*L-y_phys.y_phys) is how to be understand. I canot understand means "ux = 4 * uMax / (LL) * (y_phys.*L-y_phys.*y_phys)". can initial condition and boundary condition that is pressure and pressure different be for example. thank a lot in advance.

jmazo · April 12, 2010, 3:28pm

ydt

that formula generates de poiseuille speed profile at the inlet. not an linear speed lets say 0,04. it generates a parabolic profile. that depends on the width of the channel.

jmazo · April 14, 2010, 4:10pm

I just want to mention that if you implement the corners as stated by GoncaolSilva, it works
it matches the analytical solution.