本文同步发表于知乎 和CSDN
我在阅读其他文献时，看到有文献提到 Kupershtokh 的一篇关于 LBM 外力项的文章，他提出的格式在其他文章中又被称为 Exact difference method（EDM） 。
本想搜索这篇文章的PDF进行阅读，但最后只找到一个网站 介绍了部分信息。这份粗略笔记的内容，都是基于其他文献中对该方法的介绍而整理的。
密度分布函数 f f f
∂ f ∂ t + ξ ⋅ ∇ f + a ⋅ ∇ ξ f = − 1 τ f [ f − f ( e q ) ] (1) \frac{\partial f}{\partial t} + \bold{\xi} \cdot \nabla f + \bold{a} \cdot \nabla_{\bold{\xi}} f = -\frac{1}{\tau_f} \left[ f - f^{(eq)}\right]
\tag{1}
∂ t ∂ f + ξ ⋅ ∇ f + a ⋅ ∇ ξ f = − τ f 1 [ f − f ( e q ) ] ( 1 )
其中：
f ( e q ) = ρ ( 2 π R g T ) D / 2 exp ( − ( ξ − u ) 2 2 R g T ) (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}
f ( e q ) = ( 2 π R g T ) D / 2 ρ exp ( − 2 R g T ( ξ − u ) 2 ) ( 2 )
符号含义分别为
ξ \bold{\xi} ξ ρ \rho ρ u \bold{u} u R g R_g R g T T T a = d u / d t \bold{a}=\mathrm{d}\bold{u}/\mathrm{d}t a = d u / d t
当Knudsen数较小 时，Boltzmann方程的体力项可以基于f ( e q ) f^{(eq)} f ( e q )
∇ ξ f ≈ ∇ ξ f ( e q ) \nabla_{\bold{\xi}} f \approx \nabla_{\bold{\xi}} f^{(eq)}
∇ ξ f ≈ ∇ ξ f ( e q )
根据式(2)， f ( e q ) f^{(eq)} f ( e q ) C = ξ − u \bold{C} = \bold{\xi} - \bold{u} C = ξ − u
∂ f ( e q ) ∂ u = − ∂ f ( e q ) ∂ ξ = − ∇ ξ f ( e q ) \frac{\partial f^{(eq)}}{\partial \bold{u}} = -\frac{\partial f^{(eq)}}{\partial \bold{\xi}} = -\nabla_{\bold{\xi}} f^{(eq)}
∂ u ∂ f ( e q ) = − ∂ ξ ∂ f ( e q ) = − ∇ ξ f ( e q )
若假设密度 ρ \rho ρ ，则 d ρ d t = 0 \frac{\mathrm{d} \rho}{\mathrm{d} t}=0 d t d ρ = 0 f ( e q ) f^{(eq)} f ( e q )
d f ( e q ) d t = ∂ f ( e q ) ∂ u ⋅ d u d t + ∂ f ( e q ) ∂ ρ ⋅ d ρ d t ≈ ∂ f ( e q ) ∂ u ⋅ d u d t = − a ⋅ ∇ ξ f ( e q ) \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} = -\bold{a} \cdot \nabla_{\bold{\xi}} f^{(eq)}
\end{aligned}
d t d f ( e q ) = ∂ u ∂ f ( e q ) ⋅ d t d u + ∂ ρ ∂ f ( e q ) ⋅ d t d ρ ≈ ∂ u ∂ f ( e q ) ⋅ d t d u = − a ⋅ ∇ ξ f ( e q )
也就是说，Kupershtokh 等近似的 Boltzmann-BGK方程写作：
∂ f ∂ t + ξ ⋅ ∇ f = − 1 τ f [ f − f ( e q ) ] + d f ( e q ) d t (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}
∂ t ∂ f + ξ ⋅ ∇ f = − τ f 1 [ f − f ( e q ) ] + d t d f ( e q ) ( 3 )
将式(3)在[ t , t + δ t ] [t,t+\delta_t] [ t , t + δ t ] e α \bold{e}_\alpha e α
f α ( x + e α δ t , t + δ t ) − f α ( x , t ) = − 1 τ [ f α ( x , t ) − f α ( e q ) ( ρ , u ∗ ) ] + [ f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( 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}
f α ( x + e α δ t , t + δ t ) − f α ( x , t ) = − τ 1 [ f α ( x , t ) − f α ( e q ) ( ρ , u ∗ ) ] + [ f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) ] ( 4 )
其中：
δ t \delta_t δ t τ = τ f / δ t \tau = \tau_f / \delta_t τ = τ f / δ t x \bold{x} x ρ , u ∗ \rho, \bold{u^*} ρ , u ∗ t t t u \bold{u} u f α , f α ( e q ) f_\alpha, f_\alpha^{(eq)} f α , f α ( e q ) e α \bold{e}_\alpha e α 分布函数 和平衡态 ；
f α ( e q ) = w α ρ ⋅ [ 1 + e α ⋅ u c s 2 + ( e α ⋅ u ) 2 2 c s 4 − ∣ u ∣ 2 2 c s 2 ] = w α ρ ⋅ [ 1 + e α ⋅ u c s 2 + 1 2 c s 4 ( u u ) : ( e α e α − c s 2 [ 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}
f α ( e q ) = w α ρ ⋅ [ 1 + c s 2 e α ⋅ u + 2 c s 4 ( e α ⋅ u ) 2 − 2 c s 2 ∣ u ∣ 2 ] = w α ρ ⋅ [ 1 + c s 2 e α ⋅ u + 2 c s 4 1 ( u u ) : ( e α e α − c s 2 [ I ] ) ] ( 4 - 1 )
这里 w α w_\alpha w α e α \bold{e}_\alpha e α [ I ] [\bold{I}] [ I ] [ A ] [\bold{A}] [ A ] [ B ] [\bold{B}] [ B ] [ A ] : [ B ] = ∑ i = 1 m ∑ j = 1 n A i j B i j \displaystyle [\bold{A}] : [\bold{B}] = \sum^{m}_{i=1}\sum^{n}_{j=1} A_{ij} B_{ij} [ A ] : [ B ] = i = 1 ∑ m j = 1 ∑ n A i j B i j
F α F_\alpha F α F = ρ d u ∗ d t = ρ a \bold{F} = \rho \dfrac{\mathrm{d} \bold{u^*}}{\mathrm{d} t} = \rho \bold{a} F = ρ d t d u ∗ = ρ a
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) F_\alpha = f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t))
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) )
并且
u ∗ ( t ) = ( ∑ α f α e α ) / ( ∑ α f α ) u ∗ ( t + δ t ) = u ∗ ( t ) + ∫ t t + δ t d u ∗ d t d t = u ∗ ( t ) + ∫ t t + δ t F ρ d t ≈ u ∗ ( 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}
u ∗ ( t ) u ∗ ( t + δ t ) = ( α ∑ f α e α ) / ( α ∑ f α ) = u ∗ ( t ) + ∫ t t + δ t d t d u ∗ d t = u ∗ ( t ) + ∫ t t + δ t ρ F d t ≈ u ∗ ( t ) + ρ F δ t ( 5 )
这意味着加速度 a = F / ρ \bold{a}=\bold{F}/\rho a = F / ρ [ t , t + δ t ] [t,t+\delta_t] [ t , t + δ t ]
Kupershtokh等在这篇文章说明，在使用该源项时，流体宏观量的计算应表示为：
ρ = ∑ α f α , ρ u = ∑ i e α f α + δ t F 2 (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}
ρ = α ∑ f α , ρ u = i ∑ e α f α + 2 δ t F ( 6 )
Li 等通过 Chapman-Enskog 展开证明：该格式向宏观方程引入的误差项为：
E r r o r E D M = − δ t 2 ∇ ⋅ ( F F 4 ρ ) (7) \mathbf{Error}_{\mathrm{EDM}} = -\delta_t^2 \nabla\cdot\left( \frac{\bold{FF}}{4 \rho} \right)
\tag{7}
E r r o r E D M = − δ t 2 ∇ ⋅ ( 4 ρ F F ) ( 7 )
为便于分析，这里采用如下书写形式的平衡态：
f α ( e q ) = w α ρ ⋅ [ 1 + e α ⋅ u c s 2 + 1 2 c s 4 ( u u ) : ( e α e α − c s 2 [ 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]
f α ( e q ) = w α ρ ⋅ [ 1 + c s 2 e α ⋅ u + 2 c s 4 1 ( u u ) : ( e α e α − c s 2 [ I ] ) ]
将其带入 EDM 源项，记 Δ u = F ρ δ t \Delta\bold{u} = \frac{\bold{F}}{\rho} \delta_t Δ u = ρ F δ t
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) = ρ w α { e α ⋅ Δ u c s 2 + 1 2 c s 4 [ ( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ ] : ( e α e α − c s 2 [ 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}
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) = ρ w α { c s 2 e α ⋅ Δ u + 2 c s 4 1 [ ( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ ] : ( e α e α − c s 2 [ I ] ) }
这里涉及矩阵 ( u ∗ + Δ u ) ( u ∗ + Δ u ) (\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) ( u ∗ + Δ u ) ( u ∗ + Δ u ) u ∗ u ∗ \bold{u}^* \bold{u}^* u ∗ u ∗
( u ∗ + Δ u ) ( u ∗ + Δ u ) = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) , u ∗ u ∗ = u i ∗ u j ∗ (\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}
( u ∗ + Δ u ) ( u ∗ + Δ u ) = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) , u ∗ u ∗ = u i ∗ u j ∗
其中 u i ∗ u_i^* u i ∗ Δ u i = F i δ t / ρ \Delta u_i = F_i \delta_{t} / \rho Δ u i = F i δ t / ρ u ∗ \bold{u}^* u ∗ Δ u \Delta\bold{u} Δ u
( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) − u i ∗ u j ∗ = δ t ρ ( u i ∗ F j + u j ∗ F i ) + δ t 2 ρ 2 F i F j = δ t ρ ( u i ∗ + F i δ t 2 ρ ) F j + δ t ρ ( u j ∗ + F j δ t 2 ρ ) F i \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}
( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) − u i ∗ u j ∗ = ρ δ t ( u i ∗ F j + u j ∗ F i ) + ρ 2 δ t 2 F i F j = ρ δ t ( u i ∗ + 2 ρ F i δ t ) F j + ρ δ t ( u j ∗ + 2 ρ F j δ t ) F i
这里令 v E D M = u ∗ + F 2 ρ δ t \bold{v}_{\rm{EDM}} = \bold{u}^{*} + \frac{\bold{F}}{2 \rho} \delta_{t} v E D M = u ∗ + 2 ρ F δ t
F α = δ t w α [ e α ⋅ F c s 2 + 1 2 c s 4 ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ 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}
F α = δ t w α [ c s 2 e α ⋅ F + 2 c s 4 1 ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ I ] ) ] ( 8 )
上式即为 EDM 源项的展开结果。