趙偉濤 吳九匯?

1)(西安交通大學(xué)機(jī)械學(xué)院，西安 710049)

2)(西安交通大學(xué)，機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710049)

(2013年4月2日收到;2013年5月4日收到修改稿)

1 引言

自從17世紀(jì)Fourier建立了導(dǎo)熱的數(shù)學(xué)模型，F(xiàn)ourier定律隨之被廣泛應(yīng)用于導(dǎo)熱問(wèn)題分析的各個(gè)領(lǐng)域.對(duì)于熱作用時(shí)間較長(zhǎng)的穩(wěn)態(tài)傳熱過(guò)程以及熱傳播速度較快的非穩(wěn)態(tài)常規(guī)傳熱過(guò)程，采用Fourier定律來(lái)描述熱流密度與溫度梯度之間的關(guān)系是可以滿足精度要求的.但是，F(xiàn)ourier定律不涉及傳熱時(shí)間項(xiàng)，隱含了熱擾動(dòng)傳播速度為無(wú)限大的假設(shè).隨著科學(xué)技術(shù)的進(jìn)步，超短激光脈沖的出現(xiàn)和制冷水平的提高，存在著極高(低)溫條件下的傳熱問(wèn)題和超急速傳熱問(wèn)題，使得Fourier定律中的準(zhǔn)平衡條件假設(shè)不再成立[1-4].

為了克服Fourier定律的局限性，Cattaneo[5]和Vernotte[6]分別獨(dú)立地提出了具有熱流延遲相的非Fourier熱傳導(dǎo)模型，計(jì)及熱流變化率對(duì)熱傳導(dǎo)的影響，其修正后的熱傳導(dǎo)方程為

其中，q是熱流矢量，T為溫度，k為熱傳導(dǎo)率，t是時(shí)間，τ0為熱松弛時(shí)間.方程(1)結(jié)合能量守恒方程得到以溫度T描述的雙曲線型熱傳導(dǎo)方程：

其中，a=k/(ρc)為熱擴(kuò)散系數(shù)，ρ為密度，c為常應(yīng)變比熱.

為描述熱以有限速度傳播這一超常規(guī)熱傳導(dǎo)現(xiàn)象，研究者們對(duì)方程(2)在不同初始條件和邊界條件下進(jìn)行了求解.Tao等[7]采用數(shù)值模擬研究了固體激光器的非Fourier導(dǎo)熱問(wèn)題.李世榮等[8]采用單相延遲模型，研究了薄板在受周期熱流邊界條件下板內(nèi)的溫度響應(yīng).Sarkar和Haji-Sheikh[9]采用Laplace變換技術(shù)研究了由介電材料制成的有限平板的雙曲型熱傳導(dǎo)問(wèn)題.Tang和Araki[10]求解了表面周期性加熱條件下有限介質(zhì)的非Fourier熱傳導(dǎo)問(wèn)題.Moosaie[11]在Tang等的工作基礎(chǔ)上，運(yùn)用Fourier積分表達(dá)式，求解了有限介質(zhì)在周期表面熱擾動(dòng)條件下的非Fourier熱傳導(dǎo)問(wèn)題.Barletta和Zanchini[12]分析了無(wú)限圓柱體存在內(nèi)熱源以及與外界流體有熱對(duì)流的情況下的雙曲型熱傳導(dǎo)問(wèn)題.Ate fi和Talaee[13]使用分離變量法求解了邊界條件不隨時(shí)間變化的無(wú)限長(zhǎng)圓筒的非Fourier溫度場(chǎng).Mishra和Sahai[14]運(yùn)用格子Boltzmann法研究了一維圓柱和球體的非Fourier熱傳導(dǎo)問(wèn)題.Jiang[15]運(yùn)用Laplace變換法研究了空心球體在內(nèi)外兩個(gè)表面溫度突然變化時(shí)的雙曲型熱傳導(dǎo)問(wèn)題.Shirmohammadi和Moosaie[16]采用分離變量法得到空心球體在周期表面熱流條件下的解析解.

本文主要運(yùn)用Duhamel積分和Fourier級(jí)數(shù)展開(kāi)法得到平板前表面熱流任意周期變化時(shí)雙曲型熱傳導(dǎo)方程的解析解.首先采用Duhamel積分和分離變量法分析了方程在平板前表面遭受突變熱流和簡(jiǎn)諧熱流這兩類特殊情況下解的形式，在此基礎(chǔ)上應(yīng)用Fourier級(jí)數(shù)展開(kāi)法和疊加原理研究了平板前表面熱流任意周期變化時(shí)解的形式.按照這些表達(dá)式，不同周期邊界條件下平板的雙曲線熱傳導(dǎo)行為得到分析和研究.這為工程應(yīng)用和數(shù)值計(jì)算的驗(yàn)證提供了便利.

