Apostila MEF - Método dos Elementos Finitos (Rubens - PUC)

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

calfem_beam2crd.m
 function [excd,eycd]=beam2crd(ex,ey,ed,mag)
%-------------------------------------------------------------
% PURPOSE
% Calculate the element continous displacements for a 
% number of identical 2D Bernoulli beam elements. 
% 
% INPUT: ex,ey,
% ed,
% mag 
%
% OUTPUT: excd,eycd 
%-------------------------------------------------------------
% LAST MODIFIED: P-E AUSTRELL 1993-10-15
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
% 
 [nie,ned]=size(ed);
 
 for i=1:nie
 
 b=[ex(i,2)-ex(i,1) ey(i,2)-ey(i,1)]; L=sqrt(b*b'); n=b/L;
 
 G=[n(1) n(2) 0 0 0 0;
 -n(2) n(1) 0 0 0 0;
 0 0 1 0 0 0;
 0 0 0 n(1) n(2) 0;
 0 0 0 -n(2) n(1) 0; 
 0 0 0 0 0 1];
 
 d=ed(i,:)'; dl=G*d ;
 %
 xl=[0:L/20:L]'; one=ones(size(xl));
 %
 Cis=[-1 1;
 L 0]/L; ds=[dl(1);dl(4)];
 % 
 ul=([xl one]*Cis*ds)'; 
 % 
 Cib=[ 12 6*L -12 6*L;
 -6*L -4*L^2 6*L -2*L^2;
 0 L^3 0 0;
 L^3 0 0 0]/L^3; 
 % 
 db=[dl(2);dl(3);dl(5);dl(6)];
 % 
 vl=([xl.^3/6 xl.^2/2 xl one]*Cib*db)';
 % 
 cld=[ ul ;
 vl ]; A=[n(1) -n(2);
 n(2) n(1)]; cd=A*cld;
 %
 xyc=A(:,1)*xl'+ [ex(i,1);ey(i,1)]*one';
 % 
 excd(i,:)=xyc(1,:)+mag*cd(1,:); 
 eycd(i,:)=xyc(2,:)+mag*cd(2,:); 
 end
%--------------------------end--------------------------------
calfem_beam2d.m
function [Ke,Me,Ce]=beam2d(ex,ey,ep)
% [Ke,Me]=beam2d(ex,ey,ep)
% [Ke,Me,Ce]=beam2d(ex,ey,ep)
%-------------------------------------------------------------
% PURPOSE
% Calculate the stiffness matrix Ke, the mass matrix Me 
% and the damping matrix Ce for a 2D elastic Bernoulli
% beam element. 
%
% INPUT: ex = [x1 x2]
% ey = [y1 y2] element node coordinates 
%
% ep = [E A I m (a b)] 
% E: Young's modulus
% A: cross section area
% I: moment of inertia
% m: mass per unit length
% a,b: damping coefficients,
% Ce=aMe+bKe 
% 
% OUTPUT: Ke : element stiffness matrix (6 x 6)
% Me : element mass matrix 
% Ce : element damping matrix, optional
%-------------------------------------------------------------
 
% LAST MODIFIED: K Persson 1995-08-23
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
 b=[ ex(2)-ex(1); ey(2)-ey(1) ];
 L=sqrt(b'*b); n=b/L;
 
 E=ep(1); A=ep(2); I=ep(3); m=ep(4);
 a=0 ; b=0 ;
 if length(ep)==6 ; a=ep(5) ; b=ep(6) ; end
%
 Kle=[E*A/L 0 0 -E*A/L 0 0 ;
 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2;
 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;
 -E*A/L 0 0 E*A/L 0 0 ;
 0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;
 0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L];
%
 Mle=m*L/420*[140 0 0 70 0 0 ;
 0 156 22*L 0 54 -13*L ;
 0 22*L 4*L^2 0 13*L -3*L^2 ;
 70 0 0 140 0 0 ;
 0 54 13*L 0 156 -22*L ;
 0 -13*L -3*L^2 0 -22*L 4*L^2];
%
 Cle=a*Mle+b*Kle;
%
 G=[n(1) n(2) 0 0 0 0;
 -n(2) n(1) 0 0 0 0;
 0 0 1 0 0 0;
 0 0 0 n(1) n(2) 0;
 0 0 0 -n(2) n(1) 0;
 0 0 0 0 0 1];
%
 Ke=G'*Kle*G; Me=G'*Mle*G; Ce=G'*Cle*G;
%--------------------------end--------------------------------
 
calfem_beam2d_t.m
function [Ke,Me,Ce]=beam2d(ex,ey,ep)
% [Ke,Me]=beam2d(ex,ey,ep)
% [Ke,Me,Ce]=beam2d(ex,ey,ep)
%-------------------------------------------------------------
% PURPOSE
% Calculate the stiffness matrix Ke, the mass matrix Me 
% and the damping matrix Ce for a 2D elastic Bernoulli
% beam element. 
%
% INPUT: ex = [x1 x2]
% ey = [y1 y2] element node coordinates 
%
% ep = [EI mm (a b)] 
% EI: Rigidez da viga
% mm: mass per unit length
% a,b: damping coefficients,
% Ce=aMe+bKe 
% 
% OUTPUT: Ke : element stiffness matrix (6 x 6)
% Me : element mass matrix 
% Ce : element damping matrix, optional
%-------------------------------------------------------------
 
% LAST MODIFIED: K Persson 1995-08-23
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
%% -------------------------------------------------------------------
%% Ultima atualizacao: THIAGO RITTO - JAN/05
%%--------------------------------------------------------------------
 b=[ ex(2)-ex(1); ey(2)-ey(1) ];
 L=sqrt(b'*b); n=b/L;
 
 EI=ep(1); mm=ep(2);
 a=0 ; b=0 ;
 if length(ep)==4 ; a=ep(3) ; b=ep(4) ; end
%
 Kle=[ 12*EI/L^3 6*EI/L^2 -12*EI/L^3 6*EI/L^2;
 6*EI/L^2 4*EI/L -6*EI/L^2 2*EI/L;
 -12*EI/L^3 -6*EI/L^2 12*EI/L^3 -6*EI/L^2;
 6*EI/L^2 2*EI/L -6*EI/L^2 4*EI/L];
%
 Mle=mm*L/420*[ 156 22*L 54 -13*L ;
 22*L 4*L^2 13*L -3*L^2 ;
 54 13*L 156 -22*L ;
 -13*L -3*L^2 -22*L 4*L^2];
%
 Cle=a*Mle+b*Kle;
% soh aceita elementosna na horizontal - nao faz transformacao de
% coordenadas
 G=[1 0 0 0;
 0 1 0 0;
 0 0 1 0;
 0 0 0 1];
%
 Ke=G'*Kle*G;
 Me=G'*Mle*G;
 Ce=G'*Cle*G;
%--------------------------end--------------------------------
 
calfem_beam2e.m
function [Ke,fe]=beam2e(ex,ey,ep,eq);
% Ke=beam2e(ex,ey,ep)
% [Ke,fe]=beam2e(ex,ey,ep,eq)
%---------------------------------------------------------------------
% PURPOSE
% Compute the stiffness matrix for a two dimensional beam element. 
% 
% INPUT: ex = [x1 x2]
% ey = [y1 y2] element node coordinates
%
% ep = [E A I] element properties
% E: Young's modulus
% A: Cross section area
% I: Moment of inertia
%
% eq = [qx qy] distributed loads, local directions
% 
% OUTPUT: Ke : element stiffness matrix (6 x 6)
%
% fe : element load vector (6 x 1)
%--------------------------------------------------------------------
% LAST MODIFIED: K Persson 1995-08-23
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
 b=[ ex(2)-ex(1); ey(2)-ey(1) ];
 L=sqrt(b'*b); n=b/L;
 E=ep(1); A=ep(2); I=ep(3);
 
 qx=0; qy=0; if nargin>3; qx=eq(1); qy=eq(2); end
 Kle=[E*A/L 0 0 -E*A/L 0 0 ;
 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2;
 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;
 -E*A/L 0
0 E*A/L 0 0 ;
 0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;
 0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L];
 
 fle=L*[qx/2 qy/2 qy*L/12 qx/2 qy/2 -qy*L/12]';
 G=[n(1) n(2) 0 0 0 0;
 -n(2) n(1) 0 0 0 0;
 0 0 1 0 0 0;
 0 0 0 n(1) n(2) 0;
 0 0 0 -n(2) n(1) 0;
 0 0 0 0 0 1];
 Ke=G'*Kle*G; fe=G'*fle; 
%--------------------------end--------------------------------
calfem_coordxtr.m
function [Ex,Ey,Ez]=coordxtr(Edof,Coord,Dof,nen)
%[Ex,Ey,Ez]=coordxtr(Edof,Coord,Dof,nen)
%-------------------------------------------------------------
% PURPOSE
% Extract nodal coordinate data from the global coordinate 
% matrix for a number of elements with equal number of 
% element nodes and dof's. 
%
% INPUT: Edof : topology matrix , dim(t)= nie x ned+1
% nie= number of identical elements
% ned= number of element dof's 
%
% Coord: global coordinate matrix 
%
% Dof: global nodal dof matrix 
%
% nen: number of element nodes
%
% OUTPUT: Ex,Ey,Ez : element coordinate matrices
% Ex=[x1 x2 ...xnen; one row for each element
% ... ... ;
% nel ... ] 
% dim= nel x nen ; nel:number of elemnts 
%-------------------------------------------------------------
% LAST MODIFIED: P-E Austrell 1993-10-14 
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
 [nel,dum]=size(Edof);
 ned=dum-1;
 [n,nsd]=size(Coord);
 [n,nd]=size(Dof);
 nend=ned/nen;
 %
 for i = 1:nel
 nodnum=zeros(1,nen);
 for j = 1:nen
 check=Dof(:,1:nend)-ones(n,1)*Edof(i,(j-1)*nend+2:j*nend+1);
 [indx,dum]=find(check==0);
 nodnum(j)=indx(1);
 end
 % 
 Ex(i,:)=Coord(nodnum,1)';
 if nsd>1
 Ey(i,:)=Coord(nodnum,2)';
 end
 if nsd>2 
 Ez(i,:)=Coord(nodnum,3)';
 end
 end 
%--------------------------end--------------------------------
calfem_eigen.m
 function [L,X]=eigen(K,M,b)
% [L]=eigen(K,M)
% [L]=eigen(K,M,b)
% [L,X]=eigen(K,M)
% [L,X]=eigen(K,M,b)
%-------------------------------------------------------------
% PURPOSE
% Solve the generalized eigenvalue problem
% [K-LM]X = 0, considering boundary conditions.
%
% INPUT:
% K : global stiffness matrix, dim(K)= nd x nd
% M : global mass matrix, dim(M)= nd x nd
% b : boundary condition matrix
% dim(b)= nb x 1
% OUTPUT:
% L : eigenvalues stored in a vector with length (nd-nb) 
% X : eigenvectors dim(X)= nd x nfdof, nfdof : number of dof's
%-------------------------------------------------------------
% LAST MODIFIED: H Carlsson 1993-09-21
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
 [nd,nd]=size(K);
 fdof=[1:nd]';
%
 if nargin==3
 pdof=b(:);
 fdof(pdof)=[]; 
 if nargout==2
 [X1,D]=eig(K(fdof,fdof),M(fdof,fdof));
 [nfdof,nfdof]=size(X1);
 for j=1:nfdof;
 mnorm=sqrt(X1(:,j)'*M(fdof,fdof)*X1(:,j));
%warning off MATLAB:divideByZero
 X1(:,j)=X1(:,j)/mnorm;
 end
 d=diag(D);
 [L,i]=sort(d);
 X2=X1(:,i);
 X=zeros(nd,nfdof);
 X(fdof,:)=X2;
 else
 d=eig(K(fdof,fdof),M(fdof,fdof));
 L=sort(d);
 end
 else
 if nargout==2
 [X1,D]=eig(K,M);
 for j=1:nd;
 mnorm=sqrt(X1(:,j)'*M*X1(:,j));
 X1(:,j)=X1(:,j)/mnorm;
 end
 d=diag(D);
 [L,i]=sort(d);
 X=X1(:,i);
 else
 d=eig(K,M);
 L=sort(d);
 end
 end
%--------------------------end--------------------------------
calfem_eldisp2.m
function [magnfac]=eldisp2(ex,ey,ed,plotpar,magnfac)
%eldisp2(ex,ey,ed,plotpar,magnfac)
%[magnfac]=eldisp2(ex,ey,ed,plotpar)
%[magnfac]=eldisp2(ex,ey,ed)
%-------------------------------------------------------------
% PURPOSE 
% Draw the deformed 2D mesh for a number of elements of 
% the same type. Supported elements are:
% 
% 1) -> bar element 2) -> beam el. 
% 3) -> triangular 3 node el. 4) -> quadrilateral 4 node el. 
% 5) -> 8-node isopar. element
% INPUT
% ex,ey:.......... nen: number of element nodes
% nel: number of elements 
% ed: element displacement matrix
%
% plotpar=[ linetype, linecolor, nodemark] 
%
% linetype=1 -> solid linecolor=1 -> black
% 2 -> dashed 2 -> blue
% 3 -> dotted 3 -> magenta
% 4 -> red
% nodemark=1 -> circle 
% 2 -> star 
% 0 -> no mark 
%
% magnfac: magnification factor for displacements 
% 
% Rem. Default if magnfac and plotpar is left out is auto magnification 
% and dashed white lines with circles at nodes -> plotpar=[2 1 1] 
%-------------------------------------------------------------
% LAST MODIFIED: J Lindemann 1999-01-29
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
%
 if ~((nargin==3)|(nargin==4)|(nargin==5))
 disp('??? Wrong number of input arguments!')
 return
 end
 
 a=size(ex); b=size(ey);
 
 if (a-b)==[0 0]
 nen=a(2); 
 else
 disp('??? Check size of coordinate input arguments!')
 return
 end
 
 c=size(ed);
 
 if ~(c(1)==a(1))
 disp('??? Check size of displacement input arguments!')
 return 
 end
 
 ned=c(2);
 
 dxmax=max(max(ex')-min(ex')); dymax=max(max(ey')-min(ey'));
 dlmax=max(dxmax,dymax);
 edmax=max(max(abs(ed)));
 krel=0.1;
 
 if nargin==3; 
 plotpar=[2 1 1]; magnfac=krel*dlmax/edmax;
 elseif nargin==4
 magnfac=krel*dlmax/edmax;
 end
 
 [s1,s2]=calfem_pltstyle(plotpar);
 k=magnfac;
% ********** Bar or Beam elements *************
 if nen==2 
 if ned==4 % ----------- Bar elements -------------
 x=(ex+k*ed(:,[1 3]))'; 
 y=(ey+k*ed(:,[2 4]))';
 xc=x;
 yc=y;
 elseif ned==6 % -------- Beam elements ------------
 x=(ex+k*ed(:,[1 4]))'; 
 y=(ey+k*ed(:,[2 5]))';
 [exc,eyc]=calfem_beam2crd(ex,ey,ed,k);
 xc=exc';
 yc=eyc';
 end 
% ********* 2D triangular elements ************
 elseif nen==3
 x=(ex+k*ed(:,[1 3 5]))'; 
 y=(ey+k*ed(:,[2 4 6]))';
 xc=[x; x(1,:)];
 yc=[y; y(1,:)];
 
% ********* 2D quadrilateral elements *********
 elseif nen==4 
 x=(ex+k*ed(:,[1 3 5 7]))'; 
 y=(ey+k*ed(:,[2 4 6 8]))';
 xc=[x; x(1,:)];
 yc=[y; y(1,:)]; 
 % ********* 2D 8-node quadratic elements *********
 elseif nen==8 
 x=(ex+k*ed(:,[1 3 5 7 9 11 13 15])); 
 y=(ey+k*ed(:,[2 4 6 8 10 12 14 16]));
% xc=[x(1); x(5); x(2); x(6); x(3); x(7); x(4); x(8);x(1)];
% yc=[y(1); y(5); y(2); y(6); y(3); y(7); y(4); y(8);y(1)]; 
%
% isoparametric elements
%
 t=-1;
 n=0;
 for s=-1:0.4:1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3
,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%
 s=1;
 n=0;
 for t=-1:0.4:1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3 ,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%
 t=1;
 n=0;
 for s=1:-0.4:-1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3 ,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%
 s=-1;
 n=0;
 for t=1:-0.4:-1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3 ,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%
%********************************************************** 
 else
 disp('Sorry, this element is currently not supported!') 
 return
 end
% ************* plot commands *******************
 axis('equal')
 hold on
 plot(xc,yc,s1)
 plot(x,y,s2)
 hold off
%--------------------------end--------------------------------
calfem_eldraw2.m
function eldraw2(ex,ey,plotpar,elnum)
%eldraw2(ex,ey,plotpar,elnum)
%eldraw2(ex,ey,plotpar)
%eldraw2(ex,ey)
%-------------------------------------------------------------
% PURPOSE 
% Draw the undeformed 2D mesh for a number of elements of 
% the same type. Supported elements are:
%
% 1) -> bar element 2) -> beam el. 
% 3) -> triangular 3 node el. 4) -> quadrilateral 4 node el. 
% 5) -> 8-node isopar. elemen
%
% INPUT 
% ex,ey:.......... nen: number of element nodes
% nel: number of elements
% plotpar=[ linetype, linecolor, nodemark]
% 
% linetype=1 -> solid linecolor=1 -> black
% 2 -> dashed 2 -> blue
% 3 -> dotted 3 -> magenta
% 4 -> red
%
% nodemark=1 -> circle 
% 2 -> star 
% 0 -> no mark 
% 
% elnum=edof(:,1) ; i.e. the first column in the topology matrix
% 
% Rem. Default is solid white lines with circles at nodes.
% 
%-------------------------------------------------------------
% LAST MODIFIED: J Lindemann 1999-01-29
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
%
 if ~((nargin==2)|(nargin==3)|(nargin==4))
 disp('??? Wrong number of input arguments!')
 return
 end
 a=size(ex); b=size(ey);
 
 if (a-b)==[0 0]
 nel=a(1);nen=a(2);
 else
 disp('??? Check size of coordinate input arguments!')
 return
 end
 if nargin==2; 
 plotpar=[1 1 1]; 
 end
 [s1,s2]=calfem_pltstyle(plotpar);
 
% ************************************************** 
% ************* plot coordinates *******************
% ************************************************** 
 x0=sum(ex')/nen; y0=sum(ey')/nen;
 
% ********** Bar or Beam elements *************
 if nen==2 
 x=ex' ;
 y=ey'; 
 xc=x ;yc=y; 
 
% ********* 2D triangular elements ************
 elseif nen==3 
 x=ex' ;
 y=ey';
 xc=[x ; x(1,:)]; yc=[y ; y(1,:)];
 
% ********* 2D quadrilateral elements *********
 elseif nen==4 
 x=ex' ;
 y=ey'; 
 xc=[x ; x(1,:)]; yc=[y ; y(1,:)];
% ********* 2D 8 node quadratic elements *********
 elseif nen==8
 x=ex;
 y=ey; 
% xc=[x(1);x(5);x(2);x(6);x(3);x(7);x(4);x(8);x(1)];
% yc=[y(1);y(5);y(2);y(6);y(3);y(7);y(4);y(8);y(1)];
%
% isoparametric elements
%
 t=-1;
 n=0;
 for s=-1:0.4:1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3 ,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%
 s=1;
 n=0;
 for t=-1:0.4:1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3 ,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%
 t=1;
 n=0;
 for s=1:-0.4:-1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3 ,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%
 s=-1;
 n=0;
 for t=1:-0.4:-1
 n=n+1;
 N1=-1/4*(1-t)*(1-s)*(1+t+s);
 N2=-1/4*(1+t)*(1-s)*(1-t+s);
 N3=-1/4*(1+t)*(1+s)*(1-t-s);
 N4=-1/4*(1-t)*(1+s)*(1+t-s);
 N5=1/2*(1-t*t)*(1-s);
 N6=1/2*(1+t)*(1-s*s);
 N7=1/2*(1-t*t)*(1+s);
 N8=1/2*(1-t)*(1-s*s);
 N=[ N1, N2, N3 ,N4, N5, N6, N7, N8 ];
 x1(n,:)=N*x';
 y1(n,:)=N*y';
 end;
 xc=[xc x1];
 yc=[yc y1];
 clear x1
 clear y1
%********************************************************** 
 else
 disp('!!!! Sorry, this element is currently not supported!') 
 return
 end
% ************************************************** 
% **************** plot commands *******************
% ************************************************** 
 axis('equal')
 hold on
 plot(xc,yc,s1) 
 plot(x,y,s2)
 if nargin==4
 for i=1:nel
 h=text(x0(i),y0(i),int2str(elnum(i)));
 set(h,'fontsize',8);
 end
 end
 xlabel('x'); ylabel('y');
 hold off 
%--------------------------end--------------------------------
calfem_extract.m
 function [ed]=extract(edof,a)
% ed=extract(edof,a)
%-------------------------------------------------------------
% PURPOSE
% Extract element displacements from the global displacement
% vector according to the topology matrix edof.
%
% INPUT: a: the global displacement vector
%
% edof: topology matrix
%
% OUTPUT: ed: element displacement matrix
%-------------------------------------------------------------
% LAST MODIFIED: M Ristinmaa 1993-08-24
% Copyright (c) 1991-94 by Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
 [nie,n]=size(edof);
%
 t=edof(:,2:n);
%
 for i = 1:nie
 ed(i,1:(n-1))=a(t(i,:))';
 end
%
%--------------------------end--------------------------------
calfem_freqresp.m
 function [Resp,Freq]=freqresp(d,dt)
% [Resp,Freq]=freqresp(d,dt)
%----------------------------------------------------------
% PURPOSE
% Calculate the fourier transform of a function
% and plots the result in a loglog scale.
%
% INPUT:
% d : time history function to be transformed
% dt: time step used in the creation of d
%
% OUTPUT:
% Resp: response in frequncy domain
% Freq: corresponding frequencies
%
%----------------------------------------------------------
% LAST MODIFIED: G Sandberg 1993-11-07
% THIAGO RITTO FEV/05
%
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%----------------------------------------------------------
nr=length(d);
if (nr<=80); f2=64; end
if (nr>80 & nr<=150); f2=128; end
if (nr>150 & nr<=350); f2=256; end
if (nr>350 & nr<=850); f2=512; end 
if (nr>850 & nr<=1550); f2=1024; end 
if (nr>1550); f2=2048; end 
ff2=f2/2;
Y=fft(d,f2);
Py=Y.*conj(Y)/f2;
Resp=Py(1:ff2);
Freq=(1/dt)*(0:ff2-1)/f2;
plot(Freq,Resp,'*-b')
%--------------------------end--------------------------------
calfem_pltstyle.m
function [s1,s2]=pltstyle(plotpar)
%-------------------------------------------------------------
% PURPOSE 
% Define define linetype,linecolor and markertype character codes. 
%
% INPUT 
% plotpar=[ linetype, linecolor, nodemark ]
% 
% linetype=1 -> solid linecolor=1 -> black
% 2 -> dashed 2 -> blue
% 3 -> dotted 3 -> magenta
% 4 -> red
%
% nodemark=1 -> circle 
% 2 -> star 
% 0 -> no mark 
% OUTPUT
% s1: linetype and color for mesh lines
% s2: type and color for node markers
%-------------------------------------------------------------
% LAST MODIFIED: Jonas Lindemann 1999-01-29 
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
%
 if plotpar(1)==1 ; s1='-';
 elseif plotpar(1)==2 ; s1='--';
 elseif plotpar(1)==3 ; s1=':';
 else disp('??? Error in variable plotpar(1)!');
 return;
 end
 
 if plotpar(2)==1 ; s1=[s1,'k'];
 elseif plotpar(2)==2 ; s1=[s1,'b'];
 elseif plotpar(2)==3 ; s1=[s1,'m'];
 elseif plotpar(2)==4 ; s1=[s1,'r'];
 else disp('??? Error in variable plotpar(2)!');
 return;
 end
 if plotpar(3)==1 ; s2='ko';
 elseif plotpar(3)==2 ; s2='k*';
 elseif plotpar(3)==0 ; s2='k.';
 else disp('??? Error in variable plotpar(3)!');
 return;
 end
%--------------------------end-------------------------------- 
 
modal_barra_ef.m
% -------------------------------------------------------------------
% ANALISE MODAL DE UMA BARRA
% - APROXIMAÇÃO POR ELEMENTOS FINITOS - 
% -------------------------------------------------------------------
% PROPOSITO
% Calcular modos de vibraçao e frequencias naturais de uma barra,
% discretizada pelo MEF, dada uma certa precisao:
%
% ------------------
%
% |------------------
%
% |------------------|
%
% |------------------ -/\/\/- |
%
% ||||
% |------------------ ||||
% ||||
%
% ENTRADAS:
% L : comprimento da barra (m)
% h :			altura da barra (m)
% b :			espessura da barra (m)
%
% E :			modulo de elasticidade (Pa)
% rho :		densidade (Kg/m3)
% kk : cte da mola da extremidade (N/m)
% mm : massa na extremidade (Ns/m)
%
% cc : condicoes de contorno:
% cc=0 --> livre-livre
% cc=1 --> fixa-livre
% cc=2 --> fixa-fixa
%
% e :			precisão especificada (%)
% N : 		Número de modos que se quer obter com a precisão especificada
%
% NE_ini :	número de elementos iniciais (número PAR!!)
% passo :		passo do número de elementos (número PAR!!)
%
% SAIDAS:
% R_Freq : frequencias naturais (Hz)
% RMat_modo : modos de vibraçao
% NE_prec : numero de elementos necessarios para a precisao requerida
%-------------------------------------------------------------
% -------------------------------------------------------------------
% Autor: THIAGO RITTO
% Ultima atualizacao: MARCO/07
%--------------------------------------------------------------------
clear
clc
close all
% -------------------------------------------------------------------
% INÍCIO da Entrada de Dados
% -------------------------------------------------------------------
L=1;						% comprimento da barra
h = 5e-2;					% altura da barra
b = 10e-2;					% espessura da barra
A = b*h;
E = 2e11;					% modulo de elasticidade
rho = 7850;					% densidade
kk = 0; % cte da mola da extremidade
mm = 0; % massa na extremidade
cc=2; % condicoes de contorno:
 % cc=0 --> livre-livre
 % cc=1 --> fixa-livre
 % cc=2 --> fixa-fixa
N = 6; 						% Número de modos em que se quer obter com a precisão especificada
e = .2;		 		 % precisão especificada --> erro percentual admitido para o N-esimo modo
NE_ini=10;					% número de elementos iniciais (número PAR!!)
passo=10;					% passo do número de elementos (número PAR!!)
% -------------------------------------------------------------------
% FIM da Entrada de Dados
% -------------------------------------------------------------------
fig=0;						% auxiliar para numeração das figuras
 
if NE_ini < 4
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Número de Elementos Iniciais muito pequeno!')
 disp(' ')
 break
end
if rem(NE_ini,2) > 0
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Entrar com número de elementos iniciais PAR!')
 disp(' ')
 break
elseif rem(passo,2) > 0
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Entrar com passos dos elementos PAR!')
 disp(' ')
 break
end
% No loop a seguir aumenta-se o número de elemetos e calcula-se o erro entre as freq nat
% calculadas analiticamente com as calculadas por EF
flag=1;
NE=NE_ini; % número de elementos
erro=100;
while erro > e
 	% ------ topologia -----------------------------------------------
	m=0;
	clear RMat_Edof1 RMat_Edof2
	for j=1:NE
	 R_Edof1(j,1)=j;
	 for k=1:2;	% 2 representa o número de graus de liberdade do elemento
	 RMat_Edof2(j,k)=k+m;
	 end
	 m=m+1;
	end
	RMat_Edof=[R_Edof1 RMat_Edof2];
	% ------ lista de coordenadas -----------------------------------
 clear RMat_Coord
	RMat_Coord=[[0;R_Edof1/NE*L] zeros(NE+1,1)];
 
	% ------ lista de graus de liberdade -----------------------------
	clear RMat_Dof
	m=0;
	for j=1:NE+1
	 for k=1:1 % 2 graus de liberdade por noh
	 RMat_Dof(j,k)=k+m;
	 end
	 m=m+1;
	end
	% ------ gerando matrizes dos elementos, montando nas matrizes globais
	clear RMat_K RMat_M RMat_Ex RMat_Ey
	RMat_K=zeros(RMat_Edof(NE,3)); RMat_M=zeros(RMat_Edof(NE,3));
	[RMat_Ex,RMat_Ey]=calfem_coordxtr(RMat_Edof,RMat_Coord,RMat_Dof,2);		% coordxtr.m --> arq do CALFEM
	
	R_ep=[E A rho];
	
	clear RMat_k RMat_m
	for i=1:NE
	 [RMat_k,RMat_m]=calfem_bar2d(RMat_Ex(i,:),RMat_Ey(i,:),R_ep);		% beam2d.m --> arq do CALFEM
 RMat_K=calfem_assem(RMat_Edof(i,:),RMat_K,RMat_k); RMat_M=calfem_assem(RMat_Edof(i,:),RMat_M,RMat_m); 	% assem.m --> arq do CALFEM
	end
 % ----- Graus de liberdade restritos --------------------
[nd,nd]=size(RMat_M);
 % Condicoes de contorno
 if cc==0
 R_b=[]';
 elseif cc==1
 R_b=[1]';
 elseif cc==2
 R_b=[1 nd]';
 end
 % Corrigindo matrizes M e K para incluir as condicoes de contorno (mola
 % ou massa) na extremidade.
 RMat_K(length(RMat_K),length(RMat_K))=RMat_K(length(RMat_K),length(RMat_K))+kk;
 RMat_M(length(RMat_M),length(RMat_M))=RMat_M(length(RMat_M),length(RMat_M))+mm;
 
 % ----------- PROBLEMA DE AUTOVALOR -----------------------------------
 [R_La,RMat_modo]=calfem_eigen(RMat_K,RMat_M,R_b);		% eigen.m --> arq do CALFEM
 R_Freq1=sqrt(R_La)/(2*pi); % mudanca de unidades rd/s -> Hz
 R_Freq(:,flag)=R_Freq1(1:N);
 
 % ----------- Erro da aproximação ----------
 if flag > 1
 erro(flag-1) = 100*abs(R_Freq(N,flag)-R_Freq(N,flag-1))/R_Freq(N,flag);
 end
 R_Ns(flag)=NE;
 flag=flag+1;
 NE=NE+passo;
 
end
% ---------------
% POS-PROCESSAMENTO - VISUALIZACAO DOS RESULTADOS
% ---------------
disp('Frequências Naturais EF (Hz)')	% N primeiras frequencias naturais calculadas por EF
disp(R_Freq(1:N,flag-1))	% N primeiras frequencias naturais calculadas por EF
disp(' ')
disp('Número de elementos necessários para a precisão requerida')	% número de elementos
NE_prec=NE-passo;
disp(NE_prec)	% número de elementos
disp(' ')
disp('Erro da aproximacao (%)')
disp(erro(flag-2))
MM=round(N/2);
R_X = linspace(0, L, (NE-passo)+1);
fig=fig+1;
figure(fig)
for i=1:N
 subplot(MM,2,i)
 plot(R_X,RMat_modo(:,i),'linewidth',2)
 xlabel('x','fontsize',16)
 text(-.2,0,'u(x)','fontsize',10)
 if i==1
 title('Modos de Vibração','fontsize',16)
 end
 grid on
end
if flag > 2
 fig=fig+1;
 figure(fig)
 plot(R_Ns(2:length(R_Ns)),erro,'linewidth',2)
 ylabel('Erro (%)','fontsize',16)
 xlabel('Numero de elementos','fontsize',16)
 grid on
end
modal_L_ef.m
% -------------------------------------------------------------------
% ANALISE MODAL DE UMA ESTRUTURA EM L
% - APROXIMAÇÃO POR ELEMENTOS FINITOS - 
% -------------------------------------------------------------------
%
% PROPOSITO
% Calcular os modos de vibraçao e frequencias naturais de uma estrutura em
% L discretizada pelo MEF, dada uma certa precisao.
%
% -
% |
% |
% |
% |
% |
% ------------
%
%
% ENTRADAS:
% L : comprimento da barra (m)
% h :			altura da barra (m)
% b :			espessura da barra (m)
%
% E :			modulo de elasticidade (Pa)
% rho :		densidade (Kg/m3)
%
% e :			precisão especificada (%)
% N : 		Número de modos que se quer obter com a precisão especificada
% NE_ini :	número de elementos iniciais (número PAR!!)
%
% SAIDAS:
% R_Freq : frequencias naturais (Hz)
% RMat_modo : modos de vibraçao
% NE_prec : numero de elementos necessarios para a precisao requerida
%-------------------------------------------------------------
% -------------------------------------------------------------------
% Autor: THIAGO RITTO
% Ultima atualizacao: MARCO/07
%--------------------------------------------------------------------
clear
clc
close all
% -------------------------------------------------------------------
% INÍCIO da Entrada de Dados
% -------------------------------------------------------------------
%-------------------------------------------------------------------------
% Definição das variáveis geométricas e do material
%-------------------------------------------------------------------------
L=3;						% comprimento da viga
h = 5e-2;					% altura da viga
b = 10e-2;					% espessura da viga
E=2e11;						% modulo de elasticidade
rho=7850;					% densidade
A = b*h;					% área da seção transversal
I = (b*h^3)/12;			 % momento de inércia da viga
N=6; 						% Número de modos que se quer obter com a precisão especificada
e=.005;					 % precisão especificada --> erro percentual admitido para o N-esimo modo
NE_ini=10;					% número de elementos iniciais (número PAR!!)
%obs. a cada iteracao dobra-se o numero de elementos
% -------------------------------------------------------------------
% FIM da Entrada de Dados
% -------------------------------------------------------------------
if NE_ini < 4
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Número de Elementos Iniciais muito pequeno!')
 disp(' ')
 break
end
if rem(NE_ini,2) > 0
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Entrar com número de elementos iniciais PAR!')
 disp(' ')
 break
end
fig=0;						% auxiliar para numeração das figuras
% No loop a seguir aumenta-se o número de elemetos e calcula-se o erro entre as freq nat
% calculadas analiticamente com as calculadas por EF
flag=1;
NE=NE_ini; % número de elementos
erro=100;
while erro > e
 	% ------ topologia -----------------------------------------------
	m=0;
	clear RMat_Edof1 RMat_Edof2
	for j=1:NE
	 RMat_Edof1(j,1)=j;
	 for k=1:6;	% 6 representa o número de graus de liberdade do elemento
	 RMat_Edof2(j,k)=k+m;
	 end
	 m=m+3;
	end
	RMat_Edof=[RMat_Edof1 RMat_Edof2];
	 
	% ------ lista de coordenadas -----------------------------------
 clear RMat_Coord RMat_Edof3
 RMat_Edof3=flipud(RMat_Edof1(1:NE/2));
 RMat_Coord=[zeros(NE/2,1) [RMat_Edof3(1:NE/2)/NE*L*2]
 	[0; RMat_Edof1(1:NE/2)/(NE)*L*2] zeros(NE/2+1,1)] ;
 
	
	% ------ lista de graus de liberdade -----------------------------
	clear RMat_Dof
	m=0;
	for j=1:NE+1
	 for k=1:3
	 RMat_Dof(j,k)=k+m;
	 end
	 m=m+3;
	end
	
	% ------ gerando matrizes dos elementos, montando nas matrizes globais
	clear RMat_K RMat_M RMat_Ex RMat_Ey
	RMat_K=zeros(RMat_Edof(NE,7)); RMat_M=zeros(RMat_Edof(NE,7));
	[RMat_Ex,RMat_Ey]=calfem_coordxtr(RMat_Edof,RMat_Coord,RMat_Dof,2);		% coordxtr.m --> arq do CALFEM
	
	R_ep=[E A I rho*A 0 0]; % 0 0 corresponde a alfa=0 e beta=0, C=alfa.M + beta.K
	
	clear RMat_k RMat_m RMat_c
	for i=1:NE
	 [RMat_k,RMat_m,RMat_c]=calfem_beam2d(RMat_Ex(i,:),RMat_Ey(i,:),R_ep);		% beam2d.m --> arq do CALFEM
 RMat_K=calfem_assem(RMat_Edof(i,:),RMat_K,RMat_k); RMat_M=calfem_assem(RMat_Edof(i,:),RMat_M,RMat_m);
	end
 
 % ----- Graus de liberdade restritos --------------------
 [nd,nd]=size(RMat_M);
 R_b=[1 2 3]';
 
 % ----------- PROBLEMA DE AUTOVALOR -----------------------------------
 [R_La,RMat_modo]=calfem_eigen(RMat_K,RMat_M,R_b);		% eigen.m --> arq do CALFEM
 R_Freq1=sqrt(R_La)/(2*pi); % mudanca de unidades rd/s -> Hz
 R_Freq(:,flag)=R_Freq1(1:N);
 
 % ----------- Erro da aproximação ----------
 if flag > 1
 erro(flag-1) = 100*abs(R_Freq(N,flag)-R_Freq(N,flag-1))/R_Freq(N,flag);
 end
 R_Ns(flag)=NE;
 flag=flag+1;
 passo = NE;
 NE = NE+passo;
 clear fdof 
end
clc
disp('Frequências Naturais')	% cinco primeiras frequencias naturais calculadas por EF
disp(R_Freq(1:N,flag-1))	% cinco primeiras frequencias naturais calculadas por EF
disp(' ')
disp('Número de elementos necessários para a precisão requerida')	% número de elementos
disp(NE-passo)	% número de elementos
disp(' ')
disp('Erro da aproximacao (Hz)')
disp(erro(flag-2))
fig=fig+1;
figure(fig)
title('Modos de Vibração')
MM=round(N/2);
for i=1:N
 subplot(MM,2,i)
 calfem_eldraw2(RMat_Ex,RMat_Ey,[2 3 1]); % eldraw2.m e eldisp2.m --> arq do CALFEM
 Edb=calfem_extract(RMat_Edof,RMat_modo(:,i));
 calfem_eldisp2(RMat_Ex,RMat_Ey,Edb,[1 2 2],2);
 if i==1
 title('Modos de Vibração','fontsize',16)
 xlabel('x','fontsize',16)
 end
 grid on
end
if flag > 2
 fig=fig+1;
 figure(fig)
 plot(R_Ns(2:length(R_Ns)),erro,'linewidth',2)
 ylabel('Erro (%)','fontsize',16)
 xlabel('Numero de elementos','fontsize',16)
 grid on
end
modal_triang_ef.m
 % -------------------------------------------------------------------
% ANALISE MODAL DE UMA ESTRUTURA COM CONFIGURACAO TRIANGULO
% - APROXIMAÇÃO POR ELEMENTOS FINITOS - 
% -------------------------------------------------------------------
% PROPOSITO
% Calcular os modos de vibraçao e frequencias naturais de uma estrutura com,
% configuraçao triangulo discretizada pelo MEF, dada uma certa precisao.
%
% \/
% /\ 
% / \
% / \
% / \
% / \
% / \
% / ---------
%
%
% ENTRADAS:
% L : comprimento da barra (m)
% h :			altura da barra (m)
% b :			espessura da barra (m)
%
% E :			modulo de elasticidade (Pa)
% rho :		densidade (Kg/m3)
%
% e :			precisão especificada (%)
% N : 		Número de modos que se quer obter com a precisão especificada
% NE_ini :	número de elementos iniciais (múltiplo de TRÊS!!)
%
% SAIDAS:
% R_Freq : frequencias naturais (Hz)
% RMat_modo : modos de vibraçao
% NE_prec : numero de elementos necessarios para a precisao requerida
%-------------------------------------------------------------
% -------------------------------------------------------------------
% Autor: THIAGO RITTO
% Ultima atualizacao: MARCO/07
%--------------------------------------------------------------------
clear
clc
close all
% -------------------------------------------------------------------
% INÍCIO da Entrada de Dados
% -------------------------------------------------------------------
%-------------------------------------------------------------------------
% Definição das variáveis geométricas e do material
%-------------------------------------------------------------------------
L=3;						% comprimento da viga
h = 5e-2;					% altura da viga
b = 10e-2;					% espessura da viga
E=2e11;						% modulo de elasticidade
rho=7850;					% densidade
A = b*h;					% área da seção transversal
I = (b*h^3)/12; 			% momento de inércia da viga
N=6; 						% Número de modos que se quer obter com a precisão especificada
e=1e-1;						% precisão especificada --> erro percentual admitido para o N-esimo modo
NE_ini=15;					% número de elementos iniciais (múltiplo de TRÊS!!)
%obs. a cada iteracao dobra-se o numero de elementos
% -------------------------------------------------------------------
% FIM da Entrada de Dados
% -------------------------------------------------------------------
if NE_ini < 6
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Número de Elementos Iniciais muito pequeno!')
 disp(' ')
 break
end
if rem(NE_ini,3) ~= 0
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Entrar com número de elementos iniciais MÚLTIPLO de 3!')
 disp(' ')
 break
