提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法与流程

1.本发明属于核反应堆系统安全分析计算技术领域，具体涉及一种提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法。


背景技术：

2.钠冷快堆由于冷却剂载热性能好、热惰性大、温度水平高，因此一般采用直流蒸汽发生器来减小反应堆尺寸，提高经济性。三回路直流蒸汽发生器管侧为水，经二回路液态钠加热，与壁面的换热模式包括过冷水对流、过冷沸腾、核态沸腾、膜态沸腾和过热蒸汽对流。各模式换热系数相差很大，尤其是chf点前后的核态沸腾和膜态沸腾。固定网格内可能出现两种不同的换热模式，若只选取其中一种模式进行换热计算，会引起换热量的计算误差。同时，如果网格的换热模式发生变化，那么由于差分方程的变化会人为造成解的不连续。这种固定网格内换热模式转变引起的热量剧烈变换，将导致直流蒸汽发生器模拟的震荡问题。
3.目前直流蒸汽发生器换热计算方法有两种，一种是固定网格法，另一种是移动边界法。固定网格法是指在计算前将系统划分多个网格，在瞬态计算过程中网格不变，利用守恒方程计算流体的热力学状态。固定网格法适用于系统分析程序，但会造成解的不连续。可以通过增加网格的数量和减少时间不长的方法降低不连续性，但会降低计算效率。移动边界法中，网格边界随换热模式变化，因此网格界面是时间的连续函数，避免计算结果的突跳。同时，移动边界法可以减少网格数目，提高计算效率，但移动边界法不适用于系统分析程序。
4.钠冷快堆系统分析程序中直流蒸汽发生器仍采用传统固定网格法计算换热，计算效率和计算结果稳定性难以两全。亟待为钠冷快堆系统分析程序开发一种快速准确求解直流蒸汽发生器换热的方法，提升钠冷快堆系统程序瞬态计算的精度和速度。


技术实现要素：

