單噴管火箭自由噴流噪聲數(shù)值模擬

陳勁松，何冠杰，吳新躍，賀建華，賈延奎

(北京航天發(fā)射技術(shù)研究所，北京100076)

0 引 言

大型火箭發(fā)射過程中，發(fā)射中心5 m范圍內(nèi)的噴流噪聲聲壓級(jí)(參考聲壓2×10-5Pa，文章下同)會(huì)超過170 dB，這種強(qiáng)噴流噪聲會(huì)造成火箭、發(fā)射系統(tǒng)的控制儀器、電器元器件性能顯著下降甚至失效，對(duì)火箭、發(fā)射系統(tǒng)的薄壁結(jié)構(gòu)件、焊縫也會(huì)發(fā)生聲疲勞破壞作用。為確保火箭安全發(fā)射，目前正在發(fā)展噴水降噪技術(shù)之類噴流噪聲控制技術(shù)。

發(fā)展火箭噴流噪聲控制技術(shù)，需要采取必要的手段預(yù)示噴流噪聲。過去限噴流噪聲理論限制，主要借助實(shí)物試驗(yàn)提煉半經(jīng)驗(yàn)式理論預(yù)示方法[1-4]。近十年來，受大規(guī)模并行計(jì)算技術(shù)帶動(dòng)，氣動(dòng)噪聲領(lǐng)域開始采用數(shù)值模擬技術(shù)預(yù)示氣動(dòng)噪聲，主要發(fā)展了兩類預(yù)示方法：一是以大渦模擬為代表的近似或直接數(shù)值模擬方法[5]；二是以弗克斯-威廉姆斯(Ffowcs-Williams)和豪金斯(Hawkings)發(fā)展的FWH方法為代表的混合數(shù)值模擬方法[6]，該方法在復(fù)雜工程領(lǐng)域應(yīng)用廣泛。

直接數(shù)值模擬方法完全依托高精度湍流流場(chǎng)數(shù)值模擬結(jié)果[7]，當(dāng)前計(jì)算資源、耗時(shí)特別嚴(yán)重，難以應(yīng)用于規(guī)模巨大、結(jié)構(gòu)極其復(fù)雜的運(yùn)載火箭發(fā)射噴流噪聲預(yù)示領(lǐng)域。例如，燃?xì)饬餍u尺度lv與來流雷諾數(shù)Re關(guān)系為：

lv≈DeRe-0.75

(1)

(2)

聲波尺度λ與噴口馬赫數(shù)Mae關(guān)系為：

(3)

則直接數(shù)值模擬需要的數(shù)值模擬網(wǎng)格量Nt為：

(4)

式(1)～(4)中：De為火箭發(fā)動(dòng)機(jī)噴口直徑；Ste為噴口斯特努哈爾數(shù)；ue，ρe，μe分別為火箭發(fā)動(dòng)機(jī)噴口燃?xì)饬魉俣取⒚芏取ば韵禂?shù)。式(1)～(4)相關(guān)參數(shù)量綱依托國(guó)際單位制確定，文章下同。

由式(1)～(4)不難推斷噴口直徑De=1 m、噴口馬赫數(shù)Mae=3.0的大型火箭，噴流噪聲直接數(shù)值模擬網(wǎng)格總量將不低于2×1012，這對(duì)大型火箭噴流噪聲預(yù)示而言，仍然是難以實(shí)現(xiàn)的目標(biāo)。

FWH方法主要用于求解遠(yuǎn)場(chǎng)噪聲輻射特性，該方法無法計(jì)及噪聲傳播途徑中的聲波相互干擾問題，對(duì)于火箭發(fā)射時(shí)多聲源產(chǎn)生的噴流噪聲傳播干擾，以及復(fù)雜發(fā)射系統(tǒng)結(jié)構(gòu)擾動(dòng)時(shí)造成的噴流噪聲傳播干擾問題無能為力。

