許嘯峰，張純禹，劉 洋

(中山大學(xué) 中法核工程與技術(shù)學(xué)院，廣東 珠海 519082)

核安全與人類(lèi)生活息息相關(guān)，2011年的福島核事故[1]由于缺乏良好的應(yīng)急系統(tǒng)，沒(méi)能快速地獲取和預(yù)測(cè)場(chǎng)外污染狀況，嚴(yán)重影響了事故處理的及時(shí)性和有效性。使用計(jì)算機(jī)進(jìn)行數(shù)值模擬是對(duì)污染物擴(kuò)散進(jìn)行預(yù)報(bào)的最常用手段，根據(jù)所求解守恒方程的不同，可將其分為采用解析或半解析模型的模擬方法和直接求解大氣湍流輸運(yùn)方程的模擬方法。前者常采用多種高斯模型，如經(jīng)典高斯模型[2]、高斯煙團(tuán)模型[3]和分段高斯煙羽模型[4]。這類(lèi)方法的優(yōu)點(diǎn)是計(jì)算量小，缺點(diǎn)是為獲得解析解需對(duì)復(fù)雜的邊界條件進(jìn)行高度簡(jiǎn)化，這會(huì)影響模擬的真實(shí)性。而采用有限元、有限體積等數(shù)值方法直接求解大氣湍流輸運(yùn)方程的模擬方法，盡管能較真實(shí)地描述源項(xiàng)條件、氣象條件和地理?xiàng)l件，但由于計(jì)算量大、計(jì)算速度慢而難以在實(shí)際的工業(yè)應(yīng)急研究中應(yīng)用。在求解大氣湍流輸運(yùn)方程進(jìn)行污染物的輸運(yùn)和擴(kuò)散模擬時(shí)，將變化的系數(shù)或條件視為參數(shù)，則可用參數(shù)化的控制方程描述該類(lèi)問(wèn)題。加速求解這類(lèi)問(wèn)題的最直接手段是應(yīng)用各種降階法(ROM)降低模型的階次，大幅減少需求解的未知量數(shù)。在已發(fā)展的多種降階法中，有精度檢驗(yàn)的縮減基有限元方法(CRB-FEM)由于兼?zhèn)涓咝Ш蜏?zhǔn)確的優(yōu)點(diǎn)而在近十年得到快速發(fā)展[5]。

針對(duì)參數(shù)化偏微分方程[6]的求解，縮減基有限元方法的基本思路是預(yù)先求解少量有代表性的經(jīng)典有限元解，然后利用這組解構(gòu)造整個(gè)解空間的基函數(shù)[7]，再采用經(jīng)典有限元法將方程轉(zhuǎn)化為低階次線性方程組進(jìn)行求解。本文采用縮減基有限元方法進(jìn)行模擬求解，并計(jì)算污染物擴(kuò)散分布的情況。

1 模型和方法

1.1 污染物大氣湍流擴(kuò)散模型

氣載污染物在大氣中的遷移擴(kuò)散[8]包括大氣輸送和大氣擴(kuò)散。大氣輸送是氣載污染物在風(fēng)場(chǎng)作用下發(fā)生遷移的過(guò)程，大氣擴(kuò)散是污染物通過(guò)湍流作用與各方向的大氣進(jìn)行混合。借鑒經(jīng)典的湍流理論[9]，大氣湍流擴(kuò)散理論中污染物輸運(yùn)的脈動(dòng)量正比于污染物濃度的梯度，將脈動(dòng)量與平均量進(jìn)行關(guān)聯(lián)處理，以滿(mǎn)足方程組的閉合性。

(1)

綜合流體輸運(yùn)方程和上述湍流擴(kuò)散方程，可得大氣湍流擴(kuò)散的控制方程為：

(2)

式(2)左邊的第2項(xiàng)為輸送項(xiàng)，右邊第1項(xiàng)為湍流擴(kuò)散項(xiàng)，f為污染物源項(xiàng)。不失一般性，本文假定污染源滿(mǎn)足高斯分布：

(3)

式中：t為時(shí)間參數(shù)；Q0為源強(qiáng)；x0、y0、z0為污染源的位置。需要指出的是，σx、σy、σz僅為各方向上的標(biāo)準(zhǔn)差，與傳統(tǒng)的大氣湍流擴(kuò)散模型中的擴(kuò)散系數(shù)無(wú)關(guān)。假定湍動(dòng)擴(kuò)散為各向同性，即Kx=Ky=Kz=k，式(2)可重新寫(xiě)為：

(4)

式中：k為湍動(dòng)擴(kuò)散系數(shù)；b為x、y、z方向速度分量組成的向量。

1.2 縮減基有限元方法

1) 仿射分解與縮減基的構(gòu)造

設(shè)所求解問(wèn)題的有效區(qū)域?yàn)棣?，其邊界?Ω。在求解域上定義Hilbert空間[10]為：

V=V(Ω)={v∈H1(Ω),v|?Ω=0}

(5)

式中：H1為Sobolev空間；v為Hilbert空間的子元素。

將封閉的參數(shù)空間記為G，在該區(qū)域內(nèi)參數(shù)μ∈G，對(duì)應(yīng)的待求量為u(μ):G→V。其中參數(shù)μ可能包含多個(gè)分量，可寫(xiě)為μ=(μ1，μ2，μ3,…,μp)，p為參數(shù)的個(gè)數(shù)。假設(shè)在參數(shù)空間中對(duì)于任一參數(shù)，式(4)左邊項(xiàng)滿(mǎn)足連續(xù)性條件和強(qiáng)制性條件，右邊項(xiàng)滿(mǎn)足連續(xù)性條件，利用Lax-Milgram定理[11]可知方程存在唯一解。該方程的近似解可寫(xiě)為基函數(shù)的線性組合：

(6)

式中：ω為純數(shù)學(xué)意義上的一階參數(shù)；φj(ω)為基函數(shù)，共有M個(gè)；uj為未知的待求系數(shù)。

將式(4)由強(qiáng)形式轉(zhuǎn)化為弱形式為：

(7)

式中，vδ為任意測(cè)試函數(shù)。利用經(jīng)典的Galerkin有限元方法[12]可知，測(cè)試函數(shù)vδ采用與uδ相同的基函數(shù)進(jìn)行展開(kāi)，可表示為：

(8)

式中，ci為任意常數(shù)。將式(6)、(8)代入式(7)可得：

(9)

因?yàn)闇y(cè)試函數(shù)要滿(mǎn)足任意性，所以式(9)左右兩邊可同時(shí)消除系數(shù)ci，可得：

(10)

將式(10)展開(kāi)后可得：

(11)

(12)

將式(6)和式(12)代入式(11)可得：

(13)

將式(13)寫(xiě)為：

(14)

經(jīng)典有限元方法一般采用形式簡(jiǎn)單但通用的多項(xiàng)式作為基函數(shù)，這導(dǎo)致只有采用數(shù)量龐大的基函數(shù)才能較好地逼近復(fù)雜問(wèn)題的真實(shí)解。

為了降低求解問(wèn)題的自由度，縮減基有限元方法采用典型參數(shù)對(duì)應(yīng)的經(jīng)典有限元解構(gòu)造近似解空間，即：

VRB=span{uδ(μ1),uδ(μ2),…,uδ(μN(yùn))}

(15)

式中，N為構(gòu)造近似解空間里的縮減基數(shù)目。

為了使新的解空間里各基函數(shù)相互線性無(wú)關(guān)，可利用Gram-Schemidt正交化方法來(lái)對(duì)這些基函數(shù)進(jìn)行轉(zhuǎn)換處理：

VRB=span{ξ1,ξ2,…,ξN}

(16)

則采用這種縮減基的近似解可表示為：

(17)

式中，ξj(ω)為縮減基的基函數(shù)，各基函數(shù)之間相互線性無(wú)關(guān)。

采用與經(jīng)典有限元方法同樣的步驟，可將偏微分方程轉(zhuǎn)化為線性方程組：

(18)

(19)

其中：

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

源項(xiàng)中的參數(shù)和空間坐標(biāo)可采用經(jīng)驗(yàn)插值法(EIM)[13]進(jìn)行近似但足夠精確的仿射分解：

(28)

Qf為達(dá)到精度要求所需的疊加次數(shù)，縮減基所采用的基函數(shù)是經(jīng)典有限元的解，二者之間存在轉(zhuǎn)換映射關(guān)系：

(29)

(30)

(31)

(32)

將式(19)、(27)～(28)代入式(30)～(32)可得：

(33)

(34)

(35)

2) 源項(xiàng)的分解方法

以式(3)和(13)為基礎(chǔ)，利用EIM可將源項(xiàng)寫(xiě)為如下形式：

(36)

式中：e為達(dá)到要求精度時(shí)所需的疊加次數(shù)；gq(x,y,z)為參數(shù)無(wú)關(guān)項(xiàng)；hq(μ)為參數(shù)有關(guān)項(xiàng)；rq(x0,y0,z0)為源項(xiàng)位置參數(shù)相關(guān)項(xiàng)；sq(σx,σy,σz)為標(biāo)準(zhǔn)差參數(shù)相關(guān)項(xiàng)。由強(qiáng)形式轉(zhuǎn)換為弱形式后可得：

(37)

故可得式(28)。

2 算例

選取大亞灣核電站為事故模擬點(diǎn)，發(fā)生核反應(yīng)堆事故時(shí)，近地面的空氣里包含有多種放射性核素。131I、137Cs等為釋放份額相對(duì)較高的核素，以131I和137Cs污染物為源項(xiàng)，具體源項(xiàng)數(shù)據(jù)列于表1。

表1 源項(xiàng)數(shù)據(jù)Table 1 Source data

取以污染物泄漏點(diǎn)為中心且邊長(zhǎng)為60 km的正方形區(qū)域?yàn)槟M區(qū)域，每km設(shè)置3個(gè)網(wǎng)格點(diǎn)，時(shí)間步長(zhǎng)為0.001 h，總計(jì)算步數(shù)為3 000，計(jì)算核事故發(fā)生后3 h內(nèi)污染物在空氣中的濃度分布情況。各參數(shù)的范圍和用于驗(yàn)證的參數(shù)列于表2。

表2 高斯模型參數(shù)Table 2 Parameter of Gauss model

算出近地面污染物131I、137Cs在不同條件下的濃度分布，每個(gè)狀態(tài)時(shí)間間隔為0.5 h，如圖1、2所示。

與非縮減基模型結(jié)果對(duì)比，需90個(gè)縮減基可使得計(jì)算濃度相對(duì)誤差低于10-3。其中，離線階段耗時(shí)20 h左右，在線階段計(jì)算時(shí)間僅需4 s左右(Intel i5-2520M CPU@2.5 GHz，包括文件讀寫(xiě)的時(shí)間)，若不考慮文件讀寫(xiě)，在線階段僅需2 s左右。根據(jù)不同的環(huán)境條件更改參數(shù)文件里的可變參數(shù)即可進(jìn)行在線階段的計(jì)算，極大提高了計(jì)算效率。

3 結(jié)論

運(yùn)用縮減基有限元法求解污染物在大氣中的運(yùn)輸擴(kuò)散模型是一種高效精確的求解方法，利用仿射分解將所求解的問(wèn)題分解為參數(shù)相關(guān)和參數(shù)無(wú)關(guān)部分，并因此劃分離線階段和在線階段。在線階段計(jì)算速度快，滿(mǎn)足實(shí)時(shí)性要求，因此該方法可運(yùn)用于核事故污染物大氣擴(kuò)散的近實(shí)時(shí)模擬。

圖1 131I污染物濃度分布圖Fig.1 Distribution of 131I pollutant

圖2 137Cs污染物濃度分布Fig.2 Distribution of 137Cs pollutant