高藝航，段靜波，雷勇軍

(1. 國(guó)防科技大學(xué) 空天科學(xué)學(xué)院， 湖南 長(zhǎng)沙 410073； 2. 中國(guó)運(yùn)載火箭技術(shù)研究院 北京宇航系統(tǒng)工程研究所， 北京 100076；3. 石家莊鐵道大學(xué) 工程力學(xué)系， 河北 石家莊 050043)

隨著航天器朝著高速、強(qiáng)突防、低成本、短周期方向發(fā)展，對(duì)結(jié)構(gòu)構(gòu)件的材料與成型提出了更高的要求。熱塑性樹脂基復(fù)合材料由于其耐高溫、高強(qiáng)度、高韌性、可重復(fù)使用等諸多優(yōu)異性能，成為可重復(fù)使用航天器首選的先進(jìn)復(fù)合材料。然而，熱塑性復(fù)合材料基體韌性大，受熱塑性顯著，在航天器極端復(fù)雜的氣動(dòng)熱環(huán)境中，熱塑性復(fù)合材料制成的壁板類結(jié)構(gòu)的氣動(dòng)彈性力學(xué)行力及機(jī)理十分值得研究。

關(guān)于壁板熱顫振問(wèn)題的研究最早可以追溯到 20 世紀(jì)50年代，Houbolt[1]最先開始研究溫度場(chǎng)均勻分布的二維壁板顫振邊界及其屈曲失穩(wěn)特性。隨后，Dowell[2]針對(duì)金屬壁板幾何非線性顫振問(wèn)題進(jìn)行了大量研究。Kouchakzadeh等[3]采用經(jīng)典板理論和von Karman非線性位移應(yīng)變關(guān)系進(jìn)行結(jié)構(gòu)建模，研究復(fù)合材料層合板在超音速氣流中的非線性氣動(dòng)彈性問(wèn)題。如果飛行馬赫數(shù)較高，在強(qiáng)烈的氣動(dòng)加熱效應(yīng)作用下，壁板將面臨極端嚴(yán)酷的熱環(huán)境。Zhou等[4]發(fā)現(xiàn)受熱壁板在氣流中的運(yùn)動(dòng)形式極其復(fù)雜，具有從低速氣流中的衰減振動(dòng)、熱屈曲振動(dòng)形式到高速氣流中各種類型的振動(dòng)。楊智春等[5]采用分步分析方法研究了復(fù)合材料鋪層方式對(duì)壁板熱顫振特性的影響。李凱倫等[6]對(duì)高超聲速環(huán)境中功能梯度薄板的熱氣動(dòng)彈性問(wèn)題進(jìn)行了研究，發(fā)現(xiàn)氣動(dòng)加熱效應(yīng)能夠使薄板發(fā)生熱屈曲或者提前進(jìn)入振動(dòng)狀態(tài)。Yang等[7-8]針對(duì)高超聲速飛行器一體化防熱結(jié)構(gòu)，建立了泡沫填充復(fù)合材料波紋夾芯板的高階分層氣動(dòng)彈性模型，討論了幾何參數(shù)和材料性能對(duì)顫振臨界動(dòng)壓的影響。Li等[9]研究了三角形柵格芯夾層板在超音速氣流中的顫振和屈曲，并采用位移反饋的方法設(shè)計(jì)了主動(dòng)控制器，為夾層結(jié)構(gòu)在飛行器設(shè)計(jì)中的應(yīng)用提供了理論依據(jù)。

關(guān)于熱塑性復(fù)合材料，?zen等[10]分別采用實(shí)驗(yàn)和仿真分析方法研究了熱塑性蜂窩夾芯板低能量沖擊響應(yīng)特性。Chen等[11]設(shè)計(jì)和制備了多層熱塑性復(fù)合材料(ThermoPlastic Composite, TPC)波紋夾芯板，采用平板壓縮試驗(yàn)研究了芯材結(jié)構(gòu)的失效機(jī)理。目前，公開的文獻(xiàn)資料大多是關(guān)于熱塑性復(fù)合材料的制造工藝、本構(gòu)關(guān)系、彈塑性變形、損傷測(cè)試方面的研究[12]，而關(guān)于熱塑性復(fù)合材料壁板的顫振問(wèn)題報(bào)道還比較少[13-15]?；诖?，本文以熱塑性復(fù)合材料層合板為研究對(duì)象，研究其顫振特性，分析關(guān)鍵參數(shù)對(duì)壁板熱顫振特性的影響規(guī)律。

