基于元胞自动机的井巷火灾仿真

D0L:10.13374.issn1001-053x.2013.12.006 第35卷第12期 北京科技大学学报 Vol.35 No.12 2013年12月 Journal of University of Science and Technology Beijing Dec.2013 基于元胞自动机的井巷火灾仿真 李翠平☒，胡磊，侯定勇，张佳 北京科技大学金属矿山高效开采与安全教育部重点实验室，北京100083 ☒通信作者，E-mai让：cpli@ustb.edu.cm 摘要提出了一种基于元胞自动机的井巷火灾可视化仿真方法.在矿井巷道可视化的基础上，通过对火灾元胞进行表 征，综合考虑可燃物类型与投放密度、井巷通风、井巷坡度等因素对井巷火源引燃效果的影响及双扩散作用、井巷通风、 浮力作用和节流作用等因素对火灾烟气蔓延效果的影响，采用概率函数进行元胞自动机建模，构建了表达元胞温度的井 巷火源燃烧模型和表达元胞浓度的井巷火灾烟气蔓延模型.基于火源元胞燃烧演化规则和烟气元胞蔓延演化规则，通过 可视化手段展示了井巷火灾火源燃烧和有害气体浓度的时空发展变化.同时以矿山实际数据进行检验，说明了基于元胞 自动机的井巷火灾仿真的可行性与有效性. 关键词矿山火灾：元胞自动机：计算机仿真：可视化 分类号TD75+2 Cellular automata-based tunnel fire simulation LI Cui--ping☒，HU Lei,,HOU Ding-yomg,ZHANG Jia Key Laboratory of High Efficiency Mining and Safety for Metal Mines (Ministry of Education of China),University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:cpli@ustb.edu.cn ABSTRACT This article introduces a mine fire visualization simulation method based on cellular automata.On the basis of mine tunnel visualization,a tunnel fire combustion model for expressing cellular temperature and a tunnel fire smoke spread model for expressing cellular concentration were built by fire cell characterization.In these models the influence of fuel type and input density,ventilation and roadway slope on mine fire ignition and the effect of double diffusion,ventilation,buoyancy and throttling on fire smoke spread were taken into account,and a probability function was used for cellular automata modelling.Then according to the fire cellular combustion evolution rules and smoke cellular spreading evolution rules,the spatial changes of tunnel fire burning and harmful gas concentration were demonstrated through visualization means.Actual data from a mine proves the feasibility and effectiveness of tunnel fire simulation based on cellular automata. KEY WORDS mine fires;cellular automata;computer simulation;visualization 地下矿山火灾仍是我国矿山典型的灾害，对矿好地解释矿井火灾本质，实验研究又受实验费用高 山生产的正常运行和井下工作人员的生命安全造成昂、实验条件有限、实验过程不具可重复性等客观 了威胁.为控制和降低火灾造成的影响，进行矿井条件的限制，而计算机模拟技术因其计算能力强、 火灾规律的研究具有重要意义-4刂.虽然从20世 模拟条件灵活、可重复性强等特点在矿井火灾研究 纪50年代起就开始了矿井火灾的大量研究，但以中具有一定优势.为此，本文探讨了基于元胞自动 热力学、流体力学为基础的理论研究目前仍不能很 机的井巷火灾仿真模型的构建与可视化的研究. 收稿日期：2013-03-10 基金项目：因家自然科学基金资助项目(51174032)：教育部新世纪优秀人才支持计划资助项目(NCET-10-0225):中央高校基本科 研业务费专项资金资助项目(FRF-TP-09-001A)

第 35 卷 第 12 期 北 京 科 技 大 学 学 报 Vol. 35 No. 12 2013 年 12 月 Journal of University of Science and Technology Beijing Dec. 2013 基于元胞自动机的井巷火灾仿真 李翠平 ，胡 磊，侯定勇，张 佳 北京科技大学金属矿山高效开采与安全教育部重点实验室，北京 100083 通信作者，E-mail: cpli@ustb.edu.cn 摘 要 提出了一种基于元胞自动机的井巷火灾可视化仿真方法. 在矿井巷道可视化的基础上，通过对火灾元胞进行表 征，综合考虑可燃物类型与投放密度、井巷通风、井巷坡度等因素对井巷火源引燃效果的影响及双扩散作用、井巷通风、 浮力作用和节流作用等因素对火灾烟气蔓延效果的影响，采用概率函数进行元胞自动机建模，构建了表达元胞温度的井 巷火源燃烧模型和表达元胞浓度的井巷火灾烟气蔓延模型. 基于火源元胞燃烧演化规则和烟气元胞蔓延演化规则，通过 可视化手段展示了井巷火灾火源燃烧和有害气体浓度的时空发展变化. 同时以矿山实际数据进行检验，说明了基于元胞 自动机的井巷火灾仿真的可行性与有效性. 关键词 矿山火灾；元胞自动机；计算机仿真；可视化 分类号 TD75+2 Cellular automata-based tunnel fire simulation LI Cui-ping , HU Lei, HOU Ding-yong, ZHANG Jia Key Laboratory of High Efficiency Mining and Safety for Metal Mines (Ministry of Education of China), University of Science and Technology Beijing, Beijing 100083, China Corresponding author, E-mail: cpli@ustb.edu.cn ABSTRACT This article introduces a mine fire visualization simulation method based on cellular automata. On the basis of mine tunnel visualization, a tunnel fire combustion model for expressing cellular temperature and a tunnel fire smoke spread model for expressing cellular concentration were built by fire cell characterization. In these models the influence of fuel type and input density, ventilation and roadway slope on mine fire ignition and the effect of double diffusion, ventilation, buoyancy and throttling on fire smoke spread were taken into account, and a probability function was used for cellular automata modelling. Then according to the fire cellular combustion evolution rules and smoke cellular spreading evolution rules, the spatial changes of tunnel fire burning and harmful gas concentration were demonstrated through visualization means. Actual data from a mine proves the feasibility and effectiveness of tunnel fire simulation based on cellular automata. KEY WORDS mine fires; cellular automata; computer simulation; visualization 地下矿山火灾仍是我国矿山典型的灾害，对矿 山生产的正常运行和井下工作人员的生命安全造成 了威胁. 为控制和降低火灾造成的影响，进行矿井 火灾规律的研究具有重要意义 [1−4] . 虽然从 20 世 纪 50 年代起就开始了矿井火灾的大量研究，但以 热力学、流体力学为基础的理论研究目前仍不能很 好地解释矿井火灾本质，实验研究又受实验费用高 昂、实验条件有限、实验过程不具可重复性等客观 条件的限制，而计算机模拟技术因其计算能力强、 模拟条件灵活、可重复性强等特点在矿井火灾研究 中具有一定优势. 为此，本文探讨了基于元胞自动 机的井巷火灾仿真模型的构建与可视化的研究. 收稿日期：2013-03-10 基金项目：国家自然科学基金资助项目 (51174032)；教育部新世纪优秀人才支持计划资助项目 (NCET-10-0225)；中央高校基本科 研业务费专项资金资助项目 (FRF-TP-09-001A) DOI:10.13374/j.issn1001-053x.2013.12.006

第12期 李翠平等：基于元胞自动机的井巷火灾仿真 1547· 1井巷火灾的元胞表征 S,=3,且结1=0，表示元胞（位，）燃烧熄灭后变 1.1元胞自动机简介 为空元胞 元胞自动机(cellular automata,简称CA),是 2井巷火源燃烧模型 时间和空间离散且物理参量为有限数值集的物理系 2.1火源元胞引燃模型 统模型).元胞自动机不由方程或函数描述，而由 井巷中火灾火源燃烧是可燃物被引燃及燃烧 一系列模型构造的规则构成，其基本组成包括元胞 的过程，燃烧是在可燃物被热流加热至自身温度达 (cell)、元胞空间(lattice)、邻居(neighbor)及规则 到甚至超过着火点温度、所需氧气充足的情况下发 (le).它的特点为时间、空间、状态均离散，每个 生的.为分析井巷内可燃物的燃烧过程，本文采用 变量取有限个状态，每个状态有一确定值：状态改 变规则在时间和空间上是局部的，在离散时间上元 概率函数Pg实现对可燃元胞引燃的控制，同时设 定概率阈值Ph,当Pn≥Ph时，元胞将被引燃， 胞状态同步更新，更新规则确定而统一：且元胞演 化只依赖于本元胞及相邻元胞的状态.元胞自动机 元胞燃烧状态S,值变为2 最初用于生物系统，目前在化学自组织系统、通信、 若设定T为可燃物的着火点温度，T为t 城市发展规划、交通系统和森林火灾等领域应用广 时刻元胞(i,)内可燃物的温度，P为t时刻元胞 泛6-16 被引燃的概率，则当可燃物温度超过着火点温度， 1.2火灾元胞状态 即T,>Tgm时，有 (1)元胞可燃状态.元胞可燃状态用于描述 P=1 (1) t时刻元胞(，)的可燃性，取值0或1.取值0时 为不可燃元胞，即空元胞，空元胞经热流和烟气流 从式(1)可见，元胞内温度达到着火点温度后， 入，成为火灾传热传质作用发生的空间：取值1时 元胞的引燃概率随温度升高、加热时间增长而逐渐 为可燃元胞 增大，温度越高、加热时间越长其被引燃概率越大 (2)元胞燃烧状态.对于可燃元胞即.，=1 巷道火灾发生过程中，火源引燃受可燃物类型 时，不同时刻t的元胞燃烧状态用S表示，取 与投放密度、井巷坡度与形态、井巷通风等因素影 值0、1、2和3，分别表示t时刻元胞（亿）处于未 响，故设定t时刻受上述因素影响的元胞被引燃的 燃烧、正被引燃、正在燃烧和熄灭四个状态 概率为 1.3火灾元胞属性 P。=P(1+P)PwPter. (2) (1)元胞物理属性.元胞物理属性用于描述井 式中，P、Pw和Per分别表示可燃物类型与投放密 巷火灾各空间位置的温度T和烟气浓度C,其中 度、井巷通风和井巷坡度对引燃的影响. T表示元胞(i,)在t时刻的温度，C,表示元胞 (1)可燃物类型与投放密度对引燃的影响.由 (i,》)在t时刻的烟气浓度. 于不同类别可燃物燃烧产生的热量不同，故对可燃 (2)元胞生命属性.火灾元胞不仅具有各自的 物引燃的影响亦不同.可将井巷可燃物分为木材、 状态和物理属性，还具有生命属性，即元胞的状态 胶带、电缆类聚合物及机电设备四类，将可燃物投 与属性是时间的函数.对于不可燃元胞，火灾热流 放密度分为低、中、高三个等级7-20.通过组织 和烟气的涌入代表元胞的产生，当元胞温度或烟气 可燃物类型和可燃物投放密度的矩阵，测算各概率 浓度与空元胞持平时表示元胞的死亡.对于可燃元 P,具体测算结果为木材、胶带、电缆类聚合物及 胞，其元胞状态被引燃、燃烧到熄灭是一个动态的 机电设备四类可燃物在低密度时对应的影响概率均 过程，不断地有元胞产生、发展和死亡，其经历的 为0：中密度时的概率分别为0.2、0.1、0.2和0：高 时间tc称为可燃元胞的燃烧生命周期. 密度时的概率分别为0.4、0.2、0.3和0.1. 1.4元胞演化规则 (②)井巷通风的影响.风是影响火灾燃烧及烟气 基于元胞的生命特性，其元胞状态和属性随时 流动的重要要素之一，井巷中主要是井巷通风速度 间更新变化.为表征井巷火灾，特将元胞约束如下 和通风方向对引燃的影响21-2，其影响概率表示 规则.当元胞(，)的初始状态=0，且=0， 为 表示元胞(i,)未燃烧：当元胞(i,)的燃烧状态 S=1,且5S0=3,表示元胞（位，）在经过生 Pw e(ciCOS0:ca)v (3) 命周期tc时间后熄灭；当元胞(i,)的燃烧状态 式中：V为巷道风速；01为风流方向与火源燃烧发

第 12 期 李翠平等：基于元胞自动机的井巷火灾仿真 1547 ·· 1 井巷火灾的元胞表征 1.1 元胞自动机简介 元胞自动机 (cellular automata，简称 CA)，是 时间和空间离散且物理参量为有限数值集的物理系 统模型 [5] . 元胞自动机不由方程或函数描述，而由 一系列模型构造的规则构成，其基本组成包括元胞 (cell)、元胞空间 (lattice)、邻居 (neighbor) 及规则 (rule). 它的特点为时间、空间、状态均离散，每个 变量取有限个状态，每个状态有一确定值；状态改 变规则在时间和空间上是局部的，在离散时间上元 胞状态同步更新，更新规则确定而统一；且元胞演 化只依赖于本元胞及相邻元胞的状态. 元胞自动机 最初用于生物系统，目前在化学自组织系统、通信、 城市发展规划、交通系统和森林火灾等领域应用广 泛 [6−16] . 1.2 火灾元胞状态 (1) 元胞可燃状态. 元胞可燃状态 I t i,j 用于描述 t 时刻元胞 (i,j) 的可燃性，取值 0 或 1. 取值 0 时 为不可燃元胞，即空元胞，空元胞经热流和烟气流 入，成为火灾传热传质作用发生的空间；取值 1 时 为可燃元胞. (2) 元胞燃烧状态. 对于可燃元胞即 I t i,j = 1 时，不同时刻 t 的元胞燃烧状态用 S t i,j 表示，取 值 0、1、2 和 3，分别表示 t 时刻元胞 (i,j) 处于未 燃烧、正被引燃、正在燃烧和熄灭四个状态. 1.3 火灾元胞属性 (1) 元胞物理属性. 元胞物理属性用于描述井 巷火灾各空间位置的温度 T 和烟气浓度 C，其中 T t i,j 表示元胞 (i,j) 在 t 时刻的温度，C t i,j 表示元胞 (i,j) 在 t 时刻的烟气浓度. (2) 元胞生命属性. 火灾元胞不仅具有各自的 状态和物理属性，还具有生命属性，即元胞的状态 与属性是时间的函数. 对于不可燃元胞，火灾热流 和烟气的涌入代表元胞的产生，当元胞温度或烟气 浓度与空元胞持平时表示元胞的死亡. 对于可燃元 胞，其元胞状态被引燃、燃烧到熄灭是一个动态的 过程，不断地有元胞产生、发展和死亡，其经历的 时间 tC 称为可燃元胞的燃烧生命周期. 1.4 元胞演化规则 基于元胞的生命特性，其元胞状态和属性随时 间更新变化. 为表征井巷火灾，特将元胞约束如下 规则. 当元胞 (i,j) 的初始状态 I t i,j = 0，且 I t+1 i,j = 0， 表示元胞 (i,j) 未燃烧；当元胞 (i,j) 的燃烧状态 S t i,j = 1，且 S t+tC i,j = 3，表示元胞 (i,j) 在经过生 命周期 tC 时间后熄灭；当元胞 (i,j) 的燃烧状态 S t i,j = 3，且 I t+1 i,j = 0，表示元胞 (i,j) 燃烧熄灭后变 为空元胞. 2 井巷火源燃烧模型 2.1 火源元胞引燃模型 井巷中火灾火源燃烧是可燃物被引燃及燃烧 的过程，燃烧是在可燃物被热流加热至自身温度达 到甚至超过着火点温度、所需氧气充足的情况下发 生的. 为分析井巷内可燃物的燃烧过程，本文采用 概率函数 Pign 实现对可燃元胞引燃的控制，同时设 定概率阈值 Pth，当 Pign > Pth 时，元胞将被引燃， 元胞燃烧状态 S t i,j 值变为 2. 若设定 Tign 为可燃物的着火点温度，T t i,j 为 t 时刻元胞 (i,j) 内可燃物的温度，P t 0 为 t 时刻元胞 被引燃的概率，则当可燃物温度超过着火点温度， 即 T t i,j > Tign 时，有 P t 0 = 1 − Ã Tign T t i,j !t . (1) 从式 (1) 可见，元胞内温度达到着火点温度后， 元胞的引燃概率随温度升高、加热时间增长而逐渐 增大，温度越高、加热时间越长其被引燃概率越大. 巷道火灾发生过程中，火源引燃受可燃物类型 与投放密度、井巷坡度与形态、井巷通风等因素影 响，故设定 t 时刻受上述因素影响的元胞被引燃的 概率为 P t ign = P t 0 (1 + Pf)PwPter. (2) 式中，Pf、Pw 和 Pter 分别表示可燃物类型与投放密 度、井巷通风和井巷坡度对引燃的影响. (1) 可燃物类型与投放密度对引燃的影响. 由 于不同类别可燃物燃烧产生的热量不同，故对可燃 物引燃的影响亦不同. 可将井巷可燃物分为木材、 胶带、电缆类聚合物及机电设备四类，将可燃物投 放密度分为低、中、高三个等级 [17−20] . 通过组织 可燃物类型和可燃物投放密度的矩阵，测算各概率 Pf，具体测算结果为木材、胶带、电缆类聚合物及 机电设备四类可燃物在低密度时对应的影响概率均 为 0；中密度时的概率分别为 0.2、0.1、0.2 和 0；高 密度时的概率分别为 0.4、0.2、0.3 和 0.1. (2) 井巷通风的影响. 风是影响火灾燃烧及烟气 流动的重要要素之一，井巷中主要是井巷通风速度 和通风方向对引燃的影响 [21−23]，其影响概率表示 为 Pw = e(c1cosθ1 + c2)V . (3) 式中：V 为巷道风速；θ1 为风流方向与火源燃烧发

.1548 北京科技大学学报 第35卷 展方向的夹角，0°≤01≤360°：C1和c2为调节系 生的热烟气从中心元胞流入相邻元胞的可能性与元 数 胞间的温度差有关，温差越大流向该相邻元胞的可 (3)井巷坡度的影响.由于坡度产生的浮力效 能性越大.本文采用摩尔型元胞，即中心元胞的上、 应，呈现温度上坡方向升温明显、下坡方向减慢的 下、左、右、左上、右上、右下、左下相邻八个元 特征24-2可.故本文考虑井巷坡度对引燃的影响概 胞为该中心元胞的邻居，则中心元胞内的热烟气共 率为 有9种可能流动路径选择.设每条路径的选择概率 Pter=ec0 (4) 为Pk,可得热流流动概率矩阵P=P]3×3,即 式中，0：为井巷坡度，ct为调整系数 P1 P2 P3 2.2火源元胞燃烧演化模型 P= P4 Po P5 (7) 基于上述元胞引燃模型，可得出元胞燃烧的演 P6 P7 化规则为：若元胞（⑦，）的燃烧状态为正在燃烧，其 其中，Pk由下式确定： 相邻元胞从中获取能量使其加热温度上升，在t+1 时刻元胞（位，）的相邻元胞（位±1，j±1）被引燃的概 ∑T-l Pk= ∑k∑T-T可 (8) 率为P+1,一旦邻居元胞被引燃，将不断向外传递 热量，相继引燃其他邻居元胞，从而使得火源继续 式中，飞，l=0,1,·,8,用于标识中心元胞分别与其 燃烧下去 自身及上、下、左、右、左上、右上、右下、左下相 至此，基于元胞引燃模型和元胞燃烧演化模 邻元胞的烟气流动路径；Tk表示第k条路径的烟 型，便建立了火源燃烧模型 气温度 根据热力学第二定律，通过对每个中心元胞的 3井巷火灾烟流蔓延模型 邻居元胞内的温度进行判定，取温度低于中心元胞 3.1烟气元胞蔓延模型 的邻居作为烟气可能流动路径，故得表达式为 火源燃烧产生的烟气随风流在井巷中蔓延的 To-Tk 过程实质是烟流在巷道中的传递过程，该蔓延过程 Pk=∑(-T) (9) 受双扩散作用、井巷通风、浮力作用和节流作用的 影响2.故烟气蔓延的概率函数P。可表示为 (2)菲克扩散作用.菲克(Fick)于1855年提出 Pg=PaPw Pb P (5) 了菲克第一定律，将扩散通量和浓度梯度联系起来 认为在稳定扩散传质过程中，扩散组分在单位时间 式中，P、Pw、乃和P分别为双扩散作用、井巷 内通过垂直于扩散方向的单位截面积的扩散物质流 通风、浮力作用和节流作用对烟气蔓延的影响. 量（称为扩散通量）与该截面处的浓度梯度成正比， 3.1.1双扩散作用影响 有 烟流在井巷中蔓延时由于存在浓度梯度和温 度梯度，烟流蔓延定会有热量和质量的扩散，形成 J=-Ddc (10) dx 双扩散效应.本文从傅里叶热传导和菲克扩散作用 式中，J为扩散通量，D为扩散组分的扩散系数，C 两方面探讨双扩散的作用效果 (1)傅里叶热传导作用.因热传导服从傅里叶 为扩散组分的体积浓度，正 为浓度梯度 同理可确定菲克扩散作用下的热烟气蔓延路 定律，即在不均匀的温度场中，因导热所形成的某 径选择概率为 点热流密度正比于该时刻同一地点的温度梯度.在 Co-Ck 一维温度场中表现为 Pk= ∑(C0-Ck) (11) dT qx=-k (6) dx 根据Boussinesq假设，烟气密度p是温度梯度 式中：q红为热流密度，又称热通量，表示单位时间 与浓度梯度的函数，即 内经单位面积传递的热量： d亚为沿x方向的温度 d 梯度：k为热导率.负号表示热通量方向与温度梯度 p=po(1-0r△T+Bc△C). (12) 方向相反 式中，p0为初始风流密度，△T为温度差，△C为浓 井巷火灾的烟流蔓延遵从傅里叶定律。故在基 度差，1和B。分别为密度对于温度和浓度的体积 于元胞自动机的烟流蔓延模型构建中，火源燃烧产 膨胀系数

· 1548 · 北 京 科 技 大 学 学 报 第 35 卷 展方向的夹角，0 ◦ 6 θ1 6 360◦；c1 和 c2 为调节系 数. (3) 井巷坡度的影响. 由于坡度产生的浮力效 应，呈现温度上坡方向升温明显、下坡方向减慢的 特征 [24−27] . 故本文考虑井巷坡度对引燃的影响概 率为 Pter=ectθt . (4) 式中，θt 为井巷坡度，ct 为调整系数. 2.2 火源元胞燃烧演化模型 基于上述元胞引燃模型，可得出元胞燃烧的演 化规则为：若元胞 (i,j) 的燃烧状态为正在燃烧，其 相邻元胞从中获取能量使其加热温度上升，在 t+ 1 时刻元胞 (i,j) 的相邻元胞 (i ± 1,j ± 1) 被引燃的概 率为 P t+1 ign ，一旦邻居元胞被引燃，将不断向外传递 热量，相继引燃其他邻居元胞，从而使得火源继续 燃烧下去. 至此，基于元胞引燃模型和元胞燃烧演化模 型，便建立了火源燃烧模型. 3 井巷火灾烟流蔓延模型 3.1 烟气元胞蔓延模型 火源燃烧产生的烟气随风流在井巷中蔓延的 过程实质是烟流在巷道中的传递过程，该蔓延过程 受双扩散作用、井巷通风、浮力作用和节流作用的 影响 [28] . 故烟气蔓延的概率函数 Pg 可表示为 Pg=PdPwPbPt. (5) 式中，Pd、Pw、Pb 和 Pt 分别为双扩散作用、井巷 通风、浮力作用和节流作用对烟气蔓延的影响. 3.1.1 双扩散作用影响 烟流在井巷中蔓延时由于存在浓度梯度和温 度梯度，烟流蔓延定会有热量和质量的扩散，形成 双扩散效应. 本文从傅里叶热传导和菲克扩散作用 两方面探讨双扩散的作用效果. (1) 傅里叶热传导作用. 因热传导服从傅里叶 定律，即在不均匀的温度场中，因导热所形成的某 点热流密度正比于该时刻同一地点的温度梯度. 在 一维温度场中表现为 qx= −k dT dx . (6) 式中：qx 为热流密度，又称热通量，表示单位时间 内经单位面积传递的热量；dT dx 为沿 x 方向的温度 梯度；k 为热导率. 负号表示热通量方向与温度梯度 方向相反. 井巷火灾的烟流蔓延遵从傅里叶定律. 故在基 于元胞自动机的烟流蔓延模型构建中，火源燃烧产 生的热烟气从中心元胞流入相邻元胞的可能性与元 胞间的温度差有关，温差越大流向该相邻元胞的可 能性越大. 本文采用摩尔型元胞，即中心元胞的上、 下、左、右、左上、右上、右下、左下相邻八个元 胞为该中心元胞的邻居，则中心元胞内的热烟气共 有 9 种可能流动路径选择. 设每条路径的选择概率 为 pk，可得热流流动概率矩阵 P = [pk] 3×3，即 P =     p1 p2 p3 p4 p0 p5 p6 p7 p8     . (7) 其中，pk 由下式确定： pk = P l |Tk − Tl | P k P l |Tk − Tl | . (8) 式中，k, l=0,1,· · · ,8，用于标识中心元胞分别与其 自身及上、下、左、右、左上、右上、右下、左下相 邻元胞的烟气流动路径；Tk 表示第 k 条路径的烟 气温度. 根据热力学第二定律，通过对每个中心元胞的 邻居元胞内的温度进行判定，取温度低于中心元胞 的邻居作为烟气可能流动路径，故得 pk 表达式为 pk = P T0 − Tk k (T0 − Tk) . (9) (2) 菲克扩散作用. 菲克 (Fick) 于 1855 年提出 了菲克第一定律，将扩散通量和浓度梯度联系起来. 认为在稳定扩散传质过程中，扩散组分在单位时间 内通过垂直于扩散方向的单位截面积的扩散物质流 量 (称为扩散通量) 与该截面处的浓度梯度成正比， 有 J = −D dC dx . (10) 式中，J 为扩散通量，D 为扩散组分的扩散系数，C 为扩散组分的体积浓度，dC dx 为浓度梯度. 同理可确定菲克扩散作用下的热烟气蔓延路 径选择概率为 pk = P C0 − Ck k (C0 − Ck) . (11) 根据 Boussinesq 假设，烟气密度 ρ 是温度梯度 与浓度梯度的函数，即 ρ=ρ0 (1−βT∆T+βC∆C). (12) 式中，ρ0 为初始风流密度，∆T 为温度差，∆C 为浓 度差，βT 和 βC 分别为密度对于温度和浓度的体积 膨胀系数

第12期 李翠平等：基于元胞自动机的井巷火灾仿真 ·1549· 为分析温度梯度导致的热传导作用和浓度梯 热膨胀和火灾燃烧气体产物的大量加入使得火区通 度导致的扩散作用对烟气蔓延的影响程度，引入N 风阻力增加.常用节流系数μ来描述节流程度29， 表示双扩散效应对烟流蔓延的相对影响，有 表示随压力、温度的变化，有 r△T T-To N= μ=8T/8p≈△T/△p= (18) BcAC (13) P-Po 又因火区空气密度ρ是空气温度和浓度的函数，故 当N≥1时，表明热传导为主要作用，此时 T-To 烟流蔓延路径选择概率取式(9)：当N<1时，表明 μ=pV2-po昭 (19) 扩散为主要作用，此时烟流蔓延路径选择概率取式 式中，T为火区当前空气温度，T为初始空气温 (11). 度，p为当前巷道压力，p0为初始巷道压力，P%为 3.1.2井巷通风影响 初始风流密度，V为当前风速，%为初始风速. 井巷通风对热烟气蔓延过程的影响主要体现 因节流力的作用方向与风流方向相反，故中心 在烟气流动随风速变化而变化.故此建立通风对烟 元胞向左上、左、左下邻居流动路径的烟气受到节 流蔓延的影响概率Pw同式(3)，只是此时1为风 流力的正向作用，即推动力：而中心元胞向右上、 流方向与烟气蔓延方向的夹角 右、右下邻居流动路径的烟气受到节流力的反向作 3.1.3浮力作用的影响 用，即阻力.同理引入角度0表示烟气流动路径与 在烟流蔓延过程中，浮力主要影响流动参数 风流方向的夹角，0°≤0≤360°.则浮力对于烟气 (速度、温度和组分等)的空间分布及其随时间的 蔓延的影响概率P为 变化.在热传导和扩散共同作用下，浮力的影响可 P=ecsucost (20) 通过格拉晓夫数(Grashof number,Gr)与雷诺数 (Reynolds number,Re)平方之比来判别.其中格拉 式中，4为烟气蔓延的节流力影响系数；c4为调节 系数. 晓夫数为 Gr=9BATL3 3.2烟气元胞蔓延演化模型 (14) v2 基于上述作用对烟流蔓延的影响效果，可得 式中：g为重力加速度：L为特征尺度：v为运动黏 巷道火灾蔓延演化规则：燃烧元胞内的烟气蔓延 滞系数：△T为温差：B为热胀系数.雷诺数为 方向选择为与其相邻元胞蔓延影响概率最大的第 Re=pV L/n (15) kmax条路径，即符合P-kmx=max{Pg}= max{PaPwPP}.可见，选定概率最大的路径作为 式中：V为流体流速：p为流体密度：刀为黏性系 烟气流动路径，烟气从中心燃烧元胞0流动至kmax 数，=p.则烟流蔓延的浮力影响系数为 路径上的相邻元胞，两个元胞间能质交换平衡后， 元胞内温度和浓度参数发生变化.则t+1时刻元 V2 (16) Y值越小，气体浮力越小，则浮力对于烟气蔓 胞0的温度、浓度属性更新为 延的影响越小. To1=(1-Pokm)To+PokmTax (21) 因浮力方向竖直向上，故中心元胞向左上、上、 右上邻居流动路径的烟气受到浮力的正向作用，即 Ct+1=(1-Po-km)C6+乃6一maxCkmax;(22) 推动力：而中心元胞向左下、下、右下邻居流动路 t+1时刻kmax路径上的相邻元胞的温度、浓度属 径的烟气受到浮力的反向作用，即阻力.为分析浮 性更新为 力对各路径烟气的作用力，引入一角度表示烟气 流动路径与竖直向上方向的夹角，0°≤≤360°， T+=(1-乃6-km)Tn十B一kaT6, (23) 则浮力对于烟流蔓延的影响概率。为 Cki=(1-Po-kmCk+Po-kma Co.(24) Pb=ecs ycos0b (17) 至此，建立了基于元胞自动机的井巷火灾烟气 式中，Y为烟气蔓延的浮力影响系数，c3为调节系 蔓延模型 数 3.1.4节流作用的影响 4 井巷火灾模型可视化 火灾烟流蔓延除受浮力作用即火风压的影响， 基于上述建立的井巷火源燃烧及火灾烟流蔓 节流作用也是影响烟流蔓延的主要因素，因气体受 延的元胞自动机模型，从火灾温度属性与烟流浓度

第 12 期 李翠平等：基于元胞自动机的井巷火灾仿真 1549 ·· 为分析温度梯度导致的热传导作用和浓度梯 度导致的扩散作用对烟气蔓延的影响程度，引入 N 表示双扩散效应对烟流蔓延的相对影响，有 N= βT∆T βC∆C . (13) 当 N > 1 时，表明热传导为主要作用，此时 烟流蔓延路径选择概率取式 (9)；当 N<1 时，表明 扩散为主要作用，此时烟流蔓延路径选择概率取式 (11). 3.1.2 井巷通风影响 井巷通风对热烟气蔓延过程的影响主要体现 在烟气流动随风速变化而变化. 故此建立通风对烟 流蔓延的影响概率 Pw 同式 (3)，只是此时 θ1 为风 流方向与烟气蔓延方向的夹角. 3.1.3 浮力作用的影响 在烟流蔓延过程中，浮力主要影响流动参数 (速度、温度和组分等) 的空间分布及其随时间的 变化. 在热传导和扩散共同作用下，浮力的影响可 通过格拉晓夫数 (Grashof number，Gr ) 与雷诺数 (Reynolds number，Re) 平方之比来判别. 其中格拉 晓夫数为 Gr= gβ∆T L3 v 2 . (14) 式中：g 为重力加速度；L 为特征尺度；v 为运动黏 滞系数；∆T 为温差；β 为热胀系数. 雷诺数为 Re=ρV L/η. (15) 式中：V 为流体流速；ρ 为流体密度；η 为黏性系 数，η=ρv. 则烟流蔓延的浮力影响系数为 γ = Gr Re2 = gβ∆T L V 2 . (16) γ 值越小，气体浮力越小，则浮力对于烟气蔓 延的影响越小. 因浮力方向竖直向上，故中心元胞向左上、上、 右上邻居流动路径的烟气受到浮力的正向作用，即 推动力；而中心元胞向左下、下、右下邻居流动路 径的烟气受到浮力的反向作用，即阻力. 为分析浮 力对各路径烟气的作用力，引入一角度 θb 表示烟气 流动路径与竖直向上方向的夹角，0 ◦ 6 θb 6 360◦， 则浮力对于烟流蔓延的影响概率 Pb 为 Pb=ec3γcosθb . (17) 式中，γ 为烟气蔓延的浮力影响系数，c3 为调节系 数. 3.1.4 节流作用的影响 火灾烟流蔓延除受浮力作用即火风压的影响， 节流作用也是影响烟流蔓延的主要因素，因气体受 热膨胀和火灾燃烧气体产物的大量加入使得火区通 风阻力增加. 常用节流系数 µ 来描述节流程度 [29]， 表示随压力、温度的变化，有 µ = ∂T/∂p ≈ ∆T/∆p = T − T0 p − p0 , (18) 又因火区空气密度 ρ 是空气温度和浓度的函数，故 µ = T − T0 ρV 2 − ρ0V 2 0 . (19) 式中，T 为火区当前空气温度，T0 为初始空气温 度，p 为当前巷道压力，p0 为初始巷道压力，ρ0 为 初始风流密度，V 为当前风速，V0 为初始风速. 因节流力的作用方向与风流方向相反，故中心 元胞向左上、左、左下邻居流动路径的烟气受到节 流力的正向作用，即推动力；而中心元胞向右上、 右、右下邻居流动路径的烟气受到节流力的反向作 用，即阻力. 同理引入角度 θt 表示烟气流动路径与 风流方向的夹角，0 ◦ 6 θt 6 360◦ . 则浮力对于烟气 蔓延的影响概率 Pt 为 Pt=ec4µcosθt . (20) 式中，µ 为烟气蔓延的节流力影响系数；c4 为调节 系数. 3.2 烟气元胞蔓延演化模型 基于上述作用对烟流蔓延的影响效果，可得 巷道火灾蔓延演化规则：燃烧元胞内的烟气蔓延 方向选择为与其相邻元胞蔓延影响概率最大的第 kmax 条路径， 即符合 P0→kmax = max{Pg} = max{PdPwPbPt}. 可见，选定概率最大的路径作为 烟气流动路径，烟气从中心燃烧元胞 0 流动至 kmax 路径上的相邻元胞，两个元胞间能质交换平衡后， 元胞内温度和浓度参数发生变化. 则 t + 1 时刻元 胞 0 的温度、浓度属性更新为 T t+1 0 = (1 − P0→kmax )T t 0 + P0→kmax T t kmax , (21) C t+1 0 = (1 − P0→kmax )C t 0 + P0→kmax C t kmax ; (22) t + 1 时刻 kmax 路径上的相邻元胞的温度、浓度属 性更新为 T t+1 kmax = (1 − P0→kmax )T t kmax + P0→kmax T t 0 , (23) C t+1 kmax = (1 − P0→kmax )C t kmax + P0→kmax C t 0 . (24) 至此，建立了基于元胞自动机的井巷火灾烟气 蔓延模型. 4 井巷火灾模型可视化 基于上述建立的井巷火源燃烧及火灾烟流蔓 延的元胞自动机模型，从火灾温度属性与烟流浓度

.1550 北京科技大学学报 第35卷 属性两个方面，以时间为序对每一元胞状态和属性 自30930至31110方向)，火源处燃烧温度为800℃. 变化进行动态展示，进而实现火灾燃烧与蔓延过程 燃烧产生的烟气成分本文仅考虑C0.该矿-140中 的动态可视化. 段图如图1所示. 4.1可视化流程 基于元胞自动机的井巷火灾可视化仿真流程 为： (1)在井巷可视化的基础上，进行井巷空间网 格的划分： (2)确定火源位置，设定仿真时间步长，同时通 图1某矿-140m中段图 过相关数据计算火灾火源热释放量和烟气产生量： Fig.1 -140 m center section drawing of a mine (3)在通风解算基础上，按时间步长动态更新 火灾所引起的各风路压力、风速、风量等参数，并 4.3井巷火灾温度分布可视化 以此为基础计算巷道通风对火源燃烧和烟气蔓延的 图2为运输平巷中心面的温度分布变化可视 影响： 化效果图.当t=5min时，温度500℃以上的火 (4)基于井巷火源燃烧模型和烟流蔓延模型， 源区长度为10m,火灾影响的巷道长度为60m. 按时间步长动态更新巷道元胞状态参数，编程绘制 随着火源燃烧持续，燃烧产生的热烟气不断涌入巷 以不同颜色与灰度等级标识的火灾温度与烟流浓度 道，使得整个巷道温度逐渐升高.t=l0min时，温度 的动态发展变化，从而实现井巷火灾的可视化仿真. 500℃以上的火源区长度为23m,火灾影响的巷道 4.2数据准备 长度为89m:t=20min时，温度500℃以上的火源 以中国某典型地下和矿山为例，进行火灾元胞自 区长度为38m,火灾影响的巷道长度为135m:t= 动机模型的可视化实现.火灾发生前，井巷风流温 30mim时，温度500℃以上的火源区长度 度20℃.巷道岩壁设为绝热边界，即风流不与巷道 为47m,火灾影响的巷道长度为169m,且 壁发生热交换.火灾由-140m中段运输平巷中电 火源上游影响长度58m,火源下游影响长度 气引火导致运输巷道中胶带燃烧造成.火源位于空 111m.可见火灾发生30min内，人员在火源上 间位置(30990m,22070m,-140m)处.火源强度 游58m外、火源下游111m外不会受到火灾温度 为500kW,在通风情况下（通风风流为x正向，即 造成的伤害 温度/℃ 巷道轴/m 800 30930309503097030990310103103031050 310703109031110 700 t5 min 600 30930309503097030990310103103031050310703109031110 500 10 min 400 30930309503097030990310103103031050310703109031110 300 200 -20 min 30930 3095030970309903101031030310503107031090 31110 100 0 =30 min 图2火灾巷道温度分布 Fig.2 Temperature distribution in the fire tunnel

