網格對聲爆近場預測影響的數值研究

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

(西北工業大學 航空學院, 陜西 西安 710072)

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

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

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

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

從現有研究結果來看,目前國外學者針對網格量和網格拓撲結構對近場波形的影響開展了一系列的研究工作[20],但對結構和非結構網格、不同網格類型和網格分布等因素對于近場聲爆預測結果的影響缺少系統研究。為此,本文選取旋成體雙錐和帶翼錐體構型采用數值方法分別就網格對齊方式、網格類型、網格生成方法和分布對聲爆近場預測的精度的影響進行了分析和研究。最后,使用實驗數據齊全的LM1021標模驗證了所得到的結論。

1 數值方法

本文使用自研非結構混合網格復雜流場模擬程序HUNS3D[21]數值求解三維Euler/雷諾平均Navier-Stokes方程來實現近場聲爆的數值預測。

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

(1)

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

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

2 網格生成方法及其影響

2.1 網格對齊方式的影響

聲爆是飛行器在超聲速飛行過程中產生的特定聲學現象,其本質是超聲速飛行時飛行器周圍產生的激波,當其傳播至地面時就形成了聲爆。因此,要準確的對聲爆進行預測,關鍵在于精確的捕捉流場中產生的激波信息。對于定常超聲速流動,激波角μ計算公式為

μ=sin-1(1/M)

(2)

式中,M為來流馬赫數。與傳統正交網格相比,馬赫錐構型的網格計算時相對前者更能夠減小計算耗散,因而可以更為精準的捕捉激波信息,進而能夠更為精準的預測聲爆。

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

(3)

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

圖1雙錐模型示意圖 圖2 雙錐模型網格拓撲示意圖 圖3 雙錐模型網格示意圖

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

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

2.2 結構和非結構網格的影響

網格是進行CFD流場計算的基礎,合理的網格分布是正確模擬一些復雜流動的基本前提。網格主要分為結構網格和非結構網格。結構網格雖然具有可以控制網格拓撲的劃分和網格的正交性的優點,但是對于復雜構型來說,網格生成周期長甚至無法生成。而非結構網格則具有很強的描述復雜幾何形狀構型的能力,易于得到高質量的網格,同時,非結構網格生成過程快速且自動化程度高。

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

圖5 DWB模型、2種計算域及近壁面網格示意圖

圖5a)～c)所示為DWB模型、計算域示意圖。根據拓撲結構可將網格分為模型附近帶圓柱區域和不帶圓柱區域。所謂圓柱區域,即為網格拓撲中模型附近包裹模型的部分。對于帶圓柱區域網格,為減小激波信號傳播過程中的耗散,圓柱應盡可能小地包圍模型。圓柱區域以外采用馬赫錐構型,即網格沿馬赫角由圓柱表面往外拉伸,本算例中拉伸部分網格為六面體單元。圖5d)～f)分別為采用非結構混合網格近場、帶圓柱內場結構網格近場和不帶圓柱結構網格近場示意圖。網格數據如表1所示。

表1 不同網格類型網格量示意圖

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

圖6 使用結構和非結構網格計算的DWB近場過壓預測結果與實驗結果的比較

除此之外,由圖6a)～d)可以明顯看出隨著方位角的變化,尾部激波從單個激波演變為多個激波,數值模擬捕捉到的演變程與實驗曲線表現的變化規律基本一致。

2.3 網格類型的影響

為研究圓柱區域外網格類型對聲爆波形傳播的影響,本文在DWB楔形拓撲結構基礎上在圓柱區域以外分別使用了四面體、六面體混合網格和全四面體網格,如圖7a)～b)所示。其中,四面體、六面體混合網格網格節點數為5 303 747,單元數為8 734 621,四面體混合網格網格節點數為3 172 122,單元數為9 755 505。

圖7 DWB模型不同網格類型示意圖

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

2.4 網格分布的影響

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

圖9 網格分布參數定義

表2 不同展弦比網格信息

2.4.1 網格展弦比的影響

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

圖10 DWB模型不同展弦比網格計算的近場聲爆過壓預測結果與實驗結果對比

2.4.2 周向網格間距的影響

表3給出了不同周向間距的3套測試網格詳細信息。

表3 不同周向距離網格信息

圖11 DWB模型不同周向距離網格計算的近場聲爆過壓預測結果與實驗結果對比

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

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

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

3 復雜構型近場聲爆預測

3.1 模型和網格生成

為驗證前文所述的網格生成策略在復雜構型中的應用效果,本文選取了SBPW-1中作為可選算例的標模LM1021全機方案進行測試。采用帶圓柱區域的混合網格拓撲結構,其中圓柱內周向角度為2,遠場網格展弦比為1。圓柱區域以外采用沿馬赫角對齊六面體網格以便于捕捉聲爆傳播。

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

3.2 計算分析

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

圖12 LM1021構型及網格示意圖

圖13 LM1021構型CFD計算的近場聲爆過壓預測結果與實驗結果對比

4 結 論

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