张希祥， 张 强， 张华兵， 程万洲， 魏然然，杨玉锋， 王 新， 高冬东， 李明阅

(1.中国石油天然气股份有限公司管道分公司管道科技研究中心油气管道输送安全国家工程实验室，河北廊坊065000；2.中国石油天然气股份有限公司管道分公司长沙输油气分公司，湖南长沙430000；3.中国石油天然气股份有限公司管道分公司山东中油天然气有限公司，山东济南250000)

1 风险评价方法简述

工业生产中的风险评价方法主要有定性风险评价法、半定量风险评价法[1]和定量风险评价法。定量风险评价通过对每个失效事件事故发生频率和后果严重程度进行风险分析[2-4]，从而得到具有单位的风险值来表征风险，比较有代表性的是荷兰应用科学院(TNO)提出的网格法[5]10。定量风险评价法一般依据历史失效数据进行失效频率调整，同时通过扩散、燃烧等过程的模拟进行后果分析，结果的精确度和准确性很高。

2 基于网格法的个人风险计算

文献[5]介绍了毒性气体扩散基于网格法的个人风险定量计算方法，该方法通过将目标区域划分为一个一个的网格，独立地计算各个网格的个人风险值，最终获得计算区域内各处的个人风险值。

2.1 网格的划分

采用笛卡尔直角坐标系，选取区域中的某点作为起始点对目标区域进行网格划分，划分的网格应该足够小，从而保证计算结果准确性，即个人风险在每个网格内没有特别大的变化。网格大小可以这样确定：

① 如果失效事件的影响距离小于等于300 m，选取的网格为25 m×25 m；

② 如果失效事件的影响距离大于300 m，选取的网格为100 m×100 m；

③ 根据实际情况，还可以采用混合网格，用小网格计算影响距离不超过300 m的区域，而对300 m以上区域用大网格进行计算。

确定网格大小的目的是为了保证网格内的属性基本一致，这样该网格中心点的个人风险值能充分代表该网格的个人风险值。通过计算每个网格的个人风险值，得到整个目标区域各网格的个人风险值，最终就可以获得该区域个人风险等值线图，以便进行下一步风险分析。

2.2 计算流程

① 选择失效事件i，如发生燃烧、爆炸、有毒物质泄漏等事件，进而对失效事件展开研究，如果仅考虑单个失效事件，则i=1。

② 选择气象等级M和一定角度范围的风向φ，确定联合分布概率。气象等级M包括大气稳定度和代表性风速，出现的概率为PM。气象等级数据可以参考气象行业的Rotterdam数据库或Pasquill数据库，也可使用蒙特卡洛法模拟某地区气象等级出现的概率。我国一般推荐使用Pasquill数据库，将大气稳定度等级分为强不稳定、不稳定、弱不稳定、中性、较稳定和稳定6级，分别用大写字母A、B、C、D、E、F表示。在每个大气稳定度等级后面加上代表性风速即代表一个气象等级，如C 3.0 m/s表示在弱不稳定条件下3 m/s风速的气象等级。该气象等级条件下某个角度范围的风向的概率为Pφ，φ为风向的编号，风向角度以正北为0°，顺时针旋转一周360°重合到0°结束，例如某一风向φ对应的风向角度为60°～90°，代表我们考虑60°～90°角度范围内的风向。不管使用经典气象数据还是使用蒙特卡洛法模拟，通过确定气象等级M和一定角度范围的风向φ，查询上文提到的数据库即可得到该地区气象等级和一定角度范围风向的出现概率，用PM×Pφ的值来表示。

③ 根据选择的失效事件，展开研究，选择引发事件S并查询其发生频率fS(单位为a-1)。例如阀门失效导致存放有毒气体的储罐发生泄漏，则阀门失效就是一个引发事件，可以从文献[5]78中查询各种引发事件的发生频率。如果仅考虑一个引发事件，则S=1。

④ 在确定失效事件、气象等级和风向、引发事件情况下，计算网格中心点的死亡概率Pd。

网格中心点的死亡概率用式(1)计算。

Pd=PclPci

(1)

式中Pd——该网格中心点的死亡概率

Pcl——该网格中心点处于有毒气体范围内导致人员死亡的概率

Pci——该网格中心点在有毒气体范围内概率

处于有毒气体范围内导致人员死亡的概率Pcl与死亡概率函数偏差值Pr相关，在文献[5]中有二者直接的对应表，根据Pr可确定Pcl，部分对应关系数据见表1[5]113。

表1 Pcl与死亡概率函数偏差值Pr的对应关系

死亡概率函数偏差值Pr通过式(2)确定[5]114。

Pr=a+bln(tρc)

(2)

式中Pr——死亡概率函数偏差值

a、b、c——描述不同物质毒性的常数

t——人员处于有毒气体中的时间，min

ρ——有毒气体质量浓度，mg/m3

有毒气体质量浓度ρ可通过高斯经典扩散模型获得[5]63，不同物质毒性的常数a、b和c可以通过表2查询，表2取自文献中数据[5]116。

表2 常见有毒物质的常数a、b和c的值

网格中心点处在有毒气体范围内的概率Pci的计算方法如下。

有毒气体的扩散是一个复杂的过程，本文采用一种简化的近似模型，即只考虑从泄漏点到网格中心点方向的扩散形式。以有毒气体泄漏点为坐标原点O，正北方向为y轴正向，正东方向为x轴正向，建立直角坐标系。本文不考虑有毒气体扩散高度，仅考虑平面扩散，结合文献[5]提到的近似模型，近似地认为有毒气体在水平面上扩散到网格中心点的形状为一个等腰三角形，即等腰三角形区域内为有毒气体扩散后覆盖的区域，见图1。其中，网格中心点C为底边的中点，坐标原点O到网格中心点C的直线为等腰三角形的高，则△OAB即为扩散区域。将等腰三角形底边AB的长度称为有毒气体宽度，用LE来表示，LE与气体扩散性质和当地气象条件相关。

网格中心点处于有毒气体范围内的概率Pci可按式(3)近似计算[5]237。

(3)

式中U——风向φ的总数

LE——有毒气体宽度，m

d——网格中心点距泄漏点距离，m

U表示在计算过程中考虑的一定角度范围风向的总数，例如，如果仅考虑1种角度范围的风向，则U=1；如果考虑2种角度范围的风向，则U=2，以此类推。

有毒气体宽度LE通过该处引起死亡的概率代表值PT来查询，PT表示在模拟工况下人员死亡的概率，模拟工况包括有毒气体的种类、有毒气体的浓度、气象等级、风速等条件。例如，PT=10%表示模拟引起10%的人员死亡。不同的气体由于扩散性质和对人体的毒性不同，其概率代表值PT也不尽相同，需要大量的模型函数和实际对照才可以准确地进行模拟。通常选择PT=10%进行计算。CO在气象等级为D 5.0 m/s，风向为196°～225°，CO质量浓度为21.3 g/m3条件下，有毒气体宽度LE与该处引起死亡的概率代表值PT的对应关系来自文献[5]149，见图2。

确定了Pcl和Pci的值后，通过式(1)确定网格中心点上的死亡概率Pd。

⑤ 计算网格中心点在失效事件i、气象等级M、风向φ和引发事件S单一条件下的个人风险值ΔIr。

ΔIr=PdfSPMPφ

(4)

式中 ΔIr——网格中心点单一条件下个人风险值，a-1

fS——引发事件S的发生频率，a-1

S——引发事件编号

PM——气象等级M出现的概率

M——气象等级编号

Pφ——风向φ出现的概率

φ——风向编号

⑥ 对于所有失效事件重复以上计算步骤。则该网格中心点总的个人风险值Ir计算如下：

(5)

式中Ir——该网格中心点总的个人风险值，a-1

J——失效事件总数

i——失效事件编号

K——气象等级总数

Q——引发事件总数

3 毒性气体泄漏个人风险计算实例

本文根据文献[5]提到的计算例子进行个人风险计算。某煤化工厂使用煤气作为燃料，煤气的主要成分是具有毒性的CO气体(我们将此处煤气近似为纯CO气体进行分析)。假设运输煤气的管道发生破裂引起煤气泄漏(全孔径破裂，泄漏质量流量为100 kg/s)[5]148，则失效事件为煤气泄漏，引发事件为管道破裂。设引发事件的发生频率fS为5×10-7a-1 [5]148，使用Rotterdam气象观测站的数据(气象等级为D 5.0 m/s条件，仅考虑出现1个角度范围的风向，即196°～225°角度范围的风向，PM和Pφ的乘积为0.036 8)[5]148。以煤气管道泄漏位置作为坐标原点，计算此次失效事件下，人员办公区域的网格中心点 (200 m,300 m)的个人风险值(CO扩散到该中心点有毒气体质量浓度为21.3 g/m3，计算使用的引起死亡的概率代表值PT为10 %，人员处于有毒气体的时间按30 min考虑)[5]148。

具体计算步骤如下：

① 选择失效事件，为管道发生煤气泄漏。

② 查询气象等级D代表性风速5.0 m/s，196°～225°角度范围的风向条件下，PM和Pφ的乘积为0.036 8。

③ 选择引发事件，为管道发生破裂，其发生频率为fS=5×10-7a-1。

④ 以泄漏点为坐标原点，正北方向为y轴正向，正东方向为x轴正向，建立笛卡尔坐标系。网格中心点(200 m，300 m)和原点(0，0)的距离d=361 m，网格中心点处有毒气体质量浓度ρ=21.3 g/m3。查表2，CO气体的毒性常数为a=-7.4,b=1,c=1。根据式(2)得，死亡概率函数偏差值Pr=5.97。根据概率偏差值Pr，查表1得Pcl=0.835。

根据图2，当PT=10%时，得LE=52 m。则由式(4)得有毒气体宽度LE=62.28 m。本例仅考虑出现1个风向，因此，U=1，根据式(3)计算，则网格中心点(200 m，300 m)在有毒气体范围内的概率Pci=0.023。

根据式(1)计算，则网格中心点(200 m，300 m)的死亡概率Pd=0.019。

⑤ 根据式(4)计算，则该失效事件对网格中心点(200 m，300 m)处给定条件下的个人风险值为ΔIr=3.5×10-10a-1。