CN113177346B - 一种锅炉煤粉运输安全性判断方法及系统 - Google Patents

Abstract

本发明涉及一种锅炉煤粉运输安全性判断方法及系统。该方法包括：获取煤粉运输过程中的煤颗粒数据，煤颗粒数据为煤颗粒堆积体的尺寸数据；根据尺寸数据，确定煤粉运输过程中煤颗粒堆积的粒子模型；建立描述煤颗粒介质全部流动状态的本构理论模型；对本构理论模型采用数值方法进行数值离散，得到离散方程组；根据离散方程组和煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果；对计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；根据煤颗粒流的相关信息，判断煤粉运输过程是否安全。本发明能够快速精确地判断煤粉颗粒堆是否存在安全隐患。

Description

一种锅炉煤粉运输安全性判断方法及系统

技术领域

本发明涉及锅炉煤粉运输安全性判断领域，特别是涉及一种锅炉煤粉运输安全性判断方法及系统。

背景技术

锅炉煤粉运输过程中，受外界的干扰如运输设备的震动、碰撞或地震、台风、降雨等恶劣环境，经常发生滑动、铺展和迁移，使得煤粉远离初始的堆积状态，造成煤粉与外界空气接触的面积增加，极易发生自燃甚至爆炸；同时堆积高度较高、体积较大的煤粉在滑动的过程中还会对周围的物体造成冲击损坏，部分煤粉颗粒也会上扬到空中，造成空气污染。这些都对外界环境产生了极大的安全隐患，因此，亟需一种锅炉煤粉运输安全性判断方法及系统。

发明内容

本发明的目的是提供一种锅炉煤粉运输安全性判断方法及系统，能够快速精确地判断煤粉颗粒是否会出现滑动以及出现滑动后的铺展范围。

为实现上述目的，本发明提供了如下方案：

一种锅炉煤粉运输安全性判断方法，包括：

获取煤粉运输过程中的煤颗粒数据，所述煤颗粒数据为煤颗粒堆积体的尺寸数据，所述尺寸数据包括煤颗粒堆积的长度数据、宽度数据和高度数据；

根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积的粒子模型；

建立描述煤颗粒介质全部流动状态的本构理论模型，所述煤颗粒介质全部流动状态包括类固态、类液态、类气态和惯性态；

对所述本构理论模型采用数值方法进行数值离散，得到离散方程组；

根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，所述计算结果包括场变量和位移；

对所述计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，所述相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；

根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全。

可选地，所述根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积的粒子模型，具体包括：

根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积三维几何模型；

采用网格划分软件对所述三维几何模型进行网格划分，得到粒子模型。

可选地，所述建立描述煤颗粒介质全部流动状态的本构理论模型，具体包括：

根据稀疏程度的不同，将类固态和类液态称为浓密煤颗粒流，将类气态称为稀疏煤颗粒流，将惯性态称为超稀疏煤颗粒流；

分别建立描述浓密煤颗粒流区域的理论模型、稀疏煤颗粒流区域的理论模型以及超稀疏煤颗粒流区域的理论模型；

分别建立浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则以及稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则。

可选地，所述根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，具体包括：

将煤颗粒介质初始速度设为0，初始拟温度值设为0，根据所述煤颗粒介质的粒子模型确定煤颗粒介质的初始位置；

从0时刻开始，根据所述离散方程组逐个时间步计算场变量单位时间内的变化量以及位移单位时间内的变化量；

根据所述场变量单位时间内的变化量以及位移单位时间内的变化量采用蛙跳时间更新获得所有时刻的场变量以及位移。

可选地，所述根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全，具体包括：

获取煤粉滑落的历史数据；

根据所述煤粉滑落的历史数据，得到煤粉运输过程安全的相关参数的阈值，所述煤粉运输过程安全的相关参数参数包括：煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；

判断所述煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息是否大于煤粉运输过程安全的相关参数的阈值；

若大于煤粉运输过程安全的相关参数的阈值，则煤粉运输过程存在不安全隐患；

若不大于煤粉运输过程安全的相关参数的阈值，则煤粉运输过程安全。

一种锅炉煤粉运输安全性判断系统，包括：

煤颗粒数据获取模块，用于获取煤粉运输过程中的煤颗粒数据，所述煤颗粒数据为煤颗粒堆积体的尺寸数据，所述尺寸数据包括煤颗粒堆积的长度数据、宽度数据和高度数据；

粒子模型确定模块，用于根据所述尺寸数据，确定煤粉运输过程中煤颗粒介质的粒子模型；

本构理论模型建立模块，用于建立描述煤颗粒介质全部流动状态的本构理论模型，所述煤颗粒介质全部流动状态包括类固态、类液态、类气态和惯性态；

数值离散模块，用于对所述本构理论模型采用数值方法进行数值离散，得到离散方程组；

煤颗粒介质运动过程计算模块，用于根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，所述计算结果包括场变量和位移；

后处理模块，用于对所述计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，所述相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；

判断模块，用于根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全。

可选地，所述粒子模型确定模块，具体包括：

三维几何模型建立单元，用于根据所述尺寸数据，确定煤粉运输过程中煤颗粒三维几何模型；

网格划分单元，用于采用网格划分软件对所述三维几何模型进行网格划分，得到粒子模型。

可选地，所述本构理论模型建立模块，具体包括：

稀疏程度划分单元，用于根据稀疏程度的不同，将类固态和类液态称为浓密煤颗粒流，将类气态称为稀疏煤颗粒流，将惯性态称为超稀疏煤颗粒流；

理论模型建立单元，用于分别建立描述浓密煤颗粒流区域的理论模型、稀疏煤颗粒流区域的理论模型以及超稀疏煤颗粒流区域的理论模型；

转变原则建立单元，用于分别建立浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则以及稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则。

可选地，所述煤颗粒介质运动过程计算模块，具体包括：

初始条件设置单元，用于将煤颗粒介质初始速度设为0，初始拟温度值设为0，根据所述煤颗粒介质的粒子模型确定煤颗粒介质的初始位置；

变化量第一确定单元，用于从0时刻开始，根据所述离散方程组逐个时间步计算场变量单位时间内的变化量以及位移单位时间内的变化量；

变化量第二确定单元，用于根据所述场变量单位时间内的变化量以及位移单位时间内的变化量采用蛙跳时间更新获得所有时刻的场变量以及位移。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供一种锅炉煤粉运输安全性判断方法，该判断方法包括：获取煤粉运输过程中的煤颗粒数据，煤颗粒数据为煤颗粒堆积体的尺寸数据，尺寸数据包括煤颗粒堆积的长度数据、宽度数据和高度数据；根据尺寸数据，确定煤粉运输过程中煤颗粒堆积的粒子模型；建立描述煤颗粒介质全部流动状态的本构理论模型，煤颗粒介质全部流动状态包括类固态、类液态、类气态和惯性态；对本构理论模型采用数值方法进行数值离散，得到离散方程组；根据离散方程组和煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，计算结果包括场变量和位移；对计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；根据煤颗粒流的相关信息，判断煤粉运输过程是否安全。采用本发明的上述方法能够快速精确地判断煤粉颗粒是否会出现滑动以及出现滑动后的铺展范围。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明锅炉煤粉运输安全性判断方法流程图；

图2为所建的圆柱型煤颗粒柱在自由表面上坍塌算例的粒子模型；

