Ives's Pages

发表于 分类于

双松弛时间（Two-relaxation-time）格子 Boltzmann 方法是多松弛时间（Multi-Relaxation-Time）LBM 的一种特殊形式。它仅具有两个松弛时间，分为对称和非对称两部分，同时保持了与 LBGK 模型一致的简单算法和计算效率。

基本方程

对于 tt 时刻的点 x\boldsymbol{x} 沿 ci\boldsymbol{c}_i 方向的分布函数 fi(x,t)f_i(\boldsymbol{x}, t)，其在 Δt\Delta t 时间步长内的演化为：

fi(x+ciΔt,t+Δt)−fi(x,t)=−1τ+[fi+−fi(eq)+]−1τ−[fi−−fi(eq)−]+FiΔt(1)\begin{aligned} f_i(\boldsymbol{x}+\boldsymbol{c}_i \Delta t, t + \Delta t) - f_i(\boldsymbol{x}, t) &= -\frac{1}{\tau^{+}} \left[ f_i^{+} - f_i^{(\mathrm{eq})+} \right] \\ &\quad -\frac{1}{\tau^{-}} \left[ f_i^{-} - f_i^{(\mathrm{eq})-} \right] \\ &\quad + \mathbb{F}_i \Delta t \end{aligned} \tag{1}

其中：

  • τ+\tau^{+} 和 τ−\tau^{-} 分别为两部分的松弛时间（均为正数），其数值由流体运动粘度 ν\nu 和参数 Λ=(τ+−12)(τ−−12)\Lambda = (\tau^{+} - \frac{1}{2}) (\tau^{-} - \frac{1}{2}) 共同决定。d’Humières 等[1]提出 Λ\Lambda 控制了 TRT 方程的稳态场。并且得出 Λ=1/4\Lambda=1/4 可使 TRT 的运算保持稳定。

τ+=νcs2+12,τ−=2Λ2τ+−1+12\tau^{+} = \frac{\nu}{c_s^2} + \frac{1}{2} ,\quad \tau^{-} = \frac{2 \Lambda}{2 \tau^{+} - 1} + \frac{1}{2}

  • f−i(x,t)f_{-i}(\boldsymbol{x}, t) 和 f−i(eq)(x,t)f_{-i}^{(\mathrm{eq})}(\boldsymbol{x}, t) 为 c−i=−ci\boldsymbol{c}_{-i}=-\boldsymbol{c}_{i} 方向的分布函数和平衡态分布。

  • fi+f_i^{+} 和 fi−f_i^{-} 分别为分布函数的对称部分和非对称部分，即

fi+=12(fi(x,t)+f−i(x,t)),fi−=12(fi(x,t)−f−i(x,t))f_i^{+} = \frac{1}{2}\left( f_i(\boldsymbol{x}, t) + f_{-i}(\boldsymbol{x}, t) \right) ,\quad f_i^{-} = \frac{1}{2}\left( f_i(\boldsymbol{x}, t) - f_{-i}(\boldsymbol{x}, t) \right)

  • 同理， fi(eq)+f_i^{(\mathrm{eq})+} 和 fi(eq)−f_i^{(\mathrm{eq})-} 分别为平衡态分布 fi(eq)f_i^{(\mathrm{eq})} 的对称部分和非对称部分，即

fieq=wiρ[1+ci⋅ucs2+(ci⋅u)22cs4−u22cs2],fi(eq)+=12(fi(eq)(x,t)+f−i(eq)(x,t))=wiρ[1+(ci⋅u)22cs4−u22cs2],fi(eq)−=12(fi(eq)(x,t)−f−i(eq)(x,t))=wiρci⋅ucs2\begin{aligned} f^{\mathrm{eq}}_i &= w_i \rho \left[ 1 + \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u})^2}{2 c_s^4} - \frac{u^2}{2 c_s^2} \right] ,\\ f_i^{(\mathrm{eq})+} &= \frac{1}{2}\left( f_i^{(\mathrm{eq})}(\boldsymbol{x}, t) + f_{-i}^{(\mathrm{eq})}(\boldsymbol{x}, t) \right) =w_i \rho \left[ 1 + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u})^2}{2 c_s^4} - \frac{u^2}{2 c_s^2} \right] ,\\ f_i^{(\mathrm{eq})-} &= \frac{1}{2}\left( f_i^{(\mathrm{eq})}(\boldsymbol{x}, t) - f_{-i}^{(\mathrm{eq})}(\boldsymbol{x}, t) \right) = w_i \rho \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} \end{aligned}

  • wiw_i 为 ci\boldsymbol{c}_{i} 方向的权重， ρ\rho 和 u\boldsymbol{u} 为宏观密度和速度， csc_s 为DnQb模型的格子声速。

  • Fi\mathbb{F}_i 为源项 F\boldsymbol{F} 在 ii 方向上的函数。

TRT 的外力项

下面记 ti=wi/cs2t_i = {w_i} / {c_s^2}，以及 Fi+=12(Fi+F−i)F_{i}^{+} = \frac{1}{2} (\mathbb{F}_i + \mathbb{F}_{-i}) 和 Fi−=12(Fi−F−i)F_{i}^{-} = \frac{1}{2} (\mathbb{F}_i - \mathbb{F}_{-i}) 。

Postma 等[2]的研究对几种主流的源项模型进行了研究。不过他原文里列出的源项方程我还有些看不懂，我就照着Bawazeer等的文章抄了两个简单的源项：

（1）Buick and Greated [3]:

Fi=(1−12τ)tqci⋅F\mathbb{F}_i = \left(1 - \frac{1}{2 \tau} \right) t_q \boldsymbol{c}_{i} \cdot \boldsymbol{F}

转写为正负两项后表示为：

Fi+=0,Fi−=(1−12τ−)tici⋅FF_{i}^{+} = 0, \quad F_{i}^{-} = \left(1 - \frac{1}{2 \tau^{-}} \right) t_i \boldsymbol{c}_{i} \cdot \boldsymbol{F}

（2）Guo et al. [4]:

Fi=(1−12τ)wi⋅[ci−ucs2+(ci⋅u)cics4]⋅F\mathbb{F}_{i} = \left(1 - \frac{1}{2 \tau} \right) w_i \cdot \left[\frac{\boldsymbol{c}_i - \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u}) \boldsymbol{c}_i}{c_s^4}\right] \cdot \boldsymbol{F}

转写为正负两项后表示为：

Fi+=(1−12τ+)ti⋅[−u+(ci⋅u)cics2]⋅FFi−=(1−12τ−)tici⋅F\begin{aligned} F_i^{+} =& \left(1 - \frac{1}{2 \tau^{+}} \right) t_i \cdot \left[-\boldsymbol{u} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u}) \boldsymbol{c}_i}{c_s^2}\right] \cdot \boldsymbol{F} \\ F_i^{-} =& \left(1 - \frac{1}{2 \tau^{-}} \right) t_i \boldsymbol{c}_i \cdot \boldsymbol{F} \\ \end{aligned}

宏观密度 ρ\rho 和动量 j\boldsymbol{j} 的计算均表示为：

ρ=∑i=0b−1fi,j=ρu=∑i=0b−1fici+F2\rho = \sum\limits_{i=0}^{b-1} f_i ,\quad \boldsymbol{j} = \rho \boldsymbol{u} = \sum\limits_{i=0}^{b-1} f_i \boldsymbol{c}_{i} + \frac{\boldsymbol{F}}{2}

单步TRT-LBM的计算流程

(1) 计算 fi+f_i^{+}、 fi−f_i^{-} 以及 fi(eq)+f_i^{(\mathrm{eq})+}、 fi(eq)−f_i^{(\mathrm{eq})-}。

(2) 计算碰撞步

fi(pc)(x,t)=fi(x,t)−1τ+[fi+−fi(eq)+]−1τ−[fi−−fi(eq)−]+FiΔt(2)f_i^{(\mathrm{pc})}(\boldsymbol{x}, t) = f_i(\boldsymbol{x}, t) -\frac{1}{\tau^{+}} \left[ f_i^{+} - f_i^{(\mathrm{eq})+} \right] -\frac{1}{\tau^{-}} \left[ f_i^{-} - f_i^{(\mathrm{eq})-} \right] + \mathbb{F}_i \Delta t \tag{2}

(3) 计算流动步（与常规LBM一致）

(4) 更新宏观量。

TRT-LBM 与其他 LBM 的联系

TRT 与单松弛 LBGK

将式(1)整理，可得式(3)：

fi(x+ciΔt,t+Δt)−fi(x,t)=−(12τ++12τ−)[fi(x,t)−fi(eq)]−(12τ+−12τ−)[f−i(x,t)−f−i(eq)]+FiΔt(3)\begin{aligned} f_i(\boldsymbol{x}+\boldsymbol{c}_i \Delta t, t + \Delta t) - f_i(\boldsymbol{x}, t) &= - \left(\frac{1}{2 \tau^{+}} + \frac{1}{2 \tau^{-}}\right) [f_i(\boldsymbol{x}, t) - f^{(\mathrm{eq})}_i] \\ &\quad -\left(\frac{1}{2 \tau^{+}} - \frac{1}{2 \tau^{-}}\right) [f_{-i}(\boldsymbol{x}, t) - f^{(\mathrm{eq})}_{-i}] \\ &\quad + \mathbb{F}_i \Delta t \end{aligned} \tag{3}

所以，当 τ+=τ−\tau^{+} = \tau^{-} 时，TRT 方程可以退化为单松弛 LBGK 方程。此时： Λ=ν2cs4\Lambda = \frac{\nu^2}{c_s^4}，τ+\tau^{+} 也回归到单松弛 LBGK 方程的定义。

TRT 与 MRT

下文为方便起见，记 w1=(1τ++1τ−)/2\mathtt{w}_1 = (\frac{1}{\tau^{+}} + \frac{1}{\tau^{-}})/2 和 w2=(1τ+−1τ−)/2\mathtt{w}_2 = (\frac{1}{\tau^{+}} - \frac{1}{\tau^{-}})/2 ，并以D2Q9模型为例。

先定义：

  • f=[f0(x,t),...,f8(x,t)]T\boldsymbol{f} = [f_0(\boldsymbol{x}, t), ..., f_8(\boldsymbol{x}, t)]^\mathrm{T} 分布函数向量；
  • f(pc)=[f0(pc)(x,t),...,f8(pc)(x,t)]T\boldsymbol{f}^{(\mathrm{pc})} = [f_0^{(\mathrm{pc})}(\boldsymbol{x}, t), ..., f_8^{(\mathrm{pc})}(\boldsymbol{x}, t)]^\mathrm{T} 碰撞后的分布函数向量；
  • f(eq)=[f0(eq)(x,t),...,f8(eq)(x,t)]T\boldsymbol{f}^{(\mathrm{eq})} = [f_0^{(\mathrm{eq})}(\boldsymbol{x}, t), ..., f_8^{(\mathrm{eq})}(\boldsymbol{x}, t)]^\mathrm{T} 平衡态分布函数向量。

则 TRT 方程的碰撞步在D2Q9模型下的矩阵形式为：

f(pc)=[H](f−f(eq))(4)\boldsymbol{f}^{(\mathrm{pc})} = [\bold{H}] (\boldsymbol{f} - \boldsymbol{f}^{(\mathrm{eq})}) \tag{4}

其中矩阵 [H][\bold{H}] 为：

[H]=[w1+w2w1w2w1w2w2w1w2w1w1w2w1w2w2w1w2w1][\bold{H}] = \begin{bmatrix} \mathtt{w}_1+\mathtt{w}_2 &&&&&&&\\ & \mathtt{w}_1 & & \mathtt{w}_2 &&&&&\\ && \mathtt{w}_1 & & \mathtt{w}_2 &&&&\\ & \mathtt{w}_2 & & \mathtt{w}_1 &&&&&\\ && \mathtt{w}_2 & & \mathtt{w}_1 &&&&\\ &&&&& \mathtt{w}_1 & & \mathtt{w}_2 &\\ &&&&&& \mathtt{w}_1 & & \mathtt{w}_2 \\ &&&&& \mathtt{w}_2 & & \mathtt{w}_1 &\\ &&&&&& \mathtt{w}_2 & & \mathtt{w}_1 \end{bmatrix}

对比 MRT 方程的碰撞步在D2Q9模型下的矩阵形式：

f(pc)=[M−1][S][M](f−f(eq))(5)\boldsymbol{f}^{(\mathrm{pc})} = [\bold{M}^{-1}][\bold{S}][\bold{M}] (\boldsymbol{f} - \boldsymbol{f}^{(\mathrm{eq})}) \tag{5}

其中 [S][\bold{S}] 为松弛系数的对角矩阵。变换矩阵 [M][\bold{M}] 为：

[M]=[111111111−4−1−1−1−122224−2−2−2−21111010−101−1−110−20201−1−110010−111−1−100−20211−1−101−11−10000000001−11−1][\bold{M}] = \begin{bmatrix} 1&1&1&1&1&1&1&1&1\\ -4&-1&-1&-1&-1&2&2&2&2\\ 4&-2&-2&-2&-2&1&1&1&1\\ 0&1&0&-1&0&1&-1&-1&1\\ 0&-2&0&2&0&1&-1&-1&1\\ 0&0&1&0&-1&1&1&-1&-1\\ 0&0&-2&0&2&1&1&-1&-1\\ 0&1&-1&1&-1&0&0&0&0\\ 0&0&0&0&0&1&-1&1&-1 \end{bmatrix}

求解方程 [H]=[M−1][S][M][\bold{H}] = [\bold{M}^{-1}][\bold{S}][\bold{M}]，可得：

[S]=[M][H][M−1]=diag(w1+w2,w1+w2,w1+w2,w1−w2,w1−w2,w1−w2,w1−w2,w1+w2,w1+w2)\begin{aligned} [\bold{S}] =& [\bold{M}][\bold{H}][\bold{M}^{-1}] \\ =&\mathrm{diag}( \mathtt{w}_1+\mathtt{w}_2, \mathtt{w}_1+\mathtt{w}_2, \mathtt{w}_1+\mathtt{w}_2,\\&\quad \mathtt{w}_1-\mathtt{w}_2, \mathtt{w}_1-\mathtt{w}_2, \mathtt{w}_1-\mathtt{w}_2,\\&\quad \mathtt{w}_1-\mathtt{w}_2, \mathtt{w}_1+\mathtt{w}_2, \mathtt{w}_1+\mathtt{w}_2 ) \end{aligned}

可见 TRT 模型就是 MRT 模型的特例，将 MRT 中所使用的诸多松弛系数减少至到 2 个（1τ+\frac{1}{\tau^{+}} 和 1τ−\frac{1}{\tau^{-}}）。

  1. d’Humières, D., & Ginzburg, I. (2009). Viscosity independent numerical errors for Lattice Boltzmann models: From recurrence equations to “magic” collision numbers. Computers & Mathematics with Applications, 58(5), 823–840. https://doi.org/10.1016/j.camwa.2009.02.008 ↩︎

  2. Postma B, Silva G. Force methods for the two-relaxation-times lattice Boltzmann[J]. Physical Review E, 2020, 102(6): 063307. ↩︎

  3. Buick, J. M., Greated, C. A. Gravity in a lattice Boltzmann model[J]. Physical Review E, 2000, 61(5): 5307–5320. ↩︎

  4. Guo, Z., Zheng, C., & Shi, B. (2002). Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys. Rev. E, 65, 046308. ↩︎

发表于 更新于 分类于 ，

在 “【笔记】格子Boltzmann方法的相场模型” 中，我把LBM里的相场模型简单记了一下。笔记里有提到，序参数 ϕ\phi 的梯度可以通过下面的差分进行计算。

{∇ϕ(x)=∑i≠0wiciϕ(x+ciδt)cs2δt,∇2ϕ(x)=∑i≠02wics2δt2[ϕ(x+ciδ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}

其中 ii 为 DnQb 离散速度模型中的第 ii 个离散速度 ci\boldsymbol{c}_i 的下标， wiw_i 是 ci\boldsymbol{c}_i 的权重系数。

序参数 Φ 的一阶导数

对于点 x\boldsymbol{x} 及其邻近格子 x+ciδt\boldsymbol{x}+\boldsymbol{c}_i \delta_t，考虑序参数 ϕ\phi 在点 x\boldsymbol{x} 的Taylor展开，得：

ϕ(x+ciδt)=ϕ(x)+(ciδt)⋅∇ϕ(x)+δt22(ci⋅∇)2ϕ(x)+O(δt3)(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}

式(1)两边同时乘上 wiciw_i \boldsymbol{c}_i，并对所有方向求和，得到：

∑iwiciϕ(x+ciδt)=∑iwiciϕ(x)+∑iwici[(ciδt)⋅∇ϕ(x)]+∑iwici[δt22(ci⋅∇)2ϕ(x)]+O(δt3)(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}

下面分析式(2)右侧的三个项。

（1）对第1项，ϕ(x)\phi(\boldsymbol{x}) 是定值，所以

∑iwiciϕ(x)=ϕ(x)⋅∑i(wici)\sum_i w_i \boldsymbol{c}_i \phi(\boldsymbol{x}) = \phi(\boldsymbol{x}) \cdot \sum_i (w_i \boldsymbol{c}_i )

由于 ∑iwici=0⃗\sum_i w_i \boldsymbol{c}_i = \vec{0}，所以第1项为 0 。

（2）对第2项，

∑iwici[(ciδt)⋅∇ϕ(x)]=δt∇ϕ(x)⋅∑i(wicici)=δtcs2∇ϕ(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(wicici)=cs2[I]\sum_i (w_i \boldsymbol{c}_i \boldsymbol{c}_i) = c_s^2 [\boldsymbol{I}] 。

这个关系可以由常规 DnQb 模型中的 ∑icicifieq=ρ(uu+cs2[I])\sum_i \boldsymbol{c}_i \boldsymbol{c}_i f_i^{eq} = \rho (\boldsymbol{u}\boldsymbol{u} + c_s^2 [\boldsymbol{I}]) 推出。（取 u=0⃗\boldsymbol{u}=\vec{0}）

（3）对第3项，为简化考虑，这里推导其中一个分量，其形式为：

δt22∑iwiciγ[ciαciβ∂2ϕ(x)∂xα∂xβ]=δt22∂2ϕ(x)∂xα∂xβ∑iwiciαciβciγ\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}

由于在常规 DnQb 模型中，

∑iciαciβciγfieq=cs2ρ[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}]

这里的 δ\delta 为克罗内克函数。同样取 u=0⃗\boldsymbol{u}=\vec{0} ，得：

∑iwiciαciβciγ=cs2[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

所以式(2)被化简为

∑iwiciϕ(x+ciδt)=δtcs2∇ϕ(x)+O(δt3)\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)

即：

∇ϕ(x)=1δtcs2∑iwiciϕ(x+ciδt)+O(δt2)(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}

序参数 Φ 的拉普拉斯算子

将式(1)改写成

ϕ(x+ciδt)−ϕ(x)=(ciδt)⋅∇ϕ(x)+δt22(ci⋅∇)2ϕ(x)+O(δt3)(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}

式(4)两边同时乘上 wiw_i，并对所有方向求和，得到：

∑iwi[ϕ(x+ciδt)−ϕ(x)]=∑iwi[(ciδt)⋅∇ϕ(x)]+∑iwi[δt22(ci⋅∇)2ϕ(x)]+O(δt3)(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}

参照上一节的推导，易得式(5)等式右侧的两项分别满足

∑iwiciαδt∇αϕ(x)=δt∇αϕ(x)∑iwiciα=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

以及

∑iwiδt22ciαciβ⋅∇αβϕ(x)=δt22∇αβϕ(x)∑iwiciαciβ=δt2cs2δαβ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})

其中 ∇α=∂∂xα\nabla_{\alpha} = \frac{\partial}{\partial x_{\alpha}} 和 ∇αβ=∂2∂xα∂xβ\nabla_{\alpha\beta} = \frac{\partial^2}{\partial x_{\alpha}\partial x_{\beta}} 。而 δαβ\delta_{\alpha\beta} 同样是克罗内克函数。

所以式(5)表示为

∑iwi[ϕ(x+ciδt)−ϕ(x)]=δt2cs22∇2ϕ(x)+O(δt3)\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)

即：

∇2ϕ(x)=2δt2cs2∑iwi[ϕ(x+ciδt)−ϕ(x)]+O(δt2)(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}

关于中心差分

假设我们采用中心差分形式，需要对 x−ciδt\boldsymbol{x}-\boldsymbol{c}_i \delta_t 进行 Taylor 展开：

ϕ(x−ciδt)=ϕ(x)−(ciδt)⋅∇ϕ(x)+δt22(ci⋅∇)2ϕ(x)+O(δt3)(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}

（A）将式(1)与式(7)相减，偶数阶项（零阶项 ϕ(x)\phi(\boldsymbol{x}) 和二阶项）会被消去，得到：

ϕ(x+ciδt)−ϕ(x−ciδt)=2(ciδt)⋅∇ϕ(x)+O(δt3)\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)

两边同时乘上 wiciw_i \boldsymbol{c}_i 并对所有方向求和：

∑iwici[ϕ(x+ciδt)−ϕ(x−ciδt)]=2δt∇ϕ(x)⋅∑i(wicici)+O(δt3)\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)

利用前面提到的二阶矩性质 ∑iwicici=cs2I\sum_i w_i \boldsymbol{c}_i \boldsymbol{c}_i = c_s^2 \boldsymbol{I}，等式右侧化简为 2δtcs2∇ϕ(x)2\delta_t c_s^2 \nabla\phi(\boldsymbol{x})。因此：

∇ϕ(x)=12δtcs2∑iwici[ϕ(x+ciδt)−ϕ(x−ciδt)]+O(δt2)(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}

式(8)就是 Lee 等在文献[1]的式(64)的格式。

（B）将式(1)与式(7)相加，奇数阶项（一阶项和三阶项）会被消去，得到：

ϕ(x+ciδt)−2ϕ(x)+ϕ(x−ciδt)=δt2(ci⋅∇)2ϕ(x)+O(δt3)\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)

两边同时乘上 wiw_i 并对所有方向求和：

∑iwi[ϕ(x+ciδt)−2ϕ(x)+ϕ(x−ciδt)]=∑iδt2(ci⋅∇)2ϕ(x)+O(δt3)\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)

同理可得：

∇2ϕ(x)=1cs2δt2∑iwi[ϕ(x+ciδt)−2ϕ(x)+ϕ(x−ciδt)]+O(δt2)(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}

式(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.

0%