利用小样品获取大尺度周期结构声散射系数的可行性分析

王海涛, 曾向阳, 刘延善

(西北工业大学 航海学院, 陕西 西安 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