1 超音速流場(chǎng)中復(fù)合壁板顫振有限元方程

1.1 結(jié)構(gòu)模型

如圖1所示，矩型熱塑性復(fù)合材料壁板的長(zhǎng)度和寬度分別為L(zhǎng)和H，厚度為h。由于復(fù)合材料壁板屬于厚板，因而，本文采用Mindlin厚板理論，壁板內(nèi)任意一點(diǎn)的位移場(chǎng)可寫為：

(1)

圖1 熱塑性復(fù)合材料壁板示意圖Fig.1 Sketch of thermoplastic composite panel

其中，u0、υ0、w0分別為中面上的點(diǎn)沿x軸、y軸和z軸方向的位移，θx和θy分別為中面法向轉(zhuǎn)角。

考慮高速氣流作用下壁板可能產(chǎn)生的幾何非線性變形，采用Von-Karman大變形理論，給出熱塑性壁板的應(yīng)變-位移關(guān)系如下：

(2)

其中：ε0包括壁板中面面內(nèi)位移產(chǎn)生的應(yīng)變分量和壁板大變形時(shí)撓度引起的非線性應(yīng)變分量；κ為彎曲時(shí)壁板的曲率向量；γ為橫向剪切應(yīng)變向量；z為壁板厚度方向的坐標(biāo)。

在熱塑性壁板應(yīng)力-應(yīng)變關(guān)系方面，考慮溫度對(duì)熱塑性壁板熱應(yīng)力的影響以及材料力學(xué)性能隨溫度的變化，暫不引入應(yīng)力應(yīng)變的塑性本構(gòu)關(guān)系，但在分析中討論屈服極限對(duì)壁板顫振應(yīng)力的影響。對(duì)于第k鋪層，并忽略法向正應(yīng)力，則壁板應(yīng)力-應(yīng)變關(guān)系為：

(3)

根據(jù)復(fù)合材料層合理論，n層復(fù)合變剛度壁板本構(gòu)關(guān)系可寫為：

(4)

其中，N、M、Fs分別為復(fù)合變剛度壁板的膜力、彎矩、橫向剪力等內(nèi)力，矩陣A、B、D、As見文獻(xiàn)[16]，矩陣NT、MT的表達(dá)式為：

(5)

1.2 氣動(dòng)力模型

在Ma>1.6的高速氣流中，壁板顫振可采用一階準(zhǔn)定常活塞理論計(jì)算氣動(dòng)力。 當(dāng)壁板表面來(lái)流沿x方向時(shí)，氣動(dòng)力表達(dá)式[17]為：

(6)

1.3 非線性顫振有限元方程

根據(jù)Hamilton原理，建立復(fù)合壁板的運(yùn)動(dòng)微分方程，即

(7)

其中，δT和δU分別為虛動(dòng)能和虛應(yīng)變能，δW為氣動(dòng)力和結(jié)構(gòu)阻尼所做的虛功。

基于復(fù)合壁板本構(gòu)關(guān)系、幾何關(guān)系以及氣動(dòng)力模型，分別寫出壁板體積域V內(nèi)的虛應(yīng)變能和虛動(dòng)能，以及其表面S上的外力虛功。 其中，虛應(yīng)變能包括壁板振動(dòng)產(chǎn)生的應(yīng)變能δUM和熱載荷產(chǎn)生的應(yīng)變能δUΔT兩部分，具體如下：

(8)

運(yùn)用有限元方法求解壁板顫振，采用四邊形四節(jié)點(diǎn)板單元對(duì)壁板進(jìn)行網(wǎng)格離散，總體位移列陣可記為如下形式：

(9)

其中

(10)

將總體位移列陣代入式(8)，通過(guò)變分則可獲得壁板的有限元顫振方程：

(11)

(12)

1.4 復(fù)合材料屈服準(zhǔn)則

基于Hill正交各向異性屈服準(zhǔn)則，Chen和Sun提出了適用于各向異性纖維增強(qiáng)復(fù)合材料的三維塑性勢(shì)函數(shù)。Weeks和Sun進(jìn)一步假設(shè)復(fù)合材料在纖維方向1上是線彈性的，在纖維方向2～3 面內(nèi)是橫觀各向同性的，從而將三維塑性勢(shì)函數(shù)簡(jiǎn)化為[14]：

