[注] 本文同步发表在 知乎.

宇宙安全声明:
1.这份文档来自我的一份课程作业, 由于笔者水平有限, 所以在内容上写得会比较粗糙, 或者出现一些错误.如有错误, 敬请谅解并指正.
2. 后续可能随缘更新?【如果又有空去读什么其他文章的话?】

Boltzmann方程和格子Boltzmann方法

Boltzmann方程

Boltzmann方程源自气体动理论, 它通过建模粒子系统的概率分布函数, 描述其集体行为.速度分布函数 f(x,ξ,t)f\left( \mathbf{x},\mathbf{\xi},t \right) 的含义通常表示为:在 tt 时刻位于 x\mathbf{x} 位置的球体微元 dV=[x,x+dx]dV = [\mathbf{x}, \mathbf{x}+ d\mathbf{x}] 内, 粒子速度位于区间 [ξ,ξ+dξ][\mathbf{\xi},\mathbf{\xi} + d\mathbf{\xi}] 的分子数量为 f(x,ξ,t)dξdVf\left( \mathbf{x},\mathbf{\xi},t \right)d\mathbf{\xi}dV .Boltzmann方程可以视作是 f(x,ξ,t)f\left( \mathbf{x},\mathbf{\xi},t \right) 的雷诺输运方程, 即:

DfDt=ft+ξfx+afξ=Ω(1.1)\frac{\mathrm{D}f}{\mathrm{D}t} = \frac{\partial f}{\partial t} + \mathbf{\xi} \cdot \frac{\partial f}{\partial\mathbf{x}} + \mathbf{a} \cdot \frac{\partial f}{\partial\mathbf{\xi}} = \Omega \tag{1.1}

其中 ξ=xt\mathbf{\xi =}\frac{\partial\mathbf{x}}{\partial t} 为粒子速度, a=ξt=Fm\mathbf{a =}\frac{\partial\mathbf{\xi}}{\partial t}\mathbf{=}\frac{\mathbf{F}}{m} 为粒子的外力加速度, F\mathbf{F} 为粒子所受外力, mm 为粒子质量, Ω\Omega 为描述分子碰撞作用的项.分布函数 ff 的速度矩可导出不同宏观量, 即.

n=fdξ,ρ=mn,ρu=mfξdξ,ρE=mξ22fdξ\begin{aligned} n = \int f d\mathbf{\xi} ,&\qquad \rho = m \cdot n,\\ \rho\mathbf{u} = m\int {f \mathbf{\xi}} d\mathbf{\xi} ,&\qquad \rho E = m\int {\frac{\xi^{2}}{2} f} d\mathbf{\xi} \end{aligned}

nn 为分子数密度, ρ\rho 为流体密度, u\mathbf{u} 为流体速度, EE为总能的密度.

对分子的碰撞进行建模并非易事.为简化式(1.1)的计算, Bhatnagar等提出了BGK碰撞项.这一思路主张将复杂的碰撞视作当前状态 ff 向Maxwell平衡态分布 f(eq)f^{(eq)} 演化的过程, 即:

Ω=1τ0(ff(eq))\Omega = - \frac{1}{\tau_{0}}\left( f - f^{(eq)} \right)

其中 τ0\tau_{0} 为松弛时间.DD 维空间下的Maxwell平衡态分布表示为:

f(eq)=n(2πRgT)D2exp[(ξu)22RgT].(1.2)f^{(eq)} = \dfrac{n}{\left( 2\pi R_g T \right)^{\frac{D}{2}}}\exp\left\lbrack - \frac{\left( \mathbf{\xi} - \mathbf{u} \right)^{2}}{2R_g T} \right]. \tag{1.2}

$R_g $ 为气体常数, TT 为温度.

在格子Boltzmann方法(Lattice Boltzmann Method, LBM)中, 使用的分布函数是密度分布函数 mfmf , 下文为方便起见会将密度分布函数记作 ff .密度分布函数的演化方程与式(1.1)一致, 其平衡态分布为:

f(eq)=ρ(2πRgT)D2exp[(ξu)22RgT].(1.3)f^{(eq)} = \frac{\rho}{\left( 2\pi R_g T \right)^{\frac{D}{2}}}\exp\left\lbrack - \frac{\left( \mathbf{\xi} - \mathbf{u} \right)^{2}}{2R_g T} \right]. \tag{1.3}

宏观量则表示为:

ρ=fdξ,j=ρu=fξdξ,ρE=ξ22fdξ.\rho = \int f d\mathbf{\xi},\quad \mathbf{j} = \rho\mathbf{u} = \int {f \mathbf{\xi}}d\mathbf{\xi} ,\quad \rho E = \int_{}^{}{\frac{\xi^{2}}{2}f}d\mathbf{\xi}.

格子Boltzmann方法

广义上说, LBM是离散Boltzmann方法(Discrete Boltzmann Method, DBM)的一种, 是式(1.1)的一种数值格式.本节介绍从式(1.1)到格子Boltzmann方程的离散思路, 推导中仅考虑外力加速度 a=F/m0\mathbf{a} =\mathbf{F}/ m \equiv 0 的情况.

对式(1.1)沿着特征线 ξ=x/t\mathbf{\xi} = \mathbf{x}/tΔt\Delta t 时间段内作积分, 略去 O(Δt2)O({\Delta t}^{2}) , 得

f(x+ξΔt,ξ,t+Δt)f(x,ξ,t)=1τ[f(x,ξ,t)f(eq)(x,ξ,t)](1.4)f( \mathbf{x} + \mathbf{\xi}\Delta t,\mathbf{\xi},t + \Delta t ) - f( \mathbf{x},\mathbf{\xi,}t) = -\frac{1}{\tau} \left\lbrack f(\mathbf{x},\mathbf{\xi},t) - f^{(eq)}(\mathbf{x},\mathbf{\xi},t) \right] \tag{1.4}

其中 τ=τ0/Δt\tau = \tau_{0}/\Delta t 为无量纲松弛时间.为了数值地计算这一连续方程, 需要将其在 ξ\mathbf{\xi} 的空间中进行离散, 并使其平衡态满足速度矩的约束, 即:

f(eq)(x,ξ,t)Ψ(ξ)dξ=iwiΨ(ξi)f(eq)(x,ξi,t)(1.5)\int f^{(eq)} \left( \mathbf{x},\mathbf{\xi}, t \right) \cdot \Psi\left( \mathbf{\xi} \right) d\mathbf{\xi} = \sum_{\mathbf{i}} {w_i\Psi\left( \mathbf{\xi}_i \right) f^{(eq)} \left( \mathbf{x},\mathbf{\xi}_i,t \right)} \tag{1.5}

这里 Ψ(ξ)\Psi\left( \mathbf{\xi} \right)ξ\mathbf{\xi} 的多项式, ξi\mathbf{\xi}_i 为离散速度, wiw_i 为积分的权系数.具体到Boltzmann方程中, 则至少应满足质量和动量守恒:

ρ=ifi=ifi(eq),j=ρu=iξifi=iξifi(eq)\rho = \sum_i f_i = \sum_i f_i^{(eq)},\qquad \mathbf{j} = \rho\mathbf{u} = \sum_i {\mathbf{\xi}_if_i} = \sum_i {\mathbf{\xi}_if_i^{(eq)}}

以及内能 ee 的方程:

ρe=12i(ξiu)2fi=12i(ξiu)2fi(eq)\rho e = \frac{1}{2}\sum_i^{}{\left( \mathbf{\xi}_i\mathbf{- u} \right)^{2}f_i} = \frac{1}{2}\sum_i^{}{\left( \mathbf{\xi}_i\mathbf{- u} \right)^{2}f_i^{(eq)}}

其中 fi=Wif(x,ξi,t)f_i = W_if\left( \mathbf{x},\mathbf{\xi}_i\mathbf{,}t \right)fi(eq)=Wif(eq)(x,ξi,t)f_i^{(eq)} = W_if^{(eq)}\left( \mathbf{x},\mathbf{\xi}_i\mathbf{,}t \right) .选定合适的离散速度集 {ci}(i=0,1,,b1)\{\mathbf{c}_i\}(i = 0,1,\ldots,b - 1) 及其对应的权重后, Boltzmann-BGK方程即可被数值计算.一般地, n 维空间下 b 速度的离散速度模型简称 DnQb 模型.

在低马赫数下, 平衡态可对 u\mathbf{u} 进行Taylor展开, 并截断至二阶:

f(eq)(x,ξ,t)=ρ(2πRgT)D2exp(ξ22RgT)exp(2(ξu)u22RgT)=ρexp(ξ22RgT)(2πRgT)D2[1+ξuRgT+(ξu)22(RgT)2u22RgT]+O(u3)(1.6)\begin{aligned} f^{(eq)}\left( \mathbf{x},\mathbf{\xi,}t \right) &= \frac{\rho}{\left( 2\pi R_g T \right)^{\frac{D}{2}}} \cdot \exp\left( - \frac{\mathbf{\xi}^{2}}{2R_g T} \right) \cdot \exp\left( \frac{2\left( \mathbf{\xi \cdot u} \right) - \mathbf{u}^{2}}{2R_g T} \right) \\ &= \frac{\rho\exp\left( - \frac{\mathbf{\xi}^{2}}{2R_g T} \right)}{\left( 2\pi R_g T \right)^{\frac{D}{2}}} \cdot \left[ 1 + \frac{\mathbf{\xi} \cdot \mathrm{u}}{R_g T} + \frac{\left( \mathbf{\xi \cdot u} \right)^{2}}{2 \left( R_g T \right)^{2}} - \frac{\mathbf{u}^{2}}{2R_g T} \right] + O\left( \mathbf{u}^{3} \right) \end{aligned} \tag{1.6}

因此, 当使用 截断的平衡态 计算 f(eq)Ψ(ξ)dξ\int {f^{(eq)} \cdot \Psi\left( \mathbf{\xi} \right)}d\mathbf{\xi} 时, 会导出形如 exp(x2)Ψ(x)dx\int {\exp\left( - x^{2} \right) \cdot \Psi(x)}dx 的积分, 它可以用高斯积分进行数值计算.根据上式的具体情况, 记

Im=+exp(ζ2)ζmdζI_{m} = \int_{- \infty}^{+ \infty}{\exp\left( - \zeta^{2} \right) \cdot \zeta^{m}} d\zeta

其中 ζ=ξ/2RgT\zeta = \xi/\sqrt{2 R_g T} .ImI_{m} 是权函数 exp(ζ2)\exp\left(-\zeta^{2}\right) 在某一实轴上的 mm 阶矩.

这里以D2Q9离散速度集为例介绍权重的计算, 其离散速度集为

