[1]:
from lbmpy.session import *
from lbmpy.relaxationrates import *

Tutorial 06: Modifying a LBM method: Smagorinsky model

In this demo, we show how to modify a lattice Boltzmann method. As example we are going to add a simple turbulence model, by introducing a rule that locally computes the relaxation parameter dependent on the local strain rate tensor. The Smagorinsky model is implemented directly in lbmpy as well, however here we take the manual approach to demonstrate how a LB method can be changed in lbmpy.

1) Theoretical background

Since we have sympy available, we want to start out with the basic model equations and derive the concrete equations ourselves. This approach is less error prone, since the calculations are done by the computer algebra system, and oftentimes this approach is also more general and easier to understand.

a) Smagorinsky model

The basic idea of the Smagorinsky turbulence model is to safe compute time, by not resolving the smallest eddies of the flow on the grid, but model them by an artifical dissipation term. The energy dissipation of small scale vortices is taken into account by introducing a “turbulent viscosity”. This additional viscosity depends on local flow properties, namely the local shear rates. The larger the local shear rates the higher the turbulent viscosity and the more artifical dissipation is added.

The total viscosity is

\[\nu_{total} = \nu_0 + \underbrace{(C_S \Delta)^2 |S|}_{\nu_{t}}\]

where \(\nu_0\) is the normal viscosity, \(C_S\) is the Smagorinsky constant, not to be confused with the speed of sound! Typical values of the Smagorinsky constant are between 0.1 - 0.2. The filter length \(\Delta\) is chosen as 1 in lattice coordinates.

The quantity \(|S|\) is computed from the local strain rate tensor \(S\) that is given by

\[S_{ij} = \frac{1}{2} \left( \partial_i u_j + \partial_j u_i \right)\]

and

\[|S| = \sqrt{2 S_{ij} S_{ij}}\]

b) LBM implementation of Smagorinsky model

To add the Smagorinsky model to a LB scheme one has to first compute the strain rate tensor \(S_{ij}\) in each cell, and compute the turbulent viscosity \(\nu_t\) from it. Then the local relaxation rate has to be adapted to match the total viscosity \(\nu_{total}\) instead of the standard viscosity \(\nu_0\).

A fortunate property of LB methods is, that the strain rate tensor can be computed locally from the non-equilibrium part of the distribution function. This is somewhat surprising, since the strain rate tensor contains first order derivatives. The strain rate tensor can be obtained by

\[S_{ij} = - \frac{3 \omega_s}{2 \rho_{(0)}} \Pi_{ij}^{(neq)}\]

where \(\omega_s\) is the relaxation rate that determines the viscosity, \(\rho_{(0)}\) is \(\rho\) in compressible models and \(1\) for incompressible schemes. \(\Pi_{ij}^{(neq)}\) is the second order moment tensor of the non-equilibrium part of the distribution functions \(f^{(neq)} = f - f^{(eq)}\) and can be computed as

\[\Pi_{ij}^{(neq)} = \sum_q c_{qi} c_{qj} \; f_q^{(neq)}\]

We first have to find a closed form for \(S_{ij}\) since in the formula above, it depends on \(\omega\), which should be adapated according to \(S_{ij}\). So we compute \(\omega\) and insert it into the formula for \(S\):

[2]:
τ_0, ρ, ω, ω_total, ω_0 = sp.symbols("tau_0 rho omega omega_total omega_0", positive=True, real=True)
ν_0, C_S, S, Π = sp.symbols("nu_0, C_S, |S|, Pi", positive=True, real=True)

Seq = sp.Eq(S, 3 * ω / 2 * Π)
Seq
[2]:
$\displaystyle |S| = \frac{3 \Pi \omega}{2}$

Note that we left of the minus, since we took the absolute value of both tensor. The absolute value is defined as above, with the factor of two inside the square root. The \(\rho_{(0)}\) has been left out, remembering that \(\Pi^{(neq)}\) has to be divided by \(\rho\) in case of compressible models|.

Next, we compute \(\omega\) from the total viscosity as given by the Smagorinsky equation:

[3]:
relaxation_rate_from_lattice_viscosity(ν_0 + C_S ** 2 * S)

[3]:
$\displaystyle \frac{2}{6 C_{S}^{2} |S| + 6 \nu_{0} + 1}$

and insert it into the equation for \(|S|\)

[4]:
Seq2 = Seq.subs(ω, relaxation_rate_from_lattice_viscosity(ν_0 + C_S **2 * S ))
Seq2
[4]:
$\displaystyle |S| = \frac{3 \Pi}{6 C_{S}^{2} |S| + 6 \nu_{0} + 1}$

This equation contains only known quantities, such that we can solve it for \(|S|\). Additionally we substitute the lattice viscosity \(\nu_0\) by the original relaxation time \(\tau_0\). The resulting equations get simpler using relaxation times instead of rates.

[5]:
solveRes = sp.solve(Seq2, S)
assert len(solveRes) == 1
SVal = solveRes[0]
SVal = SVal.subs(ν_0, lattice_viscosity_from_relaxation_rate(1 / τ_0)).expand()
SVal
[5]:
$\displaystyle - \frac{\tau_{0}}{6 C_{S}^{2}} + \frac{\sqrt{72 C_{S}^{2} \Pi + 4 \tau_{0}^{2}}}{12 C_{S}^{2}}$

Knowning \(|S|\) we can compute the total relaxation time using

\[\nu_{total} = \nu_0 +C_S^2 |S|\]
[6]:
τ_val = 1 / (relaxation_rate_from_lattice_viscosity(lattice_viscosity_from_relaxation_rate(1/τ_0) + C_S**2 * SVal)).cancel()
τ_val
[6]:
$\displaystyle \frac{\tau_{0}}{2} + \frac{\sqrt{18 C_{S}^{2} \Pi + \tau_{0}^{2}}}{2}$

To compute \(\Pi^{(neq)}\) we use the following functions:

[7]:
def second_order_moment_tensor(function_values, stencil):
    assert len(function_values) == len(stencil)
    dim = len(stencil[0])
    return sp.Matrix(dim, dim, lambda i, j: sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))


def frobenius_norm(matrix, factor=1):
    return sp.sqrt(sum(i*i for i in matrix) * factor)

In the next cell we construct equations that take an standard relaxation rate \(\omega_0\) and compute a new relaxation rate \(\omega_{total}\) according to the Smagorinksy model, using τ_val computed above

[8]:
def smagorinsky_equations(ω_0, ω_total, method):
    f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
    return [sp.Eq(τ_0, 1 / ω_0),
            sp.Eq(Π, frobenius_norm(second_order_moment_tensor(f_neq, method.stencil), factor=2)),
            sp.Eq(ω_total, 1 / τ_val)]


smagorinsky_equations(ω_0, ω_total, create_lb_method())
[8]:
$\displaystyle \left[ \tau_{0} = \frac{1}{\omega_{0}}, \ \Pi = \sqrt{4 \left(- f_{5} + f_{6} + f_{7} - f_{8} - u_{0} u_{1}\right)^{2} + 2 \left(- \frac{\delta_{\rho}}{3} + f_{1} + f_{2} + f_{5} + f_{6} + f_{7} + f_{8} - u_{1}^{2}\right)^{2} + 2 \left(- \frac{\delta_{\rho}}{3} + f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8} - u_{0}^{2}\right)^{2}}, \ \omega_{total} = \frac{1}{\frac{\tau_{0}}{2} + \frac{\sqrt{18 C_{S}^{2} \Pi + \tau_{0}^{2}}}{2}}\right]$

2) Application: Channel flow

Next we modify a lbmpy scenario to use the Smagorinsky model. We create a MRT method, where we fix all relaxation rates except the relaxation rate that controls the viscosity.

[9]:
lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, force=(1e-6, 0),
                       force_model=ForceModel.LUO, relaxation_rates=[ω, 1.9, 1.9, 1.9])

method = create_lb_method(lbm_config=lbm_config)
method
[9]:
Moment-Based Method Stencil: D2Q9 Zero-Centered Storage: ✓ Force Model: Luo
Continuous Hydrodynamic Maxwellian Equilibrium $f (\rho, \left( u_{0}, \ u_{1}\right), \left( v_{0}, \ v_{1}\right)) = \frac{3 \delta_{\rho} e^{- \frac{3 v_{0}^{2}}{2} - \frac{3 v_{1}^{2}}{2}}}{2 \pi} - \frac{3 e^{- \frac{3 v_{0}^{2}}{2} - \frac{3 v_{1}^{2}}{2}}}{2 \pi} + \frac{3 e^{- \frac{3 \left(- u_{0} + v_{0}\right)^{2}}{2} - \frac{3 \left(- u_{1} + v_{1}\right)^{2}}{2}}}{2 \pi}$
Compressible: ✗ Deviation Only: ✓ Order: 2
Relaxation Info
Moment Eq. Value Relaxation Rate
$1$ $\delta_{\rho}$ $0.0$
$x$ $u_{0}$ $0.0$
$y$ $u_{1}$ $0.0$
$x^{2} - y^{2}$ $u_{0}^{2} - u_{1}^{2}$ $\omega$
$x y$ $u_{0} u_{1}$ $\omega$
$3 x^{2} + 3 y^{2} - 2$ $3 u_{0}^{2} + 3 u_{1}^{2}$ $1.9$
$3 x^{2} y - y$ $0$ $1.9$
$3 x y^{2} - x$ $0$ $1.9$
$9 x^{2} y^{2} - 3 x^{2} - 3 y^{2} + 1$ $0$ $1.9$

Only the collision rule has to be changed. Thus we first construct the collision rule, add the Smagorinsky equations and create a normal scenario from the modified collision rule. To avoid that the macroscopic quantity symbols in the Smagorinsky equations fall prey to optimization, we must disable simplification:

[10]:
optimization = {'simplification' : False}
collision_rule = create_lb_collision_rule(lb_method=method, optimization=optimization)
collision_rule = collision_rule.new_with_substitutions({ω: ω_total})

collision_rule.subexpressions += smagorinsky_equations(ω, ω_total, method)
collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
collision_rule
[10]:
Subexpressions:
$$rr_{0} \leftarrow 0.0$$
$$rr_{1} \leftarrow 1.9$$
$$F_{x} \leftarrow 1.0 \cdot 10^{-6}$$
$$F_{y} \leftarrow 0$$
$$partial_{m m1 e 0} \leftarrow f_{3} + f_{5} + f_{7}$$
$$partial_{m 0 e 0} \leftarrow f_{0} + f_{1} + f_{2}$$
$$partial_{m 1 e 0} \leftarrow f_{4} + f_{6} + f_{8}$$
$$partial_{m m1 e 1} \leftarrow f_{5} - f_{7}$$
$$partial_{m 0 e 1} \leftarrow f_{1} - f_{2}$$
$$partial_{m 1 e 1} \leftarrow f_{6} - f_{8}$$
$$partial_{m m1 e 2} \leftarrow f_{5} + f_{7}$$
$$partial_{m 0 e 2} \leftarrow f_{1} + f_{2}$$
$$partial_{m 1 e 2} \leftarrow f_{6} + f_{8}$$
$$m_{00} \leftarrow partial_{m 0 e 0} + partial_{m 1 e 0} + partial_{m m1 e 0}$$
$$m_{10} \leftarrow partial_{m 1 e 0} - partial_{m m1 e 0}$$
$$m_{01} \leftarrow partial_{m 0 e 1} + partial_{m 1 e 1} + partial_{m m1 e 1}$$
$$m_{20} \leftarrow partial_{m 1 e 0} + partial_{m m1 e 0}$$
$$m_{02} \leftarrow partial_{m 0 e 2} + partial_{m 1 e 2} + partial_{m m1 e 2}$$
$$m_{11} \leftarrow partial_{m 1 e 1} - partial_{m m1 e 1}$$
$$m_{21} \leftarrow partial_{m 1 e 1} + partial_{m m1 e 1}$$
$$m_{12} \leftarrow partial_{m 1 e 2} - partial_{m m1 e 2}$$
$$m_{22} \leftarrow partial_{m 1 e 2} + partial_{m m1 e 2}$$
$$\delta_{\rho} \leftarrow m_{00}$$
$$u_{0} \leftarrow m_{10}$$
$$u_{1} \leftarrow m_{01}$$
$$\rho \leftarrow \delta_{\rho} + 1$$
$$M_{0} \leftarrow m_{00}$$
$$M_{1} \leftarrow m_{10}$$
$$M_{2} \leftarrow m_{01}$$
$$M_{3} \leftarrow - m_{02} + m_{20}$$
$$M_{4} \leftarrow m_{11}$$
$$M_{5} \leftarrow - 2 m_{00} + 3 m_{02} + 3 m_{20}$$
$$M_{6} \leftarrow - m_{01} + 3 m_{21}$$
$$M_{7} \leftarrow - m_{10} + 3 m_{12}$$
$$M_{8} \leftarrow m_{00} - 3 m_{02} - 3 m_{20} + 9 m_{22}$$
$$M_{post 0} \leftarrow M_{0} + rr_{0} \left(- M_{0} + \delta_{\rho}\right)$$
$$M_{post 1} \leftarrow F_{x} + M_{1} + rr_{0} \left(- M_{1} + u_{0}\right)$$
$$M_{post 2} \leftarrow F_{y} + M_{2} + rr_{0} \left(- M_{2} + u_{1}\right)$$
$$M_{post 5} \leftarrow 6 F_{x} u_{0} + 6 F_{y} u_{1} + M_{5} + rr_{1} \left(- M_{5} + 3 u_{0}^{2} + 3 u_{1}^{2}\right)$$
$$M_{post 6} \leftarrow - M_{6} rr_{1} + M_{6}$$
$$M_{post 7} \leftarrow - M_{7} rr_{1} + M_{7}$$
$$M_{post 8} \leftarrow - M_{8} rr_{1} + M_{8}$$
$$sub_{k to f 18} \leftarrow \frac{1}{2}$$
$$sub_{k to f 19} \leftarrow \frac{1}{3}$$
$$sub_{k to f 20} \leftarrow \frac{1}{6}$$
$$sub_{k to f 21} \leftarrow \frac{1}{9}$$
$$sub_{k to f 22} \leftarrow \frac{1}{4}$$
$$sub_{k to f 0} \leftarrow M_{post 0}$$
$$sub_{k to f 1} \leftarrow M_{post 1}$$
$$sub_{k to f 2} \leftarrow M_{post 2}$$
$$sub_{k to f 5} \leftarrow M_{post 5}$$
$$sub_{k to f 6} \leftarrow M_{post 6}$$
$$sub_{k to f 7} \leftarrow M_{post 7}$$
$$sub_{k to f 8} \leftarrow M_{post 8}$$
$$m_{post 00} \leftarrow sub_{k to f 0}$$
$$m_{post 10} \leftarrow sub_{k to f 1}$$
$$m_{post 01} \leftarrow sub_{k to f 2}$$
$$m_{post 21} \leftarrow sub_{k to f 19} \left(sub_{k to f 2} + sub_{k to f 6}\right)$$
$$m_{post 12} \leftarrow sub_{k to f 19} \left(sub_{k to f 1} + sub_{k to f 7}\right)$$
$$m_{post 22} \leftarrow sub_{k to f 21} \left(sub_{k to f 0} + sub_{k to f 5} + sub_{k to f 8}\right)$$
$$sub_{k to f 11} \leftarrow sub_{k to f 18} \left(m_{post 01} - m_{post 21}\right)$$
$$sub_{k to f 13} \leftarrow sub_{k to f 18} \left(- m_{post 10} + m_{post 12}\right)$$
$$sub_{k to f 15} \leftarrow sub_{k to f 22} \left(- m_{post 12} + m_{post 21}\right)$$
$$sub_{k to f 17} \leftarrow sub_{k to f 22} \left(m_{post 12} + m_{post 21}\right)$$
$$\tau_{0} = \frac{1}{\omega}$$
$$\Pi = \sqrt{4 \left(- f_{5} + f_{6} + f_{7} - f_{8} - u_{0} u_{1}\right)^{2} + 2 \left(- \frac{\delta_{\rho}}{3} + f_{1} + f_{2} + f_{5} + f_{6} + f_{7} + f_{8} - u_{1}^{2}\right)^{2} + 2 \left(- \frac{\delta_{\rho}}{3} + f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8} - u_{0}^{2}\right)^{2}}$$
$$\omega_{total} = \frac{1}{\frac{\tau_{0}}{2} + \frac{\sqrt{18 C_{S}^{2} \Pi + \tau_{0}^{2}}}{2}}$$
$$M_{post 3} \leftarrow 2 F_{x} u_{0} - 2 F_{y} u_{1} + M_{3} + \omega_{total} \left(- M_{3} + u_{0}^{2} - u_{1}^{2}\right)$$
$$M_{post 4} \leftarrow F_{x} u_{1} + F_{y} u_{0} + M_{4} + \omega_{total} \left(- M_{4} + u_{0} u_{1}\right)$$
$$sub_{k to f 3} \leftarrow M_{post 3}$$
$$sub_{k to f 4} \leftarrow M_{post 4}$$
$$m_{post 20} \leftarrow sub_{k to f 0} sub_{k to f 19} + sub_{k to f 18} sub_{k to f 3} + sub_{k to f 20} sub_{k to f 5}$$
$$m_{post 02} \leftarrow sub_{k to f 0} sub_{k to f 19} - sub_{k to f 18} sub_{k to f 3} + sub_{k to f 20} sub_{k to f 5}$$
$$m_{post 11} \leftarrow sub_{k to f 4}$$
$$sub_{k to f 10} \leftarrow sub_{k to f 18} \left(m_{post 02} - m_{post 22}\right)$$
$$sub_{k to f 12} \leftarrow sub_{k to f 18} \left(m_{post 20} - m_{post 22}\right)$$
$$sub_{k to f 14} \leftarrow sub_{k to f 22} \left(- m_{post 11} + m_{post 22}\right)$$
$$sub_{k to f 16} \leftarrow sub_{k to f 22} \left(m_{post 11} + m_{post 22}\right)$$
Main Assignments:
$$d_{0} \leftarrow m_{post 00} - m_{post 02} - m_{post 20} + m_{post 22}$$
$$d_{1} \leftarrow sub_{k to f 10} + sub_{k to f 11}$$
$$d_{2} \leftarrow sub_{k to f 10} - sub_{k to f 11}$$
$$d_{3} \leftarrow sub_{k to f 12} + sub_{k to f 13}$$
$$d_{4} \leftarrow sub_{k to f 12} - sub_{k to f 13}$$
$$d_{5} \leftarrow sub_{k to f 14} + sub_{k to f 15}$$
$$d_{6} \leftarrow sub_{k to f 16} + sub_{k to f 17}$$
$$d_{7} \leftarrow sub_{k to f 16} - sub_{k to f 17}$$
$$d_{8} \leftarrow sub_{k to f 14} - sub_{k to f 15}$$

