湯 珂 唐文濤 張 玙 金 滔 林小鋼 雷 田

(1浙江大學(xué)制冷與低溫研究所 杭州 310027)

(2中國科學(xué)院理化技術(shù)研究所低溫工程學(xué)重點實驗室 北京 100190)

1 引言

回?zé)崾降蜏刂评錂C(例如:斯特林制冷機、GM制冷機、脈管制冷機等)以流體工質(zhì)交變流動為基本特征。在交變流動中，流體工質(zhì)反復(fù)流經(jīng)傳熱表面與之換熱，與通常單一方向流動的傳熱特性顯著不同。交變流流動的傳熱特性是深入理解回?zé)崾降蜏刂评錂C工作原理，以及實現(xiàn)其精確設(shè)計的重要基礎(chǔ)。而目前，對于管內(nèi)交變流動傳熱的研究不足。

美國麻省理工學(xué)院Dean，Smith等人針對斯特林熱機中出現(xiàn)的交變流動及傳熱現(xiàn)象，開展了管內(nèi)交變流動傳熱的實驗研究?？疾炝私蛔冾l率、測試管道兩端所連接的壓縮機活塞之間的相位角對于傳熱的影響，給出第一階和第二階努塞爾數(shù)與貝克利數(shù)的經(jīng)驗關(guān)聯(lián)式［1－3］?？死锓蛱m州立大學(xué) Kannapareddy 等針對斯特林熱機中回?zé)崞鳌⒗鋮s器以及加熱器內(nèi)的交變流動傳熱現(xiàn)象進行了數(shù)值模擬研究［4］。研究結(jié)果中，以管壁溫度和流體混合平均溫度的差值作為溫差來求得無量綱換熱系數(shù)努塞爾數(shù)，在流體流動反向的時刻會出現(xiàn)努塞爾數(shù)趨于無窮大的情況。香港科技大學(xué)Zhao等研究圓管內(nèi)不可壓縮層流交變流動傳熱，采用時空平均努塞爾數(shù)來表征傳熱的結(jié)果，主要考察無量綱流體位移振幅以及動態(tài)雷諾數(shù)對于時空平均努塞爾數(shù)的影響，并給出了相應(yīng)的經(jīng)驗關(guān)聯(lián)式［5］。

目前，對于管內(nèi)交變流動和傳熱特性的研究，主要是針對不可壓縮流體［4－6］。而以色列理工學(xué)院Garaway等的研究指出氣體工質(zhì)的可壓縮性將顯著影響壓力、速度間的相位角與交變頻率間的相關(guān)性［7］。為此，本文通過數(shù)值模擬研究圓管內(nèi)可壓縮層流交變流動傳熱現(xiàn)象。采用離散傅立葉級數(shù)的表達形式表征管內(nèi)交變流動傳熱特性。根據(jù)前人的研究工作［1－3］以及對管內(nèi)可壓縮交變流動傳熱相似準則體系的分析［8］可知，貝克利數(shù)其中:ω為角頻率，Dh為管道水力直徑，α為熱擴散系數(shù)。)是描述交變流動傳熱的關(guān)鍵相似準則。因此本文重點考察貝克利數(shù)對于管內(nèi)可壓縮層流交變流動傳熱特性的影響，計算分析了0階努塞爾數(shù)(即時均努塞爾數(shù))隨貝克利數(shù)的變化關(guān)系。

2 控制方程

可壓縮黏性流體非定常流動(忽略體積力、無內(nèi)熱源)的一般控制方程［9］如下:

連續(xù)性方程:

動量方程:

能量方程:

3 計算模型

計算模型為內(nèi)徑1 mm的圓管流道。文獻［10］等針對管內(nèi)不可壓縮層流交變流，推導(dǎo)了熱發(fā)展段的長度Ld的表達式。根據(jù)該研究結(jié)果，在計算模型中，將管長選定為3倍熱發(fā)展段長度，以使得所關(guān)注的管道軸向中心位置處于充分發(fā)展區(qū)。

流道左端設(shè)為隨時間做正弦變化的質(zhì)量流入口，溫度恒定為330 K。右端為隨時間做正弦變化的壓力出口，操作壓力為2.0 MPa，流體工質(zhì)為氦氣，溫度恒定為330 K。通過調(diào)節(jié)質(zhì)量流入口和壓力出口的相位和幅值使得觀測點(流道軸向中心位置處)的速度和壓力的相位差為0，最大雷諾數(shù)為650，壓比為1.2。管壁溫度保持300 K。采用二維軸對稱計算模型，整個計算區(qū)域軸向網(wǎng)格情況為50/mm，徑向網(wǎng)格數(shù)為50，壁面附近的流體網(wǎng)格加密。時間步長取為交變周期的1/180。利用FLUENT軟件進行求解，壓力/速度耦合采用SIMPLE算法，時間迭代采用一階隱式算法。

根據(jù)脈管制冷機、熱聲熱機等回?zé)崾綗釞C換熱器內(nèi)的實際工況，確定貝克利數(shù)的研究范圍為2到300。研究對象是管內(nèi)層流交變流，依據(jù)Hino等對于管內(nèi)交變流湍流轉(zhuǎn)變的研究結(jié)果［11］，將最大雷諾數(shù)(Remax=umaxDh/ν，umax為最大速度，對于正弦規(guī)律變化情況即為速度振幅，ν為運動粘度)選定為650，以保證流動狀態(tài)為層流。

對模擬計算進行了穩(wěn)定性判定。當(dāng)計算進行到第11個周期及以后，流道中心點速度振幅、溫度振幅及壓力振幅的相對變化率小于1%，認為數(shù)據(jù)已基本穩(wěn)定，取第11個周期的數(shù)據(jù)。同樣針對網(wǎng)格獨立性也進行了相應(yīng)的考察。對于計算中所采用的網(wǎng)格，分別將網(wǎng)格在徑向和軸向加倍后，監(jiān)測研究中所關(guān)心的物理量，如流道中心點速度振幅、壓力振幅、溫度振幅等參數(shù)，各物理量的相對變化率均小于3%。

4 計算結(jié)果與討論

4.1 管內(nèi)交變流動傳熱性能的表征

管內(nèi)交變流動傳熱過程中，流道中心流體與管壁的溫差以及壁面處的熱流密度均為周期變化，如果簡單地根據(jù)牛頓冷卻定律將壁面處的熱流密度除以中心流體與管壁的溫差來獲取瞬時的傳熱系數(shù)可能出現(xiàn)由于中心流體與管壁的溫差趨近于零而導(dǎo)致?lián)Q熱系數(shù)無窮大的情況。關(guān)于管內(nèi)交變流動傳熱性能的表征，香港科技大學(xué)Zhao等采用時空平均努塞爾數(shù)，其中流體的溫度采用測試管段兩端的混合室中氣體的周期平均溫度，而壁面溫度取4個軸向位置測點溫度時均值的平均值，無法反映周期內(nèi)的變化情況［12］。法國南特大學(xué)Bouvier等通過各階諧次累加的熱流密度和壁面溫度與截面平均溫度的差值來計算的無量綱換熱系數(shù)［13］，能夠反映周期內(nèi)的變化情況，但測量截面流體平均溫度及基于反傳熱方法獲得壁溫等都非常困難，精度難以保證。

參考麻省理工學(xué)院Dean，Smith等的研究思路，采用多階努塞爾數(shù)表征管內(nèi)交變流動傳熱性能［1－3］。具體而言，將模擬計算獲得的流道中心流體與管壁的溫差和壁面處的熱流密度的周期變化均表達為離散傅立葉級數(shù)形式，根據(jù)牛頓冷卻定律得到各階對流換熱系數(shù)，然后將其表達為無量綱的各階努塞爾數(shù)。其中第0階表示時均傳熱情況，而從一階開始的以后各階共同反映周期內(nèi)的波動情況。例如:首先將流道中心流體與管壁的溫差和壁面處的熱流密度的周期變化表達為離散傅立葉級數(shù)形式，如下所示:

其中:N為一個流動周期內(nèi)流道中心流體與管壁的溫差和壁面處的熱流密度的數(shù)據(jù)點數(shù)，然后求解各階傳熱系數(shù)，即 hm=qm/ΔTm(m=0，1，…，N/2)，最后由 Num=hmDh/k(m=0，1，…，N/2)求得各階努塞爾數(shù)。

圖1和圖2給出了流道中心流體與管壁的溫差和壁面處的熱流密度分別采用三階、四階以及五階級數(shù)表達式擬合的數(shù)據(jù)。采用式(6)計算一個周期內(nèi)擬合值與模擬計算的原始數(shù)據(jù)的相對誤差:

式中:xi為一個周期內(nèi)第i個時刻的級數(shù)擬合值;xri表示由模擬計算得到的一個周期內(nèi)第i個時刻的原始數(shù)據(jù);x—r則表示原始數(shù)據(jù)的周期平均值。通過將各階擬合值與模擬計算所得的原始數(shù)據(jù)對比發(fā)現(xiàn)，當(dāng)采用五階級數(shù)擬合時，溫差和熱流密度的擬合值與計算所得的原始數(shù)據(jù)的相對誤差均可小于10%，能夠較好地反映原始數(shù)據(jù)的變化情況。

圖1 傅立葉級數(shù)擬合值與溫差原始計算數(shù)據(jù)的比較Fig.1 Comparison between the fitting values with discrete Fourier series and the computed temperature difference

圖2 傅立葉級數(shù)擬合值與熱流密度原始計算數(shù)據(jù)的比較Fig.2 Comparison between the fitting values with discrete Fourier series and the computed heat flux

4.2 管內(nèi)交變流動傳熱特性

首先，給出貝克利數(shù)對一個周期內(nèi)交變流動傳熱過程影響的有代表性計算結(jié)果。圖3和圖4給出了貝克利數(shù)等于2時，計算管段軸向中心位置一個流動周期內(nèi)流道截面中心點的靜壓和速度以及流道截面中心點的溫度和壁面處熱流密度的變化情況。圖5和圖6對應(yīng)給出的是貝克利數(shù)為300時的情況。結(jié)果顯示，在兩種貝克利數(shù)下，流道中心點的溫度和壁面處的熱流密度隨時間的變化呈現(xiàn)不同的變化規(guī)律。如圖4所示，當(dāng)貝克利數(shù)較小時，在一個流動周期內(nèi)，流道中心點的溫度和壁面處的熱流密度會出現(xiàn)兩個峰值，與2005年法國南特大學(xué)Bouvier等［13］報道的實驗測量觀察到流體溫度信號2倍于流體運動頻率的情況相符;而當(dāng)貝克利數(shù)較大時，流道中心點的溫度和壁面處的熱流密度則與靜壓和速度是同頻率變化的，如圖6所示?？梢?，貝克利數(shù)對于管內(nèi)交變流動傳熱會存在顯著的影響。貝克利數(shù)的物理意義是:熱擴散的時間尺度與交變周期之比。當(dāng)貝克利數(shù)較小時，熱擴散的時間相對于交變周期來說很短，壁面處的傳熱足以通過徑向擴散影響到流道中心的溫度，因此在一個流動周期中，流體往復(fù)兩次反向的流動傳熱過程導(dǎo)致流道中心處溫度信號出現(xiàn)2倍于交變流動頻率的情況。而當(dāng)貝克利數(shù)較大時，熱擴散的時間相對于交變周期而言很長，也即熱量還來不及擴散到管道中心處，流體已經(jīng)流過該傳熱面，流體雖然也是兩次流經(jīng)傳熱面，但是流道中心溫度受壁面?zhèn)鳠岬挠绊懶?，其變化主要源于壓力的波動，因此保持與交變流動頻率一致。

圖3 流道中心點的靜壓和速度，Peω=2Fig.3 Pressure and velocity at the channel center，Peω =2

在實際換熱器的設(shè)計中，主要關(guān)注的是周期平均換熱效果，因此進一步考察時均努塞爾數(shù)(即第0階努塞爾數(shù))隨著貝克利數(shù)的改變而呈現(xiàn)的變化規(guī)律。計算結(jié)果示于圖7。從圖中可以看出，在設(shè)定最大雷諾數(shù)為650的條件下，當(dāng)貝克利數(shù)小于20時，時均努塞爾數(shù)基本不隨貝克利數(shù)變化;而當(dāng)貝克利數(shù)在20到300范圍內(nèi)時，時均努塞爾數(shù)則隨貝克利數(shù)增大呈現(xiàn)增大的趨勢?？梢?，當(dāng)貝克利數(shù)較大時，其對交變流動傳熱的影響更為顯著。

圖4 流道中心點的溫度和壁面處熱流密度，Peω=2Fig.4 Temperature at the channel center and heat flux at wall，Peω =2

