From b7b3b0caaab41f5fe093f7c9a9ce11bc2efe1b8b Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Thu, 25 Jun 2020 18:56:16 +0200
Subject: [PATCH] Schiller force model
---
lbmpy/creationfunctions.py | 1 +
lbmpy/forcemodels.py | 32 +++++++++++++++++++++++++++++++-
lbmpy/relaxationrates.py | 14 +++++++++++++-
3 files changed, 45 insertions(+), 2 deletions(-)
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index dfdd5c4c..c36f867c 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -479,6 +479,7 @@ def force_model_from_string(force_model_name, force_values):
'buick': forcemodels.Buick,
'silva': forcemodels.Buick,
'edm': forcemodels.EDM,
+ 'schiller': forcemodels.Schiller,
}
if force_model_name.lower() not in force_model_dict:
raise ValueError("Unknown force model %s" % (force_model_name,))
diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index 017927c3..0b96e844 100644
--- a/lbmpy/forcemodels.py
+++ b/lbmpy/forcemodels.py
@@ -88,7 +88,7 @@ second moment nonzero :class:`Luo` :class:`Guo`
import sympy as sp
-from lbmpy.relaxationrates import get_shear_relaxation_rate
+from lbmpy.relaxationrates import get_bulk_relaxation_rate, get_shear_relaxation_rate
class Simple:
@@ -148,6 +148,7 @@ class Guo:
luo = Luo(self._force)
shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
+ assert len(set(lb_method.relaxation_rates)) == 1, "Guo only works for SRT, use Schiller instead"
correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
return [correction_factor * t for t in luo(lb_method)]
@@ -158,6 +159,35 @@ class Guo:
return default_velocity_shift(density, self._force)
+class Schiller:
+ r"""
+ Force model by Schiller :cite:`schiller2008thermal`
+ Equivalent to Guo but not restricted to SRT.
+ """
+ def __init__(self, force):
+ self._force = force
+
+ def __call__(self, lb_method):
+ u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
+ force = sp.Matrix(self._force)
+
+ uf = u.dot(force) * sp.eye(len(force))
+ omega = get_shear_relaxation_rate(lb_method)
+ omega_bulk = get_bulk_relaxation_rate(lb_method)
+ G = (u * force.transpose() + force * u.transpose() - uf * sp.Rational(2, 3)) * sp.Rational(1, 2) * \
+ (2 + omega) + uf * sp.Rational(1, 3) * (2 + omega_bulk)
+
+ result = []
+ for direction, w_i in zip(lb_method.stencil, lb_method.weights):
+ direction = sp.Matrix(direction)
+ tr = sp.trace(G * (direction * direction.transpose() - sp.Rational(1, 3) * sp.eye(len(force))))
+ result.append(3 * w_i * (force.dot(direction) + sp.Rational(3, 2) * tr))
+ return result
+
+ def macroscopic_velocity_shift(self, density):
+ return default_velocity_shift(density, self._force)
+
+
class Buick:
r"""
This force model :cite:`buick2000gravity` has a force term with zero second moment.
diff --git a/lbmpy/relaxationrates.py b/lbmpy/relaxationrates.py
index 70c9bdd5..adcdfcb6 100644
--- a/lbmpy/relaxationrates.py
+++ b/lbmpy/relaxationrates.py
@@ -1,6 +1,6 @@
import sympy as sp
-from lbmpy.moments import is_shear_moment
+from lbmpy.moments import is_bulk_moment, is_shear_moment
def relaxation_rate_from_lattice_viscosity(nu):
@@ -46,6 +46,18 @@ def get_shear_relaxation_rate(method):
"Can not determine their relaxation rate automatically")
+def get_bulk_relaxation_rate(method):
+ """
+ The bulk moment is x^2 + y^2 + z^2, plus a constant for orthogonalization.
+ If this moment does not exist, the bulk relaxation is part of the shear relaxation.
+ The bulk relaxation rate determines the bulk viscosity in hydrodynamic LBM schemes.
+ """
+ for moment, relax_info in method.relaxation_info_dict.items():
+ if is_bulk_moment(moment, method.dim):
+ return relax_info.relaxation_rate
+ return get_shear_relaxation_rate(method)
+
+
def relaxation_rate_scaling(omega, level_scale_factor):
"""Computes adapted omega for refinement.
--
GitLab