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

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

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

引言

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

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

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

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

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

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

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

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

Navier-Stokes方程中對(duì)流通量Fc的空間離散采用AUSM+格式［11］，黏性通量Fv的空間離散采用中心差分格 式，湍 流模擬采用Menter′s SSTk-ω兩方程模型［12］，時(shí)間推進(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 交錯(cuò)迭代耦合分析方法

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

圖1 耦合模型Fig.1 Coupled model

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

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

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

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

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

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

1）首先建立CFD，NHT和CSD數(shù)值分析模型，通過(guò)CFD分析獲得初始?xì)鈩?dòng)力和壁面熱流；

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

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

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

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

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

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

3 翼面熱靜氣動(dòng)彈性分析

3.1 計(jì)算模型

筆者采用小展弦比翼面進(jìn)行熱靜氣動(dòng)彈性的分析，其平面和剖面如圖4所示。小展弦比翼面來(lái)源于F104戰(zhàn)斗機(jī)機(jī)翼，圖中參數(shù)C為翼根弦長(zhǎng)。假設(shè)高超聲速自由來(lái)流馬赫數(shù)從6增加至9，飛行高度為10 km，攻角為0.5°，且整個(gè)翼面結(jié)構(gòu)的初始溫度為300 K。劃分了流場(chǎng)的結(jié)構(gòu)網(wǎng)格，網(wǎng)格量約為80萬(wàn)，為了保證翼面熱流的計(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)溫度場(chǎng)，并研究了在熱環(huán)境下彈性模量和熱應(yīng)力對(duì)模態(tài)的影響。筆者考慮了氣動(dòng)熱、氣動(dòng)力、結(jié)構(gòu)傳熱和結(jié)構(gòu)變形，是典型的熱靜氣動(dòng)彈性問(wèn)題。翼面外表面有高輻射率的涂層，其發(fā)射率為0.85。此外流場(chǎng)計(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)境及熱靜氣動(dòng)彈性變形

筆者通過(guò)流固熱交錯(cuò)迭代耦合方法獲得了馬赫數(shù)從6增加至9的計(jì)算結(jié)果。翼面結(jié)構(gòu)溫度場(chǎng)和翼面變形均在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)溫度場(chǎng)Fig.6 Temperature field of wing for different Mach number

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

筆者通過(guò)流固熱交錯(cuò)迭代耦合方法分析獲得了馬赫數(shù)從6增加至9的翼面法向（垂直于翼面方向）的熱靜氣動(dòng)彈性變形，其熱環(huán)境下翼面法向位移如圖7所示，結(jié)果表明在不同馬赫數(shù)下翼面熱靜氣動(dòng)彈性變形情況相似。圖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)時(shí)翼面法向位移明顯大于300 K時(shí)的分析結(jié)果。馬赫數(shù)為6，7，8和9時(shí)，熱環(huán)境下翼面法向最大位移分別比300 K時(shí)大41.05%，51.10%，57.64%和63.28%，即馬赫數(shù)越大，考慮熱環(huán)境和不考慮熱環(huán)境獲得的翼面變形量之差越大。這是由于熱環(huán)境下翼面結(jié)構(gòu)剛度降低了，且馬赫數(shù)越大翼面結(jié)構(gòu)剛度下降越快，造成變形量迅速增加。

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

4 結(jié)論

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

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