MATH920: Lattice Boltzmann Equilibrium Distribution

Overview

• Boltzmann equation, conservation laws

• Hermite polynomial expansions

• Equilibrium distribution

Boltzmann equation

• $f(𝒙,𝝃,t)$ probability distribution function

• Force-free Boltzmann equation (dimensionless)

  ${\displaystyle \frac{\partial f}{\partial t}+{\xi }_{\alpha }\frac{\partial f}{\partial x_{\alpha }}=\Omega \left(f\right)}$

• Conservation laws

  ${\displaystyle \begin{array}{llll}\mathrm{Mass}& \int f(𝒙,𝝃,t)d𝝃=& \int f^{\mathrm{eq}}(𝒙,𝝃,t)d𝝃=& \rho (𝒙,t)\\ \mathrm{Momentum}& \int f(𝒙,𝝃,t)𝝃d𝝃=& \int f^{\mathrm{eq}}(𝒙,𝝃,t)𝝃d𝝃=& \rho 𝒖(𝒙,t)\\ \mathrm{Total}\mathrm{energy}& \frac{1}{2}\int f(𝒙,𝝃,t){\left|𝝃\right|}^{2}d𝝃=& \frac{1}{2}\int f^{\mathrm{eq}}(𝒙,𝝃,t){\left|𝝃\right|}^{2}d𝝃=& \rho E(𝒙,t)\\ \mathrm{Internal}\mathrm{energy}& \frac{1}{2}\int f(𝒙,𝝃,t){|𝝃-𝒖|}^{2}d𝝃=& \frac{1}{2}\int f^{\mathrm{eq}}(𝒙,𝝃,t){|𝝃-𝒖|}^{2}d𝝃=& \rho e(𝒙,t)\end{array}}$

Hermite Polynomials in 1D

• Definition from generating (weight) function

  ${\displaystyle \omega \left(x\right)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{x^2}{2}),H^{\left(n\right)}\left(x\right)={(-1)}^n\frac{1}{\omega \left(x\right)}\frac{d^n}{d{x}^n}\omega \left(x\right)}$

  In[1]:=

  omega[x_] = Exp[-x^2/2]/Sqrt[2Pi]

  ${\displaystyle \frac{e^{-\frac{x^2}{2}}}{\sqrt{2\pi }}}$

  In[2]:=

  H[n_,x_]:=Expand[(-1)^n/omega[x] D[omega[x],{x,n}]]; Table[H[n,x],{n,0,5}]

  ${\displaystyle \{1,x,x^2-1,x^3-3x,x^4-6x^2+3,x^5-10x^3+15x\}}$

  In[3]:=

  HHd1=Table[Integrate[omega[x] H[i,x] H[j,x],{x,-Infinity,Infinity}], {i,0,4},{j,0,4}] // MatrixForm

  ${\displaystyle \left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 2& 0& 0\\ 0& 0& 0& 6& 0\\ 0& 0& 0& 0& 24\end{array}\right)}$

  In[4]:=

Multidimensional Hermite polynomials

• Define multivariate generating function

  ${\displaystyle \omega \left(𝒙\right)=\frac{1}{{\left(2\pi \right)}^{d/2}}\exp (-\frac{{𝒙}^2}{2}),{𝑯}^{\left(n\right)}\left(𝒙\right)={(-1)}^n\frac{1}{\omega \left(𝒙\right)}{\nabla }^{\left(n\right)}\omega \left(𝒙\right),{\nabla }^{\left(n\right)}=\nabla \otimes \dots \otimes \nabla }$ ${\displaystyle {\nabla }^{\left(1\right)}=\nabla =\left(\begin{array}{l}\frac{\partial }{\partial x_1}\\ ⋮\\ \frac{\partial }{\partial x_d}\end{array}\right),{\nabla }^{\left(2\right)}=\nabla \otimes \nabla =\left(\begin{array}{lll}\frac{{\partial }^2}{\partial x_1}& \dots & \frac{{\partial }^2}{\partial x_1\partial x_d}\\ ⋮& \ddots & ⋮\\ \frac{{\partial }^2}{\partial x_d\partial x_1}& \dots & \frac{{\partial }^2}{\partial x_d\partial x_d}\end{array}\right)}$

  In[4]:=

  omega[x_,y_,z_] = Exp[-(x^2+y^2+z^2)/2]/(2Pi)^(3/2)

  ${\displaystyle \frac{e^{\frac{1}{2}(-x^2-y^2-z^2)}}{2\sqrt{2}{\pi }^{3/2}}}$

  In[5]:=

