C/SiC蓋板式熱防護結構動力學建模方法

孫 兵，蔡巧言，解海鷗，郭鵬飛，石 偉

（中國運載火箭技術研究院研究發展中心，北京，100076）

0 引 言

近年來，可重復使用運載器（Reusable Launch Vehicle，RLV）成為國內外研究熱點，其中對熱防護系統的研究是可重復使用運載技術研究中非常重要和關鍵的一項。RLV在工作時遭遇嚴峻的氣動加熱，熱防護系統（Thermal Protection System，TPS）的材料和結構需要能夠抵抗極端惡劣的氣動加熱情況。

第1代可重復使用運載器航天飛機使用粘接的剛性隔熱瓦[1]作為其在迎風面的熱防護系統，構型如圖1所示。下一代低成本可重復使用運載器對熱防護系統提出了質量輕、成本低、健康管理能力高、可維護性好等需求，粘接的剛性隔熱瓦已不能滿足這些需求。歐洲借鑒并評估了航天飛機等以往TPS方案設計的經驗，提出一種采用C/SiC先進復合材料的蓋板式TPS方案，構型如圖2所示[2]。蓋板式TPS一般安裝在可重復使用運載器的迎風面以抵抗大熱載，其在運輸、飛行等環境中都會受到動載荷作用，對 C/SiC蓋板式TPS進行模態參數預示是十分必要的。在 PRE-X項目[3]中，歐洲對C/SiC蓋板式TPS的設計做了諸多研究工作，包括基于有限元模型對單塊 C/SiC蓋板式熱防護結構裝配體進行模態參數預示并開展了一系列結構動力學試驗研究。

圖1 航天飛機剛性隔熱瓦[1]Fig.1 Rigid Insulation Ceramic of Space Shuttle

圖2 蓋板式熱防護方案概念[2]Fig.2 Conceptual Demonstration of Shingle TPS

中國也開展了關于C/SiC蓋板式TPS結構動力學方面的研究。文獻[4]為分析C/SiC面板的熱氣彈響應，先對C/SiC面板進行了模態分析。郝秉磊、殷小瑋等[5]建立了 C/SiC陶瓷復合材料螺栓連接有限元模型，并通過試驗驗證了模型的正確性。但是，針對 C/SiC蓋板式熱防護結構裝配體的結構動力學研究還處在起步階段，相關方面的文獻資料有限。本文將對單塊C/SiC蓋板式熱防護結構裝配體的結構動力學建模方法開展研究，不僅可為其結構動態設計提供參考，還可為整個飛行器的模態分析、模態試驗等相關工作提供支撐。

1 理論模型的建立

C/SiC蓋板式TPS的固有振動特性需先建立合理的結構動力學模型，由圖 2知單塊蓋板式熱防護結構為夾層結構，建立的分析模型如圖3所示，此結構模型是由兩塊平板和彈簧組成的二維系統。C/SiC蓋板和機身蒙皮是兩塊彈性板，而隔熱層簡化為彈簧系統。因為隔熱層相對較厚而模量相對較低，所以當 C/SiC蓋板和機身蒙皮都有撓曲變形時，隔熱層結構能被壓縮。模型假設上下兩塊板沿長度方向的每一點都被彈簧連接。此假設是合理的，因為在實際情況下無論通過緊固件連接還是膠接，均使這3層結構的面接觸在一起。

圖3 單塊蓋板式熱防護結構模型Fig.3 Structural Model of a Single Shingle TPS a，b—單塊蓋板式熱防護結構的長和寬

不考慮外力和阻尼，利用經典的拉格朗日函數能量法推導此系統的運動方程。系統的總勢能為機身蒙皮勢能C/SiC蓋板勢能和隔熱層勢能的總和：

根據板的振動理論可知，假設板在振動時的撓曲變形是符合胡克定律的微小彈性變形，這種變形明顯小于板厚。因此，這里只考慮了板的彎曲勢能，即假設板在發生撓曲變形時，沒有明顯伸長。

系統的總動能 TTOT為機身蒙皮動能、C/SiC蓋板動能和隔熱層動能的總和：

式中

將式（1）～（9）代入式（11）中得到無阻尼系統自由振動運動方程：

將式（12）寫成矩陣形式有：

式中 M為質量矩陣；K為剛度矩陣；q為時間。

可以用戈遼金法、瑞利-瑞茲法或有限元法等數值方法求解方程的近似解。

2 有限元模型的建立

在忽略阻尼影響的情況下，有限元模型要求能夠準確反映結構的質量分布特性和剛度分布特性。材料參數、單元類型和幾何構型等均將影響結構的剛度分布特性，從而影響模態分析結果。此外，邊界條件、連接件的預緊力等也將對結構的固有振動特性產生影響。以典型蓋板式TPS的簡化模型為例進行研究，構型如圖4所示。

圖4 典型蓋板式TPS構型Fig.4 Typical Shingle TPS Structure

圖4中模型主要由4個部分組成：C/SiC蓋板、隔熱層、機身蒙皮、連接件。具體尺寸及材料參數見表 1[6]，其中彈性模量均為室溫下的參數，所用材料均為均質各向同性材料。

表1 尺寸及材料參數Tab.1 Dimensions and Material Physical Properties

2.1 單元類型

為建立合理的蓋板式TPS有限元模型，下面將從賦予結構不同單元類型的角度開展分析。隔熱層選用實體單元，機身蒙皮和蓋板由于厚度非常薄可用實體單元也可用殼單元建模。分析模型如圖5a所示，其中模型a中的單元全部為8節點六面體實體單元，單元總數為4800個（圖6a）；在四邊自由的邊界條件下，利用Lanczos迭代法計算得到模態基頻為263.68 Hz。模型b中蓋板和機身蒙皮采用殼單元，隔熱層選用實體單元，單元總數為4800個（見圖6b）；邊界條件為四邊自由，經 Lanczos迭代法計算得到模態基頻為254.26 Hz。

圖5 模型示意Fig.5 Concepts of Model

圖6 模型的有限元模型Fig.6 Illustration of Finite Element Model

表2為模型質量統計結果，由表2可知，模型a和模型b的質量與實際質量相比存在誤差，誤差分別為 1%和 0.4%，這是由于結構被離散導致的。單元的一致質量矩陣為

式中 ρ為密度；N為單元位移插值函數矩陣；Ω表示體積；下角標e表示單元。

表2 模型質量統計Tab.2 Model Mass Details

而在計算模型a和模型b的模態基頻時，單元質量矩陣為集中質量矩陣。集中質量矩陣是假定單元的質量集中在節點上，因而得到的質量矩陣是對角線矩陣，這將簡化計算提高效率，并且對低階頻率的精度影響很小。

實體單元將單元一致質量矩陣eM轉換為單元集中質量矩陣的一種常用方法如式（15）所示。

式中eV為單元的體積。

式中 下標 d為單元在幾何空間的維數；W為單元的質量；dI為單位矩陣；eV為單元的體積。

殼單元將單元一致質量矩陣eM轉換為單元集中質量矩陣的常用方法是忽略 M 中對應于轉動自由

e度的元素。對于中與位移自由度相關的元素可采用式（15）進行處理。實際上，實體單元和殼單元集中質量矩陣的不同，本質上是由于節點自由度不同引起的，從而使得模型a和模型b的質量誤差不同。本例中機身蒙皮和蓋板采用殼單元質量分布特性更接近真實情況。

提取模型a和模型b中蓋板和隔熱層交界面相同位置角節點的剛度矩陣元素，位置在圖5圓點處，節點編號為i?？傮w坐標下，由于模型a中全部為實體單元，每個節點有3個方向的自由度，而模型b由于有殼單元，其節點自由度為6。按節點劃分結構的總剛度矩陣，提取的iiK如式（17）和式（18）所示。

式中 上標a和b為模型編號；下標i表示角節點編號。

因為節點位置為交界面的角點上，所以此節點只為交界面上下兩個單元共用的節點。根據直接剛度法可知，的值是2個實體單元在角節點i處相應剛度元素的疊加；的值是1個實體單元和1個殼單元在角節點 i處相應剛度元素的疊加。觀察矩陣可以看到和主對角線上的值基本在同一量級上，并且的剛度元素稍大。

