Calculating the density on velocity Dirichlet boundaries


I have a question about the way that the density is updated in VelocityDirichletBoundaryDynamics.

The function VelocityDirechletBoundaryDynamics::computeRhoBar() finally returns

T ​rhoBar =((T)2*rhoNormal+rhoOnWall-Descriptor<T>::SkordosFactor()*velNormal)
                     ​/ ((T)1+velNormal);

Considering the following case (copied for Kruger’s book),

it should give,
rho = [f0+f1+f3 + 2(f2+f5+f6)] / (1+velNormal)
where f0+f1+3 = rhoOnWall and f2+f5+f6 = rhoNormal, which leads to
rho = (2*rhoNormal + rhoOnWall) / (1+velNormal)

Based on that rhoBar = rho - 1.

rhoBar = (2*rhoNormal + rhoOnWall - 1 - velNormal) / (1+velNormal)
which is different with the one implemented in the code.
rhoBar = (2*rhoNormal + rhoOnWall - SkordosFactor * velNormal) / (1+velNormal)

So here I have two main questions,

  1. What is SkordosFactor()? It seems always to be one as set in src/latticeBoltzmann/roundOffPolicy.h.
  2. Where does the difference come from?

Thank you!


I happen to find out the answer to question 2. It is related to the difference between rho and rhoBar in Palabos. As explained in the link:

The “populations” in Palabos are defined as f_i^bar = f_i-t_i, in order to gain accuracy,

where t_i is the weight.

To calculate rho, one needs to calculate FF = [f0+f1+f3 + 2(f2+f5+f6)]. Let’s define FF_bar = f0^bar+f1^bar+f3^bar+2(f2^bar+f5^bar+f6^bar), so that FF = FF_bar + 1.

In Palabos, (T)2*rhoNormal+rhoOnWall equals FF_bar intead of FF. So, rho = (FF_bar + 1)/(1+velNormal), which leads to the returned value in computeRhoBar(), taking rhoBar = rho - 1:
rhoBar = (FF_bar - velNormal) / (1+velNormal).