From 99c0ae37c7dd14e4573d78a2e2824aae2cfafa01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20M=2E=20Requena=20Plens?= Date: Sun, 19 Jul 2026 16:58:55 +0200 Subject: [PATCH 1/5] feat(electroacoustics): radiating circular piston in an infinite baffle Add the canonical rigid circular piston radiator (Beranek & Mellow 2e, sections 4.19 and 13.7): the radiation impedance from the piston resistance R1 = 1 - 2 J1(2ka)/2ka and reactance X1 = 2 H1(2ka)/2ka, the low-frequency radiation mass 8 rho a^3 / 3, the far-field directivity 2 J1(ka sin t)/(ka sin t) and the directivity index. radiating_piston returns a frozen RadiatingPistonResult with .plot(); piston_resistance, piston_reactance and piston_directivity are exposed as building blocks. Documented in the electroacoustics guide. --- docs/electroacoustics.md | 38 +++ .../docs/es/guides/electroacoustics.mdx | 45 ++- .../content/docs/guides/electroacoustics.mdx | 44 ++- .../reference/api/electroacoustics/piston.md | 209 +++++++++++++ src/phonometry/_plot/electroacoustics.py | 32 ++ src/phonometry/electroacoustics/__init__.py | 12 + src/phonometry/electroacoustics/piston.py | 289 ++++++++++++++++++ tests/electroacoustics/test_piston.py | 127 ++++++++ 8 files changed, 794 insertions(+), 2 deletions(-) create mode 100644 site/src/content/docs/reference/api/electroacoustics/piston.md create mode 100644 src/phonometry/electroacoustics/piston.py create mode 100644 tests/electroacoustics/test_piston.py diff --git a/docs/electroacoustics.md b/docs/electroacoustics.md index a0466b48..0b55b3b5 100644 --- a/docs/electroacoustics.md +++ b/docs/electroacoustics.md @@ -289,6 +289,44 @@ hide in that sentence: "1 m" figure is then a referred quantity, not what a microphone placed at 1 m would read. +## 6. Radiating piston: radiation impedance and directivity + +The rigid circular **piston in an infinite baffle** is the canonical radiator +behind a loudspeaker cone, the open end of a duct and the radiation efficiency +of any finite vibrating surface (Beranek & Mellow §4.19, §13.7). Its mechanical +radiation impedance is `Z_r = rho c S (R1 + j X1)` with `S = pi a^2` and the +dimensionless resistance and reactance functions (Eqs. (13.117), (13.118)) + +``` +R1(x) = 1 - 2 J1(x) / x , X1(x) = 2 H1(x) / x , x = 2ka, +``` + +with `J1` the Bessel function and `H1` the Struve function of order one. At low +frequency `R1 -> (ka)^2 / 2` and the reactance is mass-like with the radiation +mass `M_r = 8 rho a^3 / 3` (Eq. (4.151)); at high frequency `R1 -> 1`, `X1 -> 0` +and the piston radiates as into an infinite tube. The far field follows the +directivity `D(theta) = 2 J1(ka sin theta) / (ka sin theta)`, whose first null +is at `ka sin theta = 3.8317`. + +```python +import numpy as np +from phonometry import radiating_piston + +res = radiating_piston(radius=0.1, frequencies=np.geomspace(20, 20000, 200), + angles=np.linspace(0.0, np.pi / 2, 91)) +print(round(res.radiation_mass, 4)) # 8 rho a^3 / 3, kg +print(round(float(res.directivity_index[0]), 2)) # 3.01 dB half-space limit +res.plot() # R1 and X1 vs ka +``` + +`radiating_piston` returns a `RadiatingPistonResult` with the normalized +`resistance`/`reactance`, the mechanical `radiation_resistance`/`radiation_reactance`, +the `radiation_mass`, the `directivity_index`, the far-field `directivity` +pattern (when `angles` are given) and `.plot()`. The building blocks +`piston_resistance`, `piston_reactance` and `piston_directivity` are also +callable directly. The piston is the companion radiator of the +[industrial noise-control silencers](noise-control.md). + ## References - Beranek, L. L., & Mellow, T. J. (2012). *Acoustics: Sound fields and diff --git a/site/src/content/docs/es/guides/electroacoustics.mdx b/site/src/content/docs/es/guides/electroacoustics.mdx index fd407f46..0adae0fa 100644 --- a/site/src/content/docs/es/guides/electroacoustics.mdx +++ b/site/src/content/docs/es/guides/electroacoustics.mdx @@ -355,6 +355,49 @@ esconden dos normalizaciones: campo cercano aún no haya terminado, y la cifra «a 1 m» es entonces una magnitud referida, no lo que leería un micrófono colocado a 1 m. +## 6. Pistón radiante: impedancia de radiación y directividad + +El **pistón circular rígido en pantalla infinita** es el radiador canónico tras +un cono de altavoz, el extremo abierto de un conducto y la eficiencia de +radiación de cualquier superficie vibrante finita (Beranek y Mellow §4.19, +§13.7). Su impedancia mecánica de radiación es $Z_r = \rho c S\,(R_1 + jX_1)$ con +$S = \pi a^2$ y las funciones adimensionales de resistencia y reactancia +(Ecs. (13.117), (13.118)) + +$$ +R_1(x) = 1 - \frac{2 J_1(x)}{x}, \qquad X_1(x) = \frac{2 H_1(x)}{x}, +\qquad x = 2ka, +$$ + +donde $J_1$ es la función de Bessel y $H_1$ la de Struve de orden uno. A baja +frecuencia $R_1 \to (ka)^2/2$ y la reactancia es de tipo másico con la **masa de +radiación** $M_r = 8\rho a^3/3$ (Ec. (4.151)); a alta frecuencia $R_1 \to 1$, +$X_1 \to 0$ y el pistón radia como en un tubo infinito. El campo lejano sigue la +directividad $D(\theta) = 2 J_1(ka\sin\theta)/(ka\sin\theta)$, cuyo primer nulo +está en $ka\sin\theta = 3.8317$. + +```python +import numpy as np +from phonometry import radiating_piston + +res = radiating_piston(radius=0.1, frequencies=np.geomspace(20, 20000, 200), + angles=np.linspace(0.0, np.pi / 2, 91)) +print(round(res.radiation_mass, 4)) # 8 rho a^3 / 3, kg +print(round(float(res.directivity_index[0]), 2)) # 3.01 dB, límite de semiespacio +res.plot() # R1 y X1 vs ka +``` + +`radiating_piston` devuelve un `RadiatingPistonResult` con la `resistance` y +`reactance` normalizadas, la `radiation_resistance`/`radiation_reactance` +mecánicas, la `radiation_mass`, el `directivity_index`, el patrón de campo lejano +`directivity` (cuando se dan `angles`) y `.plot()`. Los bloques +`piston_resistance`, `piston_reactance` y `piston_directivity` también son +invocables directamente. El pistón es el modelo de radiador complementario de los +[silenciadores de control de ruido industrial](/phonometry/es/guides/noise-control/). + ## Véase también -- Referencia de la API: [`electroacoustics.distortion`](/phonometry/es/reference/api/electroacoustics/distortion/) y [`electroacoustics.frequency_response`](/phonometry/es/reference/api/electroacoustics/frequency-response/). +- [Control de ruido industrial](/phonometry/es/guides/noise-control/): + silenciadores reactivos, métodos HVAC de conducto y cerramientos de máquina, + que comparten el pistón radiante como radiador canónico. +- Referencia de la API: [`electroacoustics.distortion`](/phonometry/es/reference/api/electroacoustics/distortion/), [`electroacoustics.frequency_response`](/phonometry/es/reference/api/electroacoustics/frequency-response/) y [`electroacoustics.piston`](/phonometry/es/reference/api/electroacoustics/piston/). diff --git a/site/src/content/docs/guides/electroacoustics.mdx b/site/src/content/docs/guides/electroacoustics.mdx index 8424023e..2f808ea4 100644 --- a/site/src/content/docs/guides/electroacoustics.mdx +++ b/site/src/content/docs/guides/electroacoustics.mdx @@ -346,6 +346,48 @@ hide in that sentence: "1 m" figure is then a referred quantity, not what a microphone placed at 1 m would read. +## 6. Radiating piston: radiation impedance and directivity + +The rigid circular **piston in an infinite baffle** is the canonical radiator +behind a loudspeaker cone, the open end of a duct and the radiation efficiency +of any finite vibrating surface (Beranek & Mellow §4.19, §13.7). Its mechanical +radiation impedance is $Z_r = \rho c S\,(R_1 + jX_1)$ with $S = \pi a^2$ and the +dimensionless resistance and reactance functions (Eqs. (13.117), (13.118)) + +$$ +R_1(x) = 1 - \frac{2 J_1(x)}{x}, \qquad X_1(x) = \frac{2 H_1(x)}{x}, +\qquad x = 2ka, +$$ + +where $J_1$ is the Bessel function and $H_1$ the Struve function of order one. +At low frequency $R_1 \to (ka)^2/2$ and the reactance is mass-like with the +**radiation mass** $M_r = 8\rho a^3/3$ (Eq. (4.151)); at high frequency +$R_1 \to 1$, $X_1 \to 0$ and the piston radiates as if into an infinite tube. +The far field follows the directivity $D(\theta) = 2 J_1(ka\sin\theta)/(ka\sin\theta)$, +whose first null is at $ka\sin\theta = 3.8317$. + +```python +import numpy as np +from phonometry import radiating_piston + +res = radiating_piston(radius=0.1, frequencies=np.geomspace(20, 20000, 200), + angles=np.linspace(0.0, np.pi / 2, 91)) +print(round(res.radiation_mass, 4)) # 8 rho a^3 / 3, kg +print(round(float(res.directivity_index[0]), 2)) # 3.01 dB half-space limit +res.plot() # R1 and X1 vs ka +``` + +`radiating_piston` returns a `RadiatingPistonResult` with the normalized +`resistance`/`reactance`, the mechanical `radiation_resistance`/`radiation_reactance`, +the `radiation_mass`, the `directivity_index`, the far-field `directivity` +pattern (when `angles` are given) and `.plot()`. The building blocks +`piston_resistance`, `piston_reactance` and `piston_directivity` are also +callable directly. The piston is the companion radiator model of the +[industrial noise-control silencers](/phonometry/guides/noise-control/). + ## See also -- API reference: [`electroacoustics.distortion`](/phonometry/reference/api/electroacoustics/distortion/) and [`electroacoustics.frequency_response`](/phonometry/reference/api/electroacoustics/frequency-response/). +- [Industrial noise control](/phonometry/guides/noise-control/): reactive + silencers, HVAC duct methods and machine enclosures, sharing the radiating + piston as the canonical radiator. +- API reference: [`electroacoustics.distortion`](/phonometry/reference/api/electroacoustics/distortion/), [`electroacoustics.frequency_response`](/phonometry/reference/api/electroacoustics/frequency-response/) and [`electroacoustics.piston`](/phonometry/reference/api/electroacoustics/piston/). diff --git a/site/src/content/docs/reference/api/electroacoustics/piston.md b/site/src/content/docs/reference/api/electroacoustics/piston.md new file mode 100644 index 00000000..a8a01296 --- /dev/null +++ b/site/src/content/docs/reference/api/electroacoustics/piston.md @@ -0,0 +1,209 @@ +--- +title: "electroacoustics.piston" +description: "Public API of phonometry.electroacoustics.piston (auto-generated)." +sidebar: + label: "piston" +--- + +> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand. + +Radiation of a rigid circular piston set in an infinite baffle. + +The baffled circular piston is the canonical acoustic radiator: a flat rigid +disc of radius `a` vibrating with a uniform normal velocity in an otherwise +rigid infinite plane. It is the model behind a loudspeaker cone in a large +cabinet, the open end of a duct, and the reference source for the radiation +efficiency of any finite vibrating surface, so its two results -- the +**radiation impedance** the air presents to the piston and the **directivity** +of the far field -- are the base of the electroacoustics domain (Beranek & +Mellow, *Acoustics: Sound Fields, Transducers and Vibration* 2nd ed., §4.4; +Bies, Hansen & Howard, *Engineering Noise Control* 5th ed.). + +**Radiation impedance.** The reaction force of the air on the piston is +`F = Z_r u` with the mechanical radiation impedance + + Z_r = rho c S ( R1(2ka) + j X1(2ka) ), S = pi a^2, + +where `k = omega / c` is the wavenumber, `rho c` the characteristic +impedance of air and `S` the piston area. The dimensionless **piston +resistance** and **reactance** functions are (Beranek & Mellow Eq. (4.30)) + + R1(x) = 1 - 2 J1(x) / x, X1(x) = 2 H1(x) / x, + +with `J1` the Bessel function of the first kind and `H1` the Struve +function, both of order one, evaluated at `x = 2ka`. + +* **Low frequency** (`ka << 1`): `R1 -> (ka)^2 / 2` so the radiated power + rises as `f^2`, and `X1 -> (8 / 3 pi) ka`. The reactance is mass-like, + `X_r = rho c S X1 = omega M_r` with the **radiation (accreted) mass** + + M_r = 8 rho a^3 / 3 + + (Beranek & Mellow Eq. (4.32)): the piston drags an extra `8 rho a^3 / 3` of + air, equivalent to a layer `8a / 3 pi` thick over its face. +* **High frequency** (`ka >> 1`): `R1 -> 1` and `X1 -> 0`, so + `Z_r -> rho c S` -- the piston radiates as if into an infinite tube and the + air loads it purely resistively. + +**Directivity.** The far-field pressure of the baffled piston varies with the +polar angle `theta` from the axis as (Beranek & Mellow Eq. (4.42)) + + D(theta) = 2 J1(ka sin theta) / (ka sin theta), D(0) = 1. + +The main lobe narrows as `ka` grows; its first null is at +`ka sin theta = 3.8317` (the first zero of `J1`), which exists only once +`ka > 3.8317`. The **directivity factor** `Q` (on-axis intensity over the +intensity of a point source of equal power radiating into the full sphere) and +the **directivity index** `DI = 10 log10 Q` follow from integrating +`|D|^2` over the radiating hemisphere, + + Q = 2 / integral_0^(pi/2) |D(theta)|^2 sin theta d theta, + +which tends to `Q = 2` (`DI = 3.01 dB`, the half-space baffle gain) at low +`ka` and to `Q ~ (ka)^2` (`DI ~ 20 log10 ka`) at high `ka`. + +## piston_directivity + +```python +piston_directivity(ka: ArrayLike, theta: ArrayLike) -> np.ndarray | float +``` + +Far-field directivity `D = 2 J1(ka sin theta) / (ka sin theta)`. + +The pressure amplitude of a baffled circular piston relative to its on-axis +value (Beranek & Mellow Eq. (4.42)), normalized so `D(0) = 1`. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `ka` | Wavenumber-radius product `ka` (scalar or array). | +| `theta` | Polar angle from the axis, rad (scalar or array). Broadcast against `ka`. | + +**Returns:** `D` (float for scalar inputs, else an array). + +## piston_reactance + +```python +piston_reactance(x: ArrayLike) -> np.ndarray | float +``` + +Piston reactance function `X1(x) = 2 H1(x) / x` (`H1` Struve order 1). + +The imaginary part of the normalized radiation impedance of a baffled +circular piston (Beranek & Mellow Eq. (4.30)). It rises as +`(8 / 3 pi) ka` (mass-like) at low `x = 2ka` and decays to 0 at high +`x`. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `x` | Argument `x = 2ka` (scalar or array), dimensionless. | + +**Returns:** `X1(x)` (float for scalar input, else an array). + +## piston_resistance + +```python +piston_resistance(x: ArrayLike) -> np.ndarray | float +``` + +Piston resistance function `R1(x) = 1 - 2 J1(x) / x`. + +The real part of the normalized radiation impedance of a baffled circular +piston, as a function of `x = 2ka` (Beranek & Mellow Eq. (4.30)). It +rises as `x^2 / 8 = (ka)^2 / 2` at low `x` and tends to 1 at high `x`. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `x` | Argument `x = 2ka` (scalar or array), dimensionless. | + +**Returns:** `R1(x)` (float for scalar input, else an array). + +## radiating_piston + +```python +radiating_piston( + radius: float, + frequencies: ArrayLike, + *, + speed_of_sound: float = 343.0, + density: float = 1.206, + angles: ArrayLike | None = None, +) -> RadiatingPistonResult +``` + +Radiation impedance and directivity of a rigid baffled circular piston. + +Evaluates the piston resistance `R1(2ka)` and reactance `X1(2ka)`, the +mechanical radiation impedance `rho c S (R1 + j X1)`, the low-frequency +radiation mass `8 rho a^3 / 3` and the directivity index over the given +frequencies (Beranek & Mellow §4.4). Pass `angles` to also sample the +far-field directivity pattern `D(theta)`. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `radius` | Piston radius `a`, m. | +| `frequencies` | Frequencies `f`, Hz (scalar or 1-D array), all > 0. | +| `speed_of_sound` | Speed of sound `c`, m/s (default 343). | +| `density` | Air density `rho`, kg/m3 (default 1.206). | +| `angles` | Optional polar angles `theta` from the axis, rad, at which to sample the directivity pattern. | + +**Returns:** A [`RadiatingPistonResult`](/phonometry/reference/api/electroacoustics/piston/#radiatingpistonresult). + +## RadiatingPistonResult + +```python +RadiatingPistonResult( + frequencies: np.ndarray, + ka: np.ndarray, + resistance: np.ndarray, + reactance: np.ndarray, + radiation_resistance: np.ndarray, + radiation_reactance: np.ndarray, + radiation_mass: float, + directivity_index: np.ndarray, + angles: np.ndarray | None, + directivity: np.ndarray | None, + radius: float, + speed_of_sound: float, + density: float, +) +``` + +Radiation impedance and directivity of a baffled circular piston. + +**Attributes** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz. | +| `ka` | Wavenumber-radius product `ka` at each frequency. | +| `resistance` | Normalized piston resistance `R1(2ka)` (real part of `Z_r / (rho c S)`). | +| `reactance` | Normalized piston reactance `X1(2ka)` (imaginary part of `Z_r / (rho c S)`). | +| `radiation_resistance` | Mechanical radiation resistance `rho c S R1`, N s/m. | +| `radiation_reactance` | Mechanical radiation reactance `rho c S X1`, N s/m. | +| `radiation_mass` | Low-frequency accreted air mass `M_r = 8 rho a^3/3`, kg (a single value; the mass limit of `radiation_reactance / omega`). | +| `directivity_index` | Directivity index `DI = 10 log10 Q`, dB. | +| `angles` | Polar angles of `directivity`, rad, or `None` if not requested. | +| `directivity` | Far-field directivity `D(theta)` as a `(n_freq, n_angle)` array, or `None` if `angles` was not given. | +| `radius` | Piston radius `a`, m. | +| `speed_of_sound` | Speed of sound `c`, m/s. | +| `density` | Air density `rho`, kg/m3. | + +### RadiatingPistonResult.plot() + +```python +RadiatingPistonResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes +``` + +Plot the normalized piston resistance and reactance against `ka`. + +Reproduces the classic Beranek & Mellow figure: `R1` rising to 1 and +`X1` peaking then decaying, over the `ka` range of the result. +Requires matplotlib (`pip install phonometry[plot]`). diff --git a/src/phonometry/_plot/electroacoustics.py b/src/phonometry/_plot/electroacoustics.py index 304f85b3..c5ac2b23 100644 --- a/src/phonometry/_plot/electroacoustics.py +++ b/src/phonometry/_plot/electroacoustics.py @@ -22,6 +22,7 @@ from matplotlib.axes import Axes from ..electroacoustics.distortion import HarmonicDistortionResult from ..electroacoustics.frequency_response import FrequencyResponseResult + from ..electroacoustics.piston import RadiatingPistonResult from ..electroacoustics.swept_sine import SweptSineDistortionResult #: Shared frequency-axis label of the electroacoustics renderers. @@ -200,3 +201,34 @@ def _thd_panel(axt: Axes) -> None: axes[1].set_xlabel("Excitation frequency [Hz]") format_frequency_axis(axes[1]) return axes + + +def plot_piston_impedance( + result: "RadiatingPistonResult", ax: Axes | None = None, **kwargs: Any +) -> Axes: + """Normalized radiation resistance and reactance of a baffled piston. + + Draws the piston resistance ``R1(2ka)`` and reactance ``X1(2ka)`` against + ``ka`` on a logarithmic ``ka`` axis (the classic Beranek & Mellow figure): + ``R1`` rising to 1, ``X1`` peaking near ``ka ~ 0.5`` then decaying. + + :param result: A :class:`~phonometry.electroacoustics.piston.RadiatingPistonResult`. + :param ax: Existing axes, or ``None`` to create a figure. + :param kwargs: Forwarded to the resistance ``Axes.plot`` (its primary curve). + :return: The axes. + """ + ax = ax if ax is not None else _new_axes() + ka = np.asarray(result.ka, dtype=np.float64) + kwargs.setdefault("color", _C_PRIMARY) + kwargs.setdefault("label", r"$R_1$ (resistance)") + ax.semilogx(ka, np.asarray(result.resistance), lw=1.8, **kwargs) + ax.semilogx(ka, np.asarray(result.reactance), color=_C_SECONDARY, lw=1.8, + ls="--", label=r"$X_1$ (reactance)") + ax.axhline(1.0, color=_C_TERTIARY, ls=":", lw=1.0, + label=r"$R_1 \to 1$ ($ka \gg 1$)") + ax.set_xlabel(r"$ka$") + ax.set_ylabel(r"Normalized radiation impedance $Z_r / \rho c S$") + ax.set_title("Baffled circular piston radiation impedance") + ax.grid(True, which="both", alpha=0.3) + ax.legend(loc="best", fontsize="small") + return ax diff --git a/src/phonometry/electroacoustics/__init__.py b/src/phonometry/electroacoustics/__init__.py index 02659bd7..96093d87 100644 --- a/src/phonometry/electroacoustics/__init__.py +++ b/src/phonometry/electroacoustics/__init__.py @@ -21,6 +21,13 @@ weighted_thd, ) from .frequency_response import FrequencyResponseResult, coherence, transfer_function +from .piston import ( + RadiatingPistonResult, + piston_directivity, + piston_reactance, + piston_resistance, + radiating_piston, +) from .swept_sine import ( SweptSineDistortionResult, swept_sine_distortion, @@ -31,6 +38,7 @@ "FrequencyResponseResult", "HarmonicDistortionResult", "ModulationDistortionResult", + "RadiatingPistonResult", "SweptSineDistortionResult", "coherence", "difference_frequency_distortion", @@ -41,6 +49,10 @@ "idle_channel_noise", "itu_r_468_weighting", "modulation_distortion", + "piston_directivity", + "piston_reactance", + "piston_resistance", + "radiating_piston", "sinad", "swept_sine_distortion", "synchronized_sweep_signal", diff --git a/src/phonometry/electroacoustics/piston.py b/src/phonometry/electroacoustics/piston.py new file mode 100644 index 00000000..dd7f3212 --- /dev/null +++ b/src/phonometry/electroacoustics/piston.py @@ -0,0 +1,289 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +""" +Radiation of a rigid circular piston set in an infinite baffle. + +The baffled circular piston is the canonical acoustic radiator: a flat rigid +disc of radius ``a`` vibrating with a uniform normal velocity in an otherwise +rigid infinite plane. It is the model behind a loudspeaker cone in a large +cabinet, the open end of a duct, and the reference source for the radiation +efficiency of any finite vibrating surface, so its two results -- the +**radiation impedance** the air presents to the piston and the **directivity** +of the far field -- are the base of the electroacoustics domain (Beranek & +Mellow, *Acoustics: Sound Fields, Transducers and Vibration* 2nd ed., §4.4; +Bies, Hansen & Howard, *Engineering Noise Control* 5th ed.). + +**Radiation impedance.** The reaction force of the air on the piston is +``F = Z_r u`` with the mechanical radiation impedance + + Z_r = rho c S ( R1(2ka) + j X1(2ka) ), S = pi a^2, + +where ``k = omega / c`` is the wavenumber, ``rho c`` the characteristic +impedance of air and ``S`` the piston area. The dimensionless **piston +resistance** and **reactance** functions are (Beranek & Mellow Eq. (4.30)) + + R1(x) = 1 - 2 J1(x) / x, X1(x) = 2 H1(x) / x, + +with ``J1`` the Bessel function of the first kind and ``H1`` the Struve +function, both of order one, evaluated at ``x = 2ka``. + +* **Low frequency** (``ka << 1``): ``R1 -> (ka)^2 / 2`` so the radiated power + rises as ``f^2``, and ``X1 -> (8 / 3 pi) ka``. The reactance is mass-like, + ``X_r = rho c S X1 = omega M_r`` with the **radiation (accreted) mass** + + M_r = 8 rho a^3 / 3 + + (Beranek & Mellow Eq. (4.32)): the piston drags an extra ``8 rho a^3 / 3`` of + air, equivalent to a layer ``8a / 3 pi`` thick over its face. +* **High frequency** (``ka >> 1``): ``R1 -> 1`` and ``X1 -> 0``, so + ``Z_r -> rho c S`` -- the piston radiates as if into an infinite tube and the + air loads it purely resistively. + +**Directivity.** The far-field pressure of the baffled piston varies with the +polar angle ``theta`` from the axis as (Beranek & Mellow Eq. (4.42)) + + D(theta) = 2 J1(ka sin theta) / (ka sin theta), D(0) = 1. + +The main lobe narrows as ``ka`` grows; its first null is at +``ka sin theta = 3.8317`` (the first zero of ``J1``), which exists only once +``ka > 3.8317``. The **directivity factor** ``Q`` (on-axis intensity over the +intensity of a point source of equal power radiating into the full sphere) and +the **directivity index** ``DI = 10 log10 Q`` follow from integrating +``|D|^2`` over the radiating hemisphere, + + Q = 2 / integral_0^(pi/2) |D(theta)|^2 sin theta d theta, + +which tends to ``Q = 2`` (``DI = 3.01 dB``, the half-space baffle gain) at low +``ka`` and to ``Q ~ (ka)^2`` (``DI ~ 20 log10 ka``) at high ``ka``. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import numpy as np +from numpy.typing import ArrayLike, NDArray +from scipy import special + +from .._internal.validation import require_positive + +if TYPE_CHECKING: + from matplotlib.axes import Axes + +#: First zero of the Bessel function ``J1``; the first directivity null of the +#: baffled piston sits at ``ka sin theta = J1_FIRST_ZERO`` (Beranek & Mellow). +J1_FIRST_ZERO = 3.8317059702075125 + +#: Reference characteristic impedance factors of air at 20 degC, 101.325 kPa. +_RHO_AIR = 1.206 +_C_AIR = 343.0 + + +def piston_resistance(x: ArrayLike) -> np.ndarray | float: + """Piston resistance function ``R1(x) = 1 - 2 J1(x) / x``. + + The real part of the normalized radiation impedance of a baffled circular + piston, as a function of ``x = 2ka`` (Beranek & Mellow Eq. (4.30)). It + rises as ``x^2 / 8 = (ka)^2 / 2`` at low ``x`` and tends to 1 at high ``x``. + + :param x: Argument ``x = 2ka`` (scalar or array), dimensionless. + :return: ``R1(x)`` (float for scalar input, else an array). + """ + arr = np.asarray(x, dtype=np.float64) + if np.any(arr < 0.0) or not np.all(np.isfinite(arr)): + raise ValueError("'x' must be non-negative and finite.") + with np.errstate(invalid="ignore", divide="ignore"): + out = 1.0 - 2.0 * special.j1(arr) / arr + # J1(x)/x -> 1/2 as x -> 0, so R1 -> 0; fill the removable singularity. + out = np.where(arr == 0.0, 0.0, out) + return out[()] if out.ndim == 0 else out + + +def piston_reactance(x: ArrayLike) -> np.ndarray | float: + """Piston reactance function ``X1(x) = 2 H1(x) / x`` (``H1`` Struve order 1). + + The imaginary part of the normalized radiation impedance of a baffled + circular piston (Beranek & Mellow Eq. (4.30)). It rises as + ``(8 / 3 pi) ka`` (mass-like) at low ``x = 2ka`` and decays to 0 at high + ``x``. + + :param x: Argument ``x = 2ka`` (scalar or array), dimensionless. + :return: ``X1(x)`` (float for scalar input, else an array). + """ + arr = np.asarray(x, dtype=np.float64) + if np.any(arr < 0.0) or not np.all(np.isfinite(arr)): + raise ValueError("'x' must be non-negative and finite.") + with np.errstate(invalid="ignore", divide="ignore"): + out = 2.0 * special.struve(1, arr) / arr + # H1(x)/x -> 0 as x -> 0 (H1 ~ 2 x^2 / 3 pi), so X1 -> 0. + out = np.where(arr == 0.0, 0.0, out) + return out[()] if out.ndim == 0 else out + + +def piston_directivity(ka: ArrayLike, theta: ArrayLike) -> np.ndarray | float: + """Far-field directivity ``D = 2 J1(ka sin theta) / (ka sin theta)``. + + The pressure amplitude of a baffled circular piston relative to its on-axis + value (Beranek & Mellow Eq. (4.42)), normalized so ``D(0) = 1``. + + :param ka: Wavenumber-radius product ``ka`` (scalar or array). + :param theta: Polar angle from the axis, rad (scalar or array). Broadcast + against ``ka``. + :return: ``D`` (float for scalar inputs, else an array). + """ + ka_arr = np.asarray(ka, dtype=np.float64) + theta_arr = np.asarray(theta, dtype=np.float64) + if np.any(ka_arr < 0.0) or not np.all(np.isfinite(ka_arr)): + raise ValueError("'ka' must be non-negative and finite.") + if not np.all(np.isfinite(theta_arr)): + raise ValueError("'theta' must be finite.") + u = ka_arr * np.sin(theta_arr) + with np.errstate(invalid="ignore", divide="ignore"): + out = 2.0 * special.j1(u) / u + # 2 J1(u)/u -> 1 as u -> 0 (on-axis, or ka = 0). + out = np.where(u == 0.0, 1.0, out) + return out[()] if out.ndim == 0 else out + + +def _directivity_index(ka: NDArray[np.float64]) -> NDArray[np.float64]: + """Directivity index ``DI = 10 log10 Q`` per ``ka`` by hemisphere quadrature. + + ``Q = 2 / integral_0^(pi/2) |D|^2 sin theta d theta`` (radiation into the + half-space in front of the baffle). Uses a fine trapezoid grid; ``ka -> 0`` + gives the half-space baffle gain ``DI = 10 log10 2 = 3.01 dB`` and high + ``ka`` tends to ``Q = (ka)^2`` (``DI = 20 log10 ka``). + """ + theta = np.linspace(0.0, 0.5 * np.pi, 2001) + sin_t = np.sin(theta) + di = np.empty_like(ka) + for i, kai in enumerate(ka): + d = np.asarray(piston_directivity(float(kai), theta), dtype=np.float64) + integrand = d**2 * sin_t + integral = float(np.trapezoid(integrand, theta)) + di[i] = 10.0 * np.log10(2.0 / integral) + return di + + +@dataclass(frozen=True) +class RadiatingPistonResult: + """Radiation impedance and directivity of a baffled circular piston. + + :ivar frequencies: Frequencies ``f``, Hz. + :ivar ka: Wavenumber-radius product ``ka`` at each frequency. + :ivar resistance: Normalized piston resistance ``R1(2ka)`` (real part of + ``Z_r / (rho c S)``). + :ivar reactance: Normalized piston reactance ``X1(2ka)`` (imaginary part of + ``Z_r / (rho c S)``). + :ivar radiation_resistance: Mechanical radiation resistance + ``rho c S R1``, N s/m. + :ivar radiation_reactance: Mechanical radiation reactance + ``rho c S X1``, N s/m. + :ivar radiation_mass: Low-frequency accreted air mass ``M_r = 8 rho a^3/3``, + kg (a single value; the mass limit of ``radiation_reactance / omega``). + :ivar directivity_index: Directivity index ``DI = 10 log10 Q``, dB. + :ivar angles: Polar angles of ``directivity``, rad, or ``None`` if not + requested. + :ivar directivity: Far-field directivity ``D(theta)`` as a + ``(n_freq, n_angle)`` array, or ``None`` if ``angles`` was not given. + :ivar radius: Piston radius ``a``, m. + :ivar speed_of_sound: Speed of sound ``c``, m/s. + :ivar density: Air density ``rho``, kg/m3. + """ + + frequencies: np.ndarray + ka: np.ndarray + resistance: np.ndarray + reactance: np.ndarray + radiation_resistance: np.ndarray + radiation_reactance: np.ndarray + radiation_mass: float + directivity_index: np.ndarray + angles: np.ndarray | None + directivity: np.ndarray | None + radius: float + speed_of_sound: float + density: float + + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": + """Plot the normalized piston resistance and reactance against ``ka``. + + Reproduces the classic Beranek & Mellow figure: ``R1`` rising to 1 and + ``X1`` peaking then decaying, over the ``ka`` range of the result. + Requires matplotlib (``pip install phonometry[plot]``). + """ + from .._plot.electroacoustics import plot_piston_impedance + + return plot_piston_impedance(self, ax=ax, **kwargs) + + +def radiating_piston( + radius: float, + frequencies: ArrayLike, + *, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, + angles: ArrayLike | None = None, +) -> RadiatingPistonResult: + """Radiation impedance and directivity of a rigid baffled circular piston. + + Evaluates the piston resistance ``R1(2ka)`` and reactance ``X1(2ka)``, the + mechanical radiation impedance ``rho c S (R1 + j X1)``, the low-frequency + radiation mass ``8 rho a^3 / 3`` and the directivity index over the given + frequencies (Beranek & Mellow §4.4). Pass ``angles`` to also sample the + far-field directivity pattern ``D(theta)``. + + :param radius: Piston radius ``a``, m. + :param frequencies: Frequencies ``f``, Hz (scalar or 1-D array), all > 0. + :param speed_of_sound: Speed of sound ``c``, m/s (default 343). + :param density: Air density ``rho``, kg/m3 (default 1.206). + :param angles: Optional polar angles ``theta`` from the axis, rad, at which + to sample the directivity pattern. + :return: A :class:`RadiatingPistonResult`. + """ + a = require_positive(radius, "radius") + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + f = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) + if f.ndim != 1 or f.size == 0: + raise ValueError("'frequencies' must be a non-empty 1-D array.") + if np.any(f <= 0.0) or not np.all(np.isfinite(f)): + raise ValueError("'frequencies' must be positive and finite.") + + omega = 2.0 * np.pi * f + k = omega / c + ka = k * a + area = np.pi * a**2 + r1 = np.asarray(piston_resistance(2.0 * ka), dtype=np.float64) + x1 = np.asarray(piston_reactance(2.0 * ka), dtype=np.float64) + rho_c_s = rho * c * area + radiation_mass = 8.0 * rho * a**3 / 3.0 + di = _directivity_index(ka) + + directivity: np.ndarray | None = None + angle_arr: np.ndarray | None = None + if angles is not None: + angle_arr = np.atleast_1d(np.asarray(angles, dtype=np.float64)) + if angle_arr.ndim != 1 or angle_arr.size == 0: + raise ValueError("'angles' must be a non-empty 1-D array.") + if not np.all(np.isfinite(angle_arr)): + raise ValueError("'angles' must be finite.") + directivity = np.asarray( + piston_directivity(ka[:, None], angle_arr[None, :]), + dtype=np.float64, + ) + + return RadiatingPistonResult( + frequencies=f, + ka=ka, + resistance=r1, + reactance=x1, + radiation_resistance=rho_c_s * r1, + radiation_reactance=rho_c_s * x1, + radiation_mass=radiation_mass, + directivity_index=di, + angles=angle_arr, + directivity=directivity, + radius=a, + speed_of_sound=c, + density=rho, + ) diff --git a/tests/electroacoustics/test_piston.py b/tests/electroacoustics/test_piston.py new file mode 100644 index 00000000..e2bcef6f --- /dev/null +++ b/tests/electroacoustics/test_piston.py @@ -0,0 +1,127 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +"""Tests for the baffled circular piston (Beranek & Mellow 2e, §4.19 / §13.7). + +Oracles are the exact closed forms (``R1 = 1 - 2 J1(x)/x``, ``X1 = 2 H1(x)/x`` +with ``x = 2ka``), their low- and high-frequency limits (Eqs. (13.117), +(13.118), radiation mass Eq. (4.151)), the first directivity null at the first +zero of ``J1`` (3.8317), and the half-space baffle directivity index of 3.01 dB +at ``ka -> 0``. Point values were cross-checked against scipy in the source +extraction. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest +from scipy import special + +from phonometry import ( + RadiatingPistonResult, + piston_directivity, + piston_reactance, + piston_resistance, + radiating_piston, +) +from phonometry.electroacoustics.piston import J1_FIRST_ZERO + + +def test_resistance_reactance_closed_form() -> None: + x = np.array([0.1, 0.5, 1.0, 2.0, 3.0, 5.0, 8.0, 10.0]) + assert np.allclose(piston_resistance(x), 1.0 - 2.0 * special.j1(x) / x) + assert np.allclose(piston_reactance(x), 2.0 * special.struve(1, x) / x) + + +def test_digitized_reference_points() -> None: + # Points verified from Eqs. (13.117)/(13.118) in the source extraction. + ref = { + 0.1: (0.001249, 0.042413), + 1.0: (0.119899, 0.396915), + 3.0: (0.773961, 0.680073), + 5.0: (1.131032, 0.323125), + 10.0: (0.991305, 0.178366), + } + for x, (r1, x1) in ref.items(): + assert piston_resistance(x) == pytest.approx(r1, abs=1e-5) + assert piston_reactance(x) == pytest.approx(x1, abs=1e-5) + + +def test_low_frequency_limits() -> None: + # R1 -> (ka)^2/2 = x^2/8 and X1 -> (8/3pi) ka = (4/3pi) x as x -> 0. + x = 1e-3 + assert piston_resistance(x) == pytest.approx(x**2 / 8.0, rel=1e-4) + assert piston_reactance(x) == pytest.approx(4.0 / (3.0 * math.pi) * x, rel=1e-4) + + +def test_high_frequency_limit_resistive() -> None: + # R1 -> 1 and X1 -> 0 at large x: Z_r -> rho c S (purely resistive). + assert piston_resistance(80.0) == pytest.approx(1.0, abs=0.02) + assert abs(piston_reactance(80.0)) < 0.05 + + +def test_zero_argument_regular() -> None: + assert piston_resistance(0.0) == 0.0 + assert piston_reactance(0.0) == 0.0 + assert not np.any(np.isnan(piston_resistance(np.array([0.0, 1.0])))) + + +def test_directivity_onaxis_and_first_null() -> None: + assert piston_directivity(5.0, 0.0) == pytest.approx(1.0) + ka = 6.0 + theta_null = math.asin(J1_FIRST_ZERO / ka) + assert abs(piston_directivity(ka, theta_null)) < 1e-9 + # No null exists when ka < first zero of J1. + theta = np.linspace(0.0, math.pi / 2, 200) + assert np.min(np.abs(piston_directivity(3.0, theta))) > 0.1 + + +def test_radiation_mass_coefficient() -> None: + a, rho = 0.12, 1.206 + res = radiating_piston(a, [100.0], density=rho) + assert res.radiation_mass == pytest.approx(8.0 * rho * a**3 / 3.0) + + +def test_radiation_reactance_tends_to_mass_times_omega() -> None: + # At low ka the mechanical reactance rho c S X1 = omega M_r. + a = 0.05 + res = radiating_piston(a, [20.0], density=1.206, speed_of_sound=343.0) + omega = 2.0 * math.pi * 20.0 + assert res.radiation_reactance[0] == pytest.approx( + omega * res.radiation_mass, rel=2e-3 + ) + + +def test_directivity_index_half_space_limit() -> None: + # DI -> 10 log10(2) = 3.01 dB as ka -> 0 (radiation into the half-space). + res = radiating_piston(0.01, [1.0]) + assert res.directivity_index[0] == pytest.approx(3.0103, abs=1e-3) + + +def test_directivity_index_high_ka() -> None: + # DI -> 20 log10(ka) at high ka (Q = (ka)^2). + a, f, c = 1.0, 5000.0, 343.0 + res = radiating_piston(a, [f], speed_of_sound=c) + ka = 2.0 * math.pi * f * a / c + assert res.directivity_index[0] == pytest.approx(20.0 * math.log10(ka), abs=0.3) + + +def test_result_shapes_and_plot() -> None: + res = radiating_piston(0.1, [100.0, 1000.0], angles=[0.0, 0.3, 0.6]) + assert isinstance(res, RadiatingPistonResult) + assert res.directivity is not None and res.directivity.shape == (2, 3) + assert res.ka.shape == (2,) + import matplotlib + + matplotlib.use("Agg") + ax = res.plot() + assert ax.get_ylabel() + + +def test_validation() -> None: + with pytest.raises(ValueError): + radiating_piston(-1.0, [100.0]) + with pytest.raises(ValueError): + radiating_piston(0.1, [0.0]) + with pytest.raises(ValueError): + piston_resistance(-1.0) From fc82d4be4af8641d0cbf063529df4d7b44960088 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20M=2E=20Requena=20Plens?= Date: Sun, 19 Jul 2026 16:59:09 +0200 Subject: [PATCH 2/5] feat(noise_control): silencers, HVAC duct methods and machine enclosures New phonometry.noise_control domain (Bies, Hansen & Howard 5e; Munjal): - Reactive silencers by the four-pole transmission-matrix method: expansion chamber (matches the closed form 10 lg[1 + (1/4)(m - 1/m)^2 sin^2 kL] exactly), Helmholtz and quarter-wave side-branch resonators, and extended-tube chambers, with transmission and insertion loss for configurable source/radiation impedances. The four-pole primitives (duct_matrix, shunt_matrix, cascade, transmission_loss, insertion_loss, helmholtz_impedance, quarter_wave_impedance) are exposed for composition. - HVAC duct methods: end reflection (ASHRAE Table 8.14), bend insertion loss (Table 8.11), plenum transmission loss (Wells' method) and flow-generated noise of straight ducts and mitred bends (VDI 2081). - Machine-enclosure insertion loss IL = R - C, where the panel R is supplied by the caller (array or callable) and the interior correction reuses the room constant of phonometry.room. Validated against the closed forms, digitised table nodes and an independent 2D FDTD cross-check of a single expansion chamber. Documented with a guide in both languages and a concept figure of the expansion-chamber TL. --- .github/images/silencer_expansion_chamber.svg | 1278 +++++++++++++++++ .../silencer_expansion_chamber_dark.svg | 1278 +++++++++++++++++ .../images/silencer_expansion_chamber_es.svg | 1278 +++++++++++++++++ .../silencer_expansion_chamber_es_dark.svg | 1278 +++++++++++++++++ docs/noise-control.md | 173 +++ scripts/generate_graphs.py | 34 + site/astro.config.mjs | 1 + .../content/docs/es/guides/noise-control.mdx | 276 ++++ .../src/content/docs/guides/noise-control.mdx | 272 ++++ .../reference/api/noise_control/enclosures.md | 106 ++ .../docs/reference/api/noise_control/hvac.md | 217 +++ .../reference/api/noise_control/silencers.md | 242 ++++ src/phonometry/_plot/noise_control.py | 126 ++ src/phonometry/noise_control/__init__.py | 51 + src/phonometry/noise_control/enclosures.py | 155 ++ src/phonometry/noise_control/hvac.py | 343 +++++ src/phonometry/noise_control/silencers.py | 564 ++++++++ tests/noise_control/test_enclosures.py | 94 ++ tests/noise_control/test_fdtd_crosscheck.py | 93 ++ tests/noise_control/test_hvac.py | 126 ++ tests/noise_control/test_silencers.py | 148 ++ tests/test_package_architecture.py | 3 + 22 files changed, 8136 insertions(+) create mode 100644 .github/images/silencer_expansion_chamber.svg create mode 100644 .github/images/silencer_expansion_chamber_dark.svg create mode 100644 .github/images/silencer_expansion_chamber_es.svg create mode 100644 .github/images/silencer_expansion_chamber_es_dark.svg create mode 100644 docs/noise-control.md create mode 100644 site/src/content/docs/es/guides/noise-control.mdx create mode 100644 site/src/content/docs/guides/noise-control.mdx create mode 100644 site/src/content/docs/reference/api/noise_control/enclosures.md create mode 100644 site/src/content/docs/reference/api/noise_control/hvac.md create mode 100644 site/src/content/docs/reference/api/noise_control/silencers.md create mode 100644 src/phonometry/_plot/noise_control.py create mode 100644 src/phonometry/noise_control/__init__.py create mode 100644 src/phonometry/noise_control/enclosures.py create mode 100644 src/phonometry/noise_control/hvac.py create mode 100644 src/phonometry/noise_control/silencers.py create mode 100644 tests/noise_control/test_enclosures.py create mode 100644 tests/noise_control/test_fdtd_crosscheck.py create mode 100644 tests/noise_control/test_hvac.py create mode 100644 tests/noise_control/test_silencers.py diff --git a/.github/images/silencer_expansion_chamber.svg b/.github/images/silencer_expansion_chamber.svg new file mode 100644 index 00000000..d922ee6c --- /dev/null +++ b/.github/images/silencer_expansion_chamber.svg @@ -0,0 +1,1278 @@ + + + + + + + + image/svg+xml + + + Matplotlib v3.11.0, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 31.5 + + + + + + + + + + + + + 63 + + + + + + + + + + + + + 125 + + + + + + + + + + + + + 250 + + + + + + + + + + + + + 500 + + + + + + + + + + + + + 1k + + + + + + + + + + + + + 2k + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Frequency [Hz] + + + + + + + + + + + + + + + + + 0.0 + + + + + + + + + + + + + 2.5 + + + + + + + + + + + + + 5.0 + + + + + + + + + + + + + 7.5 + + + + + + + + + + + + + 10.0 + + + + + + + + + + + + + 12.5 + + + + + + + + + + + + + 15.0 + + + + + + + + + + + + + 17.5 + + + + + + + + + + + + + 20.0 + + + + Transmission loss [dB] + + + + + + + + + + + + + + + + + + + + + + + + + + + + Expansion-chamber transmission loss (Bies Eq. 8.111) + + + + + + + Area ratio m = Sexp/Sduct + + + + + + m = 2 → 1.9 dB + + + + + + m = 4 → 6.5 dB + + + + + + m = 8 → 12.2 dB + + + + + + m = 16 → 18.1 dB + + + + + + + + + + diff --git a/.github/images/silencer_expansion_chamber_dark.svg b/.github/images/silencer_expansion_chamber_dark.svg new file mode 100644 index 00000000..0bb18367 --- /dev/null +++ b/.github/images/silencer_expansion_chamber_dark.svg @@ -0,0 +1,1278 @@ + + + + + + + + image/svg+xml + + + Matplotlib v3.11.0, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 31.5 + + + + + + + + + + + + + 63 + + + + + + + + + + + + + 125 + + + + + + + + + + + + + 250 + + + + + + + + + + + + + 500 + + + + + + + + + + + + + 1k + + + + + + + + + + + + + 2k + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Frequency [Hz] + + + + + + + + + + + + + + + + + 0.0 + + + + + + + + + + + + + 2.5 + + + + + + + + + + + + + 5.0 + + + + + + + + + + + + + 7.5 + + + + + + + + + + + + + 10.0 + + + + + + + + + + + + + 12.5 + + + + + + + + + + + + + 15.0 + + + + + + + + + + + + + 17.5 + + + + + + + + + + + + + 20.0 + + + + Transmission loss [dB] + + + + + + + + + + + + + + + + + + + + + + + + + + + + Expansion-chamber transmission loss (Bies Eq. 8.111) + + + + + + + Area ratio m = Sexp/Sduct + + + + + + m = 2 → 1.9 dB + + + + + + m = 4 → 6.5 dB + + + + + + m = 8 → 12.2 dB + + + + + + m = 16 → 18.1 dB + + + + + + + + + + diff --git a/.github/images/silencer_expansion_chamber_es.svg b/.github/images/silencer_expansion_chamber_es.svg new file mode 100644 index 00000000..147741c7 --- /dev/null +++ b/.github/images/silencer_expansion_chamber_es.svg @@ -0,0 +1,1278 @@ + + + + + + + + image/svg+xml + + + Matplotlib v3.11.0, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 31,5 + + + + + + + + + + + + + 63 + + + + + + + + + + + + + 125 + + + + + + + + + + + + + 250 + + + + + + + + + + + + + 500 + + + + + + + + + + + + + 1k + + + + + + + + + + + + + 2k + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Frecuencia [Hz] + + + + + + + + + + + + + + + + + 0 + + + + + + + + + + + + + 2,5 + + + + + + + + + + + + + 5 + + + + + + + + + + + + + 7,5 + + + + + + + + + + + + + 10 + + + + + + + + + + + + + 12,5 + + + + + + + + + + + + + 15 + + + + + + + + + + + + + 17,5 + + + + + + + + + + + + + 20 + + + + Pérdida por transmisión [dB] + + + + + + + + + + + + + + + + + + + + + + + + + + + + Pérdida de transmisión de cámara de expansión (Bies Ec. 8,111) + + + + + + + Relación de áreas m = Sexp/Sduct + + + + + + m = 2 → 1,9 dB + + + + + + m = 4 → 6,5 dB + + + + + + m = 8 → 12,2 dB + + + + + + m = 16 → 18,1 dB + + + + + + + + + + diff --git a/.github/images/silencer_expansion_chamber_es_dark.svg b/.github/images/silencer_expansion_chamber_es_dark.svg new file mode 100644 index 00000000..c0431285 --- /dev/null +++ b/.github/images/silencer_expansion_chamber_es_dark.svg @@ -0,0 +1,1278 @@ + + + + + + + + image/svg+xml + + + Matplotlib v3.11.0, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 31,5 + + + + + + + + + + + + + 63 + + + + + + + + + + + + + 125 + + + + + + + + + + + + + 250 + + + + + + + + + + + + + 500 + + + + + + + + + + + + + 1k + + + + + + + + + + + + + 2k + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Frecuencia [Hz] + + + + + + + + + + + + + + + + + 0 + + + + + + + + + + + + + 2,5 + + + + + + + + + + + + + 5 + + + + + + + + + + + + + 7,5 + + + + + + + + + + + + + 10 + + + + + + + + + + + + + 12,5 + + + + + + + + + + + + + 15 + + + + + + + + + + + + + 17,5 + + + + + + + + + + + + + 20 + + + + Pérdida por transmisión [dB] + + + + + + + + + + + + + + + + + + + + + + + + + + + + Pérdida de transmisión de cámara de expansión (Bies Ec. 8,111) + + + + + + + Relación de áreas m = Sexp/Sduct + + + + + + m = 2 → 1,9 dB + + + + + + m = 4 → 6,5 dB + + + + + + m = 8 → 12,2 dB + + + + + + m = 16 → 18,1 dB + + + + + + + + + + diff --git a/docs/noise-control.md b/docs/noise-control.md new file mode 100644 index 00000000..9fb65023 --- /dev/null +++ b/docs/noise-control.md @@ -0,0 +1,173 @@ +← [Documentation index](README.md) + +# Industrial noise control: silencers, HVAC and enclosures + +Three passive measures dominate applied noise control, and the +`noise_control` domain covers all three with the engineering theory of Bies, +Hansen & Howard, *Engineering Noise Control* (5th ed., CRC Press 2017): +**reactive silencers** in a duct (the four-pole transmission-matrix method), +the passive attenuations and regenerated noise of an **HVAC** run, and the +insertion loss of a **machine enclosure**. The radiating piston of the +[electroacoustics](electroacoustics.md) domain is the companion radiator +model. + +## 1. Reactive silencers (four-pole method) + +A reactive silencer attenuates by *reflecting* sound with impedance +discontinuities. Each acoustic element is a 2×2 **transfer (four-pole) +matrix** relating the sound pressure `p` and the volume velocity `Su` at its +two ends (Bies Eq. (8.133); Munjal, *Acoustics of Ducts and Mufflers*), and a +compound silencer is the ordered matrix product of its elements. A straight +duct of length `L` and area `S` is (Bies Eq. (8.143), no flow) + +``` +[ cos(kL) j (rho c / S) sin(kL) ] +[ j (S / rho c) sin(kL) cos(kL) ] , k = omega / c, +``` + +and a side branch of acoustic impedance `Z_b` is the shunt +`[[1, 0], [1/Z_b, 1]]` (Eq. (8.144)). The **transmission loss** follows from +the compound matrix `T` (Bies Eq. (8.141); for equal inlet/outlet areas, +Eq. (8.148)) + +``` +TL = 20 log10( (1/2) | T11 + T12/Zc + Zc T21 + T22 | ) , Zc = rho c / S, +``` + +and the **insertion loss** for a source impedance `Z_s` and a radiation +impedance `Z_r` is the extra attenuation over a direct (zero-length) +connection, so a through connection gives `IL = 0`. + +### Expansion chamber + +A chamber of area `S_exp` and length `L` between pipes of area `S_duct` has the +closed-form transmission loss (Bies Eq. (8.111)) with area ratio +`m = S_exp / S_duct`: + +``` +TL = 10 log10[ 1 + (1/4) (m - 1/m)^2 sin^2(kL) ] , +``` + +peaking at `10 log10[1 + (1/4)(m - 1/m)^2]` at `kL = pi/2, 3 pi/2, ...` +(1.94 dB for `m = 2`, 6.55 dB for `m = 4`, 12.18 dB for `m = 8`, 18.10 dB for +`m = 16`) and dropping to 0 at `kL = n pi`, where the chamber is a +half-wavelength long and transparent. The four-pole product reproduces this +exactly. + +Expansion-chamber transmission loss against frequency for area ratios m = 2, 4, 8 and 16, showing periodic peaks rising with m at odd multiples of the quarter-wave frequency and troughs returning to 0 dB at every half-wavelength of the chamber length + +```python +import numpy as np +from phonometry import expansion_chamber + +freqs = np.linspace(20.0, 2000.0, 2000) +res = expansion_chamber(freqs, length=0.3, chamber_area=0.04, pipe_area=0.01) +print(round(res.transmission_loss.max(), 2)) # 6.55 dB peak (m = 4) +res.plot() # TL (and IL) vs frequency +``` + +### Side-branch and extended-tube resonators + +A **Helmholtz resonator** (`helmholtz_resonator`) and a closed **quarter-wave +tube** (`quarter_wave_resonator`) each short the duct at their tuning +frequency, `f_0 = (c / 2 pi) sqrt(S_neck / (l_e V))` (Bies Eq. (8.46)) and +`f = c / 4 l_e` (Eq. (8.44)), giving a sharp transmission-loss spike there. An +**extended-tube chamber** (`extended_tube_chamber`) buries quarter-wave side +branches in an expansion chamber to fill its troughs; with zero extensions it +reduces exactly to the plain chamber. Advanced layouts chain elements directly +with `duct_matrix`, `shunt_matrix`, `cascade`, `transmission_loss` and +`insertion_loss`. + +Each device returns a `ReactiveSilencerResult` with `transmission_loss`, +`insertion_loss` (when source/radiation impedances are given), the compound +`transfer_matrix`, the tuning `resonances` and `.plot()`. + +## 2. HVAC duct attenuation and flow noise + +`noise_control.hvac` gathers the Bies Chapter 8 duct methods: + +- `end_reflection_loss` — the low-frequency reflection back up an open duct end + (ASHRAE Table 8.14, interpolated over diameter and frequency; it passes + exactly through the tabulated nodes). +- `elbow_insertion_loss` — the insertion loss per bend for square/round, + vaned/unvaned and lined/unlined elbows keyed by `W / lambda` (ASHRAE + Table 8.11). +- `plenum_attenuation` — the plenum-chamber transmission loss by Wells' method + (Eq. (8.275)), whose reverberant term uses the plenum + [room constant](room-image-sources.md). +- `flow_noise_straight_duct`, `flow_noise_bend` — the flow-generated (self) + noise sound power of straight ducts and mitred bends (VDI 2081, Eqs. (8.251), + (8.254)). + +```python +from phonometry.noise_control import hvac + +bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0] +er = hvac.end_reflection_loss(bands, diameter=0.30, termination="flush") +tl = hvac.plenum_attenuation(0.1, 1.0, 20.0, 0.2) # Wells' method, dB +fn = hvac.flow_noise_straight_duct(bands, flow_velocity=10.0, area=0.04) +``` + +Rectangular ducts use the equivalent diameter `D = sqrt(4 S / pi)`. Bies 5th +ed. gives the duct end reflection only as the ASHRAE table (no closed form in +that edition); this module reproduces and interpolates it. + +## 3. Machine enclosures + +A sealed enclosure reduces the radiated noise by its panel transmission loss +`R`, minus a penalty `C` for the reverberant build-up inside the small, hard +cavity (Bies Eqs. (7.103), (7.111)): + +``` +IL = R - C , C = 10 log10( 0.3 + S_E / R_i ) , +``` + +with the external area `S_E` and the interior room constant +`R_i = S_i alpha_i / (1 - alpha_i)` (the same `room_constant` as the +steady-state room field). A hard interior wastes much of the panel `R`; lining +it drives `C` toward its floor `10 log10 0.3 = -5.2 dB`. + +**The panel transmission loss `R` is supplied by the caller** — measured, or +predicted by a panel model — as a per-band array or a callable of frequency. +This module never predicts `R` itself; it combines a given `R` with the +interior absorption. + +```python +import numpy as np +from phonometry import enclosure_insertion_loss + +bands = np.array([125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0]) +panel_R = np.array([18.0, 24.0, 30.0, 36.0, 42.0, 46.0]) # measured, dB +enc = enclosure_insertion_loss(panel_R, external_area=6.0, internal_area=5.0, + internal_absorption=0.3, frequencies=bands) +print(np.round(enc.insertion_loss, 1)) # net IL = R - C per band +enc.plot() +``` + +`enclosure_insertion_loss` returns an `EnclosureResult` with the panel +`panel_transmission_loss`, the interior `correction`, the net `insertion_loss`, +the interior `room_constant` and `.plot()`. + +## Cross-check against the FDTD solver + +The four-pole expansion chamber is cross-checked against the independent 2D +[FDTD wave solver](fdtd-simulation.md): a plane-wave duct that widens into a +chamber and narrows back transmits far less at the four-pole TL peak +(`kL = pi/2`) than at the transparent trough (`kL = pi`), and the measured +amplitude ratio reproduces the closed-form peak transmission loss to a fraction +of a decibel (test `tests/noise_control/test_fdtd_crosscheck.py`). + +## References + +- Bies, D. A., Hansen, C. H., & Howard, C. Q. (2017). *Engineering noise + control* (5th ed.). CRC Press. + [doi:10.1201/9781351228152](https://doi.org/10.1201/9781351228152). The + muffler four-pole method and expansion-chamber TL (§8.8–8.9), the HVAC duct + methods (§8.11–8.17) and the machine-enclosure noise reduction (§7.4). +- Munjal, M. L. (2014). *Acoustics of ducts and mufflers* (2nd ed.). Wiley. + [doi:10.1002/9781118443767](https://doi.org/10.1002/9781118443767). The + transfer-matrix formulation behind the element matrices. +- Vér, I. L., & Beranek, L. L. (2006). *Noise and vibration control + engineering* (2nd ed.). Wiley. + [doi:10.1002/9780470172568](https://doi.org/10.1002/9780470172568). The + companion treatment of mufflers, ducts and enclosures. diff --git a/scripts/generate_graphs.py b/scripts/generate_graphs.py index b0453ade..3114111b 100644 --- a/scripts/generate_graphs.py +++ b/scripts/generate_graphs.py @@ -65,6 +65,9 @@ "Aircraft Atmospheric Absorption (SAE ARP 5534)": "Absorción atmosférica aeronáutica (SAE ARP 5534)", "Attenuation [dB]": "Atenuación [dB]", + "Expansion-chamber transmission loss (Bies Eq. 8.111)": + "Pérdida de transmisión de cámara de expansión (Bies Ec. 8.111)", + "Area ratio m = Sexp/Sduct": "Relación de áreas m = Sexp/Sduct", "Noise-Power-Distance Curves (ECAC Doc 29)": "Curvas nivel-potencia-distancia (ECAC Doc 29)", "Aircraft Departure SEL Contour (ECAC Doc 29)": @@ -7454,6 +7457,36 @@ def generate_panel_insulation_concept(output_dir: str) -> None: plt.close() +def generate_silencer_expansion_chamber(output_dir: str) -> None: + """Expansion-chamber transmission loss for four area ratios (Bies 8.111).""" + print("Generating silencer_expansion_chamber.svg...") + from phonometry import expansion_chamber + + freqs = np.linspace(20.0, 2000.0, 2000) + pipe_area, length = 0.01, 0.3 + ratios = (2.0, 4.0, 8.0, 16.0) + colors = (COLOR_PRIMARY, COLOR_SECONDARY, COLOR_TERTIARY, "#9467bd") + + fig, ax = plt.subplots(figsize=(9.0, 5.2)) + for m, color in zip(ratios, colors): + res = expansion_chamber(freqs, length, m * pipe_area, pipe_area) + peak = 10.0 * np.log10(1.0 + 0.25 * (m - 1.0 / m) ** 2) + ax.plot(freqs, res.transmission_loss, color=color, lw=1.8, + label=f"m = {int(m)} → {peak:.1f} dB") + ax.set_xlabel("Frequency [Hz]") + ax.set_ylabel("Transmission loss [dB]") + ax.set_title("Expansion-chamber transmission loss (Bies Eq. 8.111)", + fontweight="bold", pad=10) + ax.set_xlim(20.0, 2000.0) + ax.set_ylim(0.0, 20.0) + format_frequency_axis(ax, 20.0, 2000.0) + ax.grid(True, which="both", alpha=0.4) + ax.legend(loc="upper right", fontsize="small", title="Area ratio m = Sexp/Sduct") + plt.tight_layout() + save_figure(output_dir, "silencer_expansion_chamber.svg") + plt.close() + + # Every documentation figure, in the order the sequential generator has always # produced them. This registry is the single source of truth for both the # sequential path (``generate_all``) and the parallel runner (``--jobs``), @@ -7633,6 +7666,7 @@ def generate_panel_insulation_concept(output_dir: str) -> None: generate_fdtd_simulation, # Theoretical panel sound insulation (single/double wall, radiation, slit). generate_panel_insulation_concept, + generate_silencer_expansion_chamber, ) diff --git a/site/astro.config.mjs b/site/astro.config.mjs index f2cc0d6f..5bf5407d 100644 --- a/site/astro.config.mjs +++ b/site/astro.config.mjs @@ -518,6 +518,7 @@ export default defineConfig({ 'guides/sound-power', 'guides/electroacoustics', 'guides/swept-sine-distortion', + 'guides/noise-control', 'guides/program-loudness', ], }, diff --git a/site/src/content/docs/es/guides/noise-control.mdx b/site/src/content/docs/es/guides/noise-control.mdx new file mode 100644 index 00000000..4f924fe4 --- /dev/null +++ b/site/src/content/docs/es/guides/noise-control.mdx @@ -0,0 +1,276 @@ +--- +title: "Control de ruido industrial: silenciadores, HVAC y cerramientos" +description: "Silenciadores reactivos por el método de matrices de cuatro polos (cámaras de expansión, resonadores de Helmholtz y de cuarto de onda, tubos extendidos), atenuación y ruido de flujo en conductos HVAC, y pérdida de inserción de cerramientos de máquina, según Bies, Hansen y Howard." +references: + - type: book + authors: ["Bies, D. A.", "Hansen, C. H.", "Howard, C. Q."] + year: 2017 + title: "Engineering noise control" + edition: "5th ed." + publisher: "CRC Press" + doi: "10.1201/9781351228152" + note: "El método de cuatro polos de silenciadores y la TL de cámara de expansión (§8.8-8.9), los métodos HVAC de conducto (§8.11-8.17) y la reducción de ruido del cerramiento de máquina (§7.4) de esta página." + - type: book + authors: ["Munjal, M. L."] + year: 2014 + title: "Acoustics of ducts and mufflers" + edition: "2nd ed." + publisher: "Wiley" + doi: "10.1002/9781118443767" + note: "La formulación por matrices de transferencia y las matrices de los elementos del §1." + - type: book + authors: ["Vér, I. L.", "Beranek, L. L."] + year: 2006 + title: "Noise and vibration control engineering: Principles and applications" + edition: "2nd ed." + publisher: "Wiley" + doi: "10.1002/9780470172568" + note: "El tratamiento complementario de silenciadores, conductos y cerramientos contrastado en esta página." +--- + +import ThemeImage from '../../../../components/ThemeImage.astro'; + +Tres medidas pasivas dominan el control de ruido aplicado: los **silenciadores** +en un conducto, las atenuaciones pasivas y el ruido autogenerado de una +instalación **HVAC**, y el **cerramiento de máquina**. `phonometry.noise_control` +cubre las tres con la teoría de ingeniería de Bies, Hansen y Howard. Los +silenciadores se construyen sobre el método unidimensional de cuatro polos +(matrices de transferencia), que encadena elementos acústicos reutilizables; los +métodos HVAC son las tablas ASHRAE y las formas cerradas de codos, reflexión de +extremo, plenums y ruido de flujo; y el modelo de cerramiento combina una pérdida +de transmisión del panel aportada por el usuario con la acumulación reverberante +en el interior. El [pistón](/phonometry/es/guides/electroacoustics/) radiante del +dominio de electroacústica es el modelo de radiador complementario. + +## 1. Silenciadores reactivos (método de cuatro polos) + +Un silenciador reactivo atenúa *reflejando* el sonido con discontinuidades de +impedancia. Cada elemento acústico es una matriz de transferencia 2x2 que +relaciona la presión `p` y la velocidad de volumen `Su` en sus dos extremos, y un +silenciador compuesto es el producto ordenado de las matrices de sus elementos +(Bies §8.9). Un conducto recto de longitud `L` y área `S` es + +$$ +\begin{bmatrix} \cos kL & j\,\tfrac{\rho c}{S}\sin kL \\[2pt] +j\,\tfrac{S}{\rho c}\sin kL & \cos kL \end{bmatrix}, +\qquad k = \omega/c, +$$ + +y una rama lateral de impedancia $Z_b$ es la derivación +$\left[\begin{smallmatrix} 1 & 0 \\ 1/Z_b & 1 \end{smallmatrix}\right]$. La +**pérdida de transmisión** se obtiene de la matriz compuesta $T$ (Bies Ec. (8.141); +para áreas de entrada y salida iguales $S$, Ec. (8.148)) + +$$ +\mathrm{TL} = 20\log_{10}\!\left(\tfrac{1}{2}\left|\,T_{11} ++ \tfrac{T_{12}}{Z_c} + Z_c\,T_{21} + T_{22}\right|\right), +\qquad Z_c = \frac{\rho c}{S}, +$$ + +y la **pérdida de inserción** para una impedancia de fuente $Z_s$ y una de +radiación $Z_r$ es la atenuación adicional frente a una conexión directa. + +### Cámara de expansión + +El silenciador más simple, una cámara de área $S_\text{exp}$ y longitud $L$ entre +tubos de área $S_\text{duct}$, tiene la pérdida de transmisión en forma cerrada +(Bies Ec. (8.111)) con la relación de áreas $m = S_\text{exp}/S_\text{duct}$ + +$$ +\mathrm{TL} = 10\log_{10}\!\left[1 + \tfrac{1}{4}\left(m - \tfrac{1}{m}\right)^2 +\sin^2 kL\right], +$$ + +con máximos en $10\log_{10}[1 + \tfrac14(m-1/m)^2]$ cuando $kL = \pi/2, 3\pi/2, \dots$ +y caídas a $0$ en $kL = n\pi$, donde la cámara mide media longitud de onda y es +transparente. El producto de cuatro polos lo reproduce exactamente. + +```python +import numpy as np +from phonometry import expansion_chamber + +freqs = np.linspace(20.0, 2000.0, 2000) +res = expansion_chamber(freqs, length=0.3, chamber_area=0.04, pipe_area=0.01) +print(round(res.transmission_loss.max(), 2)) # máximo 6.55 dB (m = 4) +# Las caídas en f = n c / 2L son exactamente 0 dB (sin disipación). +print(round(float(res.transmission_loss[np.argmin(res.transmission_loss)]), 6)) +``` + +La figura muestra la pérdida de transmisión para cuatro relaciones de áreas: un +desajuste $m$ mayor eleva todos los máximos, pero las caídas permanecen en 0 dB y +los máximos permanecen en las mismas frecuencias, fijadas solo por la longitud de +la cámara. + + + +
+Mostrar el código de esta figura + +```python +import matplotlib.pyplot as plt +import numpy as np +from phonometry import expansion_chamber + +freqs = np.linspace(20.0, 2000.0, 2000) + +# Una línea para una cámara: TL vs frecuencia (con la pérdida de inserción si +# se dan las impedancias de fuente y radiación). +expansion_chamber(freqs, 0.3, 0.04, 0.01, + source_impedance=4e4, radiation_impedance=5e3).plot() +plt.show() + +# A mano: la familia de relaciones de áreas de la figura de concepto. +fig, ax = plt.subplots() +for m in (2.0, 4.0, 8.0, 16.0): + res = expansion_chamber(freqs, 0.3, m * 0.01, 0.01) + ax.plot(freqs, res.transmission_loss, label=f"m = {int(m)}") +ax.set_xlabel("Frecuencia [Hz]"); ax.set_ylabel("Pérdida de transmisión [dB]") +ax.legend() +plt.show() +``` + +
+ +### Resonadores laterales y de tubo extendido + +Un **resonador de Helmholtz** (área de cuello $S_n$, longitud efectiva $l_e$, +volumen de cavidad $V$) y un **tubo de cuarto de onda** cerrado (longitud $l_e$) +cortocircuitan el conducto en su frecuencia de sintonía, dando un pico agudo de +pérdida de transmisión: $f_0 = \tfrac{c}{2\pi}\sqrt{S_n/(l_e V)}$ (Bies Ec. (8.46)) +y $f = c/4l_e$ (Ec. (8.44)). Una **cámara de tubo extendido** aloja ramas de +cuarto de onda dentro de una cámara de expansión para rellenar sus caídas. + +```python +import numpy as np +from phonometry import ( + helmholtz_resonator, quarter_wave_resonator, extended_tube_chamber, +) + +f = np.linspace(20.0, 600.0, 4000) + +hr = helmholtz_resonator(f, duct_area=0.01, neck_area=1e-4, + neck_length=0.02, cavity_volume=1e-3) +print(round(float(hr.resonances[0]), 1)) # frecuencia de sintonía, Hz + +qw = quarter_wave_resonator(f, duct_area=0.01, length=1.516, branch_area=2e-3, + speed_of_sound=343.24) +print(round(float(qw.resonances[0]), 1)) # 56.6 Hz (Bies Ejemplo 8.1) + +# Una extensión de entrada de L/4 rellena la primera caída de la cámara. +et = extended_tube_chamber(f, length=0.4, chamber_area=0.04, pipe_area=0.01, + inlet_extension=0.1) +``` + +Cada uno devuelve un `ReactiveSilencerResult` con `transmission_loss`, +`insertion_loss` (cuando se dan las impedancias), la matriz compuesta +`transfer_matrix` y `.plot()`. Las configuraciones avanzadas encadenan elementos +directamente con `duct_matrix`, `shunt_matrix`, `cascade`, `transmission_loss` e +`insertion_loss`. + +## 2. Atenuación y ruido de flujo en conductos HVAC + +Una instalación de ventilación atenúa el ruido del ventilador en los codos, en +el extremo abierto del conducto y en los plenums, y regenera ruido allí donde el +flujo se perturba. `phonometry.noise_control.hvac` reúne los métodos del +Capítulo 8 de Bies. + +```python +from phonometry.noise_control import hvac + +bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0] + +# Reflexión de baja frecuencia en el extremo abierto (ASHRAE Tabla 8.14). +er = hvac.end_reflection_loss(bands, diameter=0.30, termination="flush") + +# Pérdida de inserción de un codo cuadrado a 90 grados con revestimiento (Tabla 8.11). +el = hvac.elbow_insertion_loss(bands, width=0.3, bend_type="square", lined=True) + +# TL de plenum por el método de Wells (forma cerrada). +tl = hvac.plenum_attenuation(exit_area=0.1, line_of_sight=1.0, + wall_area=20.0, mean_absorption=0.2) +print(round(tl, 1)) # dB + +# Ruido autogenerado (de flujo) de un conducto recto (VDI 2081). +fn = hvac.flow_noise_straight_duct(bands, flow_velocity=10.0, area=0.04) +``` + +Los métodos de reflexión de extremo y de codos interpolan las tablas ASHRAE +(pasan exactamente por los nodos tabulados); los métodos de plenum (Wells) y de +ruido de flujo (VDI 2081) son formas cerradas. `end_reflection_loss`, +`elbow_insertion_loss`, `flow_noise_straight_duct` y `flow_noise_bend` devuelven +un `HvacSpectrumResult` (atenuación o nivel de potencia sonora regenerada) con +`.plot()`; `plenum_attenuation` devuelve la pérdida de transmisión directamente. +Los conductos rectangulares usan el diámetro equivalente $D = \sqrt{4S/\pi}$. + +:::note +Bies 5.ª ed. da la reflexión de extremo solo como la consulta de la Tabla 8.14 +de ASHRAE; no hay forma cerrada en esa edición. Este módulo reproduce e +interpola la tabla. +::: + +## 3. Cerramientos de máquina + +Un cerramiento sellado reduce el ruido radiado en la pérdida de transmisión del +panel $R$, menos una penalización $C$ por la acumulación reverberante en la +cavidad pequeña y dura (Bies Ecs. (7.103), (7.111)): + +$$ +\mathrm{IL} = R - C,\qquad C = 10\log_{10}\!\left(0.3 + \frac{S_E}{R_i}\right), +$$ + +con el área externa $S_E$ y la **constante de sala** interior +$R_i = S_i \alpha_i/(1-\alpha_i)$ (la misma +[`room_constant`](/phonometry/es/guides/room-image-sources/) del campo estacionario +de sala). Un interior duro desperdicia gran parte de $R$; revestirlo lleva $C$ +hacia su límite inferior $10\log_{10}0.3 = -5.2$ dB. + +**La pérdida de transmisión del panel $R$ la aportas tú** como un array por banda +(medido, o predicho por un modelo de panel), o como una función de la frecuencia. +Este módulo nunca predice $R$; combina un $R$ dado con la absorción interior. + +```python +import numpy as np +from phonometry import enclosure_insertion_loss + +bands = np.array([125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0]) +panel_R = np.array([18.0, 24.0, 30.0, 36.0, 42.0, 46.0]) # medido, dB + +enc = enclosure_insertion_loss( + panel_R, external_area=6.0, internal_area=5.0, + internal_absorption=0.3, frequencies=bands, +) +print(np.round(enc.insertion_loss, 1)) # IL neta = R - C por banda +enc.plot() # panel R, corrección C e IL +``` + +`enclosure_insertion_loss` devuelve un `EnclosureResult` con la +`panel_transmission_loss` del panel, la `correction` interior, la `insertion_loss` +neta, la `room_constant` interior y `.plot()`. + +## Contraste con el solucionador FDTD + +La cámara de expansión de cuatro polos se contrasta con el +[solucionador de ondas FDTD 2D](/phonometry/es/guides/fdtd-simulation/) +independiente: un conducto de onda plana que se ensancha en una cámara y vuelve a +estrecharse transmite mucho menos en el máximo de TL de cuatro polos ($kL = \pi/2$) +que en la caída transparente ($kL = \pi$), y la razón de amplitudes medida +reproduce la pérdida de transmisión máxima en forma cerrada con una fracción de +decibelio. + +## Véase también + +- [Electroacústica](/phonometry/es/guides/electroacoustics/): el pistón radiante + (impedancia de radiación y directividad), el modelo de radiador complementario. +- [Potencia sonora](/phonometry/es/guides/sound-power/): la `Lw` de fuente que + alimenta un conducto o un cerramiento. +- [Fuentes imagen y campo estacionario](/phonometry/es/guides/room-image-sources/): + la `room_constant` reutilizada por la corrección interior del cerramiento. +- [Simulación de ondas FDTD 2D](/phonometry/es/guides/fdtd-simulation/): el + solucionador independiente que contrasta la cámara de expansión. +- [Informe de conformidad](https://github.com/jmrplens/phonometry/blob/main/docs/CONFORMANCE.md): + las formas cerradas y los anclajes resueltos con los que se validan estas + implementaciones. +- Referencia de la API: + [`noise_control.silencers`](/phonometry/reference/api/noise_control/silencers/), + [`noise_control.hvac`](/phonometry/reference/api/noise_control/hvac/) y + [`noise_control.enclosures`](/phonometry/reference/api/noise_control/enclosures/). diff --git a/site/src/content/docs/guides/noise-control.mdx b/site/src/content/docs/guides/noise-control.mdx new file mode 100644 index 00000000..bf1701f3 --- /dev/null +++ b/site/src/content/docs/guides/noise-control.mdx @@ -0,0 +1,272 @@ +--- +title: "Industrial Noise Control: Silencers, HVAC and Enclosures" +description: "Reactive silencers by the four-pole transmission-matrix method (expansion chambers, Helmholtz and quarter-wave resonators, extended tubes), HVAC duct attenuation and flow noise, and machine-enclosure insertion loss, from Bies, Hansen & Howard." +references: + - type: book + authors: ["Bies, D. A.", "Hansen, C. H.", "Howard, C. Q."] + year: 2017 + title: "Engineering noise control" + edition: "5th ed." + publisher: "CRC Press" + doi: "10.1201/9781351228152" + note: "The muffler four-pole method and expansion-chamber TL (§8.8-8.9), the HVAC duct methods (§8.11-8.17) and the machine-enclosure noise reduction (§7.4) of this page." + - type: book + authors: ["Munjal, M. L."] + year: 2014 + title: "Acoustics of ducts and mufflers" + edition: "2nd ed." + publisher: "Wiley" + doi: "10.1002/9781118443767" + note: "The transfer-matrix formulation and the element matrices behind §1." + - type: book + authors: ["Vér, I. L.", "Beranek, L. L."] + year: 2006 + title: "Noise and vibration control engineering: Principles and applications" + edition: "2nd ed." + publisher: "Wiley" + doi: "10.1002/9780470172568" + note: "The companion treatment of mufflers, ducts and enclosures cross-checked in this page." +--- + +import ThemeImage from '../../../components/ThemeImage.astro'; + +Three passive measures dominate applied noise control: **silencers** in a duct, +the passive attenuations and regenerated noise of an **HVAC** run, and a +**machine enclosure**. `phonometry.noise_control` covers all three with the +engineering theory of Bies, Hansen & Howard. The silencers are built on the +one-dimensional four-pole (transfer-matrix) method, which chains reusable +acoustic elements; the HVAC methods are the ASHRAE tables and closed forms for +bends, end reflection, plenums and flow noise; and the enclosure model combines +a user-supplied panel transmission loss with the reverberant build-up inside +the enclosure. The radiating [piston](/phonometry/guides/electroacoustics/) of +the electroacoustics domain is the companion radiator model. + +## 1. Reactive silencers (four-pole method) + +A reactive silencer attenuates by *reflecting* sound with impedance +discontinuities. Each acoustic element is a 2x2 **transfer matrix** relating the +pressure `p` and volume velocity `Su` at its two ends, and a compound silencer +is the ordered matrix product of its elements (Bies §8.9). A straight duct of +length `L` and area `S` is + +$$ +\begin{bmatrix} \cos kL & j\,\tfrac{\rho c}{S}\sin kL \\[2pt] +j\,\tfrac{S}{\rho c}\sin kL & \cos kL \end{bmatrix}, +\qquad k = \omega/c, +$$ + +and a side branch of impedance $Z_b$ is the shunt +$\left[\begin{smallmatrix} 1 & 0 \\ 1/Z_b & 1 \end{smallmatrix}\right]$. The +**transmission loss** follows from the compound matrix $T$ (Bies Eq. (8.141); +for equal inlet/outlet areas $S$, Eq. (8.148)) + +$$ +\mathrm{TL} = 20\log_{10}\!\left(\tfrac{1}{2}\left|\,T_{11} ++ \tfrac{T_{12}}{Z_c} + Z_c\,T_{21} + T_{22}\right|\right), +\qquad Z_c = \frac{\rho c}{S}, +$$ + +and the **insertion loss** for a source impedance $Z_s$ and radiation impedance +$Z_r$ is the extra attenuation over a direct connection. + +### Expansion chamber + +The simplest silencer, a chamber of area $S_\text{exp}$ and length $L$ between +pipes of area $S_\text{duct}$, has the closed-form transmission loss (Bies +Eq. (8.111)) with area ratio $m = S_\text{exp}/S_\text{duct}$ + +$$ +\mathrm{TL} = 10\log_{10}\!\left[1 + \tfrac{1}{4}\left(m - \tfrac{1}{m}\right)^2 +\sin^2 kL\right], +$$ + +peaking at $10\log_{10}[1 + \tfrac14(m-1/m)^2]$ when $kL = \pi/2, 3\pi/2, \dots$ +and dropping to $0$ at $kL = n\pi$ where the chamber is a half-wavelength long +and transparent. The four-pole product reproduces this exactly. + +```python +import numpy as np +from phonometry import expansion_chamber + +freqs = np.linspace(20.0, 2000.0, 2000) +res = expansion_chamber(freqs, length=0.3, chamber_area=0.04, pipe_area=0.01) +print(round(res.transmission_loss.max(), 2)) # 6.55 dB peak (m = 4) +# The troughs at f = n c / 2L are exactly 0 dB (no dissipation). +print(round(float(res.transmission_loss[np.argmin(res.transmission_loss)]), 6)) +``` + +The figure below shows the transmission loss for four area ratios: a larger +mismatch $m$ lifts every peak, but the troughs stay at 0 dB and the peaks stay +at the same frequencies, set only by the chamber length. + + + +
+Show the code for this figure + +```python +import matplotlib.pyplot as plt +import numpy as np +from phonometry import expansion_chamber + +freqs = np.linspace(20.0, 2000.0, 2000) + +# One line for one chamber: TL vs frequency (with the insertion loss too if +# source/radiation impedances are given). +expansion_chamber(freqs, 0.3, 0.04, 0.01, + source_impedance=4e4, radiation_impedance=5e3).plot() +plt.show() + +# By hand: the family of area ratios of the concept figure. +fig, ax = plt.subplots() +for m in (2.0, 4.0, 8.0, 16.0): + res = expansion_chamber(freqs, 0.3, m * 0.01, 0.01) + ax.plot(freqs, res.transmission_loss, label=f"m = {int(m)}") +ax.set_xlabel("Frequency [Hz]"); ax.set_ylabel("Transmission loss [dB]") +ax.legend() +plt.show() +``` + +
+ +### Side-branch and extended-tube resonators + +A **Helmholtz resonator** (neck area $S_n$, effective length $l_e$, cavity +volume $V$) and a closed **quarter-wave tube** (length $l_e$) each short the +duct at their tuning frequency, giving a sharp transmission-loss spike there: +$f_0 = \tfrac{c}{2\pi}\sqrt{S_n/(l_e V)}$ (Bies Eq. (8.46)) and $f = c/4l_e$ +(Eq. (8.44)). An **extended-tube chamber** buries quarter-wave side branches in +an expansion chamber to fill the plain chamber's troughs. + +```python +import numpy as np +from phonometry import ( + helmholtz_resonator, quarter_wave_resonator, extended_tube_chamber, +) + +f = np.linspace(20.0, 600.0, 4000) + +hr = helmholtz_resonator(f, duct_area=0.01, neck_area=1e-4, + neck_length=0.02, cavity_volume=1e-3) +print(round(float(hr.resonances[0]), 1)) # tuning frequency, Hz + +qw = quarter_wave_resonator(f, duct_area=0.01, length=1.516, branch_area=2e-3, + speed_of_sound=343.24) +print(round(float(qw.resonances[0]), 1)) # 56.6 Hz (Bies Example 8.1) + +# An inlet extension of L/4 fills the first expansion-chamber trough. +et = extended_tube_chamber(f, length=0.4, chamber_area=0.04, pipe_area=0.01, + inlet_extension=0.1) +``` + +Each returns a `ReactiveSilencerResult` with `transmission_loss`, +`insertion_loss` (when impedances are given), the compound `transfer_matrix` +and `.plot()`. Advanced layouts chain elements directly with `duct_matrix`, +`shunt_matrix`, `cascade`, `transmission_loss` and `insertion_loss`. + +## 2. HVAC duct attenuation and flow noise + +A ventilation run attenuates fan noise at bends, at the open duct end and in +plenums, and regenerates noise wherever the airflow is disturbed. `phonometry.noise_control.hvac` +gathers the Bies Chapter 8 methods. + +```python +from phonometry.noise_control import hvac + +bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0] + +# Low-frequency reflection back up an open duct end (ASHRAE Table 8.14). +er = hvac.end_reflection_loss(bands, diameter=0.30, termination="flush") + +# Insertion loss of a lined square 90-degree elbow (ASHRAE Table 8.11). +el = hvac.elbow_insertion_loss(bands, width=0.3, bend_type="square", lined=True) + +# Plenum-chamber TL by Wells' method (closed form). +tl = hvac.plenum_attenuation(exit_area=0.1, line_of_sight=1.0, + wall_area=20.0, mean_absorption=0.2) +print(round(tl, 1)) # dB + +# Flow-generated (self) noise of a straight duct (VDI 2081). +fn = hvac.flow_noise_straight_duct(bands, flow_velocity=10.0, area=0.04) +``` + +The end-reflection and elbow methods interpolate the ASHRAE tables (they pass +exactly through the tabulated nodes); the plenum (Wells) and flow-noise (VDI +2081) methods are closed forms. `end_reflection_loss`, `elbow_insertion_loss`, +`flow_noise_straight_duct` and `flow_noise_bend` return an `HvacSpectrumResult` +(attenuation or regenerated sound power level) with `.plot()`; `plenum_attenuation` +returns the transmission loss directly. Rectangular ducts use the equivalent +diameter $D = \sqrt{4S/\pi}$. + +:::note +Bies 5th ed. gives the duct end reflection only as the ASHRAE Table 8.14 +look-up; there is no closed form in that edition. This module reproduces and +interpolates the table. +::: + +## 3. Machine enclosures + +A sealed enclosure reduces the radiated noise by its panel transmission loss +$R$, minus a penalty $C$ for the reverberant build-up inside the small, hard +cavity (Bies Eqs. (7.103), (7.111)): + +$$ +\mathrm{IL} = R - C,\qquad C = 10\log_{10}\!\left(0.3 + \frac{S_E}{R_i}\right), +$$ + +with the external area $S_E$ and the interior **room constant** +$R_i = S_i \alpha_i/(1-\alpha_i)$ (the same +[`room_constant`](/phonometry/guides/room-image-sources/) as the steady-state +room field). A hard interior wastes much of the panel $R$; lining it drives $C$ +toward its floor $10\log_{10}0.3 = -5.2$ dB. + +**The panel transmission loss $R$ is supplied by you** as a per-band array +(measured, or predicted by a panel model), or as a callable of frequency. This +module never predicts $R$ itself; it combines a given $R$ with the interior +absorption. + +```python +import numpy as np +from phonometry import enclosure_insertion_loss + +bands = np.array([125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0]) +panel_R = np.array([18.0, 24.0, 30.0, 36.0, 42.0, 46.0]) # measured, dB + +enc = enclosure_insertion_loss( + panel_R, external_area=6.0, internal_area=5.0, + internal_absorption=0.3, frequencies=bands, +) +print(np.round(enc.insertion_loss, 1)) # net IL = R - C per band +enc.plot() # panel R, correction C and IL +``` + +`enclosure_insertion_loss` returns an `EnclosureResult` with the panel +`panel_transmission_loss`, the interior `correction`, the net `insertion_loss`, +the interior `room_constant`, and `.plot()`. + +## Cross-check against the FDTD solver + +The four-pole expansion chamber is cross-checked against the independent 2D +[FDTD wave solver](/phonometry/guides/fdtd-simulation/): a plane-wave duct that +widens into a chamber and narrows back transmits far less at the four-pole TL +peak ($kL = \pi/2$) than at the transparent trough ($kL = \pi$), and the +measured amplitude ratio reproduces the closed-form peak transmission loss to a +fraction of a decibel. + +## See also + +- [Electroacoustics](/phonometry/guides/electroacoustics/): the radiating + piston (radiation impedance and directivity), the companion radiator model. +- [Sound Power](/phonometry/guides/sound-power/): the source `Lw` that feeds a + duct or an enclosure. +- [Room image sources and steady field](/phonometry/guides/room-image-sources/): + the `room_constant` reused by the enclosure interior correction. +- [2D FDTD wave simulation](/phonometry/guides/fdtd-simulation/): the + independent solver that cross-checks the expansion chamber. +- [Conformance report](https://github.com/jmrplens/phonometry/blob/main/docs/CONFORMANCE.md): + the closed forms and worked anchors these implementations are validated + against. +- API reference: + [`noise_control.silencers`](/phonometry/reference/api/noise_control/silencers/), + [`noise_control.hvac`](/phonometry/reference/api/noise_control/hvac/) and + [`noise_control.enclosures`](/phonometry/reference/api/noise_control/enclosures/). diff --git a/site/src/content/docs/reference/api/noise_control/enclosures.md b/site/src/content/docs/reference/api/noise_control/enclosures.md new file mode 100644 index 00000000..986ce7e9 --- /dev/null +++ b/site/src/content/docs/reference/api/noise_control/enclosures.md @@ -0,0 +1,106 @@ +--- +title: "noise_control.enclosures" +description: "Public API of phonometry.noise_control.enclosures (auto-generated)." +sidebar: + label: "enclosures" +--- + +> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand. + +Insertion loss of a close or free-standing machine enclosure. + +Wrapping a machine in a sealed enclosure reduces the radiated noise by the +transmission loss of its panels, *minus* a penalty for the reverberant build-up +inside the small, hard cavity. Bies, Hansen & Howard, *Engineering Noise +Control* 5th ed., §7.4.2 (Eqs. (7.103), (7.111)) write the net reduction as + + IL = R - C, C = 10 log10[ 0.3 + S_E (1 - alpha_i) / (S_i alpha_i) ], + +where `R` is the field-incidence transmission loss of the enclosure panels, +`S_E` the external surface area, `S_i` the internal surface area (including +the machine) and `alpha_i` the mean absorption of the enclosure interior. The +reverberant term is exactly `S_E` over the interior **room constant** +`R_i = S_i alpha_i / (1 - alpha_i)` ([`phonometry.room.room_constant`](/phonometry/reference/api/rooms/steady-field/#room_constant)), so + + C = 10 log10( 0.3 + S_E / R_i ). + +A hard interior (`alpha_i` small) makes `C` large and wastes much of the +panel `R`; lining the enclosure drives `C` toward its floor +`10 log10 0.3 = -5.2 dB` (a fully absorbing interior, where `IL = R + 5.2`). +Bies terms this net reduction the enclosure *noise reduction*; it is the +insertion loss of the enclosure. + +**The panel transmission loss `R` is supplied by the caller** -- measured, or +predicted by a panel model -- as a per-band array or a callable of frequency. +This module never predicts `R` itself; it combines a given `R` with the +interior absorption. The interior room constant reuses +[`phonometry.room.room_constant`](/phonometry/reference/api/rooms/steady-field/#room_constant). + +## enclosure_insertion_loss + +```python +enclosure_insertion_loss( + panel_transmission_loss: ArrayLike | Callable[[NDArray[np.float64]], ArrayLike], + external_area: float, + internal_area: float, + internal_absorption: ArrayLike, + *, + frequencies: ArrayLike | None = None, +) -> EnclosureResult +``` + +Net insertion loss of a machine enclosure (Bies Eqs. (7.103), (7.111)). + +`IL = R - C` with `C = 10 log10(0.3 + S_E / R_i)` and the interior room +constant `R_i = S_i alpha_i / (1 - alpha_i)`. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `panel_transmission_loss` | Panel transmission loss `R` per band, dB. Either a per-band array (measured or predicted elsewhere) or a callable mapping a frequency array to per-band `R` (then `frequencies` is required). This function does not predict `R`. | +| `external_area` | External enclosure surface area `S_E`, m2. | +| `internal_area` | Internal surface area `S_i` (including the machine), m2. | +| `internal_absorption` | Mean interior absorption `alpha_i` in `(0, 1)` (scalar or per-band). | +| `frequencies` | Band centre frequencies, Hz; required when `panel_transmission_loss` is a callable, optional otherwise (used to label the result and the plot). | + +**Returns:** An [`EnclosureResult`](/phonometry/reference/api/noise_control/enclosures/#enclosureresult). + +## EnclosureResult + +```python +EnclosureResult( + frequencies: np.ndarray | None, + panel_transmission_loss: np.ndarray, + correction: np.ndarray, + insertion_loss: np.ndarray, + external_area: float, + internal_area: float, + room_constant: np.ndarray, +) +``` + +Insertion loss of a machine enclosure over frequency (Bies §7.4.2). + +**Attributes** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz, or `None` if the panel `R` was given as a bare per-band array with no frequency labels. | +| `panel_transmission_loss` | The supplied panel transmission loss `R` per band, dB. | +| `correction` | The interior-build-up correction `C` per band, dB. | +| `insertion_loss` | The net enclosure insertion loss `IL = R - C`, dB. | +| `external_area` | External enclosure surface area `S_E`, m2. | +| `internal_area` | Internal surface area `S_i`, m2. | +| `room_constant` | Interior room constant `R_i` per band, m2. | + +### EnclosureResult.plot() + +```python +EnclosureResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes +``` + +Plot the panel `R`, correction `C` and net insertion loss. + +Requires matplotlib (`pip install phonometry[plot]`); returns the +`Axes`. diff --git a/site/src/content/docs/reference/api/noise_control/hvac.md b/site/src/content/docs/reference/api/noise_control/hvac.md new file mode 100644 index 00000000..55877be3 --- /dev/null +++ b/site/src/content/docs/reference/api/noise_control/hvac.md @@ -0,0 +1,217 @@ +--- +title: "noise_control.hvac" +description: "Public API of phonometry.noise_control.hvac (auto-generated)." +sidebar: + label: "hvac" +--- + +> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand. + +HVAC duct acoustics: end reflection, bends, plenums and flow-generated noise. + +A ventilation duct network attenuates fan noise through several mechanisms +that add up along the path, and it *regenerates* noise wherever the airflow is +disturbed. This module gathers the engineering methods of Bies, Hansen & +Howard, *Engineering Noise Control* 5th ed., Chapter 8, for the passive +attenuations -- **duct end reflection** (§8.13, Table 8.14), **bends/elbows** +(§8.11, Table 8.11) and **plenum chambers** (§8.17, Wells' method) -- and for +the **flow-generated (self) noise** of straight ducts and bends (§8.15). + +The end-reflection and elbow methods are empirical look-up tables (ASHRAE); +they are interpolated over the duct size and, for the elbows, over the +frequency-to-width ratio `W / lambda`. The plenum and flow-noise methods are +closed forms evaluated directly. + +:::note +Bies 5th ed. gives the duct end reflection only as the ASHRAE Table 8.14 +look-up (there is no closed form in this edition); this module reproduces +that table and interpolates it. Rectangular ducts use the equivalent +diameter `D = sqrt(4 S / pi)`. +::: + +## elbow_insertion_loss + +```python +elbow_insertion_loss( + frequencies: ArrayLike, + width: float, + *, + bend_type: str = 'square', + vanes: bool = False, + lined: bool = False, + speed_of_sound: float = 343.0, +) -> HvacSpectrumResult +``` + +Duct bend/elbow insertion loss per bend (Bies Table 8.11, ASHRAE). + +Indexed by the frequency-to-width ratio `W / lambda` (`lambda = c / f`). +Lined bends assume the lining extends at least three duct diameters up- and +downstream. Round bends are treated as unlined with no vanes. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz (1-D array). | +| `width` | Duct width `W` in the plane of the bend, m. | +| `bend_type` | `"square"` or `"round"`. | +| `vanes` | Turning vanes fitted (square bends only). | +| `lined` | Acoustically lined bend (square bends only). | +| `speed_of_sound` | Speed of sound `c`, m/s. | + +**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the insertion loss, dB per bend. + +## end_reflection_loss + +```python +end_reflection_loss( + frequencies: ArrayLike, + diameter: float, + *, + termination: str = 'flush', + speed_of_sound: float = 343.0, +) -> HvacSpectrumResult +``` + +Duct end reflection loss (Bies Table 8.14, ASHRAE). + +The low-frequency reflection of sound back up a duct at its open +termination into a room. Interpolated over `log` diameter and `log` +frequency from Table 8.14; it passes exactly through the tabulated +`(diameter, octave band)` nodes. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz (1-D array). | +| `diameter` | Duct internal diameter `D`, m (use `D = sqrt(4 S / pi)` for a rectangular duct of area `S`). | +| `termination` | `"flush"` (duct flush with a wall/ceiling) or `"free"` (free space / suspended in the room). | +| `speed_of_sound` | Speed of sound `c`, m/s (kept for signature symmetry; the table is indexed by frequency directly). | + +**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the reflection loss, dB. + +## flow_noise_bend + +```python +flow_noise_bend( + frequencies: ArrayLike, + flow_velocity: float, + area: float, + height: float, + *, + density: float = 1.206, +) -> HvacSpectrumResult +``` + +Flow-generated octave-band sound power of a mitred bend (Bies Eqs. (8.252), (8.254)). + +`L_WB = L_Ws - 10 log10(1 + 0.165 N_s^2) + 30 log10(U) - 103` with the +stream power level `L_Ws = 30 log10(U) + 10 log10(S) + 10 log10(rho) + 117` +and the Strouhal number `N_s = f H / U` (`H` the duct height in the +plane of the bend). The regeneration scales as `U^3` at low `N_s` and +`U^5` at high `N_s`. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Octave-band centre frequencies `f`, Hz (1-D array). | +| `flow_velocity` | Mean flow speed `U`, m/s. | +| `area` | Duct cross-sectional area `S`, m2. | +| `height` | Duct height `H` in the plane of the bend, m. | +| `density` | Air density `rho`, kg/m3. | + +**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the band sound power level, dB re 1e-12 W. + +## flow_noise_straight_duct + +```python +flow_noise_straight_duct( + frequencies: ArrayLike, + flow_velocity: float, + area: float, +) -> HvacSpectrumResult +``` + +Flow-generated octave-band sound power of a straight duct (Bies Eq. (8.251)). + +`L_WB = 7 + 50 log10(U) + 10 log10(S) - 26 log10(1.14 + 0.02 f / U)` in +dB re 1e-12 W (VDI 2081-1), for airflow speed `U` in a duct of area +`S`. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Octave-band centre frequencies `f`, Hz (1-D array). | +| `flow_velocity` | Mean flow speed `U`, m/s. | +| `area` | Duct cross-sectional area `S`, m2. | + +**Returns:** A [`HvacSpectrumResult`](/phonometry/reference/api/noise_control/hvac/#hvacspectrumresult) of the band sound power level, dB re 1e-12 W. + +## HvacSpectrumResult + +```python +HvacSpectrumResult( + frequencies: np.ndarray, + values: np.ndarray, + quantity: str, + label: str, +) +``` + +A per-frequency HVAC quantity (attenuation or regenerated power level). + +**Attributes** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz. | +| `values` | The quantity per frequency (dB, or dB re 1e-12 W for a sound power level). | +| `quantity` | What `values` holds (`"attenuation"` or `"sound_power_level"`). | +| `label` | A short human label of the element. | + +### HvacSpectrumResult.plot() + +```python +HvacSpectrumResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes +``` + +Plot the quantity against a continuous log-frequency axis. + +Requires matplotlib (`pip install phonometry[plot]`). + +## plenum_attenuation + +```python +plenum_attenuation( + exit_area: float, + line_of_sight: float, + wall_area: float, + mean_absorption: ArrayLike, + *, + angle: float = 0.0, +) -> np.ndarray | float +``` + +Plenum-chamber transmission loss by Wells' method (Bies Eq. (8.275)). + +`TL = -10 log10[ S_out ( cos(theta) / (pi r^2) + (1 - alpha) / (S_w alpha) ) ]`, +where the reverberant term uses the plenum room constant +`R = S_w alpha / (1 - alpha)` ([`phonometry.room.room_constant`](/phonometry/reference/api/rooms/steady-field/#room_constant)). The +method holds above the inlet cut-on and when the plenum is large compared +with the wavelength; it underpredicts the low-frequency loss by 5-10 dB. + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `exit_area` | Outlet-opening area `S_out`, m2. | +| `line_of_sight` | Straight-line inlet-to-outlet distance `r`, m. | +| `wall_area` | Total internal wall area `S_w`, m2. | +| `mean_absorption` | Mean Sabine wall absorption `alpha` in `(0, 1)` (scalar or per-band). | +| `angle` | Angle `theta` between the inlet axis and the line to the outlet, rad (default 0). | + +**Returns:** The transmission loss, dB (float for scalar absorption, else a per-band array). diff --git a/site/src/content/docs/reference/api/noise_control/silencers.md b/site/src/content/docs/reference/api/noise_control/silencers.md new file mode 100644 index 00000000..789a34d5 --- /dev/null +++ b/site/src/content/docs/reference/api/noise_control/silencers.md @@ -0,0 +1,242 @@ +--- +title: "noise_control.silencers" +description: "Public API of phonometry.noise_control.silencers (auto-generated)." +sidebar: + label: "silencers" +--- + +> Auto-generated from the source docstrings by `scripts/generate_api_docs.py` (`make api-docs`). Do not edit by hand. + +Reactive silencers by the four-pole (transmission-matrix) method. + +A reactive silencer controls noise by *reflecting* it back to the source with +impedance discontinuities -- sudden area changes and side branches -- rather +than by dissipating it in absorptive material. The one-dimensional plane-wave +theory represents each acoustic element by a 2x2 **transfer (four-pole) +matrix** relating the sound pressure `p` and volume velocity `S u` at its +two ends, and a compound silencer is the ordered matrix product of its +elements (Bies, Hansen & Howard, *Engineering Noise Control* 5th ed., §8.8-8.9; +Munjal, *Acoustics of Ducts and Mufflers*). + +**Transfer matrix** (Bies Eq. (8.133)), state vector `[p, S u]` with the +characteristic acoustic impedance `Z = rho c / S`. The plane-wave element for +a straight duct of length `L` and area `S` is (Bies Eq. (8.143), no flow) + + [[ cos(kL), j (rho c / S) sin(kL) ], + [ j (S / rho c) sin(kL), cos(kL) ]], k = omega / c, + +and a **side branch** of acoustic impedance `Z_b` is the shunt element +(Bies Eq. (8.144)) + + [[ 1, 0 ], + [ 1 / Z_b, 1 ]]. + +**Transmission loss** from the compound matrix `T` (Bies Eq. (8.141), no +flow; reduces to Eq. (8.148) for equal inlet/outlet areas): + + TL = 10 log10[ (Z1 / Zn) (1/4) | T11 + T12 / Zn + Z1 T21 + (Z1 / Zn) T22 |^2 ] + +with `Z1 = rho c / S_in` and `Zn = rho c / S_out`. `TL` is the intrinsic +attenuation for an anechoic termination. The **insertion loss** for a source of +internal impedance `Z_s` radiating into a termination impedance `Z_r` is +the extra attenuation of inserting the silencer in place of a direct +connection, + + IL = 20 log10 | (Z_s + Z_r) / (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) |, + +which is `0` when the silencer reduces to a through connection (`T = I`). + +**Simple expansion chamber.** A chamber of area `S_exp` and length `L` +between pipes of area `S_duct` has the closed-form transmission loss (Bies +Eq. (8.111)) with area ratio `m = S_exp / S_duct` + + TL = 10 log10[ 1 + (1/4) (m - 1/m)^2 sin^2(kL) ], + +peaking at `10 log10[1 + (1/4)(m - 1/m)^2]` when `kL = pi/2, 3pi/2, ...` and +dropping to `0` at `kL = n pi` (no dissipation). The four-pole product +reproduces this exactly, and the machinery extends to side-branch (Helmholtz, +quarter-wave) and extended-tube resonators that the closed form cannot cover. + +## expansion_chamber + +```python +expansion_chamber( + frequencies: ArrayLike, + length: float, + chamber_area: float, + pipe_area: float, + *, + speed_of_sound: float = 343.0, + density: float = 1.206, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult +``` + +Simple expansion-chamber silencer (Bies Eq. (8.111) / four-pole). + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz (1-D array). | +| `length` | Chamber length `L`, m. | +| `chamber_area` | Chamber cross-sectional area `S_exp`, m2. | +| `pipe_area` | Inlet/outlet pipe area `S_duct`, m2. | +| `speed_of_sound` | Speed of sound `c`, m/s. | +| `density` | Air density `rho`, kg/m3. | +| `source_impedance` | Optional source impedance `Z_s` for the insertion loss, Pa s/m3. | +| `radiation_impedance` | Optional radiation impedance `Z_r` for the insertion loss, Pa s/m3. | + +**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult) (its `transmission_loss` equals the closed form `10 log10[1 + (1/4)(m - 1/m)^2 sin^2(kL)]`). + +## extended_tube_chamber + +```python +extended_tube_chamber( + frequencies: ArrayLike, + length: float, + chamber_area: float, + pipe_area: float, + *, + inlet_extension: float = 0.0, + outlet_extension: float = 0.0, + speed_of_sound: float = 343.0, + density: float = 1.206, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult +``` + +Extended-inlet/outlet expansion chamber (Bies §8.9.7). + +The inlet and outlet pipes extend a distance into the chamber, forming +annular quarter-wave side branches (of area `S_exp - S_duct` and lengths +equal to the extensions, Bies Eq. (8.156)) at the two junctions. Tuning the +extensions (classically `L/4` and `L/2`) places quarter-wave peaks that +fill the `kL = n pi` troughs of the plain expansion chamber. With both +extensions `0` the result reduces exactly to [`expansion_chamber`](/phonometry/reference/api/noise_control/silencers/#expansion_chamber). + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz (1-D array). | +| `length` | Chamber length `L`, m. | +| `chamber_area` | Chamber cross-sectional area `S_exp`, m2. | +| `pipe_area` | Inlet/outlet pipe area `S_duct`, m2. | +| `inlet_extension` | Inlet pipe extension into the chamber `L_a`, m. | +| `outlet_extension` | Outlet pipe extension into the chamber `L_b`, m. | +| `speed_of_sound` | Speed of sound `c`, m/s. | +| `density` | Air density `rho`, kg/m3. | +| `source_impedance` | Optional source impedance `Z_s`, Pa s/m3. | +| `radiation_impedance` | Optional radiation impedance `Z_r`, Pa s/m3. | + +**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult). + +## helmholtz_resonator + +```python +helmholtz_resonator( + frequencies: ArrayLike, + duct_area: float, + neck_area: float, + neck_length: float, + cavity_volume: float, + *, + resistance: float = 0.0, + speed_of_sound: float = 343.0, + density: float = 1.206, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult +``` + +Side-branch Helmholtz resonator on a duct (Bies Eqs. (8.144), (8.152)). + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz (1-D array). | +| `duct_area` | Main-duct cross-sectional area `S_d`, m2. | +| `neck_area` | Resonator neck area `S_neck`, m2. | +| `neck_length` | Effective neck length `l_e`, m. | +| `cavity_volume` | Cavity volume `V`, m3. | +| `resistance` | Neck acoustic resistance `R`, Pa s/m3 (default 0). | +| `speed_of_sound` | Speed of sound `c`, m/s. | +| `density` | Air density `rho`, kg/m3. | +| `source_impedance` | Optional source impedance `Z_s`, Pa s/m3. | +| `radiation_impedance` | Optional radiation impedance `Z_r`, Pa s/m3. | + +**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult); `resonances` holds `f_0 = (c / 2 pi) sqrt(S_neck / (l_e V))`. + +## quarter_wave_resonator + +```python +quarter_wave_resonator( + frequencies: ArrayLike, + duct_area: float, + length: float, + branch_area: float, + *, + speed_of_sound: float = 343.0, + density: float = 1.206, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult +``` + +Closed quarter-wave side-branch tube on a duct (Bies Eqs. (8.144), (8.146)). + +**Parameters** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz (1-D array). | +| `duct_area` | Main-duct cross-sectional area `S_d`, m2. | +| `length` | Effective branch length `l_e`, m. | +| `branch_area` | Branch tube area `S`, m2. | +| `speed_of_sound` | Speed of sound `c`, m/s. | +| `density` | Air density `rho`, kg/m3. | +| `source_impedance` | Optional source impedance `Z_s`, Pa s/m3. | +| `radiation_impedance` | Optional radiation impedance `Z_r`, Pa s/m3. | + +**Returns:** A [`ReactiveSilencerResult`](/phonometry/reference/api/noise_control/silencers/#reactivesilencerresult); `resonances` holds the odd multiples of `f = c / (4 l_e)` within the frequency range. + +## ReactiveSilencerResult + +```python +ReactiveSilencerResult( + frequencies: np.ndarray, + transmission_loss: np.ndarray, + insertion_loss: np.ndarray | None, + transfer_matrix: np.ndarray, + kind: str, + resonances: np.ndarray | None = None, +) +``` + +Transmission and insertion loss of a reactive silencer over frequency. + +**Attributes** + +| Name | Description | +| :--- | :--- | +| `frequencies` | Frequencies `f`, Hz. | +| `transmission_loss` | Transmission loss per frequency, dB. | +| `insertion_loss` | Insertion loss per frequency, dB, or `None` when no source/radiation impedance was supplied. | +| `transfer_matrix` | The compound `(n_freq, 2, 2)` four-pole matrix. | +| `kind` | A short label of the device (e.g. `"expansion chamber"`). | +| `resonances` | Notable resonance frequencies, Hz (e.g. the resonator tuning frequency), or `None`. | + +### ReactiveSilencerResult.plot() + +```python +ReactiveSilencerResult.plot(ax: Axes | None = None, **kwargs: Any) -> Axes +``` + +Plot the transmission (and insertion) loss against frequency. + +Requires matplotlib (`pip install phonometry[plot]`); returns the +`Axes`. diff --git a/src/phonometry/_plot/noise_control.py b/src/phonometry/_plot/noise_control.py new file mode 100644 index 00000000..026bce27 --- /dev/null +++ b/src/phonometry/_plot/noise_control.py @@ -0,0 +1,126 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +"""Plot renderers for the noise_control domain (lazy imports from result .plot()).""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, Any + +import numpy as np + +from .common import ( + _C_PRIMARY, + _C_REFERENCE, + _C_SECONDARY, + _C_TERTIARY, + _new_axes, + format_frequency_axis, +) + +if TYPE_CHECKING: + from matplotlib.axes import Axes + + from ..noise_control.enclosures import EnclosureResult + from ..noise_control.hvac import HvacSpectrumResult + from ..noise_control.silencers import ReactiveSilencerResult + +_FREQ_LABEL = "Frequency [Hz]" + + +def plot_reactive_silencer( + result: "ReactiveSilencerResult", ax: "Axes | None" = None, **kwargs: Any +) -> "Axes": + """Transmission (and insertion) loss of a reactive silencer over frequency. + + :param result: A + :class:`~phonometry.noise_control.silencers.ReactiveSilencerResult`. + :param ax: Existing axes, or ``None`` to create a figure. + :param kwargs: Forwarded to the transmission-loss ``Axes.plot``. + :return: The axes. + """ + ax = ax if ax is not None else _new_axes() + f = np.asarray(result.frequencies, dtype=np.float64) + kwargs.setdefault("color", _C_PRIMARY) + kwargs.setdefault("label", "Transmission loss") + ax.plot(f, np.asarray(result.transmission_loss), lw=1.8, **kwargs) + if result.insertion_loss is not None: + ax.plot(f, np.asarray(result.insertion_loss), color=_C_SECONDARY, + lw=1.4, ls="--", label="Insertion loss") + if result.resonances is not None: + for i, fr in enumerate(np.atleast_1d(np.asarray(result.resonances))): + if f.min() <= fr <= f.max(): + ax.axvline(fr, color=_C_TERTIARY, ls=":", lw=1.0, + label="Resonance" if i == 0 else None) + ax.set_xlabel(_FREQ_LABEL) + ax.set_ylabel("Loss [dB]") + ax.set_title(f"Reactive silencer: {result.kind}") + ax.grid(True, which="both", alpha=0.3) + format_frequency_axis(ax) + ax.legend(loc="best", fontsize="small") + return ax + + +def plot_hvac_spectrum( + result: "HvacSpectrumResult", ax: "Axes | None" = None, **kwargs: Any +) -> "Axes": + """Per-frequency HVAC attenuation or regenerated sound power level. + + :param result: A :class:`~phonometry.noise_control.hvac.HvacSpectrumResult`. + :param ax: Existing axes, or ``None`` to create a figure. + :param kwargs: Forwarded to ``Axes.plot``. + :return: The axes. + """ + ax = ax if ax is not None else _new_axes() + f = np.asarray(result.frequencies, dtype=np.float64) + is_power = result.quantity == "sound_power_level" + kwargs.setdefault("color", _C_SECONDARY if is_power else _C_PRIMARY) + kwargs.setdefault("label", result.label) + ax.plot(f, np.asarray(result.values), lw=1.8, marker="o", ms=3, **kwargs) + ax.set_xlabel(_FREQ_LABEL) + ax.set_ylabel( + "Sound power level [dB re 1 pW]" if is_power else "Attenuation [dB]" + ) + ax.set_title(result.label) + ax.grid(True, which="both", alpha=0.3) + format_frequency_axis(ax) + ax.legend(loc="best", fontsize="small") + return ax + + +def plot_enclosure( + result: "EnclosureResult", ax: "Axes | None" = None, **kwargs: Any +) -> "Axes": + """Panel R, interior correction C and net insertion loss of an enclosure. + + :param result: An + :class:`~phonometry.noise_control.enclosures.EnclosureResult`. + :param ax: Existing axes, or ``None`` to create a figure. + :param kwargs: Forwarded to the insertion-loss ``Axes.plot``. + :return: The axes. + """ + ax = ax if ax is not None else _new_axes() + n = np.asarray(result.insertion_loss).size + if result.frequencies is not None: + x = np.asarray(result.frequencies, dtype=np.float64) + continuous = True + else: + x = np.arange(n, dtype=np.float64) + continuous = False + ax.plot(x, np.asarray(result.panel_transmission_loss), color=_C_REFERENCE, + lw=1.3, ls="--", marker="s", ms=3, label="Panel R") + ax.plot(x, np.asarray(result.correction), color=_C_TERTIARY, lw=1.3, + ls=":", marker="^", ms=3, label="Interior correction C") + kwargs.setdefault("color", _C_PRIMARY) + kwargs.setdefault("label", "Insertion loss (R - C)") + ax.plot(x, np.asarray(result.insertion_loss), lw=1.9, marker="o", ms=3, + **kwargs) + ax.set_ylabel("Level [dB]") + ax.set_title("Machine enclosure insertion loss") + ax.grid(True, which="both", alpha=0.3) + if continuous: + ax.set_xlabel(_FREQ_LABEL) + format_frequency_axis(ax) + else: + ax.set_xlabel("Band") + ax.set_xticks(x) + ax.legend(loc="best", fontsize="small") + return ax diff --git a/src/phonometry/noise_control/__init__.py b/src/phonometry/noise_control/__init__.py new file mode 100644 index 00000000..ba13487d --- /dev/null +++ b/src/phonometry/noise_control/__init__.py @@ -0,0 +1,51 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +"""noise_control domain of phonometry (see module docstrings).""" + +from __future__ import annotations + +from .enclosures import EnclosureResult, enclosure_insertion_loss +from .hvac import ( + HvacSpectrumResult, + elbow_insertion_loss, + end_reflection_loss, + flow_noise_bend, + flow_noise_straight_duct, + plenum_attenuation, +) +from .silencers import ( + ReactiveSilencerResult, + cascade, + duct_matrix, + expansion_chamber, + extended_tube_chamber, + helmholtz_impedance, + helmholtz_resonator, + insertion_loss, + quarter_wave_impedance, + quarter_wave_resonator, + shunt_matrix, + transmission_loss, +) + +__all__ = [ + "EnclosureResult", + "HvacSpectrumResult", + "ReactiveSilencerResult", + "cascade", + "duct_matrix", + "elbow_insertion_loss", + "enclosure_insertion_loss", + "end_reflection_loss", + "expansion_chamber", + "extended_tube_chamber", + "flow_noise_bend", + "flow_noise_straight_duct", + "helmholtz_impedance", + "helmholtz_resonator", + "insertion_loss", + "plenum_attenuation", + "quarter_wave_impedance", + "quarter_wave_resonator", + "shunt_matrix", + "transmission_loss", +] diff --git a/src/phonometry/noise_control/enclosures.py b/src/phonometry/noise_control/enclosures.py new file mode 100644 index 00000000..aaf06591 --- /dev/null +++ b/src/phonometry/noise_control/enclosures.py @@ -0,0 +1,155 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +""" +Insertion loss of a close or free-standing machine enclosure. + +Wrapping a machine in a sealed enclosure reduces the radiated noise by the +transmission loss of its panels, *minus* a penalty for the reverberant build-up +inside the small, hard cavity. Bies, Hansen & Howard, *Engineering Noise +Control* 5th ed., §7.4.2 (Eqs. (7.103), (7.111)) write the net reduction as + + IL = R - C, C = 10 log10[ 0.3 + S_E (1 - alpha_i) / (S_i alpha_i) ], + +where ``R`` is the field-incidence transmission loss of the enclosure panels, +``S_E`` the external surface area, ``S_i`` the internal surface area (including +the machine) and ``alpha_i`` the mean absorption of the enclosure interior. The +reverberant term is exactly ``S_E`` over the interior **room constant** +``R_i = S_i alpha_i / (1 - alpha_i)`` (:func:`phonometry.room.room_constant`), so + + C = 10 log10( 0.3 + S_E / R_i ). + +A hard interior (``alpha_i`` small) makes ``C`` large and wastes much of the +panel ``R``; lining the enclosure drives ``C`` toward its floor +``10 log10 0.3 = -5.2 dB`` (a fully absorbing interior, where ``IL = R + 5.2``). +Bies terms this net reduction the enclosure *noise reduction*; it is the +insertion loss of the enclosure. + +**The panel transmission loss ``R`` is supplied by the caller** -- measured, or +predicted by a panel model -- as a per-band array or a callable of frequency. +This module never predicts ``R`` itself; it combines a given ``R`` with the +interior absorption. The interior room constant reuses +:func:`phonometry.room.room_constant`. +""" + +from __future__ import annotations + +from collections.abc import Callable +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +from .._internal.validation import require_positive +from ..room.steady_field import room_constant + +if TYPE_CHECKING: + from matplotlib.axes import Axes + + +@dataclass(frozen=True) +class EnclosureResult: + """Insertion loss of a machine enclosure over frequency (Bies §7.4.2). + + :ivar frequencies: Frequencies ``f``, Hz, or ``None`` if the panel ``R`` + was given as a bare per-band array with no frequency labels. + :ivar panel_transmission_loss: The supplied panel transmission loss ``R`` + per band, dB. + :ivar correction: The interior-build-up correction ``C`` per band, dB. + :ivar insertion_loss: The net enclosure insertion loss ``IL = R - C``, dB. + :ivar external_area: External enclosure surface area ``S_E``, m2. + :ivar internal_area: Internal surface area ``S_i``, m2. + :ivar room_constant: Interior room constant ``R_i`` per band, m2. + """ + + frequencies: np.ndarray | None + panel_transmission_loss: np.ndarray + correction: np.ndarray + insertion_loss: np.ndarray + external_area: float + internal_area: float + room_constant: np.ndarray + + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": + """Plot the panel ``R``, correction ``C`` and net insertion loss. + + Requires matplotlib (``pip install phonometry[plot]``); returns the + :class:`~matplotlib.axes.Axes`. + """ + from .._plot.noise_control import plot_enclosure + + return plot_enclosure(self, ax=ax, **kwargs) + + +def enclosure_insertion_loss( + panel_transmission_loss: ArrayLike | Callable[[NDArray[np.float64]], ArrayLike], + external_area: float, + internal_area: float, + internal_absorption: ArrayLike, + *, + frequencies: ArrayLike | None = None, +) -> EnclosureResult: + """Net insertion loss of a machine enclosure (Bies Eqs. (7.103), (7.111)). + + ``IL = R - C`` with ``C = 10 log10(0.3 + S_E / R_i)`` and the interior room + constant ``R_i = S_i alpha_i / (1 - alpha_i)``. + + :param panel_transmission_loss: Panel transmission loss ``R`` per band, dB. + Either a per-band array (measured or predicted elsewhere) or a callable + mapping a frequency array to per-band ``R`` (then ``frequencies`` is + required). This function does not predict ``R``. + :param external_area: External enclosure surface area ``S_E``, m2. + :param internal_area: Internal surface area ``S_i`` (including the machine), + m2. + :param internal_absorption: Mean interior absorption ``alpha_i`` in + ``(0, 1)`` (scalar or per-band). + :param frequencies: Band centre frequencies, Hz; required when + ``panel_transmission_loss`` is a callable, optional otherwise (used to + label the result and the plot). + :return: An :class:`EnclosureResult`. + """ + s_e = require_positive(external_area, "external_area") + s_i = require_positive(internal_area, "internal_area") + + freqs: NDArray[np.float64] | None = None + if frequencies is not None: + freqs = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) + if freqs.ndim != 1 or freqs.size == 0: + raise ValueError("'frequencies' must be a non-empty 1-D array.") + if np.any(freqs <= 0.0) or not np.all(np.isfinite(freqs)): + raise ValueError("'frequencies' must be positive and finite.") + + if callable(panel_transmission_loss): + if freqs is None: + raise ValueError( + "'frequencies' is required when 'panel_transmission_loss' " + "is a callable." + ) + r = np.atleast_1d(np.asarray(panel_transmission_loss(freqs), dtype=np.float64)) + else: + r = np.atleast_1d(np.asarray(panel_transmission_loss, dtype=np.float64)) + if r.ndim != 1 or r.size == 0: + raise ValueError("'panel_transmission_loss' must be a non-empty 1-D array.") + if not np.all(np.isfinite(r)): + raise ValueError("'panel_transmission_loss' must be finite.") + + alpha = np.asarray(internal_absorption, dtype=np.float64) + if np.any(alpha <= 0.0) or np.any(alpha >= 1.0) or not np.all(np.isfinite(alpha)): + raise ValueError("'internal_absorption' must lie strictly in (0, 1).") + + r_i = np.atleast_1d(np.asarray(room_constant(s_i, alpha), dtype=np.float64)) + r_i_b, r_b = np.broadcast_arrays(r_i, r) + if freqs is not None and freqs.shape != r_b.shape: + raise ValueError( + "'frequencies' must match the number of panel-R / absorption bands." + ) + correction = 10.0 * np.log10(0.3 + s_e / r_i_b) + il = r_b - correction + return EnclosureResult( + frequencies=freqs, + panel_transmission_loss=np.array(r_b, dtype=np.float64), + correction=np.array(correction, dtype=np.float64), + insertion_loss=np.array(il, dtype=np.float64), + external_area=s_e, + internal_area=s_i, + room_constant=np.array(r_i_b, dtype=np.float64), + ) diff --git a/src/phonometry/noise_control/hvac.py b/src/phonometry/noise_control/hvac.py new file mode 100644 index 00000000..4db4c9f3 --- /dev/null +++ b/src/phonometry/noise_control/hvac.py @@ -0,0 +1,343 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +""" +HVAC duct acoustics: end reflection, bends, plenums and flow-generated noise. + +A ventilation duct network attenuates fan noise through several mechanisms +that add up along the path, and it *regenerates* noise wherever the airflow is +disturbed. This module gathers the engineering methods of Bies, Hansen & +Howard, *Engineering Noise Control* 5th ed., Chapter 8, for the passive +attenuations -- **duct end reflection** (§8.13, Table 8.14), **bends/elbows** +(§8.11, Table 8.11) and **plenum chambers** (§8.17, Wells' method) -- and for +the **flow-generated (self) noise** of straight ducts and bends (§8.15). + +The end-reflection and elbow methods are empirical look-up tables (ASHRAE); +they are interpolated over the duct size and, for the elbows, over the +frequency-to-width ratio ``W / lambda``. The plenum and flow-noise methods are +closed forms evaluated directly. + +.. note:: + Bies 5th ed. gives the duct end reflection only as the ASHRAE Table 8.14 + look-up (there is no closed form in this edition); this module reproduces + that table and interpolates it. Rectangular ducts use the equivalent + diameter ``D = sqrt(4 S / pi)``. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +from .._internal.validation import require_positive +from ..room.steady_field import room_constant + +if TYPE_CHECKING: + from matplotlib.axes import Axes + +_C_AIR = 343.0 + +# --------------------------------------------------------------------------- +# Bies Table 8.14 -- duct end reflection loss (dB), ASHRAE. +# Rows: internal diameter (mm). Columns: octave band centre (Hz). +# Two termination conditions: "flush" (duct flush with a wall/ceiling) and +# "free" (free space / suspended in the room). +# --------------------------------------------------------------------------- +_END_REFLECTION_BANDS: NDArray[np.float64] = np.array( + [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0] +) +_END_REFLECTION_DIAMETERS_MM: NDArray[np.float64] = np.array( + [150, 200, 250, 300, 400, 510, 610, 710, 810, 910, 1220, 1830], dtype=float +) +_END_REFLECTION_FLUSH: NDArray[np.float64] = np.array( + [ + [18, 12, 7, 3, 1, 0], + [15, 10, 5, 2, 1, 0], + [14, 8, 4, 1, 0, 0], + [12, 7, 3, 1, 0, 0], + [10, 5, 2, 1, 0, 0], + [8, 4, 1, 0, 0, 0], + [7, 3, 1, 0, 0, 0], + [6, 2, 1, 0, 0, 0], + [5, 2, 1, 0, 0, 0], + [4, 2, 0, 0, 0, 0], + [3, 1, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 0], + ], + dtype=float, +) +_END_REFLECTION_FREE: NDArray[np.float64] = np.array( + [ + [20, 14, 9, 5, 2, 1], + [18, 12, 7, 3, 1, 0], + [16, 11, 6, 2, 1, 0], + [14, 9, 5, 2, 1, 0], + [12, 7, 3, 1, 0, 0], + [10, 6, 2, 1, 0, 0], + [9, 5, 2, 1, 0, 0], + [8, 4, 1, 0, 0, 0], + [7, 3, 1, 0, 0, 0], + [6, 3, 1, 0, 0, 0], + [5, 2, 0, 0, 0, 0], + [3, 1, 0, 0, 0, 0], + ], + dtype=float, +) + +# --------------------------------------------------------------------------- +# Bies Table 8.11 -- elbow/bend insertion loss (dB per bend) vs W / lambda. +# The five columns are the supported (bend_type, vanes, lined, round) cases; +# each row is the value for the W/lambda band whose upper edge is the key. +# --------------------------------------------------------------------------- +_ELBOW_WL_UPPER: NDArray[np.float64] = np.array( + [0.14, 0.28, 0.55, 1.11, 2.22, np.inf] +) +_ELBOW_TABLE: dict[str, NDArray[np.float64]] = { + # square, no vanes, unlined + "square": np.array([0, 1, 5, 8, 4, 3], dtype=float), + # square, no vanes, lined + "square_lined": np.array([0, 1, 6, 11, 10, 10], dtype=float), + # square, with vanes, unlined + "square_vanes": np.array([0, 1, 4, 6, 4, 4], dtype=float), + # square, with vanes, lined + "square_vanes_lined": np.array([0, 1, 4, 7, 7, 7], dtype=float), + # round, no vanes, unlined + "round": np.array([0, 1, 2, 3, 3, 3], dtype=float), +} + + +def _frequencies(frequencies: ArrayLike) -> NDArray[np.float64]: + f = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) + if f.ndim != 1 or f.size == 0: + raise ValueError("'frequencies' must be a non-empty 1-D array.") + if np.any(f <= 0.0) or not np.all(np.isfinite(f)): + raise ValueError("'frequencies' must be positive and finite.") + return f + + +@dataclass(frozen=True) +class HvacSpectrumResult: + """A per-frequency HVAC quantity (attenuation or regenerated power level). + + :ivar frequencies: Frequencies ``f``, Hz. + :ivar values: The quantity per frequency (dB, or dB re 1e-12 W for a + sound power level). + :ivar quantity: What ``values`` holds (``"attenuation"`` or + ``"sound_power_level"``). + :ivar label: A short human label of the element. + """ + + frequencies: np.ndarray + values: np.ndarray + quantity: str + label: str + + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": + """Plot the quantity against a continuous log-frequency axis. + + Requires matplotlib (``pip install phonometry[plot]``). + """ + from .._plot.noise_control import plot_hvac_spectrum + + return plot_hvac_spectrum(self, ax=ax, **kwargs) + + +def end_reflection_loss( + frequencies: ArrayLike, + diameter: float, + *, + termination: str = "flush", + speed_of_sound: float = _C_AIR, +) -> HvacSpectrumResult: + """Duct end reflection loss (Bies Table 8.14, ASHRAE). + + The low-frequency reflection of sound back up a duct at its open + termination into a room. Interpolated over ``log`` diameter and ``log`` + frequency from Table 8.14; it passes exactly through the tabulated + ``(diameter, octave band)`` nodes. + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param diameter: Duct internal diameter ``D``, m (use + ``D = sqrt(4 S / pi)`` for a rectangular duct of area ``S``). + :param termination: ``"flush"`` (duct flush with a wall/ceiling) or + ``"free"`` (free space / suspended in the room). + :param speed_of_sound: Speed of sound ``c``, m/s (kept for signature + symmetry; the table is indexed by frequency directly). + :return: A :class:`HvacSpectrumResult` of the reflection loss, dB. + """ + f = _frequencies(frequencies) + d_mm = require_positive(diameter, "diameter") * 1000.0 + if termination == "flush": + table = _END_REFLECTION_FLUSH + elif termination == "free": + table = _END_REFLECTION_FREE + else: + raise ValueError("'termination' must be 'flush' or 'free'.") + require_positive(speed_of_sound, "speed_of_sound") + + log_d = np.log(_END_REFLECTION_DIAMETERS_MM) + log_f_band = np.log(_END_REFLECTION_BANDS) + # Interpolate the table in log-diameter for each band, then in log-freq. + per_band = np.array( + [np.interp(np.log(d_mm), log_d, table[:, j]) for j in range(table.shape[1])] + ) + values = np.interp(np.log(f), log_f_band, per_band) + return HvacSpectrumResult( + frequencies=f, values=values, quantity="attenuation", + label=f"End reflection ({termination}, D = {diameter * 1000:.0f} mm)", + ) + + +def elbow_insertion_loss( + frequencies: ArrayLike, + width: float, + *, + bend_type: str = "square", + vanes: bool = False, + lined: bool = False, + speed_of_sound: float = _C_AIR, +) -> HvacSpectrumResult: + """Duct bend/elbow insertion loss per bend (Bies Table 8.11, ASHRAE). + + Indexed by the frequency-to-width ratio ``W / lambda`` (``lambda = c / f``). + Lined bends assume the lining extends at least three duct diameters up- and + downstream. Round bends are treated as unlined with no vanes. + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param width: Duct width ``W`` in the plane of the bend, m. + :param bend_type: ``"square"`` or ``"round"``. + :param vanes: Turning vanes fitted (square bends only). + :param lined: Acoustically lined bend (square bends only). + :param speed_of_sound: Speed of sound ``c``, m/s. + :return: A :class:`HvacSpectrumResult` of the insertion loss, dB per bend. + """ + f = _frequencies(frequencies) + w = require_positive(width, "width") + c = require_positive(speed_of_sound, "speed_of_sound") + if bend_type == "round": + if vanes or lined: + raise ValueError("round bends take neither vanes nor lining.") + key = "round" + elif bend_type == "square": + key = "square" + ("_vanes" if vanes else "") + ("_lined" if lined else "") + else: + raise ValueError("'bend_type' must be 'square' or 'round'.") + col = _ELBOW_TABLE[key] + wl = w * f / c + idx = np.searchsorted(_ELBOW_WL_UPPER, wl, side="left") + idx = np.clip(idx, 0, col.size - 1) + values = col[idx] + return HvacSpectrumResult( + frequencies=f, values=values, quantity="attenuation", + label=f"Elbow ({key.replace('_', ', ')}, W = {w * 1000:.0f} mm)", + ) + + +def plenum_attenuation( + exit_area: float, + line_of_sight: float, + wall_area: float, + mean_absorption: ArrayLike, + *, + angle: float = 0.0, +) -> np.ndarray | float: + """Plenum-chamber transmission loss by Wells' method (Bies Eq. (8.275)). + + ``TL = -10 log10[ S_out ( cos(theta) / (pi r^2) + (1 - alpha) / (S_w alpha) ) ]``, + where the reverberant term uses the plenum room constant + ``R = S_w alpha / (1 - alpha)`` (:func:`phonometry.room.room_constant`). The + method holds above the inlet cut-on and when the plenum is large compared + with the wavelength; it underpredicts the low-frequency loss by 5-10 dB. + + :param exit_area: Outlet-opening area ``S_out``, m2. + :param line_of_sight: Straight-line inlet-to-outlet distance ``r``, m. + :param wall_area: Total internal wall area ``S_w``, m2. + :param mean_absorption: Mean Sabine wall absorption ``alpha`` in ``(0, 1)`` + (scalar or per-band). + :param angle: Angle ``theta`` between the inlet axis and the line to the + outlet, rad (default 0). + :return: The transmission loss, dB (float for scalar absorption, else a + per-band array). + """ + s_out = require_positive(exit_area, "exit_area") + r = require_positive(line_of_sight, "line_of_sight") + s_w = require_positive(wall_area, "wall_area") + alpha = np.asarray(mean_absorption, dtype=np.float64) + if np.any(alpha <= 0.0) or np.any(alpha >= 1.0) or not np.all(np.isfinite(alpha)): + raise ValueError("'mean_absorption' must lie strictly in (0, 1).") + r_const = np.asarray(room_constant(s_w, alpha), dtype=np.float64) + direct = np.cos(angle) / (np.pi * r**2) + reverberant = 1.0 / r_const + tl = -10.0 * np.log10(s_out * (direct + reverberant)) + return float(tl) if tl.ndim == 0 else tl + + +def flow_noise_straight_duct( + frequencies: ArrayLike, + flow_velocity: float, + area: float, +) -> HvacSpectrumResult: + """Flow-generated octave-band sound power of a straight duct (Bies Eq. (8.251)). + + ``L_WB = 7 + 50 log10(U) + 10 log10(S) - 26 log10(1.14 + 0.02 f / U)`` in + dB re 1e-12 W (VDI 2081-1), for airflow speed ``U`` in a duct of area + ``S``. + + :param frequencies: Octave-band centre frequencies ``f``, Hz (1-D array). + :param flow_velocity: Mean flow speed ``U``, m/s. + :param area: Duct cross-sectional area ``S``, m2. + :return: A :class:`HvacSpectrumResult` of the band sound power level, + dB re 1e-12 W. + """ + f = _frequencies(frequencies) + u = require_positive(flow_velocity, "flow_velocity") + s = require_positive(area, "area") + lw = ( + 7.0 + + 50.0 * np.log10(u) + + 10.0 * np.log10(s) + - 26.0 * np.log10(1.14 + 0.02 * f / u) + ) + return HvacSpectrumResult( + frequencies=f, values=lw, quantity="sound_power_level", + label=f"Straight-duct flow noise (U = {u:.1f} m/s)", + ) + + +def flow_noise_bend( + frequencies: ArrayLike, + flow_velocity: float, + area: float, + height: float, + *, + density: float = 1.206, +) -> HvacSpectrumResult: + """Flow-generated octave-band sound power of a mitred bend (Bies Eqs. (8.252), (8.254)). + + ``L_WB = L_Ws - 10 log10(1 + 0.165 N_s^2) + 30 log10(U) - 103`` with the + stream power level ``L_Ws = 30 log10(U) + 10 log10(S) + 10 log10(rho) + 117`` + and the Strouhal number ``N_s = f H / U`` (``H`` the duct height in the + plane of the bend). The regeneration scales as ``U^3`` at low ``N_s`` and + ``U^5`` at high ``N_s``. + + :param frequencies: Octave-band centre frequencies ``f``, Hz (1-D array). + :param flow_velocity: Mean flow speed ``U``, m/s. + :param area: Duct cross-sectional area ``S``, m2. + :param height: Duct height ``H`` in the plane of the bend, m. + :param density: Air density ``rho``, kg/m3. + :return: A :class:`HvacSpectrumResult` of the band sound power level, + dB re 1e-12 W. + """ + f = _frequencies(frequencies) + u = require_positive(flow_velocity, "flow_velocity") + s = require_positive(area, "area") + h = require_positive(height, "height") + rho = require_positive(density, "density") + lws = 30.0 * np.log10(u) + 10.0 * np.log10(s) + 10.0 * np.log10(rho) + 117.0 + ns = f * h / u + lw = lws - 10.0 * np.log10(1.0 + 0.165 * ns**2) + 30.0 * np.log10(u) - 103.0 + return HvacSpectrumResult( + frequencies=f, values=lw, quantity="sound_power_level", + label=f"Mitred-bend flow noise (U = {u:.1f} m/s)", + ) diff --git a/src/phonometry/noise_control/silencers.py b/src/phonometry/noise_control/silencers.py new file mode 100644 index 00000000..9dee2942 --- /dev/null +++ b/src/phonometry/noise_control/silencers.py @@ -0,0 +1,564 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +""" +Reactive silencers by the four-pole (transmission-matrix) method. + +A reactive silencer controls noise by *reflecting* it back to the source with +impedance discontinuities -- sudden area changes and side branches -- rather +than by dissipating it in absorptive material. The one-dimensional plane-wave +theory represents each acoustic element by a 2x2 **transfer (four-pole) +matrix** relating the sound pressure ``p`` and volume velocity ``S u`` at its +two ends, and a compound silencer is the ordered matrix product of its +elements (Bies, Hansen & Howard, *Engineering Noise Control* 5th ed., §8.8-8.9; +Munjal, *Acoustics of Ducts and Mufflers*). + +**Transfer matrix** (Bies Eq. (8.133)), state vector ``[p, S u]`` with the +characteristic acoustic impedance ``Z = rho c / S``. The plane-wave element for +a straight duct of length ``L`` and area ``S`` is (Bies Eq. (8.143), no flow) + + [[ cos(kL), j (rho c / S) sin(kL) ], + [ j (S / rho c) sin(kL), cos(kL) ]], k = omega / c, + +and a **side branch** of acoustic impedance ``Z_b`` is the shunt element +(Bies Eq. (8.144)) + + [[ 1, 0 ], + [ 1 / Z_b, 1 ]]. + +**Transmission loss** from the compound matrix ``T`` (Bies Eq. (8.141), no +flow; reduces to Eq. (8.148) for equal inlet/outlet areas): + + TL = 10 log10[ (Z1 / Zn) (1/4) | T11 + T12 / Zn + Z1 T21 + (Z1 / Zn) T22 |^2 ] + +with ``Z1 = rho c / S_in`` and ``Zn = rho c / S_out``. ``TL`` is the intrinsic +attenuation for an anechoic termination. The **insertion loss** for a source of +internal impedance ``Z_s`` radiating into a termination impedance ``Z_r`` is +the extra attenuation of inserting the silencer in place of a direct +connection, + + IL = 20 log10 | (Z_s + Z_r) / (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) |, + +which is ``0`` when the silencer reduces to a through connection (``T = I``). + +**Simple expansion chamber.** A chamber of area ``S_exp`` and length ``L`` +between pipes of area ``S_duct`` has the closed-form transmission loss (Bies +Eq. (8.111)) with area ratio ``m = S_exp / S_duct`` + + TL = 10 log10[ 1 + (1/4) (m - 1/m)^2 sin^2(kL) ], + +peaking at ``10 log10[1 + (1/4)(m - 1/m)^2]`` when ``kL = pi/2, 3pi/2, ...`` and +dropping to ``0`` at ``kL = n pi`` (no dissipation). The four-pole product +reproduces this exactly, and the machinery extends to side-branch (Helmholtz, +quarter-wave) and extended-tube resonators that the closed form cannot cover. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +from .._internal.validation import require_non_negative, require_positive + +if TYPE_CHECKING: + from matplotlib.axes import Axes + +#: Reference air properties at 20 degC, 101.325 kPa. +_C_AIR = 343.0 +_RHO_AIR = 1.206 + +_Complex = NDArray[np.complex128] + + +def _frequencies(frequencies: ArrayLike) -> NDArray[np.float64]: + """Validate a strictly positive, finite 1-D frequency grid (Hz).""" + f = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) + if f.ndim != 1 or f.size == 0: + raise ValueError("'frequencies' must be a non-empty 1-D array.") + if np.any(f <= 0.0) or not np.all(np.isfinite(f)): + raise ValueError("'frequencies' must be positive and finite.") + return f + + +def duct_matrix( + frequencies: ArrayLike, + length: float, + area: float, + *, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, +) -> _Complex: + """Four-pole matrix of a straight duct (Bies Eq. (8.143), no flow). + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param length: Duct length ``L``, m. + :param area: Cross-sectional area ``S``, m2. + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :return: A ``(n_freq, 2, 2)`` complex transfer-matrix array. + """ + f = _frequencies(frequencies) + length = require_non_negative(length, "length") + area = require_positive(area, "area") + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + k = 2.0 * np.pi * f / c + z = rho * c / area + cos = np.cos(k * length) + sin = np.sin(k * length) + t = np.empty((f.size, 2, 2), dtype=np.complex128) + t[:, 0, 0] = cos + t[:, 0, 1] = 1j * z * sin + t[:, 1, 0] = 1j * sin / z + t[:, 1, 1] = cos + return t + + +def shunt_matrix(branch_impedance: ArrayLike) -> _Complex: + """Four-pole matrix of a side branch of impedance ``Z_b`` (Bies Eq. (8.144)). + + :param branch_impedance: Acoustic impedance ``Z_b`` of the branch, + Pa s/m3 (1-D complex array over frequency). + :return: A ``(n_freq, 2, 2)`` complex transfer-matrix array. + """ + zb = np.atleast_1d(np.asarray(branch_impedance, dtype=np.complex128)) + if zb.ndim != 1 or zb.size == 0: + raise ValueError("'branch_impedance' must be a non-empty 1-D array.") + t = np.zeros((zb.size, 2, 2), dtype=np.complex128) + t[:, 0, 0] = 1.0 + t[:, 1, 0] = 1.0 / zb + t[:, 1, 1] = 1.0 + return t + + +def cascade(*matrices: _Complex) -> _Complex: + """Cascade element four-pole matrices from inlet to outlet. + + The compound matrix is the ordered product ``T1 @ T2 @ ... @ Tn`` (the + state at the inlet equals the compound matrix times the state at the + outlet), broadcast over the frequency axis. + + :param matrices: One or more ``(n_freq, 2, 2)`` arrays sharing ``n_freq``. + :return: The compound ``(n_freq, 2, 2)`` array. + """ + if not matrices: + raise ValueError("cascade() needs at least one matrix.") + total = matrices[0] + for m in matrices[1:]: + total = np.matmul(total, m) + return total + + +def transmission_loss( + transfer_matrix: _Complex, + *, + inlet_area: float, + outlet_area: float, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, +) -> NDArray[np.float64]: + """Transmission loss of a four-pole element (Bies Eq. (8.141), no flow). + + :param transfer_matrix: A ``(n_freq, 2, 2)`` compound matrix. + :param inlet_area: Inlet pipe area ``S_in``, m2. + :param outlet_area: Outlet pipe area ``S_out``, m2. + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :return: The transmission loss per frequency, dB. + """ + s_in = require_positive(inlet_area, "inlet_area") + s_out = require_positive(outlet_area, "outlet_area") + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + z1 = rho * c / s_in + zn = rho * c / s_out + t11 = transfer_matrix[:, 0, 0] + t12 = transfer_matrix[:, 0, 1] + t21 = transfer_matrix[:, 1, 0] + t22 = transfer_matrix[:, 1, 1] + term = t11 + t12 / zn + z1 * t21 + (z1 / zn) * t22 + return np.asarray( + 10.0 * np.log10((z1 / zn) * 0.25 * np.abs(term) ** 2), + dtype=np.float64, + ) + + +def insertion_loss( + transfer_matrix: _Complex, + *, + source_impedance: ArrayLike, + radiation_impedance: ArrayLike, +) -> NDArray[np.float64]: + """Insertion loss of a four-pole element for given end impedances. + + The attenuation from inserting the element in place of a direct (zero + length) connection between a source of internal impedance ``Z_s`` and a + radiation (termination) impedance ``Z_r``: + + IL = 20 log10 | (Z_s + Z_r) / + (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) |. + + :param transfer_matrix: A ``(n_freq, 2, 2)`` compound matrix. + :param source_impedance: Source internal acoustic impedance ``Z_s``, + Pa s/m3 (scalar or per-frequency, real or complex). + :param radiation_impedance: Termination/radiation acoustic impedance + ``Z_r``, Pa s/m3 (scalar or per-frequency). + :return: The insertion loss per frequency, dB. + """ + n = transfer_matrix.shape[0] + zs = np.broadcast_to( + np.asarray(source_impedance, dtype=np.complex128), (n,) + ) + zr = np.broadcast_to( + np.asarray(radiation_impedance, dtype=np.complex128), (n,) + ) + t11 = transfer_matrix[:, 0, 0] + t12 = transfer_matrix[:, 0, 1] + t21 = transfer_matrix[:, 1, 0] + t22 = transfer_matrix[:, 1, 1] + denom = t11 * zr + t12 + zs * zr * t21 + zs * t22 + return np.asarray( + 20.0 * np.log10(np.abs((zs + zr) / denom)), dtype=np.float64 + ) + + +def helmholtz_impedance( + frequencies: ArrayLike, + neck_area: float, + neck_length: float, + cavity_volume: float, + *, + resistance: float = 0.0, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, +) -> _Complex: + """Acoustic impedance of a Helmholtz side branch (Bies Eq. (8.152)). + + ``Z = R + j( rho omega l_e / S_neck - rho c^2 / (omega V) )`` with acoustic + mass ``rho l_e / S_neck`` and compliance ``V / (rho c^2)``; the resonance + ``f_0 = (c / 2 pi) sqrt(S_neck / (l_e V))`` (Bies Eq. (8.46)) is where the + reactance vanishes and the branch shorts the duct. + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param neck_area: Neck cross-sectional area ``S_neck``, m2. + :param neck_length: Effective neck length ``l_e`` (with end corrections), m. + :param cavity_volume: Cavity volume ``V``, m3. + :param resistance: Acoustic resistance ``R``, Pa s/m3 (default 0, lossless). + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :return: The complex branch impedance per frequency, Pa s/m3. + """ + f = _frequencies(frequencies) + s_neck = require_positive(neck_area, "neck_area") + le = require_positive(neck_length, "neck_length") + vol = require_positive(cavity_volume, "cavity_volume") + r = require_non_negative(resistance, "resistance") + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + omega = 2.0 * np.pi * f + reactance = rho * omega * le / s_neck - rho * c**2 / (omega * vol) + return np.asarray(r + 1j * reactance, dtype=np.complex128) + + +def quarter_wave_impedance( + frequencies: ArrayLike, + length: float, + area: float, + *, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, +) -> _Complex: + """Acoustic impedance of a closed quarter-wave side branch (Bies Eq. (8.146)). + + ``Z = -j (rho c / S) cot(k l_e)``; the reactance vanishes at + ``l_e = lambda / 4`` (``f = c / 4 l_e``), where the closed tube presents a + pressure node and shorts the duct. + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param length: Effective tube length ``l_e`` (with end correction), m. + :param area: Tube cross-sectional area ``S``, m2. + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :return: The complex branch impedance per frequency, Pa s/m3. + """ + f = _frequencies(frequencies) + le = require_positive(length, "length") + area = require_positive(area, "area") + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + k = 2.0 * np.pi * f / c + z = rho * c / area + # At the half-wave frequencies (k l_e = n pi) the closed tube is transparent + # (Z -> infinity); the divide-by-zero is expected and its shunt 1/Z is 0. + with np.errstate(divide="ignore", invalid="ignore"): + zb = -1j * z / np.tan(k * le) + return np.asarray(zb, dtype=np.complex128) + + +@dataclass(frozen=True) +class ReactiveSilencerResult: + """Transmission and insertion loss of a reactive silencer over frequency. + + :ivar frequencies: Frequencies ``f``, Hz. + :ivar transmission_loss: Transmission loss per frequency, dB. + :ivar insertion_loss: Insertion loss per frequency, dB, or ``None`` when no + source/radiation impedance was supplied. + :ivar transfer_matrix: The compound ``(n_freq, 2, 2)`` four-pole matrix. + :ivar kind: A short label of the device (e.g. ``"expansion chamber"``). + :ivar resonances: Notable resonance frequencies, Hz (e.g. the resonator + tuning frequency), or ``None``. + """ + + frequencies: np.ndarray + transmission_loss: np.ndarray + insertion_loss: np.ndarray | None + transfer_matrix: np.ndarray + kind: str + resonances: np.ndarray | None = None + + def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": + """Plot the transmission (and insertion) loss against frequency. + + Requires matplotlib (``pip install phonometry[plot]``); returns the + :class:`~matplotlib.axes.Axes`. + """ + from .._plot.noise_control import plot_reactive_silencer + + return plot_reactive_silencer(self, ax=ax, **kwargs) + + +def _result( + f: NDArray[np.float64], + t: _Complex, + *, + inlet_area: float, + outlet_area: float, + c: float, + rho: float, + source_impedance: ArrayLike | None, + radiation_impedance: ArrayLike | None, + kind: str, + resonances: NDArray[np.float64] | None = None, +) -> ReactiveSilencerResult: + """Assemble a :class:`ReactiveSilencerResult` from a compound matrix.""" + tl = transmission_loss( + t, inlet_area=inlet_area, outlet_area=outlet_area, + speed_of_sound=c, density=rho, + ) + il: NDArray[np.float64] | None = None + if source_impedance is not None and radiation_impedance is not None: + il = insertion_loss( + t, source_impedance=source_impedance, + radiation_impedance=radiation_impedance, + ) + return ReactiveSilencerResult( + frequencies=f, + transmission_loss=tl, + insertion_loss=il, + transfer_matrix=t, + kind=kind, + resonances=resonances, + ) + + +def expansion_chamber( + frequencies: ArrayLike, + length: float, + chamber_area: float, + pipe_area: float, + *, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult: + """Simple expansion-chamber silencer (Bies Eq. (8.111) / four-pole). + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param length: Chamber length ``L``, m. + :param chamber_area: Chamber cross-sectional area ``S_exp``, m2. + :param pipe_area: Inlet/outlet pipe area ``S_duct``, m2. + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :param source_impedance: Optional source impedance ``Z_s`` for the + insertion loss, Pa s/m3. + :param radiation_impedance: Optional radiation impedance ``Z_r`` for the + insertion loss, Pa s/m3. + :return: A :class:`ReactiveSilencerResult` (its ``transmission_loss`` + equals the closed form ``10 log10[1 + (1/4)(m - 1/m)^2 sin^2(kL)]``). + """ + f = _frequencies(frequencies) + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + require_positive(chamber_area, "chamber_area") + require_positive(pipe_area, "pipe_area") + t = duct_matrix(f, length, chamber_area, speed_of_sound=c, density=rho) + return _result( + f, t, inlet_area=pipe_area, outlet_area=pipe_area, c=c, rho=rho, + source_impedance=source_impedance, + radiation_impedance=radiation_impedance, kind="expansion chamber", + ) + + +def helmholtz_resonator( + frequencies: ArrayLike, + duct_area: float, + neck_area: float, + neck_length: float, + cavity_volume: float, + *, + resistance: float = 0.0, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult: + """Side-branch Helmholtz resonator on a duct (Bies Eqs. (8.144), (8.152)). + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param duct_area: Main-duct cross-sectional area ``S_d``, m2. + :param neck_area: Resonator neck area ``S_neck``, m2. + :param neck_length: Effective neck length ``l_e``, m. + :param cavity_volume: Cavity volume ``V``, m3. + :param resistance: Neck acoustic resistance ``R``, Pa s/m3 (default 0). + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :param source_impedance: Optional source impedance ``Z_s``, Pa s/m3. + :param radiation_impedance: Optional radiation impedance ``Z_r``, Pa s/m3. + :return: A :class:`ReactiveSilencerResult`; ``resonances`` holds + ``f_0 = (c / 2 pi) sqrt(S_neck / (l_e V))``. + """ + f = _frequencies(frequencies) + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + s_d = require_positive(duct_area, "duct_area") + zb = helmholtz_impedance( + f, neck_area, neck_length, cavity_volume, resistance=resistance, + speed_of_sound=c, density=rho, + ) + t = shunt_matrix(zb) + f0 = c / (2.0 * np.pi) * np.sqrt(neck_area / (neck_length * cavity_volume)) + return _result( + f, t, inlet_area=s_d, outlet_area=s_d, c=c, rho=rho, + source_impedance=source_impedance, + radiation_impedance=radiation_impedance, + kind="Helmholtz resonator", + resonances=np.array([f0], dtype=np.float64), + ) + + +def quarter_wave_resonator( + frequencies: ArrayLike, + duct_area: float, + length: float, + branch_area: float, + *, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult: + """Closed quarter-wave side-branch tube on a duct (Bies Eqs. (8.144), (8.146)). + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param duct_area: Main-duct cross-sectional area ``S_d``, m2. + :param length: Effective branch length ``l_e``, m. + :param branch_area: Branch tube area ``S``, m2. + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :param source_impedance: Optional source impedance ``Z_s``, Pa s/m3. + :param radiation_impedance: Optional radiation impedance ``Z_r``, Pa s/m3. + :return: A :class:`ReactiveSilencerResult`; ``resonances`` holds the odd + multiples of ``f = c / (4 l_e)`` within the frequency range. + """ + f = _frequencies(frequencies) + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + s_d = require_positive(duct_area, "duct_area") + le = require_positive(length, "length") + zb = quarter_wave_impedance( + f, le, branch_area, speed_of_sound=c, density=rho + ) + t = shunt_matrix(zb) + f_fundamental = c / (4.0 * le) + odds = np.arange(1, int(2.0 * f.max() / f_fundamental) + 2, 2) + res = np.asarray(odds * f_fundamental, dtype=np.float64) + res = res[res <= f.max()] + if res.size == 0: + res = np.array([f_fundamental], dtype=np.float64) + return _result( + f, t, inlet_area=s_d, outlet_area=s_d, c=c, rho=rho, + source_impedance=source_impedance, + radiation_impedance=radiation_impedance, + kind="quarter-wave resonator", + resonances=res, + ) + + +def extended_tube_chamber( + frequencies: ArrayLike, + length: float, + chamber_area: float, + pipe_area: float, + *, + inlet_extension: float = 0.0, + outlet_extension: float = 0.0, + speed_of_sound: float = _C_AIR, + density: float = _RHO_AIR, + source_impedance: ArrayLike | None = None, + radiation_impedance: ArrayLike | None = None, +) -> ReactiveSilencerResult: + """Extended-inlet/outlet expansion chamber (Bies §8.9.7). + + The inlet and outlet pipes extend a distance into the chamber, forming + annular quarter-wave side branches (of area ``S_exp - S_duct`` and lengths + equal to the extensions, Bies Eq. (8.156)) at the two junctions. Tuning the + extensions (classically ``L/4`` and ``L/2``) places quarter-wave peaks that + fill the ``kL = n pi`` troughs of the plain expansion chamber. With both + extensions ``0`` the result reduces exactly to :func:`expansion_chamber`. + + :param frequencies: Frequencies ``f``, Hz (1-D array). + :param length: Chamber length ``L``, m. + :param chamber_area: Chamber cross-sectional area ``S_exp``, m2. + :param pipe_area: Inlet/outlet pipe area ``S_duct``, m2. + :param inlet_extension: Inlet pipe extension into the chamber ``L_a``, m. + :param outlet_extension: Outlet pipe extension into the chamber ``L_b``, m. + :param speed_of_sound: Speed of sound ``c``, m/s. + :param density: Air density ``rho``, kg/m3. + :param source_impedance: Optional source impedance ``Z_s``, Pa s/m3. + :param radiation_impedance: Optional radiation impedance ``Z_r``, Pa s/m3. + :return: A :class:`ReactiveSilencerResult`. + """ + f = _frequencies(frequencies) + c = require_positive(speed_of_sound, "speed_of_sound") + rho = require_positive(density, "density") + s_exp = require_positive(chamber_area, "chamber_area") + s_duct = require_positive(pipe_area, "pipe_area") + la = require_non_negative(inlet_extension, "inlet_extension") + lb = require_non_negative(outlet_extension, "outlet_extension") + if s_exp <= s_duct: + raise ValueError("'chamber_area' must exceed 'pipe_area'.") + annulus = s_exp - s_duct + + elements = [] + resonances = [] + if la > 0.0: + elements.append(shunt_matrix( + quarter_wave_impedance(f, la, annulus, speed_of_sound=c, density=rho) + )) + resonances.append(c / (4.0 * la)) + elements.append(duct_matrix(f, length, s_exp, speed_of_sound=c, density=rho)) + if lb > 0.0: + elements.append(shunt_matrix( + quarter_wave_impedance(f, lb, annulus, speed_of_sound=c, density=rho) + )) + resonances.append(c / (4.0 * lb)) + t = cascade(*elements) + res = np.array(resonances, dtype=np.float64) if resonances else None + return _result( + f, t, inlet_area=s_duct, outlet_area=s_duct, c=c, rho=rho, + source_impedance=source_impedance, + radiation_impedance=radiation_impedance, + kind="extended-tube chamber", resonances=res, + ) diff --git a/tests/noise_control/test_enclosures.py b/tests/noise_control/test_enclosures.py new file mode 100644 index 00000000..439abaa3 --- /dev/null +++ b/tests/noise_control/test_enclosures.py @@ -0,0 +1,94 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +"""Tests for machine-enclosure insertion loss (Bies §7.4.2). + +Oracles: the closed form ``IL = R - C`` with +``C = 10 log10(0.3 + S_E / R_i)`` and ``R_i`` the interior room constant +(Eqs. (7.103), (7.111)); the fully absorbing limit ``C -> 10 log10 0.3`` +(``IL = R + 5.23``); consistency of ``R_i`` with +:func:`phonometry.room.room_constant`; and the panel ``R`` being supplied by +the caller (array or callable), never predicted here. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from phonometry import enclosure_insertion_loss, room_constant +from phonometry.noise_control.enclosures import EnclosureResult + + +def test_closed_form() -> None: + r = np.array([20.0, 30.0, 40.0]) + s_e, s_i, alpha = 6.0, 5.0, 0.3 + res = enclosure_insertion_loss(r, s_e, s_i, alpha) + r_i = float(room_constant(s_i, alpha)) + c = 10.0 * math.log10(0.3 + s_e / r_i) + assert np.allclose(res.correction, c) + assert np.allclose(res.insertion_loss, r - c) + + +def test_room_constant_consistency() -> None: + res = enclosure_insertion_loss([30.0], 6.0, 5.0, np.array([0.3])) + assert res.room_constant[0] == pytest.approx(float(room_constant(5.0, 0.3))) + + +def test_fully_absorbing_interior_floor() -> None: + # alpha -> 1: R_i -> inf, C -> 10 log10(0.3) = -5.23 dB, IL = R + 5.23. + res = enclosure_insertion_loss([40.0], 6.0, 5.0, 0.999999) + assert res.correction[0] == pytest.approx(10.0 * math.log10(0.3), abs=1e-3) + assert res.insertion_loss[0] == pytest.approx(40.0 - 10.0 * math.log10(0.3), abs=1e-3) + + +def test_harder_interior_wastes_more_panel() -> None: + live = enclosure_insertion_loss([40.0], 6.0, 5.0, 0.05) + dead = enclosure_insertion_loss([40.0], 6.0, 5.0, 0.5) + assert dead.insertion_loss[0] > live.insertion_loss[0] + + +def test_per_band_absorption() -> None: + r = np.array([20.0, 25.0, 30.0, 35.0]) + alpha = np.array([0.1, 0.2, 0.3, 0.4]) + res = enclosure_insertion_loss(r, 6.0, 5.0, alpha) + assert res.insertion_loss.shape == (4,) + r_i = np.asarray(room_constant(5.0, alpha)) + assert np.allclose(res.correction, 10.0 * np.log10(0.3 + 6.0 / r_i)) + + +def test_callable_panel_r() -> None: + freqs = np.array([125.0, 250.0, 500.0, 1000.0]) + + def mass_law(f: np.ndarray) -> np.ndarray: + return 20.0 * np.log10(f) - 10.0 # arbitrary user-supplied R(f) + + res = enclosure_insertion_loss(mass_law, 6.0, 5.0, 0.3, frequencies=freqs) + assert np.allclose(res.panel_transmission_loss, mass_law(freqs)) + assert res.frequencies is not None + + +def test_callable_requires_frequencies() -> None: + with pytest.raises(ValueError, match="frequencies"): + enclosure_insertion_loss(lambda f: f, 6.0, 5.0, 0.3) + + +def test_result_type_and_plot() -> None: + res = enclosure_insertion_loss( + [20.0, 30.0, 40.0], 6.0, 5.0, 0.3, + frequencies=[125.0, 250.0, 500.0], + ) + assert isinstance(res, EnclosureResult) + import matplotlib + + matplotlib.use("Agg") + assert res.plot().get_ylabel() + # No-frequency plot uses band indices. + assert enclosure_insertion_loss([20.0, 30.0], 6.0, 5.0, 0.3).plot() is not None + + +def test_validation() -> None: + with pytest.raises(ValueError): + enclosure_insertion_loss([20.0], -1.0, 5.0, 0.3) + with pytest.raises(ValueError): + enclosure_insertion_loss([20.0], 6.0, 5.0, 1.0) diff --git a/tests/noise_control/test_fdtd_crosscheck.py b/tests/noise_control/test_fdtd_crosscheck.py new file mode 100644 index 00000000..35bcf85e --- /dev/null +++ b/tests/noise_control/test_fdtd_crosscheck.py @@ -0,0 +1,93 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +"""Independent FDTD cross-check of a single 2D expansion chamber. + +The four-pole expansion-chamber transmission loss (Bies Eq. (8.111)) predicts, +for an area ratio ``m = S_exp / S_duct``, a ``TL`` peak of +``10 log10[1 + (1/4)(m - 1/m)^2]`` at ``kL = pi/2`` (``f = c / 4L``) and a +``TL = 0`` trough (fully transparent) at ``kL = pi`` (``f = c / 2L``). + +This test drives an *independent* solver -- the 2D FDTD engine of +:mod:`phonometry.simulation` -- with a plane-wave duct that widens into a +chamber and narrows back, and compares the transmitted amplitude at the peak +and trough frequencies. The measured ``20 log10`` ratio of the transmitted +amplitudes (trough over peak) must reproduce the closed-form peak ``TL``. The +run is deterministic and small (a 16x200 grid, ~3000 steps, well under a +second). +""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from phonometry.noise_control import silencers as sl +from phonometry.simulation.fdtd import CWSource, fdtd_simulation + +_C = 343.0 +_DX = 0.02 +_NY, _NX = 16, 200 +_DUCT_ROWS = (6, 7, 8, 9) # 4-cell duct; chamber is the full 16 rows -> m = 4 +_X1, _X2 = 80, 100 # chamber columns +_L = (_X2 - _X1) * _DX # chamber length, 0.4 m +_M = _NY / len(_DUCT_ROWS) # area (height) ratio = 4 + + +def _chamber_mask() -> np.ndarray: + mask = np.ones((_NY, _NX), dtype=bool) + for r in _DUCT_ROWS: + mask[r, :] = False # open the narrow inlet/outlet duct + mask[:, _X1:_X2] = False # open the full-height chamber + return mask + + +def _transmitted_rms(frequency: float) -> float: + mask = _chamber_mask() + boundaries = { + "left": "absorbing", "right": "absorbing", + "top": "rigid", "bottom": "rigid", + } + src = CWSource(ix=25, iy=7, frequency=frequency, amplitude=1.0, ramp_cycles=6.0) + steps = 3000 + dt = 0.6 * _DX / (_C * math.sqrt(2.0)) + res = fdtd_simulation( + _C, _DX, steps * dt, sources=[src], shape=(_NY, _NX), + probes=[(170, 7)], boundaries=boundaries, + absorbing_layer_cells=10, obstacle_mask=mask, + ) + p = res.pressures[0] + tail = p[int(p.size * 0.6):] # steady-state window + return float(np.sqrt(np.mean(tail**2))) + + +def test_fdtd_matches_expansion_chamber_peak_tl() -> None: + f_peak = _C / (4.0 * _L) # kL = pi/2, four-pole TL maximum + f_trough = _C / (2.0 * _L) # kL = pi, four-pole TL = 0 (transparent) + + down_peak = _transmitted_rms(f_peak) + down_trough = _transmitted_rms(f_trough) + + # The chamber transmits more at the trough than at the peak. + assert down_trough > down_peak + + # The measured attenuation difference reproduces the closed-form peak TL. + measured_delta = 20.0 * math.log10(down_trough / down_peak) + peak_tl = 10.0 * math.log10(1.0 + 0.25 * (_M - 1.0 / _M) ** 2) + assert peak_tl == pytest.approx(6.55, abs=0.05) + assert measured_delta == pytest.approx(peak_tl, abs=1.2) + + +def test_fdtd_run_is_deterministic() -> None: + f = _C / (4.0 * _L) + assert _transmitted_rms(f) == _transmitted_rms(f) + + +def test_four_pole_peak_frequency_and_trough() -> None: + # The four-pole model itself: TL peak at f = c/4L, TL = 0 at f = c/2L. + f = np.array([_C / (4.0 * _L), _C / (2.0 * _L)]) + res = sl.expansion_chamber(f, _L, _M * 0.01, 0.01) + assert res.transmission_loss[0] == pytest.approx( + 10.0 * math.log10(1.0 + 0.25 * (_M - 1.0 / _M) ** 2), abs=1e-6 + ) + assert res.transmission_loss[1] == pytest.approx(0.0, abs=1e-9) diff --git a/tests/noise_control/test_hvac.py b/tests/noise_control/test_hvac.py new file mode 100644 index 00000000..878f3ded --- /dev/null +++ b/tests/noise_control/test_hvac.py @@ -0,0 +1,126 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +"""Tests for the HVAC duct methods (Bies Chapter 8). + +Oracles: exact points of the ASHRAE end-reflection table (Bies Table 8.14) and +elbow table (Bies Table 8.11), Wells' plenum closed form (Eq. (8.275)) with its +room-constant reverberant term and limits, and the VDI 2081 flow-noise formulas +(Eqs. (8.251), (8.254)) evaluated directly. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from phonometry.noise_control import hvac + + +def test_end_reflection_table_nodes_flush() -> None: + # Bies Table 8.14, flush: D = 200 mm -> [15, 10, 5, 2, 1, 0] dB. + bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0] + res = hvac.end_reflection_loss(bands, 0.200, termination="flush") + assert np.allclose(res.values, [15, 10, 5, 2, 1, 0], atol=1e-6) + + +def test_end_reflection_table_nodes_free() -> None: + # Bies Table 8.14, free space: D = 150 mm -> [20, 14, 9, 5, 2, 1] dB. + bands = [63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0] + res = hvac.end_reflection_loss(bands, 0.150, termination="free") + assert np.allclose(res.values, [20, 14, 9, 5, 2, 1], atol=1e-6) + + +def test_end_reflection_free_exceeds_flush() -> None: + res_flush = hvac.end_reflection_loss([125.0], 0.3, termination="flush") + res_free = hvac.end_reflection_loss([125.0], 0.3, termination="free") + assert res_free.values[0] > res_flush.values[0] + + +def test_elbow_table_bands() -> None: + # Square, no vanes, unlined W = 0.3 m; W/lambda selects the table step. + c = 343.0 + # f such that W/lambda = 0.4 (in 0.28-0.55) -> 5 dB + f = 0.4 * c / 0.3 + res = hvac.elbow_insertion_loss([f], 0.3, bend_type="square") + assert res.values[0] == pytest.approx(5.0) + # W/lambda < 0.14 -> 0 dB + res0 = hvac.elbow_insertion_loss([0.1 * c / 0.3], 0.3) + assert res0.values[0] == 0.0 + + +def test_elbow_lined_beats_unlined_high_freq() -> None: + c = 343.0 + f = 1.5 * c / 0.3 # W/lambda = 1.5 (1.11-2.22 band) + unlined = hvac.elbow_insertion_loss([f], 0.3, bend_type="square").values[0] + lined = hvac.elbow_insertion_loss([f], 0.3, bend_type="square", lined=True).values[0] + assert lined > unlined # 10 vs 4 dB + + +def test_elbow_round_rejects_options() -> None: + with pytest.raises(ValueError): + hvac.elbow_insertion_loss([500.0], 0.3, bend_type="round", lined=True) + + +def test_plenum_wells_closed_form() -> None: + # TL = -10 log10[S_out(cos theta/(pi r^2) + (1-alpha)/(S_w alpha))]. + s_out, r, s_w, alpha = 0.1, 1.0, 20.0, 0.2 + tl = hvac.plenum_attenuation(s_out, r, s_w, alpha) + direct = 1.0 / (math.pi * r**2) + reverb = (1.0 - alpha) / (s_w * alpha) + expected = -10.0 * math.log10(s_out * (direct + reverb)) + assert tl == pytest.approx(expected) + + +def test_plenum_more_absorption_more_loss() -> None: + a = hvac.plenum_attenuation(0.1, 1.0, 20.0, 0.1) + b = hvac.plenum_attenuation(0.1, 1.0, 20.0, 0.5) + assert b > a + + +def test_plenum_per_band() -> None: + tl = hvac.plenum_attenuation(0.1, 1.0, 20.0, np.array([0.1, 0.3, 0.5])) + assert isinstance(tl, np.ndarray) and tl.shape == (3,) + + +def test_flow_noise_straight_duct_formula() -> None: + f = np.array([250.0]) + u, s = 10.0, 0.04 + res = hvac.flow_noise_straight_duct(f, u, s) + expected = ( + 7.0 + 50.0 * math.log10(u) + 10.0 * math.log10(s) + - 26.0 * math.log10(1.14 + 0.02 * 250.0 / u) + ) + assert res.values[0] == pytest.approx(expected) + + +def test_flow_noise_scales_with_velocity() -> None: + # Regenerated noise rises steeply with flow speed (dominant 50 log10 U term, + # tempered by the frequency term which also depends on U). + f = np.array([250.0]) + levels = [hvac.flow_noise_straight_duct(f, u, 0.04).values[0] + for u in (5.0, 10.0, 15.0)] + assert levels[0] < levels[1] < levels[2] + assert levels[2] - levels[0] > 20.0 + + +def test_flow_noise_bend_formula() -> None: + f = np.array([500.0]) + u, s, h, rho = 12.0, 0.04, 0.2, 1.206 + res = hvac.flow_noise_bend(f, u, s, h, density=rho) + lws = 30.0 * math.log10(u) + 10.0 * math.log10(s) + 10.0 * math.log10(rho) + 117.0 + ns = 500.0 * h / u + expected = lws - 10.0 * math.log10(1.0 + 0.165 * ns**2) + 30.0 * math.log10(u) - 103.0 + assert res.values[0] == pytest.approx(expected) + + +def test_plot_and_validation() -> None: + import matplotlib + + matplotlib.use("Agg") + hvac.end_reflection_loss([125.0, 250.0], 0.3).plot() + hvac.flow_noise_straight_duct([250.0, 500.0], 10.0, 0.04).plot() + with pytest.raises(ValueError): + hvac.end_reflection_loss([125.0], 0.3, termination="bad") + with pytest.raises(ValueError): + hvac.plenum_attenuation(0.1, 1.0, 20.0, 1.0) diff --git a/tests/noise_control/test_silencers.py b/tests/noise_control/test_silencers.py new file mode 100644 index 00000000..46fbacb8 --- /dev/null +++ b/tests/noise_control/test_silencers.py @@ -0,0 +1,148 @@ +# Copyright (c) 2026. Jose M. Requena-Plens +"""Tests for reactive silencers (Bies §8.8-8.9, four-pole method). + +Oracles: the simple-expansion-chamber closed form ``TL = 10 log10[1 + (1/4) +(m - 1/m)^2 sin^2(kL)]`` (Bies Eq. (8.111)) with its published peak values +(m = 2 -> 1.94, 4 -> 6.55, 8 -> 12.18, 16 -> 18.10 dB), the four-pole TL +reducing to the side-branch closed form (Bies Eq. (8.73)), the quarter-wave +tube tuning ``f = c/(4 l_e)`` (Example 8.1: 56.6 Hz), the Helmholtz resonance +``f_0 = (c/2 pi) sqrt(S/(l_e V))``, matrix reciprocity, and limiting cases. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from phonometry.noise_control import silencers as sl + + +def _k(f: np.ndarray, c: float = 343.0) -> np.ndarray: + return 2.0 * np.pi * f / c + + +def test_expansion_chamber_matches_closed_form() -> None: + f = np.linspace(20.0, 2000.0, 2000) + for m in (2.0, 4.0, 8.0, 16.0): + s_duct, length = 0.01, 0.3 + res = sl.expansion_chamber(f, length, m * s_duct, s_duct) + cf = 10.0 * np.log10(1.0 + 0.25 * (m - 1.0 / m) ** 2 * np.sin(_k(f) * length) ** 2) + assert np.allclose(res.transmission_loss, cf, atol=1e-9) + + +@pytest.mark.parametrize( + ("m", "peak"), [(2.0, 1.94), (4.0, 6.55), (8.0, 12.18), (16.0, 18.10)] +) +def test_expansion_chamber_peaks(m: float, peak: float) -> None: + f = np.linspace(20.0, 2000.0, 6000) + res = sl.expansion_chamber(f, 0.3, m * 0.01, 0.01) + assert res.transmission_loss.max() == pytest.approx(peak, abs=0.02) + + +def test_expansion_chamber_troughs_zero() -> None: + # TL = 0 at kL = n pi (chamber transparent, no dissipation). + c, length, s_duct = 343.0, 0.3, 0.01 + f_trough = c / (2.0 * length) # kL = pi + res = sl.expansion_chamber(np.array([f_trough]), length, 0.04, s_duct) + assert res.transmission_loss[0] == pytest.approx(0.0, abs=1e-9) + + +def test_side_branch_reduces_to_closed_form() -> None: + # A shunt branch Z_b between equal-area ducts: TL = 20 log10|1 + rho c/(2 Sd Zb)|. + f = np.linspace(50.0, 500.0, 200) + zb = sl.quarter_wave_impedance(f, 0.5, 0.002) + t = sl.shunt_matrix(zb) + tl = sl.transmission_loss(t, inlet_area=0.01, outlet_area=0.01) + closed = 20.0 * np.log10(np.abs(1.0 + 1.206 * 343.0 / (2.0 * 0.01 * zb))) + assert np.allclose(tl, closed, atol=1e-9) + + +def test_quarter_wave_tuning_frequency() -> None: + # Bies Example 8.1: l_e = 1.516 m, c = 343.24 -> first peak at 56.6 Hz. + # Restrict below the third harmonic (169.8 Hz) so the fundamental is the + # only resonance in range (all odd harmonics peak equally for a QWT). + f = np.linspace(10.0, 120.0, 8000) + area = np.pi * 0.05**2 / 4.0 + res = sl.quarter_wave_resonator(f, area, 1.516, area, speed_of_sound=343.24) + peak_f = f[np.argmax(res.transmission_loss)] + assert peak_f == pytest.approx(56.6, abs=0.2) + assert res.resonances[0] == pytest.approx(56.6, abs=0.1) + + +def test_helmholtz_resonance_frequency() -> None: + c, s_neck, le, vol = 343.0, 1e-4, 0.02, 1e-3 + res = sl.helmholtz_resonator(np.linspace(50.0, 400.0, 50), 0.01, s_neck, le, vol) + f0 = c / (2.0 * np.pi) * np.sqrt(s_neck / (le * vol)) + assert res.resonances[0] == pytest.approx(f0) + + +def test_helmholtz_peak_at_resonance() -> None: + # Lossless resonator: TL peaks sharply at f_0. + c, s_neck, le, vol = 343.0, 1e-4, 0.02, 1e-3 + f0 = c / (2.0 * np.pi) * np.sqrt(s_neck / (le * vol)) + f = np.linspace(0.6 * f0, 1.4 * f0, 4000) + res = sl.helmholtz_resonator(f, 0.01, s_neck, le, vol) + assert f[np.argmax(res.transmission_loss)] == pytest.approx(f0, rel=0.02) + + +def test_insertion_loss_identity_is_zero() -> None: + f = np.linspace(50.0, 500.0, 50) + ident = sl.duct_matrix(f, 0.0, 0.01) # zero-length duct = identity + il = sl.insertion_loss(ident, source_impedance=400.0, radiation_impedance=50.0) + assert np.allclose(il, 0.0, atol=1e-12) + + +def test_transfer_matrix_reciprocity() -> None: + # Reciprocal passive elements have det(T) = 1. + f = np.linspace(50.0, 500.0, 20) + t = sl.expansion_chamber(f, 0.3, 0.04, 0.01).transfer_matrix + dets = np.linalg.det(t) + assert np.allclose(dets, 1.0, atol=1e-9) + + +def test_cascade_order() -> None: + f = np.linspace(50.0, 500.0, 10) + a = sl.duct_matrix(f, 0.1, 0.01) + b = sl.duct_matrix(f, 0.2, 0.01) + # Two concatenated ducts equal one duct of the summed length. + both = sl.cascade(a, b) + whole = sl.duct_matrix(f, 0.3, 0.01) + assert np.allclose(both, whole, atol=1e-9) + + +def test_extended_tube_reduces_to_expansion_chamber() -> None: + f = np.linspace(20.0, 2000.0, 500) + ext = sl.extended_tube_chamber(f, 0.3, 0.04, 0.01) + plain = sl.expansion_chamber(f, 0.3, 0.04, 0.01) + assert np.allclose(ext.transmission_loss, plain.transmission_loss, atol=1e-9) + + +def test_extended_tube_fills_trough() -> None: + # An inlet extension tuned to L/4 raises the TL at the first chamber trough. + c, length = 343.0, 0.4 + f = np.array([c / (2.0 * length)]) # first plain-chamber trough (TL = 0) + plain = sl.expansion_chamber(f, length, 0.04, 0.01) + tuned = sl.extended_tube_chamber( + f, length, 0.04, 0.01, inlet_extension=length / 4.0 + ) + assert plain.transmission_loss[0] == pytest.approx(0.0, abs=1e-9) + assert tuned.transmission_loss[0] > 3.0 + + +def test_insertion_loss_present_when_impedances_given() -> None: + f = np.linspace(50.0, 500.0, 20) + res = sl.expansion_chamber( + f, 0.3, 0.04, 0.01, source_impedance=4e4, radiation_impedance=5e3 + ) + assert res.insertion_loss is not None + plain = sl.expansion_chamber(f, 0.3, 0.04, 0.01) + assert plain.insertion_loss is None + + +def test_validation() -> None: + with pytest.raises(ValueError): + sl.expansion_chamber([0.0], 0.3, 0.04, 0.01) + with pytest.raises(ValueError): + sl.expansion_chamber([100.0], 0.3, -0.04, 0.01) + with pytest.raises(ValueError, match="must exceed"): + sl.extended_tube_chamber([100.0], 0.3, 0.01, 0.02) diff --git a/tests/test_package_architecture.py b/tests/test_package_architecture.py index 8a9439e9..299b0a93 100644 --- a/tests/test_package_architecture.py +++ b/tests/test_package_architecture.py @@ -35,6 +35,9 @@ ("building", "vibration"), # double-wall cavity fill uses the porous equivalent-fluid model ("building", "materials"), + # HVAC plenum and machine enclosures reuse the room constant + # R = S*alpha/(1 - alpha) of the steady-state room field + ("noise_control", "room"), } From e8f5f222ffc5dec69a4bd3eb44ba6a6d107d77cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20M=2E=20Requena=20Plens?= Date: Sun, 19 Jul 2026 16:59:18 +0200 Subject: [PATCH 3/5] feat: expose piston and noise-control APIs, conformance and indexes Re-export the piston and high-level noise-control functions at the top level, add the noise_control API-reference section and the piston module to the taxonomy, and register 13 additive numerical conformance checks (piston resistance/reactance/mass/directivity and the silencer, plenum, end-reflection and enclosure closed forms). Regenerate the conformance report, the API reference pages/sidebar and the quick-reference table, and note the additions in the changelog. --- CHANGELOG.md | 33 ++++ docs/CONFORMANCE.md | 34 ++-- docs/README.md | 1 + docs/api-reference.md | 14 ++ scripts/api_taxonomy.py | 12 ++ scripts/conformance_report.py | 171 +++++++++++++++++++ site/src/content/docs/reference/api/index.md | 9 + site/src/generated/api-sidebar.mjs | 10 ++ src/phonometry/__init__.py | 45 +++++ 9 files changed, 312 insertions(+), 17 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 434342c6..820dd566 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -67,6 +67,39 @@ and this project adheres to [Semantic Versioning](https://semver.org/). `W = 0.5 |F|^2 Re{Y}` (Cremer, Heckl & Petersson 2005, *Structure-Borne Sound* 3e, Table 5.1), the theoretical companions of the measured ISO 7626 mobilities. +- Industrial noise control: silencers, HVAC and machine enclosures, a new + `phonometry.noise_control` domain (Bies, Hansen & Howard, *Engineering Noise + Control* 5th ed.; Munjal, *Acoustics of Ducts and Mufflers*). `expansion_chamber`, + `helmholtz_resonator`, `quarter_wave_resonator` and `extended_tube_chamber` + build **reactive silencers** by the one-dimensional four-pole + (transmission-matrix) method, returning a `ReactiveSilencerResult` with the + transmission loss, the insertion loss for configurable source/radiation + impedances, the compound transfer matrix and `.plot()`; the expansion-chamber + transmission loss matches the closed form `10 lg[1 + (1/4)(m - 1/m)^2 sin^2 kL]` + (Bies Eq. (8.111)) exactly, and the four-pole machinery is available under + `phonometry.noise_control` (`duct_matrix`, `shunt_matrix`, `cascade`, + `transmission_loss`, `insertion_loss`, `helmholtz_impedance`, + `quarter_wave_impedance`). The + `noise_control.hvac` module adds the duct end reflection (`end_reflection_loss`, + ASHRAE Table 8.14), bend insertion loss (`elbow_insertion_loss`, Table 8.11), + the plenum transmission loss by Wells' method (`plenum_attenuation`) and the + flow-generated (self) noise of straight ducts and mitred bends + (`flow_noise_straight_duct`, `flow_noise_bend`, VDI 2081). `enclosure_insertion_loss` + gives the machine-enclosure insertion loss `IL = R - C` (Bies Eqs. (7.103), + (7.111)) from a **caller-supplied** panel transmission loss `R` (array or + callable; never predicted here) and the interior room constant reused from + `phonometry.room.room_constant`. The four-pole expansion chamber is + cross-checked against the independent 2D FDTD solver. +- Radiation of a rigid circular piston in an infinite baffle + (`phonometry.electroacoustics.piston`, re-exported at the top level). `radiating_piston` + computes the radiation impedance `rho c S (R1 + j X1)` from the piston + resistance `R1 = 1 - 2 J1(2ka)/2ka` and reactance `X1 = 2 H1(2ka)/2ka` + (Beranek & Mellow, *Acoustics* 2nd ed., Eqs. (13.117), (13.118)), the + low-frequency radiation mass `8 rho a^3 / 3` (Eq. (4.151)), the far-field + directivity `2 J1(ka sin theta)/(ka sin theta)` and the directivity index, + returning a `RadiatingPistonResult` with `.plot()`; the building blocks + `piston_resistance`, `piston_reactance` and `piston_directivity` are exposed + directly. - Objective speech intelligibility STOI and ESTOI (`phonometry.hearing.objective_intelligibility.stoi`, also `phonometry.stoi`). `stoi(clean, degraded, fs)` computes the short-time objective intelligibility diff --git a/docs/CONFORMANCE.md b/docs/CONFORMANCE.md index 836c2784..63cb2bdb 100644 --- a/docs/CONFORMANCE.md +++ b/docs/CONFORMANCE.md @@ -15,7 +15,7 @@ ## Numerical conformance report -✅ **338/338 conformance checks pass** across 42 domains and 210 standards - filters class 1 - weightings within IEC 61672-1 class 1. +✅ **338/338 conformance checks pass** across 42 domains and 211 standards - filters class 1 - weightings within IEC 61672-1 class 1. Each row pins a standard clause to its expected normative value and the value the library computes. Every section below is collapsible and stays collapsed while all of its rows pass; a section with any failing row opens automatically. @@ -690,32 +690,32 @@ Only **Butterworth** (the library default) and **Chebyshev-II** are class-compli
-✅ Panel & aperture sound insulation (Bies / Hopkins / Cremer): 100% (11/11) +✅ Electroacoustics: 100% (6/6) | Standard | Quantity | Expected (norm) | Computed | Δ | Status | |:---|:---|:---|:---|:---|:---:| -| Bies 5e Eq. 7.40 (mass law) | 6 dB per octave (500 -> 1000 Hz) | 6.0206 dB (+/-0.01 dB) | 6.02 dB | -0.001 dB | ✅ | -| Bies 5e Eq. 7.40 (mass law) | 6 dB per doubling of mass | 6.0206 dB (+/-0.01 dB) | 6.02 dB | -0.001 dB | ✅ | -| Bies 5e Eq. 7.42 (field incidence) | One-third-octave correction 5.5 dB | 5.5 dB (+/-0.001 dB) | 5.5 dB | 0 dB | ✅ | -| Hopkins Eq. 2.201 / Bies Eq. 7.3 | Coincidence frequency, 6 mm glass | 2079 Hz (+/-3%) | 2107.3639 Hz | 28.364 Hz | ✅ | -| Cremer Table 5.1 | Thin-plate point impedance Z = 8 sqrt(B' m'') | 2529.8221 N.s/m (+/-0 N.s/m) | 2529.8221 N.s/m | 0 N.s/m | ✅ | -| Cremer Table 5.1 | Infinite-beam mobility phase -45 deg | -45 deg (+/-0 deg) | -45 deg | 0 deg | ✅ | -| Hopkins Eq. 2.229 (Leppington/Maidanik) | Radiation efficiency at f = 2 fc | 1.4142 (+/-0) | 1.4142 | 0 | ✅ | -| Bies Eq. 7.62 / Hopkins Eq. 4.73 | Mass-air-mass resonance f0, empty cavity | 76.9484 Hz (+/-0.5%) | 76.8521 Hz | -0.096 Hz | ✅ | -| Bies Eq. 7.64 (double wall) | Below f0 = mass law of the combined mass | 11.6144 dB (+/-0 dB) | 11.6144 dB | 0 dB | ✅ | -| Hopkins Eq. 4.92 (composite) | 1 % open area caps R at 10 lg(S/Sa) | 20 dB (+/-0.05 dB) | 19.9996 dB | 0 dB | ✅ | -| Hopkins Eq. 4.99/4.101 (Gomperts slit) | Transmission maximum at first resonance | 1544.9615 Hz (+/-15 Hz) | 1542.9615 Hz | -2 Hz | ✅ | +| Beranek & Mellow 2e Eq. (13.117) | Piston resistance R1(x) = 1 - 2 J1(x)/x at x = 2ka = 2 | 0.423275 (+/-0.00001) | 0.423275 | 0 | ✅ | +| Beranek & Mellow 2e Eq. (13.118) | Piston reactance X1(x) = 2 H1(x)/x at x = 2ka = 2 | 0.646764 (+/-0.00001) | 0.646764 | 0 | ✅ | +| Beranek & Mellow 2e Eq. (13.117) (low-frequency limit) | R1 -> (ka)^2/2 as ka -> 0 (x = 0.02, ka = 0.01) | 0.00005 (+/-0.01%) | 0.00005 | 0 | ✅ | +| Beranek & Mellow 2e Eq. (4.151) | Radiation mass M = 8 rho a^3 / 3 (a = 0.1 m, rho = 1.206) | 0.003216 kg (+/-0 kg) | 0.003216 kg | 0 kg | ✅ | +| Beranek & Mellow 2e Eq. (13.102), Table 14.1 | First directivity null at ka sin(theta) = 3.8317 (first zero of J1) | 0 (+/-0.000001) | 0 | 0 | ✅ | +| Beranek & Mellow 2e §4.19 (half-space baffle) | Directivity index DI -> 10 lg 2 = 3.01 dB as ka -> 0 | 3.0103 dB (+/-0.001 dB) | 3.0103 dB | 0 dB | ✅ |
-✅ Atmospheric refraction (Salomons rays / GFPE): 100% (3/3) +✅ Industrial noise control: 100% (8/8) | Standard | Quantity | Expected (norm) | Computed | Δ | Status | |:---|:---|:---|:---|:---|:---:| -| Salomons Sec. 4.4 (ray turning height, linear profile) | Turning height of a 10 deg ray vs Rc(1 - cos theta0) (circular arc), m | 26.457 m (+/-0.1 m) | 26.457 m | 0 m | ✅ | -| Salomons Eq. (3.4) (GFPE vs spherical-wave ground effect, homogeneous) | PE relative level at 500 m over grassland vs Weyl-Van der Pol, dB | -13.919 dB (+/-0.5 dB) | -13.842 dB | 0.077 dB | ✅ | -| Salomons Eq. (3.4) (GFPE hard ground vs two-ray, homogeneous) | PE relative level at 500 m over a rigid ground vs the coherent two-ray, dB | 5.997 dB (+/-0.6 dB) | 5.593 dB | -0.405 dB | ✅ | +| Bies 5e Eq. (8.111) | Expansion-chamber peak TL = 10 lg[1 + (1/4)(m - 1/m)^2], m = 4 at kL = pi/2 | 6.5472 dB (+/-0 dB) | 6.5472 dB | 0 dB | ✅ | +| Bies 5e Eq. (8.111) | Expansion-chamber trough TL = 0 at kL = pi (chamber transparent) | 0 dB (+/-0 dB) | 0 dB | 0 dB | ✅ | +| Bies 5e Eq. (8.44) / Example 8.1 | Quarter-wave tube tuning f = c/(4 l_e), l_e = 1.516 m -> 56.6 Hz | 56.6 Hz (+/-0.1 Hz) | 56.6 Hz | 0.003 Hz | ✅ | +| Bies 5e Eq. (8.46) | Helmholtz resonance f0 = (c/2pi) sqrt(S/(l_e V)) (S=1e-4, l_e=0.02, V=1e-3) | 122.067 Hz (+/-0 Hz) | 122.067 Hz | 0 Hz | ✅ | +| Bies 5e Eq. (8.73) | Side-branch TL = 20 lg|1 + rho c/(2 Sd Zb)| (QWT branch, closed form) | 0.1638 dB (+/-0 dB) | 0.1638 dB | 0 dB | ✅ | +| Bies 5e Eq. (8.275) (Wells' plenum method) | Plenum TL = -10 lg[S_out(cos0/pi r^2 + (1-a)/(Sw a))] (S_out=.1,r=1,Sw=20,a=.2) | 12.8541 dB (+/-0 dB) | 12.8541 dB | 0 dB | ✅ | +| Bies 5e Table 8.14 (ASHRAE end reflection, flush) | Duct end reflection D = 200 mm at 125 Hz = 10 dB (table node) | 10 dB (+/-0 dB) | 10 dB | 0 dB | ✅ | +| Bies 5e Eqs. (7.103), (7.111) (enclosure, fully absorbing limit) | Enclosure correction C -> 10 lg 0.3 = -5.23 dB as alpha_i -> 1 | -5.2288 dB (+/-0.001 dB) | -5.2288 dB | 0 dB | ✅ |
diff --git a/docs/README.md b/docs/README.md index 6b1f7ab0..f75b38e1 100644 --- a/docs/README.md +++ b/docs/README.md @@ -23,6 +23,7 @@ Full documentation for phonometry. Also available as a website: - [Objective intelligibility (STOI & ESTOI)](objective-intelligibility.md): the correlation-based intelligibility measures for time-frequency weighted noisy speech from a clean/degraded pair: STOI (Taal et al. 2011), the clipped per-band envelope correlation, and ESTOI (Jensen & Taal 2016), the row- and column-normalised spectral correlation that tracks modulated maskers - [Electroacoustics: distortion & frequency response](electroacoustics.md): the IEC 60268-3 distortion set (THD, nth-order harmonic, THD+N and SINAD via AES17, SMPTE and CCIF intermodulation, dynamic intermodulation and weighted THD), the AES17 dynamic range and idle channel noise, and the Bendat & Piersol H1/H2 frequency-response estimators with the ordinary coherence γ² - [Swept-sine distortion and phase utilities](swept-sine-distortion.md): harmonic separation from one exponential sweep (Farina 2000) with the synchronized swept-sine of Novak et al. 2015 for coherent harmonic phases, THD as a function of the excitation frequency, and minimum phase from |H| (real cepstrum), group delay and excess phase +- [Industrial noise control: silencers, HVAC & enclosures](noise-control.md): reactive silencers by the four-pole transmission-matrix method (expansion chambers, Helmholtz and quarter-wave resonators, extended tubes) with transmission and insertion loss, the HVAC duct methods (end reflection, elbows, plenums and flow-generated noise) and machine-enclosure insertion loss from a supplied panel R and the interior room constant (Bies, Hansen & Howard) - [Programme loudness & true peak](program-loudness.md): the ITU-R BS.1770-5 programme loudness (K-weighting, gated 400 ms blocks, channel weights including the Annex 3 positions) and the oversampled true-peak level in dBTP, with the EBU R 128 −23 LUFS practice, the Tech 3341 EBU Mode momentary/short-term/integrated meters and the Tech 3342 loudness range, validated against the official EBU test signals - [Underwater acoustics: radiated noise & pile driving](underwater-acoustics.md): the ISO 18405 reference levels (SPL, SEL, peak re 1 µPa), the ISO 17208 ship radiated noise level and equivalent monopole source level via the Lloyd's-mirror correction, and the ISO 18406 single-strike, peak and cumulative pile-driving sound exposure - [Underwater sound propagation](underwater-propagation.md): closed-form transmission loss (geometrical spreading plus volume absorption by Francois-Garrison, Ainslie-McColm or Thorp), the speed of sound in sea water (UNESCO/Chen-Millero, Del Grosso, Mackenzie) with the sound-speed profile, the passive/active sonar equation, seabed reflection loss (Rayleigh), the ocean ambient-noise spectrum (Wenz wind/thermal plus JOMOPANS-ECHO ship-traffic source levels) and numerical solvers (normal modes, ray tracing, parabolic equation) diff --git a/docs/api-reference.md b/docs/api-reference.md index dca8129d..57bf7922 100644 --- a/docs/api-reference.md +++ b/docs/api-reference.md @@ -356,6 +356,20 @@ cycle and warn on use; they are removed in 4.0. | `synchronized_sweep_signal` | `function` | **Synchronized exponential sweep (Novak et al. 2015, Eqs. 47/49).**
• `fs`, `f1` < `f2` [Hz], `seconds` (quantized by the rounding)
• `amplitude` (Default: 1.0), `fade` (Default: 0.0)
Rate L = round(f1·T̃/ln(f2/f1))/f1 → f1·L integer, so a L·ln(n) shift equals the nth harmonic and harmonic phases become system properties | `x = synchronized_sweep_signal(48000, 20.0, 6000.0, 4.0)`

• sweep samples (duration L·ln(f2/f1)) | | `swept_sine_distortion` | `function` | **Harmonic separation and THD(f) from one sweep (Farina 2000 / Novak et al. 2015).**
• `recorded`, `fs`, plus the generator parameters `f1`, `f2`, `seconds`
• `method`: 'synchronized' (Default; analytic deconvolution, meaningful phases) / 'farina' (classical ESS, magnitudes only)
• `n_harmonics` (Default: 5), `ir_length` (Default: largest power of two fitting the closest arrivals)
• `amplitude` (Default: 1.0), `fade` (Default: generator default), `remove_dc` (Default: True) | `res = swept_sine_distortion(y, fs, 20.0, 6000.0, 4.0)`

• `SweptSineDistortionResult` | | `SweptSineDistortionResult` | `dataclass` | **Swept-sine harmonic separation.**
• `frequencies` [Hz], `harmonic_responses`: complex H₁..H_N, `harmonic_irs` (windowed, centred), `delays` L·ln(n) [s]
• `thd_frequencies` [Hz], `thd`, `distortion_ratios`: \|Hₙ(n·f)\|/\|H₁(f)\| per order
• `fs`, `f1`, `f2`, `duration`, `rate` L, `method`, `n_harmonics`
• `.plot()`: \|Hₙ\| magnitudes + THD(f) | `res.thd, res.harmonic_responses` | +| `radiating_piston` | `function` | **Radiation of a rigid baffled circular piston (Beranek & Mellow §4.19/§13.7).**
• `radius` a [m], `frequencies` [Hz]
• `speed_of_sound` (Default: 343), `density` (Default: 1.206)
• `angles`: directivity polar angles [rad] (Default: None) | `res = radiating_piston(0.1, freqs)`

• `RadiatingPistonResult` | +| `piston_resistance` / `piston_reactance` | `function` | **Normalized piston radiation impedance R1 = 1−2J1(x)/x, X1 = 2H1(x)/x, x = 2ka (Beranek & Mellow Eqs. 13.117/13.118).** | `piston_resistance(2.0) # 0.4233` | +| `piston_directivity` | `function` | **Far-field piston directivity 2·J1(ka·sinθ)/(ka·sinθ).**
• `ka`, `theta` [rad] | `piston_directivity(5.0, 0.3)` | +| `RadiatingPistonResult` | `dataclass` | **Piston radiation result.**
• `ka`, `resistance`/`reactance`, `radiation_resistance`/`radiation_reactance` [N·s/m], `radiation_mass` = 8ρa³/3 [kg], `directivity_index` [dB], `directivity`
• `.plot()`: R1/X1 vs ka | `res.radiation_mass` | +| `expansion_chamber` | `function` | **Expansion-chamber silencer TL/IL (Bies Eq. 8.111, four-pole).**
• `frequencies`, `length` L, `chamber_area`, `pipe_area`
• `source_impedance`/`radiation_impedance` (Default: None → no IL) | `res = expansion_chamber(f, 0.3, 0.04, 0.01)`

• `ReactiveSilencerResult` | +| `helmholtz_resonator` / `quarter_wave_resonator` | `function` | **Side-branch resonator silencers (Bies Eqs. 8.46/8.44).**
• Helmholtz: `duct_area`, `neck_area`, `neck_length`, `cavity_volume`
• Quarter-wave: `duct_area`, `length`, `branch_area` | `quarter_wave_resonator(f, 0.01, 1.5, 2e-3)` | +| `extended_tube_chamber` | `function` | **Extended-inlet/outlet expansion chamber (Bies §8.9.7); reduces to `expansion_chamber` at zero extension.**
• `length`, `chamber_area`, `pipe_area`, `inlet_extension`, `outlet_extension` | `extended_tube_chamber(f, 0.4, 0.04, 0.01, inlet_extension=0.1)` | +| `ReactiveSilencerResult` | `dataclass` | **Reactive-silencer result.**
• `frequencies`, `transmission_loss`, `insertion_loss` (or None), `transfer_matrix`, `kind`, `resonances`
• `.plot()`: TL/IL vs frequency | `res.transmission_loss` | +| `end_reflection_loss` / `elbow_insertion_loss` | `function` | **HVAC duct end reflection & bend insertion loss (Bies Tables 8.14/8.11, ASHRAE).**
• end: `frequencies`, `diameter`, `termination` 'flush'/'free'
• elbow: `frequencies`, `width`, `bend_type`, `vanes`, `lined` | `end_reflection_loss(bands, 0.3)`

• `HvacSpectrumResult` | +| `plenum_attenuation` | `function` | **Plenum-chamber TL by Wells' method (Bies Eq. 8.275).**
• `exit_area`, `line_of_sight`, `wall_area`, `mean_absorption`, `angle` | `plenum_attenuation(0.1, 1.0, 20.0, 0.2) # dB` | +| `flow_noise_straight_duct` / `flow_noise_bend` | `function` | **HVAC flow-generated (self) noise sound power (VDI 2081, Bies Eqs. 8.251/8.254).**
• `frequencies`, `flow_velocity`, `area` (+ `height` for the bend) | `flow_noise_straight_duct(bands, 10.0, 0.04)`

• `HvacSpectrumResult` | +| `HvacSpectrumResult` | `dataclass` | **Per-frequency HVAC attenuation or regenerated Lw.**
• `frequencies`, `values`, `quantity`, `label`
• `.plot()` | `res.values` | +| `enclosure_insertion_loss` | `function` | **Machine-enclosure insertion loss IL = R − C (Bies Eqs. 7.103/7.111); panel R supplied by the caller, never predicted.**
• `panel_transmission_loss` (per-band array or callable of f)
• `external_area`, `internal_area`, `internal_absorption`
• `frequencies` (required for a callable R) | `enclosure_insertion_loss(R, 6.0, 5.0, 0.3)`

• `EnclosureResult` | +| `EnclosureResult` | `dataclass` | **Enclosure insertion-loss result.**
• `panel_transmission_loss`, `correction` C, `insertion_loss`, `room_constant`, `external_area`, `internal_area`
• `.plot()` | `res.insertion_loss` | | `program_loudness` | `function` | **EBU R 128 measurement set of a programme (BS.1770-5 + Tech 3341/3342).**
• `x`: signal (1D mono or 2D [channels, samples], 1.0 = 0 dBFS)
• `fs` [Hz]
• `weights`: per-channel Gi (Default: Table 3 by channel count)
• `momentary_step` [s] (Default: 0.01), `short_term_step` [s] (Default: 0.1)
• `oversample`: true-peak factor (Default: None → reaches 192 kHz) | `res = program_loudness(x, fs)`

• `ProgramLoudnessResult` | | `integrated_loudness` | `function` | **Programme (integrated) loudness with the two-stage gate (BS.1770-5 Annex 1).**
• `x`, `fs` [Hz]
• `weights` (Default: Table 3 by channel count)
Absolute gate −70 LKFS, relative −10 LU | `integrated_loudness(x, fs) # LUFS` | | `loudness_range` | `function` | **Loudness range LRA (EBU Tech 3342).**
• `short_term_loudness`: S readings [LUFS] at ≥ 10 Hz
Cascaded gate (−70 LUFS abs, −20 LU rel), 10th-95th percentile spread | `loudness_range(res.short_term) # LU` | diff --git a/scripts/api_taxonomy.py b/scripts/api_taxonomy.py index 732a5abe..6b4d2d3a 100644 --- a/scripts/api_taxonomy.py +++ b/scripts/api_taxonomy.py @@ -232,6 +232,17 @@ class Section: "phonometry.electroacoustics.distortion", "phonometry.electroacoustics.frequency_response", "phonometry.electroacoustics.swept_sine", + "phonometry.electroacoustics.piston", + ), + ), + Section( + key="noise_control", + label_en="Industrial noise control", + label_es="Control de ruido industrial", + modules=( + "phonometry.noise_control.silencers", + "phonometry.noise_control.hvac", + "phonometry.noise_control.enclosures", ), ), Section( @@ -295,6 +306,7 @@ class Section: "underwater": ("underwater",), "power": ("emission",), "electroacoustics": ("electroacoustics",), + "noise_control": ("noise_control",), "broadcast": ("broadcast",), "metrology": ("metrology",), "spectra": ("metrology",), diff --git a/scripts/conformance_report.py b/scripts/conformance_report.py index d0e3393f..4818bec0 100644 --- a/scripts/conformance_report.py +++ b/scripts/conformance_report.py @@ -5743,6 +5743,177 @@ def _chk_atm_pe_hard_ground() -> Outcome: k = 2.0 * math.pi * freq / 343.0 two_ray = 20.0 * math.log10(abs(1.0 + (r1 / r2) * cmath.exp(1j * k * (r2 - r1)))) return numeric(two_ray, float(pe.level_at_height(zr)[i]), 0.6, unit="dB", places=3) +# Domain - Electroacoustics (baffled piston, Beranek & Mellow 2e) +# =========================================================================== +_ELECTROACOUSTICS = "Electroacoustics" + + +@register( + _ELECTROACOUSTICS, + "Beranek & Mellow 2e Eq. (13.117)", + "Piston resistance R1(x) = 1 - 2 J1(x)/x at x = 2ka = 2", +) +def _chk_piston_resistance() -> Outcome: + return numeric(0.423275, float(ph.piston_resistance(2.0)), 1e-5, places=6) + + +@register( + _ELECTROACOUSTICS, + "Beranek & Mellow 2e Eq. (13.118)", + "Piston reactance X1(x) = 2 H1(x)/x at x = 2ka = 2", +) +def _chk_piston_reactance() -> Outcome: + return numeric(0.646764, float(ph.piston_reactance(2.0)), 1e-5, places=6) + + +@register( + _ELECTROACOUSTICS, + "Beranek & Mellow 2e Eq. (13.117) (low-frequency limit)", + "R1 -> (ka)^2/2 as ka -> 0 (x = 0.02, ka = 0.01)", +) +def _chk_piston_resistance_limit() -> Outcome: + ka = 0.01 + return numeric(ka**2 / 2.0, float(ph.piston_resistance(2.0 * ka)), 1e-4, + rel=True, places=8) + + +@register( + _ELECTROACOUSTICS, + "Beranek & Mellow 2e Eq. (4.151)", + "Radiation mass M = 8 rho a^3 / 3 (a = 0.1 m, rho = 1.206)", +) +def _chk_piston_radiation_mass() -> Outcome: + res = ph.radiating_piston(0.1, [100.0], density=1.206) + return numeric(8.0 * 1.206 * 0.1**3 / 3.0, res.radiation_mass, 1e-9, + unit="kg", places=8) + + +@register( + _ELECTROACOUSTICS, + "Beranek & Mellow 2e Eq. (13.102), Table 14.1", + "First directivity null at ka sin(theta) = 3.8317 (first zero of J1)", +) +def _chk_piston_directivity_null() -> Outcome: + # ka sin(theta) = 3.8317 at ka = 3.8317, theta = pi/2. + d = float(ph.piston_directivity(3.8317059702075125, math.pi / 2.0)) + return numeric(0.0, d, 1e-6, places=8) + + +@register( + _ELECTROACOUSTICS, + "Beranek & Mellow 2e §4.19 (half-space baffle)", + "Directivity index DI -> 10 lg 2 = 3.01 dB as ka -> 0", +) +def _chk_piston_directivity_index() -> Outcome: + res = ph.radiating_piston(0.01, [1.0]) + return numeric(10.0 * math.log10(2.0), float(res.directivity_index[0]), + 1e-3, unit="dB") + + +# =========================================================================== +# Domain - Industrial noise control (Bies 5e; silencers, HVAC, enclosures) +# =========================================================================== +_NOISE_CONTROL = "Industrial noise control" + + +@register( + _NOISE_CONTROL, + "Bies 5e Eq. (8.111)", + "Expansion-chamber peak TL = 10 lg[1 + (1/4)(m - 1/m)^2], m = 4 at kL = pi/2", +) +def _chk_expansion_chamber_peak() -> Outcome: + c, length, s_duct = 343.0, 0.3, 0.01 + f = np.array([c / (4.0 * length)]) # kL = pi/2 + res = ph.expansion_chamber(f, length, 4.0 * s_duct, s_duct, speed_of_sound=c) + expected = 10.0 * math.log10(1.0 + 0.25 * (4.0 - 0.25) ** 2) + return numeric(expected, float(res.transmission_loss[0]), 1e-6, unit="dB") + + +@register( + _NOISE_CONTROL, + "Bies 5e Eq. (8.111)", + "Expansion-chamber trough TL = 0 at kL = pi (chamber transparent)", +) +def _chk_expansion_chamber_trough() -> Outcome: + c, length, s_duct = 343.0, 0.3, 0.01 + f = np.array([c / (2.0 * length)]) # kL = pi + res = ph.expansion_chamber(f, length, 4.0 * s_duct, s_duct, speed_of_sound=c) + return numeric(0.0, float(res.transmission_loss[0]), 1e-9, unit="dB") + + +@register( + _NOISE_CONTROL, + "Bies 5e Eq. (8.44) / Example 8.1", + "Quarter-wave tube tuning f = c/(4 l_e), l_e = 1.516 m -> 56.6 Hz", +) +def _chk_quarter_wave_tuning() -> Outcome: + area = math.pi * 0.05**2 / 4.0 + res = ph.quarter_wave_resonator([100.0], area, 1.516, area, speed_of_sound=343.24) + assert res.resonances is not None + return numeric(56.6, float(res.resonances[0]), 0.1, unit="Hz", places=2) + + +@register( + _NOISE_CONTROL, + "Bies 5e Eq. (8.46)", + "Helmholtz resonance f0 = (c/2pi) sqrt(S/(l_e V)) (S=1e-4, l_e=0.02, V=1e-3)", +) +def _chk_helmholtz_resonance() -> Outcome: + c = 343.0 + res = ph.helmholtz_resonator([100.0], 0.01, 1e-4, 0.02, 1e-3, speed_of_sound=c) + assert res.resonances is not None + expected = c / (2.0 * math.pi) * math.sqrt(1e-4 / (0.02 * 1e-3)) + return numeric(expected, float(res.resonances[0]), 1e-6, unit="Hz", places=3) + + +@register( + _NOISE_CONTROL, + "Bies 5e Eq. (8.73)", + "Side-branch TL = 20 lg|1 + rho c/(2 Sd Zb)| (QWT branch, closed form)", +) +def _chk_side_branch_closed_form() -> Outcome: + from phonometry.noise_control import silencers as _sl + + f = np.array([120.0]) + zb = _sl.quarter_wave_impedance(f, 0.5, 0.002) + t = _sl.shunt_matrix(zb) + tl = float(_sl.transmission_loss(t, inlet_area=0.01, outlet_area=0.01)[0]) + closed = 20.0 * math.log10(abs(1.0 + 1.206 * 343.0 / (2.0 * 0.01 * zb[0]))) + return numeric(closed, tl, 1e-9, unit="dB") + + +@register( + _NOISE_CONTROL, + "Bies 5e Eq. (8.275) (Wells' plenum method)", + "Plenum TL = -10 lg[S_out(cos0/pi r^2 + (1-a)/(Sw a))] (S_out=.1,r=1,Sw=20,a=.2)", +) +def _chk_plenum_wells() -> Outcome: + tl = float(ph.plenum_attenuation(0.1, 1.0, 20.0, 0.2)) + direct = 1.0 / (math.pi * 1.0**2) + reverb = (1.0 - 0.2) / (20.0 * 0.2) + expected = -10.0 * math.log10(0.1 * (direct + reverb)) + return numeric(expected, tl, 1e-6, unit="dB") + + +@register( + _NOISE_CONTROL, + "Bies 5e Table 8.14 (ASHRAE end reflection, flush)", + "Duct end reflection D = 200 mm at 125 Hz = 10 dB (table node)", +) +def _chk_end_reflection_table() -> Outcome: + res = ph.end_reflection_loss([125.0], 0.200, termination="flush") + return numeric(10.0, float(res.values[0]), 1e-6, unit="dB") + + +@register( + _NOISE_CONTROL, + "Bies 5e Eqs. (7.103), (7.111) (enclosure, fully absorbing limit)", + "Enclosure correction C -> 10 lg 0.3 = -5.23 dB as alpha_i -> 1", +) +def _chk_enclosure_floor() -> Outcome: + res = ph.enclosure_insertion_loss([40.0], 6.0, 5.0, 0.999999) + return numeric(10.0 * math.log10(0.3), float(res.correction[0]), 1e-3, + unit="dB") # =========================================================================== diff --git a/site/src/content/docs/reference/api/index.md b/site/src/content/docs/reference/api/index.md index 93a41e24..99040fe6 100644 --- a/site/src/content/docs/reference/api/index.md +++ b/site/src/content/docs/reference/api/index.md @@ -183,6 +183,15 @@ La referencia de la API se genera a partir de los docstrings del código (en ing | [`electroacoustics.distortion`](/phonometry/reference/api/electroacoustics/distortion/) | Distortion metrics for electroacoustic equipment (IEC 60268-3 / AES17). | | [`electroacoustics.frequency_response`](/phonometry/reference/api/electroacoustics/frequency-response/) | Frequency-response and coherence estimators (Bendat & Piersol). | | [`electroacoustics.swept_sine`](/phonometry/reference/api/electroacoustics/swept-sine/) | Harmonic-distortion separation with exponential sweeps (Farina / Novak). | +| [`electroacoustics.piston`](/phonometry/reference/api/electroacoustics/piston/) | Radiation of a rigid circular piston set in an infinite baffle. | + +## Industrial noise control + +| Module | Summary | +| :--- | :--- | +| [`noise_control.silencers`](/phonometry/reference/api/noise_control/silencers/) | Reactive silencers by the four-pole (transmission-matrix) method. | +| [`noise_control.hvac`](/phonometry/reference/api/noise_control/hvac/) | HVAC duct acoustics: end reflection, bends, plenums and flow-generated noise. | +| [`noise_control.enclosures`](/phonometry/reference/api/noise_control/enclosures/) | Insertion loss of a close or free-standing machine enclosure. | ## Program loudness diff --git a/site/src/generated/api-sidebar.mjs b/site/src/generated/api-sidebar.mjs index bba59195..c1025a42 100644 --- a/site/src/generated/api-sidebar.mjs +++ b/site/src/generated/api-sidebar.mjs @@ -179,6 +179,16 @@ export const apiSidebar = { 'reference/api/electroacoustics/distortion', 'reference/api/electroacoustics/frequency-response', 'reference/api/electroacoustics/swept-sine', + 'reference/api/electroacoustics/piston', + ], + }, + { + label: 'Industrial noise control', + translations: { es: 'Control de ruido industrial' }, + items: [ + 'reference/api/noise_control/silencers', + 'reference/api/noise_control/hvac', + 'reference/api/noise_control/enclosures', ], }, { diff --git a/src/phonometry/__init__.py b/src/phonometry/__init__.py index 0b30156b..bf8fa342 100644 --- a/src/phonometry/__init__.py +++ b/src/phonometry/__init__.py @@ -309,6 +309,13 @@ coherence, transfer_function, ) +from .electroacoustics.piston import ( + RadiatingPistonResult, + piston_directivity, + piston_reactance, + piston_resistance, + radiating_piston, +) from .electroacoustics.swept_sine import ( SweptSineDistortionResult, swept_sine_distortion, @@ -918,6 +925,22 @@ SignalSource, fdtd_simulation, ) +from .noise_control.silencers import ( + ReactiveSilencerResult, + expansion_chamber, + extended_tube_chamber, + helmholtz_resonator, + quarter_wave_resonator, +) +from .noise_control.hvac import ( + HvacSpectrumResult, + elbow_insertion_loss, + end_reflection_loss, + flow_noise_bend, + flow_noise_straight_duct, + plenum_attenuation, +) +from .noise_control.enclosures import EnclosureResult, enclosure_insertion_loss from ._version import __version__ # Public methods @@ -993,6 +1016,11 @@ "swept_sine_distortion", "synchronized_sweep_signal", "SweptSineDistortionResult", + "radiating_piston", + "piston_resistance", + "piston_reactance", + "piston_directivity", + "RadiatingPistonResult", "program_loudness", "integrated_loudness", "loudness_range", @@ -1694,6 +1722,23 @@ "circular_aperture_transmission_coefficient", "composite_transmission_loss", "transmission_loss_from_coefficient", + # Industrial noise control (silencers, HVAC, enclosures). The low-level + # four-pole primitives (duct_matrix, shunt_matrix, cascade, transmission_loss, + # insertion_loss, ...) stay under phonometry.noise_control to keep the top + # level clean and avoid clashing with the underwater transmission_loss. + "expansion_chamber", + "helmholtz_resonator", + "quarter_wave_resonator", + "extended_tube_chamber", + "ReactiveSilencerResult", + "end_reflection_loss", + "elbow_insertion_loss", + "plenum_attenuation", + "flow_noise_straight_duct", + "flow_noise_bend", + "HvacSpectrumResult", + "enclosure_insertion_loss", + "EnclosureResult", ] From bedcd830dd11e34a910fbf58a4b2843467b1c5d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20M=2E=20Requena=20Plens?= Date: Sun, 19 Jul 2026 17:15:31 +0200 Subject: [PATCH 4/5] fix(noise_control): correct insertion-loss sign and straight-duct flow-noise constant Adversarial review and a source re-check turned up two numerical errors: - insertion_loss returned the reciprocal of the correct ratio, so a real silencer reported a negative insertion loss. Use IL = 20 lg|(T11 Zr + T12 + Zs Zr T21 + Zs T22)/(Zs + Zr)|, which is 0 for a through connection and equals the transmission loss for the anechoic reference Zs = rho c/S_in, Zr = rho c/S_out. Added a test pinning IL = TL for that reference (the previous identity-only test could not see the sign) and a conformance check. - flow_noise_straight_duct dropped the -2 constant of Bies Eq. (8.251) (VDI 2081-1); added it and pinned the numeric value. Also harden shunt_matrix against the divide-by-zero at an exact lossless resonance, and correct the flow_noise_bend docstring: the radiated power grows as U^6/U^8 (dipole/quadrupole) while only the stream-power efficiency grows as U^3/U^5 (the formula itself already matched Bies Eqs. (8.252)/(8.254)). --- docs/CONFORMANCE.md | 5 +++-- scripts/conformance_report.py | 18 ++++++++++++++++++ .../docs/reference/api/noise_control/hvac.md | 13 ++++++++----- .../reference/api/noise_control/silencers.md | 6 ++++-- src/phonometry/noise_control/hvac.py | 14 +++++++++----- src/phonometry/noise_control/silencers.py | 19 ++++++++++++------- tests/noise_control/test_hvac.py | 5 ++++- tests/noise_control/test_silencers.py | 13 +++++++++++++ 8 files changed, 71 insertions(+), 22 deletions(-) diff --git a/docs/CONFORMANCE.md b/docs/CONFORMANCE.md index 63cb2bdb..29ebd668 100644 --- a/docs/CONFORMANCE.md +++ b/docs/CONFORMANCE.md @@ -15,7 +15,7 @@ ## Numerical conformance report -✅ **338/338 conformance checks pass** across 42 domains and 211 standards - filters class 1 - weightings within IEC 61672-1 class 1. +✅ **339/339 conformance checks pass** across 42 domains and 212 standards - filters class 1 - weightings within IEC 61672-1 class 1. Each row pins a standard clause to its expected normative value and the value the library computes. Every section below is collapsible and stays collapsed while all of its rows pass; a section with any failing row opens automatically. @@ -704,7 +704,7 @@ Only **Butterworth** (the library default) and **Chebyshev-II** are class-compli
-✅ Industrial noise control: 100% (8/8) +✅ Industrial noise control: 100% (9/9) | Standard | Quantity | Expected (norm) | Computed | Δ | Status | |:---|:---|:---|:---|:---|:---:| @@ -713,6 +713,7 @@ Only **Butterworth** (the library default) and **Chebyshev-II** are class-compli | Bies 5e Eq. (8.44) / Example 8.1 | Quarter-wave tube tuning f = c/(4 l_e), l_e = 1.516 m -> 56.6 Hz | 56.6 Hz (+/-0.1 Hz) | 56.6 Hz | 0.003 Hz | ✅ | | Bies 5e Eq. (8.46) | Helmholtz resonance f0 = (c/2pi) sqrt(S/(l_e V)) (S=1e-4, l_e=0.02, V=1e-3) | 122.067 Hz (+/-0 Hz) | 122.067 Hz | 0 Hz | ✅ | | Bies 5e Eq. (8.73) | Side-branch TL = 20 lg|1 + rho c/(2 Sd Zb)| (QWT branch, closed form) | 0.1638 dB (+/-0 dB) | 0.1638 dB | 0 dB | ✅ | +| Bies 5e Eqs. (8.141)/(8.148) (four-pole insertion loss) | Insertion loss = transmission loss for the anechoic reference Zs=Zr=rho c/S | 6.2498 dB (= TL) | 6.2498 dB | 0 dB | ✅ | | Bies 5e Eq. (8.275) (Wells' plenum method) | Plenum TL = -10 lg[S_out(cos0/pi r^2 + (1-a)/(Sw a))] (S_out=.1,r=1,Sw=20,a=.2) | 12.8541 dB (+/-0 dB) | 12.8541 dB | 0 dB | ✅ | | Bies 5e Table 8.14 (ASHRAE end reflection, flush) | Duct end reflection D = 200 mm at 125 Hz = 10 dB (table node) | 10 dB (+/-0 dB) | 10 dB | 0 dB | ✅ | | Bies 5e Eqs. (7.103), (7.111) (enclosure, fully absorbing limit) | Enclosure correction C -> 10 lg 0.3 = -5.23 dB as alpha_i -> 1 | -5.2288 dB (+/-0.001 dB) | -5.2288 dB | 0 dB | ✅ | diff --git a/scripts/conformance_report.py b/scripts/conformance_report.py index 4818bec0..cafee6fc 100644 --- a/scripts/conformance_report.py +++ b/scripts/conformance_report.py @@ -5882,6 +5882,24 @@ def _chk_side_branch_closed_form() -> Outcome: return numeric(closed, tl, 1e-9, unit="dB") +@register( + _NOISE_CONTROL, + "Bies 5e Eqs. (8.141)/(8.148) (four-pole insertion loss)", + "Insertion loss = transmission loss for the anechoic reference Zs=Zr=rho c/S", +) +def _chk_insertion_loss_equals_tl() -> Outcome: + from phonometry.noise_control import silencers as _sl + + c, rho, s = 343.0, 1.206, 0.01 + z = rho * c / s + f = np.array([232.0]) # a frequency with a substantial, positive TL + t = _sl.expansion_chamber(f, 0.3, 0.04, s).transfer_matrix + tl = float(_sl.transmission_loss(t, inlet_area=s, outlet_area=s)[0]) + il = float(_sl.insertion_loss(t, source_impedance=z, radiation_impedance=z)[0]) + return numeric(tl, il, 1e-9, unit="dB", + expected_label=f"{tl:.4f} dB (= TL)") + + @register( _NOISE_CONTROL, "Bies 5e Eq. (8.275) (Wells' plenum method)", diff --git a/site/src/content/docs/reference/api/noise_control/hvac.md b/site/src/content/docs/reference/api/noise_control/hvac.md index 55877be3..4b9123ee 100644 --- a/site/src/content/docs/reference/api/noise_control/hvac.md +++ b/site/src/content/docs/reference/api/noise_control/hvac.md @@ -109,9 +109,12 @@ Flow-generated octave-band sound power of a mitred bend (Bies Eqs. (8.252), (8.2 `L_WB = L_Ws - 10 log10(1 + 0.165 N_s^2) + 30 log10(U) - 103` with the stream power level `L_Ws = 30 log10(U) + 10 log10(S) + 10 log10(rho) + 117` -and the Strouhal number `N_s = f H / U` (`H` the duct height in the -plane of the bend). The regeneration scales as `U^3` at low `N_s` and -`U^5` at high `N_s`. +(Bies Eq. (8.252)) and the Strouhal number `N_s = f H / U` (`H` the duct +height in the plane of the bend). The radiated sound power grows as the +sixth power of the stream speed at low `N_s` (the inner-corner drag +dipole) and the eighth power at high `N_s` (the outer-corner shear +quadrupole); equivalently, the *efficiency* referenced to the stream power +grows as `U^3` and `U^5` respectively. **Parameters** @@ -137,8 +140,8 @@ flow_noise_straight_duct( Flow-generated octave-band sound power of a straight duct (Bies Eq. (8.251)). -`L_WB = 7 + 50 log10(U) + 10 log10(S) - 26 log10(1.14 + 0.02 f / U)` in -dB re 1e-12 W (VDI 2081-1), for airflow speed `U` in a duct of area +`L_WB = 7 + 50 log10(U) + 10 log10(S) - 2 - 26 log10(1.14 + 0.02 f / U)` +in dB re 1e-12 W (VDI 2081-1), for airflow speed `U` in a duct of area `S`. **Parameters** diff --git a/site/src/content/docs/reference/api/noise_control/silencers.md b/site/src/content/docs/reference/api/noise_control/silencers.md index 789a34d5..506c48d4 100644 --- a/site/src/content/docs/reference/api/noise_control/silencers.md +++ b/site/src/content/docs/reference/api/noise_control/silencers.md @@ -42,9 +42,11 @@ internal impedance `Z_s` radiating into a termination impedance `Z_r` is the extra attenuation of inserting the silencer in place of a direct connection, - IL = 20 log10 | (Z_s + Z_r) / (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) |, + IL = 20 log10 | (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) / (Z_s + Z_r) |, -which is `0` when the silencer reduces to a through connection (`T = I`). +which is `0` when the silencer reduces to a through connection (`T = I`) +and equals the transmission loss for the anechoic reference +`Z_s = rho c / S_in`, `Z_r = rho c / S_out`. **Simple expansion chamber.** A chamber of area `S_exp` and length `L` between pipes of area `S_duct` has the closed-form transmission loss (Bies diff --git a/src/phonometry/noise_control/hvac.py b/src/phonometry/noise_control/hvac.py index 4db4c9f3..9717a5c8 100644 --- a/src/phonometry/noise_control/hvac.py +++ b/src/phonometry/noise_control/hvac.py @@ -280,8 +280,8 @@ def flow_noise_straight_duct( ) -> HvacSpectrumResult: """Flow-generated octave-band sound power of a straight duct (Bies Eq. (8.251)). - ``L_WB = 7 + 50 log10(U) + 10 log10(S) - 26 log10(1.14 + 0.02 f / U)`` in - dB re 1e-12 W (VDI 2081-1), for airflow speed ``U`` in a duct of area + ``L_WB = 7 + 50 log10(U) + 10 log10(S) - 2 - 26 log10(1.14 + 0.02 f / U)`` + in dB re 1e-12 W (VDI 2081-1), for airflow speed ``U`` in a duct of area ``S``. :param frequencies: Octave-band centre frequencies ``f``, Hz (1-D array). @@ -297,6 +297,7 @@ def flow_noise_straight_duct( 7.0 + 50.0 * np.log10(u) + 10.0 * np.log10(s) + - 2.0 - 26.0 * np.log10(1.14 + 0.02 * f / u) ) return HvacSpectrumResult( @@ -317,9 +318,12 @@ def flow_noise_bend( ``L_WB = L_Ws - 10 log10(1 + 0.165 N_s^2) + 30 log10(U) - 103`` with the stream power level ``L_Ws = 30 log10(U) + 10 log10(S) + 10 log10(rho) + 117`` - and the Strouhal number ``N_s = f H / U`` (``H`` the duct height in the - plane of the bend). The regeneration scales as ``U^3`` at low ``N_s`` and - ``U^5`` at high ``N_s``. + (Bies Eq. (8.252)) and the Strouhal number ``N_s = f H / U`` (``H`` the duct + height in the plane of the bend). The radiated sound power grows as the + sixth power of the stream speed at low ``N_s`` (the inner-corner drag + dipole) and the eighth power at high ``N_s`` (the outer-corner shear + quadrupole); equivalently, the *efficiency* referenced to the stream power + grows as ``U^3`` and ``U^5`` respectively. :param frequencies: Octave-band centre frequencies ``f``, Hz (1-D array). :param flow_velocity: Mean flow speed ``U``, m/s. diff --git a/src/phonometry/noise_control/silencers.py b/src/phonometry/noise_control/silencers.py index 9dee2942..480dae7d 100644 --- a/src/phonometry/noise_control/silencers.py +++ b/src/phonometry/noise_control/silencers.py @@ -35,9 +35,11 @@ the extra attenuation of inserting the silencer in place of a direct connection, - IL = 20 log10 | (Z_s + Z_r) / (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) |, + IL = 20 log10 | (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) / (Z_s + Z_r) |, -which is ``0`` when the silencer reduces to a through connection (``T = I``). +which is ``0`` when the silencer reduces to a through connection (``T = I``) +and equals the transmission loss for the anechoic reference +``Z_s = rho c / S_in``, ``Z_r = rho c / S_out``. **Simple expansion chamber.** A chamber of area ``S_exp`` and length ``L`` between pipes of area ``S_duct`` has the closed-form transmission loss (Bies @@ -127,7 +129,10 @@ def shunt_matrix(branch_impedance: ArrayLike) -> _Complex: raise ValueError("'branch_impedance' must be a non-empty 1-D array.") t = np.zeros((zb.size, 2, 2), dtype=np.complex128) t[:, 0, 0] = 1.0 - t[:, 1, 0] = 1.0 / zb + # A lossless branch at exact resonance has Z_b -> 0 (it shorts the duct); + # 1/Z_b -> infinity gives the ideal infinite attenuation there. + with np.errstate(divide="ignore", invalid="ignore"): + t[:, 1, 0] = 1.0 / zb t[:, 1, 1] = 1.0 return t @@ -196,8 +201,8 @@ def insertion_loss( length) connection between a source of internal impedance ``Z_s`` and a radiation (termination) impedance ``Z_r``: - IL = 20 log10 | (Z_s + Z_r) / - (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) |. + IL = 20 log10 | (T11 Z_r + T12 + Z_s Z_r T21 + Z_s T22) / + (Z_s + Z_r) |. :param transfer_matrix: A ``(n_freq, 2, 2)`` compound matrix. :param source_impedance: Source internal acoustic impedance ``Z_s``, @@ -217,9 +222,9 @@ def insertion_loss( t12 = transfer_matrix[:, 0, 1] t21 = transfer_matrix[:, 1, 0] t22 = transfer_matrix[:, 1, 1] - denom = t11 * zr + t12 + zs * zr * t21 + zs * t22 + num = t11 * zr + t12 + zs * zr * t21 + zs * t22 return np.asarray( - 20.0 * np.log10(np.abs((zs + zr) / denom)), dtype=np.float64 + 20.0 * np.log10(np.abs(num / (zs + zr))), dtype=np.float64 ) diff --git a/tests/noise_control/test_hvac.py b/tests/noise_control/test_hvac.py index 878f3ded..7c447ce2 100644 --- a/tests/noise_control/test_hvac.py +++ b/tests/noise_control/test_hvac.py @@ -84,14 +84,17 @@ def test_plenum_per_band() -> None: def test_flow_noise_straight_duct_formula() -> None: + # Bies Eq. (8.251), VDI 2081-1: note the -2 constant term. f = np.array([250.0]) u, s = 10.0, 0.04 res = hvac.flow_noise_straight_duct(f, u, s) expected = ( - 7.0 + 50.0 * math.log10(u) + 10.0 * math.log10(s) + 7.0 + 50.0 * math.log10(u) + 10.0 * math.log10(s) - 2.0 - 26.0 * math.log10(1.14 + 0.02 * 250.0 / u) ) assert res.values[0] == pytest.approx(expected) + # Pinned numeric anchor so the constant cannot silently drift. + assert res.values[0] == pytest.approx(35.4347, abs=1e-3) def test_flow_noise_scales_with_velocity() -> None: diff --git a/tests/noise_control/test_silencers.py b/tests/noise_control/test_silencers.py index 46fbacb8..e98b5918 100644 --- a/tests/noise_control/test_silencers.py +++ b/tests/noise_control/test_silencers.py @@ -92,6 +92,19 @@ def test_insertion_loss_identity_is_zero() -> None: assert np.allclose(il, 0.0, atol=1e-12) +def test_insertion_loss_equals_tl_for_anechoic_reference() -> None: + # The insertion loss with Z_s = rho c / S_in and Z_r = rho c / S_out equals + # the (anechoic) transmission loss; a positive, not negative, quantity. + c, rho, s = 343.0, 1.206, 0.01 + z = rho * c / s + f = np.linspace(50.0, 1500.0, 400) + t = sl.expansion_chamber(f, 0.3, 0.04, s).transfer_matrix + tl = sl.transmission_loss(t, inlet_area=s, outlet_area=s) + il = sl.insertion_loss(t, source_impedance=z, radiation_impedance=z) + assert np.allclose(il, tl, atol=1e-9) + assert il.max() > 3.0 # a real silencer gives a positive insertion loss + + def test_transfer_matrix_reciprocity() -> None: # Reciprocal passive elements have det(T) = 1. f = np.linspace(50.0, 500.0, 20) From 82b8402964a80770f345ae9e51caa5869c59ac18 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20M=2E=20Requena=20Plens?= Date: Sun, 19 Jul 2026 19:02:20 +0200 Subject: [PATCH 5/5] fix(noise_control,electroacoustics): address Sonar and bot review findings Sonar: - S1244 (piston): replace exact float == 0 tests with finite-masking at the removable 0/0 singularities of R1, X1 and the directivity (physics and all oracles unchanged). - S3776 (enclosures): split enclosure_insertion_loss into _resolve_frequencies and _resolve_panel_r helpers with identical semantics. - S5863 (FDTD cross-check test): compare two independently computed runs instead of a self-comparison, keeping the determinism check meaningful. Review bots (verified against Bies / Beranek & Mellow before applying): - Vectorise the piston directivity-index quadrature over frequency (bit-identical results, confirmed by the byte-stable conformance report). - Validate that a silencer chain shares one frequency grid (cascade) and that an extended tube's inlet+outlet extensions do not exceed the chamber length. - Reject a multi-dimensional interior absorption in enclosures. - Use kwargs.setdefault for the plot line styles so a caller-supplied lw/marker/ms no longer raises, and keep the resonance legend label when the first resonance is out of range. - Rephrase the side-branch conformance description with abs(...) so the pipes no longer break the Markdown table, and add a text language hint to the formula code blocks of the noise-control doc. Rejected: converting the enclosures docstring equation blocks to inline literals (the indented-equation-block style is the library-wide docstring convention, used by ~120 modules; changing one would be inconsistent). --- docs/CONFORMANCE.md | 34 +++++++++++- docs/noise-control.md | 8 +-- llms-full.txt | 14 +++++ scripts/conformance_report.py | 2 +- src/phonometry/_plot/noise_control.py | 20 ++++--- src/phonometry/electroacoustics/piston.py | 28 +++++----- src/phonometry/noise_control/enclosures.py | 60 +++++++++++++-------- src/phonometry/noise_control/silencers.py | 10 ++++ tests/noise_control/test_fdtd_crosscheck.py | 5 +- 9 files changed, 133 insertions(+), 48 deletions(-) diff --git a/docs/CONFORMANCE.md b/docs/CONFORMANCE.md index 29ebd668..e535f374 100644 --- a/docs/CONFORMANCE.md +++ b/docs/CONFORMANCE.md @@ -15,7 +15,7 @@ ## Numerical conformance report -✅ **339/339 conformance checks pass** across 42 domains and 212 standards - filters class 1 - weightings within IEC 61672-1 class 1. +✅ **353/353 conformance checks pass** across 44 domains and 224 standards - filters class 1 - weightings within IEC 61672-1 class 1. Each row pins a standard clause to its expected normative value and the value the library computes. Every section below is collapsible and stays collapsed while all of its rows pass; a section with any failing row opens automatically. @@ -689,6 +689,36 @@ Only **Butterworth** (the library default) and **Chebyshev-II** are class-compli
+
+✅ Panel & aperture sound insulation (Bies / Hopkins / Cremer): 100% (11/11) + +| Standard | Quantity | Expected (norm) | Computed | Δ | Status | +|:---|:---|:---|:---|:---|:---:| +| Bies 5e Eq. 7.40 (mass law) | 6 dB per octave (500 -> 1000 Hz) | 6.0206 dB (+/-0.01 dB) | 6.02 dB | -0.001 dB | ✅ | +| Bies 5e Eq. 7.40 (mass law) | 6 dB per doubling of mass | 6.0206 dB (+/-0.01 dB) | 6.02 dB | -0.001 dB | ✅ | +| Bies 5e Eq. 7.42 (field incidence) | One-third-octave correction 5.5 dB | 5.5 dB (+/-0.001 dB) | 5.5 dB | 0 dB | ✅ | +| Hopkins Eq. 2.201 / Bies Eq. 7.3 | Coincidence frequency, 6 mm glass | 2079 Hz (+/-3%) | 2107.3639 Hz | 28.364 Hz | ✅ | +| Cremer Table 5.1 | Thin-plate point impedance Z = 8 sqrt(B' m'') | 2529.8221 N.s/m (+/-0 N.s/m) | 2529.8221 N.s/m | 0 N.s/m | ✅ | +| Cremer Table 5.1 | Infinite-beam mobility phase -45 deg | -45 deg (+/-0 deg) | -45 deg | 0 deg | ✅ | +| Hopkins Eq. 2.229 (Leppington/Maidanik) | Radiation efficiency at f = 2 fc | 1.4142 (+/-0) | 1.4142 | 0 | ✅ | +| Bies Eq. 7.62 / Hopkins Eq. 4.73 | Mass-air-mass resonance f0, empty cavity | 76.9484 Hz (+/-0.5%) | 76.8521 Hz | -0.096 Hz | ✅ | +| Bies Eq. 7.64 (double wall) | Below f0 = mass law of the combined mass | 11.6144 dB (+/-0 dB) | 11.6144 dB | 0 dB | ✅ | +| Hopkins Eq. 4.92 (composite) | 1 % open area caps R at 10 lg(S/Sa) | 20 dB (+/-0.05 dB) | 19.9996 dB | 0 dB | ✅ | +| Hopkins Eq. 4.99/4.101 (Gomperts slit) | Transmission maximum at first resonance | 1544.9615 Hz (+/-15 Hz) | 1542.9615 Hz | -2 Hz | ✅ | + +
+ +
+✅ Atmospheric refraction (Salomons rays / GFPE): 100% (3/3) + +| Standard | Quantity | Expected (norm) | Computed | Δ | Status | +|:---|:---|:---|:---|:---|:---:| +| Salomons Sec. 4.4 (ray turning height, linear profile) | Turning height of a 10 deg ray vs Rc(1 - cos theta0) (circular arc), m | 26.457 m (+/-0.1 m) | 26.457 m | 0 m | ✅ | +| Salomons Eq. (3.4) (GFPE vs spherical-wave ground effect, homogeneous) | PE relative level at 500 m over grassland vs Weyl-Van der Pol, dB | -13.919 dB (+/-0.5 dB) | -13.842 dB | 0.077 dB | ✅ | +| Salomons Eq. (3.4) (GFPE hard ground vs two-ray, homogeneous) | PE relative level at 500 m over a rigid ground vs the coherent two-ray, dB | 5.997 dB (+/-0.6 dB) | 5.593 dB | -0.405 dB | ✅ | + +
+
✅ Electroacoustics: 100% (6/6) @@ -712,7 +742,7 @@ Only **Butterworth** (the library default) and **Chebyshev-II** are class-compli | Bies 5e Eq. (8.111) | Expansion-chamber trough TL = 0 at kL = pi (chamber transparent) | 0 dB (+/-0 dB) | 0 dB | 0 dB | ✅ | | Bies 5e Eq. (8.44) / Example 8.1 | Quarter-wave tube tuning f = c/(4 l_e), l_e = 1.516 m -> 56.6 Hz | 56.6 Hz (+/-0.1 Hz) | 56.6 Hz | 0.003 Hz | ✅ | | Bies 5e Eq. (8.46) | Helmholtz resonance f0 = (c/2pi) sqrt(S/(l_e V)) (S=1e-4, l_e=0.02, V=1e-3) | 122.067 Hz (+/-0 Hz) | 122.067 Hz | 0 Hz | ✅ | -| Bies 5e Eq. (8.73) | Side-branch TL = 20 lg|1 + rho c/(2 Sd Zb)| (QWT branch, closed form) | 0.1638 dB (+/-0 dB) | 0.1638 dB | 0 dB | ✅ | +| Bies 5e Eq. (8.73) | Side-branch TL = 20 lg abs(1 + rho c/(2 Sd Zb)) (QWT branch, closed form) | 0.1638 dB (+/-0 dB) | 0.1638 dB | 0 dB | ✅ | | Bies 5e Eqs. (8.141)/(8.148) (four-pole insertion loss) | Insertion loss = transmission loss for the anechoic reference Zs=Zr=rho c/S | 6.2498 dB (= TL) | 6.2498 dB | 0 dB | ✅ | | Bies 5e Eq. (8.275) (Wells' plenum method) | Plenum TL = -10 lg[S_out(cos0/pi r^2 + (1-a)/(Sw a))] (S_out=.1,r=1,Sw=20,a=.2) | 12.8541 dB (+/-0 dB) | 12.8541 dB | 0 dB | ✅ | | Bies 5e Table 8.14 (ASHRAE end reflection, flush) | Duct end reflection D = 200 mm at 125 Hz = 10 dB (table node) | 10 dB (+/-0 dB) | 10 dB | 0 dB | ✅ | diff --git a/docs/noise-control.md b/docs/noise-control.md index 9fb65023..f399cfda 100644 --- a/docs/noise-control.md +++ b/docs/noise-control.md @@ -20,7 +20,7 @@ two ends (Bies Eq. (8.133); Munjal, *Acoustics of Ducts and Mufflers*), and a compound silencer is the ordered matrix product of its elements. A straight duct of length `L` and area `S` is (Bies Eq. (8.143), no flow) -``` +```text [ cos(kL) j (rho c / S) sin(kL) ] [ j (S / rho c) sin(kL) cos(kL) ] , k = omega / c, ``` @@ -30,7 +30,7 @@ and a side branch of acoustic impedance `Z_b` is the shunt the compound matrix `T` (Bies Eq. (8.141); for equal inlet/outlet areas, Eq. (8.148)) -``` +```text TL = 20 log10( (1/2) | T11 + T12/Zc + Zc T21 + T22 | ) , Zc = rho c / S, ``` @@ -44,7 +44,7 @@ A chamber of area `S_exp` and length `L` between pipes of area `S_duct` has the closed-form transmission loss (Bies Eq. (8.111)) with area ratio `m = S_exp / S_duct`: -``` +```text TL = 10 log10[ 1 + (1/4) (m - 1/m)^2 sin^2(kL) ] , ``` @@ -118,7 +118,7 @@ A sealed enclosure reduces the radiated noise by its panel transmission loss `R`, minus a penalty `C` for the reverberant build-up inside the small, hard cavity (Bies Eqs. (7.103), (7.111)): -``` +```text IL = R - C , C = 10 log10( 0.3 + S_E / R_i ) , ``` diff --git a/llms-full.txt b/llms-full.txt index 71e7ac2e..56db3ab5 100644 --- a/llms-full.txt +++ b/llms-full.txt @@ -10246,6 +10246,20 @@ cycle and warn on use; they are removed in 4.0. | `synchronized_sweep_signal` | `function` | **Synchronized exponential sweep (Novak et al. 2015, Eqs. 47/49).**
• `fs`, `f1` < `f2` [Hz], `seconds` (quantized by the rounding)
• `amplitude` (Default: 1.0), `fade` (Default: 0.0)
Rate L = round(f1·T̃/ln(f2/f1))/f1 → f1·L integer, so a L·ln(n) shift equals the nth harmonic and harmonic phases become system properties | `x = synchronized_sweep_signal(48000, 20.0, 6000.0, 4.0)`

• sweep samples (duration L·ln(f2/f1)) | | `swept_sine_distortion` | `function` | **Harmonic separation and THD(f) from one sweep (Farina 2000 / Novak et al. 2015).**
• `recorded`, `fs`, plus the generator parameters `f1`, `f2`, `seconds`
• `method`: 'synchronized' (Default; analytic deconvolution, meaningful phases) / 'farina' (classical ESS, magnitudes only)
• `n_harmonics` (Default: 5), `ir_length` (Default: largest power of two fitting the closest arrivals)
• `amplitude` (Default: 1.0), `fade` (Default: generator default), `remove_dc` (Default: True) | `res = swept_sine_distortion(y, fs, 20.0, 6000.0, 4.0)`

• `SweptSineDistortionResult` | | `SweptSineDistortionResult` | `dataclass` | **Swept-sine harmonic separation.**
• `frequencies` [Hz], `harmonic_responses`: complex H₁..H_N, `harmonic_irs` (windowed, centred), `delays` L·ln(n) [s]
• `thd_frequencies` [Hz], `thd`, `distortion_ratios`: \|Hₙ(n·f)\|/\|H₁(f)\| per order
• `fs`, `f1`, `f2`, `duration`, `rate` L, `method`, `n_harmonics`
• `.plot()`: \|Hₙ\| magnitudes + THD(f) | `res.thd, res.harmonic_responses` | +| `radiating_piston` | `function` | **Radiation of a rigid baffled circular piston (Beranek & Mellow §4.19/§13.7).**
• `radius` a [m], `frequencies` [Hz]
• `speed_of_sound` (Default: 343), `density` (Default: 1.206)
• `angles`: directivity polar angles [rad] (Default: None) | `res = radiating_piston(0.1, freqs)`

• `RadiatingPistonResult` | +| `piston_resistance` / `piston_reactance` | `function` | **Normalized piston radiation impedance R1 = 1−2J1(x)/x, X1 = 2H1(x)/x, x = 2ka (Beranek & Mellow Eqs. 13.117/13.118).** | `piston_resistance(2.0) # 0.4233` | +| `piston_directivity` | `function` | **Far-field piston directivity 2·J1(ka·sinθ)/(ka·sinθ).**
• `ka`, `theta` [rad] | `piston_directivity(5.0, 0.3)` | +| `RadiatingPistonResult` | `dataclass` | **Piston radiation result.**
• `ka`, `resistance`/`reactance`, `radiation_resistance`/`radiation_reactance` [N·s/m], `radiation_mass` = 8ρa³/3 [kg], `directivity_index` [dB], `directivity`
• `.plot()`: R1/X1 vs ka | `res.radiation_mass` | +| `expansion_chamber` | `function` | **Expansion-chamber silencer TL/IL (Bies Eq. 8.111, four-pole).**
• `frequencies`, `length` L, `chamber_area`, `pipe_area`
• `source_impedance`/`radiation_impedance` (Default: None → no IL) | `res = expansion_chamber(f, 0.3, 0.04, 0.01)`

• `ReactiveSilencerResult` | +| `helmholtz_resonator` / `quarter_wave_resonator` | `function` | **Side-branch resonator silencers (Bies Eqs. 8.46/8.44).**
• Helmholtz: `duct_area`, `neck_area`, `neck_length`, `cavity_volume`
• Quarter-wave: `duct_area`, `length`, `branch_area` | `quarter_wave_resonator(f, 0.01, 1.5, 2e-3)` | +| `extended_tube_chamber` | `function` | **Extended-inlet/outlet expansion chamber (Bies §8.9.7); reduces to `expansion_chamber` at zero extension.**
• `length`, `chamber_area`, `pipe_area`, `inlet_extension`, `outlet_extension` | `extended_tube_chamber(f, 0.4, 0.04, 0.01, inlet_extension=0.1)` | +| `ReactiveSilencerResult` | `dataclass` | **Reactive-silencer result.**
• `frequencies`, `transmission_loss`, `insertion_loss` (or None), `transfer_matrix`, `kind`, `resonances`
• `.plot()`: TL/IL vs frequency | `res.transmission_loss` | +| `end_reflection_loss` / `elbow_insertion_loss` | `function` | **HVAC duct end reflection & bend insertion loss (Bies Tables 8.14/8.11, ASHRAE).**
• end: `frequencies`, `diameter`, `termination` 'flush'/'free'
• elbow: `frequencies`, `width`, `bend_type`, `vanes`, `lined` | `end_reflection_loss(bands, 0.3)`

• `HvacSpectrumResult` | +| `plenum_attenuation` | `function` | **Plenum-chamber TL by Wells' method (Bies Eq. 8.275).**
• `exit_area`, `line_of_sight`, `wall_area`, `mean_absorption`, `angle` | `plenum_attenuation(0.1, 1.0, 20.0, 0.2) # dB` | +| `flow_noise_straight_duct` / `flow_noise_bend` | `function` | **HVAC flow-generated (self) noise sound power (VDI 2081, Bies Eqs. 8.251/8.254).**
• `frequencies`, `flow_velocity`, `area` (+ `height` for the bend) | `flow_noise_straight_duct(bands, 10.0, 0.04)`

• `HvacSpectrumResult` | +| `HvacSpectrumResult` | `dataclass` | **Per-frequency HVAC attenuation or regenerated Lw.**
• `frequencies`, `values`, `quantity`, `label`
• `.plot()` | `res.values` | +| `enclosure_insertion_loss` | `function` | **Machine-enclosure insertion loss IL = R − C (Bies Eqs. 7.103/7.111); panel R supplied by the caller, never predicted.**
• `panel_transmission_loss` (per-band array or callable of f)
• `external_area`, `internal_area`, `internal_absorption`
• `frequencies` (required for a callable R) | `enclosure_insertion_loss(R, 6.0, 5.0, 0.3)`

• `EnclosureResult` | +| `EnclosureResult` | `dataclass` | **Enclosure insertion-loss result.**
• `panel_transmission_loss`, `correction` C, `insertion_loss`, `room_constant`, `external_area`, `internal_area`
• `.plot()` | `res.insertion_loss` | | `program_loudness` | `function` | **EBU R 128 measurement set of a programme (BS.1770-5 + Tech 3341/3342).**
• `x`: signal (1D mono or 2D [channels, samples], 1.0 = 0 dBFS)
• `fs` [Hz]
• `weights`: per-channel Gi (Default: Table 3 by channel count)
• `momentary_step` [s] (Default: 0.01), `short_term_step` [s] (Default: 0.1)
• `oversample`: true-peak factor (Default: None → reaches 192 kHz) | `res = program_loudness(x, fs)`

• `ProgramLoudnessResult` | | `integrated_loudness` | `function` | **Programme (integrated) loudness with the two-stage gate (BS.1770-5 Annex 1).**
• `x`, `fs` [Hz]
• `weights` (Default: Table 3 by channel count)
Absolute gate −70 LKFS, relative −10 LU | `integrated_loudness(x, fs) # LUFS` | | `loudness_range` | `function` | **Loudness range LRA (EBU Tech 3342).**
• `short_term_loudness`: S readings [LUFS] at ≥ 10 Hz
Cascaded gate (−70 LUFS abs, −20 LU rel), 10th-95th percentile spread | `loudness_range(res.short_term) # LU` | diff --git a/scripts/conformance_report.py b/scripts/conformance_report.py index cafee6fc..1376e9e5 100644 --- a/scripts/conformance_report.py +++ b/scripts/conformance_report.py @@ -5869,7 +5869,7 @@ def _chk_helmholtz_resonance() -> Outcome: @register( _NOISE_CONTROL, "Bies 5e Eq. (8.73)", - "Side-branch TL = 20 lg|1 + rho c/(2 Sd Zb)| (QWT branch, closed form)", + "Side-branch TL = 20 lg abs(1 + rho c/(2 Sd Zb)) (QWT branch, closed form)", ) def _chk_side_branch_closed_form() -> Outcome: from phonometry.noise_control import silencers as _sl diff --git a/src/phonometry/_plot/noise_control.py b/src/phonometry/_plot/noise_control.py index 026bce27..1179142a 100644 --- a/src/phonometry/_plot/noise_control.py +++ b/src/phonometry/_plot/noise_control.py @@ -41,15 +41,18 @@ def plot_reactive_silencer( f = np.asarray(result.frequencies, dtype=np.float64) kwargs.setdefault("color", _C_PRIMARY) kwargs.setdefault("label", "Transmission loss") - ax.plot(f, np.asarray(result.transmission_loss), lw=1.8, **kwargs) + kwargs.setdefault("lw", 1.8) + ax.plot(f, np.asarray(result.transmission_loss), **kwargs) if result.insertion_loss is not None: ax.plot(f, np.asarray(result.insertion_loss), color=_C_SECONDARY, lw=1.4, ls="--", label="Insertion loss") if result.resonances is not None: - for i, fr in enumerate(np.atleast_1d(np.asarray(result.resonances))): + labeled = False + for fr in np.atleast_1d(np.asarray(result.resonances)): if f.min() <= fr <= f.max(): ax.axvline(fr, color=_C_TERTIARY, ls=":", lw=1.0, - label="Resonance" if i == 0 else None) + label="Resonance" if not labeled else None) + labeled = True ax.set_xlabel(_FREQ_LABEL) ax.set_ylabel("Loss [dB]") ax.set_title(f"Reactive silencer: {result.kind}") @@ -74,7 +77,10 @@ def plot_hvac_spectrum( is_power = result.quantity == "sound_power_level" kwargs.setdefault("color", _C_SECONDARY if is_power else _C_PRIMARY) kwargs.setdefault("label", result.label) - ax.plot(f, np.asarray(result.values), lw=1.8, marker="o", ms=3, **kwargs) + kwargs.setdefault("lw", 1.8) + kwargs.setdefault("marker", "o") + kwargs.setdefault("ms", 3) + ax.plot(f, np.asarray(result.values), **kwargs) ax.set_xlabel(_FREQ_LABEL) ax.set_ylabel( "Sound power level [dB re 1 pW]" if is_power else "Attenuation [dB]" @@ -111,8 +117,10 @@ def plot_enclosure( ls=":", marker="^", ms=3, label="Interior correction C") kwargs.setdefault("color", _C_PRIMARY) kwargs.setdefault("label", "Insertion loss (R - C)") - ax.plot(x, np.asarray(result.insertion_loss), lw=1.9, marker="o", ms=3, - **kwargs) + kwargs.setdefault("lw", 1.9) + kwargs.setdefault("marker", "o") + kwargs.setdefault("ms", 3) + ax.plot(x, np.asarray(result.insertion_loss), **kwargs) ax.set_ylabel("Level [dB]") ax.set_title("Machine enclosure insertion loss") ax.grid(True, which="both", alpha=0.3) diff --git a/src/phonometry/electroacoustics/piston.py b/src/phonometry/electroacoustics/piston.py index dd7f3212..55e92a75 100644 --- a/src/phonometry/electroacoustics/piston.py +++ b/src/phonometry/electroacoustics/piston.py @@ -94,8 +94,9 @@ def piston_resistance(x: ArrayLike) -> np.ndarray | float: raise ValueError("'x' must be non-negative and finite.") with np.errstate(invalid="ignore", divide="ignore"): out = 1.0 - 2.0 * special.j1(arr) / arr - # J1(x)/x -> 1/2 as x -> 0, so R1 -> 0; fill the removable singularity. - out = np.where(arr == 0.0, 0.0, out) + # J1(x)/x -> 1/2 as x -> 0, so R1 -> 0; the input is validated finite, so + # the 0/0 at x = 0 is the only non-finite point and is set to that limit. + out = np.where(np.isfinite(out), out, 0.0) return out[()] if out.ndim == 0 else out @@ -115,8 +116,9 @@ def piston_reactance(x: ArrayLike) -> np.ndarray | float: raise ValueError("'x' must be non-negative and finite.") with np.errstate(invalid="ignore", divide="ignore"): out = 2.0 * special.struve(1, arr) / arr - # H1(x)/x -> 0 as x -> 0 (H1 ~ 2 x^2 / 3 pi), so X1 -> 0. - out = np.where(arr == 0.0, 0.0, out) + # H1(x)/x -> 0 as x -> 0 (H1 ~ 2 x^2 / 3 pi), so X1 -> 0; the 0/0 at x = 0 + # is the only non-finite point (the input is validated finite). + out = np.where(np.isfinite(out), out, 0.0) return out[()] if out.ndim == 0 else out @@ -140,8 +142,9 @@ def piston_directivity(ka: ArrayLike, theta: ArrayLike) -> np.ndarray | float: u = ka_arr * np.sin(theta_arr) with np.errstate(invalid="ignore", divide="ignore"): out = 2.0 * special.j1(u) / u - # 2 J1(u)/u -> 1 as u -> 0 (on-axis, or ka = 0). - out = np.where(u == 0.0, 1.0, out) + # 2 J1(u)/u -> 1 as u -> 0 (on-axis, or ka = 0): the 0/0 there is the only + # non-finite point and is set to that limit. + out = np.where(np.isfinite(out), out, 1.0) return out[()] if out.ndim == 0 else out @@ -155,13 +158,12 @@ def _directivity_index(ka: NDArray[np.float64]) -> NDArray[np.float64]: """ theta = np.linspace(0.0, 0.5 * np.pi, 2001) sin_t = np.sin(theta) - di = np.empty_like(ka) - for i, kai in enumerate(ka): - d = np.asarray(piston_directivity(float(kai), theta), dtype=np.float64) - integrand = d**2 * sin_t - integral = float(np.trapezoid(integrand, theta)) - di[i] = 10.0 * np.log10(2.0 / integral) - return di + d = np.asarray( + piston_directivity(ka[:, None], theta[None, :]), dtype=np.float64 + ) + integrand = d**2 * sin_t[None, :] + integral = np.trapezoid(integrand, theta, axis=-1) + return np.asarray(10.0 * np.log10(2.0 / integral), dtype=np.float64) @dataclass(frozen=True) diff --git a/src/phonometry/noise_control/enclosures.py b/src/phonometry/noise_control/enclosures.py index aaf06591..574926e3 100644 --- a/src/phonometry/noise_control/enclosures.py +++ b/src/phonometry/noise_control/enclosures.py @@ -80,6 +80,41 @@ def plot(self, ax: "Axes | None" = None, **kwargs: Any) -> "Axes": return plot_enclosure(self, ax=ax, **kwargs) +def _resolve_frequencies( + frequencies: ArrayLike | None, +) -> NDArray[np.float64] | None: + """Validate an optional 1-D positive frequency grid (Hz).""" + if frequencies is None: + return None + freqs = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) + if freqs.ndim != 1 or freqs.size == 0: + raise ValueError("'frequencies' must be a non-empty 1-D array.") + if np.any(freqs <= 0.0) or not np.all(np.isfinite(freqs)): + raise ValueError("'frequencies' must be positive and finite.") + return freqs + + +def _resolve_panel_r( + panel_transmission_loss: ArrayLike | Callable[[NDArray[np.float64]], ArrayLike], + freqs: NDArray[np.float64] | None, +) -> NDArray[np.float64]: + """Resolve the panel transmission loss ``R`` into a validated 1-D array.""" + if callable(panel_transmission_loss): + if freqs is None: + raise ValueError( + "'frequencies' is required when 'panel_transmission_loss' " + "is a callable." + ) + r = np.atleast_1d(np.asarray(panel_transmission_loss(freqs), dtype=np.float64)) + else: + r = np.atleast_1d(np.asarray(panel_transmission_loss, dtype=np.float64)) + if r.ndim != 1 or r.size == 0: + raise ValueError("'panel_transmission_loss' must be a non-empty 1-D array.") + if not np.all(np.isfinite(r)): + raise ValueError("'panel_transmission_loss' must be finite.") + return r + + def enclosure_insertion_loss( panel_transmission_loss: ArrayLike | Callable[[NDArray[np.float64]], ArrayLike], external_area: float, @@ -110,29 +145,12 @@ def enclosure_insertion_loss( s_e = require_positive(external_area, "external_area") s_i = require_positive(internal_area, "internal_area") - freqs: NDArray[np.float64] | None = None - if frequencies is not None: - freqs = np.atleast_1d(np.asarray(frequencies, dtype=np.float64)) - if freqs.ndim != 1 or freqs.size == 0: - raise ValueError("'frequencies' must be a non-empty 1-D array.") - if np.any(freqs <= 0.0) or not np.all(np.isfinite(freqs)): - raise ValueError("'frequencies' must be positive and finite.") - - if callable(panel_transmission_loss): - if freqs is None: - raise ValueError( - "'frequencies' is required when 'panel_transmission_loss' " - "is a callable." - ) - r = np.atleast_1d(np.asarray(panel_transmission_loss(freqs), dtype=np.float64)) - else: - r = np.atleast_1d(np.asarray(panel_transmission_loss, dtype=np.float64)) - if r.ndim != 1 or r.size == 0: - raise ValueError("'panel_transmission_loss' must be a non-empty 1-D array.") - if not np.all(np.isfinite(r)): - raise ValueError("'panel_transmission_loss' must be finite.") + freqs = _resolve_frequencies(frequencies) + r = _resolve_panel_r(panel_transmission_loss, freqs) alpha = np.asarray(internal_absorption, dtype=np.float64) + if alpha.ndim > 1: + raise ValueError("'internal_absorption' must be a scalar or a 1-D array.") if np.any(alpha <= 0.0) or np.any(alpha >= 1.0) or not np.all(np.isfinite(alpha)): raise ValueError("'internal_absorption' must lie strictly in (0, 1).") diff --git a/src/phonometry/noise_control/silencers.py b/src/phonometry/noise_control/silencers.py index 480dae7d..5244c679 100644 --- a/src/phonometry/noise_control/silencers.py +++ b/src/phonometry/noise_control/silencers.py @@ -149,6 +149,11 @@ def cascade(*matrices: _Complex) -> _Complex: """ if not matrices: raise ValueError("cascade() needs at least one matrix.") + n = matrices[0].shape[0] + if any(m.shape[0] != n for m in matrices[1:]): + raise ValueError( + "cascade() matrices must share the same frequency grid (n_freq)." + ) total = matrices[0] for m in matrices[1:]: total = np.matmul(total, m) @@ -544,6 +549,11 @@ def extended_tube_chamber( lb = require_non_negative(outlet_extension, "outlet_extension") if s_exp <= s_duct: raise ValueError("'chamber_area' must exceed 'pipe_area'.") + if la + lb > length: + raise ValueError( + "the inlet and outlet extensions cannot together exceed the " + "chamber length (they would overlap inside the chamber)." + ) annulus = s_exp - s_duct elements = [] diff --git a/tests/noise_control/test_fdtd_crosscheck.py b/tests/noise_control/test_fdtd_crosscheck.py index 35bcf85e..aeacb03c 100644 --- a/tests/noise_control/test_fdtd_crosscheck.py +++ b/tests/noise_control/test_fdtd_crosscheck.py @@ -79,8 +79,11 @@ def test_fdtd_matches_expansion_chamber_peak_tl() -> None: def test_fdtd_run_is_deterministic() -> None: + # Two independent runs of the (RNG-free) solver must be bit-identical. f = _C / (4.0 * _L) - assert _transmitted_rms(f) == _transmitted_rms(f) + first = _transmitted_rms(f) + second = _transmitted_rms(f) + assert first == second def test_four_pole_peak_frequency_and_trough() -> None: