In this tutorial, we test a very classical lattice Boltzmann scheme \(\DdQq{2}{5}\) on the heat equation.
The problem reads
where \(\mu\) is a constant scalar.
The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximatation of the solution on discret points of \((0,1)^2\) at discret instants.
To simulate this system of equations, we use the \(\DdQq{2}{5}\) scheme given by
five velocities \(v_0=(0,0)\), \(v_1=(1,0)\), \(v_2=(0,1)\), \(v_3=(-1,0)\), and \(v_4=(0,-1)\) with associated distribution functions \(\fk{i}\), \(0\leq i\leq 4\),
a space step \(\dx\) and a time step \(\dt\), the ration \(\lambda=\dx/\dt\) is called the scheme velocity,
five moments
and their equilibrium values \(\mke{k}\), \(0\leq k\leq 4\). * two relaxation parameters \(s_1\) and \(s_2\) lying in \([0,2]\) (\(s_1\) for the odd moments and \(s_2\) for the odd ones).
In order to use the formalism of the package pyLBM, we introduce the five polynomials that define the moments: \(P_0 = 1\), \(P_1=X\), \(P_2=Y\), \(P_3=(X^2+Y^2)/2\), and \(P_4=(X^2-Y^2)/2\), such that
The transformation \((\fk{0}, \fk{1}, \fk{2}, \fk{3}, \fk{4})\mapsto(\mk{0},\mk{1}, \mk{2}, \mk{3}, \mk{4})\) is invertible if, and only if, the polynomials \((P_0,P_1,P_2,P_3,P_4)\) is a free set over the stencil of velocities.
The lattice Boltzmann method consists to compute the distribution functions \(\fk{i}\), \(0\leq i\leq 4\) in each point of the lattice \(x\) and at each time \(t^n=n\dt\). A step of the scheme can be read as a splitting between the relaxation phase and the transport phase:
relaxation:
m2f:
transport:
f2m:
The moment of order \(0\), \(\mk{0}\), being conserved during the relaxation phase, a diffusive scaling \(\dt=\dx^2\), yields to the following equivalent equation
if \(\mke{1}=0\). In order to be consistent with the heat equation, the following choice is done:
pyLBM uses Python dictionary to describe the simulation. In the following, we will build this dictionary step by step.
In pyLBM, the geometry is defined by a box and a label for the boundaries. We define here a square \((0, 1)^2\).
Geometry informations
spatial dimension: 2
bounds of the box:
[[ 0. 1.]
[ 0. 1.]]
pyLBM provides a class stencil that is used to define the discret velocities of the scheme. In this example, the stencil is composed by the velocities \(v_0=(0,0)\), \(v_1=(1,0)\), \(v_2=(-1,0)\), \(v_3=(0,1)\), and \(v_4=(0,-1)\) numbered by \([0,1,2,3,4]\).
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: 5
- velocities: (0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -1),
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).
We construct a domain with \(N=10\) points in space.
Domain informations
spatial dimension: 2
space step: dx= 1.000e-01
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:
(see the documentation for more details)
The scheme velocity could be taken to \(1/\dx\) and the inital value of \(u\) to
[0] WARNING pyLBM.scheme in function __init__ line 229
The value 'space_step' is not given or wrong.
The scheme takes default value: dx = 1.
WARNING:pyLBM.scheme:The value 'space_step' is not given or wrong.
The scheme takes default value: dx = 1.
Scheme informations
spatial dimension: dim=2
number of schemes: nscheme=1
number of velocities:
Stencil.nv[0]=5
velocities value:
v[0]=(0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -1),
polynomials:
P[0]=Matrix([[1], [X], [Y], [X**2/2 + Y**2/2], [X**2/2 - Y**2/2]])
equilibria:
EQ[0]=Matrix([[u], [0.0], [0.0], [0.5*u], [0.0]])
relaxation parameters:
s[0]=[0.0, 0.4, 0.4, 1.0, 1.0]
moments matrices
M = [Matrix([
[1, 1, 1, 1, 1],
[0, 1, 0, -1, 0],
[0, 0, 1, 0, -1],
[0, 1/2, 1/2, 1/2, 1/2],
[0, 1/2, -1/2, 1/2, -1/2]])]
invM = [Matrix([
[1, 0, 0, -2, 0],
[0, 1/2, 0, 1/2, 1/2],
[0, 0, 1/2, 1/2, -1/2],
[0, -1/2, 0, 1/2, 1/2],
[0, 0, -1/2, 1/2, -1/2]])]
A simulation is built by defining a correct dictionary.
We combine the previous dictionaries to build a simulation. In order to impose the homogeneous Dirichlet conditions in \(x=0\), \(x=1\), \(y=0\), and \(y=1\), the dictionary should contain the key ‘boundary_conditions’ (we use pyLBM.bc.Anti_bounce_back function).
Simulation informations:
Domain informations
spatial dimension: 2
space step: dx= 1.000e-01
Scheme informations
spatial dimension: dim=2
number of schemes: nscheme=1
number of velocities:
Stencil.nv[0]=5
velocities value:
v[0]=(0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -1),
polynomials:
P[0]=Matrix([[1], [X], [Y], [X**2/2 + Y**2/2], [X**2/2 - Y**2/2]])
equilibria:
EQ[0]=Matrix([[u], [0.0], [0.0], [0.5*u], [0.0]])
relaxation parameters:
s[0]=[0.0, 0.4, 0.4, 1.0, 1.0]
moments matrices
M = [Matrix([
[1, 1, 1, 1, 1],
[0, 1, 0, -1, 0],
[0, 0, 1, 0, -1],
[0, 1/2, 1/2, 1/2, 1/2],
[0, 1/2, -1/2, 1/2, -1/2]])]
invM = [Matrix([
[1, 0, 0, -2, 0],
[0, 1/2, 0, 1/2, 1/2],
[0, 0, 1/2, 1/2, -1/2],
[0, -1/2, 0, 1/2, 1/2],
[0, 0, -1/2, 1/2, -1/2]])]
Once the simulation is initialized, one time step can be performed by using the function one_time_step.
We compute the solution of the heat equation at \(t=0.1\). On the same graphic, we plot the initial condition, the exact solution and the numerical solution.