CN107145698A - 核与辐射突发事件后果大尺度模拟方法及模拟*** - Google Patents

Abstract

本发明涉及一种核与辐射突发事件后果大尺度模拟方法及***，该模拟方法包括对核爆源项的估算、对气载放射性核素长距离迁移轨迹的计算、对扩散的计算、对剂量的估算；模拟***包括用于对核爆源项估算的源项估算模块、用于对气载放射性核素长距离迁移轨迹计算的轨迹计算模块、用于对气载放射性核素长距离迁移拉格朗日粒子扩散计算的粒子扩散计算模块、用于对气载放射性物质剂量估算的剂量估算模块。本发明的模拟方法和***，能够在全球范围内核与辐射突发事件下，对气载放射性物质的迁移、扩散及其环境安全后果进行评估，为决策人员面对类似事件发生时，能够快速、科学、有效地提出决策建议与方案，最终达到保护我国公众与环境安全的目的。

Description

核与辐射突发事件后果大尺度模拟方法及模拟***

技术领域

本发明属于核技术领域，具体涉及一种核与辐射突发事件后果大尺度模拟方法及模拟***，用于对全球范围的核与辐射突发事件下，对气载放射性物质的迁移、扩散及其环境安全后果的评估。

背景技术

一直以来，世界各国均非常重视和关注放射性核素大尺度大气迁移问题，不仅是军事、国防方面的需要，而且是核事故应急、环境危害问题的需要。放射性物质在大气中的长距离迁移和后果评价与分析技术涉及众多的科学领域，属于多学科交叉的技术，模拟技术和评价方法也在不断地更新与发展，因此，包括美国、欧盟、日本在内的发达国家、地区也在不断加强该技术的研发。另外，放射性核素大尺度迁移数值模拟技术已经作为一种国家的技术资源，随着相关领域技术的不断发展与更新，从战略角度考虑，一个国家必须有技术储备和可持续发展的能力。当前国际形势变化多端，需要我们不断建立与完善大气输运模拟技术，并提供技术支撑。

近些年来，在朝核试验、日本福岛核电泄漏等事件接连发生的情况下，建立并不断完善一套可应对我国周边数百至数千公里范围内气载放射性物质释放情形应急响应***的必要性已经成为共识。本发明提供的模拟***及模拟方法，将为决策人员能够面对类似事件，快速、科学、有效地提出决策建议与方案，最终达到保护我国公众与环境安全的目标。

为此，中国辐射防护研究院自主开发、建立了核与辐射突发事件后果大尺度模拟***，包括放射性核素大尺度大气迁移数值模拟技术，可应对全球范围的核与辐射突发事件下，气载放射性物质的迁移、扩散及其环境安全后果的评价***。

发明内容

针对现有技术中存在的问题，本发明提供一种核与辐射突发事件后果大尺度模拟方法，能够在核与突发事件后对核爆源项、气载放射性核素长距离迁移轨迹、气载放射性核素长距离迁移拉格朗日粒子扩散及气载放射性物质剂量进行估算，为决策人员面对类似事件发生时，能够快速、科学、有效地提出决策建议与方案，最终达到保护我国公众与环境安全的目的。

为达到以上目的，本发明采用的技术方案是：提供一种核与辐射突发事件后果大尺度模拟方法，包括如下步骤：

1)对模拟参数的准备；即对天气预报数据及模拟方向、模拟起止时间、释放起止时间、轨迹时间间隔、释放位置、模拟区域范围的准备；

2)对核与辐射事件后果进行模拟计算，包括：

对核爆源项估算；即通过TNT当量对***烟云的高度和半径的估算；根据裂变方程，对核***后裂变产物、放射性活度、份额随时间变化的估算；

对气载放射性核素长距离迁移轨迹的计算；即通过解二阶差分形式的轨迹运动方程，确定空气质点的运行轨迹；

对气载放射性核素长距离迁移拉格朗日粒子扩散的计算；即通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；通过统计粒子数、估算网格浓度；

对气载放射性物质剂量的估算；即对烟云外照射剂量估算；对地面沉积外照射剂量估算；对烟云吸入内照射剂量估算；对总剂量水平估算。

进一步，步骤2)中，所述***烟云高度和半径通过下列公式进行估算：

其中：W为***当量；LW＝log₁₀(W/kT)；

进一步，步骤2)中，所述裂变方程为：

式中：i—单线链中核素的序号；λ_i为第i个核素的衰变常数；y_i—单线链中i核素的绝对独立产额；p(t)—单位时间内核燃料的裂变次数；A_i—同量异位链中i核素的放射性。

进一步，步骤2)中，所述二阶差分形式的轨迹运动方程如下：

式中，x_i(t)＝[x(t),y(t),z(t)]为t时刻质点的坐标(x,y,z)；为质点在位置处x_i(t)的风速；分别为三维风分量。

进一步，如果计算时间步长Δt内气团运动加速度恒定，那么运动轨迹的坐标可由下列方程表述：

当给定初始条件v₀[x_0i(t),t],即u(x₀,y₀,z₀,t)、v(x₀,y₀,z₀,t)和w(x₀,y₀,z₀,t),通过迭代计算确定空气气团运动的轨迹；x₀,y₀,z₀分别为笛卡尔坐标系中的初始坐标。

进一步，步骤2)中，所述积分粒子运动的拉格朗日方程为：

x_i(t+Δt)＝x_i(t)+v_i(x_i(t)，t)Δt+v′_i(x_i(t)，t)Δt+v″_i(x_i(t)，t)Δt i＝1，2，3

其中，x_i为粒子的三维坐标分量；v_i为平均风速分量v′_i为湍流脉动速度分量(u′,v′,w′)；v″_i为中尺度风脉动速度分量(u"，v"，w″)；t为时间序列；Δt为时间步长。

进一步，步骤2)中，所述网格浓度通过下列公式计算：

其中，Q，释放的放射性活度；N，释放的粒子总数；T，所有粒子扩散时间的总和；N_ik，第i个粒子在k网格中穿行时间内的时间步数；Δt，时间步长；V_k，网格k的体积。

进一步，步骤2)中，所述烟云外照射剂量通过下列公式计算：

其中：t₁、t₂代表某时段的起始和结束时刻；org代表某种器官；DC为在网格(x,y)处，器官org在t₁至t₂时段内所受到的烟云外照射剂量；AC为在网格(x,y)处，在t₁至t₂时段内核素nuc的积分空气浓度；DFC为核素nuc对器官org的半无限烟云剂量转换因子；SFC为烟云外照射的屏蔽因子，由t₁至t₂时段公众所采取的防护行动的种类决定。

