基于Lagrangian-Eulerian耦合的放射性物质大气扩散模型研究

邰 扬，申世飞

(清华大学 工程物理系,北京 100084)

近年来，我国核电建设持续增长，核电在建规模长期保持全球第一，按目前已形成的核电发展条件预测，2021—2035年每年将核准建设6～8台核电机组，2035年我国核电装机规模将达到1.5～1.8亿kW，核电发电量占比将提升至10%以上，核能在国民经济中将发挥越来越重要的作用[1]。保证核能的安全使用是开发和利用核能的前提，对核电来说，除了通过合理选址进行前期安全规划和建立并维护行之有效的安全运营体系之外，对核事故发生后放射性物质的扩散过程进行模拟并根据结果制定安全预案，也是保障核电安全性的重要方面。以2011年福岛核事故为例，大量的放射性物质持续性地排放到大气中，并在事故发生后的一个多月内随大气运动而大面积扩散，将周围大面积区域都拖入放射性污染的威胁中。因此，针对核能发电设施在核事故发生时放射性物质的扩散过程、时间与空间分布、衰减规律进行模拟研究并用于核事故后的安全疏散等应急决策，具有十分重要的现实意义。

在对放射性物质大气扩散过程的评估中，观测站的观测结果是重要数据来源，然而仅依靠事故后的观测结果会存在覆盖范围有限、无法重现、无法提前进行安全预案等问题，而数值模拟研究则相对具有很大的优势。在常用的放射性物质大气扩散的数值模拟中，基于稳态条件下烟羽在竖直和水平维度上呈现出的高斯分布特征的高斯烟羽模型(Gaussian plume model)是早期最广泛使用的扩散模型，但是高斯烟羽模型存在一些严重的缺陷，例如采取稳态近似而忽略了污染物传播到探测器的时间，且忽略了烟羽之间相互作用关系等假设造成在小尺度、近源区域和大尺度上都存在明显缺陷[2]。因此，常用于模拟放射性物质大气扩散过程的模型大多为拉格朗日模型(Lagrangian model)和欧拉模型(Eulerian model)[3]。

拉格朗日模型将放射性物质在大气中的扩散过程视作粒子团在大气中的运动，利用扩散物理过程中不同途径产生的概率函数来追踪粒子团的轨迹，再将大量粒子团的扩散结果叠加，从而得到放射性物质的时空分布。因为拉格朗日模型的关注点在粒子团本身，所以对大气条件的要求较低，在非稳态和非正态分布条件下依然可以得到较好的准确度，适合用来模拟核事故时近源区域的放射性物质扩散。但是由于粒子团会随着模拟空间的扩大而迅速分散，因此当模拟尺度增大时，为了保证远源条件下的准确度和平滑性，模拟所需要的粒子团数目将大幅上升，计算量也会大幅增大。相比之下，欧拉模型将模拟区域栅格化处理后只考虑放射性物质在相邻栅格间的传输，由于栅格尺寸的灵活性，欧拉模型可以灵活控制计算消耗。但欧拉模型也采取了稳态假设和正态分布假设，且在栅格化的过程中将一般为点排放的排放源近似成一个栅格，因此在非稳态条件和近源区域存在着明显的误差[4]。

因此，将拉格朗日模型和欧拉模型结合起来的耦合模型成为一种可以兼顾拉格朗日模型和欧拉模型优点的模拟方法。国外已有学者对此做了一些研究，如OETTL等[5]建立了可以在非稳态大气条件下使用的拉格朗日-欧拉耦合模型(GRAL模型)，但该模型不包含化学变化和物理衰变模块；HURLEY等[6]开发了用于非静态大气扩散的耦合模型(TAPM模型)，但该模型未考虑气溶胶扩散的情况，且湿沉降部分仅能处理可溶性气体；还有一种广泛使用的ARIA模型[7]，可以支持拉格朗日模式和欧拉模式，但未对两种模型进行耦合。考虑到大多数模型未考虑放射性物质的沉降与衰变特点，且在中尺度和大尺度上普遍存在准确性和易用性上的不足[8]，同时对于耦合判据和耦合过程缺乏深入的探讨，笔者将基于两款成熟的大气扩散模型，即采取拉格朗日模型的FLEXPART[9]和欧拉模型的WRF-CHEM[10]开发一款适用于核事故条件下放射性物质在中尺度和大尺度范围扩散的拉格朗日-欧拉耦合模型，并对耦合方法及其效率进行初步讨论。

1 耦合模型

笔者所使用的拉格朗日-欧拉模型在接近扩散源的位置使用拉格朗日模块FLEXPART[11]进行计算，而在远离扩散源的位置，则根据一定的标准将拉格朗日粒子转化为欧拉浓度利用WRF-CHEM进行计算。本质上来说，耦合模式其实是先用拉格朗日模块对排放源进行预处理，再将结果导入欧拉模块中作为欧拉模式的初始条件，在这个过程中通过耦合判据来判断放射性物质分布是否足够分散，以减少欧拉模式中栅格平均化所带来的误差。

具体来说，耦合模式需要跟踪并分别控制拉格朗日模块中排放源所释放出来的每个粒子团[12]，将粒子团中的某一个物理量作为判据，如扩散时间、标准差等，将其与预设判据基准做比较，一旦符合转化条件则将该粒子团的状态量(输运时间、位置、物质质量等)记录下来，经栅格化处理后转入欧拉模式继续计算。耦合模式的耦合判断过程如图1所示。

