碰撞与聚结

Simcenter STAR-CCM+ 使用宿主网格单元或碰撞网格方法检测碰撞。仅当两个颗粒包位于同一网格单元或同一碰撞网格中时才会发生碰撞。Simcenter STAR-CCM+ 用以下碰撞算法中速度较快的一种检测任何特定网格单元中的碰撞。以下算法可用：

• NTC（无计时器）碰撞检测算法 [706]
• 奥罗克算法 [685]

在两个液滴之间检测到碰撞后，可能会产生不同的碰撞结果：反弹、聚结、反射分离和延伸分离（擦边碰撞）[665]：

这种实施会忽略在高碰撞韦伯数下可能形成的卫星液滴。

碰撞结果由三个无量纲参数表示：

• 碰撞韦伯数 ${\text{We}}_{\text{coll}}$ 如下：

  $${\text{We}}_{\text{coll}}=\frac{{\rho }_l{\left(v_{1,2}\right)}^2(r_1+r_2)}{2\sigma }$$

  (3178)

  其中， $\sigma =\frac{1}{2}({\sigma }_1+{\sigma }_2)$ 为液滴 1 和 2 的平均表面张力， ${\rho }_l=\frac{1}{2}({\rho }_1+{\rho }_2)$ 为平均密度。相对速度 $v_{1,2}$ 如下：

  $$v_{1,2}=\left|{\mathbf{v}}_1-{\mathbf{v}}_2\right|$$

  (3179)

  其中， ${\mathbf{v}}_1$ 和 ${\mathbf{v}}_2$ 为碰撞液滴的速度， $r_1$ 和 $r_2$ 为其对应的半径。

  碰撞韦伯数还可以根据较小液滴的直径 $d_s$ 定义，相对于以下公式中给出的 ${\text{We}}_{\text{coll}}$ ：

  $${\text{We}}_d=\frac{{\rho }_l{\left(v_{1,2}\right)}^2{d}_s}{\sigma }=\frac{4\mathrm{\Delta }}{(1+\mathrm{\Delta })}{\text{We}}_{coll}$$

  (3180)
• 冲击参数 $B$ 如下：

  $$B=\frac{b}{r_1+r_2}$$

  (3181)

  其中， $b$ 计算为从一个液滴的中心到放置在另一个液滴的中心上的相对速度矢量的距离。当 $B=0$ 时，碰撞为正面碰撞。当 $B=1$ 时，液滴彼此几乎不擦边。

• 液滴尺寸比率 $\gamma$ 如下：

  $$\gamma =\frac{r_2}{r_1},r_2>r_1$$

  (3182)
  文献中常用的附加尺寸比参数是 $\mathrm{\Delta }=\frac{1}{\gamma },0<\mathrm{\Delta }<1$ 。
• 碰撞 Ohnesorge 数 $Oh$ 定义为：
  $$Oh=\frac{\mu }{\sqrt{\rho \sigma d_s}}$$
  其中 $\mu$ 为碰撞液滴的平均动力粘度。

碰撞结果可以通过碰撞结果图表示，具体取决于影响参数 $B$ 和碰撞韦伯数 $\text{We}$ 。根据它们在碰撞结果图上的位置，液滴相互反弹（简单反弹或反射分离）、聚结或相互擦边（延伸分离）。曲线的常规形状可以随不同尺寸比参数以及液滴粘度等其他物理属性而变化。

Simcenter STAR-CCM+ 中实施了三个碰撞结果图。O'Rourke 图 ([684]、[677])、Ashgriz 图 ([651]) 和复合图。

O'Rourke 和 Ashgiz 图绘制碰撞效率 $E$ 与碰撞韦伯数 ${\text{We}}_{\text{coll}}$ ，通过 Eqn. (3178) 根据平均液滴直径为每对相互作用液滴计算。碰撞效率 $E$ 是由获得给定结果的概率定义的 ( $E=B\cdot B$ )。这两个图中的任何一个均不包含粘性耗散导致的能量损失效应。

O'Rourke 图

O'Rourke 图考虑了反弹、聚结和延伸/擦边分离模态。

Ashgriz 图

Ashgriz 图考虑了聚结、反射分离和延伸/擦边分离模态。

聚结区域的边界由曲线 $E_{\text{coal}}$ 和 $E_{\text{bounce}}$ 或 $E_{\text{graze}}$ 和 $E_{\text{reflex}}$ 表示。

