CN104834773B - 一种直管式直流蒸汽发生器换热性能的仿真方法 - Google Patents

Abstract

本发明公开了一种直管式直流蒸汽发生器换热性能的仿真计算模型及其算法。本发明建立了基于传热分区的直流蒸汽发生器一维均相流数学模型，采用适当的判别条件对二次侧换热过程进行划分，在此基础上自主开发了基于MATLAB软件的直流蒸汽发生器换热仿真程序，进行直流蒸汽发生器在不同工况下稳态换热性能的仿真。该发明不仅简化并降低了编程难度，还可实现自动判断分段计算，准确预测不同工况下壁温飞升的幅度、位置及各传热区域传热管的长度，得到直流蒸汽发生器工作恶劣区域。采用该发明可为预防直流蒸汽发生器超温爆管等事故提供理论支持。

Description

一种直管式直流蒸汽发生器换热性能的仿真方法

技术领域

本发明属于核电技术领域，涉及一种直管式直流蒸汽发生器换热性能的仿真方法。

背景技术

直流蒸汽发生器是核反应堆中的重要设备，近年来多被应用于一体化压水堆。由于二次侧工质从过冷水被加热为过热蒸汽的过程中会发生蒸干现象，易造成传热管破裂，进而导致整个反应堆停堆、核辐射泄露等严重问题。因此在核动力装置运行过程中直流蒸汽发生器的安全可靠性引起了大家的高度重视。

1982年Walter应用一阶一维守恒方程，采用漂移流模型描述两相流部分，建立了蒸汽发生器集总参数模型；1997年解衡等人采用可移动边界差分法及节点划分法编制了热工水力程序MOFS，可应用于直管式直流蒸汽发生器的稳态及瞬态计算；2002年付明玉等人采用以集总参数代替分布参数的方法，把每个换热段都看作一个热交换器，串联构成整个蒸发器，以此来研究直流蒸汽发生器静态和动态特性。2013年李娜等人基于集总参数法建立了蒸汽发生器一维均相流动态数学模型，对蒸汽发生器不同工况进行了稳态和动态仿真；

目前对一维模型的分析大多采用集总参数法，将一维问题简化为零维问题能更加简便的进行运算，研究蒸汽发生器的热工水力参数的响应特性，或采用类集总参数法，将系统分为若干个控制体，分别研究各个控制体进、出口的热工水力参数变化规律，且在处理二次侧计算模型时简化了换热区域，无法准确探知直流蒸汽发生器传热管温度的变化规律。模型建立完成之后，既可选择某一软件或语言自行编程进行仿真，例如C语言、Fortran、Matlab/Simulink，也可选择专业仿真软件来进行仿真，例如CFX、ATHOS、RELAP5等。前者编写程序有一定难度，但具有运算速度快等优点；后者得到的仿真结果虽然精确，但却计算速度慢、建模复杂。

发明内容

本发明的目的在于提供一种直管式直流蒸汽发生器换热性能的仿真方法，解决了目前的直流蒸汽发生器的传热计算方法复杂，传热管壁温计算结果不准确的问题。

本发明所采用的技术方案是按照以下步骤：

步骤1：建立直流蒸汽发生器物理模型，设定直流蒸汽发生器结构参数和边界条件；

步骤2：划分二次侧流体的换热区域，得到分段判别条件；

步骤3：计算各个换热区域的对流传热系数；

步骤4：建立直流蒸汽发生器数学模型，以一、二次侧流体、传热管内、外壁为研究对象，以直流蒸汽发生器两个固定管板间距离为计算长度，对微元dx应用能量守恒定律，建立一维均相流传热数学模型；运用有限差分法建立直流蒸汽发生器数学模型的常微分方程；建立五个传热区域的系数矩阵；

步骤5：根据步骤1中的结构参数、边界条件；步骤2中的分段判别条件；步骤3中的各区域对流传热系数；步骤4中的系数矩阵，进行基于MATLAB的仿真。

