这里所说的 “简化的格子 Boltzmann 方法” 其实是 Simplified and Highly Stable Lattice Boltzmann Method (SHSLBM)(1,2,3,4)。该方法的特色是仅使用平衡态分布函数、宏观密度、宏观速度进行 LBM 计算,并保证无条件稳定。
下面的内容只是我对相关文献进行简单阅读后的部分整理
Simplified and Highly Stable Lattice Boltzmann Method
单松弛格子 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α(eq)](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α(eq) 为 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α(eq)=ρwα[1+cs2cα⋅u+2cs4(cα⋅u)2−2cs2u⋅u](2)
式 (2) 中 ρ \rho ρ 为流体密度, u = ( u x , u y , u z ) T \boldsymbol{u} = (u_x, u_y, u_z)^T u=(ux,uy,uz)T 为流体宏观速度, w α w_\alpha wα 为 c α \boldsymbol{c}_\alpha cα 的权重系数, c s c_s cs 为格子声速。 w α w_\alpha wα 、 c α \boldsymbol{c}_\alpha cα 和 c s c_s cs 的数值由选取的离散速度模型决定。
对于 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α(eq),ρu=α∑cαfα=α∑cαfα(eq)(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α(eq)+ϵfα(1)+ϵ2fα(2)+...,∂t∂=ϵ∂t0∂+ϵ2∂t1∂,∂r∂=ϵ∂r0∂
上述展开保留至
ϵ
1
\epsilon^{1}
ϵ1 ,即可还原为5
∂
ρ
∂
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α(eq))=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α(eq)+(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α(eq)=fα(neq) 近似为非平衡态分布函数。此时
Π
β
γ
=
∑
α
(
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α(eq)+(1−2τ1)fα(neq)](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α(neq)≈ϵfα(1)=−τ[∂t∂+cα⋅∇]fα(eq)Δ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α(neq)≈−τ⋅(Δt)⋅Dfα(eq)。
流体运动粘度
ν
\nu
ν 和松弛时间
τ
\tau
τ 的关系为
ν
=
(
τ
−
1
2
)
c
s
2
Δ
t
\nu = \left( \tau - \frac{1}{2} \right) c_s^2 \Delta t
ν=(τ−21)cs2Δt
SHSLBM
式 (7) 说明:在保留至 ϵ 1 \epsilon^{1} ϵ1 的情况下, f α ( n e q ) f_{\alpha}^{(neq)} fα(neq) 和 f α ( e q ) f_{\alpha}^{(eq)} fα(eq) 存在联系。考虑到式 (1) 中由 − 1 τ \frac{-1}{\tau} τ−1 松弛的部分即为 f α ( n e q ) f_{\alpha}^{(neq)} fα(neq) ,因此式 (1) 理应可以被写作完全由 f α ( e q ) f_{\alpha}^{(eq)} fα(eq) 进行表示的形式。这便是 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α(eq)(r−cαΔt,t−Δt),ρ∗u∗=α∑cαfα(eq)(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α(neq)(r+2Δtcα,t−2Δt)−fα(neq)(r−2Δtcα,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α(neq)(r−2Δtcα,t−2Δt)=−τ[fα(eq)(r,t)−fα(eq)(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α(neq)(r+2Δtcα,t−2Δt)=−τ[fα(eq)(r+cαΔt,t)−fα(eq)(r,t−Δt)]+O((Δt)3)
在这两条式子中, f α ( e q ) ( r , t ) f_{\alpha}^{(eq)}(\boldsymbol{r}, t) fα(eq)(r,t) 和 f α ( e q ) ( r + c α Δ t , t ) f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) fα(eq)(r+cαΔt,t) 都是未知的。参照Shu等6的做法,对上面两条公式中的这两个量使用预报步数据进行近似,即
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α(eq)(r,t)=fα(eq)(ρ∗(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α(neq)(r−2Δtcα,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α(neq)(r−2Δtcα,t−2Δt)=−τ[α∑cαfα(eq)(r,t)−α∑cαfα(eq)(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α(neq)(r+2Δtcα,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α(eq))=0∂t∂ρu+∇⋅[α∑(cα)β(cα)γfα(eq)+2τ1α∑(cα)β(cα)γfα(neq)]=0(10)
首先,对 f α ( e q ) ( r − c α Δ t , t − Δ t ) f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t) fα(eq)(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{align} 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{align} \tag{11} fα(eq)(r−cαΔt,t−Δt)=fα(eq)(r,t)−(Δt)Dfα(eq)(r,t)+2(Δt)2D2fα(eq)(r,t)+O((Δt)3)=fα(eq)(r,t)−(Δt)Dfα(eq)(r,t)−2τΔtDfα(neq)(r,t)+O((Δt)3)(11)
将式 (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{align} \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{align} \tag{12}
ρ∗=α∑fα(eq)(r,t)−(Δt)α∑Dfα(eq)(r,t)−2τΔtα∑Dfα(neq)(r,t)+O((Δt)3)(12)
式 (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{align} \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{align} \tag{13}
ρ∗u∗=α∑cαfα(eq)(r,t)−(Δt)α∑cαDfα(eq)(r,t)−2τΔtα∑cαDfα(neq)(r,t)+O((Δt)3)=α∑cαfα(eq)(r,t)−(Δt)[∂t∂α∑cαfα(eq)(r,t)+∇⋅α∑(cα)β(cα)γ(fα(eq)(r,t)+2τ1fα(neq)(r,t))]−2τΔtα∑∂t∂cαfα(neq)(r,t)+O((Δt)3)(13)
式 (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α(neq)(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α(neq)]=0(14)
这里主要介绍一下动量方程的还原。对 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α(neq)(r±2Δtcα,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α(neq)(r±2Δtcα,t−2Δt)=fα(neq)(r,t−2Δt)±2Δt(cα⋅∇)fα(neq)(r,t−2Δt)+8(Δt)2(cα⋅∇)2fα(neq)(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)Δt1α∑cα[fα(neq)(r+2Δtcα,t−2Δt)−fα(neq)(r−2Δtcα,t−2Δt)]=(1−τ1)α∑cα(cα⋅∇)fα(neq)(r,t−2Δt)+O((Δt)2)=∇⋅[(1−τ1)α∑(cα)β(cα)γfα(neq)(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=(y1,y2,y3,y4)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 y1=ρ,y2=ρux,y3=ρuy,y4=ρuz
将这些符号代换到式 (2) 的 f α ( e q ) f_{\alpha}^{(eq)} fα(eq) 中,得:
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{align} 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{align} \tag{15} fα(eq)=y1αwα+cs2wα(cαxy2α+cαyy3α+cαzy4α)+2y1αcs4wα(cαxy2α+cαyy3α+cαzy4α)2−2y1αcs2wα(y2α2+y3α2+y4α2)(15)
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).2
预报步
对于第 n n n 步的结果 Y n \mathbf{Y}^{n} Yn,式 (8) 的预报步写作 Y ∗ = G ( Y n ) \mathbf{Y}^{*} = \mathbf{G}(\mathbf{Y}^{n}) Y∗=G(Yn)。
为便于进行诺依曼稳定性分析,对上式进行求导,有
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
dY∗=[∂y1∂Y∗]ndy1n+[∂y2∂Y∗]ndy2n+[∂y3∂Y∗]ndy3n+[∂y4∂Y∗]ndy4n
上标 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 dY∗=(dy1∗,dy2∗,dy3∗,dy4∗)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 dyλ∗=∑κ=14[∂yκ∂yλ∗]ndyκ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 [∂y1∂y1∗]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{align} \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{align} \tag{16} [∂y1∂y1∗]n=∂y1∂α∑fα(eq)=α∑[wα−2(y1αn)2cs4wα(cαxy2αn+cαyy3αn+cαzy4αn)2+2(y1αn)2cs2wα((y2αn)2+(y3αn)2+(y4αn)2)](16)
假设式 (16) 中的 d y i ∗ \mathrm{d} y^{*}_{i} dyi∗, d y i n \mathrm{d} y^{n}_{i} dyin, y i α n y^{n}_{i \alpha} yiα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}) dyi∗=dyi∗exp(jk⋅r),dyin=dyinexp(jk⋅r),yiαn=yinexp(jk⋅r)
代入到式 (16) 和 d Y ∗ \mathrm{d \mathbf{Y}^{*}} dY∗ 表达式,得
[ ∂ 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{align} \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{align} \tag{17} [∂y1∂y1∗]n=∂y1∂α∑fα(eq)=α∑[wα−2(y1n)2cs4wα(cαxy2n+cαyy3n+cαzy4n)2+2(y1n)2cs2wα((y2n)2+(y3n)2+(y4n)2)](17)
格子张量的一些性质2
∑ α 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αγ=cs2δβγ.α∈[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 [∂y1∂y1∗]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} dY∗=dYn 。
校正步
令 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}) dyin+1=dyin+1exp(jk⋅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} dy1n+1=dy1∗=dy1n
根据式 (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α(neq)(r+2Δtcα,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} y2n+1=y2∗−(τ−1)y2n+(τ−1)α∑cαxfα(eq)(r+cαΔt,t)y3n+1=y3∗−(τ−1)y3n+(τ−1)α∑cαyfα(eq)(r+cαΔt,t)y4n+1=y4∗−(τ−1)y4n+(τ−1)α∑cαzfα(eq)(r+cαΔt,t)(18)
这里同样以 y 2 n + 1 y_2^{n+1} y2n+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αxfα(eq)(r+cαΔt,t)=α∑cαxwαy1(−α)∗+α∑cs2cαxwα(cαxy2(−α)∗+cαyy3(−α)∗+cαzy4(−α)∗)+α∑2cs4y1(−α)∗cαxwα(cαxy2(−α)∗+cαyy3(−α)∗+cαzy4(−α)∗)2−α∑2cs2y1(−α)∗cαxwα((y2(−α)∗)2+(y3(−α)∗)2+(y4(−α)∗)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} y2n+1=y2∗+(τ−1)(S−y2n)(19)
所以其差分表示为
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} dy2n+1=dy2∗+(τ−1)(dS−dy2n)(20)
对于 d S \mathrm{d} S dS,有
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^{*} dS=(∂y1∂S)∗dy1∗+(∂y2∂S)∗dy2∗+(∂y3∂S)∗dy3∗+(∂y4∂S)∗dy4∗
以 ( ∂ S ∂ y 1 ) ∗ \left( \frac{\partial S}{\partial y_1} \right)^{*} (∂y1∂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} (∂y1∂S)∗=α∑cαxwα−α∑2cs4(y1(−α)∗)2cαxwα(cαxy2(−α)∗+cαyy3(−α)∗+cαzy4(−α)∗)2+α∑2cs2(y1(−α)∗)2cαxwα((y2(−α)∗)2+(y3(−α)∗)2+(y4(−α)∗)2)2
与预报步的做法类似,将 y i ( − α ) ∗ y_{i(-\alpha)}^{*} yi(−α)∗ 近似为 y i ∗ y_{i}^{*} yi∗,得到:
( ∂ 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} (∂y1∂S)∗=α∑cαxwα−α∑2cs4(y1∗)2cαxwα(cαxy2∗+cαyy3∗+cαzy4∗)2+α∑2cs2(y1∗)2cαxwα((y2∗)2+(y3∗)2+(y4∗)2)2
此时所有 y i ∗ y_{i}^{*} yi∗ 均可被视为常数提取出来。根据格子张量的性质,可得: ( ∂ S ∂ y 1 ) ∗ = 0 \left( \frac{\partial S}{\partial y_1} \right)^{*} = 0 (∂y1∂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
(∂y2∂S)∗=1,(∂y3∂S)∗=(∂y4∂S)∗=0
此时 d S ‾ = d y 2 ‾ ∗ \overline{\mathrm{d} S } = \overline{\mathrm{d} y_2}^{*} dS=dy2∗ ,则式 (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} dy2n+1=dy2∗=dy2n
同理,对 y 3 n + 1 y_3^{n+1} y3n+1 和 y 4 n + 1 y_4^{n+1} y4n+1 计算可得: d y 3 ‾ n + 1 = d y 3 ‾ n \overline{\mathrm{d} y_3}^{n+1} = \overline{\mathrm{d} y_3}^{n} dy3n+1=dy3n 和 d y 4 ‾ n + 1 = d y 4 ‾ n \overline{\mathrm{d} y_4}^{n+1} = \overline{\mathrm{d} y_4}^{n} dy4n+1=dy4n 。
记 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 dYn=(dy1n,dy2n,dy3n,dy4n)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 dYn+1=(dy1n+1,dy2n+1,dy3n+1,dy4n+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} dYn+1=IdYn
其中 I \mathbf{I} I 为单位矩阵。显然,特征矩阵的特征值全是 1 ,代表着该格式无条件稳定。
Chen, Z., Shu, C., Wang, Y., Yang, L. M., & Tan, D. (2017). A Simplified Lattice Boltzmann Method without Evolution of Distribution Function. Advances in Applied Mathematics and Mechanics, 9(1), 1–22. DOI:10.4208/aamm.OA-2016-0029 ↩︎
Z. Chen, C. Shu, D. Tan; Three-dimensional simplified and unconditionally stable lattice Boltzmann method for incompressible isothermal and thermal flows. Physics of Fluids. 1 May 2017; 29 (5): 053601. DOI:10.1063/1.4983339 ↩︎ ↩︎ ↩︎
Chen Z, Shu C, Tan D, Wu C. On improvements of simplified and highly stable lattice Boltzmann method: Formulations, boundary treatment, and stability analysis. Int J Numer Meth Fluids. 2018; 87: 161–179. DOI:10.1002/fld.4485 ↩︎
Chen, Z., & Shu, C. (2020). Simplified and Highly Stable Lattice Boltzmann Method. WORLD SCIENTIFIC. DOI:10.1142/12047 ↩︎
本文中所有 ∇ \nabla ∇ 均表示空间偏导 ∂ ∂ r \frac{\partial}{\partial \boldsymbol{r}} ∂r∂。 ↩︎
C. Shu, Y. Wang, C. Teo, and J. Wu. Development of lattice Boltzmann flux solver for simulation of incompressible flows. Advances In Applied Mathematics And Mechanics. 6, 436 (2014). ↩︎