最近自己也是在学相场格子Boltzmann方法，这篇笔记主要是在整理刚学到的内容。
相场采用局部守恒性 Allen-Cahn 方程描述[1,2]，其核心量为序参数 ϕ \phi ϕ
∂ ϕ ∂ t + ∇ ⋅ ( ϕ u ) = M ϕ ( ∇ 2 ϕ − ∇ ⋅ ( λ n ) ) \frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \boldsymbol{u}) =
M_{\phi} (\nabla^2\phi - \nabla\cdot(\lambda\boldsymbol{n}))
∂ t ∂ ϕ + ∇ ⋅ ( ϕ u ) = M ϕ ( ∇ 2 ϕ − ∇ ⋅ ( λ n ) )
其中
M ϕ M_{\phi} M ϕ n = ∇ ϕ / ∥ ∇ ϕ ∥ \boldsymbol{n}=\nabla\phi/\|\nabla\phi\| n = ∇ ϕ / ∥ ∇ ϕ ∥ λ = 4 ( 1 − ϕ ) ϕ / W \lambda = 4(1-\phi)\phi/W λ = 4 ( 1 − ϕ ) ϕ / W W W W
表面张力通常采用如下形式：
F s = μ ϕ ∇ ϕ \boldsymbol{F}_s = \mu_\phi \nabla \phi
F s = μ ϕ ∇ ϕ
化学势 μ ϕ \mu_\phi μ ϕ
μ ϕ = 4 β ϕ ( ϕ − 1 ) ( ϕ − 0.5 ) − k ∇ 2 ϕ . \mu_\phi = 4\beta\phi(\phi-1)(\phi-0.5) - k \nabla^2\phi.
μ ϕ = 4 β ϕ ( ϕ − 1 ) ( ϕ − 0 . 5 ) − k ∇ 2 ϕ .
k k k β \beta β W W W σ \sigma σ
k = 3 σ W / 2 , β = 12 σ / W k = 3 \sigma W /2, \quad
\beta = 12 \sigma / W
k = 3 σ W / 2 , β = 1 2 σ / W
注： Open Access的文献[1] 非常详尽地介绍了这一部分的推导，这里不做赘述。
流场的控制方程组分为两个式子，这里我们简单讨论两种不可压流体的情况。
连续性方程为：
∇ ⋅ u = 0 , \nabla \cdot \boldsymbol{u} = 0,
∇ ⋅ u = 0 ,
Navier-Stokes方程为：
∂ ( ρ u ) ∂ t + ∇ ⋅ ( ρ u u ) = − ∇ p + ∇ ⋅ [ μ ( ∇ u + ( ∇ u ) T ) ] + F \begin{aligned}
\frac{\partial (\rho\boldsymbol{u})}{\partial t} &+ \nabla \cdot (\rho \boldsymbol{u} \boldsymbol{u}) = - \nabla p\\
& + \nabla \cdot [\mu (\nabla\boldsymbol{u} + (\nabla\boldsymbol{u})^\mathrm{T})] + \mathbb{F}
\end{aligned}
∂ t ∂ ( ρ u ) + ∇ ⋅ ( ρ u u ) = − ∇ p + ∇ ⋅ [ μ ( ∇ u + ( ∇ u ) T ) ] + F
其中
ρ , u , p \rho, \boldsymbol{u}, p ρ , u , p μ , ν \mu, \nu μ , ν μ = ρ ν \mu = \rho\nu μ = ρ ν 外力 F = F s + G \mathbb{F} = \boldsymbol{F}_s + \boldsymbol{G} F = F s + G F s \boldsymbol{F}_s F s G \boldsymbol{G} G
分布函数 f i f_i f i
f i ( x + c i δ t , t + δ t ) − f i ( x , t ) = − 1 τ f [ f i ( x , t ) − f i e q ( x , t ) ] + F i ( x , t ) δ t (1) \begin{aligned}
f_i (\boldsymbol{x} &+ \boldsymbol{c}_i \delta_t, t + \delta_t) - f_i (\boldsymbol{x}, t) =\\
&-\frac{1}{\tau_f}
\left[f_i (\boldsymbol{x}, t) - f_i^{eq} (\boldsymbol{x}, t)\right] +
F_i (\boldsymbol{x}, t) \delta_t
\end{aligned}
\tag{1}
f i ( x + c i δ t , t + δ t ) − f i ( x , t ) = − τ f 1 [ f i ( x , t ) − f i e q ( x , t ) ] + F i ( x , t ) δ t ( 1 )
其中，相场中的平衡态 f i e q f_i^{eq} f i e q
f i e q = ω i ϕ ( 1 + c i ⋅ u c s 2 ) . f_i^{eq} = \omega_i \phi \left( 1 + \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} \right).
f i e q = ω i ϕ ( 1 + c s 2 c i ⋅ u ) .
式(1)里的作用力项 F i F_i F i
F i = ( 1 − 1 2 τ f ) w i c s 2 c i ⋅ [ ∂ ( ϕ u ) ∂ t + c s 2 λ n ] . F_i = \left( 1 - \frac{1}{2 \tau_f} \right)
\frac{w_i}{c_s^2}
\boldsymbol{c}_i \cdot \left[
\frac{\partial (\phi \boldsymbol{u})}{\partial t} + c_s^2 \lambda \boldsymbol{n}
\right].
F i = ( 1 − 2 τ f 1 ) c s 2 w i c i ⋅ [ ∂ t ∂ ( ϕ u ) + c s 2 λ n ] .
松弛时间 τ f \tau_f τ f M ϕ M_{\phi} M ϕ
M ϕ = ( τ f − 1 2 ) c s 2 δ t M_{\phi} = (\tau_f - \frac{1}{2}) c_s^2 \delta_t
M ϕ = ( τ f − 2 1 ) c s 2 δ t
序参数 (order parameter) ϕ \phi ϕ ρ \rho ρ ρ l \rho_l ρ l ρ g \rho_g ρ g
ρ = ϕ ρ l + ( 1 − ϕ ) ρ g . \rho = \phi \rho_l + (1 - \phi) \rho_g .
ρ = ϕ ρ l + ( 1 − ϕ ) ρ g .
而式(1)的模型中， ϕ \phi ϕ ϕ = ∑ i f i \displaystyle\phi=\sum_{i}f_{i} ϕ = i ∑ f i
混合流体的运动粘度 ν \nu ν ϕ \phi ϕ
（1）阶跃函数 ν ( ϕ ) = { ν g , ϕ < 0.5 ν l , ϕ ≥ 0.5 \displaystyle\nu(\phi)=\begin{cases}\nu_g,&\phi < 0.5\\ \nu_l,&\phi \ge 0.5\end{cases} ν ( ϕ ) = { ν g , ν l , ϕ < 0 . 5 ϕ ≥ 0 . 5 ν g \nu_g ν g ν l \nu_l ν l
（2）线性函数 ν ( ϕ ) = ϕ ν l + ( 1 − ϕ ) ν g \displaystyle\nu(\phi)=\phi\nu_{l}+(1-\phi)\nu_{g} ν ( ϕ ) = ϕ ν l + ( 1 − ϕ ) ν g
（3）倒数的线性形式：1 ν ( ϕ ) = ϕ ν l + 1 − ϕ ν g \displaystyle\frac{1}{\nu(\phi)}=\frac{\phi}{\nu_{l}}+\frac{1-\phi}{\nu_{g}} ν ( ϕ ) 1 = ν l ϕ + ν g 1 − ϕ
分布函数 g i g_i g i
g i ( x + c i δ t , t + δ t ) − g i ( x , t ) = − 1 τ g [ g i ( x , t ) − g i e q ( x , t ) ] + G i ( x , t ) ⋅ δ t (2) \begin{aligned}
g_i (\boldsymbol{x} &+ \boldsymbol{c}_i \delta_t, t + \delta_t) - g_i (\boldsymbol{x}, t) =\\
&-\frac{1}{\tau_g}
\left[g_i (\boldsymbol{x}, t) - g_i^{eq} (\boldsymbol{x}, t)\right] +
G_i (\boldsymbol{x}, t) \cdot \delta_t
\end{aligned}
\tag{2}
g i ( x + c i δ t , t + δ t ) − g i ( x , t ) = − τ g 1 [ g i ( x , t ) − g i e q ( x , t ) ] + G i ( x , t ) ⋅ δ t ( 2 )
局部平衡态 g i e q g_i^{eq} g i e q
g i e q = { p c s 2 ( ω i − 1 ) + ρ s i ( u ) , i = 0 p c s 2 ω i + ρ s i ( u ) , i ≠ 0 (3) g_i^{eq} =
\begin{cases}
\frac{p}{c_s^2} (\omega_i - 1) + \rho s_{i}(\boldsymbol{u}), & i = 0\\
\frac{p}{c_s^2} \omega_i + \rho s_{i}(\boldsymbol{u}), & i \neq 0\\
\end{cases}
\tag{3}
g i e q = { c s 2 p ( ω i − 1 ) + ρ s i ( u ) , c s 2 p ω i + ρ s i ( u ) , i = 0 i = 0 ( 3 )
其中
s i ( u ) = [ c i ⋅ u c s 2 + ( c i ⋅ u ) 2 2 c s 4 − u ⋅ u 2 c s 2 ] ⋅ ω i s_{i}(\boldsymbol{u}) = \left[
\frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} +
\frac{(\boldsymbol{c}_i \cdot \boldsymbol{u})^2}{2 c_s^4} -
\frac{\boldsymbol{u} \cdot \boldsymbol{u}}{2 c_s^2}
\right] \cdot \omega_i
s i ( u ) = [ c s 2 c i ⋅ u + 2 c s 4 ( c i ⋅ u ) 2 − 2 c s 2 u ⋅ u ] ⋅ ω i
源项为
G i ( x , t ) = ( 1 − 1 2 τ g ) ω i [ u ⋅ ∇ ρ + c i ⋅ F c s 2 + u ∇ ρ : ( c i c i − c s 2 I ) c s 2 ] = ( 1 − 1 2 τ g ) ω i [ c i ⋅ F c s 2 + ( ρ l − ρ g ) u ∇ ϕ : ( c i c i ) c s 2 ] (4) \begin{aligned}
G_i (\boldsymbol{x}, t) &=
\left(1 - \frac{1}{2 \tau_g}\right) \omega_i
\left[
\boldsymbol{u} \cdot \nabla \rho +
\frac{\boldsymbol{c}_i \cdot \mathbb{F}}{c_s^2} +
\frac{\boldsymbol{u} \nabla \rho: (\boldsymbol{c}_i \boldsymbol{c}_i - c_s^2 \boldsymbol{I})}{c_s^2}
\right]\\
&=
\left(1 - \frac{1}{2 \tau_g}\right) \omega_i
\left[
\frac{\boldsymbol{c}_i \cdot \mathbb{F}}{c_s^2} +
\frac{(\rho_l - \rho_g) \boldsymbol{u} \nabla\phi: (\boldsymbol{c}_i \boldsymbol{c}_i)}{c_s^2}
\right]
\end{aligned}
\tag{4}
G i ( x , t ) = ( 1 − 2 τ g 1 ) ω i [ u ⋅ ∇ ρ + c s 2 c i ⋅ F + c s 2 u ∇ ρ : ( c i c i − c s 2 I ) ] = ( 1 − 2 τ g 1 ) ω i [ c s 2 c i ⋅ F + c s 2 ( ρ l − ρ g ) u ∇ ϕ : ( c i c i ) ] ( 4 )
式(4)用到了 ∇ ρ = ( ρ l − ρ g ) ∇ ϕ \nabla \rho = (\rho_l - \rho_g) \nabla\phi ∇ ρ = ( ρ l − ρ g ) ∇ ϕ
所以，宏观量的计算为：
ρ u = δ t F 2 + ∑ i c i g i , p = c s 2 1 − ω 0 [ ∑ i ≠ 0 g i + δ t 2 u ⋅ ∇ ρ + ρ s 0 ( u ) ] = c s 2 1 − ω 0 [ ∑ i ≠ 0 g i + δ t 2 ( ρ l − ρ g ) u ⋅ ∇ ϕ + ρ s 0 ( u ) ] . (5) \begin{aligned}
\rho \boldsymbol{u} &= \frac{\delta_t \mathbb{F}}{2} + \sum_{i} \boldsymbol{c}_i g_i ,\\
p &= \frac{c_s^2}{1 - \omega_0} \left[ \sum_{i \neq 0}{g_i} + \frac{\delta_t}{2} \boldsymbol{u} \cdot \nabla \rho + \rho s_{0}(\boldsymbol{u}) \right]\\
&= \frac{c_s^2}{1 - \omega_0} \left[ \sum_{i \neq 0}{g_i} + \frac{\delta_t}{2} (\rho_l - \rho_g)\boldsymbol{u} \cdot \nabla\phi + \rho s_{0}(\boldsymbol{u}) \right].
\end{aligned}
\tag{5}
ρ u p = 2 δ t F + i ∑ c i g i , = 1 − ω 0 c s 2 ⎣ ⎢ ⎡ i = 0 ∑ g i + 2 δ t u ⋅ ∇ ρ + ρ s 0 ( u ) ⎦ ⎥ ⎤ = 1 − ω 0 c s 2 ⎣ ⎢ ⎡ i = 0 ∑ g i + 2 δ t ( ρ l − ρ g ) u ⋅ ∇ ϕ + ρ s 0 ( u ) ⎦ ⎥ ⎤ . ( 5 )
记 A = − c s 2 τ f δ t \displaystyle\mathcal{A}=-c_s^2 \tau_f \delta_t A = − c s 2 τ f δ t B = M ϕ λ δ t \displaystyle\mathcal{B} = M_{\phi} \lambda \delta_t B = M ϕ λ δ t C = δ t 2 ϕ u + ∑ i c i ( f i − f i e q ) \displaystyle\mathcal{C} = \frac{\delta_t}{2}\phi\boldsymbol{u} + \sum_{i}{\boldsymbol{c}_i (f_i - f_i^{eq})} C = 2 δ t ϕ u + i ∑ c i ( f i − f i e q )
∥ ∇ ϕ ∥ = − ∥ C ∥ − B A , ∇ ϕ = C A + B / ∥ ∇ ϕ ∥ . (6) \|\nabla\phi\| = \frac{-\|\mathcal{C}\|-\mathcal{B}}{\mathcal{A}} ,
\quad
\nabla\phi = \frac{\mathcal{C}}{\mathcal{A} + \mathcal{B}/\|\nabla\phi\|}.
\tag{6}
∥ ∇ ϕ ∥ = A − ∥ C ∥ − B , ∇ ϕ = A + B / ∥ ∇ ϕ ∥ C . ( 6 )
化学势 μ ϕ \mu_\phi μ ϕ ∇ 2 ϕ \nabla^2\phi ∇ 2 ϕ ∇ ϕ \nabla\phi ∇ ϕ
当然，也有完全使用差分的计算公式
{ ∇ ϕ ( x ) = ∑ i ≠ 0 ω i c i ϕ ( x + c i δ t ) / c s 2 δ t , ∇ 2 ϕ ( x ) = ∑ i ≠ 0 2 ω i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] / c s 2 δ t 2 \begin{cases}
\nabla\phi(\boldsymbol{x}) &= \sum_{i \neq 0} {\omega_i \boldsymbol{c}_i \phi(\boldsymbol{x} + \boldsymbol{c}_i \delta_t)} / {c_s^2 \delta_t} ,\\
\nabla^2\phi(\boldsymbol{x}) &= \sum_{i \neq 0} {2\omega_i [\phi(\boldsymbol{x} + \boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x})]} / {c_s^2 \delta_t^2}
\end{cases}
{ ∇ ϕ ( x ) ∇ 2 ϕ ( x ) = ∑ i = 0 ω i c i ϕ ( x + c i δ t ) / c s 2 δ t , = ∑ i = 0 2 ω i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] / c s 2 δ t 2
值得一提的是，Tim Reis [3] 提出了一种不需要计算 ∇ 2 ϕ \nabla^2\phi ∇ 2 ϕ g i g_i g i
定义表面张力为 F s = ∇ ⋅ T \boldsymbol{F}_s = \nabla\cdot\mathrm{T} F s = ∇ ⋅ T T \mathrm{T} T
T = σ ( I − n n ) δ S \mathrm{T}=\sigma (\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n}) \delta_{\mathrm{S}}
T = σ ( I − n n ) δ S
而 δ S ≃ ∥ ∇ ϕ ∥ \delta_{\mathrm{S}}\simeq\|\nabla\phi\| δ S ≃ ∥ ∇ ϕ ∥
Tim Reis 基于多尺度展开结果，把 T \mathrm{T} T
g i e q = ω i [ p c s 2 + ρ c i ⋅ u c s 2 + 1 2 c s 4 ( ρ u u + T ) : ( c i c i − c s 2 I ) ] , g_{i}^{eq}=\omega_i \left[ \frac{p}{c_s^2} + \rho \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} + \frac{1}{2c_s^4}(\rho\boldsymbol{u}\boldsymbol{u}+\mathrm{T}):(\boldsymbol{c}_i \boldsymbol{c}_i - c_s^2 \boldsymbol{I}) \right],
g i e q = ω i [ c s 2 p + ρ c s 2 c i ⋅ u + 2 c s 4 1 ( ρ u u + T ) : ( c i c i − c s 2 I ) ] ,
这使得化学势 μ ϕ \mu_\phi μ ϕ G \boldsymbol{G} G G i G_i G i
G i = ( 1 − 1 2 τ g ) [ R i ( G ) + S i ] G_i = \left(1 - \frac{1}{2 \tau_g}\right) [R_i(\boldsymbol{G}) + S_i]
G i = ( 1 − 2 τ g 1 ) [ R i ( G ) + S i ]
其中 R i ( G ) R_i(\boldsymbol{G}) R i ( G )
R i = ω i ( c i − u c s 2 + c i ⋅ u c s 4 c i ) ⋅ G , R_i = \omega_i \left(
\frac{\boldsymbol{c}_i - \boldsymbol{u}}{c_s^2} +
\frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^4} \boldsymbol{c}_i
\right) \cdot \boldsymbol{G},
R i = ω i ( c s 2 c i − u + c s 4 c i ⋅ u c i ) ⋅ G ,
而 S i S_i S i
{ S i = ( Γ i − ω i ) ( c i − u ) ⋅ ∇ ρ Γ i = ( 1 + c i − u c s 2 + u u : ( c i c i − c s 2 I ) 2 c s 4 ) ω i \begin{cases}
S_i = (\Gamma_i - \omega_i) (\boldsymbol{c}_i - \boldsymbol{u}) \cdot \nabla\rho \\
\Gamma_i =\left(
1 + \frac{\boldsymbol{c}_i - \boldsymbol{u}}{c_s^2} +
\frac{\boldsymbol{u}\boldsymbol{u}:(\boldsymbol{c}_i \boldsymbol{c}_i - c_s^2 \boldsymbol{I})}{2c_s^4}
\right) \omega_i
\end{cases}
{ S i = ( Γ i − ω i ) ( c i − u ) ⋅ ∇ ρ Γ i = ( 1 + c s 2 c i − u + 2 c s 4 u u : ( c i c i − c s 2 I ) ) ω i
宏观量的计算也要变成
{ ∑ i g i = p c s 2 − δ t 2 ( u ⋅ ∇ ) ρ ∑ i g i c i = ρ u + δ t 2 G . \begin{cases}
\sum_i g_i = \frac{p}{c_s^2} - \frac{\delta_t}{2} (\boldsymbol{u} \cdot \nabla)\rho \\
\sum_i g_i \boldsymbol{c}_i = \rho \boldsymbol{u} + \frac{\delta_t}{2} \boldsymbol{G}
\end{cases}.
{ ∑ i g i = c s 2 p − 2 δ t ( u ⋅ ∇ ) ρ ∑ i g i c i = ρ u + 2 δ t G .
带有 '※' 标记的文章意味着可通过免费渠道获得，如：开放获取 (Open Access)；Research gate；arXiv预印本；Sci hub；作者主页等。
[1] ※ WANG H, YUAN X, LIANG H, et al . A brief review of the phase-field-based lattice Boltzmann method for multiphase flows [J]. Capillarity, 2019, 2(3): 33-52. DOI:10.26804/capi.2019.03.01 .
[2] ※ LIANG H, XU J, CHEN J, et al . Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows [J]. Physical Review E, 2018, 97(3): 033309. DOI:10.1103/PhysRevE.97.033309 .
[3] ※ REIS T. A lattice Boltzmann formulation of the one-fluid model for multiphase flow [J]. Journal of Computational Physics, 2022, 453: 110962. DOI:10.1016/j.jcp.2022.110962 .
