屈佑文，安效民，鄧 斌，閆 浩，周 悅

(西北工業(yè)大學 航天學院，西安 710072)

0 引 言

為了減輕結構重量，改善結構性能，各種新型復合材料被大量應用于高超聲速飛行器結構中。飛行器高超聲速飛行時，會出現(xiàn)氣動加熱，復合材料薄壁結構在氣動載荷和熱載荷作用下發(fā)生變形，導致流場的邊界發(fā)生改變，使得作用在結構上的載荷發(fā)生變化，進一步影響薄壁結構的響應，使其產(chǎn)生新的變形或者振蕩。這種結構、 氣動和熱之間的相互耦合作用，造成了壁板的熱氣動彈性響應問題。與傳統(tǒng)翼、 舵等升力面構型不同的是，薄壁結構在非定常載荷作用下會表現(xiàn)出非線性流-固耦合特征： 其位移響應通常與厚度為相同量級，面內(nèi)應力存在，彎曲和拉伸之間會發(fā)生耦合，在周圍受約束條件下，形變一般不會引發(fā)結構的迅速破壞，表現(xiàn)為幾何非線性。

國內(nèi)外一些學者在該領域有一些很有意義的研究，Nydick等基于Marguerre薄殼理論和Galerkin方法求解結構變形方程，分析了高超聲速氣流中三維正交各向異性薄板的極限環(huán)幅值和振蕩形態(tài)，指出壁板的氣動彈性響應受到初始曲率、 氣動熱、 激波運動等的影響。Shiau等研究了高超聲速氣流中層合板的非線性顫振響應，指出材料各向異性對壁板的極限環(huán)和混沌運動形態(tài)等有顯著影響。Gray結合Von-Karman板理論和三階活塞理論分析了復合材料薄壁板的非線性氣動彈性響應，分析了支承條件、 系統(tǒng)參數(shù)對顫振特性的影響。Kouchakzadeh等基于Galerkin方法計算了超聲速氣流中復合材料層合板的顫振特性，并討論了面內(nèi)應力、 纖維方向、 氣動阻尼等對非線性響應的影響。

梁路等研究了具有復合材料壁板結構的機翼氣動彈性問題，結果表明復合材料不同鋪層的單層厚度比對機翼的剛度影響很大。歐陽小穂等研究了高速氣流下變剛度復合材料壁板的非線性顫振響應，分析了支承條件和纖維方向等對顫振邊界的影響。林華剛研究了非均勻溫度場中復合材料圓柱殼的顫振特性，并進行了氣動熱-氣動彈性耦合的計算分析。苑凱華等分析了復合材料壁板的非線性顫振特性，結果表明高超聲速氣流中氣動力的非線性和熱載荷對壁板顫振的振動幅值的影響更為明顯。楊智春等利用有限元方法分析了不同鋪層方向和鋪層順序?qū)雍蠌秃喜牧媳诎孱澱裉匦缘挠绊?，發(fā)現(xiàn)隨溫度升高，發(fā)生頻率重合型顫振的耦合模態(tài)會發(fā)生演變，使得壁板有可能在較低的速壓下發(fā)生顫振。

對于大尺度高超聲速飛行器而言，在持續(xù)穩(wěn)定飛行狀態(tài)下，飛行器的表面存在近似均勻升溫區(qū)域。本文基于均勻溫度場針對復合材料薄壁結構在高超聲速氣流中的非線性氣動彈性響應展開分析，建立了考慮熱效應的復合材料壁板的非線性有限元求解模型，并將其應用于氣動-結構的二階松耦合方法中，用于分析考慮非線性非定常氣動載荷和溫度載荷共同作用下的非線性氣動彈性問題，揭示非線性耦合的基本關系。

1 復合材料薄壁結構的非線性結構動力學建模

1.1 復合材料薄壁殼元的切線剛度矩陣

