CN107133455B - 利用耦合蒙特卡罗方法模拟ads系统瞬态问题的方法 - Google Patents

Abstract

一种利用耦合蒙特卡罗方法模拟ADS系统瞬态问题的方法，包括如下步骤：1、利用预估校正准静态方法得到预估通量方程；2、使用转置矩阵计算0时刻共轭通量；3、引入裂变源系数，进行迭代加速求解预估通量方程；4、将前一步得到的通量分布用共轭通量归一化得到目标时刻的通量；采用转置矩阵和构造迭代加速的操作，解决共轭通量计算和低次临界度下模拟不稳定的问题；与全确定论的瞬态计算方法相比，具有计算精度高的优点，同时解决了在耦合蒙卡方法中的问题；该方法将动力学部分与中子输运部分进行外耦合，方便使用各种蒙卡程序进行替换，具有通用性。

Description

利用耦合蒙特卡罗方法模拟ADS系统瞬态问题的方法

技术领域：

本发明涉及核反应堆设计和反应堆物理计算领域，具体涉及一种利用耦合蒙特卡罗方法模拟ADS系统瞬态问题的方法。

背景技术：

加速器驱动次临界反应堆系统(ADS)是一种新型核能反应堆系统，通过加速器产生高能质子流，轰击散裂靶产生散裂中子，作为外中子源驱动次临界堆芯，嬗变其中的核废料。该系统有中子通量密度高和能谱较硬的特点，其U/Pu核素装载份额的减少会导致堆芯的有效缓发中子份额减少，也会导致燃料的多普勒温度反馈效应急剧减弱，降低了超功率过程中的负反馈机制。该型反应堆的瞬态特性特殊，需要对其进行精确的计算以保证堆芯设计的安全。

现有的反应堆瞬态分析方法中，中子动力学的计算为确定论方法，将时空中子动力学方程经全隐式向后差分后离散为求解各时间步上的固定源形式的中子学方程，得到每个时间步上的通量。现有方法中，多采用求解扩散方程的形式得到通量，精确度上还有待提高。

蒙特卡罗方法利用模拟粒子在堆内的行为进行计算，最后统计得到通量参数。由于其属于统计学方法，根据数学上的大数定律和中心极限定律可确定，当样本数量足够大时，能够收敛于真解。将蒙特卡罗方法用于求解中子输运方程可得到更为精确的通量值，用于动力学可得到更为精确的动力学计算结果。

针对加速器驱动次临界问题，在实际将蒙特卡罗用于动力学各时间步下的中子学求解时，由于需要在计算初始给出共轭通量作为归一化预估通量的权重，需要对蒙特卡罗计算进行处理，另外由于问题处于低次临界度下，模拟时中子有接近于1的增殖比，难以截断，造成模拟时间长且计算不稳定，容易发散。

从以上分析可看出，要针对加速器驱动次临界系统进行精确的瞬态分析计算，需要利用耦合蒙特卡罗方法进行分析，而耦合蒙特卡罗方法需要解决共轭通量计算和低次临界度下模拟不稳定的问题。

发明内容：

为了克服上述现有技术存在的问题，本发明的目的在于提供一种利用耦合蒙特卡罗方法模拟ADS系统瞬态问题的方法，采用转置矩阵和构造迭代加速的操作，解决共轭通量计算和低次临界度下模拟不稳定的问题；与全确定论的瞬态计算方法相比，具有计算精度高的优点，同时解决了在耦合蒙卡方法中的问题；该方法将动力学部分与中子输运部分进行外耦合，方便使用各种蒙卡程序进行替换，具有通用性。

为达到上述目的，本发明采取了以下技术方案予以实施：

一种利用耦合蒙特卡罗方法模拟ADS系统瞬态问题的方法，步骤如下：

步骤1：采用预估校正改进准静态方法处理带缓发中子先驱核的中子时空动力学方程，得到完全的预估通量方程，如公式(1)：

式中：

n--时间步序号；

Ω--空间角度；

r--空间位置；

E--中子能量；

--n+1时刻的预估通量；

