Alright.

Today I have solved the problem, and I would like to explain it in some detail here so that other people can use my enlightenment to save a lot of time.

First, about the correct definition of the deviatoric shear stress in physical units and LB units:

I am using Jonas’ PhD thesis as reference. In eq. (2.66) Jonas writes that

```
Pi^(1) = - 2 * c_s^2 * tau * rho * S
```

where c_s is the sound speed, tau is the relaxation parameter, rho is the density and S is the strain rate tensor (1.6),

```
S = 1 / 2 * (nabla u + (nabla u)^T)
```

On the other hand, we know that the physical deviatoric shear stress tensor is (1.5)

```
sigma = 2 * rho * nu * S
```

where nu is the kinematic viscosity. Moreover the total momentum tensor is (1.3)

```
Pi = rho * u * u + p * 1 - sigma
```

The tensor Pi^(1) can be calculated using the LBM (3.23),

```
Pi^(1) = sum_i f_i * c_i * c_i - rho * u * u - c_s^2 * 1
```

Now, when I compare everything, then I find the important result

```
sigma_ij = - (1 - 1 / (2 * tau)) * [sum_k f_k * c_ki * c_kj - rho * u_i * u_j - c_s^2 * delta_ij]
```

since

```
nu = c_s^2 (tau - 1 / 2)
```

This means that I need an additional factor of (1 - 1 / (2 * tau)) and a minus sign for the calculation of the shear stress tensor sigma.

Second, I have wasted a lot of time, because I was really stupid. I tried to implement this factor using the code

```
stress *= (1 - 1 / 2 / tau);
```

As you know, 1 and 2 are integers, which means that c++ thinks that 1 / 2 is also an integer, i. e. 0 in this case. 0 divided by tau is always zero, so this factor is always 1. Now I have changed the code to

```
stress *= (1 - 1 / tau / 2);
```

and the simulated results for the shear stress match the analytical results. So much about explaining a computer how to compute numbers.

Timm