邰 揚(yáng)，申世飛

(清華大學(xué) 工程物理系,北京 100084)

近年來(lái)，我國(guó)核電建設(shè)持續(xù)增長(zhǎng)，核電在建規(guī)模長(zhǎng)期保持全球第一，按目前已形成的核電發(fā)展條件預(yù)測(cè)，2021—2035年每年將核準(zhǔn)建設(shè)6～8臺(tái)核電機(jī)組，2035年我國(guó)核電裝機(jī)規(guī)模將達(dá)到1.5～1.8億kW，核電發(fā)電量占比將提升至10%以上，核能在國(guó)民經(jīng)濟(jì)中將發(fā)揮越來(lái)越重要的作用[1]。保證核能的安全使用是開(kāi)發(fā)和利用核能的前提，對(duì)核電來(lái)說(shuō)，除了通過(guò)合理選址進(jìn)行前期安全規(guī)劃和建立并維護(hù)行之有效的安全運(yùn)營(yíng)體系之外，對(duì)核事故發(fā)生后放射性物質(zhì)的擴(kuò)散過(guò)程進(jìn)行模擬并根據(jù)結(jié)果制定安全預(yù)案，也是保障核電安全性的重要方面。以2011年福島核事故為例，大量的放射性物質(zhì)持續(xù)性地排放到大氣中，并在事故發(fā)生后的一個(gè)多月內(nèi)隨大氣運(yùn)動(dòng)而大面積擴(kuò)散，將周圍大面積區(qū)域都拖入放射性污染的威脅中。因此，針對(duì)核能發(fā)電設(shè)施在核事故發(fā)生時(shí)放射性物質(zhì)的擴(kuò)散過(guò)程、時(shí)間與空間分布、衰減規(guī)律進(jìn)行模擬研究并用于核事故后的安全疏散等應(yīng)急決策，具有十分重要的現(xiàn)實(shí)意義。

在對(duì)放射性物質(zhì)大氣擴(kuò)散過(guò)程的評(píng)估中，觀測(cè)站的觀測(cè)結(jié)果是重要數(shù)據(jù)來(lái)源，然而僅依靠事故后的觀測(cè)結(jié)果會(huì)存在覆蓋范圍有限、無(wú)法重現(xiàn)、無(wú)法提前進(jìn)行安全預(yù)案等問(wèn)題，而數(shù)值模擬研究則相對(duì)具有很大的優(yōu)勢(shì)。在常用的放射性物質(zhì)大氣擴(kuò)散的數(shù)值模擬中，基于穩(wěn)態(tài)條件下煙羽在豎直和水平維度上呈現(xiàn)出的高斯分布特征的高斯煙羽模型(Gaussian plume model)是早期最廣泛使用的擴(kuò)散模型，但是高斯煙羽模型存在一些嚴(yán)重的缺陷，例如采取穩(wěn)態(tài)近似而忽略了污染物傳播到探測(cè)器的時(shí)間，且忽略了煙羽之間相互作用關(guān)系等假設(shè)造成在小尺度、近源區(qū)域和大尺度上都存在明顯缺陷[2]。因此，常用于模擬放射性物質(zhì)大氣擴(kuò)散過(guò)程的模型大多為拉格朗日模型(Lagrangian model)和歐拉模型(Eulerian model)[3]。

