本发明公开了一种描述堆积散体运动的数学模型的构建方法，涉及散体运动技术领域，解决了现有离散运动模型对散体运动的分析计算不够精准技术问题，其技术方案要点是根据其构建的数学模型对于宽筛分、形状不一、数量巨大的工业中的散体的运动进行计算。利用剪切系数描述粒径、形状、速度、密度、摩擦系数对散体运动的影响，没有复杂的颗粒之间相互作用的计算。该数学模型计算理论性强，数据处理快速便捷，克服了工业中大规模散体流动计算的难点，减少计算时间和硬件的要求。
1.一种描述堆积散体运动的数学模型的构建方法，其特征在于，包括：S1：通过实验测得散体颗粒内摩擦角与水分的关系，获得经验公式，根据所述经验公式得到摩擦系数与水分的关系；
S2：利用量纲分析法得到散体的粒径、形状、速度、密度、摩擦系数与剪切系数的关系式；
S3：将散体区域分布在欧拉网格下，获取各个网格内散体颗粒受到的重力、摩擦力和剪切力；其中，网格之间的剪切力表达为剪切系数和网格之间浓度梯度的函数；
S4：根据各个网格内散体受到的重力、摩擦力和剪切力，获得网格之间散体颗粒的质量交换，构建描述堆积散体运动的数学模型，根据所述数学模型得到散体颗粒的运动情况以及散体在计算区域内的浓度分布；
其中，φ表示散体颗粒内摩擦角；MC表示水分；λ1、λ2都表示实际测得的参数，根据煤粉粒径的大小而不同；
步骤S2中，利用量纲分析法得到散体的粒径、形状、速度、密度、摩擦系数与剪切系数的关系式，包括：μc＝tanφ； (3)；
其中，μ表示剪切系数；μc表示摩擦系数；ρ表示散体的堆积密度；v表示网格内散体颗粒2
的平均速度；g表示重力加速度，取9.8kg/s ；K表示经验常数且1
其中，τc表示两个网格之间的剪切力；Cl和Cq分别表示网格l和网格q内散体颗粒的体积分数；σl,q表示网格l和网格q之间的距离；
m(i,j,k)‑m(i,j,k)＝±Δmz(i,j,k±1)±Δmx(i±1,j,k)±Δmy(i,j±1,k)；
ax(i±1,j,k)＝Fx(i±1,j,k)/mx(i±1,j,k)； (6)；
ay(i,j±1,k)＝Fy(i,j±1,k)/my(i,j±1,k)； (8)；
az(i,j,k±1)＝Fz(i,j,k±1)/mz(i,j,k±1)； (10)；
Fx(i±1,j,k)＝τc(i,i±1)+f(i±1,j,k)； (12)；
Fy(i,j±1,k)＝τc(j,j±1)+f(i,j±1,k)； (13)；
Fz(i,j,k±1)＝mi,j,k+1g+f(i,j,k±1)； (14)；
mi,j,k+1＝ραi,j,k+1Vi,j,k+1； (15)；
f(i,j,k+1)＝‑μcmi,j,k+1g； (16)；
其中，i、j、k分别表示网格W(i,j,k)在x、y、z方向的标号；
Fx(i±1,j,k)、Fy(i,j±1,k)、FZ(i,j,k±1)分别表示网格W(i,j,k)在x、y、z方向的受力；Δmx(i±1,j,k)、Δmy(i,j±1,k)、ΔmZ(i,j,k±1)分别表示网格W(i,j,k)在x、y、z方向的质量变化；mx(i±1,j,k)、my(i,j±1,k)、mz(i,j,k±1)分别表示网格W(i±1,j,k)、W(i,j
±1,k)、W(i,j,k±1)内散体颗粒的质量；v x(i±1,j,k)、vy(i,j±1,k)、vz(i,j,k±1)分别表示网格W(i,j,k)内散体颗粒在x、y、z方向的初始速度；ax(i±1,j,k)、ay(i,j±1,k)、aZ(i,j,k±1)分别表示网格W(i,j,k)在x、y、z方向的加速度；mi,j,k+1表示网格W(i,j,k+1)内0
散体颗粒的质量；mi,j,k表示当前时间步网格W(i,j,k)内散体颗粒的质量，mi,j,k表示初始时间步网格W(i,j,k)内散体颗粒的质量；αi,j,k+1表示网格W(i,j,k+1)内堆积体的体积分数；
Vi±1,j,k、Vi,j±1,k、Vi,j,k±1分别表示网格W(i±1,j,k)、W(i,j±1,k)、W(i,j,k±1)的体积；f(i±1,j,k)、f(i,j±1,k)、f(i,j,k±1)分别表示网格W(i,j,k)内散体颗粒在x、y、z方向的摩擦力；τc(i,i±1)表示x方向上两个网格之间的剪切力；τc(j,j±1)表示y方向上两个网格之间的剪切力。
描述堆积散体运动的数学模型的构建方法
技术领域
[0001] 本公开涉及散体运动技术领域，尤其涉及一种描述堆积散体运动的数学模型的构建方法。
背景技术
[0002] 工业中的散体颗粒，不仅数量远远大于一般实验研究的数量级，而且粒径分布广，属于宽筛分颗粒，颗粒形状不一，表面存在不同程度的粗糙。工业中的煤仓体积庞大，常见3
规模有约9000m ，数值化煤仓管理需要对煤仓内煤料的运动进行计算，现阶段普遍用于计算颗粒运动是离散单元模型，但是该离散单元模型计算大规模的颗粒运动需要消耗大量的计算成本。并且离散单元模型需要对宽筛分的颗粒进行粒径分级处理，对形状不一的颗粒进行简化处理，从而引起了计算结果与实际情况的差异。所以传统模拟方法并不能满足工业中散体颗粒运动的计算研究，迫切需要新的数学模型满足工业中的需求。
发明内容
[0003] 本公开提供了一种描述堆积散体运动的数学模型的构建方法，其技术目的是将散体的运动描述为重力、摩擦力和剪切力的综合作用，而相邻网格之间的剪切力表达为剪切系数和网格之间浓度梯度的函数，从而得到散体大小、密度、形状、水分、表面粗糙度等对剪切力的影响，最终实现对散体运动的精准分析。
[0004] 本公开的上述技术目的是通过以下技术方案得以实现的：
[0005] 一种描述堆积散体运动的数学模型的构建方法，包括：
[0006] S1：通过实验测得散体颗粒内摩擦角与水分的关系，获得经验公式，根据所述经验公式得到摩擦系数与水分的关系；
[0007] S2：利用量纲分析法得到散体的粒径、形状、速度、密度、摩擦系数与剪切系数的关系式；
[0008] S3：将散体区域分布在欧拉网络下，获取各个网格内散体颗粒受到的重力、摩擦力和剪切力；其中，网格之间的剪切力表达为剪切系数和网格之间浓度梯度的函数；
[0009] S4：根据各个网格内散体受到的重力、摩擦力和剪切力，获得网格之间散体颗粒的质量交换，构建描述堆积散体运动的数学模型，根据所述数学模型得到散体颗粒的运动情况以及散体在计算区域内的浓度分布。
[0010] 本公开的有益效果在于：本申请所述的描述堆积散体运动的数学模型的构建方法，根据其构建的数学模型对于宽筛分、形状不一、数量巨大的工业中的散体的运动进行计算。利用剪切系数描述粒径、形状、速度、密度、摩擦系数对散体运动的影响，没有复杂的颗粒之间相互作用的计算。该数学模型计算理论性强，数据处理快速便捷，克服了工业中大规模散体流动计算的难点，减少计算时间和硬件的要求。
附图说明
[0011] 图1为本申请所述方法的流程图；
[0012] 图2为区域内的网格示意图；
[0013] 图3为网格W(i,j,k)内散体质量交换的示意图；
[0014] 图4为本申请实施例中料仓及网格划分情况示意图。
具体实施方式
[0015] 下面将结合附图对本公开技术方案进行详细说明。
[0016] 图1为本申请所述方法的流程图，如图1所示，该方法包括：
[0017] 步骤S1：通过实验测得散体颗粒内摩擦角与水分的关系，获得经验公式，根据所述经验公式得到摩擦系数与水分的关系。
[0018] 具体地，经验公式表示为：
[0019] φ＝λ1·MC+λ2； (1)；
[0020] 其中，φ表示散体颗粒内摩擦角；MC表示水分；λ1、λ2都表示实际测得的参数，根据煤粉粒径的大小而不同。
[0021] 例如，以平均粒径为41μm的煤粉为例，λ1＝204.76，λ2＝14.01。
[0022] 步骤S2：利用量纲分析法得到散体的粒径、形状、速度、密度、摩擦系数与剪切系数的关系式。
[0023] 具体地，包括：
[0024]
[0025] μc＝tanφ； (3)；
[0026] 其中，μ表示剪切系数，单位为N·m；μc表示摩擦系数；ρ表示散体的堆积密度，单位3 2
为kg/m；v表示网格内散体颗粒的平均速度，单位为m/s；g表示重力加速度，取9.8kg/s ；K表示经验常数且13
为m(米)；V表示网格体积，单位为m。
[0027] 步骤S3：将散体区域分布在欧拉网络下，获取各个网格内散体颗粒受到的重力、摩擦力和剪切力；其中，网格之间的剪切力表达为剪切系数和网格之间浓度梯度的函数。
[0028] 具体地，剪切力的计算公式为：
[0029]
[0030] 其中，τc表示两个网格之间的剪切力；Cl和Cq分别表示网格l和网格q内散体颗粒的体积分数；σl,q表示网格l和网格q之间的距离。
[0031] 步骤S4：根据各个网格内散体受到的重力、摩擦力和剪切力，获得网格之间散体颗粒的质量交换，构建描述堆积散体运动的数学模型，根据所述数学模型得到散体颗粒的运动情况以及散体在计算区域内的浓度分布。
[0032] 具体地，网格之间散体颗粒的质量交换如图3所示，具体过程包括：
[0033] m(i，j，k)‑m0(i，j，k)＝±Δmz(i，j，k±1)±Δmx(i±1，j，k)±Δmy(i，j±1，k)； (5)；
[0034] ax(i±1，j，k)＝Fx(i±1，j，k)/mx(i±1，j，k)； (6)；
[0035]
[0036] ay(i，j±1，k)＝Fy(i，j±1，k)/my(i，j±1，k)； (8)；
[0037]
[0038] az(i，j，k±1)＝Fz(i，j，k±1)/mz(i，j，k±1)； (10)；
[0039]
[0040] Fx(i±1，j，k)＝τc(i，i±1)+f(i±1，j，k)； (12)；
[0041] Fy(i，j±1，k)＝τc(j，j±1)+f(i，j±1，k)； (13)；
[0042] Fz(i，j，k±1)＝mi，j，k+1g+f(i，j，k±1)； (14)；
[0043] mi，j，k+1＝ραi，j，k+1Vi，j，k+1； (15)；
[0044] f(i，j，k+1)＝‑μcmi，j，k+1g； (16)；
[0045] 其中，i、j、k分别表示网格w(i，j，k)在x、y、z方向的标号，网格w(i，j，k)的示意图如图2所示。
[0046] Fx(i±1，j，k)、Fy(i，j±1，k)、FZ(i，j，k±1)分别表示网格W(i，j，k)在x、y、z方向的受力；Δmx(i±1，j，k)、Δmy(i，j±1，k)、ΔmZ(i，j，k±1)分别表示网格W(i，j，k)在x、y、z方向的质量变化；mx(i±1，j，k)、my(i，j±1，k)、mz(i，j，k±1)分别表示网格W(i±1，j，k)、W0 0 0
(i，j±1，k)、W(i，j，k±1)内散体颗粒的质量；vx(i±1，j，k)、vy(i，j±1，k)、v z(i，j，k±
1)分别表示网格w(i，j，k)内散体颗粒在x、y、z方向的初始速度；ax(i±1，j，k)、ay(i，j±1，k)、aZ(i,j,k±1)分别表示网格W(i,j,k)在x、y、z方向的加速度；mi,j,k+1表示网格W(i,j,k+
0
1)内散体颗粒的质量；mi,j,k表示当前时间步网格W(i,j,k)内散体颗粒的质量，mi,j,k表示初始时间步网格W(i,j,k)内散体颗粒的质量；αi,j,k+1表示网格W(i,j,k+1)内堆积体的体积分数；Vi±1,j,k、Vi,j±1,k、Vi,j,k±1分别表示网格W(i±1,j,k)、W(i,j±1,k)、W(i,j,k±1)的体积；
f(i±1,j,k)、f(i,j±1,k)、f(i,j,k±1)分别表示网格W(i,j,k)内散体颗粒在x、y、z方向的摩擦力；τc(i,i±1)表示x方向上两个网格之间的剪切力；τc(j,j±1)表示y方向上两个网格之间的剪切力。
[0047] 作为具体实施例地，以煤粉在料仓中卸料为例进行说明，如图4所示。首先对该料仓进行网格划分，在欧拉网络中计算散体的卸料过程，该料仓高1m，宽0.55m，深0.15m，卸料口大小为0.14×0.15的矩形。
[0048] 根据本申请构建的数学模型，对该料仓中散体的运动情况以及散体在计算区域内的浓度分布进行计算，具体步骤不再赘述。对于该案例，可以看到煤粉颗粒在该料仓内的卸料过程，料面呈凹形逐渐下降，直至所有散体全部流出料仓。网格尺寸为0.02m，时间步长为0.005s，使用vb软件计算整个卸料过程仅需12s，远远快于传统的离散单元模型的计算速度。
[0049] 与离散单元模型不同，本申请能够大大节约计算时间，减小对硬件的要求，不用对实际的工业中的散体进行简化，能够满足工业中大规模散体流动计算的需要。由于用重力、摩擦力、剪切力描述散体的运动，使得数据处理快速便捷，能够对工业中的散体运动进行预测。
[0050] 应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。
