聚合物固体火蔓延是火灾和燃烧科学的重点研究领域之一。若不考虑材料的相变过程，对聚合材料的传热和热解模型已经有较为深入的研究。基于数值模拟对材料逆流火蔓延的研究，Bhattacharjee等^[1-5]开展了气相流动、辐射、氧浓度和材料厚度等因素作用下的逆流火蔓延研究。另外，Fernandez-Pello等^[6-9]和Wichman等^[10]对火蔓延过程所涉及到的传热与流动偏微分方程模型进行适当假设和简化，得到一定条件下的简化模型。对于聚甲基丙烯酸甲酯(PMMA)逆流火蔓延的数值模拟研究，除上述经典文献外，近年来也开展了不少的工作。Tolejko等^[11]和Johnston等^[12]研究了Lewis数对材料的火蔓延影响，以三维方式对热薄型材料的点火到火蔓延的转变过程进行了模拟，并研究了火蔓延过程的熄火条件。Blasi等^[13-14]通过实验和数值模拟研究了热厚型PMMA逆流火蔓延，及不同风速条件下气固两相的温度场和浓度场变化特征。

PMMA等聚合材料的逆流火蔓延不仅包含了普通固体材料火蔓延过程中的传热模型，还包括了材料的熔融相变过程，另外还可能存在熔融流变和黏塑性等力学性质的变化。由于此类问题的复杂性，目前鲜有相关的研究。Zheng等^[15]最先以数值方法研究材料的相变和火蔓延，结果表明火焰尺寸随着液相导热系数增加而增大，随液相比热容的增加而减小。

本文通过数值模拟的方法研究了PMMA的逆流火蔓延，考虑凝聚相存在糊状区的熔融相变过程，建立了逆流火蔓延传热、熔融相变和热解模型，研究了气相火焰结构和熔融界面形态，并建立了逆流火蔓延简化模型以分析熔融界面形态特征。

1 数值计算模型 1.1 固相传热和熔融相变模型

在前人相变材料的火蔓延模型^[15]中，固相的熔化和相变过程发生在确定的熔融温度，当固体聚合物的温度大于熔融温度时，“固-液”相变过程立即发生。然而，对于聚合物的熔化，相变通常发生在一定的温度范围内而不是在某确切温度^[16]，并且在此过程中材料的导热系数随温度变化。即存在一个所谓的糊状区域(muzhy zone)，材料在此区域既不完全是固体也不完全是液体，而是固相和液相共同存在。

图 1a给出了材料右侧受热融化的一维模型示意图，图中糊状区与固体接触的等温度线为固相线，与液体接触的等温度线为液相线，糊状区温度介于固化温度和液化温度之间。针对聚合材料的燃烧和火蔓延，在固相模型中，考虑存在糊状区的熔融相变过程。图 1b给出了此类聚合材料逆流火蔓延模型示意图。为方便建模表示，以材料表面左端为原点建立坐标系，x轴正方向平行于材料表面向右，y轴方向正垂直于材料表面向上，y>0表示气相区域，y < 0表示材料区域，y=0表示气固交界面。模型计算域长度为L_sx，固相区域高度为L_sy，气相区域高度为L_gy。

针对此类聚合材料熔融相变过程，固相能量方程为

$ \frac{\partial }{{\partial t}}({\rho _{\rm{c}}}{C_p}T) = k{\nabla ^2}T + {q_{\rm{c}}}{w_{\rm{c}}} - {\rho _{\rm{c}}}L\frac{{\partial {f_1}}}{{\partial t}}. $

(1)

其中：ρ_c、C_p和k分别为材料密度、比热容和导热系数；q_c为热解单位质量燃料吸收的热量(单位为J/kg)；L为熔融潜热(单位为J/kg)；w_c为热解反应速率；T表示温度场；材料的含液率用f_l表示，当f_l为1时代表纯液相，当f_l为0时代表纯固相，当f_l在0和1之间代表糊状区。而含液率可近似为温度的函数：

$ {f_1} = \left\{ {\begin{array}{*{20}{l}} {0,}&{T \le {T_{\rm{s}}};}\\ {\frac{{T - {T_{\rm{s}}}}}{{{T_1} - {T_{\rm{s}}}}},}&{{T_{\rm{s}}} < T < {T_1};}\\ {1,}&{T \ge {T_1}.} \end{array}} \right. $

(2)

其中：T_s和T_l分别为固化和液化温度。由式(2)可知，温度小于或等于T_s和温度大于或等于T_l分别代表液相和固相，温度介于T_s和T_l之间代表固液两相共存的糊状区。凝聚相中热导率k随材料的组分而变，可以表示为含液率的函数：

$ k = {f_1}{k_1} + (1 - {f_1}){k_{\rm{s}}}. $

(3)

其中k_s和k_l分别表示纯固态和纯液态的导热系数。假设为单步一阶反应，w_c可以表示为

$ {w_{\rm{c}}} = - {\rho _{\rm{c}}}{A_{\rm{c}}}{{\rm{e}}^{ - \frac{{{E_{\rm{c}}}}}{{RT}}}}. $

(4)

其中：A_c为热解反应的指前因子，E_c为热解反应的化学反应活化能，R为摩尔气体常数。

材料的热解反应一方面产生裂解可燃气体，另一方面使固体材料损失，材料表面因燃料消耗而形变和收缩。这里对模型作一定简化，通过固体密度的连续变化描述消耗过程，材料密度随时间的变化，即∂ρ_c/∂t=w_c。

1.2 气相传热和组分扩散

气相的物理过程主要是热解产生的可燃气体与来流空气中的氧气之间的反应。反应生成热量使气体温度升高，高温气体热量向燃料传递使固体热解，以维持材料的持续燃烧。此过程中伴随可燃气体和氧气组分的扩散和对流、化学反应和能量的传递过程。另外，由于本文模拟的尺寸较小，因此不考虑火焰辐射作用。建立气相的可燃气体和氧气的组分扩散方程，且可燃气体与氧气的反应速率用Arrihenius公式近似。

热解可燃气体和氧气的组分扩散方程为

$ \begin{array}{*{20}{c}} {\frac{{\partial ({\rho _{\rm{g}}}{Y_i})}}{{\partial t}} + u\frac{{\partial ({\rho _{\rm{g}}}{Y_i})}}{{\partial x}} + v\frac{{\partial ({\rho _{\rm{g}}}{Y_i})}}{{\partial y}} = }\\ {{\rho _{\rm{g}}}D\left( {\frac{{{\partial ^2}{Y_i}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{Y_i}}}{{\partial {y^2}}}} \right) + {\mu _i}{w_{\rm{g}}}.} \end{array} $

(5)

气相能量方程为

$ \begin{array}{*{20}{c}} {\frac{{\partial ({\rho _{\rm{g}}}{C_{p{\rm{g}}}}T)}}{{\partial t}} + u\frac{{\partial ({\rho _{\rm{g}}}{C_{p{\rm{g}}}}T)}}{{\partial x}} + v\frac{{\partial ({\rho _{\rm{g}}}{C_{p{\rm{g}}}}T)}}{{\partial y}} = }\\ {{k_{\rm{g}}}\left( {\frac{{{\partial ^2}T}}{{\partial {x^2}}} + \frac{{{\partial ^2}T}}{{\partial {y^2}}}} \right) + {q_{\rm{g}}}{w_{\rm{g}}}.} \end{array} $

(6)

其中：Y_i代表燃料或者氧气的组分质量分数, i=f代表燃料，i=o代表氧气；μ_i=(M_iν_i)/(M_fν_f)为气相反应物的化学计量比，M_i代表相对分子质量，ν_i为化学反应式方程式中的计量系数，显然μ_f=1, 参考文^[22]，μ_o=4；C_pg为气体比热容；ρ_g为气体密度；D为组分扩散系数；q_g为气相燃烧热。

w_g为气相燃料与氧气的化学反应速率，可以用Arrhenius公式近似。设气相反应为单步二阶不可逆反应，则气相化学反应速率为

$ {w_{\rm{g}}} = - \rho _{\rm{g}}^2{Y_{\rm{o}}}{Y_{\rm{f}}}{A_{\rm{g}}}{{\rm{e}}^{ - {E_{\rm{g}}}/RT}}. $

(7)

其中：A_g为气相反应的指前因子，E_g为气相反应的化学反应活化能。

式(5)和(6)中气相的流动(u, v)受气体流动方程控制。部分已有的表面火蔓延模型^{[15, 17-19]}中，通常将气相的流动近似为Oseen流，无需耦合求解速度场与能量场，可以简化求解过程，同时降低了数值求解的不稳定性^[15]，因此本文将Oseen流作为气相流动近似。

本文考虑竖直向下火蔓延的情况，此时外界气流速率由高温火焰诱导的自然对流引起^[20]，其来流速率U_∞可近似为

$ {U_\infty } = {C_{{\rm{bc}}}}{[{a_{\rm{g}}}g({T_{\rm{f}}} - {T_{\rm{v}}})/{T_\infty }]^{1/3}}. $

(8)

其中：a_g和g分别为气体扩散系数(k_g/(C_pgρ_g))和重力常数；T_f为绝热火焰温度(本文取1 810 K)；T_v和T_∞分别为蒸发热解温度和环境温度；C_bc为向下火蔓延的浮力羽流系数，本文取0.57^[21]。

1.3 边界条件设置

1) 气相和凝聚相界面的边界条件。

在气相的下边界，固体热解产生的可燃气体向气相扩散。下边界的传热主要包括气相和固相热传导及固体表面向环境的辐射。组分和温度方程的边界条件可以表示为

$ \begin{array}{l} \left\{ {\begin{array}{*{20}{l}} {{\rho _{\rm{g}}}D{{\left. {\frac{{\partial {Y_{\rm{f}}}}}{{\partial y}}} \right|}_{0 + }} = (1 - {Y_{\rm{f}}}{|_{0 + }}){{\dot m}^{\prime \prime }};}\\ {{\rho _{\rm{g}}}D{{\left. {\frac{{\partial {Y_{\rm{o}}}}}{{\partial y}}} \right|}_{0 + }} = (0 - {Y_{\rm{o}}}{|_{0 + }}){{\dot m}^{\prime \prime }};}\\ { - {k_{\rm{s}}}{{\left. {\frac{{\partial T}}{{\partial y}}} \right|}_{0 - }} = - {k_g}{{\left. {\frac{{\partial T}}{{\partial y}}} \right|}_{0 + }} + \varepsilon \sigma ({T^4} - T_\infty ^4);} \end{array}} \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} y = 0,\quad 0 < x < {L_{{\rm{s}}x}}. \end{array} $

(9)

其中：ε为材料表面发射系数，σ为辐射常数，$\dot m″$为气相下边界的固体热解挥发的可燃气质量流率。

若假设固体热解产生的可燃气体全部通过固体上边界扩散到气相，则沿气相下边界各位置的质量流率可以表示为$\int_{-L_{\mathrm{s} y}}^{0} w_{\mathrm{c}} \mathrm{d} y$，L_sy代表材料厚度。而实际情况下，固体材料产生的可燃气体不能立即到达材料表面，并且可能残留在固体内部，因此实际流入气相下边界可燃气体流率要小。使用一个(0，1)范围内的渗透系数K_e, 对气相下边界的可燃气体质量流率进行修正，因此$\dot m″$的最终结果为

$ {\dot m^{\prime \prime }} = {K_{\rm{e}}}\int_{ - {L_{{\rm{s}}y}}}^0 {{w_{\rm{c}}}} {\rm{d}}y. $

(10)

2) 针对计算区域的其他边界条件。

来流氧气由气相区域的左边界进入求解控制区域，因此气相左边界的氧气浓度为固定值，且与来流空气中的氧气浓度相同。气相左边界的可燃气浓度设定为0。另外，左边界气相和固相的温度均固定为环境温度；氧气和可燃气体组分方程的上边界和温度方向的上边界均采用零梯度边界条件；氧气和可燃气体组分方程的右边界和温度方向的右边界均采用零梯度边界条件，固相温度方程右边界和下边界均采用零梯度边界条件。

1.4 计算方法与模型参数

本文将通过数值模拟的方法，研究材料的逆流火蔓延过程，并考虑材料熔融相变过程。以PMMA的物性参数作为数值模型的材料输入参数^[22]，输入参数如表 1所示。

计算区间的x和y方向尺寸分别为24和21 mm。采用有限差分算法对以上偏微分方程进行求解，并采用均匀正交网格，网格尺寸为0.1 mm×0.1 mm。其中，凝聚相y方向网格数量为27(对应厚度为2.7 mm)，气相y方向网格数为183，气相与凝聚相在x方向共用同一网格尺寸，x方向网格个数均为240，时间步长Δt取值为0.2×10^-4 s。

气相和凝聚相方程中，扩散项使用二阶中心差分，气相对流项使用迎风格式，采用全隐式离散方法，方程中的非线性化学反应源项用同一时间步的上一次迭代值替代。通过气相下边界温度梯度，获得固相上边界的温度梯度边界，求解凝聚相温度方程；获得凝聚相热解速率后，进而计算气相下边界的可燃气体质量流率，求解气相浓度方程；将固相上边界的温度代入气相下边界，求解气相温度方程。在同一时间步内迭代以上步骤，待气相温度场、浓度场、凝聚相温度场均达到收敛条件后，进入下一个时刻的求解过程。在每一时间步的迭代求解过程中，气相下边界与凝聚相上边界互为边界条件，从而耦合凝聚相和气相方程。

2 数值计算结果 2.1 火蔓延瞬态传输过程

为实现材料点火到稳定的火蔓延过程，需在材料表面设置一个较大的外界热流$\dot{q}_{\text {ig }}^{\prime \prime}=2 \times 10^{5} \mathrm{W} / \mathrm{m}^{2}$，以提供初始的点火能量。点火源设置在材料表面x=20~21 mm范围内，当点火源所在位置的气相向固体材料表面的热流大于点火源所施加的热流时，则判定为材料已经被点燃，点火成功后初始点火源将被撤除。

图 2a给出了t=1.0、1.6和6.0 s时材料表面温度分布，分别代表点火、非稳定火蔓延和稳定火蔓延阶段的材料表面温度分布。在点火阶段，由于外加的辐射热流较强，热源位置出现较大的温度峰值，且峰值温度约800 K；在非稳定火蔓延阶段，随着外加辐射热源的消失，原点火位置的表面温度下降，且峰值温度向火焰传播方向移动，这说明此时材料已经被成功点燃。

图 2a标记了稳定火蔓延(t=6.0 s)的火焰前锋位置x_f，本文的火焰前锋位置定义为垂直于材料表面方向的温度梯度最大时所对应的横坐标x。

火焰前锋的未燃区温度初始有一个缓慢的上升，随着与火焰前锋距离的拉近，温升速率出现陡增，这主要受气相高温火焰热量向下传递的影响。Fernandez-Pello等^[23]针对热厚型PMMA(厚2.7 mm)逆流火蔓延的材料表面的温度进行了测量。将图 2a的稳定火蔓延表面温度曲线(t=6 s)横坐标作(x-x_f)的平移，转化为火焰前锋位置为坐标原点的温度数据，结果如图 2b所示。将本文模型结果与Fernandez-Pello实验所测的材料表面的温度进行对比，两者一致性较好。

图 3给出了t=1.0、1.6和6.0 s时沿材料表面的y方向温度梯度分布。由图可见，在稳定火蔓延阶段材料表面在火焰前锋位置的温度梯度出现一个峰值。由图 3的稳定阶段(t=6 s)表面温度梯度分布，在火焰前锋两侧，表面温度梯度出现的急剧下降，这说明火焰向材料的传热主要集中于火焰前锋附近。

图 4a和4b分别给出了t=6.0 s时稳定火蔓延阶段材料内部的温度和温度梯度分布。由各深度y位置温度梯度曲线可见，峰值温度梯度从表面向内部减小，表明材料背面附近的传热小于样品表面附近的传热，凝聚相中的传热主要集中在靠近固体表面的区域。图 4b中同时显示了熔融糊状区以及对应的温度梯度值。在表面附近峰值温度梯度的x位置位于熔化区间的下游，如曲线y=0(材料表面)，而在距离表面的较深位置，峰值温度梯度的x位置相对地向固相区的位置移动, 如曲线y=－1.0L_sy(材料背面)。这是由于火焰向材料的传热主要集中于火焰前锋附近，从材料表面附近的燃烧释放的热量可通过凝聚相传递到材料内部，温度梯度反映了热量传递的强度，所以结果显示内部的温度梯度峰值x位置趋于向火焰前锋位置集中，而熔融区主要受温度值的影响，因此在图 4b中呈现出熔融区x位置的变化范围比温度梯度峰值x位置变化范围大。

2.2 凝聚相温度传输过程

图 5a和5b分别给出t为1.6和6.0 s时的固相温度场。为了表示方便，将原坐标系x坐标值与当前时刻的火焰前锋位置x_f相减，得到以火焰前锋为坐标原点的坐标系。因此，在图 5中，横坐标x为0表示火焰前锋位置，x轴坐标大于0表示火焰前锋下游的火焰区域，x轴坐标小于0表示未燃材料区域，本文的后续的结果采用相同的处理方法。

图 5a中，当材料刚被点燃时，火焰前锋位置的等温线相对于火焰后方的等温线密集，这是由于气相燃烧反应开始放热，并向材料表面传递，因此火焰前锋附近出现较大的温度梯度。根据图 5b的稳定的火蔓延阶段t=6.0 s时的温度场，对于熔融温度区间内的糊状区，火焰前缘附近的糊状区最薄，火焰前锋下游糊状区有缓慢增厚，且固相线和液相线具有相似的形态。由图 5b可见，在火焰前锋附近(x坐标为0的位置)液相熔融界面较浅，这是由于此时材料受热时间较短。在火焰区(x坐标小于0的区域)，随着与火焰前锋距离的增大，材料受火焰加热的时间越长，熔融深度越大。

