考慮氣動加熱瞬態效應的復合材料壁板超音速LCO特性研究

劉雅婷, 段靜波,2, 徐步青,2, 高藝航

(1. 石家莊鐵道大學 工程力學系,石家莊 050043; 2. 石家莊鐵道大學 河北省工程力學基礎學科研究中心,石家莊 050043; 3. 國防科技大學 空天科學學院,長沙 410073; 4. 北京航天系統工程研究所,北京 100076)

超音速飛行器在高速飛行過程中會產生大量的熱量。嚴重的氣動加熱會顯著改變飛行器蒙皮結構溫度場,導致惡劣的熱環境,從而可能導致飛行器力學性能、氣動彈性大變形和失穩邊界發生變化,進一步降低飛行器的安全性[1]。為了解決這些問題,許多研究者對飛行器在熱環境下的氣動彈性行為進行了研究。通常,文獻研究中主要采用穩態溫度場,忽略了氣動加熱產生的瞬態溫度場對力學性能和結構穩定性的影響。自吳志剛等[2]提出分層求解的方法之后,研究者們廣泛采用這種單向耦合的方法來研究瞬態氣動加熱對氣動彈性行為的影響[3]。Tran等[4]研究了F-16二維機翼在氣動加熱和氣動力作用下的溫度分布和時域響應。Mcnamara等[5]將氣流與結構之間的傳熱結合起來,討論了氣動加熱對結構穩定性的影響。Miller等[6]提出了一種流固熱松耦合模型,并證明了壓力、位移插值和熱流密度的估算保持了壁板運動參數和溫度的二階精度。Huo等[7]采用基于活塞理論和薄激波層的時適應松耦合方法研究了C/SiC復合材料板的氣熱彈性行為。Chen等[8]采用時間自適應高超聲速熱-流-固松耦合模型用于熱氣動彈性響應分析,重點研究了低展弦比機翼的自振和熱模態變化。Khalafi等[9]利用三階活塞理論和廣義微分求積法分析了功能梯度板的熱氣動彈性特性。Ye等[10]研究了高超聲速機翼熱氣動彈性變形對機翼扭轉角和下表面壓力的變化規律的影響。Kamali等[11]采用有限元求解器計算熱解和結構解模擬了熱氣動彈性耦合,并將計算結果與試驗數據進行了對比。

隨著高超聲速飛行器的出現,氣動和結構的耦合作用愈加顯著,氣動彈性變形對瞬態氣動加熱的影響不容忽視[12]。Thornton等[13]采用有限元法分析了不銹鋼薄板的流-熱-固耦合作用。結果表明,氣動彈性變形會影響加熱速率和溫度分布。Li等[14]采用熱氣動彈性雙向耦合方案,研究了功能梯度曲板不同拱高對顫振特性的影響。Wan等[15]采用雙向耦合的靜態熱氣動彈性分析方法,得到了高超聲速機翼的熱流密度、氣動力、熱場和變形。研究發現,由溫度引起的變形對氣動力分布有很大影響。Ye等[16]開發了松耦合熱氣動彈性分析系統,采用單向耦合和雙向耦合計算方法研究了熱氣動彈性變形對高超聲速進氣道性能和流場的影響。Huang等[17]利用模態降階法方法和緊耦合模式進行了熱氣動彈性的計算。

基于此,本文建立了考慮氣動加熱瞬態溫度效應的復合材料壁板熱氣動彈性模型。在模型驗證后,重點討論了氣動加熱瞬態溫度場中斜激波、氣動壓力、傳熱系數、初始干擾力等關鍵因素對復合材料壁板非線性熱氣動彈性LCO響應的影響規律,為超音速飛行器復合材料壁板的工程設計提供技術支持。

1 熱氣動彈性建模

以超音速氣流中的復合材料層合板為研究對象,如圖1所示。壁板沿長度x方向為a,沿長度y方向為b,壁板的厚度為h。建模過程中,考慮超音速氣流經過楔形體產生的斜激波、氣動加熱產生的壁板內瞬態溫度場,以及壁板結構的大撓度變形等。

圖1 復合材料壁板示意圖Fig.1 Schematic diagram of the composite panel

1.1 激波參數和氣動力

超音速氣流經過楔形體會被壓縮,在楔形體前緣產生斜激波,因此,復合材料壁板邊界層氣體參數會發生變化。壁板表面激波后氣體參數與來流氣體參數關系如下[18]

(1a)

(1b)

(1c)

(1d)

式中:β為斜激波角;ρ∞,V∞,Ma∞和T∞分別為來流密度、來流速度、來流馬赫數和溫度;ρe,Ve和Te為激波后氣流參數。

當來流馬赫數Ma∞>1.6時,壁板表面沿x方向的氣動力分布根據一階活塞理論獲得[19]

(2)

(3)

式中,D110為鋪層角為0°時壁板彎曲剛度矩陣的第一個元素。

1.2 壁板的本構關系

復合材料壁板在高速氣流中可能會產生非線性變形,因此,基于Von karman大變形假設,壁板的位移-應變關系可以表示為

(5)

式中:u0,v0和w0分別為中性層在x,y和z方向的位移;ε0為中面的總面內應變向量;κ為板的曲率;γ為剪切應變;θx和θy分別為中性面繞x軸和y軸的轉角。

同時,考慮瞬態非均勻溫度場對復合材料壁板熱應力和力學性能的影響,結構的應力-應變關系定義為

(6)

根據層合板的層合效應,n層對稱復合材料壁板的內力可以由下式進行計算

(7a)

Fs=Asγ

(7b)

(7c)

式中:N,M和Fs為壁板的膜力、彎扭矩和剪力;A,D,As分別為剛度矩陣、彎曲剛度矩陣和剪切剛度矩陣;NΔT和MΔT為溫度變化產生的膜力和彎矩。

1.3 氣動加熱模型

考慮壁板上表面對外界的氣體輻射效應,根據牛頓冷卻公式,可以得到氣動加熱產生的復合材料壁板表面的熱流

(8)

式中:α為熱傳遞系數;Tw為壁溫;Taw為絕熱壁溫可以表示為

(9)

α=ρ*VeSt*cp

(10)

式中:Ve為邊界層氣流速度;cp為定壓比熱容;ρ*為參考溫度下的邊界層氣流密度;St*為參考溫度下的斯坦頓數。參考溫度可以由下式進行計算

(11)

1.4 壁板熱傳導模型

由于氣動加熱產生的熱流會在氣流方向發生變化,因此,考慮壁板熱傳導沿厚度方向(z方向)和氣流方向(x方向)的溫度差異,采用二維熱傳導模型分析,其控制方程為

(12)

式中:ρ和c分別為材料密度和比熱容;kx和kz為材料在x和z方向的熱傳遞系數。

根據壁板上表面與空氣的對流換熱邊界條件可得

(13)

式中:k為壁板與空氣之間的對流換熱系數;q為氣動加熱的熱流密度,可以由式(8)求得。

將邊界條件代入控制方程,采用有限差分法計算瞬態溫度場。沿x方向和z方向進行網格劃分,z方向網格數與壁板層數一致。具體差分格式可以表示為

(14)

式中,i和j分別為x和z方向的節點編號。

1.5 氣動彈性模型

復合材料壁板結構變形產生的虛應變能δUM、溫度變化產生的虛應變能δUT、壁板的虛動能δT,氣動力和阻尼力產生的外力虛功δW,具體如下

(15a)

(15c)

δW=?AΔpδwdA

(15d)

式中,ρ為材料密度。

然后,利用四邊形四節點單元對復合材料層合板進行離散。單元的位移向量可表示為

u=Nuδe

(16)

式中,u,δe和Nu分別為各單元的未知位移向量、節點位移向量和形函數矩陣。具體表達式為

u=[u,v,w,θx,θy]

(17a)

δe=[u1,v1,w1,θx1,θy1…u4,v4,w4,θx4,θy4]

(17b)

(17c)

根據哈密頓變分原理

(18)

