球?qū)ΨQ瞬態(tài)非傅里葉熱傳導(dǎo)分析1)

劉承斌 龔順風(fēng) 王惠明

?(浙江大學(xué)土木工程學(xué)院，杭州310058)

?(浙江大學(xué)工程力學(xué)系，杭州310027)

球殼結(jié)構(gòu)因其受力合理在工程中有廣泛的應(yīng)用，大至儲(chǔ)油罐，小至壓電傳感器等。當(dāng)球殼結(jié)構(gòu)工作于溫差突然有較大變化的環(huán)境中時(shí)，急劇的溫度變化使得結(jié)構(gòu)內(nèi)部產(chǎn)生瞬態(tài)熱應(yīng)力。通常將這種短促時(shí)間內(nèi)出現(xiàn)的熱應(yīng)力數(shù)值上非常大的現(xiàn)象稱為熱沖擊，熱沖擊有可能產(chǎn)生很大的拉壓應(yīng)力，甚至可能達(dá)到脆性材料的破壞應(yīng)力，從而造成材料破壞。以往常規(guī)熱分析大多基于經(jīng)典傅里葉熱傳導(dǎo)定律，這種經(jīng)典傅里葉熱傳導(dǎo)的描述可以滿足絕大部分工程理論分析的精度要求，但存在熱傳播速度無(wú)窮大這一與實(shí)際物理規(guī)律不符合的弊端。隨著研究的不斷深入和應(yīng)用范圍的不斷擴(kuò)大，特別是在極端環(huán)境下服役時(shí)，經(jīng)典理論已不能滿足分析要求[1]。為克服經(jīng)典傅里葉熱傳導(dǎo)定律的局限性，Cattaneo[2]和 Vernotte[3]給出了能夠描述熱以有限速度傳播的雙曲型熱傳導(dǎo)方程，被稱為C-V 方程或非傅里葉熱傳導(dǎo)方程。這種廣義熱傳導(dǎo)問(wèn)題可采用積分變換、差分法、有限元以及解析法等多種方法進(jìn)行求解。積分變換方法將時(shí)間域的偏微分方程變換為頻域的常微分方程，再進(jìn)行求解，然后進(jìn)行積分反變換，求得時(shí)域響應(yīng)。Jiang等[4]利用拉普拉斯變換，針對(duì)某些邊界條件，得到了空心球殼受熱沖擊的解析解。Babaei 等[5]利用拉普拉斯變換和常規(guī)數(shù)值反變換得到功能梯度空心球殼的廣義熱傳導(dǎo)響應(yīng)。Conker 等[6]將拉普拉斯變換和修正的數(shù)值反變換方法相結(jié)合，探討了材料特性按指數(shù)變化的功能梯度厚壁圓柱殼和球殼雙曲型非傅里葉熱傳導(dǎo)問(wèn)題。Wang 等[7]利用積分變換研究了考慮界面裂紋的層狀復(fù)合材料的非傅里葉熱傳導(dǎo)現(xiàn)象。有限差分方法根據(jù)不同的邊界條件和初始條件選擇合適的差分格式，將問(wèn)題離散化后求解，所采用的差分格式?jīng)Q定了求解精度，且必須滿足收斂性和穩(wěn)定性要求[8]。廣義熱彈性有限元方法，可避免數(shù)值積分反變換造成的截?cái)嗾`差，直接在時(shí)間域求解，但推導(dǎo)相對(duì)復(fù)雜[9-11]。分離變量法[12-15]在球殼表面受熱沖擊作用的熱傳導(dǎo)分析中也非常有效。此外，Rahideh等[16]利用微分求積法研究了疊層功能梯度平面域非傅里葉熱傳導(dǎo)問(wèn)題。張浙等[17]對(duì)非傅里葉熱傳導(dǎo)的研究進(jìn)展做了評(píng)述，從實(shí)質(zhì)、模型的建立和求解及應(yīng)用于實(shí)驗(yàn)等多個(gè)方面進(jìn)行了論述，并指出了今后的研究方向。蔣方明等[18]進(jìn)一步介紹了非傅里葉導(dǎo)熱的研究進(jìn)展。近年來(lái)，Moosaie[19]研究了空心球的球?qū)ΨQ瞬態(tài)熱傳導(dǎo)問(wèn)題，為便于求解，將解分為兩部分之和：一部分穩(wěn)態(tài)解對(duì)應(yīng)非齊次邊界條件，另一部分瞬態(tài)解對(duì)應(yīng)齊次邊界條件。Zhao 等[20]通過(guò)分離變量法結(jié)合杜哈梅原理，分析了實(shí)心球外受任意表面熱沖擊條件下的非傅里葉導(dǎo)熱響應(yīng)。許光映等[21]利用積分變換方法探討了非傅里葉導(dǎo)熱邊界條件對(duì)激光輻射生物組織熱傳導(dǎo)的影響。劉承斌等[22]從三維壓電彈性理論出發(fā)，利用拉普拉斯變換結(jié)合狀態(tài)空間法，分析了壓電球殼受表面球?qū)ΨQ熱沖擊作用下的熱彈性問(wèn)題。張彥博等[23]基于雙曲型單相延遲非傅里葉熱傳導(dǎo)方程，采用有限元法研究了含裂紋厚壁圓筒結(jié)構(gòu)在熱沖擊載荷下的熱力學(xué)數(shù)值解。陳豪龍等[24]提出微分轉(zhuǎn)換雙重互易邊界元法求解功能梯度材料非傅里葉瞬態(tài)熱傳導(dǎo)問(wèn)題。結(jié)果表明時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果的影響較小，該方法可以有效求解功能梯度材料非傅里葉瞬態(tài)熱傳導(dǎo)問(wèn)題。郭松林[25]基于雙相延遲非傅里葉熱傳導(dǎo)模型對(duì)工程中常見(jiàn)的圓柱、板、涂層以及夾芯板等結(jié)構(gòu)在熱沖擊載荷作用下的裂紋問(wèn)題進(jìn)行了分析，系統(tǒng)分析了結(jié)構(gòu)熱彈性響應(yīng)和裂紋尖端應(yīng)力強(qiáng)度因子。

綜上所述，盡管廣義熱傳導(dǎo)問(wèn)題已有很多有效的求解方法，但隨著新材料技術(shù)的不斷發(fā)展和服役環(huán)境的更加嚴(yán)苛，需要尋求計(jì)算效率高且準(zhǔn)確的方法。本文針對(duì)球?qū)ΨQ瞬態(tài)非傅里葉熱傳導(dǎo)問(wèn)題，給出了兩種不同的求解方法，并通過(guò)具體算例分析比較了兩種方法的計(jì)算效果。

1 球殼瞬態(tài)熱傳導(dǎo)基本方程

考慮內(nèi)外半徑分別為a和b的空心球殼受球?qū)ΨQ熱沖擊作用 (圖 1)，球?qū)ΨQ情形的 L–S 型廣義熱傳導(dǎo)方程為

式中，k33為徑向熱傳導(dǎo)系數(shù)，τ為熱松弛時(shí)間，T為溫度變化值，ρ為質(zhì)量密度，cv為比熱容。由非傅里葉定理可知

式中,?2=r?/?r,Θr=rqr，qr為徑向熱流密度,符號(hào)上面的·表示對(duì)時(shí)間t求導(dǎo)。

圖1 球坐標(biāo)及球殼示意圖

存在三種不同的溫度邊界條件：給定邊界處的溫度、熱流密度或者兩者的線性組合。為便于推導(dǎo)，本文考慮溫度邊界條件，其他邊界條件采用同樣的步驟可以類推。球殼內(nèi)外表面的溫度邊界條件和初始條件可表示為

式中,φ0和φ1為某特定值，H(t) 為階躍函數(shù)。

2 頻域求解

2.1 狀態(tài)方程的建立

將式(2) 代入式(1)，整理后得到

為消除時(shí)間變量，將偏微分方程組轉(zhuǎn)化為常微分方程組，引入拉普拉斯變換

將式(2) 和式(5) 利用式(6) 變換后，聯(lián)合得到狀態(tài)方程

球殼邊界條件式(3) 經(jīng)過(guò)拉普拉斯變換后變?yōu)?/p>

2.2 狀態(tài)方程的求解

為求解方便，對(duì)變量進(jìn)行無(wú)量綱化

其中c33為材料的彈性參數(shù)，c為彈性波速，η=ρcv/k33為熱黏度系數(shù)。為簡(jiǎn)便考慮，省略后續(xù)公式符號(hào)中的上標(biāo)符號(hào)“*”。將式(9) 代入式(7)，可得

式中V= [T,Θr]T。由于矩陣M中的元素含有與半徑r相關(guān)的量，可采用層合模型將其轉(zhuǎn)化為常系數(shù)矩陣，當(dāng)分層的模型足夠細(xì)時(shí)，這種近似方法有足夠的精度。將球殼沿厚度方向等分為p層，矩陣中的坐標(biāo)變量r取為子層中間厚度處坐標(biāo)值，即r=a+(2i ?1)(b ?a)/2p，由矩陣?yán)碚摰檬?(10) 的解為

相鄰子層間的連續(xù)條件要求狀態(tài)變量連續(xù)，由式(11)出發(fā)遞推可得

為二階方陣。引入球殼內(nèi)外溫度邊界條件式(8)，代入式 (12)，得

從而可求得球殼內(nèi)表面的熱流密度狀態(tài)量

式中Tij為矩陣T的元素。利用式(11) 可以求得沿徑向任意位置的狀態(tài)量溫度值。

對(duì)于大部分邊界條件, 無(wú)法直接將頻域解解析為時(shí)域解，必須采用近似方法進(jìn)行拉普拉斯數(shù)值反變換[26-29]。本文采用Durbin[26]的傅里葉級(jí)數(shù)法

其中U為反變換的周期，合理的參數(shù)取值為 5αU10，50NUSM5000。

3 時(shí)域求解

3.1 分解

將式(9) 代入式(5) 化簡(jiǎn)后為

由于邊界條件式(3)非齊次，不能直接利用分離變量法來(lái)求解方程 (16)，因而需要進(jìn)行邊界條件的齊次化，將溫度場(chǎng)解分解成兩部分：對(duì)應(yīng)非齊次邊界條件的解TS(r) 和對(duì)應(yīng)齊次邊界條件的解TD(r,t)，即

3.2 穩(wěn)態(tài)解 TS(r)

非齊次邊界條件下的熱傳導(dǎo)邊值問(wèn)題可以表示為

式 (18) 的解為

其中A0和B0為待定常數(shù)，將式 (20) 代入邊界條件式 (19)，得到

3.3 瞬態(tài)解 TD(r)

滿足齊次邊界條件及初始條件的球?qū)ΨQ熱傳導(dǎo)問(wèn)題為

假設(shè)TD(r,t) =RT(r)L(t)，將其代入式 (22)，從而得到

考慮式(25) 兩端只能等于常數(shù)?ω2，分解成兩個(gè)2階常微分方程

相應(yīng)的邊界條件為

式(26) 為貝塞爾方程，其解析解為

其中C0和D0為待定系數(shù)。將式 (29) 代入邊界條件式(28)，可以得到

為保證式(30) 存在非零解，系數(shù)矩陣行列式必須為零，即

求解上述超越方程，可以得到無(wú)窮個(gè)特征根ωk(k=1,2,···)，對(duì)于任意ωk，式 (29) 的解為

二階常微分方程(27) 的解為

由初始條件式(24) 第1 式得

對(duì)應(yīng)不同本征值的本征函數(shù)RTk(r)和RTj(r)正交性證明推導(dǎo)過(guò)程及模值詳見(jiàn)文獻(xiàn)[30]，這里直接給出模值

利用級(jí)數(shù)正交展開技術(shù)，由式(35) 可得

由初始條件式(24) 第2 式得

3.4 全解

將穩(wěn)態(tài)解和瞬態(tài)解兩部分疊加，得到總的溫度場(chǎng)

4 算例

為驗(yàn)證本文解的正確性和有效性，將兩種方法與文獻(xiàn)[31] 計(jì)算結(jié)果進(jìn)行對(duì)比。球殼無(wú)量綱參數(shù)為：內(nèi)半徑a= 1，外半徑b= 2，波速c1= 1，c2=0.535，ε=0.02。溫度邊界條件為

求解時(shí)，頻域法采用層合近似法將球殼沿厚度方向細(xì)分60 等分，時(shí)域法采用前70 階特征根，和文獻(xiàn)[31] 得到的球殼內(nèi)表面溫度隨時(shí)間變化的結(jié)果吻合得非常好，如圖 2 所示，表明本文理論推導(dǎo)和數(shù)值計(jì)算的正確性。

圖2 球殼內(nèi)表面溫度隨時(shí)間變化

下面以硒化鎘球殼受內(nèi)外溫差變化為例，材料參數(shù)為：c33= 83.6 GPa,ρ= 5.68 t/m3,k33=12.9 W/(m·K),cv= 0.46×103J·kg/K。無(wú)量綱熱弛豫時(shí)間τ= 0.1，球殼內(nèi)外徑之比a/b取為0.5，內(nèi)/外表面分別同時(shí)突加溫升 200?C/400?C。圖3 為采用頻域法所求得對(duì)應(yīng)不同時(shí)刻內(nèi)部溫度變化沿著球殼徑向分布情況，計(jì)算時(shí)，將球殼沿厚度方向80 等分。觀察后可以發(fā)現(xiàn)，內(nèi)外表面熱沖擊所產(chǎn)生的熱波，以無(wú)量綱速度向球殼對(duì)向傳播。在溫度波前尚未到達(dá)的區(qū)域內(nèi)，溫度保持初始溫度不變，當(dāng)熱波波前到達(dá)某位置，該處形成一個(gè)明顯的溫度階躍現(xiàn)象，當(dāng)雙向的熱波交匯后，球體內(nèi)部的溫度會(huì)超過(guò)邊界處設(shè)定的溫度。產(chǎn)生的原因在于熱傳導(dǎo)時(shí)間與熱馳豫時(shí)間相當(dāng)時(shí)，材料的局部熱平衡已經(jīng)不能維持[1]。

圖3 對(duì)應(yīng)不同時(shí)間t 溫度沿球殼厚度方向變化(頻域法)

圖4 為時(shí)域法求得的結(jié)果，計(jì)算時(shí)采用了前70階頻率。對(duì)比圖3 和圖4 可以看出，兩種不同方法所求得的溫度分布具有大致相似的結(jié)果，表明兩種計(jì)算方法的有效性。頻域法所求得結(jié)果不存在明顯的波動(dòng)和階躍現(xiàn)象，原因?yàn)閿?shù)值反變換近似方法所固有的截?cái)嗾`差和數(shù)值震蕩所導(dǎo)致；時(shí)域方法求解的溫度結(jié)果存在更明顯的階躍，可以避免拉普拉斯反變換所帶來(lái)的數(shù)值誤差，顯示其可清晰捕捉熱波傳播的前沿波陣面，表明采用該方法求解具有更高的精度。

圖4 對(duì)應(yīng)不同時(shí)間t 溫度沿球殼厚度方向變化(時(shí)域法)

5 結(jié)論

本文針對(duì)球?qū)ΨQ非傅里葉瞬態(tài)熱傳導(dǎo)問(wèn)題，給出了兩種不同的求解方法，并采用具體算例進(jìn)行比較，經(jīng)分析，可以得到如下結(jié)論：

(1) 頻域法求解過(guò)程清晰統(tǒng)一，計(jì)算簡(jiǎn)便，容易推廣至材料參數(shù)沿徑向任意變化的功能梯度球殼熱傳導(dǎo)問(wèn)題的求解。

(2) 時(shí)域法求解時(shí)結(jié)合了分離變量法和疊加法，避免了拉普拉斯數(shù)值反變換帶來(lái)的誤差，可更準(zhǔn)確分析球殼的瞬態(tài)非傅里葉熱傳導(dǎo)問(wèn)題。

(3)頻域法求解時(shí)僅需考慮層合近似、矩陣計(jì)算及拉普拉斯數(shù)值反變換，計(jì)算效率較高。時(shí)域法需要采用二分法求解特征根，然后進(jìn)行累加，計(jì)算效率相對(duì)較低。