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壓縮機的設計和優化提供高效手段。