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]
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.