楊 玉，吳家榮，韓煜航，吳帥帥，蔣世希，李凱倫，張一帆，李紅智

（西安熱工研究院有限公司，陜西西安 710054）

1 引言

以超臨界二氧化碳（CO2）為工質(zhì)的布雷頓循環(huán)相比于傳統(tǒng)蒸汽朗肯循環(huán)具有更高的效率、更好的靈活性和更緊湊簡(jiǎn)單的機(jī)組結(jié)構(gòu)［1-3］，將在我國(guó)以新能源為主體的新型電力系統(tǒng)中起到“壓艙石”的托底作用［4］。

國(guó)內(nèi)外主流研究機(jī)構(gòu)對(duì)超臨界CO2循環(huán)開展了大量的研究工作。美國(guó)Sandia實(shí)驗(yàn)室搭建了240 kW的超臨界CO2再壓縮實(shí)驗(yàn)平臺(tái)［5-6］，采用了2個(gè)透平-電機(jī)-壓縮機(jī)模塊，開展了深入細(xì)致的實(shí)驗(yàn)測(cè)試工作。Bechtel海軍推進(jìn)公司開發(fā)了一個(gè)雙軸的簡(jiǎn)單回?zé)岢R界CO2發(fā)電系統(tǒng)［7］，采用了一個(gè)可變轉(zhuǎn)速的透平和定轉(zhuǎn)速的壓縮機(jī)，設(shè)計(jì)輸出功率為100 kW。日本的應(yīng)用能源研究所（IAE）與東京工業(yè)大學(xué)（TIT）合作建設(shè)了一臺(tái)10 kW的超臨界CO2實(shí)驗(yàn)裝置，用于小尺寸的透平機(jī)械和發(fā)電的研究［8］。韓國(guó)原子能研究院（KEARI）、韓國(guó)先進(jìn)科技大學(xué)（KAIST）和浦項(xiàng)科技大學(xué)（POSTECH）聯(lián)合設(shè)計(jì)了一臺(tái)80 kW 的超臨界CO2循環(huán)試驗(yàn)裝，并分階段建設(shè)［9］。韓國(guó)能源研究院（KIER）也正在建設(shè)用于余熱回收領(lǐng)域的超臨界CO2試驗(yàn)裝置，加熱器功率約為648 kW［10］。我國(guó)西安熱工研究院有限公司于2021 年成功投運(yùn)了5 MWe 超臨界CO2布雷頓循環(huán)發(fā)電機(jī)組，實(shí)現(xiàn)23.72%的循環(huán)效率［11］。中國(guó)科學(xué)院工程熱物理研究所搭建了超臨界CO2壓縮機(jī)測(cè)試系統(tǒng)，并完成了對(duì)1 MW 超臨界CO2壓縮機(jī)的性能測(cè)試工作［12］。綜上所述，目前國(guó)際上美國(guó)、日本、韓國(guó)、法國(guó)均在系統(tǒng)研究、機(jī)組設(shè)計(jì)、零部件加工及系統(tǒng)示范驗(yàn)證方面開展了較多的工作，美國(guó)預(yù)期在近年內(nèi)完成商業(yè)示范，在3～5年將形成1～10 MW級(jí)超臨界二氧化碳機(jī)組完整產(chǎn)業(yè)鏈和成功商業(yè)案例。我國(guó)西安熱工研究院下一步將開展50 MW超臨界CO2機(jī)組的商業(yè)示范工程的建設(shè)。

