溫浩，史愛明，鄢榮

西北工業(yè)大學(xué) 航空學(xué)院 NPU-Duke空氣動力與氣動彈性聯(lián)合實驗室，西安 710072

超聲速飛行器氣動設(shè)計需要提高升阻比和降低音爆分貝值[1]。這兩個問題的核心研究對象則是超聲速飛行時必然會產(chǎn)生的激波。減弱甚至消除激波不僅可使超聲速飛行的阻力減小[2]，還可通過降低激波強度(本文以激波后壓強比激波前壓強表征)直接減小近場音爆分貝值[3-4]。

激波強度存在隨飛行馬赫數(shù)變化的規(guī)律。正激波的強度由波前來流馬赫數(shù)唯一確定[5]。對飛機設(shè)計來說，一維正激波強度規(guī)律是不夠的。史愛明和Dowell[6]研究了二維斜激波中的斜激波極值規(guī)律，表明相同流動偏角下，必然存在流線穿過斜激波時激波強度最小的來流馬赫數(shù)，且激波角與流動偏角在最優(yōu)馬赫數(shù)下呈現(xiàn)為線性關(guān)系。Emanuel[7]研究了斜激波后掠的情況，當(dāng)氣流偏角固定時，激波強度參數(shù)會隨后掠角變化出現(xiàn)極值點，并得出極值點位置馬赫數(shù)和后掠與否無關(guān)。Elaichi和Zebbiche[8]研究了總溫對圓錐流動的影響，從其研究結(jié)果的分析中可以得出：當(dāng)圓錐錐角固定時，激波強度參數(shù)隨馬赫數(shù)變化也存在極值現(xiàn)象。因研究對象的不同，Emanuel和Elaichi都未能指出斜激波極值規(guī)律的存在。

為對斜激波極值規(guī)律在黏性流動中的適用性進行研究，嘗試將無黏流動中的斜激波規(guī)律拓展到黏性流動中。在黏性流動中，激波的最終狀態(tài)不僅取決于幾何外形，還會受到邊界層的影響。對于可壓縮流動中邊界層流動參數(shù)，一般會結(jié)合Blasius不可壓層流平板邊界層解[9]進行理論計算。可壓縮層流邊界層的理論分析發(fā)展產(chǎn)生了變換方法[10-12]和參考溫度法[13-15]，建立了可壓縮邊界層與不可壓縮邊界層參數(shù)間的關(guān)系。變換法對可壓縮邊界層方程進行簡化，往往應(yīng)用于邊界層外流動參數(shù)與來流參數(shù)相差不大的情況[16]，對于復(fù)雜的工程應(yīng)用來說精度較差。工程上則常用參考溫度法進行計算，高超聲速快速計算方法[17]、乘波體設(shè)計[18]均據(jù)此對邊界層影響進行分析計算。湍流流動較為復(fù)雜，往往通過近似的邊界層速度型分布得到邊界層參數(shù)[19]；可壓縮湍流邊界層的理論計算主要是參考溫度法[15]，變換方法則沒有邊界層厚度參數(shù)與Blasius解的關(guān)系研究。

本文建立了超聲速流動中的黏性楔面激波模型，分別使用Eckert參考溫度法[18](ERT)和Illingworth-Stewartson變換法[16](IST)進行了邊界層影響的計算，并與計算流體力學(xué)(CFD)[20]方法得到的結(jié)果進行比較。結(jié)果表明，單一的ERT方法結(jié)果較CFD偏低，而IST方法則恰好較CFD偏高，且二者的這種表現(xiàn)規(guī)律在湍流情況下更加明顯。因此建立基于ERT和IST的加權(quán)均值方法(WAM)。ERT的層流、湍流模型區(qū)別在于普朗特數(shù)導(dǎo)致的壁面溫度恢復(fù)系數(shù)的差別；因湍流模型較為復(fù)雜，IST的湍流模型設(shè)置為與層流一致，忽略湍流項導(dǎo)致的差異。邊界層理論模型的結(jié)果表明，黏性楔面激波同樣存在著極值規(guī)律。

1 楔面激波理論模型

1.1 可壓縮邊界層及楔面激波

圖1 超聲速楔面黏性流動的激波結(jié)構(gòu)Fig.1 Shock wave structure of supersonic wedge in viscous flow

楔面流動模型如圖1所示，在超聲速流動中存在一個二維楔面，楔面前不存在擾動。在特定的條件下，流體經(jīng)過楔面會產(chǎn)生附體的斜激波；同時楔面附近因無滑移條件而產(chǎn)生邊界層，使得邊界層外的流線偏離楔面，偏離量可以用位移厚度表示。楔面模型邊界層的發(fā)展與平板類似，在層流邊界層和湍流邊界層處，邊界層位移厚度隨楔面長度增加，其增長率降低，流線偏離自身，產(chǎn)生膨脹波；在邊界層轉(zhuǎn)捩處，位移厚度先減小后增加，根據(jù)位移厚度增長率的變化產(chǎn)生壓縮波、膨脹波。因位移厚度變化率而產(chǎn)生的波與上游的激波相互作用，使激波略微彎曲：壓縮波使得激波強度增加，激波角增大；膨脹波則反之。

斜激波理論中，氣流方向的變化是激波產(chǎn)生的原因，經(jīng)過激波的氣流偏角θ是斜激波結(jié)構(gòu)的一個重要參數(shù)。無黏楔面流動中，經(jīng)過激波的氣流偏角等于楔面與來流的夾角——楔面角θs，即無黏楔面激波結(jié)構(gòu)具有θ=θs。黏性楔面流動中，邊界層會導(dǎo)致氣流偏角發(fā)生改變，全層流和全湍流時均具有θ>θs。1.2、1.3節(jié)將構(gòu)建受邊界層影響的氣流偏角的計算方法，從而在1.4節(jié)中應(yīng)用斜激波、膨脹波理論計算激波強度。

1.2 邊界層外流動參數(shù)

平板流動中，使用邊界層理論計算邊界層位移厚度等參數(shù)時，主要基于自由來流參數(shù)。而在楔面激波模型中，為使得計算更加準確，考慮以無黏理論確定的斜激波后參數(shù)進行計算。主要計算激波下游的單位雷諾數(shù)Reds和激波下游溫度T2。根據(jù)Sutherland公式計算黏性系數(shù)

(1)

式中:T為溫度；參考溫度Tv=273.2 K；常數(shù)c=110.4 K；黏性系數(shù)μ(Tv) = 1.72×10-5kg/(m·s)。

根據(jù)無黏理論，在已知來流馬赫數(shù)Ma1和無黏氣流偏角θs的條件下，可得激波后溫度與激波前溫度的比值T2/T1，從而得到激波后溫度T2。

(2)

式中：T1=300 K為設(shè)置的激波上游溫度；Man1為激波上游法向馬赫數(shù)；γ=1.4為比熱比；下標1表示激波前參數(shù)，下標2表示無黏理論確定的激波后參數(shù)，后文與此相同。根據(jù)式 (1)和式(2)，可以得到激波后黏性系數(shù)的比值μ2/μ1，進而由ρ1u1=ρ2u2(ρ、u為密度和速度)得到單位雷諾數(shù)之比及激波下游邊界層外的單位雷諾數(shù)Reds，表達式為

(3)

其中:Reus為斜激波上游的單位雷諾數(shù)。

1.3 可壓縮邊界層計算方法

1.3.1 層流邊界層計算方法

參考溫度法模型以參考雷諾數(shù)來代替不可壓縮Blasius解中的雷諾數(shù)，基于此進行邊界層的計算，其模型的發(fā)展更多關(guān)注于參考溫度的算法研究。ERT方法確定了參考溫度T*，其表達式為

T*=0.5(Tw+T2)+0.22(Tr-T2)

(4)

根據(jù)參考溫度法，得到參考雷諾數(shù)Re*=ρ*ued/μ*，其中：ρ*為參考密度；ue為邊界層外速度；d=1 m為特征長度；μ*為參考黏性系數(shù)。根據(jù)理想氣體狀態(tài)方程ρ*=p*/RT*(p*=p2表示激波下游壓強，R=287 J/(kg·K)表示氣體常數(shù))，可以得到

(5)

式中：μ*由T*確定。取距離楔面頂點水平距離x位置為層流邊界層位移厚度的計算位置，則ERT方法確定的位移厚度為

(6)

求導(dǎo)可得x位置處流動方向相對楔面方向的偏轉(zhuǎn)角θΔ，其表達為

(7)

那么最終x位置處流線方向相對來流方向的偏轉(zhuǎn)角θef=θs+θΔ，如圖2(a)所示。

變換法也可建立與不可壓縮流動相似的可壓縮層流邊界層參數(shù)表示方法，該方法中Pr=1。對于ERT方法，Pr對最終理論模型結(jié)果的影響遠小于模型與CFD結(jié)果間的差距，可以忽略不計；這表明ERT方法和IST方法間的差距不在于Pr的不同，而在于自身模型的構(gòu)建過程。IST方法確定的位移厚度為

(8)

圖2 考慮邊界層影響的激波參數(shù)計算模型
Fig.2 Computation model of shock parameters with boundary layer effect considered

1.3.2 湍流邊界層計算方法

IST法中，式(8)相關(guān)的不可壓縮湍流邊界層厚度為

其余參數(shù)不變。之后計算可壓縮湍流邊界層造成的流線方向改變量為

1.4 可壓縮邊界層流場影響

對于全層流或全湍流流動，在來流馬赫數(shù)Ma1、楔面角θs、來流單位雷諾數(shù)Reus確定的情況下，隨著距楔面頂點水平位置x增大，δx逐漸增大，但θΔ卻逐漸減小。楔面流動中，楔面頂點附近壁面上會存在流動駐點，該位置處流動壓強等于總壓，之后隨流動發(fā)展，壓強會減小。在距離駐點足夠遠處，駐點總壓的影響可以忽略不計。邊界層理論可以確定不同x處的流線方向：x較小時，流動偏角增量較大；隨x增大，流動偏角增量降低。

對于邊界層影響量的計算，本文建立了斜激波和膨脹波的組合理論模型，用以計算x位置的無量綱壓強px/p1，也稱為激波強度。如圖2(b)所示，在位置b處，使用斜激波理論確定位置b的流動參數(shù)，以之模擬流動駐點的影響。依據(jù)可壓縮邊界層理論計算初始位置的氣流偏角θef,b；根據(jù)無黏斜激波理論，使用來流馬赫數(shù)Ma1和氣流偏角θef,b，確定初始斜激波強度px,b/p1等流場參數(shù)。假設(shè)b到e為等熵膨脹過程，根據(jù)b處流場信息，通過膨脹波理論和Prandtl-Meyer函數(shù)，計算e處對應(yīng)的激波強度px,e/p1等參數(shù)。

膨脹波關(guān)系表示為

θef,b-θef,e=ν(Mae)-ν(Mab)

(9)

式中：Mab和Mae分別為b處和e處的馬赫數(shù)；ν分Prandtl-Meyer函數(shù)，有

(10)

其中：α=Ma2-1。根據(jù)等熵過程總參數(shù)不變即可求得

(11)

進而得到e處的激波強度px,e/p1。

至此，在楔面激波模型中，將斜激波和邊界層理論發(fā)展到了邊界層外流場參數(shù)的計算中。邊界層理論模型——ERT和IST方法會得出不同的激波強度計算結(jié)果。本節(jié)模型適用于邊界層為全層流或全湍流的情況。對于流動轉(zhuǎn)捩，需考慮轉(zhuǎn)捩過程中邊界層位移厚度的變化規(guī)律，同時在轉(zhuǎn)捩前后分別應(yīng)用層流和湍流理論模型。

2 理論模型預(yù)測精度分析

2.1 數(shù)值方法與網(wǎng)格

使用CFD的方法進行評估時，分別進行無黏Euler方程、層流Navier-Stokes方程、湍流基于雷諾平均Navier-Stokes(RANS)方程的計算。本文Euler方程空間離散為二階守恒的單調(diào)迎風(fēng)格式MUSCL[21]，與國內(nèi)提出的無波動、無自由參數(shù)的耗散差分格式NND[22]相比，二者在數(shù)值計算的應(yīng)用層面完全一致；限制器類型為Venkatakrishnan[23]，對流項處理的數(shù)值格式為Harten-Lax-van Leer-Contact (HLLC)[24]；時間離散為歐拉隱式(Euler-Implicit)；采用當(dāng)?shù)貢r間步、CFL數(shù)自適應(yīng)、多重網(wǎng)格等方法加速收斂。層流Navier-Stokes方程無黏部分的離散格式與Euler方程一致，黏性項離散采用中心差分。湍流RANS方程湍流模型為Spalart-Allmaras(S-A)，湍流模型方程空間離散為一階標量迎風(fēng)(Scalar Upwind)格式，時間離散為歐拉隱式[20]。

二維超聲速楔面黏性流動的計算網(wǎng)格如圖3所示，采用非等距四邊形網(wǎng)格。網(wǎng)格拓撲包括：與壁面平行的邊界層網(wǎng)格，與激波平行的邊界層外區(qū)域網(wǎng)格。沿壁面方向網(wǎng)格尺度最小為2×10-5m，線性增長，增長率為1.2，數(shù)據(jù)采集點處約為2×10-3m；垂直壁面方向網(wǎng)格尺度最小為4×10-7m，線性增長，增長率為1.2，數(shù)據(jù)采集點處約為7×10-4m。網(wǎng)格區(qū)域大小為0.4 m×1.2 m，總數(shù)為177×170；加密區(qū)域為激波和壁面附近，大小為0.03 m×0.02 m，數(shù)目為51×49。所有黏性CFD計算中，單位雷諾數(shù)Reus=3.5×107/m；所有CFD計算中，激波強度采集位置xe≈0.085 m。

圖3 楔面激波數(shù)值模擬網(wǎng)格
Fig.3 Numerical mesh of wedge shock

2.2 理論模型中膨脹波部分的優(yōu)化

本節(jié)及2.3節(jié)的CFD計算中，楔面角為3°，馬赫數(shù)計算點取Ma1=1.195 7、1.234 9、1.319 0、1.368 6、1.469 1、1.600 0、1.700 0、2.000 0、2.300 0，其中Ma1=1.469 1產(chǎn)生的無黏激波強度最小。

理論模型中，對流動駐點影響的考慮在于xb處定義的激波強度計算，模型中膨脹波理論描述了邊界層變化對激波強度的影響。模型中尚有xb未能確定，為分析其對xe處激波強度預(yù)測結(jié)果的影響，取Ma1=1.195 7、1.469 1、2.300 0結(jié)果進行研究。圖4為層流ERT結(jié)果的理論模型計算值分布，可以看到在xb>0.01 m(約xe/8位置)后，預(yù)測結(jié)果隨xb變化不大。在xb<0.01 m時，xb的影響較為明顯。特別是當(dāng)xb接近0時，與其他xb相比，理論模型預(yù)測結(jié)果會產(chǎn)生較大變化。這是因為邊界層理論模型在xb為0時出現(xiàn)奇異值，使得預(yù)測的θΔ趨于90°。此時起始位置激波造成的壓強增大，之后膨脹波造成的壓強降低，二者在量值上都很大，造成最終計算精度的降低。xb>0.01 m后，xb影響的標準差為1.471 3×10-5，遠小于理論模型與CFD結(jié)果間差值的標準差1.581 8×10-4；IST結(jié)果、湍流結(jié)果也是如此。這表明在處理簡單的超聲速楔面流動時，單一的斜激波理論可以達到較好的精度。因此，后續(xù)理論模型對于激波強度的預(yù)測，可以直接在xe處應(yīng)用斜激波理論，采集點處壓強寫為px。

圖4 初始激波位置xb對層流ERT預(yù)測結(jié)果的影響
Fig.4 Influence of initial shock locationxbon prediction results of laminar ERT model

2.3 單一ERT和IST理論方法結(jié)果分析

下面根據(jù)層流和湍流的數(shù)值結(jié)果對理論模型的預(yù)測結(jié)果進行分析。圖5為理論模型預(yù)測結(jié)果與CFD間差距分析，圖中CFD表示層流或湍流計算結(jié)果，WAM方法將在2.4節(jié)具體敘述?？梢钥吹剑瑹o論是層流還是湍流，在激波強度的預(yù)測上，ERT、IST兩種理論模型得到了隨馬赫數(shù)變化的極值現(xiàn)象，均與無黏理論表現(xiàn)一致?？梢钥吹紼RT的預(yù)測結(jié)果小于CFD結(jié)果，而IST的預(yù)測結(jié)果大于CFD結(jié)果，這在湍流結(jié)果中表現(xiàn)得更加明顯。這說明兩種邊界層方法具有各自的內(nèi)在特點，且分別表現(xiàn)了對結(jié)果的低估和高估。因此考慮兩種模型的加權(quán)平均結(jié)果，以取得較好的黏性結(jié)果預(yù)測精度。

圖5θs=3°理論模型和CFD計算結(jié)果的對比
Fig.5 Comparison between theoretical model and CFD results atθs=3°

2.4 WAM方法及其結(jié)果分析

WAM方法的表達式為

(12)

式中:λ為加權(quán)因子。ERT和IST方法均采用無黏激波后參數(shù)作為邊界層參數(shù)來源。根據(jù)二者的理論方法構(gòu)建過程，自變量參數(shù)可單一取為Ma2。

在λ的擬合上，考慮如表1所示的楔面角和來流馬赫數(shù)范圍來產(chǎn)生訓(xùn)練集和測試集。θs=3°、5°、10°在相應(yīng)的馬赫數(shù)范圍內(nèi)每組取9個馬赫數(shù)點；而θs=20°、30°、35°時，每組取5個馬赫數(shù)點。各楔面角下馬赫數(shù)的設(shè)置基于激波強度極值規(guī)律，保證激波強度最小值點的存在。

Euler方程的計算結(jié)果表明，隨楔面角的增大，數(shù)值結(jié)果與無黏理論值間的差距增大，即大楔面角下，CFD結(jié)果不夠準確，如表2所示。層流和湍流的數(shù)值計算也是如此，故λ的擬合考慮角度θs=3°、5°、10°。因CFD結(jié)果與方程計算的精確值相比有一定的差距，但相較于理論與CFD間的誤差，這一差距并不大。因此，采用層流結(jié)果難以訓(xùn)練出較為準確的λ值。而湍流結(jié)果中，這一誤差相比之下卻足夠小，得到的擬合前的λ分布也較為穩(wěn)定。且湍流與層流使用的是相同的可壓縮處理方法，所以本文使用湍流計算結(jié)果來擬合λ，減少CFD本身誤差的影響。

表1 CFD狀態(tài)點范圍Table 1 Condition range of CFD

λ經(jīng)三次多項式擬合得到

(13)

下面對WAM模型的預(yù)測結(jié)果使用測試集進行分析，測試集選為全部計算黏性結(jié)果。WAM與CFD間的相對誤差為

(14)

式中：i=無黏理論，層流WAM，湍流WAM；j=無黏CFD，層流CFD，湍流CFD。對同一楔面角下的所有馬赫數(shù)結(jié)果相對誤差求標準差，得到相對誤差標準差分布如表2所示。隨楔面角增大，WAM與CFD間的相對誤差增大，在楔面角小于20°時，此時馬赫數(shù)范圍約為1.23.1，湍流WAM預(yù)測結(jié)果誤差與其他楔面角和層流結(jié)果匹配較差，這表明WAM模型在楔面角約大于30°的湍流結(jié)果的可信度開始迅速下降，預(yù)測相對誤差大于0.01；而層流結(jié)果則仍能保持一定精度。

表2 WAM與CFD間的相對誤差標準差

在無黏斜激波規(guī)律中，隨著流動偏角θ的增大，最優(yōu)馬赫數(shù)與激波恰好脫體時的馬赫數(shù)間的差距逐漸減小，這意味著斜激波極值規(guī)律在大流動偏角，或者說高馬赫數(shù)時，其適用性較差。因此對于所要研究的斜激波極值規(guī)律，WAM方法可以在合適的范圍內(nèi)取得較好的預(yù)測精度。

3 邊界層對斜激波極值規(guī)律的影響

3.1 無黏斜激波極值規(guī)律討論

無黏流動中的斜激波極值規(guī)律表明，在相同的氣流偏轉(zhuǎn)角θ下，具有最小激波強度的波前馬赫數(shù)并不是激波恰好脫體的馬赫數(shù)值，反而大于該值；稱該最小激波強度下的來流馬赫數(shù)為最優(yōu)馬赫數(shù)。該激波強度可以用壓強比、總壓損失率、溫度比衡量，在激波強度于最優(yōu)馬赫數(shù)下取得極值時，壓強比等參數(shù)均會取得極值。該規(guī)律基于斜激波理論得到的法向馬赫數(shù)關(guān)系式為[10]

(15)

式中：θ為經(jīng)過激波的氣流偏角，在無黏的超聲速楔面模型中，其等于楔面角θs；β為激波角，即激波面與來流方向的夾角。此關(guān)系式建立了斜激波空間結(jié)構(gòu)參數(shù)θ、β與激波強度唯一變量Man1的關(guān)系?；诖?，得到最優(yōu)馬赫數(shù)為

(16)

