1 CoordenadasElemento1=[1,1;0,1;0,0;1.5,0]; CoordenadasElemento2=[2.5,1;1,1;1.5,0;2.5,0]; NodosDatos=[1,2,5,6]; ValorNodosDatos=[273,273,278,278]; VectorFuerzas=zeros(6,1); NODES=[3,5;1,3;2,4;4,6]; datos=[400,0]; %CalcK(XY,tipodeproblema,tipodeelemento, puntosdecuadratura,datos) KE1=CalcK(CoordenadasElemento1,1,3,2,datos); KE2=CalcK(CoordenadasElemento2,1,3,2,datos); KG=zeros(6,6); KG=ensamble(KG,KE1,NODES,1); KG=ensamble(KG,KE2,NODES,2); Solucion=reducir(KG,VectorFuerzas,NodosDatos, ValorNodosDatos); Datos: Coordenadas de cada nodo Condiciones de contorno Relación entre las C locales y globales Datos de la Matriz constitutiva Calculo de las matrices elementales Ensamble de las matrices elementales en la matriz global Incorporación de las CC y solución del sistema de ecuaciones Estructura básica del programa de elementos finitos

2 CalcK if puntosdecuadratura==1 R=[0]; P=[2]; elseif puntosdecuadratura==2 R=[-sqrt(1/3),sqrt(1/3)]; P=[1,1]; elseif puntosdecuadratura==3 R=[-sqrt(3/5),0,sqrt(3/5)]; P=[5./9.,8./9.,5./9.]; elseif puntosdecuadratura==4 R=[-0.861136311594053,-0.339981043584856,0.339981043584856,0.861136311594053]; P=[0.347854845137454,0.652145154862546,0.652145154862546,0.347854845137454]; end

3 CalcK 2 K=zeros(4,4); k=datos(1); Ce=k*eye(2); for i=1:puntosdecuadratura for j=1:puntosdecuadratura r=R(j); s=R(i); p=P(i)*P(j); DH=Calc_DH_2D_4N(r,s); %Derivadas de la forma: % Dh1/Dr Dh2/Dr.... Dh4/Dr % Dh1/Ds Dh2/Ds.... Dh4/Ds Je=DH*CoordenadasElemento; %genera el Jacobiano InvJe=Je^-1; %Calcula la inversa detJe=det(Je); %Calcula el determinante del jacobiano Be=InvJe*DH; % Be es la matriz de las derivadas de las funciones de forma, de la forma: % Dh1/Dx Dh2/Dx.... Dh4/Dx % Dh1/Dy Dh2/Dy.... Dh4s/Dy K=K+p*Be'*Ce*Be*detJe; end

4 Ensamble function a=ensamble(KG,KE,NODES,elemento) Nodosporelemento=length(KE); for il=1:Nodosporelemento for jl=1:Nodosporelemento ig=NODES(il,elemento); jg=NODES(jl,elemento); KG(ig,jg)=KG(ig,jg)+K(il,jl); end a=KG;

5 Reducir function Solucion=reducir(KG,VectorFuerzas,NodosDatos,ValorNodosDatos) Ncondicionesdecontorno=length(NodosDatos); Nnodosglobales=length(KG); for ig=1:Ncondicionesdecontorno nodo=NodosDatos(ig); for jg=1:Nnodosglobales if jg~=nodo VectorFuerzas(jg)=VectorFuerzas(jg)-KG(jl,nodo)*ValorNodosDatos(ig); KG(jg,nodo)=0; end for jg=1:Nnodosglobales if jg~=nodo KG(nodo,jg)=0; else KG(nodo,nodo)=1; VectorFuerzas(nodo)=ValorNodosDatos(ig); end %[V,D] = eig(KG); %V %D Solucion=KG\VectorFuerzas;

6 Problema de elasticidad CoordenadasElemento1=[ 10,10; 0,10; 0,0; 10,0; 5,10; 0,5; 5,0; 10,5]; NodosDatos=[6,7,8,13,14]; ValorNodosDatos=zeros(1,5); NODES=[1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16]; VectorFuerzas=zeros(16,1); VectorFuerzas(2,1)=-100; VectorFuerzas(4,1)=100; VectorFuerzas(10,1)=100; E=2.1E6; v=0.3; datos=[E,v]; %CalcK(XY,tipodeproblema,tipodeelemento,pu ntosdecuadratura,datos) KE1=CalcK(CoordenadasElemento1,2,4,3,datos); KG=zeros(16,16); KG=ensamble(KG,KE1,NODES,1); Solucion=reducir(KG,VectorFuerzas,NodosDatos,V alorNodosDatos); Datos: Coordenadas de cada nodo Condiciones de contorno Relación entre las C locales y globales Datos de la Matriz constitutiva Calculo de las matrices elementales Ensamble de las matrices elementales en la matriz global Incorporación de las CC y solución del sistema de ecuaciones

7 CalcK elasticidad [E,v]=datos; Ce2D=Calc_Ce2D(E,v); K=zeros(16,16); Be2D=zeros(3,16); for i=1:puntosdecuadratura for j=1:puntosdecuadratura r=R(j); s=R(i); p=P(i)*P(j); DH=Calc_DH_2D_8N(r,s); NNodos=8; Je=DH*CoordenadasElemento; detJe=det(Je); InvJe=Je^-1; Be=InvJe*DH; % Be es la matriz de las derivadas de las funciones de forma, de la forma: % Dh1/Dx Dh2/Dx.... Dh8/Dx % Dh1/Dy Dh2/Dy.... Dh8/Dy for m=1:NNodos Be2D(1,(m*2-1))=Be(1,m); Be2D(2,(m*2))=Be(2,m); Be2D(3,(m*2-1))=Be(2,m)/2; Be2D(3,(m*2))=Be(1,m)/2; end % Con eso creo una matriz de 3x2NNodos, de la forma % Dh1/Dx 0 Dh2/Dx 0.... % 0 Dh1/Dy 0 Dh2/Dy.... % 1/2*Dh1/Dx 1/2*Dh1/Dy 1/2*Dh2/Dx 1/2*Dh2/Dy.... K=K+p*Be2D'*Ce2D*Be2D*detJe; end End

8 Recuperar deformaciones y tensiones ValoresNodalesElemento=zeros(16,1); for jl=1:16 jg=NODES(jl,1); ValoresNodalesElemento(jl)=Solucion(jg); end VectorTensiones=zeros(21,3); VectorDeformaciones=zeros(21,3); for m=1:21 r=1; s=(m-1)/10-1; Be=CalcB(CoordenadasElemento1,2,4,datos,r,s); Ce2D=Calc_Ce2D(E,v); deformaciones=Be*ValoresNodalesElemento; tensiones=Ce2D*deformaciones; VectorDeformaciones(m,1:3)=deformaciones'; VectorTensiones(m,1:3)=tensiones'; end save VectorDeformaciones.dat VectorDeformaciones -Ascii; save VectorTensiones.dat VectorTensiones -Ascii;