Modeling of a Variable-Mass Nonholonomic Gyrostatic Rocket Car Using Extended Kane’s Equations¶
Note
You can download this example as a Python script:
rocket-car.py
or Jupyter notebook:
rocket-car.ipynb
.
This is an implementation of the nonholonomic rocket-engine-powered jet racing car example from [Ge1982]. This example provides insight into analytically modeling variable-mass systems using Kane’s method and the extended Kane’s equations for variable-mass systems. It also demonstrates the efficacy of Kane’s method in the modeling of complex dynamical systems.
An idealized model of a jet racing car that is propelled by a rocket engine at point P is considered. The rocket engine is treated as a variable-mass particle at P.
Here, \(\{a_x: g_3, a_y: g_1, a_z: g_2\}\)
import sympy as sm
import sympy.physics.mechanics as me
from pydy.system import System
import numpy as np
from sympy.simplify.fu import TR2
import matplotlib.pyplot as plt
from scipy.integrate import odeint
me.init_vprinting()
from IPython.display import Latex, display
def print_answer(x, x_res):
for xi, xri in zip(x, x_res):
display(Latex(sm.latex(sm.Eq(xi, xri), mode='inline')))
Reference Frames, Generalized Coordinates, and Generalized Speeds¶
N = me.ReferenceFrame('N')
q1, q2, q3, q4, q5, q6, q7, q8 = me.dynamicsymbols('q1:9')
A2 = N.orientnew('A_2', 'Axis', [q3, N.z])
A1 = A2.orientnew('A_1', 'Axis', [q4, A2.z])
B1 = A1.orientnew('B_1', 'Axis', [q5, A1.y])
B2 = A1.orientnew('B_2', 'Axis', [q6, A1.y])
B3 = A2.orientnew('B_3', 'Axis', [q7, A2.y])
B4 = A2.orientnew('B_4', 'Axis', [q8, A2.y])
t = me.dynamicsymbols._t
O = me.Point('O') # fixed point in the inertial reference frame
O.set_vel(N, 0)
L, l , a, b, r1, r2 = sm.symbols('L, l , a, b, r_1, r_2')
Q = O.locatenew('Q', q1 * N.x + q2 * N.y)
P = Q.locatenew('P', L * -A2.x)
C = P.locatenew('C', l * A2.x)
Q.set_vel(N, Q.pos_from(O).dt(N))
Q.vel(N)
P.v2pt_theory(Q, N, A2)
P.vel(N)
C.v2pt_theory(P, N, A2)
# C.vel(N)
A1.ang_vel_in(A2).express(A1)
u1, u2 = me.dynamicsymbols('u_1:3')
z1 = sm.Eq(u1, A1.ang_vel_in(A2).dot(A1.z))
z2 = sm.Eq(u2, Q.vel(N).dot(A1.x))
u = sm.trigsimp(sm.solve([z1, z2], A1.ang_vel_in(A2).dot(A1.z), Q.vel(N).dot(A1.x)))
u
Formulation of the Constraint Equations¶
Nonholonomic Constraints: \(B_1\)¶
B1_center = Q.locatenew('B_1_center', a * A1.y)
B1_center.pos_from(Q)
B1_center.v2pt_theory(Q, N, A1)
B1_center.vel(N).express(A1).simplify()
B1_ground = B1_center.locatenew('B_1_ground', r1 * -A1.z)
B1_ground.pos_from(B1_center)
B1_ground.v2pt_theory(B1_center, N, B1)
B1_ground.vel(N).simplify()
B1_cons = [me.dot(B1_ground.vel(N).simplify(), uv) for uv in A1]
for i in range(len(B1_cons)):
display(sm.trigsimp(B1_cons[i]))
eq1 = sm.Eq(B1_cons[0].simplify().subs(u), 0)
eq1
eq2 = sm.Eq(B1_cons[1].simplify().subs(u), 0)
eq2
Nonholonomic Constraints: \(B_2\)¶
B2_center = Q.locatenew('B_1_center', a * -A1.y)
B2_center.pos_from(Q)
B2_center.v2pt_theory(Q, N, A1)
B2_center.vel(N).express(A1).simplify()
B2_ground = B2_center.locatenew('B_2_ground', r1 * -A1.z)
B2_ground.pos_from(B2_center)
B2_ground.v2pt_theory(B2_center, N, B2)
B2_ground.vel(N).simplify()
B2_cons = [me.dot(B2_ground.vel(N).simplify(), uv) for uv in A1]
for i in range(len(B2_cons)):
display(sm.trigsimp(B2_cons[i]))
eq3 = sm.Eq(B2_cons[0].simplify().subs(u), 0)
eq3
eq4 = sm.Eq(B2_cons[1].simplify().subs(u), 0)
eq4
Nonholonomic Constraints: \(B_3\)¶
B3_center = P.locatenew('B_3_center', b * A2.y)
B3_center.pos_from(P)
B3_center.v2pt_theory(P, N, A2)
B3_center.vel(N).express(A2).simplify()
B3_ground = B3_center.locatenew('B_3_ground', r2 * -A2.z)
B3_ground.pos_from(B3_center)
B3_ground.v2pt_theory(B3_center, N, B3)
B3_ground.vel(N).simplify()
B3_cons = [me.dot(B3_ground.vel(N).simplify(), uv) for uv in A2]
for i in range(len(B3_cons)):
display(sm.trigsimp(B3_cons[i]))
eq5 = sm.Eq(B3_cons[0].simplify().subs(u), 0)
eq5
eq6 = sm.Eq(B3_cons[1].simplify().subs(u), 0)
eq6
Nonholonomic Constraints: \(B_4\)¶
B4_center = P.locatenew('B_4_center', b * -A2.y)
B4_center.pos_from(P)
B4_center.v2pt_theory(P, N, A2)
B4_center.vel(N).express(A2).simplify()
B4_ground = B4_center.locatenew('B_4_ground', r2 * -A2.z)
B4_ground.pos_from(B4_center)
B4_ground.v2pt_theory(B4_center, N, B4)
B4_ground.vel(N).simplify()
B4_cons = [me.dot(B4_ground.vel(N).simplify(), uv) for uv in A2]
for i in range(len(B4_cons)):
display(sm.trigsimp(B4_cons[i]))
eq7 = sm.Eq(B4_cons[0].simplify().subs(u), 0)
eq7
eq8 = sm.Eq(B4_cons[1].simplify().subs(u), 0)
eq8
\(\text{LHS} \Longleftrightarrow \text{RHS}\ \text{in}\ z_1, z_2 \rightarrow \text{Equations}\ 9, 10\)
eq9 = sm.Eq(A1.ang_vel_in(A2).dot(A1.z), u1)
eq9
eq10 = sm.Eq(Q.vel(N).dot(A1.x), u2)
eq10
Solving the System of Linear Equations¶
The system of equations is linear in \(\dot{q}_1, \dot{q}_2,...\)
Note: eq4
\(\equiv\) eq2
; eq8
\(\equiv\) eq6
solution = sm.linsolve([eq1, eq2, eq3, eq5, eq6, eq7, eq9, eq10], q1.diff(), q2.diff(), q3.diff(), q4.diff(), q5.diff(), q6.diff(), q7.diff(), q8.diff())
sollist_keys = [q1.diff(), q2.diff(), q3.diff(), q4.diff(), q5.diff(), q6.diff(), q7.diff(), q8.diff()]
sollist_keys
sollist_values = list(solution.args[0])
sollist_values_simple = []
for i in range(len(sollist_values)):
sollist_values_simple.append(sm.factor(TR2(sollist_values[i]).simplify()))
soldict = dict(zip(sollist_keys, sollist_values_simple))
print_answer(sollist_keys, sollist_values_simple)
Reformulated Velocity and Angular Velocity Expressions¶
N_v_Q = Q.vel(N).subs(soldict).express(A1).simplify()
N_v_Q
N_v_P = P.vel(N).subs(soldict).express(A2).simplify()
N_v_P
N_v_C = C.vel(N).subs(soldict).express(A2).simplify()
N_v_C
N_w_A1 = A1.ang_vel_in(N).subs(soldict).express(A1).simplify()
N_w_A1
N_w_A2 = A2.ang_vel_in(N).subs(soldict).express(A2).simplify()
N_w_A2
Partial Velocities and Partial Angular Velocities¶
V_1_Q = N_v_Q.diff(u1, N)
V_1_Q
V_2_Q = N_v_Q.diff(u2, N)
V_2_Q
V_1_C = N_v_C.diff(u1, N)
V_1_C
V_2_C = N_v_C.diff(u2, N)
V_2_C
V_1_P = N_v_P.diff(u1, N)
V_1_P
V_2_P = N_v_P.diff(u2, N)
V_2_P
w_1_A1 = N_w_A1.diff(u1, N)
w_1_A1
w_2_A1 = N_w_A1.diff(u2, N)
w_2_A1
w_1_A2 = N_w_A2.diff(u1, N)
w_1_A2
w_2_A2 = N_w_A2.diff(u2, N)
w_2_A2
Accelerations and Angular Accelerations¶
a_1__P, a_2__P, a_3__P, a_1__C, a_2__C, a_3__C, a__Q, alpha__A1, alpha__A2 = sm.symbols('a_1__P, a_2__P, a_3__P, a_1__C, a_2__C, a_3__C, a__Q, alpha__A1, alpha__A2')
N_a_P = N_v_P.dt(N).subs(soldict)
N_a_P
N_a_C = N_v_C.dt(N).subs(soldict)
N_a_C
N_a_Q = N_v_Q.dt(N).subs(soldict)
N_a_Q
N_aa_A1 = N_w_A1.dt(N).subs(soldict)
N_aa_A1
N_aa_A2 = N_w_A2.dt(N).subs(soldict)
N_aa_A2
Forces and Torques¶
\((F_r^*)_G = (F_r^*)_{GR} + (F_r^*)_{GI}\)
where,
\((F_r^*)_{GR} = {V_r}^G \cdot {F_G}^* + \omega_r^A \cdot {T_G}^*\)
\(F_G^* = -m_G {a^G}^*\)
\(T_G^* \overset{\Delta}{=} -[\alpha_A \cdot I_G + \omega_r^A \times (I_G \cdot \omega_r^A)]\)
\(({F_r}^*)_{GI} = -J\{\omega_r^A [\ddot{q_k} g_1 + \dot{q_k} (\omega_3^A g_2 - \omega_2^A g_3)] + C_{kr} (\dot{\omega_1^A} + \ddot{q_k}) \}\)
Naming Convention:
\(({F_r}^*)_{G_n R}\) (rigid)
\(({F_r}^*)_{G_n I}\) (internal)
Masses and Moments of Inertia¶
M1, M2 = sm.symbols('M_1, M_2')
m = me.dynamicsymbols('m')
I1x, I1y, I1z = sm.symbols('I_{1_x}, I_{1_y}, I_{1_z}')
I2x, I2y, I2z = sm.symbols('I_{2_x}, I_{2_y}, I_{2_z}')
J1, J2 = sm.symbols('J_1, J_2')
I1 = me.inertia(A1, I1x, I1y, I1z)
I1
I2 = me.inertia(A2, I2x, I2y, I2z)
I2
Gyrostat \(G_1\)¶
\(\rightarrow {F_G}^* = -m_G {a^G}^*\)
Fstar_G1 = -M1 * N_a_Q
Fstar_G1
\(\rightarrow {T_G}^* \overset{\Delta}{=} -[\alpha_A \cdot I_G + {\omega_r}^A \times (I_G \cdot {\omega_r}^A)]\)
Tstar_G1 = -(N_aa_A1.dot(I1) + me.cross(N_w_A1, I1.dot(N_w_A1)))
Tstar_G1
\(\rightarrow ({F_r}^*)_{GR} = {V_r}^G \cdot {F_G}^* + {\omega_r}^A \cdot {T_G}^*\)
Fstar_1_G1_R = V_1_Q.dot(Fstar_G1) + w_1_A1.dot(Tstar_G1).subs(soldict)
Fstar_1_G1_R.subs({N_w_A1.dt(N).subs(soldict).dot(A1.z): alpha__A1})
Fstar_2_G1_R = V_2_Q.dot(Fstar_G1) + w_2_A1.dot(Tstar_G1).subs(soldict)
Fstar_2_G1_R.subs({N_w_A1.dt(N).subs(soldict).dot(A1.z): alpha__A1})
\(\rightarrow (F_r^*)_{GI} = -J\{\omega_r^A \cdot [\ddot{q_k} g_1 + \dot{q_k} (\omega_3^A g_2 - \omega_2^A g_3)] + C_{kr} (\dot{\omega}_1^A + \ddot{q_k}) \} \qquad (r=1,...,n-m)\)
Here, \(\{\omega_1^A: \omega_2^A,\ \omega_2^A: \omega_3^A,\ \omega_3^A: \omega_1^A\}\)
\(\rightarrow \dot{q_k} = \sum_{s = 1}^{n - m} C_{ks} u_s + D_k\) (Generalized Speeds)
\(\omega_i^A \overset{\Delta}{=} \omega^A \cdot \hat{g}_i \quad (i = 1, 2, 3)\)
# C_kr
C51, C61 = sm.symbols('C_51, C_61')
C_51 = soldict[q5.diff()].diff(u1)
C_61 = soldict[q6.diff()].diff(u1)
Fstar_1_G1_I = -J1 * (N_w_A1.dot(q5.diff().diff() * A1.y + q5.diff()*(N_w_A1.dot(A1.x)*A1.z - N_w_A1.dot(A1.z)*A1.x)) + C_51 * (N_w_A1.dot(A1.y).diff() + q5.diff().diff())) \
-J1 * (N_w_A1.dot(q6.diff().diff() * A1.y + q6.diff()*(N_w_A1.dot(A1.x)*A1.z - N_w_A1.dot(A1.z)*A1.x)) + C_61 * (N_w_A1.dot(A1.y).diff() + q6.diff().diff())) # B1 \ B2
Fstar_1_G1_I, C_51, C_61, Fstar_1_G1_I.subs({-C_51: -C51, -C_61: -C61}).simplify()
# C_kr
C52, C62 = sm.symbols('C_52, C_62')
C_52 = soldict[q5.diff()].diff(u2)
C_62 = soldict[q6.diff()].diff(u2)
Fstar_2_G1_I = -J1 * (N_w_A1.dot(q5.diff().diff() * A1.y + q5.diff()*(N_w_A1.dot(A1.x)*A1.z - N_w_A1.dot(A1.z)*A1.x)) + C_52 * (N_w_A1.dot(A1.y).diff() + q5.diff().diff())) \
-J1 * (N_w_A1.dot(q6.diff().diff() * A1.y + q6.diff()*(N_w_A1.dot(A1.x)*A1.z - N_w_A1.dot(A1.z)*A1.x)) + C_62 * (N_w_A1.dot(A1.y).diff() + q6.diff().diff())) # B1 \ B2
display(Fstar_2_G1_I),
display(C_52)
display(C_62)
display(Fstar_2_G1_I.subs({-C_52: -C52, -C_62: -C62}).simplify())
\(\rightarrow (F_r^*)_G = (F_r^*)_{GR} + (F_r^*)_{GI}\)
Fstar_1_G1 = Fstar_1_G1_R + Fstar_1_G1_I
Fstar_1_G1.subs({N_w_A1.dt(N).subs(soldict).dot(A1.z): alpha__A1}).subs({-C_51: -C51, -C_61: -C61}).simplify()
Fstar_2_G1 = Fstar_2_G1_R + Fstar_2_G1_I
Fstar_2_G1.subs({N_w_A1.dt(N).subs(soldict).dot(A1.z): alpha__A1}).subs({-C_52: -C52, -C_62: -C62}).simplify()
Gyrostat \(G_2\)¶
\(\rightarrow {F_G}^* = -m_G {a^G}^*\)
Fstar_G2 = -M2 * N_a_C
Fstar_G2
\(\rightarrow {T_G}^* \overset{\Delta}{=} -[\alpha_A \cdot I_G + {\omega_r}^A \times (I_G \cdot {\omega_r}^A)]\)
Tstar_G2 = -(N_aa_A2.dot(I2) + me.cross(N_w_A2, I2.dot(N_w_A2)))
Tstar_G2
\(\rightarrow ({F_r}^*)_{GR} = {V_r}^G \cdot {F_G}^* + {\omega_r}^A \cdot {T_G}^*\)
Fstar_1_G2_R = V_1_C.dot(Fstar_G2) + w_1_A2.dot(Tstar_G2).subs(soldict)
Fstar_1_G2_R.subs({N_w_A2.dt(N).subs(soldict).dot(A2.z): alpha__A2})
Fstar_2_G2_R = V_2_C.dot(Fstar_G2) + w_2_A1.dot(Tstar_G2).subs(soldict)
Fstar_2_G2_R.subs({N_w_A2.dt(N).subs(soldict).dot(A2.z): alpha__A2})
\(\rightarrow (F_r^*)_{GI} = -J\{\omega_r^A \cdot [\ddot{q_k} g_1 + \dot{q_k} (\omega_3^A g_2 - \omega_2^A g_3)] + C_{kr} (\dot{\omega}_1^A + \ddot{q_k}) \} \qquad (r=1,...,n-m)\)
Here, \(\{\omega_1^A: \omega_2^A,\ \omega_2^A: \omega_3^A,\ \omega_3^A: \omega_1^A\}\)
\(\rightarrow \dot{q_k} = \sum_{s = 1}^{n - m} C_{ks} u_s + D_k\) (Generalized Speeds)
\(\omega_i^A \overset{\Delta}{=} \omega^A \cdot \hat{g}_i \quad (i = 1, 2, 3)\)
# C_kr
C71, C81 = sm.symbols('C_71, C_81')
C_71 = soldict[q7.diff()].diff(u1)
C_81 = soldict[q8.diff()].diff(u1)
Fstar_1_G2_I = -J2 * (N_w_A2.dot(q7.diff().diff() * A2.y + q7.diff()*(N_w_A2.dot(A2.x)*A2.z - N_w_A2.dot(A2.z)*A2.x)) + C_71 * (N_w_A2.dot(A2.y).diff() + q7.diff().diff())) \
-J2 * (N_w_A2.dot(q8.diff().diff() * A2.y + q8.diff()*(N_w_A2.dot(A2.x)*A2.z - N_w_A2.dot(A2.z)*A2.x)) + C_81 * (N_w_A2.dot(A2.y).diff() + q8.diff().diff())) # B1 \ B2
Fstar_1_G2_I, C_71, C_81, # Fstar_1_G2_I.subs({-C_71: -C71, -C_81: -C81}).simplify()
# C_kr
C72, C82 = sm.symbols('C_72, C_82')
C_72 = soldict[q7.diff()].diff(u2)
C_82 = soldict[q8.diff()].diff(u2)
Fstar_2_G2_I = -J2 * (N_w_A2.dot(q7.diff().diff() * A2.y + q7.diff()*(N_w_A2.dot(A2.x)*A2.z - N_w_A2.dot(A2.z)*A2.x)) + C_72 * (N_w_A2.dot(A2.y).diff() + q7.diff().diff())) \
-J2 * (N_w_A2.dot(q8.diff().diff() * A2.y + q8.diff()*(N_w_A2.dot(A2.x)*A2.z - N_w_A2.dot(A2.z)*A2.x)) + C_82 * (N_w_A2.dot(A2.y).diff() + q8.diff().diff())) # B1 \ B2
display(Fstar_2_G2_I)
display(C_72)
display(C_82)
display(Fstar_2_G2_I.subs({-C_72: -C72, -C_82: -C82}).simplify())
\(\rightarrow (F_r^*)_G = (F_r^*)_{GR} + (F_r^*)_{GI}\)
Fstar_1_G2 = Fstar_1_G2_R + Fstar_1_G2_I
# Fstar_1_G2.subs({N_w_A2.dt(N).subs(soldict).dot(A2.z): alpha__A2}) # .subs({-C_71: -C71, -C_81: -C81}).simplify()
Fstar_1_G2 = 0
Here, \(\{a_1^C: a_2^C,\ a_2^C: a_3^C,\ a_3^C: a_1^C\}\)
Fstar_2_G2 = Fstar_2_G2_R + Fstar_2_G2_I
Fstar_2_G2.subs({N_w_A2.dt(N).subs(soldict).dot(A2.z): alpha__A2}).subs({N_v_C.dt(N).subs(soldict).dot(A2.x): a_3__C}).subs({N_v_C.dt(N).subs(soldict).dot(A2.y): a_1__C}).subs({-C_72: -C72, -C_82: -C82}).simplify()
Variable-Mass Particle, \(P\)¶
\(\rightarrow {F_G}^* = -m_G {a^G}^*\)
Fstar_P = -m * N_a_P
Fstar_P
\(\rightarrow ({F_r}^*)_{GR} = {V_r}^G \cdot {F_G}^*\)
Fstar_1_P_R = V_1_P.dot(Fstar_P)
Fstar_1_P_R
Fstar_2_P_R = V_2_P.dot(Fstar_P)
Fstar_2_P_R
\(\rightarrow (F_r^*)_G = (F_r^*)_{GR}\)
Fstar_1_P = Fstar_1_P_R
Fstar_1_P
Here, \(\{a_1^P: a_2^P,\ a_2^P: a_3^P,\ a_3^P: a_1^P\}\)
Fstar_2_P = Fstar_2_P_R
Fstar_2_P.subs({N_v_P.dt(N).subs(soldict).dot(A2.x): a_3__P}).subs({N_v_P.dt(N).subs(soldict).dot(A2.y): a_1__P}).simplify()
Generalized Inertia Forces¶
\(\rightarrow F_r^* = (F_r^*)_{G_1} + (F_r^*)_{G_2} + (F_r^*)_{P} \quad (r = 1, 2)\)
Fstar_1 = Fstar_1_G1 + Fstar_1_G2 + Fstar_1_P
Fstar_1.subs(soldict).simplify()
Fstar_2 = Fstar_2_G1 + Fstar_2_G2 + Fstar_2_P
Fstar_2.subs(soldict).simplify()
Velocity of material ejected at \(P\) relative to \(A_2 \rightarrow -C(t)g_3^{'}\)
\(C(t) \rightarrow\) positive
C = me.dynamicsymbols('C')
C_t = -C*A2.x
C_t
Generalized Thrust¶
\(\rightarrow F_r^{\prime} \triangleq \sum_{i=1}^{N} \mathbf{V}_{r}^{P i} \cdot \mathbf{C}^{P i} \dot{m}_{i} \quad (r=1, \ldots, k)\)
Fprime_1 = V_1_P.dot(C_t)*m.diff()
Fprime_1
Fprime_2 = V_2_P.dot(C_t)*m.diff()
Fprime_2
Extended Kane’s Equations for Variable-Mass Systems¶
\(\rightarrow F_r + F_r^* + F_r^{\prime} = 0 \quad (r = 1,..., k)\)
Here, \(F_r = 0 \rightarrow\) no forces contributing to generalized active forces
kane_1 = Fstar_1.simplify() + Fprime_1.simplify()
kane_1.subs(soldict).simplify()
kane_2 = Fstar_2 + Fprime_2
kane_2.subs(soldict).simplify()
kane_1_eq = sm.Eq(kane_1.simplify().subs(soldict).simplify().subs(u).simplify(), 0)
kane_1_eq
kane_2_eq = sm.Eq(kane_2.simplify().subs(soldict).simplify().subs(u).simplify(), 0)
kane_2_eq
References¶
- Ge1982
Ge, Z., and Cheng, Y. (June 1, 1982). “Extended Kane’s Equations for Nonholonomic Variable Mass System.” ASME. J. Appl. Mech. June 1982; 49(2): 429–431. https://doi.org/10.1115/1.3162105
- Kane1978
Kane, T.R., 1978. Nonholonomic multibody systems containing gyrostats. In Dynamics of Multibody Systems (pp. 97-107). Springer, Berlin, Heidelberg.