Hermite expansion of equilibrium distribution

• Seek Hermite expansion of equilibrium function ($d=3$)

  ${\displaystyle f^{\mathrm{eq}}(𝝃;\rho ,𝒖,\theta )=\frac{\rho }{{\left(2\pi \theta \right)}^{d/2}}\exp [-\frac{{(𝝃-𝒖)}^2}{2\theta }]=\omega \left(𝝃\right)\sum _{n=0}^{\infty }\frac{1}{n!}{𝒂}^{\left(n\right)}(\rho ,𝒖,\theta )\cdot {𝑯}^{\left(n\right)}\left(𝝃\right)}$ ${\displaystyle {𝒂}^{\left(n\right)}(\rho ,𝒖,\theta )=\int f^{\mathrm{eq}}(𝝃;\rho ,𝒖,\theta ){𝑯}^{\left(n\right)}\left(𝝃\right)d𝝃}$

• $f^{\mathrm{eq}}\left(𝝃\right)$ of same form as $\omega \left(𝝃\right)=\exp (-{𝒙}^2/2)/{\left(2\pi \right)}^{d/2}$

  ${\displaystyle f^{\mathrm{eq}}(𝝃;\rho ,𝒖,\theta )=\frac{\rho }{{\left(2\pi \theta \right)}^{d/2}}\exp [-\frac{{(𝝃-𝒖)}^2}{2\theta }]=\frac{\rho }{{\theta }^{d/2}}\omega \left(\frac{𝝃-𝒖}{\sqrt{\theta }}\right)}$

• Compute coefficients ${𝒂}^{\left(n\right)}(\rho ,𝒖,\theta )$, $𝜼=(𝝃-𝒖)/\sqrt{\theta }$ (scale to thermal velocity)

  ${\displaystyle {𝒂}^{\left(n\right),\mathrm{eq}}(\rho ,𝒖,\theta )=\frac{\rho }{{\theta }^{d/2}}\int \omega \left(\frac{𝝃-𝒖}{\sqrt{\theta }}\right){𝑯}^{\left(n\right)}\left(𝝃\right)d𝝃=\rho \int \omega \left(𝜼\right){𝑯}^{\left(n\right)}(\sqrt{\theta }𝜼+𝒖)d𝜼}$

Evaluate additional equilibrium distribution coefficients (fluxes)

In[12]:= 

AppendTo[Aeq,aeq[3]];
Aeq[[1]]

In[13]:= 

Aeq[[2]]

In[14]:= 

Aeq[[3]]

In[15]:= 

Aeq[[4]]

In[16]:= 

First three terms of Hermite expansion correspond to conserved quantities

Hermite coefficient recursion

• Verify ${𝒂}^{\left(n\right),\mathrm{eq}}={𝒂}^{\left(n\right),\mathrm{eq}}\otimes 𝒖$, $a^{\left(0\right),\mathrm{eq}}=\rho$

  In[29]:=

  V={u,v,w}; a0[1] = aeq[1] /. rho -> Subscript[rho,0]

  ${\displaystyle \{{\rho }_{0}u,{\rho }_{0}v,{\rho }_{0}w\}}$

  In[30]:=

  a0[2] = Outer[Times,a0[1],V]

  ${\displaystyle \left(\begin{array}{ccc}{\rho }_0{u}^2& {\rho }_{0}uv& {\rho }_{0}uw\\ {\rho }_{0}uv& {\rho }_0{v}^2& {\rho }_{0}vw\\ {\rho }_{0}uw& {\rho }_{0}vw& {\rho }_0{w}^2\end{array}\right)}$

  In[31]:=

  aeq[2] /. theta->1

  ${\displaystyle \{\rho u^2,\rho uv,\rho v^2,\rho uw,\rho vw,\rho w^2\}}$

  In[32]:=

Hermite coefficient recursion

