基于動網格的喉栓式推力可調噴管內流場數值模擬①

唐金蘭，宋慧敏，李進賢，馮喜平

(西北工業大學燃燒、熱結構與內流場重點實驗室，西安 710072)

0 引言

喉栓式推力可調噴管能夠實現固體火箭發動機的推力連續可調，從而提升了導彈武器的機動性和突防能力，具有強烈的軍事需求背景。因此，國外對喉栓式推力可調噴管發動機的相關技術在理論方面、數值模擬、試驗研究以及相關的基礎學科研究領域都開展了較多的研究［1］。國內對于喉栓式推力可調噴管的相關研究［2］起步較晚，且研究多集中在性能分析與流場的穩態數值模擬方面，基于動網格的喉栓式推力可調噴管內流場的動態特性的研究鮮見報道。

本文利用動網格技術，對喉栓調節運動過程中噴管的內流場進行了數值模擬和分析，分析結果對喉栓式推力可調噴管的實驗器設計及實驗研究具有一定的指導意義。

1 數理模型與計算方法

1．1 物理模型

由文獻［3］的研究結果可知，燃氣進入噴管的角度對噴管內流場的影響不大。因此，本文對軸向進氣的喉栓式推力可調噴管進行軸對稱二維內流場分析，以探討喉栓調節運動過程對噴管內流場特性的影響。喉栓式可調噴管的物理模型如圖1所示。

圖1 可調噴管示意圖Fig．1 Diagram of pintle-controlled nozzle

式(1)給出了喉栓式推力可調的基本原理。

式中 ρp、C*、a、n分別為固體推進劑的密度、特征速度、燃速系數和壓強指數;Ab為推進劑藥柱燃面面積;At為發動機噴管喉部面積。

對于選定的推進劑，ρp、C*、a、n可視為常量。假定推進劑藥柱等面燃燒，即Ab不變，則At是唯一的變量，只要改變At就可有效地改變燃燒室的壓強，從而改變了發動機的推力(F=CFpcAt)。

1．2 數學模型

1．2．1 控制方程

采用N-S方程，將雷諾守恒定律得到的質量、動量和能量守恒方程寫成微分形式，就得到了考慮氣體粘性和熱傳導等流體屬性的、描述燃氣流動的非線性偏微分方程組，式(2)為N-S方程的張量形式:

式中 φ 為通用變量(φ=1，u，v，w，T分別對應連續方程、動量方程和能量方程);Γφ是與φ對應的廣義擴散系數;湍流脈動附加項(Γi是與 φ 對應的湍流擴散系數);S為廣義源項。

1．2．2 湍流模型

采用RNG k-ε湍流模型，湍流動能k方程和湍流耗散率 ε 方程［5］如下:

湍動能k方程:

湍流耗散率ε方程:

式中 Gk是由平均速度梯度引起的湍流動能;Gb是由浮力影響引起的湍流動能;YM可壓縮湍流脈動膨脹對總的耗散率的響應;C1ε、C2ε、C3ε是常量(一般取 C1ε=1．44，C2ε=1．92，C3ε=0．09);αk、αε是 k、ε 方程的湍流普朗特數(一般取 αk=1．0，αε=1．3)。

1．2．3 動網格技術

采用動網格技術后，計算區域是變化的。所以，在任意的控制體對任意廣義標量φ，有積分守恒方程［6］:

式中 Vs是控制體體積;Ls是控制體的表面邊界;u是流體時均速度;ug是動網格邊界移動速度;n是表面邊界的外法線單位向量;Γ是擴散系數;qφ是源項。

對于喉栓調節運動區域的網格，采用動態層更新方法，即根據相鄰運動邊界網格層的高度變化，合并或者分割動態網格層，如圖2所示。

圖2 動態網格層的合并與分割Fig．2 Diagram of dynamic layering update by merging and splitting

如果分割網格層，網格單元高度的變化有一臨界值:

式中 hmin為網格單元的最小高度;h0為理想網格單元高度;αs為網格層的分割因子。

在滿足式(5)的情況下，采用常值高度法(即網格單元分割的結果是產生相同高度的網格)進行網格分割。本文中，喉栓表面是運動邊界，而其他部分保持靜止。因此，在喉栓表面的運動邊界上運用動網格技術。

1．3 初始與邊界條件

(1)初始條件。初始壓強、初始質量流率、總溫分別為2 MPa、0．077 34 kg/s、1 500 K 的條件下的定常計算結果作為本文的初始條件。

(2)入口邊界條件。定義為自適應的質量流率入口。通過UDF程序讀入燃燒室壓力，自動調整入口質量流量，使其滿足質量流量變化規律。喉栓從完全打開位置向噴管幾何喉部運動時，噴管質量流率的變化由式(6)給出。

式中 n為固體推進劑壓強指數;pc、At、m·為喉栓未介入時的燃燒室壓強、噴管幾何喉部面積和初始質量流率;pc，x、At，x、m·x為喉栓相應調節位置下的燃燒室壓強、噴管等效喉部面積和質量流率。

由式(6)可見，只要計算出喉栓運動到某位置時的噴管等效喉部面積，即可得出該位置下相應的質量流率和燃燒室壓強。

(3)出口邊界條件。定義出口邊界壓強為101 325 Pa。由于噴管內流場出口為超聲速流動，壓強出口邊界各參數外推得到。

(4)壁面邊界。采用無滑移邊界條件，近壁區采用標準壁面函數。

(5)對稱邊界。對稱軸上法向速度為零，其他所有流動變量的梯度為零。

1．4 計算分析方法

1．4．1 網格劃分

數理模型建立后，計算分析前，應對物理模型進行網格劃分［4］。網格劃分如圖3所示，區域1和區域3采用四邊形結構化網格，區域2采用三角形非結構化網格，并進行網格細化，以防止出現由于喉栓的調節運動而可能出現的負網格。

圖3 喉栓式推力可調噴管內流場網格劃分Fig．3 Meshing of pintle-controlled nozzle flow field

1．4．2 計算分析流程

計算分析流程主要包括流場壓強、質量流率等流場參數的更新、流場計算分析、網格更新等幾個過程，如圖4所示。計算時，需編寫喉栓頭部運動函數，用函數DEFINE_CG_MOTION定義喉栓頭部運動速度。

圖4 計算分析流程圖Fig．4 Flow chart of CFD program

2 計算結果分析與討論

2．1 計算結果與分析

圖5 給出了 0．1、0．5、1．0 m/s調節速度下，喉栓從未介入到完全介入過程中，噴管入口靜壓與時間的變化。從圖5可見，調節速度越大，調節喉栓所需的時間越短、壓強的變化趨勢也越快，但在調節的初始位置，調節速度越大，噴管的入口靜壓就越小(小于初始條件給定的壓強值2 MPa)，導致這一現象的原因主要是喉栓的快速運動，使壓強建立過程出現了一定的延遲。

圖5 入口靜壓隨調節時間的變化關系Fig．5 Inlet static pressure with different regulation time

對推力可調噴管喉栓處于不同調節速度下的內流場進行了動態特性分析，調節速度為1．0 m/s時，3個位置的分析結果如圖6所示。

圖6 調節喉栓處于不同位置時噴管內的壓強和馬赫數分布Fig．6 Pressure and Mach number distribution of pintle-controlled at different positions

由圖6(a)可見，由于喉栓尚未介入噴管喉部，喉栓對噴管喉部處的內流場幾乎沒有影響。此時，噴管喉部的燃氣膨脹均勻、達到聲速，噴管等效喉部面積最大。為50．27 mm2，噴管入口靜壓最低，為1．8 MPa，噴管出口處的平均馬赫數為1．34。

由圖6(b)可見，當喉栓處于調節介入的中間階段時，除等效喉部處外，在喉栓頭部尖端附近也形成了一個馬赫數為1的小區域。分析該區域產生的原因，是由于在喉栓頭部尖端出現了氣流分離現象，使氣流在喉栓頭部尖端區域的壓強(約1．61～2 MPa)明顯比周圍區域的壓強(0．79～1．34 MPa)要高，而馬赫數則比周圍區域略低。由于此時喉栓的調節介入，噴管等效喉部面積減小為27．47 mm2、噴管入口靜壓上升達到5．35 MPa。因此，噴管出口處的平均馬赫數也隨之提高到約 2．31。

由圖6(c)可見，當喉栓完全介入噴管幾何喉部位置，氣流明顯受到喉栓的影響，而在喉栓頭部尖端變化劇烈，并在擴張段中心區域形成了一個受喉栓影響產生的低馬赫數區域。喉栓處于完全介入狀態。此時，噴管等效喉部面積達到最小為21．99 mm2，噴管入口靜壓也達到最大為11．29 MPa。因此，噴管出口處的平均馬赫數也隨之達到最大值2．9。

