劉鑫鵬，許曉勇，程圣清

0 引 言

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

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

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

1 計算方法

1.1 控制方程

實際的火箭發(fā)動機噴管內(nèi)的燃?xì)饬鲃邮謴?fù)雜，而附面層只是靠近壁面極薄的一層。為便于分析，將流動簡化為二維軸對稱流動，忽略壓力在垂直于壁面方向上的變化，并不計徹體力和源項，則對應(yīng)的控制方程為

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

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

圖1 計算坐標(biāo)系Fig.1 Calculation Coordinate

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

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

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

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

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

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

式中 mI''，mE''分別為內(nèi)、外邊界的質(zhì)量轉(zhuǎn)移率，與內(nèi)、外邊界流函數(shù)在流向上的梯度值有關(guān)；下標(biāo)I，E分別為內(nèi)外邊界

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

式中 a，b，c，d為與計算位置、工質(zhì)物性等相關(guān)的常數(shù)。

1.2 數(shù)值方法

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

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

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

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

式中 A，B，C為與上游物理量相關(guān)的常數(shù)。

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

1.3 近壁面處理

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

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

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

將壁面切應(yīng)力公式轉(zhuǎn)換到壁面坐標(biāo)系可表示為

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

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

Fortran程序?qū)⒆羁拷诿娴膬蓚€網(wǎng)格點連線的中點作為結(jié)合點，其速度可以取為臨近兩點的算術(shù)平均，并根據(jù) R e= U+y+和式（13）聯(lián)立求解，進(jìn)而可表示出壁面切應(yīng)力和摩擦阻力系數(shù)，再求出新的摩擦速度。如此反復(fù)迭代可得到近壁面較準(zhǔn)確的速度分布。對能量方程的處理方法與此類似，只是參照的無量綱數(shù)不同。

為了準(zhǔn)確計算靠近壁面的各物理量，本文還引入了滑移點[9]，其主要思想是假定物理量在靠近壁面區(qū)域內(nèi)滿足指數(shù)規(guī)律，利用插值和Taylor級數(shù)展開寫出結(jié)合點處的變量值和一階導(dǎo)數(shù)值，從而求出滑移點處的變量值。

在對壁面區(qū)域的處理中，F(xiàn)ortran程序還引入了一個子函數(shù)用于對壁面網(wǎng)格進(jìn)行調(diào)整，其主要作用是通過在壁面刪去或插入網(wǎng)格點，保證結(jié)合點處的壁面y+在一個合理的范圍內(nèi)。對于本文采用的壁面處理方法，一般控制y+在0～1范圍內(nèi)，保證在Couette流方程中的湍流粘性項可以忽略。

1.4 程序說明

根據(jù)以上數(shù)值方法，參考STAN5程序的結(jié)構(gòu)，編制了采用積分方法計算噴管附面層的Fortran程序。該程序由一個主程序驅(qū)動6個子程序分別完成相應(yīng)計算，其結(jié)構(gòu)框圖如圖3所示。

圖3 積分方法程序結(jié)構(gòu)Fig.3 Structure Diagram of the Program Using Integral Method

由圖3可知，求解函數(shù)STEP內(nèi)包含滑移點計算、坐標(biāo)變換、附面層參數(shù)計算、方程求解和初始化 5個模塊，根據(jù)調(diào)用時輸入?yún)?shù)的不同跳轉(zhuǎn)到相應(yīng)模塊運行。

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

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

2 計算示例與結(jié)果分析

2.1 計算示例

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

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

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

2.2 結(jié)果分析

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

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

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

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

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

續(xù)圖6

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

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

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

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

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

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

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

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

3 結(jié) 論

經(jīng)過以上計算分析，可以得到以下結(jié)論：

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

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

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