李佳偉，王江峰，楊天鵬，李龍飛，王丁

南京航空航天大學(xué) 非定?？諝鈩?dòng)力學(xué)與流動(dòng)控制工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室，南京 210016

吸氣式高超聲速飛行器的研究已成為當(dāng)今世界航空航天發(fā)展的一個(gè)重要方向，但同時(shí)，與之伴隨的氣動(dòng)熱防護(hù)問題也成為一項(xiàng)亟待解決的重要技術(shù)難題[1]。如圖1所示，以超燃沖壓發(fā)動(dòng)機(jī)為動(dòng)力的高速飛行器在大氣層中以高馬赫數(shù)飛行時(shí)，將面臨十分嚴(yán)峻的氣動(dòng)熱力學(xué)環(huán)境。特別是，其前體產(chǎn)生的激波系與進(jìn)氣道唇緣弓形激波發(fā)生干擾，形成激波-邊界層、激波-剪切層、激波-激波相互作用等非線性物理現(xiàn)象[2]。激波干擾產(chǎn)生的復(fù)雜流場(chǎng)波系結(jié)構(gòu)，在高速飛行器表面產(chǎn)生局部高溫高壓載荷，給飛行器的性能和結(jié)構(gòu)帶來極大安全隱患，同時(shí)也給高超聲速飛行器熱防護(hù)問題提出了很大挑戰(zhàn)。

Edney[3]根據(jù)激波干擾點(diǎn)的位置，將激波干擾分為6種類型(Ⅰ～Ⅵ型)，其中Ⅳ型干擾形成超聲速射流沖擊壁面，產(chǎn)生極大的壓強(qiáng)和熱流[4-5]，對(duì)飛行器安全影響最大，因此備受關(guān)注，也是本文主要研究?jī)?nèi)容。另外，高速飛行器外部流場(chǎng)的氣動(dòng)加熱與熱防護(hù)結(jié)構(gòu)內(nèi)的熱傳導(dǎo)相互作用不可分割，因此要準(zhǔn)確預(yù)測(cè)“Ⅳ”型激波干擾下飛行器的氣動(dòng)熱環(huán)境，必須考慮氣動(dòng)加熱與結(jié)構(gòu)傳熱的相互耦合作用[6]。

圖1 前體-進(jìn)氣道激波干擾示意圖
Fig.1 Schematic diagram of forebody-inlet shock wave interference

目前，國(guó)內(nèi)外學(xué)者針對(duì)高超聲速流-固-熱多場(chǎng)耦合問題已經(jīng)開展了大量研究工作。例如，Thornton和Dechaumphai[7]基于有限元方法提出了一種雙向耦合計(jì)算方法來求解不銹鋼平板在高超聲速流動(dòng)下的氣動(dòng)加熱傳熱與變形問題，計(jì)算結(jié)果證明了多場(chǎng)耦合計(jì)算的重要性。Kazemi-Kamyab等[8]在分區(qū)雙向耦合的基礎(chǔ)上發(fā)展了一種強(qiáng)耦合迭代計(jì)算方法，采用高階隱式時(shí)間迭代方法來求解流場(chǎng)與結(jié)構(gòu)非穩(wěn)態(tài)共軛傳熱計(jì)算問題，從而提高時(shí)域上的計(jì)算精度。Chen等[9]發(fā)展了一種高超聲速流場(chǎng)-熱-結(jié)構(gòu)耦合分析平臺(tái)(HyCCD)，分析了高超聲速前緣在高溫高壓條件下的熱力學(xué)行為，同時(shí)引入自適應(yīng)時(shí)間步長(zhǎng)計(jì)算方法，提高計(jì)算效率。Qin等[10]利用MPCCI多平臺(tái)耦合計(jì)算軟件將FLUENT與ABAQUS兩種商用軟件進(jìn)行聯(lián)合使用，計(jì)算模擬了高超聲速條件鈍體外形在不同支桿外形下流-熱-固多場(chǎng)耦合特性，分析了不同外形支桿對(duì)高超聲速鈍頭體的減阻防熱效果。Zhao等[11]依據(jù)結(jié)構(gòu)熱響應(yīng)特征時(shí)間較小的特點(diǎn)，發(fā)展了一種新的多物理場(chǎng)松耦合計(jì)算方法，數(shù)值模擬了馬赫數(shù)為6.47下三維非定常圓管流-熱-固耦合算例，得到的熱壁熱流、冷壁熱流以及圓管溫度與試驗(yàn)值吻合。Zhang等[12]基于自動(dòng)控制理論發(fā)展了一種基于PID(Proportional-Integral-Differential)控制的自適應(yīng)時(shí)間步長(zhǎng)的松耦合高效計(jì)算方法，用于計(jì)算分析高超聲速流中的共軛傳熱問題。

然而需要指出的是，上述計(jì)算方法都?xì)w屬于多場(chǎng)分區(qū)耦合迭代計(jì)算方法，即將流場(chǎng)和固體結(jié)構(gòu)獨(dú)立分開建模，基于“準(zhǔn)穩(wěn)態(tài)耦合”策略，在耦合交界面進(jìn)行數(shù)據(jù)傳遞與交換，同時(shí)在時(shí)間域上進(jìn)行耦合交替迭代。后來，NASA蘭利研究中心的Wieting等[13]認(rèn)為多物理場(chǎng)分區(qū)耦合計(jì)算方法將連續(xù)的物理場(chǎng)進(jìn)行分割建模，耦合界面需要額外的數(shù)據(jù)插值方法和交換策略，這種計(jì)算方法雖然較容易實(shí)現(xiàn)，但無法準(zhǔn)確模擬流-熱-固多物理場(chǎng)耦合特性。李鵬飛和吳頌平[14]在對(duì)高速飛機(jī)機(jī)身結(jié)構(gòu)進(jìn)行氣動(dòng)加熱與結(jié)構(gòu)傳熱的耦合數(shù)值計(jì)算分析后，認(rèn)為為了準(zhǔn)確分析熱壁對(duì)傳熱計(jì)算的影響，很有必要開展氣動(dòng)加熱與結(jié)構(gòu)傳熱一體化數(shù)值計(jì)算方法研究。姜貴慶等[15]為了滿足多場(chǎng)耦合計(jì)算問題中連續(xù)變化的物理?xiàng)l件，提出了基于有限元方法的流-熱-固一體化計(jì)算方法，流場(chǎng)與結(jié)構(gòu)溫度傳熱計(jì)算采用統(tǒng)一網(wǎng)格和計(jì)算程序，計(jì)算結(jié)果表明一體化計(jì)算可以保證熱結(jié)構(gòu)的計(jì)算精度。

本文基于多物理場(chǎng)連續(xù)變化的物理?xiàng)l件，發(fā)展了一種高超聲速流-熱-固一體化計(jì)算方法。該方法基于有限體積法，將高超聲速流場(chǎng)與結(jié)構(gòu)溫度場(chǎng)統(tǒng)一到同一控制方程中，并進(jìn)行統(tǒng)一離散與統(tǒng)一求解，避開了傳統(tǒng)分區(qū)耦合方法所需的數(shù)據(jù)傳遞策略。同時(shí)，提出一種新的雙溫阻模型以保證交界面上物理場(chǎng)的連續(xù)性。選取典型高超聲速繞流二維鈍體算例驗(yàn)證本文方法的可靠性與正確性，最后采用本文計(jì)算方法對(duì)高超聲速“Ⅳ型”激波干擾開展流-熱-固一體化計(jì)算與分析研究。

1 一體化數(shù)值計(jì)算方法

1.1 一體化控制方程及求解方法

首先推導(dǎo)一體化控制方程。固體結(jié)構(gòu)傳熱控制方程積分形式(不含熱源，暫忽略熱輻射換熱)可寫為

(1)

對(duì)于結(jié)構(gòu)外部流場(chǎng)計(jì)算，采用積分形式的可壓縮Navier-Stokes方程作為控制方程，方程形式與式(1)相似，因此，流場(chǎng)計(jì)算方程與固體結(jié)構(gòu)傳熱方程可寫成如下統(tǒng)一控制方程：

(2)

式中：W為守恒量；Fc為對(duì)流通量；Fv為黏性通量。W、Fc、Fv的表達(dá)式為

上述控制方程用于流場(chǎng)計(jì)算時(shí)，u、v、w分別為控制體3個(gè)方向的速度；V=nxu+nyv+nzw為面法向速度，其中nx、ny、nz為面法矢分量；τij為應(yīng)力張量分量；H為總焓，H=E+p/ρ，其中E為流體控制體總能,p為壓力。另外控制方程滿足理想氣體狀態(tài)方程p=ρRT，R為理想氣體常數(shù)。對(duì)于固體結(jié)構(gòu)傳熱計(jì)算時(shí)，結(jié)構(gòu)無變形滿足u=v=w=0，式(2)中對(duì)流通量項(xiàng)Fc=0，黏性通量項(xiàng)中，Θx=k?T/?x，Θy=k?T/?y，Θz=k?T/?z;E=CsT為固體結(jié)構(gòu)單元控制體內(nèi)能。

采用有限體積法對(duì)上述一體化控制方程進(jìn)行離散求解。其中采用AUSM+[16]迎風(fēng)格式求解對(duì)流通量項(xiàng)，同時(shí)基于Venkatakrishnan[17]限制器采用MUSCL[18]插值方法進(jìn)行線性重構(gòu)，以實(shí)現(xiàn)二階精度計(jì)算。采用二階中心差分格式計(jì)算黏性通量項(xiàng)。時(shí)間推進(jìn)采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)隱式時(shí)間迭代計(jì)算方法，其中流場(chǎng)當(dāng)?shù)貢r(shí)間步長(zhǎng)計(jì)算公式為

(3)

式中:Ncfl為CFL數(shù)；Λc和Λv分別為控制體I單元所有邊界的對(duì)流譜半徑和黏性譜半徑的累加；ΩI為第I個(gè)單元的體積；C為常數(shù)。

針對(duì)結(jié)構(gòu)傳熱方程的計(jì)算時(shí)間步長(zhǎng)，需要注意的是，結(jié)構(gòu)傳熱過程中產(chǎn)生的熱脈沖需要限制在當(dāng)?shù)毓腆w結(jié)構(gòu)單元中，才能保證計(jì)算的穩(wěn)定性。因此，結(jié)構(gòu)傳熱方程當(dāng)?shù)貢r(shí)間步長(zhǎng)計(jì)算公式為

(4)

式中：δx為單元最小邊長(zhǎng);α為熱擴(kuò)散率。

另外，為保證壁面熱流計(jì)算的準(zhǔn)確性，需要選取合適的湍流模型。對(duì)多種湍流模型進(jìn)行對(duì)比分析后，確定選取Menter[19]提出的SST(Shear Stress Transport)k-ω兩方程湍流模型，該模型不需要阻尼函數(shù)，能較好計(jì)算邊界層近壁面熱通量。

1.2 流-固耦合界面參數(shù)定義

由于外部流場(chǎng)與固體結(jié)構(gòu)熱力學(xué)特性差異較大，本文一體化計(jì)算方法為了保證流-固耦合界面熱通量與溫度的連續(xù)性，需要對(duì)交界面的熱力學(xué)參數(shù)進(jìn)行特別計(jì)算。

圖2給出了流-固交界面左右單元意圖，Ωf和Ωs分別表示流體單元(左邊)和固體結(jié)構(gòu)單元(右邊)。耦合界面的溫度計(jì)算式為

(5)

式中：下標(biāo)l、r代表耦合界面左、右單元。

(6)

式中:Llr為左右單元中心之間的距離;rlr為單位向量。

值得注意的是，流-固交界面的熱傳導(dǎo)系數(shù)的計(jì)算對(duì)一體化傳熱計(jì)算起到關(guān)鍵作用。固體結(jié)構(gòu)內(nèi)部導(dǎo)熱系數(shù)較大，交界面的熱傳導(dǎo)系數(shù)采用簡(jiǎn)單加權(quán)平均會(huì)對(duì)計(jì)算產(chǎn)生很大誤差，為了提高計(jì)算的準(zhǔn)確性與合理性，本文提出了一種雙-熱阻抗模型用于計(jì)算交界面熱傳導(dǎo)系數(shù)。其中熱阻抗表征的物理意義為熱量在介質(zhì)中傳導(dǎo)時(shí)所遇到的阻力，其定義為

圖2 流-固交界面左右單元意圖
Fig.2 Schematic diagram of left and right cells at fluid-solid interface

(7)

式中:L為熱通量傳導(dǎo)的距離;A為熱量傳導(dǎo)的截面面積。因此，可以得到交界面的熱阻關(guān)系式為

Rt,bnd=Rt,l+Rt,r

(8)

其中：Rt,l和Rt,r分別為交界面左右單元的熱阻。將式(7)代入式(8)中，得到

(9)

因此，可得交界面當(dāng)量導(dǎo)熱系數(shù)k為

(10)

1.3 非定常時(shí)間步長(zhǎng)自適應(yīng)計(jì)算方法

為了提高非定常計(jì)算效率，本文非定常計(jì)算中，真實(shí)物理時(shí)間步長(zhǎng)采用于經(jīng)典PID控制器[22]調(diào)節(jié)的自適應(yīng)控制算法進(jìn)行實(shí)時(shí)計(jì)算，以取代傳統(tǒng)固定時(shí)間步長(zhǎng)，偽時(shí)間步長(zhǎng)選取當(dāng)?shù)貢r(shí)間步長(zhǎng)。

圖3為自適應(yīng)時(shí)間步長(zhǎng)PID控制器示意圖，整個(gè)自適應(yīng)控制器由處理器和控制器2個(gè)部分組成。自適應(yīng)控制基本原理為：處理器將時(shí)間步長(zhǎng)Δt作為輸入控制參數(shù)，進(jìn)行數(shù)值迭代輸出誤差值e；控制器根據(jù)輸出誤差值e和給定輸入的容錯(cuò)值tol進(jìn)行反饋控制，計(jì)算合適的新時(shí)間步長(zhǎng)Δt。

通過PID控制器得到的真實(shí)物理時(shí)間步長(zhǎng)Δt計(jì)算為

(11)

圖3 時(shí)間步長(zhǎng)PID控制器示意圖
Fig.3 Schematic diagram of PID control system in time step-size

式中:n為計(jì)算迭代時(shí)間步；kP、kI、kD為PID控制器參數(shù),分別取值為kP=0.075，kI=0.175,kD=0.015；en為溫度誤差值，其計(jì)算式為

(12)

其中:tol為給定值，本文給定tol=1.0。

為了保證PID控制器的穩(wěn)定性和計(jì)算迭代的收斂性，需要對(duì)時(shí)間步長(zhǎng)的大小和增長(zhǎng)率進(jìn)行一定約束：

Δtmin≤Δtn≤Δtmax

(13)

(14)

(15)

式中：Δtmin和Δtmax分別為設(shè)定的最小和最大的非定常計(jì)算時(shí)間步;αref為時(shí)間步長(zhǎng)控制率，一般0.2<αref<0.4，理想取值αref≈0.2;α的表達(dá)為

(16)

2 驗(yàn)證算例

2.1 計(jì)算模型及網(wǎng)格

選取經(jīng)典不銹鋼圓管前緣氣動(dòng)加熱與傳熱試驗(yàn)[23]來驗(yàn)證一體化計(jì)算方法的準(zhǔn)確性，該算例已多次被用于驗(yàn)證高超聲速流-熱-固多場(chǎng)耦合計(jì)算。圓管外徑尺寸為Rout=38.1 mm,內(nèi)徑尺寸為Rin=25.4 mm,結(jié)構(gòu)材料為不銹鋼，表1給出了不銹鋼材料熱性能參數(shù)[24]。

圖4給出了驗(yàn)證算例的計(jì)算網(wǎng)格和邊界條件。流場(chǎng)與結(jié)構(gòu)傳熱計(jì)算采用統(tǒng)一網(wǎng)格，圖中綠色部分為流場(chǎng)計(jì)算網(wǎng)格，紅色部分為圓管結(jié)構(gòu)傳熱計(jì)算網(wǎng)格。為了更好捕捉流場(chǎng)激波結(jié)構(gòu)，流場(chǎng)網(wǎng)格在激波位置附近進(jìn)行加密。圓管的熱邊界條件為：兩端為絕熱壁，內(nèi)壁為等溫壁，初始溫度為294.4 K。表2給出了初始來流計(jì)算條件，Ma∞、T∞、p∞、Re∞分別為來流馬赫數(shù)、溫度、壓力和單位雷諾數(shù)。

表1 不銹鋼材料熱力學(xué)參數(shù)[24]Table 1 Material thermal properties of stainless steel[24]

圖4 計(jì)算網(wǎng)格與邊界條件
Fig.4 Computational grids and boundary conditions

表2 初始來流條件Table 2 Initial freestream condition

由于近壁面溫度梯度的準(zhǔn)確計(jì)算直接依賴于近壁面網(wǎng)格單元的尺寸，因此為了更加準(zhǔn)確計(jì)算壁面熱流，采用當(dāng)?shù)鼐W(wǎng)格雷諾數(shù)Relocal=ρUΔh/μ[25]來確定物面法向第1層網(wǎng)格的合適高度，其中ρ、U、μ為流體密度、速度和動(dòng)力黏性系數(shù)，Δh為網(wǎng)格高度。文獻(xiàn)[25]指出：Relocal應(yīng)不大于3才能保證壁面熱流計(jì)算結(jié)果的準(zhǔn)確性。因此，流場(chǎng)網(wǎng)格在壁面法向第1層高度取1×10-6m，當(dāng)?shù)鼐W(wǎng)格雷諾數(shù)約為1.21，滿足網(wǎng)格計(jì)算要求。同時(shí)，為保證熱流計(jì)算的連續(xù)性，圓管網(wǎng)格在流固交界面處第1層高度也取1×10-6m。

2.2 計(jì)算結(jié)果與分析

為了與參考文獻(xiàn)和試驗(yàn)結(jié)果更好進(jìn)行對(duì)比，采用提出的一體化計(jì)算方法對(duì)不銹鋼材質(zhì)圓管開展高超聲速氣動(dòng)加熱與結(jié)構(gòu)傳熱非定常數(shù)值模擬，初始真實(shí)時(shí)間步長(zhǎng)取Δt0=1×10-4s。

圖5 密度等值線圖與試驗(yàn)紋影圖[26]對(duì)比
Fig.5 Comparison of predicted density contours with experimental schlieren photograph[26]

圖6 圓管表面壓力分布對(duì)比
Fig.6 Comparison of surface pressure distributions on cylinder

圖5直觀地給出了初始時(shí)刻計(jì)算得到的圓管外部流場(chǎng)的密度云圖與風(fēng)洞試驗(yàn)紋影圖[26]的對(duì)比結(jié)果，從圖中可以看出計(jì)算預(yù)測(cè)的激波位置和流場(chǎng)結(jié)構(gòu)與試驗(yàn)結(jié)果吻合較好，說明本算例暫不考慮高溫空氣化學(xué)非平衡效應(yīng)是合理的。另外，圖6為歸一化后的圓管表面壓力(p/p0，p0為駐點(diǎn)壓力)分布與試驗(yàn)值的對(duì)比，從圖中可以看出壓力沿圓管的分布和試驗(yàn)值吻合較好。

圖7給出了初始時(shí)刻(t=0 s)圓管外表面歸一化熱流(q/q0,q0為駐點(diǎn)熱流)分布與風(fēng)洞試驗(yàn)值的對(duì)比曲線，結(jié)果顯示計(jì)算預(yù)測(cè)的熱流分布與試驗(yàn)結(jié)果保持一致。另外，t=0 s時(shí)刻，圓管表面駐點(diǎn)處熱流值最大，本文方法計(jì)算得到的駐點(diǎn)處熱流值為49.21×104W/m2，與采用Fay-Riddell[27]工程經(jīng)驗(yàn)公式計(jì)算得到的駐點(diǎn)熱流值48.27×104W/m2以及采用黏性激波層[28]計(jì)算方法得到的47.02×104W/m2基本一致。但是與風(fēng)洞試驗(yàn)值67.0×104W/m2差別較大，這是由于數(shù)值計(jì)算中暫未考慮風(fēng)洞試驗(yàn)中來流湍流度對(duì)氣動(dòng)加熱熱流的影響所導(dǎo)致。

圖7t=0 s圓管表面熱流分布對(duì)比
Fig.7 Comparison of surface heat flux distributions on cylinder att=0 s

圖8展示了不同時(shí)刻下流場(chǎng)與圓管結(jié)構(gòu)沿中心線(y=0 m)溫度變化曲線。從曲線圖可以看出，高速自由來流在圓管前端產(chǎn)生脫體弓形激波，溫度從初始的241.5 K急劇躍升到2 162 K，與文獻(xiàn)[12]計(jì)算結(jié)果2 163.7 K基本吻合。另外，本文計(jì)算得到的弓形激波的位置約在-54.5 mm處，與激波經(jīng)驗(yàn)公式[29]計(jì)算得到的-54.4 mm基本吻合，說明了本文計(jì)算方法可以較準(zhǔn)確捕捉高速流場(chǎng)激波位置。初始時(shí)刻，在圓管駐點(diǎn)處，流場(chǎng)溫度急劇降低至294.4 K，在駐點(diǎn)區(qū)域會(huì)產(chǎn)生一個(gè)明顯的薄熱邊界層，其厚度約為激波層厚度的2.9%。從曲線圖可以看出，在薄熱邊界層中流場(chǎng)溫度變化非常明顯，產(chǎn)生很大溫度梯度，因此在該區(qū)域會(huì)產(chǎn)生很高的壁面熱流。另外，圖8也給出了圓管溫度沿中心線隨時(shí)間變化的放大曲線圖，可以觀察到隨著時(shí)間的推移，熱量逐漸傳入不銹鋼圓管內(nèi)部，圓管溫度逐漸升高且趨勢(shì)明顯。

圖9直觀地給出了不銹鋼圓管內(nèi)部溫度分布隨時(shí)間變化的云圖(t=0.1,0.5,1.0,2.0 s)。可以觀察到，隨著時(shí)間的變化，外部高速氣流產(chǎn)生的熱量不斷傳入圓管結(jié)構(gòu)內(nèi)部，圓管內(nèi)部溫度逐漸升高，高溫區(qū)域也不斷擴(kuò)大，駐點(diǎn)區(qū)域溫度始終最高且不斷升高。t=2.0 s時(shí)刻，本文計(jì)算得到駐點(diǎn)溫度為389.2 K，與Dechaumphai等[26]的計(jì)算結(jié)果388.8 K吻合很好，誤差在0.1%左右。

為了更加準(zhǔn)確對(duì)比驗(yàn)證本文計(jì)算結(jié)果的可靠性，圖10給出了圓管駐點(diǎn)溫度(T0)隨時(shí)間變化曲線與文獻(xiàn)[30-31]的對(duì)比結(jié)果，可以看出本文駐點(diǎn)溫度變化趨勢(shì)與參考文獻(xiàn)結(jié)果相同，t=2.0 s時(shí)，最大溫度差值約為3.9 K，計(jì)算相對(duì)誤差在0.9%以內(nèi)。另外，圖中還給出了駐點(diǎn)溫度變化率隨時(shí)間變化曲線，可以觀察到初始時(shí)刻駐點(diǎn)溫度升高速度最快，溫升速度逐漸減小趨于平穩(wěn)。

圖8 不同時(shí)刻流場(chǎng)與圓管沿中心線溫度分布
Fig.8 Temperature distribution along centerline within fluid and cylinder domains at different time

圖9 不同時(shí)刻圓管內(nèi)部溫度云圖
Fig.9 Temperature contours within cylinder domain at different time

圖10 圓管駐點(diǎn)溫度隨時(shí)間的變化
Fig.10 Temporal variation of stagnation point temperature on cylinder

圖11為圓管駐點(diǎn)熱流與熱流變化率隨時(shí)間的變化曲線。隨著時(shí)間的推進(jìn)，圓管壁面溫度逐漸升高，駐點(diǎn)區(qū)域熱邊界層內(nèi)的溫度梯度逐漸減小，導(dǎo)致駐點(diǎn)熱流逐漸下降。t=2.0 s時(shí)，熱流下降約6.5%，與Dechaumphai等[26]計(jì)算的2 s內(nèi)駐點(diǎn)熱流下降8%接近。同時(shí)從虛線可以看出駐點(diǎn)熱流下降速度隨時(shí)間變化也逐漸減小趨于平穩(wěn)。

為了更清楚地驗(yàn)證本文方法計(jì)算結(jié)果，表3給出了t=2 s時(shí)刻駐點(diǎn)溫度和初始時(shí)刻駐點(diǎn)熱流值與不同參考文獻(xiàn)的對(duì)比。從表中可以看出本文計(jì)算結(jié)果與各參考文獻(xiàn)值都吻合較好，計(jì)算誤差在允許范圍內(nèi)。綜上，本算例驗(yàn)證了本文流-熱-固一體化計(jì)算方法的可靠性和可行性，并可用于下文的計(jì)算與分析。

圖11 圓管駐點(diǎn)熱流隨時(shí)間的變化
Fig.11 Temporal variation of stagnation point heat flux on cylinder

表3 駐點(diǎn)溫度、熱流值對(duì)比

另外，為了對(duì)比本文一體化計(jì)算方法與傳統(tǒng)分區(qū)耦合計(jì)算方法的計(jì)算效率，表4列出了不同計(jì)算方法下，迭代步數(shù)與計(jì)算時(shí)間統(tǒng)計(jì)結(jié)果。本算例在圖形工作站(主頻 2.3 Hz，內(nèi)存16 GB)上進(jìn)行串行計(jì)算。

從表4中可以看出，相比傳統(tǒng)分區(qū)傳遞計(jì)算方法，一體化計(jì)算方法無數(shù)據(jù)交換計(jì)算過程，可以節(jié)省計(jì)算時(shí)間。采用自適應(yīng)時(shí)間步長(zhǎng)算法，可以進(jìn)一步大幅度提高計(jì)算效率。

表4 迭代步數(shù)與計(jì)算時(shí)間統(tǒng)計(jì)結(jié)果

3 “Ⅳ型”激波干擾算例

3.1 計(jì)算模型及網(wǎng)格

選取與2.1節(jié)中尺寸與材質(zhì)完全相同的不銹鋼圓管作為計(jì)算模型，開展高超聲速“Ⅳ型”激波干擾流-熱-固一體化計(jì)算分析研究。為了獲得高超聲速斜激波與正激波干擾產(chǎn)生的“Ⅳ型”激波干擾，本算例計(jì)算模型還包含了斜激波發(fā)生器，圖12給出了計(jì)算模型尺寸和激波發(fā)生器的相對(duì)位置。表5給出了初始流場(chǎng)來流條件，U∞為來流速度。

圖13給出了“Ⅳ型”激波干擾算例的計(jì)算網(wǎng)格和邊界條件。同樣，流場(chǎng)與結(jié)構(gòu)傳熱計(jì)算采用統(tǒng)一四邊形網(wǎng)格，圖中綠色部分為流場(chǎng)計(jì)算網(wǎng)格，紅色部分為圓管結(jié)構(gòu)傳熱計(jì)算網(wǎng)格。流場(chǎng)網(wǎng)格單元數(shù)約176 013，圓管內(nèi)部結(jié)構(gòu)網(wǎng)格單元約為49 900。流場(chǎng)與結(jié)構(gòu)耦合交界面法向第1層網(wǎng)格高度都取值為1×10-6m，當(dāng)?shù)鼐W(wǎng)格雷諾數(shù)約為1.98，圓管的熱邊界條件與2.1節(jié)中相同。

圖12 計(jì)算模型
Fig.12 Computational model

表5 “Ⅳ型”激波干擾算例的初始來流條件

圖13 “Ⅳ型”激波干擾算例的計(jì)算網(wǎng)格與邊界條件
Fig.13 Computational grids and boundary conditions of “Type Ⅳ” shock wave interference case

3.2 計(jì)算結(jié)果與分析

對(duì)不銹鋼材質(zhì)圓管分別開展高超聲速“Ⅳ型”激波干擾氣動(dòng)加熱與結(jié)構(gòu)傳熱一體化的定常與非定常數(shù)值計(jì)算，給出計(jì)算結(jié)果并進(jìn)行分析。

3.2.1 定常計(jì)算結(jié)果

在進(jìn)行氣動(dòng)加熱與結(jié)構(gòu)傳熱一體化分析之前，首先要了解高超聲速“Ⅳ型”激波干擾的流動(dòng)特征。圖14顯示了圓管前緣的“Ⅳ型”激波干擾流動(dòng)的穩(wěn)態(tài)壓力云圖與流線結(jié)構(gòu)?？梢杂^察到，流場(chǎng)具有明顯的弓形激波和透射激波的特征。計(jì)算發(fā)現(xiàn)當(dāng)入射斜激波與圓管前方的弓形激波接近正激波的位置相交時(shí)，發(fā)生“Ⅳ型”激波干擾現(xiàn)象。激波發(fā)生器產(chǎn)生的斜激波與圓管前端產(chǎn)生的弓形激波干擾產(chǎn)生一道明顯的透射激波，將弓形激波分為上下兩個(gè)部分，且激波干擾產(chǎn)生的高溫高壓使得弓形脫體激波的脫體距離明顯增大。同時(shí)，上、下兩道弓形激波波后的亞聲速區(qū)中產(chǎn)生上、下兩道剪切層，在這兩層剪切層之間產(chǎn)生一股超聲速“噴流”，超聲速“噴流”沖擊圓管壁面產(chǎn)生一道較短的弓形滯止激波，形成一個(gè)小的滯止加熱區(qū)域，引起壁面熱流、壓力急劇上升。由此看出，“Ⅳ型”激波干擾流場(chǎng)結(jié)構(gòu)是非常復(fù)雜的。

