王海濤, 曾向陽, 劉延善

(西北工業(yè)大學(xué) 航海學(xué)院, 陜西 西安 710072)

界面聲反射模擬是室內(nèi)復(fù)雜聲學(xué)現(xiàn)象模擬的關(guān)鍵。早期的模型只考慮鏡面反射,近年來,散射聲受到了普遍重視,因為它與室內(nèi)擴散緊密相關(guān),也對混響時間計算、可聽化、聲學(xué)舒適度等諸多因素存在影響。周期結(jié)構(gòu)(periodic structure)是一種特殊的界面散射結(jié)構(gòu),它是指材料的幾何形狀在空間上按照周期規(guī)則排列的結(jié)構(gòu)。周期結(jié)構(gòu)在形式上早已存在且十分常見,在音樂廳、錄音室、電影院等對聲場擴散要求較高的場所經(jīng)?？梢钥吹?與傳統(tǒng)無規(guī)散射體相比,周期結(jié)構(gòu)具有散射性質(zhì)可預(yù)測、可控等特點,雖然有時因其存在也會產(chǎn)生不利音質(zhì)的因素,如聲染色效應(yīng)(coloration effect)[1],但是周期結(jié)構(gòu)仍然以散射性質(zhì)可控、易于設(shè)計等特點而得到廣泛應(yīng)用。

散射系數(shù)是描述周期結(jié)構(gòu)散射性質(zhì)的主要參數(shù),它定義為非鏡面反射能量與全部反射能量之比[2]。獲取散射系數(shù)主要有實驗測量及數(shù)值計算2種方法。2004年ISO組織制定了散射系數(shù)測量標(biāo)準(zhǔn)[3],但需要的硬件條件較高,過程復(fù)雜,且測量精度受環(huán)境因素影響較大,因而近年來數(shù)值計算方法受到的研究更多,且有多種可選的算法,可分為2類:①解析方法,包括kirchhoff approximation(KA)[4]法、holford-urusovskii(HU)[5]法等,它們適用于某些具有特定輪廓的周期結(jié)構(gòu),例如矩形輪廓的周期結(jié)構(gòu)[6],或者周期結(jié)構(gòu)的輪廓可用數(shù)學(xué)函數(shù)精確表達的情況,顯然,對于某些形狀不規(guī)則的周期結(jié)構(gòu),此類方法受到的局限性較高。②以boundary element method (BEM)[7],finite element-plane wave decomposition(FE-PWD)[8]為代表的數(shù)值算法,這類方法結(jié)合插值算法,可對任意形狀的周期結(jié)構(gòu)進行計算,是目前應(yīng)用最多的散射系數(shù)計算方法。

實際中的周期結(jié)構(gòu)通常附著在建筑內(nèi)部的墻壁或頂部,一般尺度較大,而在前述數(shù)值方法中,周期結(jié)構(gòu)的幾何模型需要被離散化以進行插值運算,因而會導(dǎo)致計算量太大而無法對原始尺度進行研究。現(xiàn)有研究通常利用小尺度周期結(jié)構(gòu)(為方便起見,本文中稱為小樣品)的散射系數(shù)來等效大尺度的周期結(jié)構(gòu)。但是由于缺乏不同尺度周期結(jié)構(gòu)散射性質(zhì)之間聯(lián)系的定量研究,使實際選取樣品時顯得無據(jù)可依。目前僅能確定的是圓形周期結(jié)構(gòu)與方形周期結(jié)構(gòu)得到的散射系數(shù)是一致的,這是Vorl?nder等[9]和Lin等[10]在研究邊緣效應(yīng)問題時利用實驗測量得出的結(jié)論,Kosaka等[7]利用數(shù)值計算也給出了相同的結(jié)論。但關(guān)于不同尺度周期結(jié)構(gòu)散射系數(shù)的關(guān)系尚沒有研究報道。

