How to compute the eigenvalue problem of a truss

Question

Rita Akonobi 2020-10-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/613336-how-to-compute-the-eigenvalue-problem-of-a-truss

编辑： Rita Akonobi 2020-10-14

Hi guys, I'm trying to implement the natural frequency constraint in y code but i'm not sure, the function 'eig' solves the eigenvalue problem of a truss as I am getting negative values. Also, when I try to extract the minimum value of the frequency, I notice i get a row vector and not a scalar. Please can someone help. I attached the codes below. The eigenvalue problem is det[K -ω^2M] = 0

ELEMENTS = 11
1   1   2  0.00361  1.79e-6   
2   2   3  0.00361  1.79e-6
3   3   4  0.00361  1.79e-6
4   5   6  0.00361  1.79e-6
5   6   7  0.00361  1.79e-6
6   2   5  0.00203  0.90e-6
7   3   6  0.00203  0.90e-6
8   4   7  0.00203  0.90e-6
9   1   5  0.00361  1.79e-6
10  3   5  0.00361  1.79e-6
11  4   6  0.00203  1.79e-6
NODE_COORDINATES = 7
1    0.0    0.0
2    4.0    0.0
3    8.0    0.0
4   12.0    0.0
5   4.0     3.0
6    8.0    3.0
7   12.0    3.0
NODES_WITH_PRESCRIBED_DISPLACEMENTS = 2
1  11 0.000 0.000
7  10 0.000 0.000
YOUNG_MODULUS =100000000.0
YIELD STRESS =100000.0
DENSITY =8050
NATURAL FREQUENCY =2.5
NODES_WITH_POINT_LOAD = 3
2 0 -40
3 0 -40
4 0 -20
PLOTTING_AMPLIFICATION_FACTOR =30
NODES_WITH_VARIABLE_COORDINATES = 3
5   11 0.5  1.0  6.0  6.0
6   11 4.0  1.0  10.0 6.0
7   01 5.0  1.0  16.0 6.0
DEFLECTION_LIMIT = 6
2   11  0.01 0.03
3   11  0.01 0.03
4   11  0.01 0.03
5   11  0.01 0.03
6   11  0.01 0.03
7   11  0.01 0.03
%**********************************************************************
% Reads all data for a 2D linear elastic truss
%
% HISTORY
% October   2018, :   Initial coding (from "truss2D.m")
%***********************************************************************
function [coorddofvar,fixeddata] = truss2D_readdata(filename)
%
% Set some basic variables
% ========================
% Number of degrees of freedom per node
ndofn=2;
% Number of nodes of the element
nnode=2;
%
% Read information from data file
% ===============================
%
fid = fopen(filename, 'r');
title = fscanf(fid, 'TITLE = %s',1);
%
% Total number of elements in the mesh
nelem = fscanf(fid, '\nELEMENTS = %d', 1);
%
%Read  Table of connectivities
%
lnods = fscanf(fid, '\n%d %d %d %f %f', [nnode+3,nelem]);
lnods = lnods';
sortrows(lnods,1);
%...cross-sectional areas
csarea = lnods(:,4);
%...second moment of areas
sndmoa = lnods(:,5);
%...store only connectivities in lnods
lnods = lnods(:,2:nnode+1);
%...create table of element global degrees of freedom
eldofX = lnods*ndofn-1;
eldofY = lnods*ndofn;
eldofs = [ eldofX(:,1) eldofY(:,1) eldofX(:,2) eldofY(:,2) ];
%
% Read Nodal coordinates
%
npoin = fscanf(fid, '\nNODE_COORDINATES = %d', 1);
coord = fscanf(fid, '\n%d %f %f', [3, npoin]);
coord = coord';sortrows(coord,1);nodnumbers=coord(:,1);
coord = coord(:,2:3);
%...create table of nodal degrees of freedom
nodofs = [ nodnumbers*ndofn-1 nodnumbers*ndofn ];
%
%
% Read Prescribed displacements
%
nnodefix = fscanf(fid,'\nNODES_WITH_PRESCRIBED_DISPLACEMENTS = %d',1);
ipresc   = fscanf(fid, '\n%d %d %f %f', [2+ndofn, nnodefix]);
ipresc   = ipresc';fixednodes=ipresc(:,1);
%...create tables of fixed dofs and corresponding prescribed values
ifdoffix = zeros(npoin*ndofn,1);
icount = 0;
for inodefix = 1:nnodefix
    ipoin=ipresc(inodefix,1);
    dofX=nodofs(ipoin,1);
    dofY=nodofs(ipoin,2);
    if  ipresc(inodefix,2)==11
        ifdoffix(dofX)=1;
        ifdoffix(dofY)=1 ;
        icount = icount+1;
        valuedoffix(icount)   = ipresc(inodefix,3);
        fixeddoftable(icount) = dofX;
        icount = icount+1;
        valuedoffix(icount)   = ipresc(inodefix,4);
        fixeddoftable(icount) = dofY;
    elseif  ipresc(inodefix,2)==1
        ifdoffix(dofY)=1;
        icount = icount+1;
        valuedoffix(icount)   = ipresc(inodefix,4);
        fixeddoftable(icount) = dofY;
    elseif  ipresc(inodefix,2)==10
        ifdoffix(dofX)=1;
        icount = icount+1;
        valuedoffix(icount)   = ipresc(inodefix,3);
        fixeddoftable(icount) = dofX;
    elseif  ipresc(inodefix,2)==0
    else
        error('Wrong displacement prescription code in data file')
    end
end
%...create table of free dofs by subtracting the set of fixed dofs
%   from the set of all dofs
ngdof=npoin*ndofn;
alldoftable=[1:ngdof]';
freedoftable = setxor(alldoftable,fixeddoftable);
%
% Read Material properties
%
matprop.young = fscanf(fid, '\nYOUNG_MODULUS = %f', 1);
matprop.yield = fscanf(fid, '\nYIELD STRESS = %f', 1);
%added by Ritamatprop
matprop.density = fscanf(fid, '\nDENSITY = %f', 1);
matprop.frequency = fscanf(fid,'\nNATURAL FREQUENCY =%f', 1);
%
%Read  Load vector (point loads only)
%
npload = fscanf(fid,'\nNODES_WITH_POINT_LOAD = %d',1);
pload  = fscanf(fid, '\n%d %f %f', [1+ndofn, npload]);
pload  = pload';
%...add point loads to global load vector
F=zeros(ngdof,1);
loadednodes=pload(:,1);
loadeddofs=nodofs(loadednodes,:);
F(loadeddofs) = pload(:,2:3);
%
% Reads deformation amplification factor (for plotting only)
amplfact = fscanf(fid,'\nPLOTTING_AMPLIFICATION_FACTOR = %d',1);
%
% Reads upper and lower bounds of variable coordinates to be optimised
%
nnodevar = fscanf(fid, '\nNODES_WITH_VARIABLE_COORDINATES = %d', 1);
ivar   = fscanf(fid, '\n%d %d %f %f %f %f', [2+2*ndofn, nnodevar]);
ivar   = ivar';
varnodes=ivar(:,1);%...create tables of variable coordinates and corresponding prescribed
%   upper and lower bounds
ifcoordvar = zeros(npoin*ndofn,1); icount = 0;
upperbound=[];
lowerbound=[];
dofcoordvar=[];
for inodevar = 1:nnodevar
    ipoin=ivar(inodevar,1);
    dofX=nodofs(ipoin,1);
    dofY=nodofs(ipoin,2);
    if  ivar(inodevar,2)==11
        % both x- and y-coordinates are variable
        ifcoordvar(dofX)=1;
        ifcoordvar(dofY)=1;
        icount = icount+1;
        lowerbound(icount) = ivar(inodevar,3);
        upperbound(icount) = ivar(inodevar,5);
        dofcoordvar(icount) = dofX;
        icount = icount+1;
        lowerbound(icount) = ivar(inodevar,4);
        upperbound(icount) = ivar(inodevar,6);
        dofcoordvar(icount) = dofY;
    elseif  ivar(inodevar,2)==1
        % only y-coordinate is variable
        ifcoordvar(dofY)=1;
        icount = icount+1;
        lowerbound(icount) = ivar(inodevar,4);
        upperbound(icount) = ivar(inodevar,6);
        dofcoordvar(icount) = dofY;
    elseif  ivar(inodevar,2)==10
        % only x-coordinate is variable
        ifcoordvar(dofX)=1;
        icount = icount+1;
        lowerbound(icount) = ivar(inodevar,3);
        upperbound(icount) = ivar(inodevar,5);
        dofcoordvar(icount) = dofX;
    elseif  ivar(inodevar,2)==0
        disp('WARNING: Node listed with no variable coordinates for optimisation process')
    else
        error('Wrong variable coordinate prescription code in data file')
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RITA CODE STARTS FROM HERE
%Read deflection limit from data file
nnodelimit = fscanf(fid,'\nDEFLECTION_LIMIT = %d',1);
ilimit   = fscanf(fid, '\n%d %d %f %f', [2+ndofn, nnodelimit]); 
ilimit   = ilimit';limitednodes=ilimit(:,1);   
%...create tables of fixed dofs and corresponding prescribed values
ifdoflimit = zeros(npoin*ndofn,1); iicount = 0;
for inodelimit = 1:nnodelimit    
    ipoin=ilimit(inodelimit,1);     
    dofX=nodofs(ipoin,1);
    dofY=nodofs(ipoin,2);     
    if  ilimit(inodelimit,2)==11        
        ifdoflimit(dofX)=1;   
        ifdoflimit(dofY)=1;        
        iicount = iicount+1;         
        valuedoflimit(iicount)= ilimit(inodelimit,3);         
        limiteddoftable(iicount) = dofX;         
        iicount = iicount+1;         
        valuedoflimit(iicount) = ilimit(inodelimit,4);        
        limiteddoftable(iicount)= dofY;    
    elseif  ilimit(inodelimit,2)==1        
        ifdoflimit(dofY)=1;         
        iicount = iicount+1;         
        valuedoflimit(iicount)= ilimit(inodelimit,4);         
        limiteddoftable(iicount) = dofY;    
    elseif  ilimit(inodefix,2)==10         
        ifdoflimit(dofX)=1;        
        iicount = iicount+1;         
        valuedoflimit(iicount) = ilimit(inodelimit,3);         
        limiteddoftable(iicount) = dofX;     
    elseif  ilimit(inodelimit,2)==0     
    else         error('Deflection limitaion value is not inserted')    
    end
end
%************************************************************************
%
%                           IMPORTANT
%                 (DEFINIFION OF DESIGN VARIABLES)
%
% set of dofs associated with variable coordinates
dofcoordfix=setxor(alldoftable,dofcoordvar);
if isempty(dofcoordvar)
    disp('No variables in optimisation problem - all nodes have fixed coordinates')
    error('No variables in optimisation problem - all nodes have fixed coordinates')
end
%
% arrange coordinates following dof ordering (in one-dimensional array, as
% in global displacement vector)
coorddof=zeros(ngdof,1);
for ipoin=1:npoin
    dofx=ipoin*2-1;
    dofy=ipoin*2;
    coorddof(dofx)=coord(ipoin,1);
    coorddof(dofy)=coord(ipoin,2);
end
% split coordinates array into variable and fixed coordinates
coorddofvar=coorddof(dofcoordvar);
coorddoffix=coorddof(dofcoordfix);
%AppRita
coorddoflimit=coorddof(alldoftable); 
%x = coorddofvar
%************************************************************************
%
% Store fixed data in structure
fixeddata.amplfact=amplfact;
fixeddata.coorddoffix=coorddoffix;
fixeddata.csarea=csarea;
fixeddata.sndmoa=sndmoa;
fixeddata.dofcoordfix=dofcoordfix;
fixeddata.dofcoordvar=dofcoordvar;
fixeddata.eldofs=eldofs;
fixeddata.F=F;
fixeddata.fixeddoftable=fixeddoftable;
fixeddata.freedoftable=freedoftable;
fixeddata.lnods=lnods;
fixeddata.lowerbound=lowerbound;
fixeddata.matprop=matprop;
fixeddata.ndofn=ndofn;
fixeddata.nelem=nelem;
fixeddata.ngdof=ngdof;
fixeddata.npoin=npoin;
fixeddata.upperbound=upperbound;
fixeddata.valuedoffix=valuedoffix;
%fixeddata.coorddoflimit=coorddoflimit; %appended by Rita
%fixeddata.coorddoflimit=coorddoflimit; 
fixeddata.ifdoflimit=ifdoflimit;
fixeddata.limiteddoftable=limiteddoftable;
fixeddata.valuedoflimit=valuedoflimit; 
fixeddata.csarea=csarea; 
fixeddata.sndmoa=sndmoa; 
% Close file(s)
status = fclose(fid);
%
end