拉格朗日模型將放射性物質(zhì)在大氣中的擴(kuò)散過(guò)程視作粒子團(tuán)在大氣中的運(yùn)動(dòng)，利用擴(kuò)散物理過(guò)程中不同途徑產(chǎn)生的概率函數(shù)來(lái)追蹤粒子團(tuán)的軌跡，再將大量粒子團(tuán)的擴(kuò)散結(jié)果疊加，從而得到放射性物質(zhì)的時(shí)空分布。因?yàn)槔窭嗜漳Ｐ偷年P(guān)注點(diǎn)在粒子團(tuán)本身，所以對(duì)大氣條件的要求較低，在非穩(wěn)態(tài)和非正態(tài)分布條件下依然可以得到較好的準(zhǔn)確度，適合用來(lái)模擬核事故時(shí)近源區(qū)域的放射性物質(zhì)擴(kuò)散。但是由于粒子團(tuán)會(huì)隨著模擬空間的擴(kuò)大而迅速分散，因此當(dāng)模擬尺度增大時(shí)，為了保證遠(yuǎn)源條件下的準(zhǔn)確度和平滑性，模擬所需要的粒子團(tuán)數(shù)目將大幅上升，計(jì)算量也會(huì)大幅增大。相比之下，歐拉模型將模擬區(qū)域柵格化處理后只考慮放射性物質(zhì)在相鄰柵格間的傳輸，由于柵格尺寸的靈活性，歐拉模型可以靈活控制計(jì)算消耗。但歐拉模型也采取了穩(wěn)態(tài)假設(shè)和正態(tài)分布假設(shè)，且在柵格化的過(guò)程中將一般為點(diǎn)排放的排放源近似成一個(gè)柵格，因此在非穩(wěn)態(tài)條件和近源區(qū)域存在著明顯的誤差[4]。

因此，將拉格朗日模型和歐拉模型結(jié)合起來(lái)的耦合模型成為一種可以兼顧拉格朗日模型和歐拉模型優(yōu)點(diǎn)的模擬方法。國(guó)外已有學(xué)者對(duì)此做了一些研究，如OETTL等[5]建立了可以在非穩(wěn)態(tài)大氣條件下使用的拉格朗日-歐拉耦合模型(GRAL模型)，但該模型不包含化學(xué)變化和物理衰變模塊；HURLEY等[6]開(kāi)發(fā)了用于非靜態(tài)大氣擴(kuò)散的耦合模型(TAPM模型)，但該模型未考慮氣溶膠擴(kuò)散的情況，且濕沉降部分僅能處理可溶性氣體；還有一種廣泛使用的ARIA模型[7]，可以支持拉格朗日模式和歐拉模式，但未對(duì)兩種模型進(jìn)行耦合。考慮到大多數(shù)模型未考慮放射性物質(zhì)的沉降與衰變特點(diǎn)，且在中尺度和大尺度上普遍存在準(zhǔn)確性和易用性上的不足[8]，同時(shí)對(duì)于耦合判據(jù)和耦合過(guò)程缺乏深入的探討，筆者將基于兩款成熟的大氣擴(kuò)散模型，即采取拉格朗日模型的FLEXPART[9]和歐拉模型的WRF-CHEM[10]開(kāi)發(fā)一款適用于核事故條件下放射性物質(zhì)在中尺度和大尺度范圍擴(kuò)散的拉格朗日-歐拉耦合模型，并對(duì)耦合方法及其效率進(jìn)行初步討論。

1 耦合模型

筆者所使用的拉格朗日-歐拉模型在接近擴(kuò)散源的位置使用拉格朗日模塊FLEXPART[11]進(jìn)行計(jì)算，而在遠(yuǎn)離擴(kuò)散源的位置，則根據(jù)一定的標(biāo)準(zhǔn)將拉格朗日粒子轉(zhuǎn)化為歐拉濃度利用WRF-CHEM進(jìn)行計(jì)算。本質(zhì)上來(lái)說(shuō)，耦合模式其實(shí)是先用拉格朗日模塊對(duì)排放源進(jìn)行預(yù)處理，再將結(jié)果導(dǎo)入歐拉模塊中作為歐拉模式的初始條件，在這個(gè)過(guò)程中通過(guò)耦合判據(jù)來(lái)判斷放射性物質(zhì)分布是否足夠分散，以減少歐拉模式中柵格平均化所帶來(lái)的誤差。

具體來(lái)說(shuō)，耦合模式需要跟蹤并分別控制拉格朗日模塊中排放源所釋放出來(lái)的每個(gè)粒子團(tuán)[12]，將粒子團(tuán)中的某一個(gè)物理量作為判據(jù)，如擴(kuò)散時(shí)間、標(biāo)準(zhǔn)差等，將其與預(yù)設(shè)判據(jù)基準(zhǔn)做比較，一旦符合轉(zhuǎn)化條件則將該粒子團(tuán)的狀態(tài)量(輸運(yùn)時(shí)間、位置、物質(zhì)質(zhì)量等)記錄下來(lái)，經(jīng)柵格化處理后轉(zhuǎn)入歐拉模式繼續(xù)計(jì)算。耦合模式的耦合判斷過(guò)程如圖1所示。

