非絕熱單孔容腔瞬態(tài)響應的零維建模

丁水汀，于航，邱天，*

（1.飛機/發(fā)動機綜合系統(tǒng)安全性北京市重點實驗室，北京 100083； 2.先進航空發(fā)動機協(xié)同創(chuàng)新中心，北京 100083；3.北京航空航天大學 能源與動力工程學院，北京 100083）

隨著現(xiàn)代航空發(fā)動機不斷發(fā)展，對高機動性的需求在不斷提高［1］，發(fā)動機經歷從慢車到最大轉速的時間縮短，其表現(xiàn)出來的瞬變特征也更加顯著，這加劇了發(fā)動機狀態(tài)瞬變演化規(guī)律的復雜程度，發(fā)動機空氣系統(tǒng)氣路參數(shù)也將經歷復雜的瞬變過程［2］?？諝庀到y(tǒng)在航空發(fā)動機中發(fā)揮著冷卻、防冰、封嚴和軸向力控制等重要作用，掌握空氣系統(tǒng)瞬變演化規(guī)律極為重要。容腔類元件屬于航空發(fā)動機空氣系統(tǒng)中一種重要元件，在過渡過程中對于氣體的儲存和釋放作用使得空氣系統(tǒng)內部各元件響應速度有著不同程度的滯后，直接影響空氣系統(tǒng)瞬態(tài)分析，因此在強瞬變過程中要重點考慮容腔。

Dutton和 Coverdill［3］對容腔充放氣過程進行了試驗研究，對比了在不同的容腔體積和進出口噴嘴面積下的試驗數(shù)據(jù)與容腔絕熱模型計算結果，并分析了產生差異的原因。Thorncroft等［4］通過試驗修正容腔內氣體多變過程的多變指數(shù)來實現(xiàn)容腔內壓力、溫度響應規(guī)律與試驗數(shù)據(jù)吻合，說明容腔內氣體歷經的是非確定過程，既不是等溫過程也不是絕熱過程，但未研究氣體與容腔壁面換熱的影響作用。文獻［5-7］則主要分析容腔壁面帶有離散熱源的封閉容腔內的自然對流過程中的流動與換熱特征，缺乏對于帶有換熱過程的容腔充放氣過程的研究。Owczarek和 Progelhof［8］對帶有噴嘴的圓柱筒容器快速放氣過程進行了試驗研究，表明在不考慮容器內換熱和耗散的情況下，試驗和計算結果差異主要來自噴嘴內邊界層的影響，說明噴嘴對于容腔的瞬態(tài)響應也有影響。Chow等［9］應用了速矢變換的方法研究了帶有軸向對稱可控閥門的容腔放氣問題。上述對于容腔的理論和試驗研究中沒有考慮容腔的熱邊界，缺乏對換熱影響的分析，未針對容腔建立真實響應模型。傳統(tǒng)容腔建模方法一般采用絕熱假設，然而在航空發(fā)動機中不存在完全絕熱的容腔，而且數(shù)據(jù)表明容腔換熱使得充放氣響應時間較絕熱過程增加［3］，說明容腔換熱的影響顯著，需要在發(fā)動機瞬態(tài)空氣系統(tǒng)分析中充分考慮，但目前用于分析瞬態(tài)空氣系統(tǒng)的流體網絡法中缺少容腔換熱模型的支撐。

在航空發(fā)動機中存在多種通氣形式的容腔，如單孔、單進單出、單進多出和多進多出等。其中，單孔形式的容腔在航空發(fā)動機 CFM56、D30KY-154等型號中均存在［10］，單孔作為卸壓孔，可防止容腔內壓力變化劇烈造成周圍結構失穩(wěn)。單孔形式容腔無論對發(fā)動機空氣系統(tǒng)瞬態(tài)演化過程，還是在維持結構穩(wěn)定性方面，都起到很大作用。此外，單孔容腔結構簡單，可在此基礎上對其他形式容腔展開研究。本文重點研究單孔形式的容腔，建立非絕熱單孔容腔零維瞬態(tài)模型。

1 理論分析

1.1 基本假設

Horlock和 Woods［11］全面分析了容腔充放氣過程的熱力學特征，如果容腔充放氣過程迅速，可近似認為是絕熱過程，緩慢則可近似認為是等溫過程。在實際過程中，絕熱模型或等溫模型都不能準確地描述腔內氣體的變化過程。

在容腔充放氣的過程中，實際腔內的壓力和溫度分布不均勻，但是容腔三維特征并不明顯，本文針對非絕熱單孔容腔充放氣過程建模，提出如下假設：

1）容腔時間尺度參數(shù)滿足 τtotal?tchar，其中τtotal為總充放氣時長；tchar為特征時間，表示壓力波傳遍整個容腔所需時長，定義為tchar＝l/a，l為容腔特征尺度（對于方體容腔為邊長），a為容腔內氣體聲速。當滿足τtotal?tchar時，可認為容腔發(fā)生質量存儲的時間尺度遠大于腔內建立穩(wěn)定流場的時間尺度，即在任一時刻下，腔內流動與換熱僅與當前時刻進口氣體流動參數(shù)和腔內氣體狀態(tài)有關，與瞬變過程無關。

2）腔內氣體狀態(tài)可集總假設，即假設在瞬變過程中腔內氣體速度足夠低，容腔內氣體參數(shù)分布足夠均勻，容腔狀態(tài)只與時間有關，而與空間無關，可用平均狀態(tài)參數(shù)表示［12］。

3）腔內氣體常物性假設，即容腔內氣體物性不隨氣體熱力狀態(tài)變化［12］。

1.2 控制方程

對于理想氣體從上游總壓到下游靜壓的一維等熵可壓縮膨脹過程，理論流量可以表示為［13］

式中：W為氣體流量；A為流通截面積；ρ為氣體密度；p*和p分別為氣體上游總壓和下游靜壓；κ為絕熱指數(shù)。式（1）適用范圍為亞聲速流動工況。

非定常流動的質量方程［14］為

式中：Wi為容腔進出口氣體流量；Hi為進出口氣體焓；對于本文研究的單孔容腔，n＝1為容腔體積；e為氣體內能；Q為氣體與容腔壁面對流換熱熱流；t為時間。

本文重點關注如何處理氣體與容腔壁面的對流換熱熱流Q，而傳統(tǒng)容腔建模方法忽略了換熱的影響，對流換熱熱流表達式為

式中：hav為平均對流換熱系數(shù)；Aw為氣體與容腔壁面換熱面積；Tref和Tw分別為氣體參考溫度及壁面溫度，這里氣體參考溫度為腔內氣體平均溫度。換熱控制方程（4）無法表達成解析的形式，因為平均對流換熱系數(shù)hav通常難以解析表達。

為了解決非絕熱單孔容腔零維瞬態(tài)模型中換熱熱流如何處理這一問題，本文綜合考慮了影響容腔換熱的因素，采用量綱分析推導了單孔容腔換熱特性的特征數(shù)方程，進而確定熱流的大小。

1.3 換熱特性

非絕熱單孔容腔零維瞬態(tài)建模重點是獲得壁面對流換熱特性，其準確與否直接影響模型的精確性。氣體與容腔壁面換熱和多個參數(shù)相關，為了將高維問題簡化，本文針對單孔的容腔形式，采用無量綱分析法［15］推導單孔容腔換熱的相似準則及準則方程：

1）影響單孔容腔壁面對流換熱系數(shù)h的主要獨立參數(shù)有：容腔內氣體的物性參數(shù)（導熱系數(shù)k、動力黏度 μ、氣體密度 ρ、定壓比熱容cp）、孔徑d（特征長度）、壁面處氣流速度u。依據(jù)假設認為容腔內建立穩(wěn)定流場的時間尺度足夠小，則壁面處氣流速度u與進氣速度ui是一一對應的，氣體密度ρ是容腔內壓力Pv和溫度Tv的函數(shù)，同時Pv受進氣流量 ρi Aui影響（ρi為進氣密度）。本模型幾何形式固定，不考慮幾何形狀參數(shù)的影響。綜上所述，影響壁面對流換熱系數(shù)的函數(shù)關系式可以寫成

2）基本量綱：質量 M、長度 L、時間 T和溫度Θ。

3）基本物理量：ui、d、k、μ。

4）各獨立參數(shù)的基本量綱見表1。

5）分別對定壓比熱容、進口密度和對流換熱系數(shù)建立量綱方程，分別得到無量綱數(shù)：普朗特數(shù)Pr＝cpμ/k、進氣雷諾數(shù)Re＝ρi ui d/μ和努賽爾數(shù)Nu＝hd/k，對腔內溫度建立量綱方程得到無量綱數(shù)：速度系數(shù) λ和溫比Tv/T0，其中T0為進口總溫。

至此，式（5）可轉化為

式（6）即為在處理氣體與容腔壁面間的換熱時所采用的換熱特性準則關系式。

