Commit 43a63d2d by Jorge Pérez Zerpa

Initial commit.

parents
%% FMTrusS - Force's Method 2D Truss Solver
% The FMTrusS is a simple GNU-Octave script that allows to solve 2D truss
% analysis problems using the Force's Method. The current version
% was generated for educational purposes. NO WARRANTY is provided and
% all errors and/or suggestions are welcome https://www.fing.edu.uy/~jorgepz/ .
addpath( [ pwd '/sources'] );
%% 1- Input data
% Input parameters for geometry and material definition of all the trusses.
% Applied loads and fixed degrees of freedom (supports) are also set here.
FMTruss_input
%% 2- Pre-process
% Computation and storage of length, angle, E and A of all truss elements.
% Static indetermination computation. Equilibrium matrix assembly and
% virtual states definition.
FMTruss_preprocess
%% 3- Process
% Equilibrium equations solved, flexibility matrix assembly and
% virtual forces determined.
FMTruss_process
%% 4- Output
% Plot of truss structure with reactions, normal forces and external loads.
FMTruss_plots
clear all, close all, format compact g
% section properties: vector with the cross-section areas
As = [ 0.001 0.01 ]' ;
% material properties: vector with young moduli values
Es = [ 210e9 ] ;
% --- defintions ---
NodsCoord = [ ...
-10 0 ;
0 10 ;
0 0 ;
10 0 ] ;
% inode jnode material section
ElemConec = [ ...
1 2 1 1 ;
2 3 1 2 ; ...
2 4 1 1 ...
] ;
fixeddofs = [ 1 2 5 6 7 8 ] ;
NodalLoads = [ 2 10000 0 ];
# Deformed scale factor
scalefactor = 1e3;
# row vector with the dofs related to supports which are replaced by virtual forces
virtualforcessupports = [ 7 ] ;
# row vector with the truss elements which are replaced by virtual forces
virtualforceselements = [ ] ;
# degree of freedom whose displacement must be determined (leave empty if none)
unkndispdof = 3 ;
% -----------------------------
LW = 2;
MS = 5;
loadscfactor = max( max(NodsCoord(:,1)) - min(NodsCoord(:,1)) , ...
max(NodsCoord(:,2)) - min(NodsCoord(:,2)) ) ...
/ max( abs(Fext) ) * 0.2 ;
figure, hold on
quiver( NodsCoord(:,1), NodsCoord(:,2) , Fext(1:2:end)*loadscfactor ...
, Fext(2:2:end)*loadscfactor ,0,'g',"filled")
quiver( NodsCoord(:,1), NodsCoord(:,2) , Rext(1:2:end)*loadscfactor ...
, Rext(2:2:end)*loadscfactor ,0,'r',"filled")
for i=1:nelems
xselem = NodsCoord( ElemConec(i,1:2) , 1 ) ;
yselem = NodsCoord( ElemConec(i,1:2) , 2 ) ;
plot( xselem, yselem, 'k--o', 'linewidth', LW, 'markersize', MS );
elemdofs = nodes2dofs( ElemConec(i,1:2), 2 ) ;
if ResultNormalForces(i) >0, colornormalforce='b';
elseif ResultNormalForces(i) <0, colornormalforce='r';
else colornormalforce='k'; end
text( sum(xselem)*0.5, sum(yselem)*0.5, sprintf( '%8.2e', ResultNormalForces(i) ) ...
,'color',colornormalforce, 'fontsize', 14);
end
axis equal, xlabel('x'), ylabel('y')
title('Forces results: green: external loads, red: support reactions and blue normal forces.')
% --- compute lengths and inclination of undeformed elements ---
Lengths = sqrt ( sum( ( NodsCoord( ElemConec(:,2),:) - ...
NodsCoord( ElemConec(:,1),:) ).^2 , 2 ) ) ;
Angles = atan2( ( NodsCoord( ElemConec(:,2),2) - NodsCoord( ElemConec(:,1),2) ) , ...
( NodsCoord( ElemConec(:,2),1) - NodsCoord( ElemConec(:,1),1) ) ) ;
Areas = As( ElemConec(:,4) ) ;
Youngs = Es( ElemConec(:,3) ) ;
nnodes = size( NodsCoord,1);
nelems = size( ElemConec,1);
nfixeddofs = length(fixeddofs) ;
hiperdegree = ( nfixeddofs + nelems ) - 1.0*(2*nnodes) ;
freedofs = 1:(2*nnodes);
freedofs(fixeddofs)=[];
# calculates the total number of force unknowns (element stress and reactions)
nforceunknowns = nfixeddofs + nelems ;
# assembles the vector of virtual force states
if length(virtualforcessupports)>0
virtualforces = [ find( fixeddofs == virtualforcessupports) ...
virtualforceselements+nfixeddofs ] ;
else
virtualforces = [ virtualforceselements+nfixeddofs ] ;
end
# row vector with the indexes of the dofs of the supports left in the isostatic structure
isostaticsupports = 1:nfixeddofs ;
if length(virtualforcessupports)>0
# deletes the entries corresponding to the supports which are replaced by virtual forces
isostaticsupports( find( fixeddofs == virtualforcessupports) ) = [] ;
end
isostaticforceselem = (1:nelems) ;
isostaticforceselem(virtualforceselements) = [] ;
isostaticforceselem = isostaticforceselem + nfixeddofs ;
# row vector with the dofs and elements present in the isostatic structure
isostaticforces = [ isostaticsupports isostaticforceselem ] ;
# Equilibrium matrix assembly
# each row contains the coefficients of each equilibrium equation for the nodes
# each column contains the coefficients corresponding to an unknown reaction
# or normal force
Meq = zeros( 2*nnodes , nforceunknowns ) ;
for i=1:nfixeddofs
Meq(fixeddofs(i),i) = 1.0 ;
end
# stress
for i=1:nelems
l = Lengths(i) ;
ang = Angles(i);
nodi = ElemConec(i,1);
nodj = ElemConec(i,2);
elemdofs = nodes2dofs ( [ nodi nodj ]' ,2) ;
ca = cos(ang); sa = sin(ang);
auxproj = [ -ca; -sa; ca; sa] ;
Meq( elemdofs , i+nfixeddofs ) = - auxproj ;
end
# assemble the external nodal loads vector
Fext = zeros(2*nnodes,1);
for i=1:size(NodalLoads,1)
aux = nodes2dofs ( NodalLoads(i,1), 2 ) ;
Fext( aux ) = Fext( aux ) + NodalLoads(i,2:3)' ;
end
# force states independent terms
ForceIndepTerm = zeros( 2*nnodes,1+hiperdegree);
ForceIndepTerm(:,1) = - Fext ;
Meqred = Meq;
for j=1:hiperdegree
ForceIndepTerm(:,1+j) = -Meq(:,virtualforces(j)) ;
Meqred(:,virtualforces(j) ) = [] ;
end
x = Meqred \ ForceIndepTerm ;
supportreactions = zeros(nfixeddofs, hiperdegree+1) ;
supportreactions( isostaticsupports , : ) = x( 1:length(isostaticsupports) ,:) ;
if length(virtualforcessupports)>0
aux = find( fixeddofs == virtualforcessupports) ;
for j=1:hiperdegree
supportreactions( aux(j),1+j ) = 1 ;
end
end
fprintf('The support reactions in the real and virtual states are:\n');
supportreactions
Ns = zeros(nelems,hiperdegree+1) ;
%~ Ns( isostaticforceselem-length(isostaticsupports) , : ) = x( isostaticforceselem, : ) ;
Ns( isostaticforceselem- nfixeddofs , : ) = x( (length(isostaticsupports)+1):end, : ) ;
if length( virtualforceselements ) > 0
for i=1:length(virtualforceselements)
Ns(virtualforceselements,i+1) = 1 ;
end
end
fprintf('The normal forces of the bars in the real and virtual states are:\n');
Ns
Kf = zeros( hiperdegree, hiperdegree) ;
Ff = zeros( hiperdegree, 1) ;
for i=1:hiperdegree
for j = 1:hiperdegree
Kf(i,j) = sum( Ns(:,1+i).*Ns(:,1+j) ./ ( Youngs .* Areas ) .* Lengths ) ;
end
%
Ff(i) = - sum( Ns(:,1) .* Ns(:,1+i) ./ ( Youngs .* Areas ) .* Lengths ) ;
end
X = Kf \ Ff
ResultReactions = supportreactions * [ 1; X]
ResultNormalForces = Ns * [ 1 ; X ]
Rext = zeros( 2*nnodes,1) ;
Rext(fixeddofs) = ResultReactions ;
Fextaux = zeros(2*nnodes,1);
Fextaux(unkndispdof) = 1 ;
xaux = Meqred \ - Fextaux ;
Nsaux = zeros( nelems ,1) ;
Nsaux( isostaticforceselem- nfixeddofs ) = xaux ( (length(isostaticsupports)+1):end )
disp = sum( ResultNormalForces .* Nsaux ./ ( Youngs .* Areas ) .* Lengths )
## Auxiliar function - nodes2dofs
# Auxiliar function to perform conversion from node number to correponding degree of freedom
##
function v= nodes2dofs (u,degree)
v=zeros(degree*length(u),1);
for i=1:length(u)
v( [ (degree*(i-1) + 1):(degree*i) ] ) = [ (degree*(u(i)-1) + 1):(degree*u(i)) ] ;
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment