注:该文章同步发在知乎专栏 和CSDN
这里所说的 “简化的格子 Boltzmann 方法” 其实是 Simplified and Highly Stable Lattice Boltzmann Method (SHSLBM) (, , , )。该方法的特色是仅使用平衡态分布函数、宏观密度、宏观速度进行 LBM 计算,并保证无条件稳定。
下面的内容只是我对相关文献进行简单阅读后的部分整理
单松弛格子 Boltzmann 方法
基本方程
针对流场的 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 τ 为 BGK 碰撞项的松弛时间。 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 )
Chapman-Enskog 展开
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
SHSLBM
式 (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) 中由 − 1 τ \frac{-1}{\tau} τ − 1 松弛的部分即为 f α ( n e q ) f_{\alpha}^{(neq)} f α ( n e q ) ,因此式 (1) 理应可以被写作完全由 f α ( e q ) f_{\alpha}^{(eq)} f α ( e q ) 进行表示的形式。这便是 SHSLBM 的核心思想。
基本算法
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 ) 都是未知的。参照Shu等的做法,对上面两条公式中的这两个量使用预报步数据进行近似,即
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 ) 进行展开,并代入式 (7) ,得到:
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^{*} ρ ∗ 相等,右侧第三项根据非平衡态的守恒律可知其为 0 。所以在消去式 (12) 中的 Δ t \Delta t Δ t 后, ρ ∗ \rho^{*} ρ ∗ 的计算与式 (10) 中第一条公式对应,并具有 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 ) 项根据非平衡态的守恒律可知其为 0 。所以在消去式 (13) 中的 Δ t \Delta t Δ t 后, ρ ∗ u ∗ \rho^{*} \boldsymbol{u}^{*} ρ ∗ u ∗ 的计算与式 (10) 中第二条公式对应,并具有 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 ) 进行 Taylor 展开,得到
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) 中的动量方程。
SHSLBM 的稳定性分析
这里令 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 α denotes the quantities at a space level of ( r − c α Δ t ) (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t) ( r − c α Δ t ) .
预报步
对于第 n n n 步的结果 Y n \mathbf{Y}^{n} Y n ,式 (8) 的预报步写作 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 多达 16 个。因此下文仅拿 [ ∂ 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 ] )。由于校正步(式 (9-2) )中并未对 ρ \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 ∗ ,则式 (20) 可化简为:
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 为单位矩阵。显然,特征矩阵的特征值全是 1 ,代表着该格式无条件稳定。