本发明涉及核反应堆严重事故领域，具体涉及一种核反应堆蒸汽爆炸综合分析方法。
背景技术：
日本福岛核电事故引发了对反应堆严重事故新的关注和研究，其中蒸汽爆炸因涉及高温高压下不同时间尺度、空间尺度现象的耦合问题而成为严重事故研究的热点与难点。堆内蒸汽爆炸会对压力容器下封头的早期失效及下封头能否有效冷却有很大影响。堆外蒸汽爆炸则可能造成安全壳失效，放射性物质外泄。因此研究蒸汽爆炸现象，对于熔融燃料与冷却剂接触反应的认识及评估现存压水反应堆严重事故进程中安全屏障完整性及放射性物质外泄与否有重要作用，同时对于新型反应堆的设计及安全评估也有重要意义。
目前国际上对于轻水堆严重事故的研究分析，基本都会考虑蒸汽爆炸的可能性及其影响。国际上已经开展了大量的蒸汽爆炸实验研究，这些实验研究主要分为两类，一类是采用控制变量法研究单一现象或机理的小规模实验，小规模实验主要用于研究fci的触发、碎裂和膜态沸腾等机理现象，为fci理论模型的开发提供支持，第二类实验在于研究和评估蒸汽爆炸后果的大规模实验研究这类实验主要研究触发蒸汽爆炸所需要的熔融物冷却剂粗混合条件、蒸汽爆炸的传播过程、蒸汽爆炸的载荷及能量转换率和尺度效应以及基于实验结果评估真实反应堆蒸汽爆炸的潜在后果。针对蒸汽爆炸也开展了大量理论研究，开发了多套蒸汽爆炸模型及程序，大部分蒸汽爆炸程序能分析粗混合过程及蒸汽爆炸过程，但没有对蒸汽爆炸发生后对反应堆压力容器和堆内结构产生的影响进行分析。为了确定其对压力容器以及堆腔结构完整性的影响，蒸汽爆炸的发生、压力波的传播及对压力容器以及堆腔结构的冲击影响均需要计算分析，从而为更加全面、有效地评估反应堆的安全性提供依据。
技术实现要素：
为了克服上述现有技术存在的问题，本发明的目的是为了综合分析蒸汽爆炸现象对反应堆的影响，提供一种综合分析方法。该方法先对蒸汽爆炸现象进行分析，计算粗混合阶段和爆炸膨胀阶段的熔融物碎裂情况，随后计算计算区域内的质量交换、动量交换和能量交换，通过压力迭代法对质量守恒方程、动量守恒方程和能量守恒方程计算，根据质量守恒方程的误差进行压力修正，根据修正后的压力进行下一步迭代，直至质量守恒方程的误差低于允许误差时结束迭代，并输出结果，主要包括蒸汽爆炸产生的压力波的幅值沿轴向的分布。在蒸汽爆炸计算完成后，采用商用有限元程序对反应堆压力容器、堆内结构和堆腔结构进行几何建模，并输入蒸汽爆炸计算得到的压力分布，计算得到反应堆压力容器、堆内结构和堆腔结构的应力变化，有助于综合评估反应堆的安全性。
为了实现上述目的，本发明采取了以下技术方案予以实施：
一种核反应堆蒸汽爆炸综合分析方法，步骤如下：
步骤1：获取核反应堆蒸汽爆炸初始参数，包括堆芯熔融物入射速度、堆芯熔融物入射直径、堆芯熔融物温度、堆芯熔融物密度、水池深度、水池温度和水池直径；
步骤2：计算核反应堆蒸汽爆炸过程熔融物的变化；计算中将蒸汽爆炸过程分为粗混合阶段和爆炸膨胀阶段，将堆芯熔融物简化成一连串的熔融物球状液滴，总共分为k组；
采用公式(1)计算粗混合过程熔融物直径随时间的变化，
式中：
we——韦伯数；
u——熔融物速度/m·s-1；
tn+1——第n+1时层；
tn——第n时层；
ρc——冷却剂密度/kg·m-3；
ρf——熔融物密度/kg·m-3；
采用公式(2)、公式(3)和公式(4)计算爆炸膨胀阶段中发生碎裂的熔融物质量，
式中：
mfr——发生碎裂的熔融物质量/kg；
c——碎裂常数；
vlim——空泡份额限制因子；
ρf——熔融物密度/kg·m-3；
rk,f——熔融物粒子的直径/m；
nk,p——第k组熔融物的粒子数；
mk,f——第k组熔融物的质量/kg；
vjet——冷却剂入射速度/m·s-1；
p——压力/pa；
pth——压力阈值/pa；
ρc——冷却剂密度/kg·m-3；
步骤3：计算相界面的能量交换，采用欧拉方法描述液相和蒸汽相，采用拉格朗日方法描述连续熔融物相和分散熔融物相；将相界面的能量交换分为连续液相、连续汽相内的换热；
连续液相内的能量交换主要考虑熔融物与汽-液相界面的对流换热，熔融物与汽-液相界面的辐射换热，主流液体与汽-液相界面及第k组熔融物的换热；连续汽相内的能量交换主要考虑第k组熔融物与液滴的辐射换热，熔融物与周围壁面的辐射换热；具体计算如下：
(1)连续液相内的换热模型
第k组熔融物与汽-液相界面的对流换热由公式(5)和公式(6)计算
hfilm＝max(hn,hf)公式(6)
式中：
qfilm——第k组熔融物与汽-液相界面的对流换热量/w；
dk,f——第k组熔融物直径/m；
hfilm——第k组熔融物表面蒸汽膜的对流换热系数/w·(m2·k)-1；
tf,s——第k组固体熔融物的温度/k；
tsat——液体的饱和温度/k；
hn——蒸汽膜自然对流换热系数/w·(m2·k)-1；
hf——蒸汽膜强制对流换热系数/w·(m2·k)-1；
第k组熔融物与汽-液相界面的辐射换热由公式(7)计算
式中：
qrad——第k组熔融物与汽-液相界面的辐射换热量/w；
dk,f——第k组熔融物直径/m；
hrad——第k组熔融物与汽-液相界面的辐射换热系数/w·(m2·k)-1；
tf——第k组熔融物的温度/k；
tsat——液体的饱和温度/k；
主流液体与汽-液相界面及第k组熔融物的换热由公式(8)计算
qliquid＝π(df+2δfilm)2hliquid(tsat-tl)+cradqrad公式(8)
式中：
qliquid——主流液体与汽-液相界面及第k组熔融物的换热量/w；
df——第k组熔融物直径/m；
δfilm——蒸汽膜厚度/m；
hliquid——主流液体与汽-液相界面的对流换热系数/w·(m2·k)-1；
tf——第k组熔融物的温度/k；
tsat——液体的饱和温度/k；
crad——辐射常数；
qrad——第k组熔融物与汽-液相界面的辐射换热量/w；
(2)连续汽相内的换热模型
第k组熔融物与液滴的辐射换热由公式(9)计算
式中：
tsat——液体的饱和温度；
第k组熔融物与周围壁面的辐射换热由公式(10)计算
式中：
tw——壁面温度；
步骤4：计算相界面的质量传递，采用欧拉方法描述液相和蒸汽相，采用拉格朗日方法描述连续熔融物相和分散熔融物相；将相界面质量传递分为熔融物周围蒸汽的蒸发，连续液相中的汽泡冷凝和连续汽相中的液滴蒸发；具体质量传递计算如下：
(1)熔融物周围蒸汽产生率由公式(11)和公式(12)计算
qnet,f＝qfilm+(1-crad)qrad-qliquid公式(12)
式中：
ms,f——蒸汽产生率/kg·s-1；
qnet,f——汽-液相界面上的净热流量/w；
hl,v——液体汽化潜热/j·kg-1；
vcell——控制体体积；
crad——辐射常数；
qfilm——熔融物与汽-液相界面的对流换热量/w；
qrad——熔融物与汽-液相界面的辐射换热量/w；
(2)连续液相中汽泡冷凝率由公式(13)和公式(14)计算
qnet,b＝ql,vs-ql,sl公式(14)
式中：
ms,b——汽泡冷凝率/kg·s-1；
qnet,b——主流液体内汽泡表面的的热流量/w；
hl,v——液体汽化潜热/j·kg-1；
ql,vs——汽泡向汽-液相界面传递的热量/w；
ql,sl——主流液体向汽-液相界面传递的热量/w；
(3)连续汽相中的液滴蒸发率由公式(15)和公式(16)计算
qnet,d＝(1-crad)qrad,l+qv,vl-qv,sl公式(16)
式中：
ms,d——液滴蒸发率/kg·s-1；
qnet,d——汽相内液滴表面的的热流量/w；
hl,v——液体汽化潜热/j·kg-1；
crad——辐射常数；
qrad,l——第k组熔融物与液滴间的辐射换热量/w；
qv,vl——连续蒸汽相向液滴的对流热量/w；
qv,sl——汽-液相界面向液滴的传热量/w；
单位体积内相变引起的质量传递由公式(17)计算，
ms＝ms,f+ms,b+ms,d公式(17)
式中：
ms——相变引起的质量传递率/kg·s-1；
ms,f——蒸汽产生率/kg·s-1；
ms,b——汽泡冷凝率/kg·s-1；
ms,d——液滴蒸发率/kg·s-1；
步骤5：计算相界面的动量交换，采用欧拉方法描述液相和蒸汽相，采用拉格朗日方法描述连续熔融物相和分散熔融物相；将相界面上的动量交换分为汽-液相间动量交换、汽-液相与熔融物的动量交换和与壁面的动量交换；
(1)汽-液相间动量交换主要由相间的曳力引起，两相间的曳力由稳态的相间曳力和虚拟质量力组成，由公式(18)、公式(19)和公式(20)计算，
fl＝fd+fv公式(18)
fd＝kv,l(uv-ul)公式(19)
式中：
fi——两相间的曳力/n·m-3；
fd——标准稳态相间曳力/n·m-3；；
fv——相间虚拟质量力/n·m-3；
kv,l——汽-液相间宏观阻力系数；
uv——液滴蒸发率/m·s-1；
ul——液滴蒸发率/m·s-1；
am——虚拟质量力系数；
(2)汽-液相冷却剂作用在堆芯熔融物上的拖曳力将影响堆芯熔融物的碎裂速率及下降的速度；汽-液相与熔融物的曳力由公式(21)计算，
式中：
uk,f——第k组熔融物的速度/m·s-1；
umix——汽-液冷却剂的速度/m·s-1；
(3)汽-液两相与壁面间的摩擦力由公式(22)和公式(23)计算
fv,w＝kv,wuv公式(22)
fl,w＝kl,wul公式(23)
式中：
fv,w——汽相与壁面间的摩擦力/n·m-3；
kv,w——汽相与壁面间的宏观阻力系数；
uv——汽相的速度/m·s-1；
fl,w——液相与壁面间的摩擦力/n·m-3；
kl,w——液相与壁面间的宏观阻力系数；
ul——液相的速度/m·s-1；
步骤6：计算中采用一维三流体数学物理模型建立质量守恒方程、动量、动量守恒方程和能量守恒方程，采用欧拉方法描述汽-液两相冷却剂，拉格朗日方法描述熔融物；采用交错网格对质量守恒方程、动量守恒方程和能量守恒方程进行离散，交错网格将质量、压力、温度、内能存储于质量能量控制体中心，速度存储于质量能量控制体边界，即动量控制体中心，相邻两节点的压力差构成了动量守恒方程中的压力梯度；采用压力迭代法对质量守恒方程、动量守恒方程和能量守恒方程进行求解；压力迭代循环过程中，先半隐式求解能量和动量方程得到汽相和液相的内能；
动量守恒方程的离散形式如公式(24)和公式(25)，
式中，
i——第i个质量能量控制体；
i+1/2——第i+1/2个动量控制体；
i+1——第i+1个质量能量控制体；
n——第n时层；
n+1——第n+1时层；
ρv——汽相的密度/kg·m-3；
uv——汽相的速度/m·s-1；
δuv——汽相的速度变化量/m·s-1；
δt——时间间隔/s；
δx——控制体间隔/m；
g——重力加速度/m·s-2；
αv——汽相加权空泡份额；
p——压力/pa；
kv,l——汽-液相间宏观阻力系数；
ul——液相的速度/m·s-1；
kv,w——汽相与壁面间的宏观阻力系数；
vv——汽相粘性耗散/n·m-1；
am——虚拟质量力系数；
γe——蒸汽产生率/kg·s-1；
mv,f——汽相与熔融物间阻力/n·m-1；
ρl——液相的密度/kg·m-3；
δul——液相的速度变化量/m·s-1；
αl——液相加权空泡份额；
kl,v——液汽相间宏观阻力系数；
kv,w——液相与壁面间的宏观阻力系数；
vls——液相粘性耗散/n·m-1；
γc——蒸汽凝固率/kg·s-1；
ml,f——液相与熔融物间阻力/n·m-1；
能量守恒方程的离散形式如公式(26)和公式(27)
式中，
i——第i个质量能量控制体；
n——第n时层；
n+1/2——第n+1/2时层；
n+1——第n+1时层；
δt——时间间隔/s；
δx——控制体间隔/m；
p——压力/pa；
ρv——汽相的密度/kg·m-3；
uv——汽相的速度/m·s-1；
iv——汽相比内能/j·kg-1；
αv——汽相加权空泡份额；
wv——单位体积内汽相粘性功/w·m-3；
qv,w——汽相与壁面间的换热量/w·m-3；
qv,f——汽相与熔融物间的换热量/w·m-3；
qv,int——汽相与汽-液相界面的换热量/w·m-3；
jv——汽相间的导热项/w·m-3；
sv——汽相间的内热源项/w·m-3；
γe——蒸汽产生率/kg·s-1；
γc——蒸汽凝固率/kg·s-1；
ρl——液相的密度/kg·m-3；
ul——液相的速度/m·s-1；
il——液相比内能/j·kg-1；
αl——液相加权空泡份额；
wl——单位体积内液相粘性功/w·m-3；
ql,w——液相与壁面间的换热量/w·m-3；
ql,f——液相与熔融物间的换热量/w·m-3；
ql,int——液相与汽-液相界面的换热量/w·m-3；
jl——液相间的导热项/w·m-3；
sl——液相间的内热源项/w·m-3；
步骤7：全隐式求解质量守恒方程，根据质量守恒方程的误差对蒸汽爆炸产生的压力进行修正，直至质量守恒方程的误差低于允许误差时结束迭代，并输出结果；输出结果主要包括蒸汽爆炸产生的压力波沿轴向的分布；
质量守恒方程的离散形式如公式(28)和公式(29)，
式中，
i——第i个质量能量控制体；
n——第n时层；
n+1——第n+1时层；
δt——时间间隔/s；
δx——控制体间隔/m；
ρv——汽相的密度/kg·m-3；
uv——汽相的速度/m·s-1；
ρl——液相的密度/kg·m-3；
ul——液相的速度/m·s-1；
γe——蒸汽产生率/kg·s-1；
γc——蒸汽凝固率/kg·s-1；
通过修正每个动量控制体内的压力，使连续方程的误差函数等于0，误差函数是其对应动量控制体内以及相邻动量控制体内压力的函数，第k+1次迭代的误差函数由公式(30)、公式(31)、公式(32)和公式(33)计算，
式中，
i——第i个质量能量控制体；
i-1/2——第i-1/2个动量控制体；
i-1——第i-1个质量能量控制体；
i+1/2——第i+1/2个动量控制体；
i+1——第i+1个质量能量控制体；
k——第k次迭代；
k+1——第k+1次迭代；
δt——时间间隔/s；
δx——控制体间隔/m；
ρ——密度/kg·m-3；
dik——第k次迭代第i个质量能量控制体的误差函数；
dik+1——第k+1次迭代第i个质量能量控制体的误差函数；
pi-1——第i-1个质量能量控制体的压力；
pi——第i个质量能量控制体的压力；
pi+1——第i+1个质量能量控制体的压力；
δpi-1——第i-1个质量能量控制体的压力修正值；
δpi——第i个质量能量控制体的压力修正值；
δpi+1——第i+1个质量能量控制体的压力修正值；
ρk——第k次迭代的密度/kg·m-3；
uk——第k次迭代的速度/m·s-1；
要使第k+1次迭代的误差函数为0，可通过公式(34)、公式(35)和公式(36)计算得到δp，
βk·δp＝-dk公式(34)
pk+1＝pk+δp公式(36)
式中，
i——第i个质量能量控制体；
i-1——第i-1个质量能量控制体；
i+1——第i+1个质量能量控制体；
k——第k次迭代；
k+1——第k+1次迭代；
d——质量能量控制体的误差函数；
p——质量能量控制体的压力；
βk——第k次迭代的系数矩阵；
pk+1——第k+1次迭代的压力；
pk——第k次迭代的压力；
δp——质量能量控制体的压力修正值；
步骤8：使用有限元程序abaqus在abaqus-cae界面下part模块对核反应堆压力容器、堆内结构和堆腔结构进行几何建模，并在abaqus-cae界面下使用assembly模块将上述在part模块中创建的核反应堆压力容器、堆内结构和堆腔结构装配起来；
步骤9：使用有限元程序abaqus在abaqus-cae界面下使用property模块创建混凝土、不锈钢和水的材料物性，创建截面属性、给部件赋予截面属性，在abaqus-cae界面下使用step模块和load模块将创建分析步骤，load模块中的载荷采用步骤7的蒸汽爆炸产生的压力波沿轴向的分布；
步骤10：使用有限元程序abaqus在abaqus-cae界面下使用mesh模块对步骤8中建立的反应堆压力容器、堆内结构和堆腔结构划分网格；
步骤11：使用有限元程序abaqus在abaqus-cae界面下使用job模块提交分析，输出反应堆压力容器、堆内结构和堆腔结构的应力变化。
和现有技术相比较，本发明具备如下优点：
本发明的综合分析方法能对核反应堆蒸汽爆炸事故下的压力波分布进行求解，并结合有限元程序分析该压力波对核反应堆压力容器、堆内结构和堆腔结构的冲击作用，从而为更加全面、有效地评估反应堆的安全性提供依据。
附图说明
图1是蒸汽爆炸综合分析方法的流程图。
图2是交错网格示意图。
图3是堆外蒸汽爆炸示意图。
图4是距离堆腔底部0.2m处的压力变化曲线。
图5是反应堆堆坑结构示意图。
具体实施方式
下面结合附图具体实施方式对本发明方法作进一步详细说明：
如图1所示，本发明一种核反应堆蒸汽爆炸综合分析方法，步骤如下：
步骤1：获取核反应堆蒸汽爆炸初始参数，包括堆芯熔融物入射速度、堆芯熔融物入射直径、堆芯熔融物温度、堆芯熔融物密度、水池深度、水池温度和水池直径；
步骤2：计算核反应堆蒸汽爆炸过程熔融物的变化；计算中将蒸汽爆炸过程分为粗混合阶段和爆炸膨胀阶段，将堆芯熔融物简化成一连串的熔融物球状液滴，总共分为k组；
采用公式(1)计算粗混合过程熔融物直径随时间的变化，
式中：
we——韦伯数；
u——熔融物速度/m·s-1；
tn+1——第n+1时层；
tn——第n时层；
ρc——冷却剂密度/kg·m-3；
ρf——熔融物密度/kg·m-3；
采用公式(2)、公式(3)和公式(4)计算爆炸膨胀阶段中发生碎裂的熔融物质量，
式中：
mfr——发生碎裂的熔融物质量/kg；
c——碎裂常数；
vlim——空泡份额限制因子；
ρf——熔融物密度/kg·m-3；
rk,f——熔融物粒子的直径/m；
nk,p——第k组熔融物的粒子数；
mk,f——第k组熔融物的质量/kg；
vjet——冷却剂入射速度/m·s-1；
p——压力/pa；
pth——压力阈值/pa；
ρc——冷却剂密度/kg·m-3；
步骤3：计算相界面的能量交换，采用欧拉方法描述液相和蒸汽相，采用拉格朗日方法描述连续熔融物相和分散熔融物相；将相界面的能量交换分为连续液相、连续汽相内的换热；
连续液相内的能量交换主要考虑熔融物与汽-液相界面的对流换热，熔融物与汽-液相界面的辐射换热，主流液体与汽-液相界面及第k组熔融物的换热；连续汽相内的能量交换主要考虑第k组熔融物与液滴的辐射换热，熔融物与周围壁面的辐射换热；具体计算如下：
(1)连续液相内的换热模型
第k组熔融物与汽-液相界面的对流换热由公式(5)和公式(6)计算
hfilm＝max(hn,hf)公式(6)
式中：
qfilm——第k组熔融物与汽-液相界面的对流换热量/w；
dk,f——第k组熔融物直径/m；
hfilm——第k组熔融物表面蒸汽膜的对流换热系数/w·(m2·k)-1；
tf,s——第k组固体熔融物的温度/k；
tsat——液体的饱和温度/k；
hn——蒸汽膜自然对流换热系数/w·(m2·k)-1；
hf——蒸汽膜强制对流换热系数/w·(m2·k)-1；
第k组熔融物与汽-液相界面的辐射换热由公式(7)计算
式中：
qrad——第k组熔融物与汽-液相界面的辐射换热量/w；
dk,f——第k组熔融物直径/m；
hrad——第k组熔融物与汽-液相界面的辐射换热系数/w·(m2·k)-1；
tf——第k组熔融物的温度/k；
tsat——液体的饱和温度/k；
主流液体与汽-液相界面及第k组熔融物的换热由公式(8)计算
qliquid＝π(df+2δfilm)2hliquid(tsat-tl)+cradqrad公式(8)式中：
qliquid——主流液体与汽-液相界面及第k组熔融物的换热量/w；
df——第k组熔融物直径/m；
δfilm——蒸汽膜厚度/m；
hliquid——主流液体与汽-液相界面的对流换热系数/w·(m2·k)-1；
tf——第k组熔融物的温度/k；
tsat——液体的饱和温度/k；
crad——辐射常数；
qrad——第k组熔融物与汽-液相界面的辐射换热量/w；
(2)连续汽相内的换热模型
第k组熔融物与液滴的辐射换热由公式(9)计算
式中：
tsat——液体的饱和温度；
第k组熔融物与周围壁面的辐射换热由公式(10)计算
式中：
tw——壁面温度；
步骤4：计算相界面的质量传递，采用欧拉方法描述液相和蒸汽相，采用拉格朗日方法描述连续熔融物相和分散熔融物相；将相界面质量传递分为熔融物周围蒸汽的蒸发，连续液相中的汽泡冷凝和连续汽相中的液滴蒸发；具体质量传递计算如下：
(1)熔融物周围蒸汽产生率由公式(11)和公式(12)计算
qnet,f＝qfilm+(1-crad)qrad-qliquid公式(12)
式中：
ms,f——蒸汽产生率/kg·s-1；
qnet,f——汽-液相界面上的净热流量/w；
hl,v——液体汽化潜热/j·kg-1；
vcell——控制体体积；
crad——辐射常数；
qfilm——熔融物与汽-液相界面的对流换热量/w；
qrad——熔融物与汽-液相界面的辐射换热量/w；
(2)连续液相中汽泡冷凝率由公式(13)和公式(14)计算
qnet,b＝ql,vs-ql,sl公式(14)
式中：
ms,b——汽泡冷凝率/kg·s-1；
qnet,b——主流液体内汽泡表面的的热流量/w；
hl,v——液体汽化潜热/j·kg-1；
ql,vs——汽泡向汽-液相界面传递的热量/w；
ql,sl——主流液体向汽-液相界面传递的热量/w；
(3)连续汽相中的液滴蒸发率由公式(15)和公式(16)计算
qnet,d＝(1-crad)qrad,l+qv,vl-qv,sl公式(16)
式中：
ms,d——液滴蒸发率/kg·s-1；
qnet,d——汽相内液滴表面的的热流量/w；
hl,v——液体汽化潜热/j·kg-1；
crad——辐射常数；
qrad,l——第k组熔融物与液滴间的辐射换热量/w；
qv,vl——连续蒸汽相向液滴的对流热量/w；
qv,sl——汽-液相界面向液滴的传热量/w；
单位体积内相变引起的质量传递由公式(17)计算，
ms＝ms,f+ms,b+ms,d公式(17)
式中：
ms——相变引起的质量传递率/kg·s-1；
ms,f——蒸汽产生率/kg·s-1；
ms,b——汽泡冷凝率/kg·s-1；
ms,d——液滴蒸发率/kg·s-1；
步骤5：计算相界面的动量交换，采用欧拉方法描述液相和蒸汽相，采用拉格朗日方法描述连续熔融物相和分散熔融物相；将相界面上的动量交换分为汽-液相间动量交换、汽-液相与熔融物的动量交换和与壁面的动量交换；
(1)汽-液相间动量交换主要由相间的曳力引起，两相间的曳力由稳态的相间曳力和虚拟质量力组成，由公式(18)、公式(19)和公式(20)计算，
fl＝fd+fv公式(18)
fd＝kv,l(uv-ul)公式(19)
式中：
fi——两相间的曳力/n·m-3；
fd——标准稳态相间曳力/n·m-3；；
fv——相间虚拟质量力/n·m-3；
kv,l——汽-液相间宏观阻力系数；
uv——液滴蒸发率/m·s-1；
ul——液滴蒸发率/m·s-1；
am——虚拟质量力系数；
(2)汽-液相冷却剂作用在堆芯熔融物上的拖曳力将影响堆芯熔融物的碎裂速率及下降的速度；汽-液相与熔融物的曳力由公式(21)计算，
式中：
uk,f——第k组熔融物的速度/m·s-1；
umix——汽-液冷却剂的速度/m·s-1；
(3)汽-液两相与壁面间的摩擦力由公式(22)和公式(23)计算
fv,w＝kv,wuv公式(22)
fl,w＝kl,wul公式(23)
式中：
fv,w——汽相与壁面间的摩擦力/n·m-3；
kv,w——汽相与壁面间的宏观阻力系数；
uv——汽相的速度/m·s-1；
fl,w——液相与壁面间的摩擦力/n·m-3；
kl,w——液相与壁面间的宏观阻力系数；
ul——液相的速度/m·s-1；
步骤6：计算中采用一维三流体数学物理模型建立质量守恒方程、动量、动量守恒方程和能量守恒方程，采用欧拉方法描述汽-液两相冷却剂，拉格朗日方法描述熔融物；采用交错网格对质量守恒方程、动量守恒方程和能量守恒方程进行离散，如图2所示，交错网格将质量、压力、温度、内能存储于质量能量控制体中心，速度存储于质量能量控制体边界，即动量控制体中心，相邻两节点的压力差构成了动量守恒方程中的压力梯度；采用压力迭代法对质量守恒方程、动量守恒方程和能量守恒方程进行求解；压力迭代循环过程中，先半隐式求解能量和动量方程得到汽相和液相的内能；
动量守恒方程的离散形式如公式(24)和公式(25)，
式中，
i——第i个质量能量控制体；
i+1/2——第i+1/2个动量控制体；
i+1——第i+1个质量能量控制体；
n——第n时层；
n+1——第n+1时层；
ρv——汽相的密度/kg·m-3；
uv——汽相的速度/m·s-1；
δuv——汽相的速度变化量/m·s-1；
δt——时间间隔/s；
δx——控制体间隔/m；
g——重力加速度/m·s-2；
αv——汽相加权空泡份额；
p——压力/pa；
kv,l——汽-液相间宏观阻力系数；
ul——液相的速度/m·s-1；
kv,w——汽相与壁面间的宏观阻力系数；
vv——汽相粘性耗散/n·m-1；
am——虚拟质量力系数；
γe——蒸汽产生率/kg·s-1；
mv,f——汽相与熔融物间阻力/n·m-1；
ρl——液相的密度/kg·m-3；
δul——液相的速度变化量/m·s-1；
αl——液相加权空泡份额；
kl,v——液汽相间宏观阻力系数；
kv,w——液相与壁面间的宏观阻力系数；
vls——液相粘性耗散/n·m-1；
γc——蒸汽凝固率/kg·s-1；
ml,f——液相与熔融物间阻力/n·m-1；
能量守恒方程的离散形式如公式(26)和公式(27)
式中，
i——第i个质量能量控制体；
n——第n时层；
n+1/2——第n+1/2时层；
n+1——第n+1时层；
δt——时间间隔/s；
δx——控制体间隔/m；
p——压力/pa；
ρv——汽相的密度/kg·m-3；
uv——汽相的速度/m·s-1；
iv——汽相比内能/j·kg-1；
αv——汽相加权空泡份额；
wv——单位体积内汽相粘性功/w·m-3；
qv,w——汽相与壁面间的换热量/w·m-3；
qv,f——汽相与熔融物间的换热量/w·m-3；
qv,int——汽相与汽-液相界面的换热量/w·m-3；
jv——汽相间的导热项/w·m-3；
sv——汽相间的内热源项/w·m-3；
γe——蒸汽产生率/kg·s-1；
γc——蒸汽凝固率/kg·s-1；
ρl——液相的密度/kg·m-3；
ul——液相的速度/m·s-1；
il——液相比内能/j·kg-1；
αl——液相加权空泡份额；
wl——单位体积内液相粘性功/w·m-3；
ql,w——液相与壁面间的换热量/w·m-3；
ql,f——液相与熔融物间的换热量/w·m-3；
ql,int——液相与汽-液相界面的换热量/w·m-3；
jl——液相间的导热项/w·m-3；
sl——液相间的内热源项/w·m-3；
步骤7：全隐式求解质量守恒方程，根据质量守恒方程的误差对蒸汽爆炸产生的压力进行修正，直至质量守恒方程的误差低于允许误差时结束迭代，并输出结果；输出结果主要包括蒸汽爆炸产生的压力波沿轴向的分布；
质量守恒方程的离散形式如公式(28)和公式(29)，
式中，
i——第i个质量能量控制体；
n——第n时层；
n+1——第n+1时层；
δt——时间间隔/s；
δx——控制体间隔/m；
ρv——汽相的密度/kg·m-3；
uv——汽相的速度/m·s-1；
ρl——液相的密度/kg·m-3；
ul——液相的速度/m·s-1；
γe——蒸汽产生率/kg·s-1；
γc——蒸汽凝固率/kg·s-1；
通过修正每个动量控制体内的压力，使连续方程的误差函数等于0，误差函数是其对应动量控制体内以及相邻动量控制体内压力的函数，第k+1次迭代的误差函数由公式(30)、公式(31)、公式(32)和公式(33)计算，
式中，
i——第i个质量能量控制体；
i-1/2——第i-1/2个动量控制体；
i-1——第i-1个质量能量控制体；
i+1/2——第i+1/2个动量控制体；
i+1——第i+1个质量能量控制体；
k——第k次迭代；
k+1——第k+1次迭代；
δt——时间间隔/s；
δx——控制体间隔/m；
ρ——密度/kg·m-3；
pi-1——第i-1个质量能量控制体的压力；
pi——第i个质量能量控制体的压力；
pi+1——第i+1个质量能量控制体的压力；
δpi-1——第i-1个质量能量控制体的压力修正值；
δpi——第i个质量能量控制体的压力修正值；
δpi+1——第i+1个质量能量控制体的压力修正值；
ρk——第k次迭代的密度/kg·m-3；
uk——第k次迭代的速度/m·s-1；
要使第k+1次迭代的误差函数为0，可通过公式(34)、公式(35)和公式(36)计算得到δp，
βk·δp＝-dk公式(34)
pk+1＝pk+δp公式(36)
式中，
i——第i个质量能量控制体；
i-1——第i-1个质量能量控制体；
i+1——第i+1个质量能量控制体；
k——第k次迭代；
k+1——第k+1次迭代；
d——质量能量控制体的误差函数；
p——质量能量控制体的压力；
βk——第k次迭代的系数矩阵；
pk+1——第k+1次迭代的压力；
pk——第k次迭代的压力；
δp——质量能量控制体的压力修正值；
步骤8：使用有限元程序abaqus在abaqus-cae界面下part模块对核反应堆压力容器、堆内结构和堆腔结构进行几何建模，并在abaqus-cae界面下使用assembly模块将上述在part模块中创建的核反应堆压力容器、堆内结构和堆腔结构装配起来；
步骤9：使用有限元程序abaqus在abaqus-cae界面下使用property模块创建混凝土、不锈钢和水的材料物性，创建截面属性、给部件赋予截面属性，在abaqus-cae界面下使用step模块和load模块将创建分析步骤，load模块中的载荷采用步骤7的蒸汽爆炸产生的压力波沿轴向的分布；
步骤10：使用有限元程序abaqus在abaqus-cae界面下使用mesh模块对步骤8中建立的反应堆压力容器、堆内结构和堆腔结构划分网格；
步骤11：使用有限元程序abaqus在abaqus-cae界面下使用job模块提交分析，输出反应堆压力容器、堆内结构和堆腔结构的应力变化。
下面结合具体的计算对象说明本发明的效果，以图3所示的堆外蒸汽爆炸为例，先获取核反应堆蒸汽爆炸初始参数，包括堆芯熔融物入射速度、堆芯熔融物入射直径、堆芯熔融物温度、堆芯熔融物密度、水池深度、水池温度和水池直径，堆芯熔融物落入水中后会发生蒸汽爆炸，根据上述步骤2-步骤7可以对这一过程进行计算，包括熔融物在粗混合阶段和爆炸膨胀阶段的碎裂情况，相界面的交换，最后输出蒸汽爆炸的压力波沿轴向的分布，距离底部0.2m处的压力变化曲线如图4所示。根据步骤8采用abaqus对图5所示的堆腔结构进行建模，根据步骤9采用abaqus设置混凝土和水的物性，创建截面属性、给部件赋予截面属性，同时创建分析步骤，采用的载荷为图4的压力波曲线，根据步骤10对图5所示的堆腔结构进行网格划分，最后提交分析，输出结果。输出的结果主要为堆腔结构的应力变化曲线，从中可以得到堆腔结构应力最大位置变化，有助于对堆腔结构的完整性进行分析，从而为更加全面、有效地评估反应堆的安全性提供依据。