榮吉利，阿尼蘇，程修妍，范博超

( 北京理工大學(xué) 宇航學(xué)院，北京 100081)

近年來，隨著火箭導(dǎo)彈以及航天飛行助推器的廣泛應(yīng)用[1]，火箭發(fā)動機的高溫超音速噴流噪聲問題成為了不斷被研究的對象，噴流噪聲對工作人員、地面設(shè)施及有效載荷極具威脅，因此深入研究噴流噪聲具有重要意義.

1954年，Lighthill[2]提出了聲類比理論，從理論上完成了對噴流發(fā)聲問題的研究. 隨著計算機技術(shù)的發(fā)展，目前數(shù)值模擬已成為噴流噪聲問題的主要研究方法. Gamet和Estivalezes[3]于1998年對馬赫數(shù)2，雷諾數(shù)6×104的超音速高溫噴流進(jìn)行了數(shù)值模擬，計算結(jié)果與試驗結(jié)果相差不大. 2006年，劉耀峰和徐文燦等[4]數(shù)值模擬了噴口附近的流場結(jié)構(gòu)及渦系結(jié)構(gòu)等流場特性，計算結(jié)果與實驗結(jié)果吻合較好. 2017年，馬紅鵬等[5]對設(shè)計的某微型燃燒室的穩(wěn)態(tài)燃燒過程進(jìn)行數(shù)值模擬，得到燃燒室穩(wěn)態(tài)工作時的速度場、溫度場及各組分濃度分布，仿真與實驗數(shù)據(jù)契合度高. 孫博等[6]于2018年采用3維黏彈性有限元方法分析了固體火箭發(fā)動機點火發(fā)射時藥柱界面裂紋的穩(wěn)定性. Troyes和Vuillot[7]于2016年采用大渦模擬方法和Ffowcs Williams-Hawkings方程(簡稱FW-H方程)分別計算了馬赫數(shù)為3.1的高溫超音速噴流的流場和聲場，計算結(jié)果誤差較小. 2019年，Langenais[8]改進(jìn)了Troyes和Vuillot[7]的遠(yuǎn)場噪聲求解方法，大大提高了計算結(jié)果的精度(幾乎所有測點總聲壓級的誤差小于1 dB). 目前廣泛使用的固體火箭發(fā)動機會產(chǎn)生凝相固體顆粒，但幾乎所有關(guān)于固體火箭發(fā)動機噴流噪聲的研究都忽略了固體顆粒對噴流流場和聲場的影響.

文中針對Seiner&Ponton噴管模型[9]，選取了6種不同的湍流模型，包括Spalart-Allmatas模型、Standardk-ε模型、RNGk-ε模型、Realizablek-ε模型、Standardk-ω模型和SSTk-ω模型，進(jìn)行了穩(wěn)態(tài)計算，研究了這6種湍流模型對噴流流場穩(wěn)態(tài)結(jié)果的影響；又以RNGk-ε模型的穩(wěn)態(tài)計算結(jié)果為初始流場，使用大渦模擬方法(large Eddy simulation，簡稱LES)計算了該噴管模型的非穩(wěn)態(tài)結(jié)果，并使用FW-H方程進(jìn)行了聲場分析，計算結(jié)果與試驗結(jié)果進(jìn)行了對比；最后在此模型的基礎(chǔ)上，添加了固體顆粒，研究固體顆粒對噴流流場及聲場的影響.

1 數(shù)值模擬方法

1.1 控制方程

在笛卡爾坐標(biāo)系下，3維可壓縮流的N-S方程為

(1)

式中，

(2)

式中：ρ為氣體密度；u，v，w分別為x，y，z方向的速度；p為壓強；qi為熱傳導(dǎo)引起的能量通量；τij為應(yīng)力項.

(6)

式中μ為黏性系數(shù).

單位氣體總能量e為

(7)

式中γ為氣體比熱.

最后再引入理想氣體狀態(tài)方程，

p=ρRT,

(8)

式中:R為氣體常數(shù);T為溫度. 即構(gòu)成了封閉的控制方程組.

1.2 湍流模型

Spalart-Allmatas模型是求解湍流黏性輸運方程的單方程模型. 此模型適用于低雷諾數(shù)模型，不適用于自由剪切流. Spalart-Allmatas模型的輸運方程形式為

(9)

Standardk-ε模型對湍動能(k)和耗散率(ε)建立了輸運方程. Standardk-ε模型不適用于強分離流、包含大曲率的流動和強壓力梯度流動，適用于高雷諾數(shù)湍流，是工程計算中應(yīng)用十分廣泛的湍流模型. 湍動能(k)和耗散率(ε)從下列輸運方程中得出，

Gk+Gb-ρε-YM，

(10)

(11)

式中:Gk為平均速度梯度產(chǎn)生的湍動能;Gb為由浮力產(chǎn)生的湍動能;YM為可壓縮湍流中脈動膨脹對總耗散率的貢獻(xiàn)量;C1ε、C2ε和C3ε是常數(shù);σk和σε分別是k和ε的湍流Prandtl數(shù).

RNGk-ε模型是Standardk-ε模型的改進(jìn)模型，RNGk-ε模型在ε方程中增加了修正項，提高了計算快速應(yīng)變流的精度；考慮了湍流漩渦；為湍流Prandtl數(shù)提供了一個解析公式；使得RNGk-ε模型比Standardk-ε模型在計算更廣泛的流動中有更高的可信度和精度.

Realizablek-ε模型也是Standardk-ε模型的改進(jìn)模型，Realizablek-ε模型為湍流黏度提供了可選公式；修正了耗散率ε的輸運方程. Realizablek-ε模型的不足是：在計算同時包含旋轉(zhuǎn)和靜態(tài)流動區(qū)域時不能提供自然的湍流黏度.

Standardk-ω模型基于Wilcoxk-ω模型，考慮了低雷諾數(shù)、可壓縮性和剪切流傳播；該模型的缺點是對來流因素過于敏感；該模型適用于墻壁束縛流動和自由剪切流動；其湍動能(k)和比耗散率(ω)的輸運方程為

(12)

(13)

式中Γk和Γω分別為k和ω的有效擴散率.

SSTk-ω模型在近壁面使用Standardk-ω模型，而在主流區(qū)域使用k-ε模型；在ω方程中加入了阻尼交叉擴散導(dǎo)數(shù)項并且在湍流黏度的定義中包含了湍流剪切應(yīng)力的傳遞；這些改進(jìn)使得SSTk-ω模型比Standardk-ω模型有更高的精度和可靠性.

LES是近幾十年才發(fā)展起來的一個流體力學(xué)中重要的數(shù)值模擬研究方法；其基本思想是使用過濾方程濾掉小尺度渦，使用統(tǒng)一模型計算小尺度渦，再使用N-S方程直接求解大尺度渦. 假設(shè)用uf來表示流場速度，

(14)

(15)

為了使方程組封閉，需要使用Smagorinsky模型來建立亞格子應(yīng)力的數(shù)學(xué)模型.

1.3 數(shù)值格式

本文采用密度基求解器來求解控制方程，差分格式選用了Roe通量差分格式.

密度基求解器首先同時求解連續(xù)方程、動量方程、能量方程及組分輸運方程，再逐一求解其他標(biāo)量方程. 密度基通常用于高速可壓縮流動的求解.

為減少數(shù)值耗散造成的誤差，準(zhǔn)確捕捉激波現(xiàn)象需要采用高階或者高分辨率的差分格式進(jìn)行離散，Roe格式是通量差分分裂最具代表性的格式，它以其優(yōu)秀的間斷分辨率，得到了廣泛應(yīng)用[11].

1.4 FW-H方程

1969年，F(xiàn)fowcs Williams和Hawkings[12]將噴流發(fā)聲的理論擴展到了任意運動固體邊界在流動中發(fā)聲的問題，即為FW-H方程. 其方程如下,

(16)

式中：ui為xi方向的流速大小;un為積分面法向流速大小;vi為xi方向面速度大小;vn為積分面法向面速度大小;δ(f)為狄拉克(Dirac delta)函數(shù);H(f)為亥維賽(Heaviside)函數(shù);Tij為萊特希爾(Lighthill)應(yīng)力張量;Pij為可壓縮應(yīng)力張量.

FW-H方程可以計算在流體中流體與運動的物體相互作用而發(fā)出聲音的問題.

1.5 離散相模型

固體發(fā)動機噴流中含有固體顆粒，在流場計算時應(yīng)考慮固體顆粒的輻射作用，但目前幾乎所有關(guān)于固體火箭發(fā)動機噴流噪聲的研究都忽略了固體顆粒對噴流流場的影響. 在進(jìn)行氣-固兩相流場的耦合計算時，氣相的控制方程相較于只有氣相的控制方程，僅多了氣固兩相相互作用的源項，

式中:Ep為粒子的輻射率;a為介質(zhì)的吸收系數(shù);ap為粒子的吸收系數(shù);n為介質(zhì)的折射率;σ為玻爾茲曼常數(shù);G為入射輻射.

通過對粒子受力平衡等式的積分來求解粒子軌跡線

(18)

1.6 湍流強度

湍流強度是一種度量氣流脈動程度的針對湍流統(tǒng)計規(guī)律的描述，用脈動速度均方根與時均速度之比表示為

I=urms/utm.

(19)

2 數(shù)值仿真結(jié)果分析及驗證

2.1 噴管噴流數(shù)值仿真模型

Seiner&Ponton噴管模型[9]的出口直徑Dj=91.44 mm，出口馬赫數(shù)為2. 由于文獻(xiàn)中沒有詳細(xì)介紹此噴管的幾何尺寸，因此參照王克印等[13]提出的噴管設(shè)計方法設(shè)計了最簡單的縮擴型噴管. 具體模型如圖1所示，為了超音速噴流可以充分發(fā)展以及選取合理的聲源積分面，因此計算模型的軸向長度必須足夠長；同樣徑向也需要足夠?qū)捯苑乐惯吔鐥l件會對流場造成影響，影響計算結(jié)果. 結(jié)合以上因素的考慮，將計算模型的軸向長度設(shè)為50Dj，徑向長度設(shè)為30Dj.

模型中包含入口邊界兩個，即燃?xì)馊肟诤涂諝馊肟? 按照實驗條件，給定燃?xì)馊肟?即噴管入口)處總溫1 375 K，總壓780 000 Pa；空氣入口總溫300 K，總壓101 325 Pa. 出口均設(shè)置為壓力出口；聲源面設(shè)置為內(nèi)部邊界，聲源面形狀為逐漸擴張的圓管，沿計算域軸向的長度為25Dj，如圖1所示，聲源面左側(cè)距噴管出口面0.5Dj，左側(cè)面半徑為2Dj，右側(cè)面半徑為4Dj. 具體邊界條件劃分如圖1所示.

網(wǎng)格圖如圖2所示. 整個計算模型均采用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格最小尺度為0.5 mm，位于噴管壁面處，以便可以準(zhǔn)確模擬噴管內(nèi)部流動；在聲源面內(nèi)部加密網(wǎng)格，因為此處為湍流劇烈區(qū)域；聲源面外使用了逐漸粗糙的網(wǎng)格，可以減小計算量并減小出口界面?zhèn)畏瓷洳ǖ挠绊? 總網(wǎng)格數(shù)為283萬.

2.2 不同湍流模型穩(wěn)態(tài)計算結(jié)果對比

以此計算模型為基礎(chǔ)，分別使用6種不同的湍流模型，包括Spalart-Allmatas模型、Standardk-ε模型、RNGk-ε模型、Realizablek-ε模型、Standardk-ω模型及SSTk-ω模型，計算了Seiner&Ponton噴管模型的穩(wěn)態(tài)流場.

圖3展示了不同湍流模型的穩(wěn)態(tài)數(shù)值結(jié)果中軸向?qū)ΨQ界面流場馬赫數(shù)分布云圖，幾種湍流模型都能計算出噴流的壓縮膨脹波系，且膨脹、壓縮的位置幾乎相同，但是噴流的膨脹壓縮波系的波節(jié)數(shù)以及湍流核心區(qū)域的長度存在著明顯的不同；幾種湍流模型中，RNGk-ε模型得到的湍流核心區(qū)域最長，且波節(jié)最多最為明顯；Spalart-Allmatas模型核心區(qū)域最短，波節(jié)最少.

圖4分別對比了不同湍流模型計算結(jié)果軸線上的馬赫數(shù)、壓強以及溫度，不同湍流模型的計算結(jié)果總體趨勢大致相同，馬赫數(shù)及壓強在約10Dj前呈現(xiàn)震蕩變化趨勢，隨后快速下降. 其中，Spalart-Allmatas模型的衰減最快，經(jīng)歷3個周期后便開始衰減，而RNGk-ε模型衰減最慢，約經(jīng)歷5個周期后衰減.

Varnier[14]通過一系列試驗得出的超音速段長度經(jīng)驗公式為

(20)

式中Maj為出口馬赫數(shù).

表1展示了6種湍流模型計算結(jié)果流場的超音速段長度，針對此模型，式(20)計算結(jié)果為18.21Dj，RNGk-ε與Standardk-ε的計算結(jié)果與式(20)計算結(jié)果最接近. 再考慮到RNGk-ε為Standardk-ε的改進(jìn)模型，接下來將以RNGk-ε計算的穩(wěn)態(tài)流場為初始流場來繼續(xù)進(jìn)行計算.

