Trabajo Final Elementos Finitos

Powered by TCPDF (www.tcpdf.
org)
Universidad Nacional de San Cristóbal de Huamanga
Facultad de Ingeniería Minas, Geología y Civil
Escuela de Formación Profesional de Ingeniera Civil
CURSO
MÉTODOS NUMÉRICOS (IC-343)
INFORME SOBRE EL PROGRAMA Cálculos por

Elemetos Finitos (Programa hecho en MATLAB y GID)
CATEDRÁTICO :
Ing. CASTRO PEREZ Cristhian
INTEGRANTES:
ACHALMA GALINDO, Melvin código: 16141315

ARESTEGUI TACURI, Renzo Jair código: 16162101
CASTRO CÁRDENAS, Renzo Jair código: 16140411
CUSICHE MENDOZA, Ever J. código: 16140177
HUICHO AUCCATOMA, Michael código: 16147269
LEANDRO PEREZ, Cristhian R. código: 16142832
Ayacucho, 30 de julio del 2018

CONTENIDO
Introducción 2
Presentación 3
Enfoque 4
Objetivos 5
1 FUNDAMENTO TEÓRICO 6
1.1 ELASTICIDAD BIDIMENSIONAL 6

1.2 DEFINICIÓN 6
1.3 CARACTERÍSTICAS 6
1.4 MECANISMO Y CONFORMACIÓN DE FUNCIONAMIENTO
6
1.5 CARACTERÍSTICA DE TRABAJO-TRANSMICIÓN DE CAR-
GAS 7
1.6 CÁLCULO DE LA RESISTENCIA/ROTURA DE UNA VIGA
EN VOLADIZO 7
1.6.1 CÁLCULO DEL MOMENTO . ............................................................................ .7
1.6.2 CÁLCULO DEL MOMENTO DE INERCIA . ........................................ .8
1.7 VENTAJAS DE LA ESTRUCTURA 8
1.8 DESVENTAJAS DE LA ESTRUCTURA 8
1.9 MEJORAS DE LA ESTRUCTURA 8
1.9.1 MODULACIÓN ESTRUCTURAL . ................................................................ .8
1.9.2 CONTROL DE LAS FUERZAS DE REACCIÓN EN SOPORTES
8
1.9.3 SISTEMA POR ZETEO . .............................................................................................. .9
1.9.4 SISTEMA SUPLEX. .......................................................................................................... .9
1.10 ANÁLISIS ESTRUCTURAL 10
1.10.1FUERZAS DE DISEÑO. ........................................................................................... .10
1.10.2DEFORMACIÓN EN VIGAS. ............................................................................ .10
1.10.3MÉTODOS DE CÁLCULO . ............................................................................... .12
1.11 Ejemplo de MEF 14
2 GID 22
1
2.1 Discretización 22
2.2 PROCESOS PARA LA DISCRETIZACION 23
3 ANÁLISIS, PROCEDIMIENTO Y RESULTADOS 33
3.1 Ecuaciones Y Proceso De Uso 33

3.1.1 La Matriz de Propiedades o Matriz Constitutiva. ......................... .33
3.1.2 Funciones de Interpolacion . ............................................................................... .33
3.1.3 Coordenadas de Elemento . .................................................................................. .33
3.1.4 La Matriz de Funciones de Forma . ............................................................. .34
3.1.5 La Matriz de Deformacion Unitaria . .......................................................... .34
3.1.6 La Matriz de Rigidez. ................................................................................................. .34
3.2 Formulación del problema 34
3.3 Propiedades del material 35
3.4 Procedimiento y Calculo mediante GUIDE y MATLAB 35
3.4.1 Identificación de Nudos y Elementos . .................................................... .35
3.4.2 Matriz de Nudos . .......................................................................................................... .35
3.4.3 Matriz de Elementos . ................................................................................................. .35
3.4.4 Coordenadas de los Elementos . ...................................................................... .36
3.4.5 Área de un elemento . ................................................................................................. .36
3.4.6 Matriz de Propiedades . ........................................................................................... .36
3.4.7 Matriz de Restricciones. ........................................................................................... .37
3.4.8 Vector de Fuerzas en cada nudo considerando Restricciones 37
3.4.9 Proceso de Análisis . .................................................................................................... .38
3.4.10Análisis Grafico con el programa GID. .................................................... .41
4 Programa 44
4.1 Código fuente 44

4.1.1 programa Mat-FEM. .................................................................................................... .44
4.1.2 Función TrStif v1.3. ....................................................................................................... .49
4.1.3 Función TrStrs v1.3 . .................................................................................................... .49
4.1.4 Función QdStif v1.3. .................................................................................................... .50
4.1.5 Función QdStrs v1.3 . ................................................................................................. .51
4.1.6 Función constt . ................................................................................................................ .52
Conclusión 54
BIBLIOGRAFÍA 55
INTRODUCCIÓN
El Método de los Elementos Finitos (MEF) constituyen una muy eficiente técnica para resolver problemas
de caracter muy general, existiendo multitud de programas basados en este método.
Debido al extraordinario desarrollo de este método, ciertos problemas que presenta, deben ser analizados
minuciosamente. En general, estos problemas están relacionados con el grado de exactitud de los
resultados. De hecho, durante la pasada década, se han realizado esfuerzos considerables para
determinar la bondad de los resultados y obtener algunas estimaciones prácticas del error producido
en los resultados de un cálculo por Elementos Finitos. La decisión del usuario en la colocación de los
nodos influye de manera directa en la magnitud del error.
Actualmente, existen ciertas investigaciones, dirigidos conseguir métodos automáticos de refinamiento
de una malla existente introduciendo nuevos grados de libertad (más nodos o polinomios de mayor
grado en algunos elementos). Se considera una malla mejorada, aquella en la que la distribución de los
errores es más uniforme. En este trabajo, se tratará de obtener una malla óptima conservando el mismo
número de grados de libertad, siendo la posición de cada nodo la incógnita a resolver.
PRESENTACIÓN
La solución al problema de encontrar la malla óptima en Elementos Finitos (EF), con un determinado
número de nudos de libertad, presenta un indudable interés en la aplicación del método.
En la actualidad, el problema se plantea en términos de un proceso que permite obtener una mejor malla
de elementos finitos a partir de una inicial. La nueva malla se diseña matemáticamente (remallado) de
forma que el error del método sea lo más uniforme posible en todo el dominio de cálculo. Sin embargo,
esta técnica de indudable interés y aplicación, al aumentar el número de grados de libertad (gdl) de la
aproximación, no permite deducir de un modo directo el problema de la malla óptima condicionada
a un número demasiado grande de grados de libertad. Con la solución de este problema se podrán
deducir algunos criterios y recomendaciones para el diseño de una malla de elementos finitos, que le
exigirá, en general, en un proceso de remallado, modificaciones menores.
Para problemas unidimensionales (barras y pilares simples), se pueden encontrar soluciones analiticas.
Para problemas más complicados (tensión y deformación plana), se han utilizado métodos numéricos
para obtener la malla óptima.
ENFOQUE
onsideremos la viga en voladizo que aparece en la figura 1. Se plantea la pregunta por una
C ecuación general que relacione la carga p con el desplazamiento vertical de cualquier punto de
la viga, u( x ). Para ello se debe resolver la ecuación diferencial de la deflexión de la viga
d2 u
EI = −M
dx2
donde M es el momento flector. El resultado es
lx2 x3

1
u( x ) = −
EI 2 6
Figure 1: Viga en Voladizo.
En este cálculo se ha establecido una ecuación de la deflexión de la viga que es válida en toda su
longitud. Para ello se ha empleado una modelación de la estructura como un sistema continuo y, por
tanto, se habla en casos así de modelos estructurales continuos.
OBJETIVOS
Objetivos
Poder determinar y editar el código fuente del programa Cálculo de tensiones y deformaciones de
vigas bidimensionales en voladizo con cargas puntuales o distribuidos, muy usados en Resistencia
de Materiales I utilizando el MÉTODO DE LOS ELEMENTOS FINITOS (FEM) .
1 1
FUNDAMENTO TEÓRICO
ELEMENTOS FINITOS BIDIMENSIONALES

1.1 ELASTICIDAD BIDIMENSIONAL
Los problemas de elasticidad bidimensional son muy frecuentes en Ingeniería, y son asimismo los
primeros en los que se aplicó el MEF. En este caso el medio continuo que se analiza es plano, y se
considera situado en el plano XY. Se denomina t al espesor del dominio en su dirección transversal, el
cual se considera despreciable frente a las dimensiones del dominio en el plano XY.
La posición de un punto está definida por dos coordenadas (x,y), y su deformación tiene dos compo-
nentes u(x,y), v(x,y) en las direcciones x,y respectivamente. El campo de deformaciones es por lo tanto
un vector:
u( x, y)
u=
v( x, y)
Dentro de la elasticidad en dos dimensiones existen dos problemas diferentes:
1. Tensión plana: cuando la tensión σz en sentido perpendicular al plano XY es cero,ya que el sólido
puede dilatarse libremente en el sentido de su espesor. Por lo tanto existe una deformación
unitaria ε z no nula en dicha dirección.
2. Deformación plana: cuando en el sentido del espesor del sólido no hay posibilidad de deformación,
es decir ε z = 0 por lo que se genera una tensión en dicha dirección σz no nula.
En ambos casos la tensión y la deformación en la dirección z no contribuyen a la energía elástica del

sistema.
Las ecuaciones diferenciales que rigen el problema son de orden m=2 en las deformaciones, pues
contienen la derivada primera de las tensiones, que a su vez son derivadas de las deformaciones. En la
expresión del potencial total del sistema aparecen las deformaciones unitarias ε, que son las derivadas
primeras de las deformaciones u, luego el orden de derivación de las incógnitas primarias u en el
potencial es n=1. Se requieren por lo tanto funciones de interpolación con continuidad C0 para asegurar
la convergencia del método.
ESTRUCTURA EN VOLADIZO
1.2 DEFINICIÓN
El termino ménsula refiere a un tipo de viga, denominado más frecuentemente “voladizo”, que se
caracteriza por estar apoyado en solo uno de sus extremos mediante un empotramiento. Debido a
la necesidad de dicho empotramiento, los voladizos suelen ser prolongaciones de vigas continuas de
varios apoyos.
1.3 CARACTERÍSTICAS
Un efecto mecánico que se produce en los voladizos es la deflexión, afectada por:
- Longitud - Peralte y ancho - Material - Localización de la carga - Forma de la sección transversal
1.4 MECANISMO Y CONFORMACIÓN DE FUNCIONAMIENTO

Una carga sobre una viga en voladizo produce dos reacciones en el soporte: la fuerza de corte vertical,
que contrarresta el peso del objeto; y el momento de flexión, que evita que la viga rote.
Ingeniería Civil 6
CAPÍTULO 1 FUNDAMENTO TEÓRICO
1.5 CARACTERÍSTICA DE TRABAJO-TRANSMICIÓN DE CARGAS

Lo que hace el voladizo es hacer contrapeso al claro interior, mientras mayor sea el contrapeso, mas se
reduce la deformación. Lo recomendable seria la igualdad absoluta de los momentos mencionados.
1.6 CÁLCULO DE LA RESISTENCIA/ROTURA DE UNA VIGA EN VOLADIZO

Las vigas son elementos estructurales que han de soportar esfuerzos de flexión. Para el cálculo de una
viga, considerando los esfuerzos de flexión que ha de soportar, se puede usar la fórmula siguiente: M =
R * Z Donde: - M es el momento de las fuerzas que producen la flexión en una sección determinada de
la viga. - Z es el módulo de la sección. - R es el valor de la resistencia del material a tracción.
1.6.1 CÁLCULO DEL MOMENTO

Este momento depende del peso que ha de soportar la viga, de como está situado a lo largo de la misma
y de su longitud. Mientras más larga es la viga y mayor es el peso que ha de soportar, mayor será el
valor del momento M. A continuación se muestran las fórmulas que nos permiten calcular el momento
de las fuerzas exteriores M(máximo momento flexor) para distintos tipos de vigas en voladizo:
Ingeniería Civil 7
1.6.2 CÁLCULO DEL MOMENTO DE INERCIA

Depende de la forma geométrica que tenga la sección de la viga y de sus dimensiones.
1.7 VENTAJAS DE LA ESTRUCTURA

- Es una de las formas de controlar la deformación excesiva de las vigas.
- Utilizar voladizos es una forma eficiente de aumentar la superficie sin ocupar demasiado dentro de un
terreno.
- Se utiliza con frecuencia en la construcción de puentes hasta la aparición de la técnica del puente
colgante.
1.8 DESVENTAJAS DE LA ESTRUCTURA

- Sin protección alguna, los voladizos suelen ser la parte que se degrada con rapidez.
- El sobrepeso en el voladizo provoca deformaciones.
- La humedad y la temperatura son factores que fácilmente afectan a los voladizos.
1.9 MEJORAS DE LA ESTRUCTURA

1.9.1 MODULACIÓN ESTRUCTURAL
Esta forma trapezoidal recta es la más eficiente para una viga en voladizo, ya que el esfuerzo de flexión
permanece relativamente constante en toda la longitud.
1.9.2 CONTROL DE LAS FUERZAS DE REACCIÓN EN SOPORTES

Las fuerzas de levantamiento en las bases tienen que ser evitadas. Estas pueden ser eliminadas mediante
el aumento de la carga hacia abajo en los soportes.
Ingeniería Civil 8
1.9.3 SISTEMA POR ZETEO

Uso del acero como un solo elemento de disipación de cargas, el cual ayuda a transmitir esfuerzos de
tracción , desde el extremo externo superior del voladizo, hacia el centro inferior de la viga en apoyo, lo
cual causa mayor rigidez en la estructura, disminuyendo el esfuerzo flector del segmento en apoyos por
tracción de el alma.
1.9.4 SISTEMA SUPLEX

El de mayor empleo en la construcción de artículos en voladizo. Consiste en colocar el acero mayor de
manera convencional (Superior abajo e inferior arriba) a lo largo de todo el tramo, localizando en la
parte superior y central entre el encuentro de las 2 vigas, con longitud total de uniones entre centros de
las mismas, una barra de acero de mayor grosor, a fin de retener el diagrama de vuelco
Ingeniería Civil 9
1.10 ANÁLISIS ESTRUCTURAL

El predimensionado de las vigas consiste en determinar lasdimensiones necesarias para que el elemento
sea capaz de resistir laflexión y el corte, así como también debe tener dimensiones tales que laflecha no
sea excesiva. Así, el esquema para cumplir con los requisitos deuna viga consiste en:
-Determinar las cargas
-Cuantificar las fuerzas de diseño
-Predimensionar mediante criterio de Resistencia
-Comprobar las dimensiones por rigidez
1.10.1 FUERZAS DE DISEÑO

Los efectos que producen las cargas sobre una viga son de dos tipos: Fuerza Cortante (V) y Momento
Flector (M ). La magnitud de estas fuerzas son variables a lo largo de la longitud de la viga, siendo así
el objetivo principal de determinar la magnitud de la fuerza cortante y el momento flector máximo
aplicado en la viga (Vmax ; Mmax). El procedimiento básico para cuantificar las fuerzas de diseño
consiste en:
1.Aislar el elemento del sistema estructural
2.Determinar las reacciones por las ecuaciones estáticas o de lascondiciones de apoyos.
3.Realizar un corte en la sección donde se desea conocer la magnitudde las fuerzas internas con un
plano perpendicular al eje del elemento.
4.Las fuerzas internas se obtienen de aplicar el equilibrio sobre una de lasdos porciones obtenidas por
el corte.
1.10.2 DEFORMACIÓN EN VIGAS

LINEA ELASTICA
Denominaremos linea elástica a la curva que forma la fibra neutra una vez cargada la viga,considerando
que ésta se encontraba inicialmente recta
SUPUESTO BASE
Para establecer una serie de relaciones al interior de la sección,indicamos que se trata de una viga,cuyo
material se encuentra solicitado dentro del rango de proporcionalidad entre tensiones y de formaciones,
y en donde seadmite la conservacion de lascaras planas. Dicho en otra forma,donde se cumplen la ley
de hookey la hipótesis de Bernouilli-Navier.
Ingeniería Civil 10
LEY DE HOOKE
Establece que la relación entre la tensión y la deformación unitaria es una constante y se denomina
modulo de elasticidad.
τ
E=
ε
Donde:
E=Elasticidad
τ=Tensión
ε=Deformación Unitaria
o otra forma de expresarlo
τ=E×ε
FÓRMULA DE FLEXIÓN
ε = MV
EI
Donde:
τ=Tensión
ε=Deformación Unitaria
E=Elasticidad
M=Momento flector
V=Distancia desde la fibra neutra más traccionada o comprimida
I=Inercia
ANÁLISIS DE LA SECCIÓN
La sección cc’ tt’, inicialmente recta, se curva con un radio R como indica la figura.
La fibra cc’ se acorta a cc”
La fibra tt’ se alarga a tt” y
La fibra nn’ permanece del mismo largo.
Mediante semejanza de triángulos e igualando expresiones se obtiene:
Mdx
dφ =
EI
1.10.3 MÉTODOS DE CÁLCULO

Nos enfocaremos al caso de VIGA EN VOLADIZO CON CARGA UNIFORME REPARTIDA
MÉTODO DE ÁREA DEL MOMENTO

Establecemos el equilibrio externo:
Ra = qL
Determinamos la ecuación general del momento flector:
qx2
Mx = −
2
El ángulo entre las tangentes trazadas en ambos extremos de la viga lo obtenemos aplicando el Primer
Teorema de Mohr.
1 R L qx2
φOL = − dx
EI 0 2
−qL3
φA = −
EI
Calculando la desviación tangencial en 0(extremo libre de la viga) con respecto a la tangente trazada en
el otro extremo,determinamos la flecha máxima.
1 R L qx2
tO − L = − xdx
EI 0 2
qL4
Ymax = tO− L = −
8EI
MÉTODO DE DOBLE INTEGRACIÓN

Con la ecuación general de momento flector establecemos la ecuación diferencial de la elástica.
d2 y qx2
EI 2 = − 2
dx
Integrando la ecuación diferencial dos veces se obtiene:
qx4
EIy = − + C1 x + C2
24
Tomando condiciones iniciales(flecha y pendiente nula):
qL3
C1 =
6
qL4
C2 = −
8
Finalmente:
Ecuación general de la pendiente
dy qx3 qL3
EI dx = − 6 + 6
Ecuación general de la flecha
qx4 qL3 x qL4
EIy = − + −
24 6 8
MÉTODO DE LA VIGA CONJUGADA
Con el gráfico de momento flector y los valores característicos generamos la viga ficticia.
qx2
Mx = −
2
A la viga ficticia se le aplica como carga el momento flector de la viga dada divido por EI.
La relación establecida entre la viga ficticia y la viga real es que los valores de cortante y momento de la
viga ficticia equivalen a la pendiente y a la flecha de la viga real.
Pero en el caso particular de las vigas en voladizo, la pendiente en el apoyo es nula, así como su
descenso. En este punto no deberían existir R ni M por lo tanto para la aplicación de este método, es
necesario invertir el apoyo de la viga ficticia al otro extremo de la viga, de manera de encontrar R Y
M(maxima) en el punto correspondiente.
Mmax qL2
q0 = =
EI 2EI
qL2 L
φ A = Ra0 =
2EI 3
qL3
φA =
6EI
qL2 L 3L
Ymax = M0max =
2EI 3 4
qL 4
Ymax =
8EI
1.11 Ejemplo de MEF

Como ya se menciono el método de elementos finitos nos facilita la resolución de Ecuaciones Diferen-
ciales Parciales, por ello a continuación se describirá como se desarrolla el procedimiento para dicha
solución.
Problema 1
Suponga que para resolver la ecuación de Poisson, ∇2 u = p( x, y),en un medio bidimensional se

utilizan mallas de elementos finitos triangulares con 3 nudos, como se muestra en las figuras, con
∆x = ∆y = a. Determine las ecuaciones típicas que resultan en cada caso (correspondientes al nudo
E en la primera primera malla y a los nudos F y G en el caso de la segunda).
SOLUCIÓN
SOLUCIÓN
En la ecuación diferencial ∇ T K (∇u) + λu + p = 0 en Ω. Empleando el método de Galerkin se obtuvo
que:
Z Z Z Z Z
( B T KBa − N T λNa)dΩ = N T pdΩ + N T q̄dS
Ω Ω Sq
Ha=b
Z Z
H= B T KBdΩ − N T λNdΩ
Ω Ω
Z Z
b= N T pdΩ − N T q̄dS
Ω Sq
Empleando Elemento Finito Triangular (CST)

u = N1 ( x, y)u1 + N2 ( x, y)u2 + N3 ( x, y)u3

 ai = x j y m − x m y j


1
Ni = 2A ( ai + bi x + ci y) bi = y j − ym

c = x − x

i m j
i, j, m −→ son permutaciones cíclicas de los índices 1,2,3

1 1 1
1
A = xi xj xm
2
yi yj xm
u = N1 u1 + N2 u2 + N3 u3 = Na(e)
! ! u  
u1

∂u ∂N1 ∂N2 ∂N3 1
∇u = ∂x
∂u = ∂x
∂N1
∂x
∂N2
∂x
∂N3
 u2  = 1 b1 b2 b3  u2  = Ba(e)
∂y 2A c1 c2 c3
∂y ∂y ∂y u3 u3
Z
H = ∑ H (e) = ∑ ( B T KB − N T λN )dΩ
e e Ωe
 
Z 1
pA 
Z
b = ∑ b(e) = ∑ N T pdΩ + N T q̄dq = be = 1  + t.b.
e e Ωe Sqe 3
1
 2   2 
b b1 b2 b1 b3 c1 c1 c2 c1 c3
Kxx  1 K yy
H= b1 b2 b22 b2 b3  +  c1 c2 c2
2 c2 c3  +
4A 4A
b1 b3 b2 b3 b32 c1 c3 c2 c3 c23
 
2b1 c1 b1 c2 + b2 c1 b1 c3 + b3 c1
Kxy
 b1 c2 + b2 c1 2b2 c2 b2 c3 + b3 c2 
4A
b1 c3 + b3 c1 b2 c3 + b3 c2 2b3 c3
Para el caso de materiales homogéneos Kxy = 0. Así mismo para simplificar la matriz de los elementos,
asumimos que k = Kxx = Kyy .
Finalmente la matriz H se reduce a:
b12 + c21
 
b1 b2 + c1 c2 b1 b3 + c1 c3
k 
H= b1 b2 + c1 c2 b22 + c22 b2 b3 + c2 c3 
2A
b1 b3 + c1 c3 b2 b3 + c2 c3 b32 + c23
SOLUCIÓN FIGURA 01
En la figura se muestran 6 triángulos que conectan al punto E, de los cuales los triángulos 1,3 y 5
presentan los mismos valores para bi y ci , así mismo para los triángulos 2,4 y 6.
Cálculo para los elementos 1,3 y 5


 bi = y j − y m = a − a = 0





 b j = y m − yi = a − 0 = a

 bm = y i − y j = 0 − a = − a



H (1) = H (3) = H (5)

ci = x m − x j = 0 − a = − a









 c j = xi − x m = 0 − 0 = 0

c m = x j − xi = a − 0 = a

Por lo tanto:
 2
− a2 0 −1
  
a 0 1
k
H (1) = H (3) = H (5) =  0 a2 − a2  = k  0 1 −1 
2A
− a2 − a2 2a 2 −1 −1 2



 bi = y j − y m = 0 − a = − a





 b j = y m − yi = a − 0 = a

 bm = y i − y j = 0 − 0 = 0



H (2) = H (4) = H (6)

 ci = x m − x j = a − a = 0





 c j = xi − x m = 0 − a = − a




c m = x j − xi = a − 0 = a

Por lo tanto:
a2 − a2 −1 0
   
0 1
(2) (4) (6) k 
H =H =H = − a2 2a2 − a2  = k  −1 2 −1 
2A
0 − a2 a2 0 −1 1
Ensamblando las matrices con los correspondientes Nudos

 
1+1 −1 −1 0 0 0 0 A

 −1 1+2 0 −1 − 1 0 0 0 
 B

 −1 0 2+1 −1 − 1 0 0 0 
 D
H = k 0 −1 − 1 −1 − 1 1 + 2 + 1 + 1 + 2 + 1 −1 − 1 −1 − 1 0 E
 

−1 − 1 −1

1+2

 0 0 0 0  F
 
 0 0 0 −1 − 1 0 2+1 −1  H
0 0 0 0 −1 −1 1+1 I
 
2 −1 −1 0 0 0 0 A

 −1 3 0 −2 0 0 0  B


 −1 0 3 −2 0 0  D
0 
H = k 0 −2 −2 8 −2 −2 0  E
 
0 −2 3 0 −1  F
 
 0 0
 
 0 0 0 −2 0 3 −1  H
0 0 0 0 −1 −1 2 I
También el vector b ensamblado se tiene:
 
1
pA
be =  1 
3
1
;  
2

 2 

2
 2 
pa  
b= 6
 
6 
 

 2 
 
 2 
2
Formando la ecuación final en Hu = b
    
2 −1 −1 0 0 0 0 uA 2
 −1 3 0 −2 0 0 0  uB   2 
    
 −1 0 3 −2 0 0 0  uD   2 
   pa2  
k  0 −2 −2 8 −2 −2 0 uE = 6
    
6 
  
0 −2 −1
   
 0 0 3 0  uF   2 
    
 0 0 0 −2 0 3 −1  uH   2 
0 0 0 0 −1 −1 2 uI 2
Finalmente la Ecuación en el Nudo E:
pa2
−2u B − 2u D + 8u E − 2u F − 2u H =
k
Donde:
pa2
k + 2u B + 2u D + 2u F + 2u H
uE =
8
SOLUCIÓN FIGURA 02
En la figura 02 se muestran a todos los elementos que relacionan a los nudos F y G. Los elementos
2,3,6,9 y 10 se repiten en la figura 01 analizada anteriormente.
Por lo tanto:
−1
 
1 0
H (2) = H (9) = k  0 1 −1 
−1 −1 2
−1 0
 
1
H (3) = H (6) = H (10) = k  −1 2 −1 
0 −1 1


 bi = y j − y m = a − a = 0



 b j = y m − yi = a − 0 = a




 bm = y i − y j = 0 − a = − a



H (1) = H (5) = H (8)

ci = x m − x j = 0 − a = − a









 c j = xi − x m = a − 0 = a

