CN111859752B - 一种核反应堆蒸汽发生器流致振动计算方法 - Google Patents

Abstract

本发明公开了一种核反应堆蒸汽发生器流致振动计算方法，步骤如下：1、建立蒸汽发生器传热管湿模态分析模型，开展湿模态分析并提取传热管的质量矩阵[M],刚度矩阵[K]和阻尼矩阵[C]；2、将求解瞬态动力学方程的Newmark‑β方法和动网格模型写成用户自定义函数文件；3、建立传热管内外一、二次侧流体域流体动力学计算模型，并开展第一时间步迭代计算；4、调用并执行第2步中用户自定义函数文件，流体域网格进行更新，网格更新后的流体域进行下一时间步的迭代计算；5、循环执行第4步，直到计算到设定好的计算终止时间，计算停止。通过计算获得核反应堆蒸汽发生器的流致振动特性，对核反应堆蒸汽发生器的设计和安全分析具有重要意义。

Description

一种核反应堆蒸汽发生器流致振动计算方法

技术领域

本发明属于核反应堆蒸汽发生器流致振动技术领域，具体涉及一种核反应堆蒸汽发生器流致振动计算方法。

背景技术

蒸汽发生器是压水堆核电站连接一、二回路的核心部件，其有效性关乎核反应堆的安全、稳定运行。压水堆核电站目前使用的蒸汽发生器主要有立式U型传热管自然循环蒸汽发生器(UTSG)和管式直流蒸汽发生器(OTSG)。对立式U型传热管自然循环蒸汽发生器而言，一回路的高压热流体由热腿经进口腔室流入U型传热管，在U型传热管内向二回路传热并通过出口腔室经冷腿流回压力容器。U型传热管作为蒸汽发生器中的关键部件，承担一回路压力边界和传热边界的作用，是蒸汽发生器事故的多发位置。在高温高压的环境中，二次侧流体冲刷U型传热管管束，会引起传热管振动，即流致振动现象，还会导致传热管冲击、磨损和应力腐蚀，这些现象会破坏传热管的完整性，还可能进一步导致传热管的破裂，造成冷却剂泄露，破坏蒸汽发生器的安全性和完整性，严重影响压水堆核电站的安全。传热管损坏的原因有很多，大多数损坏是由冷却剂的流致振动现象引起，包括流弹性失稳、非定常压力波动激励、漩涡脱落，以及湍流激振，引起传热管恶化的主要振动原因来自于二次侧流体的激振作用。传热管的流致振动现象涉及流体和固体之间的相互作用，伴随着流固耦合面的数据传递。随着计算机技术的发展，研究人员越来越多的采用数值模拟的方法来研究蒸汽发生器的热工水力特性以及传热管的流致振动现象，因此，利用流固耦合方法，对蒸汽发生器传热管的流致振动现象进行研究，具有重要意义。流固耦合方法分为单向、双向流固耦合方法，目前，研究人员大多数采用中间数据交换平台耦合流体动力学计算程序和结构动力学分析程序的方法，实现流体域和固体域实时双向的数据传递。这种方法伴随着计算量大，耦合面数据交换速度慢且容易出错的弊端，且对计算程序有一定的适用性要求和版本要求，尤其当计算对象几何结构复杂的时候，双向流固耦合工作的开展将受限于计算资源。单向流固耦合虽然不能像双向流固耦合那样实现流固耦合交界面实时双向的数据传递，但是当固体形变较小，且我们更加关注流体对固体的影响时，采用单向流固耦合是更加合适的选择。并且，如果能够在一个计算程序中同时实现流体动力学方程的求解和瞬态动力学方程的求解，则会大大减少整个计算程序的依赖性。

发明内容

本发明提供了一种核反应堆蒸汽发生器流致振动计算方法，该方法通过用户自定义函数文件的方式，实现了在流体动力学计算程序中，求解瞬态动力学方程，实现流固耦合面的数据传递，模拟蒸汽发生器U型传热管流致振动特性。

为达到上述目的，本发明采用了一下技术方案：