近些年，采用有限元法數(shù)值模擬復(fù)雜聲源條件或復(fù)雜結(jié)構(gòu)擾動(dòng)的氣動(dòng)噪聲傳播問題取得很大進(jìn)展[8]。其中，比利時(shí)FFT(Free Field Technology)公司與空客公司合作，利用伽遼金(Galerkin)有限元方法開展了民航氣動(dòng)噪聲研究，取得了開創(chuàng)性研究結(jié)果[9]；我國(guó)南京航空航天大學(xué)Lv等[10]利用伽遼金有限元方法取得了機(jī)翼擾動(dòng)聲場(chǎng)研究結(jié)果。本文利用伽遼金有限元方法開展了單噴管火箭自由噴流噪聲數(shù)值模擬研究，下面介紹具體方法及研究結(jié)果。

1 噴流噪聲數(shù)值模擬研究方法

噴流噪聲數(shù)值模擬研究方法如圖1所示。

圖1 噴流噪聲數(shù)值模擬研究方法框圖Fig.1 Method diagram of numerical simulation study on jet noise

圖1所示噴流噪聲數(shù)值模擬研究方法中，非定常燃?xì)饬鲌?chǎng)數(shù)值模擬是噴流噪聲數(shù)值模擬的基礎(chǔ)和前提，但其依據(jù)的模型和噴流噪聲數(shù)值模擬模型必須聯(lián)合開發(fā)，如圖2、圖3所示。

圖2 燃?xì)饬鲌?chǎng)數(shù)值模擬模型Fig.2 Numerical simulation model of jet flow field

圖3 噴流噪聲場(chǎng)數(shù)值模擬模型Fig.3 Numerical simulation model of jet noise field

從圖2、圖3可以看出，噴流噪聲數(shù)值模擬模型與燃?xì)饬鲌?chǎng)數(shù)值模擬模型的主要差別在于是否包含燃?xì)饬骱诵膮^(qū)，噴流噪聲數(shù)值模擬模型不包含燃?xì)饬骱诵膮^(qū)，主要原因是當(dāng)前超聲速氣流噪聲預(yù)示方法尚不成熟，發(fā)展的伽遼金有限元方法、FWH方法適用于亞聲速氣流噪聲傳播、聲輻射特性研究，也具有統(tǒng)一、規(guī)范的理論基礎(chǔ)。燃?xì)饬鲌?chǎng)及噴流噪聲場(chǎng)實(shí)際建模過程中，規(guī)劃圖2、圖3所示燃?xì)饬骱诵膮^(qū)域包絡(luò)了全部超聲速流動(dòng)區(qū)域，也包絡(luò)了剪切速率極高的自由噴流邊界層以及下游部分高亞聲速(Mae≥0.8)區(qū)域，主要目的是立足亞聲速燃?xì)饬髁鲃?dòng)區(qū)域數(shù)值模擬研究噴流噪聲分布特性及形成機(jī)理，使數(shù)值模擬結(jié)果具有一定可信度。

自由噴流非定常燃?xì)饬鲌?chǎng)數(shù)值模擬程序主要基于計(jì)算流體力學(xué)發(fā)展的相對(duì)成熟的Roe-FDS隱式算法實(shí)現(xiàn)，該程序給出燃?xì)饬骱诵膮^(qū)包絡(luò)邊界上壓力、速度、密度脈動(dòng)信息，作為聲源邊界條件，驅(qū)動(dòng)后續(xù)噴流噪聲聲傳播數(shù)值模擬。噴流噪聲聲傳播數(shù)值模擬基于歐拉(Euler)方程開展，歐拉方程微分形式為：

(5)

采用高階伽遼金法數(shù)值模擬噴流噪聲聲傳播規(guī)律時(shí)，將式(5)轉(zhuǎn)換為有限元積分形式：

(6)

其中,Γα標(biāo)序α的網(wǎng)格單元型函數(shù)，應(yīng)用時(shí)采用拉格朗日插值函數(shù)，其表述形式為：

(7)

式(5)～(7)中：φ為流場(chǎng)參數(shù)，代表密度、動(dòng)量或能量參數(shù)；Fi為反映經(jīng)過網(wǎng)格單元的流動(dòng)特性函數(shù)，具體表述由歐拉方程展開形式得到；t為時(shí)間；xi為標(biāo)序i的坐標(biāo)方向；ni為網(wǎng)格單元法向單位矢量沿i向坐標(biāo)分量；V為標(biāo)序α的網(wǎng)格單元的體積；S為標(biāo)序α的網(wǎng)格單元的表面積；r為標(biāo)序α的網(wǎng)格單元內(nèi)空間任意位置坐標(biāo)的矢量距離；rk，rj分別為標(biāo)序k，j節(jié)點(diǎn)位置坐標(biāo)的矢量距離；x，y，z為直角坐標(biāo)系的三個(gè)坐標(biāo)分量；m為標(biāo)序α的網(wǎng)格單元節(jié)點(diǎn)數(shù)。