所述地面沉积外照射剂量通过下列公式计算：

其中：DG为在网格(x,y)处，器官org在t₁至t₂时段内所受到的地面沉积外照射剂量；GC为在网格(x,y)处，t时刻核素nuc的地面沉积浓度；GRF为核素nuc对器官org的地面沉积外照射剂量转换因子；SFG为地面沉积外照射的屏蔽因子，由t₁至t₂时段公众所采取的防护行动的种类决定。

所述烟云吸入内照射剂量通过下列公式计算：

其中：DI为在网格(x,y)处，器官org在t₁至t₂时段内吸入放射性物质所造成的待积吸入剂量；AC为在网格(x,y)处，在t₁至t₂时段内核素nuc的积分空气浓度；IF为核素nuc对器官org的待积吸入内照射剂量转换因子；BR为公众在t₁至t₂时段内的呼吸率；IPF为公众在t₁至t₂时段内的吸入防护因子。

进一步，步骤2)中，所述总剂量水平的估算，是将烟云外照射剂量估算的结果、地面沉积外照射剂量估算的结果以及烟云吸入内照射剂量估算的结果进行累计。

本发明还提供一种核与辐射突发事件后果大尺度模拟***，包括模拟参数准备模块、源项估算模块、轨迹计算模块、粒子扩散计算模块、剂量计算模块；

模拟参数准备模块：用于对天气预报数据及轨迹方向、时间长度、起止时间、轨迹时间间隔、轨迹初始坐标、模拟区域范围的准备；

源项估算模块：用于对核爆源项的估算，包括对***烟云的高度和半径的估算；根据裂变方程，对核***后裂变产物、放射性活度、份额随时间变化的估算；

源项估算模块：用于对气载放射性核素长距离迁移轨迹的计算，即通过解二阶差分形式的轨迹运动方程，确定空气质点的运行轨迹；

粒子扩散计算模块：用于对气载放射性核素长距离迁移拉格朗日粒子扩散的计算，即通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；统计粒子数、计算网格浓度。

剂量计算模块：用于对气载放射性物质剂量的估算，包括对烟云外照射剂量估算、对地面沉积外照射剂量估算、对烟云吸入内照射剂量估算、对总剂量水平估算。

本发明的有益技术效果在于：(1)本发明的模拟方法，能够对核爆源项、对气载放射性核素长距离迁移轨迹、对气载放射性核素长距离迁移拉格朗日粒子扩散及对气载放射性物质剂量进行估算，从而在核与辐射突发事件下，提供源项数据、轨迹结果数据、扩散模拟结果数据及剂量数据，便于为决策人员采取相应的应急响应决策提供技术支持，保护我国公众和环境的安全。

(2)本发明的模拟***，具有可靠性、开放性、兼容性、可扩展性、可维护性及可管理性。

附图说明

图1是本发明核与辐射突发事件后果大尺度模拟方法的流程图；

图2是本发明源项估算方法中核爆当量与烟云高度和半径的关系图；

图3是本发明源项估算方法中各同位素随时间剂量贡献份额的变化曲线；

图4是源项估算模块的结构图；

图5是本发明气载放射性核素长距离迁移轨迹计算方法的流程图；

图6是本发明气载放射性核素长距离迁移拉格朗日粒子扩散方法的流程图。

具体实施方式

下面结合附图，对本发明的具体实施方式作进一步详细的描述。

如图1所示，是本发明提供的核与辐射突发事件后果大尺度模拟方法，包括如下步骤：

1)对模拟参数准备；即对天气预报数据及模拟方向、模拟起止时间、释放起止时间、轨迹时间间隔、释放位置、模拟区域范围的准备；

2)对核与辐射事件后果进行模拟计算，包括对核爆源项估算、对气载放射性核素长距离迁移轨迹的计算、对气载放射性核素长距离迁移拉格朗日粒子扩散计算及对气载放射性物质剂量的估算；

3)将上述各项估算数值生成图表，运行结束。

其中，在步骤1)中，如果模拟参数准备不正常，需返回用户，重新准备。

下面分别对步骤2)中的核爆源项的估算、气载放射性核素长距离迁移轨迹的计算、气载放射性核素长距离迁移拉格朗日粒子扩散的计算及气载放射性物质剂量的估算进行说明：

(一)源项估算

1、源项估算方法

本发明主要针对地面或空中核***，根据***的TNT当量给出确定***烟云抬升高度、***生成的稳定态(零后0.5天)的放射性核素及其活度，主要考虑粒径小于100μm部分，该部分物质会经长距离迁移的。

1)***烟云几何尺寸的确定

空中或地面***后约1秒钟时，火球达到最大尺寸。短时间“升腾”后，炽热的气泡开始加速上升，几秒钟内达到最大上升速度(约300英尺/秒，即91.44m/s)。只要由辐射、周围空气的卷入和绝热膨胀仍未使气团冷却到周围空气的温度，气团就一直上升。

在上升期间内，烟云发展成熟知的蘑菇状。蘑菇顶初始由强烈的旋转旋涡构成。随云不断上升旋转逐渐减缓。对标称弹来说，在烟云稳定后蘑菇顶厚度约达10000英尺(3048m)。这种情况一般出现在碎片云进入同温层底部(在中纬度地区约为30000—40000英尺，即9144—12192m)时，但如果可用的初始能量较少时，也可能出现在某一较低高度。由***到稳定的时间约为5—15分钟。

烟云低部(蘑菇茎)特征与***类型有关，高空***时，实际上没有茎，因为地面物质几乎未被卷到云中。在这种情况下，茎的可见部分虽然也有大量的放射性，但主要是在气流上升时形成的水滴云所组成的。低空***时，冲击波和热辐射的加热作用，把大量泥土卷到上升气流中。云中心部分密度较低，这表明蘑菇中有一旋转环。细茎与蘑菇顶分开，这也是空中***的特征。低部含尘云的高度约只占茎高的三分之一，其余可见部分主要是水滴云。但是在地面***时或塔上***时，由下沉的尘埃和碎片构成的连续熏烟柱，由地表面直达蘑菇顶，蘑菇和茎的构形取决于气象学参数。稳定的气层虽不能充分地阻止上升，但可以引起更迅速的减缓，而使云的颗粒部分拉宽。更重要的气象因素是：风的切变起着“撕裂”作用。并在有些情况下使云的某些部分整片断裂。根据核爆试验，拟合出***当量与源高和半径的关系。

