LBM-IBM-DEM进阶之路3----DEM中的软球碰撞模型

上节我们了解到了流场中的固体颗粒可能会受到哪些力，比如流场作用力、重力、碰撞力、DLVO力等等。详情见：LBM-IBM-DEM相关论文分享2-CSDN博客

那么本节我们介绍颗粒的软球碰撞模型，让我们对LBM-IBM-DEM耦合方法逐渐加深理解。

软球模型是DEM中常用的用于分析两个离散固体颗粒之间或固体颗粒与墙壁之间的相互作用（接触力）的方法。如果两个离散固体颗粒相互碰撞，会有轻微的重叠，从而产生变形力。两个颗粒之间的相互作用被简化为减震器压缩反弹的过程，我们称之为弹簧阻尼模型，如下图所示。

source: Fo et al., 2022

当颗粒  和 碰撞时，会产生法向和切向接触力 和 。 

---

1. 法向力 根据线性弹簧模型计算：

其中， 是颗粒  和  接触的法向重叠（因采用的是软球模型，可以理解为两个小颗粒碰撞之后会有一部分形变。） 是依赖位移的弹性法向系数，主要取决于接触系统的几何和物理特性，如接触颗粒的半径、杨氏模量和泊松比。  是从颗粒 到颗粒  的单位矢量。法向阻尼系数 如下：

上式中系数的计算公式为：；

上式中系数的计算公式为：;

其中，​ 是通过实验确定的恢复系数（通常依据具体的case来定）。

以下解答一些疑惑：

（1）的来源

或许大家很疑惑的来源，在这里做一个简单介绍。

弹性法向系数由Hertz接触理论来定（可参考《颗粒物质力学导论》孙其诚&王光谦，p.33):

其中，和是颗粒材料的弹性模量和泊松比；是颗粒半径；下标分别代表发生接触的颗粒和。

如果颗粒和同质且粒径相等，则简化为：

在颗粒接触过程中，接触条件（如法向重叠、接触面积、材料属性等）可能会发生变化。接触刚度 ​ 的计算通常基于赫兹接触理论或其他复杂的接触模型，这些模型考虑了颗粒材料的弹性模量、泊松比等参数。由于这些参数在接触过程中可能会变化，因此严格上来说需要在每个模拟时步实计算。但这显然增加了计算量。因此在具体的case中，我们常常会普遍采用一个较小的弹性系数来获得较大的时间步长，进而加快计算过程。

（2）的来源

代表两个颗粒的有效质量 (reduced mass)，它是一个标准的物理量，常用于描述两个物体在碰撞或相互作用中的质量影响。这一概念在经典力学中被广泛应用，尤其是在描述双体问题以及粒子碰撞问题中。

（3）模拟中我们怎么确定两个小球碰撞所产生的法相重叠量

假设两个小球 和  的位置分别为 和，它们的半径分别为 和 。法向重叠量 ​ 可以通过以下几何关系确定：

其中：表示两个小球中心之间的距离； 和 是两个小球的半径。

首先计算两小球中心之间的距离：

其中 和 是小球 iii 和 jjj 的位置向量。

如果 大于 ，则说明两个小球没有接触，法向重叠量 为零。

以下是一个简单的C/C++实现：

#include <iostream>
#include <cmath>

// 定义小球结构
struct Sphere {
    double x, y, z;  // 小球中心的位置
    double radius;   // 小球的半径
};

// 计算两个小球之间的法向重叠量
double calculateOverlap(const Sphere& sphere1, const Sphere& sphere2) {
    // 计算两个小球中心之间的距离
    double dx = sphere1.x - sphere2.x;
    double dy = sphere1.y - sphere2.y;
    double dz = sphere1.z - sphere2.z;
    double distance = std::sqrt(dx * dx + dy * dy + dz * dz);
    
    // 计算法向重叠量
    double overlap = sphere1.radius + sphere2.radius - distance;
    
    // 如果没有重叠，则重叠量为零
    if (overlap < 0) {
        overlap = 0;
    }
    
    return overlap;
}

int main() {
    // 定义两个小球
    Sphere sphere1 = {1.0, 2.0, 3.0, 0.5};
    Sphere sphere2 = {1.5, 2.5, 3.5, 0.6};
    
    // 计算法向重叠量
    double overlap = calculateOverlap(sphere1, sphere2);
    
    std::cout << "The overlap between the two spheres is: " << overlap << std::endl;
    
    return 0;
}

请注意，上述代码仅是一个非常简单的case, 对于静态的两个小球。在实际的LBM-DEM模拟，小球是移动的，我们需要在每个时步中捕捉到小球的几何/质心坐标。另外一点需要提的是，如代码中的calculateOverlap()函数，通常我们在编写大规模C++并行程序时，我们会将这些函数在封装到相应的头文件中（.h)，然后在具体的case中(.cpp)来引用该头文件，并按需调用相应的函数。

---

2. 切向接触力 的计算公式如下：

其中，是切向刚度系数。 ，和 分别表示切向位移、阻尼系数和切向摩擦系数。 是接触点的滑动速度：

；

其中，和  是颗粒 和 的半径 和  是颗粒  和  的角速度。 是切向单位矢量，表示如下：

case1: 当时，切向接触力由弹性力和阻尼力决定。这种情况通常发生在接触初期或滑动速度较小的情况下，颗粒间的相对运动较小，摩擦力未达到最大静摩擦力。

case2: 当  时，切向接触力由库仑摩擦定律决定，切向接触力等于摩擦系数与法向接触力的乘积，方向由切向单位矢量决定。这种情况通常发生在滑动速度较大或相对运动较大的情况下，摩擦力达到了最大静摩擦力并转化为动摩擦力。

以下解答一些疑惑：

1. 刚体的角速度：

如果已知刚体上某点的速度 和该点到旋转轴的距离 ，则角速度 可以通过以下公式计算：

其中  是线速度，是该点到旋转轴的矢量， 是该矢量的模。

2. 离散元方法（DEM）中的角速度：

在DEM模拟中，角速度  通过积分旋转运动方程计算。旋转运动方程如下：

其中是刚体的转动惯量， 是作用在刚体上的总转矩。

通过数值积分方法（如欧拉法或龙格-库塔法），可以计算每个时间步上的角速度变化：

3. 颗粒 的接触力 和转矩 ​ 详述如下:

上式中，表示力臂矢量，它是从颗粒质心到接触点的距离矢量。

以下是一个简单的C/C++实现计算的函数：

// 计算切向接触力
Vector3 calculateTangentialForce(const Sphere& sphere1, const Sphere& sphere2, const Vector3& n_ij,
                                 double k_t, double delta_t, double eta_t, double mu, double F_n_ij, double theta) {
    // 计算切向相对速度
    Vector3 u_t_ij = (sphere1.velocity - sphere2.velocity) - (sphere1.velocity.dot(n_ij)) * n_ij
                   + (sphere1.radius * sphere1.angularVelocity + sphere2.radius * sphere2.angularVelocity).cross(n_ij);
    
    // 计算弹性力和阻尼力部分
    Vector3 F_t_ij = -k_t * delta_t * u_t_ij - eta_t * u_t_ij;

    // 判断切向接触力的大小
    if (F_t_ij.norm() <= theta * std::abs(F_n_ij)) {
        return F_t_ij;
    } else {
        return -mu * std::abs(F_n_ij) * u_t_ij.normalize();
    }
}

---

参考文献：

Fo, Bin, et al. "Numerical simulation of fine particle liquid–solid flow in porous media based on LBM‐IBM‐DEM." The Canadian Journal of Chemical Engineering 101.6 (2023): 3576-3591.

孙其诚, and 王光谦. 颗粒物质力学导论. 科学出版社, 2009.

---

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

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