Implementation of the Momentum Exchange (ME) method

Lattice Boltzmann Method LBM: Algorithms
#1

Hello,

I am a novice programmer and trying to use the LBM on the classic circular cylinder unsteady problem. I have written a code in the Julia language using a D2Q9 lattice. I have the code working and see the alternate vortex shedding phenomena one would expect at low Reynolds numbers. I also have managed (with much difficulty!) to get the Momentum Exchange method working to estimate my Cl and Cd on the cylinder. I am getting close to the expected results but have the following question:

I am summing all the forces in the x and y (separately) by identifying the fluid nodes directly adjacent to a solid node and on said nodes taking the momentum exchange in thee directions (three bordering points per node, one horizontal or vertical and two diagonals). I can see that my code is picking up any node adjacent, but I think there is a small force missing in a node not adjacent but diagonal. These diagonal nodes will have 1 and only 1 link from fluid to solid and as such, I would assume are imparting (or being imparted) a small force – perhaps this is only a very small amount but I want to know if I am thinking correctly and if I need to modify my code?

I have included a diagram and an excerpt of my code to show where the issue may be.

Any reviews would be greatly appreciated.

Regards,

Jordan.

http://imgur.com/a/JRWRO

FORCE ESTIMATION STEP

for i = 2:(Nx - 1)
  for j = 2:(Ny - 1)
    for k = 1:9
      if (Cylinder[i, j] == 0) && (Cylinder[i + Cx[k], j] == 1) # HORIZONTAL AXIS
        Momentum_Parameter_x[i, j, k] =  Cx[k]*(f_OUT[i, j, k] + f_IN[i + Cx[k], j, Opp[k]])
      end
      if (Cylinder[i, j] == 0) && (Cylinder[i, j + Cy[k]] == 1) # VERTICAL AXIS
        Momentum_Parameter_y[i, j, k] =  Cy[k]*(f_OUT[i, j, k] + f_IN[i, j + Cy[k], Opp[k]])
      end
      #if (k > 5) && (Cylinder[i, j] == 0) && (Cylinder[i + Cx[k], j + Cy[k]] == 1) # DIAGONAL???
        #Momentum_Parameter_xy[i, j, k] =  Cy[k]*(f_OUT[i, j, k] + f_IN[i + Cx[k], j + Cy[k], Opp[k]])
      #end
    end
  end
end

Fx_Cylinder[Cycle] = sum(Momentum_Parameter_x)
Fy_Cylinder[Cycle] = sum(Momentum_Parameter_y)

#2

Hi All,

I have managed to find other info on the above and ironed out some of the issues.

I now have predictions of the Coefficient of Lift and Coefficient of Drag with good precision to the experimental results of that particular Reynolds Number. What is interesting for me is that I am getting about one order of magnitude more ‘noise’ than with other LBM predictions. My Cd values for example are oscillating a bit too much (even before instabilities occur in a circular or square cylinder example).

Is there common fixes for this in LBM? Any ideas on where to start looking in my method / code? I have done tests at various Reynolds Numbers, different shapes, with and without boundary walls, all have the same noise patterns.

For reference my code is a simple SRT model, bounce back BCs and Zou He inflow/outflow BCs.

Regards,

Jordan.