[发明专利]核电厂破口失水事故水力载荷分析方法

权利要求书

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-时间步长，-管道段进口截面流体压力，-管道段出口截面流体压力，-管道段进口截面流体质量流密度，-管道段出口截面流体质量流密度，——管道段进口截面的流速，——管道段出口截面的流速；

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