vibration_analysis.tex

\documentclass[11pt,a4paper]{article}

% Document Setup
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amssymb}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{float}
\usepackage{caption}
\usepackage{hyperref}

% PythonTeX Setup
\usepackage[makestderr]{pythontex}

\title{Vibration Analysis: SDOF Systems and Modal Analysis}
\author{Mechanical Engineering Laboratory}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This report presents computational analysis of mechanical vibrations including free and forced response of single degree of freedom (SDOF) systems, damping effects, frequency response, and modal analysis. Python-based computations provide quantitative analysis with dynamic visualization.
\end{abstract}

\tableofcontents
\newpage

\section{Introduction to Mechanical Vibrations}

Vibration analysis is essential for:
\begin{itemize}
    \item Machine design and reliability
    \item Structural dynamics and earthquake engineering
    \item Noise and vibration control
    \item Condition monitoring and diagnostics
\end{itemize}

% Initialize Python environment
\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.signal import lti, step, bode

plt.rcParams['figure.figsize'] = (8, 5)
plt.rcParams['font.size'] = 10
plt.rcParams['text.usetex'] = True

def save_fig(filename):
    plt.savefig(filename, dpi=150, bbox_inches='tight')
    plt.close()
\end{pycode}

\section{SDOF System Fundamentals}

\subsection{Equation of Motion}

For a mass-spring-damper system:
\begin{equation}
m\ddot{x} + c\dot{x} + kx = F(t)
\end{equation}

Natural frequency and damping ratio:
\begin{equation}
\omega_n = \sqrt{\frac{k}{m}}, \quad \zeta = \frac{c}{2\sqrt{km}} = \frac{c}{2m\omega_n}
\end{equation}

\begin{pycode}
# SDOF system parameters
m = 1.0    # kg
k = 100    # N/m
c = 4.0    # Ns/m

omega_n = np.sqrt(k/m)
zeta = c / (2*np.sqrt(k*m))
omega_d = omega_n * np.sqrt(1 - zeta**2)
f_n = omega_n / (2*np.pi)

# Free vibration response
def free_vibration(y, t, m, c, k):
    x, v = y
    return [v, (-c*v - k*x)/m]

t = np.linspace(0, 5, 1000)
y0 = [0.1, 0]  # Initial displacement, zero velocity
sol = odeint(free_vibration, y0, t, args=(m, c, k))

# Analytical solution for underdamped
x_analytical = y0[0] * np.exp(-zeta*omega_n*t) * \
               (np.cos(omega_d*t) + zeta*omega_n/omega_d * np.sin(omega_d*t))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Time response
ax1.plot(t, sol[:, 0]*1000, 'b-', linewidth=2, label='Numerical')
ax1.plot(t, x_analytical*1000, 'r--', linewidth=1, label='Analytical')
ax1.plot(t, y0[0]*1000*np.exp(-zeta*omega_n*t), 'k--', linewidth=1, alpha=0.5, label='Envelope')
ax1.plot(t, -y0[0]*1000*np.exp(-zeta*omega_n*t), 'k--', linewidth=1, alpha=0.5)
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Displacement (mm)')
ax1.set_title('Free Vibration Response')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Phase portrait
ax2.plot(sol[:, 0]*1000, sol[:, 1]*1000, 'b-', linewidth=1)
ax2.plot(sol[0, 0]*1000, sol[0, 1]*1000, 'go', markersize=10, label='Start')
ax2.plot(sol[-1, 0]*1000, sol[-1, 1]*1000, 'ro', markersize=10, label='End')
ax2.set_xlabel('Displacement (mm)')
ax2.set_ylabel('Velocity (mm/s)')
ax2.set_title('Phase Portrait')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
save_fig('free_vibration.pdf')
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{free_vibration.pdf}
\caption{Free vibration: time response with envelope decay and phase portrait.}
\end{figure}

\begin{table}[H]
\centering
\caption{SDOF System Parameters}
\begin{tabular}{lcc}
\toprule
Parameter & Value & Units \\
\midrule
Natural frequency $\omega_n$ & \py{f"{omega_n:.2f}"} & rad/s \\
Natural frequency $f_n$ & \py{f"{f_n:.2f}"} & Hz \\
Damping ratio $\zeta$ & \py{f"{zeta:.3f}"} & -- \\
Damped frequency $\omega_d$ & \py{f"{omega_d:.2f}"} & rad/s \\
\bottomrule
\end{tabular}
\end{table}

\section{Effect of Damping}

\begin{pycode}
# Effect of damping ratio on response
zeta_values = [0.05, 0.1, 0.2, 0.5, 1.0, 2.0]
colors = plt.cm.viridis(np.linspace(0, 1, len(zeta_values)))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

for zeta_i, color in zip(zeta_values, colors):
    c_i = 2*zeta_i*np.sqrt(k*m)
    sol_i = odeint(free_vibration, y0, t, args=(m, c_i, k))
    ax1.plot(t, sol_i[:, 0]*1000, color=color, linewidth=1.5, label=f'$\\zeta$ = {zeta_i}')

ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Displacement (mm)')
ax1.set_title('Effect of Damping on Free Response')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)

# Logarithmic decrement
# For underdamped systems
delta = 2*np.pi*zeta / np.sqrt(1 - zeta**2) if zeta < 1 else 0

# Peak amplitudes for calculating log decrement
peaks_idx = []
for i in range(1, len(sol)-1):
    if sol[i, 0] > sol[i-1, 0] and sol[i, 0] > sol[i+1, 0]:
        peaks_idx.append(i)

if len(peaks_idx) >= 2:
    x1 = sol[peaks_idx[0], 0]
    x2 = sol[peaks_idx[1], 0]
    delta_measured = np.log(x1/x2)
    zeta_measured = delta_measured / np.sqrt(4*np.pi**2 + delta_measured**2)
else:
    delta_measured = 0
    zeta_measured = 0

# Show peak decay
ax2.semilogy(t[peaks_idx], np.abs(sol[peaks_idx, 0])*1000, 'ro-', linewidth=2, markersize=8)
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Peak Amplitude (mm)')
ax2.set_title(f'Logarithmic Decrement: $\\delta$ = {delta_measured:.3f}')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
save_fig('damping_effect.pdf')
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{damping_effect.pdf}
\caption{Effect of damping ratio on free vibration and logarithmic decrement measurement.}
\end{figure}

\section{Forced Vibration Response}

\subsection{Harmonic Excitation}

For $F(t) = F_0\sin(\omega t)$, the steady-state response is:
\begin{equation}
X = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}, \quad r = \frac{\omega}{\omega_n}
\end{equation}

\begin{pycode}
# Frequency response function
r = np.linspace(0.01, 3, 500)  # Frequency ratio
zeta_values = [0.05, 0.1, 0.2, 0.5, 0.7, 1.0]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Magnitude
for zeta_i in zeta_values:
    X_ratio = 1 / np.sqrt((1 - r**2)**2 + (2*zeta_i*r)**2)
    ax1.plot(r, X_ratio, linewidth=1.5, label=f'$\\zeta$ = {zeta_i}')

ax1.axvline(1, color='k', linestyle='--', alpha=0.5)
ax1.set_xlabel('Frequency Ratio $r = \\omega/\\omega_n$')
ax1.set_ylabel('Amplitude Ratio $X/(F_0/k)$')
ax1.set_title('Frequency Response - Magnitude')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 3)
ax1.set_ylim(0, 10)

# Phase
for zeta_i in zeta_values:
    phi = np.arctan2(2*zeta_i*r, 1 - r**2)
    ax2.plot(r, np.degrees(phi), linewidth=1.5, label=f'$\\zeta$ = {zeta_i}')

ax2.axvline(1, color='k', linestyle='--', alpha=0.5)
ax2.axhline(90, color='k', linestyle=':', alpha=0.5)
ax2.set_xlabel('Frequency Ratio $r = \\omega/\\omega_n$')
ax2.set_ylabel('Phase Angle (degrees)')
ax2.set_title('Frequency Response - Phase')
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
save_fig('frequency_response.pdf')
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{frequency_response.pdf}
\caption{Frequency response function showing magnitude and phase for various damping ratios.}
\end{figure}

\subsection{Resonance Analysis}

\begin{pycode}
# Resonance characteristics
# Peak frequency and amplitude for underdamped systems
zeta_range = np.linspace(0.01, 0.7, 100)
r_peak = np.sqrt(1 - 2*zeta_range**2)
Q_factor = 1 / (2*zeta_range)  # Quality factor

