3688ceacymsim: Python port of cymbalsim with a modal linear model and a PyTorch von Karman modelJeremy Magland 1"""Per-mode decay rates of the linear model from the mode shapes: the radiation
2into the air and the thermoelastic loss. Both are exact for the discrete model
3given the shape, so that the decay of each mode is a result of the geometry and
4the material rather than a prescribed curve.
6Radiation. A baffled plate whose nodes move with the velocities v_k radiates
7the power P = v^T R v with the radiation resistance matrix
8R_jk = (rho_air omega^2 / (4 pi c)) dA_j dA_k sinc(k r_jk) (Elliott and Johnson
91993), one face; both faces double it. For a mass-normalized mode of amplitude
10A the velocity is omega A phi, the energy is omega^2 A^2 / 2, and the amplitude
11decays at sigma = P / (2 E), which gives
12 sigma_rad = (rho_air omega^2 / (2 pi c)) sum_jk a_j a_k sinc(k r_jk),
13with a_k = dA_k phi(k). On a regular grid the kernel depends only on the offset
14j - k, so the double sum is a convolution, and by Parseval it is a sum over the
15spectrum: two 2-D FFTs on a zero-padded grid per mode instead of N^2 terms.
17Thermoelastic. Zener's loss factor applies to the bending energy and depends on
18the local thickness, so a mode's rate is the sum over the nodes of the local
19loss factor times the local bending energy, divided by omega: for a
20mass-normalized mode the total strain energy is omega^2 / 2, and a flat uniform
21plate recovers eta_te omega / 2.
22"""
24from __future__ import annotations
26import math
28import numpy as np
29import scipy.fft
31from ..constants import RHO_AIR, SPEED_OF_SOUND
32from ..loss import thermoelastic_loss_factor
33from ..plate import Plate
34from .modes import Modes
37def radiation_rates(modes: Modes, batch: int = 64) -> np.ndarray:
38 """The radiation decay rates [1/s] of all modes by the FFT convolution."""
39 pl = modes.plate
40 g = pl.grid
41 nx, ny, hx, hy = g.nx, g.ny, g.hx, g.hy
42 Px, Py = scipy.fft.next_fast_len(2 * nx - 1, real=True), scipy.fft.next_fast_len(2 * ny - 1, real=True)
43 c = SPEED_OF_SOUND
44 dx = np.arange(Px)
45 dx = np.where(dx <= Px // 2, dx, dx - Px)
46 dy = np.arange(Py)
47 dy = np.where(dy <= Py // 2, dy, dy - Py)
48 r = np.hypot(dx[None, :] * hx, dy[:, None] * hy)
49 reach = (np.abs(dx)[None, :] <= nx - 1) & (np.abs(dy)[:, None] <= ny - 1)
50 area = g.node_area
51 out = np.zeros(modes.count)
52 for start in range(0, modes.count, batch):
53 stop = min(modes.count, start + batch)
54 omega = 2 * math.pi * modes.freqs[start:stop]
55 k = omega / c
56 a = np.zeros((stop - start, Py, Px))
57 a[:, :ny, :nx] = area[None] * modes.w[start:stop]
58 A = scipy.fft.rfft2(a, axes=(-2, -1))
59 x = k[:, None, None] * r[None]
60 kern = np.where(x < 1e-9, 1.0, np.sin(x) / np.where(x < 1e-9, 1.0, x)) * reach[None]
61 S = scipy.fft.rfft2(kern, axes=(-2, -1)).real
62 # Parseval over the half spectrum: interior columns count twice.
63 wcol = np.full(Px // 2 + 1, 2.0)
64 wcol[0] = 1.0
65 if Px % 2 == 0:
66 wcol[-1] = 1.0
67 total = (S * (A.real**2 + A.imag**2) * wcol[None, None, :]).sum(axis=(-2, -1)) / (Px * Py)
68 out[start:stop] = RHO_AIR * omega**2 / (2 * math.pi * c) * total
69 return out
72def radiation_rate_direct(pl: Plate, shape: np.ndarray, freq: float) -> float:
73 """The same rate by the direct double sum, for the tests and small plates."""
74 g = pl.grid
75 sel = g.node_area > 0
76 a = (g.node_area * shape)[sel]
77 xs, ys = g.X[sel], g.Y[sel]
78 omega = 2 * math.pi * freq
79 k = omega / SPEED_OF_SOUND
80 x = k * np.hypot(xs[:, None] - xs[None, :], ys[:, None] - ys[None, :])
81 sinc = np.where(x < 1e-9, 1.0, np.sin(x) / np.where(x < 1e-9, 1.0, x))
82 return float(RHO_AIR * omega**2 / (2 * math.pi * SPEED_OF_SOUND) * (a[:, None] * a[None, :] * sinc).sum())
85def bending_energy_per_node(pl: Plate, w: np.ndarray) -> np.ndarray:
86 """The bending energy of transverse shapes ``w`` ``(..., ny, nx)`` on the grid,
87 split over the nodes (the mixed term of each cell is shared by its four
88 corners), as in the discrete energy of the full model."""
89 g, m = pl.grid, pl.material
90 hx, hy, nu = g.hx, g.hy, m.nu
91 left, right, down, up = g.neighbor_masks()
92 has_a = left & right
93 has_b = down & up
94 wp = np.pad(w, [(0, 0)] * (w.ndim - 2) + [(1, 1), (1, 1)])
95 a = np.where(has_a, (wp[..., 1:-1, 2:] - 2 * w + wp[..., 1:-1, :-2]) / hx**2, 0.0)
96 b = np.where(has_b, (wp[..., 2:, 1:-1] - 2 * w + wp[..., :-2, 1:-1]) / hy**2, 0.0)
97 U = 0.5 * hx * hy * g.node_weight * m.bending_stiffness(pl.h) * (a * a + b * b + 2 * nu * a * b)
98 cc = (w[..., 1:, 1:] - w[..., 1:, :-1] - w[..., :-1, 1:] + w[..., :-1, :-1]) / (hx * hy)
99 Uc = 0.5 * hx * hy * 2 * (1 - nu) * m.bending_stiffness(pl.cell_thickness) * cc * cc * g.cell_mask / 4
100 U[..., :-1, :-1] += Uc
101 U[..., :-1, 1:] += Uc
102 U[..., 1:, 1:] += Uc
103 U[..., 1:, :-1] += Uc
104 return U
107def thermoelastic_rates(modes: Modes, batch: int = 64) -> np.ndarray:
108 """The thermoelastic decay rates [1/s]: sum_k eta_te(omega, h_k) U_k / omega with U_k the bending energy at node k."""
109 pl = modes.plate
110 out = np.zeros(modes.count)
111 if pl.material.thermal is None:
112 return out
113 sel = pl.grid.node_area > 0
114 h = pl.h[sel]
115 for start in range(0, modes.count, batch):
116 stop = min(modes.count, start + batch)
117 omega = 2 * math.pi * modes.freqs[start:stop]
118 U = bending_energy_per_node(pl, modes.w[start:stop])[:, sel]
119 eta = thermoelastic_loss_factor(pl.material, h[None, :], omega[:, None], pl.temperature)
120 out[start:stop] = (eta * U).sum(axis=1) / omega
121 return out