Baixe o app para aproveitar ainda mais
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
Compartilhar