c m = x j − xi = a − a = 0

Por lo tanto:
 2
− a2 −1 0
  
a 0 1
k
H (1) = H (5) = H (8) =  − a2 2a2 − a2  = k  −1 2 −1 
4A
0 − a2 a2 0 −1 1
Cálculo para los elementos 4 y 7



 bi = y j − y m = 0 − a = − a





 b j = y m − yi = a − 0 = a

 bm = y i − y j = 0 − 0 = 0



H (4) = H (7)

ci = x m − x j = 0 − a = − a









 c j = xi − x m = 0 − 0 = 0

c m = x j − xi = a − 0 = a

Por lo tanto:
2a2 − a2 − a2 2 −1 −1
   
(4) (7) k
H =H =  − a2 a2 0  = k  −1 1 0 
4A
− a2 0 a 2 −1 0 1
Ensamblando las matrices con los correspondientes Nudos
 
3 −1 0 0 −2 0 0 0 0 0 B
 −1 4 −1 0 0 −2 0 0 0 0  C
 
 0 −1 2 0 0 0 −1 0 0 0 D
 

2 −2 0
 
 0 0 0 0 0 0 0  E
 
 −2 0 0 − 2 8 −2 0 −2 0 0  F
H = k  0 −2 0

 0 −2 7 −2 0 −2 0 
 G
 0
 0 −1 0 0 −2 4 0 0 −1 
 H
 0 0 0 0 −2 0 0 3 −1 0 J
 

0 −2 −1 −1
 
 0 0 0 0 0 4  K
0 0 0 0 0 0 −1 0 −1 2 L
También el vector b ensamblado se tiene:
 
1
pA
be =  1 
3
1
;  
3

 2 

2
 
 
 
 2 
2
 
pa  4 
b=  
6  8 


 2 

3
 
 
 
 2 
2
Formando la ecuación final en Hu = b
    
3 −1 0 0 −2 0 0 0 0 0 uB 3
 −1 4 −1 0 0 −2 0 0 0 0  uC   2 
    
 0 −1 2 0 0 0 −1 0 0 0 uD 2
    
   
2 −2
    
 0 0 0 0 0 0 0 0  uE   2 
 pa2 
    
 −2 0 0 − 2 8 −2 0 −2 0 0  uF 4 
k  =  
 0 −2 0 0 −2 7 −2 0 −2 0  uG  6  8 
    
 0
 0 −1 0 0 −2 4 0 0 −1 
 uH 


 2 

 0 0 0 0 −2 0 0 3 −1 0 uJ 3
    
   
−2 −1 −1
    
 0 0 0 0 0 0 4  uK   2 
0 0 0 0 0 0 −1 0 −1 2 uL 2
La Ecuación en el Nudo F:
2pa2
−2u B − 2u E + 8u F − 2uG − 2u J =
3k
Donde:
2pa2
3k + 2u B + 2u E + 2uG + 2u J
uF =
8
La Ecuación en el Nudo G:
4pa2
−2uC − 2u F + 7uG − 2u H − 2uK =
3k
Donde:
4pa2
3k + 2uC + 2u F + 2u H + 2uK
uG =
7
2
1. La ecuación diferencial D ddxT2 = 0 gobierna la distribución de temperaturas en el estado estacionario en los
muros compuestos mostrados en las figuras. D es la conductividad térmica. Utilizando la hoja de cálculo en
xls, determine en cada caso la variación de temperatura en el interior del muro y el flujo de calor por unidad
de área.
SOLUCIÓN
Utilizando la hoja de cálculo proporcionada en clase, se tienen los resultados:

Solución a:
n x X u e i j p pu’
1 0.0 x 20.000 1 1 2 0.020 -0.015
2 1.3 19.044 2 2 3 0.005 -0.015
3 9.3 -4.493 3 3 4 0.004 -0.015
4 11.8 x -15.000
Solución b:
n x X u e i j p pu’
1 0.0 x 30.000 1 1 2 2.0 -0.784
2 1.0 29.608 2 2 3 1.0 -0.784
3 6.0 25.686 3 3 4 0.2 -0.784
4 10.0 x 10.000
2
2 1
GID
ELEMENTOS FINITOS BIDIMENSIONALES

2.1 Discretización
La discretización es un proceso de modelación que consiste en la división de un cuerpo, de manera
equivalente, en un sistema que se conforma por cuerpos de menor tamaño a los que le damos el nombre
de “Elementos Finitos”, que son conectados entre sí por medio de puntos comunes llamados nodos,
estos a su vez pueden llegar a formar superficies y volúmenes de control completamente independientes
entre sí, y que son afectados por las condiciones de frontera que afectan el cuerpo estudiado.
Un modelo de elementos finitos es una representación discreta de una parte continua, que se va a
analizar, como se muestra en la figura .Esto se lleva mediante la utilización de nodos, los cuales se van
conectando para formar elementos. Un nodo es la representación discreta de esa parte continua que se
analizara y donde se predice el resultado de acuerdo a las cargas aplicadas.
La malla comprende a todo el conjunto de elementos que une a los nodos, es debido a ella que se
podrán transferir las cargas a través de los elementos, pero para que estas cargas se puedan transferir
de una parte a otra de manera satisfactoria se necesita que todos los nodos sean comunes donde estas
se junten .La respuesta que el análisis arroje estará ligada a los grados de libertad a los que los nodos
estarán sujetos.
Después de someter al análisis la viga modelada a partir de las suposiciones de los materiales y el
entorno ,el diseñador se dará cuenta en términos simples cual es el punto más débil de la viga, cual es
el comportamiento de la misma y bajo qué condiciones de trabajo se presenta la falla.
CAPÍTULO 2 GID
En el presente trabajo realizaremos el proceso de discretización mediante el programa GID, el cual se

explicara a detalle a continuación.
MANEJO DEL SOFTWARE GID

2.2 PROCESOS PARA LA DISCRETIZACION
Es necesario usar este programa para la discretizacion de la viga, esto debido a que la discretizacion
manual es muy tedioso. Con el GID podremos discretizar y esos datos lo usaremos en el MATLAB, a
continuacion se mostraran los pasos a seguir para poder discretizar la viga:
1. Dibujaremos la viga en voladizo(estara representado por tres apoyos fijos en dicha seccion), utlizando las
coordenadas de cada vertice y de cada apoyo fijo.
CAPÍTULO 2 GID
2. Creamos la superficie, ya que anteriormente solo eran lineas unidas.
CAPÍTULO 2 GID
3. Cargamos el programa "Mat-femv1.3"
CAPÍTULO 2 GID
4. Dibujamos los apoyos.
CAPÍTULO 2 GID
5. Dibujamos la carga distribuida.
CAPÍTULO 2 GID
6. Introducimos las propiedades de los elementos.
7. Ponemos las unidades para los esfuerzos y deformaciones.
CAPÍTULO 2 GID
8. Discretizamos la superficie de la viga.
CAPÍTULO 2 GID
9. Guardamos el archivo donde se encuentra el programa "Mat-femv1.3" con el nombre que queremos con la
extension ".m" para guardarlo en el formato de MATLAB.
CAPÍTULO 2 GID
CAPÍTULO 2 GID
2
3
3 1
ANÁLISIS, PROCEDIMIENTO Y RESULTADOS
3.1 Ecuaciones Y Proceso De Uso

3.1.1 La Matriz de Propiedades o Matriz Constitutiva
La matriz se ve reflejada mediante las propiedades del material, siendo esta constante para todo el
proceso iterativo que se quiera utilizar.Y se expresa de la siguiente manera:
 υ υ 
1 (1− υ )
0 (1− υ )
 υ υ 
E (1− υ )  (1− υ ) 1 0 (1− υ ) 
D = (1+υ)(1−2υ)  
(1−2υ)

 0 0 2(1− υ )
0  
υ υ
(1− υ ) (1− υ )
0 1
Donde:
E=Modulo de elasticidad
υ=Coeficiente de Poisson
Cabe recalcar que esta matriz es único para una estructura de análisis aximetrico,sin embargo ca-
sos con tension plana o deformación plana llevan una diferente constitución de la matriz D.
3.1.2 Funciones de Interpolacion

Las funciones de interpolacion depende de la forma y nudos que se quiera adoptar a momento de
transformar un elemento cualquiera a uno simétrico. Para este caso con elementos triangulares de tres
nudos, emplearemos:
1
N1 = 2A (r2 z3 − r3 z2 + R(z2 − z3 ) + Z (r3 − r2 ))
1
N2 = 2A (r3 z1 − r1 z3 + R(z3 − z1 ) + Z (r1 − r3 ))
1
N3 = 2A (r1 z2 − r2 z1 + R(z1 − z2 ) + Z (r2 − r1 ))
Donde:
 
1 r1 z1
2A = 1 r2 z2 
1 r3 z3
Ni =Función de interpolacion i
A=área del elemento
X=coordenadas locales de las abscisas en el elemento transformado
Y=coordenadas locales de las ordenadas en el elemento transformado
3.1.3 Coordenadas de Elemento

