CIE_4140_Lecture_3_2_Python

Free vibration of an SDOF with viscous damping under harmonic excitation

import sympy as sp

u = sp.symbols('u',cls=sp.Function)
t, omega_n, f_0, m, Omega, zeta = sp.symbols('t, omega_n, f_0, m, Omega, zeta',real=True,positive=True)
u_0, v_0 = sp.symbols('u_0, v_0')

Equation_of_Motion= sp.Eq(sp.diff(u(t),t,2)+2*zeta*omega_n*sp.diff(u(t),t)+omega_n **2 * u(t),f_0*sp.cos(Omega*t)/m)
display(Equation_of_Motion)

\[\displaystyle \omega_{n}^{2} u{\left(t \right)} + 2 \omega_{n} \zeta \frac{d}{d t} u{\left(t \right)} + \frac{d^{2}}{d t^{2}} u{\left(t \right)} = \frac{f_{0} \cos{\left(\Omega t \right)}}{m}\]

u_sol_general = sp.dsolve(Equation_of_Motion, u(t)).rhs
display(u_sol_general)

\[\displaystyle C_{1} e^{\omega_{n} t \left(- \zeta + \sqrt{\zeta - 1} \sqrt{\zeta + 1}\right)} + C_{2} e^{- \omega_{n} t \left(\zeta + \sqrt{\zeta - 1} \sqrt{\zeta + 1}\right)} - \frac{\Omega^{2} f_{0} \cos{\left(\Omega t \right)}}{m \left(\Omega^{4} + 4 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} + \omega_{n}^{4}\right)} + \frac{2 \Omega f_{0} \omega_{n} \zeta \sin{\left(\Omega t \right)}}{m \left(\Omega^{4} + 4 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} + \omega_{n}^{4}\right)} + \frac{f_{0} \omega_{n}^{2} \cos{\left(\Omega t \right)}}{m \left(\Omega^{4} + 4 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} + \omega_{n}^{4}\right)}\]

The expresion above can be written in the following real-valued form:

A, B = sp.symbols('A, B',real=True)
u_sol_general = (sp.exp(-omega_n*t*zeta) * (A*sp.cos(omega_n*t*sp.sqrt(-zeta**2 + 1)) +
                                   B*sp.sin(omega_n*t*sp.sqrt(-zeta**2 + 1))) + 
         f_0*(2*sp.sin(Omega*t)*Omega*zeta*omega_n -
              sp.cos(Omega*t)*Omega**2 +
              sp.cos(Omega*t)*omega_n**2)
         /(omega_n**4 + (4*zeta**2 - 2)*Omega**2*omega_n**2 + Omega**4))
display(u_sol_general)

\[\displaystyle \frac{f_{0} \left(- \Omega^{2} \cos{\left(\Omega t \right)} + 2 \Omega \omega_{n} \zeta \sin{\left(\Omega t \right)} + \omega_{n}^{2} \cos{\left(\Omega t \right)}\right)}{\Omega^{4} + \Omega^{2} \omega_{n}^{2} \cdot \left(4 \zeta^{2} - 2\right) + \omega_{n}^{4}} + \left(A \cos{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)} + B \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)}\right) e^{- \omega_{n} t \zeta}\]

IC1 = sp.Eq(u_sol_general.subs(t,0),u_0)
display(IC1)

\[\displaystyle A + \frac{f_{0} \left(- \Omega^{2} + \omega_{n}^{2}\right)}{\Omega^{4} + \Omega^{2} \omega_{n}^{2} \cdot \left(4 \zeta^{2} - 2\right) + \omega_{n}^{4}} = u_{0}\]

IC2 = sp.Eq(u_sol_general.diff(t).subs(t,0),v_0)
display(IC2)

\[\displaystyle - A \omega_{n} \zeta + B \omega_{n} \sqrt{1 - \zeta^{2}} + \frac{2 \Omega^{2} f_{0} \omega_{n} \zeta}{\Omega^{4} + \Omega^{2} \omega_{n}^{2} \cdot \left(4 \zeta^{2} - 2\right) + \omega_{n}^{4}} = v_{0}\]

sol = sp.solve((IC1,IC2),(A,B))
sol[A] = sp.simplify(sol[A])
sol[B] = sp.simplify(sol[B])
display(sol[A])
display(sol[B])

\[\displaystyle \frac{\Omega^{4} u_{0} + \Omega^{2} f_{0} + 4 \Omega^{2} \omega_{n}^{2} u_{0} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} u_{0} - f_{0} \omega_{n}^{2} + \omega_{n}^{4} u_{0}}{\Omega^{4} + 4 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} + \omega_{n}^{4}}\]

\[\displaystyle \frac{\Omega^{4} \omega_{n} u_{0} \zeta + \Omega^{4} v_{0} - \Omega^{2} f_{0} \omega_{n} \zeta + 4 \Omega^{2} \omega_{n}^{3} u_{0} \zeta^{3} - 2 \Omega^{2} \omega_{n}^{3} u_{0} \zeta + 4 \Omega^{2} \omega_{n}^{2} v_{0} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} v_{0} - f_{0} \omega_{n}^{3} \zeta + \omega_{n}^{5} u_{0} \zeta + \omega_{n}^{4} v_{0}}{\omega_{n} \sqrt{1 - \zeta^{2}} \left(\Omega^{4} + 4 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} + \omega_{n}^{4}\right)}\]

u_sol = u_sol_general.subs(sol)
display(u_sol)