对于每条 $E$ 曲线：

其中 $A$、$a$ 和 ${\text{We}}_{\text{c}}$ 是用于将 $E$ 曲线拟合到数据的经验常数。

这些常数的默认值为：

液滴直径比率修正定义如下：

其中，$\gamma$ 为液滴大小比。常数的默认值为：

• $a_3=2.4$
• $a_2=-5.76$
• $a_1=6.48$
• $a_0=0$

使用上表中列出的 O'Rourke 的默认值，对于碰撞和聚结的这些公式的 Eqn. (3183) 的 $E_{\text{coal}}$ 和 $E_{\text{bounce}}$ 减少：

以下流程图描述了整体碰撞动力学 O'Rourke 算法 [677]：

$Y_{r,2}$ 为 [0,1] 范围内的随机数。

对于 O'Rourke 和 Ashgriz 方法，如果发生聚结，可以通过查找 $m$ 的值来确定每个液滴发生的聚结次数 $m$，且

其中，$P_n$ 为发生 $n$ 次碰撞的概率，由 Eqn. (3214) 给出。

发生聚结时，较稀疏颗粒包中的液滴数量将从较密集的颗粒包中减去。如果结果为擦边碰撞或反弹，则每个液滴仅计算一次碰撞。液滴的碰撞后速度、温度和质量将基于质量、动量和能量守恒进行计算。

复合图

复合图提供三条边界曲线的模型。用于反弹的 Sommerfeld 模型、用于延伸分离的 Suo-Jia 模型和用于反射分离的 Ashgriz 模型。

Sommerfeld

此模型基于 Hu 等方法 [670]，其中在临界碰撞韦伯数中考虑了能量平衡中的粘性耗散效应。

$${\text{We}}_d=\frac{g\left(\mathrm{\Delta }\right)(\mathrm{\Phi }-3)}{(1-B^2)(1-\beta )}$$

(3188)

其中：

• $g\left(\mathrm{\Delta }\right)=\frac{4(1+{\mathrm{\Delta }}^2)(1+{\mathrm{\Delta }}^3)}{{\mathrm{\Delta }}^2}$
• $\beta$ 为能量耗散因子
• $\mathrm{\Phi }$ 为与碰撞过程中液滴变形程度相关的形状因子。

Sommerfeld ([716]、[712]) 方法通过适应具有不同粘度的可用数据，自动适应 $\beta$ 和 $\mathrm{\Phi }$ ，通过 $Oh<0.5$ （对于纯液体）和 $Oh<0.4$ （对于溶液）进行限制。参数 $\beta$ 和 $\mathrm{\Phi }$ 假设为 $B$ 的线性函数。

$$\mathrm{\Phi }\left(B\right)=-k_{\mathrm{\Phi }}Oh+{\mathrm{\Phi }}_{init}Oh$$

(3189)

$$\beta \left(B\right)=(1-{\beta }_{init})Oh+{\beta }_{init}Oh$$

(3190)

$$\begin{array}{c}\begin{array}{c}{k}_{\mathrm{\Phi }}=\max (k_3{Oh}^3+k_2{Oh}^2+k_{1}Oh+k_0,k_l)\\ {\mathrm{\Phi }}_{init}=\max ({\mathrm{\Phi }}_3{Oh}^3+{\mathrm{\Phi }}_2{Oh}^2+{\mathrm{\Phi }}_{1}Oh+{\mathrm{\Phi }}_0,{\mathrm{\Phi }}_l)\\ {\beta }_{init}=\min ({\beta }_3{Oh}^3+{\beta }_2{Oh}^2+{\beta }_{1}Oh+{\beta }_0,{\beta }_l)\end{array}\end{array}$$

(3191)

这些常数的默认值为：

常数	$k_0$	$k_1$	$k_2$	$k_3$	$k_l$
纯液体	0.82	-4.03	10.93	-9.6	0.31
溶液	0.725	-10.2	52.5	-87.7	0.0834

常数	${\mathrm{\Phi }}_0$	${\mathrm{\Phi }}_1$	${\mathrm{\Phi }}_2$	${\mathrm{\Phi }}_3$	${\mathrm{\Phi }}_l$
纯液体	3.9	-4.32	12.4	-11.7	3.4
溶液	3.76	-6.05	18.31	-18.25	3.1

