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

吳雨桐

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

0 引 言

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

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

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

1 基本方程

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

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

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

(1)

(2)

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

在不計體力時,圓筒平衡方程為:

(3)

熱彈性應力和位移關系為

(4)

(5)

式中:σr和σφ分別表示圓筒徑向和環向應力。

初始條件設為:

θ(r,0)=0

(6)

邊界條件設為:

σ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分別代表圓筒在與外界進行對流換熱時的傳熱系數。

2 狀態方程的建立與求解

由式(4)可得:

(8)

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

(9)

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

(10)

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

(11)

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

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

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

(12)

對式(12)積分可得到狀態向量

[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)可以求出任意時刻在任意位置處的狀態變量。

3 算例和結果分析

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

(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分別表示圓筒內、外壁的導熱系數;rm為參考半徑,rm=(a+b)/2。

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

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*)為無量綱溫度函數。

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

f(T*)=1-ηT*

(16)

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

圖1給出當時間τ=0.002、0.005和達到穩態后不同初始溫度相關參數ξ0下溫度Θ沿徑向R的分布。從圖1可以看出,在靠近圓筒內壁附近區域,隨著ξ0的增加,同一時刻和位置的圓筒溫度會升高,而在其他區域,隨著ξ0的增加,同一時刻和位置的圓筒溫度會降低。但總的來說,在該算例中,溫度相關參數對圓筒溫度場的影響很小。

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

圖2給出當時間τ=0.002、0.005和達到穩態后不同初始溫度相關參數ξ0下徑向位移Ur沿徑向R的分布。從圖2可以看出,溫度相關參數對圓筒徑向位移有較大的影響。初始溫度相關參數ξ0越大,同一時刻圓筒徑向位移越小,隨著時間的增加,圓筒內部溫度升高,溫度相關參數ξ也會增加,材料參數的數值會隨ξ的增加而減小,進而使得圓筒徑向位移相對減小幅度更大。從圖中還可以看出,同一時刻在不同的ξ0下,圓筒內部各點的徑向位移變化幅度相近。

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

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

圖3 功能梯度圓筒徑向應力沿徑向分布

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

圖4 ξ0=1.05時對溫度函數不同的處理方法得到的徑向應力分布差異對比

4 結 論

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