CN111125972B - 核电厂破口失水事故水力载荷分析方法 - Google Patents

Abstract

一种核电厂破口失水事故水力载荷分析方法，包括：考虑了控制体面积变化随时间变化的两流体水力学模型；基于半隐差分的数值解法；基于逐步分段解析法的结构动力学分析方法；实现了水力学与力学之间的流固耦合计算。与现有技术相比，本发明方法考虑了控制体面积随时间的变化，实现了水力学与力学的耦合计算，能够更加精确地进行核电厂失水事故下结构受力分析，为反应堆的设计提供标准。

Description

核电厂破口失水事故水力载荷分析方法

技术领域

本发明属于核电厂安全分析技术领域，具体涉及一种核电厂失水 事故水力载荷分析方法。

背景技术

作为压水堆核电厂设计基准事故之一，冷却剂丧失事故一直是国 内外核反应堆设计时关注的重点。破口失水事故发生时，瞬间产生的 卸压载荷会引起反应堆堆内结构的振动，威胁反应堆压力容器的安全 性和完整性，所以反应堆的破口失水事故动力分析一直受到核电发达 国家的高度重视，开发了一系列的分析方法进行数值模拟。但水力载 荷分析方法在我国依旧处于发展阶段，设计先进的破口失水事故水力 载荷分析方法对于我国的核电自主化具有重要意义。

目前国际主流的核反应堆水力载荷分析方法已经得到了广泛的 应用和验证，其发展已相对成熟。但是依旧存在一些问题亟待解决和 优化。

对于常见的水力载荷专用分析程序，如：MULTIFLUX或者ATHIS 程序，其使用均相流模型模拟两相流动，而均相流模型将两相流体均 匀处理，没有考虑相间的相互作用，模拟两相流动存在一定的误差；

当使用反应堆事故分析程序分析水力载荷时，由于程序全都认为 管道时刚性的，流体流动面积不会随时间变化，但是破口失水事故发 生后，瞬间产生的泄压波会冲击反应堆压力容器内部结构，导致结构 震动，流道发生变化，因此这些程序计算水力载荷时存在偏差；

一些计算流体力学程序(如：ANSYS或者ABAQUS)也可以用来 分析反应堆的破口失水事故水力载荷，但是这些程序建模复杂，计算 耗费大。

总之，目前已存在的水力载荷分析方法存在一些模型方面的问题， 会导致核电厂水力载荷分析的结果不够准确和真实，因此需要针对模 型和算法进行更适用于破口失水事故水力载荷分析的研究。

发明内容

本发明中开展了相关研究以解决现有水力载荷分析方法存在的 问题。首先建立了考虑控制体面积随时间变化的两流体水力学分析模 型，计算一回路各个节点的水力学参数；其次建立了基于模态分析的 结构力学分析方法，使用逐步分段解析法求解结构动力学方程；随后 建立了流固耦合模型，实现了水力学与力学的耦合计算；最后基于计 算结果计算结构受力，得到了核电厂破口失水事故水力载荷分析方法。

本发明采用如下技术方案：

一种核电厂破口失水事故水力载荷分析方法，其特征在于，包括 如下步骤：

第一步：根据核电厂一回路的运行参数和管道结构参数建立数学 模型；

第二步：确定反应堆吊篮的节点划分与周围管道的对应关系，并 对吊篮结构进行模态振型分析，得到吊篮结构的固有振动频率和模态 振型向量参数；

由于吊篮结构复杂，因此将圆筒状的模型当作板状进行简化处理， 并划分成若干个小的振荡单元，而每个小的振荡单元均存在结构动力 学方程：

式中：m—吊篮结构的质量；c—吊篮结构的阻尼系数；k—吊篮结 构的刚度系数；v(t)—吊篮结构的位移；f(t)—吊篮结构随时间的受 到的载荷；

——吊篮结构的速度；

——吊篮结构的加速度；

由于各个吊篮结构振荡单元之间存在相互影响，所以m、k、c都 不是对角矩阵，为了消除这种相互影响，简化计算，引入模态振型分 析；首先，通过坐标变换，将物理坐标转化为模态坐标：

v(t)＝φ_ny_n(t) (2)

式中：φ_n—吊篮结构的模态矩阵，每一列各代表一种模态的振型；y_n(t)—第n种模态的模态位移；

