馬博平, 王剛, 雷知錦, 葉正寅

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072)

聲爆的抑制問題是新一代環(huán)保型超聲速客機研發(fā)過程中亟需解決的關(guān)鍵技術(shù)瓶頸之一[1]。精確地評估超聲速客機飛行時地面上感受到的聲爆噪聲水平是通過低聲爆設(shè)計降低聲爆的先決條件。因此,聲爆的預(yù)測問題一直是全球航空航天研究機構(gòu)關(guān)注的熱點之一。2008年,美國航空航天管理局(the national aeronautics and space administration,NASA)發(fā)布了用于聲爆近場預(yù)測技術(shù)驗證與確認(rèn)的多組標(biāo)準(zhǔn)算例,并組織了研討會論證超聲速客機近場聲爆預(yù)測技術(shù)的發(fā)展?fàn)顩r與未來趨勢[2]。在此基礎(chǔ)上,美國航空航天學(xué)會(American Institute of Aeronautics and Astronautics, AIAA)又分別于2014年和2017年組織了第一、二屆國際聲爆預(yù)測研討會(The 1stand 2ndSonic Boom Prediction Workshop, SBPW-1[3]、SBPW-2[4-5]),并發(fā)布了更加系統(tǒng)的測試算例,包括中國在內(nèi)的多國學(xué)者參會交流了其在近場和遠(yuǎn)場聲爆預(yù)測方面的最新進(jìn)展。

主流的聲爆預(yù)測方法分為2步:首先通過風(fēng)洞實驗或計算流體力學(xué)(computational fluid danymics,CFD)方法計算飛行器周圍較小范圍的過壓(over-pressure)分布,再以此為輸入根據(jù)Whitham[6]和Hayes[7]發(fā)展的準(zhǔn)線性理論或非線性Burgers方程[8]計算地面聲爆波形,最后,計算地面聲爆的感覺噪聲級(PLdB)。因此,提高近場聲爆的預(yù)測精度是進(jìn)行全場聲爆預(yù)測的關(guān)鍵。

近場聲爆預(yù)測時往往需要考慮超聲速流動特有的一些特點以便保證預(yù)測結(jié)果的可靠性,這其中就包括特殊的網(wǎng)格生成技術(shù)。國外針對近場聲爆預(yù)測的網(wǎng)格生成技術(shù)已有一些研究。Cliff 和 Thomas在其研究中分別使用了結(jié)構(gòu)網(wǎng)格和各向同性非結(jié)構(gòu)網(wǎng)格[9]。Choi,Alonso和van der Weide在流場中壓力梯度大的區(qū)域使用了各向同性自適應(yīng)網(wǎng)格技術(shù),在一倍機體以內(nèi)表現(xiàn)出良好的效果[10]。Carter等研究了非結(jié)構(gòu)網(wǎng)格分辨率對聲爆近場預(yù)測的影響,并指出采用沿來流馬赫角對齊的各向異性網(wǎng)格可以顯著改善計算結(jié)果[11]。Campbell等對上述方法進(jìn)行了驗證與確認(rèn)[12]。Meredith在其Shaped Sonic Boom Demonstrator(SSBD)項目研究中在笛卡爾坐標(biāo)系下采用了嵌套網(wǎng)格技術(shù)及沿激波角對齊的網(wǎng)格策略,取得了良好的效果[13]。Laflin等在模型周圍采用自適應(yīng)各向同性網(wǎng)格,在此之外則采用沿來流馬赫角對齊的結(jié)構(gòu)網(wǎng)格捕捉流場激波,二者之間通過插值實現(xiàn)流場變量的傳遞,由此提出了一種對于復(fù)雜構(gòu)型快速網(wǎng)格生成方法[14]。

國內(nèi)對于這方面的研究相對較少。馮曉強采用與文獻(xiàn)[11]實質(zhì)相同的混合網(wǎng)格方法使用CFD分析得到的近場壓力分布代替原始的F函數(shù)和等效截面積分布作為匹配目標(biāo)值,提出了改進(jìn)的F函數(shù)方法[15],并在此基礎(chǔ)上基于線化理論SGD(seebass-george-darden)反設(shè)計方法對聲爆/氣動一體化設(shè)計方法進(jìn)行了研究,但由于計算資源的限制,其采用了楔形網(wǎng)格進(jìn)行CFD計算[16]。徐悅等采用馬赫錐構(gòu)型的結(jié)構(gòu)化網(wǎng)格研究了黏性對聲爆近場預(yù)測結(jié)果的影響[17]。王剛等采用與文獻(xiàn)[11]一致的非結(jié)構(gòu)混合網(wǎng)格研究了其自研求解器對于聲爆標(biāo)模近場聲爆預(yù)測的可靠性及精度[18],在此基礎(chǔ)上研究了數(shù)值離散格式和是否考慮黏性對近場聲爆預(yù)測結(jié)果的影響和以此為輸入得到的遠(yuǎn)場波形的差異,并定性分析了其對人的感受性的影響[19]。

從現(xiàn)有研究結(jié)果來看,目前國外學(xué)者針對網(wǎng)格量和網(wǎng)格拓?fù)浣Y(jié)構(gòu)對近場波形的影響開展了一系列的研究工作[20],但對結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格、不同網(wǎng)格類型和網(wǎng)格分布等因素對于近場聲爆預(yù)測結(jié)果的影響缺少系統(tǒng)研究。為此,本文選取旋成體雙錐和帶翼錐體構(gòu)型采用數(shù)值方法分別就網(wǎng)格對齊方式、網(wǎng)格類型、網(wǎng)格生成方法和分布對聲爆近場預(yù)測的精度的影響進(jìn)行了分析和研究。最后,使用實驗數(shù)據(jù)齊全的LM1021標(biāo)模驗證了所得到的結(jié)論。

1 數(shù)值方法

本文使用自研非結(jié)構(gòu)混合網(wǎng)格復(fù)雜流場模擬程序HUNS3D[21]數(shù)值求解三維Euler/雷諾平均Navier-Stokes方程來實現(xiàn)近場聲爆的數(shù)值預(yù)測。

直角坐標(biāo)系下,采用Spalart-Allmars(SA)一方程湍流模型[22]封閉的積分形式的三維RANS方程為:

(1)

在HUNS3D中,流動控制方程通過格心格式有限體積法進(jìn)行離散。對流通量Fc可以采用包括ROE[23]、AUSM(advection upstream splitting method)系列格式[24]、熵守恒格式(entropy consistent,EC)[25]在內(nèi)的多種迎風(fēng)型格式和中心格式JST(jameson-schmidt-turkel)[26]進(jìn)行數(shù)值計算。在超聲速流動數(shù)值模擬中上述迎風(fēng)格式均需要和Venkatakrishnan斜率限制器[27]配合使用以保證總變差減小(TVD,total-variation-diminishing)特性。

在HUNS3D,黏性通量Fv采用中心格式進(jìn)行離散;湍流雷諾應(yīng)力可以選用包括SA一方程模型、Menter′s MSST(shear stress transport)兩方程模型[28]和DES(detached eddy simulation)模型等湍流模擬方法進(jìn)行封閉。時間離散采用LU-SGS(lower-upper symmetric Gauss-seidel)方法,該算法利用矩陣的LU分解避免了復(fù)雜的矩陣求逆運算,具備較高的計算效率。

2 網(wǎng)格生成方法及其影響

2.1 網(wǎng)格對齊方式的影響

聲爆是飛行器在超聲速飛行過程中產(chǎn)生的特定聲學(xué)現(xiàn)象,其本質(zhì)是超聲速飛行時飛行器周圍產(chǎn)生的激波,當(dāng)其傳播至地面時就形成了聲爆。因此,要準(zhǔn)確的對聲爆進(jìn)行預(yù)測,關(guān)鍵在于精確的捕捉流場中產(chǎn)生的激波信息。對于定常超聲速流動,激波角μ計算公式為

μ=sin-1(1/M)

(2)

式中,M為來流馬赫數(shù)。與傳統(tǒng)正交網(wǎng)格相比,馬赫錐構(gòu)型的網(wǎng)格計算時相對前者更能夠減小計算耗散,因而可以更為精準(zhǔn)的捕捉激波信息,進(jìn)而能夠更為精準(zhǔn)的預(yù)測聲爆。

