Laval喷管中CO2凝结特性数值模拟研究

杨 玉，吴家荣，韩煜航，吴帅帅，蒋世希，李凯伦，张一帆，李红智

（西安热工研究院有限公司，陕西西安 710054）

1 引言

以超临界二氧化碳（CO2）为工质的布雷顿循环相比于传统蒸汽朗肯循环具有更高的效率、更好的灵活性和更紧凑简单的机组结构［1-3］，将在我国以新能源为主体的新型电力系统中起到“压舱石”的托底作用［4］。

国内外主流研究机构对超临界CO2循环开展了大量的研究工作。美国Sandia实验室搭建了240 kW的超临界CO2再压缩实验平台［5-6］，采用了2个透平-电机-压缩机模块，开展了深入细致的实验测试工作。Bechtel海军推进公司开发了一个双轴的简单回热超临界CO2发电系统［7］，采用了一个可变转速的透平和定转速的压缩机，设计输出功率为100 kW。日本的应用能源研究所（IAE）与东京工业大学（TIT）合作建设了一台10 kW的超临界CO2实验装置，用于小尺寸的透平机械和发电的研究［8］。韩国原子能研究院（KEARI）、韩国先进科技大学（KAIST）和浦项科技大学（POSTECH）联合设计了一台80 kW 的超临界CO2循环试验装，并分阶段建设［9］。韩国能源研究院（KIER）也正在建设用于余热回收领域的超临界CO2试验装置，加热器功率约为648 kW［10］。我国西安热工研究院有限公司于2021 年成功投运了5 MWe 超临界CO2布雷顿循环发电机组，实现23.72%的循环效率［11］。中国科学院工程热物理研究所搭建了超临界CO2压缩机测试系统，并完成了对1 MW 超临界CO2压缩机的性能测试工作［12］。综上所述，目前国际上美国、日本、韩国、法国均在系统研究、机组设计、零部件加工及系统示范验证方面开展了较多的工作，美国预期在近年内完成商业示范，在3～5年将形成1～10 MW级超临界二氧化碳机组完整产业链和成功商业案例。我国西安热工研究院下一步将开展50 MW超临界CO2机组的商业示范工程的建设。

压缩机是超临界CO2布雷顿循环的关键核心设备之一，压缩机入口参数越接近临界点，则功耗越小，系统效率越高。但是当压缩机入口接近临界点时，压缩机中的局部加速现象会引发两相区的出现［13］。压缩机内的局部凝结会对压缩机的高效安全运行产生不可忽视的影响。美国Sandia实验室的试验测试结果表明，当压缩机入口温度和压力进入两相区时，质量流量和密度开始脉动，脉动周期与工质循环呈强耦合关系［14］。朱玉铭等［12］的试验测试结果表明，两相流动在压缩机内的局部分层以及流动掺混会导致压缩机效率的降低。针对这一问题，国内外学者开展了大量的数值模拟研究，以期揭示压缩机内的CO2凝结规律及其对压缩机性能的影响机制。Bao 等［13］采用NUMEDA Fine/turbo 软件对压缩机内的CO2流动情况进行了模拟研究，基于均相平衡模型评价了压缩机内的两相流动和冷凝情况。Romei等［15］基于ANSYS Fluent 软件采用RANS 方法对一台50 MW 离心压缩机的流动进行进行了模拟研究，获得了压缩机进口叶片附近的两相区分布。Xu等［16］采用ANSYS CFX软件对超临界CO2压缩机内的流动特性进行了研究，获得了两相区域的分布。上述数值模拟研究都把两相区的流体当作单相来处理，通过划定凝结温度或压力来获得相变区域的分布，这对研究超临界CO2压缩机的相变规律和提出改进措施起到了很大作用。Lettieri 等的试验观测表明，超临界CO2进入两相区的凝结是呈现雾状液滴分布的［17］，上述方法都不能描述CO2的凝结成雾特点，也无法追踪CO2液滴演变过程，难以满足进一步提高超临界CO2压缩机性能的研究需求。

鉴于此，本文结合经典成核和生长模型将离散颗粒模型（discrete particle model，DPM）用于描述CO2雾化液滴的运行规律，开发出超临界CO2凝结的CFD-DPM 数值模拟方法，研究了Laval 喷管中的CO2凝结特性并进行验证，提出的模拟方法为超临界CO2的压缩机设计和优化奠定基础。

2 研究方法

2.1 物理模型

本文研究的Laval 喷管模型如图1 所示，喷管进口尺寸为12.7 mm×10.0 mm，总长为98.4 mm，其具体的几何参数取自Lettieri 等［17］的试验工作。该Laval喷管进口CO2参数为8 MPa/311 K，从超临界状态膨胀到两相区。进入到两相区的CO2并不会立即凝结为液态，而是当其参数超过Wilson 线后才会发生相变。过饱和状态的CO2工质物性采用基于气相区域的物性向两相区外推获得，图2 所示为工质焓值在两相区的外推计算过程。图中的离散点为基于NIST 物性数据库的焓值，采用2 元5 次多项式拟合，获得外推物性曲面。本文的模拟中，进入两相区的工质物性采用外推的曲面上参数表示。

