随着风电行业的迅猛发展,风电机组得到广泛应用,风机火灾事故随之增多,由于风机安装在高空,一旦发生火灾,地面无法及时有效灭火,给风机带来毁灭性破坏,还可能引燃附近的草场林场,带来更大的社会危害。由于风电机组火灾事故不断,世界各国对风机的消防越来越重视,德国制订了DIN EN 50308-2005《风力发电机组防护措施的设计、操作和维修要求》、欧洲消防协会联合会也制定了CFPA-ENo22:2012《风力发电机防火导则》。
由于风力发电机组具有内部结构复杂、工况恶劣、通风换气迅速、火灾风险程度高、火灾扑救难度大等特点,目前国内对于风电机组的火灾研究较少。风机机舱发生火灾后依赖火灾自动探测报警系统,很多学者对大空间的火灾特性开展了相关研究,李玉峰等[1] 通过数值手段对有通风弹药舱内慢速火灾特性进行了研究;张宏等[2]对地下扁平型空间火灾特性进行了研究;Zhang等[3]对停车场火灾进行了数值模拟;邹高万等[4]对船舶机舱火灾热流场特性进行了数值仿真研究。风机火灾会表现出一些独有的规律特性,若能较好掌握风机机舱内的火灾特性,则可采取有针对性的多传感器融合探测等技术进行火灾初期的准确自动探测报警。本文利用数值模拟方法对风机机舱的典型火灾特性进行研究,重点考察了风机机舱内温度及烟雾的变化规律,给火灾报警阈值的设定提供参考,为风机消防设计提供依据,对此类空间内的火灾自动探测技术改进具有指导意义。
1 数学模型 1.1 控制方程组热浮力是风机机舱内火灾热流场的主要驱动力,火场燃烧速度一般远低于1 Ma,火源的热释放速率远远大于耗散,因此忽略耗散项,进而受重力作用的火场的控制方程组:
质量方程
$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho {{u}}} \right) = 0 \text{,}$ | (1) |
动量方程
$ \frac{{\partial \left( {\rho {{u}}} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho {{u}}{{u}}} \right) + \nabla P = \rho {{g}} + \nabla \cdot {{\tau}}\text{,} $ | (2) |
能量方程
$ \frac{\partial }{{\partial t}}\left( {\rho h} \right) + \nabla \cdot \left( {\rho h{{u}}} \right) = \frac{{{\rm D}P}}{{{\rm D}t}} + {\dot q^m} - \nabla \cdot {q_r}\text{,} $ | (3) |
组分方程
$ \frac{\partial }{{\partial t}}\left( {\rho {Y_i}} \right) + \nabla \cdot \left( {\rho {Y_i}{{u}}} \right) = \nabla \cdot \left( {\rho {D_i}\nabla {Y_i}} \right) + {\dot m_i}^m\text{,} $ | (4) |
状态方程
$ P = \frac{{\rho RT}}{M} = \rho RT\sum\limits_i {\left( {\frac{{{Y_i}}}{{{M_i}}}} \right)} \text{。} $ | (5) |
式中:
计算过程中,根据火灾流场特性和计算求解速度精度的要求,对方程组处理的方法比较多,常用的方法具体处理过程文献[4-5]中有详细的推导过程与阐述。通过变形并简化速度散度项
气体连续性方程
$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho {{V}}} \right) = 0 \text{,}$ | (6) |
气体组分方程
$ \frac{{\partial \left( {\rho {Y_l}} \right)}}{{\partial t}} + {{V}} \cdot \left( {\nabla \rho {Y_l}} \right) = - \rho {Y_l}\nabla \cdot {{V}} + \nabla \cdot \left( {\rho {D_l}\nabla {Y_l}} \right) + {\dot m'''_l}\text{,} $ | (7) |
气体动量方程
$ \frac{{\partial {{V}}}}{{\partial t}} - {{V}} \times {{\varOmega }} + \nabla H = \frac{1}{\rho }\left[ {(\rho - {\rho _\infty }){{g}} + \nabla \cdot {\bf{\tau }}} \right] - \left( {\frac{1}{\rho } - \frac{1}{{{\rho _\infty }}}} \right)\nabla \tilde P \text{,} $ | (8) |
气体速度散度方程
$ \begin{split} & \nabla \cdot {{V}} = \frac{1}{{\rho {C_P}T}}\left[ {\dot q + \nabla \cdot \left( {k\nabla T} \right) + \nabla \cdot \left( {\sum\limits_l {{h_{s,l}}\rho {D_l}\nabla {Y_l}} } \right) - \nabla \cdot {{\bf{q}}_r}} \right. - \\ & \left. { \sum\limits_l {{h_l}{{\dot m'''}_l} - \sum\limits_l {{h_l}\nabla \cdot \left( {\rho {D_l}\nabla {Y_l}} \right)} } } \right] + \frac{{\bar M}}{\rho }\sum\limits_l {\nabla \cdot \rho {D_l}\nabla \left( {{{{Y_l}} / {{M_l}}}} \right)} +\\ & \frac{1}{\rho }\sum\limits_l {\left( {\frac{{\bar M}}{{{M_l}}} - \frac{{{h_l}}}{{{C_P}T}}} \right){{\dot m'''}_l}} + \left( {\frac{1}{{\rho {C_p}T}} - \frac{1}{{{P_0}}}} \right)\frac{{D{P_0}}}{{Dt}} \text{,} \\[-23pt] \end{split} $ | (9) |
气体状态方程
$ {P_0}\left( t \right) = \rho TR\sum\limits_l {\frac{{{Y_l}}}{{{M_l}}}}\text{,} $ | (10) |
气体压力方程
$ \begin{split} & {\nabla ^2}H = - \frac{{\partial \left( {\nabla \cdot {\bf{V}}} \right)}}{{\partial t}} - \nabla \cdot {\bf{F}};\\ & {{F}} = - {{V}} \times {{\varOmega }} + \left( {\frac{{\rm{1}}}{\rho } - \frac{1}{{{\rho _\infty }}}} \right)\nabla \tilde P - \frac{1}{\rho }\left[ {\left( {\rho - {\rho _\infty }} \right){{g}} + \nabla \cdot {\bf{\tau }}} \right]\text{。}\!\!\!\!\!\!\!\! \end{split} $ | (11) |
由于对三维几何模型划分的网格为长方体或立方体结构,一般采用交错网格技术以便计算机求解控制方程组,具体做法如下:网格节点上存储标量(如
$ \left\{ \begin{aligned} & {\varPhi ^{{{{{(n}} + {\rm{1)}}}_{\rm{e}}}}} = {\varPhi ^{{{(n)}}}} + \delta t\;f({t^{{{(n)}}}},{\varPhi ^{{{(n)}}}})\text{,}\\ & {\varPhi ^{{{(n}} + {\rm{1)}}}} = \;{\varPhi ^{{{(n)}}}} + \frac{{\delta t}}{2}\;\left[ {f({t^{{{(n)}}}},{\varPhi ^{{{(n)}}}}) + f({t^{{{(n)}}}} + \delta t,{\varPhi ^{{{{{(n}} + {\rm{1)}}}_{\rm{e}}}}})} \right]\text{,}\\ & (n = 0,1, \cdots ,N - 1)\text{。} \end{aligned} \right. $ | (12) |
通过在物理空间上对与时间相关的Navier-Stokes方程组进行滤波后所得到的控制方程组,就可以用大涡模拟技术进行求解计算,对于某一个变量
$ \phi \left( x \right) = \overline \phi \left( x \right) + \phi '\left( x \right)\text{。} $ | (13) |
顶盖滤波、高斯滤波和傅里叶滤波是目前普遍使用的滤波函数。对于用于求解火灾流场大涡模拟,一般采用顶盖滤波,具体表示为:
$ G(x) = \frac{1}{\Delta }H\left( {\frac{1}{2}\Delta - \left| x \right|} \right) \text{。} $ | (14) |
在大空间火灾数值仿真计算中,通常在对火灾燃烧过程的简化处理时,燃烧模型一般不考虑燃烧过程中燃料和助燃氧气混合物燃烧速率上升的过程,一般认为足够快的达到热释放速率的稳定,也即燃烧的稳定,燃烧过程只和混合物的组分占比相关,可以简化计算量,可以得到燃料与助燃氧气混合物瞬时进行的化学反应如下:
$ \begin{split} & {{\rm{C}}_{\rm{x}}}{{\rm{H}}_{\rm{y}}}{{\rm{O}}_{\rm{z}}}{{\rm{N}}_{\rm{v}}} + {v_{{{\rm{O}}_{\rm{2}}}}}{{\rm{O}}_{\rm{2}}} \to {v_{{\rm{C}}{{\rm{O}}_{\rm{2}}}}}{\rm{C}}{{\rm{O}}_{\rm{2}}} +\\ & {v_{{{\rm{H}}_{\rm{2}}}{\rm{O}}}}{{\rm{H}}_{\rm{2}}}{\rm{O}} + {v_{{\rm{co}}}}{\rm{CO}} + {v_{\rm{s}}}{\rm{Soot}} + {v_{{{\rm{N}}_{\rm{2}}}}}{{\rm{N}}_{\rm{2}}}\text{。} \end{split} $ | (15) |
风机机舱实际结构模型比较复杂,主要由机舱平台和平台下部的底部夹层组成,如图1所示。其中机舱平台的主要设备有主轴、齿轮箱、发电机、润滑散热系统、刹车系统、各类电气机柜等;机舱平台底部夹层主要有偏航齿轮箱、传输控制电缆等。在仿真计算时风机机舱简化为4 m×10 m×3.4 m的立方空间,风机机舱内其他设备同样进行合理的简化,简化后火源及火源周围适当区域采用边长为0.025 m尺寸的网格以便尽可能准确地捕捉燃烧信息参数,火源面网格总数量为400个,为20×20[6],火源外的其他区域网格尺寸为边长0.05 m,整个风机机舱空间总共约115万计算网格,计算过程中忽略风机机舱表面和机舱内设备表面的热量传递,假定其均为固体绝热壁面边界,没有热量传递和流失。
数值计算过程中,对湍流进行求解主要是采用Deardorff模型[7]里面的LES技术,而对热辐射的求解则基于有限体积法(FVM),时间微分以显性2阶龙格-库塔法预测/校正进行离散,CFL条件维持计算结果的收敛。
2.2 火源位置和大小根据风机机舱内设备实际布置情况,联轴器附近的润滑油残留或泄漏发生火灾可能性较大,联轴器上端面设置为火源,综合考虑各个因素,最终确定将模拟火源面积拟设为0.5 m×0.5 m尺寸进行仿真。
实际风机火灾燃烧工况非常复杂,以当前的数值模拟手段无法准确地对火源的真实燃烧过程进行完全描述,只能在尽量真实描述整个燃烧反应过程、捕捉燃烧基本参数和规律的基础上,将火源进行有效简化,大量对煤油燃烧实验测试结果显示,0.25 m2燃烧面积油池火焰的稳态热释放率约为200 kW,结合国际通用标准ISO/TS16733的规定,火灾计算中普遍使用
在风机机舱距舱顶壁0.05 m平面上布置15个温度测点(编号为T1~T15))和烟雾高度测线(编号为L1~L5),对计算过程中流场特征参数进行实时记录,温度和烟雾的采样周期均为0.2 s,频率均为5 Hz,各测点的布置如图2所示。
本文用数值仿真计算方法,根据风机机舱火灾流场特性构建合适的控制方程组特性,实现了计算速度和计算精度的有效平衡,计算基本流程如下:
1)依据机舱初始环境参数,赋给流场各特征参数合适的初始值,同时确定适当的计算时间步长;
2)预测并计算下一时刻的参数值;
3)对比步骤1求得的速度散度值和预测的速度散度值,若一致性较好则判断为收敛;
4)依据CFL稳定条件来判定是否收敛,需要说明的是CFL条件是收敛的必要而不充分条件;
5)进而求速度散度的修正值;
6)根据步骤4中解得的散度值与速度修正值散度的一致性来判定收敛性;
7)再次回到步骤2继续计算下时刻的流场参数;
8)经过多次迭代后计算设定时间段内的整个流场。
3 计算结果和分析图3为燃烧过程中热释放速率和总释放热量的积累随时间变化曲线。火源的热释放率按照设定的形式变化,约在50 s达到最大值,此后保持最大值不变,后期的稳定阶段燃烧阶段火源热释放率有一定的波动,这是由于火源燃烧过程的不稳定性、舱内空气流动扰动以及通风口作用等原因共同耦合造成的。整个200 s燃烧过程热时放总量为33 315 kW·s,在50 s之后热量释放总量的增加值和时间基本成线形关系。
图4为风机机舱内烟气和颗粒的分布。在t=50 s时,火灾热烟气及燃烧颗粒已经基本蔓延至整个风机机舱的整个上壁面区域,在t=100 s时烟气层厚度和烟气浓度均不断加大,悬浮的颗粒物越来越多,在火焰正上方区域悬浮颗粒浓度最大。
图5为火灾探测器所在的高度平面上50 s,100 s时刻的温度分布云图。随着火灾的发展100 s时的高温区域明显增大大于50 s时的高温区域。在50 s时在风机长度方向上高温区域基本对称,在100 s时高温区域由火源向通风口下游方向蔓延,蔓延至风机机舱长度方向的尽头。可以看出,火灾发生后温度变化速率较快的区域先发生在火源附近的垂向区域,随后向通风口下游方向蔓延,最后达到末端后由舱壁两边向通风上游扩大,通风口端的中心区域是火灾温度场中最后受影响的区域。
由图6可知,在垂向上,除火焰附近区域外,整个高温区域在风机机舱上层一定厚度的区域,100 s时高温区域大于50 s时,起重机的阻碍作用下,在火灾发展初期起重器一侧的温度远小于另一侧的温度,起重机起到类似防火分区的部分作用,尤其是温度没有充分蔓延时,起重机的隔离作用更明显。
对于风机机舱内温度的变化,可近似看做为温度沿x轴或y轴方向一维变化的过程,根据文献[8-10]中对大空间热流场推导得到的公式来表示风机机舱内温度变化:
$ T^* = {e^{ - Kd*}} {\text{。}} $ | (16) |
式中:温度
图7为火灾发生后在风机机舱上壁面附近的5个温度探测器T1-T5温度随时间变化结果。靠火源附近且处于通风口下游的T4最先温升,升温速率也最快,稳定燃烧后能达到的最高温度也最大。T4测点大约在火灾发生后10 s开始出现温升,在72 s时达到195 ℃,随后在波动中温度基本稳定,波动幅度在20 ℃之内。随着火灾的进行,其余4个测点温度陆续上升,在所有测点温度变化过程中,在同一时刻,测点温度有高到低依次为T4>T3>T5>T2>T1。在火灾发展稳定阶段,5个测点温度基本稳定后,T1测点的波动幅度最小,波动幅度的大小依次为T4>T3>T5>T2>T1。
烟气的遮光性对火灾探测有较大影响,机舱内烟气高度的变化规律需要把握。图8为烟气厚度的变化。随着火灾的发展烟气厚度逐渐增加,在50 s左右各测试线上烟气厚底基本稳定,在火灾发展初始阶段,烟气最先蔓延的区域和温度变化的区域一致,顺序依次为L4>L3>L5>L2>L1;随着火灾的发展,烟气厚度的变化和温度变化的顺序不一致,在35 s时L1烟气厚度在1 m左右。整个火灾过程烟气厚度的波动要大于温度的波动,火灾稳定阶段,靠近通风口测烟气厚度大于通风口下游。
此外对温度测点T1~T5温度梯度进行计算,计算结果显示T4测点的温度梯度变化波动较大,幅值范围为±12 ℃/s,其他测点T1,T2,T3和T5的温度梯度变化幅值远小于T4;距离火源最远的T1测点温度梯度最小,温度梯度幅值小于±2 ℃/s。烟气厚度梯度的变化波动频率更高,在个别时刻波动幅度更大。
图9为通风口热量流出速率,在t=28 s之前通风口的热量流出速率为0,说明在此之前烟气热流场没有运动至通风口,在t=28 s之后通风口热量流出速率快速上升,在32 s时热量流出速率达到最大值约为0.6 kW,随后热量流出速率在波动中下降,在火源燃烧过程的不稳定性和舱内空气流动扰动作用下有一定幅度的波动。热量流出积累显示,整个燃烧过程200 s流出热量为11.8 kW·s,占整个燃烧过程火源热释总放量的0.000 35左右。尽管通风口流出热量占热源总释放量的比例很小,但对于舱内流场的影响还是很大的。
图10为质量流出速率。在整个燃烧过程中,最大质量流出速率为0.18 kg/s,质量流出总量为12 kg。在t=0–30 s时间段内,也即燃烧过程初始阶段,由于燃烧部分缺氧周围空气自发的向燃烧处集聚,从而将外界空气的卷吸进来,主要以质量流入为主,最大流入速率为0.073 kg/s,随着燃烧的持续进行,燃烧进入稳定阶段,烟气弥漫并沉积于开口附近,成为在波动中质量流出。
通过对有通风风机机舱进行的热释放率为200 kW,燃烧时间为200 s的火灾场景数值模拟,考虑通风和风机机舱内设备对温度场和烟雾场的影响,总结如下:
1)火灾发生后温度变化速率较快的区域在火源附近的垂向区域,随后向通风口下游方向蔓延,最后达到末端后由舱壁两边向通风口上游蔓延。在温度探测过程中,火源附近且偏向通风口下游的温度测点最先响应,且温升速率最快,在8 s达到最高温度198 ℃,在燃烧不稳定性及舱内气流扰动的耦合作用下温度波动幅度也最大。远离火源区域且靠近通风口区域的温升最慢,能达到的最高温度也最低。
2)在火灾发展初始阶段,烟气最先蔓延的区域和温升区域一致,烟气的蔓延速率比温度蔓延提前3 s左右。
3)通风口热量的流出占火源热释放总量的比例很低,约为0.035%,但通风口的存在对于机舱内流场特性分布影响较大。
[1] |
李玉峰, 张宏, 霍岩. 有通风弹药舱内慢速火灾特性数值研究[J]. 舰船科学技术, 2015, 37(12): 160-165. LI Yu-feng, ZHANG Hong, HUO Yan. Research on characteristics of slow fires in the ventilation ammunition cabin[J]. Ship Science and Technology, 2015, 37(12): 160-165. DOI:10.3404/j.issn.1672-7649.2015.12.033 |
[2] |
张宏, 李玉锋, 霍岩, 等. 扁平型空间内火灾顶棚射流温度特性[J]. 哈尔滨工程大学学报, 2017(10): 79-85. ZHANG Hong, HUO Yan, et al. Ceiling jet flow temperature characteristics of a flat space in case of fire[J]. Journal of Harbin Engineering University, 2017(10): 79-85. |
[3] |
X. G. ZHANG, Y. C. GUO, C. K. CHAN, W. Y. LIN. Numerical simulations on fire spread and smoke movement in an underground car park[J]. Building and Environment, 2007, 42(10): 3466-3475. DOI:10.1016/j.buildenv.2006.11.002 |
[4] |
邹高万. 船舶机舱火灾热流场特性研究[D]. 哈尔滨: 哈尔滨工程大学, 2005: 86-97 ZOU Gao-Wan. A Study on Thermal Flow Field Characteristics of Plant Room Fires on Ship Vessels[D]. Harbin: Harbin Engineering University, 2005: 86-97 |
[5] |
NAKAMURA Y, KASHIWAGI T, MCGRATTAN K B, et al. Enclosure effects on flame spread over solid fuels in microgravity[J]. Combustion & Flame, 2002, 130(4): 307-321. |
[6] |
WANG H Y. Numerical study of under-ventilated fire in medium-scale enclosure[J]. Building and Environment, 2009, 44(6): 1215-1227. DOI:10.1016/j.buildenv.2008.09.011 |
[7] |
DEARDORFF J. W.. Numerical investigation of neutral and unstable planetary boundary layers[J]. Journal of Atmospheric Sciences, 1972, 29: 91-115. DOI:10.1175/1520-0469(1972)029<0091:NIONAU>2.0.CO;2 |
[8] |
胡隆华. 隧道火灾烟气蔓延的热物理特性研究[D]. 合肥: 中国科学技术大学, 2006. HU Long-hua. Stuidies on thermal physics of smoke movement in tunnel fires[D]. Hefei: University of Science and Technology of China, 2006. |
[9] |
DELICHATSIOS M. A.. The flow of fire gases under a beamed ceiling[J]. Combustion and Flame, 1981, 43: 1-10. DOI:10.1016/0010-2180(81)90002-X |
[10] |
EVERS E., WATERHOUSE A. A complete model for analyzing smoke movement in buildings[R]. Building Research Establishment, BRE CP 69/78, 1981.
|