雙向渦數值模擬雷諾應力各向異性研究與分析

黨進鋒，郜 冶，劉加寧

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

?

雙向渦數值模擬雷諾應力各向異性研究與分析

黨進鋒，郜 冶，劉加寧

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

以實驗數據和理論結果作為參考，探討基于布茨涅克各向同性假設的RNGk-ε模型、SST模型和各向異性RSM湍流模型在雙向渦流數值計算中的適用性，并對不同入射速度下的雙向渦流結構進行了比較分析。計算結果表明，基于布茨涅克各向同性假設的RNG模型與SST模型在不同入射速度下均無法得到內部受迫渦、外部自由渦的雙向渦流結構，而RSM模型則可以準確得到雙向渦流結構。內、外渦分界面處雷諾正應力各向異性特征顯著，且隨入射速度的提高而增強，最大差值達到了198%，軸線附近流場雷諾正應力各向異性特征相對較弱，但是最大差值依然達到46%，因此基于布茨涅克各向同性假設的湍流模型不適用于渦流冷壁發動機數值模擬研究。

雙向渦流；渦流冷壁發動機；雷諾應力；各向異性

0 引言

渦流是自然界中普遍存在的一種流動形式，大到天體星云的運動，小到湍流中的粘性耗散渦，均是不同形式的渦流。由于渦流具有豐富多樣的結構形式，近年來被越來越多的應用在流體機械中，如渦流燃燒器、渦流固液組合發動機、渦流分離器等。2002年，美國Orbital科技公司Chiaverini等提出了一種新概念液體火箭發動機形式[1]，氧化劑從發動機噴管上方沿切向噴入發動機內形成外部渦流，外部渦流向發動機頂端旋轉流動，在發動機頂端改變流動方向，并與頂端噴入的燃料混合形成內部燃燒渦流，燃燒的內部渦流從噴管流出產生推力，該發動機被稱為渦流冷壁發動機(Vortex Combustion Cold-Wall Chamber VCCWC)。VCCWC的優點在于氧化劑所形成的外部渦流對于發動機壁面起到了冷卻作用，大大降低了發動機內壁熱防護層厚度，核心流內氧化劑與燃料旋轉燃燒，增加了反應時間從而縮短了反應距離，最終可減小發動機長度，該兩方面優點均有助于降低發動機成本，因此美國NASA將VCCWC確定為未來廉價發動機的候選方案[2-3]。由于VCCWC內部渦與外部渦流動方向相反，因此該組合渦稱為雙向渦流(Bidirectional Vortex)。

從VCCWC中雙向渦流問題的提出至今十余年時間里，有限研究工作主要分為兩部分，一部分為Chiaverini等所做VCCWC冷流與燃燒實驗[4]，實驗驗證了雙向渦流的存在并得到內渦旋轉燃燒圖像。另一部分為Majdalani等開展的雙向渦流理論研究，文獻[5]中推導得到不可壓無粘流雙向渦流結構解析式，文獻[6]中推導得到不可壓粘性核心雙向渦流結構解析式，文獻[7]得到定常剪切應力核心流模型，文獻[8-9]中分別通過Bragg-Hawthorne流函數法和Beltramian流場近似法得到可壓縮雙向渦流模型。

與之相比，雙向渦流CFD研究進展較慢。由于雙向渦流三維效應顯著，且涉及到VCCWC化學反應數值計算，計算量巨大，因此LES、DNS等方法的數值研究還未開展。Fang在文獻[10]中結合Standardk-ε模型和RSM模型數值計算模擬了VCCWC冷流與燃燒過程，從其冷流計算結果來看速度曲線與理論曲線趨勢相似，能夠得到內、外雙向渦流結構。文獻[11]中通過數值計算指出外部渦流厚度約占發動機半徑的25%～29%，并提出冷卻劑噴嘴傾斜一定的角度可增強內部渦流強度、改善推進劑摻混程度提高燃燒效率。文獻[12]中數值計算研究了以GH2/GO2為燃料的VCCWC燃燒過程。

文獻[13]中通過二維物理模型比較了RNG 模型和RSM模型在VCCWC數值模擬中的準確性，指出RSM模型更適合于強旋流場。本文在文獻[13]工作基礎上，對VCCWC內流場進行了三維完整物理模型、3種不同湍流模型、3種不同入射速度下的數值計算研究，定量探究RSM模型更適合VCCWC數值計算的內在原因。計算發現雙向渦流結構對于湍流模型的選擇性非常大，不合適的湍流模型計算所得結果與實驗結果差別非常明顯，甚至可以說是完全錯誤的，分析原因認為：雙向渦流與普通流動的區別在于，雙向渦流在遠離壁面的核心流區域內湍流脈動雷諾正應力各向異性特性依然異常顯著，而基于布茨涅克各向同性假設(Boussinesq Isotropic Hypothesis)的湍流模型則認為核心流區域各向雷諾正應力均相等[14]，這就導致了傳統基于布茨涅克各向同性假設的湍流模型無法得到正確結果。

本文以文獻[15]中實驗數據及文獻[5，16]中理論結果作為參考，探討基于布茨涅克各向同性假設的RNGk-ε模型、SST模型和各向異性RSM湍流模型在雙向渦流數值計算中的適用性，并對不同入射速度下的雙向渦流結構進行了比較分析，結論可用于指導雙向渦流數值模擬湍流模型的選擇。

1 物理模型及數值計算方法

圖1(a)、(b)分別為雙向渦流物理模型示意圖及入射噴口分配示意圖。物理模型直徑d=100 mm，長度L=100 mm，即長徑比α=1，模型頂部為平面。入射噴口橫截面為5 mm×4 mm矩形，模型底部12個入射噴口周向均勻排列，采用切向速度入口形式。本文中物理模型與文獻[10，13]中類似，因此采用與文獻[10，13]相同的Outflow出口邊界條件，出口直徑為30 mm，其余均為無滑移壁面。

(a)VCCWC物理模型 (b)入射噴口分配示意圖

物理模型網格示意圖如圖2所示。計算中，3種湍流模型(RNGk-ε模型、SST模型及RSM模型)均采用SIMPLE算法求解，梯度差分采用G-S節點基格式，壓力差分采用二階格式，其余通量差分均為二階迎風格式。Majdalani在文獻[5]中指出雙向渦流在不可壓縮與可壓縮流動中均存在，因此計算中采用低速不可壓縮流形式，流體介質為理想不可壓空氣，3種入口速度分別為10、15、20 m/s。計算工作在64位服務器上完成，耗時320 h，所使用的CFD軟件為ANSYS FLUENT11.0。

圖2 物理模型網格示意圖Fig.2 Geometric mesh of physical model

2 計算可信性驗證

圖3為入射速度vin=20 m/s時L=50 mm處無量綱壓力差徑向分布曲線，3條曲線分別為RNG、SST以及RSM模型計算所得，并與文獻[15]中實驗數據進行了比較，文獻[15]中只給出了0.00

圖4所示為本文中RSM模型計算所得速度矢量與Majdalani在文獻[5]中所得理論速度矢量比較。

圖3 L=50 mm無量綱壓力徑向分布與文獻[15]實驗數據Fig.3 Non-dimensional pressure difference distribution inradial direction at L=50 mm and experimentalresults in Reference[15]

可看到在不考慮出口對流場的影響時，CFD結果與理論結果非常一致，Majdalani在文獻[5]推導過程中忽略了出口直徑對于流場的影響，假設出口直徑與內渦半徑相等，因此內、外渦流動方向相反，不存在回流區。當出口直徑小于內渦半徑時就會存在如本文中數值計算中所示的回流區。

綜合以上壓力與速度矢量圖比較可知，本文中RSM計算結果與實驗數據及理論結果非常符合，說明本文中所采用物理模型、計算網格、邊界條件及數值求解方法具有足夠的可信度。

3 計算結果及討論

3.1 不同湍流模型速度及雙向渦流線分析

圖5為不同入射速度下L=50 mm處切向速度vt徑向分布曲線(圖5(a)～(c)所示)。

vin分別為10、15、20 m/s，以及文獻[16]雙向渦vt理論曲線圖如5(d)所示，理論曲線使用模型半徑和入口速度對橫、縱軸進行了無量綱化。RSM模型所得速度曲線：vin=10 m/s時vt，max位于r≈±0.045 m處，vin=15 m/s時vt，max向軸線位置處移動位于r≈±0.017 m處，vin=20 m/s時vt，max位于r≈±0.015 m處。同時vt，max隨vin的提高而增大，說明內部受迫渦強度隨vin提高而增大。RNG模型和SST模型在不同入射速度情況下vt，max始終位于r≈±0.045 m處。