${\beta }_0$	${\beta }_1$	${\beta }_2$	${\beta }_3$	${\beta }_l$
0.07	5.24	-19.6	22.6	0.5

对于迎头碰撞 (B = 0)，临界韦伯数为 ${\text{We}}_d=\frac{g\left(\mathrm{\Delta }\right)(\mathrm{\Phi }-3)}{(1-\beta )}$ 。

Suo-Jia

用于延伸分离的 Suo 和 Jia ([717]) 模型仅包括无量纲参数，由以下公式给出：

$$\sqrt{{\text{We}}_d}=\frac{1}{B}(a+\frac{h}{d_s}Oh)$$

(3192)

其中 a 为 Ohnesorge 数和液滴尺寸比的函数，由以下公式给出：

$$a=(a_0-a_1\mathrm{\Delta })Oh+a_2$$

(3193)

常数的默认值为：

• $a_2=2.58$
• $a_1=180$
• $a_0=183$

此相关性适用于：

$\begin{array}{ccc}\begin{array}{c}1mPas<\mu <47.2mPas\end{array}& for& \mathrm{\Delta }=1\\ 1mPas<\mu <1.2mPas& for& \mathrm{\Delta }<1\end{array}$

Ashgriz

Ashgriz 方法分离临界韦伯数的反射分离和聚结模态。

$${\text{We}}_d=\frac{f\left(\mathrm{\Delta }\right)}{{\mathrm{\Delta }}^6{\eta }_s+{\eta }_l}$$

(3194)

其中：

$$\begin{array}{c}{\eta }_s=2{(1-\xi )}^2{(1-{\xi }^2)}^{\frac{1}{2}}-1\\ {\eta }_l=2{(\mathrm{\Delta }-\xi )}^2{({\mathrm{\Delta }}^2-{\xi }^2)}^{\frac{1}{2}}-{\mathrm{\Delta }}^3\end{array}$$

(3195)

且 $\xi =\frac{1}{2}B(1+\mathrm{\Delta })$

$$f\left(\mathrm{\Delta }\right)=3\mathrm{\Delta }{(1+{\mathrm{\Delta }}^3)}^2[7{(1+{\mathrm{\Delta }}^3)}^{\frac{2}{3}}-4(1+{\mathrm{\Delta }}^2\left)\right]$$

(3196)

对于迎头碰撞 (B = 0)，临界韦伯数为 ${\text{We}}_d=\frac{f\left(\mathrm{\Delta }\right)}{{\mathrm{\Delta }}^3(1+{\mathrm{\Delta }}^3)}$ 。对于大小相等的液滴的迎头碰撞 ${\text{We}}_d=18.67$ 。该数值随着液滴尺寸差异的增加而增加。

碰撞网格生成

碰撞网格的初始尺寸是域中最大网格单元尺寸的两倍。当网格中的颗粒包数增加到足以评估给定碰撞网格中的颗粒平均自由路径时，将对网格进行重构，网格尺寸由最大颗粒平均自由路径确定。

每个碰撞网格生成事件都会围绕域中心生成五个网格大小相同但旋转角度不同的碰撞网格。每个时间级别随机使用这些碰撞网格之一。这样，就可以通过按时间级别重新排列碰撞网格边界来移除碰撞网格边界的障碍效应。旋转碰撞网格示意图如下所示：

本例中的求解域是用手绘封闭曲线绘制的，而碰撞网格边界则是用笛卡尔框绘制的。上图是给定时间级别的配置，下图则是不同时间级别的另一种配置。两种配置的碰撞网格大小相同，但碰撞网格的边界线不同。此方法可以防止碰撞网格的边界成为碰撞的障碍。

并非始终需要在每个时间级别使用不同的碰撞网格大小。可以指定使用不同网格大小重构碰撞网格的时间间隔（更新频率）。在模拟达到下一个更新时间级别之前，它将使用现有的碰撞网格。注意，即使未重构碰撞网格，不同时间级别中使用的网格也不相同；网格大小相同，但旋转角度不同。

对于域连通性或网格单元连通性随时间变化的模拟，例如移动网格或自适应网格细化，即使更新频率不需要，也会重构碰撞网格。

碰撞网格与拉格朗日双向耦合求解器的体积和壳体源项光顺法的网格单元群集完全兼容。它们可以在同一模拟中共存，并且彼此完全独立。

颗粒平均自由路径

颗粒的碰撞频率由 [690] 给出：

