0 引 言
巷道局部阻力类型主要为多条巷道岔路区域、巷道拐弯区域、巷道断面改变区域，其中巷道转弯区域是井下最常见的局部阻力区域，研究巷道转弯区域局部阻力损失对于进行新掘巷道局部通风阻力预测工作，以及为矿井风网解算提供至关重要的准确的风网风阻参数具有重要意义。不同于巷道摩擦阻力，巷道转弯区域局部阻力产生原因是风流场突然发生急剧变化，出现复杂的耗能漩涡。虽然前人建立了巷道转弯区域局部阻力计算公式，但公式中作为关键参数的局部阻力系数依靠经验取值，无法保证阻力计算结果的准确性。
巷道通风阻力损失与巷道风流场密切相关，从巷道转弯区域风流场分布规律研究出发，揭示巷道转弯区域局部阻力产生原因及局部阻力损失规律,从而能够准确计算巷道转弯区域局部通风阻力。现阶段计算机数值模拟方法的兴起为风流场精准描述提供了技术支持，目前越来越多的研究人员开始采用数值模拟方法对复杂结构空间内风流场进行研究，张浪[1]借助CFD方法模拟测风站巷道风流场结构，对巷道测风站风速传感器平均风速测定位置优化；李艳昌等[2]模拟研究了人员存在条件下巷道风速场分布规律；张京兆等[3]通过数值模拟方法研究了矩形转弯巷道风速场分布规律；支学艺等[4]通过数值模拟方法研究了风流经转弯后风速场分布规律；张宏杰等[5]通过数值模拟方法研究了防爆胶轮车行驶过程巷道风流场分布规律；尹斌等[6]通过数值模拟方法研究了通道式巷道调节风窗局部风阻计算公式模型；文献[7-8]通过巷道断面风流场模拟确定巷道断面平均风速监测位置区域，李曼等[9]借助CFD模拟方法对比分析巷道中有无测风支架条件下巷道风速场分布规律，根据模拟结果改进巷道测风方法。众多研究人员借助数值模拟方法对井巷风流场、通风构筑物所在巷道风流场进行定量分析研究，揭示通风阻力损失，并取得了一定成果。
在前人研究基础上，笔者采用数值模拟方法对巷道转弯区域这类最常见巷道的局部阻力进行深入研究，以山西天地王坡煤矿南集中回风巷与上寺头回风巷的巷道转弯区域为研究对象，利用FLUENT流体动力学数值模拟软件对一定巷道尺寸及弯角结构条件下巷道转弯区域风流场进行模拟分析，通过对比现场实测结果、CFD模拟计算结果、经验公式计算结果，考察CFD模拟方法的可靠性，在此基础上，推导建立以巷道风量、转弯角度为自变量的南集中回风巷与上寺头回风巷转弯区域局部阻力系数计算公式。
1 巷道转弯区域风流场CFD计算模型
1.1 几何模型建立
山西天地王坡煤矿南集中回风巷与上寺头回风巷距离较长，且巷道沿程存在多个转弯区域，南集中回风巷转弯区域采用直线弯角结构，而上寺头回风巷转弯区域采用曲线弯角结构，以南集中回风巷与上寺头回风巷为研究对象，利用GAMBIT几何建模软件分别构建相应的几何模型，并实测巷道未转弯区域摩擦阻力系数，几何参数及摩擦阻力系数取值见表1，图1为转弯角度为90°的南集中回风巷与上寺头回风巷转弯区域几何模型。
表1 南集中回风巷与上寺头回风巷几何建模参数
Table 1 Geometric modeling parameters of South Centralized return lane and Shangsitou return lane
巷道参数南集中回风巷上寺头回风巷断面形状矩形矩形支护类型锚喷锚喷长度/m2 3411 572高度/m3.213.08宽度/m4.293.87弯角结构直线弯角结构曲线弯角结构断面积/m213.7711.92周长/m15.0013.90摩擦阻力系数/(N·s2·m-4)0.007 7290.010 857
图1 巷道转弯区域几何模型
Fig.1 Geometric model of turning area of the lane
1.2 数学模型建立
巷道内空气流动过程采用N-S方程组求解，满足质量守恒、动量守恒、能量守恒，巷道流场方程组见式(1)，巷道内风流为紊流状态，描述巷道风流场需要在N-S方程组基础上补充湍流模型方程组，由于巷道转弯区域风流场会发生急剧变化，REG k-ε模型可以很好地处理高应变率及流线弯曲程度大的流动，N-S方程组与REG k-ε模型联立就是巷道流场计算模型的微分方程组，在此基础上再加上边界条件，就构成了完整的可求解巷道风流场的CFD计算模型，计算模型边界条件如下：
巷道壁面Wall(壁面)巷道入口Velocity-in(风速入口)巷道出口Flow-Out(风速出口)
(1)
式中:ρf为风流密度，kg/m3；x为计算区域某一位置x轴方向坐标值，m；y为计算区域某一位置y轴方向坐标值，m；z为计算区域某一位置z轴方向坐标值，m；；u为计算区域内某一位置x轴方向风速，m/s；v为计算区域内某一位置y轴方向风速，m/s；w为计算区域内某一位置z轴方向风速，m/s；μ为空气动力黏度，取值为17.9×10-6 Pa·s；p为气体压力，Pa；Su为x轴方向的动量源项，kg/(m2·s2)；Sv为y轴方向的动量源项，kg/(m2·s2)；Sw为z轴方向的动量源项，kg/(m2·s2)；T为气体温度，K；k为空气传热系数，取值为0.023 W/(m·K)；CP为空气比热容，取值为717 J/(kg·K)；ST为黏性耗散项，kg·K/(m3·s)。
1.3 巷道转弯区域流场分布模拟结果分析
1.3.1 巷道转弯区域流场分布
图2、图3分别为转弯角度为45°与90°、风量为48 m3/s条件下上寺头回风巷与南集中回风巷转弯区域静压场分布模拟结果，由图2、图3可得：
1)当巷道弯角较小时，低静压区域集中于弯角拐弯区域，上寺头回风巷与南集中回风巷转弯区域静压分布无明显差异。
2)当巷道弯角较大时，低静压区域集中于巷道转弯后方区域，巷道转弯角度对巷道风流场影响巨大，上寺头回风巷与南集中回风巷转弯区域静压分布存在明显差异。
图2 巷道转弯区域静压场数值模拟结果(转弯角度45°)
Fig.2 Numerical simulation results of static pressure field distribution in corner area of roadway (turning angle was 45°)
图3 巷道转弯区域静压场数值模拟结果(转弯角度90°)
Fig.3 Numerical simulation results of static pressure field distribution in corner area of roadway (The turning angle was 90°)
1.3.2 巷道转弯区域沿程气压及阻力损失变化规律
图4为南集中回风巷转弯区域沿程气压变化及阻力损失变化曲线，由图4可得：
1)静压变化曲线呈现“平缓减小—大幅度下降—小幅度增大—平缓减小”趋势变化，动压变化曲线呈现“平稳—大幅度增大—大幅度减小-平稳”趋势变化，全压变化曲线呈现“缓慢减小—大幅度减小-平缓减小—小幅度减小-平缓减小”趋势变化，静压、动压、全压在巷道转弯处发生同步剧烈变化。
2)巷道沿程阻力损失出现2个波峰变化，首个波峰出现在静压、动压、全压大幅度骤变位置，第2个波峰出现在静压、动压、全压小幅度骤变位置，且首个波峰明显大于第2个波峰，巷道转弯区域阻力损失主要为首个波峰区域，巷道转弯区域通风阻力损失具有明显的突变特性，属于典型的局部阻力损失。
图4 巷道转弯区域沿程气压变化及通风阻力损失变化曲线(巷道转弯角度90°、巷道风量48 m3/s)
Fig.4 Variation curves of on-way ventilation resistance loss and pressure of turning area of roadway (The turning angle was 90°,roadway air volume was 48 m3/s)
2 巷道转弯区域通风阻力CFD计算方法误差分析
通过对比巷道转弯区域通风阻力现场实测结果、经验计算公式计算结果、CFD模拟计算结果的差异，考察巷道转弯区域风流场CFD计算方法的可靠性以及经验计算公式的计算误差。分别采用经验计算公式方法、现场实测方法和CFD模拟计算方法对王坡煤矿南集中回风巷3处巷道转弯区域和上寺头回风巷4处巷道转弯区域通风阻力进行测试或求解。
采用巷道转弯区域通风阻力经验计算公式计算南集中回风巷和上寺头回风巷转弯区域通风阻力[10]，见式(4)，计算中各参数取值见表2。
(2)
式中：hf为巷转弯区域通风阻力，Pa；ξ0为巷道壁面完全光滑、巷道弯角为90°条件下的局部阻力系数，该参数取决于转弯处内角曲率半径与巷道宽度之比，(N·s2)/m4；α为巷道摩擦阻力系数，(N·s2)/m4；H为巷道高度，m；B为巷道宽度，m；β为巷道弯角影响系数；ρ为风流密度，kg/m3；S为巷道转弯区域风流入口位置断面面积，m2；Q为巷道风量，m3/s。
巷道转弯区域通风阻力现场实测方法实施步骤如下：
1)利用机械式风表测试巷道风量，利用干湿温度计实测巷道风流温度、湿度，计算巷道风流密度，测试数据见表2。
2)利用倾斜压差计与皮托管配合实测巷道转弯区域风流入口位置与出口位置之间气流总压差，计算得到巷道转弯区域通风阻力，计算结果见表3。
表2 巷道转弯区域实际尺寸、实测阻力特性参数及实测风流参数汇总
Table 2 Summary of laneway size parameters, resistance characteristics parameters and air flow parameters in turning area of roadway
巷道转弯区域转弯角度/(°)实测风量/(m3·s-1)巷道宽/m巷道高/m实测温度/℃转弯处内角曲率半径/m实测湿度/%ρ/(kg·m-3)巷道弯角影响系数ξ0/(N·s2·m-4)南集中回风巷12059.334.293.2119.744.29821.0060.120.33南集中回风巷24559.334.293.2119.934.29830.9960.330.33南集中回风巷38559.334.293.2119.864.29821.0080.900.33上寺头回风巷13078.283.873.0819.713.87861.0110.190.33上寺头回风巷25278.283.873.0819.263.87860.9960.410.33上寺头回风巷37078.283.873.0819.563.87850.9910.630.33上寺头回风巷49078.283.873.0819.813.87850.9721.000.33
巷道转弯区域通风阻力CFD计算方法实施步骤如下：
1)以实际巷道尺寸、实测巷道风量、实测巷道转弯角度进行CFD几何建模与物理建模，构建CFD计算模型，求解巷道转弯区域风流场，建模参数见表2。
2)根据CFD风流场模拟结果计算巷道转弯区域通风阻力，巷道转弯区域局部阻力为巷道转弯区域风流入口位置气流总压减去巷道转弯区域出口位置气流总压，见式(3)，CFD方法计算结果见表3。
(3)
式中：n、m分别为巷道转弯区域风流入口位置和出口位置巷道断面划分单元数量；Vi、Vj分别为巷道转弯区域风流入口位置断面i单元与出口断面j单元处风速，m/s；Pi、Pj分别为巷道转弯区域风流入口断面i单元与巷道出口断面j单元处气压，Pa；hi、hj分别为巷道转弯区域风流入口断面i单元与巷道出口断面j单元处标高，m；g为重力加速度，m/s2。
南集中回风巷和上寺头回风巷转弯区域通风阻力现场实测结果与经验计算公式计算结果之间相对误差在60%以上，而现场实测结果与CFD模拟计算结果二者之间相对误差小于5%，表明以实际巷道尺寸、风量、密度参数为建模参数的巷道转弯区域风流场CFD计算方法能够用于准确计算巷道转弯区域风流场，而经验计算公式无法满足计算要求。
表3 巷道转弯区域通风阻力现场实测结果、经验公式计算结果、CFD计算结果及相对误差汇总
Table 3 Summary of field measurement results, empirical formula calculation results, CFD calculation results and relative errors of ventilation resistance in turning area of roadway
巷道转弯角度/(°)实测风阻/Pa经验公式计算风阻/PaCFD计算风阻/Pa经验公式计算风阻与实测风阻相对误差/%CFD计算风阻与实测风阻相对误差/%南集中回风巷205.000.734.7685.364.80南集中回风巷458.191.998.5375.674.15南集中回风巷8520.455.5019.4573.104.89上寺头回风巷3010.123.039.8770.082.47上寺头回风巷5216.726.4417.5461.504.90上寺头回风巷7027.669.8426.8764.422.86上寺头回风巷9045.6215.3247.6266.414.80
3 风量与转弯角度对巷道转弯区域局部通风阻力的影响规律定量分析
在分析了巷道转弯区域通风阻力损失特征以及论证了CFD模拟方法可靠性的基础上，笔者以南集中回风巷与上寺头回风巷为研究对象，因研究对象的巷道尺寸及巷道弯角结构已确定，采用CFD模拟方法重点研究巷道风量与弯角角度这2个主要因素对巷道转弯区域通风阻力的影响规律。
图5 巷道转弯区域通风阻力随巷道转弯角度变化曲线
Fig.5 Change curves of ventilation resistance in laneway turning area with laneway turning angle
巷道转弯区域通风阻力CFD计算模拟中巷道风量分别取12、24、36、48、60、72、84、96 m3/s八种情况，巷道转弯角度分别取15°、30°、45°、60°、75°、90°六种情况，南集中回风巷与上寺头回风巷转弯区域通风阻力随转弯角度的变化曲线如图5所示，南集中回风巷与上寺头回风巷转弯区域通风阻力随巷道风量的变化曲线如图6所示，由图5、图6可得，巷道转弯区域通风阻力随巷道转弯角度增大呈现“缓慢增大—迅速增大—缓慢增大”S型曲线趋势增大；巷道转弯区域通风阻力随巷道风量增大以“斜率逐渐增大”的趋势增大，相比于巷道风量因素，转弯角度对巷道转弯区域通风阻力的影响表现得更为复杂。
图6 巷道转弯区域通风阻力随巷道风量变化曲线
Fig.6 Change curves of ventilation resistance changing with roadway air volume in turning area of roadway
4 基于CFD模拟方法的巷道转弯区域局部阻力系数计算公式推导
巷道转弯区域通风阻力属于局部阻力范畴，以风阻为定值的通风阻力定律不能直接应用于巷道转弯区域通风阻力计算，局部通风阻力按式(4)计算，式中的局部阻力系数为关键参数，该参数直接决定局部阻力计算结果准确性，而以往的局部阻力经验计算公式的计算结果与实际情况误差偏大。
为获得准确的巷道转弯区域局部通风阻力系数计算公式，笔者以巷道转弯区域风流场研究为基础，根据CFD模拟结果推导建立以巷道风量、转弯角度为自变量的南集中回风巷与上寺头回风巷转弯区域局部通风阻力系数计算公式
(4)
式中：ξ为巷道转弯区域局部阻力系数，(N·s2)/m4。
4.1 计算公式中引入巷道风量
首先定量考察巷道风量对巷道转弯区域局部阻力的影响程度，将巷道风量参数引入南集中回风巷与上寺头回风巷转弯区域局部通风阻力系数计算公式。南集中回风巷与上寺头回风巷转弯区域局部阻力随巷道风量变化拟合曲线如图7所示，拟合结果表明：巷道转弯区域局部阻力与巷道风量之间满足幂函数关系，初步构建了以巷道风量为自变量的巷道转弯区域局部阻力计算公式，见式(5)，拟合参数见表4，将式(5)转化为式(6)形式。联立式(4)与式(6)，初步得到以巷道风量为自变量的巷道转弯区域局部风阻系数计算公式，见式(7)。
图7 巷道转弯区域局部阻力随巷道风量变化拟合曲线
Fig.7 Fitting curves of local resistance in laneway turning area with laneway air volume change
hf=aQb
(5)
hf=(aQb-2)Q2
(6)
(7)
式中：a、b均为拟合参数，取值与巷道转弯角度、巷道转弯结构、巷道尺寸等参数相关，其中b参数取值范围为[1,2]，笔直巷道内风流处于层流状态条件下b取值为1，笔直巷道内风流处于紊流状态条件下b取值为2。
表4 巷道转弯区域通风阻力随风量变化拟合曲线参数
Table 4 Parameters of fitting curves of ventilation resistance varying with air volume in turning area of roadway
巷道转弯区域转弯角度/(°)abR2上寺头回风巷1150.005 21.674 70.999 8上寺头回风巷2300.005 61.708 40.999 4上寺头回风巷3450.006 21.756 60.999 4上寺头回风巷4600.007 41.803 10.998 5上寺头回风巷5750.009 11.843 40.998 7上寺头回风巷6900.009 91.956 20.998 1南集中回风巷1150.005 41.669 50.999 7南集中回风巷2300.005 51.689 80.999 9南集中回风巷3450.005 61.750 10.999 9南集中回风巷4600.007 31.808 10.999 1南集中回风巷5750.008 21.863 80.997 7南集中回风巷6900.010 11.870 70.998 1
4.2 计算公式中引入巷道转弯角度
在考虑巷道风量对巷道转弯区域局部阻力影响的基础上，定量考察巷道转弯角度对巷道转弯区域局部阻力的影响程度，将巷道弯角角度引入南集中回风巷与上寺头回风巷转弯区域局部通风阻力系数计算公式，对式(7)中a值、b值与巷道转弯角度之间关系进行非线性拟合，拟合曲线如图8、图9所示。
图8 上寺头回风巷局部阻力系数计算公式a、b值随巷道风量变化拟合曲线
Fig.8 Fitting curves of a and b values in the calculation formula of local resistance coefficient of Shangsitou air return roadway with changes in roadway air volume
图9 南集中回风巷局部阻力系数计算公式a、b值随巷道风量变化拟合曲线
Fig.9 Fitting curves of values a and b in the calculation formula of local resistance coefficient of South Central air return roadway with changes in roadway air volume
a值、b值随转弯角度均呈指数函数形式变化，见式(8)、式(9)，q、r、s、p为拟合参数，具体拟合参数见表5，将式(8)、式(9)代入式(7)，得到以巷道转弯角度和巷道风量为自变量的巷道转弯区域局部风阻系数计算公式，见式(10)。
a=qerθ
(8)
b=sepθ
(9)
式中：θ为巷道转弯角度，(°)；q、r、s、p均为拟合参数，参数取值与巷道转弯结构、尺寸相关。
(10)
表5 局部阻力系数计算公式q、r、s、p参数取值
Table 5 Values of q,r,s and p parameters in the calculation formula of local resistance coefficient
巷道qrsp上寺头回风巷0.004 30.009 21.612 80.002 0南集中回风巷0.004 40.008 61.615 90.001 8
式(10)是在充分研究了风量与转弯角度对南集中回风巷与上寺头回风巷转弯区域局部通风阻力影响规律的基础上所得的局部通风阻力系数计算公式，该式能够预测计算任意风量与转弯角度条件下南集中回风巷与上寺头回风巷巷道转弯局部通风阻力，由于该公式是在巷道尺寸及弯角结构已确定的前提下所得，公式中q、r、s、p参数取值与巷道转弯结构、巷道断面尺寸之间的关系尚未明确，需进一步深入研究。
5 结 论
1)以山西天地王坡煤矿南集中回风巷与上寺头回风巷转弯区域为研究对象，构建巷道转弯区域风流场CFD计算模型，分析巷道转弯区域气压场分布特征，分析沿程静压变化规律、动压变化规律、全压变化规律、通风阻力损失变化规律，研究表明巷道转弯区域具有“前大波峰-后小波峰”的双波峰通风阻力损失特征。
2)以山西天地王坡煤矿南集中回风巷与上寺头回风巷的巷道转弯区域为研究对象，通过对比通风阻力经验计算公式计算结果、现场实测结果、CFD模拟计算结果，验证了CFD模拟巷道转弯区域风流场的可靠性。在此基础上，考察了巷道风量与巷道转弯角度对巷道转弯区域局部通风阻力损失的影响规律，巷道转弯区域局部通风阻力随转弯角度增大呈S型增大变化趋势，巷道转弯区域局部通风阻力随风量呈现以斜率逐渐增大的增大趋势变化。
3)针对目前巷道转弯区域通风阻力经验计算公式计算结果与实际偏差大的问题，在CFD模拟计算方法的基础上，采用数值分析方法构建了以巷道风量、巷道转弯角度为自变量，南集中回风巷和上寺头回风巷巷道转弯区域局部通风阻力系数计算公式，该计算公式能够预测计算南集中回风巷与上寺头回风巷在任意风量与转弯角度条件下巷道转弯区域的局部通风阻力。尽管该研究所得局部阻力系数计算公式仅能够用于计算南集中回风巷和上寺头回风巷巷道转弯区域局部阻力，但笔者论证了采用“CFD模拟分析+数值分析函数曲线拟合”方式能够建立不同尺寸结构巷道转弯区域局部阻力系数计算公式，且计算结果具有较高的可靠性，该研究可用于新建矿井或新掘弯曲巷道的通风阻力预测工作和矿井通风阻力精细化测试工作，同时也能为矿井风网解算提供准确的风网风阻参数。
[1] 张 浪.巷道测风站风速传感器平均风速测定位置优化研究[J].煤炭科学技术,2018,46(3):96-102.
ZHANG Lang. Optimized study on location to measure average air velocity with air velocity sensor in wind measuring station of underground mine[J]. Coal Science and Technology,2018,46(3):96-102.
[2] 李艳昌,张 箭,杨雨濛,等.人员对巷道风速测量影响的数值模拟分析[J].金属矿山，2016(1):167-170.
LI Yanchang, ZHANG Jian, YANG Yumeng，et al. Numerical simulation analysis of pedestrian’s influence on airflow velocity measuring in roadway[J].Metal Mine, 2016 (1):167-170.
[3] 张京兆,张红记,姜红敏,等.矩形巷道转弯后测风站位置的数值分析[J].中国矿业，2015,24(1):142-145.
ZHANG Jingzhao, ZHANG Hongji, JIANG Hongmin，et al. Numerical analysis of air measuring station location in a rectangular air way after turning[J].China Mining Magazine, 2015, 24(1):142-145.
[4] 支学艺,刘建华,赖春明.风流经转弯后稳定段长度数值模拟及应用研究[J].有色金属科学与工程,2016,7(1):96-101.
ZHI Xueyi, LIU Jianhua, LAI Chunming. Numerical simulation and application of airflow’s stable length after turning[J].Nonferrous Metals Science and Engineering, 2016,7(1):96-101.
[5] 张宏杰,刘彦青,姚鹏飞,等.井下防爆胶轮车运行对矿井通风阻力影响的定量分析[J].煤矿安全,2017,48(7):184-187.
ZHANG Hongjie, LIU Yanqing, YAO Pengfei，et al. Quantitative research on the effect of running processes of explosion-proof rubber-tyred car on the resistance of mine ventilation system[J]. Safety in Coal Mines,2017,48(7):184-187.
[6] 尹 斌,刘彦青,李 伟,等.矿用自动控制风窗局部风阻数值模拟[J]. 煤矿安全,2016,47(10):172-175.
YIN Bin，LIU Yanqing， LI Wei，et al. Numerical simulation for local wind resistance of mine automatic control wind window[J]. Safety in Coal Mines,2016, 47(10): 172-175.
[7] 鹿广利,武赞龙,赵剑锋.不同拐弯角度下巷道内风流变化规律的数值模拟[J].矿业研究与开发,2019,39(12):116-121.
LU Guangli,WU Zanlong,ZHAO Jianfeng.Numerical simulation on the change law of air flow in roadway with different turning angles [J].Mining R & D,2019,39(12):116-121.
[8] 宋 莹,刘 剑,李雪冰,等.矿井巷道风流平均风速分布规律的试验与模拟研究[J]. 中国安全科学学报,2016,26(6):146-151.
SONG Ying, LIU Jian, LI Xuebing，et al.Experiment and numerical simulation of average wind speed distribution law of airflow in mine tunnel[J]. China Safety Science Journal,2016,26(6): 146-151.
[9] 李 曼,马 欢.矿井主通风机风量测试方法的模拟研究[J].煤炭科学技术,2017,45(2):151-155.
LI Man, MA Huan. Simulation research on airflow measuring method of mine main ventilator[J]. Coal Science and Technology,2017,45(2):151-155.
[10] 张国枢.通风安全学[M].徐州:中国矿业大学出版社,2007.