Las coordenadas de elementos están expresadas en funcion de las Funciones de Interpolacion y Las
coordenadas de los elementos en estado real.
R( x, y) = N1 ( x, y)r1 + N2 ( x, y)r2 + N3 ( x, y)r3
Z ( x, y) = N1 ( x, y)z1 + N2 ( x, y)z2 + N3 ( x, y)z3
CAPÍTULO 3 ANÁLISIS, PROCEDIMIENTO Y RESULTADOS
3.1.4 La Matriz de Funciones de Forma

Esta dado por la siguiente expresión:

N1 0 N2 0 N3 0
N=
0 N1 0 N2 0 N3
Que es la constitucion de las funciones de forma.
3.1.5 La Matriz de Deformacion Unitaria

La matriz de deformación es la que proviene de la matriz de función de forma y la transformación
mediante la matriz Jacobiana; dando como resultado para elementos triangulares de tres nudos de la
siguiente manera:
  −1
1 r1 z1 0 0 0
0 1 0 0 0 0 0 0 0 1 r1 z1 
   
1 z
1 r 0 0 0 1 r2 z2 0 0 0 
 
1 r

B = 2A ∗
 0 0 0 0 0 1 0 0 0 1 r2 z2 


 
0 0 1 0 1 0 1 r3 z3 0 1 0 
0 0 0 1 r3 z3
3.1.6 La Matriz de Rigidez

La matriz de rigidez elemental para un elemento de tras nudos esta compuesta de los parámetros
anteriores como son: La matriz constitutiva, la matriz de deformación unitaria, y el volumen diferencial
de cada elemento.
Ke = B T DBdυ
H
y caso especifico para un solido de revolución analizado por elementos triangulares de tres nudos, se
realiza mediante una integración de métodos numéricos, dando resultados convenientes:
Ke = A ∗ θ ∗ R g ∗ B T ∗ D ∗ B
Donde:
A:Área de un elemento
B:Matriz de deformación unitaria que depende de Rg y Zg que son los centros de
gravedad de cada elemento
D:Matriz constitutiva
θ:Angulo de revolución del solido
r:Coordenada radial del Solido de revolución
3.2 Formulación del problema
Problema 2
Resolver mediante el método de elementos finitos bidimensionales una viga en voladizo. con las
siguientes características.
W Kg/m
SOLUCIÓN
3.3 Propiedades del material
3.4 Procedimiento y Calculo mediante GUIDE y MATLAB

3.4.1 Identificación de Nudos y Elementos
En este caso se ha programado de modo que todos los datos generalmente utilizados se generen
automáticamente a partir de los datos básicos ingresados como:
cantidad de nudos = nn = (ix + 1) ∗ (iy + 1)
Cantidad de elementos = ne = 2 ∗ ix ∗ iy
3.4.2 Matriz de Nudos

Como se muestra, se generan las coordenadas de los nudos a partir del algoritmo siguiente.
3.4.3 Matriz de Elementos

De igual manera, se ha programado de modo que el programa genere automaticamente la matriz de
elementos:
y los valores para una división de 12 nudos con 12 elementos:
 
11 10 8

 10 9 7 


 2 5 3 

8 10 7
 
 
 
 8 7 6 
 

 7 9 4 


 7 4 6 


 2 3 4 

2 4 1
 
 
 
 3 5 6 
 
 3 6 4 
8 6 11
3.4.4 Coordenadas de los Elementos

Se realiza la programación de obtención de la matriz de coordenadas de cada elemento:
 
0 0.03

 0 0.015 

0.01558 0.0075
 
 
0.03 0.03
 
 
 

 0 0 


 0.03 0 


 0.04442 0.0225 

0.04861 0.00937
 
 
0.06 0.03
 
 
 
 0.06 0.015 
0.06 0
3.4.5 Área de un elemento

El área de un elemento triangular cualquiera se obtiene mediante determinantes:
 
1 x1(m) y1(m)
1
A(m) = det 1 x2(m) y2(m)
2
1 x3(m) y3(m)
3.4.6 Matriz de Propiedades

Numéricamente es igual a:
Elemento 1
 
493.1358 −184.7491 −413.0847 252.7811 −80.0511 −68.0321

 −184.7491 323.2693 254.745 −467.2367 −69.9959 143.9674 


 −413.0847 254.745 450.8284 −256.2401 −37.7437 1.4951


252.7811 −467.2367 −256.2401 835.1885 3.4589 −367.9517
 
 
 
 −80.0511 −69.9959 −37.7437 3.4589 117.7948 66.537 
−68.0321 143.9674 1.4951 −367.9517 66.537 223.9843
Elemento 2
 
263.1579 0 0 −131.5789 263.1579 131.5789

 0 785.546 −129.6151 0 129.6151 −785.546 

 0 −129.6151 196.3865 0 −196.3865 129.6151 

−131.5789 0 0 65.7895 131.5789 −65.7895 
 

 
 −263.1579 129.6151 −196.3865 131.5789 459.5444 −261.194 
131.5789 −785.546 129.6151 −65.7895 −261.194 851.3354
Elemento 12
 
172.7814 −98.0348 21.3636 −33.5441 −194.145 131.5789

 −98.0348 322.6568 −31.5803 −257.6182 129.6151 −65.0386 

 21.3636 −31.5803 301.7516 163.1592 −323.1153 −131.5789 

−33.5441 −257.6182 163.1592 365.8618 −129.6151 −108.2436 
 

 
 −194.145 129.6151 −323.1153 −129.6151 517.2603 0 
131.5789 −65.0386 −131.5789 −108.2436 0 173.2822
3.4.7 Matriz de Restricciones

Representa a los nudos que se han sujeto a cierta restricción. En este caso representa a los nudos que se
encuentran en la base del muro de sostención.
Las restricciones en sus respectivos nudos.
 
1 1 0

 1 2 0 

2 1 0
 
 
2 2 0
 
 
 
 5 1 0 
5 2 0
3.4.8 Vector de Fuerzas en cada nudo considerando Restricciones

Según lo desarrollado en la parte teórica, tenemos:
Y las fuerzas generadas incluyendo la restricción:
0 0
 

 0 0 


 0 0 


 0 0 

0 0
 
 
 
 0 0 
 

 : : 

 −10 −10 
: :
3.4.9 Proceso de Análisis

Grados de Libertad
Se ensamblan teniendo en cuenta los movimientos permitidos en cada elemento, para el primer
elementos del ejemplo tenemos:
Para el elemento 5 por ejemplo:
Matriz de Rigidez Total Restringido

Se ensambla sumando todas las matrices de rigideces elementales y colocando con un arreglo que
depende del vector de grados de libertad. Finalmente queda así:
Obteniendo la matriz de rigidez total, considerando las restricciones para una discretización de 12
elementos tenemos:
Que se convierte en una matriz de 22x22:
 
−0.4595 −0.2612 0.2632 0.1296 0 0 0.1964 0.1316 0 0 0 0 0 0 0 0 0 0 0 0 0 0

 −0.2612 0.8513 0.1316 −0.7855 0 0 0.1296 −0.0658 0 0 0 0 0 0 0 0 0 0 0 0 0 0 


 −0.2632 0.1316 0.9875 −0.3153 −0.6022 0.3844 −0.0801 −0.1996 −0.0421 −0.0421 0 0 0 0 0 0 0 0 0 0 0 0 

0.1296 −0.7855 −0.3153 1.5484 0.3844 −0.5306 −0.1996 0.144 0.001 0.001 0 0 0 0 0 0 0 0 0 0 0 0
 
 
−0.6022 0.3844 −0.2562 −0.242 −0.1281 −0.4421 −0.4421 −0.8862
 
 0 0 2.1724 0.2612 0 0 0 0 0 0 0 0 0 0 
 

 0 0 0.3844 −0.5306 −0.2562 2.8067 −0.1281 −0.4364 −0.2612 −0.2612 0.2612 −1.0212 0 0 0 0 0 0 0 0 0 0 


 −0.1964 0.1296 −0.0801 −0.1996 −0.242 −0.1281 1.1915 −0.1997 0 0 0.1799 0.1306 −0.8862 0.2612 0 0 0.0332 0.006 0 0 0 0 


 0.1316 −0.0658 −0.1996 0.144 −0.1281 −0.4364 −0.1997 1.2951 0 0 0.1306 −0.2749 0.2612 −1.0212 0 0 0.0041 0.3593 0 0 0 0 

