沙云東,朱 林,欒孝馳,張國治,馮飛飛
(沈陽航空航天大學 遼寧省航空推進系統先進測試技術重點實驗室,沈陽 110136)
高超聲速飛行器在飛行過程中承受著嚴酷的航空動力加熱。為了維持高溫環境下結構的強度,高超飛行器的結構設計要求不同于低馬赫數飛行器。依據工作的溫度范圍,飛行器結構可被稱之為熱結構。由高溫合金設計而成的熱結構能夠在1 000至1 500高溫范圍內工作。若熱結構由C/C復合材料制造,則工作溫度可高達3 000。典型的高超飛行器熱結構部件包括:由單片鈦合金或金屬基復合材料制造而成的加筋板;由超級合金制造而成的蜂窩夾芯板以及C/C復合材料所制造而成的升降舵補助翼(或者襟翼)[1]。承受外部加熱的薄金屬板在其厚度方向存在較小的溫度梯度。然而,即便是在均勻的外部加熱條件下,由于支撐結構充當散熱片的角色,壁板表面也會存在空間溫度梯度分布。因此,典型的結構壁板在其中心承受
著較高的溫度,而在與支撐結構相連處溫度較低。這些溫度梯度結合壁板的初始缺陷,并通過壓縮薄膜應力的作用從而引發結構發生熱屈曲現象[2]。Ko[3-5]廣泛地研究了熱結構壁板(無散熱片作用)在均勻熱載荷條件下的熱屈曲問題。Thornton等[2]用有限單元法描述薄板由空間溫度梯度作用而引發的屈曲現象。并研究了穩態與非穩態溫度分布下薄板的彈性熱屈曲特性。Javaheri等[6]研究了不同溫度分布下薄板的熱屈曲的特性,并利用PDE/Galerkin法推導出了薄板在一維線性溫度梯度熱載荷作用下的臨界熱屈曲溫度差。不僅是受熱載荷的影響,高超飛行器在飛行過程中需要巨大的推力,從而引起巨大的強噪聲載荷,其數量級可達到150 dB到170 dB[7-8]。熱聲載荷聯合作用會使得飛行器壁板造成嚴重的疲勞破壞。Mei等[8-13]對帶有均勻穩態熱載荷分布的薄壁結構在隨機聲載荷作用下的動態響應進行了詳細的分析,表明薄壁板在熱聲載荷作用下會產生復雜的非線性動態響應,其中包括圍繞初始平衡位置的非線性隨機振動、圍繞熱屈曲后兩個平衡位置的非線性跳變運動以及圍繞熱后屈曲某一平衡位置的非線性隨機振動。但很少資料論述帶有溫度梯度熱載荷的薄板在隨機噪聲作用下的動態響應。因此,本文選用帶有一維線性溫度梯度熱載荷分布形式的薄板為例,利用有限單元法計算出了定常聲壓級下帶有溫度梯度熱載荷的四邊簡支矩形板在熱屈曲前后應力動態響應,并對應力響應進行了詳細地統計分析(概率密度 PDF、功率譜密度 PSD、有效值RMS)。
結構由于溫度變化所引起的熱變形受到約束時就會產生熱應力。分析溫度變化引起的應力、應變,則需要在廣義胡克定律中增加溫度項。由于下文所研究的結構為薄壁結構,忽略厚度方向溫度梯度的影響。因此,可將彈性力學中的應力-應變關系表示為(假設參考溫度T ref=0℃)

其中,α為熱膨脹系數,T(x,y)為薄板溫度分布,

根據直法線假設,距薄板中面距離為z的平面上任意一點的應變可用橫向位移w與中面應變分量表示如下

其中,中面應變分量有如下表示

將式(2)作相應微分處理,可以轉化為利用中面應變分量與橫向位移表示的變形協調方程

由式(1)的三個應力分量沿薄板厚度方向積分,薄板內力與彎矩有如下形式

將中面應變帶入其中,可得到面內力與彎矩的矩陣表達形式

其中,κ表示曲率,NT與MT分別表示熱內力與熱彎矩。而熱內力與熱彎矩可用如下方程表示

薄膜矩陣A與彎曲剛度矩陣D表示如下

根據M,N及橫向位移w可以給出薄板的平衡方程

引入應力函數ψ,

式(7)滿足力的平衡方程。將式(7)帶入式(6)可得到薄板的大撓度方程如下

假設四邊簡支服從不可動邊界條件,即四邊既無位移又無剪切力,可有如下方程

同時,四邊簡支矩形平板板邊界上的彎矩為0,即
根據式(8),結合邊界條件,可導出四邊簡支矩形平板的模態頻率計算公式[14]

將式(7)帶入式(5),導出利用應力函數所表示的中面應變分量,再帶入到變形協調式(3),最后可推出應力函數與位移共同表示的變形協調方程

式(8)與式(11)共同構成薄板大撓度控制方程。
若T為均勻穩態溫度熱載荷,則根據文獻[14],通過變形協調方程可導出四邊簡支板的臨界屈曲溫度
若僅考慮矩形板平面溫度梯度,設左右兩邊溫度分別維持在T1、T2,其物理參數不隨溫度而變化,并且無內熱源,根據穩態傳熱過程,可將傳熱學[15]中的三維導熱微分方程轉化為一維問題,轉化后的一維導熱微分方程形式如下

邊界條件為 T|x=0,T|x=a;對式(13)兩次微分可求出矩形平板的線性溫度分布如下

式中溫度差ΔT=T2-T1.結合文獻[6],可導出四邊簡支矩形板在該溫度梯度分布下的臨界熱屈曲溫度差計算公式

式中Tcr為四邊簡支矩形板在均勻穩態熱載荷下的第一階臨界熱屈曲溫度,由式(12)給出;此時,當 ΔT=Tcru時,四邊簡支矩形平板發生熱屈曲現象。
本文主要采用有限元數值計算方法進行計算求解,因此,根據Hamilton虛功原理,給出薄板熱聲載荷下的有限元運動方程如下[16]

式中M為質量矩陣;C為阻尼矩陣;K為線性剛度矩陣;KT是由于溫度變化引起的熱應力剛度矩陣;K1,K2為運動方程的第一階與第二階非線性剛度矩陣,且與位移w有關;而f與fT分別表示隨機聲載荷以及由溫度引起的等效熱載荷。
在給定邊界條件下,式(16)忽略慣性項與隨機載荷項,可得到非線性熱屈曲方程如下

利用Newton-Raphson迭代法對式(17)進行迭代,可以確定指定溫度分布下的靜態熱屈曲撓度{w}s。一旦獲得{w}s,可作為初始條件,通過添加慣性項,忽略非線性項,得到關于熱屈曲平衡位置{w}s的線性運動方程如下

其中,{w}t表示動態位移,結構總位移響應{w}=t;切線剛度矩陣 Ktan(ws)可在上述 Newton-Raphson迭代法計算中得到[17]。
根據式(18),熱屈曲薄板的模態頻率ω及振型n可通過下列方程進行求解

本文選取四邊簡支矩形鈦合金板作為研究對象,幾何尺寸與材料屬性如下表1所示。針對溫度梯度與聲載荷聯合作用,采用有限元數值計算方法對其進行了仿真計算,得到了四邊簡支矩形鈦合金板在帶有溫度梯度的熱載荷與聲聯合作用下的非線性應力動態響應。

表1 幾何尺寸與材料屬性Tab.1 The dimensions and material properties
根據式(10)與有限元法分別計算出四邊簡支矩形鈦合金板在常溫下模態頻率的解析解與數值解,如表2所示。利用式(12)以及有限元數值計算方法,計算出了四邊簡支矩形鈦合金板在均勻穩態熱載荷下第一階臨界熱屈曲溫度的解析解與數值解,如表3所示。解析解與數值解經過對比,極為近似,從而可驗證數值解的正確性。

表2 簡支板模態頻率(無熱預應力)Tab.2 Themodal frequencies of simply-supported p late(Hz)(No therm al prestress)

表3 簡支板臨界熱屈曲溫度(均溫)Tab.3The critical buckling temperature of simply-supported pate(Uniform tem perature)
上述式(14)給出了矩形板在x方向的溫度梯度分布,根據式(15)可以求出簡支板在該溫度梯度熱載荷分布下發生熱屈曲時的第一階臨界熱屈曲溫度差Tcru。結合表3中臨界熱屈曲溫度的數值解,利用式(15)可以擬合出該溫度梯度熱載荷分布下簡支板第一階臨界熱屈曲溫度差隨T1的變化規律曲線,如圖1所示,從圖中發現,隨著T1不斷增加,臨界熱屈曲溫度差線性遞減,直到T1=42.4℃時,臨界熱屈曲溫度差等于0,即簡支板在均勻溫度42.4℃下,矩形板發生熱屈曲失穩,此時,可認為均溫導致的失穩為溫度梯度作用下的一種特殊情況。
本文以圖1中所標注點的位置為例,所研究的溫度梯度是以溫度差ΔT=30℃為基準,且T1線性變化??梢园l現,當T1=27.4℃時,矩形板發生熱屈曲現象。為了后文表達方便,可令 S=T1/27.4,即保持 ΔT=30℃,當S=1時,矩形板發生熱屈曲失穩現象。本文計算所選取的S范圍從0到4,間隔0.2。

圖1 臨界熱屈曲溫度差(溫度梯度)Fig.1 The critical buckling temperature difference(temperature gradient)
計算得到了不同S(ΔT=30℃)下四邊簡支板的熱預應力模態頻率,如下表4所示。結合表2與表4,熱屈曲前,溫度的升高使模態頻率逐漸下降,屈曲后階段,溫度致使簡支板基頻呈現上升趨勢。以S=1(熱屈曲時)為例,給出四邊簡支矩形鈦合金板的第一階Von Mises應力振型云圖,如下圖2所示。從圖中發現,溫度梯度導致簡支板一彎振動最大Von Mises應力有向溫度較高一側偏移的趨勢。
本文所考慮的隨機聲激勵載荷是以總聲壓級160 dB的限帶高斯白噪聲,如圖3所示。其截止頻率為1 500 Hz,且均勻分布在矩形平板上,根據上述模態頻率計算,可以覆蓋多階模態頻率;載荷信號時長為1 s時的聲壓時間歷程如下圖3(a)所示,其聲壓的概率密度服從正態分布;相應的聲壓功率譜如下圖3(b)所示;

圖2 簡支板Von Mises應力分布(S=1,ΔT=30℃)Fig.2 Von Mises stress distribution for simple supported plate(S=1,ΔT=30℃)

圖3 限帶高斯白噪聲(總聲壓級160 dB)Fig.3 Band limited white Gauss noise(the overall sound pressure level 160 dB)
本節利用有限元數值計算方法計算了簡支板在不同熱屈曲系數S(溫度梯度不變,ΔT=30℃)、總聲壓級為160 dB下的應力動態響應。
由于所研究模型為對稱模型,且溫度梯度分布沿板x向中心線對稱,因此,以簡支板x向中心線為基準,提取薄板中心線及下側部分結點進行分析,所要計算的結點位置及編號如下圖4所示。

圖4 簡支板結點編號及位置Fig.4 Node number and position of simple supported plate