进一步，所述步骤1中直流蒸汽发生器结构参数包括传热管长度、传热管外径、传热管内径、传热管数目；100％工况的边界条件：一次侧包括入口温度、质量流量、入口压力，二次侧包括入口温度、质量流量、入口压力。

进一步，所述步骤2中将二次侧流体分为预热段、过冷沸腾段、饱和核态沸腾段、强制对流蒸发段、缺液段、过热段，其中由于对流传热系数计算公式的限制，饱和核态沸腾段和强制对流蒸发段作为一个计算区域，称为核态沸腾段。

进一步，所述步骤3中，一次侧流体和二次侧预热段、过热段单相对流换热区采用Dittus-Boelter关联式：

对于流体被加热：

对于流体被冷却：

式中，α为对流传热系数，λ₁为一次侧流体导热系数，Re_f为一次侧雷诺数，Pr_f为一次侧普朗特数；

过冷沸腾换热区采用Pradanovic et al提出的关联式：

Ja^*＝c_pΔT_sub/r (5)

式中，α_tp、α_l分别为两相、液相流体对流传热系数，Bo为沸腾数，Ja*为Jakob数，ρ_g、ρ_f分别为汽体和液体的密度，Pr为普朗特数，

为二次侧质量流速，r为汽化潜热，c_p为定压比热容，ΔT_sub为过冷度；

核态沸腾换热区采用Chen关系式：

α_chen＝α_c+α_NCB (6)

式中，α_chen为强迫对流蒸发和饱和核态沸腾的传热系数，α_NCB为饱和核态沸腾传热系数，α_c为强迫对流传热系数；

式中，σ是汽、液分界面上水的表面张力，k_f为饱和水导热系数，c_pf为饱和水定压比热容，μ_f为饱和水动力粘度，t_w为外壁温度，t_s为饱和温度，p_w是对应壁温t_w的饱和压力，d_e为二次侧当量直径，参数S是核态沸腾抑制因子；

S和F的计算公式如下：

对于S：

对于F：

式中，X_u是Martinelli参数：

而：

缺液换热区采用米罗波利斯基关联式：

其中：

式中，Re_g为二次侧饱和蒸汽的雷诺数，ν_g、ν_l分别为饱和蒸汽和饱和水的运动粘度，Pr_w是以壁温为定性温度的普朗特数。

进一步，所述预热段、过冷沸腾段、核态沸腾段、缺液段、过热段的判断方法为：

在二次侧预热段的对流传热中，是否发生过冷沸腾用Jens-Lottes公式进行判断：

其中，ΔT_s为壁面过热度，p为二次侧压力，q为热流密度；

过冷沸腾段到核态沸腾段判断依据为二次侧流体温度达到对应压力下的饱和温度；

用质量含汽率判断是否进入缺液段：X_DO＝0.3+0.7e^-45ω，

其中X_OD为计算的临界含汽率，ω为计算的无量纲参量，

为二次侧流体质量流速，μ_l为饱和水动力粘度，ρ_l、ρ_g分别为饱和水和饱和蒸汽的密度，σ为汽、液分界面上水的表面张力；

当二次侧流体总焓值达到饱和蒸汽的焓值时进入过热段。

本发明的有益效果是可快速、准确的预测直流蒸汽发生器的传热特性。

附图说明

图1为直流蒸汽发生器物理模型图；

图2为本发明仿真具体流程图；

图3为100％工况下温度的变化。

具体实施方式

下面结合具体实施方式对本发明进行详细说明。

本发明采用分布参数法建立了基于五个计算区域的直流蒸汽发生器数学模型，并在仿真过程中实现自动分段计算以准确分析直流蒸汽发生器的换热特性。

步骤1，建立直流蒸汽发生器物理模型：

设定直流蒸汽发生器结构参数和边界条件，如下表1、表2。

表1 直流蒸汽发生器结构参数

表2 100％工况的边界条件

结合直流蒸汽发生器实际工作工程做出以下假设：1)将直流蒸汽发生器的管束简化成一根单管，认为直流蒸汽发生器中每根单管的性质都相同；2)二次侧流体为均相流，在同一截面汽、液两相速度相同，温度相同；3)工质的物性参数沿横截面方向不变，即采用沿管长方向的一维模型；4)忽略传热管的轴向导热，认为管壁轴向温度仅取决于一、二次侧流体温度变化；5)假设一次侧流体密度不变，忽略一次侧的压降。

根据假设得到直流蒸汽发生器简图如图1所示，图中沿竖直管二次侧流动方向为x轴，在轴向有效长度L中任意选取微元段dx，作为计算单元。二次侧流体温度T₂分别对应一次侧流体温度T₁、管内壁温度T_wi、管外壁温度T_wo。T_1,in、c_p1,in、m₁为一次侧流体入口参数，分别为温度、定压比热容、质量流量；T_1,out、c_p1,out、m₁为一次侧流体出口参数；T_2,in、c_p2,in、m₂为二次侧流体入口参数；T_2,out、c_p2,out、m₂为二次侧流体出口参数。这些参数将作为边界条件在仿真程序中给出。图2为本发明的仿真步骤。图3为100％工况下温度的变化情况。

步骤2，根据步骤1中所得参数，划分二次侧流体的换热区域，得到分段判别条件；

为准确研究直流蒸汽发生器一、二次侧流体以及传热管管壁的温度分布规律，针对直流蒸汽发生器二次侧在换热过程中会发生沸腾放热危机，需采用适当的判别条件将其换热区域进行划分，才能准确得到壁温在干涸点处的飞升幅度。

将二次侧流体分为预热段、过冷沸腾段、饱和核态沸腾段、强制对流蒸发段、缺液段、过热段，其中由于对流传热系数计算公式的限制，饱和核态沸腾段和强制对流蒸发段作为一个计算区域，称为核态沸腾段。

在二次侧预热段的对流传热中，是否发生过冷沸腾可用Jens-Lottes公式进行判断：

其中，ΔT_s为壁面过热度，℃；p为二次侧压力，MPa；q为热流密度，MW/m²。

过冷沸腾段到核态沸腾段判断依据为二次侧流体温度达到对应压力下的饱和温度。

接下来用质量含汽率判断是否到达干涸点，是否进入缺液段：X_DO＝0.3+0.7e^-45ω，

其中X_OD为计算的临界含汽率；ω为计算的无量纲参量；

为二次侧流体质量流速，kg/m²s；μ_l为饱和水动力粘度，Pa·s；ρ_l、ρ_g分别为饱和水和饱和蒸汽的密度，kg/m³；σ为汽、液分界面上水的表面张力，N/m。

当二次侧流体总焓值达到饱和蒸汽的焓值时进入过热段。

根据划分二次侧流体换热区域的计算公式所得判别条件如表3所示。

表3 各个传热区域分段的判别条件

步骤3，根据步骤1中所得参数，计算5个区域的对流传热系数：

一次侧流体和二次侧预热段、过热段单相对流换热区采用Dittus-Boelter关联式：

对于流体被加热：

对于流体被冷却：

式中，α为对流传热系数，W/(m²·℃)；λ₁为一次侧流体导热系数，W/(m·℃)；Re_f为一次侧雷诺数；Pr_f为一次侧普朗特数，各准则数的物性参数取流体的算术平均温度作为定性温度。

过冷沸腾换热区采用Pradanovic et al提出的关联式：

Ja^*＝c_pΔT_sub/r (5)

式(3)～(5)中，α_tp、α_l分别为两相、液相流体对流传热系数，W/(m²·℃)；Bo为沸腾数；Ja*为Jakob数；ρ_g、ρ_f分别为汽体和液体的密度，kg/m³；Pr为普朗特数；

为二次侧质量流速，kg/m²s；r为汽化潜热，J/kg；c_p为定压比热容，J/(kg·℃)；ΔT_sub为过冷度，℃。

核态沸腾换热区采用Chen关系式：

α_chen＝α_c+α_NCB (6)

式中，α_chen为强迫对流蒸发和饱和核态沸腾的传热系数(简称核态沸腾传热系数)，W/(m²·℃)；α_NCB为饱和核态沸腾传热系数，W/(m²·℃)；α_c为强迫对流传热系数，W/(m²·℃)。