图1 Laval喷管的物理模型Fig.1 Physical model of the Laval Nozzle

图2 两相区物性外推Fig.2 Temperature extrapolation for the two-phase zone

2.2 数学模型

2.2.1 液滴成核和生长模型

成核速率J定义为［18-19］：

式中：σ为液滴表面张力，Mm为单个水分子质量，ρl为液滴密度，ρv为蒸汽密度，k为玻尔兹曼常数，T为当地蒸汽温度，ΔG为吉布斯自由能变化，其计算公式为：

式中S为过饱和度，其表达式为：

式中：P为蒸汽压力，Psat（T）为平衡相变的平衡压力。

Kelvin-Helmholtz液滴的临界半径计算公示为：

式中Rv为二氧化碳的气体常数。

液滴一经形成，其核心就会以一定的速率长大，该生长速率取决于蒸汽与液滴之间的热质平衡。Young［20-21］提出的液滴生长速率模型为：

式中：hvl，γ，Cp分别为潜热、比热容比和比热，Td为液滴温度，其计算公式为：

式中r为液滴的半径。

2.2.2 连续相和离散相控制方程

连续相的质量、动量和能量方程如式（7）～（9）所示：

式中：ρ，u，p和H分别为密度、速度、压力和总焓，τ为应力张量，λeff为有效导热系数，Fint为蒸汽与液滴之间的相互作用力。Sm，Su和Sh分别为由于凝结产生的质量、动量和能量源项。

液滴的运动采用基于牛顿第二定律的拉格朗日方法计算：

式中：mp为液滴质量，up为液滴速度，FD为连续相与离散相之间的曳力。

冷凝过程产生的液滴质量增长速率为：

式中N为单位体积的液滴数量。式（11）右侧的第一项表示由于蒸汽凝结产生的液滴质量增长，第二项表示成核后液滴的长大速率。连续相的质量、动量和能量源项与液滴质量变化之间的关系为：

本文采用SSTk-ω湍流模型封闭连续相的动量方程，采用SIMPLE 方法计算速度和压力耦合，采用两步耦合方法计算连续相与离散的液滴相之间的相互影响。两相区以外的二氧化碳物性计算基于NIST数据库，过饱和二氧化碳的物性通过将气相物性外延至两相区内获取。

3 结果与讨论

为了验证提出的数值模拟方法，将计算结果与实验数据对比，如图3和图4所示。图3所示为喷管中心的相对压力沿着流动方向的变化规律，从图中可以看出，进入喷管后的CO2工质压力呈单调下降趋势。在x/L=-0.3位置的下游，压力急剧下降，在x/L=0 的喉口位置前后下降趋势变缓，然后在x/L=0.2下游以后压力大致呈线性下降趋势。从x/L=-0.3到喉口位置的压力急剧下降对应于Laval 喷管的收缩段，但是喉口位置出现了工质凝结现象，凝结过程有热量释放，导致工质的压力下降速度变缓。同时由于部分工质凝结释放的热量使工质参数向饱和气相线偏移，凝结工质比率降低，因而随着工质的继续加速下游的压力降低曲线的斜率大致相同。

图3 实验和数值模拟的压力分布对比Fig.3 Comparison of the pressure distribution between the experimental and numerical results

图4 实验与数值模拟的液滴分布对比Fig.4 Comparison of the droplet distribution between the experimental and numerical results

图4所示为数值模拟的液滴分布与可视化实验拍摄的照片的对比。图4（a）所示为实验拍摄的CO2在Laval喷管中的凝结情况，图中白色的雾状区域对应于凝结的CO2液滴在喷管中的分布，这表明凝结的CO2液滴是弥散在流体中的，且跟随流体一起流动。图4（b）所示为本文提出的CFD-DPM 方法模拟出来的CO2液滴分布，可以看出本方法很好滴再现了凝结的CO2液滴在喷管中的雾状分布情况。与双欧拉模型相比，该新方法可以更好的展现CO2凝结液滴的运动变化规律。此外，从图4中还可以看出，CO2凝结的位置发生在喉口前，数值模拟结果与实验测试结果吻合。图3 和图4 中的数值模拟与实验测试的结果对比表明，本文提出的CFD-DPM 方法不仅是可靠的，还具有详细描述凝结液滴演化规律的功能，如果应用在超临界CO2的模拟中，能更好的展示压缩机叶轮进口的凝结液滴的发展演化规律，为高性能超临界CO2压缩机的设计和优化提供高效手段。

