Ffowcs Williams-Hawkings 模型

Ffowcs Williams-Hawkings (FW-H) 声学积分公式是预测高空飞机噪声等远场噪声的首选策略。此模型可计算远场声信号，这些信号由 CFD 求解得出的近场流数据扩展得到。目的是预测每个接收器位置处的小振幅声压波动。

索引

流向
位于声音体表面（例如，圆柱表面）的源上的不渗透表面
封闭音源的渗透表面
计算域边界/直接噪声模拟近场求解的极限
远场中的接收器
四极噪声可视化

适用于固定运动的 FW-H 声学模型可模拟具有固定源（例如飞机起落架或翼面尾缘）的风洞气动声学测量。FW-H 声学模型只用于预测自由空间中的声音传播。它不包括声音反射、折射或材料改性等效果。

Ffowcs Williams-Hawkings 公式基于 Farassat 的公式 1A（[897]、[905]）和 Dunn Farassat Padula 1A ([915])，二者均适用于 FW-H 的非对流形式。因此，无需十分关注数值精度。通过 FW-H 方法（一种空间与时间二阶离散格式）以及明智选择的网格分辨率和时间步长，可准确预测远场噪声。

Ffowcs Williams-Hawkings 模型

Ffowcs Williams 和 Hawkings 的研究 [909] 将 Lighthill [918] 和 Curle 的研究 [901] 扩展为由任意运动的表面生成的空气动力声公式。

与 Lighthill 方程相似，他们以不均匀波方程的形式改写了流体运动方程。定义 FW-H 方程的一个重要方面是将特定的流体变量编写为广义函数。

Ffowcs Williams 和 Hawkings（请参见 [909]）获得的结果是针对不渗透表面正式提出的，但是他们明白其公式中的数学表面不必与固体重合。

不渗透和渗透表面均适用于 Ffowcs Williams-Hawkings 声学类比模型，相关公式如下。

Ffowcs Williams-Hawkings 公式

Ffowcs Williams-Hawkings (FW-H) 公式基于 Farassat 的公式 1A 和 Dunn Farassat Padula 公式 1A。

Farassat 的公式 1A（请参见 [896]）是默认且首选的公式（用于通用亚音速源区域的 FW-H 的非对流形式），其中包括用于整个域和周期域中的瞬态旋转运动的不渗透公式。
保持 Dunn Farassat Padula 1A（请参见 [907]）的目的是确保亚音速旋转截面（表示整个域和周期域中的瞬态旋转运动的不渗透公式）的向后兼容性（FW-H 的非对流形式）。

FW-H 方程是将连续性方程和动量方程精确地重新整理为不均匀波方程形式的结果。即使在积分表面位于非线性流体区域中的情况下，FW-H 方程的结果同样精确。它根据自由空间格林函数来计算观察器位置 $x$ 处的声压。对于区域或一组表面中的流体辐射到静止介质中的压力，FW-H 方程如下：

p′ (x, t) = {p′}_{T} (x, t) + {p′}_{L} (x, t) + {p′}_{Q} (x, t)

(4884)

单极项如下：

{p′}_{T} (x, t) = \frac{1}{4 π} ((\frac{\partial}{\partial t}) \int_{S}^{} {[\frac{Q}{(r (1 - M_{r}))}]}_{r e t} d S)

(4885)

偶极项如下：

{p′}_{L} (x, t) = \frac{1}{4 π} ((- \frac{\partial}{\partial x_{i}}) \int_{S}^{} {[\frac{L_{i}}{(r (1 - M_{r}))}]}_{r e t} d S)

(4886)

四极项如下：

{p′}_{Q} (x, t) = \frac{1}{4 π} ((\frac{\partial^{2}}{(\partial x_{i}) (\partial x_{j})}) \int_{V}^{} {[\frac{T_{i j}}{(r (1 - M_{r}))}]}_{r e t} d V)

(4887)

且：

Q = ρ_{0} U_{i} n_{i}

(4888)

U_{i} = (1 - \frac{ρ}{ρ_{0}}) v_{i} + \frac{ρ u_{i}}{ρ_{0}}