%**********************************************************************
% Constraint function for optimization of linear elastic 2D trusses
%
% HISTORY
% October 2018, EA de Souza Neto:   Initial coding
%***********************************************************************
function [C,Ceq] = truss2D_nlcon(coorddofvar)
%
global fixeddata
global R
global U
%
% Retrieve fixed data from structure
coorddoffix=fixeddata.coorddoffix;
csarea=fixeddata.csarea;
dofcoordfix=fixeddata.dofcoordfix;
dofcoordvar=fixeddata.dofcoordvar;
eldofs=fixeddata.eldofs;
F=fixeddata.F;
fixeddoftable=fixeddata.fixeddoftable;
freedoftable=fixeddata.freedoftable;
lnods=fixeddata.lnods;
matprop=fixeddata.matprop;
ndofn=fixeddata.ndofn;
nelem=fixeddata.nelem;
ngdof=fixeddata.ngdof;
npoin=fixeddata.npoin;
valuedoffix=fixeddata.valuedoffix;
sndmoa=fixeddata.sndmoa;
%fixeddata.coorddoflimit=coorddoflimit; %appended by Rita
ifdoflimit=fixeddata.ifdoflimit;
limiteddoftable=fixeddata.limiteddoftable; 
valuedoflimit=fixeddata.valuedoflimit; 
%coorddofvar=fixeddata.coorddofvar; 
%
% arrange coordinates in usual 2-dimensional array format (node by node)
coorddof(dofcoordfix)=coorddoffix;
coorddof(dofcoordvar)=coorddofvar;
%...compute elements length and angle sine/cosine
elength=zeros(nelem,1);
for ielem=1:nelem
    elvec=coorddof(eldofs(ielem,3:4))-coorddof(eldofs(ielem,1:2));
    elength(ielem)=sqrt(elvec*elvec');
    
end
% arrange coordinates in usual 2-dimensional array format (node by node)
coorddof(dofcoordfix)=coorddoffix;
coorddof(dofcoordvar)=coorddofvar;
% coord=zeros(npoin,2);
% for ipoin=1:npoin
%     dofx=ipoin*2-1;dofy=ipoin*2;
%     coord(ipoin,1)=coorddof(dofx);coord(ipoin,2)=coorddof(dofy);
% end
%...compute elements length and angle sine/cosine
elength=zeros(nelem,1);
ecos=zeros(nelem,1);
esin=zeros(nelem,1);
for ielem=1:nelem
    elvec=coorddof(eldofs(ielem,3:4))-coorddof(eldofs(ielem,1:2));
    elength(ielem)=sqrt(elvec*elvec');
    ecos(ielem)=elvec(1)/elength(ielem);
    esin(ielem)=elvec(2)/elength(ielem);
end
%
% Solve FE truss equilibrium problem
% ----------------------------------
%
% Compute global stiffness and assemble system of equations
%
K = zeros(npoin*ndofn,npoin*ndofn);
for ielem = 1:nelem
%...compute element stiffness
    ke = stiffTRUSS2D(csarea(ielem),elength(ielem),ecos(ielem),esin(ielem),matprop.young);  
%...add element contribution to global stiffness
    gpos = eldofs(ielem,:);
    K(gpos,gpos) = K(gpos,gpos) + ke;
end
% Compute global mass matrix by Rita
M = zeros(npoin*ndofn,npoin*ndofn);
for ielem =1:nelem
    %compute individual element mass matrix
   me = element_mass_matrix(matprop.density,csarea(ielem),elength(ielem));
%...add element contribution to global mass matrix
    gpos = eldofs(ielem,:);
   M(gpos,gpos) = M(gpos,gpos) + me;
end
K;
%...assemble reduced stifness matrix and load vector
%   (remove prescibed d.o.f's)
Kstar = K(freedoftable,freedoftable);
Fstar = F(freedoftable);
%...add contributions from prescribed displacements to
%   right hand side (Fstar)
Fstar=Fstar-K(freedoftable,fixeddoftable)*valuedoffix';
%
% Solve system of equations for unknown displacement
% vector Ustar
%
Ustar = Kstar \ Fstar;
%...assemble full global displacement vector
U=zeros(ngdof,1);
U(freedoftable)=Ustar;
U(fixeddoftable)=valuedoffix;
%
% Compute nodal reactions
%
R=zeros(ngdof,1);
R(fixeddoftable)=K(fixeddoftable,:)*U-F(fixeddoftable);
%
%
% Post-processing of results
% --------------------------
%
% Compute element stresses
%
young=matprop.young;
beall=[ecos esin -ecos -esin]./elength;
Bg=zeros(nelem,ngdof);
for ielem = 1:nelem
%...compute element B-matrix
    be = beall(ielem,:);
%...compute global B-matrix
    Bg(ielem,eldofs(ielem,:))=be;
end
% element stresses
estress=young*Bg*U;
%creating an array of tension=0 and compression=estress
%estress
% element von Mises stresses
vmestress=abs(estress);
%Compute natural frequency(Rita)
[V,D] =eig(K,M);
eig_vec = V;
nfreq = D;
%Extract minimum frequency (Rita)
min_freq=abs(min(nfreq));
% Compute stress constraint array
% ------------------------------
yield=matprop.yield;
%Compute frequency constraint (Rita)
frequency=matprop.frequency;
C=(vmestress-yield)/yield; % array of inequality constraints
Ceq=[]; % no equality constraints in present case
%Calculate buckling stress of each element
bstress=zeros (nelem,1);
bstress=((pi^2)*young *sndmoa)./(csarea.*(elength.*elength));
%Extract the compressive stress in the estress
compstress=zeros(length(estress),1);
compstress(estress<0)=estress(estress<0);
compstress=abs(compstress);
%compare compstress to bstress. If compstress in greater than bstress, the
%element will buckle
Cbuck=(compstress-bstress)./bstress;
C=[C;Cbuck];
%Calculation deflection limit
dof=zeros(ngdof,1); icount=0;
for idof=1:ngdof 
    if ifdoflimit(idof)==1
        icount=icount+1;
        dof=limiteddoftable(icount);
        maxdefl=valuedoflimit(icount);
        elemdefl=(U(dof));
        elemdefl=abs(elemdefl);
        Cdef=(elemdefl-maxdefl)/maxdefl;
        C=[C;Cdef];
    end
end
%Calculation of Natural Frequency
C_freq=(min_freq -frequency)/frequency;
C_freq=C_freq';
C=[C;C_freq];
end