1"""The losses of the plate, derived from physical constants rather than prescribed.
3Three mechanisms drain the energy of a vibrating plate in this model. Internal
4friction of the metal is a hysteretic loss with the material's loss factor eta:
5a mode of angular frequency omega decays at the rate eta omega / 2. Thermoelastic
6damping (Zener) comes from the heat that flows across the thickness between the
7compressed and the stretched face of a bent plate; it follows from the thermal
8constants of the metal and the local thickness, with a peak at the frequency set
9by the thermal diffusion time across the thickness. Radiation into the air is
10the energy the plate loses as sound; the linear model computes it exactly for
11each mode from the mode's shape (``linear.damping``), and this module gives the
12analytic average of Leppington for a plate of given area and perimeter, which
13the full model's two-parameter loss operator is fitted to.
14"""
16from __future__ import annotations
18import math
19from dataclasses import dataclass
21import numpy as np
23from .constants import LN1000, RHO_AIR, SPEED_OF_SOUND
24from .grid import Grid
25from .material import Material
26from .plate import Plate
28__all__ = [
29 "LN1000",
30 "material_rate",
31 "thermoelastic_relaxation",
32 "thermoelastic_loss_factor",
33 "outline_perimeter",
34 "coincidence_frequency",
35 "average_radiation_efficiency",
36 "radiation_rate_from_efficiency",
37 "PlateGeometry",
38 "plate_geometry",
39 "LossEstimate",
40 "estimated_rate",
41 "in_plane_rate",
42 "SigmaFit",
43 "fit_sigma",
44 "fit_line",
45 "fitted_rate",
46 "t60",
47]
50def t60(sigma):
51 """Decay time [s] of a rate sigma [1/s]: the time in which a partial falls by 60 dB."""
52 with np.errstate(divide="ignore"):
53 return LN1000 / np.asarray(sigma, dtype=float)
56# ---------- Internal friction ----------
59def material_rate(eta: float, omega):
60 """Decay rate of a mode of angular frequency omega from a hysteretic loss
61 factor eta: the energy lost per cycle is 2 pi eta times the stored energy."""
62 return eta * omega / 2
65# ---------- Thermoelastic damping ----------
68def thermoelastic_relaxation(mat: Material, h, T: float):
69 """Zener's relaxation strength ``delta`` and time ``tau`` [s] of a plate of
70 thickness ``h``. The strength carries the factor (1 + nu) / (1 - nu) of
71 the biaxial stress state of a bent plate (Zener's beam value is
72 E alpha^2 T / (rho cp)); the time is the thermal diffusion time across the
73 thickness, h^2 / (pi^2 chi) with the diffusivity chi = k / (rho cp)."""
74 th = mat.thermal
75 if th is None:
76 return 0.0, np.ones_like(np.asarray(h, dtype=float))
77 delta = (mat.E * th.alpha**2 * T) / (mat.rho * th.cp) * ((1 + mat.nu) / (1 - mat.nu))
78 tau = np.asarray(h, dtype=float) ** 2 * mat.rho * th.cp / (math.pi**2 * th.k)
79 return delta, tau
82def thermoelastic_loss_factor(mat: Material, h, omega, T: float):
83 """Thermoelastic loss factor of bending at angular frequency ``omega`` in a
84 plate of thickness ``h``: delta x / (1 + x^2) with x = omega tau, largest
85 (delta / 2) at x = 1. ``h`` and ``omega`` broadcast."""
86 delta, tau = thermoelastic_relaxation(mat, h, T)
87 x = np.asarray(omega, dtype=float) * tau
88 return delta * x / (1 + x * x)
91# ---------- Geometry of the outline ----------
93_DIAG = 0
94_HORIZ = 1
95_VERT = 2
96# Segments of the marching-squares contour by code, as counts of (diagonal, horizontal, vertical) pieces.
97_SEGMENTS = {1: (1, 0, 0), 14: (1, 0, 0), 2: (1, 0, 0), 13: (1, 0, 0), 3: (0, 1, 0), 12: (0, 1, 0), 4: (1, 0, 0), 11: (1, 0, 0), 6: (0, 0, 1), 9: (0, 0, 1), 7: (1, 0, 0), 8: (1, 0, 0), 5: (2, 0, 0), 10: (2, 0, 0)}
100def outline_perimeter(g: Grid) -> float:
101 """Length [m] of the outline of the plate, from the marching-squares contour
102 of the cell mask at level one half (cells outside the grid count as outside).
103 A staircase count would overstate the perimeter of a disc by 4 / pi; the
104 contour, whose segments run in eight directions, overstates it by a few
105 percent, which is well within the accuracy of the average radiation
106 efficiency it feeds."""
107 c = np.pad(g.cell_mask, 1).astype(np.int64)
108 a, b, cc, d = c[:-1, :-1], c[:-1, 1:], c[1:, 1:], c[1:, :-1]
109 code = a | (b << 1) | (cc << 2) | (d << 3)
110 counts = np.bincount(code.ravel(), minlength=16)
111 diag = math.hypot(0.5 * g.hx, 0.5 * g.hy)
112 length = 0.0
113 for k, (nd, nh, nv) in _SEGMENTS.items():
114 length += counts[k] * (nd * diag + nh * g.hx + nv * g.hy)
115 return float(length)
118# ---------- Radiation ----------
121def coincidence_frequency(mat: Material, h) -> float:
122 """Frequency at which the bending wavelength equals the wavelength of sound, c^2 / (2 pi kappa) [Hz]."""
123 return mat.coincidence_frequency(h)
126def average_radiation_efficiency(f, fc: float, area: float, perimeter: float, a: float, b: float, free_edges: bool = True):
127 """The radiation efficiency of a baffled plate averaged over its modes near
128 the frequency ``f``. The simply supported case is that of Leppington,
129 Broadbent and Heron (1982): below coincidence the edges radiate and the
130 efficiency grows about linearly with frequency, at coincidence it peaks
131 at a value set by the size of the plate in coincidence wavelengths, and
132 above it tends to one. A plate with free edges radiates much less below
133 coincidence (its modes have no monopole or dipole moment, and at a free
134 edge the bending moment and the shear vanish); the per-mode computation
135 on the presets follows the simply supported curve times (f / fc)^2 within
136 a factor of two from 100 Hz to 8 kHz, and that is what ``free_edges``
137 applies. ``a`` and ``b`` are the shorter and the longer side."""
138 f = np.asarray(f, dtype=float)
139 c = SPEED_OF_SOUND
140 lambda_c = c / fc
141 cap = max(1.0, (0.5 - 0.15 * (a / b)) * math.sqrt(a / lambda_c))
142 above = f >= fc
143 fb = np.where(above, 0.5 * fc, f) # safe stand-in where the other branch applies
144 fa = np.where(above, f, 2 * fc)
145 with np.errstate(divide="ignore", invalid="ignore"):
146 s_above = np.minimum(cap, 1 / np.sqrt(np.maximum(1e-12, 1 - fc / fa)))
147 mu = np.sqrt(fc / fb)
148 s_below = (perimeter * c) / (4 * math.pi**2 * area * fc) / np.sqrt(mu * mu - 1) * (np.log((mu + 1) / (mu - 1)) + 2 * mu / (mu * mu - 1))
149 s_below = np.minimum(cap, s_below) * ((fb / fc) ** 2 if free_edges else 1.0)
150 out = np.where(above, s_above, s_below)
151 return float(out) if out.ndim == 0 else out
154def radiation_rate_from_efficiency(efficiency, mat: Material, h):
155 """Decay rate from a radiation efficiency: the two faces of a baffled plate
156 of surface mass rho h radiate the power 2 rho_air c sigma_rad <v^2> A, so
157 sigma = rho_air c sigma_rad / (rho h)."""
158 return RHO_AIR * SPEED_OF_SOUND * efficiency / (mat.rho * h)
161# ---------- The estimate for a whole plate ----------
165class PlateGeometry:
166 """What the loss estimate needs of a plate."""
168 area: float
169 perimeter: float
170 h_mean: float
171 kappa: float
172 fc: float
173 a: float
174 """Shorter side of the bounding box [m]."""
175 b: float
176 """Longer side [m]."""
179def plate_geometry(pl: Plate) -> PlateGeometry:
180 hm = pl.mean_thickness
181 return PlateGeometry(pl.area, outline_perimeter(pl.grid), hm, float(pl.material.bending_constant(hm)), coincidence_frequency(pl.material, hm), min(pl.Lx, pl.Ly), max(pl.Lx, pl.Ly))
185class LossEstimate:
186 """The three contributions to the decay rate of a bending partial [1/s], and their sum."""
188 material: float
189 thermoelastic: float
190 radiation: float
191 total: float
194def estimated_rate(pl: Plate, geom: PlateGeometry, f: float) -> LossEstimate:
195 """The decay rate of a bending partial at frequency ``f`` estimated without
196 the modes: the material rate, the thermoelastic rate averaged over the
197 plate area with the local thickness, and the radiation rate from the
198 average efficiency at the mean thickness. The modal computation replaces
199 the last two by exact per-mode values; this estimate serves the full
200 model's loss operator and the summary of a design."""
201 omega = 2 * math.pi * f
202 material = material_rate(pl.material.loss_factor, omega)
203 te = 0.0
204 if pl.material.thermal is not None:
205 area = pl.grid.node_area
206 sel = area > 0
207 eta = thermoelastic_loss_factor(pl.material, pl.h[sel], omega, pl.temperature)
208 te = float((area[sel] * eta).sum() / area[sel].sum()) * omega / 2
209 radiation = radiation_rate_from_efficiency(average_radiation_efficiency(f, geom.fc, geom.area, geom.perimeter, geom.a, geom.b), pl.material, geom.h_mean)
210 return LossEstimate(material, te, float(radiation), material + te + float(radiation))
213def in_plane_rate(pl: Plate) -> float:
214 """The in-plane motion of the full model has no radiation and no thermoelastic
215 loss in this model; it decays by internal friction at the rate of the lowest
216 in-plane resonance, c_p / (2 L) of the longer side."""
217 f_ip = pl.material.in_plane_speed / (2 * max(pl.Lx, pl.Ly))
218 return material_rate(pl.material.loss_factor, 2 * math.pi * f_ip)
221# ---------- The two-parameter operator of the full model ----------
225class SigmaFit:
226 """The two-parameter loss model of the full model: a bending mode of angular
227 frequency omega decays at the rate sigma0 + sigma1 k^2 with k^2 = omega / kappa."""
229 sigma0: float
230 """[1/s]"""
231 sigma1: float
232 """[m^2/s]"""
233 kappa: float
234 """Bending constant used for the fit [m^2/s]."""
237def fit_sigma(pl: Plate, geom: PlateGeometry | None = None, f_lo: float = 50.0, f_hi: float = 10000.0, n: int = 40) -> SigmaFit:
238 """Non-negative least-squares fit of sigma0 + sigma1 k^2 to the estimated
239 rates on ``n`` log-spaced frequencies from ``f_lo`` to ``f_hi``, at the mean thickness."""
240 geom = geom or plate_geometry(pl)
241 fs = f_lo * (f_hi / f_lo) ** (np.arange(n) / (n - 1))
242 x = 2 * math.pi * fs / geom.kappa
243 y = np.array([estimated_rate(pl, geom, f).total for f in fs])
244 return fit_line(x, y, geom.kappa)
247def fit_line(x, y, kappa: float) -> SigmaFit:
248 """Non-negative least squares of y = s0 + s1 x: with two parameters the
249 solution is the unconstrained one when it is feasible, else the best of
250 the fits with one parameter at zero, else zero."""
251 x = np.asarray(x, float)
252 y = np.asarray(y, float)
253 n = len(x)
254 sx, sy, sxx, sxy = x.sum(), y.sum(), (x * x).sum(), (x * y).sum()
255 det = n * sxx - sx * sx
256 candidates = [(0.0, 0.0), (sy / n, 0.0)]
257 if det > 0:
258 s1 = (n * sxy - sx * sy) / det
259 candidates.append(((sy - s1 * sx) / n, s1))
260 if sxx > 0:
261 candidates.append((0.0, sxy / sxx))
262 best, best_res = (0.0, 0.0), math.inf
263 for s0, s1 in candidates:
264 if not (s0 >= 0 and s1 >= 0):
265 continue
266 res = float(((y - s0 - s1 * x) ** 2).sum())
267 if res < best_res:
268 best, best_res = (s0, s1), res
269 return SigmaFit(best[0], best[1], kappa)
272def fitted_rate(fit: SigmaFit, f):
273 """Decay rate [1/s] of the fitted model at a frequency."""
274 return fit.sigma0 + fit.sigma1 * 2 * math.pi * np.asarray(f, float) / fit.kappa