对于吊篮结构无阻尼自由振动方程组—式(3)，存在振荡单元振 动的位移解析解—式(4),其中w—吊篮的固有振动频率，φ—吊篮的 振幅,即模态振型矩阵，α—吊篮的振动初始角频率；

v(t)＝φsin(w·t+α) (4)

将式(4)代入式(3)中，得：

M^-1Kφ＝w²φ (5)

式中：M^-1—吊篮质量矩阵的逆矩阵，是n×n阶矩阵，K—吊篮模态刚 度矩阵，是n×n阶矩阵，φ—吊篮的模态振型矩阵，是n×h阶矩阵， 其中第l行第j列的元素表示第l个振荡单元第j种模态下的振幅，h 为预设的模态数，n为振荡单元的个数；

式(5)是一个特征值问题，能够求解得出每一阶模态对应的特征值 和特征向量，特征值是吊篮结构的固有振动频率w²，特征向量是该频 率下的模态振型向量φ_i，其中φ_i是第i个振动频率对应的模态向量， 是n×1阶矩阵；

每个小的振荡单元的结构动力学方程—式(1)组成的方程组转变 为一组独立的、以模态坐标和模态参数描述的方程—式(6)，从而实 现方程组的独立解耦；

式中：M—吊篮模态质量矩阵，

—吊篮第n种模态的模态加速度， C—吊篮模态阻尼矩阵，

—吊篮第n种模态的模态速度，K—吊 篮模态刚度矩阵，F(t)—吊篮模态受力矩阵；

第三步：选用Henry-Fauske模型计算破口处的临界流量，根据 流体流动状态判断流型，选用相应流型对应的换热关系式，并计算流 体受到的摩擦力；

第四步：使用半隐差分算法求解水力学守恒方程；水力学模型是 两流体模型，在水力学守恒方程中考虑了管道的面积变化随时间的变 化,如式(7)的汽相质量守恒方程；

式中：A—控制体截面积，α_g—控制体中流体的汽相空泡份额，ρ_g— 流体的汽相密度；V_g—流体的汽相速度；Γ_g—控制体中产生的汽相的 质量；

水力学模块采用交错网格，在网格中心采用标量网格，网格边界 采用矢量网格；对于水力学守恒方程的求解，首先求解两相的动量守 恒方程，得到新时刻控制体速度和压力的关系式—式(8)；

式中：

—接管j处i相新时刻的速度，

—接管j处i相旧时刻的 速度，

—接管j上游的控制体k的压力，P_l ⁿ⁺¹—接管j下游的控制 体k的压力，A_i—动量方程中接管j处i相旧时刻速度的系数，B_i— 动量方程中接管j上游的控制体k压力的系数，C_i—接管j下游的控 制体k压力的系数，D_i—i相动量方程中除旧时刻速度项和上下游压 力项外的其他项；其中i相为汽相或液相；

其次求解质量和能量守恒方程，将汽相内能、液相内能、汽相空 泡份额和压力四个变量作为求解变量，将汽相质量守恒方程、液相质 量守恒方程、汽相能量守恒方程、液相能量守恒方程四个守恒方程的 方程组，利用高斯消去法求解方程组，先得到控制体压力和速度的关 系式—式(9)；

式中：P_i ⁿ⁺¹—控制体i新时刻的压力，P_i ⁿ—控制体i旧时刻的压力，

—控制体i进口的汽相新时刻的速度，

—控制体i出口的汽相新 时刻的速度，

—控制体i进口的液相新时刻的速度，

—控 制体i出口的液相新时刻的速度，a、b、c、d分别表示方程组中控制 体i液相出、进口速度的系数和方程组中控制体i汽相出、进口速度 的系数，e—方程组中除控制体i气液相进出口速度以外的其他项；

然后联立压力和速度的关系式(8)和式(9)，求解得到新时刻的压 力和速度；最后通过质量和能量守恒方程方程组，回代新时刻的压力 求得气液相内能、汽相空泡份额变量；

第五步：利用水力学守恒方程中求得的控制体的气、液相速度、 空泡份额和密度计算每个控制体的流体总质量，总质量除以控制体的 体积即该控制体中的流体混合密度ρ_m,1；同时根据控制体的压力和水 的物性表查找水和水蒸汽的密度，然后根据式(10)计算控制体中的流 体混合密度ρ_m,2；