In the next cell the collision rule is simplified by extracting common subexpressions

[11]:
from pystencils.simp import sympy_cse
#collision_rule = sympy_cse(collision_rule)

A channel scenario can be created from a modified collision rule:

[12]:
ch = create_channel((300, 100), force=1e-6, collision_rule=collision_rule,
                    kernel_params={"C_S": 0.12, "omega": 1.999})
[13]:
#show_code(ch.ast)
[14]:
ch.run(5000)
[15]:
plt.figure(dpi=200)
plt.vector_field(ch.velocity[:, :])
np.max(ch.velocity[:, :])
[15]:
$\displaystyle 0.00504266401703371$
../_images/notebooks_06_tutorial_modifying_method_smagorinsky_26_1.png

Appendix: Strain rate tensor formula from Chapman Enskog

The connection between \(S_{ij}\) and \(\Pi_{ij}^{(neq)}\) can be seen using a Chapman Enskog expansion. Since lbmpy has a module that automatically does this expansions we can have a look at it:

[16]:
from lbmpy.chapman_enskog import ChapmanEnskogAnalysis, CeMoment
from lbmpy.chapman_enskog.chapman_enskog import remove_higher_order_u
compressible_model = create_lb_method(stencil=Stencil.D2Q9, compressible=True, zero_centered=False)
incompressible_model = create_lb_method(stencil=Stencil.D2Q9, compressible=False, zero_centered=False)

ce_compressible = ChapmanEnskogAnalysis(compressible_model)
ce_incompressible = ChapmanEnskogAnalysis(incompressible_model)

The Chapman Enskog analysis yields expresssions for the moment

