Introducción al método de LOS ELEMENTOS FINITOS...
Transcript of Introducción al método de LOS ELEMENTOS FINITOS...
Publicaciones científicasUniversidad de las Fuerzas Armadas ESPE
José Fernando, Olmedo Salazar
Introducción al método de
APLICANDO MATHCADLOS ELEMENTOS FINITOS
campo unidimensional
Introducción al método de
APLICANDO MATHCADLOS ELEMENTOS FINITOS
campo unidimensional
Créditos
Campo unidimensionalFernando Olmedo
ISBN:978-9942-765-13-0
Marcelo TulioRommy Pérez
Universidad de las Fuerzas Armadas ESPECrnl. Ramiro Pazmiño (Rector)
Publicación autorizada por:Comisión Editorial de la Universidad de las Fuerzas Armadas ESPE
Edición y producción:David Andrade [email protected]
Diseño:David Cabrera [email protected]
Derechos reservados. Se prohíbe la reproducción de esta obra por cualquier
exclusiva responsabilidad del autor.
Universidad de las Fuerzas Armadas ESPEAv. General Rumiñahui s/n, Sangolquí, Ecuadorwww.espe.edu.ec
Los derechos de esta edición electrónica son de la Universidad de las Fuerzas Armadas ESPE, para consulta de profesores y estudiantes de la
universidad e investigadores en www.repositorio.espe.edu.ec.
Indice
Cápitulo IGeneralidadesIntroduccion1.3.- Problemás típicos1.4.- Conceptos1.5.- Procedimiento1.6.- Notación matricial1.7. - Simplificación1.8.-Aproximación directa en base de las ecuaciones de equilibrio1.9.- Ejemplos de ensamble1.10. - Ejemplos de aplicación1.11.- Aproximación mediante energía potencial1.12.- Solución mediante MathCAD1.13.- Ejercicios propuestos
Capítulo II Armaduras Planas
2.1.- Definición2.2.- Determinar la matriz de rigidez del elemento barra2.3.- Ecuaciones de transformación2.4.- Matriz global de rigidez2.5.-Ejemplos2.6.- Calculo de tensiones2.7.- Ejercicios2.8.- Resolución con Matlab según A.J.M. FERREIRA2.9.- Resolución con ANSYS APDL
Capítulo IIIAnálisis por el método de elementos finitos de elementos a flexión. Vigas
3.1.- Determinar la matriz de rigidez del elemento viga3.2.- Ejemplos cargas puntuales3.3.- Ejemplos cargas puntuales longitud y áreas de inercia variables3.4.- Problemas propuestos3.5.- Carga distribuida3.6.- Carga distribuida triangular3.7.-Calculo de vigas con elementos de alto orden (hP-FEM3.8.- Viga de Timoshenko3.9.- Resolución con ANSYS APDL
Capítulo IVAnálisis por el método de elementos finitos, Pórticos planos
4.1.- Determinar la matriz de rigidez del pórtico4.2.- Ejercicios con MATHCAD4.3.- Resolución de pórticos con ANSYS APDL
Pag.
Pag.Lista de Referencias
Lista de tablas
Tabla 2.1. Coordenadas en los nodos
Tabla 3.1. Coordenadas keypoints
Bibliografia
Acerca del autor
Dedicatoria
Dedicamos este manual a todos los estudiantes del DECEM que recibieron o recibirán la asignatura de elementos finitos aplicados.
6
Agradecimientos
Agradezco al Dr. David Andrade Aguirre por su ayuda e interés, manifestado en que esta obra se culmine, como también al Dr. Marcelo Tulio Piovan por sus consejos para mejorar la misma y por sus cono-cimientos entregados sobre el método de elementos finitos, en el curso impartido por él, durante su estancia como Prometeo en el Ecuador. También agradezco al Ing. Romy Pérez Moreno por su tiempo de-dicado a la revisión, como también al director del DECEM Ing. Carlos Naranjo G. y a todos los compañeros del departamento.
7
Resumen
El método de elementos finitos es hoy por hoy una técnica suma-mente extendida a una innumerable cantidad de aplicaciones ingenie-riles, aplicar los principios a través de programación es trascendental para entender el método. Este es un primer manual que trata de ex-poner en una forma accesible y pedagógica las complejidades inheren-tes del método de elementos finitos, considerando que el estudiante que toma el curso de elementos finitos, es un estudiante de último ni-vel, habido de poner en práctica los conocimientos teóricos impartidos y deseoso de experimentar con retos ingenieriles más que lidiar con matemática. Para el presente trabajo nos inspiramos en algunas obras recientes como por ejemplo “Finite Element Modeling and Simulation with ANSYS Workbench de Xialin Chen”, que hacen un tránsito por el método de los elementos finitos conjugando el desarrollo teórico con aplicaciones de uno de los programas más potentes del mercado. Exis-ten muchos textos que abordan el método con aplicaciones de Matlab, debido a que en el Departamento de Energía y Mecánica se trabaja en las asignaturas de mecanismos y vibraciones con MathCAD y puesto que los estudiantes están familiarizados con este programa, se lo aplico con bastante éxito desde el punto de vista pedagógico que es el aspecto que nos motivaba principalmente.
Trabajar con el MathCAD significo priorizar la comprensión del método sobre la efectividad y rapidez de calculo que podría ofrecer el Matlab, sin embargo en el capítulo I, únicamente para establecer una comparativa, se aplica para el cálculo de armaduras el programa del profesor A. J. M. Ferreira. El presente libro de texto cubre solamente la parte unidimensional a través de la aproximación directa, es decir sis-temas de resortes, armaduras, vigas y pórticos planos, en una futura edición o volumen se abordara el campo bidimensional.
11
Introducción al método de los elementos finitos aplicando MATHCAD
Objetivo de la obra:
El presente texto fue realizado para brindar a los estudiantes el adecuado apoyo en la asignatura de Elementos Finitos Aplicados, transparentando los principios que existen detrás de los programas comerciales que son tan extendidos en la actualidad y de uso cotidiano en las oficinas de diseño. Básicamente se han buscado tres objetivos:1. Que el estudiante entienda y aplique los conceptos fundamentales
del modelamiento por elementos finitos para resolver e implemen-tar soluciones en problemas simples.
2. Una vez alcanzado esto, obtener la competencia en el manejo de un software de propósito general como ANSYS APDL sin perder de vista la teoría.
3. Modelar problemas ingenieriles de forma profesional con el soft-ware de elementos finitos ANSYS APDL y ANSYS WORKBENCH, licencias adquiridas por la ESPE en el año 2012.
Introducción:
El método de Análisis por Elementos Finitos (comúnmente llamado FEA, del inglés Finite Element Analysis), tiene su génesis en el diseño estructural y fue introducido por Argyris, Turner, Clough y Zienkiewicz1
Esencialmente es un método numérico que genera soluciones de tipo aproximado para cualquier tipo de problemas de la ingeniería, in-calculables con métodos matemáticos convencionales.
El método de elementos finitos “FEM” (del inglés Finite Element Method) ha llegado a ser un paso esencial en el diseño y modelado en varias disciplinas de ingeniería.
La base del FEM se establece en la descomposición del dominio en un número finito de subdominios (elementos) a los cuales se aplica leyes constitutivas y se crean un sistema de ecuaciones algebraicas simultáneas por medio de aplicar las siguientes aproximaciones prin-cipales para construir una solución basada en el concepto de FEM:
• Aproximación directa: esta estrategia se utiliza en problemas rela-tivamente simples, y sirve generalmente como medio para explicar el concepto de FEM y sus pasos importantes.
<?> A brief history of the beginning of the finite element method, http://ed.Iitm.Ac.In/~palramu/ed403_2012/files/fehis-tory_gupta.Pdf
12
Generalidades
• Residuos ponderados: este es un método versátil, permitiendo el uso de FEM a los problemas cuyo funcional (energía potencial) no puede ser construido. Esta aproximación utiliza directamente las ecuaciones diferenciales de gobierno, tales como trasferencia de calor, mecánica de fluidos y torsión.
• Aproximación variacional: este acercamiento confía en el cálculo de variaciones, que implica el extremar un funcional. Este funcio-nal corresponde a la representación energética del sistema (si es una estructura será energía potencial, si se trata de un problema de temperatura, Energía calórica, etc)
Problemás típicos
El método de elementos finitos “FEM” permite resolver problemas de • Análisis estructural• Transferencia de calor• Fluidos• Acústica• Transporte de masa• Potencial electromagnético
Los siguientes son algunos ejemplos:
Figura 1.1. Visualización de esfuerzos de Von Mises en un modelo Formula SAE Student.
Fuente propia
13
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 1.2. Visualización de la velocidad en una camioneta.
Fuente propia
Figura 1.3. Determinación de la estabilidad de un cilindro.
Fuente propia
14
Generalidades
Figura 1.4. Campo Magnético con Flexpde.
Fuente propia
Conceptos
El método de elementos finitos aborda el problema ingenieril median-te la división de un dominio complejo en elementos no intersecantes y expresa la variable de campo desconocida (desplazamientos, temperatura, velocidad) en términos de una función arbitraria de aproximación (aproxi-mar la solución), que se asume dentro de cada elemento. Estas variables se reemplazan en la ecuación constitutiva del problema físico que luego será incorporada a una ecuación de energía potencial que se minimiza mediante derivación parcial para obtener la ecuación de equilibrio y la matriz de rigidez de un elemento. La Fig. 1.5 es un mapa conceptual que ayuda a entender el método de elementos finitos.
Esta función de aproximación (también llamada función de inter-polación) es definida en términos de los valores de la variable de campo en puntos específicos referidos como nodos. Los nodos son usualmente localizados a lo largo de los límites de los elementos y ellos conectan elementos adyacentes.
La característica que ha hecho al método de los elementos finitos tan popular es que puede ser dividido en un conjunto de pasos lógicos y puede ser usado para analizar un amplio rango de problemas solo cambiando los datos de entrada del dominio. Estos pasos se explicarán en el próximo apartado.
15
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 1.5. Mapa Conceptual
Fuente propia
Procedimiento
El método de los elementos finitos requiere los siguientes pasos principales:
• Discretización del dominio en un número finito de subdominios (elementos), los elementos deben ser lo suficientemente pequeños para dar resultados correctos y lo suficientemente grande para re-ducir el esfuerzo computacional. Estos elementos son conectados el uno al otro por sus nodos comunes (Fig. 1.6). Un nodo especifica la localización coordinada en el espacio donde existen los grados de libertad (desplazamientos) y las acciones del problema físico. Las incógnitas nodales en el sistema de la matriz de ecuaciones (desplazamientos) representan una (o más) de las variables primarias del campo. Las variables nodales asignadas a un elemento se llaman los grados de libertad del elemento. Dependiendo del problema se utilizan los siguientes tipos de ele-mentos. (Fig. 1.7).
16
Generalidades
Figura 1.6. Discretización de dominio y de fuerzas
Fuente propia
Figura 1.7. Discretización de dominio y de fuerzas
Fuente propia
17
Introducción al método de los elementos finitos aplicando MATHCAD
• Selección de las funciones de desplazamiento. Para un ele-mento lineal es función de una magnitud y si es bidimensional es función de x, y. El campo de desplazamiento desconocido dentro de un elemento finito puede ser interpolado por una distribución aproximada. Esta distribución se vuelve más exacta conforme se consideran más elementos en el modelo y pueden ser expresadas como funciones polinomiales que pueden ser fácilmente derivadas e integradas. Por ejemplo, para el caso más simple de elemento fi-nito que es el elemento barra, los desplazamientos son expresados en términos de los desplazamientos nodales {u1, u2} por medio del polinomio de primer grado.
(1.1)
Donde l es la longitud del elemento finito, x es la variable indepen-diente y las funciones de interpolación N1 y N2 son:
(1.2)La figura 1.8 representa la interpolación lineal del campo de des-
plazamiento dentro del elemento barra:
Figura 1.8. Representación de la interpolación lineal
Fuente propia
18
Generalidades
• Desarrollo de la matriz de elementos para el subdominio, utilizan-do las ecuaciones constitutivas
• Ensamblaje de la matriz de elementos para cada subdominio para obtener la matriz global del dominio entero
• Imposición de las condiciones de frontera• Solución de ecuaciones • Postprocesado
Notación matricial
En notación matricial, el sistema de ecuaciones global puede ser definido como
(1.3)
Donde la matriz K representa la matriz de rigidez del sistema, u es el vector del campo desconocido, y F es el vector de la fuerza. Dependiendo de la naturaleza del problema, K puede ser dependiente en u, es decir, K = K (u) y F pueden ser dependientes del tiempo, es decir, F = F (t).
Generalizando se obtiene:
(1.4)
Si asumimos que el desplazamiento de u1 es 1 y de v1 hasta wn es 0, se puede resolver fácilmente el sistema
F1x = K11 F1y = K21;. . .; Fnz = Kn1
19
Introducción al método de los elementos finitos aplicando MATHCAD
Simplificación
Figura 1.9. Simplificaciones geométricas
Fuente propia
La competencia en la teoría de elementos finitos proporciona crite-rios para simplificar las geometrías para optimizar recursos computa-cionales sin perder precisión.
Aproximación directa en base de las ecuaciones de equilibrio
Determinar la matriz de rigidez del elemento resorte. Usando la aproximación directa en base de las ecuaciones de equi-
librio, se define la matriz de rigidez para un resorte lineal de una di-mensión, el resorte Fíg.1.10 obedece la ley de Hooke y resiste fuerzas únicamente en la dirección axial.
El resorte de la figura está sometido a las fuerzas locales: f1x, f2x y sufre los desplazamientos locales: u1, u2 por lo tanto existen 2 grados de libertad
20
Generalidades
Figura 1.10. Resorte dos grados de libertad
Fuente propia
1. Selección del tipo de elementoElemento resorte lineal con 2 grados de libertad, sujeto a tensio-
nes nodales T y con longitud L
Figura 1.11. Resortes tensionadosFuente propia
Como se puede observar un extremo se deforma U2 y el otro U1 donde U2<U1 y por tanto la deformación relativa total es U2-U1
2. Seleccionar una función de desplazamiento Se debe escoger la función matemática que represente el des-
plazamiento bajo carga del resorte, se selecciona en forma arbitraria una función lineal de 2 constantes, ya que el elemento tiene 2 gra-dos de libertad:
u = a1 + a2 x, en forma matricial se obtiene:
(1.5)
Para determinar u se debe hallar las constantes a1 y a2 utilizando las condiciones de frontera u1, u2
21
Introducción al método de los elementos finitos aplicando MATHCAD
Reordenando se tiene
En forma matricial
(1.6)
Donde N1 y N2 son las funciones de forma o interpolación
Figura 1.12. Grafica de la función de interpolación
Fuente propia
Cuando x= 0 u= u1; u1 = a1 + a2 (0); por tanto a1 = u1
Cuando x = L u= u2; u2 = u1 + a2 (L); por tanto a2 = (u2-u1)/L
22
Generalidades
3. Definir la relación tensión / deformación La fuerza tensora T produce una deformación δ en el resorte:
δ = u (L) – u (0) = u2 – u1
T = k δ = k (u2 – u1) (1.7)
4. Determinar la matriz de rigidez
f1x = - T = k (u1 – u2)
f2x = T = k (-u1 + u2) (1.8)
En forma matricial
(1.9)
La matriz de rigidez local es por tanto:
(1.10)
5. Ensamblar las ecuaciones de los elementos para obtener las ecuaciones globales
(1.11)
6. Resuelva para los desplazamientos nodales
(1.12)
Los desplazamientos son determinados imponiéndose condiciones de frontera
7. Resuelva las fuerzas de los elementos Las fuerzas son determinadas por sustitución inversa aplicando
nuevamente (1.12)
23
Introducción al método de los elementos finitos aplicando MATHCAD
Ejemplos de ensamble
Ensamblaje por medio de ecuaciones de equilibrioLos pórticos, puentes y otras estructuras se componen de com-
ponentes estructurales básicos, el primer método para determinar la matriz de rigidez global es mediante el desarrollo de las ecuaciones de equilibrio a través de generar el diagrama de cuerpo libre de cada ele-mento2
Para el elemento 1 y usando (1.9) se tiene:
Donde los superíndices se refieren al elemento
Para el elemento 2:
Debido a que el nodo 3 es común se establece la siguiente ecuación de compatibilidad
Basado en las ecuaciones de equilibrio de los nodos establecemos las siguientes ecuaciones
Fuerza global en 1 es igual a fuerza local del nodo 1
2 Daryl L. Logan, “A First Course in the Finite Element Method”, Pág. 35
24
Generalidades
Por lo queF1x = k1 u1 – k1 u3 F2x = -k2 u3 + k2 u2 F3x = (- k1 u1 + k1 u3) + (k2 u3 – k2 u2)
Re ensamblando nuevamente la matriz se tiene:
Dónde:
Es el vector de fuerzas nodales globales
Es la matriz de rigidez del sistema
Es el vector de desplazamiento global
Ensamblaje por superposición
Como se puede inferir si es que se aumentan los elementos, el método anterior es impráctico. Un método más rápido es el método de superposición donde las ecuaciones locales matriciales se expanden de la siguiente manera3:
Para el elemento 1 se tiene:
Para el elemento 2 se tiene:
3 Daryl L. Logan, “A First Course in the Finite Element Method”, Pág. 38
25
Introducción al método de los elementos finitos aplicando MATHCAD
La sumatoria genera:
Ensamblaje por expansión directa
Un paso más en la automatización es el método de expansión directa donde las ecuaciones locales matriciales se expanden uti-lizando un identificador de acuerdo al nodo correspondiente, de la siguiente manera: Para el elemento 1 Para el elemento 2
Luego se genera la matriz de rigidez global llenando las filas y co-lumnas de acuerdo a lo establecido anteriormente. Lo que se ha hecho es formar fila y columna numerando los nodos y sumamos
Condiciones de frontera homogeneas
Sin condiciones apropiadas la estructura no resistiría las fuerzas aplicadas
Figura 1.13. Condición de frontera homogénea
Daryl L. Logan, “A First Course in the Finite Element Method”, Pág. 35
26
Generalidades
Debido a que el resorte esta empotrado en la pared se tiene que u1 = 0
Si se resuelve el sistema se obtiene el mismo resultado que elimi-nar la primera fila y la primera columna puesto que es cero el primer término del vector:
Substituyendo en la ecuación general se obtiene F1x, F2x y F3x
Condiciones de frontera no homogeneas
Uno o más desplazamientos no son cero, la tensión es debida al estiramiento que recibe el resorte para ensamblarlo
Figura 1.14. Condición de frontera no homogénea4
Daryl L. Logan, “A First Course in the Finite Element Method”
4 Daryl L. Logan, “A First Course in the Finite Element Method”, Pág. 46
27
Introducción al método de los elementos finitos aplicando MATHCAD
Se resuelve en forma convencional y puesto que F1x es desconocida
Ejemplos de aplicación1. Determine la matriz de rigidez global, los desplazamientos de los
nodos 2 y 3, las reacciones en 1 y 4, y las fuerzas en cada resorte. Una fuerza de 5000 N es aplicada en el nodo 3 en la dirección x.
1. Definimos las matrices de rigidez individuales
2. Ensamblamos la matriz de rigidez global y establecemos la ecua-ción de elasticidad, 4 grados de libertad de los cuatro nodos la ma-triz es 4x4
La ecuación de rigidez es:
28
Generalidades
Debido a que los desplazamientos son cero se puede eliminar la primera fila, primera columna, así como la cuarta fila y cuarta columna
3. Determinamos los desplazamientos y las fuerzas
u2 = 0.909, u3 = 1.364
Mediante substitución inversa:F1x = -909.091
F2x = 0
F3x = 5000
F4x = -4091
4. Determinar las fuerzas de cada elemento usando las ecuaciones locales
Elemento 1:
Elemento 2:
Elemento 3:
29
Introducción al método de los elementos finitos aplicando MATHCAD
2. Para el ensamblaje de resortes obtenga la matriz de rigidez global, los desplazamientos de los nodos 2, 3, 4, las fuerzas noda-les globales y las fuerzas locales de cada resorte. El nodo 1 es fijo mientras que al nodo 5 se le da un desplazamiento conocido de 20 mm. Las constantes de los resortes son todas iguales a k = 200 kN/mm
La matriz de rigidez global es de 5 x 5 y relaciona fuerzas a despla-zamientos como:
0 0 0 200 200 0 0 02 0 0 200 400 200 0 0
0 2 0 0 200 400 200 00 0 2 0 0 200 400 2000 0 0 0 0 0 200 200
k kk k k
k k kk k k
k k
� �� � � �� � � �� � � �� � � �� � � �=� � � �� � � �
� � � �� � � �� � � �� �� � � �
1 200 200 0 0 0 12 200 400 200 0 0 23 0 200 400 200 0 34 0 0 200 400 200 45 0 0 0 200 200 5
F x uF x uF X uF X uF X u
�� � � �� �� � � �� �� �� � � �� �� � � �� �= � �� � � �
� �� � � �� �� �� � � �� ��� � � �� � � �� �
Utilizando el formato de solución numérica de MathCAD
30
Generalidades
Aproximación mediante energía potencial
Otro método utilizado para obtener la matriz de rigidez es el prin-cipio de energía potencial mínima. Este método general se utiliza para tensión y deformación plana, tensión axisimétrico, flexión en placas y tensión en elementos sólidos. Este método es aplicable únicamente a materiales con elasticidad lineal. El principio dice que: “De todas las formas geométricamente posibles que un cuerpo puede asumir, la que corresponde a la satisfacción de equilibrio estable del cuerpo, es identifi-cada por un valor mínimo de la energía potencial total.”
Por lo tanto si se obtiene la energía potencial total y la minimizamos con respecto a los desplazamientos obtendremos la matriz de rigidez.
La energía potencial total se define como la suma de la energía de deformación interna U (relacionado a fuerzas internas) y la energía potencial de las fuerzas externas Ω, es decir,
(1.13)
Usando integración
(1.14)
31
Introducción al método de los elementos finitos aplicando MATHCAD
La energía potencial de las fuerzas externas que es perdida es representada por:
(1.15)
Por lo tanto:
(1.16)
Para aplicar el principio de energía potencial mínima debemos minimizar la función
(1.17)
Ejemplo: Para el resorte dado evalúe la energía potencial en fun-ción de x y demuestre que la mínima energía potencial corresponde a la posición de equilibrio
K := 500 F := 1000 x := 50, 49.9..50 E x( )k x2
×
2F x�:=
Figura 1.15. Energía Potencial de un resorte
3� 2.1� 1.2� 0.3� 0.6 1.5 2.4 3.3 4.2 5.1 6
1.5� 103�1.15� 103�
800�450�100�250600950
1.3 103�1.65 103�
2 103�
Energía
E x( )
2
x
La energía potencial mínima es -1000 N mmUsando la expresión (1.17):
32
Generalidades
Donde x = 2 mm, este valor es reemplazado en la expresión de la energía (1.16)
Por lo tanto se demuestra que la energía potencial mínima corres-ponde con la posición de equilibrio.
Derivación de la matriz de rigidez en un resorte utilizando el principio de mínima energía potencial
(1.18)
La minimización de la energía potencial requiere tomar derivadas parciales con respecto a cada desplazamiento nodal
Simplificando se obtiene:
Y en forma matricial
33
Introducción al método de los elementos finitos aplicando MATHCAD
Solución mediante Mathcad
Para el ensamblaje de resortes obtenga la matriz de rigidez global, los desplazamientos de los nodos 2, las fuerzas nodales glo-bales y las fuerzas locales de cada resorte.
k1=1000 N/mm, k2=2000 N/mm, k3 = 3000 N/mm
Figura 1.15. Sistema de resortes en serie y en paralelo
Se ha propuesto el siguiente algoritmo de MathCAD para resol-ver los problemas de resortes (utilizar el programa resorte general que acompaña al texto).
ORIGIN := 1
Elementos:= 3 nodos := 4
conectividad :=
K k( ) k1
1
1
1
k
1000
2000
3000
K k1( )1 103
�
1� 103�
1� 103�
1 103�
����
�÷÷�
= K k2( )2 103
�
2� 103�
2� 103�
2 103�
����
�÷÷�
= K k3( )3 103
�
3� 103�
3� 103�
3 103�
����
�÷÷�
=
Knodos nodos, 0:=
34
Generalidades
K
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
������
�÷÷÷÷�
=
K
Kconectividadn i, conectividadn j, , Kconectividadn i, conectividadn j, ,
Kin( )i j,
+¬
j 1 2..�for
K
i 1 2..�for
K
n 1 elementos..�for
K
:=
K
1 103
1 103
0
0
1 103
6 103
2 103
3 103
0
2 103
2 103
0
0
3 103
0
3 103
u1 0:= u3 0:= u4 0:= f2 10:=
K
u1
u2
u3
u4
������
�÷÷÷÷�
×
f1
f2
f3
f4
������
�÷÷÷÷�
solve u2, f1, f3, f4, 1
60053
�103
� 5����
�÷�
®
La solución indica que el desplazamiento del nodo 2 es de 1/600 mm y las respec-tivas fuerzas internas son -5/3,-10/3 y -5
35
Introducción al método de los elementos finitos aplicando MATHCAD
Ejercicios propuestos
1. Para el ensamblaje de resortes obtenga la matriz de rigidez glo-bal. Resuelva con MathCAD, k = 100 N/mm
2. Resuelva tanto con MathCAD, k = 100 N/mm y F = 1000 N
39
Introducción al método de los elementos finitos aplicando MATHCAD
Definición: Una armadura es una construcción reticulada conformada gene-
ralmente por triángulos formados por elementos rectos y que se utiliza para soportar cargas. Las armaduras pueden ser planas o espaciales. El análisis de armaduras dio pie al desarrollo de la teoría de elementos finitos y su aplicación es amplia: torres, puentes, cubiertas se fabrican a través de esta conjunción de barras que solo se pueden cargar axial-mente. Ver Fig. 2.1
Figura 2.1. Algunos tipos de armaduras planas
http//es.slideshare.net/guest1f9b03a/cap6r
Determinar la matriz de rigidez del elemento barra:
Para determinar la ecuación de rigidez debemos analizar el ele-mento barra que tiene las siguientes características:
Figura 2.2. Elemento barra horizontal
Fuente propia
40
Armaduras planas
Fuerzas locales: f1x, f2x
Desplazamientos locales: u1, u2 por lo tanto existen 2 grados de libertadDe la ley de Hooke se obtiene las ecuaciones constitutivas:
(2.1)
(2.2)
De las condiciones de equilibrio se tiene:
(2.3)
Por lo tanto la derivada será cero:
(2.4)
Y está es la ecuación diferencial que gobierna el sistema y que pue-de ser utilizada en los métodos de residuos ponderados.
Los pasos a seguir para determinar la matriz de rigidez son:
1. Selección del tipo de elementoElemento barra con 2 grados de libertad, sujeto a tensiones noda-les T y con longitud L
2. Seleccionar una función de desplazamiento Se debe escoger la función matemática que represente la deforma-da bajo carga del resorte, se selecciona la función lineal de 2 cons-tantes porque se tiene dos grados de libertad:
Se representa u en función de las condiciones de frontera u1, u2
41
Introducción al método de los elementos finitos aplicando MATHCAD
En x=0 u=u1;por tanto u1=a1+a2(0); a1=u1En x=L u=u2; por tanto u2=u1+a2(L); a2=(u2-u1)/L
u=u1+(u2-u1)(x/L)
Figura 2.3. Grafico funciones de interpolación
Fuente propia
En forma matricial las funciones de interpolación se representan de la siguiente manera:
(2.5)
N1, N2 son funciones de forma o de interpolación
3. Definir la relación tensión / deformación La relación deformación unitaria vs desplazamiento es:
(2.6)
La relación deformación unitaria vs esfuerzo es:
σx=E ϵx (2.7)
42
Armaduras planas
4. Determinar la matriz de rigidez del elemento De igual manera que se realizó para el resorte
En forma matricial
La matriz de rigidez local es por tanto: (2.8)
5. Ensamblar las ecuaciones de los elementos para obtener las ecua-ciones globalesSe efectúan las sumatorias respectivas
(2.9)
6. Resuelva para los desplazamientos nodales
[F]=[K][d] (2.10)
Los desplazamientos son determinados imponiéndose condiciones de frontera
7. Resuelva las fuerzas de los elementos Las fuerzas son determinadas por sustitución inversa usando (2.10)
43
Introducción al método de los elementos finitos aplicando MATHCAD
Ejercicio de ensamble Para el ensamblaje siguiente determine a) La matriz de rigidez global, b) los desplazamientos de los nodos 2 y 3, c) Las reacciones en los nodos 1 y 4, Una fuerza de 15000 N es aplicada en el nodo 2, E = 200000 MPa y D = 25 mm para el elemento 1 y 3 y E = 100000 MPa y D = 50 mm para el elemento 2
ORIGIN := 1
Elementos:= 3
nodos := 4
conectividad := k
k1
k2
k3
2.945 105
5.89 105
2.945 105
K k( ) k1
1
1
1
n := 1.. elemento
K k1( )2.945 105
2.945 105
2.945 105
2.945 105
K k2( )
5.89 105
5.89 105
5.89 105
5.89 105
K k3( )2.945 105
2.945 105
2.945 105
2.945 105
44
Armaduras planas
Ki1 K k1
Ki2 K k2
Ki3 K k3
Knodos nodos0
K
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
K
Kconectividadn i conectividadn jKconectividadn i conectividadn j
Kin i j
j 1 2for
K
i 1 2for
K
n 1 elementosfor
K
K
2.945 105
2.945 105
0
0
2.945 105
8.836 105
5.89 105
0
0
5.89 105
8.836 105
2.945 105
0
0
2.945 105
2.945 105
A continuación colocamos los datos conocidos y resolvemos sim-bólicamente
u1 0 u4 0 f2 15000 f3 0
K
u1
u2
u3
u4
f1
f2
f3
f4
solve
f1
u2
u3
f4
9000.0 0.0305577490736439051010.020371832715762603401 6000.0( )
u2 = 0.03 mm u3 = 0.02 mmf1 = -9000 N f4 = -6000 N
45
Introducción al método de los elementos finitos aplicando MATHCAD
Ecuaciones de transformación:
Como se puede observar en la Fig. 2.1. Las barras de las armadu-ras están orientadas en el plano con un ángulo definido, es convenien-te por tanto, transformar coordenadas locales a lo largo de la barra a coordenadas globales referenciadas a un sistema de coordenadas abso-luto o viceversa. En este caso se van a relacionar los desplazamientos del sistema de coordenadas d global (x, y) a coordenadas locales (x´, y´), fig. 2.4.
Figura 2.4. Barra general
Fuente propia
Las proyecciones del vector en el nodo i con respecto al sistema global son:
Matricialmente
(2.11)
46
Armaduras planas
Por lo tanto si despejamos resolviendo el sistema de ecuaciones o invirtiendo la matriz se obtiene:
Donde son los desplazamientos locales y dx son los desplaza-mientos globales
En forma matricial
(2.12)
c
s
s
c
1
c
c2 s2
s
c2 s2
s
c2 s2
c
c2 s2
Relación de transformación para el vector de desplazamientosPartiendo de la relación demostrada en (2.12) y considerando am-
bos nodos en una barra, se debe expandir la matriz
(2.13)
47
Introducción al método de los elementos finitos aplicando MATHCAD
Relación de transformación para el vector de fuerzas:Las fuerzas se transforman de igual manera que los desplazamientos:
(2.14)
Matriz global de rigidezPartiendo de la matriz local de un elemento genérico, se ensambla-
ra la matriz global de la estructura. Para una barra en el sistema local de coordenadas se tiene:
(2.15)
Puesto que la barra solo soporta cargas axiales los componentes en y son nulos por tanto la matriz local de rigidez se expande también según:
(2.16)
Reemplazando (2.13) y (2.14) en (2.16)
O lo que es lo mismo
48
Armaduras planas
Despejando:
La matriz global de rigidez es por tanto:
Simétrica (2.17)Se ha obtenido la relación fuerzas nodales globales con los despla-
zamientos nodales globales
(2.18)Ejemplos
El desplazamiento global fue determinado en d2x = 10 mm y en d2y = 20mm. Determine el desplazamiento local d2x y d2y
49
Introducción al método de los elementos finitos aplicando MATHCAD
Determine la matriz de rigidez global con respecto a los ejes x & y de la barra inclinada
50
Armaduras planas
Calculo de tensiones:
En una barra los desplazamientos se relacionan con las fuerzas según la expresión conocida (2.8) de la cual se toma únicamente cualquiera de las dos componentes de la matriz de rigidez local, en este caso .
(2.19)
De la relación siguiente:
Tomando solo las componentes en x
(2.20)
La definición del esfuerzo está dada por:
(2.21)
Combinando (2.20) en (2.19)
(2.22)
Donde C´ es
(2.23)
51
Introducción al método de los elementos finitos aplicando MATHCAD
La operación se la efectúa mediante:
(2.24)
Ejercicios:
1. Para la barra indicada en la figura determine el esfuerzo axial, si A = 4 E-4 m2, E = 210 GPa., L = 2m, Los desplazamientos glo-bales fueron determinados d1x = 0.25 mm, d1y = 0.0 mm, d2x = 0.50 mm y d2y = 0.75 mm
52
Armaduras planas
2. Para la estructura plana indicada sujeta a una fuerza de 10000 N aplicada en el nodo 1 determine la matriz de rigidez, las de-formaciones, las tensiones
Figura 2.5. Armadura de tres barras
Ensamblaje ManualPara cada barra se debe calcular la matriz simétrica
Si son 4 nodos y dos grados de libertad por nodo, la matriz global es de 8x8 y por lo tanto debemos expandir las matrices individuales
d1x d1y d2x d2y d1x d1y d3x d3y
53
Introducción al método de los elementos finitos aplicando MATHCAD
d1x d1y d4x d4y
A continuación se procede a sumar las matrices:
Si la armadura está sujeta en los nodos 2,3 y 4 se tiene la siguiente ecuación de elasticidad:
Por lo que se pueden eliminar desde la segunda a la cuarta fila y columna mediante la sentencia:
Kr submatrixKtotal 0 1 0 1( )
Kr2.256 105
5.893 104
5.893 104
2.256 105
Los desplazamientos se los calcula utilizando la matriz inversa:
d Kr 1 0
10000
d0.012
0.048
54
Armaduras planas
Ensamblaje con MathCADEl siguiente es el algoritmo desarrollado en el archivo “armadura2.mcd”
Datos de elementos, nodos y conectividad
Grados Longitud
90
45
0
180
1.571
0.785
0
Lo
1
2
1
Matriz de rigidez
K ( )
cos ( )2
cos ( ) sin ( )
cos ( )2
cos ( ) sin ( )( )
cos ( ) sin ( )
sin ( )2
cos ( ) sin ( )( )
sin ( )2
cos ( )2
cos ( ) sin ( )( )
cos ( )2
cos ( ) sin ( )
cos ( ) sin ( )( )
sin ( )2
cos ( ) sin ( )
sin ( )2
n 1 elementos
K 1
0
0
0
0
0
1
0
1
0
0
0
0
0
1
0
1
K 2
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
K 3
1
0
1
0
0
0
0
0
1
0
1
0
0
0
0
0
55
Introducción al método de los elementos finitos aplicando MATHCAD
Matriz de cerosKt1 2 nodos 0
Kt
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
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
0
0
0
0
fuerza
0
10000
f2x
f2y
f3x
f3y
f4x
f4y
f2x
d
d1x
d1y
0
0
0
0
0
0
d1x
Algoritmo de adición matricialKt
p conectividadn 1
q conectividadn 2
r q p
Ktp p Ktp pK n 1 1
Lon
Ktp p 1 Ktp p 1
K n 1 2 Lon
Ktp q r Ktp q r
K n 1 3 Lon
Ktp q r 1 Ktp q r 1
K n 1 4 Lon
Ktp 1 p Ktp 1 pK n 2 1
Lon
Ktp 1 p 1 Ktp 1 p 1
K n 2 2 Lon
Ktp 1 q r Ktp 1 q r
K n 2 3 Lon
Ktp 1 q r 1 Ktp 1 q r 1
K n 2 4 Lon
Ktq r p Ktq r pK n 3 1
Lon
Ktq r p 1 Ktq r p 1
K n 3 2 Lon
Ktq r q r Ktq r q r
K n 3 3 Lon
Ktq r q r 1 Ktq r q r 1
K n 3 4 Lon
Ktq r 1 p Ktq r 1 pK n 4 1
Lon
Ktq r 1 p 1 Ktq r 1 p 1
K n 4 2 Lon
Ktq r 1 q r Ktq r 1 q r
K n 4 3 Lon
Ktq r 1 q r 1 Ktq r 1 q r 1
K n 4 4 Lon
Kt
n 1 elementosfor
Kt Kt EArea
Longitud
56
Armaduras planas
Se despliega la matriz global
Kt
2.256 105
5.893 104
0
1.021 10 11
5.893 104
5.893 104
1.667 105
0
5.893 104
2.256 105
1.021 10 11
1.667 105
5.893 104
5.893 104
0
0
0
1.021 10 11
0
1.021 10 11
0
0
0
0
1.021 10 11
1.667 105
1.021 10 11
1.667 105
0
0
0
0
5.893 104
5.893 104
0
0
5.893 104
5.893 104
0
0
5.893 104
5.893 104
0
0
5.893 104
5.893 104
0
0
1.667 105
0
0
0
0
0
1.667 105
0
0
0
0
0
0
0
0
0
Se plantea solución simbólica mediante
Kt d fuerza solve
d1x
d1y
f2x
f2y
f3x
f3y
f4x
f4y
0.012426 0.047574 1.5248e-13 7928.9 2071.1 2071.1 2071.1 0( )
Los desplazamientos son:
Kt
d1x
d1y
0
0
0
0
0
0
0
10 103
0
7.929 103
2.071 103
2.071 103
2.071 103
0
dd1x
d1y
12.426 10 3
47.574 10 3
57
Introducción al método de los elementos finitos aplicando MATHCAD
Parte del postproceso es graficar la deformada
Figura 2.6. Deformada
0 500 1 103�
0
500
1 103�
DEFORMADA DE LA ESTRUCTURA
Ay
By
Cy
A1y
B1y
C1y
Ax Bx, Cx, A1x, B1x, C1x,
Fuente propia. porograma MathCad
58
Armaduras planas
Cálculo de tensiones Para determinar las tensiones en cada barra utilizamos la expre-
sión (2.24):
59
Introducción al método de los elementos finitos aplicando MATHCAD
3. Para la estructura plana indicada sujeta a las fuerzas indicadas aplicada en el nodo 3, determine la matriz de rigidez, las deforma-ciones, las tensiones
Figura 2.7. Armadura, números en negro nodos, números en rojo elementos
Fuente propia.
Datos de elementos, nodos y conectividad
Grados Longitud
0
45
90
0
180
0
0.785
1.571
0
Lo
1
2
1
1
Area 300 E 200000 Longitud 1500
60
Armaduras planas
Matriz de rigidez
K ( )
cos ( )2
cos ( ) sin ( )
cos ( )2
cos ( ) sin ( )( )
cos ( ) sin ( )
sin ( )2
cos ( ) sin ( )( )
sin ( )2
cos ( )2
cos ( ) sin ( )( )
cos ( )2
cos ( ) sin ( )
cos ( ) sin ( )( )
sin ( )2
cos ( ) sin ( )
sin ( )2
n 1 elementos
K 1
1
0
1
0
0
0
0
0
1
0
1
0
0
0
0
0
K 2
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
K 3
1
0
1
0
0
0
0
0
1
0
1
0
0
0
0
0
Matriz de ceros
Kt1 2 nodos 0 Kt2 nodos 1 0
Kt
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
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
0
0
0
0
fuerza
f1x
f1y
0
0
5000
3000
f4x
f4y
f1x
d
0
0
d2x
d2y
d3x
d3y
0
0
d2x
61
Introducción al método de los elementos finitos aplicando MATHCAD
Kt
p conectividadn 1
q conectividadn 2
r q p
Ktp p Ktp pK n 1 1
Lon
Ktp p 1 Ktp p 1
K n 1 2 Lon
Ktp q r Ktp q r
K n 1 3 Lon
Ktp q r 1 Ktp q r 1
K n 1 4 Lon
Ktp 1 p Ktp 1 pK n 2 1
Lon
Ktp 1 p 1 Ktp 1 p 1
K n 2 2 Lon
Ktp 1 q r Ktp 1 q r
K n 2 3 Lon
Ktp 1 q r 1 Ktp 1 q r 1
K n 2 4 Lon
Ktq r p Ktq r pK n 3 1
Lon
Ktq r p 1 Ktq r p 1
K n 3 2 Lon
Ktq r q r Ktq r q r
K n 3 3 Lon
Ktq r q r 1 Ktq r q r 1
K n 3 4 Lon
Ktq r 1 p Ktq r 1 pK n 4 1
Lon
Ktq r 1 p 1 Ktq r 1 p 1
K n 4 2 Lon
Ktq r 1 q r Ktq r 1 q r
K n 4 3 Lon
Ktq r 1 q r 1 Ktq r 1 q r 1
K n 4 4 Lon
Kt
n 1 elementosfor
Kt Kt EArea
Longitud
62
Armaduras planas
Kt
5.414 104
1.414 104
4 104
0
1.414 104
1.414 104
0
0
1.414 104
5.414 104
2.449 10 12
4 104
1.414 104
1.414 104
0
0
4 104
2.449 10 12
8 104
2.449 10 12
4 104
0
0
0
0
4 104
2.449 10 12
4 104
2.449 10 12
0
0
0
1.414 104
1.414 104
4 104
2.449 10 12
5.414 104
1.414 104
0
0
1.414 104
1.414 104
0
0
1.414 104
1.414 104
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Kt d fuerza solve
d2x
d2y
d3x
d3y
f1x
f1y
f4x
f4y
0.049999999999999995558 9.6156608458240956913e-19 0.099999999999999991117 0.11213203435596427285 5000.0 3000.0 0 0( )
d2x 0.05 d2y 0 d3x 0.0999 d3y 0.11213
Kt
0
0
d2x
d2y
d3x
d3y
0
0
4.999 103
2.999 103
4
1.222 10 13
4.995 103
2.999 103
0
0
d
0
0
d2x
d2y
d3x
d3y
0
0
0
0
0.05
0
0.1
0.112
0
0
Nodeforx
0
Longitud
Longitud
0
0
Longitud
Nodefory
0
0
Longitud
Longitud
0
Longitud
deforx
0
Longitud df3
Longitud df5
0
0
Longitud df5
defory
0
0 df4
Longitud df6
Longitud
0
Longitud df6
63
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 2.7. Deformada de la armadura de la Fig. 2.6
Fuente propia.
Resolución con Matlab según A.J.M. FerreiraDebido a la opción de evaluación simbólica de MathCAD, problemas
con mayor número de barras requieren una programación más elabo-rada por lo que se ha optado por la alternativa presentada por A.J.M. Ferreira, “Matlab Codes for Finite Element Analysis”, los archivos .m que se necesitan para resolver problemas sencillos de armaduras pla-nas se presentan. Los archivos se descargan de la página: http://ex-tras.springer.com/2009/978-1-4020-9199-5/MATLABsoftware
Figura 2.8. Armadura de tres barras
Fuente: A.J.M. Ferreira.
64
Armaduras planas
Las propiedades físicas se dan con los siguientes comandos:
%................................................................
% MATLAB codes for Finite Element Analysis
% problem4.m
% Antonio Ferreira 2008
% clear memory
clear all
clc
% E; modulus of elasticity
% A: area of cross section
% L: length of bar
E=200000; A=1000; EA=E*A;
La generación de coordenadas se obtiene ingresandoNúmero de elementos = 3Número de Nodos = 4Conexiones = [1 2; 1 3; 1 4]Coordenadas de los nodos = [0 0; 0 1200; 1200 1200; 1200 0];
% generation of coordinates and connectivities
numberElements=3;
numberNodes=4;
%conectividad
elementNodes=[1 2;1 3;1 4];
nodeCoordinates=[ 0 0;0 1200;1200 1200;1200 0];
xx=nodeCoordinates(:,1);
yy=nodeCoordinates(:,2);
xx = 0
0
1200
1200
yy = 0
1200
1200
0
65
Introducción al método de los elementos finitos aplicando MATHCAD
Aplicación de la fuerza según el siguiente esquema
Figura 2.9.Ordenamiento de fuerzas
Fuente: A.J.M. Ferreira.
Por lo tanto corresponde una fuerza en (2) de -10000
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
GDof=2*numberNodes; % GDof: total number of degrees of freedom
displacements=zeros(GDof,1);
force=zeros(GDof,1);
% applied load at node 2
force(2)=-10000.0
El vector fuerza queda:force = 0 -10000 0 0 0 0 0 0
66
Armaduras planas
Llamado a subrutina para la evaluación de la matriz de rigidez formStiffness2Dtruss , en este punto de dan los grados de libertad res-tringidos y se direcciona a la subrutina solution
% computation of the system stiffness matrix
[stiffness]=...
formStiffness2Dtruss(GDof,numberElements,elementNodes,numberNodes,nodeCoordinates,xx,yy,EA);
% boundary conditions and solution
prescribedDof= [3 4 5 6 7 8]’;
% solution
displacements=solution(GDof,prescribedDof,stiffness,force);
Postprocesado
% drawing displacements
us=1:2:2*numberNodes-1;
vs=2:2:2*numberNodes;
figure
L=xx(2)-xx(1);
%L=node(2,1)-node(1,1);
XX=displacements(us);YY=displacements(vs);
dispNorm=max(sqrt(XX.^2+YY.^2));
scaleFact=15000*dispNorm;
clf
hold on
drawingMesh(nodeCoordinates+scaleFact*[XX YY],elementNodes,’L2’,’k.-’);
drawingMesh(nodeCoordinates,elementNodes,’L2’,’k.--’);
% stresses at elements
stresses2Dtruss(numberElements,elementNodes,...
xx,yy,displacements,E)
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,GDof,prescribedDof)
67
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 2.9.Postprocesado
Subrutinas formstiffness2dtruss:
function [stiffness]=...
formStiffness2Dtruss(GDof,numberElements,elementNodes,numberNodes,nodeCoordinates,xx,yy,EA);
stiffness=zeros(GDof);
% computation of the system stiffness matrix
for e=1:numberElements;
% elementDof: element degrees of freedom (Dof)
indice=elementNodes(e,:) ;
elementDof=[ indice(1)*2-1 indice(1)*2 indice(2)*2-1 indice(2)*2] ;
xa=xx(indice(2))-xx(indice(1));
ya=yy(indice(2))-yy(indice(1));
length_element=sqrt(xa*xa+ya*ya);
C=xa/length_element;
S=ya/length_element;
k1=EA/length_element*...
[C*C C*S -C*C -C*S; C*S S*S -C*S -S*S;
-C*C -C*S C*C C*S;-C*S -S*S C*S S*S];
stiffness(elementDof,elementDof)=...
stiffness(elementDof,elementDof)+k1;
end
68
Armaduras planas
Subrutina solution:
%................................................................
function displacements=solution(GDof,prescribedDof,stiffness,force)
% function to find solution in terms of global displacements
activeDof=setdiff([1:GDof]’, ...
[prescribedDof]);
U=stiffness(activeDof,activeDof)\force(activeDof);
displacements=zeros(GDof,1);
displacements(activeDof)=U;
Subrutina tensiones
function stresses2Dtruss(numberElements,elementNodes,...
xx,yy,displacements,E)
% stresses at elements
for e=1:numberElements
indice=elementNodes(e,:);
elementDof=[ indice(1)*2-1 indice(1)*2 indice(2)*2-1 indice(2)*2] ;
xa=xx(indice(2))-xx(indice(1));
ya=yy(indice(2))-yy(indice(1));
length_element=sqrt(xa*xa+ya*ya);
C=xa/length_element;
S=ya/length_element;
sigma(e)=E/length_element* ...
[-C -S C S]*displacements(elementDof);
end
disp('stresses')
sigma'
Considere la armadura plana mostrada en la figura dado:Stresses ans = 7.9289 2.9289 -2.0711
69
Introducción al método de los elementos finitos aplicando MATHCAD
Subrutina desplazamientos y reacciones
%................................................................
function outputDisplacementsReactions...
(displacements,stiffness,GDof,prescribedDof)
% output of displacements and reactions in
% tabular form
% GDof: total number of degrees of freedom of
% the problem
% displacements
disp('Displacements')
%displacements=displacements1;
jj=1:GDof; format
[jj' displacements]
% reactions
F=stiffness*displacements;
reactions=F(prescribedDof);
disp('reactions')
[prescribedDof reactions]
70
Armaduras planas
Resolución con ANSYS APDL
Determine las deflexiones nodales, fuerzas de reacción y tensiones para el sistema de armadura indicado: (E = 200GPa, A = 1021mm2).
Figura 2.10.Ejemplo de Armadura
Fuente propia
Las unidades en que se trabajará son mm y N/mm2
1. 1.- Cambie el título del programa a tutorial de armadura para puente:En el Utility menu bar seleccione File < Change Title:
Luego ir a Plot < Replot para visualizar el titulo
71
Introducción al método de los elementos finitos aplicando MATHCAD
2. Crear la geometríaIngresar Keypoints de acuerdo a la tabla siguiente:
Tabla 2.1. Coordenadas en los nodos.
Keypoint x y
1 0 0
2 3650 0
3 7300 0
4 10950 0
5 14600 0
6 0 6000
7 3650 5250
8 7300 4500
9 10950 3750
10 14600 3000
Usar Preprocessor < Modeling < Create < Keypoints < In Active CS
Clic en Apply para aceptar el resultado, El resultado es el siguiente:
Figura 2.11.Visualización de keypoints
Fuente. Software Ansys APDL
72
Armaduras planas
Si se desea cambiar el fondo de pantalla mediante:Plot Ctrls < Style < Colors < Reverse VideoPara dibujar las líneas se selecciona:Preprocessor < Modeling < Create < Lines < Lines < In Active Coord
Con el mouse se selecciona señalando los Keypoints hasta obtener la armadura, no se debe saltar ningún Keypoint:
Figura 2.12.Visualización de la armadura
Fuente. Software Ansys APDL
73
Introducción al método de los elementos finitos aplicando MATHCAD
En caso de que desaparezcan las líneas las volver a activar con plot < lines
3. Defina el tipo de elementoSeleccione: Element Type < Add/Edit/Delete y 3D finit stn 180
Este elemento finito LINK180 puede modelar armaduras bidimen-sionales o tridimensionales, cables y resortes.
4. Definir propiedades geométricasCon: Real Constants < Add/Edit/Delete, En área colocar 1021 mm2 correspondiente a un miembro estructural rectángulo de 120 x 60 x 3
74
Armaduras planas
5. Definir propiedades del materialEn Preprocessor seleccionar Material Props < Material ModelsSeleccionar material elástico, lineal e isotrópico, llenar en el campo de Módulo de elasticidad el valor de E = 200000 MPa, y coeficiente de Poisson, 0.3.
Ok y Material < Exit
75
Introducción al método de los elementos finitos aplicando MATHCAD
6. MalladoEn este caso los elementos finitos son las barras por lo tanto el nú-mero de divisiones a utilizarse es 1Seleccione Meshing < MeshTool
Lines Set < Pick All
76
Armaduras planas
Luego Meshing < MeshTool < MeshPara numerar los elementos ir a Utility Menu seleccionar PlotCtrls < Numbering...Y numerar según nodos o elementos
Figura 2.13.Mallado de la armadura
Fuente. Software Ansys APDL
77
Introducción al método de los elementos finitos aplicando MATHCAD
7. Fase de soluciónDefinimos el tipo de análisisSe selecciona: Solution Menu < Analysis Type < New Analysis <Static.
Y se verifica que este selecciuonada la opción StaticEl nodo 2 es totalmente restringido, por lo tanto se selecciona a:Define Loads < Apply < Structural < Displacement < On Keypoints
78
Armaduras planas
El nodo 1 está sujeto a un rodillo por lo que no se permite el des-plazamiento en x
Las restricciones se observan en la Fig. 2.14
Figura 2.14.Restricciones de la armadura
Fuente. Software Ansys APDL
79
Introducción al método de los elementos finitos aplicando MATHCAD
8. Aplicar cargasAplicar las cargas especificadas en la Fig.2.10 respectivamente.
Figura 2.15.Sistema de cargas
Fuente. Software Ansys APDL
Resolver el sistema de ecuaciones mediante: Solution < Solve < Cu-rrent LS. Aparecerá el mensaje: Solution is done
80
Armaduras planas
9. Postprocesado revisión de resultadosGeneral Postproc < List Results < Reaction Solution, escoger All ItemsPlot Results < Deformed Shape
Figura 2.16. Deformada del sistema
Fuente. Software Ansys APDL
81
Introducción al método de los elementos finitos aplicando MATHCAD
La deflexión puede ser obtenida también como una lista con: Gene-ral Postproc < List Results < Nodal Solution
Figura 2.17.Esfuerzos de Von Mises
Fuente. Software Ansys APDL
82
Armaduras planas
10. Determinación de la tensión axialPara definir el comando correspondiente a la tensión axial efectuar el siguiente procedimiento: Escribir en la línea de comandos HELP LINK180 para obtener la información del elemento en particular.
Figura 2.17.Help Ansys APDL
Fuente. Software Ansys APDL
83
Introducción al método de los elementos finitos aplicando MATHCAD
Se despliega la información del elemento y debemos buscar las va-riables de salida
La variable Sxx es la que se quiere graficar. Este parámetro está asociado con el Item LS como se observa en la tabla 180.2
84
Armaduras planas
Se selecciona entonces a:General Postprocessor < Element Table < Define TableSeleccionar Add y escribir en los cuadros de dialogo Tensión axial, seleccionarBy sequence num y ponga 1 después de LS
Grafique las tensiones seleccionando General Postproc < Plot Re-sults < Contour Plot Element Table < Plot Elem Table
Figura 2.18.Tensiones en las barras
Fuente. Software Ansys APDL
85
Introducción al método de los elementos finitos aplicando MATHCAD
Donde se especifican las tensiones axiales en MPa y los elementos en rojo están en tensión y los azules en compresión. Con estos datos y el perfil estructural se podría determinar el diseño de la estructura. Igual se pueden listar las tensionesGeneral Postproc < List Results < Element Table data < TENSION_AXIAL
87
Introducción al método de los elementos finitos aplicando MATHCAD
Capítulo III
Análisis por el método de elementos finitos
de elementos a flexión:Vigas
89
Introducción al método de los elementos finitos aplicando MATHCAD
Determinar la matriz de rigidez del elemento viga
Una viga es un elemento largo y delgado, sujeto a carga transver-sal que produce significativos efectos de flexión, esta flexión es medida como un desplazamiento transversal y una rotación, por lo tanto los grados de libertad por nodo son dos. La viga es un elemento fundamen-tal en la construcción de edificaciones, máquinas, etc.
En un elemento finito para una viga se presentan los siguientes parámetros• Fuerzas locales: f i• Momentos flexores locales: m i, positivo en la dirección anti horaria• Desplazamientos locales: d i • Rotaciones: Фi
Figura 3.1.Viga IPN
Figura 3.2.Elemento viga
Fuente propia
90
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Ecuaciones constitutivasLa ecuación diferencial que gobierna una viga se la obtiene del
desarrollo de la viga de Euler - Bernoulli. Las secciones planas inicial-mente perpendiculares al eje de la viga, siguen siendo perpendiculares al eje de la viga una vez que se ha curvado.
Figura 3.3.Elementodiferencial viga
Efectuando el análisis en un elemento diferencial de viga usando las ecuaciones de equilibrio se obtiene:
(3.1)
(3.2)
Figura 3.4. Desplazamiento en una viga
De la figura 3.4, ν(x) es la función que representa el desplaza-miento transversal de la viga, por otro lado el ángulo φ o pendiente está dado por:
(3.3)
91
Introducción al método de los elementos finitos aplicando MATHCAD
Del análisis de la teoría de vigas se tiene que la curvatura de una viga es:
Figura 3.5. Ecuaciones constitutivas
De la ley de Hooke, la deformación unitaria en un arco es directa-mente proporcional a la variación radial
(3.4)
La curvatura se calcula mediante la expresión conocida:
(3.5)
Y para pequeña curvatura se tiene simplemente:
(3.6)
Por lo tanto igualando (3.4) con (3.5) el Momento es igual a:
(3.7)
92
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Puesto que la fuerza cortante es la derivada del momento, según (3.2) se obtiene:
(3.8)
Y la carga distribuida w(x) la derivada del cortante, según (3.1):
(3.9)
Si se elige un polinomio arbitrario para v(x) y se lo reemplaza en las ecuaciones constitutivas se hallaría la matriz del elemento finito:
A continuación se establecen los pasos para hallar la matriz que relaciona fuerzas con desplazamientos
1. Selección del tipo de elementoElemento viga con 4 grados de libertad, 2 por nodo, sujeto a fuer-zas transversales F y momentos M
2. Seleccionar una función de desplazamiento Se debe escoger por lo tanto un polinomio con cuatro constantes y de grado 3 que satisfaga la ecuación diferencial.
ν(x)= a1 x3+ a2 x2+ a3 x+a4 (3.10)
Figura 3.6.Elemento viga
Fuente propia
93
Introducción al método de los elementos finitos aplicando MATHCAD
Se precisa determinar a1, a2, a3 y a4 en función de las condiciones de frontera que son:
Resolviendo se obtiene la función de desplazamiento:
a1 y a2, se obtienen al resolver el siguiente sistema de ecuaciones lineales:
(3.11)
Para resolverlas se utilizara la resolución simbólica de MathCAD.
Given
a1 L3 a2 L2
1 L d1y d2y 0
3 a1 L2 2 a2 L 1 2 0
94
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Find a1 a2( )
2 d1y 2 d2y L 1 L 2
L3
3 d1y 3 d2y 2 L 1 L 2
L2
Por lo tanto ν(x) queda:
Reagrupando y factorando en función de los desplazamientos o variables de campo
(3.12)
En forma matricial se obtiene:
(3.13)
Donde los términos de la matriz de funciones de forma N están dados por
(3.14)
95
Introducción al método de los elementos finitos aplicando MATHCAD
En forma matricial queda:
� x( ) N1 x( ) N2 x( ) N3 x( ) N4 x( )( )
d1y
�1
d2y
�2
������
�÷÷÷÷�
×:=
(3.15)
N1, N2, N3, N4 son las funciones de forma para elementos viga y se denominan polinomios de interpolación cúbica de Hermite
3. Definir las relaciones Deformación/Desplazamiento y Esfuer-zo/Deformación Se utilizan las relaciones anteriores (3.1) y (3.2):
4. Definir la matriz de rigidez Se deriva la función de forma y se la reemplaza en (3.1) y (3.2) eva-luándose en los extremos del elemento finito viga, a continuación se aprecia las sucesivas derivadas del vector fila, obtenidas con la herramienta de derivada simbólica de MathCAD (3.15):
1
L32 x3 3 L x2
L3
1
L3x3 L 3 L2
x2 x L3
1
L32 x3 3 L x2
1
L3x3 L L2 x2
6 x2 6 L x
L3L3 6 L2
x 3 L x2
L36 x2 6 L x
L3
2 L2 x 3 L x2
L3
6 L 12 x
L3
6 L2 6 L x
L3
6 L 12 x
L32 L2 6 L x
L3
(3.16)
96
Análisis por el método de elementos finitos de elementos a flexión. Vigas
(3.17)
Reordenando en forma matricial se obtiene finalmente:
(3.18)
Ejemplos cargas puntuales1. Resuelva la deflexión, rotación y momentos de una viga en can-tiléver, compare con la ecuación clásica
Figura 3.7.Elemento viga empotrado en un extremo y cargado en el otro
97
Introducción al método de los elementos finitos aplicando MATHCAD
Como se puede apreciar los resultados generados por el método de elementos finitos coinciden con los resultados obtenidos por los méto-dos de mecánica de materiales, esto es integrando dos veces la ecua-ción diferencial de la viga y utilizando las condiciones de frontera.
(3.19)
A continuación se grafica la deflexión de los puntos intermedios de la viga obtenida tanto por la función de interpolación como por la solu-ción de la ecuación diferencial.
x 0 0.1 L L 1000 E 200000 P 1000
d2yL3 P3 E I
2L2 P2 E I
b 10 h 100 I
112
b h3
FEA x( )1
L32 x3 3 x2
L L3 d1y
1
L3x3 L 2 L2
x2 x L3
11
L32 x3 3 x2
L d2y1
L3x3 L x2 L2
2
Se obtiene la deflexión por el método clásico integrando dos veces la ecuación diferencial:
Figura 3.8.Deflexión viga
0 100 200 300 400 500 600 700 800 900 1 103�
2.5�2.25�
2�1.75�1.5�
1.25�1�
0.75�0.5�
0.25�0
�FEA x( )�
�EDO x( )�
500
x
Fuente propia
98
Análisis por el método de elementos finitos de elementos a flexión. Vigas
El método de los elementos finitos coincide perfectamente para los puntos intermedios. Para hallar el momento aplicamos la ecuación (3.1),
Momento x( ) E I 2xFEA x( )d
d
2
Figura 3.9.Momento viga
0 100 200 300 400 500 600 700 800 900 1 103�
01 105�
2 105�
3 105�
4 105�
5 105�
6 105�
7 105�
8 105�
9 105�
1 106�
Momento x( )
x
Fuente propia
El momento teórico se lo obtiene multiplicando la fuerza en el ex-tremo por la longitud de la viga: M teórico P (L) = 1000 (1000) = 1 e6
2. Considere la viga de la figura 3.10 para ilustrar el procedimiento de ensamblaje, E = 200000 MPa, Sección b = 10mm, L = 1000mm, h = 100mm y un límite de fluencia de 250 MPa
Figura 310.Viga con dos elementos
Fuente propia
99
Introducción al método de los elementos finitos aplicando MATHCAD
Debido a la carga intermedia se precisa discretizar la viga con dos elementos
Nodo 1 a 2 Nodo 2 a 3
Expandiendo las matrices se obtiene:
d1 1 d2 2 d3 3( ) d1 1 d2 2 d3 3( )
K1E I
L3
12
6 L
12
6 L
0
0
6 L
4 L2
6 L
2 L2
0
0
12
6 L
12
6 L
0
0
6 L
2 L2
6 L
4 L2
0
0
0
0
0
0
0
0
0
0
0
0
0
0
E
K2E I
L3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
12
6 L
12
6 L
0
0
6 L
4 L2
6 L
2 L2
0
0
12
6 L
12
6 L
0
0
6 L
2 L2
6 L
4 L2
E
Sumando
K1 K2
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
(3.20)
100
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Una manera de automatizar el ensamblaje es mediante el siguiente algoritmo de MathCAD
k1E I
L3
12
6 L
12
6 L
6 L
4 L2
6 L
2 L2
12
6 L
12
6 L
6 L
2 L2
6 L
4 L2
E
n 2
K2 n 1 2 n 1 0
K
K2 i 2 m 2 i 2 p K2 i 2 m 2 i 2 p k1m p
p 0 3for
K
m 0 3for
i 1 nfor
K
k1
K
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
Las ecuaciones de gobierno por lo tanto están dadas por:
(3.21)
Las condiciones de frontera sond1y = 0, Ф1 =0, d3y = 0Adicionalmente se conoceF2Y = -1000 N, M2 = + 1000 N mm, M3= 0
101
Introducción al método de los elementos finitos aplicando MATHCAD
Se debe determinard2y, F1y, M1 y F3y
En MathCAD esto puede resolverse automáticamente mediante las herramientas de algebra simbólica. Donde en primer lugar se plantean los parámetros conocidos:
d1y 0 1 0 d3y 0 F2y 1000 M2 1000 M3 0
La ecuación que relaciona fuerzas con desplazamientos es:
Se activan las herramientas de evaluación simbólica
Figura 3.11.Symbolic Keyword toolbar
Programa MathCAD
102
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Se escoge solve y separadas por coma “,” las incógnitas a encon-trarse en formato matricial
F1y
M1
F2y
M2
F3y
M3
K
d1y
1
d2y
2
d3y
3
solve
d2y
2
3
F1y
M1
F3y
7 L3 P 3 L2
M
96 E I
L2 P 5 L M
32 E I
L2 P L M
8 E I9 M 11 L P
16 LM8
3 L P8
9 M 5 L P
16 L
Los resultados son dados automáticamente:
Para evaluar las funciones se utilizan los datos numéricos siguientes:Módulo de elasticidad E 200000:= Altura de la viga h 100:= Ancho de la viga b 10:=
Inercia I112
b× h3×:=
Longitud del elemento L 1000:=
Carga P 1000:=
Momento M 1000:=
Para graficar la deformada reemplazando los desplazamientos y rotaciones en las funciones de interpolación:
d2y7 L3 P 3 L2
M
96 E I
103
Introducción al método de los elementos finitos aplicando MATHCAD
�FEA1 x( )1
L32 x3
× 3 x2× L×� L3
+( )× d1y×1
L3x3 L× 2 L2
× x2×� x L3
×+( )× �1×+1
L32� x3
× 3 x2× L×+( )× d2y×+
1
L3x3 L× x2 L2
×�( )× �2×+:=
�FEA2 x( )1
L32 x L�( )3
× 3 x L�( )2× L×� L3
+�� ��× d2y×1
L3x L�( )3 L× 2 L2
× x L�( )2×� x L�( ) L3
×+�� ��× �2×+1
L32� x L�( )3
× 3 x L�( )2× L×+�� ��× d3y×+
1
L3x L�( )3 L× x L�( )2 L2
×��� ��× �3×+:=
FEA x( ) FEA1 x( ) 0 x Lif
FEA2 x( ) L x 2 Lif
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�0.45�
0.405�0.36�
0.315�0.27�
0.225�0.18�
0.135�0.09�
0.045�0
Deformada
�FEA x( )
x
La deflexión máxima es: -0.45 mmLa rotación es:
x( )xFEA x( )d
d
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�8� 10 4��
6.4� 10 4��4.8� 10 4��3.2� 10 4��1.6� 10 4��
01.6 10 4��3.2 10 4��4.8 10 4��6.4 10 4��
8 10 4��
Rotación
� x( )
x
El momento se lo calcula mediante: M x( ) E I 2xFEA x( )d
d
2
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�4� 105�
3.2� 105�2.4� 105�1.6� 105�
8� 104�0
8 104�
1.6 105�2.4 105�3.2 105�
4 105�
Momento
0M x( )
x
104
Análisis por el método de elementos finitos de elementos a flexión. Vigas
El cortante se lo calcula con: Vcortante x( ) E I 3xFEA x( )d
d
3
0 200 400 600 800 1 103 1.2 103 1.4 103 1.6 103 1.8 103 2 1034002801604080
200320440560680800
Cortante
Vcortante x( )
x
La tensión se calcula mediante:
� x( )h�
2E× 2x
�FEA x( )d
d
2×:=
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�20�
15.5�11�
6.5�2�
2.57
11.516
20.525
Tensión
0� x( )
x
El factor de seguridad mínimo es:
N x( )Sy x( )
25022.5
11.111=
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�10
14
18
22
26
30
34
38
42
46
50Factor de seguridad
N x( )
x
105
Introducción al método de los elementos finitos aplicando MATHCAD
3. Considere la viga de la figura, E = 200000 MPa, Sección b = 10mm, L = 1000mm, h = 100mm y un límite de fluencia de 250 MPa
Figura 3.12.Viga simplemente apoyada
K2 n 1 2 n 1 0 n 2
K
K2 i 2 m 2 i 2 p K2 i 2 m 2 i 2 p km p
p 0 3for
K
m 0 3for
i 1 nfor
K
k
K
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
Condiciones de frontera:
d1y 0 d3y 0 M1 0 F2y P P M2 0 M3 0
106
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Solución
F1y
M1
F2y
M2
F3y
M3
K
d1y
1
d2y
2
d3y
3
solve
1
2
3
d2y
F1y
F3y
L2 P4 E I
0L2 P4 E I
L3 P6 E I
P2
P2
Parámetros
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�1.5� 10 3��1.2� 10 3��
9� 10 4��6� 10 4��3� 10 4��
03 10 4��6 10 4��9 10 4��
1.2 10 3��1.5 10 3��
Rotación
� x( )
x
107
Introducción al método de los elementos finitos aplicando MATHCAD
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�0
5 104�1 105�
1.5 105�2 105�
2.5 105�3 105�
3.5 105�4 105�
4.5 105�5 105�
Momento
0
M x( )
x
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�600�480�360�240�120�
0120240360480600
Cortante
Vcortante x( )
x
x( )h2
E 2xFEA x( )d
d
2
108
Análisis por el método de elementos finitos de elementos a flexión. Vigas
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�30�27�24�21�18�15�12�9�6�3�0
Tensión0
� x( )
x
Sy 250
N x( )Sy x( )
0 200 400 600 800 1 103� 1.2 103� 1.4 103� 1.6 103� 1.8 103� 2 103�5
9.5
14
18.5
23
27.5
32
36.5
41
45.5
50Factor de seguridad
N x( )
x
4. Determinar los desplazamientos nodales y las rotaciones, las fuerzas globales nodales, de la siguiente viga:
Figura 3.13.Viga empotrada tres elementos
109
Introducción al método de los elementos finitos aplicando MATHCAD
En este caso se utilizaran tres elementos y la matriz de rigidez queda:
K
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
Condiciones de frontera:
Solución:
F1y
M1
F2y
M2
F3y
M3
F4y
M4
K
d1y
1
d2y
2
d3y
3
d4y
4
solve
F1y
M1
d2y
2
d3y
3
F4y
M4
20 P27
4 L P
9
8 L3 P
81 E I2 L2 P
27 E I11 L3
P162 E I
5 L2 P
54 E I
7 P27
2 L P
9
110
Análisis por el método de elementos finitos de elementos a flexión. Vigas
0 300 600 900 1.2 103� 1.5 103� 1.8 103� 2.1 103� 2.4 103� 2.7 103� 3 103�0.7�
0.63�0.56�0.49�0.42�0.35�0.28�0.21�0.14�0.07�
0Deformada
�FEA x( )
x
Ejemplos cargas puntuales longitud y areas de inercia variables
Determinar los desplazamientos nodales y las rotaciones, del si-guiente eje:
Figura 3.13.Eje de máquina
Fuente propia
111
Introducción al método de los elementos finitos aplicando MATHCAD
II1
I2
I2
LL1
L2
L1
k I L( )E I
L3
12
6 L
12
6 L
6 L
4 L2
6 L
2 L2
12
6 L
12
6 L
6 L
2 L2
6 L
4 L2
E
Número de elementos finitos
n 2K2 n 1 2 n 1 0
K
K2 i 2 m 2 i 2 p K2 i 2 m 2 i 2 p k Ii 1 Li 1 m p
p 0 3for
K
m 0 3for
i 1 nfor
K
k
La matriz de riguidez toma la siguiente forma:
K
12 E I1
L13
6 E I1
L12
12 E I1
L13
6 E I1
L12
0
0
6 E I1
L12
4 E I1L1
6 E I1
L12
2 E I1L1
0
0
12 E I1
L13
6 E I1
L12
12 E I1
L1312 E I2
L23
6 E I2
L226 E I1
L12
12 E I2
L23
6 E I2
L22
6 E I1
L12
2 E I1L1
6 E I2
L226 E I1
L12
4 E I1L1
4 E I2L2
6 E I2
L22
2 E I2L2
0
0
12 E I2
L23
6 E I2
L22
12 E I2
L23
6 E I2
L22
0
0
6 E I2
L22
2 E I2L2
6 E I2
L22
4 E I2L2
112
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Condiciones de frontera:
Solución
Datos:
113
Introducción al método de los elementos finitos aplicando MATHCAD
0 27 54 81 108 135 162 189 216 243 2701.8�
1.62�1.44�1.26�1.08�0.9�
0.72�0.54�0.36�0.18�
0Deformada
�FEA x( )
x
0 27 54 81 108 135 162 189 216 243 2700.015�0.011�
7� 10 3��3� 10 3��1 10 3��5 10 3��9 10 3��
0.0130.0170.0210.025
Rotación
� x( )
x
0 27 54 81 108 135 162 189 216 243 2700
7 103�1.4 104�2.1 104�2.8 104�3.5 104�4.2 104�4.9 104�5.6 104�6.3 104�
7 104�
Momento
0
M x( )
x
Problemas propuestos
Efectuar el diseño de un eje que debe transmitir 2.5 hp a 1725 rpm, el factor de seguridad mínimo es 2.5, considere mínimo 4 elemen-tos finitos, se debe trabajar en 2 planos.
114
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Las medidas tentativas del eje son:
Figura 3.14.Eje de máquina
Carga distribuida
Si se observa en las edificaciones constataremos que las vigas a diferencia de los ejes deben soportar carga distribuida así como cargas concentradas. La estrategia que se emplea es discretizar la carga distribuida reemplazándola por cargas y momentos en los ex-tremos que generen el mismo desplazamiento que el producido por la carga distribuida.
115
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 3.15.Discretización de cargas
Se puede utilizar el método del trabajo equivalente para reempla-zar una carga distribuida por un grupo de cargas concentradas
(3.22)
Donde w(x) es la función de carga distribuida y ν(x) es el despla-zamiento transversal. El trabajo debido a las fuerzas discretas nodales está dado por:
(3.23)
Se puede determinar los momentos y fuerzas nodales igualando
W distribuida=Wdiscreto (3.24)
En el caso de la viga con carga uniformemente distribuida se eva-lúa la integral que equipara los trabajos según (3.24)
(3.25)
Donde ν(x) ya fue desarrollado anteriormente según (3.12)
116
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Para hallar m1 se reemplaza en la integral los siguientes valores, note que F1 es uno
1 1 2 0 d1y 0 d2y 0
0
L
xw2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
dL2 w
12
(3.26)
Para hallar m2 se reemplaza en la integral los siguientes valores, note ahora que F2 es uno:
1 0 2 1 d1y 0 d2y 0
0
L
xw2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
dL2 w
12
(3.27)
1 0 2 0 d1y 1 d2y 0
0
L
xw2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
dL w
2
(3.28)
1 0 2 0 d1y 0 d2y 1
0
L
xw2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
dL w
2
(3.29)
117
Introducción al método de los elementos finitos aplicando MATHCAD
Por la tanto de esta manera se hallaron las fuerzas nodales equi-valentes:
Figura 3.15.Fuerzas equivalentes
En formato matricial las fuerzas quedan de esta manera
( 3.30)
A continuación se presentan una serie de problemas resueltos que involucran carga distribuida
Para la viga en cantiléver sujeta a la carga distribuida determinar los desplazamientos y rotaciones involucrados
Figura 3.15.Viga carga distribuida
La ecuación general que se utiliza para calcular vigas concentradas o distribuida es la siguiente
(3.31)
118
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Donde F0 son las fuerzas nodales equivalentes:Para determinar los desplazamientos nodales, se parte de la ecua-
ción general:
Donde en la matriz de fuerzas se colocan las fuerzas equivalentes:
La solución se la hace por el método conocido
d1y 0 1 0 f2yw L2
w m2
w L2
12
w
La solución coincide con la solución analítica que se tiene para la viga empotrada
(3.32)
A continuación se halla las fuerzas globales nodales efectivas que es el producto de la matriz de rigidez por el vector desplazamiento, pro-ceso conocido como substitución inversa.
E I
L3
12
6 L
12
6 L
6 L
4 L2
6 L
2 L2
12
6 L
12
6 L
6 L
2 L2
6 L
4 L2
0
0
L4 w8 E I
L3 w6 E I
L w2
5 L2 w12
L w2
L2 w12
119
Introducción al método de los elementos finitos aplicando MATHCAD
Las fuerzas reales se las encuentra restando las obtenidas menos las equivalentes:
(3.33)L w
2
5 L2 w12
L w2
L2 w12
L w2
w L2
12
L w2
L2 w12
L w
L2 w2
0
0
La deformada se la determina utilizando las funciones de interpo-lación y los datos siguientes:
L 1000 w 100 E 200000 x 0 0.1 L
b 10 h 100 I112
b h3 I 8.333 105
d2yL4 w8 E I
2L3 w6 E I
FEA x( )1
L32 x3 3 x2
L L3 d1y
1
L3x3 L 2 L2
x2 x L3
11
L32 x3 3 x2
L d2y1
L3x3 L x2 L2
2
El desplazamiento obtenido resolviendo la ecuación diferencial corresponde a la fórmula:
ymax x( )1
E Iw x4
24w L x3
6
w L2 x2
4
120
Análisis por el método de elementos finitos de elementos a flexión. Vigas
0 100 200 300 400 500 600 700 800 900 1 103�80�72�64�56�48�40�32�24�16�8�0
Deformada ambos métodos
�FEA x( )
ymax x( )
500
x
Figura 3.16.Deformada
Como se puede notar el desplazamiento es exacto en los nodos, pero existe un error en la parte intermedia. Ahora graficamos el momento.
Momento x( ) E I 2xFEA x( )d
d
2
0 100 200 300 400 500 600 700 800 900 1 103�4.5� 107�
3.95� 107�3.4� 107�
2.85� 107�2.3� 107�
1.75� 107�1.2� 107�6.5� 106�
1� 106�4.5 106�
1 107�
Momento
0
Momento x( )
x
Figura 3.16.Deformada
Con un elemento finito la solución claramente no converge para el momento, y adicionalmente es una función parabólica. En un extremo debe ser cero y en el otro debe dar – 5 e7. Porque sucede esto?, Este es un aspecto intrínseco del método de los elementos finitos, puesto que si la siguiente ecuación de interpolación seleccionada es de grado 3.
121
Introducción al método de los elementos finitos aplicando MATHCAD
No sería compatible con la siguiente ecuación diferencial
Si queremos por tanto que coincida el momento debemos incre-mentar el grado del polinomio o aumentar el número de elementos.
Efectuar el mismo análisis para 2 elementos.Las fuerzas equivalentes se las construyen bajo el siguiente esquema:
Figura 3.17.Fuerzas equivalentes para 2 elementos
La matriz de rigidez es por tanto:
K1 K2
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
Las condiciones de frontera son:
d1y 0 1 0
122
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Las fuerzas equivalentes de los nodos 2 y 3 son:
F2yw L2
w L2
www
M2 0
F3yw L2
w
M3
w L2
12
w
Resolviendo simbólicamente se obtiene:
d 0 0 d2y �2 d3y �3( ):= d2y F F1y M1 w L 0w L2
w L2
12
F1y
FT K dT solve
F1y
M1
d2y
2
d3y
3
3 L w
223 L2
w12
17 L4 w
24 E I
7 L3 w6 E I
2 L4 wE I
4 L3 w3 E I
La substitución inversa y las fuerzas reales se obtienen restando lo siguiente:
K
0
0
17 L4 w
24 E I
7 L3 w6 E I
2 L4 wE I
4 L3 w3 E I
w L2
w L2
12
w L
0
w L2
w L2
12
simplify
2 L w
2 L2 w
0
0
0
0
Los datos para graficar son los siguientes:
E 200000 b 10 h 100 I
112
b h3
L 500 w 100
123
Introducción al método de los elementos finitos aplicando MATHCAD
Y las variables de campo se las toma de la solución
d2y17 L4 w
24 E I
2
7 L3 w6 E I
d3y2 L4 wE I
34 L3 w3 E I
x 0 0.01 2 L
FEA1 x( )1
L32 x3 3 x2
L L3 d1y
1
L3x3 L 2 L2
x2 x L3
11
L32 x3 3 x2
L d2y1
L3x3 L x2 L2
2
FEA2 x( )1
L32 x L( )3 3 x L( )2
L L3 d2y
1
L3x L( )3 L 2 L2
x L( )2 x L( ) L3
21
L32 x L( )3 3 x L( )2
L d3y1
L3x L( )3 L x L( )2 L2
3
FEA x( ) FEA1 x( ) 0 x Lif
FEA2 x( ) L x 2 Lif
ymax x( )
1E I
w x4
24w 2 L x3
6
w 2 L( )2 x2
4
Figura 3.18. Deformada de la viga
0 100 200 300 400 500 600 700 800 900 1 10380726456484032241680
Deformada ambos métodos0
FEA x( )
ymax x( )
x
Como se constata en la Figura la diferencia entre las curvas es im-perceptible.
El momento se lo construye mediante:
M x( ) E I× 2x�FEA x( )d
d
2×:= Mreal x( ) E I× 2x
ymax x( )d
d
2×:=
124
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Como se ve existe una discrepancia entre las funciones. El esfuerzo es:
x( )h2
E 2xFEA x( )d
d
2
Sy 250
El factor de seguridad en este caso usando la relación N x( )Sy
� x( ):= es
0.08 por lo tanto el área de inercia es insuficiente, se propone a los alumnos hallar la Inercia adecuada.
125
Introducción al método de los elementos finitos aplicando MATHCAD
Para la viga en cantiléver sujeta a la carga distribuida y carga concentrada en el extremo resolver los desplazamientos y rotaciones involucrados
Figura 3.15.Viga carga distribuida
Nuevamente se utiliza
E I
L3
12
6 L
6 L
4 L2
d2y
2
wL2 P
wL2
12
solve d2y 23 w L4
8 P L3
24 E I
w L3 3 P L2
6 E I
126
Análisis por el método de elementos finitos de elementos a flexión. Vigas
A continuación se hallan las fuerzas globales nodales efectivas que es el producto
E I×
L3
12
6 L×
12�
6 L×
6 L×
4 L2×
6� L×
2 L2×
12�
6� L×
12
6� L×
6 L×
2 L2×
6� L×
4 L2×
�������
�÷÷÷÷÷�
×
0
0
3 w× L4× 8 P× L3
×+
24 E× I×�
w L3× 3 P× L2
×+
6 E× I×�
���������
�÷÷÷÷÷÷÷�
×
3 w× L4× 8 P× L3
×+
2 L3×
w L3× 3 P× L2
×+
L2�
3 w× L4× 8 P× L3
×+
4 L2×
w L3× 3 P× L2
×+
3 L�
w L3× 3 P× L2
×+
L2
3 w× L4× 8 P× L3
×+
2 L3×
�
3 w× L4× 8 P× L3
×+
4 L2×
2 w L3× 3 P× L2
×+( )×
3 L�
���������������
���������������
® simplify
PL w×
2+
5 w× L2×
12P L×+
P�L w×
2�
L2 w×12
������������
�÷÷÷÷÷÷÷÷÷÷�
®
Finalmente se restan la fuerza nodal equivalente
PL w
2
5 w L2
12P L
PL w
2
L2 w12
L w2
w L2
12
L w2
L2 w12
P L w
w L2
2P L
P
0
La deformada se obtiene:
0 100 200 300 400 500 600 700 800 900 1 103�
6�5.4�4.8�4.2�3.6�
3�2.4�1.8�1.2�0.6�
0
�FEA x( )
500
x
127
Introducción al método de los elementos finitos aplicando MATHCAD
Para la viga en cantiléver sujeta a la carga distribuida y empotrada en los dos extremos
Figura 3.16.Viga carga distribuida
En este caso se utilizarán dos elementos finitos, sería imposible resolver con uno solo. La matriz de rigidez con las condiciones son las siguientes.
128
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Se eliminan las filas y columnas correspondientes y se calcula la deflexión
w L
0
24 E I
L3
0
0
8 E IL
d2y
2
solve d2y 2L4 w24 E I
0
A continuación se halla las fuerzas globales nodales efectivas que es el producto
k d
f1y
m1
f2y
m2
f3y
m3
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
0
0
L4 w24 E I
0
0
0
f1y
m1
f2y
m2
f3y
m3
L w2
L2 w4
L w
0
L w2
L2 w4
Las fuerzas reales se obtiene restando K d – F0:
f1y
m1
f2y
0m3
f3y
m3
L w2
L2 w4
L w
0
L w2
L2 w4
w L2
L2 w12
L w
0
w L2
L2 w12
f1y
m1
f2y
0
f3y
m3
L w
L2 w3
0
0
L w
L2 w3
129
Introducción al método de los elementos finitos aplicando MATHCAD
Con las condiciones se valora:
d2yL4 w24 E I
2 0 d3y 0 3 0
Obteniéndose los siguientes gráficos de la deformada:
0 100 200 300 400 500 600 700 800 900 1 103�
1.6�1.44�1.28�1.12�0.96�0.8�
0.64�0.48�0.32�0.16�
0 0
�FEA x( )
x
Momento:
0 100 200 300 400 500 600 700 800 900 1 103�
8� 106�
6.4� 106�
4.8� 106�
3.2� 106�
1.6� 106�
01.6 106
�
3.2 106�
4.8 106�
6.4 106�
8 106�
0M x( )
x
Cortante:
0 100 200 300 400 500 600 700 800 900 1 103�
3� 104�
2.4� 104�
1.8� 104�
1.2� 104�
6� 103�
06 103
�1.2 104
�1.8 104
�2.4 104
�3 104
�
Vcortante x( )
x
130
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Para la misma viga en cantiléver sujeta a la carga distribuida y empotrada en los dos extremos, hacer el cálculo para tres elementos
Figura 3.17.Viga carga distribuida
K1 K2 K3
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
d1y 0 1 0 F2yw L2
w L2
w m2 0
131
Introducción al método de los elementos finitos aplicando MATHCAD
F1y
M1
F2y
M2
F3y
M3
F4y
M4
K
d1y
1
d2y
2
d3y
3
d4y
4
solve
F1y
M1
d2y
2
d3y
3
F4y
M4
L w2 L2 w
3L4 w6 E I
L3 w6 E I
L4 w6 E I
L3 w6 E I
L w2 L2 w
3
La substitución inversa se obtiene de:
K
d1y
1
d2y
2
d3y
3
d4y
4
L w2
L2 w
12
L w
0
L w
0
L w2
L2 w12
simplify
3 L w2
3 L2 w
4
0
0
0
0
3 L w2
3 L2 w
4
La deflexión es por tanto
FEA x( ) FEA1 x( ) 0 x Lif
FEA2 x( ) L x 2 Lif
FEA3 x( ) 2 L x 3 Lif
0 300 600 900 1.2 103� 1.5 103� 1.8 103� 2.1 103� 2.4 103� 2.7 103� 3 103�140�126�112�98�84�70�56�42�28�14�0
Deformada
�FEA x( )
x
132
Análisis por el método de elementos finitos de elementos a flexión. Vigas
0 300 600 900 1.2 103� 1.5 103� 1.8 103� 2.1 103� 2.4 103� 2.7 103� 3 103�8� 107�
6.8� 107�5.6� 107�4.4� 107�3.2� 107�
2� 107�8� 106�4 106�
1.6 107�2.8 107�
4 107�
Momento
0M x( )
MAB x( )
x
Vcortante x( ) E I× 3x�FEA x( )d
d
3×:=De igual manera se puede proceder con más elementos
0 100 200 300 400 500 600 700 800 900 1 1038 106
6.4 106
4.8 106
3.2 106
1.6 106
0
1.6 106
3.2 106
4.8 106
6.4 106
8 106
POLINÓMIO CÚBICO: 2, 3 y 4 ELEMENTOS
0
M2 x1( )
M3 x( )
M4 x4( )
x1 x x4
133
Introducción al método de los elementos finitos aplicando MATHCAD
Carga distribuida triangular
En este punto se va a revisar cómo tratar un problema que presen-te carga distribuida triangular, esta distribución de carga se lo encuen-tra en los elementos sometidos a presión hidrostática como por ejemplo en puertas de exclusa.
Figura 3.18.Viga carga distribuida triangular
La estrategia que se sigue es utilizar las integrales de igualación de trabajo y empezamos calculando las fuerzas y momentos para el ele-mento finito 1
La carga distribuida es: (3.34)
2 0 d2y 0 1 1 d1y 0
0
L
xwx
2 L
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
dL2 w60
(3.35)
134
Análisis por el método de elementos finitos de elementos a flexión. Vigas
1 0 2 1 d1y 0 d2y 0
0
L
xwx
2 L
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
dL2 w40
(3.36)
d1y 1:= �2 0:= d2y 0:= �1 0:=
0
L
xwx
2 L
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
d3 L w
40
(3.37)
�1 0:= �2 0:= d1y 0:= d2y 1:=
0
L
xwx
2 L
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
d7 L w
40
(3.38)
Calculo de las fuerzas y momentos para el elemento finito 2La carga distribuida para el tramo 2 es: (3.39)
�1 1:= �2 0:= d1y 0:= d2y 0:=
0
L
xwx
2 L
w2
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
d7 L2 w120
(3.40)
�1 0:= �2 1:= d1y 0:= d2y 0:=
0
L
xwx
2 L
w2
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
dL2 w15
(3.41)
135
Introducción al método de los elementos finitos aplicando MATHCAD
�1 0:= �2 0:= d1y 1:= d2y 0:=
0
L
xwx
2 L
w2
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
d13 L w
40
(3.42)
�1 0:= �2 0:= d1y 0:= d2y 1:=
0
L
xwx
2 L
w2
2
L3d1y d2y( )
1
L21 2( )
x3
3
L2
d1y d2y( )1L
2 1 2( )
x2 1 x d1y
d17 L w
40
(3.43)
En base del siguiente diagrama se adicionaran las fuerzas y mo-mentos equivalentes por tanto:
Efectuando la suma correspondiente se obtiene:
136
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Para dos elementos iguales se genera la matriz de rigidez:
K
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
4 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
24 E I
L3
0
12 E I
L3
6 E I
L2
6 E I
L2
2 E IL
0
8 E IL
6 E I
L2
2 E IL
0
0
12 E I
L3
6 E I
L2
12 E I
L3
6 E I
L2
0
0
6 E I
L2
2 E IL
6 E I
L2
4 E IL
Las condiciones de frontera son por lo tanto:
d 0 0 d2y 2 d3y 3( ) d2y
F F1y M1 wL2
w L2
3017 w L
40w L2
15
F1yF1y
La ecuación a resolver
FT K dT solve
F1y
M1
d2y
2
d3y
3
37 L w
4079 L2
w60
121 L4 w
240 E I
41 L3 w
48 E I
22 L4 w
15 E I
L3 wE I
137
Introducción al método de los elementos finitos aplicando MATHCAD
La sustitución inversa da los verdaderos valores de la reacción:
K
0
0
121 L4 w
240 E I
41 L3 w
48 E I
22 L4 w
15 E I
L3 wE I
3 w L40
w L2
60
w L2
w L2
30
17 w L40
w L2
15
simplify
L w
4 L2 w
3
0
0
0
0
La deformada se la efectúa con las condiciones conocidas:
d1y 0:= �1 0:= d2y121 L4
× w×240 E× I×
�:= �241 L3
× w×48 E× I×
�:= �3L3 w×E I×
�:= d3y22 L4
× w×15 E× I×
�:=
FEA x( ) FEA1 x( ) 0 x Lif
FEA2 x( ) L x 2 Lif
ymax x( )
1E I
w x5
240 Lw L x3
6
2 w L( )2 x2
3
Figura 3.18.Deformada carga distribuida triangular
0 100 200 300 400 500 600 700 800 900 1 103�60�54�48�42�36�30�24�18�12�6�0
Deformada0
�FEA x( )
ymax x( )
x
Fuente propia, programa MathCAD
En la deformada se expresa una coincidencia total
138
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Cálculo de vigas con elementos de alto orden (hP-FEM)5
Los textos que abordan la teoría de elementos finitos en el aparta-do correspondiente a vigas, utiliza una función cúbica para describir el desplazamiento. Esta función genera resultados exactos cuando las cargas son puntuales, pero cuando la carga es distribuida existe una incompatibilidad en el momento ya que al ser este igual a la segunda derivada del desplazamiento se generara una representación lineal del mismo. Esta dificultad es solventada conforme se aumentan el núme-ro de nodos o lo que es lo mismo discretizando la viga, conforme se incremente el número de elementos finitos los resultados obtenidos serán más cercanos a los resultados que predice la teoría de vigas. La alternativa que se propone es incrementar el orden de la función de desplazamiento mediante la inclusión de un nodo interno. Por lo tanto la propuesta es desarrollar una función de alto orden de modo que con apenas uno o dos elementos finitos dependiendo del proble-ma, se alcance mejor exactitud que con el método tradicional de so-lución (polinomio cúbico). Los pasos a seguir serán básicamente pro-poner una nueva función de desplazamiento incrementando el grado del polinomio, derivar las funciones de interpolación , desarrollando una matriz de rigidez que incrementará la exactitud de cualquier pro-blema de vigas. Adicionalmente estos problemas se pueden resolver ventajosamente con las herramientas de cálculo simbólico del software MathCAD proporcionando un contexto pedagógico muy amplio para la explicación conceptual del método de los elementos finitos.
Obtención de las funciones de interpolación Se propone ahora utilizar un elemento de alto orden con un nodo
interno. Cuando se incrementa el orden del polinomio se dice que el método de refinamiento es “p-refinement”. Se sabe por lo tanto que el elemento finito tendrá 6 grados de libertad, 2 por nodo, sujeto a fuerzas transversales y momentos (Figura 4) y que el nodo interno no se conec-ta con ningún elemento.
5 José F. Olmedo, ANÁLISIS DE VIGAS POR EL MÉTODO DE ELEMENTOS FINITOS CON ELEMENTOS DE ALTO ORDEN PARA RESOLVER PROBLEMAS DE CARGA DISTRIBUIDA, Congreso COLIM 2014.
139
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 3.19.Elemento finito propuesto
Fuente propia
Si se tienen 6 grados de libertad el polinomio de interpolación de-berá ser de grado 5 de la forma siguiente:
(3.44)
(3.45)
Las condiciones de frontera de acuerdo a la Figura 4 son:
140
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Y reemplazando y ordenando en forma matricial se obtiene:
0
0
L2
5
5L2
4
L5
5 L4
0
0
L2
4
4L2
3
L4
4 L3
0
0
L2
3
3L2
2
L3
3 L2
0
0
L2
2
2L2
L2
2 L
0
1
L2
1
L
1
1
0
1
0
1
0
a1
a2
a3
a4
a5
a6
d1y
1
d2y
2
d3y
3
(3.46)
Se logra obtener el vector columna de incógnitas mediante:
[K]*[a]=[d] → [a]=[K]^(-1)*[d] (3.47)
24
L5
68
L4
66
L3
23
L2
0
1
4
L4
12
L3
13
L2
6L
1
0
0
16
L4
32
L3
16
L2
0
0
16
L4
40
L3
32
L2
8L
0
0
24
L5
52
L4
34
L3
7
L2
0
0
4
L4
8
L3
5
L2
1L
0
0
d1y
1
d2y
2
d3y
3
24 d1y
L5
24 d3y
L5
4 1
L4
16 2
L4
4 3
L4
16 d2y
L468 d1y
L4
52 d3y
L4
12 1
L3
40 2
L3
8 3
L3
66 d1y
L3
32 d2y
L3
34 d3y
L3
13 1
L2
32 2
L2
5 3
L2
16 d2y
L2
23 d1y
L2
7 d3y
L2
6 1L
8 2
L
3L
1
d1y
El polinomio de interpolación se puede ensamblar multiplicando el vector fila (x5 x4 x3 x2 x 1) Por el resultado anterior [a]
(3.48)
141
Introducción al método de los elementos finitos aplicando MATHCAD
El resultado de esta operación genera un polinomio sumamente grande para ser desplegadas en el documento. Este polinomio se re-agrupa de tal manera que las funciones de interpolación Ni puedan mostrarse en forma explícita y por tanto la función de desplazamiento se representara en la siguiente forma:
(3.49)
Obteniendo el siguiente resultado para cada función de interpolación:
(3.50)
Para L=1 se pueden graficar las funciones de interpolación para verificar que tenga el valor de 1 en su nodo asociado y cero en los otros nodos, ver figura 3.20 y 3.21.
Figura 3.20. Funciones de interpolación de alto orden N1, N3, N5
Fuente propia
142
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Figura 3.21: Derivadas de las funciones de interpolación de alto orden N2, N4, N6
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 11.2
1
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
Derivas de las funciones de interpolación asociadas a las rotaciones
xN2 x( )d
d
xN4 x( )d
d
xN6 x( )d
d
x
Fuente propia
Obtención de la matriz de rigidezEl paso más importante del procedimiento es determinar la matriz
de rigidez para lo cual se seguirá la formulación dada por el método de Galerkin para vigas donde:
(3.51)
Para N1 Para N2 Para N3....
(3.52)
143
Introducción al método de los elementos finitos aplicando MATHCAD
Afortunadamente el programa MathCAD posee herramientas de cálculo simbólico que permite obtener fácilmente estas integrales, como ejemplo se despliega k11 y k12
0
L
x2x24
xL
5 68
xL
4- 66
xL
3+ 23
xL
2- 1+
d
d
2
2x24
xL
5 68
xL
4- 66
xL
3+ 23
xL
2- 1+
d
d
2
d5092
35 L3
0
L
x2x24
xL
5 68
xL
4- 66
xL
3+ 23
xL
2- 1+
d
d
2
2x4- x
xL
4 12 x
xL
3+ 13 x
xL
2- 6 x
xL
+ x-
d
d
2
d1138
35 L2
-
Obteniendo finalmente la matriz de rigidez y con ella la formula-ción de elemento finito dada por:
Carga distribuidaComo anteriormente se manifestó, los efectos de las cargas pun-
tuales son exactamente evaluados cuando se utiliza la función de des-plazamiento cúbica, mientras que cuando se utiliza carga distribuida necesariamente se debe discretizar la viga, por tal razón el enfoque será con respecto a esta carga. Las restricciones y las cargas son apli-cadas únicamente en los nodos. La aproximación usual es reemplazar
144
Análisis por el método de elementos finitos de elementos a flexión. Vigas
la carga distribuida con fuerzas nodales y momentos tal que el trabajo mecánico efectuado por las cargas nodales (19) sea igual al generado por la carga distribuida (18). El trabajo mecánico debido a la carga distribuida w está dado por:
(3.53)
Donde es la función desplazamiento. El trabajo generado por fuer-zas discretas para tres nodos estaría dado por:
(3.54)
Se igualan ambos trabajos según la expresión:
(3.55)
Y se evalúa la integral para cada desplazamiento arbitrario
(3.56)
(3.57)
145
Introducción al método de los elementos finitos aplicando MATHCAD
(3.58)
(3.59)
(3.60)
(3.61)
Las fuerzas nodales equivalentes obtenidas están representadas en la siguiente figura 3.22.
Figura 3.22: Fuerzas discretizadas equivalentes
Fuente propia
146
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Ejemplo de implementación en MathcadSeguidamente se realiza un ejemplo con el fin obtener la solución
por el método de elementos finitos en base a la nueva matriz obtenida, así como comparar la solución tanto con la teoría de vigas como con los resultados obtenidos en base del polinomio cúbico. Se realiza por lo tanto el cálculo de una viga de acero biempotrada con un solo elemento finito, sometida a una carga distribuida w =-100 N/mm, con una longi-tud L=1 m, de medidas 100 mm x 10 mm. Las propiedades físicas son: E = 200000 MPa, b=10, h=100, , I= 1/12 b h3 ver figura 3.23
Figura 3.23: Viga biempotrada sometida a carga distribuida
Fuente propia
Donde las condiciones de contorno son:
d1y=0 ϕ1=0 m2=0 f2y=-8 (L w)/15 d3y=0 ϕ3=0
F1y
M1
F2y
M2
F3y
M3
E I
5092
35 L3
1138
35 L2
512
5 L3
384
7 L2
1508
35 L3
242
35 L2
1138
35 L2
33235 L
128
5 L2
647 L
242
35 L2
3835 L
512
5 L3
128
5 L2
1024
5 L3
0
512
5 L3
128
5 L2
384
7 L2
647 L
0
2567 L
384
7 L2
647 L
1508
35 L3
242
35 L2
512
5 L3
384
7 L2
5092
35 L3
1138
35 L2
242
35 L2
3835 L
128
5 L2
647 L
1138
35 L2
33235 L
d1y
1
d2y
2
d3y
3
= solve d2y 2 F1y M1 F3y M3L4 w
384 E I 0
4 L w15
L2 w15
4 L w
15L2 w15
Como se puede observar la deflexión en el nodo 2 coincide con el valor teórico para la viga biempotrada. Se efectúa la substi-tución inversa.
147
Introducción al método de los elementos finitos aplicando MATHCAD
E I
5092
35 L3
1138
35 L2
512
5 L3
384
7 L2
1508
35 L3
242
35 L2
1138
35 L2
332
35 L
128
5 L2
64
7 L
242
35 L2
38
35 L
512
5 L3
128
5 L2
1024
5 L3
0
512
5 L3
128
5 L2
384
7 L2
64
7 L
0
256
7 L
384
7 L2
64
7 L
1508
35 L3
242
35 L2
512
5 L3
384
7 L2
5092
35 L3
1138
35 L2
242
35 L2
38
35 L
128
5 L2
64
7 L
1138
35 L2
332
35 L
0
0
L4 w
384 E I
0
0
0
4 L w
15
L2 w
15
8 L w
15
0
4 L w
15
L2 w
15
Para determinar las fuerzas y momentos se resta del producto de la substitución inversa las fuerzas equivalentes obtenidas en el punto. {f}= [K]{d}- {f_0 }
4 L w
15
L2 w
15
8 L w
15
0
4 L w
15
L2 w
15
7 L w
30
L2 w
60
8 L w
15
0
7 L w
30
L2 w
60
L w
2
L2 w
12
0
0
L w
2
L2 w
12
La deformada se la obtiene a través de la función de desplazamien-to transversal con las condiciones siguientes:
x 0 0.01 LN1 x( ) 24
x
L
5 68
x
L
4 66
x
L
3 23
x
L
2 1
N2 x( ) 4 xx
L
4 12 x
x
L
3 13 x
x
L
2 6 x
x
L
x
148
Análisis por el método de elementos finitos de elementos a flexión. Vigas
N3 x( ) 16x
L
4 32
x
L
3 16
x
L
2
N4 x( ) 16 xx
L
4 40 x
x
L
3 32 x
x
L
2 8 x
x
L
N5 x( ) 24x
L
5 52
x
L
4 34
x
L
3 7
x
L
2
N6 x( ) 4 xx
L
4 8 x
x
L
3 5 x
x
L
2
x2
L
νFEA(x)=N1 (x) d1y+N2 (x) ϕ1+N3 (x) d2y+N4 (x) ϕ2+N5 (x) d3y+N6 (x) ϕ3
Adicionalmente en el mismo gráfico se dibujara la curva elástica obtenida de la ecuación general de vigas: (3.62)
Pudiendo observarse en la figura 3.24 la total concordancia de ambas curvas
Figura 3.24: Comparativa curvas
Fuente propia
Igualmente se obtiene el Momento, figura 3.25 a lo largo del eje mediante:
149
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 3.25: Momento obtenido para la teoría propuesta
0 100 200 300 400 500 600 700 800 900 1 103�
1� 107�
8.4� 106�
6.8� 106�
5.2� 106�
3.6� 106�
2� 106�
4� 105�
1.2 106�
2.8 106�
4.4 106�
6 106�
MOMENTO
0
M x( )
x
Fuente propia
Y el diagrama de cortante figura 3.26 con
0 100 200 300 400 500 600 700 800 900 1 103�
6� 104�
4.8� 104�
3.6� 104�
2.4� 104�
1.2� 104�
01.2 104
�
2.4 104�
3.6 104�
4.8 104�
6 104�
CORTANTE
0Vcortante x( )
x
Figura 3.26: Cortante obtenido para la teoría propuestaFuente propia
Ejemplo de implementación IIEl segundo ejemplo, figura 3.27 consistirá en calcular una viga
de acero biempotrada con un elemento finito, sometida a una carga distribuida w =-100N/mm hasta la mitad de la longitud de la viga, de una longitud L=1 m, de medidas 100 mm x 10 mm. Las propiedades físicas son: E = 200000 MPa, b=10mm, h=100 mm, I= 1/12 b h3. Este problema puede realizarse de dos formas, utilizando 2 elementos fini-
150
Análisis por el método de elementos finitos de elementos a flexión. Vigas
tos, en ese caso la matriz de rigidez será de 10 x 10 o determinar por el método de igualación del trabajo generado por las cargas y utilizar un solo elemento finito, se va a optar por este último método
Figura 3.27: Viga biempotrada sometida a carga distribuida hasta la mitad de la viga
Fuente propia
Efectuando la igualación del trabajo mecánico según el punto se obtienen las fuerzas equivalentes siguientes según la figura 3.28
Figura 3.28: Fuerzas discretizadas equivalentesFuente propia
Donde las condiciones de contorno son:
F1y
M1
F2y
M2
F3y
M3
E I
5092
35 L3
1138
35 L2
512
5 L3
384
7 L2
1508
35 L3
242
35 L2
1138
35 L2
332
35 L
128
5 L2
64
7 L
242
35 L2
38
35 L
512
5 L3
128
5 L2
1024
5 L3
0
512
5 L3
128
5 L2
384
7 L2
64
7 L
0
256
7 L
384
7 L2
64
7 L
1508
35 L3
242
35 L2
512
5 L3
384
7 L2
5092
35 L3
1138
35 L2
242
35 L2
38
35 L
128
5 L2
64
7 L
1138
35 L2
332
35 L
d1y
1
d2y
2
d3y
3
= solve d2y 2 F1y M1 F3y M3L4 w
768 E I
7 L3 w
6144 E I
47 L w
240
7 L2 w
160
17 L w
240
11 L2 w
480
151
Introducción al método de los elementos finitos aplicando MATHCAD
Como se puede observar la deflexión en el nodo 2 en este caso es
E I
5092
35 L3
1138
35 L2
512
5 L3
384
7 L2
1508
35 L3
242
35 L2
1138
35 L2
332
35 L
128
5 L2
64
7 L
242
35 L2
38
35 L
512
5 L3
128
5 L2
1024
5 L3
0
512
5 L3
128
5 L2
384
7 L2
64
7 L
0
256
7 L
384
7 L2
64
7 L
1508
35 L3
242
35 L2
512
5 L3
384
7 L2
5092
35 L3
1138
35 L2
242
35 L2
38
35 L
128
5 L2
64
7 L
1138
35 L2
332
35 L
0
0
L4 w
768 E I
7 L3 w
6144 E I
0
0
47 L w
240
7 L2 w
160
4 L w
15
L2 w
24
17 L w
240
11 L2 w
480
Para determinar las fuerzas y momentos se resta del producto de la substitución inversa las fuerzas equivalentes obtenidas en el punto. {f}= [K]{d}- {f_0 }
4 L w
15
L2 w
15
8 L w
15
0
4 L w
15
L2 w
15
7 L w
30
L2 w
60
8 L w
15
0
7 L w
30
L2 w
60
L w
2
L2 w
12
0
0
L w
2
L2 w
12
Finalmente generamos las curvas figura 3.29 mediante:
vFEA(x)=N1 (x) d1y+N2 (x) ϕ1+N3 (x) d2y+N4 (x) ϕ2+N5 (x) d3y+N6 (x) ϕ3
152
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Figura 3.28: Curva elástica
Fuente propia
Igualmente se obtiene el Momento a lo largo del eje mediante:
Figura 3.29: Momento
0 100 200 300 400 500 600 700 800 900 1 103�
6� 106�
5.1� 106�
4.2� 106�
3.3� 106�
2.4� 106�
1.5� 106�
6� 105�
3 105�
1.2 106�
2.1 106�
3 106�
MOMENTO
0
M x( )
x
Fuente propia
Y el diagrama de cortante con , Figura 3.30
Figura 3.30: Cortante
0 100 200 300 400 500 600 700 800 900 1 103�
2� 104�
1.3� 104�
6� 103�
1 103�
8 103�
1.5 104�
2.2 104�
2.9 104�
3.6 104�
4.3 104�
5 104�
CORTANTE
0
Vcortante x( )
x
Fuente propia
153
Introducción al método de los elementos finitos aplicando MATHCAD
Viga de timoshenko
A medio camino entre una viga y una placa, existen un tipo de viga corta, ver Figura 3.31 en las cuales la teoría de Bernoulli genera resultados inexactos ya que en este caso no se pueden despreciar las deformaciones debida al cortante. Se considera principalmente que las secciones transversales normales permanecen planas pero no necesa-riamente perpendiculares al eje de la viga.6
Figura 3.31: Viga corta
Fuente propia
Para ejemplificar la resolución de este tipo de viga se procederá de la siguiente manera:
Datos:d1y 0:= �1 0:= m2 0:= f2y P�:= P
La matriz de rigidez para viga de Timoshenko es:
f1y
m1
f2y
m2
E I
L3 1 ( )
12
6 L
12
6 L
6 L
4 ( ) L2
6 L
2 ( ) L2
12
6 L
12
6 L
6 L
2 ( ) L2
6 L
4 ( ) L2
d1y
1
d2y
2
solve
f1y
m1
d2y
2
P L P4 L3 P L3 P
12 E I
L2 P2 E I
6 Marcelo Tulio Piovan, TENSIONES Y DEFORMACIONES REVISIÓN DE PRINCIPIOS FÍSICOS
154
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Donde F en la matriz de rigidez es el término de corrección por cor-te. Los parámetros de la viga son:
L 500:= P 1000:= E 200000:= x 0 0.1, L..:=
b 10:= h 160:= � 0.3:= I1
12b× h3
×:= ks 0.83:= A b h×:= G
E2 1 �+( )×
:=
Siendo ks el factor de corrección de Área.
gE I
ks A G6.683 103
12 g
L20.321 (3.63)
De la solución los desplazamientos son:
2L2 P2 E I
d2y4 L3 P L3 P
12 E I (3.64)
Las funciones de forma son
N1 x( ) 2 x3 3 x2
L L3 12 g
N2 x( ) x3 L 2 L2 6 g x2
x L3 6 g L
N3 x( ) 2 x3 3 x2
L 12 g
N4 x( ) x3 L x2 L2 6 g 6 g L x
N4 x( ) x3 L× x2 L2� 6 g×+( )×+ 6 g× L× x×�:= (3.65)
Y el desplazamiento por tanto es:
FEA x( )1
L L2 12 g N1 x( ) d1y
1
L L2 12 g N2 x( ) 1
1
L L2 12 g N3 x( ) d2y
1
L L2 12 g N4 x( ) 2
EDO x( )P L
2 E Ix2 x3
3 L
155
Introducción al método de los elementos finitos aplicando MATHCAD
Figura 3.32: Comparativa viga Timoshenko vs Bernoulli
0 100 200 300 400 5000.08�
0.06�
0.04�
0.02�
0
0.02Viga Timoshenko
�FEA x( )
�EDO x( )
500
x
Fuente propia
Resolución de vigas con ansys apdl
Determine las deflexiones y tensiones
Figura 3.33: Viga continúa 7
7 Amar Khennane, Introduction to Finite Element Analysis Using MATLAB® and Abaqus, pág. 90
156
Análisis por el método de elementos finitos de elementos a flexión. Vigas
1. En preferences seleccione structural
2. Crear la geometriaUsar: Preprocessor < Modeling < Create < Keypoints < In Active CSIngresar Keypoints de acuerdo a la tabla siguiente:
Tabla 3.2. Coordenadas keypoints
Keypoints x y
1 0 0
2 2000 0
3 4000 0
4 9000 0
5 16000 0
3. Clic en Apply para aceptar el resultado, El resultado es el siguiente:
157
Introducción al método de los elementos finitos aplicando MATHCAD
Para dibujar las líneas vamos a:Preprocessor < Modeling < Create < Lines < Lines < In Active Coord
Con el mouse vamos señalando los Keypoints hasta obtener el si-guiente dibujo:
158
Análisis por el método de elementos finitos de elementos a flexión. Vigas
4. Defina el tipo de elementoSeleccione: Element Type < Add/Edit/DeleteY seleccione Beam, 3node 189, correspondiente a una viga con ele-mento de alto orden
Definir la sección de la viga mediante: Sections < Beam < Common Sections
159
Introducción al método de los elementos finitos aplicando MATHCAD
5. Definir propiedades del materialEn Preprocessor seleccionar Material Props < Material Models. Se-leccionar material elástico, lineal e isotrópico, llenar en el campo de Módulo de elasticidad el valor de E = 200000 MPa y 0.3 para Poisson que corresponde al acero.
6. MalladoSeleccionar Meshing < Mesh Tool < Size Controls < Lines < Set
160
Análisis por el método de elementos finitos de elementos a flexión. Vigas
Seleccionar Meshing < Mesh Tool < Mesh < Pick All
Para verificar la dirección de las fuerzas en el command prompt escribir eshape,1
Luego con Plot < Elements
161
Introducción al método de los elementos finitos aplicando MATHCAD
7. Se observa que la dirección de las fuerzas debe ser en z, y se acti-van los Keypoints para localizar las fuerzas
En las líneas entre el Keypoints 3 y 4 colocamos una presión de -4 N/mm
Igualmente entre los keypoints de 5 a 6 colocamos -10 N/mm
162
Análisis por el método de elementos finitos de elementos a flexión. Vigas
El nodo 1 es totalmente restringido, por lo tanto vamos a:Define Loads < Apply < Structural < Displacement < On Keypoints
En los Keypoints 3,4,5 se colocan restricciones en Z
163
Introducción al método de los elementos finitos aplicando MATHCAD
Resolver el sistema de ecuaciones mediante:Solution < Solve < Current LSAparecerá el mensaje: Solution is done
8. Postprocesado revisión de resultadosGeneral Postproc < Plot Results < Deformed ShapePlot Results < Deformed Shape
La deformación es de 5.7444 mm y las tensiones son:
165
Introducción al método de los elementos finitos aplicando MATHCAD
Determinar la matriz de rigidez del portico
Las estructuras más complejas que involucran elementos a flexión se denominan pórticos. Los cuales al igual que las armaduras se re-suelven utilizando matrices de transformación, en este caso se debe hallar la matriz de rigidez del elemento viga orientada en el plano con un ángulo arbitrario. Luego se incluye el grado de desplazamiento axial en la matriz de rigidez local.
En este caso el elemento finito corresponde a una viga de inclina-ción arbitraria con tres grados de libertad por nodo, ver figura 4.2.
Figura 4.1: Pórtico de taller
Figura 4.2: Elemento finito viga inclinada
168
Análisis por el método de elementos finitos. Pórticos planos
Dada la ecuación de transformación conocida:
En las ecuaciones de transformación se integra el término corres-pondiente a la rotación, el cual no es influenciado por la orientación de la viga.
(4.1)
Donde la matriz de transformación T es por tanto
(4.2)
La matriz local de rigidez correspondiente a la viga se expande se-gún la expresión
(4.3)
169
Introducción al método de los elementos finitos aplicando MATHCAD
Esta matriz se debe combinar con la matriz de efectos axiales:
(4.4)
Cuya combinación genera el siguiente arreglo matricial:
k
C1
0
0
C1
0
0
0
12 C2
6 C2 L
0
12 C2
6 C2 L
0
6 C2 L
4 C2 L2
0
6 C2 L
2 C2 L2
C1
0
0
C1
0
0
0
12 C2
6 C2 L
0
12 C2
6 C2 L
0
6 C2 L
2 C2 L2
0
6 C2 L
4 C2 L2
C1C1
(4.5)
Donde C1 = A E/L y C2 = E I/L3. Utilizando la fórmula conocida:
(4.6)
Se obtiene la matriz de rigidez global:
K TT k T T
K simplify
E A C2 L2
12 I S2
L3
C E S 12 I A L2
L3
6 E I S
L2
E 12 I S2 A C2
L2
L3
C E S A L2 12 I
L3
6 E I S
L2
C E S 12 I A L2
L3
E 12 I C2 A L2
S2
L3
6 C E I
L2
C E S A L2 12 I
L3
E A L2 S2
12 C2 I
L3
6 C E I
L2
6 E I S
L2
6 C E I
L2
4 E IL
6 E I S
L2
6 C E I
L2
2 E IL
E A C2 L2
12 I S2
L3
C E S 12 I A L2
L3
6 E I S
L2
E A C2 L2
12 I S2
L3
C E S 12 I A L2
L3
6 E I S
L2
C E S 12 I A L2
L3
E 12 I C2 A L2
S2
L3
6 C E I
L2
C E S 12 I A L2
L3
E 12 I C2 A L2
S2
L3
6 C E I
L2
6 E I S
L2
6 C E I
L2
2 E IL
6 E I S
L2
6 C E I
L2
4 E IL
(4.7)
170
Análisis por el método de elementos finitos. Pórticos planos
Ejercicios con mathcad
Para el pórtico mostrado resolver los desplazamientos y rotaciones involucrados
Figura 4.3: Elemento finito viga inclinada
Donde la longitud es igual para todas las barras y mide 1000 mm, módulo de elasticidad 200000 MPa, A = 1000 mm2, I = 100000 mm4.
R1A EL
R2E I
L3
T
C
S
0
0
0
0
S
C
0
0
0
0
0
0
1
0
0
0
0
0
0
C
S
0
0
0
0
S
C
0
0
0
0
0
0
1
SS
k
R1
0
0
R1
0
0
0
12 R2
6 R2 L
0
12 R2
6 R2 L
0
6 R2 L
4 R2 L2
0
6 R2 L
2 R2 L2
R1
0
0
R1
0
0
0
12 R2
6 R2 L
0
12 R2
6 R2 L
0
6 R2 L
2 R2 L2
0
6 R2 L
4 R2 L2
K TT k T T
K simplify
200000 C2 240 S2
199760 C S
120000 S
200000 C2 240 S2
199760 C S
120000 S
199760 C S
240 C2 200000 S2
120000 C
199760 C S
240 C2 200000 S2
120000 C
120000 S
120000 C
80000000
120000 S
120000 C
40000000
200000 C2 240 S2
199760 C S
120000 S
200000 C2 240 S2
199760 C S
120000 S
199760 C S
240 C2 200000 S2
120000 C
199760 C S
240 C2 200000 S2
120000 C
120000 S
120000 C
40000000
120000 S
120000 C
80000000
171
Introducción al método de los elementos finitos aplicando MATHCAD
Elemento 1: La matriz de rigidez del elemento vertical Θ= 90
90
180
C cos ( ) S sin ( )
K1
400000 C2 1200 S2
398800 C S
600000 S
400000 C2 1200 S2
398800 C S
600000 S
398800 C S
1200 C2 400000 S2
600000 C
398800 C S
1200 C2 400000 S2
600000 C
600000 S
600000 C
400000000
600000 S
600000 C
200000000
400000 C2 1200 S2
398800 C S
600000 S
400000 C2 1200 S2
398800 C S
600000 S
398800 C S
1200 C2 400000 S2
600000 C
398800 C S
1200 C2 400000 S2
600000 C
600000 S
600000 C
200000000
600000 S
600000 C
400000000
Con las siguientes sentencias augment y stack se expande la ma-triz de rigidez
zeros
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
z1 augment zerosT zerosT z2 augment z1 z1( )
K1 stack K1 z2( ) K1 augment K1 z1( )
K1
1200
0
600000
1200
0
600000
0
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
600000
0
400000000
600000
0
200000000
0
0
0
0
0
0
1200
0
600000
1200
0
600000
0
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
600000
0
200000000
600000
0
400000000
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
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
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Elemento 2: La matriz de rigidez del elemento horizontal Θ= 0
0
180
C cos ( ) S sin ( )
172
Análisis por el método de elementos finitos. Pórticos planos
K2
400000 C2 1200 S2
398800 C S
600000 S
400000 C2 1200 S2
398800 C S
600000 S
398800 C S
1200 C2 400000 S2
600000 C
398800 C S
1200 C2 400000 S2
600000 C
600000 S
600000 C
400000000
600000 S
600000 C
200000000
400000 C2 1200 S2
398800 C S
600000 S
400000 C2 1200 S2
398800 C S
600000 S
398800 C S
1200 C2 400000 S2
600000 C
398800 C S
1200 C2 400000 S2
600000 C
600000 S
600000 C
200000000
600000 S
600000 C
400000000
zeros
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
z1 augment zeros zeros( )
K2 augment zerosT K2 zerosT K2 stack z1 K2 z1( )
K2
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
0
0
0
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
0
0
1200
600000
0
1200
600000
0
0
0
0
0
0
0
600000
400000000
0
600000
200000000
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
0
0
1200
600000
0
1200
600000
0
0
0
0
0
0
0
600000
200000000
0
600000
400000000
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
0
0
0
0
0
0
0
0
0
Elemento 3: La matriz de rigidez del elemento horizontal Θ= 270
270
180
S sin ( ) C cos ( )
K3
400000 C2 1200 S2
398800 C S
600000 S
400000 C2 1200 S2
398800 C S
600000 S
398800 C S
1200 C2 400000 S2
600000 C
398800 C S
1200 C2 400000 S2
600000 C
600000 S
600000 C
400000000
600000 S
600000 C
200000000
400000 C2 1200 S2
398800 C S
600000 S
400000 C2 1200 S2
398800 C S
600000 S
398800 C S
1200 C2 400000 S2
600000 C
398800 C S
1200 C2 400000 S2
600000 C
600000 S
600000 C
200000000
600000 S
600000 C
400000000
zeros
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
z1 augment zerosT zerosT z2 augment z1 z1( )
173
Introducción al método de los elementos finitos aplicando MATHCAD
K3 augment z1 K3( )
K3 stack z2 K3( )
K3
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
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
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1200
0
600000
1200
0
600000
0
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
600000
0
400000000
600000
0
200000000
0
0
0
0
0
0
1200
0
600000
1200
0
600000
0
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
600000
0
200000000
600000
0
400000000
La suma de las matrices individuales da:
K
1200
0
600000
1200
0
600000
0
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
600000
0
400000000
600000
0
200000000
0
0
0
0
0
0
1200
0
600000
401200
0
600000
400000
0
0
0
0
0
0
400000
0
0
401200
600000
0
1200
600000
0
0
0
600000
0
200000000
600000
600000
800000000
0
600000
200000000
0
0
0
0
0
0
400000
0
0
401200
0
600000
1200
0
600000
0
0
0
0
1200
600000
0
401200
600000
0
400000
0
0
0
0
0
600000
200000000
600000
600000
800000000
600000
0
200000000
0
0
0
0
0
0
1200
0
600000
1200
0
600000
0
0
0
0
0
0
0
400000
0
0
400000
0
0
0
0
0
0
0
600000
0
200000000
600000
0
400000000
Condiciones de Frontera
d1x 0 d1y 0 1 0
d4x 0 d4y 0 4 0
f2x 10000 f2y 0 m2 0
f3x 0 f3y 0 m3 5000
La submatriz se extrae
KS submatrixK 3 8 3 8( )
174
Análisis por el método de elementos finitos. Pórticos planos
KS
401200
0
600000
400000
0
0
0
401200
600000
0
1200
600000
600000
600000
800000000
0
600000
200000000
400000
0
0
401200
0
600000
0
1200
600000
0
401200
600000
0
600000
200000000
600000
600000
800000000
Los desplazamientos se determinan de acuerdo a la siguiente expresión:
d KS 1
f2x
f2y
m2
f3x
f3y
m3
d
5.966
0.011
3.597 10 3
5.954
0.011
3.576 10 3
Se extraen las soluciones mediante:
d2x d0 d3x d3
d2y d1 d3y d4
2 d2 3 d5
Postprocesado: Las fuerzas se obtienen mediante la sustitución inversa
FTOTAL K
d1x
d1y
1
d2x
d2y
2
d3x
d3y
3
d4x
d4y
4
175
Introducción al método de los elementos finitos aplicando MATHCAD
Obteniéndose los siguientes resultados:
Los desplazamientos se los determina mediante: factor := 30
Lx
0
0
1000
1000
Ly
0
1000
1000
0
Lxc
0 d1x factor
0 d2x factor
1000 d3x factor
1000 d4x factor
Lyc
0 d1y factor
1000 d2y factor
1000 d3y factor
0 d4y factor
Figura 4.3: Elemento finito viga inclinada
100� 60 220 380 540 700 860 1.02 103� 1.18 103
� 1.34 103� 1.5 103
�100�
60
220
380
540
700
860
1.02 103�
1.18 103�
1.34 103�
1.5 103�
Desplazamiento
Ly
Lyc
Lx Lxc,
Fuente propia
176
Análisis por el método de elementos finitos. Pórticos planos
Resolución de pórticos con ANSYS APDL
Diseñe una grua de portico para 10 toneladas y una luz de 6000 mm
Grua de pórtico8
Crear Keypoints en 0,0 y 6000,0Crear la linea respectivaSeleccionar una primera alternativa
8 http://pro-grua.blogspot.com/2010/05/gruas-portico.html
177
Introducción al método de los elementos finitos aplicando MATHCAD
Que podría ser la viga IPN 450Cuyas medidas seria h= 450, b= 170, tw = 16,2 y tf = 24.3 mmEscogemos Elemento beam 188 con options k3 = cuadratic
179
Introducción al método de los elementos finitos aplicando MATHCAD
Propiedades del acero
Se malla con una magnitud de 250 mm
180
Análisis por el método de elementos finitos. Pórticos planos
En la línea de comandos ponemos /ESHAPE,1 y plot elements
Se puede ver la orientación para los otros Keypoints.
183
Introducción al método de los elementos finitos aplicando MATHCAD
Dibujar las 8 lineas restantes
184
Análisis por el método de elementos finitos. Pórticos planos
Para los miembros inferiores se escoge tubo sin costura de 4 plg ( 114.3 x 8.56 mm)
Cargas:10 tons = 100000 N
185
Introducción al método de los elementos finitos aplicando MATHCAD
RestriccionesUn apoyo fijo y los 3 restantes solo restricción en
186
Análisis por el método de elementos finitos. Pórticos planos
17 mm de desplazamiento Tensiones da: 71.9 MPa, N = 210/71.9 = 3
187
Introducción al método de los elementos finitos aplicando MATHCAD
Lista de referencias
Logan, Daryl, (2012), Fifth Edition, CENGAGE LEARNING, Stamford. A FIRST COURSE IN THE FINITE ELEMENT METHOD, [UN PRIMER CURSO EN EL MÉTODO DE ELEMENTOS FINITOS].
Hutton, David, (2002), Second. Edition, McGraw-Hill, New York. FUNDAMENTAL OF FINITE ELEMENT ANALYSIS, ,[FUNDAMENTOS DEL ANÁLISIS POR ELE-MENTOS FINITOS]
Reedy, J, (2005), Third. Edition, McGraw-Hill, New York, AN INTRODUCTION TO THE FINITE ELEMENT METHOD, [INTRODUCCIÓN AL MÉTODO DE ELE-MENTOS FINITOS].
Amar Khennane, (2013), CRC Press, INTRODUCTION TO FINITE ELEMENT ANALY-SIS USING MATLAB® AND ABAQUS