ρ_m,2＝α_gρ_g+α_fρ_f (10)

第六步：判断流体速度和物性参数即压力、温度、内能、密度的 计算是否收敛，若不收敛则时间步长减半并重新进行第三步计算，直 至满足收敛准则；

收敛准则为：1、比较利用水力学守恒方程计算得到的控制体中 的流体混合密度ρ_m,1和状态方程求得的控制体中的流体混合密度ρ_m,2， 二者的相对误差需要满足预设的限值；2、由于水力学守恒方程求解 采用半隐算法，因此所有的控制体需要满足声速库伦特限值条件；

第七步：判断外迭代是否收敛；如果不收敛，则进行第八步结构 受力计算和第九步位移计算，再通过第九十步计算管道新的面积和面 积变化的相对变化率，然后重复第四步、第五步和第六步；如果收敛， 则进行第十一步计算反应堆堆内结构部件的受力计算；收敛准则为前 后两次外迭代计算得到的管道面积变化的相对变化率小于规定的限 值；因为需要比较前后两次迭代的面积相对变化值，因此在进行每一 个时间步长内的首次外迭代计算时，由于没有前一次迭代的计算值， 无法进行外迭代收敛判断，因此外迭代默认不收敛，即每个时间步内 首次进行外迭代时不进行第七步；

第八步：传递水力学的压力到力学计算部分，求解吊篮结构受到 的载荷；将圆筒状的吊篮当作板状，将吊篮周围的控制体的压力投影 到垂直于板的方向，将所有的投影值相加即得到吊篮结构的受力；

第九步：利用逐步分段解析法求解吊篮结构的位移；

将时间划分成若干个不大于10^-5秒的时间步长，并且假设一个时 间步长内结构受到的力是不变的，由于结构动力学方程是二阶常微分 方程，因此存在解析解：

式中：w_D—阻尼圆周频率，w—无阻尼圆周频率，ξ——阻尼比，P₀— 结构受到的力，k—刚度系数，y(τ)—吊篮结构的模态位移；

由于阻尼系数对吊篮结构的影响小，因此忽略阻尼的影响，即阻 尼比为0，阻尼圆周频率等于无阻尼圆周频率即自由振动频率，则解 析解变为：

将式(10)求一阶导数，得到模态速度随时间变化的关系式：

根据初始条件τ＝0和刚度系数k：

k＝mw² (16)

由式(12)、式(13)、式(14)、式(15)和式(16)联立求解，得到模 态位移随时间的变化；

将当前计算时间步吊篮结构中每一个振荡单元的受力初值p₀、速 度初值

和位移初值y₀代入式(17)得到该时间步的模态位移，通过式 (2)将模态位移转化为物理坐标下的位移，即得到吊篮各部分的位移；

第十步：将吊篮结构的位移转化为管道的新面积，传递给水力学 计算部分；由于水力学守恒方程中求解需要面积变化项dA/A，所以需 要将吊篮结构的位移转化为面积；转化方法是通过式(18)将投影到一 维方向上的位移重新转化到真实的方向上，则吊篮结构相邻管道的直 径也会变化响应的值，然后计算管道的新面积；

式中：A_i—管道i的横截面积；s_i—管道i对应振荡单元的面积；l_i— 管道i的长度；x_j—振荡单元j的位移；θ_j—振荡单元j的方位角；

第十一步：计算反应堆压力容器及其内部结构和一回路管道的受 力，并分析结构的受力是否满足安全要求；反应堆堆内结构的受力分 为水平方向的力和竖直方向的力，水平方向的力是将所有反应堆内吊 篮结构水平方向上相邻管道的压力投影到x和y方向上然后计算合力， 竖直方向上的力是将反应堆内吊篮结构受到的压力和流体施加的摩 擦力和局部阻力相加；管道受到的力是将管道分为很多段，利用动量 方程—式19计算每一段管道的受力；

式中：F—管道段收到的力，

—管道段i当前时间步的流体质量流 密度，

—管道段i上一时间步的流体质量流密度，A—管道截面积， Δx—管道段长度，Δt—时间步长，

—管道段进口截面流体压力，

—管道段出口截面流体压力，

—管道段进口截面流体质量流 密度，

—管道段出口截面流体质量流密度，

——管道段进口 截面的流速，

——管道段出口截面的流速；

第十二步：判断当前计算时间是否到达结束时间；若没有结束， 则返回到第三步，重新计算；若结束，则结束计算。

和现有技术相比较，本发明具备如下优点：

1：由于常见的核电厂事故分析程序没有考虑控制体面积随时间 的变化，导致计算水力载荷时结果不可靠。本发明中使用了考虑控制 体面积变化的水力学模型，该模型可以模拟由于瞬间的泄压波导致的 压力容器内部流道面积发生变化的现象，并通过结构动力学分析计算 结构位移，修正管道面积，结构的水力载荷计算更加准确。

2：由于传统的专用水力载荷分析方法使用均相流模型，导致计 算两相流的计算结果不准确。本发明中使用两流体模型模拟两相流动 现象，该模型分别对气相和液相求解守恒方程，且考虑了两相之间的 质量、动量和能量交换，可以较为真实的反映两相流各种物理现象内 在机理的实际过程；

3：计算流体力学程序建模复杂，计算耗时耗力，效率较低。本 发明中建模简单，无需划分大量的网格也可以精确的计算，计算量小， 计算速度快。

附图说明

图1为本发明提出的适用于核电厂破口失水事故水力载荷分析 方法的流程示意图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

如图1所示，本发明一种核电厂破口失水事故水力载荷分析方法， 包括如下步骤：

第一步：根据核电厂一回路的运行参数和管道结构参数建立数学 模型；

第二步：确定反应堆吊篮的节点划分与周围管道的对应关系，并 对吊篮结构进行模态振型分析，得到吊篮结构的固有振动频率和模态 振型向量参数；

由于吊篮结构复杂，因此将圆筒状的模型当作板状进行简化处理， 并划分成若干个小的振荡单元，而每个小的振荡单元均存在结构动力 学方程：

式中：m—吊篮结构的质量；c—吊篮结构的阻尼系数；k—吊篮结 构的刚度系数；v(t)—吊篮结构的位移；f(t)—吊篮结构随时间的受 到的载荷；

——吊篮结构的速度；

——吊篮结构的加速度；

由于各个吊篮结构振荡单元之间存在相互影响，所以m、k、c都 不是对角矩阵，为了消除这种相互影响，简化计算，引入模态振型分 析；首先，通过坐标变换，将物理坐标转化为模态坐标：

v(t)＝φ_ny_n(t) (2)

式中：φ_n—吊篮结构的模态矩阵，每一列各代表一种模态的振型； y_n(t)—第n种模态的模态位移；

对于吊篮结构无阻尼自由振动方程组—式(3)，存在振荡单元振 动的位移解析解—式(4),其中w—吊篮的固有振动频率，φ—吊篮的 振幅,即模态振型矩阵，α—吊篮的振动初始角频率；

v(t)＝φsin(w·t+α) (4)

将式(4)代入式(3)中，得：

M^-1Kφ＝w²φ (5)

式中：M^-1—吊篮质量矩阵的逆矩阵，是n×n阶矩阵，K—吊篮模态刚 度矩阵，是n×n阶矩阵，φ—吊篮的模态振型矩阵，是n×h阶矩阵， 其中第l行第j列的元素表示第l个振荡单元第j种模态下的振幅，h 为预设的模态数，n为振荡单元的个数；

式(5)是一个特征值问题，能够求解得出每一阶模态对应的特征值 和特征向量，特征值是吊篮结构的固有振动频率w²，特征向量是该频 率下的模态振型向量φ_i，其中φ_i是第i个振动频率对应的模态向量， 是n×1阶矩阵；

每个小的振荡单元的结构动力学方程—式(1)组成的方程组转变 为一组独立的、以模态坐标和模态参数描述的方程—式(6)，从而实 现方程组的独立解耦；

式中：M—吊篮模态质量矩阵，

—吊篮第n种模态的模态加速度， C—吊篮模态阻尼矩阵，

—吊篮第n种模态的模态速度，K—吊 篮模态刚度矩阵，F(t)—吊篮模态受力矩阵；

第三步：选用Henry-Fauske模型计算破口处的临界流量，根据 流体流动状态判断流型，选用相应流型对应的换热关系式，并计算流 体受到的摩擦力；

