李映坤，韓珺禮，2，陳 雄，沈振華

(1．南京理工大學機械工程學院，南京 210094;2．北京機電研究所，北京 100012)

0 引言

準確計算固體火箭發(fā)動機的內流場，對于固體火箭發(fā)動機的結構設計與優(yōu)化有著非常重要的作用。固體火箭發(fā)動機工作過程中的燃氣流動大多呈湍流，推進劑表面的燃燒反應會受到湍流流動的影響，燃燒室內的湍流流動不僅對推進劑的燃燒產生直接影響，而且對噴管內跨聲速流動及噴管外羽流產生影響［1］。

為研究固體火箭發(fā)動機燃燒室內流場，1990年Dunlap R等對直徑為101．6 mm、長度為1 453 mm的燃燒室進行了冷流實驗［2］，利用氮氣流入含有444個圓孔的圓柱表面模擬發(fā)動機推進劑表面燃燒時的側向加質，對諸如徑向速度、脈動速度等重要物理參量進行了測量，為燃燒室內流場的數(shù)值研究［3－6］提供了重要參考，但尚未將 Menter F R提出的 k-ω SST(shearstress-transport)湍流模型［7－8］應用在結構網格固體火箭發(fā)動機燃燒室內流場數(shù)值模擬中的報道。

Menter F R提出的k-ω SST(shear-stress-transport)湍流模型，能充分發(fā)揮k-ε模型對自由流和k-ω模型對壁面受限流動的處理優(yōu)勢，在空氣動力學［8］、航空發(fā)動機［9－10］等領域有著廣泛的應用。

本文基于有限體積法求解Navier-Stokes流動控制方程組，將該湍流模型應用于固體火箭發(fā)動機內流場的數(shù)值模擬，并將計算結果與實驗值及Spalart-Allmaras湍流模型、Wilcox的k-ω兩方程湍流模型進行了對比，為固體火箭發(fā)動機內流場仿真過程中湍流模型的選擇提供了參考。

1 控制方程

1．1 流體控制方程

積分形式的可壓縮非定常Navier-Stokes方程為

式中 ρ為氣體密度;u、v為流體運動速度矢量的2個分量;T為氣體的溫度;E為單位體積氣體的總能量;τ表示應力張量，其具體形式參考文獻［11］。

1．2 湍流模型

為便于對比分析，本文共采用了4種湍流模型。Spalart-Allmaras湍流模型［12］基于經驗和量綱分析，考慮了自由剪切流、高雷諾數(shù)時的近壁區(qū)流動、有限雷諾數(shù)的近壁區(qū)流動等，直接針對渦粘性建立的方程。

Wilcox的k-ω兩方程湍流模型［13］是應用最廣的兩方程渦粘性模式之一，為積分到壁面的不可壓縮或可壓縮湍流兩方程渦粘性模式，該模式不需要顯示的壁面衰減函數(shù)。

Menter提出的 k-ω SST(shear-stress-transport)剪切應力輸運模型(以下簡稱為SST k-ω湍流模型)，該模型通過混合函數(shù)F1將k-ε模型和k-ω模型結合起來，這樣充分發(fā)揮了k-ε模型對自由流和k-ω模型對壁面受限流動的處理優(yōu)勢。具體描述如下:

其中

式中 k為湍動能;ω為比耗散率;μt為湍流粘性系數(shù)，其他參數(shù)的具體形式見參考文獻［7］。

目前，經性傳播已成為我國艾滋病傳播的主要方式，而家庭內配偶間經性傳播已成為艾滋病進一步蔓延的重要因素之一，我國2011年估計的78萬艾滋病患者中經異性傳播占46.5%，其中約1/4為配偶間性傳播［1］。因此，了解配偶間人類免疫缺陷病毒（human immunodeficiency virus，HIV）傳播狀況及其相關影響因素，采取相應措施降低配偶間HIV傳播尤為重要，現(xiàn)將相關研究進展綜述如下。

2003年，為減小原湍流模型對近壁處網格尺寸的依賴性［8］，Menter F R對原始的湍流模型進行了改進(以下簡稱為SST k-ω-2003湍流模型)，其改進后形式如下:

2 數(shù)值求解方法

本文自編程序，采用基于格心的多塊結構網格迎風型有限體積法，求解上述Navier-Stokes流動控制方程組，對流項的離散是計算的關鍵，對于圖1所示的控制體ABCD，流過AB邊通量的計算步驟為

