As an example we will study a model of washing machine vibrating. A washing machine typically has a structure that sits on the ground, usually on padded feet that provide some damping and then there is a drum that spins relative to the machine's main housing structure. Inside the drum, clothes can end up clumped in ways that cause the mass center of the contents of the drum to lie somewhere other than the spin axis. See this video for an example of what can happen:
from IPython.display import YouTubeVideo
YouTubeVideo('JqdAKPzdHDo')
For nonlinear analysis:
For linear analysis:
import sympy as sm
Here is a hack that will let $\frac{d y}{dt}$ appear as $\dot{y}$.
from sympy.printing.latex import LatexPrinter
class TimeDerivativePrinter(LatexPrinter):
def _print_Derivative(self, expr):
dim = len(expr.variables)
t = sm.Symbol('t')
if dim == 1 and expr.variables == (t,):
return r"\dot{{{}}}".format(self._print(expr.args[0]))
elif dim == 2 and expr.variables == (t, t):
return r"\ddot{{{}}}".format(self._print(expr.args[0]))
else:
return self._print(expr)
def latex(expr, **settings):
return TimeDerivativePrinter(settings).doprint(expr)
sm.init_printing(latex_printer=latex)
I want to make the model (at least the first one) to be as simple as I can such that we can try to show the vibrations. Here are some assumptions:
With those in mind, here is a basic free body diagram of the system:
mw, mc, md, Icd, k, c, e, w, g, y0 = sm.symbols('m_w, m_c, m_d, I_cd, k, c, e, omega, g, y_0')
All generalized coordinates should be functions of time.
t = sm.symbols('t')
y = sm.Function('y')(t)
y
y.diff(t)
y.diff(t, 2)
T = (mw * y.diff()**2 / 2 + # vertical motion of washing machine
md * y.diff()**2 / 2 + # vertical motion of the drum
Icd * w**2 / 2 + # rotation of the drum
mc * ((y.diff() + w * e * sm.cos(w * t))**2 + (w * e * sm.sin(w * t))**2) / 2) # motion of the clothes mass center
T
U = (mw * g * y + # due to raising height of washing machine
md * g * y + # due to raising height of drum
mc * g * (y + e * sm.sin(w * t)) + # due to raising height of clothes
k * (y - y0)**2 / 2) # due to stretch away from free length of spring
U
R = c * y.diff()**2 / 2 # due to dashpot
R
L = T - U
L
There is only one generalized coordinate, $y$, so there will be one equation of motion for this single degree of freedom system.
$$f(y, \dot{y}, \ddot{y}) = \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{y}}\right) - \frac{\partial L}{\partial y} + \frac{\partial R}{\partial\dot{y}} = 0$$f = L.diff(y.diff()).diff(t) - L.diff(y) + R.diff(y.diff())
f
.expand()
and .collect()
can be used to show the equation in terms of each $y$ term.
f = f.expand().collect((y, y.diff(), y.diff(t, 2)))
f
Notice that the rotational centroidal inertia of the drum plays no role in the dynamics because it can't affect the vertical motion of the machine. We could have ignored that from the beginning if we made that realization and further simplified the model description. Note also that the left hand side of this equation is linear in $y,\dot{y},\ddot{y}$. Lastly, note the constant term $$-m_T g + k y_0$$. This term has to do with the free length and the total weight.
$$m_T \ddot{y} + c \dot{y} + k y = em_c\omega^2\sin(\omega t) - m_T g + k y_0$$where
$$m_T = m_c + m_d + m_w$$The above equations are valid, but with linear systems, we are typically interested in vibration about the equilibrium. The value of $y$ for equilibrium should be the value of the free length of the spring minus the distance compressed under the static load of the machine, i.e.:
$$y_{eq} = y_0 - \frac{m_T g}{k}$$If you set the velocities and accelerations to zero in the equation of motion, this leaves the static force balance and solve it for the generalized coordinate will give the equilibrium value.
yeq = sm.solve(f.subs({y.diff(t, 2): 0, y.diff(t): 0, w: 0}), y)[0] # note that solve returns a list, so [0] is needed
yeq
Now that we know the value of $y$ the system is at a equilibrium, we can write the equation of motion in terms it's motion with respect to the motion $\Delta y$ about the equilibrium, i.e. $y = y_{eq} + \Delta y$.
dely = sm.Function('\Delta y')(t)
dely
f_eq = f.subs({y.diff(t, 2): dely.diff(t, 2)}).subs({y.diff(t): dely.diff(t)}).subs({y: yeq + dely}).simplify()
f_eq
Notice how the gravity related force terms and the need for the spring's free length were eliminated.
The equation above also is linear wrt $\Delta y$, $\Delta\dot{y}$, and $\Delta\ddot{y}$. The $m$, $c$, and $k$ coefficients can be extracted:
m_coeff = f_eq.coeff(dely.diff(t, 2))
m_coeff
c_coeff = f_eq.coeff(dely.diff(t))
c_coeff
k_coeff = f_eq.coeff(dely)
k_coeff
The remaining term is the right hand side forcing term:
forcing = -(f_eq - m_coeff * dely.diff(t, 2) - c_coeff * dely.diff(t) - k_coeff * dely)
forcing
Now we have a linear single degree of freedom system with a sinusoidal forcing term that takes this form:
$$ m \Delta\ddot{y}(t) + c \Delta\dot{y}(t) + k \Delta y(t) = f_0 \sin(\omega t) $$sm.Eq(m_coeff * dely.diff(t, 2) + c_coeff * dely.diff(t) + k_coeff * dely, forcing)
Notice that the equation of motion of vibration about the equilibrium does not involve the free length or the weight of the machine. If you originally define your origin with respect to the static equilibrium you can write $U$ without the gravity potential energy terms related to the total weight because the energy stored in the compressed spring always equates to the energy gained from raising or lowering the mass. The free body diagram would look more like:
So the equation of motion can be formed more easily like:
U = (mc * g * e * sm.sin(w * t) + # potential from height raise of mass of the clothes only
k * y**2 / 2) # potential due to stretch of spring away from equilibrium
U
L = T - U
f = L.diff(y.diff()).diff(t) - L.diff(y) + R.diff(y.diff())
sm.Eq(f.expand().collect(y.diff(t, 2)), 0)
Same equation of motion as above, but with fewer steps!