[FEATURE] add pytta.absorbers with layered absorbers

docs/pyttadoc/absorbers.rst

@@ -0,0 +1,7 @@
+.. currentmodule:: pytta
+
+Absorbers
+=========
+
+.. automodule:: pytta.absorbers
+   :members:

docs/pyttadoc/index.rst

@@ -57,11 +57,15 @@ Modules:
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.absorbers:
+        surface impedance and absorption of layered absorbers;
+
 .. toctree::
    classes/index
    utils/index
    apps/index
    functions
    generate
    rooms
+   absorbers
    properties

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

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

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()