黨進鋒，郜 冶，劉平安

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，哈爾濱 150001)

可壓縮流渦流冷壁發(fā)動機數(shù)值模擬

黨進鋒，郜 冶，劉平安

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，哈爾濱 150001)

為深入探索入射傾角、入射壓降、發(fā)動機圓柱段長度等參數(shù)與渦流冷壁發(fā)動機內(nèi)雙向渦流結(jié)構(gòu)特征量之間的內(nèi)在關(guān)系，以三維物理模型為基礎(chǔ)，進行了可壓縮流渦流冷壁發(fā)動機冷流數(shù)值模擬。在前期湍流模型選擇中，比較了RSM模型、帶旋流修正RNGk-ε模型、RNGk-ε模型3種模型的可信度，從準確度和經(jīng)濟性兩方面綜合考慮選擇帶旋流修正RNGk-ε模型作為計算模型。研究發(fā)現(xiàn)，切向速度及最大切向速度均隨入射傾角的增大而減小，隨入射壓降增大而增大，隨長徑比的增大而減小。不同入射傾角、不同入射壓降、不同長徑比下，整個發(fā)動機內(nèi)最大切向速度的無量綱徑向位置均恒定在0.19附近。對于不同長徑比工況，長徑比為1.0時，最大切向速度從發(fā)動機頂端到入射口附近逐漸增大，長徑比為1.5時，最大切向速度從發(fā)動機頂端到入射口附近先增大、后減小。切向速度及最大切向速度的軸向衰減率維持在3%以內(nèi)。無量綱渦幔半徑從發(fā)動機頂端開始線性增大到入射口附近，變化范圍為0.71～0.82。入射傾角相同時，隨入射壓降的增大，燃燒室長徑比對最大切向速度大小的影響將隨之增大。

渦流冷壁發(fā)動機；可壓縮；湍流模型；入流條件；長徑比；數(shù)值模擬

0 引言

渦流冷壁發(fā)動機(Vortex Combustion Cold-Wall Chamber,VCCWC)由于內(nèi)部強旋轉(zhuǎn)流動的存在，數(shù)值計算難度相較于無旋轉(zhuǎn)或弱旋轉(zhuǎn)流動增大許多，其計算所需的湍流模型、網(wǎng)格分布、數(shù)值方法等均有別于普通流動。因此，數(shù)值計算研究進展非常有限。美國Majdalani等針對渦流冷壁發(fā)動機內(nèi)雙向渦流進行了大量理論研究[1-3]，其研究成果對于驗證渦流冷壁發(fā)動機可行性、認識雙向渦流結(jié)構(gòu)特點等方面具有非常重要的意義。其與Orbital科技公司Chiaverini等人所做的渦流冷壁發(fā)動機燃燒實驗，為渦流冷壁發(fā)動機研究奠定了基礎(chǔ)[4]。在Majdalani與Chiaverini等人帶動下，國內(nèi)從2010年左右，也開始對渦流冷壁發(fā)動機進行了相關(guān)理論、實驗與數(shù)值計算研究。孫得川等[5]通過理論推導(dǎo)，得到渦流冷壁推力室傳熱模型，通過數(shù)值計算研究，指出一次進氣有利于燃料與空氣的摻混燃燒[6]，文獻[7]通過實驗研究了雙向渦流速度分布規(guī)律。李家文、俞南嘉等人在文獻[8-9]中，通過實驗和數(shù)值計算，驗證了渦流冷卻技術(shù)的可行性，在文獻[10]中指出，冷卻劑噴嘴傾斜一定角度，可改善內(nèi)渦強度，改進燃料摻混、提高燃燒效率。文獻[11]中采用三維全尺寸可壓縮流體模型，驗證了推力室內(nèi)存在同軸反向的雙層渦旋結(jié)構(gòu)，并指出在后續(xù)研究中，有必要明確渦旋速度場環(huán)形區(qū)域?qū)ν屏κ胰紵Ч木唧w影響；文獻中同時還指出，有必要針對渦旋來流條件、圓柱段幾何尺寸、燃料噴注方式等，開展優(yōu)化推力室結(jié)構(gòu)的仿真研究。可看到，現(xiàn)階段對于VCCWC的研究，仍停留在可行性驗證與定性的規(guī)律描述上。因此，針對VCCWC內(nèi)雙向渦流結(jié)構(gòu)特征影響因素及其定量描述的研究就顯得尤為迫切。

本文在文獻[11]研究基礎(chǔ)上，針對文獻[11]中提出的渦旋來流條件、圓柱段幾何尺寸等內(nèi)容，進行了三維可壓縮流數(shù)值模擬研究，所得結(jié)論對于深入認識VCCWC內(nèi)渦流結(jié)構(gòu)隨模型幾何尺寸及來流條件的變化具有重要作用，對于未來VCCWC設(shè)計具有一定指導(dǎo)作用。

1 物理模型、計算工況及邊界條件

由于VCCWC具有軸對稱特征，文中采用與文獻[8]中相同的三維周期性邊界，取全模型周向1/6作為計算模型，該1/6模型上周向均布2個入射口，即全模型為12個入射口。本文為亞音速冷流數(shù)值模擬，流動未達到臨界狀態(tài)，因此采用收縮噴管形式，模型具體尺寸如圖1所示。圖2為三維模型網(wǎng)格示意圖。

表1為本文中所完成的16種計算工況。發(fā)動機圓柱段長度L分別為100 mm和150 mm，對應(yīng)發(fā)動機長徑比αL=L/D分別為1和1.5，發(fā)動機頂部無量綱軸向位置記為α=0，則切向入射處無量綱軸向位置分別為α=1和α=1.5。無量綱徑向位置記為β=r/R，則壁面處β=1，對稱軸處β=0。入射傾斜角θ從0°增大30°，入射壓降Δpin分別為41 700 Pa和62 550 Pa，操作壓力為101 325 Pa，出口為壓力出口，計算工質(zhì)為300 K可壓縮理想氣體。

2 計算可信性驗證

2.1 湍流模型選擇及可信性驗證

文獻[6]中指出，求解渦流冷壁推力室內(nèi)的強旋轉(zhuǎn)流動，采用RSM(雷諾應(yīng)力)模型較臺適，文獻[8]采用RNGk-ε模型進行了燃燒數(shù)值模擬，文獻[11]采用RNGk-ε模型，并激活旋流修正進行了燃燒數(shù)值模擬， Fang在文獻[12]中結(jié)合Standardk-ε模型和RSM模型數(shù)值計算，模擬了VCCWC冷流與燃燒過程。筆者在前期研究中也發(fā)現(xiàn)，RSM湍流模型相較于兩方程模型具有一定優(yōu)勢。兩方程模型中，RNGk-ε模型由于可加入旋流修正，大大改善了其旋轉(zhuǎn)流動數(shù)值計算的準確性。本文在選擇計算所采用的湍流模型時，對RSM、RNGk-ε模型、帶旋流修正的RNGk-ε模型進行了比較，從數(shù)值結(jié)果準確性與計算經(jīng)濟性兩方面，綜合決定最終所采用的湍流模型。計算中，壓力-速度耦合算法為SIMPLE算法，梯度差分采用G-S節(jié)點基格式，壓力差分采用二階格式，其余通量差分均為二階迎風(fēng)格式，計算采用Ansys Fluent 11.0軟件。

表1 計算工況

圖3為3種驗證模型計算所得壓力曲線與文獻[13]中實驗數(shù)據(jù)對比曲線，圖4為3種驗證模型計算所得壓力云圖。從圖3可看到，各湍流模型計算所得燃燒室內(nèi)壓力值變化趨勢并不相同。實驗數(shù)據(jù)在0.6<β<1.0范圍內(nèi)，壓力下降緩慢，β<0.6后，壓力下降速度逐漸增大。3種模型中，RSM模型曲線與實驗數(shù)據(jù)吻合度最好，其次是帶旋流修正的RNGk-ε模型，不帶旋流的RNGk-ε模型相差較多。圖4分別為RSM模型、帶旋流的RNGk-ε模型及不帶旋流的RNGk-ε模型壓力云圖。從圖4同樣可看到，3種模型所得壓力云圖變化規(guī)律與圖3的曲線相同，二者相互驗證。

