[AdvancedSimulation]

Dear colleagues,


I am relatively new to LBM and I have tried for more than a week to compute the lift and drag coeff in the 2D problem of a Re=100 round cylinder in a channel with a Poiseuille flow (to be validated with Schafer & Turek https://link.springer.com/chapter/10...322-89849-4_39). The reported benchmark values in Schafer & Turek are periodic oscillations witch Cd_max = 3.22 with Cd_avg =3.20 approx; Cl=+/- 0.96 approx.


For the calculation I have used LBM with SRT and the momentum exchange method with the easy-to-code bounce-back method (for cylinder and walls). I am aware the method detects my cylinder boundaries properly and does not include the channel walls in the calculation. I attach a dynamic gif which shows how each solid node is properly identified from each node of the boundary (from the i-th yellow node position, neighbouring solid nodes are plotted sequentially).



The problem is that every time I change the number of lattices in x or y direction, the lift and drag plots vary (especially in terms of the period of the oscillation, which i.e. seems to be kinda doubled if number of y-lattices is doubled, as in the attached figure). The flow velocity seems to be correct (changes may be due to the different frequency in the oscillations that change the "time" at which velocity profile is being visualised in my figure).



My deltax=deltay=deltat=1, and the proportion of the domain is kept (image attached) as so I do not understand why my Cl and Cd have such changes.


May I have a bug in my code? May I be missing anything?

May I be missing anything particularly in the post-processing? Perhaps I have to make changes in the coefficient since Fx and Fy increase with the number of lattices (my guess is no because D also increases)?



Possible relevant facts:
ly (number of nodes in y-direction) determines the channel height (H=ly-2, since the top and bottom walls nodes must be not counted). Since H/D=4.1, I compute the diameter from ly as D=H/4.1 = (ly - 2)/4.1. Thus, dimensionalisation depends on ly.

I compute force coefficients by (note U=U_mean at inlet and D is the one calculated above in the beginning of the code):
Fx./(0.5.*(U^2).*ceil(D))
Fy./(0.5.*(U^2).*ceil(D))



Any support is pretty much appreciated!!!


Many thanks in advance