图3为所建煤颗粒沿滑槽坍塌滑动算例的粒子模型；

图4为煤粉颗粒介质稀疏程度划分示意图；

图5为建立整个煤颗粒介质全流态本构理论模型的步骤之间的逻辑关系示意图；

图6为煤颗粒介质全流态理论示意图；

图7为采用新发明技术计算获得的三维圆柱型颗粒柱在长径比为0.55工况下，不同时刻的煤颗粒介质运动及最终沉积形态与实验结果进行对比情况图；

图8为采用本发明计算获得的煤颗粒沿滑槽滑动过程结果图；

图9为本发明锅炉煤粉运输安全性判断系统结构图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种锅炉煤粉运输安全性判断方法及系统，能够快速精确地判断煤粉颗粒是否会出现滑动以及出现滑动后的铺展范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

传统对于煤颗粒介质的描述大多仅停留在单一流动状态的描述上，与自然界和工业中煤颗粒介质多流动状态共存并且相互之间频繁转化的实际现象存在较大差别。仅有少部分研究在现有单一流态理论的基础上进行结合获得多个流态下的理论并采用数值方法进行模拟但存在着方程迭代的问题，计算复杂且计算量大，同时这种理论也只是结合了煤颗粒介质的多个流动状态，而非全部状态，状态之间的过渡连接还不成熟，本发明就是针对存在的这些不足，创造性地提出了描述煤颗粒介质类固、类液、类气和惯性态四种流动状态的新的理论，不仅状态之间的转化更加顺畅，不需要方程迭代，同时状态之间的转化严格遵循质量、动量和能量的守恒。现有的对于煤颗粒介质多流态的计算采用了物质点法进行离散求解，计算过程中需要使用背景网格，依赖于背景网格求解动量方程，不可避免的存在网格重分、大范围布置背景网格、网格与物质点之间需要不断反复插值的缺陷，同时物质点法也无法实现对煤颗粒介质全部流动状态的模拟。本发明克服了现有数值方法的不足，采用完全无网格粒子方法进行求解，针对不同的煤颗粒流态采用不同的与之相适应的粒子方法(SPH、SDPH、DEM)，既保证对各个流态的准确描述和动态再现，又大大降低计算量，从而能够快速精确地判断锅炉煤粉运输是否安全。如图1所示，本发明提供的一种锅炉煤粉运输安全性判断方法包括以下步骤：

步骤101：获取煤粉运输过程中的煤颗粒数据，所述煤颗粒数据为煤颗粒堆积体的尺寸数据，所述尺寸数据包括煤颗粒堆积的长度数据、宽度数据和高度数据。

步骤102：根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积的粒子模型，具体包括：

根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积三维几何模型；

采用网格划分软件对所述三维几何模型进行网格划分，得到粒子模型。

煤颗粒堆积三维几何模型建立采用商业软件完成，通过商业软件建立三维几何模型，再导入网格划分软件，进行细致均匀的网格划分，最后网格文件导入到程序中按照一个网格转化成一个粒子的转化原则，全部转化为粒子，具体过程为：

1)首先建立实际煤颗粒介质材料的几何模型以及该煤颗粒介质材料所处的环境条件下的几何模型：采用达索公司旗下的SolidWorks子公司负责开发的三维软件SolidWorks软件完成。

2)在步骤1)的基础上，对于步骤1)建立的几何模型进行网格划分：对于模型网格生成，采用功能强大的CAE应用软件包—Hypermesh软件完成。

3)将步骤2)划分后形成的网格文件导入到自编程序中进行网格到粒子的转化：按照一个网格对应一个粒子的原则，采用任意六面体体积计算方法计算六面体网格的体积即为SPH粒子的体积，采用任意六面体质心计算方法计算六面体网格的质心即为SPH粒子的质心也就是SPH粒子的初始位置，这样就获得了结构的粒子模型。SPH单个粒子的体积决定了SPH单个粒子的质量，在下面步骤(5)中会涉及，是算法计算的核心数据，SPH粒子的质心，也就是SPH粒子的初始位置，直接决定了物质的初始位置。图2为所建的圆柱型煤颗粒柱在自由表面上坍塌算例的粒子模型，颜色较深的粒子表征的是处于初始堆积状态的颗粒柱，下方较浅颜色的平面粒子表征的是自由表面，提供颗粒柱坍塌的边界条件。煤颗粒柱的尺寸为半径75mm，高度为41.25mm，底部边界的尺寸为300mm×300mm，高度为一层粒子的高度，整个结构总粒子数为182484。

图3为所建煤颗粒沿滑槽坍塌滑动算例的粒子模型，颜色较深的粒子表征的是初始处于上方堆积状态的煤颗粒堆积体，下方较浅颜色的粒子表征的是滑槽的边界。初始煤颗粒堆积体的尺寸为长500mm，高为500mm，宽500mm，滑槽边界斜边长度为3500mm，高为500mm，宽500mm，下平板尺寸为长3000mm，宽2500mm，整个结构总的粒子数为46100。

步骤103：建立描述煤颗粒介质全部流动状态的本构理论模型，所述煤颗粒介质全部流动状态包括类固态、类液态、类气态和惯性态，具体包括：

根据稀疏程度的不同，将类固态和类液态称为浓密煤颗粒流，将类气态称为稀疏煤颗粒流，将惯性态称为超稀疏煤颗粒流。如图4所示。

分别建立描述浓密煤颗粒流区域的理论模型、稀疏煤颗粒流区域的理论模型以及超稀疏煤颗粒流区域的理论模型。

分别建立浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则以及稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则。如图5所示。

如图4所示，煤颗粒介质的全部流动状态共包括类固态、类液态、类气态和惯性态四种状态，根据稀疏程度的不同，将类固态和类液态的区域称之为浓密煤颗粒流区；将类气态区域称之为稀疏煤颗粒流区域；将惯性态区域称之为超稀疏煤颗粒流区域，因此建立煤颗粒介质的全部流动状态就需要根据他所处的不同状态分别进行建模，同时要建立不同态之间的联系，这也是下面每个步骤开展的原因。图5为建立整个煤颗粒介质全流态本构理论模型的步骤之间的逻辑关系示意图。

建立煤颗粒介质全部流动状态的本构理论，需要建立根据煤颗粒介质的四种流动状态所划分的三个煤颗粒流区域的本构理论，包括浓密煤颗粒流区域的理论、稀疏煤颗粒流区域的理论以及超稀疏煤颗粒流区域的理论，建立完三个流动区域的理论之后，针对不同区域之间的转变原则，再分别建立浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则以及稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则。

在整个步骤103建立本构理论获得本构方程组的基础上，采用数值方法对本构方程组进行离散获得离散方程组(步骤104)，便可以在计算机上进行程序的编写和实现了。

煤颗粒介质全流态理论示意图如图6所示，将煤颗粒介质所处的流动状态共分为三个大的区域，分别为浓密煤颗粒流区域(0.5≤ξ_p≤ξ_max，ξ_p为颗粒体积分数，ξ_max为最大体积分数，下同)、稀疏煤颗粒流区域(0.02≤ξ_p＜0.5)和超稀疏煤颗粒流区域(ξ_p＜0.02)，其中浓密煤颗粒流区域主要是指煤颗粒流处于不可压缩状态，颗粒之间以长时接触为主，颗粒的体积分数基本保持不变，主要包括颗粒的类固态和类液态两种状态；稀疏煤颗粒流区域主要是指煤颗粒流处于可压缩状态，颗粒的体积分数变化较为明显，颗粒之间以二体碰撞为主，主要包括颗粒的类液-类气之间的过渡态和类气态；超稀疏煤颗粒流区域主要是指颗粒的宏观流动特征时间尺度明显小于微观碰撞特征时间尺度，颗粒之间的二体碰撞假设不再满足，不遵循颗粒的宏观连续描述，主要包括颗粒的惯性态。

那么，对于煤颗粒介质不同的流动状态采用的本构理论是不同的：

a：对于浓密煤颗粒流区域状态采用粘弹塑性本构理论描述，具体针对浓密煤颗粒流的类固态采用弹塑性理论描述，对于浓密煤颗粒流的类液态采用流变学理论描述；

b：对于稀疏煤颗粒流区域状态采用颗粒动理学理论和摩擦动力学理论描述，具体针对稀疏煤颗粒流的类液-类气过渡态采用颗粒动理学和摩擦动力学相结合的方式描述，对于稀疏煤颗粒流的类气态采用颗粒动理学描述；

c：对于超稀疏煤颗粒流区域状态采用质点动力学理论描述，具体就是指针对超稀疏煤颗粒流的惯性态采用质点动力学理论进行描述。

下面就分别介绍这些不同的理论，以及这些理论之间的转变原则。

1)描述浓密煤颗粒流区域状态的粘-弹-塑性理论(0.5≤ξ_p≤ξ_max)

首先建立描述浓密煤颗粒流区域状态的动力学控制方程组：

α、β为自由指标，分别表示空间中的方向x,y,z，ρ为煤颗粒介质的有效密度，即空间平均密度，ρ＝ρ_p·ξ_p，ρ_p为煤颗粒介质的实际密度，ξ_p为煤颗粒介质的体积分数；v为煤颗粒介质的速度；σ^αβ为煤颗粒介质的总应力，f^α为煤颗粒介质受到的除内部应力之外的其他作用力，如重力、浮力等。

总应力张量σ^αβ写为各向同性的静水压力与应力偏张量s之和的形式：

σ^αβ＝-pδ^αβ+s^αβ (2)

δ^αβ是克罗内克函数，当α＝β时δ^αβ＝1，当α≠β时δ^αβ＝0。

静水压力p采用总应力张量分量的形式表示：

因此，求得总应力张量便可求得静水压力和应力偏张量，便可求解动力学方程组(1)。

下面就阐述总应力张量σ^αβ的求解过程。首先将总应变率划分为弹性应变率和塑性应变率两个部分：

式中，ε代表应变，

为应变率，e表示弹性分量的下标，p表示塑性分量的下标，总应变率

用公式可以直接计算出来。

v为颗粒的速度。煤颗粒介质在达到塑性之前完全处于弹性状态，遵循胡克定律，忽略准静态下的流动行为。

是剪切应力率的偏量；υ是泊松比；E是杨氏模量；G是剪切模量

σ为应力，

为应力变化率，

是3个法向应力分量的和，

δ^αβ是克罗内克函数，当α＝β时δ^αβ＝1，当α≠β时δ^αβ＝0。

在煤颗粒介质达到塑性屈服之后，按照塑性位势流动理论，塑性应变率张量采用塑性流动规律计算。

是非负的塑性标量因子，表示塑性应变增量的大小，λ为塑性乘子，Q是塑性势函数，指定了塑性应变的发展方向。由公式(3)和(5)可获得总的应变率张量以应力率的形式表示出来。

使用总应力张量的标准定义，

即

代入公式(6)中反向求解

可以获得以下增量形式的应力-应变关系式：

式中，

是偏剪切应变率张量，m、n是伪索引，在公式(7)的基础上，选取屈服准则为μ(I)流变学本构理论公式。

将公式(8)写成应力张量的第一、二不变量的形式为：

塑性流动势函数取μ(I)流变学公式(9)：

同时屈服函数f也取公式(9)，按照全微分展开可以得到：

将式(11)、(12)代入式(7)中，得到关联塑性流动规律下的应力-应变关系式：

公式(13)中的塑性流动乘子变化率公式采用以下方式获得：将式(7)代入一致性条件公式(14)中：

得到塑性流动乘子变化率公式为:

塑性势函数选取与屈服准则相同的形式(10)得到关联塑性流动规律下的塑性乘子变化率公式：

当考虑大变形问题时，本构关系必须采用与刚体转动无关的应力率。在本发明中，Jaumann应力率

采用如下形式：

为自旋速率张量，公式为：

最终获得关联型塑性流动规则下的浓密颗粒材料粘弹塑性本构理论的应力-应变关系式：

其中塑性乘子变化率公式采用(16)对应关联型塑性流动规则。使用该公式计算获得的应力应变关系应该始终保持在材料达到屈服后一直位于公式(9)所表征的屈服面上，为了实现这一要求，采用以下的退回原则执行：

如果材料在时间步长n处的应力状态超过屈服面顶点，满足以下条件时:

将法向应力分量调整为新值，以便静水压力与顶点处的静水压力相对应，调整公式如下：

式中，n代表时刻，剪切应力

保持不变。

当理想弹塑性材料发生塑性变形时，在塑性加载过程中，应力状态必须始终位于屈服面上，采用应力重标度程序将应力状态恢复到屈服面：引入了比例因子r。对于μ(I)流变学屈服准则，时间步长n处的比例因子定义为：

因此，当浓密煤颗粒介质的应力状态超过屈服面时，根据μ(I)流变学屈服准则，对应于以下条件：

偏差剪切应力分量按比例因子r减小，而静水应力分量I₁保持不变，根据以下关系：

其中

2)描述稀疏煤颗粒流区域状态的颗粒动理学理论和摩擦动力学理论(0.02≤ξ_p＜0.5)。

a：对于稀疏煤颗粒流区域中的类气态，采用颗粒动理学理论(简称KTGF)进行描述：一共包括质量守恒方程、动量守恒方程、能量守恒方程以及拟温度守恒方程四个，如下：

其中，下标p表示的是煤颗粒介质的标志，ξ_p,ρ_p和v_p分别为煤颗粒介质的体积分数、密度和速度，▽P_p为煤颗粒介质压力梯度，ξ_pρ_pf为煤颗粒介质所受外部其他作用力，本发明中该作用力主要指重力作用。h_p为煤颗粒介质的能量焓值及q_p为煤颗粒介质内部的热传导量。τ_p为煤颗粒介质粘性应力张量。

其中，μ_p和λ_p为煤颗粒介质的剪切粘度和体粘度，I为单位张量。引入对于煤颗粒介质压力P_p和粘性应力τ_p的描述，从而实现对煤颗粒介质控制方程(27)(28)的封闭。采用颗粒动理学理论获知，颗粒相压力P_p和粘性应力τ_p与颗粒速度脉动的最大值相关，而颗粒的速度脉动由颗粒拟温度描述，颗粒拟温度守恒方程如式(29)，其中(-p_pI+τ_p):▽v_p为由颗粒相应力产生的能量，k_p▽θ_p为能量耗散项，k_p为能量耗散系数，N_cθ_p为颗粒间碰撞产生的能量耗散项，具体公式为：

