Forcing is handled incorrectly when using a compressible equilibrium

Forcing terms are currently used incorrectly when using the compressible form of the PDF's equilibrium in lbmpy.

In LBM, we want the user to specify the variable force in terms of a body force density vector

\boldsymbol{F}

. In the compressible case, is should be formally

\boldsymbol{F}(\boldsymbol{x},t) = \rho(\boldsymbol{x},t) \boldsymbol{a}

, where

\rho(\boldsymbol{x},t)

is the density and

\boldsymbol{a}

is an acceleration. I assume that

\boldsymbol{a}

is constant in the whole domain, as I could not find a counter example where this is not the case. In contrast, I assume

\rho(\boldsymbol{x},t)

to depend on the position

\boldsymbol{x}

and time

t

because that's what one observes in compressible fluids.

In lbmpy, such a compressible case is currently treated as follows:

1. A user can only specify the body force density
   \boldsymbol{F} = \text{const.}
   , which is a single vector and independent of
   \boldsymbol{x}
   or
   t
   . Let us assume that
   \rho(\boldsymbol{x},0) = 1
   . During the simulation,
   \rho(\boldsymbol{x},t)
   will eventually change because of the flow dynamics.
2. In conservedquantitycomputation.py, line 240, the macroscopic velocity for constructing the equilibrium, is shifted by
   \boldsymbol{F}/2/\rho(\boldsymbol{x},t)
   . The density is indeed taken dependent on
   \boldsymbol{x}
   and
   t
   here. Since the user has specified a force density independent of
   \boldsymbol{x}
   and
   t
   , lbmpy essentially assumes that the acceleration
   \boldsymbol{a}
   varies with
   \boldsymbol{x}
   and
   t
   : We have
   \boldsymbol{F} = \text{const.} = \rho(\boldsymbol{x},0) \boldsymbol{a}
   . Therefore, if
   \boldsymbol{F} / \rho(\boldsymbol{x},t)
   with
   \rho(\boldsymbol{x},0) \neq \rho(\boldsymbol{x},t)
   , the acceleration
   \boldsymbol{a}
   must change in time
   t
   and likewise, it could also change in
   \boldsymbol{x}
   , so that
   \boldsymbol{a}(\boldsymbol{x}, t) \neq \text{const.}
3. Then, in momentbasedmethod.py, line 397, lbmpy directly uses
   \boldsymbol{F}
   as specified by the user. Here, the change in
   \boldsymbol{a}(\boldsymbol{x}, t)
   is not taken into account.

How it should be implemented: When using a compressible equilibrium with a force model, the user should be allowed to either specify

• a whole force density field, where the force density may change, so that
  \boldsymbol{F}(\boldsymbol{x},t)
• the acceleration
  \boldsymbol{a} = \text{const.}
  rather than a force density
  \boldsymbol{F}
  . This acceleration can be directly used in the computation of the macroscopic velocity, without needing to take the density into account. Finally, in the LBM collision, lbmpy should use the force density
  \rho(\boldsymbol{x},t) \boldsymbol{a}
  .

To illustrate the problem, I have reprinted the relevant equations from the book of T. Krüger, p. 233:

\boldsymbol{u}(\boldsymbol{x},t) = 1/\rho(\boldsymbol{x},t) \sum_{i}f_{i}(\boldsymbol{x},t) + \boldsymbol{F} \Delta t / 2 / \rho(\boldsymbol{x},t)

, where I assumed that a user specifies a constant force density

\boldsymbol{F}

for the whole domain.

\boldsymbol{S_{i}}(\boldsymbol{x},t) \sim \boldsymbol{F}

, where

\boldsymbol{S_{i}}

is a source term added to the collision operator. Note that the index

i

should not be bold, but the LaTeX interpreter here does not allow to write it correctly.

added Bug label

assigned to @holzer