\[\displaystyle \frac{f_{0} \left(- \Omega^{2} \cos{\left(\Omega t \right)} + 2 \Omega \omega_{n} \zeta \sin{\left(\Omega t \right)} + \omega_{n}^{2} \cos{\left(\Omega t \right)}\right)}{\Omega^{4} + \Omega^{2} \omega_{n}^{2} \cdot \left(4 \zeta^{2} - 2\right) + \omega_{n}^{4}} + \left(\frac{\left(\Omega^{4} u_{0} + \Omega^{2} f_{0} + 4 \Omega^{2} \omega_{n}^{2} u_{0} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} u_{0} - f_{0} \omega_{n}^{2} + \omega_{n}^{4} u_{0}\right) \cos{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)}}{\Omega^{4} + 4 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} + \omega_{n}^{4}} + \frac{\left(\Omega^{4} \omega_{n} u_{0} \zeta + \Omega^{4} v_{0} - \Omega^{2} f_{0} \omega_{n} \zeta + 4 \Omega^{2} \omega_{n}^{3} u_{0} \zeta^{3} - 2 \Omega^{2} \omega_{n}^{3} u_{0} \zeta + 4 \Omega^{2} \omega_{n}^{2} v_{0} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} v_{0} - f_{0} \omega_{n}^{3} \zeta + \omega_{n}^{5} u_{0} \zeta + \omega_{n}^{4} v_{0}\right) \sin{\left(\omega_{n} t \sqrt{1 - \zeta^{2}} \right)}}{\omega_{n} \sqrt{1 - \zeta^{2}} \left(\Omega^{4} + 4 \Omega^{2} \omega_{n}^{2} \zeta^{2} - 2 \Omega^{2} \omega_{n}^{2} + \omega_{n}^{4}\right)}\right) e^{- \omega_{n} t \zeta}\]

Development of resoncance in time

solution_in_resonance = u_sol.subs([(Omega,omega_n*sp.sqrt(1-zeta**2))])
sp.plot(solution_in_resonance.subs([(u_0,0.1),(v_0,0.1),(f_0,1),(omega_n,1),(zeta,0.05)]), (t , 0 , 152));

%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

u_resonance_func = sp.lambdify((u_0,v_0,omega_n,zeta,f_0,t),solution_in_resonance)

fig, ax = plt.subplots()
tdata = np.linspace(0,100,500)
line, = ax.plot([], [])
ax.set_xlim(0, 100)
ax.set_ylim(-30, 30)

def update(frame):
    ydata = u_resonance_func(u_0=0.1,v_0=1,omega_n=1,f_0=1,zeta=frame,t=tdata)
    ax.set_title("Displacement versus time for $\zeta$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(0.01,0.8,100),interval = 100)
plt.show()

The steady-state solution

u_steady = u_sol.args[1]
display(u_steady)

\[\displaystyle \frac{f_{0} \left(- \Omega^{2} \cos{\left(\Omega t \right)} + 2 \Omega \omega_{n} \zeta \sin{\left(\Omega t \right)} + \omega_{n}^{2} \cos{\left(\Omega t \right)}\right)}{\Omega^{4} + \Omega^{2} \omega_{n}^{2} \cdot \left(4 \zeta^{2} - 2\right) + \omega_{n}^{4}}\]

p0 = sp.plotting.plot(u_steady.subs([(u_0,0.1),(v_0,0.1),(f_0,1),(omega_n,1),(zeta,0.05),(Omega,2.3)]), (t , 0 , 152),label='Steady-state' ,legend=True,show=False,adaptive=False,nb_of_points=3000)
p1 = sp.plotting.plot(   u_sol.subs([(u_0,0.1),(v_0,0.1),(f_0,1),(omega_n,1),(zeta,0.05),(Omega,2.3)]), (t , 0 , 152),label='Full solution',legend=True,show=False,adaptive=False,nb_of_points=3000)
p0.append(p1[0])
p0.show()

u_func = sp.lambdify((u_0,v_0,omega_n,zeta,f_0,t,Omega),u_sol)

fig, ax = plt.subplots()
tdata = np.linspace(0,100,500)
line, = ax.plot([], [])
ax.set_xlim(0, 100)
ax.set_ylim(-0.5, 0.5)

def update(frame):
    ydata = u_func(u_0=0.1,v_0=0.1,f_0=1,omega_n=1,Omega=2.3,zeta=frame,t=tdata)
    ax.set_title("Displacement versus time for $\zeta$ = "+str(np.round(frame,2)))
    line.set_data(tdata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(0.01,.8,100),interval = 100)
plt.show()

The magnification factor

Omega_dl = sp.symbols('Omega_dl')
magnification_factor = 1 / sp.sqrt((1-Omega_dl**2)**2+4*zeta**2*Omega_dl**2)

magnification_factor_func = sp.lambdify((Omega_dl,zeta), magnification_factor)

fig, ax = plt.subplots()
frequencydata = np.linspace(0,3,100)
line, = ax.plot([], [])
ax.set_xlim(0, 3)
ax.set_ylim(0,25)

def update(frame):
    ydata = magnification_factor_func(Omega_dl = frequencydata , zeta = frame)
    ax.set_title("Magnification factor versus frequency for $\zeta$ = "+str(np.round(frame,2)))
    line.set_data(frequencydata, ydata)

ani = FuncAnimation(fig, update, frames=np.linspace(0.02,1.3,100),interval = 100)
plt.show()