氣動力熱作用下飛行器結構分析的分區耦合CG/DG有限元方法

趙建寧，魏 東，杜雁霞，劉冬歡，*，桂業偉

(1. 北京科技大學 數理學院應用力學系，北京 100083；2. 中國空氣動力研究與發展中心，綿陽 621000)

0 引 言

近空間高超聲速飛行器是當前國內外的研究熱點之一，其發展水平也是一個國家綜合國力的重要體現[1]。近空間高超聲速飛行器的研發涉及到氣動布局、氣動加熱環境模擬、結構熱防護、發動機研制、航空航天材料、高能燃料、全球定位、矢量控制等多學科技術問題，其中熱防護結構是起到決定性作用的一項關鍵技術[2-4]。對熱防護結構的設計及安全評估涉及到氣動熱力耦合分析，即高超聲速飛行引起的氣動力和氣動熱對飛行器結構的傳熱防熱及熱力耦合行為的影響[5-6]，以及由此引起的熱氣彈問題[7-8]。

國內外研究者針對氣動熱力耦合問題分析開展了大量的研究工作，取得了一系列豐富的研究成果。桂業偉等給出了氣動力/熱與結構熱響應多場耦合問題的數據流程，深入探討了多物理場耦合的策略與方法[9]。郭亮等建立了一種針對流場-熱-結構的多場耦合分析方法，實現了對固體隔絕內外流場溫度動態變化問題的仿真分析[10]。張章等利用流固及熱固單向耦合的方法分析了考慮高超聲速流場氣動壓力和氣動熱作用下空間再入充氣結構的特性變化[11]。張勝濤等基于松耦合分析策略框架分析了高超聲速流場-熱-結構耦合問題，采用自適應耦合計算時間步長、混合插值策略和復雜外形網格變形等方法，實現了多場耦合分析平臺[12]。Huang等建立了一體化熱氣彈耦合的計算框架，基于松耦合方法分析了復合材料雙曲淺殼的溫度場及熱應力分布[13]。Li等基于有限體積法建立了高超聲速飛行器圓柱形前緣的空氣動力-結構-傳熱一體化分析方法[14]。Huang等建立了針對熱防護系統的流-固-熱緊耦合數值方法，與松耦合相比，緊耦合需要增加內部迭代過程。研究結果表明，緊耦合方法可以降低時滯效應并增加時間步長，但是由于迭代量的增加實際計算時間花費并未降低[15-16]。Chen等基于松耦合方法建立了時間步長可自適應改變的高超聲速飛行器前緣結構的氣動-熱-力耦合分析方法[17]。

特別的，為了提高結構熱防護的效率，不同的飛行器結構部位采用的防熱結構和使用的防熱材料體系也往往不同，因此產生了大量的異種材料裝配界面。裝配界面會導致兩方面的問題，一是非完好接觸的裝配界面會產生阻礙熱量流動的接觸熱阻，從而對結構溫度場產生影響，二是由于界面兩側不同材料的熱物性參數可能相差很大，由此帶來了嚴重的熱失配問題。由于界面接觸熱阻受界面溫度和界面接觸應力的影響，因此形成了復雜的耦合關系，如圖1所示。

圖1 考慮接觸熱阻的熱力耦合問題Fig. 1 Thermomechanical coupling problem with thermal contact resistance

對于含接觸熱阻的多場耦合問題，采用間斷伽遼金有限元方法進行處理是一種很自然的選擇[18-19]。間斷伽遼金方法最早是為了求解中子輸運的雙曲型方程[20]，間斷伽遼金有限元法既能夠像經典有限元方法一樣處理復雜邊界問題，又吸收了有限差分法的優點，適用于顯式求解。間斷伽遼金有限元方法允許單元的插值函數存在間斷，可以采用不一致的網格劃分以及利用間斷的分片多項式構造近似函數和權函數空間，在構造高階單元時具有很大的靈活性，可以較好地進行網格加密及升階和并行化[21]；間斷伽遼金有限元方法不再強制單元之間必須協調連續，因此可以避免諸如薄膜自鎖、剪切自鎖、體積自鎖等現象，且提高了解的計算精度和收斂性[22]；數值通量的引入進一步消除了間斷處的虛假數值振動，相當于計算流體力學中的人工黏性，因此在處理高梯度以及間斷問題上很有優勢[23]；同時由于采用間斷的插值函數，因此非常適合含間斷問題的處理。但是間斷伽遼金有限元方法的最大弱點在于計算量大，特別是對于三維問題[24]。因此，將間斷伽遼金有限元方法和傳統連續伽遼金方法結合起來就成為一種更好的選擇[24]。Nguyen等研究了基于間斷伽遼金有限元方法的氣動熱力耦合問題中界面的自適應處理方法[25]。劉冬歡等研究了接觸熱阻對疏導式熱防護結構傳熱防熱效果的影響[26]。關于氣動熱力多場耦合問題雖然已經有了大量的研究，但是針對考慮接觸熱阻的氣動熱力耦合問題，國內外相關的研究并不多見。

