Model and numerical solution method the failure of air volume and air flow components in the ventilation network during the period of external fire
-
摘要:
为了尽可能复现火灾时期通风网络内风量波动、风流逆转、风温异常、有毒有害气体蔓延、氧气体积分数降低等灾变过程，构建了外因火灾通风网络风质风量失效物理模型，采用时序化分析方法进行了解耦处理。通过风阻、火风压、风流密度关键通风参量计算模型，将基于质量流量守恒的通风网络回路风压平衡方程、通风网络风流−围岩热量交换−传播控制方程、有毒有害气体弥散−扩散控制方程、氧气消耗−运移控制方程进行联立，构建了外因火灾通风网络“风量风向−风温−有毒有害气体体积分数−氧气体积分数”非稳定多场耦合数学模型。采用回路风量牛顿法求解通风网络风量风向变化，采用迎风离散格式的有限差分法求解通风网络内风流传热−传质过程，建立了火灾时期全风网“节点−分支−网络”风量风质失效多物理场耦合解算方法，构建了多物理场间接耦合解算流程，编写了解算程序。通过在角联结构试验巷道内开展外因火灾试验，初步验证了模型的准确性与适用性，在此基础上分别模拟了进风区分支与用风区分支火源位置情况下复杂通风网络内风量变化、风流逆转、风温异常、有毒有害气体蔓延、氧气体积分数波动的动态变化过程。模拟结果表明：火源位置在进风井条件下风量显著变化分支达到通风网络分支总数的61%，风流逆转分支数量占通风网络分支总数的11%，风温显著变化分支数量达到通风网络分支总数的76%，CO侵入分支数量达到通风网络分支总数的84%，氧气体积分数显著变化分支数量达到通风网络分支总数的41%，进风井分支发生外因火灾造成的灾害影响范围和影响强度均明显大于其他分支。
Abstract:In order to eeproduced as much as possible the disaster processes such as air volume fluctuations, airflow reversal, abnormal airflow temperature, the spread of toxic and harmful smoke, and the decrease in oxygen volume fraction in the ventilation network during the period of fire occurrence, the airflow pressure balance equation of the ventilation network circuit based on the conservation of airflow mass flow rate, the control equations for heat exchange and propagation in ventilation networks, the control equations for the diffusion and dispersion of toxic and harmful gases in ventilation networks, and the control equation of the consumption and transport of oxygen component in ventilation network were constructed. Based on the theory of multi physics field coupling, the above equations were combined using key ventilation parameters such as airflow resistance, fire pressure, and air flow density to construct the multi field coupling mathematical model for an unstable ventilation network during the period of external fire, which included “air volume and airflow direction, airflow temperature, toxic and harmful gases volume fraction, and oxygen volume fraction”. Based on the Newton method of loop air volume, the variation of air volume and airflow direction in the ventilation network was solved. The finite difference method of upwind discrete format was used to solve the heat and mass transfer process of airflow in the ventilation network. The indirect coupled multi physical field solution method based on temporal logic was adopted to construct the solution process. The accuracy and applicability of the model were preliminarily verified by carrying out external fire experiments in the test roadway of the angular structure, and on this basis, the dynamic change process of air volume change, wind flow reversal, abnormal wind temperature, spread of toxic and harmful gases, and oxygen volume fraction fluctuation in the complex ventilation network under the location of the fire source of the branch in the air inlet area and the branch branch in the air consumption area were simulated respectively.The simulation results showed that under the condition that the ignition source was located in the air inlet well, the number of branches with significant changes in air volume reached 61% of the total number of branches in the ventilation network, the number of branches with significant changes in wind flow reached 11% of the total number of branches in the ventilation network, the number of branches with significant changes in wind temperature reached 76% of the total number of branches in the ventilation network, the number of branches with CO intrusion reached 84% of the total number of branches in the ventilation network, and the number of branches with significant changes in oxygen volume fraction reached 41% of the total number of branches in the ventilation network. The impact scope and intensity of disasters caused by fire in the air inlet shaft branch were significantly greater than those in other branches.
-
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）文中研究对象为干煤散料，对于含水煤散料还需要进一步深入研究。尽管文章提出了双尺度系统，但目前还没有理论方法来预测粗粒化和原始系统之间的计算误差，未来需要研究解决这个问题。此外，还需要进一步研究与旋转相关的物理量尺度规律，这在文中没有进行详尽探究。
-
![]()
图 1 外因火灾通风网络风量风质失效耦合作用机制
Figure 1. Coupling mechanism of failure of air volume and air flow components in ventilation network during the period of external fires
![]()
图 2 外因火灾通风网络风量风质失效耦合流程
Figure 2. Coupling process of failure of air volume and air flow components in ventilation network during the period of external fires
![]()
图 3 外因火灾通风网络风质风量失效解算方法
Figure 3. Calculation method for sub process of failure of air volume and air flow components in ventilation network during the period of external fires
![]()
图 4 通风网络分支元素的风流传热解算网格划分示意
Figure 4. Schematic meshing of the wind transfer heat transfer of the branch elements of the ventilation network
![]()
图 5 通风网络节点元素的风流传热解算网格划分示意
Figure 5. Schematic diagram of the calculation of the air heat transfer of the node elements of the ventilation network
![]()
图 6 巷道围岩内部节点控制体划分示意
Figure 6. Schematic diagram of the division of the internal node control body in the surrounding rock of the roadway
![]()
图 7 通风网络分支元素的风流传质解算网格划分示意
Figure 7. Schematic diagram of the solving meshing of the airflow quality of the branch elements of the ventilation network
![]()
图 8 通风网络节点元素的风流传质解算网格划分示意
Figure 8. Schematic diagram of the solving meshing of the airflow quality of the node elements of the ventilation network
![]()
图 9 火灾时期通风网络内风质风量失效耦合解算流程
Figure 9. Coupled calculation process for failure of air volume and air flow components in ventilation network during theperiod of fire occurrence
![]()
图 10 角联结构试验巷道通风系统图与通风网络图
Figure 10. Ventilation system diagram and ventilation network diagram of corner joint structure test tunnel
![]()
图 11 角联结构试验巷道各分支风量变化情况
Figure 11. Changes in air volume and airflow direction of branches of the corner joint structure test tunnel
![]()
图 12 角联结构试验巷道各节点风流温度变化情况
Figure 12. Changes in temperature of airflow at each node of the corner joint structure test tunnel
![]()
图 13 角联结构试验巷道各节点风流CO体积分数变化情况
Figure 13. Changes in CO volume fraction of airflow at each node of the corner joint structure test tunnel
![]()
图 14 角联结构试验巷道各节点风流O2体积分数变化情况
Figure 14. Changes in O2 volume fraction of airflow at each node of the corner joint structure test tunnel
![]()
图 16 火灾发生30 min之后通风网络内分支风量风向变化模拟结果
注：粉色分支代表分支风量增大；蓝色分支代表分支风量减小；红色分支代表分支风流逆转。
Figure 16. Simulation results of changes in air volume and airflow direction of branches in the ventilation network after 30 minutes of fire
![]()
图 17 火灾发生30 min之后通风网络分支风温变化模拟结果
注：红色分支代表风温变化分支
Figure 17. Simulation results of airflow temperature changes of branches in ventilation network after 30 minutes of fire
![]()
图 18 外因火灾发生30 min后通风网络内分支CO蔓延模拟结果
注：红色分支代表CO侵入分支
Figure 18. Simulation results of CO spread of branches in the ventilation network after 30 minutes of fire
![]()
图 19 外因火灾发生30 min后通风网络内分支氧气体积分数变化模拟结果
注：红色分支代表氧浓度降低分支
Figure 19. Simulation results of oxygen volume fraction changes of branches in the ventilation network after 30 minutes of fire
表 1 外因火灾模拟参数取值
Table 1 Parameter values for external fire simulation
序号 参数类型 参数名称 参数取值 1 物性参数 空气导热系数/(W·m−1·K−1) 2.31 2 巷道围岩比热容/(J·kg−1·K−1) 878 3 围岩与空气对流换热系数/(W·m−2·K−1) 14.75 4 氧气自由扩散系数/(m2·s−1) 1.78×10−5 5 一氧化碳自由扩散系数/(m2·s−1) 1.43×10−5 6 风机设备 风机特性参数A0 500 7 风机特性参数A1 0.1 8 风机特性参数A2 −0.7 9 巷道气候和巷道岩温条件参数 矿井进风流温度/℃ 19.5 10 进风流氧气体积分数/% 20 11 巷道壁面温度/℃ 19.0 巷道原始岩温/℃ 45 12 巷道围岩散热圈半径/m 32 14 火源特征参数 火源生成热量速率/(kJ·s−1) 5.45 15 火源生成CO速率/(m3·s−1) 0.34 16 火源消耗氧气速率/(m3·s−1) 0.86 17 数值解算设置参数 模拟控制单元尺寸/m 0.5 18 风网风量解算残差/% 0.1 19 非稳态模拟时间步长/s 2 20 模拟时长/s 1800 
下载: 导出CSV
-
[1] 王德明,邵振鲁,朱云飞. 煤矿热动力重大灾害中的几个科学问题[J]. 煤炭学报,2021,46(1):57−64. WANG Deming,SHAO Zhenlu,ZHU Yunfei. Several scientific issues on major thermodynamic disasters in coal mines[J]. Journal of China Coal Society,2021,46(1):57−64.
[2] 王德明,张伟,王和堂,等. 煤矿热动力重大灾害的不确定性风险特性研究[J]. 采矿与安全工程学报,2023,40(4):826−837. WANG Deming,ZHANG Wei,WANG Hetang,et al. Uncertain risk characteristics of major thermodynamic disasters in underground coal mines[J]. Journal of Mining & Safety Engineering,2023,40(4):826−837.
[3] 周福宝,辛海会,魏连江,等. 矿井智能通风理论与技术研究进展[J]. 煤炭科学技术,2023,51(1):313−328. ZHOU Fubao,XIN Haihui,WEI Lianjiang,et al. Research progress of mine intelligent ventilation theory and technology[J]. Coal Science and Technology,2023,51(1):313−328.
[4] 张浪,刘彦青. 矿井智能通风与关键技术研究[J]. 煤炭科学技术,2024,52(1):178−195. doi: 10.12438/cst.2023-1987 ZHANG Lang,LIU Yanqing. Research on technology of key steps of intelligent ventilation in mines[J]. Coal Science and Technology,2024,52(1):178−195. doi: 10.12438/cst.2023-1987
[5] 李祥春,蒋颖,李梅生. 巷道火灾时期流场及瓦斯浓度变化规律数值模拟研究[J]. 煤炭科学技术,2019,47(5):119−125. LI Xiangchun,JIANG Ying,LI Meisheng. Study on numerical simulation of variations of airflow field and gas concentration during roadway fire[J]. Coal Science and Technology,2019,47(5):119−125.
[6] 张春华,康璇,申嘉辉. 矿井L形巷道火灾蔓延规律模拟研究[J]. 安全与环境学报,2022,22(6):3111−3118. ZHANG Chunhua,KANG Xuan,SHEN Jiahui. Simulation study on fire spread law of L-shaped roadway in mine[J]. Journal of Safety and Environment,2022,22(6):3111−3118.
[7] 田水承,窦培谦,张成镇. 基于Pyrosim的风速对矿井火灾蔓延规律影响研究[J]. 金属矿山,2020(2):199−204. TIAN Shuicheng,DOU Peiqian,ZHANG Chengzhen. Impact of different wind speeds on the law of mine fire spread based on pyrosim[J]. Metal Mine,2020(2):199−204.
[8] 刘雨晴,张培红. 大坡度倾斜巷道火灾烟气温度分布特征研究[J]. 中国安全科学学报,2021,31(4):156−162. LIU Yuqing,ZHANG Peihong. Study on fire-induced smoke temperature distribution characteristics in largely inclined roadway[J]. China Safety Science Journal,2021,31(4):156−162.
[9] 李宗翔,王海文,李腾,等. 下行风流火灾管道试验与烟流动力特征研究[J]. 安全与环境学报,2023,23(2):391−396. LI Zongxiang,WANG Haiwen,LI Teng,et al. Downstream airflow fire dust test and smoke flow dynamic characteristics study[J]. Journal of Safety and Environment,2023,23(2):391−396.
[10] 李晴,康建宏,周福宝,等. 全尺寸巷/隧道火灾风烟流温度预测模型与验证[J]. 中国安全生产科学技术,2022,18(8):5−12. LI Qing,KANG Jianhong,ZHOU Fubao,et al. Prediction model and verification of smoke flow temperature in full-scale roadway/tunnel fires[J]. Journal of Safety Science and Technology,2022,18(8):5−12.
[11] 马砺,张鹏宇,李超华,等. 不同封堵条件下巷道火灾烟气温度场预测模型[J]. 西安科技大学学报,2021,41(2):213−220. MA Li,ZHANG Pengyu,LI Chaohua,et al. Prediction model of roadway fire temperature field under different plugging conditions[J]. Journal of Xi’an University of Science and Technology,2021,41(2):213−220.
[12] 王海燕,徐祚卉,徐良伟,等. 矿井火灾时期火源分支巷道风阻变化趋势试验研究[J]. 安全与环境学报,2023,23(8):2661−2668. WANG Haiyan,XU Zuohui,XU Liangwei,et al. Experimental study on the variation law of wind resistance in branch roadway of fire source during mine fire[J]. Journal of Safety and Environment,2023,23(8):2661−2668.
[13] 白志鹏,王南,杜珮颖,等. 矿井复杂管网内有害气体的非稳态运移模型与试验研究[J]. 煤炭工程,2021,53(11):135−139. BAI Zhipeng,WANG Nan,DU Peiying,et al. Unsteady migration model of harmful gas in complex pipe network of mine and the experiment verification[J]. Coal Engineering,2021,53(11):135−139.
[14] 周福宝,王德明. 矿山通风网络火灾的计算机仿真模拟[J]. 金属矿山,2008(2):115−117,151. doi: 10.3321/j.issn:1001-1250.2008.02.029 ZHOU Fubao,WANG Deming. Computer simulation of mine fire of ventilation network[J]. Metal Mine,2008(2):115−117,151. doi: 10.3321/j.issn:1001-1250.2008.02.029
[15] 李贺,田丽,曾钢,等. 基于FDS的风速对矿井火灾蔓延规律的影响研究[J]. 中国安全生产科学技术,2022,18(5):143−149. LI He,TIAN Li,ZENG Gang,et al. Study on influence of wind speed on spread law of mine fire based on FDS[J]. Journal of Safety Science and Technology,2022,18(5):143−149.
[16] 戚宜欣,王省身,鲍庆国. 矿井火灾时期风流流动及通风系统变化的动态模拟[J]. 中国矿业大学学报,1995,24(3):19−23. QI Yixin,WANG Xingshen,BAO Qingguo. Dynamic simulation of the flow of mine air and the change of ventilation systems during mine fire[J]. Journal of China University of Mining & Technology,1995,24(3):19−23.
[17] 李宗翔,王雅迪,李林. 上行风流火灾3D矿井通风系统灾变过程仿真[J]. 煤炭学报,2015,40(1):115−121. LI Zongxiang,WANG Yadi,LI Lin. 3D simulation of disaster process in mine ventilation system during fire period[J]. Journal of China Coal Society,2015,40(1):115−121.
[18] 李宗翔,张慧博,路宝生,等. 矿井系统下行风流火灾试验与TF1M3D平台仿真研究[J]. 中国安全生产科学技术,2018,14(1):30−34. LI Zongxiang,ZHANG Huibo,LU Baosheng,et al. Research of test and TF1M3D platform simulation on descending airflow fire in mine[J]. Journal of Safety Science and Technology,2018,14(1):30−34.
[19] 张明乾,李宗翔,李腾,等. 矿井网域系统上行火灾灾情演变特征研究[J]. 煤炭学报,2021,46(S2):785−792. ZHANG Mingqian,LI Zongxiang,LI Teng,et al. Experiments of mine fire in ascensional airflow and the mine fire control simulation[J]. Journal of China Coal Society,2021,46(S2):785−792.
[20] 张景钢,孙春峰,张海洋,等. 矿井火灾模拟解算软件开发研究[J]. 华北科技学院学报,2015,12(1):30−35. ZHANG Jinggang,SUN Chunfeng,ZHANG Haiyang,et al. Development & research of simulation calculating software of mine fire simulation[J]. Journal of North China Institute of Science and Technology,2015,12(1):30−35.
[21] 裴晓东,郝海清,王凯,等. 矿井复杂风网火灾风烟流应急调控技术及应用[J]. 煤炭科学技术,2023,51(5):124−132. PEI Xiaodong,HAO Haiqing,WANG Kai,et al. Research and application of fire air and smoke flow emergency control technology for mine complex ventilation network[J]. Coal Science and Technology,2023,51(5):124−132.
[22] 郝海清,王凯,张春玉,等. 矿井皮带巷火灾风烟流场-区-网演化与调控规律[J]. 中国矿业大学学报,2021,50(4):716−724. HAO Haiqing,WANG Kai,ZHANG Chunyu,et al. Evolution and regulation law of wind and smoke flow field area network in mine belt roadway fire[J]. Journal of China University of Mining & Technology,2021,50(4):716−724.
[23] 李翠平,曹志国,钟媛. 矿井火灾的场量模型构建及其可视化仿真[J]. 煤炭学报,2015,40(4):902−908. LI Cuiping,CAO Zhiguo,ZHONG Yuan. Field variables modeling and visualization simulation of fire disaster in underground mine[J]. Journal of China Coal Society,2015,40(4):902−908.
[24] 傅培舫. 实际巷道火灾过程热物理参数变化规律与计算机仿真的研究[D]. 徐州:中国矿业大学,2002. FU Peifang. Study on variation law of thermophysical parameters and computer simulation in actual roadway fire process[D]. Xuzhou:China University of Mining and Technology,2002.
[25] 高建良,杨明. 巷道围岩温度分布及调热圈半径的影响因素分析[J]. 中国安全科学学报,2005,15(2):73−76. doi: 10.3969/j.issn.1003-3033.2005.02.017 GAO Jianliang,YANG Ming. Analysis of the factors influencing temperature distribution of surrounding rock and cooled zone radius[J]. China Safety Science Journal,2005,15(2):73−76. doi: 10.3969/j.issn.1003-3033.2005.02.017
[26] 周心权,吴兵. 矿井火灾救灾理论与实践[M]. 北京:煤炭工业出版社,1996. [27] 姚勇征,张文明,吴兵,等. 巷道火灾对通风系统影响的全尺寸试验与模拟[J]. 中国矿业大学学报,2021,50(4):709−715. YAO Yongzheng,ZHANG Wenming,WU Bing,et al. Full-scale experimental and simulation study of the influences of laneway fire on ventilation system[J]. Journal of China University of Mining & Technology,2021,50(4):709−715.
-
期刊类型引用(1)
1. 徐乔木，赵星杰，卜侃侃. 旋转多斗装车系统受限空间内复杂粒度散料堆积特性研究. 煤. 2025(04): 86-91 . 
百度学术
其他类型引用(1)
计量
- 文章访问数: 29
- HTML全文浏览量: 2
- PDF下载量: 12
- 被引次数: 2