基于式(7)，同樣的網(wǎng)格單元，伽遼金有限元方法可以利用更多節(jié)點(diǎn)，使得噴流噪聲數(shù)值模擬模型網(wǎng)格單元總數(shù)可以減少很多。實(shí)踐發(fā)現(xiàn)，這種單元數(shù)也必須合理控制，當(dāng)網(wǎng)格單元數(shù)過少造成網(wǎng)格分辨率不足時(shí)，噴流噪聲聲壓分布數(shù)值模擬結(jié)果存在局部失真，如圖4所示。圖4中為方便說明，將核心區(qū)燃?xì)饬鲌?chǎng)分布云圖與噴流噪聲場(chǎng)聲壓云圖一并顯示。

后續(xù)研究總結(jié)提出了網(wǎng)格分辨率控制方法，主要是控制相同空間區(qū)域內(nèi)噴流噪聲數(shù)值模擬模型網(wǎng)格單元數(shù)Nt_FEA與燃?xì)饬鲌?chǎng)數(shù)值模擬網(wǎng)格單元數(shù)Nt_CFD關(guān)系：

Nt_CFD=(5～50)Nt_FEA

(8)

2 噴流噪聲分布規(guī)律

利用圖2所示模型開展的非定常燃?xì)饬鲌?chǎng)數(shù)值模擬結(jié)果顯示：高溫燃?xì)饬髯試姽苤饾u向下游推進(jìn)，伴隨推進(jìn)過程燃?xì)饬骷げㄏ抵饾u形成，推進(jìn)前鋒受環(huán)境氣流阻滯呈現(xiàn)球冠狀抽卷狀態(tài)，如圖5所示。

圖5 燃?xì)饬魍七M(jìn)過程靜溫分布云圖Fig.5 Static temperature cloud map of advancing combustion gas flow

基于瞬態(tài)燃?xì)饬鲌?chǎng)，噴流噪聲數(shù)值模擬給出了如圖6、圖7所示的噴流噪聲分布規(guī)律。

圖6 50.0 Hz頻段噴流噪聲聲壓分布云圖Fig.6 Sound pressure cloud map of jet noise at 50.0 Hz

圖7 150.0 Hz頻段噴流噪聲聲壓分布云圖Fig.7 Sound pressure cloud map of jet noise at 150.0 Hz

圖5～圖7顯示：伴隨燃?xì)饬魍七M(jìn)，燃?xì)饬髑颁h、燃?xì)饬鞯饶軈^(qū)域末稍以及兩者之間的燃?xì)饬鲃?dòng)區(qū)域形成了三個(gè)強(qiáng)噪聲源，強(qiáng)噪聲源驅(qū)動(dòng)下，圖6、圖7中下方區(qū)域發(fā)展成了略傾斜向下的強(qiáng)噪聲傳播區(qū)域；三個(gè)強(qiáng)噪聲傳播區(qū)域中，尤以位置靠上的兩個(gè)強(qiáng)噪聲區(qū)域表現(xiàn)強(qiáng)勁，它們?cè)谙蜻h(yuǎn)處傳播過程中，50 Hz以下頻段逐漸融合成一體，而150 Hz以上頻段相對(duì)獨(dú)立。

隨著燃?xì)饬髑颁h越過圖6、圖7所示的流場(chǎng)包絡(luò)底邊界并遠(yuǎn)離燃?xì)饬鞯饶軈^(qū)域末稍，圖6、圖7所示強(qiáng)噪聲傳播區(qū)域相應(yīng)變化，如圖8所示。

圖8 250.0 Hz頻段噴流噪聲聲壓分布云圖Fig.8 Sound pressure cloud map of jet noise at 250.0 Hz

