含相变压力振荡改进算法及其验证

刘 培 启，李 想，冯 明 宇，于 洋，胡 大 鹏

( 大连理工大学 化工学院, 辽宁 大连 116024 )

0 引 言

波转子技术通过气体自身压力能来实现能量的传递．该技术的核心部件为波转子，是一种由若干双开口方形直管，即压力振荡管(pressure oscillation tube)组成的转毂[1]．目前波转子技术已广泛应用于增压和制冷领域[2]，该技术在制冷领域的应用便是压力振荡制冷技术．

压力振荡制冷技术依靠高压气体与低压气体间产生的压力波完成能量的直接交换．低压气体经激波压缩为高温气体，而高压气体经膨胀波膨胀为低温气体．目前，压力振荡制冷技术已经广泛应用在天然气脱水、低温风洞等冷凝分离工艺领域．起初，压力振荡制冷技术依靠单开口压力振荡管实现制冷．然而当进气含有水蒸气等可凝组分时，压力振荡管内凝结形成的液体无法彻底排出并不断地积聚在振荡管的封闭端，最终导致制冷性能下降[3]．Dai等提出用波转子取代单开口压力振荡管应用在压力振荡制冷技术中[4]．由于波转子内压力振荡管两端开口，管内的液滴不会形成堆积而是通过任一管口排出，这样既提高了压力振荡制冷技术的带液能力，也减小了设备的占地面积[5]．在此基础上，赵家权对压力振荡管内高压气体的膨胀功进行回收，进一步提高了压力振荡制冷技术的制冷效率[6]．

在压力振荡制冷的过程中，如果进气是含湿天然气等含有可凝组分的混合气，压力振荡管内会存在包含凝结与蒸发在内的相变行为．而可凝组分的凝结过程会产生压力波，对压力振荡管内原有波系造成干扰[7]．由于压力振荡制冷技术依靠压力波的运动实现制冷且压力振荡管与各端口之间的匹配往往依据忽略相变的理想波系进行设计，凝结过程会使得内部实际波系偏离理想波系，从而使得压力振荡制冷的效率无法达到最大值．近几年来，科研工作者开始重视压力振荡管内的相变过程．赵家权等忽略蒸发现象，建立了理想气体凝结数值模型，并模拟了压力振荡管内的凝结过程，探究了凝结对于流场的影响[8]．徐思远搭建了压力振荡管相变可视化实验平台，并通过高速摄影捕捉到管内的凝结行为[9]．吴柯含将蒸发模型与凝结模型相结合，对压力振荡管内理想气体的凝结与蒸发过程进行了数值分析[10]．然而作为可凝组分的水是极性物质，水分子相互缔合形成氢键，使得其实际热力学性质与并没有考虑到氢键作用的理想气体偏差较大[11]．因此，利用理想气体状态方程无法准确地对压力振荡管内水的相变行为进行数值计算．

本文利用CPA气体状态方程对压力振荡管内凝结蒸发数值模型进行改进，并通过实验对改进后的数值模型进行验证．在此基础上，本文分别利用改进后的数值模型与理想气体模型对湿空气压力振荡制冷过程进行数值计算，并对计算结果进行对比总结．

1 压力振荡制冷技术

压力振荡制冷技术的工作过程可通过将波转子与各端口展开得到的二维波系图进行描述，如图1所示．当压力振荡管与高压入口接通时，由于二者之间存在压差，激波S1会在二者交界处形成并向管内运动，压缩管内气体．同时相对地，在高压入口与管内气体之间也会形成一系列膨胀波使得高压进气温度和压力降低．当向上移动的压力振荡管与高温出口接通时，被激波加热的气体通过高温出口排出，同时激波到达压力振荡管的右端并反射形成一系列左向移动的膨胀波E1进一步膨胀管内气体．随后，在压力振荡管与高压入口错开时，由于管内气体的惯性，振荡管的左侧会形成一系列右向膨胀波E2再一次膨胀管内气体．管内气体经膨胀波多次膨胀后温度和压力下降至某一较低值．从高温出口流出的气体经过换热至室温后从中压入口流回压力振荡管，管内低温低压气体会在压力振荡管两侧压差的推动下由低温出口排出．一个周期的制冷过程结束．在整个制冷过程中高压进气与管内原有气体之间会形成较为明显且随着时间运动的交界面．两股气体之间存在通过交界面换热的现象，但是不存在较为明显的混合现象．