圖14 “Ⅳ型”激波干擾穩(wěn)態(tài)流場(chǎng)結(jié)構(gòu)
Fig.14 Steady flow field structure of “Type Ⅳ” shock wave interference

圖15展示了 “Ⅳ型”激波干擾流場(chǎng)與圓管結(jié)構(gòu)內(nèi)部的穩(wěn)態(tài)溫度分布。自由來流經(jīng)過弓形激波第1次加熱后，超聲速“噴流”產(chǎn)生的二次弓形滯止激波對(duì)其進(jìn)一步加熱，流場(chǎng)溫度從初始的226.5 K急劇增加，最高溫度達(dá)到3 722.1 K。不銹鋼圓管由于外部流場(chǎng)的持續(xù)氣動(dòng)加熱，內(nèi)部結(jié)構(gòu)溫度整體升高。其中，二次弓形激波在圓管壁面形成了一個(gè)局部的高溫高壓滯止加熱區(qū)域，由此導(dǎo)致圓管結(jié)構(gòu)在該區(qū)域形成一個(gè)高溫區(qū)，最高溫度可達(dá)3 406.2 K。

圖15 流場(chǎng)與圓管穩(wěn)態(tài)溫度云圖
Fig.15 Steady temperature contours within fluid and cylinder domains

圖16 圓管表面穩(wěn)態(tài)熱流與壓力系數(shù)分布曲線
Fig.16 Curves of steady surface heat flux and pressure coefficient distributions on cylinder

為了定量分析“Ⅳ型”激波干擾產(chǎn)生的高溫高壓效應(yīng)，圖16給出了圓管表面穩(wěn)態(tài)熱流與壓力系數(shù)(Cp)分布曲線。可以明顯看出，“Ⅳ型”激波干擾產(chǎn)生的超聲速“噴流”沖擊圓管壁面導(dǎo)致表面熱流與壓力系數(shù)出現(xiàn)顯著峰值。超聲速“噴流”沖擊壁面的作用位置大致位于120°的地方，在此處，表面熱流與壓力系數(shù)分別高達(dá)4 348.1 kW/m2和15.8。

值得一提的是，本定常算例也展示了本文一體化計(jì)算方法的優(yōu)勢(shì)：可以高效求解復(fù)雜高超聲速流場(chǎng)下的流-熱-固耦合定常計(jì)算問題。傳統(tǒng)分區(qū)耦合計(jì)算方法在求解定常問題時(shí)，需要在時(shí)間域上進(jìn)行大量的耦合交替迭代直至流場(chǎng)與固體結(jié)構(gòu)溫度場(chǎng)定常收斂，因此必將導(dǎo)致計(jì)算效率的大幅度降低，相比之下，本文一體化方法可以很好解決這一問題。

3.2.2 非定常計(jì)算結(jié)果

圖17給出了非定常計(jì)算中不銹鋼圓管內(nèi)部溫度分布隨時(shí)間變化的云圖(t=0.1,0.5,1.0,2.0 s)。從圖中可以直觀地觀察內(nèi)部結(jié)構(gòu)的溫度分布隨時(shí)間的演變歷程。由于“Ⅳ型”激波干擾的作用，產(chǎn)生的超聲速“噴流”不斷沖擊壁面，產(chǎn)生的高熱量首先在沖擊區(qū)域快速積累，導(dǎo)致該區(qū)域的圓管內(nèi)部結(jié)構(gòu)溫度急劇上升，形成一個(gè)顯著的高溫區(qū)。隨著時(shí)間的推移，“噴流”沖擊產(chǎn)生的高熱流持續(xù)傳入圓管結(jié)構(gòu)內(nèi)部且不斷延伸，導(dǎo)致結(jié)構(gòu)內(nèi)部整體溫度不斷升高，同時(shí)沖擊區(qū)域的高溫區(qū)域也不斷擴(kuò)大。由此可見，“Ⅳ型”激波干擾的熱沖擊作用導(dǎo)致結(jié)構(gòu)內(nèi)部產(chǎn)生集中的高溫區(qū)，在此高溫區(qū)域必然會(huì)出現(xiàn)較大的熱應(yīng)力集中現(xiàn)象，這也將是后續(xù)研究的重要內(nèi)容。

圖18和圖19分別給出了圓管表面的壓力系數(shù)分布和熱流分布隨時(shí)間的變化歷程。同時(shí)，為了更好地對(duì)比分析“Ⅳ型”激波干擾產(chǎn)生的高溫高壓效應(yīng)，圖中也給出了無斜激波干擾時(shí)，初始時(shí)刻圓管表面的壓力與熱流分布，并用此時(shí)的駐點(diǎn)壓力系數(shù)(Cp0)和熱流值對(duì)預(yù)測(cè)的“Ⅳ型”激波干擾壁面壓力系數(shù)和熱通量進(jìn)行統(tǒng)一歸一化處理。

圖17 圓管內(nèi)部溫度云圖隨時(shí)間變化歷程
Fig.17 Temporal variation of temperature contours within cylinder