(13)

本文研究的熱塑性C/PPS材料為緞紋機(jī)織復(fù)合材料，材料在纖維方向1、2上是線彈性的，且不考慮纖維方向3(厚度方向)的正應(yīng)力σ33，因而，上述塑性勢(shì)函數(shù)可最終退化為：

(14)

2 算例驗(yàn)證與分析

2.1 正確性驗(yàn)證

根據(jù)文獻(xiàn)[18]中T300/5208型環(huán)氧樹脂復(fù)合材料壁板算例，壁板幾何尺寸為0.3 m×0.12 m，厚度為1 mm，具體材料參數(shù)取值及顫振動(dòng)壓無(wú)量綱化與文獻(xiàn)[18]一致，采用本文程序進(jìn)行了固有頻率及顫振速度計(jì)算，結(jié)果見表1和表2。從表1中可以看出，本文計(jì)算的固有頻率與文獻(xiàn)結(jié)果吻合較好。采用有限元法分析時(shí)，通常需要通過(guò)細(xì)密網(wǎng)格才能保證分析結(jié)果的精度。因此，需要進(jìn)行網(wǎng)格收斂性驗(yàn)證。

表1 環(huán)氧樹脂復(fù)合材料壁板固有頻率

從表2中可以看到，隨著網(wǎng)格劃分?jǐn)?shù)量增多，結(jié)果逐漸收斂。當(dāng)網(wǎng)格劃分?jǐn)?shù)為40×16時(shí)，本文顫振速度計(jì)算結(jié)果與文獻(xiàn)吻合得很好。因此，綜合考慮線性壁板顫振的計(jì)算精度和效率，本文后續(xù)計(jì)算均采用40×16的網(wǎng)格數(shù)量。

表2 環(huán)氧樹脂復(fù)合材料壁板無(wú)量綱顫振動(dòng)壓

2.2 頻域線性顫振分析

以四邊簡(jiǎn)支的C/PPS熱塑性復(fù)合材料矩形壁板為對(duì)象，采用有限元法求解研究其熱固有特性和頻域熱顫振特性。熱塑性復(fù)合材料密度為0.9～1.6 g/cm3，本文計(jì)算中取1.25 g/cm3，泊松比取0.08，熱膨脹系數(shù)取1.6×10-6℃，其隨溫度變化的材料力學(xué)性能參數(shù)見表3[15]。壁板幾何尺寸為0.4 m×0.2 m，厚度為3 mm，鋪層方式為[0]7，來(lái)流方向平行于x軸。

表3 熱塑性材料力學(xué)性能參數(shù)

由于溫度載荷不僅產(chǎn)生熱應(yīng)力，而且溫度還改變材料的力學(xué)性能，這兩方面均對(duì)壁板顫振產(chǎn)生影響。熱塑性材料受溫度變化影響較大，需要同時(shí)考慮兩方面的影響。本文采用V-g法求解熱塑性復(fù)合材料壁板頻域顫振特性，并分析材料力學(xué)性能隨溫度變化對(duì)熱塑性復(fù)合材料壁板熱顫振特性的影響。

圖2分別給出了60 ℃和90 ℃下壁板前兩階模態(tài)阻尼隨無(wú)量綱動(dòng)壓的變化規(guī)律。壁板無(wú)量綱顫振動(dòng)壓定義見文獻(xiàn)[19]。從圖2中可以看到，兩種溫度下壁板熱顫振均由1階模態(tài)和2階模態(tài)耦合引起。在60 ℃下，由阻尼曲線在來(lái)流動(dòng)壓為211.0時(shí)過(guò)零點(diǎn)，這說(shuō)明壁板的無(wú)量綱顫振動(dòng)壓為211.0。同樣，在90 ℃下，由阻尼曲線過(guò)零點(diǎn)可得知，壁板無(wú)量綱顫振動(dòng)壓為149.5。

(a) 60 ℃下壁板的模態(tài)耦合(a) Modal coupling of the panel at 60 ℃

圖3給出了23 ℃、60 ℃、90 ℃、120 ℃四個(gè)溫度下考慮與不考慮材料性能隨溫度變化兩種情況下熱塑性壁板的熱顫振特性?？傮w上看，溫度的升高顯著降低了壁板的顫振動(dòng)壓。在常溫23 ℃下，壁板顫振無(wú)量綱動(dòng)壓最大?？紤]熱塑性材料溫變后得到的壁板無(wú)量綱動(dòng)壓低于僅考慮熱應(yīng)力時(shí)的顫振動(dòng)壓，從圖3中可以看到，在60 ℃、90 ℃、120 ℃下，壁板無(wú)量綱動(dòng)壓分別降低了5.1、17.0、34.8。由此可得出，考慮材料性能隨溫度變化后，得出的壁板顫振動(dòng)壓更低，而且溫度越高，材料性能隨溫度變化對(duì)壁板顫振特性影響越明顯。

圖3 不同溫度下熱塑性壁板熱顫振動(dòng)壓Fig.3 Thermal flutter dynamic pressure of thermoplastic panels at different temperatures

2.3 時(shí)域非線性顫振分析

以四邊簡(jiǎn)支的C/PPS熱塑性復(fù)合材料壁板為對(duì)象，采用有限元法從時(shí)域求解壁板非線性熱顫振特性。

圖4～6給出了60 ℃時(shí)熱塑性復(fù)合材料壁板特征點(diǎn)(0.5a，0.5b)在不同來(lái)流速度下的顫振時(shí)間歷程和相平面圖。從圖4可以看到，當(dāng)來(lái)流動(dòng)壓λ=188.6時(shí)，壁板受到初始擾動(dòng)后，振動(dòng)是收斂的，這表明動(dòng)壓未達(dá)到顫振臨界動(dòng)壓。當(dāng)來(lái)流動(dòng)壓λ=213.7時(shí)，壁板受到初始擾動(dòng)后，進(jìn)入極限環(huán)振動(dòng)，振幅數(shù)量級(jí)為1.7×10-3mm，此時(shí)動(dòng)壓超過(guò)顫振臨界動(dòng)壓，壁板做周期振動(dòng)，發(fā)生顫振，如圖5所示。當(dāng)動(dòng)壓進(jìn)一步增大，取來(lái)流動(dòng)壓λ=234.7時(shí)，壁板受到初始擾動(dòng)后，仍處于顫振極限環(huán)振動(dòng)，但是壁板做周期振動(dòng)的振幅增大為4.3×10-1mm，如圖6所示。

(a) 板中點(diǎn)位置撓度-時(shí)間曲線(a) Time history plot of deflection at the centre of the plate

(a) 板中點(diǎn)位置撓度時(shí)間曲線(a) Time history plot of deflection at the centre of the plate

(a) 板中點(diǎn)位置撓度時(shí)間曲線(a) Time history plot of deflection at the centre of the plate

圖7給出了熱塑性復(fù)合材料壁板在23 ℃、60 ℃、90 ℃、120 ℃下壁板極限環(huán)振幅隨無(wú)量綱動(dòng)壓的變化情況。從圖7中可以看出，不同溫度下壁板的顫振極限環(huán)振幅隨動(dòng)壓的增大而增大，增大趨勢(shì)基本一致。但隨著溫度從23 ℃逐漸增加到120 ℃，壁板進(jìn)入顫振的動(dòng)壓越來(lái)越低。此外，從圖7中還可以看出考慮與不考慮材料性能隨溫度變化對(duì)壁板極限環(huán)振幅的影響。在常溫23 ℃下，兩者沒有差異。在60 ℃下，考慮材料力學(xué)性能隨溫度變化時(shí)，壁板進(jìn)入顫振的動(dòng)壓低于不考慮材料力學(xué)性能隨溫度變化的情況。與此同時(shí)，相同動(dòng)壓下，考慮材料力學(xué)性能隨溫度變化時(shí)的極限環(huán)振幅要高于不考慮材料力學(xué)性能隨溫度變化的情況。在90 ℃、120 ℃下，規(guī)律相似于60 ℃的情形，只是壁板進(jìn)入顫振的動(dòng)壓差值隨著溫度從60 ℃到120 ℃逐漸增大，相同動(dòng)壓下，極限環(huán)振幅差異也越來(lái)越大。

圖7 不同溫度下壁板極限環(huán)振幅隨無(wú)量綱動(dòng)壓的變化Fig.7 Change of limit cycle amplitude of the panel along with dimensionless dynamic pressure at different temperatures

2.4 顫振應(yīng)力分析

文獻(xiàn)[15]給出的不同溫度下剪切試驗(yàn)數(shù)據(jù)表明，C/PPS熱塑性材料剪切強(qiáng)度都隨著溫度的增加而出現(xiàn)明顯的降低，在120 ℃時(shí)剪切強(qiáng)度比常溫時(shí)降低了63.64%。對(duì)于熱塑性壁板，在顫振發(fā)生過(guò)程中，如果應(yīng)力響應(yīng)達(dá)到材料塑性屈服強(qiáng)度將引起不可恢復(fù)的永久變形，這是熱塑性壁板熱顫振設(shè)計(jì)中需要考慮的問(wèn)題。下面仍以四邊簡(jiǎn)支的C/PPS熱塑性復(fù)合壁板為例，研究壁板顫振亞臨界過(guò)程中的應(yīng)力變化情況。

圖8給出了熱塑性復(fù)合材料壁板分別在23 ℃、60 ℃、90 ℃、120 ℃下，壁板極限環(huán)振蕩下等效應(yīng)力幅值隨動(dòng)壓變化的情況。從圖8中可以看到，在同一溫度下，隨著動(dòng)壓增大，壁板的應(yīng)力水平在增大。例如在23 ℃時(shí)，隨著壁板顫振動(dòng)壓由280.5增加到322.1，顫振等效應(yīng)力由0.399 MPa增加到1.27 MPa；在60 ℃時(shí)，隨著壁板顫振動(dòng)壓由220.2增加到216.7，顫振等效應(yīng)力由0.462 MPa增加到0.840 MPa。比較不同溫度下壁板的顫振等效應(yīng)力，可以看出，隨著壁板溫度的升高，在溫度應(yīng)力和材料溫變的共同作用下，壁板顫振等效應(yīng)力是逐漸降低。原因在于，溫度產(chǎn)生的結(jié)構(gòu)溫度應(yīng)力是壓應(yīng)力，而壁板顫振產(chǎn)生的應(yīng)力是反向的拉應(yīng)力，兩者呈現(xiàn)相互抵消的效應(yīng)，且溫度越高，抵消效應(yīng)越明顯。此外，從圖8中可以發(fā)現(xiàn)，壁板極限環(huán)振蕩下等效應(yīng)力幅值總體較小，均沒有達(dá)到各溫度下C/PPS熱塑性復(fù)合材料的屈服極限(見表3)，壁板沒有出現(xiàn)屈服區(qū)域。這主要是由于壁板顫振過(guò)程中產(chǎn)生的主要是雙向拉或壓應(yīng)力狀態(tài)，切應(yīng)力較小，而等效應(yīng)力則是反映切應(yīng)力水平，因而壁板極限環(huán)振蕩下等效應(yīng)力幅值總體較小。

圖8 不同溫度下壁板極限環(huán)振蕩等效應(yīng)力變化Fig.8 Change of limit cycle oscillation equivalent stress of the panel at different temperatures

3 結(jié)論

本文建立了超音速流場(chǎng)中熱塑性復(fù)合材料壁板熱顫振的有限元模型，與相應(yīng)文獻(xiàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證了本文模型及算法的正確性。進(jìn)而采用V-g法、Newmark法分別從頻域、時(shí)域求解復(fù)合變剛度壁板顫振特性，得出的主要結(jié)論如下：

1)熱塑性壁板熱顫振由1階模態(tài)(縱向一彎)和2階模態(tài)(縱向二彎)耦合產(chǎn)生，隨著溫度的升高，壁板的顫振動(dòng)壓顯著降低，而且考慮熱塑性材料溫變特性后得到的壁板顫振動(dòng)壓要低于僅考慮熱應(yīng)力時(shí)的顫振動(dòng)壓。

2)不同溫度下，當(dāng)考慮材料力學(xué)性能隨溫度變化時(shí)，壁板進(jìn)入顫振的動(dòng)壓低于不考慮材料力學(xué)性能隨溫度變化的情況。與此同時(shí)，相同來(lái)流動(dòng)壓下，考慮材料力學(xué)性能隨溫度變化時(shí)的極限環(huán)振幅要高于不考慮材料力學(xué)性能隨溫度變化的情況。

3)隨著壁板溫度的升高，壁板顫振等效應(yīng)力是逐漸降低的，而且壁板極限環(huán)振蕩下等效應(yīng)力水平總體較低，均沒有達(dá)到各溫度下復(fù)合材料的屈服極限。