diff --git a/creationfunctions.py b/creationfunctions.py
index ae5bbcaf2c3cc15816b618f429db92795e4c6131..74ab264142477a974a9ea43815007842134f9dd7 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -40,7 +40,7 @@ General:
compressible.
- ``equilibriumAccuracyOrder=2``: order in velocity, at which the equilibrium moment/cumulant approximation is
truncated. Order 2 is sufficient to approximate Navier-Stokes
-- ``forceModel=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'``. For details see
+- ``forceModel=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'`` or ``'buick'``. For details see
:mod:`lbmpy.forcemodels`
- ``force=(0,0,0)``: either constant force or a symbolic expression depending on field value
- ``useContinuousMaxwellianEquilibrium=True``: way to compute equilibrium moments/cumulants, if False the standard
@@ -422,8 +422,8 @@ def createLatticeBoltzmannMethod(**params):
forceModel = forceModels.Luo(params['force'][:dim])
elif forceModelName.lower() == 'guo':
forceModel = forceModels.Guo(params['force'][:dim])
- elif forceModelName.lower() == 'silva':
- forceModel = forceModels.Silva(params['force'][:dim])
+ elif forceModelName.lower() == 'silva' or forceModelName.lower() == 'buick':
+ forceModel = forceModels.Buick(params['force'][:dim])
else:
raise ValueError("Unknown force model %s" % (forceModelName,))
else:
diff --git a/forcemodels.py b/forcemodels.py
index ad0539e49c37cd8054e8395956f3a01a4c3d4b2e..39cc4e0425c21623d0dd33ca3609ccd7481ddce2 100644
--- a/forcemodels.py
+++ b/forcemodels.py
@@ -1,8 +1,90 @@
-"""
+r"""
.. module:: forcemodels
:synopsis: Collection of forcing terms for hydrodynamic LBM simulations
+
+Get started:
+-----------
+
+This module offers different models to introduce a body force in the lattice Boltzmann scheme.
+If you don't know which model to choose, use :class:`lbmpy.forcemodels.Guo`.
+For incompressible collision models the :class:`lbmpy.forcemodels.Buick` model can be better.
+
+
+Detailed information:
+---------------------
+
+Force models add a term :math:`C_F` to the collision equation:
+
+.. math ::
+
+ f(\pmb{x} + c_q \Delta t, t + \Delta t) - f(\pmb{x},t) = \Omega(f, f^{(eq)})
+ + \underbrace{F_q}_{\mbox{forcing term}}
+
+The form of this term depends on the concrete force model: the first moment of this forcing term is equal
+to the acceleration :math:`\pmb{a}` for all force models.
+
+.. math ::
+
+ \sum_q \pmb{c}_q F_q = \pmb{a}
+
+
+The second order moment is different for the forcing models - if it is zero the model is suited for
+incompressible flows. For weakly compressible collision operators a force model with a corrected second order moment
+should be chosen.
+
+.. math ::
+
+ \sum_q c_{qi} c_{qj} f_q = F_i u_j + F_j u_i \hspace{1cm} \mbox{for Guo, Luo models}
+
+ \sum_q c_{qi} c_{qj} f_q = 0 \hspace{1cm} \mbox{for Simple, Buick}
+
+Models with zero second order moment have:
+
+.. math ::
+
+ F_q = \frac{w_q}{c_s^2} c_{qi} \; a_i
+
+Models with nonzero second moment have:
+
+.. math ::
+
+ F_q = \frac{w_q}{c_s^2} c_{qi} \; a_i + \frac{w_q}{c_s^4} (c_{qi} c_{qj} - c_s^2 \delta_{ij} ) u_j \, a_i
+
+
+For all force models the computation of the macroscopic velocity has to be adapted (shifted) by adding a term
+:math:`S_{macro}` that we call "macroscopic velocity shift"
+
+ .. math ::
+
+ \pmb{u} = \sum_q \pmb{c}_q f_q + S_{macro}
+
+ S_{macro} = \frac{\Delta t}{2} \sum_q F_q
+
+
+Some models also shift the velocity entering the equilibrium distribution.
+
+Comparison
+----------
+
+Force models can be distinguished by 2 options:
+
+**Option 1**:
+ :math:`C_F = 1` and equilibrium velocity is not shifted, or :math:`C_F=(1 - \frac{\omega}{2})` and equilibrum is shifted.
+
+**Option 2**:
+ second velocity moment is zero or :math:`F_i u_j + F_j u_i`
+
+
+===================== ==================== =================
+ Option2 \\ Option1 no equilibrium shift equilibrium shift
+===================== ==================== =================
+second moment zero :class:`Simple` :class:`Buick`
+second moment nonzero :class:`Luo` :class:`Guo`
+===================== ==================== =================
+
"""
+
import sympy as sp
from lbmpy.methods.relaxationrates import getShearRelaxationRate
@@ -10,7 +92,7 @@ from lbmpy.methods.relaxationrates import getShearRelaxationRate
class Simple(object):
r"""
A simple force model which introduces the following additional force term in the
- collision process :math:`3 w_i e_i f_i` (often: force = rho * acceleration)
+ collision process :math:`\frac{w_q}{c_s^2} c_{qi} \; a_i` (often: force = rho * acceleration)
Should only be used with constant forces!
Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
"""
@@ -30,31 +112,9 @@ class Simple(object):
return defaultVelocityShift(density, self._force)
-class Silva(object):
- def __init__(self, force):
- self._force = force
-
- def __call__(self, lbMethod, **kwargs):
- simple = Simple(self._force)
-
- shearRelaxationRate = getShearRelaxationRate(lbMethod)
- correctionFactor = (1 - sp.Rational(1, 2) * shearRelaxationRate)
- return [correctionFactor * t for t in simple(lbMethod)]
-
- def macroscopicVelocityShift(self, density):
- return defaultVelocityShift(density, self._force)
-
- def equilibriumVelocityShift(self, density):
- return defaultVelocityShift(density, self._force)
-
-
class Luo(object):
r"""
- Force model by Luo with the following forcing term
-
- .. math ::
-
- F_i = w_i \left( \frac{c_i - u}{c_s^2} + \frac{c_i \left<c_i, u\right>}{c_s^4} \right) F
+ Force model by Luo :cite:`luo1993lattice`
Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
"""
@@ -77,12 +137,7 @@ class Luo(object):
class Guo(object):
r"""
- Force model by Guo with the following term:
-
- .. math ::
-
- F_i = w_i ( 1 - \frac{\omega}{2} ) \left( \frac{c_i - u}{c_s^2} + \frac{c_i \left<c_i, u\right>}{c_s^4} \right) F
-
+ Force model by Guo :cite:`guo2002discrete`
Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity (both shifted by F/2)!
"""
def __init__(self, force):
@@ -102,6 +157,31 @@ class Guo(object):
return defaultVelocityShift(density, self._force)
+class Buick(object):
+ r"""
+ This force model :cite:`buick2000gravity` has a force term with zero second moment.
+ It is suited for incompressible lattice models.
+
+
+ """
+
+ def __init__(self, force):
+ self._force = force
+
+ def __call__(self, lbMethod, **kwargs):
+ simple = Simple(self._force)
+
+ shearRelaxationRate = getShearRelaxationRate(lbMethod)
+ correctionFactor = (1 - sp.Rational(1, 2) * shearRelaxationRate)
+ return [correctionFactor * t for t in simple(lbMethod)]
+
+ def macroscopicVelocityShift(self, density):
+ return defaultVelocityShift(density, self._force)
+
+ def equilibriumVelocityShift(self, density):
+ return defaultVelocityShift(density, self._force)
+
+
# -------------------------------- Helper functions ------------------------------------------------------------------