一种压水堆核热耦合的计算方法

1.本发明涉及核反应堆堆芯技术领域，具体而言，涉及一种压水堆核热耦合的计算方法。


背景技术：

2.反应核热耦合技术能够将反应堆中复杂的中子扩散方程组和热工水利方程组联立，以完成对堆芯中子和热工水力等物理场的仿真。当前，传统的耦合技术主要是通过算符分解(operator splitting，os)方法或picard迭代方法。os方法将中子扩散方程组模块和热工水利方程组模块分别求解后的部分解作为对方边界条件的新估计，在每个时间步内交替计算直到仿真结束。该方法由于在瞬态计算中各物理场参数的更新存在滞后，要求耦合时间步长较小，否则将引起较大误差。而picard迭代方法在每个时间步内，要求中子扩散方程组和热工水力方程组迭代收敛后，再推进下一个时间步。picard迭代方法是一种支持全隐式耦合的方法，能够确保每个时间步内完成收敛，可以适当增大时间步长，但是存在计算速度慢、计算稳定性较弱的问题。
3.无雅可比克雷洛夫子空间方法(jacobian-free newton
–
krylov,jfnk)方法是一种强耦合算法，其将待耦合的扩散方程组和热工水利方程组置于大型矩阵中统一求解，以保证迭代求解过程中参数的统一更新，从而提高收敛速度并减小迭代误差。一般而言，jfnk方法相较于传统的算符分解法或picard迭代法，具有更快的收敛速度、更高的收敛精度、更好的计算稳定性。因此，采用jfnk方法有利于提高压水堆核热耦合仿真的计算效率。但是在瞬态仿真中，压水堆冷却剂相变时，若统一使用两相流体方程兼容单相流体方程，则可能会出现病态，甚至退化为不定方程组而造成失真的现象。而如果在出现相变时，根据具体情况列出相应方程，再联立求解可以解决病态方程组和不定方程组的问题。但是jfnk方法本身在大型矩阵联立求解时，其计算流程较为复杂，且扩展能力不强。当物理场方程组数目和格式均发生变化时，使得整个待解矩阵发生变化，对jfnk算法而言是一项挑战。因此，如何利用jfnk方法稳定高效地求解压水堆两相核热耦合是一个待解决的技术问题。


技术实现要素：

