声学方程的有限元求解

1 声波方程

空气中的压力变动传达到人耳中表现为声音。所以声波方程其实就是空气的振动方程，可以从Navier-Stokes方程推导出来。由于空气粘性极小，因此空气的运动为无旋的纵波。此时的Navier-Stokes方程为 

再考虑空气的本构方程

在这里c为声速。

把方程（2）带入方程（1-1）得

对方程（3）两边取 D/Dt ，对方程（1-2）两边取 ∇ 然后相加得

由于速度值v的变化Dv/Dt一般较小，我们忽略左边第二项。得

再进一步忽略D/Dt中的移流项（由于Dv/Dt一般较小），可简化为

该方程为最为常见的声波控制方程。将上式作傅立叶变换得

2. 声学解析中的边界条件

粒子通过边界面的速度和声压的比定义为音响impendance。 在频域内可写为

但在时域内该值并不好确定。

3. 方程的变分形式

3.1 Helmholtz方程的变分形式

考虑下述带边界条件的Helmholtz方程

3.2 波动方程（4）的变分形式

4. 算例

4.1 标准benchmark测试问题AHLV100

该问题为Code_Aster导入的测试问题（1）。如下图所示，该问题在长梁的一端指定速度，另一端设置为吸收面 Z=ρ0c0 .计算在指定频率下(f=500Hz）下方程式（6）的解。由于该问题有理论解，数值计算结果可以与该理论解比较，用于验证计算软件实装的正确性。

Benchmark测试问题

下图为各种软件的计算结果。最下方一个是本软件计算结果，上面两个出自参考文献(2)

计算结果(声压的实数部分)

计算结果(声压的虚数部分)

4.2 居室内的声波传输

该算例通过求解方程（8）来模拟居室内的声音传播。在时间t=0时在下图的红点处施加一个Dirac函数格式的点声源，居室内各个壁面均为全反射。

有限元网格暨点声源位置

居室内音响: 声压的时间演化

本计算采用的是Newmark- β 法时间积分格式。上图显示的是室内声压随时间的变化。

5. 参考文献

1）

AHLV100 - Guide of wave at anechoic exitwww.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiakduFmaTyAhUFwosBHUFlCt8QFnoECAQQAQ&url=https%3A%2F%2Fwww.code-aster.org%2FV2%2Fdoc%2Fv13%2Fen%2Fman_v%2Fv8%2Fv8.22.100.pdf&usg=AOvVaw2sUyFOVu4fW-URza8rjrEq

2)

各種オープンソースFEMコードによる周波数領域音響伝搬解析精度のベンチマークテストwww.opencae.or.jp/wp-content/uploads/2015/06/Open