图1 压力振荡制冷波图

制冷效率，即等熵效率η，是用来衡量压力振荡制冷技术制冷性能的重要参数，其表达式为

(1)

式中：T1为高压进气温度，K；T2为低温排气温度，K；p1为高压入口压力，kPa；p2为低温出口压力，kPa；γ为绝热系数．

当压力振荡制冷技术应用在天然气脱水等冷凝分离技术中时，液化率也是衡量其工作性能的重要参数．在本文中，液化率是低温出口液态水的质量流量与高压进气中水蒸气的质量流量的比值．液化率利用进出口水蒸气质量分数可表达为

(2)

式中：α为液化率，w1为高压入口的水蒸气质量分数，w2为低温出口的水蒸气质量分数．

2 数值模拟及流场分析

含湿压力振荡制冷实验平台以及实验流程如图2、3所示．经压缩机排出的高压空气首先流向两条管路．在一条管路中，高压空气在经过雾化器与过滤器后获得饱和水蒸气，而在另一条管路中，空气经吸收式干燥机后成为干空气．通过控制两条管路的气体流量并将二者混合，就可以控制压力振荡管进气的相对湿度．在压力振荡管进气端安装湿度仪可以对其湿度进行监测．压力振荡管的尺寸参数如表1所示．

图2 实验平台

图3 实验流程

表1 压力振荡管结构参数

3 数值模型

3.1 凝结蒸发模型

在建立凝结蒸发模型前，本文作出如下假设：液滴与气体之间无滑移，液滴体积为零，液滴之间无相互作用，液滴内部温度分布均匀，流场中没有离子、灰尘等外界颗粒．

凝结过程分为成核与液滴生长两部分．由于流场已被假定没有外界颗粒，压力振荡管内的成核过程为均相成核，即蒸汽分子自身聚集形成凝结核．本文采用由经典成核模型发展而来，经由Feder修正后的Frenkel均相成核率公式[12]：

(3)

式中：J为成核率，1/(kg·s)；qc为凝结系数；ρv为可凝气体密度，kg/m3；ρl为液滴密度，kg/m3；Mm为可凝组分的摩尔质量，kg/mol；φ为Feder绝热修正系数；r*为凝结临界半径，m；σ为液滴表面张力，N/m；k为玻尔兹曼常数，1.38×10-23J/K．

对于液滴生长模型，采用经实验验证有效的Hertz-Knudsen方程[13]：

(4)

式中：r为液滴半径，m；t为时间，s；p为压力，kPa；T为温度，K；R为气体常数，8.314 J/(mol·K)．

对于滴液蒸发模型，经实验验证，采用逆向描述Hertz-Knudsen方程的方法可有效地描述液滴的蒸发行为[14]，因此本文采用同种方式对蒸发行为进行描述．

3.2 数值模型改进

为解决理想气体相变模型计算准确性不足的问题，本文采用可以准确描述极性物质性质的CPA气体状态方程对数值模型进行改进．CPA气体状态方程由SRK方程作为物理项和Wertheim 理论提出的缔合项组成[15]，可以较为准确地描述极性气体的逸度、比热容等参数．

(5)

式中：V为摩尔体积，L/mol；a为SRK方程物理项能量参数，kPa·L2/mol2；b为SRK方程物理项的体积参数，L/mol；g为径向分布函数；xi为组分i的摩尔分数；xAi为分子i上没有与其他端位缔合的A端的摩尔分数．

如果计算介质为混合气，需要采用混合规则对混合物的物理参数进行定义．本文采用经典的范德华混合规则[16]．

(1)密度计算

在理想气体相变模型中，气体密度通过理想气体状态方程求得，液体密度由于可压缩性较差，则默认不变．在本文的改进模型中，气体与液体的密度通过CPA气体状态方程在压力p、温度T下求得．

(2)相平衡计算

