Stress analysis of bicycle wheels using the Mode Matrix method.

bike-wheel-calc is a Python module for simulating the quasistatic structural behavior of wire-spoked bicycle wheels. Using the ModeMatrix object, the deformation of the wheel can be calculated given a set of external forces applied to the rim. In addition, predefined routines in the theory submodule may be used to calculate the wheel stiffness, buckling tension, and other properties.

Some example calculations that can be performed with bike-wheel-calc are:

• calculating change in spoke tension under rider weight, acceleration, braking, etc
• calculating the wheel stiffness, including rotational, lateral, and radial
• determining the effect of breaking a spoke
• comparing the performance of different spoke lacing patterns
• comparing the response of drive-side vs. non-drive-side spokes

You can try out bike-wheel-calc with an interactive, graphical, online app for simulating wheels based on the code here.

1. Download the entire repository.
2. Open a console (with Python) and navigate to the root directory
3. Install the package using

pip install .

The BicycleWheel objects defines the geometry, material properties, and spoke arrangement of the wheel. First, create an empty wheel object:

wheel = BicycleWheel()

Next, define the hub using the subclass Hub:

wheel.hub = wheel.Hub(diameter=0.040, width=0.050)

This creates a hub with a flange diameter of 40 mm and a total width of 50 mm (flange to flange). Optionally, you can specify diameter_ds, diameter_nds, width_ds, and width_nds. to create a hub with different flange diameters and/or rim dish.

Next, define the rim.

wheel.rim = Rim(radius=0.3,        # [m] radius at the beam axis
                area=100.0e-6,     # [m^2] cross-sectional area of the rim
                I_lat=1500e-12,    # [m^4] Second moment of area for lateral bending
                I_rad=3000e-12,    # [m^4] Second moment of area for radial bending
                J_tor=500e-12,     # [m^4] Torsion constant
                Iw=0.0,            # [m^6] Warping constant
                young_mod=69.0e9,  # [N/m^2] Young's modulus
                shear_mod=26.0e9)  # [N/m^2] Shear modulus

The torsion constant, and second moments of area are geometric properties of the rim cross-section that determine its stiffness. Iw is the warping constant (set to zero if you are unsure), young_mod is the rim material Young's Modulus, and shear_mod is the rim material Shear Modulus.

Next, create spokes by either using a predefined spoke configuration:

wheel.lace_cross(n_spokes=36, n_cross=3, diameter=2.0e-3, young_mod=210e9)

or by defining spokes manually:

rim_pt = (r_rim, theta_rim, z_rim)
hub_pt = (r_hub, theta_hub, z_hub)
wheel.spokes.append(wheel.Spoke(rim_pt, hub_pt, diameter=2.0e-3, young_mod=210e9))
...

The first two arguments are tuples defining the position of the spoke nipple and hub eyelet, respectively, in polar coordinates (R, theta, z). The spoke nipple does not need to lie on the rim centroid line. For example, a "fat-bike" wheel typically has spoke nipples which are offset from the centerline of the rim to provide torsional stability to the rim.

(Optional) Finally, apply spoke tension

wheel.apply_tension(800.)  # Apply 800 Newtons average radial tension.

The Mode Matrix Method, described in [1], is a numerical technique for solving the equations of motion of the wheel to arbitrary precision. In this method, the deformations functions are approximated by sine and cosine functions, e.g.

u(t) = u_0 + SUM(u_n_c*cos(n*theta) + u_n_s*sin(n*theta))
v(t) = v_0 + SUM(v_n_c*cos(n*theta) + v_n_s*sin(n*theta))
w(t) = ...
phi(t) = ...

where n goes from 1 to N in the SUM. u,v,w,phi are the lateral, radial, and tangential displacements and twist angle of the rim, respectively. The highest mode number N is chosen to achieve the desired precision. More modes requires more computational power, but results in a more precise solution.

