貝克利數(shù)對(duì)圓管流道交變流動(dòng)傳熱的影響分析

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

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

(2中國(guó)科學(xué)院理化技術(shù)研究所低溫工程學(xué)重點(diǎn)實(shí)驗(yàn)室 北京 100190)

1 引言

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

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

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

2 控制方程

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

連續(xù)性方程:

動(dòng)量方程:

能量方程:

3 計(jì)算模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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 湯 珂，張 玙，林小鋼，等，交變流動(dòng)傳熱的相似準(zhǔn)則［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.