2.2 裝配建模方法

由圖 4知此裝配體為夾層結構，載荷主要通過上下板面接觸及連接件傳遞。接觸過程中2個物體在接觸界面上的相互作用是復雜的力學現象，接觸界面高度非線性。這種非線性來源于2個方面：a）接觸界面的區域大小和相互位置以及接觸狀態事先未知，并且隨時間變化，需在求解過程中確定；b）接觸條件非線性，接觸條件的內容包括：接觸物體不可相互侵入；接觸力的法向分量只能是壓力；切向接觸的摩擦條件。這些條件區別于一般約束條件，其特點是單邊性的不等式約束，具有強烈的非線性[7]。

假設預緊力足夠大使得連接件在振動過程中可以忽略非線性影響，按以往的工程經驗這種連接屬于緊密連接[8]，連接面兩端位移可以認為是協調的，建模時可采用共節點方法或綁定接觸模式進行裝配建模。以模型b為例分別使用上述兩種裝配建模方法進行模態計算，結果如表3所示。共節點方法和綁定接觸模式均是使裝配體組件在接觸面間不存在切向相對滑動或法向相對分離，但是在交界面上模擬力的傳遞方式不一樣。

表3 不同裝配建模方法計算結果Tab.3 Simulation Results of Modeling

共節點是通過直接剛度法，將下標相同的剛度元素疊加在一起，從而實現不同類型單元的結合。其物理意義是：整個結構抵抗某種變形的能力等于各單元抵抗該變形能力之和。

綁定接觸模式則通過一定的接觸算法實現力的傳遞，但計算時只需要迭代一次，屬于線性接觸類型，從而可在模態分析這種線性計算中使用。利用有限元分析軟件ANSYS Workbench，在模型b的交界面生成接觸單元CONTA 174和TARGE 170從而建立接觸對，根據罰函數法進行接觸計算。

罰函數法采用一個接觸“彈簧”在兩個面間建立關系，如圖7所示，這個彈簧的剛度稱為接觸剛度[9]。在力F的作用下，彈簧被壓縮Δ（侵入數值），且滿足平衡方程

式中 k為接觸剛度。

圖7 罰函數法的接觸彈簧Fig.7 The Contact Spring of Penalty Function Method

式中nK為法向接觸剛度；nu為接觸間隙大小。根據庫倫法則，摩擦應力如下（zτ類同）：

式中sK為切向接觸剛度；yu為y向接觸滑動距離；μ為摩擦系數。

本算例中初始法向接觸剛度和切向接觸剛度分別設置為10 N/m2、1 N/m2。采用這種綁定接觸模式的裝配建模方法，計算出模型 b的模態基頻結果為262.76 Hz。

在本算例中采用共節點和綁定接觸方法對模型進行裝配建模，由表 3可知采用這兩種裝配建模方法模態基頻的計算結果相差不大。使用共節點方法進行裝配時，要求相鄰組件網格節點在交界面上一一對應，這將給各組件的網格劃分帶來限制，而使用綁定接觸方法則不會帶來這樣的問題。但是采用綁定接觸方法會產生接觸單元，需額外進行接觸計算，在求解全飛行器的動力學問題時將降低計算效率。

3 結 論

本文分別通過解析方法和有限元法建立了蓋板式熱防護結構的理論模型和有限元模型。建立理論模型的優點是物理量的意義簡單直觀，缺點是需進行大量簡化與假設。雖然有限元法是高效且精確的工具，但是針對不同的結構需要重新建立不同的計算模型。對于文中算例，機身蒙皮和蓋板采用殼單元模擬質量分布特性更接近真實情況；采用共節點和綁定接觸方法對模型進行裝配建模，模態基頻的計算結果相差不大。文中對 C/SiC蓋板式熱防護結構動力學建模方法的研究，可為建立真實飛行器上的 C/SiC蓋板式熱防護結構動力學模型提供參考。