diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index dfdd5c4c2c5dbd246d24a112a52089f3c4f4909e..c36f867c937aad8c56ff97b89657673ad73c0066 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 017927c37c8cb27f962497ab24770f35a11329de..0b96e844434556d409471b10c0f15f1998bfb8b4 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 70c9bdd5c99149546d55c2e01cc28321183c520b..adcdfcb65bdf6f1aac690fdc610935abafb0c8f1 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.