圖8顯示：燃?xì)饬髑颁h推進(jìn)較遠(yuǎn)距離后，原先跟隨高溫燃?xì)饬髑颁h的高聲強(qiáng)噪聲源及其傳播區(qū)域繼續(xù)向下游推進(jìn)，強(qiáng)度依然較高；而燃?xì)饬鞯饶軈^(qū)末稍的噪聲源強(qiáng)度相對(duì)下降，變成次強(qiáng)噪聲源，由此形成的傳播區(qū)域變成次強(qiáng)噪聲傳播區(qū)域。圖8還顯示計(jì)算域聲源邊界下方，受燃?xì)饬髑颁h推進(jìn)驅(qū)動(dòng)生成新的強(qiáng)噪聲源，由此產(chǎn)生沿燃?xì)饬髑颁h推進(jìn)方向傳播的強(qiáng)噪聲區(qū)域。

圖9 100.0 Hz頻段噴流噪聲聲壓分布云圖Fig.9 Sound pressure cloud map of jet noise at 100.0 Hz

燃?xì)饬鲌?chǎng)逐漸穩(wěn)定后，噴流噪聲分布如圖9所示。圖9顯示：穩(wěn)定燃?xì)饬鲌?chǎng)條件下，燃?xì)饬鞯饶軈^(qū)末稍附近強(qiáng)噪聲源及其傳播區(qū)域再次凸顯，并且還伴隨有次強(qiáng)噪聲源及其傳播區(qū)域。圖8中聲源計(jì)算域邊界下方的強(qiáng)噪聲源及其傳播區(qū)域在圖9中仍保持存在，但強(qiáng)度較燃?xì)饬鞯饶苣┥詤^(qū)域附近噪聲強(qiáng)度要低些。

仔細(xì)分析圖7～圖9，還可以發(fā)現(xiàn)強(qiáng)噪聲傳播區(qū)域之間還存在其它強(qiáng)度略低的噪聲傳播區(qū)域。針對(duì)噴流噪聲原理性試驗(yàn)[4]專題數(shù)值模擬結(jié)果顯示，小尺度試驗(yàn)中這些強(qiáng)度較低的噴流噪聲有時(shí)會(huì)表現(xiàn)得比較突出。例如噴口馬赫數(shù)2.5的小尺度試驗(yàn)噴流噪聲分布如圖10所示。

圖10 小尺度試驗(yàn)100.0 Hz頻段噪聲聲壓云圖Fig.10 Sound pressure cloud map of scaled jet test noise at 100.0 Hz

對(duì)比9、圖10可以看出：原先大尺度噴管下方表現(xiàn)不明顯的噪聲傳播區(qū)域在小尺度試驗(yàn)中表現(xiàn)得比較突出;同時(shí)，小尺度試驗(yàn)燃?xì)饬鞯饶軈^(qū)末稍附近強(qiáng)噪聲源相對(duì)傳播區(qū)域也增加較多。

3 噴流噪聲機(jī)理

圖6～圖8中強(qiáng)噪聲主要位于高溫燃?xì)饬髑颁h和燃?xì)饬鞯饶軈^(qū)末稍之間的空間區(qū)域，該空間區(qū)域燃?xì)鈭?chǎng)流線分布如圖11所示。

圖11 局部燃?xì)饬鲌?chǎng)流線分布圖Fig.11 Streamline distribution map of local jet field

圖11顯示：受高溫燃?xì)饬髑颁h推進(jìn)過程中卷吸帶動(dòng)作用，在高溫燃?xì)饬髑颁h和燃?xì)饬鞯饶軈^(qū)末稍之間的空間區(qū)域形成了十分明顯的卷吸大渦。非定常燃?xì)饬鲌?chǎng)數(shù)值模擬結(jié)果表明這種大渦是不穩(wěn)定的，它隨著燃?xì)饬髑颁h繼續(xù)推進(jìn)不斷翻滾、脫落，又不斷形成，正是這種不穩(wěn)定卷吸渦驅(qū)動(dòng)作用形成了圖6～圖8所示的強(qiáng)噪聲源及其傳播區(qū)域。

圖6～圖8所示燃?xì)饬鞯饶軈^(qū)末稍附近的強(qiáng)噪聲在燃?xì)饬鲌?chǎng)穩(wěn)定后仍然得到保持，如圖9、圖10所示。形成燃?xì)饬鞯饶軈^(qū)末稍附近強(qiáng)噪聲機(jī)理可結(jié)合圖12、圖13所示燃?xì)饬鲃?dòng)及分布特性說明。

