摘要：采用ANSYS有限元软件计算得到减压阀膜片的反力时程曲线和运动部件的动态特性，在此基础上将结构简化为两个单自由度质量弹簧阻尼动力学模型，采用Newmark算法计算。流体控制方程为三维非定常积分形式的ALE(Arbitrary Lagrangian Eulerian)方程，采用基于弹簧近似的动网格有限体积格式求解，应用了新的离散几何守恒律和高精度界面算法；同时采用虚拟网格通气技术实现阀门部件运动过程所引起的拓扑变化。计算表明，在较宽的上游增压速率范围内减压阀出口压力存在振荡，均值接近按照静态性能设计的理论值。通过对计算流场进行分析，确定了造成两种开启故障的主要机理，修改模型参数可以排除异常。
关键词：流固耦合； 减压阀； 非结构动网格； 压力振荡； 动态特性
1 概 述
减压阀通过其内部构件调节管路系统的流量和压力，使上游的高压气体在下游出口处维持相对稳定的低压。传统设计流程中主要考虑静态性能，目前已有比较成熟的计算方法能够预测减压阀稳定工作状态下的出口压力和流量特性[1-2]。由于分析手段有限，通常情况下只能依靠样品实验方法进行动态特性研究，由于阀门内部结构复杂、体积小，很难进行数据测量工作，动态试验往往只能给出总体性能是否满足设计指标的定性结论[3]。
近年来有许多学者采用数值仿真手段开展阀门研究，根据前期调研情况看，对液体阀门管路系统中产生的水锤、空化、气蚀等现象进行仿真研究相对较多；而气体管路系统中也会出现所谓的“气柱振荡”现象，在上世纪70年代以后逐渐引起学术界关注，利用这种共振产生激波进行能量传导可用于制冷[4-6]。减压阀的减压功能主要依靠调整流道最小截面产生壅塞来实现，由于气体可压缩性强、声速较低，在阀门内部复杂的结构中流动时常会形成激波等超声速流动现象。文献[7-8]通过数值模拟得到的结论可以看出，由于阀门的限流作用，内部流场中形成的激波导致了气动载荷剧烈变化，这是美国Stennes航天中心(SSC)在火箭发动机试车过程中稳压阀出现故障的主要原因。文献[9-10]进行固体火箭发动机调压阀的数值模拟时也发现存在马赫数高达4的局部区域。由于流场内存在激波，常用的商业软件中SIMPLE或PISO等基于不可压缩流动建立的算法不再适用，需要在弱解理论指导下发展能够捕捉激波的新方法。除此之外，阀门内部复杂结构网格生成、启闭引起的空间拓扑变化、描述运动部件的网格技术、流固耦合的界面算法、计算结果的验证等问题也是数值模拟应用于阀门动态性能分析遇到的挑战。
图1是一种膜片式减压阀的结构示意图。低压腔与出口管路相连接，其下方密封膜片盒内的主弹簧K1使膜片K2向上变形，推动阀芯M02向上运动，在阀芯和阀口(下限位K5)之间形成流通面，高压腔和低压腔之间的压差在该处形成壅塞，实现降压功能。阀芯和阀座顶盖之间的副弹簧K3和上限位K6控制阀芯运动规律，因此流通面大小受很多因素所决定。如果出口压力超过预定值，作用在膜片的气动力破坏弹簧K1，K2和K3之间受力平衡，推动膜片向下运动，带动阀芯使缝隙变窄，减少进入低压腔的气体流量，致使压力降低；反之，如果出口压力小于预定值，气动力推动阀芯使缝隙变宽，更多气体进入低压腔导致压力升高，这是实现其稳压的原理。为了避免高压冲击载荷引起的下游波动，采用阀芯与膜片中心落座可以分离的结构设计，如果膜片向下运动位移大于下限位K5，膜片脱离阀芯，副弹簧K3使得阀芯和阀口完全贴合，阀门完全关闭，该流通面的流量为0。同时，为了避免膜片回弹过量，膜片盒内设置了下限位K4。阀芯和上限位K6之间的空腔称为卸荷腔，阀芯和阀座之间缝隙使得高压气体进入卸荷腔进而影响阀芯受到的气动载荷，为此阀芯中间开孔与低压腔相通，即使阀门完全关闭也存在一定流量。
图1 减压阀结构示意图
Fig.1 Schematic of pressure relief valve
基于“虚拟现实”的优势，本文通过数值模拟的方法预测这种减压阀开启过程的动态特性，分析设计参数的影响机理，进而为工艺流程控制和产品检验提供依据。
2 结构有限元模拟
尽管在结构分析中常用ANSYS、FLUENT等商用软件，这些软件中也包含流体单元的计算模型，但这些模型都是针对不可压流动方程建立，不适用于气体减压阀。由于ANSYS软件和流场求解器之间存在耦合问题，本文把减压阀内部运动简化为如下两个单自由度的质量弹簧模型。
如图1所示，在阀芯与膜片的分离面处将阀门结构分为上下两部分。膜片系统的位移变量采用x1表示，方向向上为正，运动方程如下
F+F1-Fm(x1)
(1)
式中 M01为膜片质量，K1为主弹簧刚度，M1为质量，L1为预压缩量，下限位行程L4=1．574 1 mm。膜片运动中，调整块与下限位K4碰撞过程采用“接触→压缩→回弹→脱离”非线性弹簧模型，弹簧刚度K4表示为
(2)
在阀芯与上限位K6接触的全开状态作为初始位置x1=0，运动到x=-1.8 mm位置时下限位K4发挥作用。F是阀芯和膜片相互间作用力；F1是作用在膜片上的气动力，可以根据低压腔内流场压力积分得到，在与上部接触时去除阀芯对应的面积；Fm(x1)是膜片变形产生的反力，由于膜片厚度很薄，当位移相对其厚度较大时，需要考虑位移大变形效应，膜片内部应力很大，可能超出其屈服限，因此，膜片变形提供的力和位移之间的关系是非线性的。对于图1中刚度K2进行预处理建模，在流固耦合模拟前先用ANSYS软件对主阀膜片进行有限元分析，得到位移变形与反力的关系Fm(x1)供本文计算时使用。
阀芯系统的位移变量采用x2表示，运动方程如下
(3)
式中 M01为阀芯质量，K3为副弹簧刚度，M3为质量，L3为预压缩量。阀芯碰撞过程用非线性弹簧K5和K6近似，下限位行程为L5=0.874 1 mm，初始位置x2=0是全开状态，上限位行程为L3=0；f是阀芯和阀座之间涨圈的运动摩擦力。如果阀芯脱离膜片中心落座块的情况下，相互间作用力F=0；如果阀芯接触膜片一起运动，以上2个动力学模型合并成1个方程，这时F为内力，相互抵消；F2是作用在阀芯上的气动力，受低压腔、高压腔和卸荷腔内压力影响。
由于模型中受力均是非保守力，构成非保守系统，采用Newmark算法计算。在本文数值模拟过程中，由于三维流动的复杂性，每个迭代步内流场计算耗时远大于固体位移场计算耗时。
3 流体数值方法
气体的粘性比液体小，在气体阀门中，局部空间会形成超声速流动，由激波和膨胀波引起的载荷变化远远大于粘性的贡献，表现为NS方程中时间导数项比粘性项的影响要大。为了提高计算效率，从三维Euler方程出发进行数值模拟，本文的算例也验证了这种简化方式是可行的。为了描述网格变形，必须把流体控制方程变换为ALE形式
(4)
式中 Ω为控制体，∂Ω为控制体边界；n为控制体边界外法向单位向量n=(nx,ny,nz,)；守恒变量Q和对流项F(Q)为
式中
为流体相对于网格的速度矢量，V=(u,v,w)为流体的速度矢量
为网格运动速度矢量，at=xtnx+ytny+ztnz为网格运动的法向速度。
采用守恒变量存贮在单元中心的格心格式进行空间离散，内部复杂的空间结构均采用四面体非结构网格，在进出口和低压腔布置测点如图2所示。
图2 模型网格
Fig.2 The grid of model
假设物理量在单元内线性分布，可以构造出空间二阶精度的格式。在跨声速和超声速绕流计算中为了捕捉激波，依据Riemann问题提出了多种守恒型的通量计算格式。本文采用Van Leer矢通量分裂格式计算通量，时间离散采用二阶时间精度的四步Runge-Kutta方法，采用弹簧模型处理网格变形，在流场边界运动以后带动内部网格节点。采用这种动网格技术除了需要发展可靠性高的网格变形算法，还要注意采用离散几何守恒律克服网格变形给流场计算引入的误差。关于时空二阶的有限体积法推导过程、抑制非物理振荡的限制器、弹簧模型网格变形方法、虚拟挡板通气技术和新的离散几何守恒律的具体细节见参考文献[11]。
虚拟挡板通气技术的基本原理[11-12]如图3所示，生成初始网格时在接触/分离位置或小尺度缝隙处设计一层网格，一般情况下这层网格是全场最小尺度的网格，根据运动物体在网格内部设置虚拟档板，即图3中阴影部分作壁面处理，其他部分作通气处理。当运动物体的位移等于虚拟档板的初始网格以后，采用正常的变形网格拉伸来描述宽度增加，当拉伸距离较大导致网格质量变差，采用局部重构进行网格更新。阀门尺度为100 mm量级，最小网格尺度大约在0.5 mm，运动行程有限，一般进行2～4次重构可以满足要求。这种新技术可以计算包括接近零的无限小距离，模拟从接触到逐步远离的真实物理过程，同时也不会因网格尺度太小而影响计算时间推进步长。
本文在3个地方设置了虚拟挡板：(1)阀芯和下限位K5之间缝隙设置模拟开度小于0.404 1 mm，直到完全关闭状态的开度变化；(2)阀芯和阀座的涨圈结构存在缝隙使得高压气体进入阀芯上方的卸荷腔，影响阀芯受到的气动载荷，卸荷腔与低压腔相通，也影响出口压力，因此，在缝隙进口位置设置网格宽度为0.45 mm的网格来模拟缝隙宽度为0.009 mm的微小流量；(3)高压气流在进入高压腔之前经过了一个多孔网状结构的过滤器，根据前期实验数据有明显压降，也需要设置虚拟挡板。
图3 虚拟档板示意图
Fig.3 Sketch of virtual baffle faces
4 流固耦合界面算法和程序验证
流固耦合系统的求解通常有两种方法：全耦合方法(Monolithic Method)和分解方法(Partitioned Method)。全耦合算法是指流体控制方程和结构方程作为一个系统方程同时求解，公式推导复杂，建立高精度算法困难，求解过程涉及到矩阵运算，难以应用于实际工程中的多自由度结构系统，目前使用很少。分解算法是把耦合问题分解为流场和结构两部分，分别进行求解，流场和结构之间的耦合通过流固耦合界面上的信息交换来实现，又称为交错迭代算法(The Conventional Serial Staggered —CSS)。相对于全耦合算法，分解算法具有以下优点：(1)可以充分利用已有的CFD和CSD计算方法和程序；(2)可以减少程序开发的难度；(3)可以保持程序的模块化。所以，自从这一算法提出以来得到就受到重视。为提高精度和计算效率，在流体和结构方程交错求解中又引入内迭代，内迭代达到收敛后全系统统一时间。这样一来，根据是否存在内迭代，进一步把CSS分为紧耦合算法(Tightly-Coupled Method)和松耦合算法(Loosely-Coupled Method)。实际应用表明，紧耦合算法在有些情况下达到收敛标准的内迭代计算量很大，计算效率比时间小步长推进的松耦合算法还低。在文献[13]中结合新的几何守恒律算法，构建了高精度的松耦合算法。
如图4示意，高度H=1、长度为L=10的二维空间均匀分布无量纲密度ρ0=1、压力p0=1、声速a0=1的静止气体，右侧连接质量m=1的活塞，活塞和固定点之间有刚度系数k=1的弹簧。根据非定常一维均熵流理论可以得到活塞运动速度与压力的关系式
(5)
图4 验证模型及其结果比较
Fig.4 Comparison of displacement for piston problem
采用Runge-Kutta方法求解以上常微分方程可得到活塞运动特性，用来验证前面介绍的流固耦合计算方法的有效性。采用2 844个网格点、5 260个单元的三角形网格离散初始空间，动网格模拟流场，在t=3时两种方法计算得到的活塞速度和位移比较如图4所示，计算结果非常一致，表明本文建立的流固耦合计算程序具有良好的计算精度，为下面的应用研究打下很好的基础。
5 数值计算与结果分析
减压阀上游通过管路连接高压气源，下游管路有孔板限流器。工作介质为常温条件下的空气。为了避免减压阀内部压力波对进出口边界影响，在三维计算域的入口和出口处连接足够长的一维等截面管路模型，一维模型与三维模型之间通过延拓两排虚拟网格点来传递流场物理量，一维模型的自由端分别定义为无反射入口和出口边界条件。入口长度为10 m，出口长度为2.5 m。
按照静态性能设计模型，主弹簧的预压力为1 050 N，出口压力1.7 MPa为理论值。标准工况的上游开启总压为21 MPa，初始流场为标准大气参数0.1 MPa。在一维模型的入口处总压在给定时间内从0.1 MPa线性增加到21 MPa，开启时间大约在10～100 ms之间。计算得到的压力随着时间的变化规律如图5所示。尽管低压腔测点位置不同，曲线差异很小，表明在低压腔内主要进行亚声速流动。算例结果基本围绕比理论值略高的均值波动，考虑到黏性影响，计算值高于理论值也是合理的。
图5 开启压力为21 MPa时不同监测点压力变化时程
Fig.5 Pressure curves of observation points at the case of 21 MPa upstream pressure
计算过程中在过滤器虚拟挡板两侧形成压力梯度，涨圈结构也有明显的压力差存在，说明本文采用虚拟挡板成功模拟滤网结构的阻挡作用和高压气体进入卸荷腔的流动过程。如果主弹簧预压力为678 N，出口压力调整到1.2 MPa为理论值，数值模拟的低压腔测点压力规律与图5类似，也是围绕理论值波动。以上这些论据验证本文计算是合理有效的。
减压器在正常工作过程中，由于流固耦合相互干扰，膜片脱离阀芯后向下运动，与限位器碰撞后反向运动，逐渐在限位器和原点之间发生高频振荡；阀芯振幅均值较小，除了接触膜片一起向上运动，在脱离膜片后也有高频小幅振动。
同时，在数值模拟过程中出现了减压器失效的故障状态，对多种工况计算结果进行总结分析，主要的故障状态可以分为两种，第1种如图6(a)所示，在工作一段时间以后，膜片的运动位移逐渐增大直至发散，这样可能导致膜片损坏；第2种如图6(b)所示，减压阀可以稳定工作，但是这种工况下的出口压力低于理论值。
图6 故障状态下阀芯与膜片位移历程
Fig.6 Displacement curves of valve core and diaphragm
图6(a)为第1种故障状态的计算结果，当开启压力为15 MPa时，减压阀阀芯和膜片以小幅振动的形式工作了一段时间后便振动发散，低压腔压力均值从1.7 MPa变化到7.5 MPa，如图7所示。由这个故障工况的结果推测，采用原来设计的低压腔外形，当近距离布置出口孔板时，低压腔连通孔板前管路的容积相对较小，而高压腔压力很高，低压腔压力对阀门开度过于敏感，每一次波动都会导致低压腔压力出现一次跃升，如果压力的振动与结构的运动产生共振就会导致颤振，最终必然导致结构损坏。
图7 压力变化时程曲线
Fig.7 Pressure curves of low-pressure cavity
分析图6(b)的第2种故障原因，由于阀芯有下限位K5，因此在图中坐标原点处以下是膜片脱离阀芯后的单独运动状态，原点以上可以是阀芯单独运动也可以包含膜片的重合运动。从图中可以看出，膜片与限位器K5碰撞后继续向上运动，并与阀芯顶杆碰撞，碰撞时两者的动量方向相反但大小接近，因此碰撞后两者并保持接触，直到关闭状态并保持速度不变。低压腔压力最终稳定在一定值，但不是设计值。由于出口孔板一直在向外界排气，低压腔压力之所以保持稳定，是因为涨圈结构的漏气进入卸荷腔再经由阀芯中心孔进入低压腔与出口排气达到了平衡。此时弹性元件的受力也达到平衡，图8为稳定工作时的马赫数云图。这种非正常工作状态是由于建模时设定的涨圈漏气量过大造成的，说明涨圈漏气量过大会导致减压阀失效。
图8 开启压力为2 MPa时阀门结构的马赫数云图
Fig.8 Mach contours in failure state at the case of 2 MPa upstream pressure
上文给出的模拟结果都是近距离布置出口孔板的情况，有些工况出现故障状态。增大低压腔容积虽有助于改善弹性元件的工作状态，但效果并不明显。对原始模型进行改进，模拟远距离布置出口孔板的工况，设定的孔板放置位置距离低压腔2 m，但对于2 m长的出口管路真实建模将会导致计算域不同方向的尺度相差太大，因此本文将2 m长的出口管路等效成一个体积相同的容腔，置于减压阀出口附近，其模型及网格如图9所示。
图9 新模型网格
Fig.9 The new grid of model
采用这种新模型重新进行计算，除了计算外形不同外，计算模型中还重新设定了膜片硬芯与下限位的碰撞能量损失，由原来的无损失改为损失20%；另外卸荷腔虚拟挡板的通气量也根据设计值重新进行了设定。图10给出了开启压力为15 MPa时采用新模型进行计算得到的稳定工作状态下的马赫数云图。
图10 开启压力为15 MPa新模型稳定工作时的马赫数云图
Fig.10 Mach contours within the new model at the case of 2 MPa upstream pressure
分别针对两种故障工况，采用新模型重新进行模拟，计算结果得到明显改善，说明修改模型以后，可以抑制导致故障状态的振动发生。
从图10可以看出，更改卸荷腔虚拟挡板通气量后，顶杆中心孔内的流速非常小，高压气体几乎不通过涨圈结构进入卸荷腔；从图11可以看出，开启压力为15 MPa和8 MPa的工况在稳定工作时低压腔压力没有振动而是接近理论值1.7 MPa左右，开启压力为21 MPa的工况低压腔压力呈现出有限幅值的波动。低压腔压力趋于稳定的工作状态，将不会再出现颤振现象，有效地防止了减压器的结构出现损坏。
图11 压力变化时程曲线
Fig.11 Pressure curves of low-pressure cavity
6 结 论
以上研究表明：
(1)在两种主弹簧预压力的大部分算例中，计算得到的出口压力的均值接近按照静态性能设计出的理论值，基于非结构动网格的流固耦合计算方法有效开展减压阀的动态特性研究，采用虚拟网格技术实现了减压阀从完全关闭到开启的真实工作过程。
(2)在进口气体的冲击载荷作用下，这款减压阀正常工作的出口压力表现为周期振荡，而且对于进口气体的不同增压速率有很好的适应性。
(3)减压阀会出现气动力与动力学系统相互干扰形成的颤振现象，定性分析这种流固耦合现象是导致颤振的主要机理，可采取相应措施可以消除颤振。
(4)由于涨圈结构的设计问题，导致涨圈漏气量过大造成减压阀失效，修改模型参数可以消除这种故障，在工艺控制时需要关注缝隙宽度的检测。
