1. 引言
能源短缺已成为现在社会的一大难题,核能作为一种高效能源,带给了人们的巨大利益。但由于一些技术风险和管理纰漏,核泄漏问题屡屡发生。如日本福岛核电站,在发生特大地震时由于备用系统的不充分和急救措施的不及时,引起了核泄漏,对人类和自然环境都造成了极大的伤害,引起了人们对核问题的深思。
由于重大的突发性核泄漏紧急灾害事件具有爆发性、空间分布不连续性、对周边地形和气象条件的敏感性等特点,研究核事故所释放的物质的时空分布需要高度精确的技术。该技术的研究对于更好地减轻核辐射危害、保护环境和降低损失等方面有着极其重要的意义。
2. 核辐射粉尘扩散模型
2.1. 模型概况
国内有许多专家学者对于核辐射粉尘的扩散问题进行研究,如乔方利 [1] 等人对日本福岛核泄漏扩散路径进行模拟和预测,但其没有考虑粉尘的沉降问题;张斌才 [2] 等人用高斯烟羽模型来研究大气污染扩散,虽然可以快速模拟出扩散过程但模型的参数和实用性需要随着温度地形等因素不断修正,可靠性不高;黄文宏 [3] 等人用三维传递现象来模拟气体扩散,但模型的方法复杂计算困难,张振华 [4] 等用有限差分计算与模拟放射性气体在介质中扩散,但在进行隐式有限差分计算时,若所选的时间差分增量较大,计算的误差相对较大,误差不可忽略,导致计算和模拟不够精确,陈坤 [5] 等人研究了高斯烟羽模型的山区含硫天然气泄漏的扩散,由于在地势起伏的山地、丘壑等地面特征较复杂的地区,将地形因子引入计算模型来降低预测值与实际值的差距,才能使普适性增强,Denise Hertwig [6] 等建立的高斯模型仅考虑了下风速度的影响,影响因素考虑不全面。基于此,本文通过改进的高斯烟团模型来考察,不但充分考虑到沉降、风力、地势等条件,能够对放射线粉尘的扩散进行模拟,而且克服了高斯烟羽模型对流动条件的局限性使对粉尘的时空分布的模拟结果更加真实可靠。
2.2. 粉尘扩散模型
假设核泄漏速率恒定,且放射性粉尘在无风情况下向各个方向的扩散速度相等且扩散过程中不发生化学反应。令放射性气体发生泄漏的时刻
,在上述理想状态下粉尘扩散模型如下:
(1)
其中
为高斯误差函数;
为初始浓度;
为积分时间;D为扩散强度;r为扩散系数;t为时间。
上述模型仅是一个理想状态下的预测浓度模型,没有考虑实际天气条件,以及有效泄漏源等实际因素的影响,预测结果十分粗糙。
基于实际情况,粉尘扩散受外界因素影响,综合考虑这些因素,构建高斯烟团模型的粉尘扩散修正模型.假设泄漏源有效高度为H,取其地面投影为坐标原点,x轴指向风向,由高斯烟团模型 [4] 得到点源泄漏的浓度分布:
(2)
其中,Q为源强;c为污染物质量浓度;k为环境风速;
分别为用浓度标准差表示的
轴上的扩散参数;H为泄漏有效高度。
如考虑风速对辐射的影响,令风速
,假设风速风向在采取数据阶段恒定不变,那么
(3)
风速随高度升高的变化 [7]:
(4)
上述各式中,
为扩散物温度;
为环境温度,K;P为气压;
为离地面10 m处的平均速度;
为核反应堆泄漏点距地面的几何高度;
为核反应堆泄漏点高度处的平均风速;D为核反应堆泄漏处直径;
为粉尘扩散速度;m为风速廓线系数,其值与大气稳定度有关。
由Matlab计算得,当
时,
。
由放射性物质的衰变规律 [1] [3],其浓度随时间的变化公式为:
(5)
其中
为初始浓度;
衰变常数;t经过的时间。
记
为放射性核素的半衰期,有
(6)
其中
为衰变下的源强。
由于131I的半衰期为8.04天 [8] [9],上述方程为:
含碘放射性核素的干沉积速度 [10] 为
,其沉积轨迹如图1。
从而沉积量为:
(7)
模型(2)变化为:
(8)
其中
为反射系数,放射性核素一般取0.5 [11] [12]。
假设在t = 0时刻,在原点释放的一个气团将随风飘动,并因扩散不断胀大。得到经粒子衰变与重力沉降(干沉积)修正后的浓度变化模型:
Figure 1. Depositional trajectories of 131I
图1. 131I的沉积轨迹图
2.3. 实际观测位置辐射污染模型及数值求解
取核辐射源
,直径
,环境温度为
,温差
。
选取观测地点为距核电站62 km的福岛市,其合扩散速度为:
(9)
假设
,则
。由于假设了x轴与平行风向平行,因而可以得到
,继而,
。
记开始监测放射性的时刻为0时刻,考虑风速作用,由反射倾斜云模式 [13] [14],得出瞬时点源浓度变化方程
其中
表示水平和竖直方向;
分别表示合扩散速度在x轴和z轴的分速度。
取放射源强
[9],故实际观测位置辐射污染模型(数值)可表示如下:
由Matlab软件进行编程计算,得
(约204天后)。考虑到实际情况,核电站在福岛县内距福岛市并不远,因此受辐射强度大,此值可以接受。沿风向方向初始放射源强度为Q辐射污染扩散仿真图,如图2。
如考虑137Cs等放射性元素的衰变,则有:
(其中
均为常数)
Figure 2. The initial radiation source intensity along the direction of wind direction is the simulation diagram of Q radiation pollution diffusion, where (A) is the equipotential diagram and (B) is the diffusion diagram along the direction of wind direction
图2. 沿风向方向初始放射源强度为Q辐射污染扩散仿真图,其中(A)为等势图,(B)为沿风向方向扩散图
3. 基于层次分析法的处理措施优化模型
核泄漏事故不仅会造成经济损失,更会危害到附近居民的健康和产生较大的社会影响。对核电站抢修可最大程度的降低人群受放射污染的计量,保护人民群众的安全。人们对核污染的恐慌心理会造成股市波动,土地废弃等,从而直接或间接地对各个产业造成影响。因此有关部门应及时公布相关信息降低群众的恐慌。
为了在核事故中最大限度地减少辐射伤害和降低损失,在评估处理措施时考虑三个最重要的因素,即时间,距离,庇护场所。
综合上述因素,大致有以下4种基本方案:
方案一:以疏散人群为首要任务(尽快将核电站工人及附近居民疏散到安全地带);
方案二:以阻断污染源为首要任务(对核泄漏威胁进行全面客观分析,不惜一切代价防止核泄漏);
方案三:以及时准确的公布信息为首要任务(公布放射性物质变化真实情况,消除不必要的恐慌,也让避难的人们提前做好心理准备,以便紧急情况下的避难行动);
方案四:以切断放射性物质的传播途径及控制传染源为首要任务(对已受核污染的人群及牲畜进行去污处理,防止二次污染)。
通过以上的方案及准则,构造决策层次图如下图3。
Figure 3. Nuclear leakage incident management decision hierarchy diagram
图3. 核泄漏事故处理决策层次图
根据Saaty等人提出的1~9尺度 [7],分别确定各层对上层每一因素的成对比较阵以及方案。
子准则层对父准则层的属性判断:
健康影响
;经济代价
;社会影响
。
父准则层对目标层的属性判断:
子准则层对目标层的判断矩阵
,
,
,
,
,
经过Matlab编程,得到一致性指标 [9] 比例数CR如下表1:
Table 1. The proportion of consistency index of each parameter
表1. 各参数的一致性指标比例数
由上表可以看出矩阵O, A1, A3的CR值大于0.1,即不满足一致性,采用最大方向改进方法进行一致性调整后,得到方案对总目标的权重如下表2:
Table 2. Weight of each alternative to the target
表2. 各备选方案对目标的权重
由上表看出各方案对总目标的权重排序为B > D > C > A。
根据上述模型,得到的解决方案如下:
第一:为了将核辐射危害降到最低,将损失减到最小,阻断核污染源是当务之急;
第二:对于已受污染及可能受污染影响的人群,牲畜及食品作物等进行全面的核辐射检查及去污处理,并进行有效地控制,防止二次污染;
第三:要及时向公众公开核泄漏相关信息,以及当局和东电公司的核处理计划,尽量降低公众的恐慌以及避免因为焦虑等引起的愤怒情绪对社会稳定造成的影响;
第四:疏散人群也是不容忽视的,尽可能的将辐射区人群安置到相对安全的地方,将人们的经济损失及健康危害降到最低。
基金项目
国家自然科学基金(11761025,11901114),广东省教育厅青年创新人才类(2017KQNCX081),广州市科技创新一般项目(201904010010),中山大学广东省计算科学重点实验室开放课题基金资助(2018001),海南省研究生创新科研课题项目(Hys2019-59)。
附录
附录1:
Briggs扩散参数
附录2:层次分析法的权重比较表
1. 日本核泄漏事故处理决策 判断矩阵一致性比例:0.1701;对总目标的权重:1.0000
日本核泄漏事故
处理决策 健康影响 经济代价 社会影响 Wi
健康影响 1.0000 1.0000 0.5000 0.2130
经济代价 1.0000 1.0000 0.1429 0.1403
社会影响 2.0000 7.0000 1.0000 0.6467
2. 健康影响 判断矩阵一致性比例:0.0000;对总目标的权重:0.2130
健康影响 个人可避免计量 集体可避免计量 Wi
个人可避免计量 1.0000 8.0000 0.8889
集体可避免计量 0.1250 1.0000 0.1111
3. 经济代价 判断矩阵一致性比例:0.0000;对总目标的权重:0.1403
经济代价 直接费用 损失费用 Wi
直接费用 1.0000 9.0000 0.9000
损失费用 0.1111 1.0000 0.1000
4. 社会影响 判断矩阵一致性比例:0.0000;对总目标的权重:0.6467
社会影响 公众心理影响 社会稳定性 Wi
公众心理影响 1.0000 0.2500 0.2000
社会稳定性 4.0000 1.0000 0.8000
5. 个人可避免计量 判断矩阵一致性比例:0.0451;对总目标的权重:0.1893
个人可避免计量 备选方案A 备选方案B 备选方案C 备选方案D Wi
备选方案A 1.0000 1.0000 0.1429 1.0000 0.1039
备选方案B 1.0000 1.0000 0.1429 0.3333 0.0789
备选方案C 7.0000 7.0000 1.0000 5.0000 0.6685
备选方案D 1.0000 3.0000 0.2000 1.0000 0.1487
6. 集体可避免计量 判断矩阵一致性比例:0.0345;对总目标的权重:0.0237
集体可避免
计量
备选方案A 备选方案B 备选方案C 备选方案D Wi
备选方案A 1.0000 1.0000 0.1667 0.2000 0.0806
备选方案B 1.0000 1.0000 0.2500 0.2500 0.0943
备选方案C 6.0000 4.0000 1.0000 0.5000 0.3509
备选方案D 5.0000 4.0000 2.0000 1.0000 0.4742
7. 直接费用 判断矩阵一致性比例:0.0361;对总目标的权重:0.1262
直接费用备选方案A 备选方案B 备选方案C 备选方案D Wi
备选方案A 1.0000 1.0000 0.3333 2.0000 0.1910
备选方案B 1.0000 1.0000 0.2500 2.0000 0.1777
备选方案C 3.0000 4.0000 1.0000 3.0000 0.5177
备选方案D 0.5000 0.5000 0.3333 1.0000 0.1136
8. 损失费用 判断矩阵一致性比例:0.0449;对总目标的权重:0.0140
损失费用备选方案A 备选方案B 备选方案C 备选方案D Wi
备选方案A 1.0000 7.0000 1.0000 3.0000 0.4223
备选方案B 0.1429 1.0000 0.2500 0.2000 0.0574
备选方案C 1.0000 4.0000 1.0000 2.0000 0.3318
备选方案D 0.3333 5.0000 0.5000 1.0000 0.1885
9. 公众心理影响 判断矩阵一致性比例:0.0327;对总目标的权重:0.1293
公众心理
影响 备选方案A 备选方案B 备选方案C 备选方案D Wi
备选方案A 1.0000 0.3333 3.0000 1.0000 0.2091
备选方案B 3.0000 1.0000 4.0000 2.0000 0.4628
备选方案C 0.3333 0.2500 1.0000 0.2500 0.0794
备选方案D 1.0000 0.5000 4.0000 1.0000 0.2487
10. 社会稳定性 判断矩阵一致性比例:0.0327;对总目标的权重:0.5174
社会稳定性备选方案A 备选方案B 备选方案C 备选方案D Wi
备选方案A 1.0000 0.3333 3.0000 0.5000 0.1799
备选方案B 3.0000 1.0000 4.0000 1.0000 0.3981
备选方案C 0.3333 0.2500 1.0000 0.3333 0.0873
备选方案D 2.0000 1.0000 3.0000 1.0000 0.3347
附录3:
%%%%%%%%%通过Matlab模拟辐射污染扩散
Q=7.7*10^18;%输入强源
u=9.475;%输入风速
d=1;%迭代步长
Z=0.5;%放射性核素强度
x=20:d:100;%沿风向距离
y=-100:d:100;%垂直风向距离
by0=0.08*x.*(1+0.0001*x).^(-1/2);
bz0=0.06*x.*(1+0.0015*x).^(-1/2);
by=by0.*(1+0.38*Z);
fz=(2.53-0.13*log(x)).*(0.55+0.042*log(x)).^(-1).*Z.^(0.35-0.03*log(x));
bz=bz0.*fz;
tempy1=-y.*y./by./by./2;
tempy2=2.718282.^(tempy1);
c=Q/pi/u*((by.*bz).^(-1)).*tempy2
Cs=40;
contour(x,y,c,Cs);
shading interp;
colorbar;
grid off;
附录4:
%%%%%%%%%通过Matlab模拟辐射污染扩散2
Q=input('泄露源:');%%泄露源强度为7.7*10^18
d=input('计算机精度:');%%计算精度为1
Z0=input('地面粗糙长度:');%%%取0.5
H=input('源高度:');%%%源高度为20
t=100;
z=0;
by0=0.08.*x.*(1+0.0001*x).^(-1/2);
bz0=0.06.*x.*(1+0.0015*x).^(-1/2);
by=by0.*(1+0.38*Z0);
bx0=0.08.*x.*(1+0.0001*x).^(-1/2);
bx=bx0.*(1+0.38*Z0);
fz=(2.53-0.13*log(x)).*(0.55+0.042*log(x)).^(-1).*Z0.^(0.35-0.03*log(x));
bz=bz0.*fz;
n=Q/((4*pi*t).^3/2).*(bx.*bz.*by).^(1/2);
m=exp(-x.^2./(4.*bx.*t)-y.^2./(4.*by.*t)-(z-H).^2./(4.*bz.*t));
c=n.*m;
l=input('输入条数:');%%条数为20
%contour(x,y,c,l)
mesh(x,y,c)
附录5:层次分析法Matlab程序
function [ W,ahpResult] = ahp(C)
n=length(C);
ahpResult=cell(n,1);
for k=1:n
m_n=size(C{k},2);%k层的元素个数
ahpResult{k}=cell(m_k,1);
for kk=1:m_k%求k+1各元素对k层kk元素的成对比较矩阵的特征值和特征向量
%为存储第k+1层所有元素相对k层kk元素的权重预留出空间,长度应等于C{k}{2,kk}的长度
ahpResult{k}{1,kk}=zeros(length(C{k}{2,kk},1));
%将相应正互反矩阵属于最大特征值的特征向量归一化后赋给ahpResult{k}{1,kk}中相应位置,
%这些位置由逻辑术组C{k}{2,kk}决定
ahpResult{k}{1,kk}(C{k}{2,kk})=V(:,ind)/sum(V(:,ind));
ahpResult{k}{2,kk}=maxD;%正互反矩阵的最大特征值
nn=size(C{k}{1,kk},1);%C{k}{1,kk}的阶数
ahpResult{k}{3,kk}=(maxD-nn)/(nn-1)/RI(nn);%相应的一致性比例
end
end
W=ahpResult{1}{1,1};
for k=2:n
W=cat(2,ahpResult{k}{1,:})*W;%cat(2,ahpResult{k}{1,:})是把k+1层所有元素相对k层各个元素的权重向量横向排在一起生成的权重矩阵
end
NOTES
*通讯作者。