賀正堯+陳文振+歐陽可漢

摘 要：目前對海洋環(huán)境中反應堆嚴重事故的研究較少，為給海上核應急提供參考，采用歐拉模型，針對海洋環(huán)境中放射性核素大氣擴散問題，基于GPU編寫了大氣擴散模擬程序，選取133Xe和133I兩種核素作為研究對象，分析核素在30公里范圍內(nèi)的擴散過程，并將計算數(shù)據(jù)實時顯示。結(jié)果表明，133I的積分濃度在下風方向約為133Xe的一半，干沉降是133I總量減少的主導因素?；贕PU的程序較非GPU版本加速約15%。

關(guān)鍵詞：反應堆事故；海洋大氣擴散；歐拉模型；放射性核素；GPU

中圖分類號：TL732 文獻標志碼：A 文章編號：2095-2945（2017）33-0001-04

引言

目前國內(nèi)外對位于沿?；騼?nèi)陸核電站的放射性物質(zhì)大氣擴散研究已經(jīng)比較成熟，建立了完善的應急處理機制。

近年，小型反應堆逐漸成為熱門話題，我國研發(fā)出了多種應用于海洋核動力平臺的堆型。但反應堆并不是絕對安全的，在嚴重事故的情況下，裂變產(chǎn)生的放射性核素有可能突破安全屏障大量進入到大氣中。

本文假想位于海上的小型反應堆發(fā)生嚴重事故，對事故中放射性核素的大氣擴散過程進行計算機模擬，期望為海上核應急提供參考。

為了能直觀準確地反映大氣擴散的詳細過程，本文選用了污染物大氣擴散模型中的歐拉模型?？紤]到數(shù)值求解歐拉模型時計算量大、并行度高的特點，本文將GPU通用計算技術(shù)應用到大氣擴散模擬研究中，一方面加速計算，降低計算成本；另一方面利于計算結(jié)果的實時顯示，數(shù)據(jù)的直觀性更強。

1 歐拉模型介紹

歐拉模型是一種網(wǎng)格模型，將需要研究的空間區(qū)域劃分成小的網(wǎng)格，應用數(shù)值方法離散化求解關(guān)于濃度的對流-擴散方程。

相較于工業(yè)中較常使用的高斯模型，歐拉模型沒有將煙羽視為一個整體，因此精度更高，能反映濃度隨時間變化的整個過程，而不是直接通過公式得到時間積分濃度。

1.1 對流-擴散方程

歐拉模型中的對流-擴散方程是根據(jù)質(zhì)量守恒定律和菲克定律導出的。這里定義在任意時刻、位置單位體積空氣中某種核素的活度為該核素的濃度，記為c（x，y，z，t），單位Bq/m3。暫時不考慮核素活度的變化，并將空氣視為不可壓縮流體，則在直角坐標系中有以下方程描述核素的濃度變化[1]：

（1）

式中S為該點的源強，單位Bq/（m3·s）；u、v、w分別代表風速在x、y、z方向上的分量向量；Kx、Ky、Kz分別是核素在三個方向上的擴散參數(shù)，單位m2/s。

在網(wǎng)格節(jié)點上求解對流-擴散方程時，不直接將（1）式離散化，而是先忽略z方向上的物質(zhì)交換，將Kx、Ky視為與橫縱坐標無關(guān)，在高度相同的每一層內(nèi)部做一次迭代：

（2）

之后計算第k層網(wǎng)格和上下兩層網(wǎng)格在豎直方向的物質(zhì)交換：

（3）

最底層和最上層網(wǎng)格的計算較為特殊，最底層網(wǎng)格需要考慮干沉降；最上層網(wǎng)格處于大氣混合層頂部，通常認為污染物不會再向上傳輸，因此這兩層網(wǎng)格節(jié)點的計算式為：

（2）～（5）式中下標i，j，k表示節(jié)點由東向西、由南向北、由低向高的序號，上標n表示經(jīng)過n個時間步長Δt之后的值；上標n'表示迭代過程的中間量；式（1）中的Kz由第k層與第k+1層之間的擴散參數(shù)Kz，k代替；vd為干沉降速度。

1.2 擴散參數(shù)的確定

