From 73e88a51e9fd1015adaa4e0e18641c3d8cc9bd67 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 01:23:58 +0100 Subject: [PATCH] [FEATURE] add pytta.absorbers with layered absorbers pytta.absorbers computes the surface impedance and absorption of layered absorbers by the transfer matrix method, as the Impedance Calculator of the ITA-Toolbox: layers of air, porous material (Delany-Bazley, Miki with the constants of the paper, or Komatsu), microperforated panels (Maa) and limp masses on a rigid backing, at normal or oblique incidence, and the random-incidence absorption by Paris' formula. Tests check that a rigid-backed air layer is purely reactive, that the transfer matrix gives Zc coth(Gamma d) for a porous layer and Zc for a very thick one, that the models tend to air, the Poiseuille limit of the microperforated panel, and Paris' formula against its closed form for a real impedance. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/absorbers.rst | 7 ++ docs/pyttadoc/index.rst | 4 + pytta/__init__.py | 5 + pytta/absorbers.py | 191 ++++++++++++++++++++++++++++++++++++ tests/test_absorbers.py | 96 ++++++++++++++++++ 5 files changed, 303 insertions(+) create mode 100644 docs/pyttadoc/absorbers.rst create mode 100644 pytta/absorbers.py create mode 100644 tests/test_absorbers.py diff --git a/docs/pyttadoc/absorbers.rst b/docs/pyttadoc/absorbers.rst new file mode 100644 index 0000000..25a559b --- /dev/null +++ b/docs/pyttadoc/absorbers.rst @@ -0,0 +1,7 @@ +.. currentmodule:: pytta + +Absorbers +========= + +.. automodule:: pytta.absorbers + :members: diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..4c74601 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -57,6 +57,9 @@ Modules: * pytta.iso3741: calculations according to the standard; + * pytta.absorbers: + surface impedance and absorption of layered absorbers; + .. toctree:: classes/index utils/index @@ -64,4 +67,5 @@ Modules: functions generate rooms + absorbers properties diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..bc91934 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -72,6 +72,9 @@ * pytta.iso3741: calculations according to the standard; + * pytta.absorbers: + surface impedance and absorption of layered absorbers; + You can find usage examples in the toolbox's examples folder. @@ -108,6 +111,7 @@ from . import utils from . import rooms from . import iso3741 +from . import absorbers from .functions import read_wav, write_wav, merge, split, \ list_devices, print_devices, get_device_from_user,\ @@ -130,6 +134,7 @@ 'generate', 'utils', 'iso3741', + 'absorbers', 'rooms', # Functions diff --git a/pytta/absorbers.py b/pytta/absorbers.py new file mode 100644 index 0000000..16fa376 --- /dev/null +++ b/pytta/absorbers.py @@ -0,0 +1,191 @@ +""" +PyTTa absorbers: +---------------- + +Surface impedance and absorption of layered absorbers by the transfer +matrix method, as the Impedance Calculator of the ITA-Toolbox +(applications/Analytic/Impedance_Calculator): layers of air, porous +material (Delany-Bazley, Miki or Komatsu), microperforated panels (Maa) and +limp masses on a rigid backing, at normal or oblique incidence, and the +random-incidence absorption of a locally reacting surface. + +Available functions: + + >>> Zc, k = pytta.absorbers.porous(f, sigma, model='miki') + >>> Zs = pytta.absorbers.surface_impedance(f, layers) + >>> alpha = pytta.absorbers.absorption(Zs) + >>> alpha_d = pytta.absorbers.random_incidence_absorption(f, layers) + +Time dependence exp(j w t). Impedances in Pa s/m, flow resistivity sigma in +Pa s/m^2, lengths in m. + +References +---------- + +Delany, M. E., & Bazley, E. N. (1970). Acoustical properties of fibrous +absorbent materials. Applied Acoustics, 3(2), 105-116. + +Miki, Y. (1990). Acoustical properties of porous materials. Modifications of +Delany-Bazley models. Journal of the Acoustical Society of Japan (E), 11(1), +19-24. + +Komatsu, T. (2008). Improvement of the Delany-Bazley and Miki models for +fibrous sound-absorbing materials. Acoustical Science and Technology, 29(2), +121-129. + +Maa, D.-Y. (1998). Potential of microperforated panel absorber. The Journal +of the Acoustical Society of America, 104(5), 2861-2866. + +Cox, T. J., & D'Antonio, P. (2009). Acoustic Absorbers and Diffusers (2nd +ed.). Taylor & Francis. +""" + +import numpy as np + +__all__ = ['porous', 'mpp_impedance', 'surface_impedance', 'absorption', + 'random_incidence_absorption'] + +_RHO, _C, _ETA = 1.204, 343.2, 1.81e-5 + + +def porous(f, sigma, model='miki', rho=_RHO, c=_C): + """ + Characteristic impedance Zc (Pa s/m) and wavenumber k (1/m, with a + negative imaginary part) of a porous material of flow resistivity sigma + (Pa s/m^2), by the empirical models of Delany and Bazley (1970), Miki + (1990) or Komatsu (2008): + + * 'delany_bazley': X = rho f / sigma, Zc = rho c (1 + 0.0571 X^-0.754 + - j 0.087 X^-0.732), k = w/c (1 + 0.0978 X^-0.700 - j 0.189 X^-0.595); + * 'miki': X = 1000 f / sigma, Zc = rho c (1 + 5.50 X^-0.632 - j 8.43 + X^-0.632), k = w/c (1 + 7.81 X^-0.618 - j 11.41 X^-0.618), the + constants of Miki's paper, Eqs. (27) and (28) (often rounded in SI + units to 0.0699, 0.107, 0.109 and 0.160 with f / sigma); + * 'komatsu': u = 2 - lg(f / sigma), Zc = rho c (1 + 0.00027 u^6.2 - + j 0.0047 u^4.1), k = w/c (1 + 0.0004 u^6.2 - j 0.0069 u^4.1), for + f / sigma up to 100 (u >= 0), where it reaches air; above, nan. + """ + f = np.asarray(f, dtype=float) + k0 = 2 * np.pi * f / c + if model == 'delany_bazley': + X = rho * f / sigma + Z = 1 + 0.0571 * X ** -0.754 - 1j * 0.087 * X ** -0.732 + k = 1 + 0.0978 * X ** -0.700 - 1j * 0.189 * X ** -0.595 + elif model == 'miki': + X = 1000 * f / sigma # f / sigma in Hz per kPa s/m^2 + Z = 1 + 5.50 * X ** -0.632 - 1j * 8.43 * X ** -0.632 + k = 1 + 7.81 * X ** -0.618 - 1j * 11.41 * X ** -0.618 + elif model == 'komatsu': + u = 2 - np.log10(f / sigma) + u = np.where(u >= 0, u, np.nan) + Z = 1 + 0.00027 * u ** 6.2 - 1j * 0.0047 * u ** 4.1 + k = 1 + 0.0004 * u ** 6.2 - 1j * 0.0069 * u ** 4.1 + else: + raise ValueError("model must be 'delany_bazley', 'miki' or 'komatsu'.") + return rho * c * Z, k0 * k + + +def mpp_impedance(f, diameter, thickness, porosity, rho=_RHO, c=_C, eta=_ETA): + """ + Series impedance (Pa s/m) of a microperforated panel, Maa (1998): + + Z / (rho c) = 32 eta t k_r / (p rho c d^2) + j w t k_m / (p c), + k_r = sqrt(1 + x^2 / 32) + sqrt(2) x d / (32 t), + k_m = 1 + 1 / sqrt(9 + x^2 / 2) + 0.85 d / t, x = d sqrt(w rho / (4 eta)), + + with hole diameter d, panel thickness t and perforation ratio p. + """ + f = np.asarray(f, dtype=float) + w = 2 * np.pi * f + d, t, p = diameter, thickness, porosity + x = d * np.sqrt(w * rho / (4 * eta)) + kr = np.sqrt(1 + x ** 2 / 32) + np.sqrt(2) / 32 * x * d / t + km = 1 + 1 / np.sqrt(9 + x ** 2 / 2) + 0.85 * d / t + return rho * c * (32 * eta * t / (p * rho * c * d ** 2) * kr + 1j * w * t / (p * c) * km) + + +def _layer(f, layer, sin0, rho, c): + """Transfer matrix (2, 2, nf) of one layer along the normal, for an + incidence with sin(theta0) = sin0 in air (laterally reacting fluids).""" + kind = layer[0] + nf = np.size(f) + if kind in ('air', 'porous'): + if kind == 'air': + Zc = np.full(nf, rho * c, dtype=complex) + k = 2 * np.pi * np.asarray(f, dtype=float) / c + 0j + d = layer[1] + else: + _, model, sigma, d = layer + Zc, k = porous(f, sigma, model, rho, c) + k0 = 2 * np.pi * np.asarray(f, dtype=float) / c + kx = np.sqrt(k ** 2 - (k0 * sin0) ** 2) # wavenumber along the normal + Zn = Zc * k / kx + cs, sn = np.cos(kx * d), np.sin(kx * d) + return np.array([[cs, 1j * Zn * sn], [1j * sn / Zn, cs]]) + if kind == 'mpp': + _, dia, t, p = layer + Z = mpp_impedance(f, dia, t, p, rho, c) + elif kind == 'mass': + Z = 1j * 2 * np.pi * np.asarray(f, dtype=float) * layer[1] + else: + raise ValueError(f'unknown layer {kind!r}.') + one, zero = np.ones(nf, dtype=complex), np.zeros(nf, dtype=complex) + return np.array([[one, Z], [zero, one]]) + + +def surface_impedance(f, layers, angle=0., rho=_RHO, c=_C): + """ + Surface impedance (Pa s/m) of layers on a rigid backing, by the transfer + matrix method, at an angle of incidence (degrees) in air. + + Parameters + ---------- + f : array_like + Frequencies (Hz). + layers : list of tuples + From the surface down to the backing: + ('air', thickness), ('porous', model, sigma, thickness), + ('mpp', hole diameter, thickness, perforation ratio) or + ('mass', mass per unit area in kg/m^2). + angle : float + Angle of incidence in degrees. Air and porous layers react + laterally (the trace wavenumber is kept); the panels react locally. + """ + f = np.asarray(f, dtype=float) + sin0 = np.sin(np.radians(angle)) + T = np.broadcast_to(np.eye(2, dtype=complex)[:, :, None], (2, 2, f.size)).copy() + for layer in layers: + T = np.einsum('ijn,jkn->ikn', T, _layer(f, layer, sin0, rho, c)) + with np.errstate(divide='ignore'): + return T[0, 0] / T[1, 0] # rigid backing: velocity 0 at the bottom + + +def absorption(Zs, angle=0., rho=_RHO, c=_C): + """ + Absorption coefficient 1 - |R|^2 of a surface impedance at an angle of + incidence (degrees), R = (Zs cos(theta) - rho c) / (Zs cos(theta) + rho c). + """ + Zs = np.asarray(Zs, dtype=complex) + a = Zs * np.cos(np.radians(angle)) + R = (a - rho * c) / (a + rho * c) + return 1 - np.abs(R) ** 2 + + +def random_incidence_absorption(f, layers=None, Zs=None, n_angles=180, + rho=_RHO, c=_C): + """ + Random-incidence absorption by Paris' formula, the integral of + alpha(theta) sin(2 theta) over 0 to 90 degrees. + + Give layers, whose surface impedance is computed at each angle, or the + impedance Zs of a locally reacting surface, the same at all angles. + """ + edges = np.linspace(0, 90, n_angles + 1) + th = 0.5 * (edges[1:] + edges[:-1]) + dth = np.radians(edges[1] - edges[0]) + out = 0 + for a in th: + Z = (surface_impedance(f, layers, a, rho, c) if Zs is None + else np.broadcast_to(np.asarray(Zs, dtype=complex), np.shape(f))) + out = out + absorption(Z, a, rho, c) * np.sin(2 * np.radians(a)) * dth + return out diff --git a/tests/test_absorbers.py b/tests/test_absorbers.py new file mode 100644 index 0000000..dfd97c5 --- /dev/null +++ b/tests/test_absorbers.py @@ -0,0 +1,96 @@ +""" +Tests for pytta.absorbers, layered absorbers by the transfer matrix method +after the Impedance Calculator of the ITA-Toolbox. Each test pins one rule +of the models or of the method. +""" +import unittest + +import numpy as np + +from pytta import absorbers as A + +RHO, C = 1.204, 343.2 +F = np.array([63, 125, 250, 500, 1000, 2000, 4000, 8000.]) + + +class TestAbsorbers(unittest.TestCase): + + def test_air_layer_on_rigid_backing_does_not_absorb(self): + # Zs = -j rho c cot(k d), purely reactive, so alpha = 0. + d = 0.1 + Zs = A.surface_impedance(F, [('air', d)]) + k = 2 * np.pi * F / C + e = np.abs(Zs - (-1j * RHO * C / np.tan(k * d))).max() / (RHO * C) + print(f'\n air layer: Zs against -j rho c cot(kd) {e:.1e}, max alpha ' + f'{np.abs(A.absorption(Zs)).max():.1e}') + self.assertLess(e, 1e-9) + self.assertLess(np.abs(A.absorption(Zs)).max(), 1e-12) + + def test_porous_layer_by_the_transfer_matrix(self): + # One porous layer on a rigid backing: Zs = -j Zc cot(k d) (Aretz, + # eq. 7.14, with Gamma = j k); a very thick layer tends to Zc. + Zc, k = A.porous(F, 15000, 'miki') + Zs = A.surface_impedance(F, [('porous', 'miki', 15000, 0.08)]) + thick = A.surface_impedance(F[3:], [('porous', 'miki', 15000, 5.0)]) + e = np.abs(Zs - (-1j * Zc / np.tan(k * 0.08))).max() / np.abs(Zc).max() + e_thick = np.abs(thick / Zc[3:] - 1).max() + print(f'\n porous layer {e:.1e}; 5 m layer against Zc {e_thick:.1e}') + self.assertLess(e, 1e-12) + self.assertLess(e_thick, 1e-6) + + def test_models_tend_to_air_and_miki_constants(self): + # At f / sigma large Delany-Bazley and Miki tend to air, and Komatsu + # reaches it at f / sigma = 100, the end of its range; Miki's constants in + # the model are those of his paper, Eqs. (27) and (28), and its usual + # SI form (0.0699, 0.107, 0.109, 0.160 with f / sigma) agrees to its + # rounding. + f = np.array([1e6]) + for model in ('delany_bazley', 'miki'): + Zc, k = A.porous(f, 100, model) + self.assertAlmostEqual(abs(Zc[0] / (RHO * C) - 1), 0, delta=2e-3) + self.assertAlmostEqual(abs(k[0] / (2 * np.pi * f[0] / C) - 1), 0, delta=2e-3) + Zc, k = A.porous(np.array([1e4]), 100, 'komatsu') + self.assertAlmostEqual(abs(Zc[0] / (RHO * C) - 1), 0, places=12) + self.assertTrue(np.isnan(A.porous(np.array([1e5]), 100, 'komatsu')[0][0])) + f, sigma = np.array([500.]), 10000. + Zc, k = A.porous(f, sigma, 'miki') + X = f / sigma + Zp = 1 + 0.0699 * X ** -0.632 - 1j * 0.107 * X ** -0.632 + kp = 1 + 0.109 * X ** -0.618 - 1j * 0.160 * X ** -0.618 + e = max(abs(Zc[0] / (RHO * C) - Zp[0]), abs(k[0] / (2 * np.pi * 500 / C) - kp[0])) + print(f'\n Miki against its rounded SI form {e:.1e}') + self.assertLess(e, 3e-3) + + def test_mpp_tends_to_poiseuille(self): + # Maa: as w -> 0, k_r -> 1 (its second term is linear in x, hence a + # very low frequency) and k_m -> 4/3 + 0.85 d / t, so the + # resistance is the Poiseuille resistance 32 eta t / (p d^2) and the + # mass rho t (4/3 + 0.85 d / t) / p. + d, t, p, eta = 0.4e-3, 0.5e-3, 0.01, 1.81e-5 + f = np.array([1e-6]) + Z = A.mpp_impedance(f, d, t, p)[0] + R0 = 32 * eta * t / (p * d ** 2) + M0 = RHO * t * (4 / 3 + 0.85 * d / t) / p + print(f'\n MPP at 1 uHz: R {Z.real:.4f} against {R0:.4f}, ' + f'm {Z.imag / (2 * np.pi * f[0]):.6f} against {M0:.6f}') + self.assertAlmostEqual(Z.real / R0, 1, delta=1e-5) + self.assertAlmostEqual(Z.imag / (2 * np.pi * f[0]) / M0, 1, delta=1e-4) + + def test_random_incidence_for_a_real_impedance(self): + # For a locally reacting surface of real normalised impedance z, + # Paris' integral is 8/z (1 + 1/(1 + z) - 2/z ln(1 + z)). + for z in (0.5, 1.0, 3.0, 10.0): + a = A.random_incidence_absorption(np.array([1000.]), Zs=z * RHO * C, n_angles=2000)[0] + ref = 8 / z * (1 + 1 / (1 + z) - 2 / z * np.log(1 + z)) + print(f'\n z = {z}: {a:.6f} against {ref:.6f}', end='') + self.assertAlmostEqual(a, ref, delta=1e-5) + + def test_rejects_unknown_input(self): + with self.assertRaises(ValueError): + A.porous(F, 1e4, 'biot') + with self.assertRaises(ValueError): + A.surface_impedance(F, [('foam', 0.1)]) + + +if __name__ == '__main__': + unittest.main()