針對單孔容腔幾何形式，分析了在瞬態(tài)響應過程中適用的換熱特征數(shù)方程，為研究考慮換熱影響的非絕熱容腔真實響應過程提供基礎。對于單進單出等多種其他幾何形式容腔均可采用量綱分析方法獲取適用不同幾何形式容腔的換熱特征數(shù)方程進行建模。

表1 參數(shù)量綱Table 1 Dimension of parameters

2 建模方法

首先分析換熱對非絕熱單孔容腔瞬態(tài)響應規(guī)律的影響程度，基于此分析討論如何反映換熱影響，顯式表達換熱項，進而引入換熱建立非絕熱單孔容腔零維瞬態(tài)模型，最后通過CFD數(shù)值模擬結果對模型進行對比驗證。

第1節(jié)給出了非絕熱單孔容腔零維瞬態(tài)建模的理論分析，明確了建模需求，通過量綱分析推導出與換熱特性相關聯(lián)的準則數(shù)，為建模提供了基礎。換熱特性作為非絕熱單孔容腔零維瞬態(tài)模型的重要輸入，其準確性直接影響模型的準確性。本文從CFD數(shù)值模擬中獲取影響換熱的相關參數(shù)，通過整理得到準則數(shù)并借助非線性擬合出具體的準則關系式，利用該式顯示表達換熱項，聯(lián)立求解質量方程和能量方程得到非絕熱單孔容腔零維瞬態(tài)響應規(guī)律，建模方法如圖1所示。

通過與CFD數(shù)值模擬結果進行對比驗證非絕熱單孔容腔零維瞬態(tài)模型的準確性。CFD數(shù)值模擬本身涵蓋了動量方程和能量方程的交互機理，兩者耦合求解。在處理容腔問題時，在時間的維度上是非穩(wěn)態(tài)的，在空間的維度上是不均勻的，反映了物理本質和容腔內氣體的真實變化，利用CFD數(shù)值模擬可以驗證在本文研究的時間和空間尺度上模型的準確性。

圖1 非絕熱單孔容腔零維瞬態(tài)建模方法Fig.1 Zero-dimensional transientmodeling method of non-adiabatic cavity with single opening

3 模型計算方法及驗證

3.1 計算方法

3.1.1 非絕熱單孔容腔零維瞬態(tài)模型計算

基于第1節(jié)和第2節(jié)分析，在Visual C++環(huán)境下編制計算程序，建立了非絕熱零維瞬態(tài)模型（Non-adiabatic Zero-dimensional Transient Model，NZTM）。NZTM首先需要獲取幾何參數(shù)、邊界條件、時間步長和模擬總時長，根據(jù)初始條件，在不同時刻通過迭代求解控制方程組得到當前時刻容腔內平均壓力和溫度的響應規(guī)律。NZTM具體計算流程如圖2所示。

為建立NZTM，需要在能量方程中添加氣體與容腔壁面間的對流換熱項。通過CFD數(shù)值模擬獲取與換熱相關的物理參數(shù)，依據(jù)本文中給出的處理方式整理出換熱特性的具體表達形式，從而支撐模型的建立。

圖2 非絕熱零維瞬態(tài)模型計算流程Fig.2 Calculation flowchart of non-adiabatic zero-dimensional transient model

3.1.2 CFD數(shù)值模擬

CFD數(shù)值模擬一方面提供與換熱相關的準則數(shù)數(shù)據(jù)，一方面用于驗證模型的準確性。本文針對單孔容腔的充氣過程進行驗證，因為在該過程中進口氣體對腔內氣體具有沖擊作用，換熱較放氣過程更明顯，選擇充氣過程更有利于驗證非絕熱單孔容腔零維瞬態(tài)模型的準確性。

容腔物理模型是邊長為0.5 m的立方體容腔，壁面邊界厚度為 0，容腔一側壁面有直徑0.025m的孔，如圖3所示。

圖3 容腔物理模型Fig.3 Physical model of cavity

采用Fluent計算，模型網格劃分如圖4所示，為六面體結構化網格，數(shù)量20萬，計算采用基于密度的算法進行瞬態(tài)求解，控制方程在網格系統(tǒng)上采用有限體積法離散，對流項采取二階迎風格式，時間項采取一階隱式差分格式，湍流模型為標準k-ε湍流模型。初始邊界（時刻t＝0）和計算邊界（時刻t＞0）條件參數(shù)設置見表2。

數(shù)值計算邊界進口為壓力進口，壓力從101 325 Pa階躍至191 504 Pa，容腔壁面為等壁溫條件。

圖4 數(shù)值計算模型網格劃分Fig.4 Meshing of numerical calculation model

表2 初始邊界和計算邊界條件參數(shù)Tab le 2 Initial boundary and computational boundary condition parameters

