Dear all,
I am experimenting with a simple stationary Poiseuille flow created from a D2Q9 LBGK model, where flow is induced by an external force along the x-Axis with periodic boundary conditions in the same direction and bounce-back boundaries at top and bottom (Actually my setup is the same as the channel flow in He,Zou,Luo,Dembo 1997). However, I also tried to evaluate the stress tensor recovered from the non-equilibrium part of the pdfs. There I am getting non-zero entries on the diagonal of the stress tensor.
E.g., at a center velocity of 0.00228 and an error in the velocity profile below 1e-16 (everywhere!), I have the following the following profile for the s[sub]xx[/sub] entries of my stress tensor along the y-axis (channel hight is exactly 6 nodes):
-1.3e-6 2.2e-7 1.0e-6 1.0e-6 2.2e-7 -1.3e-6.
Computing the stress tensor by central differences from the velocity field results in a trace free stress tensor - which is what i really want as of course the x-velocity is constant along x, and the y-velocity is constant along y. Is this then a bug in my implementation or just another LBGK / Bounce Back artifact? I have used a relaxation time with vanishing slip velocities at boundaries for this case!
have you checked if the implementation for the stress evaluation is correct? Please try a larger channel and plot the stress profile (both xx- and xy-components).
The first entry of the stress tensor should be zero - but this does not seem to be the case.
Could this be caused by the body force term?
Could one easily try this with OpenLB or Palabos (I mean setting up a 2d poiseuille flow and evaluating the stress tensor).
Okay, so here is how I compute the stress tensor:
I take the second order moment of my populations field:
S := Sum[ c[sub]i[/sub][sup]2[/sup] (f[sub]i[/sub] - f[sub]i[/sub][sup]eq[/sup]) ],
where c[sub]i[/sub] are the lattice velocities and f[sub]i[/sub] are the population densities at each site. I subtract the equilibrium part (which corresponds to the normal stresses and momentum flux tensor) f[sub]i[/sub][sup]eq[/sup].
I would expect that, in steady state, the xx-component should indeed vanish. Since your flow is gravity driven along the x-axis, there should be a force correction term in the stress tensor evaluation. This may cause the deviations.
Do you use a force correction when you compute the velocity from the populations?
Hi!
Good point - but this should not make a difference here. I guess if that was the mistake, I would also have deviations in the velocity profile - but here I have only a truncation error!
-Simon
The setup is the same as the channel flow in the paper (He,Zou,Luo,Dembo 1997) I was mentioning. The force term is added after the BGK relaxation. See equations (1) and (5), respectively.
This may well just be a problem with bounce back. We know it gets the velocity wrong at the wall, unless you use an MRT/TRT scheme (or rather, it gives an error of second order) so it could also get the stress wrong. It does sound a little strange that you predict non-zero diagonal components of the stress tensor in a simple channel flow where the velocity field is one dimensional, but then again the bounce back method doesn’t say anything about the stress field.
It is certainly worth trying a different forcing scheme. The basic method, as you quote from the He at al paper, is not consistent with the Chapman-Enskog expansion. In fact, it alters the momentum flux tensor, which gives the viscous stress. To stop this, a body force R_i should have the following moments:
\sum R_i=0
\sum R_ic_i=F
\sum R_ic_i*c_i=uF+Fu (where u and F are both vectors).
In Poiseuille flow the basic formulation should be sufficient because u_x is a function of y only and u_y=0. However, any variation could, in theory alter the stress. Note that if the third moment above is violated you do not recover the Navier-Stokes equations (at least not to the desired order of accuracy). You can use the explicit form for R_i given in Luo PRL 1998 (title has something about non-ideal gasses in it). If that doesn’t sort things out, and you’re convinced your code is bug free, then I’d be inclined to say that the bounce back gives some higher order (in Knundsen number, ie Burnet-type) contribution (or error) to the stress field.
ps
If you are using a small tau (or have a high reynolds number) you may also want to try using an MRT scheme because BGK collisions over-relax f_i (and this is what causes a lot of the instabilities). This gives rise to oscillations in the higher non-conserved moments, which should die out over longer timescales, but perhaps they are still present? (you could, for example, plot the non-equilibrium moments over time and see what they look like) Remember also that the stress tensor comes from moments of f_i^(1), so getting your hands on it with f_i^{non-eq)=epsilonf_i^(1)+epsilon^2f_i^(2)+… could give you additional (non-Navier-Stokes) contributions.
Thanks for sharing your thoughts! I tried the improved force term fomr Luo 1998 but the results were the same in this test case.
Would be interesting to know whether this is some BGK-oddity. Since I have tried two different matlab-implementations and both showed the same error, I am pretty confident, that it is not a programming issue.
Maybe I will try Zou-He instead of Bounce Back or even MRT in the next few weeks and then post the results…
Before doing that, have a look again at those non-zero numbers for the stress. I wouldn’t be surprised if they are proportional to the shear stress (the constant of proportionality being something like tau^2 ish)? If so, this is probably a Burnett-type effect.
interesting discussion. I am not sure if I remember correctly: the last time I checked, the diagonal elements in steady Poiseuille flow (BB at walls, periodicty along flow direction, Guo forcing) vanished (up to machine precision). If I find some time, I will run a small simulation with my code and let you know.