一种核反应堆蒸汽发生器流致振动计算方法，针对处于核反应堆一、二次侧流体环境中的蒸汽发生器U型传热管的流致振动现象，通过湿模态分析求取U型传热管的质量矩阵[M]，刚度矩阵[K]和阻尼矩阵[C]，通过流体动力学计算程序调用并执行用户自定义函数文件求解瞬态动力学方程，通过动网格模型更新每一时间步之后的网格，实现蒸汽发生器U型传热管在一、二次侧流体共同作用下的流致振动计算：包括以下步骤：

步骤1：获取核反应堆蒸汽发生器U型传热管几何参数，建立单根U型传热管固体域、管内一次侧流体域、管外二次侧流体域的简化几何模型装配体，在有限元分析程序中，赋予U型传热管Inconel690材料物性，利用简支和固支方式，在U型传热管两端根部施加固支约束，在U型传热管弯头30°、60°、120°、150°方向施加简支约束，模拟防振条对U型传热管的作用，根据核反应堆一、二次侧流体的实际运行条件，施加流体域的进出口边界条件，建立U型传热管湿模态分析模型，利用有限元分析程序开展湿模态分析并提取传热管的质量矩阵[M],刚度矩阵[K]和阻尼矩阵[C]；

步骤2：根据步骤1中提取的质量矩阵[M]、刚度矩阵[K]和阻尼矩阵[C]，将求解瞬态动力学方程的Newmark-β方法、动网格模型添加进用户自定义函数文件，具体步骤如下：

步骤2-1：根据步骤1中提取的质量矩阵[M]、刚度矩阵[K]和阻尼矩阵[C]，在用户自定义函数文件中定义这三者；

步骤2-2：给定U型传热管的初始位移{u}₀、初始速度

初始加速度

初始外力{F}均为零，即{u}₀＝0,

{F}₀＝0；

步骤2-3：根据Newmark-β方法假设：

{u}，

{F}分别表示U型传热管的位移、速度、加速度和外力，下标表示不同时刻，Δt表示时间步长，β和γ为自定义参数，选择β＝1/2，γ＝1/4，此时Newmark-β方法无条件稳定，由公式(1)、公式(2)推导，并代入β和γ的值，得到：

t+Δt时刻的瞬态动力学方程：

将公式(3)和公式(4)代入公式(5)得到

其中：

公式(5)中的[M]、[C]、[K]分别为质量矩阵、阻尼矩阵、刚度矩阵，{F}_t+Δt为t+Δt时刻的外力，公式(6)为t+Δt时刻的瞬态动力学方程的等效表达式，其中

为等效刚度矩阵，其表完整达式为公式(7)，

为t+Δt时刻的等效外力，其完整表达式为公式(8)；在用户自定义函数文件中，定义Δt＝0.001s，定义获取U型传热管表面外力的函数，定义公式(3)、(4)、(6)、(7)、(8)，根据这些定义，可以根据t时刻的{u}_t，

可求得t+Δt时刻{u}_t+Δt，

步骤2-4：添加动网格模型到用户自定义函数文件中，用来实现网格更新功能；

步骤3：建立单根U型传热管内一次侧、管外二次侧流体域计算模型，在一、二次侧流体域计算模型中标识出和固体U型传热管内、外壁接触的耦合面，定义动网格边界，选择动网格更新方法，利用流体动力学计算程序开展第一时间步迭代计算至收敛(计算残差小于10^-5)；

步骤4：调用步骤2中用户自定义函数文件，流体域网格进行更新，再对网格更新后的流体域进行下一时间步的迭代计算；

步骤5：循环执行第4步，确保每一时间步收敛(计算残差小于10^-5)直到计算到设定好的计算终止时间，计算停止。

本发明方法可以实现在流体动力学计算程序中，同时求解流体动力学方程和瞬态动力学方程，实现蒸汽发生器传热管流致振动特性的流固耦合方法，此方法可以为核反应堆蒸汽发生器的设计和安全分析提供参考，具有重要意义。和现有技术相比较，本发明具有以下优点和有益效果：

1.该方法提供了一种核反应堆蒸汽发生器流致振动计算方法，该方法适用于现有的大多数流体动力学计算程序，如FLUENT、STAR-CCM+、OPENFOAM等；

2.该方法基于流固耦合方法和动网格模型，通过压力场、位移场之间的交互耦合，在模拟流体侧三维流动状态的同时，可以模拟固体侧流致振动现象；