cα={(0,0),α=0c(cosπ(α1)2,sinπ(α1)2),α=142c(cosπ(2α1)4,sinπ(2α1)4),α=58(1.7)\mathbf{c}_{\alpha} = \begin{cases} (0,0) & ,\alpha = 0 \\ c\left( \cos\frac{\pi(\alpha - 1)}{2},\sin\frac{\pi(\alpha - 1)}{2} \right) & ,\alpha = 1\ldots 4 \\ \sqrt{2}c\left( \cos\frac{\pi(2\alpha - 1)}{4},\sin\frac{\pi(2\alpha - 1)}{4} \right) & , \alpha = 5 \ldots 8 \\ \end{cases} \tag{1.7}

其中 c=Δx/Δtc = \Delta x/\Delta t .二维情况下的 Ψ\Psi 可定义为 Ψm,n(ξ)=ξxmξyn\Psi_{m,n}\left( \mathbf{\xi} \right) = \xi_{x}^{m}\xi_{y}^{n}, 则式(1.6)被展开为[1]

I=f(eq)Ψm,n(ξ)dξ=ρπ(2RgT)m+n2[(1u22RgT)ImIn+22RgT(uxIm+1In+uyImIn+1) +1RgT(ux2Im+2In+2uxuyIm+1In+1+uy2ImIn+2)](1.8)\begin{aligned} I & = \int f^{(eq)} \cdot \Psi_{m,n} \left( \mathbf{\xi} \right) d\mathbf{\xi} \\ &= \frac{\rho}{\pi} \left( 2 R_g T \right)^{\frac{m + n}{2}} \left[ \left( 1 - \frac{u^{2}}{2R_{g}T} \right)I_{m}I_{n} \right. \\ & + \frac{2}{\sqrt{2R_{g}T}}\left( u_{x}I_{m + 1}I_{n} + u_{y}I_{m}I_{n + 1} \right) \\ & \left. \ + \frac{1}{R_{g}T}\left( u_{x}^{2}I_{m + 2}I_{n} + 2u_{x}u_{y}I_{m + 1}I_{n + 1} + u_{y}^{2}I_{m}I_{n + 2} \right) \right] \end{aligned} \tag{1.8}

对D2Q9模型, ImI_{m} 可被替换为三点Gauss-Hermite数值积分:

Im=+exp(ζ2)ζmdζ=α=13ωαζαmI_{m} = \int_{- \infty}^{+ \infty}{\exp\left( - \zeta^{2} \right) \cdot \zeta^{m}}d\zeta = \sum_{\alpha = 1}^{3}{\omega_{\alpha}\zeta_{\alpha}^{m}}

[NOTE]
积分点 ζi\zeta_i 是Hermite多项式 H3(x)=8x312xH_{3}(x) = 8x^{3} - 12x 的零点, 对应权重 ωi=2313!H2(ζi)232π\omega_i = \frac{2^{3 - 1} \cdot 3!}{H_{2}\left( \zeta_i \right)^{2} \cdot 3^{2}} \cdot \sqrt{\pi}.
所以积分点和权重为:ζ1=3/2,w1=π/6\zeta_{1} = - \sqrt{3/2},w_{1} = \sqrt{\pi}/6ζ2=0,w2=2π/3\zeta_{2} = 0,w_{2} = 2\sqrt{\pi}/3ζ3=3/2,ω3=π/6\zeta_{3} = \sqrt{3/2},\omega_{3} = \sqrt{\pi}/6.

则式(1.8)可被数值离散为 [1:1]

I=ρi,j=13Ψ(ξi,j)ωiωjπ[1+ξi,juRgT+(ξi,ju)22(RgT)2u22RgT](1.9)I = \rho\sum_{i,j = 1}^{3} \Psi(\mathbf{\xi}_{i,j}) \cdot \frac{\omega_i \omega_j}{\pi} \cdot \left[ 1 + \frac{\mathbf{\xi}_{i,j} \cdot \mathbf{u}}{R_g T} + \frac{\left( \mathbf{\xi}_{i,j} \cdot \mathbf{u} \right)^{2}}{2 (R_g T)^{2}} - \frac{\mathbf{u}^{2}}{2 R_g T} \right] \tag{1.9}

此处 ξi,j=(ξi,ξj)=2RgT(ζi,ζj)\mathbf{\xi}_{i,j}\mathbf{=}\left( \xi_i\mathbf{,}\xi_{j} \right)\mathbf{=}\sqrt{2R_g T}\left( \zeta_i\mathbf{,}\zeta_{j} \right) .也就是说, 若令 RgT=cs2=c2/3R_g T = c_{s}^{2} = c^{2}/3 , 则 ξi,j\mathbf{\xi}_{i,j} 便对应D2Q9模型的 cα\mathbf{c}_{\alpha} .根据式(1.9), cα\mathbf{c}_{\alpha} 方向的平衡态分布 fα(eq)f_{\alpha}^{(eq)} 定义为

fα(eq)=ρwα[1+cαucs2+(cαu)22cs4u22cs2](1.10)f_{\alpha}^{(eq)} = \rho w_{\alpha} \cdot \left\lbrack 1 + \frac{\mathbf{c}_{\alpha} \cdot \mathbf{u}}{c_{s}^{2}} + \frac{\left( \mathbf{c}_{\alpha} \cdot \mathbf{u} \right)^{2}}{2c_{s}^{4}} - \frac{\mathbf{u}^{2}}{2c_{s}^{2}} \right\rbrack \tag{1.10}

其中

wα={4/9,i=j=2,α=01/9,(i=2)(j=2),α=1,,41/36,(i2)(j2),α=5,,8w_{\alpha} = \begin{cases} 4/9 ,& i = j = 2 ,& \alpha = 0 \\ 1/9 ,& (i = 2)\bigoplus(j = 2) ,& \alpha = 1, \ldots, 4 \\ 1/36 ,& (i \neq 2)\bigwedge(j \neq 2) ,& \alpha = 5, \ldots, 8 \\ \end{cases}

Qian等[2]已给出常用的DnQb模型参数, 大量学者也给出了其他离散速度模型的推导方式 [1:2][3][4].

在前文中, RgTR_g T 为常数一直是LBGK方程的推导前提, 这使得LBGK模型的最佳适用范围被限制在等温流动或无温流动.LBGK在模拟不可压缩流体时, 格子系统的马赫数应小于 0.15 [5].而LBGK框架虽然可通过部分手段实现亚音速[6]和高雷诺数[7]的模拟, 但仍旧无法直接解决热流问题.

就目前来说, 热流体的LBM模拟可以分为多速LBE模型 [8] [9] [10] [11] [12] [13] [14] 、双分布函数LBE模型 [15] [16] [17] [18] 、基于简化LBM的实现[19]、热流格子Boltzmann通量求解器[20]、多松弛模型[21] [22] [23] 以及与其他数值方法结合的混合模型等多种形式 [24] .

多速LBE模型

多速LBE模型(即Multi-Speed LBE)是一种完全基于LBM经典算法的改进.这类方法演化方程与常规LBM一致, 在保持算法不变的基础上通过使用比常规LBM更为复杂的离散速度模型还原高阶速度矩.它因只使用单一分布函数,物理描述上更接近事实而受到了广泛关注.并且使用更多的离散速度达到高阶精度的思路亦已经被广泛拓展到各种各样的LBM模拟中. 虽然Qian[11:1]和Alexander等[12:1]才是早期的提出者, 但本章首先以Chen等 [8:1] 的高阶MS-LBGK模型的二维形式为例引入多速LBM的概念, 随后介绍几种不同思路下的MS-LBM改进形式.

高阶MS-LBGK模型

Chen等 [8:2] 提出的多速LBM的平衡态分布函数定义为

fσi(eq)=Aσ+Bσ(cσiu)+Cσu2+Dσ(cσiu)2+Eσu2(cσiu)+Fσ(cσiu)3+Gσu2(cσiu)2+Hσu4(2.1)\begin{aligned} f_{\sigma i}^{(eq)} = & A_{\sigma} + B_{\sigma} \cdot \left( \mathbf{c}_{\sigma i} \cdot \mathbf{u} \right) + C_{\sigma}\left\| \mathbf{u} \right\|^{2} + D_{\sigma} \cdot \left( \mathbf{c}_{\sigma i} \cdot \mathbf{u} \right)^{2} \\ & + E_{\sigma}\left\| \mathbf{u} \right\|^{2} \cdot \left( \mathbf{c}_{\sigma i} \cdot \mathbf{u} \right) + F_{\sigma} \cdot \left( \mathbf{c}_{\sigma i} \cdot \mathbf{u} \right)^{3} \\ & + G_{\sigma}\left\| \mathbf{u} \right\|^{2} \cdot \left( \mathbf{c}_{\sigma i} \cdot \mathbf{u} \right)^{2} + H_{\sigma}\left\| \mathbf{u} \right\|^{4} \\ \end{aligned} \tag{2.1}

其中 AσA_{\sigma}HσH_{\sigma} 均为待定参数, 它们均可被表示为内能 ee 和密度 ρ\rho 的复合函数, 即 Xσ=(xσ0+xσ1e+xσ2e2)ρX_{\sigma} = \left( x_{\sigma 0} + x_{\sigma 1}e + x_{\sigma 2}e^{2} \right)\rho .并且 fσi(eq)f_{\sigma i}^{(eq)} 的各阶速度矩需要满足下列约束条件:

σifσi(eq)=ρ,σicσifσi(eq)=ρu,σicσicσifσi(eq)=ρuu+pI\sum_{\sigma i} f_{\sigma i}^{(eq)} = \rho, \qquad \sum_{\sigma i} {\mathbf{c}_{\sigma i}f_{\sigma i}^{(eq)}} = \rho\mathbf{u}, \qquad \sum_{\sigma i} {\mathbf{c}_{\sigma i}\mathbf{c}_{\sigma i}f_{\sigma i}^{(eq)}} = \rho\mathbf{uu} + p\mathbf{I}

σicσicσicσifσi(eq)=ρuuu+p [uδ]\sum_{\sigma i} {\mathbf{c}_{\sigma i}\mathbf{c}_{\sigma i}\mathbf{c}_{\sigma i}f_{\sigma i}^{(eq)}} = \rho\mathbf{uuu} + p\ [\mathbf{u \delta} ]

σicσicσicσi2fσi(eq)=ρu2uu+p(D+4)uu+pu2I+2(D+2)peD\sum_{\sigma i} {\mathbf{c}_{\sigma i}\mathbf{c}_{\sigma i} \vert{c_{\sigma i}}\vert ^{2}f_{\sigma i}^{(eq)}} = \rho |u|^{2} \mathbf{uu}\mathbf{+}p(D + 4)\mathbf{uu}\mathbf{+}p|u|^{2}\mathbf{I}\mathbf{+}\frac{2(D + 2)pe}{D}

其中压强 p=2ρeDp = \frac{2\rho e}{D} , 三阶张量 [uδ]= uαδβγ+uβδαγ+uγδαβ\left\lbrack \mathbf{u\delta} \right\rbrack = \ u_{\alpha}\delta_{\beta\gamma} + u_{\beta}\delta_{\alpha\gamma} + u_{\gamma}\delta_{\alpha\beta} .

为满足这五个守恒方程的约束, 张量 icσicσi...cσin\displaystyle \sum_i \overset{n\text{个}}{\overbrace{\mathbf{c}_{\sigma i}\mathbf{c}_{\sigma i}...\mathbf{c}_{\sigma i}}}\qquad 需要在 0~6 阶均为各向同性.因此Chen等 [8:3] 使用的离散速度集表示为

cσi=cpki=Perm{k(±1,,±1p,0,...,0Dp)}\mathbf{c}_{\sigma i} = \mathbf{c}_{pki}^{'} = \mathbf{Perm}\left\{ k\left( \underbrace{\pm 1,\ldots, \pm 1}_{p} , \overset{D - p}{\overbrace{0,...,0}} \right) \right\}

这里 Perm\mathbf{Perm} 是由 ()(*) 的所有排列构成的集合, kk 为格子的缩放系数, 且 σ=cσi2=k2p\sigma = c_{\sigma i}^{2} = k^{2}p. 因此, 高阶MS-LBGK模型的运动粘度为 μ=2ρeD(τ12)\mu = \frac{2\rho e}{D}\left( \tau - \frac{1}{2} \right) , 体积粘度 λ=2μD\lambda = - \frac{2\mu}{D} , 热扩散系数 κ=(D+2)ρeD(τ12)\kappa = \frac{(D + 2)\rho e}{D}\left( \tau - \frac{1}{2} \right) .对于气体常数为 $R_g $ 的单原子分子, 定压和定容比热容分别为 cp=(D+2)Rg/2c_{p} = (D + 2)R_g /2cv=DRg/2c_{v} = D R_g /2 .由于无量纲格子系统中 Rg=1R_g = 1 , 该模型的Prandtl数为 Pr=μcp/κ=1Pr = \mu c_{p}/\kappa = 1 .这也是MS-LBGK模型这类建模方式的常见问题.

目前, 学界已有数类方式在多速LBE模型中实现对Prandtl数的自由调节, 如:双松弛时间方法[9:1];引入熵格式[10:1] [13:1] [14:1] ;修改网格速度的定义 [3:1] [25] [26] 等.由于已有文献指出若不局限于常规晶格模型则可获得更高精度的离散格式, 并且大多数方法都涉及对DnQb模型的修改, 因此这里仅简要介绍前两种方式的实现思路.

对多速LBE模型的改进

基于双松弛时间的实现

基于双松弛时间的方法将碰撞分为由 fif_ifif_{-i} 控制的两个部分, 并分别进行松弛.以Chen等[9:2]提出的双松弛模型为例, 其演化方程为(不考虑外力项):

fi(x+ciΔt,t+Δt)fi(x,t)=Ωi+Ωi(2.2)f_i\left( \mathbf{x} + \mathbf{c}_i\Delta t,t + \Delta t \right) - f_i\left( \mathbf{x},t \right) = \Omega_i + \Omega_{\mathbb{i}} \tag{2.2}

其中

Ωi=1τ1(fifi(eq)),Ωi=1τ2(fifi(eq))\Omega_i = - \frac{1}{\tau_{1}}\left( f_i - f_i^{(eq)} \right),\qquad \Omega_{\mathbb{i}} = - \frac{1}{\tau_{2}}\left( f_{- i} - f_{- i}^{(eq)} \right)

下标 i-i 表示 ci=ci\mathbf{c}_{- i} = - \mathbf{c}_i 方向.τv=τ1τ2/(τ1+τ2)\tau_{v} = \tau_{1}\tau_{2}/(\tau_{1} + \tau_{2}) , τk=τ1τ2/(τ2τ1)\tau_{k} = \tau_{1}\tau_{2}/(\tau_{2} - \tau_{1}) .Chen等[9:3]指出其导出的宏观方程组为

DρDt=ρu(ρu)t+(ρuu)=p+Tv(ρe)t+(ρeu)=pu+(κT)+Tk:u\begin{aligned} \frac{D\rho}{Dt} &= \rho\nabla \cdot \mathbf{u} \\ \frac{\partial\left( \rho\mathbf{u} \right)}{\partial t} + \nabla \cdot \left( \rho\mathbf{uu} \right) &= - \nabla p + \nabla \cdot \mathbf{T}_{v} \\ \frac{\partial(\rho e)}{\partial t} + \nabla \cdot \left( \rho e\mathbf{u} \right) &= p\nabla \cdot \mathbf{u} + \nabla \cdot (\kappa\nabla T) + \mathbf{T}_{k}:\nabla\mathbf{u} \end{aligned}

其中

Tv=2μvS+λv(u)I,Tk=2μkS+λk(u)I\mathbf{T}_{v} = 2\mu_{v}\mathbf{S} + \lambda_{v}\left( \nabla \cdot \mathbf{u} \right)\mathbf{I},\mathbf{\qquad }\mathbf{T}_{k} = 2\mu_{k}\mathbf{S} + \lambda_{k}\left( \nabla \cdot \mathbf{u} \right)\mathbf{I}

S\mathbf{S} 为应变率张量.各输运系数为

μv=2Dρe(τv12),λv=4D2ρe(τv12)μk=2Dρe(τk12),λk=4D2ρe(τk12)\begin{aligned} \mu_{v} = \frac{2}{D}\rho e\left( \tau_{v} - \frac{1}{2} \right),& \lambda_{v} = \frac{- 4}{D^{2}}\rho e\left( \tau_{v} - \frac{1}{2} \right) \\ \mu_{k} = \frac{2}{D}\rho e\left( \tau_{k} - \frac{1}{2} \right),& \lambda_{k} = \frac{- 4}{D^{2}}\rho e\left( \tau_{k} - \frac{1}{2} \right) \end{aligned}

温度 T=De/2T = De/2 , 热传导系数 κ=ρe(D+2)(τk1/2)/D\kappa = \rho e(D + 2)\left( \tau_{k} - 1/2 \right)/D. 综上所述, 该模型的Prandtl数为

Pr=μvcpκ=2τv12τk1Pr = \frac{\mu_{v}c_{p}}{\kappa} = \frac{2\tau_{v} - 1}{2\tau_{k} - 1}

这种方法通过使用两个松弛时间 τ1,τ2\tau_{1},\tau_{2} 分别对动量和能量方程进行松弛(τ1,τ2>12\tau_{1},\tau_{2} > \frac{1}{2}), 缓解LBGK单个松弛时间的参数依赖性.但由于能量方程的粘性应力项 TkTv\mathbf{T}_{k} \neq \mathbf{T}_{v} , 因而其能量方程是与动量方程存在冲突的.

基于熵函数的实现

熵格子Boltzmann方程(Entropic Lattice Boltzmann Equation, ELBE)是另一种基于传统LBE的修改.这里以单松弛的ELBE为例, 它将LBE的演化方程修改为

fi(x+ciΔt,t+Δt)fi(x,t)=αβ(fifi(eq))(2.3)f_i\left( \mathbf{x} + \mathbf{c}_i\Delta t,t + \Delta t \right) - f_i\left( \mathbf{x},t \right) = - \alpha\beta\left( f_i - f_i^{(eq)} \right) \tag{2.3}

其中 β=1/(1+2ν/cs2)\beta = 1/(1 + 2\nu/c_{s}^{2}) , ν\nu 为运动粘度.α\alpha 为方程 H(f+α(ff(eq)))=H(f)H\left( f + \alpha\left( f - f^{(eq)} \right) \right) = H(f) 的解, 熵函数 HH 定义为[27]

H=ifiln(fiwi)(2.4)H = \sum_i^{}{f_i\ln\left( \frac{f_i}{w_i} \right)} \tag{2.4}

α\alpha 的解可通过迭代法或近似公式计算[28].由于部分常规DnQb模型在使用 f(eq)f^{(eq)} 的低阶展开(如式(1.10))时无法满足熵方程[29], 因此该方法通常需要对DnQb模型进行修改, 或搭配多速DnQb模型一同使用.

Pransianakis等[10:2] [14:2] 在ELBE的总体框架上, 将碰撞项视为从 fif_i 到中间态 fif_i^{*} 再到平衡态 fi(eq)f_i^{(eq)} 的两步松弛, 即:

1τ1(fifi)1τ2(fifi(eq))\frac{-1}{\tau_{1}}\left( f_i - f_i^{*} \right) - \frac{1}{\tau_{2}}\left( f_i^{*} - f_i^{(eq)} \right)

中间态 fif_i^{*} 是平衡态的扰动, 即:fi=fi(eq)+δfif_i^{*} = f_i^{(eq)} + \delta f_i^{*}.校正量δfi\delta f_i^{*} 的计算可基于 “中间态的热通量不变” 假设构造 δfi\delta f_i^{*} 的方程组进行求解[10:3]. 松弛时间为 τ1=μ/ρT0\tau_{1} = \mu/\rho T_0τ2=2κ/ρT0\tau_{2} = 2\kappa/\rho T_0 , 分别控制应力张量和热流量的松弛, T0T_0 是格子系统中的参考温度.因此, Prandtl数为 Pr=4τ1/τ2\Pr = 4\tau_{1}/\tau_{2} .

[NOTE]
需要强调的是, 在Pransianakis等[10:4]在文中使用的D2Q9模型里, 虽然特征速度方向与式(1.7)一致, 但各个ci\mathbf{c}_i的权重 WiW_i 是温度TT的函数, 即: Wi=(1T2)T2(1T)ci2\displaystyle W_i = \frac{\left( 1 - T^{2} \right)T}{2(1 - T)}c_i^{2}

而平衡态分布写作: fi(eq)=ρWi{1+cijρT+jjT2(ρT)2:[ciciT4T2+ci2(13T)2(1T)I]}\displaystyle f_i^{(eq)} = \rho W_i\left\{ 1 + \frac{\mathbf{c}_i \cdot \mathbf{j}}{\rho T} + \frac{\mathbf{j}\mathbf{j}^{T}}{2(\rho T)^{2}}:\left\lbrack \mathbf{c}_i\mathbf{c}_i^{T} - \frac{4T^{2} + c_i^{2}(1 - 3T)}{2(1 - T)}\mathbf{I} \right\rbrack \right\}

Frapolli等[13:2]提出了将多速LBE模型与熵LBE框架进行结合的思路.以D2Q25-ZOT模型为例, 其演化方程为:

fi(x+ciΔt,t+Δt)fi(x,t)=ω1(fifi)+ω2(fifi(eq))(2.5)f_i\left( \mathbf{x} + \mathbf{c}_i\Delta t,t + \Delta t \right) - f_i\left( \mathbf{x},t \right) = \omega_{1}\left( f_i - f_i^{*} \right) + \omega_{2}\left( f_i^{*} - f_i^{(eq)} \right) \tag{2.5}

在不考虑ELBE的校正时, 动力粘度 μ\mu 和热扩散系数 κ\kappa 的表达式为:

μ={(1ω112)ρT,Pr1(1ω212)ρT,Pr1,κ={(1ω212)ρTcp,Pr1(1ω112)ρTcp,Pr1\mu = \begin{cases} \left( \frac{1}{\omega_{1}} - \frac{1}{2} \right)\rho T,& \Pr \leq 1 \\ \left( \frac{1}{\omega_{2}} - \frac{1}{2} \right)\rho T,& \Pr \geq 1 \end{cases}, \qquad \kappa = \begin{cases} \left( \frac{1}{\omega_{2}} - \frac{1}{2} \right)\rho Tc_{p},& \Pr \leq 1 \\ \left( \frac{1}{\omega_{1}} - \frac{1}{2} \right)\rho Tc_{p},& \Pr \geq 1 \end{cases}

Frapolli等[13:3]认为, 上述计算中仅需要对涉及 μ\mu 的部分进行校正.当 Pr<1\Pr < 1 时, ω1=αβ1,ω2=β2\omega_{1} = \alpha\beta_{1},\omega_{2} = \beta_{2} ;当 Pr>1\Pr > 1 时, ω1=β1,ω2=αβ2\omega_{1} = \beta_{1},\omega_{2} = \alpha\beta_{2} . α\alpha 仍为方程 H(f+α(ff(eq)))=H(f)H\left( f + \alpha\left( f - f^{(eq)} \right) \right) = H(f) 的解.因此, μ\muκ\kappa 的计算变为:

μ={12(1β11)ρT,Pr112(1β21)ρT,Pr1,κ={(1β212)ρT,Pr1(1β112)ρT,Pr1 .\mu = \left\{ \begin{matrix} \frac{1}{2}\left( \frac{1}{\beta_{1}} - 1 \right)\rho T,\qquad \Pr \leq 1 \\ \frac{1}{2}\left( \frac{1}{\beta_{2}} - 1 \right)\rho T,\qquad \Pr \geq 1 \\ \end{matrix} \right., \qquad \kappa = \left\{ \begin{matrix} \left( \frac{1}{\beta_{2}} - \frac{1}{2} \right)\rho T,\qquad \Pr \leq 1 \\ \left( \frac{1}{\beta_{1}} - \frac{1}{2} \right)\rho T,\qquad \Pr \geq 1 \\ \end{matrix} \right.\ .

[NOTE]
Pr=1\Pr = 1 时模型退化为LBGK方程, 这里不做讨论.

上述所介绍的基于ELBE的热流模拟能够成功还原Fourier-Navier-Stokes方程组, 实现低马赫数下各种热流的数值模拟, 但也没有实现热流和流场的耦合.此外, 由于 HH 函数的引入, 这类方法需要额外求解 H(f+α(ff(eq)))=H(f)H\left( f + \alpha\left( f - f^{(eq)} \right) \right) = H(f) , 在数值计算上导致一定不便.

双分布函数LBE模型

双分布函数(double distribution function, DDF)的LBE模型的主要思路是将速度场和温度场分别用两套LBE进行计算.DDF-LBE下的热流分为两类:其一是将温度视作被动标量, 这要求温度变化对流动的影响较小;另一类考虑的温度对流场的影响(如可压缩热流).

不考虑黏性热耗散和压缩功的DDF-LBE

基本框架

Bartoloni等 [15:1] 指出:在压力做功和黏性热耗散项均可被忽略的场景中, 可将流场和温度场分别采用两套分布函数体系进行对流扩散.一般来说, 无源项的DDF-LBE演化方程组为:

fi(x+ciΔt,t+Δt)fi(x,t)=1τ1(fifi(eq))(3.1)f_i\left( \mathbf{x} + \mathbf{c}_i\Delta t,t + \Delta t \right) - f_i\left( \mathbf{x},t \right) = - \frac{1}{\tau_{1}}\left( f_i - f_i^{(eq)} \right)\tag{3.1}

gi(x+ciΔt,t+Δt)gi(x,t)=1τ2(gigi(eq))(3.2)g_i\left( \mathbf{x} + \mathbf{c}_i\Delta t,t + \Delta t \right) - g_i\left( \mathbf{x},t \right) = - \frac{1}{\tau_{2}}\left( g_i - g_i^{(eq)} \right)\tag{3.2}

其中 fif_igig_i 分别用于计算速度和温度, 并通过如下速度矩还原宏观物理量

ρ=ifi,j=ρu=icifi,ρT=igi(3.3)\rho = \sum_i f_i, \qquad \mathbf{j} = \rho\mathbf{u} = \sum_i {\mathbf{c}_if_i}, \qquad \rho T = \sum_i g_i \tag{3.3}

需要强调的是, fif_igig_i 分别使用不同的DnQb模型, 因此它们的格子声速分别为 cs1c_{s1}cs2c_{s2} .两个松弛时间( τ1,τ2\tau_{1},\tau_{2} )和运动粘度 ν\nu 、热扩散系数 κ\kappa 的关系为:

τ1=νcs12+12,τ2=κcs22+12(3.4)\tau_{1} = \frac{\nu}{c_{s1}^{2}} + \frac{1}{2},\qquad \tau_{2} = \frac{\kappa}{c_{s2}^{2}} + \frac{1}{2} \tag{3.4}

这种思路仅额外添加了一组分布函数及其平衡态的定义, 并且由于其易于实现, 目前是低速小Mach数场景下热流的一种主流LBM计算方法.

含黏性热耗散和压缩功的DDF-LBE

前面介绍的DDF-LBE都是在忽略黏性热耗散和压缩功的前提下, 单独另开一个分布函数 gig_i 来计算温度 TT 这一标量的变化.虽然基于温度的建模方式足够简单, 但其应用场景限制较大.为更全面地计算低马赫数热流中流场和温度的相互作用, 学界已经提出了一些方案.

基于能量的DDF-LBE

既然可以为温度TT构建一套LBE, 那么为其他可以表示能量的宏观量(如内能 ρe\rho e [16:1] 、总能 ρE=ρe+ρu2/2\rho E = \rho e + \rho u^{2}/2 [17:1] 、总焓 h=e+p/ρh = e + p/\rho [18:1] )构建一套LBE也是可行的.通过模拟能量场的演化, 更准确地描述流场和温度场之间的作用.由于这套思路的建模方式总体上较为相似, 下面仅以基于总能ρE\rho E的DDF-LBE为例对这种思路进行介绍.

单位体积内的总能可以由分布函数的速度矩表示为

ρE=ρe+ρu22=ρcvT+ρu22=12ξ2fdξ(3.5)\rho E = \rho e + \frac{\rho u^{2}}{2} = \rho c_{v}T + \frac{\rho u^{2}}{2} = \int_{}^{}{\frac{1}{2}\xi^{2}f}d\mathbf{\xi} \tag{3.5}

因此可以将总能的分布函数定义为 h=12ξ2fh = \frac{1}{2}\xi^{2}f , 则 ρE=hdξ\rho E = \int h d\mathbf{\xi} .

根据 ff 的演化方程(即式(1.1)), 可得到 hh 的演化方程为

ht+ξf+aξhfξa=Ωh(3.6)\frac{\partial h}{\partial t} + \mathbf{\xi} \cdot \nabla f + a \cdot \nabla_{\mathbf{\xi}}h - f\mathbf{\xi \cdot a} = \Omega_{h} \tag{3.6}

其中 Ωh=12ξ2Ωf\Omega_{h} = \frac{1}{2}\xi^{2}\Omega_{f} 为总能的碰撞算子.Guo等 [17:2]Ωh\Omega_{h} 视作机械能和内能两部分作用的组合, 分别采用 τf\tau_{f}τh\tau_{h} 控制, Ωh\Omega_{h} 则表示为

Ωh=1τh(hh(eq))+Zτhf(ff(eq))(3.7)\Omega_{h} = - \frac{1}{\tau_{h}}\left( h - h^{(eq)} \right) + \frac{Z}{\tau_{hf}}\left( f - f^{(eq)} \right) \tag{3.7}

其中 Z=ξuu22Z = \mathbf{\xi \cdot u} - \frac{u^{2}}{2}, 1τhf=1τh1τf\frac{1}{\tau_{hf}} = \frac{1}{\tau_{h}} - \frac{1}{\tau_{f}}.

多原子分子气体的 ff 也受到分子转动和振动的影响, 因此其平衡态 f(eq)f^{(eq)} 为:

f~(eq)(x,ξ,η,t)=ρ(2πRgT)D+K2exp[(ξu)2+η22RgT](3.8){\widetilde{f}}^{(eq)}\left( \mathbf{x,\xi,}\mathbf{\eta},t \right) = \frac{\rho}{\left( 2\pi R_g T \right)^{\frac{D + K}{2}}} \cdot \exp\left\lbrack - \frac{\left( \mathbf{\xi} - \mathbf{u} \right)^{2} + \mathbf{\eta}^{2}}{2R_g T} \right\rbrack \tag{3.8}

其中 η\mathbf{\eta} 是一个含有K个元素的内部自由度矢量.引入 f=fdη\overline{f} = \int f d\mathbf{\eta}, h=ξ2+η22fdη \overline{h}\mathbf{=}\int \frac{\mathbf{\xi}^{2}\mathbf{+}\mathbf{\eta}^{2}}{2}f d\mathbf{\eta\ } , 则可得到一组新的输运方程:

ft+ξf+aξf=1τf(ff(eq)),(3.9)\frac{\partial\overline{f}}{\partial t} + \mathbf{\xi} \cdot \nabla\overline{f} + \mathbf{a} \cdot \nabla_{\mathbf{\xi}}\overline{f} = - \frac{1}{\tau_{f}}\left( \overline{f} - {\overline{f}}^{(eq)} \right), \tag{3.9}

ht+ξh+aξhfξa=1τh(hh(eq))+Zτhf(ff(eq)).(3.10)\frac{\partial\overline{h}}{\partial t} + \mathbf{\xi} \cdot \nabla\overline{h} + a \cdot \nabla_{\mathbf{\xi}}\overline{h} - \overline{f}\mathbf{\xi \cdot a} = - \frac{1}{\tau_{h}}\left( \overline{h} - {\overline{h}}^{(eq)} \right) + \frac{Z}{\tau_{hf}}\left( \overline{f} - {\overline{f}}^{(eq)} \right). \tag{3.10}

它们的平衡态分别为

f(eq)=f~(eq)dη=ρ(2πRgT)D2exp[(ξu)22RgT],(3.11){\overline{f}}^{(eq)} = \int_{}^{}{\widetilde{f}}^{(eq)}d\mathbf{\eta =}\frac{\rho}{\left( 2\pi R_g T \right)^{\frac{D}{2}}}\exp\left\lbrack - \frac{\left( \mathbf{\xi} - \mathbf{u} \right)^{2}}{2R_g T} \right\rbrack, \tag{3.11}

h(eq)=ξ2+η22f~(eq)dη=ρ(ξ2+KRgT)2(2πRgT)D2exp[(ξu)22RgT].(3.12){\overline{h}}^{(eq)} = \int_{}^{}{\frac{\mathbf{\xi}^{2}\mathbf{+}\mathbf{\eta}^{2}}{2}{\widetilde{f}}^{(eq)}}d\mathbf{\eta =}\frac{\rho\left( \mathbf{\xi}^{2} + KR_g T \right)}{2 \cdot \left( 2\pi R_g T \right)^{\frac{D}{2}}}\exp\left\lbrack - \frac{\left( \mathbf{\xi} - \mathbf{u} \right)^{2}}{2R_g T} \right\rbrack. \tag{3.12}

可以看到多原子分子的流场与单原子分子的区别主要体现在能量上, 而 f(eq){\overline{f}}^{(eq)} 与式(1.10)一致.流场宏观量表示为:

ρ=fdξ,j=ρu=ξfdξ,ρE=hdξ.\rho = \int \overline{f}d\mathbf{\xi},\qquad \mathbf{j} =\rho\mathbf{u} = \int {\mathbf{\xi}\overline{f}}d\mathbf{\xi},\qquad \rho E = \int \overline{h}d\mathbf{\xi}.

定容比热比 cvc_{v} 和定压比热比 cpc_{p} 与内部自由度K相关:

cv=D+K2Rg,cp=D+K+22Rg.c_{v} = \frac{D + K}{2}R_g ,\qquad c_{p} = \frac{D + K + 2}{2}R_g .

对单原子分子而言, η=0 \mathbf{\eta = 0\ }K=0K = 0.

式(3.9)和式(3.10)在$\mathbf{c}_i $方向上的方程为

fit+cifi=1τf(fifi(eq))+Fi,(3.13)\frac{\partial{\overline{f}}_i}{\partial t} + \mathbf{c}_i \cdot \nabla{\overline{f}}_i = - \frac{1}{\tau_{f}}\left( {\overline{f}}_i - {\overline{f}}_i^{(eq)} \right) + F_i, \tag{3.13}

hit+cihi=1τh(hihi(eq))+Ziτhf(fifi(eq))+qi.(3.14)\frac{\partial{\overline{h}}_i}{\partial t} + \mathbf{c}_i \cdot \nabla{\overline{h}}_i = - \frac{1}{\tau_{h}}\left( {\overline{h}}_i - {\overline{h}}_i^{(eq)} \right) + \frac{Z_i}{\tau_{hf}}\left( {\overline{f}}_i - {\overline{f}}_i^{(eq)} \right) + q_i. \tag{3.14}

其中 Zi=ciuu22Z_i = \mathbf{c}_i \mathbf{\cdot u} - \frac{u^{2}}{2} .将其采用二阶积分格式离散后可得:

fi(x+ciΔt,t+Δt)fi(x,t)=ωf(fifi(eq))+(1ωf2)FiΔt,(3.15){\overline{f}}_i\left( \mathbf{x} + \mathbf{c}_i\Delta t,t + \Delta t \right) - {\overline{f}}_i\left( \mathbf{x},t \right) = - \omega_{f}\left( {\overline{f}}_i - {\overline{f}}_i^{(eq)} \right) + \left( 1 - \frac{\omega_{f}}{2} \right)F_i\Delta t, \tag{3.15}

hi(x+ciΔt,t+Δt)hi(x,t)=Ωh(hihi(eq))+Ziτhf(fifi(eq))+qi.(3.16){\overline{h}}_i\left( \mathbf{x} + \mathbf{c}_i\Delta t,t + \Delta t \right) - {\overline{h}}_i\left( \mathbf{x},t \right) = - \Omega_{h}\left( {\overline{h}}_i - {\overline{h}}_i^{(eq)} \right) + \frac{Z_i}{\tau_{hf}}\left( {\overline{f}}_i - {\overline{f}}_i^{(eq)} \right) + q_i .\tag{3.16}

Guo等 [17:3] 将平衡态使用Hermite展开对 f(eq){\overline{f}}^{(eq)}h(eq){\overline{h}}^{(eq)} 进行离散, 从而使低阶速度矩不会因高阶项的舍去而产生影响.由于涉及局部温度 TT , 因此展开过程是基于参考温度 T0T_0 进行的.离散平衡态为:

fi(eq)=ρwi[1+ciuRgT0+12(ciuRgT0)2u22RgT0],(3.17)\overline{f}_i^{(eq)} = \rho w_i \left[ 1 + \frac{\mathbf{c}_i \mathbf{\cdot u}}{R_g T_0} + \frac{1}{2}\left( \frac{\mathbf{c}_i \mathbf{\cdot u}}{R_g T_0} \right)^{2} - \frac{u^{2}}{2R_g T_0} \right], \tag{3.17}

hi(eq)=p0wi[ciuRgT0+(ciuRgT0)2u2RgT0+12(ci2RgT0D)]+fi(eq)E.(3.18)\overline{h}_i^{(eq)} = p_0 w_i \left[ \frac{\mathbf{c}_i \cdot \mathbf{u}}{R_g T_0} + \left( \frac{\mathbf{c}_i \cdot \mathbf{u}}{R_g T_0} \right)^{2} - \frac{u^{2}}{R_g T_0} + \frac{1}{2}\left( \frac{\mathbf{c}_i ^{2}}{R_g T_0} - D \right) \right] + \overline{f}_i^{(eq)}E. \tag{3.18}

源项为:

Fi=ρwi[ciaRgT0+(ciu)(cia)(RgT0)2auRgT0],qi=(ρwiERgT0+fi)(cia).(3.19)F_i = \rho w_i\left[ \frac{\mathbf{c}_{\mathbf{i}}\mathbf{\cdot}\mathbf{a}}{R_g T_{0}} + \frac{\left( \mathbf{c}_{\mathbf{i}}\mathbf{\cdot u} \right)\left( \mathbf{c}_{\mathbf{i}}\mathbf{\cdot}\mathbf{a} \right)}{\left( R_g T_{0} \right)^{2}} - \frac{\mathbf{a}\mathbf{\cdot u}}{R_g T_{0}} \right],\, q_i = \left( \frac{\rho w_iE}{R_g T_{0}} + f_i \right)\left( \mathbf{c}_{\mathbf{i}}\mathbf{\cdot}\mathbf{a} \right). \tag{3.19}

各宏观量则为:

ρ=ifi,j=ρu=icifi,ρE=ihi.\rho = \sum_i^{}{\overline{f}}_i,\, \mathbf{j} = \rho\mathbf{u} = \sum_i^{}{\mathbf{c}_i{\overline{f}}_i},\, \rho E = \sum_i^{}{\overline{h}}_i.

运动粘度为μ=τfp0\mu = \tau_{f}p_{0}, 热扩散系数为 κ=cvτhp0=D+K2Rgτhp0\kappa = c_{v}\tau_{h}p_{0} = \frac{D + K}{2}R_g \tau_{h}p_{0} .因此该模型的Prandtl数为 Pr=μcp/κ=γτf/τh\Pr = \mu c_{p}/\kappa = \gamma\tau_{f}/\tau_{h} , 其中 γ\gamma 为比热比.

[NOTE]
对应的参考压强 p0=ρRgT0p_{0} = \rho R_g T_{0}.

耦合DDF-LBE

何雅玲团队在Guo等 [17:4] 的理论基础上, 提出了一种耦合的DDF-LBE[30], 实现了热流和流场的耦合模拟.这里主要介绍耦合DDF-LBE与基于总能的DDF-LBE的区别.

为简化讨论, 考虑无外力项耦合DDF-LBE在 ci\mathbf{c}_i 方向上的演化方程组:

fit+cifi=1τf(fifi(eq)),(3.20)\frac{\partial f_i}{\partial t} + \mathbf{c}_{\mathbf{i}} \cdot \nabla f_i = - \frac{1}{\tau_{f}}\left( f_i - f_i^{(eq)} \right),\tag{3.20}

hit+cihi=1τh(hihi(eq))+(ciu)τhf(fifi(eq)).(3.21)\frac{\partial h_i}{\partial t} + \mathbf{c}_{\mathbf{i}} \cdot \nabla h_i = - \frac{1}{\tau_{h}}\left( h_i - h_i^{(eq)} \right) + \frac{\left( \mathbf{c}_{\mathbf{i}}\mathbf{\cdot u} \right)}{\tau_{hf}}\left( f_i - f_i^{(eq)} \right).\tag{3.21}

式(3.21)中将Zi=ciuu22Z_i = \mathbf{c}_{\mathbf{i}}\mathbf{\cdot u} - \frac{u^{2}}{2}换为(ciu)\left( \mathbf{c}_{\mathbf{i}}\mathbf{\cdot u} \right)并不会影响Navier-Stokes方程中的能量表示. 式(3.20)和式(3.21)中的两个平衡态分布(fα(eq)f_{\alpha}^{(eq)}hα(eq)h_{\alpha}^{(eq)})可以用两类方式进行求解.第一种方式是将其写作式(3.22)和式(3.23)的截断形式:

fα(eq)=Ai+Bσ(cαu)+Cαu2+Dα(cαu)2+Eαu2(cαu)+Fα(cαu)3,(3.22)\begin{aligned} f_{\alpha}^{(eq)} = & A_i + B_{\sigma} \cdot \left( \mathbf{c}_{\alpha} \cdot \mathbf{u} \right) + C_{\alpha}\left\| \mathbf{u} \right\|^{2} + D_{\alpha} \cdot \left( \mathbf{c}_{\alpha} \cdot \mathbf{u} \right)^{2} \\ & + E_{\alpha}\left\| \mathbf{u} \right\|^{2} \cdot \left( \mathbf{c}_{\alpha} \cdot \mathbf{u} \right) + F_{\alpha} \cdot \left( \mathbf{c}_{\alpha} \cdot \mathbf{u} \right)^{3} \end{aligned}, \tag{3.22}

hα(eq)=Kα+Lα(cαu)+Mαu2+Nα(cαu)2.(3.23)h_{\alpha}^{(eq)} = K_{\alpha} + L_{\alpha} \cdot \left( \mathbf{c}_{\alpha} \cdot \mathbf{u} \right) + M_{\alpha}\left\| \mathbf{u} \right\|^{2} + N_{\alpha} \cdot \left( \mathbf{c}_{\alpha} \cdot \mathbf{u} \right)^{2}. \tag{3.23}

则各系数可通过式(3.11)和式(3.12)的Chapman-Enskog展开得出[30:1].

第二种方式则采用其他 f(eq)f^{(eq)} 及其基于拉格朗日多项式的离散形式 fα(eq)f_{\alpha}^{(eq)} .令

h(eq)=[E+(ξu)u]f(eq)+ω(ξ,T)h^{(eq)} = \left[ E + \left( \mathbf{\xi} - \mathbf{u} \right) \cdot \mathbf{u} \right] f^{(eq)} + \omega(\mathbf{\xi},T)

其中 ω(ξ,T)\omega\left( \mathbf{\xi},T \right) 满足

ω(ξ,T)dξ=0,ω(ξ,T)ξdξ=0,ω(ξ,T)ξξTdξ=pRgTI\int {\omega(\mathbf{\xi},T)}d\mathbf{\xi =}0, \qquad \int {\omega(\mathbf{\xi},T)\mathbf{\xi}}d\mathbf{\xi =}0, \qquad \int {\omega\left( \mathbf{\xi},T \right)\mathbf{\xi}\mathbf{\xi}^{T}}d\mathbf{\xi =}pR_g T\mathbf{I}

因此可通过构建等价离散形式(式(3.24)):

hα(eq)=wαpc2RgT+[E+(cαu)u]fα(eq)(3.24)h_{\alpha}^{(eq)} = \frac{w_{\alpha}p}{c^{2}}R_g T + \left[ E + \left( \mathbf{c}_{\alpha} - \mathbf{u} \right) \cdot \mathbf{u} \right] \cdot f_{\alpha}^{(eq)} \tag{3.24}

计算总能的离散平衡态.

对于二维问题, 基于Qu等提出的圆函数[31]进行修改:

f(eq)={ρ2πcξu=cD(γ1)e and λ=ep,0otherwise(3.25)f^{(eq)} = \begin{cases} \frac{\rho}{2\pi c} & \parallel \mathbf{\xi} - \mathbf{u} \parallel = c \equiv \sqrt{D(\gamma - 1)e} \text{ and } \lambda = e_{p}, \\ 0 & \text{otherwise} \end{cases} \tag{3.25}

其中 λ\lambda 为静能, ep=[1D2(γ1)]ee_{p} = \left[ 1 - \frac{D}{2}(\gamma - 1) \right] e .对于三维问题, Li等[32]采用球面函数表示 f(eq)f^{(eq)}

f(eq)={ρ4πr2ξu=r=r,0otherwise.(3.26)f^{(eq)} = \begin{cases} \frac{\rho}{4\pi r^{2}} & \left\| \xi - \mathbf{u} \right\| = \left\| \mathbf{r} \right\| = r, \\ 0 & \text{otherwise}. \\ \end{cases} \tag{3.26}

以Li等[32:1]使用的D3Q27模型为例, 可得ρr2/3=p\rho r^{2}/3 = p, 因此对于理想气体有 r=3RgTr = \sqrt{3 R_g T}.

此外, Li等[30:2]指出:空间离散采用五阶WENO和TVD(total variation diminishing scheme)格式捕捉可压缩流中的不连续性;时间离散可采用二阶(或三阶)IMEX Runge-Kutta格式推进.虽然式(3.20)和式(3.21)同样可被离散为如下形式(ttt+Δtt + \Delta t和上标表示时间):

fαt+Δt=fαtΔt(cαfαt)+Δtτft(fα(eq),tfαt)(3.27)f_{\alpha}^{t + \Delta t} = f_{\alpha}^{t} - \Delta t \cdot \left( \mathbf{c}_{\alpha} \cdot \mathbf{\nabla}f_{\alpha}^{t} \right) + \frac{\Delta t}{\tau_{f}^{t}}\left( f_{\alpha}^{(eq),t} - f_{\alpha}^{t} \right) \tag{3.27}

hαt+Δt=hαt+Δt[1τht(hα(eq),thαt)(cαut)τhft(fα(eq),tfαt)(cαhαt)](3.28)h_{\alpha}^{t + \Delta t} = h_{\alpha}^{t} + \Delta t\left[ \frac{1}{\tau_{h}^{t}}\left( h_{\alpha}^{(eq),t} - h_{\alpha}^{t} \right) - \frac{\left( \mathbf{c}_{\alpha} \cdot \mathbf{u}^{t} \right)}{\tau_{hf}^{t}}\left( f_{\alpha}^{(eq),t} - f_{\alpha}^{t} \right) - \left( \mathbf{c}_{\alpha} \cdot \mathbf{\nabla}h_{\alpha}^{t} \right) \right] \tag{3.28}

但这种简单的形式仅适用于流场中无激波、无间断的情况, 且在时间上仅有一阶精度.


  1. He, X., & Luo, L.-S. (1997). Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. Physical Review E, 56(6), 6811–6817. DOI:10.1103/PhysRevE.56.6811 ↩︎ ↩︎ ↩︎

  2. Qian, Y. H., D’Humières, D., & Lallemand, P. (1992). Lattice BGK Models for Navier-Stokes Equation. Europhysics Letters, 17(6), 479–484. DOI:10.1209/0295-5075/17/6/001 ↩︎

  3. 杨鲲 & 单肖文. (2020). 多层速度格子Boltzmann方法进展及展望. 空气动力学学报, 40(3), 23–45. DOI:10.7638/kqdlxxb-2021.0348 ↩︎ ↩︎

  4. Koelman, J. M. V. A. (1991). A Simple Lattice Boltzmann Scheme for Navier-Stokes Fluid Flow. Europhysics Letters, 15(6), 603–607. DOI:10.1209/0295-5075/15/6/007 ↩︎

  5. He, X., & Luo, L.-S. (1997). Lattice Boltzmann Model for the Incompressible Navier–Stokes Equation. Journal of Statistical Physics, 88(3/4), 927–944. DOI:10.1023/B:JOSS.0000015179.12689.e4 ↩︎

  6. Zaki Abiza, Miguel Chavez, David M. Holman, & Ruddy Brionnaud. Prediction of Finned Projectile Aerodynamics using a Lattice-Boltzmann Method CFD Solution (Vol:10, No:05, 2016). World Academy of Science, Engineering and Technology. ↩︎

  7. 邵菲, 韩端锋, 刘强, & 谢伟. (2016). 熵格子Boltzmann方法的亚格子尺度模型. 中国舰船研究, 11(3), 43–47. DOI:10.3969/j.issn.1673-3185.2016.03.008 ↩︎

  8. Chen, Y., Ohashi, H., & Akiyama, M. (1994). Thermal lattice Bhatnagar-Gross-Krook model without nonlinear deviations in macrodynamic equations. Phys. Rev. E, 50, 2776–2783. DOI: 10.1103/PhysRevE.50.2776 ↩︎ ↩︎ ↩︎ ↩︎

  9. Chen, Y., Ohashi, H., & Akiyama, M. (1997). Two-Parameter Thermal Lattice BGK Model with a Controllable Prandtl Number. Journal of Scientific Computing, 12(2), 169–185. DOI:10.1023/A:1025621832215 ↩︎ ↩︎ ↩︎ ↩︎

  10. Nikolaos I. Prasianakis, & Konstantinos Boulouchos (2007). Lattice Boltzmann Method For Simulation Of Weakly Compressible Flows At Arbitrary Prandtl Number. International Journal of Modern Physics C, 18, 602-609. DOI:10.1142/S012918310701084X ↩︎ ↩︎ ↩︎ ↩︎ ↩︎

  11. Qian, Y.H (1993). Simulating thermohydrodynamics with lattice BGK models. Journal of Scientific Computing, 8, 231–242. DOI:10.1007/BF01060932 ↩︎ ↩︎

  12. Alexander, F., Chen, S., & Sterling, J. (1993). Lattice Boltzmann thermohydrodynamics. Physical Review E, 47, R2249–R2252. DOI:10.1103/PhysRevE.47.R2249 ↩︎ ↩︎

  13. Frapolli, N., Chikatamarla, S., & Karlin, I. (2014). Multispeed entropic lattice Boltzmann model for thermal flows. Physical Review E, 90, 043306. DOI:10.1103/PhysRevE.90.043306 ↩︎ ↩︎ ↩︎ ↩︎

  14. N.I. Prasianakis, S.S. Chikatamarla, I.V. Karlin, S. Ansumali, & K. Boulouchos (2006). Entropic lattice Boltzmann method for simulation of thermal flows. Mathematics and Computers in Simulation, 72(2), 179-183. DOI:10.1016/j.matcom.2006.05.012 ↩︎ ↩︎ ↩︎

  15. A. Bartoloni, Claudia Battista, Simone Cabasino, P. S. Paolucci, J. Pech, Renata Sarno, Gian Marco Todesco, Mario Torelli, Walter Tross, Piero Vicini, Roberto Benzi, Nicola Cabibbo, Federico Massaioli, & Raffaele Tripiccione (1993). LBE simulations of Rayleigh-Bénard convection on the APE100 parallel processor. International Journal of Modern Physics C, 04, 993-1006. DOI:10.1142/S012918319300077X ↩︎ ↩︎

  16. Xiaoyi He, Shiyi Chen, & Gary D. Doolen (1998). A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit. Journal of Computational Physics, 146(1), 282-300. DOI:10.1006/jcph.1998.6057 ↩︎ ↩︎

  17. Guo, Z., Zheng, C., Shi, B., & Zhao, T. (2007). Thermal lattice Boltzmann equation for low Mach number flows: Decoupling model. Physical Review E, 75, 036704. DOI:10.1103/PhysRevE.75.036704 ↩︎ ↩︎ ↩︎ ↩︎ ↩︎

  18. Dipankar Chatterjee (2009). An enthalpy-based thermal lattice Boltzmann model for non-isothermal systems. Europhysics Letters, 86(1), 14004. DOI:10.1209/0295-5075/86/14004 ↩︎ ↩︎

  19. Chen, Z., Shu, C., & Tan, D. (2017). Three-dimensional simplified and unconditionally stable lattice Boltzmann method for incompressible isothermal and thermal flows. Physics of Fluids, 29(5), 053601. DOI:10.1063/1.4983339 ↩︎

  20. Y. Wang, C. Shu, & C.J. Teo (2014). Thermal lattice Boltzmann flux solver and its application for simulation of incompressible thermal flows. Computers & Fluids, 94, 98-111. DOI: 10.1016/j.compfluid.2014.02.006 ↩︎

  21. Zheng, L., Shi, B., & Guo, Z. (2008). Multiple-relaxation-time model for the correct thermohydrodynamic equations. Physical Review E, 78, 026705. DOI:10.1103/PhysRevE.78.026705 ↩︎

  22. Ai-Guo Xu, Guang-Cai Zhang, Yan-Biao Gan, Feng Chen & Xi-Jun Yu (2012). Lattice Boltzmann modeling and simulation of compressible flows. Frontiers of Physics. 7, 582–600. DOI:10.1007/s11467-012-0269-5 ↩︎

  23. Feng Chen, Aiguo Xu, Guangcai Zhang, Yingjun Li, & Sauro Succi (2010). Multiple-relaxation-time lattice Boltzmann approach to compressible flows with flexible specific-heat ratio and Prandtl number. Europhysics Letters, 90(5), 54003. DOI:10.1209/0295-5075/90/54003 ↩︎

  24. Sharma, K. V., Straka, R., & Tavares, F. W. (2020). Current status of Lattice Boltzmann Methods applied to aerodynamic, aeroacoustic, and thermal flows. Progress in Aerospace Sciences, 115, 100616. DOI:10.1016/j.paerosci.2020.100616 ↩︎

  25. Kataoka, T., & Tsutahara, M. (2004). Lattice Boltzmann model for the compressible Navier-Stokes equations with flexible specific-heat ratio. Physical Review E, 69, 035701. DOI:10.1103/PhysRevE.69.035701 ↩︎

  26. Gan, Y.-B., Xu, A.-G., Zhang, G.-C., & Li, Y.-J. (2011). Flux Limiter Lattice Boltzmann Scheme Approach to Compressible Flows with Flexible Specific-Heat Ratio and Prandtl Number. Communications in Theoretical Physics, 56(3), 490–498. DOI:10.1088/0253-6102/56/3/18 ↩︎

  27. Hosseini, S. A., Atif, M., Ansumali, S., & Karlin, I. V. (2023). Entropic lattice Boltzmann methods: A review. Computers & Fluids, 259, 105884. DOI:10.1016/j.compfluid.2023.105884 ↩︎

  28. Anirudh Jonnalagadda; Atul Sharma; Amit Agrawal (2021). Single relaxation time entropic lattice Boltzmann methods: A developer’s perspective for stable and accurate simulations. Computers & Fluids, 215, 104792. DOI:10.1016/j.compfluid.2020.104792 ↩︎

  29. Yong, Wen-an., Luo, Li-Shi (2005). Nonexistence of H Theorem for Some Lattice Boltzmann Models. Journal of Statistical Physics, 121, 91–103. DOI:10.1007/s10955-005-5958-9 ↩︎

  30. Li, Q., He, Y. L., Wang, Y., & Tao, W. Q. (2007). Coupled double-distribution-function lattice Boltzmann method for the compressible Navier-Stokes equations. Physical Review E, 76(5), 056705. DOI:10.1103/PhysRevE.76.056705 ↩︎ ↩︎ ↩︎

  31. Qu, K., Shu, C., & Chew, Y. (2007). Alternative method to construct equilibrium distribution functions in lattice-Boltzmann method simulation of inviscid compressible flows at high Mach number. Physical Review E, 75, 036706. ↩︎

  32. Q. Li, Y.L. He, Y. Wang, & G.H. Tang (2009). Three-dimensional non-free-parameter lattice-Boltzmann model and its application to inviscid compressible flows. Physics Letters A, 373(25), 2101-2108. DOI:10.1016/j.physleta.2009.04.036 ↩︎ ↩︎

工位上换了新电脑,于是需要重新安装 RDPWrap 让机子上的 Windows 11 家庭版能够打开 RDP 远程连接。

  1. 安装 RDPWrap

RDPWrap 的 Github 网址上下载它的最新 Release 【其实已经很久没更新了】。 解压后运行 install.bat 进行安装。

  1. 出现 Listener state: Not-listening 的解决方案。

1)首先在 sebaxakerhtc/rdpwrap.ini 这个仓库里下一份最新的 rdpwrap.ini,一般来说大多数Windows OS 版本号它都有配置。
2)然后将自行下载的 rdpwrap.ini 替换掉原始文件[1],原始文件位于 C:\Program Files\RDP Wrapper 文件夹。
3)替换完成后,回到解压 RDPWrap 的文件夹并在此处打开终端,运行 .\RDPWInst.exe -r [2]。这之后我的电脑就设置完成了。


  1. 执行此操作前先在具有管理员权限的终端里把对应远程服务关了:net stop TermService↩︎

  2. 这里参考了 Github 上的 issue 999↩︎

工位上换了新电脑,于是需要给新电脑安装 Linux 子系统(WSL 2)用来跑代码。
这里简单记录一下我的安装过程。

  1. 安装 WSL 2

这里主要是按照微软官方的教程[1][2]来做。

先在 powershell 里启用 WSL 和虚拟机功能

1
2
dism.exe /online /enable-feature /featurename:Microsoft-Windows-Subsystem-Linux /all /norestart  
dism.exe /online /enable-feature /featurename:VirtualMachinePlatform /all /norestart

然后安装适用于 x64 计算机的 WSL2 Linux 内核更新包,并在 powershell 里将 WSL 2 设置为默认。

1
wsl --set-default-version 2

WSL 2 的自定义安装位置参考的是CSDN上给的做法:我下载了Ubuntu 24.04 的 AppxBundle,从包里提取 Ubuntu_2404.0.5.0_x64.Appx ,然后在我想要的位置将它当成 zip 解压。最后运行一下 ubuntu2404.exe 即可。

  1. 导入Windows系统的字体

这里主要参照知乎上给的方法[3],做了个软链接。

1
2
3
sudo mkdir /usr/share/fonts/win11 
sudo ln -s /mnt/c/Windows/Fonts/* /usr/share/fonts/win11
fc-cache -fv
  1. 安装CUDA

Nvidia 的安装教程 和它给的下载链接主要是针对 WSL-Ubuntu 的。

我用的网络下载,在shell里就应该是:

1
2
3
4
5
wget https://developer.download.nvidia.com/compute/cuda/repos/wsl-ubuntu/x86_64/cuda-keyring_1.1-1_all.deb
sudo dpkg -i cuda-keyring_1.1-1_all.deb
sudo apt-get update
sudo apt-get -y install cuda-toolkit-12-6
sudo apt-get -y install cuda

即使是这么安装完,还是会有非常幽默的 nvcc 找不到、 libcuda.so 找不到等问题。我就参照网上的做法[4][5]~/.bashrc 里补充了这么一段:

1
2
3
4
5
6
# >>> Nvidia CUDA initialize >>>
export CUDA_HOME=/usr/local/cuda
export PATH=${CUDA_HOME}/bin:${PATH}
export LD_LIBRARY_PATH=${CUDA_HOME}/lib64:${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}
export LD_LIBRARY_PATH=/usr/lib/wsl/lib/:${LD_LIBRARY_PATH}
# <<< Nvidia CUDA initialize <<<

搞定。


  1. https://learn.microsoft.com/zh-cn/windows/wsl/install ↩︎

  2. https://learn.microsoft.com/zh-cn/windows/wsl/install-manual ↩︎

  3. https://zhuanlan.zhihu.com/p/683470601 ↩︎

  4. https://stackoverflow.com/questions/68221962/nvcc-not-found-but-cuda-runs-fine ↩︎

  5. https://github.com/microsoft/WSL/issues/8587 ↩︎

题 1

一无限流体被单位质量力 μr3/2\mu r^{-3/2} 的力作用,方向指向原点。若开始时流体静止且其中有一个空腔 r=cr=c,证明:在过 (25μ)1/2c5/4(\frac{2}{5 \mu})^{1/2} c^{5/4} 时间之后,空腔被流体填满。

这里假设流体不可压缩。取球坐标系 (r,θ,ϕ)(r, \theta, \phi) 进行计算,并以空腔中心为原点。

思路1:功能关系

假设当 t=t0t=t_0 时,空腔边界的球壳形薄液面位于 r=r0r=r_0 处。该薄液面厚度为 dr0\mathrm{d}r_0, 液面速度 vr(r0,t0)=0v_r (r_0, t_0) = 0

t0t_0 时刻,位于 r=r0r=r_0 处且厚度为 dr0\mathrm{d}r_0 的液面向内推动了 dr\mathrm{d} r, 由于该液面的质量积为 dm=ρ(4πr02)dr0\mathrm{d}m = \rho (4 \pi r_0^2) \mathrm{d}r_0,则这一过程中体积力(Fbr=μr3/2F_{br} = -\mu r^{-3/2})所作的功为

dW=dmr0+Fbrdr=dm(2μr01/2)=8μρπr03/2dr0.\mathrm{d} W = \mathrm{d}m \cdot \int_{r_0}^{+\infty} F_{br} \mathrm{d}r = \mathrm{d}m \cdot (-2 \mu r_0^{-1/2}) = -8 \mu \rho \pi r_0^{3/2} \mathrm{d}r_0 .

t1=t0+δtt_1 = t_0 + \delta_t 时,该液面前进到 r=r1r=r_1 位置(r1r0r_1 \le r_0),此时液面速度为 vr(r1,t1)=dr1dtv_r(r_1, t_1) = \dfrac{\mathrm{d} r_1}{\mathrm{d} t}

所以薄液面从 r0r_0 运动到 r1r_1 过程中, FbrF_{br} 做的总功为:

W=r0r1dW=165μρπ(r05/2r15/2)W = \int_{r_0}^{r_1} \mathrm{d} W = \frac{16}{5} \mu \rho \pi (r_0^{5/2} - r_1^{5/2})

记速度 v=[vr,0,0]T\bold{v} = [v_r, 0, 0]^T。代入球坐标系下的连续性方程,化简得:

(r2vr)r=0\frac{\partial (r^2 v_r)}{\partial r} = 0

rr 积分,可得 r2vr=f(t)r^2 v_r = f(t) ,其中 f(t)f(t) 为待定函数。根据这一方程,可得在 t1t_1 时刻满足 r02vr(r0,t1)=r12vr(r1,t1)r_0^2 v_r(r_0, t_1) = r_1^2 v_r(r_1, t_1)

t1t_1 时刻, r=r0r=r_0 层液体的动能表示为

dE=12dmvr2(r0,t1)=2ρπr14r02vr2(r1,t1)dr0\mathrm{d} E = \frac{1}{2} \cdot \mathrm{d}m \cdot v_r^2(r_0, t_1) = 2 \rho \pi \cdot \frac{r_1^4}{r_0^2} \cdot v_r^2(r_1, t_1) \mathrm{d}r_0

这一时刻下整个系统的动能为:

E=r1+dE=2ρπr13vr2(r1,t1)E = \int_{r_1}^{+\infty} \mathrm{d} E = 2 \rho \pi r_1^3 \cdot v_r^2(r_1, t_1)

据此,我们令 t0t_0 时刻即为初始时刻(t0=0t_0=0),满足 r0=cr_0 = cvr(t=0)=0v_r(t=0)=0。所以 t0t_0 时刻的动能为 0。对 t0t_0t1t_1 时刻的过程, FbrF_{br} 所做的功 WW 全部转化为流场动能 EE,即:W=EW=E

165μρπ(c5/2r15/2)=2ρπr13vr2(r1,t1)0\frac{16}{5} \mu \rho \pi (c^{5/2} - r_1^{5/2}) = 2 \rho \pi r_1^3 \cdot v_r^2(r_1, t_1) - 0

也就有:

vr(r1,t1)=dr1dt=8μ5r13/2(c5/2r15/2)v_r(r_1, t_1) = \frac{\mathrm{d} r_1}{\mathrm{d} t} = \sqrt{\frac{8 \mu}{5}} \cdot \sqrt{r_1^{-3/2} (c^{5/2} - r_1^{5/2})}

分离变量后解得:

1625(c5/2r15/2)=(8μ5t+c1)2\frac{16}{25} (c^{5/2} - r_1^{5/2}) = \left( \sqrt{\frac{8 \mu}{5}} t + c_1 \right)^2

其中 c1c_1 为待定常数。当 t=0t=0 时, r1=cr_1 = c,所以 c1=0c_1 = 0

假设当 t=Tmt=T_{m}r1=0r_1=0, 则可解得: Tm=(25μ)1/2c5/4T_m = (\frac{2}{5 \mu})^{1/2} c^{5/4}

思路2:联立动量方程和连续性方程

记速度 v=[vr,0,0]T\bold{v} = [v_r, 0, 0]^T。代入连续性方程,化简得:

(r2vr)r=0\frac{\partial (r^2 v_r)}{\partial r} = 0

rr 积分,可得 r2vr=f(t)r^2 v_r = f(t) ,其中 f(t)f(t) 为待定函数。

该场景下的运动方程为

ρDvDt=ρFb\rho \frac{\mathrm{D} \bold{v}}{\mathrm{D} t} = \rho \bold{F}_b

所以,动量方程写作:

vrt+vrvrr=Fbr\frac{\partial v_r}{\partial t} + v_r \frac{\partial v_r}{\partial r} = F_{br}

由题目可知 Fbr=μr3/2F_{br} = - \mu r^{-3/2}。将 r2vr=f(t)r^2 v_r = f(t) 代入到运动方程中,得:

f(t)r2+vrvrr=μr3/2\frac{f'(t)}{r^2} + v_r \frac{\partial v_r}{\partial r} = -\mu r^{-3/2}

R(t)R(t) 为真空球面的半径变化,则 vr(r=R(t))=R(t)v_r (r=R(t)) = R'(t)。对 rr 积分,得到:

R(t)f(t)r2+vrvrrdr=μR(t)r3/2dr[f(t)r+(vr)22]R(t)=2μrR(t)\begin{aligned} \int_{\infty}^{R(t)} \frac{f'(t)}{r^2} + v_r \frac{\partial v_r}{\partial r} \mathrm{d}r &= -\mu \int_{\infty}^{R(t)} r^{-3/2} \mathrm{d}r \\ \left.\left[ -\frac{f'(t)}{r} + \frac{(v_r)^2}{2} \right]\right|_{\infty}^{R(t)} &= \left. \frac{2 \mu}{\sqrt{r}} \right|_{\infty}^{R(t)} \end{aligned}

代入 vr(r=)=0v_r (r=\infty)=0,整理得到:

f(t)R(t)+(R(t))22=2μR(t)-\frac{f'(t)}{R(t)} + \frac{(R'(t))^2}{2} = \frac{2 \mu}{\sqrt{R(t)}}

在边界上有 f(t)=R2(t)vr(r=R(t))=R2Rf(t) = R^2(t) \cdot v_r (r=R(t)) = R^2 \cdot R' , 所以 f(t)=2RR2+R2Rf'(t) = 2 R R'^{2} + R^2 \cdot R{''} 。上式化简为:

32R2RR=2μR-\frac{3}{2} R'^{2} -R \cdot R{''} = \frac{2 \mu}{\sqrt{R}}

由于 d(R1/2(R)2)dR=R1/22(R)2+2R1/2R\frac{\mathrm{d} (R^{1/2} (R')^{2})}{\mathrm{d} R} =\frac{R^{-1/2}}{2} (R')^{2} + 2 R^{1/2} R{''},上式整理得:

d(R1/2(R)2)8μ+5(R1/2(R)2)=dR2R-\frac{\mathrm{d}(R^{1/2} (R')^{2})}{8\mu + 5 (R^{1/2} (R')^{2})} = \frac{\mathrm{d}R}{2R}

分离变量后并整理出 RR' ,可表示为:

R=R1/25(c1R5/28μ)R' = \sqrt{\frac{R^{-1/2}}{5} (c_1 R^{-5/2} - 8 \mu) }

R(t=0)=cR(t=0)=cR(t=0)=0R'(t=0) = 0 代入,解得: c1=8μc5/2c_1 = 8 \mu c^{5/2} 。即:

R=R1/258μ(c5/2R5/21)=8μ5R1/2(c5/2R5/21)R' = \sqrt{\frac{R^{-1/2}}{5} 8 \mu (c^{5/2} R^{-5/2} - 1)} \\= \sqrt{\frac{8 \mu}{5}} \cdot \sqrt{R^{-1/2} (c^{5/2} R^{-5/2} - 1)}

综上,真空球面的总运动时长表示为:

T=0c1RdR=58μ0c[R1/2((cR)5/21)]1/2dRT = \int_{0}^{c} \frac{1}{R'} \mathrm{d}R = \sqrt{\frac{5}{8 \mu}} \int_{0}^{c} \left[ R^{-1/2} \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}R

因为:

0c[R1/2((cR)5/21)]1/2dR=0cR1/4[((cR)5/21)]1/2dR=450c[((cR)5/21)]1/2d(R5/4)=t=R5/4450c5/4[(c5/2t21)]1/2dt=450c5/4t[(c5/2t2)]1/2dt=45(12)0c5/4[(c5/2t2)]1/2d(c5/2t2)=45(12)2(c5/2t2)1/20c5/4=45c5/4\begin{aligned} & \int_{0}^{c} \left[ R^{-1/2} \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}R \\ =& \int_{0}^{c} R^{1/4} \left[ \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}R \\ =& \frac{4}{5} \int_{0}^{c} \left[ \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}(R^{5/4}) \\ \xlongequal{t=R^{5/4}} & \frac{4}{5} \int_{0}^{c^{5/4}} \left[ \left( c^{5/2} t^{-2} - 1 \right) \right]^{-1/2} \mathrm{d}t \quad= \frac{4}{5} \int_{0}^{c^{5/4}} t \left[ \left( c^{5/2} - t^2 \right) \right]^{-1/2} \mathrm{d}t \\ =& \frac{4}{5} \cdot \left(-\frac{1}{2}\right) \int_{0}^{c^{5/4}} \left[ \left( c^{5/2} - t^2 \right) \right]^{-1/2} \mathrm{d}(c^{5/2} - t^2) \\ =& \frac{4}{5} \cdot \left(-\frac{1}{2}\right) \left.2 (c^{5/2} - t^2)^{1/2}\right|_{0}^{c^{5/4}} = \frac{4}{5} c^{5/4} \end{aligned}

所以 T=58μ45c5/4=25μc5/4T = \sqrt{\frac{5}{8 \mu}} \cdot \frac{4}{5} c^{5/4} = \sqrt{\frac{2}{5 \mu}} c^{5/4} .

证毕。

题 2

体积为 43πa3\frac{4}{3} \pi a^3 的液体充满于两个同心球面之间,外球面有压强 π\pi 作用,无质量力,内球面压强为 0 。开始时液体静止,内球面半径为 2a2a 。证明:当内球面半径变为 aa 时,其速度为 (14π3ρ21/321/31)12\displaystyle \left(\frac{14 \pi}{3 \rho} \frac{2^{1/3}}{2^{1/3} - 1}\right)^{\frac{1}{2}}

这里假设流体不可压缩。取球坐标系 (r,θ,ϕ)(r, \theta, \phi) 进行计算,并以空腔中心为原点。

记初始状态下(即 t=0t=0 时刻),内球面半径为 R0=2aR_0=2 a,外球面半径为 R1R_1;当内球面半径变为 aa 时(即 t1t_1 时刻),此时的内、外球面半径为分别为 R01=aR_{01} = aR11R_{11}。由于流体体积 V=43πa3V = \frac{4}{3} \pi a^3,可得:

V=43π(R13R03)=43π(R113R013)\begin{aligned} V &= \frac{4}{3} \pi (R_1^3 - R_0^3) \\ &= \frac{4}{3} \pi (R_{11}^3 - R_{01}^3) \end{aligned}

解得: R1=91/3aR_1 = 9^{1/3} aR11=21/3aR_{11} = 2^{1/3} a

记速度 u=[vr,0,0]T\bold{u} = [v_r, 0, 0]^T。代入连续性方程,化简得:

(r2vr)r=0\frac{\partial (r^2 v_r)}{\partial r} = 0

rr 积分,可得 r2vr=f(t)r^2 v_r = f(t) ,其中 f(t)f(t) 为待定函数。

t[0,t1]t \in [0,t_1] 时间段,外界压强在外表面做的功为:

W=R1R11(π)S(r)dr=283π2a3W = \int_{R_{1}}^{R_{11}} (-\pi) S(r) \mathrm{d} r = \frac{28}{3} \pi^2 a^3

其中 S(r)=4πr2S(r) = 4 \pi r^2 为球的表面积公式。

t1t_1 时刻, r2vr=f(t1)r^2 v_r = f(t_1) 为定值。若记此时刻的内球面速度为 u0=[v0,0,0]T\bold{u}_0 = [v_0, 0, 0]^Trr 处的流体速度为 ur=[vr,0,0]T\bold{u}_r = [v_r, 0, 0]^Tr[R01,R11]r \in [R_{01}, R_{11}]),则 R012v0=r2vrR_{01}^2 v_0 = r^2 v_r 。 所以 vr=(R01r)2v0v_r = \left( \frac{R_{01}}{r} \right)^2 v_0

t1t_1 时刻的流场总动能表示为:

Et1=R01R1112vr2(ρS(r)dr)=2πρv02a321/3121/3E_{t_1} = \int_{R_{01}}^{R_{11}} \frac{1}{2} v_r^2 \cdot (\rho S(r) \mathrm{d}r) = 2 \pi \rho v_0^2 a^3 \frac{2^{1/3} - 1}{2^{1/3}}

由于初始时刻的动能 E0=0E_0 = 0,根据功能关系可知:

W=Et1E0W = E_{t_1} - E_0

解得 t1t_1 时刻的内球面速度为: v0=(14π3ρ21/321/31)12\displaystyle v_0= - \left(\frac{14 \pi}{3 \rho} \frac{2^{1/3}}{2^{1/3} - 1}\right)^{\frac{1}{2}} ,方向指向原点。

本文同步发表于知乎CSDN

我在阅读其他文献时,看到有文献提到 Kupershtokh 的一篇关于 LBM 外力项的文章[1],他提出的格式在其他文章中又被称为 Exact difference method(EDM)

本想搜索这篇文章的PDF进行阅读,但最后只找到一个网站介绍了部分信息。这份粗略笔记的内容,都是基于其他文献[2][3][4]中对该方法的介绍而整理的。

介绍

Boltzmann-BGK方程的近似

密度分布函数 ff 的 Boltzmann-BGK方程写作:

ft+ξf+aξf=1τf[ff(eq)](1)\frac{\partial f}{\partial t} + \bold{\xi} \cdot \nabla f + \boldsymbol{a} \cdot \nabla_{\bold{\xi}} f = -\frac{1}{\tau_f} \left[ f - f^{(eq)}\right] \tag{1}

其中:

f(eq)=ρ(2πRgT)D/2exp((ξu)22RgT)(2)f^{(eq)} = \frac{\rho}{(2 \pi R_g T)^{D/2}} \exp{\left(- \frac{(\bold{\xi} - \bold{u})^2}{2 R_g T}\right)} \tag{2}

符号含义分别为

  • ξ\bold{\xi} : 粒子速度;
  • ρ\rho : 流体密度;
  • u\bold{u} : 流场宏观速度;
  • RgR_g : 气体常数;
  • TT : 温度;
  • a=du/dt\boldsymbol{a}=\mathrm{d}\bold{u}/\mathrm{d}t : 流场加速度。

Knudsen数较小时,Boltzmann方程的体力项可以基于f(eq)f^{(eq)}进行近似,即

ξfξf(eq)\nabla_{\bold{\xi}} f \approx \nabla_{\bold{\xi}} f^{(eq)}

根据式(2), f(eq)f^{(eq)} 是分子随机速度 C=ξu\bold{C} = \bold{\xi} - \bold{u} 的指数函数,则:

f(eq)u=f(eq)ξ=ξf(eq)\frac{\partial f^{(eq)}}{\partial \bold{u}} = -\frac{\partial f^{(eq)}}{\partial \bold{\xi}} = -\nabla_{\bold{\xi}} f^{(eq)}

假设密度 ρ\rho 为常数,则 dρdt=0\frac{\mathrm{d} \rho}{\mathrm{d} t}=0,那么 f(eq)f^{(eq)} 的全导数为:

df(eq)dt=f(eq)ududt+f(eq)ρdρdtf(eq)ududt=aξf(eq)\begin{aligned} \frac{\mathrm{d} f^{(eq)}}{\mathrm{d} t} &= \frac{\partial f^{(eq)}}{\partial \bold{u}} \cdot \frac{\mathrm{d} \bold{u}}{\mathrm{d} t} + \frac{\partial f^{(eq)}}{\partial \rho} \cdot \frac{\mathrm{d} \rho}{\mathrm{d} t} \\ &\approx \frac{\partial f^{(eq)}}{\partial \bold{u}} \cdot \frac{\mathrm{d} \bold{u}}{\mathrm{d} t} = -\boldsymbol{a} \cdot \nabla_{\bold{\xi}} f^{(eq)} \end{aligned}

也就是说,Kupershtokh 等[1:1][2:1]近似的 Boltzmann-BGK方程写作:

ft+ξf=1τf[ff(eq)]+df(eq)dt(3)\frac{\partial f}{\partial t} + \bold{\xi} \cdot \nabla f = -\frac{1}{\tau_f} \left[ f - f^{(eq)}\right] + \frac{\mathrm{d} f^{(eq)}}{\mathrm{d} t} \tag{3}

近似方程的离散化

将式(3)在[t,t+δt][t,t+\delta_t]内沿着eα\bold{e}_\alpha方向积分,可得

fα(x+eαδt,t+δt)fα(x,t)=1τ[fα(x,t)fα(eq)(ρ,u)]+[fα(eq)(u(t+δt))fα(eq)(u(t))](4)\begin{aligned} f_\alpha (\bold{x} + \bold{e}_\alpha \delta_t, t+\delta_t) - f_\alpha (\bold{x}, t) =& -\frac{1}{\tau} \left[ f_\alpha (\bold{x}, t) - f_\alpha^{(eq)} (\rho, \bold{u^*}) \right] \\ & + \left[ f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t)) \right] \end{aligned} \tag{4}

其中:

  • δt\delta_t 为时间步长;
  • τ=τf/δt\tau = \tau_f / \delta_t 为无量纲松弛时间;
  • x\bold{x} 为网格点位置,ρ,u\rho, \bold{u^*} 分别为该点在 tt 时刻的密度和速度(与宏观速度 u\bold{u} 有区别);
  • fα,fα(eq)f_\alpha, f_\alpha^{(eq)} 分别为 eα\bold{e}_\alpha 方向的分布函数平衡态

fα(eq)=wαρ[1+eαucs2+(eαu)22cs4u22cs2]=wαρ[1+eαucs2+12cs4(uu):(eαeαcs2[I])](4-1)\begin{aligned} f_\alpha^{(eq)} &= w_\alpha \rho \cdot \left[ 1 + \frac{\bold{e}_\alpha \cdot \bold{u}}{c_s^2} + \frac{(\bold{e}_\alpha \cdot \bold{u})^2}{2 c_s^4} - \frac{|\bold{u}|^2}{2 c_s^2} \right] \\ &= w_\alpha \rho \cdot \left[ 1 + \frac{\bold{e}_\alpha \cdot \bold{u}}{c_s^2} + \frac{1}{2 c_s^4}(\bold{u}\bold{u}):(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \\ \end{aligned} \tag{4-1}

这里 wαw_\alphaeα\bold{e}_\alpha 的权重系数, [I][\bold{I}] 为单位矩阵;
两个同样为m*n大小的矩阵 [A][\bold{A}][B][\bold{B}] ,则 [A]:[B]=i=1mj=1nAijBij\displaystyle [\bold{A}] : [\bold{B}] = \sum^{m}_{i=1}\sum^{n}_{j=1} A_{ij} B_{ij}

  • FαF_\alpha 为外力项,其表达式基于外力 F=ρdudt=ρa\bold{F} = \rho \dfrac{\mathrm{d} \bold{u^*}}{\mathrm{d} t} = \rho \bold{a} 进行构建。

Fα=fα(eq)(u(t+δt))fα(eq)(u(t))F_\alpha = f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t))

并且

u(t)=(αfαeα)/(αfα)u(t+δt)=u(t)+tt+δtdudtdt=u(t)+tt+δtFρdtu(t)+Fρδt(5)\begin{aligned} \bold{u^*}(t) &= \left(\sum_{\alpha} f_\alpha \bold{e}_\alpha \right) / \left(\sum_{\alpha} f_\alpha \right) \\ \bold{u^*}(t+\delta_t) &= \bold{u^*}(t) + \int_{t}^{t+\delta_t} \frac{\mathrm{d} \bold{u^*}}{\mathrm{d} t} \mathrm{d} t \\ &= \bold{u^*}(t) + \int_{t}^{t+\delta_t} \frac{\bold{F}}{\rho} \mathrm{d} t \approx \bold{u^*}(t) + \frac{\bold{F}}{\rho} \delta_t \end{aligned} \tag{5}

这意味着加速度 a=F/ρ\bold{a}=\bold{F}/\rho[t,t+δt][t,t+\delta_t] 时间段内被视为常数。

Kupershtokh等在这篇文章[2:2]说明,在使用该源项时,流体宏观量的计算应表示为:

ρ=αfα,ρu=ieαfα+δtF2(6)\rho = \sum_{\alpha} f_\alpha ,\quad \rho \bold{u} = \sum_{i} \bold{e}_\alpha f_\alpha + \frac{\delta_t \bold{F}}{2} \tag{6}

Li 等[4:1]通过 Chapman-Enskog 展开证明:该格式向宏观方程引入的误差项为:

ErrorEDM=δt2(FF4ρ)(7)\mathbf{Error}_{\mathrm{EDM}} = -\delta_t^2 \nabla\cdot\left( \frac{\bold{FF}}{4 \rho} \right) \tag{7}

EDM 源项的展开式

为便于分析,这里采用如下书写形式的平衡态:

fα(eq)=wαρ[1+eαucs2+12cs4(uu):(eαeαcs2[I])]f_\alpha^{(eq)} = w_\alpha \rho \cdot \left[ 1 + \frac{\bold{e}_\alpha \cdot \bold{u}}{c_s^2} + \frac{1}{2 c_s^4}(\bold{u}\bold{u}):(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right]

将其带入 EDM 源项,记 Δu=Fρδt\Delta\bold{u} = \frac{\bold{F}}{\rho} \delta_t,整理得:

Fα=fα(eq)(u(t+δt))fα(eq)(u(t))=ρwα{eαΔucs2+12cs4[(u+Δu)(u+Δu)uu]:(eαeαcs2[I])}\begin{aligned} F_{\alpha} &= f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t)) \\&= \rho w_\alpha \left\{ \frac{\bold{e}_\alpha \cdot \Delta\bold{u}}{c_s^2} + \newline \frac{1}{2 c_s^4} \left[ (\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) \right.\right. \\ &\quad \left. - \bold{u}^* \bold{u}^* \right] : \left.(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right\} \end{aligned}

这里涉及矩阵 (u+Δu)(u+Δu)(\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u})uu\bold{u}^* \bold{u}^* 的相减,由于:

(u+Δu)(u+Δu)=(ui+Δui)(uj+Δuj),uu=uiuj(\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) = (u^*_{i} + \Delta u_i)(u^*_{j} + \Delta u_j), \quad \bold{u}^* \bold{u}^* = u^*_{i} u^*_{j}

其中 uiu_i^*Δui=Fiδt/ρ\Delta u_i = F_i \delta_{t} / \rho 分别为 u\bold{u}^*Δu\Delta\bold{u} 在 i 轴上的分量,所以:

(u+Δu)(u+Δu)uu=(ui+Δui)(uj+Δuj)uiuj=δtρ(uiFj+ujFi)+δt2ρ2FiFj=δtρ(ui+Fiδt2ρ)Fj+δtρ(uj+Fjδt2ρ)Fi\begin{aligned} (\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) - \bold{u}^* \bold{u}^{*} &= (u^*_{i} + \Delta u_i)(u^*_{j} + \Delta u_j) - u^*_{i} u^*_{j} \\ &= \frac{\delta_t}{\rho} (u_i^* F_j + u_j^* F_i) + \frac{\delta_t^2}{\rho^2} F_i F_j \\ &= \frac{\delta_t}{\rho} (u_i^* + \frac{F_i \delta_t}{2 \rho}) F_j + \frac{\delta_t}{\rho} (u_j^* + \frac{F_j \delta_t}{2 \rho}) F_i \end{aligned}

这里令 vEDM=u+F2ρδt\bold{v}_{\rm{EDM}} = \bold{u}^{*} + \frac{\bold{F}}{2 \rho} \delta_{t} ,则 EDM 源项的展开式化简为:

Fα=δtwα[eαFcs2+12cs4(vEDMF+FvEDM):(eαeαcs2[I])](8)F_\alpha = \delta_t w_\alpha \left[ \frac{\bold{e}_\alpha \cdot \bold{F}}{c_s^2} + \frac{1}{2 c_s^4} \left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \tag{8}

上式即为 EDM 源项的展开结果。

源项误差的 Chapman-Enskog 分析

这里的推导主要参考 Li 等[4:2]的文章,并补充了一些自己的笔记。

需要提醒的是,下文的速度记号和前文不同。

源项的各阶速度矩

根据 LBM 格子张量对称性,有前 5 阶表达式:

αwα=1,αwαeαi=0,αwαeαieαj=cs2δij,αwαeαieαjeαk=0αwαeαieαjeαkeαl=cs4(δijδkl+δikδjl+δilδkj)αwαeαieαjeαkeαleαm=0\sum_{\alpha} w_{\alpha} = 1 ,\quad \sum_{\alpha} w_{\alpha} e_{\alpha i} = 0 ,\\ \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} = c_s^2 \delta_{ij} ,\quad \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} e_{\alpha k} = 0 \\ \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} e_{\alpha k} e_{\alpha l} = c_s^4 (\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{kj}) \\ \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} e_{\alpha k} e_{\alpha l} e_{\alpha m} = 0

其中 eαie_{\alpha i} 表示 D 维向量 eα\bold{e}_{\alpha} 的第 ii 个分量, wαw_\alphaeα\bold{e}_{\alpha} 的权重系数。 δij\delta_{ij} 为克罗内克符号。

所以,对于 EDM 格式的源项 Fα,EDMF_{\alpha,\mathrm{EDM}} ,其各阶速度矩为:

αFα,EDM=0,αeαFα,EDM=F,αeαeαFα,EDM=δt((uF+Fu)+δtFFρ)(9)\sum_{\alpha} F_{\alpha,\mathrm{EDM}} = 0 ,\quad \sum_{\alpha} \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} = \bold{F} ,\quad\\ \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} = \delta_t \left( (\bold{u}^* \bold{F} + \bold{F}\bold{u}^*) + \delta_t \frac{\bold{FF}}{\rho} \right) \tag{9}

技术细节可见附录(A)。

LBM 的常规 Chapman-Enskog 展开

令展开参数 ϵ\epsilon 是与克努森数同阶的小量,则将分布函数展开为

fα=fα(0)+ϵfα(1)+ϵ2fα(2)+...,f_{\alpha} = f_{\alpha}^{(0)} + \epsilon f_{\alpha}^{(1)} + \epsilon^2 f_{\alpha}^{(2)} + ...,

源项展开为 Fα=ϵF1αF_{\alpha} = \epsilon F_{1 \alpha},时间和空间偏导展开为

t=ϵt1+ϵ2t2,x=ϵx1.\frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_1} + \epsilon^2 \frac{\partial}{\partial t_2} ,\quad \frac{\partial}{\partial \bold{x}} = \epsilon \frac{\partial}{\partial \bold{x}_1}.

将式(4)中的 fα(x+eαδt,t+δt)f_\alpha (\bold{x} + \bold{e}_\alpha \delta_t, t+\delta_t)fα(x,t)f_\alpha (\bold{x}, t) 处进行 Taylor 展开(至2阶项),并把上述展开式代入,可得 ϵ\epsilon 各阶项的方程为:

O(ϵ0):fα(0)=fα(eq),(101)O(ϵ1):D1αfα(0)=1τδtfα(1)+F1αδt,(102)O(ϵ2):fα(0)t2+D1αfα(1)+δt2D1α2fα(0)=1τδtfα(2).(103)\begin{aligned} O(\epsilon^0):\quad & f_{\alpha}^{(0)} = f_{\alpha}^{(eq)} , &(10-1)\\ O(\epsilon^1):\quad & \mathrm{D}_{1 \alpha} f_{\alpha}^{(0)} = -\frac{1}{\tau \delta_t} f_{\alpha}^{(1)} + \frac{F_{1 \alpha}}{\delta_t} , &(10-2)\\ O(\epsilon^2):\quad & \frac{\partial f_{\alpha}^{(0)}}{\partial t_2} + \mathrm{D}_{1 \alpha} f_{\alpha}^{(1)} + \frac{\delta_t}{2} \mathrm{D}_{1 \alpha}^2 f_{\alpha}^{(0)} = -\frac{1}{\tau \delta_t} f_{\alpha}^{(2)} . &(10-3)\\ \end{aligned}

其中 D1α=t1+eαx1\displaystyle \mathrm{D}_{1 \alpha} = \frac{\partial}{\partial t_1} + \bold{e}_{\alpha} \cdot \frac{\partial}{\partial \bold{x}_1}。将式(10-2)代入式(10-3)中,式(10-3)可被化简为式(10-4):

O(ϵ2):fα(0)t2+D1α(112τ)fα(1)+12D1αF1α=1τδtfα(2).(104)\begin{aligned} O(\epsilon^2):\quad & \frac{\partial f_{\alpha}^{(0)}}{\partial t_2} + \mathrm{D}_{1 \alpha} \left(1 - \frac{1}{2 \tau}\right) f_{\alpha}^{(1)} + \frac{1}{2} \mathrm{D}_{1 \alpha} F_{1 \alpha} = -\frac{1}{\tau \delta_t} f_{\alpha}^{(2)} . &(10-4) \end{aligned}

(1)考虑对式(10-3)和式(10-4)计算零阶矩,得:

ρt1+1(ρu)=0ρt2+121(δtF1)=0\frac{\partial \rho}{\partial t_1} + \nabla_1 \cdot (\rho \bold{u}^*) = 0 \\ \frac{\partial \rho}{\partial t_2} + \frac{1}{2} \nabla_1 \cdot (\delta_t \bold{F}_1) = 0 \\

此处 u=(ieαfα)/(ifα)\bold{u}^* = (\sum_{i} \bold{e}_\alpha f_\alpha) / (\sum_{i} f_\alpha) 并非真实流场速度,而 u^=u+Fδt/(2ρ)\hat{\bold{u}} = \bold{u}^* + \bold{F} \delta_t / (2 \rho) 为真实流场速度(也就是式(6)给出的形式)。

基于前面假设的偏微分关系式,我们可以得到 ρt+(ρu^)=0\displaystyle \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \hat{\bold{u}}) = 0

(2)同理,对式(10-3)和式(10-4)计算一阶矩,得:

ρut1+1(ρuu)=1p+F1ρut2+1[(112τ)αeαeαfα(1)]+δt2F1t1+121(αeαeαF1α)=0\frac{\partial \rho \bold{u}^*}{\partial t_1} + \nabla_1 \cdot (\rho \bold{u}^* \bold{u}^*) = -\nabla_1 p + \bold{F}_1 \\ \frac{\partial \rho \bold{u}^*}{\partial t_2} + \nabla_1 \cdot \left[\left(1 - \frac{1}{2 \tau}\right) \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(1)} \right] + \frac{\delta_t}{2} \frac{\partial \bold{F}_1}{\partial t_1} + \frac{1}{2} \nabla_1 \cdot \left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{1 \alpha} \right) = 0

其中 p=ρcs2p = \rho c_s^2

根据式(10-2),有:

αeαeαfα(1)=τ[δtt1(αeαeαfα(0))+δt1(αeαeαeαfα(0))αeαeαF1α]\begin{aligned} \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(1)} &= -\tau \left[ \delta_t \frac{\partial}{\partial t_1}\left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) \right. +\\ &\quad \left. \delta_t \nabla_1 \cdot \left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) - \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{1 \alpha} \right] \end{aligned}

Π(1)=δt(τ12)[t1(αeαeαfα(0))+1(αeαeαeαfα(0))]=δt(τ12){ρcs2[1u+(1u)T]+uF1+F1u}\begin{aligned} \Pi^{(1)} &= -\delta_t \left( \tau - \frac{1}{2} \right) \left[ \frac{\partial}{\partial t_1}\left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) + \nabla_1 \cdot \left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) \right] \\ &= -\delta_t \left( \tau - \frac{1}{2} \right) \{\rho c_s^2 [\nabla_1 \bold{u}^* + (\nabla_1 \bold{u}^*)^T] + \bold{u}^*\bold{F}_1 + \bold{F}_1 \bold{u}^* \} \end{aligned}

式(10-4)的一阶矩表示为: ρut2+1Π(1)+δt2F1t1+1(ταeαeαF1α)=0\displaystyle \frac{\partial \rho \bold{u}^*}{\partial t_2} + \nabla_1 \Pi^{(1)} + \frac{\delta_t}{2} \frac{\partial \bold{F}_1}{\partial t_1} + \nabla_1 \cdot \left( \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{1 \alpha} \right) = 0

那么同样根据假设的微分关系,即可还原出一阶矩对应的方程。这里直接写出使用真实流场速度 u^=u+Fδt/(2ρ)\hat{\bold{u}} = \bold{u}^* + \bold{F} \delta_t / (2 \rho) 表示的形式:

(ρu^)t+(ρu^u^)=p+F+{ρν[1u+(1u)T}{ρν[1a+(1a)T]}+δt2ϵ2Ft2+[τδt(uF+Fu)+δt2FF4ρταeαeαFα]\begin{aligned} \frac{\partial (\rho \hat{\bold{u}})}{\partial t} + \nabla \cdot (\rho \hat{\bold{u}} \hat{\bold{u}}) &= -\nabla p + \bold{F} + \nabla \cdot \{ \rho \nu [\nabla_1 \bold{u}^* + (\nabla_1 \bold{u}^*)^T \} \\ &\quad - \nabla \cdot \{ \rho \nu [\nabla_1 \bold{a} + (\nabla_1 \bold{a})^T] \} + \frac{\delta_t}{2} \epsilon^2 \frac{\partial \bold{F}}{\partial t_2} \\ &\quad+ \nabla \cdot \left[ \tau \delta_t (\bold{u^* F} + \bold{F u^*}) + \delta_t^2 \frac{\bold{FF}}{4 \rho} - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha} \right] \end{aligned}

其中 a=δtF/2ρ\bold{a} = \delta_t \bold{F} / 2 \rhoν=(τ12)cs2δt\nu = (\tau - \frac{1}{2}) c_s^2 \delta_t 为流体运动粘度。 通过对比N-S方程,可见其误差项为:

Err={ρν[1a+(1a)T]}+δt2ϵ2Ft2+[τδt(uF+Fu)+δt2FF4ρταeαeαFα]\text{Err} = - \nabla \cdot \{ \rho \nu [\nabla_1 \bold{a} + (\nabla_1 \bold{a})^T] \} + \frac{\delta_t}{2} \epsilon^2 \frac{\partial \bold{F}}{\partial t_2} + \nabla \cdot \left[ \tau \delta_t (\bold{u^* F} + \bold{F u^*}) + \delta_t^2 \frac{\bold{FF}}{4 \rho} - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha} \right]

但上述的展开方式和Wagner等[5]使用Taylor展开的结果并不一致。 Li 等[4:3] 指出这种错误的来源是:上文将 fα(0)=fα(eq)(ρ,u)f_{\alpha}^{(0)} = f_{\alpha}^{(eq)}(\rho, \bold{u}^*),从而得到了 αeαfα(1)=0\sum_{\alpha} \bold{e}_{\alpha} f_{\alpha}^{(1)} = 0。 然而实际上 fα(0)=fα(eq)(ρ,u^)f_{\alpha}^{(0)} = f_{\alpha}^{(eq)}(\rho, \hat{\bold{u}}) ,从而有 αeαfα(1)=δtF1/2\sum_{\alpha} \bold{e}_{\alpha} f_{\alpha}^{(1)} = -\delta_t \bold{F}_1 / 2

修正后的 Chapman-Enskog 展开

因此,需要将原本的 LBE 等效地修改为:

fα(x+eαδt,t+δt)fα(x,t)=F^α1τ[fα(x,t)fα(eq)(ρ,u^)](11)f_\alpha (\bold{x} + \bold{e}_\alpha \delta_t, t+\delta_t) - f_\alpha (\bold{x}, t) = \hat{F}_{\alpha} -\frac{1}{\tau} \left[ f_\alpha (\bold{x}, t) - f_\alpha^{(eq)} (\rho, \hat{\bold{u}}) \right] \tag{11}

其中 F^α=Fα1τ(fα(eq)(ρ,u^)fα(eq)(ρ,u))\hat{F}_{\alpha} = F_{\alpha} - \frac{1}{\tau} (f_\alpha^{(eq)} (\rho, \hat{\bold{u}}) - f_\alpha^{(eq)} (\rho, \bold{u}^*)) 为基于原始源项 FαF_{\alpha} 修改而来的等效源项。并且需要强调的是,式(11)中的平衡态采用真实流场速度 u^\hat{\bold{u}} 计算。 在式(11)这一等效公式上,执行与上一节一致的 Chapman-Enskog 展开后,可得到还原的N-S方程为:

(ρu^)t+(ρu^u^)=p+F+{ρν[1u^+(1u^)T}+[(τ12)δt(u^F+Fu^)ταeαeαF^α](12)\begin{aligned} \frac{\partial (\rho \hat{\bold{u}})}{\partial t} + \nabla \cdot (\rho \hat{\bold{u}} \hat{\bold{u}}) &= -\nabla p + \bold{F} + \nabla \cdot \{ \rho \nu [\nabla_1 \hat{\bold{u}} + (\nabla_1 \hat{\bold{u}})^T \} \\ &\quad + \nabla \cdot \left[ \left(\tau - \frac{1}{2}\right) \delta_t (\hat{\bold{u}}\bold{F} + \bold{F}\hat{\bold{u}}) - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \hat{F}_{\alpha} \right] \end{aligned} \tag{12}

式(12)的误差项为:

Err=[(τ12)δt(u^F+Fu^)ταeαeαF^α](13)\text{Err} = \nabla \cdot \left[ \left(\tau - \frac{1}{2}\right) \delta_t (\hat{\bold{u}}\bold{F} + \bold{F}\hat{\bold{u}}) - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \hat{F}_{\alpha} \right] \tag{13}

将 EDM 的源项代入式(13)中,得到: Err=δt2(FF4ρ)\text{Err}=-\delta_t^2 \nabla\cdot\left( \frac{\bold{FF}}{4 \rho} \right)

附录

(A)源项各阶速度矩的推导

这里仅以 αeαeαFα,EDM\displaystyle \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} 的计算举例来简单说明。

αeαeαFα,EDM=αeαeαδtwα[eαFcs2+12cs4(vEDMF+FvEDM):(eαeαcs2[I])]=δtcs2α[eαeα(wαeαF)]+δt2cs4αwαeαeα[(vEDMF+FvEDM):(eαeαcs2[I])]\begin{aligned} \sum_{\alpha} \bold{e}_{\alpha} \bold{e}_{\alpha} F_{\alpha,\mathrm{EDM}} &= \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \cdot \delta_t w_\alpha \left[ \frac{\bold{e}_\alpha \cdot \bold{F}}{c_s^2} + \right.\\ &\quad\left. \frac{1}{2 c_s^4} \left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \\ &= \frac{\delta_t}{c_s^2} \sum_{\alpha} \left[ \bold{e}_\alpha \bold{e}_\alpha (w_\alpha \bold{e}_\alpha \cdot \bold{F}) \right] + \\ &\quad \frac{\delta_t}{2 c_s^4} \sum_{\alpha} w_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \cdot [\left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}])] \end{aligned}

一方面,对于 α[eαeα(wαeαF)]\displaystyle \sum_{\alpha} \left[ \bold{e}_\alpha \bold{e}_\alpha (w_\alpha \bold{e}_\alpha \cdot \bold{F}) \right] 这个矩阵,它的第 ii 行第 jj 列写作

kFkα(wαeαieαjeαk)=kFk0=0\sum_{k} F_{k} \sum_{\alpha}(w_{\alpha} e_{\alpha i} e_{\alpha j} e_{\alpha k}) = \sum_{k} F_{k} \cdot 0 = 0

另一方面,已知 vEDM=u+F2ρδt\bold{v}_{\rm{EDM}} = \bold{u}^{*} + \frac{\bold{F}}{2 \rho} \delta_{t},由于 (vEDMF+FvEDM)\left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) 为矩阵。所以,下文为方便起见,我们记

[B]=vEDMF+FvEDM=uF+Fu+δtρFF[\bold{\mathbb{B}}] = \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} = \bold{u}^{*} \bold{F} + \bold{F} \bold{u}^{*} + \frac{\delta_t}{\rho}\bold{F}\bold{F}

[B]:(eαeαcs2[I])[\bold{\mathbb{B}}] : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) 是一个标量。
对于矩阵 [A]=αwαeαeα{[B]:(eαeαcs2[I])}\displaystyle [\bold{\mathbb{A}}] = \sum_{\alpha} w_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \cdot \{[\bold{\mathbb{B}}] : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}])\} ,它的第 ii 行第 jjAij\mathbb{A}_{ij} 写作

Aij=αwαeαieαj[klBkl(eαkeαlcs2δkl)]=klBklα{wαeαieαj(eαkeαlcs2δkl)}=klBklcs4(δilδjk+δikδjl)=(Bij+Bji)cs4\begin{aligned} \mathbb{A}_{ij} &= \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} \left[ \sum_{k} \sum_{l} \mathbb{B}_{kl} (e_{\alpha k} e_{\alpha l} - c_s^2 \delta_{kl}) \right] \\&= \sum_{k} \sum_{l} \mathbb{B}_{kl} \sum_{\alpha} \left\{ w_{\alpha} e_{\alpha i} e_{\alpha j} (e_{\alpha k} e_{\alpha l} - c_s^2 \delta_{kl}) \right\} \\&= \sum_{k} \sum_{l} \mathbb{B}_{kl} \cdot c_s^4 (\delta_{il} \delta_{jk} + \delta_{ik} \delta_{jl}) \\&= (\mathbb{B}_{ij} +\mathbb{B}_{ji}) c_s^4 \end{aligned}

因为 [B][\mathbb{B}] 是对称矩阵,所以 Aij=2cs4Bij\mathbb{A}_{ij} =2 c_s^4 \mathbb{B}_{ij}。也就有:

αeαeαFα,EDM=0+δt2cs4αwαeαeα[(vEDMF+FvEDM):(eαeαcs2[I])]=δt(vEDMF+FvEDM)=δt(uF+Fu+δtρFF)\begin{aligned} \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} &= \vec{\bold{0}} + \frac{\delta_t}{2 c_s^4} \sum_{\alpha} w_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \cdot [\left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}])] \\&= \delta_t \left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) \\&= \delta_t \left( \bold{u}^{*} \bold{F} + \bold{F} \bold{u}^{*} + \frac{\delta_t}{\rho}\bold{F}\bold{F} \right) \end{aligned}

综上, αeαeαFα,EDM=δt(uF+Fu+δtFFρ)\displaystyle \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} = \delta_t \left( \bold{u}^* \bold{F} + \bold{F}\bold{u}^* + \delta_t \frac{\bold{FF}}{\rho} \right) 成立。其他低阶矩均采用类似方法进行计算。当然,这也同样可以拓展至其他源项的速度矩计算。

(B)修正源项的计算

由于 a=δtF/2ρ=u^u\bold{a} = \delta_t \bold{F} / 2 \rho = \hat{\bold{u}} - \bold{u}^* ,则 ρa=δtF/2\rho \bold{a} = \delta_t \bold{F} / 2 。因此:

1τ(fα(eq)(ρ,u^)fα(eq)(ρ,u))=ρwατ[eαacs2+12cs4(u^u^uu):(eαeαcs2[I])]=ρwατ[eαacs2+12cs4(ua+au+aa):(eαeαcs2[I])]=wαδtτ[eαF2cs2+12cs4(12(uF+Fu)+δtFF4ρ):(eαeαcs2[I])]\begin{aligned} & \frac{1}{\tau} (f_\alpha^{(eq)} (\rho, \hat{\bold{u}}) - f_\alpha^{(eq)} (\rho, \bold{u}^*)) \\=& \frac{\rho w_{\alpha}}{\tau} \left[ \frac{\bold{e}_{\alpha} \cdot \bold{a}}{c_s^2} + \frac{1}{2 c_s^4} (\hat{\bold{u}}\hat{\bold{u}} - \bold{u^* u^*}):(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \\=& \frac{\rho w_{\alpha}}{\tau} \left[ \frac{\bold{e}_{\alpha} \cdot \bold{a}}{c_s^2} + \frac{1}{2 c_s^4} ( \bold{u^* a} + \bold{a u^*} + \bold{aa}):(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \\=& \frac{w_{\alpha} \delta_t}{\tau} \left[ \frac{\bold{e}_{\alpha} \cdot \bold{F}}{2 c_s^2} + \frac{1}{2 c_s^4} \left( \frac{1}{2} (\bold{u^* F} + \bold{F u^*}) + \delta_t \frac{\bold{FF}}{4 \rho} \right) : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \end{aligned}

由于 EDM 源项的表达式为:

Fα,EDM=δtwα[eαFcs2+12cs4(u^F+Fu^):(eαeαcs2[I])]F_{\alpha,\mathrm{EDM}} = \delta_t w_\alpha \left[ \frac{\bold{e}_\alpha \cdot \bold{F}}{c_s^2} + \frac{1}{2 c_s^4} \left( \hat{\bold{u}} \bold{F} + \bold{F} \hat{\bold{u}} \right) : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right]

因此,其对应修正源项为

Fα,EDM1τ(fα(eq)(ρ,u^)fα(eq)(ρ,u))=eαFcs2(112τ)+wαδt2cs4{[(12τ1)(uF+Fu)(aF+Fa)δt4ρτFF]:(eαeαcs2[I])}\begin{aligned} & F_{\alpha,\mathrm{EDM}} - \frac{1}{\tau} (f_\alpha^{(eq)} (\rho, \hat{\bold{u}}) - f_\alpha^{(eq)} (\rho, \bold{u}^*)) \\=& \frac{\bold{e}_\alpha \cdot \bold{F}}{c_s^2} (1 - \frac{1}{2 \tau}) + \\& \frac{w_{\alpha} \delta_t}{2 c_s^4} \left\{ \left[ (\frac{1}{2 \tau} - 1) (\bold{u^* F} + \bold{F u^*}) - (\bold{a F} + \bold{F a}) \right.\right. \\& \left.\left. - \frac{\delta_t}{4 \rho \tau} \bold{FF} \right] : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right\} \end{aligned}


  1. Kupershtokh, A. L. New method of incorporating a body force term into the lattice Boltzmann equation. in Proceeding of the 5th international EHD workshop 241–246 (Poitiers, France, 2004). ↩︎ ↩︎

  2. A.L. Kupershtokh, D.A. Medvedev, & D.I. Karpov (2009). On equations of state in a lattice Boltzmann method. Computers & Mathematics with Applications, 58(5), 965-974. DOI:10.1016/j.camwa.2009.02.024. ↩︎ ↩︎ ↩︎

  3. Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, & Q. Liu (2016). Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Progress in Energy and Combustion Science, 52, 62-105. DOI:10.1016/j.pecs.2015.10.001. ↩︎

  4. Li, Q., Zhou, P., & Yan, H. (2016). Revised Chapman-Enskog analysis for a class of forcing schemes in the lattice Boltzmann method. Physical Review E, 94, 043313. DOI:10.1103/PhysRevE.94.043313. ↩︎ ↩︎ ↩︎ ↩︎

  5. Wagner, A. (2006). Thermodynamic consistency of liquid-gas lattice Boltzmann simulations. Physical Review E, 74, 056703. DOI: 10.1103/PhysRevE.74.056703. ↩︎

基于这篇的内容,我用 Python 写了一份对单种颗粒系统进行粗粒化的代码。

其具体使用流程为:

  1. obj = CG_mono_3D(...): 初始化。
  2. obj.set_domain(...): 设置矩形的计算域。
  3. obj.get_CG_data(R, CG_width): 对 R 点进行粗粒化,其粗粒化的长度为 CG_width
DEM_CG_mono.py 这里使用Github Gist展示代码。如果无法查看下面内容,青科学上网。

离散单元法(Discrete Element Method,简称DEM)是颗粒流数值模拟的重要数值方法,被用于各种各样的颗粒流定量研究中。颗粒流的粗粒化(Coarse Graining),便是将颗粒流流场的大量数据转化为连续场数据的数值方法。

在开始介绍之前,先简要列出主要的参考资料:
  • Breard, E. C. P., Dufek, J., Fullard, L., & Carrara, A. (2020). The Basal Friction Coefficient of Granular Flows With and Without Excess Pore Pressure: Implications for Pyroclastic Density Currents, Water‐Rich Debris Flows, and Rock and Submarine Avalanches. Journal of Geophysical Research: Solid Earth, 125(12), e2020JB020203. DOI:10.1029/2020JB020203.
  • Lacaze, L., & Kerswell, R. R. (2009). Axisymmetric Granular Collapse: A Transient 3D Flow Test of Viscoplasticity. Physical Review Letters, 102(10), 108305. DOI:10.1103/PhysRevLett.102.108305.
  • Weinhart, T., Labra, C., Luding, S., & Ooi, J. Y. (2016). Influence of coarse‐graining parameters on the analysis of DEM simulations of silo flow. Powder Technology, 293, 138–148. DOI:10.1016/j.powtec.2015.11.052

单种颗粒的粗粒化

为简单起见,先从仅存在单种球形颗粒的颗粒流开始说起,并且只考虑获取空间内一点 r\boldsymbol{r}tt 时刻的粗粒化数据。

粗粒化函数

在这里,引入高斯函数

W(r)={1Vwexp(r22w),r<c0,otherwiseW(\boldsymbol{r}) = \begin{cases} \frac{1}{V_w} \exp{\left(- \frac{\vert\boldsymbol{r}\vert^2}{2w} \right)} ,& {\vert\boldsymbol{r}\vert} < c \\ 0 ,& \text{otherwise} \end{cases}

其中: ww 为粗粒化宽度; c=3wc = 3 w 为粗粒化截断长度; VwV_w 是确保其密度积分等于总质量的系数

Vw=22w3π32erf(c22w)4cw2πexp(c22w2)V_w = 2 \sqrt{2} \cdot w^3 \pi^{\frac{3}{2}} \mathrm{erf}\left( \frac{c \sqrt{2}}{2 w} \right) - 4 c w^2 \pi \exp{\left( \frac{- c^2}{2 w^2} \right)}

Weinhart 等建议 ww 的取值范围应在 [0.75dmean,1.25dmean][0.75 d_{\mathrm{mean}}, 1.25 d_{\mathrm{mean}}] 之间,其中 dmeand_{\mathrm{mean}} 为平均颗粒直径。

密度和动量的计算

tt 时刻,单种颗粒的集合为 qq 。对于其中的第 ii 个球,其位置矢量为 ri\boldsymbol{r}_i ,质量为 mim_i ,速度为 ui\boldsymbol{u}_i

密度的粗粒化表达式为:

ρ(r,t)=iqmiW(rri(t))\rho(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))

动量的粗粒化表达式为:

j(r,t)=iqmiuiW(rri(t))\boldsymbol{j}(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i \boldsymbol{u}_i W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))

所以,速度的计算需要通过 j\boldsymbol{j}ρ\rho 完成:

u(r,t)=j(r,t)/ρ(r,t)\boldsymbol{u}(\boldsymbol{r}, t) = \boldsymbol{j}(\boldsymbol{r}, t) / \rho(\boldsymbol{r}, t)

根据 W(r)W(\boldsymbol{r}) 的定义,当 rri(t)c\vert \boldsymbol{r} - \boldsymbol{r}_{i}(t) \vert \ge c 时, W(rri(t))=0W(\boldsymbol{r} - \boldsymbol{r}_{i}(t)) = 0
因此在实际计算中可根据这一点跳过某些球的计算。

应力张量的计算

应力张量 [σ]=σαβ[\boldsymbol{\sigma}] = \sigma_{\alpha\beta} 被表示为由接触作用产生的部分 σαβ(c)\sigma_{\alpha\beta}^{(c)} ,以及由动能作用产生的部分 σαβ(k)\sigma_{\alpha\beta}^{(k)} 。即:

σαβ=σαβ(k)+σαβ(c)\sigma_{\alpha\beta} = \sigma_{\alpha\beta}^{(k)} + \sigma_{\alpha\beta}^{(c)}

对于 σαβ(k)\sigma_{\alpha\beta}^{(k)} ,有:

σαβ(k)(r,t)=iqmivi,αvi,βW(rri(t))\sigma_{\alpha\beta}^{(k)}(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i v_{i,\alpha}^{'} v_{i,\beta}^{'} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))

其中 vi,α=ui,αuα(r,t)v_{i,\alpha}^{'} = u_{i,\alpha} - u_{\alpha}(\boldsymbol{r}, t) 为第 ii 个球的速度波动值。

对于 σαβ(c)\sigma_{\alpha\beta}^{(c)} ,考虑第 ii 个球的所有接触:

σαβ(c)(r,t)=iqjQˉjiFij,αlij,β01W(rri(t)+slij)ds\sigma_{\alpha\beta}^{(c)}(\boldsymbol{r}, t) = \sum\limits_{i \in q} \sum\limits_{j \in \bar{Q}}^{j \neq i} F_{ij,\alpha} l_{ij,\beta} \int_{0}^{1} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t) + s \boldsymbol{l}_{ij}) \mathrm{d}s

式中

  • jj 是整个颗粒系统 Qˉ\bar{Q} 中的一个球(jij \neq i),并且球 ii 与球 jj 发生接触;
  • Fij\boldsymbol{F}_{ij} 是球 ii 所受到的来自球 jj 接触的力, Fij,αF_{ij,\alpha} 为其 α\alpha 方向分量;
  • lij\boldsymbol{l}_{ij} 是球 ii 与球 jj 的接触矢量(从 iijj), lij,αl_{ij,\alpha} 为其 α\alpha 方向分量。

在实际操作中,对于

01W(rri(t)+slij)ds\int_{0}^{1} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t) + s \boldsymbol{l}_{ij}) \mathrm{d}s

我选择使用复化 Simpson 积分计算其数值,以简化操作。

多种颗粒的粗粒化

假设一个颗粒流系统 Qˉ\bar{Q} 中有 QQ 类颗粒,对于第 qq 类颗粒 qQq \in Q, 可以使用单类颗粒粗粒化的公式计算出:第 qq 类颗粒的粗粒化密度 ρq(r,t)\rho^{q}(\boldsymbol{r}, t) 、速度 uq(r,t)\boldsymbol{u}^{q}(\boldsymbol{r}, t) 和应力张量 σαβq(r,t)\sigma_{\alpha\beta}^{q}(\boldsymbol{r}, t)

所以整个颗粒流系统 Qˉ\bar{Q} 的粗粒化结果为:

ρ(r,t)=qQρq(r,t)u(r,t)=qQuq(r,t)σαβ(r,t)=qQσαβq(r,t)\begin{aligned} \rho(\boldsymbol{r}, t) &= \sum_{q \in Q} \rho^{q}(\boldsymbol{r}, t) \\ \boldsymbol{u}(\boldsymbol{r}, t) &= \sum_{q \in Q} \boldsymbol{u}^{q}(\boldsymbol{r}, t) \\ \sigma_{\alpha\beta}(\boldsymbol{r}, t) &= \sum_{q \in Q} \sigma_{\alpha\beta}^{q}(\boldsymbol{r}, t) \\ \end{aligned}

使用粗粒化数据提取颗粒流流变参数

这里以 μ(I)\mu(I) 流变关系为例子,以上文的文献作为参照。

记速度剪切率张量 γij=uixj+ujxi\gamma_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} ,其中 uiu_i 源自粗粒化速度场。令:

γijd=γijδij3tr([γ])\gamma_{ij}^{d} = \gamma_{ij} - \frac{\delta_{ij}}{3} \mathrm{tr}([\gamma])

其中 tr([γ])=γxx+γyy+γzz\mathrm{tr}([\gamma]) = \gamma_{xx} + \gamma_{yy} + \gamma_{zz}[γ][\gamma] 的迹。

所以剪切应变率定义为:

γ˙=12γijdγijd|\dot{\gamma}| = \sqrt{\frac{1}{2} {\gamma}_{ij}^{d} {\gamma}_{ij}^{d} }

定义压强为 p=13tr([σ])=(σxx+σyy+σzz)/3p = \frac{1}{3} \mathrm{tr}([\sigma]) = (\sigma_{xx} + \sigma_{yy} + \sigma_{zz}) / 3,偏应力张量为: [τ]=τij=σijpδij[\tau] = \tau_{ij} = \sigma_{ij} - p \cdot \delta_{ij}, 所以切应力为:

τ=12τijτij|\tau| = \sqrt{\frac{1}{2} \tau_{ij} \tau_{ij} }

简单的粗粒化代码示意

我用 Python 写了一份对单种颗粒系统进行粗粒化的代码,详情见此网页

我的这个 Github Page 大约花了我一天时间从完全不懂到凑活着用。
这里简单分享一下我使用 hexo + hexo-next-theme + giscus 实现这个 Github Page 的参考资料。

  1. 建站
    从零开始用Hexo+GithubPage搭建个人网站(保姆级) - 哔哩哔哩

  2. 使用 Github Action
    如何优雅的使用Github Action服务来将Hexo部署到Github Pages - CSDN博客
    GitHub Pages | Hexo

  3. Katex 渲染公式
    Hexo博客渲染KaTeX数学公式 | Feng’s Blog

  4. Giscus 评论系统
    Hexo 使用 giscus 评论系统 | 404频道
    Hexo 评论系统(Giscus)

注:该文章同步发在知乎专栏CSDN
这里所说的 “简化的格子 Boltzmann 方法” 其实是 Simplified and Highly Stable Lattice Boltzmann Method (SHSLBM)([1], [2], [3], [4])。该方法的特色是仅使用平衡态分布函数、宏观密度、宏观速度进行 LBM 计算,并保证无条件稳定。
下面的内容只是我对相关文献进行简单阅读后的部分整理

单松弛格子 Boltzmann 方法

基本方程

针对流场的 LBGK 演化方程为

fα(r+cα+Δt,t+Δt)fα(r,t)=1τ[fα(r,t)fα(eq)](1)f_{\alpha} (\boldsymbol{r} + \boldsymbol{c}_\alpha + \Delta t, t + \Delta t) - f_{\alpha} (\boldsymbol{r}, t) = -\frac{1}{\tau} \left[ f_{\alpha} (\boldsymbol{r}, t) - f_{\alpha}^{(eq)} \right] \tag{1}

其中 fα(r,t)f_{\alpha} (\boldsymbol{r}, t)tt 时刻在 r\boldsymbol{r} 处沿着 cα\boldsymbol{c}_\alpha 方向的密度分布函数, τ\tau 为 BGK 碰撞项的松弛时间。 fα(eq)f_{\alpha}^{(eq)}cα\boldsymbol{c}_\alpha 方向的平衡态分布

fα(eq)=ρwα[1+cαucs2+(cαu)22cs4uu2cs2](2)f_{\alpha}^{(eq)} = \rho w_\alpha \left[ 1 + \frac{\boldsymbol{c}_\alpha \cdot \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_\alpha \cdot \boldsymbol{u})^2}{2 c_s^4} - \frac{\boldsymbol{u} \cdot \boldsymbol{u}}{2 c_s^2} \right] \tag{2}

式 (2) 中 ρ\rho 为流体密度, u=(ux,uy,uz)T\boldsymbol{u} = (u_x, u_y, u_z)^T 为流体宏观速度, wαw_\alphacα\boldsymbol{c}_\alpha 的权重系数, csc_s 为格子声速。 wαw_\alphacα\boldsymbol{c}_\alphacsc_s 的数值由选取的离散速度模型决定。

对于 LBM 而言,宏观量表示为

ρ=αfα=αfα(eq),ρu=αcαfα=αcαfα(eq)(3)\rho = \sum_{\alpha} f_{\alpha} = \sum_{\alpha} f_{\alpha}^{(eq)}, \quad \rho\boldsymbol{u} = \sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha} = \sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \tag{3}

Chapman-Enskog 展开

Boltzmann 方程和 LBGK 方程都可以通过多尺度展开还原为 Navier-Stokes 方程。这里以 Chapman-Enskog 为例:

fα=fα(eq)+ϵfα(1)+ϵ2fα(2)+...,t=ϵt0+ϵ2t1,r=ϵr0f_{\alpha} = f_{\alpha}^{(eq)} + \epsilon f_{\alpha}^{(1)} + \epsilon^2 f_{\alpha}^{(2)} + ... , \quad \frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_0} + \epsilon^2 \frac{\partial}{\partial t_1}, \quad \frac{\partial}{\partial \boldsymbol{r}} = \epsilon \frac{\partial}{\partial \boldsymbol{r}_0}

上述展开保留至 ϵ1\epsilon^{1} ,即可还原为[5]

ρt+(αcαfα(eq))=0(4)\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \right) = 0 \tag{4}

ρut+Π=0(5)\frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \Pi = \boldsymbol{0} \tag{5}

式 (5) 中张量 Π\Pi 表示为

Πβγ=α(cα)β(cα)γ[fα(eq)+(112τ)ϵfα(1)]\Pi_{\beta \gamma} = \sum_{\alpha} (\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} \left[ f_{\alpha}^{(eq)} + \left(1 - \frac{1}{2 \tau}\right) \epsilon f_{\alpha}^{(1)} \right]

由于只保留到 ϵ1\epsilon^{1} ,所以 ϵfα(1)fαfα(eq)=fα(neq)\epsilon f_{\alpha}^{(1)} \approx f_{\alpha} - f_{\alpha}^{(eq)} = f_{\alpha}^{(neq)} 近似为非平衡态分布函数。此时

Πβγ=α(cα)β(cα)γ[fα(eq)+(112τ)fα(neq)](6)\Pi_{\beta \gamma} = \sum_{\alpha} (\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} \left[ f_{\alpha}^{(eq)} + \left(1 - \frac{1}{2 \tau}\right) f_{\alpha}^{(neq)} \right] \tag{6}

并且存在如下关系

fα(neq)ϵfα(1)=τ[t+cα]fα(eq)Δt(7)f_{\alpha}^{(neq)} \approx \epsilon f_{\alpha}^{(1)} = -\tau \left[ \frac{\partial}{\partial t} + \boldsymbol{c}_{\alpha} \cdot \nabla \right] f_{\alpha}^{(eq)} \Delta t \tag{7}

D=t+cα\mathrm{D} = \frac{\partial}{\partial t} + \boldsymbol{c}_{\alpha} \cdot \nabla ,则 fα(neq)τ(Δt)Dfα(eq)f_{\alpha}^{(neq)} \approx -\tau \cdot (\Delta t) \cdot \mathrm{D}f_{\alpha}^{(eq)}

流体运动粘度 ν\nu 和松弛时间 τ\tau 的关系为

ν=(τ12)cs2Δt\nu = \left( \tau - \frac{1}{2} \right) c_s^2 \Delta t

SHSLBM

式 (7) 说明:在保留至 ϵ1\epsilon^{1} 的情况下, fα(neq)f_{\alpha}^{(neq)}fα(eq)f_{\alpha}^{(eq)} 存在联系。考虑到式 (1) 中由 1τ\frac{-1}{\tau} 松弛的部分即为 fα(neq)f_{\alpha}^{(neq)} ,因此式 (1) 理应可以被写作完全由 fα(eq)f_{\alpha}^{(eq)} 进行表示的形式。这便是 SHSLBM 的核心思想。

基本算法

SHSLBM 的单个时间步迭代分为预报和校正两部分,具体写作:

  1. 预报步:计算 (r,t)(\boldsymbol{r},t) 的预估密度 ρ\rho^{*} 和预估速度 u\boldsymbol{u}^{*}

ρ=αfα(eq)(rcαΔt,tΔt),ρu=αcαfα(eq)(rcαΔt,tΔt)(8)\rho^{*} = \sum_{\alpha} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t), \quad \rho^{*} \boldsymbol{u}^{*} = \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t) \tag{8}

  1. 校正步:根据 ρ\rho^{*}u\boldsymbol{u}^{*} 得出校正值

{ρ=ρρu=ρu(11τ)αcα[fα(neq)(r+Δt2cα,tΔt2)fα(neq)(rΔt2cα,tΔt2)](9-1)\begin{cases} \rho = \rho^{*} \\ \rho \boldsymbol{u} = \rho^{*} \boldsymbol{u}^{*} - \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} [f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) - f_{\alpha}^{(neq)}(\boldsymbol{r}-\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2})] \end{cases} \tag{9-1}

式 (9-1) 中非平衡态的计算通过对式 (7) 的差分实现,即:

fα(neq)(rΔt2cα,tΔt2)=τ[fα(eq)(r,t)fα(eq)(rcαΔt,tΔt)]+O((Δt)3)f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) = -\tau \left[ f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - f_{\alpha}^{(eq)}(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t - \Delta t) \right] + O((\Delta t)^3)

同理:

fα(neq)(r+Δt2cα,tΔt2)=τ[fα(eq)(r+cαΔt,t)fα(eq)(r,tΔt)]+O((Δt)3)f_{\alpha}^{(neq)}(\boldsymbol{r} + \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) = -\tau \left[ f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) - f_{\alpha}^{(eq)}(\boldsymbol{r}, t - \Delta t) \right] + O((\Delta t)^3)

在这两条式子中, fα(eq)(r,t)f_{\alpha}^{(eq)}(\boldsymbol{r}, t)fα(eq)(r+cαΔt,t)f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) 都是未知的。参照Shu等[6]的做法,对上面两条公式中的这两个量使用预报步数据进行近似,即

fα(eq)(r,t)=fα(eq)(ρ(r,t),u(r,t))f_{\alpha}^{(eq)}(\boldsymbol{r}, t) = f_{\alpha}^{(eq)}(\rho^{*}(\boldsymbol{r}, t), \boldsymbol{u}^{*}(\boldsymbol{r}, t))

式 (9-1) 中 ρu\rho \boldsymbol{u} 的计算还可以进一步化简。注意到,对于 fα(neq)(rΔt2cα,tΔt2)f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) ,有

αcαfα(neq)(rΔt2cα,tΔt2)=τ[αcαfα(eq)(r,t)αcαfα(eq)(rcαΔt,tΔt)]+O((Δt)3)=τ(ρuρu)+O((Δt)3)=O((Δt)3)\begin{aligned} \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) &= -\tau \left[ \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \right.\\ &\quad\left. -\sum\limits_{\alpha} \boldsymbol{c}_{\alpha}f_{\alpha}^{(eq)}(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t - \Delta t)\right] + O((\Delta t)^3) \\ &= -\tau (\rho^{*} \boldsymbol{u}^{*} - \rho^{*} \boldsymbol{u}^{*}) + O((\Delta t)^3) \\ &= O((\Delta t)^3) \end{aligned}

所以式 (9-1) 可进一步化简为式 (9-2)

{ρ=ρρu=ρu(11τ)αcαfα(neq)(r+Δt2cα,tΔt2)(9-2)\begin{cases} \rho = \rho^{*} \\ \rho \boldsymbol{u} = \rho^{*} \boldsymbol{u}^{*} - \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) \end{cases} \tag{9-2}

对应的宏观流程

事实上,式 (8) 和式 (9-1), (9-2) 的预报和校正存在着对应的宏观偏微分方程。下面简要说明一下。

预报步

式 (8) 对应的宏观方程为

{ρt+(αcαfα(eq))=0ρut+[α(cα)β(cα)γfα(eq)+12τα(cα)β(cα)γfα(neq)]=0(10)\begin{cases} \frac{\partial \rho}{\partial t} + \nabla \cdot \left(\sum\limits_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \right) = 0 \\ \frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \left[ \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(eq)} + \frac{1}{2 \tau} \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)} \right] = \boldsymbol{0} \end{cases} \tag{10}

首先,对 fα(eq)(rcαΔt,tΔt)f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t)(r,t)(\boldsymbol{r},t) 进行展开,并代入式 (7) ,得到:

fα(eq)(rcαΔt,tΔt)=fα(eq)(r,t)(Δt)Dfα(eq)(r,t)+(Δt)22D2fα(eq)(r,t)+O((Δt)3)=fα(eq)(r,t)(Δt)Dfα(eq)(r,t)Δt2τDfα(neq)(r,t)+O((Δt)3)(11)\begin{aligned} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t) &= f_{\alpha}^{(eq)} (\boldsymbol{r}, t) - (\Delta t) \mathrm{D} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) \\ &\quad+ \frac{(\Delta t)^2}{2} \mathrm{D}^{2} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) + O((\Delta t)^3) \\ &= f_{\alpha}^{(eq)} (\boldsymbol{r}, t) - (\Delta t) \mathrm{D} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) \\ &\quad - \frac{\Delta t}{2 \tau} \mathrm{D} f_{\alpha}^{(neq)} (\boldsymbol{r}, t) + O((\Delta t)^3) \\ \end{aligned} \tag{11}

将式 (11) 代入式 (8) 的第一条公式,有:

ρ=αfα(eq)(r,t)(Δt)αDfα(eq)(r,t)Δt2ταDfα(neq)(r,t)+O((Δt)3)(12)\begin{aligned} \rho^{*} &= \sum_{\alpha}f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \sum_{\alpha} \mathrm{D} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \\ &\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \mathrm{D} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3) \end{aligned} \tag{12}

式 (12) 右侧第一项与 ρ\rho^{*} 相等,右侧第三项根据非平衡态的守恒律可知其为 0 。所以在消去式 (12) 中的 Δt\Delta t 后, ρ\rho^{*} 的计算与式 (10) 中第一条公式对应,并具有 O((Δt)2)O((\Delta t)^2) 精度。

将式 (11) 代入式 (8) 的第二条公式,整理得:

ρu=αcαfα(eq)(r,t)(Δt)αcαDfα(eq)(r,t)Δt2ταcαDfα(neq)(r,t)+O((Δt)3)=αcαfα(eq)(r,t)(Δt)[tαcαfα(eq)(r,t)+α(cα)β(cα)γ(fα(eq)(r,t)+12τfα(neq)(r,t))]Δt2ταtcαfα(neq)(r,t)+O((Δt)3)(13)\begin{aligned} \rho^{*} \boldsymbol{u}^{*} &= \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \sum_{\alpha} \boldsymbol{c}_{\alpha} \mathrm{D} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \\ &\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \boldsymbol{c}_{\alpha} \mathrm{D} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3) \\ &= \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \left[ \frac{\partial}{\partial t} \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) + \right. \\ &\quad \left. \nabla \cdot \sum_{\alpha} (\boldsymbol{c}_{\alpha})_\beta (\boldsymbol{c}_{\alpha})_\gamma \left( f_{\alpha}^{(eq)}(\boldsymbol{r}, t) + \frac{1}{2 \tau} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) \right) \right] \\ &\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \frac{\partial}{\partial t} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3) \end{aligned} \tag{13}

式 (13) 右侧第一项与 ρu\rho^{*} \boldsymbol{u}^{*} 相等,右侧 αtcαfα(neq)(r,t)\sum\limits_{\alpha} \frac{\partial}{\partial t} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) 项根据非平衡态的守恒律可知其为 0 。所以在消去式 (13) 中的 Δt\Delta t 后, ρu\rho^{*} \boldsymbol{u}^{*} 的计算与式 (10) 中第二条公式对应,并具有 O((Δt)2)O((\Delta t)^2) 精度。

校正步

式 (9) 对应的宏观方程为

{ρt=0ρut+[(112τ)α(cα)β(cα)γfα(neq)]=0(14)\begin{cases} \frac{\partial \rho}{\partial t} = 0 \\ \frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \left[ \left( 1 - \frac{1}{2 \tau} \right) \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)} \right] = \boldsymbol{0} \end{cases} \tag{14}

这里主要介绍一下动量方程的还原。对 fα(neq)(r±Δt2cα,tΔt2)f_{\alpha}^{(neq)}(\boldsymbol{r} \pm \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) 进行 Taylor 展开,得到

fα(neq)(r±Δt2cα,tΔt2)=fα(neq)(r,tΔt2)±Δt2(cα)fα(neq)(r,tΔt2)+(Δt)28(cα)2fα(neq)(r,tΔt2)+O((Δt)3)\begin{aligned} f_{\alpha}^{(neq)}(\boldsymbol{r} \pm \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) &= f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \\ &\quad \pm \frac{\Delta t}{2} (\boldsymbol{c}_{\alpha} \cdot \nabla) f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \\ &\quad + \frac{(\Delta t)^2}{8} (\boldsymbol{c}_{\alpha} \cdot \nabla)^2 f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) + O((\Delta t)^3) \end{aligned}

这里我们以未化简的式 (9-1) 为例。在式 (9-1) 的动量校正步中,有:

(11τ)1Δtαcα[fα(neq)(r+Δt2cα,tΔt2)fα(neq)(rΔt2cα,tΔt2)]=(11τ)αcα(cα)fα(neq)(r,tΔt2)+O((Δt)2)=[(11τ)α(cα)β(cα)γfα(neq)(r,tΔt2)]+O((Δt)2)\begin{aligned} &\quad \left(1 - \frac{1}{\tau}\right) \frac{1}{\Delta t} \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} [f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) - f_{\alpha}^{(neq)}(\boldsymbol{r}-\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2})] \\ &= \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} (\boldsymbol{c}_{\alpha} \cdot \nabla) f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) + O((\Delta t)^2) \\ &= \nabla \cdot \left[ \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \right] + O((\Delta t)^2) \\ \end{aligned}

此时联立式 (13) ,即可还原式 (14) 中的动量方程。

SHSLBM 的稳定性分析

这里令 Y=(y1,y2,y3,y4)T\mathbf{Y} = (y_1, y_2, y_3, y_4)^T。其中

y1=ρ,y2=ρux,y3=ρuy,y4=ρuzy_1 = \rho,\quad y_2 = \rho u_x ,\quad y_3 = \rho u_y ,\quad y_4 = \rho u_z

将这些符号代换到式 (2) 的 fα(eq)f_{\alpha}^{(eq)} 中,得:

fα(eq)=y1αwα+wαcs2(cαxy2α+cαyy3α+cαzy4α)+wα2y1αcs4(cαxy2α+cαyy3α+cαzy4α)2wα2y1αcs2(y2α2+y3α2+y4α2)(15)\begin{aligned} f_{\alpha}^{(eq)} &= y_{1 \alpha} w_{\alpha} + \frac{w_\alpha}{c_s^2} ( c_{\alpha x} y_{2 \alpha}+ c_{\alpha y} y_{3 \alpha} + c_{\alpha z} y_{4 \alpha} ) \\ &\quad + \frac{w_\alpha}{2 y_{1 \alpha} c_s^4} ( c_{\alpha x} y_{2 \alpha}+ c_{\alpha y} y_{3 \alpha} + c_{\alpha z} y_{4 \alpha} )^2 \\ &\quad - \frac{w_\alpha}{2 y_{1 \alpha} c_s^2} (y_{2 \alpha}^2 + y_{3 \alpha}^2 + y_{4 \alpha}^2) \end{aligned} \tag{15}

where the subscript α\alpha denotes the quantities at a space level of (rcαΔt)(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t).[2:1]

预报步

对于第 nn 步的结果 Yn\mathbf{Y}^{n},式 (8) 的预报步写作 Y=G(Yn)\mathbf{Y}^{*} = \mathbf{G}(\mathbf{Y}^{n})

为便于进行诺依曼稳定性分析,对上式进行求导,有

dY=[Yy1]ndy1n+[Yy2]ndy2n+[Yy3]ndy3n+[Yy4]ndy4n\mathrm{d \mathbf{Y}^{*}} = \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_1} \right]^n \mathrm{d}y_1^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_2} \right]^n \mathrm{d}y_2^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_3} \right]^n \mathrm{d}y_3^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_4} \right]^n \mathrm{d}y_4^n

上标 nn 代表第 nn 个时间步。由于 dY=(dy1,dy2,dy3,dy4)T\mathrm{d \mathbf{Y}^{*}} = (\mathrm{d} y_1^{*}, \mathrm{d} y_2^{*}, \mathrm{d} y_3^{*}, \mathrm{d} y_4^{*})^T ,所以,对于预报步结果的每个分量而言有: dyλ=κ=14[yλyκ]ndyκn\mathrm{d} y_{\lambda}^{*} = \sum_{\kappa = 1}^{4}\left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n \mathrm{d}y_{\kappa}^nλ,κ[1,2,3,4]\lambda, \kappa \in [1,2,3,4])。

对于本文列出的方程组, [yλyκ]n\left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n 多达 16 个。因此下文仅拿 [y1y1]n\left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n 举例进行计算。

[y1y1]n=y1αfα(eq)=α[wαwα2(y1αn)2cs4(cαxy2αn+cαyy3αn+cαzy4αn)2+wα2(y1αn)2cs2((y2αn)2+(y3αn)2+(y4αn)2)](16)\begin{aligned} \left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n &= \frac{\partial}{\partial y_{1}} \sum_{\alpha} f_{\alpha}^{(eq)} \\ &= \sum_{\alpha} \left[ w_{\alpha} - \frac{w_\alpha}{2 (y_{1\alpha}^n)^2 c_s^4} ( c_{\alpha x} y_{2 \alpha}^n + c_{\alpha y} y_{3 \alpha}^n + c_{\alpha z} y_{4 \alpha}^n )^2 \right. \\ &\quad\left. + \frac{w_\alpha}{2 (y_{1\alpha}^n)^2 c_s^2} ((y_{2 \alpha}^n)^2 + (y_{3 \alpha}^n)^2 + (y_{4 \alpha}^n)^2) \right] \end{aligned} \tag{16}

假设式 (16) 中的 dyi\mathrm{d} y^{*}_{i}, dyin\mathrm{d} y^{n}_{i}, yiαny^{n}_{i \alpha} 均可被展开为如下形式(i=[1,2,3,4]i = [1,2,3,4]j=1j=\sqrt{-1}

dyi=dyiexp(jkr),dyin=dyinexp(jkr),yiαn=yinexp(jkr)\mathrm{d} y^{*}_{i} = \overline{\mathrm{d} y_{i}}^{*} \exp(j \boldsymbol{k \cdot r}) ,\quad \mathrm{d} y^{n}_{i} = \overline{\mathrm{d} y_{i}}^{n} \exp(j \boldsymbol{k \cdot r}) ,\quad y^{n}_{i \alpha} = y_{i}^{n} \exp(j \boldsymbol{k \cdot r})

代入到式 (16) 和 dY\mathrm{d \mathbf{Y}^{*}} 表达式,得

[y1y1]n=y1αfα(eq)=α[wαwα2(y1n)2cs4(cαxy2n+cαyy3n+cαzy4n)2+wα2(y1n)2cs2((y2n)2+(y3n)2+(y4n)2)](17)\begin{aligned} \left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n &= \frac{\partial}{\partial y_{1}} \sum_{\alpha} f_{\alpha}^{(eq)} \\ &= \sum_{\alpha} \left[ w_{\alpha} - \frac{w_\alpha}{2 (y_{1}^n)^2 c_s^4} ( c_{\alpha x} y_{2}^n + c_{\alpha y} y_{3}^n + c_{\alpha z} y_{4}^n )^2 \right. \\ &\quad\left. + \frac{w_\alpha}{2 (y_{1}^n)^2 c_s^2} ((y_{2}^n)^2 + (y_{3}^n)^2 + (y_{4}^n)^2) \right] \end{aligned} \tag{17}

格子张量的一些性质[2:2]

αwα=1.α[0,1,...,q1]αwαcαβ=0.β[x,y,z]αwαcαβcαγ=cs2δβγ.β,γ[x,y,z]\begin{aligned} & \sum_{\alpha} w_{\alpha} = 1 .&\qquad \alpha \in [0,1, ..., q-1]\\ & \sum_{\alpha} w_{\alpha} c_{\alpha \beta} = 0 .&\qquad \beta \in [x,y,z] \\ & \sum_{\alpha} w_{\alpha} c_{\alpha \beta} c_{\alpha \gamma} = c_s^2 \delta_{\beta \gamma} .&\qquad \beta, \gamma \in [x,y,z] \end{aligned}

注: qq 为离散速度数目

基于格子张量的性质,式 (17) 右侧仅剩下 αwα\sum\limits_{\alpha} w_{\alpha} ,即:

[y1y1]n=1\left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n = 1

同理,其他项也可进行类似分析,并得出

[yλyκ]n=δλκ\left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n = \delta_{\lambda \kappa}

其中 δλκ\delta_{\lambda \kappa} 为狄拉克函数,且 λ,κ[1,2,3,4]\lambda,\kappa \in [1,2,3,4] 。所以我们可以得到 dY=dYn\overline{\mathrm{d} \mathbf{Y}}^{*} = \overline{\mathrm{d} \mathbf{Y}}^{n}

校正步

dyin+1=dyin+1exp(jkr){\mathrm{d} y_i}^{n+1} = \overline{\mathrm{d} y_i}^{n+1} \exp(j \boldsymbol{k \cdot r})i[1,2,3,4]i \in [1,2,3,4])。由于校正步(式 (9-2) )中并未对 ρ\rho 做修改,易得:

dy1n+1=dy1=dy1n\overline{\mathrm{d} y_1}^{n+1} = \overline{\mathrm{d} y_1}^{*} = \overline{\mathrm{d} y_1}^{n}

根据式 (9-2) 的校正步,以及 fα(neq)(r+Δt2cα,tΔt2)f_{\alpha}^{(neq)}(\boldsymbol{r} + \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) 的差分表示,有

y2n+1=y2(τ1)y2n+(τ1)αcαxfα(eq)(r+cαΔt,t)y3n+1=y3(τ1)y3n+(τ1)αcαyfα(eq)(r+cαΔt,t)y4n+1=y4(τ1)y4n+(τ1)αcαzfα(eq)(r+cαΔt,t)(18)\begin{aligned} y_2^{n+1} = y_2^{*} - (\tau - 1) y_2^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha x} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\ y_3^{n+1} = y_3^{*} - (\tau - 1) y_3^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha y} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\ y_4^{n+1} = y_4^{*} - (\tau - 1) y_4^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha z} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \end{aligned} \tag{18}

这里同样以 y2n+1y_2^{n+1} 的计算为例,记

S=αcαxfα(eq)(r+cαΔt,t)=αcαxwαy1(α)+αcαxwαcs2(cαxy2(α)+cαyy3(α)+cαzy4(α))+αcαxwα2cs4y1(α)(cαxy2(α)+cαyy3(α)+cαzy4(α))2αcαxwα2cs2y1(α)((y2(α))2+(y3(α))2+(y4(α))2)2\begin{aligned} S &= \sum\limits_{\alpha} c_{\alpha x} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\ &= \sum\limits_{\alpha} c_{\alpha x} w_{\alpha} y_{1(-\alpha)}^{*} \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{c_s^2} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*}) \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 y_{1(-\alpha)}^{*}} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*})^2 \\ &\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 y_{1(-\alpha)}^{*}} ( (y_{2(-\alpha)}^{*})^2 + (y_{3(-\alpha)}^{*})^2 + (y_{4(-\alpha)}^{*})^2 )^2 \\ \end{aligned}

其中下标 (α)(-\alpha) 代表来自 r+cαΔt\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t 的值。则

y2n+1=y2+(τ1)(Sy2n)(19)y_2^{n+1} = y_2^{*} + (\tau - 1) (S - y_2^{n}) \tag{19}

所以其差分表示为

dy2n+1=dy2+(τ1)(dSdy2n)(20)\mathrm{d} y_2^{n+1} = \mathrm{d} y_2^{*} + (\tau - 1) (\mathrm{d} S - \mathrm{d} y_2^{n}) \tag{20}

对于 dS\mathrm{d} S,有

dS=(Sy1)dy1+(Sy2)dy2+(Sy3)dy3+(Sy4)dy4\mathrm{d} S = \left( \frac{\partial S}{\partial y_1} \right)^{*} \mathrm{d}y_1^{*} + \left( \frac{\partial S}{\partial y_2} \right)^{*} \mathrm{d}y_2^{*} + \left( \frac{\partial S}{\partial y_3} \right)^{*} \mathrm{d}y_3^{*} + \left( \frac{\partial S}{\partial y_4} \right)^{*} \mathrm{d}y_4^{*}

(Sy1)\left( \frac{\partial S}{\partial y_1} \right)^{*} 为例,有

(Sy1)=αcαxwααcαxwα2cs4(y1(α))2(cαxy2(α)+cαyy3(α)+cαzy4(α))2+αcαxwα2cs2(y1(α))2((y2(α))2+(y3(α))2+(y4(α))2)2\begin{aligned} \left( \frac{\partial S}{\partial y_1} \right)^{*} &= \sum\limits_{\alpha} c_{\alpha x} w_{\alpha} \\ &\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 (y_{1(-\alpha)}^{*})^2} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*})^2 \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 (y_{1(-\alpha)}^{*})^2} ( (y_{2(-\alpha)}^{*})^2 + (y_{3(-\alpha)}^{*})^2 + (y_{4(-\alpha)}^{*})^2 )^2 \end{aligned}

与预报步的做法类似,将 yi(α)y_{i(-\alpha)}^{*} 近似为 yiy_{i}^{*},得到:

(Sy1)=αcαxwααcαxwα2cs4(y1)2(cαxy2+cαyy3+cαzy4)2+αcαxwα2cs2(y1)2((y2)2+(y3)2+(y4)2)2\begin{aligned} \left( \frac{\partial S}{\partial y_1} \right)^{*} &= \sum\limits_{\alpha} c_{\alpha x} w_{\alpha} \\ &\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 (y_{1}^{*})^2} (c_{\alpha x} y_{2}^{*} + c_{\alpha y} y_{3}^{*} + c_{\alpha z} y_{4}^{*})^2 \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 (y_{1}^{*})^2} ( (y_{2}^{*})^2 + (y_{3}^{*})^2 + (y_{4}^{*})^2 )^2 \end{aligned}

此时所有 yiy_{i}^{*} 均可被视为常数提取出来。根据格子张量的性质,可得: (Sy1)=0\left( \frac{\partial S}{\partial y_1} \right)^{*} = 0

同理

(Sy2)=1,(Sy3)=(Sy4)=0\left( \frac{\partial S}{\partial y_2} \right)^{*} = 1,\quad \left( \frac{\partial S}{\partial y_3} \right)^{*} = \left( \frac{\partial S}{\partial y_4} \right)^{*} = 0

此时 dS=dy2\overline{\mathrm{d} S } = \overline{\mathrm{d} y_2}^{*} ,则式 (20) 可化简为:

dy2n+1=dy2=dy2n\overline{\mathrm{d} y_2}^{n+1} = \overline{\mathrm{d} y_2}^{*} = \overline{\mathrm{d} y_2}^{n}

同理,对 y3n+1y_3^{n+1}y4n+1y_4^{n+1} 计算可得: dy3n+1=dy3n\overline{\mathrm{d} y_3}^{n+1} = \overline{\mathrm{d} y_3}^{n}dy4n+1=dy4n\overline{\mathrm{d} y_4}^{n+1} = \overline{\mathrm{d} y_4}^{n}

dYn=(dy1n,dy2n,dy3n,dy4n)T\mathrm{d \mathbf{Y}^{n}} = (\mathrm{d} y_1^{n}, \mathrm{d} y_2^{n}, \mathrm{d} y_3^{n}, \mathrm{d} y_4^{n})^TdYn+1=(dy1n+1,dy2n+1,dy3n+1,dy4n+1)T\mathrm{d \mathbf{Y}^{n+1}} = (\mathrm{d} y_1^{n+1}, \mathrm{d} y_2^{n+1}, \mathrm{d} y_3^{n+1}, \mathrm{d} y_4^{n+1})^T ,那么:

dYn+1=IdYn\overline{\mathrm{d \mathbf{Y}}}^{n+1} = \mathbf{I} \, \overline{\mathrm{d \mathbf{Y}}}^{n}

其中 I\mathbf{I} 为单位矩阵。显然,特征矩阵的特征值全是 1 ,代表着该格式无条件稳定。


  1. Chen, Z., Shu, C., Wang, Y., Yang, L. M., & Tan, D. (2017). A Simplified Lattice Boltzmann Method without Evolution of Distribution Function. Advances in Applied Mathematics and Mechanics, 9(1), 1–22. DOI:10.4208/aamm.OA-2016-0029 ↩︎

  2. Z. Chen, C. Shu, D. Tan; Three-dimensional simplified and unconditionally stable lattice Boltzmann method for incompressible isothermal and thermal flows. Physics of Fluids. 1 May 2017; 29 (5): 053601. DOI:10.1063/1.4983339 ↩︎ ↩︎ ↩︎

  3. Chen Z, Shu C, Tan D, Wu C. On improvements of simplified and highly stable lattice Boltzmann method: Formulations, boundary treatment, and stability analysis. Int J Numer Meth Fluids. 2018; 87: 161–179. DOI:10.1002/fld.4485 ↩︎

  4. Chen, Z., & Shu, C. (2020). Simplified and Highly Stable Lattice Boltzmann Method. WORLD SCIENTIFIC. DOI:10.1142/12047 ↩︎

  5. 本文中所有 \nabla 均表示空间偏导 r\frac{\partial}{\partial \boldsymbol{r}}↩︎

  6. C. Shu, Y. Wang, C. Teo, and J. Wu. Development of lattice Boltzmann flux solver for simulation of incompressible flows. Advances In Applied Mathematics And Mechanics. 6, 436 (2014). ↩︎

这是我之前的一次期末大报告作业,里面的内容也算是从其他文章里搬运的,在此重新誊抄出来。内容如有错误,欢迎各位指正。

[注] 1. 下文提到的 xxuuξ\xi 无论有没有加粗均为向量。
2. 本文已经同步发表在本人的知乎CSDN上,可在其他网站获得同样的阅读效果。

一、 单松弛格子Boltzmann方程

1.1 Boltzmann方程

格子Boltzmann方法基于Boltzmann方程(即(1)式),该方程引入速度分布函数,描述了非热力学平衡状态的热力学系统统计行为:

ft+ξfx+Fmfξ=Ω(f)(1)\frac{\partial f}{\partial t} + \mathbf{\xi} \cdot \frac{\partial f}{\partial \boldsymbol{x}} + \frac{F}{m} \cdot \frac{\partial f}{\partial \mathbf{\xi}}=\Omega(f) \tag{1}

(1)式中 Ω(f)\Omega(f) 为碰撞项, ξ=xt\mathbf{\xi}=\frac{\partial x}{\partial t} 。速度分布函数的定义是:在时刻tt,以xx为中心的空间微元 dx\mathrm{d}\mathit{x} 内,速度在 ξ\xi(ξ+dξ)(\xi+\mathrm{d}\it{\xi}) 之间的分子的数量为 fdxdξf \mathrm{d}\it{x} \mathrm{d}\it{\xi} 。记分子质量为 mm ,流体宏观速度为uu,单位体积内能为ee,分子数密度为nn,则:

ρ=mn=mfdξ\rho = mn = m \int{f} \mathrm{d}\xi

ρu=mξfdξ\rho u = m \int{\xi f} \mathrm{d}\xi

ρE=ρe+12ρu2=mfξ22dξ\rho E = \rho e + \frac{1}{2} \rho u^2 = m \int{ f \frac{\xi^2}{2} } \mathrm{d} \xi

碰撞项 描述了分子碰撞对整个系统运动产生的影响,它具有两个特点:①满足质量、动量和能量守恒;②能够反映系统趋于平衡态的趋势。Bhatnagur、Gross和Krook三人在此基础上提出了最为简单且著名的BGK模型,即(2)式。

Ω(f)=1τc(ff(eq))(2)\Omega(f) = -\frac{1}{\tau_c} \left ( f - f^{(eq)} \right ) \tag{2}

其中τc\tau_c为松弛时间, f(eq)f^{(eq)} 为平衡态分布函数。 可根据Boltzmann H定理计算得:

f(eq)=n(2πRT)3/2exp((ξu)22RT)(3)f^{(eq)} = \frac{n}{(2 \pi RT)^{3/2}} \mathrm{exp}\left( -\frac{(\mathbf{\xi-u})^2}{2RT} \right) \tag{3}

1.2 格子Boltzmann方程

为了将Boltzmann方程应用于实际的数值模拟,需要将流体视为大量的离散粒子。离散粒子于流体分子不同,它们驻留在网格上,并按一定的规则在网格上发生碰撞和迁移。因此,1.2节将 mfmf 记为 ff,此时 ff 为密度分布函数。对密度分布函数 ff ,不考虑 FF 的Boltzmann-BGK方程为式。

ft+ξfx=1τc[ff(eq)](4)\frac{\partial f}{\partial t} + \mathbf{\xi} \cdot \frac{\partial f}{\partial \mathit{x}} = - \frac{1}{\tau_c} \left[ f - f^{(eq)} \right] \tag{4}

其宏观量按如下方式确定

\rho = \int { f } \mathrm{d} \xi ,\rho u = \int { \xi f } \mathrm{d} \xi ,\rho E = \int { f \frac{\xi^2}{2} } \mathrm{d} \xi

ff 和其对应的 f(eq)f^{(eq)} 可进行如下离散:fi=wiff_i = w_i ffi(eq)=wif(eq)f_i^{(eq)} = w_i f^{(eq)}。其中下标 ii 表示特征方向 ci\boldsymbol{c}_ifif_ifi(eq)f_i^{(eq)}均沿着ci\boldsymbol{c}_i方向运动。据此,沿着特征方向离散(4)式得到的差分格式为

fi(x+ξiΔt,t+Δt)fi(x,t)=1τ[fi(x,t)fi(eq)(x,t)](5)f_i (x+ξ_i Δt,t+Δt)-f_i (x,t)=-\frac{1}{τ} \left[f_i (x,t)-f_i^{(eq)} (x,t) \right] \tag{5}

其中 τ=τcΔt\tau = \frac{\tau_c}{\Delta t} 为无量纲松弛时间。由于对应的 f(eq)f^{(eq)} 为:

f(eq)=ρexp((ξu)22RT)(6-1)f^{(eq)} = \rho \cdot \mathrm{exp} ( - \frac{(\xi - \it{u})^2}{2\it{RT}} ) \tag{6-1}

所以将 f(eq)f^{(eq)} 展开至2阶并离散,可得:

fi(eq)=wiρ[1+ξiucs2+(ξiu)22cs4uu2cs2](6-2)f^{(eq)}_i = w_i \rho \left[ 1 + \frac{\xi_i \cdot u}{c_s^2} + \frac{(\xi_i \cdot u)^2}{2 c_s^4} - \frac{u \cdot u}{2 c_s^2} \right] \tag{6-2}

其中 cs=RTc_s = \sqrt{RT} 在LBM模拟中一般被称作格子声速。

二、Chapman-Enskog展开分析

将基本碰撞不变量 φi=1,ξ,ξ22φ_i=1,ξ,\frac{ξ^2}{2} 乘以(1)式两端,并对 ξ\xi 积分有:

ρt+(ρu)=0(7-1)\frac{∂ρ}{\partial t}+∇⋅(ρ \mathbf{u})=0 \tag{7-1}

(ρu)t+(ρuu+P)=ρa(7-2)\frac{∂(ρu)}{\partial t}+∇⋅(ρ \mathbf{uu}+P) = ρ \mathit{\bf{a}} \tag{7-2}

(ρE)t+(ρuE+Pu+Q)=ρau(7-3)\frac{∂(\rho E)}{\partial t} + ∇⋅(\rho \mathbf{u}E + P⋅\mathbf{u} + \mathbf{Q})=\rho \mathbf{a⋅u} \tag{7-3}

Chapman-Enskog展开(下文简称C-E展开)是处理Boltzmann方程的一种分析方法,该方法引入小量ε\varepsilon作为展开因子(一般与Knudsen数同阶),将BGK方程在不同的尺度上展开,即:

f=n=0εnf(n)(8)f = \sum_{n=0}^{\infty} \varepsilon^n f^{(n)} \tag{8}

并将分布函数、导数、物理量等都按照ε\varepsilon的不同阶次展开。该方法对 f(n)f^{(n)} 进行了限制,使其满足:f(n)φidξ=0\int f^{(n)} φ_i \mathrm{d}ξ=0φi(i=1,..,5)φ_i (i=1,..,5) 为基本碰撞不变量。将 ff 视为 ρ\rhou\mathbf{u}TT 及其各阶导数的泛函,而时间偏导项则展开为

t=n=0nt(9)\frac{\partial }{\partial t} = \sum_{n=0}^{\infty } \frac{\partial _n}{\partial t} \tag{9}

将(8)式和(9)式代入到(7)式中。得 O(ε0)O({\varepsilon}^0) 方程组为:

0ρt+(ρu)=0(10-1)\frac{∂_0 ρ}{\partial t}+∇⋅(ρ \mathbf{u})=0 \tag{10-1}

0(ρu)t+(ρuu+P(0))=ρa(10-2)\frac{∂_0 (ρu)}{\partial t}+∇⋅(ρ \mathbf{uu}+P^{(0)}) = ρ \mathbf{a} \tag{10-2}

0(ρE)t+(ρuE+P(0)u+Q(0))=ρau(10-3)\frac{∂_0(\rho E)}{\partial t} + ∇⋅(\rho \mathbf{u}E + P^{(0)} ⋅ \mathbf{u} + \mathbf{Q}^{(0)})=\rho \mathbf{a⋅u} \tag{10-3}

以及 O(εn)O({\varepsilon}^n) 方程组为 (n>0)(n>0)

nρt=0(11-1)\frac{∂_n ρ}{\partial t}=0 \tag{11-1}

n(ρu)t+P(n)=0(11-2)\frac{∂_n (ρu)}{\partial t}+∇ ⋅ P^{(n)} = 0 \tag{11-2}

n(ρE)t+(P(n)u+Q(n))=0(11-3)\frac{∂_n (\rho E)}{\partial t} + ∇⋅( P^{(n)} ⋅ \mathbf{u} + \mathbf{Q}^{(n)} ) = 0 \tag{11-3}

其中

C=ξu,P(n)=mCCf(n)dξ,Q(n)=mCC22f(n)dξC=\boldsymbol{ξ-u}, P^{(n)}=m∫CCf^{(n)} \mathrm{d}ξ , Q^{(n)}=m∫C \frac{C^2}{2} f^{(n)} \mathrm{d}ξ

只需取精确到O(ε1)O({\varepsilon}^1)的宏观方程组进行分析,即可还原出Naiver-Stokes方程。且具有如下关系:

P=P(0)+εP(1)=pI2εμ(STr(S)3I)P=P^{(0)} + \varepsilon P^{(1)} = p \mathbf{I} - 2 \varepsilon \mu (\mathbf{S} - \frac{\mathrm{Tr}(S)}{3} \mathbf{I})

Q=Q(0)+εQ(1)=εκTQ = Q^{(0)} + \varepsilon Q^{(1)} = - \varepsilon \kappa \nabla{T}

其中:κ\kappa为热传导系数, μ\mu 为粘性系数。

离散的LBGK方程也可采用类似的方法展开,还原出宏观流体运动方程(见附录A)。展开方式如下:

fi=fi(0)+εfi(1)+ε2fi(2)+...f_i = f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)} + ...

x=εx1,x1=εx\frac{\partial}{\partial x} = \varepsilon \frac{\partial}{\partial x_1} , x_1 = \varepsilon x

略去O(ε3)O({\varepsilon}^3)项,并对O(ε0)O({\varepsilon}^0)项、O(ε1)O({\varepsilon}^1)项和O(ε2)O({\varepsilon}^2)项求0阶和1阶速度矩,即可得到不可压缩流体的Navier-Stokes方程。由于在分析过程中略去了O(ε3)O({\varepsilon}^3)项和速度的3阶项,导致结果缺失了部分物理信息。

三、 Grad-13矩方法及其发展

3.1 Grad的矩分析方法

Grad -13矩是一种由数学家Harold Grad提出的分析方法,其基于Hermite多项式进行分析。对于一个长度为 NN 的向量 x\boldsymbol{x} ,若定义权函数

ω(x)=1(2π)N/2exp(xx2)(12)\omega( \boldsymbol{x}) = \frac{1}{(2 \pi)^{N/2}} \mathrm{exp} \left(- \frac{ \boldsymbol{x} \cdot \boldsymbol{x}}{2} \right) \tag{12}

则其对应的nn阶Hermite多项式写作

H(n)(x)=(1)nω(x)nω(x)(13)\mathcal{H} ^{(n)} ( \boldsymbol{x}) = \frac{(-1)^n}{\omega ( \boldsymbol{x})} \nabla^{n} \omega ( \boldsymbol{x}) \tag{13}

由Hermite多项式的正交性,可证明得到:

+ω(x)Hi(m)(x)Hj(n)(x)dx=δmnδij(14)∫_{-∞}^{+∞} \omega ( \boldsymbol{x}) H_i^{(m)}( \boldsymbol{x}) H_j^{(n)}( \boldsymbol{x}) \mathrm{d} \boldsymbol{x} = \delta_{mn} \delta_{ij} \tag{14}

据此,可将任意一个函数 f(x)f(x) 使用Hermite多项式进行展开[1]。即:

f(x)=ω12n=0i=1Nai(n)Hi(n)(x)(15)f( \boldsymbol{x})=ω^{\frac{1}{2}} \sum_{n=0}^∞ \sum_{i=1}^N a_i^{(n)} H_i^{(n)}( \boldsymbol{x}) \tag{15}

所以

ω1/2f(x)Hj(m)(x)dx=i=1Nai(m)δijm=m!aim(16)\int_{-\infty}^{-\infty} \omega^{1/2} f(x) \mathcal{H}_j^{(m)} (x) \mathrm{d} \it{x} = \sum_{i=\mathrm{1}}^{N} a_i^{\mathrm{(}m\mathrm{)}} \delta_{ij}^m = m\mathrm{!} a_i^{m} \mathrm{\tag{16}}

其中 ai(m)a_i^{(m)} 可借由Hermite多项式的正交性进行求解。 δijm\delta_{ij}^{m} 是一个记号,表示 mmδij\delta_{ij}的乘积之和, iijj 分别来自下标集 (i1,i2,...,im)(i_1 , i_2 , ... , i_m)(j1,j2,...,jm)(j_1 , j_2 , ... , j_m)
在此数学基础之上,Harold Grad将分布函数 ff 展开为[2]

f(x,ξ,t)=ωn=01n!a(n)H(n)(ξ)(17)f(x, \xi, t) = \omega \sum_{n=0}^{\infty} \frac{1}{n!} a^{(n)} H^{(n)}(\xi) \tag{17}

根据式(16), a(n)a^{(n)} 的表达式为:

a(n)=+f(x,ξ,t)H(n)(ξ)dξ(18)a^{(n)} = \int_{-\infty}^{+\infty} f(x, \xi, t) H^{(n)}(\xi) d\xi \tag{18}

可得 a(n)a^{(n)} 的前4项为:

a(0)=+f(x,ξ,t)dξ=ρ(19-1)a^{(0)} = \int_{-\infty}^{+\infty} f(x, \xi, t) d\xi= \rho \tag{19-1}

a(1)=+f(x,ξ,t)ξdξ=ρu(19-2)a^{(1)} = \int_{-\infty}^{+\infty} f(x, \xi, t) \xi d\xi = \rho u \tag{19-2}

a(2)=+f(x,ξ,t)(ξ2δ)dξ=ρ(u2δ)+P(19-3)a^{(2)} = \int_{-\infty}^{+\infty} f(x, \xi, t) (\xi^2 - \delta) d\xi = \rho (\mathbf{u}^2 - \mathbf{\delta}) + \mathbf{P} \tag{19-3}

a(3)=+f(x,ξ,t)(ξ3ξδ)dξ=(D1)ρu3+ua(2)+Q(19-4)a^{(3)} = \int_{-\infty}^{+\infty} f(x, \xi, t) (\xi^3 - \xi \delta) d\xi = (D-1) \rho \bf{u}^3 + \bf{u} \it{a}^{\mathrm{(2)}} \mathrm{+} \mathbf{Q} \tag{19-4}

注意到从 a(0)a^{(0)}a(3)a^{(3)} 都描述了流场的宏观量,所以通过上述表达式的值可以解出流场宏观量。与C-E分析方法不同的是分析过程中没有舍去高阶量。因此这种展开方法得到的式子合理地引入了更高阶的矩,物理意义更加明确。

3.2 Gauss-Hermite积分的引入和离散化BGK方程

Grad的思路开创了描述的新方式,而单肖文、袁学锋、陈沪东三位学者在此基础上进一步发展,于2006年在 Journal of Fluid Mechanics 上提出了将原思路拓展到离散的Boltzmann-BGK方程中的方法[3],简化了求解过程。
首先,基于式(16)将 ff 截断至前 NN 阶,得:

ffN=ωn=0N1n!a(n)H(n)(ξ)(20)f \approx f^N = \omega \sum_{n=0}^{N} \frac{1}{n!} a^{(n)} H^{(n)}(\xi) \tag{20}

fNf^Nff 有着完全相同的前 NN 阶矩。通过截断操作将BGK方程对宏观流动方程的近似保持在动力学层面而不是热力学层面。此时fNf^N可表示为:

fN(x,ξ,t)H(n)(ξ)=ω(ξ)p(x,ξ,t)(21)f^{N}(x, \xi , t) H^{(n)}(\xi) = \omega(\xi) p(x, \xi, t) \tag{21}

其中 p(x,ξ,t)p(x, \xi, t) 为阶数不超过 2N2N 的多项式。使用Gauss-Hermite积分可将 a(n)a^{(n)} 精确地展开为 p(x,ξ,t)p(x, \xi, t) 的加权和。

a(n)=ω(ξ)p(x,ξ,t)dξ=a=1dwap(x,ξa,t)=a=1dwaω(ξa)fN(x,ξa,t)H(n)(ξa)(22)a^{(n)} =\int \omega(\xi) p(x, \xi, t) d \xi = \sum_{a=1}^{d} w_a p(x, \xi_a, t) = \sum_{a=1}^{d} \frac{w_a}{\omega(\xi_a)} f^{N}(x, \xi_a, t) H^{(n)}(\xi_a) \tag{22}

式(22)中 waw_a 和 $ \boldsymbol{\xi}_a$ 是Gauss-Hermite积分的权重和特征向量( a=1,2,...,da=1,2,...,d )。所以使用 fN(x,ξa,t)f^N (x, \boldsymbol{\xi}_a, t) 表示 fNf^N 的前 NN 阶矩在数学上可行。
据此可定义离散化的分布函数 faf_aa=1,2,...,da=1,2,...,d ),其对应的权重和特征方向为 waw_aξa\xi_a

fa(x,t)=waω(ξa)f(x,ξa,t)(23)f_a (x,t) = \frac{w_a}{\omega( \boldsymbol{\xi}_a)} f(x, \boldsymbol{\xi}_a, t) \tag{23}

则宏观量可表达为:

ρ=a=1dfa,ρu=a=1dfaξa,ρuu+P=a=1dfaξaξa(24)\rho = \sum_{a=1}^{d} f_a , \rho \mathbf{u} = \sum_{a=1}^{d} f_a \mathbf{\xi}_a , \rho \mathbf{uu} + \mathbf{P} = \sum_{a=1}^{d} f_a \mathbf{\xi}_a \mathbf{\xi}_a \tag{24}

下面将(1)式的碰撞项换为(2)式,并把 Fmfξ\frac{\mathbf{F}}{m} \frac{\partial f}{\partial \xi} 直接记作外力项 F(ξ)-\mathbf{F} (\xi) ,得

ft+ξfx=1τ(ff(eq))+F(ξ)(25)\frac{\partial f}{\partial t} + \mathbf{\xi} \cdot \frac{\partial f}{\partial \mathit{x}} = \frac{-1}{\tau} ( f - f^{(eq)} ) + \mathbf{F} (\xi) \tag{25}

对式(25),令 ξ=ξa\boldsymbol{\xi} = \boldsymbol{\xi}_a ,并在两端乘上 waω(ξa)\frac{w_a}{\omega( \boldsymbol{\xi}_a)} ,即可得到使用 faf_a 离散的方程( a=1,2,...,da=1,2,...,d ):

fat+ξafa=1τ(fafa(eq))+Fa(26-1)\frac{\partial f_a}{\partial t} + \xi_a \cdot \nabla{f_a} = -\frac{1}{\tau} ( f_a - f_a^{(eq)} ) + F_a \tag{26-1}

fa(eq)=waω(ξa)f(eq)(ξa)(26-2)f_a^{(eq)} = \frac{w_a}{\omega (\xi_a)} f^{(eq)} (\xi_a) \tag{26-2}

Fa=waω(ξa)F(ξa)(26-3)F_a = \frac{w_a}{\omega (\xi_a)} \mathbf{F} (\xi_a) \tag{26-3}

此时只需将 f(eq)f^{(eq)}FF 进行类似的展开,即可完成对整个系统的离散。

附录A. 使用C-E展开方法将离散的LBGK方程还原为Navier-Stokes方程

LBE方程为(A-1.1)式。源项 F\boldsymbol{F}f(eq)f^{(\mathrm{eq})} 为由郭照立、郑楚光和施保昌提出的格式,即(A-1.2)式和(A-1.3)式。[注意,这里的 ff 是密度分布函数]

fi(x+ξiΔt,t+Δt)fi(x,t)=1τ[fi(x,t)fi(eq)(x,t)]+ΔtFi(x,t)(A-1.1)f_i (x+ξ_i Δt,t+Δt)-f_i (x,t)=-\frac{1}{τ} \left[f_i (x,t)-f_i^{(eq)} (x,t) \right]+Δt⋅F_i (x,t) \tag{A-1.1}

Fi=(1Δt2τ)wi[ξiucs2+(ξiu)ξics4]F(A-1.2)F_i=(1 - \frac{Δt}{2τ}) w_i \left[ \frac{\mathbf{\xi_i-u} }{c_s^2 } + \frac{\mathbf{(\xi_i \cdot u) \xi_i} }{c_s^4} \right] ⋅ F \tag{A-1.2}

fi(eq)(x,t)=ρwi[1+ξiucs2+uu:(ξiξics2I)2cs4](A-1.3)f_i^{(eq)} (x,t) = \rho w_i \left[ 1 + \frac{\mathbf{ξ_i⋅u}}{c_s^2} +\frac{ \mathbf{uu} : ( \mathbf{ξ_i ξ_i} - c_s^2 \mathbf{I}) }{2c_s^4} \right] \tag{A-1.3}

并且对于 fif_ifi(eq)f_i^{(eq)} 有:

{ρ=ifiρu=iξifi+Δt2F\begin{cases} \rho = \sum_i {f_i} \\ \rho \mathbf{u} = \sum_i {\mathbf{\xi_i} f_i} + \frac{\Delta t}{2} F \end{cases}

{ρ=ifi(eq)ρu=iξifi(eq)ρ(uu+cs2I)=iξiξifi(eq)\begin{cases} \rho = \sum_i {f_i^{(eq)}} \\ \rho \mathbf{u} = \sum_i {\mathbf{\xi_i} f_i^{(eq)}} \\ \rho (\mathbf{uu} + c_s^2 \mathbf{I}) = \sum_i {\mathbf{\xi_i \xi_i} f_i^{(eq)}} \end{cases}

将(A-1.1)式在 (x,t)(\mathbf{x} , t)处展开,取至 2 阶项得:

Δtfit+ξiΔtfix+12[Δt22fit2+2Δt(ξiΔt)2fitx+Δt2(ξiξi):2fix2]+1τ[fifi(eq)]ΔtFi=0(A-2)\begin{aligned} \Delta t \frac{\partial f_i}{\partial t} + ξ_i \Delta t ⋅ \frac{\partial f_i}{∂x} + &\\ \frac{1}{2} \left[\Delta t^2 \frac{\partial^2 f_i}{\partial t^2} + 2\Delta t(\xi_i \Delta t) ⋅ \frac{\partial^2 f_i}{\partial t∂x} + \Delta t^2 (\xi_i \xi_i): \frac{\partial^2 f_i}{∂x^2} \right] + &\\ \frac{1}{\tau } \left[f_i - f_i^{(eq)} \right] - \Delta t F_i &= 0 \tag{A-2} \end{aligned}

代入多尺度展开等式,并舍去 ε3\varepsilon^3 项和更高阶的项,可得:

(εfi(0)t1+ε2fi(1)t1+ε2fi(0)t2)+ξi(εfi(0)x1+ε2fi(0)x1)+Δt2[ε22fi(0)t12+2ξiε22fi(0)t1x1+(ξiξi):(ε22fi(0)x2)]+1τΔt[fi(0)+εfi(1)+ε2fi(2)fi(eq)]εFi(1)+O(ε3)=0(A-3)\begin{aligned} (\varepsilon \frac{\partial f_i^{(0)}}{\partial t_1} + \varepsilon^2 \frac{\partial f_i^{(1)}}{\partial t_1} + \varepsilon^2 \frac{\partial f_i^{(0)}}{\partial t_2} ) + \xi_i ⋅ (\varepsilon \frac{\partial f_i^{(0)}}{\partial x_1} + \varepsilon^2 \frac{\partial f_i^{(0)}}{\partial x_1}) + &\\ \frac{Δt}{2} \left[ \varepsilon^2 \frac{\partial^2 f_i^{(0)}}{\partial t_1^2} + 2 \xi_i \varepsilon^2 ⋅ \frac{\partial^2 f_i^{(0)}}{\partial t_1 \partial x_1} + (\xi_i \xi_i) : (\varepsilon^2 \frac{\partial^2 f_i^{(0)}}{\partial x^2} ) \right] &\\ + \frac{1}{τΔt} \left[ f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)} -f_i^{(eq)} \right] - \varepsilon F_i^{(1)} + O(\varepsilon^3 ) &=0 \tag{A-3} \end{aligned}

从(A-3)式中提取O(ε0)O({\varepsilon}^0)项、O(ε1)O({\varepsilon}^1)项和O(ε2)O({\varepsilon}^2)项,并令 Dk=tk+ξkxkD_k = \frac{\partial}{\partial t_k} + \xi_k \cdot \frac{\partial}{\partial \mathbf{x}_k} 。(A-3)式对任意 ε\varepsilon 都成立,则关于各阶的系数均为0,即:

O(ε0)=1τΔt(fi(0)fi(eq))=0(A-4.1)O({\varepsilon}^0) = \frac{1}{\tau \Delta t} (f_i^{(0)} - f_i^{(eq)}) = 0 \tag{A-4.1}

O(ε1)=D1fi(0)+fi(1)τΔtFi(1)=0(A-4.2)O(\varepsilon^1) = D_1 f_i^{(0)} + \frac{ f_i^{(1)} }{\tau \Delta t} - F_i^{(1)} = 0 \tag{A-4.2}

O(ε2)=fi(0)t2(112τ)D1fi(1)+fi(2)τΔt+Δt2D1Fi(1)=0(A-4.3)O(\varepsilon^2) = \frac{\partial f_i^{(0)}}{\partial t_2} - (1 - \frac{1}{2\tau}) D_1 f_i^{(1)} + \frac{ f_i^{(2)} }{\tau \Delta t} + \frac{\Delta t}{2} D_1 F_i^{(1)} = 0 \tag{A-4.3}

根据多尺度展开的变量关系,可得:
ifi(1)=ifi(2)=0\sum_i {f_i^{(1)}} = \sum_i {f_i^{(2)}} = 0iξifi(2)=0\sum_i { \xi_i f_i^{(2)}} = 0iξifi(1)=Δt2Fi(1)\sum_i { \xi_i f_i^{(1)}} = - \frac{\Delta t}{2} F_i^{(1)}

{iFi=0iξiFi=(112τ)FiξiξiFi=(112τ)2uF\left\{\begin{matrix} \sum_{i} {F_i} = 0 \\ \sum_{i} {\xi_i F_i} = (1 - \frac{1}{2 \tau}) F \\ \sum_{i} {\xi_i \xi_i F_i} = (1 - \frac{1}{2 \tau}) 2 u F \end{matrix}\right.

{iFi(1)=0iξiFi(1)=(112τ)F(1)iξiξiFi(1)=(112τ)2uF(1)\left\{\begin{matrix} \sum_{i} {F_i^{(1)}} = 0 \\ \sum_{i} {\xi_i F_i^{(1)}} = (1 - \frac{1}{2 \tau}) F^{(1)} \\ \sum_{i} {\xi_i \xi_i F_i^{(1)}} = (1 - \frac{1}{2 \tau}) 2 u F^{(1)} \end{matrix}\right.

为方便下文分析,记: Π(k)=iξiξifi(k)\Pi^{(k)} = \sum_i { \xi_i \xi_i f_i^{(k)} }Γ(k)=iξiξiξifi(k)\Gamma^{(k)} = \sum_i { \xi_i \xi_i \xi_i f_i^{(k)} }
O(ε1)O({\varepsilon}^1) 项求1阶和2阶速度矩,得:

ρt1+(ρu)x1=0(A-5.1)\frac{\partial \rho}{\partial t_1} + \frac{\partial (\rho u)}{\partial x_1} = 0 \tag{A-5.1}

(ρu)t1+Π(0)x1=0(A-5.2)\frac{\partial (\rho u)}{\partial t_1} + \frac{\partial \Pi^{(0)}}{\partial x_1} = 0 \tag{A-5.2}

Π(0)=ρ(uu+cs2I)(A-5.3)\Pi^{(0)} = \rho (\it{\bf{uu}} \mathrm{+} c_s^2 \bf{I}) \tag{A-5.3}

O(ε2)O({\varepsilon}^2) 项求1阶和2阶速度矩,得:

ρt2=0(A-6.1)\frac{\partial {\rho}}{\partial t_{2}} = 0 \tag{A-6.1}

(ρu)t2+(112τ)Π(1)x1+Δt(112τ)F(1)x1=0(A-6.2)\frac{\partial (\rho \boldsymbol{u})}{\partial t_{2}} + (1 - \frac{1}{2 \tau}) \frac{\partial \Pi^{(1)}}{\partial \boldsymbol{x}_{1}} + \Delta t (1 - \frac{1}{2 \tau}) \frac{\partial \boldsymbol{F}^{(1)}}{\partial \boldsymbol{x}_{1}} = 0 \tag{A-6.2}

由(A-4.2)式,得: fi(1)=τΔt{D1fi(0)Fi(1)}f_i^{(1)} = \tau \Delta t \{ D_1 f_i^{(0)} - F_i^{(1)}\} ,所以

Π(1)=iξiξifi(1)=τΔt[(ρuu)t1+cs2(ρI)t1+Γ(0)x1(112τ)2uF(1)](A-7)\Pi^{(1)} = \sum_{i} \boldsymbol{\xi}_{i} \boldsymbol{\xi}_{i} f_{i}^{(1)} = -\tau \Delta t [\frac{\partial (\rho \boldsymbol{u} \boldsymbol{u})}{\partial t_{1}} + c_s^2 \frac{\partial (\rho \boldsymbol{I})}{\partial t_{1}} + \frac{\partial \Gamma^{(0)}}{\partial \boldsymbol{x}_{1}} - (1 - \frac{1}{2 \tau}) 2\boldsymbol{uF^{(1)}} ] \tag{A-7}

uF(1)\boldsymbol{uF^{(1)}}是Grad的论文中的记号, (uF(1))lm=12(ulFm(1)+umFl(1))(\bf{\it{uF}}^{(1)} )_{\it{lm}} = \mathrm{\frac{1}{2}} (\it{ u_l F_m^{(1)} + u_m F_l^{(1)} } \mathrm{)} 。下文为方便分析,取 Πlm(1)\Pi_{lm}^{(1)} 分析。

对第1项,有: (ρulum)t1=ρulumt1+um(ρul)t1\frac{\partial (\rho u_l u_m)}{\partial t_1} = \rho u_l \frac{\partial u_m}{\partial t_1} + u_m \frac{\partial (\rho u_l)}{\partial t_1} 。由(A-5.2)式,可得[4]

ρumt1+umρt1+x1n(ρumun+ρcs2δmn)=Fm(1)(A-8.1)\rho \frac{\partial u_m}{\partial t_1} + u_m \frac{\partial \rho}{\partial t_1} + \frac{\partial}{\partial x_{1n}} \left( \rho u_m u_n + \rho c_s^2 \delta_{mn} \right) = F_m^{(1)} \tag{A-8.1}

(ρul)t1+x1n(ρulun+ρcs2δln)=Fl(1)(A-8.2)\frac{\partial (\rho u_l)}{\partial t_1}+ \frac{\partial}{\partial x_{1n}} \left( \rho u_l u_n + \rho c_s^2 \delta_{ln} \right) = F_l^{(1)} \tag{A-8.2}

并且

(ρulumun)x1n=ulum(ρun)x1n+ul(ρumun)x1n+um(ρulun)x1n(A-8.3)\frac{\partial (\rho u_l u_m u_n)}{\partial x_{1n}} = - u_l u_m \frac{\partial (\rho u_n)}{\partial x_{1n}} + u_l \frac{\partial (\rho u_m u_n)}{\partial x_{1n}} + u_m \frac{\partial (\rho u_l u_n)}{\partial x_{1n}} \tag{A-8.3}

所以

(ρulum)x1n=(ulFm(1)+umFl(1))(ρulumun)x1ncs2(ulρx1m+umρx1l)(A-8.4)\frac{\partial (\rho u_l u_m)}{\partial x_{1n}} = (u_l F_m^{(1)} + u_m F_l^{(1)}) - \frac{\partial (\rho u_l u_m u_n)}{\partial x_{1n}} - c_s^2 \left( u_l \frac{\partial \rho}{\partial x_{1m}} + u_m \frac{\partial \rho}{\partial x_{1l}} \right) \tag{A-8.4}

所以第2项的分量为: cs2δlmρt1c_s^2 \delta_{lm} \frac{\partial \rho}{\partial t_1}

对第3项的分量,由于

Γlmn(0)=ρcs2(ulδmn+unδlm+umδln)\Gamma_{lmn}^{(0)} = \rho c_s^2 ( u_l \delta_{mn} + u_n \delta_{lm} + u_m \delta_{ln} )

,所以:

Γlmn(0)x1n=cs2(δmnρulx1n+δlmρunx1n+δlnρumx1n)(A-8.5)\frac{\partial \Gamma_{lmn}^{(0)}}{\partial x_{1n}} = c_s^2 ( \delta_{mn} \frac{\partial \rho u_l}{\partial x_{1n}} + \delta_{lm} \frac{\partial \rho u_n}{\partial x_{1n}} + \delta_{ln} \frac{\partial \rho u_m}{\partial x_{1n}} ) \tag{A-8.5}

所以

Πlm(1)=τΔt[(ρulumun)x1n+12τ(ulFm(1)+umFl(1))+ρcs2(ulx1m+umx1l)](A-8.6)\Pi_{lm}^{(1)} = - \tau \Delta t \left[ - \frac{\partial (\rho u_l u_m u_n)}{\partial x_{1n}} + \frac{1}{2 \tau} (u_l F_m^{(1)} + u_m F_l^{(1)}) + \rho c_s^2 ( \frac{\partial u_l}{\partial x_{1m}} + \frac{\partial u_m}{\partial x_{1l}} ) \right] \tag{A-8.6}

略去Πlm(1)\Pi_{lm}^{(1)} 中的速度的3阶项(即 x1n{ρulumun}\frac{\partial}{\partial x_{1n}} \{ \rho u_l u_m u_n \} ),有

Π(1)=ΔtuF(1)τρcs2Δt[ux1+(ux1)T](A-9)\Pi^{(1)} = - \Delta t \cdot \boldsymbol{u F}^{(1)} - \tau \rho c_s^2 \Delta t \left[ \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} + \left( \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} \right)^{\mathrm{T}} \right] \tag{A-9}

回代(A-6.2)式,得等价形式为

(ρu)t2+cs2Δt(τ12)x1{ρ[ux1+(ux1)T]}=0(A-10)\frac{\partial (\rho \boldsymbol{u})}{\partial t_2} + c_s^2 \Delta t (\tau - \frac{1}{2}) \frac{\partial}{\partial \boldsymbol{x_1}} \left\{ \rho \left[ \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} + \left( \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} \right)^{\mathrm{T}} \right] \right\} = \boldsymbol{0} \tag{A-10}

将(A-10)式和(A-5.2)式联立,得

(ρu)t+(ρuu)x=(ρcs2)x+cs2Δt(τ12)x{ρ[ux+(ux)T]}+F(A-11)\frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u u)}{\partial x} = - \frac{\partial (\rho c_s^2)}{\partial x} + c_s^2 \Delta t (\tau - \frac{1}{2}) \frac{\partial}{\partial x}\{\rho [\frac{\partial u}{\partial x} + (\frac{\partial u}{\partial x})^T ]\} + F \tag{A-11}

此时令 p=cs2ρp = c_s^2 \rhoν=cs2Δt(τ12)\nu = c_s^2 \Delta t (\tau - \frac{1}{2}) ,即可还原回Navier-Stokes方程

附录B.Gauss-Hermite积分

对于任意一维函数 f(ξ)f(\xi) ,高斯积分通过寻找积分点 ξk\xi_k 和其对应权重 wkw_k 实现对 abω(ξ)f(ξ)dξ\int_{a}^{b} \omega(\xi) f(\xi) d \xi 的近似,即

abω(ξ)f(ξ)dξk=1nwkξk(B-1)\int_{a}^{b} \omega(\xi) f(\xi) d \xi \approx \sum_{k=1}^{n} w_k \xi_k \tag{B-1}

其中 ω(ξ)\omega (\xi) 是任意的权重函数。n 点高斯求积的积分点正是对应的第n个正交多项式的根,此时 ξk\xi_k 对应的权重 wkw_k

wk=<Pn1,Pn1>Pn1(ξk)Pn(ξk)(B-2)w_k = \frac{ <P_{n-1}, P_{n-1}> }{ P_{n-1}(\xi_k) P_{n}^{'}(\xi_k) } \tag{B-2}

其中 Pn=dPndxP_{n}^{'} = \frac{d P_n}{d x} 。(B-2)式的精度为 2n12n-1
所以在一维Gauss-Hermite积分中,取 Pn=H(n)P_n = \mathcal{H}^{(n)}H(n)\mathcal{H}^{(n)} 的权函数和表达式见(12)式和(13)式(下文同理)。 n点高斯求积的积分点即为 H(n)\mathcal{H}^{(n)} 的根。因为 dH(n)dξ=ξH(n)H(n+1)=nH(n1)\frac{d \mathcal{H}^{(n)}}{d \xi} = \xi \mathcal{H}^{(n)} - \mathcal{H}^{(n+1)} = n \mathcal{H}^{(n-1)} ,所以有

wk=n!/[nH(n1)(ξk)]2(B-3)w_k = n! / [n \mathcal{H}^{(n-1)} (\xi_k)]^2 \tag{B-3}

对于更高维的情况,可做类似构造[3:1]。分析如下积分

1(2π)D/2exp(ξξ2)p(ξ)dξ(B-4)\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi \tag{B-4}

其中 ξ\xi 是长度为D的向量。 p(ξ)p(\xi) 是D维n次多项式,可写为(B-5)式

p(ξ)=n1+n2+...+nDnan1n2...nDj=1Dξjnj(B-5)p(\xi) = \sum\limits_{n_1 + n_2 + ... +n_D \le n} {a_{n_1 n_2 ... n_D}} \prod\limits_{j=1}\limits^{D}{\xi_{j}^{n_j}} \tag{B-5}

其中 an1n2...nDa_{n_1 n_2 ... n_D}为常数。

考虑(B-5)式中的任意一项 j=1Dξjnj\prod\limits_{j=1}\limits^{D}{\xi_{j}^{n_j}} (由于an1n2...nDa_{n_1 n_2 ... n_D} 为常数所以可暂不处理,在最后一步被归并到常数项中),回代到式(B-4)中,得

1(2π)D/2exp(ξξ2)p(ξ)dξ=j=1D(12πexp(ξj22)ξjnjdξj)=j=1D(k=1nwkξk)(B-6.1)\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi = \prod\limits_{j=1}\limits^{D} \left( \frac{1}{\sqrt{2 \pi}} \int \mathrm{exp}(- \frac{\xi_{j}^2}{2}) \xi_j^{n_j} \mathrm{d}\xi_j \right) = \prod\limits_{j=1}\limits^{D} ( \sum\limits_{k=1}\limits^{n} w_k \xi_k ) \tag{B-6.1}

其中 wkw_kξk\xi_k是一维n次多项式高斯积分的权重和积分点。积分结果被转化为多个一维Gauss-Hermite积分的连乘。注意到(B-6.1)式等价于

1(2π)D/2exp(ξξ2)p(ξ)dξ=k1=1n...kD=1n[(wk1wk2...wkD)(ξk1n1ξk2n2...ξkDnD)](B-6.2)\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi = \sum\limits_{k_1=1}\limits^{n} ... \sum\limits_{k_D=1}\limits^{n} [(w_{k_1} w_{k_2} ... w_{k_D}) (\xi_{k_1}^{n_1} \xi_{k_2}^{n_2} ... \xi_{k_D}^{n_D})] \tag{B-6.2}

如果定义 ξk1...kD=(ξk1,ξk2,...,ξkD)\xi_{k_1 ... k_D} = (\xi_{k_1}, \xi_{k_2}, ... , \xi_{k_D})wk1...kD=wk1wk2...wkDw_{k_1 ... k_D} = w_{k_1} w_{k_2} ... w_{k_D} ,那么(B-4)式等价为

1(2π)D/2exp(ξξ2)p(ξ)dξ=wk1...kDp(ξk1...kD)(B-7)\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi = \sum w_{k_1 ... k_D} p(\xi_{k_1 ... k_D}) \tag{B-7}

此时则将(B-4)式转换为与一维Gauss-Hermite积分类似的形式,即(B-7)式。

参考


  1. GRAD H. Note on N-Dimensional Hermite Polynomials[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 325–330. http://doi.org/10.1002/cpa.3160020402 ↩︎

  2. GRAD H. On the Kinetic Theory of Rarefied Gases[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 331–407. http://doi.org/10.1002/cpa.3160020403 ↩︎

  3. SHAN X, YUAN X-F, CHEN H. Kinetic Theory Representation of Hydrodynamics: A Way beyond the Navier–Stokes Equation[J]. Journal of Fluid Mechanics, 2006, 550(1): 413. http://doi.org/10.1017/S0022112005008153 ↩︎ ↩︎

  4. CAO W. Investigation of the Applicability of the Lattice Boltzmann Method to Free-Surface Hydrodynamic Problems in Marine Engineering[D/OL]. Laboratoire de recherche en Hydrodynamique, Énergétique et Environnement Atmosphérique (LHEEA): École centrale de Nantes, 2019. https://tel.archives-ouvertes.fr/tel-02383174 ↩︎

0%