2 數(shù)學(xué)模型

考慮一厚度為L(zhǎng)，邊界絕緣，熱物性為常數(shù)的平板，假設(shè)其初始溫度T(x,0)=T0，從時(shí)間t=0時(shí)起，平板一表面x=0處(稱該表面為平板前表面，另一面稱平板后表面)遭受一熱流為q0·q(t)的作用.在這種情況下，由方程(2)可得一維平板的非Fourier熱傳導(dǎo)方程：

為了獲得方程(3)—(5)的無(wú)量綱形式，特引入以下無(wú)量綱量：

式中，F(xiàn)o是 Fourier數(shù)，Ve是Vernotte數(shù)，1/Ve代表溫度波傳播的無(wú)量綱速度.

則方程(3)—(5)無(wú)量綱化之后為

3 分析求解

由于邊界條件(8)式中的f(Fo)是任意函數(shù)，這就使得直接求解方程(7)變得不可能.因此，首先假定f(Fo)為一時(shí)間無(wú)關(guān)量f，求解此條件下的溫度場(chǎng);其次運(yùn)用Duhamel積分，求解在任意函數(shù)f(Fo)下的溫度場(chǎng).

設(shè)在邊界條件f下方程(7)的解為

當(dāng)Ve→0時(shí)，方程(7)對(duì)應(yīng)經(jīng)典Fourier導(dǎo)熱的情形.此時(shí)βn為實(shí)數(shù)，(35)式變?yōu)?/p>

把(37)式代入(34)式可得經(jīng)典Fourier導(dǎo)熱時(shí)平板內(nèi)部的溫度場(chǎng).

4 結(jié)果分析與討論

為了得到任意周期邊界條件下平板溫度場(chǎng)的解析解，將做如下遞進(jìn)式推導(dǎo).4.1給出了邊界條件為常數(shù)的情況;4.2的邊界條件為簡(jiǎn)諧周期函數(shù);4.3在前兩步的基礎(chǔ)上，計(jì)算與討論了邊界條件可用Fourier級(jí)數(shù)展開(kāi)的任意周期函數(shù)的情況.

4.1 突變邊界條件

當(dāng)厚度為L(zhǎng)的平板前表面遭受大小為q0的突變熱擾動(dòng)(q(t)=1)時(shí)，無(wú)量綱化后

此時(shí)(30)式就是方程(7)在邊界條件為常數(shù)下的解

對(duì)于Fourier導(dǎo)熱，當(dāng)Ve→0時(shí)，方程的解為

根據(jù)(39)式，可以計(jì)算得到厚度為L(zhǎng)的平板前面遭受一突變熱擾動(dòng)時(shí)的溫度分布.為了驗(yàn)證結(jié)果的可靠性，直接對(duì)方程(7)—(9)進(jìn)行數(shù)值求解.圖1給出了平板前表面遭受突變熱擾動(dòng)時(shí)，在Ve=0.8下，平板后表面溫度響應(yīng)的數(shù)值解與解析解的對(duì)比.可以看出，數(shù)值解和解析解的誤差非常小，幾乎看不到差異，從而驗(yàn)證了解析解的正確性.

圖1 平板前表面遭受突變熱擾動(dòng)時(shí)后表面溫度響應(yīng)的數(shù)值解與解析解的比較

圖2 平板前表面遭受突變熱擾動(dòng)時(shí)不同熱松弛時(shí)間Ve下的表面溫度響應(yīng) (a)前表面溫度響應(yīng);(b)后表面溫度響應(yīng)

圖2給出了平板前表面遭受突變熱擾動(dòng)時(shí)不同無(wú)量綱熱松弛時(shí)間Ve下，平板兩表面無(wú)量綱溫度隨無(wú)量綱時(shí)間的變化曲線.可以看出當(dāng)Ve=0.3，0.6，0.8時(shí)，其溫度分布曲線與典型的Fourier溫度分布曲線已不一致，尤其是當(dāng)Ve=0.6，0.8時(shí)，其曲線有了凸凹點(diǎn)，熱傳導(dǎo)的非Fourier效應(yīng)變得更加明顯.從圖2(b)可以看出其溫度響應(yīng)曲線均存在一明顯的延遲，由于熱波的傳播速度為1/Ve，所以在平板后表面(無(wú)量綱化后X=1)，其延遲時(shí)間為Fo0=Ve.

4.2 簡(jiǎn)諧變化邊界條件

當(dāng)厚度為L(zhǎng)的平板前表面遭受一幅值為q0，角頻率為α的簡(jiǎn)諧熱擾動(dòng)(q(t)=cos(αt))時(shí)

式中，F(xiàn)o1=a/αL2.

把(41)式代入(34)和(35)式得平板前表面熱擾動(dòng)簡(jiǎn)諧變化時(shí)平板內(nèi)部的溫度場(chǎng)分布

式中

對(duì)于Fourier導(dǎo)熱，當(dāng)Ve→0時(shí)，把(41)式代入(34)和(37)式可得

使用(42)式，可以模擬得到厚度為L(zhǎng)的平板在其前表面受到簡(jiǎn)諧熱擾動(dòng)時(shí)，平板前后兩表面的溫度響應(yīng).

圖3給出了平板前表面在簡(jiǎn)諧熱流作用下，當(dāng)Fo1=0.25，Ve=0.8時(shí)，平板前后兩表面的溫度響應(yīng)曲線，反映出了非Fourier熱傳導(dǎo)的瞬態(tài)溫度特性.從圖3可以看出兩條曲線都存在一系列階躍點(diǎn)，它們表示溫度波波前經(jīng)過(guò)邊界的傳播和反射到達(dá)的時(shí)刻.由于熱波的傳播速度為1/Ve，所以在兩表面X(X=0,1)處，階躍點(diǎn)位置為Fo=Ve·X，Ve(2+X)，Ve(4+X),···.由于擴(kuò)散的作用，前表面溫度波的階躍值和幅值要比后表面的大.從圖3還可以看出溫度波曲線的階躍點(diǎn)隨著時(shí)間的增加而減小并逐漸消失，曲線最終變成光滑的周期曲線.

圖3 平板前表面遭受簡(jiǎn)諧熱擾動(dòng)時(shí)前后兩表面的溫度響應(yīng)(Fo1=0.25，Ve=0.8)

圖4 平板前表面遭受簡(jiǎn)諧熱擾動(dòng)時(shí)不同熱松弛時(shí)間Ve下前后兩表面的溫度響應(yīng)(Fo1=0.25) (a)前表面溫度響應(yīng);(b)后表面溫度響應(yīng)

圖4給出了在前表面遭受簡(jiǎn)諧熱擾動(dòng)和Fo1=0.25的情況下，不同無(wú)量綱熱松弛時(shí)間Ve下，平板前后兩表面處無(wú)量綱溫度隨無(wú)量綱時(shí)間的變化曲線，以及采用(44)式得到的Fourier熱傳導(dǎo)模型下的溫度響應(yīng)曲線.從圖中可以看出，隨著熱松弛時(shí)間Ve的減小，熱波的傳播速度增大，熱量能夠迅速地傳遞，從而導(dǎo)致平板內(nèi)溫度響應(yīng)的幅值也隨之減小，非Fourier效應(yīng)逐漸減弱.Ve越小，非Fourier和Fourier的溫度響應(yīng)曲線越接近.

4.3 任意周期邊界條件

當(dāng)厚度為L(zhǎng)的平板前表面遭受任意周期熱流變化時(shí)，假設(shè)該周期函數(shù)q0·q(t)的周期為2l.在實(shí)際問(wèn)題中，q(t)只能在t≥0上有定義，故我們可以在t＜0的區(qū)間內(nèi)將函數(shù)q(t)進(jìn)行偶延拓，使q(t)=q(-t)，即延拓后的函數(shù)為偶函數(shù).這樣周期函數(shù)q(t)可以展開(kāi)成余弦級(jí)數(shù)的形式

當(dāng)i=0時(shí)，方程(7)的解可以由(39)式得到：

把(50)與(52)式進(jìn)行疊加，可得平板表面溫度任意周期變化時(shí)的非Fourier溫度場(chǎng)：

當(dāng)Ve→0時(shí)，根據(jù)疊加原理與(40)，(44)和(45)式，可得在任意周期邊界條件下的非Fourier導(dǎo)熱溫度場(chǎng)：

4.4 應(yīng)用舉例

設(shè)q(t)是周期為2π的三角波函數(shù)(如圖5所示)，它在[0,2π)上的表達(dá)式為

將函數(shù)q(t)進(jìn)行偶延拓后，展開(kāi)成余弦級(jí)數(shù)

將(57)式分別代入(53)式和(54)式得三角波條件下的非Fourier和Fourier導(dǎo)熱溫度場(chǎng).

圖5 三角波熱擾動(dòng)

圖6給出了在前表面遭受三角波熱擾動(dòng)和a/L2=0.25的情況下，不同無(wú)量綱熱松弛時(shí)間Ve下，平板前后兩表面處無(wú)量綱溫度隨無(wú)量綱時(shí)間的變化曲線.從圖6可以看出，該溫度響應(yīng)曲線除具有前面介紹的規(guī)律外，還具有隨著Ve的增大，特別是Ve=0.50時(shí)，溫度響應(yīng)曲線形狀越接近前表面輸入的三角波熱擾動(dòng)信號(hào)這樣的現(xiàn)象.這是因?yàn)?，Ve越大，熱波的傳播速度越小，其延遲效應(yīng)越明顯，隨著時(shí)間的增加，當(dāng)階躍消失時(shí)，其溫度響應(yīng)曲線形狀將越接近輸入信號(hào).

圖6 平板前表面遭受三角波熱擾動(dòng)時(shí)不同熱松弛時(shí)間Ve下前后兩表面的溫度響應(yīng)(a/L2=0.25) (a)前表面溫度響應(yīng);(b)后表面溫度響應(yīng)

5 結(jié)論

本文首先通過(guò)分離變量法和Duhamel積分原理，對(duì)平板前表面遭受突變熱擾動(dòng)和簡(jiǎn)諧熱擾動(dòng)這兩類邊界條件下的非Fourier熱傳導(dǎo)問(wèn)題進(jìn)行了分析求解，在此基礎(chǔ)上，利用Fourier級(jí)數(shù)展開(kāi)法和疊加原理，得到了雙曲線型熱傳導(dǎo)方程在平板前表面遭受任意周期變化熱流時(shí)這個(gè)最一般情況下的解析解.通過(guò)理論計(jì)算與數(shù)值模擬，展示了非Fourier熱傳導(dǎo)模型所給出的溫度響應(yīng)與Fourier熱傳導(dǎo)模型的如下差別：

1)溫度響應(yīng)曲線存在一系列的階躍點(diǎn)，在兩表面X(X=0,1)處，它們依次出現(xiàn)的時(shí)間為Fo=Ve·X，Ve(2+X)，Ve(4+X)，···，但隨著時(shí)間的延長(zhǎng)，溫度響應(yīng)曲線逐漸變得光滑;

2)同一位置處，溫度響應(yīng)的幅值隨著熱松弛時(shí)間Ve的減小而減小，最終趨近于Fourier熱傳導(dǎo)時(shí)的溫度響應(yīng)幅值;

3)具有有限的傳播速度，其大小與1/Ve成線性關(guān)系;

4)熱松弛時(shí)間Ve越大，溫度響應(yīng)曲線的形狀越和前表面遭受的熱流形狀相似.

[1]Wang Y Z，Song X N 2012Acta Phys.Sin.61 234601(in Chinese)[王穎澤，宋新南2012物理學(xué)報(bào)61 234601]

[2]Guo Z Y，Cao B Y 2008Acta Phys.Sin.57 4273(in Chinese)[過(guò)增元，曹炳陽(yáng)2008物理學(xué)報(bào)57 4273]

[3]Tian X G，Shen Y P 2012Adv.Mech.42 18(in Chinese)[田曉耕，沈亞鵬2012力學(xué)進(jìn)展42 18]

[4]Tung T L，F(xiàn)ong E 2011Int.J.Heat Mass Transfer54 4796

[5]Cattaneo C 1948Atti.Sem.Mat.Fis.Univ.Modena3 83

[6]Vernotte P 1958C.R.Acad.Sci.246 3154

[7]Tao Y J，Huai X L，Li Z G 2006Chin.Phys.Lett.23 2487

[8]Li S R，Zhou F X，Wu H M 2007Engineer.Mech.24 48(in Chinese)[李世榮，周鳳璽，吳紅梅2007工程力學(xué)24 48]

[9]Sarkar D，Haji-Sheikh A 2012Int.Commun.Heat Mass Transfer39 1009

[10]Tang D W，Araki N 1996Int.J.Heat Mass Transfer39 1585

[11]Moosaie A 2007Int.Commun.Heat Mass Transfer34 996

[12]Barletta A，Zanchini E 1997Heat and Mass Transfer32 285

[13]Ate fiG，Talaee M R 2011Arch.Appl.Mech.81 569

[14]Mishra S C，Sahai H 2012Int.J.Heat Mass Transfer55 7015

[15]Jiang F M 2006Heat and Mass Transfer42 1083

[16]Shirmohammadi R，Moosaie A 2009Int.Commun.Heat Mass Transfer36 827