4.本发明提供一种压水堆核热耦合的计算方法，用以克服现有技术中存在的至少一个技术问题。
5.本发明实施例提供一种压水堆核热耦合的计算方法，包括：
6.步骤s1、建立待分析的压水堆堆芯的中子扩散方程组和热工水利方程组，将所述中子扩散方程组和所述热工水利方程组进行耦合，得到核热耦合方程组；利用扩展节块法对所述中子扩散方程组进行离散化处理，得到第一网格集和，以及利用子通道法对所述热工水利方程组进行离散化处理，得到第二网格集和；建立所述第一网格集和所述第二网格集和之间的映射关系；所述第一网格集和和所述第二网格集和基于对所述压水堆堆芯采用相同的网格划分方式而得到；
7.根据所述压水堆堆芯的物理场的边界条件和初始条件，初始化所述堆芯的物理场的参数，利用少群截面库初始化包括输运截面、散射截面、吸收截面、裂变截面、中子产生截面等反应截面在内的物理量；
8.根据下述公式计算初始的时间步长，计算公式如下：
[0009][0010]
其中，符号δt表示允许的最大时间步长；符号δxi表示单元网格尺寸，单位为m；符号ui表示单元网格中的最大速度，单位为m/s；符号η表示大于0经验系数；
[0011]
步骤s2、动态调整预处理矩阵及残差方程组格式，具体包括：
[0012]
步骤s21、读取相对于待分析时刻的上一时刻的各物理参数，计算出过程参数，将所述过程参数存入计算机内存；
[0013]
步骤s22、判断所述第一网格集和和所述第二网格集和中各个网格的相态，对于其中的第一任意一个网格，若所述第一任意一个网格为液相，则该网格对应的两相残差方程组退化为液相残差方程组，对应的待解向量中汽相物理量删除；若所述第一任意一个网格为汽相，则该网格对应的两相残差方程组退化为汽相残差方程组，对应的待解向量中液相物理量删除；若所述第一任意一个网格为两相，则两相残差方程组和待解向量不变；遍历所述第一网格集和和所述第二网格集和中所有网格的相态参数，统计待解方程和待解物理量的数目，确定对应的预处理矩阵和残差方程组规模及生成形式；
[0014]
步骤s3、动态更新预处理矩阵及残差方程组，具体包括：
[0015]
步骤s31、从内存中读取所述上一时刻的各物理参数及过程参数；
[0016]
步骤s32、判断所述第一网格集和和所述第二网格集和中各个网格的相态，对于其中的第二任意一个网格，若所述第二任意一个网格为两相，则基于两相物理参数数据库填写该网格的残差方程组和待解向量，若所述第二任意一个网格为液相，基于液相物理参数数据库填写该网格的残差方程组和待解向量，若所述第二任意一个网格为汽相，基于汽相物理参数数据库填写该网格的残差方程组和待解向量，直至完成残差方程组和待解向量；
[0017]
步骤s4、结合预处理矩阵的jfnk求解核热耦合残差方程组；
[0018]
步骤s5、判断jfnk求解收敛性并调节仿真步长；
[0019]
步骤s6、判断仿真过程是否结束。
[0020]
优选的，所述若所述第二任意一个网格为两相，则基于两相物理参数数据库填写该网格的残差方程组和待解向量，具体包括：
[0021]
基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及两相的压力、温度、流速、截面含气率在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及两相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。
[0022]
优选的，所述若所述第二任意一个网格为液相，基于液相物理参数数据库填写该网格的残差方程组和待解向量，具体包括：
[0023]
基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及液相的压力、温度、流速在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及液相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。
[0024]
优选的，所述若所述第二任意一个网格为汽相，基于汽相物理参数数据库填写该网格的残差方程组和待解向量，具体包括：
[0025]
基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及汽相的压力、温度、流速在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及汽相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。
[0026]
优选的，当相态为两相时，所述核热耦合方程组为：
[0027][0028][0029]
x＝[φ,c,q,p,αg,ug,u
l
,wg,w
l
,hg,h
l
]
t
[0030]
其中，f
φ
表示中子通量残差方程组；fc表示缓发中子浓度残差方程组；fq表示裂变能量残差方程组；f
p
表示汽相质量残差方程组；表示液相质量残差方程组；表示汽相轴向动量残差方程组；表示液相轴向动量残差方程组；表示汽相横向动量残差方程组；表示液相横向动量残差方程组；表示汽相能量残差方程组；表示液相能量残差方程组；上标n+1为本时间步长的变量值，上标n为上一时间步长的变量值；δt是时间步长，单位为s；φ是中子通量密度，单位为cm-2
s-1
；v是中子运动速度，单位为cm/s；c是缓发中子先驱核的密度，单位为cm-3
；j是中子流密度，单位为cm-2
s-1
；a1、a2、a3、a4是扩展节块法相关系数；λd是缓发中子先驱核的衰变常数，单位为t-1
；σf是裂变截面，单位为cm-1
；β、k
eff
分别是缓发中子份额和有效增殖因子；ef是单次裂变产生能量，单位为j；q是裂变产生能量，单位为w；v是所述网格的体积，单位为m3；au、aw分别是正交于轴向和横向速度的流通面积，单位为m2；ρg、ρ
l
分别是汽相和液相的密度，单位为kg/m3；αg、α
l
分别是截面含汽率和截面含液率；hg、h
l
分别是汽相和液相的比焓，单位为kj/kg；ug、u
l
分别是汽相和液相的轴向流速，单位为m/s；wg、w
l
分别是汽相和液相的横向流速，m/s；γ是液相汽化为汽相的相变质量流量，单位为kg/s；pj、p
j+1
分别是轴向上游和下游网格的压力，单位为pa；p
ii
、p
jj
分别是横向上游和下游网格
的压力，单位为pa；τ
wg,u
、τ
wl,u
分别是轴向汽相和液相的壁面阻力，单位为kg
·
m/s2；τ
wg,w
、τ
wl,w
分别是横向汽相和液相的壁面阻力，单位为kg
·
m/s2；τ
gl
是汽液两相间阻力，单位为kg
·
m/s2；qg、q
l
分别是网格内传递给汽相和液相的总热量，单位为w；nj表示轴向方向相邻网格总数；n
gap
表示横向方向相邻网格总数。
[0031]
本说明书一个实施例至少能够达到以下有益效果：
[0032]
(1)本发明的压水堆核热耦合计算方法能够高效稳定地计算压水堆两相核热耦合问题，避免出现病态、不定解等问题。
[0033]
(2)本发明利用jfnk方法对压水堆堆芯进行两相核热耦合联立求解，各物理场变量在jfnk迭代中同时更新，相比于传统的os方法、picard等方法，具有更高的收敛精度、更快的收敛速度和更好的稳定性。
[0034]
(3)本发明的动态jacobi预处理方法和动态残差方程组方法，依据各离散网格中的相态以动态调整预处理矩阵和残差方程组规模及格式，规避无穷小量，减小病态，消除不定方程组，避免出现不定解。
[0035]
(4)本发明提出智能的计算时间步长方案，充分利用jfnk算法支持的全隐式迭代特性，提高仿真时间步长，减少仿真次数，大幅提升计算效率。
[0036]
(5)本发明不需要限定中子扩散方程和热工水力方程组的离散形式，仅需要用离散方程描述单元格间物理变量关系，以构建相应的jacobi预处理矩阵和差分方程组即可。
附图说明
[0037]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0038]
图1表示本发明提供的一种压水堆核热耦合的计算方法的方法流程示意图；
[0039]
图2表示本发明提供的一种压水堆核热耦合的计算方法中动态调整预处理矩阵和待解方程组格式的流程图；
[0040]
图3表示本发明提供的一种压水堆核热耦合的计算方法中动态更新预处理矩阵和待解方程组的流程图。
具体实施方式
[0041]
为使本说明书一个或多个实施例的目的、技术方案和优点更加清楚，下面将结合本说明书具体实施例及相应的附图对本说明书一个或多个实施例的技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本说明书的一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本说明书一个或多个实施例保护的范围。
[0042]
应当理解，尽管在本技术文件中可能采用术语第一、第二、第三等来描述各种信息，但这些信息不应限于这些术语。这些术语仅用来将同一类型的信息彼此区分开。
[0043]
如前文陈述，当前系统仿真由于传统的管路节点和热力构件建模方法，难以对高维、多物理场对象扩展。于是，如何高效地将各对象以各自合适的方式建模并耦合，以形成
多对象、多尺度、多物理场系统是当前的研究热点。经典的耦合方法有算符分解法(os)和picard迭代法等。采用传统的os耦合方法可能导致时间步长内各物理场参数不收敛，使得精细的模型算出来不精确甚至错误的结果。而采用picard方法可以使得各步长内物理场参数收敛，但是计算量较大、收敛稳定性较差、计算效率较低。因此，兼容多对象、多尺度、多物理场的高收敛速度高稳定性的系统仿真方法值得研究。本发明提供一种压水堆核热耦合的计算方法，用以解决现有技术中的技术问题。
[0044]
如图1所示，该仿真方法流程可以包括如下几个步骤：
[0045]
步骤s1：初始化方程组、时间步长和物理场参数：
[0046]
步骤s11：初始化方程组。本发明可以将中子扩散方程组和热工水力方程组合并为一个大型核热耦合方程组。以压水堆典型堆芯建模为例，可以利用扩展节块法对中子扩散方程组进行离散化处理，得到第一网格集和，以及，利用子通道法对热工水利方程组进行离散化处理，得到第二网格集和，然后建立该第一网格集和该第二网格集和之间的映射关系。在进行具体的网格划分时，可以基于对压水堆堆芯采用相同的网格划分方式而得到第一网格集和和第二网格集和，从而便于建立扩展节块法的离散网格和子通道的离散网格之间的映射关系。其中，当相态为两相的所述核热耦合方程如下所示：
[0047][0048][0049]
x＝[φ,c,q,p,αg,ug,u
l
,wg,w
l
,hg,h
l
]
t
[0050]
其中，f
φ
表示中子通量残差方程组；fc表示缓发中子浓度残差方程组；fq表示裂变能量残差方程组；f
p
表示汽相质量残差方程组；表示液相质量残差方程组；表示汽相轴向动量残差方程组；表示液相轴向动量残差方程组；表示汽相横向动量残差方程组；表示液相横向动量残差方程组；表示汽相能量残差方程组；表示液相能量残差方程组；上标n+1为本时间步长的变量值，上标n为上一时间步长的变量值；δt是时间步长，单位
为s；φ是中子通量密度，单位为cm-2
s-1
；v是中子运动速度，单位为cm/s；c是缓发中子先驱核的密度，单位为cm-3
；j是中子流密度，单位为cm-2
s-1
；a1、a2、a3、a4是扩展节块法相关系数；λd是缓发中子先驱核的衰变常数，单位为t-1
；σf是裂变截面，单位为cm-1
；β、k
eff
分别是缓发中子份额和有效增殖因子；ef是单次裂变产生能量，单位为j；q是裂变产生能量，单位为w；v是所述网格的体积，单位为m3；au、aw分别是正交于轴向和横向速度的流通面积，单位为m2；ρg、ρ
l
分别是汽相和液相的密度，单位为kg/m3；αg、α
l
分别是截面含汽率和截面含液率；hg、h
l
分别是汽相和液相的比焓，单位为kj/kg；ug、u
l
分别是汽相和液相的轴向流速，单位为m/s；wg、w
l
分别是汽相和液相的横向流速，m/s；γ是液相汽化为汽相的相变质量流量，单位为kg/s；pj、p
j+1
分别是轴向上游和下游网格的压力，单位为pa；p
ii
、p
jj
分别是横向上游和下游网格的压力，单位为pa；τ
wg,u
、τ
wl,u
分别是轴向汽相和液相的壁面阻力，单位为kg
·
m/s2；τ
wg,w
、τ
wl,w
分别表示横向汽相和液相的壁面阻力，单位为kg
·
m/s2；τ
gl
是汽液两相间阻力，单位为kg
·
m/s2；qg、q
l
分别是网格内传递给汽相和液相的总热量，单位为w；nj表示轴向方向相邻网格总数；n
gap
表示横向方向相邻网格总数。
[0051]
步骤s12：初始化物理场参数。根据堆芯物理场的边界条件和初始条件，初始化堆芯物理场的参数，再利用少群宏观截面库初始化反应截面等物理量。堆芯物理场边界条件和初始条件一般为流速、压力、温度、比焓等物理量，本例选为堆芯入口处流速、比焓和堆芯出口出为压力。少群宏观反应截面一般包括输运截面、散射截面、吸收截面、裂变截面、中子产生截面等。少群宏观截面库能够根据流体的密度、温度、硼浓度以及燃料的温度等信息查询出相应的截面数据。
[0052]
步骤s13：初始化时间步长。根据库朗数计算初始的时间步长，并存于内存之中，以便后续计算调用。时间步长的计算公式如下：
[0053][0054]
其中，δt是允许的最大时间步长；δxi是单元网格尺寸，m；ui是单元网格中的最大速度，m/s；η是经验系数，大于0，当方程组为全隐格式时，可以取值为3甚至更高。
[0055]
步骤s2：动态调整预处理矩阵及残差方程组格式，如图2所示，具体可以包括：
[0056]
步骤s21：读取上一时刻的解向量参数，并计算出过程参数，并存入内存方便别的步骤读取。因为核热耦合是多物理场的仿真，考虑到网格数目和方程组待解量数目，若将所有的物理量作为待解量并写成大型待解向量，则残差方程和预处理矩阵将非常巨大，计算困难。而本发明将燃料棒、控制棒导管、定位格架等固体结构件的导热方程组作为过程参数方程组，有利于保证计算精度同时减少计算量，提高计算效率。因为过程参数方程组将在每次jfnk核热耦合的非精确牛顿迭代后完成更新求解，并参与下一次非精确牛顿迭代直至收敛，待解向量和过程参数都完成了更新，满足全隐式耦合的要求，计算精度较小。而将过程参数方程组独立求解，有利于减少核热耦合方程组规模，大幅减少jfnk核热耦合计算量，有利于提高计算效率。
[0057]
步骤s22：确定预处理矩阵及残差方程组的组成格式。由于瞬态计算中各单元格中相态会发生变化而导致残差方程组数目及格式变化，以及对应的待解向量发生变化。从而使得预处理矩阵也发生变化。因此，需要遍历所有网格的相态参数，统计待解方程和待解物
理量的数目，确定对应的预处理矩阵和残差方程组规模及生成形式。
[0058]
本步骤中，可以判断第一网格集和和第二网格集和中各个网格的相态，对于其中的第一任意一个网格，若该第一任意一个网格为液相，则该网格对应的两相残差方程组退化为液相残差方程组，对应的待解向量中汽相物理量删除；若该第一任意一个网格为汽相，则该网格对应的两相残差方程组退化为汽相残差方程组，对应的待解向量中液相物理量删除；若该第一任意一个网格为两相，则两相残差方程组和待解向量不变；遍历所述第一网格集和和所述第二网格集和中所有网格的相态参数，统计待解方程和待解物理量的数目，确定对应的预处理矩阵和残差方程组规模及生成形式。
[0059]
更具体的，本发明可以根据上一步每个网格的截面含汽率数值结果，判断该网格处于何种相态，并动态调整相应的方程组形式。若为两相，则两相残差方程组和待解向量不变。若为液相，则该网格对应的两相残差方程组退化为液相残差方程组，对应待解向量中汽相物理量删除，如下式所示：
[0060][0061]
同理，若为汽相，则该网格对应的两相残差方程组退化为汽相残差方程组，对应待解向量中液相物理量删除，如下式所示：
[0062][0063]
本发明采用jacobi矩阵作为jfnk预处理矩阵，该预处理矩阵由所有网格的离散方程组和待解变量经过求偏导生成，其块矩阵形式如下：
[0064][0065]
当某网格相态为单相时，其残差方程仅为单相的方程。在形成jacobi矩阵时，使用所述单相残差方程对其自变量求偏导以得到偏导项，并将偏导项置于所述jacobi矩阵相应位置即可。
[0066]
由此可见，本发明在存在单相状态的网格时，由于两相方程组退化成了单相方程，整体的残差方程组和预处理矩阵规模会减小，因此减小了计算量，能够提升计算速度。
[0067]
步骤s3：动态更新预处理矩阵及残差方程组，如图3所示，具体可以包括：
[0068]
步骤s31：从内存中读取所述上一时刻的解向量参数、过程参数以及时间步长。
[0069]
步骤s32：根据第2步中生成的残差方程组格式，更新残差方程组参数。利用所述从内存中读取的上一时刻的解向量参数、过程参数和时间步长，更新残差方程组。
[0070]
具体的，如图3所示，若第一网格集和和第二网格集和和中的第二任意一个网格为两相，则基于两相物理参数数据库填写该网格的残差方程组和待解向量，具体可以包括：
[0071]
基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及两相的压力、温度、流速、截面含气率在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及两相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。
[0072]
若该第二任意一个网格为液相，基于液相物理参数数据库填写该网格的残差方程组和待解向量，具体可以包括：
[0073]
基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及液相的压力、温度、流速在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及液相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。
[0074]
若该第二任意一个网格为汽相，基于汽相物理参数数据库填写该网格的残差方程组和待解向量，具体可以包括：
[0075]
基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及汽相的压力、温度、流速在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及汽相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。
[0076]
步骤s33：根据第2步中生成的预处理矩阵格式，更新预处理矩阵。利用所述从内存中读取的上一时刻的解向量参数、过程参数和时间步长，更新预处理矩阵。本例采用jacobi矩阵作为预处理矩阵。jacobi矩阵是由各残差方程组对待解向量中所有元素求偏导形成的大型矩阵。其格式如下：
[0077][0078]
该jacobi矩阵用于jfnk中krylov子空间的预处理计算，以提高krylov子空间法的计算效率。如果使用左预处理，则将该jacobi矩阵经不完全lu分解后取逆运算，再左乘于待解矩阵上。
[0079]
步骤s4：结合预处理矩阵，利用jfnk求解核热耦合残差方程组包括：
[0080]
步骤s41：非精确牛顿迭代。非精确牛顿迭代法是jfnk迭代求解的外循环，非精确牛顿迭代包括：构建jacobi矩阵向量积线性方程组、迭代解列的更新和迭代收敛的判断等过程。其中，迭代解列更新包括对krylov子空间解的更新以及过程参数的更新。迭代收敛判断指在外循环中通过判断方程残差是否满足收敛要求，若满足要求则退出循环，若不满足要求则继续外循环。构建jacobi矩阵向量积线性方程组指将非线性方程组按照非精确牛顿方程的形式变形，并对jacobi矩阵向量积差分，形成线性方程组，以便krylov子空间迭代求解。非精确牛顿方程如下式所示：
[0081][0082]
其中，f(xk)是第k次迭代函数值，f
′
(xk)是第k次迭代的jacobi矩阵值，是需要求解的迭代步长，ηk是强制项。但是，由于jacobi矩阵求解较为困难，且非精确牛顿方程中jacobi矩阵是与迭代步长以矩阵向量积的形式出现，所以jfnk方法采用差分近似处理如下：
[0083][0084]
其中，ε是较小的常数。由此便将复杂的非精确牛顿方程转换为线性方程组，便于下一步利用krylov子空间方法高效求解方程组。
[0085]
步骤s42：krylov子空间迭代求解线性方程组。在jfnk求解过程中，krylov子空间迭代法是jfnk迭代求解的内循环算法，用于求解大型线性方程组。本发明使用第3步的预处理矩阵并结合krylov子空间迭代法便可以稳定高效地求解出迭代步长。记左预处理矩阵为m，经左预处理后的非精确牛顿方程为：
[0086][0087]
上式在一个迭代步中，除了迭代步长是未知数外，其他数均已知。使用krylov子空间迭代法便可以以较高的收敛速度完成迭代步的计算。
[0088]
步骤s5：判断jfnk求解收敛性并调节仿真步长包括：
[0089]
jfnk求解收敛性的判断。在jfnk求解过程中，若预设的迭代次数(如在一种可能的实施例技术方案中，可以将该迭代次数设置为1千次)并未收敛，则认定为方程求解不收敛。当方程计算不收敛时，将当前仿真时间步长除以2，并更新内存中的时间步长值，重新执行第3步(动态更新预处理矩阵及残差方程组)，以更新预处理矩阵及残差方程组。若方程计算收敛，则认为完成本仿真时刻的物理场参数计算，进行下一步。
[0090]
步骤s6：对仿真结束时间的判断，并做相应仿真进程调整，具体可以包括：
[0091]
对仿真结束时间的判断：比较当前仿真时间和需求的仿真时间，若未达到需求的仿真时间，则根据库朗数计算仿真步长，执行第2步(确定预处理矩阵及残差方程组的组成形式)。若达到需求的仿真时间则程序终止。
[0092]
本技术技术方案中在使用jfnk方法求解前述大型核热耦合方程组过程时，中子场、温度场、压力场、流速场等物理量将在每个非精确牛顿迭代步中同时更新，相比于将物理场分开求解的os方法及picard方法，有更高的核热耦合计算精度。而且相比于os方法及
picard方法的一阶收敛速度，本发明采用的jfnk方法具有二阶收敛速度，有更快的计算速度。此外，本发明根据离散网格的相态动态调整对应方程组，相比对各网格相态不加区分地使用两相方程建模，本发明方程个数更少，计算量更小，且能够避免两相方程计算单相物理场时产生的病态问题。
[0093]
本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
[0094]
本领域普通技术人员可以理解：实施例中的装置中的模块可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块。
[0095]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。

技术特征：
1.一种压水堆核热耦合的计算方法，其特征在于，所述仿真方法包括：步骤s1、建立待分析的压水堆堆芯的中子扩散方程组和热工水利方程组，将所述中子扩散方程组和所述热工水利方程组进行耦合，得到核热耦合方程组；利用扩展节块法对所述中子扩散方程组进行离散化处理，得到第一网格集和，以及利用子通道法对所述热工水利方程组进行离散化处理，得到第二网格集和；建立所述第一网格集和所述第二网格集和之间的映射关系；所述第一网格集和和所述第二网格集和基于对所述压水堆堆芯采用相同的网格划分方式而得到；根据所述压水堆堆芯的物理场的边界条件和初始条件，初始化所述堆芯的物理场的参数，利用少群截面库初始化包括输运截面、散射截面、吸收截面、裂变截面、中子产生截面等反应截面在内的物理量；根据下述公式计算初始的时间步长，计算公式如下：其中，符号δt表示允许的最大时间步长；符号δx
i
表示单元网格尺寸，单位为m；符号u
i
表示单元网格中的最大速度，单位为m/s；符号η表示大于0经验系数；步骤s2、动态调整预处理矩阵及残差方程组格式，具体包括：步骤s21、读取相对于待分析时刻的上一时刻的各物理参数，计算出过程参数，将所述过程参数存入计算机内存；步骤s22、判断所述第一网格集和和所述第二网格集和中各个网格的相态，对于其中的第一任意一个网格，若所述第一任意一个网格为液相，则该网格对应的两相残差方程组退化为液相残差方程组，对应的待解向量中汽相物理量删除；若所述第一任意一个网格为汽相，则该网格对应的两相残差方程组退化为汽相残差方程组，对应的待解向量中液相物理量删除；若所述第一任意一个网格为两相，则两相残差方程组和待解向量不变；遍历所述第一网格集和和所述第二网格集和中所有网格的相态参数，统计待解方程和待解物理量的数目，确定对应的预处理矩阵和残差方程组规模及生成形式；步骤s3、动态更新预处理矩阵及残差方程组，具体包括：步骤s31、从内存中读取所述上一时刻的各物理参数及过程参数；步骤s32、判断所述第一网格集和和所述第二网格集和中各个网格的相态，对于其中的第二任意一个网格，若所述第二任意一个网格为两相，则基于两相物理参数数据库填写该网格的残差方程组和待解向量，若所述第二任意一个网格为液相，基于液相物理参数数据库填写该网格的残差方程组和待解向量，若所述第二任意一个网格为汽相，基于汽相物理参数数据库填写该网格的残差方程组和待解向量，直至完成残差方程组和待解向量；步骤s4、结合预处理矩阵的jfnk求解核热耦合残差方程组；步骤s5、判断jfnk求解收敛性并调节仿真步长；步骤s6、判断仿真过程是否结束。2.根据权利要求1所述的压水堆核热耦合的计算方法，其特征在于，所述若所述第二任意一个网格为两相，则基于两相物理参数数据库填写该网格的残差方程组和待解向量，具体包括：
基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及两相的压力、温度、流速、截面含气率在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及两相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。3.根据权利要求1所述的压水堆核热耦合的计算方法，其特征在于，所述若所述第二任意一个网格为液相，基于液相物理参数数据库填写该网格的残差方程组和待解向量，具体包括：基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及液相的压力、温度、流速在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及液相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。4.根据权利要求1所述的压水堆核热耦合的计算方法，其特征在于，所述若所述第二任意一个网格为汽相，基于汽相物理参数数据库填写该网格的残差方程组和待解向量，具体包括：基于从所述内存中读取出的关于所述上一时刻的包括中子通量密度以及汽相的压力、温度、流速在内的物理参数，计算出包括中子流密度、宏观反应截面、总热功率密度以及汽相的换热系数、壁面热功率、壁面阻力系数在内的过程参数。5.根据权利要求1所述的压水堆核热耦合的计算方法，其特征在于，当相态为两相时，所述核热耦合方程组为：所述核热耦合方程组为：x＝[φ,c,q,p,α
g
,u
g
,u
l
,w
g
,w
l
,h
g
,h
l
]
t
其中，f
φ
表示中子通量残差方程组；f
c
表示缓发中子浓度残差方程组；f
q
表示裂变能量残差方程组；f
p
表示汽相质量残差方程组；表示液相质量残差方程组；表示汽相轴向动量残差方程组；表示液相轴向动量残差方程组；表示汽相横向动量残差方程组；
表示液相横向动量残差方程组；表示汽相能量残差方程组；表示液相能量残差方程组；上标n+1为本时间步长的变量值，上标n为上一时间步长的变量值；δt是时间步长，单位为s；φ是中子通量密度，单位为cm-2
s-1
；v是中子运动速度，单位为cm/s；c是缓发中子先驱核的密度，单位为cm-3
；j是中子流密度，单位为cm-2
s-1
；a1、a2、a3、a4是扩展节块法相关系数；λ
d
是缓发中子先驱核的衰变常数，单位为t-1
；σ
f
是裂变截面，单位为cm-1
；β、k
eff
分别是缓发中子份额和有效增殖因子；e
f
是单次裂变产生能量，单位为j；q是裂变产生能量，单位为w；v是所述网格的体积，单位为m3；a
u
、a
w
分别是正交于轴向和横向速度的流通面积，单位为m2；ρ
g
、ρ
l
分别是汽相和液相的密度，单位为kg/m3；α
g
、α
l
分别是截面含汽率和截面含液率；h
g
、h
l
分别是汽相和液相的比焓，单位为kj/kg；u
g
、u
l
分别是汽相和液相的轴向流速，单位为m/s；w
g
、w
l
分别是汽相和液相的横向流速，m/s；γ是液相汽化为汽相的相变质量流量，单位为kg/s；p
j
、p
j+1
分别是轴向上游和下游网格的压力，单位为pa；p
ii
、p
jj
分别是横向上游和下游网格的压力，单位为pa；τ
wg,u
、τ
wl,u
分别是轴向汽相和液相的壁面阻力，单位为kg
·
m/s2；τ
wg,w
、τ
wl,w
分别是横向汽相和液相的壁面阻力，单位为kg
·
m/s2；τ
gl
是汽液两相间阻力，单位为kg
·
m/s2；q
g
、q
l
分别是网格内传递给汽相和液相的总热量，单位为w；n
j
表示轴向方向相邻网格总数；n
gap
表示横向方向相邻网格总数。

技术总结
本发明公开一种压水堆核热耦合的计算方法，方案可以包括：建立待分析的压水堆堆芯的中子扩散方程组和热工水利方程组，将所述中子扩散方程组和所述热工水利方程组进行耦合，得到核热耦合方程组；动态调整预处理矩阵及残差方程组格式；动态更新预处理矩阵及残差方程组；结合预处理矩阵的JFNK求解核热耦合残差方程组；判断JFNK求解收敛性并调节仿真步长；判断仿真过程是否结束。断仿真过程是否结束。断仿真过程是否结束。


技术研发人员：李磊 张宇航 田兆斐
受保护的技术使用者：哈尔滨工程大学
技术研发日：2023.04.28
技术公布日：2023/8/13