Club Fitting Simulation: Forward Dynamics, Twist Counterfactuals, and Reproducible Reporting
In Layman’s Terms
Changing a club changes the mechanical system the golfer moves. It can also change how the golfer moves it. Simulation helps separate specified versions of these questions: what happens under the same applied torque history, under the same feedback policy, or under the same handle motion? Those are different experiments, with different quantities held fixed.
A repeatable simulation is valuable because its assumptions and calculations can be inspected. It does not predict a person’s response exactly. Camera errors, uncertain shaft properties, unknown grip loads, imperfect contact models and adaptation can all change the observed result. Good fitting therefore combines controlled measurements with models whose limits are explicit.
The chain is geometry and material identification → flexible response → delivery state → contact and ball flight → a fitting decision. An error at an earlier link can appear as an apparently convincing explanation at a later one. This article develops the mechanical links and a proposed reporting structure; it does not report a validated fitting product or human trial.
Why This Article Exists
A lower shot with a stiffer shaft could reflect different shaft deformation, changed handle motion, a different strike location, measurement variation, or several of these together. Launch data alone do not identify which explanation is responsible. A model makes alternatives explicit and suggests additional measurements that can distinguish them.
The parameters include mass, center of mass, rotational inertia, shaft geometry, bending and torsional stiffness, damping, grip boundary conditions and the contact model. Identifying all of them from a few launch outcomes is generally underdetermined. Measure accessible properties independently, examine parameter sensitivity and correlation, and reserve trials that were not used for calibration to assess predictions.
For a declared effective plant, a control-affine representation may be written \(\dot x=f_p(x,t)+G_p(x,t)u\). State \(x\) includes the chosen rigid, flexible and actuator states; \(p\) identifies model parameters. Explicit time permits a declared prescribed boundary motion. Within a smooth mode, control-affine structure depends on the chosen input: muscle activation, muscle force and joint torque are not interchangeable signals. Contact transitions need separately declared event and reset laws.
A Zero-Torque Counterfactual sets the declared applied control channels to zero. Gravity, passive elastic loads, retained actuators, contact reactions and prescribed-boundary work do not vanish merely because those channels are zero. Replaying a nonzero input on different equipment is an equipment counterfactual, not a zero-input trajectory. The Drift-Control Ratio compares projected drift acceleration with bounded control authority under a stated metric and admissible input set; it does not rank clubs or identify adaptation on its own. See the notation contract.
Part I: Delivery Kinematics and Equipment Counterfactuals
A Twist Requires Frames and a Reference Point
Immediately before contact, a rigid-head approximation has a pose and a velocity. A twist supplies the velocity part; it does not contain the head’s position, orientation, face shape or contact point. Here all vectors are measured relative to an inertial world frame and expressed in its axes. For a head-fixed point \(A\), use angular-first ordering:
\[ \boldsymbol{\xi}_A= \begin{pmatrix}\boldsymbol{\omega}\\\mathbf v_A\end{pmatrix}. \]
This six-vector representation of rigid-body velocity is often called a twist; its relation to an \(\mathfrak{se}(3)\) matrix depends on the selected spatial/body convention. For two points fixed in the same rigid head, let \(\mathbf r_{B/A}=\mathbf r_B-\mathbf r_A\) and define \([\mathbf r]_\times\mathbf a=\mathbf r\times\mathbf a\). Then
\[ \begin{aligned} \mathbf v_B&=\mathbf v_A+\boldsymbol\omega\times\mathbf r_{B/A},\\ \boldsymbol\xi_B&= \begin{pmatrix}I_3&0\\-[\mathbf{r}_{B/A}]_\times&I_3\end{pmatrix} \boldsymbol\xi_A. \end{aligned} \]
The minus sign follows from \(\boldsymbol\omega\times\mathbf r=-\mathbf r\times\boldsymbol\omega\). For \(\boldsymbol\omega=(0,0,2)\) rad/s and \(\mathbf r_{B/A}=(0.3,0,0)\) m, the added velocity is \((0,0.6,0)\) m/s. This is a point shift with unchanged axes, not a general rotation-and-translation frame change. Drake explicitly distinguishes the measured-in frame, expressed-in frame and reference point (Drake Developers, n.d.).
Wrenches must use the dual shift. With force \(\mathbf F\) and moment \(\boldsymbol\tau_A\), \(\boldsymbol\tau_B=\boldsymbol\tau_A-\mathbf r_{B/A}\times\mathbf F\). Consequently, \(\boldsymbol\tau_A\cdot\boldsymbol\omega+\mathbf F\cdot\mathbf v_A\) equals the corresponding power evaluated at \(B\). Changing the point in a velocity record without changing the moment reference can create a false energy-transfer claim.
What Is Held Fixed?
Let the baseline model use parameters \(p_0\), input history \(u_0(t)\) and state \(x_0(t)\). At a declared branch time \(t_b\), define an admissible initial-state map \(\Phi_{01}\) into candidate model \(p_1\). Equal numerical coordinates are not necessarily the same physical state when length, reference configuration, mass distribution or flexible coordinates change. Position, velocity, strain and stored energy cannot always all be held fixed simultaneously; record the choice and enforce the candidate’s constraints.
An open-loop input replay is
\[ \begin{aligned} \dot x_1&=f_{p_1}(x_1,t)+G_{p_1}(x_1,t)u_0(t),\\ x_1(t_b)&=\Phi_{01}(x_0(t_b)). \end{aligned} \]
Compare it with a baseline simulation using the same declared integration and initialization policy. Baseline fit error against measured motion must be reported separately. The replay isolates a parameter intervention within that model, conditional on the frozen input, initial-state mapping and retained loads. The input may be estimated rather than measured; inverse dynamics identifies net generalized loads only under its own model and external-force assumptions. It does not uniquely recover muscle activations.
Other questions require other interventions:
- Same feedback policy: \(u_1(t)=\pi(x_1(t),t)\). The policy is fixed, but its numerical output generally changes as the state changes. This differs from open-loop replay.
- Same handle motion: prescribe a measured pose trajectory at the grip and solve for shaft/head response and reaction loads. Input work can change. Prescribing all handle motions and all conjugate loads independently generally overconstrains the boundary.
- Zero declared control: integrate the specified zero-input model. A retained moving grip can still supply energy, and initial shaft strain can still release energy.
At a fixed diagnostic time \(t_f\), compare delivery twists at the same named geometric reference convention and in common axes. Club-specific centers of mass are different physical points, so their speeds are not automatically a face-center comparison. For actual impact, each branch must instead locate its own valid contact event. The candidate can arrive earlier, later, at a different face location, or miss the ball; evaluating both at the old impact time does not ensure that both are impacting.
An observed-minus-predicted residual is not an identified neural adaptation. It also contains measurement error, parameter uncertainty, initial-state mismatch, unmodeled loads and numerical error. Repeated, appropriately controlled trials and additional measurements are needed to study adaptation. Determinism makes a model replay auditable; it does not remove those alternatives.
Part II: Flexible Shaft Dynamics
Bending and Torsion Have Different Equations
Let \(s\in[0,L]\) denote the undeformed shaft coordinate. For one small transverse displacement \(w(s,t)\) about a declared reference motion, a bounded Euler–Bernoulli model with prescribed axial tension \(N(s,t)\) is
\[ \lambda w_{tt}+(EIw_{ss})_{ss}-(Nw_s)_s=q. \]
Here \(\lambda\) is mass per length in kg/m, \(EI\) is bending rigidity in N·m², \(N\) is axial force in N, and \(q\) is transverse load per length in N/m. Every term is a force density. A second bending direction needs its own rigidity and coordinates; anisotropic composite shafts may also require coupling terms. Euler–Bernoulli theory neglects shear deformation. Timoshenko theory introduces an independent section rotation, shear stiffness and rotary inertia; these names do not describe one interchangeable equation (Delft University of Technology, n.d.).
For a small Saint-Venant torsional angle \(\theta(s,t)\), a separate simple model is
\[ \mu_\theta\theta_{tt}-(GJ\theta_s)_s=m_s. \]
\(\mu_\theta=\int_A\rho r_\perp^2\,dA\) is polar mass moment per length in kg·m; \(GJ\) is effective torsional rigidity in N·m²; \(m_s\) is applied torque per length in N·m/m. For a noncircular or anisotropic section, the appropriate effective torsional rigidity is not obtained by blindly substituting the geometric polar second moment. Applying the bending operator and translational mass density to \(\theta\) is dimensionally wrong.
These equations are explanatory linearizations. A moving, rotating base requires consistent transport/inertial terms; large bending, bend–twist coupling, damping, gravity, grip compliance and head boundary inertia may require a richer formulation. The tip is not a free unloaded end when a clubhead is attached. Initial conditions, grip conditions and force/moment balance at the hosel close the model. Boundary excitation must not be counted twice in \(q\) and the prescribed motion.
A Centrifugal-Load Check
Consider only a straight radial shaft rotating steadily at angular speed \(\Omega\), with its root at radius \(R_0\) and a point head mass \(m_h\) at \(R_0+L\). Its tensile preload is
\[ \begin{aligned} N(s)&=\Omega^2\int_s^L\lambda(\sigma)(R_0+\sigma)\,d\sigma\\ &\quad+\Omega^2m_h(R_0+L). \end{aligned} \]
The acceleration is proportional to radius, not radius squared. Integrating \(\lambda v^2\) over length gives energy units, not force. For the synthetic values \(\lambda=0.06\) kg/m, \(L=1.1\) m, \(R_0=0.2\) m, \(\Omega=40\) rad/s and \(m_h=0.2\) kg, the base load is \(79.2+416=495.2\) N. This is a dimensional/numerical example, not a measured swing load or a formula for an arbitrarily oriented, accelerating golf shaft.
Lead/lag, droop and twist describe components in a declared shaft or delivery frame, not fixed world directions throughout the swing. A force whose line of action passes through the head center of mass contributes \(\mathbf M=\mathbf r_{CG/hosel}\times\mathbf F\) about the hosel. Other forces use their actual application points; any free couples must be added. The component in the selected bending plane depends on both vectors. An offset alone does not establish a toe-down angle. Lead need not always add loft, and torsion need not always delay face rotation: mapping shaft deformation into delivered face orientation depends on pose, mounting and the excitation history. No universal acceleration, droop or twist ranges are established here.
Discretization Must Preserve Coupling
A collection of diagonal springs is not automatically equivalent to a beam. For a uniform Euler–Bernoulli segment of length \(\ell\), clamp one end and use tip displacement \(w\) and slope \(\phi=dw/ds\) with their conjugate transverse force and moment. The remaining stiffness is
\[ \begin{pmatrix}F\\M\end{pmatrix} =\frac{EI}{\ell^3} \begin{pmatrix}12&-6\ell\\-6\ell&4\ell^2\end{pmatrix} \begin{pmatrix}w\\\phi\end{pmatrix}. \]
The off-diagonal terms are essential. With zero tip moment, solving both equations gives \(w=F\ell^3/(3EI)\) and \(\phi=F\ell^2/(2EI)\). Keeping only the translational diagonal gives \(F\ell^3/(12EI)\), four times too small for that boundary condition. If tip rotation is instead constrained, the diagonal translation result answers that different experiment. The complete two-node beam element retains both ends and their coupling (Woudenberg and Rocha, n.d.).
Similarly, uniform static torsion gives \(\theta=ML/(GJ)\) and stored energy \(GJ\theta^2/(2L)\). These limits are useful tests before simulating a swing. Select element count through convergence of delivery quantities, modes and energy balance; a fixed count such as 10–20 does not certify accuracy. Compare time-step refinement and damping models as well. Equivalent static stiffness alone does not ensure equivalent dynamic response.
Part III: Geometry, Density and Inertia
What a Watertight Mesh Can Establish
A consistently oriented, closed, non-self-intersecting mesh can define a polyhedral material region. An exterior scan alone usually does not identify wall thickness, cavities, internal weights or density. A hollow multimaterial head needs those regions and material assignments; watertightness alone does not reveal its mass distribution. Exact integration of a specified polyhedron is also distinct from exact knowledge of the manufactured part or curved CAD surface.
The divergence theorem reduces polynomial volume moments to surface integrals. Eberly’s primary derivation explains this reduction and its constant-density assumptions (Eberly 2009). The following signed-tetrahedron form provides a compact independent check. For each outward-oriented triangular face with vertices \(\mathbf a_k,\mathbf b_k,\mathbf c_k\) expressed from a common origin, define
\[ \begin{aligned} \Delta V_k=\frac{1}{6}\mathbf a_k\cdot(\mathbf b_k\times\mathbf c_k),\\ V=\sum_k\Delta V_k,\qquad m=\rho V,\\ \mathbf r_{CG}=\frac{1}{4V}\sum_k\Delta V_k (\mathbf a_k+\mathbf b_k+\mathbf c_k). \end{aligned} \]
Thus the centroid numerator written directly in scalar triple products carries \(1/24\), not \(1/4\). Individual signed volumes may be negative if the origin lies outside the body; taking their absolute values separately corrupts the integral. Cavity boundaries must have the orientation appropriate to the material region. Check manifold closure, units and orientation before interpreting the result; nearly cancelling terms may require a better origin or numerical summation strategy.
The unit tetrahedron with vertices \(0,\mathbf e_x,\mathbf e_y,\mathbf e_z\) has \(V=1/6\) and centroid \((1/4,1/4,1/4)\). These identities immediately detect the earlier factor errors. They also remain valid after translating the whole tetrahedron and its centroid.
For a density field \(\rho(\mathbf r)\) and an origin \(O\),
\[ \begin{aligned} I_O&=\int_\Omega\rho(\mathbf r) \left(\|\mathbf r\|^2I_3-\mathbf r\mathbf r^T\right)dV,\\ I_{CG}&=I_O-m\left(\|\mathbf r_{CG}\|^2I_3 -\mathbf r_{CG}\mathbf r_{CG}^T\right). \end{aligned} \]
Off-diagonal tensor entries are the negatives of the corresponding product integrals. For the unit tetrahedron at unit density, \(I_{CG}\) has diagonal entries \(1/80\) and off-diagonal entries \(1/480\), in the consistent chosen length/mass units. A tensor must be symmetric and physically realizable; nonnegative eigenvalues alone are insufficient because principal moments must also satisfy the triangle inequalities.
Part IV: Proposed Equipment Data
The next three JSON blocks are proposed data examples, not published wire standards. The earlier golf-club/fitting-document/v1, swing-sim/delivery-trajectory/v1 and golf-club/fitting-report/v1 names did not have matching schemas or implementations in this repository. These examples therefore omit schema URLs and compatibility guarantees. All identities and numbers below are synthetic; they are not OEM specifications, measured player records, qualified engine runs or purchase recommendations.
A future equipment contract should identify geometry/material artifacts, measurement methods, uncertainty, cut versus uncut shaft state, mounting transforms, stiffness/damping profiles and versioned unit conventions. Here a short SI example illustrates consistent axis naming:
{
"example_status": "synthetic_not_a_validated_run",
"proposal": "equipment_record",
"equipment_id": "example_club_A",
"head_frame": {
"origin": "head_center_of_mass",
"axes": "orthonormal right-handed; x face normal, y selected heel-toe direction, z=x cross y"
},
"head": {
"mass_kg": 0.1985,
"cg_from_hosel_expressed_in_head_m": [0.0382, 0.0145, -0.0128],
"inertia_about_cg_expressed_in_head_kg_m2": [
[0.000485, -0.0000085, 0.000012],
[-0.0000085, 0.000312, -0.0000045],
[0.000012, -0.0000045, 0.000541]
]
},
"shaft": {
"state": "cut_and_prepared_for_this_example",
"length_m": 1.143,
"mass_kg": 0.0621,
"stiffness_and_damping_profile": null,
"profile_status": "not_supplied_not_simulation_ready"
}
}The earlier g·cm² unit converts by \(1\ \mathrm{g\,cm^2}=10^{-7}\ \mathrm{kg\,m^2}\). Here the physical direction of the positive heel-toe axis must be selected consistently with the stated right-handed construction and the particular club. Do not independently label an inconsistent third axis “up.” An actual schema must specify missing-value behavior and reject incomplete records for simulation. A single shaft frequency does not identify its full mass, stiffness and damping distributions.
Part V: Proposed Delivery Records and Engine Comparisons
A trajectory record supports comparison but cannot reconstruct the governing model by itself. Exchange the model, actuator and constraint definitions, flexible discretization, initialization, input policy, contact parameters, integrator and solver settings as well. Record engine version, platform and numerical tolerances. MuJoCo’s model compiler and contact/constraint options are examples of backend-specific choices that a label such as “compatible” does not settle (MuJoCo Developers, n.d.). Cross-engine agreement needs executable comparison cases and reported errors; bitwise identity is a separate, usually much stricter requirement.
{
"example_status": "synthetic_not_a_validated_run",
"proposal": "delivery_excerpt",
"equipment_ref": "example_club_A",
"engine_run": null,
"world_axes": "orthonormal right-handed; x lateral, y target, z up",
"pose_convention": "head-to-world rotation; head origin at center of mass",
"velocity_measured_in": "world",
"velocity_expressed_in": "world",
"velocity_reference_point": "head_center_of_mass",
"time_reference": "declared diagnostic instant; contact not established",
"sampling_status": "two illustrative samples; no integration step implied",
"trajectory_samples": [
{
"time_s": -0.05,
"position_world_m": [0.124, -0.45, 1.25],
"orientation_quaternion_wxyz": [1.0, 0.0, 0.0, 0.0],
"angular_velocity_world_rad_s": [12.4, 5.8, 18.2],
"cg_velocity_world_m_s": [14.2, 28.5, -8.4],
"shaft_coordinates": {"lead_m": -0.0184, "droop_m": -0.0121, "twist_rad": -0.0314}
},
{
"time_s": 0.0,
"position_world_m": [0.0, 0.0, 0.025],
"orientation_quaternion_wxyz": [0.8, 0.0, 0.6, 0.0],
"angular_velocity_world_rad_s": [22.8, 8.4, 48.6],
"cg_velocity_world_m_s": [1.4, 48.2, 4.1],
"shaft_coordinates": {"lead_m": 0.0246, "droop_m": -0.0282, "twist_rad": 0.0070}
}
]
}The quaternions are unit length and their mapping direction is named. Their samples and velocities have not been integrated into a dynamically consistent trajectory. A production validator must distinguish such a structural example from usable motion: verify temporal ordering, quaternion/velocity consistency, reference-point transforms, interpolation and contact events. The shaft coordinates also need a fully specified moving basis, sign convention and reference configuration before use. Sparse reporting times are not solver time steps.
Part VI: Reproducible Reports and Statistical Meaning
A deterministic forward run produces an output conditional on inputs and model choices. Dispersion requires an ensemble: repeated measured swings, a declared input distribution, parameter uncertainty or some combination. Report which variation was sampled and which uncertainties were omitted. Separate validation data from calibration data.
A 95% confidence interval for a population mean is not a 95% interval for individual shots. For an illustrative independent normal sample with known standard deviation \(\sigma=4.2\) m and \(n=25\), the mean interval half-width is approximately \(1.96\sigma/\sqrt n=1.6464\) m. The central interval about the known population mean for individual outcomes has half-width \(1.96\sigma=8.232\) m. With estimated parameters, dependence or non-normal outcomes, use and describe an appropriate method; do not attach “95%” to an unexplained dispersion number.
{
"example_status": "synthetic_not_a_validated_run",
"proposal": "fitting_report",
"equipment_ref": "example_club_A",
"candidate_delivery_and_launch": null,
"uncertainty_example": {
"model": "independent normal observations with known standard deviation",
"sample_count": 25,
"known_carry_standard_deviation_m": 4.2,
"normal_quantile_approximation": 1.96,
"mean_confidence_interval_half_width_m": 1.6464,
"population_central_interval_half_width_m": 8.232
},
"sensitivity_estimates": null,
"physical_validation": "not_performed",
"recommendation": "withheld_pending_identification_and_validation"
}Sensitivity gradients must name their input units and intervention. A derivative with respect to “shaft EI per unit” is ambiguous for a spatial profile: define, for example, a dimensionless scale multiplying a fixed profile, and state which other parameters are held fixed. Perturbing head mass at fixed geometry also changes inertia unless a different artificial intervention is explicitly intended. Contact-event changes can make a local derivative discontinuous or inapplicable. Check step-size dependence before interpreting it as a design direction.
Delivery twist alone cannot determine launch. Relative contact velocity, face normal, strike location, ball properties, friction, face/ball compliance and the chosen contact law also enter. A scalar restitution parameter is an effective model input under defined conditions, not a universal material constant or equipment-conformance certificate. Carry adds aerodynamic and environmental assumptions. Optimality is therefore conditional on a declared objective, player population, constraints and uncertainty, rather than the largest simulated carry number.
Part VII: Verification Architecture
Start with independent checks at each link. Test signed-volume integration on analytic solids, verify frame shifts with power invariance, and check static bending/torsion limits before dynamic simulation. Establish convergence and energy balance, then compare predicted shaft motion and head delivery with measurements that were not used to fit the parameters. Test the impact and flight models against their own evidence rather than allowing their error to be absorbed into shaft stiffness. Finally assess repeatability and player response in a controlled fitting study.
These checks connect the article’s central questions. A center-of-mass error changes both the inertial loading and the conversion from center-of-mass velocity to face-point velocity. An omitted beam coupling changes deformation and reaction work. A shifted impact time changes which delivery state enters the collision model. Treating the resulting prediction residual as neural adaptation would hide all three mechanical alternatives.
The useful outcome is an inspectable chain of conditional predictions and tests. Empirical fitting supplies essential evidence for that chain; simulation helps identify what to measure and which explanations remain distinguishable.
References
References
Citation
@online{olson2026,
author = {Olson, Dieter},
title = {Club {Fitting} {Simulation:} {Forward} {Dynamics,} {Twist}
{Counterfactuals,} and {Reproducible} {Reporting}},
date = {2026-09-28},
url = {https://affinedrift.com/articles/technology-club-fitting.html},
langid = {en}
}