(1)采用具有保單調性的三階MUSCL(Monotone Upstream-centered schemes for Conservation Laws)數(shù)值格式，利用界面左右兩側4個點的信息，計算出控制體界面AB處的物理量(密度、速度、壓力等)，并使用Van abada限制器，以避免間斷處的非物理震蕩，具體形式如下:

其中，S為Van Albada限制器，具體形式為

式中 ε=10－6為一小量，防止上式分母為0。

(2)利用控制體界面處的物理量，采用AUSM-PW通量技術［14］，計算出單位時間內通過控制體邊界的流通量，以i方向上的通量F為例，I+1/2界面上的通量可寫為

式中 c為單元界面聲速;Φ為守恒通量;p為壓力項。

(3)利用積分型的守恒方程，計算出下一個時間步控制體內物理量的平均值。

圖1 控制體示意圖系Fig．1 The diagram of a control volume

粘性項采用Jameson中心差分法離散，時間推進采用三階三步TVD型Runge-Kutta顯式方法，此方法計算過程簡單，內存需求小，但時間步長受穩(wěn)定性條件限制，而必須取得很小。因此，本文采用局部時間步長加速收斂技術來加快定常流場計算的收斂速度。該方法的基本思想是定常問題在求解過程中，時間和空間離散是非耦合的，流場最終的計算結果不受各點的發(fā)展歷程的影響，最終的定常計算結果與所采用的時間步長無關。因此，可適當?shù)馗淖冇嬎愕闹虚g過程，以加快得到最終定常解的目的。

3 結果與分析

為說明SST湍流模型模擬固體火箭發(fā)動機內流場的能力，本文選擇Dunlap R的圓柱冷流實驗模型作為對比，該實驗利用氮氣流入含有444個圓孔的圓柱表面，模擬發(fā)動機推進劑表面典型的入射速度和雷諾數(shù)，Dunlap R對此模型進行了大量的研究，得到了許多有用的實驗數(shù)據，本文將這些實驗數(shù)據與數(shù)值計算結果進行了對比。

該實驗模型包括一個等截面的圓柱和一個收斂擴張的噴管，如圖2所示，圓柱段的長度L=1 453 mm，直徑D=101．6 mm，且L/D=14．3。噴管的收斂擴張段之間有一段半角為20°的圓弧，且該圓弧半徑等于噴管喉部半徑，擴張段下游的形狀由一段拋物線段確定，噴管出口直徑與圓柱段直徑相等。因此，該噴管的擴張比為 2．27。

圖2 冷流實驗模型及邊界條件Fig．2 of the cylindrical-port cold flow model and boundary conditions

根據文獻［2］的實驗，按照文獻［4－6］的處理方法，忽略圓柱壁面上圓孔之間的間隙，將整個圓柱表面作為入口邊界，圓柱表面注入氮氣的速度 uinj=0．001 8 Ma，溫度Tinj=303 K，基于圓柱半徑rw和注入速度的雷諾數(shù)為Reinj=180 00。文獻［4］的計算結果表明，入口處湍動能k和湍動能比擴散率ω的取值對計算結果有很大的影響。因此，本文根據文獻［15］的方法，采用式(24)計算入口處的k和ω:

式中 σv=2．5 ×10－6;lw=1．0。

固體壁面假設為絕熱壁，采用無滑移邊界條件;出口邊界根據馬赫數(shù)判定，當出口為超聲速時，此時所有物理量外推，當出口為亞聲速時，給定環(huán)境反壓101 325 Pa，其他參數(shù)由內向外插值。

3．1 網格影響分析

計算區(qū)域網格劃分的疏密對數(shù)值計算有非常大的影響。因此，有必要開展網格相關性研究，排除網格對計算結果的影響。本文采用多塊對接網格，如圖3所示。

圖3 模擬固體發(fā)動機計算網格Fig．3 The computational grid of model SRM

網格分為粗網格(454×60)、中等網格(539×84)和細網格(654×120)3種，壁面第一層網格的間距為1．0×10－5m。圖4是采用 SST k-ω 湍流模型計算的x/D=10．30處軸向無量綱速度沿著徑向的分布(圖中軸向速度u以中心處速度uc無量綱化，燃燒室半徑rw=50．8 mm)，D為燃燒室直徑，x為距離燃燒室頭部的距離。從圖4可看出，在中心線附近粗網格與實驗值吻合得很好，但在靠近壁面附近，粗網格與實驗值相差很大，而中等網格與細網格整體上都與實驗值吻合得很好，尤其是在靠近中心線處，中等網格計算的結果幾乎與細網格一致。因此，本文后面的計算均采用中等網格。