# Maximum amplitude ratio at resonance
X_peak = 1 / (2*zeta_range*np.sqrt(1 - zeta_range**2))

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Peak frequency ratio
axes[0, 0].plot(zeta_range, r_peak, 'b-', linewidth=2)
axes[0, 0].set_xlabel('Damping Ratio $\\zeta$')
axes[0, 0].set_ylabel('Peak Frequency Ratio $r_{peak}$')
axes[0, 0].set_title('Resonance Peak Location')
axes[0, 0].grid(True, alpha=0.3)

# Quality factor
axes[0, 1].semilogy(zeta_range, Q_factor, 'r-', linewidth=2)
axes[0, 1].set_xlabel('Damping Ratio $\\zeta$')
axes[0, 1].set_ylabel('Quality Factor $Q$')
axes[0, 1].set_title('Quality Factor vs Damping')
axes[0, 1].grid(True, alpha=0.3)

# Half-power bandwidth
bandwidth = 2*zeta_range
axes[1, 0].plot(zeta_range, bandwidth, 'g-', linewidth=2)
axes[1, 0].set_xlabel('Damping Ratio $\\zeta$')
axes[1, 0].set_ylabel('Half-Power Bandwidth $\\Delta r$')
axes[1, 0].set_title('Half-Power Bandwidth')
axes[1, 0].grid(True, alpha=0.3)

# Time domain at resonance
omega_f = omega_n  # Forcing at natural frequency
t_res = np.linspace(0, 10, 2000)

def forced_vibration(y, t, m, c, k, F0, omega):
    x, v = y
    F = F0 * np.sin(omega * t)
    return [v, (F - c*v - k*x)/m]

F0 = 10  # N
sol_res = odeint(forced_vibration, [0, 0], t_res, args=(m, c, k, F0, omega_f))

axes[1, 1].plot(t_res, sol_res[:, 0]*1000, 'b-', linewidth=1)
axes[1, 1].set_xlabel('Time (s)')
axes[1, 1].set_ylabel('Displacement (mm)')
axes[1, 1].set_title(f'Response at Resonance ($\\omega = \\omega_n$)')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
save_fig('resonance.pdf')

# Peak amplitude at resonance
X_resonance = F0/k / (2*zeta)
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{resonance.pdf}
\caption{Resonance characteristics: peak location, quality factor, bandwidth, and time response.}
\end{figure}

\section{Transmissibility}

Force transmitted to foundation:
\begin{equation}
TR = \frac{F_T}{F_0} = \sqrt{\frac{1 + (2\zeta r)^2}{(1-r^2)^2 + (2\zeta r)^2}}
\end{equation}

\begin{pycode}
# Transmissibility
r = np.linspace(0.01, 4, 500)
zeta_values = [0.05, 0.1, 0.2, 0.5]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Force transmissibility
for zeta_i in zeta_values:
    TR = np.sqrt((1 + (2*zeta_i*r)**2) / ((1 - r**2)**2 + (2*zeta_i*r)**2))
    ax1.plot(r, TR, linewidth=2, label=f'$\\zeta$ = {zeta_i}')

ax1.axhline(1, color='k', linestyle='--', alpha=0.5)
ax1.axvline(np.sqrt(2), color='r', linestyle=':', alpha=0.5, label='$r = \\sqrt{2}$')
ax1.set_xlabel('Frequency Ratio $r$')
ax1.set_ylabel('Transmissibility $TR$')
ax1.set_title('Force Transmissibility')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_ylim(0, 5)

# Isolation efficiency
isolation = (1 - 1/TR) * 100  # When TR < 1 (r > sqrt(2))
for zeta_i in [0.1, 0.2]:
    TR_i = np.sqrt((1 + (2*zeta_i*r)**2) / ((1 - r**2)**2 + (2*zeta_i*r)**2))
    isolation_i = np.maximum(0, (1 - TR_i) * 100)
    ax2.plot(r, isolation_i, linewidth=2, label=f'$\\zeta$ = {zeta_i}')

ax2.axvline(np.sqrt(2), color='r', linestyle=':', alpha=0.5)
ax2.set_xlabel('Frequency Ratio $r$')
ax2.set_ylabel('Isolation Efficiency (\\%)')
ax2.set_title('Vibration Isolation')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
save_fig('transmissibility.pdf')
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{transmissibility.pdf}
\caption{Force transmissibility and vibration isolation efficiency.}
\end{figure}

\section{Two-DOF System (Modal Analysis)}

\begin{pycode}
# Two-DOF system
m1 = m2 = 1.0  # kg
k1 = k2 = k3 = 100  # N/m

# Mass and stiffness matrices
M = np.array([[m1, 0], [0, m2]])
K = np.array([[k1+k2, -k2], [-k2, k2+k3]])

# Eigenvalue problem
eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(M) @ K)
omega_natural = np.sqrt(eigenvalues)
omega_sorted = np.sort(omega_natural)

# Mode shapes
idx = np.argsort(eigenvalues)
modes = eigenvectors[:, idx]
modes = modes / modes[0, :]  # Normalize to first component

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Mode shape visualization
x = [0, 1, 2]  # Node positions
mode1 = [0, modes[0, 0], modes[1, 0]]
mode2 = [0, modes[0, 1], modes[1, 1]]

axes[0, 0].plot(x, mode1, 'bo-', linewidth=2, markersize=10, label=f'Mode 1: $\\omega$ = {omega_sorted[0]:.2f} rad/s')
axes[0, 0].axhline(0, color='k', linestyle='-', linewidth=0.5)
axes[0, 0].set_xlabel('DOF')
axes[0, 0].set_ylabel('Mode Shape Amplitude')
axes[0, 0].set_title('First Mode Shape')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].plot(x, mode2, 'ro-', linewidth=2, markersize=10, label=f'Mode 2: $\\omega$ = {omega_sorted[1]:.2f} rad/s')
axes[0, 1].axhline(0, color='k', linestyle='-', linewidth=0.5)
axes[0, 1].set_xlabel('DOF')
axes[0, 1].set_ylabel('Mode Shape Amplitude')
axes[0, 1].set_title('Second Mode Shape')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Free response of 2-DOF system
def two_dof_system(y, t, M, K, C):
    x = y[:2]
    v = y[2:]
    a = np.linalg.inv(M) @ (-C @ v - K @ x)
    return list(v) + list(a)

C = 0.1 * K  # Proportional damping
y0_2dof = [0.1, 0, 0, 0]  # Initial conditions
t_2dof = np.linspace(0, 5, 1000)
sol_2dof = odeint(two_dof_system, y0_2dof, t_2dof, args=(M, K, C))

axes[1, 0].plot(t_2dof, sol_2dof[:, 0]*1000, 'b-', linewidth=1.5, label='Mass 1')
axes[1, 0].plot(t_2dof, sol_2dof[:, 1]*1000, 'r-', linewidth=1.5, label='Mass 2')
axes[1, 0].set_xlabel('Time (s)')
axes[1, 0].set_ylabel('Displacement (mm)')
axes[1, 0].set_title('2-DOF Free Response')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# FFT to identify natural frequencies
from scipy.fft import fft, fftfreq
N = len(t_2dof)
dt = t_2dof[1] - t_2dof[0]
yf = fft(sol_2dof[:, 0])
xf = fftfreq(N, dt)

axes[1, 1].plot(xf[:N//2], 2/N * np.abs(yf[:N//2]), 'b-', linewidth=1.5)
axes[1, 1].axvline(omega_sorted[0]/(2*np.pi), color='r', linestyle='--', alpha=0.5)
axes[1, 1].axvline(omega_sorted[1]/(2*np.pi), color='r', linestyle='--', alpha=0.5)
axes[1, 1].set_xlabel('Frequency (Hz)')
axes[1, 1].set_ylabel('Amplitude')
axes[1, 1].set_title('Frequency Content (FFT)')
axes[1, 1].set_xlim(0, 5)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
save_fig('modal_analysis.pdf')
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{modal_analysis.pdf}
\caption{Modal analysis of 2-DOF system: mode shapes, time response, and FFT.}
\end{figure}

Natural frequencies: $\omega_1 = \py{f"{omega_sorted[0]:.2f}"}$ rad/s, $\omega_2 = \py{f"{omega_sorted[1]:.2f}"}$ rad/s

\section{Conclusions}

This analysis demonstrates key aspects of vibration analysis:
\begin{enumerate}
    \item SDOF systems are characterized by natural frequency and damping ratio
    \item Logarithmic decrement provides experimental damping measurement
    \item Resonance occurs near natural frequency with phase shift of 90 degrees
    \item Vibration isolation requires $r > \sqrt{2}$
    \item MDOF systems have multiple natural frequencies and mode shapes
    \item FFT analysis identifies frequency content from time domain data
\end{enumerate}

\end{document}