蒸汽浸沒射流與凝結(jié)數(shù)值仿真分析

汪廣懷，孫露，王暢，郝銳，王建軍

1中國艦船研究設(shè)計中心，湖北武漢430064

2哈爾濱工程大學核科學與技術(shù)學院，黑龍江哈爾濱150001

0 引 言

蒸汽浸沒射流是汽水直接接觸冷凝的一種形式，因其良好的傳熱傳質(zhì)特性，在制冷、石油化工、航天航空及核能等工業(yè)領(lǐng)域的應(yīng)用較為廣泛［1］。船用蒸汽動力系統(tǒng)進行蒸汽排放時，也是通過射流冷凝的方式實現(xiàn)高效換熱。高速蒸汽射入冷水而形成的穩(wěn)定連續(xù)蒸汽區(qū)即汽羽，在穩(wěn)定射流中形成的汽羽能夠反映浸沒射流的傳熱特性，所以國內(nèi)外開展了大量的研究工作。曹巖等［2］針對0.2～0.5 MPa飽和蒸汽在20～70℃過冷水中超音速浸沒射流凝結(jié)而形成的汽羽形狀開展了實驗研究，其根據(jù)汽羽膨脹的次數(shù)提出了5類汽羽形狀，得到了汽羽的穿透深度隨著蒸汽壓力增大和過冷水溫度上升而逐漸增加的結(jié)論。Zhou等［3］運用CFD方法對超音速浸沒射流冷凝進行數(shù)值分析，得到了汽羽在凝水中的貫穿長度和膨脹比將隨著環(huán)境壓力的增加而降低的結(jié)論，同時還發(fā)現(xiàn)噴管出口處汽相的體積波動將隨著環(huán)境壓力的下降而上升。袁方等［4］針對音速與超音速蒸汽浸沒射流核心汽羽區(qū)的穿透長度開展實驗研究，得到了汽羽的無量綱穿透長度將隨著蒸汽質(zhì)量流率的增大和過冷水溫度的上升而逐漸增加的結(jié)論，同時發(fā)現(xiàn)超音速蒸汽射流的無量綱穿透長度小于音速蒸汽射流。王超等［5］研究了超音速氣體射流在水中的噴射過程，并基于流體體積函數(shù)（Volume of Fluid，VOF）方法對浸沒射流過程開展二維仿真分析，發(fā)現(xiàn)了射流出口處氣泡生長的周期性及伴隨的脹鼓現(xiàn)象。Isaev等［6］采用修正的SST湍流模型，基于 Zapryagaev等［7］的實驗結(jié)果，精確模擬了浸沒射流汽羽，發(fā)現(xiàn)與實驗結(jié)果高度吻合。

綜上所述，壓力、凝水溫度、噴管口徑、浸沒深度等因素對浸沒射流冷凝的傳熱傳質(zhì)特性影響較大。目前的研究成果主要與汽羽形態(tài)、長度等空間分布有關(guān)，然而，對其形成過程中汽羽內(nèi)部速度場、壓力場等流體細節(jié)卻沒有形成清晰的認識。本文擬基于前期的蒸汽浸沒射流實驗成果，利用CFD仿真軟件ANSYS Fluent進行仿真模擬，分析蒸汽進口壓力、凝水溫度等因素對浸沒射流汽羽的影響，同時研究汽羽內(nèi)部的速度場和壓力場分布情況，用以為浸沒射流的高效換熱應(yīng)用場景設(shè)計提供指導。

1 CFD模型

ANSYS Fluent中的Mixture兩相流模型是基于Euler-Euler兩相流方程開發(fā)的應(yīng)用，其遵守的基本控制方程如下。

1）連續(xù)性方程：

式中：ρm為混合相密度；t為時間；vm為混合相基于質(zhì)量平均的速度。

其中，

式中：k=1，2，…，n，為混合相中的不同組分；αk為k相的體積分數(shù)；ρk為k相密度；vk為k相速度。

2）動量方程：

式中：P為壓力；μm為混合相的粘度；g為重力加速度；F為體積力；vm，k為k相相對于混合相的漂移速度。

其中，

式中：μk為k相的粘度。

3）漂移速度：

使用漂移流模型時，1.2節(jié)動量方程的漂移速度可按照如下模型定義：

