Ives's Pages

发表于 分类于

最近在网上刷到了这篇 arxiv 预印本，

Zhou J G. Single MRT-featured Lattice Boltzmann Method (SmrtLBM)[PP/OL]. arXiv（2024-09-02）. http://arxiv.org/abs/2409.01076.

加上休息时间有空，就简单看了这篇文章并做了笔记。

基本原理

这篇文章以 D2Q9 离散速度模型为例，对其多松弛形式进行了简化。

多松弛模型 (MRT) 的简单介绍

定义：

  • 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} 平衡态分布函数向量。

这里的平衡态用的仍旧是 LBGK 模型的形式

fi(eq)=wiρ[1+ci⋅ucs2+(ci⋅u)22cs4−u22cs2]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]

其中 ci\boldsymbol{c}_i 为离散速度方向， wiw_i 为 ci\boldsymbol{c}_i 方向的权重系数， csc_s 为格子声速。 ρ\rho 和 u\boldsymbol{u} 分别为流场密度和速度。

在 D2Q9 模型中，cs=1/3c_s=1/\sqrt{3}，

ci={(0,0)i=0(±1,0),(0,±1)ci=1−4(±1,±1)ci=5−8\boldsymbol{c}_{i} = \begin{cases} (0,0) & i=0 \\ (\pm 1, 0), (0, \pm 1) c & i=1-4 \\ (\pm 1, \pm 1) c & i=5-8 \end{cases}

权重系数 w0=4/9w_0 = 4/9 , w1−4=1/9w_{1-4} = 1/9 , w5−8=1/36w_{5-8} = 1/36 .

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

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

松弛步则写作：

fi(x+ciδt,t+δt)=fi(pc)(x,t+δt)(2)f_i (\boldsymbol{x}+\boldsymbol{c}_i \delta_t, t + \delta_t) = f_i^{(\mathrm{pc})} (\boldsymbol{x}, t + \delta_t) \tag{2}

其中 δt\delta_t 为时间步长， [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}

流场密度 ρ\rho 和速度 u\boldsymbol{u} 分别为

ρ=∑ifi,ρu=∑icifi\rho = \sum_i f_i , \quad \rho\boldsymbol{u} = \sum_i \boldsymbol{c}_i f_i

郭照立等[1]记松弛时间的对角矩阵为：

[S]=diag(0,se,sε,0,sq,0,sq,sν,sν)[\bold{S}] = \text{diag}(0, s_e, s_{\varepsilon}, 0, s_q, 0, s_q, s_{\nu}, s_{\nu})

则运动粘度 ν=cs2(1sν−12)δt\nu = c_s^2 (\frac{1}{s_{\nu}} - \frac{1}{2}) \delta_t ，以及体黏性系数 ζ=cs2(1se−12)δt\zeta = c_s^2 (\frac{1}{s_e} - \frac{1}{2}) \delta_t 。

Zhou 的简化

Zhou[2]将矩阵 [S][\bold{S}] 改写成

[S]=[S0]+[St]=diag(1/τs,1/τs,1/τs,1/τs,1/τs,1/τs,1/τs,1/τ,1/τ),(3)\begin{aligned} [\bold{S}] &= [\bold{S}_0] + [\bold{S}_t] \\ &= \text{diag}( 1/\tau_s, 1/\tau_s, 1/\tau_s, 1/\tau_s, 1/\tau_s, 1/\tau_s, 1/\tau_s, 1/\tau, 1/\tau ), \end{aligned} \tag{3}

其中：

[S0]=(1τ−1τs)diag(0,0,0,0,0,0,0,1,1),[St]=1τs diag(1,1,1,1,1,1,1,1,1).\begin{aligned} [\bold{S}_0] &= \left( \frac{1}{\tau} - \frac{1}{\tau_s} \right) \text{diag}(0, 0, 0, 0, 0, 0, 0, 1, 1), \\ [\bold{S}_t] &= \frac{1}{\tau_s} \, \text{diag}(1, 1, 1, 1, 1, 1, 1, 1, 1). \end{aligned}

可见在这种构造形式下，τ=1/sν\tau= 1 / s_{\nu}，与LBGK模型一致。

为了实现化简，Zhou[2]令 τs=1\tau_s=1，并把式(3)代入式(1)中展开，整理得：

fi(pc)(x,t)=fi(eq)(x,t)+(−1)i(1−1τ)×∑k=18(−1)k4[1[1,2,3,4](i,k)+1[5,6,7,8](i,k)][fk(x,t)−fk(eq)(x,t)](4)\begin{aligned} f_i^{(\mathrm{pc})}(\boldsymbol{x},t) &= f_{i}^{(\mathrm{eq})}(\boldsymbol{x},t) + (-1)^i \left( 1 - \frac{1}{\tau} \right) \times \\ & \sum_{k=1}^{8} \frac{(-1)^k}{4} \left[ \boldsymbol{1}_{[1,2,3,4]}(i,k) + \boldsymbol{1}_{[5,6,7,8]}(i,k) \right] \left[ f_{k}(\boldsymbol{x},t) - f_{k}^{(\mathrm{eq})}(\boldsymbol{x},t) \right] \end{aligned} \tag{4}

式(4)中的函数 1[a,b,c,d](i,k)\boldsymbol{1}_{[a,b,c,d]}(i,k) 定义为

1[a,b,c,d](i,k)={1,if i in [a,b,c,d], and k in [a,b,c,d]0,otherwise\boldsymbol{1}_{[a,b,c,d]}(i,k) = \begin{cases} 1, & \text{if i in [a,b,c,d], and k in [a,b,c,d]}\\ 0, & \text{otherwise} \end{cases}

实例代码

这里顺手做了个双周期剪切流的简单算例（拿上次双松弛LBM的代码改的），设置如下。

  • 计算域： [0,L]×[0,L][0,L] \times [0,L]，其中 L=2πL=2\pi。
  • 雷诺数： Re=1/ν=1000Re = 1 / \nu = 1000。 其中 ν\nu 为流体运动粘度 (Physical)。
  • 初始速度场

{ux=U0tanh⁡(σ (14−∣yL−12∣))uy=U0δsin⁡(2π (xL+12))\begin{cases} u_x =& U_0 \tanh (\sigma\, (\frac{1}{4} - |\frac{y}{L} - \frac{1}{2}| )) \\ u_y =& U_0 \delta \sin (2 \pi \, (\frac{x}{L} + \frac{1}{2}) ) \end{cases}

  • 压强场 pp ： 用伪谱法求解压力泊松方程。
  • LBM松弛时间 τ\tau: 0.515
  • 网格大小： δx=L/N\delta_x = L/N (N=200N=200)。
  • 时间步 δt\delta_t ： 用 LBM 中 τ\tau 的关系反推。
  • 初始速度场参数： U0=1U_0 = 1, σ=30\sigma=30, δ=0.05\delta=0.05

代码我放在了Github Gist上，计算结果如下图所示

Reference

[1] 郭照立，郑楚光. 格子Boltzmann方法的原理及应用[M]. 科学出版社，2009.

[2] Zhou J G. Single MRT-featured Lattice Boltzmann Method (SmrtLBM)[PP/OL]. arXiv（2024-09-02）. http://arxiv.org/abs/2409.01076.

发表于 更新于 分类于

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

Concepts

基本方程

对于 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^{-}}）。

实例代码

为了练手，我写了一个双周期剪切流的代码。

  • 计算域： [0,L]×[0,L][0,L] \times [0,L]，其中 L=2πL=2\pi。
  • 网格大小： δx=L/N\delta_x = L/N (N=128N=128)。
  • 雷诺数： Re=1/ν=100Re = 1 / \nu = 100。 其中 ν\nu 为流体运动粘度 (Physical)。
  • 时间步 δt\delta_t ： 用 LBM 中 τ+\tau^{+} 的关系反推。
  • 速度场

{ux=U0tanh⁡(σ (14−∣yL−12∣))uy=U0δsin⁡(2π (xL+12))\begin{cases} u_x =& U_0 \tanh (\sigma\, (\frac{1}{4} - |\frac{y}{L} - \frac{1}{2}| )) \\ u_y =& U_0 \delta \sin (2 \pi \, (\frac{x}{L} + \frac{1}{2}) ) \end{cases}

  • 压强场 pp ： 求解压力泊松方程。由于这个场是周期的，我就用伪谱法算。

代码我放在了Github Gist上，欢迎直接查看

TRT-LBM 的碰撞算子会对全部对称非平衡分量进行松弛，隐式保留了非流体力学行为的伪模态（ghost modes）。这会在高雷诺数工况下降低计算稳定性，比如上面的代码在高 Re 下会崩溃。 Yu 等的文章中对此有相关说明，并给出了基于矩展开的修正 [5]。

  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. ↩︎

  5. Yu Y，Yu Y，Zhou Y，et al. Purified Two-Relaxation-Time Lattice Boltzmann Method: Removing Ghost Modes from TRT for Enhanced Stability[PP/OL]. arXiv（2026-02-06）. http://arxiv.org/abs/2602.06686. DOI:10.48550/arXiv.2602.06686. ↩︎

0%