邰 揚(yáng)，申世飛

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

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

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

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

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

1 耦合模型

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

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

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

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

(1)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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