• Third-order

  In[32]:=

  a0[3] = DeleteDuplicates[Flatten[Outer[Times,a0[2],V]]]

  ${\displaystyle \{{\rho }_0{u}^3,{\rho }_0{u}^{2}v,{\rho }_0{u}^{2}w,{\rho }_{0}u{v}^2,{\rho }_{0}uvw,{\rho }_{0}u{w}^2,{\rho }_0{v}^3,{\rho }_0{v}^{2}w,{\rho }_{0}v{w}^2,{\rho }_0{w}^3\}}$

  In[33]:=

  aeq[3] /. theta->1

  ${\displaystyle \{\rho u^3,\rho u^{2}v,\rho u{v}^2,\rho v^3,\rho u^{2}w,\rho uvw,\rho v^{2}w,\rho u{w}^2,\rho v{w}^2,\rho w^3\}}$

  In[34]:=

Hermite quadrature

• ${𝒂}^{\left(n\right),\mathrm{eq}}(\rho ,𝒖,\theta )=\int f^{\mathrm{eq}}\left(𝝃\right){𝑯}^{\left(n\right)}\left(𝝃\right)d𝝃$, $f^{\mathrm{eq},N}\left(𝝃\right)=\omega \left(𝝃\right)\rho Q(𝒖,\theta ,𝝃)\Rightarrow$

  ${\displaystyle {𝒂}^{\left(n\right),\mathrm{eq}}(\rho ,𝒖,\theta )=\rho \int \omega \left(𝝃\right)Q(𝒖,\theta ,𝝃){𝑯}^{\left(n\right)}\left(𝝃\right)d𝝃=\sum _{i=0}^{[(n+2)/2]}{w}_{i}Q(𝒖,\theta ,{𝝃}_i){𝑯}^{\left(n\right)}\left({𝝃}_i\right)}$

  In[34]:=

  H[2,x]

  ${\displaystyle x^2-1}$

  In[35]:=

Discrete equilibrium functions

• Consider a discrete velocity set $\{{𝝃}_0,\dots ,{𝝃}_q\}$. The Hermite expansion of the Maxwell-Boltzmann distribution up to second order $N=2$ is

  ${\displaystyle f^{\mathrm{eq}}(𝝃;\rho ,𝒖,\theta )\cong f^{\mathrm{eq},N}(𝝃;\rho ,𝒖,\theta )=\omega \left(𝝃\right)\rho Q(𝝃;𝒖,\theta )}$

• Evaluate the $N=2$ approximation for discrete velocity ${𝝃}_i$

  ${\displaystyle f_i^{\mathrm{eq}}=\rho w_{i}Q({𝝃}_i;𝒖,\theta ),w_i\equiv \omega \left({𝝃}_i\right)}$

  In[39]:=

  Q[xi_,eta_,zeta_]=Simplify[feqH[xi,eta,zeta]/omega[xi,eta,zeta]/rho]

  ${\displaystyle \frac{1}{2}({\eta }^2(\theta +v^2-1)+\eta v(u\xi +w\zeta +2)+\theta ({\xi }^2+{\zeta }^2-3)+u^2{\xi }^2-u^2+uw\xi \zeta +2u\xi -v^2+w^2{\zeta }^2-w^2+2w\zeta -{\xi }^2-{\zeta }^2+5)}$

  In[47]:=

  Q[xi,0,0] /. theta->1

  ${\displaystyle \frac{1}{2}(u^2{\xi }^2-u^2+2u\xi -v^2-w^2+2)}$

  In[48]:=

Discrete Boltzmann equation

• For a discrete set of velocities ${\left\{{𝒄}_i\right\}}_{i=0,1,\dots ,Q}$ the discrete Boltzmann equation is

  ${\displaystyle {\partial }_t{f}_i+c_{i\alpha }{\partial }_{\alpha }{f}_i=\Omega \left(f_i\right),i=0,1,\dots ,Q,}$ ${\displaystyle {\partial }_t{f}_i+({𝒄}_i\cdot \nabla )f_i=\Omega \left(f_i\right)}$

• Step 1 of operator splitting, ${\partial }_t{f}_i+({𝒄}_i\cdot \nabla )f_i=0$,

  ${\displaystyle f_i(t+\delta t,𝒙)=f_i(t,𝒙-{𝒄}_i\delta t)}$

• Step 2 of operator splitting

  ${\displaystyle {\partial }_t{f}_i=\Omega \left(f_i\right)}$