压力容器-管路-安全阀系统变保真度仿真模型研究

宗超勇，李清野，周威豪，宋学官，颜少华

(1.大连理工大学 机械工程学院，辽宁 大连 116000；2.中广核惠州核电有限公司，广东 惠州 516000)

引言

压力容器-管路-安全阀(以下简称：压力系统)是现代工业，特别是压力、能源系统中的常用配置，被广泛应用于军工、石化、油气以及核电等众多领域[1-3]。其中，压力容器和管路分别主要用于能量的存储和传输，而安全阀(本研究特指弹簧式安全阀)则为整个压力系统提供超压保护。

在针对压力系统的研究中，建立系统级仿真模型，然后基于仿真模型进行压力系统的早期设计、过程分析和优化设计是有效的途径之一。在压力系统的模型研究方面，国内外学者开展了大量的研究，如DARBY[4-6]和ERDÖDI[7]分别采用偏微分方程组联立的方式建立了压力系统的系统级解析模型，并通过实验测试对解析模型的计算精度进行了验证。张仅等[8]利用AMESim软件建立飞机燃油管路的数学模型对不同弹簧刚度的安全阀所造成的压力脉动现象进行了仿真分析，通过对弹簧优化得到比较理想的压力脉动抑制模型。郭崇志等[9]结合ANSYS CFX动网格技术和CEL编译语言对安全阀开启瞬态过程进行了数值模拟，分析了流场参数对安全阀动态特性的影响，并采用自行设计的安全阀测试实验装置对CFD模型进行了验证。杨留等[10]使用全CFD方法完成了高压管路-阀门系统瞬态仿真计算模型的构建，并基于瞬态仿真的结果进行了弹簧式安全阀的参数分析和启闭压差的优化设计。BOUZIDI等[11-12]对工作在低压、亚临界工况下的弹簧式安全阀及其所属的压力系统进行了解析建模和实验分析，并试图据此探究弹簧式安全阀动态不稳定性产生的机理。

综合分析现有的涉及压力系统建模研究的文献发现，目前针对压力系统进行建模的方法主要有解析法和CFD法两种。其中，解析法建模具有高效、通用等优点，但在建模时需要引入大量的简化和假设，这不可避免的会导致模型精度的损失，因此一般仅适用于几何结构简单且流体流动不太复杂的工况。而CFD方法具有计算准确、结果直观等优点，具备处理复杂流动情况的能力，但同时对计算资源提出了更高的要求，在进行精细仿真的过程中往往难以兼顾计算效率。

综合考虑上述解析模型和CFD模型的特点，本研究提出了一种兼顾仿真效率和仿真精度的压力系统变保真度建模方案。其基本原则为：对于流体流动相对简单的压力容器和连接管路，分别使用等效压力点方法(零维)和特征线法(一维)进行容器动态压力的等效和管路压力波动的捕捉，以实现更高的仿真效率；对于具有复杂流动，诸如激波演化、流动分离现象存在的弹簧式安全阀，使用精度更高的CFD方法进行建模，以保证足够的仿真精度。

1 压力系统的基本结构和工作原理

本研究的压力系统特指压力容器-管路-安全阀系统，其基本结构如图1所示，其中压力容器主要用于工作介质的产生或储存，是系统的压力来源；管路与压力容器相连，将压力容器内的压力输送到各执行元件；弹簧式安全阀一般安装在压力容器或管路附近，为系统提供超压保护。当压力系统运行时，系统内各组成部件的动作相互影响，流体介质与阀门固件发生强烈耦合。

图1 压力系统的基本结构

图2 弹簧式安全阀的工作原理图

弹簧式安全阀的主要功能部件包括阀门喷嘴、阀盘、弹簧、调整圈等。阀门喷嘴一般通过螺纹或者法兰与被保护的压力系统连接，是安全阀的压力引入部件，其内部压力一般等于系统压力。当系统压力处于正常范围时，安全阀保持关闭，阀盘在弹簧力的作用下被挤压在阀门喷嘴(阀座)的端部，系统压力仅可抵达阀盘的前部端面。当系统压力增加至安全阀设定压力时，作用在阀盘上的流体力超过弹簧力，阀盘在向上合力的驱动下启跳，系统内流体开始被泄放；随着系统内超压介质的泄放，安全阀内的压力以及作用在阀盘上的流体力同步降低。当阀盘流体力小于弹簧弹力时，阀盘合力向下，阀盘开始回座；当系统压力低至安全阀的回座压力以下时，阀盘与阀座再次接触，阀门密封重新建立，至此完成了安全阀的动作过程。

2 压力系统变保真度模型构建

建模工作主要涉及2个方面，分别为压力系统的元件级模型构建和压力系统的系统级模型耦合。

2.1 压力容器零维模型

压力系统正常工作时，压力容器作为系统边界，不仅遵循理想气体状态方程，还要受介质的流入和流出2个因素的影响，实时压力变化如式(1)所示：

(1)

式中，p——压力容器(管路入口)的实时压力

p0——上一时刻的容器压力

v——系统充压速率(取决于压力容器介质流入的速度)

Δt——瞬态仿真的时间步长

Q——系统入口(容器出口)处的质量流量，代表了压力容器的介质流出速率

R——空气的摩尔气体常数

T——温度

V——压力容器的体积

M——空气的摩尔质量

式(1)既是压力容器的压力变化情况，也代表着连接管路的入口总压，公式等号右边第一项代表了上一个时刻的容器压力，第二项v·Δt代表了单位时间内容器压力的补充效果，第三项代表了由于系统泄放导致的容器压力降低效果。

2.2 连接管道MOC一维模型

MOC方法是一种用于近似求解双曲型偏微分方程的近似解法，其基本求解思路是通过求解空间的转换将适用于全流场的偏微分控制方程转化为仅在特征线上成立的常微分方程，然后使用常规方法对常微分方程进行求解。在本研究的压力系统中，使用MOC方法对连接管路进行一维模型构建的基本思路如下文所述。

同所有的流体力学一样，连接管路内的流场同样遵循流体力学基本控制方程，即连续性方程和动量方程(对于可压缩介质还需补充能量方程和理想气体状态方程)。将连接管路内的流体假设为一维单向流动，2个控制方程的具体形式分别如式(2)、式(3)所示：

(2)

(3)

f——Darcy-Weisbach摩阻系数

D——管道直径

a——管路中的声速

q——流量

g——重力加速度

A——管路截面积

p——管路压力

t——时间

x——管路位置

式(2)、式(3)均为偏微分方程，为了便于计算，对其等号左边项目进行线性组合，得到：

(4)

其中，χ为未知算子，进一步整理后得：

(5)

其中，流量q和压力p均为管路位置x和时间t的函数，其全导数分别为：

(6)

(7)

全导数式(6)、式(7)为控制方程式(5)由偏微分形式向常微分形式转换提供了可能，此时需要定义合适的χ值，令：

(8)

结合式(6)、式(7)，式(5)可重写为：

(9)

(10)

其中，式(9)、式(10)分别仅在左右特征线上成立，如式(11)、式(12)以及图3所示：

(11)

(12)

至此，完成了管路流体控制方程从偏微分方程式(5)向常微分方程组式(9)、式(10)的转换。此方程转换可有效降低方程的求解难度，但同时也缩小了控制方程的适用性。因为式(5)在管路的空间-时间空间内任何一点均适用，而简化后的方程仅沿着左特征线(图3中点划线)和右特征线(图3中虚线)才分别成立。

图3 连接管路MOC模型的基本结构

图3是管路系统的x-t空间，在初始状态下(时间t=0)，点A和B的流场参数，包括当地压力和流量通过初始化或者依据稳态计算结果可以直接给定，点C为过点A和B的2个特征线的交点，即相容方程式(11)、式(12)在C点是同时成立的。因此，通过方程联立即可求解出C点位置在Δt时刻的流场信息，再通过对x-t空间进行上述过程的遍历即可实现不同时刻不同管路位置处的流场信息的求解和更新。

2.3 弹簧式安全阀二维CFD模型

如前文所述，安全阀动作时，其阀门附近会出现激波迁移、流动分离等复杂流动状态。由于这些流动的高度非线性和随机性，使用解析法或特征线法均难以实现精确模拟。为了保证压力系统的整体仿真精度，本研究使用了CFD方法对弹簧式安全阀进行了动态模型的构建，其计算域的划分和计算网格如图4所示。

图4为弹簧式安全阀的二维CFD网格模型，为了便于网格密度和质量控制，本研究对弹簧式安全阀的计算域进行了分解，将其拆分为3个子域，分别为阀内子域、喷嘴子域和外场子域，子域之间以Interface接口进行模型连接和数据交换，如图4所示。在3个子域中，喷嘴子域最为关键，几乎所有的复杂流动均在此区域内发生，因而也具有最高的网格密度。另外，弹簧式安全阀的阀盘运动受阀盘流体力和弹簧力共同控制，当作用在阀盘上的合力不为0时，阀盘将开始运动，这就要求CFD模型中的阀盘边界同步更新。为了实现在瞬态CFD仿真中对阀盘的运动进行控制，本研究将喷嘴子域中阀盘和喷嘴部分使用Interface交界面分开，具体方式如图4中的放大图框所示。

图4 弹簧式安全阀二维轴对称CFD模型

其中，图4中的A区域为初始状态下的网格分布，由于Fluent不能直接求解分离的流场，为了保持求解域的连续，本研究将安全阀全开度的1%(0.036 mm)作为关闭时的阀门状态。

