談文姬,王德忠,馬元巍,吉志龍
(上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院,上海 200240)
自1986年發(fā)生切爾諾貝利事故以來,國(guó)際上對(duì)核事故后果評(píng)價(jià)和決策支持系統(tǒng)的研究有了很大進(jìn)展,特別是在2011 年福島核事故之后,放射性核素釋放到環(huán)境中的影響引起了國(guó)際社會(huì)的廣泛關(guān)注。核事故后果決策支持系統(tǒng)能對(duì)事故發(fā)生后放射性核素泄漏到空氣中的濃度分布進(jìn)行實(shí)時(shí)和預(yù)測(cè)性計(jì)算,為保護(hù)周圍公眾的安全提供數(shù)據(jù)依據(jù)。
目前,用于核事故應(yīng)急決策支持系統(tǒng)中的大氣擴(kuò)散模型主要有3類:高斯煙羽模型、拉格朗日煙團(tuán)模型和蒙特卡羅粒子模型[1]。高斯煙羽模型適用于穩(wěn)定均勻流動(dòng)及平坦地形的條件,雖模型簡(jiǎn)單、運(yùn)算快捷,但難以適用于復(fù)雜流場(chǎng)下的實(shí)際核事故。蒙特卡羅粒子模型能較為真實(shí)地模擬粒子在大氣中的輸運(yùn)過程,但計(jì)算時(shí)間較長(zhǎng),并存在邊界稀疏的問題,往往不能滿足實(shí)時(shí)核事故應(yīng)急響應(yīng)[2]。煙團(tuán)模型能適用于非穩(wěn)定流動(dòng)及復(fù)雜地形的條件,在局地尺度上有較好的精度,應(yīng)用較為廣泛,如歐洲RODOS等核事故應(yīng)急決策支持系統(tǒng)均采用該模型[3]。
描述煙團(tuán)擴(kuò)散能力的大氣擴(kuò)散系數(shù),是影響拉格朗日煙團(tuán)模型的重要參數(shù)。傳統(tǒng)的拉格朗日煙團(tuán)模型采用的擴(kuò)散參數(shù)是由Pasquill-Gifford(P-G)曲線確定的。P-G 曲線是根據(jù)理想條件下的實(shí)驗(yàn)結(jié)果擬合得到的,與實(shí)際情況的不符將給煙團(tuán)模型帶來偏差。本文針對(duì)拉格朗日煙團(tuán)模型不能實(shí)時(shí)反映局地?cái)U(kuò)散能力的缺點(diǎn),提出大氣擴(kuò)散系數(shù)自適應(yīng)的改進(jìn)方法,并使用實(shí)際大氣擴(kuò)散示蹤實(shí)驗(yàn)的數(shù)據(jù),對(duì)大氣擴(kuò)散系數(shù)自適應(yīng)修正方法進(jìn)行檢驗(yàn)。
煙團(tuán)模型是將釋放的污染物假想成離散的煙團(tuán),對(duì)每個(gè)煙團(tuán)中心的輸運(yùn)過程進(jìn)行模擬。每個(gè)煙團(tuán)在某一特定時(shí)間間隔內(nèi)固定不動(dòng),根據(jù)此刻固定住的煙團(tuán)計(jì)算濃度。然后煙團(tuán)在下一個(gè)時(shí)間步長(zhǎng)內(nèi)繼續(xù)移動(dòng),大小和強(qiáng)度繼續(xù)變化,直到下次采樣時(shí)間再次固定。在基本時(shí)間步長(zhǎng)內(nèi),接受點(diǎn)的濃度為周圍所有煙團(tuán)采樣時(shí)間內(nèi)的平均濃度總和。有風(fēng)的條件下,第i個(gè)煙團(tuán)在某一點(diǎn)的濃度[4]為:
式中:Ci(x,y,z)為第i個(gè)煙團(tuán)(x,y,z)點(diǎn)的濃度,g/m3;Q(i)為第i個(gè)煙團(tuán)的強(qiáng)度,g/s;σx、σy、σz分別為沿風(fēng)向、水平方向、垂直方向的擴(kuò)散參數(shù),m;(xc(i),yc(i),zc(i))為第i個(gè)煙團(tuán)的中心位置,m。
計(jì)算污染物濃度及分布時(shí),必須知道源強(qiáng)、風(fēng)場(chǎng)、有效源高和大氣擴(kuò)散參數(shù)。其中,源強(qiáng)往往是通過測(cè)量獲得或由設(shè)計(jì)基準(zhǔn)事故給出,氣象信息由當(dāng)?shù)貧庀鬁y(cè)量站獲取。大氣擴(kuò)散參數(shù)對(duì)擴(kuò)散后果有較大的影響,但目前仍為理想條件下的示蹤實(shí)驗(yàn)數(shù)據(jù)計(jì)算得到,這些數(shù)據(jù)在核事故后果評(píng)價(jià)中與實(shí)際情況有較大偏差,應(yīng)提出新的方法對(duì)其進(jìn)行修正。
大氣擴(kuò)散參數(shù)(σx、σy、σz)表示污染物在隨機(jī)湍流場(chǎng)中的擴(kuò)散能力和散布范圍,它與下風(fēng)向距離x、大氣穩(wěn)定度和下墊面粗糙度等相關(guān)。大氣擴(kuò)散參數(shù)隨x 及下墊面粗糙度的增大而增大,隨大氣穩(wěn)定度的增大而減小。
大多數(shù)模型均通過Pasquill氣象穩(wěn)定度與P-G 曲線來確定大氣擴(kuò)散參數(shù),該曲線是由大量擴(kuò)散實(shí)驗(yàn)和理論分析得到的,具體確定擴(kuò)散參數(shù)的方法為:首先根據(jù)云況和日射及地面風(fēng)速選擇大氣穩(wěn)定度,然后根據(jù)擴(kuò)散曲線讀出不同下風(fēng)向距離處的擴(kuò)散參數(shù)。RIMPUFF模型中使用的Karlsruhe-Jülich(K-J)大氣擴(kuò)散系數(shù)即采用此方法,K-J大氣擴(kuò)散系數(shù)是根據(jù)大氣穩(wěn)定度及源項(xiàng)的釋放高度將大氣擴(kuò)散系數(shù)進(jìn)行分類,共分為18套大氣擴(kuò)散系數(shù)[5]。但由于P-G曲線是依據(jù)平坦理想條件下的實(shí)驗(yàn)得到的,對(duì)于不同源高、下墊面等條件下的擴(kuò)散曲線還有待進(jìn)一步修正。
本文提出一種根據(jù)當(dāng)?shù)氐膶?shí)測(cè)數(shù)據(jù)(源項(xiàng)釋放數(shù)據(jù)、觀測(cè)點(diǎn)位置和測(cè)得的濃度數(shù)據(jù)、氣象數(shù)據(jù))及煙團(tuán)模型的濃度計(jì)算公式,對(duì)大氣擴(kuò)散系數(shù)進(jìn)行實(shí)時(shí)估計(jì)。這種將上一步大氣擴(kuò)散系數(shù)的估計(jì)結(jié)果作為修正的大氣擴(kuò)散系數(shù)進(jìn)行下一個(gè)時(shí)間步長(zhǎng)運(yùn)算的方法,稱為大氣擴(kuò)散系數(shù)的自適應(yīng)修正。將得到的適用于當(dāng)?shù)氐膶?shí)時(shí)大氣擴(kuò)散系數(shù)作為修正后的系數(shù)繼續(xù)進(jìn)行擴(kuò)散計(jì)算,較使用理想條件下的實(shí)驗(yàn)得到的擴(kuò)散參數(shù)更符合實(shí)際情況,能提高拉格朗日煙團(tuán)模型的準(zhǔn)確性。
假設(shè)擴(kuò)散參數(shù)為下風(fēng)向距離x 的冪函數(shù),即:
式中:py、qy、pz、qz為大氣擴(kuò)散系數(shù)。
若確定py、qy、pz、qz,就 可 得 到 適 合 當(dāng) 時(shí)當(dāng)?shù)氐匦螝庀髼l件的σy和σz。這4個(gè)大氣擴(kuò)散系數(shù)可利用最小二乘法確定。為了權(quán)衡觀測(cè)數(shù)據(jù)的精度,引入一標(biāo)志測(cè)量精度的權(quán)重g,使地面濃度的計(jì)算值Cmi(Cmi=C(xi,yi,0))與實(shí)際觀測(cè)值COi之差的平方和S 最小,則S 可表示為:
式中:n為實(shí)際觀測(cè)的數(shù)據(jù)個(gè)數(shù);gi為采樣點(diǎn)的權(quán)重,gi=COi/CO,max,CO,max為觀測(cè)點(diǎn)采集數(shù) 據(jù)中的最大濃度。
采用梯度迭代法計(jì)算求解使S 為最小的非線性方程組,可得到實(shí)時(shí)的大氣擴(kuò)散系數(shù)。本文根據(jù)上述方法,針對(duì)拉格朗日煙團(tuán)模型進(jìn)行修正,開發(fā)了相應(yīng)的大氣擴(kuò)散系數(shù)自適應(yīng)模型。
MVK 是荷蘭國(guó)家公共健康與環(huán)境研究院開發(fā)的大氣擴(kuò)散模型驗(yàn)證對(duì)比工具[6],包括Kincaid、Indianapolis、Copenhagen、Lillestrom 4個(gè)實(shí)驗(yàn)數(shù)據(jù)集,定量的統(tǒng)計(jì)分析工具(BOOT程序)和定性的圖形分析工具(SIGPLOT 程序),以及可用于評(píng)價(jià)模型的準(zhǔn)確性、適用性和揭示模型的錯(cuò)誤和數(shù)據(jù)的壞點(diǎn)。目前,MVK 已廣泛用于驗(yàn)證和比較自主開發(fā)的大氣擴(kuò)散模型[7]。
Kincaid實(shí)驗(yàn)是EPRI煙羽模型驗(yàn)證和改進(jìn)項(xiàng)目的一部分。EPRI空氣質(zhì)量數(shù)據(jù)中心提供了實(shí)驗(yàn)時(shí)的氣象數(shù)據(jù),包括氣壓、溫度、濕度、風(fēng)速和風(fēng)向。Kincaid電站坐落在美國(guó)伊利諾斯州(北緯39.59°,西經(jīng)89.49°),平均海拔約為180m,其周圍是平坦的農(nóng)田及湖泊。圖1為以Kincaid電站為中心的80km×80km 范圍的等高地形圖。圖1中為UTM 坐標(biāo),x正方向?yàn)闁|,y正方向?yàn)楸保W(wǎng)格劃分為100m×100m。
圖1 以Kincaid電站為中心的80km×80km 范圍的等高地形圖Fig.1 80km×80km wide topographic contour centered at Kincaid Power Plant
實(shí)驗(yàn)時(shí),六氟化硫?qū)碾娬局械囊蛔?87m、直徑9 m 的煙囪頂部釋放。與釋放點(diǎn)距離為1、2、3、5、7和10km 的圓弧上設(shè)有觀測(cè)點(diǎn),采樣并記錄數(shù)據(jù)。實(shí)驗(yàn)數(shù)據(jù)集包含約171h的數(shù)據(jù)。
BOOT程序是MVK中的統(tǒng)計(jì)工具,能將模型預(yù)測(cè)值與實(shí)驗(yàn)觀測(cè)值進(jìn)行比較,且可計(jì)算出各種統(tǒng)計(jì)參數(shù)。除了模型預(yù)測(cè)與實(shí)驗(yàn)觀測(cè)的平均值MEAN、標(biāo)準(zhǔn)差SIGMA 外,其他數(shù)據(jù)還包括平均偏差BIAS、標(biāo)準(zhǔn)化均方差NMSE、相關(guān)系數(shù)CORR、比例偏差FB、FA2等,這些統(tǒng)計(jì)參數(shù)能為模型評(píng)價(jià)提供依據(jù)[8],其計(jì)算公式為:
式中:CM為模型預(yù)測(cè)濃度,μg/m3;CO為觀測(cè)濃度,μg/m3;σM為模型預(yù)測(cè)濃度的標(biāo)準(zhǔn)差,μg/m3;σO為觀測(cè)濃度的標(biāo)準(zhǔn)差,μg/m3。FA2為0.5≤CM/CO≤2的幾率,即模型計(jì)算值誤差在50%以內(nèi)的幾率。
擬合大氣擴(kuò)散系數(shù)所需數(shù)據(jù)包括:觀測(cè)站數(shù)據(jù)(觀測(cè)站相對(duì)于源的位置坐標(biāo)和測(cè)得的濃度數(shù)據(jù))、氣象數(shù)據(jù)(當(dāng)時(shí)的風(fēng)速、風(fēng)向,由氣象站提供)、源項(xiàng)釋放率(實(shí)驗(yàn)設(shè)計(jì)提供)。1980年4月20日14:00的測(cè)量數(shù)據(jù)列于表1。其中源項(xiàng)釋放處的高度為源高處,源高處風(fēng)向以正北方向?yàn)?°。
表1 Kincaid實(shí)驗(yàn)1980年4月20日測(cè)量數(shù)據(jù)Table 1 Measurement data of 1980-04-20in Kincaid experiment
大氣擴(kuò)散系數(shù)初值的設(shè)置會(huì)對(duì)計(jì)算得到的大氣擴(kuò)散系數(shù)帶來較大影響。不同的初值求解得到的大氣擴(kuò)散系數(shù)不同,需對(duì)解得的大氣擴(kuò)散系數(shù)進(jìn)行篩選。本文提出一種大氣擴(kuò)散系數(shù)最優(yōu)擬合方法,根據(jù)表1中的數(shù)據(jù),當(dāng)大氣擴(kuò)散系數(shù)py、qy、pz、qz初值為(1,1,1,1)時(shí),求解得到大氣擴(kuò)散系數(shù)A=(0.739,0.881,0.838,0.959);當(dāng)初值為(0.8,0.8,0.8,0.8)時(shí),解得大氣擴(kuò)散系數(shù)B=(3.84,0.681,0.553,1.007)。將A 和B代入拉格朗日煙團(tuán)模型時(shí),兩者計(jì)算得到的觀測(cè)點(diǎn)坐標(biāo)處的濃度相同,無法立即判斷這兩組大氣擴(kuò)散系數(shù)的適用性。在這種情況下,將A、B 兩組大氣擴(kuò)散系數(shù)作為初值代入下一個(gè)時(shí)間步長(zhǎng)的迭代計(jì)算。下一步計(jì)算結(jié)果表明,使用A 計(jì)算的濃度較使用B 計(jì)算的濃度更接近觀測(cè)值,則選擇A 代入下一步的計(jì)算。重復(fù)該過程即可得到大氣擴(kuò)散系數(shù)的最優(yōu)解。
用MVK 中的Kincaid數(shù)據(jù)集對(duì)大氣擴(kuò)散系數(shù)自適應(yīng)修正方法進(jìn)行檢驗(yàn),并對(duì)大氣擴(kuò)散系數(shù)自適應(yīng)修正方法的煙團(tuán)模型與K-J大氣擴(kuò)散系數(shù)煙團(tuán)模型的計(jì)算結(jié)果進(jìn)行比較,以判斷該方法的可行性。根據(jù)Kincaid實(shí)驗(yàn)記錄,從K-J大氣擴(kuò)散系數(shù)中選擇與源項(xiàng)釋放高度為180m 和D 類大氣穩(wěn)定度相對(duì)應(yīng)的一組大氣擴(kuò)散系數(shù)(0.208,0.903,0.307,0.734)[9]。
對(duì)實(shí)驗(yàn)值(OBS)、K-J大氣擴(kuò)散系數(shù)煙團(tuán)模型的計(jì)算值和大氣擴(kuò)散系數(shù)自適應(yīng)修正的煙團(tuán)模型(APM)的計(jì)算值進(jìn)行對(duì)比的結(jié)果列于表2~6。
表2 1980年4月20日3h18個(gè)有效觀測(cè)站數(shù)據(jù)檢驗(yàn)結(jié)果Table 2 Result of comparison of 18valid experiment data in 3hon 1980-04-20
表3 1980年4月25日6h45個(gè)有效觀測(cè)站數(shù)據(jù)檢驗(yàn)結(jié)果Table 3 Result of comparison of 45valid experiment data in 6hon 1980-04-25
表4 1980年5月1日6h20個(gè)有效觀測(cè)站數(shù)據(jù)檢驗(yàn)結(jié)果Table 4 Result of comparison of 20valid experiment data in 6hon 1980-05-01
表5 1980年5月4日6h19個(gè)有效觀測(cè)站數(shù)據(jù)檢驗(yàn)結(jié)果Table 5 Result of comparison of 19valid experiment data in 6hon 1980-05-04
由表2、3可看出,與大氣擴(kuò)散系數(shù)煙團(tuán)模型相比,大氣擴(kuò)散系數(shù)自適應(yīng)修正方法計(jì)算得到的平均值、標(biāo)準(zhǔn)差更接近于觀測(cè)值,比例偏差FB更接近于0、FA2更接近于1。采取大氣擴(kuò)散系數(shù)自適應(yīng)方法的煙團(tuán)模型的各項(xiàng)統(tǒng)計(jì)參數(shù)均優(yōu)于KJ大氣擴(kuò)散系數(shù)煙團(tuán)模型的。就模型計(jì)算值與觀測(cè)值誤差在50%以內(nèi)的幾率而言,準(zhǔn)確性提高約20%。表2、3、6的數(shù)據(jù)修正效果優(yōu)于表4、5的,這是由于風(fēng)向變化較大對(duì)修正效果產(chǎn)生了負(fù)影響。1980年4月20日、4月25日及5月7日3天的風(fēng)向變化較?。?°~40°),而5月1日、4日兩天的風(fēng)向變化較大,有時(shí)達(dá)180°。
表6 1980年5月7日9h36個(gè)有效觀測(cè)站數(shù)據(jù)檢驗(yàn)結(jié)果Table 6 Result of comparison of 36valid experiment data in 9hon 1980-05-07
圖2、3分別為大氣擴(kuò)散系數(shù)自適應(yīng)修正方法模型及K-J大氣擴(kuò)散系數(shù)模型的預(yù)測(cè)值與觀測(cè)值的分位數(shù)圖(Q-Q 圖)。
圖2 APM 模型的Q-Q 圖Fig.2 Q-Q plot of APM model
圖3 K-J模型的Q-Q 圖Fig.3 Q-Q plot of K-J model
從圖2和3可看出,兩個(gè)模型的預(yù)測(cè)值較實(shí)驗(yàn)值均偏大,這種保守的估計(jì),在核事故濃度預(yù)測(cè)上是可被接受的。APM 模型點(diǎn)的個(gè)數(shù)多于K-J模型的,這是由于APM 模型預(yù)測(cè)的觀測(cè)點(diǎn)處濃度的非零值較多。修正模型的預(yù)測(cè)值的Q-Q 圖更接近y=x 直線,表明APM 模型的預(yù)測(cè)值較K-J模型的更符合觀測(cè)值。
由于傳統(tǒng)煙團(tuán)模型中的擴(kuò)散參數(shù)是根據(jù)大氣穩(wěn)定度及釋放高度的條件選取的,可供選取的值是以理想條件下的實(shí)驗(yàn)結(jié)果為依據(jù),因此,擴(kuò)散參數(shù)在復(fù)雜流動(dòng)和地形條件下,適用性受到了很大局限。為改善擴(kuò)散參數(shù)的適用性,利用當(dāng)時(shí)當(dāng)?shù)氐臍庀髷?shù)據(jù)、源項(xiàng)數(shù)據(jù)及觀測(cè)數(shù)據(jù),以最小二乘法求解大氣擴(kuò)散系數(shù),將符合實(shí)際情況的大氣擴(kuò)散系數(shù)代入模型中,以達(dá)到系數(shù)自適應(yīng)修正的效果。通過計(jì)算,將傳統(tǒng)擴(kuò)散參數(shù)的煙團(tuán)模型與系數(shù)自適應(yīng)修正的煙團(tuán)模型的計(jì)算結(jié)果進(jìn)行比較,后者與觀測(cè)數(shù)據(jù)一致性的統(tǒng)計(jì)結(jié)果得到明顯的提高,證明大氣擴(kuò)散系數(shù)自適應(yīng)修正方法是有效的。
在擬合求解大氣擴(kuò)散系數(shù)時(shí),一個(gè)時(shí)間步長(zhǎng)中風(fēng)速風(fēng)向的影響會(huì)累積到下一個(gè)計(jì)算時(shí)間步長(zhǎng)內(nèi)。縮短計(jì)算時(shí)間步長(zhǎng)可減少風(fēng)速風(fēng)向變化帶來的影響,更準(zhǔn)確地求解大氣擴(kuò)散系數(shù)。但在計(jì)算過程中要求有與計(jì)算時(shí)間步長(zhǎng)相匹配的氣象數(shù)據(jù)。
由計(jì)算過程可發(fā)現(xiàn),大氣擴(kuò)散系數(shù)自適應(yīng)修正方法需要有較為準(zhǔn)確的觀測(cè)值和釋放數(shù)據(jù)。由于對(duì)觀測(cè)值的數(shù)量和精度均有一定的要求,因此,在廠區(qū)內(nèi)如何合理有效地布置監(jiān)測(cè)站是需要考慮的問題。事故發(fā)生時(shí),破口出現(xiàn)的位置及泄漏的數(shù)據(jù)應(yīng)盡可能地詳細(xì),因此安全殼內(nèi)各部件的監(jiān)測(cè)系統(tǒng)需靈敏、準(zhǔn)確,且需發(fā)展源項(xiàng)反演技術(shù),以便獲得準(zhǔn)確的釋放數(shù)據(jù)。
本文針對(duì)拉格朗日煙團(tuán)模型提出了一種大氣擴(kuò)散系數(shù)自適應(yīng)的方法,建立了動(dòng)態(tài)的大氣擴(kuò)散系數(shù)自適應(yīng)修正的拉格朗日煙團(tuán)模型。通過Kincaid實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證表明,大氣擴(kuò)散系數(shù)自適應(yīng)修正方法的煙團(tuán)模型優(yōu)于傳統(tǒng)P-G 曲線擴(kuò)散參數(shù)的煙團(tuán)模型。基于新方法建立的動(dòng)態(tài)大氣擴(kuò)散系數(shù)自適應(yīng)修正煙團(tuán)模型可提高大氣擴(kuò)散模型計(jì)算結(jié)果的準(zhǔn)確性。
[1] 胡二邦,姚仁太,宣儀仁,等.核事故后果評(píng)價(jià)系統(tǒng)的進(jìn)展與比較[J].輻射防護(hù)通訊,2003,23(1):6-12.HU Erbang,YAO Rentai,XUAN Yiren,et al.Comparison and advance of consequences assessment system for nuclear accident[J].Radiation Protection Bulletin,2003,23(1):6-12(in Chinese).
[2] 曲靜原,曹建主,劉磊,等.我國(guó)核應(yīng)急決策支持系統(tǒng)研究開發(fā)的現(xiàn)狀與展望[J].原子能科學(xué)技術(shù),2001,35(3):283-288.QU Jingyuan,CAO Jianzhu,LIU Lei,et al.Status and perspective on the research and development of the Chinese decision support system for nuclear emergency management[J].Atomic Energy Science and Technology,2001,35(3):283-288(in Chinese).
[3] EHRHARDT J,WEIS A.Development of a comprehensive decision support system for nuclear emergency in europe following an accidental release to atmosphere,F(xiàn)ZKA-5772[R].Karlsruhe:Forshungszentrum Karlsruhe,1996.
[4] 蔣維楣.空氣污染氣象學(xué)[M].南京:南京大學(xué)出版社,2003:102-111.
[5] 曲靜原,曹建主.RODOS4.0與RODOS的未來發(fā)展[J].輻射防護(hù),2002,22(4):198-206.QU Jingyuan,CAO Jianzhu.The future development of RODOS4.0and RODOS[J].Radiation Protection,2002,22(4):198-206(in Chinese).
[6] OLESEN H R.Model validation kit-status and outlook[J].International Journal of Environment and Pollution,2000,14(12):65-76.
[7] OLESEN H R.User’s guide to the model validation kit[R].Denmark:National Environment Research Institute,2005.
[8] ELEVELD H.Application of a methodology to validate atmospheric dispersion models[C]∥Proceedings of the 7th International Conference on Harmonization Within Atmospheric Dispersion Modelling for Regulatory Purposes.Italy:[s.n.],2001.
[9] THYKIER-NIELSEN S,DEME S,MIKKELSEN T.Description of the atmospheric dispersion module RIMPUFF[R].Denmark:Risoe National Laboratory,1999.