材料參數(shù)隨溫度變化的功能梯度圓筒熱彈性分析

吳雨桐

(合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，安徽 合肥 230009)

0 引 言

功能梯度材料(functionally graded materials,FGM)是指構(gòu)成材料的組成、結(jié)構(gòu)沿厚度方向由一側(cè)向另一側(cè)呈連續(xù)梯度變化,從而使材料性質(zhì)和功能也呈梯度變化的一種新型材料[1],因其良好的力學(xué)與熱學(xué)性能適用于高溫環(huán)境。功能梯度空心圓筒和圓柱結(jié)構(gòu)是工程中的常用結(jié)構(gòu),而溫度荷載也是其在工作過程中常見的荷載類型。因此,研究功能梯度圓筒熱彈性問題對(duì)該類材料結(jié)構(gòu)的設(shè)計(jì)和應(yīng)用具有重要的理論指導(dǎo)意義。

在功能梯度圓筒的熱彈性分析中,由于材料性質(zhì)隨梯度改變,且熱傳導(dǎo)方程中還包含了溫度對(duì)時(shí)間的偏導(dǎo)數(shù),使得求解過程非常復(fù)雜。Shao等[2]采用Laplace變換,借助Galerkin法和級(jí)數(shù)法,得到了功能梯度圓筒熱彈性場(chǎng)的半解析解。Darabseh和Bani-Salameh[3]采用一種基于隱式差分格式的數(shù)值方法,計(jì)算得到了功能梯度圓柱的熱彈性場(chǎng)分布。通過使用一種攝動(dòng)解法,果立成等[4]研究了功能梯度圓柱殼在熱沖擊載荷下的瞬態(tài)溫度場(chǎng)和瞬態(tài)熱應(yīng)力場(chǎng)。盛宏玉等[5-6]將狀態(tài)空間理論與有限差分法相結(jié)合,對(duì)層合材料和組合實(shí)心圓柱進(jìn)行了瞬態(tài)熱傳導(dǎo)分析。

目前在求解功能梯度圓筒瞬態(tài)熱彈性問題時(shí)以積分變換法居多,求解過程煩瑣,且材料參數(shù)隨溫度變化屬于復(fù)雜的非線性問題,很多相關(guān)研究都沒有將其考慮在內(nèi)。馬永斌等[7]和何天虎等[8]雖考慮了溫度對(duì)材料參數(shù)的影響,但為了簡化計(jì)算,在材料參數(shù)與溫度相關(guān)的函數(shù)中用材料初始溫度代替了變化的溫度場(chǎng),這在溫度變化較大的條件下計(jì)算結(jié)果有較大的誤差。本文在空間上使用狀態(tài)空間技術(shù),時(shí)域內(nèi)應(yīng)用有限差分格式,可將某時(shí)刻求得的溫度場(chǎng)應(yīng)用到下一時(shí)刻材料參數(shù)與溫度相關(guān)的函數(shù)中,從而得到該熱彈性問題的準(zhǔn)非線性解,精度能得到進(jìn)一步提升,且方法簡便,計(jì)算效率高。陳新宇[9]應(yīng)用本文方法對(duì)材料參數(shù)只沿徑向坐標(biāo)變化的功能梯度圓筒熱彈性問題進(jìn)行了分析,并與文獻(xiàn)[2]進(jìn)行了對(duì)比,驗(yàn)證了本文方法的正確性和有效性。

1 基本方程

考慮一內(nèi)、外半徑分別為a和b的無限長彈性圓筒,其橫截面與所建坐標(biāo)系如圖1所示。材料的物理參數(shù)κ、λ、α和彈性參數(shù)E和ν分別表示其熱擴(kuò)散率、熱傳導(dǎo)系數(shù)、熱膨脹系數(shù)、彈性模量和泊松比。除泊松比外，其他幾項(xiàng)均為關(guān)于徑向坐標(biāo)r以及溫度場(chǎng)T的函數(shù)。假設(shè)圓筒初始時(shí)刻的溫度T0均勻分布,內(nèi)、外環(huán)境溫度分別為Ta和Tb,且保持不變,則該情況屬于軸對(duì)稱平面應(yīng)變問題。圓筒在環(huán)境溫度作用下,產(chǎn)生一維瞬態(tài)溫度場(chǎng)、位移場(chǎng)和應(yīng)力場(chǎng)。

圖1 功能梯度圓筒的橫截面及坐標(biāo)和分層示意圖

傅里葉方程和能量守恒方程在該問題下的表達(dá)式分別為:

(1)

(2)

式中:q為熱流密度;θ=T-T0,為圓筒的溫度場(chǎng)相對(duì)于初始溫度的溫差。

在不計(jì)體力時(shí),圓筒平衡方程為:

(3)