表1 不同湍流模型的超音速段長度

圖5給出了以RNGk-ε模型計算結(jié)果的軸向?qū)ΨQ截面流場各參數(shù)分布云圖，噴流的流場內(nèi)部是由一系列膨脹壓縮波組成的.

2.3 LES非穩(wěn)態(tài)計算結(jié)果

以RNGk-ε模型的計算的穩(wěn)態(tài)結(jié)果為初始流場，使用LES湍流模型計算了該噴管模型的非穩(wěn)態(tài)流場，使用的時間步長為5×10-6s.

圖6給出了非穩(wěn)態(tài)計算至0.2 s時，模型軸向?qū)ΨQ截面流場的各參數(shù)分布云圖. 與穩(wěn)態(tài)計算結(jié)果(如圖5所示)對比可得，前5個波系結(jié)構(gòu)，流場結(jié)構(gòu)幾乎一致，但隨后氣流變得雜亂無章.

對非穩(wěn)態(tài)計算結(jié)果軸線上的速度進(jìn)行了統(tǒng)計平均，再將軸線上平均速度分布以噴管出口速度vj為參考得到了軸線上量綱一速度分布，以便與試驗數(shù)據(jù)進(jìn)行對比. 圖7展示了軸線上的量綱一速度分布，與Seiner等[9]的試驗數(shù)據(jù)進(jìn)行了對比，數(shù)值計

算結(jié)果與試驗數(shù)據(jù)十分接近，證明了數(shù)值計算結(jié)果的合理性.

2.4 聲場計算結(jié)果

非穩(wěn)態(tài)計算至噴流充分發(fā)展(湍流特性完全顯現(xiàn))后，在非穩(wěn)態(tài)流場的基礎(chǔ)上繼續(xù)進(jìn)行計算，并保存所選取的聲源面上的流場脈動信息(每10個時間步保存一次信息)，再使用FW-H方程預(yù)測遠(yuǎn)場聲學(xué)特性. 將聲場計算結(jié)果與Seiner等[9]的試驗數(shù)據(jù)進(jìn)行對比，驗證計算結(jié)果的正確性.

參照Seiner等[9]試驗方案中布置噪聲測點的位置，在此計算模型中設(shè)置了噪聲測點，設(shè)置噪聲測點的具體方式是：如圖8所示，以噴管出口中心點為圓心，在模型的軸對稱平面上取一半徑為40Dj(3.66 m)的半圓，設(shè)噴流方向為0°，在此半圓上，從0°～180°每隔10°便設(shè)置一個噪聲測點，共設(shè)置19個噪聲測點.

將20°和50°噪聲測點的1/3倍頻程譜與試驗數(shù)據(jù)進(jìn)行對比，如圖9所示. 計算結(jié)果與試驗數(shù)據(jù)的變化趨勢大體一致，最大誤差約為5 dB，在工程可接受的范圍內(nèi).

計算各個噪聲測點處的總聲壓級，并將計算結(jié)果與Seiner等[9]的試驗數(shù)據(jù)對比，如圖10所示，在遠(yuǎn)場半徑為40Dj的半圓上，聲壓級最大值為149.94 dB，位于50°的測點處，試驗數(shù)據(jù)的最大值也大約出現(xiàn)在50°測點處，計算結(jié)果與試驗數(shù)據(jù)的變化趨勢相同，60°～100°之間，數(shù)值計算結(jié)果與試驗數(shù)據(jù)基本一致，在20°～60°之間，數(shù)值計算結(jié)果存在一定的誤差，最大誤差約為6 dB. 但在20°～60之間數(shù)值計算的結(jié)果與Biancherin[15]使用FW-H計算的結(jié)果十分接近，差別在2 dB以內(nèi)，且比Biancherin[15]的計算結(jié)果更為接近試驗數(shù)據(jù).

誤差出現(xiàn)原因主要有：①數(shù)值計算模型的環(huán)境條件與試驗的環(huán)境條件存在差別，因此會出現(xiàn)一定的誤差；②此噴管出口噴流馬赫數(shù)較高，對網(wǎng)格的質(zhì)量要求較高，可以進(jìn)一步加密網(wǎng)格來減小誤差；③由于沒有Seiner等[9]試驗噴管的詳細(xì)結(jié)構(gòu)，數(shù)值計算所用噴管模型是自行設(shè)計的構(gòu)型簡單的噴管，會對流場計算結(jié)果產(chǎn)生影響，也導(dǎo)致了一定的誤差.

但綜合流場對比、單點聲壓級頻程譜對比以及測點總聲壓級對比來說，數(shù)值計算的流場、聲場結(jié)果還是比較準(zhǔn)確的.

3 固體顆粒對噴流流場、聲場的影響

3.1 固體顆粒對噴流流場的影響

以RNGk-ε模型的計算結(jié)果為初始流場，以噴管出口面為固體顆粒釋放面，進(jìn)行穩(wěn)態(tài)計算，獲取固體顆粒對噴流流場的影響，固體顆粒相關(guān)參數(shù)如表2所示.

表2 固體顆粒參數(shù)Tab.2 Solid particle parameters

表中:ρp為顆粒密度;Dp為顆粒直徑;Cp為顆粒的定壓比熱;ep為輻射率;qmrp為顆粒的質(zhì)量流量.

對固體顆粒做出如下假設(shè)[16]：①氣體與顆粒之間不考慮化學(xué)反應(yīng)；②顆粒間無碰撞等相互作用；③顆粒形狀為球形.

含固體顆粒的噴流流場明顯與無顆粒噴流流場不同，含固體顆粒的噴流的超音速段長度為16.51Dj，較無固體顆粒的超音速段長度減小了2.6Dj，圖11分別對比了有無固體顆粒仿真模型的計算結(jié)果軸線上的馬赫數(shù)、壓強以及溫度，加入固體顆粒后，燃?xì)馀c固體顆粒相互作用，使得燃?xì)鈬娏髋蛎涀龉υ龃?，?dǎo)致噴流的膨脹壓縮波數(shù)減少，馬赫數(shù)迅速衰減，溫度升高.

3.2 固體顆粒對噴流聲場的影響

以固體顆粒的穩(wěn)態(tài)流場為初始流場，使用LES計算該模型的非穩(wěn)態(tài)流場，計算足夠時間后，保存所選取的聲源面上的流動參數(shù)，再使用FW-H方程預(yù)測遠(yuǎn)場聲學(xué)特性.

在此計算模型中設(shè)置了與無固體顆粒非穩(wěn)態(tài)計算模型位置完全相同的聲源面及噪聲測點以方便對比，圖12對比了有無固體顆粒時遠(yuǎn)場噪聲測點總聲壓級分布的情況.

0°～30°范圍內(nèi)，有固體顆粒噪聲測點的總聲壓級大于無固體顆粒噪聲測點的總聲壓級. 0°～30°處于噴流的下游，噴流下游噪聲的主要成分為湍流混合噪聲[17]. 湍流混合噪聲與湍流強度成正比例增加，因此對比了有無固體顆粒噴流軸線的湍流強度，如圖13所示，在27Dj處，有固體顆粒噴流的湍流強度達(dá)到了0.8，遠(yuǎn)大于無固體顆粒噴流的湍流強度，因此導(dǎo)致了在噴流下游處有固體顆粒聲壓級較大.

40°～180°范圍內(nèi)，有固體顆粒噪聲測點的總聲壓級低于無固體顆粒噪聲測點的總聲壓級. 40°～

180°為噴流的中上游區(qū)域，噪聲的主要成分為寬頻激波噪聲[17]，寬頻激波噪聲與激波強度成正比例增加[18]. 有固體顆粒噴流的馬赫數(shù)小于無固體顆粒噴流的馬赫數(shù)，且超音速段長度較短，即有固體顆粒噴流的激波強度較小，導(dǎo)致了噴流中上游區(qū)域噪聲較小.

4 結(jié) 論

針對Seiner&Ponton噴管模型，建立了3維噴流數(shù)值模擬方法，進(jìn)行了遠(yuǎn)場噪聲分析，并考慮了固體顆粒對噴流流場和聲場的影響，得到了以下結(jié)論：

① RNGk-ε模型最適合用于穩(wěn)態(tài)噴流數(shù)值模擬；

② 使用FW-H方程計算遠(yuǎn)場噪聲的誤差有一定的準(zhǔn)確性；

③ 固體顆粒的存在減小了噴流流場的超音速段長度，增大了膨脹壓縮半徑，提高了噴流溫度；

④ 在噴流下游區(qū)域，由于固體顆粒的存在，增大了噴流的湍流強度，導(dǎo)致噴流下游區(qū)域噪聲增大，在噴流中上游區(qū)域，減小了噴流的激波強度較小，導(dǎo)致中上游區(qū)域噪聲降低.