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

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

?

計(jì)算氣動(dòng)聲學(xué)中的伽遼金玻爾茲曼方法研究

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

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

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

計(jì)算氣動(dòng)聲學(xué);伽遼金玻爾茲曼方法;無反射邊界條件

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

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

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

1 伽遼金玻爾茲曼方法

1.1 無黏速度空間BGK方程

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

α=0,1,…,Nα

(1)

式中:fα為沿格子速度eα=(eαx,eαy)的粒子分布函數(shù);λ為弛豫時(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。

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

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

(4)

對(duì)式(4)右端進(jìn)行積分有

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

(5)

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

(6)

于是式(5)可重寫為

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

(7)

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

碰撞步

(8)

對(duì)流步

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

(9)

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

(10)

(11)

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

(12)

1.2 時(shí)空離散格式

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

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

(13)

(14)

對(duì)式(14)中通量進(jìn)行分步積分得

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

(15)

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

(16)

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

(17)

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

(18)

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

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

(19)

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

(20)

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

(21)

(22)

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

(23)

式(23)中的導(dǎo)數(shù)項(xiàng)有

(24)

(25)

式(24)、(25)中幾何項(xiàng)可逐點(diǎn)計(jì)算

(26)

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

(27)

式(27)中的每一項(xiàng)有

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

(28)

(29)

(30)

(31)

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

(32)

顯式五級(jí)四階龍格庫塔法[14]的穩(wěn)定區(qū)域和存儲(chǔ)空間都優(yōu)于經(jīng)典龍格庫塔法,對(duì)式(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為文獻(xiàn)優(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)

對(duì)于聲傳播的硬壁面,物理量滿足

(36)

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

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

(37)

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

根據(jù)流體黏性對(duì)聲的吸收和衰減原理,構(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方向?yàn)橹芷谛赃吔鐥l件,這里取ρ0=1.0,Δρ=0.017 321,Λ=50,υ=10υair=0.044,υair表示空氣黏性系數(shù)。采用伽遼金玻爾茲曼方法求解上述算例,圖2給出了聲粒子速度在t=3 000時(shí)模擬值和解析解的比較。二者吻合較好,說明了該方法可以模擬聲的吸收和衰減作用。

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

圖3 運(yùn)動(dòng)黏性系數(shù)對(duì)聲粒子速度的影響

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

圖4 波長對(duì)聲粒子速度的影響

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

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

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

(39)

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

(40)

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

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

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

(41)

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

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

(42)

(a)脈動(dòng)聲源傳播的計(jì)算域

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

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

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

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

圖7 4個(gè)不同時(shí)刻的聲壓分布

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

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

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

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

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

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

3 結(jié) 論

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

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

[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方程的氣動(dòng)聲學(xué)計(jì)算研究 [J]. 航空學(xué)報(bào), 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ì)算氣動(dòng)聲學(xué)進(jìn)展與展望 [J]. 中國科學(xué): 物理學(xué)力學(xué)天文學(xué), 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ì)算氣動(dòng)聲學(xué)中的高精度緊致差分格式研究 [J]. 航空動(dòng)力學(xué)報(bào), 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)用于氣動(dòng)聲學(xué)研究 [J]. 計(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—),男,博士生;李軍(通信作者),男,教授,博士生導(dǎo)師。 基金項(xiàng)目:國家自然科學(xué)基金資助項(xiàng)目(51376144)。

時(shí)間: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