液體火箭發動機噴管附面層計算

劉鑫鵬，許曉勇，程圣清

0 引 言

噴管是火箭發動機推力室的重要組成部分，其功能是將高溫燃氣工質加速、熱能轉化為動能、高速排出以產生推力。燃氣在噴管內流動會有各種損失，主要為化學不平衡損失、非軸向流損失、氣流與壁面的摩擦損失及通過壁面的傳熱損失[1]。其中，化學不平衡損失是噴管的主要損失，原因是噴管內的復合反應滯后于流動狀態的變化，難以控制和補償；非軸向流損失和傳熱損失在采用了再生冷卻的特型噴管中影響極小。摩擦損失對發動機比沖有一定的影響，它與貼近壁面的附面層直接相關。此外，由于附面層內流速迅速降低至壁面處的靜止，使得壁面處的流通能力降低，導致按非粘性流設計的推力最優噴管型面實際工作時偏離設計狀態，造成比沖性能損失。因此，對噴管附面層的計算以及基于附面層的噴管型面修正對于摩擦損失的估算和控制以及比沖性能的提升具有重要意義。

火箭發動機噴管中的流動包括亞聲速、跨聲速和超聲速等多種情況，在附面層中存在著化學反應、物質交換、對流傳熱等過程，計算十分困難。Evans[2]、Nickerson[3]對發動機噴管性能進行了準確的計算，但程序規模大且輸入輸出困難，因此在計算精度要求不高時可對各種影響因素作一定取舍，同樣可以得到在工程研制中有用的結果；Coats等[4]開發的VIPER程序是對TDK中BLM模塊的改進，重點關注工質的流動，適用于大面積比和縱向曲率較大的噴管附面層計算；黃崇錫[5]提出的大縱向曲率的附面層理論主要考慮了型面的縱向曲率對流動的影響，在計算中忽略其它過程；Crawford等[6]開發的STAN5程序中考慮了壁面附近存在的物質交換，并能求解對應的物質交換方程；于勝春等[7]采用積分方法求解湍流附面層控制方程組，考慮了噴管內的對流換熱問題，得到的結果與 Bartz公式的計算結果較為吻合。

本文采用 Patankar-Spalding積分方法求解動量方程和能量方程，以STAN5程序為基礎，編制了用于計算火箭發動機噴管附面層的Fortran程序，并對某型號發動機型面進行了附面層計算，得到的結果與CFD軟件計算結果接近。

1 計算方法

1.1 控制方程

實際的火箭發動機噴管內的燃氣流動十分復雜，而附面層只是靠近壁面極薄的一層。為便于分析，將流動簡化為二維軸對稱流動，忽略壓力在垂直于壁面方向上的變化，并不計徹體力和源項，則對應的控制方程為

式中 U，V分別為x，y方向的速度；r為到對稱軸的距離；P為靜壓；h，h?分別為工質的靜焓和總焓；ρ，k，c， μ分別為工質的密度、熱導率、比熱和動力粘度；u′，v′，h′分別為速度和總焓由于湍流造成的波動項。

本文中采用的是曲線坐標系，x方向為沿壁面方向，y方向為垂直壁面方向，如圖1所示。

圖1 計算坐標系Fig.1 Calculation Coordinate

動量方程中的速度脈動項可以通過引入渦粘性系數Mε和湍流粘度tμ加以描述，即：

渦粘性系數Mε可由混合長度理論計算得到，即：

式中 l為混合長度，可由下式計算：

式中 κ為馮·卡門常數，取值為 0.4～0.41；D為范德列斯特衰減函數； y+(νw/ν)為到壁面的無量綱距離；A+為等效粘性底層厚度，取值為25。

為表述方便，引入等效粘度μeff，即分子動力粘度和湍流粘度之和，則動量方程中粘性力項可表示為μeff(? U/?y)。同理，能量方程中的總焓脈動項則引入湍流普朗特數 P rt表示，其值通過經驗公式計算，引入等效普朗特數 P reff可對方程進一步化簡。

本文引入了流函數ψ使得 ρ r U =?ψ / ? y成立，通過流函數可使連續方程自動成立并將方程從x?y坐標系轉換到x?ψ坐標系。通過Pantankar-Spalding坐標轉換，引入 ω = ( ψ ?ψI) /(ψE? ψI)將方程進一步轉換到x?ω坐標系。經過上述變換后的動量方程和能量方程如下：

式中 mI''，mE''分別為內、外邊界的質量轉移率，與內、外邊界流函數在流向上的梯度值有關；下標I，E分別為內外邊界

從式（7）、式（8）可以看出，整理后的動量方程和能量方程形式上一致，僅兩式的未知數不同，因此可以將兩式統一表示為變量為φ的形式：

式中 a，b，c，d為與計算位置、工質物性等相關的常數。

1.2 數值方法

本文采用有限體積法求解上述得到的偏微分方程，如圖2所示。

圖2 計算選取的控制體積Fig.2 Control Volume for Calculation

由圖2可以看出，xu，xd表示壁面上相鄰計算截面的位置，垂直壁面方向劃分為若干節點。本文采用沿x方向推進的計算方法，xu截面的所有物理量作為已知，將 xd截面的物理量作為未知量。

選取圖2所示的控制體積，對式（9）中的各項進行從xu到xd、ωi?1/2到ωi+1/2的二重積分，并采用文獻[8]中的方法建立對應的差分方程。這種通過積分得到差分方程的方法，比采用Taylor級數展開進行差分在防止誤差和滿足守恒方程物理概念方面具有明顯的區別，所得到的差分方程是一個六點隱式差分格式，對于任意的推進步長都是穩定的[8]。對差分方程進行整理，經過化簡可以得到如下的有限差分形式：

式中 A，B，C為與上游物理量相關的常數。

對dx截面所有節點進行計算，可以得到一個線性方程組，未知數的系數組成一個三對角線矩陣，在給定邊界值的條件下通過迭代法可以求解。

1.3 近壁面處理

以上構造有限差分方程的過程中假定所求變量在ω方向上是線性變化的，而噴管流動中靠近壁面的區域，速度、總焓等梯度較大，這一假定是不準確的，需要另作處理。本文在近壁面設置一個結合點，在靠近主流區域內求解上述有限差分方程，在近壁面區域單獨求解，并使兩者的解在結合點處相匹配，即對應的無量綱相似準則數相等。

以動量方程為例，本文將噴管近壁面處的流動按Couette流動處理，將式（2）轉化為二維平板流動的動量方程。引入摩擦速度 Uτ=，將各物理量均轉化為對應的無量綱量，將控制方程沿垂直壁面方向積分，并忽略徹體力影響和積分后的小量，可以得到：

式中+τ，wV+，U+，P+，y+分別為無量綱剪切應力、壁面 y向速度、沿壁面方向速度、壓力和垂直壁面距離。

將壁面切應力公式轉換到壁面坐標系可表示為

式中 μ+為無量綱動力粘度，為有效動力粘度與壁面分子動力粘度之比。

在靠近壁面的區域內，流體的分子粘性作用占主導，因此可以假設μ+=1，由式（11）、式（12）可得到關于U+的一階線性常微分方程，利用邊界條件y+= 0 時 U+= 0 對其進行求解并化簡得到：

Fortran程序將最靠近壁面的兩個網格點連線的中點作為結合點，其速度可以取為臨近兩點的算術平均，并根據 R e= U+y+和式（13）聯立求解，進而可表示出壁面切應力和摩擦阻力系數，再求出新的摩擦速度。如此反復迭代可得到近壁面較準確的速度分布。對能量方程的處理方法與此類似，只是參照的無量綱數不同。

為了準確計算靠近壁面的各物理量，本文還引入了滑移點[9]，其主要思想是假定物理量在靠近壁面區域內滿足指數規律，利用插值和Taylor級數展開寫出結合點處的變量值和一階導數值，從而求出滑移點處的變量值。

在對壁面區域的處理中，Fortran程序還引入了一個子函數用于對壁面網格進行調整，其主要作用是通過在壁面刪去或插入網格點，保證結合點處的壁面y+在一個合理的范圍內。對于本文采用的壁面處理方法，一般控制y+在0～1范圍內，保證在Couette流方程中的湍流粘性項可以忽略。

1.4 程序說明

根據以上數值方法，參考STAN5程序的結構，編制了采用積分方法計算噴管附面層的Fortran程序。該程序由一個主程序驅動6個子程序分別完成相應計算，其結構框圖如圖3所示。

圖3 積分方法程序結構Fig.3 Structure Diagram of the Program Using Integral Method

由圖3可知，求解函數STEP內包含滑移點計算、坐標變換、附面層參數計算、方程求解和初始化 5個模塊，根據調用時輸入參數的不同跳轉到相應模塊運行。

主程序的流程如圖4所示。從圖4中可以看出，由于動量方程和能量方程具有相同的形式，在求解過程中只需分別對相應系數進行計算即可對兩個方程共同求解。此外，積分程序中不存在迭代過程，完成一個截面的計算后判斷下一截面是否仍在計算范圍內，若在則進行推進求解計算，一直到噴管出口。

圖4 積分方法主程序流程Fig.4 Flow Diagram of the Main Program using Integral Method

2 計算示例與結果分析

2.1 計算示例

本文采用以上的計算方法，以STAN5程序為基礎，編制了用于計算火箭發動機噴管附面層的Fortran程序，對某型號發動機的噴管型面進行了計算，所計算的噴管擴張比為80，入口總壓為11.5 MPa。為計算簡便，工質的物性由多項式擬合給出，計算中根據變截面一維定常絕能流模型給定了壁面上82個點分別對應的主流速度等參數，從噴管入口處開始推進計算，推進過程中的主流速度等參數通過已知點參數線性插值得到。

