翼面熱靜氣動彈性的流固熱交錯迭代耦合分析*

常斌，黃杰，2，姚衛(wèi)星，2

（1.南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室 南京，210016）

（2.南京航空航天大學(xué)飛行器先進(jìn)設(shè)計(jì)技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室 南京，210016）

引言

翼面的靜氣動彈性是飛行器設(shè)計(jì)必須考慮的問題，它涉及到氣動力與結(jié)構(gòu)彈性變形之間的耦合。隨著飛行器速度越來越快，氣動加熱引起翼面結(jié)構(gòu)溫度升高［1-3］，導(dǎo)致翼面結(jié)構(gòu)剛度發(fā)生變化［4-5］，并由此提出了熱靜氣動彈性的問題，熱靜氣動彈性的研究對高超聲速飛行器的設(shè)計(jì)至關(guān)重要。

最早的高超聲速翼面熱靜氣動彈性分析方法假設(shè)翼面結(jié)構(gòu)具有均勻的溫度場分布，分析不同溫度下的翼面熱剛度，并在此熱剛度下進(jìn)行翼面的靜氣動彈性分析。隨著氣動熱分析技術(shù)研究的深入，一些學(xué)者采用以Eckert參考焓法［6］為代表的工程近似方法分析翼面氣動熱問題，將分析獲得的翼面熱流作為邊界條件，并通過有限差分法等求解熱傳導(dǎo)方程獲得翼面溫度場分布，再研究翼面的熱剛度和熱靜氣動彈性［7-8］。以上傳統(tǒng)翼面熱靜氣動彈性分析方法的計(jì)算精度較差，已無法滿足現(xiàn)代及未來高超聲速飛行器熱靜氣動彈性的分析要求。傳統(tǒng)方法最主要的缺點(diǎn)是未考慮結(jié)構(gòu)溫度場對氣動熱的反饋效應(yīng)。實(shí)際上氣動熱與結(jié)構(gòu)溫度場之間存在著強(qiáng)烈的耦合效應(yīng)，但傳統(tǒng)分析方法僅考慮了氣動熱引起的結(jié)構(gòu)溫度場升高，而未考慮翼面溫度升高后對壁面熱流的影響。

隨著計(jì)算流體力學(xué)（computational fluid dynamics，簡稱CFD）、數(shù)值傳熱學(xué)（numerical heat transfer，簡稱NHT）及計(jì)算結(jié)構(gòu)動力學(xué)（computational structural dynamics，簡稱CSD）的發(fā)展，近年來一些學(xué)者通過CFD和NHT之間的耦合分析求解結(jié)構(gòu)溫度場，研究熱環(huán)境對結(jié)構(gòu)熱剛度的影響，最后通過CFD和CSD之間的耦合分析求解熱靜氣動彈性響應(yīng)［9-10］。該方法分2步計(jì)算，能獲得精確的壁面熱流和結(jié)構(gòu)溫度場，并準(zhǔn)確評價(jià)翼面結(jié)構(gòu)熱剛度，具有計(jì)算精度高的特點(diǎn)。由于CFD與NHT的耦合以及CFD與CSD的耦合需要進(jìn)行2次CFD計(jì)算，故該方法還存在計(jì)算效率較低的缺點(diǎn)。此外該方法也無法考慮翼面變形對氣動熱帶來的影響，因此需要一種兼顧計(jì)算精度和計(jì)算效率的耦合分析方法。

筆者提出了一種針對高超聲速翼面熱靜氣動彈性的流固熱交錯迭代耦合方法，其主要特點(diǎn)是能考慮翼面變形對氣動熱的影響，且不需要像傳統(tǒng)熱靜氣動彈性耦合方法分2步求解，并且利用該方法進(jìn)行了典型高超聲速翼面熱靜氣動彈性問題的分析。

1 流體和結(jié)構(gòu)傳熱控制方程

流體動力學(xué)的控制方程為Navier-Stokes方程，

其在直角坐標(biāo)系下的積分形式為

其中：W為守恒向量；Fc為 對 流 通量；Fv為黏性通量；dS為控制體V的邊界面；n為邊界面dS的 外法線單位向量。

Navier-Stokes方程中對流通量Fc的空間離散采用AUSM+格式［11］，黏性通量Fv的空間離散采用中心差分格 式，湍 流模擬采用Menter′s SSTk-ω兩方程模型［12］，時間推進(jìn)采用LU-SGS格式［13］，且采用聚合多重網(wǎng)格法［14］加速收斂。

翼面結(jié)構(gòu)傳熱分析的控制方程為

其中：kx，ky和kz分別為材料在x，y和z方向的導(dǎo)熱系數(shù)；ρ和c分別為結(jié)構(gòu)材料的密度和比熱容。

翼面結(jié)構(gòu)傳熱分析的邊界條件包括壁面熱流Qw和壁面輻射熱流Qr，其表達(dá)式為

其中：ε為表面發(fā)射率；σ為玻爾茲曼常數(shù)；Tw為壁面溫度。

筆者采用有限元法進(jìn)行翼面結(jié)構(gòu)傳熱及變形的分析。

2 交錯迭代耦合分析方法

高超聲速飛行器做定常飛行時，翼面所受的氣動力和氣動熱均穩(wěn)定，其可通過數(shù)值方法計(jì)算同時得到，因此統(tǒng)稱為氣動環(huán)境。翼面熱靜氣動彈性分析中主要存在2種耦合效應(yīng)：①氣動環(huán)境中的氣動熱與結(jié)構(gòu)傳熱之間的耦合效應(yīng)，當(dāng)氣動熱作用于翼面時，翼面結(jié)構(gòu)溫度升高，翼面結(jié)構(gòu)溫度的升高導(dǎo)致邊界層內(nèi)氣體與壁面溫度的梯度減小，即氣動熱效應(yīng)將減弱；②氣動環(huán)境中的氣動力與結(jié)構(gòu)變形之間的耦合效應(yīng)，當(dāng)氣動力作用于翼面時，翼面產(chǎn)生變形，變形會造成氣動力重分布。此外，翼面變形也會造成壁面熱流的重分布，且熱環(huán)境還會影響翼面結(jié)構(gòu)剛度（材料剛度和熱應(yīng)力引起的附加幾何剛度），進(jìn)而影響結(jié)構(gòu)變形。筆者考慮以上所述的耦合效應(yīng)，而忽略其他次要因素，建立了如圖1所示的翼面熱靜氣動彈性分析耦合模型。

圖1 耦合模型Fig.1 Coupled model

根據(jù)圖1中的耦合模型，筆者采用流固熱交錯迭代耦合方法進(jìn)行高超聲速翼面熱靜氣動彈性的研究，交錯迭代耦合方法如圖2所示，圖中的i為耦合迭代步。采用CFD方法進(jìn)行翼面氣動力和氣動熱的計(jì)算，采用NHT方法進(jìn)行翼面結(jié)構(gòu)溫度場的計(jì)算，采用CSD方法進(jìn)行翼面結(jié)構(gòu)變形的計(jì)算，其中CFD方法基于有限體積法，而NHT和CSD方法均基于有限元法。流固熱交錯迭代耦合方法的基本假設(shè)和特點(diǎn)如下：

圖2 交錯迭代耦合方法Fig.2 Staggered iterative coupled method

1）在流場計(jì)算過程中（氣動力和氣動熱）翼面結(jié)構(gòu)參數(shù)（結(jié)構(gòu)變形和壁面溫度）保持不變，同理在結(jié)構(gòu)計(jì)算過程中（結(jié)構(gòu)變形和壁面溫度）翼面氣動環(huán)境（氣動力和氣動熱）保持不變；

2）只需進(jìn)行一次耦合分析即可獲得翼面氣動力、壁面熱流、溫度場和結(jié)構(gòu)變形的全部穩(wěn)態(tài)結(jié)果。

翼面熱靜氣動彈性的交錯迭代耦合方法分析流程如圖3所示，其主要步驟為：

圖3 交錯迭代耦合分析流程Fig.3 Staggered iterative coupled analysis process

1）首先建立CFD，NHT和CSD數(shù)值分析模型，通過CFD分析獲得初始?xì)鈩恿捅诿鏌崃鳎?/p>

2）進(jìn)行第i步的CFD計(jì)算，將計(jì)算獲得的氣動力Fi和壁面熱流Qi分別傳遞給CSD模型和NHT模型；

3）進(jìn)行第i步的NHT計(jì)算，將計(jì)算獲得的翼面結(jié)構(gòu)溫度場Ti傳遞給CSD模型，并更新翼面結(jié)構(gòu)剛度；

4）進(jìn)行第i步的CSD計(jì)算，將獲得的翼面表面位移Uwi及表面溫度Twi傳遞給CFD模型，并進(jìn)行流體網(wǎng)格變形；

5）若翼面結(jié)構(gòu)變形和溫度場均收斂，則結(jié)束迭代。否則重復(fù)步驟2～4，直至結(jié)果收斂，結(jié)束分析。

由于流場和結(jié)構(gòu)的網(wǎng)格尺寸相差較大，耦合變量需要在耦合面上通過插值算法實(shí)現(xiàn)數(shù)據(jù)的傳遞。筆者采用文獻(xiàn)［15］中的基于控制面的雙向映射插值方法進(jìn)行壁面熱流、壁面溫度、氣動力以及翼面變形的數(shù)據(jù)傳遞。此外為了加速收斂，在流場計(jì)算中引入了子迭代。

目前還沒有熱靜氣動彈性的標(biāo)準(zhǔn)試驗(yàn)，且翼面熱靜氣動彈性分析的核心是精確計(jì)算氣動熱和結(jié)構(gòu)溫度場，故筆者采用NASA的高超聲速空心圓管模型［16］進(jìn)行氣動熱和結(jié)構(gòu)溫度場計(jì)算精度的驗(yàn)證。文獻(xiàn)［17］已經(jīng)進(jìn)行了空心圓管模型的耦合計(jì)算，外壁面熱流和溫度分布的計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好，驗(yàn)證了計(jì)算精度。

3 翼面熱靜氣動彈性分析

3.1 計(jì)算模型

筆者采用小展弦比翼面進(jìn)行熱靜氣動彈性的分析，其平面和剖面如圖4所示。小展弦比翼面來源于F104戰(zhàn)斗機(jī)機(jī)翼，圖中參數(shù)C為翼根弦長。假設(shè)高超聲速自由來流馬赫數(shù)從6增加至9，飛行高度為10 km，攻角為0.5°，且整個翼面結(jié)構(gòu)的初始溫度為300 K。劃分了流場的結(jié)構(gòu)網(wǎng)格，網(wǎng)格量約為80萬，為了保證翼面熱流的計(jì)算精度，壁面第1層網(wǎng)格高度小于1×10-5m，圖5為CFD計(jì)算網(wǎng)格及結(jié)構(gòu)傳熱和變形分析的翼面結(jié)構(gòu)表面計(jì)算網(wǎng)格。筆者假設(shè)翼面前緣為耐高溫碳/碳復(fù)合材料結(jié)構(gòu)，其導(dǎo)熱系數(shù)、彈性模量和熱膨脹系數(shù)分別為42 W/（m·K），95 GPa和4×10-6K-1。翼面其他部位采用TA7耐高溫鈦合金結(jié)構(gòu)，材料密度為4.4×103kg/m3，但材料導(dǎo)熱系數(shù)、彈性模量和熱膨脹系數(shù)均與溫度相關(guān)，可參考文獻(xiàn)［15］。文獻(xiàn)［15］進(jìn)行了翼面的CFD和NHT耦合計(jì)算，獲得了結(jié)構(gòu)溫度場，并研究了在熱環(huán)境下彈性模量和熱應(yīng)力對模態(tài)的影響。筆者考慮了氣動熱、氣動力、結(jié)構(gòu)傳熱和結(jié)構(gòu)變形，是典型的熱靜氣動彈性問題。翼面外表面有高輻射率的涂層，其發(fā)射率為0.85。此外流場計(jì)算中引入了子迭代，子迭代步數(shù)設(shè)置為30。

圖4 小展弦比翼面的平面和剖面Fig.4 Platform and cross-sectional views of the low aspect ratio wing

圖5 計(jì)算網(wǎng)格Fig.5 Computational grids

3.2 翼面熱環(huán)境及熱靜氣動彈性變形

筆者通過流固熱交錯迭代耦合方法獲得了馬赫數(shù)從6增加至9的計(jì)算結(jié)果。翼面結(jié)構(gòu)溫度場和翼面變形均在25步以內(nèi)達(dá)到收斂狀態(tài)。圖6給出了其表面溫度云圖，結(jié)果表明在不同馬赫數(shù)下計(jì)算獲得的翼面結(jié)構(gòu)溫度分布情況相似。由于翼面前緣當(dāng)量半徑很小，故該部位的壁面熱流和溫度最高，且往下游翼面結(jié)構(gòu)溫度逐漸降低。此外馬赫數(shù)從6增加到9，翼面結(jié)構(gòu)最高溫度從1 525.4 K上升到2 163.3 K，翼面結(jié)構(gòu)最低溫度從824.9 K上升到996.2 K。

