LBM-IBM-DEM相关论文分享2

CFD升级打怪的龙猫

已于 2024-06-26 21:01:12 修改

阅读量1k

点赞数 43

分类专栏： LBM-IBM-DEM初级进阶文章标签：算法学习方法

于 2024-06-26 20:55:02 首次发布

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

版权

LBM-IBM-DEM初级进阶专栏收录该内容

4 篇文章 1 订阅

订阅专栏

上一篇博文已经介绍了LBM-IBM耦合的其中一种方式，即在碰撞迁移公式中加入一个力项。并介绍了求该力项的方法。

在正式讲将DEM耦合进来之前，我想花一些篇幅讲颗粒在流体中可能受到的力。

本次分享的论文来源：Lu, Taijia, et al. "Pore‐scale simulation of the influence of temperature on particle motion based on LBM‐IBM‐DEM method." Hydrological Processes 37.9 (2023): e14981.

在多孔介质流体模拟中，悬浮颗粒会受到多种力的共同作用，这些力影响颗粒的运动状态和沉积行为。根据2.2节的内容，主要的力包括：

1. 流体力（Fluid Forces）：

阻力（Drag Force）：流体流动相对于颗粒表面产生的阻力，使颗粒减速。

升力（Lift Force）：流体流动在颗粒上下表面产生的压差，导致颗粒垂直于流动方向的运动。

由以下公式计算

流体对颗粒施加的力，影响颗粒的线速度和加速度

$F_f = \left( \frac{\Delta x}{\Delta t} \right)^2 \sum_{k=1}^{n} \left( B_{sk} \sum_{i=0}^{8} \Omega_i^s c_i \right)$

流体对颗粒施加的力矩，影响颗粒的角速度和角加速度。

$T_f = \left( \frac{\Delta x}{\Delta t} \right)^2 \sum_{k=1}^{n} \left( \left( x_k - x_p \right) \times B_{s_k} \sum_{i=0}^{8} \Omega_i^s c_i \right)$

$T_f$ ：流体对颗粒施加的力矩。
$\Delta x$ ：格子间距。
$\Delta t$ ：时间步长。
$n$ ：与颗粒重叠的LBM节点总数。
$x_k$ ：第k个节点的位置。
$x_p$ ：颗粒的质心位置，本论文中是球形颗粒，也为几何中心。
$B_{s_k}$ ：额外碰撞项的权重函数，用于调整颗粒节点的碰撞项。
$\Omega_i^s$ ：额外的碰撞项，表示粒子在流体中的影响。
$c_i$ ：离散速度方向的速度。
( $x_k$ − $x_p$ )：从颗粒质心到流体格子节点的位置向量。

再解释一下 $x_k$ 和 $x_p$ 的关系：

如上图所示，在LBM-IBM-DEM模型中，通常固体颗粒是拉格朗日物体，而流体网格是欧拉网格， $x_p$ 是固体颗粒的质心位置。拉格朗日物体映射到欧拉网格中会覆盖掉一些流体网格点，称这些被覆盖的流体网格点的位置为 $x_k$ 。在代码中，我们会将lagrangian points的属性映射到Euler grid上，这一点我将来会分享。

我们继续看 $B_{s_k}$ 和 $\Omega_i^s$ ：

$B_s (\epsilon_s, \tau_{BGK}) = \frac{\epsilon_s \left( \frac{\tau_{BGK}}{\Delta t} - \frac{1}{2} \right)}{(1 - \epsilon_s) + \epsilon_s \left( \frac{\tau_{BGK}}{\Delta t} - \frac{1}{2} \right)}$

式中， $B_s$ ：附加碰撞项的权重函数； $\epsilon_s$ ：表示固体体积分数的参数； $\tau_{BGK}$ : BGK模型中的弛豫时间; $\Delta t$ : 时间步长。

$\Omega_i^s = \left[ f_{-i}(x,t) - f_i^{eq}(\rho, u) \right] - \left[ f_i(x,t) - f_i^{eq}(\rho, v_p) \right]$

式中， $\Omega_i^s$ ：附加的碰撞项，用于表示固体对流体的影响； $f_{-i}(x,t)$ ：在位置 x 和时间 t 上沿着-i方向的分布函数； $f_{i}(x,t)$ ：在位置 x 和时间 t 上沿着方向 -i 的分布函数； $f_i^{eq}(\rho, u)$ ：基于流体密度 $\rho$ 和流体速度 $u$ 的平衡分布函数; $f_i^{eq}(\rho, v_p)$ ：基于流体密度 $\rho$ 和颗粒速度 $v_p$ 的平衡分布函数。

2. 重力（Gravity Force）：

作用：由于重力作用，颗粒受到向下的力，影响颗粒的沉降速度和分布，

$F_g = \left( 1 - \frac{\rho}{\rho_p} \right) g$

其中 $\rho_p$ 是颗粒密度， $g$ 是重力加速度。

该公式表示颗粒在流体中的重力，考虑了浮力的影响。公式中的 $( 1 - \frac{\rho}{\rho_p})$ 部分用来计算颗粒在流体中的有效重力。

重力的实际作用力是颗粒自身的重力减去浮力，浮力等于排开流体的重量。

Fg=(颗粒的重量)−(浮力)

将颗粒的重量和浮力代入公式，我们得到：

$F_g = \rho_p V g - \rho V g$

$F_g = (\rho_p - \rho) V g$

$F_g = m \left( 1 - \frac{\rho}{\rho_p} \right) g$

这里，我们假设 m=1m = 1m=1 ，则得到公式：

$F_g = \left( 1 - \frac{\rho}{\rho_p} \right) g$

3. 布朗力（Brownian Force）：

作用：由于流体分子的随机热运动，特别是对微小颗粒，布朗力引起颗粒的随机运动。

$F_B = \xi \sqrt{\frac{12 \pi k_B T \mu}{\Delta t}}$

$\xi$ ：一个随机向量，具有高斯随机数的特性，代表布朗运动的随机性。

$k_B$ ：玻尔兹曼常数（Boltzmann constant），其值为 1.3807×10^−23 J/K，用于将温度转换为能量。

T：流体温度，以开尔文（Kelvin）为单位。温度越高，布朗力越大。

μ：流体的动力粘度（dynamic viscosity），表征流体的内部摩擦力。

Δt：时间步长，表示模拟中每一步的时间间隔。时间步长越小，布朗力的变化越频繁。

对于微小颗粒（通常在1微米及以下），布朗力是显著的。这是因为这些颗粒的质量很小，流体分子的随机热运动对它们的影响较大，导致颗粒表现出显著的随机运动。对于孔隙尺度的研究，通常固体颗粒都在微米以内，故而应考虑布朗力。

对于较大颗粒，布朗力的影响可以忽略不计。

4. 碰撞力（Collision Force）：

作用：颗粒与颗粒、颗粒与多孔介质骨架之间的碰撞力，包含法向碰撞力和切向碰撞力。

$F_c = F_c^n + F_c^t$

（1）法向碰撞力 $F_c^n$ 的计算：

$F_c^n = \begin{cases} -k_n \delta_{ij} n_{ij} - \zeta_n u_{r,ij}^n, & \delta_{ij} > 0 \\ 0, & \delta_{ij} \leq 0 \end{cases}$

法向碰撞力是根据重叠量 $\delta_{ij}$ 和法向相对速度 $u_{r,ij}^n$ 计算得出的。

当 $\delta_{ij}$ 大于零（即颗粒接触或重叠时），法向碰撞力由弹簧力（ $-k_n \delta_{ij} n_{ij}$ ）和阻尼力（ $- \zeta_n u_{r,ij}^n$ ）组成。

当 $\delta_{ij} \leq 0$ （即颗粒没有接触或重叠时），法向碰撞力为零。

（2）切向碰撞力 $F_c^t$ 的计算：

$F_c^t = \begin{cases} - \min \left( \mu_c \left| F_c^n \right|, \zeta_t \left| u_{r,ij}^t \right| \right) t_{ij}, & \left| u_{r,ij}^t \right| > 0 \\ 0, & \left| u_{r,ij}^t \right| \leq 0 \end{cases}$

切向碰撞力是根据法向碰撞力的大小、切向相对速度和摩擦系数计算得出的。

当切向相对速度 $\left| u_{r,ij}^t \right| > 0$ 大于零时，切向碰撞力由摩擦力决定。当切向相对速度 $\left| u_{r,ij}^t \right| \leq 0$ 时，切向碰撞力为零。

（3）碰撞力矩 $T_c$ 的计算：

$T_c = \sum_{j=1}^{N} \left( \frac{d_p}{2} F_{c,ij}^t \times n_{ij} \right)$

碰撞力矩是由所有切向碰撞力对颗粒质心的力矩之和。计算时需要考虑每个切向碰撞力 $F_{c,ij}^t$ 和接触点到质心的向量 $n_{ij}$ 。

5. DLVO力（DLVO Forces）：

范德华力（Van der Waals Force）：是短程吸引力。

（1）当两个颗粒之间的距离较近时起主要作用。通过Hamaker常数和颗粒的几何特性来计算。

$F_{ad} = - \frac{A_H (4r_i r_j)^3}{6} \left( \frac{S_{ij}}{(S_{ij}^2 - (r_i + r_j)^2)^2\cdot\ (S_{ij}^2 - (r_i + r_j)^2 + 4r_i r_j)^2} \right)$

其中 $A_H$ 是Hamaker常数， $S_{ij}$ 是两个颗粒的中心距离， $r_i$ 和 $r_j$ 是两个颗粒的半径。

（2）当两个颗粒接触时，使用Johnson-Kendall-Roberts (JKR)理论来计算吸附力：

$F_{ad} = \begin{cases} \frac{3}{2} \pi r w_{ij}, & \text{particle-skeleton} \\ \frac{3}{4} \pi r w_{ij}, & \text{particle-particle} \end{cases}$

其中， $w_{ij}$ 是每单位面积的有效附着功（它的具体来源和计算取决于颗粒表面材料的性质以及颗粒之间的相互作用类型。）， $l$ 表示液相。请注意：骨架（多孔介质的固相）被视为一个平面或者具有较大半径的固体结构，相对于颗粒半径 $r$ 来说，可以简化为使用颗粒半径 rrr 来计算。

静电排斥力（Electrostatic Repulsive Force）：需要注意得是颗粒之静电排斥力是由于颗粒表面的电荷相互排斥产生的力。在带电颗粒悬浮液中，这种力在防止颗粒聚集和稳定悬浮液中起关键作用。它是长程力，随着颗粒之间的距离增加而显著减弱。间的电荷相互作用产生的排斥力。

$F_{Re} = 64 \pi \epsilon_r \epsilon_0 \kappa \left( \frac{RT}{zF} \right)^2 \cdot \tanh \left( \frac{zF \psi_i}{4RT} \right) \cdot \tanh \left( \frac{zF \psi_j}{4RT} \right) \cdot \frac{r_i r_j}{r_i + r_j} \cdot e^{-\kappa h_{ij}}$

$\kappa = \sqrt{\frac{2 e_l^2 N_A I}{\epsilon_r \epsilon_0 k_B T}}$

其中 $\epsilon_r$ 是流体的相对介电常数， $\epsilon_0$ 是真空的介电常数， $R$ 是气体常数， $z$ 是电解质的价态， $F$ 是法拉第常数， $\psi$ 是表面Zeta电位， $h_{ij}$ 是颗粒的表面距离。 $\kappa$ 是Debye-Hückel常数，表示静电力的作用范围。 $e_l$ 是电子的基本电荷， $N_A$ 是阿伏加德罗常数， $I$ 是溶液的离子浓度。