利用小樣品獲取大尺度周期結構聲散射系數的可行性分析

王海濤, 曾向陽, 劉延善

(西北工業大學 航海學院, 陜西 西安 710072)

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

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

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

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

1 理論分析

在混響環境中,根據散射系數的定義可得周期結構散射系數的計算公式:

(1)

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

(2)

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

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

(3)

圖1 測量實驗中不同尺度周期結構及轉臺示意圖

(4)

(5)

小樣品的散射系數可表示為:

(6)

可見,對于不同尺度的周期結構,無規入射吸聲系數是相同的,計算散射系數時僅有無規入射鏡面吸聲系數不同。

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

2 數值計算

2.1 邊界無網格型散射系數計算模型BMM

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

δ(θ,φ)=

(7)

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

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

(8)

因此,散射體及參考板接收點聲壓的計算是計算散射系數的關鍵。

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

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

(9)

式中：k=ω/c,為波數。

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

(10)

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

(11)

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

(12)

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

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

(13)

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

(14)

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

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

(C+D)·p=F

(15)

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

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

pω(R)=-B·p

(16)

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

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

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

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

圖2 取自大尺度周期結構不同局部區域的小樣品形式

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

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

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

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

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

(17)

式中：δs,i、δf,i分別為小樣品及大尺度周期結構的對應于入射方向i的方向入射散射系數,n表示聲源的數量。

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

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

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

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

圖8 不同形狀周期結構大小樣品散射系數對比

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

3 結 論

本文對不同尺度周期結構的散射性質之間的聯系進行了理論分析和數值計算研究,首次給出了利用小尺度周期結構的散射系數等效大尺度周期結構的準則。文中首先對不同尺度周期結構散射系數之間的關系進行了理論分析,給出了利用小尺度周期結構的散射系數來等效大尺度周期結構的理論依據;然后推導了計算散射系數的邊界無網格數值計算模型BMM,對不同形式小樣品的散射系數進行了計算,結果表明只有那些周期數與大尺度周期結構相同的小樣品的散射系數具有等效性,通過對方向入射散射系數的進一步對比研究,發現當小樣品的面積至少為大尺度周期結構的一半時,散射系數具有很好的等效性。

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

參考文獻：

[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