陽林峰，楊宏偉，駱志平，吳建平，程衛(wèi)亞

（1．中國原子能科學(xué)研究院輻射安全研究所，北京 102413；2．成都理工大學(xué)核技術(shù)與自動(dòng)化學(xué)院，四川成都 610059）

利用核函數(shù)提高隨機(jī)游走模式中污染物濃度計(jì)算的效率

陽林峰1，2，楊宏偉1，駱志平1，吳建平2，程衛(wèi)亞1

（1．中國原子能科學(xué)研究院輻射安全研究所，北京 102413；2．成都理工大學(xué)核技術(shù)與自動(dòng)化學(xué)院，四川成都 610059）

本文在拉格朗日隨機(jī)游走模式中引入核函數(shù)法代替質(zhì)點(diǎn)法進(jìn)行濃度計(jì)算，并提出了對(duì)核函數(shù)的頻寬進(jìn)行變頻寬處理，將其與對(duì)流擴(kuò)散時(shí)間相關(guān)聯(lián)，形成變頻寬核函數(shù)。通過對(duì)固定頻寬高斯核函數(shù)、固定頻寬伊番科尼可夫核函數(shù)、變頻寬高斯核函數(shù)和變頻寬伊番科尼可夫核函數(shù)的濃度計(jì)算結(jié)果與質(zhì)點(diǎn)法以及解析解的濃度計(jì)算結(jié)果進(jìn)行比較，發(fā)現(xiàn)變頻寬伊番科尼可夫核函數(shù)在不增加模擬粒子數(shù)的情況下能有效提高濃度計(jì)算結(jié)果的精確性。

核事故后果評(píng)價(jià)；隨機(jī)游走；濃度計(jì)算；核函數(shù)法；變頻寬

隨著福島第一核電站泄漏事故的發(fā)生，污染物在海洋中的對(duì)流擴(kuò)散特性愈發(fā)受到外界關(guān)注［13］。在事故后，通過數(shù)值模擬模式及時(shí)、準(zhǔn)確地給出污染物的濃度分布成為當(dāng)前研究的熱點(diǎn)。目前，國內(nèi)外用于研究污染物擴(kuò)散的數(shù)值模擬模式主要有歐拉模式和拉格朗日模式兩大類［4－8］。因拉格朗日模式在數(shù)值模擬計(jì)算過程中避免了處理歐拉模式的閉合問題而被廣泛采用。在拉格朗日隨機(jī)游走模式濃度計(jì)算部分中常采用質(zhì)點(diǎn)法來計(jì)算污染物的濃度分布［9－10］，即在統(tǒng)計(jì)粒子最終位置分布的過程中將每個(gè)數(shù)值粒子看成帶有質(zhì)量的粒子，因此污染物的最終濃度分布可通過相應(yīng)網(wǎng)格內(nèi)的粒子數(shù)乘以每個(gè)粒子所攜帶的質(zhì)量除以該網(wǎng)格的體積得到，但這種方法在粒子數(shù)較少時(shí)常因粒子對(duì)流擴(kuò)散距離增加，導(dǎo)致計(jì)算網(wǎng)格內(nèi)粒子數(shù)稀疏，使得最終的濃度計(jì)算結(jié)果出現(xiàn)大幅漲落，增加了數(shù)值計(jì)算結(jié)果的不確定度。為提高該方法計(jì)算濃度結(jié)果的精確性，可通過采取增加粒子數(shù)的方式，但增加粒子數(shù)的同時(shí)也相應(yīng)地增加了數(shù)值模擬計(jì)算時(shí)間。因此，為在不增加數(shù)值模擬計(jì)算時(shí)間的基礎(chǔ)上，提高濃度計(jì)算結(jié)果的精確性，本工作考慮引入核函數(shù)法用于拉格朗日隨機(jī)游走模式中的濃度計(jì)算。

1 核函數(shù)計(jì)算方法

核函數(shù)法計(jì)算污染物濃度，即將對(duì)流擴(kuò)散后每個(gè)粒子對(duì)最終污染物濃度分布的貢獻(xiàn)看成是某一函數(shù)的分布［11］。因此，采用核函數(shù)后拉格朗日隨機(jī)游走模式中粒子對(duì)空間各位置總的濃度分布貢獻(xiàn)可表示為：

高斯核函數(shù)：

伊番科尼可夫核函數(shù)：

式中，νd為維度d的單位球體積，νd＝2πd／2／dΓ（d／2），Γ（x）為伽馬函數(shù)。

