Hello everyone,
I need a simple code (Poiseuil 2D) with the boundary conditions Halfway. I was able to implement (Poiseuil 2D, 3D) with the boundary condition Fullway, but for the boundary conditions Halfway I still have a problem. I count on your kindness. Thank you in Avence.
Hi zarbaya,
Have you had a look at the following post?
http://www.palabos.org/forum/read.php?9,7184
The half-way boundary is explained in simple terms.
If you want, you can post your code with full-way bounce-back and we’ll see how to change it for half-way.
Regards,
ndel
Hello nicolasD and everyone in this forum
I’m struggling with the exact same problem right now, but somehow the velocities on the solid nodes won’t be negative.
The only thing I get is zero velocity in flow direction on the walls just like in a fullway bounceback. Can you please help me turn this to a halfway bounceback?
I also wonder if my body force implementation is correct.
Thanks in advance,
Totomo
P.s.: Here’s an extract from my MATLAB code:
% define matrices and variables
...
% initialization
for a = 1 : 9
uc = (cx(a)*ux(1) + cy(a)*uy(1));
f_prop(:,:,a) = rho * w(a) * ( 1 + uc/cs^2 + uc.^2/(2*cs^4) - (ux.^2+uy.^2)/(2*cs^2) );
end
for t = 1 : runtime
% adding body force
ux = ux + dt * Fx ./ (2*density);
uy = uy + dt * Fy ./ (2*density);
% calculate equilibrium
for a = 1 : 9
uc = ( cx(a)*ux(:,:) + cy(a)*uy(:,:) );
f_eq(:,:,a) = rho * w(a) * ( 1 + uc/cs^2 + uc.^2/(2*cs^4) - (ux(:,:).^2+uy(:,:).^2)/(2*cs^2) );
end
% collision with body force
for i = 1 : ly
for j = 1 : lx
for a = 1 : 9
F_a(i,j,a) = (1 - 1/(2*tau)) * w(a) * ( ( (c(1,a) - ux(i,j)) / cs^2 + ...
(c(1,a) * ux(i,j) + c(2,a) * uy(i,j)) * c(1,a) / cs^4 ) * Fx + ...
( (c(2,a) - uy(i,j)) / cs^2 + (c(1,a) * ux(i,j) + c(2,a) * uy(i,j)) * ...
c(2,a) / cs^4 ) * Fy );
end
end
end
for a = 1 : 9
f_c(:,:,a) = f_prop(:,:,a) - 1/tau * (f_prop(:,:,a) - f_eq(:,:,a)) + dt * F_a(a);
end
% bounceback
for g = 1 : length(row)
temp = f_c(row(g),col(g),7); f_prop(row(g),col(g),7) = f_c(row(g),col(g),5); f_prop(row(g),col(g),5) = temp;
temp = f_c(row(g),col(g),4); f_prop(row(g),col(g),4) = f_c(row(g),col(g),2); f_prop(row(g),col(g),2) = temp;
temp = f_c(row(g),col(g),8); f_prop(row(g),col(g),8) = f_c(row(g),col(g),6); f_prop(row(g),col(g),6) = temp;
end
% propagation / streaming
for i = 1 : ly
if i > 1
in = i-1;
else
in = ly;
end
if i < ly
ip = i+1;
else
ip = 1;
end
for j = 1 : lx
if j > 1
jn = j-1;
else
jn = lx;
end
if j < lx
jp = j+1;
else
jp = 1;
end
ftemp(i,j,9) = f_c(i,j,9);
ftemp(i,jp,1) = f_c(i,j,1);
ftemp(ip,j,4) = f_c(i,j,4);
ftemp(i,jn,3) = f_c(i,j,3);
ftemp(in,j ,2) = f_c(i,j,2);
ftemp(ip,jp,8) = f_c(i,j,8);
ftemp(ip,jn,7) = f_c(i,j,7);
ftemp(in,jn,6) = f_c(i,j,6);
ftemp(in,jp,5) = f_c(i,j,5);
end
end
f_prop(:,:,:)=ftemp(:,:,:);
% density and velocity
density = sum(f_prop,3);
ux = (f_prop(:,:,1)-f_prop(:,:,3)+f_prop(:,:,5)-f_prop(:,:,6)+f_prop(:,:,8)-f_prop(:,:,7))./density;
uy = (f_prop(:,:,2)-f_prop(:,:,4)+f_prop(:,:,5)-f_prop(:,:,8)+f_prop(:,:,6)-f_prop(:,:,7))./density;
% visualization
...
Hi Zahrbaya,
could you share it with us your code?
thank you