圖1 拉格朗日-歐拉耦合模式的耦合判斷過(guò)程

在耦合判據(jù)的選取方面，筆者以單個(gè)粒子團(tuán)擴(kuò)散過(guò)程中的彌散程度為判據(jù)，即粒子團(tuán)本身的不確定度S[13]，其表達(dá)式為：

(1)

其中，xi、xj分別表示第i個(gè)和第j個(gè)粒子的位置。本模型在拉格朗日模塊中每隔若干時(shí)間步長(zhǎng)對(duì)粒子團(tuán)的不確定度進(jìn)行一次檢驗(yàn)，若該粒子團(tuán)的不確定度S達(dá)到預(yù)設(shè)閾值S0，則將該粒子團(tuán)的狀態(tài)量[14]轉(zhuǎn)入歐拉模塊中繼續(xù)計(jì)算。

2 模型驗(yàn)證與計(jì)算

利用搭建的拉格朗日-歐拉耦合模型，對(duì)我國(guó)東南沿海某核電站進(jìn)行模擬運(yùn)算(虛擬算例)?？紤]到臺(tái)風(fēng)過(guò)境期間在固定時(shí)間點(diǎn)上較大空間范圍的風(fēng)向、風(fēng)速分布會(huì)比較穩(wěn)定，故選擇臺(tái)風(fēng)過(guò)境期間模擬核電站泄漏事故，這樣不僅可以利用大尺度上的放射性物質(zhì)擴(kuò)散與沉降來(lái)驗(yàn)證模型的有效性，還可以利用中小尺度上的擴(kuò)散與沉降來(lái)比較耦合模型與單一模型的差異性。

2011年8月上旬第9號(hào)超強(qiáng)臺(tái)風(fēng)“梅花”曾掠過(guò)我國(guó)東南沿海并引起嚴(yán)重的氣象災(zāi)害，假設(shè)該核電站于2011年8月4日0點(diǎn)發(fā)生泄漏，參考福島核事故數(shù)據(jù)，設(shè)置泄漏劑量為一次性排放131I為2.3×1013Bq，137Cs為2.3×1012Bq，模擬時(shí)間為北京時(shí)間2011年8月1日0時(shí)至8月11日0時(shí)。因氣象場(chǎng)為GFS(global forecasting system)數(shù)據(jù)插值后由WRF(weather research and forecasting)模型模擬得出，所以需提前計(jì)算3天的氣象數(shù)據(jù)來(lái)建立穩(wěn)定的氣象場(chǎng)。為保證模擬的有效性，拉格朗日模塊和歐拉模塊均選擇相同的微物理參數(shù)化方案[15]和干濕沉降模型[16]，并考慮放射性物質(zhì)隨時(shí)間的衰變。為了進(jìn)一步說(shuō)明模型的有效性，將僅采用歐拉模型進(jìn)行計(jì)算的模擬算例作為參照組。

首先對(duì)模擬的氣象條件進(jìn)行了驗(yàn)證，選擇了泄漏點(diǎn)200 km范圍內(nèi)7個(gè)氣象觀測(cè)站的觀測(cè)數(shù)值進(jìn)行驗(yàn)證。因超強(qiáng)臺(tái)風(fēng)帶來(lái)的強(qiáng)降水幅度巨大，大部分觀測(cè)站未能公布降水的分時(shí)詳細(xì)數(shù)據(jù)，故選取觀測(cè)數(shù)據(jù)較為詳實(shí)的風(fēng)向和風(fēng)速數(shù)據(jù)作為參照。觀測(cè)站點(diǎn)分別位于櫟社、蕭山、浦東、嵊泗、虹橋、定海和石浦(以下分別稱為觀測(cè)站1～觀測(cè)站7)。風(fēng)速和風(fēng)向的模擬值與觀測(cè)值的對(duì)比結(jié)果分別如圖2和圖3所示。

圖2 風(fēng)速的模擬值與觀測(cè)值

