ME3255_group9
##1) Analytically solve for the shape of the beam if q(x)=cst, P=0, and ω=0 and create a function called shape_simple_support.m that returns the displacement w(x) given q and x
function w =shape_simple_support(x, q)
%The Differential equation for the shape of the beam was integrated to
%solve for the equation of the beam. The input is the postition on the beam
%and the constant distributed load and the output is the shape of the beam.
%The plot shows the maximum deflection that occurs at loads ranging from 1 to 1000
E=70e9;
b=0.1;
h=0.01;
I=(b*h^3)/12;
w=q*((x^4)/(24*E*I))-((q*(x^3))/(12*E*I))+((q*x)/(24*E*I));
q1=1:1:1000;
for n= 1:1:1000
d= (n/(384*E*I))-(n/(96*E*I))+(n/(48*E*I));
D(n)=d;
end
plot(D,q1,'.');
xlabel 'Maximum Deflection'
ylabel 'Distributed Load'
title('q vs Maximum Deflection')
#1a
As the load increases the maximum deflection of the beam increases proportionally and is a linear trend.
#1b Use a Monte Carlo model to determine the mean and standard deviation for the maximum deflection δx if b and h are normally distributed random variables with 0.1 % standard deviations at q=50 N/m.
function [mean1, stand_dev ]=montecarlo_support(N)
%Uses monte carlo model to input number of randomn varaibles to generate varying variables b and
%h so that they have a standard deviation
%The maximum deflection is then solved for and the mean and standard deviation of that is returned.
E=70e9;
q=50;
for i=1:N
b=normrnd(0.1,0.0001);
h=normrnd(0.01,0.00001);
I=(b*h^3)/12;
D= (q/(384*E*I))-(q/(96*E*I))+(q/(48*E*I));
d(i)=D;
end
mean1=mean(d);
stand_dev=std(d);
end
This function takes the number of random variables as an inout and computes a b and h with standard deviation of 0.1% and those are used to calculate the maximum deflection and the mean and standard deviation of the maximum deflections are computed.
#2) Now use the central difference approximation to set up a system of equations for the beam for q(x)=cst, P=0, and ω = 0. Use the boundary conditions with a numerical differentiation to determine the values of the end points.
#2 a-c In the following script, the matrices were solved by hand using the central difference formulas, the fourth and second derivative were used (found in the textbook page 529) for segments of 6, 10 and 20, the max deflection was then calulcated using for loops for each segment at each resepective distributed load. The values of the max deflcetion were then plotted vs the distributed load.
seg6 = [ 5 -4 1 0 0
-4 6 -4 1 0
1 -4 6 -4 1
0 1 -4 6 -4
0 0 1 -4 5];
seg10= [5 -4 1 0 0 0 0 0 0
-4 6 -4 1 0 0 0 0 0
1 -4 6 -4 1 0 0 0 0
0 1 -4 6 -4 1 0 0 0
0 0 1 -4 6 -4 1 0 0
0 0 0 1 -4 6 -4 1 0
0 0 0 0 1 -4 6 -4 1
0 0 0 0 0 1 -4 6 -4
0 0 0 0 0 0 1 -4 5];
seg20= [5 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-4 6 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 -4 6 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 -4 6 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 -4 6 -4 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 -4 6 -4 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 -4 6 -4 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 -4 6 -4 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 -4 6 -4 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 -4 6 -4 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 -4 6 -4 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 -4 6 -4 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 -4 6 -4 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 -4 6 -4 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 -4 6 -4 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 6 -4 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 6 -4 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 6 -4
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 5];
b=0.1;
h=0.01;
I=b*(h^3)/12;
E=70e9;
q=[1,10,20,30,50];
n=[6,10,20];
for i=1:5
for N = 1:3
b=(((1/n(N))^4)*q(i))/(E*I);
if n(N)==6
B=b*ones((n(N)-1),1);
W=seg6\B;
dx(i,N)=max(W);
elseif n(N)==10
B=b*ones((n(N)-1),1);
W=seg10\B;
dx(i,N)=max(W);
elseif n(N)==20
B=b*ones((n(N)-1),1);
W=seg20\B;
dx(i,N)=max(W);
end
end
end
plot(dx,q,'o')
xlabel('dx')
ylabel('q')
title('Load vs Max Deflection')
2d
As the number of segments increase the slope of the linear relationship between load and maximum deflection decreases.
#3) Now set up the system of equations using a central difference method if P>0 and ω = 0
3 a-c
function centraldifference2(n)
%Takes the number of segments as the input argument and then runs a for
%loop of the P value and a for loop of the q value to find the maximum deflection
%with those given parameters, the calculated max deflcetion valus are then
%plotted with respect to the given distributed loads.
%Within the for loop the matrix is solved for using the central difference
%formulas and patterns that it follows and then the backslash operator is
%used
b=0.1;
h=0.01;
I=b*(h^3)/12;
E=70e9;
P=[0,100,200,300];
q=[1,10,20,30,50];
for p=1:4
for i=1:5
D=-4-((P(p)*((1/n)^2))/(E*I));
e=ones((n-2),1)*D;
F= 6+(P(p)*(2*(1/n)^2)/(E*I));
f=ones(n-1,1)*F;
g=ones((n-3),1);
A=diag(e,-1)+ diag(f)+diag(e,1)+diag(g,2)+diag(g,-2);
A(1,1)=5+(P(p)*(2*(1/n)^2)/(E*I));
A((n-1),(n-1))=5+(P(p)*(2*(1/n)^2)/(E*I));
b=(((1/n)^4)*q(i))/(E*I);
if n==6
B=b*ones((n-1),1);
W=A\B;
dx(i)=max(W);
elseif n==10
B=b*ones((n-1),1);
W=A\B;
dx(i)=max(W);
elseif n==20
B=b*ones((n-1),1);
W=A\B;
dx(i)=max(W);
end
end
DX(p,:)=dx;
plot(q,DX)
xlabel('q (force/length)')
ylabel('dx (max deflection)')
title('Load vs Max Deflection')
legend('P=0','P=100','P=200','P=300')
end
end
3d
An input Argument of 6 yields this graph:
An input of 10 yields this graph:
An input of 20 yields this graph:
The graphs above show that as the transverse loading of the beam increases the maximum deflction decreases which makes sense when a free body diagram is drawn.
#4) Now set up an eigenvalue problem to solve for the natural frequencies of the simply supported beam if P=0 and q=0.
% 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')
function [ w1 ] = Coefficients( n,EIpa,L )
%Coefficients calculates the values of the w matrix where the values
%represent the coefficients of each node in the computational domain. The
%coefficient matrix is of the form: Ax = b
%The problem is posed as a Ax = (lambda)x = 0 so this makes sense to write
%the matrix as a coefficient matrix
% Change in x for central difference formula for 6 segments, dx6
dx = L/n;
% Precalculate denominator in central difference equation
dx4 = dx^4;
% Precalculate coefficients' constant (see derivation)
K = (EIpa)/dx4; % Central points coefficient
Kf = (EIpa)/(2*dx); % boundary points coefficients (i = 2, i = n)
% Preallocate w matrix
w1 = zeros(n+1,n+1);
% Generate coefficient matrix for computational domain
% Boundary nodes w(1) and w(n+1) should equal zero
for i = 1:n+1
if i == 1
w1(i,i) = 6*K;
w1(i,i+1) = -4*K;
w1(i,i+2) = K;
elseif i == 2
w1(i,i-1) = -4*K;
w1(i,i) = 6*K;
w1(i,i+1) = -4*K;
w1(i,i+2) = K;
elseif i == n
w1(i,i-2) = K;
w1(i,i-1) = -4*K;
w1(i,i) = 6*K;
w1(i,i+1) = -4*K;
elseif i == n+1
w1(i,i-2) = K;
w1(i,i-1) = -4*K;
w1(i,i) = 6*K;
else
w1(i,i-2) = K;
w1(i,i-1) = -4*K;
w1(i,i) = 6*K;
w1(i,i+1) = -4*K;
w1(i,i+2) = K;
end
end
%% Code Debug: NOT NECESSARY (CONFIRMED VIA EMAIL)
% w2 = w1;
% % If Boundary nodes are different:
% % Written using forward and backward difference equations
% % First Node
% i = 1;
% w2(i,i) = -3*Kf;
% w2(i,i+1) = 4*Kf;
% w2(i,i+2) = -Kf;
% % Last Node
% i = n+1;
% w2(i,i) = 3*Kf;
% w2(i,i-1) = -4*Kf;
% w2(i,i-2) = Kf;
%
end
function [ W ] = LHS ( w,EIpa,n,omega,t )
%LHS rewrites the coefficient matrix w to the form denoted in
%equation 12 in the derivation, after the e^iwt coefficient was expanded to
%cos(wt)-isin(wt)
% Rewrite coefficient matrix in to follow form of equation 4 in derivation
wp = w./EIpa;
% Preallocate W for speed
W = wp;
% Precalculate coefficient to w term in equation 4
A = ((omega)^2)/EIpa;
% Rewrite eigenvalue problem such that the w matrix becomes (Aw - lambda w)
% which is the same as the terms in paranthesis in equation 4 (see issue in
% our project)
for i = 1:n+1
W(i,i) = wp(i,i) - A;
end
% Element by element multiplication of term in paranthesis by cos(omega*t)
W = W.*cos(omega*t);
end
function [ RHS ] = RHS( w,q,E,I )
%RHS writes the RHS of the equation (w(x))cos(wt) = q/EI
dimen = length(w);
RHS = zeros(dimen,1);
for i = 1:dimen
RHS(i) = q/(E*I);
end
end
Plot of the beam for the first natural frequency
Plot of the beam for the Second natural Frequency
Plot of the beam for the third natural frequency
#5 Extra Credit
% ME 3255 Computational Mechanics
% Peter Joseph Damian, Group 9
% Part 4 Rewritten for bonus problem
clear; clc;
%% Please see issue "Derivation" on github for the derivations of these equations
%% Equation 16: Setting up eigenvalue problem to solve for lowest natural frequency
%% Constants and Properties
% 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
%% 20 Segment System of Equations
% Number of Segments
n20 = 20;
% Write Values of P
P = linspace(-1000,1000,2000);
omega = zeros(length(P),1);
for i = 1:length(P)
Q = P(i);
% Write w coefficient matrix (see above for details)
[ w20 ] = Coefficients( n20,EIpa,L,Q,E,I );
% Calculate eigenvalues of system
Lambda20 = eig(w20);
% Calculate natural frequencies of system
omega20 = Lambda20.^(1/2);
% Store lowest natural frequency for use outside of loop
omega(i) = omega20(1);
end
%% Plotting
setdefaults
figure(1)
plot(P,omega)
title('Lowest Nat Freq vs Load P')
xlabel('Load P,[N]')
ylabel('Omega, [1/s]')
function [ w1 ] = Coefficients( n,EIpa,L,P,E,I )
%Coefficients calculates the values of the w matrix where the values
%represent the coefficients of each node in the computational domain. The
%coefficient matrix is of the form: Ax = b
%The problem is posed as a Ax = (lambda)x = 0 so this makes sense to write
%the matrix as a coefficient matrix
% Change in x for central difference formula for 6 segments, dx6
dx = L/n;
% Precalculate denominator in central difference equation
dx2 = dx^2;
dx4 = dx^4;
% Precalculate coefficients' constant (see derivation)
K = (EIpa)/dx4; % Central points coefficient
% Coefficient for Second Derivative term (new term in this equation)
SD = EIpa*(P/(E*I*dx2));
% Preallocate w matrix
w1 = zeros(n+1,n+1);
% Generate coefficient matrix for computational domain
% Boundary nodes w(1) and w(n+1) should equal zero
for i = 1:n+1
if i == 1
w1(i,i) = 6*K+2*SD;
w1(i,i+1) = -4*K-SD;
w1(i,i+2) = K;
elseif i == 2
w1(i,i-1) = -4*K-SD;
w1(i,i) = 6*K+2*SD;
w1(i,i+1) = -4*K-SD;
w1(i,i+2) = K;
elseif i == n
w1(i,i-2) = K;
w1(i,i-1) = -4*K-SD;
w1(i,i) = 6*K+2*SD;
w1(i,i+1) = -4*K-SD;
elseif i == n+1
w1(i,i-2) = K;
w1(i,i-1) = -4*K-SD;
w1(i,i) = 6*K+2*SD;
else
w1(i,i-2) = K;
w1(i,i-1) = -4*K;
w1(i,i) = 6*K;
w1(i,i+1) = -4*K;
w1(i,i+2) = K;
end
end
%% Code Debug: NOT NECESSARY (CONFIRMED VIA EMAIL)
% w2 = w1;
% % If Boundary nodes are different:
% % Written using forward and backward difference equations
% % First Node
% i = 1;
% w2(i,i) = -3*Kf;
% w2(i,i+1) = 4*Kf;
% w2(i,i+2) = -Kf;
% % Last Node
% i = n+1;
% w2(i,i) = 3*Kf;
% w2(i,i-1) = -4*Kf;
% w2(i,i-2) = Kf;
%
end
