羅 錕，汪振國

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

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

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

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

1 CA理論及其聲學模擬

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

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

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

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

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

(1)

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

2 一維聲學CA模型

2.1 局部演化規(guī)則

聲傳播的一維聲波方程為

(2)

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

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

(3)

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

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

(4)

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

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時刻的聲壓值與其左、右兩相鄰位置x+Δx，x-Δx元胞t時刻聲壓和該元胞t-Δt時刻聲壓相關，因此鄰居類型為R=1型，如圖3所示。同時，因為CA模型中聲傳播速度ca等于聲速c，所以局部演化系數(shù)

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

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

(6)

那么式(5)可進一步簡化為

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

(7)

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

2.2 管道聲學CA模擬

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

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

在邊界條件①的聲管中，選擇x=1 m處為測點，入射聲波p在t1時刻到達該位置，則該測點聲壓解析解的表達式為

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

(8)

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

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

(9)

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

圖5 聲管測點聲壓時程曲線Fig.5 Sound pressure time history curve of measuring points in acoustic tube

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

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

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

f=mc/2l

(10)

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

從圖6可知：聲學CA模型能夠模擬聲共振現(xiàn)象，當在自然頻率附近激發(fā)正弦波時，合成聲場會出現(xiàn)駐波，圖5(b)中合成波幅值不等于10 Pa，這是因為入射波頻率為2 000 Hz，介于共振頻率第23階(1 972.25 Hz)～第24階(2 058.25 Hz)。

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

3 二維聲學CA模型

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

聲傳播的二維聲波方程為

(11)

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

定義元胞單元尺寸為Δx=Δy=Δl，時間步為Δt，局部演化系數(shù)φ=(c·Δt/Δl)2，在任意t時刻，式(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時刻的聲壓值與其上、下、左、右4處相鄰位置(x，y+Δy)，(x，y-Δy)，(x-Δx，y)，(x+Δx，y)的元胞t時刻聲壓和該元胞t-Δt時刻聲壓相關，因此鄰居類型為Von Neumann型，如圖7所示。

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

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

(13)

將式(13)代入式(12)，可得Von Neumann型二維聲學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求解點聲源自由空間的輻射聲場，聲場中過點源的某一切片結(jié)果，如圖8所示。由圖8可知：在聲源附近和離聲源較遠位置處，解析解和CA解相差不大，但是在離聲源一定距離的位置上，解析解與CA解相差很大，最大誤差在50%以上，這是由CA模型各向異性的更新規(guī)則所致。

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

3.2 脈動球源輻射聲場

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

(15)

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

(16)

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

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

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

(17)

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

θ=arctan(1/kr0)

(18)

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

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

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

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

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

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

(19)

其中，

(20)

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

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

由圖12與圖13可知：解析解與CA解吻合良好，表明本文所建立的二維CA模型能準確求解脈動球源半自由空間聲輻射問題；同時也說明了對于雙球源組合聲場，該模型依舊具有較高的計算精度和良好的適用性。

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

4 結(jié) 論

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

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

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

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

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