Poiseuille flow¶

\[\renewcommand{\DdQq}[2]{{\mathrm D}_{#1}{\mathrm Q}_{#2}} \renewcommand{\drondt}{\partial_t} \renewcommand{\drondx}{\partial_x} \renewcommand{\drondy}{\partial_y} \renewcommand{\drondtt}{\partial_{tt}} \renewcommand{\drondxx}{\partial_{xx}} \renewcommand{\drondyy}{\partial_{yy}} \renewcommand{\dx}{\Delta x} \renewcommand{\dt}{\Delta t} \renewcommand{\grandO}{{\mathcal O}} \renewcommand{\density}[2]{\,f_{#1}^{#2}} \renewcommand{\fk}[1]{\density{#1}{\vphantom{\star}}} \renewcommand{\fks}[1]{\density{#1}{\star}} \renewcommand{\moment}[2]{\,m_{#1}^{#2}} \renewcommand{\mk}[1]{\moment{#1}{\vphantom{\star}}} \renewcommand{\mke}[1]{\moment{#1}{e}} \renewcommand{\mks}[1]{\moment{#1}{\star}}\]

In this tutorial, we consider the classical \(\DdQq{2}{9}\) to simulate a Poiseuille flow modeling by the Navier-Stokes equations.

The \(\DdQq{2}{9}\) for Navier-Stokes¶

The \(\DdQq{2}{9}\) is defined by:

a space step \(\dx\) and a time step \(\dt\) related to the scheme velocity \(\lambda\) by the relation \(\lambda=\dx/\dt\),
nine velocities \(\{(0,0), (\pm1,0), (0,\pm1), (\pm1, \pm1)\}\), identified in pyLBM by the numbers \(0\) to \(8\),
nine polynomials used to build the moments

\[\{1, \lambda X, \lambda Y, 3E-4, (9E^2-21E+8)/2, 3XE-5X, 3YE-5Y,X^2-Y^2, XY\},\]

where \(E = X^2+Y^2\).

three conserved moments \(\rho\), \(q_x\), and \(q_y\),
nine relaxation parameters (three are \(0\) corresponding to conserved moments): \(\{0,0,0,s_\mu,s_\mu,s_\eta,s_\eta,s_\eta,s_\eta\}\), where \(s_\mu\) and \(s_\eta\) are in \((0,2)\),
equilibrium value of the non conserved moments

\[\begin{split}\begin{aligned}\mke{3} &= -2\rho + 3(q_x^2+q_y^2)/(\rho_0\lambda^2), \\ \mke{4} &= \rho-3(q_x^2+q_y^2)/(\rho_0\lambda^2), \\ \mke{5} &= -q_x/\lambda, \\ \mke{6} &= -q_y/\lambda, \\ \mke{7} &= (q_x^2-q_y^2)/(\rho_0\lambda^2), \\ \mke{8} &= q_xq_y/(\rho_0\lambda^2),\end{aligned}\end{split}\]

where \(\rho_0\) is a given scalar.

This scheme is consistant at second order with the following equations (taken \(\rho_0=1\)) \begin{gathered} :raw-latex:`\drondt`:raw-latex:`\rho `+ :raw-latex:`drondx q_x + :raw-latex:drondy q_y = 0,\ :raw-latex:drondt q_x + :raw-latex:drondx (q_x^2+p) + :raw-latex:drondy (q_xq_y) = :raw-latex:mu :raw-latex:drondx (:raw-latex:drondx q_x + :raw-latex:drondy q_y ) + :raw-latex:eta (:raw-latex:drondxx`+:raw-latex:drondyy)q_x, \ :raw-latex:`\drondt `q\_y + :raw-latex:`drondx (q_xq_y) + :raw-latex:drondy (q_y^2+p) = :raw-latex:mu :raw-latex:drondy (:raw-latex:drondx q_x + :raw-latex:drondy q_y ) + :raw-latex:eta (:raw-latex:drondxx`+:raw-latex:drondyy)q_y,\end{gathered} with \(p=\rho\lambda^2/3\).

Build the simulation with pyLBM¶

In the following, we build the dictionary of the simulation step by step.

The geometry¶

The simulation is done on a rectangle of length \(L\) and width \(W\). We can use \(L=W=1\).

We propose a dictionary that build the geometry of the domain. The labels of the bounds can be specified to different values for the moment.

Geometry informations
     spatial dimension: 2
     bounds of the box:
[[ 0.   1. ]
 [-0.5  0.5]]

The stencil¶

The stencil of the \(\DdQq{2}{9}\) is composed by the nine following velocities in 2D:

\[\begin{split}\begin{gathered}v_0=(0,0),\\ v_1=(1,0), \quad v_2=(0,1), \quad v_3=(-1,0), \quad v_4=(0,-1), \\ v_5=(1,1), \quad v_6=(-1,1), \quad v_7=(-1,-1), \quad v_8=(1,-1).\end{gathered}\end{split}\]

Stencil informations
     * spatial dimension: 2
     * maximal velocity in each direction: [1 1]
     * minimal velocity in each direction: [-1 -1]
     * Informations for each elementary stencil:
            stencil 0
             - number of velocities:  9
             - velocities: (0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -1), (5: 1, 1), (6: -1, 1), (7: -1, -1), (8: 1, -1),

The domain¶

In order to build the domain of the simulation, the dictionary should contain the space step \(\dx\) and the stencils of the velocities (one for each scheme).

Domain informations
     spatial dimension: 2
     space step: dx= 1.000e-01

The scheme¶

In pyLBM, a simulation can be performed by using several coupled schemes. In this example, a single scheme is used and defined through a list of one single dictionary. This dictionary should contain:

‘velocities’: a list of the velocities
‘conserved_moments’: a list of the conserved moments as sympy variables
‘polynomials’: a list of the polynomials that define the moments
‘equilibrium’: a list of the equilibrium value of all the moments
‘relaxation_parameters’: a list of the relaxation parameters (\(0\) for the conserved moments)
‘init’: a dictionary to initialize the conserved moments

(see the documentation for more details)

In order to fix the bulk (\(\mu\)) and the shear (\(\eta\)) viscosities, we impose

\[s_\eta = \frac{2}{1+\eta d}, \qquad s_\mu = \frac{2}{1+\mu d}, \qquad d = \frac{6}{\lambda\rho_0\dx}.\]

The scheme velocity could be taken to \(1\) and the inital value of \(\rho\) to \(\rho_0=1\), \(q_x\) and \(q_y\) to \(0\).

In order to accelerate the simulation, we can use another generator. The default generator is Numpy (pure python). We can use for instance Cython that generates a more efficient code. This generator can be activated by using ‘generator’: pyLBM.generator.CythonGenerator in the dictionary.

Scheme informations
     spatial dimension: dim=2
     number of schemes: nscheme=1
     number of velocities:
    Stencil.nv[0]=9
     velocities value:
    v[0]=(0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -1), (5: 1, 1), (6: -1, 1), (7: -1, -1), (8: 1, -1),
     polynomials:
    P[0]=Matrix([[1], [LA*X], [LA*Y], [3*X**2 + 3*Y**2 - 4], [-21*X**2/2 - 21*Y**2/2 + 9*(X**2 + Y**2)**2/2 + 4], [3*X*(X**2 + Y**2) - 5*X], [3*Y*(X**2 + Y**2) - 5*Y], [X**2 - Y**2], [X*Y]])
     equilibria:
    EQ[0]=Matrix([[rho], [qx], [qy], [-2*rho + 3.0*qx**2/LA**2 + 3.0*qy**2/LA**2], [rho - 3.0*qx**2/LA**2 - 3.0*qy**2/LA**2], [-qx/LA], [-qy/LA], [1.0*qx**2/LA**2 - 1.0*qy**2/LA**2], [1.0*qx*qy/LA**2]])
     relaxation parameters:
    s[0]=[0.0, 0.0, 0.0, 1.1312217194570136, 1.1312217194570136, 0.025706940874035987, 0.025706940874035987, 0.025706940874035987, 0.025706940874035987]
     moments matrices
M = [Matrix([
[ 1,  1,  1,   1,   1,  1,   1,   1,   1],
[ 0, LA,  0, -LA,   0, LA, -LA, -LA,  LA],
[ 0,  0, LA,   0, -LA, LA,  LA, -LA, -LA],
[-4, -1, -1,  -1,  -1,  2,   2,   2,   2],
[ 4, -2, -2,  -2,  -2,  1,   1,   1,   1],
[ 0, -2,  0,   2,   0,  1,  -1,  -1,   1],
[ 0,  0, -2,   0,   2,  1,   1,  -1,  -1],
[ 0,  1, -1,   1,  -1,  0,   0,   0,   0],
[ 0,  0,  0,   0,   0,  1,  -1,   1,  -1]])]
invM = [Matrix([
[1/9,         0,         0,  -1/9,   1/9,     0,     0,    0,    0],
[1/9,  1/(6*LA),         0, -1/36, -1/18,  -1/6,     0,  1/4,    0],
[1/9,         0,  1/(6*LA), -1/36, -1/18,     0,  -1/6, -1/4,    0],
[1/9, -1/(6*LA),         0, -1/36, -1/18,   1/6,     0,  1/4,    0],
[1/9,         0, -1/(6*LA), -1/36, -1/18,     0,   1/6, -1/4,    0],
[1/9,  1/(6*LA),  1/(6*LA),  1/18,  1/36,  1/12,  1/12,    0,  1/4],
[1/9, -1/(6*LA),  1/(6*LA),  1/18,  1/36, -1/12,  1/12,    0, -1/4],
[1/9, -1/(6*LA), -1/(6*LA),  1/18,  1/36, -1/12, -1/12,    0,  1/4],
[1/9,  1/(6*LA), -1/(6*LA),  1/18,  1/36,  1/12, -1/12,    0, -1/4]])]

Run the simulation¶

For the simulation, we take

The domain \(x\in(0, L)\) and \(y\in(-W/2,W/2)\), \(L=2\), \(W=1\),
the viscosities \(\mu=10^{-2}=\eta=10^{-2}\),
the space step \(\dx=1/128\), the scheme velocity \(\lambda=1\),
the mean density \(\rho_0=1\).

Concerning the boundary conditions, we impose the velocity on all the edges by a bounce-back condition with a source term that reads

\[ \begin{align}\begin{aligned}q_x(x, y) = \rho_0 v_{\text{max}} \Bigl( 1 - \frac{4y^2}{W^2} \Bigr), \qquad q_y(x, y) = 0,\\with :math:`v_{\text{max}}=0.1`.\end{aligned}\end{align} \]

We compute the solution for \(t\in(0,50)\) and we plot several slices of the solution during the simulation.

This problem has an exact solution given by

\[ \begin{align}\begin{aligned}q_x = \rho_0 v_{\text{max}} \Bigl( 1 - \frac{4y^2}{W^2} \Bigr), \qquad q_y = 0, \qquad p = p_0 + K x,\\where the pressure gradient :math:`K` reads\end{aligned}\end{align} \]

\[K = -\frac{8 v_{\text{max}} \eta}{W^2}.\]

We compute the exact and the numerical gradients of the pressure.

Exact pressure gradient    : -8.000e-03
Numerical pressure gradient: -7.074e-03