由于非完好接觸界面的存在，產生了由界面接觸熱阻誘導的非線性熱力耦合、由高溫引起的材料熱物性的非線性等復雜的非線性效應，傳統的三維有限元計算方法在求解此類問題時遇到了嚴重的挑戰，不僅計算規模大、計算效率低而且計算精度有限，給工程設計帶來了極大的不便，成為制約我國新型高超聲速飛行器結構設計及安全評估的瓶頸問題之一。在國家數值風洞工程項目的支持下，本文建立了分區耦合的間斷/連續伽遼金有限元方法及其算法框架，并編制相應的有限元計算軟件，形成適應于工程大規模計算的含多材料搭接界面結構的三維傳熱、熱力耦合問題研究的計算分析能力，這對提高我國高超聲速飛行器結構性能預測水平及優化設計等具有重要意義。

1 有限元格式及計算流程

1.1 氣動節點數據向有限元節點的轉換方法

雖然國內外針對流固界面數據交互方法進行了大量的研究，但是很多往往過于理論化，采用非常復雜的插值方法來重構界面物理場。Pidaparti通過等參映射實現結構分析和流體分析界面的數據轉換[27]，三點挑方法將氣動網格節點上的氣動力按照靜力等效原則分配到鄰近的三個結構點上，三點挑方法可以推廣到多點挑方法，本質上是要求離氣動點近的結構點多分配一些載荷而遠的少分配一些[28]。Liu等通過一個基于表面樣條理論的轉換矩陣將流場節點位移列向量用結構有限元節點位移列向量來表示，同時基于虛功原理將流場網格上的壓力轉換到結構網格上的有限元節點力[29]。劉深深等提出一種基于幾何尺度進行歸一化的高超聲速飛行器多場耦合數據傳遞新方法[30]。張建剛等采用樣條曲面擬合的方法得到翼面壓力分布曲面，由該曲面得到有限元節點上的壓力值，再在有限元模型單元上積分得到有限元節點載荷供強度設計使用[31]。這種數據插值擬合效果的好壞很大程度上取決于數據點的分布情況，若某個區域數據點分布較少，則很難通過插值的方法恢復該區域，因為這個區域已知的信息量不足以高精度地恢復數據。由于多場耦合分析涉及多次迭代計算，要求氣動節點數據向有限元節點轉換的效率一定要高，因此本文從工程實用的角度出發，開發流固界面數據轉換軟件。基本的思想是確定待定有限元節點在氣動網格內的歸屬單元，即有限元節點被哪個氣動單元所包圍。這可以通過比較節點坐標范圍的方法進行判斷，接下來基于有限元節點在氣動單元內部的不同位置，基于多點挑或插值的方法來確定有限元節點的場量值。

1.2 溫度場分析的連續/間斷有限元格式

考慮邊界為 ? Ω 的區域 Ω上的橢圓型熱傳導方程：

其中，▽為梯度算子，k是熱傳導系數矩陣，一般是對稱的，T為溫度，f是熱源函數。問題的Dirichlet邊界條件和Neumann邊界條件分別為：

其中， ?ΩD和 ?ΩN分別表示Dirichlet邊界和Neumann邊界，邊界 ?Ω 的單位外法線為n， 熱流密度為q，邊界上給定的溫度和熱流密度分別為和。

等效積分弱形式的熱傳導方程可寫成：

其中，熱流密度為qe=?ke▽Te。通過分部積分可以得到其等效形式：

