郭靜 王建民 張忠 李海波 高博

（北京強度環(huán)境研究所 可靠性與環(huán)境工程技術(shù)重點實驗室，北京100076）

0 引言

高馬赫數(shù)飛行的航天飛行器，氣動力、熱、結(jié)構(gòu)多物理場相互耦合，易引發(fā)熱氣動彈性穩(wěn)定性問題[1]。新型復合材料的廣泛應用，使得壁板結(jié)構(gòu)易出現(xiàn)極限環(huán)振蕩、混沌等復雜非線性熱氣動彈性響應形式，影響飛行安全，亟需開展壁板結(jié)構(gòu)非線性熱氣動彈性分析方法研究[2-3]。

大變形引起的幾何非線性和非定常氣動力的高效預示是非線性熱氣動彈性分析應考慮的主要問題。壁板顫振、超聲速顫振和熱顫振廣泛采用簡單、高效的活塞理論預示非定常氣動力[4-8]，活塞理論也在應用中不斷改進和發(fā)展[9]，先后提出了Van Dyke活塞理論[10]、當?shù)亓骰钊碚揫11]等。楊智春等[12]以超音速氣流中二維受熱壁板為對象，開展了非線性熱顫振分析，給出了其基本型和二次失穩(wěn)型兩種顫振模式。夏巍等[13]用一階活塞理論模擬壁板受到的氣動力，分析了超音速氣流中受熱壁板后屈曲狀態(tài)下的穩(wěn)定性。楊翊仁等[14-16]研究了二維、三維壁板在超音速氣流作用下各種復雜的熱彈性顫振動力學現(xiàn)象，給出了二維壁板發(fā)生動態(tài)Hopf分岔的條件[17]。夏巍[18]、Mei[19]推導了考慮熱效應的顫振有限元方程，研究了幾何非線性引起的壁板分岔與混沌響應。

本文充分利用Nastran軟件SOL400中的非線性有限元求解程序，采用DMAP語言對其進行二次開發(fā)，將溫度場分析和Van Dyke修正活塞理論非定常氣動力計算引入SOL400非線性瞬態(tài)響應直接求解模塊，編制非線性熱氣動彈性分析程序，并采用PCL語言在Patran軟件中開發(fā)嵌入式前后處理界面，將Van Dyke修正活塞理論非定常氣動力計算的高效性與SOL400非線性瞬態(tài)響應計算的高效性相結(jié)合，建立了工程實用的非線性熱氣動彈性分析工具。以復合材料壁板為研究對象，進行典型溫度工況下的非線性熱氣動彈性分析。結(jié)果表明，本文方法可對壁板結(jié)構(gòu)LCO，周期，非周期和混沌等復雜非線性熱氣動彈性響應進行有效預示。

1 Van Dyke修正活塞理論

根據(jù)經(jīng)典活塞理論，翼面上某點壓力與該點下洗速度之間滿足如下關(guān)系[20]

式中，p為翼面上某點的壓力，p∞、ρ∞、a∞分別為無窮遠來流壓力，密度和聲速，W為翼面下洗速度，γ為氣體絕熱比，空氣通常取 =1.4γ。當下洗速度小于音速時，活塞理論可根據(jù)麥克勞林級數(shù)保留項階次稱為一階、二階或三階活塞理論。

針對馬赫數(shù)較小的超聲速工況，采用Van Dyke二階理論修正后的壓力系數(shù)Cp可表示為

式中，M為馬赫數(shù)，a為當?shù)匾羲?，則活塞理論可轉(zhuǎn)化為

2 基于Patran/Nastran的非線性熱氣動彈性分析程序

矩陣形式的熱氣動彈性運動方程表達式為

式中，[M ]，[K]分別為質(zhì)量矩陣和剛度矩陣，A(u)為結(jié)點力列向量，由Van Dyke修正活塞理論求得。剛度矩陣[K]由小位移線性剛度矩陣[k0]，幾何剛度矩陣[kσ]和初位移矩陣[kL]疊加表示，[K ]= [k0]+ [kσ]+[kL]。剛度矩陣[K]和結(jié)點力列向量A(u)均為位移向量{u}的非線性函數(shù)，須采用直接時域積分方法進行非線性熱氣動彈性分析。本文基于Nastran軟件中的非線性瞬態(tài)響應直接求解模塊SOL400中的非線性有限元求解程序，采用DMAP語言對其進行二次開發(fā)，將溫度場分析和Van Dyke修正活塞理論非定常氣動力計算引入SOL400非線性瞬態(tài)響應直接求解模塊，編制了非線性熱氣動彈性分析程序，將Van Dyke修正活塞理論求解非定常氣動力的高效性和SOL400求解結(jié)構(gòu)非線性瞬態(tài)響應的高效性充分結(jié)合，避免了非線性瞬態(tài)響應復雜有限元程序的編制，技術(shù)路線如圖 1所示。

