The Yeadon software was developed as a component in our tool chain for performing simulations of bicycle dynamics that include the effects of a human rider. We choose to model the bicycle rider such that they were rigidly attached to the rear frame of the bicycle thus we required estimates of the inertial properties of a human in a seated position. Yeadon's popular human inertia model is a great candidate for this purpose. This notebook describes the process of solving for the joint angles needed to configure a general Yeadon human model to sit on a typical bicycle with hands on the handlebars, arms hanging down, feet at the crank axle, and butt on the seat. Then, we express the inertial properties of the rider in this configuration in terms that allow them to be combined with the inertial properties of the rear frame of the bicycle.
from IPython.display import Image
import sympy as sym
import sympy.physics.mechanics as me
import numpy as np
from scipy.optimize import fsolve
# This makes sure the mayavi gui isn't blocking in IPython.
%gui wx
from mayavi import mlab
import yeadon
sym.init_printing(print_builtin=False)
%precision 5
u'%.5f'
We describe the basic bicycle geometry as follows:
Image('bicycle-geometry.png')
The human is described by 95 geometrical measurements using Yeadon's method and we will need a subset of them for the configuration computations:
Image('http://yeadon.readthedocs.org/en/latest/_images/measurements.png')
The human's configuration is defined by 21 joint angles, also described in Yeadon's method.
Image('http://yeadon.readthedocs.org/en/latest/_images/configuration.png')
Our first objective is to solve for the somersault
angle; the leg angles PJ1extension
and J1J2flexion
; and the arm angles CA1extension
, CA1adduction
, CA1rotation
and A1A2extension
that place the rider in the proper configuration on a bicycle. We assume the rider's configuration is symmetric so that the rider's right side configuration is a mirror of the left, for example, PJ1extension
is equal to PK1extension
.
We first configure the rider by adding geometric points of the bicycle to the main Yeadon reference frame that represent the desired location of the butt, hands, and feet. We set the center of the base of the thumb stadia $L_{a5}$ and $L_{b5}$ to lie on the handlebar grips, the origin of the Yeadon model to be at the top of the seat, and the centers of the heel of the foot stadia $L_{j6}$ and $L_{k6}$ to lie on the bottom bracket axis.
We first define a reference frame that represents the Yeadon model's base frame. The $\hat{x}$ unit vector is directed from right to left, the $\hat{y}$ unit vector directed anterior to posterior, and the $\hat{z}$ unit vector is directed inferior to superior. This is consistent with the configuration figure above.
yeadon_rf = me.ReferenceFrame('Y')
The origin, $O$, of the Yeadon model is at the base of pelvis.
origin = me.Point('O')
The pelvis, $P$, thorax, $T$, and chest, $C$, segements in Yeadon's model are assumed to be fixed relative to each other. The pelvis is then rotated forward through a somersault angle such that the rider's trunk ($P$,$T$, & $C$) is hunched over with respect to the bicycle. Experimentally, we determine the somersault angle, $\phi_P$, by photographing the seated bicyclist and estimating the angle at which the rider is leaned forward.
somersault = sym.symbols('phi_P')
pelvis = yeadon_rf.orientnew('P', 'Axis', (somersault, yeadon_rf.x))
This rotation and all remaining rotations are body-fixed $X$-$Y$-$Z$ rotations.
To place the heels of the rider, $L_{j6}$ and $L_{k6}$, on the bike pedals, we must solve for the required hip extension and knee flexion angles. The thigh is defined to extend at the hip through the angle PJ1extension
and the lower leg to flex at the knee through the angle J1J2flexion
. We choose hip adduction to be zero, so the centers of the heels are the same distance from the medial plane as the hip centers are. First, we'll describe the position of the heel from the pelvis. Next, we will describe the position of the bottom bracket from the seat. Then, we will constrain the two positions to be equal.
PJ1extension, J1J2flexion = sym.symbols('phi_J1, phi_J2')
We create two new reference frames for the left thigh, $J1$ and the left lower leg, $J2$.
left_thigh = pelvis.orientnew('J1', 'Axis', (PJ1extension, pelvis.x))
left_lower_leg = left_thigh.orientnew('J2', 'Axis', (J1J2flexion, left_thigh.x))
The heel center is located with respect to the seat via measurements of the subject when standing upright:
The heel point, $H$, is thuse defined:
thigh_length, Lj5L, Lj6L = sym.symbols('L_j3L, L_j5L, L_j6L')
lower_leg_length = (Lj5L - thigh_length) + Lj6L
heel = origin.locatenew('H', -thigh_length * left_thigh.z - lower_leg_length * left_lower_leg.z)
heel.pos_from(origin)
The seat point, $S$, is the point at the top of the seat that lies on the seat tube axis, and it is coincident with the origin, $O$, of the Yeadon model.
seat = origin.locatenew('S', 0)
The location from the origin to the bottom bracket can be described by the geometrical measurements:
seat_tube_length, seat_post_length, seat_tube_angle = sym.symbols('l_st, l_sp, lambda_st')
seat_distance = seat_tube_length + seat_post_length # relative to the bottom bracket
bottom_bracket = me.Point('B')
bottom_bracket.set_pos(origin, -seat_distance * (sym.cos(seat_tube_angle) * yeadon_rf.y +
sym.sin(seat_tube_angle) * yeadon_rf.z))
We find the hip extension and knee flexion by enforcing that the distance from the heel to the bottom bracket is zero. This gives two nonlinear equations in two unknowns.
leg_zero_vector = bottom_bracket.pos_from(heel)
zero1 = me.dot(leg_zero_vector, yeadon_rf.y)
zero2 = me.dot(leg_zero_vector, yeadon_rf.z)
zero1.trigsimp()
zero2.trigsimp()
These two equations can be solved for the hip extension angle, $\phi_{J1}$, and the knee flexion angle, $\phi_{J2}$, given all other variables (remember, we measure $\phi_{P}$).
The relationship between the bicycle geometry and the arm's coinfguration angles can be found in a similar manner, although now we require full three dimensional rotations.
The left shoulder, $A1$, is oriented with respect to the trunk through body-fixed $X$-$Y$-$Z$ rotations $\phi_{A1}$, $\theta_{A1}$, and $\psi_{A}$.
CA1extension, CA1adduction, CA1rotation = sym.symbols('phi_A1, theta_A1, psi_A')
left_upper_arm = pelvis.orientnew('A1', 'Body',
(CA1extension, CA1adduction, CA1rotation), 'XYZ')
The lower arm is oriented with respect to the upper arm, $A1$, through a simple rotation: A1A2extension
or $\phi_{A2}$.
A1A2extension = sym.symbols('phi_A2')
left_lower_arm = left_upper_arm.orientnew('A2', 'Axis', (A1A2extension, left_upper_arm.x))
To define the position of the hand, we require the following measurements, all of which come from the measurements figure above.
shoulder_width = sym.symbols('L_s4w')
torso_length, upper_arm_length, La4L, La5L = sym.symbols('L_s4L, L_a2L, L_a4L, L_a5L')
shoulder = origin.locatenew('La0', shoulder_width / 2 * pelvis.x + torso_length * pelvis.z)
lower_arm_length = (La4L - upper_arm_length) + La5L
hand = shoulder.locatenew('H', -upper_arm_length * left_upper_arm.z - lower_arm_length * left_lower_arm.z)
Now, we locate the handlebar grip points with respect to the seat. But several new bicycle geometric quantities are needed. We choose quantities that are simple to measure in practice, and will use these to solve for quantities that are more natural to use in our model.
We define four unknowns that give the $Y$ and $Z$ measures of the vectors from the front and rear wheel centers to the handlebar grips: $\alpha_{y}$, $\alpha_{z}$, $\beta_{y}$, $\beta_{z}$.
chain_stay_length = sym.symbols('l_cs')
front_wheel_radius, rear_wheel_radius = sym.symbols('r_F, r_R')
bottom_bracket_height = sym.symbols('h_bb')
handle_width = sym.symbols('w_hb')
wheelbase = sym.symbols('w')
alpha_y, alpha_z, beta_y, beta_z = sym.symbols('alpha_y, alpha_z, beta_y, beta_z')
rear_wheel_center = me.Point('cR')
measure_z = rear_wheel_radius - bottom_bracket_height
measure_y = sym.sqrt(chain_stay_length**2 - (measure_z)**2)
rear_wheel_center.set_pos(bottom_bracket,
measure_y * yeadon_rf.y +
measure_z * yeadon_rf.z)
front_wheel_center = me.Point('cF')
front_wheel_center.set_pos(rear_wheel_center,
-rear_wheel_radius * yeadon_rf.z -
wheelbase * yeadon_rf.y +
front_wheel_radius * yeadon_rf.z)
grip_defined_from_front = me.Point('hbF')
grip_defined_from_front.set_pos(front_wheel_center,
handle_width / 2.0 * yeadon_rf.x +
alpha_y * yeadon_rf.y +
alpha_z * yeadon_rf.z)
grip_defined_from_rear = me.Point('hbR')
grip_defined_from_rear.set_pos(rear_wheel_center,
handle_width / 2.0 * yeadon_rf.x +
beta_y * yeadon_rf.y +
beta_z * yeadon_rf.z)
The last two vectors are expressed in terms of unknowns. Now, we'll solve for those unknowns, using measured values for the distances between the wheel centers and the handle grip, $L_{hbF}, L_{hbR}$. We require the seat
point we defined in a previous section. We know that the magnitudes of the vectors from the wheel centers to the grip must equate to the measurements.
front_to_grip_distance, rear_to_grip_distance = sym.symbols('L_hbF, L_hbR')
zero3 = grip_defined_from_front.pos_from(front_wheel_center).magnitude() - front_to_grip_distance
zero4 = grip_defined_from_rear.pos_from(rear_wheel_center).magnitude() - rear_to_grip_distance
zero3.trigsimp()
zero4.trigsimp()
Furthermore, the two handlebar vectors relative to the seat must be equal.
bike_param_zero_vector = (grip_defined_from_front.pos_from(seat) -
grip_defined_from_rear.pos_from(seat))
zero5 = me.dot(bike_param_zero_vector, yeadon_rf.y)
zero6 = me.dot(bike_param_zero_vector, yeadon_rf.z)
zero5
zero6
These four equations can be solved for $\alpha_{y}$, $\alpha_{z}$, $\beta_{y}$, $\beta_{z}$. Now, we know the position of the handlebar grip point, through grip_defined_from_front
or grip_defined_from_rear
in terms of the bicycle's geometry.
Finally, we generate equations to solve for the arm angles with the additional assumption that the plane defined by axis of the upper arm is always parallel to the ground, $XY$ plane.
arm_zero_vector = hand.pos_from(origin) - grip_defined_from_rear.pos_from(origin)
zero7 = me.dot(arm_zero_vector, yeadon_rf.x)
zero8 = me.dot(arm_zero_vector, yeadon_rf.y)
zero9 = me.dot(arm_zero_vector, yeadon_rf.z)
zero7.trigsimp()
zero8.trigsimp()
zero9.trigsimp()
The arm x unit vector is orthogonal to Yeadon z unit vector (i.e. parallel to the ground).
zero10 = me.dot(left_upper_arm.x, yeadon_rf.z)
zero10.trigsimp()
These four equations allow us to solve for $\phi_{A1}$, $\theta_{A1}$, $\psi_{A}$, and $\phi_{A2}$.
We now have 10 equations that we can solve simultaneously for the left leg and left arm angles. This will generally require give good guesses in the numerical solver.
equations = [zero1, zero2, zero3, zero4, zero5, zero6, zero7, zero8, zero9, zero10]
We need to solve the ten equations for $\phi_{J1},\phi_{J2},\phi_{A1},\theta_{A1},\psi_{A1},\phi_{A2},\alpha_{y},\alpha_{z},\beta_{y},\beta_{z}$ given values for all of the measurements. The first two equations for the leg angles are not coupled to the remaining equations for the arms. Also, these are nonlinear equations; we will need to use a solver such as scipy.optimize.fsolve
to find the solutions to the equations.
Here we list all the known measurements, as well as all the unknowns for which we are solving. Note that the order is important.
unknown_symbols = (PJ1extension, J1J2flexion, CA1extension, CA1adduction, CA1rotation, A1A2extension,
alpha_y, alpha_z, beta_y, beta_z)
bicycle_measurement_symbols = (seat_tube_length, seat_post_length, seat_tube_angle, chain_stay_length,
front_wheel_radius, rear_wheel_radius, bottom_bracket_height, handle_width,
wheelbase, front_to_grip_distance, rear_to_grip_distance)
human_measurement_symbols = (thigh_length, Lj5L, Lj6L, torso_length, shoulder_width, upper_arm_length,
La4L, La5L, somersault)
all_symbols = unknown_symbols + bicycle_measurement_symbols + human_measurement_symbols
Now we generate a function that is based on the SymPy expressions and that can be evaluated numerically.
numerical_function = sym.lambdify(all_symbols, equations, modules="numpy")
These are typical values of the measurements, expressed in meters except otherwise noted (from the Browser bicycle: http://moorepants.github.io/dissertation/physicalparameters.html#parameter-tables).
n_lst = 0.53
n_lsp = 0.24
n_lamst = 1.195550538 # radians
n_lcs = 0.46
n_rF = 0.34352982332
n_rR = 0.340958858855
n_hbb = 0.295
n_whb = 0.58
n_w = 1.121
n_LhbF = 0.8930
n_LhbR = 0.9213
args = [n_lst, n_lsp, n_lamst, n_lcs, n_rF, n_rR, n_hbb, n_whb, n_w, n_LhbF, n_LhbR]
These typical measurements in meters from an adult rider for the Yeadon model (Jason's measurements) can be loaded from a YAML file using the Yeadon package.
human = yeadon.Human('JasonYeadonMeas.txt')
We then add the needed measurements to the args
variable.
for measurement in ['Lj3L', 'Lj5L', 'Lj6L', 'Ls4L', 'Ls4w', 'La2L', 'La4L', 'La5L']:
args.append(human.meas[measurement])
#n_somersault = 0.123918377 # radians
n_somersault = 0.175 # try 10 degrees instead of actual measurement because of incompatible geometry
args.append(n_somersault)
Now we must give initial guesses for all of the unknowns. These should be close to the correct solution:
g_PJ1extension = -np.deg2rad(90.0)
g_J1J2flexion = np.deg2rad(75.0)
g_CA1extension = -np.deg2rad(15.0)
g_CA1adduction = np.deg2rad(2.0)
g_CA1rotation = np.deg2rad(2.0)
g_A1A2extension = -np.deg2rad(40.0)
g_alpha_y = n_LhbF * np.cos(np.deg2rad(45.0))
g_alpha_z = n_LhbF * np.sin(np.deg2rad(45.0))
g_beta_y = -n_LhbR * np.cos(np.deg2rad(30.0))
g_beta_z = n_LhbR * np.sin(np.deg2rad(30.0))
guess = [g_PJ1extension, g_J1J2flexion, g_CA1extension, g_CA1adduction, g_CA1rotation, g_A1A2extension,
g_alpha_y, g_alpha_z, g_beta_y, g_beta_z]
Now we use scipy.optimize.fsolve
to hone in on the exact solution.
def my_nonlinear_function(unknowns, knowns):
"""Returns the value of the vector nonlinear function given the measured quantities.
Parameters
----------
unknowns : array_like, shape(10,)
The numerical values of the unknown quantities,
i.e. joint angles and d_ry, d_rz, d_fy, d_fz.
knowns : array_like, shape(19,)
The numerical values of the known quantities: human and
bicycle geometry.
Returns
-------
sol : ndarray, shape(10,)
"""
all_values = np.hstack((unknowns, knowns))
return np.array(numerical_function(*all_values))
solution = fsolve(my_nonlinear_function, guess, args)
We can see that the solution evaluates the functions to zero.
my_nonlinear_function(solution, args)
array([ -7.74700e-13, -1.63214e-12, 3.69027e-12, -6.43707e-12, 0.00000e+00, 5.55112e-17, -3.92657e-11, 4.28030e-12, 1.54883e-11, -1.83722e-10])
for s, v in zip(unknown_symbols, solution):
print('{} = {:1.2f}'.format(s, v))
phi_J1 = -1.16 phi_J2 = 1.19 phi_A1 = -0.62 theta_A1 = -0.16 psi_A = 0.32 phi_A2 = -0.25 alpha_y = 0.54 alpha_z = 0.65 beta_y = -0.58 beta_z = 0.65
The magnitude of the vectors from the front and rear wheel centers to the grip should be equal to the measured distances, $L_{hbf}$ and $L_{hbR}$.
print('{} == {}'.format(np.sqrt((n_whb / 2) ** 2 + solution[-3] ** 2 + solution[-4] ** 2), n_LhbF))
0.893000000004 == 0.893
print('{} == {}'.format(np.sqrt((n_whb / 2) ** 2 + solution[-1] ** 2 + solution[-2] ** 2), n_LhbR))
0.921299999994 == 0.9213
The equations can also be printed as Python code.
for i, eq in enumerate(equations):
print('zero[{}] = {}'.format(i, sym.printing.lambdarepr.lambdarepr(eq)))
print('')
zero[0] = L_j3L*(-sin(phi_J1)*cos(phi_P) - sin(phi_P)*cos(phi_J1)) + (-l_sp - l_st)*cos(lambda_st) + (-(-sin(phi_J1)*sin(phi_P) + cos(phi_J1)*cos(phi_P))*sin(phi_J2) + (-sin(phi_J1)*cos(phi_P) - sin(phi_P)*cos(phi_J1))*cos(phi_J2))*(-L_j3L + L_j5L + L_j6L) zero[1] = L_j3L*(-sin(phi_J1)*sin(phi_P) + cos(phi_J1)*cos(phi_P)) + (-l_sp - l_st)*sin(lambda_st) + ((-sin(phi_J1)*sin(phi_P) + cos(phi_J1)*cos(phi_P))*cos(phi_J2) - (sin(phi_J1)*cos(phi_P) + sin(phi_P)*cos(phi_J1))*sin(phi_J2))*(-L_j3L + L_j5L + L_j6L) zero[2] = -L_hbF + sqrt(alpha_y**2 + alpha_z**2 + 0.25*w_hb**2) zero[3] = -L_hbR + sqrt(beta_y**2 + beta_z**2 + 0.25*w_hb**2) zero[4] = alpha_y - beta_y - w zero[5] = alpha_z - beta_z + r_F - r_R zero[6] = -L_a2L*sin(theta_A1) + L_s4w/2 - 0.5*w_hb + (sin(phi_A2)*sin(psi_A)*cos(theta_A1) + sin(theta_A1)*cos(phi_A2))*(L_a2L - L_a4L - L_a5L) zero[7] = -L_a2L*(-sin(phi_A1)*cos(phi_P)*cos(theta_A1) - sin(phi_P)*cos(phi_A1)*cos(theta_A1)) - L_s4L*sin(phi_P) - beta_y - sqrt(l_cs**2 - (-h_bb + r_R)**2) - (-l_sp - l_st)*cos(lambda_st) + (-(-(sin(phi_A1)*cos(psi_A) + sin(psi_A)*sin(theta_A1)*cos(phi_A1))*sin(phi_P) + (-sin(phi_A1)*sin(psi_A)*sin(theta_A1) + cos(phi_A1)*cos(psi_A))*cos(phi_P))*sin(phi_A2) + (-sin(phi_A1)*cos(phi_P)*cos(theta_A1) - sin(phi_P)*cos(phi_A1)*cos(theta_A1))*cos(phi_A2))*(L_a2L - L_a4L - L_a5L) zero[8] = -L_a2L*(-sin(phi_A1)*sin(phi_P)*cos(theta_A1) + cos(phi_A1)*cos(phi_P)*cos(theta_A1)) + L_s4L*cos(phi_P) - beta_z + h_bb - r_R - (-l_sp - l_st)*sin(lambda_st) + (-((sin(phi_A1)*cos(psi_A) + sin(psi_A)*sin(theta_A1)*cos(phi_A1))*cos(phi_P) + (-sin(phi_A1)*sin(psi_A)*sin(theta_A1) + cos(phi_A1)*cos(psi_A))*sin(phi_P))*sin(phi_A2) + (-sin(phi_A1)*sin(phi_P)*cos(theta_A1) + cos(phi_A1)*cos(phi_P)*cos(theta_A1))*cos(phi_A2))*(L_a2L - L_a4L - L_a5L) zero[9] = (sin(phi_A1)*sin(psi_A) - sin(theta_A1)*cos(phi_A1)*cos(psi_A))*cos(phi_P) + (sin(phi_A1)*sin(theta_A1)*cos(psi_A) + sin(psi_A)*cos(phi_A1))*sin(phi_P)
We have completed our first objective! We have oriented the rider to sit on the bicycle. Now, we need to determine the rider's inertia properties? First, we tell the Yeadon human
object about the configuration in which we want the rider. The default dictionary can be extracted from the human object with:
cfg_dict = human.CFG.copy()
The left side can be set with the solution:
cfg_dict['PJ1extension'] = solution[0]
cfg_dict['J1J2flexion'] = solution[1]
cfg_dict['CA1extension'] = solution[2]
cfg_dict['CA1adduction'] = solution[3]
cfg_dict['CA1rotation'] = solution[4]
cfg_dict['A1A2extension'] = solution[5]
cfg_dict['somersault'] = n_somersault
And the right side can just be a mirror of the left.
cfg_dict['PK1extension'] = cfg_dict['PJ1extension']
cfg_dict['K1K2flexion'] = cfg_dict['J1J2flexion']
cfg_dict['CB1extension'] = cfg_dict['CA1extension']
cfg_dict['CB1abduction'] = -cfg_dict['CA1adduction']
cfg_dict['CB1rotation'] = -cfg_dict['CA1rotation']
cfg_dict['B1B2extension'] = cfg_dict['A1A2extension']
The values of the configuration are:
cfg_dict
{'A1A2extension': -0.2482939191873583, 'B1B2extension': -0.2482939191873583, 'CA1adduction': -0.15756629042407669, 'CA1extension': -0.61930340001927042, 'CA1rotation': 0.31840490580456943, 'CB1abduction': 0.15756629042407669, 'CB1extension': -0.61930340001927042, 'CB1rotation': -0.31840490580456943, 'J1J2flexion': 1.1850079456424161, 'K1K2flexion': 1.1850079456424161, 'PJ1adduction': 0.0, 'PJ1extension': -1.1645088784332249, 'PK1abduction': 0.0, 'PK1extension': -1.1645088784332249, 'PTbending': 0.0, 'PTsagittalFlexion': 0.0, 'TCsagittalSpinalFlexion': 0.0, 'TCspinalTorsion': 0.0, 'somersault': 0.175, 'tilt': 0.0, 'twist': 0.0}
The configuration can the be set with:
human.set_CFG_dict(cfg_dict)
Before visualizing the human in the correct configuration we will add the basic geometry of the bicycle for reference.
def numerical_vector(point, origin, ref_frame):
"""Returns an (3,) length array with the x, y, z coordinates of the
given point with respect to the given origin in the given reference frame."""
matrix = point.pos_from(origin).to_matrix(ref_frame)
f = sym.lambdify(all_symbols, matrix, modules='numpy')
return np.array(f(*np.hstack((solution, args)))).flatten()
numerical_vec_to_fc = numerical_vector(front_wheel_center, origin, yeadon_rf)
numerical_vec_to_rc = numerical_vector(rear_wheel_center, origin, yeadon_rf)
numerical_vec_to_left_grip = numerical_vector(grip_defined_from_front, origin, yeadon_rf)
numerical_vec_to_bottom_bracket = numerical_vector(bottom_bracket, origin, yeadon_rf)
def add_torus(outer_radius, cross_section_radius, center, orientation):
"""Adds a torus to represent the bicycle wheeels to the scene."""
torus = mlab.pipeline.parametric_surface()
torus.function = 'torus'
torus.parametric_function.ring_radius = outer_radius - cross_section_radius
torus.parametric_function.cross_section_radius = cross_section_radius
mlab.pipeline.surface(torus)
torus.children[0].children[0].actor.actor.orientation = orientation
torus.children[0].children[0].actor.actor.position = center
return torus
front_torus = add_torus(n_rF - 0.02, 0.02, numerical_vec_to_fc, np.array([0.0, 90.0, 0.0]))
rear_torus = add_torus(n_rR - 0.02, 0.02, numerical_vec_to_rc, np.array([0.0, 90.0, 0.0]))
Show some important points:
points = np.vstack((numerical_vec_to_left_grip,
numerical_vec_to_bottom_bracket,
numerical_vec_to_fc,
numerical_vec_to_rc))
points
array([[ 0.29 , -0.40641, -0.01772], [ 0. , -0.28221, -0.71642], [ 0. , -0.94551, -0.66789], [ 0. , 0.17549, -0.67046]])
null = mlab.points3d(points[:, 0], points[:, 1], points[:, 2],
0.125 * np.ones_like(points[:, 0]), scale_factor=1.0)
Draw the Yeadon model:
human.draw(mlabobj=mlab)
Take some orthographic screen shots:
mlab.view(azimuth=0.0, elevation=90.0)
mlab.get_engine().scenes[0].scene.camera.parallel_projection = True
mlab.orientation_axes()
mlab.savefig('side-view.png')
Image('side-view.png') #, width='800px')
mlab.view(azimuth=-90.0)
mlab.savefig('front-view.png')
Image('front-view.png') #, width='800px')
The inertial properties of the human in this configuration can be expressed in the Yeadon reference frame with respect to the center of mass of the human:
human.mass
83.500000000000014
human.center_of_mass
array([[ 0. ], [-0.14028], [ 0.08757]])
human.inertia
matrix([[ 11.2571 , 0. , 0. ], [ 0. , 10.98524, -1.70211], [ 0. , -1.70211, 2.26922]])
Miejaard et. al 2007 uses the SAE standard for defining the reference frame of a bicycle. The $\hat{x}_b$ unit vector is directed from the back to the front of the bicycle, the $\hat{y}_b$ is directed from the left side of the bicycle to the right side, and the $\hat{z}_b$ is directed from top to bottom. So to rotate a reference frame relative to the Yeadon base reference frame into alignment with the bicycle reference frame, you can first rotate through $\pi$ about $\hat{x}$ and then through $-\frac{\pi}{2}$ about $\hat{z}$. The direction cosine matrix that transforms a vector in blank to a vector in blank is then:
bicycle_rf = yeadon_rf.orientnew('B', 'Space', (sym.pi, -sym.pi / 2, 0), 'XZY')
bicycle_rf.dcm(yeadon_rf)
The inertia of the human with respect to the bicycle reference frame and about the human center of mass is then found by:
R = np.array([[0.0, -1.0, 0.0],
[-1.0, 0.0, 0.0],
[0.0, 0.0, -1.0]])
human.inertia_transformed(rotmat=R)
matrix([[ 10.98524, 0. , -1.70211], [ 0. , 11.2571 , 0. ], [ -1.70211, 0. , 2.26922]])
The human center of mass position with respect to the bicycle seat is expressed in the bicycle reference frame is:
R.dot(human.center_of_mass)
array([[ 0.14028], [ 0. ], [-0.08757]])
Now if we assume the human is rigidly affixed to the rear bicycle frame we can combine the inertia of the two rigid bodies into a single rigid body by using the parallel axis thereom.
First we form the vector from the bicycle origin (rear wheel contact point) to the Yeadon origin:
bicycle_origin_to_yeadon_origin = origin.pos_from(rear_wheel_center) + rear_wheel_radius * yeadon_rf.z
We can express this vector in the bicycle reference frame:
bicycle_origin_to_yeadon_origin.express(bicycle_rf)
And use this expression to compute the location given our bicycle measurements:
compute_yeadon_origin = sym.lambdify(bicycle_measurement_symbols[:7],
bicycle_origin_to_yeadon_origin.to_matrix(bicycle_rf), modules='numpy')
yeadon_origin = np.matrix(compute_yeadon_origin(*args[:7])).reshape(3, 1)
yeadon_origin
matrix([[ 0.17549], [ 0. ], [-1.01142]])
The mass center of the human can then be computed with respect to the bicycle origin and in the bicycle reference frame:
human_mass_center_in_bicycle_rf = yeadon_origin + R.dot(human.center_of_mass)
human_mass_center_in_bicycle_rf
matrix([[ 0.31577], [ 0. ], [-1.09899]])
Now we need some parameters for a the rear frame of the bicycle, inertia, center of mass, and mass. These are taken from the Browser bicycle, see http://moorepants.github.io/dissertation/physicalparameters.html.
n_IBxx = 0.52962890621
n_IBxz = -0.116285607878
n_IByy = 1.3163960125
n_IBzz = 0.756786895402
n_xB = 0.275951285677
n_zB = -0.537842424305
n_mB = 9.86
rear_frame_mass_center = np.matrix([[n_xB],
[0.0],
[n_zB]])
rear_frame_mass_center
matrix([[ 0.27595], [ 0. ], [-0.53784]])
rear_frame_inertia_tensor = np.matrix([[n_IBxx, 0.0, n_IBxz],
[0.0, n_IByy, 0.0],
[n_IBxz, 0.0, n_IBzz]])
rear_frame_inertia_tensor
matrix([[ 0.52963, 0. , -0.11629], [ 0. , 1.3164 , 0. ], [-0.11629, 0. , 0.75679]])
Total mass:
total_mass = n_mB + human.mass
total_mass
93.360000000000014
The total center of mass of the bicycle rear frame and the rider is:
tmp = human.mass * human_mass_center_in_bicycle_rf
total_center_of_mass = np.matrix([[n_mB * n_xB + tmp[0, 0]],
[n_mB * 0.0 + tmp[1, 0]],
[n_mB * n_zB + tmp[2, 0]]]) / total_mass
total_center_of_mass
matrix([[ 0.31156], [ 0. ], [-1.03972]])
Now use the parallel axis thereom to compute the inertias about the new center of mass and sum them.
tmp = (total_center_of_mass - rear_frame_mass_center)
a = tmp[0, 0]
b = tmp[1, 0]
c = tmp[2, 0]
rear_frame_inertia_about_total_CoM = \
rear_frame_inertia_tensor + n_mB * np.matrix([[b**2 + c**2, -a * b, -a * c],
[-a * b, c**2 + a**2, -b * c],
[-a * c, -b * c, a**2 + b**2]])
rear_frame_inertia_about_total_CoM
matrix([[ 3.01322, 0. , 0.05994], [ 0. , 3.81249, 0. ], [ 0.05994, 0. , 0.76929]])
tmp = (total_center_of_mass - human_mass_center_in_bicycle_rf)
a = tmp[0, 0]
b = tmp[1, 0]
c = tmp[2, 0]
human_inertia_about_total_CoM = \
human.inertia_transformed(rotmat=R) + human.mass * np.matrix([[b**2 + c**2, -a * b, -a * c],
[-a * b, c**2 + a**2, -b * c],
[-a * c, -b * c, a**2 + b**2]])
human_inertia_about_total_CoM
matrix([[ 11.27852, 0. , -1.6813 ], [ 0. , 11.55185, 0. ], [ -1.6813 , 0. , 2.2707 ]])
Finally, we have the combined inertia of the rider and the bicycle rear frame expressed in the bicycle frame about the combined center of mass.
total_inertia = rear_frame_inertia_about_total_CoM + human_inertia_about_total_CoM
total_inertia
matrix([[ 14.29174, 0. , -1.62136], [ 0. , 15.36434, 0. ], [ -1.62136, 0. , 3.03999]])
%install_ext version_information.py
Installed version_information.py. To use it, type: %load_ext version_information
%load_ext version_information
%version_information numpy, scipy, sympy, mayavi, yeadon
Software | Version |
---|---|
Python | 2.7.5+ (default, Feb 27 2014, 19:37:08) [GCC 4.8.1] |
IPython | 2.0.0-dev |
OS | posix [linux2] |
numpy | 1.7.1 |
scipy | 0.12.0 |
sympy | 0.7.5-git |
mayavi | 4.1.0 |
yeadon | 1.2.0 |
Wed Sep 03 15:43:06 2014 PDT |