在相变过程中需要对可凝组分的气液相平衡进行计算．理想气体相变模型并没有考虑到可凝气的相平衡．由于处于相平衡状态的组分在各相中的逸度相等，在改进模型中，逸度将作为判断相平衡的基本参数．CPA气体状态方程中各组分的逸度系数可由剩余亥姆霍兹能求出[17]，即

(6)

式中：Ar为剩余亥姆霍兹能，J；φi为组分i的逸度系数；ni为组分i的物质的量，mol．

CPA气体状态方程的剩余亥姆霍兹能由SRK物理项和缔合项两部分组成[18]，即

(7)

(8)

(9)

(3)其他热力学参数计算

由于数值模拟中涉及能量的传递，需要对介质的热力学参数进行改进．这些热力学参数包括定压比热容cp、定容比热容cV、绝热系数γ以及声速c．

CPA气体的cp可通过下式进行计算：

cp=(∂H/∂T)p

(10)

(11)

(12)

计算得出cp后，cV可通过下式进行计算：

(13)

据此便可求出声速c与绝热系数γ[19-20]：

(14)

γ=cp/cV

(15)

3.3 几何模型及求解器设置

压力振荡管的三维模型可以准确地模拟其内部的流场，存在的不足就是对电脑性能要求较高且计算时间较长．由于压力振荡管的长径比大于10，且其转速较低，离心力对流场的影响可以被忽略．除此之外，在过去的研究中将压力振荡管沿周向展开得到的二维模型已经被证实可以准确地描述内部流场[21]．因此，本文采用与实验设备尺寸相同的二维模型进行数值计算．计算模型网格划分如图4所示．压力振荡管与各个端口之间采用滑移网格，振荡管的上下边界设为周期性边界．高压入口设置为压力入口，高温、低温出口设置为压力出口．中压入口设置为质量流量入口且通过自定义函数(UDF)将其质量流量设置为与高温出口相等．

图4 几何模型与网格划分

气波制冷机的内部流场为可压缩非定常流场，该流场遵循动量守恒定律、质量守恒定律和能量守恒定律．控制方程为

(16)

(17)

式中：μ为黏度系数，Pa·s；λ为流动介质的传热系数，W/(m·K)；Su和Sv均为源项；ST为黏性耗散项：

(18)

(19)

(20)

湍流模型采用Realizablek-ε方程．模拟介质采用湿空气．采用AUSM格式进行离散，其广泛用于不同音速下气体流场的计算．求解格式选用基于密度求解，时间步长为1×10-6s，计算过程中，保证各时间步下连续性、速度、能量、湍流模型参数、液滴数量、气相质量、液相质量各参数的最终残差值均小于1×10-3．该数值计算需要对计算流域赋予初值，初值不会影响计算结果，为保证计算的收敛性，对全场流域赋予温度298 K、总压0.1 MPa的初值．

3.4 理论模型验证

膨胀激波管的内部流动与压力振荡管相似，因此本文采用埃因霍芬大学进行的水蒸气-氦气混合气在脉冲膨胀激波管内相变的实验数据[22]分别对改进前后的凝结蒸发数值模型进行验证．如图5所示为脉冲膨胀激波管高压区末端处液滴半径变化的实验值与数值模拟结果．

图5 膨胀激波管内液滴半径分布

可见，从液滴半径的大小与变化趋势来看，采用CPA气体状态方程改进过的凝结蒸发模型更能准确描述水蒸气的相变行为．

4 结果与讨论

实验的操作条件：高压进气温度为298 K，总压设置在0.2～0.4 MPa，进气中水蒸气的相对湿度固定不变为100%．低温出口压力设为0.1 MPa，中压回气温度设置为298 K．对于数值模拟，在与实验相同的操作条件下，分别对CPA气体与理想气体在压力振荡管内的制冷过程进行数值模拟并将模拟结果与实验结果进行对比．

