LBM-IBM-DEM进阶之路4----Two-relaxation time (TRT) LBM 双松弛时间LBM（之1）

CFD升级打怪的龙猫

已于 2024-07-03 21:53:11 修改

阅读量607

点赞数 17

文章标签：算法学习

于 2024-07-03 21:42:13 首次发布

本文链接：https://blog.csdn.net/2401_82853022/article/details/140152653

版权

本节我们来看双松弛时间LBM，TRT LBM的执行方法。

双松弛时间（TRT）模型是对传统的单松弛时间（SRT）和多松弛时间（MRT）模型的一种改进。TRT的核心思想在于：

松弛时间的分离：TRT模型将粒子分布函数分解为对称和反对称分量，并分别使用不同的松弛时间（ $\tau^+$ 和 $\tau^-$ ）进行松弛。这种方法允许更精细地控制扩散和对流过程的数值稳定性和准确性。
更高的数值稳定性：通过分别调节对称和反对称分量的松弛时间，TRT模型可以在保持高数值精度的同时，提高系统的数值稳定性。

1. 背景及初始化

TRT LBM中溶质的热力学状态由一个Q维粒子分布函数 $f_q(\mathbf{r}, t)$ 定义，其中q=0,…,Q−1。该函数在每个格点( $\mathbf{r}$ )和每个离散时间( $t$ )上定义。格点通过一组与格子轴和对角线对齐的离散速度 $\mathbf{c}_q$ 相互连接。对所有离散速度求和粒子分布得到溶质浓度 $C$ ：

$C = \sum_{q=0}^{Q-1} f_q$

在每个节点上，粒子分布函数可以分解为对称和反对称分量，

$f_q = f_q^+ + f_q^-$

其中， $f_q^+ = \frac{f_q + f_{\bar{q}}}{2}$ ， $f_q^- = \frac{f_q - f_{\bar{q}}}{2}$

TRT LBM的物理解释:

对称分量反映扩散：对称分量 $f_q^+$ 反映了粒子分布的平均行为，可以与经典的扩散现象相对应。
反对称分量反映对流：反对称分量 $f_q^-$ 反映了粒子分布的偏移行为，可以与经典的对流现象相对应。

整数 $\bar{q}$ 是指向 $\mathbf{c}_q$ 反方向的速度索引 $(\mathbf{c}_{\bar{q}} = -\mathbf{c}_q)$ 。

对于零速度的其余粒子 $(\mathbf{c}_0 = 0)$ ， $f_0^+ = f_0 = C - \sum_{q=1}^{Q-1} f_q$ ， $f_0^- = 0$

为了简化对所有速度 $\mathbf{c}_q$ 的求和，假设速度按以下顺序排列：

$\bar{q} = q + \frac{(Q-1)}{2}$ ， $1 \leq q \leq \frac{(Q-1)}{2}$ ，这里举例：对于D2Q5, 3 = 1 + (5-1)/2; 4 = 1 + (5-1)/2;

$f_0^+ = C - 2 \sum_{q=1}^{(Q-1)/2} f_q^+$

2. 碰撞迁移方程

一个节点上的粒子通过非零离散速度 $\mathbf{c}_q$ $(\mathbf{c}_q \neq 0 \ldots, Q-1)$ 移动到邻近节点，或者对于 $\mathbf{c}_0$ 保持在节点上。一旦粒子到达邻近节点，结果分布通过TRT碰撞算子放松到平衡状态：

$f_q(\mathbf{r} + \mathbf{c}_q, t + 1) = f_q(\mathbf{r}, t) - \frac{1}{\tau^+}(f_q^+ - e_q^+) - \frac{1}{\tau^-}(f_q^- - e_q^-)$

其中 $e_q^+$ 和 $e_q^-$ 是平衡粒子分布的对称和反对称分量 $(e_q = e_q^+ + e_q^-)$ ； $\tau^+$ 和 $\tau^-$ 分别是对称和反对称松弛参数； $f_q(\mathbf{r}, t)$ 表示在时间 $t$ 时，位置 $\mathbf{r}$ 处的粒子分布函数； $\frac{1}{\tau^+}(f_q^+ - e_q^+)$ 是对称分量的碰撞项， $\frac{1}{\tau^-}(f_q^- - e_q^-)$ 是反对称分量的碰撞项。 $\frac{1}{\tau^+}(f_q^+ - e_q^+) + \frac{1}{\tau^-}(f_q^- - e_q^-)$ 为碰撞项，是粒子在碰撞过程中趋向于平衡态的变化量。将当前时刻的粒子分布函数减去对称和反对称分量的碰撞项，其目的是为了模拟流体中粒子在碰撞过程中的松弛效应，使其趋向于平衡态。

记上述迁移方程的整个右侧称为碰撞后粒子分布 $\tilde{f}_q(\mathbf{r}, t)$ 。

