張 強，楊 永

（西北工業(yè)大學(xué) 翼型、葉柵空氣動力學(xué)國防科技重點實驗室，陜西 西安 710072）

0 引 言

湍流是自然界中普遍存在的一種流動形態(tài)。在宏觀上，湍流表現(xiàn)為流體微團的不規(guī)則運動，湍流運動產(chǎn)生的質(zhì)量和能量的輸運遠大于分子熱運動所產(chǎn)生的輸運，湍流運動引起的質(zhì)量和能量的平均擴散是不能忽略的。因此湍流研究一直是流體力學(xué)研究中的一個重要方向。

國內(nèi)外針對平板、翼型、機翼等外流標(biāo)準(zhǔn)算例，開展的數(shù)值模擬研究已經(jīng)很多，但針對噴管內(nèi)流的湍流模型研究主要以Spalart－Allmaras（SA）一方程模型或者改進的k－ε兩方程模型為主，其它湍流模型的驗證研究還比較少見?，F(xiàn)今兩方程湍流模型的研究按研究角度主要分為三類：構(gòu)造k－ε和k－ω的混合模型；構(gòu)造更為精確的剪應(yīng)力構(gòu)成模型；構(gòu)造新的更準(zhǔn)確的耗散方程。其代表模型分別是：SSTk－ω兩方程模型，EASMk－ω模型和基于新的耗散方程的k－ζ模型。

國外對射流推力矢量噴管已進行了系統(tǒng)的試驗研究［1］，形成了較為完整的數(shù)據(jù)庫，為CFD數(shù)值模擬研究提供了良好的條件。因此本文的研究目的就是針對激波控制射流推力矢量噴管的內(nèi)外流場、噴管性能參數(shù)等，分別采用SSTk－ω，EASMk－ω和k－ζ模型，進行系統(tǒng)的驗證分析，并與SA模型的計算結(jié)果和試驗數(shù)據(jù)進行對比，考察不同湍流模型對噴管內(nèi)外流場的流場結(jié)構(gòu)、噴管性能參數(shù)以及主噴管內(nèi)壁面的壓力分布等的預(yù)測能力。之所以選擇激波控制射流推力矢量噴管作為驗證算例，主要是因為其流場中由于二次射流的干擾，會在主噴流中誘導(dǎo)產(chǎn)生斜激波，存在較強的激波／附面層干擾，這樣的流動對目前湍流模型的研究也是一種極大的挑戰(zhàn)，因此其數(shù)值模擬結(jié)果也更具說服力。

1 計算模型及網(wǎng)格

研究激波控制推力矢量噴管時，所采用的模型是NASA LaRC試驗用的噴管，設(shè)計主流壓比為NPRD＝8.78，其幾何形狀及尺寸見參考文獻［2］，單縫引射對應(yīng)外形1，雙縫引射對應(yīng)外形5。

單縫引射的推力矢量噴管，計算網(wǎng)格如圖1所示，圖中生成的計算網(wǎng)格共包括5塊，包圍主噴管的箱體上下各有一塊網(wǎng)格，主噴管和次流噴管各生成一塊網(wǎng)格，主噴管的出口及箱體上下兩塊網(wǎng)格的下游，點對點生成一塊網(wǎng)格，以計算噴管的尾跡流場。雙縫引射噴管外部網(wǎng)格拓?fù)浣Y(jié)構(gòu)與單縫引射噴管相同，只是在主噴管擴張段上多了另外一個小噴口的網(wǎng)格，如圖2所示，并且兩個噴口總的開口面積和單縫時的開口面積相同。計算中參考長度取為從主噴管出口到其喉道的長度c，由于噴管總長度為2c，因此后面的噴管無量綱化的長度為2，這樣處理以方便與試驗值的比較。

圖1 單縫引射噴管的計算網(wǎng)格及噴口局部圖Fig.1 Computational mesh and close－up of one injection port

圖2 雙縫引射噴管的噴口局部圖Fig.2 Close－up of two injection ports

2 數(shù)值計算方法

針對射流推力矢量噴管的數(shù)值模擬，采用的是課題組自有的研究型軟件。假定射流推力矢量噴管內(nèi)外流動是定常的，通過求解不同湍流模型封閉的RANS方程，模擬射流推力矢量噴管的內(nèi)外流場。RANS流動控制方程的無粘項采用Roe－FDS格式，粘性項采用中心差分，湍流模型的空間項離散也采用二階迎風(fēng)格式，流動主控方程和湍流模型的控制方程均采用近似因子隱式時間推進，以提高數(shù)值計算效率和穩(wěn)定性，關(guān)于算法的詳細介紹可以參考文獻［3－5］。

計算條件取為：單縫引射時外流Ma＝0.05，雙縫引射時Ma＝0.1，雷諾數(shù)Re＝1.65×105，分別采用SA、SSTk－ω、EASMk－ω和k－ζ四種湍流模型求解了主流壓比NPR＝4.6，7.0，8.78，10.0和SPR＝0.7、1.0時的射流推力矢量噴管的內(nèi)外流場。靜壓固定為標(biāo)準(zhǔn)海平面大氣壓，總溫固定為300K。

3 單縫引射噴管的數(shù)值模擬

利用文獻［2］的試驗結(jié)果首先考察了外流馬赫數(shù)對數(shù)值模擬結(jié)果的影響，在兩個外流馬赫數(shù)Ma＝0.05及Ma＝0.1下，采用SA湍流模型，模擬了NPR＝4.6、10.0，SPR＝0.7時單縫引射的噴管內(nèi)外流場。計算的主噴管內(nèi)壁面壓力分布與在靜止大氣中進行的試驗值的比較如圖3（a、b）所示，總的來看，在兩個馬赫數(shù)下的數(shù)值模擬結(jié)果與試驗值比較吻合。在兩個較小的外流馬赫數(shù)下，只在NPR＝4.6時，噴管出口截面處的壓強略有微小差別，并且我們的數(shù)值驗證也表明，在這兩個馬赫數(shù)下，噴管性能參數(shù)的差別也可以忽略不計，考慮到自由來流馬赫數(shù)比較小時（例如Ma＝0.05或0.01），流動控制方程不易收斂，為了提高數(shù)值計算的穩(wěn)定性和效率，自由來流馬赫數(shù)也可以取0.1，對數(shù)值模擬結(jié)果的影響不大，可以認(rèn)為它們是相同的。

Ma＝0.05，NPR＝4.6，SPR＝0.7，SA和SSTk－ω模型給出的噴管流場及二次射流局部流場如圖4所示，由于二次射流的干擾，在次流噴口之前誘導(dǎo)出現(xiàn)斜激波，在次流噴口上游和下游區(qū)域均會出現(xiàn)較大的分離區(qū)，兩種湍流模型都給出了與試驗很接近的激波結(jié)構(gòu)（圖5），從EASMk－ω和k－ζ模型的計算結(jié)果也可以得出相同的結(jié)論。四種模型計算給出的主噴管上壁面，次流引起的斜激波位置都比試驗值稍微靠后，約在x／c＝1.62處，而試驗測量則是在x／c＝1.543處。

圖3 SPR＝0.7，Ma＝0.05和0.1時，SA模型計算的主噴管內(nèi)壁面壓力分布與試驗值的比較Fig.3 Experimental and computational pressure distribution at Ma＝0.05，0.1with SPR＝0.7，SA model

圖4 SA和SST k－ω模型計算馬赫云圖和次噴口局部流場（NPR＝4.6，SPR＝0.7，Ma＝0.05）Fig.4 Mach counters and flow pattern near the injection port for NPR＝4.6，SPR＝0.7，Ma＝0.05with SA and SST model

圖5 NPR＝4.6，SPR＝0.7時試驗陰影照片F(xiàn)ig.5 Experimental shadowgraph for one injection port，NPR＝4.6，SPR＝0.7

Ma＝0.05，SPR＝0.7，NPR分別為4.6、7.0、8.78、10.0時，主噴管內(nèi)壁面的壓力分布與試驗值的對比如圖6（a～d）所示。在主噴管下壁面，當(dāng)NPR＝4.6時，有一段斜激波在噴管內(nèi)部，SA模型計算的激波位置最靠近出口，其它三種模型給出的激波位置差別不大，但都優(yōu)于SA模型的計算結(jié)果，更接近實驗值，從喉道到主噴管出口方向依次是k－ζ，SSTk－ω，EASMk－ω模型，在其他NPR下，下壁面不存在斜激波，四種湍流模型計算給出的壓力分布相同，與試驗值吻合良好。