\(\Pi = \Pi^{(eq)} + \epsilon \Pi^{(1)} + \epsilon^2 \Pi^{(2)} \cdots\) and the strain rate tensor is related to \(\Pi^{(1)}\). However the best approximation we have for \(\Pi^{(1)}\) is \(\Pi^{(neq)}\). For details, see the paper “Shear stress in lattice Boltzmann simulations” by Krüger, Varnik and Raabe from 2009.

Lets look at the values of \(\Pi^{(1)}\) obtained from the Chapman enskog expansion:

[17]:
Π_1_xy = CeMoment("\\Pi", moment_tuple=(1,1), superscript=1)
Π_1_xx = CeMoment("\\Pi", moment_tuple=(2,0), superscript=1)
Π_1_yy = CeMoment("\\Pi", moment_tuple=(0,2), superscript=1)
components = (Π_1_xx, Π_1_yy, Π_1_xy)

Π_1_xy_val = ce_compressible.higher_order_moments[Π_1_xy]
Π_1_xy_val
[17]:
$\displaystyle \frac{\rho u_{0}^{2} {\partial^{(1)}_{0} u_{1}} + 2 \rho u_{0} u_{1} {\partial^{(1)}_{0} u_{0}} + 2 \rho u_{0} u_{1} {\partial^{(1)}_{1} u_{1}} + \rho u_{1}^{2} {\partial^{(1)}_{1} u_{0}} - \frac{\rho {\partial^{(1)}_{1} u_{0}}}{3} - \frac{\rho {\partial^{(1)}_{0} u_{1}}}{3} + u_{0}^{2} u_{1} {\partial^{(1)}_{0} \rho} + u_{0} u_{1}^{2} {\partial^{(1)}_{1} \rho}}{\omega}$

This term has lots of higher order error terms in it. We assume that \(u\) is small in lattice coordinates, so if we neglect all terms in \(u\) that are quadratic or higher we get:

[18]:
remove_higher_order_u(Π_1_xy_val.expand())
[18]:
$\displaystyle - \frac{\rho {\partial^{(1)}_{1} u_{0}}}{3 \omega} - \frac{\rho {\partial^{(1)}_{0} u_{1}}}{3 \omega}$

Putting these steps together into a function, we can display them for the different cases quickly:

[19]:
def get_Π_1(ce_analysis, component):
    val = ce_analysis.higher_order_moments[component]
    return remove_higher_order_u(val.expand())

Compressible case:

[20]:
tuple(get_Π_1(ce_compressible, Pi) for Pi in components)
[20]:
$\displaystyle \left( - \frac{2 \rho {\partial^{(1)}_{0} u_{0}}}{3 \omega}, \ - \frac{2 \rho {\partial^{(1)}_{1} u_{1}}}{3 \omega}, \ - \frac{\rho {\partial^{(1)}_{1} u_{0}}}{3 \omega} - \frac{\rho {\partial^{(1)}_{0} u_{1}}}{3 \omega}\right)$

Incompressible case:

[21]:
tuple(get_Π_1(ce_incompressible, Pi) for Pi in components)
[21]:
$\displaystyle \left( \frac{2 u_{0} {\partial^{(1)}_{0} \rho}}{3 \omega} - \frac{2 {\partial^{(1)}_{0} u_{0}}}{3 \omega}, \ \frac{2 u_{1} {\partial^{(1)}_{1} \rho}}{3 \omega} - \frac{2 {\partial^{(1)}_{1} u_{1}}}{3 \omega}, \ \frac{u_{0} {\partial^{(1)}_{1} \rho}}{3 \omega} + \frac{u_{1} {\partial^{(1)}_{0} \rho}}{3 \omega} - \frac{{\partial^{(1)}_{1} u_{0}}}{3 \omega} - \frac{{\partial^{(1)}_{0} u_{1}}}{3 \omega}\right)$

In the incompressible case has some terms \(\partial \rho\) which are zero, since \(\rho\) is assumed constant.

Leaving out the error terms we finally obtain:

\[\Pi_{ij}^{(neq)} \approx \Pi_{ij}^{(1)} = -\frac{2 \rho_{(0)}}{3 \omega_s} \left( \partial_i u_j + \partial_j u_i \right)\]