陳以勒,俞凱凱,徐驚雷
南京航空航天大學(xué) 能源與動力學(xué)院,南京 210016
超燃沖壓發(fā)動機(jī)是高超聲速飛行器的主要動力選擇,具有結(jié)構(gòu)簡單、無需自身攜帶氧化劑、有效載荷大等優(yōu)點(diǎn),作為高超聲速技術(shù)實(shí)現(xiàn)的關(guān)鍵,是近年來國際研究的熱點(diǎn)之一[1-4]。超燃沖壓發(fā)動機(jī)一般主要由4部分構(gòu)成,分別為進(jìn)氣道、隔離段、燃燒室和尾噴管[5-6]。其中尾噴管是提供推力的重要部件,特別是當(dāng)馬赫數(shù)超過6時(shí),噴管可以產(chǎn)生約占整個(gè)推進(jìn)系統(tǒng)70%的凈推力[7]。這說明噴管的性能將直接影響超燃沖壓發(fā)動機(jī)乃至高超聲速飛行器的性能,因此研究者們嘗試各種手段以提高噴管的氣動性能[8-11]。此外,由于超高聲速飛行器通常要求機(jī)體/推進(jìn)系統(tǒng)一體化,進(jìn)而對噴管的幾何尺寸產(chǎn)生了一定限制,因此噴管設(shè)計(jì)方法除了需要考慮其氣動性能的優(yōu)化,還需要兼顧幾何約束的要求。
為提升噴管的氣動性能,Rao[12]采用變分法推導(dǎo)了最大推力理論,在此基礎(chǔ)上提出了一種最大推力的火箭噴管設(shè)計(jì)方法。最大推力理論在后來的噴管設(shè)計(jì)中得到廣泛的應(yīng)用。Mo等[13]基于最大推力理論提出了一種非對稱噴管的設(shè)計(jì)方法,并應(yīng)用該方法考慮非均勻來流設(shè)計(jì)噴管,最終可使噴管推力系數(shù)、升力和俯仰力矩分別提高1.75 %、 6.51 %和 6.25 %。Lv等[14]提出了一種基于特征線法的噴管設(shè)計(jì)方法,討論了下壁面長度對噴管性能的影響,該方法可以使噴管軸向推力系數(shù)、升力以及俯仰力矩分別增加5.5%、1 098.2% 和20.3%。Lu等[15]將Rao理論與流線追蹤技術(shù)相結(jié)合,設(shè)計(jì)了三維超燃沖壓發(fā)動機(jī)噴管,結(jié)果表明該方法可使噴管的推力和升力分別提升2.7%和69.5%。
為滿足噴管的幾何約束,Argrow和Emanuel[16]發(fā)展了最短噴管理論并用于噴管設(shè)計(jì),但采用該方法設(shè)計(jì)得到的理想膨脹噴管仍然較長,難以滿足推力噴管的幾何約束。Shyne和Keith[17]提出可以利用截短方法使噴管滿足尺寸約束并達(dá)到減重的目的,結(jié)果表明需要在噴管的重量和性能損失上做出權(quán)衡。Hoffman[18]提出了一種壓縮型截短噴管的設(shè)計(jì)方法,應(yīng)用該設(shè)計(jì)方法可有效縮短噴管長度,且壓縮型截短噴管的性能與 Rao 噴管性能相差僅0.04%~0.34%。
在現(xiàn)有的噴管設(shè)計(jì)方法中,盡管有些方法可以滿足部分尺寸約束[16-18],但卻未能充分利用幾何空間來提升噴管氣動性能。因此,本文提出一種幾何尺寸約束的超燃沖壓發(fā)動機(jī)推力噴管設(shè)計(jì)方法,實(shí)現(xiàn)對噴管長度和高度的約束,同時(shí)優(yōu)化其氣動性能。并且通過2個(gè)設(shè)計(jì)因子實(shí)現(xiàn)對噴管設(shè)計(jì)的主動控制,以滿足實(shí)際應(yīng)用過程中針對噴管的不同性能需求。本文首先詳細(xì)介紹了提出的噴管設(shè)計(jì)原理與過程;其次校核了所使用的數(shù)值模擬方法的有效性,并確定了網(wǎng)格分辨率;然后,采用數(shù)值模擬對本文提出的設(shè)計(jì)方法進(jìn)行了有效性驗(yàn)證;在此基礎(chǔ)上,開展了2個(gè)設(shè)計(jì)因子對噴管性能影響的研究;最后,通過與典型截短設(shè)計(jì)方法進(jìn)行對比分析,驗(yàn)證了本文提出的設(shè)計(jì)方法的優(yōu)越性。
特征線法(Method of Characteristic, MOC)是一種用于求解雙曲型偏微分方程的通用數(shù)學(xué)方法,其在超聲速領(lǐng)域得到了廣泛的應(yīng)用[12-15]。定常二維、平面或軸對稱的等熵超聲速流動的控制方程如式(1)所示。通過有旋特征線法可以對目標(biāo)流場構(gòu)造特征線網(wǎng)格,如圖1所示。其中任意節(jié)點(diǎn)均會發(fā)出3條特征線,分別為左行特征線C+,右行特征線C-和流線C0,分別對應(yīng)特征方程式(2)和式(3)。
圖1 特征線網(wǎng)格
(1)
(2)
(3)
式中:θ為流動角;α為馬赫角;ρ為密度;p為靜壓;a為聲速;u為x方向的速度分量;v為y方向的速度分量;ux為u在x方向的導(dǎo)數(shù),其余類似表達(dá)含義類似;δ表示流動類型(對于平面流取0,軸對稱流取1)。特征方程對應(yīng)的相容性方程為
ρVdV+dp=0
(4)
dp-a2dρ=0
(5)
(6)
其中:Ma為馬赫數(shù);V為速度。通過數(shù)值程序求解式(2)~式(6),就可以獲得目標(biāo)流場內(nèi)各節(jié)點(diǎn)的氣動參數(shù)。此外,根據(jù)節(jié)點(diǎn)的類型可以分為幾種不同的單元過程,本文提出的設(shè)計(jì)方法主要包括3種類型的節(jié)點(diǎn):內(nèi)點(diǎn)、壁面點(diǎn)和變向區(qū)域點(diǎn),單元點(diǎn)的具體求解方法可以參考文獻(xiàn)[19]。
為提升噴管氣動性能,本文提出的設(shè)計(jì)方法采用了最大推力理論。Rao提出的最大推力理論目前已得到廣泛應(yīng)用[13-15],研究結(jié)果表明最大推力理論具有顯著的優(yōu)越性。在此基礎(chǔ)上發(fā)展出了二維最大推力理論,二維情況時(shí)最后左行特征線上的氣動參數(shù)需滿足式(7)和式(8),具體的理論推導(dǎo)過程可見文獻(xiàn)[14]。
(7)
ρV2sin2θtanα=C2
(8)
式中:C1、C2均為常數(shù)。圖2為幾何尺寸約束的推力噴管設(shè)計(jì)方法示意圖。噴管內(nèi)各區(qū)域的設(shè)計(jì)細(xì)節(jié)如圖3所示。因此,為實(shí)現(xiàn)噴管推力最大化,沿最后左行特征線KuM上的氣動參數(shù)均需滿足式(7)和式(8),其中C1、C2的值可由上位關(guān)鍵點(diǎn)Ku的氣動參數(shù)計(jì)算得到。同理,沿最后右行特征線KdN上的氣動參數(shù)可采用相似的方法確定[20]。值得注意的是,本文提出的設(shè)計(jì)方法根據(jù)幾何約束確定最大推力理論中的壓力條件。
圖2 噴管設(shè)計(jì)過程示意圖
圖3 各區(qū)域設(shè)計(jì)過程
為克服以往方法未能充分利用幾何空間提升噴管氣動性能的問題,本文提出了一種直接從幾何約束(包括長度約束和高度約束)出發(fā)的設(shè)計(jì)方法,該設(shè)計(jì)方法結(jié)合了最大推力理論,可以使噴管在滿足幾何約束的同時(shí)優(yōu)化其氣動性能。一般來說,進(jìn)口條件和高度約束可以基本確定噴管內(nèi)氣流膨脹程度。但是為了進(jìn)一步有效控制噴管幾何尺寸并提升氣動性能,本設(shè)計(jì)方法引入2個(gè)設(shè)計(jì)因子,即比例因子和非對稱因子,用于實(shí)現(xiàn)對噴管設(shè)計(jì)的主動控制,以滿足實(shí)際應(yīng)用過程中的不同性能參數(shù)需求。其中,比例因子用于確定核心點(diǎn)的設(shè)計(jì)參數(shù),從而實(shí)現(xiàn)噴管內(nèi)氣流膨脹程度可控,進(jìn)一步提升噴管氣動性能。非對稱因子則用于控制噴管上、下壁面的非對稱膨脹程度,以平衡推升力。圖4為本文提出的設(shè)計(jì)方法的流程圖。具體設(shè)計(jì)過程如下。
圖4 噴管設(shè)計(jì)流程
步驟1根據(jù)噴管的進(jìn)口參數(shù)和高度約束,利用一維流理論確定完全膨脹下噴管的出口馬赫數(shù),并將其作為噴管設(shè)計(jì)過程中的初始馬赫數(shù)Ma0:
q(Ma0)Ae=q(Main)Ain
(9)
(10)
式中:Ma0為初始馬赫數(shù);Main為噴管進(jìn)口馬赫數(shù);Ae為噴管出口截面面積;Ain為噴管進(jìn)口截面面積;k為氣體比熱比。
步驟2確定比例因子β,利用式(11)和等熵關(guān)系式,獲得噴管的設(shè)計(jì)馬赫數(shù)Mad以及其他氣動參數(shù)。
Mad=β·Ma0
(11)
步驟3確定上壁面初始膨脹角δu,并利用式(12)得到下壁面初始膨脹角δd。根據(jù)噴管進(jìn)口參數(shù)分布,初始膨脹段以及下壁面直線EE′,利用特征線法便可以求解區(qū)域①、②、③和④,并得到上下膨脹段發(fā)出的最后特征線的交點(diǎn),即核心點(diǎn)I。其中,區(qū)域①的求解過程見圖3(a)。由特征線法可知,對于一條非特征線的初值線,可以采用單元過程中的內(nèi)點(diǎn)法求解其影響域。在本設(shè)計(jì)方法中,進(jìn)口條件確定了初值線AB,因此利用內(nèi)點(diǎn)法對區(qū)域①進(jìn)行求解,并獲得線AC上的流動參數(shù)。然后,根據(jù)式(13)和δu確定上初始膨脹段AD,并進(jìn)一步采用壁面點(diǎn)法即可獲得線AD上各節(jié)點(diǎn)的流動參數(shù)。在此基礎(chǔ)上,利用內(nèi)點(diǎn)法求解區(qū)域②,如圖3(b)所示。區(qū)域③和④的求解同理。需要注意的是,線EE′為一條可變直線,其長度在核心點(diǎn)I的迭代過程中確定。
(12)
y=b1x2+b2
(13)
式中:φ為非對稱因子,表征上、下壁面膨脹程度的非對稱性;x和y分別為節(jié)點(diǎn)的橫、縱坐標(biāo);b1和b2為常數(shù)。
步驟4為了使核心點(diǎn)I的氣動參數(shù)滿足步驟2中的設(shè)計(jì)要求,需要對步驟3進(jìn)行迭代,得到滿足式(14)的上初始膨脹角δu和下壁面直線EE′的長度ld。進(jìn)一步獲得區(qū)域①、②、③、④的流動參數(shù)。
(14)
步驟5在上膨脹段AD發(fā)出的最后特征線DI上構(gòu)建上位關(guān)鍵點(diǎn)Ku的求解空間。區(qū)域⑤的求解過程如圖3(c)所示。利用最大推力理論確定沿著上位關(guān)鍵點(diǎn)Ku發(fā)出的左行特征線KuM上的流動參數(shù),進(jìn)一步借助變向點(diǎn)法求解區(qū)域⑤。此外根據(jù)流量守恒原理,流經(jīng)特征線KuM的流量與通過特征線DKu的流量相等,并以此逐個(gè)確定構(gòu)成型線DM的各點(diǎn)參數(shù),獲得壁面型線AM。
需要注意的是,此時(shí)的壁面AM未必滿足長度約束,因此在求解空間DI中,采用二分法尋找到合適的Ku點(diǎn)以保證噴管的長度在約束范圍內(nèi),并求解出待設(shè)計(jì)的推力噴管上壁面AM。
步驟6類似步驟5,采用相同方法,根據(jù)高度約束確定下位關(guān)鍵點(diǎn)Kd并獲得待設(shè)計(jì)的推力噴管下壁面BN,完成噴管設(shè)計(jì)。
本文采用計(jì)算流體力學(xué)方法(Computational Fluid Dynamics, CFD)來獲得噴管的流場結(jié)構(gòu)和氣動性能,數(shù)值方法中的控制方程為二維、穩(wěn)態(tài)、可壓縮的雷諾平均Navier-stokes方程,利用有限體積法對控制方程進(jìn)行離散,控制面上的無黏通量采用Roe-FDS通量差分格式獲取,湍流模型使用Shear Stress Transport(SST)k-ω兩方程模型。為了準(zhǔn)確獲取壁面處的流動信息,對壁面附近網(wǎng)格進(jìn)行加密,壁面處y+約為1。在數(shù)值模擬過程中,采用了壓力入口、壓力出口和壓力遠(yuǎn)場等邊界條件。
為驗(yàn)證上述數(shù)值方法的有效性,對文獻(xiàn)[13]中的噴管利用上述方法進(jìn)行數(shù)值模擬校核。該噴管的幾何構(gòu)型和計(jì)算網(wǎng)格如圖5所示,其中模型的幾何尺寸使用噴管的喉部高度ht進(jìn)行無量綱化。根據(jù)文獻(xiàn)[13]中的實(shí)驗(yàn)工況確定數(shù)值模擬中進(jìn)出口邊界條件。數(shù)值模擬得到的噴管壁面壓力與實(shí)驗(yàn)結(jié)果(EXP)對比如圖6所示,并利用進(jìn)口總壓p*對壓力分布進(jìn)行無量綱化。從圖中可以看到該數(shù)值方法能夠準(zhǔn)確地獲得噴管內(nèi)的壓力變化,壁面壓力的平均相對偏差僅為4.9%,因此可以證明上述數(shù)值方法的有效性和準(zhǔn)確性。
圖5 噴管計(jì)算網(wǎng)格
圖6 噴管上壁面壓力分布
為消除計(jì)算網(wǎng)格分辨率對數(shù)值模擬準(zhǔn)確性的影響,本節(jié)開展了網(wǎng)格無關(guān)性研究。噴管構(gòu)型和遠(yuǎn)場設(shè)置分別如圖7中的紅色區(qū)域和綠色區(qū)域所示。對該幾何模型使用商用軟件ICEM進(jìn)行網(wǎng)格劃分,獲得了網(wǎng)格尺寸為粗(Coarse)、中(Middle)、細(xì)(Fine)尺度的計(jì)算網(wǎng)格(網(wǎng)格量分別為7萬、15萬和30萬),其中,中尺度網(wǎng)格見圖7。相同工況下,不同尺度網(wǎng)格得到的上壁面壓力分布如圖8所示,可以看到整體上3種尺度網(wǎng)格的壓力分布基本一致。但通過局部放大,發(fā)現(xiàn)相比中尺度和細(xì)尺度的網(wǎng)格,粗尺度網(wǎng)格在x≈0.3, 0.75 m 處壓力分布出現(xiàn)較大偏差,最大偏差可達(dá)1.72%。因此,說明中尺度網(wǎng)格可以滿足計(jì)算精度的要求,在下文中均采用中等尺度的計(jì)算網(wǎng)格開展研究。
圖7 目標(biāo)噴管的計(jì)算域及網(wǎng)格劃分
圖8 不同尺度網(wǎng)格的噴管上壁面壓力分布
為驗(yàn)證本文提出的幾何尺寸約束的超燃沖壓發(fā)動機(jī)推力噴管設(shè)計(jì)方法的有效性,本節(jié)運(yùn)用該方法在某典型工況下對噴管進(jìn)行設(shè)計(jì),并對設(shè)計(jì)得到的噴管進(jìn)行無黏模擬驗(yàn)證。設(shè)計(jì)工況和幾何約束如表1所示,其中Lc為長度約束,Hc為高度約束,pin為進(jìn)口壓力,Tin為進(jìn)口溫度,Main為進(jìn)口馬赫數(shù),H為飛行高度,Mac為巡航馬赫數(shù)。此外本驗(yàn)證算例中β取1.02,φ取0.9。設(shè)計(jì)得到的噴管構(gòu)型見圖9,可以看到噴管的幾何尺寸能夠滿足約束條件。圖9為分別通過MOC和CFD獲得的噴管在設(shè)計(jì)工況下運(yùn)行時(shí)的壓力云圖,可以看到MOC的壓力分布與CFD結(jié)果基本一致,同時(shí)MOC很好地捕捉到了噴管出口附近下壁面的反射膨脹波。此外,MOC獲得的壁面壓力相對CFD壁面壓力結(jié)果的平均誤差僅為1.36%。由此說明該方法能夠有效地設(shè)計(jì)噴管內(nèi)流場結(jié)構(gòu),并獲得滿足幾何約束的噴管型線。
表1 設(shè)計(jì)工況
圖9 通過MOC和CFD獲得的噴管壓力云圖
為研究在設(shè)計(jì)過程中比例因子對噴管型線以及氣動性能的影響,本節(jié)使用3個(gè)不同的比例因子在相同工況和幾何約束下對噴管進(jìn)行設(shè)計(jì),并獲得噴管的流場結(jié)構(gòu)和氣動性能。設(shè)計(jì)得到的噴管及其流場結(jié)構(gòu)如圖10所示,可以看到比例因子會顯著影響噴管構(gòu)型,進(jìn)而改變噴管內(nèi)流場結(jié)構(gòu)。隨著比例因子增加,噴管的擴(kuò)張程度逐漸增加,導(dǎo)致上壁面高度增加,同時(shí)由于高度約束的存在,使下壁面逐漸縮短。此外,噴管出口馬赫數(shù)也隨著比例因子的增加而增加。
圖10 不同比例因子的噴管馬赫數(shù)云圖
表2給出了由不同比例因子設(shè)計(jì)得到的噴管的氣動性能,表中Fx為噴管軸向推力,L為升力,M為俯仰力矩,俯仰力矩的參考中心取在噴管進(jìn)口下壁面處,Cfx為噴管推力系數(shù),具體計(jì)算式為
(5)地表氧化鐵帽標(biāo)志:礦體經(jīng)風(fēng)化后,呈黃褐色,具蜂窩狀構(gòu)造,可見有大量的褐鐵礦斑點(diǎn)(鐵帽),是找礦的直接標(biāo)志。
(15)
(16)
(17)
圖11 不同比例因子的噴管氣動性能分布
為研究非對稱因子對噴管設(shè)計(jì)及性能的影響,本節(jié)使用了3個(gè)不同的非對稱因子分別設(shè)計(jì)噴管,比例因子均取本研究中推力最優(yōu)的比例因子,即β=1.02。通過數(shù)值模擬可以得到相應(yīng)的噴管性能和流場結(jié)構(gòu)。由圖12可以發(fā)現(xiàn)噴管下壁面會隨著非對稱因子的增加而縮短,下壁面的縮短意味著推力和負(fù)升力的有效作用面積減小,這將導(dǎo)致推力降低同時(shí)升力增加;此外,噴管上壁面擴(kuò)張程度出現(xiàn)輕微減小,導(dǎo)致氣流膨脹程度降低,從而使上壁面壓力增加,這有利于推、升力的增加。
圖12 不同非對稱因子的噴管馬赫數(shù)云圖
不同非對稱因子的噴管氣動性能見表3。噴管氣動性能隨非對稱因子的變化趨勢如圖13 所示。可以看到隨著非對稱因子逐漸增加,噴管的推力系數(shù)逐漸降低,在φ為1.0處達(dá)到最小值0.741, 同時(shí)噴管的升力和俯仰力矩隨著非對稱因子的增加而增大,在φ為1.0處達(dá)到最大值,分別為3 009.97 N和2 154.46 N·m,這與上述分析結(jié)果一致。
表3 不同非對稱因子的噴管氣動性能
圖13 不同非對稱因子的噴管氣動性能分布
以上研究表明,通過調(diào)整比例因子和非對稱因子能夠?qū)崿F(xiàn)對噴管氣動性能的主動控制,在現(xiàn)實(shí)應(yīng)用中,可以根據(jù)對噴管性能的實(shí)際需求選取合適的比例因子和非對稱因子。
為證明本文提出的設(shè)計(jì)方法的優(yōu)越性,本節(jié)將其與典型的噴管截短設(shè)計(jì)方法[13]進(jìn)行對比分析。典型截短設(shè)計(jì)方法示意圖如圖14所示,其根據(jù)進(jìn)出口條件設(shè)計(jì)完全膨脹噴管,并利用等比例截短方法使噴管滿足幾何約束,具體設(shè)計(jì)過程可見文獻(xiàn)[13]。
圖14 典型截短設(shè)計(jì)方法示意圖
應(yīng)用截短設(shè)計(jì)方法設(shè)計(jì)噴管時(shí),由于無法通過截短方法一次性獲得滿足高度和長度約束的噴管,因此首先根據(jù)長度約束進(jìn)行等比例截短獲得噴管6,其幾何構(gòu)型和流場結(jié)構(gòu)可見圖15,具體尺寸在表4中給出,其中Ln為噴管長度,Hn為噴管高度??梢钥吹絿姽?的高度并不能滿足約束條件。其次,根據(jù)高度約束進(jìn)行等比例截短獲得噴管7,可以發(fā)現(xiàn)此時(shí)噴管7的長度會小于長度約束。這從側(cè)面再次說明截短設(shè)計(jì)方法并不能總是使噴管完全滿足幾何約束,并且可能會使部分尺寸過度截短,導(dǎo)致噴管性能惡化。
圖15 傳統(tǒng)噴管的馬赫數(shù)云圖
根據(jù)2種截短條件得到的噴管性能參數(shù)在表4 中給出,同時(shí)給出了噴管2(由本文設(shè)計(jì)方法設(shè)計(jì)得到)的性能參數(shù)進(jìn)行對比,表中ΔXn的定義為
表4 設(shè)計(jì)工況下噴管性能對比
(18)
式中:X為噴管的氣動性能參數(shù),可分別取Cfx、L和M,下標(biāo)“n”為噴管序號??梢钥吹礁鶕?jù)長度約束得到的噴管6的推力性能高于根據(jù)高度約束得到的噴管7的推力性能,這說明隨著截短程度的增加,噴管的性能損失也隨之增加。此外,相對噴管6,本文提出的設(shè)計(jì)方法得到的噴管2的推力系數(shù)和升力分別提高了11.95%和209.50%,俯仰力矩下降了3.02%。需要注意的是,噴管6的高度并不滿足約束條件。而相對可以滿足幾何約束的噴管7,噴管2的性能提升更為顯著,其推力系數(shù)、升力和俯仰力矩的增幅分別可達(dá)33.36%、265.75%和37.21%。
因此,本文提出的噴管設(shè)計(jì)方法不僅可以使噴管滿足幾何約束條件,并且可以有效提升噴管的氣動性能。此外,通過調(diào)整比例因子和非對稱因子可以實(shí)現(xiàn)對噴管推力和升力性能的控制,以滿足實(shí)際應(yīng)用中的需求。
本文提出了一種基于最大推力理論和特征線方法的幾何尺寸約束的超燃沖壓發(fā)動機(jī)噴管設(shè)計(jì)方法,該設(shè)計(jì)方法能夠在滿足尺寸約束的前提下優(yōu)化噴管氣動性能,并且可以通過調(diào)節(jié)比例因子和非對稱因子來實(shí)現(xiàn)對噴管推、升力的控制。本文的主要結(jié)論如下:
1) 隨著比例因子增加,噴管推力系數(shù)先增加再降低,升力和俯仰力矩逐漸增加,對于特定的設(shè)計(jì)工況和幾何約束存在推力最佳的比例因子。隨著非對稱因子增加,噴管推力系數(shù)降低,升力和俯仰力矩逐漸增加,在非對稱因子為1.0時(shí),噴管推力系數(shù)達(dá)到最小值0.741,升力和俯仰力矩達(dá)到最大值,分別為3 009.97 N和2 154.46 N·m。
2) 相比典型的噴管截短設(shè)計(jì)方法,本文提出的設(shè)計(jì)方法不僅可以保證噴管滿足幾何約束,同時(shí)還能使其性能得到顯著提升,噴管的推力系數(shù)、升力和俯仰力矩最高可提升33.36%、265.75%和37.21%。此外,在實(shí)際應(yīng)用過程中可以根據(jù)不同的氣動性能需求來調(diào)整比例因子和非對稱因子,以獲得滿足不同性能要求的噴管。