圖12 噴管下方空間區(qū)域流線分布圖Fig.12 Streamline distribution map of local jet field below the nozzle

圖13 燃?xì)饬鲝?qiáng)湍流區(qū)附近噪聲聲壓云圖Fig.13 Sound pressure cloud map of jet noise near the strong turbulence

圖12顯示：一級(jí)發(fā)動(dòng)機(jī)燃?xì)饬鲃?dòng)接近穩(wěn)定狀態(tài)后，燃?xì)饬鞯饶軈^(qū)域附近并不存在類似圖11所示的流動(dòng)大渦，因此，圖9、圖10中位置靠上的兩個(gè)強(qiáng)噪聲源及強(qiáng)噪聲傳播區(qū)并不是由燃?xì)饬髁鲃?dòng)大渦引起。圖13顯示：強(qiáng)噪聲源形成與燃?xì)饬魍牧鲝?qiáng)度關(guān)系密切，圖9、圖10中位置靠上的兩個(gè)強(qiáng)噪聲源分別緊鄰燃?xì)饬鲝?qiáng)湍流區(qū)末稍附近以及湍流強(qiáng)度梯度劇烈變化區(qū)域。湍流強(qiáng)度指標(biāo)反映了燃?xì)饬髁鲃?dòng)的脈動(dòng)特性，湍流強(qiáng)度越強(qiáng)表示燃?xì)饬髅}動(dòng)越強(qiáng)，據(jù)此說明燃?xì)饬髅}動(dòng)較強(qiáng)區(qū)域以及燃?xì)饬髅}動(dòng)變化劇烈區(qū)域均是燃?xì)饬鲝?qiáng)噪聲生成區(qū)域。

對(duì)于圖10所示小尺度試驗(yàn)噴管下方強(qiáng)噪聲區(qū)域形成機(jī)理，可結(jié)合燃?xì)饬鲌?chǎng)核心區(qū)流動(dòng)特性說明，參考圖14、圖15所示。

圖14 噴管下方附近噪聲聲壓及燃?xì)饬黛o溫分布云圖Fig.14 Sound pressure and static temperature cloud map below the nozzle

圖15 燃?xì)饬鲃?dòng)對(duì)稱軸線上靜壓變化曲線Fig.15 Static pressure curve along the symmetric axis of jet

圖14是圖10局部放大圖。圖15為燃?xì)饬魍七M(jìn)軸線上靜壓變化曲線圖，圖中橫坐標(biāo)表示燃?xì)饬髯試娍谙蛳掠瓮七M(jìn)的無量綱距離，縱坐標(biāo)表示無量綱燃?xì)饬黛o壓，無量綱距離、燃?xì)饬黛o壓分別相對(duì)噴口直徑、噴口燃?xì)饬鲃?dòng)壓計(jì)及。圖14、圖15顯示燃?xì)饬鲃?dòng)核心區(qū)存在靜溫與靜壓突變。從靜壓分布曲線能夠分辨三處壓力突變情況，說明這種小尺度試驗(yàn)存在三組激波系，這三組激波系附近空間正是圖10中噴管下方強(qiáng)噴流噪聲源形成空間，可以推斷：燃?xì)饬鞔┻^激波系過程中，流動(dòng)狀態(tài)發(fā)生的劇烈擾動(dòng)變化是形成噴管下方強(qiáng)噪聲的主要原因之一。圖13已經(jīng)說明噴管下方激波系附近湍流脈動(dòng)強(qiáng)度也比較強(qiáng)，因此燃?xì)饬鞯膹?qiáng)湍流脈動(dòng)因素也對(duì)噴管下方強(qiáng)噴流噪聲源形成起了重要貢獻(xiàn)。

4 數(shù)值模擬試驗(yàn)驗(yàn)證

噴流噪聲數(shù)值模擬方法及具體數(shù)值模擬結(jié)果依據(jù)研發(fā)的噴流噪聲模擬試驗(yàn)系統(tǒng)進(jìn)行了驗(yàn)證研究，文獻(xiàn)[4]已經(jīng)對(duì)此進(jìn)行了系統(tǒng)介紹，本文限篇幅不再說明，僅給出利用該試驗(yàn)開展的自由噴流噪聲試驗(yàn)結(jié)果和數(shù)值模擬結(jié)果對(duì)比典型情況。由于試驗(yàn)?zāi)Ｐ椭袦y(cè)點(diǎn)和數(shù)值模擬模型中網(wǎng)格單元節(jié)點(diǎn)難以保持一一對(duì)應(yīng)關(guān)系，為此在擇取一個(gè)試驗(yàn)測(cè)點(diǎn)數(shù)據(jù)時(shí)，相應(yīng)擇取了該測(cè)點(diǎn)上下兩個(gè)相鄰網(wǎng)格單元節(jié)點(diǎn)預(yù)示數(shù)據(jù)。擇取過程中，自由噴流噪聲試驗(yàn)?zāi)Ｐ秃蛿?shù)值模擬模型均利用直角坐標(biāo)系x向、y向(高度方向)、z向等三向無量綱坐標(biāo)標(biāo)識(shí)空間具體位置，噴口中心設(shè)定為坐標(biāo)原點(diǎn)。實(shí)例中，自由噴流噪聲試驗(yàn)測(cè)點(diǎn)坐標(biāo)為(8.78,-1.22,0.00)，數(shù)值模擬模型中上部相鄰網(wǎng)格單元節(jié)點(diǎn)坐標(biāo)為(8.80,-1.15,0.00)，下部相鄰網(wǎng)格單元節(jié)點(diǎn)坐標(biāo)為(8.80,-1.25,0.00)。數(shù)值模擬監(jiān)測(cè)點(diǎn)與試驗(yàn)監(jiān)測(cè)點(diǎn)噴流噪聲聲壓級(jí)對(duì)比曲線如圖16所示。

圖16 噴流噪聲數(shù)值模擬結(jié)果與測(cè)試結(jié)果對(duì)比曲線Fig.16 The contrast SPL curve between numerical simulation and test of jet noise

圖16顯示：在600 Hz～2500 Hz范圍內(nèi)，噴流噪聲聲壓級(jí)數(shù)值模擬曲線與試驗(yàn)測(cè)試曲線趨勢(shì)吻合，數(shù)值模擬結(jié)果與試驗(yàn)測(cè)試結(jié)果相對(duì)接近，差別大致在±5.0 dB范圍內(nèi)，數(shù)值模擬結(jié)果可以作為分析噴流噪聲的依據(jù)；在2500 Hz～5000 Hz范圍內(nèi)，數(shù)值模擬結(jié)果嚴(yán)重偏離測(cè)試值。造成這種嚴(yán)重偏離現(xiàn)象的原因與噴流噪聲數(shù)值模擬模型網(wǎng)格單元分辨率仍然不足有關(guān)，后續(xù)仍需持續(xù)探索改進(jìn)數(shù)值模擬方法。

5 結(jié) 論

1) 文章綜合非定常燃?xì)饬鲌?chǎng)數(shù)值模擬方法、伽遼金有限元聲傳播數(shù)值模擬方法以及FWH聲輻射數(shù)值模擬方法，實(shí)現(xiàn)并完成了單噴管火箭自由噴流噪聲數(shù)值模擬。

2) 數(shù)值模擬結(jié)果顯示單噴管火箭燃?xì)饬魍七M(jìn)初期，卷吸前鋒帶動(dòng)產(chǎn)生的卷吸渦附帶產(chǎn)生了強(qiáng)噴流噪聲；強(qiáng)噴流噪聲會(huì)跟隨燃?xì)饬髑颁h持續(xù)向下游移動(dòng)；燃?xì)饬鲌?chǎng)相對(duì)穩(wěn)定后，強(qiáng)噴流噪聲位于等能區(qū)末稍附近，主要由燃?xì)饬鲝?qiáng)湍流脈動(dòng)及湍流脈動(dòng)劇烈變化引起；一些小尺度試驗(yàn)中，激波擾動(dòng)以及強(qiáng)湍流脈動(dòng)共同作用還會(huì)造成緊鄰噴管下方空間區(qū)域存在強(qiáng)噴流噪聲傳播區(qū)域。

3) 受網(wǎng)格分辨率限制，當(dāng)前中低頻段噪聲聲壓級(jí)數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果總體比較接近，高頻段偏差較大，數(shù)值模擬方法仍需改進(jìn)。