將式(16)代入式(18),通過變分可以得到復合材料壁板的熱氣動彈性分析的非線性運動微分方程

(19)

式中:M為質量矩陣;Ca為氣動阻尼矩陣;K0,Knon,Ka和KT分別為線性剛度矩陣、非線性剛度矩陣、氣動力矩陣和熱應變剛度矩陣;fΔT為溫度引起的壁板總體載荷列陣。對應矩陣的具體表達式見附錄A。

在對式(19)求解時,首先進行模態降階,縮減后的運動方程可以表示為

(20)

[F({q}i+1)]{q}i+1-{f}i+1=0

(21)

其中,

(22b)

最后,由于運動方程中具有非線性項,采用牛頓迭代法對式(21)進行求解。

2 氣動加熱和氣動彈性振動耦合計算

由于計算中考慮氣動加熱的瞬態效應,因此,氣動加熱和氣動彈性振動需在時間增量步內分別求解,再進行瞬態溫度響應和氣動彈性振動響應的數據交換,依此進行氣動加熱與氣動彈性振動耦合計算。具體而言,首先,根據t時刻溫度得到材料參數,進而計算t～t+Δt時刻壁板的氣動彈性響應。其次,根據壁板的變形量可以獲得壁板表面t+Δt時刻的氣動力分布。同時,可以計算t+Δt時刻與壁板變形有關的氣動參數。在氣動加熱分析中,利用更新后的氣動參數可以得到t+Δt時刻的溫度場。接著,用同樣的方法可以得到下一時間步壁板的氣動彈性響應,直至最后。計算流程圖如圖2所示。

圖2 熱氣動彈性耦合計算流程圖Fig.2 Flow chart of aerothermoelastic coupling calculation

3 算例驗證與結果分析

3.1 算例驗證

為了驗證現有算法和程序,首先采用本文的程序計算壁板的瞬態溫度場,并與Culler等研究進行對比,驗證壁板中心點溫度隨時間的變化。壁板幾何參數和材料物性參數均與Culler等研究中的一致,計算結果如圖3所示。可以看到,本文計算結果與文獻結果吻合很好。這說明本文程序關于壁板瞬態溫場計算的程序是正確的。

圖3 壁板中心點瞬態溫度場的變化Fig.3 Variation of transient temperature field at the center point of the composite panel

其次,通過穩態溫度場中四邊簡支矩形壁板的顫振驗證本文關于氣動彈性模型計算程序的正確性。復合材料壁板的幾何尺寸和材料參數與歐陽小穗等研究中的相同。采用歐陽小穗等研究中的兩種鋪層結構。將本程序得到的結果與歐陽小穗等的研究結果進行比較,如圖4所示。這說明本文關于壁板氣動彈性計算程序是正確的。

圖4 LCO振幅隨無量綱動壓的變化Fig.4 Variation of LCO amplitude with dimensionless dynamic pressure

3.2 氣動加熱瞬態溫度場中壁板時域響應特性

以四邊固定的矩形復合材料層合壁板為研究對象,研究其在考慮氣動加熱瞬態溫度效應下的顫振響應特性。壁板尺寸及材料參數,如表1所示。考慮壁板上表面氣動加熱以及壁板上表面對外界的氣體輻射效應,氣動加熱產生的熱流為q,壁板左右兩側為絕熱邊界條件,下表面為恒溫邊界條件,T=300 K,如圖5所示。

表1 壁板尺寸和材料參數Tab.1 Panel dimensions and material parameters

圖5 熱邊界條件示意圖Fig.5 Schematic diagram of thermal boundary conditions

為了便于比較,本文同時給出了瞬態溫度場和穩態溫度場中復合材料壁板的顫振時域響應和溫度響應,如圖6所示。由圖6(a)中可以看出,復合材料壁板在考慮氣動加熱瞬態溫度效應下的LCO幅值與穩態溫度下的計算結果基本一致。這說明通過穩態溫度場近似預測瞬態熱環境下壁板顫振特性是合理的。然而,從壁板整個時間響應歷程來看,受溫度場瞬態特性的影響,復合材料壁板的氣動彈性振動時域響應不同于穩態溫度場情形,主要體現在振動的早期階段。從圖6(a)可以看到,穩態溫度場中的復合材料壁板振動響應很快就會進入LCO狀態,時間過程相對較短。但是考慮了氣動加熱產生的溫度場瞬態效應后,復合材料壁板在進入LCO狀態前有一個相對較長的振動歷程,而且壁板在進入LCO狀態前其振動位移是先收斂后增大,見圖6(a)所示。主要原因是氣動加熱下壁板內溫度場存在一個溫度逐漸升高過程,如圖6(b)所示。值得注意的是,計算得到的瞬態溫度場最終將是隨時間持續微幅振蕩變化的,這主要是由壁板最終處于LCO狀態導致的。此外,圖6(c)和圖6(d)分別給出了復合材料壁板上一個周期內不同時刻的位移振型圖和溫度場分布圖。可以看到,壁板的周期變形影響其表面的氣動力分布,進而影響氣動加熱熱流,最終導致壁板溫度場在長度方向以圖6(d)所示形態微幅波動。

圖6 復合材料壁板氣動彈性振動時域響應與瞬態溫度響應(λ=166)Fig.6 Aeroelastic vibration time-domain and transient temperature responses of the composite panel(λ=166)

3.3 主要參數影響分析

進一步討論斜激波、傳熱系數、初始擾動力等因素對考慮氣動加熱瞬態溫度效應下壁板氣動彈性LCO響應和瞬態溫度響應的影響。

首先研究斜激波對復合材料壁板在瞬態溫度場中LCO響應的影響。如圖7所示給出了楔形體半錐角為10°時復合材料壁板的振動時域響應和瞬態溫度響應。從圖7(a)中可以發現,斜激波會使壁板的LCO幅值顯著增加,也就是說,斜激波的產生會大大降低壁板的臨界顫振動壓,降低壁板的顫振穩定性邊界。圖7(b)給出了斜激波作用下壁板的瞬態溫度場。從圖中可以看出,斜激波的產生也會使壁板的溫度場顯著增大,這也是導致壁板臨界顫振動壓變化的原因所在。而且,從圖8中不同動壓下瞬態溫度場響應可以看出,動壓越大,斜激波對壁板瞬態溫度場影響越明顯,進而對壁板的顫振影響也會越大。

圖7 激波前后復合材料壁板氣動彈性振動時域響應與瞬態溫度響應(λ=157)Fig.7 Aeroelastic vibration time-domain and transient temperature responses of the composite panel ahead of and behind shock wave(λ=157)

圖8 不同動壓下激波前后復合材料壁板瞬態溫度響應Fig.8 Transient temperature responses of the composite panel ahead of and behind shock waves at different dynamic pressures

圖9～圖11分別給出了復合材料壁板在長度方向和厚度方向的溫度場分布。從圖9中不難發現,壁板前緣受氣動加熱的影響最嚴重,溫度梯度較大,而沿著壁板的長度方向溫度梯度逐漸減緩。另外,瞬態溫度場會隨動壓的增大而增大,但相比于激波前,激波后動壓對瞬態溫度場的影響變得更為明顯。圖10所示動壓為λ=164時,復合材料壁板內不同位置的瞬態溫度場隨時間的變化情況,可以看出,壁板內瞬態溫度場隨著時間變化逐漸趨于等幅波動,但不同位置處的溫度場波動幅值不同。圖11進一步給出了壁板厚度方向的溫度場分布。可以看到,壁板沿厚度方向會存在非均勻的溫度梯度。而且,通過動壓λ=155和λ=164兩組溫度分布結果對比,可以得出,隨著動壓增大,壁板厚度方向的溫度變化明顯增大,這將會影響壁板的臨界顫振動壓。因此,考慮壁板厚度方向的溫度非均勻分布是十分必要的。

圖9 不同動壓下激波前后復合材料壁板上表面溫度場分布Fig.9 Temperature field distribution on the upper surface of the composite panel ahead of and behind shock waves under different dynamic pressures

圖10 激波后復合材料壁板不同位置的瞬態溫度響應(λ=164)Fig.10 Transient temperature responses of the composite panel at different positions behind shock wave(λ=164)

圖11 激波后復合材料壁板沿厚度方向的溫度場分布Fig.11 Temperature field distribution along thickness of the composite panel behind shock wave

圖12為傳熱系數對氣動加熱瞬態溫度場中復合材料壁板LCO響應的影響。其中,傳熱系數k分別取值為20 W/m/K,40 W/m/K,60 W/m/K三種情況。從圖中可以看出,隨著傳熱系數從60 W/m/K降低到20 W/m/K,壁板的LCO幅值顯著增大,這意味著壁板的臨界顫振動壓會明顯下降。原因在于傳熱系數較大時,氣動加熱會在復合材料壁板內產生變化相對較小的瞬態溫度場,如圖12(b)所示,進而影響了壁板的臨界顫振動壓。

圖12 不同傳熱系數下復合材料壁板氣動彈性振動時域響應與溫度響應(λ=180)Fig.12 Aeroelastic vibration time-domain and transient temperature responses of the composite panel at different heat transfer coefficients(λ=180)

圖13給出了不同初始擾動力作用下考慮氣動加熱瞬態溫度效應的復合材料壁板振動時域響應和瞬態溫度響應。從圖13(a)中可以看出,初始擾動力不會影響復合材料壁板的顫振臨界動壓,但是在時間歷程上會影響復合材料壁板初始階段的收斂振動,以及壁板進入LCO的時間。當初始擾動力較大時,壁板初始階段的收斂振動幅值變化顯著,且壁板進入LCO的時間較早。當初始擾動力較小時,壁板初始階段的收斂振動幅值變化不明顯,壁板進入LCO的時間相對較晚。從圖13(b)中可以看到,不同初始擾動力作用下,壁板的瞬態溫度場也會發生變化,這其中的原因主要是初始擾動影響了壁板的氣動彈性振動響應,進而影響了壁板的瞬態溫度場。

圖13 不同初始擾動力下復合材料壁板氣動彈性振動時域響應與瞬態溫度響應(λ=157)Fig.13 Aeroelastic vibration time-domain and transient temperature responses of the composite panel at different initial disturbance forces(λ=157)

4 結 論

本文建立了氣動加熱瞬態溫度效應的復合材料壁板氣動彈性模型,給出了壁板的氣動彈性時域響應和瞬態溫度響應特性,并在考慮氣動加熱瞬態溫度效應下,重點討論分析了斜激波、來流動壓、傳熱系數、初始干擾力等對復合材料壁板的氣動彈性LCO響應的影響。主要結論如下:

(1)考慮了氣動加熱產生的溫度場瞬態效應后,復合材料壁板的氣動彈性振動時域響應早期階段不同于穩態溫度場情形,壁板在進入LCO狀態前有一個相對較長的振動歷程,且壁板振動位移是先收斂后增大,計算得到的溫度場最終是隨時間持續微幅振蕩變化的。

(2)斜激波會使壁板的溫度場升高,并且相比于激波前,激波后動壓對瞬態溫度場的影響更為明顯,這進而導致復合材料壁板的LCO幅值顯著增加,大大降低壁板的臨界顫振動壓,降低壁板的顫振穩定性邊界。

(3)傳熱系數的變化對壁板的LCO振幅有顯著的影響,傳熱系數減小,會導致壁板內瞬態溫度場變化加劇,從而通過影響壁板剛度使壁板的LCO振動幅值增大,降低壁板的臨界顫振動壓。

(4)初始擾動力不會影響復合材料壁板的顫振臨界動壓,但是在時間歷程上會影響復合材料壁板初始階段的收斂振動,以及壁板進入LCO的時間,也會影響瞬態溫度場振動發生的時間。

附錄A

(A.1)

(A.2)

(A.3)

(A.4)

(A.5)

(A.6)

(A.7)

其中,

(A.8)

(A.9)

(A.12)