非牛顿流体管道流动方程及其近似解

2020-12-14兰文涛吴爱祥

哈尔滨工业大学学报 2020年12期

兰文涛，吴爱祥

(1.北京科技大学土木与资源工程学院，北京 100083； 2.金属矿山高效开采与安全教育部重点试验室(北京科技大学)，北京 100083)

如何将地表堆积固体废弃物安全高效的回填至充填采空区，是充填采矿法的核心工程问题之一.在长期的充填实践中，人们逐渐认识到水力管道输送模式与机械搬运模式相比，具有成本低、效率高和工艺简单等诸多经济技术优势[1-2].管道输送工艺也已成为充填采矿技术的最重要的核心工艺环节[3].

管道输送工艺应用过程中，由于充填原材料物化性质的不同，所制备的充填浆体有不同的流变特性和流变参数.具备不同流变特性和流变参数的充填浆体在不同压差、管径和管路长度条件下的充填流量和沿程阻力计算结果，是充填配比优化、充填系统设计、泵送设备选型和充填管道选型的必要依据.只有准确确定充填料浆的流变参数、利用准确的流动方程才能准确计算不同工况下充填料浆的管道流量和沿程阻力，保证充填系统的安全运行.

目前关于充填料浆管道输送基础理论的研究严重滞后于技术发展[4].充填料浆一般为固液两相流，在固体质量分数低时表现为非均质流体，固体质量分数高时表现为均质流体[5-6].对于低固体质量分数的非均质充填料浆，在管道流动时绝大多数情况下为湍流，只能通过大量的试验得到经验方程，无法获取准确的数学模型[7].随着膏体充填技术的发展和推广，矿山企业逐渐意识到充填固体质量分数的提高对降低充填材料成本和保护环境的显著作用，充填浆体也由非均质流体向均质膏体转变，均质充填浆体管道输送的问题越来越多被广大矿山企业和研究单位所关注，其中以牛顿流体和非牛顿流体中的宾汉体作为本构模型的充填浆体管道流动方程应用最为广泛，上述两个本构模型所推导的流动方程分别为Poiseuille定律和Buckingham流动方程[8].由于均质充填料浆工况的复杂性[9-10]，牛顿流体和宾汉体本构模型通常并不精确[11-12]，而较为精确的幂率模型所得出的流动方程形式复杂、计算量大，不易进行工程应用，而盲目采用牛顿流体和宾汉体本构模型所得出的流动方程计算结果与现场实际测量结果常常存在较大误差，进而导致在充填过程中常常发生堵管事故，影响井下采矿效率[13].

针对上述问题，本文力图推导出具有普适性的非牛顿流体管道流动方程，使之适用于所有流型的非牛顿流体，从而能够有效指导充填料浆管道输送工程实践.为了便于现场工作人员应用，本文进一步推导出非牛顿流体沿程阻力和管道流量方程的近似表达式，大幅度的降低了计算量，简化了计算过程.

1 核心内容及研究路线

本文核心内容是根据流变学以及流体力学理论，推导非牛顿流体管道流动方程与近似管道流量方程，建立非牛顿流体充填料浆流变参数与管道压差、管径、管路长度及充填流量的函数关系，并通过管道试验验证近似管道流量方程的准确性.具体研究路线如图1所示.

图1 本文的研究路线

2 非牛顿流体管道层流流动方程

2.1 流速分布方程

充填料浆的本构模型采用可涵盖绝大多数常见非牛顿流体的幂率模型[14-21]，即

式中：τ为剪应力；τB为流变参数中的屈服应力；K为流变参数中的稠度系数；n为流变参数中的幂率指数；du/dr为速度梯度.

由幂率模型的定义可得：

(1)

流体在一定压差作用下在管道内的剪应力分布如图2所示.

图2 流体在圆管内的剪应力分布

对管壁内半径为R的流体部分，进行静力平衡分析可得：

ΔP·πR2-2πRL·τw=0.

式中：ΔP为管道两端压差；R为管道半径；L为管路长度；τw为管壁四周的剪切应力.

同理，对任意半径为r的流体部分，进行静力平衡分析可得:

ΔP·πr2-2πrL·τ=0.

(2)

式中：ΔP为圆柱体两端压差；r为圆柱体半径；L为圆柱体长度；τ为圆柱体的剪切应力.

将式(2)代入式(1)可得:

(3)

对式(3)两边分离变量，并进行积分可得:

(4)

式中C为常数.根据边界条件，r=R时，u=0，可求得常数C的数值为

将常数C代入式(4)可得流体在管道内的任一流层的速度分布u为

当考虑τB存在情况时，流体在管道内的速度分布如图3所示.

图3 流体在圆管内的层流流动

当τB不为0时，管道中柱塞流区和剪切流区同时存在，柱塞流区和剪切流区的特征为：在柱塞流区r0内，各流层的速度相等且最大；在剪切流区各流层的速度不等，且越靠近管壁速度越小.

2.2 管道流量方程

管道流量即管道流速与管道截面积的乘积，因此管道流量可用管道内流体速度分布方程进行积分求得，即

积分后可得流量Q的表达式为

(5)

式(5)即非牛顿流体管道流动方程.可见只要能够准确确定充填料浆的流变参数，就通过该方程能准确计算充填浆体在不同压差、管径和管路长度条件下的充填流量.反之，如果只要能够通过简单的管道试验(L管、倾斜管或环管等)准确测定充填浆体在不同压差、管径或管路长度条件下的充填流量，也就可以通过以该方程所构建的非牛顿流体管道流动方程组准确的反向求解流变参数.