在主噴管上壁面，由于二次射流的影響，在次流噴口的上游和下游位置都會出現(xiàn)較大的分離區(qū)，二次噴流誘導(dǎo)產(chǎn)生的斜激波位于上游分離區(qū)之前，更接近喉道處，在此激波之前主噴管的上壁面，不同NPR下各湍流模型給出的壓力分布相同，與試驗值完全吻合。在斜激波和次流噴口之間區(qū)域，不同NPR下，四種湍流模型給出的激波位置都比試驗值略靠后，但壓力分布差別不大，僅略低于試驗值。在次流噴口以后的區(qū)域，受外部回流分離區(qū)的影響，不同湍流模型給出的壓力分布略有不同，NPR越大，回流分離區(qū)越小，各湍流模型的預(yù)測越準(zhǔn)確。在圖6（c）中，當(dāng)NPR＝8.78時，在主噴管上壁面次流噴口以后的區(qū)域，EASMk－ω的結(jié)果與試驗值差別較大，甚至還不如SA模型，這可能與我們沒有考慮流動的非定常效應(yīng)有關(guān)。可見，在有較大分離的區(qū)域，四種湍流模型的預(yù)測能力還是有差別的。從總的壓力分布來看，相比較而言，SSTk－ω模型最 準(zhǔn)確，其次是k－ζ和EASMk－ω，SA模型的誤差最大。

在SPR＝0.7時，分別采用SA、SSTk－ω、EASMk－ω、k－ζ四種湍流模型計算的射流推力矢量噴管性能參數(shù)隨NPR的變化如圖7所示，不同湍流模型計算給出的噴管性能參數(shù)差別不大，隨NPR的增加，流量系數(shù)幾乎保持一個常數(shù)，推力系數(shù)是先增加后減小，矢量偏角和矢量效率都是一直減小的，與試驗得出的規(guī)律一致［2］。矢量偏角和矢量效率隨NPR的增加減小的原因是：由于二次射流與主噴管流動相互干擾，誘導(dǎo)生成斜激波，隨著NPR的增加，受主噴流的影響，斜激波末端逐漸外移，遠離主噴管出口，并且強度變?nèi)?，?dǎo)致產(chǎn)生的矢量偏角變小和推力矢量效率降低。

圖6 不同NPR下各湍流模型計算的主噴管壓力分布（SPR＝0.7）Fig.6 Pressure distribution of the primary nozzle wall for different NPR（SPR＝0.7）

圖7 4種湍流模型計算的噴管性能參數(shù)（SPR＝0.7）Fig.7 Computational nozzle performance with four turbulence models（SPR＝0.7）

圖8是SSTk－ω模型在NPR＝4.6，SPR＝0.7和1.0時計算給出的主噴管上壁面的壓力分布，可以看出，在相同NPR下，SPR越大，次流噴口前的激波強度越強，位置越靠近喉道，即隨SPR的增加，激波前移，強度增加。不同湍流模型的計算結(jié)果均表明，從次噴口進入主噴口的流量與NPR無關(guān)，只隨SPR的增加而增加。在不同NPR下，SPR＝0.7時，四種湍流模型計算的次噴流流量約是主噴流的4.2%，試驗是4%，SPR＝1.0時，次噴流流量約是主噴流的5.7%，試驗是6%。

圖8 不同SPRSST k－ω模型計算的上壁面壓力分布（NPR＝4.6）Fig.8 Pressure distribution of the primary nozzle upper wall with SST k－ω model，NPR＝4.6，SPR＝0.7，1.0

此外，四種湍流模型在SPR＝0.7時，數(shù)值計算給出的推力矢量偏角隨NPR的變化與試驗值的比較如圖9所示，各湍流模型的計算結(jié)果與試驗值吻合良好。在NPR大于4之后，隨著NPR的增加，推力矢量偏角變小，與試驗得出的結(jié)論一致。

4 雙縫引射噴管的數(shù)值模擬

圖9 推力矢量偏角隨NPR變化與試驗值的比較（SPR＝0.7）Fig.9 SPR＝0.7，Computational thrust vector angle and experiments varying with NPR

