技术领域：
本发明涉及核反应堆设计和反应堆物理计算领域，具体涉及一种利用耦合蒙特卡罗方法模拟ads系统瞬态问题的方法。
背景技术：
：
加速器驱动次临界反应堆系统(ads)是一种新型核能反应堆系统，通过加速器产生高能质子流，轰击散裂靶产生散裂中子，作为外中子源驱动次临界堆芯，嬗变其中的核废料。该系统有中子通量密度高和能谱较硬的特点，其u/pu核素装载份额的减少会导致堆芯的有效缓发中子份额减少，也会导致燃料的多普勒温度反馈效应急剧减弱，降低了超功率过程中的负反馈机制。该型反应堆的瞬态特性特殊，需要对其进行精确的计算以保证堆芯设计的安全。
现有的反应堆瞬态分析方法中，中子动力学的计算为确定论方法，将时空中子动力学方程经全隐式向后差分后离散为求解各时间步上的固定源形式的中子学方程，得到每个时间步上的通量。现有方法中，多采用求解扩散方程的形式得到通量，精确度上还有待提高。
蒙特卡罗方法利用模拟粒子在堆内的行为进行计算，最后统计得到通量参数。由于其属于统计学方法，根据数学上的大数定律和中心极限定律可确定，当样本数量足够大时，能够收敛于真解。将蒙特卡罗方法用于求解中子输运方程可得到更为精确的通量值，用于动力学可得到更为精确的动力学计算结果。
针对加速器驱动次临界问题，在实际将蒙特卡罗用于动力学各时间步下的中子学求解时，由于需要在计算初始给出共轭通量作为归一化预估通量的权重，需要对蒙特卡罗计算进行处理，另外由于问题处于低次临界度下，模拟时中子有接近于1的增殖比，难以截断，造成模拟时间长且计算不稳定，容易发散。
从以上分析可看出，要针对加速器驱动次临界系统进行精确的瞬态分析计算，需要利用耦合蒙特卡罗方法进行分析，而耦合蒙特卡罗方法需要解决共轭通量计算和低次临界度下模拟不稳定的问题。
技术实现要素：
：
为了克服上述现有技术存在的问题，本发明的目的在于提供一种利用耦合蒙特卡罗方法模拟ads系统瞬态问题的方法，采用转置矩阵和构造迭代加速的操作，解决共轭通量计算和低次临界度下模拟不稳定的问题；与全确定论的瞬态计算方法相比，具有计算精度高的优点，同时解决了在耦合蒙卡方法中的问题；该方法将动力学部分与中子输运部分进行外耦合，方便使用各种蒙卡程序进行替换，具有通用性。
为达到上述目的，本发明采取了以下技术方案予以实施：
一种利用耦合蒙特卡罗方法模拟ads系统瞬态问题的方法，步骤如下：
步骤1：采用预估校正改进准静态方法处理带缓发中子先驱核的中子时空动力学方程，得到完全的预估通量方程，如公式(1)：
式中：
n--时间步序号；
ω--空间角度；
r--空间位置；
e--中子能量；
σ'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时刻的时间相关的缓发中子产生截面；
sn--n时刻的形式外源，包括实际外源，缓发中子先驱核相关源和时间相关源；
同时在预估校正改进准静态方法中，需要对预估通量进行归一化，其权重因子为共轭通量；
归一化公式为：
式中：
v--中子速度；
ψ--形状函数分布；
步骤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--迭代次数序号；
km--迭代第m步时的裂变源系数；
km-1--迭代第m-1步时的裂变源系数；
s0--归一化后的外源；
φm--迭代第m步时通量
φm-1--迭代第m-1步时的通量
l--泄漏、吸收、散射项算子；
b--裂变产生项算子；
<>--关于全相空间的积分算子；
在蒙特卡罗中子模拟中，将裂变反应视为吸收，在第m-1次计算时模拟计算相应通量，并通过式(6)计算迭代系数km-1，然后通过式(5)构造新的源项，不改变模拟问题的材料组成，进行第m次计算，反复迭代，直到两次迭代中的km-1和km的差小于10-6时，视为收敛，此时的φm即为n-1时刻的预估通量
步骤4：根据步骤3求得的n+1时刻的预估通量
式中：
ψn+1--第n+1时刻的形状函数；
由预估校正改进准静态方法，在形状函数已知的情况下，能够求出n+1时刻的点堆参数，再由点堆动力学方程能够得到n+1时刻的通量幅度tn+1；最终可得到n+1时刻的通量：
φn+1＝tn+1ψn+1公式(8)
重复步骤3步骤4继续计算n+2时刻的通量，依次类推，直到计算到目标时刻。进而得到通量随时间和空间的分布φ(r,e,t)，完成反应堆瞬态计算。
本发明与现有技术相比具有如下优点：
1、采用耦合蒙特卡罗方法，且为外部耦合，可以提高瞬态计算中各时间步上中子输运计算的准确度，另外还能灵活选取蒙特卡罗求解程序。
2、将裂变源迭代加速技巧用于蒙特卡罗求解，解决了低次临界度下固定源问题模拟缓慢且不稳定的问题，使蒙特卡罗方法更好地适用于加速器驱动次临界系统的瞬态分析计算。
附图说明
图1为利用耦合蒙特卡罗方法模拟ads系统瞬态问题的方法的计算流程图。
具体实施方式
下面结合附图和具体实施例对本发明作进一步详细说明：
本发明在原有的确定论的预估校正准静态计算方法中，将瞬态分析中分时间步的中子输运计算用蒙特卡罗方法进行耦合，通过转置裂变截面和散射截面处理共轭通量计算，通过引入裂变源迭代因子，加速并稳定蒙特卡罗中子模拟，达到预期效果。
如图1所示，本发明的具体实施步骤如下：
步骤1：采用预估校正改进准静态方法处理带缓发中子先驱核的中子时空动力学方程，得到完全的预估通量方程：
式中：
n--时间步序号；
ω--空间角度；
r--空间位置；
e--中子能量；
σ'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时刻的瞬发中子产生截面和时间相关的缓发中子产生截面；
sn--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)进行蒙特卡罗输运计算，其中外源sn由n时刻的输运计算的通量φn得到，n＝0时即为初始时刻实际外源；截面以蒙特卡罗指定材料的形式给出；通过目标为计算n+1时刻的预估通量
式中：
m--迭代次数序号；
km--迭代第m步时的裂变源系数；
km-1--迭代第m-1步时的裂变源系数；
s0--归一化后的外源；
φm--迭代第m步时通量
φm-1--迭代第m-1步时的通量
l--泄漏、吸收、散射项算子；
b--裂变产生项算子；
<>--关于全相空间的积分算子；
在蒙特卡罗中子模拟中，将裂变反应视为吸收，在第m-1次计算时模拟计算相应通量φm-1，并通过式(6)计算迭代系数km-1，然后通过式(5)构造新的源项，不改变模拟问题的材料组成，将源项改为以上计算的源进行第m次计算，反复迭代，直到两次迭代中的km-1和km的相对误差小于10-6时，视为收敛，此时的即为n+1时刻的预估通量
步骤4：
根据步骤3求得的n+1时刻的预估通量
ψn+1--第n+1时刻的形状函数；
由预估校正改进准静态方法，在形状函数已知的情况下，可以求出n+1时刻的点堆参数，由此可以得到n+1时刻的通量幅度tn+1；最终可得到n+1时刻的通量：
φn+1＝tn+1ψn+1公式(8)
重复步骤3可继续计算n+2时刻的通量，依次类推，直到计算到目标时刻；进而可得到通量随时间和空间的分布φ(r,e,t)，完成反应堆瞬态计算。