基于GPU的放射性核素海洋大气扩散研究

贺正尧+陈文振+欧阳可汉

摘 要：目前对海洋环境中反应堆严重事故的研究较少，为给海上核应急提供参考，采用欧拉模型，针对海洋环境中放射性核素大气扩散问题，基于GPU编写了大气扩散模拟程序，选取133Xe和133I两种核素作为研究对象，分析核素在30公里范围内的扩散过程，并将计算数据实时显示。结果表明，133I的积分浓度在下风方向约为133Xe的一半，干沉降是133I总量减少的主导因素。基于GPU的程序较非GPU版本加速约15%。

关键词：反应堆事故；海洋大气扩散；欧拉模型；放射性核素；GPU

中图分类号：TL732 文献标志码：A 文章编号：2095-2945（2017）33-0001-04

引言

目前国内外对位于沿海或内陆核电站的放射性物质大气扩散研究已经比较成熟，建立了完善的应急处理机制。

近年，小型反应堆逐渐成为热门话题，我国研发出了多种应用于海洋核动力平台的堆型。但反应堆并不是绝对安全的，在严重事故的情况下，裂变产生的放射性核素有可能突破安全屏障大量进入到大气中。

本文假想位于海上的小型反应堆发生严重事故，对事故中放射性核素的大气扩散过程进行计算机模拟，期望为海上核应急提供参考。

为了能直观准确地反映大气扩散的详细过程，本文选用了污染物大气扩散模型中的欧拉模型。考虑到数值求解欧拉模型时计算量大、并行度高的特点，本文将GPU通用计算技术应用到大气扩散模拟研究中，一方面加速计算，降低计算成本；另一方面利于计算结果的实时显示，数据的直观性更强。

1 欧拉模型介绍

欧拉模型是一种网格模型，将需要研究的空间区域划分成小的网格，应用数值方法离散化求解关于浓度的对流-扩散方程。

相较于工业中较常使用的高斯模型，欧拉模型没有将烟羽视为一个整体，因此精度更高，能反映浓度随时间变化的整个过程，而不是直接通过公式得到时间积分浓度。

1.1 对流-扩散方程

欧拉模型中的对流-扩散方程是根据质量守恒定律和菲克定律导出的。这里定义在任意时刻、位置单位体积空气中某种核素的活度为该核素的浓度，记为c（x，y，z，t），单位Bq/m3。暂时不考虑核素活度的变化，并将空气视为不可压缩流体，则在直角坐标系中有以下方程描述核素的浓度变化[1]：

（1）

式中S为该点的源强，单位Bq/（m3·s）；u、v、w分别代表风速在x、y、z方向上的分量向量；Kx、Ky、Kz分别是核素在三个方向上的扩散参数，单位m2/s。

在网格节点上求解对流-扩散方程时，不直接将（1）式离散化，而是先忽略z方向上的物质交换，将Kx、Ky视为与横纵坐标无关，在高度相同的每一层内部做一次迭代：

（2）

之后计算第k层网格和上下两层网格在竖直方向的物质交换：

（3）

最底层和最上层网格的计算较为特殊，最底层网格需要考虑干沉降；最上层网格处于大气混合层顶部，通常认为污染物不会再向上传输，因此这两层网格节点的计算式为：

（2）～（5）式中下标i，j，k表示节点由东向西、由南向北、由低向高的序号，上标n表示经过n个时间步长Δt之后的值；上标n'表示迭代过程的中间量；式（1）中的Kz由第k层与第k+1层之间的扩散参数Kz，k代替；vd为干沉降速度。

1.2 扩散参数的确定

x、y、z三个方向的扩散参数对于大气扩散模型是至关重要的，选取合适的计算公式能提高模拟的精度。本文采用Seinfeld和Pandis提出的公式计算x、y方向的扩散参数：

（6）

当大气稳定度为近中性时，Kz，k可由下式计算：

（7）

