流域网格中嵌入式复杂声源的时域传播仿真

申肖雪，卢浩，肖友刚



流域网格中嵌入式复杂声源的时域传播仿真

申肖雪，卢浩，肖友刚

(中南大学交通运输工程学院，湖南长沙 410075)

针对声学商业软件较难模拟任意形状、时变复杂声源的声辐射问题，使用有限体积法在时域内求解无声源项亥姆霍兹(Helmholtz)方程，将复杂声源嵌入到有限体积单元节点，推导了由给定声源表面声压或振动位移得到速度势公式，提高了声源处理的灵活性和计算效率。该方法允许对初始场问题及复杂时变声源声辐射进行仿真。对常见声源及二阶圆柱体声源声辐射进行了数值模拟，结果与解析解及商业软件结果进行了对比，误差均小于15%。程序具有良好的封装性及通用性，可以灵活地对不同声源进行组合，得出任意复杂声源时域的传播规律，为复杂声源声辐射等线性声学问题的研究提供了一个可靠的平台。

复杂声源；时域；有限体积法；声辐射

0 引言

要研究复杂声源在复杂密闭空间内的辐射声波特性，如何容易、方便地设置声源便成为必须要面对的问题。目前声学商业软件中声源类型大多为单/偶/四极子声源、平面波声源、混响声源等，而对复杂形状、任意变化的时变声源则较难进行仿真，在研究论文中也比较少见。对于声源的仿真及声场计算，MILLAN[1]使用了有限元素法(Finite Element Method，FEM)和边界元法(Boundary Element Method，BEM)对时域内声波传播进行仿真，对单极子声源的处理是在给定球形网格并在球形表面设定位移。Takuya Oshima[2]开发了一款时域内声波传输的通用求解器potentialWaveFoam。Schmalz J[3]在OpenFOAM[4-5]下实现了Lighthill和Curle声类比，声源作用分别体现在Lighthill方程及FW-H方程。Nilsson[6]在OpenFOAM和CALFEM中实现声类比，对声源的处理是在求解器每一个时间步计算完成后根据公式计算声源项。Ji H[7]使用有限体积法计算声传播，对声源的处理是将声源项加入Helmholtz方程，然后对整个方程进行离散求解。

以上文献中无声源项的Helmholtz方程适用于初始场问题的求解，而对于时变声源，则是将声源项加入到声波方程，然后对其进行离散、矩阵求解等处理，大大占用了计算资源。另外，大部分声学商业软件可以对简易声源进行处理，但对复杂形状或复杂时变声源的处理技术较为缺乏。本文将声源以边界条件的形式巧妙地嵌入到数值计算的网格中，将求解时域内含声源项的Helmholtz方程的求解问题转化为求解无声源项Helmholtz方程，将声源以边界条件的形式加入到计算过程。在计算初期对流域划分网格时即考虑声源作用及位置，然后根据声源类型的不同在现有计算域网格中“掏去”不同形状的网格，使之形成边界，然后再在此边界上设定速度势，以模拟声源。在提高对复杂声源处理的灵活性的同时，提高了计算效率。

1 声波控制方程

理想气体的Helmholtz方程[8]为

将式(2)～(4)代入式(5)，可得

由于本文侧重于声源的实现方式，不对其方程的离散方法及修正算法进行细述。采用有限体积法计算声波的传播，对方程(1)进行空间离散，对式(6)进行时间离散，离散格式及修正算法参考文献[2]，算法迭代过程、误差分析见文献[9]。

2 嵌入式网格表面设置速度势作声源

空气噪声声源可以看作是单极子、偶极子、四极子声源的叠加。此节对时域内单极子、偶极子和四极子声源的声辐射进行计算，将结果进行FFT变换，并与商业软件的频域计算结果进行对比。然后对复杂时变二阶圆柱声源声辐射进行了时域仿真，结果与理论值吻合。

2.1 单极子声源

2.1.1 单极子声源

如图1所示，在边长为1 m的立方体声腔中心存在一个直径为0.01 m的脉动球形声源。脉动球源是进行均匀涨缩振动的球面声源，也就是在球源表面上各点沿着径向做同振幅、同相位的振动。

图1 在立方体声腔中心的半径为0.01m的球形声源

其中：为声源的强度，定义=0.2 N/m，=0.01 m，则球形声源的表面上的声压为20 Pa。

2.1.2 仿真模型

在立方体中心有一球形声源，由于整个模型关于立方体中心点对称，所以取立方体的1/4和球形声源的1/4进行仿真即可，如图2所示。原立方体的3个表面都设置为吸声边界，吸声系数为0.01，根据吸声系数的定义有[9]

图2 计算区域网格，原尺寸的1/4

Fig.2 The meshes in computation region, 1/4 of origin geometry

使用OpenFOAM自带的groovyBC边界条件设置工具，设置随时间变化的边界条件，在声腔左侧采用Dirichlet边界条件，根据式(2)的离散形式：

其中：为声波振幅，为声波频率。这样就能保证在球形表面上的声压以正弦方式变化。

2.1.3 计算结果

在初始时刻，声场中各点的声压值为0，声源声压在初始时刻开始以正弦方式变化，模拟时间为1.5 s。P点声压变化如图3所示。从图3可以看出，在1.0 s时间后，声压趋于稳定。

计算时间步长为10-6s，每20个时间步长取一次数据，总共有7 500个数据。取其中最后的2 048个时域数据进行FFT变换，得到其频域谱如图4所示。从图4可以看到，该点的声压在200 Hz时为0.17 Pa，商业软件中有限元法得到的该点声压为0.18 Pa，误差为5.6%。考虑到信号的采样频率及采样长度，频谱中幅值最大值所在的频率在200 Hz左右存在偏移。由于截断信号并非整周期，频谱在200 Hz以外发生频率泄漏。根据单极子声源的声压变化规律，设置边界上速度势，模拟单极子声源的声辐射。

图3 1.5 s内P点的声压变化

图4 时域数据的FFT变换

2.2 偶极子声源、四极子声源

本文对偶极子和四极子声源的仿真同样采用在流域声腔网格内嵌入球形网格的方式进行处理。在边长为1 m的立方体声腔中心，偶极子声源由两个声压幅值相同、相位相反的关于中心对称的脉动球源组成，两个声源的球心位置为(0.5，0.5，0.487 5)和(0.5，0.5，0.512 5)，半径为0.006 25 m，球形表面设置声压幅值均为20 Pa，频率为200 Hz，但是相位相反。声腔表面设为无反射边界，在时域内对声场进行计算。在声压稳定后，两声源的中垂面上由于两声源幅值相等而相位相反，声压相互抵消为0。在两球源中心连线上的声压叠加增强，声压最大。四极子声源由一对极性相反的偶极子组成。同单极子声源一样，使用snappyHexMesh命令指定四个球源的球心为(0.5，0.5，0.487 5)、(0.5，0.5，0.512 5)和(0.5，0.487 5，0.5)、(0.5，0.512 5，0.5)，半径为0.006 25 m及相关划分网格参数，划分好后施加边界条件即可将声源以边界条件的方式容易地加入到数值计算。同之前单极子声源的处理相似，将偶极子声源和四极子声源辐射声压时域内的结果经FFT变换后，同解析解及商业软件频域结果进行对比，结果吻合。使用嵌入式声源方法大大简化了计算过程。

对比MATLAB中k-wave工具箱中对偶极子、四极子声源的处理，将声源项加入二阶Helmholtz方程[10]，得

2.3 复杂时变声源

对简易声源的声辐射理论研究有不少[11]，但对复杂时变声源的数值计算研究并不多。在数值计算中，大多数将声源项加入声波方程然后将整个方程进行格式离散。对于单、偶极子等简单声源可以将其在网格点上设置为点源，将数据导入到计算过程。但是对于复杂声源(尤其近场辐射问题)，不能在方程及求解中设定几个网格点以代表整个声源，其表面每个点的变化规律也可能不同，因此需要尽量细致地描述声源。另外，对于复杂声源的声辐射仿真大多使用边界元、无限元或格林函数方法，要求仿真者对声辐射理论公式等有一定的了解，相关的数值计算方法不多。本文将任意形式的声源辐射仿真问题简化为求解无声源项的Helmholtz方程问题，将求解过程中的声源处理问题转化到前处理过程，使用者只需关注前处理过程，这样使得求解器具有良好的通用性。二阶圆柱体声源变化规律较单极子、偶极子等声源复杂，本节推导了极坐标系下圆柱体表面上速度势与位移关系式，对给定位移变化的二阶振动圆柱体声辐射案例进行仿真，并对理论值和数值解进行了对比。

2.3.1 问题描述与解析解

本小节对圆柱体振动辐射声场进行了仿真，对数值计算结果和理论值进行了对比。假设有一无限长圆柱体[9]，在柱状坐标系下，半径随角度的变化规律为

径向速度的表达式可写为

其中：为振动圆柱体的平均半径，为圆柱体径向振动幅值。根据文献[12]和文献[13]，阶圆柱体声源辐射声场的表达式为

2.3.2 圆柱体声源辐射仿真

已知圆柱体表面的位移变化，推导声源表面速度势变化。当圆柱体振动时，圆柱体表面速度与空气粒子速度相同。根据式(3)，圆柱体面上各点径向振动速度等于速度势的负梯度，在极坐标中可以表示为

在OpenFOAM中使用blockMesh工具对边长为1 m的立方体划分为20*20*20=8 000个六面体单元。然后使用cellSet和refineMesh工具对圆柱体所在的中心区域进行4层网格细化，最后使用snappyHexMesh工具划出半径为0.01 m、高度为1 m的圆柱体网格，并将其移去从而形成圆柱体边界。划分好后，将圆柱面的速度势使用groovyBC设置成为时变量，程序代码对应为

outCylinder_region0

{ type groovyBC;

valueExpression "-(cos(2*pi*200*time())-1)*2*pi*200*

(pow((pos().x-0.5),2)-pow((pos().y-0.5),2))";

//后半段代表cos(2)

value uniform 0;

gradientExpression "0";

}

其他边界设置为无反射边界，以模拟圆柱体声源在无限空间的声辐射。

2.3.3 数值计算结果与解析值对比

(a)= 0.01 m

(b)= 0.32 m

图5 距离圆柱体轴心= 0.01 m和= 0.32 m处的声压变化

Fig.5 The sound pressures at distance= 0.01 m and= 0.32 m in time domain

由图5可见，在= 0.01 m处，声压幅值为186 Pa。在= 0.32 m的地方，声压幅值为0.21 Pa，理论值为0.18 Pa，相差14.3%，吻合度较高。符合声压与距离平方成反比的规律。

3 结论

(1)本文提出了通过在边界上设置速度势的方式将声源引入数值计算的方法，推导了由给定的声源表面声压或振动位移得到速度势的公式。该方法在建立网格初期即考虑到声源的作用，声源以边界条件的方式出现，求解器只需要求解Helmholtz方程，不需要将声源项加入到Helmholtz方程，使得声场仿真变得简单、直观。

(2)在时域内对单/偶/四极子声源等简易声源及复杂时变二阶圆柱体声源的声辐射进行了数值计算，单极子声源声辐射结果与有限元结果相差5.6%，圆柱体声源的声辐射结果与解析解相差14.3%，误差较小。为更复杂声源的声辐射及声场仿真提供了思路。

(3)该方法将前处理、求解、后处理三个过程全部分开，程序具有良好的封装性及通用性，可用于研究任意形状、任意变化的物体的声辐射及传播等线性声学问题。

(4) 该方法对于远场仿真具有较好效果，而对近场仿真尚有一定的不准确性。方程(7)和方程(15)是自由声场下的声源声辐射公式。假设空间中存在反射波，即非自由声场下，以算例1为例，=0.01 m的球形表面声压为20 Pa，而商业软件得到的=0.01 m的结果并非20 Pa，而是有一定的误差。说明存在壁面反射时，近场声场仿真并不精确。

[1] MILLAN BEV. Acoustic time-domain simulation with BEM and FEM[D]. Stuttgart: University of Stuttgart, 2011: 36-39.