其中，單元e的邊界 ? Ωe上的外法線為ne。

上式可見，單元邊界? Ωe上的積分項出現在單元平衡方程中，在間斷伽遼金有限元方法中單元內部邊界上允許場變量出現間斷，該積分項不能隨意消去，可通過定義在單元邊界上的Bassi-Rebay型數值通量來考慮這種間斷效應，間斷伽遼金有限元方法就是通過將單元邊界熱流qe用其對應的熱流數值通量進行替換來實現的，即單元的間斷伽遼金有限元方程為：

其中，熱流數值通量為：

這里，上標“eb”代表單元e的鄰居單元。而在單元外邊界上，數值通量取為：

同時必須在邊界上加入穩定項，這里將單元內部邊界基本場變量的跳變引入到平衡方程中，稱為穩定項，此時單元間斷伽遼金有限元方程可寫成：

其中，溫度穩定化參數 τT是O(‖ke‖/h)量級的正數，這里 ‖ke‖是熱傳導系數矩陣的某種范數，h是單元的特征尺寸。

在引入界面的非完全熱接觸條件時，首先需要將單元間斷伽遼金有限元方程中的溫度穩定化參數賦零，同時將界面的熱流數值通量由非完全接觸的定量描述即接觸熱阻的定義式給出：

同時為了體現這種強的間斷效應，在存在接觸熱阻的單元內邊界上要將穩定化參數賦零。

假設單元溫度場可以表示為：

這里，NTe為溫度場插值函數，Te為單元節點的溫度列向量。權函數取為ve=(NTe)T，代入單元有限元方程，并引入氣動熱載荷和輻射邊界條件，經過推導并組裝可得用溫度節點列向量為唯一未知量的結構總體傳熱有限元方程組：

在引入強制溫度邊界條件并求解此方程組即可得到結構的溫度場，如果考慮材料屬性的溫度相關性或者接觸熱阻的溫度相關性，則需要進行迭代求解，進一步可得到熱流密度場等相關場量。這樣就自然地將接觸熱阻引入到單元平衡方程中來，避免了采用界面單元的麻煩。而在不存在接觸熱阻的地方，可以采用經典的連續伽遼金有限元方法，無需穩定項和數值通量項，從而大大提高了計算效率。

1.3 位移場分析的連續/間斷有限元格式

結構位移場分析的平衡方程可寫成：

其中，σ為應力向量，f為體力向量，▽為微分算子。

平衡方程的等效積分弱形式可寫成：

其中，we為 定義在區域 Ω上矢量形式的權函數。

通過分部積分可以得到其等效形式：

通過將單元邊界應力 σe用應力數值通量進行替換即可得到單元應力分析的間斷伽遼金有限元方程：

和溫度場分析時引入接觸界面條件類似，應力場分析時需要引入位移的不可穿透條件。不可穿透條件意味著接觸界面兩側的節點不能互相穿透對方的表面而進入其內部，這在數值計算中是一個很難處理的不等式約束問題。一般有兩種方法可用于引入不可穿透條件，一種是采用界面單元來耦合接觸區域，通過給定界面單元的剛度與節點相對位置之間的關系來描述各種各樣的接觸狀態；另一種是在計算格式中直接施加位移約束，這需要首先進行接觸搜索確定接觸對，進而采用諸如罰函數方法、拉格朗日乘子法、約束函數方法等對接觸對施加位移約束。與此同時，在許多工程實際問題中罰函數方法已經被證明是很高效的。另一方面，這也同間斷伽遼金有限元方法實現無接觸單元界面位移連續的穩定項是一致的，因此這里我們采用罰函數方法來引入接觸界面的不可穿透條件。此時有：

其中，gN表示單元e與其鄰居單元eb之間的初始預留間隙。

假設單元e的位移場和應變場通過插值可以表示為：

其中，Nue為 單元位移場的插值函數矩陣，Ue為單元節點位移列向量，Be=▽TNue為幾何矩陣。

進行應力場分析時，需要考慮溫度產生的變形，此時的本構方程為：

其中，De是彈性矩陣，T0是計算熱變形的參考溫度，αe是材料的熱膨脹系數向量。

引入強制位移邊界條件即可求解整體有限元方程組，得到結構位移場，進一步經過應力磨平等處理即可得到高精度的應力場。

值得注意的是，本文假設結構的變形處于小變形狀態，因此熱應力對結構的影響僅體現在等效載荷列向量中，并未考慮其對結構幾何剛度的影響。

1.4 計算流程及軟件開發

求解多場耦合問題的基本方法有緊耦合和松耦合兩種，本文采用松耦合的方法通過迭代得到耦合問題的解，具體的計算流程圖如圖2所示。

圖2 含接觸熱阻的熱力耦合問題計算流程圖Fig. 2 Computational algorithm for solving thermomechanical coupling problems with thermal contact resistance

首先進行流固界面載荷轉換，輸入氣動熱力計算結果和結構有限元網格信息，利用界面載荷轉換程序得到有限元節點或者單元上的氣動熱流和氣動力信息。

其次輸入必須的各種數據，包括幾何信息，比如有限元網格的單元組成和節點坐標；材料屬性信息，比如隨溫度變化的熱學參數（熱導率、熱膨脹系數），力學參數（彈性模量、泊松比、塑性參數），接觸熱阻模型；載荷信息（經過界面場量轉換軟件得到的有限元節點上的氣動力和氣動熱流密度、結構位移約束類型及所在的位置）；控制參數（比如分析類型、收斂準則、自適應策略、初始化參數等）。

基于初始化的參數進行溫度場的試算，并判斷當前溫度場是否收斂：如不收斂，則對溫度場插值單元進行升階或者加密網格并重新計算；如收斂，則繼續進行位移場和應力場的計算。得到位移場和應力場結果之后，需要對位移場結果進行收斂性判斷：如不收斂，則進行單元插值函數升階或網格加密再計算；如收斂，那么基于當前的界面接觸應力場和接觸熱阻模型，對界面接觸熱阻進行更新，并與之前的界面接觸熱阻進行收斂性判斷：如不收斂，則需要基于當前的新界面接觸熱阻重新回到溫度場計算；如收斂，則針對計算得到的結果進行后處理并輸出結果，計算結束。

基于上述有限元格式和計算流程圖編制了相應的計算程序，開發了相配套的軟件，并作為一個功能模塊集成到國家數值風洞軟件群中。

2 算法驗證

2.1 問題描述

為驗證算法的準確性，構建了如圖3所示由四面體單元構成的組合桿件，將本文結果同有限元軟件仿真結果對比分析。桿件1長度為0.4 m，材料熱導率為30 W/m·K，熱膨脹系數為10×10?6/K，彈性模量為200 GPa；桿件2長度為0.2 m，材料熱導率為100 W/m·K，熱膨脹系數為15×10?6/K，彈性模量為200 GPa。泊松比均為0.3。組合桿件左右兩端分別給定恒溫300 K和700 K，界面熱阻為1×10?3m2·K/W。左右兩端位移固定，熱變形參考溫度300 K。

圖3 組合桿件模型有限元網格圖Fig. 3 The mesh of the composite rod model for finite element analyses

2.2 結果對比

圖4～圖6分別給出了組合桿件溫度場、總位移場(USUM)和等效應力場(SEQV)基于本文有限元程序結果和通用有限元軟件結果的對比。同時選取了組合桿件軸線方向截面中心處的一條路徑，通過場量沿該路徑的分布圖比較不同界面接觸熱阻R下兩種方法的計算結果誤差，見圖7，并在組合桿件中心選擇同一位置比較等效應力的結果誤差，見表1。

表1 組合桿件中點位置的等效應力Table 1 Equivalent stress at the center of the composite rod

圖4 組合桿件溫度場的對比Fig. 4 Comparisons of the temperature field of the composite rod

圖6 組合桿件等效應力場的對比Fig. 6 Comparisons of the equivalent stress of the composite rod

圖7 不同接觸熱阻時軸向溫度及總位移的對比Fig. 7 Comparisons of the axial temperatures and total displacements with different thermal contact resistances

從以上結果對比可以看出，接觸熱阻的存在導致溫度場分布在界面處產生跳變，接觸熱阻越大，溫度的跳變越大。沿軸線路徑上的溫度場和位移場以及桿件中點處的等效應力數值，本文程序結果和理論解及通用軟件結果一致，桿件中點處的等效應力的相對誤差不超過0.4%。除此算例外，本文有限元程序在設計開發中，基于軟件工程的思想，對所有的邊界材料、材料屬性、結構形式、載荷形式進行了全方位的測試，本程序結果均與通用有限元軟件結果或者理論解完全吻合，這充分說明了本文程序的精度和可靠性，可進一步用于實際工程結構的數值仿真。

圖5 組合桿件總位移場的對比Fig. 5 Comparisons of the total displacement of the composite rod

3 工程算例

3.1 問題描述

圖8為高速飛行器的水平舵面模型圖，整體結構由翼前緣和翼后舵結構（蒙皮、蜂窩夾心、結構骨架組成）構成，材料均為GH99。翼根處弦長約537 mm，翼梢處弦長約380 mm，翼展約317 mm，翼面由三個折面組成。針對飛行器強烈的氣動加熱現象，首先基于給定的冷壁熱流和恢復焓對水平舵結構進行熱分析，進而同時考慮氣動力和氣動加熱的作用，對結構的熱強度進行計算分析。溫度場分析時需考慮翼前緣與翼后舵之間的界面接觸熱阻的影響。

圖8 舵面模型結構圖Fig. 8 The geometry of the rudder model

3.2 流固界面氣動熱力轉換

基于1.1節的轉換方法得到的轉換結果如圖9～圖11所示。結構溫度場和應力場分析時需要給定氣動熱流、氣動力邊界條件，相關邊界條件由CFD計算給出，提供的環境數據為翼面上各個氣動節點的氣動力、冷壁熱流和恢復焓。由于氣動網格節點和有限元分析網格的節點是不重合的，因此首先需要將CFD氣動點的氣動力、冷壁熱流以及恢復焓轉換到有限元分析時的單元結構節點上。

圖9 氣動和結構分析時計算網格的對比Fig. 9 Comparisons of meshes for CFD and FEA computations

圖10 轉換前后氣動和結構分析氣動壓強的對比Fig. 10 Comparisons of aerodynamic pressures for the CFD and FEA computations after conversion

圖11 轉換前后氣動和結構分析氣動熱流的對比Fig. 11 Comparisons of aerodynamic heat fluxes for the CFD and FEA computations after conversion

可以看出，飛行器飛行過程中，水平舵翼前緣承受較大的氣動熱和氣動壓強，基本在2 000 kW/m2和170 kPa左右，而舵面數值偏小，與前緣相差約為一個數量級。

由于本算例主要針對面分布的熱流和氣動壓強進行轉換，從快速確定有限元節點場量的角度出發，采用簡化的轉換算法，將距離有限元節點最近的氣動點的場量值作為該有限元節點的場量值。

從轉換前后氣動和有限元的場量云圖可以看出，在大部分區域的轉換效果都很好，而在翼前緣區域由于轉換前后網格差異較大，導致轉換誤差稍微大一些，同時由于該區域范圍較小，因此轉換誤差對結構變形和應力的影響是可以忽略的。

3.3 邊界條件以及熱阻模型

進行熱分析時對翼舵外表面施加給定熱流邊界條件，根據CFD計算得到的冷壁熱流qc和恢復焓hr計算得到施加在結構上的熱壁熱流qn為：

這里考慮了舵面輻射散熱效果，ε為Stefan-Boltzmann常數，取為5.67×10?8W/m2K4， σ為輻射系數，取為0.8，T為待求的壁面溫度。焓值計算時的參考溫度為280 K，壁面溫度和參考溫度條件下的焓值hw和hw0分別為：

其中，cp為空氣的定壓比熱容，在250～2500 K之間可通過數據擬合為：

溫度場分析時翼舵底部安裝面給定溫度280 K，輻射環境參考溫度280 K，由于給定的熱流邊界條件與待求溫度T相關，因此這是一個非線性溫度場問題，依據上述程序進行迭代求解。

在對翼舵強度分析時，外表面施加給定的氣動力載荷，翼舵安裝面位移固定，熱變形參考溫度280 K，而界面接觸熱阻與界面應力密切相關[32-33]：

其中， σc為界面接觸應力，MPa。當應力較小時，界面接觸熱阻較大，隨界面應力的增大，界面接觸熱阻隨應力指數型減小。