結(jié)合圖3、圖4結(jié)果可知，RSM湍流模型精確度最高，與實驗值吻合最好，其次為帶旋流修正的RNGk-ε模型，所得壓力變化趨勢與實驗結(jié)果相一致，關(guān)鍵位置處誤差在可接受范圍內(nèi)，不帶旋流修正的RNGk-ε模型則與實驗結(jié)果相差較遠。從經(jīng)濟性角度出發(fā)，RSM模型相較于帶旋流修正的RNGk-ε模型，增加了約1倍的計算量，且RSM模型的精細化使其難以收斂。因此，綜合考量準確度與經(jīng)濟性因素，最終選擇帶旋流修正的RNGk-ε湍流模型作為計算模型。

2.2 網(wǎng)格依賴性驗證

圖5、圖6分別為不同網(wǎng)格數(shù)計算所得的壓力、切向速度分布曲線。3種網(wǎng)格mesh-1、mesh-2、mesh-3網(wǎng)格數(shù)分別為707 362、451 938、207 816，分別對應(yīng)于致密網(wǎng)格、中等網(wǎng)格及稀疏網(wǎng)格，計算所用湍流模型為2.1節(jié)結(jié)論中所得帶旋流修正的RNGk-ε湍流模型。

從圖5可看到，不同網(wǎng)格所得壓力徑向分布基本重合，變化趨勢以及關(guān)鍵位置處與實驗數(shù)據(jù)吻合得非常好。圖6為不同網(wǎng)格所得切向速度vt徑向分布，3種網(wǎng)格所得vt分布在整個半徑方向基本重合。

綜合圖5、圖6中結(jié)果可知，所選湍流模型及計算設(shè)置對網(wǎng)格的依賴性非常低，即致密網(wǎng)格、中等網(wǎng)格以及稀疏網(wǎng)格均可得較可信的計算結(jié)果。但同樣考慮到計算精確度和經(jīng)濟性因素，致密網(wǎng)格計算較慢，稀疏

網(wǎng)格捕捉流場信息較少。因此，綜合考慮選擇中等網(wǎng)格設(shè)置作為最終計算網(wǎng)格。

3 計算結(jié)果及討論

通常研究者對VCCWC中特別關(guān)心的流場信息包括壓力p分布、軸向速度va分布及切向速度vt分布。從第2章中看到，VCCWC內(nèi)p分布易于常規(guī)固體或液體火箭發(fā)動機，其推力計算有別于常規(guī)火箭發(fā)動機[4,9,11]。va分布決定了渦幔半徑rm的大小，對VCCWC燃燒及冷壁作用均有重要意義。vt分布決定了燃料與氧化劑的摻混程度及壓力梯度力與離心力所形成的力學(xué)平衡。

3.1 入射傾角θ對內(nèi)流場渦結(jié)構(gòu)影響

入射傾角數(shù)值研究中，Δpin=41 700 Pa，αL=1.5。計算所得α處vt徑向分布曲線見圖7。從圖7可見，vt曲線在0<β<0.05范圍內(nèi)接近0 m/s，從β=0.05開始，vt迅速增大，在β=0.19左右達到vt,max，隨后逐漸減小，靠近β=1時，迅速減小為0 m/s，以匹配無滑移壁面條件，該速度分布曲線與Majdalani所得理論曲線變化趨勢相一致[13-17]。隨著θ的增大，vt及vt,max均減小，說明入射傾角θ的增大，減小了來流總動量中切向動量比例。θ=0°與θ=10° 2條曲線基本重合，當(dāng)θ>10°時，vt及vt,max均迅速減小，說明vt減小的速率隨著θ的增大而增大。對于相同θ值，當(dāng)α從0.3增大到0.9時，vt分布曲線幾乎無變化。

圖8為vt,max的無量綱徑向位置βmax沿α變化曲線，看到在整個發(fā)動機內(nèi)，βmax幾乎保持恒定，以上兩方面共同說明，vt軸向衰減率非常小。

圖9為軸向速度va分布曲線，α分別為0.3、0.6、0.9，軸向速度va曲線與0 m/s曲線交點的徑向位置，即為渦幔(Mantle)[1]半徑rm(無量綱渦幔半徑βm定義為rm/R，如圖9中所示)。相同α處，va隨θ增大而增大，這是由于來流總動量中軸向動量所占比例增大導(dǎo)致的，與圖7結(jié)果相一致。當(dāng)α從0.9減小到0.3時，va曲線同樣減小，說明va沿軸向具有顯著衰減變化。

圖10為入射傾角θ下βm的軸向分布，從入射口附近(α=1.2)到發(fā)動機頂端(α=0.2)，βm以線性規(guī)律減小，線性斜率約為0.093 7。βm從α=1.2處的0.81±0.01減小至α=0.2處的0.73±0.01，該變化范圍與文獻[4,6,10]所得結(jié)果相一致。

3.2 入射壓降Δpin對內(nèi)流場渦結(jié)構(gòu)影響

入射壓降Δpin研究中，模型長徑比αL=1.0，即發(fā)動機長度L=100 mm，入射壓降分別為Δpin-1=41 700 Pa和Δpin-2=62 550 Pa，入射傾角θ分別為0°、10°、20°、30°。圖11為α=0.6處無量綱壓力Δp/Δpin徑向分布。可看到，在0.19<β<1范圍內(nèi)，各曲線基本重合，在0<β<0.19范圍內(nèi)，壓力曲線出現(xiàn)差異，Δpin-2壓力曲線在β≈0.19變陡，斜率增大。圖12為θ=0°、10°、20°、30°時，不同Δpin下，α=0.6處vt分布?？煽吹?，從θ=0°增大到30°時，vt均隨Δpin增大而增大，vt,max隨θ增大而減小，這與3.1節(jié)中結(jié)論相一致。比較圖11、圖12可發(fā)現(xiàn)，壓力曲線斜率增大位置與vt,max位置均出現(xiàn)在β≈0.19，這是由渦流理論決定的。對于渦流，存在式(1)所示力學(xué)平衡關(guān)系[18-19]：

(1)

以模型半徑R對式(1)進行無量綱化，并進行變形，可得：

(2)

3.3 長徑比αL對內(nèi)流場渦結(jié)構(gòu)影響

本節(jié)討論中，對表1中所有工況進行了匯總。圖13為Δpin=41 700 Pa所得結(jié)果，圖14為Δpin=62 550 Pa所得結(jié)果。當(dāng)Δpin=41 700 Pa時，各入射傾角下，不同αL所得βm的軸向分布并不相同。αL=1.0所得βm在整個發(fā)動機內(nèi)均大于αL=1.5時，但二者變化趨勢基本相同，均從入射口處附近(α=0.2)至發(fā)動機頂端處附近(對于αL=1.0，發(fā)動機頂端處附近α=0.8；對于αL=1.5，發(fā)動機頂端處附近α=1.2)以線性規(guī)律減小，且線性斜率近似相同。

對于vt,max，從圖13(b)看到，各入射傾角下，不同αL所得βmax在整個發(fā)動機內(nèi)均保持在0.19附近，αL對βmax幾乎無影響。圖13(c)為各入射傾角下，不同αL所得vt,max軸向分布?？煽吹?，αL=1.0所得vt,max在整個發(fā)動機內(nèi)均大于αL=1.5，且當(dāng)αL=1.0時，vt,max在發(fā)動機內(nèi)從入射口處附近至發(fā)動機頂端處附近逐漸減小。當(dāng)αL=1.5時，vt,max在發(fā)動機內(nèi)從入射口處附近至發(fā)動機頂端處附近，呈現(xiàn)出先增大、后減小的變化趨勢。但對于αL=1.0和αL=1.5而言，vt,max變化范圍均在3%以內(nèi)?？芍?，vt,max的軸向衰減率非常小，與3.1節(jié)所得結(jié)論一致。

Δpin=62 550 Pa時，從圖14(a)可知，βm軸向分布與Δpin=41 770 Pa相似，即αL=1.0所得βm在整個發(fā)動機內(nèi)均大于αL=1.5時，且二者線性變化斜率也基本相同，βm變化范圍為0.71±0.01～0.82±0.01。從圖14(b)可知，βmax軸向分布也與Δpin=41 770 Pa時相似。從圖14(c)可知，對于αL=1.0和αL=1.5兩種工況，當(dāng)θ相同時，αL=1.0所得vt,max與αL=1.5所得vt,max在數(shù)值上的差值明顯要大于Δpin=41 770 Pa時，說明隨著Δpin的增大，長徑比αL對于vt,max大小的影響將增大。

4 結(jié)論

(1)切向速度vt及vt,max均隨θ增大而減小，隨Δpin增大而增大，隨αL增大而減小。不同θ、不同Δpin、不同αL下，整個發(fā)動機內(nèi)vt,max的徑向位置βmax均恒定在0.19附近。αL=1.0時，vt,max從發(fā)動機頂端到入射口附近逐漸增大；αL=1.5時，vt,max從發(fā)動機頂端到入射口附近先增大、后減小。vt及vt,max的軸向衰減率非常小，維持在3%以內(nèi)。

(2)軸向速度va隨θ增大而增大，渦幔半徑βm從發(fā)動機頂端開始線性增大到入射口附近，變化范圍為0.71<βm<0.82，不同αL下，該線性斜率近似相同，va軸向衰減率要顯著大于vt及vt,max。

(3)發(fā)動機內(nèi)雙向渦力學(xué)特性符合渦流理論，增大Δpin在增大vt及vt,max的同時，會使壓力曲線在βmax附近梯度值增大。

(4)入射傾角θ相同時，隨入射壓降Δpin的增大，長徑比αL對vt,max的影響將隨之增強。

[1] Vyas A B,Majdalani J,Chiaverini M J.The bidirectional vortex part 1 an exact inviscid solution[C]//39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit 20-23 July 2003,Huntsville,Alabama.

[2] Vyas A B,Majdalani J,Chiaverini M J.The bidirectional vortex part 2 viscous core corrections[C]//39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit,20-23 July 2003,Huntsville,Alabama.

[3] Maicke B A,Majdalani J.Heuristic representation of the swirl velocity in the core of the bidirectional vortex[C]//37th AIAA Fluid Dynamics Conference and Exhibit,25-28 June 2007,Miami,FL.

[4] Chiaverini M J,Malecki M J,Majdalani J,et al.Vortex thrust chamber testing and analysis for O2-H2propulsion applications[C]//39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit,20-23 July 2003,Huntsville,Alabama.

[5] 孫得川,白榮博,劉上.渦流冷壁推力室傳熱模型分析計算[J].計算機仿真,2011,28(4): 87-90.

[6] 孫得川，白榮博，劉上.渦流燃燒發(fā)動機燃燒室數(shù)值模擬[J].彈箭與制導(dǎo)學(xué)報，2011,31(2): 111-116.

[7] Sun De-chuan,Liu Shang.Experimental research on bidirectional vortices in cold wall rocket thruster[J].Aerospace Science and Technology,2011,18(2012): 56-62.

[8] 吳東波，李家文，?？擞?GH2/GO2渦流冷卻推力室設(shè)計與數(shù)值計算[J].火箭推進，2010,36(5): 17-22.

[9] 李家文，唐飛，俞南嘉.推力室渦流冷卻技術(shù)試驗研究[J].推進技術(shù),2012,33(6): 956-960.

[10] 唐飛，李家文，?？擞?渦流冷卻推力室中渦流結(jié)構(gòu)的分析與優(yōu)化[J].推進技術(shù),2010,31(2): 165-169.

[11] 李恭楠，俞南嘉，路強.渦流冷卻推力室流場特征與性能仿真[J].航空動力學(xué)報，2014,29(2): 420-426.

[12] Fang D,Majdalani J,Chiaverini M J.Simulation of the cold-wall swirl driven combustion chamber[C]//39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit,20-23 July 2003,Huntsville,Alabama.

[13] Maicke B A,Majdalani J.A constant shear stress core flow model of the bidirectional vortex[J].Proceedings of Royal Society Association,2009,465(3): 915-935.

[14] Maicke B A,Majdalani J.The compressible cidirectional vortex[C]//44th AIAA/ASME /SAE/ASEE Joint Propulsion Conference & Exhibit,21-23 July 2008,Hartford,CT.

[15] Saad T,Majdalani J.Energy based solutions of the bidirectional vortex[C]//44th AIAA/ASME /SAE/ASEE Joint Propulsion Conference & Exhibit,21-23 July 2008,Hartford,CT.

[16] Maicke B A,Majdalani J.On the compressible bidirectional vortex part 1 a bragg-hawthorne stream function formulation[C]//50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition,09-12 January 2012,Nashville,Tennessee.

[17] Maicke B A,Majdalani J.On the compressible bidirectional vortex part 2 a beltramian flowfield approximation[C]//50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition,09-12 January 2012,Nashville,Tennessee.

[18] Basina I P,Tonkonogii A V,Kukueev B N.Motion of burning particles in cyclone chambers[J].Thermal Engineering,1974,21(3): 99-103.

[19] Georgios H V.Tangential velocity and static pressure distributions in vortex chambers[J].AIAA Journal,1987,25(8): 1139-1140.

(編輯：崔賢彬)

Numerical research of compressible vortex combustion cold-wall chamber

DANG Jin-feng,GAO Ye,LIU Ping-an

(College of Aerospace and Civil Engineering,Harbin Engineering University,Harbin 150001,China)

In order to deeply investigate the inner relationship between the inlet parameters and the structure characteristics of bidirectional vortex in the Vortex Combustion Cold-Wall Chamber(VCCWC),the compressible cold simulation research of the VCCWC was carried out based on three-dimensional physical model.The inlet parameters included injection angle,injection pressure drop and chamber length.Three models,including RSM model,RNGk-εmodel with swirl factor and RNGk-εmodel,were compared with the experimental data in the pre-research.The RNGk-εmodel with swirl factor was finally selected as the numerical turbulent model based on the performance of the accuracy and the cost.Results show that: the tangential velocity and the maximum tangential velocity decrease with the injection angle and length to diameter ration but increase with the injection pressure drop.The dimensionless radial positon of the maximum tangential velocity is kept around 0.19 in the whole chamber regardless of injection angle,injection pressure drop and the length to diameter ratio.For the different length to diameter ratio,the maximum tangential velocity increases for the length to diameter ratio of 1.0 but increases first and then decreases for the length to diameter ratio of 1.5 from the top of the chamber to the injection section.The axial decay rates of the tangential velocity and the maximum tangential velocity are kept within 3%.The dimensionless mantle radius is not kept constant but linearly increases from the top of the chamber to the injection section.The range of the dimensionless mantle radius is from 0.71 corresponding to 0.82.For the same injection angle,the impact of the length to diameter ratio to the magnitude of the maximum tangential velocity will strengthen with the increase of the injection pressure drop.

vortex combustion cold-wall chamber；compressible；turbulent models；injection conditions；length to diameter ratio；simulation

2015-09-06；

2015-12-21。

自然科學(xué)基金面上項目(11372079)；中央高?；究蒲谢?HEUCFD1404)。

黨進鋒(1987—)，男，博士生，研究方向為渦流旋轉(zhuǎn)燃燒發(fā)動機數(shù)值研究。E-mail:dangjinfeng1987@163.com

V434.13

A

1006-2793(2017)01-0016-08

10.7673/j.issn.1006-2793.2017.01.003