由特征數(shù)方程式（6）可知，直接反映換熱特性的努賽爾數(shù)Nu與4個準則數(shù)相關，在一定條件下可以對該關系式進行簡化。

氣體密度的靜參數(shù)和總參數(shù)之比的氣動函數(shù)為

式中：ρ和ρ*分別為密度的靜參數(shù)和總參數(shù)。

絕能等熵滯止條件下的氣體狀態(tài)方程為

式中：T*為總溫；R為氣體常數(shù)。

雷諾數(shù)Re及速度系數(shù)λ的表達式分別為

式中：ccr為臨界聲速。

將式（9）和式（10）整理可得如下關系式：

由上述分析可知，在確定的進口初始邊界（包括進口總壓和總溫）及幾何條件下，雷諾數(shù)與速度系數(shù)關系是一一對應的。同時依據(jù)常物性假設，普朗特數(shù)不變，所以準則關系式（6）可簡化為

將不同雷諾數(shù)Re、溫比Tv/T0和對應的壁面平均努塞爾數(shù)Nu擬合成函數(shù)關系式，各無量綱數(shù)關系及擬合數(shù)據(jù)如圖5所示，擬合數(shù)據(jù)最大誤差不超過2%，擬合得到的特征數(shù)方程為

式（13）即對應該容腔幾何形式和初始邊界條件的換熱特性準則關系式。在不同時刻下，根據(jù)腔內氣體狀態(tài)及進口氣流參數(shù)確定唯一的努賽爾數(shù)。至此，處理了氣體與容腔壁面間的換熱，滿足了建模需求。

3.2 計算結果對比分析

1）將NZTM與CFD數(shù)值模擬以及絕熱模型計算結果進行對比，如圖6所示。

圖5 無量綱數(shù)關系及擬合數(shù)據(jù)Fig.5 Relation of dimensionless number and fitting data

圖6 CFD數(shù)值模擬和NZTM、絕熱模型結果對比Fig.6 Comparison of results calculated by CFD simulation，NZTM and adiabatic model

可以看出：①在相同的邊界條件下，NZTM與CFD數(shù)值模擬結果無論在趨勢上還是數(shù)值上都吻合得很好，兩者最大相對誤差小于0.8%，如圖7所示。②NZTM與CFD數(shù)值模擬結果存在微小偏差，產生偏差的原因：一方面在于NZTM處理換熱特性的方式，即將換熱特性擬合成準則關系式，而數(shù)據(jù)擬合本身存在誤差；另一方面對于容腔的基本假設也與真實情況存在稍許差別。③非絕熱單孔容腔氣體壓力和溫度響應與絕熱模型計算結果存在很大差異，其中壓力響應的最大相對誤差達6%，如圖8所示，這也充分說明了換熱對于容腔的瞬態(tài)響應有重要影響。同時絕熱容腔內氣體達到穩(wěn)定壓力所需要的時間較非絕熱容腔短，這是因為在容腔充氣過程中絕熱容腔內氣體溫度上升更快，在充氣過程完成時腔內壓力相同，則穩(wěn)定后絕熱容腔內氣體密度小，說明絕熱過程充氣量較少，達到穩(wěn)定的時間也較短。從溫度響應的對比圖（見圖6（b））中可以看到非絕熱容腔內的氣體溫度更低，在狀態(tài)穩(wěn)定后，溫度還有下降趨勢，這主要是因為腔內氣體與壁面有溫差，換熱過程還在進行。

2）將NZTM與絕熱模型和等溫模型壓力響應進行對比，如圖9所示。

可以看出：NZTM壓力響應曲線介于絕熱和等溫模型曲線之間。由于響應過程較緩慢，氣體與壁面換熱充分，所以氣體歷經過程更接近等溫，偏離絕熱。三線對比說明了本文建立的NZTM較絕熱和等溫模型能夠更精確地反映單孔容腔內氣體真實瞬態(tài)響應規(guī)律。

圖8 NZTM與絕熱模型壓力的相對誤差Fig.8 Relative error of pressure between NZTM and adiabatic model

圖9 NZTM與等溫、絕熱模型壓力計算結果對比Fig.9 Comparison of pressure calculated by NZTM，isothermal and adiabatic model

4 結 論

容腔內氣體與腔壁換熱對于容腔瞬態(tài)響應有重要影響，絕熱容腔模型不能精確反映容腔真實響應歷程，本文基于此需求建立了非絕熱單孔容腔零維瞬態(tài)模型（NZTM）。

1）NZTM與CFD數(shù)值模擬結果對比最大誤差不超過0.8%（絕熱模型為6%），表明本文提出的NZTM較絕熱模型更好地反映了容腔內氣體狀態(tài)的真實變化。

2）雖然非絕熱單孔容腔的換熱特性是通過數(shù)值模擬建立的，但是對比結果首先驗證了本文對于時間尺度分析的正確性，其次說明了依據(jù)本文處理容腔換熱的方法，可以復現(xiàn)容腔內氣體的真實瞬態(tài)響應，體現(xiàn)了準確的換熱特性可以保證模型的準確性，精確反映容腔真實瞬態(tài)響應。

3）非絕熱零維模型在縮減了3個維度之后，較CFD數(shù)值模擬計算量少了數(shù)萬倍，而精度較高，誤差不超過1%，表明方法是可行且有效的。

（References）

［1］SEHRA A K，JR WHITLOW W.Propulsion and power for 21st century aviation［J］.Progress in Aerospace Sciences，2004，40（4-5）：199-235.

［2］劉傳凱，李圓圓，李艷茹，等.渦輪軸斷裂條件下空氣系統(tǒng)強瞬變過程分析［J］.北京航空航天大學學報，2016，42（1）：47-53.LIU C K，LIY Y，LIY R，et al.Dynamic analysis of air system with fast transients in shaft failure event［J］.Journal of Beijing University of Aeronautics and Astronautics，2016，42（1）：47-53（in Chinese）.

［3］DUTTON JC，COVERDILL R E.Experiments to study the gaseous discharge and filling of vessels［J］.International Journal of Engineering Education，1997，13（2）：123-134.

［4］THORNCROFT G，PATTON J S，GORDON R.Modeling compressible air flow in a charging or discharging vessel and assessment of polytropic exponent［C］∥ASEE Annual Conference.Honolulu：ASEE，2007：1-18.

［5］NARDINIG，PARONCINIM.Heat transfer experiment on natural convection in a square cavity with discrete sources［J］.Heat and Mass Transfer，2012，48（11）：1855-1865.

［6］DENG Q H.Fluid flow and heat transfer characteristics of natural convection in square cavities due to discrete source-sink pairs［J］.International Journal of Heat and Mass Transfer，2008，51（25-26）：5949-5957.

［7］NASSAR N N，MEHROTRA A K.Design of a laboratory experiment on heat transfer in an agitated vessel［J］.Education for Chemical Engineers，2011，6（3）：83-89.

［8］OWCZAREKJA，PROGELHOF R C.Rapid discharge of a gas from a cylindrical vessel through a nozzle［J］.American Institute of Aeronautics and Astronautics，1963，1（9）：2182-2183.

［9］CHOW W L，TING A A L，TSAIP H.Discharge from a vessel through an axisymmetric control valve［J］.Journal of App lied Mechanics，1987，54（2）：447-452.

［10］陳光，洪杰.航空燃氣渦輪發(fā)動機結構［M］.北京：北京航空航天大學出版社，2010：170-171.CHEN G，HONG J.Structure of aircraft gas turbine engine［M］.Beijing：Beihang University Press，2010：170-171（in Chinese）.

［11］HORLOCKJH，WOODSWA.The thermodynamics of charging and discharging processes［C］∥Proceedings of the Institution of Mechanical Engineers，London：IMechE，1965，180（310）：16-24.

［12］劉傳凱，劉海明，李艷茹，等.強瞬變空氣系統(tǒng)的模塊化仿真建模［J］.航空動力學報，2015，30（8）：1826-1833.LIU C K，LIU H M，LIY R，et al.Modularized simulation of air system with fast transient［J］.Journal of Aerospace Power，2015，30（8）：1826-1833（in Chinese）.

［13］潘錦珊，單鵬.氣體動力學基礎［M］.北京：國防工業(yè)出版社，2011：87-90.PAN JS，SHAN P.Fundamentals of gasdynamics［M］.Beijing：National Defense Industry Press，2011：87-90（in Chinese）.

［14］駱廣琦，桑增產，王如根，等.航空燃氣渦輪發(fā)動機數(shù)值仿真［M］.北京：國防工業(yè)出版社，2007：84-86.LUO G Q，SANG Z C，WANG R G，et al.Numerical methods for aviation gas turbine engine simulation［M］.Beijing：National Defense Industry Press，2007：84-86（in Chinese）.

［15］鄒滋祥.相似理論在葉輪機械模型研究中的應用［M］.北京：科學出版社，1984：17-24.ZOU Z X.Application of similarity theory in turbo machinery model［M］.Beijing：Science Press，1984：17-24（in Chinese）.