邵衛(wèi)東,李軍,2

(1.西安交通大學能源與動力工程學院,710049,西安;2.先進航空發(fā)動機協(xié)同創(chuàng)新中心,100191,北京)

?

計算氣動聲學中的伽遼金玻爾茲曼方法研究

邵衛(wèi)東1,李軍1,2

(1.西安交通大學能源與動力工程學院,710049,西安;2.先進航空發(fā)動機協(xié)同創(chuàng)新中心,100191,北京)

為獲得氣動聲學的高精度和低耗散特性的數(shù)值方法,發(fā)展了伽遼金玻爾茲曼方法和相應(yīng)的無反射邊界條件。首先,引入新粒子分布函數(shù)到格子玻爾茲曼BGK方程中并重構(gòu)歐拉方程;然后,在空間上采用高精度的交點間斷伽遼金有限元方法,在時間上采用顯式五級四階龍格庫塔離散方法對解耦得到的對流步方程進行離散求解;最后,通過數(shù)值通量構(gòu)造速度邊界、聲學硬壁面邊界和無反射邊界條件。采用包含聲反射和多普勒效應(yīng)的數(shù)值算例進行驗證,可得模擬值與解析解吻合一致,從而證明了伽遼金玻爾茲曼方法和無反射邊界條件用于氣動聲學計算的有效性和準確性。

計算氣動聲學;伽遼金玻爾茲曼方法;無反射邊界條件

氣動聲學問題一般有兩類數(shù)值解法:基于Lighthill聲類比理論[1]及拓展的聲比擬預測方法;直接模擬研究噪聲的產(chǎn)生機理及傳播特性的計算氣動聲學(CAA)方法。通過計算流體動力學或試驗得到流場信息并結(jié)合聲比擬方法預測噪聲在工程中有廣泛的應(yīng)用[2-3],但聲比擬方法既不能考慮流場與聲場的相互作用,又不能預測非線性氣動噪聲,故要從更本質(zhì)的角度研究流動發(fā)聲機理就必須采用CAA方法。

在CAA中兩大關(guān)鍵技術(shù)是高精度的時空離散格式和無反射邊界條件(NRBC)[4-5]。Tam等提出了色散相關(guān)保持差分格式[4],Hu等提出了低耗散低色散龍格庫塔格式[6],柳占新等發(fā)展了五對角緊致差分格式[7]。這些高階格式能夠顯著減少每個波長所需的網(wǎng)格點數(shù),但處理無反射邊界時仍需重新構(gòu)造且精度降低。與Navier-Stokes方程的高階格式相比,格子玻爾茲曼方法(LBM)[8]具有更低的耗散特性和更快的計算速度,但采用LBM直接模擬[9-10]要有規(guī)則的幾何和強大的計算資源。鄧義求等用LBM研究了點源輻射等聲學現(xiàn)象[11],但與NRBC的結(jié)合并不理想。Min等將間斷伽遼金有限元與LBM結(jié)合起來研究了黏性流動而沒有考慮其在無黏聲傳播中的應(yīng)用[12-13]。

為兼具高精度格式和LBM的優(yōu)點,本文發(fā)展了時空上高精度離散求解無黏速度空間BGK方程的伽遼金玻爾茲曼方法,建立了相應(yīng)的NRBC,并用聲反射和多普勒效應(yīng)的算例進行了驗證。

1 伽遼金玻爾茲曼方法

1.1 無黏速度空間BGK方程

玻爾茲曼方程在離散速度空間下并考慮BGK近似的方程為

α=0,1,…,Nα

(1)

式中:fα為沿格子速度eα=(eαx,eαy)的粒子分布函數(shù);λ為弛豫時間;Nα為離散速度數(shù)。本文采用D2Q9模型為

(2)

平衡態(tài)粒子分布函數(shù)定義為

(3)

式中:ρ、u為宏觀密度和速度;cs=3-1/2為格子聲速;tα為權(quán)系數(shù),t0=4/9,t1,2,3,4=1/9,t5,6,7,8=1/36。

對式(1)沿特征線積分有

fα(x+eαΔt,t+Δt)-fα(x,t)=

(4)

對式(4)右端進行積分有

fα(x+eαΔt,t+Δt)-fα(x,t)=

(5)

式中:τ=λ/Δt為歸一化弛豫時間;?為隱式所占的比例,本文中取為0.5。

(6)

于是式(5)可重寫為

gα(x+eαΔt,t+Δt)-gα(x,t)=

(7)

式(7)的求解又可分為:

碰撞步

(8)

對流步

gα(x+eαΔt,t+Δt)=gα(x,t)

(9)

則宏觀物理量密度ρ、速度u=(u,v)、壓力p可由粒子分布函數(shù)表示為

(10)

(11)

式(9)準確描述了對流過程,但是要求均勻網(wǎng)格和單位CFL數(shù),計算時間長且對不規(guī)則幾何無法求解。為此,求解可替代的純對流方程

(12)

1.2 時空離散格式

交點間斷伽遼金有限元方法具有局部守恒性和選取數(shù)值通量的靈活性,從而獲得在一般非規(guī)則網(wǎng)格上的高精度并保證波占優(yōu)問題的穩(wěn)定性,故對式(12)應(yīng)用此法求解。

引入通量Gα(g)=eαgα,式(12)變?yōu)?/p>

(13)

(14)

對式(14)中通量進行分步積分得

-∫Γkφn·Gα(g)dΓ

(15)

式中:Γk為Ωk的邊界;n=(nx,ny)為外法向量。

(16)

對式(16)再次分步積分,單元弱形式方程為

(17)

為引入迎風特性,對式(17)右端中數(shù)值通量采用Lax-Friedrichs通量,即

(18)

式中:C=max(n·?G/?g)=n·eα。

則式(17)右端中數(shù)值通量為

(19)

在單元Ωk中定義新粒子分布函數(shù)

(20)

采用雙線性四邊形等參單元將單元內(nèi)的點(x,y)∈Ωk映射到參考域內(nèi)的點(ξ,η)∈I=[-1,1]2上。采用高斯GLL積分點的一維拉格朗日插值基函數(shù)為

(21)

(22)

N+1個高斯GLL積分點的集合{ξ0,ξ1,…,ξN}是方程(22)的解,則二維張量基可定義為ψij(ξ,η)=li(ξ(x))lj(η(y))。將式(20)代入式(17),并選取試驗函數(shù)為張量基φ=ψi*j*,則離散弱形式為

(23)

式(23)中的導數(shù)項有

(24)

(25)

式(24)、(25)中幾何項可逐點計算

(26)

對式(23)采用高斯積分法則,可得單元Ωk上的半離散形式

(27)

式(27)中的每一項有

Mk=(ψij,ψi*j*)Ωk=

(28)

(29)

(30)

(31)

對式(27)兩邊左乘質(zhì)量矩陣Mk的逆矩陣,得

(32)

顯式五級四階龍格庫塔法[14]的穩(wěn)定區(qū)域和存儲空間都優(yōu)于經(jīng)典龍格庫塔法,對式(32)應(yīng)用該方法,得

k(i)=aik(i-1)+ΔtL(p(i-1),t+ciΔt),i∈[1,…,5]

p(i)=p(i-1)+bik(i)

(33)

式中:ai、bi、ci為文獻優(yōu)化的常數(shù)。

本文取參考長度Lref、速度Vref和密度ρref為10-5m、343.4 m/s和1.2 kg/m3,文中變量均以參考量作歸一化處理。

1.3 邊界條件的構(gòu)造

速度邊界條件通過數(shù)值通量施加,將式(19)簡化為

(34)

定義交界面鄰近單元的粒子分布函數(shù)為

(35)

對于聲傳播的硬壁面,物理量滿足

(36)

式中:u′=(u′,v′)為聲粒子速度;p′為聲壓。

以圖1為例,鄰近單元的粒子分布函數(shù)為

(37)

注:0～8表示粒子速度的方向圖1 壁面粒子分布函數(shù)示意圖

根據(jù)流體黏性對聲的吸收和衰減原理,構(gòu)造無反射邊界條件?？紤]沿x軸正方向波數(shù)k=2π/Λ和圓頻率ω=csk的平面聲波從源x=0處開始傳播,產(chǎn)生的正弦波的幅值U=csΔρ/ρ0,相應(yīng)的密度ρ=ρ0+Δρsin(ωt)。上述算例的解析解為

u(x,t)=Uexp(-asx)sin(ωt-kx)

(38)

取x方向足夠長和y方向為周期性邊界條件,這里取ρ0=1.0,Δρ=0.017 321,Λ=50,υ=10υair=0.044,υair表示空氣黏性系數(shù)。采用伽遼金玻爾茲曼方法求解上述算例,圖2給出了聲粒子速度在t=3 000時模擬值和解析解的比較。二者吻合較好,說明了該方法可以模擬聲的吸收和衰減作用。

圖2 聲粒子速度的模擬值與解析解的比較

圖3 運動黏性系數(shù)對聲粒子速度的影響

考察3種不同運動黏性系數(shù)25υair、50υair和100υair對聲吸收的影響,圖3給出了t=3 000時的聲粒子速度沿x軸分布。聲波沿著傳播方向衰減,因此其幅值越來越小。當運動黏性系數(shù)增加,聲粒子速度幅值減小且衰減較為明顯。

圖4 波長對聲粒子速度的影響

考察3種不同波長25、50和100在運動黏性系數(shù)υ=25υair時對聲傳播的影響,圖4給出了t=3 000時的聲粒子速度沿x軸正方向分布。聲粒子速度幅值沿著波的傳播方向衰減,而隨著波長的增加衰減變得緩慢。

由傅里葉級數(shù)可知,任意聲波可由許多正弦波疊加而成。對聲波的衰減可轉(zhuǎn)換為對每一個正弦波的衰減,聲粒子速度幅值滿足

|u(x,t)|/U=σ

(39)

式中:σ=1%為給定的衰減程度。由式(38)可得,在給定最大波長下衰減聲波到給定范圍所需的最小距離為

(40)

基于式(40)構(gòu)造無反射邊界條件,為盡可能減小計算域的大小,采用高運動黏性系數(shù)來建立高黏性吸聲區(qū)(HVAZ)。在HVAZ和無黏區(qū)(IZ)之間需要添加過渡吸聲區(qū)(TAZ)來緩和物理量的劇烈變化,從而提高數(shù)值計算的穩(wěn)定性。對于對流步3個區(qū)域求解方程相同,而對于碰撞步IZ、TAZ和HVAZ 3個區(qū)域弛豫時間對應(yīng)的黏性系數(shù)取為0、2υair和100υair。

2 數(shù)值算例與驗證

為驗證上述方法和邊界條件的可行性,本文研究了包含聲反射和多普勒效應(yīng)的算例:初始脈動源在馬赫數(shù)Ma為0.5的亞聲速平均流中與硬壁面相互作用,計算網(wǎng)格如圖5所示。初始t=0脈動源為

(41)

式中:xs=0,ys=25為初始源位置;β=ln(2/25)為源形狀因子。為了長距離傳播而不產(chǎn)生激波,選取源壓力脈動幅值ε=0.01,并給出解析解[15]

[J0(ζξ)+J0(ζη)]ζdζ

(42)

(a)脈動聲源傳播的計算域

(b)脈動聲源傳播的計算網(wǎng)格圖5 脈動聲源傳播的計算域和網(wǎng)格

HVAZ的邊界采用速度邊界條件,壁面采用聲學硬壁面邊界條件,初始粒子分布函數(shù)采用麥克斯韋平衡函數(shù)。為驗證網(wǎng)格的無關(guān)性,選取Ma=0.5并使用12 221和17 061個節(jié)點來對比4個不同時刻聲壓的模擬值和解析解。圖6給出了兩種不同網(wǎng)格數(shù)目計算得到的聲壓與解析解的比較,可知兩種網(wǎng)格數(shù)目獲得的數(shù)值解均與解析解吻合一致。采用12 221個節(jié)點足夠描述該物理現(xiàn)象且該方法能夠捕捉聲波的變化。隨著傳播距離的增大,聲壓的幅值逐漸減小。

圖6 沿著直線x=y,s=(x2+y2)1/2的聲壓分布

圖7給出了4個不同時刻的聲壓分布云圖。t=15時刻之前聲波各向同性地向四周傳播且在t=15時刻幾乎到達壁面;在t=45時刻波前撞擊壁面而形成反射波;隨著時間推進,反射波跟隨著原始波傳播并在壁面附近相互干涉。隨著與壁面距離的增加,聲壓減弱。兩種波在t=75、100時刻都光滑地離開右邊界而沒有受到HVAZ邊界的干擾,證明了NRBC構(gòu)造的有效性。由于平均流作用,向右傳播的波前速度是向左波前速度的3倍。

圖7 4個不同時刻的聲壓分布

圖8給出了4個不同時刻計算聲壓和解析解沿著壁面的分布比較,二者吻合一致,從而證明了聲學硬壁面邊界條件構(gòu)造的準確性。在t=15時刻聲波剛到達壁面,只能在壁面上形成單波峰且聲壓值較小;隨著時間推移,聲波撞擊壁面區(qū)域變大,故能形成雙波峰且聲壓值較大;在t=75、100時刻右波峰已經(jīng)傳出了右邊界。

圖8 聲壓沿著y=0的分布比較

圖9 t=45時刻不同馬赫數(shù)下的聲壓分布

進一步研究了4種不同馬赫數(shù)下聲波與壁面的干涉情況,圖9給出了t=45時刻不同馬赫數(shù)下的聲壓分布。仔細觀察發(fā)現(xiàn),隨著馬赫數(shù)的增大,左波前向左傳播的距離逐漸減小,右波前向右傳播的距離逐漸增大,這一現(xiàn)象也可從聲波與壁面干涉的位置上發(fā)現(xiàn)。盡管聲波的位置不同,但聲壓大小基本一致。

圖10 t=45時刻聲壓沿著y=0的分布

圖10給出了不同馬赫數(shù)下t=45時刻聲壓沿著y=0的分布,隨著馬赫數(shù)的增大,聲壓分布整體向右平移。這是由于多普勒效應(yīng)的作用,聲波右行速度為左行速度的(1+Ma)/(1-Ma)倍,與圖7中Ma=0.5的結(jié)論一致。但是,聲壓分布的形狀并未改變,波峰與波谷之間的距離也沒有變化,說明在相對流體的坐標系中聲波的傳播速度是絕對的,也就是說以聲速傳播。

3 結(jié) 論

本文發(fā)展了可用于氣動聲學計算的伽遼金玻爾茲曼方法。將格子玻爾茲曼BGK方程解耦為碰撞步和對流步,調(diào)節(jié)運動黏性系數(shù)以獲得模擬聲傳播的歐拉方程;對對流步方程應(yīng)用空間上的交點間斷伽遼金有限元方法和時間上的顯式五級四階龍格庫塔方法離散求解;通過數(shù)值通量構(gòu)造相應(yīng)的聲學邊界條件。

本文研究了包含聲反射和多普勒效應(yīng)的算例,數(shù)值解與解析解的良好吻合證明了伽遼金玻爾茲曼方法和邊界條件的有效性和準確性。隨著聲波與源距離的增大,聲壓幅值減小;隨著馬赫數(shù)增大,多普勒效應(yīng)顯著但聲波相對于流體仍以聲速傳播。

[1] LIGHTHILL M J. On sound generated aerodynamically: I General theory [C]∥Proceedings of the Royal Society of London: A Mathematical, Physical and Engineering Sciences. London, UK: The Royal Society, 1952, 211(1107): 564-587.

[2] 余雷, 宋文萍, 韓忠華, 等. 基于混合RANS/LES方法與FW-H方程的氣動聲學計算研究 [J]. 航空學報, 2013, 34(8): 1795-1805. YU Lei, SONG Wenping, HAN Zhonghua, et al. Aeroacoustic noise prediction using hybrid RANS/LES method and FW-H equation [J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(8): 1795-1805.

[3] MAO Y, XU C, QI D. Analytical solution for sound radiated from the rotating point source in uniform subsonic axial flow [J]. Applied Acoustics, 2015, 92: 6-11.

[4] TAM C K W. Computational aeroacoustics: issues and methods [J]. AIAA Journal, 1995, 33(10): 1788-1796.

[5] 李曉東, 旻江, 高軍輝. 計算氣動聲學進展與展望 [J]. 中國科學: 物理學力學天文學, 2014, 44(3): 234-248. LI Xiaodong, MIN Jiang, GAO Junhui. Progress and prospective of computational aeroacoustics [J]. Sci Sin: Phys Mech Astron, 2014, 44(3): 234-248.

[6] HU F Q, HUSSAINI M Y, MANTHEY J L. Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics [J]. Journal of Computational Physics, 1996, 124(1): 177-191.

[7] 柳占新, 黃其柏, 胡溧, 等. 計算氣動聲學中的高精度緊致差分格式研究 [J]. 航空動力學報, 2009, 24(1): 83-90. LIU Zhanxin, HUANG Qibai, HU Li, et al. Study on the high-accuracy compact-finite-difference schemes for computational aeroacoustics [J]. Journal of Aerospace Power, 2009, 24(1): 83-90.

[8] MARIé S, RICOT D, SAGAUT P. Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics [J]. Journal of Computational Physics, 2009, 228(4): 1056-1070.

[9] LI X M, LEUNG R C, SO R M. One-step aeroacoustics simulation using lattice Boltzmann method [J]. AIAA Journal, 2006, 44(1): 78-89.

[10]TSUTAHARA M, KATAOKA T, SHIKATA K, et al. New model and scheme for compressible fluids of the finite difference lattice Boltzmann method and direct simulations of aerodynamic sound [J]. Computers & Fluids, 2008, 37(1): 79-89.

[11]鄧義求, 唐政, 董宇紅. 格子Boltzmann方法應(yīng)用于氣動聲學研究 [J]. 計算物理, 2013, 30(6): 808-814. DENG Yiqiu, TANG Zheng, DONG Yuhong. Lattice Boltzmann method for simulating propagating acoustic waves [J]. Chinese Journal of Computational Physics, 2013, 30(6): 808-814.

[12]MIN M, LEE T. A spectral-element discontinuous Galerkin lattice Boltzmann method for nearly incompressible flows [J]. Journal of Computational Physics, 2011, 230(1): 245-259.

[13]ZADEHGOL A, ASHRAFIZAADEH M, MUSAVI S H. A nodal discontinuous Galerkin lattice Boltzmann method for fluid flow problems [J]. Computers & Fluids, 2014, 105: 58-65.

[14]CARPENTER M H, KENNEDY C A. Fourth-order 2N-storage Runge-Kutta schemes [J/OL]. [2015-06-06]. http: ∥www.ece.uvic.ca/～bctill/papers/numacoust/Carpenter_Kennedy_1994.pdf.

[15]HARDIN J C, RISTORCELLI J R, TAM C K W. ICASE/LARC workshop on benchmark problems in computational aeroacoustics [R]. Washington, DC, USA: NASA, 1995: 10-11.

(編輯 趙煒 苗凌)

Study on the Galerkin Boltzmann Method for Computational Aeroacoustics

SHAO Weidong1,LI Jun1,2

(1. School of Energy & Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China;2. Collaborative Innovation Center of Advanced Aero-Engine, Beijing 100191, China)

To get high-accuracy and low dissipative numerical method in aeroacoustics, Galerkin Boltzmann method and corresponding nonreflecting boundary condition (NRBC) were developed. A modified particle distribution function was introduced to lattice Boltzmann BGK equation in order to reconstruct the Euler equation. After decoupling the collision step from the streaming step, we implemented the high-accuracy nodal discontinuous Galerkin spatial discretization and fourth-order, five-stage Runge-Kutta time marching scheme to solve the resulting advection equation. Velocity boundary condition, acoustically hard wall boundary condition and NRBC were constructed through numerical flux. A benchmark problem including acoustic reflection and Doppler effects was used to test the functionality and accuracy of this method and NRBC. Computational results are in good agreement with the analytical solution, implying that it is a promising method for computational aeroacoustics.

computational aeroacoustics; Galerkin Boltzmann method; nonreflecting boundary condition

10.7652/xjtuxb201603021

2015-08-06。 作者簡介:邵衛(wèi)東(1989—),男,博士生;李軍(通信作者),男,教授,博士生導師。 基金項目:國家自然科學基金資助項目(51376144)。

時間:2015-12-28

http:∥www.cnki.net/kcms/detail/61.1069.T.20151228.1956.006.html

V2111 3

:A

:0253-987X(2016)03-0134-07