基于Modelica的并联多通道系统流动不稳定性仿真研究

童秋实，邓康杰，吕 星，曾小康，王艳林，杨 浩，陈 路

（1.苏州同元软控信息技术有限公司，江苏 苏州 215123；2.中国核动力研究设计院中核核反应堆热工水力技术重点实验室，四川 成都 610213）

0 引言

流动不稳定性（flow instability）指系统受到瞬时扰动后偏离原来稳定状态，热工参数发生非周期性漂移或周期性振荡现象［1］。对于存在多个并联加热通道的系统，加热通道间可能存在并联通道之间的流动不稳定性现象，表现为各加热通道间的流量发生周期性脉动，此时将导致临界热流密度（critical heat flux）明显下降或对设备产生机械振荡，从而危及系统、设备的安全运行。因此，在两相流研究中并联加热通道流动不稳定性现象一直是重要的研究课题。

1 相关研究

针对并联加热多通道系统的流动不稳定性现象，许多学者开展了实验研究。Cheng 等［2］通过实验研究自然循环棒束形并联通道流动不稳定，主要研究稳压器连接在加热通道上游时，出现的压力降流动不稳定性的实验现象和产生机理，并研究入口过冷度和热流密度对流动不稳定性的影响。唐瑜等［3］在强迫循环条件下采用矩形并联通道开展流动不稳定性的实验研究，考虑了系统压力、质量流量、入口过冷度等热工参数对流动不稳定性界限的影响，通过实验数据拟合得到了流动不稳定性无量纲关系式。Taleyarhan 等［4］建立平行加热通道系统模型考虑过冷沸腾、热通量分布和管形阻力系数等因素，描述沸水堆中的热工水力现象，仿真分析表明在平行加热通道中功率、流量会对流动不稳定性发生存在较大影响，实验数据证实了该模型对流动不稳定性边界的预测能力。

目前，主要采用实验方式研究并联加热通道流动不稳定性现象［5］，但实验研究周期长、实验工况有限、无法较好捕捉流动不稳定性内部机理现象等。因此，有必要采用仿真手段进行流动不稳定性现象模拟，通过不同工况条件下并联加热通道的快速仿真模拟，研究流动不稳定机理与确定流动不稳定性边界，为实验研究提供参考依据。

在工程实践中，通常采用一维系统仿真程序对并联加热通道的流动不稳定性现象进行仿真分析，由于强迫循环下并联加热通道内流体状态变化剧烈，因此一般将流道分为若干个控制体，该方法相较于研究流道内流场温度场的三维数值模拟方法既降低了计算复杂性，又表达了流动不稳定性热工水力特性，现已广泛应用于流动不稳定性的建模仿真分析中［6-7］。常用的仿真程序一般采用Fortran、C语言、Matlab/Simulink 等语言进行开发，或直接采用RELAP5、THEATRE 等核工程专业仿真软件进行计算［8-9］。然而，上述程序和软件采用卡片式建模，建模过程缺乏图形化界面，用户难以直观获得系统设备的构成、连接关系、参数信息等。此外，美国MMS、芬兰APROS、中国STAR-90 等图形化仿真软件已广泛应用于动力系统仿真领域，通过拖拽式建模及设备参数可视化配置的方式极大降低了建模工作量，但商业软件代码相对封闭，不便于用户根据自身需求进行二次开发。

为此，苏州同元软控信息技术有限公司联合中国核动力研究设计院构建了两相热工水力特性模型架构，采用两流体六方程数学模型开发了基于Modelica 多专业、多层次、多系统的核反应堆统一建模与分析平台NUMAP（Nuclear reactor Unified Modeling and Analysis Platform）［10］。具体为，基于Modelica 语言开发的NUMAP 可参考系统拓扑图进行拖拽式建模，相较于传统热工水力软件具备强大、快捷的可视化交互功能，能极大提升建模效率［11-12］。本文参考某实验台架，在NUMAP 平台上建立高温高压并联加热双通道系统模型，通过分析在均匀、非均匀加热工况下的流动不稳定性、实验台架和NUMAP 系统模型的出口质量含气率值，验证了NUMAP 平台可用于模拟分析并联加热通道流动不稳定性瞬态特性。

