离散单元法（D iscrete E lement M ethod，简称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
为简单起见，先从仅存在单种球形颗粒的颗粒流开始说起，并且只考虑获取空间内一点 r \boldsymbol{r} r t t t
在这里，引入高斯函数
W ( r ) = { 1 V w exp ( − ∣ r ∣ 2 2 w ) , ∣ r ∣ < c 0 , otherwise W(\boldsymbol{r}) = \begin{cases}
\frac{1}{V_w} \exp{\left(- \frac{\vert\boldsymbol{r}\vert^2}{2w} \right)} ,&
{\vert\boldsymbol{r}\vert} < c \\
0 ,& \text{otherwise}
\end{cases}
W ( r ) = { V w 1 exp ( − 2 w ∣ r ∣ 2 ) , 0 , ∣ r ∣ < c otherwise
其中： w w w c = 3 w c = 3 w c = 3 w V w V_w V w
V w = 2 2 ⋅ w 3 π 3 2 e r f ( c 2 2 w ) − 4 c w 2 π exp ( − c 2 2 w 2 ) V_w = 2 \sqrt{2} \cdot w^3 \pi^{\frac{3}{2}} \mathrm{erf}\left( \frac{c \sqrt{2}}{2 w} \right) - 4 c w^2 \pi \exp{\left( \frac{- c^2}{2 w^2} \right)}
V w = 2 2 ⋅ w 3 π 2 3 e r f ( 2 w c 2 ) − 4 c w 2 π exp ( 2 w 2 − c 2 )
Weinhart 等建议 w w w [ 0.75 d m e a n , 1.25 d m e a n ] [0.75 d_{\mathrm{mean}}, 1.25 d_{\mathrm{mean}}] [ 0 . 7 5 d m e a n , 1 . 2 5 d m e a n ] d m e a n d_{\mathrm{mean}} d m e a n
在 t t t q q q i i i r i \boldsymbol{r}_i r i m i m_i m i u i \boldsymbol{u}_i u i
密度的粗粒化表达式为：
ρ ( r , t ) = ∑ i ∈ q m i W ( r − r i ( t ) ) \rho(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))
ρ ( r , t ) = i ∈ q ∑ m i W ( r − r i ( t ) )
动量的粗粒化表达式为：
j ( r , t ) = ∑ i ∈ q m i u i W ( r − r i ( t ) ) \boldsymbol{j}(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i \boldsymbol{u}_i W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))
j ( r , t ) = i ∈ q ∑ m i u i W ( r − r i ( t ) )
所以，速度的计算需要通过 j \boldsymbol{j} j ρ \rho ρ
u ( r , t ) = j ( r , t ) / ρ ( r , t ) \boldsymbol{u}(\boldsymbol{r}, t) = \boldsymbol{j}(\boldsymbol{r}, t) / \rho(\boldsymbol{r}, t)
u ( r , t ) = j ( r , t ) / ρ ( r , t )
根据 W ( r ) W(\boldsymbol{r}) W ( r ) ∣ r − r i ( t ) ∣ ≥ c \vert \boldsymbol{r} - \boldsymbol{r}_{i}(t) \vert \ge c ∣ r − r i ( t ) ∣ ≥ c W ( r − r i ( t ) ) = 0 W(\boldsymbol{r} - \boldsymbol{r}_{i}(t)) = 0 W ( r − r i ( t ) ) = 0
应力张量 [ σ ] = σ α β [\boldsymbol{\sigma}] = \sigma_{\alpha\beta} [ σ ] = σ α β σ α β ( c ) \sigma_{\alpha\beta}^{(c)} σ α β ( c ) σ α β ( k ) \sigma_{\alpha\beta}^{(k)} σ α β ( k )
σ α β = σ α β ( k ) + σ α β ( c ) \sigma_{\alpha\beta} = \sigma_{\alpha\beta}^{(k)} + \sigma_{\alpha\beta}^{(c)}
σ α β = σ α β ( k ) + σ α β ( c )
对于 σ α β ( k ) \sigma_{\alpha\beta}^{(k)} σ α β ( k )
σ α β ( k ) ( r , t ) = ∑ i ∈ q m i v i , α ′ v i , β ′ W ( r − r i ( t ) ) \sigma_{\alpha\beta}^{(k)}(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i v_{i,\alpha}^{'} v_{i,\beta}^{'} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))
σ α β ( k ) ( r , t ) = i ∈ q ∑ m i v i , α ′ v i , β ′ W ( r − r i ( t ) )
其中 v i , α ′ = u i , α − u α ( r , t ) v_{i,\alpha}^{'} = u_{i,\alpha} - u_{\alpha}(\boldsymbol{r}, t) v i , α ′ = u i , α − u α ( r , t ) i i i
对于 σ α β ( c ) \sigma_{\alpha\beta}^{(c)} σ α β ( c ) i i i
σ α β ( c ) ( r , t ) = ∑ i ∈ q ∑ j ∈ Q ˉ j ≠ i F i j , α l i j , β ∫ 0 1 W ( r − r i ( t ) + s l i j ) d s \sigma_{\alpha\beta}^{(c)}(\boldsymbol{r}, t) = \sum\limits_{i \in q} \sum\limits_{j \in \bar{Q}}^{j \neq i} F_{ij,\alpha} l_{ij,\beta} \int_{0}^{1} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t) + s \boldsymbol{l}_{ij}) \mathrm{d}s
σ α β ( c ) ( r , t ) = i ∈ q ∑ j ∈ Q ˉ ∑ j = i F i j , α l i j , β ∫ 0 1 W ( r − r i ( t ) + s l i j ) d s
式中
j j j Q ˉ \bar{Q} Q ˉ j ≠ i j \neq i j = i i i i j j j F i j \boldsymbol{F}_{ij} F i j i i i j j j F i j , α F_{ij,\alpha} F i j , α α \alpha α l i j \boldsymbol{l}_{ij} l i j i i i j j j i i i j j j l i j , α l_{ij,\alpha} l i j , α α \alpha α
在实际操作中，对于
∫ 0 1 W ( r − r i ( t ) + s l i j ) d s \int_{0}^{1} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t) + s \boldsymbol{l}_{ij}) \mathrm{d}s
∫ 0 1 W ( r − r i ( t ) + s l i j ) d s
我选择使用复化 Simpson 积分 计算其数值，以简化操作。
假设一个颗粒流系统 Q ˉ \bar{Q} Q ˉ Q Q Q q q q q ∈ Q q \in Q q ∈ Q q q q ρ q ( r , t ) \rho^{q}(\boldsymbol{r}, t) ρ q ( r , t ) u q ( r , t ) \boldsymbol{u}^{q}(\boldsymbol{r}, t) u q ( r , t ) σ α β q ( r , t ) \sigma_{\alpha\beta}^{q}(\boldsymbol{r}, t) σ α β q ( r , t )
所以整个颗粒流系统 Q ˉ \bar{Q} Q ˉ
ρ ( r , t ) = ∑ q ∈ Q ρ q ( r , t ) u ( r , t ) = ∑ q ∈ Q u q ( r , t ) σ α β ( r , t ) = ∑ q ∈ Q σ α β q ( r , t ) \begin{aligned}
\rho(\boldsymbol{r}, t) &= \sum_{q \in Q} \rho^{q}(\boldsymbol{r}, t) \\
\boldsymbol{u}(\boldsymbol{r}, t) &= \sum_{q \in Q} \boldsymbol{u}^{q}(\boldsymbol{r}, t) \\
\sigma_{\alpha\beta}(\boldsymbol{r}, t) &= \sum_{q \in Q} \sigma_{\alpha\beta}^{q}(\boldsymbol{r}, t) \\
\end{aligned}
ρ ( r , t ) u ( r , t ) σ α β ( r , t ) = q ∈ Q ∑ ρ q ( r , t ) = q ∈ Q ∑ u q ( r , t ) = q ∈ Q ∑ σ α β q ( r , t )
这里以 μ ( I ) \mu(I) μ ( I )
记速度剪切率张量 γ i j = ∂ u i ∂ x j + ∂ u j ∂ x i \gamma_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} γ i j = ∂ x j ∂ u i + ∂ x i ∂ u j u i u_i u i
γ i j d = γ i j − δ i j 3 t r ( [ γ ] ) \gamma_{ij}^{d} = \gamma_{ij} - \frac{\delta_{ij}}{3} \mathrm{tr}([\gamma])
γ i j d = γ i j − 3 δ i j t r ( [ γ ] )
其中 t r ( [ γ ] ) = γ x x + γ y y + γ z z \mathrm{tr}([\gamma]) = \gamma_{xx} + \gamma_{yy} + \gamma_{zz} t r ( [ γ ] ) = γ x x + γ y y + γ z z [ γ ] [\gamma] [ γ ]
所以剪切应变率定义为：
∣ γ ˙ ∣ = 1 2 γ i j d γ i j d |\dot{\gamma}| = \sqrt{\frac{1}{2} {\gamma}_{ij}^{d} {\gamma}_{ij}^{d} }
∣ γ ˙ ∣ = 2 1 γ i j d γ i j d
定义压强为 p = 1 3 t r ( [ σ ] ) = ( σ x x + σ y y + σ z z ) / 3 p = \frac{1}{3} \mathrm{tr}([\sigma]) = (\sigma_{xx} + \sigma_{yy} + \sigma_{zz}) / 3 p = 3 1 t r ( [ σ ] ) = ( σ x x + σ y y + σ z z ) / 3 [ τ ] = τ i j = σ i j − p ⋅ δ i j [\tau] = \tau_{ij} = \sigma_{ij} - p \cdot \delta_{ij} [ τ ] = τ i j = σ i j − p ⋅ δ i j
∣ τ ∣ = 1 2 τ i j τ i j |\tau| = \sqrt{\frac{1}{2} \tau_{ij} \tau_{ij} }
∣ τ ∣ = 2 1 τ i j τ i j
我用 Python 写了一份对单种颗粒系统 进行粗粒化的代码，详情见此网页 。