趙偉濤 吳九匯?

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

2)(西安交通大學(xué)，機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國家重點(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)熱問題分析的各個(gè)領(lǐng)域.對(duì)于熱作用時(shí)間較長的穩(wěn)態(tài)傳熱過程以及熱傳播速度較快的非穩(wěn)態(tài)常規(guī)傳熱過程，采用Fourier定律來描述熱流密度與溫度梯度之間的關(guān)系是可以滿足精度要求的.但是，F(xiàn)ourier定律不涉及傳熱時(shí)間項(xiàng)，隱含了熱擾動(dòng)傳播速度為無限大的假設(shè).隨著科學(xué)技術(shù)的進(jìn)步，超短激光脈沖的出現(xià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)熱問題.李世榮等[8]采用單相延遲模型，研究了薄板在受周期熱流邊界條件下板內(nèi)的溫度響應(yīng).Sarkar和Haji-Sheikh[9]采用Laplace變換技術(shù)研究了由介電材料制成的有限平板的雙曲型熱傳導(dǎo)問題.Tang和Araki[10]求解了表面周期性加熱條件下有限介質(zhì)的非Fourier熱傳導(dǎo)問題.Moosaie[11]在Tang等的工作基礎(chǔ)上，運(yùn)用Fourier積分表達(dá)式，求解了有限介質(zhì)在周期表面熱擾動(dòng)條件下的非Fourier熱傳導(dǎo)問題.Barletta和Zanchini[12]分析了無限圓柱體存在內(nèi)熱源以及與外界流體有熱對(duì)流的情況下的雙曲型熱傳導(dǎo)問題.Ate fi和Talaee[13]使用分離變量法求解了邊界條件不隨時(shí)間變化的無限長圓筒的非Fourier溫度場.Mishra和Sahai[14]運(yùn)用格子Boltzmann法研究了一維圓柱和球體的非Fourier熱傳導(dǎo)問題.Jiang[15]運(yùn)用Laplace變換法研究了空心球體在內(nèi)外兩個(gè)表面溫度突然變化時(shí)的雙曲型熱傳導(dǎo)問題.Shirmohammadi和Moosaie[16]采用分離變量法得到空心球體在周期表面熱流條件下的解析解.

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

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

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

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

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

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

3 分析求解

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

設(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)部的溫度場.

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

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

4.1 突變邊界條件

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

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

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

根據(jù)(39)式，可以計(jì)算得到厚度為L的平板前面遭受一突變熱擾動(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í)不同無量綱熱松弛時(shí)間Ve下，平板兩表面無量綱溫度隨無量綱時(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，所以在平板后表面(無量綱化后X=1)，其延遲時(shí)間為Fo0=Ve.

4.2 簡諧變化邊界條件

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

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

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

式中

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

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

圖3給出了平板前表面在簡諧熱流作用下，當(dāng)Fo1=0.25，Ve=0.8時(shí)，平板前后兩表面的溫度響應(yīng)曲線，反映出了非Fourier熱傳導(dǎo)的瞬態(tài)溫度特性.從圖3可以看出兩條曲線都存在一系列階躍點(diǎn)，它們表示溫度波波前經(jīng)過邊界的傳播和反射到達(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 平板前表面遭受簡諧熱擾動(dòng)時(shí)前后兩表面的溫度響應(yīng)(Fo1=0.25，Ve=0.8)

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

圖4給出了在前表面遭受簡諧熱擾動(dòng)和Fo1=0.25的情況下，不同無量綱熱松弛時(shí)間Ve下，平板前后兩表面處無量綱溫度隨無量綱時(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的平板前表面遭受任意周期熱流變化時(shí)，假設(shè)該周期函數(shù)q0·q(t)的周期為2l.在實(shí)際問題中，q(t)只能在t≥0上有定義，故我們可以在t＜0的區(qū)間內(nèi)將函數(shù)q(t)進(jìn)行偶延拓，使q(t)=q(-t)，即延拓后的函數(shù)為偶函數(shù).這樣周期函數(shù)q(t)可以展開成余弦級(jí)數(shù)的形式

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

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

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

4.4 應(yīng)用舉例

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

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

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

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

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

本文首先通過分離變量法和Duhamel積分原理，對(duì)平板前表面遭受突變熱擾動(dòng)和簡諧熱擾動(dòng)這兩類邊界條件下的非Fourier熱傳導(dǎo)問題進(jìn)行了分析求解，在此基礎(chǔ)上，利用Fourier級(jí)數(shù)展開法和疊加原理，得到了雙曲線型熱傳導(dǎo)方程在平板前表面遭受任意周期變化熱流時(shí)這個(gè)最一般情況下的解析解.通過理論計(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í)間的延長，溫度響應(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)[過增元，曹炳陽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