3.4 計算結果

3.4.1 結構溫度場

基于上述有限元方法，在不考慮界面接觸熱阻時得到的結構溫度場如圖12所示。

圖12 無接觸熱阻時舵面溫度場Fig. 12 The temperature field of the rudder without the thermal contact resistance

從以上計算結果可以看出，由于水平舵兩面氣動熱邊界條件的差異導致了溫度場分布的不同，其中前緣最高溫度達到1 276 K，舵面整體溫度大致在600～1 100 K之間。實際工程應用中，溫度過高對結構變形影響較為顯著，因此針對升溫引起的結構應力變化的研究很有必要。

不同接觸熱阻時舵面溫度場如圖13所示，翼舵結構界面兩側節點溫度隨接觸熱阻的變化趨勢如圖14所示。

圖13 不同接觸熱阻Rc時的舵面溫度場Fig. 13 Temperature fields of the rudder with different thermal contact resistances Rc

圖14 界面接觸熱阻對接觸界面兩側節點溫度的影響Fig. 14 Variations of interface temperatures with the interface thermal contact resistance

界面接觸熱阻在組合結構熱傳導過程中起到阻礙熱量傳遞的作用，會造成界面兩側溫度的跳變，從本翼舵算例結果可以看出，界面接觸熱阻越大界面兩側溫度跳變越明顯，但對于整體溫度場分布的影響不明顯，這是因為本算例中熱流邊界條件與輻射邊界同時施加在翼舵外表面上，氣動熱和輻射相互作用的影響對于翼舵表面溫度場的分布占據主導作用，因此界面接觸熱阻效果對于翼舵界面處熱傳導影響效果甚微，而對于非氣動熱外表面的接觸界面，界面接觸熱阻會占據主導作用，極大地影響結構溫度場和熱防護效果。

3.4.2 結構位移場和應力場

在同時考慮氣動壓力、氣動熱以及界面接觸熱阻的情況下，利用前述有限元方法對翼舵結構進行熱強度計算，得到的總體位移場和等效應力場如圖15、圖16所示。

圖15 翼舵結構總位移場Fig. 15 The total displacement field of the rudder

圖16 翼舵結構等效應力場Fig. 16 The equivalent stress field of the rudder

從計算結果可以看出，如果同時考慮氣動力與氣動熱的作用，由于GH99熱膨脹系數較大，結構位移的變化主要由受熱變形的主導，氣動力對結構應力的變化影響較小。在翼舵與底盤裝配區域等效應力是最大的，最大值可以達到600 MPa，在界面處，由于涉及到不同結構的裝配，因此接觸應力相對較高，接觸熱阻較小，因此界面對結構溫度場的影響較小，界面接觸熱阻可忽略不計，無明顯溫度跳躍現象，界面兩側的應力場分布連續，整體結構仍處于材料的強度極限范圍內。一般來講，結構溫度場對結構應力的分布影響較大，計算結果也表明了彈塑性變形對結構強度分析的重要性。因此，通過調控界面接觸熱阻去實現結構溫度場的重分布，從而影響結構的等效應力場，是一種提高結構安全性的有效手段。本算例充分展示了本文發展的分區耦合連續/間斷伽遼金有限元方法和及編制的計算軟件在處理氣動-熱-力耦合分析問題上的可行性。

4 結 論

在國家數值風洞工程項目的支持下，采用松耦合的方法建立了氣動熱力耦合問題的連續/間斷伽遼金有限元計算格式，充分考慮搭接結構在界面處廣泛存在的非完好接觸現象，引入溫度和應力相關的界面接觸熱阻，在界面處使用間斷伽遼金有限元方法，在連續區域采用經典的連續伽遼金有限元方法，極大提高了計算精度和計算效率，并編制了具有自主知識產權的氣動-熱-力耦合計算分析軟件。數值算例驗證了本文建立的方法和編制的計算軟件的精度和可靠性，形成了大規模工程熱力耦合問題的計算分析能力，為我國新型飛行器熱防護結構的研發提供重要的技術支撐。下一步將建立瞬態溫度場和結構動力學的連續/間斷伽遼金有限元方法，將計算能力從靜態擴展到動態分析，實現飛行全過程的實時高精度數值仿真。