其中w？鄢为对流速度尺度，h为混合层高度，z为两层网格交界面的高度。

1.3 放射性核素的衰变

放射性核素在大气扩散过程中，每时每刻都在发生衰变，其浓度会进一步降低，因此需要在计算大气扩散的过程中加入对衰变量的计算。

对于衰变常数为λ的放射性核素，其在t时刻活度记为A（t），经过Δt时间的衰变后，其活度变为：

（8）

若该核素是另一核素衰变的子体，且母核素在严重事故中也被释放到大气中，情况就会变得复杂。设子核素的衰变常数为λ0，t时刻浓度为A0（t）；母核素衰变常数为λ1，t时刻浓度为A1（t）。若不考虑多级衰变，则母核素按照（8）式的规律衰变，而子核素经过Δt时间后浓度为：

（9）

可見子核素的浓度计算需要用到上一时刻母核素的浓度，因此在计算时要同时跟踪两种核素。另外数值计算中需要注意精度的问题，为了保证收敛，时间步长不能取得太长，而较短的时间步长内核素的衰变量及其微小，因此本文设定衰变的计算和扩散的计算不同步，在一种核素衰变了1%后计算一次衰变量，该时间间隔为：

（10）

2 基于GPU的模拟程序设计

GPU芯片中处理单元的数量在数百至数千量级，非常适合于并行计算，基于GPU的编程工具日渐成熟，GPU通用计算在工业领域的应用范围也越来越广。前文提到的欧拉模型中，求解对流-扩散方程的计算量很大，迭代式（2）～（5）的形式适合并行运算，因此本文基于GPU，采用CUDA语言编写大气扩散计算程序。

2.1 计算区域的设置

本文将坐标原点设置在海面上，放射性核素释放点的正下方。坐标系x轴指向东，y轴指向北，z轴竖直向上。设定x、y、z三个方向的网格数量分别为127、127、20，则海面上的空间被划分为322580个长方形体积单元，每一个体积的长宽均为500m，高度为h/20，控制点位于体积单元的中心，控制点坐标对应的核素浓度c代表该核素在体积单元内的平均浓度。由于海上核设施的高度不高，反应堆功率不大，则释放源位于最底层中心的体积单元内，不考虑烟羽被抬升的情况。endprint

由于GPU不能直接读取主机内存中的数据，因此需要同时在主机的内存和GPU显存上划分浓度数据的储存空间，并在需要的时候进行数据交换，代码如图2。

2.2 并行线程的分配

从前文可知，计算过程中每一次迭代都是对众多控制点上数值的一次更新，对于每一个控制点上的浓度数据来说，下一时刻的浓度只依赖于上一时刻的数据，因此可以将一个控制点的数值计算任务分配给GPU的一个线程（thread），让线程按图3的流程处理数据。

本文使用的GT720芯片可以同时执行1024个线程，编程人员可以给GPU安排不同的线程组织形式，以适应不同的计算任务。本文程序中的代码如下：

此段代码安排的线程组织形式如图5，其中每一个小格代表一个线程。每一个平面内最后一行和最后一列线程闲置，因此每一个线程都和空间中的体积单元一一对应。

在分配完线程之后，只需在代码中引用线程的编号来对相应控制点的数据进行读取和写入，代码如图6。

2.3 数据的可视化

由于GPU的基本工作就是处理图像信息，因此使用GPU程序显示显存中的数据会非常方便快捷。但CUDA中没有直接控制显示的函数，而需要借助OpenGL绑定相应显存区域显示输出，代码如图7。

3 放射性核素在海洋环境下的大气扩散过程

本文在分析海洋环境特点的基础上拟定了较为合适的大气扩散条件，借助GPU程序对133Xe和133I两种核素的大气扩散过程进行分析，其中133Xe是133I的衰变子体，衰变常数分别为1.52×10-6s-1、9.26×10-6s-1。

3.1 海洋环境下大气扩散条件分析

相较于陆地上或海陆间的污染物大气扩散，海洋环境的大气扩散有以下特点：

下垫面平坦。与陆地上地形的高低起伏不同，海平面的高度在没有大浪的情况下变化量不大，因此研究海洋环境中的大气扩散不需要考虑地形对风场的影响，也不需要考虑释放源的建筑物尾流效应。

海表粗糙度长度远低于陆地表面，通常陆地表面的粗糙度长度在0.01～1m的范围内，而海表为10-4～10-3m。这导致海洋大气风廓线和陆地不同，见图8（假设100m高空风速为10m/s），使得近地面污染物向下风向运动得更快。

海洋环境中干沉降的机制不同，影响因素更多。水面的物理特性与干燥的地面明显不同，许多对水面干沉降的研究都描述了湿度和波浪对干沉降速度的影响[2]，认为水面附近的空气湿度大导致气溶胶粒径增长，波浪的作用也会增加水体与气溶胶的接触，加速沉降。实测结果[3][4]也表明海洋环境中的干沉降速度要比陆地上的大数倍。本文算例中133I的干沉降速度设为0.06m/s，133Xe不发生干沉降。

本文根据海洋环境的特点拟定了一组气象数据，大气稳定度为中性，没有降水，较重要的数据列于表1内。其中t为距离放射性物质开始释放的时间。设定事故进程中放射性物质仅释放一次，每种核素的释放速率为1.0×1011Bq/s，释放时间为7200秒，释放高度为10m，事故进程中的大气条件根据表中数据线性插值得到。

3.2 扩散过程分析

经过基于GPU的程序计算，得到了133I和133Xe两种核素的大气扩散过程，见图9；以及在下垫面的积分浓度分布，见图10。时间步长为0.1s，不采用GPU加速的程序迭代速度为每步18.0毫秒，加速后为每步15.7毫秒，加速约15%。

从两种核素的扩散过程来看，风向的改变使得原本往ESE方向运动的放射性煙团最终从SE方向离开计算区域，从开始释放到放射性物质的主体离开计算区域经过了约4.2小时；反应堆释放放射性物质时，其附近的单种核素浓度超过了5×105Bq/m3，停止释放后核素浓度迅速降低，很快被稀释到103Bq/m3的水平，由于取了浓度的对数，对比两种核素的瞬时浓度图像很难发现其数值上的差别。另外，分析图10中数量级为107的等值线，发现133I的积分浓度在干沉降和衰变的共同作用下，明显低于133Xe的积分浓度，经过具体计算，t=7200s时133I的干沉降量与衰变量分别占总减少量的76.1%和23.9%。

4 结束语

本文采用欧拉模型对海洋环境下放射性核素大气扩散进行研究，借助GPU的并行计算能力将速度提高了15%，程序显示的图像能直观实时地反映大气扩散模拟过程。经过计算分析，海洋环境下非惰性气体核素133I在下风方向的积分浓度约为惰性气体核素133Xe的一半，干沉降在133I总量减少的过程中起主导作用。

参考文献：

[1]Stockie J M. The Mathematics of Atmospheric Dispersion Modeling[J]. Society for Industrial and Applied Mathematics，2011，53（2）：349-372.

[2]Calec N， Boyer P， Anselmet F， et al. Dry deposition velocities of submicron aerosols on water surfaces： Laboratory experimental data and modelling approach[J]. Journal of Aerosol Science，2017，105：179-192.

[3]丁敏霞，苏玲玲，刘国卿，等.深圳市大气7Be沉降通量与气溶胶沉积速率[J].地球化学，2017，46（1）：81-86.

[4]詹建琼，陈立奇，张远辉，等.台湾海峡海表气溶胶干沉降通量研究[J].台湾海峡，2010，29（2）：258-264.

[5]宫钊，樊海燕.浅谈放射性污染源及监测方法[J].科技创新与应用，2012（02）：79.endprint