Estructura bsica del programa de elementos finitos Coordenadas

Estructura básica del programa de elementos finitos Coordenadas. Elemento 1=[1, 1; 0, 0; 1. 5, 0]; Datos: Coordenadas de cada nodo Condiciones de contorno Relación entre las C locales y globales Datos de la Matriz constitutiva Coordenadas. Elemento 2=[2. 5, 1; 1. 5, 0; 2. 5, 0]; Nodos. Datos=[1, 2, 5, 6]; Valor. Nodos. Datos=[273, 278, 278]; Vector. Fuerzas=zeros(6, 1); NODES=[3, 5; 1, 3; 2, 4; 4, 6]; datos=[400, 0]; %Calc. K(XY, tipodeproblema, tipodeelemento, Calculo de las matrices elementales puntosdecuadratura, datos) KE 1=Calc. K(Coordenadas. Elemento 1, 1, 3, 2, datos); KE 2=Calc. K(Coordenadas. Elemento 2, 1, 3, 2, datos); KG=zeros(6, 6); Ensamble de las matrices elementales en la matriz global KG=ensamble(KG, KE 1, NODES, 1); KG=ensamble(KG, KE 2, NODES, 2); Solucion=reducir(KG, Vector. Fuerzas, Nodos. Datos, Incorporación de las CC y solución del sistema de ecuaciones Valor. Nodos. Datos);

Calc. K 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. 861136311594053]; P=[0. 347854845137454, 0. 652145154862546, 0. 347854845137454]; end

Calc. K 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_2 D_4 N(r, s); %Derivadas de la forma: % Dh 1/Dr Dh 2/Dr. . Dh 4/Dr % Dh 1/Ds Dh 2/Ds. . Dh 4/Ds Je=DH*Coordenadas. Elemento; %genera el Jacobiano Inv. Je=Je^-1; %Calcula la inversa det. Je=det(Je); %Calcula el determinante del jacobiano Be=Inv. Je*DH; % Be es la matriz de las derivadas de las funciones de forma, de la forma: % Dh 1/Dx Dh 2/Dx. . Dh 4/Dx % Dh 1/Dy Dh 2/Dy. . Dh 4 s/Dy K=K+p*Be'*Ce*Be*det. Je; end

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;

Reducir function Solucion=reducir(KG, Vector. Fuerzas, Nodos. Datos, Valor. Nodos. Datos) Ncondicionesdecontorno=length(Nodos. Datos); Nnodosglobales=length(KG); for ig=1: Ncondicionesdecontorno nodo=Nodos. Datos(ig); for jg=1: Nnodosglobales if jg~=nodo Vector. Fuerzas(jg)=Vector. Fuerzas(jg)-KG(jl, nodo)*Valor. Nodos. Datos(ig); KG(jg, nodo)=0; end for jg=1: Nnodosglobales if jg~=nodo KG(nodo, jg)=0; else KG(nodo, nodo)=1; Vector. Fuerzas(nodo)=Valor. Nodos. Datos(ig); end end %[V, D] = eig(KG); %V %D Solucion=KGVector. Fuerzas;

Problema de elasticidad Coordenadas. Elemento 1=[ 10, 10; 0, 0; 10, 0; 5, 10; 0, 5; 5, 0; Datos: Coordenadas de cada nodo Condiciones de contorno Relación entre las C locales y globales Datos de la Matriz constitutiva 10, 5]; Nodos. Datos=[6, 7, 8, 13, 14]; Valor. Nodos. Datos=zeros(1, 5); NODES=[1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16]; Vector. Fuerzas=zeros(16, 1); Vector. Fuerzas(2, 1)=-100; Vector. Fuerzas(4, 1)=100; Vector. Fuerzas(10, 1)=100; E=2. 1 E 6; v=0. 3; datos=[E, v]; %Calc. K(XY, tipodeproblema, tipodeelemento, pu ntosdecuadratura, datos) Calculo de las matrices elementales KE 1=Calc. K(Coordenadas. Elemento 1, 2, 4, 3, datos); KG=zeros(16, 16); Ensamble de las matrices elementales en la matriz global KG=ensamble(KG, KE 1, NODES, 1); Solucion=reducir(KG, Vector. Fuerzas, Nodos. Datos, Valor. Nodos. Datos); Incorporación de las CC y solución del sistema de ecuaciones

Calc. K elasticidad [E, v]=datos; Ce 2 D=Calc_Ce 2 D(E, v); K=zeros(16, 16); Be 2 D=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_2 D_8 N(r, s); NNodos=8; Je=DH*Coordenadas. Elemento; det. Je=det(Je); Inv. Je=Je^-1; Be=Inv. Je*DH; % Be es la matriz de las derivadas de las funciones de forma, de la forma: % Dh 1/Dx Dh 2/Dx. . Dh 8/Dx % Dh 1/Dy Dh 2/Dy. . Dh 8/Dy for m=1: NNodos Be 2 D(1, (m*2 -1))=Be(1, m); Be 2 D(2, (m*2))=Be(2, m); Be 2 D(3, (m*2 -1))=Be(2, m)/2; Be 2 D(3, (m*2))=Be(1, m)/2; end % Con eso creo una matriz de 3 x 2 NNodos, de la forma % % Dh 1/Dx 0 0 Dh 2/Dx Dh 1/Dy 0 0 Dh 2/Dy . . . . % 1/2*Dh 1/Dx 1/2*Dh 1/Dy 1/2*Dh 2/Dx 1/2*Dh 2/Dy K=K+p*Be 2 D'*Ce 2 D*Be 2 D*det. Je; end End . .

Recuperar deformaciones y tensiones Valores. Nodales. Elemento=zeros(16, 1); for jl=1: 16 jg=NODES(jl, 1); Valores. Nodales. Elemento(jl)=Solucion(jg); end Vector. Tensiones=zeros(21, 3); Vector. Deformaciones=zeros(21, 3); for m=1: 21 r=1; s=(m-1)/10 -1; Be=Calc. B(Coordenadas. Elemento 1, 2, 4, datos, r, s); Ce 2 D=Calc_Ce 2 D(E, v); deformaciones=Be*Valores. Nodales. Elemento; tensiones=Ce 2 D*deformaciones; Vector. Deformaciones(m, 1: 3)=deformaciones'; Vector. Tensiones(m, 1: 3)=tensiones'; end save Vector. Deformaciones. dat Vector. Deformaciones -Ascii; save Vector. Tensiones. dat Vector. Tensiones -Ascii;