5.本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法，既保留系统程序的所有特点，又提升钠冷快堆系统程序瞬态计算的精度和速度。
6.本发明采用以下技术方案：
7.提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法，包括以下步骤：
8.s1、根据蒸汽发生器水力部件控制体划分方式，采用固定网格法确定蒸汽发生器热构件网格划分，计算chf位置；
9.s2、判定各网格换热模式，针对单相液换热、核态沸腾、膜态沸腾和过热蒸汽工况，分别计算各网格换热系数h；
10.s3、基于固定网格，使用chf可移动边界模型对发生chf现象的控制体再次网格划分，重新计算上、下热构件的换热系数；
11.s4、设置边界条件，重复以上步骤，分别计算每个网格下一时间步换热系数，最终
求得每一时刻蒸汽发生器的换热量。
12.具体的，步骤s1中，chf点在网格中的具体位置为：
[0013][0014]
其中，p表示钠侧，δl为chf位置距前一个网格的距离，t
p,chf
为发生chf时钠侧控制体温度，t
s
为水侧饱和温度，t
p,i
为钠侧控制体i的温度，u为综合传热系数，r
i
为传热管内半径，w
p
为钠侧流体质量流量，c
p
为钠的比热容。
[0015]
具体的，步骤s2中，固定网格水侧换热模式判定采用平衡含汽率xe和chf位置，逻辑如下：
[0016]
(1)xe≤0，网格为单相对流换热；
[0017]
(2)0＜xe＜1，且网格在chf发生前，网格为核态沸腾换热；
[0018]
(3)0＜xe＜1，且网格在chf发生后，网格为膜态沸腾换热；
[0019]
(4)xe≥1，网格为过热蒸汽换热。
[0020]
进一步的，水侧单相对流换热关系式根据流型选择，具体为：
[0021]
对于层流，即re<2300，选择sellars关系式进行计算：
[0022][0023]
对于湍流，即re≥2300，选取dittus-boelter公式进行计算：
[0024]
nu＝0.023re
0.8
pr
0.4
[0025]
其中，k
f
为液相热导率；d
h
为管道的水力直径，h为内表面传热系数，nu为努塞尔数，re为雷诺数，pr为普朗特数。
[0026]
进一步的，水侧核态沸腾换热采用chen关系式，具体为：
[0027]
h
tp
＝fh
mac
+sh
mic
[0028]
其中，h
mac
为宏观对流项换热系数；f为雷诺数因子；h
mic
为微观沸腾项换热系数；s为抑制因子。
[0029]
更进一步的，对于饱和核态沸腾，雷诺数因子f计算如下：
[0030][0031][0032]
对于过冷核态沸腾，雷诺数因子f计算如下：
[0033][0034]
其中，g为流体质量流速，ρ
f
为液相密度，ρ
g
为汽相密度，g
g
为液相流量，g
f
为汽相流量，μ
g
为液相粘性系数，μ
f
为汽相粘性系数，t
s
为水侧饱和温度，t
f
为液相温度。
[0035]
进一步的，水侧膜态沸腾换热关系式根据流量进行选择，具体为：
[0036]
g≤271.24kg/(m2·
s)，采用bromely关系式计算如下：
[0037][0038]
g＞271.24kg/(m2·
s)，采用groneveld关系式计算如下：
[0039][0040]
其中，h为内表面传热系数，λ
g
为汽相导热系数，g为重力加速度，ρ
g
为汽相密度，μ
g
为液相粘性系数，ρ
f
为液相密度，t
w
为壁面温度，t
s
为水侧饱和温度，h'
fg
为汽化比焓，d为管道直径，g为流体质量流速，c
pg
为汽相比热容，y为修正系数，x为平衡含汽率，d
e
为当量直径。
[0041]
进一步的，水侧过热蒸汽换热采用sider-tate关系式，具体为：
[0042][0043]
其中，nu为努塞尔数，re为雷诺数，pr为普朗特数，μ为粘性系数，μ
w
为壁面粘性系数。
[0044]
具体的，步骤s2中，钠侧换热系数采用westinghouse棒束关系式：
[0045][0046]
适用于任意棒束排布，且
[0047][0048]
适用于任意棒束排布，且
[0049][0050]
适用于任意棒束排布，且
[0051]
具体的，步骤s3中，采用移动边界法对发生chf的网格再划分，将随时间变化的chf位置作为移动边界，网格边界下为核态沸腾，边界上为膜态沸腾，并重新计算边界上、下区域各自的换热。
[0052]
与现有技术相比，本发明至少具有以下有益效果：
[0053]
本发明一种提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法，使用固定网格加移动边界法计算直流蒸汽发生器的换热，保留了固定网格在系统分析程序中的优势，可以通过差分形式计算热工水力特性。以chf位置作为移动边界，并只对发生chf网格重新计算换热，极大地简化了单一移动边界法在系统程序中的编程困难问题，同时有效减弱传统固定网格法造成的计算结果不稳定性，对移动边界上下区域重新计算换热，提高粗网格的计算精度，同时提高了计算效率。
[0054]
进一步的，确定chf点在网格中的具体位置，热构件的网格划分依据水力部件控制体划分方式，满足钠冷快堆系统分析程序划分方式。
[0055]
进一步的，水侧换热逻辑依据平衡含汽率xe和chf位置，简化了固定网格水侧换热模式判定过程。
[0056]
进一步的，通过水侧单相对流换热关系式考虑了层流和湍流的影响。
[0057]
进一步的，水侧核态沸腾采用更精确的chen关系式。
[0058]
进一步的，水侧膜态沸腾关系式分别考虑了高低流量的影响。
[0059]
进一步的，水侧过热蒸汽关系式可通过粘性系数比、re和pr直接求得。
[0060]
进一步的，钠侧换热系数关系式考虑了棒束，且棒束可以任意排布。
[0061]
进一步的，在固定网格下再次划分移动边界，可以增加网格的连续性。
[0062]
综上所述，本发明采用固定网格划分蒸汽发生器热构件，并在发生chf的网格划分可移动边界，针对不同的工况，壁面换热关系式采用了不同的计算方式，最后可以更快地得到更精确的钠冷快堆蒸汽发生器的壁面换热量。
[0063]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0064]
图1为本发明示意图，其中(a)为固定网格划分，(b)为固定网格加移动边界划分；
[0065]
图2为本发明流程图；
[0066]
图3为chf轴向位置随时间变化的曲线图。
具体实施方式
[0067]
本发明提供了一种提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法，采用固定网格法划分钠冷快堆直流蒸汽发生器，采用固定网格法确定蒸汽发生器热构件网格划分，计算chf发生的位置。判定各网格换热模式，针对单相液换热、核态沸腾、膜态沸腾和过热蒸汽工况，分别计算各网格换热系数。基于固定网格，使用chf可移动边界模型对发生chf现象的控制体再次网格划分，重新计算上、下热构件的换热系数。设置边界条件，重复以上步骤，分别计算每个网格下一时间步换热系数，最终求得每一时刻蒸汽发生器的换热量。既保证钠冷快堆系统分析程序中直流蒸汽发生器换热计算精度，又提升反应堆数值计算程序的求解速度。
[0068]
请参阅图2，本发明一种提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法，包括以下步骤：
[0069]
s1、根据蒸汽发生器水力部件控制体划分方式，采用固定网格法确定蒸汽发生器热构件网格划分，计算chf位置；
[0070]
首先建立蒸汽发生器正常运行时一、二次侧工质换热量计算公式：
[0071][0072]
其中，下标i为管道内工质参数；下标o为管道外工质参数；t为流体温度/k；h为内表面传热系数/j
·
m-2
·
k-1
；l为控制体长度/m；λ为传热管热导率/j
·
m-1
·
k-1
。
[0073]
对上式进行进一步简化：
[0074][0075][0076]
其中，
[0077][0078]
根据上式写出发生chf时，蒸汽发生器一、二次侧工质换热量的计算关系式：
[0079][0080]
其中，u为综合传热系数/j
·
m-2
·
k-1
；t
p.chf
为发生chf时钠侧控制体温度/k；t
s
为水侧饱和温度/k。
[0081]
其中，综合传热系数u为：
[0082][0083]
当蒸汽发生器达到稳定状态时，单位时间内钠侧工质能量的减少等于水侧工质能量增加，如下式所示：
[0084]
c
p
w
p
(t
p.chf-t
pi
)＝w
s
(h
s,chf-h
si
)
[0085]
其中，c
p
为钠的比热容/j
·
kg-1
·
k-1
；t
p.chf
为发生chf时钠侧温度/k；w
p
为钠侧流体质量流量/kg
·
s-1
；t
pi
为钠侧控制体i的温度/k；w
s
为钠侧工质质量流量/kg
·
s-1
；h
s,chf
为发生chf时水侧工质焓值/j
·
kg-1
；h
si
为水侧控制体i的焓值/j
·
kg-1
。
[0086]
水在两相区的焓值变化还可以用含汽率来表示：
[0087]
w
s
(h
chf-h
i
)＝w
s
h
fg
(x
chf-x
i
)
[0088]
水侧含汽率与钠侧温度的关系：
[0089]
w
p
c
p
(t
p,chf-t
p,i
)＝w
s
h
fg
(x
chf-x
i
)
[0090]
其中，h
fg
为水的汽化潜热/j
·
kg-1
；x
chf
为发生chf时水侧的含汽率；x
i
为水侧控制体i的含汽率。
[0091]
对上式进行化简：
[0092][0093]
可以得到发生chf时，水侧含汽率关于钠侧温度的关系式：
[0094][0095]
根据chf查询表，临界热流密度q
chf
与含汽率x
chf
，流体压力和流体质量流速g有关，如下式所示：
[0096]
q
chf
＝f(g,x
chf
,p)
[0097]
迭代求解得到临界热流密度q
chf
，计算出chf点处的钠侧温度t
p.chf
。
[0098]
将换热计算关系式转成微分形式：
[0099]
dq＝w
p
c
p
dt＝-2πu(t
p-t
s
)r
i
dl
[0100]
并沿蒸汽发生器轴向进行积分，得：
[0101][0102]
其中，t1、t2为蒸汽发生器钠侧两端温度/k；t
s
为水侧饱和温度/k；w
p
为钠侧质量流量/kg
·
s-1
；r
i
为传热管内半径/m；u为蒸汽发生器传热管一、二次侧之间的综合传热系数/w
·
(m2·
k)-1
；c
p
为钠侧比热容/j
·
(kg
·
k)-1
；δl为水侧两端距离/m。
[0103]
如果t1＝t
p.chf
，t2＝t
p.i
，那么δl就是chf发生位置与前一个控制体i之间的距离，这样就可以计算出chf实际位置。接下来将发生chf的控制体进行进一步网格划分，细化该控制体的换热过程，从而提高计算精度。
[0104]
请参阅图1，(a)图为采用固定网格划分的节点，chf边界只能在控制体边界上，(b)图采用固定网格和可移动边界划分节点，chf边界可以出现在控制体内。
[0105]
s2、针对不同工况，换热关系式采用不同的计算方式；
[0106]
对钠侧采用westinghouse棒束关系式：
[0107][0108]
适用于任意棒束排布，且
[0109][0110]
适用于任意棒束排布，且
[0111][0112]
适用于任意棒束排布，且
[0113]
对水侧，则进行换热模式判断，再根据不同工况选取换热系数计算公式。
[0114]
(1)xe≤0，网格为单相对流换热；
[0115]
对于层流，即re<2300，选择sellars关系式进行计算：
[0116][0117]
其中，k
f
为液相热导率/j
·
(m
·
k)-1；d
h
为管道的水力直径/m。
[0118]
对于湍流，即re≥2300，选取dittus-boelter公式进行计算：
[0119]
nu＝0.023re
0.8
pr
0.4
[0120]
(2)0＜xe＜1，且网格在chf发生前，网格为核态沸腾换热；
[0121]
核态沸腾换热采用chen关系式：
[0122]
h
tp
＝fh
mac
+sh
mic
[0123]
其中，h
mac
为宏观对流项换热系数/j
·
m-2
·
k-1
；f为雷诺数因子；h
mic
为微观沸腾项换热系数/j
·
m-2
·
k-1
；s为抑制因子。
[0124]
其中，宏观对流换热系数h
mac
计算公式为：
[0125][0126]
微观微观沸腾项换热系数为h
mic
计算公式为：
[0127][0128]
其中，cp
f
为液相比热容/j
·
(kg
·
k)-1
；ρ
f
为液相密度/kg
·
m-3
；σ为流体表面张力/n
·
m-1
；ρ
g
为汽相密度/kg
·
m-3
。
[0129]
对于饱和核态沸腾，雷诺数因子f计算关系式如下：
[0130][0131][0132]
其中，g为流体质量流速/kg
·
m-2
·
s-1
。
[0133]
对于过冷核态沸腾，雷诺数因子f计算关系式如下：
[0134][0135]
抑制因子计算关系式如下：
[0136][0137]
其中，
[0138]
re
tp
＝min(10-4
re
f
f
1.25
,70)
[0139][0140]
(3)0＜xe＜1，且网格在chf发生后，网格为膜态沸腾换热；
[0141]
膜态沸腾换热关系式根据流量进行选择。
[0142]
对小流量(g≤271.24kg/(m2·
s))采用bromely关系式：
[0143][0144]
对大流量(g＞271.24kg/(m2·
s))采用groneveld关系式：
[0145][0146]
(4)xe≥1，网格为过热蒸汽换热。
[0147]
水侧过热蒸汽换热采用sider-tate关系式如下：
[0148][0149]
s3、基于固定网格，使用chf可移动边界模型对发生chf现象的控制体再次网格划分，将随时间变化的chf位置作为移动边界，网格边界下为核态沸腾，边界上为膜态沸腾，并重新计算边界上、下区域各自的换热，关系式选择方式与步骤s2中相同。
[0150]
s4、设置边界条件，重复以上步骤，分别计算每个网格下一时间步换热系数，最终求得每一时刻蒸汽发生器的换热量。
[0151]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0152]
根据etec实验的相关参数，选取了a2、a5、a7和a11实验工况采用nusol-lmr程序进行了模拟。表1给出了水侧出口温度在是否考虑chf可移动边界模型前后的对比。表2为考虑chf可移动边界模型前、后程序计算etec稳态工况的水侧进出口焓差与实验测量值的对比表。表3为本程序分别使用固定网格和chf可移动边界模型计算etec蒸汽发生器时在不同网格数目下计算精度和计算时间的对比表。
[0153]
表1考虑chf可移动边界模型前、后程序计算的水侧出口温度对比
[0154][0155]
表2考虑chf可移动边界模型前、后程序计算的水侧进出口焓差与实验值的对比
[0156][0157]
表3在不同网格数目下使用两种网格划分方法计算etec蒸汽发生器的结果对比表
[0158][0159][0160]
通过将本发明方法的结果与实验结果进行对比，使用了chf可移动边界法之后，程序计算结果与实验测量值之间的误差得到减小，尤其在控制体中chf发生的位置靠后，使用chf可移动边界模型后，控制体的换热量得到较大提升，蒸汽发生器的计算结果有较大改进。增加蒸汽发生器轴向节点数目和使用chf可移动边界模型都可以提高程序的计算结果的准确性，但chf可移动边界模型可以提高程序的计算速度。
[0161]
图3为固定网格和chf可移动边界模型下程序计算etec蒸汽发生器一次侧工质入口温度由450k升高到770k时chf轴向位置随时间变化曲线。该瞬态工况下固定网格中chf位置在控制体边界上发生移动时会在较大范围内发生反复振荡。而使用了chf可移动边界模型之后，因为把chf边界定义在了控制体内部，所以当蒸汽发生器工质参数发生变化时，程序可以在每一个时间步长上都计算出chf发生时其边界在控制体中的具体位置，chf边界在控制体中的位置变化是连续的，不会发生突跳。因此chf可移动边界模型会增加程序的稳定性，不易发生振荡。
[0162]
综上所述，本发明一种提高钠冷快堆系统程序中蒸汽发生器换热计算性能的方法，结合了系统程序固定网格划分方式，采用chf可移动边界模型计算壁面换热，既保留了系统程序建模方式，同时又能提高蒸汽发生器换热计算的精度、稳定性和速度。
[0163]
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。