乏燃料公路運輸事故條件下放射性核素彌散模型研究

陳黎偉，陶龍龍，周百昌，汪 進，龍鵬程，王 芳

(1.中國科學院核能安全技術研究所 中子輸運理論與輻射安全重點實驗室，安徽 合肥230031；2.中國科學技術大學，安徽 合肥 230027；3.安徽省核應急專業技術支持中心，安徽 合肥 230031；4.中國人民解放軍96901部隊23分隊，北京 100095)

截至2018年9月，我國有44個核反應堆在運行，容量為40.6 kMW，13個在建，容量為14 kMW，并計劃增建36 kMW瓦的核反應堆，到2020年，中國將擁有58 kMW的容量，乏燃料累積量將達到近9 000 tHM[1, 2]。2017年，我國出臺的《核安全與放射性污染防治“十三五”規劃及2025遠景目標》文件強調要推進乏燃料貯存和處理能力建設，建立保障機制，優化運行管理。可見，乏燃料管理已成為我國核燃料循環技術發展的重點問題[2, 3]。

乏燃料管理有三種主要方式：直接處理、后處理和乏燃料出口。無論是哪種方式，乏燃料運輸在乏燃料的管理中都起著關鍵作用[4]。然而，乏燃料在運輸過程中可能會受到不可抗力因素的影響，如：災害天氣、車禍、人為誤操作等，導致乏燃料組件及運輸容器受損，從而發生放射性核素泄漏[5]。泄漏的放射性核素會從燃料組件到運輸容器內腔再釋放到外部環境，對環境造成放射性污染[6]。有統計顯示，美國在1974年至1993年間與放射性運輸容器表面污染相關事故共53起[7]；2013年，歐盟一輛裝載乏燃料容器的火車因脫軌而導致有放射性污染源泄漏，且有25萬人居住在距離事故發生地僅2.5公里的范圍內[7]；我國在1988年至1998年有13起與放射性材料運輸相關的事故，且造成了57人受照[8]。李越等人系統分析了國內外乏燃料運輸現狀及存在的問題，指出目前我國主要的運輸方式為公路運輸，并提出現階段應加強乏燃料運輸事故情景及后果評價方面研究[2]。相對于乏燃料存貯泄漏事故[9]，由于事故地點的不確定性，乏燃料運輸事故具有更大的輻射風險[10]，因此，事故條件下乏燃料公路運輸核素彌散后果評價技術對于當前我國核事故早期應急是必不可少的。

乏燃料運輸事故后果評價關鍵技術中，大氣彌散模型是計算個人劑量與集體劑量的基礎，美國基于高斯模型開發的RADTRAN軟件實現了放射性材料運輸輻射影響的計算和評價[11]；我國王任澤等人基于高斯模型開發了放射性物質運輸后果評價軟件CRAMTRA1.0，并進行了對比測試[10]，馬文娟等人利用CRAMTRA1.0軟件分析了在三種運輸方案下乏燃料從中國廣東某核電站到西北地區后處理廠的輻射影響[12]。然而，運輸事故過程中的放射性物質擴散與遷移具有室外復雜環境的特點，如沉降、植被效應、復雜地形與建筑等，尤其當泄漏源高度位于植被或建筑物高度以下時，其泄漏的核素受植被效應與建筑物壁面效應影響較大，而高斯模型對于此類情況下并不能很好的模擬[13]。

計算流體力學(CFD)方法是解決復雜條件下大氣污染物彌散分布的有效工具[14]。謝清芳等人利用CFD方法，分析了不同來流方向條件下植被對鈾尾礦庫下風向氡的擴散分布及規律的影響[15]。Ghasemian等人利用CFD方法模擬了植被屏障下的流場和交通排放污染物濃度，為城市規劃者提供了一些有效建議[16]。Wang CH等人基于CFD方法系統分析了樹木與道路兩邊建筑物高度差異及道路兩邊建筑間隔距離差異對交通排放大氣污染物擴散規律影響機制[17]。相比于交通污染源近地面排放，乏燃料公路運輸泄漏事故還具有泄漏源位置不確定性的特點，從而為放射性污染物擴散與遷移帶來了不確定性。