在进气压力为0.2 MPa时的温度云图如图6所示．可见在压力振荡管的左侧低温区，CPA气体的温度高于理想气体．提取数值模拟与实验结果并计算制冷效率，可得到制冷效率随进气压力的变化曲线如图7所示．从曲线可知，虽然CPA气体与理想气体的制冷效率变化趋势相同，都是随着进气压力的增加而逐渐上升且上升趋势逐渐平缓，但是CPA气体的压力振荡制冷效率与实验值最为接近，二者最大差值仅为0.5%，最小差值为0.35%．且随着进气压力的增加，CPA气体与理想气体制冷效率差值Δη也逐渐增大．

图6 高压进气0.2 MPa下温度云图

图7 制冷效率与高压进气压力关系

造成CPA气体与理想气体状态方程数值模拟结果存在差距的主要原因是压力振荡管内实际波系与理想波系存在偏差[23]．在压力振荡管中，激波是由多数压缩波叠加形成并且管中压缩波与膨胀波的移动速度为声速．提取计算结果中的声速云图如图8所示．可见在标示出的激波形成区中，CPA气体中的声速与理想气体不同，这导致了激波的形成与运动存在偏差，进一步导致了后续压力波系的偏差．而压力振荡制冷技术依靠压力波的运动实现制冷，并且目前压力振荡管与各端口的匹配基于理想气体内的波系运动进行设计．因此CPA气体中压力波无法与各端口很好地匹配，导致了制冷效率的下降．除此之外，相较于CPA气体，压力波在理想气体中能量传递效率较高，因此经CPA气体状态方程对气体的热力学参数计算进行改进后，压力波在压力振荡管内的能量传递效率有所降低，这也是导致CPA气体制冷效率较理想气体低的原因之一．二者制冷效率的差别随着进气压力的升高而增加，这是因为CPA气体的压力越高，其状态偏离理想气体状态的程度越大．

分别利用改进前后的数值模型对水蒸气在压力振荡管内的液化率进行计算．图9为进气压力0.2 MPa下水蒸气质量分数云图．可见在CPA气体中水蒸气的质量分数较高，说明其凝结成液态水的量较少．提取进出口水蒸气的质量分数计算液化率可得到水蒸气液化率随着进气压力变化的曲线以及CPA气体与理想气体液化率差值Δα的变化曲线，如图10所示．可见CPA气体与理想气体在压力振荡管中液化率的变化趋势相同，即随着进气压力的提升液化率先上升后下降并在进气压力为0.30 MPa时达到最大值．这一变化趋势是由压力振荡管内波系与各端口之间的匹配引起的．对比CPA气体与理想气体中水蒸气的液化率，发现CPA气体中水蒸气的液化率低于理想气体，且两者之间的差值也随着进气压力的增加而逐渐增加．这是因为水蒸气的凝结主要取决于其分压与温度．由于水蒸气在凝结过程中相对湿度保持为100%，其分压仅取决于温度．因此水蒸气的温度决定了其凝结过程．由于CPA气体的制冷效率低于理想气体，压力振荡管内CPA气体的温度高于理想气体，这导致了CPA气体中水蒸气的凝结受到了抑制，液态水的质量减少．随着进气压力的增加，CPA气体与理想气体之间制冷效率差距的增加也导致了液化率差距的增加．

图8 高压进气0.2 MPa下声速云图

图9 高压进气0.2 MPa下水蒸气质量分数云图

图10 水蒸气液化率与高压进气压力的关系

5 结 论

(1)通过与实验数值对比发现，CPA气体状态方程改进后的凝结蒸发模型相较于理想气体模型可以更准确地描述含湿气体凝结蒸发行为，是一种有效的含相变压力振荡数值计算方法，为含湿天然气压力振荡制冷提供了理论基础．

(2)随着高压进气压力从0.2 MPa增加至0.4 MPa，CPA改进模型和理想气体模型得到的制冷效率的差距逐渐增大，其中CPA改进模型数值偏低，更接近实验结果，计算得到的制冷效率变化趋势也更符合实际．

(3)随着高压进气压力从0.2 MPa增加至0.4 MPa，CPA改进模型获得的压力振荡管内含湿空气液化率的变化规律与理想混合气相似．在数值上，CPA改进模型得到的含湿气体液化率低于理想气体，且两者差距随着高压进气压力的升高而增加．该结论可为含湿天然气制冷后的气液分离过程提供设计依据．