核武器***蘑菇云最终抬升高度H与半径R是其***当量W的函数：

其中，LW＝log₁₀(W/kT)。图2给出了上述关系式的曲线，H为烟云高度、R为烟云半径。对于低当量的核爆可以近似表示为：

H≈3000m(W/kT)^1/3

对于LW<1.5，例如W<32kT采用上述公式内差。上述公式可计算近似到50MT的当量(“Tsar Bomba”试验，1961)。可以外推到100MT当量的***。

2)核***后源项考虑的核素

核***后裂变产物、剩余的裂变物质和结构材料在高温火球中气化，形成具有很强放射性的气溶胶颗粒物。较大的放射性气溶胶颗粒物在距爆心几百公里范围内的重力沉降，称为局部放射性沉降；注入对流层的较小放射性气溶胶颗粒物主要在同一半球范围内围绕地球沉降，称为对流层放射性沉降；注入平流层的微小放射性气溶胶颗粒物形成的世界范围内的沉降，称为全球性放射性沉降。核***产生的引起外照射的主要核素有：¹³⁷Cs和短寿命核素⁹⁵Zr、¹⁰⁶Ru、¹⁴⁰Ba、¹⁴⁴Ce、¹⁰³Ru和¹⁴¹Ce等；引起外照射的主要核素有：¹⁴C、¹³⁷Cs、⁹⁰Sr、¹⁰⁶Ru、¹⁴⁴Ce、³H、¹³¹I、^239,240,241Pu、⁵⁵Fe、²⁴¹Am和⁸⁹Sr。

根据核素随时间迁移产生的剂量份额的不同，以²³⁵U的裂变产物为例，最常见的裂变产物有同位素碘、铯、锶、氙和钡，随时间的迁移，其剂量贡献的份额发生变化。如图3所示，为各同位素随时间的剂量贡献份额的变化曲线，因此在大尺度大气迁移扩散中核爆中核素的选择还应考虑依据中长时间尺度较大范围内最终导致的剂量贡献。根据文献“Hunter and Ballon，²³⁵U的瞬时慢中子裂变，ADC-65,P73”中描述，考虑到²³⁵U裂化核的相对放射性，主要考虑下述核素，其中各同位素的相对放射性份额见表1，该结果为零后冷却时间为2.5天及以后若干时间的结果。

表1 ²³⁵U裂化核的相对放射性产额(％)

对于国内有关单位提供的资料，按0.5d冷却时间结果给出释放源项，该结果中列有300多种核素，考虑以下两方面因素：释放活度占总贡献的大小和剂量转换因子的大小，排选出72种核素，本***按72中核素计算。

3)核爆放射性烟云的粒径特征

核爆后，裂变产物、未裂变的核材料与被火球加热汽化的残留物产生的粒径在10nm～20μm直径的精细微粒，这些粒子将被迅速带入同温层，尤其是***当量超过10kt时，这些辐射微粒在大气中的扩散将是全球尺度的。早期核爆产生的放射性微粒对于全球尺度的影响很少被关注，美国原子能委员会假定微粒将是在全球范围内扩散的，受大气中风的影响而扩散，随着几周、月或年逐渐沉降到地球表面。对于⁹⁰Sr与¹³⁷Cs等长半衰期核素的潜在累积影响是长期的，但危害要小于局地的放射性粒子的影响。

4)放射性核素活度的估算

(1)核***释放能量

重核(²³⁵U，or²³⁹Pu)+n→Σfpi+υ+2-3n+200MeV；1kt TNT释放能量为1×10¹²Cal＝4.2×10¹⁹erg；这时，核***释放能量中，以瞬发能量释放的180MeV/次fp；因此，1kt TNT当量的核***，共有1.45×10²³次核裂变；共计56克²³⁵U全部裂变；即1kt TNT当量＝1.45×10²³次裂变。

(2)核***生成的放射性核素活度

重核裂变后，裂变碎片变成裂变产物，裂变产物的质量数A在72～165范围内，这些产物是不稳定的β放射性核素，构成了同量异位链。平均而言，大部分初始裂变产物在达到稳定状态之前要经过3次或4次放射性蜕变。

裂变生成的放射性核素活度A_io为A_io＝λ_iN_iD＝λ_in·η_i，式中：λ_i为第i个核素的衰变常数；n为总裂变次数；η_i是第i种核素的裂变产额；A_io是裂变生成的初始活度。

裂变产额η分为绝对独立产额，相对独立产额和链产额，上述表1给出不同冷却时间的相对产额。对于***5～15min后各核素的裂变产额，一方面，今后还需要进一步查阅文献获得；另一方面，如果未能获得该数据，则可根据目前获得的各核素裂变产额随冷却时间的变化关系外推。

(3)瞬发裂变产物的放射性特性

核***当量乃是原子核裂变或聚变时释放的能量。1克梯恩梯(TNT)******释放出4.18×10¹⁰尔格(1000卡)的能量。如前所述，当一个铀核、钚核或者钍核裂变时释放出约200兆电子伏的能量，这相当于3.2×10^-4尔格。

考虑到在瞬发裂变过程中仅有大约90％的能量释放出来，其余10％的能量在裂变以后伴随裂变产物的β和γ蜕变而滞后释放。可以计算出与1kt TNT当量的核***相对应的裂变物质中发生裂变的核数为

这个发生裂变的核数大约相当于60克铀或钚。严格说来，梯恩梯当量和裂变核数之间的关系式仅仅对***才是正确的，它所释放的能量完全是核燃料裂变链式反应的结果。

在热核***时，利用能量为14兆电子伏的聚变中子引起²³⁸U的瞬发裂变，梯恩梯当量由三个过程的能量构成：²³⁵U、²³⁹Pu或²³³U的“引爆”裂变，热核装置中轻核聚变和14兆电子伏的快中子引起的²³⁸U裂变。这样，在其中发生1.45×10²³次²³⁸U裂变的热核装置的额定梯恩梯当量要比1kt TNT***时释放的能量大。

在不同的热核装置内，由裂变反应释放的能量大不相同，这与产品的结构特点有关。这里给出的所有放射特性均对应于1kt TNT当量裂变，即1.45×10²³个发生裂变的核数。这里没有考虑裂变的瞬发γ光子以及产品构件和周围介质(空气、土壤)的原子俘获中子而产生的γ辐射。关于单线同量异位链瞬发裂变产物的放射性A的积累过程，在裂变时，j核素的放射性A_j由下式确定：

