注：该文章同步发在知乎专栏 和CSDN Simplified and Highly Stable Lattice Boltzmann Method (SHSLBM) 下面的内容只是我对相关文献进行简单阅读后的部分整理
针对流场的 LBGK 演化方程为
f α ( r + c α + Δ t , t + Δ t ) − f α ( r , t ) = − 1 τ [ f α ( r , t ) − f α ( e q ) ] (1) f_{\alpha} (\boldsymbol{r} + \boldsymbol{c}_\alpha + \Delta t, t + \Delta t) - f_{\alpha} (\boldsymbol{r}, t) =
-\frac{1}{\tau} \left[ f_{\alpha} (\boldsymbol{r}, t) - f_{\alpha}^{(eq)} \right]
\tag{1}
f α ( r + c α + Δ t , t + Δ t ) − f α ( r , t ) = − τ 1 [ f α ( r , t ) − f α ( e q ) ] ( 1 )
其中 f α ( r , t ) f_{\alpha} (\boldsymbol{r}, t) f α ( r , t ) t t t r \boldsymbol{r} r c α \boldsymbol{c}_\alpha c α τ \tau τ f α ( e q ) f_{\alpha}^{(eq)} f α ( e q ) c α \boldsymbol{c}_\alpha c α
f α ( e q ) = ρ w α [ 1 + c α ⋅ u c s 2 + ( c α ⋅ u ) 2 2 c s 4 − u ⋅ u 2 c s 2 ] (2) f_{\alpha}^{(eq)} = \rho w_\alpha \left[ 1 + \frac{\boldsymbol{c}_\alpha \cdot \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_\alpha \cdot \boldsymbol{u})^2}{2 c_s^4} - \frac{\boldsymbol{u} \cdot \boldsymbol{u}}{2 c_s^2} \right] \tag{2}
f α ( e q ) = ρ w α [ 1 + c s 2 c α ⋅ u + 2 c s 4 ( c α ⋅ u ) 2 − 2 c s 2 u ⋅ u ] ( 2 )
式 (2) 中 ρ \rho ρ u = ( u x , u y , u z ) T \boldsymbol{u} = (u_x, u_y, u_z)^T u = ( u x , u y , u z ) T w α w_\alpha w α c α \boldsymbol{c}_\alpha c α c s c_s c s w α w_\alpha w α c α \boldsymbol{c}_\alpha c α c s c_s c s
对于 LBM 而言，宏观量表示为
ρ = ∑ α f α = ∑ α f α ( e q ) , ρ u = ∑ α c α f α = ∑ α c α f α ( e q ) (3) \rho = \sum_{\alpha} f_{\alpha} = \sum_{\alpha} f_{\alpha}^{(eq)}, \quad \rho\boldsymbol{u} = \sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha} = \sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)}
\tag{3}
ρ = α ∑ f α = α ∑ f α ( e q ) , ρ u = α ∑ c α f α = α ∑ c α f α ( e q ) ( 3 )
Boltzmann 方程和 LBGK 方程都可以通过多尺度展开还原为 Navier-Stokes 方程。这里以 Chapman-Enskog 为例：
f α = f α ( e q ) + ϵ f α ( 1 ) + ϵ 2 f α ( 2 ) + . . . , ∂ ∂ t = ϵ ∂ ∂ t 0 + ϵ 2 ∂ ∂ t 1 , ∂ ∂ r = ϵ ∂ ∂ r 0 f_{\alpha} = f_{\alpha}^{(eq)} + \epsilon f_{\alpha}^{(1)} + \epsilon^2 f_{\alpha}^{(2)} + ... , \quad \frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_0} + \epsilon^2 \frac{\partial}{\partial t_1}, \quad \frac{\partial}{\partial \boldsymbol{r}} = \epsilon \frac{\partial}{\partial \boldsymbol{r}_0}
f α = f α ( e q ) + ϵ f α ( 1 ) + ϵ 2 f α ( 2 ) + . . . , ∂ t ∂ = ϵ ∂ t 0 ∂ + ϵ 2 ∂ t 1 ∂ , ∂ r ∂ = ϵ ∂ r 0 ∂
上述展开保留至 ϵ 1 \epsilon^{1} ϵ 1
∂ ρ ∂ t + ∇ ⋅ ( ∑ α c α f α ( e q ) ) = 0 (4) \frac{\partial \rho}{\partial t} + \nabla \cdot \left(\sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \right) = 0
\tag{4}
∂ t ∂ ρ + ∇ ⋅ ( α ∑ c α f α ( e q ) ) = 0 ( 4 )
和
∂ ρ u ∂ t + ∇ ⋅ Π = 0 (5) \frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \Pi = \boldsymbol{0} \tag{5}
∂ t ∂ ρ u + ∇ ⋅ Π = 0 ( 5 )
式 (5) 中张量 Π \Pi Π
Π β γ = ∑ α ( c α ) β ( c α ) γ [ f α ( e q ) + ( 1 − 1 2 τ ) ϵ f α ( 1 ) ] \Pi_{\beta \gamma} = \sum_{\alpha} (\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} \left[ f_{\alpha}^{(eq)} + \left(1 - \frac{1}{2 \tau}\right) \epsilon f_{\alpha}^{(1)} \right]
Π β γ = α ∑ ( c α ) β ( c α ) γ [ f α ( e q ) + ( 1 − 2 τ 1 ) ϵ f α ( 1 ) ]
由于只保留到 ϵ 1 \epsilon^{1} ϵ 1 ϵ f α ( 1 ) ≈ f α − f α ( e q ) = f α ( n e q ) \epsilon f_{\alpha}^{(1)} \approx f_{\alpha} - f_{\alpha}^{(eq)} = f_{\alpha}^{(neq)} ϵ f α ( 1 ) ≈ f α − f α ( e q ) = f α ( n e q )
Π β γ = ∑ α ( c α ) β ( c α ) γ [ f α ( e q ) + ( 1 − 1 2 τ ) f α ( n e q ) ] (6) \Pi_{\beta \gamma} = \sum_{\alpha} (\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} \left[ f_{\alpha}^{(eq)} + \left(1 - \frac{1}{2 \tau}\right) f_{\alpha}^{(neq)} \right]
\tag{6}
Π β γ = α ∑ ( c α ) β ( c α ) γ [ f α ( e q ) + ( 1 − 2 τ 1 ) f α ( n e q ) ] ( 6 )
并且存在如下关系
f α ( n e q ) ≈ ϵ f α ( 1 ) = − τ [ ∂ ∂ t + c α ⋅ ∇ ] f α ( e q ) Δ t (7) f_{\alpha}^{(neq)} \approx \epsilon f_{\alpha}^{(1)} = -\tau \left[ \frac{\partial}{\partial t} + \boldsymbol{c}_{\alpha} \cdot \nabla \right] f_{\alpha}^{(eq)} \Delta t
\tag{7}
f α ( n e q ) ≈ ϵ f α ( 1 ) = − τ [ ∂ t ∂ + c α ⋅ ∇ ] f α ( e q ) Δ t ( 7 )
令 D = ∂ ∂ t + c α ⋅ ∇ \mathrm{D} = \frac{\partial}{\partial t} + \boldsymbol{c}_{\alpha} \cdot \nabla D = ∂ t ∂ + c α ⋅ ∇ f α ( n e q ) ≈ − τ ⋅ ( Δ t ) ⋅ D f α ( e q ) f_{\alpha}^{(neq)} \approx -\tau \cdot (\Delta t) \cdot \mathrm{D}f_{\alpha}^{(eq)} f α ( n e q ) ≈ − τ ⋅ ( Δ t ) ⋅ D f α ( e q )
流体运动粘度 ν \nu ν τ \tau τ
ν = ( τ − 1 2 ) c s 2 Δ t \nu = \left( \tau - \frac{1}{2} \right) c_s^2 \Delta t
ν = ( τ − 2 1 ) c s 2 Δ t
式 (7) 说明：在保留至 ϵ 1 \epsilon^{1} ϵ 1 f α ( n e q ) f_{\alpha}^{(neq)} f α ( n e q ) f α ( e q ) f_{\alpha}^{(eq)} f α ( e q ) − 1 τ \frac{-1}{\tau} τ − 1 f α ( n e q ) f_{\alpha}^{(neq)} f α ( n e q ) f α ( e q ) f_{\alpha}^{(eq)} f α ( e q )
SHSLBM 的单个时间步迭代分为预报和校正两部分，具体写作：
预报步：计算 ( r , t ) (\boldsymbol{r},t) ( r , t ) ρ ∗ \rho^{*} ρ ∗ u ∗ \boldsymbol{u}^{*} u ∗
ρ ∗ = ∑ α f α ( e q ) ( r − c α Δ t , t − Δ t ) , ρ ∗ u ∗ = ∑ α c α f α ( e q ) ( r − c α Δ t , t − Δ t ) (8) \rho^{*} = \sum_{\alpha} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t), \quad
\rho^{*} \boldsymbol{u}^{*} = \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t)
\tag{8}
ρ ∗ = α ∑ f α ( e q ) ( r − c α Δ t , t − Δ t ) , ρ ∗ u ∗ = α ∑ c α f α ( e q ) ( r − c α Δ t , t − Δ t ) ( 8 )
校正步：根据 ρ ∗ \rho^{*} ρ ∗ u ∗ \boldsymbol{u}^{*} u ∗
{ ρ = ρ ∗ ρ u = ρ ∗ u ∗ − ( 1 − 1 τ ) ∑ α c α [ f α ( n e q ) ( r + Δ t 2 c α , t − Δ t 2 ) − f α ( n e q ) ( r − Δ t 2 c α , t − Δ t 2 ) ] (9-1) \begin{cases}
\rho = \rho^{*} \\
\rho \boldsymbol{u} = \rho^{*} \boldsymbol{u}^{*} - \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} [f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) - f_{\alpha}^{(neq)}(\boldsymbol{r}-\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2})]
\end{cases}
\tag{9-1}
⎩ ⎪ ⎨ ⎪ ⎧ ρ = ρ ∗ ρ u = ρ ∗ u ∗ − ( 1 − τ 1 ) α ∑ c α [ f α ( n e q ) ( r + 2 Δ t c α , t − 2 Δ t ) − f α ( n e q ) ( r − 2 Δ t c α , t − 2 Δ t ) ] ( 9 - 1 )
式 (9-1) 中非平衡态的计算通过对式 (7) 的差分实现，即：
f α ( n e q ) ( r − Δ t 2 c α , t − Δ t 2 ) = − τ [ f α ( e q ) ( r , t ) − f α ( e q ) ( r − c α Δ t , t − Δ t ) ] + O ( ( Δ t ) 3 ) f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) = -\tau \left[ f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - f_{\alpha}^{(eq)}(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t - \Delta t) \right] + O((\Delta t)^3)
f α ( n e q ) ( r − 2 Δ t c α , t − 2 Δ t ) = − τ [ f α ( e q ) ( r , t ) − f α ( e q ) ( r − c α Δ t , t − Δ t ) ] + O ( ( Δ t ) 3 )
同理：
f α ( n e q ) ( r + Δ t 2 c α , t − Δ t 2 ) = − τ [ f α ( e q ) ( r + c α Δ t , t ) − f α ( e q ) ( r , t − Δ t ) ] + O ( ( Δ t ) 3 ) f_{\alpha}^{(neq)}(\boldsymbol{r} + \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) = -\tau \left[ f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) - f_{\alpha}^{(eq)}(\boldsymbol{r}, t - \Delta t) \right] + O((\Delta t)^3)
f α ( n e q ) ( r + 2 Δ t c α , t − 2 Δ t ) = − τ [ f α ( e q ) ( r + c α Δ t , t ) − f α ( e q ) ( r , t − Δ t ) ] + O ( ( Δ t ) 3 )
在这两条式子中， f α ( e q ) ( r , t ) f_{\alpha}^{(eq)}(\boldsymbol{r}, t) f α ( e q ) ( r , t ) f α ( e q ) ( r + c α Δ t , t ) f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) f α ( e q ) ( r + c α Δ t , t )
f α ( e q ) ( r , t ) = f α ( e q ) ( ρ ∗ ( r , t ) , u ∗ ( r , t ) ) f_{\alpha}^{(eq)}(\boldsymbol{r}, t) = f_{\alpha}^{(eq)}(\rho^{*}(\boldsymbol{r}, t), \boldsymbol{u}^{*}(\boldsymbol{r}, t))
f α ( e q ) ( r , t ) = f α ( e q ) ( ρ ∗ ( r , t ) , u ∗ ( r , t ) )
式 (9-1) 中 ρ u \rho \boldsymbol{u} ρ u f α ( n e q ) ( r − Δ t 2 c α , t − Δ t 2 ) f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) f α ( n e q ) ( r − 2 Δ t c α , t − 2 Δ t )
∑ α c α f α ( n e q ) ( r − Δ t 2 c α , t − Δ t 2 ) = − τ [ ∑ α c α f α ( e q ) ( r , t ) − ∑ α c α f α ( e q ) ( r − c α Δ t , t − Δ t ) ] + O ( ( Δ t ) 3 ) = − τ ( ρ ∗ u ∗ − ρ ∗ u ∗ ) + O ( ( Δ t ) 3 ) = O ( ( Δ t ) 3 ) \begin{aligned}
\sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) &= -\tau \left[ \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \right.\\
&\quad\left. -\sum\limits_{\alpha} \boldsymbol{c}_{\alpha}f_{\alpha}^{(eq)}(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t - \Delta t)\right] + O((\Delta t)^3) \\
&= -\tau (\rho^{*} \boldsymbol{u}^{*} - \rho^{*} \boldsymbol{u}^{*}) + O((\Delta t)^3) \\
&= O((\Delta t)^3)
\end{aligned}
α ∑ c α f α ( n e q ) ( r − 2 Δ t c α , t − 2 Δ t ) = − τ [ α ∑ c α f α ( e q ) ( r , t ) − α ∑ c α f α ( e q ) ( r − c α Δ t , t − Δ t ) ] + O ( ( Δ t ) 3 ) = − τ ( ρ ∗ u ∗ − ρ ∗ u ∗ ) + O ( ( Δ t ) 3 ) = O ( ( Δ t ) 3 )
所以式 (9-1) 可进一步化简为式 (9-2)
{ ρ = ρ ∗ ρ u = ρ ∗ u ∗ − ( 1 − 1 τ ) ∑ α c α f α ( n e q ) ( r + Δ t 2 c α , t − Δ t 2 ) (9-2) \begin{cases}
\rho = \rho^{*} \\
\rho \boldsymbol{u} = \rho^{*} \boldsymbol{u}^{*} - \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2})
\end{cases}
\tag{9-2}
⎩ ⎪ ⎨ ⎪ ⎧ ρ = ρ ∗ ρ u = ρ ∗ u ∗ − ( 1 − τ 1 ) α ∑ c α f α ( n e q ) ( r + 2 Δ t c α , t − 2 Δ t ) ( 9 - 2 )
事实上，式 (8) 和式 (9-1), (9-2) 的预报和校正存在着对应的宏观偏微分方程。下面简要说明一下。
式 (8) 对应的宏观方程为
{ ∂ ρ ∂ t + ∇ ⋅ ( ∑ α c α f α ( e q ) ) = 0 ∂ ρ u ∂ t + ∇ ⋅ [ ∑ α ( c α ) β ( c α ) γ f α ( e q ) + 1 2 τ ∑ α ( c α ) β ( c α ) γ f α ( n e q ) ] = 0 (10) \begin{cases}
\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\sum\limits_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \right) = 0 \\
\frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \left[ \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(eq)} + \frac{1}{2 \tau} \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)} \right] = \boldsymbol{0}
\end{cases}
\tag{10}
⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ ∂ t ∂ ρ + ∇ ⋅ ( α ∑ c α f α ( e q ) ) = 0 ∂ t ∂ ρ u + ∇ ⋅ [ α ∑ ( c α ) β ( c α ) γ f α ( e q ) + 2 τ 1 α ∑ ( c α ) β ( c α ) γ f α ( n e q ) ] = 0 ( 1 0 )
首先，对 f α ( e q ) ( r − c α Δ t , t − Δ t ) f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t) f α ( e q ) ( r − c α Δ t , t − Δ t ) ( r , t ) (\boldsymbol{r},t) ( r , t )
f α ( e q ) ( r − c α Δ t , t − Δ t ) = f α ( e q ) ( r , t ) − ( Δ t ) D f α ( e q ) ( r , t ) + ( Δ t ) 2 2 D 2 f α ( e q ) ( r , t ) + O ( ( Δ t ) 3 ) = f α ( e q ) ( r , t ) − ( Δ t ) D f α ( e q ) ( r , t ) − Δ t 2 τ D f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) (11) \begin{aligned}
f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t) &= f_{\alpha}^{(eq)} (\boldsymbol{r}, t) - (\Delta t) \mathrm{D} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) \\
&\quad+ \frac{(\Delta t)^2}{2} \mathrm{D}^{2} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) + O((\Delta t)^3) \\
&= f_{\alpha}^{(eq)} (\boldsymbol{r}, t) - (\Delta t) \mathrm{D} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) \\
&\quad - \frac{\Delta t}{2 \tau} \mathrm{D} f_{\alpha}^{(neq)} (\boldsymbol{r}, t) + O((\Delta t)^3) \\
\end{aligned}
\tag{11}
f α ( e q ) ( r − c α Δ t , t − Δ t ) = f α ( e q ) ( r , t ) − ( Δ t ) D f α ( e q ) ( r , t ) + 2 ( Δ t ) 2 D 2 f α ( e q ) ( r , t ) + O ( ( Δ t ) 3 ) = f α ( e q ) ( r , t ) − ( Δ t ) D f α ( e q ) ( r , t ) − 2 τ Δ t D f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) ( 1 1 )
将式 (11) 代入式 (8) 的第一条公式，有：
ρ ∗ = ∑ α f α ( e q ) ( r , t ) − ( Δ t ) ∑ α D f α ( e q ) ( r , t ) − Δ t 2 τ ∑ α D f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) (12) \begin{aligned}
\rho^{*} &= \sum_{\alpha}f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \sum_{\alpha} \mathrm{D} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \\
&\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \mathrm{D} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3)
\end{aligned}
\tag{12}
ρ ∗ = α ∑ f α ( e q ) ( r , t ) − ( Δ t ) α ∑ D f α ( e q ) ( r , t ) − 2 τ Δ t α ∑ D f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) ( 1 2 )
式 (12) 右侧第一项与 ρ ∗ \rho^{*} ρ ∗ Δ t \Delta t Δ t ρ ∗ \rho^{*} ρ ∗ O ( ( Δ t ) 2 ) O((\Delta t)^2) O ( ( Δ t ) 2 )
将式 (11) 代入式 (8) 的第二条公式，整理得：
ρ ∗ u ∗ = ∑ α c α f α ( e q ) ( r , t ) − ( Δ t ) ∑ α c α D f α ( e q ) ( r , t ) − Δ t 2 τ ∑ α c α D f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) = ∑ α c α f α ( e q ) ( r , t ) − ( Δ t ) [ ∂ ∂ t ∑ α c α f α ( e q ) ( r , t ) + ∇ ⋅ ∑ α ( c α ) β ( c α ) γ ( f α ( e q ) ( r , t ) + 1 2 τ f α ( n e q ) ( r , t ) ) ] − Δ t 2 τ ∑ α ∂ ∂ t c α f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) (13) \begin{aligned}
\rho^{*} \boldsymbol{u}^{*} &= \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \sum_{\alpha} \boldsymbol{c}_{\alpha} \mathrm{D} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \\
&\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \boldsymbol{c}_{\alpha} \mathrm{D} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3) \\
&= \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \left[ \frac{\partial}{\partial t} \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) + \right. \\
&\quad \left. \nabla \cdot \sum_{\alpha} (\boldsymbol{c}_{\alpha})_\beta (\boldsymbol{c}_{\alpha})_\gamma \left( f_{\alpha}^{(eq)}(\boldsymbol{r}, t) + \frac{1}{2 \tau} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) \right) \right] \\
&\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \frac{\partial}{\partial t} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3)
\end{aligned}
\tag{13}
ρ ∗ u ∗ = α ∑ c α f α ( e q ) ( r , t ) − ( Δ t ) α ∑ c α D f α ( e q ) ( r , t ) − 2 τ Δ t α ∑ c α D f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) = α ∑ c α f α ( e q ) ( r , t ) − ( Δ t ) [ ∂ t ∂ α ∑ c α f α ( e q ) ( r , t ) + ∇ ⋅ α ∑ ( c α ) β ( c α ) γ ( f α ( e q ) ( r , t ) + 2 τ 1 f α ( n e q ) ( r , t ) ) ] − 2 τ Δ t α ∑ ∂ t ∂ c α f α ( n e q ) ( r , t ) + O ( ( Δ t ) 3 ) ( 1 3 )
式 (13) 右侧第一项与 ρ ∗ u ∗ \rho^{*} \boldsymbol{u}^{*} ρ ∗ u ∗ ∑ α ∂ ∂ t c α f α ( n e q ) ( r , t ) \sum\limits_{\alpha} \frac{\partial}{\partial t} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) α ∑ ∂ t ∂ c α f α ( n e q ) ( r , t ) Δ t \Delta t Δ t ρ ∗ u ∗ \rho^{*} \boldsymbol{u}^{*} ρ ∗ u ∗ O ( ( Δ t ) 2 ) O((\Delta t)^2) O ( ( Δ t ) 2 )
式 (9) 对应的宏观方程为
{ ∂ ρ ∂ t = 0 ∂ ρ u ∂ t + ∇ ⋅ [ ( 1 − 1 2 τ ) ∑ α ( c α ) β ( c α ) γ f α ( n e q ) ] = 0 (14) \begin{cases}
\frac{\partial \rho}{\partial t} = 0 \\
\frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \left[ \left( 1 - \frac{1}{2 \tau} \right) \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)} \right] = \boldsymbol{0}
\end{cases}
\tag{14}
⎩ ⎪ ⎨ ⎪ ⎧ ∂ t ∂ ρ = 0 ∂ t ∂ ρ u + ∇ ⋅ [ ( 1 − 2 τ 1 ) α ∑ ( c α ) β ( c α ) γ f α ( n e q ) ] = 0 ( 1 4 )
这里主要介绍一下动量方程的还原。对 f α ( n e q ) ( r ± Δ t 2 c α , t − Δ t 2 ) f_{\alpha}^{(neq)}(\boldsymbol{r} \pm \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) f α ( n e q ) ( r ± 2 Δ t c α , t − 2 Δ t )
f α ( n e q ) ( r ± Δ t 2 c α , t − Δ t 2 ) = f α ( n e q ) ( r , t − Δ t 2 ) ± Δ t 2 ( c α ⋅ ∇ ) f α ( n e q ) ( r , t − Δ t 2 ) + ( Δ t ) 2 8 ( c α ⋅ ∇ ) 2 f α ( n e q ) ( r , t − Δ t 2 ) + O ( ( Δ t ) 3 ) \begin{aligned}
f_{\alpha}^{(neq)}(\boldsymbol{r} \pm \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) &= f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \\
&\quad \pm \frac{\Delta t}{2} (\boldsymbol{c}_{\alpha} \cdot \nabla) f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \\
&\quad + \frac{(\Delta t)^2}{8} (\boldsymbol{c}_{\alpha} \cdot \nabla)^2 f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) + O((\Delta t)^3)
\end{aligned}
f α ( n e q ) ( r ± 2 Δ t c α , t − 2 Δ t ) = f α ( n e q ) ( r , t − 2 Δ t ) ± 2 Δ t ( c α ⋅ ∇ ) f α ( n e q ) ( r , t − 2 Δ t ) + 8 ( Δ t ) 2 ( c α ⋅ ∇ ) 2 f α ( n e q ) ( r , t − 2 Δ t ) + O ( ( Δ t ) 3 )
这里我们以未化简的式 (9-1) 为例。在式 (9-1) 的动量校正步中，有：
( 1 − 1 τ ) 1 Δ t ∑ α c α [ f α ( n e q ) ( r + Δ t 2 c α , t − Δ t 2 ) − f α ( n e q ) ( r − Δ t 2 c α , t − Δ t 2 ) ] = ( 1 − 1 τ ) ∑ α c α ( c α ⋅ ∇ ) f α ( n e q ) ( r , t − Δ t 2 ) + O ( ( Δ t ) 2 ) = ∇ ⋅ [ ( 1 − 1 τ ) ∑ α ( c α ) β ( c α ) γ f α ( n e q ) ( r , t − Δ t 2 ) ] + O ( ( Δ t ) 2 ) \begin{aligned}
&\quad \left(1 - \frac{1}{\tau}\right) \frac{1}{\Delta t} \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} [f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) - f_{\alpha}^{(neq)}(\boldsymbol{r}-\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2})] \\
&= \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} (\boldsymbol{c}_{\alpha} \cdot \nabla) f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) + O((\Delta t)^2) \\
&= \nabla \cdot \left[ \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \right] + O((\Delta t)^2) \\
\end{aligned}
( 1 − τ 1 ) Δ t 1 α ∑ c α [ f α ( n e q ) ( r + 2 Δ t c α , t − 2 Δ t ) − f α ( n e q ) ( r − 2 Δ t c α , t − 2 Δ t ) ] = ( 1 − τ 1 ) α ∑ c α ( c α ⋅ ∇ ) f α ( n e q ) ( r , t − 2 Δ t ) + O ( ( Δ t ) 2 ) = ∇ ⋅ [ ( 1 − τ 1 ) α ∑ ( c α ) β ( c α ) γ f α ( n e q ) ( r , t − 2 Δ t ) ] + O ( ( Δ t ) 2 )
此时联立式 (13) ，即可还原式 (14) 中的动量方程。
这里令 Y = ( y 1 , y 2 , y 3 , y 4 ) T \mathbf{Y} = (y_1, y_2, y_3, y_4)^T Y = ( y 1 , y 2 , y 3 , y 4 ) T
y 1 = ρ , y 2 = ρ u x , y 3 = ρ u y , y 4 = ρ u z y_1 = \rho,\quad y_2 = \rho u_x ,\quad y_3 = \rho u_y ,\quad y_4 = \rho u_z
y 1 = ρ , y 2 = ρ u x , y 3 = ρ u y , y 4 = ρ u z
将这些符号代换到式 (2) 的 f α ( e q ) f_{\alpha}^{(eq)} f α ( e q )
f α ( e q ) = y 1 α w α + w α c s 2 ( c α x y 2 α + c α y y 3 α + c α z y 4 α ) + w α 2 y 1 α c s 4 ( c α x y 2 α + c α y y 3 α + c α z y 4 α ) 2 − w α 2 y 1 α c s 2 ( y 2 α 2 + y 3 α 2 + y 4 α 2 ) (15) \begin{aligned}
f_{\alpha}^{(eq)} &= y_{1 \alpha} w_{\alpha} + \frac{w_\alpha}{c_s^2} ( c_{\alpha x} y_{2 \alpha}+ c_{\alpha y} y_{3 \alpha} + c_{\alpha z} y_{4 \alpha} ) \\
&\quad + \frac{w_\alpha}{2 y_{1 \alpha} c_s^4} ( c_{\alpha x} y_{2 \alpha}+ c_{\alpha y} y_{3 \alpha} + c_{\alpha z} y_{4 \alpha} )^2 \\
&\quad - \frac{w_\alpha}{2 y_{1 \alpha} c_s^2} (y_{2 \alpha}^2 + y_{3 \alpha}^2 + y_{4 \alpha}^2)
\end{aligned}
\tag{15}
f α ( e q ) = y 1 α w α + c s 2 w α ( c α x y 2 α + c α y y 3 α + c α z y 4 α ) + 2 y 1 α c s 4 w α ( c α x y 2 α + c α y y 3 α + c α z y 4 α ) 2 − 2 y 1 α c s 2 w α ( y 2 α 2 + y 3 α 2 + y 4 α 2 ) ( 1 5 )
where the subscript α \alpha α ( r − c α Δ t ) (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t) ( r − c α Δ t )
对于第 n n n Y n \mathbf{Y}^{n} Y n Y ∗ = G ( Y n ) \mathbf{Y}^{*} = \mathbf{G}(\mathbf{Y}^{n}) Y ∗ = G ( Y n )
为便于进行诺依曼稳定性分析，对上式进行求导，有
d Y ∗ = [ ∂ Y ∗ ∂ y 1 ] n d y 1 n + [ ∂ Y ∗ ∂ y 2 ] n d y 2 n + [ ∂ Y ∗ ∂ y 3 ] n d y 3 n + [ ∂ Y ∗ ∂ y 4 ] n d y 4 n \mathrm{d \mathbf{Y}^{*}} = \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_1} \right]^n \mathrm{d}y_1^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_2} \right]^n \mathrm{d}y_2^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_3} \right]^n \mathrm{d}y_3^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_4} \right]^n \mathrm{d}y_4^n
d Y ∗ = [ ∂ y 1 ∂ Y ∗ ] n d y 1 n + [ ∂ y 2 ∂ Y ∗ ] n d y 2 n + [ ∂ y 3 ∂ Y ∗ ] n d y 3 n + [ ∂ y 4 ∂ Y ∗ ] n d y 4 n
上标 n n n n n n d Y ∗ = ( d y 1 ∗ , d y 2 ∗ , d y 3 ∗ , d y 4 ∗ ) T \mathrm{d \mathbf{Y}^{*}} = (\mathrm{d} y_1^{*}, \mathrm{d} y_2^{*}, \mathrm{d} y_3^{*}, \mathrm{d} y_4^{*})^T d Y ∗ = ( d y 1 ∗ , d y 2 ∗ , d y 3 ∗ , d y 4 ∗ ) T d y λ ∗ = ∑ κ = 1 4 [ ∂ y λ ∗ ∂ y κ ] n d y κ n \mathrm{d} y_{\lambda}^{*} = \sum_{\kappa = 1}^{4}\left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n \mathrm{d}y_{\kappa}^n d y λ ∗ = ∑ κ = 1 4 [ ∂ y κ ∂ y λ ∗ ] n d y κ n λ , κ ∈ [ 1 , 2 , 3 , 4 ] \lambda, \kappa \in [1,2,3,4] λ , κ ∈ [ 1 , 2 , 3 , 4 ]
对于本文列出的方程组， [ ∂ y λ ∗ ∂ y κ ] n \left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n [ ∂ y κ ∂ y λ ∗ ] n [ ∂ y 1 ∗ ∂ y 1 ] n \left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n [ ∂ y 1 ∂ y 1 ∗ ] n
[ ∂ y 1 ∗ ∂ y 1 ] n = ∂ ∂ y 1 ∑ α f α ( e q ) = ∑ α [ w α − w α 2 ( y 1 α n ) 2 c s 4 ( c α x y 2 α n + c α y y 3 α n + c α z y 4 α n ) 2 + w α 2 ( y 1 α n ) 2 c s 2 ( ( y 2 α n ) 2 + ( y 3 α n ) 2 + ( y 4 α n ) 2 ) ] (16) \begin{aligned}
\left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n
&= \frac{\partial}{\partial y_{1}} \sum_{\alpha} f_{\alpha}^{(eq)} \\
&= \sum_{\alpha} \left[ w_{\alpha} - \frac{w_\alpha}{2 (y_{1\alpha}^n)^2 c_s^4} ( c_{\alpha x} y_{2 \alpha}^n + c_{\alpha y} y_{3 \alpha}^n + c_{\alpha z} y_{4 \alpha}^n )^2 \right. \\
&\quad\left. + \frac{w_\alpha}{2 (y_{1\alpha}^n)^2 c_s^2} ((y_{2 \alpha}^n)^2 + (y_{3 \alpha}^n)^2 + (y_{4 \alpha}^n)^2) \right]
\end{aligned}
\tag{16}
[ ∂ y 1 ∂ y 1 ∗ ] n = ∂ y 1 ∂ α ∑ f α ( e q ) = α ∑ [ w α − 2 ( y 1 α n ) 2 c s 4 w α ( c α x y 2 α n + c α y y 3 α n + c α z y 4 α n ) 2 + 2 ( y 1 α n ) 2 c s 2 w α ( ( y 2 α n ) 2 + ( y 3 α n ) 2 + ( y 4 α n ) 2 ) ] ( 1 6 )
假设式 (16) 中的 d y i ∗ \mathrm{d} y^{*}_{i} d y i ∗ d y i n \mathrm{d} y^{n}_{i} d y i n y i α n y^{n}_{i \alpha} y i α n i = [ 1 , 2 , 3 , 4 ] i = [1,2,3,4] i = [ 1 , 2 , 3 , 4 ] j = − 1 j=\sqrt{-1} j = − 1
d y i ∗ = d y i ‾ ∗ exp ( j k ⋅ r ) , d y i n = d y i ‾ n exp ( j k ⋅ r ) , y i α n = y i n exp ( j k ⋅ r ) \mathrm{d} y^{*}_{i} = \overline{\mathrm{d} y_{i}}^{*} \exp(j \boldsymbol{k \cdot r}) ,\quad
\mathrm{d} y^{n}_{i} = \overline{\mathrm{d} y_{i}}^{n} \exp(j \boldsymbol{k \cdot r}) ,\quad
y^{n}_{i \alpha} = y_{i}^{n} \exp(j \boldsymbol{k \cdot r})
d y i ∗ = d y i ∗ exp ( j k ⋅ r ) , d y i n = d y i n exp ( j k ⋅ r ) , y i α n = y i n exp ( j k ⋅ r )
代入到式 (16) 和 d Y ∗ \mathrm{d \mathbf{Y}^{*}} d Y ∗
[ ∂ y 1 ∗ ∂ y 1 ] n = ∂ ∂ y 1 ∑ α f α ( e q ) = ∑ α [ w α − w α 2 ( y 1 n ) 2 c s 4 ( c α x y 2 n + c α y y 3 n + c α z y 4 n ) 2 + w α 2 ( y 1 n ) 2 c s 2 ( ( y 2 n ) 2 + ( y 3 n ) 2 + ( y 4 n ) 2 ) ] (17) \begin{aligned}
\left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n
&= \frac{\partial}{\partial y_{1}} \sum_{\alpha} f_{\alpha}^{(eq)} \\
&= \sum_{\alpha} \left[ w_{\alpha} - \frac{w_\alpha}{2 (y_{1}^n)^2 c_s^4} ( c_{\alpha x} y_{2}^n + c_{\alpha y} y_{3}^n + c_{\alpha z} y_{4}^n )^2 \right. \\
&\quad\left. + \frac{w_\alpha}{2 (y_{1}^n)^2 c_s^2} ((y_{2}^n)^2 + (y_{3}^n)^2 + (y_{4}^n)^2) \right]
\end{aligned}
\tag{17}
[ ∂ y 1 ∂ y 1 ∗ ] n = ∂ y 1 ∂ α ∑ f α ( e q ) = α ∑ [ w α − 2 ( y 1 n ) 2 c s 4 w α ( c α x y 2 n + c α y y 3 n + c α z y 4 n ) 2 + 2 ( y 1 n ) 2 c s 2 w α ( ( y 2 n ) 2 + ( y 3 n ) 2 + ( y 4 n ) 2 ) ] ( 1 7 )
格子张量的一些性质
∑ α w α = 1. α ∈ [ 0 , 1 , . . . , q − 1 ] ∑ α w α c α β = 0. β ∈ [ x , y , z ] ∑ α w α c α β c α γ = c s 2 δ β γ . β , γ ∈ [ x , y , z ] \begin{aligned}
& \sum_{\alpha} w_{\alpha} = 1 .&\qquad \alpha \in [0,1, ..., q-1]\\
& \sum_{\alpha} w_{\alpha} c_{\alpha \beta} = 0 .&\qquad \beta \in [x,y,z] \\
& \sum_{\alpha} w_{\alpha} c_{\alpha \beta} c_{\alpha \gamma} = c_s^2 \delta_{\beta \gamma} .&\qquad \beta, \gamma \in [x,y,z]
\end{aligned}
α ∑ w α = 1 . α ∑ w α c α β = 0 . α ∑ w α c α β c α γ = c s 2 δ β γ . α ∈ [ 0 , 1 , . . . , q − 1 ] β ∈ [ x , y , z ] β , γ ∈ [ x , y , z ]
注: q q q
基于格子张量的性质，式 (17) 右侧仅剩下 ∑ α w α \sum\limits_{\alpha} w_{\alpha} α ∑ w α
[ ∂ y 1 ∗ ∂ y 1 ] n = 1 \left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n = 1
[ ∂ y 1 ∂ y 1 ∗ ] n = 1
同理，其他项也可进行类似分析，并得出
[ ∂ y λ ∗ ∂ y κ ] n = δ λ κ \left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n = \delta_{\lambda \kappa}
[ ∂ y κ ∂ y λ ∗ ] n = δ λ κ
其中 δ λ κ \delta_{\lambda \kappa} δ λ κ λ , κ ∈ [ 1 , 2 , 3 , 4 ] \lambda,\kappa \in [1,2,3,4] λ , κ ∈ [ 1 , 2 , 3 , 4 ] d Y ‾ ∗ = d Y ‾ n \overline{\mathrm{d} \mathbf{Y}}^{*} = \overline{\mathrm{d} \mathbf{Y}}^{n} d Y ∗ = d Y n
令 d y i n + 1 = d y i ‾ n + 1 exp ( j k ⋅ r ) {\mathrm{d} y_i}^{n+1} = \overline{\mathrm{d} y_i}^{n+1} \exp(j \boldsymbol{k \cdot r}) d y i n + 1 = d y i n + 1 exp ( j k ⋅ r ) i ∈ [ 1 , 2 , 3 , 4 ] i \in [1,2,3,4] i ∈ [ 1 , 2 , 3 , 4 ] ρ \rho ρ
d y 1 ‾ n + 1 = d y 1 ‾ ∗ = d y 1 ‾ n \overline{\mathrm{d} y_1}^{n+1} = \overline{\mathrm{d} y_1}^{*} = \overline{\mathrm{d} y_1}^{n}
d y 1 n + 1 = d y 1 ∗ = d y 1 n
根据式 (9-2) 的校正步，以及 f α ( n e q ) ( r + Δ t 2 c α , t − Δ t 2 ) f_{\alpha}^{(neq)}(\boldsymbol{r} + \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) f α ( n e q ) ( r + 2 Δ t c α , t − 2 Δ t )
y 2 n + 1 = y 2 ∗ − ( τ − 1 ) y 2 n + ( τ − 1 ) ∑ α c α x f α ( e q ) ( r + c α Δ t , t ) y 3 n + 1 = y 3 ∗ − ( τ − 1 ) y 3 n + ( τ − 1 ) ∑ α c α y f α ( e q ) ( r + c α Δ t , t ) y 4 n + 1 = y 4 ∗ − ( τ − 1 ) y 4 n + ( τ − 1 ) ∑ α c α z f α ( e q ) ( r + c α Δ t , t ) (18) \begin{aligned}
y_2^{n+1} = y_2^{*} - (\tau - 1) y_2^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha x} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\
y_3^{n+1} = y_3^{*} - (\tau - 1) y_3^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha y} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\
y_4^{n+1} = y_4^{*} - (\tau - 1) y_4^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha z} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t)
\end{aligned}
\tag{18}
y 2 n + 1 = y 2 ∗ − ( τ − 1 ) y 2 n + ( τ − 1 ) α ∑ c α x f α ( e q ) ( r + c α Δ t , t ) y 3 n + 1 = y 3 ∗ − ( τ − 1 ) y 3 n + ( τ − 1 ) α ∑ c α y f α ( e q ) ( r + c α Δ t , t ) y 4 n + 1 = y 4 ∗ − ( τ − 1 ) y 4 n + ( τ − 1 ) α ∑ c α z f α ( e q ) ( r + c α Δ t , t ) ( 1 8 )
这里同样以 y 2 n + 1 y_2^{n+1} y 2 n + 1
S = ∑ α c α x f α ( e q ) ( r + c α Δ t , t ) = ∑ α c α x w α y 1 ( − α ) ∗ + ∑ α c α x w α c s 2 ( c α x y 2 ( − α ) ∗ + c α y y 3 ( − α ) ∗ + c α z y 4 ( − α ) ∗ ) + ∑ α c α x w α 2 c s 4 y 1 ( − α ) ∗ ( c α x y 2 ( − α ) ∗ + c α y y 3 ( − α ) ∗ + c α z y 4 ( − α ) ∗ ) 2 − ∑ α c α x w α 2 c s 2 y 1 ( − α ) ∗ ( ( y 2 ( − α ) ∗ ) 2 + ( y 3 ( − α ) ∗ ) 2 + ( y 4 ( − α ) ∗ ) 2 ) 2 \begin{aligned}
S &= \sum\limits_{\alpha} c_{\alpha x} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\
&= \sum\limits_{\alpha} c_{\alpha x} w_{\alpha} y_{1(-\alpha)}^{*} \\
&\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{c_s^2} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*}) \\
&\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 y_{1(-\alpha)}^{*}} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*})^2 \\
&\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 y_{1(-\alpha)}^{*}} ( (y_{2(-\alpha)}^{*})^2 + (y_{3(-\alpha)}^{*})^2 + (y_{4(-\alpha)}^{*})^2 )^2 \\
\end{aligned}
S = α ∑ c α x f α ( e q ) ( r + c α Δ t , t ) = α ∑ c α x w α y 1 ( − α ) ∗ + α ∑ c s 2 c α x w α ( c α x y 2 ( − α ) ∗ + c α y y 3 ( − α ) ∗ + c α z y 4 ( − α ) ∗ ) + α ∑ 2 c s 4 y 1 ( − α ) ∗ c α x w α ( c α x y 2 ( − α ) ∗ + c α y y 3 ( − α ) ∗ + c α z y 4 ( − α ) ∗ ) 2 − α ∑ 2 c s 2 y 1 ( − α ) ∗ c α x w α ( ( y 2 ( − α ) ∗ ) 2 + ( y 3 ( − α ) ∗ ) 2 + ( y 4 ( − α ) ∗ ) 2 ) 2
其中下标 ( − α ) (-\alpha) ( − α ) r + c α Δ t \boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t r + c α Δ t
y 2 n + 1 = y 2 ∗ + ( τ − 1 ) ( S − y 2 n ) (19) y_2^{n+1} = y_2^{*} + (\tau - 1) (S - y_2^{n}) \tag{19}
y 2 n + 1 = y 2 ∗ + ( τ − 1 ) ( S − y 2 n ) ( 1 9 )
所以其差分表示为
d y 2 n + 1 = d y 2 ∗ + ( τ − 1 ) ( d S − d y 2 n ) (20) \mathrm{d} y_2^{n+1} = \mathrm{d} y_2^{*} + (\tau - 1) (\mathrm{d} S - \mathrm{d} y_2^{n}) \tag{20}
d y 2 n + 1 = d y 2 ∗ + ( τ − 1 ) ( d S − d y 2 n ) ( 2 0 )
对于 d S \mathrm{d} S d S
d S = ( ∂ S ∂ y 1 ) ∗ d y 1 ∗ + ( ∂ S ∂ y 2 ) ∗ d y 2 ∗ + ( ∂ S ∂ y 3 ) ∗ d y 3 ∗ + ( ∂ S ∂ y 4 ) ∗ d y 4 ∗ \mathrm{d} S = \left( \frac{\partial S}{\partial y_1} \right)^{*} \mathrm{d}y_1^{*} + \left( \frac{\partial S}{\partial y_2} \right)^{*} \mathrm{d}y_2^{*} + \left( \frac{\partial S}{\partial y_3} \right)^{*} \mathrm{d}y_3^{*} + \left( \frac{\partial S}{\partial y_4} \right)^{*} \mathrm{d}y_4^{*}
d S = ( ∂ y 1 ∂ S ) ∗ d y 1 ∗ + ( ∂ y 2 ∂ S ) ∗ d y 2 ∗ + ( ∂ y 3 ∂ S ) ∗ d y 3 ∗ + ( ∂ y 4 ∂ S ) ∗ d y 4 ∗
以 ( ∂ S ∂ y 1 ) ∗ \left( \frac{\partial S}{\partial y_1} \right)^{*} ( ∂ y 1 ∂ S ) ∗
( ∂ S ∂ y 1 ) ∗ = ∑ α c α x w α − ∑ α c α x w α 2 c s 4 ( y 1 ( − α ) ∗ ) 2 ( c α x y 2 ( − α ) ∗ + c α y y 3 ( − α ) ∗ + c α z y 4 ( − α ) ∗ ) 2 + ∑ α c α x w α 2 c s 2 ( y 1 ( − α ) ∗ ) 2 ( ( y 2 ( − α ) ∗ ) 2 + ( y 3 ( − α ) ∗ ) 2 + ( y 4 ( − α ) ∗ ) 2 ) 2 \begin{aligned}
\left( \frac{\partial S}{\partial y_1} \right)^{*} &=
\sum\limits_{\alpha} c_{\alpha x} w_{\alpha} \\
&\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 (y_{1(-\alpha)}^{*})^2} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*})^2 \\
&\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 (y_{1(-\alpha)}^{*})^2} ( (y_{2(-\alpha)}^{*})^2 + (y_{3(-\alpha)}^{*})^2 + (y_{4(-\alpha)}^{*})^2 )^2
\end{aligned}
( ∂ y 1 ∂ S ) ∗ = α ∑ c α x w α − α ∑ 2 c s 4 ( y 1 ( − α ) ∗ ) 2 c α x w α ( c α x y 2 ( − α ) ∗ + c α y y 3 ( − α ) ∗ + c α z y 4 ( − α ) ∗ ) 2 + α ∑ 2 c s 2 ( y 1 ( − α ) ∗ ) 2 c α x w α ( ( y 2 ( − α ) ∗ ) 2 + ( y 3 ( − α ) ∗ ) 2 + ( y 4 ( − α ) ∗ ) 2 ) 2
与预报步的做法类似，将 y i ( − α ) ∗ y_{i(-\alpha)}^{*} y i ( − α ) ∗ y i ∗ y_{i}^{*} y i ∗
( ∂ S ∂ y 1 ) ∗ = ∑ α c α x w α − ∑ α c α x w α 2 c s 4 ( y 1 ∗ ) 2 ( c α x y 2 ∗ + c α y y 3 ∗ + c α z y 4 ∗ ) 2 + ∑ α c α x w α 2 c s 2 ( y 1 ∗ ) 2 ( ( y 2 ∗ ) 2 + ( y 3 ∗ ) 2 + ( y 4 ∗ ) 2 ) 2 \begin{aligned}
\left( \frac{\partial S}{\partial y_1} \right)^{*} &=
\sum\limits_{\alpha} c_{\alpha x} w_{\alpha} \\
&\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 (y_{1}^{*})^2} (c_{\alpha x} y_{2}^{*} + c_{\alpha y} y_{3}^{*} + c_{\alpha z} y_{4}^{*})^2 \\
&\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 (y_{1}^{*})^2} ( (y_{2}^{*})^2 + (y_{3}^{*})^2 + (y_{4}^{*})^2 )^2
\end{aligned}
( ∂ y 1 ∂ S ) ∗ = α ∑ c α x w α − α ∑ 2 c s 4 ( y 1 ∗ ) 2 c α x w α ( c α x y 2 ∗ + c α y y 3 ∗ + c α z y 4 ∗ ) 2 + α ∑ 2 c s 2 ( y 1 ∗ ) 2 c α x w α ( ( y 2 ∗ ) 2 + ( y 3 ∗ ) 2 + ( y 4 ∗ ) 2 ) 2
此时所有 y i ∗ y_{i}^{*} y i ∗ ( ∂ S ∂ y 1 ) ∗ = 0 \left( \frac{\partial S}{\partial y_1} \right)^{*} = 0 ( ∂ y 1 ∂ S ) ∗ = 0
同理
( ∂ S ∂ y 2 ) ∗ = 1 , ( ∂ S ∂ y 3 ) ∗ = ( ∂ S ∂ y 4 ) ∗ = 0 \left( \frac{\partial S}{\partial y_2} \right)^{*} = 1,\quad
\left( \frac{\partial S}{\partial y_3} \right)^{*} = \left( \frac{\partial S}{\partial y_4} \right)^{*} = 0
( ∂ y 2 ∂ S ) ∗ = 1 , ( ∂ y 3 ∂ S ) ∗ = ( ∂ y 4 ∂ S ) ∗ = 0
此时 d S ‾ = d y 2 ‾ ∗ \overline{\mathrm{d} S } = \overline{\mathrm{d} y_2}^{*} d S = d y 2 ∗
d y 2 ‾ n + 1 = d y 2 ‾ ∗ = d y 2 ‾ n \overline{\mathrm{d} y_2}^{n+1} = \overline{\mathrm{d} y_2}^{*} = \overline{\mathrm{d} y_2}^{n}
d y 2 n + 1 = d y 2 ∗ = d y 2 n
同理，对 y 3 n + 1 y_3^{n+1} y 3 n + 1 y 4 n + 1 y_4^{n+1} y 4 n + 1 d y 3 ‾ n + 1 = d y 3 ‾ n \overline{\mathrm{d} y_3}^{n+1} = \overline{\mathrm{d} y_3}^{n} d y 3 n + 1 = d y 3 n d y 4 ‾ n + 1 = d y 4 ‾ n \overline{\mathrm{d} y_4}^{n+1} = \overline{\mathrm{d} y_4}^{n} d y 4 n + 1 = d y 4 n
记 d Y n = ( d y 1 n , d y 2 n , d y 3 n , d y 4 n ) T \mathrm{d \mathbf{Y}^{n}} = (\mathrm{d} y_1^{n}, \mathrm{d} y_2^{n}, \mathrm{d} y_3^{n}, \mathrm{d} y_4^{n})^T d Y n = ( d y 1 n , d y 2 n , d y 3 n , d y 4 n ) T d Y n + 1 = ( d y 1 n + 1 , d y 2 n + 1 , d y 3 n + 1 , d y 4 n + 1 ) T \mathrm{d \mathbf{Y}^{n+1}} = (\mathrm{d} y_1^{n+1}, \mathrm{d} y_2^{n+1}, \mathrm{d} y_3^{n+1}, \mathrm{d} y_4^{n+1})^T d Y n + 1 = ( d y 1 n + 1 , d y 2 n + 1 , d y 3 n + 1 , d y 4 n + 1 ) T
d Y ‾ n + 1 = I d Y ‾ n \overline{\mathrm{d \mathbf{Y}}}^{n+1} = \mathbf{I} \, \overline{\mathrm{d \mathbf{Y}}}^{n}
d Y n + 1 = I d Y n
其中 I \mathbf{I} I