A two-scale coarse-grained discrete element methodand experimental verification of bulk coal
-
摘要:
离散元法（DEM）常用于工程尺度下颗粒系统的模拟，但计算效率仍是大规模颗粒系统运行的制约因素。现有的粗粒化处理方法适用性有限且缺乏普适的理论基础。为解决此问题，研究利用量纲分析描述精确缩尺系统中物理量的缩放定律，通过代表性体积单元 (RVE)，在粗粒化系统和原始系统之间建立了质量、动量和能量的近似守恒关系，并获得了2个不同尺度（宏观和细观）下相应物理量的缩放关系。为验证所提出的双尺度粗粒化离散元方法的正确性，将该方法应用于煤散料下落和回转试验中，真实试验中煤颗粒尺寸为4 mm，仿真中采用5组不同缩放系数的煤颗粒，其尺寸范围为4～12 mm，试验与仿真中煤散料总体积均为0.001 m3。煤散料下落试验采用平均冲击力和堆积休止角作为对比指标，结果表明，随着煤颗粒尺寸增大，两者的相对误差总体呈现增大趋势（尺寸为12 mm时，冲击力误差为14.36%，休止角误差为19.05%）。煤散料回转试验采用上试样受力和煤散料堆积轮廓曲线相关系数作为对比指标，结果表明，随着颗粒尺寸增大，相对误差随之增大，相关系数总体降低（尺寸为12 mm时，受力误差为39.29%，相关系数为0.957 4）。仿真所需时长比预测时长还要少，结果充分证明了双尺度粗粒化离散元法的有效性，粗粒化系统在可接受的误差范围内可显著提高计算效率。
Abstract:The discrete element method (DEM) is commonly used for simulating particle systems at the engineering scale, but computational efficiency remains a limiting factor for large-scale particle system simulations. Existing coarse-grained methods have limited applicability and lack a universal theoretical basis. To address this issue, this study utilizes dimensional analysis to describe the scaling laws of physical quantities in an exact scaled system. By using representative volume elements (RVE), approximate conservation relationships for mass, momentum, and energy are established between the coarse-grained system and the original system. Scaling relationships for corresponding physical quantities at two different scales (global and particle level) are obtained.To validate the correctness of the proposed two-scale coarse-grained DEM method, the method is applied to the falling and rotation tests of bulk coal. The particle size in the actual tests is 4mm, while in the simulations, five groups of coal particles with different scaling coefficients ranging from 4 mm to 12 mm are used. The total volume of bulk coal in both the tests and simulations is 0.001 m3. The falling test of bulk coal uses the average impact force and the reposeangle as comparative indicators. The results show that as the particle size of coal increases, the relative errors of both indicators generally increase (at a particle size of 12 mm, the error in impact force is 14.36%, and the error in the reposeangle is 19.05%).The rotational test of bulk coal uses the correlation coefficient between the applied force on the upper sample and the profile curve of the bulk coal heap as a comparative indicator. The results show that as the particle size increases, the relative error increases, and the correlation coefficient generally decreases (at a particle size of 12 mm, the force error is 39.29%, and the correlation coefficient is
0.9574 ).The simulations require less time compared to the predicted duration, fully demonstrating the effectiveness of the dual-scale coarse-graining DEM method. The coarse-grained system can significantly improve computational efficiency within an acceptable range of errors.-
Keywords:
- two-scale coarse-grained /
- exact scaling /
- DEM /
- coal falling test /
- coal rotation test
-
0. 引 言
离散元法广泛应用于不同领域的颗粒力学行为模拟，包括化学与制药工业[1]、农业[2]以及土木工程[3]等。在煤炭开采[4]、运输[5]与筛分[6]过程中，离散元法被用于分析煤散料与刚性结构间的相互作用，优化机械参数从而提高作业性能。煤散料形状不规则且数量巨大，模拟实际尺寸的煤颗粒系统在计算上非常耗时。普通计算机硬件无法满足计算需求，当前的离散元模拟颗粒数量最多可以达到数千万的水平[7]，然而，工业系统通常拥有数万亿颗粒，即使采用GPU加速也无法满足模拟要求。因此，在工业尺度上使用离散元法模拟大规模煤散料系统仍是一个挑战。煤散料的堆积与运输在煤炭工业现场十分常见，以往学者在进行煤散料系统模拟时，通常使用比实际尺寸大得多的煤颗粒，却并没有明确说明颗粒尺寸选择的依据。
近年来，研究人员在粗粒化理论方面进行了较多研究，以期解决在工程应用中大规模颗粒系统计算的难题。有学者[8-9]提出相似颗粒组合（SPA）模型，并被用于大规模离散元法模拟。在该方法中，具有相似物理或化学特性的颗粒被放大h倍得到粗粒化颗粒。该方法假设粗粒化系统中颗粒的排列方式与原始系统相似，然而，该方法对于控制方程的缩放定律缺乏严格的理论推导，通过假定颗粒大小在颗粒动力行为中起决定性作用，直接将h3作为接触力、液桥力及拖曳力的放大系数。
粗粒化方法更加关注由大多数颗粒引起的集体力学行为，而非个体颗粒轨迹。SAKAI等[10-12]提出了粗粒化CG模型（Corse Grain Model），其中粗粒化系统中的颗粒尺寸是原始系统中颗粒尺寸的l倍。该方法假设粗粒化颗粒的动能与原始颗粒动能相同，粗粒化颗粒的平移和旋转运动假设与原始颗粒相同。在粗粒化颗粒的二元碰撞中，作用在颗粒上的接触力为原始颗粒的l3倍。粗粒化系统中的拖曳力和重力按原始颗粒系统的l3倍进行缩放。对于范德华力，粗粒化模型根据能量守恒原理得到其缩放系数为l2。
BIERWISCH等[13]提出了一种利用量纲分析法构建粗粒化系统的方法，通过构建具有缩放颗粒尺寸的有效介质，使其能够包含与原始系统未缩放颗粒相同的能量密度和能量密度演化。一些学者[14-15]提出了用于模拟颗粒剪切流动的粗粒化方法CGSF（Corse-grained method for granular shear flow），其中滑动摩擦因数、线性刚度系数和恢复系数被缩放以满足能量守恒关系。然而，该方法缩放定律复杂，且模型仅适用于颗粒剪切混合运动，并不适用于所有情况。LU等[16]提出了一种基于能量最小化多尺度EMMS（energy-minimization multi-scale）的粗粒化模型，它的计算速度比普通离散元计算快几个数量级，并且能充分利用GPU。采用两相分解的EMMS模型确定粗粒化颗粒的尺寸及固体浓度，以及它们与气流的相互作用力，根据颗粒流动力学理论确定了粗粒化颗粒之间的相互作用。CHU等[17]在粗粒化系统与真实系统的总能量（包括势能和动能）保持一致的条件下建立了CFD-DEM耦合下的粗粒化模型用于模拟重介质旋流器中的旋转多相流，通过与原始颗粒的模拟结果进行比较，发现在流体密度较大时，两模型之间的差异很小。Tausendschön等[18]研究了基于包块的CFD-DPM流化气固流模拟中颗粒间黏聚力和颗粒间液体传递的粗粒化方法。考虑了4种不同的粗粒化策略，并通过量纲分析确定了如何缩放与范德华内聚力及颗粒间液桥相关的参数。基于无量纲重叠、应力和有效恢复系数策略产生的粗粒化规则相同，而基于匹配粘结键数的策略则产生不同的规则。程宏旸等[19]利用粗粒化（CG）推导有限元－离散元（FEM-DEM）表面和体积耦合的一般性表达式，对于表面耦合，CG将耦合力分布到颗粒−单元接触点以外的位置，如相邻的积分点；对于体积耦合，CG将颗粒尺度的运动均匀化到耦合单元上。
FENG等[20-21]从普适性角度出发，提出了精确缩尺系统，并推导出精确缩尺系统中颗粒集合各物理量之间应满足的缩放关系。利用该模型，赵婷婷等[22]通过多尺度方法建立了粗粒化系统与原始系统之间的缩放关系，并得到了离散元接触模型中相关参数的缩放定律，但该文仅考虑了圆球形颗粒，且没有与真实试验进行对比。
研究旨在基于精确缩尺系统，建立粗粒化系统与原始系统之间的双尺度缩放关系，并得出离散元接触模型中相关参数的缩放定律，将该方法扩展到三维不规则颗粒系统，为验证粗粒化方法对形状不规则煤散料的适用性，进行了煤散料下落试验以及回转试验，在离散元仿真中采用不同缩放系数的煤颗粒进行了同等条件下的模拟。
1. 双尺度粗粒化方法
1.1 精确缩尺系统
图1描述了原始系统、精确缩尺系统和粗粒化系统之间的关系，其中原始系统为一个由4×4×4个球形颗粒组成的规则排列，粗粒化系统占据与原始系统相同的几何域。在精确缩尺系统中，相比原始系统，颗粒尺寸和几何域均按照2倍比例放大。这种精确缩尺使我们能够得到原始系统和缩放颗粒组之间不同物理量的比例关系。然而，精确缩尺同时也会对计算域进行缩放，导致精确缩尺系统中的颗粒数与原始系统中颗粒数相同，无法提高计算效率。因此，在工业规模的颗粒系统模拟中，必须使用粗粒化系统。在推导粗粒化系统的缩放定律前，首先需要了解精确缩尺系统的缩放定律，其主要的推导过程可参考这2篇文献[20-21]，接下来将会简要介绍其推导过程。
1.2 精确缩尺系统中颗粒物理量的缩放
若只考虑颗粒的机械运动，那么任何物理量q的大小均可由国际单位制（SI）中3个基本变量长度[L]、质量[M]和时间[T]的组合推导出来。
然而，对于文献[20-21]中描述的离散元模型，选择的3个基本单位变量为：u1—长度[L]，u2—时间[T]和u3—密度[ρ]，其中只有u3与SI单位的选择不同。密度的量纲可以表示为[ρ]=ML−3。在SI单位下，基本单位组合{L, T, ρ}和基本单位组合{L, M, T}的转换矩阵及其逆矩阵为：
$$ {{\boldsymbol{R} = }}\left[ {\begin{array}{*{20}{c}} 1&0&{ - 3} \\ 0&0&1 \\ 0&1&0 \end{array}} \right] , {{{{\boldsymbol{R}}}}^{ - 1}}{\text{ = }}\left[ {\begin{array}{*{20}{c}} 1&3&0 \\ 0&0&1 \\ 0&1&0 \end{array}} \right] $$ (1) 在精确缩尺系统中，通常可以独立选择与3个基本变量相对应的缩放系数，分别为长度h，时间h，密度1，即基本单元转换系数为Hb={h, h, 1}。然后，与精确缩尺系统中任意物理量q对应的缩放系数hq可从基本单元转换系数、转换矩阵的逆R−1推导出来：
$$ {h}_{q}={{\boldsymbol{H}}}_{b}\widehat{·}\langle {q}^{\text{'}}\rangle \equiv {h}_{1}^{{a}_{1}^{\text{'}}}{h}_{2}^{{a}_{2}^{\text{'}}}{h}_{3}^{{a}_{3}^{\text{'}}}={\displaystyle \prod _{i=1}^{3}{h}_{i}^{{a}_{i}^{\text{'}}}}={{\boldsymbol{H}}}_{b}\widehat{·}\left({{\boldsymbol{R}}}^{-1}\langle q\rangle \right) $$ (2) 其中，$ \left\langle q \right\rangle $=[a1, a2, a3]T为精确缩尺系统中量纲[q]的向量表示，即$\left[ q \right]{\rm{ = }}{L^{{{\rm{a}}_{\rm{1}}}}}{M^{{{\rm{a}}_{\rm{2}}}}}{T^{{{\rm{a}}_{\rm{3}}}}} $。
以物理量力F为例，在由SI单位表示的粒子系统中，其量纲为
$$ [F]=N=LMT^{-2}, \left\langle F \right\rangle = {\left[ {1,1, - 2} \right]^{\mathrm{T}}} $$ (3) 在基本单位为{L, T, ρ}的精确缩尺系统中，力F的量纲由以下向量表示：
$$ \left\langle {{F'}} \right\rangle = {R^{ - 1}}\left\langle F \right\rangle = {\left[ {4, - 2,1} \right]^{\mathrm{T}}} $$ (4) 根据公式(2)，可以得到精确缩尺系统中力F的缩放系数hF，计算公式如下：
$$ {h}_{{\mathrm{F}}}={{\boldsymbol{H}}}_{{\mathrm{b}}}\widehat{·}\langle {F}^{\text{'}}\rangle ={h}^{4}·{h}^{-2}·{1}^{1}={h}^{2} $$ (5) 根据上述推导过程，可以计算精确缩尺系统中每个物理量的缩放系数。表1列出了一些主要的物理量，而其他物理量可以在文献[21]中找到。选择基本缩放系数集合Hb={h, h, 1}可确保精确缩尺系统中的应力、应变、杨氏模量和能量相关性质与原始系统中的相等。
表 1 精确缩尺模型中主要物理量的缩放系数Table 1. Scaling factors of main physical quantities in exact scaled model物理量 符号 量纲 缩放系数 长度 L [L] h 时间 t [T] h 密度 ρ [ρ] 1 体积 V [L]3 h3 质量 M [L]3[ ρ] h3 力 F [L]4[T]−2[ρ] h2 应力 σ [L]2[T]−2[ ρ] 1 应变 ε [1] 1 弹性模量 E [L]2[T]−2[ ρ] 1 动能密度 ek [L]2[T]−2[ ρ] 1 力矩 MO [L]5[T]−2[ ρ] h3 从表1可看出，在精确缩尺模型中，时间变量比原始系统大h倍，导致相应的时步放大h倍。因此，当采用中心差分方法进行时间积分时，使用精确缩尺系统并不能提高整体的计算效率，因为两个系统所需的时步数保持不变。另外，对于无量纲系数如摩擦因数、泊松比和阻尼系数等，不需要进行缩放处理。
1.3 粗粒化模型
粗粒化系统与原始系统所占据的几何区域相同。这样可以减少粒子数量，同时仍保留在精确缩尺系统中获得的整体相似规律。以图2中立方体内包含的颗粒集合作为代表性体积单元（RVE）。为确保粗粒化系统中的离散元计算结果与原始系统的物理力学性质相一致，RVE在两个系统中需要满足一定条件。具体来说，RVE必须满足几何一致，并且在质量、动量和能量密度方面满足近似守恒。
![]()
图 2 原始系统及粗粒化系统中的代表性单元Figure 2. Representative units in original system and coarse-grained system假设粗粒化系统和原始系统之间存在几何缩放系数h，那么粗粒化RVE中的粒子数$ \overline N $和接触数$ {\overline N _{\rm{c}}} $与原始RVE中的粒子数N和接触数Nc之间的关系满足：
$$ N \approx {h^3}\overline N \text{，} {N_{\rm{c}}} \approx {h^3}{\overline N _{\rm{c}}} $$ (6) RVE边界上的颗粒接触数满足下式，其中$ {N_{\rm{b}}} $和$ {\overline N _{\rm{b}}} $分别为原始RVE和粗粒化RVE边界上的颗粒接触数：
$$ {N_{\rm{b}}} \approx {h^2}{\overline N _{\rm{b}}} $$ (7) 由于满足质量守恒条件，单个粗粒化粒子的质量$ \overline m $可以与原始系统中的粒子质量m相关联，其关系为
$$ \overline m = {h^3}m $$ (8) 动量守恒条件要求：
$$ \sum\limits_{i = 1}^N {{m_i}{v_i}} \approx \sum\limits_{i = 1}^{\overline N } {{{\overline m }_i}{{\overline v }_i}} $$ (9) 能量密度守恒包括动能密度守恒、应变能密度守恒和能量耗散率守恒。2个RVE粒子系统的平均柯西应力表达式应大致相同，如下所示，其中$ {x_k} $为原始系统中的接触点位置，$ {\overline x _k} $为粗粒化系统中的接触点位置：
$$ \sigma = \frac{1}{{\left| \varOmega \right|}}\sum\limits_{k = 1}^{{N_b}} {{f_k} \otimes {x_k} \approx } \overline \sigma = \frac{1}{{\left| \varOmega \right|}}\sum\limits_{k = 1}^{{{\overline N }_b}} {{{\overline f }_k} \otimes {{\overline x }_k}} $$ (10) 根据方程（7），可推导出两个系统中RVE边界上接触力的缩放关系：
$$ {\overline f _k} \approx {h^2}{f_k} $$ (11) 原始系统中RVE的整体平动方程为：
$$ \begin{split} &\sum\limits_{i = 1}^N {m{}_i{a_i}} = \sum\limits_{i = 1}^N {\left( {\sum\limits_j {f_{ij}^{\rm{c}} + {m_i}g + f_i^{\rm{a}} + f_i^{\rm{f}}} } \right)} = \\ &\qquad\sum\limits_{k = 1}^{{N_{\rm{b}}}} {{f_{\rm{k}}}} + \sum\limits_{i = 1}^N {\left( {{m_i}g + f_i^{\rm{a}} + f_i^{\rm{f}}} \right)} \end{split} $$ (12) 粗粒化系统中RVE的整体平动方程为，其中$ \overline f _{ij}^c $为颗粒间的接触力：
$$ \begin{split} &\sum\limits_{i = 1}^N {\overline m {}_i{{\overline a }_i}} = \sum\limits_{i = 1}^N {\left( {\sum\limits_j {\overline f _{ij}^{\rm{c}} + {{\overline m }_i}g + \overline f _i^{\rm{a}} + \overline f _i^{\rm{f}}} } \right)} = \\ &\qquad\sum\limits_{k = 1}^{{{\overline N }_{\rm{b}}}} {{f_{\rm{k}}}} + \sum\limits_{i = 1}^{\overline N } {\left( {{{\overline m }_i}g + \overline f _i^{\rm{a}} + \overline f _i^{\rm{f}}} \right)} \end{split} $$ (13) 根据式（9），可得到$ \displaystyle\sum\limits_{i = 1}^N {{m_i}{a_i}} \approx \displaystyle\sum\limits_{i = 1}^{\overline N } {{{\overline m }_i}{{\overline a }_i}} $。根据式（6）和式（7），控制方程中各力分量的缩放系数不同：细观尺度接触力的缩放系数是h2，即$ {\overline f _k} \approx {h^2}{f_k} $；宏观尺度力（包括重力、外加力和拖曳力）的缩放系数是h3，即$ {\overline m _i}g = {h^3}{m_i}g $, $ \overline f _i^{\mathrm{a}} = {h^3}f_i^{\mathrm{a}} $, $ \overline f _i^{\mathrm{f}} = {h^3}f_i^{\mathrm{f}} $。
从以上分析可以看出，在精确缩尺系统中，接触力的缩放系数与粗粒化系统完全相同。因此，精确缩尺模型中提出的对于不同类型的离散元接触模型的缩放定律可以完全应用于粗粒化系统中的离散元计算。
所提出的“双尺度”为原始系统和粗粒化系统之间存在的两种尺度的缩放关系，粒子间接触力代表了“细观”尺度力，重力、外加力、拖曳力等代表“宏观”尺度力。在粗粒化系统中，细观尺度力的缩放关系与精确缩尺模型的缩放关系一致（h2），而宏观尺度力的缩放系数为h3。
对于时间变量，与细观尺度相关的量（例如时间步长）的缩放系数是h，与宏观尺度相关的量（例如总模拟时间）保持不变。这样就形成了一个双尺度的粗粒化框架。在计算成本方面，将粒子尺寸扩大h倍将使粒子数量减少h3倍，从而使每个时间步的计算效率提高相同的倍数。此外，由于时间步长也扩大了h倍，时间积分步数的总数减少了h倍。因此，粗粒化系统的总时间理论上是原始系统的1/h4。
1.4 煤散料模型
为验证以上双尺度粗粒化方法的正确性，进行了煤散料下落和回转试验，这些试验常用于校准和修正与冲击和输运相关的煤散料参数[23-24]。在仿真过程中，粒子集合的总体积保持不变，只改变煤颗粒尺寸。煤散料的不规则形状用具有相同半径的圆球形颗粒拼接表示，在仿真中考虑采用3种不同形状的煤颗粒组合用于表示真实煤散料，包括扁平状、类锥状和类块状，如图3所示。颗粒组合的最大长度用于表示单个煤颗粒尺寸。仿真中采用5组不同缩放系数的颗粒，见表2，以C-1为例，煤颗粒模型的最长边尺寸设定为4 mm，总体积为0.001 m3，每种煤颗粒形状占1/3（质量分数）。在其他组的模拟中，保证煤散料总体积为0.001 m3不变，C-2到C-5中的煤颗粒通过将C-1中的等比例放大而得（放大系数分别为1.5、2、2.5和3）。颗粒尺寸的选取主要是依据进行试验的设备尺寸来定的，颗粒尺寸过大会导致颗粒无法通过相应设备，且颗粒过少无法形成堆积，而颗粒尺寸过小会导致仿真时间急剧增加。
根据笔者课题组对中国陕西省榆林市神木县长焰煤的标定结果，设置了相应的参数见表3[25]。使用Hertz-Mindlin（no slip）模型来模拟煤散料之间以及煤散料与刮板输送机之间的相互作用。从文献[21]可知，在三维离散元计算中，赫兹模型是尺度无关模型，原始系统中使用的接触参数可应用于不同尺度的粗粒化系统。
表 2 仿真计算条件Table 2. Calculation conditions参数 C-1 C-2 C-3 C-4 C-5 仿真煤颗粒尺寸/mm 4 6 8 10 12 缩放系数 1 1.5 2 2.5 3 表 3 验证试验所需仿真参数Table 3. Simulation input parameters for verification test材料本征参数 数值 接触参数 数值 煤剪切模量 /Pa 4.7×108 煤−煤恢复因数 6.4×10−1 煤泊松比 3×10−1 煤−煤静摩擦因数 3.29×10−1 煤密度/（kg·m−3） 1.229×103 煤−煤滚动摩擦因数 3.6×10−2 钢剪切模量 /Pa 8×1010 煤−钢恢复因数 6.5×10−1 钢泊松比 3×10−1 煤−钢静摩擦因数 4.6×10−1 钢密度/（kg·m−3） 75×103 煤−钢滚动摩擦因数 3.2×10−2 2. 煤散料下落试验
2.1 煤散料下落试验装置及仿真模型
研究使用的煤散料下落试验装置如图4a所示，使用了来自神木县的长焰煤进行试验。试验过程如图4b所示。压力传感器（DYZ-102）安装在底板下方，并与测力仪表（D800）相连接，其由24V直流供电，采用RS485标准与计算机进行通信（图4c）。压力传感器用于测量煤散料下落到底板上的冲击力，而高速摄像机（GC-P100）用于记录煤散料流动形态及试验结束时煤散料堆积形成的休止角。底板上的平均冲击力和休止角度作为表征参数以验证粗粒化离散元模型的准确性。试验进行了3次，取平均值作为最终测试结果。
仿真中的煤下落试验装置三维模型与真实试验相一致（图4d），煤散料参数与表3中的参数一致，使用EDEM软件进行离散元模拟。料斗中的颗粒工厂生成速率设置为5 kg/s，料斗下方的上挡板在t=0.4 s时（煤散料稳定后）以10 m/s的速度被水平抽出。表2中显示的煤颗粒由料斗中的颗粒工厂生成。将从仿真中获得的平均冲击力和休止角与实际试验结果进行比较，以验证粗粒化离散元模型的准确性。
2.2 仿真与试验结果对比
1）不同时刻煤散料流动形态对比。在试验过程中，使用高速摄像机记录了4 mm煤颗粒的流动形态。完全脱离漏斗的时刻被视为时刻0。选择了0.1～0.6 s不同时刻的煤流动形态，并将其与相同尺寸颗粒仿真中对应时刻进行比较。如图5所示，可以观察到试验和仿真中相同时刻的煤流动形态非常相似。
2）冲击力对比。图6a所示为使用4 mm煤散料进行真实试验和仿真的冲击力对比。在底板上的最大冲击力发生在煤落下的初始阶段，随后波动逐渐减小直至冲击结束。为比较平均冲击力，对两条曲线进行积分，并除以积分区间的长度得到相应的平均值。真实试验测试中的平均冲击力为5.855 N，而仿真结果显示平均冲击力为6.066 N。仿真和试验值之间的相对误差为3.60%，两者结果非常接近。从图6还可发现，随着颗粒尺寸的增大，冲击力的波动性变得愈加严重，这是煤散料颗粒尺寸增加导致下落过程中颗粒离散程度变大，从而使底板受到不连续冲击造成的。从图6中可发现，冲击结束后底板上的真实作用力大于仿真中的，这是由于仿真中使用的煤散料由球形颗粒拼接而成，因此会有更好的流动性，而真实试验中的煤散料形状不规则，流动性偏差，导致底板上堆积的煤散料更多。详细的试验与仿真结果见表4，总体来讲，随着颗粒尺寸的增加，真实试验与仿真之间的平均冲击力相对误差增加，当仿真中煤颗粒尺寸为12 mm时，误差达到最大（14.36%）。
表 4 煤散料下落试验的仿真条件与结果Table 4. Calculation conditions and results of coal falling test仿真条件 煤颗粒
尺寸/mm缩放系数 颗粒数 计算时长 /h 平均冲击力 /N 平均冲击力的
相对误差休止角/(°) 休止角的
相对误差/%真实试验 4 — — — 5.855 — 28.5285 — C-1 4 1 49183 91.18 6.066 3.60% 30.6340 7.38 C-2 6 1.5 13676 11.83 5.986 2.24% 29.2924 2.68 C-3 8 2 5548 2.67 5.171 11.68% 32.6396 14.41 C-4 10 2.5 2687 0.98 5.240 10.50% 24.3230 14.74 C-5 12 3 1509 0.40 5.014 14.36% 23.0934 19.05 3）休止角对比。在煤散料下落完成后，使用MATLAB软件对仿真和试验中获得的煤散料堆轮廓进行拟合。如图7a和表4所示，真实试验中的煤散料休止角为
28.5285 °，仿真中的煤散料休止角为30.6340 °，相对误差为7.38%。图7对真实试验与仿真堆积轮廓进行了线性拟合（分别为线A与线B）。从表4可以看出，随着粒子尺寸的增加，休止角并没有明显的变化趋势。从休止角的相对误差可以看出，总体来讲，随着煤颗粒尺寸的增大，相对误差变大。在仿真中，C-2到C-5的总计算时间甚至少于预测值（C-1的1/h4），其原因是模拟中使用的煤颗粒由小的球形颗粒组合而成，为确保每次模拟料斗中的煤散料总体积为0.001 m3，C-2到C-5中生成的粒子数量少于计算出的粒子数量（C-1的1/h3），使得模拟时间减少，计算速度比预期更快。3. 煤散料回转试验
3.1 煤散料回转试验装置及仿真模型
煤散料回转试验装置由多个零部件组成，包括机架、料槽、上试样和夹具等（图8a）。机架固定在地面上，料槽由电机驱动旋转。上试样及其附属夹具固定在机架上。上试样的下端倾角为75°，为消除外部作用力的潜在干扰，试验前上试样被抬高1 mm。三向压力传感器（T501）通过夹具连接到机架上，传感器组成如图8b所示。为创建回转试验装置简化模型，使用UG软件进行建模（图8c）。为模拟试验过程中煤散料运动行为，将运动子系统添加到动力学软件RecurDyn中，保持料槽旋转速度恒定（4 rad/s），然后将动力学模型导入EDEM软件中作为Wall文件，实现颗粒和刚体运动状态的实时共享。模拟中使用的煤散料体积为0.001 m3，煤散料参数与表2中一致。总模拟时长设定为2 s，其中前0.5 s为煤散料生成时间，后1.5 s为模型运动时间。
在旋转过程中，煤散料随着料槽的运动与上试样碰撞，在上试样的内侧和外侧形成堆积。模型坐标系如图9a所示，X轴与上试样的径向对齐，Y轴平行于槽的内侧线。为测量煤散料的横截面曲线，料槽初始位置定为0°，将90°、180°和270°这3个位置（图9b）的轮廓坐标取平均值，利用MATLAB进行图像处理和曲线拟合，得到相应的轮廓曲线。在实际测量中，使用直尺测量轮廓位置坐标，每隔10 mm选择一个点，共9个点用以拟合轮廓曲线。
3.2 仿真与试验结果对比
1）上试样受力。图10所示为真实试验及仿真中上试样受力的结果对比。图10a为上试样在4 mm煤散料作用下的平均力，将曲线积分并除以区间长度，得到相应的力平均值。试验中施加于上试样的平均力为0.560 N，仿真结果为0.596 N，相对误差为6.43%，仿真与试验结果非常接近。如图10所示，上试样的受力波动随煤颗粒尺寸增加而增加（尤其对于10 mm和12 mm的颗粒），这是由于煤颗粒尺寸的增大导致煤散料的离散性变大，导致作用力波动程度变大。仿真中，随着煤颗粒尺寸的增加，作用于上试样的平均力相对误差逐渐增大，具体数据见表5，在误差允许范围内，粗粒化方法可以进行相应过程的模拟。
表 5 煤散料旋转试验的仿真条件与结果Table 5. Calculation conditions and results of the rotation test仿真条件 煤颗粒尺寸/mm 缩放系数 颗粒数 计算时长 /h 平均力/N 平均力的相对误差/% 轮廓曲线的相关系数 真实试验 4 — — — 0.560 — — C-1 4 1 49183 102.45 0.596 6.43 0.9898 C-2 6 1.5 13676 14.00 0.625 11.61 0.9920 C-3 8 2 5548 4.58 0.633 13.04 0.9841 C-4 10 2.5 2687 1.33 0.707 26.25 0.9761 C-5 12 3 1509 0.40 0.780 39.29 0.9574 2）煤散料剖面曲线。图11a为4 mm煤散料的试验及仿真剖面堆积轮廓的对比图，轮廓位置通过将90°、180°和270°位置处的数值进行平均得到。水平和垂直坐标是相对于原点的位移（基于图9a中的坐标系）。通过Matlab对煤散料剖面曲线进行拟合，试验和仿真数据的相关系数R分别为0.950和0.936 3，比较试验和仿真剖面曲线，得到其相关系数为
0.9898 ，表明两者之间的煤颗粒堆积轮廓非常相似。图11为煤颗粒尺寸从4 mm到12 mm的煤散料剖面轮廓曲线的比较。总体来讲，随着煤颗粒尺寸逐渐增大，相关系数变小，这可以归因于煤散料的流动性随尺寸增大而减小，相关系数数值见表5。另外，C-2到C-5的总模拟时间均少于预测值（C-1的1/h4），与煤散料下落试验类似。4. 结 论
1）利用量纲分析，获得了颗粒系统每个物理量在原始系统和缩尺系统之间的尺度关系，这为离散元接触模型中接触参数的处理提供了理论基础。采用多尺度描述方法，在粗粒化和原始系统之间建立了RVE单元。通过应用不同系统的RVE单元之间的质量、动量和能量守恒关系，获得了两个不同尺度（宏观和细观）上的关系。从精确缩尺模型得到的缩放定律在细观尺度上完全适用于粗粒化系统。
2）通过对煤散料下落和回转试验中的冲击力、煤散料休止角、上试样受力和煤堆剖面曲线进行离散元计算，验证了该方法的有效性。在可接受的误差范围内（依据具体问题确定误差范围），粗粒化方法可有效模拟相应的过程，为工程尺度上的大规模离散元计算提供了高效的解决方案。
3）文中研究对象为干煤散料，对于含水煤散料还需要进一步深入研究。尽管文章提出了双尺度系统，但目前还没有理论方法来预测粗粒化和原始系统之间的计算误差，未来需要研究解决这个问题。此外，还需要进一步研究与旋转相关的物理量尺度规律，这在文中没有进行详尽探究。
-
![]()
图 2 原始系统及粗粒化系统中的代表性单元
Figure 2. Representative units in original system and coarse-grained system
表 1 精确缩尺模型中主要物理量的缩放系数
Table 1 Scaling factors of main physical quantities in exact scaled model
物理量 符号 量纲 缩放系数 长度 L [L] h 时间 t [T] h 密度 ρ [ρ] 1 体积 V [L]3 h3 质量 M [L]3[ ρ] h3 力 F [L]4[T]−2[ρ] h2 应力 σ [L]2[T]−2[ ρ] 1 应变 ε [1] 1 弹性模量 E [L]2[T]−2[ ρ] 1 动能密度 ek [L]2[T]−2[ ρ] 1 力矩 MO [L]5[T]−2[ ρ] h3 
下载: 导出CSV
表 2 仿真计算条件
Table 2 Calculation conditions
参数 C-1 C-2 C-3 C-4 C-5 仿真煤颗粒尺寸/mm 4 6 8 10 12 缩放系数 1 1.5 2 2.5 3
下载: 导出CSV
表 3 验证试验所需仿真参数
Table 3 Simulation input parameters for verification test
材料本征参数 数值 接触参数 数值 煤剪切模量 /Pa 4.7×108 煤−煤恢复因数 6.4×10−1 煤泊松比 3×10−1 煤−煤静摩擦因数 3.29×10−1 煤密度/（kg·m−3） 1.229×103 煤−煤滚动摩擦因数 3.6×10−2 钢剪切模量 /Pa 8×1010 煤−钢恢复因数 6.5×10−1 钢泊松比 3×10−1 煤−钢静摩擦因数 4.6×10−1 钢密度/（kg·m−3） 75×103 煤−钢滚动摩擦因数 3.2×10−2
下载: 导出CSV
表 4 煤散料下落试验的仿真条件与结果
Table 4 Calculation conditions and results of coal falling test
仿真条件 煤颗粒
尺寸/mm缩放系数 颗粒数 计算时长 /h 平均冲击力 /N 平均冲击力的
相对误差休止角/(°) 休止角的
相对误差/%真实试验 4 — — — 5.855 — 28.5285 — C-1 4 1 49183 91.18 6.066 3.60% 30.6340 7.38 C-2 6 1.5 13676 11.83 5.986 2.24% 29.2924 2.68 C-3 8 2 5548 2.67 5.171 11.68% 32.6396 14.41 C-4 10 2.5 2687 0.98 5.240 10.50% 24.3230 14.74 C-5 12 3 1509 0.40 5.014 14.36% 23.0934 19.05
下载: 导出CSV
表 5 煤散料旋转试验的仿真条件与结果
Table 5 Calculation conditions and results of the rotation test
仿真条件 煤颗粒尺寸/mm 缩放系数 颗粒数 计算时长 /h 平均力/N 平均力的相对误差/% 轮廓曲线的相关系数 真实试验 4 — — — 0.560 — — C-1 4 1 49183 102.45 0.596 6.43 0.9898 C-2 6 1.5 13676 14.00 0.625 11.61 0.9920 C-3 8 2 5548 4.58 0.633 13.04 0.9841 C-4 10 2.5 2687 1.33 0.707 26.25 0.9761 C-5 12 3 1509 0.40 0.780 39.29 0.9574
下载: 导出CSV
-
[1] ROUABAH M,BOURGEOIS S,BRIANON S,et al. A numerical tool to predict powder behaviour for pharmaceutical handling and processing[J]. Journal of Drug Delivery Science and Technology,2022,70:103258. doi: 10.1016/j.jddst.2022.103258
[2] ZENGY,MAOB Q,JIAF G,et al. Modelling of grain breakage of in a vertical rice mill based on DEM simulation combining particle replacement model[J]. Biosystems Engineering,2022,215:32−48. doi: 10.1016/j.biosystemseng.2021.12.022
[3] XU WJ,WANG L,CHENG K. The failure and river blocking mechanism of large-scale anti-dip rock landslide induced by earthquake[J]. Rock Mechanics and Rock Engineering,2022,55(8):4941−4961. doi: 10.1007/s00603-022-02903-x
[4] 张美晨,赵丽娟,李明昊,等. 基于双向耦合法的采煤机螺旋滚筒振动特性分析及实验研究[J]. 煤炭科学技术,2024,52(3):200−216 ZHANG Meichen,ZHAO Lijuan,LIM inghao, et al. Analysis and experimental study on the vibration characteristics of the spiral drum of a shearer based on two-way coupling methods[J]. Journal of China Coal Society,2024,52(3):200−216
[5] ZHAO Y,SUN X X,MENG W J. Research on the axial velocity of the raw coal particles in vertical screw conveyor by using the discrete element method[J]. Journal of Mechanical Science and Technology,2021,35(6):2551−2560. doi: 10.1007/s12206-021-0526-z
[6] WANG Z Q,WANG X J,ZHUANG J J,et al. Multiple parameter collaborative optimization of a particle separation equipment for coal cleaning production[J]. Journal of Environmental Chemical Engineering,2021,9(4):105646. doi: 10.1016/j.jece.2021.105646
[7] GAN J Q,ZHOU Z Y,YU A B. A GPU-based DEM approach for modelling of particulate systems[J]. Powder Technology,2016,301:1172−1182. doi: 10.1016/j.powtec.2016.07.072
[8] KUWAGI K,MOKHTAR M A,OKADA H, et al. Numerical experiment of thermoset particles in surface modification system with discrete element method (Quantization of cohesive force between particles by agglomerates analysis)[J]. Numerical Heat Transfer,Part A:Applications,2009,56(8):647−664.
[9] MOKHTAR M A,KUWAGI K,TAKAMI T,et al. Validation of the similar particle assembly (SPA) model for the fluidization of Geldart's group A and D particles[J]. AIChE Journal,2012,58(1):87−98. doi: 10.1002/aic.12568
[10] SAKAI M,ABE M,SHIGETO Y,et al. Verification and validation of a coarse grain model of the DEM in a bubbling fluidized bed[J]. Chemical Engineering Journal,2014,244:33−43. doi: 10.1016/j.cej.2014.01.029
[11] SAKAIM,KOSHIZUKA S. Large-scale discrete element modeling in pneumatic conveying[J]. Chemical Engineering Science,2009,64(3):533−539. doi: 10.1016/j.ces.2008.10.003
[12] TAKABATAKEK,MORIY,KHINASTJG,et al. Numerical investigation of a coarse-grain discrete element method in solid mixing in a spouted bed[J]. Chemical Engineering Journal,2018,346:416−426. doi: 10.1016/j.cej.2018.04.015
[13] BIERWISCHC,KRAFTT,RIEDEL H,et al. Three-dimensional discrete element models for the granular statics and dynamics of powders in cavity filling[J]. Journal of the Mechanics and Physics of Solids,2009,57(1):10−31. doi: 10.1016/j.jmps.2008.10.006
[14] NAKAMURAH,TAKIMOTO H,KISHIDA N,et al. Coarse-grained discrete element method for granular shear flow[J]. Chemical Engineering Journal Advances,2020,4:100050. doi: 10.1016/j.ceja.2020.100050
[15] KISHIDAN,NAKAMURA H,TAKIMOTOH,et al. Coarse-grained discrete element simulation of particle flow and mixing in a vertical high-shear mixer[J]. Powder Technology,2021,390:1−10. doi: 10.1016/j.powtec.2021.05.028
[16] LU L Q,XU J,GE W,et al. EMMS-based discrete particle method (EMMS–DPM) for simulation of gas–solid flows[J]. Chemical Engineering Science,2014,120:67−87. doi: 10.1016/j.ces.2014.08.004
[17] CHU K W,CHEN J,YU A B. Applicability of a coarse-grained CFD–DEM model on dense medium cyclone[J]. Minerals Engineering,2016,90:43−54. doi: 10.1016/j.mineng.2016.01.020
[18] TAUSENDSCHÖNJ,KOLEHMAINENJ,SUNDARESAN S,et al. Coarse graining Euler-Lagrange simulations of cohesive particle fluidization[J]. Powder Technology,2020,364:167−182. doi: 10.1016/j.powtec.2020.01.056
[19] 程宏旸,THOMAS Weinhart. 关于采用粗粒化提高颗粒材料多尺度模拟守恒特性的研究[J]. 计算力学学报,2022,39(3):373−380. doi: 10.7511/jslxCMGM202214 CHENG Hongyang,THOMAS Weinhart. On the conservation properties of CG-enriched concurrent coupling methods for muti-scale modeling of granular materials[J]. Chinese Journal of Computational Mechanics,2022,39(3):373−380. doi: 10.7511/jslxCMGM202214
[20] FENG Y T,HAN K,Owen D R J,et al. On upscaling of discrete element models:similarity principles[J]. Engineering Computations,2009,26:599−609. doi: 10.1108/02644400910975405
[21] FENG Y T,OWEN D R J. Discrete element modelling of large scale particle systems—I:exact scaling laws[J]. Computational Particle Mechanics,2014,1:159−168. doi: 10.1007/s40571-014-0010-y
[22] 赵婷婷,冯云田. 大规模颗粒系统的精确缩尺和粗粒化离散元方法[J]. 计算力学学报,2022,39:365−372. ZHAO Tingting,FENG Yuntian. Exact scaling laws and coarse-grained discrete element modellingof largescale granular systems[J]. Chinese Journal of Computational Mechanics,2022,39:365−372.
[23] XIA R,LI B,WANG X W,et al. Measurement and calibration of the discrete element parameters of wet bulk coal[J]. Measurement,2019,142:84−95. doi: 10.1016/j.measurement.2019.04.069
[24] WANG Z S,LI B,LIANG C,et al. Response analysis of a scraper conveyor under chain faults based on MBD-DEM-FEM[J]. Strojniski Vestnik-Journal of Mechanical Engineering,2021,67(10):501−516. doi: 10.5545/sv-jme.2021.7300
[25] 李铁军,王学文,李 博,等. 基于离散元法的煤颗粒模型参数优化[J]. 中国粉体技术,2018,24(5):6−12. LI Tiejun,WANG Xuewen,LI Bo, et al. Optimization method for coal particle model parameters based on discrete element method [J]. China Powder Science and Technology,2018,24(5):6−12.
计量
- 文章访问数: 87
- HTML全文浏览量: 10
- PDF下载量: 27