· 1550 · 北 京 科 技 大 学 学 报 第 35 卷 属性两个方面，以时间为序对每一元胞状态和属性 变化进行动态展示，进而实现火灾燃烧与蔓延过程 的动态可视化. 4.1 可视化流程 基于元胞自动机的井巷火灾可视化仿真流程 为： (1) 在井巷可视化的基础上，进行井巷空间网 格的划分； (2) 确定火源位置，设定仿真时间步长，同时通 过相关数据计算火灾火源热释放量和烟气产生量； (3) 在通风解算基础上，按时间步长动态更新 火灾所引起的各风路压力、风速、风量等参数，并 以此为基础计算巷道通风对火源燃烧和烟气蔓延的 影响； (4) 基于井巷火源燃烧模型和烟流蔓延模型， 按时间步长动态更新巷道元胞状态参数，编程绘制 以不同颜色与灰度等级标识的火灾温度与烟流浓度 的动态发展变化，从而实现井巷火灾的可视化仿真. 4.2 数据准备 以中国某典型地下矿山为例，进行火灾元胞自 动机模型的可视化实现. 火灾发生前，井巷风流温 度 20 ℃. 巷道岩壁设为绝热边界，即风流不与巷道 壁发生热交换. 火灾由 –140 m 中段运输平巷中电 气引火导致运输巷道中胶带燃烧造成. 火源位于空 间位置 (30990 m，22070 m，–140 m) 处. 火源强度 为 500 kW，在通风情况下 (通风风流为 x 正向，即 自 30930 至 31110 方向)，火源处燃烧温度为 800 ℃. 燃烧产生的烟气成分本文仅考虑 CO. 该矿 –140 中 段图如图 1 所示. 图 1 某矿 –140 m 中段图 Fig.1 –140 m center section drawing of a mine 4.3 井巷火灾温度分布可视化 图 2 为运输平巷中心面的温度分布变化可视 化效果图. 当 t=5 min 时，温度 500 ℃以上的火 源区长度为 10 m，火灾影响的巷道长度为 60 m. 随着火源燃烧持续，燃烧产生的热烟气不断涌入巷 道，使得整个巷道温度逐渐升高.t=10 min 时，温度 500 ℃以上的火源区长度为 23 m，火灾影响的巷道 长度为 89 m；t=20 min 时，温度 500 ℃以上的火源 区长度为 38 m，火灾影响的巷道长度为 135 m；t= 30 min 时， 温度 500 ℃以上的火源区长度 为 47 m， 火灾影响的巷道长度为 169 m， 且 火源上游影响长度 58 m， 火源下游影响长度 111 m. 可见火灾发生 30 min 内，人员在火源上 游 58 m 外、火源下游 111 m 外不会受到火灾温度 造成的伤害. 图 2 火灾巷道温度分布 Fig.2 Temperature distribution in the fire tunnel