式中，σ是汽、液分界面上水的表面张力，N/m；k_f为饱和水导热系数，W/(m·℃)；c_pf为饱和水定压比热容，J/(kg·℃)；μ_f为饱和水动力粘度，Pa·s；t_w为外壁温度，℃；t_s为饱和温度，℃；p_w是对应壁温t_w的饱和压力，Pa；d_e为二次侧当量直径，m；参数S是核态沸腾抑制因子，它和参数F都是图解函数。S和F的计算公式如下：

对于S：

对于F：

式中，X_u是Martinelli参数：

而：

缺液换热区采用米罗波利斯基关联式：

其中：

式中，Re_g为二次侧饱和蒸汽的雷诺数；ν_g、ν_l分别为饱和蒸汽和饱和水的运动粘度，m²/s；Pr_w是以壁温为定性温度的普朗特数。

根据公式(1)～(14)所得各个区域的对流传热系数如表4所示。

表4 各个计算区域的对流传热系数

步骤4，建立直流蒸汽发生器数学模型，分别以一、二次侧流体、传热管内、外壁为研究对象，以直流蒸汽发生器两个固定管板间距离为计算长度，对微元dx应用能量守恒定律，建立一维均相流传热数学模型。

(1)一次侧流体：

(2)二次侧流体：

①预热段：

②过冷沸腾段：

③核态沸腾段：

④缺液段：

⑤过热段：

(3)传热管壁：

传热管内壁：

由于二次侧流体的换热计算被分为五个部分。因此相对应的传热管外壁的热平衡方程也被分为五个部分，其中K₂随传热区域的不同而不同。

式(15)～(22)中，x为沿传热管轴向的位置坐标；t为时间；n为传热管根数；T₁为一次侧流体温度，℃；T₂为二次侧流体温度，℃；T_2s为二次侧工作压力下对应的饱和温度，℃；T_wi、T_wo为传热管内、外壁温度，℃；K₁、K₂为一二次侧对流换热系数，W/(m²·℃)；c_p1为一次侧流体定压比热容，J/(kg·℃)；c_pw为直管壁的定压比热容，J/(kg·℃)；c_p2y、c_p2l、c_p2g分别为二次侧预热段、过冷沸腾段、过热段定压比热容，J/(kg·℃)；d_i、d_o为传热管内、外径，m；λ_w为管壁的导热系数，W/(m·℃)；m₁、m₂一、二次侧流体的质量流量，kg/s。M₁为一次侧流体沿轴向单位长度的质量，kg/m；M_2y、M_2h、M_2m、M_2g、M_w为二次侧预热段、核态沸腾段、缺液段、过热段流体及管壁沿轴向单位长度的质量，kg/m；h₂为二回路沸腾段的焓值，J/kg。

步骤5：运用有限差分法对式(15)～(22)因变量的空间项进行离散，将偏微分方程组转化为常微分方程组。即将传热管沿x轴方向离散为N份，可得到一组与时间相关的变量T_i(t)，(i＝1，2，3......，N)，各段断点用x(i)表示。第i段长度用Δx(i)表示或统一表示成Δx。运用有限差分法，采用对i点的向后差分则有：

可得一次侧流体常微分方程：

同理可得传热管内、外壁及二次侧流体常微分方程。

步骤6，根据步骤5中所得常微分方程组，建立五个传热区域的系数矩阵：

根据最终得到4N个常微分方程，将得到的微分方程组等号右边项表达成状态空间的形式：

等式(24)中，X代表一、二次侧温度，传热管内、外壁温度或二次侧饱和沸腾段焓值项，即X＝[T₁,T_wi,T_wo,T₂(h₂)]；A(4×4)为X的系数矩阵，B(4×2)为U的系数矩阵。

U为边界条件项，即U＝[T_1,i-1,T_2,i-1(h_2,i-1)]^T。根据实际情况，预热段的边界条件可表示为：T_1,i-1＝T_1,0＝T_1,out，T_2,i-1＝T_2,0＝T_2,in，h_2,i-1＝h_2s由T_2s确定。为了便于Matlab仿真编程，将各个传热区域的系数矩阵整理如下：