表4 簡支板模態頻率 (ΔT=30℃)Tab.4 Themodal frequencies of sim p ly-supported p late(Hz)(ΔT=30℃)
上圖2指出,一彎振動最大Von Mises應力雖然向高溫一側便宜,但偏移量較小,因此,為了說明帶有溫度梯度熱載荷的簡支板在隨機聲載荷作用下的動態特性,以編號為57的中心結點為例,給出其x向應力響應時間歷程,如圖5所示,相應的概率密度如圖6所示。從中可以清晰地揭示薄板在熱聲載荷作用下的三種運動形式,即屈曲前,薄板圍繞初始平衡位置做隨機振動;屈曲后薄板的跳變現象以及圍繞熱后屈曲某一平衡位置做非線性隨機振動。
屈曲前,隨著溫度的增大,應力響應時間歷程的對稱中心逐漸遠離x軸,說明熱應力在上升,應力均值絕對值增大,而應力幅值也在增加,如圖5(a),圖5(b),圖5(c)所示。從下述概率密度分布圖6(a),圖6(b)中可以看出,S=0且ΔT=0℃時,應力均值為0,S=0.6且ΔT=30℃時,應力均值增加到約40 MPa,最大應力幅值大約從20 MPa增加到60 MPa。
圖6(a)指出,無熱載荷影響時,其概率密度基本服從正態分布。而圖6(b),圖6(c)表明由于溫度的影響,其概率密度不服從正態分布。
在熱后屈曲跳變過程中,薄板中點x向應力時間歷程處于拉伸與壓縮兩種狀態,如圖5(d)所示,拉伸時,應力幅值大于壓縮應力幅值,而拉伸應力均值卻小于壓縮應力均值,這一點在圖6(d)概率密度中體現的更加顯著。同時,在下述概率密度分布圖6(d)中了解到,由于跳變的產生,應力概率密度出現雙峰值狀態,且不服從正態分布。

圖5 簡支板中點應力響應 (SPL=160 dB)Fig.5 x component of stress for themidpoint of simple supported plate
圖 5(e),圖 5(f)指出薄板在 S=1.6,2(ΔT=30℃)時完全處于熱后屈曲拉伸階段的應力時間歷程,隨著溫度的增加,由于熱屈曲撓度的作用,拉伸應力均值逐漸增大,但應力振幅卻逐漸減小。在下述概率分布圖6(e)與圖 6(f)中,應力均值約從 60 MPa(S=1.6)增加到大約 128 MPa(S=2),而最大應力幅值約從60 MPa(S=1.6)減小到20 MPa(S=2),并且,在熱后屈曲狀態,應力概率密度恢復為正態分布。

圖6 簡支板中點應力概率密度(SPL=160 dB)Fig.6 Probability density for themidpoint of simple supported plate

圖7 簡支板中點應力功率譜密度(ΔT=30℃)Fig.7 Stress power spectrum density for themidpoint of simple supported plate
對于總聲壓級為160 dB,熱屈曲系數 S=0,0.6,1,1.2,1.6(ΔT=30℃)時薄板中心結點 x向應力的時間歷程通過自相關函數計算,然后經過傅里葉變換,可轉化為頻域下的應力功率譜密度,如下圖7所示。由于結構振動基頻最為主要,因此,僅考慮基頻變化情況??砂l現,定常聲壓級下,屈曲前,隨著溫度的升高,薄板基頻下降;屈曲后,溫度升高,導致薄板基頻開始上升;這是由于在熱屈曲前,薄板隨著溫度的增加出現軟化過程,導致薄板剛度降低;而熱屈曲后,溫度的增加致使薄板逐漸硬化,導致其基率升高;以上基頻變化情況與表2與表4中熱屈曲前后基頻變化規律基本一致。這種穩定、失穩、再穩定的過程,在文獻[9,11]中均溫熱載荷對矩形平板模態頻率的影響論述研究中也有所體現。
本文通過應力有效值(RMS)統計分析,得到了結點57在總聲壓級160 dB、不同S(ΔT=30℃)下的Von Mises應力有效值(RMS),如圖下8所示。
熱屈曲前(S=0)直到進入熱后屈曲的頻繁跳變階段(S=1.2),簡支板應力有效值隨著溫度的增加而線性遞增;從頻繁跳變(S=1.2)開始直到跳變結束階段(S=1.5),由于跳變逐漸減弱致使簡支板的應力有效值隨著溫度的增加呈現下降趨勢;而后,應力有效值逐漸上升,此時,簡支板圍繞熱后屈曲某一平衡位置振動,由于隨機性,結點57的應力將會出現拉伸或壓縮兩種狀態,這里采用拉伸應力對其有效值進行計算。而由于熱屈曲前與熱后屈曲跳變階段,拉伸與壓縮應力概率基本相同,計算應力有效值時,可不必對其進行區分。值得說明的一點,簡支板在受到熱聲載荷的聯合作用時,存在彎曲拉伸應力、彎曲壓縮應力以及壓縮熱應力。根據不同的溫度以及簡支板結點位置的不同,三種應力在數值上大小各不相同,因而三種應力相互疊加的結果會出現單純意義上的拉伸與壓縮狀態。
考慮到溫度梯度(ΔT=30℃)對簡支板的影響,結合圖4,下面給出了簡支板不同位置結點Von Mises應力有效值比較示意圖,如下圖9所示。圖9(b),圖9(c)中結點應力在熱后屈曲階會呈現出拉伸或壓縮狀態,因此,簡支板在熱后屈曲階段均采用拉伸應力進行有效值計算。而圖9(a),9(d)所對應的結點應力在熱屈曲前后均為壓縮狀態,即為壓縮應力有效值。圖9(a)與圖9(d)中簡支板三邊中點及角點應力有效值曲線表明,溫度較高一側,其應力有效值相對較大。圖9(b)所給出的三點應力有效值經比較發現,中點應力有效值最大,其余對稱兩結點溫度較高一側應力有效值較大。圖9(c)中指出y=0.3×1/8 m處的結點應力有效值,由于溫度梯度的存在,在熱屈曲前后數值的大小出現交替現象,熱屈曲前到跳變結束階段,溫度較高一側應力有效值較大,熱后屈曲圍繞上凸平衡位置振動時,溫度較低一側應力有效值較大。
圖9指出,簡支板上所求得應力有效值的結點中,除了角點1和2以外,其余結點應力有效值隨溫度的變化趨勢基本相同,同時受到頻繁跳變、間歇跳變直至跳變結束這段過程的影響,應力有效值在這段過程隨溫度的增加均出現下降現象。然而,簡支板角點1和2的應力有效值在熱屈曲前后隨著溫度的增加卻線性遞增。圖9(a)中三個邊界中點均為壓縮應力,雖然彎曲壓縮應力與壓縮熱應力的總和大于彎曲拉伸應力,但在跳變過程中,彎曲拉伸應力卻在其中占有重要位置。彎曲拉伸應力、彎曲壓縮應力以及壓縮熱應力之間相互作用,導致應力有效值出現下降趨勢。角點1、2所處的位置其彎曲拉伸應力遠小于壓縮熱應力,因而壓縮熱應力致使角點1、2的應力有效值隨著溫度的增加線性遞增。
本文利用有限元數值計算方法對帶有溫度梯度熱載荷的四邊簡支矩形鈦合金板進行熱聲激振計算,得到結論如下:
(1)根據簡支板在線性熱梯度作用下的第一階臨界熱屈曲溫度差計算公式可知,隨著簡支板低溫一側邊界線上溫度的增加,臨界熱屈曲溫度差線性降低。
(2)簡支板在熱屈曲前,隨著溫度的遞增,應力均值增大,應力幅值增加;跳變過程中,應力處于拉伸與壓縮兩種狀態,拉伸應力幅值大于壓縮應力幅值,但拉伸應力均值卻小于壓縮時應力均值;隨著溫度繼續增加,拉伸應力均值逐漸增大,但應力振幅卻逐漸減小。
(3)在熱屈曲前直到熱后屈曲頻繁跳變階段,除角點以外,簡支板其余結點應力有效值隨著溫度的增加而線性遞增;從熱后屈曲頻繁跳變到跳變結束階段,簡支板應力有效值呈下降趨勢;溫度繼續提高,拉伸應力有效值逐漸上升。而角點應力有效值隨著溫度的增加線性遞增。
[1]Ko W L.Thermal buckling analysis of rectangular panels subjected to humped temperature profile heating[R].NASA/TP-2004-212041,2004.
[2]Thornton E A,Kolenski J D,Marino R P.Finite element study of plate buckling induced by spatial temperature gradients[C] //AIAA/ASME/ASCE/AHS/ASC 34th Structures,Structural Dynamics,and Materials Conference.1993,1:2313-2326.
[3]Ko W L.Mechanical and thermal buckling analysis of rectangular sandwich panels under different edge conditions[J]. NASA STI/Recon Technical Report N, 1994,94:33704.
[4]Ko W L.Predictions of thermal buckling strengths of hypersonic aircraft sandwich panels using minimum potential energy and finite elementmethods[M].National Aeronautics and Space Administration,Office of Management,Scientific and Technical Information Program,1995.
[5]Ko W L.Thermostructural behavior of a hypersonic aircraft sandwich panel subjected to heating on one side[M].National Aeronautics and Space Administration,Office of Management,Scientific and Technical Information Program,1997.
[6]Javaheri R,Eslami M R.Thermal buckling of functionally graded plates[J].AIAA Journal,2002,40:162-169.
[7]Dhainaut J M,Guo X,Mei C,et al.Nonlinear random response of panels in an elevated thermal-acoustic environment[J].Journal of Aircraft,2003,40(4):683-691.
[8]Mei C,Dhainaut JM,Duan B,et al.Nonlinear random response of composite panels in an elevated thermal environment[R].Old Dominion Univ Norfolk,VA,2000.
[9]Sha Y D,Li J Y,Gao Z J.Dynamic response of pre/post buckled thin-walled structure under thermo-acoustic loading[J].Applied Mechanics and Materials,2011;80-81:536-541.
[10]Przekop A,Rizzi SA.Dynamic snap-through of thin-walled structures by a reduced-order method[J].AIAA Journal,2007,45(10):2510-2519.
[11]Sha Y D,Gao Z J,Xu F,et al.Influence of thermal loading on the dynamic response of thin-walled structure under thermo-acoustic loading[J].Advanced Engineering Forum,2011,2-3:876-881.
[12]沙云東,魏靜,高志軍,等.熱聲激勵下金屬薄壁結構的隨機疲勞壽命估算[J].振動與沖擊,2013,32(10):162-166.SHA Yun-dong,WEI Jing,GAO Zhi-jun,et al.Random fatigue life prediction ofmetallic thin-walled structures under thermo-acoustic excitation[J]. Journal of Vibration and Shock,2013,32(10):162-166.
[13]沙云東,魏靜,高志軍,等.熱聲載荷作用下薄壁結構的非線性響應特性[J].航空學報,2013,34(6):1336-1346.SHA Yun-dong,WEI Jing,GAO Zhi-jun,et al.Nonlinear responses characteristics of thin-walled structures under thermo-acoustic loadings[J]. Acta Aeronautica et Astronautica Sinica,2013,34(6):1336-1346.
[14]Sha Y D,Gao Z J,Xu F.Influences of thermal loads on nonlinear response of thin-walled structures in thermo-acoustic environment[J].Applied Mechanics and Materials,2011,105-107:220-226.
[15]楊世銘,陶文銓.傳熱學(第4版)[M].北京:高等教育出版社,2006.8.
[16]Guo X Y,Przekop A,Mei C,Nonlinear random response of shallow shells at elevated temperatures using finite element modal method[C]. Structural Dynamics& Materials Conference,2004,19-22.
[17]Guo X Y, Lee Y, Mei C, et al. Thermal buckling suppression of supersonic vehicle surface panels using shape memory alloy[J].Journal of Aircraft,2004,41(6):1498-1504.