百叶式风窗流场分布及局部阻力特征研究*

王 宁，李雨成，张 欢，张 静，李博伦

(太原理工大学 安全与应急管理工程学院，山西 太原 030024)

0 引言

通风系统是矿山地下工程建设与优化的子系统之一，而风量的精准调控是实现通风优化这一目标的基石[1-5]。风量精准调控的实质是对有效过风面积的调节。与无构筑物平直巷道不同，当风窗存在时，通风系统中的局部阻力系数会明显增加[6-7]。之前对风窗两端的风阻多以挡板式调节风窗作为研究对象，随着百叶式风窗的出现，传统风窗的局部阻力计算公式不再适用[8]。百叶式风窗局部阻力与流场分布特征密切相关，为此，有必要对百叶式风窗局部阻力及流场分布特征进行深入研究。

许多学者曾采用数值模拟对复杂结构空间内风流场分布特征进行研究[9]：1995年，高建良等[10]对梯形巷道两扇风窗的局部阻力进行了分析，发现理论公式中未考虑断面形状、风窗开口形状和位置、断面收缩变化等因素；2003年胡创义[11]对调节风窗面积的计算法进行研究,提出间接计算风窗面积的计算公式；2008年张迎新等[12]从理论上分析了矿井风量调节中调节风窗的作用机理；2014年黄光球等[13]分析井下风窗安装前后通风系统中风流流动特征，给出了度量风量和阻力反应敏感性强弱的方法；2015年邵和等[14]分别采用局部阻力测定实验系统和数值模拟的方法研究了不同雷诺数对风窗局部阻力系数的影响；2018年许克南等[15]基于压能守恒方程和风压特性曲线，推导出增设风窗前后压能变化的表达式；2018年刘永红等[16]利用粒子图像测速仪(PIV)对调节风窗模型流场进行实验测试，获得巷道流场纵向截面风流分布特征及调节风窗相关参数的关系式。虽然前人建立了传统风窗的局部通风阻力计算公式，但未考虑叶片角度参数与局部阻力的函数关系，因此，无法保证百叶式风窗局部阻力计算结果的准确性[17]。

本文运用CFD流体力学软件建立巷道-百叶式风窗-空气模型，对百叶式风窗区域流场分布特征进行数值模拟，并采用专业后处理软件Tecplot对计算结果进行绘制等值线图、迹线图等处理，分析风窗叶片开启不同角度下的局部阻力系数特征，准确计算出百叶式自动风窗区域局部通风阻力系数，进而实现风量的精准调控和矿井通风系统的优化。

1 数值模型建立与参数设置

1.1 百叶式风窗几何模型建立

为研究百叶式风窗前后风流运动过程、流场分布特征和局部阻力大小，利用CFD软件建立百叶式风窗的三维模型。如图1所示，箭头所指为风流流动方向，风窗所在巷道长100 m，风窗建立在Z=40 m的XY面上，窗框和门框厚度0.5 m，其他尺寸见表1。

表1 百叶式风窗几何模型关键参数Table 1 Key parameters in geometric model of louvered windshield

图1 百叶式风窗几何模型Fig.1 Geometric model of louvered windshield

过风面积直接影响百叶式风窗前后流场分布特征及局部阻力的大小。实验中通过变换叶片开启角度以改变过风面积，几何模型如图2所示。

图2 叶片不同开启角度示意Fig.2 Schematic diagram of blades with different opening angles

1.2 模型参数设置

为提高百叶式风窗流场模拟结果的准确性，采用四面体与六面体混合方式对物理模型进行网格划分，对百叶式风窗处的网格进行局部加密，从而降低网格质量对计算精度和稳定性的影响。结合矿井百叶式风窗的实际情况，在不影响模拟结果规律性的前提下，设置模型相关参数见表2。假设巷道壁面和百叶窗叶片表面光滑，空气为不可压缩流体，密度取1.225 kg·m-3；巷道中空气温度恒定为298 K；流场中的风流为稳态流动，动力黏度取1.789 4e-5kg/(m·s)；在基于压力的稳态求解器中，选择通过重整化群理论[18]的统计方法推导出来的两方程的湍流模型，采用SIMPLEC算法，对百叶式风窗流场分布特征进行数值模拟计算。

表2 百叶式风窗巷道风流场数值计算模型参数设置Table 2 Parameter settings of numerical calculation model for wind flow field of louvered windshield roadway

2 百叶式风窗流场空间分布特征

百叶式风窗叶片角度改变后，会使风窗局部流场的运动及分布特征发生急剧变化，这种变化包括风流汇集或分散的非稳定流动和高低风速接触处扰动形成贴壁涡流等。

2.1 风速场分布特征

以百叶式风窗叶片开启45°、入口风速2.5 m/s的工况，说明百叶式风窗速度变化规律。分别选取相邻叶片的中线做水平截面，得到Y=1.25,1.75,2.25,2.75 m高度下的风速分布特征。

图3 水平截面风窗前后速度Fig.3 Front and rear speed traces at horizontal sections of windshield

Y=1.25 m水平截面中，距离风窗入口1 m(Z=39 m)处的风速呈均匀分布，大小为4.5 m/s,风窗区域叶片开启45°后引起风窗入口风速突增到9 m/s，出口风速最大为15 m/s，随着与出口距离的增加，风速在减小；Y=1.75 m水平截面中，风窗入口风速为10.5 m/s，风窗出口风速为16.5 m/s，风窗区域及风窗出口1 m以内的风速均在12 m/s以上，形成主流的高风速区域，高风速区域面积扩大了近1倍；Y=2.25 m水平截面中，最高风速为16.5 m/s，随风流移动到了风窗后方0.5 m处，高风速区域面积延伸到了2.5 m，有向左右两侧扩展的趋势；Y=2.75 m水平截面中，高风速区域面积由风窗中心向四周扩散达到了最大。

整体来看，在百叶窗入口开始出现高风速流场，其流场区域面积和风速峰值随着底板高度逐渐增加，形成主流风速区；这种现象的产生原因与百叶窗叶片开启状态相关，叶片角度改变导致风窗区域局部阻力方向发生变化。

2.2 迹线特征

由图4所示，风流在百叶式风窗前的平直巷道中基本呈平滑直线分布，在风窗后方形成了增速减压区和减速增压区，这是由于风窗叶片向上开启45°后局部阻力突变，导致风流重新分配造成的。

图4 垂直截面风窗前后速度迹线Fig.4 Front and rear speed traces at vertical section of windshield

在增速减压区，紊流流体通过突变部位时，由于惯性力的作用，不能随从边壁突然转折，形成了主流风速区与边壁脱离的现象，在主流风速边界与边壁间形成上涡流区。产生的大尺度涡流不断被主风流带走，补充进去的流体又形成新的涡流，因而增加了局部阻力造成的能量损失。

在减速增压区，风窗后方巷道边壁也会产生涡流，流速沿程减小，静压不断增加，压差的作用与流动方向相反，使边壁附近本来很小的流速逐渐减小到零，在这里主流与边壁开始脱离，出现与主流方向相反的流动，形成下涡流区。

2.3 静压场特征

以百叶式风窗叶片开启45°，巷道入口风速2.5 m/s的工况，说明百叶式风窗前后压差变化规律。与风速场工况相同，分别选取风窗相邻叶片的中线做水平截面，得到Y=1.25，1.75，2.25，2.75 m高度下的静压分布特征，如图5所示。

图5 水平截面风窗前后局部放大压力云图Fig.5 Front and rear local enlarged pressure cloud diagrams at horizontal sections of windshield

Y=1.25 m水平截面中，百叶窗入口和出口静压为130，-20 Pa，相同0.5 m的距离内，风流在风窗前(Z=39.5～40 m)和风窗区域(Z=40～40.5 m)的压差为60，150 Pa；Y=1.75 m水平截面中，风窗进出口前后压差与Y=1.25 m截面基本相似，不同点在于风窗入口两侧靠近窗框位置的压差梯度变化更加明显；Y=2.25 m水平截面中，风流在风窗入口前0.5 m区域压差为30 Pa，风窗0.5 m区域内压差为180 Pa，风窗后方的负压区域面积由中心向两侧缩小明显；Y=2.75 m水平截面中，风窗入口前0.5 m区域压差为30 Pa，风窗区域压差不到150 Pa，风窗区域局部阻力在减小。

整体来看，风窗入口前的局部阻力小于风窗区域。随着距离底板高度的增加，风窗前0.5 m区域的局部阻力呈减小趋势，风窗区域局部阻力先增大后减小；风窗出口的负压区域面积不断减小，负压区域面积减小的方向是由中心向两侧移动。

3 百叶式风窗局部阻力系数及影响因素

3.1 百叶式风窗局部阻力特征

为了研究风窗局部阻力变化特征，对叶片开启角度为30°，45°，60°，90°，风速为1，2，2.5，3 m/s的不同工况条件进行数值模拟。根据色带值中颜色的差异，计算风流经过风窗前后的压差；根据公式(1)将压差值代入，得到风窗区域的局部阻力系数，不同角度下的风窗两端压差和局部阻力系数模拟计算结果见表3。

表3 百叶式风窗局部压差和阻力系数模拟计算结果汇总Table 3 Summary on simulated calculation results of local pressure difference and resistance coefficient of louvered windshield

(1)

式中：h为风窗两端压差，Pa；ρ为空气密度，kg·m-3；v为巷道风速，m·s-1；ξ为局部阻力系数。

由表3可知，巷道入口风速3 m/s，叶片开启30°，风窗两端压差最大为391.09 Pa；入口风速1 m/s，叶片开启90°时，风窗两端压差最小为20.95 Pa；整体来看，百叶窗局部阻力随风速的增大而增大，随着叶片开启角度的增大而减小；其原因是风流强度随叶片开启角度不同而急剧变化，并伴随有复杂的耗能漩涡出现，局部阻力的能量耗散主要和涡流区的存在相关。

3.2 局部阻力系数影响因素

百叶式风窗局部阻力与巷道流场分布特征密切相关，局部阻力系数可以作为局部阻力大小的评判标准，首先从风窗本身的特性规律进行研究，风窗有效过风面积随叶片开启角度的变化特征十分明显，根据其几何关系可得公式(2)。

S=5bl(1-cosθ)

(2)

式中：S为百叶式风窗有效开启面积，m2；b为单页百叶窗的宽度，m；l为单页百叶窗的长度，m；θ为百叶式风窗开启角度，(°)。

紊流状态百叶式风窗的局部阻力系数主要取决于阻力物的形状。局部风阻系数随叶片开启角度的增大而减小的规律极为明显，局部风阻系数的大小取决于有效过风面积大小。利用Origin对表4中所有相关数据进行非线性迭代拟合，如图6所示，以风窗过风面积为自变量的百叶式风窗局部阻力系数的计算见式(3)，拟合度达到0.985 2。

图6 百叶窗局部阻力系数与风窗开启面积关系Fig.6 Relationship between local resistance coefficient and opening area of louvered windshield

ξ=1.02S-0.369

(3)

将式(2)带入式(3)得到以风窗开启角度为自变量的百叶式风窗局部阻力系数计算见式(4)。

ξ=0.56[bl(1-cosθ)]-0.369(0<θ≤90°)

(4)

3.3 局部阻力系数计算公式的误差分析

为了验证百叶窗局部阻力系数计算方法的控制精度和可靠性，在井下宽4 m，高3 m的百米回风巷增设百叶窗，选取过风面积影响较稳定的风窗前后5～30 m处共布置6个传感器。为避免偶然误差，每个测试点进行6次测风，6次测风数据平均值即为巷道平均风速，根据局部阻力定律，计算得到风窗局部阻力系数，然后利用公式(4)计算相应叶片不同开启角度的局部阻力系数，对2种方法所得到的计算结果进行对比，实测与计算对比数据见表4，由表可知，实测结果与公式计算结果之间的相对误差均在5%以内，说明利用公式(4)计算得到的百叶窗局部阻力系数值可用于指导现场调风工作。

表4 百叶窗局部阻力系数实测与计算数据Table 4 Measured and calculated local resistance coefficients of louvered windshield

4 结论

1)主流风速区域在风窗入口处逐渐形成，其区域面积随巷道高度上升呈增大趋势。

2)百叶窗局部阻力与风速变化趋势一致，负压区域面积随底板高度上升呈减小趋势，压力梯度由中心向两端变小。

3)巷道边壁形成的涡流区是百叶窗局部阻力能量损失的主要影响因素。

4)百叶窗局部阻力系数取决于有效过风面积，与叶片开启角度成幂函数关系，其关系式为ξ=0.56[bl(1-cosθ)]-0.369(0<θ≤90°)。