Contents
- CÀLCUL D'INESTABILITAT PER AL PROGRAMA CME2.
- SUBMATRIUS GEOMETRIQUES EN EIXOS LOCALS.
- SUBMATRIUS DE REGIDESA EN EIXOS LOCALS.
- CÀLCUL DE LES SOL·LICITACIONS ALS NUSOS PER CÀRREGUES A LES BARRES.
- APLICACIÓ CONDICIÓ DE CONTORN RODET INCLINAT A LES MATRIUS DE RIGIDESA.
- RESOLUCIÓ DEL PROBLEMA NO LINEAL MITJANÇANT EL MÈTODE D'EULER
- RESOLUCIÓ DEL PROBLEMA NO LINEAL MITJANÇANT EL MÈTODE DE RUNGE_KUTTA
- FUNCIÓ PER REORDENAR ELS MODES DE VINCLAMENT
CÀLCUL D'INESTABILITAT PER AL PROGRAMA CME2.
MANUEL GIL RAMOS 2.013-2.042 GRAU EN ENGINYERIA DE LA CONSTRUCCIÓ UNIVERSITAT POLITÈCNICA DE CATALUNYA
Versió data 10/06/2.014
Entrada: n = sencer nombre de nusos de l'estructura. N = sencer nombre de barres de l'estructura. Ncoord = matriu n*3 amb les coordenades x, y i z dels nusos. Ncond = condicions de contorn dels nusos: nus lliure =0; rodet inclinat =1; recolzament simple =2; encastament =3; desplaçaments elàstics i girs lliures =4; girs elàstics i desplaçaments nuls =5; desplaçaments i girs elàstics =6; desplaçaments imposats i girs lliures =7; girs imposats i desplaçaments lliures =8; desplaçaments i girs imposats =9. Nextfor = vector n*6 de forces aplicades als nusos. Lnods = matriu N*2 amb nombre de nus inferior i superior de cada barra. Lartyp = matriu N*2 amb tipus d'extrem de barra inferior i superior, rígid =0; articulat =1. Lmecp = matriu N*9 de propietats mecaniques de cada barra: E, G, A, Iz, Iy, zx, zy, zz, J. barra_carrega_rep = matriu N*6 amb la càrrega repartida sobre les barres, mòdul, direcció x, y i z, distancia a extrem inferior i superior. barra_carrega_rep_proj = matriu N*6 amb la càrrega repartida projectada sobre les barres (longitut de la càrrega = longitut de la barra), mòdul, direcció x, y i z, distància a extrem inferior i superior. barra_carrega_punt = matriu N*5 amb la càrrega puntual sobre les barres, mòdul, direcció x, y i z i distància al nus inferior. barra_carrega_mom = matriu N*5 amb el moment aplicat sobre les barres, mòdul, direcció x, y i z i distància al nus inferior. iteracions = nombre d'iteracions a realitzar pel mètode iteratiu de resolució de 2on ordre alpha1 = primera variable del mètode de resolució Runge-Kutta alpha2 = segona variable del mètode de resolució Runge-Kutta mu = tercera variable del mètode de resolució Runge-Kutta metode = mètode de resolució a utilitzar pel càlcul de 2on ordre
Sortida: coef = vector amb els coeficients de vinclament ordenats per valor absolut. mode_vincl = matriu amb els moviments del nusos segons el mode de vinclament associat. TL2G = matriu 3*3*N amb les matrius de canvi de base local a global de cada barra. ierr = comptador d'errors que retorna un numero entre el 0 i el 4 segons el tipus d'error
function [coef,mode_vincl,TL2G,ierr] = calculinestabilitat(n,N,Ncoord,... Ncond,Nextfor,Lnods,Lartyp,Lmecp,barra_carrega_rep,... barra_carrega_rep_proj,barra_carrega_punt,barra_carrega_mom,... iteracions,alpha1,alpha2,mu,metode) coef=zeros(1); % Si tots els nusos tenen la mateixa coordenada x dimx=0. dimx=any(Ncoord(:,1)~=Ncoord(1,1)); % Si tots els nusos tenen la mateixa coordenada y dimy=0. dimy=any(Ncoord(:,2)~=Ncoord(1,2)); % Si tots els nusos tenen la mateixa coordenada z dimz=0. dimz=any(Ncoord(:,3)~=Ncoord(1,3)); % Si existeix una barra col·locada en una direcció... if (dimx&~dimy&~dimz)|(~dimx&dimy&~dimz)|(~dimx&~dimy&dimz) % Si la barra té les mateixes coordenades x i z (barra en eix y) if ~any(Ncoord(:,1)~=Ncoord(1,1)) && ~any(Ncoord(:,3)~=Ncoord(1,3)) Ncoord(n+1,1)=Ncoord(1,1)+1; % Creació d'un punt imaginari fora de % l'eix de la barra per evitar % eliminar certs moviments % Si la barra té les mateixes coordenades y i z (barra en eix x) elseif ~any(Ncoord(:,2)~=Ncoord(1,2)) && ~any(Ncoord(:,3)~=Ncoord(1,3)) Ncoord(n+1,2)=Ncoord(1,2)+1; % Creació d'un punt imaginari fora de % l'eix de la barra per evitar % eliminar certs moviments % Si la barra té les mateixes coordenades z (barra en eix x o y) elseif ~any(Ncoord(:,3)~=Ncoord(1,3)) ierr=5; % Com les estructures bidimensionals han d'estar als % eixos x-y, retorna un error TL2G=[]; % Dimensiona la matriu de canvi de base buida mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament % buida end end % Si existeixen càrregues en barra (repartides, moments, puntuals...) if any(barra_carrega_rep(:,1))|any(barra_carrega_rep_proj(:,1))|... any(barra_carrega_punt(:,1))|any(barra_carrega_mom(:,1)) ierr=3; % Error que indica càrregues en barra TL2G=[]; % Dimensiona la matriu de canvi de base buida mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament % buida return; % Retorna l'execució a la funció cme2.m end % Cridem a la funció del càlcul estàtic amb els arguments de sortida del % nostre interès [~,~,~,TL2G,ierr,K_Glob,Vector_esforcos,Vector_esforcos_prescrits,... sistema,ns,nn] = cme2_calcul_estatic(n,N,Ncoord,Ncond,Nextfor,... Lnods,Lartyp,Lmecp,barra_carrega_rep,barra_carrega_rep_proj,... barra_carrega_punt,barra_carrega_mom); K_Geom_Glob_cte=zeros(ns,ns); % Dimensionem la matriu geomètrica global % perquè el programa s'executi mès ràpid if ierr==1 % Si l'execució del càlcul estàtic ens retorna un error... mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament % buida return % Retorna l'execució a la funció cme2.m end kriggeom_cte_tot=zeros(12,12,N); % Dimensiona la hipermatriu de matrius % geomètriques locals en eixos locals kriggeom_cte=zeros(12,12,N); % Dimensiona la hipermatriu de matrius % geomètriques locals en eixos locals for i=1:N % Per a cada barra... n1=Lnods(i,1); % node inferior n2=Lnods(i,2); % node superior extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0 % rígid-articulat =1 % articulat-rígid =2 % articulat-articulat =3 vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra. L=norm(vector_barra); % Longitut de la barra. T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base % completa 6x6 [kriggeom_cte(:,:,i)]=calcul_submatrius_geometriques(Lmecp(i,:),... extrems,L); % Càlcul de les submatrius geomètriques en eixos % locals. kriggeom_cte_tot(:,:,i)=[T*kriggeom_cte(1:6,1:6,i)*T',T*... kriggeom_cte(1:6,7:12,i)*T';T*kriggeom_cte(7:12,1:6,i)*T'... ,T*kriggeom_cte(7:12,7:12,i)*T']; % Submatriu geomètrica en eixos % globals. kriggeom_cte_tot(:,:,i)=contorn_rodet_inclinat(kriggeom_cte_tot... (:,:,i),n1,n2,Ncond); % Aplicació de la condició de contorn de % rodet inclinat a la submatriu geometrica, % si s'escau. end Vector_moviments=zeros(ns,1); % Dimensiona el vector de moviments per % optimitzar l'execució incr=iteracions; % Li dona a la variable incr el valor de les iteracions % introduides a l'interfície gràfica switch metode % Switch per escollir mètode case 1 % Cas del mètode 1 (Euler) [K_Geom_Glob_cte,errmet]=Euler_iterative... (Vector_esforcos_prescrits,K_Glob,TL2G,K_Geom_Glob_cte,... sistema,Vector_moviments,Vector_esforcos,n,nn,ns,Ncond,N,... barra_carrega_rep,barra_carrega_rep_proj,barra_carrega_punt,... barra_carrega_mom,kriggeom_cte_tot,kriggeom_cte,Lnods,... Lartyp,Lmecp,Ncoord,incr); case 2 % Cas del mètode 2 (Runge-Kutta) [K_Geom_Glob_cte,errmet]=Runge_Kutta(Vector_esforcos_prescrits,... K_Glob,TL2G,K_Geom_Glob_cte,sistema,Vector_moviments,... Vector_esforcos,n,nn,ns,Ncond,N,barra_carrega_rep,... barra_carrega_rep_proj,barra_carrega_punt,barra_carrega_mom,... kriggeom_cte_tot,kriggeom_cte,Lnods,Lartyp,Lmecp,Ncoord,... incr,alpha1,alpha2,mu); end if errmet==1 % Si l'execució dels mètodes ens ha retornat l'error de canvi % de signe ierr=4; % Dona a la variable ierr el valor pertinent TL2G=[]; % Dimensiona la matriu de canvi de base buida mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament buida return; % Retorna l'execució a la funció cme2.m end B=K_Geom_Glob_cte(1:nn,1:nn); % Dimensiona la variable B amb la part de % la matriu de rigidesa que ens interessa A=K_Glob(1:nn,1:nn); % Dimensiona la variable A amb la part de la matriu % de rigidesa global que ens interessa [vect,mat_coefcoef]=eig(A,B); % Resolució del problema d'autovalors amb % diverses variables de retorn pels % coeficients i els modes de vinclament dimvec=size(vect); % Es troba el nombre de moviments determinats pel mode % de vinclament if dimvec(1)<6*n % Si el nombre de moviments determinats pel mode de % vinclament és inferior al total... vect((dimvec(1)+1):(6*n),:)=0; % La resta de moviments són igual a 0 end dimvec=size(vect); % Es troba el nombre de moviments totals vect_reord=zeros(dimvec(1),1); for i=1:nn % Bucle que reordena els vectors de la matriu amb els modes % de vinclament. Col·loca primer el desplaçament en x del node % 1, el desplaçament en y, etc. vect_reord(:,i)=reordenacio_vec(vect(:,i),sistema,n,Ncond); end i=1; while i<=nn % Bucle per extreure la diagonal de la matriu de coeficients vect_coef(i,1)=mat_coefcoef(i,i); i=i+1; end vect_coef_abs=abs(vect_coef); % Vector amb els valors absoluts [~,ord]=sort(vect_coef_abs); % Busquem l'ordre del vector absolut for i=1:max(ord) % Bucle que ordena els coeficients i modes de vinclament % segons el valor absolut dels coeficients coeford(i,1)=vect_coef(ord(i),1); vect_vincl(:,i)=vect_reord(:,ord(i)); end if all(isinf(coeford))|all(isnan(coeford)) % Si tots els coeficients % són infinit... ierr=2; % Dona a la variable ierr el valor pertinent mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament buida return; % Retorna l'execució a la funció cme2.m end k=1; for i=1:max(ord) % Bucle que elimina els coeficients que són infinit % (no són d'interès) if ~isinf(coeford(i)) coef(k)=coeford(i); mode_vincl(:,k)=vect_vincl(:,i); k=k+1; end end end
Input argument "Ncoord" is undefined. Error in ==> calculinestabilitat at 58 dimx=any(Ncoord(:,1)~=Ncoord(1,1));
SUBMATRIUS GEOMETRIQUES EN EIXOS LOCALS.
% Funció que calcula les matrius geomètriques locals segons els tipus % d'extrems de les barres. function [kriggeom_cte]=calcul_submatrius_geometriques(Lmecp,extrems,l) switch extrems case 0 %si la barra es rigida en tots dos extrems kriggeom_cte(1,1)=0; %k11 kriggeom_cte(2,2)=36; kriggeom_cte(2,6)=3*l; kriggeom_cte(3,3)=36; kriggeom_cte(3,5)=-3*l; kriggeom_cte(4,4)=0; kriggeom_cte(5,3)=-3*l; kriggeom_cte(5,5)=4*l^2; kriggeom_cte(6,2)=3*l; kriggeom_cte(6,6)=4*l^2; kriggeom_cte(1,7)=0; %k12 kriggeom_cte(2,8)=-36; kriggeom_cte(2,12)=3*l; kriggeom_cte(3,9)=-36; kriggeom_cte(3,11)=-3*l; kriggeom_cte(4,10)=0; kriggeom_cte(5,9)=-3*l; kriggeom_cte(5,11)=-l^2; kriggeom_cte(6,8)=-3*l; kriggeom_cte(6,12)=-l^2; kriggeom_cte(7,1)=0; %k21 kriggeom_cte(8,2)=-36; kriggeom_cte(8,6)=-3*l; kriggeom_cte(9,3)=-36; kriggeom_cte(9,5)=3*l; kriggeom_cte(10,4)=0; kriggeom_cte(11,3)=-3*l; kriggeom_cte(11,5)=-l^2; kriggeom_cte(12,2)=3*l; kriggeom_cte(12,6)=-l^2; kriggeom_cte(7,7)=0; %k22 kriggeom_cte(8,8)=36; kriggeom_cte(8,12)=-3*l; kriggeom_cte(9,9)=36; kriggeom_cte(9,11)=3*l; kriggeom_cte(10,10)=0; kriggeom_cte(11,9)=3*l; kriggeom_cte(11,11)=4*l^2; kriggeom_cte(12,8)=-3*l; kriggeom_cte(12,12)=4*l^2; case 1 %si la barra es articulada en l'extrem inferior kriggeom_cte(1,1)=0; kriggeom_cte(2,2)=36; kriggeom_cte(2,6)=0; kriggeom_cte(3,3)=36; kriggeom_cte(3,5)=0; kriggeom_cte(4,4)=0; kriggeom_cte(5,3)=-0; kriggeom_cte(5,5)=0; kriggeom_cte(6,2)=0; kriggeom_cte(6,6)=0; kriggeom_cte(1,7)=0; kriggeom_cte(2,8)=-36; kriggeom_cte(2,12)=6*l; kriggeom_cte(3,9)=-36; kriggeom_cte(3,11)=-6*l; kriggeom_cte(4,10)=-0; kriggeom_cte(5,9)=0; kriggeom_cte(5,11)=0; kriggeom_cte(6,8)=-0; kriggeom_cte(6,12)=0; kriggeom_cte(7,1)=0; kriggeom_cte(8,2)=-36; kriggeom_cte(8,6)=0; kriggeom_cte(9,3)=-36; kriggeom_cte(9,5)=0; kriggeom_cte(10,4)=-0; kriggeom_cte(11,3)=-6*l; kriggeom_cte(11,5)=0; kriggeom_cte(12,2)=6*l; kriggeom_cte(12,6)=0; kriggeom_cte(7,7)=0; kriggeom_cte(8,8)=36; kriggeom_cte(8,12)=-6*l; kriggeom_cte(9,9)=36; kriggeom_cte(9,11)=6*l; kriggeom_cte(10,10)=0; kriggeom_cte(11,9)=6*l; kriggeom_cte(11,11)=6*l^2; kriggeom_cte(12,8)=-6*l; kriggeom_cte(12,12)=6*l^2; case 2 %si la barra es articulada en l'extrem superior kriggeom_cte(1,1)=0; kriggeom_cte(2,2)=36; kriggeom_cte(2,6)=6*l; kriggeom_cte(3,3)=36; kriggeom_cte(3,5)=-6*l; kriggeom_cte(4,4)=0; kriggeom_cte(5,3)=-6*l; kriggeom_cte(5,5)=6*l^2; kriggeom_cte(6,2)=6*l; kriggeom_cte(6,6)=6*l^2; kriggeom_cte(1,7)=0; kriggeom_cte(2,8)=-36; kriggeom_cte(2,12)=0; kriggeom_cte(3,9)=-36; kriggeom_cte(3,11)=-0; kriggeom_cte(4,10)=-0; kriggeom_cte(5,9)=6*l; kriggeom_cte(5,11)=0; kriggeom_cte(6,8)=-6*l; kriggeom_cte(6,12)=0; kriggeom_cte(7,1)=0; kriggeom_cte(8,2)=-36; kriggeom_cte(8,6)=-6*l; kriggeom_cte(9,3)=-36; kriggeom_cte(9,5)=-6*l; kriggeom_cte(10,4)=-0; kriggeom_cte(11,3)=-0; kriggeom_cte(11,5)=0; kriggeom_cte(12,2)=0; kriggeom_cte(12,6)=0; kriggeom_cte(7,7)=0; kriggeom_cte(8,8)=36; kriggeom_cte(8,12)=-0; kriggeom_cte(9,9)=36; kriggeom_cte(9,11)=0; kriggeom_cte(10,10)=0; kriggeom_cte(11,9)=0; kriggeom_cte(11,11)=0; kriggeom_cte(12,8)=-0; kriggeom_cte(12,12)=0; case 3 %si la barra es articulada en tots dos extrems kriggeom_cte(1,1)=0; kriggeom_cte(2,2)=30; kriggeom_cte(2,6)=0; kriggeom_cte(3,3)=30; kriggeom_cte(3,5)=0; kriggeom_cte(4,4)=0; kriggeom_cte(5,3)=0; kriggeom_cte(5,5)=0; kriggeom_cte(6,2)=0; kriggeom_cte(6,6)=0; kriggeom_cte(1,7)=0; kriggeom_cte(2,8)=-30; kriggeom_cte(2,12)=0; kriggeom_cte(3,9)=-30; kriggeom_cte(3,11)=-0; kriggeom_cte(4,10)=-0; kriggeom_cte(5,9)=0; kriggeom_cte(5,11)=0; kriggeom_cte(6,8)=0; kriggeom_cte(6,12)=0; kriggeom_cte(7,1)=0; kriggeom_cte(8,2)=-30; kriggeom_cte(8,6)=0; kriggeom_cte(9,3)=-30; kriggeom_cte(9,5)=0; kriggeom_cte(10,4)=-0; kriggeom_cte(11,3)=-0; kriggeom_cte(11,5)=0; kriggeom_cte(12,2)=0; kriggeom_cte(12,6)=0; kriggeom_cte(7,7)=0; kriggeom_cte(8,8)=30; kriggeom_cte(8,12)=-0; kriggeom_cte(9,9)=30; kriggeom_cte(9,11)=0; kriggeom_cte(10,10)=0; kriggeom_cte(11,9)=0; kriggeom_cte(11,11)=0; kriggeom_cte(12,8)=-0; kriggeom_cte(12,12)=0; end end
SUBMATRIUS DE REGIDESA EN EIXOS LOCALS.
% Funció que calcula les matrius de rigidesa en eixos locals. Al cridar a % la funció de càlcul estàtic nomès obtenim la matriu de rigidesa global. function KL=calcul_submatrius_rigidesa(Lmecp,extrems,L) E=Lmecp(1); % Mòdul de Young Iz=Lmecp(4); % Inèrcia eix z if any(Lmecp(6:8)) % Si existeix vector direcció d'Iz Iy=Lmecp(5); % Inèrcia eix y else Iy=Iz; end GItL=Lmecp(9)/L; % Inèrcia a torsió A=Lmecp(3); % Àrea de la secció EAL=E*A/L; EIzL=E*Iz/L; EIyL=E*Iy/L; EIzL2=E*Iz/L^2; EIyL2=E*Iy/L^2; EIzL3=E*Iz/L^3; EIyL3=E*Iy/L^3; % Sub-matrius de rigidesa segons els extrems siguin rígids o articulats. switch extrems case 0 % rígid-rígid KL(1:6,1:6)=[EAL,0,0,0,0,0;0,12*EIzL3,0,0,0,6*EIzL2;0,0,... 12*EIyL3,0,-6*EIyL2,0;0,0,0,GItL,0,0;0,0,-6*EIyL2,0,4*... EIyL,0;0,6*EIzL2,0,0,0,4*EIzL]; KL(1:6,7:12)=[-EAL,0,0,0,0,0;0,-12*EIzL3,0,0,0,6*EIzL2;0,0,... -12*EIyL3,0,-6*EIyL2,0;0,0,0,-GItL,0,0;0,0,6*EIyL2,0,2*EIyL,... 0;0,-6*EIzL2,0,0,0,2*EIzL]; KL(7:12,1:6)=KL(1:6,7:12)'; KL(7:12,7:12)=[EAL,0,0,0,0,0;0,12*EIzL3,0,0,0,-6*EIzL2;0,0,... 12*EIyL3,0,6*EIyL2,0;0,0,0,GItL,0,0;0,0,6*EIyL2,0,4*EIyL,0;... 0,-6*EIzL2,0,0,0,4*EIzL]; case 1 % rígid-articulat KL(1:6,1:6)=[EAL,0,0,0,0,0;0,3*EIzL3,0,0,0,3*EIzL2;0,0,3*EIyL3,... 0,-3*EIyL2,0;0,0,0,0,0,0;0,0,-3*EIyL2,0,3*EIyL,0;0,3*EIzL2,... 0,0,0,3*EIzL]; KL(1:6,7:12)=[[-EAL,0,0;0,-3*EIzL3,0;0,0,-3*EIyL3;0,0,0;0,0,3*... EIyL2;0,-3*EIzL2,0],zeros(6,3)]; KL(7:12,1:6)=KL(1:6,7:12)'; KL(7:12,7:12)=[[EAL,0,0;0,3*EIzL3,0;0,0,3*EIyL3],zeros(3,3);... zeros(3,6)]; case 2 % articulat-rígid KL(1:6,1:6)=[[EAL,0,0;0,3*EIzL3,0;0,0,3*EIyL3],zeros(3);... zeros(3,6)]; KL(1:6,7:12)=[[-EAL,0,0,0,0,0;0,-3*EIzL3,0,0,0,3*EIzL2;0,0,... -3*EIyL3,0,-3*EIyL2,0];zeros(3,6)]; KL(7:12,1:6)=KL(1:6,7:12)'; KL(7:12,7:12)=[EAL,0,0,0,0,0;0,3*EIzL3,0,0,0,-3*EIzL2;0,0,... 3*EIyL3,0,3*EIyL2,0;0,0,0,0,0,0;0,0,3*EIyL2,0,3*EIyL,0;0,... -3*EIzL2,0,0,0,3*EIzL]; case 3 % articulat-articulat KL(1:6,1:6)=[EAL,0,0,0,0,0;zeros(5,6)]; KL(1:6,7:12)=[-EAL,0,0,0,0,0;zeros(5,6)]; KL(7:12,1:6)=KL(1:6,7:12); KL(7:12,7:12)=KL(1:6,1:6); end end
CÀLCUL DE LES SOL·LICITACIONS ALS NUSOS PER CÀRREGUES A LES BARRES.
% Per a cada tipus d'extrem de barra, rígid o articulat, en calcula els % esforços d'empotrament perfecte, per a càrregues repartides projectades % (longitut de la càrrega = longitut de la barra) o no sobre la barra, % càrregues puntuals i moments aplicats a la bara. Per a càrregues % repartides permet definir distàncies al nus inferior i superior per % delimitar la zona de la barra que rep la càrrega. Igualment amb la % distància al nus inferior es situa la càrrega puntual i el moment a la % barra. Totes les càrregues s'introdueixen, i els resultats es retornen, % en eixos globals, però es calculen en eixos locals (eix de la barra % coincideix amb l'eix X local). function reaccions=calcul_carregues_barra_local(carrega_rep,... carrega_rep_proj,carrega_punt,carrega_mom,TL2G,L,extrems) reaccions=[0,0,0,0,0,0,0,0,0,0,0,0]; if carrega_rep_proj(1) % Càlcul de les forces equivalents als nusos a la % càrrega repartida projectada sobre la barra. vbarra=[L,0,0]; % Vector barra en eixos locals; carrega=TL2G'*carrega_rep_proj(2:4)'; % Vector direcció de la càrrega % en eixos locals. carrega=carrega_rep_proj(1)*carrega/norm(carrega); % Vector càrrega en % eixos locals. vperp=cross(vbarra,carrega); vperp=cross(vperp,carrega); vperp=vperp/norm(vperp); Lc=abs(dot(vperp,vbarra)); % Longitut de la càrrega. reaccions=-calcul_repartida(carrega,carrega_rep_proj(5),... carrega_rep_proj(6),L,Lc,extrems); end if carrega_rep(1) % Càlcul de les forces equivalents als nusos a la % càrrega repartida a la barra. carrega=TL2G'*carrega_rep(2:4)'; % Vector direcció de la càrrega % repartida en eixos locals. carrega=carrega_rep(1)*carrega'/norm(carrega); % Vector càrrega reaccions=reaccions-calcul_repartida(carrega,carrega_rep(5),... carrega_rep(6),L,L,extrems); end if carrega_punt(1) % Càlcul de les forces equivalents als nusos a la % càrrega puntual a la barra. carrega=TL2G'*carrega_punt(2:4)'; % Vector direcció de la carrega % puntual carrega=carrega'*carrega_punt(1)/norm(carrega); % Vector carrega. reaccions=reaccions-calcul_puntual(carrega,carrega_punt(5),L,extrems); end if carrega_mom(1) % Càlcul de les forces equivalents als nusos al moment % aplicat a la barra. carrega=TL2G'*carrega_mom(2:4)'; % Vector direcció del moment. carrega=carrega'*carrega_mom(1)/norm(carrega); % Vector moment. reaccions=reaccions-calcul_moment(carrega,carrega_mom(5),L,extrems); end end
APLICACIÓ CONDICIÓ DE CONTORN RODET INCLINAT A LES MATRIUS DE RIGIDESA.
% Condició de contorn rodet inclinat. % Si algún dels extrems connecta amb un rodet inclinat, multiplica les % submatrius per la matriu de canvi de base als eixos del rodet, els % moviments i esforços del nus corresponent queden en eixos locals del % rodet. function KG=contorn_rodet_inclinat(KG,n1,n2,Ncond) if Ncond(n1,1)==1 % Si el nus inferior és un rodet inclinat Tr=calcul_canvi_base(Ncond(n1,2:4)); % Calcula matriu canvi de base % del rodet. T=[Tr,zeros(3);zeros(3),Tr]; KG(1:6,1:6)=T'*KG(1:6,1:6)*T; % Sub-matriu en la base del rodet. KG(1:6,7:12)=T'*KG(1:6,7:12); % Sub-matriu en la base del rodet. KG(7:12,1:6)=KG(7:12,1:6)*T; % Sub-matriu en la base del rodet. end if Ncond(n2,1)==1 % Si el nus superior és un rodet inclinat Tr=calcul_canvi_base(Ncond(n2,2:4)); % Calcula matriu canvi de base % del rodet. T=[Tr,zeros(3);zeros(3),Tr]; KG(1:6,7:12)=KG(1:6,7:12)*T; % Sub-matriu en la base del rodet. KG(7:12,1:6)=T'*KG(7:12,1:6); % Sub-matriu en la base del rodet. KG(7:12,7:12)=T'*KG(7:12,7:12)*T; % Sub-matriu en la base del rodet. end end
RESOLUCIÓ DEL PROBLEMA NO LINEAL MITJANÇANT EL MÈTODE D'EULER
% Funció que retorna la matriu geomètrica global calculada pel mètode % d'Euler function [K_Geom_Glob_cte,errmet]=Euler_iterative... (Vector_esforcos_prescrits,K_Glob,TL2G,K_Geom_Glob_cte,... sistema,Vector_moviments,Vector_esforcos,n,nn,ns,Ncond,N,... barra_carrega_rep,barra_carrega_rep_proj,barra_carrega_punt,... barra_carrega_mom,kriggeom_cte_tot,kriggeom_cte,Lnods,Lartyp,Lmecp,... Ncoord,incr) errmet=0; dland=1/incr; % La inversa dels increments ens servirà per trobar el % diferencial de càrrega Ax=zeros(N,1); % Dimensionament d'una matriu que emmagatzemarà les dàdes % dels axils dp=dland*Vector_esforcos_prescrits; % Diferencial de càrrega (utilitzant % la inversa dels increments) % Dimensiona a zero totes les dades que utlitzarem per optimitzar % l'execució K_Geom_Glob_cte=K_Geom_Glob_cte*0; % Matriu geomètrica global sum_dif_mov=Vector_moviments*0; % Moviments totals d_vector_moviments=Vector_moviments*0; % Moviments diferencials kriggeom_loc_ev=kriggeom_cte*0; % Matriu gemètrica local for j=1:incr; % Troba el diferencial de moviment d'aquest pas d_vector_moviments(1:nn,1)=(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte... (1:nn,1:nn))\dp(1:nn,1); while det(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte(1:nn,1:nn))<0 % Si el % determinant de la matriu geomètrica % global es superior al determinant % de la matriu de rigidesa global K_Geom_Glob_cte=[]; % Dimensionament de la matriu goemètrica global % buida errmet=1; % Valor pertinent per la variable d'error del mètode return % Retorna l'execució a la funció calculinestabilitat.m end sum_dif_mov=sum_dif_mov+d_vector_moviments; % Sumatori de moviments % A continuació e calcularàn els esforços en les barres per trobar % l'axil de casascuna d'elles i introduir-lo a les matrius geomètriques % locals % Calcula les reaccions als recolzaments de l'estructura. Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos; % Retorn de reaccions i moviments en l'ordre original dels nusos i % eixos globals. [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos... (sistema,Vector_reaccions,sum_dif_mov,n,Ncond); % Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems % de barra. K_Geom_Glob_cte=zeros(ns); % Dimensiona a zero la matriu geomètrica % global per evitar sumar l'actual a % l'anterior. for i=1:N n1=Lnods(i,1); % node inferior n2=Lnods(i,2); % node superior extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0 % rígid-articulat =1 % articulat-rígid =2 % articulat-articulat =3 vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra. L=norm(vector_barra); % Longitut de la barra. T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base % completa 6x6 % Càlcul de la submatriu de rigidesa en eixos locals. KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L); vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*... Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals. Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep... (i,:),barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),... barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems); Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*... vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra % en eixos locals. Ax(i,j+1)=Retorn_esforcos_barres(i,1); % Emmagatzema els axils de barra % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu % geomètrica global. Aquesta matriu geomètrica serà la que utilitzarem % al pas següent. kriggeom_cte_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte_tot(:,:,i); kriggeom_loc_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte(:,:,i); K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa... (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2); end end % Un cop s'han realitzat tots els increments trobarem la matriu % geomètrica global a partir dels desplaçaments finals dels nussos. % Calcula les reaccions als recolzaments de l'estructura. Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos; % Retorn de reaccions i moviments en l'ordre original dels nusos i % eixos globals. [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos(sistema,... Vector_reaccions,sum_dif_mov,n,Ncond); Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems % de barra. K_Geom_Glob_cte=zeros(ns); for i=1:N n1=Lnods(i,1); % node inferior n2=Lnods(i,2); % node superior extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0 % rígid-articulat =1 % articulat-rígid =2 % articulat-articulat =3 vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra. L=norm(vector_barra); % Longitut de la barra. T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base % completa 6x6 % Càlcul de la submatriu de rigidesa en eixos locals. KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L); vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*... Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals. Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),... barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),... barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems); Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*... vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra % en eixos locals. Ax_fin(i,1)=Retorn_esforcos_barres(i,1); % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu % geomètrica global. Aquesta matriu geomètrica serà la que la funció % retornarà. kriggeom_cte_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte_tot(:,:,i); kriggeom_loc_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte(:,:,i); K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa... (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2); end end
RESOLUCIÓ DEL PROBLEMA NO LINEAL MITJANÇANT EL MÈTODE DE RUNGE_KUTTA
% Funció que retorna la matriu geomètrica global calculada pel mètode % de Runge-Kutta function [K_Geom_Glob_cte,errmet]=Runge_Kutta(Vector_esforcos_prescrits,... K_Glob,TL2G,K_Geom_Glob_cte,sistema,Vector_moviments,... Vector_esforcos,n,nn,ns,Ncond,N,barra_carrega_rep,... barra_carrega_rep_proj,barra_carrega_punt,barra_carrega_mom,... kriggeom_cte_tot,kriggeom_cte,Lnods,Lartyp,Lmecp,Ncoord,incr,alpha1,... alpha2,mu) errmet=0; dland=1/incr; % La inversa dels increments ens servirà per trobar el % diferencial de càrrega Ax=zeros(N,1); % Dimensionament d'una matriu que emmagatzemarà les dàdes % dels axils Ax_prim=zeros(N,1); % Dimensionament d'una matriu que emmagatzemarà les % dàdes dels axils dels passos intermedis dp=dland*Vector_esforcos_prescrits; % Diferencial de càrrega (utilitzant % la inversa dels increments) % Dimensiona a zero totes les dades que utlitzarem per optimitzar % l'execució K_Geom_Glob_cte=K_Geom_Glob_cte*0; K_Geom_Glob_cte_prim=K_Geom_Glob_cte*0; sum_dif_mov=Vector_moviments*0; d_vector_moviments=Vector_moviments*0; d_vector_moviments_prim=Vector_moviments*0; kriggeom_loc_ev=kriggeom_cte*0; for j=1:incr; % Busca la matriu de rigidesa tangent a l'inici del pas a partir % dels moviments totals dels nussos al pas anteior. Si és el primer % pas, la matriu geomètrica serà igual a 0 i per tant la matriu de % rigidesa tangent serà igual a la matriu de rigidesa global. % Calcula les reaccions als recolzaments de l'estructura. Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos; % Retorn de reaccions i moviments en l'ordre original dels nusos i % eixos globals. [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos... (sistema,Vector_reaccions,sum_dif_mov,n,Ncond); Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems % de barra. K_Geom_Glob_cte=zeros(ns); for i=1:N n1=Lnods(i,1); % node inferior n2=Lnods(i,2); % node superior extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0 % rígid-articulat =1 % articulat-rígid =2 % articulat-articulat =3 vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra. L=norm(vector_barra); % Longitut de la barra. T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base % completa 6x6 % Càlcul de la submatriu de rigidesa en eixos locals. KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L); vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*... Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals. Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),... barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),... barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems); Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*... vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra % en eixos locals. Ax(i,j+1)=Retorn_esforcos_barres(i,1); % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu % geomètrica global. Aquesta matriu geomètrica serà la que el mètode % utilitzarà pel pas intermedi (K1). kriggeom_cte_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte_tot(:,:,i); kriggeom_loc_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte(:,:,i); K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa... (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2); end % Diferencial dels moviments al pas intermedi. Aquest moviments % s'hauràn de sumar al sumatori total de desplaçaments. d_vector_moviments_prim(1:nn,1)=(K_Glob(1:nn,1:nn)-... K_Geom_Glob_cte(1:nn,1:nn))\(mu*dp(1:nn,1)); sum_dif_mov_prim=sum_dif_mov+d_vector_moviments_prim; while det(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte(1:nn,1:nn))<0 % Si el % determinant de la matriu geomètrica % global es superior al determinant % de la matriu de rigidesa global K_Geom_Glob_cte=[]; % Dimensionament de la matriu goemètrica global % buida errmet=1; % Valor pertinent per la variable d'error del mètode return % Retorna l'execució a la funció calculinestabilitat.m end % A continuació e calcularàn els esforços en les barres per trobar % l'axil de casascuna d'elles i introduir-lo a les matrius geomètriques % locals % Calcula les reaccions als recolzaments de l'estructura. Vector_reaccions=K_Glob*sum_dif_mov_prim-Vector_esforcos; % Retorn de reaccions i moviments en l'ordre original dels nusos i % eixos globals. [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos(sistema,... Vector_reaccions,sum_dif_mov_prim,n,Ncond); Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems % de barra. K_Geom_Glob_cte_prim=zeros(ns); for i=1:N n1=Lnods(i,1); % node inferior n2=Lnods(i,2); % node superior extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0 % rígid-articulat =1 % articulat-rígid =2 % articulat-articulat =3 vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra. L=norm(vector_barra); % Longitut de la barra. T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base % completa 6x6 % Càlcul de la submatriu de rigidesa en eixos locals. KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L); vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*... Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals. Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),... barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),... barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems); Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*... vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra % en eixos locals. Ax_prim(i,j+1)=Retorn_esforcos_barres(i,1); % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu % geomètrica global. Aquesta matriu geomètrica (K2) serà la que el % mètode utilitzarà per convinar-la amb l'anterior (K1). kriggeom_cte_ev(:,:,i)=(Ax_prim(i,j+1)/(30*L))*kriggeom_cte_tot(:,:,i); kriggeom_loc_ev(:,:,i)=(Ax_prim(i,j+1)/(30*L))*kriggeom_cte(:,:,i); K_Geom_Glob_cte_prim=K_Geom_Glob_cte_prim+... ensamblatge_matriu_rigidesa(kriggeom_cte_ev(:,:,i),ns,sistema,... n1,n2); end % Diferencial dels moviments al final del pas. Aquest moviments % serviran per obtenir la matriu de rigidesa tangen a l'inici del pas % següent.. d_vector_moviments(1:nn,1)=(K_Glob(1:nn,1:nn)-(K_Geom_Glob_cte... (1:nn,1:nn)*alpha1+K_Geom_Glob_cte_prim(1:nn,1:nn)*alpha2))... \dp(1:nn,1); sum_dif_mov=sum_dif_mov+d_vector_moviments; while det(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte(1:nn,1:nn))<0 % Si el % determinant de la matriu geomètrica % global es superior al determinant % de la matriu de rigidesa global K_Geom_Glob_cte=[]; % Dimensionament de la matriu goemètrica global % buida errmet=1; % Valor pertinent per la variable d'error del mètode return % Retorna l'execució a la funció calculinestabilitat.m end end % Un cop s'han realitzat tots els increments trobarem la matriu % geomètrica global a partir dels desplaçaments finals dels nussos. % Calcula les reaccions als recolzaments de l'estructura. Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos; % Retorn de reaccions i moviments en l'ordre original dels nusos i % eixos globals. [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos(sistema,... Vector_reaccions,sum_dif_mov,n,Ncond); Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems % de barra. K_Geom_Glob_cte=zeros(ns); for i=1:N n1=Lnods(i,1); % node inferior n2=Lnods(i,2); % node superior extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0 % rígid-articulat =1 % articulat-rígid =2 % articulat-articulat =3 vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra. L=norm(vector_barra); % Longitut de la barra. T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base % completa 6x6 % Càlcul de la submatriu de rigidesa en eixos locals. KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L); vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*... Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals. Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),... barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),... barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems); Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*... vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra % en eixos locals. Ax_fin(i,1)=Retorn_esforcos_barres(i,1); % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu % geomètrica global. Aquesta matriu geomètrica serà la que la funció % retornarà. kriggeom_cte_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte_tot(:,:,i); kriggeom_loc_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte(:,:,i); K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa... (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2); end end
FUNCIÓ PER REORDENAR ELS MODES DE VINCLAMENT
% Funció que ordena els vectors del mode de vinclament perquè aquest % tinguin l'ordre dels nussos function v_reord=reordenacio_vec(v,sistema,n,Ncond) v_reord=zeros(n*6,1); for i=1:n for j=1:6 k=sistema(6*i+j-6); if k v_reord(6*i+j-6)=v(k); end end if Ncond(i,1)==1 % Si és rodet inclinat torna desplaçaments del sistema % d'eixos locals del rodet al global Tr=calcul_canvi_base(Ncond(i,2:4)); % Calcula matriu canvi de base % del rodet. T=[Tr,zeros(3);zeros(3),Tr]; v_reord(6*i-5:6*i)=T*v_reord(6*i-5:6*i); end end end