1function H = plate_hertz(params)
2% PLATE_HERTZ Hertz contact constants of the hammer tip on the plate.
3%
4% H = plate_hertz(params) returns
5% estar effective modulus, 1 / ((1 - nu_h^2) / E_h + (1 - nu^2) / E) [Pa]
6% contact_k Hertz stiffness 4/3 E* sqrt(R): force = k delta^(3/2) [N/m^1.5]
7% duration contact time of a sphere on a rigid half-space,
8% 2.94 (m^2 / (R E*^2 v))^(1/5) [s]
9% impulse momentum transferred, (1 + e) m v [N s]
10% shape_integral mean of the sin^(3/2) pulse shape over its duration,
11% Gamma(5/4) / (sqrt(pi) Gamma(7/4))
12% peak_force impulse / (duration * shape_integral) [N]
13% hunt_crossley the damping factor of the coupled model, 3 (1 - e) / (2 v) [s/m]
14%
15% The Rayleigh-Ritz method applies the pulse of this duration and impulse
16% directly; the finite-element method integrates the Hunt-Crossley force
17% k delta^(3/2) (1 + hunt_crossley delta_dot) against the moving plate and
18% uses the duration only to size its time step.
20hm = params.hammer;
21mat = params.material;
22inv_estar = (1 - hm.nu^2) / hm.E + (1 - mat.nu^2) / mat.E;
23H.estar = 1 / inv_estar;
24H.contact_k = 4 / 3 * H.estar * sqrt(hm.tip_radius);
25speed = max(hm.speed, 1e-9);
26H.duration = 2.94 * (hm.mass^2 / (hm.tip_radius * H.estar^2 * speed))^0.2;
27H.impulse = (1 + hm.restitution) * hm.mass * hm.speed;
28H.shape_integral = gamma(1.25) / (sqrt(pi) * gamma(1.75));
29H.peak_force = H.impulse / (H.duration * H.shape_integral);
30H.hunt_crossley = 1.5 * (1 - hm.restitution) / speed;
31end