圖18 圓管表面壓力系數(shù)隨時(shí)間變化
Fig.18 Pressure coefficient variation with time along cylinder surface

圖19 圓管表面熱流隨時(shí)間變化
Fig.19 Heat flux variation with time along cylinder surface

從圖18和圖19可明顯觀察到，“Ⅳ型”激波干擾的確導(dǎo)致圓管表面壓力和熱流急劇增大，且峰值都出現(xiàn)在超聲速“噴流”沖擊壁面區(qū)域。相比于無斜激波干擾的流動(dòng)，“Ⅳ型”激波干擾的射流沖擊作用使得壁面最大壓力系數(shù)增大約9倍，壁面最大熱流增大約4.7倍，計(jì)算結(jié)果與文獻(xiàn)[32-33]吻合。由此可見，此類“Ⅳ型”激波干擾將極大放大外部高速流場(chǎng)對(duì)飛行器結(jié)構(gòu)的氣動(dòng)熱與氣動(dòng)力載荷作用，給高速飛行器飛行安全造成極大安全隱患。另外，從圖中可以看出，隨著時(shí)間的推移，圓管表面壓力分布幾乎不隨時(shí)間變化，這說明圓管壁面壓力受圓管整體結(jié)構(gòu)溫度變化的影響很小。然而，壁面熱通量的峰值隨著計(jì)算時(shí)間的推進(jìn)逐漸減小。這是由于初始時(shí)刻壁面熱流峰值最大，圓管結(jié)構(gòu)吸收氣動(dòng)加熱熱量后，壁面溫度快速上升。壁溫的升高導(dǎo)致溫度邊界層內(nèi)的溫度梯度減小，壁面熱流峰值減小，削弱了外部氣動(dòng)加熱效應(yīng)從而導(dǎo)致圓管結(jié)構(gòu)溫升速率的降低。同時(shí)，結(jié)構(gòu)溫升速率的降低反過來也會(huì)導(dǎo)致壁面熱通量的變化率逐漸減小。這些現(xiàn)象正是反映了流-熱-固相互作用的耦合特性。

綜上分析，高速飛行器在進(jìn)行長(zhǎng)航時(shí)飛行時(shí)，“Ⅳ型”激波干擾作用產(chǎn)生的高溫高壓沖擊載荷，極有可能造成飛行器熱防護(hù)結(jié)構(gòu)的熱力學(xué)破壞從而形成嚴(yán)重的飛行安全隱患。因此，“Ⅳ型”激波干擾作用將給高速飛行器的熱防護(hù)設(shè)計(jì)與選材帶來嚴(yán)峻挑戰(zhàn)，需引起重要關(guān)注。

4 結(jié) 論

本文針對(duì)高超聲速進(jìn)氣道前緣“Ⅳ型”激波干擾產(chǎn)生的氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合計(jì)算問題，提出一種基于有限體積法的流-熱-固一體化計(jì)算方法，并采用經(jīng)典高超聲速二維圓管繞流算例驗(yàn)證了一體化方法的正確性與穩(wěn)定性。采用所提出的一體化計(jì)算方法對(duì)高超聲速前緣“Ⅳ型”激波干擾中的流-熱-固耦合問題進(jìn)行一體化定常/非定常計(jì)算分析，得到以下結(jié)論：

1) 采用經(jīng)典高超聲速二維圓管繞流非定常算例對(duì)本文發(fā)展的一體化計(jì)算方法進(jìn)行驗(yàn)證分析，計(jì)算結(jié)果與參考文獻(xiàn)和風(fēng)洞試驗(yàn)結(jié)果吻合較好，證明了本文高超聲速流-熱-固一體化計(jì)算方法的可行性與可靠性，同時(shí)也證明了自適應(yīng)時(shí)間步長(zhǎng)控制方法的正確性。

2) 本文所提出的高超聲速流-熱-固一體化求解方法可以高效解決流場(chǎng)與結(jié)構(gòu)傳熱穩(wěn)態(tài)求解問題，較快計(jì)算出穩(wěn)態(tài)結(jié)構(gòu)與流場(chǎng)的溫度分布。計(jì)算方法進(jìn)行全物理場(chǎng)統(tǒng)一迭代，可以很好地改善傳統(tǒng)分區(qū)耦合算法多次迭代計(jì)算效率低的不足。

3) 對(duì)高超聲速前緣“Ⅳ型”激波干擾流-熱-固耦合進(jìn)行定常/非定常一體化計(jì)算分析，計(jì)算發(fā)現(xiàn)，“Ⅳ型”激波干擾作用產(chǎn)生的超聲速“噴流”不斷沖擊壁面，使得壁面最大壓力系數(shù)增大約9倍，壁面最大熱流增大約4.7倍，給高速飛行器的熱防護(hù)設(shè)計(jì)與選材帶來嚴(yán)峻挑戰(zhàn)。

為進(jìn)一步探索高超聲速流-熱-固一體化求解方法的應(yīng)用范圍，后續(xù)將研究“Ⅳ型”激波干擾作用產(chǎn)生的高溫高壓對(duì)結(jié)構(gòu)熱力學(xué)變形產(chǎn)生的影響。