接触力

DEM 中的接触力公式通常是弹簧-阻尼器模型的一种变体。弹簧产生排斥力将颗粒推开,而阻尼器表示粘性阻尼并允许模拟除完全弹性以外的碰撞类型。

接触点上的力建模为一对弹簧-阻尼器振荡器。并行线性弹簧-阻尼器模型表示法向力,用滑块串联的并行线性弹簧-阻尼器表示力的切向(相对接触平面法向矢量)。在这两者中,弹簧考虑响应的弹性部分,而阻尼器则考虑碰撞期间的能量耗散。



Simcenter STAR-CCM+ 提供了三种接触模型:

  • Hertz Mindlin
  • 线性弹簧
  • 沃尔顿布劳恩

Hertz-Mindlin 无滑移接触模型

Hertz-Mindlin 接触模型是非线性弹簧-阻尼器接触模型(基于 Hertz-Mindlin 接触理论 [726]、[722])的一种变体。A 和 B 这两个球体之间的力由以下方程组表示:

图 1. EQUATION_DISPLAY
Fcontact=Fnn+Ftt
(3250)

其中, Fn 和 Ft 分别为法向和切向分量的幅值。

法向由以下各项定义:

  • 法向力:
    图 2. EQUATION_DISPLAY
    Fn=-Kndn-Nnvn
    (3251)
  • 法向弹簧刚度:
    图 3. EQUATION_DISPLAY
    K n = 4 3 E e q d n R e q
    (3252)
  • 法向阻尼:
    图 4. EQUATION_DISPLAY
    Nn=(5KnMeq)Nn damp
    (3253)

    其中,Nn damp 为法向阻尼系数,如 Eqn. (3257) 中给定。

切向由以下各项定义:

  • 如果 |Ktdt|<|Kndn|Cfs ,其中, C f s 为静摩擦系数,则切向力由 - K t d t - N t v t 定义。否则:
    图 5. EQUATION_DISPLAY
    Ft=|Kndn|Cfsdt|dt|
    (3254)
  • 切向弹簧刚度:
    图 6. EQUATION_DISPLAY
    K t = 8 G e q d n R e q
    (3255)
  • 切向阻尼:
    图 7. EQUATION_DISPLAY
    Nt=(5KtMeq)Nt damp
    (3256)

    其中,Nt damp 为切向阻尼系数,如 Eqn. (3258) 中给定。

图 8. EQUATION_DISPLAY
Nn⁢  damp=-ln(Cn rest)π2+⁡ln⁡(Cn rest)2
(3257)
图 9. EQUATION_DISPLAY
Nt  damp=-ln(Ct  rest)π2+⁡ln⁡(Ct  rest)2
(3258)

此处,Cn  rest 和 Ct  rest 为法向和切向恢复系数,属于用户设置的模型参数。如果 Cn  rest=0,则 Nn damp=1;如果 Ct  rest=0,则 Nt damp=1。

等效半径定义为:

图 10. EQUATION_DISPLAY
R e q = 1 1 R A + 1 R B
(3259)

等效颗粒质量为:

图 11. EQUATION_DISPLAY
M e q = 1 1 M A + 1 M B
(3260)

等效杨氏模量表示为:

图 12. EQUATION_DISPLAY
E e q = 1 1 - ν A 2 E A + 1 - ν B 2 E B
(3261)

等效剪切模量为:

图 13. EQUATION_DISPLAY
Geq=12(2-νA)(1+νA)EA+2(2-νB)(1+νB)EB
(3262)

其中:

  • M A 和 M B 为球体 A 和 B 的质量。
  • d n 和 d t 为接触点处的法向和切向重叠。
  • R A 和 R B 为球体的半径。
  • E A 和 E B 为球体的杨氏模量。
  • ν A 和 ν B 为泊松比。
  • ν n 和 ν t 为接触点处相对球体表面速度的法向和切向速度分量。

对于颗粒-壁面碰撞,这些公式保持不变,但假设壁面半径和质量为 R w a l l = ∞ 和 M w a l l = ∞ ,因此等效半径减少至 R e q = R p a r t i c l e ,等效质量减少至 M w a l l = M p a r t i c l e 。