式中：vpq為p，q兩相間的滑移速度；ck=αkρk/ρm，為k相對應(yīng)的質(zhì)量份額；vqk為主相q相對于混合相中k相的速度。

在Fluent中，相間滑移速度的Manninen模型定義如下：

其中，

式中：dp為次相氣泡的直徑；μq為主相的粘度。

基于Schiller-Naumann模型，fdrag的表達式為

式中，Re為雷諾數(shù)。

4）能量方程：

需要說明的是，本文的相間冷凝換熱系數(shù)需通過實驗結(jié)果擬合的經(jīng)驗關(guān)系式進行自定義。

5）兩熱阻模型：

在相間傳熱過程中，鑒于相界面兩側(cè)分別對應(yīng)不同的換熱系數(shù)，本文將相間換熱模型轉(zhuǎn)換為兩熱阻模型。定義主相（q）與次相（p）之間相界面的溫度為Ts，則相界面?zhèn)鬟f到q相和p相的熱量Qq和Qp分別為：

式中：hq和hp分別為q相和p相的換熱系數(shù)；Ai為相界面面積；Tq和Tp為q相和p相的溫度；m˙pq為由p相轉(zhuǎn)變?yōu)閝相的質(zhì)量速率；Hqs和 Hps分別為q相和p相的的焓值。

由于相間傳遞的熱量之和為0，則相間溫度為

相間傳熱的兩熱阻模型如下：

式中：hpq為p，q兩相的相間傳熱系數(shù)。

6）冷凝換熱關(guān)系式

本文利用試驗數(shù)據(jù)擬合得到的冷凝換熱試驗關(guān)系式為

其中，

式中：cp，w為水的定壓比熱，kJ/（kg·℃）；hfg為蒸汽潛熱；ΔT為冷凝前后水的溫差。

通過數(shù)據(jù)對比可知，冷凝換熱系數(shù)的計算結(jié)果（式（16））與實驗結(jié)果的最大偏差在±12%以內(nèi)。

2 幾何模型及仿真模型驗證

2.1 實驗裝置

在前期研究工作中，為分析蒸汽浸沒射流與凝結(jié)換熱現(xiàn)象，搭建了如圖1所示的實驗系統(tǒng)，該系統(tǒng)主要包括電加熱鍋爐、蒸汽噴嘴、水箱、配置了壓力和溫度探針的三維測量支架、可視化系統(tǒng)和數(shù)據(jù)采集系統(tǒng)［8］。

圖1 蒸汽浸沒射流實驗系統(tǒng)Fig.1 Experiment system for steam jet in water

通過調(diào)節(jié)鍋爐功率，即可獲得不同參數(shù)的蒸汽；基于水箱的充水排水系統(tǒng)，即可調(diào)節(jié)水箱內(nèi)部的溫度；利用高速攝影儀，即可對不同工況下的射流冷凝現(xiàn)象進行可視化攝像。該系統(tǒng)的實驗參數(shù)如表1所示。

表1 實驗參數(shù)Table 1 Experimental parameters

2.2 計算模型

結(jié)合實驗過程，將實驗裝置簡化為如圖2所示的計算模型。其中，噴管的上端進口設(shè)置為給定壓力邊界的壓力進口，水箱的水面設(shè)置為壓力出口邊界（壓力設(shè)為1個標準大氣壓），水箱的其余壁面則設(shè)置為壁面無滑移邊界條件。

水箱的幾何尺寸（長×寬×高）為1.0 m×0.5 m×2.0 m，噴管直徑為10 mm，浸沒深度為0.5 m。針對簡化的三維幾何模型，本文采用了結(jié)構(gòu)化網(wǎng)格劃分方法，其空間離散網(wǎng)格如圖3所示，其中噴管中心部分采用O-Block技術(shù)進行了加密處理。

圖2 計算模型幾何結(jié)構(gòu)示意圖Fig.2 Schematic of geometry structure of calculation model

圖3 網(wǎng)格劃分示意圖Fig.3 Display of mesh generation

