最近看到了 Shi et al. [1] 的新文章很有意思: MesoFD 格式 —— 从平衡态就能去构建各种偏微分方程的有限差分LBM格式 (并包括其对应的源项和平衡态).
我这里主要是把自己感兴趣的 “Navier-Stokes 方程的构造” 部分内容摘录了出来, 并做分享 (没想到用Taylor展开来分析真的不难啊).
至于文章里的稳定性分析、DnQb构造等等, 还请各位阅读施保昌老师的原文.
[1] Shi B, Du R, Chai Z. Equilibrium-distribution-function-based mesoscopic finite-difference methods for partial differential equations: modeling and analysis . Physical Review E, 2026, 113(4)：045305. DOI:10.1103/61zg-dmhh. [Also available at arXiv:2505.17057 ]Multiple-relaxation-time finite-difference lattice Boltzmann model for the nonlinear convection-diffusion equation . Physical Review E, 2021, 104:035308. DOI:10.1103/physreve.104.035308
Shi et al. [1] 提出了一个 MesoFD 格式, 核心在于将 c j \boldsymbol{c}_j c j f j f_j f j f j e q f_j^{eq} f j e q j = 0 , 1 , … b − 1 j=0,1, \dots b-1 j = 0 , 1 , … b − 1
f j ( x , t ) = ∑ q ∑ k ( a k q f j e q ( x + k c j δ t , t + q δ t ) + b k q δ t F j ( x + k c j δ t , t + q δ t ) ) (1) f_j (\boldsymbol{x},t)=\sum_{q} \sum_{k} \left( a_{kq} f_j^{eq}(\boldsymbol{x}+k\boldsymbol{c}_j \delta_t, t+q\delta_t) + b_{kq} \delta_t F_j (\boldsymbol{x}+k\boldsymbol{c}_j \delta_t, t+q\delta_t) \right)
\tag{1}
f j ( x , t ) = q ∑ k ∑ ( a k q f j e q ( x + k c j δ t , t + q δ t ) + b k q δ t F j ( x + k c j δ t , t + q δ t ) ) ( 1 )
其中 a k q , b k q a_{kq}, b_{kq} a k q , b k q F j F_j F j − m q ≤ k ≤ m q -m_q \le k \le m_q − m q ≤ k ≤ m q − n ≤ q ≤ 0 -n \le q \le 0 − n ≤ q ≤ 0 n n n m q m_q m q t − δ t t-\delta_t t − δ t
为便于书写, Shi et al. 用了下面的简记符号
f j e q ( x + k c j δ t , t + q δ t ) = f j e q , [ k , q ] , F j ( x + k c j δ t , t + q δ t ) = F j [ k , q ] f_j^{eq}(\boldsymbol{x}+k\boldsymbol{c}_j \delta_t, t+q\delta_t) = f_j^{eq,[k,q]}
,\quad
F_j (\boldsymbol{x}+k\boldsymbol{c}_j \delta_t, t+q\delta_t) = F_j^{[k,q]}
f j e q ( x + k c j δ t , t + q δ t ) = f j e q , [ k , q ] , F j ( x + k c j δ t , t + q δ t ) = F j [ k , q ]
所以式(1)写为： f j ( x , t ) = ∑ k , q [ a k q f j e q , [ k , q ] + b k q δ t F j [ k , q ] ] f_j (\boldsymbol{x},t)=\sum_{k,q} [ a_{kq} f_j^{eq,[k,q]} + b_{kq} \delta_t F_j^{[k,q]}] f j ( x , t ) = ∑ k , q [ a k q f j e q , [ k , q ] + b k q δ t F j [ k , q ] ]
ρ ( x , t ) = ∑ j ∑ k , q [ a k q f j e q , [ k , q ] + b k q δ t F j [ k , q ] ] (2a) \rho(\boldsymbol{x},t) = \sum_{j} \sum_{k,q} \left[ a_{kq} f_j^{eq,[k,q]} + b_{kq} \delta_t F_j^{[k,q]} \right]
\tag{2a}
ρ ( x , t ) = j ∑ k , q ∑ [ a k q f j e q , [ k , q ] + b k q δ t F j [ k , q ] ] ( 2 a )
以及
( ρ u ) ( x , t ) = ∑ j ∑ k , q c j [ a ˉ k q f j e q , [ k , q ] + b ˉ k q δ t F j [ k , q ] ] (2b) (\rho\boldsymbol{u})(\boldsymbol{x},t) = \sum_{j}\sum_{k,q} \boldsymbol{c}_j \left[ \bar{a}_{kq} f_j^{eq,[k,q]} + \bar{b}_{kq} \delta_t F_j^{[k,q]} \right]
\tag{2b}
( ρ u ) ( x , t ) = j ∑ k , q ∑ c j [ a ˉ k q f j e q , [ k , q ] + b ˉ k q δ t F j [ k , q ] ] ( 2 b )
需要强调, 计算动量 ρ u \rho\boldsymbol{u} ρ u { a ˉ k q , b ˉ k q } \{\bar{a}_{kq}, \bar{b}_{kq}\} { a ˉ k q , b ˉ k q }
若对一个空间内的函数 g ( x ) g(\boldsymbol{x}) g ( x ) c j \boldsymbol{c}_j c j c j ˉ = − c j \boldsymbol{c}_{\bar{j}} = - \boldsymbol{c}_j c j ˉ = − c j
一阶空间中心差分算符： Δ j , k ( 1 ) g ( x ) = g ( x + k c j δ t ) − g ( x − k c j δ t ) \Delta_{j,k}^{(1)} g(\boldsymbol{x}) = g(\boldsymbol{x} + k \boldsymbol{c}_j \delta_t) - g(\boldsymbol{x} - k \boldsymbol{c}_j \delta_t) Δ j , k ( 1 ) g ( x ) = g ( x + k c j δ t ) − g ( x − k c j δ t ) Δ j ( 1 ) g ( x ) = Δ j , 1 ( 1 ) g ( x ) \Delta_{j}^{(1)} g(\boldsymbol{x}) = \Delta_{j,1}^{(1)} g(\boldsymbol{x}) Δ j ( 1 ) g ( x ) = Δ j , 1 ( 1 ) g ( x )
二阶空间中心差分算符： Δ j , k ( 2 ) g ( x ) = g ( x + k c j δ t ) − g ( x ) + g ( x − k c j δ t ) \Delta_{j,k}^{(2)} g(\boldsymbol{x}) = g(\boldsymbol{x} + k \boldsymbol{c}_j \delta_t) - g(\boldsymbol{x}) + g(\boldsymbol{x} - k \boldsymbol{c}_j \delta_t) Δ j , k ( 2 ) g ( x ) = g ( x + k c j δ t ) − g ( x ) + g ( x − k c j δ t ) Δ j ( 2 ) g ( x ) = Δ j , 1 ( 2 ) g ( x ) \Delta_{j}^{(2)} g(\boldsymbol{x}) = \Delta_{j,1}^{(2)} g(\boldsymbol{x}) Δ j ( 2 ) g ( x ) = Δ j , 1 ( 2 ) g ( x )
我们便可发现： g ( x ± k c j δ t ) = ( Δ j , k ( 2 ) ± Δ j , k ( 1 ) + g ( x ) ) / 2 g(\boldsymbol{x} \pm k \boldsymbol{c}_j \delta_t) = (\Delta_{j,k}^{(2)} \pm \Delta_{j,k}^{(1)} + g(\boldsymbol{x})) / 2 g ( x ± k c j δ t ) = ( Δ j , k ( 2 ) ± Δ j , k ( 1 ) + g ( x ) ) / 2
那么式(1)要进一步改写为：
f j ( x , t ) = ∑ k , q a k q [ 1 + ( Δ j , k ( 2 ) + Δ j , k ( 1 ) ) / 2 ] f j e q ( x , t + q δ t ) + δ t ∑ k , q b k q [ 1 + ( Δ j , k ( 2 ) + Δ j , k ( 1 ) ) / 2 ] F j ( x , t + q δ t ) \begin{aligned}
f_j (\boldsymbol{x}, t) =& \sum_{k,q} a_{kq} [1 + (\Delta_{j,k}^{(2)} + \Delta_{j,k}^{(1)})/2] f_j^{eq}(\boldsymbol{x}, t+q\delta_t) + \\&
\delta_t \sum_{k,q} b_{kq} [1 + (\Delta_{j,k}^{(2)} + \Delta_{j,k}^{(1)})/2] F_j (\boldsymbol{x}, t+q\delta_t)
\end{aligned}
f j ( x , t ) = k , q ∑ a k q [ 1 + ( Δ j , k ( 2 ) + Δ j , k ( 1 ) ) / 2 ] f j e q ( x , t + q δ t ) + δ t k , q ∑ b k q [ 1 + ( Δ j , k ( 2 ) + Δ j , k ( 1 ) ) / 2 ] F j ( x , t + q δ t )
这个算法其实也比较容易. 假定已有 t t t t + δ t t+\delta_t t + δ t
(1) 对网格点 x \boldsymbol{x} x f j e q ( x , t ) f_j^{eq}(\boldsymbol{x}, t) f j e q ( x , t ) F j ( x , t ) F_j(\boldsymbol{x}, t) F j ( x , t )
(2) 用式(1)计算所有点的 f j f_j f j
(3) 用式(2)更新宏观物理量.
(4) 施加边界条件.
下面讨论对式(1)直接 Taylor 展开的分析, 并最高展开至 N ≥ 1 N \geq 1 N ≥ 1
记 Taylor 展开算子
D k , q = q ∂ t + k ∇ ⋅ c j D_{k,q} = q \partial_{t} + k \nabla \cdot \boldsymbol{c}_{j}
D k , q = q ∂ t + k ∇ ⋅ c j
那么 f j e q , [ k , q ] f_j^{eq,[k,q]} f j e q , [ k , q ] F j [ k , q ] F_j^{[k,q]} F j [ k , q ]
f j e q ( x + k c j δ t , t + q δ t ) = ∑ l = 0 N ( δ t ) l l ! D k , q l f j e q ( x , t ) + O ( ( δ t ) N + 1 ) f_j^{eq}(\boldsymbol{x}+k\boldsymbol{c}_j \delta_t, t+q\delta_t) = \sum_{l=0}^{N} \frac{(\delta_{t})^l}{l!} D_{k,q}^l f_j^{eq}(\boldsymbol{x}, t) + O((\delta_{t})^{N+1})
f j e q ( x + k c j δ t , t + q δ t ) = l = 0 ∑ N l ! ( δ t ) l D k , q l f j e q ( x , t ) + O ( ( δ t ) N + 1 )
和
F j ( x + k c j δ t , t + q δ t ) = ∑ l = 0 N ( δ t ) l l ! D k , q l F j ( x , t ) + O ( ( δ t ) N + 1 ) F_j(\boldsymbol{x}+k\boldsymbol{c}_j \delta_t, t+q\delta_t) = \sum_{l=0}^{N} \frac{(\delta_{t})^l}{l!} D_{k,q}^l F_j(\boldsymbol{x}, t) + O((\delta_{t})^{N+1})
F j ( x + k c j δ t , t + q δ t ) = l = 0 ∑ N l ! ( δ t ) l D k , q l F j ( x , t ) + O ( ( δ t ) N + 1 )
考虑到算子 D k , q D_{k,q} D k , q
D k , q l = ( q ∂ t + k ∇ ⋅ c j ) l = ∑ m = 0 l C l m ( q ∂ t ) l − m ( k ∇ ⋅ c j ) m , C l m = l ! / ( m ! ( l − m ) ! ) \begin{aligned}
D_{k,q}^{l} &= (q \partial_{t} + k \nabla \cdot \boldsymbol{c}_{j})^l = \sum_{m=0}^{l} C_l^m \, (q \partial_{t})^{l-m} \, (k\nabla \cdot \boldsymbol{c}_j)^m ,\\
C_l^m &= l! / (m! (l-m)!)
\end{aligned}
D k , q l C l m = ( q ∂ t + k ∇ ⋅ c j ) l = m = 0 ∑ l C l m ( q ∂ t ) l − m ( k ∇ ⋅ c j ) m , = l ! / ( m ! ( l − m ) ! )
那么, 我们把和空间、时间有关的系数部分提取出来
A l m = C l m l ! ∑ k , q a k q q l − m k m , B l m = C l m l ! ∑ k , q b k q q l − m k m A_{lm} = \frac{C_l^m}{l!} \sum_{k,q} a_{kq} q^{l-m} k^{m} ,\quad
B_{lm} = \frac{C_l^m}{l!} \sum_{k,q} b_{kq} q^{l-m} k^{m}
A l m = l ! C l m k , q ∑ a k q q l − m k m , B l m = l ! C l m k , q ∑ b k q q l − m k m
则式(1)的展开要写作：
f j ( x , t ) = ∑ l = 0 N ( δ t ) l ∑ m = 0 l [ A l m ∂ t l − m ( ∇ ⋅ c j ) m f j e q + δ t B l m ∂ t l − m ( ∇ ⋅ c j ) m F j ] + O ( ( δ t ) N + 1 ) (3) f_j (\boldsymbol{x},t) = \sum_{l=0}^{N} (\delta_t)^l \sum_{m=0}^{l} \left[
A_{lm} \partial_{t}^{l-m} (\nabla \cdot \boldsymbol{c}_j)^m f_j^{eq} +
\delta_t B_{lm} \partial_{t}^{l-m} (\nabla \cdot \boldsymbol{c}_j)^m F_j
\right] + O((\delta_{t})^{N+1})
\tag{3}
f j ( x , t ) = l = 0 ∑ N ( δ t ) l m = 0 ∑ l [ A l m ∂ t l − m ( ∇ ⋅ c j ) m f j e q + δ t B l m ∂ t l − m ( ∇ ⋅ c j ) m F j ] + O ( ( δ t ) N + 1 ) ( 3 )
其中 f j e q = f j e q ( x , t ) f_j^{eq} = f_j^{eq}(\boldsymbol{x},t) f j e q = f j e q ( x , t ) F j = F j ( x , t ) F_j = F_j(\boldsymbol{x},t) F j = F j ( x , t )
取 N = 2 N=2 N = 2 F j F_j F j l = 2 l=2 l = 2 O ( δ t 3 ) O(\delta_t^3) O ( δ t 3 )
f j ( x , t ) = A 00 f j e q + δ t B 00 F j ⏟ l = 0 + δ t ( A 10 ∂ t + A 11 ∇ ⋅ c j ) f j e q + δ t 2 ( B 10 ∂ t + B 11 ∇ ⋅ c j ) F j ⏟ l = 1 + δ t 2 [ A 20 ∂ t 2 + A 21 ∂ t ( ∇ ⋅ c j ) + A 22 ( ∇ ⋅ c j ) 2 ] f j e q ⏟ l = 2 ( f e q part ) + O ( δ t 3 ) \begin{aligned}
f_j (\boldsymbol{x},t) =& \underbrace{A_{00} f_j^{eq} + \delta_t B_{00} F_j}_{l=0}
\\& + \underbrace{\delta_t (A_{10} \partial_t + A_{11} \nabla \cdot \boldsymbol{c}_j) f_j^{eq} + \delta_t^2 (B_{10} \partial_t + B_{11} \nabla \cdot \boldsymbol{c}_j) F_j}_{l=1}
\\& +\underbrace{\delta_t^2 [A_{20} \partial_t^2 + A_{21} \partial_t \, (\nabla \cdot \boldsymbol{c}_j) + A_{22} (\nabla \cdot \boldsymbol{c}_j)^2] f_j^{eq} }_{l=2 \,(f^{eq} \text{ part})}
\\& + O(\delta_t^3)
\end{aligned}
f j ( x , t ) = l = 0 A 0 0 f j e q + δ t B 0 0 F j + l = 1 δ t ( A 1 0 ∂ t + A 1 1 ∇ ⋅ c j ) f j e q + δ t 2 ( B 1 0 ∂ t + B 1 1 ∇ ⋅ c j ) F j + l = 2 ( f e q part ) δ t 2 [ A 2 0 ∂ t 2 + A 2 1 ∂ t ( ∇ ⋅ c j ) + A 2 2 ( ∇ ⋅ c j ) 2 ] f j e q + O ( δ t 3 )
这里我们写列出 k 阶速度矩的表达式：
M 0 = ∑ j f j e q = ∑ j f j , M 0 F = ∑ j F j , M 1 = ∑ j c j f j e q , M 1 F = ∑ j c j F j , M k = ∑ j c j … c j ⏟ k f j e q , M k F = ∑ j c j … c j ⏟ k F j , ( k > 1 ) \begin{aligned}
M_0 &= \sum_j f_j^{eq} = \sum_j f_j ,& M_{0F} &= \sum_j F_j, \\
M_1 &= \sum_j \boldsymbol{c}_j f_j^{eq},& M_{1F} &= \sum_j \boldsymbol{c}_j F_j, \\
M_k &= \sum_j \underbrace{\boldsymbol{c}_j \dots \boldsymbol{c}_j}_{k} f_j^{eq},& M_{kF} &= \sum_j \underbrace{\boldsymbol{c}_j \dots \boldsymbol{c}_j}_{k} F_j, (k>1)\\
\end{aligned}
M 0 M 1 M k = j ∑ f j e q = j ∑ f j , = j ∑ c j f j e q , = j ∑ k c j … c j f j e q , M 0 F M 1 F M k F = j ∑ F j , = j ∑ c j F j , = j ∑ k c j … c j F j , ( k > 1 )
按理来说, 只要把这些速度矩也匹配上就能去还原宏观方程.
毕竟我们知道, 对于常用的正方形+对称 DnQb 模型, 如果是还原Navier-Stokes方程, 就有:
M 0 = ρ M_0 = \rho M 0 = ρ M 1 = ρ u M_1 = \rho \boldsymbol{u} M 1 = ρ u M 2 = ρ u u + p I M_2 = \rho \boldsymbol{u} \boldsymbol{u} + p \boldsymbol{I} M 2 = ρ u u + p I
不过, 先别急.
让我们注意一下式(3)中的 ( ∇ ⋅ c j ) m f j e q (\nabla \cdot \boldsymbol{c}_j)^m f_j^{eq} ( ∇ ⋅ c j ) m f j e q j j j
∑ j ( ∇ ⋅ c j ) m f j e q = ( ∇ … ∇ ) ⊗ ∑ j ( c j … c j ⏟ k f j e q ) = ∇ m ⊗ M m \begin{aligned}
\sum_{j} (\nabla \cdot \boldsymbol{c}_j)^m f_j^{eq} &= (\nabla \dots \nabla) \otimes \sum_{j} (\underbrace{\boldsymbol{c}_j \dots \boldsymbol{c}_j}_{k} f_j^{eq}) \\
&= \nabla^{m} \otimes M_{m}
\end{aligned}
j ∑ ( ∇ ⋅ c j ) m f j e q = ( ∇ … ∇ ) ⊗ j ∑ ( k c j … c j f j e q ) = ∇ m ⊗ M m
也就是说, M 0 M_0 M 0 M m M_{m} M m m = [ 0 , 1 , . . , l ] m=[0,1,..,l] m = [ 0 , 1 , . . , l ] l = [ 0 , 1 , . . . N ] l=[0,1,...N] l = [ 0 , 1 , . . . N ] N ≥ 1 N \geq 1 N ≥ 1 M k M_k M k
将式(3)代入 M 0 M_0 M 0
M 0 = A 00 M 0 + ∑ l = 1 N δ t l ∑ m = 0 l [ A l m ∂ t l − m ∇ m ⊗ M m ] + ∑ l = 0 N − 1 δ t l + 1 ∑ m = 0 l [ B l m ∂ t l − m ∇ m ⊗ M m F ] + O ( δ t l + 1 ) \begin{aligned}
M_0 =& A_{00} M_0 + \sum_{l=1}^{N} \delta_{t}^l \sum_{m=0}^{l} [A_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{m}] \\
& + \sum_{l=0}^{N-1} \delta_{t}^{l+1} \sum_{m=0}^{l} [B_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{mF}] + O(\delta_{t}^{l+1})
\end{aligned}
M 0 = A 0 0 M 0 + l = 1 ∑ N δ t l m = 0 ∑ l [ A l m ∂ t l − m ∇ m ⊗ M m ] + l = 0 ∑ N − 1 δ t l + 1 m = 0 ∑ l [ B l m ∂ t l − m ∇ m ⊗ M m F ] + O ( δ t l + 1 )
令 A 00 = 1 A_{00}=1 A 0 0 = 1 M 0 M_0 M 0
∑ l = 1 N δ t l − 1 ∑ m = 0 l [ A l m ∂ t l − m ∇ m ⊗ M m ] + ∑ l = 0 N − 1 δ t l ∑ m = 0 l [ B l m ∂ t l − m ∇ m ⊗ M m F ] = O ( δ t l ) (4) \sum_{l=1}^{N} \delta_{t}^{l-1} \sum_{m=0}^{l} [A_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{m}] + \sum_{l=0}^{N-1} \delta_{t}^l \sum_{m=0}^{l} [B_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{mF}] = O(\delta_{t}^l)
\tag{4}
l = 1 ∑ N δ t l − 1 m = 0 ∑ l [ A l m ∂ t l − m ∇ m ⊗ M m ] + l = 0 ∑ N − 1 δ t l m = 0 ∑ l [ B l m ∂ t l − m ∇ m ⊗ M m F ] = O ( δ t l ) ( 4 )
同理, 将式(3)代入 M 1 M_1 M 1
M 1 = A ˉ 00 M 1 + ∑ l = 1 N δ t l ∑ m = 0 l [ A ˉ l m ∂ t l − m ∇ m ⊗ M m ] + ∑ l = 0 N − 1 δ t l + 1 ∑ m = 0 l [ B ˉ l m ∂ t l − m ∇ m ⊗ M m F ] + O ( δ t l + 1 ) \begin{aligned}
M_1 =& \bar{A}_{00} M_1 + \sum_{l=1}^{N} \delta_{t}^l \sum_{m=0}^{l} [\bar{A}_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{m}] \\
& + \sum_{l=0}^{N-1} \delta_{t}^{l+1} \sum_{m=0}^{l} [\bar{B}_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{mF}] + O(\delta_{t}^{l+1})
\end{aligned}
M 1 = A ˉ 0 0 M 1 + l = 1 ∑ N δ t l m = 0 ∑ l [ A ˉ l m ∂ t l − m ∇ m ⊗ M m ] + l = 0 ∑ N − 1 δ t l + 1 m = 0 ∑ l [ B ˉ l m ∂ t l − m ∇ m ⊗ M m F ] + O ( δ t l + 1 )
其中 A ˉ l m , B ˉ l m \bar{A}_{lm}, \bar{B}_{lm} A ˉ l m , B ˉ l m
A ˉ l m = C l m l ! ∑ k , q a ˉ k q q l − m k m , B ˉ l m = C l m l ! ∑ k , q b ˉ k q q l − m k m \bar{A}_{lm} = \frac{C_l^m}{l!} \sum_{k,q} \bar{a}_{kq} q^{l-m} k^{m} ,\quad
\bar{B}_{lm} = \frac{C_l^m}{l!} \sum_{k,q} \bar{b}_{kq} q^{l-m} k^{m}
A ˉ l m = l ! C l m k , q ∑ a ˉ k q q l − m k m , B ˉ l m = l ! C l m k , q ∑ b ˉ k q q l − m k m
令 A ˉ 00 = 1 \bar{A}_{00}=1 A ˉ 0 0 = 1 M 1 M_1 M 1
∑ l = 1 N δ t l − 1 ∑ m = 0 l [ A ˉ l m ∂ t l − m ∇ m ⊗ M m + 1 ] + ∑ l = 0 N − 1 δ t l ∑ m = 0 l [ B ˉ l m ∂ t l − m ∇ m ⊗ M ( m + 1 ) F ] = O ( δ t l ) (5) \sum_{l=1}^{N} \delta_{t}^{l-1} \sum_{m=0}^{l} [\bar{A}_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{m+1}] + \sum_{l=0}^{N-1} \delta_{t}^l \sum_{m=0}^{l} [\bar{B}_{lm} \, \partial_{t}^{l-m} \, \nabla^{m} \otimes M_{(m+1)F}] = O(\delta_{t}^l)
\tag{5}
l = 1 ∑ N δ t l − 1 m = 0 ∑ l [ A ˉ l m ∂ t l − m ∇ m ⊗ M m + 1 ] + l = 0 ∑ N − 1 δ t l m = 0 ∑ l [ B ˉ l m ∂ t l − m ∇ m ⊗ M ( m + 1 ) F ] = O ( δ t l ) ( 5 )
【1】对于式(4), 当 N = 1 N=1 N = 1
A 10 ∂ t M 0 + A 11 ∇ ⋅ M 1 = − B 00 M 0 F + O ( δ t ) (6a) A_{10} \partial_{t} M_0 + A_{11} \nabla \cdot M_1 = - B_{00} M_{0F} + O(\delta_t)
\tag{6a}
A 1 0 ∂ t M 0 + A 1 1 ∇ ⋅ M 1 = − B 0 0 M 0 F + O ( δ t ) ( 6 a )
当 N = 2 N=2 N = 2
( A 10 ∂ t M 0 + A 11 ∇ ⋅ M 1 ) + δ t [ A 20 ∂ t 2 M 0 + A 21 ∂ t ∇ ⋅ M 1 + A 22 ∇ ∇ : M 2 ] = − B 00 M 0 F − δ t ( B 10 ∂ t M 0 F + B 11 ∇ ⋅ M 1 F ) + O ( δ t 2 ) (6b) \begin{aligned}
( A_{10} \partial_{t} M_0 + A_{11} \nabla \cdot M_1 ) &+ \delta_t [
A_{20} \partial_{t}^2 M_0 + A_{21} \partial_{t} \nabla \cdot M_1 + A_{22} \nabla\nabla : M_2
] \\
& = -B_{00} M_{0F} - \delta_t (B_{10} \partial_{t} M_{0F} + B_{11} \nabla \cdot M_{1F}) + O(\delta_t^2)
\end{aligned}
\tag{6b}
( A 1 0 ∂ t M 0 + A 1 1 ∇ ⋅ M 1 ) + δ t [ A 2 0 ∂ t 2 M 0 + A 2 1 ∂ t ∇ ⋅ M 1 + A 2 2 ∇ ∇ : M 2 ] = − B 0 0 M 0 F − δ t ( B 1 0 ∂ t M 0 F + B 1 1 ∇ ⋅ M 1 F ) + O ( δ t 2 ) ( 6 b )
【2】对于式(5), 当 N = 1 N=1 N = 1
A ˉ 10 ∂ t M 1 + A ˉ 11 ∇ ⋅ M 2 = − B ˉ 00 M 1 F + O ( δ t ) (7a) \bar{A}_{10} \partial_{t} M_1 + \bar{A}_{11} \nabla \cdot M_2 = -\bar{B}_{00} M_{1F} + O(\delta_t)
\tag{7a}
A ˉ 1 0 ∂ t M 1 + A ˉ 1 1 ∇ ⋅ M 2 = − B ˉ 0 0 M 1 F + O ( δ t ) ( 7 a )
当 N = 2 N=2 N = 2
( A ˉ 10 ∂ t M 1 + A ˉ 11 ∇ ⋅ M 2 ) + δ t [ A ˉ 20 ∂ t 2 M 1 + A ˉ 21 ∂ t ∇ ⋅ M 2 + A ˉ 22 ∇ ∇ : M 3 ] = − B ˉ 00 M 1 F − δ t ( B ˉ 10 ∂ t M 1 F + B ˉ 11 ∇ ⋅ M 2 F ) + O ( δ t 2 ) (7b) \begin{aligned}
( \bar{A}_{10} \partial_{t} M_1 + \bar{A}_{11} \nabla \cdot M_2 ) &+ \delta_t [
\bar{A}_{20} \partial_{t}^2 M_1 + \bar{A}_{21} \partial_{t} \nabla \cdot M_2 + \bar{A}_{22} \nabla\nabla : M_3
] \\
& = -\bar{B}_{00} M_{1F} - \delta_t (\bar{B}_{10} \partial_{t} M_{1F} + \bar{B}_{11} \nabla \cdot M_{2F}) + O(\delta_t^2)
\end{aligned}
\tag{7b}
( A ˉ 1 0 ∂ t M 1 + A ˉ 1 1 ∇ ⋅ M 2 ) + δ t [ A ˉ 2 0 ∂ t 2 M 1 + A ˉ 2 1 ∂ t ∇ ⋅ M 2 + A ˉ 2 2 ∇ ∇ : M 3 ] = − B ˉ 0 0 M 1 F − δ t ( B ˉ 1 0 ∂ t M 1 F + B ˉ 1 1 ∇ ⋅ M 2 F ) + O ( δ t 2 ) ( 7 b )
Shi et al. 做出下面代换来简化式(6-7), 这些代换在下文会用到:
A 10 = A 11 = − B 00 A_{10} = A_{11} = - B_{00} A 1 0 = A 1 1 = − B 0 0 A ˉ 10 = A ˉ 11 = − B ˉ 00 \bar{A}_{10} = \bar{A}_{11} = - \bar{B}_{00} A ˉ 1 0 = A ˉ 1 1 = − B ˉ 0 0 S = − B 00 A 10 M 0 F = M 0 F \displaystyle S = -\frac{B_{00}}{A_{10}} M_{0F} = M_{0F} S = − A 1 0 B 0 0 M 0 F = M 0 F F = − B ˉ 00 A ˉ 10 M 1 F = M 1 F \displaystyle \boldsymbol{F} = -\frac{\bar{B}_{00}}{\bar{A}_{10}} M_{1F} = M_{1F} F = − A ˉ 1 0 B ˉ 0 0 M 1 F = M 1 F
这里我们先把非线性对流扩散方程写出来, 作为参照:
∂ t ϕ + ∇ ⋅ B = S + α ∇ ∇ : D + O ( δ t 2 ) = S + ∇ ⋅ [ α ∇ ⋅ D ] + O ( δ t 2 ) (8) \begin{aligned}
\partial_{t} \phi + \nabla \cdot \bold{B} &= \boldsymbol{S} + \alpha \nabla\nabla:\bold{D} + O(\delta_t^2)\\
&= \boldsymbol{S} + \nabla \cdot [\alpha \nabla \cdot \bold{D}] + O(\delta_t^2)
\end{aligned}
\tag{8}
∂ t ϕ + ∇ ⋅ B = S + α ∇ ∇ : D + O ( δ t 2 ) = S + ∇ ⋅ [ α ∇ ⋅ D ] + O ( δ t 2 ) ( 8 )
若我们只考虑 “各向同性”+“不可压缩流” 来降低心智负担, 则
ϕ \phi ϕ u \boldsymbol{u} u ϕ \phi ϕ α \alpha α B = ϕ u \bold{B}=\phi \, \boldsymbol{u} B = ϕ u ∇ ⋅ u = 0 \nabla\cdot\boldsymbol{u}=0 ∇ ⋅ u = 0 ∇ ⋅ B = u ⋅ ∇ ϕ \nabla\cdot\bold{B}=\boldsymbol{u}\cdot\nabla\phi ∇ ⋅ B = u ⋅ ∇ ϕ D = ϕ [ I ] \bold{D}=\phi \, [\bold{I}] D = ϕ [ I ] S \boldsymbol{S} S
式(6)可化简为
∂ t M 0 + ∇ ⋅ M 1 = S + O ( δ t ) (9a) \partial_{t} M_0 + \nabla \cdot M_1 = S + O(\delta_t)
\tag{9a}
∂ t M 0 + ∇ ⋅ M 1 = S + O ( δ t ) ( 9 a )
( ∂ t M 0 + ∇ ⋅ M 1 ) + δ t A 10 [ A 20 ∂ t 2 M 0 + A 21 ∂ t ∇ ⋅ M 1 + A 22 ∇ ∇ : M 2 ] = S − δ t A 10 ( B 10 ∂ t M 0 F + B 11 ∇ ⋅ M 1 F ) + O ( δ t 2 ) (9b) \begin{aligned}
( \partial_{t} M_0 + \nabla \cdot M_1 ) &+ \frac{\delta_t}{A_{10}} [
A_{20} \partial_{t}^2 M_0 + A_{21} \partial_{t} \nabla \cdot M_1 + A_{22} \nabla\nabla : M_2
] \\
& = S - \frac{\delta_t}{A_{10}} (B_{10} \partial_{t} M_{0F} + B_{11} \nabla \cdot M_{1F}) + O(\delta_t^2)
\end{aligned}
\tag{9b}
( ∂ t M 0 + ∇ ⋅ M 1 ) + A 1 0 δ t [ A 2 0 ∂ t 2 M 0 + A 2 1 ∂ t ∇ ⋅ M 1 + A 2 2 ∇ ∇ : M 2 ] = S − A 1 0 δ t ( B 1 0 ∂ t M 0 F + B 1 1 ∇ ⋅ M 1 F ) + O ( δ t 2 ) ( 9 b )
而 ∂ t 2 M 0 \partial_{t}^2 M_0 ∂ t 2 M 0
∂ t 2 M 0 = − ∂ t ∇ ⋅ M 1 + ∂ t S + O ( δ t ) \partial_{t}^2 M_0 = - \partial_{t} \nabla \cdot M_1 + \partial_{t} S + O(\delta_t)
∂ t 2 M 0 = − ∂ t ∇ ⋅ M 1 + ∂ t S + O ( δ t )
则式(9b)可以进一步化简为
( ∂ t M 0 + ∇ ⋅ M 1 ) + δ t A 10 ( B 10 + A 20 ) ∂ t S = S − δ t A 10 [ ( A 21 − A 20 ) ∂ t ∇ ⋅ M 1 + B 11 ∇ ⋅ M 1 F + A 22 ∇ ∇ : M 2 ] + O ( δ t 2 ) (9c) \begin{aligned}
( \partial_{t} M_0 + \nabla \cdot M_1 ) &+ \frac{\delta_t}{A_{10}} (B_{10} + A_{20}) \partial_{t}S \\
& = S - \frac{\delta_t}{A_{10}} [(A_{21} - A_{20}) \partial_{t} \nabla \cdot M_1 \\
&\quad + B_{11} \nabla \cdot M_{1F} + A_{22} \nabla\nabla : M_2] \\
&\quad + O(\delta_t^2)
\end{aligned}
\tag{9c}
( ∂ t M 0 + ∇ ⋅ M 1 ) + A 1 0 δ t ( B 1 0 + A 2 0 ) ∂ t S = S − A 1 0 δ t [ ( A 2 1 − A 2 0 ) ∂ t ∇ ⋅ M 1 + B 1 1 ∇ ⋅ M 1 F + A 2 2 ∇ ∇ : M 2 ] + O ( δ t 2 ) ( 9 c )
令 ( B 10 + A 20 ) ∂ t S = 0 (B_{10} + A_{20}) \partial_{t}S = 0 ( B 1 0 + A 2 0 ) ∂ t S = 0
M 0 = ϕ , M 1 = B , M 2 = β c s 2 D + C , α = − A 22 A 10 δ t β c s 2 , M 1 F = − [ ( A 21 − A 20 ) ∂ t B + A 22 ∇ ⋅ C ] \begin{aligned}
M_0 = \phi, \,
M_1 &= \bold{B}, \,
M_2 = \beta c_s^2 \bold{D} + \bold{C}, \,
\alpha = -\frac{A_{22}}{A_{10}} \delta_t \beta c_s^2 ,\\
M_{1F} &= -[(A_{21} - A_{20}) \partial_{t}\bold{B} + A_{22} \nabla\cdot\bold{C}]
\end{aligned}
M 0 = ϕ , M 1 M 1 F = B , M 2 = β c s 2 D + C , α = − A 1 0 A 2 2 δ t β c s 2 , = − [ ( A 2 1 − A 2 0 ) ∂ t B + A 2 2 ∇ ⋅ C ]
这里和 C \bold{C} C A 22 ∇ ∇ : M 2 A_{22} \nabla\nabla : M_2 A 2 2 ∇ ∇ : M 2
首先, 式(6)的二阶结果即为式(9c). 这里只需要结合式(7)的二阶结果, 再和Navier-Stokes方程组进行参数匹配即可. 所以下面主要讨论式(7).
式(7)可化简为
∂ t M 1 + ∇ ⋅ M 2 = F + O ( δ t ) (10a) \partial_{t} M_1 +\nabla \cdot M_2 = \boldsymbol{F} + O(\delta_t)
\tag{10a}
∂ t M 1 + ∇ ⋅ M 2 = F + O ( δ t ) ( 1 0 a )
( ∂ t M 1 + ∇ ⋅ M 2 ) + δ t A ˉ 10 [ A ˉ 20 ∂ t 2 M 1 + A ˉ 21 ∂ t ∇ ⋅ M 2 + A ˉ 22 ∇ ∇ : M 3 ] = F − δ t A ˉ 10 ( B ˉ 10 ∂ t F + B ˉ 11 ∇ ⋅ M 2 F ) + O ( δ t 2 ) (10b) \begin{aligned}
( \partial_{t} M_1 + \nabla \cdot M_2 ) &+ \frac{\delta_t}{\bar{A}_{10}} [
\bar{A}_{20} \partial_{t}^2 M_1 + \bar{A}_{21} \partial_{t} \nabla \cdot M_2 + \bar{A}_{22} \nabla\nabla : M_3
] \\
& = \boldsymbol{F} - \frac{\delta_t}{\bar{A}_{10}} (\bar{B}_{10} \partial_{t} \boldsymbol{F} + \bar{B}_{11} \nabla \cdot M_{2F}) + O(\delta_t^2)
\end{aligned}
\tag{10b}
( ∂ t M 1 + ∇ ⋅ M 2 ) + A ˉ 1 0 δ t [ A ˉ 2 0 ∂ t 2 M 1 + A ˉ 2 1 ∂ t ∇ ⋅ M 2 + A ˉ 2 2 ∇ ∇ : M 3 ] = F − A ˉ 1 0 δ t ( B ˉ 1 0 ∂ t F + B ˉ 1 1 ∇ ⋅ M 2 F ) + O ( δ t 2 ) ( 1 0 b )
同理, 用式(10a)表示 ∂ t 2 M 1 \partial_{t}^2 M_1 ∂ t 2 M 1
∂ t 2 M 1 = − ∂ t ∇ ⋅ M 2 + ∂ t F + O ( δ t ) \partial_{t}^2 M_1 = -\partial_{t} \nabla \cdot M_2 + \partial_{t}\boldsymbol{F} + O(\delta_t)
∂ t 2 M 1 = − ∂ t ∇ ⋅ M 2 + ∂ t F + O ( δ t )
再次代入式(10b)进行化简, 得到:
( ∂ t M 1 + ∇ ⋅ M 2 ) + δ t A ˉ 10 ( B ˉ 10 ∂ t F + A ˉ 20 ∂ t F ) = F − δ t A ˉ 10 [ ( A ˉ 21 − A ˉ 20 ) ∂ t ∇ ⋅ M 2 + B ˉ 11 ∇ ⋅ M 2 F + A ˉ 22 ∇ ∇ : M 3 ] + O ( δ t 2 ) (10c) \begin{aligned}
( \partial_{t} M_1 + \nabla \cdot M_2 ) &+ \frac{\delta_t}{\bar{A}_{10}} (\bar{B}_{10} \partial_{t}\boldsymbol{F} + \bar{A}_{20} \partial_{t}\boldsymbol{F}) \\
& = \boldsymbol{F} - \frac{\delta_t}{\bar{A}_{10}} [(\bar{A}_{21} - \bar{A}_{20}) \partial_{t} \nabla \cdot M_2 \\
&\quad + \bar{B}_{11} \nabla \cdot M_{2F} + \bar{A}_{22} \nabla\nabla : M_3] \\
&\quad + O(\delta_t^2)
\end{aligned}
\tag{10c}
( ∂ t M 1 + ∇ ⋅ M 2 ) + A ˉ 1 0 δ t ( B ˉ 1 0 ∂ t F + A ˉ 2 0 ∂ t F ) = F − A ˉ 1 0 δ t [ ( A ˉ 2 1 − A ˉ 2 0 ) ∂ t ∇ ⋅ M 2 + B ˉ 1 1 ∇ ⋅ M 2 F + A ˉ 2 2 ∇ ∇ : M 3 ] + O ( δ t 2 ) ( 1 0 c )
【1】对于式(6)的二阶结果(即式(9c)), 再次利用式(10a), 可得
∂ t ∇ ⋅ M 1 = − ∇ ∇ : M 2 + ∇ ⋅ F + O ( δ t ) \partial_{t} \nabla \cdot M_1 = -\nabla \nabla : M_2 + \nabla \cdot \boldsymbol{F} + O(\delta_t)
∂ t ∇ ⋅ M 1 = − ∇ ∇ : M 2 + ∇ ⋅ F + O ( δ t )
再令 − B 11 F = ( A 21 − A 20 ) F -B_{11} \boldsymbol{F} = (A_{21} - A_{20}) \boldsymbol{F} − B 1 1 F = ( A 2 1 − A 2 0 ) F ( B 10 + A 20 ) ∂ t S = 0 (B_{10} + A_{20}) \partial_{t}S = 0 ( B 1 0 + A 2 0 ) ∂ t S = 0
∂ t M 0 + ∇ ⋅ M 1 = S − δ t A 10 [ ( A 22 − ( A 21 − A 20 ) ) ∇ ∇ : M 2 ] + O ( δ t 2 ) (11a) \partial_{t} M_0 + \nabla\cdot M_1 = S - \frac{\delta_t}{A_{10}} [(A_{22} - (A_{21} - A_{20}))\nabla \nabla : M_2] + O(\delta_t^2)
\tag{11a}
∂ t M 0 + ∇ ⋅ M 1 = S − A 1 0 δ t [ ( A 2 2 − ( A 2 1 − A 2 0 ) ) ∇ ∇ : M 2 ] + O ( δ t 2 ) ( 1 1 a )
式(11a)将用于匹配连续性方程.
【2】对于式(7)的二阶结果(即式(10c))
令 ( B ˉ 10 ∂ t F + A ˉ 20 ∂ t F ) = 0 (\bar{B}_{10} \partial_{t}\boldsymbol{F} + \bar{A}_{20} \partial_{t}\boldsymbol{F}) = 0 ( B ˉ 1 0 ∂ t F + A ˉ 2 0 ∂ t F ) = 0
( ∂ t M 1 + ∇ ⋅ M 2 ) = F − δ t A ˉ 10 [ ( A ˉ 21 − A ˉ 20 ) ∂ t ∇ ⋅ M 2 + B ˉ 11 ∇ ⋅ M 2 F + A ˉ 22 ∇ ∇ : M 3 ] + O ( δ t 2 ) (11b) \begin{aligned}
( \partial_{t} M_1 + \nabla \cdot M_2 ) & = \boldsymbol{F} - \frac{\delta_t}{\bar{A}_{10}} [(\bar{A}_{21} - \bar{A}_{20}) \partial_{t} \nabla \cdot M_2 \\
&\quad + \bar{B}_{11} \nabla \cdot M_{2F} + \bar{A}_{22} \nabla\nabla : M_3] \\
&\quad + O(\delta_t^2)
\end{aligned}
\tag{11b}
( ∂ t M 1 + ∇ ⋅ M 2 ) = F − A ˉ 1 0 δ t [ ( A ˉ 2 1 − A ˉ 2 0 ) ∂ t ∇ ⋅ M 2 + B ˉ 1 1 ∇ ⋅ M 2 F + A ˉ 2 2 ∇ ∇ : M 3 ] + O ( δ t 2 ) ( 1 1 b )
式(11a)将用于匹配Navier-Stokes方程.
【3】 针对式(11a)和式(11b), 我们需要给出如下参数匹配
M 0 = ρ , M 1 = ρ u , M 2 = ρ u u + p I , M 3 = c s 2 ρ u I + M 30 , M 2 F = − 1 B ˉ 11 [ ( A ˉ 21 − A ˉ 20 ) ∂ t M 2 + A ˉ 22 ( ∇ ⋅ M 30 + ( ∇ ⋅ ρ u I + u ∇ ρ + ( u ∇ ρ ) T ) c s 2 ) ] \begin{aligned}
M_0 &= \rho , \\
M_1 &= \rho \bold{u} , \\
M_2 &= \rho \bold{u} \bold{u} + p \bold{I} , \\
M_3 &= c_s^2 \rho \bold{u} \bold{I} + M_{30} , \\
M_{2F} &= \frac{-1}{\bar{B}_{11}} [(\bar{A}_{21} - \bar{A}_{20}) \partial_{t} M_2 + \bar{A}_{22} ( \nabla \cdot M_{30} +
\\&\quad (\nabla \cdot \rho \bold{u} \bold{I} + \bold{u} \nabla\rho + (\bold{u} \nabla\rho)^{T}) c_s^2 \, )] \\
\end{aligned}
M 0 M 1 M 2 M 3 M 2 F = ρ , = ρ u , = ρ u u + p I , = c s 2 ρ u I + M 3 0 , = B ˉ 1 1 − 1 [ ( A ˉ 2 1 − A ˉ 2 0 ) ∂ t M 2 + A ˉ 2 2 ( ∇ ⋅ M 3 0 + ( ∇ ⋅ ρ u I + u ∇ ρ + ( u ∇ ρ ) T ) c s 2 ) ]
且运动粘度 ν = A ˉ 22 A ˉ 10 δ t c s 2 \nu = \frac{\bar{A}_{22}}{\bar{A}_{10}} \delta_t c_s^2 ν = A ˉ 1 0 A ˉ 2 2 δ t c s 2
∂ t ρ + ∇ ⋅ ( ρ u ) = S − δ t A 10 [ ( A 22 − ( A 21 − A 20 ) ) ∇ ∇ : M 2 ] + O ( δ t 2 ) ∂ t ( ρ u ) + ∇ ⋅ ( ρ u u + p I ) = F + ν ∇ ⋅ ρ [ ∇ u + ( ∇ u ) T ] + O ( δ t 2 ) \begin{aligned}
\partial_t \rho + \nabla \cdot (\rho \boldsymbol{u}) &=
S - \frac{\delta_t}{A_{10}} [(A_{22} - (A_{21} - A_{20}))\nabla \nabla : M_2] + O(\delta_t^2) \\
\partial_t (\rho \boldsymbol{u}) + \nabla \cdot (\rho \boldsymbol{u}\boldsymbol{u} + p \bold{I}) &= \boldsymbol{F} + \nu \nabla \cdot \rho [\nabla\boldsymbol{u} + (\nabla\boldsymbol{u})^T] + O(\delta_t^2) \\
\end{aligned}
∂ t ρ + ∇ ⋅ ( ρ u ) ∂ t ( ρ u ) + ∇ ⋅ ( ρ u u + p I ) = S − A 1 0 δ t [ ( A 2 2 − ( A 2 1 − A 2 0 ) ) ∇ ∇ : M 2 ] + O ( δ t 2 ) = F + ν ∇ ⋅ ρ [ ∇ u + ( ∇ u ) T ] + O ( δ t 2 )
注:
Shi et al. 原文式(61a)写的 M 3 = c s 2 Δ ⋅ ρ u + M 30 M_3 = c_s^2 \Delta \cdot \rho \bold{u} + M_{30} M 3 = c s 2 Δ ⋅ ρ u + M 3 0 Δ \Delta Δ ρ u \rho \boldsymbol{u} ρ u
这里的 u I \bold{u} \bold{I} u I u α δ β γ + u β δ α γ + u γ δ α β u_{\alpha} \delta_{\beta\gamma} + u_{\beta} \delta_{\alpha\gamma} + u_{\gamma} \delta_{\alpha\beta} u α δ β γ + u β δ α γ + u γ δ α β
p = c s 2 ρ p = c_s^2 \rho p = c s 2 ρ