Drag and Lift Forces: Momentum exchange method

Lattice Boltzmann Method LBM: Theory
#1

Hi all,

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
#2

Hi 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:

[code=“fortran”]
Xforce = Xforce + 2.d0 * CVEL(1,ii)
& *(ftemp(i,x,y) - ftemp(ii,ix,iy))
Yforce = Yforce + 2.d0 * CVEL(2,ii)
& *(ftemp(i,x,y) - ftemp(ii,ix,iy))



I hope this helps!

Regards,

Michael
#3

Thank You Michael,

I will try your suggestion. Will let you know the results.

According to the article “Force Evaluation in the Lattice Boltzmann Method Involving Curved Geometry” by Mei, Renwei et al. (2.3)

Summation of all the contribution over all boundary nodes xb, belonging to the body, the total force (by the solid body on the fluid) is:

F = sum(all xb) : sum(for non-zero lattice) e_i( f_i (xb) + f_opp(i) (x_b + e_opp(i))

I am following this equation.

In my previous code, there was a bug that I removed and saw some results, however, not accurate one, not 100% success.

Regards
Ali

#4

Hi Ali,

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:


F = 2 * f_j (x_b + e_i) * e_i = 2 * f_j (x_f) * e_i

and therefore you could obtain the force by changing the lines for Xforce and Yforce to:

[code=“fortran”]
Xforce = Xforce + 2.d0 * CVEL(1,ii) * ftemp(ii,ix,iy)
Yforce = Yforce + 2.d0 * CVEL(2,ii) * ftemp(ii,ix,iy)



(If your distribution functions are zero for the boundary nodes, you could still stick to the changes I suggested earlier.)

Regards,

Michael
#5

Thanks Dear Michael,

I want to share my code with you if you can give me your email.

Best Regards
Ali

#6

Hello,

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:

Xforce = Xforce + 2.d0 * CVEL(1,ii) *(ftemp(i,xi,yi) - ftemp(ii,xi,yi))
Yforce = Yforce + 2.d0 * CVEL(2,ii) *(ftemp(i,xi,yi) - ftemp(ii,xi,yi))

I am wondering if I am doing it in a correct way.
Thanks