薄壁結構中面變形可以看作是薄膜變形和彎曲變形相互耦合的結果，引入一階剪切變形假設。復合材料薄壁結構如圖1所示，其中(,,)為單元參考坐標系，(,)為鋪層材料參考坐標系。

圖1 復合材料薄壁結構示意圖

考慮熱效應，單元參考坐標系下的單層應力-應變關系可寫為

(1)

單元參考坐標系內(nèi)，包括考慮熱效應后的合力和力矩在內(nèi)的所有合力和合力矩，都可以寫成如下形式：

(2)

(3)

考慮單元參考坐標系和材料坐標系之間的轉(zhuǎn)換關系，可知：

(4)

由此，可以推導得到單元參考坐標系中復合材料壁板殼元的剛度矩陣為

(5)

式中：，和分別為薄膜應變矩陣、 彎曲應變矩陣和橫向剪切應變矩陣； 由虛功原理，熱效應可以轉(zhuǎn)換為溫度載荷，其在單元參考坐標系中的表達式可以寫為

(6)

利用高斯積分可以獲得單元坐標下的剛度矩陣和溫度載荷。考慮熱效應后，載荷和單元位移的關系如下：

+Δ=

(7)

式中：和分別是單元參考系下的位移和載荷。

基于Co-rotational formulation (CR)理論，總體坐標系下的復合材料薄壁結構的切線剛度矩陣為

=+m

(8)

式中：m為幾何剛度矩陣。

1.2 熱效應下結構非線性動力學方程

基于哈密爾頓方程和虛功原理，不考慮結構阻尼的情況下，推導得到增量形式的結構動力學運動方程如下：

(9)

2 非定常氣動載荷建模

無黏非定常流場的控制方程為Euler方程，其守恒積分形式如下：

(10)

(11)

3 氣動-結構二階松耦合方法的計算流程

本文基于二階松耦合方法(如圖2所示)，步驟如下：

(1) 利用時刻的結構信息，預測+12時刻的結構狀態(tài)：

(12)

將耦合邊界上的結構信息轉(zhuǎn)換到流場網(wǎng)格中。為了在邊界上滿足應力連續(xù)，流場物面壓強梯度考慮施加慣性力影響：

(13)

(2) 利用動網(wǎng)格算法更新流場網(wǎng)格，計算新的控制體體積。

(3) 求解式(11)，獲得+12時刻流場的壓強，積分得到物面的氣動載荷，將其轉(zhuǎn)換為結構等效載荷s, +12，預估+1時刻的載荷：

s, +1=2s, +12-s,

(14)

(4) 利用Newmark方法求解式(9)，得到+1時刻的結構狀態(tài)。

圖2 流-固耦合示意圖

4 算法驗證

4.1 復合材料壁板的有限元建模驗證

考慮一個四周固支的四層[/-/-/]復合材料層合板模型，如圖3所示。材料屬性為：=1，=1100，=500，=2 GPa，=154，=03,=08，=1 500 kg/m。

圖3 復合材料層合板模型

考慮兩種鋪層結構：=45°, 90°。結構有限元網(wǎng)格數(shù)量為45×45，本文計算所得的前6階自然頻率如表1所示。本文的結果與文獻中的結果相當吻合，說明了本文的復合材料有限元模型的可靠性。

表1 前6階自然頻率計算結果

4.2 熱結構算例驗證

以一個各向同性平板模型為例對熱結構模型進行驗證，其結構為：=1 m，=0.1 m，=1/100，=03，=2 700 kg/m，=10×10。彈性模量隨溫度的變化如圖4所示。兩端固支約束，有限元網(wǎng)格數(shù)量為41×5，初始溫度為=100 ℃。

圖4 彈性模量隨溫度的變化曲線

考慮溫度場的變化[100 ℃，900 ℃]，基于本文的熱效應有限元模型進行模態(tài)分析，并與Nastran的結果進行比較。圖5給出了在不同溫度下前6階熱模態(tài)頻率的計算結果對比?？梢钥闯?，本文的計算結果與Nastran的結果頗為接近，除了第4階模態(tài)(該模態(tài)為扭轉(zhuǎn)模態(tài))有稍許差異外，其他模態(tài)的頻率吻合較好。第4階模態(tài)最大誤差不超過9%，這種差異可能與兩種算法的非線性剛度矩陣計算的不同有關。

由圖可知，隨溫度的升高，各階頻率會隨之降低，更高階的頻率下降幅度更大。這種特性很容易造成各階模態(tài)運動之間的耦合，從而誘發(fā)如壁板顫振等不利的氣動彈性現(xiàn)象。前兩階頻率隨溫度出現(xiàn)了先降低后增高的趨勢。以第1階結構模態(tài)為例，其頻率隨溫度的升高先變小，在400 ℃時達到最小，在超過429.5 ℃時，1階頻率又上升，此時結構發(fā)生了熱屈曲，熱屈曲以后會產(chǎn)生較大的變形。屈曲對板頻率的影響，其原因較為復雜，與薄板屈曲構型、 溫度梯度造成的熱應力以及振動模態(tài)有關。本文計算中，基頻模態(tài)在屈曲后頻率隨溫度升高而增大。這也與Schneider等的研究結果相符。

圖5 前6階結構固有頻率隨溫度的變化對比

4.3 非定常氣動力計算驗證

以AGARD CT5標模為例進行驗證，翼型為NACA0012，馬赫數(shù)為0.755，翼型非定常運動規(guī)律為()=0016+251·sin(2×0081 4) 。圖6顯示了本文計算的升力系數(shù)和俯仰力矩系數(shù)與實驗值的比較。兩者吻合很好，說明了本文所述方法應用于無黏非定常流場計算時的可靠性。

圖6 AGARD CT5算例結果

4.4 三維壁板顫振分析驗證

考慮一個三維各向同性壁板，如圖7所示。其幾何和材料特性為：=1，=0002，=03，=()=01，顫振分析時選取的計算狀態(tài)為： 超聲速=12, 壁板的流場網(wǎng)格為H型，包括121×121×39節(jié)點，結構有限元模型由1 152個三角殼元構成。

圖7 三維壁板的幾何構型

圖8顯示了壁板極限環(huán)振蕩幅值和頻率，并與Dowell和Gordnier的結果進行了對比。顯然，與Gordnier的結果相符較好。Gordnier的求解器中同樣采用了基于CFD的非定常載荷計算，而Dowell幅值稍高，因為其選用了線性勢流理論來計算流場。這個算例也說明了本文耦合計算方法在處理氣動力非線性和幾何大變形結構非線性引起的氣動彈性響應上的有效性。

圖8 壁板極限環(huán)振蕩幅值和頻率隨動壓的變化

5 壁板的非線性熱氣動彈性響應分析

5.1 壁板模型與參數(shù)

壁板幾何和材料屬性參數(shù)為：==1.0 m，=0.005 m，=400，=06，=05，=025，=1 500，質(zhì)量比=()=02。壁板四個邊界均固支，壁板包含5個等厚度鋪層?？紤]兩種鋪層： 正交鋪層板[0°90°0°90°0°](cross-ply)和斜交鋪層板[45°-45°45°-45°45°](angle-ply)。計算工況為=5，給定溫度條件分別為Δ=0 ℃, 50 ℃, 100 ℃, 200 ℃，分析壁板在不同動壓下的非線性熱氣動彈性響應，以=075,=05作為參考點。

5.2 網(wǎng)格無關性和時間步長收斂性分析

壁板流場網(wǎng)格為H型，考核了三種網(wǎng)格密度，沿壁板順氣流方向、 沿展向以及法向的網(wǎng)格點分別設置為： 81×81×29(壁板表面41×41)、 121×121×39(壁板表面61×61)和181×181×49(壁板表面81×81)。選取三個不同的無量綱時間步長： Δ=0015, 003, 006。采用Richardson外插方法考核網(wǎng)格和時間步長的影響，圖9顯示了以歸一化無量綱振蕩幅值為考核的收斂情況。

圖9 不同網(wǎng)格和時間步長下的GCI收斂情況

綜合考慮計算精度和效率，本文后續(xù)計算的網(wǎng)格選取為121×121×39，耦合計算時間步長Δ=003。圖10顯示了壁板流場網(wǎng)格，壁板在下表面中心位置，下表面設置為無穿透條件，流入邊界上的變量由自由流超聲速條件確定。在流出邊界上，通過外推獲得變量，其他邊界變量由遠場黎曼不變邊界條件定義。圖11顯示了/=0.5剖面處壁板發(fā)生變形后的網(wǎng)格圖。壁板流-固耦合響應計算中，不同下的流場初始值取的是剛性壁板的定常流場。

圖10 壁板流場網(wǎng)格圖

圖11 y/b=0.5處壁板變形后的網(wǎng)格剖面

5.3 正交鋪層壁板在不同溫度下的響應分析

(1) Δ= 0 ℃下的響應分析

對于正交壁板，Δ=0 ℃時，發(fā)現(xiàn)=766時出現(xiàn)了顫振，其相平面圖如圖12(a)所示。隨著動壓的增大，后顫振的幅值和頻率也會隨之增加，在較大動壓下，相圖會出現(xiàn)粉刺現(xiàn)象，與高階模態(tài)介入有關, 如圖12(b)所示。圖13～14分別顯示了=1 500、 參考點振蕩達到=0°(正向極值)和=180°(負向極值)時的壁板變形云圖和流場壓強云圖?？梢钥闯?，此時的振蕩形式表現(xiàn)除了結構的2階模態(tài)之外，還有高階形態(tài)參與其中。從壓強云圖來看，壁板參考點正向極值時前緣有激波，后緣為膨脹波； 其負向極值時前緣有膨脹波，后緣有激波出現(xiàn)。

圖12 正交鋪層在ΔT=0 ℃時不同動壓時的相平面圖

圖13 正交鋪層在ΔT=0 ℃下λ=1 500時壁板位移圖

圖14 正交鋪層在ΔT=0 ℃下λ=1 500時壓強云圖

(2) Δ=50 ℃下的響應分析

當Δ=50 ℃時，在=377時出現(xiàn)了顫振，顫振動壓相比未施加溫度載荷的結果變小約50%，其相平面圖如圖15(a)所示。隨動壓增大，顫振幅值和頻率也會隨之增加, 如圖15(b)所示。

圖15 正交鋪層在ΔT=50 ℃時不同動壓下的相平面圖

(3) Δ=100 ℃下的響應分析

當Δ=100 ℃時，壁板在亞臨界下會出現(xiàn)靜態(tài)屈曲。在不同動壓條件下，該屈曲位置由正向朝負向發(fā)生了轉(zhuǎn)變, 如圖16所示。當動壓增大到=120時出現(xiàn)了低頻顫振，此時壁板無法保持在熱屈曲的平衡位置，受到擾動后出現(xiàn)了振蕩。 其響應相平面圖如圖17(a)所示，

圖16 正交鋪層在ΔT=100 ℃時屈曲響應

出現(xiàn)了雙環(huán)狀構型。隨動壓增大，逐漸演化為單環(huán)形，如圖17(b)所示。圖18顯示=1 000，=0°和=180°時的壁板變形云圖?？梢钥闯?，此時的振蕩形式表現(xiàn)為1階自然振型。

圖17 正交鋪層在ΔT=100 ℃時各動壓下的壁板相平面圖

圖18 正交鋪層在ΔT=100 ℃下λ=1 000時壁板位移圖

(4) Δ=200 ℃下的響應分析

當Δ=200 ℃時，其亞臨界下的形態(tài)也表現(xiàn)為靜態(tài)屈曲，不同動壓下該屈曲位置也出現(xiàn)了正向到負向的轉(zhuǎn)變，并且隨著動壓的增大，其收斂過程較長。其顫振出現(xiàn)在=249，從圖19中可以看出，顫振及后顫振時壁板的形態(tài)表現(xiàn)為類似混沌運動。圖20～21分別顯示了=1 000，=0°和=180°時的壁板變形云圖和流場壓強云圖?？梢钥闯?，此時的振蕩形式中含有高階自然振型，流場的壓縮波和膨脹波分布與未施加溫度載荷時相似。

圖19 正交鋪層在ΔT=200 ℃時不同動壓下的相平面圖

圖20 正交鋪層在ΔT=200 ℃下λ=1 000時壁板位移圖

圖21 正交鋪層在ΔT=200 ℃下λ=1 000時壓強云圖

5.4 斜交鋪層壁板在不同溫度下的響應分析

(1) Δ=0 ℃下的響應分析

對于斜交鋪層壁板，當Δ=0 ℃，=1 265時才會出現(xiàn)顫振，該臨界動壓遠遠超過了正交鋪層壁板的值(約1.7倍)，其響應如圖22(a)所示，=1 500時的相平面圖如圖22(b)所示。圖23顯示了=1 500，=0°和=180°時的壁板變形云圖?？梢钥闯?，此時的振蕩形式主要表現(xiàn)為2階自然振型，含有部分高階形態(tài)。

圖22 斜交鋪層在ΔT=0 ℃時不同動壓下的相平面圖

圖23 斜交鋪層在ΔT=0 ℃下λ=1 500時壁板位移圖

(2) Δ=50 ℃下的響應分析

當Δ=50 ℃,=689時出現(xiàn)了顫振。隨動壓增大，顫振幅值和頻率也會增加，如圖24所示。

(3) Δ=100 ℃下的響應分析

當Δ=100 ℃時，斜交壁板在未發(fā)生顫振時也出現(xiàn)了靜態(tài)屈曲。隨動壓的變化，屈曲位置也會由正向朝負向轉(zhuǎn)變。當動壓增大到=284時出現(xiàn)了低頻顫振，其響應相平面圖出現(xiàn)了雙環(huán)狀構型，如圖25(a)所示。隨動壓增大，這種構型逐漸演化為單環(huán)形, 如圖25(b)所示。

圖24 斜交鋪層在ΔT=50 ℃時各動壓下的相平面圖

圖25 斜交鋪層在ΔT=100 ℃時各動壓下的相平面圖

(4)Δ=200 ℃下的響應分析

當Δ=200 ℃時，亞臨界下靜態(tài)屈曲仍然存在，并且參考點屈曲位置也隨著動壓的變化而出現(xiàn)正向到負向的轉(zhuǎn)換。顫振發(fā)生在=386，顫振及后顫振響應表現(xiàn)出擬周期特點，如圖26所示。圖27顯示了=1 000,=0°和=180°時的壁板變形云圖，可以看出，此時的振蕩形式主要表現(xiàn)為1階自然振型。

圖26 斜交鋪層在ΔT=200 ℃時各動壓時的壁板響應相平面圖

圖27 斜交鋪層在ΔT=200 ℃下λ=1 000時壁板位移圖

5.5 不同溫度載荷下復合材料壁板的響應對比

圖28比較了兩個壁板在不同溫度下的振蕩幅值與頻率隨動壓的變化，可以看出，幅值和頻率均隨動壓的增大而增大。Δ=0 ℃時，相同動壓下，正交壁板的振蕩幅值與頻率均高于斜交壁板的值； Δ=50 ℃，=2 000時，正交壁板的幅值略小于斜交壁板的值(約91%)； Δ=100 ℃時，兩個壁板在小動壓下呈現(xiàn)低頻振蕩(<10 Hz)、 振蕩極值呈對稱性； Δ=200 ℃時，斜交壁板在小動壓時也表現(xiàn)為低頻特征，且在較大動壓(≥1 000)下，振蕩極值呈現(xiàn)出非對稱性。