雙縫引射的噴管，在計算過程中取外流馬赫數(shù)為0.1，正如我們在前文中所指出的，當(dāng)馬赫數(shù)不大于0.1時，外流馬赫數(shù)對噴管性能的影響很小，可以忽略不計。由于試驗中兩個次流注入縫是從同一氣源引氣，因此在數(shù)值計算過程中兩個次流噴口的入口給定總溫、總壓邊界條件設(shè)置相同（即SPR相同），其它邊界條件的設(shè)定與單縫引射時完全一樣。

NPR＝4.6，SPR＝0.7，不同湍流模型給出的噴管及尾跡流場如圖10所示，雙縫引射會在主噴管中誘導(dǎo)產(chǎn)生兩道斜激波，第一道斜激波位于噴口2（圖2）之前，第二道斜激波在兩個噴口之間，兩個噴口左右兩側(cè)都會出現(xiàn)分離區(qū)，但噴口1兩側(cè)的分離區(qū)明顯大于噴口2，各湍流模型都捕捉到了試驗所觀測到的激波結(jié)構(gòu)（圖11）。四種模型計算給出的主噴管內(nèi)第一道斜激波的位置都比試驗值靠后，約x／c＝1.4處，而試驗測量則是在x／c＝1.24處，第二道激波位置與試驗相比，k－ζ模型計算給出的激波位置偏靠前，SSTk－ω和EASMk－ω模型與試驗值頗為吻合，SA模型給出的激波位置偏靠后，下壁面出口附近的斜激波試驗測量在x／c＝1.89處，三種兩方程湍流模型的計算結(jié)果比較接近，都約在x／c＝1.9處，SA最靠后，約在x／c＝1.94處。

雙縫引射的推力矢量噴管，SPR＝0.7時，在一系列NPR下，不同湍流模型計算的主噴管內(nèi)壁面壓力分布與試驗值的比較如圖12（a～d）所示，總的來看均與試驗值吻合良好。雙縫引射時第一道斜激波的壓力曲線比試驗值更陡，激波位置比試驗值偏靠后。第二道斜激波處的壓力分布比試驗值要低，但對激波位置的捕捉要好于第一道斜激波。總體來說，kζ模型計算給出的激波位置最為靠前，SSTk－ω模型次之，之后是EASMk－ω模型，SA模型最為靠后。

圖10 不同模型計算的馬赫云圖（NPR＝4.6，SPR＝0.7）Fig.10 Mach counters for two injection ports nozzle at NPR＝4.6，SPR＝0.7，with different turbulence models

圖11 NPR＝4.6，SPR＝0.7時試驗陰影照片F(xiàn)ig.11 Experimental shadowgraph for two injection ports，NPR＝4.6，SPR＝0.7

圖12 雙縫引射噴管上下壁面壓力分布與試驗值對比Fig.12 Pressure distribution of the primary nozzle wall with two injection ports at different NPR

雙縫引射時，不同湍流模型計算的流量系數(shù)，推力系數(shù)隨NPR、SPR的變化規(guī)律與單縫引射相同，在此不再贅述。

圖13是SPR＝0.7，不同湍流模型計算的雙縫引射的推力矢量偏角隨NPR的增加與試驗值的比較，從圖中可以看出，各湍流模型都是在NPR＝4.6時給出了最大的推力矢量偏角，然后隨NPR增加，矢量偏角減小。在NPR＝7.0時，SSTk－ω和k－ζ模型計算的推力矢量偏角稍偏大，主要是由于這兩種湍流模型過高估計了噴口1后回流分離區(qū)的壓力，如圖12（b）所示，可能是由于此狀態(tài)下，回流分離區(qū)的非定常效應(yīng)較強，采用定常方法計算，這兩種模型無法給出準(zhǔn)確的結(jié)果。圖14是NPR＝4.6，不同SPR下，SSTk－ω模型給出的噴管上壁面壓力分布，從圖中可以看出，與單縫引射類似，隨SPR的增加，兩道斜激波的強度增加，位置前移。

圖13 推力矢量偏角隨NPR的變化（SPR＝0.7）Fig.13 SPR＝0.7，Computational thrust vector angel and experiments varying with NPR