In the Mode Matrix method, we solve for the coefficients u_0, u_1_c, u_1_s, u_2_c, u_2_s, .... If the highest mode is N, there are a total of 8N + 4 modes: That's one coefficient for each sine/cosine, each degree of freedom, and each mode = 2*4*N, plus 4 coefficients for u_0, v_0, w_0, phi_0 for the zero-mode.

The governing matrix equation of the wheel is

(K_rim + K_spk)*d = F_ext + A_adj*a

where K_rim + K_spk is the wheel modal stiffness matrix, d is the vector of mode coefficients u_0, u_1_c, ..., F_ext is the modal external force vector, A_adj is the spoke adjustment matrix, and a is the vector of spoke adjustments (change in length due to spoke nipple rotation only).

The real displacements are related to the mode coefficient vector d by a linear transformation given by:

u = B_u * d
v = B_v * d
...

The B matrix is calculated by the method B_theta(theta, components) of the ModeMatrix object.

The steps to solving a basic stress-analysis problem are as follows:

The ModeMatrix class contains the methods and objects to implement the the Mode Matrix method for stress analysis of the wheel, developed by Matthew Ford in his Ph.D. thesis [1].

Create a ModeMatrix object as follows:

mm = ModeMatrix(wheel, N=24)

The mode stiffness matrix is composed of two parts: the rim stiffness matrix and the spoke stiffness matrix:

K = (mm.K_rim(tension=True, r0=True) +
     mm.K_spk(tension=True, smeared_spokes=True)

The option tension=True on the rim stiffness matrix takes into account the effect of spoke tension and the compressive stress in the rim on lateral stiffness. The option tension=True on the spoke stiffness matrix takes into account the stiffening effect of tension on the spoke system lateral stiffness. The option smeared_spokes=True indicates that the wheel should be treated as if it had infinite spokes, smeared out into an effective "disc" with the same average stiffness as the original spokes. Choose smeared_spokes=False to get all the effects of discrete spokes.

Apply forces and torques by creating a F_ext object for each force, and adding together as many forces as desired.

# 120 Newton lateral force and 500 Newton radial force from the road
F_contact = mm.F_ext(theta=0., f=[120., 500., 0., 0.])

# Pair of forces from rim braking
F_brake_1 = mm.F_ext(theta=0., f=[0., 0., -100., 0.])     # Tangential force at the road
F_brake_2 = mm.F_ext(theta=3.1415, f=[0., 0., 100., 0.])  # Tangential force at the rim brakes

F = F_contact + F_brake_1 + F_brake_2

In this scenario, we have a created a radial and lateral force at the road contact point (theta=0), and a pair of forces representing a rim brakig scenario: -100 Newtons at the road contact point and 100 Newtons at the brake pads (theta=pi)

Note, you can also add a torque to the rim (say, if one flange of the rim is loaded) by specifying f=[0., 0., 0., T].

Solve the governing equation for the mode coefficients d using

d = numpy.linalg.solve(K, F)

Once you have d, evaluate the rim displacement at arbitrary points using

lat_disp = mm.rim_def_lat(theta, d)
rad_disp = mm.rim_def_rad(theta, d)

theta is a 1-D numpy array of locations along the rim.

Forces and displacements are specified in terms of three degrees of freedom: lateral: u, radial: v, and tangential: w, and a twist angle of the rim about its own axis: phi. The positive radial direction points inwards towards the center of the hub. The positive lateral direction points towards the non-drive side (NDS) of the bike. The positive tangential direction points along the rim in the counter-clockwise direction (forwards, at the road contact point). Positive rim twist is defined as a right-handed rotation around the tangential direction (imagine pointing your thumb along the tangential direction, and curling your fingers: that's the direction of positive rim twist).

If the lace_cross() method is used, the bottom-most spoke is a non-drive-side leading (pushing) spoke. The spokes alternate NDS, DS, NDS, ... and leading, leading, trailing, trailing...

[1] Matthew Ford, Reinventing the Wheel: Stress Analysis, Stability, and Optimization of the Bicycle Wheel, Ph.D. Thesis, Northwestern University (2018)