圖28 兩種鋪層壁板的振蕩幅值和頻率隨動壓的變化

圖29顯示了Δ=100 ℃和Δ=200 ℃時亞臨界條件下參考點的屈曲位置。明顯觀察到兩個屈曲平衡點，隨著動壓的增大，壁板的位置從正向平衡點向負向平衡點轉(zhuǎn)換。相同動壓條件下，斜交壁板的屈曲位移更大。

圖30比較了正交鋪層和斜交鋪層在不同溫度下的振蕩幅值和頻率隨動壓的變化?？梢园l(fā)現(xiàn)，極限環(huán)振蕩的幅值隨著溫度的增加而增加，極值也呈對稱特點。對于斜交壁板，在較大溫度和動壓條件下，會出現(xiàn)非對稱的振蕩。隨著溫度的增大，對于正交壁板，溫度最大時，其振蕩頻率也最大，但斜交壁板對應最大溫度載荷下的振蕩頻率最小。

圖29 兩種鋪層壁板的屈曲位置隨動壓的變化

圖30 兩種鋪層壁板的振蕩幅值和頻率隨動壓的變化

圖31為正交鋪層和斜交鋪層在不同溫度下的顫振臨界動壓和頻率?？梢园l(fā)現(xiàn)，隨著溫度的增加，顫振臨界動壓先明顯減小，Δ=100 ℃時，相比常溫狀態(tài)，兩個壁板的減小量分別為84%和78%，而后隨溫度增加又增大； 相同溫度下，斜交壁板的顫振臨界動壓高于正交壁板。對于顫振頻率而言，正交壁板的顫振頻率隨溫度載荷的增加，先減小后增大，而斜交壁板的顫振頻率隨溫度載荷的增加而變小。

圖31 兩種鋪層壁板在不同溫度載荷下的顫振臨界動壓和頻率

6 結 論

本文基于氣動-結構二階耦合的求解思路，針對非線性非定常氣動載荷與非線性結構響應在熱環(huán)境下的非線性氣動彈性問題展開分析，建立了考慮熱效應的幾何非線性復合材料壁板求解模型，針對高超聲速下不同鋪層壁板的顫振特性對比研究，探討了鋪層方向、 溫度等對顫振特性的影響。結論如下：

(1) 不同鋪層的薄壁結構會呈現(xiàn)出不同的顫振特性。一般而言，斜交壁板的顫振動壓高于正交壁板(增量大于55%)，在相同動壓下，正交壁板后顫振時的頻率高于斜交壁板。在一定溫度下，會出現(xiàn)低頻顫振特性，與壁板屈曲有關。

(2) 隨著溫度的增大，顫振動壓會下降，但出現(xiàn)熱屈曲時，顫振動壓又會緩慢增大。

(3) 當溫度超過結構的臨界屈曲溫度時，壁板在亞臨界條件下會出現(xiàn)靜屈曲響應，而且存在多個屈曲平衡位置，該位置隨著不同動壓條件發(fā)生轉(zhuǎn)變。

(4) 本文分析中溫度為給定情況，流場求解考慮的Euler方程，后續(xù)研究中需要分析真實條件下氣動-熱-結構耦合下的響應情況。

(5) 本文僅針對高超聲速典型馬赫數(shù)狀態(tài)點=5開展不同均勻溫度場下的流-固耦合分析研究，所采用的方法對于其他高超聲速馬赫數(shù)狀態(tài)下的研究具有參考意義，后續(xù)研究還需增加高超聲速狀態(tài)，以完整刻畫其演變規(guī)律。