LimB: A 2D nonequilibrium entropy limiter lattice Boltzmann code

Current version: 1.06 (June 2009)

Author: Robert Brownlee

Nonequilibrium entropy limiters for lattice Boltzmann methods (LBMs) are a novel technique for the stabilisation of LBMs. This innovation was pioneered and championed by R. A. Brownlee, A. N. Gorban and J. Levesley in the series of papers:

1. Nonequilibrium entropy limiters in lattice Boltzmann methods. Physica A, 387(2–3):385–406, 2008 (arXiv:0704.0043);
2. Stable simulation of fluid flows with high-Reynolds number using Ehrenfests' steps. Numer. Algorithms, 45(1–4):389–408, 2007 (arXiv:0705.4371);
3. Stability and stabilisation of the lattice Boltzmann method. Phys. Rev. E, 75(3):036711, 2007 (arXiv:cond-mat/0611444);
4. The Ehrenfests' coarse-graining as a paradigm for stabilising the lattice Boltzmann method. In Proc. Russian National Seminar "Modeling of Nonequilibrium Systems", Krasnoyarsk, October 13-15 2006, pages 3-13, Krasnoyarsk State Technical University Press, 2006;
5. Stabilisation of the lattice Boltzmann method using the Ehrenfests' coarse-graining idea. Phys. Rev. E, 74(3):037703, 2006 (arXiv:cond-mat/0605359).

The code used to test the concept of nonequilibrium entropy limiters in the above papers is available for download: LimB.zip

During the first 6 months of 2009, the code was extensively revisited, overhauled and improved by the code's author (R. A. Brownlee) and subsequently given the code name LimB. Perhaps considered a naive implementation by some, LimB will be of interest to those concerned with the development of nonequilibrium entropy limiters for LBMs as well as to newcomers who are developing their first LBM code. LimB is free software, but only made available under the proviso that appropriate acknowledgement is provided. Limited discretionary support is available by contacting the author. Freely distributable movies and figures produced using LimB are also available (see below).

LimB is written in Fortran 90 and its features are summarised in the following list:

• 2D isothermal lattice Boltzmann;
• D2Q9 lattice;
• solves for user defined initial conditions as well as 9 fundamental example problems:
  • Taylor-Green vortex;
  • Poiseuille flow;
  • Couette flow;
  • lid-driven cavity;
  • double-driven cavity;
  • 2D shock tube;
  • flow around a square cylinder;
  • backward facing step;
  • forward facing step.
• optional rectangular object in flow;
• choice of boundary condition:
  • periodic;
  • no-slip
    • bounce-back;
    • diffusive [Kinetic boundary conditions in the lattice Boltzmann method. Phys. Rev. E, 66(2):026311, 2002 (arXiv:nlin/0206033v1)].
  • free-slip
  • inlet/outlet
• choice of equilibrium:
  • polynomial;
  • entropy maximising.
• choice of nonequilibrium entropy limiter:
  • Ehrenfests' steps [Stabilisation of the lattice Boltzmann method using the Ehrenfests' coarse-graining idea. Phys. Rev. E, 74(3):037703, 2006 (arXiv:cond-mat/0605359)];
  • median filter [Nonequilibrium entropy limiters in lattice Boltzmann methods. Physica A, 387(2–3):385–406, 2008 (arXiv:0704.0043)].
• data is output in FORTRAN unformatted format.
• other code accomplishments:
  • main loop wholly unswitched;
  • main loop totally vectorised.
• tested working using both Intel FORTRAN and gfortran. However, not guaranteed to work on any compiler because compiler-dependent implementations of subroutine pointers are used.

The code has the following modular structure:

• LimB.f90 [main loop];
• LimBparam.f90 [parameter and variable declarations];
• LimBinout.f90 [collection of subroutines used to initialise LimB and control output];
• LimBprobs.f90 [subroutines used to initialise LimB for various problems];
• LimBsubs.f90 [various fundamental subroutines used in the LimB code (including propagate, collide, evolve)];
• LimBbcs.f90 [subroutines used in setting up boundary buffer for various conditions];
• LimBneq.f90 [subroutines for computing nonequilibrium entropy];
• LimBlims.f90 [subroutines for nonequilibrium entropy limiters];

The main loop is contained within the file LimB.f90 and the code is executed by submitting a parameter file LimBprob.par with the following structure:

