[注] 本文同步发表在 知乎 .
宇宙安全声明： 1.这份文档来自我的一份课程作业, 由于笔者水平有限, 所以在内容上写得会比较粗糙, 或者出现一些错误.如有错误, 敬请谅解并指正. 2. 后续可能随缘更新？【如果又有空去读什么其他文章的话？】
Boltzmann方程源自气体动理论, 它通过建模粒子系统的概率分布函数, 描述其集体行为.速度分布函数 f ( x , ξ , t ) f\left( \mathbf{x},\mathbf{\xi},t \right) f ( x , ξ , t ) t t t x \mathbf{x} x d V = [ x , x + d x ] dV = [\mathbf{x}, \mathbf{x}+ d\mathbf{x}] d V = [ x , x + d x ] [ ξ , ξ + d ξ ] [\mathbf{\xi},\mathbf{\xi} + d\mathbf{\xi}] [ ξ , ξ + d ξ ] f ( x , ξ , t ) d ξ d V f\left( \mathbf{x},\mathbf{\xi},t \right)d\mathbf{\xi}dV f ( x , ξ , t ) d ξ d V f ( x , ξ , t ) f\left( \mathbf{x},\mathbf{\xi},t \right) f ( x , ξ , t )
D f D t = ∂ f ∂ t + ξ ⋅ ∂ f ∂ x + a ⋅ ∂ f ∂ ξ = Ω (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}
D t D f = ∂ t ∂ f + ξ ⋅ ∂ x ∂ f + a ⋅ ∂ ξ ∂ f = Ω ( 1 . 1 )
其中 ξ = ∂ x ∂ t \mathbf{\xi =}\frac{\partial\mathbf{x}}{\partial t} ξ = ∂ t ∂ x a = ∂ ξ ∂ t = F m \mathbf{a =}\frac{\partial\mathbf{\xi}}{\partial t}\mathbf{=}\frac{\mathbf{F}}{m} a = ∂ t ∂ ξ = m F F \mathbf{F} F m m m Ω \Omega Ω f f f
n = ∫ f d ξ , ρ = m ⋅ n , ρ u = m ∫ f ξ d ξ , ρ E = m ∫ ξ 2 2 f d ξ \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}
n = ∫ f d ξ , ρ u = m ∫ f ξ d ξ , ρ = m ⋅ n , ρ E = m ∫ 2 ξ 2 f d ξ
n n n ρ \rho ρ u \mathbf{u} u E E E
对分子的碰撞进行建模并非易事.为简化式(1.1)的计算, Bhatnagar等提出了BGK碰撞项.这一思路主张将复杂的碰撞视作当前状态 f f f f ( e q ) f^{(eq)} f ( e q )
Ω = − 1 τ 0 ( f − f ( e q ) ) \Omega = - \frac{1}{\tau_{0}}\left( f - f^{(eq)} \right)
Ω = − τ 0 1 ( f − f ( e q ) )
其中 τ 0 \tau_{0} τ 0 D D D
f ( e q ) = n ( 2 π R g T ) D 2 exp [ − ( ξ − u ) 2 2 R g T ] . (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}
f ( e q ) = ( 2 π R g T ) 2 D n exp [ − 2 R g T ( ξ − u ) 2 ] . ( 1 . 2 )
$R_g $ 为气体常数, T T T
在格子Boltzmann方法(Lattice Boltzmann Method, LBM)中, 使用的分布函数是密度分布函数 m f mf m f f f f
f ( e q ) = ρ ( 2 π R g T ) D 2 exp [ − ( ξ − u ) 2 2 R g T ] . (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}
f ( e q ) = ( 2 π R g T ) 2 D ρ exp [ − 2 R g T ( ξ − u ) 2 ] . ( 1 . 3 )
宏观量则表示为：
ρ = ∫ f d ξ , j = ρ u = ∫ f ξ d ξ , ρ E = ∫ ξ 2 2 f d ξ . \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}.
ρ = ∫ f d ξ , j = ρ u = ∫ f ξ d ξ , ρ E = ∫ 2 ξ 2 f d ξ .
广义上说, LBM是离散Boltzmann方法(Discrete Boltzmann Method, DBM)的一种, 是式(1.1)的一种数值格式.本节介绍从式(1.1)到格子Boltzmann方程的离散思路, 推导中仅考虑外力加速度 a = F / m ≡ 0 \mathbf{a} =\mathbf{F}/ m \equiv 0 a = F / m ≡ 0
对式(1.1)沿着特征线 ξ = x / t \mathbf{\xi} = \mathbf{x}/t ξ = x / t Δ t \Delta t Δ t O ( Δ t 2 ) O({\Delta t}^{2}) O ( Δ t 2 )
f ( x + ξ Δ t , ξ , t + Δ t ) − f ( x , ξ , t ) = − 1 τ [ f ( x , ξ , t ) − f ( e q ) ( 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}
f ( x + ξ Δ t , ξ , t + Δ t ) − f ( x , ξ , t ) = − τ 1 [ f ( x , ξ , t ) − f ( e q ) ( x , ξ , t ) ] ( 1 . 4 )
其中 τ = τ 0 / Δ t \tau = \tau_{0}/\Delta t τ = τ 0 / Δ t ξ \mathbf{\xi} ξ
∫ f ( e q ) ( x , ξ , t ) ⋅ Ψ ( ξ ) d ξ = ∑ i w i Ψ ( ξ i ) f ( e q ) ( 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}
∫ f ( e q ) ( x , ξ , t ) ⋅ Ψ ( ξ ) d ξ = i ∑ w i Ψ ( ξ i ) f ( e q ) ( x , ξ i , t ) ( 1 . 5 )
这里 Ψ ( ξ ) \Psi\left( \mathbf{\xi} \right) Ψ ( ξ ) ξ \mathbf{\xi} ξ ξ i \mathbf{\xi}_i ξ i w i w_i w i
ρ = ∑ i f i = ∑ i f i ( e q ) , j = ρ u = ∑ i ξ i f i = ∑ i ξ i f i ( e q ) \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)}}
ρ = i ∑ f i = i ∑ f i ( e q ) , j = ρ u = i ∑ ξ i f i = i ∑ ξ i f i ( e q )
以及内能 e e e
ρ e = 1 2 ∑ i ( ξ i − u ) 2 f i = 1 2 ∑ i ( ξ i − u ) 2 f i ( e q ) \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)}}
ρ e = 2 1 i ∑ ( ξ i − u ) 2 f i = 2 1 i ∑ ( ξ i − u ) 2 f i ( e q )
其中 f i = W i f ( x , ξ i , t ) f_i = W_if\left( \mathbf{x},\mathbf{\xi}_i\mathbf{,}t \right) f i = W i f ( x , ξ i , t ) f i ( e q ) = W i f ( e q ) ( x , ξ i , t ) f_i^{(eq)} = W_if^{(eq)}\left( \mathbf{x},\mathbf{\xi}_i\mathbf{,}t \right) f i ( e q ) = W i f ( e q ) ( x , ξ i , t ) { c i } ( i = 0 , 1 , … , b − 1 ) \{\mathbf{c}_i\}(i = 0,1,\ldots,b - 1) { c i } ( i = 0 , 1 , … , b − 1 ) n 维空间下 b 速度的离散速度模型简称 Dn Qb 模型.
在低马赫数下, 平衡态可对 u \mathbf{u} u
f ( e q ) ( x , ξ , t ) = ρ ( 2 π R g T ) D 2 ⋅ exp ( − ξ 2 2 R g T ) ⋅ exp ( 2 ( ξ ⋅ u ) − u 2 2 R g T ) = ρ exp ( − ξ 2 2 R g T ) ( 2 π R g T ) D 2 ⋅ [ 1 + ξ ⋅ u R g T + ( ξ ⋅ u ) 2 2 ( R g T ) 2 − u 2 2 R g T ] + O ( u 3 ) (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 ( e q ) ( x , ξ , t ) = ( 2 π R g T ) 2 D ρ ⋅ exp ( − 2 R g T ξ 2 ) ⋅ exp ( 2 R g T 2 ( ξ ⋅ u ) − u 2 ) = ( 2 π R g T ) 2 D ρ exp ( − 2 R g T ξ 2 ) ⋅ [ 1 + R g T ξ ⋅ u + 2 ( R g T ) 2 ( ξ ⋅ u ) 2 − 2 R g T u 2 ] + O ( u 3 ) ( 1 . 6 )
因此, 当使用 截断的平衡态 计算 ∫ f ( e q ) ⋅ Ψ ( ξ ) d ξ \int {f^{(eq)} \cdot \Psi\left( \mathbf{\xi} \right)}d\mathbf{\xi} ∫ f ( e q ) ⋅ Ψ ( ξ ) d ξ ∫ exp ( − x 2 ) ⋅ Ψ ( x ) d x \int {\exp\left( - x^{2} \right) \cdot \Psi(x)}dx ∫ exp ( − x 2 ) ⋅ Ψ ( x ) d x
I m = ∫ − ∞ + ∞ exp ( − ζ 2 ) ⋅ ζ m d ζ I_{m} = \int_{- \infty}^{+ \infty}{\exp\left( - \zeta^{2} \right) \cdot \zeta^{m}} d\zeta
I m = ∫ − ∞ + ∞ exp ( − ζ 2 ) ⋅ ζ m d ζ
其中 ζ = ξ / 2 R g T \zeta = \xi/\sqrt{2 R_g T} ζ = ξ / 2 R g T I m I_{m} I m exp ( − ζ 2 ) \exp\left(-\zeta^{2}\right) exp ( − ζ 2 ) m m m
这里以D2Q9离散速度集为例介绍权重的计算, 其离散速度集为
c α = { ( 0 , 0 ) , α = 0 c ( cos π ( α − 1 ) 2 , sin π ( α − 1 ) 2 ) , α = 1 … 4 2 c ( cos π ( 2 α − 1 ) 4 , sin π ( 2 α − 1 ) 4 ) , α = 5 … 8 (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 α = ⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ ( 0 , 0 ) c ( cos 2 π ( α − 1 ) , sin 2 π ( α − 1 ) ) 2 c ( cos 4 π ( 2 α − 1 ) , sin 4 π ( 2 α − 1 ) ) , α = 0 , α = 1 … 4 , α = 5 … 8 ( 1 . 7 )
其中 c = Δ x / Δ t c = \Delta x/\Delta t c = Δ x / Δ t Ψ \Psi Ψ Ψ m , n ( ξ ) = ξ x m ξ y n \Psi_{m,n}\left( \mathbf{\xi} \right) = \xi_{x}^{m}\xi_{y}^{n} Ψ m , n ( ξ ) = ξ x m ξ y n
I = ∫ f ( e q ) ⋅ Ψ m , n ( ξ ) d ξ = ρ π ( 2 R g T ) m + n 2 [ ( 1 − u 2 2 R g T ) I m I n + 2 2 R g T ( u x I m + 1 I n + u y I m I n + 1 ) + 1 R g T ( u x 2 I m + 2 I n + 2 u x u y I m + 1 I n + 1 + u y 2 I m I n + 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}
I = ∫ f ( e q ) ⋅ Ψ m , n ( ξ ) d ξ = π ρ ( 2 R g T ) 2 m + n [ ( 1 − 2 R g T u 2 ) I m I n + 2 R g T 2 ( u x I m + 1 I n + u y I m I n + 1 ) + R g T 1 ( u x 2 I m + 2 I n + 2 u x u y I m + 1 I n + 1 + u y 2 I m I n + 2 ) ] ( 1 . 8 )
对D2Q9模型, I m I_{m} I m
I m = ∫ − ∞ + ∞ exp ( − ζ 2 ) ⋅ ζ m d ζ = ∑ α = 1 3 ω α ζ α m I_{m} = \int_{- \infty}^{+ \infty}{\exp\left( - \zeta^{2} \right) \cdot \zeta^{m}}d\zeta = \sum_{\alpha = 1}^{3}{\omega_{\alpha}\zeta_{\alpha}^{m}}
I m = ∫ − ∞ + ∞ exp ( − ζ 2 ) ⋅ ζ m d ζ = α = 1 ∑ 3 ω α ζ α m
[NOTE] ζ i \zeta_i ζ i H 3 ( x ) = 8 x 3 − 12 x H_{3}(x) = 8x^{3} - 12x H 3 ( x ) = 8 x 3 − 1 2 x ω i = 2 3 − 1 ⋅ 3 ! H 2 ( ζ i ) 2 ⋅ 3 2 ⋅ π \omega_i = \frac{2^{3 - 1} \cdot 3!}{H_{2}\left( \zeta_i \right)^{2} \cdot 3^{2}} \cdot \sqrt{\pi} ω i = H 2 ( ζ i ) 2 ⋅ 3 2 2 3 − 1 ⋅ 3 ! ⋅ π ζ 1 = − 3 / 2 , w 1 = π / 6 \zeta_{1} = - \sqrt{3/2},w_{1} = \sqrt{\pi}/6 ζ 1 = − 3 / 2 , w 1 = π / 6 ζ 2 = 0 , w 2 = 2 π / 3 \zeta_{2} = 0,w_{2} = 2\sqrt{\pi}/3 ζ 2 = 0 , w 2 = 2 π / 3 ζ 3 = 3 / 2 , ω 3 = π / 6 \zeta_{3} = \sqrt{3/2},\omega_{3} = \sqrt{\pi}/6 ζ 3 = 3 / 2 , ω 3 = π / 6
则式(1.8)可被数值离散为 ：
I = ρ ∑ i , j = 1 3 Ψ ( ξ i , j ) ⋅ ω i ω j π ⋅ [ 1 + ξ i , j ⋅ u R g T + ( ξ i , j ⋅ u ) 2 2 ( R g T ) 2 − u 2 2 R g T ] (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 = ρ i , j = 1 ∑ 3 Ψ ( ξ i , j ) ⋅ π ω i ω j ⋅ [ 1 + R g T ξ i , j ⋅ u + 2 ( R g T ) 2 ( ξ i , j ⋅ u ) 2 − 2 R g T u 2 ] ( 1 . 9 )
此处 ξ i , j = ( ξ i , ξ j ) = 2 R g T ( ζ 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) ξ i , j = ( ξ i , ξ j ) = 2 R g T ( ζ i , ζ j ) R g T = c s 2 = c 2 / 3 R_g T = c_{s}^{2} = c^{2}/3 R g T = c s 2 = c 2 / 3 ξ i , j \mathbf{\xi}_{i,j} ξ i , j c α \mathbf{c}_{\alpha} c α c α \mathbf{c}_{\alpha} c α f α ( e q ) f_{\alpha}^{(eq)} f α ( e q )
f α ( e q ) = ρ w α ⋅ [ 1 + c α ⋅ u c s 2 + ( c α ⋅ u ) 2 2 c s 4 − u 2 2 c s 2 ] (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}
f α ( e q ) = ρ w α ⋅ [ 1 + c s 2 c α ⋅ u + 2 c s 4 ( c α ⋅ u ) 2 − 2 c s 2 u 2 ] ( 1 . 1 0 )
其中
w α = { 4 / 9 , i = j = 2 , α = 0 1 / 9 , ( i = 2 ) ⨁ ( j = 2 ) , α = 1 , … , 4 1 / 36 , ( i ≠ 2 ) ⋀ ( j ≠ 2 ) , α = 5 , … , 8 w_{\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}
w α = ⎩ ⎪ ⎪ ⎨ ⎪ ⎪ ⎧ 4 / 9 , 1 / 9 , 1 / 3 6 , i = j = 2 , ( i = 2 ) ⨁ ( j = 2 ) , ( i = 2 ) ⋀ ( j = 2 ) , α = 0 α = 1 , … , 4 α = 5 , … , 8
Qian等已给出常用的DnQb模型参数, 大量学者也给出了其他离散速度模型的推导方式 .
在前文中, R g T R_g T R g T
就目前来说, 热流体的LBM模拟可以分为多速LBE模型 、双分布函数LBE模型 、基于简化LBM的实现、热流格子Boltzmann通量求解器、多松弛模型 以及与其他数值方法结合的混合模型等多种形式 .
多速LBE模型(即Multi-Speed LBE)是一种完全基于LBM经典算法的改进.这类方法演化方程与常规LBM一致, 在保持算法不变的基础上通过使用比常规LBM更为复杂的离散速度模型还原高阶速度矩.它因只使用单一分布函数,物理描述上更接近事实而受到了广泛关注.并且使用更多的离散速度达到高阶精度的思路亦已经被广泛拓展到各种各样的LBM模拟中. 虽然Qian和Alexander等才是早期的提出者, 但本章首先以Chen等 的高阶MS-LBGK模型的二维形式为例引入多速LBM的概念, 随后介绍几种不同思路下的MS-LBM改进形式.
Chen等 提出的多速LBM的平衡态分布函数定义为
f σ i ( e q ) = A σ + B σ ⋅ ( c σ i ⋅ u ) + C σ ∥ u ∥ 2 + D σ ⋅ ( c σ i ⋅ u ) 2 + E σ ∥ u ∥ 2 ⋅ ( c σ i ⋅ u ) + F σ ⋅ ( c σ i ⋅ u ) 3 + G σ ∥ u ∥ 2 ⋅ ( c σ i ⋅ u ) 2 + H σ ∥ u ∥ 4 (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}
f σ i ( e q ) = A σ + B σ ⋅ ( c σ i ⋅ u ) + C σ ∥ u ∥ 2 + D σ ⋅ ( c σ i ⋅ u ) 2 + E σ ∥ u ∥ 2 ⋅ ( c σ i ⋅ u ) + F σ ⋅ ( c σ i ⋅ u ) 3 + G σ ∥ u ∥ 2 ⋅ ( c σ i ⋅ u ) 2 + H σ ∥ u ∥ 4 ( 2 . 1 )
其中 A σ A_{\sigma} A σ H σ H_{\sigma} H σ e e e ρ \rho ρ X σ = ( x σ 0 + x σ 1 e + x σ 2 e 2 ) ρ X_{\sigma} = \left( x_{\sigma 0} + x_{\sigma 1}e + x_{\sigma 2}e^{2} \right)\rho X σ = ( x σ 0 + x σ 1 e + x σ 2 e 2 ) ρ f σ i ( e q ) f_{\sigma i}^{(eq)} f σ i ( e q )
∑ σ i f σ i ( e q ) = ρ , ∑ σ i c σ i f σ i ( e q ) = ρ u , ∑ σ i c σ i c σ i f σ i ( e q ) = ρ u u + p I \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}
σ i ∑ f σ i ( e q ) = ρ , σ i ∑ c σ i f σ i ( e q ) = ρ u , σ i ∑ c σ i c σ i f σ i ( e q ) = ρ u u + p I
∑ σ i c σ i c σ i c σ i f σ i ( e q ) = ρ u u u + 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} ]
σ i ∑ c σ i c σ i c σ i f σ i ( e q ) = ρ u u u + p [ u δ ]
∑ σ i c σ i c σ i ∣ c σ i ∣ 2 f σ i ( e q ) = ρ ∣ u ∣ 2 u u + p ( D + 4 ) u u + p ∣ u ∣ 2 I + 2 ( D + 2 ) p e D \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}
σ i ∑ c σ i c σ i ∣ c σ i ∣ 2 f σ i ( e q ) = ρ ∣ u ∣ 2 u u + p ( D + 4 ) u u + p ∣ u ∣ 2 I + D 2 ( D + 2 ) p e
其中压强 p = 2 ρ e D p = \frac{2\rho e}{D} p = D 2 ρ e [ 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} [ u δ ] = u α δ β γ + u β δ α γ + u γ δ α β
为满足这五个守恒方程的约束, 张量 ∑ i c σ i c σ i . . . c σ i ⏞ n 个 \displaystyle \sum_i \overset{n\text{个}}{\overbrace{\mathbf{c}_{\sigma i}\mathbf{c}_{\sigma i}...\mathbf{c}_{\sigma i}}}\qquad i ∑ c σ i c σ i . . . c σ i n 个
c σ i = c p k i ′ = P e r m { k ( ± 1 , … , ± 1 ⏟ p , 0 , . . . , 0 ⏞ D − p ) } \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\}
c σ i = c p k i ′ = P e r m ⎩ ⎪ ⎨ ⎪ ⎧ k ⎝ ⎛ p ± 1 , … , ± 1 , 0 , . . . , 0 D − p ⎠ ⎞ ⎭ ⎪ ⎬ ⎪ ⎫
这里 P e r m \mathbf{Perm} P e r m ( ∗ ) (*) ( ∗ ) k k k σ = c σ i 2 = k 2 p \sigma = c_{\sigma i}^{2} = k^{2}p σ = c σ i 2 = k 2 p μ = 2 ρ e D ( τ − 1 2 ) \mu = \frac{2\rho e}{D}\left( \tau - \frac{1}{2} \right) μ = D 2 ρ e ( τ − 2 1 ) λ = − 2 μ D \lambda = - \frac{2\mu}{D} λ = − D 2 μ κ = ( D + 2 ) ρ e D ( τ − 1 2 ) \kappa = \frac{(D + 2)\rho e}{D}\left( \tau - \frac{1}{2} \right) κ = D ( D + 2 ) ρ e ( τ − 2 1 ) c p = ( D + 2 ) R g / 2 c_{p} = (D + 2)R_g /2 c p = ( D + 2 ) R g / 2 c v = D R g / 2 c_{v} = D R_g /2 c v = D R g / 2 R g = 1 R_g = 1 R g = 1 P r = μ c p / κ = 1 Pr = \mu c_{p}/\kappa = 1 P r = μ c p / κ = 1
目前, 学界已有数类方式在多速LBE模型中实现对Prandtl数的自由调节, 如：双松弛时间方法；引入熵格式 ；修改网格速度的定义 等.由于已有文献指出若不局限于常规晶格模型则可获得更高精度的离散格式, 并且大多数方法都涉及对Dn Qb模型的修改, 因此这里仅简要介绍前两种方式的实现思路.
基于双松弛时间的方法将碰撞分为由 f i f_i f i f − i f_{-i} f − i
f i ( x + c i Δ t , t + Δ t ) − f i ( 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}
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = Ω i + Ω i ( 2 . 2 )
其中
Ω i = − 1 τ 1 ( f i − f i ( e q ) ) , Ω i = − 1 τ 2 ( f − i − f − i ( e q ) ) \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 = − τ 1 1 ( f i − f i ( e q ) ) , Ω i = − τ 2 1 ( f − i − f − i ( e q ) )
下标 − i -i − i c − i = − c i \mathbf{c}_{- i} = - \mathbf{c}_i c − i = − c i τ v = τ 1 τ 2 / ( τ 1 + τ 2 ) \tau_{v} = \tau_{1}\tau_{2}/(\tau_{1} + \tau_{2}) τ v = τ 1 τ 2 / ( τ 1 + τ 2 ) τ k = τ 1 τ 2 / ( τ 2 − τ 1 ) \tau_{k} = \tau_{1}\tau_{2}/(\tau_{2} - \tau_{1}) τ k = τ 1 τ 2 / ( τ 2 − τ 1 )
D ρ D t = ρ ∇ ⋅ u ∂ ( ρ u ) ∂ t + ∇ ⋅ ( ρ u u ) = − ∇ p + ∇ ⋅ T v ∂ ( ρ e ) ∂ t + ∇ ⋅ ( ρ e u ) = p ∇ ⋅ u + ∇ ⋅ ( κ ∇ T ) + T k : ∇ 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}
D t D ρ ∂ t ∂ ( ρ u ) + ∇ ⋅ ( ρ u u ) ∂ t ∂ ( ρ e ) + ∇ ⋅ ( ρ e u ) = ρ ∇ ⋅ u = − ∇ p + ∇ ⋅ T v = p ∇ ⋅ u + ∇ ⋅ ( κ ∇ T ) + T k : ∇ u
其中
T v = 2 μ v S + λ v ( ∇ ⋅ u ) I , T k = 2 μ k S + λ 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}
T v = 2 μ v S + λ v ( ∇ ⋅ u ) I , T k = 2 μ k S + λ k ( ∇ ⋅ u ) I
S \mathbf{S} S
μ v = 2 D ρ e ( τ v − 1 2 ) , λ v = − 4 D 2 ρ e ( τ v − 1 2 ) μ k = 2 D ρ e ( τ k − 1 2 ) , λ k = − 4 D 2 ρ e ( τ k − 1 2 ) \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}
μ v = D 2 ρ e ( τ v − 2 1 ) , μ k = D 2 ρ e ( τ k − 2 1 ) , λ v = D 2 − 4 ρ e ( τ v − 2 1 ) λ k = D 2 − 4 ρ e ( τ k − 2 1 )
温度 T = D e / 2 T = De/2 T = D e / 2 κ = ρ e ( D + 2 ) ( τ k − 1 / 2 ) / D \kappa = \rho e(D + 2)\left( \tau_{k} - 1/2 \right)/D κ = ρ e ( D + 2 ) ( τ k − 1 / 2 ) / D
P r = μ v c p κ = 2 τ v − 1 2 τ k − 1 Pr = \frac{\mu_{v}c_{p}}{\kappa} = \frac{2\tau_{v} - 1}{2\tau_{k} - 1}
P r = κ μ v c p = 2 τ k − 1 2 τ v − 1
这种方法通过使用两个松弛时间 τ 1 , τ 2 \tau_{1},\tau_{2} τ 1 , τ 2 τ 1 , τ 2 > 1 2 \tau_{1},\tau_{2} > \frac{1}{2} τ 1 , τ 2 > 2 1 T k ≠ T v \mathbf{T}_{k} \neq \mathbf{T}_{v} T k = T v
熵格子Boltzmann方程(Entropic Lattice Boltzmann Equation, ELBE)是另一种基于传统LBE的修改.这里以单松弛的ELBE为例, 它将LBE的演化方程修改为
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − α β ( f i − f i ( e q ) ) (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}
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − α β ( f i − f i ( e q ) ) ( 2 . 3 )
其中 β = 1 / ( 1 + 2 ν / c s 2 ) \beta = 1/(1 + 2\nu/c_{s}^{2}) β = 1 / ( 1 + 2 ν / c s 2 ) ν \nu ν α \alpha α H ( f + α ( f − f ( e q ) ) ) = H ( f ) H\left( f + \alpha\left( f - f^{(eq)} \right) \right) = H(f) H ( f + α ( f − f ( e q ) ) ) = H ( f ) H H H
H = ∑ i f i ln ( f i w i ) (2.4) H = \sum_i^{}{f_i\ln\left( \frac{f_i}{w_i} \right)}
\tag{2.4}
H = i ∑ f i ln ( w i f i ) ( 2 . 4 )
α \alpha α n Qb模型在使用 f ( e q ) f^{(eq)} f ( e q ) n Qb模型进行修改, 或搭配多速Dn Qb模型一同使用.
Pransianakis等 在ELBE的总体框架上, 将碰撞项视为从 f i f_i f i f i ∗ f_i^{*} f i ∗ f i ( e q ) f_i^{(eq)} f i ( e q )
− 1 τ 1 ( f i − f i ∗ ) − 1 τ 2 ( f i ∗ − f i ( e q ) ) \frac{-1}{\tau_{1}}\left( f_i - f_i^{*} \right) - \frac{1}{\tau_{2}}\left( f_i^{*} - f_i^{(eq)} \right) τ 1 − 1 ( f i − f i ∗ ) − τ 2 1 ( f i ∗ − f i ( e q ) )
中间态 f i ∗ f_i^{*} f i ∗ f i ∗ = f i ( e q ) + δ f i ∗ f_i^{*} = f_i^{(eq)} + \delta f_i^{*} f i ∗ = f i ( e q ) + δ f i ∗ δ f i ∗ \delta f_i^{*} δ f i ∗ δ f i ∗ \delta f_i^{*} δ f i ∗ τ 1 = μ / ρ T 0 \tau_{1} = \mu/\rho T_0 τ 1 = μ / ρ T 0 τ 2 = 2 κ / ρ T 0 \tau_{2} = 2\kappa/\rho T_0 τ 2 = 2 κ / ρ T 0 T 0 T_0 T 0 Pr = 4 τ 1 / τ 2 \Pr = 4\tau_{1}/\tau_{2} Pr = 4 τ 1 / τ 2
[NOTE] c i \mathbf{c}_i c i W i W_i W i T T T W i = ( 1 − T 2 ) T 2 ( 1 − T ) c i 2 \displaystyle W_i = \frac{\left( 1 - T^{2} \right)T}{2(1 - T)}c_i^{2} W i = 2 ( 1 − T ) ( 1 − T 2 ) T c i 2
而平衡态分布写作： f i ( e q ) = ρ W i { 1 + c i ⋅ j ρ T + j j T 2 ( ρ T ) 2 : [ c i c i T − 4 T 2 + c i 2 ( 1 − 3 T ) 2 ( 1 − T ) 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\} f i ( e q ) = ρ W i { 1 + ρ T c i ⋅ j + 2 ( ρ T ) 2 j j T : [ c i c i T − 2 ( 1 − T ) 4 T 2 + c i 2 ( 1 − 3 T ) I ] }
Frapolli等提出了将多速LBE模型与熵LBE框架进行结合的思路.以D2Q25-ZOT模型为例, 其演化方程为：
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = ω 1 ( f i − f i ∗ ) + ω 2 ( f i ∗ − f i ( e q ) ) (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}
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = ω 1 ( f i − f i ∗ ) + ω 2 ( f i ∗ − f i ( e q ) ) ( 2 . 5 )
在不考虑ELBE的校正时, 动力粘度 μ \mu μ κ \kappa κ
μ = { ( 1 ω 1 − 1 2 ) ρ T , Pr ≤ 1 ( 1 ω 2 − 1 2 ) ρ T , Pr ≥ 1 , κ = { ( 1 ω 2 − 1 2 ) ρ T c p , Pr ≤ 1 ( 1 ω 1 − 1 2 ) ρ T c p , Pr ≥ 1 \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}
μ = ⎩ ⎪ ⎨ ⎪ ⎧ ( ω 1 1 − 2 1 ) ρ T , ( ω 2 1 − 2 1 ) ρ T , Pr ≤ 1 Pr ≥ 1 , κ = ⎩ ⎪ ⎨ ⎪ ⎧ ( ω 2 1 − 2 1 ) ρ T c p , ( ω 1 1 − 2 1 ) ρ T c p , Pr ≤ 1 Pr ≥ 1
Frapolli等认为, 上述计算中仅需要对涉及 μ \mu μ Pr < 1 \Pr < 1 Pr < 1 ω 1 = α β 1 , ω 2 = β 2 \omega_{1} = \alpha\beta_{1},\omega_{2} = \beta_{2} ω 1 = α β 1 , ω 2 = β 2 Pr > 1 \Pr > 1 Pr > 1 ω 1 = β 1 , ω 2 = α β 2 \omega_{1} = \beta_{1},\omega_{2} = \alpha\beta_{2} ω 1 = β 1 , ω 2 = α β 2 α \alpha α H ( f + α ( f − f ( e q ) ) ) = H ( f ) H\left( f + \alpha\left( f - f^{(eq)} \right) \right) = H(f) H ( f + α ( f − f ( e q ) ) ) = H ( f ) μ \mu μ κ \kappa κ
μ = { 1 2 ( 1 β 1 − 1 ) ρ T , Pr ≤ 1 1 2 ( 1 β 2 − 1 ) ρ T , Pr ≥ 1 , κ = { ( 1 β 2 − 1 2 ) ρ T , Pr ≤ 1 ( 1 β 1 − 1 2 ) ρ T , Pr ≥ 1 . \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.\ .
μ = ⎩ ⎪ ⎨ ⎪ ⎧ 2 1 ( β 1 1 − 1 ) ρ T , Pr ≤ 1 2 1 ( β 2 1 − 1 ) ρ T , Pr ≥ 1 , κ = ⎩ ⎪ ⎨ ⎪ ⎧ ( β 2 1 − 2 1 ) ρ T , Pr ≤ 1 ( β 1 1 − 2 1 ) ρ T , Pr ≥ 1 .
[NOTE] Pr = 1 \Pr = 1 Pr = 1
上述所介绍的基于ELBE的热流模拟能够成功还原Fourier-Navier-Stokes方程组, 实现低马赫数下各种热流的数值模拟, 但也没有实现热流和流场的耦合.此外, 由于 H H H H ( f + α ( f − f ( e q ) ) ) = H ( f ) H\left( f + \alpha\left( f - f^{(eq)} \right) \right) = H(f) H ( f + α ( f − f ( e q ) ) ) = H ( f )
双分布函数(double distribution function, DDF)的LBE模型的主要思路是将速度场和温度场分别用两套LBE进行计算.DDF-LBE下的热流分为两类：其一是将温度视作被动标量, 这要求温度变化对流动的影响较小；另一类考虑的温度对流场的影响(如可压缩热流).
Bartoloni等 指出：在压力做功和黏性热耗散项均可被忽略的场景中, 可将流场和温度场分别采用两套分布函数体系进行对流扩散.一般来说, 无源项的DDF-LBE演化方程组为：
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − 1 τ 1 ( f i − f i ( e q ) ) (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}
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − τ 1 1 ( f i − f i ( e q ) ) ( 3 . 1 )
g i ( x + c i Δ t , t + Δ t ) − g i ( x , t ) = − 1 τ 2 ( g i − g i ( e q ) ) (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}
g i ( x + c i Δ t , t + Δ t ) − g i ( x , t ) = − τ 2 1 ( g i − g i ( e q ) ) ( 3 . 2 )
其中 f i f_i f i g i g_i g i
ρ = ∑ i f i , j = ρ u = ∑ i c i f i , ρ T = ∑ i g i (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}
ρ = i ∑ f i , j = ρ u = i ∑ c i f i , ρ T = i ∑ g i ( 3 . 3 )
需要强调的是, f i f_i f i g i g_i g i c s 1 c_{s1} c s 1 c s 2 c_{s2} c s 2 τ 1 , τ 2 \tau_{1},\tau_{2} τ 1 , τ 2 ν \nu ν κ \kappa κ
τ 1 = ν c s 1 2 + 1 2 , τ 2 = κ c s 2 2 + 1 2 (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}
τ 1 = c s 1 2 ν + 2 1 , τ 2 = c s 2 2 κ + 2 1 ( 3 . 4 )
这种思路仅额外添加了一组分布函数及其平衡态的定义, 并且由于其易于实现, 目前是低速小Mach数场景下热流的一种主流LBM计算方法.
前面介绍的DDF-LBE都是在忽略黏性热耗散和压缩功的前提下, 单独另开一个分布函数 g i g_i g i T T T
既然可以为温度T T T ρ e \rho e ρ e ρ E = ρ e + ρ u 2 / 2 \rho E = \rho e + \rho u^{2}/2 ρ E = ρ e + ρ u 2 / 2 h = e + p / ρ h = e + p/\rho h = e + p / ρ ρ E \rho E ρ E
单位体积内的总能可以由分布函数的速度矩表示为
ρ E = ρ e + ρ u 2 2 = ρ c v T + ρ u 2 2 = ∫ 1 2 ξ 2 f d ξ (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}
ρ E = ρ e + 2 ρ u 2 = ρ c v T + 2 ρ u 2 = ∫ 2 1 ξ 2 f d ξ ( 3 . 5 )
因此可以将总能的分布函数定义为 h = 1 2 ξ 2 f h = \frac{1}{2}\xi^{2}f h = 2 1 ξ 2 f ρ E = ∫ h d ξ \rho E = \int h d\mathbf{\xi} ρ E = ∫ h d ξ
根据 f f f h h h
∂ h ∂ t + ξ ⋅ ∇ f + a ⋅ ∇ ξ h − f ξ ⋅ 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}
∂ t ∂ h + ξ ⋅ ∇ f + a ⋅ ∇ ξ h − f ξ ⋅ a = Ω h ( 3 . 6 )
其中 Ω h = 1 2 ξ 2 Ω f \Omega_{h} = \frac{1}{2}\xi^{2}\Omega_{f} Ω h = 2 1 ξ 2 Ω f Ω h \Omega_{h} Ω h τ f \tau_{f} τ f τ h \tau_{h} τ h Ω h \Omega_{h} Ω h
Ω h = − 1 τ h ( h − h ( e q ) ) + Z τ h f ( f − f ( e q ) ) (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}
Ω h = − τ h 1 ( h − h ( e q ) ) + τ h f Z ( f − f ( e q ) ) ( 3 . 7 )
其中 Z = ξ ⋅ u − u 2 2 Z = \mathbf{\xi \cdot u} - \frac{u^{2}}{2} Z = ξ ⋅ u − 2 u 2 1 τ h f = 1 τ h − 1 τ f \frac{1}{\tau_{hf}} = \frac{1}{\tau_{h}} - \frac{1}{\tau_{f}} τ h f 1 = τ h 1 − τ f 1
多原子分子气体的 f f f f ( e q ) f^{(eq)} f ( e q )
f ~ ( e q ) ( x , ξ , η , t ) = ρ ( 2 π R g T ) D + K 2 ⋅ exp [ − ( ξ − u ) 2 + η 2 2 R g T ] (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}
f ( e q ) ( x , ξ , η , t ) = ( 2 π R g T ) 2 D + K ρ ⋅ exp [ − 2 R g T ( ξ − u ) 2 + η 2 ] ( 3 . 8 )
其中 η \mathbf{\eta} η f ‾ = ∫ f d η \overline{f} = \int f d\mathbf{\eta} f = ∫ f d η h ‾ = ∫ ξ 2 + η 2 2 f d η \overline{h}\mathbf{=}\int \frac{\mathbf{\xi}^{2}\mathbf{+}\mathbf{\eta}^{2}}{2}f d\mathbf{\eta\ } h = ∫ 2 ξ 2 + η 2 f d η
∂ f ‾ ∂ t + ξ ⋅ ∇ f ‾ + a ⋅ ∇ ξ f ‾ = − 1 τ f ( f ‾ − f ‾ ( e q ) ) , (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}
∂ t ∂ f + ξ ⋅ ∇ f + a ⋅ ∇ ξ f = − τ f 1 ( f − f ( e q ) ) , ( 3 . 9 )
∂ h ‾ ∂ t + ξ ⋅ ∇ h ‾ + a ⋅ ∇ ξ h ‾ − f ‾ ξ ⋅ a = − 1 τ h ( h ‾ − h ‾ ( e q ) ) + Z τ h f ( f ‾ − f ‾ ( e q ) ) . (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} ∂ t ∂ h + ξ ⋅ ∇ h + a ⋅ ∇ ξ h − f ξ ⋅ a = − τ h 1 ( h − h ( e q ) ) + τ h f Z ( f − f ( e q ) ) . ( 3 . 1 0 )
它们的平衡态分别为
f ‾ ( e q ) = ∫ f ~ ( e q ) d η = ρ ( 2 π R g T ) D 2 exp [ − ( ξ − u ) 2 2 R g T ] , (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}
f ( e q ) = ∫ f ( e q ) d η = ( 2 π R g T ) 2 D ρ exp [ − 2 R g T ( ξ − u ) 2 ] , ( 3 . 1 1 )
h ‾ ( e q ) = ∫ ξ 2 + η 2 2 f ~ ( e q ) d η = ρ ( ξ 2 + K R g T ) 2 ⋅ ( 2 π R g T ) D 2 exp [ − ( ξ − u ) 2 2 R g T ] . (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}
h ( e q ) = ∫ 2 ξ 2 + η 2 f ( e q ) d η = 2 ⋅ ( 2 π R g T ) 2 D ρ ( ξ 2 + K R g T ) exp [ − 2 R g T ( ξ − u ) 2 ] . ( 3 . 1 2 )
可以看到多原子分子的流场与单原子分子的区别主要体现在能量上, 而 f ‾ ( e q ) {\overline{f}}^{(eq)} f ( e q )
ρ = ∫ f ‾ d ξ , j = ρ u = ∫ ξ f ‾ d ξ , ρ E = ∫ h ‾ d ξ . \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}.
ρ = ∫ f d ξ , j = ρ u = ∫ ξ f d ξ , ρ E = ∫ h d ξ .
定容比热比 c v c_{v} c v c p c_{p} c p
c v = D + K 2 R g , c p = D + K + 2 2 R g . c_{v} = \frac{D + K}{2}R_g ,\qquad
c_{p} = \frac{D + K + 2}{2}R_g .
c v = 2 D + K R g , c p = 2 D + K + 2 R g .
对单原子分子而言, η = 0 \mathbf{\eta = 0\ } η = 0 K = 0 K = 0 K = 0
式(3.9)和式(3.10)在$\mathbf{c}_i $方向上的方程为
∂ f ‾ i ∂ t + c i ⋅ ∇ f ‾ i = − 1 τ f ( f ‾ i − f ‾ i ( e q ) ) + F i , (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}
∂ t ∂ f i + c i ⋅ ∇ f i = − τ f 1 ( f i − f i ( e q ) ) + F i , ( 3 . 1 3 )
∂ h ‾ i ∂ t + c i ⋅ ∇ h ‾ i = − 1 τ h ( h ‾ i − h ‾ i ( e q ) ) + Z i τ h f ( f ‾ i − f ‾ i ( e q ) ) + q i . (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}
∂ t ∂ h i + c i ⋅ ∇ h i = − τ h 1 ( h i − h i ( e q ) ) + τ h f Z i ( f i − f i ( e q ) ) + q i . ( 3 . 1 4 )
其中 Z i = c i ⋅ u − u 2 2 Z_i = \mathbf{c}_i \mathbf{\cdot u} - \frac{u^{2}}{2} Z i = c i ⋅ u − 2 u 2
f ‾ i ( x + c i Δ t , t + Δ t ) − f ‾ i ( x , t ) = − ω f ( f ‾ i − f ‾ i ( e q ) ) + ( 1 − ω f 2 ) F i Δ 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}
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − ω f ( f i − f i ( e q ) ) + ( 1 − 2 ω f ) F i Δ t , ( 3 . 1 5 )
h ‾ i ( x + c i Δ t , t + Δ t ) − h ‾ i ( x , t ) = − Ω h ( h ‾ i − h ‾ i ( e q ) ) + Z i τ h f ( f ‾ i − f ‾ i ( e q ) ) + q i . (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}
h i ( x + c i Δ t , t + Δ t ) − h i ( x , t ) = − Ω h ( h i − h i ( e q ) ) + τ h f Z i ( f i − f i ( e q ) ) + q i . ( 3 . 1 6 )
Guo等 将平衡态使用Hermite展开对 f ‾ ( e q ) {\overline{f}}^{(eq)} f ( e q ) h ‾ ( e q ) {\overline{h}}^{(eq)} h ( e q ) T T T T 0 T_0 T 0
f ‾ i ( e q ) = ρ w i [ 1 + c i ⋅ u R g T 0 + 1 2 ( c i ⋅ u R g T 0 ) 2 − u 2 2 R g T 0 ] , (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}
f i ( e q ) = ρ w i [ 1 + R g T 0 c i ⋅ u + 2 1 ( R g T 0 c i ⋅ u ) 2 − 2 R g T 0 u 2 ] , ( 3 . 1 7 )
h ‾ i ( e q ) = p 0 w i [ c i ⋅ u R g T 0 + ( c i ⋅ u R g T 0 ) 2 − u 2 R g T 0 + 1 2 ( c i 2 R g T 0 − D ) ] + f ‾ i ( e q ) 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}
h i ( e q ) = p 0 w i [ R g T 0 c i ⋅ u + ( R g T 0 c i ⋅ u ) 2 − R g T 0 u 2 + 2 1 ( R g T 0 c i 2 − D ) ] + f i ( e q ) E . ( 3 . 1 8 )
源项为：
F i = ρ w i [ c i ⋅ a R g T 0 + ( c i ⋅ u ) ( c i ⋅ a ) ( R g T 0 ) 2 − a ⋅ u R g T 0 ] , q i = ( ρ w i E R g T 0 + f i ) ( c i ⋅ a ) . (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}
F i = ρ w i [ R g T 0 c i ⋅ a + ( R g T 0 ) 2 ( c i ⋅ u ) ( c i ⋅ a ) − R g T 0 a ⋅ u ] , q i = ( R g T 0 ρ w i E + f i ) ( c i ⋅ a ) . ( 3 . 1 9 )
各宏观量则为：
ρ = ∑ i f ‾ i , j = ρ u = ∑ i c i f ‾ i , ρ E = ∑ i h ‾ i . \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.
ρ = i ∑ f i , j = ρ u = i ∑ c i f i , ρ E = i ∑ h i .
运动粘度为μ = τ f p 0 \mu = \tau_{f}p_{0} μ = τ f p 0 κ = c v τ h p 0 = D + K 2 R g τ h p 0 \kappa = c_{v}\tau_{h}p_{0} = \frac{D + K}{2}R_g \tau_{h}p_{0} κ = c v τ h p 0 = 2 D + K R g τ h p 0 Pr = μ c p / κ = γ τ f / τ h \Pr = \mu c_{p}/\kappa = \gamma\tau_{f}/\tau_{h} Pr = μ c p / κ = γ τ f / τ h γ \gamma γ
[NOTE] p 0 = ρ R g T 0 p_{0} = \rho R_g T_{0} p 0 = ρ R g T 0
何雅玲团队在Guo等 的理论基础上, 提出了一种耦合的DDF-LBE, 实现了热流和流场的耦合模拟.这里主要介绍耦合DDF-LBE与基于总能的DDF-LBE的区别.
为简化讨论, 考虑无外力项耦合DDF-LBE在 c i \mathbf{c}_i c i
∂ f i ∂ t + c i ⋅ ∇ f i = − 1 τ f ( f i − f i ( e q ) ) , (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}
∂ t ∂ f i + c i ⋅ ∇ f i = − τ f 1 ( f i − f i ( e q ) ) , ( 3 . 2 0 )
∂ h i ∂ t + c i ⋅ ∇ h i = − 1 τ h ( h i − h i ( e q ) ) + ( c i ⋅ u ) τ h f ( f i − f i ( e q ) ) . (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}
∂ t ∂ h i + c i ⋅ ∇ h i = − τ h 1 ( h i − h i ( e q ) ) + τ h f ( c i ⋅ u ) ( f i − f i ( e q ) ) . ( 3 . 2 1 )
式(3.21)中将Z i = c i ⋅ u − u 2 2 Z_i = \mathbf{c}_{\mathbf{i}}\mathbf{\cdot u} - \frac{u^{2}}{2} Z i = c i ⋅ u − 2 u 2 ( c i ⋅ u ) \left( \mathbf{c}_{\mathbf{i}}\mathbf{\cdot u} \right) ( c i ⋅ u ) f α ( e q ) f_{\alpha}^{(eq)} f α ( e q ) h α ( e q ) h_{\alpha}^{(eq)} h α ( e q )
f α ( e q ) = A i + B σ ⋅ ( c α ⋅ u ) + C α ∥ u ∥ 2 + D α ⋅ ( c α ⋅ u ) 2 + E α ∥ u ∥ 2 ⋅ ( 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}
f α ( e q ) = A i + B σ ⋅ ( c α ⋅ u ) + C α ∥ u ∥ 2 + D α ⋅ ( c α ⋅ u ) 2 + E α ∥ u ∥ 2 ⋅ ( c α ⋅ u ) + F α ⋅ ( c α ⋅ u ) 3 , ( 3 . 2 2 )
h α ( e q ) = K α + L α ⋅ ( c α ⋅ u ) + M α ∥ u ∥ 2 + 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}
h α ( e q ) = K α + L α ⋅ ( c α ⋅ u ) + M α ∥ u ∥ 2 + N α ⋅ ( c α ⋅ u ) 2 . ( 3 . 2 3 )
则各系数可通过式(3.11)和式(3.12)的Chapman-Enskog展开得出.
第二种方式则采用其他 f ( e q ) f^{(eq)} f ( e q ) f α ( e q ) f_{\alpha}^{(eq)} f α ( e q )
h ( e q ) = [ E + ( ξ − u ) ⋅ u ] f ( e q ) + ω ( ξ , T ) h^{(eq)} = \left[ E + \left( \mathbf{\xi} - \mathbf{u} \right) \cdot \mathbf{u} \right] f^{(eq)} + \omega(\mathbf{\xi},T)
h ( e q ) = [ E + ( ξ − u ) ⋅ u ] f ( e q ) + ω ( ξ , T )
其中 ω ( ξ , T ) \omega\left( \mathbf{\xi},T \right) ω ( ξ , T )
∫ ω ( ξ , T ) d ξ = 0 , ∫ ω ( ξ , T ) ξ d ξ = 0 , ∫ ω ( ξ , T ) ξ ξ T d ξ = p R g T I \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}
∫ ω ( ξ , T ) d ξ = 0 , ∫ ω ( ξ , T ) ξ d ξ = 0 , ∫ ω ( ξ , T ) ξ ξ T d ξ = p R g T I
因此可通过构建等价离散形式(式(3.24))：
h α ( e q ) = w α p c 2 R g T + [ E + ( c α − u ) ⋅ u ] ⋅ f α ( e q ) (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}
h α ( e q ) = c 2 w α p R g T + [ E + ( c α − u ) ⋅ u ] ⋅ f α ( e q ) ( 3 . 2 4 )
计算总能的离散平衡态.
对于二维问题, 基于Qu等提出的圆函数进行修改：
f ( e q ) = { ρ 2 π c ∥ ξ − u ∥ = c ≡ D ( γ − 1 ) e and λ = e p , 0 otherwise (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}
f ( e q ) = { 2 π c ρ 0 ∥ ξ − u ∥ = c ≡ D ( γ − 1 ) e and λ = e p , otherwise ( 3 . 2 5 )
其中 λ \lambda λ e p = [ 1 − D 2 ( γ − 1 ) ] e e_{p} = \left[ 1 - \frac{D}{2}(\gamma - 1) \right] e e p = [ 1 − 2 D ( γ − 1 ) ] e f ( e q ) f^{(eq)} f ( e q )
f ( e q ) = { ρ 4 π r 2 ∥ ξ − u ∥ = ∥ r ∥ = r , 0 otherwise . (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}
f ( e q ) = { 4 π r 2 ρ 0 ∥ ξ − u ∥ = ∥ r ∥ = r , otherwise . ( 3 . 2 6 )
以Li等使用的D3Q27模型为例, 可得ρ r 2 / 3 = p \rho r^{2}/3 = p ρ r 2 / 3 = p r = 3 R g T r = \sqrt{3 R_g T} r = 3 R g T
此外, Li等指出：空间离散采用五阶WENO和TVD(total variation diminishing scheme)格式捕捉可压缩流中的不连续性；时间离散可采用二阶(或三阶)IMEX Runge-Kutta格式推进.虽然式(3.20)和式(3.21)同样可被离散为如下形式(t t t t + Δ t t + \Delta t t + Δ t
f α t + Δ t = f α t − Δ t ⋅ ( c α ⋅ ∇ f α t ) + Δ t τ f t ( f α ( e q ) , t − f α 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}
f α t + Δ t = f α t − Δ t ⋅ ( c α ⋅ ∇ f α t ) + τ f t Δ t ( f α ( e q ) , t − f α t ) ( 3 . 2 7 )
h α t + Δ t = h α t + Δ t [ 1 τ h t ( h α ( e q ) , t − h α t ) − ( c α ⋅ u t ) τ h f t ( f α ( e q ) , t − f α 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}
h α t + Δ t = h α t + Δ t [ τ h t 1 ( h α ( e q ) , t − h α t ) − τ h f t ( c α ⋅ u t ) ( f α ( e q ) , t − f α t ) − ( c α ⋅ ∇ h α t ) ] ( 3 . 2 8 )
但这种简单的形式仅适用于流场中无激波、无间断的情况, 且在时间上仅有一阶精度.