基于該網(wǎng)格拓撲結(jié)構(gòu)，本文將針對4套不同基準尺度的離散網(wǎng)格（基準尺寸為10～40 mm），在入口壓力為0.2 MPa、冷卻劑溫度為25℃的條件下進行蒸汽射流仿真計算，以驗證網(wǎng)格的獨立性。在4套計算網(wǎng)格下，噴管內(nèi)蒸汽的最高射流流速如圖4所示。當網(wǎng)格尺度低于30 mm時，蒸汽射流流速的最大偏差小于0.5%，從而驗證了對應(yīng)網(wǎng)格的獨立性。因此，本文將網(wǎng)格尺度設(shè)為20 mm，而網(wǎng)格整體規(guī)模約為8×105。

圖4 網(wǎng)格無關(guān)性驗證Fig.4 Verification of grid independence

蒸汽浸沒射流的直接接觸冷凝是一個瞬態(tài)的物理過程，為保證仿真模擬的計算精度，本文將采用雙精度求解器。同時，仿真計算將考慮采用兩相相間滑移的Mixture兩相流模型，而相間傳熱傳質(zhì)將通過在兩相的體積源項中使用DEFINE_PROPERTY UDF宏添加udf進行自定義。進行壓力—速度的耦合計算時，將采用壓力的隱式算子分割算法（Pressure Implicit with Splitting of Operators，PISO）；守恒方程的離散格式中，體積分數(shù)項將采用QUICK格式，其余項采用二階迎風格式。瞬態(tài)計算中：時間步長為10-5s，計算的物理時間為100 ms，每一個計算時間步內(nèi)需要完成50～150次迭代；采用ANSYS Fluent 14.0仿真程序，服務(wù)器的基本配置為 E5 2685 v3（×2）CPU。在此基礎(chǔ)上，完成一個算例約耗時96 h。

本文需保證連續(xù)性方程的收斂殘差降至10-3量級，動量方程收斂至10-6量級以下，同時需監(jiān)控蒸汽排放管出口速度在每一個瞬態(tài)時間步中迭代的變化，并以此作為下一時間步中迭代計算的輸入。綜合上述收斂評判標準，即可保證在一個瞬態(tài)計算步的迭代滿足收斂精度。

2.3 仿真模型驗證

在蒸汽浸沒射流凝結(jié)過程中，蒸汽由噴管射入水中，典型的穩(wěn)定汽羽形狀主要包括錐形和橢球形2種，如圖5所示。其中，圖5（a）～圖5（d）為豎直噴射時不同形狀汽羽的仿真結(jié)果，圖5（e）～圖5（h）為對應(yīng)的實驗結(jié)果［8］。

在低蒸汽流量、低凝水溫度條件下，蒸汽在噴管出口處充分膨脹，然后進入凝水主流中得到充分冷凝，整個汽羽呈現(xiàn)出如圖5（a）所示的錐形汽羽。隨著蒸汽流量增加，汽羽核心區(qū)的蒸汽無法充分冷凝，故有所膨脹，從而形成如圖5（b）所示的橢球形汽羽。當蒸汽流量進一步增加時，汽羽在凝水中膨脹后收縮，但未能及時冷卻的高溫汽羽將進一步膨脹，形成如圖5（c）所示的雙膨脹收縮形汽羽。在低蒸汽流量、高凝水溫度條件下，由于氣泡產(chǎn)生和脫落，噴射蒸汽將在凝水中呈現(xiàn)如圖5（d）所示的振蕩流形。

由圖5可知，兩熱阻模型作為相變模型，能夠較好地模擬不同類型的蒸汽射流凝結(jié)實驗，且仿真計算結(jié)果與實驗結(jié)果［8］吻合度較高，驗證了該模型的適用性及準確性。

3 射流汽羽影響分析

基于上述仿真模型，本文將進一步分析冷凝溫度、入口壓力等因素對當前實驗條件下浸沒射流汽羽行為的影響。

圖5 仿真結(jié)果與實驗結(jié)果對比Fig.5 Comparison between simulation results and experimental results

3.1 汽羽長度確定

一般而言，在實驗條件下可以通過調(diào)節(jié)可視化圖像中氣、液兩相之間的對比度來確定汽羽長度［9］。以噴管直徑10 mm、浸沒深度0.5 m、水溫25℃、入口壓力0.5 MPa條件下的計算結(jié)果為例，參照文獻［10］截取計算結(jié)果中的汽羽，如圖6所示，橫坐標為試驗裝置中水箱液位的不同高度。本文將以汽相體積份額為0.95作為汽羽長度度量的截斷值，如圖6所示，即以兩相區(qū)域內(nèi)汽相份額為0.95的輪廓線作為形成汽羽的邊界。

