Dear Timm,
I read the paper by Guo, and chose the first method of body force implementation. The problem to be solved is (as reported by Sukop):
yDim = 11 lu
Re = 4.4
Thaw = 1/omega = 1.0
Kinematic viscosity = 0.166667
Maximum velocity = 0.1
Gravitational acceleration = g = 1.10210^-3
the collision term is changed to :
if (image(y,x)≠ wall)
f(y,x,i) = (1.0d0 - omega) * f(y,x,i) + omega * feq(y,x,i) + Force(i)
where,
Force(i) = (1- 0.5omega)ω (i){[ 3*(e(i,0)-u(y,x,0)) + 9*e(i,0)*u(y,x,0) + e(i,1)*u(y,x,1))e(i,0) )]F(0) + [ 3(e(i,1)-u(y,x,1)) + 9e(i,0)*u(y,x,0) + e(i,1)*u(y,x,1))e(i,1) )]F(1)}
In this case, F(0)= ρ(y,x)g and F(1)=0.
The other change that must be applied is the macro properties calculation as follows:
u(y,x,0) =[f(y,x,1) - f(y,x,3) + f(y,x,5) - f(y,x,6) - f(y,x,7) + f(y,x,8) + 0.5F(0)] / rho(y,x)
u(y,x,1) =[f(y,x,2) - f(y,x,4) + f(y,x,5) + f(y,x,6) - f(y,x,7) - f(y,x,8) + 0.5F(1)] / rho(y,x)
Would you please tell me whether the above formulas are correct or not?
And the algorithm shown below is correct or not?
do tStep = 1,tMax
CALL wall_boundaries(f,image)
CALL compute_Macros(f,rho,u,uSqr)
CALL compute_Feq(fEq,rho,u,uSqr)
CALL collide(f,fEq,omega,image,rho,u)
CALL stream(f) ! with regard to periodic boundary conditions
new=sum(u(:,:,0))
t=(new-old)/old
print,tstep,abs(t)
if(abs(t).lt.eps) then
exit
end if
old=new
end do
I highly appreciate your considerations before,
Sincerely,
Narges