Σ'_t,n+1--n+1时刻的形式总截面，包括n+1时刻的总截面和时间相关总截面；

Σ_s,n+1(r,E',Ω'→E,Ω)--n+1时刻，位置r处，中子从E'能量Ω'角度散射到E能量Ω角度的散射截面；

F'_t,n+1--n+1时刻的形式裂变产生截面，包括n+1时刻的瞬发中子产生截面和n+1时刻的时间相关的缓发中子产生截面；

S_n--n时刻的形式外源，包括实际外源，缓发中子先驱核相关源和时间相关源；

同时在预估校正改进准静态方法中，需要对预估通量进行归一化，其权重因子为共轭通量；

归一化公式为：

式中：

v--中子速度；

ψ--形状函数分布；

--0时刻的共轭通量；

步骤2：根据步骤1公式推导要求，按照公式(4)计算0时刻的共轭通量

式中：

ψ^*--共轭通量；

Σ--总截面；

Σ_s(r,E,Ω→E',Ω')--在位置r处，中子从E能量Ω角度散射到E'能量Ω'角度的散射截面；

νΣ_f(r,E)--在位置r处，能量为E的中子的裂变产生截面；

χ(E')--裂变产生中子能量为E'的占比；

ψ^*(r,E',Ω')--在位置r处，能量为E'，角度为Ω'的共轭通量；

通过转置裂变矩阵和散射矩阵的方式，保持原边界条件不变，进行前向计算，所得标通量即为对应的共轭通量，其分布能够用于对动力学计算中预估通量的计算；

步骤3：根据步骤1推导得出的公式(1)，对公式(1)进行蒙特卡罗输运计算，其中外源由n时刻的输运计算得到，n＝0时即为初始时刻实际外源；截面以蒙特卡罗指定材料的形式给出；目标为计算n+1时刻的预估通量在计算过程中，进行以下迭代格式以加速和增加稳定性：

式中：

m--迭代次数序号；

k_m--迭代第m步时的裂变源系数；

k_m-1--迭代第m-1步时的裂变源系数；

S₀--归一化后的外源；

Φ_m--迭代第m步时通量

Φ_m-1--迭代第m-1步时的通量

L--泄漏、吸收、散射项算子；

B--裂变产生项算子；

<>--关于全相空间的积分算子；

在蒙特卡罗中子模拟中，将裂变反应视为吸收，在第m-1次计算时模拟计算相应通量，并通过式(6)计算迭代系数k_m-1，然后通过式(5)构造新的源项，不改变模拟问题的材料组成，进行第m次计算，反复迭代，直到两次迭代中的k_m-1和k_m的差小于10^-6时，视为收敛，此时的Φ_m即为n-1时刻的预估通量

步骤4：根据步骤3求得的n+1时刻的预估通量和步骤2求得的0时刻的共轭通量由公式(7)，经归一化得到n+1时刻的通量：

式中：

ψ_n+1--第n+1时刻的形状函数；

由预估校正改进准静态方法，在形状函数已知的情况下，能够求出n+1时刻的点堆参数，再由点堆动力学方程能够得到n+1时刻的通量幅度T_n+1；最终可得到n+1时刻的通量：

φ_n+1＝T_n+1ψ_n+1 公式(8)

重复步骤3步骤4继续计算n+2时刻的通量，依次类推，直到计算到目标时刻。进而得到通量随时间和空间的分布φ(r,E,t)，完成反应堆瞬态计算。

本发明与现有技术相比具有如下优点：

1、采用耦合蒙特卡罗方法，且为外部耦合，可以提高瞬态计算中各时间步上中子输运计算的准确度，另外还能灵活选取蒙特卡罗求解程序。

2、将裂变源迭代加速技巧用于蒙特卡罗求解，解决了低次临界度下固定源问题模拟缓慢且不稳定的问题，使蒙特卡罗方法更好地适用于加速器驱动次临界系统的瞬态分析计算。

附图说明

图1为利用耦合蒙特卡罗方法模拟ADS系统瞬态问题的方法的计算流程图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明：

本发明在原有的确定论的预估校正准静态计算方法中，将瞬态分析中分时间步的中子输运计算用蒙特卡罗方法进行耦合，通过转置裂变截面和散射截面处理共轭通量计算，通过引入裂变源迭代因子，加速并稳定蒙特卡罗中子模拟，达到预期效果。

如图1所示，本发明的具体实施步骤如下：

步骤1：采用预估校正改进准静态方法处理带缓发中子先驱核的中子时空动力学方程，得到完全的预估通量方程：

式中：

n--时间步序号；

Ω--空间角度；

r--空间位置；

E--中子能量；

--n+1时刻的预估通量；

Σ'_t,n+1--n+1时刻的形式总截面，包括n+1时刻的总截面和时间相关总截面；

Σ_s,n+1(r,E',Ω'→E,Ω)--n+1时刻，位置r处，中子从E'能量Ω'角度散射到E能量Ω角度的散射截面；

F'_t,n+1--n+1时刻的形式裂变产生截面，包括n+1时刻的瞬发中子产生截面和时间相关的缓发中子产生截面；

S_n--n时刻的形式外源，包括实际外源，缓发中子先驱核相关源和时间相关源；

同时在预估校正改进准静态方法中，需要对预估通量进行归一化，其权重因子为共轭通量；

归一化公式为：

式中：

v--中子速度；

ψ--形状函数分布；

--共轭通量分布；

步骤2：根据第1步公式推导要求，计算共轭通量

中子输运稳态前向方程表达式如下：

中子输运稳态共轭通量表达式如下：

式中：

ψ^*--共轭通量；

Σ--总截面；

Σ_s(r,E,Ω→E',Ω')--在位置r处，中子从E能量Ω角度散射到E'能量Ω'角度的散射截面；

νΣ_f(r,E)--在位置r处，能量为E的中子的裂变产生截面；

χ(E')--裂变产生中子能量为E'的占比；

ψ^*(r,E',Ω')--在位置r处，能量为E'，角度为Ω'的共轭通量；

由公式(3)和公式(4)可知，在0时刻，通过转置散射截面，和改变裂变中子产生截面与裂变中子能谱的积分顺序，再通过前向计算求解方法求得所得标通量，即为对应的共轭通量其分布可用于对动力学计算中预估通量的计算；

步骤3：

根据步骤1推导得出的公式，对式公式(1)进行蒙特卡罗输运计算，其中外源S_n由n时刻的输运计算的通量φ_n得到，n＝0时即为初始时刻实际外源；截面以蒙特卡罗指定材料的形式给出；通过目标为计算n+1时刻的预估通量在计算过程中，进行以下迭代格式以实现加速和增加稳定性：

式中：

m--迭代次数序号；

k_m--迭代第m步时的裂变源系数；

k_m-1--迭代第m-1步时的裂变源系数；

S₀--归一化后的外源；

Φ_m--迭代第m步时通量

Φ_m-1--迭代第m-1步时的通量

L--泄漏、吸收、散射项算子；

B--裂变产生项算子；

<>--关于全相空间的积分算子；

在蒙特卡罗中子模拟中，将裂变反应视为吸收，在第m-1次计算时模拟计算相应通量Φ_m-1，并通过式(6)计算迭代系数k_m-1，然后通过式(5)构造新的源项，不改变模拟问题的材料组成，将源项改为以上计算的源进行第m次计算，反复迭代，直到两次迭代中的k_m-1和k_m的相对误差小于10-⁶时，视为收敛，此时的即为n+1时刻的预估通量由于在模拟中将裂变视为吸收，即使处在较低次临界度下，也不会出现中子增殖比接近1而无法截断的问题，大大减少了模拟时间，同时保证模拟能够完成，增加了模拟的稳定性；

步骤4：

根据步骤3求得的n+1时刻的预估通量和步骤2求得的0时刻的共轭通量由式(7)，归一化得到n+1时刻的通量：

ψ_n+1--第n+1时刻的形状函数；

由预估校正改进准静态方法，在形状函数已知的情况下，可以求出n+1时刻的点堆参数，由此可以得到n+1时刻的通量幅度T_n+1；最终可得到n+1时刻的通量：

φ_n+1＝T_n+1ψ_n+1 公式(8)

重复步骤3可继续计算n+2时刻的通量，依次类推，直到计算到目标时刻；进而可得到通量随时间和空间的分布φ(r,E,t)，完成反应堆瞬态计算。

Claims (1)

1.一种利用耦合蒙特卡罗方法模拟ADS系统瞬态问题的方法，其特征在于：步骤如下：

步骤1：采用预估校正改进准静态方法处理带缓发中子先驱核的中子时空动力学方程，得到完全的预估通量方程，如公式(1)：

式中：

n--时间步序号；

Ω--空间角度；

r--空间位置；

E--中子能量；

--n+1时刻的预估通量；

∑'_t,n+1--n+1时刻的形式总截面，包括n+1时刻的总截面和时间相关总截面；

∑_s,n+1(r,E',Ω'→E,Ω)--n+1时刻，位置r处，中子从E'能量Ω'角度散射到E能量Ω角度的散射截面；

F'_t,n+1--n+1时刻的形式裂变产生截面，包括n+1时刻的瞬发中子产生截面和n+1时刻的时间相关的缓发中子产生截面；

S_n--n时刻的形式外源，包括实际外源，缓发中子先驱核相关源和时间相关源；

同时在预估校正改进准静态方法中，需要对预估通量进行归一化，其权重因子为共轭通量；

归一化公式为：

式中：

v--中子速度；

ψ--形状函数分布；

--0时刻的共轭通量；

步骤2：根据步骤1公式推导要求，按照公式(4)计算0时刻的共轭通量

式中：

ψ^*--共轭通量；

∑--总截面；

∑_s(r,E,Ω→E',Ω')--在位置r处，中子从E能量Ω角度散射到E'能量Ω'角度的散射截面；

v∑_f(r,E)--在位置r处，能量为E的中子的裂变产生截面；

χ(E')--裂变产生中子能量为E'的占比；

ψ^*(r,E',Ω')--在位置r处，能量为E'，角度为Ω'的共轭通量；

通过转置裂变矩阵和散射矩阵的方式，保持原边界条件不变，进行前向计算，所得标通量即为对应的共轭通量，其分布能够用于对动力学计算中预估通量的计算；

步骤3：根据步骤1推导得出的公式(1)，对公式(1)进行蒙特卡罗输运计算，其中外源由n时刻的输运计算得到，n＝0时即为初始时刻实际外源；截面以蒙特卡罗指定材料的形式给出；目标为计算n+1时刻的预估通量在计算过程中，进行以下迭代格式以加速和增加稳定性：

式中：

m--迭代次数序号；

k_m--迭代第m步时的裂变源系数；

k_m-1--迭代第m-1步时的裂变源系数；

S₀--归一化后的外源；

Φ_m--迭代第m步时通量

Φ_m-1--迭代第m-1步时的通量

L--泄漏、吸收、散射项算子；

B--裂变产生项算子；

< >--关于全相空间的积分算子；

在蒙特卡罗中子模拟中，将裂变反应视为吸收，在第m-1次计算时模拟计算相应通量，并通过式(6)计算迭代系数k_m-1，然后通过式(5)构造新的源项，不改变模拟问题的材料组成，进行第m次计算，反复迭代，直到两次迭代中的k_m-1和k_m的差小于10^-6时，视为收敛，此时的Φ_m即为n-1时刻的预估通量

步骤4：根据步骤3求得的n+1时刻的预估通量和步骤2求得的0时刻的共轭通量由公式(7)，经归一化得到n+1时刻的通量：

式中：

ψ_n+1--第n+1时刻的形状函数；

由预估校正改进准静态方法，在形状函数已知的情况下，能够求出n+1时刻的点堆参数，再由点堆动力学方程能够得到n+1时刻的通量幅度T_n+1；最终可得到n+1时刻的通量：

φ_n+1＝T_n+1ψ_n+1 公式(8)

重复步骤3步骤4继续计算n+2时刻的通量，依次类推，直到计算到目标时刻，进而得到通量随时间和空间的分布φ(r,E,t)，完成反应堆瞬态计算。