熱載荷作用下顆粒增強(qiáng)復(fù)合材料溫升分布數(shù)值模擬

楊萬友， 王家序， 黃彥彥， 周青華， 楊 勇

(四川大學(xué) 空天科學(xué)與工程學(xué)院，成都 610065)

衛(wèi)星、宇宙飛船等航天器常常服役于復(fù)雜空間環(huán)境中，受到高低溫、強(qiáng)輻射、微重力、真空等嚴(yán)酷空間環(huán)境的共同作用，材料性能被嚴(yán)重影響.尤其是寬范圍的高低溫交變，將導(dǎo)致材料發(fā)生熱變形和產(chǎn)生內(nèi)部熱應(yīng)力，造成材料性能退化和航天器可靠性降低.此外，針對航天器使用的高性能顆粒增強(qiáng)復(fù)合材料，熱量在傳遞的過程中還會受到增強(qiáng)體的影響，使材料表面及內(nèi)部溫度場產(chǎn)生擾動，從而造成表面熱變形不協(xié)調(diào)、內(nèi)部熱應(yīng)力集中等不利后果，進(jìn)一步影響航天器的可靠性[1-2].

顆粒增強(qiáng)復(fù)合材料作為一種典型的工程材料，其目的是通過往基體材料中添加不同材料屬性的增強(qiáng)體以提升基體材料性能.由于熱傳導(dǎo)屬性的差異，增強(qiáng)體的存在將會影響熱載荷作用下顆粒增強(qiáng)復(fù)合材料的熱平衡.對于材料熱平衡問題的理論分析方法，國內(nèi)外學(xué)者開展了許多相關(guān)研究.早期的理論研究多基于均質(zhì)材料這一假設(shè)，因此可用解析方法直接進(jìn)行求解.如Carslaw等[3]給出了移動熱源作用下半空間表面溫升解析解.其后圓形熱源[4]、正方形和橢圓形熱源[5]作用下半空間表面溫升分布的解析公式被其他學(xué)者們推導(dǎo)出來.Liu等[1,6]提出了點(diǎn)力或者點(diǎn)熱源作用下求解半空間3種狀態(tài)(暫瞬態(tài)、恒瞬態(tài)和穩(wěn)態(tài))下溫升、熱應(yīng)力[1]和表面熱位移[6]的頻率響應(yīng)函數(shù)，這避免了采用Green函數(shù)計(jì)算溫度場時函數(shù)存在奇點(diǎn)的問題[7-8].

然而顆粒增強(qiáng)復(fù)合材料是由多種性質(zhì)不同的材料復(fù)合而成，呈現(xiàn)出典型的非均質(zhì)特性.由于其復(fù)雜的內(nèi)部微觀結(jié)構(gòu)，難以找到解析解求解相關(guān)的熱平衡等問題，所以，有限元法(FEM)和半解析法(SAM)等數(shù)值方法逐漸被廣大學(xué)者采用來求解非均質(zhì)材料問題.如Bazyar等[9]利用比例邊界FEM求解了各向異性非均質(zhì)材料的二維熱傳導(dǎo)問題；Mortazavi等[10]建立了有限元模型來預(yù)測各向異性隨機(jī)雙相復(fù)合材料的有效導(dǎo)熱系數(shù)和彈性模量，并將之與其他模型進(jìn)行對比分析研究.此外，余天堂等[11]提出了以導(dǎo)熱系數(shù)為基本參數(shù)的熱傳導(dǎo)擴(kuò)展FEM來研究非均質(zhì)材料中的熱傳導(dǎo)問題；陳康等[12]給出了一種適用于梯度復(fù)合材料熱傳導(dǎo)分析的梯度單元法，結(jié)果表明，梯度單元和均勻單元得到的溫度場基本一致.類似研究還包括文獻(xiàn)[13-15].盡管FEM在各種商業(yè)軟件支持下廣泛應(yīng)用于顆粒增強(qiáng)復(fù)合材料熱傳導(dǎo)分析，但是該方法具有材料微觀結(jié)構(gòu)建模困難、計(jì)算耗時長、結(jié)果過度依賴網(wǎng)格尺寸和類型等缺點(diǎn).

為了克服上述FEM的缺點(diǎn)，學(xué)者們提出利用SAM來求解與顆粒增強(qiáng)復(fù)合材料相關(guān)的科學(xué)問題.在Eshelby[16]提出等效夾雜方法(EIM)求解非均質(zhì)材料彈性場這一開創(chuàng)性工作之后，SAM被逐漸發(fā)展、豐富和完善[17-22]，可求解含任意形狀、分布增強(qiáng)體顆粒增強(qiáng)復(fù)合材料的粗糙面彈性接觸問題.然而對于顆粒增強(qiáng)復(fù)合材料熱平衡分析，目前鮮有文獻(xiàn)報(bào)道.Hatta等[23]類比Eshelby[16]提出的彈性場等效夾雜方法，將應(yīng)力對應(yīng)熱流密度、應(yīng)變對應(yīng)溫度梯度、彈性模量對應(yīng)導(dǎo)熱系數(shù)，提出了適用于求解含圓/球形、橢圓/球形增強(qiáng)體顆粒增強(qiáng)復(fù)合材料的穩(wěn)態(tài)熱傳導(dǎo)全空間等效夾雜方法.然而，全空間解只能處理遠(yuǎn)端熱載荷問題，真實(shí)熱載荷作用下顆粒增強(qiáng)復(fù)合材料的熱平衡問題則需要在半空間條件下求解.

本文將顆粒增強(qiáng)復(fù)合材料中增強(qiáng)體設(shè)定為橢球，利用鏡像法[24]完成熱傳導(dǎo)等效夾雜方法從全空間到半空間的轉(zhuǎn)換，用于研究熱載荷作用下顆粒增強(qiáng)復(fù)合材料內(nèi)部和表面溫升分布規(guī)律.分析了不同導(dǎo)熱系數(shù)、深度和位置朝向的增強(qiáng)體對顆粒增強(qiáng)復(fù)合材料熱平衡行為的影響，并通過對含任意分布球形增強(qiáng)體顆粒增強(qiáng)復(fù)合材料溫升分布進(jìn)行求解，驗(yàn)證了算法處理真實(shí)顆粒增強(qiáng)復(fù)合材料熱平衡分析的能力.

1 基本理論和方程

(1)

圖1 穩(wěn)態(tài)熱傳導(dǎo)等效夾雜方法Fig.1 EIM for steady state heat conduction

圖2 鏡像法基本原理Fig.2 Basic principle of method of images

采用鏡像法[24]，實(shí)現(xiàn)穩(wěn)態(tài)熱傳導(dǎo)EIM從全空間到半空間的轉(zhuǎn)換.基于鏡像法原理(見圖2)，半空間夾雜引起的溫度梯度量是兩個全空間夾雜引起的溫度梯度量和對稱面法向熱流引起的溫度梯度量之和，其求解方法詳見文獻(xiàn)[1,23].

進(jìn)一步引入類似彈性場中的Eshelby張量Sij，文獻(xiàn)[25]已給出其推導(dǎo)過程及具體形式，此處不再贅述.利用Eshelby張量Sij關(guān)聯(lián)擾動溫度梯度與本征溫度梯度：

(2)

將上式代入式(1)整理得

(3)

(4)

式中：Φ,i為牛頓勢函數(shù)Φ的一階導(dǎo)數(shù)，詳情參見文獻(xiàn)[25]；xi為目標(biāo)溫度坐標(biāo)點(diǎn)到增強(qiáng)體單元中心坐標(biāo)的向量.顆粒增強(qiáng)復(fù)合材料總溫升Tt為均質(zhì)材料溫升T與增強(qiáng)體擾動溫升Td之和[26]：

Tt=T+Td

(5)

(6)

式中：Sij為關(guān)聯(lián)第j個等效夾雜相應(yīng)溫度梯度與第i個本征溫度梯度之間相互影響的系數(shù)；n為增強(qiáng)體的數(shù)目.

2 顆粒增強(qiáng)復(fù)合材料熱平衡算法驗(yàn)證

利用熱載荷作用下含單橢球增強(qiáng)體顆粒增強(qiáng)復(fù)合材料模型驗(yàn)證本文所提算法的有效性.模型示意圖如圖1(a)所示，其中橢球形增強(qiáng)體尺寸設(shè)置為a=2b=2c=r/2，并令其中心距離表面h=r.外部加載熱流密度可表示為

(7)

模型計(jì)算區(qū)域大小設(shè)置為6r×6r×3r，等距離散為128×128×64個網(wǎng)格點(diǎn)，利用上述算法對不同導(dǎo)熱系數(shù)(K*=4Km或K*=Km/4，其中，基體材料導(dǎo)熱系數(shù)Km=50.2 W/(m·K)[1])增強(qiáng)體進(jìn)行計(jì)算.同時，利用Abaqus軟件中軸對稱方式建立同樣條件的有限元驗(yàn)證模型，如圖3所示.另外值得注意的是，第1節(jié)中已直接或間接說明了本文所提模型/算法與相關(guān)的有限元分析模型基于下列3點(diǎn)假設(shè)：① 顆粒增強(qiáng)復(fù)合材料是半無限體；② 熱傳導(dǎo)分析是在穩(wěn)態(tài)時進(jìn)行；③ 顆粒增強(qiáng)體之間不重疊.

圖4所示為利用上述兩種方法計(jì)算得到的顆粒增強(qiáng)復(fù)合材料下表面溫升分布云圖，NT11表示節(jié)點(diǎn)溫度輸出變量.當(dāng)增強(qiáng)體與基體導(dǎo)熱系數(shù)關(guān)系為K*=4Km時，增強(qiáng)體內(nèi)部溫度等值線分布相較基體稀疏，溫度變化率較??；而K*=Km/4時，增強(qiáng)體內(nèi)部溫度等值線分布較基體緊湊，溫度變化率較大.圖5所示為利用兩種方法計(jì)算得到的模型結(jié)果沿z軸的溫升分布.從圖5中可以看出，增強(qiáng)體單元對基體溫度場有一處明顯的擾動，并且K*=4Km時會在增強(qiáng)體距離表面近端(z/r=0.5)以及遠(yuǎn)端(z/r=1.5)附近分別造成相對均質(zhì)材料的溫度降低和升高，而增強(qiáng)體材料屬性為K*=Km/4時造成的影響與前者相反且相對較弱.此外，利用FEM計(jì)算的結(jié)果與本文算法模擬的結(jié)果在數(shù)值上能夠較好地吻合，驗(yàn)證了所提算法的有效性.

圖3 含單一橢球增強(qiáng)體二維有限元對稱穩(wěn)態(tài)熱傳導(dǎo)模型Fig.3 2D axisymmetric FEM including an ellipsoidal reinforcement for steady state heat conduction

圖4 含單一橢球增強(qiáng)體下表面溫升分布Fig.4 Lower surface temperature rise distribution including an ellipsoidal reinforcement

圖5 含單一橢球增強(qiáng)體沿z軸的溫升分布Fig.5 Temperature rise distribution along the z axis including an ellipsoidal reinforcement

3 計(jì)算結(jié)果和討論

采用圖1(a)的模型，除了特別說明外，其余參數(shù)設(shè)置參見第2節(jié)，以橢球形增強(qiáng)體導(dǎo)熱系數(shù)、位置、朝向、數(shù)目為研究參數(shù)，分析其對熱載荷作用下顆粒增強(qiáng)復(fù)合材料溫升分布影響.

3.1 單個增強(qiáng)體熱平衡分析

設(shè)置橢球形增強(qiáng)體尺寸參數(shù)為a=1.25b=1.25c=r，增強(qiáng)體中心至表面的距離h=0.75r.由于橢球本身各軸長短不一，在空間中有著復(fù)雜的位置關(guān)系.如圖6所示，可采用歐拉角α、β和γ來描述橢球增強(qiáng)體所處的空間位置.目標(biāo)點(diǎn)從全局坐標(biāo)到局部坐標(biāo)的轉(zhuǎn)換可表示為

Rxyz(m)=Rxyz(m)+

a(m,n)[Pxyz(m)-Cxyz(m)]

(8)

m,n=1,2,3

圖6 歐拉角和坐標(biāo)轉(zhuǎn)換關(guān)系Fig.6 Relationship of Euler angles and coordinate transformation

式中：m,n分別代表x,y,z；Rxyz(m)和Pxyz(m)分別為目標(biāo)點(diǎn)的局部坐標(biāo)和全局坐標(biāo)；Cxyz(m)為橢球中心的全局坐標(biāo);a(m,n)為對應(yīng)的旋轉(zhuǎn)因子，具體可表示為

(9)

由上述歐拉角的定義，圖1(a)模型中橢球增強(qiáng)體各歐拉角分別為α=0°、β=90° 和γ=90°.不論增強(qiáng)體處于什么位置，其對溫度場或多或少都有擾動，且這種擾動直接表現(xiàn)為溫度的變化.本節(jié)利用表面溫度變化率η來分析各參數(shù)對材料熱平衡行為的影響.具體表達(dá)式為

(10)

圖7 不同導(dǎo)熱系數(shù)增強(qiáng)體表面溫升分布及溫度變化率Fig.7 Distribution and change rate of surface temperature rise including an ellipsoidal reinforcement with different thermal conduction properties

圖8 不同導(dǎo)熱系數(shù)增強(qiáng)體下表面溫升分布Fig.8 Lower surface temperature rise distribution including an ellipsoidal reinforcement with different thermal conduction properties

3.1.1導(dǎo)熱系數(shù) 導(dǎo)熱系數(shù)是影響材料熱平衡的關(guān)鍵因素，本節(jié)將通過改變增強(qiáng)體與基體導(dǎo)熱系數(shù)之比λ(λ=K*/Km或λ=Km/K*)使之從1～20之間變化，以探究顆粒增強(qiáng)復(fù)合材料內(nèi)部熱平衡規(guī)律，結(jié)果如圖7和8所示.從圖7(a)中可以看出，當(dāng)λ=K*/Km時，隨著λ的增大，表面溫升分布趨于平緩，表面溫度值逐漸減?。划?dāng)λ=Km/K*時，則得到與前者相反的結(jié)論.圖7(a)中還給出了K*=20Km、K*=Km/20和K*=Km3種情況下溫度標(biāo)尺相同的表面溫升三維等值線圖，可直觀地看出表面溫升情況.圖7(b)結(jié)果表明，隨著λ的增大，η先快速增長后增速逐漸放緩，且λ=K*/Km時相對λ=Km/K*增速更為顯著.圖8為圖7(a)情況對應(yīng)的下表面沿z軸的溫升分布.圖中曲線表明增強(qiáng)體距離表面近端和遠(yuǎn)端位置附近都會對溫度場有較大擾動，且λ越大，擾動越明顯.產(chǎn)生上述現(xiàn)象的原因在于增強(qiáng)體的導(dǎo)熱系數(shù)較基體大時，則熱量在增強(qiáng)體中集中，導(dǎo)致在增強(qiáng)體距表面近端處溫度相對均質(zhì)材料要低，而遠(yuǎn)端溫度相對均質(zhì)材料要高.增強(qiáng)體導(dǎo)熱系數(shù)較基體小時則相反.

3.1.2深度 由于材料受摩擦熱等外部熱載荷作用時，是從表面往內(nèi)部傳遞的，所以，顆粒增強(qiáng)復(fù)合材料內(nèi)部熱平衡往往與增強(qiáng)體的深度位置密切相關(guān).設(shè)置橢球增強(qiáng)體中心到表面的距離h∈[0.6r, 2r]，研究不同深度增強(qiáng)體對顆粒增強(qiáng)復(fù)合材料熱平衡的影響規(guī)律.圖9所示為含不同深度增強(qiáng)體顆粒增強(qiáng)復(fù)合材料表面溫升分布及溫度變化率.圖9(a)中表面溫升分布曲線表明，不論增強(qiáng)體導(dǎo)熱系數(shù)大于或小于基體，隨著增強(qiáng)體深度的增加，表面溫升分布都將越來越接近均質(zhì)材料.在h/r=0.6時，K*=20Km和K*=Km/20兩種情況下基體表面溫升分布可從相應(yīng)的三維等值線圖直觀看出.從圖9(b)中可以看出，隨著增強(qiáng)體深度的增加，無論其導(dǎo)熱系數(shù)大于或小于基體，顆粒增強(qiáng)復(fù)合材料表面溫度變化率都快速下降，逐漸趨近于0.上述結(jié)果表明，隨著增強(qiáng)體深度的增加，表面熱載荷作用對顆粒增強(qiáng)復(fù)合材料內(nèi)部及表面溫度場擾動逐漸減小.

圖9 不同深度位置增強(qiáng)體表面溫升分布及溫度變化率Fig.9 Distribution and change rate of surface temperature rise including an ellipsoidal reinforcement at different depth locations

3.1.3朝向角度 顆粒增強(qiáng)復(fù)合材料中類纖維狀增強(qiáng)體空間朝向具有隨機(jī)性，且增強(qiáng)體的朝向角度已被證實(shí)對應(yīng)力場會產(chǎn)生影響[17].本小節(jié)將改變其朝向角度，保持α=0° 和β=90° 不變，使歐拉角γ在0°～ 90° 之間變化，探究其對顆粒增強(qiáng)復(fù)合材料溫度場的影響規(guī)律.圖10給出了γ分別為0°、45° 和90° 時顆粒增強(qiáng)復(fù)合材料下表面溫升分布云圖，此時K*=4Km.由圖可知，增強(qiáng)體不同空間朝向?qū)Σ牧蠝厣植籍a(chǎn)生較大影響.從圖11顆粒增強(qiáng)復(fù)合材料表面溫度變化率曲線中可知，隨著γ的逐漸增大，η先減小后增大，且γ=45° 時η最小.

圖10 不同朝向角度橢球形增強(qiáng)體下表面溫升分布Fig.10 Lower surface temperature rise distribution including an ellipsoidal reinforcement with different facing angles

圖11 不同朝向角度橢球形增強(qiáng)體表面溫度變化率Fig.11 Distribution and change rate of surface temperature including an ellipsoidal reinforcement with different facing angles

3.2 多增強(qiáng)體熱平衡分析

增強(qiáng)體往往隨機(jī)分布于顆粒增強(qiáng)復(fù)合材料中.本小節(jié)考慮多個增強(qiáng)體之間的相互作用，進(jìn)一步分析其對顆粒增強(qiáng)復(fù)合材料溫度場的影響.

圖12(a)為含雙球形增強(qiáng)體顆粒增強(qiáng)復(fù)合材料模型，增強(qiáng)體中心位于Oxz平面上，沿z軸對稱分布，且其中心到表面的距離h=0.6r.球形增強(qiáng)體直徑d=r，間距為l.增強(qiáng)體的導(dǎo)熱系數(shù)為K*=4Km和K*=Km/4，其余參數(shù)保持不變.改變間距l(xiāng)，使其在r～3r之間變化，研究顆粒增強(qiáng)復(fù)合材料下表面溫升分布規(guī)律.

圖12(b)和12(c)所示為利用本文所提算法求出的含不同間距雙球形增強(qiáng)體(K*=4Km)顆粒增強(qiáng)復(fù)合材料下表面溫升分布云圖.從圖中可以直觀看出不同間距的增強(qiáng)體對基體溫度場的影響程度不同，而且當(dāng)l增大時，會使基體溫度場向周邊擴(kuò)散.

圖13所示為不同間距雙球形增強(qiáng)體顆粒增強(qiáng)復(fù)合材料表面溫度變化率.圖中曲線表明，隨著l增大，η逐漸降低，且兩種材料增強(qiáng)體造成的η差異也在逐漸減小.

圖14(a)所示為中心位置坐標(biāo)(x,y,z)服從正態(tài)分布的200個半徑為0.2r球形增強(qiáng)體顆粒增強(qiáng)復(fù)合材料模型，即

x～N(0,r2)

y～N(0,r2)

z～N(r,r2)

圖12 含雙球形增強(qiáng)體模型及下表面溫升分布Fig.12 Model including two spherical reinforcements and its lower surface temperature rise distribution

圖13 不同間距雙球形增強(qiáng)體表面溫度變化率Fig.13 Change rate of surface temperature including two spherical reinforcements with different distances

圖14 含正態(tài)分布球形增強(qiáng)體模型及下表面溫升分布Fig.14 Model including spherical reinforcements in normal distribution and its lower surface temperature rise distribution

增強(qiáng)體導(dǎo)熱系數(shù)為K*=20Km和K*=Km/20，其余參數(shù)參照前例，計(jì)算結(jié)果如圖14(b)和14(c)所示.從圖中可以看出，隨機(jī)分布的增強(qiáng)體對材料下表面溫升分布擾動較大，且位置接近表層的增強(qiáng)體對表面溫度及周圍溫度影響較為顯著，而較深的增強(qiáng)體對溫度場的擾動則相對較小.對比K*=20Km和K*=Km/20兩種情況，可發(fā)現(xiàn)后者溫度變化范圍更廣，且具有較大的溫升，說明即使是將增強(qiáng)體分散于基體中，其導(dǎo)熱系數(shù)對基體溫升分布也具有較大影響.

4 結(jié)論

(1) 本文提出了一種熱載荷作用下顆粒增強(qiáng)復(fù)合材料熱平衡數(shù)值建模方法.對于單橢球增強(qiáng)體模型，增強(qiáng)體導(dǎo)熱系數(shù)與基體差異越大，越容易造成表面溫度變化率增大.因此在顆粒增強(qiáng)復(fù)合材料制備過程中，增強(qiáng)體熱屬性不宜與基體差異太大，否則影響表面溫度變化率，從而導(dǎo)致材料表面性能退化.

(2) 單橢球形增強(qiáng)體越靠近表面高溫度區(qū)時，其對基體表面溫度擾動越大.此外，單橢球形增強(qiáng)體朝向會對溫度場產(chǎn)生較大影響，在本文設(shè)定情況下，歐拉角γ在45° 左右時對表面溫度的影響最小.工程實(shí)際中，顆粒增強(qiáng)復(fù)合材料增強(qiáng)體的分布應(yīng)盡量遠(yuǎn)離受熱載荷表面區(qū)域，且橢球/纖維狀增強(qiáng)體歐拉角γ應(yīng)在45° 左右.

(3) 雙球形增強(qiáng)體間距增大，表面溫度變化率逐漸降低且由不同導(dǎo)熱系數(shù)增強(qiáng)體造成的顆粒增強(qiáng)復(fù)合材料表面溫度變化率差異也在逐漸減小.分布于顆粒增強(qiáng)復(fù)合材料中的增強(qiáng)體會對基體溫度場造成擾動，且越靠近表面的增強(qiáng)體造成的擾動越顯著.

本文所提材料溫升分布求解模型可有效考慮材料微觀結(jié)構(gòu)對熱載荷作用下材料內(nèi)部熱平衡的影響，為顆粒增強(qiáng)復(fù)合材料及其他材料散熱、吸熱等功能需求為導(dǎo)向的內(nèi)部微觀結(jié)構(gòu)優(yōu)化設(shè)計(jì)提供理論指導(dǎo).耦合材料彈、塑、電、磁等不同場的求解模型，可進(jìn)一步求解不同外載聯(lián)合作用下非均質(zhì)材料復(fù)雜工程力學(xué)問題.