<jupyter_start><jupyter_text>Universidade Federal do Rio Grande do Sul (UFRGS)   Programa de Pós-Graduação em Engenharia Civil (PPGEC)    PEC00025: Introduction to Vibration Theory Class 04 - Free vibration in time domain[1.   Free vibration of undamped _sdof_ systems](section_1)  [2.   Energy dissipation](section_2)  [3.   Free vibration of damped _sdof_ systems](section_3)  [4.   Damping ratio estimation](section_4)  [5.   Assignment](section_5)  ---_Prof. <NAME>, Dr.techn._ [(ORCID)](https://orcid.org/0000-0001-5640-1020)  _Porto Alegre, RS, Brazil_<jupyter_code># Importing Python modules required for this notebook
# (this cell must be executed with "shift+enter" before any other Python cell)

import numpy as np
import matplotlib.pyplot as plt
from   MRPy import MRPy<jupyter_output><empty_output><jupyter_text>1. Free vibration of undamped _sdof_ systems  The figure below shows an idealization of a single degree of freedom (sdof) mechanical systemwith its basic elements:1. The total system mass, $m$,2. A restitutive force element, ideally a linear spring with stiffness coefficient $k$,3. A dissipative force element, ideally a viscous damper with damping coefficient $c$,4. The applied external force, $F(t)$, and5. The system response, the displacement $u(t)$.Although the idealization presented above may resemble some feasible device, the basic elements may also be seen as conceptual and represent a much broader concept. For instance, the restitutive force element could be the gravity acting upon a pendulum.Or the system response could also be a rotation angle instead of a displacement.The dissipative force will be introduced in the next section and are initiallydisregarded. If we also disregard the external force, $F(t) = 0$, any system response will be called _free vibration_ and that is our first case for analysis.By applying the equilibrium equation and accounting for the inertial force as statedby Newton's law gives:$$ m \ddot{u} + k u = 0 $$Applying Laplace transform over this equation gives:$$ \mathscr{L}\left\{ m \ddot{u} + k u \right\} =     m \left[ -s u(0) - \dot{u}(0) + s^2 \bar{u}(s) \right] + k \bar{u}(s) = 0 $$Recognizing that $u(0) = u_0$ is the system initial position and $\dot{u}(0) = v_0$is the system initial velocity and solving for $\bar{u}(s)$ gives:$$ \bar{u}(s) =  \left(\frac{s}{s^2 + \omega_{\rm n}^2}\right) u_0 +    \left(\frac{\omega_{\rm n}}{s^2 + \omega_{\rm n}^2}\right) \frac{v_0}{\omega_{\rm n}}$$where we have defined $\omega_{\rm n}^2 = k/m$.The inverse transform is promptly recognized from the sine and cosine transformspreviously presented:$$ u(t) = u_0 \cos \omega_{\rm n}t + \frac{v_0}{\omega_{\rm n}} \sin \omega_{\rm n}t $$It can be seen that the free vibration occurs at the _natural vibration frequency_ $\omega_{\rm n}$, and the movement amplitude depends on the initial displacement and initial velocity. Let us plot this function with some Python lines:<jupyter_code>fn = 1.0                 # natural vibration frequency (1 cycle per second = Hz)
wn = 2*np.pi*fn          # conversion to rad/s
u0 = 5.                  # initial displacement of 1 meter
v0 = 30.                 # oscillation starts from rest (try changing this)

Td = 8.                  # time series duration (s)
fs = 32.                 # sampling rate (Hz), fs = 1/dt (dt is the time step)

t  = np.linspace(0, Td, int(fs*Td))

u  = u0*np.cos(wn*t) + (v0/wn)*np.sin(wn*t)

plt.figure(1, figsize=(8, 4), clear=True)
plt.plot(t, u)

plt.xlim(  0, Td);  plt.xlabel('time (s)') 
plt.ylim(-10, 10);  plt.ylabel('u(t)') 

plt.grid(True)<jupyter_output><empty_output><jupyter_text>This is surely not an interesting result. Try changing the values for $u_0$ and $v_0$ tosee what happens. Indeed, the solution above can also be given in the form:$$ u(t) = u_{\rm max} \cos ( \omega_{\rm n}t - \theta ) $$with:\begin{align*} u_{\rm max} &=  \sqrt{u_0^2 + \left(\frac{v_0}{\omega_{\rm n}}\right)^2} \\ \\ \theta      &=  \arctan \left(\frac{v_0}{u_0 \omega_{\rm n}} \right)\end{align*}The most important observation so far is the natural vibration frequency, that holds asort of universal law for mechanical vibrations:$$ \omega_{\rm n} = \sqrt{\frac{k}{m}} $$ Hence, the larger the system stiffness the higher the frequency, while thelarger the system mass the lower the frequency, and vice-versa.This special frequency can be also obtained by energy conservation, if one recognizes that there are two kinds of energy present in the system:1. Potential elastic energy: $$V = k u^2(t)/2$$2. Kinectic energy: $$T = m \dot{u}^2(t)/2$$Assuming that if the maximum displacement of the system is $u_{\rm max}$, then the maximum velocity will be $\omega_{\rm n} u_{\rm max}$. When the displacement attains its maximum the potential energy is maximumand the kinectic energy is minimun, and vice-versa. By energy conservation:$$ \frac{k \; u_{\rm max}^2}{2} = \frac{m \; \omega_{\rm n}^2 \; u_{\rm max}^2}{2} $$what leads to the equation for $\omega_{\rm n}$.The natural frequency is also expressed in cycles per second, or _Hertz_, related to the frequency in radians per second with:$$ f_{\rm n} =\frac{\omega_{\rm n}}{2 \pi} $$Alternatively, the _natural period of vibration_:$$ T_{\rm n} = \frac{1}{f_{\rm n}} = \frac{2 \pi } {\omega_{\rm n}} $$can be used. These relations are illustrated in the figure below.      2. Energy dissipation  Besides the potential and kinectic energy, there is always some energy dissipation mechanism taking place during the system vibration (otherwise a perpetual motion would be possible). The most common causes for this dissipation are:1. Material internal hysteresis,  2. Shear deformation of embedding fluid,  3. Friction at supports or connections,   among others. Each one of these phenomena must be properly modelled to enter theequilibrium equation, but the most common approach is to use theso-called _newtonian_ or _viscous_ damping model for thedissipative force, $c \dot{u}$, where $c$ is the viscous damping constant.The great advantage of this model is that the equilibrium equation preservesits linearity and solution easiness, but one must keep in mind that it is_only a model_ and the definition of constant $c$ may be quite challenging.Most commonly the damping constant is presented as a _ratio to critical damping_, $\zeta$, defined as:$$ \zeta = \frac{c}{2 m \omega_{\rm n}} $$which is non-dimensional and for this reason more easily quantified.It can be shown that $\zeta = 1$ is the smallest damping prenventing the system from having a single oscillation (one can have an idea of this condition by releasing a pendulum under water).The Chapter 9 of the Brazilian Code _NBR6123-Forças devidas o vento nasedificações_ brings a table (_Tabela 2_) of suggested values for $\zeta$ according to the type of structural system undergoing wind induced vibration:|  Structure type                                           | $\zeta$  | | :---                                                      |   :---:  | | Reinforced concrete buildings without stiffning walls     |   0.020  | | Reinforced concrete buildings with stiffning walls        |   0.015  | | Reinforced concrete towers and chimneys, tappered         |   0.015  | | Reinforced concrete towers and chimneys, constant section |   0.010  | | Moment resisting steel buildings with welded connections  |   0.010  | | Steel towers and chimneys, constant section               |   0.008  || Timber structures                                         |   0.030  |It can be observed that for Civil Engineering the range of values arenot too wide, and a typical value $\zeta = 0.01 = 1\%$ is usually assumedwhen no further information is available.In the next section it will be shown that despite of all uncertainty aboutenergy dissipation, it is of utmost importance for the safety and goodperformance of mechanical systems.  3. Free vibration of damped _sdof_ systems  Assuming now the viscous model for energy dissipation gives the free vibration equilibrium:$$ m \ddot{u} + c \dot{u} + k u = 0$$Applying Laplace transform over this equation gives:$$ \mathscr{L}\left\{ m \ddot{u} + c \dot{u} + k u \right\} =     m \left[ -s u(0) - \dot{u}(0) + s^2 \bar{u}(s) \right] +     c \left[ -u(0) + s \bar{u}(s) \right] +    k \bar{u}(s) = 0 $$Solving for $\bar{u}(s)$ results in:$$ \bar{u}(s) =     \left[\frac{(s + \eta)}{(s + \eta)^2 + \omega_{\rm D}^2}\right] u_0 +    \left[\frac{\omega_{\rm D}}{(s + \eta)^2 + \omega_{\rm D}^2}\right]    \left(\frac{v_0 + \eta u_0}{\omega_{\rm D}} \right) $$where we have additionally defined:\begin{align*}\omega_{\rm D}^2 &= \omega_{\rm n}^2 - \eta^2 = \omega_{\rm n}^2 \left( 1 - \zeta^2 \right ) \\ \\\eta &= \frac{c}{2m} = \zeta \omega_{\rm n}\end{align*}The inverse transform can be obtained with the well known sine and cosine transforms and the translation theorem in frequency domain:$$ u(t) = e^{- \zeta \omega_{\rm n} t}   \left[ u_0 \cos \omega_{\rm D}t + \left( \frac{v_0 + \zeta \omega_{\rm n} u_0}   {\omega_{\rm D}} \right) \sin \omega_{\rm D}t  \right] $$It can be seen that the free vibration now occurs at the _natural damped vibration frequency_ $\omega_{\rm D}$, and the movement amplitude depends on the initial displacement and initial velocity. Although the natural frequency has changed, by recalling that the damping ration $\zeta$ usually is in the order of 1% implies, for instance, that:$$ \omega_{\rm D} = \omega_{\rm n} \sqrt{1 - 0.01^2}                     \approx 0.99995 \omega_{\rm n} \approx \omega_{\rm n}$$what means that, beyond the whole uncertainty about the actual damping, for usual damping values the natural frequency change may be disregarded.Let us take a look at the free vibration damped responseusing the same time discretization from the previous example.<jupyter_code>fn = 1.0                 # natural vibration frequency (1 cycle per second = Hz)
wn = 2*np.pi*fn          # conversion to rad/s
zt = 0.05                # damping ratio of 3% of the critical (try changing this)
u0 = 5.                  # initial displacement of 1 meter
v0 = 0.                  # oscillation starts from rest 

Td = 8.                  # time series duration (s)
fs = 32.                 # sampling rate (Hz), fs = 1/dt (dt is the time step)

t  = np.linspace(0, Td, int(fs*Td))

wD = wn*np.sqrt(1 - zt**2)
e  = np.exp(-zt*wn*t)
uD = e*(u0*np.cos(wD*t) + ((v0 + zt*wn*u0)/wD)*np.sin(wD*t))

plt.figure(2, figsize=(8, 4), clear=True)
plt.plot(t, uD)
plt.plot(t, u0*e, 'r:', t, -u0*e, 'r:')

plt.xlim(  0, Td);  plt.xlabel('time (s)') 
plt.ylim(-10, 10);  plt.ylabel('u(t)') 

plt.grid(True)<jupyter_output><empty_output><jupyter_text>It can be observed that the exponential envelope causes the undamped response to haveits amplitude decreased along time. Try changing the damping ratio in the script aboveto see how this affect the response (but it must always be within the range $0 < \zeta < 1$).  4. Damping ratio estimation  As we did for the undamped response, the damped response in free vibration can be put inthe form:$$ u(t) = u_{\rm max} e^{- \zeta \omega_{\rm n} t} \cos ( \omega_{\rm D}t - \theta ) $$with:\begin{align*} u_{\rm max} &=  \sqrt{u_0^2 +                  \left( \frac{v_0 + \zeta \omega_{\rm n} u_0} {\omega_{\rm D}} \right) ^2} \\ \\ \theta      &=  \arctan \left( \frac{v_0 + \zeta \omega_{\rm n} u_0}{\omega_{\rm D} u_0 }\right)\end{align*}Now we calculate the amplitude ratio at two time instants where the cosinefunction reaches its maximum, which are $N$ vibration cycles apart, what gives:$$ \frac{u(t)}{u(t + N T_{\rm n})} =    \frac{\exp \left[ - \zeta \omega_{\rm n}  t                \right]}        {\exp \left[ - \zeta \omega_{\rm n} (t + N T_{\rm n}) \right]} =          \exp \left( - \zeta \omega_{\rm n} N T_{\rm n}       \right)  =         \exp \left( - 2 \pi \zeta N  \right) $$   Making $u(t) = u_i$, $u(t + N T_{\rm n}) = u_{i+N}$, taking natural logarithms and isolating $\zeta$ finally gives:$$ \zeta = \frac {\log (u_i \, / \, u_{i+N})}{2 \pi N} $$This formula is exemplified in the figure below with the output of a real accelerometer.It is important to emphasize that the estimator presented above is valid no matterwhich kinemactic parameter is used for the amplitude ratio (displacement, velocity oracceleration), for they are related by a constant that will cancel upon division.  Applying the formula gives:$$ \zeta = \frac{\ln (0.52 \, / \, 0.31)}{2 \pi \cdot 8} \approx 1\% $$A usefull rule for estimating $\zeta$ without much effort is to count the numberof cycles necessary for the amplitude to be reduced to half. In this case:$$ \zeta = \frac{\ln (2 \, / \, 1)}{2 \pi \,N} \approx \frac{1}{9N} $$The ```MRPy``` module provides a more sophisticated tool for estimating all parameters of a free damped response function available as a time series. This is exemplifyed below using the same the $u(t)$ function from the previous example.<jupyter_code>fn = 1.0                 # natural vibration frequency (1 cycle per second = Hz)
wn = 2*np.pi*fn          # conversion to rad/s
zt = 0.01                # damping ratio of 1% of the critical (try changing this)
u0 = 1.                  # initial displacement of 1 meter
v0 = 0.                  # oscillation starts from rest 

Td = 8.                  # time series duration (s)
fs = 32.                 # sampling rate (Hz), fs = 1/dt (dt is the time step)

t  = np.linspace(0, Td, int(fs*Td))

wD = wn*np.sqrt(1 - zt**2)
e  = np.exp(-zt*wn*t)
uD = e*(u0*np.cos(wD*t) + ((v0 + zt*wn*u0)/wD)*np.sin(wD*t) + 0.1*u0*np.cos(0.7*wD*t)) 

uD = MRPy(uD, fs=fs)                # converting numpy to MRPy object
        
uD = uD.extract((1/8, 7/8), by='fraction')

ufit, p = uD.fit_decay()            # calls fit_decay() method
t       = uD.t_axis()               # prepare time axis for plotting

plt.figure(3, figsize=(12, 6), clear=True)
plt.plot(t, uD[0])
plt.plot(t, ufit[0], 'r:')
plt.legend(('original','fitted'))

plt.xlim( 0, 10);  plt.xlabel('time (s)') 
plt.ylim(-2,  2);  plt.ylabel('u(t)') 

plt.grid(True) 

# Print out fitted parameters

print('Amplitude:     {0:4.2f}m  '.format(p[0,0]))
print('Frequency:     {0:4.2f}Hz '.format(p[0,1]))
print('Damping ratio: {0:4.2f}%  '.format(p[0,2]*100))
print('Phase:         {0:4.2f}rad'.format(p[0,3]))<jupyter_output>Amplitude:     0.90m  
Frequency:     1.00Hz 
Damping ratio: 0.79%  
Phase:         -0.01rad<jupyter_text>The fitting result is usually not so accurate when the original function comes fromreal measurements. In some cases the fitting algorithm (which uses ```scipy``` module)may even fail in presenting any usefull result.<jupyter_code>data = MRPy.from_file('resources/data/iNVH001', form='invh')

data.printAttrib()
data.plot_time(fig=1);

free = data.extract((2.5, 5.5), by='time')

free.printAttrib()
free.plot_time(fig=2);

afit, par = free.fit_decay()
t         = free.t_axis()               # prepare time axis for plotting

plt.figure(4, figsize=(12, 6), clear=True)
plt.plot(t, free[2])
plt.plot(t, afit[2], 'r:')
plt.legend(('original','fitted'))

plt.xlim( 0, 3);  plt.xlabel('time (s)') 
plt.ylim(-3, 3);  plt.ylabel('u(t)') 

plt.grid(True) 

print('Amplitude:     {0:4.2f}m/s2'.format(par[2,0]))
print('Frequency:     {0:4.2f}Hz  '.format(par[2,1]))
print('Damping ratio: {0:4.2f}%   '.format(par[2,2]*100))
print('Phase:         {0:4.2f}rad '.format(par[2,3]))<jupyter_output>Amplitude:     2.10m/s2
Frequency:     8.34Hz  
Damping ratio: 1.25%   
Phase:         -0.89rad
{'fraction_non_alphanumeric': 0.13179145939660183, 'fraction_numerical': 0.02448478583887848, 'mean_word_length': 3.348314606741573, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 20, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '569111', 'license': 'MIT', 'n_tokens_mistral': 5151, 'n_tokens_neox': 4812, 'n_words': 1905, 'path': 'Class_04_FreeTimeDomain.ipynb', 'repo_name': 'mmaiarocha/PEC00025', 'size': '15489'}