Numerical stability


and another question:
I have simulated the velocity and shear stress profiles in a long pipe using different values of the relaxation parameter tau. Now I have observed a curious thing: For 0.51 < tau < 0.65 the velocity profile approximately matches the theoretical expectations, but the shear stress profile becomes more and more “coarse” the smaller tau is. This means that the code produces inexact values for the shear stress, but good values for the velocity.
I compute the shear stress using

void computeDevStress()
for(k = k1; k < k2; k++) {
  if(obstacle[k] == 0) {
    devStressXY[k] = 0;
    devStressXZ[k] = 0;
    devStressYZ[k] = 0;
    devStressXX[k] = 0;
    devStressYY[k] = 0;
    devStressZZ[k] = 0; } }

for(k = k1; k < k2; k++) {
  if(obstacle[k] == 0) {

    for(ii = 0; ii < NRVELO; ii++) {
      devStressXY[k] -= (1 - 1 / tau / 2) * (ff[ii][k] - fe[ii]) * vl[ii][0] * vl[ii][1];
      devStressXZ[k] -= (1 - 1 / tau / 2) * (ff[ii][k] - fe[ii]) * vl[ii][0] * vl[ii][2];
      devStressYZ[k] -= (1 - 1 / tau / 2) * (ff[ii][k] - fe[ii]) * vl[ii][1] * vl[ii][2];
      devStressXX[k] -= (1 - 1 / tau / 2) * (ff[ii][k] - fe[ii]) * (SQ(vl[ii][0]) * 2. / 3. - (SQ(vl[ii][1]) + SQ(vl[ii][2])) / 3.);
      devStressYY[k] -= (1 - 1 / tau / 2) * (ff[ii][k] - fe[ii]) * (SQ(vl[ii][1]) * 2. / 3. - (SQ(vl[ii][0]) + SQ(vl[ii][2])) / 3.);
      devStressZZ[k] -= (1 - 1 / tau / 2) * (ff[ii][k] - fe[ii]) * (SQ(vl[ii][2]) * 2. / 3. - (SQ(vl[ii][0]) + SQ(vl[ii][1])) / 3.); } } }

The LB code is a D3Q15, so it cannot handle values of tau very close to 0.5. But like I said: While the velocity still is very accurate, the shear stress is not. Is it possible that this effect originates in the fact that the velocity is the first and the shear stress the second moment of the distribution functions?
For now, I just make sure that tau is not too small, but is it possible to correct this problem generally (except changing the lattice structure)?

Thank you,

It may be that what you observe is a limitation of the numerical model, but it may also be that your setup violates a physical principle, in which case it is not astonishing that the numerics disagrees with the analytical prediction.

An important thing to clarify is how you approach the limit tau -> 1/2. Do you adjust the grid resolution (L)? And do you adjust the velocity (U), measured in lattice units? If both of them are kept constant, then the Reynolds number increases as tau approaches 1/2:

Re = U L / nu, where nu= 1/3(tau-1/2).

At large Reynolds number, you need to be aware of the possibility that the simulation enters a turbulent regime (although the laminar solution still exists, it becomes unstable). In this regime, the lattice resolution must be increased to embrace the size of small-scale eddies.

If this isn’t what you did already, it would be interesting to decrease tau while keeping the Reynolds number constant (thus decreasing the velocity U), and see if you still observe a decreasing accuracy on fitting the stress tensor Pi with velocity derivatives.

I want to see the effect of the Reynolds number and the lattice resolution on my accuracy.
For that purpose I changed the Reynolds number (between 2 and 100) and the lattice resolution (relative factors 1, 2, 4). My first observation of course is that tau decreases, if Re is increased or if the lattice resolution is decreased (lattice constant increased). This is clear.
If tau is too close to 1/2 I have the possibility to increase the Mach number (at constant Re). Now it seems that the effect I have described in my initial post only happens for small values of tau, independent of Re, Ma and the lattice resolution. This would indicate a problem with the numerical model. But I should check this in more detail.
If the coarsing of my solution depends on Re at fixed tau, then it seems to be a problem due to the lattice resolution. Am I right?
Is the accuracy of the second moment (shear stress) of f_i in principle always worse than of the first moment (velocity)?

I think that if you want to decide for sure whether your simulation shows an expected behavior or not, you should vary only one parameter at a time and observe how the error of a given quantity evolves. There are four variables Re, U, L, and tau which are linked through the equation in my first post. You may therefore choose two variables which you will fix (keep at the same value from one experiment to another) and one variable which you will vary. The fourth variable is free and varies according to the above-mentioned relation.