对于熔融界面液相温度线或固相温度线的形态，本文建立简化的分析模型。在火焰前锋下游，在距离火焰前锋Δx位置，用材料的竖直熔融深度Δy表征纵向传热尺度δ_sy。对于液相，假设材料热解区的表面温度近似恒定，特征传热尺度有

$ {\delta _{{\rm{s}}y}} \backsim \Delta y \backsim \sqrt {{\alpha _1}\Delta t} . $

(11)

其中：a_l为液相热扩散系数，Δt为受热时间，“~”表示两侧变量成正比关系。火蔓延速率恒定为V_f，则在距离火焰前锋Δx位置处，材料表面受到火焰加热的时间可由Δt=Δx/V_f计算，将其代入式(11)，得

$ \Delta y \backsim \sqrt {{\alpha _1}\Delta x/{V_{\rm{f}}}} \backsim \sqrt {\Delta x} \cdot \sqrt {{\alpha _1}/{V_{\rm{f}}}} . $

(12)

由式(14)可知，熔融界面的形态可近似为二次函数；图 6给出了固相线与液相线的数值模型和解析模型对比的结果，可以看到，在火焰前锋附近，结果一致性较好。在x>13的范围内曲线存在偏差，这主要是由于该位置为初始点火源，初始点火源对固体的传热与本节的火蔓延解析模型的传热作用不一致，从而造成两者的偏差。另外由式(12)可知，熔融层的深度随着热扩散系数的增加而增大。这直观上理解为热扩散系的增加使传热速率增大，则材料表面向内部热量传递增强，进而增大了固相熔融速率。

2.3 火焰前锋气相温度

图 7给出了火蔓延稳定阶段的气相结构，包括气相方程的温度场源项、温度场源项的火焰前锋局部位置的放大图、气相温度场以及本文模型结果与前人实验结果的对比。这里将气相反应项表示为(q_gw_g)/(ρ_gC_p)，单位为K/s。这样结果不仅可表征气相热解速率w_c，还可代表气相温度场源项。

图 7a中，在外界逆流风速的影响下，由交界面向气相区域的扩散燃料气体与来流气体中的氧气反应，形成图中向下弯曲的“弧形”反应区域。对比图 7c的气相温度场结果，反应区的温度最高。从图 7b可以看出，反应区有一定厚度，这也说明气相中燃烧化学反应速率的有限性，这不同于Wichman等^{[10, 24]}和DeRis等^[25]气相反应速率无限快的简化模型。当化学反应速率较大时，火焰前锋的位置仅由燃料蒸气的产生决定，而蒸发过程受材料受热条件控制，因此文^{[10, 24-25]}研究的火焰向前蔓延特征受材料表面温度分布的影响。如图 7b最大化学反应区并没有与材料表面接触，在材料表面上方约0.2 mm范围内存在一个反应速率接近于零的区域，称为淬火层^[26]，这主要受材料表面附近温度较低的影响。Ray等^[26]研究表明，在该区域产生的预混合火焰可能会在热解表面上形成锚定效应。

图 7c气相温度场的计算结果表明在火焰前锋位置气相具有较大的温度梯度，图 7d给出了Fernandez-Pello^[23]实验所测的火焰前锋不同距离的沿y方向的温度分布与本文模型的结果，可以看到本文结果能较好地模拟温度场的变化特征。需要注意的是，由于本模型中火焰前锋的定义为“固-气”界面竖直温度梯度最大的x位置，这与实验所测的火焰前锋位置有一定的差异^[22]，实验通过火焰图像观察的火焰前锋位置比计算的火焰前锋位置偏大，这会增大模拟与实验的差距，因此需将本模型计算的温度场向后作一定的平移(0.8 mm)修正后，方可与实验进行对比。

3 结论

本文通过数值模拟的方法，研究了熔融聚合材料的逆流火蔓延过程。考虑存在糊状区的熔融相变过程，建立传热、熔融相变和热解模型。主要结论如下：

1) 在材料表面峰值温度梯度的x位置位于熔化区间的下游；在背面，峰值温度梯度的x位置穿过熔化区间，并移动到熔化区间上游的固相区域。

2) 基于凝聚相温度场，获得了熔融相界面和糊状区形态，结果表明：火焰前锋附近的糊状区最薄，火焰前锋下游糊状区有缓慢增厚，并且固化等温线和液化等温线有良好一致性；火焰前锋附近的液相熔融界面较浅，随着与火焰前锋距离的增大，熔融深度逐渐加深。

3) 模拟结果和尺度分析表明：熔融界面的形态可近似为二次函数。

4) 基于气相温度、反应速率以及氧气和可燃气体浓度，分析了气相火焰形态，模型结果能较好地模拟实验火焰温度变化特征。