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

(1. 北京科技大學(xué) 數(shù)理學(xué)院應(yīng)用力學(xué)系，北京 100083；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽 621000)

0 引 言

近空間高超聲速飛行器是當(dāng)前國(guó)內(nèi)外的研究熱點(diǎn)之一，其發(fā)展水平也是一個(gè)國(guó)家綜合國(guó)力的重要體現(xiàn)[1]。近空間高超聲速飛行器的研發(fā)涉及到氣動(dòng)布局、氣動(dòng)加熱環(huán)境模擬、結(jié)構(gòu)熱防護(hù)、發(fā)動(dòng)機(jī)研制、航空航天材料、高能燃料、全球定位、矢量控制等多學(xué)科技術(shù)問題，其中熱防護(hù)結(jié)構(gòu)是起到?jīng)Q定性作用的一項(xiàng)關(guān)鍵技術(shù)[2-4]。對(duì)熱防護(hù)結(jié)構(gòu)的設(shè)計(jì)及安全評(píng)估涉及到氣動(dòng)熱力耦合分析，即高超聲速飛行引起的氣動(dòng)力和氣動(dòng)熱對(duì)飛行器結(jié)構(gòu)的傳熱防熱及熱力耦合行為的影響[5-6]，以及由此引起的熱氣彈問題[7-8]。

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

特別的，為了提高結(jié)構(gòu)熱防護(hù)的效率，不同的飛行器結(jié)構(gòu)部位采用的防熱結(jié)構(gòu)和使用的防熱材料體系也往往不同，因此產(chǎn)生了大量的異種材料裝配界面。裝配界面會(huì)導(dǎo)致兩方面的問題，一是非完好接觸的裝配界面會(huì)產(chǎn)生阻礙熱量流動(dòng)的接觸熱阻，從而對(duì)結(jié)構(gòu)溫度場(chǎng)產(chǎn)生影響，二是由于界面兩側(cè)不同材料的熱物性參數(shù)可能相差很大，由此帶來了嚴(yán)重的熱失配問題。由于界面接觸熱阻受界面溫度和界面接觸應(yīng)力的影響，因此形成了復(fù)雜的耦合關(guān)系，如圖1所示。

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

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

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

1 有限元格式及計(jì)算流程

1.1 氣動(dòng)節(jié)點(diǎn)數(shù)據(jù)向有限元節(jié)點(diǎn)的轉(zhuǎn)換方法

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

1.2 溫度場(chǎng)分析的連續(xù)/間斷有限元格式

考慮邊界為 ? Ω 的區(qū)域 Ω上的橢圓型熱傳導(dǎo)方程：

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

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

等效積分弱形式的熱傳導(dǎo)方程可寫成：

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

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

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

其中，熱流數(shù)值通量為：

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

同時(shí)必須在邊界上加入穩(wěn)定項(xiàng)，這里將單元內(nèi)部邊界基本場(chǎng)變量的跳變引入到平衡方程中，稱為穩(wěn)定項(xiàng)，此時(shí)單元間斷伽遼金有限元方程可寫成：

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

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

同時(shí)為了體現(xiàn)這種強(qiáng)的間斷效應(yīng)，在存在接觸熱阻的單元內(nèi)邊界上要將穩(wěn)定化參數(shù)賦零。

假設(shè)單元溫度場(chǎng)可以表示為：

這里，NTe為溫度場(chǎng)插值函數(shù)，Te為單元節(jié)點(diǎn)的溫度列向量。權(quán)函數(shù)取為ve=(NTe)T，代入單元有限元方程，并引入氣動(dòng)熱載荷和輻射邊界條件，經(jīng)過推導(dǎo)并組裝可得用溫度節(jié)點(diǎn)列向量為唯一未知量的結(jié)構(gòu)總體傳熱有限元方程組：

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

1.3 位移場(chǎng)分析的連續(xù)/間斷有限元格式

結(jié)構(gòu)位移場(chǎng)分析的平衡方程可寫成：

其中，σ為應(yīng)力向量，f為體力向量，▽為微分算子。

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

其中，we為 定義在區(qū)域 Ω上矢量形式的權(quán)函數(shù)。

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

通過將單元邊界應(yīng)力 σe用應(yīng)力數(shù)值通量進(jìn)行替換即可得到單元應(yīng)力分析的間斷伽遼金有限元方程：

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

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

假設(shè)單元e的位移場(chǎng)和應(yīng)變場(chǎng)通過插值可以表示為：

其中，Nue為 單元位移場(chǎng)的插值函數(shù)矩陣，Ue為單元節(jié)點(diǎn)位移列向量，Be=▽TNue為幾何矩陣。

進(jìn)行應(yīng)力場(chǎng)分析時(shí)，需要考慮溫度產(chǎn)生的變形，此時(shí)的本構(gòu)方程為：

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

引入強(qiáng)制位移邊界條件即可求解整體有限元方程組，得到結(jié)構(gòu)位移場(chǎng)，進(jìn)一步經(jīng)過應(yīng)力磨平等處理即可得到高精度的應(yīng)力場(chǎng)。

值得注意的是，本文假設(shè)結(jié)構(gòu)的變形處于小變形狀態(tài)，因此熱應(yīng)力對(duì)結(jié)構(gòu)的影響僅體現(xiàn)在等效載荷列向量中，并未考慮其對(duì)結(jié)構(gòu)幾何剛度的影響。

1.4 計(jì)算流程及軟件開發(fā)

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

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

首先進(jìn)行流固界面載荷轉(zhuǎn)換，輸入氣動(dòng)熱力計(jì)算結(jié)果和結(jié)構(gòu)有限元網(wǎng)格信息，利用界面載荷轉(zhuǎn)換程序得到有限元節(jié)點(diǎn)或者單元上的氣動(dòng)熱流和氣動(dòng)力信息。

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

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

基于上述有限元格式和計(jì)算流程圖編制了相應(yīng)的計(jì)算程序，開發(fā)了相配套的軟件，并作為一個(gè)功能模塊集成到國(guó)家數(shù)值風(fēng)洞軟件群中。

2 算法驗(yàn)證

2.1 問題描述

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

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

2.2 結(jié)果對(duì)比

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

表1 組合桿件中點(diǎn)位置的等效應(yīng)力Table 1 Equivalent stress at the center of the composite rod

圖4 組合桿件溫度場(chǎng)的對(duì)比Fig. 4 Comparisons of the temperature field of the composite rod

圖6 組合桿件等效應(yīng)力場(chǎng)的對(duì)比Fig. 6 Comparisons of the equivalent stress of the composite rod

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

從以上結(jié)果對(duì)比可以看出，接觸熱阻的存在導(dǎo)致溫度場(chǎng)分布在界面處產(chǎn)生跳變，接觸熱阻越大，溫度的跳變?cè)酱蟆Ｑ剌S線路徑上的溫度場(chǎng)和位移場(chǎng)以及桿件中點(diǎn)處的等效應(yīng)力數(shù)值，本文程序結(jié)果和理論解及通用軟件結(jié)果一致，桿件中點(diǎn)處的等效應(yīng)力的相對(duì)誤差不超過0.4%。除此算例外，本文有限元程序在設(shè)計(jì)開發(fā)中，基于軟件工程的思想，對(duì)所有的邊界材料、材料屬性、結(jié)構(gòu)形式、載荷形式進(jìn)行了全方位的測(cè)試，本程序結(jié)果均與通用有限元軟件結(jié)果或者理論解完全吻合，這充分說明了本文程序的精度和可靠性，可進(jìn)一步用于實(shí)際工程結(jié)構(gòu)的數(shù)值仿真。

圖5 組合桿件總位移場(chǎng)的對(duì)比Fig. 5 Comparisons of the total displacement of the composite rod

3 工程算例

3.1 問題描述

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

圖8 舵面模型結(jié)構(gòu)圖Fig. 8 The geometry of the rudder model

3.2 流固界面氣動(dòng)熱力轉(zhuǎn)換

基于1.1節(jié)的轉(zhuǎn)換方法得到的轉(zhuǎn)換結(jié)果如圖9～圖11所示。結(jié)構(gòu)溫度場(chǎng)和應(yīng)力場(chǎng)分析時(shí)需要給定氣動(dòng)熱流、氣動(dòng)力邊界條件，相關(guān)邊界條件由CFD計(jì)算給出，提供的環(huán)境數(shù)據(jù)為翼面上各個(gè)氣動(dòng)節(jié)點(diǎn)的氣動(dòng)力、冷壁熱流和恢復(fù)焓。由于氣動(dòng)網(wǎng)格節(jié)點(diǎn)和有限元分析網(wǎng)格的節(jié)點(diǎn)是不重合的，因此首先需要將CFD氣動(dòng)點(diǎn)的氣動(dòng)力、冷壁熱流以及恢復(fù)焓轉(zhuǎn)換到有限元分析時(shí)的單元結(jié)構(gòu)節(jié)點(diǎn)上。

圖9 氣動(dòng)和結(jié)構(gòu)分析時(shí)計(jì)算網(wǎng)格的對(duì)比Fig. 9 Comparisons of meshes for CFD and FEA computations

圖10 轉(zhuǎn)換前后氣動(dòng)和結(jié)構(gòu)分析氣動(dòng)壓強(qiáng)的對(duì)比Fig. 10 Comparisons of aerodynamic pressures for the CFD and FEA computations after conversion

圖11 轉(zhuǎn)換前后氣動(dòng)和結(jié)構(gòu)分析氣動(dòng)熱流的對(duì)比Fig. 11 Comparisons of aerodynamic heat fluxes for the CFD and FEA computations after conversion

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

由于本算例主要針對(duì)面分布的熱流和氣動(dòng)壓強(qiáng)進(jìn)行轉(zhuǎn)換，從快速確定有限元節(jié)點(diǎn)場(chǎng)量的角度出發(fā)，采用簡(jiǎn)化的轉(zhuǎn)換算法，將距離有限元節(jié)點(diǎn)最近的氣動(dòng)點(diǎn)的場(chǎng)量值作為該有限元節(jié)點(diǎn)的場(chǎng)量值。

從轉(zhuǎn)換前后氣動(dòng)和有限元的場(chǎng)量云圖可以看出，在大部分區(qū)域的轉(zhuǎn)換效果都很好，而在翼前緣區(qū)域由于轉(zhuǎn)換前后網(wǎng)格差異較大，導(dǎo)致轉(zhuǎn)換誤差稍微大一些，同時(shí)由于該區(qū)域范圍較小，因此轉(zhuǎn)換誤差對(duì)結(jié)構(gòu)變形和應(yīng)力的影響是可以忽略的。

3.3 邊界條件以及熱阻模型

進(jìn)行熱分析時(shí)對(duì)翼舵外表面施加給定熱流邊界條件，根據(jù)CFD計(jì)算得到的冷壁熱流qc和恢復(fù)焓hr計(jì)算得到施加在結(jié)構(gòu)上的熱壁熱流qn為：

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

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

溫度場(chǎng)分析時(shí)翼舵底部安裝面給定溫度280 K，輻射環(huán)境參考溫度280 K，由于給定的熱流邊界條件與待求溫度T相關(guān)，因此這是一個(gè)非線性溫度場(chǎng)問題，依據(jù)上述程序進(jìn)行迭代求解。

在對(duì)翼舵強(qiáng)度分析時(shí)，外表面施加給定的氣動(dòng)力載荷，翼舵安裝面位移固定，熱變形參考溫度280 K，而界面接觸熱阻與界面應(yīng)力密切相關(guān)[32-33]：

其中， σc為界面接觸應(yīng)力，MPa。當(dāng)應(yīng)力較小時(shí)，界面接觸熱阻較大，隨界面應(yīng)力的增大，界面接觸熱阻隨應(yīng)力指數(shù)型減小。

3.4 計(jì)算結(jié)果

3.4.1 結(jié)構(gòu)溫度場(chǎng)

基于上述有限元方法，在不考慮界面接觸熱阻時(shí)得到的結(jié)構(gòu)溫度場(chǎng)如圖12所示。

圖12 無接觸熱阻時(shí)舵面溫度場(chǎng)Fig. 12 The temperature field of the rudder without the thermal contact resistance

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

不同接觸熱阻時(shí)舵面溫度場(chǎng)如圖13所示，翼舵結(jié)構(gòu)界面兩側(cè)節(jié)點(diǎn)溫度隨接觸熱阻的變化趨勢(shì)如圖14所示。

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

圖14 界面接觸熱阻對(duì)接觸界面兩側(cè)節(jié)點(diǎn)溫度的影響Fig. 14 Variations of interface temperatures with the interface thermal contact resistance

界面接觸熱阻在組合結(jié)構(gòu)熱傳導(dǎo)過程中起到阻礙熱量傳遞的作用，會(huì)造成界面兩側(cè)溫度的跳變，從本翼舵算例結(jié)果可以看出，界面接觸熱阻越大界面兩側(cè)溫度跳變?cè)矫黠@，但對(duì)于整體溫度場(chǎng)分布的影響不明顯，這是因?yàn)楸舅憷袩崃鬟吔鐥l件與輻射邊界同時(shí)施加在翼舵外表面上，氣動(dòng)熱和輻射相互作用的影響對(duì)于翼舵表面溫度場(chǎng)的分布占據(jù)主導(dǎo)作用，因此界面接觸熱阻效果對(duì)于翼舵界面處熱傳導(dǎo)影響效果甚微，而對(duì)于非氣動(dòng)熱外表面的接觸界面，界面接觸熱阻會(huì)占據(jù)主導(dǎo)作用，極大地影響結(jié)構(gòu)溫度場(chǎng)和熱防護(hù)效果。

3.4.2 結(jié)構(gòu)位移場(chǎng)和應(yīng)力場(chǎng)

在同時(shí)考慮氣動(dòng)壓力、氣動(dòng)熱以及界面接觸熱阻的情況下，利用前述有限元方法對(duì)翼舵結(jié)構(gòu)進(jìn)行熱強(qiáng)度計(jì)算，得到的總體位移場(chǎng)和等效應(yīng)力場(chǎng)如圖15、圖16所示。

圖15 翼舵結(jié)構(gòu)總位移場(chǎng)Fig. 15 The total displacement field of the rudder

圖16 翼舵結(jié)構(gòu)等效應(yīng)力場(chǎng)Fig. 16 The equivalent stress field of the rudder

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

4 結(jié) 論

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