基于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