2 基于Modelica的两流体六方程理论模型

NUMAP 热工水力模型库基于Modelica 语言，采用两流体六方程进行描述。两流体六方程模型是核反应堆热工水力系统分析常用的数学方程模型，能较为准确地模拟核反应堆系统在正常运行工况、事故工况下汽液两相的传质传热行为，可用于核反应堆系统的设计验证、运行模拟、安全分析等场景。

基于Modelica 语言两流体六方程模型考虑了实际两相流体流动过程中，两相流体具备的不同物性、流速、温度，两相之间的质量、能量和动量交换，即针对汽相和液相分别建立质量、动量和能量守恒方程。为了使控制方程组封闭，还增加了相间界面摩擦、相间传热、相间质量交换、壁面摩擦、壁面传热等本构方程。

2.1 守恒方程模型

汽液质量守恒方程为：

汽液动量守恒方程为：

汽液能量守恒方程为：

式中：k表示汽相g或液相f；i表示界面；α表示体积百分比；μ表示流速；h表示焓值；Γk表示界面传质；Fwk表示壁面摩擦；Fik表示界面摩擦；Fa表示外力场作用；qik表示界面传热；qwk表示壁面传热；Γk、Fwk、Fik、Fa、qik、qwk均为未知量，需要通过本构封闭方程获得。

2.2 本构封闭方程模型

NUMAP 采用的热工水力守恒方程，在数学上封闭求解时需补充本构封闭方程模型，包括流动阻力系数模型、换热系数模型、流型图及流型转换判定条件、水和蒸汽物性计算等模型［13］。

2.3 方程离散求解

守恒方程与本构封闭方程构成了求解热工水力问题的非线性方程组。在方程组的求解过程中，通过数值方法进行求解模型中压力—速度耦合关系。在离散过程中，采用Patankar、Spalding 在交叉网格上速度和压力耦合求解［14］的基础上提出半隐式方法（Semi-Implicit Method for Pressure Linked Equations，SIMPLE），即求解压力耦合方程的半隐方法。

对于一维流体，交叉网格指将速度及压力分别储存在两套不同的网格上。其中，压力储存在控制体中心（节点）处，速度储存在两控制体间的界面上。在计算速度修正值时，无需考虑相邻界面上的速度修正值对该界面上速度修正值的影响，即界面上速度修正值仅由相邻控制体节点间的压力修正值所决定（见图1）。

Fig.1 Storage of pressure and velocity图1 压力与速度存储

3 并联通道系统模型建立与工程验证

3.1 并联双通道系统模型搭建与验证

并联加热双通道实验在高温高压汽水两相热工水力实验装置上开展，装置设计压力17.2 MPa，温度350 ℃，实验质量流量范围100～5 000 kg/（m2·s）［15］。参考实验台架设备，利用NUMAP 平台热工水力模型库中的管道、分支及进出口边界等模型，采用拖拽式建模方式快速建立并联加热双通道系统模型，如图2所示。

Fig.2 Thermal hydraulic model library of NUMAP图2 NUMAP热工水力模型库

搭建完毕后，参考实验台架通过参数面板设置结构参数与初始化参数，建立系统模型，如图3所示。

Fig.3 Model of parallel heating dual-channel system图3 并联加热双通道系统模型

仿真计算工况的开展方式是在系统压力及入口参数保持不变的前提下，阶梯式缓慢提升加热功率，每提升一次功率后保持2 min，直至流动不稳定现象发生时加热功率不变，并记录通道出口的质量含气率，仿真计算结果与实验结果如图4、图5所示。

Fig.4 Calculation error of uniform heating图4 均匀加热计算误差