当n=1时，τB不等于0时，为Buckingham流体，此时式(5)可简化为

(6)

浆体流速v的计算公式为

(7)

联立式(2)、(6)、(7)可得Buckingham流动方程:

当n≠1时，τB=0时，为假塑性流体或胀塑性流体，此时式(5)可进一步简化为

当n=1时，τB=0时，为牛顿流体，此时式(5)可进一步简化为Poiseuille定律:

2.3 过渡流速方程

非牛顿流体管道流动方程仅适用于管道层流，对于非层流状态的管道流动，本方程并不适用.为了对管道流动状态进行判别，需要引入雷诺数并在此基础上推导出过渡流速方程，明确非牛顿流体管道流动方程的适用条件，从而实现料浆流态的准确判别和管道流动方程的正确应用.

雷诺数Re计算公式为

(8)

式中：D(2R)为管道直径；η为料浆表观黏度；v为料浆平均流速；ρm为料浆比重.

根据式(8)可见，表观黏度越小，雷诺数越大.表观黏度是指在一定速度梯度下，用相应的剪切应力除以剪切速率所得的商.表观黏度η的定义为

(9)

联立式(3)、(9)可得非牛顿流体在管道内的表观黏度分布为

(10)

由式(10)可知r=R时，η有最小值，由于η与Re负相关，此时Re有最大值.即管壁处的雷诺数满足层流状态时，管道内任意位置也满足层流状态.

为求解管壁处的雷诺数，取r=R，同时联立式(7)、(8)、(10)可得非牛顿流体管道流动雷诺数Re的最大值为

(11)

方程(11)即非牛顿流体过渡流速方程，当非牛顿流体在管道内流动时，根据方程(11)计算得雷诺数Re的最大值小于2 300时，充填浆体在管道内才全部呈现层流状态，此时本文推导的非牛顿流体管道流动方程才能成立.

2.4 近似管道流量方程

由于本文推导的非牛顿流体管道流量方程形式复杂、计算量极大，为了便于现场工作人员应用，需要对其进行简化以减少计算量、简化计算过程.

在矿山常见工况条件下(0 Pa<τB<100 Pa、3.6 MPa<ΔP<7.2 MPa、900 m

由图4可见，均质充填浆体在常见管道输送工况条件下，其实际壁面剪切速率与管道平均流速正相关、与管道直径负相关，且du/dr≈8v/D.由于管壁处du/dr≈8v/D，即

(12)

将8v/D替换为4v/R，联立式(7)、(12)并化简后可得：

(13)

式(13)即非牛顿流体近似管道流量方程.通过该方程能近似计算已知流变参数的充填浆体在不同压差、管径和管路长度条件下的充填流量.反之，也可以通过该方程所构建的近似管道流动方程组反向近似求解流变参数.

图4 常见工况du/dr与8v/D数值对比

3 试验

3.1 试验材料

试验中采用矿山全尾砂、325普通硅酸盐水泥混合制备非牛顿流体充填料浆.矿山全尾砂取自某铜矿山尾矿库，体积质量为2 670 kg/m3，容重为1 392 kg/m3，孔隙率为47.87%.图5为该矿山铜尾砂及水泥的粒径分布图，全尾砂中-74 μm筛下颗粒总量占比为70.85%，-37 μm占比46.96%，水泥中的细颗粒含量占比尾砂更大，可以看出该矿山尾砂及水泥粒径较细，适于制备均质充填料浆.试验材料的化学成分分析可参考文献[22]，此处不再赘述.

图5 水泥和尾砂的粒径分布

按照11%的水泥掺量与水、尾砂混合，强力搅拌制备为充填料浆，所制备的充填料浆基本物理性质见表1.

表1 充填料浆基本物理性质

3.2 试验装置

试验系统是由深锥浓密机、可储存水泥的微粉称、一级搅拌机、二级搅拌机、膏体柱塞泵、充填管道、压力传感器和上位机构成的环管实验装置，如图6所示.

该环管装置设置的充填管道内径为Φ150 mm.为监测管内压力，共设置压力表2块，分别标记为P1、P2号压力表，P1与P2号压力表之间管道长度为26 m，为监测管道流量，在柱塞泵与管道接口法兰处布设OPTIFLUX 4300C型流量计.

图6 环管试验装置系统布置

3.3 试验原理

由于本文试验采用固定式环管试验装置，其管径R、充填管路长度L不变，通过膏体泵改变不同的管道两端压差ΔP1、ΔP2和ΔP3，测出口断面流量Q1、Q2、Q3，代入下式(14)，通过解方程组，即可求得料浆的屈服应力τB、稠度系数K和幂率指数n.将试验参数代入下式得到充填料浆流变参数的数值解以后，将屈服应力τB、稠度系数K和幂率指数n的数值代入式(11)判别其是否为层流，当Re<2 300时，为层流，此时流动方程成立.

以近似解方法求解流变参数近似解的原理与解析解一致，对于本文试验采用的试验装置，其管径R、充填管路长度L不变，通过膏体泵改变不同的管道两端压差ΔP1、ΔP2和ΔP3，测出口断面流量Q1、Q2、Q3，代入下式，通过解方程组，即可求得料浆的3个流变参数即屈服应力τB、稠度系数K和幂率指数n的近似解，将近似解方法所求得的流变参数数值与解析解所得流变参数数值对比，考察其差异.

(15)