Embedded Files
%Created by Jerry Tustaniwskyj
%Developed by Team 10
clear all;
%
global E I L L2 % save parameters that are used by other functions
E=2.85e9; % elastic modulus (Pa)
ro=0.0254*6.625/2; % outer radius
ri=0.0254*6.065/2; % inner radius
I=pi*(ro^4-ri^4)/4; % area moment of inertia
L=10; % length of spar
L2=5*sqrt(2); % distance to gut wire fastening location
%L2=5;
x0=[-140, .05]; % initial guesses for R1(0) and y'(0)
% the Matlab function fminsearch is used to find correct R1(0) and y'(0)
options = optimset('Tolfun',1.e-8,'TolX',1.e-8,'MaxFunEvals',1500);
[x, fval, exitflag, output]=fminsearch(@beam2,x0,options);
np=201; % number of points to plot from 0 to L2
x2=linspace (0,L2,np); % points from 0 to L2 (needed to plot)
p0=[x(1); 0; x(2); 0]; % correct initial conditions
[x2,p] = ode45(@gaccl,x2,p0); % integrate 0 to L2 to obtain points to plot
x3=linspace (L2,L,np); % points from L2 to L (needed to plot)
reaction=-L*frc(L); %
V=reaction-x(1); % calculate change in shear at L2
% initial conditions for second part
p0=[p(np,1)+V; p(np,2); p(np,3); p(np,4)];
[x3,p2] = ode45(@gaccl,x3,p0); % integrate L2 to L to obtain points to plot
% combine arrays for 1 plot, x2 & x3 into xp and p & p2 into pp
for i=1:np
xp(i)=x2(i);
xp(i+np)=x3(i);
pp(i,1)=p(i,1);
pp(i+np,1)=p2(i,1);
pp(i,2)=p(i,2);
pp(i+np,2)=p2(i,2);
pp(i,3)=p(i,3);
pp(i+np,3)=p2(i,3);
pp(i,4)=p(i,4);
pp(i+np,4)=p2(i,4);
end
%
figure(1);
plot (xp,pp(:,1));
grid on;
xlabel ('Length along spar (m)');
ylabel ('Shear (N)');
title('Shear vs. Spar Length');
%
figure(2);
plot (xp,pp(:,2));
grid on;
xlabel ('Length along spar (m)');
ylabel ('Moment (Nm)');
title('Moment vs. Spar Length');
%
figure(3);
plot (xp,pp(:,3));
grid on;
xlabel ('Length along spar (m)');
ylabel ('Slope (rad)');
title('Slope vs. Spar Length');
%
figure(4);
plot (xp,pp(:,4));
grid on;
xlabel ('Length along spar (m)');
ylabel ('Displacement (m)');
title('Deflection vs. Spar Length');
%
figure (5);
str=ro*pp(:,2)/I/1.e6; % calculate stress from moment
plot (xp,str);
grid on;
xlabel ('Length along spar (m)');
ylabel ('Stress (MPa)');
title('Stress vs. Spar Length');
%
function pend=beam2(roots)
% roots transferred in are provided by root finder
global E I L L2
np=3; % number of integration points
x=linspace (0,L2,np); % define number of integration points from 0 to L2
p0=[roots(1); 0; roots(2); 0]; % these are initial conditions
[x,p] = ode45(@gaccl,x,p0); % integrate from 0 to L2
psave=p(np,4); % save displacement at L2 for root finder
reaction=-L*frc(L);
V=reaction-roots(1); % calculate R2 to modify shear
x2=linspace (L2,L,np); % define number of integration points from 0 to L2
% initial conditions for second part
p0=[p(np,1)+V; p(np,2); p(np,3); p(np,4)];
[x2,p2] = ode45(@gaccl,x2,p0); % integrate from L2 to L
% The function fminsearch uses the value pend to find the correct roots
% that will make pend minimum (zero in our case). psave is multiplied
% EI to make the 3 zeros converge better
pend=sqrt((E*I*psave)^2+p2(np,1)^2+p2(np,2)^2); % should be zero
end
%
function pdot = gaccl(x,p)
% function used by ode45 to integrate beam equations
global E I L L2
% p(1) = shear
% p(2) = moment
% p(3) = slope
% p(4) = displacement
% get force per unit length at x
w = frc(x);
pdot = zeros(4,1); % initialize array
pdot(1)=w; % derivative of shear
pdot(2)=p(1); % derivative of moment
pdot(3)=p(2)/E/I; % derivative of slope
pdot(4)=p(3); % derivative of displacement
end
%
function w=frc(x)
% this function defines the force per length at x for the spar
w=57.45; % assumed constant
end
Page updated
Report abuse