• LimBprob [1=Taylor Green; 2=Poiseuille flow; 3=Couette flow; 4=Lid driven cavity, 5=Double driven cavity, 6=2D shock tube, 7=Square Cylinder, 8=BFS, 9=FFS];
• k1 [problem parameter (for LimBprob=1: wave number in x-direction; LimBprob=2: NA; LimBprob=3: NA; LimBprob=4/5: 0/1 constant/Leriche's poly profile; LimBprob=6: Density ratio of big to small region; LimBprob=7: 0/1 constant/parabolic inlet profile; LimBprob=8/9: Step height)];
• k2 [problem parameter (for LimBprob=1: wave number in y-direction; LimBprob=2: NA; LimBprob=3: NA; LimBprob=4/5: paramter in Leriche's poly profile); LimBprob=6: Fluid configuration parameter; LimBprob=7: NA; LimBprob=8/9: Step length)];
• BCflag [boundary condition flag (for LimBprob=0: NA; LimBprob=1: NA; LimBprob=2: 0=Periodic+BB,1=Periodic+Diff; LimBprob=3: NA; LimBprob=4/5: 0=Constant lid+BB,1=Diff; LimBprob=6: 0=BB,1=Diff; LimBprob=7: 0=BB,1=Diff feeslip; LimBprob=8: 0=BB,1=Diff)];
• nx [number of lattice sites in x-direction];
• ny [number of lattice sites in y-direction];
• Re [Reynolds number];
• rho0 [initial problem density];
• u_f [charactistic problem velocity];
• obj [flag for precence of rectangular object];
• cox [x-coordinate of lower-left corner of object];
• coy [y-coordinate of lower-left corner of object];
• ox [charactistic length of object in x-direction];
• oy [charactistic length of object in y-direction];
• oBCtype [flag for BC on object];
• tend [run simulation for tend dimensionless time units];
• eqflag [equilibrium flag (0 for entropy maximisers, 1 for polynomial equilbria)];
• limflag [limiter flag (0: no limiter; 1: Ehrenfests' (k,delta); 2: Median limiter; 2: Median limiter (k,delta))];
• limpar1 [k parameter in limiter];
• limpar2 [delta parameter in limiter];
• DSflag [nonequilibrium entropy type (0=based on perfect entropy; 1=Kullback entropy)];
• nframes [number of frames];
• ndiag [diagnostic info rate];
• outsig [not yet implimented: output signal flag];
• sigtype [not yet implimented: sigtype=0 means U is recorded; sigtype=1 means V is recorded];
• sx1 [not yet implimented: x-coord of signal point 1];
• sy1 [not yet implimented: y-coord of signal point 1];
• sx2 [not yet implimented: x-coord of signal point 2];
• sy2 [not yet implimented: y-coord of signal point 2];
• sx3 [not yet implimented: x-coord of signal point 3];
• sy3 [not yet implimented: y-coord of signal point 3];
• LimBrestart [restart code (by reloading f)];

Examples

Flow around a square cylinder

Unless otherwise stated, the following applies to all of the 2d square cylinder simulations below:

• characteristic length L=20 (side of cylinder);
• 601×501 uniformly spaced square lattice in used (301101 lattice sites in total);
• channel blockage ratio of 4%;
• free-stream velocity (u_∞,v_∞)=(0.05,0.0) (in lattice units);
• constant inlet condition, zeroth-order extrapolation at outlet;
• free-slip boundary condition on channel walls;
• Maxwell diffusive boundary condition on cylinder;
• simulation run for 125000 time steps;
• Ehrenfests' tolerances (k,δ)=(10,10^-3);
• only part of the computational domain is shown.

Square cylinder simulations:

• Re=37: vorticity movie (close to the onset of vortex shedding) [1m18s];
• Re=50: vorticity movie (vortex shedding has begun) [1m18s];
• Re=110: vorticity movie (flow still laminar, one pronounced frequency) [1m18s];
• Re=250: vorticity movie (transitional) [1m18s];
• Re=500: vorticity movie [1m18s];
• Re=1000: vorticity movie (onset of turbulence) (with Ehnrefests' steps histogram) [1m18s];
• Re=3000: vorticity movie [1m18s] (with Ehnrefests' steps histogram) [1m18s];
• Re=5000: vorticity movie [1m18s] (with Ehnrefests' steps histogram) [1m18s];
• Re=10000: vorticity movie [1m18s] (with Ehnrefests' steps histogram) [1m18s];
• Re=20000: vorticity movie [1m18s] (with Ehnrefests' steps histogram) [1m18s];
• Strouhal–Reynolds relationship;
• comparision of LBGK with LBGK-ES at Re=3000 for various values of the parameter pair (k,δ) [19s];
• comparision of LBGK with LBGK-ES at Re=3000 for various values of the parameter pair (k,δ) showing position of Ehrenfests' steps [19s].

Lid-driven cavity flow

The following applies to all of the lid-driven cavity simulations below:

• characteristic length L=320 (length of lid);
• 321×321 uniformly spaced square lattice in used (103041 lattice sites in total);
• fluid initially at rest;
• Maxwell diffusive boundary condition on all sides of the cavity;
• lid moving with constant velocity (u_∞,v_∞)=(0.075,0.0);
• simulation run for 500000 time steps;
• Ehrenfests' tolerances (k,δ)=(16,10^-3).

Lid-driven cavity simulations:

• Re=1000: stream-function movie [13s];
• Re=5000: stream-function movie [25s];
• Re=20000: vorticity movie [25s] / stream-function movie [25s].
• Re=200000: vorticity movie [62s].