I am simulating a square cylinder in LBM and want to calculate lift and drag by Momentum exchange method. I found a subroutine in this forum and try to implement but unsuccessful (zero values for FX and FY reported). Here is my Subroutine, kindly have a look and suggest if anything wrong:

[code=“fortran”]
FX = 0.d0
FY = 0.d0

DO x = 1, Xmax
DO y = 1, Ymax
Xforce = 0.d0
Yforce = 0.d0
IF (obst(x,y)) THEN !!! obst(x,y): is the logical array having the boundary points of cylinder
DO i = 1, 8
ii = ANTIALPH(i)
ix = x + CVEL(1,ii)
iy = y + CVEL(2,ii)
IF (.not. obst(ix,iy)) THEN
Xforce = Xforce + CVEL(1,ii)
& *(ftemp(i,x,y) + ftemp(ii,ix,iy))
Yforce = Yforce + CVEL(2,ii)
& *(ftemp(i,x,y) + ftemp(ii,ix,iy))
END IF
END DO
END IF
FX = FX + Xforce
FY = FY + Yforce
END DO
END DO

The time iteration works like that;
1- Inlet_Outlet_BCs
2- Streaming
3- Bounceback for obstacle
4- Collision (BGK)
5- ForceComputation
6- Data Output
7- Back to 1
Is the sequence is right or need to modify?
Thank you for any help.
Best Regards
Ali

I think you need to calculate the difference between the two distribution functions across the boundary rather than the sum. According to Tony Ladd’s papers on solid/fluid interactions, the force acting on a boundary site for a stationary object due to link i is equal to:

F = 2 * (f_i (x_b) - f_j (x_b + e_i)) * e_i

where j is the conjugate link to i.

To correct your code, you will need to change the two lines calculating Xforce and Yforce to:

I’ve subsequently taken a closer look at the papers by both Ladd and Mei et al., and realised there is another modification you can make to your original code for the force on the object.

Ladd makes the assumption that the solid object is filled with fluid and hence there is a direct exchange of momentum across the solid/fluid boundary. I believe you can still use the same formula for the force if the boundary is entirely solid, but the value of the distribution function for the boundary (x_b) would be zero.

Meanwhile, Mei et al. state that the post-collisional distribution function for the boundary node (x_b) needs to be determined from the value for the neighbouring fluid node (x_f): this is given by equation 7 in their paper. If the object is kept stationary and the boundary is halfway between the boundary and fluid nodes (as per usual for bounce-back), the post-collisional distribution function at x_b will be equal to the value at x_f (as the equilibrium term would be cancelled out).

In both cases, this would lead to a force at the boundary equal to:

I will appreciate it if you could check out my question below, I stuck in it for a couple of weeks but the result is not satisfactory. My distribution functions are zero for the solid nodes, so for each solid node I check the neighbor nodes if it is a fluid node I am using the following formula: