氢能作为一种清洁、高效、可持续的能量载体，已被广泛认为是未来能源发展的关键[1]。在现有的储氢方法中，高压储氢技术在长期储存方面最具优势[2]，中国储氢方式也正在朝着高压化和规模化方向发展[3]，目前较大规模的储氢系统包括绿电示范制氢项目的储氢罐区(约2 MPa)和加氢站的储氢设备(20~90 MPa)。一方面，氢原子具有最小的原子半径，很容易在金属中扩散，并进一步导致容器发生氢脆现象从而失去承载能力[4]；另一方面，氢气具有较低的点火能，极易点燃，且氢气燃烧极限较宽，一旦点燃可能造成严重的火灾或爆炸事故[5]。更重要的是，一旦某一储氢容器发生泄漏爆燃，可能导致邻近存储设备发生事故升级，造成比初始事故更为严重的多米诺事故。
为定量评估储氢设施泄漏的后果，研究者开发了不同的经验事故模型，包括高压氢气泄漏的虚拟喷嘴模型[6]、喷射火模型[7]、积分扩散模型[8]、气体状态参数模型[9]等。此外，2017年美国桑迪亚国家实验室开发了HyRAM平台[10]，将许多数据、模型和方法集成到氢能的定量风险分析中。但目前对储氢系统的多米诺事故定量风险分析的研究较浅，重点仍关注氢泄漏等一次事故的定量风险分析[11-13]。在传统的石油化工储罐区中，Bayes网络(BN)被证明是一种有效的多米诺事故分析方法。Khakzad等[14]将BN应用于多米诺事故的研究，并基于BN的后验概率推理确定了罐区多米诺效应最有可能的传播路径。Wu等[15]提出了一种基于领结图分析和BN的动态风险评估方法，且提出了从领结图到BN的映射算法来描述系统的动态性和不确定性，每当给出新的证据信息时，所有事件的发生概率都会立即动态更新。伍东等[16]根据多种事故模型分析了火电厂氢气储罐的火灾、爆炸事故后果，指出氢气物理爆炸和化学爆炸对邻近设施造成损坏要大于喷射火事故。
文[17]提出了一种基于BN的复杂风速风向与事件树下危险化学品储罐区多米诺事故分析方法。为完善当前储氢系统多米诺事故的定量风险评价，本文进一步提出一种储氢系统多米诺事故传播BN分析方法。首先，方法基于氢设施泄漏的事件树模型和专用后果评估模型等，获取系统内潜在的多米诺风险概率数据。然后，采用BN建模算法自动构建多个BN模型，并由BN模型计算储氢系统内包括多米诺风险的综合风险概率分布，详细分析多米诺事故传播模式，指导储氢设备的安全生产。
1 方法及模型
1.1 研究方法
本文提出的高压储氢系统的多米诺事故BN分析方法流程如图 1所示。首先构建涉氢装置泄漏的事件树模型，结合氢装置泄漏特定后果评估模型和储氢设备的基础数据等，完成一次多米诺事故风险定量计算，为BN建模提供可信数据；然后根据计算数据，以每一装置为根节点(事故源)完成BN结构、参数建模；最后通过BN模型进行先验概率计算和后验概率推理，分析系统内存储压强与综合风险间的关系、事故传播模式与泄漏场景。
1.2 氢设施泄漏事件树分析
1.2.1 事件树拓扑结构
1.2.2 事件树节点频率
表 1 涉氢装置泄漏场景频率 |
| 泄漏场景 | 氢气球罐 | 氢气瓶 | 压缩机 | 管道 |
| 小孔泄漏 | 3.5×10-5 | 2.1×10-6 | 2.0×10-1 | 1.3×10-5 |
| 中孔泄漏 | 7.1×10-6 | 6.6×10-7 | 7.9×10-3 | 1.8×10-6 |
| 大孔泄漏 | 4.3×10-7 | 3.8×10-7 | 2.1×10-4 | 9.1×10-7 |
| 完全破裂 | 9.0×10-7 | 2.1×10-7 | 3.0×10-5 | 6.4×10-7 |
注：实际管道泄漏场景频率为表格结果与管道长度之积 |
1.3 多米诺风险定量计算
表 2 事件树分支计算所用模型 |
| 事故类型 | 物理效应 | 所用模型 |
| JF | 热辐射 | 虚拟喷嘴模型、加权多点源喷射火模型、气体状态参数模型 |
| VCE | 超压 | 虚拟喷嘴模型、积分扩散模型、气体状态参数模型、TNO模型 |
| FF | — | — |
| FB | 热辐射 | 火球模型 |
| CGEE | 超压 | 热力学有效性—虚拟距离TNT模型 |
在每一事故源的事件树的每一叶节点下，调用上述模型，计算各事故源的场景(状态)频率、事故持续时间，受影响装置的升级向量、升级概率和失效时间。鉴于压力装置不同于常压浮顶罐，可以不考虑闪火造成事故升级的可能性。喷射火方向具有不确定性，本文考虑正交的4个等可能方向，代表所有方向的喷射火。扩散VCE的射流方向假设与风向一致，同时本文考虑了16种可能的风速—风向组合及其概率分布[17]，综合评估VCE事故风险。具体计算公式如下：
1) 事故源的场景定义为设备经历某种孔径泄漏最终发生某种事故类型的过程，事故源场景频率为
其中：fc, s, u为第u个事故源的第s个场景的概率，fl为某种孔径泄漏的发生概率，fmid, s, u为第u个事故源的第s个场景所在事件树分支的其他所有节点的概率之积。
2) 事故源事故持续时间为
其中：acc为事件树叶节点事故类型，tml是根据气体状态参数计算模型得到的最大泄漏时间，tfb是根据火球模型求得的火球持续时间。
3) 考虑到喷射火和扩散蒸气云爆炸的后果影响均随时间变化，因此将tml离散化，计算各时刻受影响装置的升级向量、升级概率、失效时间，随后取tml内的平均热辐射强度作为受影响装置的热辐射通量，取tml内的最大超压值作为受影响装置的超压值，最后根据平均热辐射值或最大超压值计算受影响装置的升级概率Pd和失效时间ttf。若“超压小于阈值”或“平均热辐射下的ttf>tte或ttf>tml”，则Pd为0。然而，实际储罐泄漏(尤其是小孔泄漏)的tml往往较大，离散后的计算量仍然巨大，因此可参考最大应急响应时间等现场信息，取有限的时间段进行离散化处理。
其中：P为受影响装置的瞬时升级概率，γ为概率单位变量，U为后果模型计算得到的物理效应(超压或热辐射)，V为受影响装置的体积，Pt为0~t时间段的点燃概率，Ppresent为蒸气云经过时点火源存在的概率，ω为点火效率。式(4)—(6)表示升级概率Probit计算模型。
4) 在某个风向下，第u个事故源的第s个场景对受影响装置的损伤风险为
其中fw为第w个风速—风向组合的发生概率。
1.4 Bayes网络
1.4.1 BN结构
BN结构建模流程如图 3所示。首先，确定初始事故源(根节点)；其次，筛选升级向量大于阈值且事故持续时间大于失效时间的受影响装置；再次，由于未失效装置可能受到已失效装置(父节点)的多种事故的协同影响，在所有满足升级条件的装置中，根据所有协同事故下的“平均失效时间最短原则”确定优先失效的受影响装置(子节点)，父、子节点间通过有向边连接，同一级子节点间不存在有向边，并根据已失效节点的协同效应重新计算剩余装置的升级向量、升级概率和失效时间；然后，判断BN网络层级是否超过预设的多米诺级数；最后，构建逻辑与、或门节点即确定逻辑门节点与上述装置节点间的有向边。
子装置节点m面临的协同效应最多考虑以下4种情况：
其中：节点m每一父节点的事件树分支场景可构成Descartes积形式的协同效应组合，cartj, m表示节点m面临的第j个协同效应集合；SF表示中安全场景，FFS表示最终发生闪火的场景，EPS表示最终发生爆炸的场景，k1为安全场景的总数，k2为闪火场景总数，k3为爆炸场景的总数，Pmax, j, mEPS和Pmax, j, mFFS分别表示为cartj, m中爆炸和闪火场景对节点m造成的最大升级概率，opj, m为Pmax, j, mEPS对应的超压值，thm为节点m的超压阈值，tf, j, mr为cartj, m组合中第f个父节点除闪火外的某种火灾场景对节点m的热辐射值，Cj, m为cartj, m组合中除闪火外的所有火灾场景总数。
case1和case2分别表示cartj, m中存在多个安全和闪火场景；case3表示cartj, m中存在多个爆炸场景，且爆炸超压大于节点m的阈值和超压升级概率大于闪火；case4表示cartj, m中存在除闪火外的多个火灾场景，且不存在爆炸场景或爆炸场景不会造成节点m发生事故升级。
由此，协同效应cartj, m下节点m的升级向量Uj, m、失效时间ttfj, m和升级概率Pd, j, m分别为：
其中：F1表示根据式(6)计算失效时间的函数，F2表示根据式(4)—(5)计算升级概率的函数。
1.4.2 Bayes网络参数
对于装置节点，其状态认为是装置编号与事件树各分支通过散列算法得到的散列值，并考虑安全状态的散列值为safety；对于逻辑门节点，其状态的散列值仅包括occurrence、safety，即事故发生与未发生。在任意事故源作为根节点所构建的BN结构中，各类节点的条件概率表(CPT)计算式为：
1) 根节点的CPT：
其中S为根节点的事故状态散列值总数。
2) 第m个子装置节点的CPT：
其中：q为节点m事故状态散列值的总数，sfrs, m为节点m的第s个事故状态的先验概率与该节点所有事故状态先验概率的总和之比，n为节点m面临的cartj, m的总数，若cartj, m中各事故的协同时间大于节点m的ttfj, m，则Pd, j, m由式(14)计算，反之为0。
3) 第x个逻辑或门的CPT：
其中ox为第x个或门节点面临的cartj, x的总数。
4) 第y个逻辑与门的CPT：
其中ay为第y个与门节点面临的cartj, y的总数。
2 分析与讨论
2.1 分析案例
本文具体考虑下述2个储氢系统场景。
案例1：参考新疆库车绿氢示范项目的储罐数量及布局，建立一个常温、气态存储的罐区模型，如图 4所示。罐区共有10个氢气球罐，球罐体积均为2 000 m3，存储温度等于环境温度298.15 K。以N指向S方向为0°方向，顺时针角度增加。
对上述2个案例，考虑16种可能的风速—风向组合及其概率分布，采用1.3及1.4节所述的风险计算方法和BN建模算法，自动构建各自失效设备作为根节点的BN模型。构建结果显示，描述案例1、2多米诺事故的BN模型数量分别为160和224个。随后基于BN模型开展综合风险分析、因果和诊断推理。
2.2 方法验证
以存储压强为10 MPa储氢罐区为例，考虑1 h的事故发展阶段，并将其离散化为10个时间点，且最大构建至2级多米诺事故节点(后同)。在180°风向下假设T2储罐作为初始事故源，构建的BN结构如图 6所示，其中L1、L2为或门节点，DL1、DL2为与门节点。特定风向下的多米诺风险概率为ΣuΣsfr, ws, u，分别利用本文方法和解析法[33]计算目标装置的多米诺风险概率，2种方法的计算结果如表 3所示。其中，E1~E5分别表示如下场景的散列值：“T2中孔泄漏发生90°JF” “T2中孔泄漏发生VCE” “T2完全破裂发生FB” “T2完全破裂发生CGEE” “T2中孔泄漏发生270°JF，T1、T3、T7安全”。可以看出，给定事件树中的多种事故类型证据E1~E5，2种方法的分析结果基本相同，验证了本文方法的准确性。
2.3 综合风险分析
综合风险为各设备的累计自失效风险$ \sum\limits_{u} \sum\limits_{s} f_{c, s, u}$ 与多米诺风险概率$ \sum\limits_{w} \sum\limits_{u} \sum\limits_{s} f_{r, w}^{s, u}$ 之和。文[32]规定了不同加氢站存储系统的工作压强范围，而当前暂没有相关文件指导储氢罐区的工作压强设置，因此将考虑固定压强的加氢站综合风险分析和储氢罐区的压强—综合风险敏感性分析。
2.3.1 加氢站综合风险分析
考虑案例2中的储氢瓶作为受影响设备，图 7描述了站内综合风险来源。可见：站内的综合风险主要来源于压缩机的自失效风险，概率约2.2×10-3 a-1；其次来源于各储氢瓶多米诺风险，概率约3.3×10-4 a-1；然后是管道的自失效风险，概率约8.6×10-5 a-1；而各储氢瓶的自失效风险最低。
由上述分析知，应重点关注压缩机和管道泄漏事件，减少高温下储氢系统的持续运行时间，严格控制站内火源，降低发生JF和VCE事故的概率。
2.3.2 储氢罐区敏感性分析
考虑高压化与规模化储氢的趋势，本节旨在分析案例1中不同存储压强下风险概率和各一次事故造成的多米诺风险概率变化规律。图 10说明1~15 MPa下储氢罐区综合风险概率随存储压强增大而增加；储存压强在1~2.5 MPa增大时，多米诺风险概率快速增加；当以2 MPa及以上的压强存储时，多米诺风险概率占比在25%以上，最高可达37%；而典型2 MPa存储压强的储氢罐区多米诺风险概率占比约35%，可见多米诺风险是储氢罐区重要风险来源；当存储压强在1~1.7 MPa或10~15 MPa时，多米诺风险概率占比均低于30%，且累计自失效风险概率都维持在10-5 a-1量级。
图 11说明了各一次事故造成的多米诺风险概率随存储压强的变化关系。可见爆炸事故(CGEE和VCE)造成的多米诺风险概率对存储压强的变化更加敏感，而火灾(JF和FB)造成的多米诺风险概率随存储压强变化不明显，JF在存储压强为10 MPa以上时可能造成多米诺事故，而FB的持续时间小于各设备的失效时间，其造成的多米诺风险概率约为0。
鉴于高压储氢容器存在氢脆现象导致泄漏和破裂，应重点关注设备连接部位，且严格控制火源，以降低泄漏、破裂和点燃的风险。
2.4 因果推理
固定储氢罐区存储压强为10 MPa，T8作为初始事故源；考虑加氢站内压缩机为高自失效风险设备，以C1为初始事故源。假设环境风向为180 °，构建的BN如图 12所示。
给定证据T8中孔泄漏发生VCE时，表 4描述了各节点的状态升级概率：一级节点T3的总升级概率略大于T7和T9，而二级节点T2、T4的总升级概率相同。最可能的传播模式为：首先T3或T7或T9发生中孔泄漏JF，接着T2或T4也发生中孔泄漏JF。然而，各设备的升级概率均较大，表明一次事故中氢气泄漏VCE对邻近设备的冲击较大。
表 4 储氢罐区T8中孔泄漏VCE各节点的状态升级概率 |
| 节点 | JF | VCE | FB | CGEE | ||||||||
| 小孔泄漏 | 中孔泄漏 | 大孔泄漏 | 小孔泄漏 | 中孔泄漏 | 大孔泄漏 | 完全破裂 | 完全破裂 | |||||
| 一级节点 | T3 | 2.6×10-1 | 2.9×10-1 | 3.7×10-2 | 1.3×10-2 | 1.4×10-2 | 1.8×10-3 | 3.5×10-2 | 8.1×10-2 | |||
| T7 | 2.6×10-1 | 2.9×10-1 | 3.7×10-2 | 1.3×10-2 | 1.4×10-2 | 1.8×10-3 | 3.4×10-2 | 8.0×10-2 | ||||
| T9 | 2.6×10-1 | 2.9×10-1 | 3.7×10-2 | 1.3×10-2 | 1.4×10-2 | 1.8×10-3 | 3.4×10-2 | 8.0×10-2 | ||||
| 二级节点 | T2 | 2.8×10-1 | 3.0×10-1 | 3.9×10-2 | 1.4×10-2 | 1.5×10-2 | 1.9×10-3 | 3.7×10-2 | 8.5×10-2 | |||
| T4 | 2.8×10-1 | 3.0×10-1 | 3.9×10-2 | 1.4×10-2 | 1.5×10-2 | 1.9×10-3 | 3.7×10-2 | 8.5×10-2 | ||||
给定证据C1大孔泄漏发生180 °的JF时，表 5描述了各节点的状态升级概率：一级节点中T1、T3的升级概率相同，且所有节点的升级概率均较大，最可能事故传播模式为：首先T1或T3发生中孔泄漏JF，接着T2在C1、T1(或T3)的协同作用下也发生中孔泄漏JF。
表 5 加氢站C1大孔泄漏JF各节点的状态升级概率 |
| 节点 | JF | VCE | FB | CGEE | ||||||||
| 小孔泄漏 | 中孔泄漏 | 大孔泄漏 | 小孔泄漏 | 中孔泄漏 | 大孔泄漏 | 完全破裂 | 完全破裂 | |||||
| 一级节点 | T1 | 2.3×10-1 | 3.3×10-1 | 2.2×10-1 | 1.2×10-2 | 1.6×10-2 | 1.1×10-2 | 5.4×10-2 | 1.3×10-1 | |||
| T3 | 2.3×10-1 | 3.3×10-1 | 2.2×10-1 | 1.2×10-2 | 1.6×10-2 | 1.1×10-2 | 5.4×10-2 | 1.3×10-1 | ||||
| 二级节点 | T2 | 1.5×10-1 | 2.1×10-1 | 1.4×10-1 | 7.3×10-3 | 1.0×10-2 | 6.9×10-3 | 3.4×10-2 | 8.0×10-2 | |||
2.5 诊断推理
根据图 12的BN结构进行诊断推理，分别分析储氢罐区T8和加氢站C1作为根节点时的泄漏场景。
表 6 不同证据下的T8自失效泄漏概率推理结果 |
| 证据 | 小孔泄漏 | 中孔泄漏 | 大孔泄漏 | 完全破裂 |
| E6 | 0 | 1.5×10-1 | 1.9×10-2 | 8.3×10-1 |
| E7 | 0 | 1.4×10-1 | 1.8×10-2 | 8.4×10-1 |
注：E6和E7分别表示发生一级和两级多米诺事故的散列值 |
表 7 不同证据下的C1自失效泄漏概率推理结果 |
| 证据 | 小孔泄漏 | 中孔泄漏 | 大孔泄漏 | 完全破裂 |
| E8 | 0 | 2.7×10-1 | 5.5×10-1 | 1.9×10-1 |
| E9 | 0 | 3.1×10-1 | 5.2×10-2 | 1.6×10-1 |
注：E8和E9分别表示发生一级和两级多米诺事故散列值 |
由上述分析可知，针对不同的储氢系统，当多米诺事故发生时，储氢罐区更可能是储罐完全破裂或中孔泄漏场景，而加氢站压缩机更可能是大孔或中孔泄漏场景。为降低装置泄漏概率，对于储氢罐区，可选用高可靠性的密封技术、耐腐蚀性的材料，同时加强对储罐的检查与维护，确保其能够承受预定的压力和温度，预防破裂事故的发生。对于加氢站高泄漏风险的压缩机，在采用耐高温、高压、腐蚀的密封部件的同时，应定期检查各部件的磨损情况，预防因机械故障等原因导致的大孔和中孔泄漏。
3 结论
本文提出一种高压储氢系统的多米诺事故BN分析方法。该方法集成了氢泄漏事件树、专属风险评估等模型，为储氢系统构建了考虑多种泄漏场景和事故协同的BN模型，揭示了气态储氢罐区和加氢站的综合风险分布。研究表明：当前常见的2 MPa储氢罐区的多米诺风险概率占比约35%，物理和化学蒸汽云爆炸是造成多米诺事故的主要事故类型，且其后果随存储压强变化明显。70 MPa加注压力的加氢站综合风险主要来源于压缩机的自失效风险，且喷射火事故对储氢瓶发生多米诺事故升级产生可观的贡献。在上述储氢系统中，爆炸对多米诺事故具有重要贡献，因此可采取抑爆措施减小点火概率，以降低系统的综合风险。
下一步将采用动态BN，详尽分析储氢系统的多米诺事故演变过程。
