声波模型

Ewert 和 Schröder 提出了一系列声扰动方程 (APE)，用于模拟流体引发的声场。这些方程描述通过可压缩或不可压缩问题的非稳态模拟确定的噪声源项。正如 Lighthill 根据基本纳维-斯托克斯方程推导得出声学类比一样，声波方程也可以从声扰动方程中衍生。

此模型计算整个域中的声压历史，包括对流和折射效应。

请参见 [887] 和 [901]。

根据扰动变量 $(p^{'}, u^{a})$ ，完整的声扰动方程组如下

图 1. EQUATION_DISPLAY

\frac{\partial p^{'}}{\partial t} + c^{2} \nabla\cdot (\bar{ρ} u^{a} + \bar{v} \frac{p^{'}}{c^{2}}) \approx 0

(4711)

图 2. EQUATION_DISPLAY

\frac{\partial u^{a}}{\partial t} + \nabla (\bar{v} \cdot u^{a}) + \nabla (\frac{p^{'}}{\bar{ρ}}) \approx \nabla Φ_{p}

(4712)

其中：

$p^{'}$ 为扰动压力
$u^{a}$ 为无旋扰动速度
$\bar{ρ}$ 为时间平均密度
$\bar{v}$ 为时间平均速度
$c$ 为声速
$Φ_{p}$ 为噪声源函数

使用 Ewert 和 Schroder 提出的关系 $p^{'} = \bar{ρ} Φ_{p} + p^{a}$ （其中， $p^{a}$ 为声压）并求 Eqn. (4711) 的物质导数 $(\frac{\partial}{\partial t} + \bar{v} \cdot \nabla)$ 和 Eqn. (4712) 的发散，即可获得以 $p^{a}$ 和 $u^{a}$ 为因变量的两个方程。对于不可压缩流（即， $\nabla\cdot \bar{v} = 0$ ）合并这两个方程，即可根据不可压缩流的声压和噪声源得出一个方程：

图 3. EQUATION_DISPLAY

\frac{1}{c^{2}} \frac{\partial^{2} p^{a}}{\partial t^{2}} + \frac{2 \bar{v}}{c^{2}} \cdot \nabla \frac{\partial p^{a}}{\partial t} + \frac{(\bar{v} \cdot \nabla)}{c^{2}} (\nabla \cdot \bar{v} p^{a}) - \nabla^{2} p^{a} = - (\frac{\bar{ρ}}{c^{2}} \frac{\partial^{2} Φ_{p}}{\partial t^{2}} + \frac{2 (\bar{v} \cdot \nabla)}{c^{2}} \bar{ρ} \frac{\partial Φ_{p}}{\partial t} + \frac{(\bar{v} \cdot \nabla)}{c^{2}} (\nabla \cdot \bar{v} Φ_{p}))

(4713)

对于主要由对流效应引起的波，通过假设不可压缩流的噪声源 [887]， $\bar{ρ} \frac{\partial Φ_{p}}{\partial t} \approx \frac{\partial P^{'}}{\partial t}$ 和 $\nabla Φ_{p} \approx \frac{1}{\bar{ρ}} \nabla P^{'}$ （其中 $P^{'}$ 为流体动力压力），即可从 Eqn. (4713) 衍生出一个方程。在该方程中，需要物理阻尼机制来消除源自网格粗化转换的虚假波。根据粘性介质中的波传播概念，使用物理阻尼项的最终方程是：

图 4. EQUATION_DISPLAY

\frac{1}{c^{2}} \frac{\partial^{2} p^{a}}{\partial t^{2}} + \frac{2 \bar{v}}{c^{2}} \cdot \nabla \frac{\partial p^{a}}{\partial t} + \frac{(\bar{v} \cdot \nabla)}{c^{2}} (\nabla \cdot \bar{v} p^{a}) - \nabla^{2} (p^{a} + τ \frac{\partial p^{a}}{\partial t}) = - (\frac{1}{c^{2}} \frac{\partial^{2} P^{'}}{\partial t^{2}} + \frac{2 (\bar{v} \cdot \nabla)}{c^{2}} \frac{\partial P^{'}}{\partial t} + \frac{(\bar{v} \cdot \nabla)}{c^{2}} (\nabla \cdot \bar{v} P^{'}))

(4714)

其中，在物理阻尼项中，

τ = χ \frac{Δ t}{π λ}

。在此项中，

$χ$ 为阻尼系数（= 0 表示无阻尼，= 1 表示最大阻尼）
$Δ t$ 为时间步
$λ = \frac{c Δ t}{Δ x}$ 为局部 CFL 数（其中， $Δ x$ 为网格间距）

可使用此阻尼系数在求解不当的区域（粗糙网格）中阻尼声波，从而抑制不需要的振荡，避免影响相关区域内部的预测。

没有平均流时，Eqn. (4714) 简化为

图 5. EQUATION_DISPLAY

\frac{1}{c^{2}} \frac{\partial^{2} p^{a}}{\partial t^{2}} - \nabla^{2} (p^{a} + τ \frac{\partial p^{a}}{\partial t}) = - \frac{1}{c^{2}} \frac{\partial^{2} P^{'}}{\partial t^{2}}

(4715)

此方程类似于 Lighthill 方程，不同之处在于它基于不可压缩流的声扰动方程中的噪声源。

低通时间滤波器

声波模型中的噪声源主要源自流体压力波动对时间的二阶导数。由于网格尺寸转换，这些噪声源会包括高频虚假数值信号，进而生成高频噪声。合适的低通滤波器可抑制这些不需要的噪声。

对于以下滤波器，滤波器越高阶，它提供的信号阻尼越多，但同时引入的相位误差（即，波角与预期值的偏移量）也越多。可以使用相位误差度量不同频率的相位差引入的信号失真。

声波模型直接将低通时间滤波器应用于 APE 源。滤波 APE 源为：

- \tilde{\frac{1}{c^{2}} \frac{\partial^{2}}{\partial t^{2}} P^{'}} = \frac{1}{c^{2}} \frac{\partial^{2} p^{a}}{\partial t^{2}} - \nabla \cdot \nabla (p^{a} + τ \frac{\partial p^{a}}{\partial t})

一阶 IIR 低通时间滤波器

\frac{d^{2} P^{'}}{d t^{2}}

在时间步

n

的滤波形式为：

图 6. EQUATION_DISPLAY

{\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n} = α {(\frac{d^{2} P^{'}}{d t^{2}})}^{n} + (1 - α) {\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n - 1}

(4716)

其中：

α = \frac{2 π Δ t f_{c u t}}{2 π Δ t f_{c u t} + 1}

$\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}$ 在给定时间步的滤波值取决于上一个时间步的滤波值。

二阶 Butterworth 低通时间滤波器

\frac{d^{2} P^{'}}{d t^{2}}

在时间步

n

的滤波形式为：

图 7. EQUATION_DISPLAY

{\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n} = [\frac{n_{0}}{d_{0}} {(\frac{d^{2} P^{'}}{d t^{2}})}^{n} + \frac{n_{1}}{d_{0}} {(\frac{d^{2} P^{'}}{d t^{2}})}^{n - 1} + \frac{n_{2}}{d_{0}} {(\frac{d^{2} P^{'}}{d t^{2}})}^{n - 2}] + [- \frac{d_{1}}{d_{0}} {\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n - 1} - \frac{d_{2}}{d_{0}} {\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n - 2}]

(4717)

其中：

$n_{0} = 1$ ， $n_{1} = 2$ ，而 $n_{2} = 1$
$d_{0} = e^{2} + \sqrt{2} e + 1$
$d_{1} = - 2 (e^{2} - 1)$
$d_{2} = e^{2} - \sqrt{2} e + 1$
$e = \frac{1}{\tan (\frac{2 π Δ t f_{c u t}}{2})}$

三阶 Butterworth 低通时间滤波器

\frac{d^{2} P^{'}}{d t^{2}}

在时间步

n

的滤波形式为：

图 8. EQUATION_DISPLAY

{\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n} = [\frac{n_{0}}{d_{0}} {(\frac{d^{2} P^{'}}{d t^{2}})}^{n} + \frac{n_{1}}{d_{0}} {(\frac{d^{2} P^{'}}{d t^{2}})}^{n - 1} + \frac{n_{2}}{d_{0}} {(\frac{d^{2} P^{'}}{d t^{2}})}^{n - 2} + \frac{n_{3}}{d_{0}} {(\frac{d^{2} P^{'}}{d t^{2}})}^{n - 3}] + [- \frac{d_{1}}{d_{0}} {\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n - 1} - \frac{d_{2}}{d_{0}} {\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n - 2} - \frac{d_{3}}{d_{0}} {\tilde{(\frac{d^{2} P^{'}}{d t^{2}})}}^{n - 3}]

(4718)

其中：

$n_{0} = 1$ ， $n_{1} = 3$ ， $n_{2} = 3$ ，而 $n_{3} = 1$
$d_{0} = e^{3} + 2 e^{2} + 2 e + 1$
$d_{1} = - 3 e^{3} - 2 e^{2} + 2 e + 3$
$d_{2} = 3 e^{3} - 2 e^{2} - 2 e + 3$
$d_{3} = - e^{3} + 2 e^{2} - 2 e + 1$
$e = \frac{1}{\tan (\frac{2 π Δ t f_{c u t}}{2})}$

非反射边界条件

借助声波模型中的非反射边界条件，声波可在不产生任何虚假反射的情况下穿过边界离开域。假设边界处没有源并对动量方程 Eqn. (4712) 应用运算符 " $\nabla\cdot$ "：

图 9. EQUATION_DISPLAY

\begin{array}{l} \nabla\cdot \frac{\partial u^{a}}{\partial t} + (\bar{v} \cdot \nabla) \nabla\cdot u^{a} + \frac{1}{\bar{ρ}} \nabla\cdot \nabla p^{a} = 0 \\ \nabla\cdot u^{a} = - \frac{1}{\bar{ρ} c^{2}} (\frac{\partial p^{a}}{\partial t} + \bar{v} \cdot \nabla p^{a}) \end{array}

(4719)

声压与声速的法向分量相关：

图 10. EQUATION_DISPLAY

p^{a} = \bar{ρ} c u^{a} \cdot n

(4720)

通过求 Eqn. (4720) 的时间导数并使用 Eqn. (4719)，声压梯度变为：

图 11. EQUATION_DISPLAY

\nabla p^{a} \cdot A_{f} = - \frac{(1 - \bar{v} \cdot n / c) \partial p^{a}}{(1 - {\bar{v}}^{2} / c^{2}) c \partial t} ‖ A_{f} ‖

(4721)

其中， $A_{f}$ 是面法向面积，而 $‖ A_{f} ‖$ 是其幅值。

纽马克 α 法

纽马克 Alpha 法从声压时间导数的以下表达式开始：

M {(\frac{\partial^{2} p^{a}}{\partial t^{2}})}^{n + 1} + C {(\frac{\partial^{} p^{a}}{\partial t^{}})}^{n + 1} + K {(p^{a})}^{n + 1} = F^{n + 1}

其中，M、C 和 K 为算子：

$M = \frac{1}{c^{2}} \int_{V} d V$ 为质量算子。
$C = \frac{1}{c^{2}} \int_{s} \bar{v} d s$ 为阻尼算子。
$K = \int_{s} \nabla\cdot d s + \frac{1}{c^{2}} (\bar{v} \cdot \nabla) \bar{v} \cdot d s$ 为刚度算子。
$F$ 为噪声源项。
$n$ 为上一个时间步。
$n + 1$ 为要计算求解的时间步。

然后，该方法离散化上述方程并引入 $α$ 参数，得出以下修正控制方程：

图 12. EQUATION_DISPLAY

M {(\frac{\partial^{2} p^{a}}{\partial t^{2}})}^{n + 1} + (1 + α) C {(\frac{\partial^{} p^{a}}{\partial t^{}})}^{n + 1} + (1 + α) K {(p^{a})}^{n + 1} = {(1 + α) F}^{n + 1} - α F^{n} + α C {(\frac{\partial^{} p^{a}}{\partial t^{}})}^{n} + α K {(p^{a})}^{n}

(4722)

这些算子和外力采用矩阵的形式，其中：

$M$ 为质量矩阵
$C$ 为阻尼矩阵
$K$ 为刚度矩阵
$F$ 为外力矩阵
$α$ 为加权参数， $- 0.33 \leq α \leq 0$ 。

后 FW-H 模型的 APE 场

声波模型可与 Ffowcs Williams-Hawkings 模型相结合，用来计算远场案例的以下场：

声密度 $ρ^{a}$ ：
图 13. EQUATION_DISPLAY

$ρ^{a} = \frac{p^{a}}{c^{2}} + ρ_{f}$

(4723)
声速 $u^{a}$ ：
$u^{a}$ 在第 n 个时间步中 ( $u_{n}^{a}$ ) 可通过上一个时间步的声速和声压梯度计算得出。此项通过声扰动方程的动量方程的二阶向后差分得出。
图 14. EQUATION_DISPLAY

$u_{n}^{a} = \frac{1}{3} [4 u_{n - 1}^{a} - u_{n - 2}^{a} - \frac{2 Δ t \nabla p^{a}}{ρ_{0}}]$

(4724)
APE 总密度 $Σ ρ$ ：
图 15. EQUATION_DISPLAY

$Σ ρ = \frac{P^{a}}{c^{2}} + ρ_{f}$

(4725)
APE 总速度 $Σ v$ ：
图 16. EQUATION_DISPLAY

$Σ v = v^{a} + v_{f}$

(4726)
APE 总压力 $Σ P$ ：
图 17. EQUATION_DISPLAY

$Σ P = P^{a} + P_{f}$

(4727)

其中：

$c$ 为声速。
$ρ_{f}$ 为不可压缩流的密度。
$P^{a}$ 为声压。
$P_{f}$ 为不可压缩流的压力。
$v_{f}$ 为不可压缩流的速度。

声压重构梯度的平滑算法

声压重构梯度用于计算面形心处的声压。可以平滑此重构梯度，消除局部升高值，以帮助收敛。平滑算法使用以下步骤：

计算局部平均归一化和平滑梯度 $\bar{\nabla (p_{c}^{a})}$ ：
图 18. EQUATION_DISPLAY

$\bar{\nabla (p_{c}^{a})} = \frac{\frac{\nabla (p_{c}^{a})}{{‖ \nabla (p_{c}^{a}) ‖}^{2}} + \sum_{k \in {N B}_{c}}^{} \frac{\nabla (p_{k}^{a})}{{‖ \nabla (p_{k}^{a}) ‖}^{2}}}{\frac{1}{{‖ \nabla (p_{c}^{a}) ‖}^{2}} + \sum_{k \in {N B}_{c}}^{} \frac{1}{{‖ \nabla (p_{k}^{a}) ‖}^{2}}}$

(4728)
其中， ${N B}_{c}$ 是通过面与网格单元 $c$ 直接相邻的网格单元集。
将平滑梯度 $\bar{\nabla (p_{c}^{a})}$ 复制到工作梯度 $\nabla (p_{c}^{a})$ 。
根据需要，多次重复步骤 1 和 2 以消除局部升高值。