熱彈性應(yīng)力和位移關(guān)系為

(4)

(5)

式中:σr和σφ分別表示圓筒徑向和環(huán)向應(yīng)力。

初始條件設(shè)為:

θ(r,0)=0

(6)

邊界條件設(shè)為:

σr(a,t)=0,σr(b,t)=0

q(a,t)=βa[θa-θ(a,t)],q(b,t)=βb[θ(b,t)-θb]

(7)

式中:θa=Ta-T0,θb=Tb-T0;βa和βb分別代表圓筒在與外界進(jìn)行對(duì)流換熱時(shí)的傳熱系數(shù)。

2 狀態(tài)方程的建立與求解

由式(4)可得:

(8)

將式(4)和式(5)代入到平衡方程式(3)中,再利用式(8)得:

(9)

于是,式(1)、式(2)、式(8)和式(9)構(gòu)成了一個(gè)偏微分方程組,其矩陣形式為:

(10)

由于方程的系數(shù)矩陣中材料參數(shù)不僅沿圓筒徑向坐標(biāo)變化,還與圓筒的溫度變化有關(guān),式(10)為變系數(shù)非線性偏微分方程組,求解過程十分復(fù)雜,因此需要進(jìn)行如下幾個(gè)方面的處理。

(11)

式中:[θ(r,t)]k表示t=tk時(shí)刻的溫度場(chǎng)。

將式(11)代入式(10)中,在系數(shù)矩陣中取徑向坐標(biāo)r=hi=(ri+ri-1)/2,并在材料參數(shù)與溫度相關(guān)的函數(shù)中代入tk-1時(shí)刻的溫度場(chǎng),最終得到tk時(shí)刻關(guān)于第i層子殼的一階非齊次常微分方程組,即狀態(tài)方程:

(ri-1≤r≤ri,i=1,2,3,…,n)

(12)

對(duì)式(12)積分可得到狀態(tài)向量

[Vi(r,t)]k=Di(r-ri-1)[Vi(ri-1,t)]k+[Hi(r,t)]k,

(ri-1≤r≤ri)

(13)

(14)

利用初始條件(6)和邊界條件(7),由式(12)～(14)可以求出任意時(shí)刻在任意位置處的狀態(tài)變量。

3 算例和結(jié)果分析

下面考慮材料參數(shù)在空間上沿徑向呈連續(xù)函數(shù)變化的功能梯度圓筒,為便于分析,做以下無量綱變換:

(R,A,B)=(r,a,b)/rm,E*=E/E0,α*=α/α0,

Θa,Θb)=(θ,θa,θb)/T0,

(∑r,∑θ)=(σr,σθ)/(α0T0E0),Ur=ur/α0T0rm,

Ha=rmβa/λa,Hb=rmβb/λb

式中:E0、α0、λ0、κ0為材料參考值;λa、λb分別表示圓筒內(nèi)、外壁的導(dǎo)熱系數(shù);rm為參考半徑,rm=(a+b)/2。

這里考慮一內(nèi)、外半徑分別為a=0.9 m和b=1.1 m的功能梯度圓筒。參考半徑rm=1 m,材料各參數(shù)的參考值為E0=225 GPa,α0=4.8×10-6K-1,λ0=5.9 W(mk)-1,κ0=2.8×10-6m2/s,ν=0.3。假定圓筒的初始溫度為參考溫度T0,內(nèi)、外表面突然置于溫度為Θa=1和Θb=0的環(huán)境中,且環(huán)境溫度Θa和Θb保持不變,圓筒的內(nèi)、外表面?zhèn)鳠嵯禂?shù)分別為Ha=30和Hb=55。設(shè)圓筒的材料參數(shù)分布服從以下函數(shù):

E*=em1(R-A)·f(T*),α*=em2(R-A)·f(T*)

λ*=em3(R-A)·f(T*),κ*=em4(R-A)·f(T*)

(15)

式中：取m1=2.0,m2=0.3,m3=3.0,m4=2.0；f(T*)為無量綱溫度函數(shù)。

Rishin等[10]在研究了金屬的彈性模量與溫度之間的關(guān)系后指出,材料的彈性模量隨溫度的增加而減小。為便于分析,在一般情況下,設(shè)溫度函數(shù)

f(T*)=1-ηT*

(16)

為得到足夠收斂的解,在計(jì)算中將圓筒沿厚度方向分為100層,時(shí)間步長取Δτ=0.0001。

圖1給出當(dāng)時(shí)間τ=0.002、0.005和達(dá)到穩(wěn)態(tài)后不同初始溫度相關(guān)參數(shù)ξ0下溫度Θ沿徑向R的分布。從圖1可以看出,在靠近圓筒內(nèi)壁附近區(qū)域,隨著ξ0的增加,同一時(shí)刻和位置的圓筒溫度會(huì)升高,而在其他區(qū)域,隨著ξ0的增加,同一時(shí)刻和位置的圓筒溫度會(huì)降低。但總的來說,在該算例中,溫度相關(guān)參數(shù)對(duì)圓筒溫度場(chǎng)的影響很小。

圖1 功能梯度圓筒溫度沿徑向分布

圖2給出當(dāng)時(shí)間τ=0.002、0.005和達(dá)到穩(wěn)態(tài)后不同初始溫度相關(guān)參數(shù)ξ0下徑向位移Ur沿徑向R的分布。從圖2可以看出,溫度相關(guān)參數(shù)對(duì)圓筒徑向位移有較大的影響。初始溫度相關(guān)參數(shù)ξ0越大,同一時(shí)刻圓筒徑向位移越小,隨著時(shí)間的增加,圓筒內(nèi)部溫度升高,溫度相關(guān)參數(shù)ξ也會(huì)增加,材料參數(shù)的數(shù)值會(huì)隨ξ的增加而減小,進(jìn)而使得圓筒徑向位移相對(duì)減小幅度更大。從圖中還可以看出,同一時(shí)刻在不同的ξ0下,圓筒內(nèi)部各點(diǎn)的徑向位移變化幅度相近。

圖2 功能梯度圓筒徑向位移沿徑向分布

圖3給出當(dāng)時(shí)間τ=0.002、0.005和達(dá)到穩(wěn)態(tài)后不同初始溫度相關(guān)參數(shù)ξ0下徑向應(yīng)力Σr沿徑向R的分布。從圖3可以看出,溫度相關(guān)參數(shù)對(duì)圓筒徑向應(yīng)力的影響很大。初始溫度相關(guān)參數(shù)ξ0越大,在同一時(shí)刻圓筒徑向應(yīng)力的絕對(duì)值越小,且不同ξ0之間的徑向應(yīng)力差異十分明顯。由于和徑向位移相同的原因,隨著時(shí)間的增加,徑向應(yīng)力絕對(duì)值相對(duì)減小的幅度越大。圖中還可以看出,同一時(shí)刻在不同ξ0下,靠近圓筒中部附近位置的徑向應(yīng)力變化幅度更大。

圖3 功能梯度圓筒徑向應(yīng)力沿徑向分布

為了分析本文對(duì)溫度函數(shù)的處理方法與文獻(xiàn)[7]中處理方法之間的誤差,以對(duì)圓筒徑向應(yīng)力分布的分析為例,圖4給出了當(dāng)ξ0=1.05時(shí),不同處理方法下圓筒徑向應(yīng)力分布之間的差別對(duì)比。由于文獻(xiàn)[7]中對(duì)溫度函數(shù)式(16)進(jìn)行線性化處理時(shí),用不變的初始溫度T0*代替了變化的溫度場(chǎng)T*,也就是說,溫度相關(guān)參數(shù)始終保持為ξ0不變。而從圖1可以看出,圓筒內(nèi)部的溫度升高幅度較大,則ξ會(huì)發(fā)生較大的改變,因此不能忽略溫度變化所帶來的影響。從圖4可以看出,文獻(xiàn)[7]的處理方法與本文的處理方法計(jì)算出的結(jié)果誤差較大,由于本文在計(jì)算過程中考慮了ξ會(huì)隨著溫度的增加而增加,因此得到的徑向應(yīng)力絕對(duì)值要更小(其他物理量同理)。隨著時(shí)間的增加,圓筒內(nèi)部溫度升高,兩種方法得到的徑向應(yīng)力之間差異也就更加明顯。

圖4 ξ0=1.05時(shí)對(duì)溫度函數(shù)不同的處理方法得到的徑向應(yīng)力分布差異對(duì)比

4 結(jié) 論

本文基于傅里葉定律、能量守恒方程、熱彈性方程和平衡方程,在空間上使用狀態(tài)空間技術(shù),時(shí)域內(nèi)應(yīng)用有限差分格式,建立了圓筒熱彈性問題在軸對(duì)稱平面應(yīng)變條件下的狀態(tài)空間理論,得到了材料參數(shù)隨溫度以及沿徑向任意梯度變化的圓筒熱彈性問題的狀態(tài)空間解。算例表明,溫度相關(guān)參數(shù)的增加會(huì)使材料參數(shù)的數(shù)值減小,從而使圓筒各物理量有不同程度的降低。相較于文獻(xiàn)[7]的線性化方法,本文通過迭代的方式在求解該類非線性問題時(shí)可以使計(jì)算精度得到提高,過程更簡便,計(jì)算效率更高。此外,本文方法能夠滿足軸對(duì)稱平面應(yīng)變問題下的多種邊界條件,還能退而求解均勻材料的相關(guān)問題。