Viscous Nutation Damper



In [9]:

    
from miscpy.utils.sympyhelpers import *
init_printing()
w1,w2,w3,w1d,w2d,w3d,D,O,Od,I1,Iw,Is,It,M1,M2 = \
symbols('omega_1,omega_2,omega_3,omegadot_1,omegadot_2,\
omegadot_3,D,Omega,Omegadot,I_1,I_W,I_S,I_T,M1,M2')
diffmap = {w1:w1d,w2:w2d,w3:w3d,O:Od}

Define total and wheel MOI in $\mathcal B$ frame (about respective centers of mass):



In [10]:

    
Isat_B = diag(It,Is,It)
Iwheel_B = diag(I1,I1,Iw)
Isat_B,Iwheel_B









    Out[10]:





$$\left ( \left[\begin{matrix}I_{T} & 0 & 0\\0 & I_{S} & 0\\0 & 0 & I_{T}\end{matrix}\right], \quad \left[\begin{matrix}I_{1} & 0 & 0\\0 & I_{1} & 0\\0 & 0 & I_{W}\end{matrix}\right]\right )$$

Define ${}^\mathcal{I}\boldsymbol{\omega}^\mathcal{B}$ and ${}^\mathcal{B}\boldsymbol{\omega}^\mathcal{A}$ where $\mathcal A$ is the frame fixed to the wheel ($\mathbf b_3 \equiv \mathbf a_3$):



In [19]:

    
iWb = Matrix([w1,w2,w3])
bWa = Matrix([0,0,O])
iWb,bWc









    Out[19]:





$$\left ( \left[\begin{matrix}\omega_{1}\\\omega_{2}\\\omega_{3}\end{matrix}\right], \quad \left[\begin{matrix}0\\0\\\Omega\end{matrix}\right]\right )$$

Calculate total angular momentum:



In [20]:

    
hG_sat = Isat_B*iWb + Iwheel_B*bWa
hG_sat









    Out[20]:





$$\left[\begin{matrix}I_{T} \omega_{1}\\I_{S} \omega_{2}\\I_{T} \omega_{3} + I_{W} \Omega\end{matrix}\right]$$

Angular momentum of the wheel:



In [21]:

    
hG_wheel = Iwheel_B*(bWa+iWb); hG_wheel









    Out[21]:





$$\left[\begin{matrix}I_{1} \omega_{1}\\I_{1} \omega_{2}\\I_{W} \left(\Omega + \omega_{3}\right)\end{matrix}\right]$$

Set up and solve systems of differential equations assuming no external torques and $-D\Omega$ torque internal torque about the wheel spin axis:



In [22]:

    
eq1 = difftotalmat(hG_sat,t,diffmap) + skew(iWb)*hG_sat;eq1









    Out[22]:





$$\left[\begin{matrix}- I_{S} \omega_{2} \omega_{3} + I_{T} \dot{\omega}_{1} + \omega_{2} \left(I_{T} \omega_{3} + I_{W} \Omega\right)\\I_{S} \dot{\omega}_{2} + I_{T} \omega_{1} \omega_{3} - \omega_{1} \left(I_{T} \omega_{3} + I_{W} \Omega\right)\\I_{S} \omega_{1} \omega_{2} - I_{T} \omega_{1} \omega_{2} + I_{T} \dot{\omega}_{3} + I_{W} \dot{\Omega}\end{matrix}\right]$$



In [23]:

    
eq2 = difftotalmat(hG_wheel,t,diffmap) + skew(iWb)*hG_wheel - Matrix([M1,M2,-D*O]);eq2









    Out[23]:





$$\left[\begin{matrix}- I_{1} \omega_{2} \omega_{3} + I_{1} \dot{\omega}_{1} + I_{W} \omega_{2} \left(\Omega + \omega_{3}\right) - M_{1}\\I_{1} \omega_{1} \omega_{3} + I_{1} \dot{\omega}_{2} - I_{W} \omega_{1} \left(\Omega + \omega_{3}\right) - M_{2}\\D \Omega + I_{W} \left(\dot{\Omega} + \dot{\omega}_{3}\right)\end{matrix}\right]$$



In [37]:

    
sol1 = simplify(solve((eq1,eq2[2]),(w1d,w2d,w3d,Od)))
sol1









    Out[37]:





$$\left \{ \dot{\Omega} : \frac{- D I_{T} \Omega + I_{W} \omega_{1} \omega_{2} \left(I_{S} - I_{T}\right)}{I_{W} \left(I_{T} - I_{W}\right)}, \quad \dot{\omega}_{1} : \frac{\omega_{2} \left(I_{S} \omega_{3} - I_{T} \omega_{3} - I_{W} \Omega\right)}{I_{T}}, \quad \dot{\omega}_{2} : \frac{I_{W} \Omega \omega_{1}}{I_{S}}, \quad \dot{\omega}_{3} : \frac{D \Omega - \omega_{1} \omega_{2} \left(I_{S} - I_{T}\right)}{I_{T} - I_{W}}\right \}$$



In [49]:



In [51]:









    Out[51]:





$$\omega_{2} \left(- \Omega \frac{I_W}{I_T} + \frac{I_S}{I_T} \omega_{3} - \omega_{3}\right)$$

Finally, some numerical integration:



In [27]:

    
import numpy, scipy
import scipy.integrate
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt

Integrator state is: $z = \begin{bmatrix} \omega_1, \omega_2,\omega_3, \Omega \end{bmatrix}^T$. Note that there are common factors of $I_S/I_T$, $I_W/I_T$, and $D/I_T$. Grouping:



In [66]:

    
IsIt,IwIt, IwIs,DIt = symbols('\\frac{I_S}{I_T},\\frac{I_W}{I_T},\\frac{I_W}{I_S},\\frac{D}{I_T}')



In [106]:

    
dz = Matrix([simplify(expand(sol1[w1d]).subs(Is/It,IsIt).subs(Iw/It,IwIt)),
             simplify(expand(sol1[w2d]).subs(Iw/Is,IwIs)),
             (DIt*O - (IsIt - 1)*w1*w2)/(1 - IwIt),
             ((IsIt - 1)*w1*w2 - DIt/IwIt*O)/(1 - IwIt)])
dz









    Out[106]:





$$\left[\begin{matrix}\omega_{2} \left(- \Omega \frac{I_W}{I_T} + \frac{I_S}{I_T} \omega_{3} - \omega_{3}\right)\\\Omega \frac{I_W}{I_S} \omega_{1}\\\frac{\Omega \frac{D}{I_T} - \omega_{1} \omega_{2} \left(\frac{I_S}{I_T} - 1\right)}{- \frac{I_W}{I_T} + 1}\\\frac{- \frac{\Omega \frac{D}{I_T}}{\frac{I_W}{I_T}} + \omega_{1} \omega_{2} \left(\frac{I_S}{I_T} - 1\right)}{- \frac{I_W}{I_T} + 1}\end{matrix}\right]$$



In [65]:

    
simplify(sol1[w3d] - (D/It*O - (Is/It - 1)*w1*w2)/(1 - Iw/It))









    Out[65]:





$$0$$



In [107]:

    
simplify(((Is/It - 1)*w1*w2 - D/It/(Iw/It)*O)/(1 - Iw/It) - sol1[Od])









    Out[107]:





$$0$$

Pick some numerical values and initial conditions:



In [135]:

    
IsItn = 1.5
#IsItn = 1/1.5
IwItn = 0.06
IwIsn = IwItn/IsItn
DItn = 0.5

z0n = [0.2, 2, 0, 0]; #rad/s



In [136]:

    
dz1 = dz.subs(([IsIt,IsItn],[IwIt,IwItn],[IwIs,IwIsn],[DIt,DItn]))
dz1









    Out[136]:





$$\left[\begin{matrix}\omega_{2} \left(- 0.06 \Omega + 0.5 \omega_{3}\right)\\0.04 \Omega \omega_{1}\\0.531914893617021 \Omega - 0.531914893617021 \omega_{1} \omega_{2}\\- 8.86524822695035 \Omega + 0.531914893617021 \omega_{1} \omega_{2}\end{matrix}\right]$$



In [137]:

    
f1 = lambdify((w1,w2,w3,O),dz1)
f2 = lambda z,t: f1(z[0],z[1],z[2],z[3]).flatten()



In [138]:

    
t = numpy.linspace(0,400,1000)
y = scipy.integrate.odeint(f2,z0n,t)



In [139]:

    
plt.plot(t, y[:, 1], 'b', label='$\\omega_2(t)$')
plt.legend(loc='best')
plt.xlabel('Time (s)')
plt.ylabel('radians')
plt.grid()



In [140]:

    
plt.plot(t, y[:, 0], label='$\\omega_1(t)$')
plt.plot(t, y[:, 2], label='$\\omega_3(t)$')
plt.legend(loc='best')
plt.xlabel('Time (s)')
plt.ylabel('radians')
plt.grid()



In [141]:

    
plt.plot(t, y[:, 3], label='$\\Omega(t)$')
plt.legend(loc='best')
plt.xlabel('Time (s)')
plt.ylabel('radians')
plt.grid()



In [ ]: