Mechanical Vibrations
Assignment # 2
MATLAB Code for Viscous Dampers
Submission Date: 20-11-2023
Name: Muhammad Tayyab
Roll No 035
Batch: BSME 2020-24
Submitted To Dr. Zahid Iqbal
MATLAB Code
% Define the system parameters
clear all
clc
m = input('inter mass= '); % Mass
% c = 1; % Damping coefficient
k = input('Spring Constant= '); % Spring constant
zeta=input('Zeta Value= ');
natural_frequency=(k/m)^0.5;
c_c=2*m*natural_frequency;
c=c_c*zeta;
natural_frequency=(k/m)^0.5
x0 = 2; % Initial displacement
x_dot0 = 0; % Initial velocity
% Create a function representing the differential equation
% dxdt = x_dot
% dx_dotdt = -c/m * x_dot - k/m * x
ode = @(t, y) [y(2); -(c/m) * y(2) - (k/m) * y(1)];
% Time span for the simulation
tspan = [0 10]; % Adjust the time span as needed
% Solve the ODE using ode45
[t, y] = ode45(ode, tspan, [x0, x_dot0]);
% Extract displacement, velocity, and acceleration
x = y(:, 1); % Displacement
x_dot = y(:, 2); % Velocity
x_double_dot = -(c/m) * x_dot - (k/m) * x; % Acceleration
% Plot displacement, velocity, and acceleration
figure;
%subplot(3,1,1);
plot(t, x, 'g-', 'LineWidth', 2, 'DisplayName', 'Displacement');
title('Displacement vs Time');
xlabel('Time');
ylabel('Displacement');
hold on
%plot(3,1,2);
plot(t, x_dot, 'r-', 'LineWidth', 2, 'DisplayName', 'Velocity');
title('Velocity vs Time');
xlabel('Time');
ylabel('Velocity');
%subplot(3,1,3);
hold on
plot(t, x_double_dot, 'b-', 'LineWidth', 2, 'DisplayName', 'acceleration');
title(' acceleration, Velocity and displacement vs Time');
legend('show');
xlabel('Time');
ylabel('Acceleration, Velocity and Displacement');
Undamped
m=2 kg
k=2 N/m
Zeta=0
Underdamped
m=1 kg
k=2 N/m
Zeta=0.3
Critically Damped
m=3 kg
k=2 N/m
Zeta=1
Overdamped
m=2 kg
k=3 N/m
Zeta=2
Mechanical Vibrations
Assignment # 3
Coulomb and Hysteresis Dampers
Submission Date: 20-11-2023
Name: Muhammad Tayyab
Roll No 035
Batch: BSME 2020-24
Submitted To Dr. Zahid Iqbal
MATLAB Code for Hysteretic Damping
% Define the system parameters
clear all
clc
m = input('inter mass= '); % Mass
% c = 1; % Damping coefficient
k = input('Spring Constant= '); % Spring constant
zeta=input('Zeta Value= ');
natural_frequency=(k/m)^0.5;
c_c=2*m*natural_frequency;
c=c_c*zeta;
natural_frequency=(k/m)^0.5
x0 = 2; % Initial displacement
x_dot0 = 0; % Initial velocity
% Create a function representing the differential equation
% dxdt = x_dot
% dx_dotdt = -c/m * x_dot - k/m * x
ode = @(t, y) [y(2); -(c/m) * y(2) - (k/m) * y(1)];
% Time span for the simulation
tspan = [0 10]; % Adjust the time span as needed
% Solve the ODE using ode45
[t, y] = ode45(ode, tspan, [x0, x_dot0]);
% Extract displacement, velocity, and acceleration
x = y(:, 1); % Displacement
x_dot = y(:, 2); % Velocity
x_double_dot = -(c/m) * x_dot - (k/m) * x; % Acceleration
% Plot displacement, velocity, and acceleration
figure;
%subplot(3,1,1);
plot(t, x, 'g-', 'LineWidth', 2, 'DisplayName', 'Displacement');
title('Displacement vs Time');
xlabel('Time');
ylabel('Displacement');
hold on
%plot(3,1,2);
plot(t, x_dot, 'r-', 'LineWidth', 2, 'DisplayName', 'Velocity');
title('Velocity vs Time');
xlabel('Time');
ylabel('Velocity');
%subplot(3,1,3);
hold on
plot(t, x_double_dot, 'b-', 'LineWidth', 2, 'DisplayName', 'acceleration');
title(' acceleration, Velocity and displacement vs Time');
legend('show');
xlabel('Time');
ylabel('Acceleration, Velocity and Displacement');
%plot(t, undamped_response, 'r-', 'LineWidth', 2, 'DisplayName', 'Undamped');
Undamped
m=3 kg
k=2 N/m
Zeta=0
Underdamped
m=5 kg
k=30 N/m
Zeta=0.35
Critically Damped
m=10 kg
k=10 N/m
Zeta=1
Overdamped
m=5 kg
k=7 N/m
z=1.5
MATLAB Code for Coulomb Damping
% Parameters
clc
m = input('mass= '); % mass
mu = input('coefficient of friction= '); % damping coefficient
g = 9.81; % acceleration due to gravity
k = input('coefficient of spring constant= '); % spring constant
% Define the differential equation
ode = @(t, x) [x(2); (-mu * m * g * sign(x(2)) - k * x(1)) / m];
% Time span
tspan = [0 10];
% Initial conditions
x0 = [5; 0]; % initial displacement and velocity
% Solve the differential equation using ode45
[t, x] = ode45(ode, tspan, x0);
% Plot the results
figure;
subplot(2, 1, 1);
plot(t, x(:, 1));
xlabel('Time (t)');
ylabel('Displacement (x)');
subplot(2, 1, 2);
plot(t, x(:, 2));
xlabel('Time (t)');
ylabel('Velocity');
% Add title and labels
% suptitle('Solution of the Nonlinear ODE with Damping');
Plots
Plot 1
m=1
k=10
Plot 2
m=5
k=10
Plot 3
m=1
k=10
Mechanical Vibrations
Assignment # 4
Magnification Factor
Submission Date: 20-11-2023
Name: Muhammad Tayyab
Roll No 035
Batch: BSME 2020-24
Submitted To Dr. Zahid Iqbal
MATLAB Code for Magnification Factor
% System parameters
m = 1; % mass
%c = 0.1; % damping coefficient
k = 1; % stiffness
F0 = 1; % amplitude of the force
delta=F0/k;
% Frequency range
omega_range = 0: 0.1: 10; % adjust the range as needed
Cc=2*sqrt(k*m);
% Damping ratios to be considered
zeta_values = [0.1, 0.2, 0.3, 0.4, 0.5,1, 1.2, 1.5];
% Initialize arrays to store results
max_amplitudes = zeros(length(zeta_values), length(omega_range));
% Loop over each damping ratio
for i = 1:length(zeta_values)
zeta = zeta_values(i);
c=zeta*Cc
% Loop over each forcing frequency
for j = 1:length(omega_range)
omega = omega_range(j);
% Solve the differential equation
tspan = [0, 20];
initial_conditions = [0; 0]; % initial position and velocity
[t, x] = ode45(@(t, x) [(x(2)); (F0 * sin(omega * t) - c * x(2) - k *
x(1)) / m - 2 * zeta * sqrt(k/m) * x(2)], tspan, initial_conditions);
% Find maximum amplitude
max_amplitudes(i, j) = max(x(:, 1));
end
end
% Plot results
figure;
plot(omega_range, max_amplitudes/delta);
xlabel('Frequency Ratio (\omega/\omega_n)');
ylabel('Maximum Amplitude');
title('Force Vibration with Viscous Damping');
legend(strcat('\zeta = ', cellstr(num2str(zeta_values'))));
grid on;
Results
MATLAB Code for Phase Angle
% Define the range of r values
r = linspace(0, 5, 100);
% Define different values for zeta
zeta_values = [0.1,0.2,0.3,0.4 0.5, 1, 2];
% Plot phi vs r for each zeta value
figure;
hold on;
for i = 1:length(zeta_values)
zeta = zeta_values(i);
phi_rad =atan((2 * zeta * r) ./ (1 - r.^2));
% Convert phi from radians to degrees
phi_deg = rad2deg(phi_rad);
% Ensure phi is in the range [0, 180]
phi_deg(phi_deg < 0) = phi_deg(phi_deg < 0) + 180;
plot(r, phi_deg, 'DisplayName', ['\zeta = ' num2str(zeta)]);
end
hold off;
% Add labels and legend
xlabel('r');
ylabel('\phi (degrees)');
title('Plot of \phi vs r for different \zeta values');
legend('show');
Results
Mechanical Vibrations
Assignment # 5
Phase Angle vs Frequency Ratio
Submission Date: 27-11-2023
Name: Muhammad Tayyab
Roll No 035
Batch: BSME 2020-24
Submitted To Dr. Zahid Iqbal
MATLAB Code for Displacement Transmissibility
% System parameters
m = 1; % mass
%c = 0.1; % damping coefficient
k = 1; % stiffness
F0 = 1; % amplitude of the force
Y=5;
% Frequency range
omega_range = 0: 0.1: 10; % adjust the range as needed
Cc=2*sqrt(k*m)
% Damping ratios to be considered
zeta_values = [0.1, 0.2, 0.3, 0.4, 0.5,1, 1.2, 1.5];
% Initialize arrays to store results
max_amplitudes = zeros(length(zeta_values), length(omega_range));
% Loop over each damping ratio
for i = 1:length(zeta_values)
zeta = zeta_values(i);
c=zeta*Cc
% Loop over each forcing frequency
for j = 1:length(omega_range)
omega = omega_range(j);
% Solve the differential equation
tspan = [0, 20];
initial_conditions = [0; 0]; % initial position and velocity
[t, x] = ode45(@(t, x) [x(2); -(c/m)*(x(2) - Y*omega*cos(omega*t)) -
(k/m)*(x(1) - Y*sin(omega*t))], tspan, initial_conditions);
% Find maximum amplitude
max_amplitudes(i, j) = max(x(:, 1));
end
end
% Plot results
figure;
plot(omega_range, max_amplitudes/Y);
xlabel('Frequency Ratio (\omega/\omega_n)');
ylabel('Maximum Amplitude');
title('Force Vibration with Viscous Damping');
legend(strcat('\zeta = ', cellstr(num2str(zeta_values'))));
grid on;
Results
Phase Angle Matlab Code
% Define the range of r values
r = linspace(0, 10, 100);
% Define different values for zeta
zeta_values = [0.1,0.2,0.3,0.4 0.5, 1, 2];
% Plot phi vs r for each zeta value
figure;
hold on;
for i = 1:length(zeta_values)
zeta = zeta_values(i);
phi_rad = atan((2*zeta*r.^3)./(1+(4*zeta^2-1)*r.^2));
% Convert phi from radians to degrees
phi_deg = rad2deg(phi_rad);
% Ensure phi is in the range [0, 180]
phi_deg(phi_deg < 0) = phi_deg(phi_deg < 0) + 180;
plot(r, phi_deg, 'DisplayName', ['\zeta = ' num2str(zeta)]);
end
hold off;
% Add labels and legend
xlabel('r');
ylabel('\phi (degrees)');
title('Plot of \phi vs r for different \zeta values');
legend('show');
Result
Mechanical Vibrations
Assignment # 6
Data From Accelerometer
Submission Date: 27-11-2023
Name: Muhammad Tayyab
Roll No 035
Batch: BSME 2020-24
Submitted To Dr. Zahid Iqbal
Assignment 6
Assignment 6 focuses on the collection of data from viscous dampers, hysteresis
dampers, and Coulomb dampers. The provided code is integral to the completion
of this assignment, with further details presented in the CEP report.
Arduino Code for accelerometer adxl 345:
#include <Wire.h>
#include <Adafruit_Sensor.h>
#include <Adafruit_ADXL345_U.h>
#define ADXL345_ADDRESS_1 (0x53) // I2C address for the first ADXL345 sensor
#define ADXL345_ADDRESS_2 (0x1D) // I2C address for the second ADXL345 sensor
Adafruit_ADXL345_Unified accel1 = Adafruit_ADXL345_Unified(ADXL345_ADDRESS_1);
Adafruit_ADXL345_Unified accel2 = Adafruit_ADXL345_Unified(ADXL345_ADDRESS_2);
const int arraySize = 100; // Adjust the size of the array as needed
int dummy[arraySize];
float xValues1[arraySize], yValues1[arraySize], zValues1[arraySize];
float xValues2[arraySize], yValues2[arraySize], zValues2[arraySize];
unsigned long times[arraySize];
void setup() {
[Link](2000000);
// Initialize the first accelerometer
if (![Link]()) {
[Link]("Could not find a valid ADXL345 sensor at address 0x53, check wiring!");
while (1);
}
// Set the range to +/- 16g
[Link](ADXL345_RANGE_16_G);
// Initialize the second accelerometer
if (![Link]()) {
[Link]("Could not find a valid ADXL345 sensor at address 0x1D, check wiring!");
while (1);
}
// Set the range to +/- 16g
[Link](ADXL345_RANGE_16_G);
}
void loop() {
for (int i = 0; i < arraySize; i++)
{
// Read data from the first accelerometer
sensors_event_t event1;
[Link](&event1);
xValues1[i] = [Link].x;
yValues1[i] = [Link].y;
zValues1[i] = [Link].z;
// Read data from the second accelerometer
sensors_event_t event2;
[Link](&event2);
xValues2[i] = [Link].x;
yValues2[i] = [Link].y;
zValues2[i] = [Link].z;
times[i] = millis();
}
// Store time and dummy variable
for (int i = 0; i < arraySize; i++)
{
dummy[i] = i;
// Print data
[Link]("Dummy: ");
[Link](dummy[i]);
[Link](", Accel1 - X: ");
[Link](xValues1[i]);
[Link](", Y: ");
[Link](yValues1[i]);
[Link](", Z: ");
[Link](zValues1[i]);
[Link](", Accel2 - X: ");
[Link](xValues2[i]);
[Link](", Y: ");
[Link](yValues2[i]);
[Link](", Z: ");
[Link](zValues2[i]);
[Link](", Time: ");
[Link](times[i]);
// Add a delay if needed
delay(0.1); // Adjust the delay as needed
}
}