圖4 x/D=10．30處軸向速度網格影響分析Fig．4 Grid independence study for axial mean velocity at x/D=10．30

3．2 軸線上馬赫數(shù)分布

采用不同的湍流模型計算得到的馬赫數(shù)沿著軸向的分布如圖5所示。

圖5 不同湍流模型計算的馬赫數(shù)沿著軸線的分布Fig．5 Comparison of computed and experimental axial mach number with different turbulent models

由圖5可見，馬赫數(shù)沿著軸線呈線性增加趨勢，在x/D小于7．5時，SST k-ω-2003湍流模型與實驗值吻合得非常好;其次是 SST k-ω湍流模型，Wilcox k-ω和Spalart-Allmaras湍流模型與實驗值相差較大;當x/D超過7．5時，4種湍流模型的計算結果均與實驗值有一定差別。這是因為實驗中氣體的流動狀態(tài)迅速轉變?yōu)橥牧?，而?shù)值模擬中這一過程卻相對較慢，但SST kω-2003湍流模型與實驗值的誤差是最小的，Spalart-Allmaras湍流模型次之。

3．3 不同位置處軸向速度對比

圖6是采用不同的湍流模型計算的軸向速度與實驗值的對比，沿著軸線圖中，共給出了10個位置處的軸向速度沿著徑向分布(圖中軸向速度u以中心處速度uc無量綱化)。整體上來看，本文計算的結果與實驗值的趨勢是一致的，但在x/D=0．62處，4種湍流模型計算的結果都小于實驗值，0．5≤r/rw≤1．0范圍內，Spalart-Allmaras湍流模型計算的結果最接近實驗值，但仍有很大的差距，靠近中心線附近SST k-ω-2003湍流模型計算的結果幾乎與實驗值一致。造成這一差異的原因是x/D=0．62處燃燒室前端軸向速度較小。此時，燃燒室表面注入的徑向氣流正在向軸向方向轉變，如圖7所示，燃燒室前端x/D=0．62處，最大軸向速度僅為3．07 m/s，而燃燒室尾部x/D=12．72處的最大軸向速度是前端的17倍。同時，文獻［2］的實驗結果表明，此處形成了比較強的渦流，并沿著軸線，渦流向燃燒室中心附近匯聚，湍流的效果越來越明顯。所以，這就解釋了x/D=1．80處的計算結果相比于x/D=0．62處更加接近實驗值的原因。另外，從圖7還可看出，靠近壁面附近SST k-ω-2003湍流模型計算的結果明顯優(yōu)于其他幾個湍流模型，Wilcox k-ω湍流模型的結果與實驗值相差甚遠;在中心軸線附近，靠近燃燒室尾部的位置(r/rw≥7．88)，SST k-ω湍流模型的結果與實驗值吻合得很好，而SST k-ω-2003湍流模型計算的結果相比實驗值略微偏小，在x/D=7．88位置處達，到最大誤差，約為5．1%。產生這一差異的原因為燃燒表面注入的徑向氣流往下游移動，流動狀態(tài)從層流轉變?yōu)橥牧?，相比徑向與周向運動速度分量，軸向速度分量占據主導地位;同時，文獻［2］中的實驗分析表明，靠近燃燒室尾部處，渦量被限制在中心線附近。因此，基于渦量的SST k-ω湍流模型的計算值與實驗值吻合得很好，而基于剪切率的SST k-ω-2003湍流模型的計算值與實驗值略微有點差距，但總體趨勢還是一致的。

3．4 不同位置處湍流強度對比

基于Boussinesq渦粘性假設的兩方程湍流模型很難準確預測燃燒室內的湍流強度分布［6］，而大渦模擬方法和直接數(shù)值模擬方法預測精度較高，但其計算量較大，以目前的計算機資源很難在工程上廣泛應用。本文計算的不同位置處的湍流強度I沿著徑向的分布如圖8所示。由圖8可見，幾個湍流模型計算的結果與實驗值的趨勢是一致的，但都有一定的差異，特別是Wilcox k-ω湍流模型與實驗值相差很大，而SST湍流模型計算的結果更接近于實驗值，尤其是 SST k-ω-2003湍流模型在x/D=10．30位置處，幾乎與實驗值是重合的，但在加質壁面附近，SST k-ω-2003湍流模型計算的結果與實驗值還略微有些差距。文獻［3］認為，這與加質壁面湍流邊界條件的取值有關，而在中心線附近，SST k-ω-2003湍流模型的計算結果與實驗值吻合得很好。

圖6 不同位置處軸向速度與實驗值對比Fig．6 Comparison of computed and experimental axial velocity profiles

圖7 燃燒室前、后端速度沿徑向的分布Fig．7 Velocity distribution at the head end and the aft of combustion chamber

由圖8進一步分析可知，靠近燃燒室前端湍流強度很小，隨著氣流向下游移動，湍流強度逐漸增大，本文SST k-ω-2003湍流模型也成功地模擬出了這一趨勢，相比于 x/D=3．04位置，x/D=5．46處湍流強度明顯增大，在x/D=6．64處增大得更加明顯，這說明在此位置，流動正在由層流轉變?yōu)橥牧鳌?/p>

4 結論

(1)采用不同湍流模型計算得到的馬赫數(shù)沿軸線均呈線性增加趨勢，當x/D小于7．5時，SST k-ω-2003湍流模型與實驗值吻合得非常好;當x/D超過7．5時，4種湍流模型的計算結果均與實驗值有一定差別，但SST k-ω-2003湍流模型與實驗值的誤差最小。

(2)4種湍流模型都能準確模擬出固體火箭發(fā)動機燃燒室內的徑向速度分布，且計算結果差別不大，SST k-ω-2003湍流模型的計算結果與實驗值吻合得最好，最大誤差約為5．1%。

(3)SST k-ω-2003湍流模型計算的固體火箭發(fā)動機燃燒室內湍流強度分布與實驗的規(guī)律一致，而其余湍流模型計算的結果與實驗值有很大差異。因此，與其他湍流相比，SST k-ω-2003湍流模型能較準確地模擬固體火箭發(fā)動機的內流場，具有應用于固體火箭發(fā)動機內流場仿真的可行性和相對優(yōu)勢。

圖8 不同位置處湍流強度與實驗值對比Fig．8 Turbulence intensity in vertical direction at various axial locations

［1］ 武曉松，陳軍，王棟，等．固體火箭發(fā)動機工作過程數(shù)值仿真［M］．北京:高等教育出版社，2005．

［2］ Dunlap R，Blackner A M，Waugh R C，et al．Internal flow field studies in a simulated cylindrical port rocket chamber［J］．Journal of Propulsion and Power，1990，6(6):690-704．

［3］ Sabnis J S，Gibeling H J，McDonald H．Navier-Stokes analysis of solid propellant rocket motor internal flows［J］．Journal of Propulsion and Power，1989，5(6):657-664．

［4］ Arabshahi A，Webster R S，Briley W R，et al．Numerical analysis of solid propellant rocket motor internal flows［R］．AIAA 2006-5114．

［5］ Arabshahi A，Sreenivas K，Nichols D S，et al．Computational analysis of turbulent internal flow in ballistic solid rocket motors［R］．AIAA 2007-1449．

［6］ Yumusak Mine．Analysis and design optimization of solid rocket motors in viscous flows［J］．Computers ＆ Fluids，2013，75:22-34．

［7］ Menter F R．Two-equation eddy-viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32(8):1598-1605．

［8］ Menter F R，Kuntz M，Langtry R．Ten years of industrial experience with the SST turbulence model［J］．Turbulence，Heat and Mass Transfer，2003(4):625-632．

［9］ 姚玉，張靖周，郭文．氣膜孔角度對導葉冷卻效果影響的數(shù)值研究［J］．航空動力學報，2009，24(3):507-512．

［10］ 曾軍，王彬，卿雄杰．某雙級高壓渦輪全三維計算［J］．航空動力學報，2012，27(11):2553-2561．

［11］ 張涵信，陳堅強，高樹椿．H2/O2燃燒的超聲速非平衡流動的數(shù)值模擬［J］．宇航學報，1994，15(02):14-23．

［12］ Spalart P R，Allmaras S R．A one-equation turbulence model for aerodynamic flows［R］．AIAA 92-1439．

［13］ Wilcox D C．Turbulence modeling for CFD［M］．La Canada:DCW Industries，1998．

［14］ Meng-sing Liou．Methods for the Accurate Computations of Hypersonic Flows［J］．Journal of Computational Physics，2001，174(1):38-80．

［15］ Sachdev J S．Parallel solution-adaptive method for predicting solid propellant rocket motor core flows［D］．Doctor of Philosophy Graduate Department of Aerospace Science and Engineering University of Toronto，2007．