高化超，葛義軍，仲晨華

(海軍工程大學(xué) 艦船工程系，武漢　430033)

?

基于高斯分段煙羽模型的海上核素彌散模型

高化超，葛義軍，仲晨華

(海軍工程大學(xué) 艦船工程系，武漢430033)

摘要：根據(jù)海上核事故核素彌散特點，選擇高斯模型進(jìn)行模擬，基于高斯分段煙羽模型對事故發(fā)生后較短時間內(nèi)的核素彌散情況進(jìn)行了模擬分析，通過Matlab編程模擬煙羽擴(kuò)散形式，得出了短時間內(nèi)事故艦船周圍的核素分布情況，從而對受時間和風(fēng)場影響的核素彌散過程實現(xiàn)了實時煙羽動態(tài)模擬。整個核素彌散過程可預(yù)報一定時間內(nèi)核素彌散路徑和影響范圍，為核救援方案提供參考。

關(guān)鍵詞：核素彌散；艦船；核救援；高斯分段煙羽模型

Citation format:GAO Hua-chao, GE Yi-jun, ZHONG Chen-hua.Model Analysis of Nuclide Dispersion by Segmented Gaussian Plume Model on the Ocean[J].Journal of Ordnance Equipment Engineering,2016(3):116-120.

近年來我國核事業(yè)有了快速的發(fā)展，與之相對應(yīng)的核安全問題越來越受到普遍的重視。雖然各種核動力裝置在設(shè)計、建造、調(diào)試、運行等階段嚴(yán)格遵從核安全法規(guī)，尤其是海上核設(shè)施在設(shè)計上采取了“多道屏障”、“縱深防御”等安全設(shè)計原則和防止事故發(fā)生的多重安全措施，但由于設(shè)備和人員這兩大固有因素，涉核艦船核動力裝置在大海之上仍存在發(fā)生核事故的潛在風(fēng)險。特別要注意的是，涉核艦船在海上運行，由于機(jī)動性的要求，工況變化復(fù)雜，同時受到海浪、風(fēng)暴、洋流等因素的影響，因此涉核艦船海上事故發(fā)生后搶救難度更大且核事故的環(huán)境擴(kuò)散影響更為復(fù)雜。

核素海上大氣彌散模型在核事故應(yīng)急決策支持系統(tǒng)中的作用表現(xiàn)以下方面[1]：

1) 核事故發(fā)生后制定應(yīng)急救援的海上接近路徑；

2) 在核素擴(kuò)散期間實時預(yù)報放射性煙羽的核素濃度及擴(kuò)散范圍分布情況；

3) 對核事故后核素彌散方向的船舶及人員設(shè)備等進(jìn)行疏散或作一定防護(hù)處理。

國內(nèi)的大氣彌散原理研究由3種理論體系支撐，分別為梯度輸送理論、湍流統(tǒng)計理論和相似理論[2]。這3種理論分別基于不同的物理機(jī)制和不同的假定條件[3]。

針對海上核事故的突發(fā)情況，因救援時需要事故周圍較小區(qū)域內(nèi)的核素彌散情況，因此利用高斯模型進(jìn)行較小范圍內(nèi)的煙團(tuán)擴(kuò)散模擬分析較為真實可靠；然后針對模擬得到的結(jié)果結(jié)合實際煙羽彌散特點進(jìn)行比較，從而得到適合高斯模型的煙羽彌散限制邊界條件。因此對于海上核事故的放射性泄漏的意外情況進(jìn)行核素擴(kuò)散仿真驗證的過程對實際中的核救援方案可提供參考，降低了可能遇到的輻射危險。

1計算原理

事故后部分核素易揮發(fā)形成具有一定放射性的云團(tuán)，放射性元素較易隨風(fēng)飄移。需要重點關(guān)注的是其中具有高放射性、長壽命、易沾染的核素。燃料裂變可產(chǎn)生300種以上不同核素，只有少數(shù)是穩(wěn)定的，其余大部分為β衰變的放射性核素。主要的放射性核素如表1所示。

表1　事故產(chǎn)物揮發(fā)性分組