由于噴管型面尚無附面層相關試驗數據，因此Fortran程序計算結果主要與CFD軟件計算的結果進行比較和分析。本文采用 Fluent14.5對噴管型面進行仿真計算，噴管的流動可簡化為二維軸對稱流動，湍流模型選用Spalart-Allmaras單方程模型。為了保證附面層計算的精度，該模型要求近壁面的y+控制在 3～5左右，因此需要對計算用網格進行多次加密，最終網格數量約54萬個，計算工質為氫氣和水蒸汽的混合物，兩者物質的量之比為3∶10，采用理想氣體模型，氫氣和水蒸汽各自的物性采用 CFD軟件自帶的多項式擬合。

計算完成后，根據附面層位移厚度的定義和動量厚度的定義計算各截面的位移厚度和動量厚度，其中δ為附面層厚度。在自編程序中利用循環結構進行積分直接輸出厚度結果，在CFD軟件計算中需先導出計算截面的速度、密度等參數，利用Excel數據處理軟件計算。本文從噴管喉部到出口選取了3個截面的計算結果，對Fortran程序和CFD軟件在同樣位置的計算結果進行了對比。

2.2 結果分析

圖5給出了喉部、距喉部軸向距離為0.8 m的中間截面和出口3個截面上速度分布的結果對比。

圖5 3個截面上的速度分布計算結果曲線Fig.5 Velocity Distribution on Three Calculated Cross Sections

由圖5可知，通過Fortran程序和Fluent計算得到的同一截面速度分布規律一致，從3個速度分布圖的橫坐標的變化可以看出噴管的附面層隨 x的增加逐漸變厚。

與附面層相關的3種厚度沿壁面的變化如圖6所示，選取的3個截面的厚度結果如表1所示。從圖6可知，Fortran程序和CFD軟件的計算結果差距很小，而且利用 Fortran程序可以直接求出各截面的 3種厚度，不像采用CFD軟件計算時還要進行數據的后續處理，提高了計算效率。

圖6 附面層厚度、位移厚度和動量厚度沿壁面的變化曲線Fig.6 Boundary Layer Thickness, Displacement Thickness and Momentum Thickness Along the Wall

續圖6

表1 3個截面上的附面層厚度、位移厚度和動量厚度Tab.1 Boundary Layer Thickness, Displacement Thickness and Momentum Thickness on Three Cross Sections

由圖 6可知，附面層在噴管圓柱段有所加厚，進入收斂段后迅速變薄并在喉部厚度最小，可忽略不計，進入擴張段后厚度以基本穩定的速度不斷增加。因此在對噴管附面層的計算中，可以認為型面段的附面層是從喉部開始建立，并且在擴張段逐漸加厚。

由表 1可知，對于本文計算的面積較大的噴管，在出口截面的附面層位移厚度達到厘米量級，對噴管壁面的流通能力將造成影響，應該對該設計型面進行修正，從而提高噴管的實際性能。在工程實踐中一般利用附面層的位移厚度對噴管型面進行修正。表2給出了該型面修正前后利用CFD軟件計算得到的真空推力和比沖，修正后比沖提高0.87 s。

由于進行了近壁面處理，因此采用本文程序計算不僅可以描述附面層內的速度分布和厚度變化，還可直接給出與壁面摩擦等相關的物理量沿壁面的變化，這是利用CFD軟件計算難以輕松得到的。

表2 修正前后發動機真空推力與比沖Tab.2 Vacuum Thrust and Specific Impulse of Engine Before and After Modification

壁面切應力沿壁面的變化規律如圖7所示。由圖7可知，由于附面層的存在導致的壁面摩擦在收斂段內迅速增強，在近喉部區域內壁面切應力達到峰值；進入擴張段后壁面摩擦開始減弱。在型面段后壁面切應力的下降速度逐漸平緩，這與實際情況相符。

圖7 壁面切應力沿壁面的變化曲線Fig.7 Wall Shear Stress Along the Wall

綜上所述，采用CFD軟件和Fortran程序計算得到的結果十分接近。而采用CFD軟件計算網格數量巨大，為保證網格質量需在壁面劃分近10萬個點，要得到較好的收斂精度需要近4萬步迭代，用時超過10 h；采用程序計算可以在1 min內計算出壁面近4千個截面的參數，足夠描述附面層沿壁面的發展情況，并能直接給出位移厚度、壁面切應力等相關結果，計算時間大大縮短，效率顯著提高。因此在同樣精度下可以采用本文得到的程序進行附面層計算。

3 結 論

經過以上計算分析，可以得到以下結論：

a）本文編制了用于計算火箭發動機噴管附面層的Fortran程序，計算結果與CFD軟件計算結果接近，在同樣精度要求下可采用Fortran程序計算，大大縮短計算時間；

b）對于本文計算的面積比為 80的噴管，出口截面的附面層位移厚度達到厘米量級，采用位移厚度對設計型面進行修正后噴管真空比沖有所提高；

c）本文的計算結果還需要進一步的試驗數據進行驗證，同時，在計算中對湍流源項的模擬等方面可以進行改進，為工程實際中的噴管附面層計算和型面修正提供更可靠的參考。