本文基于雙錐構(gòu)型研究網(wǎng)格對齊方式的影響。雙錐模型如圖1所示,錐體長度為0.050 8 m(L=2 inches),沿軸向半徑分布如公式(3)所示, 風(fēng)洞試驗數(shù)據(jù)取自Carlson、Mark 和Morris于1965年在NASA蘭利研究中心4×4英寸超聲速風(fēng)洞的實驗結(jié)果。

(3)

本文數(shù)值模擬選取馬赫數(shù)為2.01,迎角0°,湍流模型為SA湍流模型,雷諾數(shù)為7.38×106/m。本算例中,分別使用了3套非結(jié)構(gòu)混合網(wǎng)格,其拓?fù)浣Y(jié)構(gòu)示意圖如圖2所示。模型壁面附近圓柱區(qū)域采用非結(jié)構(gòu)網(wǎng)格,非結(jié)構(gòu)網(wǎng)格對復(fù)雜外形具有適應(yīng)性強、劃分簡單的優(yōu)點,在外部拉伸區(qū)域采用結(jié)構(gòu)化網(wǎng)格,以便更好的捕捉激波。3套網(wǎng)格除拉伸區(qū)域傾斜角度外其他因素一律相同,圖中ψ表示網(wǎng)格拉伸方向與來流夾角,3套網(wǎng)格夾角分別為μ-5,μ和μ+5。其中,沿馬赫角對齊(ψ=μ)的網(wǎng)格示意圖見圖3。所有網(wǎng)格的邊界層第一層網(wǎng)格高度均滿足y+小于1,邊界層數(shù)40層,增長率1.21。網(wǎng)格節(jié)點數(shù)8 923 645,單元數(shù)10 303 712。

圖1雙錐模型示意圖 圖2 雙錐模型網(wǎng)格拓?fù)涫疽鈭D 圖3 雙錐模型網(wǎng)格示意圖

計算結(jié)果如圖4所示,在距離模型較近區(qū)域(h/L=1)處可以看出不同錐角計算結(jié)果基本完全重合,表明此時錐角對于捕捉到的激波形狀并沒有顯著影響,隨著距離模型增大,到h/L=5時,由于錐角不同產(chǎn)生的差異逐漸明顯,錐角為馬赫角的網(wǎng)格捕捉到的激波峰值和谷值范圍最大,其余2組網(wǎng)格捕捉到的波形明顯鈍化一些。除此之外,錐角為馬赫角的網(wǎng)格在x=17處捕捉到的激波上升處也更為迅速。圖4c)為h/L=10處的結(jié)果,圖中除不同錐角網(wǎng)格計算結(jié)果外,還對比了實驗值。由圖可以看出,在離模型較遠(yuǎn)距離處,使用錐角為馬赫角的網(wǎng)格計算得到的波形與實驗結(jié)果基本完全吻合,于此同時,其余2組計算結(jié)果與實驗結(jié)果均有明顯誤差,在激波上升處僅有錐角為馬赫角網(wǎng)格捕捉到的波形與實驗結(jié)果完全吻合。由此可以證明,采用馬赫錐構(gòu)型網(wǎng)格由于可以顯著減小激波傳播時數(shù)值耗散,有利于捕捉激波,因而更有利于提高聲爆預(yù)測的精度。

圖4 雙錐模型h/L=1,5,10處不同錐角近場過壓計算結(jié)果

2.2 結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格的影響

網(wǎng)格是進(jìn)行CFD流場計算的基礎(chǔ),合理的網(wǎng)格分布是正確模擬一些復(fù)雜流動的基本前提。網(wǎng)格主要分為結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格。結(jié)構(gòu)網(wǎng)格雖然具有可以控制網(wǎng)格拓?fù)涞膭澐趾途W(wǎng)格的正交性的優(yōu)點,但是對于復(fù)雜構(gòu)型來說,網(wǎng)格生成周期長甚至無法生成。而非結(jié)構(gòu)網(wǎng)格則具有很強的描述復(fù)雜幾何形狀構(gòu)型的能力,易于得到高質(zhì)量的網(wǎng)格,同時,非結(jié)構(gòu)網(wǎng)格生成過程快速且自動化程度高。

本文選取DWB模型來研究結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格對聲爆近場預(yù)測的影響,該模型選自Hunton等[29]于1973年的進(jìn)行風(fēng)洞實驗。DWB模型為面對稱模型,包括1個機身和2個三角翼構(gòu)型,模型參考長度為0.175 2 m,三角翼前緣后掠角為69°。該布局被廣泛用于聲爆研究,在第一屆AIAA聲爆預(yù)測會議中被選為測試模型之一。本文DWB構(gòu)型計算采用的是SA湍流模型,雷諾數(shù)為7.38×106/m。

圖5 DWB模型、2種計算域及近壁面網(wǎng)格示意圖

圖5a)～c)所示為DWB模型、計算域示意圖。根據(jù)拓?fù)浣Y(jié)構(gòu)可將網(wǎng)格分為模型附近帶圓柱區(qū)域和不帶圓柱區(qū)域。所謂圓柱區(qū)域,即為網(wǎng)格拓?fù)渲心Ｐ透浇Ｐ偷牟糠?。對于帶圓柱區(qū)域網(wǎng)格,為減小激波信號傳播過程中的耗散,圓柱應(yīng)盡可能小地包圍模型。圓柱區(qū)域以外采用馬赫錐構(gòu)型,即網(wǎng)格沿馬赫角由圓柱表面往外拉伸,本算例中拉伸部分網(wǎng)格為六面體單元。圖5d)～f)分別為采用非結(jié)構(gòu)混合網(wǎng)格近場、帶圓柱內(nèi)場結(jié)構(gòu)網(wǎng)格近場和不帶圓柱結(jié)構(gòu)網(wǎng)格近場示意圖。網(wǎng)格數(shù)據(jù)如表1所示。

表1 不同網(wǎng)格類型網(wǎng)格量示意圖

圖6a)～d)分別不同網(wǎng)格拓?fù)浣Y(jié)構(gòu)在不同周向角(φ)下近場聲爆過壓計算結(jié)果。由圖可以看出,不同網(wǎng)格拓?fù)浣Y(jié)構(gòu)捕捉到的激波形狀均與實驗吻合較好,不同結(jié)果之間的區(qū)別主要在過壓峰值處。有圓柱區(qū)域結(jié)構(gòu)網(wǎng)格和混合網(wǎng)格計算結(jié)果幾乎完全一致,這得益于2組網(wǎng)格具有同樣的拓?fù)浣Y(jié)構(gòu),同時,除圓柱內(nèi)區(qū)域以外,其余部分網(wǎng)格分布也基本相同。對于2組結(jié)構(gòu)網(wǎng)格,無圓柱區(qū)域網(wǎng)格計算結(jié)果與實驗值偏差較大,這主要是由于其在模型附近網(wǎng)格拉伸明顯,引入了更大的耗散,如圖6c)所示。因此,對于用于聲爆近場預(yù)測的楔形拓?fù)浣Y(jié)構(gòu)的網(wǎng)格,用混合網(wǎng)格能夠取得更好的效果。盡管帶圓柱區(qū)域結(jié)構(gòu)網(wǎng)格具有類似的效果,但生成結(jié)構(gòu)網(wǎng)格相較于非結(jié)構(gòu)混合網(wǎng)格需要花費更長時間,因此非結(jié)構(gòu)混合網(wǎng)格可以取得快速高效的效果,應(yīng)用更為方便。

圖6 使用結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格計算的DWB近場過壓預(yù)測結(jié)果與實驗結(jié)果的比較

除此之外,由圖6a)～d)可以明顯看出隨著方位角的變化,尾部激波從單個激波演變?yōu)槎鄠€激波,數(shù)值模擬捕捉到的演變程與實驗曲線表現(xiàn)的變化規(guī)律基本一致。

2.3 網(wǎng)格類型的影響

為研究圓柱區(qū)域外網(wǎng)格類型對聲爆波形傳播的影響,本文在DWB楔形拓?fù)浣Y(jié)構(gòu)基礎(chǔ)上在圓柱區(qū)域以外分別使用了四面體、六面體混合網(wǎng)格和全四面體網(wǎng)格,如圖7a)～b)所示。其中,四面體、六面體混合網(wǎng)格網(wǎng)格節(jié)點數(shù)為5 303 747,單元數(shù)為8 734 621,四面體混合網(wǎng)格網(wǎng)格節(jié)點數(shù)為3 172 122,單元數(shù)為9 755 505。

圖7 DWB模型不同網(wǎng)格類型示意圖

圖8a)～b)分別為四面體、六面體混合網(wǎng)格和全四面體網(wǎng)格計算得到的DWB構(gòu)型近場過壓結(jié)果和實驗值在不同周向角下的對比。由圖可以看出,六面體混合網(wǎng)格和四面體混合網(wǎng)格計算得到的近場過壓結(jié)果基本一致,且與實驗結(jié)果吻合較好,表明在網(wǎng)格量近似的情況下,網(wǎng)格類型對聲爆預(yù)測的影響有限。

2.4 網(wǎng)格分布的影響

在CFD的實際應(yīng)用中,計算精度和計算量往往是相互矛盾的關(guān)系,當(dāng)要求提高計算精度時,計算量往往會增大,反之亦然。因此如果可以將網(wǎng)格節(jié)點分布到更為流場中更為“合適”的地方,對提高計算效率可能會起到事半功倍的作用。對于本文中采用的馬赫錐構(gòu)型,當(dāng)模型軸向方向網(wǎng)格點給定后,為探究圓柱體周向和垂直模型軸方向的網(wǎng)格分布對聲爆預(yù)測的精度影響,本文利用DWB帶翼錐體模型分別生成了3套不同分布的網(wǎng)格進(jìn)行了計算,并將結(jié)果與實驗值進(jìn)行了比對,網(wǎng)格分布參數(shù)示意圖和不同展弦比網(wǎng)格信息分別如圖9和表2所示。

圖9 網(wǎng)格分布參數(shù)定義

表2 不同展弦比網(wǎng)格信息

2.4.1 網(wǎng)格展弦比的影響

圖10a)～b)為DWB模型在其中心線下方h=0.629 9 m不同周向角處結(jié)果。圖中計算都采用HUNS3D在黏性條件下計算,短折線代表網(wǎng)格展弦比為10的結(jié)果,點劃線代表的是展弦比為4的結(jié)果,雙點劃線代表網(wǎng)格展弦比為1的網(wǎng)格計算的結(jié)果,實驗測量結(jié)果用黑點表示?？梢园l(fā)現(xiàn),網(wǎng)格展弦比在一定范圍內(nèi)對近場聲爆預(yù)測基本沒有影響。不同展弦比的網(wǎng)格在不同高度不同周向角的條件下捕捉到的激波形狀都能夠很好吻合,同時,和實驗值也吻合很好。但將網(wǎng)格展弦比系數(shù)進(jìn)一步增大時,可以發(fā)現(xiàn)網(wǎng)格展弦比對捕捉到的激波峰值有明顯影響,見圖10c)～d)。

圖10 DWB模型不同展弦比網(wǎng)格計算的近場聲爆過壓預(yù)測結(jié)果與實驗結(jié)果對比

2.4.2 周向網(wǎng)格間距的影響

表3給出了不同周向間距的3套測試網(wǎng)格詳細(xì)信息。

表3 不同周向距離網(wǎng)格信息

圖11 DWB模型不同周向距離網(wǎng)格計算的近場聲爆過壓預(yù)測結(jié)果與實驗結(jié)果對比

圖11a)～b)為DWB模型在其中心線下方h=0.629 9 m不同周向角處結(jié)果。圖中計算都采用HUNS3D在黏性條件下計算,短折線代表網(wǎng)格內(nèi)部圓柱區(qū)域周向角度為8的結(jié)果,點劃線代表的是周向角度為4的結(jié)果,雙點劃線代表周向角度為2的網(wǎng)格計算的結(jié)果,實驗測量結(jié)果用黑點表示。不同周向角度網(wǎng)格的網(wǎng)格量信息見表3。

可以發(fā)現(xiàn),圓柱區(qū)域周向角度對捕捉到的激波“陡峭”程度有明顯影響,采用周向角度較小的網(wǎng)格捕捉到的激波形狀更為“尖銳”,并且和實驗值吻合更好。

以上結(jié)果表明,網(wǎng)格展弦比和馬赫錐面周向網(wǎng)格間距對流場中激波峰值和激波形狀的捕捉有明顯影響,使用相對合理的網(wǎng)格展弦比和馬赫錐面周向網(wǎng)格間距有利于提高近場聲爆預(yù)測精度。

3 復(fù)雜構(gòu)型近場聲爆預(yù)測

3.1 模型和網(wǎng)格生成

為驗證前文所述的網(wǎng)格生成策略在復(fù)雜構(gòu)型中的應(yīng)用效果,本文選取了SBPW-1中作為可選算例的標(biāo)模LM1021全機方案進(jìn)行測試。采用帶圓柱區(qū)域的混合網(wǎng)格拓?fù)浣Y(jié)構(gòu),其中圓柱內(nèi)周向角度為2,遠(yuǎn)場網(wǎng)格展弦比為1。圓柱區(qū)域以外采用沿馬赫角對齊六面體網(wǎng)格以便于捕捉聲爆傳播。

LM1021方案來自美國Lockheed-Martin公司,包括3個發(fā)動機短艙和1個V型尾翼。其參考長度為0.568 9 m,實驗數(shù)據(jù)測量條件為來流馬赫數(shù)1.6,迎角2.3°(模型在網(wǎng)格中已繞前緣點旋轉(zhuǎn),因此CFD求解中迎角為0°)。此算例中網(wǎng)格節(jié)點數(shù)為24 805 540,單元數(shù)為12 144 918。面網(wǎng)格數(shù)目為252 741,附面層33層,第1層高度滿足y+小于1。計算采用MSST湍流模型,雷諾數(shù)為7.38×106/m。模型及網(wǎng)格詳細(xì)示意圖如圖12a)～b)所示。

3.2 計算分析

圖13a)～f)為不同周向角下模型正下方h=0.807 7 m處CFD計算得到的LM1021近場過壓預(yù)測結(jié)果結(jié)果和實驗結(jié)果的對比。由圖可以看出,CFD計算結(jié)果準(zhǔn)確捕捉到了LM1021構(gòu)型流場中復(fù)雜的激波系,且與實驗測量的演化過程基本一致,表明所采用的網(wǎng)格策略有效。

圖12 LM1021構(gòu)型及網(wǎng)格示意圖

圖13 LM1021構(gòu)型CFD計算的近場聲爆過壓預(yù)測結(jié)果與實驗結(jié)果對比

4 結(jié) 論

本文使用自研流場求解器對3組典型標(biāo)模的近場聲爆過壓分布進(jìn)行了數(shù)值模擬和分析。雙錐旋成體模型的算例表明網(wǎng)格對齊方式對聲爆近場過壓分布具有非常明顯的影響,即網(wǎng)格面沿馬赫角對齊可以顯著提高近場過壓預(yù)測結(jié)果捕捉的精度。DWB模型的算例表明,模型附近采用圓柱區(qū)域包裹的結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格均能很好地捕捉聲爆近場過壓預(yù)測結(jié)果,但非結(jié)構(gòu)網(wǎng)格由于可以快速生成,應(yīng)用更加方便;馬赫錐面方向采用四面體和六面體網(wǎng)格對近場聲爆預(yù)測精度影響不大;在采用帶圓柱區(qū)域的混合網(wǎng)格時,網(wǎng)格展弦比和馬赫錐面周向網(wǎng)格間距對流場中激波峰值和激波形狀的捕捉有明顯影響,使用相對合理的網(wǎng)格展弦比和馬赫錐面周向網(wǎng)格間距有利于提高近場聲爆預(yù)測精度。最后,將總結(jié)的網(wǎng)格生成策略應(yīng)用于LM1021全機構(gòu)型1.42倍機體長度下方近場過壓預(yù)測,計算結(jié)果與實驗測量結(jié)果吻合良好,驗證了其實用性。