圖14 不同SPR下SST k－ω模型計算的上壁面壓力分布Fig.14 Pressure distribution of the primary nozzle upper wall with two injection ports，SST k－ω model

5 單／雙縫引射矢量噴管性能比較

當(dāng)SPR＝0.7時，SA湍流模型計算的單／雙縫引射噴管性能參數(shù)對比如圖15所示。雙縫引射時，盡管兩個次流注入縫總的開口面積與單縫引射噴管相同，SPR也相同，仍然會導(dǎo)致噴管的流量系數(shù)降低約1%，在相同的NPR下，推力系數(shù)比單縫引射增加了大約1%，利用其它湍流模型的計算結(jié)果也可以得出同樣的結(jié)論。

盡管兩個次流注入縫的入口及主噴管入口邊界條件設(shè)置完全相同（NPR、SPR相同），但由于次流注入縫在主噴管的位置不同，進入主噴管的流量是不同的，靠近主噴管喉道的噴口（噴口2）要比靠近主噴管出口的噴口（噴口1）的流量稍小一些，進入主噴管的總流量也比單縫引射噴管少，減小了約0.7%。

NPR＝4.6，SPR＝0.7時，不同湍流模型給出了在6.4°～6.8°范圍內(nèi)的矢量偏角（試驗是6.5°），小于單次流噴口時不同湍流模型給出的介于7.3°～7.6°之間的推力矢量偏角（試驗是7.0°），但由于雙縫引射時進入主噴管的流量較小，每1%次流注入產(chǎn)生的推力矢量偏角即矢量效率反而是雙縫引射大于單縫引射的噴管。在其它NPR下，雙縫引射產(chǎn)生的推力矢量效率卻比單縫引射小，如圖15（d）所示。

圖15 SA模型計算的單／雙縫引射噴管性能參數(shù)比較（SPR＝0.7）Fig.15 Computational nozzle performance with SA turbulence model（SPR＝0.7）

6 結(jié) 論

分別采用四種不同湍流模型，通過對激波控制單／雙縫引射推力矢量噴管的數(shù)值模擬研究分析，可以得出以下結(jié)論：

（1）四種湍流模型都能比較準(zhǔn)確的預(yù)測出單／雙縫射流推力矢量噴管的性能參數(shù)，并且計算結(jié)果差別不大，與試驗得出的規(guī)律一致。

（2）各湍流模型計算給出的單／雙縫主噴管內(nèi)壁面壓力分布均與試驗值吻合良好。在激波位置捕捉方面，SSTk－ω、k－ζ和 EASMk－ω模型較好，SA 模型稍差。針對雙縫引射噴管，在回流分離區(qū)內(nèi)，SSTkω、k－ζ模型的計算結(jié)果并不理想，會過高地估計壓力分布。

（3）無論是單縫還是雙縫引射，各湍流模型的計算結(jié)果均表明，相同NPR，如果SPR增加，進入主噴管的次流流量增加，次流誘導(dǎo)的斜激波前移，主噴管流量系數(shù)損失增加，推力系數(shù)損失也增加。

（4）雙縫引射的推力矢量效率并不一定優(yōu)于單縫引射。如文中噴管構(gòu)型，NPR＝4.6時，雙縫引射的矢量效率較高，NPR大于4.6時，單縫引射的推力矢量效率較高。

［1］DEERE K A.Summary of fluidic thrust vectoring research conducted at NASA langley research center［R］.AIAA－2003－3800，2003.

［2］KENRIC A W，DEERE K A.Experimental and computational investigation of multiple injection ports in a convergent－divergent nozzle for fluidic thrust vectoring［R］.AIAA－2003－3802，2003.

［3］ZHANG Q，YANG Y。On Exploring predicting precision of an EASMk－ωtwo－equation turbulence model［J］.JournalofNorthwesternPolytechnical University，2005，23（1）：89－92

［4］ZHANG Q，YANG Y.Applying robinson's new twoequationturbulencek－ζmodel to numerically simulating turbulent flow［J］.JournalofNorthwesternPolytechnicalUniversity，2009，27（4）：466－470.

［5］ZHANG Q，YANG Y.Numerical analysis of low reynolds number flow over airfoil［J］.ActaAerodynamic Sinica，2006，24（4）：482－486.