在核函數(shù)中頻寬主要是用來度量單個(gè)粒子在進(jìn)行濃度分布貢獻(xiàn)計(jì)算過程中的有效距離［12］。污染物輸運(yùn)過程中，兩相鄰粒子間距隨對(duì)流擴(kuò)散時(shí)間的增加不斷增大，若繼續(xù)采用固定頻寬的核函數(shù)，則當(dāng)兩相鄰粒子間距超過兩倍核函數(shù)的頻寬時(shí)，其每個(gè)粒子類似于質(zhì)點(diǎn)法被看成單個(gè)帶有質(zhì)量的粒子。為避免這種情況，在固定頻寬的基礎(chǔ)上，提出對(duì)核函數(shù)中的頻寬進(jìn)行變頻寬處理，即將核函數(shù)中的頻寬與對(duì)流擴(kuò)散時(shí)間相關(guān)聯(lián)，其處理方式如下式：

其中，ci（i為x、y、z）為常數(shù)。

由此，通過上述處理，在固定頻寬高斯核函數(shù)和固定頻寬伊番科尼可夫核函數(shù)的基礎(chǔ)上，分別得到變頻寬高斯核函數(shù)和變頻寬伊番科尼可夫核函數(shù)。在拉格朗日隨機(jī)游走濃度計(jì)算中，通過這4種核函數(shù)可將單個(gè)粒子對(duì)空間各位置處的濃度分布進(jìn)行函數(shù)化。

2 基于核函數(shù)的隨機(jī)游走濃度計(jì)算

2．1 濃度計(jì)算模型構(gòu)建

為此，當(dāng)t0＝dt時(shí)，隨機(jī)游走模型無量綱化處理后為：

式中，dW′x、dW′y、dW′z分別為平均值為0、方差為1的高斯隨機(jī)數(shù)。

對(duì)僅考慮二階全反射的相應(yīng)解析解無量綱處理后為（下式中各變量省略了上標(biāo)）：

式中：h為污染源排放點(diǎn)垂直高度；（xs，ys，zs）為污染源的釋放點(diǎn)位置。

2．2 不同核函數(shù)法的濃度計(jì)算結(jié)果比較

在濃度結(jié)果比較過程中，模型的輸入?yún)?shù)均相同，其中水平流速u＝1m／s，垂直高度h＝30m，污染源釋放點(diǎn)位置為（0，0，－10）。4種核函數(shù)法的濃度計(jì)算結(jié)果與質(zhì)點(diǎn)法以及解析解所得到的濃度計(jì)算結(jié)果在中心線上時(shí)間t＝10、100、500、800個(gè)單位時(shí)的比較情況如圖1～4所示（為便于觀察，對(duì)縱坐標(biāo)污染物濃度均放大1 000倍）。

圖1 t為10個(gè)單位時(shí)中心線上的濃度分布Fig．1 Concentration distribution in center line at t＝10units

從圖1可見，當(dāng)t為10個(gè)單位時(shí)，即在污染物釋放源附近，由于粒子剛開始對(duì)流擴(kuò)散，粒子相對(duì)集中。因此采用質(zhì)點(diǎn)法、固定頻寬高斯核函數(shù)法、固定頻寬伊番科尼可夫核函數(shù)法、變頻寬高斯核函數(shù)法以及變頻寬伊番科尼可夫核函數(shù)法，其濃度計(jì)算結(jié)果曲線均未出現(xiàn)震蕩；變頻寬高斯核函數(shù)法因高斯函數(shù)本身考慮全空間區(qū)域在計(jì)算過程中有截?cái)?，因此其濃度結(jié)果計(jì)算值較解析解的計(jì)算值低。

圖2 t為100個(gè)單位時(shí)中心線上的濃度分布Fig．2 Concentration distribution in center line at t＝100units

圖3 t為500個(gè)單位時(shí)中心線上的濃度分布Fig．3 Concentration distribution in center line at t＝500units

