, 王颂 2 , 党建武 2,3 1.天津大学 软件学院, 天津 300350, 中国;
2.天津大学 计算机科学与技术学院, 天津 300350, 中国;
3.北陆先端科学技术大学院大学 信息科学学院, 石川923-1292, 日本
收稿日期: 2016-06-24
基金项目: 国家自然科学基金资助项目（51478305，61175016）；国家自然科学基金重点项目（61233009）；教育部归国留学人员科研启动基金资助项目
作者简介: 魏建国(1971-),男,副研究员
通信作者: 侯庆志,讲师,E-mail:qhou@tju.edu.cn
摘要：在人体发音过程仿真中，考虑声道边界的动态变化以及气流的流动，可以更加准确、真实地模拟声波在声道中的传播。在处理带有移动边界的气动声学问题时，相比传统声道声学研究中广泛应用的网格方法，无网格方法可以避免网格重构、网格畸变等。基于Euler体系下的气动声学波动方程，推导了Lagrange体系下声波传播的控制方程，并建立了无网格光滑粒子动力学（smoothed particle hydrodynamics，SPH）方法的数值离散格式。通过对比静止流体中声传播问题的SPH解和时域有限差分（finite difference time domain，FDTD）解，验证了SPH方法在声学计算中的准确性和可靠性。对于一维和二维流动流体中的声传播问题，通过与基于Doppler效应的理论解对比，阐明了利用SPH方法求解复杂气动声学问题的可行性。
关键词： 气动声学 声道 无网格 光滑粒子动力学 Lagrange方法
SPH simulations of aeroacoustic problems in vocal tracts
WEI Jianguo1, HAN Jiang2, HOU Qingzhi2
, WANG Song2, DANG Jianwu2,31.School of Computer Software, Tianjin University, Tianjin 300350, China;
2.School of Computer Science and Technology, Tianjin University, Tianjin 300350, China;
3.School of Information Science, Japan Advanced Institute of Science and Technology, Ishikawa 923-1292, Japan
Abstract:Simulation of human sound wave propagation need to take into account the moving boundaries and fluid flow within the vocal tract for accurate realistic models. Traditional mesh-based methods that are widely used to study human sound production have many problems due to mesh reconstruction and distortion, so they are not as effective as meshless methods. The aeroacoustic wave equations in the Eulerian framework are transformed to the governing equations for wave propagation in the Lagrangian form and discretized using the smoothed particle hydrodynamics (SPH) method. The accuracy and reliability of SPH for wave propagation in a static media are shown by comparisons with finite difference time domain (FDTD) results. This method is validated against the Doppler effect based theoretical solutions for one-and two-dimensional aeroacoustics to verify the ability of SPH to solve complex aeroacoustic problems.
Key words: aeroacousticsvocal tractmeshlesssmoothed particle hydrodynamicsLagrangian method
在语音生理学中，声音是声道中的空气质子将声带振动以波的形式传播出去产生的。 声道主要是指由唇部、 口腔、 咽腔以及喉腔共同形成的一个管状空间，从声门(声带)开始到唇部结束[1]。 在发音过程中，声道的形状决定着声音的音色以及声音的内容。 当声道形状改变时声音也会发生改变。 此外，当气流速度较高时，介质的流动也可对声波的传播产生较大影响。 因此，为了更加准确、 真实地模拟发音过程中声波在声道中的传播，移动边界条件下流动介质中的声传播问题研究尤为重要。 作为从神经控制到声学输出的语音产生过程的神经生理控制计算模型[2]的一部分，本文主要进行底层语音生成的建模与分析。 该计算模型旨在建立一个从上到下，从神经模型到生理模型再到语音模型的完整体系，并通过生成语音的反馈，自适应调节神经模型进而对底层模型进行控制[2]。
1954年，Lighthill[3]创立了声类比理论，标志着气动声学学科的诞生。 计算流体力学(computational fluid dynamics，CFD)的快速发展为气动声学的研究指明了方向，即利用数值计算的方法来求解气动声学问题，随即产生了计算气动声学(computational aeroacoustics，CAA)[4]。 针对CAA问题的数值计算，Tam等[5]提出了色散关系保持(dispersion relation preserving，DRP)差分格式。 由于DRP格式能够同时从流场中分辨出涡波、 声波和熵波，该方法获得了广泛应用。 另一种比较成熟的空间离散格式是Lele[6]提出的紧致有限差分格式。 基于波数分析理论,许多研究者对上述2种格式进行了优化,如Kim等[7]的优化紧致差分格式、 Zhong[8]的高阶迎风格式、 Zhuang等[9]的优化迎风格式、 Gaitonde等[10]的优化紧致有限体积格式、 Lee等[11]的紧致谱格式等。 其他在CAA领域广泛采用的差分格式还有Hu等[12]提出的低耗散低色散Runge-Kutta(low dissipation and low dispersion Runge-Kutta，LDDRK)差分格式。 随着有限元技术的广泛应用，有限元方法也被引入到CAA领域，比较成熟的是间断Galerkin有限元法(discontinuous Galerkin finite element method，DGFEM)[13]。
在模拟含移动边界的气动声学问题时，传统网格方法可能会引起网格的畸变，而网格畸变可造成计算精度的快速下降，甚至导致计算错误。 在一定程度上，网格重构技术可以解决这个问题，但该技术一般十分耗时，并且在频繁网格重构时的插值过程中，计算精度也会受较大影响。 本文采用的光滑粒子动力学(smoothed particle hydrodynamics，SPH)方法[14]是一种发展迅速且应用广泛的无网格方法。 相比传统网格方法，由于不产生网格，避免了网格畸变和网格重构引起的各类问题，对于模拟边界实时动态变化的声传播过程有重要意义。
1 控制方程声波的传播需要满足质量守恒、 动量守恒以及能量守恒定律，因此可以利用连续性方程、 动量方程以及能量方程共同描述声波的传播。 对于流体中的声传播问题，流体的黏性通常可以忽略，从而可采用Euler方程作为控制方程[15]。 根据观察视角的不同，Euler方程有Euler和Lagrange 形式这2种。
1.1 Euler形式Euler形式的控制方程为
| $\left\{ \begin{matrix} \frac{\partial v}{\partial t}+v\nabla v=-\frac{\nabla p}{{{\rho }_{0}}}, \\ \frac{\partial \rho }{\partial t}+v\nabla \rho =-{{\rho }_{0}}{{c}^{2}}_{0}\nabla v. \\\end{matrix} \right.$ | (1) |
1.2 Lagrange形式应用物质导数的定义
| $\frac{d\varphi }{dt}=\frac{\partial \varphi }{\partial t}+v\nabla \varphi ,$ |
| $\left\{ \begin{matrix} \frac{dv}{dt}=-\frac{\nabla p}{{{\rho }_{0}}}, \\ \frac{dp}{dt}=-{{\rho }_{0}}{{c}_{0}}^{2}\nabla v \\ \frac{dx}{dt}=v. \\\end{matrix} \right.,$ | (2) |
两种形式的主要区别在于： 式(1)的求解不需要在每个时间步更新计算点的位置，而式(2)的求解需要根据速度对计算点位置进行更新。 另外，Lagrange形式方程中不含对流项，可以避免传统网格法因对流项离散带来的各类数值问题。 本文采用式(2)作为控制方程。
2 光滑粒子动力学方法(SPH)由Lucy[16]和Gingold等[17]在1977年提出的SPH方法是公认的最早出现的一种无网格方法。 该方法最早应用于天体物理领域如原恒星或星系的形成与发展等。 20世纪90年代初，SPH方法的研究和应用变得越来越广泛，成功解决了很多问题如星系的形成、 高速碰撞、 金属成形、 铸造和工业加工、 自由表面流体、 低Reynolds数流体、 多相流、 热传导和血管流等[18]。
2.1 SPH方法的基本概念与原理SPH方法的基本思想是在计算域上任意设置有限个节点，采用节点权函数(或核函数)来表征节点及其邻域内的物理量，即利用节点权函数近似地表示其影响域内的场函数。 SPH方法采用基于点的近似，不受网格的约束，可以方便地在求解域内增加或减少节点，适用于求解具有复杂计算域的问题。
核函数近似、 粒子近似和核函数是SPH方法中的3个重要概念，以一维情况为例，分别简述如下，具体内容可参考文[14]。
2.1.1 核函数近似在SPH方法中，可用核函数积分的形式将场函数f(x)近似为
| $\left\langle f(x) \right\rangle ={{\int }_{\Omega }}f(\xi )W(x-\xi ,h)d\xi .$ | (3) |
参照上述场函数的核函数近似，其空间导数的核函数近似可以表示为
| $\begin{align} & \left\langle f(x) \right\rangle ={{\int }_{\Omega }}f(\xi )W(x-\xi ,h)d\xi . \\ & \left\langle \nabla f(x) \right\rangle ={{\int }_{\Omega }}\nabla f(\xi )W(x-\xi ,h)d\xi = \\ & {{\int }_{S}}f(\xi )W(x-\xi ,h)ndS- \\ & {{\int }_{\Omega }}f(\xi ){{\nabla }_{\xi }}W(x-\xi ,h)d\xi = \\ & {{\int }_{\Omega }}f(\xi ){{\nabla }_{x}}W(x-\xi ,h)d\xi =\text{ }\nabla \left\langle f(x) \right\rangle . \\ \end{align}$ | (4) |
2.1.2 粒子近似粒子近似是SPH方法中另一个重要的步骤。 用任意分布的N个粒子来表示计算域Ω，每个粒子的位置表示为xj，体积表示为ΔVj(j=1,2,… ,N)。
用Riemann加和的形式来近似式(3)和(4)中的积分，可得场函数及其导数的粒子近似表达式：
| $\left\langle f(x) \right\rangle =\sum\limits_{j}{f({{x}_{j}})}W(x-{{x}_{j}},h)\Delta {{V}_{j}}$ | (5) |
| $\left\langle \nabla f(x) \right\rangle =\sum\limits_{j}{f({{x}_{j}})}\nabla W(x-{{x}_{j}},h)\Delta {{V}_{j}}.$ | (6) |
2.1.3 核函数在SPH方法中，根据不同的问题可以选择不同的核函数。 这些核函数通常满足如下条件：
① 归一化：
② 对称性： W(xj-x,h)=W(x-xj,h)；
③ 非负性： W(x-xj,h)≥0；
④ 一致性： $\underset{h\to 0}{\mathop{\lim }}\,$W(x-xj,h)=δ(x-xj)；
⑤ 递减性： 随距离r=x-xj的增大而递减；
⑥ 导数归零化：$\sum\limits_{j}{\nabla W(x-{{x}_{j}},h)}\Delta {{V}_{j}}=0$.其中δ(x-xj)是Dirac δ函数。
本文选用如下的三次样条函数作为核函数：
| $W=\frac{G}{{{h}^{\lambda }}}\left\{ \begin{matrix} 1-1.5{{q}^{2}}+0.75{{q}^{3}}, & 0\le q ＜1; \\ 0.25{{\left( 2-q \right)}^{3}}, & 1\le q＜2; \\ 0, & q\ge 2. \\\end{matrix} \right.$ |
核函数的一阶导函数表示如下：
| $\frac{\partial W(x-{{x}_{j}},h)}{\partial {{x}_{\alpha }}}=\frac{1}{h}\frac{dW}{dq}\frac{{{x}_{\alpha }}-{{x}_{j,\alpha }}}{r}.$ |
2.2 空间导数的离散化应用SPH导数近似公式(6)以及核函数导数的归零化条件⑥，对式(2)中的空间导数进行离散，得到如下的半离散格式：
| $\left\{ \begin{matrix} {{\rho }_{0}}\frac{d{{v}_{i}}}{dt}=-\sum\limits_{j}{\frac{{{m}_{j}}}{{{\rho }_{j}}}}({{p}_{j}}-{{p}_{i}}){{\nabla }_{i}}{{W}_{ij}}, \\ \frac{1}{c_{0}^{2}}\frac{d{{p}_{i}}}{dt}=-{{\rho }_{0}}\sum\limits_{j}{\frac{{{m}_{j}}}{{{\rho }_{j}}}}({{v}_{j}}-{{v}_{i}}){{\nabla }_{i}}{{W}_{ij}}, \\ \frac{d{{x}_{i}}}{dt}={{v}_{i}}. \\\end{matrix} \right.$ | (7) |
对于式(7)中时间项的积分，本文采用二阶Runge-Kutta(RK2)方法。
2.3 数值边界条件在数值计算中，边界条件为一系列针对计算域边界的附加条件。 边界条件主要包括Dirichlet、 Neumann以及Robin 3种类型。 对于入口边界，本文采用Dirichlet型，即给定入口处粒子的压强和速度； 出口边界为Neumann型自由出口，即速度和压力的梯度为零。 为了施加给定的边界条件，并克服SPH方法的边界缺陷问题[14]，本文使用受到广泛验证的虚粒子方法[19-20]。
3 数值实验与讨论在本节中，通过使用上述的SPH方法求解式(3)，模拟声波在一维和二维空间中的传播，并通过与文献结果对比，评价本文方法的数值特性。
3.1 一维算例 对于一维问题，在空间[0,1]内均匀分布501个粒子并在中心处产生如下Gauss声源：
| $gp(x)={{e}^{-\left| \frac{x-{{x}_{c}}}{20\Delta x} \right|}}.$ |
在t为250和500 μs时刻的仿真结果分别如图 1和2所示。 作为对比，图中同时给出时域有限差分(finite difference time domain,FDTD)的计算结果。 FDTD是计算声学中应用最广泛的方法之一[21]。 由图 1和2可观察到，SPH结果与FDTD结果吻合较好，从而验证了SPH方法的可靠性。 同时也表明： 对于静止介质中的波传播问题，无网格SPH方法可以获得与传统网格方法相一致的结果。
 |
| 图 1 一维静止介质中250 μs时的波形 |
| 图选项 |
 |
| 图 2 一维静止介质中500 μs时的波形 |
| 图选项 |
为了表明本文方法的优越性，声波在流动气体中的传播即气动声学问题得到了研究。 与上例相比，唯一不同之处在于流体的初始速度为200 m/s，也就是说流体的Mach数为Ma=v/c0=0.59。 对于该类问题，式(1)中的对流项不容忽视。
在t为250和500 μs时刻的仿真结果分别如图 3和4所示。 为表明流体流动对声波传播的影响，上例中流体不流动条件下的波传播结果也分别标在两幅图中。 从图中可观察到，相同时间内，声波沿流体流动方向的传播距离明显增加，而沿流体流动反方向的传播距离明显减小。 从而表明，介质的流动可对声波的传播产生重要影响。 另外，计算所得波形与向流动方向平移流体不流动时的波形(平移距离vt，考虑Doppler效应的波形)吻合较好，这也从侧面验证了本文计算结果的正确性。
 |
| 图 3 一维流动介质中250 μs时的波形 |
| 图选项 |
 |
| 图 4 一维流动介质中500 μs时的波形 |
| 图选项 |
3.2 二维算例 对于二维问题，在长为1 m的方形计算域内均匀分布10 201个粒子。 初始条件下，粒子的压力和速度都为零，声源位于方形空间中央。 计算域的上下边界为自由滑移边界，采用镜像粒子法[20]予以处理。
在t=200 μs和t=400 μs时刻的仿真结果分别如图 5和6所示。 可以看出，二维空间内SPH解和FDTD解吻合较好，从而验证了SPH方法在高维空间波传播问题模拟中的可靠性。 同时这也表明，对于二维静止介质中的波传播问题，无网格SPH方法可以获得与传统网格方法相一致的结果。
 |
| 图 5 二维静止介质中200 μs时的波形 |
| 图选项 |
 |
| 图 6 二维静止介质中400 μs时的波形 |
| 图选项 |
同样，二维空间内的气动声学问题也得到了仿真。 与上例相比，流体在x轴方向的速度为200 m/s。 在t=200 μs和t=400 μs时刻的仿真结果分别如图 7和8所示。 作为参考，流体不流动条件下的波传播结果也分别标在两幅图中。 由图示结果可以获得与一维算例相类似的结论。
 |
| 图 7 二维流动介质中200 μs时的波形 |
| 图选项 |
 |
| 图 8 二维流动介质中400 μs时的波形 |
| 图选项 |
另外，与FDTD相比，SPH可基于不规则的粒子分布构建计算格式，从而具备处理不规则边界的便利性。 虽然所给算例中的边界是规则的，但是由于粒子的动态变化，实际计算是基于非规则粒子分布的(即使初始粒子分布是均匀规则的，随着粒子的移动，也将变得不规则)。
SPH与FDTD的计算精度已通过上述算例进行了对比，下面再讨论一下二者的计算特点。 由于2种方法在不同的框架体系下完成(SPH为Lagrange形式，FDTD为Euler形式)，二者的计算效率不具有严格的可比性。 FDTD利用固定网格点构造计算格式，效率较高； 而SPH基于位置动态变化的粒子，计算格式随时间的变化而变化，为了建立新的格式，每个时间步需要根据核函数半径进行粒子搜索，因此效率不如FDTD高。 需要指出的是，如果不考虑粒子的移动，SPH将退化为有限差分方法，从而具有与FDTD相近的计算效率。
4 结 论本文推导了Lagrange体系下的波传播方程，并利用无网格SPH方法对控制方程进行空间离散，利用二阶Runge-Kutta方法进行时间离散，建立了气动声学计算的无网格SPH格式。 通过对一维和二维静止介质中声波传播问题模拟结果的对比，发现SPH数值解和FDTD解吻合较好，从而验证了SPH方法在计算波传播问题上的有效性。 通过模拟一维和二维空间内流动介质中的波传播过程，并与基于Doppler效应的数值解进行对比，验证了SPH方法在模拟气动声学问题上的优越性和准确性。 下一步将研究移动边界条件下流动介质中的声波传播问题。
参考文献
| [1] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Stevens K N. Acoustic Phonetics[M].Cambridge: The MIT Press, 2000. |
| [2] | Journal of Central South University(Science and Technology), 41(2):649-654.--> CHEN Xi, DANG Jianwu, YAN Han, et al. A neural understanding of speech motor learning[C]//Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA 2013). Kaohsiung, Taiwan, China:IEEE Press, 2013:321-325. |
| [3] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Lighthill M J. On sound generated aerodynamically:II. Turbulence as a source of sound[J]. Proc R Soc Lond A, 1954, 222 : 1–32.DOI:10.1098/rspa.1954.0049 |
| [4] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Tam C K W. Computational aeroacoustics:Issues and methods[J]. AIAA J, 1995, 33(10) : 1788–1796.DOI:10.2514/3.12728 |
| [5] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Tam C K W, Webb J C. Dispersion-relation-preserving finite difference scheme for computational acoustics[J]. J Comput Phys, 1993, 107(2) : 262–281.DOI:10.1006/jcph.1993.1142 |
| [6] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Lele S K. Compact finite difference scheme with spectral-like resolution[J]. J Comput Phys, 1992, 103(1) : 16–42.DOI:10.1016/0021-9991(92)90324-R |
| [7] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Kim J W, Lee D J. Optimized compact finite difference schemes with maximum resolution[J]. AIAA J, 1996, 34(5) : 887–893.DOI:10.2514/3.13164 |
| [8] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Zhong X L. High-order finite-difference schemes for numerical simulation of hypersonic boundary-layer transition[J]. J Comput Phys, 1998, 144(2) : 662–709.DOI:10.1006/jcph.1998.6010 |
| [9] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Zhuang M, Chen R F. Optimized upwind dispersion relation preserving finite difference scheme for computational aeroacoustics[J]. AIAA J, 1998, 36(11) : 2146–2148.DOI:10.2514/2.319 |
| [10] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Gaitonde D, Shang J S. Optimized compact difference based finite-volume schemes for linear wave phenomena[J]. J Comput Phys, 1997, 138(2) : 617–643.DOI:10.1006/jcph.1997.5836 |
| [11] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Lee C, Seo Y. A new compact spectral scheme for turbulence simulation[J]. J Comput Phys, 2002, 183(2) : 438–469.DOI:10.1006/jcph.2002.7201 |
| [12] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Hu F Q, Hussaini M Y, Manthey J L. Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics[J]. J Comput Phys, 1996, 124(1) : 177–191.DOI:10.1006/jcph.1996.0052 |
| [13] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Cockburn B, Shu C W. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II:general framework[J]. Math Comput, 1989, 52 : 411–435. |
| [14] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Liu G R, Liu M B. Smoothed Particle Hydrodynamics Methed:A Meshfree Pariticle Method[M].Singapore: World Scientific Press, 2003. |
| [15] | Journal of Central South University(Science and Technology), 41(2):649-654.-->张德良. 计算流体力学教程[M].北京: 高等教育出版社, 2010.ZHANG Deliang. A Course in Computational Fluid Dynamics[M].Beijing: Higher Education Press, 2010.(in Chinese) |
| [16] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Lucy L B. A numerical approach to the testing of the fission hypothesis[J]. Astron J, 1977, 8(12) : 1013–1024. |
| [17] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Gingold R A, Monaghan J J. Smoothed particle hydrodynamics:Theory and applications to non-spherical stars[J]. MNRAS, 1977, 18 : 375–389. |
| [18] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Monaghan J J. Smoothed particle hydrodynamics and its diverse applications[J]. Ann Rev Fluid Mech, 2012, 44(44) : 323–346. |
| [19] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Lastiwka M, Basa M, Quinlan N J. Permeable and non-reflecting boundary conditions in SPH[J]. Int J Numer Methods Fluids, 2009, 61 : 709–724.DOI:10.1002/fld.v61:7 |
| [20] | Journal of Central South University(Science and Technology), 41(2):649-654.-->HOU Qingzhi, Kruisbrink A C H, Pearce F, et al. Smoothed particle hydrodynamics simulations of flow separation at bends[J]. Comput Fluids, 2014, 90(4) : 138–146. |
| [21] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Sullivan D M. Electromagnetic Simulation Using the FDTD Method[M].New York: John Wiley & Sons, 2013. |