第四步：使用半隐差分算法求解水力学守恒方程；水力学模型是 两流体模型，在水力学守恒方程中考虑了管道的面积变化随时间的变 化,如式(7)的汽相质量守恒方程；

式中：A—控制体截面积，α_g—控制体中流体的汽相空泡份额，ρ_g— 流体的汽相密度；V_g—流体的汽相速度；Γ_g—控制体中产生的汽相的 质量；

水力学模块采用交错网格，在网格中心采用标量网格，网格边界 采用矢量网格；对于水力学守恒方程的求解，首先求解两相的动量守 恒方程，得到新时刻控制体速度和压力的关系式—式(8)；

式中：

—接管j处i相新时刻的速度，

—接管j处i相旧时刻的 速度，

—接管j上游的控制体k的压力，P_l ⁿ⁺¹—接管j下游的控制 体k的压力，A_i—动量方程中接管j处i相旧时刻速度的系数，B_i— 动量方程中接管j上游的控制体k压力的系数，C_i—接管j下游的控 制体k压力的系数，D_i—i相动量方程中除旧时刻速度项和上下游压 力项外的其他项；其中i相为汽相或液相；

其次求解质量和能量守恒方程，将汽相内能、液相内能、汽相空 泡份额和压力四个变量作为求解变量，将汽相质量守恒方程、液相质 量守恒方程、汽相能量守恒方程、液相能量守恒方程四个守恒方程的 方程组，利用高斯消去法求解方程组，先得到控制体压力和速度的关 系式—式(9)；

式中：P_i ⁿ⁺¹—控制体i新时刻的压力，P_i ⁿ—控制体i旧时刻的压力，

—控制体i进口的汽相新时刻的速度，

—控制体i出口的汽相新 时刻的速度，

—控制体i进口的液相新时刻的速度，

—控 制体i出口的液相新时刻的速度，a、b、c、d分别表示方程组中控制 体i液相出、进口速度的系数和方程组中控制体i汽相出、进口速度 的系数，e—方程组中除控制体i气液相进出口速度以外的其他项；

然后联立压力和速度的关系式(8)和式(9)，求解得到新时刻的压 力和速度；最后通过质量和能量守恒方程方程组，回代新时刻的压力 求得气液相内能、汽相空泡份额变量；

第五步：利用水力学守恒方程中求得的控制体的气、液相速度、 空泡份额和密度计算每个控制体的流体总质量，总质量除以控制体的 体积即该控制体中的流体混合密度ρ_m,1；同时根据控制体的压力和水 的物性表查找水和水蒸汽的密度，然后根据式(10)计算控制体中的流 体混合密度ρ_m,2；

ρ_m,2＝α_gρ_g+α_fρ_f (10)

第六步：判断流体速度和物性参数即压力、温度、内能、密度的 计算是否收敛，若不收敛则时间步长减半并重新进行第三步计算，直 至满足收敛准则；

收敛准则为：1、比较利用水力学守恒方程计算得到的控制体中 的流体混合密度ρ_m,1和状态方程求得的控制体中的流体混合密度ρ_m,2， 二者的相对误差需要满足预设的限值；2、由于水力学守恒方程求解 采用半隐算法，因此所有的控制体需要满足声速库伦特限值条件；

第七步：判断外迭代是否收敛；如果不收敛，则进行第八步结构 受力计算和第九步位移计算，再通过第九十步计算管道新的面积和面 积变化的相对变化率，然后重复第四步、第五步和第六步；如果收敛， 则进行第十一步计算反应堆堆内结构部件的受力计算；收敛准则为前 后两次外迭代计算得到的管道面积变化的相对变化率小于规定的限 值；因为需要比较前后两次迭代的面积相对变化值，因此在进行每一 个时间步长内的首次外迭代计算时，由于没有前一次迭代的计算值， 无法进行外迭代收敛判断，因此外迭代默认不收敛，即每个时间步内 首次进行外迭代时不进行第七步；

第八步：传递水力学的压力到力学计算部分，求解吊篮结构受到 的载荷；将圆筒状的吊篮当作板状，将吊篮周围的控制体的压力投影 到垂直于板的方向，将所有的投影值相加即得到吊篮结构的受力；

第九步：利用逐步分段解析法求解吊篮结构的位移；

将时间划分成若干个不大于10^-5秒的时间步长，并且假设一个时 间步长内结构受到的力是不变的，由于结构动力学方程是二阶常微分 方程，因此存在解析解：

式中：w_D—阻尼圆周频率，w—无阻尼圆周频率，ξ——阻尼比，P₀— 结构受到的力，k—刚度系数，y(τ)—吊篮结构的模态位移；

由于阻尼系数对吊篮结构的影响小，因此忽略阻尼的影响，即阻 尼比为0，阻尼圆周频率等于无阻尼圆周频率即自由振动频率，则解 析解变为：

将式(10)求一阶导数，得到模态速度随时间变化的关系式：

根据初始条件τ＝0和刚度系数k：

k＝mw² (16)

由式(12)、式(13)、式(14)、式(15)和式(16)联立求解，得到模 态位移随时间的变化；

将当前计算时间步吊篮结构中每一个振荡单元的受力初值p₀、速 度初值

和位移初值y₀代入式(17)得到该时间步的模态位移，通过式(2)将模态位移转化为物理坐标下的位移，即得到吊篮各部分的位移；

第十步：将吊篮结构的位移转化为管道的新面积，传递给水力学 计算部分；由于水力学守恒方程中求解需要面积变化项dA/A，所以需 要将吊篮结构的位移转化为面积；转化方法是通过式(18)将投影到一 维方向上的位移重新转化到真实的方向上，则吊篮结构相邻管道的直 径也会变化响应的值，然后计算管道的新面积；

式中：A_i—管道i的横截面积；s_i—管道i对应振荡单元的面积；l_i— 管道i的长度；x_j—振荡单元j的位移；θ_j—振荡单元j的方位角；

第十一步：计算反应堆压力容器及其内部结构和一回路管道的受 力，并分析结构的受力是否满足安全要求；反应堆堆内结构的受力分 为水平方向的力和竖直方向的力，水平方向的力是将所有反应堆内吊 篮结构水平方向上相邻管道的压力投影到x和y方向上然后计算合力， 竖直方向上的力是将反应堆内吊篮结构受到的压力和流体施加的摩 擦力和局部阻力相加；管道受到的力是将管道分为很多段，利用动量 方程—式19计算每一段管道的受力；

式中：F—管道段收到的力，

—管道段i当前时间步的流体质量流 密度，

—管道段i上一时间步的流体质量流密度，A—管道截面积， Δx—管道段长度，Δt—时间步长，

—管道段进口截面流体压力，

—管道段出口截面流体压力，

—管道段进口截面流体质量流 密度，

—管道段出口截面流体质量流密度，

——管道段进口 截面的流速，

——管道段出口截面的流速；

第十二步：判断当前计算时间是否到达结束时间；若没有结束， 则返回到第三步，重新计算；若结束，则结束计算。

Claims (1)

1.一种核电厂破口失水事故水力载荷分析方法，其特征在于，包括如下步骤：

第一步：根据核电厂一回路的运行参数和管道结构参数建立数学模型；

第二步：确定反应堆吊篮的节点划分与周围管道的对应关系，并对吊篮结构进行模态振型分析，得到吊篮结构的固有振动频率和模态振型向量参数；

由于吊篮结构复杂，因此将圆筒状的模型当作板状进行简化处理，并划分成若干个小的振荡单元，而每个小的振荡单元均存在结构动力学方程：

式中：m-吊篮结构的质量；c-吊篮结构的阻尼系数；k-吊篮结构的刚度系数；v(t)-吊篮结构的位移；f(t)-施加于吊篮结构的随时间变化的载荷；

——吊篮结构的速度；

——吊篮结构的加速度；

由于各个吊篮结构振荡单元之间存在相互影响，所以m、k、c都不是对角矩阵，为了消除这种相互影响，简化计算，引入模态振型分析；首先，通过坐标变换，将物理坐标转化为模态坐标：

v(t)＝φ_ny_n(t) (2)

式中：φ_n—吊篮结构的模态矩阵，每一列各代表一种模态的振型；y_n(t)—第n种模态的模态位移；

对于吊篮结构无阻尼自由振动方程组—式(3)，存在振荡单元振动的位移解析解—式(4),其中w—吊篮的固有振动频率，φ—吊篮的振幅,即模态振型矩阵，α—吊篮的振动初始角频率；