这里λ_j和y_j分别是裂变产物j的衰变常数和绝对独立产额。

为了求得单线同量异位链所有核素在t>0时的放射性活度，需要以上述值为初始条件，求解原始微分方程组。例如同量异位链第一、第二和第三核素的放射性积累可分别用下列各解析式来描述：

式中的e表示自然对数；

对于瞬发裂变产物混合物γ辐射的其他放射特性可用上述给出的公式计算。

2、数据处理

***根据***核爆源项特点，由用户选择核弹***位置、核弹***当量或监测到的核爆产生的地震震级与核弹的类型生成核***源项所需要的信息，根据这些信息生成源项估算模块，如图4所示。

1)核爆源项放射性核素活度处理方法

对于铀或钚等具有几个核素的单线同量异位链(A₁→A₂→…A_i→…A_n→A)，描述其中各核素的比放射性随时间(t)的变化的原始微分方程组如下：

式中：i——单线链中核素的序号，j≤i≤n；y_i——链中i核素的绝对独立产额；p(t)——单位时间内核燃料的裂变次数；A_i——同量异位链中i核素的放射性。

首先研究最普遍的情况：p(t)＝p₀＝常数。对于反应堆而言，这意味着热功率恒定。由于初始核燃料是随时间而消耗的，保持功率不变就需要提高中子通量密度。在下列初始条件下试求微分方程(1)的解：

A₁(0)＝A₁₀，A₂(0)＝A₂₀，...，A_i(0)＝A_i0，...，A_n(0)＝A_n0 (2)

同量异位链第一个核素的放射性为：

对i≠1的解将以下式表示：

式中：1≤j≤i，而系数b_i和c_i待定。将(4)式代入方程(1)，并进行相应的移项整理，得到：

为了使这个等式对所有t值均成立，其必要和充分的条件是各指数项的系数均等于零。由此可导出C_ij和b_i的递推关系式：

b_i＝b_i-1+y_ip₀ (7)

对于i＝1，从(7)式可得到b₁＝y₁p₀；而在一般情况下，则为

由于j≤i，要求有个C_ij系数。关系式(6)只能定出个系数，余下的n个系数可从t＝0时的初始条件求得：

于是，公式(4)全部系数的计算式如下：

(b)C₁₁＝A₁₀-b₁；

(d)C₂₂＝A₂₀-b₂-C₂₁；

等等，直至C_nn。

同量异位链中各核素放射性的最终递推关系式可写成如下形式：

在热核***时，由热中子、裂变谱中子和14兆电子伏快中子引起的²³⁵U、²³⁸U和²³⁹Pu的瞬发裂变产物的放射性Qj(t)，见附录表2。

表2 由热中子引起的²³⁵U瞬发裂变产物的放射性Qj(t)，居里/千吨梯恩梯当量裂变或居里/1.45×10²³次裂变

(二)气载放射性核素长距离迁移轨迹计算

如图5所示，是本发明提供的气载放射性核素长距离迁移轨迹计算方法，包括如下步骤：

(1)获取并存储天气预报数据、地形数据和模式参数。

天气预报数据包括两个水平风分量u、v、一个垂直风分量w、温度T三维场、地表压力p_s；地形数据包括地球表面经度和纬度；模式参数包括轨迹类型(等压轨迹)、坐标系类型(笛卡尔坐标)、轨迹方向、时间长度、起止时间、轨迹时间间隔、轨迹初始坐标、模拟区域范围。

(2)进行坐标变换。该坐标变换是将等压层气象数据转换为随地笛卡尔坐标。

(3)根据轨迹的运动方程，确定空气质点的运动轨迹。

一般地，每一条计算的大气轨迹代表一个空气团运动的时空轨迹。我们将这些轨迹作为是对某种物质扩散烟云平均运动的估计。

TraModel轨迹模式通过解二阶差分形式的轨迹运动方程确定空气质点(或气团)运动的轨迹。

假设一个气团为无限小，看作是一个质点，那么质点在指定的一段时间内走过的路径称为该质点的轨迹。轨迹的方程为：

其中，x_i(t)＝[x(t),y(t),z(t)]为t时刻质点(或气团)的坐标(x,y,z)；

为质点(或气团)在位置处x_i(t)的风速，x为东西向坐标、y为南北向坐标、z为垂向坐标。

对于所研究的问题需要进行离散化处理，因此，需用有限差分近似的形式解方程(14)。这意味着，数值计算的精度取决于方程(14)泰勒展开式二阶或高阶项的取舍。如果轨迹计算的积分时间步长很短，则一阶形式的差分方程的计算精度足以满足要求。综合考虑计算机资源和计算精度，同时考虑到气象场离散的时间间隔，二阶形式的差分方程是较理想的。这样，当假定在计算时间步长Δt(或风观测周期)内气团运动加速度恒定，则其运动轨迹的坐标可由下列方程表述：

对于方程(15)求解，当给定初始条件v₀[x_0i(t),t](即u(x₀,y₀,z₀,t)、v(x₀,y₀,z₀,t)和w(x₀,y₀,z₀,t))，通过迭代计算就可以确定空气气团运动的轨迹，x₀、y₀、z₀分别为笛卡尔坐标系中的初始坐标。

(4)利用风场插值法，对两个水平风分量u、v和一个垂直风分量w坐标转换，即将等压层气象数据转换为随地笛卡尔坐标。

(5)计算轨迹所在的位置。粒子所在位置由该时刻风场与计算时间步长决定，即

X(2)＝X(1)*Δt*U

Y(2)＝Y(1)*Δt*V

Z(2)＝Z(1)*Δt*W

其中，X(1)、Y(1)、Z(1)分别为当前时刻放射性核素位置；X(2)、Y(2)、Z(2)分别为下一时刻该放射性核素的位置；U、V、W风别为粒子所在位置风速三维方向的风分量，Δt为时间步长。

(6)根据轨迹所在的位置，确定轨迹是否终止，如果终止，输出轨迹，结束程序；如果未终止或有新轨迹，返回步骤(3)，继续计算。

(三)气载放射性核素长距离迁移轨迹拉格朗日粒子扩散计算

如图6所示，是本发明提供的气载放射性核素长距离迁移轨迹拉格朗日粒子扩散计算方法，该方法包括如下步骤：

1)获取天气预报数据、地形数据、下垫面特征数据、模式参数；