1.预热段：

a₁₃＝0，a₁₄＝0

a₂₄＝0

a₃₁＝0，

a₄₁＝0，a₄₂＝0，

b₁₂＝0，b₂₁＝0，b₂₂＝0，b₃₁＝0，b₃₂＝0，b₄₁＝0，

2.过冷沸腾段：

a₁₃＝0，a₁₄＝0

a₂₄＝0

a₃₁＝0，

a₄₁＝0，a₄₂＝0，

b₁₂＝0，b₂₁＝0，b₂₂＝0，b₃₁＝0，b₃₂＝0，b₄₁＝0，

3.核态沸腾段：

a₁₃＝0，a₁₄＝0

a₂₄＝0

a₃₁＝0，

a₃₄＝0

a₄₁＝0，a₄₂＝0，

b₁₂＝b₁₃＝0

b₂₁＝b₂₂＝b₂₃＝0，b₃₁＝b₃₂＝0，b₄₁＝0，

4.缺液段：

a₁₃＝0，a₁₄＝0

a₂₄＝ 0a₃₁＝0，

a₃₄＝0

a₄₁＝0，a₄₂＝0，

b₁₂＝b₁₃＝0

b₂₁＝b₂₂＝b₂₃＝0，b₃₁＝b₃₂＝0，b₄₁＝0，

5.过热段：

a₁₃＝0，a₁₄＝0

a₂₄＝0

a₃₁＝0，

a₄₁＝0，a₄₂＝0，

b₁₂＝0，b₂₁＝0，b₂₂＝0，b₃₁＝0，b₃₂＝0，b₄₁＝0，

步骤7，根据步骤1中的结构参数、边界条件；步骤2中的分段判别条件；步骤3中的各区域对流传热系数；步骤6中的系数矩阵，进行基于MATLAB的仿真；

将步骤1中的结构参数、边界条件；步骤3中的各区域对流传热系数作为初始参数在程序中给出，再分段求解系数矩阵最终得到温度分布曲线。

由于直流蒸汽发生器二次侧传热过程复杂，两相流动具有多种变化，因此需要在仿真计算过程中能够自动进行分段计算，即当所得结果满足一定条件时，结束一个换热区域的循环进入下一个换热区域的循环计算。在每一段的循环中采用Jacobi迭代法求解步骤6中的矩阵，再根据步骤2所得的判别数据运用if语句对结果进行判断，若所得结果不满足判别条件则继续循环计算，若所得结果满足判别条件则采用break语句跳出循环进行下一步。这样编程既提高了仿真结果的准确性，又减短了运行时间。

由于该直流蒸汽发生器一次侧流动方向是从上到下流动，而二次侧流动方向是从下到上流动。因此所给的边界条件不在直流蒸汽发生器的同一端，这样将使所列方程无法进行分段求解。对此，本文采用假设的方法，先将二次侧流体入口温度T_2,in和假设的一次侧流体出口温度T′_1,out作为边界条件带入循环计算，再把最后计算得到的一次侧流体入口温度T′_1,in与已知的一次侧流体入口温度T_1,in作对比，若两者之间的误差大于规定误差范围，则减小或增大假设的一次侧流体出口温度T′_1,out再进行计算，若两者之间误差满足误差范围，则结束仿真得到最终结果，以此来确保结果的准确性。按以上数据所得仿真结果如图3所示。

和现有技术相比，本发明具有如下优点和有益效果：

对直流蒸汽发生器物理模型进行简化后降低了编程难度，使得运算直观、快速，同时又针对直流蒸汽发生器的特点注重了二次侧换热区域的划分，可快速、准确的预测直流蒸汽发生器的传热特性，得到壁温在干涸点处的飞升幅度和直流蒸汽发生器传热管工作条件恶劣区域。

Claims (1)

1.一种直管式直流蒸汽发生器换热性能的仿真方法，其特征在于按照以下步骤进行：

步骤1：建立直流蒸汽发生器物理模型，设定直流蒸汽发生器结构参数和边界条件；