由圖6(b)、(c)可見，隨著喉栓調節介入的不斷深入，喉栓頭部對流場的擾動越來越大，從而使喉栓頭部處氣流的壓強、馬赫數與周圍流場的差異也越來越大，在完全介入情況下，喉栓頭部處的流動分離現象也最明顯。同時，由于喉栓調節對流場的擾動，在喉栓頭部存在馬赫數為1的等值線區域，在該區域產生了斜激波，該激波在噴管擴張段下游發生反射，并隨著喉栓的調節深入，該激波及其反射區也向噴管出口移動，燃氣經過該斜激波后的壓強略有增大。

2．2 喉栓調節速度對流場特性的影響

不同調節速度下，典型的計算結果如表1所示。由表1可知，在不同的喉栓介入狀態下，隨著喉栓調節速度的增加，噴管入口靜壓的建立有一定的延遲，需要經過緩沖一段時間才能達到，調節速度越大，壓強建立的延遲情況越明顯。分析造成這種現象的原因，是由于喉栓調節速度越大，喉栓對流場的擾動越大，壓強建立的延遲程度也越大，從而導致的發動機推力的建立延遲情況越明顯。

2．3 穩態和非穩態流場特性分析

為便于比較，喉栓調節速度為1．0 m/s、處于不同調節位置時噴管內流場的穩態壓強分布如圖7所示，典型的數據結果如表2所示。

表1 不同喉栓調節速度仿真結果的對比Table 1 Comparison of simulation results with different speeds

圖7 穩態下喉栓處于不同位置的壓強變化Fig．7 Pressure distribution of steady condition

表2 穩態和非穩態下的噴管入口靜壓對比Table 2 Comparison of nozzle inlet pressure between steady and unsteady condition

對比圖6、圖7及由表2可見，在喉栓完全未介入噴管喉部的情況下，非穩態與穩態的壓強分布基本相同;在喉栓調節介入的中間位置和完全介入時，非穩態分析得出的壓強明顯小于穩態下的壓強值?？梢?，調節速度較低時，喉栓運動對流場造成的干擾較小;隨著調節速度的增加，壓強的建立會產生延遲現象，而且調節速度越大，壓強建立的延遲情況越明顯。造成延遲的原因是喉栓調節運動對流場產生擾動，使流場分布還不能及時做出相應變化時，喉栓又運動到了下一節點，從而導致壓強的建立存在延遲現象。

由圖5可知，喉栓調節速度越快，喉栓調節時間越短。喉栓對流場的擾動時間就越短，在工程實際中，對流場擾動時間越短越好，在對材料、結構等方面不會造成太大影響的前提下，選擇一個適中的速度，已達到快速高效的推力可控。而喉栓處于同一位置時，壓強的增量隨著速度的增大而減小，這就需要對壓強的增量進行修正，才能清楚地得出不同喉栓調節速度與壓強之間的變化關系，進而實現推力的隨機控制。

3 結論

(1)將動網格技術應用到喉栓式推力可調噴管內流場動態特性分析中，可準確地模擬出流場的瞬時變化。

(2)隨著喉栓的介入，喉栓頭部會產生激波，對噴管流動產生顯著影響，且當喉栓完全介入時，影響最大。

(3)喉栓調節速度很小時，流場的分布與穩態下大體相同，隨著調節速度增大，流場演化過程有一定的延遲。

(4)在不影響噴管結構前提下，喉栓調節速度越快，對流場造成擾動的時間越短，入口壓強響應時間越短，越有利于實現快速的推力可控。

［1］ 李娟，李江，王毅林，等．喉栓式變推力發動機性能研究［J］．固體火箭技術，2007，30(6):505-509．

［2］ 滑利輝，田維平，甘曉松，等．喉栓式推力可調發動機噴管流場數值模擬［J］．固體火箭技術，2008，31(4):344-349．

［3］ 胡博．非同軸喉栓式推力可調噴管內流場及熱應力分析［D］．西北工業大學，2013．

［4］ 江帆，陳維平，李元元．基于射流與兩相流的射流曝氣器研究［J］．流體機械，2005，33(6):18-21．

［5］ 楊建明，吳建華．動網格技術數值模擬挑流沖刷過程［J］．水動力學研究與進展，2001，16(2):156-161．

［6］ Moureau V，Lartigue G，Sommerer Y，et al．Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids［J］．Journal of Computational Physics，2005，202:710-736．

［7］ Kris R，Jan V，Erik D．An arbitrary Lagrangian-Eulerian finite-volume method for the simulation of rotary displacement flow［J］．Applied Numerical Mathematics，2002，32:419-433．