天气预报数据包括：(a)空间三维场数据：两个水平风分量u和v，一个垂直风分量w，温度T，特征湿度；(b)空间二维场数据：地表压力，雪厚度，海面压力，云状，10m高u、v风速，2m高温度，2m高露点温度，大尺度降水，地面感应热通量，太阳辐射，地面应力。

模式参数包括：(a)释放源项参数，即***当量、核素种类及放射性活度；(b)模拟区域范围，并对该区域范围进行网格划分；(c)坐标系类型、时间长度、起止时间。

2)空间坐标和时间坐标变换，是通过已知的方法将等压层气象数据转换为随地笛卡尔坐标。

3)调用风场；按照文件给定顺序读取风场文件名称、各风场时刻、以及各时刻各节点气象数据信息。

4)通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；利用风场插值，获得精确的风场数据。

拉格朗日粒子弥散模式是把每个污染质点当成有标志的质点，通过释放大量粒子，计算粒子的轨迹，而这些粒子描述了气载污染物在大气中的迁移扩散。粒子在流场中按平均风输送，同时又用一系列随机位移来模拟湍流扩散，这样就表达了平流和湍流扩散两种作用，最后由这些粒子在空间和时间上的总体分布估算出污染物的分布。

积分粒子运动的拉格朗日方程，将粒子的运行轨迹写成下列形式：

x_i(t+Δt)＝x_i(t)+v_i(x_i(t)，t)Δt+v′_i(x_i(t)，t)Δt i＝1，2，3. (16)

其中，x_i为粒子的三维坐标分量；v_i为平均风速分量v′_i为湍流脉动速度分量(u′,v′,w′)；t为时间序列；Δt为时间步长。每个时间步长的脉动速度通过假设运动遵从Markov假定，即

其中，式中右边第二项代表速度涨落量中的随机部分，ξ为符合高斯分布(平均值为0、标准差为σ_i)的随机数；和为脉动量v′_i的标准差；R(Δt)＝exp(-Δt/τ)为拉格朗日自相关函数；τ为Lagrangian时间尺度；w分量公式中右边的第三项是为避免质点在低能区的堆积而引入的修正项。

方程(16)考虑了平均风输送影响和大气湍流风脉动的影响，湍流脉动反映了时间尺度小于1小时，对应于较短的长度尺度。中尺度运动可以使弥散的烟羽显著增大(Gupta等，1997)，对于大尺度模拟问题，需要考虑中尺度风脉动影响。因此，考虑中尺度风脉动影响的粒子运动方程为：

x_i(t+Δt)＝x_i(t)+v_i(x_i(t)，t)Δt+v′_i(x_i(t)，t)Δt+v″_i(x_i(t)，t)Δt i＝1，2，3 (18)

其中，v″_i为中尺度风脉动速度分量(u″,v″,w″)。

本发明运用拉格朗日粒子弥散模拟方法的关键在于确定决定粒子运动平均轨迹的平均风场、三个速度分量的拉格朗日时间尺度和涨落的标准差。

对平均风场，通过解二阶差分形式的轨迹运动方程，确定空气质点的运行轨迹，式中，x_i(t)＝[x(t),y(t),z(t)]为t时刻质点的坐标(x，y，z)；为质点在位置处x_i(t)的风速。当计算步长Δt内气团运动加速度恒定，则运动轨迹的坐标由下列方程表述，当给定初始条件v₀[x_0i(t),t],即u(x₀,y₀,z₀,t)、v(x₀,y₀,z₀,t)和w(x₀,y₀,z₀,t),通过迭代计算确定空气气团运动的轨迹，其中，x₀、y₀、z₀分别为笛卡尔坐标系中的初始坐标。

对于湍流参数，采用Hanna(1982)提出的一种参数化方法，即根据边界层参数混合层高度h,莫宁-奥布霍夫长度L,对流速度尺度w_*,粗糙长度z₀和摩擦速度u_*来计算湍流参数。由于Hanna的方法在整个行星边界层并不总能得到光滑的σ_w廓线，导致粒子不能充分混合，因此采用Ryall和Maryon(1997)提出的一种修正方法确定σ_w。

对于边界层参数(如L，u_*)，利用模式第一层高度、地面10m和2m高度的风温数据然后用廓线方法计算上述参数，采用迭代的方法解下列方程：

其中，κ，Karman常数；z_i，模式第一层高度；Δu，模式第一层高度和10m高之间的风速差；ΔΘ，模式第一层高度和2m高之间的位温差；g，重力加速度；Θ_*，温度尺度；L，莫宁-奥布霍夫长度；平均地面气温；Ψ_m和Ψ_h，动量和热量的稳定度修正函数，其函数形式为

其中，φ₁和φ₂分别为风速和温度的廓线函数，其形式分别为

对于中尺度风脉动参数，采用Maryon(1998)提出的方法：通过假定中尺度风速脉动与Hanna参数化方法覆盖的湍流脉动无关，解一个独立的拉格朗日方程来解决，所用到的时间尺度和速度的方差是通过对一个测站的风观测时间序列进行谱分析得到的。假定在网格尺度上观测的风的方差也为亚网格风速标准差提供一些信息，这样，用一种简化方法确定风速标准差的取值方法是：计算粒子所在位置周围16个网格点(在时间和空间上)的速度的标准差，然后取该标准差的一半作为解中尺度拉格朗日方程时所用的风速标准差。

5)采用源耗减的概念来计算干、湿沉积

给定物质的干沉积用沉积速度V_d(m/s)来描述，可根据核素类型和下垫面特征来确定。雨水冲洗效应是造成地面高放射性污染水平的最重要的因素之一。湿沉积可以用类似干沉积的方法计算，不同之处仅在于用冲洗系数Λ(s^-1)代替干沉积速度，冲洗系数的大小取决于降雨强度。放射性衰变的估算：根据公式

N(1)＝N(0)*0.5^Δt*λ

其中，N(1)为衰减后的放射性物质量；N(0)为衰减前的放射性物质量；△t为衰变时间，s；λ为衰变常数，s^-1。

干、湿沉积量的估算。只针对气溶胶进行沉积计算。

C_d＝C*V_d

C_d为沉降后的浓度，Bq/m²；

C为空气浓度，Bq/m³；

V_d为沉降速度，m^-1。

当近地层有降水发生时，同时计算干沉积与湿沉积量；

当近地层无降水发生时，只计算干沉积。

6)统计粒子密度、计算网格浓度

