Hi.
I m trying to solve hagen poiseuille flow in matlab but not able to get the result.
Can anyone help me…how i can apply Axisymmetric Conditions in LBM.
Hi,
There are now several papers dedicated to axisymmetric LBEs. The original idea (Halliday et al) involves adding a source term to the evolution equation which, in the macroscopic limit, furnishes the “additional” derivatives that arise when the Navier Stokes equations are written in cylindrical coordinates. See also Reis and Phillips. Another approach starts with the Boltzmann equation in cylindrical coordinates and then discretises it in space and time (I forget the authors who proposed this method but a quick Google search should show several articles).
If you don’t get the analytic solution then I’d guess you have one of the following problems: a) your boundary conditions aren’t perfect or b) you’ve got a bug in your code.
Regarding the boundary conditions at the wall, remember you can not satisfy the no-slip condition exactly when using bounce-back with the BGK collision operator. So, if you’re using this approach, you should see a numerical error which decreases like dx^2 (where dx is the grid spacing) and increases with viscosity. You also need to be careful at the line of symmetry because there will be a term in your model which, when evaluated at y=0, gives 0/0. This should be treated with L’Hopital’s formula.
Good luck!
lz=40;
lr=20;
tau=1;
omega=1;
g=0.01102;
force=tau*g;
kvis=0.1667;
ez=[1 0 -1 0 1 -1 -1 1 0];%components of lattice velocities in z-direction
er=[0 1 0 -1 1 1 -1 -1 0];%components of lattice velocities in r-direction
%Initializataion of distributions
is_solid_node=zeros(lr,lz);
for z=1:lz
%is_solid_node(1,z)=1;
is_solid_node(lr,z)=1;
end
disp(is_solid_node)
rho=ones(lr,lz);
f(:,:,1) = (1./9.)*rho;
f(:,:,2) = (1./9.)*rho;
f(:,:,3) = (1./9.)*rho;
f(:,:,4) = (1./9.)*rho;
f(:,:,5) = (1./36.)*rho;
f(:,:,6) = (1./36.)*rho;
f(:,:,7) = (1./36.)*rho;
f(:,:,8) = (1./36.)*rho;
f(:,:,9) = (4./9.)*rho;
% actual time loop
%MAXIT=input(‘Number of iterations’);
for ts=1:1000
for r=1:lr
for z=1:lz
uz(r,z)=0;
ur(r,z)=0;
rho(r,z)=0;
if ~ is_solid_node(r,z)
for a=1:9
rho(r,z)=rho(r,z)+f(r,z,a);
uz(r,z)=uz(r,z)+(ez(a).*f(r,z,a));
ur(r,z)=ur(r,z)+(er(a).*f(r,z,a));
end
uz(r,z)=uz(r,z)/rho(r,z);
ur(r,z)=ur(r,z)/rho(r,z);
end
end
end
% Equilibrium Distribution Function.
f1=3.;
f2=9./2.;
f3=3./2.;
for r=1:lr
for z=1:lz
if ~ is_solid_node(r,z)
rt0 = (4./9. )*rho(r,z); %for rest particles 9
rt1 = (1./9. )*rho(r,z); %for 1 to 4
rt2 = (1./36.)*rho(r,z); %for 5 to 8
ueqz = uz(r,z)+tau*g;
ueqr = ur(r,z);
uzsq = ueqz * ueqz;
ursq = ueqr * ueqr;
uzur5 = ueqz + ueqr;
uzur6 = -ueqz + ueqr;
uzur7 = -ueqz -ueqr;
uzur8 = ueqz -ueqr;
usq = uzsq + ursq;
feq(r,z,9) = rt0*( 1 - f3*usq);
feq(r,z,1) = rt1*( 1 + f1*ueqz + f2*uzsq - f3*usq);
feq(r,z,2) = rt1*( 1 + f1*ueqr + f2*ursq - f3*usq);
feq(r,z,3) = rt1*( 1 - f1*ueqz + f2*uzsq - f3*usq);
feq(r,z,4) = rt1*( 1 - f1*ueqr + f2*ursq - f3*usq);
feq(r,z,5) = rt2*( 1 + f1*uzur5 + f2*uzur5*uzur5 - f3*usq);
feq(r,z,6) = rt2*( 1 + f1*uzur6 + f2*uzur6*uzur6 - f3*usq);
feq(r,z,7) = rt2*( 1 + f1*uzur7 + f2*uzur7*uzur7 - f3*usq);
feq(r,z,8) = rt2*( 1 + f1*uzur8 + f2*uzur8*uzur8 - f3*usq);
end
end
end
% at r=0 ,si1=0,si2=0
% compute from r=1 to lr
for r=1:19
for z=1:lz
rt0 = (4./9. )*rho(r,z); %for rest particles 9
rt1 = (1./9. )*rho(r,z); %for 1 t0 4
rt2 = (1./36.)*rho(r,z); %for 5 to 8
s91(r,z,9)=-((rt0*ur(r,z))/r);
s11(r,z,1)=-((rt1*ur(r,z))/r);
s21(r,z,2)=-((rt1*ur(r,z))/r);
s31(r,z,3)=-((rt1*ur(r,z))/r);
s41(r,z,4)=-((rt1*ur(r,z))/r);
s51(r,z,5)=-((rt2*ur(r,z))/r);
s61(r,z,6)=-((rt2*ur(r,z))/r);
s71(r,z,7)=-((rt2*ur(r,z))/r);
s81(r,z,8)=-((rt2*ur(r,z))/r);
w1=(4/9);
w2=(1/9);
w3=(1/36);
n=-(ur(r,z)/r);
s12(r,z,1)=(3*w2)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,1)-feq(r,z,1))*1-3*ur(r,z)*omega*(f(r,z,1)-feq(r,z,1))*0-(rho(r,z)*ur(r,z).^2)/r)-((6*kvis)/(6*kvis+1)*(f(r,z,1)-feq(r,z,1))*0+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))+0*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,1)-feq(r,z,1)*0)));
s22(r,z,2)=(3*w2)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,2)-feq(r,z,2))*0*1-3*ur(r,z)*omega*(f(r,z,2)-feq(r,z,2))*1-(rho(r,z)*ur(r,z).^2)/r)-0*((6*kvis)/(6*kvis+1)*(f(r,z,2)-feq(r,z,2))*0+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))+1*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,2)-feq(r,z,2)*1)));
s32(r,z,3)=(3*w2)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,3)-feq(r,z,3))*-1*1-3*ur(r,z)*omega*(f(r,z,3)-feq(r,z,3))*0-(rho(r,z)*ur(r,z).^2)/r)+((6*kvis)/(6*kvis+1)*(f(r,z,3)-feq(r,z,3))*0+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))+0*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,3)-feq(r,z,3)*0)));
s42(r,z,4)=(3*w2)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,4)-feq(r,z,4))*0*1-3*ur(r,z)*omega*(f(r,z,4)-feq(r,z,4))*-1-(rho(r,z)*ur(r,z).^2)/r)-0*((6*kvis)/(6*kvis+1)*(f(r,z,4)-feq(r,z,4))*0+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))+-1*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,4)-feq(r,z,4)*0)));
s52(r,z,5)=(3*w2)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,5)-feq(r,z,5))*1-3*ur(r,z)*omega*(f(r,z,5)-feq(r,z,5))*1-(rho(r,z)*ur(r,z).^2)/r)-((6*kvis)/(6*kvis+1)*(f(r,z,5)-feq(r,z,5))*1+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))+1*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,5)-feq(r,z,5)*1)));
s62(r,z,6)=(3*w3)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,6)-feq(r,z,6))*1-3*ur(r,z)*omega*(f(r,z,6)-feq(r,z,6))*1-(rho(r,z)*ur(r,z).^2)/r)+((6*kvis)/(6*kvis+1)*(f(r,z,6)-feq(r,z,6))*-1+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))+1*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,6)-feq(r,z,6)*1)));
s72(r,z,7)=(3*w3)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,7)-feq(r,z,7))*1-3*ur(r,z)*omega*(f(r,z,7)-feq(r,z,7))*1-(rho(r,z)*ur(r,z).^2)/r)+((6*kvis)/(6*kvis+1)*(f(r,z,7)-feq(r,z,7))*1+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))-1*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,7)-feq(r,z,7)*1)));
s82(r,z,8)=(3*w3)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,8)-feq(r,z,8))*1-3*ur(r,z)*omega*(f(r,z,8)-feq(r,z,8))*1-(rho(r,z)*ur(r,z).^2)/r)-((6*kvis)/(6*kvis+1)*(f(r,z,8)-feq(r,z,8))*-1+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))-1*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,8)-feq(r,z,8)*1)));
s92(r,z,9)=(3*w1)/r*(1/2*(rho(r,z)*uz(r,z)*n-((3*ur(r,z)*omega)/2)*(f(r,z,9)-feq(r,z,9))*0-3*ur(r,z)*omega*(f(r,z,9)-feq(r,z,9))*0-(rho(r,z)*ur(r,z).^2)/r)-0*((6*kvis)/(6*kvis+1)*(f(r,z,9)-feq(r,z,9))*0+(rho(r,z)/6)*n-rho(r,z)*uz(r,z)*ur(r,z))+0*(1-12*kvis)*(1/(2*(1+6*kvis))*(f(r,z,9)-feq(r,z,9)*0)));
end
end
%collision step
s(r,z,1)=s11(r,z,1)+s12(r,z,1);
s(r,z,2)=s21(r,z,2)+s22(r,z,2);
s(r,z,3)=s31(r,z,3)+s32(r,z,3);
s(r,z,4)=s41(r,z,4)+s42(r,z,4);
s(r,z,5)=s51(r,z,5)+s52(r,z,5);
s(r,z,6)=s61(r,z,6)+s62(r,z,6);
s(r,z,7)=s71(r,z,7)+s72(r,z,7);
s(r,z,8)=s81(r,z,8)+s82(r,z,8);
s(r,z,9)=s91(r,z,9)+s92(r,z,9);
for r=1:lr
for z=1:lz
if~is_solid_node(r,z)
for a=1:9
f(r,z,a)=f(r,z,a)-((f(r,z,a)-feq(r,z,a))*omega)+s(r,z,a);
end
end
end
end
%bounceback BCS
for r=1:lr
for z=1:lz
if is_solid_node(r,z)
temp = f(r,z,1);
f(r,z,1) = f(r,z,3);
f(r,z,3) = temp;
temp = f(r,z,2);
f(r,z,2) = f(r,z,4);
f(r,z,4) = temp;
temp = f(r,z,5);
f(r,z,5) = f(r,z,7);
f(r,z,7) = temp;
temp = f(r,z,6);
f(r,z,6) = f(r,z,8);
f(r,z,8) = temp;
end
end
end
%streaming(propagation) step
% C6 C2 C5 ^ r
% \ | / |
% C3-C0-C1 |
% / | \ |
% C7 C4 C8 -----> z
% (in,jp) (i,jp) (ip,jp) ^ r(jp)
% \ | / |
% (in,j)-(i,j)-(ip,j) |
% / | \ |
% (in,jn) (i,jn) (ip,jn) (in)<-----(i,j) -----> z(ip)
% |
% |
% |
% (jn)
for r=1:lr
if r>1
jn=r-1;
else
jn=lr;
end
if r<lr
jp=r+1;
else
jp=1;
end
for z=1:lz
if z>1
in=z-1;
else
in=lz;
end
if z<lz
ip=z+1;
else
ip=1;
end
ftemp(r,z,9) = f(r,z,9);
ftemp(r,ip,1) = f(r,z,1);
ftemp(jp,z,2) = f(r,z,2);
ftemp(r,in,3) = f(r,z,3);
ftemp(jn,z,4) = f(r,z,4);
ftemp(jp,ip,5) = f(r,z,5);
ftemp(jp,in,6) = f(r,z,6);
ftemp(jn,in,7) = f(r,z,7);
ftemp(jn,ip,8) = f(r,z,8);
end
end
f(:,:,:)=ftemp(:,:,:);
end
I hope this doesn’t sound too grumpy or rude because I like to help, where possible, but I think you’d learn more (and feel more rewarded) if you consider any advice and comments given on this forum and, with this in mind, go through your code, bit by bit, by yourself before sending an entire code to the forum.
Indeed, you could alter the number of grid points to see if your predictions are better when you use more nodes. Likewise, you can alter your relaxation parameter. Make sure all your explicit finite difference terms are correct and consistent, especially when you are at, or near, the boundaries . Also, be careful with the force term you are using to drive the flow.
Note also that the analytic solution is u®=U(1-r^2/a^2), where a is the radius of the pipe and r is the radial coordinate. The maximum velocity is U=Ga^2/(4nu*rho), where nu is the kinematic viscosity, rho is the density and G is the pressure gradient (force) driving the flow. Your inflow condition or driving force should reflect this.
Good luck!
ps
We have all had moments of extreme frustration when trying to get our codes to behave nicely - I’ve been close to throwing my laptop out of the window twice this week already!