在 “【笔记】格子Boltzmann方法的相场模型” 中，我把LBM里的相场模型简单记了一下。笔记里有提到，序参数 ϕ \phi ϕ
{ ∇ ϕ ( x ) = ∑ i ≠ 0 w i c i ϕ ( x + c i δ t ) c s 2 δ t , ∇ 2 ϕ ( x ) = ∑ i ≠ 0 2 w i c s 2 δ t 2 [ ϕ ( x + c i δ t ) − ϕ ( x ) ] \displaystyle
\begin{cases}
\nabla\phi(\boldsymbol{x}) =& \sum\limits_{i \neq 0} \dfrac{w_i \boldsymbol{c}_i \phi(\boldsymbol{x} + \boldsymbol{c}_i \delta_t)}{c_s^2 \delta_t} ,\\
\nabla^2\phi(\boldsymbol{x}) =& \sum\limits_{i \neq 0} \dfrac{2w_i}{c_s^2 \delta_t^2} [\phi(\boldsymbol{x} + \boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x})]
\end{cases}
⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ ∇ ϕ ( x ) = ∇ 2 ϕ ( x ) = i = 0 ∑ c s 2 δ t w i c i ϕ ( x + c i δ t ) , i = 0 ∑ c s 2 δ t 2 2 w i [ ϕ ( x + c i δ t ) − ϕ ( x ) ]
其中 i i i i i i c i \boldsymbol{c}_i c i w i w_i w i c i \boldsymbol{c}_i c i
对于点 x \boldsymbol{x} x x + c i δ t \boldsymbol{x}+\boldsymbol{c}_i \delta_t x + c i δ t ϕ \phi ϕ x \boldsymbol{x} x
ϕ ( x + c i δ t ) = ϕ ( x ) + ( c i δ t ) ⋅ ∇ ϕ ( x ) + δ t 2 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) (1) \phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) = \phi(\boldsymbol{x}) + (\boldsymbol{c}_i \delta_t) \cdot \nabla\phi(\boldsymbol{x}) + \frac{\delta_t^2}{2} (\boldsymbol{c}_i \cdot \nabla)^2 \phi(\boldsymbol{x}) + O(\delta_t^3)
\tag{1}
ϕ ( x + c i δ t ) = ϕ ( x ) + ( c i δ t ) ⋅ ∇ ϕ ( x ) + 2 δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) ( 1 )
式(1)两边同时乘上 w i c i w_i \boldsymbol{c}_i w i c i
∑ i w i c i ϕ ( x + c i δ t ) = ∑ i w i c i ϕ ( x ) + ∑ i w i c i [ ( c i δ t ) ⋅ ∇ ϕ ( x ) ] + ∑ i w i c i [ δ t 2 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) ] + O ( δ t 3 ) (2) \begin{aligned}
\sum_i w_i \boldsymbol{c}_i \phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) &=
\sum_i w_i \boldsymbol{c}_i \phi(\boldsymbol{x}) \\
&\quad + \sum_i w_i \boldsymbol{c}_i [(\boldsymbol{c}_i \delta_t) \cdot \nabla\phi(\boldsymbol{x})] \\
&\quad + \sum_i w_i \boldsymbol{c}_i \left[
\frac{\delta_t^2}{2} (\boldsymbol{c}_i \cdot \nabla)^2 \phi(\boldsymbol{x})
\right] + O(\delta_t^3)
\end{aligned}
\tag{2}
i ∑ w i c i ϕ ( x + c i δ t ) = i ∑ w i c i ϕ ( x ) + i ∑ w i c i [ ( c i δ t ) ⋅ ∇ ϕ ( x ) ] + i ∑ w i c i [ 2 δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) ] + O ( δ t 3 ) ( 2 )
下面分析式(2)右侧的三个项。
（1）对第1项，ϕ ( x ) \phi(\boldsymbol{x}) ϕ ( x )
∑ i w i c i ϕ ( x ) = ϕ ( x ) ⋅ ∑ i ( w i c i ) \sum_i w_i \boldsymbol{c}_i \phi(\boldsymbol{x}) = \phi(\boldsymbol{x}) \cdot \sum_i (w_i \boldsymbol{c}_i )
i ∑ w i c i ϕ ( x ) = ϕ ( x ) ⋅ i ∑ ( w i c i )
由于 ∑ i w i c i = 0 ⃗ \sum_i w_i \boldsymbol{c}_i = \vec{0} ∑ i w i c i = 0
（2）对第2项，
∑ i w i c i [ ( c i δ t ) ⋅ ∇ ϕ ( x ) ] = δ t ∇ ϕ ( x ) ⋅ ∑ i ( w i c i c i ) = δ t c s 2 ∇ ϕ ( x ) \sum_i w_i \boldsymbol{c}_i [(\boldsymbol{c}_i \delta_t) \cdot \nabla\phi(\boldsymbol{x})] = \delta_t \nabla\phi(\boldsymbol{x}) \cdot \sum_i (w_i \boldsymbol{c}_i \boldsymbol{c}_i) = \delta_t c_s^2 \nabla\phi(\boldsymbol{x})
i ∑ w i c i [ ( c i δ t ) ⋅ ∇ ϕ ( x ) ] = δ t ∇ ϕ ( x ) ⋅ i ∑ ( w i c i c i ) = δ t c s 2 ∇ ϕ ( x )
其中 ∑ i ( w i c i c i ) = c s 2 [ I ] \sum_i (w_i \boldsymbol{c}_i \boldsymbol{c}_i) = c_s^2 [\boldsymbol{I}] ∑ i ( w i c i c i ) = c s 2 [ I ]
这个关系可以由常规 DnQb 模型中的 ∑ i c i c i f i e q = ρ ( u u + c s 2 [ I ] ) \sum_i \boldsymbol{c}_i \boldsymbol{c}_i f_i^{eq} = \rho (\boldsymbol{u}\boldsymbol{u} + c_s^2 [\boldsymbol{I}]) ∑ i c i c i f i e q = ρ ( u u + c s 2 [ I ] ) u = 0 ⃗ \boldsymbol{u}=\vec{0} u = 0
（3）对第3项，为简化考虑，这里推导其中一个分量，其形式为：
δ t 2 2 ∑ i w i c i γ [ c i α c i β ∂ 2 ϕ ( x ) ∂ x α ∂ x β ] = δ t 2 2 ∂ 2 ϕ ( x ) ∂ x α ∂ x β ∑ i w i c i α c i β c i γ \frac{\delta_t^2}{2} \sum_i w_i c_{i\gamma} \left[
c_{i\alpha} c_{i\beta} \frac{\partial^2 \phi(\boldsymbol{x})}{\partial x_{\alpha} \partial x_{\beta}}
\right] =
\frac{\delta_t^2}{2} \frac{\partial^2 \phi(\boldsymbol{x})}{\partial x_{\alpha} \partial x_{\beta}} \sum_i w_i c_{i\alpha} c_{i\beta} c_{i\gamma}
2 δ t 2 i ∑ w i c i γ [ c i α c i β ∂ x α ∂ x β ∂ 2 ϕ ( x ) ] = 2 δ t 2 ∂ x α ∂ x β ∂ 2 ϕ ( x ) i ∑ w i c i α c i β c i γ
由于在常规 DnQb 模型中，
∑ i c i α c i β c i γ f i e q = c s 2 ρ [ u α δ β γ + u β δ α γ + u γ δ α β ] \sum_i c_{i\alpha} c_{i\beta} c_{i\gamma} f_i^{eq} = c_s^2 \rho [u_{\alpha} \delta_{\beta\gamma}+ u_{\beta} \delta_{\alpha\gamma} + u_{\gamma} \delta_{\alpha\beta}]
i ∑ c i α c i β c i γ f i e q = c s 2 ρ [ u α δ β γ + u β δ α γ + u γ δ α β ]
这里的 δ \delta δ u = 0 ⃗ \boldsymbol{u}=\vec{0} u = 0
∑ i w i c i α c i β c i γ = c s 2 [ u α δ β γ + u β δ α γ + u γ δ α β ] = 0 \sum_i w_i c_{i\alpha} c_{i\beta} c_{i\gamma} = c_s^2 [u_{\alpha} \delta_{\beta\gamma}+ u_{\beta} \delta_{\alpha\gamma} + u_{\gamma} \delta_{\alpha\beta}] = 0
i ∑ w i c i α c i β c i γ = c s 2 [ u α δ β γ + u β δ α γ + u γ δ α β ] = 0
所以式(2)被化简为
∑ i w i c i ϕ ( x + c i δ t ) = δ t c s 2 ∇ ϕ ( x ) + O ( δ t 3 ) \sum_i w_i \boldsymbol{c}_i \phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) =
\delta_t c_s^2 \nabla\phi(\boldsymbol{x}) + O(\delta_t^3)
i ∑ w i c i ϕ ( x + c i δ t ) = δ t c s 2 ∇ ϕ ( x ) + O ( δ t 3 )
即：
∇ ϕ ( x ) = 1 δ t c s 2 ∑ i w i c i ϕ ( x + c i δ t ) + O ( δ t 2 ) (3) \nabla\phi(\boldsymbol{x}) = \frac{1}{\delta_t c_s^2} \sum_i w_i \boldsymbol{c}_i \phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) + O(\delta_t^2)
\tag{3}
∇ ϕ ( x ) = δ t c s 2 1 i ∑ w i c i ϕ ( x + c i δ t ) + O ( δ t 2 ) ( 3 )
将式(1)改写成
ϕ ( x + c i δ t ) − ϕ ( x ) = ( c i δ t ) ⋅ ∇ ϕ ( x ) + δ t 2 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) (4) \phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x}) = (\boldsymbol{c}_i \delta_t) \cdot \nabla\phi(\boldsymbol{x}) + \frac{\delta_t^2}{2} (\boldsymbol{c}_i \cdot \nabla)^2 \phi(\boldsymbol{x}) + O(\delta_t^3)
\tag{4}
ϕ ( x + c i δ t ) − ϕ ( x ) = ( c i δ t ) ⋅ ∇ ϕ ( x ) + 2 δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) ( 4 )
式(4)两边同时乘上 w i w_i w i
∑ i w i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] = ∑ i w i [ ( c i δ t ) ⋅ ∇ ϕ ( x ) ] + ∑ i w i [ δ t 2 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) ] + O ( δ t 3 ) (5) \begin{aligned}
\sum_i w_i [\phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x})] &= \sum_i w_i [(\boldsymbol{c}_i \delta_t) \cdot \nabla\phi(\boldsymbol{x})] \\
&\quad + \sum_i w_i \left[
\frac{\delta_t^2}{2} (\boldsymbol{c}_i \cdot \nabla)^2 \phi(\boldsymbol{x})
\right] + O(\delta_t^3)
\end{aligned}
\tag{5}
i ∑ w i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] = i ∑ w i [ ( c i δ t ) ⋅ ∇ ϕ ( x ) ] + i ∑ w i [ 2 δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) ] + O ( δ t 3 ) ( 5 )
参照上一节的推导，易得式(5)等式右侧的两项分别满足
∑ i w i c i α δ t ∇ α ϕ ( x ) = δ t ∇ α ϕ ( x ) ∑ i w i c i α = 0 \sum_i w_i c_{i\alpha} \delta_t \nabla_{\alpha}\phi(\boldsymbol{x}) =
\delta_t \nabla_{\alpha}\phi(\boldsymbol{x}) \sum_i w_i c_{i\alpha} = 0
i ∑ w i c i α δ t ∇ α ϕ ( x ) = δ t ∇ α ϕ ( x ) i ∑ w i c i α = 0
以及
∑ i w i δ t 2 2 c i α c i β ⋅ ∇ α β ϕ ( x ) = δ t 2 2 ∇ α β ϕ ( x ) ∑ i w i c i α c i β = δ t 2 c s 2 δ α β 2 ∇ α β ϕ ( x ) \sum_i w_i \frac{\delta_t^2}{2} c_{i\alpha} c_{i\beta} \cdot \nabla_{\alpha\beta}\phi(\boldsymbol{x}) = \frac{\delta_t^2}{2} \nabla_{\alpha\beta}\phi(\boldsymbol{x}) \sum_i w_i c_{i\alpha} c_{i\beta} = \frac{\delta_t^2 c_s^2 \delta_{\alpha\beta}}{2} \nabla_{\alpha\beta}\phi(\boldsymbol{x})
i ∑ w i 2 δ t 2 c i α c i β ⋅ ∇ α β ϕ ( x ) = 2 δ t 2 ∇ α β ϕ ( x ) i ∑ w i c i α c i β = 2 δ t 2 c s 2 δ α β ∇ α β ϕ ( x )
其中 ∇ α = ∂ ∂ x α \nabla_{\alpha} = \frac{\partial}{\partial x_{\alpha}} ∇ α = ∂ x α ∂ ∇ α β = ∂ 2 ∂ x α ∂ x β \nabla_{\alpha\beta} = \frac{\partial^2}{\partial x_{\alpha}\partial x_{\beta}} ∇ α β = ∂ x α ∂ x β ∂ 2 δ α β \delta_{\alpha\beta} δ α β
所以式(5)表示为
∑ i w i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] = δ t 2 c s 2 2 ∇ 2 ϕ ( x ) + O ( δ t 3 ) \sum_i w_i [\phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x})] = \frac{\delta_t^2 c_s^2}{2} \nabla^2\phi(\boldsymbol{x}) + O(\delta_t^3)
i ∑ w i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] = 2 δ t 2 c s 2 ∇ 2 ϕ ( x ) + O ( δ t 3 )
即：
∇ 2 ϕ ( x ) = 2 δ t 2 c s 2 ∑ i w i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] + O ( δ t 2 ) (6) \nabla^2\phi(\boldsymbol{x}) = \frac{2}{\delta_t^2 c_s^2} \sum_i w_i [\phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x})] + O(\delta_t^2)
\tag{6}
∇ 2 ϕ ( x ) = δ t 2 c s 2 2 i ∑ w i [ ϕ ( x + c i δ t ) − ϕ ( x ) ] + O ( δ t 2 ) ( 6 )
假设我们采用中心差分形式，需要对 x − c i δ t \boldsymbol{x}-\boldsymbol{c}_i \delta_t x − c i δ t
ϕ ( x − c i δ t ) = ϕ ( x ) − ( c i δ t ) ⋅ ∇ ϕ ( x ) + δ t 2 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) (7) \phi(\boldsymbol{x}-\boldsymbol{c}_i \delta_t) = \phi(\boldsymbol{x}) - (\boldsymbol{c}_i \delta_t) \cdot \nabla\phi(\boldsymbol{x}) + \frac{\delta_t^2}{2} (\boldsymbol{c}_i \cdot \nabla)^2 \phi(\boldsymbol{x}) + O(\delta_t^3)
\tag{7}
ϕ ( x − c i δ t ) = ϕ ( x ) − ( c i δ t ) ⋅ ∇ ϕ ( x ) + 2 δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) ( 7 )
（A）将式(1)与式(7)相减，偶数阶项（零阶项 ϕ ( x ) \phi(\boldsymbol{x}) ϕ ( x )
ϕ ( x + c i δ t ) − ϕ ( x − c i δ t ) = 2 ( c i δ t ) ⋅ ∇ ϕ ( x ) + O ( δ t 3 ) \phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x}-\boldsymbol{c}_i \delta_t) = 2(\boldsymbol{c}_i \delta_t) \cdot \nabla\phi(\boldsymbol{x}) + O(\delta_t^3)
ϕ ( x + c i δ t ) − ϕ ( x − c i δ t ) = 2 ( c i δ t ) ⋅ ∇ ϕ ( x ) + O ( δ t 3 )
两边同时乘上 w i c i w_i \boldsymbol{c}_i w i c i
∑ i w i c i [ ϕ ( x + c i δ t ) − ϕ ( x − c i δ t ) ] = 2 δ t ∇ ϕ ( x ) ⋅ ∑ i ( w i c i c i ) + O ( δ t 3 ) \sum_i w_i \boldsymbol{c}_i [\phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x}-\boldsymbol{c}_i \delta_t)] = 2\delta_t \nabla\phi(\boldsymbol{x}) \cdot \sum_i (w_i \boldsymbol{c}_i \boldsymbol{c}_i) + O(\delta_t^3)
i ∑ w i c i [ ϕ ( x + c i δ t ) − ϕ ( x − c i δ t ) ] = 2 δ t ∇ ϕ ( x ) ⋅ i ∑ ( w i c i c i ) + O ( δ t 3 )
利用前面提到的二阶矩性质 ∑ i w i c i c i = c s 2 I \sum_i w_i \boldsymbol{c}_i \boldsymbol{c}_i = c_s^2 \boldsymbol{I} ∑ i w i c i c i = c s 2 I 2 δ t c s 2 ∇ ϕ ( x ) 2\delta_t c_s^2 \nabla\phi(\boldsymbol{x}) 2 δ t c s 2 ∇ ϕ ( x )
∇ ϕ ( x ) = 1 2 δ t c s 2 ∑ i w i c i [ ϕ ( x + c i δ t ) − ϕ ( x − c i δ t ) ] + O ( δ t 2 ) (8) \nabla\phi(\boldsymbol{x}) = \frac{1}{2\delta_t c_s^2} \sum_i w_i \boldsymbol{c}_i [\phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - \phi(\boldsymbol{x}-\boldsymbol{c}_i \delta_t)] + O(\delta_t^2)
\tag{8}
∇ ϕ ( x ) = 2 δ t c s 2 1 i ∑ w i c i [ ϕ ( x + c i δ t ) − ϕ ( x − c i δ t ) ] + O ( δ t 2 ) ( 8 )
式(8)就是 Lee 等在文献[1]的式(64)的格式。
（B）将式(1)与式(7)相加，奇数阶项（一阶项和三阶项）会被消去，得到：
ϕ ( x + c i δ t ) − 2 ϕ ( x ) + ϕ ( x − c i δ t ) = δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) \phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - 2 \phi(\boldsymbol{x}) + \phi(\boldsymbol{x}-\boldsymbol{c}_i \delta_t) = \delta_t^2 (\boldsymbol{c}_i \cdot \nabla)^2 \phi(\boldsymbol{x}) + O(\delta_t^3)
ϕ ( x + c i δ t ) − 2 ϕ ( x ) + ϕ ( x − c i δ t ) = δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 )
两边同时乘上 w i w_i w i
∑ i w i [ ϕ ( x + c i δ t ) − 2 ϕ ( x ) + ϕ ( x − c i δ t ) ] = ∑ i δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 ) \sum_i w_i [\phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - 2 \phi(\boldsymbol{x}) + \phi(\boldsymbol{x}-\boldsymbol{c}_i \delta_t)] = \sum_i \delta_t^2 (\boldsymbol{c}_i \cdot \nabla)^2 \phi(\boldsymbol{x}) + O(\delta_t^3)
i ∑ w i [ ϕ ( x + c i δ t ) − 2 ϕ ( x ) + ϕ ( x − c i δ t ) ] = i ∑ δ t 2 ( c i ⋅ ∇ ) 2 ϕ ( x ) + O ( δ t 3 )
同理可得：
∇ 2 ϕ ( x ) = 1 c s 2 δ t 2 ∑ i w i [ ϕ ( x + c i δ t ) − 2 ϕ ( x ) + ϕ ( x − c i δ t ) ] + O ( δ t 2 ) (9) \nabla^2 \phi(\boldsymbol{x}) = \frac{1}{c_s^2 \delta_t^2} \sum_i w_i [\phi(\boldsymbol{x}+\boldsymbol{c}_i \delta_t) - 2 \phi(\boldsymbol{x}) + \phi(\boldsymbol{x}-\boldsymbol{c}_i \delta_t)] + O(\delta_t^2)
\tag{9}
∇ 2 ϕ ( x ) = c s 2 δ t 2 1 i ∑ w i [ ϕ ( x + c i δ t ) − 2 ϕ ( x ) + ϕ ( x − c i δ t ) ] + O ( δ t 2 ) ( 9 )
式(9)就是 Lee 等在文献[1]的式(65)的格式。
[1] Lee T，Lin C-L. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio [J]. Journal of Computational Physics ，2005，206（1）：16-47. DOI:10.1016/j.jcp.2004.12.001.