圖1 非線性熱氣動彈性分析技術(shù)方案Fig.1 Technical scheme of nonlinear thermoaeroelastic analysis

SOL 400引入非定常氣動力的DMAP二次開發(fā)技術(shù)路線如圖 2所示，并采用PCL語言在Patran軟件中開發(fā)非線性熱氣動彈性分析嵌入式前后處理界面，如圖 3所示，建立了工程實用的非線性熱氣動彈性分析工具。

圖2 SOL 400引入非定常氣動力DMAP二次開發(fā)技術(shù)方案Fig.2 Technical scheme for SOL 400 secondary development by introducing unsteady aerodynamics with DMAP language

圖3 非線性熱氣動彈性分析界面Fig.3 Nonlinear thermoaeroelastic analysis interface

3 數(shù)值算例

針對石墨-環(huán)氧復合材料壁板標準算例，壁板長度、寬度、厚度分別為0.381m、0.3048m、1.3e-3m，位移邊界為四邊固支，壁板力學性能參數(shù)和有限元模型分別如表 1和圖 4所示。

表1 壁板力學性能參數(shù)Table1 Mechanical property parameters of panel

圖4 壁板有限元模型Fig.4 Finite element model of panel

本文采用與文獻標準算例相同的Van Dyke修正活塞理論計算非定常氣動力，進行了復合材料壁板典型溫度工況的非線性熱氣動彈性分析。取壁板展向50%站位3/4弦線位置處結(jié)點的位移和速度響應與標準算例進行對比驗證，maxW表示最大位移振幅。

1）常溫工況

常溫工況取空氣密度為1.225kg/m3，馬赫數(shù)為2.0，采用本程序分析了不同飛行速度下極限環(huán)幅值隨無量綱動壓的變化規(guī)律，如圖 5所示。由圖 5可知，本文計算得到的無量綱臨界顫振動壓為275，文獻中無量綱臨界顫振動壓為250，誤差不超過10%，且隨著無量綱動壓增大，極限環(huán)幅值增大，本文與文獻結(jié)果變化趨勢一致。研究結(jié)果表明，本程序適用于復合材料壁板結(jié)構(gòu)常溫工況氣動彈性分析。

圖5 壁板極限環(huán)幅值隨無量綱動壓變化規(guī)律Fig.5 Variation of limit cycle amplitude of panel with dimensionless dynamic pressure

2）均勻溫度場工況

均勻溫度場工況取空氣密度為1.225kg/m3，馬赫數(shù)為2.0，考慮四種均勻溫度場，溫升分別為ΔT=Tcr,2Tcr,3Tcr,4Tcr,屈 曲 臨 界 溫 升 為Tcr= 21.8°C。采用本程序計算得到不同飛行速度下的非線性熱氣動彈性響應如圖 6-圖 9所示，分別對應了LCO，周期，非周期和混沌響應。非線性熱氣動彈性響應無量綱最大幅值與文獻對比如表2所示。

表2 非線性熱氣動彈性響應無量綱最大幅值（ max/W h）對比 Table2 Comparison of non-dimensional maximum amplitudes of nonlinear thermoaeroelastic responses

圖6 ΔT = Tcr，LCOFig.6 ΔT = Tcr , Limit cycle oscillation LCO

圖7 ΔT = 2Tcr ，周期Fig.7 ΔT = 2Tcr, periodic vibration

圖8 ΔT = 3Tcr ，非周期Fig.8 ΔT = 3Tcr, aperiodic vibration

圖9 ΔT = 4Tcr ，混沌Fig.9 ΔT = 4Tcr, chaos

由分析可知，本程序可以對壁板結(jié)構(gòu)LCO，周期，非周期和混沌等復雜非線性熱氣動彈性響應進行有效預示，且本文與文獻結(jié)果一致性很好，誤差不超過6%，表明針對Nastran SOL 400多學科非線性求解模塊中的非線性瞬態(tài)響應求解器的DMAP二次開發(fā)流程是正確的。

4 結(jié)論

本文采用Van Dyke修正活塞理論分析非定常氣動力，針對商用軟件Nastran SOL 400多學科非線性求解模塊中的非線性瞬態(tài)響應求解器，增加溫度場和非定常氣動力求解模塊，采用PCL和DMAP語言，基于Patran/Nastran二次開發(fā)平臺，編制了非線性熱氣動彈性分析程序，以復合材料壁板為研究對象，分析了其在典型溫度工況下的非線性氣動彈性響應，標準算例對比驗證了本文方法的正確性。

研究結(jié)果表明，常溫工況下，本文計算得到的臨界顫振動壓與文獻結(jié)果的誤差不超過10%，且變化趨勢一致，且可以對復合材料壁板結(jié)構(gòu)LCO，周期，非周期和混沌等復雜非線性熱氣動彈性響應進行有效預示。