鳳麟核能團隊圍繞反應堆事故診斷、放射性核素環境影響、核應急與公共安全等迫切需求，依托中子輸運設計與安全評價軟件系統SuperMC[18]、可靠性與概率安全分析軟件系統RiskA[19]、核安全云NCloud等核心基礎，建成了數字社會環境下的虛擬核電站Virtual4DS[20]。本文通過修正放射性核素衰變項及動量源損失項，提出一種基于CFD方法的乏燃料運輸事故條件下的放射性核素彌散模型，并在Virtual4DS框架下，基于OpenFOAM 4.1軟件引擎進行求解器開發以實現事故場景模型的數值模擬，且在模擬中結合不同泄漏源高度及有無樹木影響的條件，實現了乏燃料公路運輸事故下放射性核素分布，同時對比分析了不同條件下的核素彌散分布規律。

1 模型描述

1.1 乏燃料運輸事故場景描述

根據2003年國家核安全局發布的《關于頒發大亞灣核電站乏燃料公路運輸裝運批準書的通知》文件內容[21]，從大亞灣核電站運往甘肅四○四廠乏燃料后處理廠，途經廣州、長沙、武漢、信陽、漯河、鄭州、洛陽、西安、咸陽時為高速路段，從咸陽、蘭州、武威、甘肅乏燃料后處理廠時為國道路段。圖1為谷歌衛星地圖中的信陽至漯河高速路段的乏燃料運輸路線場景圖，圖1(a)表示有樹條件下的運輸場景；圖1(b)表示無樹條件下運輸場景。此類場景中道路兩旁均有建筑民房、服務區，為居民密集區，一旦發生核泄漏事故，將對周邊環境和人身造成放射性危害。因此，本文基于類似場景作為研究樣地進行場景建模。

圖1 乏燃料公路運輸路線場景

假設一輛乏燃料運輸卡車從臨時儲存庫行駛到后處理廠，車上裝有用于儲存壓水堆(PWR)燃料組件的NAC-STC型乏燃料容器。由于運輸乏燃料金屬容器中使用的金屬墊圈受到長期的高溫影響，金屬墊片的殘余線性載荷和總回彈距離因軟金屬外層護套的蠕變變形而松弛，在受到撞擊或其他機械振動后，乏燃料容器密封條發生松動。因此，放射性核素從儲罐泄漏到環境中[22]，以氣載的形式在環境中擴散和遷移，并在植被和土壤上沉積，對周邊環境造成了輻射危害[23]。

當樹木高于建筑物時，其對峽谷氣流和污染物濃度影響最大[17]，本工作主要考慮樹木高于建筑物的情況。乏燃料運輸場景的幾何模型如圖2所示，假定計算域為64H(x方向)×70H(y方向)×8H(z方向)(其中H=18 m)，建筑物、樹木和道路位于計算域的中部，其長度為10H(180 m)。道路的寬度是2H(36 m)，從建筑物到車輛行駛道路的距離是2/3H(12 m)，樹木位于建筑物和車輛行駛道路之間，寬度是1/3H(6 m)，高度是3/2H(18 m)，從建筑物到樹木的距離是1/3H(6 m)。大氣氣流沿x正向垂直于道路方向經過。

為準確評估放射性核素對兩邊建筑的輻射危害，假設靠近來流處的建筑物且離釋放源最近的壁面為壁面一；靠近去流方向處的建筑物且離釋放源最近的壁面為壁面二。

圖2 乏燃料公路運輸事故幾何模型

1.2 數學模型

近地面氣流在地面、植被、建筑等復雜條件的影響下，較易形成湍流[24]。為準確預測乏燃料公路運輸事故中泄漏核素的分布狀態，本文引入標準k-ε湍流模型，綜合考慮核素衰變及動量源的損失，建立了基于CFD方法的修正模型，用于模擬事故條件下的放射性核素彌散過程。模型中主要包括兩個部分，即流場模擬和放射性核素濃度模擬。

1.2.1 流場模擬

本文假設計算區域內空氣流動為不可壓縮流動，且考慮重力損失，由于在處理植被影響下的大氣污染物擴散時，可將樹木區域視為多孔介質[16, 25]，以壓力損失系數近似氣流在通過植被時產生的壓力損失，可以通過增加動量源損失項來修正動量守恒方程，以模擬樹木作用下的流場分布。簡化后的質量守恒方程與動量守恒方程如下：