Di Renzo 和 Di Maio [722] 针对切向微滑移和切向力计算的细节处理提出了几个公式。Simcenter STAR-CCM+ 使用 Tsuji 提出的公式 [737],其中假设切向力为非线性,但细节微滑移追踪由解析表达式替换。生成的代码计算效率较高,且仍与实验数据匹配。

临界速度
当颗粒壁面冲击速度超过临界速度时,颗粒有足够的能量穿过壁面边界。Simcenter STAR-CCM+ 计算临界速度的方程如下:
图 14. EQUATION_DISPLAY
v crit = 4 5 E π ρ ( 1 − ν 2 )
(3263)
其中 E 为杨氏模量, ρ 为密度, ν 为泊松比。

线性弹簧接触模型

该模型基于 Cundall 和 Strack 的成果 [720]。法向力和切向力 Fn 和 Ft 如 Eqn. (3251) 和 Eqn. (3254) 中所定义,但是参数 Kn 、 Kt 、 Nn 和 Nt 的定义如下:

  • Kn 为法向弹簧常数。
  • Kt 为切向弹簧常数。
  • Nn=2Nn⁢  dampKnMeq  其中:
    • Nn⁢  damp 为法向阻尼系数,如 Eqn. (3257) 中所示。
    • Meq 为等效颗粒质量,如 Eqn. (3260) 中所示。
  • Nt=2Nt⁢  dampKtMeq  其中, Nt⁢   damp 为切向阻尼系数,如 Eqn. (3258) 中所示。

当线性弹簧模型使用基于颗粒材料方法时,基于线性模型是线性化 Hertz 模型的假设估计弹簧刚度,其中最大法向接触力为:

图 15. EQUATION_DISPLAY
Knδmax⁡=43Eeq⁢ Reqδmax⁡δmax⁡
(3264)

其中:

  • δmax⁡ 为最大重叠。

  • Eeq 为 Eqn. (3261) 中定义的等效杨氏模量。

  • Req 为 Eqn. (3259) 中定义的等效半径。

要计算弹簧常数,假设:

图 16. EQUATION_DISPLAY
λmax⁡=δmax⁡Req⁡
(3265)

由此得出法向弹簧常数 Kn 为:

图 17. EQUATION_DISPLAY
Kn=43λmax⁡Eeq⁢ Req⁢ 
(3266)

切向弹簧常数 Kt 为:

图 18. EQUATION_DISPLAY
Kt=8λmax⁡Geq⁢ Req⁢ 
(3267)

其中, Geq 为等效剪切模量,如 Eqn. (3262) 中定义。

沃尔顿布劳恩迟滞接触模型

沃尔顿布劳恩迟滞塑料弹性接触系统具有以下关系:

图 19. EQUATION_DISPLAY
F n = - K 1 δ n
(3268)

加载接触时,以及:

图 20. EQUATION_DISPLAY
Fn=-K2(δ-δdeformation)n
(3269)

卸载期间,其中:

图 21. EQUATION_DISPLAY
K 1 = 1.6 π R Y 0
(3270)
图 22. EQUATION_DISPLAY
K 2 = 1 E f K 1
(3271)
图 23. EQUATION_DISPLAY
Y0=E⋅YieldStressFraction
(3272)
图 24. EQUATION_DISPLAY
δdeformation=δmax(1-Ef)
(3273)

E 为材料杨氏模量, δ 为接触法向重叠, n 为接触法向矢量,屈服应力分数是一个用户可控模型属性,用于定义塑性变形的开始。 E f 是一个用户可控能量分数,用于定义卸载期间恢复的能量, δ max 为接触加载期间达到的最大重叠。

接触循环

形成颗粒接触后,将在每个时间步追踪接触状态,如下所示:

  • 如果 vr⋅n≥0 且 δ > δ max ,其中, v r 为相对表面速度,则接触视为加载。

    加载接触时,使用 K 1 刚度计算力,并将 δ max 更新为 δ 的值。

  • 如果不加载接触,则为卸载 ( vr⋅n<0 ) 或重新加载 ( δ ≤ δ max )。卸载期间,使用 K 2 刚度计算法向力,但不更新 δ max 的值。
  • 如果颗粒完全分离,则移除颗粒追踪信息并将后续接触视为新接触。

几何解释

加载/卸载循环的几何表示如下: