Part4_Main_Script.m

% ME 3255 Computational Mechanics
% Peter Joseph Damian, Group 9
% Part 4
clear; clc;

%% Please see issue "Derivation" on github for the derivations of these equations
%% Equation 7: Setting up eigenvalue problem to solve for natural frequencies
%% Equation 12: writing equations as an Ax = b problem for nontrivial solution of Ax = 0

% Material Properties
E = 70*(10^6);          % Young's modulus, Pa
p = 2700;               % Density, kg/m^3

% Geometry
b = 0.1;                % base width, m
h = 0.01;               % height, m
L = 1;                  % bar length, m
A = b*h;                % area, m^2
I = (b*(h^3))/12;       % second moment of inertia, m^4

% Frequency squared coefficient
EIpa = (E*I)/(p*A);     % Inverse coefficient of omega^2

%% 6 Segment System of Equations

% Number of Segments
n = 6;

% Write w matrix coefficients in form Aw = lambda w = 0 where lambda is the
% natural frequency squared and the other coefficients, rho, A, E and I are
% shifted to the other side of the Aw = lambda w = 0 equation
[ w6 ] = Coefficients( n,EIpa,L );

% Calculate eigenvalues of system
Lambda6 = eig(w6);

% Calculate natural frequencies of system
omega6 = Lambda6.^(1/2);

%% 10 Segment System of Equations

% Number of Segments
n10 = 10;

% Write w coefficient matrix (see above for details)
[ w10 ] = Coefficients( n10,EIpa,L );

% Calculate eigenvalues of system
Lambda10 = eig(w10);

% Calculate natural frequencies of system
omega10 = Lambda10.^(1/2);

%% 20 Segment System of Equations

% Number of Segments
n20 = 20;

% Write w coefficient matrix (see above for details)
[ w20 ] = Coefficients( n20,EIpa,L );

% Calculate eigenvalues of system
Lambda20 = eig(w20);

% Calculate natural frequencies of system
omega20 = Lambda20.^(1/2);


%% Beam Shape Plotting for the first 3 natural frequencies

% Selection of evalutation time
t = 10;

% Set value of q (could write a function for this and loop through code
q = 100;

% 6 Segment System of Equation
omega61 = omega6(1);
[ W61 ] = LHS ( w6,EIpa,n,omega61,t );
[ RHS61 ] = RHS( w6,q,E,I );
W61solv = W61\RHS61;

omega62 = omega6(2);
[ W62 ] = LHS ( w6,EIpa,n,omega62,t );
[ RHS62 ] = RHS( w6,q,E,I );
W62solv = W62\RHS62;

omega63 = omega6(3);
[ W63 ] = LHS ( w6,EIpa,n,omega63,t );
[ RHS63 ] = RHS( w6,q,E,I );
W63solv = W63\RHS63;

% 10 Segment System of Equation
omega101 = omega10(1);
[ W101 ] = LHS ( w10,EIpa,n,omega101,t );
[ RHS101 ] = RHS( w10,q,E,I );
W101solv = W101\RHS101;

omega102 = omega10(2);
[ W102 ] = LHS ( w10,EIpa,n,omega102,t );
[ RHS102 ] = RHS( w10,q,E,I );
W102solv = W102\RHS102;

omega103 = omega10(3);
[ W103 ] = LHS ( w10,EIpa,n,omega103,t );
[ RHS103 ] = RHS( w10,q,E,I );
W103solv = W103\RHS103;

% 20 Segment System of Equation
omega201 = omega20(1);
[ W201 ] = LHS ( w20,EIpa,n,omega201,t );
[ RHS201 ] = RHS( w20,q,E,I );
W201solv = W201\RHS201;

omega202 = omega20(2);
[ W202 ] = LHS ( w20,EIpa,n,omega202,t );
[ RHS202 ] = RHS( w20,q,E,I );
W202solv = W202\RHS202;

omega203 = omega20(3);
[ W203 ] = LHS ( w20,EIpa,n,omega203,t );
[ RHS203 ] = RHS( w20,q,E,I );
W203solv = W203\RHS203;

%% Plotting
setdefaults

% 6 Segments
dx = L/(n+1);
x6 = linspace(dx,L-dx,n+1);

figure(1)
plot(x6,W61solv)
xlabel('x distance, m')
ylabel('Deflection')
figure(2)
plot(x6,W62solv)
xlabel('x distance, m')
ylabel('Deflection')
figure(3)
plot(x6,W63solv)
xlabel('x distance, m')
ylabel('Deflection')

% 10 Segments
dx10 = L/(n10+1);
x10 = linspace(dx10,L-dx10,n10+1);

figure(4)
plot(x10,W101solv)
xlabel('x distance, m')
ylabel('Deflection')
figure(5)
plot(x10,W102solv)
xlabel('x distance, m')
ylabel('Deflection')
figure(6)
plot(x10,W103solv)
xlabel('x distance, m')
ylabel('Deflection')

% 20 Segments
dx20 = L/(n20+1);
x20 = linspace(dx20,L-dx20,n20+1);

figure(7)
plot(x20,W201solv)
xlabel('x distance, m')
ylabel('Deflection')
figure(8)
plot(x20,W202solv)
xlabel('x distance, m')
ylabel('Deflection')
figure(9)
plot(x20,W203solv)
xlabel('x distance, m')
ylabel('Deflection')