同樣地，在確定馬赫數(shù)作為飛機設(shè)計約束、要求更小的飛行能量損失的情況下，穿過激波的流動偏角會趨向于0。因此在沒有其他約束作用的情形下，無法合理地設(shè)計出具有一定體積的超聲速外形，亞聲速飛行器設(shè)計也是如此，空氣動力學(xué)特性決定了更小飛行能量損失的飛行器永遠是厚度更小的類似設(shè)計。上述斜激波極值規(guī)律表明對于楔面激波，盡管在給定馬赫數(shù)下不存在最優(yōu)的飛行器設(shè)計，然而一旦根據(jù)其他約束條件確定飛行器外形，那么便可以使用極值規(guī)律衡量飛行器外形和馬赫數(shù)需求的匹配度。

關(guān)于直接由飛行器外形變化產(chǎn)生的附體斜激波，主要有兩種形式。其一是本文研究的楔面激波類型，特點是斜激波的上游不存在邊界層[25]；其二為楔角流動類型，斜激波前存在邊界層，因此產(chǎn)生的斜激波與邊界層相互作用，使得激波振蕩，流場與無黏理論解間差距較大[26-28]。這為斜激波極值規(guī)律的應(yīng)用提供了新的思路，比如對于楔角流動所產(chǎn)生的斜激波，其最優(yōu)馬赫數(shù)如何變化、其振蕩幅值和頻率是否表現(xiàn)對稱性，這都是可以深入研究的方向。此外，經(jīng)典的入/反射激波與平板邊界層相互作用，分離泡、流動轉(zhuǎn)捩的相關(guān)特征會如何變化，這是斜激波極值規(guī)律應(yīng)用的另一個探索[29-30]。

3.2 最優(yōu)馬赫數(shù)及激波強度變化

數(shù)值驗證的結(jié)果表明，在最優(yōu)馬赫數(shù)附近，楔面角不大時，WAM方法取了得良好的精度。

首先考慮最優(yōu)馬赫數(shù)的變化。在加入了黏性邊界層的影響后，楔面角θs≈3°～20°范圍內(nèi)，關(guān)于最優(yōu)馬赫數(shù)的變化，層流增量為0.001 5～0.003 3，湍流增量比層流稍大，為0.002 8～0.006 1，且二者均隨楔面角增大而增大。圖6為加入邊界層影響后，最優(yōu)馬赫數(shù)Ma1隨θs變化的情況。這表明邊界層不會使得該Ma1產(chǎn)生較大變化；由局部放大圖可知湍流Ma1增量約為層流2倍，20°時的增量約為3°時的2倍。

根據(jù)楔面邊界層變化規(guī)律，在距楔面頂點水平位置大于參考位置0.085 m處，最優(yōu)馬赫數(shù)的增量會更小，反之則更大，這里取θs=10°進行討論。在百倍于參考位置，即約10 m處，黏性最優(yōu)馬赫數(shù)的增量相比參考位置減小，層流不大于0.000 2，湍流不大于0.000 6。在距楔面頂點約0.01 m，接近理論模型奇異點處，最優(yōu)馬赫數(shù)增量比參考位置增加，層流和湍流均不大于0.000 6。這表明在約0.01～10 m的較大范圍內(nèi)，于參考位置0.085 m處得出的最優(yōu)馬赫數(shù)的變化規(guī)律可對整個區(qū)域的結(jié)果進行表征。所以之后的討論仍是基于參考位置的數(shù)據(jù)。

圖6 邊界層對最小激波強度馬赫數(shù)的影響
Fig.6 Changes of Mach number of minimum shock strength influenced by boundary layer

其次考慮激波強度受到的影響。激波強度增量Δ(px/p1)與最優(yōu)馬赫數(shù)增量類似，其相比無黏激波強度值變化很小。圖7為層流狀態(tài)下，激波強度增量Δ(px/p1)的分布。圖中：LMS表示無黏最小激波強度線，有黏下的最小激波強度線基本與無黏結(jié)果一致；βs指無黏超聲速楔面模型中，由Ma1和θs經(jīng)無黏斜激波理論確定的激波角。圖中激波強度增量分布存在極小值點，極值點馬赫數(shù)大于最優(yōu)馬赫數(shù)。

圖7 激波強度增量分布
Fig.7 Distribution of increments of shock strength

3.3 激波強度變化的理論解釋

為理解從無黏斜激波關(guān)系到加入邊界層影響后的激波強度和最優(yōu)馬赫數(shù)的變化，將受邊界影響的激波強度px/p1表示為無黏激波強度項與激波強度增量項之和：

