高级检索

采空区球体组装模型与多孔介质流动压损预测

陈方兴, 郝小礼, 陈世强, 李轶群, 李石林, 鲁义

陈方兴,郝小礼,陈世强,等. 采空区球体组装模型与多孔介质流动压损预测[J]. 煤炭科学技术,2025,53(5):233−242. DOI: 10.12438/cst.2024-0448
引用本文: 陈方兴,郝小礼,陈世强,等. 采空区球体组装模型与多孔介质流动压损预测[J]. 煤炭科学技术,2025,53(5):233−242. DOI: 10.12438/cst.2024-0448
CHEN Fangxing,HAO Xiaoli,CHEN Shiqiang,et al. Spherical assembly model and prediction of pressure loss of porous media flow in goaf[J]. Coal Science and Technology,2025,53(5):233−242. DOI: 10.12438/cst.2024-0448
Citation: CHEN Fangxing,HAO Xiaoli,CHEN Shiqiang,et al. Spherical assembly model and prediction of pressure loss of porous media flow in goaf[J]. Coal Science and Technology,2025,53(5):233−242. DOI: 10.12438/cst.2024-0448

采空区球体组装模型与多孔介质流动压损预测

基金项目: 

国家自然科学基金资助项目(U1361118);湖南省自然科学基金资助项目(2022JJ30251,2023JJ30243)

详细信息
    作者简介:

    陈方兴: (1996—),男,湖南株洲人,博士研究生。E-mail:weiai239@163.com

    通讯作者:

    陈世强: (1978—),男,贵州遵义人,教授,博士。E-mail:shiqiangchen@hnust.edu.cn

  • 中图分类号: TD752.2

Spherical assembly model and prediction of pressure loss of porous media flow in goaf

  • 摘要:

    快速预测采空区内部流动损失,一直是煤矿精准防灭火的关键问题。围绕采空区多孔介质流动,提出球体组装的采空区多孔介质物理模型,结合模型几何特征,找到流动的最小体元,并划分出该流动的4个阶段。基于纳维−斯托克斯(N−S)方程,建立表观速度、球径、动力黏度、流体密度与压损的关系式,得到4个阶段的压损计算式,进一步,采用时间比重的加权方法,建立了单位长度压损模型JT/L,符合Forchheimer方程形式。应用该模型预测了压损数值,对比了数值模拟结果和Ergun方程的预测值,结果表明:对数值模拟结果的影响,相比出入口效应,发现边壁效应更大,经优化确定模型为16×7的多球组装体接枝组合;对于该接枝组合内的多孔介质流动,在层流和近层流时,Ergun预测值与数值模拟结果的相对误差小于8%;引入孔隙率为0.375的多孔介质试验数据,开展无量纲分析,验证了数值模拟结果的可靠性,但是,单位长度压损模型JT/L预测值的误差最大值达69.9%,实为将N−S方程简化成一维而引起的;定义修正系数K,表征三维流动的影响,雷诺数(Re)低于10,K=2.84,雷诺数为10~3 000时,K=6.56Re−0.334;最终,得到了拟三维的JT'/L计算模型,当雷诺数低于10,预测值的误差小于6%,在雷诺数10~3 000范围内,误差低于10%。结论表明,基于球体组装模型的采空区拟三维多孔介质流动压损计算模型,可快速预测出雷诺数不超过3 000的单位长度压损,有助于科学指导采空区精准防灭火的工程实践。

    Abstract:

    Rapid prediction of flow loss in goaf is always the key problem of accurate fire prevention in coal mines. Based on the porous media flow in goaf, a physical model of porous media in goaf assembled by spheres is proposed. Combined with the geometric characteristics of the model, the minimum volume element of the flow is found, and four stages of the flow are divided. Based on Navier−Stokes equation, the relationship between apparent velocity, diameter, dynamic viscosity, fluid density and pressure loss is established, and the pressure loss calculation formula of four stages is obtained. Further, the pressure loss model per unit length JT/L is established by using the time-specific weight method, which accords with Forchheimer equation form. The model is used to predict the pressure loss, and the simulation results are compared with the predicted values of Ergun equation. The results show that the boundary effect on the numerical simulation results is greater than that of the inlet and outlet effect. The optimal model is 16×7 multi-ball assembly unit grafting combination. For the porous media flow in the grafting combination, the relative error between Ergun predicted value and simulation result is less than 8% in laminar flow or near laminar flow. The reliability of the numerical simulation results was verified through dimensionless analysis using experimental data from porous media with a porosity of 0.375. However, the maximum error of JT/L predicted value of the unit length pressure loss model is 69.9%, which is actually caused by the simplification of N−S equation to one dimension. The correction coefficient K was defined to represent the influence of three-dimensional flow. When the Reynolds number (Re) is less than 10, K=2.84, and when the Re is from 10 to 3000, K=6.56Re−0.334. Finally, the pseudo-three-dimensional JT'/L model is obtained. When the Re is below 10, the error of the predicted value is less than 6%, and when the Re is from 10 to 3000, the error is less than 10%. It is concluded that the pseudo three-dimensional porous media flow pressure loss calculation model based on the sphere assembly model can quickly predict the pressure drop per unit length of the Re less than 3000, which is helpful to scientifically guide the engineering practice of accurate fire prevention in the goaf.

  • 煤炭是我国的主要能源,在未来相当长时间仍占据重要地位[1]。随着煤矿向深部开采,采空区煤自燃现象频发,防治采空区自然发火灾害对煤矿生产安全十分重要[2-3]。现今,采空区煤自燃防灭火技术已取得显著进展,而如何防止工作面新鲜风向采空区流动正是其中的关键因素之一。因此,研究采空区近工作面堆积区中流体流动规律,揭示其流动过程中压力损失机理,建立采空区单位长度压力损失模型,对采空区煤自燃防灭火具有重要意义。

    针对采空区流动问题,已有国内外学者开展了诸多研究,并取得了一定成果。采空区结构符合多孔介质的定义[4],因此开展采空区内流场研究,大多基于多孔介质渗流理论[5-6]。黎经雷等[7]开展采空区漏风三维模拟,得到采空区自燃三带分布特征。LIU等[8]引入Sigmoid函数反映采空区“三水平区”和“三垂直区”的特征,建立了反映采空区“O”形特征的采空区孔隙度测量模型。司俊鸿等[9]通过数值模拟分析采空区渗透率三维分布模型可行性,对比实测数据得出平均误差小于8%。罗聪亮等[10]研究了工作面两端压差变化对采空区内微流动的影响。LIANG和WANG[11]引入Brinkman−Forchheimer的扩展Darcy模型,建立采空区流场数学模型。QING等[12]基于达西定律建立采空区连续多孔介质流动模型,并通过模拟与现场应用验证了理论的可靠性。

    上述研究说明运用多孔介质流动理论描述采空区内部流动是可行的,为提升采空区安全性、提高灾害防治效率亟需针对多孔介质流动阻力开展研究。最早由DARCY通过大量试验得到了达西定律,其表明层流条件下压降与速度成线性关系。FORCHHEIMER与BRINKMAN在DARCY的研究上加入惯性项以及重力项,对达西定律进行了修正。BUCHWALD等[13]开展流化床试验总结发现,多孔介质流动总阻力损失为惯性阻力损失与黏性阻力损失之和,并给出Ergun方程。此方程存在2个经验常数,A=150和B=1.75,针对这2个常数有部分学者发现其具有适用范围限制,并得出了一些不同结论。张晟等[14]开展粒径与球形度等参数对压损的影响,并通过试验修正黏性、惯性系数。LAI等[15]为了将Ergun方程适用性扩展到高孔隙率计算中,提出了一个等效校正系数。CHEN等[16]基于管球模型对压损公式进行重建,得到无经验常数计算式。KURPASKA等[17]总结了影响压损的因素排序。张四宗等[18]发现Ergun方程并不适用于不规则颗粒内流动的计算,继而利用形状因子较好的完成了流动阻力的预测。而YU等[19]引入分形理论,通过分形的方法推导得出多孔介质压力损失大小分形解。SONJIA等[20]提出了一个预测模型,用于表征高质量流量下多孔介质材料两相流动压降,并通过试验发现相对误差在30%以内。陈世强等[21-22]提出一种新的采空区多孔介质流动模型,但单周期内压力损失分布尚不清晰,且未明确计算模型的适用范围。

    针对上述问题,首先,基于一种采空区球体组装多孔介质流动模型,结合模型几何特征得到流动最小体元,并划分流动阶段,应用N−S方程对流动参数与压力损失之间的量化关系进行分析,利用时间比重加权构建单位长度压力损失计算模型。其次,开展采空区多孔介质模型压力损失数值模拟,并通过Ergun方程验证数值模拟结果的准确性。然后,将计算模型、数值模拟与Ergun方程结果进行对比,确定模型适用范围,以期实现采空区近工作面范围压力损失精准、高效预测。

    煤矿采空区是由于顶板煤体与岩层垮落所形成的一种多孔介质结构,其内部结构中存在大量微小孔洞,部分孔洞之间相互连通形成通道,在一定情况下,流体可以在其中流动。由于煤矿地质构造复杂、开采方法多样,导致采空区具有以下特征:① 自然形态丰富;② 空间结构复杂;③ 孔隙小且存在动态演化。

    1)毛细管模型。毛细管模型是将多孔介质内相连微孔通道抽象为规律排布的毛细管,可一定程度上通过迂曲度解释表观速度与真实速度的关系,但是此模型结构简单,对于流动过程中惯性力损失的表征效果较差。

    2)孔喉模型与绕流模型。孔喉模型与绕流模型间存在一定关系,在二维情况下顺排堆积的球体间孔隙可抽象为一个孔喉结构。因此为解决毛细管模型对于惯性力损失的表征问题,WU等[23]将孔喉模型与绕流模型结合,基于分形理论,得到多孔介质流动中惯性力损失分形解。

    以上模型均可在一定程度上表达多孔介质流动过程压力损失,但模型结构限制了其理论计算的适用范围。自然条件下形成的多孔介质大多具备不均匀各向异性的特性,如土壤、岩石堆积等,但其无规则特性并不利于研究开展。因此,需化繁为简从无序中寻找规律,例如著名的二维的谢尔宾斯基地毯以及三维的谢尔宾斯基海绵模型。而为了便于理论计算,常见的三维模型由等尺寸球体、正方体等顺排或插排组成,且通过试验数据验证了模型可靠性[22-23]。为了更加符合采空区自然特征,多球组装模型球体排布方式既有顺排也有插排,且球体直径存在差异,相比均匀粒径颗粒且各向同性的多孔介质模型更具优势。

    陈世强等[21]建立了一种采空区多孔介质流动模型,如图1a所示,该多球组装采空区物理模型的结构固定,由8个直径为d的大球之间填充1个直径为$\left( {\sqrt 3 - 1} \right)d$的小球以及6个直径为$\left( {\sqrt 3 - 1} \right)d$的半球,因此其孔隙率与填充颗粒粒径无关,可通过改变颗粒填充方式改变孔隙率,为方便描述压力损失推导过程,仅针对此模型开展理论推导计算。

    图  1  多球组装模型
    Figure  1.  Multi-ball assembled flow model

    此采空区多孔介质流动模型结构具有对称性,因此以模型中心坐标轴为原点切割形成8个立方体,每部分由1个大球与4个1/8的小球组成,将切割后的模型作为采空区流动最小体元,如图1b所示。

    采空区最小流动体元模型具有周期流动特性,为更好地表征流动过程,根据其结构特征划分流动阶段,如图2所示。

    图  2  流动阶段
    Figure  2.  Flow stage

    首先以图2Z方向作为流动体元的主流动方向,并根据结构特征分别得出z=0入口及4个阶段过流面积关系式。

    入口,当z=0时过流面积A0

    $$ {A_0} = {d^2} - \dfrac{1}{4} {\text{π}} {\left( {\dfrac{{\sqrt 3 - 1}}{2}d} \right)^2} $$ (1)

    第1阶段,当$ 0 \leqslant z \lt \left( {\sqrt 3 - 1} \right)0.5\;d $时过流面积Az1

    $$ \begin{gathered} {A_{z1}} = {d^2} - \dfrac{{\text{π}} }{4}\left[ {{{\left( {\dfrac{{\sqrt 3 - 1}}{2}} \right)}^2}{d^2} - {z^2}} \right]- {\text{π}} \left[ {{{\left( {\dfrac{d}{2}} \right)}^2} - {{\left( {\dfrac{d}{2} - z} \right)}^2}} \right] \end{gathered} $$ (2)

    第2阶段,当$ \left( {\sqrt 3 - 1} \right)0.5\;d \lt z \leqslant 0.5\;d $时过流面积Az2

    $$ {A_{z2}} = {d^2} - \left[ {{{\left( {\dfrac{d}{2}} \right)}^2} - {{\left( {\dfrac{d}{2} - z} \right)}^2}} \right]{\text{π}} $$ (3)

    第3阶段,当$ 0.5\;d \lt z \leqslant \left( {1.5 - 0.5\sqrt 3 } \right)d $时过流面积Az3

    $$ {A_{z3}} = {d^2} - \left[ {{{\left( {\dfrac{d}{2}} \right)}^2} - {{\left( {z - \dfrac{d}{2}} \right)}^2}} \right]{\text{π}} $$ (4)

    第4阶段,当$ \left( {1.5 - 0.5\sqrt 3 } \right)d \lt z \leqslant d $时过流面积Az4

    $$ \begin{gathered} {A_{z4}} = {d^2} - \dfrac{{3{\text{π}} }}{4}\left[ {{{\left( {\dfrac{{\sqrt 3 - 1}}{2}} \right)}^2}{d^2} - {{(d - z)}^2}} \right] -\\ {\text{π}} \left[ {{{\left( {\dfrac{d}{2}} \right)}^2} - {{\left( {z - \dfrac{d}{2}} \right)}^2}} \right] \end{gathered} $$ (5)

    在流体动力学中,对于不可压缩流体遵循动量守恒,该流动符合N−S方程,且为稳态流动,因此,在Z轴方向上不考虑质量力的影响,则有:

    $$ {u_x}\dfrac{{\partial {u_z}}}{{\partial x}} + {u_y}\dfrac{{\partial {u_z}}}{{\partial y}} + {u_z}\dfrac{{\partial {u_z}}}{{\partial z}} = - \dfrac{1}{\rho }\dfrac{{\partial P}}{{\partial z}} + \upsilon {\nabla ^2}{u_z} $$ (6)

    式中:uxuyuz为速度分量,m/s;P为流体微元上的压力,Pa;ρ为流体密度,kg/m3υ为流体运动黏度,m2/s;${\nabla ^2}$为拉普拉斯算子。

    针对第1阶段流动开展压力损失计算模型推导,z=0入口处速度为u0,即本模型表观速度,过流面积为A0;设第一阶段内实时流速为uz1,过流面积为Az1,根据连续性方程则有:

    $$ {u_{z1}} = \dfrac{{{A_0}{u_0}}}{{{A_{z1}}}} $$ (7)

    式(6)中的惯性力项与黏滞力项可用速度表达,第一阶段流体单元受到惯性力和黏滞力在Z轴方向上的分力分别为$\dfrac{{{\partial ^2}{u_{z1}}}}{{\partial {z^2}}}$和${u_{z1}}\dfrac{{\partial {u_{z1}}}}{{\partial z}}$。根据质量守恒及动量守恒,通过形心式积分得出过流断面流体Z轴方向合力大小,即:

    $$ \overline {\dfrac{{{\partial ^2}{u_z}}}{{\partial {z^2}}}} = \dfrac{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 {{ - }}1} \right)d}}{2}} {\dfrac{{{\partial ^2}{u_z}}}{{\partial {z^2}}}{A_{z1}}} {\mathrm{d}}z}}{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 {{ - }}1} \right)d}}{2}} {{A_{z1}}} {\mathrm{d}}z}} $$ (8)
    $$ \overline {{u_{z1}}\dfrac{{\partial {u_{z1}}}}{{\partial z}}} = \dfrac{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 - 1} \right)d}}{2}} {{u_{z1}}\dfrac{{\partial {u_{z1}}}}{{\partial z}}{A_{z1}}{\mathrm{d}}z} }}{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 - 1} \right)d}}{2}} {{A_{z1}}{\mathrm{d}}z} }} $$ (9)

    将式(1)、式(2)、式(3)、和式(7)代入式(8)与式(9)得:

    $$ \overline {\dfrac{{{\partial ^2}{u_z}}}{{\partial {z^2}}}} = \dfrac{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 {{ - }}1} \right)d}}{2}} {\dfrac{{{\partial ^2}{u_z}}}{{\partial {z^2}}}{A_{z1}}} {\rm{d}}z}}{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 {{ - }}1} \right)d}}{2}} {{A_{z1}}} {\rm{d}}z}} = - 7.80\dfrac{{{u_0}}}{{{d^2}}} $$ (10)
    $$ \overline {{u_{z1}}\dfrac{{\partial {u_{z1}}}}{{\partial z}}} = \dfrac{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 - 1} \right)d}}{2}} {{u_{z1}}\dfrac{{\partial {u_{z1}}}}{{\partial z}}{A_{z1}}{\rm{d}}z} }}{{\displaystyle\int_0^{\tfrac{{\left( {\sqrt 3 - 1} \right)d}}{2}} {{A_{z1}}{\rm{d}}z} }} = 11.36\dfrac{{{u_0}^2}}{d} $$ (11)

    表观速度u0方向与Z轴一致,考虑各项在XY轴方向的分量十分困难,因此将式(6)简化,只考虑式(6)主流方向Z轴的流动,将式(10)、式(11)代入式(6),则第1流动阶段的单位长度压力损失为

    $$ {J_1} = - \dfrac{{7.8\mu }}{{{d^2}}}{u_0} - \dfrac{{11.36\rho }}{d}{u_0}^2 $$ (12)

    由式(12)形式可知,符合Forchheimer方程关于初速度的二次关系的形式。同理,应用相同方法可得出各阶段单位长度压力损失计算式。

    $$ {J_2} = - \dfrac{{74.05\mu }}{{{d^2}}}{u_0} - \dfrac{{24.83\rho }}{d}{u_0}^2 $$ (13)
    $$ {J_3} = - \dfrac{{74.05\mu }}{{{d^2}}}{u_0} + \dfrac{{24.83\rho }}{d}{u_0}^2 $$ (14)
    $$ {J_4} = - \dfrac{{20.68\mu }}{{{d^2}}}{u_0} - \dfrac{{13.7\rho }}{d}{u_0}^2 $$ (15)

    推导计算得出4个阶段压力损失计算式后,通过长度比重加权求和的方法可以得到单周期流动4个阶段压力损失值,各阶段长度比重分别为0.366 d、0.134 d、0.134 d和0.366 d。由于模型断面大小随Z轴方向持续变化,流速不断在改变,而流体实际流动长度比表观长度要长许多,所以仅从长度比重考虑是不合理的。4个阶段压损大小占总压损的份额并不是表观长度比重,而是实际流动长度比重,但实际流动长度与各阶段迂曲度相关,迂曲度大多定义为实际流动长度Le与多孔介质长度L的比值(τ=Le/L),难以精确描述,于是在稳态流动情况下,流量大小不变,可通过各阶段流动时间占总时间的比值反映实际流动长度比重。因此,单周期流动过程压力损失计算式应采取时间比重加权求和的方法。

    首先,通过长度比重加权得出单周期流动4阶段流动压力损失大小表达式。

    $$ \left\{ \begin{gathered} J_1^{'} = 0.366d{J_1} = - \dfrac{{2.85\mu }}{d}{u_0} - 4.15\rho {u_0}^2 \\ J_2^{'} = 0.134d{J_2} = - \dfrac{{9.92\mu }}{d}{u_0} - 3.32\rho {u_0}^2 \\ J_3^{'} = 0.134d{J_3} = - \dfrac{{9.92\mu }}{d}{u_0} + 3.32\rho {u_0}^2 \\ J_4^{'} = 0.366d{J_4} = - \dfrac{{7.57\mu }}{d}{u_0} + 4.89\rho {u_0}^2 \\ \end{gathered} \right. $$ (16)

    分别积分求得4个阶段孔隙体积,如下式:

    $$ {V_{\mathrm{s}}} = \int_n^m {{A_z}{\rm{d}}z} $$ (17)

    式中:Vs为各阶段孔隙体积,m3Az为各阶段过流面积,m2mn为各阶段区间上下限值。

    基于连续性方程,可得出体积流量Q

    $$ Q = {A_0}{u_0} $$ (18)

    已知体积流量与各阶段孔隙体积,可分别计算出4个阶段流动时间为:

    $$ {t_{\mathrm{s}}} = \dfrac{{{V_{\mathrm{s}}}}}{Q} $$ (19)

    式中:ts为各阶段流动时间,s。

    $$ \left\{ \begin{gathered} {V_1} = \int_0^{\left( {\sqrt 3 - 1} \right)d/2} {{A_{z1}}} = {\text{0}}{\text{.181}}{d^3} \\ {V_2} = \int_{\left( {\sqrt 3 - 1} \right)d/2}^{d/2} {{A_{z2}}} = 0.031{d^3} \\ {V_3} = \int_{d/2}^{\left( {3 - \sqrt 3 } \right)d/2} {{A_{z3}}} = 0.031{d^3} \\ {V_4} = \int_{\left( {3 - \sqrt 3 } \right)d/2}^d {{A_{z4}}} = 0.130{d^3} \\ \end{gathered} \right. $$ (20)

    稳态流动体积流量不变,时间比重即为4个阶段孔隙体积比重,通过式(20)可计算得出各阶段孔隙体积。根据式(19)、式(20)计算可知,第1阶段时间比重约为0.485,第2、3阶段时间比重约为0.084,第4阶段时间比重约为0.347。则单位长度压力损失计算式如下:

    $$ \dfrac{{{J_{\rm{T}}}}}{L} = \dfrac{{378.42\mu }}{{{d^2}}}{u_0} + \dfrac{{17.93\rho }}{d}{u_0}^2 $$ (21)

    式中:L为多孔介质流动长度,m。

    孙纪国等[24]通过试验得知,在多孔介质流动中,当雷诺数小于1则压损与速度呈线性关系,大于1则呈非线性关系。因此修正雷诺数小于1应选用Laminar层流模型计算,大于1则选用RNG k-ε湍流模型计算,修正雷诺数计算式如下[25]

    $$ {R_{e}} = \dfrac{{{u_0}{d_{\mathrm{c}}}\rho }}{{\mu \left( {1 - \varepsilon } \right)}} $$ (22)

    式中:ε为孔隙率,本模型孔隙率为0.374;dc为颗粒当量直径,m。根据模型颗粒粒径参数,绘制粒径分布直方图求得模型颗粒平均粒径大小为0.01366 m。设置速度入口以及压力出口边界,入口速度最大为2 m/s。流动介质选择空气,密度ρ=1.205 m3/kg,动力黏度系数μ=1.81×10−5Pa·s。壁面为无滑移边界条件,对流项与扩散项为二阶迎风格式,各项残差收敛条件为10−5,同时监测入口与出口静压变化,当数值曲线趋于稳定则认定计算收敛。

    为了降低边壁效应以及出入口效应对模型流动压力损失计算结果的影响,针对此模型结构特征,利用三维建模软件建立了多种尺寸的采空区多孔介质流动模型,以图1a所示模型为一个多球组装体单元,针对不同尺寸模型开展计算,分析不同模型当量直径对以及纵向长度对模拟结果的影响从而确定计算模型。部分模型如图3所示,其中大球直径d=0.015 m,小球直径约为0.01098 m。图4为表观速度u0=0.1 m/s情况下双层不同组装体数量的模型计算结果分析。已知多孔介质模型边壁效应对计算结果的影响程度,取决于模型当量直径以及填充颗粒物粒径的比值,因此不断增加单层多球组装体单元数量,并开展计算分析。设置入口表观速度大小为0.1 m/s,不考虑重力的影响,计算得出1×2模型单位长度压损为52 Pa,4×2模型单位长度压损为42 Pa,9×2模型单位长度压损为37.8 Pa,16×2模型单位长度压损为36.3 Pa。从压损计算结果分析可知,压损变化率从24.7%降低至3.1%,由此可见边壁效应对于数值模拟计算结果影响较大,本多孔介质三维物理模型单层应不少于16个多球组装体。

    图  3  采空区多孔介质流动模型
    Figure  3.  Porous media flow model in goaf
    图  4  边壁效应影响分析
    Figure  4.  Analysis of the influence of boundary effect

    图5为出入口效应对计算结果的影响分析,出入口效应对计算结果的影响则取决于模型沿流向的纵向长度与颗粒粒径的比值,于是分别对2~7层多球组装体(每层具有16个多球组装体)模型开展表观速度u0=0.1 m/s情况下计算,通过结果分析得知,单位长度压力损失由36.3 Pa下降至30.7 Pa,模型计算结果变化率最大为8.8%,最小为1.98%,值得注意的是,边壁效应对结果的影响明显大于出入口效应。考虑到计算性能对网格数量的限制,最终确定16×7模型作为数值模拟模型。

    图  5  出入口效应影响分析
    Figure  5.  Analysis of the impact of entrance and exit effects

    以上分析中模型球体附近的网格大小、数量均保持一致,因此模型整体网格数量并不影响图4图5的分析结果。于是,确定16×7模型后,应针对此模型整体网格数量开展计算验证,以确保数值模拟结果精度。由于本模型结构复杂,存在大量狭小空隙,网格数量过小无法开展计算,因此选取5套网格方案进行验证,网格数量由1194万至2 000万。图6给出了入口风速为0.1 m/s时不同网格数量模型计算单位长度压损结果变化情况。很明显,当网格数量超过1200万后,单位长度压力损失的偏差均小于0.5%,考虑计算效率及结果精确度,选用1351万网格数量方案。

    图  6  不同网格数量下模型压损对比
    Figure  6.  Comparison of pressure loss under different grid numbers

    首先,为验证数值模拟结果可靠性,引入试验数据[26]开展分析。试验多孔介质由2.1 mm光滑小球压实堆积形成,流动介质为水,孔隙率0.375。此试验多孔介质孔隙率与本模型相差很小,具有较强的参考价值,但颗粒直径对于结果的影响不可忽视。为此,可开展无量纲分析,从而消除颗粒粒径因素干扰。

    若将多孔介质流动通道视为许多弯曲的毛细管,通过分析符合Forchheimer方程形式的各类计算式可得无量纲式形式为

    $$ P = \dfrac{{\Delta P}}{L}\dfrac{{\rho D_{\mathrm{p}}^3{\varepsilon ^3}}}{{{\mu ^2}{{\left( {1 - \varepsilon } \right)}^3}}} = A\tau {R_{e}} + B\tau R_{e}^2 $$ (23)

    式中:τ为迂曲度;A为黏性损失系数中的常数;B为惯性损失系数中的常数。上式出现迂曲度是考虑了孔隙通道实际长度对结果的影响,一般计算式将其影响纳入经验系数中。式(23)中τ还未确定,试验多孔介质与多球组装模型迂曲度存在差异,因此开展16×7模型数值模拟,得出结果进行分析。已知层流或近层流区黏性损失占主导地位,可忽略惯性损失,则通过变换式(23)即可计算得到两多孔介质的迂曲度。

    $$ \tau = \dfrac{P}{{A{R_{e}}}} $$ (24)

    将试验数据以及数值模拟结果代入式(24)计算,发现此试验多孔介质迂曲度约为多球组装模型迂曲度的2.1倍。为消除迂曲度不同对无量纲形式计算结果的影响,将数值模拟数据代入式(23)计算后乘2.1,与试验结果对比如图7所示。

    图  7  不同雷诺数下无量纲流动阻力
    Figure  7.  Dimensionless flow resistance under different Re

    图7可知,试验与数值模拟结果具有很好的一致性,最大偏差为7%,其中还残存孔隙率略微不同而造成的结果差异。因此,在孔隙率接近的情况下,通过数值模拟方法计算多孔介质流动压损是可行的。

    对修正雷诺数Re为0~3000范围开展计算,完成计算后得到上下端静压差与表观长度L之比,即单位长度下压力损失大小。把初始速度、空气密度ρ以及动力黏度系数μ代入式(21),得出单位长度压力损失计算值JT/L。同时,为计算结果分析提供理论依据,引入Ergun方程作为参考。其方程如下:

    $$ \dfrac{{\Delta P}}{L} = 150\dfrac{{\mu {{\left( {1 - \varepsilon } \right)}^2}}}{{{\varepsilon ^3}d_c^2}}{u_0} + 1.75\dfrac{{\rho \left( {1 - \varepsilon } \right)}}{{{\varepsilon ^3}{d_{\mathrm{c}}}}}u_0^2 $$ (25)

    式中:u0为表观速度,m/s。

    图8给出了雷诺数Re[0, 10]范围各方法单位长度压力损失计算结果,其中数值模拟与Ergun方程计算结果十分接近,相对误差不大于8%,此时惯性损失可忽略,这证明Ergun方程的黏性损失项经验系数适用于本模型。而JT/L计算式则与数值模拟相差甚远,相对误差最大为69.9%。从以上分析可知,数值模拟结果具有良好的可靠性,JT/L计算式在雷诺数处于层流区或近层流区时误差较大。观察图8单位长度压力损失曲线,发现压力损失与速度几乎呈线性关系,因此可通过系数修正直接优化JT/L计算式。以数值模拟结果为实际值,将数值模拟结果与JT/L计算结果比值,约等于2.84,作为修正系数K在雷诺数Re[0, 10]范围的取值,代入式(26)计算,图9给出JT'/L计算结果以及修正后相对误差。

    图  8  Re[0, 10]单位长度压力损失及误差分析
    Figure  8.  Analysis of pressure loss per unit length and error when Re[0, 10]
    图  9  Re[0, 10]修正结果及误差分析
    Figure  9.  Correction results and error analysis when Re[0, 10]
    $$ {\dfrac{{{J_{\rm{T}}}}}{L}'} = K \left( {\dfrac{{378.42\mu }}{{{d^2}}}{u_0} + \dfrac{{17.93\rho }}{d}{u_0}^2} \right) $$ (26)

    图9可知,通过拟三维系数K修正JT/L计算式,其单位长度压力损失的趋势、幅度以及数值均符合数值模拟结果,且相对误差最小为0.8%,最大不超过5%。在雷诺数Re[0, 10]范围内,其流动损失过程黏性力占主导,力的方向应在模型与流体接触切面内,且与速度方向相反,而推导过程并未考虑速度分量,Z轴方向速度与模型表面形成的夹角导致计算分量小于实际值,从而造成图8JT/L计算结果偏小的现象。修正系数K优化JT/L计算式,解决了其中黏性损失小于实际值的问题,形成了新的单位长度压力损失拟三维计算模型JT'/L。计算模型JT'/L在层流区(Re<1)以及近层流区(1<Re<10),针对单位长度压损具备良好的预测能力,在Re>10情况下,流动损失由黏性损失主导向惯性损失主导转变,拟三维修正系数同样发生转变。

    图10为雷诺数Re[10, 3000]范围各方法单位长度压力损失计算结果,随着雷诺数Re的增大,各单位长度压损结果总体增长趋势相近,均与速度的二次方呈正相关。但是,当雷诺数大于150后,数值模拟结果小于JT/L计算式与Ergun方程,且差值持续增大。从表面来看,其原因为数值模拟曲线二阶函数系数小于二者,而从物理意义出发则存在以下原因:其一,JT/L计算式系数因推导过程忽略XY轴方向分量,由速度项数值过大而引起的系数不准确;其二,Ergun方程是基于均匀颗粒床试验得出,非均匀颗粒取得平均粒径代入计算造成计算误差出现;其三,多孔介质模型内颗粒排布方式不同则流动过程产生差异,同样将导致Ergun方程计算结果误差。由于颗粒粒径以及排布方式差异导致两模型阻力系数不同,所以Ergun方程计算误差原因可总结为方程中的经验系数不适配本模型,当异形颗粒非球体且大小不一,则两经验系数均不适配[14],而通过分析图8图10结果可知,当颗粒为粒径不一的标准球体时,黏性损失系数适配,惯性损失系数不适配。则Ergun方程适用于预测非均匀球形颗粒多孔介质层流或近层流区的流动压损,且Re超出280后误差较大。为提高JT/L计算式结果准确性,应开展计算式优化,分析数值模拟、JT/L计算结果比值与雷诺数的关系,通过数据拟合后得出修正系数K表达式,以修正雷诺数Re[10, 3000]范围JT/L计算结果。

    图  10  Re[10, 3000]单位长度压力损失
    Figure  10.  Pressure loss per unit length when Re[10, 3000]

    图11为修正系数K拟合曲线,系数最大为2.69,整体随雷诺数增大呈递减趋势,当雷诺数大于1000后递减幅度逐步降低,符合多孔介质流动雷诺数大于1000后压损与速度的平方呈线性关系规律,由此可知,本多孔介质流动在Re[10, 1000]范围属于转捩区,Re[1000, 3000]范围属于湍流区。图11表示,雷诺数10~3000情况下修正系数与雷诺数的关联式为幂函数,且R2为0.99,则修正系数K

    图  11  Re[10, 3000]修正系数拟合曲线
    Figure  11.  Fitting curve of the correction coefficient when Re[0, 10]
    $$ K = 6.56{R_{e}}^{ - 0.334} $$ (27)

    将式(27)代入式(26),图12给出JT'/L计算结果以及修正后相对误差。其中JT'/L压力损失变化规律符合数值模拟结果特征,二者随雷诺数变化趋势相同。雷诺数大于2 000,JT'/L计算式相对误差逐渐增大,在雷诺数Re[10, 3000]区间内,JT'/L计算结果与数值模拟相对误差低于10%。

    图  12  Re[10, 3000]修正结果及误差分析
    Figure  12.  Correction results and error analysis when Re[10, 3000]

    通过以上结果分析结合式(26)可知,拟三维修正系数K在层流区和近层流区是对表观速度与流动实际速度关系的修正,其中核心是关于表观长度与实际流动长度的关系。而系数K在转捩区和湍流区则描述了随着雷诺数增大流动状态由层流向湍流转变的过程,主要针对流动实际速度以及速度平方下系数进行修正。黏性项和惯性项的常数系数则是对多孔介质孔隙率、迂曲度以及阻力系数等参数的描述。

    1)基于采空区多球组装流动模型,按照时间比重加权求和,得出新的模型周期流动单位长度压力损失表达式JT/L

    2)通过定义修正系数K优化计算式,得到单位长度压力损失拟三维计算模型JT'/L。当雷诺数低于10,修正系数K=2.84,解决了层流、近层流状态下黏性损失小于实际值的问题;当雷诺数处于10~3000区间,修正系数K=6.56Re−0.334,代入模型后消除了计算式推导过程考虑一维流动产生的误差影响。

    3)JT'/L计算模型存在适用范围,在雷诺数0~10区间内计算误差小于6%,在雷诺数10~3000范围内计算误差小于10%,证明修正后的计算式JT'/L在雷诺数低于3000情况下其计算结果具备可靠性。

  • 图  1   多球组装模型

    Figure  1.   Multi-ball assembled flow model

    图  2   流动阶段

    Figure  2.   Flow stage

    图  3   采空区多孔介质流动模型

    Figure  3.   Porous media flow model in goaf

    图  4   边壁效应影响分析

    Figure  4.   Analysis of the influence of boundary effect

    图  5   出入口效应影响分析

    Figure  5.   Analysis of the impact of entrance and exit effects

    图  6   不同网格数量下模型压损对比

    Figure  6.   Comparison of pressure loss under different grid numbers

    图  7   不同雷诺数下无量纲流动阻力

    Figure  7.   Dimensionless flow resistance under different Re

    图  8   Re[0, 10]单位长度压力损失及误差分析

    Figure  8.   Analysis of pressure loss per unit length and error when Re[0, 10]

    图  9   Re[0, 10]修正结果及误差分析

    Figure  9.   Correction results and error analysis when Re[0, 10]

    图  10   Re[10, 3000]单位长度压力损失

    Figure  10.   Pressure loss per unit length when Re[10, 3000]

    图  11   Re[10, 3000]修正系数拟合曲线

    Figure  11.   Fitting curve of the correction coefficient when Re[0, 10]

    图  12   Re[10, 3000]修正结果及误差分析

    Figure  12.   Correction results and error analysis when Re[10, 3000]

  • [1] 谢和平,吴立新,郑德志. 2025年中国能源消费及煤炭需求预测[J]. 煤炭学报,2019,44(7):1949−1960.

    XIE Heping,WU Lixin,ZHENG Dezhi. Prediction on the energy consumption and coal demand of China in 2025[J]. Journal of China Coal Society,2019,44(7):1949−1960.

    [2] 邓军,李贝,王凯,等. 我国煤火灾害防治技术研究现状及展望[J]. 煤炭科学技术,2016,44(10):1−7,101.

    DENG Jun,LI Bei,WANG Kai,et al. Research status and outlook on prevention and control technology of coal fire disaster in China[J]. Coal Science and Technology,2016,44(10):1−7,101.

    [3] 邓军,杨囡囡,王彩萍,等. 采空区煤自燃“防−抑−灭” 协同防灭火关键技术[J]. 煤矿安全,2022,53(9):1−8.

    DENG Jun,YANG Nannan,WANG Caiping,et al. Key technology of“preventing-suppressing-extinguishing” coordinated fire preventing and extinguishing for coal spontaneous combustion in goaf[J]. Safety in Coal Mines,2022,53(9):1−8.

    [4] 汪腾蛟,聂朝刚,杨小彬,等. 考虑温度变化的采空区瓦斯抽采数值模拟[J]. 煤炭科学技术,2021,49(7):85−94.

    WANG Tengjiao,NIE Chaogang,YANG Xiaobin,et al. Numerical simulation of gas drainage in gob considering temperature change[J]. Coal Science and Technology,2021,49(7):85−94.

    [5] 徐树媛,张永波,时红,等. 采空区冒落带内破碎岩体的渗流特征与渗透性试验研究[J]. 安全与环境工程,2022,29(1):128−134.

    XU Shuyuan,ZHANG Yongbo,SHI Hong,et al. Flow characteristics and experimental study on the permeability of mining-induced fractured rock mass in caving zones[J]. Safety and Environmental Engineering,2022,29(1):128−134.

    [6] 丁洋,陈文彬,林海飞,等. 煤矿采空区碳封存CO2泄漏地表扩散规律研究[J]. 西安科技大学学报,2023,43(4):705−714.

    DING Yamg,CHEN Wenbin,LIN Haifei,et al. study on surface diffusion law of carbon storage CO2 leakage in coal mine goaf[J]. Jounal of Xi’an University of Science and Technology,2023,43(4):705−714.

    [7] 黎经雷,牛会永,鲁义,等. 风速对近距离煤层采空区漏风及煤自燃影响研究[J]. 煤炭科学技术,2019,47(3):156−162.

    LI Jinglei,NIU Huiyong,LU Yi,et al. Study on effect of wind speed to air leakage and spontaneous combustion in goaf of contiguous seams[J]. Coal Science and Technology,2019,47(3):156−162.

    [8]

    LIU Q,LIN B Q,ZHOU Y,et al. Porosity model of the goaf based on overlying strata movement and deformation[J]. Environmental Earth Sciences,2022,81(7):214. doi: 10.1007/s12665-022-10329-5

    [9] 司俊鸿,程根银,朱建芳,等. 采空区非均质多孔介质渗透特性三维建模及应用[J]. 煤炭科学技术,2019,47(5):220−224.

    SI Junhong,CHENG Genyin,ZHU Jianfang,et al. Three-dimensional modeling and application of permeability characteristics of heterogeneous porous media in goaf[J]. Coal Science and Technology,2019,47(5):220−224.

    [10] 罗聪亮,王海桥,陈世强,等. 压差对采空区空气微流动场影响的模拟研究[J]. 黑龙江科技大学学报,2016,26(1):1−4. doi: 10.3969/j.issn.2095-7262.2016.01.001

    LUO Congliang,WANG Haiqiao,CHEN Shiqiang,et al. Simulation on effects of micro-flow air in goaf on changing pressure differences[J]. Journal of Heilongjiang University of Science and Technology,2016,26(1):1−4. doi: 10.3969/j.issn.2095-7262.2016.01.001

    [11]

    LIANG Y T,WANG S G. Prediction of coal mine goaf self-heating with fluid dynamics in porous media[J]. Fire Safety Journal,2017,87:49−56. doi: 10.1016/j.firesaf.2016.12.002

    [12]

    QIN B T,WANG H T,YANG J Z,et al. Large-area goaf fires:A numerical method for locating high-temperature zones and assessing the effect of liquid nitrogen fire control[J]. Environmental Earth Sciences,2016,75(21):1396. doi: 10.1007/s12665-016-6173-5

    [13]

    BUCHWALD T,SCHMANDRA G,SCHÜTZENMEISTER L,et al. Gaseous flow through coarse granular beds:The role of specific surface area[J]. Powder Technology,2020,366:821−831. doi: 10.1016/j.powtec.2020.03.028

    [14] 张晟,张晓虎,赵亮,等. 基于Ergun方程的菱镁球团填充床层阻力特性实验[J]. 东北大学学报(自然科学版),2021,42(3):347−352.

    ZHANG Sheng,ZHANG Xiaohu,ZHAO Liang,et al. Experiment of resistance characteristics for magnesite pellets packed bed based on Ergun equation[J]. Journal of Northeastern University (Natural Science),2021,42(3):347−352.

    [15]

    LAI T W,LIU X Y,XUE S,et al. Extension of Ergun equation for the calculation of the flow resistance in porous media with higher porosity and open-celled structure[J]. Applied Thermal Engineering,2020,173:115262. doi: 10.1016/j.applthermaleng.2020.115262

    [16]

    CHEN Z D,ZHENG K C,WANG X,et al. Derivation and validation of a new porous media resistance formula based on a tube-sphere model[J]. Industrial & Engineering Chemistry Research,2020,59(40):18170−18179.

    [17]

    KURPASKA S,KIEŁBASA P,SOBOL Z,et al. Analysis of air flow resistance through a porous stone bed[J]. Biosystems Engineering,2020,198:323−337. doi: 10.1016/j.biosystemseng.2020.08.002

    [18] 张四宗,温治,刘训良,等. 颗粒形状对烧结矿填充床内渗透系数和阻力系数的影响[J]. 中南大学学报(自然科学版),2021,52(4):1066−1075. doi: 10.11817/j.issn.1672-7207.2021.04.003

    ZHANG Sizong,WEN Zhi,LIU Xunliang,et al. Effects of particle shape on permeability and resistance coefficients of sinter packed bed[J]. Journal of Central South University (Science and Technology),2021,52(4):1066−1075. doi: 10.11817/j.issn.1672-7207.2021.04.003

    [19]

    YU B M,LI J H. Some fractal characters of porous media[J]. Fractals,2001,9(3):365−372. doi: 10.1142/S0218348X01000804

    [20]

    WEISE S,MEINICKE S,WETZEL T,et al. Modelling the pressure drop of two-phase flow through solid porous media[J]. International Journal of Multiphase Flow,2019,112:13−26. doi: 10.1016/j.ijmultiphaseflow.2018.12.005

    [21] 陈世强,张连会,姜文,等. 采空区气体流动:Forchheimer模型改进与压力损失预测[J]. 煤炭学报,2020,45(S1):361−366.

    CHEN Shiqiang,ZHANG Lianhui,JIANG Wen,et al. Gas flow in goaf:Improvement of Forchheimer model and prediction of pressure loss[J]. Journal of China Coal Society,2020,45(S1):361−366.

    [22] 张连会,陈世强,赵利群,等. 采空区最小流动单元压力损失计算与试验[J]. 中国安全科学学报,2021,31(12):78−84.

    ZHANG Lianhui,CHEN Shiqiang,ZHAO Liqun,et al. Experiment and calculation of pressure loss of minimum flow unit in goaf[J]. China Safety Science Journal,2021,31(12):78−84.

    [23]

    WU J S,YU B M,YUN M J. A resistance model for flow through porous media[J]. Transport in Porous Media,2008,71(3):331−343. doi: 10.1007/s11242-007-9129-0

    [24] 孙纪国,王建华. 烧结多孔结构的渗透和流阻特性研究[J]. 航空动力学报,2008,23(1):130−133.

    SUN Jiguo,WANG Jianhua. Research on permeability and flow resistance characteristic of sintered porous structure[J]. Journal of Aerospace Power,2008,23(1):130−133.

    [25]

    GRACIANO-URIBE J,PUJOL T,PUIG-BARGUÉS J,et al. Assessment of different pressure drop-flow rate equations in a pressurized porous media filter for irrigation systems[J]. Water,2021,13(16):2179. doi: 10.3390/w13162179

    [26]

    SALAHI M B,SEDGHI-ASL M,PARVIZI M. Nonlinear flow through a packed-column experiment[J]. Journal of Hydrologic Engineering,2015,20(9):04015003. doi: 10.1061/(ASCE)HE.1943-5584.0001166

图(12)
计量
  • 文章访问数:  21
  • HTML全文浏览量:  2
  • PDF下载量:  12
  • 被引次数: 0
出版历程
  • 收稿日期:  2024-04-09
  • 网络出版日期:  2025-05-14
  • 刊出日期:  2025-05-24

目录

/

返回文章
返回