質量守恒方程：

(1)

式中：ui(m·s-1)為速度分量，i=1,2,3分別對應x,y,z三個方向。

動量守恒方程：

(2)

式中：t——時間，s；

ρ——流體的密度；

P——壓力，N·m-2；

μ,μi——分別為動力粘度和湍流粘度；

g——重力加速度，m·s-2；

Sui——動量源損失項;

λ——壓力損失系數,m-1。

1.2.2 放射性核素濃度模擬

在大氣環境中，放射性核素濃度方程可表示為：

(3)

式中：C——坐標點(x，y，z)處的每個時間步的濃度值，kBq·m-3;

ΓC(m2·s-1)——湍流擴散系數,m2/s;

SC——源項kBq·m-3·s-1。

在乏燃料公路運輸事故場景中，泄漏的放射性核素會因大氣流場的驅動，受樹冠效應、壁面效應、放射性衰變和沉積[26,27]的影響，逐漸向周圍環境擴散和遷移。因此，我們需要對濃度方程作沉積、核素衰變等修正，然而，碘在乏燃料公路運輸事故中主要以氣相形式遷移，且只有極小的粒子(平均粒徑在2 μm以內)從容器中釋放出來[28]。因此，對于事故早期應急(0.5 h～1 d)，放射性核素沉積量可以忽略不計。

其濃度修正方程可以表示為：

(4)

式中：λd——核素衰變因子;

T1/2——核素半衰期。

1.3 模型驗證

由于輻射危害，用實驗方法較難驗證乏燃料公路運輸過程中的核素彌散過程。本工作將采用Jeanjean APR等人工作[25]中模擬的數據作為驗證案例，該驗證模型由一個街道峽谷組成，樹木寬度與高度比例為1∶1，在樹木兩邊分別有一排10H(長)×H(寬)×H(高)的建筑物(H=18 m)，氣流從垂直于建筑物x軸正方向經過，在兩排建筑中間和模型樹的底部有四條線源發射器，用于釋放出氣體示蹤劑SF6，其釋放率為Qi。其他邊界條件及參數設置參見文獻[25]，最后將濃度值使用歸一化濃度公式進行表示：

(5)

圖3為Jeanjean APR等人在靠來流方向較近的建筑壁面的濃度模擬結果與本工作模擬結果的比較，其中圖3(a)表示為有樹條件下的空街峽谷墻壁歸一化濃度模擬，圖3(b)表示為無樹條件下的空街峽谷墻壁歸一化濃度模擬。從模擬的結果可以看出，連續釋放源在風場及建筑壁面的作用下，在建筑壁面中間產生蓄積，且由于樹木的作用，造成了示蹤氣體濃度的升高，可見兩項工作在模擬的濃度分布情況是一致的。雖然部分區域等值線存在一些較小差異，這主要是因為本文在網格處理上與Jeanjean APR等人可能有所不同，但總體濃度擴散趨勢與分布是一致的。

圖3 靠來流方向較近的建筑壁面的濃度模擬

2 數值模擬

2.1 計算方法

在OpenFOAM 4.1多孔介質求解器的基礎上添加了修正后的濃度方程，對求解器進行重新開發，以適合本工作的需求。

具體開發步驟如下：

(1)創建輸入變量場：在多孔介質求解器中創建三個標量場(濃度場C、壓力場p，湍流擴散系數alphat)和一個矢量場(速度場U)；

(2)添加濃度方程頭文件CEqn.H：頭文件中寫入本文中的濃度輸運方程；

(3)添加傳輸參數頭文件：readTransportProperties.H：添加求解過程中所需要的參數，如：層流Prandtl數，湍流Prandtl數，衰變因子，沉積因子等；

(4)與速度、壓力聯立求解。在執行文件中(.C文件)添加 #include “CEqn.H”，使其基于pimple算法進行壓力速度耦合求解；

(5)最后進行求解器編譯；

(6)在算例求解過程中，給定速度、壓力、湍流擴散參數，湍動能、湍流耗散率、源項等初始條件和邊界條件，利用topoSet命令對樹木多孔介質區域進行指定，先對大氣風場進行穩態求解，最后將求解后的穩態風場作為初始條件，使用自定義求解器進行迭代求解。

另外，輸運方程均采用有限體積法離散[14]，所有方程采用中心差分格式和迎風格式進行求解。

2.2 邊界條件

在乏燃料運輸中，假設天氣晴朗，環境溫度為20 ℃，且為等溫狀態。空氣密度應為1.2 kg·m-3。事故發生在場景模型正中心位置，泄漏引起了放射性核素131I的連續釋放，泄漏率為3.3 kBq·m-3·s-1，大約為0.72 μg·m-3·s-1質量流量，碘在空氣中的擴散系數約為0.8×10-5m2s-1[29]，來流沿x正方向的初始風速為2 m·s-1，壁面邊界采用Neumann邊界條件，且為固體壁面無滑移條件，流體通過樹木的壓力損失系數為200，湍流長度尺度(Λ)通過選擇來流直徑確定。根據公式(6)至公式(8)對雷諾數(Re)、初始湍流動能(k)和湍流耗散率(ε)進行初始值設定。

(6)

(7)

(8)

本工作主要計算泄漏事故后發生的1 h以內的放射性131I連續釋放，分別考慮泄漏源在2 m(位于建筑物高度以下)和22.5 m(位于建筑物高度與樹高度之間)的高度，以表示運輸車輛經過普通路面、橋或邊坡等不同路面高度的情形。由于空氣中的碘摩爾分數是約為4.36×10-11，因而不影響空氣混合物流動性能，為了向核應急決策者提供更多的實際信息，在結果討論中考慮工作場所的131I最大允許濃度為0.33 kBq m-3[9]。

2.3 網格無關性測試

圖4 在t=720 s時刻y=630 m平面在z=2 m高度131I從x=558 m至x=594 m區間濃度

2.4 分析與討論

為了獲取不同乏燃料公路運輸場景下的核素彌散分布規律，下面結合有無樹木影響及不同源高兩個部分，從大氣流場的分布、濃度空間分布與建筑壁面濃度時間變化等方面對乏燃料運輸事故條件下放射性核素131I彌散規律進行了對比分析。

圖5 y=630 m平面釋放源附近穩態風場分布

圖5表示y=630 m平面釋放源附近穩態風場分布，其中圖5(a)表示無樹條件下的風場分布；圖5(b)表示有樹條件下的風場分布(虛線區域內表示樹木)。從圖中可以看出，樹木對建筑物之間的湍流漩渦有明顯的減弱作用，在無樹條件下，建筑物之間產生了較大渦旋，相比之下，有樹環境下的建筑物之間則沒有大的渦旋作用。

圖6給出了z=2 m高的釋放源131I在y=630 m平面第3 600 s時的濃度分布情況，其中圖6(a)表示無樹條件下的濃度分布；圖6(b)表示有樹條件下的濃度分布。無樹條件下的放射性核素131I受湍流回流作用向來流方向逐漸擴散，對靠近來流方向的建筑物污染較大。然而在有樹條件下，因兩建筑物之間渦旋效應被減弱，131I向泄漏源四周擴散較為明顯，且在風場的作用下，擴散逐漸向去流方向偏移，對去流方向建筑物逐漸造成污染。可以看出，樹木對放射性核素擴散范圍起到了一定的抑制作用。

圖6 z=2 m高的釋放源131I 在y=630 m平面t=3 600ths的濃度分布

從圖7中可以看出，無樹條件下的壁面一的131I最大濃度和平均濃度都明顯比壁面二要高。但是有樹條件下的情況恰恰相反，壁面二的131I最大濃度和平均濃度都明顯比壁面一要高，這是因為在無樹條件下131I受到回旋渦流的影響，釋放的核素不斷向來流方向建筑壁面蓄積，從而導致壁面一的放射性核素濃度較高現象。通過模擬，我們可以得到第3 600 s在無樹條件下，壁面二的131I平均濃度為0.005 kBq·m-3，比壁面一的平均濃度(0.025 kBq·m-3)低5倍；有樹條件下，壁面二的131I平均濃度為0.008 9 kBq·m-3，比壁面一的平均濃度(0.003 2 kBq·m-3)高2.8倍。而且，在相同邊界條件下，壁面一在無樹木條件的131I平均濃度比有樹條件下高7.8倍。由此可見，建筑物之間的樹木對建筑物的放射性污染起到了一定的抑制作用，且要加強無樹條件下來流方向附近的建筑物壁面去污工作。

