火灾事故是威胁船舶安全性的主要事故类型之一，约占我国船舶海难事故的10%^[1]。由于船舶舱室布置紧密、空间相对狭小，开口多位于顶部，且可燃物较多，发生火灾时火势蔓延迅速，难以扑灭，极易造成严重的生命财产损失。此外，火灾高温还将导致过火区域的钢材力学性能退化，削弱船体刚度，严重威胁船体结构的安全性^[2]。为了提高船舶抗火设计水平、保障人员安全性和降低财产损失，有必要开展船舶火灾下结构安全性的研究。

目前，国内外对船舶火灾事故下船体结构安全性的研究尚不深入，相关工程计算评估方法和行业标准规范尚未形成。Shetty等^[3]基于统一概率模型对海洋平台的火灾安全性开展了可靠性评估和结构优化设计研究，指出温度载荷和不同温度下钢材的力学性能对于结构安全性有着显著影响。Guedes等^[4]开展了高温下矩形板的极限压缩强度的数值研究，分析不同范围和大小的温度载荷对极限强度的影响，指出当受热面积占结构总面积的50%时，结构的承载能力将迅速降低。付丹文等^[5]采用标准火灾升温曲线模拟了火灾场景下船体结构的温升，分析了10 000 TEU集装箱船在火灾中的极限强度。针对“桑吉”号油轮碰撞着火沉没事件，李陈峰等^[6]考虑火灾温度对材料力学性能的影响，分析了“桑吉”号油轮事故前后的船体承载能力和波浪载荷，并推测了事故沉没的原因。刘云山等^[7]基于火灾场景模拟的温度载荷，开展了舱室火灾下甲板板架的结构热力响应研究。

本文将基于双区大涡火灾场景数值方法和结构热力耦合响应分析方法，建立舱室火灾下结构温度载荷分析方法和船体剩余强度分析方法，为船体结构的抗火设计和火灾事故下结构安全性评估提供技术支撑。

1 基本理论与方法 1.1 火灾场景数值模拟软件与方法

Fire Dynamic Simulation(FDS)是美国国家标准与技术研究院(NIST)开发的一款基于流体动力学的火灾分析软件^[8]，该软件基于双区大涡理论和燃烧模型，通过建立火灾场景和直接求解受火灾浮力驱动的低马赫数流动的N-S方程，可以真实地模拟火灾湍流流动过程，尤其是火灾中的烟气和热传递过程，能够准确模拟火灾场景中温度及烟气的分布与蔓延。FDS的主要控制方程为^[9]：

能量方程：

$ \frac{\partial(\rho h)}{\partial t}+\frac{\partial\left(\rho U_{j} h\right)}{\partial x_{j}}=\frac{\partial}{\partial x_{j}}\left(\varGamma_{h} \frac{\partial h}{\partial x_{j}}\right)+S_{h} $

(1)

N-S方程：

$ \frac{\partial\left(\rho U_{j}\right)}{\partial t}+\frac{\partial\left(U_{j} U_{i}\right)}{\partial x_{j}}=\frac{\partial}{\partial x_{j}}\left(\varGamma_{U} \frac{\partial U_{J}}{\partial x_{j}}\right)+S_{U_{j}} $

(2)

连续性方程：

$ \frac{\partial \rho}{\partial t}+\nabla \rho \boldsymbol{u}=0 $

(3)

化学方程：

$ \frac{\partial(\rho m)}{\partial t}+\frac{\partial\left(\rho U_{j} m\right)}{\partial x_{j}}=\frac{\partial}{\partial x_{j}}\left(\varGamma_{h} \frac{\partial m}{\partial x_{j}}\right)+S $

(4)

式中：ρ为气体的密度，kg/m³；h为气体的焓，kJ/mol；m为气体质量，kg；U_i、U_j为速度分量，m/s；Γ表示交换系数，与模型大小有关；S表示火源。

1.2 基于热传导原理的结构温度载荷分析方法

FDS能够计算获得不同时刻结构和环境的温度、气体压力、流速、CO₂浓度、热通量等火灾场景特性参数。热通量，即热流密度，指单位面积的截面内单位时间通过的热流量。为了准确获得火灾下结构温度载荷，考虑到FDS模型与结构热力响应分析模型节点的不匹配，本文以热通量为环境边界条件，通过热传导分析间接获得结构热力响应分析的温度载荷，其计算原理如下^[10]：

温度场是空间坐标和时间的函数，在直角坐标系中可表示为：

$ T = f(x,y,z,t) $

(5)

物体中等温线较密集的地方温度变化率较大，且温度变化率沿不同方向往往是不同的。在各个不同方向的温度变化率中，沿等温线法线方向的变化率是最大的，其数学表达式为：

$ \operatorname{grad} T=\frac{\partial T}{\partial \boldsymbol{n}} \boldsymbol{n}=\frac{\partial T}{\partial x} \boldsymbol{i}+\frac{\partial T}{\partial y} \boldsymbol{j}+\frac{\partial T}{\partial z} \boldsymbol{k}=\nabla T $

(6)

式中：gradT为温度梯度；$\frac{\partial T}{\partial \boldsymbol{n}}$为等温线法线方向的温度变化率；n为等温面法线方向的单位矢量，指向温度增加的方向；▽为哈密尔顿算子。

导热基本定律，即傅里叶定律的表达式为：

$ q=-\lambda \operatorname{grad} T=-\lambda \nabla T $

(7)

式中：q为热通量；λ为导热系数，表征物质导热能力的一个热物性参数。

沿x、y和z方向的热通量分量分别为：

$ q_{x}=-\lambda \frac{\partial T}{\partial x}, q_{y}=-\lambda \frac{\partial T}{\partial y}, q_{z}=-\lambda \frac{\partial T}{\partial z} $

(8)

则可得物体微元中导入和导出的热量，以x向为例：

$ \mathrm{d} \varPhi_{\mathrm{in}-x}=q_{x} \mathrm{d} y \mathrm{d} z=-\lambda \frac{\partial T}{\partial x} \mathrm{d} y \mathrm{d} z $

(9)

$ \mathrm{~d} \varPhi_{\mathrm{out}-x}=\mathrm{d} \varPhi_{\mathrm{in}-x}+\frac{\partial}{\partial x}\left(-\lambda \frac{\partial T}{\partial x}\right) \mathrm{d} y \mathrm{d} z \mathrm{d} x $

(10)

根据能量守恒定律，即导入微元的热量dΦ_in和内部热源生成的热量dQ，与导出微元的热量dΦ_out和微元内能的增量dU守恒：

$ \mathrm{d} \varPhi_{\mathrm{in}}+\mathrm{d} Q=\mathrm{d} \varPhi_{\mathrm{out}}+\mathrm{d} U $

(11)

式中微元内能的增量可表示为：

$ \mathrm{d} U=\rho c \frac{\partial t}{\partial \tau} \mathrm{d} x \mathrm{d} y \mathrm{d} z $

(12)

式中c为物体比热容。

将式(9)、(10)和(12)代入式(11)处理后可得以热通量为边界，无内部热源的导热微分方程：

$ \rho c \frac{\partial T}{\partial t}=\lambda\left(\frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}+\frac{\partial^{2} T}{\partial z^{2}}\right)=\lambda \nabla^{2} T $

(13)

式中▽²为拉普拉斯算子。

1.3 结构热力耦合响应分析理论

当结构中存在温度梯度时，其内部产生的热应力和在热应力作用下结构产生的塑性变形间的相互作用时出现的位移、应力及应变间的耦合场分析称之为热力耦合分析。结构热力耦合响应分析包括结构瞬态温度场分析和结构热应力场分析^[11]。

结构瞬态温度场分析的控制方程：

$ \dot{\boldsymbol{u}}_{N}^{\mathrm{T}}(t)\left[\boldsymbol{K}_{u} \dot{u}_{N}(t)+\boldsymbol{M}_{T} \dot{T}_{N}(t)-\boldsymbol{F}(t)\right]=0 $

(14)

式中：K_u为结构刚度矩阵；M_T为热力学刚度矩阵；F(t)为受力向量。

结构热应力场分析的控制方程：

$ \boldsymbol{T}_{N}^{\mathrm{T}}(t)\left[\boldsymbol{C}_{u} \dot{\boldsymbol{T}}_{N}(t)+\boldsymbol{M}_{u} \dot{\boldsymbol{u}}_{N}(t)-\boldsymbol{D}-\boldsymbol{R}-\boldsymbol{K}_{T} \boldsymbol{T}_{N}(t)\right]=0 $

(15)

式中：T_N(t)为节点的温度矢量；u_N(t)为节点的位移矢量；C_u为热熔矩阵；M_u为热力耦合矩阵；K_T为热传导矩阵；D为耗散向量；R为热载荷向量。

结合式(14)及式(15)可以得到结构热力耦合响应分析的有限元方法：

$ {\left[\begin{array}{cc} \boldsymbol{K}_{u} & \boldsymbol{M}_{T} \\ \boldsymbol{M}_{u} & \boldsymbol{C}_{u} \end{array}\right]\left[\begin{array}{c} \dot{\boldsymbol{u}}_{N}(t) \\ \dot{\boldsymbol{T}}_{N}(t) \end{array}\right]=\left[\begin{array}{c} \boldsymbol{F}(t) \\ \boldsymbol{Z}(t) \end{array}\right]} $

(16)

$ \boldsymbol{Z}(t)=\boldsymbol{D}+\boldsymbol{R}+\boldsymbol{K}_{T} \boldsymbol{T}_{N}(t) $

(17)

2 目标船基本信息与结构布置

根据英国海军舰船火灾事故的统计数据，50%以上的火灾事故是机舱火灾^[1]。对于水面舰船而言，机舱往往位于船舯，而船舯的船体梁弯矩最大，机舱火灾将严重威胁船体结构的安全性，因此本文以一水面舰船为例，重点开展机舱火灾的模拟和火灾下船体剩余强度的分析。

目标船船长120 m，型宽14.4 m，标准吃水3.87 m，全船有限元模型如图 1所示。该船的机舱位于船舯区域，机舱段长24.15 m，宽14.4 m，高10.8 m(包含3层甲板)，结构布置如图 2所示。表 1为钢材的主要热物理参数^[12]，材料本构模型采用理想弹塑性模型。

3 机舱火灾场景模拟与分析 3.1 机舱火灾场景的构建 3.1.1 模型范围的选取

舱室火灾场景模拟的模型范围为以机舱室为核心的三舱段模型，一方面可以减少建模工作量并提高计算效率，另一方面通过前后舱室连接实现热传导的相对连续以保证计算精度。火灾场景模型中，忽略骨材、加强筋和肘板等局部结构，只计及纵横舱壁、甲板等主要隔断结构，并设置材料属性。

由于通风条件对于火灾的发展有着显著的影响，因此根据目标船的舱室布置和舱门及通风口的设置，在FDS模型中设置相应的内部通风边界。考虑船体与周围环境的热交换，将舱段外部设置为开放边界。图 3为隐去外板的机舱火灾场景模型。

3.1.2 火源的模拟与设置

火源的设置是火灾场景建模的一个重要工作，主要涉及火源类型的选择和热释放速率的确定^[13]。机舱火灾多以柴油火为主，火灾类型为B类火。庚烷的燃烧特性与柴油相近^[14]，其物理化学性质见表 2，因此选择庚烷作为燃料模拟机舱油火，火源设置在机舱内底板靠近前舱壁位置，尺寸为1 m²。

热释放速率指材料在单位时间内燃烧所释放的热量，热释放速率模型主要有t²稳定火源模型、分段平均法、分段线性法和根据质量损失率确定的方法^[15]。其中，t²稳定火源模型采用分段函数模拟整个火灾过程，热释放速率在火灾初期随时间平方增长，达到最大热释放速率后维持一段时间，然后衰减熄灭，其达式为：

$ \bar{Q}(t)=\left\{\begin{array}{ll} \alpha t^{2}, & t \leqslant t_{g} \\ \bar{Q}_{\max }, & t_{g}<t \leqslant t_{d} \\ \bar{Q}_{\max } \mathrm{e}^{-\frac{t-t_{d}}{\tau}} t, & t>t_{d} \end{array}\right. $

(18)

式中：Q_max为燃烧过程中达到的最大的热释放速率，对于汽油、柴油这类易燃物，其最大热释放速率一般为1~3 MW；t_g、t_d分别为热释放速率达到最大值和开始衰减的时刻；τ为衰减时间；α为火灾增长系数，机舱油火属于超快速增长型火灾，推荐取0.187 8 kW·s^-2[15]。

t²稳定火源模型与实际火源较为符合，因此本文采用该模型模拟机舱火灾并忽略减弱阶段，采用稳定燃烧阶段的高温载荷进行结构热力响应分析。

3.1.3 温度测点和环境切片的设置

FDS软件提供了温度测点和环境切片2种获取火灾过程中结构温度和环境温度的方法。其中，温度测点模拟热电偶可以定点测量火灾过程中结构的温度变化；环境切片可以探测切片截面上包括温度、压力、流速等在内的火灾动态特性参数。

为了掌握机舱火灾下火源和舱室环境温度的变化，在火源周围设置了4个结构温度测点，并在中纵剖面、内底板处和2甲板下方等位置设置了环境温度切片，如图 4所示。其中，测点1和2靠近船艉方向，测点3和4靠近横舱壁在火源靠船艉方向，且测点1、2和3位于2甲板开口下方(图 5)，各测点坐标如表 3所示。

3.2 机舱火灾模拟结果与分析 3.2.1 火灾温度

柴油火灾发生后将快速发展并稳定燃烧，因此火灾场景模拟时间设定为600 s。图 6是机舱内底火源周围内底板测点温度随时间的变化，可以发现：4个结构测点的温度变化趋势基本一致，在火灾初期，结构温度增加较为缓慢，150 s左右温度迅速升高，500 s左右温度趋于稳定，说明设置600 s的模拟时间是合理的。

进一步分析图 6中4个测点的温度差异，可以发现：测点1、2和3的温度较测点4温度高，这是由于火灾温度往往随着烟气蔓延，由于测点1、2和3上方有甲板开口，火灾烟气向着通风条件较好的甲板开口移动并蔓延，因此导致通风条件较差的测点4温度较低。

图 7~9是火灾不同发展阶段不同切片上的温度分布情况。其中，图 7是机舱中纵剖面的环境温度分布，图 8和图 9分别是机舱内底和舱室顶部(2甲板)的环境温度分布。结合图 6进一步分析，可以发现：

1) 火灾初期，由于烟气向上发展，火源上方环境温度急剧升高，此时近地面的火源周围的结构温度增加缓慢。