[2] OSHIMA T, IMANO M. A full finite-volume time-domain approach towards general-purpose code development for sound propagation prediction with unstructured mesh[C]//Proceedings of Inter-Noise 2008, Shanghai, 2008, 1511-1525.

[3] SCHMALZ J, KOWALCZYK W. Implementation of acoustic analogies in openFOAM for computation of sound fields[J]. Open Journal of Acoustics, 2015, 5(2): 29-44.

[4] JASAK H. OpenFOAM: open source CFD in research and industry[J]. International Journal of Naval Architecture and Ocean Engineering, 2015, 1(2): 89-94.

[5] WELLER G, TABOR G, JASAK H. A tensorial approach to computational continuum mechanics using object-oriented techniques[J]. Computers in Physics, 1998, 12(6): 620-631.

[6] NILSSON, J. Implementation of acoustical analogies in openFOAM and CALFEM[D]. Lund: Lund University, 2010: 13-14.

[7] JI H, XU X, GUO X. Direct FVM simulation for sound propagation in an ideal wedge[J]. Shock & Vibration, 2016, 2016(9): 1-9.

[8] 杜功焕, 朱哲民, 龚秀芬. 声学基础[M]. 3版, 南京: 南京大学出版社, 2012. DU Gonghuan, ZHU Zhemin, GONG Xiufen. Fundamentals of Acoustics[M]. 3rdEdition, Nanjing: Nanjing University Press, 2012, 147-152.

[9] MECHEL P. Formulas of Acoustics[M]. Springer Berlin Heidelberg, 2008, 147-152.

[10] TREEBY B E, JAROS J, ROHRBACH D. Modelling Elastic Wave Propagation Using the K-Wave MATLAB Toolbox[C]// 2014 IEEE International Ultrasonics Symposium,Chicago, 2014, 146-149.

[11] 曹景记, 唐晓明, 苏远大. 多极声源在裸眼井外的纵横波辐射特性[J]. 声学技术, 2016, 35(6): 487-492. CAO Jingji, TANG Xiaoming, SU Yuanda. Radiation characteristics of a multipole acoustic source in an open borehole[J]. Technical Acoustics, 2016, 35(6): 487-492.

[12] RUSSELL A. On the sound field radiated by a tuning fork[J]. American Journal of Physics, 2000, 68(12): 1139-1145.

[13] JACOBSEN F, JUHL P. Radiation of Sound[R]. Lyngby: Technical University of Denmark, 2011, 25-27.

Time domain propagation characteristics of complex embedded acoustic sources in fluid meshes

SHEN Xiao-xue, LU Hao, XIAO You-gang

(School of Traffic and Transportation Engineering, Central South University, Changsha 410075, Hu’nan, China)

Commercial acousticsoftware is generally hard to simulate sound radiation of complex acoustic sources, which are of arbitrary shape or time-varying. In order to solve this problem, the Helmholtz equation (without acoustic source term) is solved in time domain via Finite Volume Method with embedding acoustic sources to Finite-volume nodes, and the relationship between velocity potential and the pressure/displacement of acoustic source is acquired. This method makes it more convenient and efficient to deal with the acoustic source in numerical calculation. In addition, it allows simulating initial value problem and time-varying source problem. The sound radiation of common acoustic sources and the second order cylindrical acoustic source has been simulated, and the results agree well with analytic solutions and results of commercial code, with numerical errors smaller than 15%. This code has good encapsulation and versatility, via which different types of acoustic sources can be combined to form an arbitrary complex acoustic source and the propagation characteristics can be obtained. It is a good platform of carrying out research on linear acoustics like sound radiation of complex acoustic sources.

complex acoustic source; time domain; Finite Volume Method(FVM); sound radiation

TB533

A

1000-3630(2017)-05-0405-05

10.16300/j.cnki.1000-3630.2017.05.002

2017-01-18;

2017-02-28

国家自然科学基金资助项目(50975289、51275531)

申肖雪(1993－), 男, 河北景县人, 硕士研究生，研究方向为铁路噪声预测及控制。

肖友刚, E-mail: csuxyg@163.com