壓縮機(jī)是超臨界CO2布雷頓循環(huán)的關(guān)鍵核心設(shè)備之一，壓縮機(jī)入口參數(shù)越接近臨界點(diǎn)，則功耗越小，系統(tǒng)效率越高。但是當(dāng)壓縮機(jī)入口接近臨界點(diǎn)時(shí)，壓縮機(jī)中的局部加速現(xiàn)象會(huì)引發(fā)兩相區(qū)的出現(xiàn)［13］。壓縮機(jī)內(nèi)的局部凝結(jié)會(huì)對(duì)壓縮機(jī)的高效安全運(yùn)行產(chǎn)生不可忽視的影響。美國(guó)Sandia實(shí)驗(yàn)室的試驗(yàn)測(cè)試結(jié)果表明，當(dāng)壓縮機(jī)入口溫度和壓力進(jìn)入兩相區(qū)時(shí)，質(zhì)量流量和密度開始脈動(dòng)，脈動(dòng)周期與工質(zhì)循環(huán)呈強(qiáng)耦合關(guān)系［14］。朱玉銘等［12］的試驗(yàn)測(cè)試結(jié)果表明，兩相流動(dòng)在壓縮機(jī)內(nèi)的局部分層以及流動(dòng)摻混會(huì)導(dǎo)致壓縮機(jī)效率的降低。針對(duì)這一問題，國(guó)內(nèi)外學(xué)者開展了大量的數(shù)值模擬研究，以期揭示壓縮機(jī)內(nèi)的CO2凝結(jié)規(guī)律及其對(duì)壓縮機(jī)性能的影響機(jī)制。Bao 等［13］采用NUMEDA Fine/turbo 軟件對(duì)壓縮機(jī)內(nèi)的CO2流動(dòng)情況進(jìn)行了模擬研究，基于均相平衡模型評(píng)價(jià)了壓縮機(jī)內(nèi)的兩相流動(dòng)和冷凝情況。Romei等［15］基于ANSYS Fluent 軟件采用RANS 方法對(duì)一臺(tái)50 MW 離心壓縮機(jī)的流動(dòng)進(jìn)行進(jìn)行了模擬研究，獲得了壓縮機(jī)進(jìn)口葉片附近的兩相區(qū)分布。Xu等［16］采用ANSYS CFX軟件對(duì)超臨界CO2壓縮機(jī)內(nèi)的流動(dòng)特性進(jìn)行了研究，獲得了兩相區(qū)域的分布。上述數(shù)值模擬研究都把兩相區(qū)的流體當(dāng)作單相來處理，通過劃定凝結(jié)溫度或壓力來獲得相變區(qū)域的分布，這對(duì)研究超臨界CO2壓縮機(jī)的相變規(guī)律和提出改進(jìn)措施起到了很大作用。Lettieri 等的試驗(yàn)觀測(cè)表明，超臨界CO2進(jìn)入兩相區(qū)的凝結(jié)是呈現(xiàn)霧狀液滴分布的［17］，上述方法都不能描述CO2的凝結(jié)成霧特點(diǎn)，也無(wú)法追蹤C(jī)O2液滴演變過程，難以滿足進(jìn)一步提高超臨界CO2壓縮機(jī)性能的研究需求。

鑒于此，本文結(jié)合經(jīng)典成核和生長(zhǎng)模型將離散顆粒模型（discrete particle model，DPM）用于描述CO2霧化液滴的運(yùn)行規(guī)律，開發(fā)出超臨界CO2凝結(jié)的CFD-DPM 數(shù)值模擬方法，研究了Laval 噴管中的CO2凝結(jié)特性并進(jìn)行驗(yàn)證，提出的模擬方法為超臨界CO2的壓縮機(jī)設(shè)計(jì)和優(yōu)化奠定基礎(chǔ)。

2 研究方法

2.1 物理模型

本文研究的Laval 噴管模型如圖1 所示，噴管進(jìn)口尺寸為12.7 mm×10.0 mm，總長(zhǎng)為98.4 mm，其具體的幾何參數(shù)取自Lettieri 等［17］的試驗(yàn)工作。該Laval噴管進(jìn)口CO2參數(shù)為8 MPa/311 K，從超臨界狀態(tài)膨脹到兩相區(qū)。進(jìn)入到兩相區(qū)的CO2并不會(huì)立即凝結(jié)為液態(tài)，而是當(dāng)其參數(shù)超過Wilson 線后才會(huì)發(fā)生相變。過飽和狀態(tài)的CO2工質(zhì)物性采用基于氣相區(qū)域的物性向兩相區(qū)外推獲得，圖2 所示為工質(zhì)焓值在兩相區(qū)的外推計(jì)算過程。圖中的離散點(diǎn)為基于NIST 物性數(shù)據(jù)庫(kù)的焓值，采用2 元5 次多項(xiàng)式擬合，獲得外推物性曲面。本文的模擬中，進(jìn)入兩相區(qū)的工質(zhì)物性采用外推的曲面上參數(shù)表示。

圖1 Laval噴管的物理模型Fig.1 Physical model of the Laval Nozzle

圖2 兩相區(qū)物性外推Fig.2 Temperature extrapolation for the two-phase zone

2.2 數(shù)學(xué)模型

2.2.1 液滴成核和生長(zhǎng)模型

成核速率J定義為［18-19］：

式中：σ為液滴表面張力，Mm為單個(gè)水分子質(zhì)量，ρl為液滴密度，ρv為蒸汽密度，k為玻爾茲曼常數(shù)，T為當(dāng)?shù)卣羝麥囟?，ΔG為吉布斯自由能變化，其計(jì)算公式為：

式中S為過飽和度，其表達(dá)式為：

式中：P為蒸汽壓力，Psat（T）為平衡相變的平衡壓力。

Kelvin-Helmholtz液滴的臨界半徑計(jì)算公示為：

式中Rv為二氧化碳的氣體常數(shù)。

液滴一經(jīng)形成，其核心就會(huì)以一定的速率長(zhǎng)大，該生長(zhǎng)速率取決于蒸汽與液滴之間的熱質(zhì)平衡。Young［20-21］提出的液滴生長(zhǎng)速率模型為：

式中：hvl，γ，Cp分別為潛熱、比熱容比和比熱，Td為液滴溫度，其計(jì)算公式為：

式中r為液滴的半徑。

2.2.2 連續(xù)相和離散相控制方程

連續(xù)相的質(zhì)量、動(dòng)量和能量方程如式（7）～（9）所示：

式中：ρ，u，p和H分別為密度、速度、壓力和總焓，τ為應(yīng)力張量，λeff為有效導(dǎo)熱系數(shù)，F(xiàn)int為蒸汽與液滴之間的相互作用力。Sm，Su和Sh分別為由于凝結(jié)產(chǎn)生的質(zhì)量、動(dòng)量和能量源項(xiàng)。

液滴的運(yùn)動(dòng)采用基于牛頓第二定律的拉格朗日方法計(jì)算：

式中：mp為液滴質(zhì)量，up為液滴速度，F(xiàn)D為連續(xù)相與離散相之間的曳力。

冷凝過程產(chǎn)生的液滴質(zhì)量增長(zhǎng)速率為：

式中N為單位體積的液滴數(shù)量。式（11）右側(cè)的第一項(xiàng)表示由于蒸汽凝結(jié)產(chǎn)生的液滴質(zhì)量增長(zhǎng)，第二項(xiàng)表示成核后液滴的長(zhǎng)大速率。連續(xù)相的質(zhì)量、動(dòng)量和能量源項(xiàng)與液滴質(zhì)量變化之間的關(guān)系為：

本文采用SSTk-ω湍流模型封閉連續(xù)相的動(dòng)量方程，采用SIMPLE 方法計(jì)算速度和壓力耦合，采用兩步耦合方法計(jì)算連續(xù)相與離散的液滴相之間的相互影響。兩相區(qū)以外的二氧化碳物性計(jì)算基于NIST數(shù)據(jù)庫(kù)，過飽和二氧化碳的物性通過將氣相物性外延至兩相區(qū)內(nèi)獲取。

3 結(jié)果與討論

為了驗(yàn)證提出的數(shù)值模擬方法，將計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，如圖3和圖4所示。圖3所示為噴管中心的相對(duì)壓力沿著流動(dòng)方向的變化規(guī)律，從圖中可以看出，進(jìn)入噴管后的CO2工質(zhì)壓力呈單調(diào)下降趨勢(shì)。在x/L=-0.3位置的下游，壓力急劇下降，在x/L=0 的喉口位置前后下降趨勢(shì)變緩，然后在x/L=0.2下游以后壓力大致呈線性下降趨勢(shì)。從x/L=-0.3到喉口位置的壓力急劇下降對(duì)應(yīng)于Laval 噴管的收縮段，但是喉口位置出現(xiàn)了工質(zhì)凝結(jié)現(xiàn)象，凝結(jié)過程有熱量釋放，導(dǎo)致工質(zhì)的壓力下降速度變緩。同時(shí)由于部分工質(zhì)凝結(jié)釋放的熱量使工質(zhì)參數(shù)向飽和氣相線偏移，凝結(jié)工質(zhì)比率降低，因而隨著工質(zhì)的繼續(xù)加速下游的壓力降低曲線的斜率大致相同。

圖3 實(shí)驗(yàn)和數(shù)值模擬的壓力分布對(duì)比Fig.3 Comparison of the pressure distribution between the experimental and numerical results

圖4 實(shí)驗(yàn)與數(shù)值模擬的液滴分布對(duì)比Fig.4 Comparison of the droplet distribution between the experimental and numerical results

圖4所示為數(shù)值模擬的液滴分布與可視化實(shí)驗(yàn)拍攝的照片的對(duì)比。圖4（a）所示為實(shí)驗(yàn)拍攝的CO2在Laval噴管中的凝結(jié)情況，圖中白色的霧狀區(qū)域?qū)?yīng)于凝結(jié)的CO2液滴在噴管中的分布，這表明凝結(jié)的CO2液滴是彌散在流體中的，且跟隨流體一起流動(dòng)。圖4（b）所示為本文提出的CFD-DPM 方法模擬出來的CO2液滴分布，可以看出本方法很好滴再現(xiàn)了凝結(jié)的CO2液滴在噴管中的霧狀分布情況。與雙歐拉模型相比，該新方法可以更好的展現(xiàn)CO2凝結(jié)液滴的運(yùn)動(dòng)變化規(guī)律。此外，從圖4中還可以看出，CO2凝結(jié)的位置發(fā)生在喉口前，數(shù)值模擬結(jié)果與實(shí)驗(yàn)測(cè)試結(jié)果吻合。圖3 和圖4 中的數(shù)值模擬與實(shí)驗(yàn)測(cè)試的結(jié)果對(duì)比表明，本文提出的CFD-DPM 方法不僅是可靠的，還具有詳細(xì)描述凝結(jié)液滴演化規(guī)律的功能，如果應(yīng)用在超臨界CO2的模擬中，能更好的展示壓縮機(jī)葉輪進(jìn)口的凝結(jié)液滴的發(fā)展演化規(guī)律，為高性能超臨界CO2壓縮機(jī)的設(shè)計(jì)和優(yōu)化提供高效手段。

圖5 至圖11 所示為L(zhǎng)aval 噴管中的關(guān)鍵參數(shù)分布云圖。圖5 所示為壓力分布云圖，進(jìn)口壓力為8 MPa，出口壓力約為2.85 MPa，喉部壓力約為5.0 MPa。圖6 所示為溫度分布云圖，進(jìn)口溫度為310 K，出口溫度為251 K，喉部溫度為277 K。根據(jù)文獻(xiàn)［17］，8 MPa/310 K絕熱膨脹線與氣相線相交的參數(shù)約為6.87 MPa/301 K，與Wilson 曲線相交的參數(shù)約為299.6 K。根據(jù)圖5 和圖6，x=-0.013 位置為凝結(jié)開始點(diǎn)，即在該處工質(zhì)CO2狀態(tài)參數(shù)與Wilson 曲線相交，而在CO2達(dá)到凝結(jié)之前的狀態(tài)為亞穩(wěn)定態(tài)。圖7 所示為噴管中的速度分布云圖，進(jìn)口的速度約為26 m/s，喉部的速度為143 m/s，出口速度約為182 m/s。圖8 所示為過飽和度的分布云圖，噴管進(jìn)口的過飽和度為0.96 左右，在x=-0.015 位置處，過飽和度達(dá)到1，在凝結(jié)開始點(diǎn)過飽和度為1.045，在噴管喉部為1.285，在噴管出口過飽和度達(dá)到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 所示為臨界半徑的分布云圖，臨界半徑為某一狀態(tài)下液滴凝結(jié)時(shí)的最小液滴半徑。從圖中可以看出，凝結(jié)時(shí)液滴臨界半徑約在0.65 nm 到0.82 nm 范圍內(nèi)。開始凝結(jié)時(shí)，臨界半徑約為0.69 nm。在x=-0.05 m 到喉部，臨界半徑降低為0.65 nm左右。隨后沿著流動(dòng)方向臨界半徑逐漸增大，在出口達(dá)到0.82 nm 左右。在喉部前的臨界半徑減小的這一過程主要是由于CO2凝結(jié)放熱導(dǎo)致工質(zhì)狀態(tài)的過飽和度降低導(dǎo)致的。圖10 所示為液滴凝結(jié)速率的分布云圖，即單位體積的工質(zhì)在單位時(shí)間內(nèi)生成的液滴數(shù)量。從圖中可以看出成核率最大的位置出現(xiàn)在x=-0.012 m 位置處，最大值約為2.58×1025個(gè)/（m3·s），該位置之前和之后的成核率都降低。圖11 所示為凝結(jié)的液滴質(zhì)量生長(zhǎng)速率分布云圖，從圖中可以看出液滴的質(zhì)量生長(zhǎng)速率最大位置與成核率最大的位置重合，都在x=-0.012 m 位置前后，最大值約為26.8 kg/（m3·s）。從x=-0.014 m 到喉部是液滴質(zhì)量集中生成的區(qū)域，其他區(qū)域液滴增加的質(zhì)量較小，這說明了在Laval噴管的漸縮部分工質(zhì)壓力急劇降低，達(dá)到了很高的過飽和度狀態(tài)，穿過亞穩(wěn)態(tài)區(qū)進(jìn)入凝結(jié)區(qū)域，創(chuàng)造極為有利于凝結(jié)的環(huán)境，因此CO2工質(zhì)在這個(gè)位置集中凝結(jié)。在CO2凝結(jié)的同時(shí)，釋放出來的熱量使工質(zhì)狀態(tài)向亞穩(wěn)態(tài)區(qū)域偏移，液滴凝結(jié)的驅(qū)動(dòng)力減弱，液滴的成核速率和生長(zhǎng)速率減緩。

圖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中的質(zhì)量生長(zhǎng)速率分布Fig.11 Distribution of mass growing rate in the Laval nozzle

圖12和圖13所示為L(zhǎng)aval噴管中的凝結(jié)液滴分布情況，在噴管漸縮區(qū)域生成的液滴跟隨流體向下游流動(dòng)并長(zhǎng)大，液滴在噴管中的停留時(shí)間約為0.000 3 s，單個(gè)液滴的平均質(zhì)量從剛生成時(shí)的9×10-25kg增長(zhǎng)到出口的（2.0～3.0）×10-24kg。

圖12 Laval中的凝結(jié)液滴分布Fig.12 Droplets distribution in the Laval nozzle

圖13 Laval中的液滴質(zhì)量分布Fig.13 Mass distribution of droplets in the Laval nozzle

4 結(jié)論

本文采用歐拉-拉格朗日數(shù)值模擬方法研究了Laval噴管中的CO2凝結(jié)特性，主要結(jié)論如下：

1）將液滴凝結(jié)經(jīng)典成核和生長(zhǎng)模型以離散顆粒模型（DPM）的形式嵌入流體計(jì)算模型，開發(fā)出描述超臨界CO2凝結(jié)CFD-DPM 數(shù)值模擬方法，其可靠性得到了實(shí)驗(yàn)驗(yàn)證。

2）本研究中Laval 噴管的漸縮部分工質(zhì)壓力急劇降低，達(dá)到了很高的過飽和度狀態(tài)，穿過亞穩(wěn)態(tài)區(qū)進(jìn)入凝結(jié)區(qū)域，創(chuàng)造極為有利于凝結(jié)的環(huán)境，因此CO2工質(zhì)在這個(gè)位置集中凝結(jié)。在CO2凝結(jié)的同時(shí)，釋放出來的熱量使工質(zhì)狀態(tài)向亞穩(wěn)態(tài)區(qū)域偏移，液滴凝結(jié)的驅(qū)動(dòng)力減弱，液滴的成核速率和生長(zhǎng)速率減緩。

3）提出的CFD-DPM 方法具有詳細(xì)描述凝結(jié)液滴演化規(guī)律的功能，可應(yīng)用在超臨界CO2壓縮機(jī)的模擬中，有望展示壓縮機(jī)葉輪進(jìn)口的凝結(jié)液滴的發(fā)展演化規(guī)律，為高性能超臨界CO2壓縮機(jī)的設(shè)計(jì)和優(yōu)化提供高效手段。