空力音響解析（FW-H方程式）
理论与物理
概述 — 什么是FW-H方程
老师，FW-H方程是什么？我听说它可以从CFD结果计算噪声，它的原理是怎样的？
问得好。FW-H方程是“Ffowcs Williams-Hawkings方程”的缩写，是John Ffowcs Williams和David Hawkings于1969年提出的一种声学类比方法。粗略地说，它是将Lighthill的声学类比扩展到运动物体表面的产物。
原来如此……具体是怎么使用的呢？不是直接用CFD计算所有的声音吗？
问题就在这里。声音是流体的微小压力波动，如果试图直接用CFD求解到远处，会受到数值衰减或反射的影响，精度难以保证。而且必须将网格扩展到观测点，计算成本会爆炸式增长。
这时FW-H就派上用场了。CFD只求解物体附近的流场，远场的声压则通过FW-H积分方程进行后处理计算。得益于这种两阶段方法，飞机喷流噪声、风扇噪声、汽车后视镜风噪等实际问题才得以在现实的计算成本下求解。
也就是说，它是作为CFD的“后处理”来使用的！那样的话计算域就可以保持紧凑了……但是，为什么通过积分就能知道远处的声压呢？
因为声音的本质是波的传播。只要知道物体表面（或其周围的虚拟面）上的压力和速度波动，就可以利用线性波动方程的格林函数计算任意远场点的声压。这正是FW-H积分的原理。
Lighthill的声学类比
刚才您提到“将Lighthill的类比进行了扩展”。能先讲讲Lighthill的类比吗？
Lighthill的声学类比（1952年）是气动声学的起点。通过变换纳维-斯托克斯方程，推导出以下波动方程：
其中 $\rho' = \rho - \rho_0$ 是密度波动，$c_0$ 是均匀场的声速，$T_{ij}$ 是Lighthill应力张量，定义如下：
$T_{ij}$ 代表什么？看右边有动量、压力、粘性应力等等……
关键在于，左边是“静止均匀介质中的波动方程”的形式。右边的 $T_{ij}$ 可以解释为等效声源。也就是说，“将实际的湍流全部塞进右边的声源项，而左边则视为静止流体中的声音传播”——这就是类比（Analogy）的本质。
实际的主导项是 $\rho u_i u_j$，这是喷流噪声的主要原因。它也被称为雷诺应力张量。在低马赫数下，其他项很小，因此可以近似为 $T_{ij} \approx \rho_0 u_i u_j$。
原来如此，因为将湍流“重新解读”为声源，所以才叫“类比”啊。但是这样一来，物体表面的效应没有包含进去吧？
很敏锐。确实如此，Lighthill最初的公式是针对自由空间中的湍流的。要处理旋转风扇叶片或飞行中飞机表面的效应，需要将边界条件纳入积分。这就是FW-H方程。
FW-H方程的公式化
终于到FW-H方程本体了。它是什么形式？
FW-H方程是使用广义函数（Heaviside函数、Dirac delta函数）将Lighthill方程扩展到运动物体表面 $f=0$ 的结果。最终，远场声压波动 $p'(\mathbf{x},t)$ 表示为三个项的和：
其中 $[\cdot]_{\mathrm{ret}}$ 表示在延迟时间（retarded time）下求值，$r$ 是声源点到观测点的距离，$M_r$ 是声源马赫数在观测方向的分量，$v_n$ 是表面的法向速度，$\ell_i = (p\delta_{ij} - \tau_{ij})n_j + \rho u_i(u_n - v_n)$ 是载荷矢量。
三个声源项的物理意义
哇，有三个项啊……它们分别对应什么具体的声音呢？
好，我们一个一个用具体例子说明。
1. 厚度声源（Thickness Noise）：物体运动时，其前方的空气被推开，后方的空气被吸引。直升机叶片旋转时，叶片的“厚度”周期性地推动空气的效果。表现为转速×叶片数的整数倍频率（BPF: 叶片通过频率）上的离散音调。
2. 载荷声源（Loading Noise）：作用在物体表面的压力（升力/阻力的波动）成为声源。风扇叶片与尾流干涉导致升力波动时，就会产生声音。对于汽车后视镜，镜面表面的压力波动会作为“呼呼”的风噪传入车内。
3. 四极子声源（Quadrupole Noise）：物体外部空间中的湍流本身发出的声音。喷气发动机排气喷流产生的“轰鸣”宽带噪声是典型例子。需要体积积分，计算成本高。但在低马赫数（约 $M < 0.3$）下，与厚度/载荷声源相比很小，因此在工程实践中常被省略。
明白了！也就是说，对于低速的汽车或风扇，只用面积分的两项就大致可以了，而对于像喷气发动机这样的高马赫数问题，四极子项也变得重要起来，对吧？
没错。实际上，在Ansys Fluent或STAR-CCM+的FW-H实现中，默认也只启用厚度+载荷的面积分，四极子项是作为选项处理的。
Lighthill的八次方定律与声源缩放
Lighthill推导的重要结果之一是低马赫数喷流的声功率与喷流速度$U$的八次方成正比，即八次方定律：
$$ W_a \propto \frac{\rho_0 U^8 D^2}{c_0^5} $$其中$D$是喷流直径。这与实验非常吻合，成为气动声学的基础。物理上，这意味着四极子声源的辐射效率依赖于马赫数的八次方。另一方面，表面上的载荷声源（偶极子）按$U^6$缩放，厚度声源（单极子）按$U^4$缩放。
FW-H方程的微分形式（广义函数形式）
FW-H方程原始的微分形式，使用Heaviside函数$H(f)$，可以写成：
$$ \Box^2 [p' H(f)] = \frac{\partial^2}{\partial x_i \partial x_j}[T_{ij} H(f)] - \frac{\partial}{\partial x_i}[\ell_i \delta(f)] + \frac{\partial}{\partial t}[Q_n \delta(f)] $$其中 $\Box^2 = \frac{1}{c_0^2}\frac{\partial^2}{\partial t^2} - \nabla^2$ 是达朗贝尔算子，$Q_n = \rho_0 v_n + \rho(u_n - v_n)$ 是质量通量项（厚度声源），$\delta(f)$ 是Dirac delta函数。用自由空间的格林函数求解这个微分方程，就得到了前述的积分形式。
各声源的指向性模式
| 声源类型 | 多极子阶数 | 速度依赖性 | 指向性 | 代表例 |
|---|---|---|---|---|
| 厚度（Thickness） | 单极子（Monopole） | $\propto U^4$ | 全向 | 叶片旋转、活塞 |
| 载荷（Loading） | 偶极子（Dipole） | $\propto U^6$ | $\cos\theta$ 模式 | 风扇噪声、后视镜 |
| 四极子（Quadrupole） | 四极子 | $\propto U^8$ | $\cos^2\theta$ 模式 | 喷流噪声 |
可渗透面（Permeable Surface）公式化
我有时会听到“可渗透面”，是指不使用物体表面，而是使用虚拟的面吗？
这是一个非常重要的工程技巧。原始的FW-H方程将积分面取在物体表面（固体面），但di Francescantonio（1997）提出的可渗透面公式化，则使用包围物体的虚拟闭合曲面（permeable surface）作为积分面。
有两个优点。首先，可以避免四极子的体积积分。只要将可渗透面取得足够大以包含声源区域，仅通过面积分就能捕捉到四极子的贡献。其次，非线性流动效应（包含激波的高马赫数流动）也能作为面上的数据被纳入。
那么使用可渗透面就万事大吉了吗？没有缺点吗？
遗憾的是，穿过可渗透面的涡旋可能会产生“虚假声源（spurious noise）”。特别是当喷流穿过可渗透面的下游边缘时，面内和面外计算的声音无法完全抵消，会产生非物理的低频噪声。这被称为端盖问题（end-cap problem）。面的布置需要工匠般的技巧。
孕育FW-H方程的两位天才
John Ffowcs Williams（威尔士人，剑桥大学教授）被称为“航空声学之父”，因协和式飞机喷流噪声研究而闻名。合作研究者David Hawkings当时是研究生，1969年发表在Philosophical Transactions期刊上的论文，成为此后50多年气动声学分析的基础。另一方面，奠定理论基础的Sir James Lighthill是流体力学领域的巨人，其1952年关于声学类比的论文引用次数超过10,000次。有趣的是，Lighthill在1952年提出的理论，在当时单台喷气发动机噪声测试需要数亿日元的时代，是一个“通过计算预测噪声”的革命性想法。
数值解法与实现
时域积分法（Farassat 1A）
实际在计算机上求解FW-H方程要怎么做呢？直接数值计算那个积分吗……？
工程实践中最常用的是Farassat 1A公式化。这是NASA的Farassat在1980年代推导出的FW-H方程在时域中的闭合形式解。
厚度声源的贡献为：
载荷声源的贡献为：
延迟时间（retarded time）的计算看起来很难。要计算从每个面元到观测者的声音到达时间……
这正是FW-H代码实现的关键所在。对于每个面元，通过下式求解延迟时间$\tau$：
这个隐式方程通过牛顿法迭代求解。在均匀流中可解析求解，但存在平均流效应时则需要修正。许多商业软件采用假设均匀流的快速算法。
频域法
除了时域方法，还有其他方法吗？
也有在频域求解FW-H方程的方法。将面上的压力/速度数据通过FFT分解为频率分量，然后对每个频率使用亥姆霍兹方程的格林函数进行积分。
优点在于当只需要计算特定频带时效率高。缺点是需要沿时间方向保存大量的非定常数据。在工程实践中，它适用于旋转机械的纯音噪声（BPF分量）分析。
与CFD的耦合流程
CFD和FW-H的耦合，具体是如何进行数据交换的呢？
典型的流程是这样的：
- 执行CFD非定常计算，在每个时间步记录FW-H面上的压力$p$、密度$\rho$、速度$\mathbf{u}$
- 获得足够的统计量后（通常是几十到几百个流动通过时间），将记录的数据传递给FW-H求解器
- FW-H求解器计算各观测点的声压时间历程$p'(t)$
- 对声压时间历程进行FFT，计算频谱（SPL vs. 频率）或OASPL（总声压级）
在Ansys Fluent或STAR-CCM+中，也可以使用与CFD计算同时执行FW-H积分的“在线（on-the-fly）”模式。这可以节省内存，但缺点是无法在事后更改观测点位置。
湍流模型的选择
CFD侧的湍流模型应该用什么？k-ε模型之类的可以吗？
RANS（k-ε, k-ω SST等）不能用于气动声学分析。这是非常重要的一点。RANS求解的是时间平均的流场，因此无法解析作为声源的瞬时压力波动。
气动声学需要的是：
- LES（大涡模拟）：直接解析大尺度涡，小尺度通过SGS模型近似。可信度最高，但计算成本也最高
- DES/DDES（分离涡模拟）：壁面附近用RANS，分离区域用LES的混合方法。具有较好的工程实用性平衡
- SAS（尺度自适应模拟）：Menter-Egorov模型。基于URANS但能实现类似LES的分辨率
例如，对于汽车后视镜风噪，DDES + FW-H已成为工程实践中的标准。对于风扇噪声，推荐使用LES；对于喷流噪声，LES更理想，但由于计算成本的关系，DDES也有使用。
网格与时间步长的要求
要计算声音的话，网格要求也比普通CFD更严格吗？
没错。气动声学所需分辨率的参考标准：
| 参数 | 推荐值 | 备注 |
|---|---|---|