end
fig=0;						% auxiliar para numeração das figuras
% No loop a seguir aumenta-se o número de elemetos e calcula-se o erro entre as freq nat
% calculadas analiticamente com as calculadas por EF
flag=1;
NE=NE_ini; % número de elementos
erro=100;
while erro > e
 	% ------ topologia -----------------------------------------------
	m=0;
	clear RMat_Edof1 RMat_Edof2
	for j=1:NE
	 RMat_Edof1(j,1)=j;
	 for k=1:6;	% 6 representa o número de graus de liberdade do elemento
	 RMat_Edof2(j,k)=k+m;
	 end
	 m=m+3;
	end
	RMat_Edof=[RMat_Edof1 RMat_Edof2];
 
 % ------ lista de coordenadas -----------------------------------
 % Para fazer as coordenadas estou fazendo cada parte do triangulo com o
 %mesmo número de elementos e a parte de baixo (a que não fecha) 70% do tamanho[
 %do lado do triangulo. Além disso estou considerando o triangulo isosceles e o ângulo 
 %de cima de 50 graus.
 
 clear RMat_Coord1 RMat_Coord2 RMat_Coord RMat_Edof3
 for i=1:NE/3
 RMat_Coord1(i,1)=i*L/(NE/3)*cos(((180-45)/2)*pi/180);
 RMat_Coord1(i,2)=i*L/(NE/3)*sin(((180-45)/2)*pi/180);
 end
 aux=NE/3-1;
 for i=NE/3+1:2*NE/3
 RMat_Coord1(i,1)=i*L/(NE/3)*cos(65*pi/180);
 RMat_Coord1(i,2)=aux*L/(NE/3)*sin(65*pi/180);
 aux=aux-1;
 end
 for i=1:NE/3
 RMat_Coord2(i,1)=8*L/4*cos(65*pi/180)-i/(NE/3)*0.7*2*L*sin(25*pi/180);
 RMat_Coord2(i,2)=0;
 end
 
 RMat_Coord=[0 	0
 RMat_Coord1
 RMat_Coord2];
 
 % ------ lista de graus de liberdade -----------------------------
	clear RMat_Dof
	m=0;
	for j=1:NE+1
	 for k=1:3
	 RMat_Dof(j,k)=k+m;
	 end
	 m=m+3;
	end
	
	% ------ gerando matrizes dos elementos, montando nas matrizes globais
	clear RMat_K RMat_M RMat_Ex RMat_Ey
	RMat_K=zeros(RMat_Edof(NE,7)); RMat_M=zeros(RMat_Edof(NE,7));
	[RMat_Ex,RMat_Ey]=calfem_coordxtr(RMat_Edof,RMat_Coord,RMat_Dof,2);		% coordxtr.m --> arq do CALFEM
	
	R_ep=[E A I rho*A 0 0]; % 0 0 corresponde a alfa=0 e beta=0, C=alfa.M + beta.K
	
	clear RMat_k RMat_m RMat_c
	for i=1:NE
	 [RMat_k,RMat_m,RMat_c]=calfem_beam2d(RMat_Ex(i,:),RMat_Ey(i,:),R_ep);		% beam2d.m --> arq do CALFEM
 RMat_K=calfem_assem(RMat_Edof(i,:),RMat_K,RMat_k); RMat_M=calfem_assem(RMat_Edof(i,:),RMat_M,RMat_m);
	end
	
	% ----- Graus de liberdade restritos --------------------
 [nd,nd]=size(RMat_M);
 R_b=[NE+1 NE+2 ]';
 
 % ----------- PROBLEMA DE AUTOVALOR -----------------------------------
 [R_La,RMat_modo]=calfem_eigen(RMat_K,RMat_M,R_b);		% eigen.m --> arq do CALFEM
 R_Freq1=sqrt(R_La)/(2*pi); % mudanca de unidades rd/s -> Hz
 R_Freq(:,flag)=R_Freq1(1:N);
 
 % ----------- Erro da aproximação ----------
 if flag > 1
 erro(flag-1) = 100*abs(R_Freq(N,flag)-R_Freq(N,flag-1))/R_Freq(N,flag);
 end
 R_Ns(flag)=NE;
 flag=flag+1;
 passo = NE;
 NE = NE+passo;
 clear fdof 
end
clc
disp('Frequências Naturais')	% cinco primeiras frequencias naturais calculadas por EF
disp(R_Freq(1:N,flag-1))	% cinco primeiras frequencias naturais calculadas por EF
disp(' ')
disp('Número de elementos necessários para a precisão requerida')	% número de elementos
disp(NE-passo)	% número de elementos
disp(' ')
disp('Erro da aproximacao (Hz)')
disp(erro(flag-2))
fig=fig+1;
figure(fig)
title('Modos de Vibração')
MM=round(N/2);
for i=1:N
 subplot(MM,2,i)
 calfem_eldraw2(RMat_Ex,RMat_Ey,[2 3 1]); % eldraw2.m e eldisp2.m --> arq do CALFEM
 Edb=calfem_extract(RMat_Edof,RMat_modo(:,i));
 calfem_eldisp2(RMat_Ex,RMat_Ey,Edb,[1 2 2],2);
 if i==1
 title('Modos de Vibração','fontsize',16)
 xlabel('x','fontsize',16)
 end
 grid on
end
if flag > 2
 fig=fig+1;
 figure(fig)
 plot(R_Ns(2:length(R_Ns)),erro,'linewidth',2)
 ylabel('Erro (%)','fontsize',16)
 xlabel('Numero de elementos','fontsize',16)
 grid on
end
modal_viga_ef.m
% -------------------------------------------------------------------
% ANALISE MODAL DE UMA VIGA
% - APROXIMAÇÃO POR ELEMENTOS FINITOS - 
% -------------------------------------------------------------------
% PROPOSITO
% Calcular modos de vibraçao e frequencias naturais de uma viga,
% discretizada pelo MEF, dada uma certa precisao:
%
% ------------------
%
% |------------------
%
% |------------------|
%
% ------------------
% /\ /\
%
% ENTRADAS:
% L : comprimento da barra (m)
% h :		 altura da barra (m)
% b :		 espessura da barra (m)
%
% E :		 modulo de elasticidade (Pa)
% rho :		densidade (Kg/m3)
%
% cc : condices de contorno:
% cc=0 --> livre-livre
% cc=1 --> engastada-livre
% cc=2 --> engastada-engastada
% cc=3 --> bi-apoiada
%
% e :			precisão especificada (%)
% N : 		Número de modos que se quer obter com a precisão especificada
% NE_ini :	número de elementos iniciais (número PAR!!)
%
% SAIDAS:
% R_Freq :
frequencias naturais (Hz)
% RMat_modo : modos de vibraçao
% NE_prec : numero de elementos necessarios para a precisao requerida
%-------------------------------------------------------------
% -------------------------------------------------------------------
% Autor: THIAGO RITTO
% Ultima atualizacao: MARCO/07
%--------------------------------------------------------------------
clear
clc
close all
% -------------------------------------------------------------------
% INÍCIO da Entrada de Dados
% -------------------------------------------------------------------
L = 3;						% comprimento da viga
h = 5e-2;					% altura da viga
b = 10e-2;					% espessura da viga
E = 2e11;					% modulo de elasticidade
rho = 7850;					% densidade
cc=3; % condicoes de contorno:
 % cc=0 --> livre-livre
 % cc=1 --> engastada-livre
 % cc=2 --> engastada-engastada
 % cc=3 --> bi-apoiada
N=6; 	 % Número de modos que se quer obter com a precisão especificada
e=1e-2;					 % precisão especificada --> erro percentual admitido para o N-esimo modo
NE_ini=20;					% número de elementos iniciais
%obs. a cada iteracao dobra-se o numero de elementos
% -------------------------------------------------------------------
% FIM da Entrada de Dados
% -------------------------------------------------------------------
if N > 20
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Número modos muito grande!')
 disp(' ')
 break
end
fig=0;						% auxiliar para numeração das figuras
A = b*h;						% área da seção transversal
I = (b*h^3)/12;			 % momento de inércia da viga
delta_x=L/1000;			 % a viga é discretizada com intervalos delta_x
if NE_ini < 4
 clc
 disp(' ')
 disp(' ERRO ')
 disp(' ')
 disp('Número de Elementos Iniciais muito pequeno!')
 disp(' ')
 break
end
% No loop a seguir aumenta-se o número de elemetos e calcula-se o erro entre as freq nat
% calculadas analiticamente com as calculadas por EF
flag=1;
NE=NE_ini; % número de elementos
erro=100;
while erro > e
 
 % ------ topologia -----------------------------------------------
	m=0;
	clear RMat_Edof1 RMat_Edof2
	for j=1:NE
	 RMat_Edof1(j,1)=j;
	 for k=1:6;	% 6 representa o número de graus de liberdade do elemento
	 RMat_Edof2(j,k)=k+m;
	 end
	 m=m+3;
	end
	RMat_Edof=[RMat_Edof1 RMat_Edof2];
	 
	% ------ lista de coordenadas -----------------------------------
	clear RMat_Coord
	RMat_Coord=[[0;RMat_Edof1/NE*L] zeros(NE+1,1)];
	
	% ------ lista de graus de liberdade -----------------------------
	clear Dof
	m=0;
	for j=1:NE+1
	 for k=1:3
	 RMat_Dof(j,k)=k+m;
	 end
	 m=m+3;
	end
	
	% ------ gerando matrizes dos elementos, montando nas matrizes globais
	clear RMat_K RMat_M RMat_Ex RMat_Ey
	RMat_K=zeros(RMat_Edof(NE,7)); RMat_M=zeros(RMat_Edof(NE,7));
	[RMat_Ex,RMat_Ey]=calfem_coordxtr(RMat_Edof,RMat_Coord,RMat_Dof,2);		% calfem_coordxtr.m --> arq do CALFEM
	
	R_ep=[E A I rho*A 0 0]; % 0 0 corresponde a alfa=0 e beta=0, C=alfa.M + beta.K
 
	clear RMat_k RMat_m RMat_c
 for i=1:NE
 [RMat_k,RMat_m,RMat_c]=calfem_beam2d(RMat_Ex(i,:),RMat_Ey(i,:),R_ep);		% beam2d.m --> arq do CALFEM
 RMat_K=calfem_assem(RMat_Edof(i,:),RMat_K,RMat_k); RMat_M=calfem_assem(RMat_Edof(i,:),RMat_M,RMat_m); 
 end
 % ----- Graus de liberdade restritos --------------------
 [nd,nd]=size(RMat_M);
 % Condicoes de contorno
 if cc==0
 R_b=[]';
 elseif cc==1
 R_b=[1 2 3]';
 elseif cc==2
 R_b=[1 2 3 (nd-2) (nd-1) (nd)]';
 elseif cc==3
 R_b=[1 2 (nd-2) (nd-1)]';
 end
 
 % ----------- PROBLEMA DE AUTOVALOR -----------------------------------
 [R_La,RMat_modo]=calfem_eigen(RMat_K,RMat_M,R_b);		% eigen.m --> arq do CALFEM
 R_Freq1=sqrt(R_La)/(2*pi); % mudanca de unidades rd/s -> Hz
 R_Freq(:,flag)=R_Freq1(1:N);
 % ----------- Erro da aproximação ----------
 if flag > 1
 erro(flag-1) = 100*abs(R_Freq(N,flag)-R_Freq(N,flag-1))/R_Freq(N,flag);
 end
 R_Ns(flag)=NE;
 flag=flag+1;
 passo = NE;
 NE = NE+passo;
 clear fdof 
 
end
% ---------------
% POS-PROCESSAMENTO - VISUALIZACAO DOS RESULTADOS
% ---------------
%clc
disp('Frequências Naturais calculadas por Elementos Finitos')	% cinco primeiras frequencias naturais calculadas por EF
disp(R_Freq(1:N,flag-1))	% cinco primeiras frequencias naturais calculadas por EF
disp(' ')
disp('Número de elementos necessários para a precisão requerida')	% número de elementos
NE_prec = NE-passo;
disp(NE_prec)	% número de elementos
disp(' ')
disp('Erro da aproximacao (Hz)')
disp(erro(flag-2))
fig=fig+1;
figure(fig)
title('Modos de Vibração')
MM=round(N/2);
for i=1:N
 subplot(MM,2,i)
 calfem_eldraw2(RMat_Ex,RMat_Ey,[2 3 1]); % eldraw2.m e eldisp2.m --> arq do CALFEM
 Edb=calfem_extract(RMat_Edof,RMat_modo(:,i));
 calfem_eldisp2(RMat_Ex,RMat_Ey,Edb,[1 2 2],5);
 if i==1
 title('Modos de Vibração','fontsize',16)
 xlabel('x','fontsize',16)
 end
 grid on
end
if flag > 2
 fig=fig+1;
 figure(fig)
 plot(R_Ns(2:length(R_Ns)),erro,'linewidth',2)
 xlabel('Erro (%)','fontsize',16)
 xlabel('Numero de elementos','fontsize',16)
 grid on
end
read me.txt
Programas da apostila "Análise Modal", RUBENS SAMPAIO e THIAGO RITTO
modal_barra_ef.m
Analise modal de uma barra discretizada pelo Metodo dos Elementos Finitos
modal_viga_ef.m
Analise modal de uma viga discretizada pelo Metodo dos Elementos Finitos
modal_l_ef.m
Analise modal de uma estrutura em L discretizada pelo Metodo dos Elementos Finitos
modal_triang_ef.m
Analise modal de uma estrutura em triangulo discretizada pelo Metodo dos Elementos Finitos
calfem_assem.m
 function [K,f]=assem(edof,K,Ke,f,fe)
% K=assem(edof,K,Ke)
% [K,f]=assem(edof,K,Ke,f,fe)
%-------------------------------------------------------------
% PURPOSE
% Assemble element matrices Ke ( and fe ) into the global
% stiffness matrix K ( and the global force vector f )
% according to the topology matrix edof.
%
% INPUT: edof: dof topology matrix
% K : the global stiffness matrix
% Ke: element stiffness matrix
% f : the global force vector
% fe: element force vector
%
% OUTPUT: K : the new global stiffness matrix
% f : the new global force vector
%-------------------------------------------------------------
% LAST MODIFIED: M Ristinmaa 1993-10-06
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
 [nie,n]=size(edof);
 t=edof(:,2:n);
 for i = 1:nie
 K(t(i,:),t(i,:)) = K(t(i,:),t(i,:))+Ke;
 if nargin==5
 f(t(i,:))=f(t(i,:))+fe;
 end
 end
%--------------------------end--------------------------------
calfem_bar2d.m
function [Ke,Me]=bar2e(ex,ey,ep)
% Ke=bar2e(ex,ey,ep)
%----------------------------------------------------------------------
% PURPOSE
% Compute the element stiffness matrix for two dimensional bar element.
% E tambem a matriz de massa do elemento...
%
% INPUT: ex = [x1 x2];
% ey = [y1 y2]; element node coordinates
%
% ep = [E A rho] E: Young's modulus rho: densidade
% A: Cross section area
%
% OUTPUT: Ke : stiffness matrix, dim(Ke)= 4 x 4
%----------------------------------------------------------------------
% LAST MODIFIED: K Persson 1995-08-23
% Copyright (c) Division of Structural Mechanics and
% Department of Solid Mechanics.
% Lund Institute of Technology
%-------------------------------------------------------------
 E=ep(1); A=ep(2); rho=ep(3); 
 b=[ ex(2)-ex(1); ey(2)-ey(1) ];
 L=sqrt(b'*b);
 Kle=E*A/L*[ 1 -1; 
 -1 1];
 Mle=rho*A*L/6*[ 2 1; 
 1 2];
 
 n=b'/L;
 
 G=[ n(1) n(2); 
 -n(2) n(1) ];
 Ke=G'*Kle*G;
 Me=G'*Mle*G;
%--------------------------end--------------------------------
Apostila MEF (Rubens - PUC).pdf
Vibrac¸o˜es Mecaˆnicas
Dinaˆmica de estruturas flex´ıveis
Rubens Sampaio, Priscilla Almeida, Thiago Ritto
Rio de Janeiro, marc¸o de 2007.
Conteu´do
1 INTRODUC¸A˜O 2
1.1 Objetivo do Trabalho 3
1.2 Organizac¸a˜o do Trabalho 5
2 MODELAGEM DE PROBLEMAS DE BARRAS 7
2.1 Dinaˆmica de Barras 7
2.2 Problema Modelo: resoluc¸a˜o por separac¸a˜o de varia´veis 10
3 FORMULAC¸A˜O FRACA 18
3.1 Vantagens da Formulac¸a˜o Fraca 19
3.2 Formulac¸a˜o Fraca e Me´todo Variacional de Aproximac¸a˜o 20
3.3 Formulac¸a˜o Fraca: problemas de barras 23
4 MODELAGEM DE PROBLEMAS DE VIGAS 40
4.1 Dinaˆmica de Vigas 40
4.2 Problema Modelo: resoluc¸a˜o por separac¸a˜o de varia´veis 42
5 FORMULAC¸A˜O FRACA PARA PROBLEMAS DE VIGAS 47
5.1 Problema de uma Viga Livre-Livre 47
5.2 Problema de uma Viga Engastada-Livre 52
5.3 Problema de uma Viga Engastada-Mola 54
6 ME´TODO DE ELEMENTOS FINITOS - MEF 60
6.1 Aproximac¸a˜o do dom´ınio 66
6.2 Aproximac¸a˜o da soluc¸a˜o no dom´ınio aproximado 66
6.3 Escolha do nu´mero de elementos iniciais (N) 80
6.4 Ana´lise do erro 81
7 APLICAC¸A˜O DO MEF PROBLEMAS DE BARRAS E VIGAS 82
7.1 Problema de uma barra fixa-livre 82
7.2 Viga engastada com massa concentrada na extremidade 87
7.3 Viga engastada com mola vertical na extremidade livre 91
8 VIGAS COM CONDIC¸O˜ES INTERMEDIA´RIAS 95
8.1 Problemas de vigas com um apoio intermedia´rio 95
8.2 Problema de viga engastada-livre com dois apoios intermedia´rios 102
8.3 Viga bi-engastada com uma mola vertical em coordenada intermedia´ria105
8.4 Viga bi-apoiada com massa concentrada em coordenada intermedia´ria106
8.5 Problema de viga apoiada-livre com uma mola torcional 108
8.6 Aproximac¸a˜o da dinaˆmica de um problema de viga com apoio
intermedia´rio 109
9 CONCLUSA˜O 116
10 PRINC´IPIO DE HAMILTON: PROBLEMAS DE BARRAS E VIGAS 117
10.1 Princ´ıpio de Hamilton: Barras 117
10.2 Princ´ıpio de Hamilton: Vigas 119
Bibliografia 123
A MANUAL DE PROGRAMAS DO MATLAB 125
A.1 Simulac¸o˜es de problemas de barras 128
A.2 Simulac¸o˜es de Vigas 144
1
INTRODUC¸A˜O
O Me´todo de Elementos Finitos (MEF) pode ser definido como uma
ferramenta nume´rica para obtenc¸a˜o da aproximac¸a˜o de um problema [14];
e consiste na discretizac¸a˜o de um meio cont´ınuo, dividindo-o em elementos.
Esses elementos sa˜o descritos por equac¸o˜es diferenciais e resolvidos por
modelos matema´ticos para que sejam obtidos os resultados desejados [3].
O desenvolvimento do MEF originou-se no final do se´culo XVIII,
quando Gauss propoˆs a utilizac¸a˜o de func¸o˜es de aproximac¸a˜o para a
soluc¸a˜o de problemas matema´ticos. Durante mais de um se´culo, diversos
matema´ticos desenvolveram teorias e te´cnicas anal´ıticas para a soluc¸a˜o
de problemas, entretanto, pouco se evoluiu devido a` dificuldade no
processamento de equac¸o˜es alge´bricas. O desenvolvimento pra´tico desta
ana´lise ocorreu mais tarde, em consequ¨eˆncia dos avanc¸os tecnolo´gicos,
por volta de 1950. Isto permitiu a elaborac¸a˜o e a resoluc¸a˜o de sistemas
de equac¸o˜es complexas [5]. Em 1960, Turner, Clough, Martins e Topp,
propuseram um me´todo de ana´lise estrutural, o MEF, e descreveram-no.
Figura 1.1: Fluxograma representativo do MEF
O fluxograma de (1.1) representa, esquematicamente, o funcionamento
de um modelo baseado em simulac¸o˜es [9]. Historicamente, os me´todos
matriciais para ana´lise estrutural (Matrix Structural Analysis - MSA), que
eram utilizados antes de 1950, antecedem o MEF.
A MSA e o MEF sa˜o baseados em treˆs aspectos principais [9]:
– Modelagem matema´tica;
– Formulac¸a˜o matricial de equac¸o˜es discretas lineares;
– Ferramentas computacionais.
Figura 1.2: Principais aspectos da MSA e do MEF
A partir da de´cada de 70, o MEF passou a ser aplicado tambe´m em
problemas de mecaˆnica dos fluidos.
1.1
Objetivo do Trabalho
Em geral, a dinaˆmica de uma estrutura e´ descrita por um sistema
de equac¸o˜es diferenciais parciais em conjunto com condic¸o˜es de contorno,
condic¸o˜es intermedia´rias e condic¸o˜es iniciais. Uma vez conhecidas as
frequ¨eˆncias naturais e os modos normais de vibrac¸a˜o do sistema (exatamente
ou aproximadamente), pode-se controlar sua dinaˆmica [10]. O MEF pode
ser utilizado para aproximar, numericamente, as frequ¨eˆncias e os modos da
estrutura; e, para isso, uma formulac¸a˜o fraca do problema e´ prefer´ıvel.
Este trabalho tem por objetivo formar um material dida´tico que sera´
destinado aos alunos de graduac¸a˜o, introduzindo conceitos ba´sicos sobre o
estudo de sistemas cont´ınuos e sua discretizac¸a˜o pelo me´todo de elementos
finitos. Pretende-se explicar o que e´ esse me´todo de aproximac¸a˜o, como
funciona e qual deve ser o procedimento para calcular aproximac¸o˜es de
sistemas cont´ınuos.
Uma ide´ia e´ na˜o usar pacotes de programac¸a˜o, mas sim desenvolver
programas usando o Matlab e estimular alunos de graduac¸a˜o a
programarem. Ale´m disso, deseja-se incentivar os alunos a se preocuparem
com noc¸o˜es de aproximac¸a˜o, erro de aproximac¸a˜o e convergeˆncia.
Deseja-se construir aproximac¸o˜es para a dinaˆmica de problemas
formados por estruturas simples unidimensionais (barras ou combinac¸a˜o de
barras, vigas ou combinac¸a˜o de vigas). Isso sera´ feito nas seguintes etapas:
– Modelagem do problema na Formulac¸a˜o Fraca: parte-se da formulac¸a˜o
forte do problema [7], pore´m a modelagem pode ser feita diretamente
usando o princ´ıpio de Hamilton;
– Construc¸a˜o de uma base de aproximac¸a˜o: de forma a representar a
soluc¸a˜o como uN(x, t) =
∑N
1 ai(t)φi(x) [7].
Caso os modos de vibrac¸a˜o do problema em estudo na˜o sejam
conhecidos, o que acontece na maioria dos casos, constro´i-se uma
aproximac¸a˜o para eles, atrave´s do Me´todo de Elementos Finitos (resolvendo
o problema de valor caracter´ıstico associado); e, em seguida, constro´i-se um
modelo reduzido com a aproximac¸a˜o dos modos [7]. A maior dificuldade
concentra-se na construc¸a˜o da aproximac¸a˜o, que deve satisfazer um crite´rio
de erro pre´-fixado (por exemplo, construir os modos com um crite´rio de erro
controla´vel).
Com a aproximac¸a˜o u(x, t) =
∑
ai(t)φi(x), deseja-se aproximar a
dinaˆmica do problema; e isso pode ser dividido em dois problemas:
1. Achar uma base de aproximac¸a˜o (aproximar φi);
2. Aproximar a dinaˆmica do sistema (aproximar ai).
O processo de soluc¸a˜o [4] pode ser visualizado pela figura (1.3) e sera´
exemplificado no decorrer do trabalho.:
Figura 1.3: Processo de soluc¸a˜o
1.2
Organizac¸a˜o do Trabalho
No cap´ıtulo 2 sera´ apresentada a modelagem de um problema de barra
[8], de forma que sera´ desenvolvida a equac¸a˜o representativa da dinaˆmica de
movimento desse sistema e, em seguida, sera´ calculada a soluc¸a˜o anal´ıtica
de um problema espec´ıfico (barra livre-livre), apresentando as frequ¨eˆncias
naturais e os modos de vibrac¸a˜o correspondentes.
O cap´ıtulo 3 apresentara´ todo o procedimento necessa´rio para calcular
a formulac¸a˜o
fraca de um sistema cont´ınuo [11], mostrando as considerac¸o˜es
necessa´rias, as vantagens desta te´cnica e as caracter´ısticas do me´todo
variacional adotado. Em seguida, sera´ desenvolvida a formulac¸a˜o fraca de
um problema de barra, considerando diferentes condic¸o˜es de contorno.
Os assuntos abordados para problemas de barras, nos cap´ıtulos 2
e 3, sera˜o apresentados para problemas de vigas nos cap´ıtulos 4 e 5,
respectivamente. No cap´ıtulo 4, sera´ modelado um problema de viga para
calcular a dinaˆmica de movimento do sistema e sera´ desenvolvida a soluc¸a˜o
anal´ıtica de um problema de viga bi-engastada [8]; enquanto, no cap´ıtulo 5,
sera´ desenvolvida a formulac¸a˜o fraca de diferentes problemas de vigas, com
diferentes condic¸o˜es de contorno.
O cap´ıtulo 6 introduzira´ conceitos ba´sicos do Me´todo de Elementos
Finitos, mostrando o funcionamento do me´todo e todas as etapas necessa´rias
para a sua aplicac¸a˜o [1]. Sera˜o apresentados diferentes tipos de elementos,
que definem o tipo de aproximac¸a˜o; ale´m da construc¸a˜o de equac¸o˜es
elementares que, acopladas, formam as equac¸o˜es globais do sistema.
O cap´ıtulo 7 mostrara´ a aplicac¸a˜o do MEF a problemas espec´ıficos
de barras e vigas. Nele sera˜o apresentadas as matrizes de massa e rigidez
elementares bem como as matrizes globais correspondentes [17].
Ja´ no cap´ıtulo 8, sera˜o apresentados problemas de estruturas que
apresentam condic¸o˜es intermedia´rias, sejam apoios, molas ou massas
concentradas [?]. Para esses problemas, sera˜o desenvolvidas as formulac¸o˜es
fracas correspondentes e sera˜o aproximadas as frequ¨eˆncias naturais, os
modos de vibrac¸a˜o e a dinaˆmica do sistema.
Em anexo, sera´ apresentado um Manual dos Programas desenvolvidos
no Matlab, mostrando todos os passos necessa´rios. Esses programas
aproximam frequ¨eˆncias naturais e modos de vibrac¸a˜o de diversos problemas,
atrave´s do me´todo de Elementos Finitos; e depois aproximam a dinaˆmica
do sistema.
2
MODELAGEM DE PROBLEMAS DE BARRAS
A barra e´ uma estrutura mecaˆnica que, num dado instante, pode ser
descrita por um u´nico paraˆmetro, o deslocamento longitudinal; portanto
e´ um objeto unidimensional [17]. Este cap´ıtulo apresentara´ a dinaˆmica de
movimento e a soluc¸a˜o de problemas envolvendo barras.
2.1
Dinaˆmica de Barras
Considere um corpo unidimensional, representado pelo intervalo [0,L],
com densidade ρ, a´rea de sec¸a˜o transversal A e mo´dulo de elasticidade E
[8]. Seja u(x, t) a posic¸a˜o do ponto x no instante t e f(x, t) a forc¸a externa
ao sistema. u(•, t) e´ a configurac¸a˜o do corpo no instante t e u(x, •) e´ a
trajeto´ria do ponto x. Um exemplo de barra (fixa em uma extremidade e
livre na outra) esta´ apresentado na figura (2.1):
Figura 2.1: Barra fixa-livre
Tomando uma parte finita do domı´nio [x1x2] e adotando as leis de
conservac¸a˜o de quantidade de movimento [15], deseja-se calcular a dinaˆmica
do corpo:
Figura 2.2: Sec¸a˜o [x1x2] do domı´nio
A sec¸a˜o [x1x2] esta´ submetida a forc¸as internas P (x, t) que dependem da
deformac¸a˜o local:
lim
x2→x1
u(x2,t)−u(x1,t)
x2−x1 =
∂u
∂x
∣∣
x=x1
(2-1)
Dessa forma, pode-se escrever a dinaˆmica em [x1x2]
∂
∂t
∫ x2
x1
ρ(x)A(x)
∂u
∂t
(x, t)dx = P (x2, t)− P (x1, t) +
∫ x2
x1
f(x, t)dx (2-2)
Define-se 4x = x2 − x1 e, se x2 −→ x1, enta˜o 4x −→ 0.
Usando a equac¸a˜o constitutiva que relaciona P (x, t) com a deformac¸a˜o
∂u
∂x
(x, t), e caracteriza materiais ela´sticos, lineares (lei de Hooke), definida
por:
P (x,t)
A(x)
= E(x)ε(x, t) = E(x)∂u
∂x
(x, t)
P (x, t) = E(x)A(x)∂u
∂x
(x, t)
(2-3)
E(x) e´ o mo´dulo de elasticidade do material, no ponto x; e ε(x,t) e´ a
deformac¸a˜o.
Aproxima-se a integral
∫ b
a
g(x)dx ∼ g(a)(b − a), e substitui-se a equac¸a˜o
(2-3) em (2-2):
ρ(x1)A(x1)
∂2u
∂t2
∣∣∣∣∣
x2
x1
4x = E(x2)A(x2)∂u
∂x
(x2, t)−E(x1)A(x1)∂u
∂x
(x1, t)+f(x1)4x
(2-4)
Divide-se toda a equac¸a˜o (2-4) por 4x e faz-se lim(4x→ 0)
ρ(x)A(x)
∂2u
∂t2
(x, t)− ∂
∂x
(E(x)A(x)
∂u
∂x
(x, t)) = f(x, t) (2-5)
Considerando constantes a sec¸a˜o transversal da barra e as propriedades do
material, a equac¸a˜o (2-5) reduz-se a:
ρA
∂2u
∂t2
(x, t)− EA∂
2u
∂x2
(x, t) = f(x, t) (2-6)
Para f = 0:
∂2u
∂t2
=
E
ρ
∂2u
∂x2
(2-7)
A equac¸a˜o (2-5) fornece parte da informac¸a˜o sobre a dinaˆmica do
corpo, porque faltam as informac¸o˜es referentes a`s extremidades 0 e L, que
sa˜o denominadas condic¸o˜es de contorno [8]. Sa˜o necessa´rias, tambe´m, as
condic¸o˜es iniciais, que informam a posic¸a˜o e a velocidade de todos os pontos;
ou seja, as condic¸o˜es da dinaˆmica no in´ıcio do processo.
A equac¸a˜o (2-5), juntamente com as condic¸o˜es de contorno e as
condic¸o˜es iniciais constituem um problema que pode-se mostrar ter soluc¸a˜o
u´nica. Alguns exemplos de problemas de barras com diferentes condic¸o˜es de
contorno [8], esta˜o apresentados na figura (2.3):
Figura 2.3: Condic¸o˜es de contorno e configurac¸o˜es associadas a cada tipo
de barra
Definic¸a˜o das condic¸o˜es de contorno:
– Extremidade fixa: u(x, t) = 0;
– Extremidade livre: P (x, t) = EA∂u
∂x
(x, t) = 0;
– Extremidade com mola (ke): P (x, t) = EA
∂u
∂x
(x, t) = −keu(x, t);
– Extremidade com massa (me): P (x, t) = EA
∂u
∂x
(x, t) = me
∂2u
∂t2
(x, t).
2.2
Problema Modelo: resoluc¸a˜o por separac¸a˜o de varia´veis
Escolheu-se o problema de uma barra livre-livre como modelo a ser
apresentado. Como esse problema tem configurac¸a˜o muito simples, a soluc¸a˜o
de frequ¨eˆncias naturais e modos de vibrac¸a˜o correspondentes encontram-se
na literatura [8]. O ca´lculo da soluc¸a˜o anal´ıtica do problema pode dividido
em treˆs etapas, que sera˜o explicadas a seguir.
Atrave´s desse problema-modelo (problema de barra livre-livre),
deseja-se introduzir os conceitos de frequ¨eˆncia natural e modo de vibrac¸a˜o;
ale´m de mostrar que a aproximac¸a˜o da forma
∑N
1 ai(t)φi(x) e´ realista [8].
Considere uma barra homogeˆnea (mesmo material) e uniforme (mesma
geometria); de comprimento L, livre nas duas extremidades e com forc¸a
externa nula, f = 0.
Figura 2.4: Barra livre-livre
Conforme foi apresentado na figura (2.3), quando a extremidade da
barra e´ livre, P (x, t) = 0. Com isso, pode-se especificar o problema da barra
livre-livre atrave´s da equac¸a˜o de movimento, das condic¸o˜es de contorno e
das condic¸o˜es iniciais; que sa˜o, respectivamente:
Problema(P )

∂2u
∂t2
(x, t) = c2 ∂
2u
∂x2
(x, t) ∴ c2 = E
ρ
, em (0, L)
EA∂u
∂x
(0, t) = 0 EA∂u
∂x
(L, t) = 0
u(x, 0) = u0(x)
∂u
∂t
(x, 0) = v0(x)
(2-8)
Para cada ponto no intervalo [0, L], tem-se uma condic¸a˜o: nos pontos
internos (0, L) a condic¸a˜o e´ prescrita por uma equac¸a˜o diferencial parcial e
nos extremos, (x = 0 e x = L), a condic¸a˜o de a barra estar ”livre”, faz com
que a forc¸a exercida seja nula. As condic¸o˜es iniciais do problema informam
a posic¸a˜o e a velocidade de cada ponto, no in´ıcio do processo.
Deseja-se encontrar a soluc¸a˜o u da equac¸a˜o diferencial da barra, que
satisfac¸a as condic¸o˜es de contorno e as condic¸o˜es iniciais (2-8). Para isso,
devem ser realizadas treˆs etapas [8]: Aplicac¸a˜o do Me´todo de Separac¸a˜o
de Varia´veis, a partir do qual sera˜o obtidas duas Equac¸o˜es Diferenciais
Ordina´rias (EDOs). Uma vez
determinadas as Soluc¸o˜es das duas EDOs, que
satisfac¸am as condic¸o˜es de contorno do problema; sera´ feita a Superposic¸a˜o
das Soluc¸o˜es de forma a satisfazer as condic¸o˜es iniciais.
– Primeira etapa: Me´todo de Separac¸a˜o de Varia´veis
Parte-se da equac¸a˜o diferencial e das condic¸o˜es de contorno,
apresentadas em (2-8); e procura-se uma soluc¸a˜o u(x, t) que seja o
produto de duas func¸o˜es, uma dependente apenas da posic¸a˜o x e outra
dependente somente do instante t:
u(x, t) = X(x)T (t) (2-9)
Nessa primeira etapa, na˜o utiliza-se as condic¸o˜es iniciais do problema.
Define-se diferentes nomenclaturas adotadas para diferenciac¸a˜o em
relac¸a˜o a` coordenada x e em relac¸a˜o a` t, quando trata-se de func¸o˜es
de uma varia´vel [8]:
∂
∂x
=′ ∂
∂t
= ˙
Substitui-se (2-9) na equac¸a˜o diferencial de (2-8):
XT¨ = c2X ′′T (2-10)
a equac¸a˜o (2-10) pode ser reescrita por:
X ′′
X
(x) =
T¨
c2T
(t) (2-11)
Como o lado esquerdo da equac¸a˜o (2-11) depende apenas da varia´vel
x e o lado esquerdo depende apenas de t, conclui-se que os dois lados
devem ser iguais a uma constante que, por convenieˆncia, sera´ chamada
de −λ2.
A equac¸a˜o (2-11) pode ser separada em duas equac¸o˜es diferenciais
ordina´rias:
X ′′ + λ2X = 0 (2-12)
T¨ + c2λ2T = 0 (2-13)
– Segunda etapa: Determinac¸a˜o das soluc¸o˜es das duas EDOs
A mesma separac¸a˜o de varia´veis realizada na equac¸a˜o (2-9) deve ser
aplicada nas condic¸o˜es de contorno de (2-8):
∂u
∂x
(0, t) = 0 =⇒ X ′(0)T (t) = 0
∂u
∂x
(L, t) = 0 =⇒ X ′(L)T (t) = 0
(2-14)
Para T (t) ≡ 0, tem-se que u(x, t) = 0, o que na˜o satisfaz condic¸o˜es
iniciais na˜o-nulas. Por isso, exige-se que:
X ′(0) = 0 e X ′(L) = 0 (2-15)
Encontra-se, enta˜o o seguinte problema de valor de contorno,
conhecido como um problema de Sturm-Liouville [8]:
X ′′ + λ2X = 0
(2-16)
X ′(0) = 0 X ′(L) = 0
Deseja-se encontrar o par (λ,X) tal que X 6= 0 e que satisfac¸a o
problema de valor caracter´ıstico apresentado pela equac¸a˜o (2-12),
satisfazendo as condic¸o˜es de contorno (2-15).
λ sa˜o os autovalores e X as autofunc¸o˜es correspondentes.
Para o caso de λ2 = 0, tem-se queX ′′ = 0 e a soluc¸a˜oX(x) = A1x+A2.
Quando sa˜o impostas as condic¸o˜es de contorno de (2-15), verifica-se
que X(x) = 0; por isso a condic¸a˜o de X
′′
X
= 0 e´ descartada.
Tem-se que a soluc¸a˜o geral da equac¸a˜o (2-12) e´:
X(x) = B1 cosλx+B2 sinλx (2-17)
e a respectiva derivada e´:
X ′(x) = −B1λ sinλx+B2λ cosλx (2-18)
A primeira condic¸a˜o de contorno X ′(0) = 0 e´ substitu´ıda em (2-18).
Sendo sin(0) = 0 e cos(0) = 1, a equac¸a˜o (2-18) reduz-se a:
B2λ = 0 =⇒ B2 = 0 (2-19)
A poss´ıvel soluc¸a˜o X(x) reduz-se a:
X(x) = B1 cosλx (2-20)
Pela segunda condic¸a˜o de contorno, X ′(L) = 0:
−B1λ sinλL = 0 =⇒ sinλL = 0 (2-21)
Da equac¸a˜o (2-21), existiriam treˆs opc¸o˜es: B1 = 0, λ = 0, sinλL = 0.
Conforme foi explicado anteriormente, λ 6= 0; e se B1 = 0, pela
equac¸a˜o (2-20), X(x) tambe´m seria nulo; portanto, necessariamente,
sinλL = 0. Essa condic¸a˜o impo˜e restric¸o˜es a` constante λ. Por motivos
que sera˜o explicados posteriormente, e´ necessa´rio achar todas as
soluc¸o˜es do problema, de forma a satisfazerem as condic¸o˜es iniciais
do mesmo (terceira etapa) [8].
Existem infinitas possibilidade que satisfazem a condic¸a˜o sin(λL) = 0
e, para representar cada uma delas utiliza-se um ı´ndice n diferente;
ou seja λnL = npi, tal que n = 1,2 ... ∞. Enta˜o:
λn =
npi
L
n = 1, 2... (2-22)
Para cada λn, existira´ uma soluc¸a˜o Xn(x) a ele associada:
Xn(x) = B1n cosλnx n = 1, 2... (2-23)
Conforme foi visto pela equac¸a˜o (2-23), existem infinitas soluc¸o˜es para
a primeira func¸a˜o Xn(x). O mesmo acontece com T (t) que, a partir
de (2-13), tem-se T¨n + λ
2
nc
2Tn = 0; tal que, uma soluc¸a˜o para Tn(t) e´:
Tn(t) = C1n cos(λnct) + C2n sin(λnct) (2-24)
Partindo-se de (2-9) tem-se que a soluc¸a˜o candidata do sistema e´:
un(x, t) = Xn(x)Tn(t) ∴ n = 1, 2... (2-25)
que, juntamente com (2-23) e (2-24) pode ser reescrita por:
un(x, t) = cos(λnx)(D1n cos(λnct) +D2n sin(λnct)) (2-26)
sendo D1n = B1nC1n e D2n = B1nC2n
– Terceira etapa: Superposic¸a˜o de soluc¸o˜es
As infinitas soluc¸o˜es do problema na˜o necessariamente satisfazem as
condic¸o˜es iniciais do problema (2-8). Para garantir que essas condic¸o˜es
sejam satisfeitas, faz-se a superposic¸a˜o das soluc¸o˜es candidatas (un);
o que resulta em:
u(x, t) =
∞∑
n=0
cos(λnx)(D1n cos(λnct) +D2n sin(λnct)) (2-27)
sendo λnc a frequ¨eˆncia natural do sistema.
D1n e D2n sa˜o calculadas a partir das condic¸o˜es iniciais, apresentadas
em (2-8). Substituindo-se a primeira delas na equac¸a˜o (2-27), tem-se:
u(x, 0) =
∞∑
n=0
D1n cos(λnx) = u0(x) (2-28)
multiplica-se os dois lados da equac¸a˜o por cos(λmx) e integra-os no
domı´nio [0 L]:
∫ L
0
∞∑
n=0
D1n cos(λnx) cos(λmx)dx =
∫ L
0
u0(x) cos(λmx)dx (2-29)
o resultado da integral do somato´rio acima e´ L/2 para n = m e 0 para
n 6= m, enta˜o:
D1n =
2
L
∫ L
0
u0(x) cos(λnx)dx (2-30)
Aplica-se a segunda condic¸a˜o inicial do sistema (2-8) e tem-se:
∞∑
n=0
D2ncλn cos(λnx) = v0(x) (2-31)
A mesma te´cnica de multiplicar toda a equac¸a˜o por cos(λmx), para
enta˜o integra´-la no domı´nio [0 L] e´ adotada:
∫ L
0
∞∑
n=0
D2ncλn cos(λnx) cos(λmx)dx =
∫ L
0
v0(x) cos(λmx)dx (2-32)
D2n =
2
Lλnc
∫ L
0
v0(x) cos(λnx)dx (2-33)
A soluc¸a˜o geral da dinaˆmica de uma barra livre-livre, com as
respectivas condic¸o˜es de contorno e condic¸o˜es iniciais, consiste na
substituic¸a˜o de D1n e D2n (2-30 e 2-33) na expressa˜o de u, de (2-27).
A soluc¸a˜o u pode ser expressa de forma simplificada [14]:
u(x, t) =
∞∑
n=1
φn(x)an(t) (2-34)
sendo:
φn(x) = cos(λnx)
an(t) = D1n cos(λnct) +D2n sin(λnct)
Associado a cada frequ¨eˆncia natural, existe um Modo de Vibrac¸a˜o, que
e´ uma configurac¸a˜o do sistema quando todos os seus pontos vibram
sincronicamente, ou seja, na mesma frequ¨eˆncia [13].
As frequ¨eˆncias naturais e os modos de vibrac¸a˜o deste sistema sa˜o,
respectivamente:
Frequ¨eˆncias Naturais: ωn = cλn = c
npi
L
(2-35)
Modos de Vibrac¸a˜o: φn(x) = cos
npix
L
(2-36)
n = 1,2,3...
u pode ser aproximada por N modos de vibrac¸a˜o, sendo representada
por uN(x, t) e essa aproximac¸a˜o gera uma parcela de erro a ela
associada, que e´ εN(x, t):
u(x, t) = uN(x, t) + εN(x, t) (2-37)
∞∑
n=1
φn(x)an(t)︸ ︷︷ ︸
u(x,t)
=
N∑
n=1
φn(x)an(t)︸ ︷︷ ︸
uN (x,t)
+
∞∑
n=N+1
φn(x)an(t)︸ ︷︷ ︸
εN (x,t)
(2-38)
aproximac¸a˜o: uN(x, t) =
N∑
n=1
φn(x)an(t) (2-39)
erro de aprox: εN(x, t) =
∞∑
n=N+1
φn(x)an(t) (2-40)
sendo N o nu´mero de modos usados na aproximac¸a˜o de uN [14].
3
FORMULAC¸A˜O FRACA
Um problema pode ser representado pela sua Formulac¸a˜o Forte,
que corresponde a`s equac¸a˜o diferenciais parciais, uma para cada ponto
interno do domı´nio, em conjunto com as condic¸o˜es de contorno, as condic¸o˜es
intermedia´rias e as condic¸o˜es iniciais [7] [14]. Por exemplo, a equac¸a˜o (2-8)
e´ a formualc¸a˜o forte do problema de uma barra livre-livre. A Formulac¸a˜o
Fraca e´ uma forma alternativa de representar problemas e consta de uma
equac¸a˜o variacional e um espac¸o de func¸o˜es admiss´ıveis [14].
A equac¸a˜o variacional incorpora as equac¸o˜es diferenciais e parte das