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

* (西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院,復(fù)雜系統(tǒng)動力學(xué)與控制工信部重點實驗室,西安 710129)

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

引言

熱傳導(dǎo)理論是熱彈耦合理論的重要基礎(chǔ),經(jīng)典傅里葉熱傳導(dǎo)定律為

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

為此,國內(nèi)外學(xué)者提出了多種廣義熱傳導(dǎo)模型,并建立了相應(yīng)的廣義熱彈性理論.Cattaneo[7]和Vernotte[8]提出了CV (Cattaneo-Vernotte)熱傳導(dǎo)模型

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

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

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

通過量綱分析,引入關(guān)系式k*=k/τ,則GN 熱傳導(dǎo)方程可進一步改寫為

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

同樣地,考慮到關(guān)系式k*=k/τ,則有

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

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

綜上,廣義熱傳導(dǎo)模型的形式多樣,對應(yīng)的熱彈耦合理論也很豐富.因此,存在的問題有: (1) 如何從連續(xù)介質(zhì)力學(xué)視角,建立起相關(guān)理論模型的熱力學(xué)基礎(chǔ);(2) 各模型之間存在怎樣的關(guān)聯(lián),如何揭示之.

基于拓展熱力學(xué)原理,本文建立了廣義熱彈耦合理論的新模型,是現(xiàn)有相關(guān)模型的統(tǒng)一形式,可通過參數(shù)賦值實現(xiàn)模型的退化.此外,提出了熱學(xué)“彈性”單元和“黏性”單元,通過串并聯(lián)模式實現(xiàn)了上述模型的重構(gòu),從而揭示了已有模型間的關(guān)系.結(jié)合數(shù)值分析結(jié)果,研究了各模型中松弛時間等對瞬態(tài)響應(yīng)結(jié)果的影響規(guī)律.

1 拓展熱力學(xué)和廣義熱彈耦合理論

根據(jù)連續(xù)介質(zhì)力學(xué),熱力學(xué)第一定律為

其中,ρ 為密度,U為比內(nèi)能,r為內(nèi)熱源,σ 為應(yīng)力張量,ε 為應(yīng)變張量.熱力學(xué)第二定律為

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

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

根據(jù)拓展熱力學(xué),現(xiàn)引入高階熱流[34]

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

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

其中

代入式(16),有

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

可得

經(jīng)典熱力學(xué)的熵流為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)可更新為

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

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

2 熱學(xué)“彈性”單元和“黏性”單元及廣義熱傳導(dǎo)方程的重構(gòu)

上節(jié)中,基于拓展熱力學(xué),通過考慮熱彈耦合得到了應(yīng)力、熵本構(gòu)關(guān)系和廣義熱傳導(dǎo)方程.若將其與運動學(xué)方程、能量守恒方程和應(yīng)變-位移關(guān)系結(jié)合,則可建立廣義熱彈耦合模型.同時發(fā)現(xiàn),本文得到的廣義熱傳導(dǎo)模型可退化為CV,DPL,GN,MGT 等模型.

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

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

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

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

2.1 CV 模型

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

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

對式(39)求導(dǎo)

其中

因此

整理后可得CV 熱傳導(dǎo)模型

2.2 GN 模型

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

可得到

整理后可得GN 熱傳導(dǎo)模型

2.3 DPL 模型

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

對第1 部分有

對第2 部分有

聯(lián)立式(47)～式(49),可得

整理可得DPL 型熱傳導(dǎo)方程

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

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

2.4 MGT 模型

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

對第1 部分有

對第2 部分有

聯(lián)立式(53)～式(55)可得

其中

聯(lián)立式(56)和式(57)有

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

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

2.5 新模型——熱場Burgers 模型

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

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

聯(lián)立式(60)～式(62)并求導(dǎo)可得

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

表1 廣義熱傳導(dǎo)模型的變換和統(tǒng)一化Table 1 Transformation and unification of generalized heat conduction models

3 算例計算

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

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

考慮了熱力耦合的應(yīng)力本構(gòu)方程為

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

能量守恒方程為

其中Q為熱源.熵本構(gòu)方程為

為簡單起見,表1 中的廣義熱傳導(dǎo)方程統(tǒng)一寫為

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

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

定義無量綱變量如下

其中

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

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

其中

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

3.1 熱源函數(shù)=0

式(79)退化為

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

邊界條件經(jīng)拉普拉斯變換后為

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

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

3.2 移動熱源=L(Q)=

式(79)具體表示為

拉普拉斯域中的溫度、應(yīng)力解同樣可表示為

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

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

邊界條件經(jīng)拉普拉斯變換為

聯(lián)立式(94)和式(95),可得

其中

4 結(jié)果與分析

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

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

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

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

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

圖3 移動熱源下各模型的位移、應(yīng)力、溫度響應(yīng)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)為該作用環(huán)境下位移隨時間演化的規(guī)律.兩模型均在經(jīng)歷過移動熱源的加熱后產(chǎn)生較大位移,其中v=2熱源作用下的位移響應(yīng)更大.

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

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

5 總結(jié)

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