為此,本文首先根據(jù)散射系數(shù)的獲取原理,分析了不同尺度的周期結(jié)構(gòu)之間散射系數(shù)的內(nèi)在聯(lián)系;然后推導(dǎo)了一種無網(wǎng)格周期結(jié)構(gòu)散射系數(shù)計算模型(boundary meshless model, BMM),對代表不同大尺度周期結(jié)構(gòu)局部區(qū)域的小樣品進行了散射系數(shù)數(shù)值計算,通過比較無規(guī)入射散射系數(shù)及方向入射散射系數(shù),確定了可等效大樣品的小樣品需滿足的條件。

1 理論分析

在混響環(huán)境中,根據(jù)散射系數(shù)的定義可得周期結(jié)構(gòu)散射系數(shù)的計算公式:

(1)

式中：αspec為無規(guī)入射鏡面吸聲系數(shù);αs為無規(guī)入射吸聲系數(shù);V為混響室的體積;c表示聲速;Sp為周期結(jié)構(gòu)的表面積;T1為無試件但有底板時得到的混響時間;T2為有試件時得到的混響時間;T3和T4的獲取需要將周期結(jié)構(gòu)放置在轉(zhuǎn)臺上。T3為無試件但底板旋轉(zhuǎn)時得到的混響時間;T4為有試件底板旋轉(zhuǎn)時得到的混響時間。根據(jù)賽賓公式,轉(zhuǎn)臺不轉(zhuǎn)時的混響時間可表示為:

(2)

式中：Seαe表示除去大尺度周期結(jié)構(gòu)所在方形區(qū)域以外的所有吸聲量;Sf表示無試件時,大尺度周期結(jié)構(gòu)所在方形區(qū)域的面積。

如圖1所示,假設(shè)此時存在一個小樣品,其面積為大尺度周期結(jié)構(gòu)的n倍(0

(3)

圖1 測量實驗中不同尺度周期結(jié)構(gòu)及轉(zhuǎn)臺示意圖

(4)

(5)

小樣品的散射系數(shù)可表示為:

(6)

可見,對于不同尺度的周期結(jié)構(gòu),無規(guī)入射吸聲系數(shù)是相同的,計算散射系數(shù)時僅有無規(guī)入射鏡面吸聲系數(shù)不同。

上述推導(dǎo)過程說明,利用小樣品等效大尺度周期結(jié)構(gòu)獲取散射系數(shù)具有理論依據(jù),但到底何種尺度的小樣品具有代表性還需進行定量研究。

2 數(shù)值計算

2.1 邊界無網(wǎng)格型散射系數(shù)計算模型BMM

首先推導(dǎo)一種結(jié)合無網(wǎng)格算法與Mommertz[11]散射系數(shù)計算理論的邊界無網(wǎng)格型散射系數(shù)計算模型BMM,可對任意形狀、尺度的周期結(jié)構(gòu)的散射系數(shù)進行計算。在Mommertz的理論中,一個與散射體尺度、材料相同的純平參考板被引入到散射系數(shù)的計算之中。假設(shè)在散射體及參考板的遠場范圍內(nèi)有一半球面,上面布置著一系列均勻分布的接收點。對于某一聲源來說,只要求得散射體及參考板在接收點的聲壓,就可以根據(jù)(7)式計算出此聲源入射方向上的方向散射系數(shù)(directional scattering coefficient)

δ(θ,φ)=

(7)

式中:θ和φ分別代表聲源的俯仰角與方位角,每個聲源都對應(yīng)于一個入射方向;n為接收點的數(shù)量;θ′和φ′分別代表第i個接收點的俯仰角與方位角;p1是散射體所對應(yīng)的接收點的聲壓;p0是參考板所對應(yīng)的接收點的聲壓;*代表復(fù)共軛。

如果同時假設(shè)在散射體及參考板的遠場范圍內(nèi)有一個與接收點半球面類似的半球面,上面分布著若干聲源,那么只要求得所有這些聲源的方向散射系數(shù),就可以根據(jù)Paris′ formula[12]求得無規(guī)入射散射系數(shù)(random-incidence scattering coefficient):

(8)

因此,散射體及參考板接收點聲壓的計算是計算散射系數(shù)的關(guān)鍵。

假設(shè)聲源在單位時間內(nèi)向單位體積的空間提供了ρ0q(q為聲源的體積速度)的媒質(zhì)質(zhì)量,則此時聲場的有源Helmholtz方程為:

2pω+k2pω=-jρ0ωqω

(9)

式中：k=ω/c,為波數(shù)。

在三維問題中,Helmholtz方程存在基本解G(P,Q)=e-jkr/4πr,其中P、Q為聲場中的任意2點,r表示2點之間的距離。它表示當(dāng)聲場中某點存在單位強度的“源”時,對另外一點所產(chǎn)生的影響,結(jié)合格林第二公式可得:

(10)

假設(shè)聲源是位于r0(x0,y0,z0)處的點聲源,則(10)式可寫為:

(11)

由于采用邊界積分形式的數(shù)學(xué)推導(dǎo),(11)式在某些頻率處無法求得唯一解,這些相應(yīng)的頻率被稱為“偽頻率”。為了克服非唯一解所帶來的問題,可以首先對(11)式求一次偏導(dǎo)[13],并應(yīng)用邊界條件?p/?n=0,然后將得到的式子與(11)式進行線性組合得到在任意頻率下都可以求得唯一解的方程。

(12)

式中：β為非零耦合常數(shù),一般要求虛部非零,通常取為β=j/k。

考慮點P位于周期結(jié)構(gòu)表面上的情況,(12)式可變換為[14]:

(13)

在無網(wǎng)格方法中,周期結(jié)構(gòu)表面上任意一點處(例如Q)的聲壓可表示為:

(14)

式中:P1,P2,…,Pn表示n個節(jié)點,NQ為各節(jié)點在Q處的形函數(shù)向量,p為各節(jié)點處聲壓的向量。

假設(shè)已經(jīng)將周期結(jié)構(gòu)薄板模型用P1,P2,…,Pnn個節(jié)點進行了劃分,那么將(14)式代入(13)式,經(jīng)過化簡即可得到計算節(jié)點聲壓差的系統(tǒng)方程:

(C+D)·p=F

(15)

式中：C、D為n×n階矩陣,p、F為n×1階向量。

根據(jù)(15)式求得周期結(jié)構(gòu)薄板上所有節(jié)點處的聲壓p之后,將其帶入(13)式的變換形式,即可求得聲場中任意一點P的聲壓:

pω(R)=-B·p

(16)

對于某一聲源來說,根據(jù)上述推導(dǎo)過程求得所有接收點的聲壓之后,就可利用(7)式求得此聲源對應(yīng)入射方向的方向散射系數(shù);若進一步求得所有聲源的此系數(shù),就可以利用(8)式求得平均散射系數(shù)。

2.2 不同形式、尺度小樣品散射系數(shù)的數(shù)值計算

為研究何種尺度、形式的小樣品對大尺度周期結(jié)構(gòu)的散射系數(shù)具有代表性,本節(jié)對3種不同形式的小樣品的散射系數(shù)進行了計算,如圖2所示。

小樣品形式A、B、C各取自大尺度周期結(jié)構(gòu)的不同局部區(qū)域,其中B的周期數(shù)與大尺度結(jié)構(gòu)相同。各模型的子結(jié)構(gòu)為相同的三角形,其高為sub-h=6 cm,底邊長度為sub-l=16 cm。l為常數(shù),其值為l=15*sub-l=2.4 m。在計算中通過變化b,分別取b=0.3l,b=0.5l,b=0.7l,b=0.9l來代表不同尺度的周期結(jié)構(gòu),計算頻率為從400 Hz到2 000 Hz的1/3倍頻程的中心頻率。3種形式的小樣品在不同尺度下的散射系數(shù)如圖3～圖5所示。

圖2 取自大尺度周期結(jié)構(gòu)不同局部區(qū)域的小樣品形式

圖3 不同尺度,形式A小樣品的散射系數(shù) 圖4 不同尺度,形式B小樣品的散射系數(shù) 圖5 不同尺度,形式C小樣品的散射系數(shù)

從計算結(jié)果來看,形式B的小樣品散射系數(shù)與大尺度周期結(jié)構(gòu)的最為接近,即使在其尺度僅有大尺度周期結(jié)構(gòu)面積的0.3倍時,兩者數(shù)值依然十分接近。而形式A和C的數(shù)值則有明顯差別。若將大尺度周期結(jié)構(gòu)的散射系數(shù)視為標(biāo)準(zhǔn)參考值,各個尺度下3種形式的小樣品散射系數(shù)的相對誤差如圖6所示。

由圖6可以看出,隨著小樣品尺度的減小,形式A和形式C的相對誤差急劇增大,但是形式B的相對誤差始終保持在很低的水平,均小于5%。上述分析表明,只有如同形式B一樣,當(dāng)小樣品的周期數(shù)與大尺度周期結(jié)構(gòu)的相同時,小樣品的散射系數(shù)才具有等效性。

圖6 不同尺度時3種形式的小樣品散射系數(shù)的相對誤差

為了確定形式B的小樣品具有等效性的最小尺度,對各種尺度的形式B小樣品的方向入射散射系數(shù)做了進一步的分析。計算中共選取了從b=0.1l到b=0.9l的9種情況,以大尺度周期結(jié)構(gòu)的方向入射散射系數(shù)為參考標(biāo)準(zhǔn),按下式計算了每個頻率下所有入射方向的平均相對誤差,如圖7所示。

(17)

式中：δs,i、δf,i分別為小樣品及大尺度周期結(jié)構(gòu)的對應(yīng)于入射方向i的方向入射散射系數(shù),n表示聲源的數(shù)量。

圖7 不同頻率所有入射方向下的方向入射散射系數(shù)平均相對誤差

由圖7可以看出,隨著b對l倍數(shù)的減小,即小樣品尺度的減小,各個頻率上的入射方向散射系數(shù)平均相對誤差均有增大的趨勢,說明,小樣品的尺度與大尺度周期結(jié)構(gòu)越接近,其散射系數(shù)越具有等效性,這也與理論推導(dǎo)得到的結(jié)論一致。當(dāng)b≥0.5l時,各個頻率的平均相對誤差均小于10%,這說明對于形式B的小樣品,當(dāng)其尺度至少為大尺度周期結(jié)構(gòu)的一半時,才具有對大尺度周期結(jié)構(gòu)散射系數(shù)的等效性。

在計算效率方面,利用BMM對大尺度周期結(jié)構(gòu)的各頻率的散射系數(shù)進行計算,共耗時65 h,工作站參數(shù)為:Intel (R) Xeon (R) CPU E5-2620 @2.00 GHz 6 cores (2CPU) with 64.0GB RAM。而當(dāng)尺度減半時,即b=0.5l時,計算時間僅為約3 h。這說明尺度的減小可以使矩陣的規(guī)模急劇下降,從而大大提高計算效率。

上述計算結(jié)果表明,當(dāng)小尺度樣品的周期數(shù)與大尺度試件相同,且總面積至少為其1/2時,其散射系數(shù)可以作為大尺度試件的散射系數(shù)。為驗證此結(jié)論的通用性,本文對另外2種具有正弦型及矩形輪廓周期結(jié)構(gòu)進行了數(shù)值計算。2種周期結(jié)構(gòu)的周期數(shù)分別為20、10。在數(shù)值計算中,2種周期結(jié)構(gòu)的大尺度樣品與圖2中所示一致,小樣品的周期數(shù)與大樣品相同,但面積為其一半,計算所得的散射系數(shù)如圖8所示。

圖8 不同形狀周期結(jié)構(gòu)大小樣品散射系數(shù)對比

圖8中2種小樣品的散射系數(shù)與大尺度樣品十分接近,這說明對于不同輪廓、周期數(shù)的大尺度周期結(jié)構(gòu)來說,獲取其散射系數(shù)時均可用滿足上述條件的小樣品來等效。

3 結(jié) 論

本文對不同尺度周期結(jié)構(gòu)的散射性質(zhì)之間的聯(lián)系進行了理論分析和數(shù)值計算研究,首次給出了利用小尺度周期結(jié)構(gòu)的散射系數(shù)等效大尺度周期結(jié)構(gòu)的準(zhǔn)則。文中首先對不同尺度周期結(jié)構(gòu)散射系數(shù)之間的關(guān)系進行了理論分析,給出了利用小尺度周期結(jié)構(gòu)的散射系數(shù)來等效大尺度周期結(jié)構(gòu)的理論依據(jù);然后推導(dǎo)了計算散射系數(shù)的邊界無網(wǎng)格數(shù)值計算模型BMM,對不同形式小樣品的散射系數(shù)進行了計算,結(jié)果表明只有那些周期數(shù)與大尺度周期結(jié)構(gòu)相同的小樣品的散射系數(shù)具有等效性,通過對方向入射散射系數(shù)的進一步對比研究,發(fā)現(xiàn)當(dāng)小樣品的面積至少為大尺度周期結(jié)構(gòu)的一半時,散射系數(shù)具有很好的等效性。

概括起來,利用小尺度周期結(jié)構(gòu)來等效大尺度周期結(jié)構(gòu)散射系數(shù)需要滿足——小樣品與大尺度周期結(jié)構(gòu)有相同的周期數(shù),且面積至少為其一半。利用這一結(jié)論,將可大大提高散射系數(shù)的計算效率,同時對于工程應(yīng)用的快速估計乃至測量實驗的改進都具有重要的意義。

參考文獻：

[1] Toyoda M, Furukawa T, Takahashi D. Study on Echo Suppression Effect and Coloration due to Periodic-Type Diffusers[J]. Appl. Acoust, 2009, 70(5): 722-729

[2] Vorlander M, Mommertz E. Definition and Measurement of Random-Incidence Scattering Coefficients[J]. Appl Acoust, 2000, 60(2): 187-199

[3] ISO17497-1. Acoustics-Sound Scattering Properties of Surfaces-Measurement of the Random-Incidence Scattering Coefficient in a Reverberation Room[S]. International Organization for Standardization, 2004

[4] Ripoll J, Ntziachristos V, Carminati R, Nieto-Vesperinas M. Kirchhoff Approximation for Diffusive Waves[J]. Phys Rev E, 2001, 64(5): 051917

[5] Holford R L. Scattering of Sound Waves at a Periodic, Pressure-Release Surface: An Exact Solution[J]. J Acoust Soc Am, 1981, 70(4): 1116-1128

[6] Embrechts J J, Billon A. Theoretical Determination of the Random-Incidence Scattering Coefficients of Infinite Rigid Surfaces with a Perioidc Rectangular Roughness Profile[J]. Acta Acust Acust, 2011, 97(4): 607-617

[7] Kosaka Y, Sakumay T. Numerical Examination on Scattering Coefficients of Architectural Surfaces Using the Boundary Element Method[J]. Acoust Sci & Tech, 2005, 26(2): 136-144

[8] Embrechts J J, Geetere L D, Vermeir G, Vorlander M, Sakuma T. Calculation of the Random Incidence Scattering Coefficients of a Sine-Shaped Surface[J]. Acta Acust Acust, 2006, 92(4): 593-603

[9] Vorlander M, Embrechts J J, Geetere L D, Vermeir G, Gomes M H. Case Studies in Measurement of Random Incidence Scattering Coefficients[J]. Acta Acust Acust, 2004, 90(5): 858-867

[10] Lin F M, Hong P Y, Lee C Y. An Experimental Investigation Into the Sound-Scattering Performance of Wooden Diffusers with Different Structures[J]. Appl Acoust, 2010, 71(1): 68-78

[11] Mommertz E. Determination of Scattering Coefficients from the Reflection Directivity of Architectural Surface[J]. Appl Acoust, 2000, 60(2): 201-203

[12] Kuttruff H. Room Acoustics [M]. Fifth Edition London: Spon, 2009

[13] Burton A J, Miller G F. The Application of Integral Equation Methods to Numerical Solution of Some Exterior Boundary Problems[J]. Proc Roy Soc Lond A, 1971, 323(1553): 201-210

[14] Terai T. On Calculation of Sound Fields around Three Dimensional Objects by Integral Equation Method[J]. J Sound Vib, 1980, 69(1): 71-100