許嘯峰，張純禹，劉 洋

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

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

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

1 模型和方法

1.1 污染物大氣湍流擴散模型

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

(1)

綜合流體輸運方程和上述湍流擴散方程，可得大氣湍流擴散的控制方程為：

(2)

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

(3)

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

(4)

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

1.2 縮減基有限元方法

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

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

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

(5)

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

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

(6)

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

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

(7)

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

(8)

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

(9)

因為測試函數(shù)要滿足任意性，所以式(9)左右兩邊可同時消除系數(shù)ci，可得：

(10)

將式(10)展開后可得：

(11)

(12)

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

(13)

將式(13)寫為：

(14)

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

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

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

(15)

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

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

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

(16)

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

(17)

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

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

(18)

(19)

其中：

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

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

(28)

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

(29)

(30)

(31)

(32)

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

(33)

(34)

(35)

2) 源項的分解方法

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

(36)

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

(37)

故可得式(28)。

2 算例

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

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

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

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

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

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

3 結(jié)論

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

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

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