陳少華 陳學(xué)軍
(北京科技大學(xué)應(yīng)用力學(xué)系,北京 100083)
涂層技術(shù)廣泛應(yīng)用于工程構(gòu)件以提高材料的表面性能,例如在航空領(lǐng)域,涂層技術(shù)可以有效延長發(fā)動機(jī)零部件的使用壽命[1].隨著服役工況要求的不斷提高,熱端零部件所承受的溫度越來越高,熱載荷作用方式越來越復(fù)雜.防護(hù)涂層承受的熱應(yīng)力一旦超過其強(qiáng)度極限,表面將萌生邊緣裂紋.在循環(huán)熱載荷的作用下,邊緣裂紋將繼續(xù)擴(kuò)展,導(dǎo)致涂層剝落甚至貫穿至基體內(nèi)部[2],嚴(yán)重影響涂層的防護(hù)性能.基于斷裂力學(xué)的觀點,可采用裂紋尖端的應(yīng)力強(qiáng)度因子表征熱致裂紋擴(kuò)展驅(qū)動力,從而有效預(yù)測和評估裂紋的擴(kuò)展行為.
作為導(dǎo)熱理論的本構(gòu)方程,經(jīng)典傅里葉定律建立了熱流密度與溫度梯度之間的線性關(guān)系,在處理熱致斷裂問題上發(fā)揮了重要作用[3-4].然而,傅里葉導(dǎo)熱方程隱含了無限大熱傳播速度的假定,適用于熱作用時間較長的傳熱過程.而對于極端傳熱過程(極高/低溫、微時空尺寸等),經(jīng)典傅里葉導(dǎo)熱模型的求解結(jié)果會出現(xiàn)偏差.Maurer 等[5]對薄板進(jìn)行脈沖激光加熱后,表面溫度比傅里葉定律模型預(yù)測的溫度高約300 K;Tzou[6]的實驗表明,當(dāng)時間尺度低至ps/Hs 或空間尺度低至μm/nm 時,傅里葉定律模型不再精確;Qiu 等[7]通對金屬涂層脈沖激光沖擊,實驗得到的反射率變化與采用傅里葉定律計算所得的數(shù)據(jù)相背離.
在對傅里葉導(dǎo)熱模型進(jìn)行修正以描述上述瞬態(tài)熱傳導(dǎo)過程的研究中,考慮時滯效應(yīng)的C-V (Cattaneo-Vernotte)模型[8]得到了廣泛使用.Wang 等[9-11]基于此模型得到了溫度場和熱應(yīng)力場的解析解,評估了帶裂紋薄板和圓柱體的熱沖擊性能;Chen 等[12]基于C-V 模型獲得了熱彈性半無限板裂紋附近的瞬態(tài)熱應(yīng)力場,并研究了相應(yīng)的熱致斷裂問題.然而,C-V 模型自身也存在局限性,Ahmadikia 等[13]分析了受周期性邊界條件的C-V 傳熱模型,調(diào)節(jié)熱沖擊位置等參數(shù),發(fā)現(xiàn)弛豫時間對溫度場的分布有較大影響,并指出在某些不穩(wěn)定邊界條件下,雙曲導(dǎo)熱方程可能違背熱力學(xué)第二定律.
最近,分?jǐn)?shù)階微積分已成功應(yīng)用于眾多工程領(lǐng)域.大量文獻(xiàn)[14-20]表明: 分?jǐn)?shù)階微分算子定義中的積分項充分體現(xiàn)了系統(tǒng)函數(shù)發(fā)展的歷史依賴性,是描述極端傳熱現(xiàn)象的有力工具,相比經(jīng)典傅里葉定律和C-V 導(dǎo)熱理論,分?jǐn)?shù)階導(dǎo)熱模型可以更精確地描述極端傳熱機(jī)理,并有效解決奇異性問題.通過對分?jǐn)?shù)階的階次α在有效范圍內(nèi)的取值,可以將熱擴(kuò)散形式分為亞擴(kuò)散(0 <α<1)、正常擴(kuò)散(α=1)、超擴(kuò)散(1 <α<2)和彈道擴(kuò)散(α=2)[21].基于分?jǐn)?shù)階導(dǎo)熱模型,Povstenko 等[22-24]提出準(zhǔn)靜態(tài)非耦合的熱彈性理論,研究了各類熱載荷下薄板、圓柱體和球體等均質(zhì)材料的傳熱機(jī)理.此外,通過數(shù)值分析解釋了亞擴(kuò)散和超擴(kuò)散之間的區(qū)別,并指出當(dāng)α=1 時,分?jǐn)?shù)階導(dǎo)熱可以退化為傳統(tǒng)傅里葉導(dǎo)熱模型.文獻(xiàn)[25-28]基于分?jǐn)?shù)階導(dǎo)熱模型分析了預(yù)置裂紋均質(zhì)材料的熱致斷裂問題,并評估了分?jǐn)?shù)階導(dǎo)熱模型的可靠性.
雖然基于分?jǐn)?shù)階導(dǎo)熱模型的熱致斷裂問題已有研究,但其研究對象主要集中于均質(zhì)材料,對涂層-基體復(fù)合材料的相關(guān)研究相對較少.因此,本文擬基于分?jǐn)?shù)階導(dǎo)熱理論,以熱應(yīng)力強(qiáng)度因子表征邊緣裂紋的擴(kuò)展驅(qū)動力,研究熱流脈沖作用下涂層邊緣裂紋的擴(kuò)展行為,為涂層的安全服役及可靠性評估提供理論支持.
考慮涂層-基體系統(tǒng),涂層厚度為h1,基體厚度為h2,脈沖熱流作用于涂層表面,如圖1(a)所示.建立笛卡爾坐標(biāo)系,x軸設(shè)定為厚度方向,y軸設(shè)定在涂層和基體的界面處.為簡單起見,假定兩種材料都是均勻的、各向同性的,涂層和基體完好接觸,忽略熱彈性耦合效應(yīng)、慣性效應(yīng)以及溫度對力學(xué)性能的影響.
圖1 涂層-基體模型Fig.1 Coating-substrate model
由于對稱性,溫度場與y坐標(biāo)無關(guān),基于Caputo分?jǐn)?shù)階導(dǎo)熱理論[29],建立涂層-基體的溫度場控制方程如下
式中,θi=Ti-T0,其中Ti(i=1,2)分別為涂層和基體的溫度,T0為初始溫度.D為熱擴(kuò)散率,其量綱為m2/sα.下標(biāo)i=1,2 分別對應(yīng)涂層和基體,t為時間.α為分?jǐn)?shù)階階次,當(dāng)α=1 時,分?jǐn)?shù)階導(dǎo)熱模型退化為經(jīng)典傅里葉導(dǎo)熱模型.
Caputo 分?jǐn)?shù)階算子由下式給出
式中,Γ(·)為Gamma 函數(shù),n為(α+1)的整數(shù)部分,熱流密度與溫度的關(guān)系滿足
式中,q為熱流密度,k為導(dǎo)熱系數(shù),為R-L(Riemann-Liouville)算子,定義如下[30]
初始條件為
邊界條件由下式給出
式中,Q0為冷卻熱流密度,δ (t)為狄拉克函數(shù).
假設(shè)界面完好接觸,溫度場和熱流密度在界面處均連續(xù),即
根據(jù)分?jǐn)?shù)階導(dǎo)數(shù)的拉普拉斯積分變換法則[30],對時間t作變換(t→s)
式中,θ?(x,s)為拉普拉斯變換后的溫度場.
考慮初始條件式(5)和式(6),式(1)拉普拉斯積分變換后如下
由式(3)和式(10)得
對式(13)作拉普拉斯積分變換得
式中,為拉普拉斯變換后的熱流密度.
同理可得
利用有限余弦積分變換法則[31]對式(12)作空間變量x的變換(x→ξm)得
涂層區(qū)域上(-h1≤x≤0)變換如下
將式(7)、式(8)、式(14)和式(15)代入式(17)和式(18)可得
將式(19)和式(20)代入式(16)并整理后得
對式(21)有限余弦逆變換可得
為確定界面處的熱流密度q(0,t),采用三角級數(shù)求和公式[32],即
利用上述公式對式(22)求和可得
由式(9)可得
聯(lián)立式(24)、式(25)和式(26)得到
基于Tauberian 近似定理,對熱傳導(dǎo)求解問題中較小的時間值,拉普拉斯變換中的指數(shù)項相比常數(shù)項會迅速衰減,因此指數(shù)項可以忽略不計,以得到一個更容易數(shù)值處理的解析解[33],由此可認(rèn)為
將式(27)近似處理可得
式(24)和式(25)作同樣近似處理,并將式 (8)代入得
最后,利用Mainardi 拉普拉斯積分逆變換法則[34]對式(29)作逆變換得到溫度場如下
式中,M(α;z)是Mainardi 函數(shù),即
涂層表面無面力(見圖1(a)),因此 σxx=0.由熱彈性平衡方程,可得下式
式中,σ 為應(yīng)力,E為彈性模量,μ 為泊松比,ε 為應(yīng)變,p為熱膨脹系數(shù).
式中,ε0和 ρ 均為待定常數(shù),可根據(jù)熱應(yīng)力自平衡條件確定,即
平面應(yīng)變條件下,涂層和基體在y方向的熱應(yīng)力為[35]
式中,ε0(t)和 ρ (t)均為與時間t相關(guān)的待定參數(shù),表達(dá)式見附錄A 中的式(A1).
涂層承受的熱應(yīng)力一旦超過其強(qiáng)度極限,其表面將萌生邊緣裂紋.必須指出,多數(shù)情況下涂層表面分布群體裂紋,為簡單起見,本文僅考慮單條邊緣裂紋,該裂紋與涂層表面垂直,且裂紋尖端位于涂層之內(nèi)(a
根據(jù)疊加原理及權(quán)函數(shù)理論[36-37],應(yīng)力強(qiáng)度因子KI可由權(quán)函數(shù)與無裂紋體在假想裂紋位置應(yīng)力乘積的積分得到,即
式中,η 軸的原點位于裂紋口部,如圖1(b)所示,η=x+h1,σ (η)=-σ1yy,w(η,a)為權(quán)函數(shù),其級數(shù)表達(dá)式為[38]
式中,Mj為待定系數(shù),通常只需取上式級數(shù)中的前3 項就具有足夠精度[38].因此,M1,M2和M3可通過兩個參考載荷作用下的應(yīng)力強(qiáng)度因子和一個附加條件確定.
附加條件常取權(quán)函數(shù)在裂紋口部(η=0)曲率為0[38],即
兩種參考載荷可分別取作用在裂紋口部的一對集中力P和作用在裂紋面上的均布載荷σ,因此
式中,KIP和KIσ為兩種參考載荷作用下的應(yīng)力強(qiáng)度因子,FP和Fσ為形狀因子.
聯(lián)立式(37)~式(40),得到M1,M2和M3表達(dá)式如下[38]
式(41)中FP和Fσ為待定參數(shù),可以通過有限元方法計算兩種參考載荷下的應(yīng)力強(qiáng)度因子確定.考慮到涂層厚度遠(yuǎn)小于基體厚度,文獻(xiàn)[4]給出了h2/h1=10時FP和Fσ的擬合表達(dá)式如下
為簡化分析并得到普適性規(guī)律,本文引入下列無量綱參數(shù)
無量綱溫度沿厚度方向的分布規(guī)律如圖2 所示,圖中無量綱時間 κ 依次為0.75,1 和1.25,其他無量綱參數(shù)=1.5,=1.5,h2/h1=10.由圖可知,溫度分布曲線在界面=1 兩側(cè)呈現(xiàn)不同斜率,這是由于涂層和基體材料性能差異,兩者具有不同的熱阻所致.圖2(a)和圖2(b)分別為分?jǐn)?shù)階階次α=0.5(亞擴(kuò)散)和α=1 (經(jīng)典擴(kuò)散)情形,可以看到脈沖熱流加載后,無量綱溫度沿厚度方向降低.圖2(c)表明,當(dāng)分?jǐn)?shù)階階次α=1.5 (超擴(kuò)散)時,溫度沿厚度分布呈現(xiàn)出先增后減的變化趨勢,這是因為熱流脈沖下的熱波在較短時間內(nèi)突破界面并傳遞到了基體,涂層溫度快速降低,與文獻(xiàn)[21]中所描述的分?jǐn)?shù)階熱擴(kuò)散現(xiàn)象一致.
圖2 無量綱溫度 在不同時間 κ 下的分布Fig.2 The distribution of the dimensionless temperature for various timesκ
圖3 無量綱應(yīng)力 在不同時間 κ 下的分布Fig.3 The distribution of the dimensionless stress for various timesκ
無量綱熱應(yīng)力強(qiáng)度因子隨無量綱時間 κ 的變化規(guī)律如圖4 所示,其中無量綱參數(shù)h2/h1=10,=1,=1.2,=1.5,=1.25,=1.可以看出,熱應(yīng)力強(qiáng)度因子在熱流脈沖后的短時間內(nèi)快速增大,達(dá)到峰值后開始下降并逐漸趨于平穩(wěn).相對裂紋長度依次取0.25,0.5 和0.75 時,熱應(yīng)力強(qiáng)度因子的峰值逐漸降低,這是由于熱應(yīng)力沿厚度急劇衰減(見圖3),從而使較長裂紋裂尖附近的應(yīng)力水平降低,說明在熱流脈沖下較短裂紋更易擴(kuò)展.當(dāng)涂層和基體具有相同的材料屬性和傳熱性能時,退化為均質(zhì)薄板,此時應(yīng)力強(qiáng)度因子的變化規(guī)律與均質(zhì)板裂紋的變化趨勢一致[28].對比圖4(a)~圖4(c)可知,隨著分?jǐn)?shù)階階次的增大,熱應(yīng)力強(qiáng)度因子峰值逐漸提高.圖4(c)為超擴(kuò)散(α=1.5)情形下的應(yīng)力強(qiáng)度因子變化規(guī)律,對比圖4(b)可以看出,相同長度邊緣裂紋的應(yīng)力強(qiáng)度因子峰值均明顯提高,表明采用經(jīng)典傅里葉導(dǎo)熱理論將低估熱載荷的破壞程度.
圖4 應(yīng)力強(qiáng)度因子的時間歷程Fig.4 The time history of SIF
本文基于Caputo 分?jǐn)?shù)階導(dǎo)熱模型,采用拉普拉斯積分變換和有限余弦積分變換方法得到了涂層-基體系統(tǒng)在熱脈沖作用下溫度場的封閉解.基于非耦合熱彈性理論和權(quán)函數(shù)法確定涂層邊緣裂紋的熱應(yīng)力強(qiáng)度因子.探討了分?jǐn)?shù)階階次對溫度場、熱應(yīng)力場和熱應(yīng)力強(qiáng)度因子的影響,以及熱應(yīng)力強(qiáng)度因子隨裂紋長度、時間等的變化規(guī)律.主要結(jié)論如下:
(1)分?jǐn)?shù)階階次強(qiáng)烈影響涂層-基體系統(tǒng)的熱響應(yīng)速度以及溫度場和熱應(yīng)力場;
(2)超擴(kuò)散(1 <α<2)情形下的裂紋尖端熱應(yīng)力強(qiáng)度因子峰值高于經(jīng)典擴(kuò)散(α=1)時的峰值,而亞擴(kuò)散(0 <α<1)情形下裂尖熱應(yīng)力強(qiáng)度因子峰值低于經(jīng)典擴(kuò)散時的峰值;
(3)熱流脈沖作用下,涂層邊緣裂紋越長,裂紋尖端應(yīng)力強(qiáng)度因子峰值越小.
附錄A
ε0(t)和 ρ (t)的表達(dá)式如下[35]
式中