为了提高扩散计算精度和效率，当采用粒子***技术和核函数方法。利用核函数方法对空间关注点浓度进行计算。核函数方法认为离散化后的烟团在迁移过程中，其自身也呈现高斯扩散的趋势。对本模式而言，由于空间网格间距较大，而单纯的粒子模式如果不使用网格嵌套技术，则网格内部呈现不合理的均匀浓度，并引起监测点浓度的***误差。所以，在模式中引入核函数方法，这样不仅可以避免复杂的嵌套计算，而且可以快速合理地得到有限关注点的浓度。同时，对于一般意义的核函数方法，不强调单个质子团的空间扩散方向性，但由于大尺度模式中，质子团稀疏，对称分布的核函数仍然会引起不真实的空间浓度分布。对此，模式根据气象条件设定不同参数，使质子团的核函数扩散呈现更为真实的三维不对称分布。

放射性核素空气浓度的计算方法：每个网格中的浓度正比于质点通过该网格时所需时间的总和，因此每个网格的浓度C_k(Bq/m³)用下列公式计算：

其中，Q，释放的放射性活度，Bq；N，释放的粒子总数；T，所有粒子扩散时间的总和，s；N_ik，第i个粒子在k网格中穿行时间内的时间步数；Δt，时间步长，s；V_k，网格k的体积，m³。在浓度计算中，根据放射性核素的半衰期和粒子迁移的总时间对放射性衰变进行修正。

7)输出时间间隔

如果达到设定的浓度输出时间，计算该时刻网格的浓度，并输出空气浓度场；如果没有达到设定的浓度输出时间，考虑粒子***条件，跟踪所有粒子、检查是否有新粒子释放；如果有新粒子释放返回步骤3)，如果没有新粒子释放，结束程序。

其中：粒子***条件，是当某粒子所在空间的12个相邻空间均无粒子分布时，对粒子进行***。即，***后的粒子质量减半，总模拟粒子数增加1。跟踪所有的粒子，检查是否有新粒子释放，是判断模拟时间是否处于释放时间段内，如果处于释放段内，则新增N个粒子的初始迁移位置，并与其它已有粒子一同进行迁移扩散模拟。

对于方程(3)求解，当给定初始条件，通过迭代计算就可以确定各随机游走粒子的运动轨迹。

(四)气载放射性物质释放剂量估算

核与辐射突发事件后果大尺度气载放射性物质释放剂量估算方法，包括空气浸没外照剂量、地表沉积外照剂量、吸入内照剂量、个人有效剂量以及甲状腺剂量五类进行合计，给出该事故下气载放射性物质在不同时刻可能引起的环境潜在剂量水平。

1、烟云外照射剂量估算

对于烟云长距离迁移，采用“半无限烟云”方法。对于一给定的时段，远距离计算网格中的烟云外照射剂量通过下列公式计算：

式中：t₁、t₂代表某时段的起始和结束时刻；org代表某种器官；DC为在网格(x,y)处，器官org在t₁至t₂时段内所受到的烟云外照射剂量(Sv)；AC为在网格(x,y)处，在t₁至t₂时段内核素nuc的积分空气浓度(Bq·s/m³)；DFC为核素nuc对器官org的半无限烟云剂量转换因子[Sv·m³/(Bq·s)]；SFC为烟云外照射的屏蔽因子(无量纲)，由t₁至t₂时段公众所采取的防护行动的种类决定。

这里的AC为在网格(x,y)处，在t₁至t₂时段内核素nuc的积分空气浓度Bq·s/m³)，该核素积分空气浓度根据计算得出，计算方法如下：首先计算每个网格的浓度，每个网格中的浓度正比于质点通过该网格时所需时间的总和，因此每个网格的浓度C_k(Bq/m³)用下列公式计算：

其中，Q，释放的放射性活度，Bq；N，释放的粒子总数；T，所有粒子扩散时间的总和，s；N_ik，第i个粒子在k网格中穿行时间内的时间步数；Δt，时间步长，s；V_k，网格k的体积，m³；然后根据网格的浓度，进行累加，获得核素积分空气浓度。

Q，释放的放射性活度，该活度根据核素的衰变常数、总裂变次数以及核素裂变产额三者之间的乘积获得；核素裂变产额需要查阅文献。

N，释放的粒子总数，采用粒子随机游走的方法模拟大气扩散，该模式中把每个污染的质点当成有标志的质点，通过释放大量粒子，通过已知方法计算粒子的轨迹，而这些粒子描述了气载污染物在大气中的迁移扩散，最后由这些粒子在空间和时间上的总体分布估算出释放粒子总数；

Δt，时间步长，该时间步长的确定取决于网格的格距、风速和离散风场的时间间隔，从空间角度遵循：

其中，Δx_i为空间格距，i代表三个方向；v_i分别代表三个方向的风分量。从时间角度遵循：

其中，Δt_风场为输入风场的时间间隔，即风观测或数值天气预报风场的周期。粒子迁移计算模式的时间步长最大值要小于Δt_空间和Δt_时间的最小值，最小时间步长可为1秒。

DFC为核素nuc对器官org的半无限烟云剂量转换因子[Sv·m³/(Bq·s)]，该转换因子可在现有的文献中获得；

SFC为烟云外照射的屏蔽因子，由t₁至t₂时段公众所采取的防护行动的种类决定，如果受照射人员采取不撤离环境的方式，屏蔽因子为1；如果采取撤离环境的方式，屏蔽因子为0；对于进入房屋等障碍物内时，用户需要给出实际的屏蔽因子；屏蔽因子值应为大于0的小数。

2、地面沉积外照射剂量的估算

对于一给定的时段，地面沉积外照射剂量通过下面的公式计算：

式中，DG为在网格(x,y)处，器官org在t₁至t₂时段内所受到的地面沉积外照射剂量(Sv)；GC为在网格(x,y)处，t时刻核素nuc的地面沉积浓度(Bq/m²)。随时间变化的地面沉积浓度与瞬时空气浓度一样，由扩散模式给出。GRF为核素nuc对器官org的地面沉积外照射剂量转换因子[Sv·m²/(Bq·s)]，该剂量转换因子可在现有的文献中获得；SFG为地面沉积外照射的屏蔽因子(无量纲)，由t₁至t₂时段公众所采取的防护行动的种类决定。

3、烟云吸入内照射剂量的估算

对于给定未服用稳定碘的时段t₁至t₂内的烟云吸入内照射剂量，计算公式：

式中，DI为在网格(x,y)处，器官org在t₁至t₂时段内吸入放射性物质所造成的待积吸入剂量(Sv)；AC为在网格(x,y)处，在t₁至t₂时段内核素nuc的积分空气浓度(Bq·s/m³)；IF为核素nuc对器官org的待积吸入内照射剂量转换因子(Sv/Bq)；BR为公众在t₁至t₂时段内的呼吸率(m³/s)；IPF为公众在t₁至t₂时段内的吸入防护因子(无量纲)。

4、总剂量水平估算

考虑到核爆源项核素种类繁多，为避免对所有释放核素逐一扩散计算带来的时间费用，按照其半衰期长短(30天为界)、是否引起甲状腺剂量、是否存在地面沉积剂量等因素将所有关注核素分为五类，依次为碘组、有沉积短半衰期、有沉积长半衰期、无沉积短半衰期、无沉积长半衰期。同时，在各组中选择某核素作为各类代表性核素，对应上述五类，其代表性核素依次为I-131、Ba-140、Cs-137、Xe-133、Kr-85。

每次扩散只针对代表性核素进行模拟，而剂量估算时再根据实际释放核素量以及核素半衰期差异对其浓度进行校正，进而对各核素引起的剂量水平进行估算，最后对所有核素引起的剂量进行累计，从而得到该次释放总剂量分布。

以下是本发明在计算过程中所涉及到的数据：

(1)坐标系

模式采用混合坐标系，即x、y、η坐标，其中η是一种对气压坐标进行变换的垂直坐标，η与气压坐标的变换方法为：

p_k＝A_k+B_kp_s

η_k＝A_k/p₀+B_k

其中，η_k为模式第k层的η值；p_s为地表压力；p₀为压力常数(101325Pa)。A_k和B_k为系数，由最接近地表(随地坐标)的值和最大压力高度层的值确定，中间高度层的系数值则根据随地表层和压力高度层之间的压力梯度确定。

此外，如果纬度超过75°时考虑极地立体投影。

(2)时间、空间分辨率

模式的时间、空间分辨率取决于两方面因素：一是输入资料的时间、空间分辨率；二是模式计算的时间、空间分辨率，在输入资料的时空分辨率一定的情况下，模式计算的时空分辨率主要取决于时间步长。

①输入资料的时间、空间分辨率

在数值模式中，地形和气象数据的空间分辨率是一致的，可接受的空间分辨率从0.25°到2.5°，时间分辨率从3h到12h。

数值计算分析结果表明，随着风场网格分辨率的降低，数值模拟结果的准确程度也随之降低。数值模拟实验表明：风特征大约可以被0.5°(约45km)空间分辨、6小时时间分辨水平的预报风场解析，而利用2.5°(约225km)与12小时的时空分辨预报风场则无法满足需要。因此，可以认为6小时的风场时间间隔是大尺度模式运行的基本条件。综合考虑模拟精度和计算资源与计算速度，在本***中采用0.5°的空间分辨率、3～6小时的时间分辨率水平完全满足实际应用的需要。

②时间步长

时间步长的确定取决于网格的格距、风速和离散风场的时间间隔，从空间角度遵循：

其中，Δx_i为空间格距，i代表三个方向；v_i分别代表三个方向的风分量。从时间角度遵循：

其中，Δt_风场为输入风场的时间间隔，即风观测或数值天气预报风场的周期。粒子迁移计算模式的时间步长最大值要小于Δt_空间和Δt_时间的最小值，最小时间步长可为1秒。

本发明还提供一种核与辐射突发事件后果大尺度模拟***，包括模拟参数准备模块、源项估算模块、轨迹计算模块、粒子扩散计算模块及剂量计算模块；

模拟参数准备模块，用于对天气预报数据及轨迹方向、时间长度、起止时间、轨迹时间间隔、轨迹初始坐标、模拟区域范围的准备；

源项估算模块，用于对已知的***当量对***烟云的高度和半径的估算、根据核***后随时间迁移产生的剂量份额不同对核素种类的估算、根据核素的衰变常数、裂变次数及裂变产额对放射性核素活度的估算；

轨迹计算模块，通过解二阶差分形式的轨迹运动方程，确定空气质点的运行轨迹；

粒子扩散计算模块，通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；计算放射性衰变和干、湿沉积；统计粒子密度、计算网格浓度。

剂量计算模块，用于气载放射性物质剂量估算方法，包括对烟云外照射剂量估算、对地面沉积外照射剂量估算、对烟云吸入内照射剂量估算、对总剂量水平估算。

本发明模拟***中的四个物理模块，可交互运行，也可根据需要进行选择。本发明模拟***具有：

(1)可靠性

运行及数据(库)管理以任务ID的方式管理，使得***运行管理、用户操作有一套清晰的思路，避免了交叉，***运行更加安全、可靠。

(2)开放性和兼容性

整个的运行设计是开放的，用户可以根据需要方便地选择物理模块的运行方式。***输入、输出参数格式公开。

(3)可扩展性

用户根据需要，可以编辑源项，根据需要选择计算核素；可以更改地理信息***的图例类型；可以修改、增加、删除地理信息***的图层。

(4)可维护性

由于***设计安全可靠，***本身一般不需要太多的维护。

(5)可管理性

***采用通用的Windows操作***，使用标准的Windows界面，对***的管理简单。

本发明模拟***输入输出及运行过程：

对于本***而言，需要用户确定的输入信息可分为模拟背景信息、源项信息、模拟时间信息三类。其中，模拟背景信息包括：评价类型、运算类型、气象数据的选取、模拟区域、扩散模拟关注点名称、位置；源项信息包括：核爆类型与当量、释放水平位置、释放高度、释放总量、释放核素类型以及非核爆、非均匀核素释放总量；模拟时间信息包括：风场启始时间、释放启止时间、模拟启止时间信息。

根据操作及***运行需要，各信息的输入方式一般分三类。具体有：地图点选、手工录入、***默认。

***输出包括主要气象数据、扩散结果、轨迹结果和剂量结果。其中，气象数据包括三维风场、三维温度场、气压、降雨；扩散模拟结果包括逐时瞬时空气浓度、地面沉积浓度、空气时间积分浓度、关注点瞬时空气浓度与时间积分浓度；轨迹模拟结果包括逐时粒子位置；剂量模拟结果包括逐时浸没外照、吸入内照、地面沉积外照、个人有效剂量，以及甲状腺剂量。***输出包括图形、表格与数据文件形式。