图1 拉格朗日-欧拉耦合模式的耦合判断过程

在耦合判据的选取方面，笔者以单个粒子团扩散过程中的弥散程度为判据，即粒子团本身的不确定度S[13]，其表达式为：

(1)

其中，xi、xj分别表示第i个和第j个粒子的位置。本模型在拉格朗日模块中每隔若干时间步长对粒子团的不确定度进行一次检验，若该粒子团的不确定度S达到预设阈值S0，则将该粒子团的状态量[14]转入欧拉模块中继续计算。

2 模型验证与计算

利用搭建的拉格朗日-欧拉耦合模型，对我国东南沿海某核电站进行模拟运算(虚拟算例)。考虑到台风过境期间在固定时间点上较大空间范围的风向、风速分布会比较稳定，故选择台风过境期间模拟核电站泄漏事故，这样不仅可以利用大尺度上的放射性物质扩散与沉降来验证模型的有效性，还可以利用中小尺度上的扩散与沉降来比较耦合模型与单一模型的差异性。

2011年8月上旬第9号超强台风“梅花”曾掠过我国东南沿海并引起严重的气象灾害，假设该核电站于2011年8月4日0点发生泄漏，参考福岛核事故数据，设置泄漏剂量为一次性排放131I为2.3×1013Bq，137Cs为2.3×1012Bq，模拟时间为北京时间2011年8月1日0时至8月11日0时。因气象场为GFS(global forecasting system)数据插值后由WRF(weather research and forecasting)模型模拟得出，所以需提前计算3天的气象数据来建立稳定的气象场。为保证模拟的有效性，拉格朗日模块和欧拉模块均选择相同的微物理参数化方案[15]和干湿沉降模型[16]，并考虑放射性物质随时间的衰变。为了进一步说明模型的有效性，将仅采用欧拉模型进行计算的模拟算例作为参照组。

首先对模拟的气象条件进行了验证，选择了泄漏点200 km范围内7个气象观测站的观测数值进行验证。因超强台风带来的强降水幅度巨大，大部分观测站未能公布降水的分时详细数据，故选取观测数据较为详实的风向和风速数据作为参照。观测站点分别位于栎社、萧山、浦东、嵊泗、虹桥、定海和石浦(以下分别称为观测站1～观测站7)。风速和风向的模拟值与观测值的对比结果分别如图2和图3所示。

图2 风速的模拟值与观测值

图3 风向的模拟值与观测值

从图2可以看出，在大多数观测站，模拟的风速值与观测数据吻合较好，很好地反映了模拟时段内风速随时间的变化规律。

在图3中，y轴为风向与正北方向的夹角，由于对角度来说360°等价于0°，故图上最高点和最低点其实是等价的。从图3可以看出，在大多数观测站，模拟的风向值与观测数据吻合较好，可以较好地反映模拟时段内风向随时间的变化规律。对于放射性物质的沉降，对照组和考察组采用了相同的模拟参数设置，但模拟流程不同，即对照组只采用欧拉模型的WRF-CHEM进行模拟，而考察组先用拉格朗日模型的FLEXPART较精确地模拟点排放源释放的过程，经过耦合模型处理后再转入欧拉模块处理。

从大尺度上验证模型的有效性，笔者选取1 500 km范围的大尺度分别对对照组和考察组进行模拟，得到源排放后12 h、24 h、48 h131I的总沉降量,如图4所示。从图4可以看出，考察组在大尺度上的模拟结果与对照组基本一致，说明耦合模型可以很好地模拟放射性物质在大气中的扩散过程。

图4 大尺度下对照组(左)与考察组(右)131I的总沉降量(单位:Bq/h)

从中尺度上对模型有效性进行验证，将模拟尺度切换到500 km范围内的中尺度上分别对对照组和考察组进行模拟，得到源排放后12 h、24 h、48 h131I的总沉降量，如图5所示。从图5可以看出，在中尺度上对照组和考察组的放射性物质沉降结果表现出较明显的差异，尤其是在相对近源区域(100 km以内)差异更明显。由于强台风带来的高风速和强降水影响，甚至在距离泄漏源200 km的区域依然有很大差异，在实际指导疏散的过程中将给疏散区域的判断带来巨大影响。

图5 中尺度下对照组(左)与考察组(右)131I的总沉降量(单位:Bq/h)

中尺度下对照组和考察组137Cs的总沉降量如图6所示，可见对137Cs的模拟也验证了这一结果，即对照组与考察组差异明显。这一结果表明，在中尺度范围内，耦合模型对于解决欧拉模型计算放射性物质扩散时将排放源均匀化到整个计算格点的假设所带来误差的问题上有着非常好的效果。

图6 中尺度下对照组(左)与考察组(右)137Cs的总沉降量(单位:Bq/h)

3 结论

笔者建立了适用于模拟核事故下放射性物质大气扩散的拉格朗日-欧拉耦合模型，并对模型进行了算例验证与分析。结果表明，在中尺度和大尺度下拉格朗日-欧拉耦合模型均有较好的模拟效果，且耦合模型可以在中尺度范围内较大幅度地减小欧拉模型的初始条件所带来的系统误差。