I’m trying to compute the nondimensionalized force (Cx and Cy) around a cylinder section using the D2Q9 sheme. I know that this topic has been seen over and over again, but i’m facing an issue in which I didn’t find any information.
I’m actually wondering why the values of Cx or Cy depend on the value of the relaxation time? This is obviously wrong because Cx and Cy only depend on the Reynolds number (for a given cylinder shape), and for a given Reynolds number you can have infinite value of omega.
Here is what I do. I put the force calculus step after the collide step and before the stream step, and I use the momentum exchang method:
Where:
S is the stress, v is the lattice velocity, and f is the distribution function (of course).
(x,y) define the position of contour lattice
(xf,yf) define the position of a neighbouring fluid lattice.
j=opposite(i)
Ok, I’m not sure I understand but here are some thoughts. The relaxation time sets the Reynolds number, so effects the force. If changing the relaxation time while keeping the Reynolds number fixed give you different values for the force, then there are a number of possible explanations.
If you alter the relaxation time for a fixed Re then you must either change the number of grid points in a characteristic length (ie in radius of the cylinder in lattice units) or change the characteristic velocity. This means you are changing the resolution of the system, ie dx and dt. And also the Mach number. This is going to influence how accurate your predictions are. If you want to check convergence, you should keep the Mach number (or the velocity in lattice units) and Reynolds number constant, and adjust only dx (or the grid size). Each mesh (ie different numbers of grid points) will give different results, but all have the same Re. Hopefully the predictions will converge to the “true” solution with second order accuracy.
Most lattice boundary Boltzmann boundary conditions give an additional error term, often call a “numerical slip”. This error depends of the relaxation time. In general it is larger for lower Reynolds number.
Also, the higher order moments and errors in the LBE are also controlled by the relaxation times, so adjusting them can effect the accuracy of a simulation.
Actually, you really well understood what was my problem and you deeply shed light on it!
In my case, I’d rather remain thrifty in terms of system resolution, so that arrays dimensions don’t grow too fast. I opt then for a new computation of omega each new Reynolds number, and an adjustment of the resolution if any. Everything is based on the following relation:
Ma = Re/(dsqrt(3)) * (tau-0.5)
omega = 1/tau = Re/(sqrt(3)dMa+0.5Re)
Where d is the number of lattices along the cylinder diameter, omega must remain less than 1.8 and Ma=0.2. Do you think this strategy fits well this kind of problem?
I can’t really check the convergence in d for now because my PC is not powerfull enough, but results seem so much better! I’ll try this tomorrow at work.
Otherwise, I’m looking for documentation on numerical slip phenomenon, do you have any reference? Does an acurate wall treatment (such a second order Yu-Mei-Luo-Shyy) will decrease this slip?
Many thanks again for your answer,
JB
PS: I’m not ready yet for multi-relaxation-time methods
What you are doing sounds reasonable. Like any numerical method, LBE is an approximation. The error in LBE (for Navier–Stokes) is, in theory, dx^2 (grid spacing squared). So, if you double the mesh resolution there error should be 4 times smaller. That said, there are often other errors which can reduce the algorithm to first order. The numerical slip error should be second order for “half way” bounce back, but many implementations are only first order. If you are interested in such things, a good paper to start with is He, Zhou, Luo and Dembo, J. Stat. Phys,. 1997. For curved boundaries, I’m not sure what the best method is.