P_p＝ξ_pρ_p[1+2(1+e)ξ_pg₀]θ_p (33)

式中，d_p为煤颗粒的直径，e_pp为煤颗粒之间的碰撞恢复系数，

为由颗粒碰撞产生的颗粒相有效容积粘度，为中间变量。g₀为颗粒的径向恢复系数：

ξ_p,max为可压缩条件下，煤颗粒介质可达到的最大体积分数值。

b：对于稀疏煤颗粒流区域中的类液-类气之间的过渡态来说，采用颗粒动理学理论与摩擦动力学理论相结合的方式描述，如下对于颗粒相的正应力和剪切应力则分别可写为以下形式：

式中，p_KTGF和p_friction分别是煤颗粒压力p_p中的动理学部分和摩擦部分，μ_KTGF和μ_friction分别是煤颗粒剪切粘度μ_p的动理学部分和摩擦部分。ξ_p,min是摩擦应力开始逐渐增加时的颗粒体积分数。小于ξ_p,min值时，通过已知资料发现未观察到颗粒间的摩擦行为，因此假设，当均匀分布的粒子不再接触时，摩擦相互作用不再发生，主要由碰撞相互作用产生相互间的应力。对于煤颗粒由摩擦产生的正应力，采用Johnson等人提出的半经验模型，如下：

式中，Fr、a和b是材料的经验常数。对于玻璃珠，ξ_p,min、Fr、a和b的经验参数值分别取为0.5、0.05、2.0和5.0。摩擦粘度的计算最早由Schaeffer在研究筒仓颗粒物料在重力作用下从锥形出口流动的过程中，假定服从理想刚塑性本构理论，获得了摩擦剪应力与正应力之间存在以下关系式：

τ^αβ为剪切应力，p为正应力，φ为内摩擦角，

为应变率张量，

为

的2范数

总应力σ^αβ＝-Pδ^αβ+τ^αβ。Johnson等人也提出了类似的摩擦剪切力和法向应力之间服从库仑关系的公式。从公式(40)深入分析可以看出，该公式与μ(I)流变学剪切力计算公式具有相同的形式，差别在剪切力系数，这里的剪切力系数为

而流变学为μ(I)＝μ_s+(μ₂-μ_s)/(I₀/I+1)，变化范围从μ_s到μ₂。而通过实验更惊奇的发现，在摩擦动力学中，对于玻璃珠来说内摩擦角φ取值为28.5，而在流变学中，对于玻璃珠的典型值为μ_s＝tan(21°)，μ₂＝tan(33°)，从而计算得到

这也就说明了，摩擦动力学是μ(I)流变学在惯性数或者说剪切应变率无限大时的极限情况。

3)描述超稀疏煤颗粒流区域状态的质点动力学理论

当煤颗粒的宏观流动特征时间尺度明显小于微观碰撞特征时间尺度时，对煤颗粒类气态的宏观连续描述不再成立，换句话说就是对于煤颗粒相体积分数小于一定数值后，煤颗粒之间的二体碰撞假设不再满足，这时需要引入描述微观单颗粒运动和碰撞的质点动力学理论，把这种状态称为超稀疏煤颗粒流动状态，或者称为惯性态，颗粒以惯性运动为主，颗粒间碰撞作用为辅。

质点动力学是指针对具有一定质量但几何尺寸大小可以忽略的物体，采用牛顿三大定律进行描述的动力学过程。质点动力学基于牛顿第二定律，建立了质点的加速度与作用力之间关系的方程式，是质点动力学的基本模型。当质点受到n个力F_i(i＝1,2,...,n)作用时，公式为

其中r为质点矢径，上标“··”表示对时间的二阶导数。式(41)是矢量形式的微分方程，也称为质点动力学基本方程。在分析和计算实际问题时，可根据不同的坐标系将基本方程表示为相应形式的微分方程组，以便应用。

对于右端项中的F_i根据质点的实际受力来决定，如惯性力、重力、浮力、静电力、液桥力、颗粒间碰撞力、气体-颗粒间曳力、颗粒壁面碰撞力、分子键力以及范德华力等，本发明中考虑颗粒运动时所受到的惯性力、阻力、重力、浮力以及压力梯度力。

煤颗粒运动时所受到的惯性力为：

在实际的两相流动中，煤颗粒的阻力大小受到许多因素的影响，它不仅与颗粒的雷诺数Re_p有关，而且还和流体的湍流运动、流体的可压缩性、流体温度与颗粒温度、颗粒的形状、壁面的存在以及颗粒群的浓度等因素有关。因此，煤颗粒的阻力很难用统一的形式表达。为研究方便，引入阻力系数的概念，它的定义为：

于是，颗粒的阻力可表示为：

式中

为球形煤颗粒半径，ρ为流体密度，v为流体的速度。

颗粒的重力：

外部气体施加在颗粒上的浮力：

颗粒在有压力梯度的流场中运动时，还会受到由压力梯度引起的作用力。其表达式为：

式中，V_p表示煤颗粒的体积，负号表示压力梯度力的方向和流场中压力梯度的方向相反。一般来说，压力梯度力同惯性力相比数量级很小，因而可以忽略不计。

4)浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则

采用类液-类气之间增加一项过渡态的方式有效耦合起来。过渡态采用颗粒动理学与摩擦动力学相耦合的方式，摩擦动力学主要用于对接类液态或类固态计算中采用的颗粒长时接触、以摩擦应力为主的假设，因为有学者通过研究表明摩擦动力学虽然可以用于高体积分数的模拟，但是其对于较高体积分数下的状态计算存在不确定性，而有了浓密煤颗粒流的粘弹塑性本构理论之后完全克服了摩擦动力学的缺陷，同时有效结合该方法唯象的优势，在无法完全阐明物理机理的基础上能够采用简单的公式有效计算出该状态下的应力。因此，在类液态与过渡态转化的界面上，需要保持转化过程中作用力的守恒，一方面通过设定转化的体积分数值进行状态的转化，另一方面设定转化时，由浓密煤颗粒流本构计算得到的正应力与摩擦动力学计算得到的正应力相同，由浓密煤颗粒流本构计算得到的剪应力与摩擦动力学计算得到的剪应力相同。正应力相同，即：

上式中的a,b,ξ_p,min,ξ_p,max都为定值，而值Fr则由两态在相互转化时正压力值相同求得，这也既保证了两态转化过程中能量的守恒，同时也保证了Fr值选取的合理性。

对于剪应力保持相同，从公式(40)深入分析可以看出，该公式与μ(I)流变学剪切力计算公式具有相同的形式，差别在剪切力系数，这里的剪切力系数为

而流变学为μ(I)＝μ_s+(μ₂-μ_s)/(I₀/I+1)，变化范围从μ_s到μ₂。而通过实验更惊奇的发现，在摩擦动力学中，对于玻璃珠来说内摩擦角φ取值为28.5，而在流变学中，对于玻璃珠的典型值为μ_s＝tan(21°)，μ₂＝tan(33°)，从而计算得到

这也就说明了，摩擦动力学是μ(I)流变学在惯性数或者说剪切应变率无限大时的极限情况。因此，在计算剪应力的过程中直接采用μ(I)代替摩擦动力学中的

进行计算，更加保证计算的守恒性和合理性，当然在设定的转化体积分数下，在该条件的惯性数下，μ(I)数值基本趋于最大值μ₂，与摩擦动力学的

数值基本相同。

从类液态到过渡态转化的B点开始，煤颗粒动理学开始计算，颗粒的拟温度θ_p从零开始累加，表征颗粒间的碰撞效应逐渐发挥作用，到C点时，表征颗粒间的摩擦作用完全消失，完全转化为颗粒间碰撞作用，进入类气体状态。

从上述描述的过渡态的计算过程来看，这样建立的过渡态可以有效连接类液态和类气态，实现光滑过渡的目的，同时又完全符合物理规律。

5)稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则

对于颗粒体积分数小到一定数值后(<＝2％)，颗粒的宏观流动特征时间尺度明显小于微观碰撞特征时间尺度，对颗粒类气态的宏观连续描述不再成立，这时转化为颗粒质点进行追踪，由于在类气态计算的过程中，进行了宏观拟流体假设，一个单元表征了在空间该位置处颗粒的统计平均信息，如颗粒的有效密度、颗粒的均值粒径、颗粒的均值速度等，这时在转化的过程中需要保证转化的质量、动量和能量的守恒。对于稀疏煤颗粒流采用的是拉格朗日粒子方法进行离散求解，而超稀疏煤颗粒流的质点动力学也是采用粒子法进行模拟，因此转化的过程较为自然，可以保证物理量的守恒，同时为了更加贴近实际，还可以结合粒子分裂算法，使两者在空间位置上也对应起来。

步骤104：对所述本构理论模型采用数值方法进行数值离散，得到离散方程组。

1)描述浓密煤颗粒流区域状态的光滑粒子流体动力学(SPH)方法

首先采用光滑粒子流体动力学方法(SPH)对浓密煤颗粒流的粘-弹-塑性理论模型方程组(1)进行离散，获得：

式中的i,j分别是指的i粒子和j粒子，W_ij为i粒子和j粒子之间的核函数的数值，W为核函数，h为光滑长度。

对于总应力张量公式(19)同样需要采用SPH进行离散，每个粒子上均需要计算总应力张量数值，离散得到：

对于粒子i的塑性乘子

采用SPH离散后仍然为公式(51)形式：

公式(50)和(51)中的应变率和旋转应变率张量同样需要采用SPH进行离散，从而得到：

通过方程组(49)获得密度ρ、速度v^α的变化量dρ、dv^α之后，采用以下公式进行新的时刻上物理量的更新：

式中

表示粒子i的密度ρ，速度v^α；

为粒子i处在α方向上的位置坐标。

是通过方程组(49)计算获得的

和dρ_i(t)；δt为初始设定的恒定的时间步长。

2)描述稀疏煤颗粒流区域状态的光滑离散颗粒流体动力学(SDPH)方法

采用SDPH方法对稀疏煤颗粒介质区域的控制方程(27)(28)(29)进行离散，得到以下方程组：

Π_ij为粒子i和粒子j之间的人工粘性，用以下公式计算：

式中，

ε＝0.01用于防止粒子相互靠近时产生的数值发散，常数α在模拟激波时一般设定为1，在本发明中模拟煤颗粒介质问题时，α最小为0.02可以保证计算稳定。

在SDPH中，SDPH粒子的质量与其所代表的颗粒群的总质量相等，密度为颗粒群的有效密度，速度为颗粒群的均值速度，拟温度以及压力均为所代表的颗粒群的均值拟温度及均值压力，同时SDPH粒子携带表征颗粒群粒径分布特性的粒径均值、方差及单颗粒数量。

拟温度梯度▽θ_p的SPH离散公式为：

对上述方程组进行封闭所涉及的颗粒相压力及剪切力公式如(30)～(34)。

通过方程组(54)-(56)获得密度ρ、速度v^α、拟温度θ_p的变化量dρ、dv^α、dθ_p之后，采用以下公式进行新的时刻上物理量的更新：

式中

表示粒子i的密度ρ，速度v^α，拟温度θ_p；

为粒子i处在α方向上的位置坐标。

是通过方程组(54)-(56)计算获得的dρ_i(t)、

dθ_p；δt为初始设定的恒定的时间步长。

3)描述超稀疏煤颗粒流区域状态的离散单元(DEM)方法

采用离散单元法(DEM)对描述超稀疏煤颗粒流区域状态的微分方程(41)进行离散，获得：

式中，m_i与

分别为煤颗粒i的质量和速度。t为时间，m_if^α为颗粒所受到的外部作用力，

与

分别为煤颗粒i与j的接触力与粘性接触阻尼力，k_i为所有与煤颗粒接触的颗粒总数。

颗粒i与j之间的接触力

分解为：法向接触力与切向接触力，即：

法向接触力的计算采用Hertz模型：

式中：

δ_n为煤颗粒i与j接触时的侵入深度：

δ_n＝R_i+R_j-|R_j-R_i| (65)

切向接触力的计算采用Coulomb准则：

式中，μ_s为静摩擦系数，切向摩擦力的方向与相对滑动的趋势相反。

对于粘性接触阻尼力

也分解为法向与切向分量形式，即：

法向粘性接触阻尼力

的计算采用以下形式：

式中，c_n为法向粘性接触阻尼系数。

切向接触阻尼力

的计算采用以下形式：

式中，c_t为切向粘性接触阻尼系数。

有了以上作用力之后，对运动方程进行求解，可计算得到速度的时间变化量

采用公式(51)和(52)所述的蛙跳更新方法获得新的时刻的速度值以及新的时刻的颗粒的位置。

4)浓密煤颗粒流的SPH与稀疏煤颗粒流的SDPH之间算法的转化及相互作用

a.浓密煤颗粒流的SPH与稀疏煤颗粒流的SDPH之间算法的转化

浓密煤颗粒流的SPH与稀疏煤颗粒流的SDPH在表征颗粒的性质方面是相同的，不论是SPH还是SDPH不仅承载颗粒的质量、速度、位置、压力等传统参量，而且承载颗粒的粒径分布形态、体积分数以及SDPH还承载由颗粒动理学引入的拟温度等颗粒属性。SPH和SDPH在计算框架上是相同的，只不过SPH离散的是浓密煤颗粒流的粘弹塑性本构，而SDPH离散的是稀疏煤颗粒流的颗粒动理学本构。浓密煤颗粒流的SPH与稀疏煤颗粒流的SDPH之间转化主要由颗粒的体积分数值控制，而颗粒的体积分数主要由颗粒的有效密度即SPH或SDPH的密度决定，公式如下

α_p为煤颗粒体积分数，

为煤颗粒有效密度，也就是SPH粒子或SDPH粒子的密度，ρ_p为颗粒的实际密度。

的计算遵循质量守恒定律。

由SPH向SDPH转化的原则：

首先要保持煤颗粒粒子的物性参量均不变，包括粒子的位置、速度、密度、能量等，主要区别在粒子间相互作用力上。由于SPH计算的浓密煤颗粒流的正应力采用弹性定律计算，剪应力采用弹性剪应力与塑性剪应力加和的方式计算，在向SDPH转变时，这两个作用力也应该保持不变，即由SPH计算的弹性正应力转化为SDPH的摩擦正应力：

由公式(72)反向计算求得Fr作为不变量，而后根据体积分数α_s的变化更新由摩擦产生正应力的值；对于摩擦剪应力则继续采用μ(I)流变学剪应力公式计算

只不过这时的p开始由(72)式计算；对于由长时接触产生的弹性剪应力则置零；以上是对于SDPH采用摩擦动力学在过渡区产生的正应力与剪应力的数值计算，保证了转化的动量的守恒；同时SDPH从转化开始，拟温度的值由零开始计算，从而由碰撞产生的正应力和剪应力逐渐增大，直到过渡区完全部转化为颗粒动理学模型计算。

由SDPH向SPH转化的原则：

同样地，首先保持粒子的物性参量不变，包括粒子的位置、速度、密度、能量等，主要区别在粒子间相互作用力上。由SDPH计算的拟温度值反向从过渡区开始逐渐降低，两两相互碰撞的应力值逐渐降低，而摩擦正应力值p_friction逐渐增加，摩擦剪应力值同样增加，采用流变学剪应力公式计算，Fr数值按照之前的计算确定的数值，直到转化点开始，由类气态进入类液态。这时保证类液态中弹性正应力的值等于摩擦正应力的值，塑性流动剪应力的值继续按照流变学剪应力公式计算，由于剪切力逐渐减小，弹性剪切力逐渐增加，颗粒速度逐渐降低，体积分数逐渐增大，按照塑性流动法则计算煤颗粒介质材料的卸载过程，直至恢复到静止状态，由于颗粒经历了流动过程，因此煤颗粒介质无法恢复到初始状态，处于在另一个位置和状态下的准静态。

b.浓密煤颗粒流的SPH与稀疏煤颗粒流的SDPH之间算法的相互作用

在浓密煤颗粒流和稀疏煤颗粒流同时存在的状况下，两种流态颗粒之间存在着相互作用，计算两种流态的SPH粒子与SDPH粒子之间也存在着相互作用。SPH方法和SDPH方法计算均依赖于邻近粒子搜索，因此当SPH粒子与SDPH粒子作为邻近粒子时需要判定两者之间是否参与到对方的计算之中。首先以SDPH粒子为主粒子，SPH粒子作为被动粒子时，SDPH在计算以二体碰撞为假设的颗粒间的正应力和剪应力时，浓密煤颗粒流状态的颗粒对稀疏煤颗粒流状态的颗粒产生类似于二体碰撞的应力，只要两两粒子相互作用时，其中一个粒子处于稀疏状态，即可认为其和另一个粒子之间也处于二体碰撞假设范围之内，因此SPH粒子对SDPH粒子之间具有贡献，SPH粒子参与SDPH粒子的速度梯度计算；那么以SPH粒子为主粒子时，浓密煤颗粒流的计算主要是依赖于长时接触的理想弹粘塑性本构，而当SDPH作为SPH粒子的邻近粒子时，由于SDPH表征的颗粒的间距超出了长时接触的范围，因此SDPH粒子无法向SPH粒子提供长时接触的作用力，SDPH粒子不参与SPH粒子的速度梯度计算。另外，SDPH粒子计算得到的正应力和剪应力与SPH粒子计算得到的总应力相互参与到对方的计算之中。

5)稀疏煤颗粒流的SDPH与超稀疏颗粒的DEM之间算法的转化及相互作用

a.稀疏煤颗粒流的SDPH与超稀疏颗粒的DEM之间算法的转化

在稀疏煤颗粒流SDPH粒子的体积分数降到一定阈值后(<0.02)，其不再遵循二体碰撞假设的颗粒动理学模型，因此将SDPH粒子转化为DEM粒子进行计算。转化策略是，将一个SDPH粒子转化为一个DEM颗粒，SDPH粒子的质量、速度、刚度、位置等参数与转化后的DEM的粒子相同，DEM粒子的密度为实际表征的颗粒的密度，因此根据DEM粒子的质量和密度便可计算出DEM在转化后的粒径大小，即：

m_SDPH＝m_DEM (73a)

v_SDPH＝v_DEM (73b)

x_SDPH＝x_DEM (73c)

ρ_DEM＝ρ_particle (73d)

通过这种方法转化后的DEM粒子实际表征的是一个颗粒群，颗粒群中具体颗粒的数量与SDPH粒子表征的颗粒的数量相同，这种方法等效于将这些数量的实际颗粒汇聚到了现在DEM粒子的质心处，密度与实际颗粒密度相同，大小等于这些颗粒汇聚到一起之后的大小，与实际存在一定的偏差，但是在整个煤颗粒流尺寸较大、DEM粒子数较少的情况下去近似的代替求解该离散相具有一定的合理性，在保证计算量不增加的同时，提高物理模型的保真度。

b.稀疏煤颗粒流的SDPH与超稀疏颗粒的DEM之间算法的相互耦合

对于超稀疏颗粒模拟的DEM与SPH或SDPH粒子之间的相互作用问题，采用DEM颗粒之间的相互作用力法则计算。SDPH粒子按照SDPH粒子转化成DEM粒子的方法，将SDPH粒子隐形的转化成DEM粒子后计算SDPH与DEM(等效两个DEM粒子)之间的相互作用力，包括接触力F_c,ij＝F_cn,ij+F_ct,ij和法向接触阻尼力F_d,ij＝F_dn,ij+F_dt,ij，作用在SDPH和DEM两者之间的作用力大小相等、方向相反，作为动量方程源项加入到各自方程计算中，公式如下：

考虑DEM粒子对SDPH粒子作用的SDPH方法动量方程：

考虑SDPH粒子对DEM粒子作用的DEM方法的动量方程：

为DEM粒子作用于SDPH粒子上的作用力在α方向的分量，F_SDPH为DEM粒子作用于SDPH粒子上的作用力矢量。

在建立煤颗粒介质全流态数值方法的基础上，将步骤(3)中获得的守恒方程离散方程组(49)、(54)-(56)、(59)以及蛙跳时间更新方程组(51)-(52)采用C++程序进行编写，不同方法之间的转化和相互耦合策略采用步骤(3)中阐述的方法进行实现，编写完程序后采用Microsoft Visual Studio2017编译器进行编译，实现对程序的计算。

步骤105：根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，所述计算结果包括场变量和位移，具体包括：

步骤1051：将煤颗粒介质初始速度设为0，初始拟温度值设为0，根据所述煤颗粒介质的粒子模型确定煤颗粒介质的初始位置。在进行步骤1051之前英选择煤颗粒介质材料参数。例如选择煤颗粒的直径为0.32mm，实际密度为2600kg/m³，初始体积分数为0.6，体积密度为1560kg/m³，弹性模量20GPa，泊松比0.3，内摩擦角30°。采用SPH方法进行初始离散，SPH粒子的密度为颗粒的有效密度1560kg/m³，初始体积分数为0.6，SPH粒子的直径为5mm，粒子总数量根据工况不同而发生改变，光滑长度为6.5mm。采用势函数边界力计算方法，确定煤颗粒介质的边界条件。

步骤1052：从0时刻开始，根据所述离散方程组逐个时间步计算场变量单位时间内的变化量以及位移单位时间内的变化量。

步骤1053：根据所述场变量单位时间内的变化量以及位移单位时间内的变化量采用蛙跳时间更新获得所有时刻的场变量以及位移。