(4889)

L_{i} = P_{i j} n_{i} + ρ u_{i} (u_{n} - v_{n})

(4890)

P_{i j} = (p - p_{0}) δ_{i j} - σ_{i j}

(4891)

T_{i j} = ρ u_{i} u_{j} + δ_{i j} [(p - p_{0}) - {c_{0}}^{2} (ρ - ρ_{0})] - σ_{i j}

(4892)

其中：

$u_{i}$ 表示 $i$ 方向上的流体速度分量。
$u_{n}$ 为与表面垂直的流体速度分量。
$v_{i}$ 表示 $i$ 方向上的表面速度分量。
$v_{n}$ 为与表面垂直的表面速度分量。
$n_{i}$ 为表面法向矢量。
$p_{0}$ 为大气/远场压力，静止介质中的流体压力。
$σ_{i j}$ 为粘性应力张量。
$ρ_{0}$ 为远场密度。
$P_{i j}$ 为压缩应力张量。
$T_{i j}$ 为 Lighthill 应力张量。

下标 $r e t$ 表示阻碍时间，即发射的时间。

根据 Eqn. (4885) 和 Eqn. (4886) 得到的空间导数将转换为时间导数。随后，观察器位置处的时间导数将移入积分中。

Eqn. (4894) 使用提前时间公式来描述单极源项。Eqn. (4895) 使用提前时间公式来描述偶极源项。

自动信号修剪

在导出计算的 FW-H 声压之前，将自动修剪 FW-H 接收器数据，即从时间信号的两端剪切数据，因为这些间隔中的声源不完整。

根据固定 FW-H 表面，信号从 FW-H 表面移到接收器的最小和最大行程时间 ${T D}_{MIN}$ 和 ${T D}_{MAX}$ 为：

{T D}_{MIN} = \frac{r_{min}}{c_{0}}, {T D}_{MAX} = \frac{r_{max}}{c_{0}}

其中：

$r_{\min}$ 为 FW-H 源表面与接收器位置之间的最短距离。
$r_{\max}$ 为 FW-H 源表面与接收器位置之间的最长距离。
$c_{0}$ 为声速。

在接收器位置处预测的远场噪声从时间 $t_{0}$ 到 $t_{n}$ 进行表格化：

\begin{matrix} t_{0} = T_{0} + {T D}_{MIN} \\ t_{n} = T_{n} + {T D}_{MAX} \end{matrix}

其中：

$t_{0}$ 为修剪之前，数据集开始在接收器处累积的时间。
$t_{n}$ 为修剪之前，数据集停止在接收器处累积的时间。
$T_{0}$ 为实时 FW-H 求解器的开始时间属性，或是用后 FW-H 求解器的 .simh 文件格式记录的第一个物理时间。
$T_{n}$ 为实时 FW-H 求解器的模拟结束时间，或是用后 FW-H 求解器的 .simh 文件格式记录的最后一个物理时间。

接收器不会接收来自 FW-H 源表面 ( $r_{\max}$ ) 上最远点的声压信号，直到开始接收时间变为：

t_{b} = T_{0} + {T D}_{MAX}

从时间 $t_{0}$ 到 $t_{b}$ ，接收器寄存器上累积的声音信号不包含来自整个 FW-H 表面区域的贡献，因此在该时间内收到的声压数据不完整。从结束时间 $t_{e}$ 到 $t_{n}$ 的期间内会出现相同的修剪过程：

\begin{matrix} t_{e} = T_{n} + {T D}_{MIN} \end{matrix}

修剪时将移除不完整的端点上的信号，持续时间为从 $t_{0}$ 到 $t_{b}$ （即持续了 ${T D}_{MAX} - {T D}_{MIN}$ 的时间间隔），然后再持续 $t_{e}$ 到 $t_{n}$ （即另一个 ${T D}_{MAX} - {T D}_{MIN}$ 间隔），因此总持续时间为 $2 \cdot ({T D}_{MAX} - {T D}_{MIN})$ 。修剪后的信号跨越的时间为从修剪数据开始时的 $t_{b}$ 到修剪数据结束时的 $t_{e}$ 。