Also, I am not sure what you mean by a coarsening of the solution. A good value to compute for a quantitative analysis is the L2-norm error between numerical and analytical solution.

If the coarsing of my solution depends on Re at fixed tau, then it seems to be a problem due to the lattice resolution

Do you vary Re by varying U or by varying L? If you increase U, you increase the Mach number and your error may be related to compressibility effects.

Is the accuracy of the second moment (shear stress) of f_i in principle always worse than of the first moment (velocity)?

That’s a good question, and I would need to think more carefully to be sure to give the right answer. Somebody else may have an opinion on this? Either way, the error on the shear stress should decrease as you reach the hydrodynamic limit. It may decrease at a first-order or second-order rate, depending on its accuracy, but it should decrease.

For BGK, the incompressible Navier-Stokes solution is known to be reached in the diffusive limit, which you obtain as follows: keep Re and nu constant, then start decreasing U and adjusting the lattice resolution (L) accordingly. The L2-norm error on the velocity should then decrease at a second-order rate, i.e. err(L) \propto 1/L^2 (if you display err(L) on a log-log plot, you should observe a straight line with slope -2). In the same way, the error on Pi vs. velocity derivatives should decrease with 1/L. You can easily observe numerically (and I would be interested in knowing the answer) whether it decreases at first-order (slope -1) or second-order (slope -2).

Hello Jonas,

I have some results now, and they are strange.
I have tested the same geometry (channel with length L, width W and height H and aspect ratio L:W:H = 8:1:2) for fixed Reynolds number (Re = 2) and three different lattice resolutions (ratio 1:2:4). Moreover I have used two different boundary conditions:

  1. periodic boundary conditions with homogeneous and constant body force
  2. velocity boundary conditions without body force
    So I have two different data sets. The relaxation parameters tau are
    0.653205 (80x10x20 lattice)
    0.84641 (160x20x40 lattice)
    1.19282 (320x40x80 lattice)
    The L2-norm errors for my velocity u_x and sigma_xy, sigma_xz (flow in x-direction) are

u_x: 0.03561 and 0.004811 and 0.002164
sigma_xy: 0.04238 and 0.002546 and 0.008582
sigma_xz: 0.08343 and 0.004798 and 0.01719

for the periodic boundary conditions and

u_x: 0.01565 and 0.001226 and 0.001313
sigma_xy: 0.04701 and 0.002503 and 0.009197
sigma_xz: 0.08690 and 0.003501 and 0.01799

for the velocity driven boundary conditions.
What do we see? At first, both boundary conditions approximately result in the same errors. Fine.
But increasing the resolution does not necessarily decrease the errors, especially when dealing with the shear stresses.
Does anybody know, what happens here?


you are varying you relaxation time in order to keep the Re number constant when increasing the resolution. This means that the reference velocity remains constant. Therefore you will have a compressibility error that will appear. (You should vary U and L instead of nu and L, see the thesis of Jonas for more explanations p. 28.) The error should therefore be computed as

E=sqrt(1/N*sum_i (u_i-u_{ia})^2/u_{ref})

where N is the total number of nodes, u_i the computed value of the velocity at node i, u_{ia} its analytical value, and u_{ref} the reference value of the velocity used for your Re computation. Therefore you compute the error per node evolution.

Does it change the results? It would be indeed interesting to see what is the decrease rate of the error of sigma.

Hi Orestis,

thank you for your additional advise. At first one should note that

\tau = \sqrt{3} L \frac{Ma}{Re} + 1 / 2

with $Re = L u / \nu$ and $Ma = u / c_s$. Since tau and Re shall be fixed, I have to make sure that the product L Ma is constant. This is of course a different approach than mine (keep Re and Ma constant) and I see that the error behaves different in this case. I will tell you when I am done.


Ok, I have performed three simulations with three different lattice resolutions, but constant Re and tau. And you are of course right. The accuracy of the velocity scales with a power of 2. And the accuracy of the shear stress scales with a power of approximately 1.5! But I have to perform more simulations after the weekend.
Now, what do you think of a factor of 1.5?

Have a nice weekend,

Back again…
After five simulations I have found the result that the L2 norm of the velocity obeys

\epsilon_u  \propto Ma^2.039

This is a nice result, because this is expected. For the two shear stress components however I have found fit values of 1.523 and 1.477, which is very close to 1.5. If I understand you correctly, this is not really the expectation, is it?