v(t)＝φsin(w·t+α) (4)

将式(4)代入式(3)中，得：

M^-1Kφ＝w²φ (5)

式中：M^-1—吊篮质量矩阵的逆矩阵，是n×n阶矩阵，K—吊篮模态刚度矩阵，是n×n阶矩阵，φ—吊篮的模态振型矩阵，是n×h阶矩阵，其中第l行第j列的元素表示第l个振荡单元第j种模态下的振幅，h为预设的模态数，n为振荡单元的个数；

式(5)是一个特征值问题，能够求解得出每一阶模态对应的特征值和特征向量，特征值是吊篮的固有振动频率w的平方，特征向量是该频率下的模态振型向量φ_i，其中φ_i是第i个振动频率对应的模态向量，是n×1阶矩阵；

每个小的振荡单元的结构动力学方程—式(1)组成的方程组转变为一组独立的、以模态坐标和模态参数描述的方程—式(6)，从而实现方程组的独立解耦；

式中：M—吊篮模态质量矩阵，

—吊篮第n种模态的模态加速度，C—吊篮模态阻尼矩阵，

—吊篮第n种模态的模态速度，K—吊篮模态刚度矩阵，F(t)-吊篮模态受力矩阵；

第三步：选用Henry-Fauske模型计算破口处的临界流量，根据流体流动状态判断流型，选用相应流型对应的换热关系式，并计算流体受到的摩擦力；

第四步：使用半隐差分算法求解水力学守恒方程；水力学模型是两流体模型，在水力学守恒方程中考虑了管道的面积变化随时间的变化,如式(7)的汽相质量守恒方程；

式中：A—控制体截面积，α_g—控制体中流体的汽相空泡份额，ρ_g—流体的汽相密度；V_g—流体的汽相速度；Γ_g—控制体中产生的汽相的质量；

水力学模块采用交错网格，在网格中心采用标量网格，网格边界采用矢量网格；对于水力学守恒方程的求解，首先求解两相的动量守恒方程，得到新时刻控制体速度和压力的关系式—式(8)；

式中：

—接管j处液相新时刻的速度，

—接管j处汽相新时刻的速度，

—接管j处液相旧时刻的速度，

—接管j处汽相旧时刻的速度，

—接管j上游的控制体k的压力，P_l ⁿ⁺¹—接管j下游的控制体l的压力，A_f—液相动量方程中接管j处液相旧时刻速度的系数，A_g—汽相动量方程中接管j处汽相旧时刻速度的系数，B_f、B_g分别表示液相、汽相动量方程中接管j上游的控制体k压力的系数，C_f、C_g分别表示液相、汽相接管j下游的控制体l压力的系数，D_f、D_g分别表示液相、汽相动量方程中除旧时刻速度项和上下游压力项外的其他项；

其次求解质量和能量守恒方程，将汽相内能、液相内能、汽相空泡份额和压力四个变量作为求解变量，将汽相质量守恒方程、液相质量守恒方程、汽相能量守恒方程、液相能量守恒方程四个守恒方程的方程组，利用高斯消去法求解方程组，先得到控制体压力和速度的关系式—式(9)；

式中：P_i ⁿ⁺¹—控制体i新时刻的压力，P_i ⁿ—控制体i旧时刻的压力，

—控制体i进口的汽相新时刻的速度，

—控制体i出口的汽相新时刻的速度，

—控制体i进口的液相新时刻的速度，

—控制体i出口的液相新时刻的速度，a、b、c、d分别表示方程组中控制体i液相出、进口速度的系数和方程组中控制体i汽相出、进口速度的系数，e—方程组中除控制体i气液相进出口速度以外的其他项；

然后联立压力和速度的关系式(8)和式(9)，求解得到新时刻的压力和速度；最后通过质量和能量守恒方程方程组，回代新时刻的压力求得气液相内能、汽相空泡份额变量；

第五步：利用水力学守恒方程中求得的控制体的气、液相速度、空泡份额和密度计算每个控制体的流体总质量，总质量除以控制体的体积即该控制体中的流体混合密度ρ_m,1；同时根据控制体的压力和水的物性表查找水和水蒸汽的密度，然后根据式(10)计算控制体中的流体混合密度ρ_m,2；

ρ_m,2＝α_gρ_g+α_fρ_f (10)

