廣義熱彈模型的熱力學基礎與瞬態響應1)

吳 華 鄒紹華 徐成輝 尉亞軍,?,2) 鄧子辰,3)

* (西北工業大學力學與土木建筑學院,復雜系統動力學與控制工信部重點實驗室,西安 710129)

? (西安交通大學機械結構強度與振動國家重點實驗室,西安 710049)

引言

熱傳導理論是熱彈耦合理論的重要基礎,經典傅里葉熱傳導定律為

其中,k為熱傳導系數,q為熱流,T為溫度.傅里葉定律及其熱彈耦合理論被廣泛用于描述熱力共同作用下材料/結構的熱傳導和變形行為.同時,傅里葉定律表明[1-3]: 熱以無限大速度傳導,即一點的溫度變化能瞬間傳導到無窮遠處.解決這一理論“悖論”成為不斷推動熱傳導方程深入發展的動力.同時,飛速發展的激光、集成電路等技術領域,具有作用載體和作用環境極小、作用時間極短的特點,使得器件在極小空間、極短時間內產生大量熱量.建立在宏觀、近平衡態的經典熱傳導和熱力耦合關系不再適用于此類以微納尺度和超快熱沖擊為代表的極端非平衡態問題[4-6].

為此,國內外學者提出了多種廣義熱傳導模型,并建立了相應的廣義熱彈性理論.Cattaneo[7]和Vernotte[8]提出了CV (Cattaneo-Vernotte)熱傳導模型

其中,τq,τT分別為熱流相、溫度梯度相的松弛時間.若左右兩端分別按一階泰勒級數展開,可得到

Green 等[10-11]提出了GN (Green-Naghdi)熱傳導模型

其中,k*為導熱比,χ為熱位移:χ˙=T.當k*=0 時,該式退化為傅里葉定律.對式(5)求導,可得

通過量綱分析,引入關系式k*=k/τ,則GN 熱傳導方程可進一步改寫為

在GN 模型(5)中也考慮熱流率和松弛時間,可以得到MGT (Moore-Gibson-Thompson) 熱傳導方程[12-13]

同樣地,考慮到關系式k*=k/τ,則有

其他的廣義熱傳導模型包括: Cao 等[14]基于相對論理論,并計及聲子質量,建立的熱質熱傳導方程;Kuang[15]通過引入慣性熵的概念,建立的慣性熵理論等.

廣義熱傳導的研究直接推動了廣義熱彈耦合理論的發展.Lord 等[16]在CV 熱傳導模型的基礎上,建立了LS 廣義熱彈性模型;Bazarra 等[17]在該模型的基礎上探討了材料孔隙率對瞬態響應的影響;Green 等[11]建立了無耗散的GN 熱彈性理論;El-Karamany 等[18]考慮了滯后效應,建立了GN 模型Ⅱ型和Ⅲ型的單相滯后、雙相滯后模型;Alizadeh Hamidi 等[19]探討了GN 模型在Euler-Bernoulli 梁的運用;Tzou[20-21]由理論推導建立了DPL 熱彈耦合理論,針對瞬態響應和滯后效應開展了大量實驗研究;Singh[22]模擬了平面波在橫觀各向同性雙相滯后廣義熱彈性固體中的傳播,以圖形方式顯示了雙相滯后熱彈性和LS 廣義熱彈性情況下的傳播角度;Quintanilla[23]發展了MGT 熱彈性理論,并考慮了雙溫下的擴展,證明了解的多項式衰減;Youssef[24]基于熱質熱傳導方程,構建了雙溫度廣義熱彈性理論;文獻[25-27]建立并發展了GL 熱彈耦合理論,討論了應變率的影響,證明了解的唯一性定理.此外,考慮空間尺度效應和時間記憶效應的熱彈耦合理論也已取得較大進展,如分數階微積分[28-31]和非局部效應[32-33]的引入等.

綜上,廣義熱傳導模型的形式多樣,對應的熱彈耦合理論也很豐富.因此,存在的問題有: (1) 如何從連續介質力學視角,建立起相關理論模型的熱力學基礎;(2) 各模型之間存在怎樣的關聯,如何揭示之.

基于拓展熱力學原理,本文建立了廣義熱彈耦合理論的新模型,是現有相關模型的統一形式,可通過參數賦值實現模型的退化.此外,提出了熱學“彈性”單元和“黏性”單元,通過串并聯模式實現了上述模型的重構,從而揭示了已有模型間的關系.結合數值分析結果,研究了各模型中松弛時間等對瞬態響應結果的影響規律.

1 拓展熱力學和廣義熱彈耦合理論

根據連續介質力學,熱力學第一定律為

其中,ρ 為密度,U為比內能,r為內熱源,σ 為應力張量,ε 為應變張量.熱力學第二定律為

其中,η為比熵,Js為熵流矢量.引入比亥姆霍茲自由能ψ

由式(10)～式(12),有

根據拓展熱力學,現引入高階熱流[34]

其中一階熱流Q(1)是經典熱流q=Q(1).經典熱力學中,熵ρη 是比內能和應變的函數.對于拓展熱力學,假設其也是各階熱流、溫度梯度的函數,則拓展的吉布斯方程為

由式(12)進一步可知,比亥姆霍茲自由能 ψ 是溫度、應變、各階熱流和溫度梯度的函數,故有

其中

代入式(16),有

由式(13)和式(17)可得

可得

經典熱力學的熵流為Js=T-1q,通過考慮高階熱流項,將其拓展為

上式代入式(21),并由正定性得

其中,A1和A2滿足A1+A2=1 ,χi(i=0,1,2,3)≥0 .考慮式(14),并由式(26)可得

將上式代入式(24),有

考慮到A1+A2=1,則式(23)可改寫為

代入式(28)

若取

以及

式(30)可更新為

至此,得到了廣義熱傳導方程.雖然,本文中自由能ρψ與應變、溫度以及各階熱流相關,但由式(19)和式(20)可看出,應力和熵本構關系具有經典熱彈耦合理論的形式.廣義熱傳導方程(34)同時考慮了時空尺度效應,若忽略空間尺度效應,可得

顯然,上式可退化為CV 模型(2)、DPL 模型(4)、GN 模型(7)和MGT 模型(9).

2 熱學“彈性”單元和“黏性”單元及廣義熱傳導方程的重構

上節中,基于拓展熱力學,通過考慮熱彈耦合得到了應力、熵本構關系和廣義熱傳導方程.若將其與運動學方程、能量守恒方程和應變-位移關系結合,則可建立廣義熱彈耦合模型.同時發現,本文得到的廣義熱傳導模型可退化為CV,DPL,GN,MGT 等模型.

為了進一步闡釋各模型之間的內在關聯,類比于黏彈性力學,本節提出熱學“彈性”單元和熱學“黏性”單元,并通過串并聯方式,實現各廣義熱傳導模型的重構.為此,在一維問題中,考慮GN 模型(5)中熱位移χ的概念(=T),熱位移梯度記為 ω=?χ,則傅里葉定律(1)可寫作

GN 模型(5)中若忽略k?T項,可得

類比于力學,稱式(37)為熱學“彈性”單元,式(36)為熱學“黏性”單元.其中q為熱流,k*和k分別為熱學“彈性”單元和熱學“黏性”單元的傳熱系數,滿足關系式k*=k/τ .基于以上單元,現通過串并聯模型,揭示廣義熱傳導模型間的關聯.基本模型如圖1 所示.

圖1 熱學黏彈性單元組合模型Fig.1 Combination model of thermovisco and thermoelastic elements

2.1 CV 模型

如圖1(a)示,熱彈性單元和熱黏性單元串聯組合,其結構與黏彈性理論中的Maxwell 模型一致.模型兩端的熱流相等,即

而模型的總熱位移χ為各單元熱位移之和,則

對式(39)求導

其中

因此

整理后可得CV 熱傳導模型

2.2 GN 模型

如圖1(b)示,熱彈性單元和熱黏性單元并聯組合,其結構與黏彈性理論中的Kelvin 模型一致.模型的熱位移與各單元的熱位移相等,而總熱流為各單元熱流之和

可得到

整理后可得GN 熱傳導模型

2.3 DPL 模型

如圖1(c)所示,兩端總熱流和各部分的熱流相等,總熱位移為各單元熱位移之和

對第1 部分有

對第2 部分有

聯立式(47)～式(49),可得

整理可得DPL 型熱傳導方程

進一步,若設定第1 和2 部分的單元參數一致,即k=k1=k2,τ2=τ,方程改寫為

值得指出的是,該方法推導出的DPL 模型是特殊形式,和一般形式(具有兩個松弛時間,熱流相松弛時間 τq和溫度梯度相松弛時間 τT)存在差異,這是因為推導過程中假定了兩個黏性單元參數保持一致.

2.4 MGT 模型

如圖1(d)示,兩端總熱流和各部分的熱流相等,總熱位移為兩部分之和

對第1 部分有

對第2 部分有

聯立式(53)～式(55)可得

其中

聯立式(56)和式(57)有

同樣地,若k=k1=k2,τ2=τ,,則式(58) 整理為

雖然式(59)和MGT 方程的一般形式仍有差別,但兩式都是由等分別作為函數構成的等式,僅僅是系數的差異.而這一問題可以通過對常系數的賦值加以解決.事實上,區分熱傳導模型也可以通過對比所含的函數這一方法來完成.

2.5 新模型——熱場Burgers 模型

基于以上的組合思想,自然可以形成一種新的模型.如圖1 (f)示,該模型在結構上是由Maxwell 模型和Kelvin 模型串聯組合形成,與黏彈性理論中的Burgers 模型的結構一致.在前文推導的基礎上,可以直接得出第1 和第2 部分的方程分別為

而對于該模型,兩端總熱流和各部分兩端熱流仍然相等,總熱位移為各部分熱位移之和

聯立式(60)～式(62)并求導可得

至此,通過熱學“彈性”和“黏性”單元的串并聯組合,重構了CV,GN,DPL,MGT 熱傳導模型,并通過Burgers 模型得到了第1 節基于拓展熱力學的廣義熱傳導模型.特別指出的是,串并聯模型給出了廣義熱傳導方程中各松弛時間的關系.綜合2.1～2.5 節的推導,廣義熱傳導各模型的控制方程如表1 第2 列所示.

表1 廣義熱傳導模型的變換和統一化Table 1 Transformation and unification of generalized heat conduction models

3 算例計算

考慮一維半無限的材料受邊界熱沖擊和移動熱源的作用,假設無窮遠處不受擾動,且左側邊界無應力,因此,邊界條件可表達為

其中,θ=T-T0為溫度變化值,T0為起始溫度,H(t)為Heaviside 單位階躍函數,u為位移.一維問題下,運動平衡方程為

考慮了熱力耦合的應力本構方程為

其中,λ 和 μ 是拉梅常數,αθ是熱膨脹系數.幾何方程為

能量守恒方程為

其中Q為熱源.熵本構方程為

為簡單起見,表1 中的廣義熱傳導方程統一寫為

其中,M和N是微分算子,如表1 第2 列所示: 對于CV 模型,有M=1+τ?t,N=1 ;對本文新建模型,有M=1+3τ?t+τ2?tt,N=1+τ?t.熱傳導方程經拉普拉斯變換后的變換方程如表1 第3 列所示.

聯立式(65)～式(70),可得位移和溫度的控制方程

定義無量綱變量如下

其中

本文運用拉普拉斯變換方法求解這一問題,為此對式(73)～式(75)做拉普拉斯變換,可得

其中,變換后的M和N(即和) 如表1 第4、5 列所示.聯立式(77)和式(78)可得

其中

至此,可以通過設定不同的熱源函數求解不同條件下的式(79).

3.1 熱源函數=0

式(79)退化為

該方程的特征方程有4 個根,考慮到一維半無限問題下,無窮遠處不受擾動,故舍去兩個正根.位移、溫度、應力解的形式為

邊界條件經拉普拉斯變換后為

聯立式(82)～式(84),可得方程組

聯立式(83)和式(85),可解得溫度表達式為

3.2 移動熱源=L(Q)=

式(79)具體表示為

拉普拉斯域中的溫度、應力解同樣可表示為

將式(88)和式(89)代入式(76),可得

聯立式(76)、式(88)、式(91)和式(92),可解得

邊界條件經拉普拉斯變換為

聯立式(94)和式(95),可得

其中

4 結果與分析

在上一節中對問題進行了描述并給出了數值方法.本節將對所得數據進行討論.(1)在無熱源的情況下,將對四個基本模型和新模型的溫度、位移、應變響應進行對比;(2) 在移動熱源下的各模型溫度、位移、應變響應對比;(3)對比不同熱源移速下x=0.12處DPL 模型和新模型溫度、位移、應變的時間演化特征.材料常數如表2 示.

表2 材料常數表(銅)Table 2 Material constants (copper)

通過對M和N的不同賦值可以實現不同模型的轉化,由此可以計算得到無熱源狀態下各模型t=0.06時刻的空間響應如圖2 所示.圖2(a)為無熱源狀態下各模型位移響應的空間分布特征.新模型表現出相比于其他模型更為激進的變化,而GN 模型則更趨保守.具體為: 新模型會發生更大的峰值位移,同時位移在空間上遞減更快,至熱波邊界x=0.3處為0,變形區域更小.相反地,GN 模型產生最小的峰值位移響應,位移遞減的曲線也最為平緩,變形區域最大.圖2(b)為無熱源狀態下各模型應力響應的空間分布特征.各模型均在約產生突變,CV,MGT和新模型在約x=0.3,即熱波邊界處產生突變,x=0.3處以外的區域熱波尚未到達,不產生應力響應.而GN 和DPL 模型并沒有清晰的波前突變,反映出光滑連續的分布特點.新模型相比于CV 和MGT 模型在波前產生的應力更小.圖2(c)為無熱源狀態下各模型溫度的空間分布特征.與應力特征類似,CV,MGT 和新模型在約x=0.3,即波前處的溫度并不連續,熱邊界條件對x=0.3 以外的區域不產生溫度影響.新模型在波前的溫度響應更小,而GN 和DPL 模型的溫度分布光滑連續,逐步歸于無響應.綜合以上信息,新模型在熱波到達時,位移與應力響應峰值更大,但響應區域更小.新模型和CV,MGT 模型均有明顯的波前突變,在熱波未達區不產生響應,但新模型的波前突變更小.GN 和DPL 模型則呈現光滑連續的變化,響應區域更大.

圖2 無熱源下各模型位移、應力、溫度響應Fig.2 Displacement,stress and temperature responses of each model without heat source

在移動熱源移速v=2 的作用下,各模型位移、應力、溫度的空間分布如圖3 所示.圖3(a)為移動熱源下各模型位移響應的空間分布,相比于無熱源狀態,各模型的位移均有較大增加,尤其在x=0.12處即移動熱源作用處,位移達到最大值.新模型仍然具有各模型中最大的位移,同時位移在空間上的遞減也更快,變形區域最小.圖3(b)為移動熱源下各模型應力響應的空間分布,熱源作用下的應力峰值相比于無熱源狀態明顯增大.無熱源狀態下的應力峰值約為0.05,移動熱源下的峰值應力已躍升至約0.125,且應力峰值出現在移動熱源作用處.新模型表現出最大的峰值應力,而GN 模型應對熱源作用產生的應力變化則最小.圖3(c)為移動熱源下各模型溫度的空間分布.熱源作用下必然會引起各模型溫度峰值的上升,其中尤為明顯的是CV,MGT 和新模型.其溫度曲線也在,即波前處發生突變,回復為常溫狀態.GN 和DPL 模型的溫度曲線是光滑連續的,溫度沿空間逐步降至常溫,故溫度變化的區域更廣.由于移動熱源的作用,導致各模型無論在波前、波后均產生更大的位移、應力和溫度響應,同時在移動熱源作用點x=0.12 處相比無熱源情況明顯升高.同無熱源狀態的響應特點相似的是: CV,MGT 和新模型在波前的響應也存在突變,新模型的突變更小.與無熱源狀態的響應特點不同的是: CV,MGT 和新模型在熱源作用點x=0.12 處也產生了突變,GN 和DPL 模型無論在熱源作用點和波前處,均呈現光滑連續的變化特點.

圖3 移動熱源下各模型的位移、應力、溫度響應Fig.3 Displacement,stress and temperature responses of each model under moving heat source

DPL 模型和新模型在x=0.12處受v=2和v=4的移動熱源作用,移動熱源到達該點的時間分別為t=0.06和t=0.03,其時間演化特征如圖4 所示.圖4 (a)為該作用環境下位移隨時間演化的規律.兩模型均在經歷過移動熱源的加熱后產生較大位移,其中v=2熱源作用下的位移響應更大.

圖4 移動熱源下DPL 模型和新模型位移、應力、溫度的時間演化特征Fig.4 Evolution laws of displacement,stress and temperature of DPL model and new model under moving heat source

新模型的位移響應較DPL 模型峰值更高.圖4(b)為DPL 和新模型的應力時間演化曲線.由于移動熱源的加熱作用,DPL 和新模型均在熱源到達時取得極值.新模型具有比DPL 模型更大的峰值應力.在移動熱源離開該點后,新模型的應力響應迅速回復至與DPL 模型相近的水平.圖4(c)為該作用環境下溫度的時間演化曲線.與圖4(b)類似,DPL 和新模型經歷熱源作用后,溫度均有明顯上升.在溫度上升階段,新模型表現出比DPL 模型更大更快的變化,擁有更高的峰值溫度.新模型的溫度回復也較快,與應力時間演化曲線呈現的特點一致,新模型在受熱源作用前后產生了突變.離開移動熱源作用后,新模型的溫度響應迅速回復至與DPL 模型相近的數值,之后兩模型以較為近似的曲線緩慢下降.熱源移速的變化,使得各項數據的達峰時間隨移速的增大而提前.在時間演化特征曲線上,可以明顯看到新模型在熱源作用時產生的響應突變,而DPL 模型則呈現為光滑連續的變化,這一點和空間響應呈現的特征是一致的.

5 總結

超快熱沖擊等極端熱作用環境下,以傅里葉定律為代表的經典熱傳導理論不再適用,為發展非平衡態熱彈耦合新理論提出迫切需求.學者們通過拓展經典熱傳導方程,已提出多種廣義熱傳導模型,并通過熱力耦合建立了廣義熱彈性理論.為了闡清已有各模型間的內在關聯,并進一步發展極端熱彈耦合理論,本文基于拓展熱力學原理,建立了考慮高階項的廣義熱傳導理論,其可退化為已有相關理論.同時,類比于黏彈性力學理論,本文抽象出熱學“彈性”單元和“黏性”單元模型,基于串并聯方式實現了已有各廣義熱傳導方程的重構,并通過Burgers 模型得到了本文新建立的熱傳導方程,又利用拓展熱力學原理,在原理層面上推導支撐了新模型的理論基礎.廣義熱傳導方程的重構過程也揭示了部分模型中,熱流率相與溫度梯度率相的松弛時間的比例關系.如對DPL 模型而言,熱流率相與溫度梯度率相的松弛時間的比例為2:1.數值結果表明: 無移動熱源作用時,新模型得到的位移、應力、溫度響應峰值更大,影響范圍較小;在移動熱源作用下,新模型預測的響應峰值明顯升高,且高于其他模型.對比DPL 模型,新模型在擁有DPL 模型雙相滯后特點的同時,解決了熱波速度無限大的悖論,在響應曲線上反映出清晰的波前突變.本研究將推動遠離平衡態極端熱彈耦合理論的發展,并為相關問題的求解提供了分析方法.