步骤1051-1052设置完成后，从0时刻开始，由计算机根据所编写的程序逐个时间步计算场变量ρ、v、θ_p以及位移x的单位时间内的变化量dρ、dv、dθ_p以及位移dx，而后采用蛙跳时间更新获得下一时刻的场变量ρ、v、θ_p以及位移x的数值。由此，可获得煤颗粒介质在所有时刻的场变量的数值(步骤1053)，也就是煤颗粒介质运动变化的动力学全过程数据。

步骤106：对所述计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，所述相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息。

将步骤105计算获得的数据导入到Tecplot软件中，生成各个时刻煤颗粒介质材料在空间中的分布状况图片，给出速度矢量、流态指标、剪切应变等数据在图中的显示情况，同时选取模型中特定的粒子为研究对象，提取该粒子的相关参量随时间的变化数据，通过Tecplot软件生成该粒子的相关变量的时间历程曲线。

下面选取两个典型案例对新的模型和方法的计算结果进行分析演示。

图7为采用新发明技术计算获得的三维圆柱型颗粒柱在长径比为0.55工况下，不同时刻的煤颗粒介质运动及最终沉积形态与实验结果进行对比情况图。可以看出两者在沉积轮廓、煤颗粒介质流态分布、颗粒坍塌铺展的范围、堆积的高度等都吻合较好，验证了该发明在计算煤颗粒介质全流态问题上的有效性。

图8为采用本发明计算获得的煤颗粒沿滑槽滑动过程结果图，可以看到从煤挡板释放后，煤颗粒介质逐渐坍塌沿滑槽壁面向下方运动的详细细节过程，计算结果很好的捕捉到了煤颗粒介质的运动形态、速度矢量分布、堆积过程、堆积形态等，最终获得的颗粒体铺展范围与实验结果吻合较好，再次验证了新发明适合用于煤颗粒介质全流态过程的数值模拟。

步骤107：根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全，具体包括：

步骤1071：获取煤粉滑落的历史数据。

步骤1072：根据所述煤粉滑落的历史数据，得到煤粉运输过程安全的相关参数的阈值，所述煤粉运输过程安全的相关参数包括：煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息。

步骤1073：判断所述煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息是否大于煤粉运输过程安全的相关参数的阈值。

步骤1074：若大于煤粉运输过程安全的相关参数的阈值，则煤粉运输过程存在不安全隐患，即煤粉颗粒会出现滑动。

步骤1075：若不大于煤粉运输过程安全的相关参数的阈值，则煤粉运输过程安全。

本发明一方面用来分析煤颗粒介质内部应力变化过程、屈服流动过程、快速流动过程、堆积过程等，揭示颗粒在不同流动状态下的运动机理以及不同流动状态之间的相互作用机理等，为解释煤颗粒介质材料的复杂运动现象提供支持，为理论预测模型的建立提供支撑；另一方面，能够精确迅速地判断煤粉运输过程是否安全，并通过煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息，确定煤粉发生滑落之后运动铺展的范围，从而指导相关人员在发生煤粉滑落之前做好防范，或者指导相关人员在发生煤粉滑落的实际破坏面积做好有效救援等。

如图9所示，本发明还提供一种锅炉煤粉运输安全性判断系统，该系统包括：

煤颗粒数据获取模块201，用于获取煤粉运输过程中的煤颗粒数据，所述煤颗粒数据为煤颗粒堆积体的尺寸数据，所述尺寸数据包括煤颗粒堆积的长度数据、宽度数据和高度数据。

粒子模型确定模块202，用于根据所述尺寸数据，确定煤粉运输过程中煤颗粒介质的粒子模型。

本构理论模型建立模块203，用于建立描述煤颗粒介质全部流动状态的本构理论模型，所述煤颗粒介质全部流动状态包括类固态、类液态、类气态和惯性态。

数值离散模块204，用于对所述本构理论模型采用数值方法进行数值离散，得到离散方程组。

煤颗粒介质运动过程计算模块205，用于根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，所述计算结果包括场变量和位移。

后处理模块206，用于对所述计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，所述相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息。

判断模块207，用于根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全。

所述粒子模型建立模块202，具体包括：

三维几何模型建立单元，用于根据所述尺寸数据，确定煤粉运输过程中煤颗粒三维几何模型。

网格划分单元，用于采用网格划分软件对所述三维几何模型进行网格划分，得到粒子模型。

所述本构理论模型建立模块203，具体包括：

稀疏程度划分单元，用于根据稀疏程度的不同，将类固态和类液态称为浓密煤颗粒流，将类气态称为稀疏煤颗粒流，将惯性态称为超稀疏煤颗粒流。

理论模型建立单元，用于分别建立描述浓密煤颗粒流区域的理论模型、稀疏煤颗粒流区域的理论模型以及超稀疏煤颗粒流区域的理论模型。

转变原则建立单元，用于分别建立浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则以及稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则。

所述煤颗粒介质运动过程计算模块205，具体包括：

初始条件设置单元，用于将煤颗粒介质初始速度设为0，初始拟温度值设为0，根据所述煤颗粒介质的粒子模型确定煤颗粒介质的初始位置。

变化量第一确定单元，用于从0时刻开始，根据所述离散方程组逐个时间步计算场变量单位时间内的变化量以及位移单位时间内的变化量。

变化量第二确定单元，用于根据所述场变量单位时间内的变化量以及位移单位时间内的变化量采用蛙跳时间更新获得所有时刻的场变量以及位移。

本发明一方面实现了对煤颗粒介质全部流动状态的描述：本发明可以有效描述煤颗粒介质在静止状态下的类固态、在屈服流动状态下的类液态、在快速流动状态下的类气态以及在超稀疏流动状态下的惯性态，克服了传统技术要么仅能描述单一流动状态，要么仅是其中几个状态的叠加的不足，本发明不仅对于煤颗粒介质不同流动状态下的理论进行了定义和描述，同时对于不同流动状态之间的转化方法和原则给出了具体描述，实用性好，可操作性强。本发明更加适用于实际自然界和工业工程中存在的煤颗粒流问题数值模拟，有效解释了煤颗粒介质多流动状态共存和转化过程中存在的典型现象。

本发明另一方面计算量相对较小，占用计算资源较少：本发明建立的煤颗粒介质多流态理论是从宏观连续介质力学定律出发的理论，除颗粒的惯性态不遵守连续介质力学假设之外，其他三种流态均遵循该假设，克服了传统分子动力学理论或介观离散单元理论对于煤颗粒介质单体描述带来的计算量大的问题；同时，本发明基于这些理论，建立了适合该理论的多方法耦合技术，对于每种状态采用最适合的数值方法进行离散，如煤颗粒介质类固态和类液态的SPH数值模拟方法、类气态的SDPH数值模拟方法、惯性态的DEM模拟方法，在颗粒达到惯性态之前均是基于拉格朗日粒子流体动力学进行模拟，相比传统采用的DEM方法计算量减小10倍左右，而由于颗粒达到惯性态的数量占比非常小，所以针对颗粒惯性态再采用DEM进行计算，增加的计算量非常有限。另外，相比传统的物质点法模拟煤颗粒介质多流态问题，本发明无需建立背景网格、无需在粒子和网格之间不断反复插值，计算的复杂度降低，计算量降低。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (7)

1.一种锅炉煤粉运输安全性判断方法，其特征在于，包括：

获取煤粉运输过程中的煤颗粒数据，所述煤颗粒数据为煤颗粒堆积体的尺寸数据，所述尺寸数据包括煤颗粒堆积的长度数据、宽度数据和高度数据；

根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积的粒子模型；

建立描述煤颗粒介质全部流动状态的本构理论模型，所述煤颗粒介质全部流动状态包括类固态、类液态、类气态和惯性态；

所述建立描述煤颗粒介质全部流动状态的本构理论模型，具体包括：

根据稀疏程度的不同，将类固态和类液态称为浓密煤颗粒流，将类气态称为稀疏煤颗粒流，将惯性态称为超稀疏煤颗粒流；

对于浓密煤颗粒流区域状态采用粘弹塑性本构理论描述，具体针对浓密煤颗粒流的类固态采用弹塑性理论描述，对于浓密煤颗粒流的类液态采用流变学理论描述；对于稀疏煤颗粒流区域状态采用颗粒动理学理论和摩擦动力学理论描述，具体针对稀疏煤颗粒流的类液-类气过渡态采用颗粒动理学和摩擦动力学相结合的方式描述，对于稀疏煤颗粒流的类气态采用颗粒动理学描述；对于超稀疏煤颗粒流区域状态采用质点动力学理论描述，具体针对超稀疏煤颗粒流的惯性态采用质点动力学理论进行描述；分别建立描述浓密煤颗粒流区域的理论模型、稀疏煤颗粒流区域的理论模型以及超稀疏煤颗粒流区域的理论模型；

分别建立浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则以及稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则；

对所述本构理论模型采用数值方法进行数值离散，得到离散方程组；具体采用光滑粒子流体动力学方法对浓密煤颗粒流的粘-弹-塑性理论模型方程组进行离散；采用光滑离散颗粒流体动力学方法对稀疏煤颗粒介质区域的控制方程进行离散；采用离散单元法对描述超稀疏煤颗粒流区域状态的微分方程进行离散；

根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，所述计算结果包括场变量和位移；

对所述计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，所述相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；

根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全。

2.根据权利要求1所述的锅炉煤粉运输安全性判断方法，其特征在于，所述根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积的粒子模型，具体包括：

根据所述尺寸数据，确定煤粉运输过程中煤颗粒堆积三维几何模型；

采用网格划分软件对所述三维几何模型进行网格划分，得到粒子模型。

3.根据权利要求1所述的锅炉煤粉运输安全性判断方法，其特征在于，所述根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，具体包括：

将煤颗粒介质初始速度设为0，初始拟温度值设为0，根据所述煤颗粒介质的粒子模型确定煤颗粒介质的初始位置；

从0时刻开始，根据所述离散方程组逐个时间步计算场变量单位时间内的变化量以及位移单位时间内的变化量；

根据所述场变量单位时间内的变化量以及位移单位时间内的变化量采用蛙跳时间更新获得所有时刻的场变量以及位移。

4.根据权利要求1所述的锅炉煤粉运输安全性判断方法，其特征在于，所述根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全，具体包括：

获取煤粉滑落的历史数据；

根据所述煤粉滑落的历史数据，得到煤粉运输过程安全的相关参数的阈值，所述煤粉运输过程安全的相关参数包括：煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；

判断所述煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息是否大于煤粉运输过程安全的相关参数的阈值；

若大于煤粉运输过程安全的相关参数的阈值，则煤粉运输过程存在不安全隐患；

若不大于煤粉运输过程安全的相关参数的阈值，则煤粉运输过程安全。

5.一种锅炉煤粉运输安全性判断系统，其特征在于，包括：

煤颗粒数据获取模块，用于获取煤粉运输过程中的煤颗粒数据，所述煤颗粒数据为煤颗粒堆积体的尺寸数据，所述尺寸数据包括煤颗粒堆积的长度数据、宽度数据和高度数据；

粒子模型确定模块，用于根据所述尺寸数据，确定煤粉运输过程中煤颗粒介质的粒子模型；

本构理论模型建立模块，用于建立描述煤颗粒介质全部流动状态的本构理论模型，所述煤颗粒介质全部流动状态包括类固态、类液态、类气态和惯性态；

所述本构理论模型建立模块，具体包括：

稀疏程度划分单元，用于根据稀疏程度的不同，将类固态和类液态称为浓密煤颗粒流，将类气态称为稀疏煤颗粒流，将惯性态称为超稀疏煤颗粒流；

理论模型建立单元，用于对于浓密煤颗粒流区域状态采用粘弹塑性本构理论描述，具体针对浓密煤颗粒流的类固态采用弹塑性理论描述，对于浓密煤颗粒流的类液态采用流变学理论描述；对于稀疏煤颗粒流区域状态采用颗粒动理学理论和摩擦动力学理论描述，具体针对稀疏煤颗粒流的类液-类气过渡态采用颗粒动理学和摩擦动力学相结合的方式描述，对于稀疏煤颗粒流的类气态采用颗粒动理学描述；对于超稀疏煤颗粒流区域状态采用质点动力学理论描述，具体针对超稀疏煤颗粒流的惯性态采用质点动力学理论进行描述；分别建立描述浓密煤颗粒流区域的理论模型、稀疏煤颗粒流区域的理论模型以及超稀疏煤颗粒流区域的理论模型；

转变原则建立单元，用于分别建立浓密煤颗粒流与稀疏煤颗粒流两个状态之间的转变原则以及稀疏煤颗粒流与超稀疏煤颗粒流两个状态之间的转变原则；

数值离散模块，用于对所述本构理论模型采用数值方法进行数值离散，得到离散方程组；具体采用光滑粒子流体动力学方法对浓密煤颗粒流的粘-弹-塑性理论模型方程组进行离散；采用光滑离散颗粒流体动力学方法对稀疏煤颗粒介质区域的控制方程进行离散；采用离散单元法对描述超稀疏煤颗粒流区域状态的微分方程进行离散；

煤颗粒介质运动过程计算模块，用于根据所述离散方程组和所述煤颗粒介质的粒子模型对煤颗粒介质运动过程进行计算，得到计算结果，所述计算结果包括场变量和位移；

后处理模块，用于对所述计算结果采用后处理软件Tecplot进行作图显示，得到煤颗粒流的相关信息，所述相关信息包括煤颗粒流的空间分布、速度场分布、铺展范围和堆积高度信息；

判断模块，用于根据所述煤颗粒流的相关信息，判断煤粉运输过程是否安全。

6.根据权利要求5所述的锅炉煤粉运输安全性判断系统，其特征在于，所述粒子模型确定模块，具体包括：

三维几何模型建立单元，用于根据所述尺寸数据，确定煤粉运输过程中煤颗粒三维几何模型；

网格划分单元，用于采用网格划分软件对所述三维几何模型进行网格划分，得到粒子模型。

7.根据权利要求5所述的锅炉煤粉运输安全性判断系统，其特征在于，所述煤颗粒介质运动过程计算模块，具体包括：

初始条件设置单元，用于将煤颗粒介质初始速度设为0，初始拟温度值设为0，根据所述煤颗粒介质的粒子模型确定煤颗粒介质的初始位置；

变化量第一确定单元，用于从0时刻开始，根据所述离散方程组逐个时间步计算场变量单位时间内的变化量以及位移单位时间内的变化量；

变化量第二确定单元，用于根据所述场变量单位时间内的变化量以及位移单位时间内的变化量采用蛙跳时间更新获得所有时刻的场变量以及位移。

Images