第六步：判断流体速度和物性参数即压力、温度、内能、密度的计算是否收敛，若不收敛则时间步长减半并重新进行第三步计算，直至满足收敛准则；

收敛准则为：1、比较利用水力学守恒方程计算得到的控制体中的流体混合密度ρ_m,1和状态方程求得的控制体中的流体混合密度ρ_m,2，二者的相对误差需要满足预设的限值；2、由于水力学守恒方程求解采用半隐算法，因此所有的控制体需要满足声速库伦特限值条件；

第七步：判断外迭代是否收敛；如果不收敛，则进行第八步结构受力计算和第九步位移计算，再通过第九、十步计算管道新的面积和面积变化的相对变化率，然后重复第四步、第五步和第六步；如果收敛，则进行第十一步计算反应堆堆内结构部件的受力计算；收敛准则为前后两次外迭代计算得到的管道面积变化的相对变化率小于规定的限值；因为需要比较前后两次迭代的面积相对变化值，因此在进行每一个时间步长内的首次外迭代计算时，由于没有前一次迭代的计算值，无法进行外迭代收敛判断，因此外迭代默认不收敛，即每个时间步内首次进行外迭代时不进行第七步；

第八步：传递水力学的压力到力学计算部分，求解吊篮结构受到的载荷；将圆筒状的吊篮当作板状，将吊篮周围的控制体的压力投影到垂直于板的方向，将所有的投影值相加即得到吊篮结构的受力；

第九步：利用逐步分段解析法求解吊篮结构的位移；

将时间划分成若干个不大于10^-5秒的时间步长，并且假设一个时间步长内结构受到的力是不变的，由于结构动力学方程是二阶常微分方程，因此存在解析解：

式中：w_D—阻尼圆周频率，w—无阻尼圆周频率，ξ——阻尼比，P₀-结构受到的力，k-刚度系数，y(τ)-吊篮结构的模态位移；

由于阻尼系数对吊篮结构的影响小，因此忽略阻尼的影响，即阻尼比为0，阻尼圆周频率等于无阻尼圆周频率即自由振动频率，则解析解变为：

将式(10)求一阶导数，得到模态速度随时间变化的关系式：

根据初始条件τ＝0和刚度系数k：

k＝mw² (16)

由式(12)、式(13)、式(14)、式(15)和式(16)联立求解，得到模态位移随时间的变化；

将当前计算时间步吊篮结构中每一个振荡单元的受力初值p₀、速度初值

和位移初值y₀代入式(17)得到该时间步的模态位移，通过式(2)将模态位移转化为物理坐标下的位移，即得到吊篮各部分的位移；

第十步：将吊篮结构的位移转化为管道的新面积，传递给水力学计算部分；由于水力学守恒方程中求解需要面积变化项dA/A，所以需要将吊篮结构的位移转化为面积；转化方法是通过式(18)将投影到一维方向上的位移重新转化到真实的方向上，则吊篮结构相邻管道的直径也会变化响应的值，然后计算管道的新面积；

式中：A_i-管道i的横截面积；s_i-管道i对应振荡单元的面积；l_i-管道i的长度；x_j-振荡单元j的位移；θ_j-振荡单元j的方位角；

第十一步：计算反应堆压力容器及其内部结构和一回路管道的受力，并分析结构的受力是否满足安全要求；反应堆堆内结构的受力分为水平方向的力和竖直方向的力，水平方向的力是将所有反应堆内吊篮结构水平方向上相邻管道的压力投影到x和y方向上然后计算合力，竖直方向上的力是将反应堆内吊篮结构受到的压力和流体施加的摩擦力和局部阻力相加；管道受到的力是将管道分为很多段，利用动量方程-式19计算每一段管道的受力；

式中：F-管道段收到的力，

-管道段i当前时间步的流体质量流密度，

-管道段i上一时间步的流体质量流密度，A_p-管道截面积，Δx-管道段长度，Δt-时间步长，

-管道段进口截面流体压力，

-管道段出口截面流体压力，

-管道段进口截面流体质量流密度，

-管道段出口截面流体质量流密度，

——管道段进口截面的流速，

——管道段出口截面的流速；

第十二步：判断当前计算时间是否到达结束时间；若没有结束，则返回到第三步，重新计算；若结束，则结束计算。

Images