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

溫浩，史愛(ài)明，鄢榮

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

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

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

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

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

1 楔面激波理論模型

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

1.3 可壓縮邊界層計(jì)算方法

1.3.1 層流邊界層計(jì)算方法

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

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

(4)

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

(5)

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

(6)

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

(7)

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

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

(8)

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

1.3.2 湍流邊界層計(jì)算方法

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

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

1.4 可壓縮邊界層流場(chǎng)影響

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

對(duì)于邊界層影響量的計(jì)算，本文建立了斜激波和膨脹波的組合理論模型，用以計(jì)算x位置的無(wú)量綱壓強(qiáng)px/p1，也稱為激波強(qiáng)度。如圖2(b)所示，在位置b處，使用斜激波理論確定位置b的流動(dòng)參數(shù)，以之模擬流動(dòng)駐點(diǎn)的影響。依據(jù)可壓縮邊界層理論計(jì)算初始位置的氣流偏角θef,b；根據(jù)無(wú)黏斜激波理論，使用來(lái)流馬赫數(shù)Ma1和氣流偏角θef,b，確定初始斜激波強(qiáng)度px,b/p1等流場(chǎng)參數(shù)。假設(shè)b到e為等熵膨脹過(guò)程，根據(jù)b處流場(chǎng)信息，通過(guò)膨脹波理論和Prandtl-Meyer函數(shù)，計(jì)算e處對(duì)應(yīng)的激波強(qiá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ù)等熵過(guò)程總參數(shù)不變即可求得

(11)

進(jìn)而得到e處的激波強(qiáng)度px,e/p1。

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

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

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

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

二維超聲速楔面黏性流動(dòng)的計(jì)算網(wǎng)格如圖3所示，采用非等距四邊形網(wǎng)格。網(wǎng)格拓?fù)浒ǎ号c壁面平行的邊界層網(wǎng)格，與激波平行的邊界層外區(qū)域網(wǎng)格。沿壁面方向網(wǎng)格尺度最小為2×10-5m，線性增長(zhǎng)，增長(zhǎng)率為1.2，數(shù)據(jù)采集點(diǎn)處約為2×10-3m；垂直壁面方向網(wǎng)格尺度最小為4×10-7m，線性增長(zhǎng)，增長(zhǎng)率為1.2，數(shù)據(jù)采集點(diǎn)處約為7×10-4m。網(wǎng)格區(qū)域大小為0.4 m×1.2 m，總數(shù)為177×170；加密區(qū)域?yàn)榧げê捅诿娓浇?，大小?.03 m×0.02 m，數(shù)目為51×49。所有黏性CFD計(jì)算中，單位雷諾數(shù)Reus=3.5×107/m；所有CFD計(jì)算中，激波強(qiáng)度采集位置xe≈0.085 m。

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

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

本節(jié)及2.3節(jié)的CFD計(jì)算中，楔面角為3°，馬赫數(shù)計(jì)算點(diǎn)取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)生的無(wú)黏激波強(qiáng)度最小。

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

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

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

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

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

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

WAM方法的表達(dá)式為

(12)

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

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

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

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

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

(13)

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

(14)

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

表2 WAM與CFD間的相對(duì)誤差標(biāo)準(zhǔn)差

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

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

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

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

(15)

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

(16)

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

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

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

數(shù)值驗(yàn)證的結(jié)果表明，在最優(yōu)馬赫數(shù)附近，楔面角不大時(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變化的情況。這表明邊界層不會(huì)使得該Ma1產(chǎn)生較大變化；由局部放大圖可知湍流Ma1增量約為層流2倍，20°時(shí)的增量約為3°時(shí)的2倍。

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

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

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

圖7 激波強(qiáng)度增量分布
Fig.7 Distribution of increments of shock strength

3.3 激波強(qiáng)度變化的理論解釋

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

(17)

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

(18)

其中：

(19)

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

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

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

圖9 LMS線上湍流和層流偏角增量、激波強(qiáng)度 增量比值
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)致的激波強(qiáng)度增量和其導(dǎo)致的流動(dòng)偏角增量呈較為準(zhǔn)確的線性關(guān)系，系數(shù)為無(wú)黏激波強(qiáng)度梯度項(xiàng)。相對(duì)于激波強(qiáng)度的極值點(diǎn)馬赫數(shù)，激波強(qiáng)度增量與激波強(qiáng)度梯度的極值點(diǎn)馬赫數(shù)均較大，這一分布關(guān)系使得最優(yōu)馬赫數(shù)增大。

4 結(jié) 論

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

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

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

致 謝

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