1"""The losses derived from physical constants, and the per-mode rates."""
3import numpy as np
5from cymsim import Plate, compute_modes, make_grid
6from cymsim.constants import RHO_AIR, ROOM_TEMPERATURE, SPEED_OF_SOUND
7from cymsim.linear.damping import bending_energy_per_node, radiation_rate_direct, radiation_rates, thermoelastic_rates
8from cymsim.loss import (
9 LN1000,
10 average_radiation_efficiency,
11 coincidence_frequency,
12 estimated_rate,
13 fit_line,
14 fit_sigma,
15 fitted_rate,
16 material_rate,
17 outline_perimeter,
18 plate_geometry,
19 thermoelastic_loss_factor,
20 thermoelastic_relaxation,
21)
23from conftest import BRONZE, STEEL_LOSSY
26def test_zener_peak():
27 h = 1e-3
28 delta, tau = thermoelastic_relaxation(BRONZE, h, ROOM_TEMPERATURE)
29 th = BRONZE.thermal
30 assert abs(tau - h * h * BRONZE.rho * th.cp / (np.pi**2 * th.k)) < 1e-15
31 assert abs(delta - BRONZE.E * th.alpha**2 * ROOM_TEMPERATURE / (BRONZE.rho * th.cp) * (1 + BRONZE.nu) / (1 - BRONZE.nu)) < 1e-15
32 peak = thermoelastic_loss_factor(BRONZE, h, 1 / tau, ROOM_TEMPERATURE)
33 assert abs(peak - delta / 2) < 1e-15
34 assert thermoelastic_loss_factor(BRONZE, h, 0.1 / tau, ROOM_TEMPERATURE) < peak > thermoelastic_loss_factor(BRONZE, h, 10 / tau, ROOM_TEMPERATURE)
35 assert 1 / (2 * np.pi * tau) < 100 and 1e-4 < delta / 2 < 1e-2
36 assert abs(thermoelastic_relaxation(BRONZE, 2 * h, ROOM_TEMPERATURE)[1] / tau - 4) < 1e-12
37 from cymsim import Material
39 assert thermoelastic_loss_factor(Material(110e9, 0.34, 8800), h, 1000, ROOM_TEMPERATURE) == 0
42def test_outline_perimeter():
43 rect = make_grid(0.3, 0.2, 0.01)
44 assert abs(outline_perimeter(rect) / (2 * (0.3 + 0.2)) - 1) < 0.02
45 R = 0.2
46 disc = make_grid(2 * R, 2 * R, 0.008, lambda x, y, r, t: r <= R)
47 p = outline_perimeter(disc)
48 assert 2 * np.pi * R < p < 1.08 * 2 * np.pi * R
49 ring = make_grid(2 * R, 2 * R, 0.008, lambda x, y, r, t: (r <= R) & (r >= 0.05))
50 assert abs(outline_perimeter(ring) / (2 * np.pi * (R + 0.05)) - 1) < 0.1
53def test_average_radiation_efficiency_limits():
54 fc = coincidence_frequency(BRONZE, 1e-3)
55 assert 15000 < fc < 20000
56 A, U, a, b = 0.13, 1.28, 0.4, 0.4
57 s1 = average_radiation_efficiency(100, fc, A, U, a, b, False)
58 s2 = average_radiation_efficiency(200, fc, A, U, a, b, False)
59 assert 0 < s1 < 0.1 and abs(s2 / s1 - 2) < 0.1
60 assert average_radiation_efficiency(fc / 2, fc, A, U, a, b, False) > s2
61 assert abs(average_radiation_efficiency(4 * fc, fc, A, U, a, b) - 1 / np.sqrt(0.75)) < 1e-9
62 assert average_radiation_efficiency(fc * 1.0001, fc, A, U, a, b) <= (0.5 - 0.15) * np.sqrt(a / (SPEED_OF_SOUND / fc)) + 1e-9
63 f1, f2 = average_radiation_efficiency(100, fc, A, U, a, b), average_radiation_efficiency(200, fc, A, U, a, b)
64 assert abs(f1 / s1 - (100 / fc) ** 2) < 1e-12 and abs(f2 / f1 - 8) < 0.4
65 # Vectorized over frequency.
66 arr = average_radiation_efficiency(np.array([100.0, 200.0, 4 * fc]), fc, A, U, a, b)
67 assert np.allclose(arr[:2], [f1, f2]) and abs(arr[2] - 1 / np.sqrt(0.75)) < 1e-9
70def test_two_parameter_fit():
71 x = np.array([1, 2, 3, 5, 8.0])
72 fit = fit_line(x, 0.7 + 2e-4 * x, 1.1)
73 assert abs(fit.sigma0 - 0.7) < 1e-9 and abs(fit.sigma1 - 2e-4) < 1e-12
74 assert fit_line(x, 3 - x, 1.1).sigma1 == 0
75 pl = Plate.uniform(make_grid(0.4, 0.4, 0.02), BRONZE, 1e-3)
76 geom = plate_geometry(pl)
77 f2 = fit_sigma(pl, geom)
78 assert f2.sigma0 >= 0 and f2.sigma1 >= 0
79 e100, e5k = estimated_rate(pl, geom, 100), estimated_rate(pl, geom, 5000)
80 assert e5k.total > e100.total
81 assert abs(e100.material - material_rate(BRONZE.loss_factor, 2 * np.pi * 100)) < 1e-12
82 assert e100.radiation > 0 and e100.thermoelastic > 0
83 assert 0.3 * e100.total < fitted_rate(f2, 1000) < 3 * e5k.total
86def small_modes():
87 pl = Plate.uniform(make_grid(0.2, 0.15, 0.0125), STEEL_LOSSY, 1e-3)
88 return pl, compute_modes(pl, 4000)
91def test_piston_radiation_rate_equals_the_closed_form():
92 pl, ms = small_modes()
93 one = ms.select([0])
94 M, S = pl.mass, pl.area
95 one.w[0] = np.where(pl.grid.mask, 1 / np.sqrt(M), 0.0)
96 one.freqs[0] = 20.0
97 omega = 2 * np.pi * 20
98 expected = RHO_AIR * omega**2 * S**2 / (2 * np.pi * SPEED_OF_SOUND * M)
99 assert abs(radiation_rates(one)[0] / expected - 1) < 2e-3
100 assert abs(radiation_rate_direct(pl, one.w[0], 20.0) / expected - 1) < 2e-3
103def test_fft_convolution_equals_the_direct_double_sum():
104 pl, ms = small_modes()
105 rates = radiation_rates(ms)
106 for m in range(0, ms.count, max(1, ms.count // 12)):
107 direct = radiation_rate_direct(pl, ms.w[m], ms.freqs[m])
108 assert abs(rates[m] - direct) < 1e-9 * max(1, direct)
109 assert (rates >= 0).all() and rates[-1] > rates[0]
112def test_wallace_low_frequency_efficiencies():
113 # sin(m pi x / a) sin(n pi y / b) on a baffled rectangle at k a << 1 (Wallace 1972).
114 a, b, h, f = 0.3, 0.45, 1e-3, 30.0
115 c = SPEED_OF_SOUND
116 k = 2 * np.pi * f / c
117 g = make_grid(a, b, 0.01)
118 pl = Plate.uniform(g, STEEL_LOSSY, h)
119 rho_h = STEEL_LOSSY.rho * h
120 xs = (g.X - g.x[0]) / a
121 ys = (g.Y - g.y[0]) / b
122 for m, n in [(1, 1), (1, 2), (2, 1), (2, 2), (3, 3)]:
123 v = np.sin(m * np.pi * xs) * np.sin(n * np.pi * ys)
124 v /= np.sqrt((rho_h * g.node_area * v * v).sum())
125 efficiency = radiation_rate_direct(pl, v, f) * rho_h / (RHO_AIR * c)
126 odd = lambda q: q % 2 == 1 # noqa: E731
127 if odd(m) and odd(n):
128 wallace = 32 * k * k * a * b / (np.pi**5 * m * m * n * n)
129 elif not odd(m) and not odd(n):
130 wallace = 2 * k**6 * a**3 * b**3 / (15 * np.pi**5 * m * m * n * n)
131 elif odd(m):
132 wallace = 8 * k**4 * a * b**3 / (3 * np.pi**5 * m * m * n * n)
133 else:
134 wallace = 8 * k**4 * a**3 * b / (3 * np.pi**5 * m * m * n * n)
135 assert abs(efficiency / wallace - 1) < 0.04, (m, n)
138def test_thermoelastic_rate_of_a_uniform_plate_and_the_total():
139 pl, ms = small_modes()
140 rates = thermoelastic_rates(ms)
141 omega = 2 * np.pi * ms.freqs
142 expected = thermoelastic_loss_factor(STEEL_LOSSY, 1e-3, omega, ROOM_TEMPERATURE) * omega / 2
143 assert np.abs(rates / expected - 1).max() < 1e-4
144 U = bending_energy_per_node(pl, ms.w[3]).sum()
145 assert abs(U / (omega[3] ** 2 / 2) - 1) < 1e-4
146 assert np.allclose(ms.sigma, ms.sigma_material + ms.sigma_thermal + ms.sigma_radiation)
147 assert (ms.sigma_material > 0).all() and (ms.sigma_radiation > 0).all() and (ms.sigma_thermal > 0).all()
148 assert ((LN1000 / ms.sigma > 0.1) & (LN1000 / ms.sigma < 1000)).all()
149 by = ms.t60_by_mechanism()
150 assert set(by) == {"material", "thermal", "radiation"}
151 assert isinstance(ms.table(5), str)