2) 当高温烟气到达舱室顶部，继续向着氧气充足的2甲板开口处移动，机舱室环境温度快速升高；高温烟气顺着舱口进入2甲板和主甲板间的舱室，并继续向舱口或通风口移动。

3) 当达到稳定状态时，舱室顶部的环境温度要比下方高，且最高环境温度达到700 ℃。

因此，舱室火灾的温度分布十分复杂，舱口、通风口等通风条件对于火灾发展和温度分布有着显著的影响，结构温度与近结构的环境温度有相关性，但并不完全一致。

3.2.2 热通量与结构温度载荷

图 10是提取的稳定状态舱段结构表面的热通量，其中着火的机舱舱室壁面上的热通量分布较为密集，其他舱室由于受通风状况的影响，热通量分布较为稀疏，这也符合实际火灾场景下的热通量分布状况。

以节点坐标为索引，将结构表面的热通量映射施加到结构有限元模型(图 2)，以热通量为环境边界，通过热传导分析获得舱段温度载荷，如图 11所示，01甲板和外底板的温度接近环境温度，主甲板温度在400 ℃左右，结构局部最高温度出现在临近2甲板开口的舱壁处，温度接近700 ℃。图 12为该舱壁的结构温度与环境温度的对比，可以发现两者的温度分布基本一致，由于临近2甲板开口，随着高温烟气流动，环境温度和结构表面热通量都较大，导致结构温度也较高。

4 火灾下船体结构剩余强度分析

为了分析火灾高温对船体承载能力的影响程度，对正常状态下和舱室火灾下船体的中垂极限承载能力进行了对比计算分析。

图 13为2种工况下船体梁的弯矩-曲率曲线，曲线的峰值点即对应工况船体极限承载能力，其中正常状态下船体的中垂极限弯矩为7.04×10⁸ N·m，火灾事故下剩余承载能力为6.56×10⁸ N·m，1 m²的机舱池火燃烧600 s导致船体承载能力降低了6.8%。

图 14和图 15分别为2种工况下船体舱段极限状态的结构应力和变形，最大应力趋于室温下材料屈服限。相较于无火灾高温影响的正常状态，火灾事故状态下船体达到极限状态时，甲板的高应力区域范围更小但结构变形更大，船底的应力水平偏低。这是由于火灾高温导致过火区域的钢材材料屈服限和弹性模量等力学性能折减引起的，具体为材料力学性能的降低导致结构的刚度和承载能力降低，过火区域结构应力水平的降低进一步引起剖面中和轴向下偏移，导致中和轴远端的甲板塑性应变增加、结构变形更为显著，同时导致船底的应力水平降低，综合效应最终导致舱段的极限承载能力降低。

5 结论

1) 基于双区大涡火灾场景数值方法适用于船舶舱室火灾场景的模拟。结构测点温度和切片环境温度显示，通风条件对于火灾发展和温度分布有着显著的影响，结构温度与近结构的环境温度有相关性，但并不完全一致。

2) 结构温度载荷是开展火灾下结构热力响应分析的基础，以热通量为环境边界结合热传导分析获得结构温度场，该温度载荷确定方法可行且有效。

3) 火灾高温导致过火区域钢材的力学性能折减和结构应力的重新分布，从而造成船体承载能力的降低。为了提高我国船舶结构抗火设计和火灾下船体结构安全性评估水平，有必要针对一些典型舱室火灾工况开展更为系统的研究分析工作。