許楊健，劉文芝

(河北工程大學(xué)土木工程學(xué)院，河北邯鄲056038)

功能梯度材料(Functionally Graded Materials，簡稱FGM)是一種新型的特殊復(fù)合材料［1－2］，優(yōu)于傳統(tǒng)復(fù)合材料的特性，已成為材料領(lǐng)域研究的重點及多學(xué)科交叉領(lǐng)域的研究熱點［3－7］。因 FGM在超高溫工作環(huán)境中的應(yīng)用日益廣泛，故分析該材料組成物體的熱應(yīng)力場特別是瞬態(tài)問題更具實用價值。

國外Obata［8］等采用攝動法對加熱、冷卻過程中的常物性FGM板瞬態(tài)熱應(yīng)力進(jìn)行了研究，但是這種方法過于繁復(fù)，不便于工程應(yīng)用。在前面的研究工作基礎(chǔ)上，本文采用有限元法，對處在不同的位移邊界條件下的常物性2D－FGM板冷卻瞬態(tài)熱應(yīng)力問題進(jìn)行研究分析，其結(jié)果對FGM的生產(chǎn)和應(yīng)用具有重要的指導(dǎo)意義。

1 研究模型

如圖1所示，選定二維常物性Al1100/Ti－6Al－4V/ZrO2FGM板作為分析模型，其平面結(jié)構(gòu)長度l=40 mm、厚度b=10 mm，且假設(shè):(1)該FGM 平面區(qū)域具有長度x與厚度y方向任意分布及連續(xù)變化的材料物性性質(zhì);(2)FGM平面區(qū)域的溫度初始條件和加熱邊界條件為:初始溫度為常溫T0，當(dāng)時間t＞0時，在區(qū)域的四周表面外側(cè)分別施加第一類加熱邊界溫度Ta、Tb、Tc及Td并保持不變;(3)分析模型前后表面絕熱，內(nèi)部無內(nèi)熱源，且不考慮熱固耦合項作用影響。

圖中，k，ρ，cρ，E，ν以及α分別是FGM 的熱導(dǎo)率、密度、比熱，彈性模量和線熱膨脹系數(shù)。

1.1 位移邊界條件

為研究外部約束對常物性2D－FGM板平面區(qū)域冷卻瞬態(tài)熱應(yīng)力分布的影響，本文選取4種位移邊界條件，即簡支、一端固定、兩端固定、四周固定，分別如圖 2 中 a、b、c、d 所示。

1.2 平面結(jié)構(gòu)和加熱、冷卻邊界條件

使用網(wǎng)格自動劃分程序?qū)ι鲜鲅芯磕Ｐ推矫鎱^(qū)域進(jìn)行網(wǎng)格剖分，在x方向上劃分80層，即選定81條縱向線段，每一條線段上取21個節(jié)點，將該平面區(qū)域劃分為3 200個單元，1 701個節(jié)點。瞬態(tài)時間步長Δt取0.1，并給出換熱邊界條件:

加熱時，板的周邊第一類換熱邊界條件為

冷卻時，板的周邊第一類換熱邊界條件為

式(1)、式(2)中的Ta，Tb，Tc和Td是平面矩形區(qū)域的四周邊界溫度，T0為初始溫度。

2 熱傳導(dǎo)和熱應(yīng)力有限元基本方程

2.1 熱傳導(dǎo)有限元基本方程

首先將整個平面區(qū)域離散成M個簡單的三角形單元，單元面積為Se使其3個節(jié)點編號i、j、m按逆時針旋轉(zhuǎn)，其中對于邊界單元而言，編號i表示內(nèi)部節(jié)點，編號j和m表示邊界節(jié)點。并且時間過程劃分為n個間隔Δtn，n=2，3，…，N。全部節(jié)點在瞬時tn－1的溫度值為列矩陣。在Δtn內(nèi)，常物性 FGM板瞬態(tài)熱傳導(dǎo)有限元基本方程為［9］

式中H、Q—溫度剛度矩陣和變溫系數(shù)矩陣;T—未知溫度值的列向量。

2.2 熱應(yīng)力有限元基本方程

設(shè)圖1所示常物性2D－FGM板的熱學(xué)性質(zhì)為位置坐標(biāo)x、y的函數(shù)，且在每一坐標(biāo)平面內(nèi)材料具有各向同性的性質(zhì)，則該板的應(yīng)變分量和應(yīng)力分量分別為［10］

式中E(x，y)、a(x，y)、v(x，y)、θ(x，y，t)—彈性模量、線性熱膨脹系數(shù)、泊松比、溫變函數(shù);u和v—單元e內(nèi)任一點(x，y)的橫向位移和縱向位移。

對于熱彈性問題而言，計算應(yīng)力場得先計算溫度場，故此處的網(wǎng)格劃分與熱傳導(dǎo)網(wǎng)格劃分完全相同。則常物性FGM板瞬態(tài)熱應(yīng)力有限元基本方程為［9］

此線性代數(shù)方程的解就是整個求解區(qū)域在tn時刻，所有節(jié)點位移u1，u2，…，uL及v1，v2，…，vL。

3 本研究方法正確性檢驗

3.1 熱傳導(dǎo)正確性檢驗

為了驗證本文有限元程序的正確性，我們應(yīng)用有限元法與分離變量法求解上述特殊非均勻材料平面熱傳導(dǎo)問題的同一算例。設(shè)平面矩形區(qū)域?qū)挒閘=10 mm，厚為b=10 mm，參數(shù)D=E=1，c=d=0.1。并將該區(qū)域劃分為441個節(jié)點，800個單元，整個區(qū)域初始溫度300 K，且區(qū)域四周均作用第一類邊界條件。利用有限元程序計算該材料平面溫度場，并從中取得t=1.0 s，3.0 s，5.0 s 時的溫度場值，并將其同利用分離變量法計算所得的結(jié)果形成對比，如表1所示。

觀察表 1 可知:t=1.0 s、t=3.0 s、t=5.0 s時平面區(qū)域水平中線上11個點的溫度值的分離變量解和有限元解最大誤差分別為 0.44%、0.36%、0.25%，并分別為(7，5)、(5，5)、(9，5)的點。由上述分析可知，利用分離變量法和有限元法求得的兩種溫度解存在的最大誤差都遠(yuǎn)小于5%，是能滿足實際工程需求的。

3.2 熱應(yīng)力正確性檢驗

為檢驗熱應(yīng)力計算部分有限元程序的正確性，我們?nèi)∫缓喼Я鹤鳛樗憷?，設(shè)其模型為2l=100 mm，h/2=5 mm，上邊界作用均布荷載q，兩端作用反力ql保持平衡，且不計體力，溫變函數(shù)θ=θ0(1－4y2/h2)為已知(θ0為常數(shù))，線熱膨脹系數(shù)α=0.000 001/K，泊松比ν=0.3，彈性模量E=200 GPa。利用網(wǎng)格自動劃分程序?qū)⒛Ｐ蛣澐譃?01條縱向線段，21條橫向線段，每條縱向線段間隔為1.0 mm，每條橫向線段間隔為0.5 mm，即網(wǎng)格劃分為2 121個節(jié)點，同時獲得解析解所需要的區(qū)域節(jié)點坐標(biāo)。如圖3所示。

首先給出簡支梁應(yīng)力分量計算表達(dá)式:

表1 溫度場正確性檢驗Tab.1 Testifying temperature distribution accuracy

我們只求解簡支梁僅在溫度荷載作用下的熱應(yīng)力，令q=0，θ0=300。通過式(7)、式(8)、式(9)計算該問題熱應(yīng)力σx的解析解，通過有限元法計算得到其數(shù)值解。給出簡支梁x坐標(biāo)為0.0的縱向線(y軸)上熱應(yīng)力的解析解和有限元解，并將兩者的計算結(jié)果進(jìn)行對比。如表2所示。

表2 溫變下的簡支梁平面熱應(yīng)力Tab.2 Plane thermal stress of simple supported beam under temperature variation

從表2可以發(fā)現(xiàn):兩者解法的最大誤差為7.37%，顯然誤差較大，但如果在應(yīng)力值都只保留1位有效數(shù)字的情況下，誤差將均為0%，加之應(yīng)力值本身也很小，故仍可認(rèn)為有限元計算結(jié)果是有效的。

4 熱應(yīng)力場分布和分析

如圖1所示模型，取組分形狀分布系數(shù)mx和my為1.0，孔隙率控制參數(shù)Ax和Ay取0.0，nx、ny、zx和zy均取為1.0，其換熱邊界條件為前面所述式(1)和式(2)，按圖2所示的外部約束，即簡支、一端固定、兩端固定和四端固定4種位移邊界條件，得到t=1.0 s時的冷卻瞬態(tài)熱應(yīng)力場分布圖形，如圖4所示。

我們知道，位移邊界的不同是不會引起溫度場變化的，那么在各種計算相同，換熱邊界條件相同的情況下，常物性2D－FGM平面區(qū)域的溫度場分布在不同位移約束條件下應(yīng)當(dāng)是保持一致的。并且根據(jù)熱彈性力學(xué)理論，若FGM平面區(qū)域的左右兩個邊界是自由邊界，那么邊界上的熱應(yīng)力σx=0。

分析圖4可知:冷卻情況下，位移邊界條件對常物性2D－FGM平面區(qū)域冷卻瞬態(tài)熱應(yīng)力分布影響頗大。

觀察圖4(a)和4(b)可知，在簡支情況下，平面區(qū)域的左右兩個邊界都是自由邊界，在一端固定情況下，平面區(qū)域的右邊界是自由邊界，因此，這些邊界上的熱應(yīng)力σx=0;從圖4(a)變化到圖4(b)，即將左邊界設(shè)置為固定約束后，在平面的左邊界上形成了上中下三個應(yīng)力聚集，但應(yīng)力值幾乎沒有發(fā)生變化，整個區(qū)域既存在壓應(yīng)力又存在拉應(yīng)力。

圖4(c)中，將左右邊界均設(shè)置為固定約束后，其右邊界不僅也產(chǎn)生了相同現(xiàn)象，而且區(qū)域內(nèi)部應(yīng)力分布形狀與數(shù)值都發(fā)生了巨大變化:平面中部應(yīng)力分布曲線更加不平行于x坐標(biāo)方向，區(qū)域內(nèi)部應(yīng)力值全部成為拉應(yīng)力，且數(shù)值增長幅度頗為顯著;將四周均設(shè)置呈固定約束后，如圖4(d)所示，平面四角上的應(yīng)力聚集消失，熱應(yīng)力數(shù)值增長幅度較兩端固定情況時更為明顯。

5 結(jié)論

冷卻過程中，位移邊界條件對常物性2DFGM平面區(qū)域冷卻瞬態(tài)熱應(yīng)力分布影響頗大，將四周均設(shè)置呈固定約束后，平面四角上的應(yīng)力聚集消失，應(yīng)力凸起位置移動到平面右上部位，熱應(yīng)力數(shù)值增長幅度較兩端固定情況時更為明顯。

［1］李信，劉海昌.功能梯度材料的研究現(xiàn)狀及展望［J］.材料導(dǎo)報，2012，26(19):370 －372.

［2］馬濤，趙忠民，劉良祥，等.功能梯度材料的研究進(jìn)展及應(yīng)用前景［J］.化工科技，2012，20(1):71－75.

［3］WUXIANG LIU，ZHENG ZHONG.Three－dimensional thermoelastic analysis of functionally graded plate［J］.Acta Mechanica Solida Sinica，2011，24(3):241 －249.

［4］許楊健，涂代惠，馬士進(jìn).加熱、冷卻下變物性梯度功能材料板瞬態(tài)熱應(yīng)力［J］.機(jī)械強(qiáng)度，2005，27(4):510－517.

［5］M R GOLBAHAR HAGHIGHI，P MALEKZADEH，H RAHIDEH，et al.Inverse Transient heat conduction problems of a multilayered functionally graded cylinder［J］.Numerical Heat Transfer，Part A:Applications.2012(9):717－733.

［6］許楊健，王飛，杜海洋.初始和邊界恒溫時二維梯度板瞬態(tài)溫度場［J］.固體力學(xué)學(xué)報，2014(2):3－6.

［7］許楊健，王飛，杜海洋，等.邊界不同恒溫時功能梯度板平面穩(wěn)態(tài)溫度場［J］.河北工程大學(xué)學(xué)報:自然科學(xué)版，2013，30(2):4 －8.

［8］OBATA Y，NODA N.Unsteady thermal stresses in a functionally gradient material plate(Analysis of one－dimensional unsteady heat transfer problem)［J］.Trans.JSME，1993，59(560):1090 －1096.

［9］王洪綱.熱彈性力學(xué)概論［M］.北京:清華大學(xué)出版社，1989.

［10］嚴(yán)宗達(dá)，王洪禮.熱應(yīng)力［M］.北京:高等教育出版社，1993.