# 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 = \
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.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.grid()







In [141]:

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







In [ ]: