Day 10 Mesoscale: LBE Tutorial Exercise 1: Single component LBE simulations¶
We start our LBE tutorials by looking at comparatively simple systems, using single-relaxation-time Bhatnagar-Gross-Krook (BGK) collisions - giving LBE with BGK or LBGK - and two-dimensional simulations to keep computer program execution times sensible. 3D LBE is a straightforward extension of 2D LBE, and often the 2D results you will obtain in these sessions may be regarded as sample ‘sections’ or slices from a 3D system, which is uniform in the direction perpendicular to the plane of the simulation.
The principal learning objectives for this tutorial are:
to familiarise you with the structure of a typical LBE algorithm,
to underscore the ‘kinetic length scale’ of the description,
to illustrate the simplicity and versatility of the LBE method, and
to illustrate strengths and shortcomings.
The flow can be visualised using the scalar stream function, \(\psi\). At steady state, contours of constant ‘rectangular stream function’ values represent fluid particle or, if you prefer, dye paths in the fluid.
To get started, access the working directory for this exercise:
~/WORKSHOP/Day_10Meso
and then open the notebook Day10LBETutorial1Ex1.ipynb
, and follow the instructions within it. While this notebook is fully contained, we include some additional background information here for you to look at while running LBE calculations. For historic reasons, we also provide some alternative (deprecated) scripts to plot results generated by the code:
lbe1plots.m
to use in MATLAB or Octave with comma-separated value (CSV) fileslbe1plots.py
to use Python (with Matplotlib) with comma-separated value (CSV) fileslbe1plots.gp
to use in Gnuplot with fixed-width delimited files ending with.dat
but we have already provided more extensive plotting scripts in the notebook.
Simulation Geometry: Lid-driven cavity¶
This 2D code is designed to work around the ‘default’ periodic boundary closure, which you should be able to understand by looking at its propagate()
function alongside the propagation look-up tables that are initialised in tables()
. While this periodic boundary closure of distribution functions \(f_i\) is native to the code, a set of hydrodynamic boundary conditions is also required to define a test flow. So-called ‘closure rules’ to capture these boundary conditions are implemented in a separate function to ‘overcome’ the default periodic closure at certain locations.
Our target ‘testbench’ flow is the computation of 2D flow within a square ‘lid-driven’ cavity (see Figure 1, left). This is defined by viscous hydrodynamic no-slip (Dirichlet) boundary conditions that require the fluid moves at the same velocity as any solid boundaries (i.e. the fluid sticks to the moving boundaries). The code implements no-slip boundary conditions using LBE’s infamous ‘bounce-back’ condition: populations propagating into the boundary wall are bounced back down the link along which they approached. This introduces dissipation at the boundary and is used in place of a standard collision for the wall. Our code has \(o(1)\) accurate on-link and \(o(1.5)\) accurate ‘mid-link’ variants to select from, which you can see working in the function DirichletBoundaryApply1()
.
It is easiest to understand the working of the bounce-back lattice closure conditions using the ‘on-link’ set. Considering the ‘on-link’ set, bounce-back is implemented along the lines
x=0
(with an image set appearing atx = XLENGTH
),y=0
(with an image set aty = YLENGTH
)
In addition, a ‘lid’ at y = YLENGTH_SYS
, with horizontal velocity \(U_0\) in the \(x\) direction is required. Figure 1 below shows how the lattice closure rules allow us to simulate the target lid-driven cavity in CCP5_01.c
.
Figure 1. Target hydrodynamic boundary conditions (left) and the corresponding LBE lattice closure rules (right) to be used throughout Exercise 1 to represent a ‘lid-driven cavity’ system. The black lines denoting the hydrodynamic boundary condition \(\mathbf{v} = \mathbf{0}\) (left) are represented by bounce-back closures in LBE simulations (right, green lines), while the red lines denote fluid assumed to be moving at velocity \(U_{0}\mathbf{\hat{e}}_{x}\). Note that the LBE domain (right) is over-sized - the top portion above the right-moving lid is simply discarded.
Accordingly, the square Eulerian mesh of this LBE simulation - the so-called flow domain - extends over the \(x, y\) region:
XLENGTH
\(> x \geq 0, ~~~~~\) YLENGTH
\(> y \geq 0\).
The lid is implemented crudely using so-called equilibrium forcing, where the distribution functions are overwritten by equilibrium values characteristic of the target motion to be induced, i.e. \(f_i = f_i^{(0)} \left(\rho, U_{0}\mathbf{\hat{e}}_{x}\right)\). To start with, we shall select on-link bounce-back for the other boundaries.
One last definition before we start: the Reynolds number is the most important example of a ‘dimensionless group’ in fluid mechanics. It is defined to be
where \(U =\) characteristic problem velocity [1], \(L =\) characteristic problem length and \(\nu =\) fluid kinematic viscosity. It is essentially the ratio of inertial forces to viscous forces and can be used to characterise various types of flow, e.g. viscosity-dominated laminar flows vs. Inertia-dominated turbulent flows.