在语音生理学中，声音是声道中的空气质子将声带振动以波的形式传播出去产生的。 声道主要是指由唇部、 口腔、 咽腔以及喉腔共同形成的一个管状空间，从声门(声带)开始到唇部结束^[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)

其中： ρ₀表示介质的密度； v表示速度矢量； p表示压强； c₀表示声音在空气中的传播速度，本文取c₀=340 m/s。 传统网格方法大多直接求解式(1)。

1.2 Lagrange形式

应用物质导数的定义

$\frac{d\varphi }{dt}=\frac{\partial \varphi }{\partial t}+v\nabla \varphi ,$

可将式(1)转化为Lagrange形式：

$\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)

其中x表示计算域中任意一点的位置向量。

两种形式的主要区别在于： 式(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)

其中： x和ξ表示计算域内点的位置； 〈f(x)〉表示场函数f(x)在x处的核估计值； Ω表示计算域； f(ξ)为坐标ξ处的场量值； W(x-ξ,h)为核函数，它有2个自变量： 粒子间距离x-ξ和光滑长度h。

参照上述场函数的核函数近似，其空间导数的核函数近似可以表示为

$\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)

其中： ∇表示梯度算子，n表示表面S的外法向单位向量。 〈∇f(x)〉=∇〈f(x)〉表明函数导数的SPH核近似具有交换律的性质。

2.1.2 粒子近似

粒子近似是SPH方法中另一个重要的步骤。 用任意分布的N个粒子来表示计算域Ω，每个粒子的位置表示为x_j，体积表示为ΔV_j(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)

在SPH方法中，粒子体积ΔV_j通常定义为粒子质量与密度的比值即ΔV_j=m_j/ρ_j。 在计算过程中粒子质量保持不变，密度和体积可变。

2.1.3 核函数

在SPH方法中，根据不同的问题可以选择不同的核函数。 这些核函数通常满足如下条件：

① 归一化： $\sum\limits_{j}{W(x-{{x}_{j}},h)}\Delta {{V}_{j}}=1$;

② 对称性： W(x_j-x,h)=W(x-x_j,h)；

③ 非负性： W(x-x_j,h)≥0；

④ 一致性： $\underset{h\to 0}{\mathop{\lim }}\,$W(x-x_j,h)=δ(x-x_j)；

⑤ 递减性： 随距离r=x-x_j的增大而递减；

⑥ 导数归零化：$\sum\limits_{j}{\nabla W(x-{{x}_{j}},h)}\Delta {{V}_{j}}=0$. 其中δ(x-x_j)是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.$

其中： q=|x-x_j|/h表示相对距离，λ表示问题的维度。 对于一维和二维问题，标准化系数G的取值分别为2/3和10/(7π)。 三次样条函数获得了SPH研究者的广泛认可^[18]。

核函数的一阶导函数表示如下：

$\frac{\partial W(x-{{x}_{j}},h)}{\partial {{x}_{\alpha }}}=\frac{1}{h}\frac{dW}{dq}\frac{{{x}_{\alpha }}-{{x}_{j,\alpha }}}{r}.$

其中α=1,2,3，表示维度信息。

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)

其中： v_i,p_i和x_i分别表示待计算点的速度、 压强和位置坐标； m_j,ρ_j,v_j和p_j分表表示待计算点周围粒子的质量、 密度、 速度和压强。

对于式(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|}}.$

其中： x_c表示声源的中心位置； x表示计算域中任意一点的位置坐标； Δx表示该点与声源中心的距离。 初始条件下，粒子的压力和速度均为零。

在t为250和500 μs时刻的仿真结果分别如图 1和2所示。 作为对比，图中同时给出时域有限差分(finite difference time domain,FDTD)的计算结果。 FDTD是计算声学中应用最广泛的方法之一^[21]。 由图 1和2可观察到，SPH结果与FDTD结果吻合较好，从而验证了SPH方法的可靠性。 同时也表明： 对于静止介质中的波传播问题，无网格SPH方法可以获得与传统网格方法相一致的结果。

为了表明本文方法的优越性，声波在流动气体中的传播即气动声学问题得到了研究。 与上例相比，唯一不同之处在于流体的初始速度为200 m/s，也就是说流体的Mach数为Ma=v/c₀=0.59。 对于该类问题，式(1)中的对流项不容忽视。

在t为250和500 μs时刻的仿真结果分别如图 3和4所示。 为表明流体流动对声波传播的影响，上例中流体不流动条件下的波传播结果也分别标在两幅图中。 从图中可观察到，相同时间内，声波沿流体流动方向的传播距离明显增加，而沿流体流动反方向的传播距离明显减小。 从而表明，介质的流动可对声波的传播产生重要影响。 另外，计算所得波形与向流动方向平移流体不流动时的波形(平移距离vt，考虑Doppler效应的波形)吻合较好，这也从侧面验证了本文计算结果的正确性。

3.2 二维算例

对于二维问题，在长为1 m的方形计算域内均匀分布10 201个粒子。 初始条件下，粒子的压力和速度都为零，声源位于方形空间中央。 计算域的上下边界为自由滑移边界，采用镜像粒子法^[20]予以处理。

在t=200 μs和t=400 μs时刻的仿真结果分别如图 5和6所示。 可以看出，二维空间内SPH解和FDTD解吻合较好，从而验证了SPH方法在高维空间波传播问题模拟中的可靠性。 同时这也表明，对于二维静止介质中的波传播问题，无网格SPH方法可以获得与传统网格方法相一致的结果。

同样，二维空间内的气动声学问题也得到了仿真。 与上例相比，流体在x轴方向的速度为200 m/s。 在t=200 μs和t=400 μs时刻的仿真结果分别如图 7和8所示。 作为参考，流体不流动条件下的波传播结果也分别标在两幅图中。 由图示结果可以获得与一维算例相类似的结论。

另外，与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方法在模拟气动声学问题上的优越性和准确性。 下一步将研究移动边界条件下流动介质中的声波传播问题。