I am trying to implement a really simple periodic flow but somehow I fail to get correct results. What is really intriguing me is that there isn’t much to go wrong so really have no clue of what is happening…
The flow is the one present in the paper of Zou et al. (1995) named “An improved incompressible Lattice Boltzmann Model for time-independent flows”, Journal of Statistical Physics Vol. 81, 35-48.
Briefly the flow has the following analytical solution:
the flow takes place in a periodic square box bounded between [0,2*pi].
To solve it I used the D2Q9 model with periodic boundary conditions in all 4 boundary planes. To produce the flow I just incorporate in the LB model the up mentioned body force. (When including the force in the BGK collision step I tested several body force models and results are more or less identical so the problem does not come from that).
In fact, the problem is that the solution that I obtain is not even close to the expected analytical one…
So I ask if anyone has ever tried to solve this simple flow?
Any help/suggestion would be very much appreciated…
have you re-computed the solutions for pressure and body force? I first had similar problems with the decaying Taylor-Green vortex flow which is quite similar (but time-dependent).
After you have made sure that you got the solutions right, you should check whether you have initialized the flow correctly. In order to do this, the equilibrium and at least the first order non-equilibrium should be initialized correctly.
How exactly does your solution look like? What happens?
Thank for you reply.
I do not know if I understand what you mean by recomputing the solutions for pressure and body force. What I do is to set the body force in the LB scheme and let it evolve. When my steady-state criterion is reached, i.e. (abs(u-u_old)+abs(v-v_old))/(abs(u_old+v_old))<1e-12, the LB simulation stops and then I compare the obtained macroscopic field with the analytical solution. The flow should look like two vortices in the middle with two half vortices at each side but my numerical solution yields two huge vortices in the middle and very very small vortices at the sides (and furthermore these are not symmetrical as they should be)…
I also implemented the same flow but I have removed the periodic BCs and changed them to velocity BCs, using Zou / He BCs with proper modifications to account for the external body force presence. These BCs have been used for several other flows and yield good results, however once more in this case, although the velocity solution looks better, the results are still not correct. What I find strange in this case is that after fixing the two velocity components the pressure at the boundary does not match that of the analytical solution. This gives a totally wrong pressure field. I believed that if this problem was fixed I would be able to use this flow as benchmark at least for velocity BCs case. Nonetheless, this is other problem…
Now, and since I am stuck for almost a week I think I will follow your last suggestion and properly initializing my field.
To be honest I always thought that for steady state the converged solution was pretty much independent of the initialization values of the distribution functions (as long as they wouldn’t threat stability). Therefore I only have initialized the equilibrium part.
you are right that for a steady flows the initialization is usually not important. So, there are two approaches:
start with some arbitrary configuration and wait until you have steady state
start with the expected solution and see whether anything changes
I think that both approaches are more or less equivalent, except for the facts that the first approach will require more time steps while the second approach will require more human setup time.
By the way, if you start with an arbitrary flow field at the beginning, its total momentum should be zero. Elsewise, your fluid will move at the end since you have no walls. This you do not want.
What I meant at the beginning of my first post: Make sure that the literature solution is really correct, i.e., plug the velocity field into the NS equations and see whether you can recover the body force correctly. Maybe, there is only a wrong factor. When you are absolutely sure that the analytic solution is correct, the error will be somewhere else.
Please let me know what the result of your analytic computation was.
Thanks for your fast reply. Yes, I have substituted the analytical solutions in Navier-Stokes and they are correct.
I will follow your hints and then let you know what I have obtained. In fact I want to thank you once more because your suggestion on starting with the analytical solution and then see what (negative) effects will the LB scheme have on it is pretty cleaver :-)…
Sorry for bothering you again. But the more I think about this problem the less I understand it
So as you suggested me I tried to start with the initial analytical solution and then apply the full LB scheme to ascertain how the method works.
To initialize the distribution functions field, f(i) I used the approach proposed in the paper of Mei et al. (2005) “Consistent initial conditions for the lattice Boltzmann simulations”.
If I understood it correctly their procedure is pretty simple: In order to initialize our f(i) field one fixes the velocity to the analytical solution and uses exactly the same LB algorithm as to find rho and u but this time u is fixed and only rho is computed from the zeroth order moment of f(i). Am I thinking well? At least using this reasoning I was able to successfully initialize the rho field of a Couette and a Poiseuille flow. However when applying it to the periodic flow my results are still awkward…
So in conclusion, I want to obtain a steady-state flow solution through LB. I use two approaches: (1) starting from an arbitrary initial solution or (2) trying to compute the distribution function that provide me proper initialization and then applying the LB algorithm. In both cases my results are totally messed up.
Since the flow solution I intend to find is correct, i.e. it respects continuity and Navier-Stokes eqs. my problem also does not come from here.
I am pretty much in an dead end
Since the initialization of the Taylor-Green vortex flow is identical to the initialization of the density in my problem (of periodic flow) I ask you if you don’t mind to give some further hints or the send me the algorithm of your code, I would be very very grateful… and that would surely help me understand what is happening!
Really helpful post as I tried to do the same thing a time ago but with Taylor-Green vortex. I experienced absolutely the same picture. But right now probably there is a bit better understanding:
Let us start with the Taylor-Green vortex decay. So if you initialize the velocity field as it is given by analytical equation you would obtain the decay with factor -2nut. You can check it. What does it mean? The equation to be solved is the Navier-Stokes equation:
\partial_t \rho u_{\alpha}+\partial_{\beta} \rho u_{\alpha} u_{\beta}=-\partial_{\alpha} p +\partial_{\beta}(\partial_{\beta} u_{\alpha}+\partial_{\alpha} u_{\beta})
So in reality If you initialize the system with the initial velocity as:
u(x,y)=U0sin(x)sin(y)
v(x,y)=U0cos(x)cos(y)
p(x,y)=1/4U0^2(cos(2x)-cos(2y))
You will obtain the Taylor-Green vortex decay if you do not use any forcing: Taylor-Green solution
If you have forcing the Navier-Stokes equation is different:
\partial_t \rho u_{\alpha}+\partial_{\beta} \rho u_{\alpha} u_{\beta}=-\partial_{\alpha} p +\partial_{\beta}(\partial_{\beta} u_{\alpha}+\partial_{\alpha} u_{\beta})+F_{\alpha}
What we want to do is to correct the pressure term with forcing. That means the term which was previously produced with the pressure gradient (and is analytical solution) is substituted with forcing. However, the pressure gradient exists in the LBM scheme that’s why it is told that you need to work in the incompressible limit, when the forcing is bigger than the gradient. The equation should become:
\partial_t \rho u_{\alpha}+\partial_{\beta} \rho u_{\alpha} u_{\beta}=F_{\alpha} +\partial_{\beta}(\partial_{\beta} u_{\alpha}+\partial_{\alpha} u_{\beta})
What happens in your case - you want to have the constant flow - that means you there are no time-dependent variables:
\partial_{\beta} \rho u_{\alpha} u_{\beta}=F_{\alpha} +\partial_{\beta}(\partial_{\beta} u_{\alpha}+\partial_{\alpha} u_{\beta})
How to achieve it? I guess the good thing to do is to perform only ONE step without forcing as it’s done in the Taylor-Green case. You know the pressure initialized by analytical formula. Take the -gradient of the analytical formula:
fx=-grad p
fy = -grad p
or use the equations which are given for forcing and which you provided them.
Try to substitute forcing which will give the same result as you expect for decay in one step(!) - you initialize the density with constant in order not to have density gradients. If you achieved it - that means you know how to control forcing to achieve the same pressure gradient provided in the analytical equation. Then you can move to the next step as controlling forcing to eliminate pressure gradient - then you will obtain constant solution
Big thing is the initialization - as soon as you introduced the forcing you need to do a half shifting velocity. Which also put the headache in terms of how to initialize. Once I tried to obtain pressure waves - so I initialize the system with the equilibrium values, but velocity was minus half shifted with forcing in order to obtain the necessary macroscopic velocity which is plus half force shifted.
Hopefully it will help a bit, I don’t have time to do it by myself right now but I am extremely interested in your results. If you can achieve the required picture could you please send me the code to shurik.kuzmin@gmail.com
@Alex
Thanks for your reply. A few weeks ago, I had a similar experience. You either take the pressure or the forcing. The interesting point is that you can also take both, as long as $- \nabla p + \vec f$ is correct.
@Goncalo
I suggest the following: Absorb the entire term $- \nabla p + \vec f$ into the pressure, i.e., $\vec f = 0$. In this case, you do not have to think about forcing (if a pressure p exists which gives the right gradient). I have also sent you an email with an article draft which may be relevant for your problem.
