Due to my previous question, where I asked about flux calculation in lattice Boltzmann (LB) method here, I have more or less same question for deviatoric stress tensor calculation due to pseudo-compressibility of LB method. In fact, strain rate tensor is calculated in LB by using non-equilibrium part of distribution functions ($f_{i}^\text{neq}$) as:
$$\hat{\varepsilon}_{\alpha\beta} = -\frac{1}{2\hat{\tau} \hat{\rho} \hat{c}_{s}^{2}} \sum_{i} f_{i}^\text{neq} c_{i\alpha}c_{i\beta}$$
Where hat quantities are dimensionless quantities, $\hat{\tau}$ is dimensionless relaxation time ($\hat{\tau} > 0.5$), $\hat{\rho}$ is dimensionless instantaneous density of the fluid, $\hat{c}_{s}^{2} = \frac{1}{3}$, and $c_{i\alpha}$ is the $i$th discrete velocity in $\alpha$ direction. Deviatoric stress is defined as:
$$\hat{\sigma}_{\alpha\beta} = 2 \hat{\mu} \hat{\varepsilon}_{\alpha\beta}$$
But, we have for dimensionless viscosity: $\hat{\mu} = \hat{c}_{s}^{2} \left(1-\frac{1}{2\hat{\tau}}\right) \hat{\tau}$
So, finally:
$$\hat{\sigma}_{\alpha\beta} = -\left(1-\frac{1}{2\hat{\tau}}\right) \frac{1}{\hat{\rho}} \sum_{i} f_{i}^\text{neq} c_{i\alpha}c_{i\beta}$$
This is in contrast with what usually LB people use, such as this one from Krüger et. al. as:
$$\hat{\sigma}_{\alpha\beta} = -\left(1-\frac{1}{2\hat{\tau}}\right) \sum_{i} f_{i}^\text{neq} c_{i\alpha} c_{i\beta}$$
I understand at incompressible limit ($Mach \rightarrow 0$), $\hat{\rho}$ should be close to 1, but for my simulations where $Mach \sim 0.06$ and $Re \sim 600$, $\hat{\rho}$ may fluctuate quite a bit. So my question is which of these formulas should be used to calculate deviatoric stress? Should I assume $\hat{\rho} \sim 1$ even with my high $Mach$ number? Any suggestion is truly appreciated.