3. 平衡态的粒子分布函数

非零-离散速度的平衡粒子分布函数可表示为(Ginzburg, 2013)：

平衡分布函数 $e_q$ 分为对称分量 $e_q^+$ 和反对称分量 $e_q^-$ ，用来分别处理扩散和对流现象。

$e_q^+ = C \left( t_q^{(m)} c_s^2 + t_q^{(u)} \mathbf{V}^2 + w_q^{(u)} \|\mathbf{c}_q\|^2 \sum_{\alpha=1}^{d} (V_\alpha^2 - \bar{V}^2)c_{q\alpha}^2 \right) + \sum_{\beta \neq \alpha} \frac{V_\alpha V_\beta c_{q\alpha} c_{q\beta}}{2 \sum_{j=1}^{(Q-1)/2} c_{j\alpha}^2 c_{j\beta}^2}$

式中，

$e_q^+$ ：粒子分布函数的对称平衡分量。对称平衡分量主要反映了系统的扩散特性。
$C$ ：溶质浓度。
$t_q^{(m)}$ ：与质量相关的系数。
$c_s^2$ ：声速的平方。
$t_q^{(u)}$ ：与速度相关的系数。
$\mathbf{V}^2$ ：对流速度的平方。
$w_q^{(u)}$ ：与速度相关的权重系数。
$\|\mathbf{c}_q\|^2$ ：离散速度的平方。
$V_\alpha$ ：对流速度的第α\alphaα个分量。
$\bar{V}^2$ ：对流速度的均方值。
$c_{q\alpha}$ ：离散速度的第 $\alpha$ 个分量。
$d$ ：维度。
$\sum_{\beta \neq \alpha}$ ：表示对所有不同于 $\alpha$ 的分量进行求和。
$\sum_{j=1}^{(Q-1)/2}$ ：表示对前半部分速度的求和。

上式展示了如何计算对称平衡分量 $e_q^+$ ，它综合了质量、速度和离散速度的各种影响因素。

$e_q^- = t_q^{(a)} C \sum_{\alpha=1}^{d} V_\alpha c_{q\alpha}$

式中，

$e_q^-$ ：这是粒子分布函数的反对称平衡分量。反对称平衡分量主要反映了系统的对流特性。
$t_q^{(a)}$ ：这是一个与对称性相关的系数。
$C$ ：这是溶质浓度。
$\sum_{\alpha=1}^{d}$ ：这表示对所有维度进行求和。
$V_\alpha$ ：这是对流速度的第 $\alpha$ 个分量。
$c_{q\alpha}$ ：这是离散速度的第 $\alpha$ 个分量。

上式展示了如何计算反对称平衡分量 $e_q^-$ ，它主要依赖于对流速度和离散速度的分量。

上面两个式子中， $t_q^{(m)}$ , $t_q^{(u)}$ , $t_q^{(a)}$ 分别是与质量、速度和对称性相关的系数， $w_q^{(u)}$ 是与速度相关的权重系数，这些权重都是非负的并且都是各项同性的，满足各项同性的条件：

$\sum_q w_q^{(\cdot)} c_{q\alpha} c_{q\beta} = \delta_{\alpha \beta}$ ， $\sum_q t_q^{(\cdot)} c_{q\alpha} c_{q\beta} = \delta_{\alpha \beta}$

其中 $\alpha, \beta = 1, \ldots, d$ 。不同的权重选择会导致不同的数值稳定性（Ginzburg 等，2010）。

$\delta_{\alpha \beta}$ 是克罗内克δ符号，它在 α = β 时为1，其他情况下为0。

另一半的平衡粒子分布（对于 $q = (Q-1)/2 + 1, \ldots, Q-1$ ）可以通过对称和反对称关系计算：

粒子分布函数反方向 $\bar{q}$ 的对称平衡分量等于粒子分布函数的对称平衡分量： $e_{\bar{q}}^+ = e_q^+$ ,

这是因为在TRT模型中，对称分量表示的是粒子分布的均匀部分，它与方向无关。因此，对于反方向的速度 $\bar{q}$ ，对称分量 $e_{\bar{q}}^+$ 与原速度 $q$ 的对称分量 $e_q^+$ 是相同的。

粒子分布函数反方向 $\bar{q}$ 的反对称平衡分量等于粒子分布函数的反对称平衡分量的负数: $e_{\bar{q}}^- = -e_q^-$

这是因为反对称分量描述的是方向性特征，当速度方向取反时，这些特征也会反转。因此，对于反方向的速度 $\bar{q}$ ，反对称分量 $e_{\bar{q}}^-$ 与原速度 $q$ 的反对称分量 $e_q^-$ 的符号相反。

对于0方向的粒子： $e_0^+ = C - 2 \sum_{q=1}^{(Q-1)/2} e_q^+$ , $e_0^- = 0$

将速度索引分为两个部分后，在计算某些量时可以减少求和的范围。例如，在计算平衡分量时，可以只对一半的速度索引进行求和，然后根据对称性得到另一半的结果。

4. TRT-LBM碰撞迁移总结

（1）初始化分布函数

在格子玻尔兹曼方法中，初始状态下的粒子分布函数 $f_q(\mathbf{r}, 0)$ 通常根据宏观物理量（如密度 $\rho$ 和速度 $\mathbf{u}$ ）以及对称和反对称分量的平衡态分布函数 $e_q^+$ 和 $e_q^-$ 来初始化。

非零方向：

$e_q^- = t_q^{(a)} C \sum_{\alpha=1}^{d} V_\alpha c_{q\alpha}$

零方向： $e_0^+ = C - 2 \sum_{q=1}^{(Q-1)/2} e_q^+$ , $e_0^- = 0$

请注意上面的 $e_q^+$ 和 $e_q^-$ 的形式比较复杂，涉及到的参数比较多，文献中给定以下特例（对于D3Q15):

source: Yan et al., 2017

在模拟中使用了具有15个离散速度的3D格子结构（D3Q15）。D3Q15模型中的离散速度可以描述为：

其中， $\|\mathbf{c}_q\|^2 = 1$ 的 $\mathbf{c}_q$ 被分类为I型速度，而 $\|\mathbf{c}_q\|^2 = 3$ 的 $\mathbf{c}_q$ 被分类为II型速度。我们可以选择常用的“水动力学”各向同性权重（Ginzburg 等，2010）：

对于I型速度: $t_q^{(a)} = t_q^{(m)} = 1/3, \quad t_q^{(u)} = 0, \quad w_q^{(u)} = 1/2$ ；

对于II型速度: $t_q^{(a)} = t_q^{(m)} = 1/24, \quad t_q^{(u)} = 1/8, \quad w_q^{(u)} = 0$ .

对于I型和II型速度，选择了 $c_s^2 = 3/8$ 。因此，平衡分布变为：

$e_q = \frac{1}{8} C + \frac{1}{3} C \mathbf{c}_q \cdot \mathbf{V} + \frac{1}{2} C (\mathbf{c}_q \cdot \mathbf{V})^2 - \frac{1}{6} C \mathbf{V} \cdot \mathbf{V}, \quad \mathbf{c}_q \in I$

$e_q = \frac{1}{64} C + \frac{1}{24} C \mathbf{c}_q \cdot \mathbf{V} + \frac{1}{16} C (\mathbf{c}_q \cdot \mathbf{V})^2 - \frac{1}{48} C \mathbf{V} \cdot \mathbf{V}, \quad \mathbf{c}_q \in II$

$e_q = \frac{1}{8} C - \frac{1}{3} C \mathbf{c}_q \cdot \mathbf{V}$

$f_q^+ = \frac{e_q + e_{\bar{q}}}{2}$ (请注意，对于上面这个特例，要注意找好方向的 $e_q$ , 如 $e_3$ (0,0,1) 和 $e_{10}$ (0,0,-1)是相反方向，都使用I型速度公式来计算；如 $e_4$ (1,1,1) 和 $e_{11}$ (-1,-1,-1)是相反方向，都使用II型速度公式来计算；

$f_q^- = \frac{e_q - e_{\bar{q}}}{2}$

初始化： $f_q = f_q^+ + f_q^-$

（2）碰撞过程

在每个时间步长内，粒子分布函数 $f_q$ 首先经历碰撞过程。碰撞过程按照TRT模型的碰撞算子进行更新：

$\tilde{f}_q(\mathbf{r}, t) = f_q(\mathbf{r}, t) - \frac{1}{\tau^+}(f_q^+ - e_q^+) - \frac{1}{\tau^-}(f_q^- - e_q^-)$

（3）迁移过程

在碰撞过程之后，粒子分布函数经历迁移过程，即从当前位置 $\mathbf{r}$ 迁移到新的位置 $\mathbf{r} + \mathbf{c}_q$

$f_q(\mathbf{r} + \mathbf{c}_q, t + 1) = \tilde{f}_q(\mathbf{r}, t)$

5. 边界条件

为了使TRT方法尽可能简单，我们应用广泛使用的Standard bounce-back（SBB）边界条件：

$f_{\bar{q}}(\mathbf{r}, t + 1) = \tilde{f}_q(\mathbf{r}, t)$

检测固体边界节点：遍历所有节点，识别出固体边界节点。
反弹操作：
- 对于每个固体边界节点，检查哪些速度方向 $q$ 指向边界。
- 应用反弹回弹边界条件，将粒子反弹回边界节点。

如果存在质量源或汇（ $M$ ），碰撞后粒子分布修改为：

$\tilde{f}_q(\mathbf{r}, t) = \tilde{f}_q(\mathbf{r}, t) + t_q^{(m)} c_s^2 M, \quad q = 1, \ldots, Q-1$

$\tilde{f}_0(\mathbf{r}, t) = \tilde{f}_0(\mathbf{r}, t) + M \left( 1 - 2 c_s^2 \sum_{q=1}^{(Q-1)/2} t_q^{(m)} \right)$

其中，右侧的 $\tilde{f}_q$ 定义在碰撞迁移方程中。这一修改使得TRT LBM能够模拟反应性传输。

6. 宏观恢复

使用Chapman-Enskog展开，碰撞迁移方程的解接近以下无量纲反应对流扩散方程（ADE）的解（Ginzburg，2005a）：

$\frac{\partial C}{\partial T} + \nabla \cdot (\mathbf{V}C) = \hat{D} \nabla^2 C + M$

其中， $\hat{D}$ 是无量纲扩散系数， $\hat{D} = (\tau^- - 1/2)c_s^2$ 。这个无量纲反应ADE可以从相应的有量纲反应ADE中推导出来：

$\frac{\partial c}{\partial t} + \nabla \cdot (\mathbf{uc}) = D \nabla^2 c + I$

其中， $c$ 是溶质浓度（ $c = c_0$ ）， $t$ 是时间 $(t = T \Delta x / u_0)$ ， $u$ 是流动的对流速度 $(\mathbf{u} = u_0 \mathbf{V})$ ， $D$ 是溶质的扩散系数 $D = \Delta x u_0 \hat{D}$ ， $I$ 是反应项 $(I = c_0 u_0 M / \Delta x)$ 。 $c_0$ , $\Delta x$ 和 $u_0$ 是用于归一化有量纲ADE的特征值。

在TRT LBM应用于溶质的反应传输中，反对称松弛参数 $\tau^-$ 由溶质的扩散系数决定， $\tau^- = 1/2 + D / (c_s^2 \Delta x u_0)$ ，而对称松弛参数 $\tau^+$ 可以自由选择，但需满足 $\tau^+ > 0.5$ 。

请注意如果将 $\tau^+$ 赋值与 $\tau^-$ 相同的值，TRT方法将与SRT方法（Ginzburg, 2005a; Zou和He, 1997）一致。一个建议的取值方式是 $\tau^+ = 1/2 + 1/(4\tau^- - 2)$ ，这来自于所谓的“magic”乘积 $((\tau^+ - 1/2)(\tau^- - 1/2) = 1/4)$ ，即已知 $\tau^-$ 值可求 $\tau^+$ 同时能保持良好的稳定性（Ginzburg等，2010）。

7. 注意事项

请注意本次分享仅作为介绍TRT LBM方法的一个guide, 内容不够完整详细，其实这里面的水还是非常深的，大家看完应该感觉这里怎么这么多参数，这些参数如何来的？这些参数的具体值如何设置？

是的，LBM就是这样的，为了得到收敛性，我们不得不引入大量的参数来“制衡”和“修正”.

TRT LBM未完待续，将来博主会再次更新补充将这些疑惑，包括如何对于这些繁杂的公式进行编程。

参考文献：

Yan, Zhifeng, et al. "Two-relaxation-time lattice Boltzmann method and its application to advective-diffusive-reactive transport." Advances in water resources 109 (2017): 333-342.

Ginzburg, Irina. "Multiple anisotropic collisions for advection–diffusion Lattice Boltzmann schemes." Advances in water resources 51 (2013): 381-404.

Ginzburg, Irina, Dominique d’Humières, and Alexander Kuzmin. "Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium." Journal of Statistical Physics 139 (2010): 1090-1143.

Ginzburg, Irina. "Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation." Advances in water resources 28.11 (2005): 1171-1195.

Zou, Qisu, and Xiaoyi He. "On pressure and velocity boundary conditions for the lattice Boltzmann BGK model." Physics of fluids 9.6 (1997): 1591-1598.

如有错误，请评论区批评指正。

良心创作，转载请注明出处----CFD龙猫

CFD升级打怪的龙猫

关注

17
点赞
踩
9

收藏

觉得还不错? 一键收藏
0
评论
LBM-IBM-DEM进阶之路4----Two-relaxation time (TRT) LBM 双松弛时间LBM（之1）

双松弛时间（TRT）模型是对传统的单松弛时间（SRT）和多松弛时间（MRT）模型的一种改进。TRT的核心思想在于：松弛时间的分离：TRT模型将粒子分布函数分解为对称和反对称分量，并分别使用不同的松弛时间（和）进行松弛。这种方法允许更精细地控制扩散和对流过程的数值稳定性和准确性。更高的数值稳定性：通过分别调节对称和反对称分量的松弛时间，TRT模型可以在保持高数值精度的同时，提高系统的数值稳定性。
复制链接

扫一扫