离散单元法的颗粒数据粗粒化

注: 本文同步发表于 个人主页 和 知乎

离散单元法（Discrete Element Method，简称DEM）是颗粒流数值模拟的重要数值方法，被用于各种各样的颗粒流定量研究中。颗粒流的粗粒化（Coarse Graining），便是将颗粒流流场的大量数据转化为连续场数据的数值方法。

在开始介绍之前，先简要列出主要的参考资料：

• Breard, E. C. P., Dufek, J., Fullard, L., & Carrara, A. (2020). The Basal Friction Coefficient of Granular Flows With and Without Excess Pore Pressure: Implications for Pyroclastic Density Currents, Water‐Rich Debris Flows, and Rock and Submarine Avalanches. Journal of Geophysical Research: Solid Earth, 125(12), e2020JB020203. DOI:10.1029/2020JB020203.
• Lacaze, L., & Kerswell, R. R. (2009). Axisymmetric Granular Collapse: A Transient 3D Flow Test of Viscoplasticity. Physical Review Letters, 102(10), 108305. DOI:10.1103/PhysRevLett.102.108305.
• Weinhart, T., Labra, C., Luding, S., & Ooi, J. Y. (2016). Influence of coarse‐graining parameters on the analysis of DEM simulations of silo flow. Powder Technology, 293, 138–148. DOI:10.1016/j.powtec.2015.11.052

单种颗粒的粗粒化

为简单起见，先从仅存在单种球形颗粒的颗粒流开始说起，并且只考虑获取空间内一点 $\mathbf{r}$ 在 $t$ 时刻的粗粒化数据。

粗粒化函数

在这里，引入高斯函数

$$W(\mathbf{r})=\{\begin{array}{ll}\frac{1}{V_w}\exp \left(-\frac{|\mathbf{r}{|}^2}{2w}\right),& |\mathbf{r}|<c\\ 0,& \text{otherwise}\end{array}$$

其中： $w$ 为粗粒化宽度； $c=3w$ 为粗粒化截断长度； $V_w$ 是确保其密度积分等于总质量的系数

$$V_w=2\sqrt{2}\cdot w^3{\pi }^{\frac{3}{2}}\mathrm{erf}\left(\frac{c\sqrt{2}}{2w}\right)-4c{w}^2\pi \exp \left(\frac{-c^2}{2w^2}\right)$$

Weinhart 等建议 $w$ 的取值范围应在 $[0.75d_{\mathrm{mean}},1.25d_{\mathrm{mean}}]$ 之间，其中 $d_{\mathrm{mean}}$ 为平均颗粒直径。

密度和动量的计算

在 $t$ 时刻，单种颗粒的集合为 $q$ 。对于其中的第 $i$ 个球，其位置矢量为 ${\mathbf{r}}_i$ ，质量为 $m_i$ ，速度为 ${\mathbf{u}}_i$。

密度的粗粒化表达式为：
$$\rho (\mathbf{r},t)=\sum _{i\in q}{m}_{i}W(\mathbf{r}-{\mathbf{r}}_i(t))$$

动量的粗粒化表达式为：
$$\mathbf{j}(\mathbf{r},t)=\sum _{i\in q}{m}_i{\mathbf{u}}_{i}W(\mathbf{r}-{\mathbf{r}}_i(t))$$

所以，速度的计算需要通过 $\mathbf{j}$ 和 $\rho$ 完成：

$$\mathbf{u}(\mathbf{r},t)=\mathbf{j}(\mathbf{r},t)/\rho (\mathbf{r},t)$$

根据 $W(\mathbf{r})$ 的定义，当 $|\mathbf{r}-{\mathbf{r}}_i(t)|\ge c$ 时， $W(\mathbf{r}-{\mathbf{r}}_i(t))=0$。
因此在实际计算中可根据这一点跳过某些球的计算。

应力张量的计算

应力张量 $[\mathbf{\sigma }]={\sigma }_{\alpha \beta }$ 被表示为由接触作用产生的部分 ${\sigma }_{\alpha \beta }^{(c)}$ ，以及由动能作用产生的部分 ${\sigma }_{\alpha \beta }^{(k)}$ 。即：

$${\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^{(k)}+{\sigma }_{\alpha \beta }^{(c)}$$

对于 ${\sigma }_{\alpha \beta }^{(k)}$ ，有：

$${\sigma }_{\alpha \beta }^{(k)}(\mathbf{r},t)=\sum _{i\in q}{m}_i{v}_{i,\alpha }^{{}^{\prime }}{v}_{i,\beta }^{{}^{\prime }}W(\mathbf{r}-{\mathbf{r}}_i(t))$$

其中 $v_{i,\alpha }^{{}^{\prime }}=u_{i,\alpha }-u_{\alpha }(\mathbf{r},t)$ 为第 $i$ 个球的速度波动值。

对于 ${\sigma }_{\alpha \beta }^{(c)}$ ，考虑第 $i$ 个球的所有接触：

$${\sigma }_{\alpha \beta }^{(c)}(\mathbf{r},t)=\sum _{i\in q}\sum _{j\in \bar{Q}}^{j\ne i}{F}_{ij,\alpha }{l}_{ij,\beta }{\int }_0^{1}W(\mathbf{r}-{\mathbf{r}}_i(t)+s{\mathbf{l}}_{ij})\mathrm{d}s$$

式中

• $j$ 是整个颗粒系统 $\bar{Q}$ 中的一个球（$j\ne i$），并且球 $i$ 与球 $j$ 发生接触；
• ${\mathbf{F}}_{ij}$ 是球 $i$ 所受到的来自球 $j$ 接触的力， $F_{ij,\alpha }$ 为其 $\alpha$ 方向分量；
• ${\mathbf{l}}_{ij}$ 是球 $i$ 与球 $j$ 的接触矢量（从 $i$ 向 $j$）， $l_{ij,\alpha }$ 为其 $\alpha$ 方向分量。

在实际操作中，对于
$${\int }_0^{1}W(\mathbf{r}-{\mathbf{r}}_i(t)+s{\mathbf{l}}_{ij})\mathrm{d}s$$
我选择使用复化 Simpson 积分计算其数值，以简化操作。

多种颗粒的粗粒化

假设一个颗粒流系统 $\bar{Q}$ 中有 $Q$ 类颗粒，对于第 $q$ 类颗粒 $q\in Q$， 可以使用单类颗粒粗粒化的公式计算出：第 $q$ 类颗粒的粗粒化密度 ${\rho }^q(\mathbf{r},t)$ 、速度 ${\mathbf{u}}^q(\mathbf{r},t)$ 和应力张量 ${\sigma }_{\alpha \beta }^q(\mathbf{r},t)$。

所以整个颗粒流系统 $\bar{Q}$ 的粗粒化结果为：

$$\begin{array}{rl}\rho (\mathbf{r},t)& =\sum _{q\in Q}{\rho }^q(\mathbf{r},t)\\ \mathbf{u}(\mathbf{r},t)& =\sum _{q\in Q}{\mathbf{u}}^q(\mathbf{r},t)\\ {\sigma }_{\alpha \beta }(\mathbf{r},t)& =\sum _{q\in Q}{\sigma }_{\alpha \beta }^q(\mathbf{r},t)\end{array}$$

使用粗粒化数据提取颗粒流流变参数

这里以 $\mu (I)$ 流变关系为例子，以上文的文献作为参照。

记速度剪切率张量 ${\gamma }_{ij}=\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}$ ，其中 $u_i$ 源自粗粒化速度场。令：

$${\gamma }_{ij}^d={\gamma }_{ij}-\frac{{\delta }_{ij}}{3}\mathrm{tr}([\gamma ])$$

其中 $\mathrm{tr}([\gamma ])={\gamma }_{xx}+{\gamma }_{yy}+{\gamma }_{zz}$ 为 $[\gamma ]$ 的迹。

所以剪切应变率定义为：

$$|\dot{\gamma }|=\sqrt{\frac{1}{2}{\gamma }_{ij}^d{\gamma }_{ij}^d}$$

定义压强为 $p=\frac{1}{3}\mathrm{tr}([\sigma ])=({\sigma }_{xx}+{\sigma }_{yy}+{\sigma }_{zz})/3$，偏应力张量为： $[\tau ]={\tau }_{ij}={\sigma }_{ij}-p\cdot {\delta }_{ij}$， 所以切应力为：
$$|\tau |=\sqrt{\frac{1}{2}{\tau }_{ij}{\tau }_{ij}}$$

简单的粗粒化代码示意

我用 Python 写了一份对单种颗粒系统进行粗粒化的代码，详情见此网页。