圖3 風(fēng)向的模擬值與觀測(cè)值

從圖2可以看出，在大多數(shù)觀測(cè)站，模擬的風(fēng)速值與觀測(cè)數(shù)據(jù)吻合較好，很好地反映了模擬時(shí)段內(nèi)風(fēng)速隨時(shí)間的變化規(guī)律。

在圖3中，y軸為風(fēng)向與正北方向的夾角，由于對(duì)角度來(lái)說(shuō)360°等價(jià)于0°，故圖上最高點(diǎn)和最低點(diǎn)其實(shí)是等價(jià)的。從圖3可以看出，在大多數(shù)觀測(cè)站，模擬的風(fēng)向值與觀測(cè)數(shù)據(jù)吻合較好，可以較好地反映模擬時(shí)段內(nèi)風(fēng)向隨時(shí)間的變化規(guī)律。對(duì)于放射性物質(zhì)的沉降，對(duì)照組和考察組采用了相同的模擬參數(shù)設(shè)置，但模擬流程不同，即對(duì)照組只采用歐拉模型的WRF-CHEM進(jìn)行模擬，而考察組先用拉格朗日模型的FLEXPART較精確地模擬點(diǎn)排放源釋放的過(guò)程，經(jīng)過(guò)耦合模型處理后再轉(zhuǎn)入歐拉模塊處理。

從大尺度上驗(yàn)證模型的有效性，筆者選取1 500 km范圍的大尺度分別對(duì)對(duì)照組和考察組進(jìn)行模擬，得到源排放后12 h、24 h、48 h131I的總沉降量,如圖4所示。從圖4可以看出，考察組在大尺度上的模擬結(jié)果與對(duì)照組基本一致，說(shuō)明耦合模型可以很好地模擬放射性物質(zhì)在大氣中的擴(kuò)散過(guò)程。

圖4 大尺度下對(duì)照組(左)與考察組(右)131I的總沉降量(單位:Bq/h)

從中尺度上對(duì)模型有效性進(jìn)行驗(yàn)證，將模擬尺度切換到500 km范圍內(nèi)的中尺度上分別對(duì)對(duì)照組和考察組進(jìn)行模擬，得到源排放后12 h、24 h、48 h131I的總沉降量，如圖5所示。從圖5可以看出，在中尺度上對(duì)照組和考察組的放射性物質(zhì)沉降結(jié)果表現(xiàn)出較明顯的差異，尤其是在相對(duì)近源區(qū)域(100 km以內(nèi))差異更明顯。由于強(qiáng)臺(tái)風(fēng)帶來(lái)的高風(fēng)速和強(qiáng)降水影響，甚至在距離泄漏源200 km的區(qū)域依然有很大差異，在實(shí)際指導(dǎo)疏散的過(guò)程中將給疏散區(qū)域的判斷帶來(lái)巨大影響。

圖5 中尺度下對(duì)照組(左)與考察組(右)131I的總沉降量(單位:Bq/h)

中尺度下對(duì)照組和考察組137Cs的總沉降量如圖6所示，可見(jiàn)對(duì)137Cs的模擬也驗(yàn)證了這一結(jié)果，即對(duì)照組與考察組差異明顯。這一結(jié)果表明，在中尺度范圍內(nèi)，耦合模型對(duì)于解決歐拉模型計(jì)算放射性物質(zhì)擴(kuò)散時(shí)將排放源均勻化到整個(gè)計(jì)算格點(diǎn)的假設(shè)所帶來(lái)誤差的問(wèn)題上有著非常好的效果。

圖6 中尺度下對(duì)照組(左)與考察組(右)137Cs的總沉降量(單位:Bq/h)

3 結(jié)論

筆者建立了適用于模擬核事故下放射性物質(zhì)大氣擴(kuò)散的拉格朗日-歐拉耦合模型，并對(duì)模型進(jìn)行了算例驗(yàn)證與分析。結(jié)果表明，在中尺度和大尺度下拉格朗日-歐拉耦合模型均有較好的模擬效果，且耦合模型可以在中尺度范圍內(nèi)較大幅度地減小歐拉模型的初始條件所帶來(lái)的系統(tǒng)誤差。