3.通过用户自定义函数的形式，在流体动力学计算程序中实现了固体结构瞬态动力学方程的Newmark-β方法求解过程，省去了对有限元计算程序的依赖性，避免了中间数据交换接口数据交换的复杂性和易错性，且用户自定义函数文件具有很好的移植性和适应性；

4.使用湿模态分析方法对处于一、二次侧流体域中的蒸汽发生器U型传热管进行湿模态分析，不光可以得到核反应堆实际运行条件下，处于一、二次侧冷却剂中U型管传热管的质量矩阵[M]、刚度矩阵[K]和阻尼矩阵[C]，还可以得到U型传热管不同阶的振动频率和振动幅值，可以为U型传热管的机械力学分析，提供参考；

5.该方法可以对不同运行工况、不同高度、以及不同约束的U型传热管进行流致振动特性流固耦合模拟；

本发明提出的核反应堆蒸汽发生器流致振动计算方法适用于不同结构参数的U型传热管，还可以根据不同工况，通过修改边界条件的初始条件进行事故状态下的流致振动现象的计算。

附图说明

图1为本发明方法的流程图。

图2为单根U型传热管固体域、管内一次侧流体域、管外二次侧流体域简化几何模型装配体。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步详细描述：

本发明提供了一种如图1所示核反应堆蒸汽发生器流致振动计算方法，具体方法如下：

步骤1：获取核反应堆蒸汽发生器U型传热管几何参数，考虑到U型管排列的对称性和重复性，建立单根U型管固体域、管内一次侧流体域、管外二次侧流体域的简化几何模型装配体，如图2所示。在有限元分析程序中，赋予传热管Inconel690材料物性，利用简支和固支方式，在U型传热管两端根部施加固支约束，在U型传热管弯头30°、60°、120°、150°方向施加简支约束，模拟防振条对U型传热管的作用，根据核反应堆一、二次侧流体的实际运行条件，施加流体域的进出口边界条件，建立U型传热管湿模态分析模型，利用有限元分析程序开展湿模态分析并提取传热管的质量矩阵[M],刚度矩阵[K]和阻尼矩阵[C]；

步骤2：根据步骤2中提取的质量矩阵[M]、刚度矩阵[K]和阻尼矩阵[C]，将瞬态动力学方程的Newmark-β方法求解过程、动网格模型添加进进用户自定义函数文件，具体步骤如下：

步骤2-1：根据步骤1中提取的质量矩阵[M]、刚度矩阵[K]和阻尼矩阵[C]，在用户自定义函数文件中定义这三者；

步骤2-2：在用户自定义函数文件中定义U型传热管的初始位移{u}₀、初始速度

初始加速度

初始外力{F}₀均为零，即{u}₀＝0,

{F}₀＝0；

步骤2-3：根据Newmark-β方法假设：

{u}，

{F}分别表示U型传热管的位移、速度、加速度和外力，下标表示不同时刻，Δt表示时间步长，β和γ为自定义参数，选择β＝1/2，α＝1/4，此时Newmark-β方法无条件稳定，由公式(1)、公式(2)推导，并代入β和γ的值，得到：

t+Δt时刻的瞬态动力学方程：

将公式(3)和公式(4)代入公式(5)得

其中：

公式(5)中的[M]、[C]、[K]分别为质量矩阵、阻尼矩阵、刚度矩阵，{F}_t+Δt为t+Δt时刻的外力，公式(6)为t+Δt时刻的瞬态动力学方程的等效表达式，其中

为等效刚度矩阵，其表完整达式为公式(7)，

为t+Δt时刻的等效外力，其完整表达式为公式(8)；在用户自定义函数文件中，定义Δt＝0.001s，定义获取U型传热管表面外力的函数，定义公式公式(3)、(4)、(6)、(7)、(8)，根据这些定义，可以根据t时刻的{u}_t，

可求得t+Δt时刻{u}_t+Δt，

步骤2-4：添加动网格模型到用户自定义函数文件中，用来实现网格更新功能；

步骤3：建立单根U型传热管内一次侧、管外二次侧流体域计算模型，在一、二次侧流体域计算模型中标识出和固体U型传热管内、外壁接触的耦合面，定义动网格边界，选择动网格更新方法，利用流体动力学计算程序开展第一时间步迭代计算至收敛(计算残差小于10^-5)；

