root/OptimizingMATLABCode/Truss/trussCantilever.m @ 11
| 10 | anderm8 | function [Yr,bars,L,N,H]=trussCantilever(N,Aval)
|
|
% Yr = trussCantilever(N,Aval);
|
|||
%
|
|||
% Usage:
|
|||
% Y=trussCantilever(N,Aval)
|
|||
% [Yr,bars,L,N,H]=trussCantilever(N,Aval)
|
|||
%
|
|||
% Inputs:
|
|||
%
|
|||
% N: number of horizontal elements (recommend range of [1,400])
|
|||
%
|
|||
% Aval: cross section area for each rod (recommend range of [10 1000])
|
|||
%
|
|||
% Outputs:
|
|||
%
|
|||
% Yr:
|
|||
%
|
|||
%
|
|||
% Copyright 2015 The MathWorks, Inc
|
|||
%% Defaults
|
|||
if nargin < 2
|
|||
% Ten horizontal truss elements in bottom cord
|
|||
N = 10;
|
|||
% Cross sectional area of 100
|
|||
Aval = 100;
|
|||
end
|
|||
%% Truss Parameters
|
|||
% Physical parameters of the material and truss
|
|||
% Horizontal Length of the truss
|
|||
L = 1;
|
|||
% Vertical Height of the truss
|
|||
H = 0.025;
|
|||
% Density of truss bar material
|
|||
rho = 1;
|
|||
% Modulus of Elasticity of truss bar material
|
|||
Eval = 1e-1;
|
|||
% Rayleigh, "alpha,beta" damping coefficients
|
|||
dampingCoeffAlpha = 0.01;
|
|||
dampingCoeffBeta = 0.01;
|
|||
%% Force Parameters
|
|||
%
|
|||
% The following downward force is applied to the truss
|
|||
%
|
|||
% max(0,(applyTime-t)*downwardForceMag/applyTime)*sin(2*pi*freq*t)
|
|||
%
|
|||
% The frequency of the force is 1. The magnitude is linearly lowered from
|
|||
% _downwardForceMag_ to zero over the interval |[0,applyTime]|. The
|
|||
% simulation time is simTime. Because of the damping the vibration will die
|
|||
% out.
|
|||
downwardForceMag = 0.02;
|
|||
freq = 1;
|
|||
applyTime = 20;
|
|||
simTime = 100;
|
|||
%% Map reduced Dofs to actual Dofs
|
|||
% The truss is fixed to the wall on the left hand side. Therefore the two
|
|||
% nodes there are not mobile and can be eliminated. Create a map between
|
|||
% the "reduced" degrees of freedom and the actual degrees of freedom
|
|||
numDofs = 2*2*(N+1)-2; % -2 because the cantilever is pointed upwards at the end
|
|||
groundDofs = [1,2,2*(N+1)+1,2*(N+1)+2]; % Degrees of Freedom that will be eliminated
|
|||
%% Add bars along with geometry
|
|||
% bars are #of Bars X [Area,E,length of bar,angle of bar w.r.t. horizontal
|
|||
% axis, node1 (from node), node2 (to node)]
|
|||
bars = zeros(2*N+2*(N-1),6);
|
|||
for n = 1:N
|
|||
% upper bars
|
|||
lelem = L/N;
|
|||
bars(n,:) = [Aval,Eval,lelem,0,n,n+1];
|
|||
% diagonal bars
|
|||
lelem = sqrt((L/N)^2+H^2);
|
|||
bars(N+n,:) = [Aval,Eval,lelem,pi/4,N+1+n,n+1];
|
|||
end
|
|||
for n = 1:N-1
|
|||
% lower bars
|
|||
lelem = L/N;
|
|||
bars(2*N+n,:) = [Aval,Eval,lelem,0,N+1+n,N+1+n+1];
|
|||
% vertical bars
|
|||
lelem = H;
|
|||
bars(2*N+N-1+n,:) = [Aval,Eval,lelem,pi/2,N+1+n+1,n+1];
|
|||
end
|
|||
%% Assemble all bars into global stiffness and mass matrices
|
|||
% Use a dense matrix for insertion efficiency
|
|||
K = zeros(numDofs);
|
|||
M = zeros(numDofs);
|
|||
% Mass element matrix
|
|||
% the mass matrix is a lumped matrix with half the mass going to one
|
|||
% node and half to the other node
|
|||
unitMassMatrix = diag([1/2 1/2 1/2 1/2]);
|
|||
melem = rho*lelem*unitMassMatrix;
|
|||
for ii=1:size(bars,1)
|
|||
% extract parameters for stiffness and mass matrices
|
|||
Aelem = bars(ii,1);
|
|||
Eelem = bars(ii,2);
|
|||
lelem = bars(ii,3);
|
|||
telem = bars(ii,4);
|
|||
kelem = localStiffness(Aelem,Eelem,lelem,telem); % stiffness matrix
|
|||
n1 = bars(ii,5);
|
|||
n2 = bars(ii,6);
|
|||
% convert to reduced dofs
|
|||
ix = [n1*2-1,2*n1,n2*2-1,n2*2];
|
|||
% element "stamping"
|
|||
K(ix,ix) = K(ix,ix) + kelem; % Add in element's contribution to K
|
|||
M(ix,ix) = M(ix,ix) + melem; % Add in element's contribution to M
|
|||
end
|
|||
% Zero out rows and columns of ground DOFs in stiffness and mass matrix
|
|||
K(groundDofs,:) = 0;
|
|||
K(:,groundDofs) = 0;
|
|||
M(groundDofs,:) = 0;
|
|||
M(:,groundDofs) = 0;
|
|||
% Add ones to diagonal of ground DOFs to freeze them
|
|||
groundDofDiag = sub2ind(size(M),groundDofs,groundDofs);
|
|||
K(groundDofDiag) = 1;
|
|||
M(groundDofDiag) = 1;
|
|||
% convert stiffness and mass matrices to sparse
|
|||
Kr = sparse(K);
|
|||
Mr = sparse(M);
|
|||
% Force vector iz zeros with magnitude on end Y Dof in Y
|
|||
F = zeros(size(Kr,1),1);
|
|||
F(2*(N+1)) = downwardForceMag;
|
|||
%% Solve the system
|
|||
% Transform M*d2x/dt2 + damping*dx/dt + K*x = F into dydt = A*y + f form
|
|||
% Convert 2nd order ODE to first order ODE with y = dx/dt transformation.
|
|||
% The state space is therefore [position,velocity] and in that order
|
|||
[LM,UM] = lu(Mr);
|
|||
A = [-UM\(LM\Kr),-UM\(LM\(dampingCoeffAlpha*Kr+dampingCoeffBeta*Mr));sparse(numDofs,numDofs),speye(numDofs)];
|
|||
f = [UM\(LM\F);zeros(numDofs,1)];
|
|||
myEvaluator = @(t,y) A*y + max(0,(applyTime-t)*f/applyTime)*sin(2*pi*freq*t);
|
|||
initialCondition = zeros(numDofs*2,1); % initial conditions for ode23t
|
|||
% Yr is #Time points X [position, velocity]
|
|||
[~,Yr] = ode23t(myEvaluator,[0 simTime],initialCondition);
|
|||
% Remove the velocity components from Yr
|
|||
Yr = Yr(:,1:numDofs);
|
|||
end
|