$$f_i=\sum _{j\in V_c}\left|{\mathbf{v}}_{ij}\right|\frac{\pi }{4}{(d_i+d_j)}^2\frac{n_j}{V_C}$$

(3215)

其中：

• $V_C$ 为碰撞网格单元的体积。
• $d_i$ 为颗粒 $i$ 的直径。
• ${\mathbf{v}}_{ij}$ 为颗粒 $i$ 和 $j$ 之间的相对颗粒速度。
• $n_i$ 为颗粒包 $i$ 所代表的颗粒数。
• $\sum _{j\in V_C}$ 为碰撞网格中所有颗粒的总和。

使用 Eqn. (3215)，碰撞的时间尺度可以估算如下：

$$\mathrm{\Delta }{t}_i=\frac{1}{f_i}$$

那么颗粒的颗粒平均自由路径可以写为：

$${\lambda }_i=\text{max}\left(\left|{\mathbf{v}}_i\right|,{\left|{\mathbf{v}}_{ij}\right|}_{\text{max}}\right)\mathrm{\Delta }{t}_i$$

(3216)

其中 ${\mathbf{v}}_i$ 为颗粒的速度， $\left|{\mathbf{v}}_{ij}\right|$ 为颗粒 $i$ 与其相邻颗粒（占据同一碰撞网格的颗粒）的最大相对颗粒速度幅值。

使用此公式计算颗粒平均自由路径非常耗时，需要针对给定颗粒包检查位于同一碰撞网格单元中的所有相邻颗粒包，因此在此处使用颗粒平均自由路径的统计估算方法。

在统计估算方法中，颗粒 $i$ 的碰撞频率 $f_i$ 写为：

$$f_i\approx \frac{\pi N_p^c{\sigma }_{v_c}{d}_{\text{avg}}^2{n}_{\text{avg}}}{V_c}$$

(3217)

其中：

• $N_p^C$ 为碰撞网格 $C$ 中的颗粒包数。
• $d_{\text{avg}}$ 为碰撞网格 $C$ 中颗粒直径的算术平均值。
  $$d_{\text{avg}}=\frac{\sum _{i=1}^{N_p^C}{d}_i}{N_p^C}$$
  其中， $d_i$ 为颗粒包 $i$ 的颗粒直径。
• $n_{\text{avg}}$ 为碰撞网格 $C$ 中颗粒数的算术平均值。
  $$n_{\text{avg}}=\frac{\sum _{i=1}^{N_p^C}{n}_i}{N_p^C}$$
  其中， $n_i$ 为颗粒包 $i$ 的颗粒数。
• $V_c$ 为碰撞网格 $C$ 的体积。
• ${\sigma }_{vc}$ 为碰撞网格 $C$ 中颗粒速度的标准偏差。
  $${\sigma }_{vc}=\frac{1}{N_p^C}\sqrt{\sum _{i=1}^{N_p^C}{(v_i-v_{\text{avg}})}^2}$$
  其中， $v_{\text{avg}}$ 为碰撞网格 $C$ 中颗粒速度的算术平均值。

最后，使用所有颗粒包的最大平均自由路径计算碰撞网格尺寸：

$$L_{\text{cm}}^{mfp}=\text{max}\left(\left|{\mathbf{v}}_i\right|\mathrm{\Delta }{t}_i,{\left|{\mathbf{v}}_{ij}\right|}_{\text{max}}\mathrm{\Delta }{t}_i\right)$$

(3218)

其中， ${\left|{\mathbf{v}}_{ij}\right|}_{\text{max}}$ 为颗粒 $i$ 与其搜索范围（即碰撞网格单元）内任何相邻颗粒之间相对速度的最大幅值。

当碰撞网格单元中的颗粒包数（采样数）较小时，Eqn. (3218) 估计的碰撞网格尺寸波动较大。在此类网格中，由于统计采样不佳，网格尺寸往往过大（包裹数密度越小，平均自由路径越大）。为了避免这种情况发生，碰撞网格大小的最终值定义如下：

$$L_{cm}=\underset{k\in C}{\text{max}}\left(L_{cm}^{mfp}\right)$$

(3219)

其中， ${\text{max}}_{k\in C}\left(L\right)$ 表示碰撞网格单元中 $L$ 的最大值，这些网格单元的颗粒包数大于等于碰撞网格节点的目标颗粒包数属性中指定的值。