How to model the dynamic tension distribution of a loaded bicycle wheel?

(Sorry for the very mathematical post, but code that is executable with no mathematics knowledge is at the end.)

First some mathematics. Let us consider the rim section having three spokes, two at both ends and one in the middle, as a simply supported beam (not necessarily 100% realistic model but should be good enough for most purposes). Its equation is

(F_S(0)-F(0))*L^3/(6*E*I) = w(0) - (w(-L) + w(L))/2

where:

• F_S(x) is the dynamic load of the spoke at location x,
• F(x) is the wheel load at location x,
• L is the distance between two neighboring spokes at the rim,
• E is the Young's modulus of aluminum (rim),
• I is the second moment of area of rim,
• w(x) is the radial deflection at location x.

There is also the equation for spokes:

w(0) = -F_S(0) L_S / (E_S A_S)

where:

• L_S is the length of a spoke,
• E_S is the Young's modulus of steel (spoke),
• A_S is the area of spoke cross section.

Let us also define:

k = L^3/(6*E*I),
k_S = L_S / (E_S A_S),
C = (1+k/k_S)

Then the equations can be combined and rewritten as:

C*w(0) - 0.5*w(-L) - 0.5*w(L) = -F(0)*k

and repeated 36 times in a circularly shifted form. This results in the following Matlab / GNU Octave code (I only tested it on GNU Octave but should work on Matlab as well):

L = 600*pi/36; % mm
L_S = 300; % mm
E = 69e3; % MPa
E_S = 193e3; % MPa
I = 1124; % mm4
A_S = pi*0.8^2; % mm2
k = L^3/(6*E*I); % mm/M
k_S = L_S/(E_S*A_S); % mm/N
C = (1+k/k_S);
m = [C, -0.5, zeros(1,33), -0.5];
M = [circshift(m, 0); circshift(m, 1); circshift(m, 2); circshift(m, 3);
     circshift(m, 4); circshift(m, 5); circshift(m, 6); circshift(m, 7);
     circshift(m, 8); circshift(m, 9); circshift(m, 10); circshift(m, 11);
     circshift(m, 12); circshift(m, 13); circshift(m, 14); circshift(m, 15);
     circshift(m, 16); circshift(m, 17); circshift(m, 18); circshift(m, 19);
     circshift(m, 20); circshift(m, 21); circshift(m, 22); circshift(m, 23);
     circshift(m, 24); circshift(m, 25); circshift(m, 26); circshift(m, 27);
     circshift(m, 28); circshift(m, 29); circshift(m, 30); circshift(m, 31);
     circshift(m, 32); circshift(m, 33); circshift(m, 34); circshift(m, 35)];

results = M[k*500;zeros(35,1)]
firstspoke = results(1)/k_S

When this code is executed, it says that the first spoke deflects 0.158 mm (Jobst Brandt's simulation showed the spoke where load is placed deflecting 0.153 mm so the simulation results are very similar) and carries 203.92 N of the 500 N load (calculated from Jobst Brandt's simulation, that simulation showed it carries 198.11 N of the 500 N load).

We can calculate from this that according to this simulation, one spoke carries 41% of load (Jobst Brandt: 40% of load) in the example wheel.

The difference between this simulation and Jobst Brandt's simulation is that this simulation considered only radial deflection so it is not necessarily as accurate, and also the rim as a simply supported beam might not be the most accurate possible model.