本发明的核与辐射突发事件后果大尺度模拟***及模拟方法，并不限于上述具体实施方式，本领域技术人员根据本发明的技术方案得出其他的实施方式，同样属于本发明的技术创新范围。

Claims (10)

1.一种核与辐射突发事件后果大尺度模拟方法包括如下步骤：

1)对模拟参数的准备；即对天气预报数据及模拟方向、模拟起止时间、释放起止时间、轨迹时间间隔、释放位置、模拟区域范围的准备；

2)对核与辐射事件后果进行模拟计算，包括：

对核爆源项估算；即通过TNT当量对***烟云的高度和半径的估算；根据裂变方程，对核***后裂变产物、放射性活度、份额随时间变化的估算；

对气载放射性核素长距离迁移轨迹的计算；即通过解二阶差分形式的轨迹运动方程，确定空气质点的运行轨迹；

对气载放射性核素长距离迁移拉格朗日粒子扩散的计算；即通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；通过统计粒子数，估算网格浓度；

对气载放射性物质剂量的估算；即对烟云外照射剂量估算；对地面沉积外照射剂量估算；对烟云吸入内照射剂量估算；对总剂量水平估算。

2.如权利要求1所述的核与辐射突发事件后果大尺度模拟方法，其特征是：步骤2)中，所述***烟云高度和半径通过下列公式进行估算：

其中：W为***当量；LW＝log₁₀(W/kT)。

3.如权利要求1所述的核与辐射突发事件后果大尺度模拟方法，其特征是：步骤2)中，所述裂变方程如下：

式中：i-单线链中核素的序号；λ_i为第i个核素的衰变常数；y_i-单线链中i核素的绝对独立产额；p(t)-单位时间内核燃料的裂变次数；A_i-同量异位链中i核素的放射性。

4.如权利要求1所述的核与辐射突发事件后果大尺度模拟方法，其特征是：步骤2)中，所述二阶差分形式的轨迹运动方程如下： 式中，x_i(t)＝[x(t),y(t),z(t)]为t时刻质点的坐标(x,y,z)；为质点在位置处x_i(t)的风速； 分别为三维平均风分量。

5.如权利要求4所述的核与辐射突发事件后果大尺度模拟方法，其特征是：如果计算时间步长Δt内气团运动加速度恒定，那么运动轨迹的坐标可由下列方程表述：当给定初始条件v₀[x_0i(t),t],即u(x₀,y₀,z₀,t)、v(x₀,y₀,z₀,t)和w(x₀,y₀,z₀,t),通过迭代计算确定空气气团运动的轨迹；x₀,y₀,z₀分别为笛卡尔坐标系中的初始坐标。

6.如权利要求1所述的核与辐射突发事件后果大尺度模拟方法，其特征是：步骤2)中，所述积分粒子运动的拉格朗日方程为：

x_i(t+Δt)＝x_i(t)+v_i(x_i(t),t)Δt+v′_i(x_i(t),t)Δt+v″_i(x_i(t),t)Δt i＝1,2,3

其中，x_i为粒子的三维坐标分量；v_i为平均风速分量v′_i为湍流脉动速度分量(u′,v′,w′)；v″_i为中尺度风脉动速度分量(u″,v″,w″)；t为时间序列；Δt为时间步长。

7.如权利要求1所述的核与辐射突发事件后果大尺度模拟方法，其特征是：步骤2)中，所述网格浓度通过下列公式计算：

其中，Q，释放的放射性活度；N，释放的粒子总数；T，所有粒子扩散时间的总和；N_ik，第i个粒子在k网格中穿行时间内的时间步数；Δt，时间步长；V_k，网格k的体积。

8.如权利要求1所述的核与辐射突发事件后果大尺度模拟方法，其特征是：步骤2)中，所述烟云外照射剂量通过下列公式计算：

其中：t₁、t₂代表某时段的起始和结束时刻；org代表某种器官；DC为在网格(x,y)处，器官org在t₁至t₂时段内所受到的烟云外照射剂量；AC为在网格(x,y)处，在t₁至t₂时段内核素nuc的积分空气浓度；DFC为核素nuc对器官org的半无限烟云剂量转换因子；SFC为烟云外照射的屏蔽因子，由t₁至t₂时段公众所采取的防护行动的种类决定。

所述地面沉积外照射剂量通过下列公式计算：

其中：DG为在网格(x,y)处，器官org在t₁至t₂时段内所受到的地面沉积外照射剂量；GC为在网格(x,y)处，t时刻核素nuc的地面沉积浓度；GRF 为核素nuc对器官org的地面沉积外照射剂量转换因子；SFG为地面沉积外照射的屏蔽因子，由t₁至t₂时段公众所采取的防护行动的种类决定。

所述烟云吸入内照射剂量通过下列公式计算：

其中：DI为在网格(x,y)处，器官org在t₁至t₂时段内吸入放射性物质所造成的待积吸入剂量；AC为在网格(x,y)处，在t₁至t₂时段内核素nuc的积分空气浓度；IF为核素nuc对器官org的待积吸入内照射剂量转换因子；BR为公众在t₁至t₂时段内的呼吸率；IPF为公众在t₁至t₂时段内的吸入防护因子。

9.如权利要求1所述的核与辐射突发事件后果大尺度模拟方法，其特征是：步骤2)中，所述总剂量水平的估算，是将烟云外照射剂量估算的结果、地面沉积外照射剂量估算的结果以及烟云吸入内照射剂量估算的结果进行累计。

10.一种核与辐射突发事件后果大尺度模拟***，其特征是：包括模拟参数准备模块、源项估算模块、轨迹计算模块、粒子扩散计算模块、剂量计算模块；

模拟参数准备模块：用于对天气预报数据及轨迹方向、时间长度、起止时间、轨迹时间间隔、轨迹初始坐标、模拟区域范围的准备；

源项估算模块：用于对核爆源项的估算，包括对***烟云的高度和半径的估算；根据裂变方程，对核***后裂变产物、放射性活度、份额随时间变化的估算；

源项估算模块：用于对气载放射性核素长距离迁移轨迹的计算，即通过 解二阶差分形式的轨迹运动方程，确定空气质点的运行轨迹；

粒子扩散计算模块：用于对气载放射性核素长距离迁移拉格朗日粒子扩散的计算，即通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；统计粒子数、计算网格浓度。

剂量计算模块：用于对气载放射性物质剂量的估算，包括对烟云外照射剂量估算、对地面沉积外照射剂量估算、对烟云吸入内照射剂量估算、对总剂量水平估算。