Fig.5 Calculation error of non-uniform heating图5 非均匀加热计算误差

参考实验工况，对均匀加热和非均匀并联加热双通道流动不稳定性分别进行仿真模拟，选取流动不稳定性发生点通道出口质量含气率的仿真计算结果与实验值进行比较，均匀加热和非均匀加热工况下仿真计算结果与实验结果误差最大值分别为18.87%、15.7%，满足该实验工况下20%的误差要求，表明NUMAP 软件可用模拟分析并联通道流动不稳定性瞬态特性。

3.2 不同加热通道数对并联多通道流动不稳定边界的影响

由于并联多通道系统存在多种核电设备中，其流动不稳定性会导致临界热流密度（CHF）明显下降或设备发生机械振荡，影响核电站正常运行。因此，在核电站的运行过程中，如何避免并联多通道系统流动不稳定现象发生一直是热议的研究课题。本文在验证NUMAP 软件对并联加热双通道流动不稳定性瞬态计算能力的基础上，对并联多通道流动不稳定性边界进行仿真计算。

3.2.1 流动不稳定性边界

Ishii［16］研究发现，流动不稳定性边界可由无量纲相变数Npch和过冷度数Nsub构成。目前许多学者采用Npch-Nsub图来绘制流动不稳定边界和不稳定区域的准确性，效果较好，因此本文也沿用该方法进行实验。无量纲相变数Npch和无量纲过冷度数Nsub为：

式中：Q为总加热功率，单位为W；W为通道进口总流量，单位为kg/s；hfg为汽化潜热，单位为J/kg；Δhsub为进口欠焓，单位为J/kg；vfg为饱和气液两相比容差，单位为m3/kg；vf为饱和液相比容，单位为m3/kg。

3.2.2 多通道流动不稳定性边界

根据上述实验工况进行并联双、三、四通道系统进行流动不稳定性边界的研究，仿真曲线如图6—图8所示。

Fig.6 Flow rate curve for parallel dual-channel system with symmetric uniform heating图6 双通道对称均匀加热流量曲线

Fig.7 Flow rate curve for parallel triple-channel system with symmetric uniform heating图7 三通道对称均匀加热流量曲线

Fig.8 Flow rate curve for parallel quadruple-channel system with symmetric uniform heating图8 四通道对称均匀加热流量曲线

在并联加热多通道系统仿真中，随着加热功率上升，加热通道出口达到饱和后开始出现气液两相流动。随着气液两相流动的含气率逐渐增大，当通道出口含气率增大到一定程度后，并联多通道系统中各通道的流量出现周期性流量脉动，即发生了流动不稳定性，如图6—图8 中局部放大部分所示。

在不同工况下设置加热通道数目分别为2、3、4 时，计算得到流动不稳定发生点过冷度数Nsub和相变数Npch，并在此基础上进行并联多通道系统流动不稳定性边界的比较，如图9 所示。由此可见，流动不稳定性边界相差在±5%以内，仿真计算的结果与王艳林等［15］实验结果一致，即存在多个并联加热通道时，采用并联双通道结构即可获得其流动不稳定性边界。

Fig.9 Flow instability boundaries under different heating channel numbers图9 不同加热通道数目下流动不稳定性边界

4 结语

流动不稳定性是核反应堆热工安全的重要研究课题之一，仿真分析可作为开展流动不稳定性研究的重要手段。本文基于Modelica 语言开发的NUMAP 软件搭建了并联双通道加热模型，仿真模拟结果表明，实验误差满足实验工况下的误差要求，验证了模型对流动不稳定现象的瞬态计算分析能力。

同时，本文使用NUMAP 软件进行并联多通道加热模型的仿真研究，在加热通道数分别为2、3、4 时，计算流动不稳定发生点过冷度数Nsub和相变数Npch的值，发现其流动不稳定性边界相差均在±5%以内，证实了在多个并联加热通道时，采用并联双通道结构可获得其流动不稳定边界。