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.
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)