由圖2b可知，隨著對(duì)流擴(kuò)散的進(jìn)行，當(dāng)對(duì)流擴(kuò)散時(shí)間t增加到100個(gè)單位時(shí)，分別采用質(zhì)點(diǎn)法和固定頻寬伊番科尼可夫核函數(shù)法計(jì)算濃度時(shí)，其濃度結(jié)果曲線已開始出現(xiàn)較大的震蕩，前者主要是因?yàn)殡S粒子對(duì)流擴(kuò)散距離的增加在部分網(wǎng)格內(nèi)粒子數(shù)稀疏，因此采用質(zhì)點(diǎn)法計(jì)算濃度時(shí)會(huì)出現(xiàn)相鄰網(wǎng)格之間的濃度值大幅漲落；后者盡管在計(jì)算過程中將每個(gè)粒子對(duì)濃度的貢獻(xiàn)函數(shù)化，但隨著相鄰兩粒子間距超過兩倍固定頻寬伊番科尼可夫核函數(shù)的頻寬時(shí)，在計(jì)算過程中每個(gè)粒子被當(dāng)作單個(gè)帶有質(zhì)量的粒子來計(jì)算，因此此時(shí)采用固定頻寬伊番科尼可夫核函數(shù)，其濃度結(jié)果曲線會(huì)出現(xiàn)震蕩。此外，從圖2a、d可見，此時(shí)采用固定頻寬高斯核函數(shù)和變頻寬伊番科尼可夫核函數(shù)法計(jì)算濃度，其濃度結(jié)果曲線與解析解的濃度結(jié)果曲線符合較好。

由圖3b可知，隨著對(duì)流擴(kuò)散距離的繼續(xù)增加，當(dāng)對(duì)流擴(kuò)散時(shí)間增加到500個(gè)單位時(shí)，采用質(zhì)點(diǎn)法和固定頻寬伊番科尼可夫核函數(shù)法計(jì)算的濃度值曲線震蕩加劇。此時(shí)，如圖3a、d所示，若采用固定頻寬高斯核函數(shù)和變頻寬伊番科尼可夫核函數(shù)法計(jì)算濃度，其濃度結(jié)果曲線與解析解的濃度結(jié)果曲線依然符合較好。

當(dāng)對(duì)流擴(kuò)散時(shí)間t增加到800個(gè)單位時(shí)，由圖4a可見，采用固定頻寬高斯核函數(shù)法計(jì)算的濃度值曲線出現(xiàn)震蕩，這是因?yàn)殡S著相鄰兩粒子間距超過兩倍固定頻寬高斯核函數(shù)的頻寬時(shí)，粒子在計(jì)算過程中被當(dāng)作單個(gè)攜帶質(zhì)量的粒子，因此其相鄰濃度值之間會(huì)出現(xiàn)大幅漲落。相比之下，從圖4c可知，當(dāng)時(shí)間增加到800個(gè)單位時(shí)采用變頻寬高斯核函數(shù)法計(jì)算濃度能有效地克服對(duì)流擴(kuò)散距離增加帶來的震蕩，但高斯函數(shù)在計(jì)算過程中因質(zhì)量虧損其濃度計(jì)算值會(huì)較解析解的值偏低；與此同時(shí)如圖4d所示，若采用變頻寬伊番科尼可夫核函數(shù)法計(jì)算的濃度值曲線與解析解的濃度值曲線依然符合較好。

綜上所述，在污染物釋放源附近或短距離處，因粒子對(duì)流擴(kuò)散距離短、相對(duì)集中，因此無論采用質(zhì)點(diǎn)法還是核函數(shù)法其濃度計(jì)算結(jié)果均較理想，但隨著對(duì)流擴(kuò)散距離增加，質(zhì)點(diǎn)法所計(jì)算的濃度值曲線開始出現(xiàn)震蕩并逐步加劇。若采用固定頻寬核函數(shù)，其抗震能力主要由核函數(shù)的頻寬決定；若采用變頻寬核函數(shù)則能有效地克服震蕩，但在使用過程中因高斯函數(shù)的質(zhì)量虧損，所以變頻寬伊番科尼可夫核函數(shù)是較理想地選擇。

圖4 t為800個(gè)單位時(shí)中心線上的濃度分布Fig．4 Concentration distribution in center line at t＝800units

3 結(jié)論

在隨機(jī)游走模式濃度計(jì)算中用核函數(shù)法代替?zhèn)鹘y(tǒng)的質(zhì)點(diǎn)法計(jì)算濃度時(shí)，通過對(duì)相應(yīng)的核函數(shù)本身進(jìn)行變頻寬處理，能有效克服濃度計(jì)算過程中結(jié)果隨粒子對(duì)流擴(kuò)散距離增加帶來的漲落，這樣在不增加模擬粒子數(shù)的前提下，能有效提高拉格朗日隨機(jī)游走模式中濃度計(jì)算結(jié)果的精確性。

