Profile with the cosine-dome bell (cymbal-bell-v1); the hole is closed, documented
2 changed files+56−14
solvers/vonkarman-modal-disk/__pycache__/vk_modal_disk.cpython-313.pycmodified+0−0View file
Binary file not shown.
solvers/vonkarman-modal-disk/vk_modal_disk.pymodified+56−14View file
@@ -27,8 +27,8 @@ Because the strike is centered on the x axis, the solution is even in y, and by
2727 modes needed for a given frequency reach (the coupling cost scales with the cube of the count). --full
2828 uses the unrestricted bases instead.
2929
30-If the challenge has a shell-of-revolution geometry (cymbal-bow-v1: a paraboloidal bow z = Z(r)), the same
31-solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
30+If the challenge has a shell-of-revolution geometry (cymbal-bow-v1: a paraboloidal bow z = Z(r);
31+cymbal-bell-v1: the bow plus a cosine-dome bell), the same solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
3232 system: L(w + w0, Phi) on the right of the bending equation and 2 L(w0, w) in the compatibility equation.
3333 With G^p_j = int Psi_p L(w0, Phi_j) dA this gives a linear membrane stiffening 2 sum_p |c_p| G^p (G^p)^T,
3434 which is diagonalized once in the flat-mode basis to give the shell modes (exact in the Gautschi step,
@@ -36,7 +36,9 @@ damped at sigma0 + eta omega^2 / 2 like every mode, as the Kelvin-Voigt law of t
3636 quadratic coupling terms F_k = sum_p c_p [ (r^T H^p r + 2 G^p . r) (H^p r)_k + (r^T H^p r) G^p_k ]. The strike
3737 acts along -z on the curved top face with the normalization of the specification, the listener integral uses
3838 the distance to the top face, and the field export is the normal displacement n_z w. Shallow-shell theory
39-neglects in-plane inertia and terms of the order of the slope squared (0.06 at the rim here).
39+neglects in-plane inertia and terms of the order of the slope squared (0.06 at the rim of the bow, 0.19 on the
40+flank of the bell). The mounting hole of the bell challenge is not represented: the solver closes it, which
41+neglects the mass and stiffness of a 12 mm disk at the apex; a stated approximation.
4042
4143 With --device cuda the coupling tensor is assembled and the time stepping runs on a GPU through PyTorch
4244 (the modes and quadrature stay in numpy/scipy); the CPU path is the reference and the two agree to round-off.
@@ -58,11 +60,14 @@ sys.path.insert(0, str(Path(__file__).resolve().parent))
5860 from bfs_disk import BFSDisk, load_challenge # noqa: E402
5961
6062
61-class Paraboloid:
62- """The bow of cymbal-bow-v1: mid-surface z = Z(r) = H (1 - r^2 / R^2), uniform thickness along the normal."""
63+class Profile:
64+ """An axisymmetric mid-surface z = Z(r): the paraboloidal bow of cymbal-bow-v1, optionally with the cosine-dome
65+ bell of cymbal-bell-v1. Uniform thickness along the normal. The hole of the bell challenge is not represented
66+ (the solver closes it; see the module docstring)."""
6367
64- def __init__(self, R: float, H: float, h: float):
68+ def __init__(self, R: float, H: float, h: float, bell_radius: float = 0.0, bell_height: float = 0.0, hole_radius: float = 0.0):
6569 self.R, self.H, self.h = R, H, h
70+ self.rb, self.hb, self.rh = bell_radius, bell_height, hole_radius
6671
6772 @staticmethod
6873 def from_challenge(ch: dict):
@@ -70,19 +75,49 @@ class Paraboloid:
7075 if not g or g.get("type") != "shell-of-revolution":
7176 return None
7277 p = ch["parameters"]
73- return Paraboloid(p["disk_radius"], p["bow_height"], p["h"])
78+ return Profile(p["disk_radius"], p["bow_height"], p["h"], p.get("bell_radius", 0.0), p.get("bell_height", 0.0), p.get("hole_radius", 0.0))
79+
80+ # ---- radial profile and its derivatives
81+ def _bell(self, r):
82+ if self.hb == 0.0:
83+ return np.zeros_like(r), np.zeros_like(r), np.zeros_like(r)
84+ k = np.pi / self.rb
85+ inside = r < self.rb
86+ z = np.where(inside, 0.5 * self.hb * (1.0 + np.cos(k * r)), 0.0)
87+ dz = np.where(inside, -0.5 * self.hb * k * np.sin(k * r), 0.0)
88+ d2z = np.where(inside, -0.5 * self.hb * k * k * np.cos(k * r), 0.0)
89+ return z, dz, d2z
90+
91+ def radial(self, r):
92+ """Z, Z' and Z'' as functions of r."""
93+ r = np.asarray(r, dtype=float)
94+ c = -2.0 * self.H / self.R**2
95+ z = self.H + 0.5 * c * r * r
96+ dz = c * r
97+ d2z = c * np.ones_like(r)
98+ zb, dzb, d2zb = self._bell(r)
99+ return z + zb, dz + dzb, d2z + d2zb
74100
75101 def z(self, x, y):
76- return self.H * (1.0 - (x * x + y * y) / self.R**2)
102+ return self.radial(np.sqrt(x * x + y * y))[0]
77103
78104 def grad(self, x, y):
79- c = -2.0 * self.H / self.R**2
80- return c * x, c * y
105+ r = np.sqrt(x * x + y * y)
106+ _, dz, _ = self.radial(r)
107+ s = np.where(r > 0, dz / np.where(r > 0, r, 1.0), 0.0) # Z'/r, finite on the axis
108+ return s * x, s * y
81109
82110 def curvatures(self, x, y):
83111 """Second derivatives (Z_xx, Z_yy, Z_xy) of the mid-surface height."""
84- c = -2.0 * self.H / self.R**2
85- return c * np.ones_like(x), c * np.ones_like(x), np.zeros_like(x)
112+ r = np.sqrt(x * x + y * y)
113+ _, dz, d2z = self.radial(r)
114+ rr = np.where(r > 0, r, 1.0)
115+ s = np.where(r > 0, dz / rr, d2z) # Z'/r -> Z'' on the axis
116+ cx, cy = x / rr, y / rr
117+ zxx = d2z * cx * cx + s * cy * cy
118+ zyy = d2z * cy * cy + s * cx * cx
119+ zxy = (d2z - s) * cx * cy
120+ return zxx, zyy, zxy
86121
87122 def normal(self, x, y):
88123 zx, zy = self.grad(x, y)
@@ -94,6 +129,13 @@ class Paraboloid:
94129 nx, ny, nz = self.normal(x, y)
95130 return x + 0.5 * self.h * nx, y + 0.5 * self.h * ny, self.z(x, y) + 0.5 * self.h * nz
96131
132+ def max_slope(self):
133+ r = np.linspace(0, self.R, 20001)
134+ return float(np.abs(self.radial(r)[1]).max())
135+
136+
137+Paraboloid = Profile # former name
138+
97139
98140 class VonKarmanDisk:
99141 def __init__(self, ch: dict, spacing: float, n_w: int, n_airy: int, force_override: float | None = None, log=None, symmetric_y: bool | None = None, device: str = "cpu", duration_override: float | None = None):
@@ -116,7 +158,7 @@ class VonKarmanDisk:
116158 if symmetric_y is None:
117159 symmetric_y = p["strike_y"] == 0.0
118160 self.symmetric_y = symmetric_y
119- self.profile = Paraboloid.from_challenge(ch)
161+ self.profile = Profile.from_challenge(ch)
120162 t0 = time.perf_counter()
121163 self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
122164 pl = self.plate
@@ -174,7 +216,7 @@ class VonKarmanDisk:
174216 t0 = time.perf_counter()
175217 self._apply_profile(xq, yq, WQ, PsiQ, Pxx, Pyy, Pxy)
176218 V = self.V
177- self.log(f" shell modes: f_low={np.round(self.omega[self.omega > 0][:4] / (2 * np.pi), 2)} Hz, f_max={self.omega.max() / (2 * np.pi):.0f} Hz in {time.perf_counter() - t0:.1f}s")
219+ self.log(f" shell modes (max slope {self.profile.max_slope():.2f}): f_low={np.round(self.omega[self.omega > 0][:4] / (2 * np.pi), 2)} Hz, f_max={self.omega.max() / (2 * np.pi):.0f} Hz in {time.perf_counter() - t0:.1f}s")
178220 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
179221 _, _, self.tau, self.cq = pl.listener_quadrature(3)
180222 self.PhiQ = Bn @ V # (Q, n_w)