(17)

式中：p2/p1為無黏激波強度。根據(jù)斜激波理論，激波強度自變量取為來流馬赫數(shù)和楔面角，即p2/p1=f(Ma1,θs)。并且最終，邊界層導(dǎo)致的流線偏角變化仍是代入無黏斜激波理論，從而計算激波強度變化。根據(jù)無黏斜激波理論，Δ(px/p1)的一階展開為

(18)

其中：

(19)

圖8為無黏激波強度梯度項?(p2/p1)/?θs和層流流動偏角增量θΔ的空間分布。無黏激波強度梯度項的極值點馬赫數(shù)異于最優(yōu)馬赫數(shù)。層流流動偏角增量項不存在極值點，其隨馬赫數(shù)單調(diào)增加。

圖8 激波強度增量項的分解結(jié)果
Fig.8 Terms decomposed from increments of shock strength

圖9為LMS線上湍流和層流情況下的激波強度增量比值及流動偏角增量比值。LMS線上，湍流自約θs=37.9°開始激波全部脫體，而層流自約θs=40.4°開始脫體。圖9中，在較大的楔面角范圍內(nèi)，激波強度增量比值和流動偏角增量比值保持了較好的一致性：楔面角小于20°時，二者相對差距不超過0.5%；在楔面角小于25°時，二者相對差距不超過1%，在楔面角小于30°時，二者相對差距不超過2.5%。因此，在WAM-CFD誤差較小的楔面角范圍內(nèi)(3°～20°)，Δ(px/p1)使用式(18)預(yù)測，其誤差相比精確的斜激波理論不大。在楔面角大于35°后，二者差距已超過7%，此時不宜使用式(18)進行Δ(px/p1)的預(yù)測。

圖9 LMS線上湍流和層流偏角增量、激波強度 增量比值
Fig.9 Ratio of increments at turbulent and laminar case for deflection angle and shock strength at line LMS

至此確定了在LMS線附近的區(qū)域內(nèi)，邊界層導(dǎo)致的激波強度增量和其導(dǎo)致的流動偏角增量呈較為準確的線性關(guān)系，系數(shù)為無黏激波強度梯度項。相對于激波強度的極值點馬赫數(shù)，激波強度增量與激波強度梯度的極值點馬赫數(shù)均較大，這一分布關(guān)系使得最優(yōu)馬赫數(shù)增大。

4 結(jié) 論

1) 建立了黏性楔面激波強度計算的理論模型，初步定量地確定了邊界層厚度對楔面激波強度影響規(guī)律，為超聲速黏性楔面流動問題研究提供了精度尚可的斜激波參數(shù)理論計算方法。通過研究計入邊界層等效厚度影響的斜激波強度極值規(guī)律隨馬赫數(shù)、楔面角的變參規(guī)律，將斜激波強度極值規(guī)律拓展到了黏性楔面流動范疇。

2) 理論模型顯示：① 邊界層位移厚度使斜激波極值規(guī)律最優(yōu)馬赫數(shù)略微增大，增加量在千分位；② 當(dāng)楔面角小于20°、馬赫數(shù)小于2.3時，理論模型結(jié)果與CFD結(jié)果的相對誤差不大于0.1%；③ 層流與湍流邊界層對斜激波極值規(guī)律的影響程度不同；湍流邊界層導(dǎo)致的氣流偏角增量約是層流邊界層的2倍，最優(yōu)馬赫數(shù)增量也約是層流邊界層的2倍。

3) ERT和IST的加權(quán)模型僅適用于全層流或全湍流流動，且對于高馬赫數(shù)(Ma>2.5)和大楔面角(θ>30°)情況，預(yù)測結(jié)果開始變得不理想(相對誤差大于0.2%)。計入轉(zhuǎn)捩因素的楔面斜激波理論模型是進一步的研究方向。文章討論了楔面上游為自由來流的楔面斜激波，楔面上游為壁面的楔面斜激波則是接下來的研究內(nèi)容。這兩種激波是超聲速飛機內(nèi)、外流問題生成斜激波的主要形式。

致 謝

感謝審稿人對論文研究工作的點評和實質(zhì)性提高論文撰寫質(zhì)量的指導(dǎo)。感謝聯(lián)合實驗室合作者美國杜克大學(xué)Earl H. DOWELL院士對本文研究內(nèi)容的探討及英文摘要的修改。