Simulation study of top coal caving and conveying process based on smoothed particle hydrodynamics
-
摘要: 目前针对综放开采中顶煤放出规律的数值模拟研究中，对于顶煤运动的连续−非连续性问题需复杂的耦合算法，必须解决煤岩界面信息精确交互问题，且忽略了刮板输送机输送过程。针对该问题，基于光滑粒子动力学构建了无网格数值计算模型，通过建立连续介质力学控制方程的光滑粒子动力学离散方程，并引入弹塑性土体本构模型和Drucker−Prager屈服准则，实现了顶煤坍塌、运移、放出过程的动态模拟。考虑采场实际放煤和输煤过程，构建了刮板输送机模型，模拟沿工作面水平方向顶煤放出和底煤输送过程，得到不同刮板输送机运行速度（0~1.5 m/s)下的煤岩界面和煤流速度变化规律。仿真结果表明：弹塑性土体本构模型可有效模拟颗粒的流动行为，通过设定摩擦角、弹性模量等材料参数，避免了传统离散元法模型的参数不定问题；煤流速度稳定后，放煤口附近的顶煤应力分布呈 “双峰”形态；刮板输送机运行速度对放煤时间影响较大，但对终止的煤岩界面和放出体形状影响较小；多支架同时放煤需考虑刮板输送机的输送能力，不同支架之间的底煤输送干涉可能导致放煤口的堵塞效应； “见矸关门”准则导致不同放煤口放煤量存在差异，40个放煤口顶煤放出量的标准差（7.52 m2）高于自动放煤的标准差（1.93 m2）。Abstract: Currently, in the numerical simulation research on the release laws of top coal during fully mechanized mining, complex coupling algorithms are required to address the continuity-discontinuity issues of top coal movement and ensure precise interaction of coal-rock interface information. However, the conveying process of scraper conveyors is typically neglected in these simulations. To address this problem, a meshless numerical computation model was constructed based on smoothed particle hydrodynamics (SPH). The discrete equations of SPH, derived from the control equations of continuous medium mechanics, were established. An elastic-plastic soil constitutive model along with the Drucker-Prager yield criterion were introduced to achieve dynamic simulation of the caving, movement, and release processes of the top coal. Considering the actual coal release and conveying processes in the mining area, a scraper conveyor model was constructed to simulate the release of top coal and the conveying of bottom coal along the working face, obtaining the variations in coal-rock interface and coal flow velocity at different scraper conveyor operating speeds (0-1.5 m/s). The simulation results indicated that the elastic-plastic soil constitutive model effectively simulated the flow behavior of particles. By setting parameters such as friction angle and elastic modulus, the issue of uncertain parameters commonly found in traditional discrete element models was avoided. After stabilization of the coal flow velocity, the stress distribution of the top coal near the coal drawing outlet exhibited a "double peak" pattern. The operating speed of the scraper conveyor significantly impacted the coal drawing time, while its effect on the coal-rock interface at termination and the shape of the released body was minimal. When multiple supports released coal simultaneously, the conveying capacity of the scraper conveyor needed to be considered, as interference in bottom coal transportation between different supports could lead to blockage effects at the release port. The "gangue closing" rule resulted in variations in the amount of coal drawing at different coal drawing outlets, with the standard deviation of top coal drawing amount from 40 coal drawing outlets (7.52 m²) being greater than that of automatic coal drawing (1.93 m²).
-
0. 引言
我国厚煤层储量和产量占地下煤炭资源和产量的45%，放顶煤技术是开采厚煤层的主要方法之一[1-2]。传统放顶煤的终止条件是“见矸关门”，需人工控制放煤口。目前自动化顶煤冒落控制系统已应用[3]。由于放顶煤工艺和条件的复杂性，全面实现智能化综放还有许多难题需要攻克[4]。以往学者直接将顶煤假设为松散体[5-6]，进行散体冒放运移和成拱机理的研究，对松软易碎或裂隙发育的煤层有较好的适用性[7]，但大多忽略了刮板输送机与支架放煤的协同过程。在实际放煤工艺中，刮板输送机是保障连续放煤和输送的关键环节。因此，有必要考虑支架放煤与刮板输送机运煤的相互影响，进一步建立符合实际采场环境的放顶煤工艺优化方法，为提高顶煤采出率和放煤效率提供指导。
数值模拟是研究顶煤运动和放出规律的重要工具。常用的数值模拟方法包括离散元法（Discrete Element Method，DEM）、有限差分法（Finite Difference Method，FDM）、有限元法（Finite Element Method，FEM）等。DEM的代表性软件有PFC，EDEM等；FDM，FEM的代表性软件有FLAC3D，LS−DYNA等。张宁波等[8]使用PFC2D软件，基于鹤壁煤电股份有限公司第十煤矿的实际条件，模拟不同顶煤厚度条件下的煤矸流场，分析了煤层厚度及采放比对煤矸流场特征和顶煤放出规律的影响。王家臣等[9]利用PFC3D软件模拟了综放开采过程中煤矸三维放出体形态，分析了不同放煤方式对放出体空间形态的影响。孙强等[10]对浅埋双硬特厚煤层放煤规律进行研究，利用FLAC3D软件分析了综放工作面回采过程中的煤岩破坏规律。何欣等[11]以中煤平朔集团有限公司34201工作面为背景，结合现场地质条件，使用PFC2D软件建立了仰斜角分别为0，13，23，33°时的数值模型，分析了不同条件下的顶煤位移场、接触应力场、煤矸分界面、顶煤位移始动点及顶煤放出体形态。张锦旺等[12]基于BBR（Boundary-Body-Ratio）研究体系，采用PFC软件建立数值模型，模拟了工作面倾角为0～50°时的顶煤放出过程，分析了工作面倾角对散体顶煤放出规律的影响。冯宇峰[13]使用LS−DYNA3D软件模拟不同装药方式和炮孔位置的爆破对顶板破碎效果的影响。邹光华等[14]使用PFC软件模拟不同煤层倾角下的煤岩分界线，建立了等效于放出椭球体面积的煤岩分界线抛物线方程。姜志刚等[15]针对开滦（集团）有限责任公司唐山矿业分公司大倾角综放工作面，使用PFC2D软件研究了煤岩分界线的不对称分布特征及放煤工艺优化方法。Wang Shuai等[16]使用PFC2D软件研究了重采工作面散体顶煤和顶板运动及煤岩质界面演化规律。Liu Yang等[17]使用PFC2D软件模拟和比较不同垮落技术下的顶煤采出率和采动平衡过程，认为独立簇群垮落技术可提高厚煤层的顶煤采出率。Wei Weijie等[18]测量了顶煤块体的尺寸分布，通过离散元计算和物理实验，研究了特厚煤层中顶煤和矸石块体的混合运动机制。
在综放开采中，顶煤产生渐进破坏过程，是典型的连续−非连续性问题。而传统的DEM，FEM一般仅适用于单一的连续或非连续问题。研究连续−非连续转换过程需不同方法进行耦合来实现。张文辉等[19]使用FLAC3D和PFC3D软件，通过界面耦合方法实现连续块体和离散颗粒的动态耦合模拟，研究特厚煤层大采高综放开采中顶煤的渐进破坏行为。Zhang Qunlei等[20]采用连续介质与非连续单元耦合方法（Continuum Discrete Element Method，CDEM）来模拟顶煤放出过程，引入液压支架的本构模型模拟液压支架与周围岩石的耦合作用。李东印等[21]采用CDEM建立沿工作面走向的二维数值模型，探究放煤步距对顶煤遗失和采出率的影响。刘闯等[22]利用CDEM建立沿工作面倾向的数值模型，模拟研究不同顶煤厚度（4.0，8.0，12.0 m）下的放煤过程，分析煤矸分界面和顶煤采出率的变化规律。然而，连续−非连续耦合方法涉及复杂的耦合算法，且在转换过程中必须解决界面信息精确交互问题。近年来，连续粒子方法，如光滑粒子动力学（Smoothed Particle Hydrodynamics，SPH）[23]、物质点方法（Material Point Method，MPM）[24]发展迅速。这些方法不仅继承了无网格技术在处理大变形问题时的优势，还具备了连续介质力学模型的确定性特征，为解决连续−非连续问题提供了一种统一的方案。
此外，以往放煤过程模拟研究往往忽略了刮板输送机的输送过程，煤流流向无阻碍的自由空间，这与实际情况不符。本文采用SPH建立综放开采顶煤放出和底煤输送过程的数值模型，结合非黏聚弹塑性土体本构，模拟颗粒散煤的流动行为，研究了煤岩界面的演化、煤流速度及支架放煤−刮板输送之间的相互影响。
1. SPH基础理论
SPH是一种采用“粒子”离散计算域的数值计算方法，通过将连续介质离散化为粒子，并利用核函数在局部积分域内近似物理量，实现对物理现象的模拟。SPH中场函数通过“核积分”表示法来近似，通过计算相邻粒子的场值并进行加权求和，完成对场函数及其导数的离散化。构建SPH离散方程主要包括核近似和粒子近似2个步骤。后者通过计算核函数影响域内相邻粒子物理量的加权求和，来估计当前粒子的变量信息。
基于核近似原理，计算域内任意场函数$ f\left(\boldsymbol{x}\right) $可近似表示为
$$ f\left(\boldsymbol{x}\right)={\int }_{\varOmega }f\left({\boldsymbol{x}}'\right)W\left(\boldsymbol{x}-{\boldsymbol{x}}', h\right){\mathrm{d}}{\boldsymbol{x}}' $$ (1) 式中：x为计算域内任意点的位置；$ \varOmega $为支持域或影响域；$ W\left(\boldsymbol{x}-{\boldsymbol{x}}', h\right) $为核函数；x'为x支持域$ \varOmega $中任意点的位置；h为光滑长度，本文取初始粒子间距的1.2倍。
场函数梯度的核近似表示为
$$ \nabla f\left(\boldsymbol{x}\right)=-{\int }_{\varOmega }f\left({\boldsymbol{x}}'\right)\nabla W\left(\boldsymbol{x}-{\boldsymbol{x}}', h\right){\mathrm{d}}{\boldsymbol{x}}' $$ (2) 式中：$ \nabla f\left(\boldsymbol{x}\right) $ 为$ f\left(\boldsymbol{x}\right) $的梯度；$ \nabla W\left(\boldsymbol{x}-{\boldsymbol{x}}', h\right) $为核函数的梯度。
由式（2）可知，SPH核近似表示为连续积分形式。因计算域离散为粒子，核近似方程可进一步转换为基于粒子分布的离散形式，即“粒子近似”，其利用支持域粒子的叠加求和完成。SPH计算域离散化及粒子支持域如图1所示。在XY坐标系中，r为计算域内任意点x与其支持域内任意点x'间的距离，k为与光滑函数有关的比例系数，kh表示光滑域的半径。
![]()
图 1 SPH计算域离散化及粒子支持域Figure 1. Smoothed particle hydrodynamics (SPH) computational domain discretisation and particle support domains式（1）和式（2）可写为以下离散形式（即粒子求和近似表达式[23]）：
$$ \left\langle{f\left({\boldsymbol{x}}_{i}\right)}\right\rangle=\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}}{{\rho }_{j}}f\left({\boldsymbol{x}}_{j}\right){W}_{ij} $$ (3) $$ \left\langle{\nabla f\left({\boldsymbol{x}}_{i}\right)}\right\rangle=-\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}}{{\rho }_{j}}f\left({\boldsymbol{x}}_{j}\right){{\nabla }_{i}W}_{ij} $$ (4) 式中：i为粒子索引；$ \left\langle{f\left({\boldsymbol{x}}_{i}\right)}\right\rangle $，$ \left\langle{\nabla f\left({\boldsymbol{x}}_{i}\right)}\right\rangle $分别为$ f\left(\boldsymbol{x}\right) $和$ f\left(\boldsymbol{x}\right) $梯度的粒子近似式；j为位于粒子i支持域内的粒子编号；$ f\left({\boldsymbol{x}}_{i}\right) $，$ f\left({\boldsymbol{x}}_{j}\right) $分别为场函数$ f\left(\boldsymbol{x}\right) $在粒子i和粒子j处的值；$ {N}_{i} $为粒子i支持域内粒子总数；$ {m}_{j} $为粒子j的质量；$ {\rho }_{j} $为粒子j的密度；$ {W}_{ij} $为粒子i与粒子j之间的核函数值，$ {W}_{ij}=W\left(\left|{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}\right|,{h}_{ij}\right) $，$ {h}_{ij} $为粒子i和粒子j光滑长度的平均值；$ {{\nabla }_{i}W}_{ij} $为粒子i和粒子j之间核函数的梯度值，$ {{\nabla }_{i}W}_{ij}={\nabla }_{i}W\left(\left|{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}\right|,{h}_{ij}\right)= \dfrac{{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}}{\left|{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}\right|}\dfrac{\partial {W}_{ij}}{\partial \left|{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}\right|} $ 。
构建SPH离散方程时，可利用式（4）表示控制方程中的偏微分项，但需做必要的转换。以矢量函数的散度为例说明转换方法。SPH中散度通常近似为[23]
$$ \left\langle {\nabla \cdot {f}\left({\boldsymbol{x}}_{i}\right)} \right\rangle =\frac{1}{{\rho }_{i}}\left[\sum _{j=1}^{{N}_{i}} {m}_{j}\left(f\left({\boldsymbol{x}}_{j}\right)-f\left({\boldsymbol{x}}_{i}\right)\right) {\nabla }_{i}{W}_{ij}\right] $$ (5) 式中$ \left\langle {\nabla \cdot {f}\left({\boldsymbol{x}}_{i}\right)} \right\rangle $ 为$ f\left(\boldsymbol{x}\right) $散度的粒子近似式。
需要注意的是，核函数是一个随核距离的增加而减小的函数，以确保每个粒子与其在一定光滑域半径内的相邻粒子相互作用。
f(x)梯度的粒子近似式为[23]
$$ \left\langle{\nabla f\left({\boldsymbol{x}}_{i}\right)}\right\rangle={\rho }_{i}\left[\sum _{j=1}^{{N}_{i}}{m}_{j}\left(\frac{f\left({\boldsymbol{x}}_{j}\right)}{{\rho }_{j}^{2}}+\frac{f\left({\boldsymbol{x}}_{i}\right)}{{\rho }_{i}^{2}}\right){{\nabla }_{i}W}_{ij}\right] $$ (6) $$ \left\langle{\nabla f\left({\boldsymbol{x}}_{i}\right)}\right\rangle={\rho }_{i}\left[\sum _{j=1}^{{N}_{i}}{m}_{j}\left(\frac{f\left({\boldsymbol{x}}_{j}\right)+f\left({\boldsymbol{x}}_{i}\right)}{{\rho }_{i}{\rho }_{j}}\right){{\nabla }_{i}W}_{ij}\right] $$ (7) 本文采用式（7）来离散控制方程中的梯度项。需要说明的是，式（6）和式（7）对数值近似精度的影响没有显著差异。根据本文经验，式（6）也是可行的。
本文选用三次样条分段函数（又称光滑函数）作为核函数，它具有核函数所要求的6个基本性质，包括紧致性、归一性、非负性等。三次样条核函数表达式为
$$ {W}_{ij}=\left\{\begin{array}{*{20}{l}}{\alpha }_{{\mathrm{d}}}\left(\dfrac{2}{3}-{q}^{2}+\dfrac{1}{2}{q}^{3}\right) & 0\le q < 1 \\ \dfrac{{\alpha }_{{\mathrm{d}}} }{6}{\left(2-q\right)}^{3}& 1\leqslant q < 2\\ 0& q\geqslant 2\end{array}\right. $$ (8) 式中：$ \alpha\mathrm{_d} $为正则化系数，用于确保核函数的归一性，在三维空间中$ {\alpha }_{{\mathrm{d}}}=\dfrac{3}{2{\text{π}} {h}^{2}} $；$ q $为无量纲距离，$ q=\dfrac{r'}{h} $，$ r' $为距函数原点的距离。
SPH中两两交互的粒子通常以粒子对的形式存在和定义。一对交互的粒子对i和j的核函数值$ {W}_{ij} $可写为$ W\left(\left|{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}\right|,{h}_{ij}\right) $或$ W\left(\left|{\boldsymbol{x}}_{ij}\right|,{h}_{ij}\right) $。
2. 基于弹塑性土体本构的放煤SPH模型
2.1 连续介质控制方程
采用SPH方法模拟颗粒材料，如土壤、砂子或煤炭颗粒。虽然这些材料本质上是离散的，但被模拟为连续介质。在连续介质力学框架下，颗粒材料的控制方程包括连续性方程（式（9））和动量方程（式（10））。为了便于建立SPH离散方程，控制方程写为拉格朗日形式。
$$ \frac{{\mathrm{d}}\rho }{{\mathrm{d}}t}=-\rho \frac{\partial \boldsymbol{v}}{\partial \boldsymbol{x}} $$ (9) $$ \frac{{\mathrm{d}}\boldsymbol{v}}{{\mathrm{d}}t}=-\frac{1}{\rho }\frac{\partial \boldsymbol{\sigma }}{\partial \boldsymbol{x}}+\frac{{\boldsymbol{F}}_{{\mathrm{c}}}}{m}+\boldsymbol{a} $$ (10) 式中：ρ为材料密度；t为时间；$ \boldsymbol{v} $ 为速度矢量；$ \boldsymbol{\sigma } $为总应力张量；$ {\boldsymbol{F}}_{{\mathrm{c}}} $为介质所受的接触力；m为材料质量；$ \boldsymbol{a} $为重力引起的加速度。
颗粒材料的宏观力学特性通过本构模型来描述，以捕捉材料在外力作用下的变形和失效行为。本文将颗粒材料描述为弹性−完全塑性材料，其应力−应变关系基于以下假设：材料在达到屈服点之前以弹性方式变形，一旦超过屈服点，材料将进入塑性变形阶段。需注意的是，式（9）和式（10）组成的控制方程组是不封闭的，本构模型的引入使方程组封闭。为便于实施本构模型，将总应力张量分解为两部分，如下。
$$ \boldsymbol{\sigma }=-p\delta+\boldsymbol{\tau } $$ (11) 式中：$ p $为各向同性的静水压力；$ \delta $为张量运算的克罗内克符号；$ \boldsymbol{\tau } $为剪应力张量。
在传统的SPH模型中，静水压力一般通过状态方程来求解。静水压力可通过主应力之和计算得到。
$$ p=-\frac{1}{3}{\mathrm{tr}}\left(\boldsymbol{\sigma }\right)=-\frac{1}{3}\left({\sigma }^{XX}+{\sigma }^{YY}+{\sigma }^{ZZ}\right) $$ (12) 式中：$ {\mathrm{tr}}\left(\boldsymbol{\sigma }\right) $为总应力张量的迹；$ {\sigma }^{XX} $, $ {\sigma }^{YY} $, $ {\sigma }^{ZZ} $分别为X，Y，Z 3个方向的主应力分量。
采用弹塑性本构模型，构建应力率张量与应变率张量之间的关系，进而可实现材料瞬时应力的求解，结合屈服准则模拟颗粒材料的塑性行为。该部分内容将在2.3节介绍。
2.2 数值离散方程
基于粒子近似方程（式（5）—式（7）），控制方程（式（9）、式（10））右侧的偏微分项可以写成粒子求和近似形式：
$$ {\left(\rho \frac{\partial \boldsymbol{v}}{\partial \boldsymbol{x}}\right)}_{i}={\left(\rho \nabla \cdot \boldsymbol{v}\right)}_{i}={\rho }_{i}\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}}{{\rho }_{j}}\left({\boldsymbol{v}}_{j}-{\boldsymbol{v}}_{i}\right)\frac{\partial {W}_{ij}}{\partial {\boldsymbol{x}}_{i}} $$ (13) $$ {\left(\frac{1}{\rho }\frac{\partial \boldsymbol{\sigma }}{\partial \boldsymbol{x}}\right)}_{i}=\sum _{j=1}^{{N}_{i}}{m}_{j}\left(\frac{{\boldsymbol{\sigma }}_{i}+{\boldsymbol{\sigma }}_{j}}{{\rho }_{i}{\rho }_{j}}\right)\frac{\partial {W}_{ij}}{\partial {\boldsymbol{x}}_{i}} $$ (14) 式中$ \nabla \cdot \boldsymbol{v} $为速度矢量$ \boldsymbol{v} $的散度。
将式（13）、式（14）代入控制方程，替代偏微分项。离散化的控制方程可以表示为
$$ \frac{{\mathrm{d}}{\rho }_{i}}{{\mathrm{d}}t}={\rho }_{i}\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}}{{\rho }_{j}}\left({\boldsymbol{v}}_{i}-{\boldsymbol{v}}_{j}\right)\frac{\partial {W}_{ij}}{\partial {\boldsymbol{x}}_{i}} $$ (15) $$ \frac{{\mathrm{d}}{\boldsymbol{v}}_{i}}{{\mathrm{d}}t}=\sum _{j=1}^{{N}_{i}}{m}_{j}\left(\frac{{\boldsymbol{\sigma }}_{i}+{\boldsymbol{\sigma }}_{j}}{{\rho }_{i}{\rho }_{j}}\right)\frac{\partial {W}_{ij}}{\partial {\boldsymbol{x}}_{i}}+\frac{{\boldsymbol{F}}_{{\mathrm{c}},i}}{{m}_{i}}+{\boldsymbol{a}}_{i} $$ (16) 式中：$ {\boldsymbol{F}}_{{\mathrm{c}},i} $为粒子i所受的接触力；$ {\boldsymbol{a}}_{i} $为粒子i重力引起的加速度。
式（15）和式（16）即离散化后的控制方程。对计算域中的每一个粒子，均得到对应的控制方程离散公式。
在SPH框架内，粒子系统的控制方程遵循质量和动量守恒定律。SPH粒子被赋予固定的质量，并且在计算过程中保持不变，自动满足质量守恒。粒子的密度会变化，根据式（15）可计算密度导数，按照时间积分格式完成更新。在式（16）表达的动量方程中，粒子的受力项总是成对出现的。例如，粒子i和粒子j为一对粒子对，相互作用力为$ {m}_{j}\left(\dfrac{{\boldsymbol{\sigma }}_{i}+{\boldsymbol{\sigma }}_{j}}{{\rho }_{i}{\rho }_{j}}\right) \times \dfrac{\partial {W}_{ij}}{\partial {\boldsymbol{x}}_{i}} $。核梯度具有特性：$ \dfrac{\partial {W}_{ij}}{\partial {\boldsymbol{x}}_{i}}= -\dfrac{\partial {W}_{ij}}{\partial {\boldsymbol{x}}_{j}} $，可得任意交互的2个粒子的相互作用力大小相等、方向相反，确保了在仅受内力的条件下，整个粒子系统的动量守恒。
2.3 弹塑性土体本构模型
SPH模拟属于显式动力学计算，因此需要根据本构模型建立应力率和应变率之间的关系。类似地，材料的总应变率张量可以分解为弹性部分和塑性部分，即
$$ \dot{\boldsymbol{\varepsilon }}={\dot{\boldsymbol{\varepsilon }}}_{{\mathrm{e}}}+{\dot{\boldsymbol{\varepsilon }}}_{{\mathrm{p}}} $$ (17) 式中：$ \dot{\boldsymbol{\varepsilon }} $为总应变率张量；$ {\dot{\boldsymbol{\varepsilon }}}_{{\mathrm{e}}} $为弹性应变率张量； $ {\dot{\boldsymbol{\varepsilon }}}_{{\mathrm{p}}} $为塑性应变率张量。
根据广义胡克定律建立弹性应变率和应力率之间的关系，按照塑性流动法则确定塑性应变率与应力率的关系，即
$$ {\dot{\boldsymbol{\varepsilon }}}_{{\mathrm{e}}}=\frac{\dot{\boldsymbol{\tau }}}{2G}+K{\mathrm{tr}}\left(\dot{\boldsymbol{\sigma }}\right)\delta $$ (18) $$ {\dot{\boldsymbol{\varepsilon }}}_{{\mathrm{p}}}=\dot{\lambda }\frac{\partial g}{\partial \boldsymbol{\sigma }} $$ (19) 式中：$ \dot{\boldsymbol{\tau }} $为剪应力率张量；G 为剪切模量；K为体积模量；$ {\mathrm{tr}}\left(\dot{\boldsymbol{\sigma }}\right) $为总应力率张量$\dot{\boldsymbol{\sigma }} $的迹；$ \dot{\lambda } $为塑性乘子的变化率；$ g $为塑性势函数。
本文采用关联流动法则，即塑性势函数与材料的屈服函数一致。弹塑性土体本构模型需要与屈服准则结合，模拟特定应力条件下颗粒材料的弹性、塑性变形行为。本文采用岩土力学中常用的Drucker−Prager（D−P）屈服准则，它定义了弹性和塑性变形之间的界限。当应力状态达到由D−P准则定义的屈服面时，材料开始发生塑性流动，如图2所示。其中O为中心，$\sigma ' _1 $为最大主应力轴，$\sigma ' _2 $为中间主应力轴，$\sigma ' _3 $为最小主应力轴，αφ和kc为D−P常数，I1和J2分别为第1和第2应力不变量。
D−P屈服准则表达式为
$$ f\left(I_1,J_2\right)=\sqrt{J_2}+\alpha_{\varphi}I_1-k_{{c}}=0 $$ (20) 式中$ f\left({I}_{1},{J}_{2}\right) $为塑性势函数。
$ {\alpha }_{\varphi }$，$ {k}_{{{c}}} $可与摩尔库仑材料参数摩擦角$ \varphi $和黏聚力$ c $建立以下关系：
$$\left\{\begin{aligned} &{\alpha }_{\varphi }=\frac{\mathrm{tan}\;\varphi }{\sqrt{9+12{{\mathrm{tan}}}^{2}\varphi }}\\ &{k}_{c}=\frac{3c}{\sqrt{9+12{{\mathrm{tan}}}^{2}\;\varphi }} \end{aligned}\right.$$ (21) 塑性乘子的变化率为[25]
$$ \dot{\lambda }=\frac{3{\alpha }_{\varphi }K{\mathrm{tr}}\left(\dot{\boldsymbol{\varepsilon }}\right)+\dfrac{G}{\sqrt{{J}_{2}}}\boldsymbol{\tau }\dot{\boldsymbol{\varepsilon }}}{27{\alpha }_{\varphi }K\mathrm{sin}\;\psi +G} $$ (22) 式中：$ {\mathrm{tr}}\left(\dot{\boldsymbol{\varepsilon }}\right) $为总应变率张量$ \dot{\boldsymbol{\varepsilon }} $的迹，代表体积应变率；$ \psi $为剪胀角。
由此可得应力率（即应力张量的时间导数）表达式：
$$ \frac{\mathrm{d}\boldsymbol{\sigma}}{\mathrm{d}t}=2G\dot{\boldsymbol{\varepsilon}}+\boldsymbol{\sigma}\dot{\boldsymbol{\omega}}+\boldsymbol{\sigma}\dot{\boldsymbol{\omega}}^{\mathrm{T}}+K\mathrm{tr}(\dot{\boldsymbol{\varepsilon}})\delta-\dot{\lambda}\left(9K\mathrm{sin}\; \psi\delta+\frac{G}{\sqrt{J_2}}\boldsymbol{\tau}\right) $$ (23) 式中$ \dot{\boldsymbol{\omega }} $为旋转率张量。
式（23）建立了应力率与应变率、当前应力状态之间的关系，通过该式可根据颗粒系统当前的运动参数计算应力随时间的变化率，进而通过数值积分更新应力，获得当前时刻的应力分布。然后根据D−P准则判定是否进入塑性。在塑性加载下，颗粒的应力状态必须始终保持在塑性屈服面上，因此需要回退映射算法来调整剪应力，确保应力不超过屈服面。应力调整过程如图2（b）所示，具体的应力调整原理和步骤可参照文献[26]。
2.4 SPH模型的实施
在显式计算框架下，一旦确定了第n个时间步的粒子速度和应力，即可利用式（15）、式（16）、式（23）计算当前时间步的导数，用于更新粒子的密度、速度和应力。粒子i的密度可通过式（15）导出的密度导数进行更新。动量方程包含了作用于粒子的各种力，包括静水压力和剪应力。粒子i的速度更新基于式（16）。此外，应力的演变依赖于材料的本构模型，它在计算过程中涉及从弹性到塑性的转变及在这种转变时的应力调整。离散化的SPH方程可以通过数值积分格式求解，如龙格−库塔、预测−校正法等。本文采用计算效率更高的蛙跳算法，场变量和粒子位置通过下式更新[23]：
$$ \left\{\begin{aligned} &{\rho }_{n+\frac{1}{2}}={\rho }_{n-\frac{1}{2}}+{\left(\frac{\mathrm{d}\rho }{\mathrm{d}t}\right)}_{n}\Delta t\\ &{\boldsymbol{v}}_{n+\frac{1}{2}}={\boldsymbol{v}}_{n-\frac{1}{2}}+{\left(\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t}\right)}_{n}\Delta t\\ &{\boldsymbol{\sigma }}_{n+1}={\boldsymbol{\sigma }}_{n}+{\left(\frac{\mathrm{d}\boldsymbol{\sigma }}{\mathrm{d}t}\right)}_{n}\Delta t\\ &{\boldsymbol{x}}_{n+1}={\boldsymbol{x}}_{n}+{\boldsymbol{v}}_{n+\frac{1}{2}}\Delta t \end{aligned}\right. . $$ (24) 式中：$ {\rho }_{n+\frac{1}{2}} $，$ {\boldsymbol{v}}_{n+\frac{1}{2}} $，$ {\boldsymbol{\sigma }}_{n+1} $，$ {\boldsymbol{x}}_{n+1} $分别为对应时间步下的密度、速度矢量、应力张量、位置矢量更新值；n+1/2，n−1/2，n，n+1分别为t+Δt/2，t−Δt/2，t，t+$ \Delta t $时刻的时间步；$ \Delta t $为时间步长。
为了在计算过程中保持数值的稳定性，时间步长须满足Courant−Friedrichs−Lewy（CFL）条件：
$$ \Delta t\leqslant C\mathrm{_{cour}}\frac{h_i}{s} $$ (25) 式中：$ {C}_{{\mathrm{cour}}} $为CFL常数，设定为0.1；hi为粒子i的光滑长度；s为材料声速。
本文中声速s取值取决于体积模量K，可通过$ s=\sqrt{{K}/{\rho }} $ 来估计。
根据以上数值流程，编制了基于Fortran90的计算代码，相关SPH模型参数设定见表1，具体的程序结构可参照文献[23]。本文后续的计算均采用该程序完成，通过Visual Studio编译环境进行运算，计算结果后处理采用Tecplot 360软件完成。
表 1 SPH模型参数设定Table 1. Parameters seting in SPH model参数 取值 参数 取值 初始粒子间距
$ {d}_{{\mathrm{ini}}} $/m取决于具体算例 声速$ {s} $/(m·s−1) $ \mathrm{ }s=\sqrt{\dfrac{K}{\rho}} $ 光滑长度$ {h}_{i} $/m $ h_i=1.2d_{\mathrm{ini}} $ 人工黏性力系数
$ {\mathrm{\alpha }}_{\mathrm{\Pi }},\;{\mathrm{\beta }}_{\mathrm{\Pi }} $$ {\mathrm{\alpha }}_{\mathrm{\Pi }}=0.1,{\mathrm{\beta }}_{\mathrm{\Pi }}=0.1 $ 时间步长$ \Delta t $/s $ \Delta t\leqslant C_{\mathrm{cour}}\left(\dfrac{h_i}{s}\right) $ 3. 仿真结果与讨论
3.1 数值模型验证
本节通过SPH模型模拟铝颗粒的塌落过程，并与文献[27]中的实验观测结果进行对比分析。实验采用等长的圆柱形铝棒，采用2D平面应变模型进行模拟和对比。SPH模拟中使用的材料参数见表2。
表 2 铝颗粒崩塌模拟采用的材料参数Table 2. Material parameters used for aluminum particle collapse simulations参数 数值 参数 数值 密度$ \rho $/(kg·m3) 2 004.0 弹性模量$ E $/MPa 5.84 摩擦角$ \phi $/（°） 21.9 泊松比$ \upsilon $ 0.3 黏聚力$ C $/Pa 0 黏聚力设定为0，颗粒材料是非黏聚性的，代表散体颗粒的运动特性。在重力作用下，颗粒会发生塌落，这一过程受到材料剪切强度的影响。塌落过程中，重力势能转换为动能和塑性应变能。由于塑性变形，能量产生耗散，最终颗粒达到稳定状态。铝颗粒床坍塌过程中，最前端移动的距离定义为滑移距离。颗粒的塌落行为还与长宽比紧密相关，当长宽比超过某个值时，靠左侧边界的颗粒不发生移动。本文中模拟从重力加载完成后开始，重力加载20 000个时间步，时间步长设定为0.000 002 5 s。加载完成后，移除右侧障碍物，颗粒随即开始塌落，此时计时为0。模拟与实验结果对比如图3所示，左侧为实验观测结果，右侧为SPH模拟结果，x，y分别为水平、垂直方向位置。为便于与实验结果进行直观比较，对SPH模型的表面轮廓进行颜色编码。可看出模拟结果在不同时间的颗粒塌落过程与实验观测结果一致。在0.41 s时，颗粒达到滑移距离，实验和模拟均观察到了未扰动和扰动颗粒之间的清晰边界线。图3中黄线表示从SPH模拟中提取的表面形态剖面，可看出其与实验观测到的表面轮廓非常吻合。
3.2 煤岩界面演化规律分析
单个放煤口放煤SPH模型如图4所示。模型中包含煤层、矸石层及刮板输送机。其中刮板输送机采用虚粒子技术构建。刮板可移动，移动速度由刮板输送机速度决定。刮板输送机的底板为单层固定虚粒子。构成刮板输送机底座和刮板的虚粒子与SPH煤粒子、矸石粒子之间采用接触算法，当煤或矸石粒子靠近刮板壁时，虚粒子给出一个排斥力，对煤岩粒子产生推动作用[28]。图4中模型顶部煤层厚度为6.0 m，矸石层厚度为3.0 m。初始时刻放煤口关闭（构成放煤口的4层虚粒子起作用）；重力加载20 000步后，煤岩初始应力分布建立，此时打开放煤口（去除放煤口宽度范围内的4层虚粒子）。模型中初始粒子间距设置为0.2 m，共生成7 668个粒子，包含1 261个虚粒子，时间步长根据CFL准则确定为0.000 3 s。放煤高度（放煤口与刮板输送机的距离）设置为1.5 m，刮板间距为1.5 m，刮板输送机运行速度为0～1.5 m/s。
![]()
图 3 铝颗粒床在重力作用下坍塌过程的模拟和实验结果对比Figure 3. Comparison between simulation and experiment results of aluminum particle bed collapse process under influence of gravity![]()
图 4 单个放煤口放煤SPH模型Figure 4. SPH model of top coal caving with a single coal drawing outlet在沿工作面倾向方向，采煤过程中的液压支架结构不会影响采煤体的形状。因此，采用无支架的数值模型来研究放顶煤开采过程，如图4（b）所示。在沿工作面走向方向上，液压支架的采煤口宽度为2.0 m。当该支架放煤时，将模型底部边界中间2.0 m宽度范围内的虚粒子边界删除，以模拟单个采煤口的打开。当前采煤口采煤的终止条件是“见矸关门”，即检测到矸石粒子后，关闭该采煤口，通过重新布置2.0 m宽的虚粒子边界来实现。即在放煤口检测到放出的矸石粒子时，计算终止，以该时刻的煤岩界面作为终止的煤岩界面。
在顶煤被采出后，煤和岩石的界面对于单个采煤口呈现出漏斗形状，如图5（a）所示。为了研究采煤体的形状和顶煤运动轨迹，采用颗粒标记方法。顶煤被采出后，计算所有释放颗粒的编号，然后在采煤前在数值模型中标记所有具有相同编号的粒子。顶煤放出体形状如图5（b）所示，可看出其在与煤壁平行的剖面上呈现近似椭圆形状。对图5（b）中椭圆边界的粒子做标记，以追踪采煤过程中顶煤的运动轨迹。
![]()
图 5 放煤终止时的煤岩界面及放出体形状Figure 5. Coal-rock interface and shape of released body at the end of coal caving在不同的采煤步骤中，标记不同放煤时刻粒子在破碎顶煤中的位置，如图6所示。可看出在不同的放煤时刻，标记的粒子在空间位置上仍呈现近似椭圆形状，该椭圆形内的顶煤最终可通过采煤口被放出。此外，椭圆形边界外的部分顶煤颗粒虽然松动，但未被采出。
![]()
图 6 不同放煤时刻标记颗粒的位置Figure 6. The positions of marked particles at different coal drawing times整个采煤过程中所有标记颗粒的整体迁移轨迹和放出体轮廓如图7所示。由图7（a）可看出，由于顶煤开采过程中顶煤颗粒的碰撞和堵塞效应，标记的粒子并不严格沿着从顶煤椭圆形边界到采煤口的直线移动，但迁移轨迹基本上符合理论模型（图7（b））。因此，通过比较图5—图7，可认为在工作面平行方向上，单个支架的放煤模拟与理论模型基本相符。
3.3 刮板输送机运行速度对放煤过程的影响分析
放煤过程中不同时刻的粒子分布和静水压力分布云图如图8所示。放煤口打开后，刮板输送机启动运行，运行方向为x轴负方向，速度Vc为1.0 m/s。放煤开始后，放煤口附近的煤体在重力作用下迅速下落，对刮板输送机表面产生冲击，并向两侧扩散；在刮板的带动下，底煤向左运动。同时，煤岩界面（黑色线）不断向下演化。煤颗粒由放煤口不断放出，煤口底部堆积的煤被刮板输送机不断带走，由此形成了一个放煤的准稳态过程。静水压力沿垂直方向由上至下逐渐增大，呈线性分布；放煤开始后，在放煤口附近两侧的压力较高，呈现“双峰”形状，且在后续的放煤过程中形成基本稳定的应力分布。由于本文采用了应力光滑技术，静水压力及应力显示出较好的分布效果（数值噪声较小）。23.4 s时放煤口检测到矸石粒子，放煤过程终止。
![]()
图 8 顶煤放出和底煤输送过程的SPH模拟结果（$ {V}_{{\mathrm{c}}}=1.0 $ m/s）Figure 8. SPH simulation results of top coal drawing and bottom coal conveying process ($ {V}_{{\mathrm{c}}}=1.0 $ m/s)不同刮板输送机运行速度Vc下放煤口煤流速度和时间的关系曲线如图9所示。可看出煤流速度在放煤开始（20 000步）时迅速增大到峰值，这是因为放煤口打开瞬间，放煤口下方空隙使放煤口附近的煤颗粒自由落体，速度迅速增大；当放煤口下方堆积煤颗粒后，堵塞了煤流放出通道，煤流速度迅速减小并趋于较稳定的值。这说明随着支架放煤进行，协同刮板输送机运行，放煤口的煤流速度达到准稳态值。刮板输送机运行速度越大，稳定后的煤流速度越大。这也说明通过调整刮板输送机的运行速度可提高放煤速度，从而快速完成放煤口的放煤过程。
![]()
图 9 放煤口监测的煤流速度−时间曲线Figure 9. Coal flow velocity-time curve monitored at the coal drawing outlet截取不同刮板输送机运行速度下9.0 s时顶煤、矸石和底煤分布的模拟结果（图10）。可看出刮板输送机速度为0时，底煤输送停止，放煤口煤放出后在刮板输送机与支架之间的空隙内堆积，阻碍了剩余煤放出，且底煤在刮板输送机上方相对于放煤口对称分布。随着刮板输送机运行速度增大，相同时刻底煤向右输送的距离增大，煤岩界面向下演化的距离逐渐增大，更多的煤被放出，表示放煤速度提升。
![]()
图 10 不同刮板输送机运行速度下煤岩界面和底煤输送距离（t=9.0 s）Figure 10. Coal-rock interface and bottom coal conveying distance at different scraper conveyor operating speeds (t=9.0 s)1.5，9.0，18.0，23.4 s时刻4种刮板输送机运行速度（Vc=0，0.5，1.0，1.5 m/s）下煤岩界面轮廓模拟结果如图11所示。可看出随着刮板输送机运行速度增大，煤岩界面向下演化速度逐渐变大；较大的刮板输送机运行速度对应更快的放煤时间（即见矸时间），如18.0 s时1.5 m/s刮板输送机速度对应的煤岩界面已到达最下方，即已经见矸，而1.0 m/s时见矸时间为23.4 s。
![]()
图 11 不同时刻刮板输送机运行速度对煤岩界面形貌的影响Figure 11. The impact of scraper conveyor speed at different times on the morphology of coal-rock interface3.4 沿工作面方向顺序放煤模拟结果及对比
2个放煤口同时放煤时的模拟结果如图12所示。2个放煤口间距设置为6.0 m。放煤口1放出的煤颗粒在刮板输送下与放煤口2放出的煤干涉。2个放煤口的煤流速度曲线如图13所示。可看出采用2个放煤口时，放煤口1稳定后的煤流速度小于放煤口2，说明放煤口放出的底煤对放煤口产生了堵塞效果；放煤口2的煤流速度大于放煤口1，但小于单放煤口的放煤速度。由于放煤口2的煤流速度更大，所以压力“双峰”在放煤口2处更明显，如图12（a）所示。这导致煤岩界面的演化对应放煤口1和2产生了不同的规律，如图12（b）所示。这也说明在实际采煤过程中，2个放煤口同时工作要考虑底煤的干涉问题。
![]()
图 12 2个放煤口同时放煤时的模拟结果（t=11.4 s）Figure 12. Simulation results of coal being released simultaneously from two coal outlets （t=11.4 s）进一步模拟沿工作面平行方向的多放煤口放煤过程。设置工作面宽度为100 m，支架放煤口宽度为2.0 m。使用2D模型研究40个支架放煤口在1个采煤周期内的整体顶煤放出效果。模拟从左到右40个支架的放煤过程。传统的放顶煤技术采用“见矸关门”，没有岩石混入被放出的顶煤中；在自动放煤过程中，通过控制放煤时间来控制放煤口，岩石颗粒会混入被放出的顶煤中。为了进一步分析煤岩混合和顶煤采出率，收集40个支架放出的煤粒子和岩石粒子。对40个放煤口的煤岩放出量进行统计，对比分析传统放煤技术与自动放煤技术的顶煤放出量，如图14所示。自动放煤时40个放煤口的顶煤放出量平均值为17.25 m2，大于传统放煤时的16.25 m2，顶煤放出量的标准差为1.93 m2，小于传统放煤时的7.52 m2。经计算，自动放煤时顶煤采出率为85.8%，高于传统放煤时的80.9%。然而，自动放煤时矸石混合率为7.2%，高于传统放煤时的0。
![]()
图 14 不同放煤技术下每个放煤口的顶煤放出量Figure 14. Amount of top coal drawing from each coal drawing outlets under different coal drawing techniques在自动放煤过程中，一些岩石可能会混入放出的顶煤中。自动放煤时的顶煤放出量平均值更大，因此顶煤的采出率高于传统放煤技术。自动放煤时顶煤放出量的标准差较小，且每个放煤口的放出量比传统放煤更均匀，人工操作更简单，刮板输送机的运输效率更高。总体而言，自动放煤技术优于传统放煤技术。
4. 结论
1） 弹塑性土体本构模型能够有效模拟颗粒流动行为，且放煤口附近的顶煤应力呈现“双峰”分布特征。模型仅需设置摩擦角、弹性模量等基本材料参数，避免了传统粒子模型中参数选择的不确定性问题。
2） 刮板输送机的运行速度对放煤时间有显著影响，但对于“见矸关门”准则下终止煤岩界面的影响相对较小。多支架同时放煤时需考虑刮板输送机的输送能力，不同支架之间的底煤输送干涉可能导致放煤口堵塞效应。
3） “见矸关门”准则导致不同放煤口放煤量存在差异，40个放煤口顶煤放出量的标准差（7.52 m2）高于自动放煤的标准差（1.93 m2），表明该准则对放煤量的均匀性有一定影响。
-
![]()
图 1 SPH计算域离散化及粒子支持域
Figure 1. Smoothed particle hydrodynamics (SPH) computational domain discretisation and particle support domains
![]()
图 3 铝颗粒床在重力作用下坍塌过程的模拟和实验结果对比
Figure 3. Comparison between simulation and experiment results of aluminum particle bed collapse process under influence of gravity
![]()
图 4 单个放煤口放煤SPH模型
Figure 4. SPH model of top coal caving with a single coal drawing outlet
![]()
图 5 放煤终止时的煤岩界面及放出体形状
Figure 5. Coal-rock interface and shape of released body at the end of coal caving
![]()
图 6 不同放煤时刻标记颗粒的位置
Figure 6. The positions of marked particles at different coal drawing times
![]()
图 8 顶煤放出和底煤输送过程的SPH模拟结果（$ {V}_{{\mathrm{c}}}=1.0 $ m/s）
Figure 8. SPH simulation results of top coal drawing and bottom coal conveying process ($ {V}_{{\mathrm{c}}}=1.0 $ m/s)
![]()
图 9 放煤口监测的煤流速度−时间曲线
Figure 9. Coal flow velocity-time curve monitored at the coal drawing outlet
![]()
图 10 不同刮板输送机运行速度下煤岩界面和底煤输送距离（t=9.0 s）
Figure 10. Coal-rock interface and bottom coal conveying distance at different scraper conveyor operating speeds (t=9.0 s)
![]()
图 11 不同时刻刮板输送机运行速度对煤岩界面形貌的影响
Figure 11. The impact of scraper conveyor speed at different times on the morphology of coal-rock interface
![]()
图 12 2个放煤口同时放煤时的模拟结果（t=11.4 s）
Figure 12. Simulation results of coal being released simultaneously from two coal outlets （t=11.4 s）
![]()
图 14 不同放煤技术下每个放煤口的顶煤放出量
Figure 14. Amount of top coal drawing from each coal drawing outlets under different coal drawing techniques
表 1 SPH模型参数设定
Table 1 Parameters seting in SPH model
参数 取值 参数 取值 初始粒子间距
$ {d}_{{\mathrm{ini}}} $/m取决于具体算例 声速$ {s} $/(m·s−1) $ \mathrm{ }s=\sqrt{\dfrac{K}{\rho}} $ 光滑长度$ {h}_{i} $/m $ h_i=1.2d_{\mathrm{ini}} $ 人工黏性力系数
$ {\mathrm{\alpha }}_{\mathrm{\Pi }},\;{\mathrm{\beta }}_{\mathrm{\Pi }} $$ {\mathrm{\alpha }}_{\mathrm{\Pi }}=0.1,{\mathrm{\beta }}_{\mathrm{\Pi }}=0.1 $ 时间步长$ \Delta t $/s $ \Delta t\leqslant C_{\mathrm{cour}}\left(\dfrac{h_i}{s}\right) $ 
下载: 导出CSV
表 2 铝颗粒崩塌模拟采用的材料参数
Table 2 Material parameters used for aluminum particle collapse simulations
参数 数值 参数 数值 密度$ \rho $/(kg·m3) 2 004.0 弹性模量$ E $/MPa 5.84 摩擦角$ \phi $/（°） 21.9 泊松比$ \upsilon $ 0.3 黏聚力$ C $/Pa 0
下载: 导出CSV
-
[1] YANG Shengli,ZHANG Jinwang,CHEN Yi,et al. Effect of upward angle on the drawing mechanism in longwall top-coal caving mining[J]. International Journal of Rock Mechanics and Mining Sciences,2016,85:92-101. DOI: 10.1016/j.ijrmms.2016.03.004
[2] 王家臣. 我国放顶煤开采的工程实践与理论进展[J]. 煤炭学报,2018,43(1):43-51. WANG Jiachen. Engineering practice and theoretical progress of top-coal caving mining technology in China[J]. Journal of China Coal Society,2018,43(1):43-51.
[3] 王国法,庞义辉,马英. 特厚煤层大采高综放自动化开采技术与装备[J]. 煤炭工程,2018,50(1):1-6. WANG Guofa,PANG Yihui,MA Ying. Automated mining technology and equipment for fully-mechanized caving mining with large mining height in extra-thick coal seam[J]. Coal Engineering,2018,50(1):1-6.
[4] 王国法,范京道,徐亚军,等. 煤炭智能化开采关键技术创新进展与展望[J]. 工矿自动化,2018,44(2):5-12. WANG Guofa,FAN Jingdao,XU Yajun,et al. Innovation progress and prospect on key technologies of intelligent coal mining[J]. Industry and Mine Automation,2018,44(2):5-12.
[5] 王家臣,富强. 低位综放开采顶煤放出的散体介质流理论与应用[J]. 煤炭学报,2002,27(4):337-341. DOI: 10.3321/j.issn:0253-9993.2002.04.001 WANG Jiachen,FU Qiang. The loose medium flow field theory and its application on the longwall top-coal caving[J]. Journal of China Coal Society,2002,27(4):337-341. DOI: 10.3321/j.issn:0253-9993.2002.04.001
[6] 王家臣,杨建立,刘颢颢,等. 顶煤放出散体介质流理论的现场观测研究[J]. 煤炭学报,2010,35(3):353-356. WANG Jiachen,YANG Jianli,LIU haohao,et al. The practical observation research on loose medium flow field theory on the top-coal caving[J]. Journal of China Coal Society,2010,35(3):353-356.
[7] 许永祥,王国法,李明忠,等. 特厚坚硬煤层超大采高综放开采支架−围岩结构耦合关系[J]. 煤炭学报,2019,44(6):1666-1678. XU Yongxiang,WANG Guofa,LI Mingzhong,et al. Structure coupling between hydraulic roof support and surrounding rock in extra-thick and hard coal seam with super large cutting height and longwall top coal caving operation[J]. Journal of China Coal Society,2019,44(6):1666-1678.
[8] 张宁波,刘长友,陈玉明. 不稳定厚煤层放顶煤开采煤矸流场规律的数值模拟研究[J]. 煤炭技术,2014,33(12):1-4. ZHANG Ningbo,LIU Changyou,CHEN Yuming. Study on coal and gangue flow rule of top-coal caving with erratic hick coal seam by PFC2D[J]. Coal Technology,2014,33(12):1-4.
[9] 王家臣,张锦旺,杨胜利,等. 多夹矸近水平煤层综放开采顶煤三维放出规律[J]. 煤炭学报,2015,40(5):979-987. WANG Jiachen,ZHANG Jinwang,YANG Shengli,et al. 3-D movement law of top-coal in near horizontal coal seam with multi-gangue under caving mining technique[J]. Journal of China Coal Society,2015,40(5):979-987.
[10] 孙强,单成方,李亚锋,等. 浅埋双硬特厚煤层放煤规律分析及参数研究[J]. 工矿自动化,2022,48(2):61-69. SUN Qiang,SHAN Chengfang,LI Yafeng,et al. Analysis of coal drawing law and parameter research in shallow buried double hard and extra-thick coal seam[J]. Industry and Mine Automation,2022,48(2):61-69.
[11] 何欣,刘长友,吴锋锋,等. 仰斜综放开采煤层仰角对顶煤放出规律的影响研究[J]. 煤炭科学技术,2021,49(9):25-31. HE Xin,LIU Changyou,WU Fengfeng,et al. Effect of elevation angle of coal seam on top-coal caving law in fully-mechanized top-coal caving mining face during topple mining[J]. Coal Science and Technology,2021,49(9):25-31.
[12] 张锦旺,王家臣,魏炜杰. 工作面倾角对综放开采散体顶煤放出规律的影响[J]. 中国矿业大学学报,2018,47(4):805-814. ZHANG Jinwang,WANG Jiachen,WEI Weijie. Effect of face dip angle on the drawing mechanism in longwall top-coal caving mining[J]. Journal of China University of Mining & Technology,2018,47(4):805-814.
[13] 冯宇峰. 综放开采含硬夹矸顶煤破碎机理及控制技术研究[J]. 煤炭科学技术,2020,48(3):133-139. FENG Yufeng. Research on top-coal caving mechanism and control technology in extra-thick coal seam with hard dirt band[J]. Coal Science and Technology,2020,48(3):133-139.
[14] 邹光华,杨健男,关书方,等. 倾斜放煤口上下侧煤岩分界线方程表征及其模拟[J/OL]. 煤炭科学技术,1-13[2024-05-26]. http://kns.cnki.net/kcms/detail/11.2402.TD.20240227.1432.003.html. ZOU Guanghua,YANG Jiannan,GUAN Shufang,et al. Characterization and simulation of the coal-rock boundary equation on the upper and lower sides of the inclined coal caving opening[J/OL]. Coal Science and Technology, 1-13[2024-05-26]. http://kns.cnki.net/kcms/detail/11.2402.TD.20240227.1432.003.html.
[15] 姜志刚,关书方,王明强,等. 大倾角综放面煤岩分界线不对称分布特征及放煤工艺优化[J]. 煤炭技术,2023,42(12):104-108. JIANG Zhigang,GUAN Shufang,WANG Mingqiang,et al. Asymmetric distribution characteristics of coal-rock boundary and optimization of coal caving technology in fully mechanized caving face with large dip angle[J]. Coal Technology,2023,42(12):104-108.
[16] WANG Shuai,ZHANG Chunhua,HE Feng,et al. Numerical modelling of loose top coal and roof mass movement for a re-mined seam using the top coal caving method[J]. PLoS One,2023,18(4). DOI: 10.1371/JOURNAL.PONE.0283883.
[17] LIU Yang,WANG Jiachen,YANG Shengli,et al. Improving the top coal recovery ratio in longwall top coal caving mining using drawing balance analysis[J]. International Journal of Mining,Reclamation and Environment,2023,37(5):319-337. DOI: 10.1080/17480930.2023.2184036
[18] WEI Weijie,YANG Shengli,LI Meng,et al. Motion mechanisms for top coal and gangue blocks in longwall top coal caving (LTCC) with an extra-thick seam[J]. Rock Mechanics and Rock Engineering,2022,55(8):5107-5121. DOI: 10.1007/s00603-022-02928-2
[19] 张文辉,李东印,王伸,等. 特厚煤层顶煤渐进破坏的块体−颗粒耦合模拟研究[J]. 河南理工大学学报(自然科学版),2022,41(6):24-35. ZHANG Wenhui,LI Dongyin,WANG Shen,et al. Study on block-particle coupling approach for modeling progressive failure of top coal in extra thick coal seams[J]. Journal of Henan Polytechnic University (Natural Science),2022,41(6):24-35.
[20] ZHANG Qunlei,YUE Jinchao,LIU Chuang,et al. Study of automated top-coal caving in extra-thick coal seams using the continuum-discontinuum element method[J]. International Journal of Rock Mechanics and Mining Sciences,2019,122. DOI: 10.1016/j.ijrmms.2019.04.019.
[21] 李东印,王耀闯,王伸,等. 特厚煤层综放开采顶煤遗失机理及放煤步距优化数值模拟研究[J]. 河南理工大学学报(自然科学版),2023,42(2):1-10. LI Dongyin,WANG Yaochuang,WANG Shen,et al. Numerical simulation study of top-coal loss mechanism and reasonable cycle step length for extra-thick coal seam top-coal drawing[J]. Journal of Henan Polytechnic University (Natural Science),2023,42(2):1-10.
[22] 刘闯,李化敏,马占元,等. 不同煤厚综放工作面合理放煤工艺数值模拟研究[J]. 河南理工大学学报(自然科学版),2023,42(4):40-46. LIU Chuang,LI Huamin,MA Zhanyuan,et al. Numerical simulation study on reasonable top coal caving technology in longwall top coal caving working face with different top coal thickness[J]. Journal of Henan Polytechnic University (Natural Science),2023,42(4):40-46.
[23] 董祥伟. 无网格弱可压缩SPH数值算法及应用扩展[M]. 北京:科学出版社,2021. DONG Xiangwei. Extension of meshless weakly compressible SPH numerical algorithms and applications[M]. Beijing:Science Publishing House,2021.
[24] 张伟, 晏飞, 王兆丰, 等. 基于物质点和深度积分耦合模型的滑坡数值分析[J]. 岩土力学,2024,45(8):2515-2526. ZHANG Wei,YAN Fei,WANG Zhaofeng,et al. Numerical analysis oflandslides based on coupling modelof material point method and depth integral[J]. Rock and Soil Mechanics,2024,45(8):2515-2526.
[25] FENG Ruofeng,FOURTAKAS G,ROGERS B D,et al. Large deformation analysis of granular materials with stabilized and noise-free stress treatment in smoothed particle hydrodynamics (SPH)[J]. Computers and Geotechnics,2021,138. DOI: 10.1016/J.COMPGEO.2021.104356.
[26] 胡嫚,谢谟文,王立伟. 基于弹塑性土体本构模型的滑坡运动过程SPH模拟[J]. 岩土工程学报,2016,38(1):58-67. HU Man,XIE Mowen,WANG Liwei. SPH simulations of post-failure flow of landslides using elastic-plastic soil constitutive model[J]. Chinese Journal of Geotechnical Engineering,2016,38(1):58-67.
[27] NGUYEN N H T,BUI H H,NGUYEN G D. Effects of material properties on the mobility of granular flow[J]. Granular Matter,2020,22(3). DOI: 10.1007/s10035-020-01024-y.
[28] 高涛,谢守勇,胡嫚,等. 基于亚塑性本构模型的土壤−触土部件SPH互作模型[J]. 农业工程学报,2022,38(13):47-55. GAO Tao,XIE Shouyong,HU Man,et al. Soil-soil engaging component SPH model based on a hypoplastic constitutive model[J]. Transactions of the Chinese Society of Agricultural Engineering,2022,38(13):47-55.
计量
- 文章访问数: 107
- HTML全文浏览量: 27
- PDF下载量: 28