对于旋转 FW-H 表面，将根据旋转轴周围的多个接收器位置计算 ${T D}_{MIN}$ 和 ${T D}_{MAX}$ 。这些位置通过相对于 FW-H 表面初始位置以 1 度增量将接收器旋转 360 度来确定。对于多个接收器，接收器数据将使用所有接收器中的最小值 ${T D}_{MIN}$ 和最大值 ${T D}_{MAX}$ 进行调整。此方法可确保与旋转轴距离相等的接收器的信号延迟一致。

对于多个接收器，分别计算每个接收器的 ${T D}_{MIN}$ 和 ${T D}_{MAX}$ 。

FW-H 表面项

当积分表面与固体表面重合时，单极项 Eqn. (4885)、偶极项 Eqn. (4886) 和四极项Eqn. (4887) 称为：

厚度表面项 ${p′}_{T} (x, t)$ ，由体通过时流体的位移导致。Eqn. (4894) 中定义了用于常规流体的该项；Eqn. (4896) 中定义了具有刚体运动或移动参考坐标系的流体的该项。
负载表面项 ${p′}_{L} (x, t)$ ，由体表面上力分布的瞬态变化导致。Eqn. (4895) 中定义了用于常规流体的该项；Eqn. (4897) 中定义了具有刚体运动或移动参考坐标系的流体的该项。
总表面项 ${p′}_{S} (x, t)$ ，由对厚度表面项和负载表面项求和获得。该项定义如下：

{p′}_{S} (x, t) = {p′}_{T} (x, t) + {p′}_{L} (x, t)

(4893)

Farrassat 的公式 1A

此公式适用于常规亚音速源区域以及常规远场噪声预测。请参见 [896] 和 [897]。

{p′}_{T} (x, t) = \frac{1}{4 π} (\int_{(f = 0)}^{} {[\frac{ρ_{0} ({\dot{U}}_{n} + U_{\bar{n}})}{r {(1 - M_{r})}^{2}}]}_{r e t} d S + \int_{(f = 0)}^{} {[\frac{ρ_{0} U_{n} [r {\dot{M}}_{r} + a_{0} (M_{r} - M^{2})]}{r^{2} {(1 - M_{r})}^{3}}]}_{r e t} d S)

(4894)

{p′}_{L} (x, t) = \frac{1}{4 π} (\frac{1}{a_{0}} \cdot \int_{(f = 0)}^{} {[\frac{{\dot{L}}_{r}}{r {(1 - M_{r})}^{2}}]}_{r e t} d S + \int_{(f = 0)}^{} {[\frac{(L_{r} - L_{M})}{r^{2} {(1 - M_{r})}^{2}}]}_{r e t} d S + \frac{1}{a_{0}} \cdot \int_{(f = 0)}^{} {[\frac{L_{r} [r {\dot{M}}_{r} + a_{0} (M_{r} - M^{2})]}{r^{2} {(1 - M_{r})}^{3}}]}_{r e t} d S)

(4895)

其中：

$L_{i} = P_{i j} \cdot n_{i} + ρ \cdot u_{i} (u_{n} - v_{n})$ 其中 $P_{i j} = (p - p_{0}) δ_{i j}$
$M_{i} = v_{i} / a_{0}$
$r = x_{observer} - y_{face}$

其中：

$f = 0$ 表示用于将外部流问题 $f > 0$ 嵌入无界限空间内的数学表面。
$f = 0$ 表示发射表面，并设置为与体、不渗透表面或渗透表面重合。

如果数据表面与固体表面重合，则流体的法向速度将与表面的法向速度相同： $u_{n} = v_{n}$

在这种情况下，Eqn. (4894) 和 Eqn. (4895) 将对应于不渗透 FW-H 表面类型并且将消除其中一些项。

Dunn-Farrassat-Padula 公式 1A