圖6所示為vin=20 m/s時RNGk-ε模型、SST模型和RSM模型所得雙向渦流結構流線圖，在物理模型、入射速度及其他邊界條件均相同的情況下，3種湍流模型所得流線圖差別非常顯著。RNG模型和SST模型也可得到雙向旋轉流動結構，結合圖3、圖5曲線可知此時并沒有形成內部受迫渦，整個流場均為自由渦，所以雙向流動并沒有延伸至整個內部流場，旋流長度僅為流場長度的2/3，外部渦流半徑也僅為流場半徑的2/3。RSM模型所得流線則可以形成完整的雙向渦流結構。

從引言中可知，基于雙向渦流的VCCWC優點在于外部渦流能夠對于發動機起到冷壁作用，所以外部渦流半徑應接近于發動機半徑，內部渦流高速旋轉能夠縮短化學反應距離。因此，vt，max應位于內部受迫渦范圍內，這均與RSM計算結果相一致。

綜合圖3中壓力曲線、圖5中速度曲線及圖6流線圖可知，SST模型和RNG模型不能得到正確的雙向渦流場分布，RSM模型能夠準確計算得到雙向渦流結構。

(b)vin=15 m/s

(c)vin=20 m/s

(d)文獻[16]理論結果

3.2 RSM模型雙向渦結構分析

圖7為不同入射速度RSM模型計算所得L=50 mm處軸向速度va徑向分布曲線。雙向渦流內、外渦軸向速度方向相反，因此va曲線與0.00 m/s速度曲線相交點為內、外渦分界面位置，軸線到該分界面距離即為渦幔半徑rm(Mantle Radius)。vin=10m/s時，rm≈0.042 m，vin=15m/s時，rm≈0.035 m，vin=20 m/s時rm與vin=15 m/s時基本相同。

(a)RNG模型

(b)SST模型

(c)RSM模型

圖8所示為vin=20 m/s時RSM模型所得不同軸向位置處va分布曲線，軸向位置分別為L=25 mm、L=50 mm和L=75 mm。從圖8看到，渦幔半徑從頂端到出口處略有變化，但是基本保持在0.035 m左右，這與現有文獻結論非常吻合[5，11，17]。

圖7 RSM 模型計算所得L=50mm處軸向速度徑向分布Fig.7 Axial velocity distribution in radial direction atL=50mm using RSM model

圖8 vin=20 m/s時RSM模型計算所得不同位置軸向速度徑向分布Fig.8 Axial velocity distribution in radial direction in inlet velocity of vin=20 m/s using RSM model

3.3 雷諾應力各向異性分析

從圖3～圖6看到，基于布茨涅克各向同性假設的RNG模型與SST模型均無法得到正確的雙向渦流結構，而RSM模型則能夠準確得到完整的雙向渦流結構，壓力差分布曲線(圖3)與實驗數值吻合的非常好，對稱面速度矢量分布(圖4)與理論結果符合較好。從國內外學者關于雙向渦流數值計算文章中可看到，部分文章依然采用基于布茨涅克各向同性假設的湍流模型進行數值計算研究。為了證明基于布茨涅克各向同性假設湍流模型的不足，圖9～圖11將RSM模型計算所得不同入射速度情況下各軸向位置處雷諾正應力分布曲線進行比較分析。

(a)vin=10 m/s (b)vin=15 m/s (c)vin=20 m/s

(a)vin=10 m/s (b)vin=15 m/s (c)vin=20 m/s

(a)vin=10 m/s (b)vin=15 m/s (c)vin=20 m/s

比較3.1節與本節雷諾正應力分析可知，RNG模型與SST模型無法得到正確雙向渦流的原因在于雙向渦流結構中雷諾正應力各向異性特性顯著，內、外渦分界面附近雷諾正應力最大差值達到198%，且隨內部受迫渦強度增大還會進一步提高，此時布茨涅克各向同性假設已經無法滿足雙向渦流物理特性，因此RNG模型與SST模型無法得到正確的雙向渦流，推論可知基于布茨涅克各向同性假設的其他湍流模型，如Standardk-ε模型，Standardk-ω模型等也不適用于雙向渦流數值計算。雷諾應力各向異性決定了湍流脈動的各向異性，化學反應數值計算準確性很大程度上決定于湍流場計算的準確性，因此未來無論進行VCCWC冷流或燃燒數值模擬時均應摒棄基于布茨涅克各向同性假設的湍流模型。

4 結論

(1)基于布茨涅克各向同性假設的RNG模型與SST模型在入射速度從vin=10 m/s增大到vin=20 m/s時均無法準確得到內部受迫渦、外部自由渦的雙向渦流結構，RSM模型能夠準確得到與實驗與理論結果相符的雙向渦流結構。

(2)對于本文中所使用的物理模型尺寸，當入射速度vin=10 m/s時無法形成穩定的內部受迫渦，當vin≥15 m/s時形成穩定的雙向渦結構，在一定的速度變化范圍內渦幔半徑保持在rm=0.035 m左右。隨入射速度提高內部受迫渦強度增大，因此最大切向速度vt，max增大，其位置隨入射速度提高而靠近軸線位置。

(3)當雙向渦流形成時，內、外渦分界面附近雷諾正應力各向異性特征顯著，且隨入射速度的提高而增強，最大差值達到了198%，軸線附近流場雷諾正應力各向異性特征相對較弱，但是最大差值依然達到46%，因此基于布茨涅克各向同性假設的湍流模型不適用于VCCWC雙向渦流數值模擬研究。

[1] Chiaverini,M J,Malecki,M J,Sauer,J A,et al.Vortex combustion chamber development for future liquid rocket engine applications[R].AIAA2002-4149.

[2] Turner A E.Aquarius low cost launch main engine study[C]//3rd Responsive Space Conference April 25-28,2005 Los Angeles,CA.

[3] Munson S M,Sauer J A,Rocholl J D,et al.Development of a low-cost vortex-cooled thrust chamber using hybrid fabrication techniques[C]//47th AIAA/ASME/SA.E/ASEE Joint Propulsion Conference & Exhibit,31 July-03 August 2011,San Diego,California.

[4] Chiaverini M J,Malecki M J,Sauer A,et al.Vortex thrust chamber testing and analysis for O2-H2propulsion applications[R].AIAA 2003-4473.

[5] Vyas A B,Majdalani J.The bidirectional vortex.part 1:an exact inviscid solution[R].AIAA2003-5052.

[6] Vyas A B,Majdalani J.The bidirectional vortex.Part 2:viscous core corrections[R].AIAA2003-5053.

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

[8] 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.

[9] 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.

[10] Fang D,Majdalani J.Simulation of the cold-wall swirl driven combustion chamber[R].AIAA 2003-5055.

[11] 唐飛,李家文,?？擞?渦流冷卻推力室中渦流結構的分析與優化[J].推進技術,2010,31(2):165-169.

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

[13] 孫得川,白榮博,劉上.渦流燃燒發動機燃燒室數值模擬[J].彈箭與制導學報,2011,31(2):111-114.

[14] Schmitt F G.About boussinesq's turbulent viscosity hypothesis:historical remarks and a direct evaluation of its validity[J].Comptes,Rendus,Mecanique,2007,9-10,335:617-627.

[15] Maicke B A,Majdalani J.Characterization of the Bidirectional vortex using particle image velocimetry[D].PhD.Giovanna Cavazzini (Ed.),ISBN:978-953-51-0625-8,InTech.

[16] Majdalani J,Halpenny E K.The bidirectional vortex with sidewall injection[C]//44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,21-23 July 2008,Hartford,CT.

[17] Smith,J L.An analysis of the vortex flow in the cyclone separator[J].Journal of Basic Engineering Transactions of the ASME,1962:609-618.

(編輯：呂耀輝)

Research and analysis of reynolds stress anisotropic in bidirectional vortex simulation

DANG Jin-feng，GAO Ye，LIU Jia-ning

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

In this paper,simulation results were compared with experimental data and theoretical results.The simulation validity of RNGk-εmodel and SST model based on Boussinesq isotropic hypothesis and RSM model in bidirectional vortex were discussed and bidirectional vortex properties in different inlet velocities were analyzed as well.It is pointed out that the RNG model and SST model could not achieve the bidirectional vortex which frees in outer area and is forced in inner space under any inlet velocity based on Boussinesq isotropic hypothesis,while the RSM model simulation results fit the experiment data and theoretical solutions well.In mantle surface of inner vortex and outer vortex,Reynolds stress anisotropy characteristics is significant and is strengthened with the increase of inlet velocity.The biggest difference between normal stresses inx,y,zdirection can reach up to 198%.In core flow area close to axial,Reynolds stress anisotropy is relatively weak,but the greatest difference can still reach as large as 46%.Therefore,it is believed that turbulent models based on Boussinesq isotropic hypothesis are invalid in VCCWC simulation.

bidirectional vortex；vortex combustion cold-wall chamber；reynolds stress；anisotropic

2014-10-30；

：2014-12-11。

自然科學基金面上項目(11372079)。

黨進鋒(1987—)，男，博士生，研究方向為發動機內流場數值仿真。E-mail：xiandjf@163.com

V430

A

1006-2793(2015)05-0640-06

10.7673/j.issn.1006-2793.2015.05.007