圖7z=2 m高的釋放源131I在建筑物壁面的濃度

Fig.7131I concentration of the walls of the building at thez=2 m leakage height

圖8給出了z=22.5 m高的釋放源131I在y=630 m平面第3 600 s時的濃度分布情況，其中圖8(a)表示無樹條件下的濃度分布；圖8(b)表示有樹條件下的濃度分布。在無樹條件下泄漏的放射性核素首先向去流方向擴散，對去流方向較近的建筑物造成了污染，然后部分核素被渦旋捕捉后逐漸向來流方向建筑物擴散。在有樹條件下，釋放源逐漸在兩樹之間較為均勻的向四周擴散。相比源高z=2 m時131I濃度分布情況，無樹條件下131I擴散趨勢呈明顯不同，這主要是因為建筑物之間的渦流對擴散造成的影響。

圖8 z=22.5 m高的釋放源131I 在y=630 m平面t=3 600th s的濃度分布

圖9 z=22.5 m高的釋放源131I在建筑物壁面的濃度

從圖7和圖9中可以看出，無論釋放源在何處，有樹條件下的核素濃度增長變化曲線相對平緩。與源高z=2m情況不同的是，當釋放源高位于z=22.5 m時，無樹條件下的壁面一的131I最大濃度和平均濃度都明顯比壁面二要低。而在有樹條件下的情況相對復雜，從第0 s到1 680 s，壁面二的131I最大濃度均比壁面一要高；從1 680 s后，壁面一的131I最大濃度均比壁面二要高；但是對于兩旁建筑物的壁面平均131I濃度變化值非常接近，呈線性增長趨勢。這是因為在無樹條件下131I主要直接受到去流風的影響，釋放的核素不斷向去流方向擴散，對壁面二建筑壁面直接造成了放射性污染；有樹條件下泄漏源周圍去流方向風得到了減弱，131I不斷向四周均勻擴散。通過模擬，我們可以得到第3 600 s時刻在無樹條件下，壁面二平均131I濃度為0.006 6 kBq·m-3，比壁面一平均131I濃度(0.001 5 kBq·m-3)高4.4倍；比有樹條件下的壁面二平均131I濃度(0.002 1 kBq·m-3)要高3倍。由此可見，當源高位于建筑物以上，樹以下時，無樹條件下去流方向附近建筑放射性污染較為嚴重。

3 結論

本文基于CFD的方法研究了乏燃料公路運輸事故條件下放射性核素彌散過程，而且我們可以通過此方法得到乏燃料運輸事故后不同泄漏源位置放射性核素的具體分布情況以及對周邊建筑的輻射危害，從而可以為核應急工作者提供參考依據。此外，根據模擬，我們主要得出了以下結論：

(1)湍流漩渦的減弱效應對泄漏源擴散范圍起到了一定的抑制作用，從而影響了核素分布規律。

(2)當源高位于建筑物以下，有樹條件下，去流方向附近的建筑放射性污染較來流方向附近的建筑嚴重。相反，在無樹條件下來流方向附近的建筑放射性污染較去流方向附近的建筑嚴重。

(3)源高位于建筑物以上樹以下，有樹條件下，泄漏源均勻向四周擴散，對道路兩邊建筑造成的放射性污染相當。而在無樹條件下，去流方向附近的建筑物的放射性污染比來流方向附近的建筑物污染要嚴重。

因此，當泄漏源發生在不同高度時，泄漏源對道路兩邊的建筑所造成的放射性污染程度都是明顯不同的。在此基礎上，我們對核應急工作者提出了兩點建議。對于發生乏燃料運輸事故后處理，應首先判斷泄漏源、樹木、建筑物的相對高度，然后通過放射性核素彌散模擬的方法估計放射性危害區域，從而為核事故后的放射性處置提供更加真實的信息。

致謝

本論文得到中國科學院信息化專項項目(XXH13 506-104)，中國科學院青年創新促進會專項項目資助。同時衷心感謝鳳麟團隊全體成員的幫助和指導。