步骤2：根据分段判别条件，划分二次侧流体的换热区域；

步骤3：计算各个换热区域的对流传热系数；

步骤4：建立直流蒸汽发生器数学模型，以一、二次侧流体、传热管内、外壁为研究对象，以直流蒸汽发生器两个固定管板间距离为计算长度，对微元dx应用能量守恒定律，建立一维均相流传热数学模型；运用有限差分法建立直流蒸汽发生器数学模型的常微分方程；建立五个传热区域的系数矩阵；

步骤5：根据步骤1中的结构参数、边界条件；步骤2中的分段判别条件；步骤3中的各区域对流传热系数；步骤4中的系数矩阵，进行基于MATLAB的仿真；

所述步骤1中直流蒸汽发生器结构参数包括传热管长度、传热管外径、传热管内径、传热管数目；边界条件包括：一次侧包括入口温度、质量流量、入口压力，二次侧包括入口温度、质量流量、入口压力；

所述步骤2中将二次侧流体分为预热段、过冷沸腾段、饱和核态沸腾段、强制对流蒸发段、缺液段、过热段，其中由于对流传热系数计算公式的限制，饱和核态沸腾段和强制对流蒸发段作为一个计算区域，称为核态沸腾段；

所述步骤3中，一次侧流体和二次侧预热段、过热段单相对流换热区采用Dittus-Boelter关联式：

对于流体被加热：

对于流体被冷却：

式中，α为对流传热系数，λ₁为流体导热系数，d_i为水力直径，Re_f为雷诺数，Pr_f为普朗特数；

过冷沸腾换热区采用Pradanovic et al提出的关联式：

Ja^*＝c_pΔT_sub/r (5)

α_tp、α_l分别为两相、液相流体对流传热系数，Bo为沸腾数，q为热流密度，Ja*为Jakob数；ρ_g、ρ_f分别为汽体和液体的密度，Pr为普朗特数，

为二次侧质量流速，r为汽化潜热，c_p为定压比热容，ΔT_sub为过冷度；

核态沸腾换热区采用Chen关系式：

α_chen＝α_c+α_NCB (6)

式中，α_chen为强迫对流蒸发和饱和核态沸腾的传热系数，α_NCB为饱和核态沸腾传热系数，α_c为强迫对流传热系数；

式中，σ是汽、液分界面上水的表面张力，k_f为饱和水导热系数，c_pf为饱和水定压比热容，μ_f为饱和水动力粘度，t_w为外壁温度，t_s为饱和温度，p_w是对应壁温t_w的饱和压力，p是二次侧压力，d_e为二次侧当量直径，参数S是核态沸腾抑制因子；G为质量流率，x为质量含汽率，k_f为饱和水导热系数；

S和F的计算公式如下：

对于S：

对于F：

式中，X_tt是Martinelli参数：

而：

缺液换热区采用米罗波利斯基关联式：

其中：

式中，Re_g为二次侧饱和蒸汽的雷诺数，Re′_tp为汽液两相的雷诺数，μ_g为气体的动力粘度，ν_g、ν_l分别为饱和蒸汽和饱和水的运动粘度，Pr_w是以壁温为定性温度的普朗特数；

所述步骤2中预热段、过冷沸腾段、核态沸腾段、缺液段、过热段的判断方法为：

在二次侧预热段的对流传热中，是否发生过冷沸腾用Jens-Lottes公式进行判断：

其中，ΔT_s为壁面过热度，p为二次侧压力，q为热流密度；

过冷沸腾段到核态沸腾段判断依据为二次侧流体温度达到对应压力下的饱和温度；

用质量含汽率判断是否进入缺液段：X_DO＝0.3+0.7e^-45ω，

其中X_DO为计算的临界含汽率，ω为计算的无量纲参量，

为二次侧流体质量流速，μ_l为饱和水动力粘度，ρ_l、ρ_g分别为饱和水和饱和蒸汽的密度，σ为汽、液分界面上水的表面张力；

当二次侧流体总焓值达到饱和蒸汽的焓值时进入过热段。

Images