UNIVERSIDAD NACIONAL DE
INGENIERÍA
FACULTAD DE INGENIERÍA MECÁNICA
Especialidad de Ingeniería Mecánica-Eléctrica
CALCULO POR ELEMENTOS FINITOS
(MC516 - F)
Práctica N°1: Tracción
PRESENTADO POR:
Apellidos y nombres Código
Mayta Aliaga Diego José 20184029A
PROFESOR
Cueva Pacheco, Ronald
LIMA, 202
Página | 0
INDICE
1. DESARROLLO DEL PROBLEMA ........................................................................2
1.1 Enunciado del problema .........................................................................................2
1.2 Modelado del cuerpo................................................................................................3
1.3 Grados de libertad nodales .......................................................................................4
1.4 Vector carga .............................................................................................................5
1.5 Matriz de Rigidez .....................................................................................................6
1.6 Ecuación de Rigidez y condiciones de contorno .....................................................7
1.7 Cálculo de esfuerzos ................................................................................................8
1.8 Resultados ................................................................................................................9
2. DIAGRAMA DE FLUJO DEL PROGRAMA ......................................................10
3. DIGITACIÓN Y EJECUCIÓN (MATLAB) .........................................................13
3.1 Código de Matlab ...................................................................................................13
3.2 Vista en el “Command Window” ..........................................................................15
4. CONCLUSIONES ....................................................................................................16
Página | 1
1. DESARROLLO DEL PROBLEMA
1.1 Enunciado del problema
En la figura se muestra una placa triangular simétrica; con espesor constante igual a 120
mm; está sometida a su peso propio y a una fuerza concentrada, tal como se indica.
DATOS:
𝑔𝑟 − 𝑓
𝛾 = 8.0
𝑐𝑚3
𝑁
𝐸 = 2.00 ∙ 105
𝑚𝑚2
HALLAR:
• La distribución de esfuerzos a lo largo de la barra; modelándola con cuatro (04)
elementos finitos (unidimensionales).
• La fuerza de reacción en el apoyo.
FIGURA:
𝟕𝟓𝟎 𝐦𝐦
𝒍 = 𝟏𝟓𝟎𝟎 𝐦𝐦
𝑷 = 𝟑𝟎𝟎𝟎𝟎 𝑵
𝒉 = 𝟏𝟎𝟎𝟎 𝐦𝐦
Página | 2
1.2 Modelado del cuerpo
Se considerarán cuatro elementos finitos, todos con la misma longitud de 375 mm.
En cuanto a los espesores, estos se calculan como la semisuma de las bases de los
trapecios que se forman al hacer cada partición:
250 + 0
𝑏1 = = 125 𝑚𝑚
2
500 + 250
𝑏2 = = 375 𝑚𝑚
2
750 + 500
𝑏3 = = 625 𝑚𝑚
2
1000 + 750
𝑏4 = = 875 𝑚𝑚
2
Entonces, el modelado del cuerpo quedaría así:
(X)
(1)
𝟑𝟕𝟓 𝐦𝐦
𝒃𝟏 = 𝟏𝟐𝟓 𝐦𝐦
(2)
𝟑𝟕𝟓 𝐦𝐦
𝒃𝟐 = 𝟑𝟕𝟓 𝐦𝐦
(3)
𝑷 = 𝟑𝟎𝟎𝟎𝟎 𝑵
𝟑𝟕𝟓 𝐦𝐦
𝒃𝟑 = 𝟔𝟐𝟓 𝐦𝐦
(4)
𝟑𝟕𝟓 𝐦𝐦
𝒃𝟒 = 𝟖𝟕𝟓 𝐦𝐦
(5)
𝑹
Página | 3
Además, las áreas de cada elemento finito se calculan de la siguiente manera:
𝑨𝒊 = 𝒃𝒊 ∙ 𝒕
Tabla de Conectividad
Grados de Libertad
Nodos
(GDL) le Ae
e (mm) (mm2)
(1) (2) 1 2
1 1 2 Q1 Q2 375 15000
2 2 3 Q2 Q3 375 45000
3 3 4 Q3 Q4 375 75000
4 4 5 Q4 Q5 375 105000
1.3 Grados de libertad nodales (Vector Desplazamiento)
El vector desplazamiento será:
(Q1)
𝑄1
𝑄2
𝑄𝑗 = 𝑄3
(Q2) 𝑄4
[𝑄5 ]
Donde el valor numérico de Q5 es 0 debido
(Q3) al empotramiento de la placa. Los valores de
𝑄1 , 𝑄2 , 𝑄3 , 𝑄4 serán calculados.
𝑄1
𝑄2
(Q4)
𝑄𝑗 = 𝑄3
𝑄4
[0]
(Q5)
Página | 4
1.4 Vector Carga
(1)
F1
Debido a que la densidad es:
(2)
𝒈𝒓 − 𝒇 𝑵
F2 𝜸 = 𝟖. 𝟎 𝟑
= 𝟕𝟖. 𝟒𝟖 ∙ 𝟏𝟎−𝟔
𝒄𝒎 𝒎𝒎𝟑
(3) Se hallará el peso con este valor asumiendo que
este se distribuye de manera simétrica en cada
F3 nodo.
𝒑 ∙ 𝑳𝒆 (𝑨𝑳)𝒆
(4) 𝑭𝒆𝒔 = +𝒇∙ + 𝑷𝑺
𝟐 𝟐
F4
(5)
F5
Analizamos las fuerzas en cada elemento finito:
Elemento Finito 1:
(𝐴𝐿)1 375
𝐹11 = 𝑓 ∙ + 𝑃1 = −78.48 ∙ 10−6 ∙ 15000 ∙ + 0 = −220.725 𝑁
2 2
(𝐴𝐿)1 375
𝐹21 =𝑓∙ + 𝑃2 = −78.48 ∙ 10−6 ∙ 15000 ∙ + 0 = −220.725 𝑁
2 2
Elemento Finito 2:
(𝐴𝐿)2 375
𝐹22 =𝑓∙ + 𝑃2 = −78.48 ∙ 10−6 ∙ 45000 ∙ + 0 = −662.175 𝑁
2 2
(𝐴𝐿)2 375
𝐹32 = 𝑓 ∙ + 𝑃3 = −78.48 ∙ 10−6 ∙ 45000 ∙ + 0 = −662.175 𝑁
2 2
Elemento Finito 3:
(𝐴𝐿)3 375
𝐹33 = 𝑓 ∙ + 𝑃3 = −78.48 ∙ 10−6 ∙ 75000 ∙ − 30000 = −31103.625 𝑁
2 2
(𝐴𝐿)3 375
𝐹43 =𝑓∙ + 𝑃4 = −78.48 ∙ 10−6 ∙ 75000 ∙ + 0 = −1103.625 𝑁
2 2
Página | 5
Elemento Finito 4:
(𝐴𝐿)4 375
𝐹44 = 𝑓 ∙ + 𝑃4 = −78.45 ∙ 10−6 ∙ 105000 ∙ + 0 = −1545.075 𝑁
2 2
(𝐴𝐿)4 375
𝐹54 = 𝑓 ∙ + 𝑃5 = −78.45 ∙ 10−6 ∙ 105000 ∙ + 𝑅 = (−1545.075 + 𝑅) 𝑁
2 2
Ahora analizamos las fuerzas en cada elemento finito:
𝐹1 = 𝐹11 = −220.725 𝑁
𝐹2 = 𝐹21 + 𝐹22 = −882.9 𝑁
𝐹3 = 𝐹32 + 𝐹33 = −31765.8 𝑁
𝐹4 = 𝐹43 + 𝐹44 = −2648.7 𝑁
𝐹5 = 𝐹54 = (−1545.075 + 𝑅) 𝑁
Entonces el vector carga queda de la siguiente manera:
𝐹1 −220.725
𝐹2 −882.9
𝐹 = 𝐹3 = −31765.8 𝑁
𝐹4 −2648.7
[𝐹5 ] [−1545.075 + 𝑅 ]
1.5 Matriz de Rigidez
A continuación, calcularemos la matriz de rigidez global:
1 −1 0 0 0 0 0 0 0 0
EA 1 −1 1 0 0 0 EA 2 0 1 −1 0 0
K ij = ( ) ∙ 0 0 0 0 0 + ( ) ∙ 0 −1 1 0 0
l l
0 0 0 0 0 0 0 0 0 0
(0 0 0 0 0) (0 0 0 0 0)
0 0 0 0 0 0 0 0 0 0
EA 3 0 0 0 0 0 EA 4 0 0 0 0 0
+( ) ∙ 0 0 1 −1 0 + ( ) ∙ 0 0 0 0 0
l l
0 0 −1 1 0 0 0 0 1 −1
(0 0 0 0 0) (0 0 0 −1 1 )
Página | 6
Considerando que E = 2*105 N/mm2 para todos los elementos finitos, se
reemplazan los datos en la ecuación y se obtiene:
8000000 −8000000 0 0 0
−8000000 32000000 −24000000 0 0 𝑁
K ij = 0 −24000000 64000000 −40000000 0
𝑚𝑚
0 0 −40000000 96000000 −56000000
( 0 0 0 −56000000 56000000 )
1.6 Ecuación de Rigidez y condiciones de contorno
La ecuación de rigidez:
𝑭𝒊 = 𝑲𝒊𝒋 ∙ 𝑸𝒋
Reemplazando los valores calculados:
−220.725 8 −8 0 0 0 𝑄1
−882.9 −8 32 −24 0 0 𝑄2
−31765.8 = 106 ∙ 0 −24 64 −40 0 × 𝑄3
−2648.7 0 0 −40 96 −56 𝑄4
[−1545.075 + 𝑅 ] (0 0 0 −56 56 ) [ 0]
Para obtener los desplazamientos:
−220.725 8 −8 0 0 𝑄1
−8 32 −24 0 𝑄
[ −882.9 ] = 106 ∙ ( ) × [ 2]
−31765.8 0 −24 64 −40 𝑄3
−2648.7 0 0 −40 96 𝑄4
Resolviendo el sistema de ecuaciones se tiene:
𝑸𝟏 = −𝟏𝟓𝟐. 𝟗𝟓𝟔𝟑 × 𝟏𝟎−𝟓 𝒎𝒎
𝑸𝟐 = −𝟏𝟓𝟎. 𝟏𝟗𝟕𝟐 × 𝟏𝟎−𝟓 𝒎𝒎
𝑸𝟑 = −𝟏𝟒𝟓. 𝟓𝟗𝟖𝟖 × 𝟏𝟎−𝟓 𝒎𝒎
𝑸𝟒 = −𝟔𝟑. 𝟒𝟐𝟓𝟐 × 𝟏𝟎−𝟓 𝒎𝒎
Página | 7
Luego para obtener la reacción en el empotramiento:
𝑄1
𝑄2
[−1545.075 + 𝑅] = 106 ∙ [0 0 0 −56 56 ] × 𝑄3
𝑄4
[0]
Resolviendo obtenemos:
𝑹 = 𝟑𝟕𝟎𝟔𝟑. 𝟐 𝑵
El valor calculado mediante las condiciones de equilibrio:
1000 ∙ 1500
𝑅 = 𝑃 + 𝑊 = 30000 + ( ) (120)(78.48 ∙ 10−6 )
2
𝑅 = 37063.2 𝑁
1.7 Cálculo de esfuerzos
Para calcular los valores de los esfuerzos por elemento finito, aplicamos la siguiente
relación:
𝑬 𝒆 𝑸
𝒆
𝝈 = ( ) [−𝟏 𝟏] [𝑸 𝒊 ]
𝒍 𝒊+𝟏
Se obtiene:
2 × 105 −152.9563 𝑵
𝝈𝟏 = ( ) [−1 1] [ ] × 10−5 = 𝟎. 𝟎𝟏𝟒𝟕𝟏𝟓
375 −150.1972 𝒎𝒎𝟐
2 × 105 −150.1972 𝑵
𝝈𝟐 = ( ) [−1 1] [ ] × 10−5 = 𝟎. 𝟎𝟐𝟒𝟓𝟐𝟓
375 −145.5988 𝒎𝒎𝟐
𝟑
2 × 105 −145.5988 𝑵
𝝈 =( ) [−1 1] [ ] × 10−5 = 𝟎. 𝟒𝟑𝟖𝟐𝟓𝟗
375 −63.4252 𝒎𝒎𝟐
2 × 105 −63.4252 𝑵
𝝈𝟒 = ( ) [−1 1] [ ] × 10−5 = 𝟎. 𝟑𝟑𝟖𝟐𝟔𝟖
375 0 𝒎𝒎𝟐
Página | 8
1.8 Resultados
𝑵
𝝈𝟏 = 𝟎. 𝟎𝟏𝟒𝟕𝟏𝟓
𝒎𝒎𝟐
𝑵
𝝈𝟐 = 𝟎. 𝟎𝟐𝟒𝟓𝟐𝟓
𝒎𝒎𝟐
a.
𝑵
𝝈𝟑 = 𝟎. 𝟒𝟑𝟖𝟐𝟓𝟗
𝒎𝒎𝟐
𝑵
𝝈𝟒 = 𝟎. 𝟑𝟑𝟖𝟐𝟔𝟖
𝒎𝒎𝟐
b. 𝑹 = 𝟑𝟕𝟎𝟔𝟑. 𝟐 𝑵
Página | 9
2. DIAGRAMA DE FLUJO DEL PROGRAMA
Página | 10
Página | 11
Página | 12
3. DIGITACIÓN Y EJECUCIÓN (MATLAB)
3.1 Código de Matlab
clc
disp('................................................................
.......')
disp(' INGRESO DE DATOS')
disp('----------------------------------------------------------------
-------')
h=input('Ingrese la altura de la placa triangular en mm: ');
b=input('Ingrese la base de la placa triangular en mm: ');
t=input('Ingrese el espesor de la placa en mm: ');
P=input('Ingrese el valor de P en N: ');
nodo=input('Nodo en el que se encuentra la fuerza P: ');
E=input('Ingrese el valor del Módulo de elasticidad en N/mm2: ');
vector_divisiones=input('Ingrese las divisiones en mm (en un vector
fila): ');
densidad=input('Ingrese el valor de la densidad en gr - f/cm3: ');
% DEFINIMOS LAS VARIABLES QUE UTILIZAREMOS MÁS ADELANTE
densidad=densidad*9.81*10^(-6);
n=length(vector_divisiones);
DIV=vector_divisiones';
bases_medias=zeros(n,1);
Areas=zeros(n,1);
bases=zeros(n,1);
fuerzas=zeros(2,n);
vector_carga=zeros(n+1,1);
A=zeros(n+1);
K=zeros(n+1);
esfuerzos=zeros(n,1);
%HALLAMOS EL ANCHO Y EL ÁREA TRANSVERSAL DE CADA ELEMENTO FINITO
for i=1:n;
bases(i)=b*sum(DIV(1:i))/h;
end
bases=[0; bases];
for i=1:n;
bases_medias(i)=(bases(i)+bases(i+1))/2;
end
for i=1:n;
Areas(i)=t*bases_medias(i);
end
Areas; %Vector que contiene las áreas transversales de cada elemento
finito
%HALLAMOS EL VECTOR CARGA
for i=1:n;
fuerzas(1,i)=(-densidad*Areas(i)*DIV(i)/2);
fuerzas(2,i)=(-densidad*Areas(i)*DIV(i)/2);
end
fuerzas=[fuerzas(1,:) 0; 0 fuerzas(2,:)];
for i=1:n+1;
vector_carga(i)=sum(fuerzas(:,i));
end
vector_carga(nodo)=vector_carga(nodo)-P; %Vector Carga (sin contar R)
Página | 13
%HALLAMOS LA MATRIZ DE RIGIDEZ
for i=1:n;
A(i,i)=1;
A(i+1,i+1)=1;
A(i,i+1)=-1;
A(i+1,i)=-1;
K=K+(E*Areas(i)/DIV(i))*A;
A=zeros(n+1);
end
K; %Matriz de Rigidez
%OBTENEMOS EL VECTOR DE DESPLAZAMIENTOS
D=K(1:n,1:n)\vector_carga(1:n);
vector_desplazamientos=[D;0]; %Vector de desplazamientos
%HALLAMOS LA REACCIÓN R EN EL APOYO:
R=K(n+1,:)*[D;0]-vector_carga(n+1);
R; %Reacción en el apoyo
%HALLAMOS EL VECTOR DE ESFUERZOS:
for i=1:n;
esfuerzos(i)=(E/DIV(i))*[-1 1]*vector_desplazamientos(i:i+1);
end
esfuerzos; %Vector de esfuerzos
%RESULTADOS
disp('................................................................
.......')
disp(' RESULTADOS')
disp('----------------------------------------------------------------
-------')
disp('El vector de desplazamientos en mm:')
disp(vector_desplazamientos)
disp('El vector de esfuerzos en MPa:')
disp(esfuerzos)
disp('La reacción R en el apoyo en N: ')
disp(R)
Página | 14
3.2 Vista en el “Command Window”
Página | 15
4. CONCLUSIONES
• El Método de los Elementos finitos nos permite analizar cuerpos y sistemas complejos
a través de su seccionamiento en partes simplificadas denominadas “Elementos
finitos”. Así, obtendremos valores aproximados para las variables analizadas.
• El Método de los Elementos Finitos solo nos proporciona resultados concretos en los
nodos del cuerpo. Por lo tanto, si se desea conocer el valor de una variable en otra
parte del cuerpo, se debe interpolar con los valores conocidos en los nodos.
• Notamos que se obtuvieron valores negativos en la matriz de desplazamientos; esto
tiene sentido ya que se intuye que el cuerpo se comprime y, por lo tanto, cada nodo
se desplaza hacia abajo. Justamente, en la dirección contraria al eje positivo X que se
definió en un inicio.
• A partir de los resultados, se observa claramente como los esfuerzos soportados por
los elementos finitos 3 y 4 son notablemente mayores a los esfuerzos hallados en el
resto del cuerpo. Esto tiene sentido ya que justamente son los elementos finitos de la
parte inferior los que soportan la fuerza P (30000 N), además del peso del propio
cuerpo.
• Los valores calculados para la reacción R con el Método de los Elementos Finitos y
con la aplicación del Equilibrio de Fuerzas fueron iguales. Esto habla muy bien del
método; no obstante, lo normal es que existiera una pequeña diferencia entre esos dos
resultados dada la ligera pérdida de exactitud que implica trabajar con métodos
iterativos en softwares, en este caso el Matlab.
Página | 16
