核事故污染物大气扩散的三维近实时模拟方法研究

时间：2024-07-28

许啸峰，张纯禹，刘洋

(中山大学中法核工程与技术学院，广东珠海 519082)

核安全与人类生活息息相关，2011年的福岛核事故[1]由于缺乏良好的应急系统，没能快速地获取和预测场外污染状况，严重影响了事故处理的及时性和有效性。使用计算机进行数值模拟是对污染物扩散进行预报的最常用手段，根据所求解守恒方程的不同，可将其分为采用解析或半解析模型的模拟方法和直接求解大气湍流输运方程的模拟方法。前者常采用多种高斯模型，如经典高斯模型[2]、高斯烟团模型[3]和分段高斯烟羽模型[4]。这类方法的优点是计算量小，缺点是为获得解析解需对复杂的边界条件进行高度简化，这会影响模拟的真实性。而采用有限元、有限体积等数值方法直接求解大气湍流输运方程的模拟方法，尽管能较真实地描述源项条件、气象条件和地理条件，但由于计算量大、计算速度慢而难以在实际的工业应急研究中应用。在求解大气湍流输运方程进行污染物的输运和扩散模拟时，将变化的系数或条件视为参数，则可用参数化的控制方程描述该类问题。加速求解这类问题的最直接手段是应用各种降阶法(ROM)降低模型的阶次，大幅减少需求解的未知量数。在已发展的多种降阶法中，有精度检验的缩减基有限元方法(CRB-FEM)由于兼备高效和准确的优点而在近十年得到快速发展[5]。

针对参数化偏微分方程[6]的求解，缩减基有限元方法的基本思路是预先求解少量有代表性的经典有限元解，然后利用这组解构造整个解空间的基函数[7]，再采用经典有限元法将方程转化为低阶次线性方程组进行求解。本文采用缩减基有限元方法进行模拟求解，并计算污染物扩散分布的情况。

1 模型和方法

1.1 污染物大气湍流扩散模型

气载污染物在大气中的迁移扩散[8]包括大气输送和大气扩散。大气输送是气载污染物在风场作用下发生迁移的过程，大气扩散是污染物通过湍流作用与各方向的大气进行混合。借鉴经典的湍流理论[9]，大气湍流扩散理论中污染物输运的脉动量正比于污染物浓度的梯度，将脉动量与平均量进行关联处理，以满足方程组的闭合性。

(1)

综合流体输运方程和上述湍流扩散方程，可得大气湍流扩散的控制方程为：

(2)

式(2)左边的第2项为输送项，右边第1项为湍流扩散项，f为污染物源项。不失一般性，本文假定污染源满足高斯分布：

(3)

式中：t为时间参数；Q0为源强；x0、y0、z0为污染源的位置。需要指出的是，σx、σy、σz仅为各方向上的标准差，与传统的大气湍流扩散模型中的扩散系数无关。假定湍动扩散为各向同性，即Kx=Ky=Kz=k，式(2)可重新写为：

(4)

式中：k为湍动扩散系数；b为x、y、z方向速度分量组成的向量。

1.2 缩减基有限元方法

1) 仿射分解与缩减基的构造

设所求解问题的有效区域为Ω，其边界为∂Ω。在求解域上定义Hilbert空间[10]为：

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

(5)

式中：H1为Sobolev空间；v为Hilbert空间的子元素。

将封闭的参数空间记为G，在该区域内参数μ∈G，对应的待求量为u(μ):G→V。其中参数μ可能包含多个分量，可写为μ=(μ1，μ2，μ3,…,μp)，p为参数的个数。假设在参数空间中对于任一参数，式(4)左边项满足连续性条件和强制性条件，右边项满足连续性条件，利用Lax-Milgram定理[11]可知方程存在唯一解。该方程的近似解可写为基函数的线性组合：

(6)

式中：ω为纯数学意义上的一阶参数；φj(ω)为基函数，共有M个；uj为未知的待求系数。

将式(4)由强形式转化为弱形式为：

(7)

式中，vδ为任意测试函数。利用经典的Galerkin有限元方法[12]可知，测试函数vδ采用与uδ相同的基函数进行展开，可表示为：

(8)

式中，ci为任意常数。将式(6)、(8)代入式(7)可得：

(9)

因为测试函数要满足任意性，所以式(9)左右两边可同时消除系数ci，可得：

(10)

将式(10)展开后可得：

(11)

(12)

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

(13)

将式(13)写为：

(14)

经典有限元方法一般采用形式简单但通用的多项式作为基函数，这导致只有采用数量庞大的基函数才能较好地逼近复杂问题的真实解。

为了降低求解问题的自由度，缩减基有限元方法采用典型参数对应的经典有限元解构造近似解空间，即：

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

(15)

式中，N为构造近似解空间里的缩减基数目。

为了使新的解空间里各基函数相互线性无关，可利用Gram-Schemidt正交化方法来对这些基函数进行转换处理：

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

(16)

则采用这种缩减基的近似解可表示为：

(17)

式中，ξj(ω)为缩减基的基函数，各基函数之间相互线性无关。

采用与经典有限元方法同样的步骤，可将偏微分方程转化为线性方程组：

(18)

(19)

其中：

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

源项中的参数和空间坐标可采用经验插值法(EIM)[13]进行近似但足够精确的仿射分解：

(28)

Qf为达到精度要求所需的叠加次数，缩减基所采用的基函数是经典有限元的解，二者之间存在转换映射关系：

(29)

(30)

(31)

(32)

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

(33)

(34)

(35)

2) 源项的分解方法

以式(3)和(13)为基础，利用EIM可将源项写为如下形式：

(36)

式中：e为达到要求精度时所需的叠加次数；gq(x,y,z)为参数无关项；hq(μ)为参数有关项；rq(x0,y0,z0)为源项位置参数相关项；sq(σx,σy,σz)为标准差参数相关项。由强形式转换为弱形式后可得：

(37)

故可得式(28)。

2 算例

选取大亚湾核电站为事故模拟点，发生核反应堆事故时，近地面的空气里包含有多种放射性核素。131I、137Cs等为释放份额相对较高的核素，以131I和137Cs污染物为源项，具体源项数据列于表1。

表1 源项数据Table 1 Source data

取以污染物泄漏点为中心且边长为60 km的正方形区域为模拟区域，每km设置3个网格点，时间步长为0.001 h，总计算步数为3 000，计算核事故发生后3 h内污染物在空气中的浓度分布情况。各参数的范围和用于验证的参数列于表2。

表2 高斯模型参数Table 2 Parameter of Gauss model

算出近地面污染物131I、137Cs在不同条件下的浓度分布，每个状态时间间隔为0.5 h，如图1、2所示。

与非缩减基模型结果对比，需90个缩减基可使得计算浓度相对误差低于10-3。其中，离线阶段耗时20 h左右，在线阶段计算时间仅需4 s左右(Intel i5-2520M CPU@2.5 GHz，包括文件读写的时间)，若不考虑文件读写，在线阶段仅需2 s左右。根据不同的环境条件更改参数文件里的可变参数即可进行在线阶段的计算，极大提高了计算效率。