x、y、z三個方向的擴散參數(shù)對于大氣擴散模型是至關(guān)重要的，選取合適的計算公式能提高模擬的精度。本文采用Seinfeld和Pandis提出的公式計算x、y方向的擴散參數(shù)：

（6）

當大氣穩(wěn)定度為近中性時，Kz，k可由下式計算：

（7）

其中w？鄢為對流速度尺度，h為混合層高度，z為兩層網(wǎng)格交界面的高度。

1.3 放射性核素的衰變

放射性核素在大氣擴散過程中，每時每刻都在發(fā)生衰變，其濃度會進一步降低，因此需要在計算大氣擴散的過程中加入對衰變量的計算。

對于衰變常數(shù)為λ的放射性核素，其在t時刻活度記為A（t），經(jīng)過Δt時間的衰變后，其活度變?yōu)椋?/p>

（8）

若該核素是另一核素衰變的子體，且母核素在嚴重事故中也被釋放到大氣中，情況就會變得復雜。設(shè)子核素的衰變常數(shù)為λ0，t時刻濃度為A0（t）；母核素衰變常數(shù)為λ1，t時刻濃度為A1（t）。若不考慮多級衰變，則母核素按照（8）式的規(guī)律衰變，而子核素經(jīng)過Δt時間后濃度為：

（9）

可見子核素的濃度計算需要用到上一時刻母核素的濃度，因此在計算時要同時跟蹤兩種核素。另外數(shù)值計算中需要注意精度的問題，為了保證收斂，時間步長不能取得太長，而較短的時間步長內(nèi)核素的衰變量及其微小，因此本文設(shè)定衰變的計算和擴散的計算不同步，在一種核素衰變了1%后計算一次衰變量，該時間間隔為：

（10）

2 基于GPU的模擬程序設(shè)計

GPU芯片中處理單元的數(shù)量在數(shù)百至數(shù)千量級，非常適合于并行計算，基于GPU的編程工具日漸成熟，GPU通用計算在工業(yè)領(lǐng)域的應用范圍也越來越廣。前文提到的歐拉模型中，求解對流-擴散方程的計算量很大，迭代式（2）～（5）的形式適合并行運算，因此本文基于GPU，采用CUDA語言編寫大氣擴散計算程序。

2.1 計算區(qū)域的設(shè)置

本文將坐標原點設(shè)置在海面上，放射性核素釋放點的正下方。坐標系x軸指向東，y軸指向北，z軸豎直向上。設(shè)定x、y、z三個方向的網(wǎng)格數(shù)量分別為127、127、20，則海面上的空間被劃分為322580個長方形體積單元，每一個體積的長寬均為500m，高度為h/20，控制點位于體積單元的中心，控制點坐標對應的核素濃度c代表該核素在體積單元內(nèi)的平均濃度。由于海上核設(shè)施的高度不高，反應堆功率不大，則釋放源位于最底層中心的體積單元內(nèi)，不考慮煙羽被抬升的情況。endprint

由于GPU不能直接讀取主機內(nèi)存中的數(shù)據(jù)，因此需要同時在主機的內(nèi)存和GPU顯存上劃分濃度數(shù)據(jù)的儲存空間，并在需要的時候進行數(shù)據(jù)交換，代碼如圖2。

2.2 并行線程的分配

從前文可知，計算過程中每一次迭代都是對眾多控制點上數(shù)值的一次更新，對于每一個控制點上的濃度數(shù)據(jù)來說，下一時刻的濃度只依賴于上一時刻的數(shù)據(jù)，因此可以將一個控制點的數(shù)值計算任務分配給GPU的一個線程（thread），讓線程按圖3的流程處理數(shù)據(jù)。

本文使用的GT720芯片可以同時執(zhí)行1024個線程，編程人員可以給GPU安排不同的線程組織形式，以適應不同的計算任務。本文程序中的代碼如下：

此段代碼安排的線程組織形式如圖5，其中每一個小格代表一個線程。每一個平面內(nèi)最后一行和最后一列線程閑置，因此每一個線程都和空間中的體積單元一一對應。

在分配完線程之后，只需在代碼中引用線程的編號來對相應控制點的數(shù)據(jù)進行讀取和寫入，代碼如圖6。

2.3 數(shù)據(jù)的可視化

由于GPU的基本工作就是處理圖像信息，因此使用GPU程序顯示顯存中的數(shù)據(jù)會非常方便快捷。但CUDA中沒有直接控制顯示的函數(shù)，而需要借助OpenGL綁定相應顯存區(qū)域顯示輸出，代碼如圖7。

3 放射性核素在海洋環(huán)境下的大氣擴散過程

本文在分析海洋環(huán)境特點的基礎(chǔ)上擬定了較為合適的大氣擴散條件，借助GPU程序?qū)?33Xe和133I兩種核素的大氣擴散過程進行分析，其中133Xe是133I的衰變子體，衰變常數(shù)分別為1.52×10-6s-1、9.26×10-6s-1。

3.1 海洋環(huán)境下大氣擴散條件分析

相較于陸地上或海陸間的污染物大氣擴散，海洋環(huán)境的大氣擴散有以下特點：

下墊面平坦。與陸地上地形的高低起伏不同，海平面的高度在沒有大浪的情況下變化量不大，因此研究海洋環(huán)境中的大氣擴散不需要考慮地形對風場的影響，也不需要考慮釋放源的建筑物尾流效應。

海表粗糙度長度遠低于陸地表面，通常陸地表面的粗糙度長度在0.01～1m的范圍內(nèi)，而海表為10-4～10-3m。這導致海洋大氣風廓線和陸地不同，見圖8（假設(shè)100m高空風速為10m/s），使得近地面污染物向下風向運動得更快。

海洋環(huán)境中干沉降的機制不同，影響因素更多。水面的物理特性與干燥的地面明顯不同，許多對水面干沉降的研究都描述了濕度和波浪對干沉降速度的影響[2]，認為水面附近的空氣濕度大導致氣溶膠粒徑增長，波浪的作用也會增加水體與氣溶膠的接觸，加速沉降。實測結(jié)果[3][4]也表明海洋環(huán)境中的干沉降速度要比陸地上的大數(shù)倍。本文算例中133I的干沉降速度設(shè)為0.06m/s，133Xe不發(fā)生干沉降。

本文根據(jù)海洋環(huán)境的特點擬定了一組氣象數(shù)據(jù)，大氣穩(wěn)定度為中性，沒有降水，較重要的數(shù)據(jù)列于表1內(nèi)。其中t為距離放射性物質(zhì)開始釋放的時間。設(shè)定事故進程中放射性物質(zhì)僅釋放一次，每種核素的釋放速率為1.0×1011Bq/s，釋放時間為7200秒，釋放高度為10m，事故進程中的大氣條件根據(jù)表中數(shù)據(jù)線性插值得到。

3.2 擴散過程分析

經(jīng)過基于GPU的程序計算，得到了133I和133Xe兩種核素的大氣擴散過程，見圖9；以及在下墊面的積分濃度分布，見圖10。時間步長為0.1s，不采用GPU加速的程序迭代速度為每步18.0毫秒，加速后為每步15.7毫秒，加速約15%。

從兩種核素的擴散過程來看，風向的改變使得原本往ESE方向運動的放射性煙團最終從SE方向離開計算區(qū)域，從開始釋放到放射性物質(zhì)的主體離開計算區(qū)域經(jīng)過了約4.2小時；反應堆釋放放射性物質(zhì)時，其附近的單種核素濃度超過了5×105Bq/m3，停止釋放后核素濃度迅速降低，很快被稀釋到103Bq/m3的水平，由于取了濃度的對數(shù)，對比兩種核素的瞬時濃度圖像很難發(fā)現(xiàn)其數(shù)值上的差別。另外，分析圖10中數(shù)量級為107的等值線，發(fā)現(xiàn)133I的積分濃度在干沉降和衰變的共同作用下，明顯低于133Xe的積分濃度，經(jīng)過具體計算，t=7200s時133I的干沉降量與衰變量分別占總減少量的76.1%和23.9%。

4 結(jié)束語

本文采用歐拉模型對海洋環(huán)境下放射性核素大氣擴散進行研究，借助GPU的并行計算能力將速度提高了15%，程序顯示的圖像能直觀實時地反映大氣擴散模擬過程。經(jīng)過計算分析，海洋環(huán)境下非惰性氣體核素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]宮釗，樊海燕.淺談放射性污染源及監(jiān)測方法[J].科技創(chuàng)新與應用，2012（02）：79.endprint