此公式基于 Farassat 的公式 1A 的修正，后者在 ANOPP（飞机噪声预测程序）和 DFP-ATP（Dunn-Farassat-Padula 高级涡轮螺旋桨发动机预测）都用于亚音速旋转叶片，以便根据 Hubbard 的研究 [915] 计算声压。

对于渗透 FW-H 表面类型，将对内部交界面边界计算 FW-H 表面积分，从原位交界面开始。

{p′}_{T} (x, t) = \frac{1}{4 π} {(\int_{(f = 0)}^{} [\frac{ρ_{0} v_{n} [(r {\dot{M}}_{r}) + a_{0} (M_{r} - M^{2})]}{r^{2} {(1 - M_{r})}^{3}}] d S)}_{r e t}

(4896)

{p′}_{L} (x, t) = \frac{1}{4 π} {(\frac{1}{a_{0}} \cdot \int_{(f = 0)}^{} [\frac{{\dot{L}}_{i} r_{i}}{r {(1 - M_{r})}^{2}}] d S + \int_{(f = 0)}^{} [\frac{(L_{r} - L_{i} M_{i})}{r^{2} {(1 - M_{r})}^{2}}] d S + \frac{1}{a_{0}} \cdot \int_{(f = 0)}^{} [\frac{L_{r} [(r {\dot{M}}_{r}) + a_{0} (M_{r} - M^{2})]}{r^{2} {(1 - M_{r})}^{3}}] d S)}_{r e t}

(4897)

和

{\dot{M}}_{r} = {\dot{M}}_{i} r_{i}

(4898)

其中：

$L_{i}$ 为叶片负载矢量。 $L_{i} = P_{i j} \cdot n_{i} + ρ \cdot u_{i} (u_{n} - v_{n})$ 其中 $P_{i j} = (p - p_{0}) δ_{i j}$
${\dot{L}}_{i}$ 为叶片负载矢量对源时间的导数。
$v_{i}$ 为叶片表面相对于静止流体的局部速度。
$v_{n}$ 为表面法向速度。
$r_{i}$ 为从源点到观察器的距离。
$M_{r}$ 为源朝向观察器的马赫数。

下标 $r e t$ 表示阻碍时间，即发射的时间。

考虑时间滞后

计算 FW-H 积分时的主要任务是如何考虑发射和接收时间之间的时间滞后。

提前时间算法可在时间上提前看到观察器察觉当前生成的声波的时间：

该方法从一系列发射时间开始（为方便起见，以下称作流动时间）。
源强度在给定发射时间内在所有源单元（积分表面的面）处计算（厚度表面噪声和负载表面噪声）。
源的贡献将插值到远场时间域中，用于构建声信号。

观察器察觉的总声压包含所有源单元的贡献。通过累积时隙中的到达信号，可获得接收器处的声压。在相同的源时间范围内，通过对 FW-H 表面的每个源单元的声信号求和，可获得观察器整体声信号。

FW-H 体积项

四极噪声是考虑了流体中的非线性的源体积分布。

这些非线性分为两种类型 [919]。首先，局部声速不是恒定的，会随颗粒加速而变化。其次，体（例如叶片）附近的有限颗粒速度会影响声音传播的速度。当流体中存在强剪切层时或当马赫数增加时，四极项不可忽略。

Farassat 和 Brentner 的研究 [896] 表明，四极 ${p′}_{Q} (x, t)$ 的噪声贡献可表达为“消失球”公式。使用此公式，由 Eqn. (4887) 得到的空间导数将转换为时间导数：

{p′}_{Q} (x, t) = \frac{1}{4 π} ((\frac{1}{c}) (\frac{\partial^{2}}{{\partial t}^{2}}) \int_{- \infty}^{t} [\int_{(f > 0)}^{} \frac{T_{r r}}{r} d Ω] d τ + (\frac{\partial}{\partial t}) \int_{- \infty}^{t} [\int_{(f > 0)}^{} \frac{{3 T}_{r r} - T_{i i}}{r^{2}} d Ω] d τ) + (c \int_{- \infty}^{t} [\int_{(f > 0)}^{} \frac{{3 T}_{r r} - T_{i i}}{r^{3}} d Ω] d τ)

