How to calculate the body fixed angular velocity of the head and body fixed acceleration of the head center

The issue includes a SymPy and symbrim notebook example but names no repository files or tests. Start by reproducing the shown angular-velocity extraction and clarify whether the intended deliverable is documentation or an implementation for acceleration. Done should specify the requested calculation and where it should be documented or tested.

Here is the most basic bicycle model from the symbrim tutorial:

import sympy as sm
import sympy.physics.mechanics as me

import symbrim as sb

bicycle_def = sb.WhippleBicycle("my_bike")
bicycle_def.rear_frame = sb.RigidRearFrame.from_convention("moore",
                                                           "rear_frame")
bicycle_def.rear_wheel = sb.KnifeEdgeWheel("rear_wheel")
bicycle_def.rear_tire = sb.NonHolonomicTire("rear_tire")
bicycle_def.ground = sb.FlatGround("ground")
bicycle_def.front_frame = sb.RigidFrontFrame("front_frame")
bicycle_def.front_wheel = sb.KnifeEdgeWheel("front_wheel")
bicycle_def.front_tire = sb.NonHolonomicTire("front_tire")

bicycle_def.define_connections()
bicycle_def.define_objects()

# Change the symbol names of the rear and front wheel radius
bicycle_def.rear_wheel.symbols["r"] = sm.Symbol("Rr")
bicycle_def.front_wheel.symbols["r"] = sm.Symbol("Rf")
bicycle_def.define_kinematics()
bicycle_def.define_loads()
bicycle_def.define_constraints()

system_def = bicycle_def.to_system()

normal = bicycle_def.ground.get_normal(bicycle_def.ground.origin)
g = sm.symbols("g")
disturbance = me.dynamicsymbols("disturbance")
steer_torque = me.dynamicsymbols("steer_torque")
system_def.apply_uniform_gravity(-g*normal)
system_def.add_loads(
    me.Force(bicycle_def.rear_frame.saddle.point,
             disturbance*bicycle_def.rear_frame.wheel_hub.axis)
)
system_def.add_actuators(
    me.TorqueActuator(steer_torque,
                      bicycle_def.rear_frame.steer_hub.axis,
                      bicycle_def.rear_frame.steer_hub.frame,
                      bicycle_def.front_frame.steer_hub.frame)
)
system_def.loads, system_def.actuators

system_def.q_ind = [*bicycle_def.q[:4], *bicycle_def.q[5:]]
system_def.q_dep = [bicycle_def.q[4]]
system_def.u_ind = [bicycle_def.u[3], *bicycle_def.u[5:7]]
system_def.u_dep = [*bicycle_def.u[:3], bicycle_def.u[4], bicycle_def.u[7]]

system_def.validate_system()

eoms = system_def.form_eoms(constraint_solver="CRAMER")

In this model, the torso and head are considered rigidly fixed to the rear frame so we can extract the angular velocity of the rear frame as the most basic approximation of an IMU on the rider's head.

The body fixed angular velocity of the rear frame observed from the ground can be extracted like this:

In [20]: rear_frame_rigid_body = system_def.bodies[7]

In [21]: rear_frame_rigid_body
Out[21]: RigidBody('rear_frame_body', masscenter=rear_frame_body_masscenter, frame=rear_frame_body_frame, mass=rear_frame_body_mass, inertia=Inertia(dyadic=rear_frame_ixx*(rear_frame_body_frame.x|rear_frame_body_frame.x) + rear_frame_izx*(rear_frame_body_frame.x|rear_frame_body_frame.z) + rear_frame_iyy*(rear_frame_body_frame.y|rear_frame_body_frame.y) + rear_frame_izx*(rear_frame_body_frame.z|rear_frame_body_frame.x) + rear_frame_izz*(rear_frame_body_frame.z|rear_frame_body_frame.z), point=rear_frame_body_masscenter))

In [22]: rear_frame_rigid_body.frame.ang_vel_in(system_def.frame).express(rear_frame_rigid_body.frame)
Out[22]: (-my_bike_u3(t)*sin(my_bike_q5(t))*cos(my_bike_q4(t)) + my_bike_u4(t)*cos(my_bike_q5(t)))*rear_frame_body_frame.x + (my_bike_u3(t)*sin(my_bike_q4(t)) + my_bike_u5(t))*rear_frame_body_frame.y + (my_bike_u3(t)*cos(my_bike_q4(t))*cos(my_bike_q5(t)) + my_bike_u4(t)*sin(my_bike_q5(t)))*rear_frame_body_frame.z

If you want a nice function to evaluate numerically, then you can do something like this:

In [24]: ang_vel_mat = rear_frame_rigid_body.frame.ang_vel_in(system_def.frame).to_matrix(rear_frame_rigid_body.frame)

In [25]: ang_vel_mat
Out[25]: 
Matrix([
[-my_bike_u3(t)*sin(my_bike_q5(t))*cos(my_bike_q4(t)) + my_bike_u4(t)*cos(my_bike_q5(t))],
[                                       my_bike_u3(t)*sin(my_bike_q4(t)) + my_bike_u5(t)],
[ my_bike_u3(t)*cos(my_bike_q4(t))*cos(my_bike_q5(t)) + my_bike_u4(t)*sin(my_bike_q5(t))]])

In [26]: me.find_dynamicsymbols(ang_vel_mat)
Out[26]: {my_bike_q4(t), my_bike_q5(t), my_bike_u3(t), my_bike_u4(t), my_bike_u5(t)}

In [27]: ang_vel_mat.free_symbols
Out[27]: {t}

In [28]: system_def.q
Out[28]: 
Matrix([
[my_bike_q1(t)],
[my_bike_q2(t)],
[my_bike_q3(t)],
[my_bike_q4(t)],
[my_bike_q6(t)],
[my_bike_q7(t)],
[my_bike_q8(t)],
[my_bike_q5(t)]])

In [29]: system_def.u
Out[29]: 
Matrix([
[my_bike_u4(t)],
[my_bike_u6(t)],
[my_bike_u7(t)],
[my_bike_u1(t)],
[my_bike_u2(t)],
[my_bike_u3(t)],
[my_bike_u5(t)],
[my_bike_u8(t)]])

In [31]: eval_ang_vel = sm.lambdify((system_def.q, system_def.u), ang_vel_mat[:])

In [32]: import numpy as np

In [33]: q_vals, u_vals = np.random.random(8), np.random.random(8)

In [34]: eval_ang_vel(q_vals, u_vals)
Out[34]: 
[np.float64(0.3952612865917701),
 np.float64(0.6600391860245965),
 np.float64(0.46399481716245594)]

This goes for any reference frame extracted from any rigid body in a bicycle model.