圖5 流道中心點的靜壓和速度，Peω=300Fig.5 Pressure and velocity at the channel center，Peω =300

圖6 流道中心點的溫度和壁面處熱流密度，Peω=300Fig.6 Temperature at the channel center and heat flux at wall，Peω =300

圖7 時均努塞爾數(shù)隨貝克利數(shù)的變化Fig.7 Time-averaged Nusselt number as a function of Peclet number

5 結(jié)論

根據(jù)可壓縮黏性流體的一般控制方程，通過數(shù)值模擬研究了管內(nèi)交變流動傳熱特性，重點考察貝克利數(shù)對于交變流動傳熱的影響，得出以下結(jié)論:

(1)采用離散傅立葉級數(shù)形式表達流道中心流體與管壁的溫差以及壁面處熱流密度時，取至五階級數(shù)擬合可與原始數(shù)據(jù)的相對誤差小于10%，能夠較好地反映原始數(shù)據(jù)的變化情況。

(2)在本文分析條件下，當(dāng)貝克利數(shù)較小時(如貝克利數(shù)等于2)，在一個流動周期內(nèi)，流道中心點的溫度和壁面處的熱流密度均會出現(xiàn)2個峰值，即這2個參數(shù)的變化2倍于交變流動的頻率;當(dāng)貝克利數(shù)較大時(如貝克利數(shù)等于300)，在一個流動周期內(nèi)，流道中心點的溫度和壁面處的熱流密度則與流體的靜壓和速度保持同頻率變化。

(3)在本文分析條件下，當(dāng)貝克利數(shù)小于20時，時均努塞爾數(shù)幾乎不隨貝克利數(shù)變化;而當(dāng)貝克利數(shù)在20到300的范圍內(nèi)，時均努塞爾數(shù)則隨著貝克利數(shù)的增長而增大，即在貝克利數(shù)較大時，其對交變流動傳熱的影響更為顯著。

1 Dean C E.Establishing baseline data for an experimental apparatus that measures heat transfer under conditions of oscillating pressure and flow［D］.Cambridge:Massachusetts Institute of Technology，1993.

2 Grassmyer W A.Experiments in heat transfer under conditions of oscillating pressure and flow［D］.Cambridge:Massachusetts Institute of Technology，1994.

3 Cheng P D L.Experiments in heat transfer under conditions of oscillating pressure and flow［D］.Cambridge:Massachusetts Institute of Technology，1995.

4 Kannapareddy M R.Numerical thermal analysis of heat exchangers for the Stirling engine application［R］.NASA－CR－199120，1995.

5 Zhao T S，Cheng P.A numerical solution of laminar forced convection in a heated pipe subjected to a reciprocating flow ［ J］.International Journal of Heat and Mass Transfer，1995，38(16):3011－3022.

6 Jaworski A J，Mao X A，Mao X R，et al.Entrance effects in the channels of the parallel plate stack in oscillatory flow conditions［J］.Experimental Thermal and Fluid Science，2009，33(3):495－502.

7 Garaway I，Grossman G.Studies in high frequency oscillating compressible flow for application in a micro regenerative cryocooler［C］.Advances in Cryogenic Engineering:Transactions of the Cryogenic Engineering Conference，2006，823:1588－1595.

8 湯 珂，張 玙，林小鋼，等，交變流動傳熱的相似準則［J］，低溫工程，2010(2):4－7.

9 Anderson J D.Computational Fluid Dynamics［M］.New York:McGraw－Hill Inc，1995.

10 Lee D Y，Park S J，Ro S T.Heat transfer in the thermally developing region of a laminar oscillating pipe flow ［ J］.Cryogenics，1998，38(6):585－594.

11 Hino M，Sawamoto M，Takasu S J.Experiments on transition to turbulence in an oscillatory pipe flow ［ J］.Journal of Fluid Mechanics，1976，75(2):193－207.

12 Zhao T S，Cheng P.Oscillatory heat transfer in a pipe subjected to a laminar reciprocating flow ［ J］.Journal of Heat Transfer，1996，118(3):592－597.

13 Bouvier P，Stouffs P，Bardon J P.Experimental study of heat transfer in oscillating flow ［ J］.International Journal of Heat and Mass Transfer，2005，48(12):2473－2482.