圖6 不同馬赫數(shù)下翼面結(jié)構(gòu)溫度場Fig.6 Temperature field of wing for different Mach number

熱靜氣動彈性與常規(guī)靜氣動彈性的主要區(qū)別是熱環(huán)境下翼面結(jié)構(gòu)剛度會發(fā)生變化。首先，高溫會直接影響材料的彈性模量，針對本研究采用的TA7耐高溫鈦合金材料，溫度越高材料的彈性模量越低；其次，高溫下翼面結(jié)構(gòu)會膨脹，結(jié)構(gòu)受到約束及內(nèi)部溫度分布不均勻時均會產(chǎn)生熱應(yīng)力。熱應(yīng)力的本質(zhì)是預(yù)應(yīng)力，預(yù)應(yīng)力會造成幾何非線性現(xiàn)象，即翼面結(jié)構(gòu)會產(chǎn)生附加的幾何剛度。因此，必須綜合考慮熱環(huán)境造成的材料彈性模量的降低及熱應(yīng)力引起的附加幾何剛度才能準(zhǔn)確計(jì)算翼面變形。

筆者通過流固熱交錯迭代耦合方法分析獲得了馬赫數(shù)從6增加至9的翼面法向（垂直于翼面方向）的熱靜氣動彈性變形，其熱環(huán)境下翼面法向位移如圖7所示，結(jié)果表明在不同馬赫數(shù)下翼面熱靜氣動彈性變形情況相似。圖8為300 K和熱環(huán)境下翼面最大厚度處的法向位移沿展向的分布情況，圖9為300 K和熱環(huán)境下翼面法向最大位移隨馬赫數(shù)的變化情況。

圖7 熱環(huán)境下翼面法向位移Fig.7 Normal displacements of wing under thermal environments

圖8 翼面法向位移沿展向分布情況Fig.8 Normal displacements of wing along the spanwise direction

圖9 翼面最大法向位移Fig.9 Maximum normal displacements of wing

由以上分析結(jié)果可知，隨著馬赫數(shù)的增加，翼面變形越大。此外在相同馬赫數(shù)下，考慮熱效應(yīng)時翼面法向位移明顯大于300 K時的分析結(jié)果。馬赫數(shù)為6，7，8和9時，熱環(huán)境下翼面法向最大位移分別比300 K時大41.05%，51.10%，57.64%和63.28%，即馬赫數(shù)越大，考慮熱環(huán)境和不考慮熱環(huán)境獲得的翼面變形量之差越大。這是由于熱環(huán)境下翼面結(jié)構(gòu)剛度降低了，且馬赫數(shù)越大翼面結(jié)構(gòu)剛度下降越快，造成變形量迅速增加。

綜上所述，進(jìn)行高超聲速翼面靜氣動彈性分析時必須考慮其熱效應(yīng)，熱環(huán)境導(dǎo)致翼面結(jié)構(gòu)剛度降低，造成翼面變形高于常溫下的結(jié)果，應(yīng)給予足夠重視。

4 結(jié)論

1）提出了一種針對高超聲速翼面熱靜氣動彈性的流固熱交錯迭代數(shù)值耦合方法。該方法充分考慮了氣動環(huán)境（氣動力和氣動熱）與結(jié)構(gòu)變形之間的耦合、氣動熱與結(jié)構(gòu)溫度場之間的耦合以及溫度場對結(jié)構(gòu)剛度的影響，能克服傳統(tǒng)熱靜氣動彈性耦合方法中未考慮翼面變形對氣動熱影響的缺點(diǎn)。

2）進(jìn)行了翼面熱靜氣動彈性分析，針對本研究的翼面結(jié)構(gòu)形式和來流狀態(tài)，熱環(huán)境造成了翼面結(jié)構(gòu)剛度下降，從而導(dǎo)致熱環(huán)境下翼面法向最大位移比300 K時增加了40%以上，且馬赫數(shù)越大，兩者之差越大。因此，熱環(huán)境對高超聲速翼面的靜氣動彈性變形影響明顯，在分析時必須考慮其影響才能獲得準(zhǔn)確的分析結(jié)果。