图5 至图11 所示为Laval 喷管中的关键参数分布云图。图5 所示为压力分布云图，进口压力为8 MPa，出口压力约为2.85 MPa，喉部压力约为5.0 MPa。图6 所示为温度分布云图，进口温度为310 K，出口温度为251 K，喉部温度为277 K。根据文献［17］，8 MPa/310 K绝热膨胀线与气相线相交的参数约为6.87 MPa/301 K，与Wilson 曲线相交的参数约为299.6 K。根据图5 和图6，x=-0.013 位置为凝结开始点，即在该处工质CO2状态参数与Wilson 曲线相交，而在CO2达到凝结之前的状态为亚稳定态。图7 所示为喷管中的速度分布云图，进口的速度约为26 m/s，喉部的速度为143 m/s，出口速度约为182 m/s。图8 所示为过饱和度的分布云图，喷管进口的过饱和度为0.96 左右，在x=-0.015 位置处，过饱和度达到1，在凝结开始点过饱和度为1.045，在喷管喉部为1.285，在喷管出口过饱和度达到1.5。

图5 Laval中的压力分布Fig.5 Pressure distribution in the Laval nozzle

图6 Laval中的温度分布Fig.6 Temperature distribution in the Laval nozzle

图7 Laval中的速度分布Fig.7 Velocity distribution in the Laval nozzle

图8 Laval中的过饱和度分布Fig.8 Distribution of the supersaturation degree in the Laval nozzle

图9 所示为临界半径的分布云图，临界半径为某一状态下液滴凝结时的最小液滴半径。从图中可以看出，凝结时液滴临界半径约在0.65 nm 到0.82 nm 范围内。开始凝结时，临界半径约为0.69 nm。在x=-0.05 m 到喉部，临界半径降低为0.65 nm左右。随后沿着流动方向临界半径逐渐增大，在出口达到0.82 nm 左右。在喉部前的临界半径减小的这一过程主要是由于CO2凝结放热导致工质状态的过饱和度降低导致的。图10 所示为液滴凝结速率的分布云图，即单位体积的工质在单位时间内生成的液滴数量。从图中可以看出成核率最大的位置出现在x=-0.012 m 位置处，最大值约为2.58×1025个/（m3·s），该位置之前和之后的成核率都降低。图11 所示为凝结的液滴质量生长速率分布云图，从图中可以看出液滴的质量生长速率最大位置与成核率最大的位置重合，都在x=-0.012 m 位置前后，最大值约为26.8 kg/（m3·s）。从x=-0.014 m 到喉部是液滴质量集中生成的区域，其他区域液滴增加的质量较小，这说明了在Laval喷管的渐缩部分工质压力急剧降低，达到了很高的过饱和度状态，穿过亚稳态区进入凝结区域，创造极为有利于凝结的环境，因此CO2工质在这个位置集中凝结。在CO2凝结的同时，释放出来的热量使工质状态向亚稳态区域偏移，液滴凝结的驱动力减弱，液滴的成核速率和生长速率减缓。

图9 Laval中的临界半径分布Fig.9 Distribution of critical radius in the Laval nozzle

图10 Laval中的成核率分布Fig.10 Distribution of nucleating rate in the Laval nozzle

图11 Laval中的质量生长速率分布Fig.11 Distribution of mass growing rate in the Laval nozzle

图12和图13所示为Laval喷管中的凝结液滴分布情况，在喷管渐缩区域生成的液滴跟随流体向下游流动并长大，液滴在喷管中的停留时间约为0.000 3 s，单个液滴的平均质量从刚生成时的9×10-25kg增长到出口的（2.0～3.0）×10-24kg。

图12 Laval中的凝结液滴分布Fig.12 Droplets distribution in the Laval nozzle

图13 Laval中的液滴质量分布Fig.13 Mass distribution of droplets in the Laval nozzle

4 结论

本文采用欧拉-拉格朗日数值模拟方法研究了Laval喷管中的CO2凝结特性，主要结论如下：

1）将液滴凝结经典成核和生长模型以离散颗粒模型（DPM）的形式嵌入流体计算模型，开发出描述超临界CO2凝结CFD-DPM 数值模拟方法，其可靠性得到了实验验证。

2）本研究中Laval 喷管的渐缩部分工质压力急剧降低，达到了很高的过饱和度状态，穿过亚稳态区进入凝结区域，创造极为有利于凝结的环境，因此CO2工质在这个位置集中凝结。在CO2凝结的同时，释放出来的热量使工质状态向亚稳态区域偏移，液滴凝结的驱动力减弱，液滴的成核速率和生长速率减缓。

3）提出的CFD-DPM 方法具有详细描述凝结液滴演化规律的功能，可应用在超临界CO2压缩机的模拟中，有望展示压缩机叶轮进口的凝结液滴的发展演化规律，为高性能超临界CO2压缩机的设计和优化提供高效手段。