0 0 −0.0421 0.001 −0.4421 −0.2612 0 0 0.451 0.451 0.0332 0.0041 0 0 0 0 0 0 0 0 0 0
 
 
−0.001 −0.3763 −0.2612 −0.8185
 
 0 0 0 0 0.2561 0.2561 0.006 0.3593 0 0 0 0 0 0 0 0 0 0 
 

 0 0 0 0 −0.8862 0.2612 0.1799 0.1306 0.0332 0.0332 1.3077 −0.1164 −0.0701 −0.225 −0.541 −0.0868 0 0 0 0 −0.0235 0.0304 


 0 0 0 0 0.2612 −1.0212 0.1306 −0.2749 0.0041 0.0041 −0.1164 1.332 −0.225 −0.1193 −0.0868 −0.5308 0 0 0 0 0.0324 0.2551 


 0 0 0 0 0 0 −0.8862 0.2612 0 0 −0.0701 −0.225 2.1744 −0.2573 −0.4555 0.2568 −0.4421 −0.2612 −0.3205 0.2256 0 0 

0 0 0 0 0 0 0.2612 −1.0212 0 0 −0.225 −0.1193 −0.2573 2.8103 0.2568 −0.8399 −0.2612 −0.8185 0.2256 −0.0113 0 0
 
 
−0.541 −0.0868 −0.4555 −0.1225 −0.5961 −0.3086 −0.4555
 
 0 0 0 0 0 0 0 0 0 0 0.2568 2.0482 0 0 0.2612 
 

 0 0 0 0 0 0 0 0 0 0 −0.0868 −0.5308 0.2568 −0.8399 −0.1225 2.5471 0 0 −0.3086 −0.3312 0.2612 −0.8451 


 0 0 0 0 0 0 0.0332 0.0041 0 0 0 0 −0.4421 −0.2612 0 0 0.451 0.2561 −0.0421 0.001 0 0 


 0 0 0 0 0 0 0.006 0.3593 0 0 0 0 −0.2612 −0.8185 0 0 0.2561 0.8355 −0.001 −0.3763 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 −0.3205 0.2256 −0.5961 −0.3086 −0.0421 −0.001 0.9374 0.1156 0.0214 −0.0316
 
 
−0.0113 −0.3086 −0.3312 −0.3763 −0.0335 −0.2576
 
 0 0 0 0 0 0 0 0 0 0 0 0 0.2256 0.001 0.1156 0.9763 
 
 0 0 0 0 0 0 0 0 0 0 −0.0235 0.0324 0 0 −0.4555 0.2612 0 0 0.0214 −0.0335 0.4576 −0.2601 
0 0 0 0 0 0 0 0 0 0 0.0304 0.2551 0 0 0.2612 −0.8451 0 0 −0.0316 −0.2576 −0.2601 0.8477
Resultado: Deformaciones de los nudos libres

Resulta a partir de la matriz de rigidez y el vector de fuerza por la siguiente ecuación matricial:
q = inv(K ) · F
Obteniendo:
 
(1, 5) −0.0168

 (1, 6) −0.0446 


 (1, 7) 0.0573


(1, 8) −0.092
 
 
 

 (1, 11) −0.0606 


 (1, 12) −0.1019 


 (1, 13) 0.0292 


 (1, 14) −0.1817 

(1, 15) −0.0287
 
 
 

 (1, 16) −0.2049 


 (1, 17) 0.0859 


 (1, 18) −0.2914 


 (1, 19) −0.002 

(1, 20) −0.2739
 
 
 
 (1, 21) −0.0815 
(1, 22) −0.2709
Resultado: Reacciones Resulta a partir de la matriz de rigidez, el desplazamiento y el vector de fuerza
por la siguiente ecuaci´on matricial:
R=K·q−F
Obteniendo:
 
(1, 1) −23.3537

 (2, 1) 13.4770 

 (3, 1) 6.7073 

(4, 1) −7.4329 
 

 
 (9, 1) 16.6463 
(10, 1) 3.9559
Resultado: Las Tensiones
Nos resulta del siguiente algoritmo de desplazamientos
 
1.5002 0.4950 −0.8068

 0.3025 0.0521 −0.5134 


 −0.6471 0.0010 −0.2693 

0.3504 0.0650 −0.3181
 
 
 

 −1.0726 0.0314 −0.5662 


 −0.5357 0.0969 −0.2922 


 0.1541 −0.2004 −0.3478 


 −0.2569 −0.1094 −0.2835 

0.4392 −0.5497 −0.2218
 
 
 
 0.0517 −0.3339 −0.2462 
−0.3808 −0.1261 −0.1987
3.4.10 Análisis Grafico con el programa GID
Desplazamientos
Figure 3.1: Desplazamiento en X
Figure 3.2: Desplazamiento en Y
Figure 3.3: Deformada de la viga en voladizo
Esfuerzos
Figure 3.4: Esfuerzos generados en la dirección X
Figure 3.5: Esfuerzos generados en la dirección Y
2
3
4 1
4 Programa
4.1 Código fuente

4.1.1 programa Mat-FEM
Inicializamos el programa MAT-FEM de nuestra carpeta de archivos "PROGRAMA CÁLCULO CON
ELEMENTOS FINITOS y luego MAT − f emv 1.3"
Figure 4.1: Matfem-v1-3
CAPÍTULO 4 Programa
Figure 4.2: RUN en MATLAB
Código Fuente 4.1: MAT-fem

1 % % MAT - fem
2 clear
3 clc
4 file_name = input ( ' Enter the file name : ' , 's ') ;
5
6 tic ; % Start clock
7 ttim = 0; % Initialize time counter
8 eval ( file_name ) ; % Read input file
9
10 % Dimensiones basicas .
11 coordinates ;
12 elements ;
13 npnod = size ( coordinates ,1) ; % Number of nodes
14 nndof = 2* npnod ; % Number of total DOF
15 nelem = size ( elements ,1) ; % Number of elements
16 nnode = size ( elements ,2) ; % Number of nodes per element
17 neleq = nnode *2; % Number of DOF per element
18
19 elements = sortrows ( elements ) ;
20
21 ttim = timing ( ' Time needed to read the input file ' , ttim ) ; ←-
% Reporting time
22
23 % Dimension de las matrices globales
24 StifMat = sparse ( nndof , nndof ) ; % Create the global ←-
stiffness matrix
25 force = sparse ( nndof , 1 ) ; % Create the global force ←-

vector
26 force1 = sparse ( nndof , 1 ) ; % Create the global force ←-
vector
27 reaction = sparse ( nndof , 1 ) ; % Create the global reaction ←-
vector
28 u = sparse ( nndof , 1 ) ; % Nodal variables
29
30 % Propiedades del material ( constante sobre el dominio )
31 dmat = constt ( young , poiss , pstrs ) ;
32
33 ttim = timing ( ' Time needed to set initial values ' , ttim ) ;
34
35 % Ciclo del elemento
36 for ielem = 1 : nelem
37
38 % Recuperar propiedades del elemento
39 lnods = elements ( ielem ,:) ; % Elem . ←-
connectivity
40 coord (1: nnode ,:) = coordinates ( lnods (1: nnode ) ,:) ; % Elem . ←-
coordinates
41
42 % Evaluar la matriz de rigidez elemental y el vector de fuerza de ←-
masa
43 if ( nnode == 3) % 3 Nds . Triangle
44 [ ElemMat , ElemFor ] = TrStif_v1_3 ( coord , dmat , thick , denss ) ;
45 else % 4 Nds . Quadrilateral
46 [ ElemMat , ElemFor ] = QdStif_v1_3 ( coord , dmat , thick , denss ) ;
47 end
48
49 % Encuentre la lista de numeros de ecuaciones para el elemento ←-
i - esimo
50 eqnum = []; % Clear the list
51 for i = 1 : nnode % Node cycle
52 eqnum = [ eqnum , lnods ( i ) *2 -1 , lnods ( i ) *2]; % Build the ←-
equation
53 end % number list
54
55 % Montar el vector de fuerza y la matriz de rigidez
56 for i = 1 : neleq
57 force ( eqnum ( i ) ) = force ( eqnum ( i ) ) + ElemFor ( i ) ;
58 for j = 1 : neleq
59 StifMat ( eqnum ( i ) , eqnum ( j ) ) = StifMat ( eqnum ( i ) , eqnum ( j ) ) + ...
60 ElemMat (i , j ) ;
61 end
62 end
63
64 end
65
66 num = length ( StifMat )
67 square = zeros ( num , num ) ;
68 for i = 1: num
69 for j =1: num
70 square (i , j ) = StifMat (i , j ) ;
71 end
72 end
73 format
74 square
75
76
77 ttim = timing ( ' Time to assemble the global system ' , ttim ) ; ←-
% Reporting time
78
79 % Agregue fuerzas laterales al vector de fuerza
80 for i = 1 : size ( sideload ,1)
81 x = coordinates ( sideload (i ,1) ,:) - coordinates ( sideload (i ,2) ,:) ;
82 l = sqrt ( x * transpose ( x ) ) ; % Find the length of the side
83 ieqn = sideload (i ,1) *2; % Find eq . number for the ←-
first node
84 force ( ieqn -1) = force ( ieqn -1) + l * sideload (i ,3) /2; % Add x ←-
force
85 force ( ieqn ) = force ( ieqn ) + l * sideload (i ,4) /2; % Add y ←-
force
86
87 ieqn = sideload (i ,2) *2; % Find eq . number for the ←-
second node
88 force ( ieqn -1) = force ( ieqn -1) + l * sideload (i ,3) /2; % Add x ←-
force
89 force ( ieqn ) = force ( ieqn ) + l * sideload (i ,4) /2; % Add y ←-
force
90 end
91
92
93
94 % Agregue condiciones de carga de punto al vector de fuerza
95 for i = 1 : size ( pointload ,1)
96 ieqn = ( pointload (i ,1) -1) *2 + pointload (i ,2) ; % Find eq . ←-
number
97 force ( ieqn ) = force ( ieqn ) + pointload (i ,3) ; % and add the ←-
force
98 end
99
100 ttim = timing ( ' Time for apply side and point load ' , ttim ) ; ←-
% Reporting time
101
102 force
103
104 num = length ( force )
105 square = zeros ( num , num ) ;
106 for i = 1: num
107 for j =1: num
108 square (i , j ) = force ( i ) ;
109 end
110 end
111 format
112 square
113
114 % Aplicar las condiciones de Dirichlet y ajustar el lado derecho
115 for i = 1 : size ( fixnodes ,1)
116 ieqn = ( fixnodes (i ,1) -1) *2 + fixnodes (i ,2) ; % Find the equation ←-
number
117 u ( ieqn ) = fixnodes (i ,3) ; % and store the ←-
solution in u
118 fix ( i ) = ieqn ; % and mark the eq . as a fix ←-
value
119 end
120
121 force1 = force - StifMat * u ; % Adjust the rhs with the ←-
known values
122
123
124
125 % Calcula la solucion resolviendo StifMat * u = fuerza para el resto
126 % de valores desconocidos de u
127 FreeNodes = setdiff ( 1: nndof , fix ) ; % Find the free ←-
node list
128 % and solve for it
129 u ( FreeNodes ) = StifMat ( FreeNodes , FreeNodes ) \ force1 ( FreeNodes ) ;
130
131
132 num = length ( u )
133
134 for i = 1: num
135 for j =1: num
136 square_U ( i ) = u ( i ) ;
137 end
138 end
139 format
140 square_U
141 % Calcule las reacciones en los nodos fijos como R = StifMat * u ←-
- F
142 reaction ( fix ) = StifMat ( fix ,1: nndof ) * u (1: nndof ) - force ( fix ) ;
143
144 ttim = timing ( ' Time to solve the stiffness matrix ' , ttim ) ; ←-
% Reporting time
145
146 % Calcule las tensiones
147 Strnod = Stress_v1_3 ( dmat , poiss , thick , pstrs , u )
148
149 ttim = timing ( ' Time to solve the nodal stresses ' , ttim ) ; ←-
% Reporting time
150
151 % Representacion grafica
152
153 ToGiD_v1_3 ( file_name ,u , reaction , Strnod ) ;
154
155 ttim = timing ( ' Time used to write the solution ' , ttim ) ; ←-
% Reporting time
156 itim = toc ; % Close ←-

last tic
157 fprintf (1 , '\ nTotal running time %12.6 f \ n \ n ' , ttim ) ; % Reporting ←-
final time
4.1.2 Función TrStif v1.3
Código Fuente 4.2: TrStif v1.3

1 function [M , F ] = TrStif_v1_3 ( nodes , dmat , thick , denss )
2
3 % % TrStif Evaluates the stiffness matrix for a triangular element
4 %
5 % Parameters :
6 %
7 % Input , nodes : Contains the 2 D coordinates of the element nodes
8 % dmat : Constitutive matrix
9 % thick : Thickness
10 % denss : Density
11 %
12 % Output , M the element local stiffness matrix
13 % F the element local force vector
14
15 b (1) = nodes (2 ,2) - nodes (3 ,2) ;
16 b (2) = nodes (3 ,2) - nodes (1 ,2) ;
17 b (3) = nodes (1 ,2) - nodes (2 ,2) ;
18
19 c (1) = nodes (3 ,1) - nodes (2 ,1) ;
20 c (2) = nodes (1 ,1) - nodes (3 ,1) ;
21 c (3) = nodes (2 ,1) - nodes (1 ,1) ;
22
23 area2 = abs ( b (1) * c (2) - b (2) * c (1) ) ;
24 area = area2 / 2;
25
26 bmat = [ b (1) , 0 ,b (2) , 0 ,b (3) , 0 ;
27 0 ,c (1) , 0 ,c (2) , 0 ,c (3) ;
28 c (1) ,b (1) ,c (2) ,b (2) ,c (3) ,b (3) ];
29
30 bmat = bmat / area2 ;
31
32 M = ( transpose ( bmat ) * dmat * bmat ) * area * thick ;
33
34 force = area * denss * thick /3;
35
36 F = [ 0 , - force , 0 , - force , 0 , - force ];
4.1.3 Función TrStrs v1.3
Código Fuente 4.3: TrStrs v1.3

1 function S = TrStrs_v1_3 ( nodes , dmat , displ , poiss , thick , pstrs )
2
3 % % TrStrs Evaluates the stresses for a triangular element
4 %
5 % Parameters :
6 %
9 % displ : Nodal displacement
10 % poiss : Poisson ratio
12 % pstrs : Flag for Plane Stress
13 %
14 % Output , S the element constant stress vector
15
16
17 b (1) = nodes (2 ,2) - nodes (3 ,2) ;
18 b (2) = nodes (3 ,2) - nodes (1 ,2) ;
19 b (3) = nodes (1 ,2) - nodes (2 ,2) ;
20
21 c (1) = nodes (3 ,1) - nodes (2 ,1) ;
22 c (2) = nodes (1 ,1) - nodes (3 ,1) ;
23 c (3) = nodes (2 ,1) - nodes (1 ,1) ;
24
25 area2 = abs ( b (1) * c (2) - b (2) * c (1) ) ;
26 area = area2 / 2;
27
28 bmat = [ b (1) , 0 ,b (2) , 0 ,b (3) , 0 ;
29 0 ,c (1) , 0 ,c (2) , 0 ,c (3) ;
30 c (1) ,b (1) ,c (2) ,b (2) ,c (3) ,b (3) ];
31
32 se = ( dmat * bmat * displ ) / area2 ;
33
34 % Plane Stress
35 if ( pstrs == 1)
36 S = se ;
37 % Plane Strain
38 else
39 S = [ se (1) , se (2) ,- poiss *( se (1) + se (2) ) , se (3) ];
40 end
4.1.4 Función QdStif v1.3
Código Fuente 4.4: QdStif v1.3

1 function [M , F ] = QdStif_v1_3 ( nodes , dmat , thick , denss )
2
3 % % QdStif Evaluates the stiffness matrix for a quadrilateral element
4 %
5 % Parameters :
6 %
10 % denss : Density
11 %
12 % Output , M the element local stiffness matrix

13 % F the element local force vector
14
15 fform = ←-
@ (s , t ) [(1 - s - t + s * t ) /4 ,(1+ s -t - s * t ) /4 ,(1+ s + t + s * t ) /4 ,(1 - s +t - s * t ) /4];
16 deriv = @ (s , t ) [( -1+ t ) /4 ,( 1 - t ) /4 ,( 1+ t ) /4 ,( -1 - t ) /4;
17 ( -1+ s ) /4 ,( -1 - s ) /4 ,( 1+ s ) /4 ,( 1 - s ) /4];
18
19 pospg = [ -0.577350269189626 E +00 , 0.57 73502 691896 26 E +00 ];
20 pespg = [ 1.0 E +00 , 1.0 E +00 ];
21 M = zeros (8 ,8) ;
22 fy = zeros (1 ,4) ;
23
24 for i = 1 : 2
25 for j = 1 : 2
26 lcffm = fform ( pospg ( i ) , pospg ( j ) ) ; % SF at Gauss point
27 lcder = deriv ( pospg ( i ) , pospg ( j ) ) ; % SF Local derivatives
28 xjacm = lcder * nodes ; % Jacobian matrix
29 ctder = xjacm \ lcder ; % SF Cartesian ←-
derivates
30 darea = det ( xjacm ) * pespg ( i ) * pespg ( j ) * thick ;
31
32 bmat = [];
33 for inode = 1 : 4
34 bmat = [ bmat , [ ctder (1 , inode ) , 0 ;
35 0 , ctder (2 , inode ) ;
36 ctder (2 , inode ) , ctder (1 , inode ) ] ];
37 end
38
39 M = M + ( transpose ( bmat ) * dmat * bmat ) * darea ;
40
41 fy = fy + lcffm * denss * darea ;
42
43 end
44 end
45
46 F = [ 0 , - fy (1) , 0 , - fy (2) , 0 , - fy (3) , 0 , - fy (4) ];
4.1.5 Función QdStrs v1.3
Código Fuente 4.5: QdStrs v1.3

1 function S = QdStrs_v1_3 ( nodes , dmat , displ , poiss , thick , pstrs )
2
3 % % QdStrs Evaluates the stresses for a quadrilateral element
4 %
5 % Parameters :
6 %
9 % displ : Nodal displacement
12 % pstrs : Flag for Plane Stress

13 %
14 % Output , S the element constant stress vector
15
16 fform = ←-
@ (s , t ) [(1 - s - t + s * t ) /4 ,(1+ s -t - s * t ) /4 ,(1+ s + t + s * t ) /4 ,(1 - s +t - s * t ) /4];
17 deriv = @ (s , t ) [( -1+ t ) /4 ,( 1 - t ) /4 ,( 1+ t ) /4 ,( -1 - t ) /4;
18 ( -1+ s ) /4 ,( -1 - s ) /4 ,( 1+ s ) /4 ,( 1 - s ) /4];
19
20 pospg = [ -0.577350269189626 E +00 , 0.57 73502 691896 26 E +00 ];
21 pespg = [ 1.0 E +00 , 1.0 E +00 ];
22
23 strsg = [];
24 extrap = [];
25 order = [ 1 , 4 ; 2 , 3 ]; % Align the G - pts . with the element ←-
corners
26
27 for i = 1 : 2
28 for j = 1 : 2
29 lcder = deriv ( pospg ( i ) , pospg ( j ) ) ; % SF Local derivatives
30 xjacm = lcder * nodes ; % Jacobian matrix
31 ctder = xjacm \ lcder ; % SF Cartesian ←-
derivates
32
33 bmat = [];
34 for inode = 1 : 4
35 bmat = [ bmat , [ ctder (1 , inode ) , 0 ;
36 0 , ctder (2 , inode ) ;
37 ctder (2 , inode ) , ctder (1 , inode ) ] ];
38 end
39
40 strsg (: , order (i , j ) ) = ( dmat * bmat * displ ) ;
41
42 a = 1/ pospg ( i ) ;
43 b = 1/ pospg ( j ) ;
44
45 extrap ( order (i , j ) ,:) = fform (a , b ) ;
46 end
47 end
48
49 se = transpose ( extrap * transpose ( strsg ) ) ;
50
51 % Plane Stress
52 if ( pstrs == 1)
53 S = se ;
54 % Plane Strain
55 else
56 S = [ se (1 ,:) ; se (2 ,:) ; - poiss *( se (1 ,:) + se (2 ,:) ) ; se (3 ,:) ];
57 end
4.1.6 Función constt
Código Fuente 4.6: constt
1 function D = constt ( young , poiss , pstrs )

2
3 % % constt Evaluates the constitutive matrix
4 %
5 % Parameters :
6 %
7 % Input , young : Young modulus
9 % pstrs : flag for Plane Stress
10 %
11 % Output , D the element constitutive matrix
12
13 % Plane Sress
14 if ( pstrs ==1)
15 aux1 = young /(1 - poiss ^2) ;
16 aux2 = poiss * aux1 ;
17 aux3 = young /2/(1+ poiss ) ;
18 % Plane Strain
19 else
20 aux1 = young *(1 - poiss ) /(1+ poiss ) /(1 -2* poiss ) ;
21 aux2 = aux1 * poiss /(1 - poiss ) ;
22 aux3 = young /2/(1+ poiss ) ;
23 thick = 1.0;
24 end
25
26 D = [ aux1 , aux2 ,0; aux2 , aux1 ,0;0 ,0 , aux3 ];
CONCLUSIÓN
Se logró estudiar y analizar los principios que abarca el MÉTODO AVANZADO DE LOS ELEMENTOS
FINITOS ; y así posteriormente llegar a determinar y editar el código fuente del programa Cálculo de
tensiones y deformaciones de vigas bidimensionales en voladizo con cargas puntuales o distribuidas ,
el cual fue de uso muy productivo tanto en nuestro curso de Métodos Numéricos , por abarcar nuestro
estudio en un método avanzando ; y así también en el curso de Resistencia de Materiales , por enfocar
el objetivo del presente trabajo en la determinación de esfuerzos y deformaciones de una estructura en
voladizo.
Hablando sobre el método de elementos finitos podremos concluir y mencionar que es una herramienta
que brinda muchas ventajas al momento que uno quiere analizar o diseñar una determinada pieza o
elemento estructural, pues gracias a su alta eficacia es posible detectar problemas (si estamos diseñando),
y en nuestro caso obtener precisos resultados de deformación y esfuerzo en una estructura, cuando la
sometemos a una determinada carga; posteriormente estos resultados serán considerados para realizar
un óptimo diseño, que es uno de nuestros objetivos como futuros Ingenieros Civiles.
BIBLIOGRAFÍA
• https://www.gidhome.com/
• https://la.mathworks.com/
• https://www.educacioneningenieria.org/index.php/edi/article/viewFile/242/152
• http://bdigital.unal.edu.co/10002/6/958932276X.2002.pdf
• https://www.youtube.com/watch?v=lPSoiR3Snv8
• https://ingemecanica.com/tutoriales/prontuario.html
• https://www.tdx.cat/bitstream/handle/10803/6294/06Efv06de23.pdf;...
jsessionid=DFB536E649B2FC0F0DBD509204B5529A?sequence=6

Trabajo Final Elementos Finitos

Cargado por

Copyright:

Formatos disponibles

Trabajo Final Elementos Finitos

Cargado por

Información del documento

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Trabajo Final Elementos Finitos

Cargado por

Copyright:

Formatos disponibles

Powered by TCPDF (www.tcpdf.

MÉTODOS NUMÉRICOS (IC-343)

INFORME SOBRE EL PROGRAMA Cálculos por

Ing. CASTRO PEREZ Cristhian

ACHALMA GALINDO, Melvin código: 16141315

Ayacucho, 30 de julio del 2018

1.1 ELASTICIDAD BIDIMENSIONAL 6

3 ANÁLISIS, PROCEDIMIENTO Y RESULTADOS 33

3.1 Ecuaciones Y Proceso De Uso 33

4.1 Código fuente 44

Figure 1: Viga en Voladizo.

ELEMENTOS FINITOS BIDIMENSIONALES

En ambos casos la tensión y la deformación en la dirección z no contribuyen a la energía elástica del

1.4 MECANISMO Y CONFORMACIÓN DE FUNCIONAMIENTO

1.5 CARACTERÍSTICA DE TRABAJO-TRANSMICIÓN DE CARGAS

1.6 CÁLCULO DE LA RESISTENCIA/ROTURA DE UNA VIGA EN VOLADIZO

1.6.1 CÁLCULO DEL MOMENTO

1.6.2 CÁLCULO DEL MOMENTO DE INERCIA

1.7 VENTAJAS DE LA ESTRUCTURA

1.8 DESVENTAJAS DE LA ESTRUCTURA

1.9 MEJORAS DE LA ESTRUCTURA

1.9.2 CONTROL DE LAS FUERZAS DE REACCIÓN EN SOPORTES

1.9.3 SISTEMA POR ZETEO

1.9.4 SISTEMA SUPLEX

1.10 ANÁLISIS ESTRUCTURAL

1.10.1 FUERZAS DE DISEÑO

1.10.2 DEFORMACIÓN EN VIGAS

1.10.3 MÉTODOS DE CÁLCULO

MÉTODO DE ÁREA DEL MOMENTO

MÉTODO DE DOBLE INTEGRACIÓN

1.11 Ejemplo de MEF

Suponga que para resolver la ecuación de Poisson, ∇2 u = p( x, y),en un medio bidimensional se

Empleando Elemento Finito Triangular (CST)

i, j, m −→ son permutaciones cíclicas de los índices 1,2,3

Cálculo para los elementos 1,3 y 5

Cálculo para los elementos 2,4 y 6

Ensamblando las matrices con los correspondientes Nudos

Finalmente la Ecuación en el Nudo E:

Cálculo para los elementos 1,5 y 8

Cálculo para los elementos 4 y 7

Utilizando la hoja de cálculo proporcionada en clase, se tienen los resultados:

ELEMENTOS FINITOS BIDIMENSIONALES

En el presente trabajo realizaremos el proceso de discretización mediante el programa GID, el cual se

MANEJO DEL SOFTWARE GID

2. Creamos la superficie, ya que anteriormente solo eran lineas unidas.

3. Cargamos el programa "Mat-femv1.3"

4. Dibujamos los apoyos.

5. Dibujamos la carga distribuida.

6. Introducimos las propiedades de los elementos.

7. Ponemos las unidades para los esfuerzos y deformaciones.

8. Discretizamos la superficie de la viga.

ANÁLISIS, PROCEDIMIENTO Y RESULTADOS

3.1 Ecuaciones Y Proceso De Uso

3.1.2 Funciones de Interpolacion

3.1.3 Coordenadas de Elemento

Z ( x, y) = N1 ( x, y)z1 + N2 ( x, y)z2 + N3 ( x, y)z3

3.1.4 La Matriz de Funciones de Forma

3.1.5 La Matriz de Deformacion Unitaria

3.1.6 La Matriz de Rigidez

3.2 Formulación del problema