羅 錕，汪振國(guó)

(1.華東交通大學(xué) 鐵路環(huán)境振動(dòng)與噪聲教育部工程研究中心，南昌 330013；2.大連理工大學(xué) 土木工程學(xué)院，遼寧 大連 116024)

在建筑結(jié)構(gòu)、交通運(yùn)輸、工程機(jī)械等諸多領(lǐng)域，噪聲污染問(wèn)題一直受到廣泛關(guān)注[1]，研究準(zhǔn)確有效的聲學(xué)數(shù)值計(jì)算方法是噪聲預(yù)測(cè)和控制技術(shù)發(fā)展的基礎(chǔ)，也是解決噪聲污染問(wèn)題的關(guān)鍵步驟。當(dāng)前，已有許多成熟的數(shù)值計(jì)算方法在聲學(xué)領(lǐng)域得到應(yīng)用，例如：有限元法[2-3]、邊界元法[4-5]、統(tǒng)計(jì)能量法[6-7]、多方法聯(lián)合求解[8-10]等。然而，適于描述復(fù)雜動(dòng)力系統(tǒng)的元胞自動(dòng)機(jī)(cellular automata，CA)方法在聲學(xué)系統(tǒng)的建模和計(jì)算方面卻鮮有應(yīng)用。

CA自其誕生至今，已衍生出幾種獨(dú)立在CA模型之外的計(jì)算模型，如格子氣自動(dòng)機(jī)(lattice gas automata，LGA)模型、格子Boltzmann模型、多粒子模型等，且在各種物理問(wèn)題的建模上均有應(yīng)用[11]。在模擬波傳播方面：Leamy[12]介紹了一種在理想化線彈性介質(zhì)中利用Moore型鄰居CA方法模擬地震波的傳播行為；Nishawala等[13]結(jié)合毗域動(dòng)力學(xué)與CA方法對(duì)彈性波傳播進(jìn)行模擬，結(jié)果顯示CA計(jì)算值與解析解吻合較好；Krutar等[14-15]利用LGA模型模擬了波浪現(xiàn)象，并探討了復(fù)雜場(chǎng)景下水中聲傳播和散射問(wèn)題；Sudo等[16-17]在此基礎(chǔ)上，進(jìn)一步改進(jìn)LGA模型，提出聲傳播的一維和二維LGA模型，其中二維聲傳播模型雖求解穩(wěn)定，但計(jì)算結(jié)果存在各向異性色散；Komatsuzaki等[18]利用CA模型模擬了一維和二維聲傳播現(xiàn)象，其中二維CA模型采用Von Neumann型鄰居，這導(dǎo)致元胞狀態(tài)變量在每一計(jì)算步中出現(xiàn)各向異性更新的問(wèn)題，從而二維CA計(jì)算結(jié)果與解析解間存在較大差異?？梢园l(fā)現(xiàn)，若要應(yīng)用CA方法對(duì)聲學(xué)問(wèn)題進(jìn)行準(zhǔn)確求解，則必須解決元胞各向異性更新的問(wèn)題。

本文首先應(yīng)用CA原理，結(jié)合聲波方程導(dǎo)出一維聲學(xué)CA模型的局部演化規(guī)則，并對(duì)聲管聲場(chǎng)進(jìn)行建模，同時(shí)考慮兩種邊界條件，利用一維CA模型計(jì)算管內(nèi)聲壓分布，并將結(jié)果與解析解進(jìn)行對(duì)比驗(yàn)證；之后，以脈動(dòng)球源聲場(chǎng)為分析對(duì)象，按球坐標(biāo)下的二維聲波方程導(dǎo)出Y函數(shù)一維CA模型的局部演化規(guī)則，進(jìn)而利用Y值與聲壓值間的關(guān)系建立二維脈動(dòng)球源聲場(chǎng)的CA模型，并計(jì)算球源聲輻射性質(zhì)；最后對(duì)比分析二維CA結(jié)果與解析解。

1 CA理論及其聲學(xué)模擬

CA是一種離散化的動(dòng)力系統(tǒng)，其由元胞、元胞空間、鄰居、狀態(tài)變量、局部演化規(guī)則和邊界條件等部分組成[19]，如圖1所示。

圖1 元胞自動(dòng)機(jī)系統(tǒng)Fig.1 Cellular automata system

對(duì)于某一具體問(wèn)題，應(yīng)用CA方法往往需要將計(jì)算區(qū)域在時(shí)間和空間上進(jìn)行離散，每一離散單元稱為元胞或元胞單元，所有元胞組成的集合就是元胞空間；狀態(tài)變量通常為需求解的數(shù)或數(shù)集，每一個(gè)元胞內(nèi)都對(duì)應(yīng)著狀態(tài)變量，且在不同時(shí)刻下實(shí)時(shí)更新；對(duì)于一個(gè)元胞來(lái)說(shuō)，其周圍相鄰的且在單位時(shí)間步內(nèi)參與演化的元胞為其鄰居，如圖2所示。一維CA系統(tǒng)鄰居通常以半徑R為認(rèn)定準(zhǔn)則，即在半徑R以內(nèi)的元胞均為該元胞的鄰居，二維CA系統(tǒng)鄰居類型多樣，其中Von Neumann型與Moore型最為常見；CA系統(tǒng)的邊界條件按所求解問(wèn)題的實(shí)際邊界條件進(jìn)行模擬。

圖2 CA系統(tǒng)鄰居類型Fig.2 Neighbor types of CA system

局部演化規(guī)則是CA系統(tǒng)中最為重要的組成部分，它能夠由當(dāng)前時(shí)刻元胞及其鄰居的狀態(tài)變量確定出下一時(shí)刻元胞的狀態(tài)變量，簡(jiǎn)單來(lái)講，局部演化規(guī)則就是一組狀態(tài)轉(zhuǎn)移函數(shù)，其表達(dá)式為

(1)

應(yīng)用CA方法求解聲學(xué)問(wèn)題的一般步驟是：首先需分析聲源及其輻射聲場(chǎng)特性，明確聲學(xué)CA模型的維度和鄰居類型；之后對(duì)計(jì)算區(qū)域進(jìn)行元胞單元?jiǎng)澐?，并確定元胞尺寸和時(shí)間步長(zhǎng)，同時(shí)以聲壓為狀態(tài)變量，考慮輻射聲場(chǎng)的聲學(xué)邊界條件，并按聲波方程導(dǎo)出局部演化規(guī)則；最后構(gòu)建聲學(xué)CA模型進(jìn)行求解。

2 一維聲學(xué)CA模型

2.1 局部演化規(guī)則

聲傳播的一維聲波方程為

(2)

式中：p為x，t的聲壓函數(shù)；c為聲速。

定義元胞單元尺寸為Δx，時(shí)間步為Δt，在任意t時(shí)刻，式(2)的差分形式為

(3)

定義局部演化系數(shù)

φ=(c·Δt/Δx)2

(4)

將式(4)代入式(3)并化簡(jiǎn)可得

p(x,t+Δt)=φ{(diào)p(x+Δx,t)+p(x-Δx,t)+[(2/φ)-2]p(x,t)-(1/φ)p(x,t-Δt)}

(5)

從式(5)可以看出：x位置處元胞t+Δt時(shí)刻的聲壓值與其左、右兩相鄰位置x+Δx，x-Δx元胞t時(shí)刻聲壓和該元胞t-Δt時(shí)刻聲壓相關(guān)，因此鄰居類型為R=1型，如圖3所示。同時(shí)，因?yàn)镃A模型中聲傳播速度ca等于聲速c，所以局部演化系數(shù)

圖3 一維聲學(xué)CA模型鄰居Fig.3 Neighbors of 1D acoustic CA model

φ=(c·Δt/Δx)2=(c/ca)2=1

(6)

那么式(5)可進(jìn)一步簡(jiǎn)化為

p(x,t+Δt)=p(x+Δx,t)+p(x-Δx,t)-p(x,t-Δt)

(7)

式(7)即為一維聲學(xué)CA模型的局部演化規(guī)則。

2.2 管道聲學(xué)CA模擬

設(shè)有l(wèi)=2 m長(zhǎng)聲管，管體為剛性壁，其法向振速為0，入射壓力波與管道平行，管內(nèi)聲波以平面波的形式傳播，在傳播過(guò)程中不發(fā)生徑向粒子的振動(dòng)，所以管內(nèi)聲場(chǎng)可由一維CA模型模擬，如圖4所示。通過(guò)考慮兩種邊界條件來(lái)確定管中聲壓分布規(guī)律：①左端激勵(lì)，右端開口；②左端激勵(lì)，右端剛壁。

圖4 聲管邊界條件Fig.4 Boundary conditions of acoustic tube

在邊界條件①的聲管中，選擇x=1 m處為測(cè)點(diǎn)，入射聲波p在t1時(shí)刻到達(dá)該位置，則該測(cè)點(diǎn)聲壓解析解的表達(dá)式為

p(x,t)=p0ei(ωt-kx+θ)

(8)

式中：p0為聲壓幅值；ω為圓頻率；k為波數(shù)；θ為位相。

在邊界條件②的聲管中，選擇x=1.5 m處為測(cè)點(diǎn)，入射聲波p在t1時(shí)刻到達(dá)該位置，在t2時(shí)刻到達(dá)末端剛壁并發(fā)生反射，形成沿-x向的反射聲波p-，p-在t3時(shí)刻到達(dá)測(cè)點(diǎn)處，此時(shí)開始，測(cè)點(diǎn)將處于p與p-的合成聲場(chǎng)中。該測(cè)點(diǎn)聲壓解析解的表達(dá)式為

(9)

將管內(nèi)聲場(chǎng)沿x向離散成單元尺寸Δx=1 mm的元胞，以式(7)為局部演化規(guī)則，建立一維聲學(xué)CA模型，對(duì)兩種邊界條件下管內(nèi)聲場(chǎng)聲壓進(jìn)行求解，將計(jì)算結(jié)果與解析解對(duì)比，如圖5所示。計(jì)算參數(shù)如表1所示。

圖5 聲管測(cè)點(diǎn)聲壓時(shí)程曲線Fig.5 Sound pressure time history curve of measuring points in acoustic tube

表1 計(jì)算參數(shù)Tab.1 Calculation parameters

由圖5可知：兩種邊界條件下的CA模型計(jì)算結(jié)果與解析解吻合良好，單向平面波場(chǎng)最大誤差僅為0.614%，合成平面波場(chǎng)最大峰值誤差為1.048%，這表明聲學(xué)CA模型能夠準(zhǔn)確模擬平面波場(chǎng)聲壓分布，且計(jì)算結(jié)果可靠。

在邊界條件②下，聲管的共振頻率f為

f=mc/2l

(10)

式中：m為階數(shù)，m=1,2,…；c為聲速；l為管長(zhǎng)；當(dāng)入射聲波頻率為前3階共振頻率時(shí)，管內(nèi)聲壓分布，如圖6所示。

從圖6可知：聲學(xué)CA模型能夠模擬聲共振現(xiàn)象，當(dāng)在自然頻率附近激發(fā)正弦波時(shí)，合成聲場(chǎng)會(huì)出現(xiàn)駐波，圖5(b)中合成波幅值不等于10 Pa，這是因?yàn)槿肷洳l率為2 000 Hz，介于共振頻率第23階(1 972.25 Hz)～第24階(2 058.25 Hz)。

圖6 聲模態(tài)Fig.6 Acoustic modes

3 二維聲學(xué)CA模型

3.1 Von Neumann型CA局部演化規(guī)則

聲傳播的二維聲波方程為

(11)

式中：p為x，y，t的聲壓函數(shù)；c為聲速。

定義元胞單元尺寸為Δx=Δy=Δl，時(shí)間步為Δt，局部演化系數(shù)φ=(c·Δt/Δl)2，在任意t時(shí)刻，式(11)的差分形式可變化為

p(x,y,t+Δt)=φ{(diào)p(x+Δx,y,t)+p(x-Δx,y,t)+
p(x,y+Δy,t)+p(x,y-Δy,t)+

2[(1/φ)-2]p(x,y,t)-(1/φ)p(x,y,t-Δt)}

(12)

從式(12)可知：(x，y)位置處元胞t+Δt時(shí)刻的聲壓值與其上、下、左、右4處相鄰位置(x，y+Δy)，(x，y-Δy)，(x-Δx，y)，(x+Δx，y)的元胞t時(shí)刻聲壓和該元胞t-Δt時(shí)刻聲壓相關(guān)，因此鄰居類型為Von Neumann型，如圖7所示。

圖7 二維聲學(xué)CA模型鄰居Fig.7 Neighbors of 2D acoustic CA model

φ=(c·Δt/Δl)2=(c/ca)2=1/2

(13)

將式(13)代入式(12)，可得Von Neumann型二維聲學(xué)CA模型的局部演化規(guī)則

p(x,y,t+Δt)=0.5{p(x+Δx,y,t)+
p(x-Δx,y,t)+p(x,y+Δy,t)+p(x,y-Δy,t)}-p(x,y,t-Δt)

(14)

Komatsuzaki等的研究曾利用Von Neumann型CA求解點(diǎn)聲源自由空間的輻射聲場(chǎng)，聲場(chǎng)中過(guò)點(diǎn)源的某一切片結(jié)果，如圖8所示。由圖8可知：在聲源附近和離聲源較遠(yuǎn)位置處，解析解和CA解相差不大，但是在離聲源一定距離的位置上，解析解與CA解相差很大，最大誤差在50%以上，這是由CA模型各向異性的更新規(guī)則所致。

圖8 Komatsuzaki等的研究中點(diǎn)源輻射聲場(chǎng)聲壓分布Fig.8 Distribution of sound pressure in the radiated sound field of point source in Komatsuzaki et al research

3.2 脈動(dòng)球源輻射聲場(chǎng)

為解決二維聲學(xué)CA模型中狀態(tài)變量更新的各向異性問(wèn)題，在導(dǎo)出局部演化規(guī)則時(shí)可將直角坐標(biāo)變換為球坐標(biāo)，以脈動(dòng)球源為例，坐標(biāo)原點(diǎn)取在球心，其球坐標(biāo)下的聲傳播二維聲波方程為

(15)

式中：r為空間任一點(diǎn)至球心的距離；p為r，t的聲壓函數(shù)；S為r處波陣面面積；c為聲速。令Y=pr，則式(15)變換為

(16)

可以發(fā)現(xiàn)，式(16)與一維聲波方程類似，那么應(yīng)用CA方法求解脈動(dòng)球源二維聲場(chǎng)聲壓時(shí)，可先建立Y函數(shù)的一維CA模型，在解出二維聲場(chǎng)中任意場(chǎng)點(diǎn)的Y值之后，只需將Y值比上該場(chǎng)點(diǎn)至球心的距離r即可得到二維聲場(chǎng)中任意場(chǎng)點(diǎn)的聲壓值p。

算例1假設(shè)在二維自由空間中存在一個(gè)脈動(dòng)球源，以球心為中心，取10 m×10 m的矩形區(qū)域?yàn)橛?jì)算空間，將計(jì)算空間以邊長(zhǎng)為0.1 m的正方形網(wǎng)格劃分為100×100個(gè)像元，計(jì)算參數(shù)如表2所示。自由空間脈動(dòng)球源聲場(chǎng)聲壓解析解表達(dá)式為

表2 計(jì)算參數(shù)Tab.2 Calculation parameters

(17)

式中：ρ為介質(zhì)密度；k為波數(shù)；r0為球源半徑；u為球源振速幅值；θ按式(18)計(jì)算

θ=arctan(1/kr0)

(18)

將計(jì)算區(qū)域聲壓分布解析解與球坐標(biāo)下的CA解對(duì)比，如圖9所示。二維自由空間計(jì)算區(qū)域內(nèi)豎向和對(duì)角切片的聲壓分布，如圖10所示。

圖9 自由空間脈動(dòng)球源輻射聲壓云圖Fig.9 Nephograms of sound pressure radiated by pulsating sphere source in free space

圖10 豎向與對(duì)角方向聲壓分布Fig.10 Distribution of sound pressure in vertical and diagonal directions

由圖9和圖10可知：CA結(jié)果顯示，二維脈動(dòng)球源在自由空間內(nèi)的聲輻射形狀為圓周，且聲壓隨距離r的增大而減小，近場(chǎng)聲壓衰減快，遠(yuǎn)場(chǎng)聲壓衰減慢，這與球源輻射聲場(chǎng)性質(zhì)一致；解析解與CA解吻合良好，表明本文所建立的二維CA模型能準(zhǔn)確求解脈動(dòng)球源自由空間聲輻射問(wèn)題；采用球坐標(biāo)聲波方程導(dǎo)出Y函數(shù)一維CA模型局部演化規(guī)則，進(jìn)而求解二維聲場(chǎng)聲壓的方法，可避免二維CA模型中元胞狀態(tài)變量各向異性更新的問(wèn)題。

算例2假設(shè)在二維半自由空間中存在一個(gè)脈動(dòng)球源1，距離球源下方5 m處有一無(wú)限大剛性壁面，以球心為中心，取10 m×10 m的矩形區(qū)域?yàn)橛?jì)算空間，將計(jì)算空間以邊長(zhǎng)為0.1 m的正方形網(wǎng)格劃分為100×100個(gè)像元，如圖11所示。計(jì)算參數(shù)如表2所示。由鏡像原理可知，半自由空間脈動(dòng)球源1聲場(chǎng)等效于自由空間脈動(dòng)球源1及其沿剛壁鏡像球源2聲場(chǎng)的疊加，此時(shí)聲場(chǎng)聲壓解析解表達(dá)式為

圖11 半自由空間脈動(dòng)球源輻射聲場(chǎng)Fig.11 Sound field radiated by pulsating sphere source in half free space

(19)

其中，

(20)

將計(jì)算區(qū)域聲壓分布解析解與CA解對(duì)比，如圖12所示。二維半自由空間計(jì)算區(qū)域內(nèi)橫向、豎向和對(duì)角切片的聲壓分布，如圖13所示。

圖12 半自由空間脈動(dòng)球源輻射聲壓云圖Fig.12 Nephograms of sound pressure radiated by pulsating sphere source in half free space

由圖12與圖13可知：解析解與CA解吻合良好，表明本文所建立的二維CA模型能準(zhǔn)確求解脈動(dòng)球源半自由空間聲輻射問(wèn)題；同時(shí)也說(shuō)明了對(duì)于雙球源組合聲場(chǎng)，該模型依舊具有較高的計(jì)算精度和良好的適用性。

圖13 橫向、豎向與對(duì)角方向聲壓分布Fig.13 Distribution of sound pressure in transverse,vertical and diagonal directions

4 結(jié) 論

本文利用CA方法對(duì)聲管和脈動(dòng)球源聲學(xué)問(wèn)題進(jìn)行建模，分析了一維平面波與二維球面波聲場(chǎng)規(guī)律，并將聲學(xué)CA模型計(jì)算結(jié)果與解析解進(jìn)行對(duì)比，得到以下幾點(diǎn)結(jié)論：

(1)本文建立的聲學(xué)CA模型能準(zhǔn)確模擬平面波與球面波場(chǎng)聲輻射規(guī)律，且結(jié)果與聲波方程的解吻合程度高。

(2)一維CA模型的局部演化規(guī)則可直接由聲波方程導(dǎo)出；而直接按聲波方程導(dǎo)出二維CA模型的局部演化規(guī)則會(huì)使元胞狀態(tài)變量更新存在各向異性的問(wèn)題。

(3)采用球坐標(biāo)聲波方程導(dǎo)出脈動(dòng)球源Y函數(shù)一維CA模型的局部演化規(guī)則，進(jìn)而求解二維聲場(chǎng)聲壓的方法，可避免二維CA模型中元胞狀態(tài)變量各向異性更新的問(wèn)題。

(4)對(duì)于雙球源組合聲場(chǎng)，聲學(xué)CA模型依舊具有較高的計(jì)算精度和良好的適用性；那么對(duì)于多點(diǎn)源組合成的面聲源輻射聲場(chǎng)，利用CA方法分析其聲場(chǎng)性質(zhì)是否可行，這將在下一步工作中繼續(xù)探討。