本文同步发表于知乎 和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 + \boldsymbol{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 \boldsymbol{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} = -\boldsymbol{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 源项的展开结果。
这里的推导主要参考 Li 等的文章，并补充了一些自己的笔记。
需要提醒的是，下文的速度记号和前文不同。
根据 LBM 格子张量对称性，有前 5 阶表达式：
∑ α w α = 1 , ∑ α w α e α i = 0 , ∑ α w α e α i e α j = c s 2 δ i j , ∑ α w α e α i e α j e α k = 0 ∑ α w α e α i e α j e α k e α l = c s 4 ( δ i j δ k l + δ i k δ j l + δ i l δ k j ) ∑ α w α e α i e α j e α k e α l e α 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
α ∑ w α = 1 , α ∑ w α e α i = 0 , α ∑ w α e α i e α j = c s 2 δ i j , α ∑ w α e α i e α j e α k = 0 α ∑ w α e α i e α j e α k e α l = c s 4 ( δ i j δ k l + δ i k δ j l + δ i l δ k j ) α ∑ w α e α i e α j e α k e α l e α m = 0
其中 e α i e_{\alpha i} e α i D 维向量 e α \bold{e}_{\alpha} e α i i i w α w_\alpha w α e α \bold{e}_{\alpha} e α δ i j \delta_{ij} δ i j
所以，对于 EDM 格式的源项 F α , E D M F_{\alpha,\mathrm{EDM}} F α , E D M
∑ α F α , E D M = 0 , ∑ α e α F α , E D M = F , ∑ α e α e α F α , E D M = δ t ( ( u ∗ F + F u ∗ ) + δ t F F ρ ) (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}
α ∑ F α , E D M = 0 , α ∑ e α F α , E D M = F , α ∑ e α e α F α , E D M = δ t ( ( u ∗ F + F u ∗ ) + δ t ρ F F ) ( 9 )
技术细节可见附录（A）。
令展开参数 ϵ \epsilon ϵ
f α = f α ( 0 ) + ϵ f α ( 1 ) + ϵ 2 f α ( 2 ) + . . . , f_{\alpha} = f_{\alpha}^{(0)} + \epsilon f_{\alpha}^{(1)} + \epsilon^2 f_{\alpha}^{(2)} + ...,
f α = f α ( 0 ) + ϵ f α ( 1 ) + ϵ 2 f α ( 2 ) + . . . ,
源项展开为 F α = ϵ F 1 α F_{\alpha} = \epsilon F_{1 \alpha} F α = ϵ F 1 α
∂ ∂ t = ϵ ∂ ∂ t 1 + ϵ 2 ∂ ∂ t 2 , ∂ ∂ x = ϵ ∂ ∂ x 1 . \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}.
∂ t ∂ = ϵ ∂ t 1 ∂ + ϵ 2 ∂ t 2 ∂ , ∂ x ∂ = ϵ ∂ x 1 ∂ .
将式(4)中的 f α ( x + e α δ t , t + δ t ) f_\alpha (\bold{x} + \bold{e}_\alpha \delta_t, t+\delta_t) f α ( x + e α δ t , t + δ t ) f α ( x , t ) f_\alpha (\bold{x}, t) f α ( x , t ) ϵ \epsilon ϵ
O ( ϵ 0 ) : f α ( 0 ) = f α ( e q ) , ( 10 − 1 ) O ( ϵ 1 ) : D 1 α f α ( 0 ) = − 1 τ δ t f α ( 1 ) + F 1 α δ t , ( 10 − 2 ) O ( ϵ 2 ) : ∂ f α ( 0 ) ∂ t 2 + D 1 α f α ( 1 ) + δ t 2 D 1 α 2 f α ( 0 ) = − 1 τ δ t f α ( 2 ) . ( 10 − 3 ) \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}
O ( ϵ 0 ) : O ( ϵ 1 ) : O ( ϵ 2 ) : f α ( 0 ) = f α ( e q ) , D 1 α f α ( 0 ) = − τ δ t 1 f α ( 1 ) + δ t F 1 α , ∂ t 2 ∂ f α ( 0 ) + D 1 α f α ( 1 ) + 2 δ t D 1 α 2 f α ( 0 ) = − τ δ t 1 f α ( 2 ) . ( 1 0 − 1 ) ( 1 0 − 2 ) ( 1 0 − 3 )
其中 D 1 α = ∂ ∂ t 1 + e α ⋅ ∂ ∂ x 1 \displaystyle \mathrm{D}_{1 \alpha} = \frac{\partial}{\partial t_1} + \bold{e}_{\alpha} \cdot \frac{\partial}{\partial \bold{x}_1} D 1 α = ∂ t 1 ∂ + e α ⋅ ∂ x 1 ∂
O ( ϵ 2 ) : ∂ f α ( 0 ) ∂ t 2 + D 1 α ( 1 − 1 2 τ ) f α ( 1 ) + 1 2 D 1 α F 1 α = − 1 τ δ t f α ( 2 ) . ( 10 − 4 ) \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}
O ( ϵ 2 ) : ∂ t 2 ∂ f α ( 0 ) + D 1 α ( 1 − 2 τ 1 ) f α ( 1 ) + 2 1 D 1 α F 1 α = − τ δ t 1 f α ( 2 ) . ( 1 0 − 4 )
（1）考虑对式(10-3)和式(10-4)计算零阶矩 ，得：
∂ ρ ∂ t 1 + ∇ 1 ⋅ ( ρ u ∗ ) = 0 ∂ ρ ∂ t 2 + 1 2 ∇ 1 ⋅ ( δ t F 1 ) = 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 \\
∂ t 1 ∂ ρ + ∇ 1 ⋅ ( ρ u ∗ ) = 0 ∂ t 2 ∂ ρ + 2 1 ∇ 1 ⋅ ( δ t F 1 ) = 0
此处 u ∗ = ( ∑ i e α f α ) / ( ∑ i f α ) \bold{u}^* = (\sum_{i} \bold{e}_\alpha f_\alpha) / (\sum_{i} f_\alpha) u ∗ = ( ∑ i e α f α ) / ( ∑ i f α ) u ^ = u ∗ + F δ t / ( 2 ρ ) \hat{\bold{u}} = \bold{u}^* + \bold{F} \delta_t / (2 \rho) u ^ = u ∗ + F δ t / ( 2 ρ )
基于前面假设的偏微分关系式，我们可以得到 ∂ ρ ∂ t + ∇ ⋅ ( ρ u ^ ) = 0 \displaystyle \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \hat{\bold{u}}) = 0 ∂ t ∂ ρ + ∇ ⋅ ( ρ u ^ ) = 0
（2）同理，对式(10-3)和式(10-4)计算一阶矩 ，得：
∂ ρ u ∗ ∂ t 1 + ∇ 1 ⋅ ( ρ u ∗ u ∗ ) = − ∇ 1 p + F 1 ∂ ρ u ∗ ∂ t 2 + ∇ 1 ⋅ [ ( 1 − 1 2 τ ) ∑ α e α e α f α ( 1 ) ] + δ t 2 ∂ F 1 ∂ t 1 + 1 2 ∇ 1 ⋅ ( ∑ α e α e α F 1 α ) = 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
∂ t 1 ∂ ρ u ∗ + ∇ 1 ⋅ ( ρ u ∗ u ∗ ) = − ∇ 1 p + F 1 ∂ t 2 ∂ ρ u ∗ + ∇ 1 ⋅ [ ( 1 − 2 τ 1 ) α ∑ e α e α f α ( 1 ) ] + 2 δ t ∂ t 1 ∂ F 1 + 2 1 ∇ 1 ⋅ ( α ∑ e α e α F 1 α ) = 0
其中 p = ρ c s 2 p = \rho c_s^2 p = ρ c s 2
根据式(10-2)，有：
∑ α e α e α f α ( 1 ) = − τ [ δ t ∂ ∂ t 1 ( ∑ α e α e α f α ( 0 ) ) + δ t ∇ 1 ⋅ ( ∑ α e α e α e α f α ( 0 ) ) − ∑ α e α e α F 1 α ] \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}
α ∑ e α e α f α ( 1 ) = − τ [ δ t ∂ t 1 ∂ ( α ∑ e α e α f α ( 0 ) ) + δ t ∇ 1 ⋅ ( α ∑ e α e α e α f α ( 0 ) ) − α ∑ e α e α F 1 α ]
记
Π ( 1 ) = − δ t ( τ − 1 2 ) [ ∂ ∂ t 1 ( ∑ α e α e α f α ( 0 ) ) + ∇ 1 ⋅ ( ∑ α e α e α e α f α ( 0 ) ) ] = − δ t ( τ − 1 2 ) { ρ c s 2 [ ∇ 1 u ∗ + ( ∇ 1 u ∗ ) T ] + u ∗ F 1 + F 1 u ∗ } \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}
Π ( 1 ) = − δ t ( τ − 2 1 ) [ ∂ t 1 ∂ ( α ∑ e α e α f α ( 0 ) ) + ∇ 1 ⋅ ( α ∑ e α e α e α f α ( 0 ) ) ] = − δ t ( τ − 2 1 ) { ρ c s 2 [ ∇ 1 u ∗ + ( ∇ 1 u ∗ ) T ] + u ∗ F 1 + F 1 u ∗ }
则式(10-4)的一阶矩 表示为： ∂ ρ u ∗ ∂ t 2 + ∇ 1 Π ( 1 ) + δ t 2 ∂ F 1 ∂ t 1 + ∇ 1 ⋅ ( τ ∑ α e α e α F 1 α ) = 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 ∂ t 2 ∂ ρ u ∗ + ∇ 1 Π ( 1 ) + 2 δ t ∂ t 1 ∂ F 1 + ∇ 1 ⋅ ( τ α ∑ e α e α F 1 α ) = 0
那么同样根据假设的微分关系，即可还原出一阶矩对应的方程。这里直接写出使用真实流场速度 u ^ = u ∗ + F δ t / ( 2 ρ ) \hat{\bold{u}} = \bold{u}^* + \bold{F} \delta_t / (2 \rho) u ^ = u ∗ + F δ t / ( 2 ρ )
∂ ( ρ u ^ ) ∂ t + ∇ ⋅ ( ρ u ^ u ^ ) = − ∇ p + F + ∇ ⋅ { ρ ν [ ∇ 1 u ∗ + ( ∇ 1 u ∗ ) T } − ∇ ⋅ { ρ ν [ ∇ 1 a + ( ∇ 1 a ) T ] } + δ t 2 ϵ 2 ∂ F ∂ t 2 + ∇ ⋅ [ τ δ t ( u ∗ F + F u ∗ ) + δ t 2 F F 4 ρ − τ ∑ α 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}
∂ t ∂ ( ρ u ^ ) + ∇ ⋅ ( ρ u ^ u ^ ) = − ∇ p + F + ∇ ⋅ { ρ ν [ ∇ 1 u ∗ + ( ∇ 1 u ∗ ) T } − ∇ ⋅ { ρ ν [ ∇ 1 a + ( ∇ 1 a ) T ] } + 2 δ t ϵ 2 ∂ t 2 ∂ F + ∇ ⋅ [ τ δ t ( u ∗ F + F u ∗ ) + δ t 2 4 ρ F F − τ α ∑ e α e α F α ]
其中 a = δ t F / 2 ρ \bold{a} = \delta_t \bold{F} / 2 \rho a = δ t F / 2 ρ ν = ( τ − 1 2 ) c s 2 δ t \nu = (\tau - \frac{1}{2}) c_s^2 \delta_t ν = ( τ − 2 1 ) c s 2 δ t
Err = − ∇ ⋅ { ρ ν [ ∇ 1 a + ( ∇ 1 a ) T ] } + δ t 2 ϵ 2 ∂ F ∂ t 2 + ∇ ⋅ [ τ δ t ( u ∗ F + F u ∗ ) + δ t 2 F F 4 ρ − τ ∑ α 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]
Err = − ∇ ⋅ { ρ ν [ ∇ 1 a + ( ∇ 1 a ) T ] } + 2 δ t ϵ 2 ∂ t 2 ∂ F + ∇ ⋅ [ τ δ t ( u ∗ F + F u ∗ ) + δ t 2 4 ρ F F − τ α ∑ e α e α F α ]
但上述的展开方式和Wagner等使用Taylor展开的结果并不一致。 Li 等 指出这种错误的来源是：上文将 f α ( 0 ) = f α ( e q ) ( ρ , u ∗ ) f_{\alpha}^{(0)} = f_{\alpha}^{(eq)}(\rho, \bold{u}^*) f α ( 0 ) = f α ( e q ) ( ρ , u ∗ ) ∑ α e α f α ( 1 ) = 0 \sum_{\alpha} \bold{e}_{\alpha} f_{\alpha}^{(1)} = 0 ∑ α e α f α ( 1 ) = 0 f α ( 0 ) = f α ( e q ) ( ρ , u ^ ) f_{\alpha}^{(0)} = f_{\alpha}^{(eq)}(\rho, \hat{\bold{u}}) f α ( 0 ) = f α ( e q ) ( ρ , u ^ ) ∑ α e α f α ( 1 ) = − δ t F 1 / 2 \sum_{\alpha} \bold{e}_{\alpha} f_{\alpha}^{(1)} = -\delta_t \bold{F}_1 / 2 ∑ α e α f α ( 1 ) = − δ t F 1 / 2
因此，需要将原本的 LBE 等效地修改为：
f α ( x + e α δ t , t + δ t ) − f α ( x , t ) = F ^ α − 1 τ [ f α ( x , t ) − f α ( e q ) ( ρ , 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 α ( x + e α δ t , t + δ t ) − f α ( x , t ) = F ^ α − τ 1 [ f α ( x , t ) − f α ( e q ) ( ρ , u ^ ) ] ( 1 1 )
其中 F ^ α = F α − 1 τ ( f α ( e q ) ( ρ , u ^ ) − f α ( e q ) ( ρ , u ∗ ) ) \hat{F}_{\alpha} = F_{\alpha} - \frac{1}{\tau} (f_\alpha^{(eq)} (\rho, \hat{\bold{u}}) - f_\alpha^{(eq)} (\rho, \bold{u}^*)) F ^ α = F α − τ 1 ( f α ( e q ) ( ρ , u ^ ) − f α ( e q ) ( ρ , u ∗ ) ) F α F_{\alpha} F α u ^ \hat{\bold{u}} u ^
∂ ( ρ u ^ ) ∂ t + ∇ ⋅ ( ρ u ^ u ^ ) = − ∇ p + F + ∇ ⋅ { ρ ν [ ∇ 1 u ^ + ( ∇ 1 u ^ ) T } + ∇ ⋅ [ ( τ − 1 2 ) δ t ( u ^ F + F u ^ ) − τ ∑ α 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}
∂ t ∂ ( ρ u ^ ) + ∇ ⋅ ( ρ u ^ u ^ ) = − ∇ p + F + ∇ ⋅ { ρ ν [ ∇ 1 u ^ + ( ∇ 1 u ^ ) T } + ∇ ⋅ [ ( τ − 2 1 ) δ t ( u ^ F + F u ^ ) − τ α ∑ e α e α F ^ α ] ( 1 2 )
式(12)的误差项为：
Err = ∇ ⋅ [ ( τ − 1 2 ) δ t ( u ^ F + F u ^ ) − τ ∑ α 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}
Err = ∇ ⋅ [ ( τ − 2 1 ) δ t ( u ^ F + F u ^ ) − τ α ∑ e α e α F ^ α ] ( 1 3 )
将 EDM 的源项代入式(13)中，得到： Err = − δ t 2 ∇ ⋅ ( F F 4 ρ ) \text{Err}=-\delta_t^2 \nabla\cdot\left( \frac{\bold{FF}}{4 \rho} \right) Err = − δ t 2 ∇ ⋅ ( 4 ρ F F )
这里仅以 ∑ α e α e α F α , E D M \displaystyle \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} α ∑ e α e α F α , E D M
∑ α e α e α F α , E D M = ∑ α e α e α ⋅ δ 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 ] ) ] = δ t c s 2 ∑ α [ e α e α ( w α e α ⋅ F ) ] + δ t 2 c s 4 ∑ α w α e α e α ⋅ [ ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ 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 α F α , E D M = α ∑ e α e α ⋅ δ 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 ] ) ] = c s 2 δ t α ∑ [ e α e α ( w α e α ⋅ F ) ] + 2 c s 4 δ t α ∑ w α e α e α ⋅ [ ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ I ] ) ]
一方面，对于 ∑ α [ e α e α ( w α e α ⋅ F ) ] \displaystyle \sum_{\alpha} \left[ \bold{e}_\alpha \bold{e}_\alpha (w_\alpha \bold{e}_\alpha \cdot \bold{F}) \right] α ∑ [ e α e α ( w α e α ⋅ F ) ] i i i j j j
∑ k F k ∑ α ( w α e α i e α j e α k ) = ∑ k F k ⋅ 0 = 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
k ∑ F k α ∑ ( w α e α i e α j e α k ) = k ∑ F k ⋅ 0 = 0
另一方面，已知 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 ( v E D M F + F v E D M ) \left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) ( v E D M F + F v E D M )
[ B ] = v E D M F + F v E D M = u ∗ F + F u ∗ + δ t ρ F F [\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 ] = v E D M F + F v E D M = u ∗ F + F u ∗ + ρ δ t F F
则 [ B ] : ( e α e α − c s 2 [ I ] ) [\bold{\mathbb{B}}] : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) [ B ] : ( e α e α − c s 2 [ I ] ) [ A ] = ∑ α w α e α e α ⋅ { [ B ] : ( e α e α − c s 2 [ 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}])\} [ A ] = α ∑ w α e α e α ⋅ { [ B ] : ( e α e α − c s 2 [ I ] ) } i i i j j j A i j \mathbb{A}_{ij} A i j
A i j = ∑ α w α e α i e α j [ ∑ k ∑ l B k l ( e α k e α l − c s 2 δ k l ) ] = ∑ k ∑ l B k l ∑ α { w α e α i e α j ( e α k e α l − c s 2 δ k l ) } = ∑ k ∑ l B k l ⋅ c s 4 ( δ i l δ j k + δ i k δ j l ) = ( B i j + B j i ) c s 4 \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}
A i j = α ∑ w α e α i e α j [ k ∑ l ∑ B k l ( e α k e α l − c s 2 δ k l ) ] = k ∑ l ∑ B k l α ∑ { w α e α i e α j ( e α k e α l − c s 2 δ k l ) } = k ∑ l ∑ B k l ⋅ c s 4 ( δ i l δ j k + δ i k δ j l ) = ( B i j + B j i ) c s 4
因为 [ B ] [\mathbb{B}] [ B ] A i j = 2 c s 4 B i j \mathbb{A}_{ij} =2 c_s^4 \mathbb{B}_{ij} A i j = 2 c s 4 B i j
∑ α e α e α F α , E D M = 0 ⃗ + δ t 2 c s 4 ∑ α w α e α e α ⋅ [ ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ I ] ) ] = δ t ( v E D M F + F v E D M ) = δ t ( u ∗ F + F u ∗ + δ t ρ F F ) \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 α , E D M = 0 + 2 c s 4 δ t α ∑ w α e α e α ⋅ [ ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ I ] ) ] = δ t ( v E D M F + F v E D M ) = δ t ( u ∗ F + F u ∗ + ρ δ t F F )
综上， ∑ α e α e α F α , E D M = δ t ( u ∗ F + F u ∗ + δ t F F ρ ) \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) α ∑ e α e α F α , E D M = δ t ( u ∗ F + F u ∗ + δ t ρ F F )
由于 a = δ t F / 2 ρ = u ^ − u ∗ \bold{a} = \delta_t \bold{F} / 2 \rho = \hat{\bold{u}} - \bold{u}^* a = δ t F / 2 ρ = u ^ − u ∗ ρ a = δ t F / 2 \rho \bold{a} = \delta_t \bold{F} / 2 ρ a = δ t F / 2
1 τ ( f α ( e q ) ( ρ , u ^ ) − f α ( e q ) ( ρ , u ∗ ) ) = ρ w α τ [ e α ⋅ a c s 2 + 1 2 c s 4 ( u ^ u ^ − u ∗ u ∗ ) : ( e α e α − c s 2 [ I ] ) ] = ρ w α τ [ e α ⋅ a c s 2 + 1 2 c s 4 ( u ∗ a + a u ∗ + a a ) : ( e α e α − c s 2 [ I ] ) ] = w α δ t τ [ e α ⋅ F 2 c s 2 + 1 2 c s 4 ( 1 2 ( u ∗ F + F u ∗ ) + δ t F F 4 ρ ) : ( e α e α − c s 2 [ 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}
= = = τ 1 ( f α ( e q ) ( ρ , u ^ ) − f α ( e q ) ( ρ , u ∗ ) ) τ ρ w α [ c s 2 e α ⋅ a + 2 c s 4 1 ( u ^ u ^ − u ∗ u ∗ ) : ( e α e α − c s 2 [ I ] ) ] τ ρ w α [ c s 2 e α ⋅ a + 2 c s 4 1 ( u ∗ a + a u ∗ + a a ) : ( e α e α − c s 2 [ I ] ) ] τ w α δ t [ 2 c s 2 e α ⋅ F + 2 c s 4 1 ( 2 1 ( u ∗ F + F u ∗ ) + δ t 4 ρ F F ) : ( e α e α − c s 2 [ I ] ) ]
由于 EDM 源项的表达式为：
F α , E D M = δ t w α [ e α ⋅ F c s 2 + 1 2 c s 4 ( u ^ F + F u ^ ) : ( e α e α − c s 2 [ 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 α , E D M = δ t w α [ c s 2 e α ⋅ F + 2 c s 4 1 ( u ^ F + F u ^ ) : ( e α e α − c s 2 [ I ] ) ]
因此，其对应修正源项为
F α , E D M − 1 τ ( f α ( e q ) ( ρ , u ^ ) − f α ( e q ) ( ρ , u ∗ ) ) = e α ⋅ F c s 2 ( 1 − 1 2 τ ) + w α δ t 2 c s 4 { [ ( 1 2 τ − 1 ) ( u ∗ F + F u ∗ ) − ( a F + F a ) − δ t 4 ρ τ F F ] : ( e α e α − c s 2 [ 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}
= F α , E D M − τ 1 ( f α ( e q ) ( ρ , u ^ ) − f α ( e q ) ( ρ , u ∗ ) ) c s 2 e α ⋅ F ( 1 − 2 τ 1 ) + 2 c s 4 w α δ t { [ ( 2 τ 1 − 1 ) ( u ∗ F + F u ∗ ) − ( a F + F a ) − 4 ρ τ δ t F F ] : ( e α e α − c s 2 [ I ] ) }