(4899)

其中：

$T_{r r}$ 为 $T_{i j}$ 的双倍收缩。
$r_{i}$ 和 $r_{i}$ 为辐射方向上单位矢量的分量。
$T_{i j}$ 为 Lighthill 应力张量。 $T_{i j} = ρ u_{i} u_{j} + δ_{i j} [(p - p_{0}) - {c_{0}}^{2} (ρ - ρ_{0})] - σ_{i j}$ ，且 $σ_{i j}$ 为粘性应力张量。

Eqn. (4899) 从消失球公式转换为使用以下四个方程的提前时间公式。在这些方程中，观察器处的时间导数将移入积分，以防止积分出现数值时间差异化。根据 [898] 获得的“源时间主导”算法用于估算 FW-H 方程的 ${p′}_{Q} (x, t)$ 体积项，如下所示：

{p′}_{Q} (x, t) = \frac{1}{4 π} (\int_{(f > 0)}^{} {[\frac{K_{1}}{c^{2} r} + \frac{K_{2}}{c r^{2}} + \frac{K_{3}}{r^{3}}]}_{r e t} d V)

(4900)

K_{1} = K_{11} + K_{12} + K_{13} = [\frac{{\ddot{T}}_{r r}}{{(1 - M_{r})}^{3}}] + [\frac{{{\ddot{M}}_{r} T}_{r r} + 3 {\dot{M}}_{r}_{\bar{T}}}{{(1 - M_{r})}^{4}}] + [\frac{3 {{\dot{M}}_{r}^{2} T}_{r r}}{{(1 - M_{r})}^{5}}]

(4901)

K_{2} = K_{21} + K_{22} + K_{23} + K_{24} = [\frac{- {\dot{T}}_{i i}}{{(1 - M_{r})}^{2}}] + [\frac{4 {\dot{T}}_{M_{r}} + 2 T_{{\dot{M}}_{r}} + {\dot{M}}_{r} T_{i i}}{{(1 - M_{r})}^{3}}] + [\frac{3 [(1 - M^{2}) {\dot{T}}_{r r} - 2 {\dot{M}}_{r} T_{M r} - M_{i} {\dot{M}}_{i} T_{r r}]}{{(1 - M_{r})}^{4}}] + [\frac{6 {\dot{M}}_{r} (1 - M^{2}) T_{r r}}{{(1 - M_{r})}^{5}}]

(4902)

K_{3} = K_{31} + K_{32} + K_{33} = [\frac{2 T_{M M} - (1 - M^{2}) T_{i i}}{{(1 - M_{r})}^{3}}] - [\frac{6 (1 - M^{2}) T_{M r}}{{(1 - M_{r})}^{4}}] + [\frac{3 {(1 - M^{2})}^{2} T_{r r}}{{(1 - M_{r})}^{5}}]

(4903)

其中：

$T_{r r} = T_{i j} r_{i} r_{j}$	$T_{M M} = T_{i j} M_{i} M_{j}$	$T_{M r} = T_{i j} M_{i} r_{j}$	$T_{\dot{M r}} = T_{i j} {\dot{M}}_{i} r_{j}$
${\dot{T}}_{M r} = {\dot{T}}_{i j} M_{i} r_{j}$	${\dot{T}}_{r r} = {\dot{T}}_{i j} r_{i} r_{j}$	${\ddot{T}}_{r r} = {\ddot{T}}_{i j} r_{i} r_{j}$	$M_{i} = \frac{v_{i}}{c}$
$M_{r} = M_{i} r_{i}$	${\dot{M}}_{r} = \dot{M} r_{i}$	${\ddot{M}}_{r} = {\ddot{M}}_{i} r_{i}$	${\ddot{M}}_{r r} = \ddot{M} r_{i} r_{i}$
$r_{i} = \frac{x_{i} - y_{i}}{r}$

其中：

$r_{i}$ 表示辐射方向上的单位矢量。
变量上方的点表示相对于该变量的源时间的导数。