核电站的安全预测及分析
谷江源，黄 冲，晋超超
摘要： 本文就核电站的安全问题建立了相关模型，根据模型对核电站的情况进行分析，并给出合理的预
测。
对于问题一，查阅资料得1990-2011年我国发电量与用电量详情，得出用电量，发电量以及去核缺口量的变化趋势，用边际效应模型对问题进行分析，得出建立核电站的必要性的结论，以及最优的数量。 对于问题二，为了模拟核泄漏的扩散范围与路径,选用烟雾扩散模型进行分析并进行改进，验证得到一个预测核扩散影响范围的模型。
综合问题一二，对我国核电发展前景作出预测，得出我国核电快速发展短暂缓和，但在未来的几十年依旧会保持稳定增长趋势的结论。
关键词：边际效应模型，烟雾扩散模型
一.问题重述
2011年3月11日，日本遭受了9级大地震并引发了强烈的海啸。这次大地震及其引发的海啸不仅给日本以重创，而且由此造成的福岛核电站的核泄漏更是引起了全世界对核电站及其安全的重新思考。
请从互联网或报刊上搜集有关数据，根据这些数据建立评估核电站安全的数学模型。考虑： （1）随着人们生活水平的提高，用电量大幅增加，假设不建设核电站，用电量和发电量之间的差距有多大？建设一个某种规格的核电站能提供多少电力？建设核电站的经济成本和效益如何？目前国内有几个核电站、在建或准备建设的有几个？也就是要求建立建设核电站必要性的数学模型并分析；
（2）以秦山或大亚湾核电站为例（选一个），如果这些地方出现了严重的自然灾害造成了核泄漏（需要你自己作出合理假设），那么，在一定气象条件（一定风向、风力、下雨等）的情况下建立核扩散的数学模型，并讨论对周围多大范围的居民进行疏散以及其他的应对措施和可能的后果；
（3）综合前两问给出核电站的发展前景预测。
二．问题分析
经问题分析，研究核电站最优数量，如果最有数量等于或大于目前我国已有核电站数目+在建核电站数目+计划建设核电站数目，那么我们认为核电在有继续建设的需求；反之，则没有必要性。
建设核电站的必要性问题，对于经济效益，考虑的主要是由核电站的建设成本以及后期投资费用，运用边际效应模型进行分析，得出核电站建设的必要性结论以及我国核电发电量发展势头短暂减缓，但在未来的几十年依旧会保持稳定增长趋势的结论。
针对核泄露的扩散问题，以秦山核电站为例，应用高斯烟雾模型，首先考虑在无风，无障碍物，无反射情况下放射性物质辐射的情形，计算出空间中任一点任意时刻的放射性物质浓度。后对模型进行改进，考虑风速，并加入地面辐射以及降雨吸收因素，用matlab绘制出核扩散的三维图，模拟出核扩散的速度和范围，从而对核电站安全性进行评估。并讨论对周围多大范围的居民进行疏散以及其他的应对措施和可能的后果。
结合前两个模型得出我国核电站的发展前景会受到日本核泄漏的暂时影响，但大体仍然保持增
长趋势，几十年之后保持稳定。
三 模型假设
3.1边际效应模型的假设
1.假设用于发电的核燃料成本不变
2.把我国现运行的核电站年统计平均发电量作为标准核电站的发电量 3.假定安全运行是核电站的常规状态，既不考虑发生事故 3.2烟雾扩散模型的假设
1.假设核扩散可以看做气云在平整、无障碍物地面上空扩散，且气云中不发生化学反应，地面对气云无吸收
2 .假设核扩散不受温度影响
4 .假设核扩散服从扩散定律，即单位时间通过单位法向面积的流量与它的浓度梯度成正比 5. 将核电站泄露过程看作连续性泄露，即泄放时间相对于扩散时间一致的泄露 6. 假设风向为水平风向，风速和风向不随时间变化
四．符号说明
4.1 边际效应模型的符号说明
F：建设一座核电站前期所需投资的成本（元） L：建设一座核电站后期所需投资的成本（元） W: 一座核电站处理核废料所需投资的成本（元） C：每年核电站购入原材料的成本 N：现有核电站总数目
Qn：为一座核电站总发电量（千瓦时） Pe ：Pe为研究时间第一年的电价，
t： 考察时间
Td：投资于单个核电站促进技术进步的资金所带来的经济效益 EP:环保系数
r: 环境治理投资占GDP的比重 QC：全国每年的用电量（千瓦时） R: 所有核电站的总经济收益函数 S：真实GDP增长率 N:最优的核电站数量 d：t年内年均增长率
k：中国用于研究开发领域经费占国内生产总值的比重 Cn:核电站碳排放量 Cf:火电站碳排放量
4.2烟雾扩散模型符号说明 H：泄漏点O距地面的有效高度 T：任意扩散时刻
C（x,y,z,t）:空间任意一点的放射性物质浓度 K:核扩散系数 Ω：空间域
V:空间域其体积
S:一规则的球面面积
Q1:在(t,t+∆t)内通过Ω的流量
*
Q2: Ω内放射性物质的增量
Q:从泄漏源泄漏的放射性物质的总量 R:降雨对于放射性物质的吸收率
:单位强度的点源函数 （x,y,z）
ws:在风速影响下其扩散速度
wa;风向
Wd:风力影响范围每天减小量
五.模型建立与求解
5.1问题一的求解 5.1.1 核电站存在必要性
第一问中用电量和发电量之间的差距问题，我们选取中国全国为例，从《中国统计年鉴》中查的相应数据详细列出表1（详见附录）
图1.我国用电的生产量与消费量（1-9分别代表表1中1990-2010年）
由表1和图一可以看出，随着我国城市化、工业化推进，社会电气化水平不断提高，我国用电需求总量在直线攀升，在中长期内仍将会不断增加。目前供电量基本能满足用电量要求，且有电能输出，出口量逐年增加。但是如果除去核电的发电量，那么我国的缺口量大约是成直线上升趋势.趋势图如图2.因此发展核电势在必需。
图2.我国除去核电缺口量走势 （1-9分别代表表1中1990-2010年） 5.1.2效益分析
中国现有核电站在t年内的成本函数为：
C=(F+L+W+ct)⨯N
核电站在研究时间范围内通过发电产生的经济收益:
t-1
i=0
与此同时，引入判断核电站是否为绿色电站的指标环保系数Ep，如果Cn÷(Cn+Cf)>0.5,则Ep=0，反
Qn⨯N⨯Pe⨯[∑(1+d)]⨯(1-r)
i
-EP
之为1。
技术进步带来的经济收益:
Td
⨯N
假设国内对于发展核电站的需求是国内总发电量与所有核电站发电量得差值的货币表示:
t-1
i
t
( Qc- Qn⨯N)⨯Pe⨯[∑(1+d)]⨯(1+S)
i=0
技术进步带来的单位核电站经济产出等于技术进步的投资:
t-1
i
所有核电站的总经济收益函数R：
t-1
Td=k⨯( Qc- Qn⨯N)⨯Pe⨯[∑(1+d)]⨯(1+S)
i=0
t
R= Qn⨯Pe⨯[∑(1+d)]⨯(1-r)
i
i=0
-EP
⨯N+Td⨯N
为了计算的方便，对R做如下的变形：
t-1
令： x=Qn⨯Pe⨯[∑(1+d)]⨯(1-r)
i
i=0
-EP
t-1
+ k⨯Qc⨯Pe⨯[∑(1+d)]⨯(1+S)
i
i=0
t
t-1
i
y= k⨯Qn⨯Pe⨯[∑(1+d)]⨯(1+S)
i=0
则所有核电站的总经济收益函数R化简得：
t
成本函数对N求导得到边际成本MC： MC
R=-yN+xN
2
=F+L+W+ct
收益函数对N求导得到边际收益：
MR=-2yN+x
由微观经济学中的利润最大化原则知，当MR=MC时，利润达到最大化，即核电站得到了最大的利用并且带来了最大的利润。所以最优的核电站数量为N:
*
N=
*
x-(F+L+W+ct)
2y
5.1.3模型求解
该模型的求解是通过计算机实现，得出合理核电站数目，以得出必要的结论。通过查阅资料得到以下数据：
Qn=[1**********]（千瓦时），为单座核电厂年发电量 全国不同行业的电价不唯一，大多分布在0.2~0.8之间，如下表2：
表2
全国用电均价Pe为：0.65⨯84.3%+0.5⨯7.7%+0.5⨯6.8%+0.8⨯1.2%=0.63元/千瓦时 F=[1**********]（元人民币）建设一座核电站前期所需投资的成本（元） L=1400000000（元人民币）建设一座核电站后期所需投资的成本（元） W=3500000000（元人民币）一座核电站处理核废料所需投资成本（元） Qc=[1**********]00（千瓦时）全国年用电量
t=为考察年限，在此定30年。t年内年均增长率d=3%
S=核电站真实GDP增长率。真实GDP增长率每年略有不同，在此设为8% 用matlab编写的程序一，代入以上数据进行模型求解，具体程序见附录： 模型求解的图像如下：
图3
由图知，随着技术进步和用电需求量直线攀升，我国最优核电站数不断增加，在未来五年，最
优的核电站数最大，约为80座，这个数字远大于我国现运行以及在建和筹建的核电站数目总和。因此，我国有必要继续建设核电站。未来10年后最优核电站数目缓慢下降，最后稳定在60座左右。这是因为核电站的成本相对于火力和水电发电偏高，且核反应的主要材料铀235的产量减少，价格有明显的上涨趋势。
5.1.4考虑政府政策对模型的修正
政策对核电站发展的影响可通过改变我国投资于核电站研究开发资金占GDP的比重来实现，即
为：k1、k2、k3。通过matlab编写的程序二（参见附录）进行模拟，得出相应曲线走势图：
图4
由上图知，随着技术进步以及我国能源需求量不断增加，核电迅速发展，在5年后左右，最优的核电站的数目达到最大，约为100座。此时所有核电站的总社会经济效益达到最大。
图中在3-5年左右时，由于受2011年日本核泄露事故的重大影响，我国核电快速发展势头暂时减缓，2011年全年没有核电机组开工就是一具体表现。在未来的几年，由于我国对核电站的政策更加谨慎，对核电的投资开发相对减少，这会更加加剧我国能源紧张局面，而使得最优的核电站数目在之后几年相比无政策影响的情况有所增大，因此最优核电站数目最大值比无政策影响时要大。
在未来大约10年后，最优核电站数目有一段迅速下降，原因可能是国家对核电站的整顿，使得部分旧核电站退役等。
但随后最优核电站数目仍会慢慢减少，并趋于稳定，约为60座，因为成本的增加，以及其他新能源（如太阳能）的发展等，平衡了对核电的需求。 5.2问题二的求解
5.2.1准静态扩散模型（不考虑风力、风向、下雨等因素）
时刻t空间中任意一点（x,y,z）的放射性物质浓度记为C。单位时间通过单位法相面的流量为：
q=-k×gradC （1）
考察空间域Ω，Ω的体积为V，包围Ω的曲面为S，S的外法线向量为n,则在[t,t+Δt]内通过Ω的流量为： Q1=
⎰
t+∆t
→→
t
⎰⎰q∙ndσdt
s
（2）
而Ω内放射性物质的增量为： Q2=
⎰⎰⎰[C(x,y,z-H,t)-C(x,y,z-H,t+∆t)]dV
v
（3）
有质量守恒定律：
Q1= Q2 （4）
根据曲面面积分的奥氏公式：
→→→
⎰⎰q∙ndσ=
s
⎰⎰⎰divqdV
v
（5）
有（1）~（5）不难得到：
2
⎛∂2C∂2C∂C
=k⋅div(grandC)=k ++2 ∂x2
∂t∂y∂(z-H⎝
∂C
)2
⎫
⎪,t>0,-∞
这是无界区域的抛物型偏微分方程。根据假设，初始条件为作用在坐标原点的点源函数，可以记作：
C（x,y,z,0）=Qδ（x,y,z-H） (7) 方程（6）满足条件（7）的解为： C=
Q
e
-
x+y+(z-H
4kt
22
)2
(4πkt)
3
（8）
结果表明，对于任意时刻t，放射性物质浓度C的等值面是球面x2+y2+（z-H）2=R2,并且随着球面半径R的增加C值是连续减少的；当R—>∞或t—>∞时，C—>0 将x+y+(z-H)=R代入方程（8),
R=-4kt⋅ln
2
2
2
2
2
C(4πkt)Q
3
（9）
R为放射源总泄漏量Q的情况下的受自然扩散影响范围的半径.
通过matlab，可以实现该模型的图像生成。
特此假定：由于扩散半径很大，相对于R，H做忽略处理，即坐标轴定为地面 k=0.00001（参考了其他元素粒子的扩散系数均为此数量级） Q=1000000（发生七级核泄漏时的放射性物质泄漏量） C=100（根据国际放射委员会确定） t=1000000(秒)（扩散时间）
程序详见附录一，此次输入参数为：e1（0.00001,1000000,100,1000000,0,0,0,0）由matlab
生成图像为：
图5
由matlab得：球面半径R=8.85公里，即距离放射源8.85公里内均能受到污染。 5.2.2考虑到地面反射对模型的修正
考虑地面会扩散物质的反射作用，同时扩散的核素粒子受沉降等作用，设地面反射系数为α。这样进入大气的核辐射物质可以看成是两个部分：一是从泄漏源O直接扩散到空间的部分；二是从地面反射进入空间的部分。
于是（8）式所得模型可以修正为：
⎧H=h+∆h
⎧x2+y2⎫⎧⎪
exp⎨-Q⎨⎬⋅⎨exp
C(x,y,z)=4kt⎭⎩1.5⎩⎪()4πkt⎩
⎡(z-H)2⎤
⎢-⎥+α∙exp
4kt⎣⎦
⎡(z+H)2⎤⎫
⎢-⎥⎬
4kt⎣⎦⎭
5.2.3考虑风力和降雨因素对模型修正
降雨因素：假设降雨对于放射性物质的吸收率为r,则方程（9）变为：
R=-4kt⋅ln
2
C(4πkt)
3
(1-R)Q
（10）
风力因素：只考虑产生在辐射原点某一角度wa，某一强度ws的风力影响。 风力的影响范围从R开始每天减少wd，并且对于影响范围内每一点的作用相同。 用matlab模拟此扩散模型。
特此假定：k、Q、C、t参数同准静态扩散模型
r=0.5(假设雨水对于放射性物质的吸收率为50%)
ws=0.72(公里每小时)（考虑到广州大亚湾周边地形、建筑的影响，假定其在风速作用下其扩散速度近似为0.2米每秒）
wa=0.785（假设风向与x0y面平行，且与x轴成正45度角）
wd=0.5（假设风力影响范围从R开始，每天减弱500米） 程序详见附录一，输入参数为：e1(0.00001,1000000,100,1000000,0.5,0.72,0.785,0.5）由matlab生成图像为：
图6
由matlab得出数据：
图3的最大球的半径R=7.1184公里，最小球半径R0=3.1184公里，球心位移为L=200公里，放射性物质沿风向扩散最长距离为：7.1184+3.1184+200=210.2368（公里），其扩散面积近似为：
(2⋅R
+2⋅R0)⋅L/2+π⋅R/2+π⋅R0/2=2142.2(平方公里)
2
2
5.2.4模型分析与扩展：
此模型结果说明在给定条件下（风向、风速、下雨等），扩散的最远距离为210.2368公里，扩散面积近似为2142.2平方公里。其优点在于考虑了影响核扩散的几个主要因素，基于这几个因素得出一个合理化的建议。
对于核泄漏的预防，基于此模型考虑，核电站应考虑在背风点选址，建立相应应急系统，比如一旦检测到核泄漏，就能立即自动注水或者封堆，以减少核物质的泄露总量。
若发生核泄漏，优先考虑在扩散面积内的人员，及时采用人工降雨等措施减弱其扩散。 模型的缺点在于一些系数的选取如扩散系数k没有考虑风速的影响，还有忽略一些因素影响（如泄漏源的有效高度、放射性物质的自身沉降作用等） 模型的扩展：可用于一些有毒烟雾的扩散研究等。
六.核电站发展前景分析与预测
由一二问结果以及我国2011年全年没有核电机组开工的政策可看出日本核事故对我核电发展已经或产生重大影响，我国核电快速发展势头出现减缓态势，但我国核电发展预期前景依旧光明，经历过缓和阶段，在未来的几十年依旧会保持稳定增长趋势。原因如下：
从总量上看，随着我国城市化、工业化推进，社会电气化水平不断提高，我国电力需求总量在中长期内仍将直线攀升。2011 年我国人均电力装机为0.78 kW,与高能效发达国家日本（人均约1.25 kW）及西欧（人均1.3～1.5 kW)想比还相差甚远，我国发展核电的空间仍长期存在。
从能源消费结构上看，也必须大力发展核电等清洁能源。我国提出到2020 年实现非化石能源占比15% 的目标，及2020 年比2005 年单位GDP二氧化碳排放强度下降40～45%目标，发展核电是必须选向。
按我国核电项目表，到2015 年我国将建成核电机组41 台，装机容量达4179.3 万kW。即使今后两年内不再审批新的核电项目，到2017 年将建成核电机组47 台、装机容量达4827.3 万kW。若仅开展前期工作的12 台机组，到2020 年我国亦可建成核电机组59 台、装机容量达6220 万kW。核电机组数的变化趋势与核能发电量变化趋势一致，根据国家政策，即使采取最保守的估计，我国未来10年核电变化趋势依旧会保持稳定增长。 参考文献：
[1]姜启源，谢金星，叶俊.数学模型（第三版）[M].北京：高等教育出版社 [2]赵静，但琦.数学建模与数学实验（第三版）[M].北京：高等教育出版社 [3]中广核官网.www.cgnpc.com.cn/n2881959/index.html [4]国家电力信息网.http://www.sp.com.cn/
[5]中华人民共和国统计局官网.http://www.stats.gov.cn/ [6]肖新建.2011年中国核电发展状况_未来趋势及政策建议 2011
附录：
表1.1990-2010年我国发电量和用电量详情
单位：亿千瓦*时
程序一：
function []=e5(d,r,Qn,Pe,Ep,k,Qc,S,F,L,W,c) t=0:1:30;N=zeros(1,31); for i=1:31 a =0 ;
for t1=0:1:t(i) b=(1+d)^t1; a=a+b; end
x=Qn*Pe*a*((1-r)^(-Ep))+k*Qc*Pe*a*(1+S)^t1; y=k*Qn*Pe*a*(1+S)^t1;
N(i)=(x-(F+L+W+c*t1))/(2*y); m=N(i)
end
plot(t,N,'r')
legend('最优核电站数目与考察时间函数图')
%可用调用：e5(0.03,0.0149,3*10^10,0.63,1,0.01,[1**********]00,0.08,[1**********],140 0000000,3500000000,2*10^7)
程序二：
function []=e6 (d,r,Qn,Pe,Ep,k1,k2,k3,Qc,S,F,L,W,c) t=0:1:30;N=zeros(1,31);
for i=1:31
a=0;
for t1=0:1:t(i)
b=(1+d)^t1;
a=a+b;
end
if t1
k=k1;
elseif t1
k=k2;
else k=k3;
end
x=Qn*Pe*a*((1-r)^(-Ep))+k*Qc*Pe*a*(1+S)^t1;
y=k*Qn*Pe*a*(1+S)^t1;
N(i)=(x-(F+L+W+c*t1))/(2*y);
mm=N(i)
end
plot(t,N,'r')
legend('最优核电站数目与考察时间函数图')
%可调用：e6(0.03,0.0149,3*10^10,0.63,1,0.01,0.006,0.05,[1**********]00,0.08,2646800000 0,1400000000,3500000000,2*10^7)
程序三：
function R=e1(k,Q,C,t,r,ws,wa,wd)
R=sqrt((-4)*k*t*log((C/((1-r)*Q))*(4*pi*k*t)^(1.5)))
for i=0:3:t/3600
if (wd*(i/3))
R=R-wd;
R0=R;
t=0:0.1:0.5*pi;
p=0:0.1:2*pi;
wsx=ws*cos(wa)*i;
wsy=ws*sin(wa)*i;
[theta,phi]=meshgrid(t,p);
x=R*sin(theta).*cos(phi)+wsx;
y=R*sin(theta).*sin(phi)+wsy;
z=R*cos(theta);
hold on
surfc(x,y,z)
end
end
R0=R0
view(30,60)
title（'核泄漏没有风、雨的影响'）;
%title('核泄漏有风、雨的影响'); xlabel('x');
ylabel('y');
zlabel('z');
axis([-10,20,-10,20,0,10]) x0=ws*cos(wa)*1000000/3600; y0=ws*sin(wa)*1000000/3600; R1=sqrt((x0)^2+(y0)^2)
end
核电站的安全预测及分析
谷江源，黄 冲，晋超超
摘要： 本文就核电站的安全问题建立了相关模型，根据模型对核电站的情况进行分析，并给出合理的预
测。
对于问题一，查阅资料得1990-2011年我国发电量与用电量详情，得出用电量，发电量以及去核缺口量的变化趋势，用边际效应模型对问题进行分析，得出建立核电站的必要性的结论，以及最优的数量。 对于问题二，为了模拟核泄漏的扩散范围与路径,选用烟雾扩散模型进行分析并进行改进，验证得到一个预测核扩散影响范围的模型。
综合问题一二，对我国核电发展前景作出预测，得出我国核电快速发展短暂缓和，但在未来的几十年依旧会保持稳定增长趋势的结论。
关键词：边际效应模型，烟雾扩散模型
一.问题重述
2011年3月11日，日本遭受了9级大地震并引发了强烈的海啸。这次大地震及其引发的海啸不仅给日本以重创，而且由此造成的福岛核电站的核泄漏更是引起了全世界对核电站及其安全的重新思考。
请从互联网或报刊上搜集有关数据，根据这些数据建立评估核电站安全的数学模型。考虑： （1）随着人们生活水平的提高，用电量大幅增加，假设不建设核电站，用电量和发电量之间的差距有多大？建设一个某种规格的核电站能提供多少电力？建设核电站的经济成本和效益如何？目前国内有几个核电站、在建或准备建设的有几个？也就是要求建立建设核电站必要性的数学模型并分析；
（2）以秦山或大亚湾核电站为例（选一个），如果这些地方出现了严重的自然灾害造成了核泄漏（需要你自己作出合理假设），那么，在一定气象条件（一定风向、风力、下雨等）的情况下建立核扩散的数学模型，并讨论对周围多大范围的居民进行疏散以及其他的应对措施和可能的后果；
（3）综合前两问给出核电站的发展前景预测。
二．问题分析
经问题分析，研究核电站最优数量，如果最有数量等于或大于目前我国已有核电站数目+在建核电站数目+计划建设核电站数目，那么我们认为核电在有继续建设的需求；反之，则没有必要性。
建设核电站的必要性问题，对于经济效益，考虑的主要是由核电站的建设成本以及后期投资费用，运用边际效应模型进行分析，得出核电站建设的必要性结论以及我国核电发电量发展势头短暂减缓，但在未来的几十年依旧会保持稳定增长趋势的结论。
针对核泄露的扩散问题，以秦山核电站为例，应用高斯烟雾模型，首先考虑在无风，无障碍物，无反射情况下放射性物质辐射的情形，计算出空间中任一点任意时刻的放射性物质浓度。后对模型进行改进，考虑风速，并加入地面辐射以及降雨吸收因素，用matlab绘制出核扩散的三维图，模拟出核扩散的速度和范围，从而对核电站安全性进行评估。并讨论对周围多大范围的居民进行疏散以及其他的应对措施和可能的后果。
结合前两个模型得出我国核电站的发展前景会受到日本核泄漏的暂时影响，但大体仍然保持增
长趋势，几十年之后保持稳定。
三 模型假设
3.1边际效应模型的假设
1.假设用于发电的核燃料成本不变
2.把我国现运行的核电站年统计平均发电量作为标准核电站的发电量 3.假定安全运行是核电站的常规状态，既不考虑发生事故 3.2烟雾扩散模型的假设
1.假设核扩散可以看做气云在平整、无障碍物地面上空扩散，且气云中不发生化学反应，地面对气云无吸收
2 .假设核扩散不受温度影响
4 .假设核扩散服从扩散定律，即单位时间通过单位法向面积的流量与它的浓度梯度成正比 5. 将核电站泄露过程看作连续性泄露，即泄放时间相对于扩散时间一致的泄露 6. 假设风向为水平风向，风速和风向不随时间变化
四．符号说明
4.1 边际效应模型的符号说明
F：建设一座核电站前期所需投资的成本（元） L：建设一座核电站后期所需投资的成本（元） W: 一座核电站处理核废料所需投资的成本（元） C：每年核电站购入原材料的成本 N：现有核电站总数目
Qn：为一座核电站总发电量（千瓦时） Pe ：Pe为研究时间第一年的电价，
t： 考察时间
Td：投资于单个核电站促进技术进步的资金所带来的经济效益 EP:环保系数
r: 环境治理投资占GDP的比重 QC：全国每年的用电量（千瓦时） R: 所有核电站的总经济收益函数 S：真实GDP增长率 N:最优的核电站数量 d：t年内年均增长率
k：中国用于研究开发领域经费占国内生产总值的比重 Cn:核电站碳排放量 Cf:火电站碳排放量
4.2烟雾扩散模型符号说明 H：泄漏点O距地面的有效高度 T：任意扩散时刻
C（x,y,z,t）:空间任意一点的放射性物质浓度 K:核扩散系数 Ω：空间域
V:空间域其体积
S:一规则的球面面积
Q1:在(t,t+∆t)内通过Ω的流量
*
Q2: Ω内放射性物质的增量
Q:从泄漏源泄漏的放射性物质的总量 R:降雨对于放射性物质的吸收率
:单位强度的点源函数 （x,y,z）
ws:在风速影响下其扩散速度
wa;风向
Wd:风力影响范围每天减小量
五.模型建立与求解
5.1问题一的求解 5.1.1 核电站存在必要性
第一问中用电量和发电量之间的差距问题，我们选取中国全国为例，从《中国统计年鉴》中查的相应数据详细列出表1（详见附录）
图1.我国用电的生产量与消费量（1-9分别代表表1中1990-2010年）
由表1和图一可以看出，随着我国城市化、工业化推进，社会电气化水平不断提高，我国用电需求总量在直线攀升，在中长期内仍将会不断增加。目前供电量基本能满足用电量要求，且有电能输出，出口量逐年增加。但是如果除去核电的发电量，那么我国的缺口量大约是成直线上升趋势.趋势图如图2.因此发展核电势在必需。
图2.我国除去核电缺口量走势 （1-9分别代表表1中1990-2010年） 5.1.2效益分析
中国现有核电站在t年内的成本函数为：
C=(F+L+W+ct)⨯N
核电站在研究时间范围内通过发电产生的经济收益:
t-1
i=0
与此同时，引入判断核电站是否为绿色电站的指标环保系数Ep，如果Cn÷(Cn+Cf)>0.5,则Ep=0，反
Qn⨯N⨯Pe⨯[∑(1+d)]⨯(1-r)
i
-EP
之为1。
技术进步带来的经济收益:
Td
⨯N
假设国内对于发展核电站的需求是国内总发电量与所有核电站发电量得差值的货币表示:
t-1
i
t
( Qc- Qn⨯N)⨯Pe⨯[∑(1+d)]⨯(1+S)
i=0
技术进步带来的单位核电站经济产出等于技术进步的投资:
t-1
i
所有核电站的总经济收益函数R：
t-1
Td=k⨯( Qc- Qn⨯N)⨯Pe⨯[∑(1+d)]⨯(1+S)
i=0
t
R= Qn⨯Pe⨯[∑(1+d)]⨯(1-r)
i
i=0
-EP
⨯N+Td⨯N
为了计算的方便，对R做如下的变形：
t-1
令： x=Qn⨯Pe⨯[∑(1+d)]⨯(1-r)
i
i=0
-EP
t-1
+ k⨯Qc⨯Pe⨯[∑(1+d)]⨯(1+S)
i
i=0
t
t-1
i
y= k⨯Qn⨯Pe⨯[∑(1+d)]⨯(1+S)
i=0
则所有核电站的总经济收益函数R化简得：
t
成本函数对N求导得到边际成本MC： MC
R=-yN+xN
2
=F+L+W+ct
收益函数对N求导得到边际收益：
MR=-2yN+x
由微观经济学中的利润最大化原则知，当MR=MC时，利润达到最大化，即核电站得到了最大的利用并且带来了最大的利润。所以最优的核电站数量为N:
*
N=
*
x-(F+L+W+ct)
2y
5.1.3模型求解
该模型的求解是通过计算机实现，得出合理核电站数目，以得出必要的结论。通过查阅资料得到以下数据：
Qn=[1**********]（千瓦时），为单座核电厂年发电量 全国不同行业的电价不唯一，大多分布在0.2~0.8之间，如下表2：
表2
全国用电均价Pe为：0.65⨯84.3%+0.5⨯7.7%+0.5⨯6.8%+0.8⨯1.2%=0.63元/千瓦时 F=[1**********]（元人民币）建设一座核电站前期所需投资的成本（元） L=1400000000（元人民币）建设一座核电站后期所需投资的成本（元） W=3500000000（元人民币）一座核电站处理核废料所需投资成本（元） Qc=[1**********]00（千瓦时）全国年用电量
t=为考察年限，在此定30年。t年内年均增长率d=3%
S=核电站真实GDP增长率。真实GDP增长率每年略有不同，在此设为8% 用matlab编写的程序一，代入以上数据进行模型求解，具体程序见附录： 模型求解的图像如下：
图3
由图知，随着技术进步和用电需求量直线攀升，我国最优核电站数不断增加，在未来五年，最
优的核电站数最大，约为80座，这个数字远大于我国现运行以及在建和筹建的核电站数目总和。因此，我国有必要继续建设核电站。未来10年后最优核电站数目缓慢下降，最后稳定在60座左右。这是因为核电站的成本相对于火力和水电发电偏高，且核反应的主要材料铀235的产量减少，价格有明显的上涨趋势。
5.1.4考虑政府政策对模型的修正
政策对核电站发展的影响可通过改变我国投资于核电站研究开发资金占GDP的比重来实现，即
为：k1、k2、k3。通过matlab编写的程序二（参见附录）进行模拟，得出相应曲线走势图：
图4
由上图知，随着技术进步以及我国能源需求量不断增加，核电迅速发展，在5年后左右，最优的核电站的数目达到最大，约为100座。此时所有核电站的总社会经济效益达到最大。
图中在3-5年左右时，由于受2011年日本核泄露事故的重大影响，我国核电快速发展势头暂时减缓，2011年全年没有核电机组开工就是一具体表现。在未来的几年，由于我国对核电站的政策更加谨慎，对核电的投资开发相对减少，这会更加加剧我国能源紧张局面，而使得最优的核电站数目在之后几年相比无政策影响的情况有所增大，因此最优核电站数目最大值比无政策影响时要大。
在未来大约10年后，最优核电站数目有一段迅速下降，原因可能是国家对核电站的整顿，使得部分旧核电站退役等。
但随后最优核电站数目仍会慢慢减少，并趋于稳定，约为60座，因为成本的增加，以及其他新能源（如太阳能）的发展等，平衡了对核电的需求。 5.2问题二的求解
5.2.1准静态扩散模型（不考虑风力、风向、下雨等因素）
时刻t空间中任意一点（x,y,z）的放射性物质浓度记为C。单位时间通过单位法相面的流量为：
q=-k×gradC （1）
考察空间域Ω，Ω的体积为V，包围Ω的曲面为S，S的外法线向量为n,则在[t,t+Δt]内通过Ω的流量为： Q1=
⎰
t+∆t
→→
t
⎰⎰q∙ndσdt
s
（2）
而Ω内放射性物质的增量为： Q2=
⎰⎰⎰[C(x,y,z-H,t)-C(x,y,z-H,t+∆t)]dV
v
（3）
有质量守恒定律：
Q1= Q2 （4）
根据曲面面积分的奥氏公式：
→→→
⎰⎰q∙ndσ=
s
⎰⎰⎰divqdV
v
（5）
有（1）~（5）不难得到：
2
⎛∂2C∂2C∂C
=k⋅div(grandC)=k ++2 ∂x2
∂t∂y∂(z-H⎝
∂C
)2
⎫
⎪,t>0,-∞
这是无界区域的抛物型偏微分方程。根据假设，初始条件为作用在坐标原点的点源函数，可以记作：
C（x,y,z,0）=Qδ（x,y,z-H） (7) 方程（6）满足条件（7）的解为： C=
Q
e
-
x+y+(z-H
4kt
22
)2
(4πkt)
3
（8）
结果表明，对于任意时刻t，放射性物质浓度C的等值面是球面x2+y2+（z-H）2=R2,并且随着球面半径R的增加C值是连续减少的；当R—>∞或t—>∞时，C—>0 将x+y+(z-H)=R代入方程（8),
R=-4kt⋅ln
2
2
2
2
2
C(4πkt)Q
3
（9）
R为放射源总泄漏量Q的情况下的受自然扩散影响范围的半径.
通过matlab，可以实现该模型的图像生成。
特此假定：由于扩散半径很大，相对于R，H做忽略处理，即坐标轴定为地面 k=0.00001（参考了其他元素粒子的扩散系数均为此数量级） Q=1000000（发生七级核泄漏时的放射性物质泄漏量） C=100（根据国际放射委员会确定） t=1000000(秒)（扩散时间）
程序详见附录一，此次输入参数为：e1（0.00001,1000000,100,1000000,0,0,0,0）由matlab
生成图像为：
图5
由matlab得：球面半径R=8.85公里，即距离放射源8.85公里内均能受到污染。 5.2.2考虑到地面反射对模型的修正
考虑地面会扩散物质的反射作用，同时扩散的核素粒子受沉降等作用，设地面反射系数为α。这样进入大气的核辐射物质可以看成是两个部分：一是从泄漏源O直接扩散到空间的部分；二是从地面反射进入空间的部分。
于是（8）式所得模型可以修正为：
⎧H=h+∆h
⎧x2+y2⎫⎧⎪
exp⎨-Q⎨⎬⋅⎨exp
C(x,y,z)=4kt⎭⎩1.5⎩⎪()4πkt⎩
⎡(z-H)2⎤
⎢-⎥+α∙exp
4kt⎣⎦
⎡(z+H)2⎤⎫
⎢-⎥⎬
4kt⎣⎦⎭
5.2.3考虑风力和降雨因素对模型修正
降雨因素：假设降雨对于放射性物质的吸收率为r,则方程（9）变为：
R=-4kt⋅ln
2
C(4πkt)
3
(1-R)Q
（10）
风力因素：只考虑产生在辐射原点某一角度wa，某一强度ws的风力影响。 风力的影响范围从R开始每天减少wd，并且对于影响范围内每一点的作用相同。 用matlab模拟此扩散模型。
特此假定：k、Q、C、t参数同准静态扩散模型
r=0.5(假设雨水对于放射性物质的吸收率为50%)
ws=0.72(公里每小时)（考虑到广州大亚湾周边地形、建筑的影响，假定其在风速作用下其扩散速度近似为0.2米每秒）
wa=0.785（假设风向与x0y面平行，且与x轴成正45度角）
wd=0.5（假设风力影响范围从R开始，每天减弱500米） 程序详见附录一，输入参数为：e1(0.00001,1000000,100,1000000,0.5,0.72,0.785,0.5）由matlab生成图像为：
图6
由matlab得出数据：
图3的最大球的半径R=7.1184公里，最小球半径R0=3.1184公里，球心位移为L=200公里，放射性物质沿风向扩散最长距离为：7.1184+3.1184+200=210.2368（公里），其扩散面积近似为：
(2⋅R
+2⋅R0)⋅L/2+π⋅R/2+π⋅R0/2=2142.2(平方公里)
2
2
5.2.4模型分析与扩展：
此模型结果说明在给定条件下（风向、风速、下雨等），扩散的最远距离为210.2368公里，扩散面积近似为2142.2平方公里。其优点在于考虑了影响核扩散的几个主要因素，基于这几个因素得出一个合理化的建议。
对于核泄漏的预防，基于此模型考虑，核电站应考虑在背风点选址，建立相应应急系统，比如一旦检测到核泄漏，就能立即自动注水或者封堆，以减少核物质的泄露总量。
若发生核泄漏，优先考虑在扩散面积内的人员，及时采用人工降雨等措施减弱其扩散。 模型的缺点在于一些系数的选取如扩散系数k没有考虑风速的影响，还有忽略一些因素影响（如泄漏源的有效高度、放射性物质的自身沉降作用等） 模型的扩展：可用于一些有毒烟雾的扩散研究等。
六.核电站发展前景分析与预测
由一二问结果以及我国2011年全年没有核电机组开工的政策可看出日本核事故对我核电发展已经或产生重大影响，我国核电快速发展势头出现减缓态势，但我国核电发展预期前景依旧光明，经历过缓和阶段，在未来的几十年依旧会保持稳定增长趋势。原因如下：
从总量上看，随着我国城市化、工业化推进，社会电气化水平不断提高，我国电力需求总量在中长期内仍将直线攀升。2011 年我国人均电力装机为0.78 kW,与高能效发达国家日本（人均约1.25 kW）及西欧（人均1.3～1.5 kW)想比还相差甚远，我国发展核电的空间仍长期存在。
从能源消费结构上看，也必须大力发展核电等清洁能源。我国提出到2020 年实现非化石能源占比15% 的目标，及2020 年比2005 年单位GDP二氧化碳排放强度下降40～45%目标，发展核电是必须选向。
按我国核电项目表，到2015 年我国将建成核电机组41 台，装机容量达4179.3 万kW。即使今后两年内不再审批新的核电项目，到2017 年将建成核电机组47 台、装机容量达4827.3 万kW。若仅开展前期工作的12 台机组，到2020 年我国亦可建成核电机组59 台、装机容量达6220 万kW。核电机组数的变化趋势与核能发电量变化趋势一致，根据国家政策，即使采取最保守的估计，我国未来10年核电变化趋势依旧会保持稳定增长。 参考文献：
[1]姜启源，谢金星，叶俊.数学模型（第三版）[M].北京：高等教育出版社 [2]赵静，但琦.数学建模与数学实验（第三版）[M].北京：高等教育出版社 [3]中广核官网.www.cgnpc.com.cn/n2881959/index.html [4]国家电力信息网.http://www.sp.com.cn/
[5]中华人民共和国统计局官网.http://www.stats.gov.cn/ [6]肖新建.2011年中国核电发展状况_未来趋势及政策建议 2011
附录：
表1.1990-2010年我国发电量和用电量详情
单位：亿千瓦*时
程序一：
function []=e5(d,r,Qn,Pe,Ep,k,Qc,S,F,L,W,c) t=0:1:30;N=zeros(1,31); for i=1:31 a =0 ;
for t1=0:1:t(i) b=(1+d)^t1; a=a+b; end
x=Qn*Pe*a*((1-r)^(-Ep))+k*Qc*Pe*a*(1+S)^t1; y=k*Qn*Pe*a*(1+S)^t1;
N(i)=(x-(F+L+W+c*t1))/(2*y); m=N(i)
end
plot(t,N,'r')
legend('最优核电站数目与考察时间函数图')
%可用调用：e5(0.03,0.0149,3*10^10,0.63,1,0.01,[1**********]00,0.08,[1**********],140 0000000,3500000000,2*10^7)
程序二：
function []=e6 (d,r,Qn,Pe,Ep,k1,k2,k3,Qc,S,F,L,W,c) t=0:1:30;N=zeros(1,31);
for i=1:31
a=0;
for t1=0:1:t(i)
b=(1+d)^t1;
a=a+b;
end
if t1
k=k1;
elseif t1
k=k2;
else k=k3;
end
x=Qn*Pe*a*((1-r)^(-Ep))+k*Qc*Pe*a*(1+S)^t1;
y=k*Qn*Pe*a*(1+S)^t1;
N(i)=(x-(F+L+W+c*t1))/(2*y);
mm=N(i)
end
plot(t,N,'r')
legend('最优核电站数目与考察时间函数图')
%可调用：e6(0.03,0.0149,3*10^10,0.63,1,0.01,0.006,0.05,[1**********]00,0.08,2646800000 0,1400000000,3500000000,2*10^7)
程序三：
function R=e1(k,Q,C,t,r,ws,wa,wd)
R=sqrt((-4)*k*t*log((C/((1-r)*Q))*(4*pi*k*t)^(1.5)))
for i=0:3:t/3600
if (wd*(i/3))
R=R-wd;
R0=R;
t=0:0.1:0.5*pi;
p=0:0.1:2*pi;
wsx=ws*cos(wa)*i;
wsy=ws*sin(wa)*i;
[theta,phi]=meshgrid(t,p);
x=R*sin(theta).*cos(phi)+wsx;
y=R*sin(theta).*sin(phi)+wsy;
z=R*cos(theta);
hold on
surfc(x,y,z)
end
end
R0=R0
view(30,60)
title（'核泄漏没有风、雨的影响'）;
%title('核泄漏有风、雨的影响'); xlabel('x');
ylabel('y');
zlabel('z');
axis([-10,20,-10,20,0,10]) x0=ws*cos(wa)*1000000/3600; y0=ws*sin(wa)*1000000/3600; R1=sqrt((x0)^2+(y0)^2)
end