［1］ CHINO M，NAKAYAMA H，NAGAI H，etal．Preliminary estimation of release amounts of131I and137Cs accidentally discharged from the Fukushima Daiichi Nuclear Power Plant into the atmosphere［J］．J Nucl Sci Technol，2011，48：1 129－1 134．

［2］ BUTLER D．Radioactivity spreads in Japan［J］．Nature，2001，471：555－556．

［3］ 喬方利，王關(guān)鎖，趙偉，等．2011年3月日本福島核泄漏物質(zhì)運(yùn)輸擴(kuò)散路徑的情景模擬和預(yù)測［J］．科學(xué)通報(bào)，2011，56（12）：964－971．

QIAO Fangli，WANG Guansuo，ZHAO Wei，et al．Predicting the spread of nuclear radiation from the damaged Fukushima Nuclear Power Plant［J］．Chinese Science Bulletin，2011，56（12）：964－971（in Chinese）．

［4］ KOSHIZUKA S，TAMAKO H，OKA Y．A particle method for incompressible viscous flow with fluid fragmentation［J］．Journal of Computational Fluid Dynamics，1995，4（1）：29－46．

［5］ VALE L M，DIAS J M．Coupling of a Lagrangian particle tracking module to a numerical hydrodynamic model：Simulation of pollution events inside an estuarine port area［J］．Journal of Coastal Research，2011，64：1 609－1 613．

［6］ THOMSON D J．Criteria for the selection of stochastic models of particle trajectories in turbulent flows［J］．J Fluid Mech，1987，180：529－556．

［7］ GRIFFA A，KIRWAN A D，Jr，MARIANO A J．Lagrangian analysis and prediction of coastal and ocean dynamics［M］．Cambridge：Cambridge University Press，2007．

［8］ WANG Cui，SUN Yinglan，ZHANG Xueqing．Study of the three－dimensional Lagrangian model and its application in Jiaozhou Bay［C］∥International Workshop on Education Technology and Training ＆International Workshop on Geoscience and Remote Sensing．［S．l．］：［s．n．］，2008．

［9］ RIDDLE A M．The specification of mixing in random walk models for dispersion in the sea［J］．Cont Shelf Res，1998，18：441－456．

［10］LIU Aihua，KUAI Linping．A review on radionuclides atmospheric dispersion modes［J］．Journal of Meteorology and Environment，2011，27（4）：59－65．

［11］SPIVAKOVSKAYA D，HEEMINK A W，DELEERSNILDER E．Lagrangianmodelling of multidimensional advection－diffusion with space－varying diffusivities：Theory and idealized test cases［J］．Ocean Dynamics，2007，57：189－203．

［12］AIZERMAN M，BRAVERMAN R L．Theoretical foundations of the potentialfunction method in pattern recognitionlearning［J］．Automation and Remote Control，1964，25：821－837．

［13］WADE E H．Simplified approach to particle tracking methods for contaminant transport［J］．Hydraulic Engineering，1997，123（12）：1 157－1 160．

［14］CRISTIANINI N，SHAWE－TAYLOR J．Kernel methods for pattern recognition［M］．Cambridge：Cambridge University Press，2004．

Using Kernel Function to Improve Efficiency of Contaminated Concentration Calculation in Radom Walk Model

YANG Lin－feng1，2，YANG Hong－wei1，LUO Zhi－ping1，WU Jian－ping2，CHENG Wei－ya1
（1．China Institute of Atomic Energy，P．O．Box275－15，Beijing102413，China；
2．College of Nuclear Technology and Automation Engineering，
Chengdu University of Technology，Chengdu610059，China）

In this paper，the kernel function method was introduced into the concentration calculation of the random walk model to replace the box method．And a change to the bandwidth of the kernel function was proposed to form a frequency bandwidth kernel function by relating it to convection diffusion time．By comparing the calculated concentration of analytical solution with the fixed bandwidth Gaussian kernel function，the fixed bandwidth Epanechnikov kernel function，the frequency bandwidth Gaussian kernel function，the frequency bandwidth Epanevhnikov kernel function and the box method，the conclusion was obtained that the frequency bandwidth Epanevhnikov kernel function can efficiently improve the accuracy of the contaminated concentration calculation without increasing simulated particles．

consequence assessment of nuclear accident；random walk；concentration calculation；kernel function method；frequency bandwidth

TL32

：A

1000－6931（2015）03－0485－06

10．7538／yzk．2015．49．03．0485

2013－11－28；

2014－01－21

陽林峰（1986—），男，四川自貢人，碩士研究生，輻射防護(hù)與環(huán)境保護(hù)專業(yè)