metodo FEM para 1D

Vista previa del material en texto

Método de elementos finitos para 1D
Por: Daniela Alejandra Benavides Tapia
Elementos finitos (2021-I)
Contenido
Introducción.............................................................................................................................................................. 1
Sistema local............................................................................................................................................................ 2
Sistema global.......................................................................................................................................................... 2
Establezco condiciones de contorno (Dirichlet)....................................................................................................3
Ecuacion diferencial a resolver:
 
donde: , , , , , y además: y , y 
tiene por dominio: 
Introducción
clc, clear all,close all
syms xi
L=4;
delta=7840;
g=9.81;
u=0;
A=1;
E=200*(10^9);
k=E*A;
f=(A*delta*g)*((L-xi)/L);
N=5; %número de elementos
xmax=L; 
xmin=0;
T1=0; %condición de contorno 1
dt=-1; %condición de contorno 2
Dt=[0;0;0;0;dt]; %vector condición de contorno
Definir funciones base y derivadas
S=[(1/2)*(1-xi);(1/2)*(1+xi)];
1
dS=[diff(S(1));diff(S(2))];
A continuacion defino el paso del sistema a partir del dominio y numero de elementos
h=(xmax-xmin)/N;
Magnitud del Jacobiano
J=h/2;
Matriz de conversion de coordenadas
GtoL=[1 2; 2 3; 3 4; 4 5; 5 6];
Sistema local
Se define la matriz de coordenadas locales kij y fi 
for i=1:2
 for j=1:2
 kij(i,j)=int((S(i)*u*dS(j)/J+k*dS(i)*(1/J)*dS(j)*(1/J))*J,-1,1); 
 end 
 fi(i,1)=int(S(i)*f*J,-1,1);
end
disp('Matriz elemental de rigidez')
Matriz elemental de rigidez
kij
kij = 
disp('Vector fuerza elemental')
Vector fuerza elemental
fi
fi = 
Sistema global
Definiendo matriz de ceros, la cual tiene el tamaño: N+1xN+1 (numero de nodos)
Kij=zeros(N+1,N+1);
Fi=zeros(N+1,1);
Rellenando las matrices globales
2
for e=1:N
 for i=1:2
 for j=1:2
 I=GtoL(e,i);
 J=GtoL(e,j);
 Kij(I,J)=Kij(I,J)+kij(i,j);
 end
 I=GtoL(e,i);
 Fi(I,1)=Fi(I)+fi(i,1);
 end
end
disp('Matriz global de rigidez')
Matriz global de rigidez
Kij
Kij = 6×6
1011 ×
 2.5000 -2.5000 0 0 0 0
 -2.5000 5.0000 -2.5000 0 0 0
 0 -2.5000 5.0000 -2.5000 0 0
 0 0 -2.5000 5.0000 -2.5000 0
 0 0 0 -2.5000 5.0000 -2.5000
 0 0 0 0 -2.5000 2.5000
disp('Vector fuerza global')
Vector fuerza global
Fi
Fi = 6×1
104 ×
 3.3328
 6.1528
 6.1528
 6.1528
 6.1528
 2.8200
Establezco condiciones de contorno (Dirichlet)
Eliminando filas y columnas 1 porque vale cero, y añadiendo el valor dado para la primera derivada dada como 
condición de contorno
for i=2:6
 Fi(i,1)=Fi(i,1)-Kij(i,1)*T1;
end 
Eliminando filas y columnas 1
Fi(1,:)=[];
Kij(1,:)=[];
Kij(:,1)=[];
Kij
3
Kij = 5×5
1011 ×
 5.0000 -2.5000 0 0 0
 -2.5000 5.0000 -2.5000 0 0
 0 -2.5000 5.0000 -2.5000 0
 0 0 -2.5000 5.0000 -2.5000
 0 0 0 -2.5000 2.5000
Añadiendo la condición de contorno dad para la primera derivada
B=Fi+Dt; 
Fi=B;
Finalmente, solucionando
T=Kij\Fi
T = 5×1
10-5 ×
 0.1097
 0.1948
 0.2553
 0.2912
 0.3025
Agregando T1 a la matriz
Tij=zeros(6,1);
for i=1:5
 Tij(i+1,1)=T(i,1);
end
x_plot=[0:h:L];
Graficando
plot(x_plot,Tij,'k-s')
title('Deformación en la barra')
legend('FEM 1D','location',"south")
xlabel('x')
ylabel('T')
4
disp('Deformacion en nodo 1 [mm]')
Deformacion en nodo 1 [mm]
Nd1=Tij(1)*1000;
disp(Nd1)
 0
disp('Deformacion en nodo 2 [mm]')
Deformacion en nodo 2 [mm]
Nd2=Tij(2)*1000;
disp(Nd2)
 0.0011
disp('Deformacion en nodo 3 [mm]')
Deformacion en nodo 3 [mm]
Nd3=Tij(3)*1000;
disp(Nd3)
 0.0019
5