3.2 汽羽形狀對比

圖6 汽羽長度的計算方法Fig.6 Calculation method of steam plume length

根據(jù)唐思等［11］的研究，入口壓力和水溫將對蒸汽排入過冷水中的膨脹—壓縮過程及蒸汽的冷凝效率產(chǎn)生影響，從而導致蒸汽汽羽形狀和流場有所區(qū)別。如圖7（a）和圖7（b）所示，蒸汽入口壓力從0.3 MPa變化到0.8 MPa，當水溫為25℃和35℃時，汽羽形狀呈現(xiàn)出與實驗結(jié)果相似的橢球形。

隨著蒸汽入口處壓力的增加，汽羽的膨脹程度也隨之增加，具體表現(xiàn)為：噴管出口處的汽羽直徑增加，蒸汽的射流速度增加，并導致汽羽長度進一步增加。如圖7（c）所示，在水溫為50℃條件下，當入口壓力由0.3 MPa增加至0.8 MPa時，汽羽形狀逐漸由橢球形轉(zhuǎn)變?yōu)殡p膨脹收縮形。隨著入口壓力的增加，在汽羽內(nèi)部首先發(fā)生第1次膨脹—壓縮過程，但蒸汽凝結(jié)不充分，進而發(fā)生第2次膨脹—壓縮過程，最終呈現(xiàn)出雙膨脹收縮型汽羽。

圖7 不同溫度、壓力條件下的汽羽形狀對比Fig.7 Comparison of steam plumes under different pressures and temperatures

3.3 溫度場分布

徑向溫度場分布是蒸汽射流凝結(jié)的一個重要特征。圖8所示為入口蒸汽壓力為0.5 MPa、冷凝水溫為25℃的工況下，在噴嘴出口不同距離處截面的徑向溫度分布情況。圖8中，x為距離噴管中心軸線的徑向距離，x=0即為噴管軸線處。

如圖8（a）所示，噴管出口截面上的蒸汽溫度分布與管內(nèi)相似，其中管壁位置出現(xiàn)了較高的汽相溫度分布。這是由于蒸汽噴射到凝水界面時與凝水產(chǎn)生劇烈壓縮，從而導致了對應(yīng)位置處的蒸汽溫度上升。

如圖8（b）所示，汽羽核心處存在低溫區(qū)，由徑向溫度分布曲線可知，蒸汽由于膨脹做功而導致溫度降低。汽羽內(nèi)部沿著徑向溫度是先上升后下降，并逐漸接近凝水溫度。由可壓縮氣體的紊流射流可知，靠近氣液界面處的蒸汽首先處于壓縮狀態(tài)，蒸汽的內(nèi)能增加，溫度升高；然后，由于卷吸和凝結(jié)作用，蒸汽溫度逐漸降低，并接近周圍過冷水的溫度。

圖8 溫度場分布Fig.8 Distribution of temperature field

在圖8（c）和圖8（d）中，由于汽相受到壓縮，且對應(yīng)位置在壓縮波處，因此均在汽羽軸線達到了峰值溫度。

3.4 壓力場分布

圖9所示為10 mm管徑下水溫為50℃時的壓力分布情況。圖9（a）為蒸汽射流流場中心截面的壓力分布，可以發(fā)現(xiàn)，蒸汽射流流場中壓力場的分布大致可以分為3個區(qū)：對應(yīng)汽羽內(nèi)部汽空間范圍為壓力核心區(qū)，其壓力梯度較大；對應(yīng)包裹汽羽汽空間的氣液混合區(qū)為壓力擴散區(qū)，其表壓為負值；對應(yīng)汽羽外側(cè)的凝水區(qū)為壓力邊緣區(qū)，其接近于環(huán)境壓力。在壓力擴散區(qū)，由于蒸汽出現(xiàn)了部分冷凝，故其壓力低于凝水的環(huán)境壓力；在接近液體單相區(qū)，其壓力逐漸接近環(huán)境壓力。

圖9（b）為不同蒸汽壓力進口條件下，距離噴嘴出口0.1 m處的徑向壓力分布圖。由中心位置沿徑向方向，壓力呈現(xiàn)出先降低后升高，然后逐漸接近周圍環(huán)境壓力的變化趨勢。

圖9（b）中的壓力谷值主要源自蒸汽射流的卷吸過程，卷吸過程會強化蒸汽的凝結(jié)作用，導致局部壓力降低。結(jié)合圖9（a）和圖9（b）可知，壓力場沿噴管軸線呈軸對稱形式，且在中軸線處達到壓力峰值。隨著入口初始壓力的增加，該峰值也隨之增加，且分布曲線的影響范圍逐漸擴大，表現(xiàn)為蒸汽汽羽直徑逐漸增大。另一方面，軸線處的壓力峰值相對于蒸汽入口壓力的壓差較大，表明蒸汽射流在凝水中的傳質(zhì)過程較為強烈，從而導致動量傳遞也較強烈。

圖9（c）和圖9（d）所示為不同蒸汽壓力進口條件下，沿噴管軸線方向的壓力分布曲線?？梢钥闯?，汽羽內(nèi)部存在一定范圍的壓力波動，表現(xiàn)為壓力變化呈現(xiàn)出2個峰值和2～3個谷值。這一結(jié)果表明了汽羽內(nèi)部壓縮波和膨脹波的存在性：首先，蒸汽在噴管出口處發(fā)生膨脹過程，導致壓力過低；然后，進入壓縮過程，導致壓力升高；最后，進入壓縮—膨脹過程。隨著入口壓力的增加，第2次壓縮—膨脹過程產(chǎn)生的壓力波動逐漸增強。同時，通過對比圖9（c）和圖9（d）可知，進口壓力對汽羽內(nèi)部壓縮—膨脹波的波節(jié)數(shù)影響較小，但較大的進口壓力將導致壓縮—膨脹波的波長增加，即壓力波動將傳播更遠，其幅值也將有所增加。

圖9 10 mm管徑下水溫為50℃時的壓力分布Fig.9 Pressure distribution under 10 mm tube and 50℃condition

4 結(jié) 論

基于筆者前期開展的浸沒射流實驗，本文針對浸沒射流冷凝現(xiàn)象開展了CFD仿真分析，通過引入實驗過程中擬合的傳熱經(jīng)驗關(guān)系式，得到了與實驗結(jié)果高度吻合的仿真結(jié)果。本文以汽相體積份額為0.95作為度量汽羽長度的臨界分數(shù)，對比研究了不同凝水溫度和進口壓力條件下浸沒射流冷凝汽羽的空間分布、速度場、溫度場及壓力場，充分掌握了這類開放空間純蒸汽浸沒射流的基礎(chǔ)物理現(xiàn)象，得到結(jié)論如下：

1）隨著蒸汽入口壓力的增加，蒸汽在水域中的膨脹程度隨之增加，產(chǎn)生的汽羽直徑和汽羽長度也逐漸增加。同時，不同溫度條件下，冷凝水箱內(nèi)形成的蒸汽汽羽形狀將有所不同。

2）蒸汽進口壓力較高時，水溫對汽羽內(nèi)部蒸汽膨脹的影響較小，且不同水溫條件下沿軸線的汽羽速度分布比較接近。

3）汽羽內(nèi)部的徑向溫度分布將呈現(xiàn)出不同的特征，由于汽羽的膨脹過程，在核心區(qū)將出現(xiàn)局部低溫現(xiàn)象。

4）蒸汽射流氣羽周圍的壓力場分布主要分為3個區(qū)：較大壓力梯度的核心區(qū)、負壓的擴散區(qū)及接近環(huán)境壓力的邊緣區(qū)。

5）隨著蒸汽初始壓力的增加，核心區(qū)靜壓的峰值逐漸升高，蒸汽的汽羽直徑也隨之增加。

6）汽羽內(nèi)部沿軸向的壓力分布一般存在2組波峰—波谷，這是因為汽羽內(nèi)部存在壓縮波和膨脹波，蒸汽在噴嘴出口處先發(fā)生膨脹過程，再進行壓縮過程。進口壓力對汽羽內(nèi)部壓縮—膨脹波的波節(jié)數(shù)影響較小，但會導致壓力波的波動幅值和影響范圍有所增加。