阀盘边界的更新是通过控制整个阀盘部分网格的运动实现的，具体方式为：激活Fluent软件中的Dynamic Mesh模块和Layering选项，然后将整个阀盘网格部分设置为运动刚体(运动速度由UDF中的DEFINE_CG_MOTION宏命令控制)，并把左右两个流场边界(Interface_01和Wall_01)设置为静止墙面，这样即可实现边界网格随着阀盘运动的实时更新，更新后某一开度下的网格如图4中的B区域所示。

网格质量是保证CFD仿真精度和计算收敛性的重要因素，为了获取最优的网格模型，本研究以弹簧式安全阀的阀盘升力为评价指标对弹簧式安全阀的计算网格进行了无关性测试，测试结果如图5所示。

图5 弹簧式安全阀网格无关性验证结果

如图5所示，本研究使用了Coarse、Fine和Very Fine 3个密度等级的网格模型，配合4个两方程湍流模型对全开阀门状态下(3.6 mm阀门开度，工作压力0.3 MPa)弹簧式安全阀的阀盘升力进行了计算，并将仿真结果与已有的实验结果进行了对比。结果显示，不同密度的网格模型确实对阀盘升力有重要的影响，当网格密度从Coarse增加到Fine等级时，阀盘升力的结果变化较大，但当网格密度大于Fine等级时，阀盘升力几乎不变(变化范围在实验值的1%以内)，这说明当网格的密度高于Fine等级时，CFD的仿真结果即与网格密度无关。综合考虑计算资源和效率因素，本研究按照Fine等级的网格生成方案对弹簧式安全阀进行了模型的构建。另外，网格无关性验证结果还表明，不同湍流模型在阀盘流体力计算方面的能力是不同的，其中SSTk-ω模型的计算结果与实验测试值之间的相对偏差最小，是用于阀盘流体升力预测的最优湍流模型。

2.4 系统级动态仿真耦合模型

2.1～2.3小节分别完成了压力容器、连接管路以及弹簧式安全阀的元件级模型构建。本节主要基于上述元件级模型构建压力系统的系统级变保真度耦合模型，建模思路如图6所示。

图6 压力系统多尺度变保真度耦合建模方案

在压力系统的仿真中，不同模型在每个时间步中同步求解，其中，压力容器零维模型计算得出的瞬变压力和阀门CFD模型计算得出的瞬变质量流量分别作为管路特征线模型的左侧和右侧边界条件，边界条件的加载以及不同模型间的数据交换通过interface实现。另外，弹簧式安全阀的阀盘位置、速度和阀盘升力以及容器的压力等关键物理量都是随时间动态变化的。为了捕捉这些变化，本研究使用UDF方法对这些物理量进行实时地记录和更新。

3 压力系统变保真度仿真模型的验证

为了验证本研究构建的压力系统变保真度仿真模型在压力系统动态特性计算方面的能力，将仿真模型的计算结果与实验结果进行了对比。由于MOC建模方法的主要作用是捕捉管路中的压力波动，而在弹簧式安全阀的动作过程中，关阀动作导致的管路压力波动现象相比开阀动作更加明显。因此，本部分的验证主要针对弹簧式安全阀的关闭过程进行，模型验证的结果如表1和图7～图10所示。

表1 实验与仿真结果对比

图7～图10分别代表了弹簧式安全阀回座过程中的阀盘运动、管路末端(阀前)压力波动、容器压力变化和管路压力波动的频率；表1为实验和仿真结果对比的量化结果。由图7～图10和表1可知，与实验结果相比，系统级多尺度变保真度仿真模型能够对弹簧式安全阀的关闭过程，以及由关阀动作引起的压力容器和管路内的压力波动进行较为准确的仿真。在管路压力波动(阀前压力)方面，实验和仿真计算得出的阀前压力波的幅值和频率的相对偏差均小于4%，这说明本研究使用的压力容器等效压力点方法、MOC一维管路建模方法、弹簧式安全阀二维CFD建模方法在压力系统的动态特性计算方面具有较高的精度，可以实现保证仿真精度的前提下提高压力系统的仿真效率。

图7 阀盘位移实验与仿真结果对比

图8 阀前压力实验与仿真结果对比

图9 容器压力实验与仿真结果对比

图10 阀前压力-频率实验与仿真结果对比

4 结论

本研究针对由压力容器-管路-安全阀组成的压力系统开展了系统级变保真度仿真模型的研究。分别使用了等效压力点技术、MOC理论和CFD方法对压力容器、连接管路和弹簧式安全阀进行了元件级模型的构建和系统级模型的耦合。为了验证系统级模型的仿真精度，将弹簧式安全阀回座阶段的仿真结果与实验结果进行了对比，结果显示本研究构建的系统级变保真度仿真模型能够对弹簧式安全阀的阀盘位移、管路压力波动以及容器压力变化准确的计算，证明了本研究提出的建模方法以及所建模型的准确性和适用性。

本研究提出的压力容器-管路-安全阀系统级建模方法在保证仿真精度的同时，可以大幅提高仿真计算的效率，这为液压与气动回路中具有类似结构的系统级建模和仿真研究提供了参考。