事故發(fā)生后彌散至空氣中的核素種類有大致范圍，但具體某核素的比活度值并不確定，目前大多通過灰塵沉積放射性檢測得到部分原始數(shù)據(jù)[3]。針對彌散核素的特點，依據(jù)風(fēng)場和監(jiān)測點放射性檢測數(shù)據(jù)，高斯模型下不同煙羽軌跡的彌散模擬有助于核事故后現(xiàn)場快速救援工作和事故風(fēng)場后彌散核素軌跡還原工作的展開。高斯模型通過采用較少量的、非球形、類似高斯煙羽的長形煙團(tuán)依據(jù)時間序列模擬煙羽釋放和擴(kuò)散。與其他模型相比，分段煙羽軌跡模型最顯著的區(qū)別在于所有的污染物都是沿著一個既定的軌跡輸送[2]。

1.1高斯分段煙羽模型理論公式

(1)

式(1)中：Cn(x,y,z)為在第n個時間段內(nèi)于(x,y,z)點的放射性核素濃度(Bq/m3)；Qn為在第n個時間段內(nèi)放射性核素源的釋放速率(Bq/s)；σy,σz則分別表示水平方向和垂直方向的擴(kuò)散參數(shù)；un為在第n個時間段內(nèi)于有效高度處的煙羽速度(m/s)；yn為在第n個時間段內(nèi)計算點距離彌散煙羽中心軌跡的垂直距離(m)；Δfn為在第n個時間段里隨時間變化下的積分修正系數(shù)。

1.2高斯分段煙羽模型計算思路

高斯分段煙羽模型是在高斯模型的基礎(chǔ)上，為了對變化風(fēng)場環(huán)境條件下大氣擴(kuò)散過程進(jìn)行快速模擬計算而研究形成的[4-5]。在利用模型進(jìn)行計算的過程中，首先確定相關(guān)現(xiàn)場氣象參數(shù)的數(shù)值，如大氣穩(wěn)定度等，然后建立計算模型，代入數(shù)據(jù)得到計算結(jié)果。高斯分段煙羽模型的煙羽段彌散軌跡如圖1所示。

圖1　高斯分段煙羽模型煙羽段軌跡示意圖

1.3模型參數(shù)

利用高斯分段煙羽模型計算煙羽彌散效果時，公式中相關(guān)參數(shù)的確定方法如下：

1) 大氣穩(wěn)定度

國家標(biāo)準(zhǔn)(GB/T 13201—1991)中規(guī)定了地方大氣污染物排放標(biāo)準(zhǔn)的判別方法，進(jìn)而劃分大氣穩(wěn)定度的高低級別，一共分為6級標(biāo)準(zhǔn)，分別為A～F，其中A級為極不穩(wěn)定；F級則為極穩(wěn)定。

大氣穩(wěn)定度的級別參考表如表2所示。

表2　大氣穩(wěn)定度的級別參考表

2) 大氣擴(kuò)散參數(shù)

在高斯分段煙羽模型中，擴(kuò)散參數(shù)σy、σz是表征擴(kuò)散的范圍及擴(kuò)散的速率大小的量，計算中應(yīng)用較多的是沃格特方法和Briggs方法，在本模型建立中主要采用Briggs方法。Briggs方法對于不同地形有其相適應(yīng)的取值方案，表3中為開闊平原田野地帶的估算方法。

表3　Briggs 擴(kuò)散參數(shù)(開闊平原田野)

2計算流程

核素彌散的模擬程序工作流程如圖2所示。首先將相關(guān)參數(shù)確定完畢，依據(jù)高斯分段煙羽模型作為計算核心，然后對煙羽段中心及擴(kuò)散區(qū)域進(jìn)行確定，在空間點擬合循環(huán)過程中確定其核素濃度分布，循環(huán)結(jié)束后完成本煙羽分段相關(guān)計算，最后再將各煙羽分段根據(jù)同樣的計算流程進(jìn)行分析、匯總、輸出，從而得到所有模擬時下次分段煙羽的軌跡模型。

Matlab模型建立部分計算程序如下：

% 主循環(huán)(時間)

for t_ind=1:T_index

con_radia = zeros(Nx,Ny);

% 第t_ind個時間段內(nèi)沿風(fēng)向輸送矢量坐標(biāo)

x_V=speed_wind(t_ind)*delt_t(t_ind)*cos(dir_wind(t_ind)*pi/180);

y_V=speed_wind(t_ind)*delt_t(t_ind)*sin(dir_wind(t_ind)*pi/180);

D_V=speed_wind(t_ind)*delt_t(t_ind);

% 煙羽分段循環(huán)

for plume_index=t_ind:-1:1

% 煙羽段中心軸線起點、終點坐標(biāo)計算

if plume_index > 1

traj_plume(1,1,plume_index)=traj_plume(1,1,plume_index-1)+x_V;

traj_plume(2,1,plume_index)=traj_plume(2,1,plume_index-1)+y_V;

traj_plume(1,2,plume_index)=traj_plume(1,2,plume_index-1)+x_V;

traj_plume(2,2,plume_index)=traj_plume(2,2,plume_index-1)+y_V;

dwd_plume(1,plume_index)=dwd_plume(1,plume_index-1)+D_V;

dwd_plume(2,plume_index)=dwd_plume(2,plume_index-1)+D_V;

theta_plume(1,plume_index)=theta_plume(plume_index-1);

else

traj_plume(1,1,plume_index)=0;

traj_plume(2,1,plume_index)=0;

traj_plume(1,2,plume_index)=x_source+x_V;

traj_plume(2,2,plume_index)=y_source+y_V;

dwd_plume(1,plume_index)=1;

dwd_plume(2,plume_index)=D_V;

theta_plume(1,plume_index)=dir_wind(t_ind);

end

r=dwd_plume(1,plume_index);

coef_plume(1,plume_index)=0.22*r/sqrt(1+0.0004*r);

r=dwd_plume(2,plume_index);

coef_plume(2,plume_index)=0.22*r/sqrt(1+0.0004*r);

% 計算煙羽段擴(kuò)散范圍(梯形區(qū)域，計算四個頂點)

traj_plume(1,3,plume_index) = traj_plume(1,1,plume_index)...

+coef_plume(1,plume_index)*sin(theta_plume(1,plume_index)/57.3);

traj_plume(2,3,plume_index) = traj_plume(2,1,plume_index)...

-coef_plume(1,plume_index)*cos(theta_plume(1,plume_index)/57.3);

traj_plume(1,4,plume_index) = traj_plume(1,2,plume_index)...

+coef_plume(2,plume_index)*sin(theta_plume(1,plume_index)/57.3);

traj_plume(2,4,plume_index) = traj_plume(2,2,plume_index)...

-coef_plume(2,plume_index)*cos(theta_plume(1,plume_index)/57.3);

traj_plume(1,5,plume_index) = traj_plume(1,2,plume_index)...

-coef_plume(2,plume_index)*sin(theta_plume(1,plume_index)/57.3);

traj_plume(2,5,plume_index) = traj_plume(2,2,plume_index)...

+coef_plume(2,plume_index)*cos(theta_plume(1,plume_index)/57.3);

traj_plume(1,6,plume_index) = traj_plume(1,1,plume_index)...

-coef_plume(1,plume_index)*sin(theta_plume(1,plume_index)/57.3);

traj_plume(2,6,plume_index) = traj_plume(2,1,plume_index)...

+coef_plume(1,plume_index)*cos(theta_plume(1,plume_index)/57.3);

% 模擬區(qū)域網(wǎng)格點循環(huán)

圖2　計算流程

3數(shù)值模擬

3.1數(shù)值模擬方案

根據(jù)海上環(huán)境特點，大氣穩(wěn)定度取為C；船用輻射源的泄漏情形可能原因有多種，不易進(jìn)行估算，此處計算中源強(qiáng)釋放速率取為4×107Bq/s(此處取值僅供參考)；事故船舶核素的泄漏從船甲板附近開始彌散的可能性較大，故釋放源高度為10 m，模擬區(qū)域為二維平面[6-11]。

接下來依據(jù)圖2中的計算流程，利用軟件Matlab進(jìn)行高斯分段煙羽模型的編程，并分別模擬定常風(fēng)場和非定常風(fēng)場兩種條件下的放射性濃度分布：定常風(fēng)場(風(fēng)速2 m/s，風(fēng)向西風(fēng)，較大尺度范圍)；非定常風(fēng)場：(風(fēng)速2 m/s，風(fēng)向每個積分時間依據(jù)步長改變，小尺度范圍)。

因目前沒有相關(guān)的第一手資料，在真實場景中小范圍的煙羽數(shù)據(jù)無法進(jìn)行比較。其他相關(guān)參數(shù)的設(shè)置參考福島核電站較遠(yuǎn)距離的煙羽彌散的監(jiān)測資料，以及結(jié)合日本氣象廳的氣象數(shù)據(jù)：事故期間(2011.3.11-2013.4.11)平均溫度6.2℃,平均相對濕度為54.5%[3]。

3.2數(shù)值模擬結(jié)果

在模擬結(jié)果中，橫坐標(biāo)和縱坐標(biāo)均代表空間位置信息，整體結(jié)果以俯視圖的形式呈現(xiàn)；煙羽中核素濃度以顏色進(jìn)行區(qū)分，核素含量標(biāo)準(zhǔn)于右側(cè)顯示。定常風(fēng)場下的核素彌散其湍流擴(kuò)散系數(shù)K為常數(shù),風(fēng)場不變,核素濃度符合高斯分布，故采用高斯煙羽模型；非定常風(fēng)場下，其彌散核素方向隨風(fēng)向而變，故適宜采用高斯分段煙羽模型。

1) 定常風(fēng)場較大尺度范圍模擬結(jié)果(高斯煙羽模型)

圖3中可以看出，定常風(fēng)場下核素的彌散沿某固定方向，對于沿風(fēng)場方向的核素彌散較為可靠，可對事故發(fā)生后沿風(fēng)場方向的核素濃度進(jìn)行預(yù)報，高斯分布的特征較為明顯。在海上核事故發(fā)生的后期，煙羽末端的彌散受海上多種因素的影響，同時結(jié)合福島核電站事故后放射性的監(jiān)測數(shù)據(jù)，可知在經(jīng)過一定時間的彌散后，其煙羽分布界限將會逐漸淡化。

圖3　較大尺度范圍及定常風(fēng)場下高斯煙羽模擬

2) 非定常風(fēng)場模擬初始彌散結(jié)果(高斯分段煙羽模型)

從圖4中可以看出，采用高斯分段煙羽模型后可以模擬風(fēng)場初始變化時的煙羽濃度，在風(fēng)場變化的情況下可對核素彌漫煙羽的擴(kuò)散軌跡和分布進(jìn)行實時模擬和預(yù)先預(yù)報，較為符合復(fù)雜風(fēng)場下的情形，且風(fēng)速和風(fēng)向均可隨實際情況而變化，需要的現(xiàn)場氣象條件較少。不過因計算過程中采用煙羽分段循環(huán)擬合，存在各煙羽分段之間濃度不連續(xù)的情況，各分段間存在一定計算間隙。

在煙羽初始階段模擬狀態(tài)下，對于煙羽的彌散較為符合煙氣的彌散特點，因此在初始狀態(tài)下的煙羽模擬與真實核素?zé)熡饛浬⒌那樾屋^為接近，能夠?qū)π》秶鷧^(qū)域內(nèi)的核素濃度分布依據(jù)風(fēng)場情況進(jìn)行預(yù)報。

圖4　非定常風(fēng)場下高斯分段煙羽初始模擬

通過模型分析和對比，可以看出兩種高斯模式煙羽的分布情形特點：

定常風(fēng)場要求海上氣象環(huán)境變化較小，適合在核事故發(fā)生后較短的時間內(nèi)且事故地點周圍的風(fēng)場為定常、環(huán)境變化不明顯時的彌散模擬；非定常風(fēng)場下的高斯分段煙羽模型適合在核事故地點為風(fēng)場周期性變化且無極端氣象條件的情況模擬，但在大尺寸情形下對現(xiàn)場的氣象數(shù)據(jù)收集有一定要求。

4結(jié)論

1) 利用高斯分段煙羽模型可模擬并預(yù)報核事故后核素彌散軌跡運動情況，其結(jié)果在無極端天氣條件的情況下是可行的，且能夠?qū)司仍暗木仍桨柑峁┛蛇x接近路線，也可及時對事故地點周圍的船舶發(fā)出預(yù)警，從而避開危險區(qū)域。

2) 高斯煙羽模型輸入少量現(xiàn)場氣象數(shù)據(jù)即可得出短時間內(nèi)的核素彌散情況，整個高斯分段煙羽模型模擬預(yù)報過程方便快捷，適合對彌散過程進(jìn)行快速計算。

3) 高斯煙羽模型在事故地點的較小范圍內(nèi)，可利用非定常風(fēng)場下的分段煙羽模型和此區(qū)域內(nèi)詳細(xì)的氣象數(shù)據(jù)對煙羽的軌跡、濃度等進(jìn)行預(yù)報，在小范圍內(nèi)的預(yù)報是有效的，但大尺度下核素彌散濃度水平的數(shù)值模擬需要更詳細(xì)的氣象數(shù)據(jù)支撐和新的計算模型嵌套。

參考文獻(xiàn)：

[1]張兵.放射性物質(zhì)大氣彌散仿真研究[D].哈爾濱：哈爾濱工程大學(xué)，2004.

[2]余琦，劉原中.分段煙羽模型和煙團(tuán)模型在核事故應(yīng)急中的應(yīng)用比較[J].核科學(xué)與工程，2001，21(3)：288-292.

[3]寧莎莎.福島核事故典型氣載放射性核素的彌散及輻射劑量研究[D].上海：上海交通大學(xué)，2013.

[4]李華，鄧?yán)^勇，王旭輝,等.用高斯模型計算大氣中放射性核素云團(tuán)的擴(kuò)散[J].輻射防護(hù)，2004，24(2):92-99.

[5]陳曉秋.核事故早期應(yīng)急響應(yīng)的風(fēng)場和煙羽濃度預(yù)測模式研究[D].北京：中國原子能科學(xué)研究院，2003.

[6]姚仁太，郝宏偉，胡二邦.RODOS系統(tǒng)中兩種大氣彌散模型鏈的比較[J].輻射防護(hù)，2003，23(3)：146-155.

[7]徐志燕.隨機(jī)游走粒子-煙團(tuán)大氣模式在核事故中的應(yīng)用[D].哈爾濱:哈爾濱工程大學(xué),2011.

[8]湯彬，葛良全，方方.核輻射測量原理[M].哈爾濱:哈爾濱工程大學(xué)出版社,2011.

[9]吳建平，蔡 琦，袁 燦，等.艦船核動力裝置事故處理規(guī)程的改進(jìn)[J].四川兵工學(xué)報,2014,35(11):137-139.

[10]中國科學(xué)院工程力學(xué)研究所.γ射線屏蔽參數(shù)手冊[M].北京:原子能出版社,1976.

[11]RODOS(WG2)-TN(97)-01.Description of the Atmospheric Dispersion Model ATSTEP[Z].1997.

(責(zé)任編輯楊繼森)

Model Analysis of Nuclide Dispersion by Segmented Gaussian Plume Model on the Ocean

GAO Hua-chao, GE Yi-jun, ZHONG Chen-hua

(Department of Ship Engineering, Naval University of Engineering, Wuhan 430033, China)

Abstract:According to the characteristic of nuclide dispersion on the sea, we chose Gaussian model to simulate, and then simulated and analyzed the nuclide dispersion which happened after the nuclear accident in a short time by segmented Gaussian Plume model, and programmed and simulated smoke plume diffusion forms with the software of Matlab. At last, we got the nuclide concentration distribution around the accident ship, thus accomplished the dynamic simulation of nuclide dispersion which affected by time and wind field after the nuclear accident in real time. The whole progress of nuclide dispersion could forecast the nuclide diffusion path and influence scope in a certain time, providing references for nuclear rescue packages.

Key words:nuclide diffusion; ship; nuclear rescue package; segmented Gaussian Plume model

文章編號：1006-0707(2016)03-0116-05

中圖分類號：U661.74；X591

文獻(xiàn)標(biāo)識碼：A

doi:10.11809/scbgxb2016.03.028

作者簡介：高化超(1990—)，男，碩士研究生，主要從事艦船核防護(hù)研究。

收稿日期：2015-08-30；修回日期：2015-09-12

本文引用格式：高化超，葛義軍，仲晨華.基于高斯分段煙羽模型的海上核素彌散模型[J].兵器裝備工程學(xué)報，2016(3):116-120.

【化學(xué)工程與材料科學(xué)】