基于修正活塞理論的壁板非線性熱氣動彈性分析

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

（北京強度環境研究所 可靠性與環境工程技術重點實驗室，北京100076）

0 引言

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

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

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

1 Van Dyke修正活塞理論

根據經典活塞理論，翼面上某點壓力與該點下洗速度之間滿足如下關系[20]

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

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

式中，M為馬赫數，a為當地音速，則活塞理論可轉化為

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

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

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

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

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

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

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

3 數值算例

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

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

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

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

1）常溫工況

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

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

2）均勻溫度場工況

均勻溫度場工況取空氣密度為1.225kg/m3，馬赫數為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

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

4 結論

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

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