基于共轭传热的火灾条件下复合材料储氢瓶的传热行为

程崇律，徐宝鹏，孔凡夫

基于共轭传热的火灾条件下复合材料储氢瓶的传热行为

程崇律，徐宝鹏，孔凡夫

(大连理工大学能源与动力学院，大连 116024)

复合材料储氢瓶为氢燃料电池汽车的关键部件，其在火灾条件下的瞬态传热行为的精准预测对储氢瓶的失效判断至关重要．文中提出了一种基于多重区域的实时双向耦合共轭传热模型，用于模拟复合材料储氢瓶在火灾条件下的瞬态传热特性．共轭传热模型采用大涡模拟(LES)对瓶外流体计算区域的火灾进行了数值模拟，并对复合材料缠绕层、高密度聚乙烯内胆和瓶内高压氢气3个计算区域分别建立了基于有限体积法的三维导热方程．复合材料缠绕层的导热模型考虑了高温热解和热解气的对流及冷却作用，并采用了依赖于复合材料成分和温度的可变物性．结合某商用Ⅳ型储氢瓶的池火火灾试验，验证了所提出的共轭传热模型的可靠性，并对其瞬态传热特性进行了数值研究．

复合材料储氢瓶；大涡模拟；共轭传热；火灾；热解

氢能源具有清洁无污染、利用率高、来源广泛、容易运输等优点，是解决能源枯竭以及坏境恶化等问题的重要途径之一[1]．目前燃料电池汽车普遍采用重容比小的复合材料储氢瓶进行储氢．复合材料储氢瓶由防止高压氢气渗透的高密度聚乙烯或铝合金内胆，以及碳纤维/环氧树脂复合材料缠绕层组成[2]．储氢瓶内氢气压力极高且易燃易爆[3]，对储氢瓶进行火灾试验具有一定危险性．因此，对复合材料储氢瓶在火灾条件下的传热特性开展数值研究具有重要意义.

复合材料储氢瓶在火灾条件下的传热涉及非预混燃烧、由浮力驱动的湍流输运、碳烟的生成、共轭传热、复合材料高温热解等多物理现象，以上多物理过程可分为火灾预测以及储氢瓶热解传热两部分. 文献[4]中普遍采用均匀分布的人为规定最大表面热通量或试验测得的平均表面热流通量作为固体储氢瓶传热模型的边界条件，将火灾预测与储氢瓶传热失效分析分离开来．Saldi等[5]则采用松耦合的方法，使用大涡模拟(LES)方法提高了预测储氢罐表面热通量的准确性，并将计算得到的3个不同位置上随时间变化的表面热通量分布作为基于有限元方法的二维失效分析热边界条件，但仅实现了火灾预测模型-固体模型的单向耦合．然而，储氢瓶在火灾下的极限强度很大程度上取决于储氢瓶表面热通量，对储氢瓶表面热通量的准确预测尤为重要．实际表面热通量具有显著的时间瞬态性以及空间分布不均匀性，以上两类研究由于方法的局限性，无法对固体表面热通量进行准确预测，预测结果与实际表面热通量差异较大，进而影响失效分析的准确性．

双向耦合共轭传热是在每一时间步长内对不同区域(火灾流体区域与复合材料固体区域)内的控制方程进行耦合求解，通过能量方程平衡流体和固体交界面上的热通量，可得到更为准确的储氢瓶表面热通量瞬态分布．目前文献中已有类似研究[6]，但其未考虑复合材料在高温下的特殊热解现象，即复合材料处于高温状态时会热解生成积炭以及热解气[7]，热解气体通过对流作用向外表面进行迁移，影响复合材料层内的温度分布，同时复合材料的组成成分也将发生改变，严重影响复合材料的热力性能和机械性能，从而导致储氢瓶的强度降低[8]．

本文提出一种基于多重区域的实时双向耦合共轭传热模型，同时考虑复合材料的特殊热解行为，并基于开源软件工具包Open FOAM，使用C++编辑语言，开发了相应的传热数值分析计算程序，用于模拟复合材料储氢瓶在火灾条件下的瞬态传热特性．计算区域包括火灾流体区域、碳纤维/环氧树脂复合材料缠绕层区域、高密度聚乙烯(HDPE)内胆区域和高压氢气区域．采用LES对火灾的非定常流动进行模拟，并对复合材料缠绕层区域、HDPE内胆区域和高压氢气区域分别建立基于有限体积法的三维控制方程，复合材料缠绕层的导热模型考虑了高温热解和热解气的对流、冷却作用以及可变物性．结合某商用Ⅳ型储氢瓶的池火火灾试验[9]，对该模型进行了验证，并对其瞬态传热行为进行了详细的数值研究.

1 数学模型

本文针对不同计算区域分别建立了控制方程，在区域界面上通过保证各界面温度和热通量的连续性实现不同计算区域的双向耦合共轭传热．

1.1 火灾预测模型

本文采用大涡模拟方法对火灾的非稳态放热和流动过程进行数值预测．火灾模型求解非稳态、可压缩、包括浮力效应和经过空间滤波的Navier-Stokes (N-S)输运方程组．方程中的亚网格尺度(SGS)源项由单方程SGS湍流模型确定[10]，模拟被滤波掉的小尺度湍流涡对直接求解的大尺度湍流涡的影响．采用扩展的涡耗散模型(EEDM)[11]模拟池火的扩散燃烧过程，该模型假定湍流燃烧过程由燃料的混合过程控制，根据Kolmogorov尺度下的“微结构”与其周围流体的湍流混合速率确定湍流燃烧速率．碳氢化合物燃料在燃烧过程中会产生大量的碳烟，碳烟对辐射换热影响显著，本文采用扩展的烟点模型[11]对其进行模拟，其将碳烟的生长以及碳烟的氧化表达为混合分数和温度的函数，仅求解碳烟质量分数守恒方程，具有可观的计算精度．储氢瓶表面热通量主要来自于扩散燃烧火焰的辐射传热，采用扩展的有限体积法模型[12]对辐射传热进行模拟，碳烟的辐射换热量占总辐射换热量的主要部分，其辐射特性远大于气体辐射特性，采用灰体介质假设．由Radcal程序中的多项式表达式计算不同种类气体光谱吸收系数，由预测碳烟体积分数和温度确定碳烟的吸收系数，同时考虑储氢罐表面的辐射作用．以上模型在前人研究文献中已详细描述，本文由于篇幅有限不再详细介绍．

1.2 复合材料缠绕层传热模型

文献[8]中复合材料传热研究通常仅考虑在厚度方向上的一维传热，本文建立了三维传热模型．其中热解源项考虑了热解引起的质量和能量损失以及热解气对流引起的能量迁移，即

假设热解气对流只发生在复合材料缠绕层的厚度方向，热解气的质量流量可表示为

复合材料和热解气的焓表示为

1.3 HDPE内胆与高压氢气传热模型

由于复合材料的热惯性较大，HDPE内胆在发生失效之前的温升较低，本研究只考虑HDPE内胆区域的导热问题．此外，通过试验[9]发现在火灾条件储氢瓶爆破发生失效前，瓶内氢气压力没有明显升高，说明罐内氢气温升较小．在之前的研究中[5]，采用定值对流换热系数对罐内表面换热进行计算．本文忽略在小温升下氢气浮力效应，即假定瓶内氢气的传热问题为静止条件下的导热问题．因此HDPE内胆和高压氢气的传热控制方程采用相同形式的导热方程.

1.4 界面边界条件

在计算区域的界面处需要保证温度和热通量的连续性．

(1) 在流体区域与复合材料缠绕层区域的界面处满足：

本研究中火灾流动的雷诺数较低，壁面边界层较厚，采用直接求解壁面边界层的方法计算罐表面对流换热．

(2)在复合材料缠绕层区域与HDPE内胆区域的边界面处满足：

(3)在HDPE内胆区域与高压氢气区域的边界面处满足：

1.5 数值求解方法

将各区域的控制方程与上述界面处边界条件耦合，按顺序依次求解火灾区域、复合材料缠绕层区域、HDPE内胆区域和高压氢气区域．在迭代求解过程中，首先初步设定各边界温度分布，然后利用在边界面处得到的热通量求解相应的能量输运方程，得到新的温度分布．重复该迭代过程，直到达到期望的收敛精度．

1.6 可变物性

复合材料的热力性质由复合材料的组成和温度确定．复合材料发生热解反应时，其组成可认为是原始复合材料和热解炭产物的混合物[8]．假定忽略复合材料在火烧过程中的几何变形，剩余未发生热解反应的复合材料占初始原始材料的比例定义如下：

碳纤维/树脂复合材特性料具有各向异性的特点，其热物理性质在不同方向上差异较大[13]，由于储氢瓶在火灾条件下的传热主要发生在厚度方向上，本文假定具有各向同性的导热系数．表1为各区域的材料热物性参数．

表1 材料的热物性参数

Tab.1 Thermal properties of materials

2 问题描述

图1 计算区域及计算网格

流体区域初始温度293K，初始压力0.1MPa，其他区域的初始温度与流体区域设置相同．燃料进口设定为定质量流量边界，进气温度设定为庚烷沸点372K．计算区域的底面固体边界设为无滑移边界条件．流体区域的其他面设为自由边界条件，允许流体流入和流出．需要指出的是，在试验中使用了挡板进行防火，同时采用了强制通风，通风速率为10m3/s．然而在数值计算中，强制通风条件难以明确设置，因此本文采用自然通风条件．简化的自然通风条件会带来一定的预测误差．

储氢瓶外表面辐射发射率取为0.9，复合材料的热解属性[8]列入表2中．

表2 热解属性

Tab.2 Pyrolytic properties

3 结果与分析

3.1 火灾区域

火灾的精准预测对复合材料储氢瓶的传热特性研究具有重要影响．本研究使用的火灾预测模型以及火灾流体区域的网格分辨率在之前的研究中[5,10-11,17]已得到充分研究验证．图2为火灾区域的预测结果．

图2 火灾区域中截面预测结果

3.2 模型的验证

图3为试验中6个监测点处(如图1中所示)的试验温度与预测温度的对比图．数值模拟中忽略了液体燃料的气化过程，故火灾的预测温度相较于试验温度在初始阶段上升更迅速．由于数值计算采用自然通风条件，储氢瓶下方没有充足的空气实现燃料与空气的充分混合，罐下方为富燃燃烧条件，火灾的预测温度6略小于试验值．监测点1与5、2与4处于对称位置，其预测温度1与5、2与4几乎相同．而监测点1处的试验温度和预测温度差距较大，与对称的监测点5处试验温度同样相差较大，其原因在于强制通风条件带来的不确定性，以及试验过程中自然风对火焰的影响、碳烟沉积等不确定因素．目前试验方法存在着测试结果的可重复性问题，导致试验值与预测值存在一定差距，随着试验时间增加，两者差距逐渐减小．

储氢瓶的预测表面热通量与试验表面热通量的对比如图4所示．试验中[9]利用空气钢瓶测量表面热通量，根据钢瓶内的空气温度估算钢瓶平均表面热通量．由于钢的导热系数远大于复合材料的导热系数，试验中计算得到的表面热通量可能与实际的表面热通量有一定偏差．由图4可见，预测的平均表面热通量小于试验的平均表面热通量．这是由于试验中采用了强制通风条件，在钢瓶四周有充足的空气实现燃料与空气的充分混合，导致预测的储氢瓶底部热通量较低，试验储氢瓶表面热通量更均匀．预测最大表面热通量与试验平均表面热通量吻合较好，验证了此观点的正确性．

复合材料储氢瓶的失效往往发生在表面热通量最大处，试验得到的表面热通量为平均表面热通量，而最大表面热通量对储氢瓶的失效起着更为重要的作用．采用试验测得的平均表面热通量进行储氢瓶的失效分析，会低估储氢罐爆破或泄漏的风险．

（a）监测点1（b）监测点2（c）监测点3 （d）监测点4（e）监测点5（f）监测点6

图4 预测表面热通量与试验表面热通量的对比

3.3 外表面瞬时热通量及温度分布

图5为不同时刻下外表面的瞬时热通量和表面温度分布．图中左右视图分别是储氢瓶底部和顶部视图．表面热通量包括净辐射热通量和对流热通量，表现出显著的时间瞬态性和空间分布不均匀性．在自然通风条件下，储氢瓶下方的燃料在富燃条件下燃烧，导致火焰温度和碳烟浓度较低，因此储氢瓶底部表面的热通量明显低于顶部表面的热通量．在＝100s后，由于储氢瓶表面温度升高，表面热通量随之减小．与表面热通量相比，由于复合材料具有较大的热惯性，温度分布趋于更加均匀，并随时间的增加而增加．由于储氢瓶水平中分面周边具有较大的对流热通量，此处温度较高．

图5 外表面瞬时热通量和表面温度分布

3.4 最大及平均外表面热通量

图6为随时间变化的外表面最大和平均热通量．预测的最大表面热通量随时间减小，从初期约为200kW/m2的峰值逐渐降低至＝600s时75kW/m2左右．从图6可看出，在最大表面热通量中辐射传热热通量起主导作用，远大于对流传热热通量．平均表面热通量随时间变化的波动较小．由于表面热通量的非均匀分布，平均表面热通量明显低于最大表面热通量．在平均表面热通量中，虽然辐射换热仍大于对流换热，但是两者差距较小．

图6 最大和平均表面热通量

3.5 最大和平均界面温度

图7为随时间变化的最大和平均界面温度．复合材料缠绕层外表面的最高温度在初始阶段快速上升至约800K，然后温度升高的速度放缓，在＝600s时温度稳定在1140K左右．在＝120s之前，由于复合材料缠绕层的热惯性，内胆内外壁面的温度几乎保持不变，在＝200s后开始明显升高，在＝600s时达到540K，而HDPE的熔点在393～453K，致使瓶内高压氢气在储氢瓶爆破之前可能发生泄漏，试验中也证明了在低工作压力下会发生泄漏失效现象．

图7 随时间变化的最大和平均界面温度

3.6 瓶内氢气状态

预测的瓶内氢气平均温度和压力如图8所示.瓶内的氢气平均温度可以由导热方程直接求得，氢气的平均压力根据氢气的真实气体状态方程[18-19]计算.

式中：p为氢气压力；v为比容；n为与温度相关的系数，由文献[18-19]拟合数据得到，.氢气的平均温度在t＝300s之前几乎保持不变，然后在t＝600s时升至305K．氢气的压力变化趋势与温度变化趋势密切相关，t＝600s时压力略微增至72.8MPa．试验中的储氢瓶在t＝400s时发生爆破，在发生爆破前氢气压力几乎保持不变[9]，预测结果与试验结果相符．

研究普遍认为高压储氢瓶的爆破由外部热负荷与内部机械负荷(高压)共同作用．预测结果发现，在爆破前，储氢瓶内部机械负荷变化较小，表明储氢瓶并非因内部机械负荷急剧增大而发生爆破，其直接原因为外部热负荷产生的高温及热解条件引起的复合材料强度的降低．

3.7 复合材料热解

复合材料热解随温度和时间的变化分别如图9所示．当复合材料被加热至300℃左右时，即约在＝20s时，复合材料外表面开始发生热解反应，在＝311s时，剩余原始材料比例分数下降至热解反应最终值0.8，外表面完成热解反应．由于复合材料缠绕层的温度分布不均匀，导致在不同深度位置的升温速率和热解速率不同．为此，在复合材料表面温度最高处，选取不同深度位置的剩余原始材料比例数据进行分析．其中，较浅处位置升温速度快，温度较低时，热解速率较慢；而深处位置升温速度慢，更接近于试验中的升温速度，与试验数据[8]吻合更好．

图9 随温度和时间变化的剩余原始材料比例

复合材料的热解对数值计算结果的影响见图10.图(a)为考虑热解及不考虑热解两种情况下，监测点3处预测温度与试验温度对比．图(b)为考虑热解及不考虑热解两种情况下，储氢瓶的预测表面热通量与试验表面热通量对比．预测结果表明，本文提出的考虑复合材料热解的导热模型能够明显提高储氢瓶瞬态传热预测的准确性．

图10 热解对预测结果的影响

3.8 网格无关性验证

本研究中的池火火灾雷诺数较低，储氢瓶壁面边界层较厚，直接对边界层对流换热进行求解，因此网格必须具有较高的分辨率．数值模拟中将壁面边界层加密，计算得到Y Plus分布满足接近于1的要求．为了研究网格分辨率对数值计算的影响，使用粗网格进行验证．粗网格在外流体区域、复合材料缠绕层区域、内胆区域和高压氢区域的厚度方向上的分辨率为1/2．如图11所示，通过不同区域交界面处的热通量分析网格分辨率对不同区域数值计算的影响．由于最大热通量波动较大，两种网格分辨率下的最大热通量略有偏差，但变化的趋势几乎相同．平均热通量随时间波动较小，两种网格分辨率下的平均热通量预测值几乎重合．网格无关性验证结果表明，本文采用的网格分辨率满足对复合材料储氢罐的传热特性数值模拟要求．

（a）外表面热通量最大值（b）外表面热通量平均值 （c）中间界面及内表面热通量最大值（d）中间界面及内表面热通量平均值

4 结 论

(1) 采用LES对火灾进行了模拟，并对复合材料缠绕层、HDPE内胆区域和高压氢气区域分别建立了基于有限体积法的三维控制方程，并考虑了复合材料高温热解和热解气的对流和冷却作用，以及与温度和组分相关的材料物性．

(2)通过某商用Ⅳ型储氢瓶的池火试验验证了该模型的可靠性，并对储氢瓶的瞬态传热特性进行了研究．

(3)计算结果表明，在火灾条件下储氢瓶的表面热通量随时间波动较大，并且呈非均匀分布．采用平均表面热通量进行失效分析会低估发生失效分析的风险．在表面热通量中虽然辐射换热占主导地位，但是对流换热不可以忽略不计．

(4)研究发现，由于复合材料层具有较大的热惯性，储氢瓶内部氢气压力在失效之前压力升高不显著，发生失效的主要原因是由于高温以及热解引起的复合材料强度降级．

［1］ 赵永志，蒙 波，陈霖新，等. 氢能源的利用现状分析[J]. 化工进展，2015，34(9)：3248-3255.

Zhao Yongzhi，Meng Bo，Chen Linxin，et al. Utilization status of hydrogen energy[J].，2015，34(9)：3248-3255(in Chinese).

［2］ 杨文刚，李文斌，林 松，等. 碳纤维缠绕复合材料储氢气瓶的研制与应用进展[J]. 玻璃钢/复合材料，2015(12)：99-104.

Yang Wengang，Li Wenbin，Lin Song，et al. Research and application progess of carbon fiber composite hydrogen storage cylinder[J]./，2015(12)：99-104(in Chinese).

［3］ Xu B P，Hima L E，Wen J X，et al. Numerical study on the spontaneous ignition of pressurized hydrogen release through a tube into air [J].，2008，21(2)：205-213.

［4］ Ojha M，Dhiman A K，Guha K C. Simulation of thermally protected cylindrical container engulfed in fire [J].，2012，25(2)：391-399.

［5］ Saldi Z S，Wen J X. Modeling thermal response of polymer composite hydrogen cylinders subjected to external fires[J].，2017；42(11)：7513-7520.

［6］ 欧可升. 碳纤维全缠绕复合材料高压储氢气瓶耐局部火烧性能研究[D]. 杭州：浙江大学化学工程与生物工程学院，2014.

Ou Kesheng. Research on Fully-Wrapped Carbon Fiber Reinforced Composite High-Pressure Hydrogen Storage Cylinder Subjected to Localized Fire[D]. Hangzhou：School of Chemical and Biological Engineering，Zhejiang University，2014(in Chinese).

［7］ 张 颖. 碳纤维/环氧树脂复合材料热解动力学研究[D]. 沈阳：沈阳航空航天大学安全工程学院，2017.

Zhang Ying. Study on Pyrolysis Kinetics of Carbon Fiber/Epoxy Resin Composites[D]. Shenyang：School of Safety Engineering，Shenyang Aerospace University，2017(in Chinese).

［8］ Henderson J B，Wiebelt J A，Tant M R. A model for the thermal response of polymer composite materials with experimental verification[J].，1985，19(6)：579-595.

［9］ Ruban S，Heudier L，Jamois D，et al. Fire risk on high-pressure full composite cylinders for automotive applications[J].，2012，37(22)：17630-17638.

［10］ Menon S，Yeung P K，Kim W W. Effect of subgrid models on the computed interscale energy transfer in isotropic turbulence[J].，1996，25(2)：165-180.

［11］ Chen Zhibin，Wen J，Xu Baopeng，et al. Large eddy simulation of a medium-scale methanol pool fire using the extended eddy dissipation concept[J].，2014，70：389-408.

［12］ Dombrovsky L A，Dembele S，Wen J X，et al. Two-step method for radiative transfer calculations in a developing pool fire at the initial stage of its suppression by a water spray[J].，2018，127：717-26.

［13］ 沈蓉影. 碳纤维复合材料导热系数研究[J]. 材料工程，1993(3)：4-5.

Shen Rongying. Research on thermal conductivity of carbon fiber composites[J].，1993(3)：4-5(in Chinese).

［14］ Moroe S，Woodfield P L，Kimura K，et al. Measurements of hydrogen thermal conductivity at high pressure and high temperature[J].，2011，32(9)：1887-1917.

［15］ Woolley H W，Scott R B，Brickwedde F G. Compilation of thermal properties of hydrogen in its various isotopic and ortho-para modifications[J].，1948，41(5)：379-475.

［16］ International Organization for Standardization[ISO]. ISO/TS 15869 2009 Gaseous Hydrogen and Hydrogen Blends-Land Vehicle Fuel Tanks[S]. 2009.

［17］ Wang C J，Wen J X，Chen Z B，et al. Predicting radiative characteristics of hydrogen and hydrogen/methane jet fires using FireFOAM[J].，2014，39(35)：20560-20569.

［18］ 李 磊. 加氢站高压氢系统工艺参数研究[D]. 杭州：浙江大学材料与化学工程学院，2007.

Li Lei. Research on Technological Parameters of High Pressure System in Hydrogen Refueling Station[D]. Hangzhou：School of Materials and Chemical Engineering，Zhejiang University，2017(in Chinese).

［19］ Chen H G，Zheng J Y，Xu P，et al. Study on real-gas equations of high pressure hydrogen[J].，2010，35(7)：3100-3104.

Heat Transfer Behavior of Composite Hydrogen Cylinder Based on Conjugate Heat Transfer Subjected to Fire Impingement

Cheng Chonglü，Xu Baopeng，Kong Fanfu

(School of Energy and Power Engineering，Dalian University of Technology，Dalian 116024，China)

Composite hydrogen cylinder is a key component of hydrogen fuel cell vehicles，and the accurate prediction of its transient heat transfer behavior subjected to fire impingement is essential to its failure criterion．In this paper，a two-way coupled conjugate heat transfer model based on a multi-region approach is proposed for modelling the transient heat transfer behavior of composite hydrogen cylinders．The fire is modelled using large eddy simulation(LES)，and three-dimensional governing equations based on the finite volume method are established to model the heat transfer through the regions of composite laminate，HDPE liner and pressurized hydrogen，respectively．The effects of the gas decomposition of composite material and the volatile gas transport are taken into account，and the variable physical properties relying on composite material composition and temperature are adopted．The proposed conjugate heat transfer model is validated against a bonfire test of a commercial Type-4 cylinder，and its transient heat transfer characteristics are also studied in detail．

composite hydrogen cylinder；large eddy simulation；conjugate heat transfer；fire；decomposition

TK91

A

1006-8740(2020)01-0051-09

10.11715/rskxjs.R201812003

2018-12-03．

教育部留学归国人员科研启动基金资助项目(ZX20150008).

程崇律（1996—  ），男，博士研究生，chchlv9611@mail.dlut.edu.cn.

徐宝鹏，男，博士，教授，baopengxu@dlut.edu.cn.