步骤4：调用步骤2中用户自定义函数文件，流体域网格进行更新，再对网格更新后的流体域进行下一时间步的迭代计算；

步骤5：循环执行第4步，确保每一时间步收敛(计算残差小于10^-5)直到计算到设定好的计算终止时间，计算停止。

Claims (1)

1.一种核反应堆蒸汽发生器流致振动计算方法，其特征在于：针对处于核反应堆一、二次侧流体环境中的蒸汽发生器U型传热管的流致振动现象，通过湿模态分析求取U型传热管的质量矩阵[M]，刚度矩阵[K]和阻尼矩阵[C]，通过流体动力学计算程序调用并执行用户自定义函数文件求解瞬态动力学方程，通过动网格模型更新每一时间步之后的网格，实现蒸汽发生器U型传热管在一、二次侧流体共同作用下的流致振动计算：

该方法包括以下步骤:

步骤1：获取核反应堆蒸汽发生器U型传热管几何参数，建立单根U型传热管固体域、管内一次侧流体域、管外二次侧流体域的简化几何模型装配体，在有限元分析程序中，赋予U型传热管Inconel690材料物性，利用简支和固支方式，在U型传热管两端根部施加固支约束，在U型传热管弯头30°、60°、120°、150°方向施加简支约束，模拟防振条对U型传热管的作用，根据核反应堆一、二次侧流体的实际运行条件，施加流体域的进出口边界条件，建立U型传热管湿模态分析模型，利用有限元分析程序开展湿模态分析并提取传热管的质量矩阵[M],刚度矩阵[K]和阻尼矩阵[C]；

步骤2：根据步骤1中提取的质量矩阵[M]、刚度矩阵[K]和阻尼矩阵[C]，将求解瞬态动力学方程的Newmark-β方法、动网格模型添加进用户自定义函数文件，具体步骤如下：

步骤2-1：根据步骤1中提取的质量矩阵[M]、刚度矩阵[K]和阻尼矩阵[C]，在用户自定义函数文件中定义这三者；

步骤2-2：给定U型传热管的初始位移{u}₀、初始速度

初始加速度

初始外力{F}₀均为零，即{u}₀＝0,

{F}₀＝0；

步骤2-3：根据Newmark-β方法假设：

{u}，

{F}分别表示U型传热管的位移、速度、加速度和外力，下标表示不同时刻，Δt表示时间步长，β和γ为自定义参数，选择β＝1/2，γ＝1/4，此时Newmark-β方法无条件稳定，由公式(1)、公式(2)推导，并代入β和γ的值，得到：

t+Δt时刻的瞬态动力学方程：

将公式(3)和公式(4)代入公式(5)得到

其中：

公式(5)中的[M]、[C]、[K]分别为质量矩阵、阻尼矩阵、刚度矩阵，{F}_t+Δt为t+Δt时刻的外力，公式(6)为t+Δt时刻的瞬态动力学方程的等效表达式，其中

为等效刚度矩阵，其完整表达式为公式(7)，

为t+Δt时刻的等效外力，其完整表达式为公式(8)；在用户自定义函数文件中，定义Δt＝0.001s，定义获取U型传热管表面外力的函数，定义公式(3)、(4)、(6)、(7)、(8)，根据这些定义，能够根据t时刻的{u}_t，

可求得t+Δt时刻{u}_t+Δt，

步骤2-4：添加动网格模型到用户自定义函数文件中，用来实现网格更新功能；

步骤3：建立单根U型传热管内一次侧、管外二次侧流体域计算模型，在一、二次侧流体域计算模型中标识出和固体U型传热管内、外壁接触的耦合面，定义动网格边界，选择动网格更新方法，利用流体动力学计算程序开展第一时间步迭代计算至收敛，收敛条件为计算残差小于10^-5；

步骤4：调用步骤2中用户自定义函数文件，流体域网格进行更新，再对网格更新后的流体域进行下一时间步的迭代计算；

步骤5：循环执行第4步，确保每一时间步收敛，直到计算到设定好的计算终止时间，计算停止。

Images