0 引 言
粉尘污染是工业生产中的常见问题，长期反复接触粉尘会造成严重的肺部慢性病，在防护不当的情况下粉尘还是引起工业爆炸的原因之一。2017年4月国家环保部印发了《国家环境保护标准“十三五”发展规划》[1]，对粉尘污染物的排放提出了新的要求。而在当今传统的火力发电厂中，煤炭作为主要的发电能源，其输送系统通常由多个生产环节组成，包括卸料、运输、存放、取出、破碎、配煤至煤仓等[2]。这些生产环节在运行过程中，或在和其他生产环节交替的过程中会产生大量的煤尘，给电厂的安全生产带来巨大隐患，尤其在煤仓间等厂房工作期间[3-5]，对工作人员的身心健康造成严重影响，迫切需要对粉尘进行治理。
目前，应用最为广泛的除尘技术有袋式除尘[6]、电式除尘[7,8]和湿式除尘[9-11]。袋式除尘技术对煤种和粉尘特性的适应性强，对不同粒度的粉尘均有较高的捕集率，可同时脱除多种有害气体，但该技术除尘效率低、占地面积大、造价高、清灰复杂、滤袋寿命有限、检修维护工作量大、维护运行成本高。另外，针对输送系统的除尘来说，稳定有效地抽取粉尘气体是布袋除尘面临的主要挑战。静电除尘技术除尘效率高、阻力小、能耗低、运行费用低、方便管理、能处理高温烟尘，但该技术对外部条件变化敏感，高电阻粉尘和二次扬尘对电除尘器影响明显，有卸灰系统，产生二次污染，并且产生大量的检修维护工作量，增加维护运行成本。干雾除尘、多管冲式除尘等湿式除尘技术对无组织排放的封闭及半封闭环境中粉尘抑制作用明显，其中干雾抑尘技术效果优异，但其设备投资大，需要连续提供动力，耗电量大，用水量大，会造成二次污染，且检修维护工作量大。以上除尘技术主要用来处理粉尘气体收集充分的领域，如燃煤烟气粉尘处理等，对于厂房、车间等较大空间中的粉尘来说，起尘源多，气流压力分布不均匀，物料、气体、粉尘多相混合且相互间转化复杂，单一除尘技术无法完成有效处理。
无动力除尘技术[12]作为一种新技术，利用空气在流动时所产生势能转换成动能，使含尘空气进入多功能尘气分离装置后进行惯性、重力沉降、过滤等自动分离除尘，在一定程度上克服了以上3种传统除尘技术的不足。无动力除尘设备无需人工操作，检测控制准确，自动回旋彻底，除尘效果显著，并且投资少，无占地，免维护，无二次污染，使用安全可靠。但目前无动力除尘理论支撑有待提高，粉尘气体的运动规律有待进一步研究掌握，粉尘气体流动的合理性需要进一步加强。
因此，本文从粉尘运动的物理机理角度出发，通过建立粉尘运动的物理模型，采用数值模拟方法对粉尘的释放、扩散、运动、沉降等过程进行仿真计算，准确、全面掌握了厂房内粉尘污染物分布及其演化规律[13-16]，分析空气对流对粉尘污染物扩散及沉降的影响，为厂房内粉尘污染源布置、通风设计的模型和优化手段，以及除尘技术的改进提供参考。
1 物理模型
粉尘气体中通常包含气体、固态颗粒等物质，建立煤粉的释放模型、沉降模型和颗粒-气体两相流模型，在此基础上分析各相的参数指标，获得了多相流中各种变量的场分布规律。
1.1 粉尘释放模型
煤炭源释放的粉尘颗粒及粉尘黏附颗粒物是厂房内粉尘污染的主要原因，厂房内粉尘污染预测、评价和控制需要建模模拟粉尘颗粒物的释放过程。假设粉尘颗粒物的释放浓度各向同性，采用双指数衰减模型描述单位时间释放量：
E(t)=E01e-k11t+E02e-k12t
(1)
式中：E(t)为单位时间单位面积的释放量随时间变换的函数，mg/(h·m2)；E01为指数模型的快速初始释放参数；E02为双指数模型的慢速初始释放参数，mg/(h·m2)；k11为快速衰减常数；k12为慢速衰减常数，h-1，释放参数由实验实测数据拟合得到。
粉尘颗粒物的初始分布采用基于Rosin-Rammler粒径分布函数的喷射模型模拟：
(2)
式中：Yd为分布函数；d为粒径尺寸，
为中位半径，mm，是当颗粒物累计重量分布Yd=0.5处所对应的颗粒物粒径，此参数需通过实验数据获得，表征颗粒物粒径分布范围。同时还需要由实验数据或经验模型设置初始释放速度等条件。
1.2 粉尘沉降模型
粉尘沉降模型用于模拟悬浮微粒沉降。悬浮颗粒物从空气中降落回到皮带机或地面是一个复杂的耦合过程。在空间中运动时，粉尘颗粒物不仅会受到扩散和对流因素的影响，而且还会受到来自垂直方向上的重力(式(3))、浮力(式(4))和阻力(式(5))的作用。将颗粒物近似看成球形颗粒，其受力分析如图1所示。
重力
(3)
式中：dp为颗粒直径，m；ρp为颗粒密度，kg/m3；g为重力加速度，本文取9.8 m/s2。
浮力
(4)
式中：ρa为气体密度，kg/m3。
阻力
(5)
式中：Cξ为空气阻力系数，无量纲参量；u为颗粒与空气之间的相对速度，m/s。
图1 颗粒的空间垂直受力状态
Figure 1 Vertical stress state of particles in space
为求粉尘颗粒的沉降速度，除了考虑颗粒以及空气的物性影响外，还需要求出空气阻力系数。空气阻力系数可通过颗粒与流体相对运动时的颗粒雷诺数的函数求得，即标准阻力公式：
(6)
式中：Rep为颗粒雷诺数，
无量纲参数；μ为空气黏性系数，kg/(m·s)；Ut=|u|为颗粒与空气之间相对速度的绝对值。
1.3 粉尘颗粒物两相流物理模型
预测粉尘污染物在某一区域、某一时刻的浓度，是评估环境空气质量基础。通过设定初始流场、污染物排放源、通风口的位置、通风量等对粉尘污染物的影响规律进行分析。
气体在受限空间内的扩散问题是无化学反应的单相扩散问题，采用守恒形式的三维N-S方程和组分输运方程进行描述，如式(7)所示：
(7)
其中：
hi=
cp,idT。
式中：
为速度矢量
为密度，kg/m3；p为压强，Pa；T为温度，K；μ为黏性系数，kg/(m·s)；k为导热系数，W/(m·K)；Yi为组分i质量分数；
为应力张量；I为单位张量；E为总能量，J；h为焓，cp为比热，J/K；Tref=298.15 K，为参考温度；Sm为质量源项，Sh为能量源项，包括化学反应能量源项Sh,c与体积热源项
当为动量源项；Ri为组分i的化学反应净生成率；Si为组分源项；
为组分i的扩散通量。
当厂房格局复杂时，气流流动会呈现出复杂的湍流特征。湍流的运动存在随机变化，具有非线性和多尺度特性，通过对湍流运动进行建模可掌握厂房内湍流流场的详细信息。湍流流动的模拟需要2方面的必备要素：1)湍流运动中出现的最大涡要求被包含在计算区域内；2)计算网格的尺度要足够小，小到足以分辨最小涡的运动。通过Navier-Stokes方程可直接求出比网格尺度大的湍流运动，通过建立模型来模拟小尺度的涡对大尺度运动的影响，从而形成了目前的大涡模拟方法(large eddy simulation，LES)。
在计算时，LES方程通过有效滤掉比过滤网格小的漩涡进而得到大涡的动量方程。过滤的变量定义为：
(8)
式中：D为流场区域；G为决定过滤尺寸的函数，进行过滤操作。
(9)
式中：V为计算单元的体积。
过滤函数G(x,x′)定义为：
(10)
用LES计算可压缩流体还不现实，目前主要用于不可压缩流体，通过过滤不可压缩N-S方程进一步得到方程见式(11)—(12)：
(11)
(12)
其中τij为亚网格张力，定义为：
(13)
在此过程中可得到亚网格黏性力μeff=μ+μt，其表达式为：
ueff=μ[1+H(x)]1/3
(14)
H(x)为Heaviside函数：
(15)
其中
(16)
(17)
式中：理论给出Crng=0.157；C=100。
湍流中颗粒扩散方程如式(18)所示：
(18)
积分时间与颗粒的湍流扩散率成正比。颗粒在流动过程中处于湍流状态的时间随着T值的增大而增长。
(19)
离散相模型的模拟计算，需要通过计算足够多的颗粒轨迹情况来模拟，此过程比较复杂，为了便于操作，可以通过较少数目的颗粒轨迹来近似模拟大量气溶胶颗粒的扩散趋势。本文中的粉尘颗粒采用颗粒轨道模型(DPM)来计算。DPM模型用于计算颗粒相的运动轨迹，是在拉格朗日坐标下模拟流场中离散颗粒相，还可以用于计算颗粒引起的热量/质量传递。在拉氏坐标系下，颗粒轨道模型通过积分颗粒作用力微分方程就可得到，颗粒作用力微分方程可表示为：
(20)
式中：Fx为附加加速度项，为附加力除以颗粒质量，在不同情况下代表的力不同，文中不考虑该作用力影响。
2 数值方法及模拟
2.1 数值方法
针对以上物理模型，本文采用非稳态、分离隐式解法迭代求解，时间采用一阶隐式，空间上采用二阶迎风格式离散控制方程，采用PISO压力-速度耦合法，在计算中进行畸变修正和相邻修正。
控制体积V中标量φ的瞬态守恒型输运方程的积分形式：
(21)
式中:为速度矢量；
为表面面积矢量；Γφ为标量φ的扩散系数；
φ为标量φ的梯度；Sφ为单位体积φ的源项；ΔV为单元体积；ΔA为包围单元面的面积。
对于时间格式的离散，采用一阶离散格式:
(22)
式中:上标0表示物理量在t时刻(时间步开始时)的值，f为介于0与1的加权因子，f=0、1/2、1时分别为显式、Crank-Nicolson、全隐式时间积分方案。本文采用自编程序来实现数值模拟方法。采用面向对象的C++语言编写，采用MPI与OpenMP共同实现大规模高性能并行计算，空间离散采用有限体积方法，二阶迎风格式，采用基于压力的分离求解器。对于颗粒离散相采用基于拉格朗日体系的颗粒轨道追踪方法，相间有完整的相互作用模块。
2.2 几何物理模型
几何物理模型是根据厂房尺寸建立的，如图2所示。厂房四周有窗户，一侧有门。初始假定操作人员采用蹲或站的姿态位于厂房操作台中央部位对产生煤炭粉尘污染的装置进行操作(如皮带机、碎煤器等)，粉尘污染装置同样位于房间中间部位，屋内设施布局按正常厂房布局。
图2 厂房几何物理模型
Figure 2 Geometric physical model of the factory building
模型区域的尺寸为25 m×6.35 m×5.2 m，顶部为拱形，拱高为2.0 m，距房顶0.25 m有一送风管道，管道上有3个尺寸为0.3 m×0.2 m送风管口。网格划分选用四面体型，网格尺寸为0.1 m，采用均一化布置，边界条件和流场类型：拱形顶部3个方形进风口设置为压力入口，压力值为-10 Pa，底部排气口设置为压力出口，压力值为-10 Pa，其余边界为固体壁面边界。气体相参数采用空气的参数：密度为1.225 kg/m3，比热容为1006.43 J/(kg·K)，热传导系数为0.2042 W/(m·K)，黏度为1.789×10-5 kg/(m·s)；颗粒相参数：粉尘颗粒的平均直径为50 μm，密度为1500 kg/m3。
3 结果与讨论
目前，对于厂房内粉尘污染物的处理，大多采用设置通风口和加装送风管道的方式对其进行吹除，处理效果不理想。为了在现有厂房设备和结构的基础上，进一步提高处理效果，对现有的送风管道、通风口状态以及房门状态进行合理设置和搭配，优化厂房内粉尘气体的流动结构，进一步提高粉尘气体流动的合理性，使更多的粉尘气体通过出口排出，降低室内污染物浓度。本文选择了目前厂房经常采取的两种措施“房门处于关闭状态，屋顶通风口和下方窗子均处于打开状态”(工况1)和“房间下方窗子、屋顶通风口以及房门均处于打开状态”(工况2)进行数值模拟，研究关键因素对房间内粉尘污染的影响，验证本模型方法的有效性。
3.1 工况1条件下粉尘分布及演化分析
首先根据现有房间布局，对房间内的气体运动状态进行计算，从而得到稳态流场分布，如图3所示。
图3 稳态流场分布
Figure 3 Distribution chart of steady flow field in space
图3a为y=0截面处速度矢量分布情况。可知：气流经进气孔以入口速度喷入内部空间中，形成内部空间气流的循环流动分布，气流由上喷入后到达地面，由于壁面阻挡作用会向四周扩散，其中中心进气口正对准操作台，造成操作台表面气流方向为水平方向，离开操作台后气流与周期气体混合。图3b为y=-3.1, 3.1截面处速度矢量分布，可看到气流经前后两面出气孔排出空间。图3c为x=-12, -6.2, 0, 6.2, 12截面处速度矢量分布，可看到进气口和排气孔之间气流贯穿流通的现象。
在计算得到稳态流场分布后，加入粉尘颗粒进行颗粒粉尘污染物分布的计算，图4为不同时刻颗粒污染物在空间中的分布状况。可知：由于顶部气流的吹动作用，颗粒在释放后首先在工作台表面水平方向运动，在离开工作台进入空间后在上方气流的挤压之下，颗粒大部分向空间底部运动，造成颗粒在底部沉积，当底部颗粒沉积到一定程度后，由于附近气流漩涡的作用突然将颗粒向上扬起，造成颗粒在空间中任何位置都存在分布，但大部分颗粒仍位于空间底部。由图5不同位置颗粒物浓度分布情况可更为清晰地看出，由于粉尘源放于工作台靠左边区域，颗粒主要排出到厂房的左方空间中，左方空间中污染最为严重，建议将粉尘源置于排气孔附近，同时操作人员位于上风向方向。
图4 工况1下不同时刻粉尘颗粒在空间中的分布
Figure 4 Distribution of dust particles in space at different time under condition 1
图5 工况1下不同位置颗粒浓度分布云图
Figure 5 Particle concentration distribution at different positions under condition 1
由图5可知：由于颗粒的沉降作用，颗粒浓度大的区域主要位于空间的底部，最大颗粒浓度已经达到16.58 mg/m3，根据指导性国家标准GBZ 2.1—2007《工作场所有害因素职业接触限值 化学有害因素》中的规定，工作场所中粉尘容许浓度≤4 mg/m3，因此对于操作人员来说，呼吸的剂量严重超过阈值限制，威胁身体健康，应采取措施降低粉尘颗粒浓度。为了验证模型和算法计算的准确性，选取了工作人员正常进行装置操作一段时间后，厂房内流动处于稳定状态时，对实际情况下厂房内的颗粒浓度分布进行测量，测得厂房空间底部最大颗粒浓度为17.74 mg/m3，相对误差为7.0%(在10%以内)，结果基本吻合。表明采用新的模型和方法不仅可以得到粉尘颗粒物在空间中的实时运动状态，同时粉尘浓度等定量值也可较好捕获。
3.2 工况2条件下粉尘分布及演化分析
为能迅速降低厂房内的粉尘污染物的浓度，将厂房门打开进行通风，分析该条件下，污染物在厂房内的分布及演化规律，以便操作人员掌握更好的防护措施。
图6为不同时刻粉尘颗粒在空间中的分布情况。可知：在计算初始时刻，粉尘粒子首先在操作台部位聚集，在气流的吹动下，粒子向左方运动，在5 s时粒子运动离开操作装置后受到周围强烈气流的吹动作用迅速向附近空间扩散，在15 s时充满整个房间。同时，由于厂房门开放通风的原因，粒子向房门运动的速度明显较大，很快运动到厂房门附近，并排出厂房外，但厂房门的排放尺寸有限，在15 s时粒子充满整个空间时，粒子浓度处于高峰值，但随着粒子通过下方8个窗户排出厂房，房内的粒子浓度又很快降低，直至在空间内达到平衡状态，如25 s时粒子空间分布状况。
图6 工况2下不同时刻颗粒空间分布情况
Figure 6 Spatial distribution of particles at different time under condition 2
图7为粉尘颗粒浓度在不同时刻、位置处的分布状况。可知：粉尘颗粒在空间中初始时刻分布浓度较大，随着时间推移，粉尘颗粒通过门窗及通风口的排放，室内浓度逐渐降低至稳定状态，除在粉尘源附近浓度较大外，空间中颗粒浓度分布基本处于均衡状态，因开门打破了空间中气流的分布，增加了内部空间与外部的气流流通，便于内部粉尘的排放。同样对此种工况下厂房内的颗粒浓度分布进行测量，测得厂房空间底部最大颗粒浓度为3.6 mg/m3，相对误差为8.2%(在10%以内)，验证了模型与方法在不同工况下计算的准确性。
图7 工况2下厂房内不同时刻粉尘颗粒空间浓度分布情况
Figure 7 Spatial concentration distribution of dust particles in the plant at different time under condition 2
图8为房间内操作人员位置处的污染物浓度随时间变化关系。在工况1条件下，操作人员位置处污染物浓度一直持续增加，顶部和下部窗口通风未明显降低污染物浓度，浓度增加到16 mg/m3以上，对操作人员造成极大的伤害。随着房门的打开，房间内气流产生较大的变动，上方进入较多的气体量通过房门和下方排气孔流出室外，大量的粉尘颗粒通过此方式被带出室外，室内浓度一直保持在安全线以下。
图8 房间内操作人员位置处的污染物浓度随时间变化关系
Figure 8 The pollutant concentration changes with time in the operator′s position in the room
通过上述粉尘颗粒扩散过程计算分析可知：受上方窗户进气的影响，离地面越近粉尘颗粒浓度越大，因此，操作人员应尽量位于窗户下方操作，处于气体流动的上风向，不易受粉尘物质侵害；同时操作人员适宜采用站姿，使口鼻处于空间较高位置处，降低粉尘污染的风险；而煤炭源应处于窗户附近，释放粉尘颗粒物后直接排出厂房；避免将煤炭置于进气孔流通路径，否则不仅会加速粉尘污染物的释放速度，同时将造成整个厂房内的粉尘污染；此外将厂房门置于开放状态，便于内部气流流动及粉尘的排放。
4 结 论
针对火力发电厂物料输送过程中产生大量粉尘的问题，从产生粉尘漂浮运动的机理出发，采用多相流体动力学理论，综合考虑输送设备的空间分布、产生扬尘的源头、外部气流的影响等多方面因素。同时，采用新型计算流体力学方法，结合气体粉尘运动机理，对物料输送过程中的粉尘沉降、气体流动、颗粒分布等过程进行数值模拟，得到以下结论：
1)本文所建立的颗粒释放模型、沉降模型、颗粒-气体两相流模型以及所采用的数值方法可很好地捕捉到粉尘颗粒在空间中从释放到运动、沉积、排出等细节现象，获得不同时刻粉尘颗粒空间分布、浓度分布等参量，分析不同通风状况对室内粉尘分布的影响，结果与实测结果基本吻合，表明本文所建立的模型和方法对于燃煤电厂粉尘运动规律预测具有适用性和有效性。
2)通过对不同通风状况下室内粉尘分布浓度的对比，可知合理调整厂房通风位置，可有效降低现场粉尘浓度；同时，由室内粉尘在空间不同位置处的分布规律可知，工作人员通过改变操作位置或改变操作方式，可有效降低受污染的程度。
[1] 孙秀艳. 国家环境保护标准“十三五”发展规划印发[J]. 绿色包装, 2017(4):27.
[2] LIN C. Research on dust control of coal conveying system in coal-fired power plant[J]. Energy & Energy Conservation, 2017(11):112-113.
[3] 张忠海.电厂输煤煤粉尘成因分析及防治措施[J].神华科技,2019,17(5):94-96.
[4] 刘伟.论电厂输煤系统扬尘成因及除尘对策分析[J].内蒙古科技与经济,2019(15):108-110.
[5] 杜珩.关于燃煤电厂粉煤尘的研究[J].广东化工,2019,46(4):117-121.
[6] 吴维峰. 袋式除尘器内部流场数值模拟分析研究[D].哈尔滨：哈尔滨理工大学,2019.
[7] 王为术,张春杰,王鑫,等.静电除尘器飞灰颗粒除尘特性的数值研究[J].科学技术与工程,2017,17(23):195-200.
[8] 奇玉玺.湿式电除尘技术及其在燃煤电厂中的应用[J].电子世界,2017(16):130-132.
[9] 张帆.燃煤电厂湿式电除尘技术及应用分析[J].科技风,2018(22):185-186.
[10] 郑善江.湿式电除尘技术及在电厂应用与发展探讨[J].科技创新与应用,2019(3):157-158.
[11] 李春亮,卢艳峰,华岁喜.无动力除尘在电厂输煤系统中的应用[J].中国高新科技,2017,1(6):56-58.
[12] XIN T, HU Q, SUN S, et al. Study on technologies of dust emission control in coal fired power plant[J]. Environmental Engineering, 2015.
[13] 蔡松.燃煤电厂除尘技术的分析比较[J].科技视界,2016(26):303-314.
[14] 安连锁,王金平,郦建国,等.中国燃煤电厂电除尘技术发展及应用综述[J].中国电力,2018,51(4):115-123.
[15] AN L S, WANG J P, LI J G, et al. Development and application overview of electrostatic precipitation technology for coal-fired power plant in china[J].Electric Power,2018,51(4):115-123.
[16] 张成国,王建.火力发电厂输煤系统粉尘综合治理技术研究[J].华电技术,2019,41(4):63-65.
[17] XU H, LI X, MIN L, et al. Research on Complex Multi-pollutants Control Technology in a Large-scale Coal-fired Power Plant[M]. Electrostatic Precipitation, 2009.
[18] 杨健晖.电厂电除尘节能减排技术改造工作研究[J].科技资讯,2018,16(5):48-50.
[19] 梁俊宁,卢立栋,柳玲.低热值煤气发电厂污染物排放特征[J].科学技术与工程, 2013,13(32): 9587-9591.
[20] 王铁营. 电除尘器气流分布与粉尘运动特性的数值分析[D].吉林：东北电力大学,2015.