Ocsillatory flow

Lattice Boltzmann Method LBM: Theory
#1

Hi friends
I have simulated ocsillatory flow in channel with sinusoidal pressure gradient between inlet and outlet. BC for inlet:( from zou-He pressure BC):

rho(:,in,co)=1+0.001sin(OMEGAcycle);
ux(:,in,colin)=1 - 1 ./ (rho(:,in,colin)) .* ( …
sum(fIn([1,3,5],in,colin)) + 2*sum(fIn([4,8,7],in,colin)) );
uy(:,in,colin)=0;

and BC for outlet:
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;

where OMEGA=2*pi/1000, and Lx=200 and Ly=20
I have used incompressible model D2Q9 for feq.(He and Lou 1997)
I have used microscopic BC for inlet and outlet the same as available code " cylider.m"
But the velocity profile in the time is not axicymmetric and dose not match with analytical solution for womersley flow.

Does anyone have experience on this problem?
Is there any mistake in my BC?

#2

Hello,

if you have not read the following article yet, please do:
Wang, He, Tang, Tao. Simulation of two-dimensional oscillating flow using the lattice Boltzmann method. Int. J. Modern Phys. C 17 (5), 615-630 (2006).
Have you tried to increase the time period (reduce the oscillation frequency)? Maybe the time scales for oscillation and sound waves are not well separated.

Best,
Timm

#3

Hi
thanks Timm for your quick reply.
I use these parameters the same as article(Lattice boltzmann model for incompressile Navier stokes Eq,By Luo,He 1997)
and My doubt is on the implementing inlet BC.
I increased T but the problem remain: magnitude of velocity in (-x and x) direction(x is direction of channel centeline) are not the equal!!

#4

Then you should carefully recheck your implementation.
How does the velocity profile look like? Is your mass conserved? Any hints which could tell us where the error could be located?

#5

I cheked the velocity profile at inlet and outlet, at each time step they are the same. and rho variation between (1(±)0.001) inlet and 1 at outlet.
which other parameter should I check for mass conservation ?? I don’t know how to insert my result. but I can say profiles are nearly the same as article 'wang" that you mentioned.but not axisymetric!

#6

Please describe the geometry you use? Is it 3D? Circular or rectangular cross-section?

#7

thanks Timm
2D rectangular channel, I set parabolic velocity for initial condition,

for cycle = 1:maxT
% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux = reshape ( (cx * reshape(fIn,9,lxly)), 1,lx,ly) ./rho;
uy = reshape ( (cy * reshape(fIn,9,lxly)), 1,lx,ly) ./rho;

% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
  % Inlet: pressure BC
rho(:,in,col)=1+0.001*sin(OMEGA*cycle); 
ux(:,in,col)=1 - 1 ./ (rho(:,in,col)) .* ( ...
   sum(fIn([1,3,5],in,col)) + 2*sum(fIn([4,8,7],in,col));
uy(:,in,col)=0;

% Outlet: pressure BC
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/3rho(:,in,col).ux(:,in,col);
fIn(6,in,col) = fIn(8,in,col) + 1/2(fIn(5,in,col)-fIn(3,in,col)) …
+ 1/2rho(:,in,col).uy(:,in,col) …
+ 1/6rho(:,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/6rho(:,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(incomresible model for feq)
for i=1:9
cu = 3*(cx(i)ux+cy(i)uy);
fEq(i,:,:slight_smile: =rho.t(i)+ rho0. t(i). …
( cu + 1/2(cu.cu) - 3/2(ux.^2+uy.^2) );
fOut(i,:,:slight_smile: = fIn(i,:,:slight_smile: - 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

thanks a lot

#8

I do not understand what you mean by axisymmetry in this case. You mean that the magnitude of the velocity at t and t + T/2 (T being the oscillation period) is not the same?

#9

hi
excuse me for my unclear explain.
Yes, I Exactly mean what you mentioned: magnitude of the velocity at t and t + T/2 (T being the oscillation period) is not the same.

#10

Hi,

I am new to PALABOS. While trying to run the examples in the showcases, I found that there are a few examples which compile but do not run. Eg. Womersley. Have anyone tried this already?
Let me know.

Thanks in advance.

Sanjib