Fundamentos del Método de Elemento Finito
Fundamentos del Método de Elemento Finito
de Elemento Finito
Primera Edición
Jaime G. Molina P.
Catedrático Emérito
Ingenierı́a Mecánica
U. M. S. A.
Procesado Procesado Procesado Procesado
mediante mediante mediante mediante
LATEX 2ε LATEX 2ε LATEX 2ε LATEX 2ε
Procesado Procesado Procesado Procesado
mediante mediante mediante
La presentación y disposición en conjunto de:
“ fundamentos del método de elemento finito ”
mediante
incluye materia preliminar, cuerpo de texto y materia posterior.
LATEX 2ε LATEX 2ε LATEX 2ε
Además también incluye un cd en la contratapa, que contiene
archivos electrónicos complementarios al documento.
LATEX 2ε
La propiedad intelectual de este trabajo se transfiere a la U.M.S.A.
Procesadocopiar,La distribuir
Procesado
Institución patrocinante y el autorProcesado
conceden licencia para
y/o modificar este documento según las condiciones
Procesado
mediante mediante mediante
de Licencia de Documentación Libre GNU, Versión 1.2 o posterior;
sin las Secciones Invariantes del documento, ningún texto de las mediante
Páginas Preliminares, y ningún texto de Materia Posterior.
Jaime G. Molina P.
Verano de 2010
Fundamentos del Método
de Elemento Finito
Primera Edición
Jaime G. Molina P.
Catedrático Emérito
Ingenierı́a Mecánica
U. M. S. A.
Índice de Contenido
1. Conceptos básicos 1
1.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Concepción del método de elemento finito ? . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1. Convergencia del método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.2. El método de diferencias finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3. Procedimiento general . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1. Pre–procesamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.2. Procesamiento o Solución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.3. Post–procesamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4. Breve historia del método de elemento finito . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.5. Ejemplos de análisis de elemento finito . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6. Objetivos del documento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2. Análisis uni–dimensional 17
2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2. El resorte linealmente elástico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.1. Ensamble del sistema en coordenadas globales . . . . . . . . . . . . . . . . . . . . 21
2.3. El elemento barra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4. Energı́a de deformación elástica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.4.1. Primer teorema de Castigliano . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.5. Energı́a potencial mı́nima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.6. Resumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
Problemas propuestos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
v
vi
4. Elementos de f lexión 85
4.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.2. Teorı́a elemental de vigas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.3. El elemento viga . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.4. Matriz de rigidez de elemento viga . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.5. Vector de carga de elemento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.6. Cargas nodales equivalentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.7. El elemento marco . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
4.8. El elemento rectilı́neo general . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
4.9. Comentarios finales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Problemas propuestos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Prefacio
La presente obra: “Fundamentos del Método de Elemento Finito”, tiene pretensión de ser
libro de texto para un curso terminal de pre–grado o inicial de nivel–superior (post–grado), diseñado
para ser incluı́do en los programas de estudio de ingenierı́a. La orientación temática de su contenido
se adapta muy apropiadamente a los programas de estudio de ingenierı́a mecánica e ingenierı́a civil,
aunque podrı́a ser considerado también como libro de consulta en otros campos de la ingenierı́a o en
un curso de matemática aplicada.
El método de elemento finito es una técnica de análisis y diseño ampliamente utilizada actual-
mente, que es esencial para los estudiantes de ingenierı́a que poseen rudimentos básicos de la teorı́a y
aplicaciones de los métodos discretos de análisis. Persiguiendo el objetivo de aprender algo más sobre
el método, a finales del año 2006 en la Universidad de Calghary — en Canadá — fuı́ invitado para
elaborar y enseñar una serie de “temas especiales” como parte de un simposio de charlas complemen-
tarias al curso extraordinario de modalidad electiva del método de elemento finito, implementado el
segundo semestre de dicho año en la universidad mencionada. El curso estaba compuesto de aproxi-
madamente tres–cuartas partes de teorı́a y una cuarta parte de uso de software comercial destinado
a la resolución de los problemas asignados a lo largo del curso. Subsecuentemente, desde dicha época,
el curso se ha convertido en uno regularmente ofrecido como asignatura de especialidad de modalidad
obligatoria en el programa de ingenierı́a mecánica, y generalmente posee una demanda muy alta de
estudiantes interesados en aprender esta materia (incluso como asignatura electiva, para ser cursada
por estudiantes pertenecientes a otras carreras de ingenierı́a).
Durante el proceso de desarrollo para el curso, los responsables de impartirlo nunca estuvimos
satisfechos con cualquier texto que se usó en el mismo, y nosotros probamos muchos libros escritos
sobre el tema, por cierto. Encontramos que los textos disponibles se ubicaban a un extremo u otro
de nuestras pretensiones; a saber, en algunos esencialmente ninguna teorı́a y mucha aplicación de
software; o toda la teorı́a y ninguna aplicación de software. Los planteamientos pedagógicos anteriores,
en mi opinión, representan: el primero, entrenamiento extensivo usando programas computacionales
(adecuado para quien domina la teorı́a y se inmiscuye exclusivamente con las aplicaciones del método);
mientras que el segundo pretende desarrollo extensivo teórico sin acercamiento al software asociado al
método (adecuado para quien pretende dominar la teorı́a, sin inmiscuirse en la aplicación práctica).
xii
En esa época esperaba que la experiencia que relato en anteriores párrafos se replicara en la uni-
versidad donde imparto enseñanza. Y afortunadamente, en la actualidad, por la apertura del ciclo de
Postgrado con el Programa de Maestrı́a en Ciencias de la Ingenierı́a Mecánica, tengo la oportunidad
de hacer el intento de replicar de alguna manera la experiencia académica antes comentada.
Es pretensión de éste documento ubicarse en lo posible en un punto medio donde ciertamente
se dé preferencia a la elaboración teórica, pero sin omitir los aspectos de tipo práctico mediante
utilización de software adecuado para ello. Esto porque consideramos este documento como un libro
de texto adecuado a un curso introductorio al método de elemento finito. En tal virtud, los problemas
propuestos para ser resueltos computacionalmente no son de un gran volumen de datos de entrada, de
modo que para resolverlos podrı́a utilizarse incluso un paquete de caracter académico (existen muchos
de descarga gratuita en diversos portales de la red Internet).
Pedagógicamente, yo creo que proporcionar entrenamiento práctico al estudiante pre–graduado en
el uso de un paquete de software particular, sin proporcionarle conocimiento de la concepción teórica
subyascente en los programas computacionales que conforman el paquete, es un perjuicio al estudiante
y puede ser hasta peligroso para sus futuros empleadores. Me considero que soy agudamente consciente
que la mayorı́a de los programas de estudio en ingenierı́a, tienen un paquete especı́fico de software de
elemento finito disponible para el uso del estudiante; y simultáneamente yo no creo que el texto que
los estudiantes usen sólo deba exclusivamente atarse en su desarrollo teórico a ese especı́fico software.
Por consiguiente, he escrito este texto para que sea considerado independiente de cualquier software
en particular. Doy énfasis a la teorı́a básica del método de elemento finito, en un contexto que pueda
entenderse por cualquier estudiante común de ingenierı́a, y dejo las porciones especı́ficas de utilización
de software a la curiosidad e inquietud del lector, para que él mismo complemente su aprendizaje
incorporando la parte práctica de utilización de ésta herramienta computacional de manera individual
y autodidacta.
Como el texto está elaborado para un curso introductorio al método de elemento finito, los requisitos
de conocimiento previo comprenden: los principios de la mecánica general en sus partes componen-
tes estática y dinámica, la mecánica de materiales, la mecánica de fluidos, la termodinámica y la
transmisión de calor, y finalmente el cálculo infinitesimal en concepción de las operaciones básicas de
integración y derivación de funciones, juntamente con la aplicación a la solución de las ecuaciones di-
ferenciales. Por necesidad, y cuando sea ineludible, se introducen ecuaciones diferenciales en derivadas
parciales al interior del desarrollo teórico presentado; pero de una manera que pueda entenderse sin
mayor dificultad, basado en los requisitos previos declarados.
Las aplicaciones del método de elemento finito a la transmisión del calor y la mecánica del medio
fluido son incluidas en Capı́tulos especı́ficos, pero las deducciones necesarias son tales que los cursos
anteriores en esos temas se requieren solamente por los conceptos fundamentales de definición que fue-
ron supuestamente elaborados en dichos cursos precedentes. Seguramente, muchos estudiantes habrán
tomado cursos regulares previos de transferencia de calor y mecánica de fluidos, y los temas aquı́ pre-
sentados pueden suponerse que resultarán ser extensión natural de aquellos temas fundamentales que
los estudiantes poseen como conocimiento previo a este curso.
Advierto que esta obra no es completamente inédita en autorı́a; mucha parte de su contenido es
traducción e interpretación libre del texto: “Fundamentals of the Finite Element Analysis” de David D.
Hutton (McGraw Hill – Higher Education Publications, 1st. Ed.,
2004 c , International Edition), con la
autorización del autor para efectuar una traducción no literal de las partes de interés. Otras partes de
desarrollo teórico sobretodo, son también traducciones no–literales de diversos documentos (artı́culos
y partes menores de otros libros descargados desde la red Internet), y una parte algo menor es de
autorı́a propia en concepción y redacción, la cual se ha incorporado en forma de párrafos intermedios
de explicación más detallada de los temas que son tratados.
El contenido del presente documento, seccionado siguiendo un formato estandar pre–establecido,
incluye los temas siguientes:
Capı́tulo 1 — es una introducción general al método del elemento finito, e incluye una descripción
del concepto básico de dividir un dominio en una serie de sub–dominios de tamaño finito. El
xiii
método de diferencias finitas se presenta para comparación con el método de elemento finito. Un
procedimiento general en la secuencia de definición, solución, y la interpretación de resultados
asociados con el modelo matemático de análisis se discute, y se relaciona a la serie generalmente
aceptada de: pre–procesamiento, solución, y post–procesamiento. También se incluye una historia
breve del método de elemento finito, como también algunos ejemplos que ilustran la aplicación
del método.
Capı́tulo 2 — introduce el concepto de matriz de rigidez de elemento finito, y la ecuación del campo
de desplazamientos asociado, en términos de las denominadas funciones de interpolación, usando
el resorte lineal como un elemento finito. Este elemento mecánico (el resorte lineal elástico) es
conocido por la mayorı́a de estudiantes pre–graduados, por lo que su comportamiento mecánico
no debe ser nuevo para ellos. Sin embargo, la representación del resorte lineal como un elemento
finito es nueva, y la virtud de éste enfoque es que proporciona un ejemplo simple y conciso
del método de elemento finito. La premisa de formulación del elemento resorte es extendida
al elemento barra, y se introducen los métodos de concepción energética. El primer teorema de
Castigliano es aplicado, como representación del principio de energı́a potencial mı́nima. El teorema
de Castigliano es un método simple para introducir al estudiante pre–graduado a los principios
de minimización funcional sin el uso explı́cito del cálculo variacional.
Capı́tulo 3 — utiliza el elemento barra del Capı́tulo anterior para ilustrar el ensamble global de las
ecuaciones de equilibrio para una estructura compuesta de muchos elementos finitos. Se desarrolla
la transformación de las ecuaciones gobernantes de comportamiento mecánico desde las coorde-
nadas del elemento hacia las coordenadas globales (de la estructura), y se ilustra el procedimiento
con ejemplos bi y tri–dimensionales. El método directo de rigidez se utiliza y se presentan dos
métodos para el ensamble de la matriz de rigidez global. Se discute la aplicación de condiciones
de borde lı́mite y la solución de las ecuaciones de restricción. Se muestra también el uso de la
solución básica del desplazamiento para obtener la tensión interna del elemento, y se muestra a
la tensión como una importante operación de la etapa de post–procesamiento.
Capı́tulo 4 — introduce el elemento viga como un elemento flexible y temáticamente como puente
a los requisitos de continuidad para los elementos de orden–superior. Se introduce el concepto
de continuidad de la derivada espacial y esto requiere un ajuste a las funciones de interpolación
asumidas, para asegurar la continuidad requerida por la solución. Se discuten los vectores de carga
nodales en el contexto de cargas discretas y distribuı́das, usando el método de equivalencia del
trabajo mecánico.
Los Capı́tulos 2, 3, y 4, introducen los procedimientos básicos de elemento–finito modelando
los sistemas considerados en el contexto de elementos estructurales simples que deben ser muy
conocidos por el estudiante desde el requisito previo de su curso de mecánica de materiales. Ası́, el
énfasis en la parte inicial del curso en que el texto es usado puede estar enmarcado en la aplicación
del método de elemento finito hacia situaciones simples, sin la introducción de nuevos conceptos
fı́sicos. Los elementos barra y viga pueden usarse para proporcionar al estudiante problemas
prácticos de estructuras reticuladas compuestas de elementos rectilı́neos, cuya solución deba ser
hallada aplicando el software disponible para el método de elemento finito. Cuando el alumno
adquiera dominio (en el contexto bidimensional) de manejo de los elementos barra y viga, este
conocimiento proporciona un ámbito adecuado para que se pueda combinar este desarrollo teórico
y generar el elemento marco, adecuado para el análisis o diseño de estructuras tipo pórtico.
Capı́tulo 5 — es el trampolı́n hacia los conceptos más avanzados de análisis por aplicación del método
de elemento finito. El método de residuos ponderados se introduce como la técnica fundamental
que será utilizada en el resto del texto. El método de Galerkin se utiliza exclusivamente desde
que particularmente he encontrado que este método es comprensible para el estudiante, y es dócil
a la formulación de una gama amplia de problemas que se presentan en los diversos campos de
la ingenierı́a. El material en este Capı́tulo repite los desarrollos de los elementos barra y viga
presentados anteriormente desde una perspectiva diferente, pero también extiende el concepto de
xiv
método básico presentado para hallar la respuesta dinámica; de modo que una cantidad conside-
rable de material del texto se consagra a la determinación de los modos naturales de vibración,
la ortogonalidad de los mismos, y la superposición modal para establecer la respuesta del sistema
debida a la perturbación aplicada. Se incluye la combinación del método de diferencias finitas y
el método de elementos finitos para resolver problemas estructurales dinámicos transitorios, o sea
variables en el tiempo.
Hacia el final del libro se incluye una serie de Apéndices, en los cuales se provee material a los
estudiantes que podrı́a ser nuevo para ellos, o que pueden tener la condición de material yá conocido
en forma previa y que de alguna manera constituyan conceptos teóricos que hayan sido yá olvidados
debido a su escasa utilización en el transcurso del tiempo hasta llegar a esta instancia de aprendizaje
del método de elemento finito. Esta serie se compone de los siguientes temas:
Apéndice A — es una revisión del álgebra matricial y debe ser material conocido por el estudiante
desde un curso de álgebra lineal, común en todos los programas de estudio de ingenierı́a.
Apéndice B — establece las relaciones constitutivas tridimensionales generales para un material
elástico, homogéneo, e isótropo. Yo he encontrado durante los años que imparto clases de mecánica
de sólidos, que los estudiantes de pre–grado no tienen un dominio firme de estas relaciones. En
general, el estudiante ha sido expuesto a tantos casos especiales que las ecuaciones tridimensionales
no son comprendidas de verdad.
Apéndice C — cubre tres métodos para resolver ecuaciones algebraicas lineales. Algunos estudian-
tes pueden usar este material como algoritmos para programar métodos de solución. Sólamente
incluı́mos éste Apéndice para que el lector sea consciente de los algoritmos que tienen presencia
debajo del software que usarán en el proceso de resolver computacionalmente los problemas de
elemento finito que les sean planteados.
Apéndice D — describe las capacidades computacionales básicas del software denominado FEPC.
El programa de elemento finito para computador personal — FEPC (Finite Element Personal
Computer) — fué desarrollado por el Dr. Charles Knight en el Instituto Politécnico de Virginia y
la Universidad Estatal; y se usa junto con este texto haciendo uso del permiso de registro de pro-
piedad. Los programas del Dr. Knight permiten el análisis de problemas bi–dimensionales usando
elementos: barra, viga y superficiales planos. El Apéndice describe en general las capacidades
y limitaciones del software. El programa FEPC está disponible para el estudiante mediante el
portal de sitio en Internet: www.mhhe.com/hutton.
Apéndice E — incluye los problemas (tipo mini–proyecto simulado) para varios capı́tulos del texto
que deben ser resueltos mediante software de elemento finito. La dimensión de los problemas
planteado no es tan grande, existiendo posibilidad de ser resueltos también mediante un pequeño
paquete de caracterı́stica académica (existen algunos que pueden descargarse desde la red Inter-
net). Adicionalmente, en el transcurrir del tiempo se agregarán problemas de esta clase a éste
Apéndice (en ediciones posteriores del libro) y también a algún sitio WEB a ser creado, en la
perspectiva de continuar la base de ejercicios propuestos que fueron puestos en este Apéndice.
El presente libro fué escrito para servir de material bibliográfico principal (libro de texto) en el
curso regular del método de elemento finito, que se dicta en el programa de Postgrado de Maestrı́a
en Ciencias de la Ingenierı́a Mecánica, de la Universidad Mayor de San Andrés (UMSA) – La Paz,
bolivia.
Agradezco a las autoridades de la carrera de Ingenierı́a Mecánica de la UMSA por depositar su
confianza y brindarme total respaldo para cumplir las funciones de docente de esta asignatura. Asimis-
mo, agradezco la orientación de mi amigo: PhD.Eng. J. Roger Saravia Luna (†), en la estructuración
del contenido temático y sus sabios consejos para abordar la redacción de éste documento.
Jaime G. Molina P.
Verano de 2010
fundamentos del método
de elemento finito
Capítulo
Conceptos básicos
1
El método de elemento finito es un método numérico general para la aproximación de soluciones de
ecuaciones diferenciales (ordinarias y parciales) muy utilizado en diversos problemas principalmente
de ingenierı́a y fı́sica–matemática aplicada.
Tı́picamente el método de los elementos finitos se programa computacionalmente para calcular el
campo de la variable dependiente del problema y, posteriormente, a través de relaciones generales y
constitutivas se calculan las variables secundarias de interés luego que se ha obtenido la solución para
el modelo de elemento finito que ha sido formulado. Esta técnica tuvo su origen en el tratamiento
de problemas planteados por la mecánica de sólidos deformables, pero en la actualidad es más gene-
ralmente aplicable a cualquier problema de la mecánica de medios continuos (que involucra a todos
los estados agregados de la materia). El método de los elementos finitos es muy usado debido a su
generalidad y a la facilidad de introducir dominios de cálculo complejos (en dos o tres dimensiones).
Además el método es fácilmente adaptable a problemas de transmisión de calor, de mecánica de fluidos
para calcular campos de velocidades y presiones (fluidodinámica) o de campo electromagnético. Dada
la imposibilidad práctica de encontrar la solución analı́tica de estos problemas, con frecuencia en la
práctica ingenieril los métodos numéricos y, en particular, los elementos finitos, se convierten en la
única alternativa práctica de cálculo.
Una importante propiedad del método es la convergencia; si se consideran particiones de elementos
finitos sucesivamente más finas, la solución numérica calculada converge rápidamente hacia la solución
exacta del sistema de ecuaciones.
1.1. Introducción
El método de elemento finito (MEF), a veces llamado análisis de elemento finito (AEF), es una
técnica computacional usada para obtener soluciones aproximadas de problemas con valores lı́mite en
diversos campos de la ingenierı́a. Declarado de manera simple, los problemas con valores lı́mite son
problemas matemáticos en los cuáles una o más variables dependientes deben satisfacer una ecuación
diferencial en todas partes dentro de un dominio conocido de variables independientes y satisfacer
1
2 CAPÍTULO 1. CONCEPTOS BÁSICOS
condiciones especı́ficas en el lı́mite del dominio de definición del problema. Los problemas de valores
lı́mite también se llaman a veces problemas de campo. El campo es el dominio de interés y a menudo
representa una estructura fı́sica. Las variables de campo son las variables dependientes de interés gober-
nadas por la ecuación diferencial. Las condiciones lı́mite son los valores especificados de las variables
de campo (o variables relacionadas como sus derivadas) sobre los lı́mites del campo. Dependiendo del
tipo de problema fı́sico que se analiza, las variables de campo pueden incluir desplazamiento fı́sico,
temperatura, flujo de calor, y velocidad de flujo, para mencionar sólamente algunas variables que se
presentan.
(x, y) (x, y)
y y
2
x x
P(x,y) 3
1
(a) Dominio de definición de la variable (b) Elemento finito triangular con nodos
de campo φ(x, y) vértice
(x, y)
Un pequeño elemento triangular que encierra un sub–dominio de tamaño–finito del área de interés
se muestra en la Figura 1.1(b). Que este elemento no sea un elemento diferencial de tamaño dx×dy
hace del mismo un elemento finito; es decir, un elemento de dimensiones finitas (nó infinitesimales !).
Cuando nosotros tratamos este ejemplo como un problema bi–dimensional, se supone que el espesor en
la dirección z es constante y la dependencia según esta dirección no se indica en la ecuación diferencial.
Los vértices del elemento triangular se numeran para indicar que estos puntos son los nodos del elemento
finito. Un nodo es un punto especı́fico perteneciente al elemento finito en el que el valor de la variable
de campo va a ser calculado explı́citamente. Los nodos exteriores son localizados en los lı́mites del
elemento finito y podrı́an usarse para conectar un elemento dado a elementos finitos adyacentes. Los
nodos que no se ubican en los lı́mites del elemento son nodos interiores y no pueden ser conectados a
los nodos de cualquier otro elemento. El elemento triangular de la Figura 1.1(b) tiene sólo nodos en
sus vértices, y por ello los mismos son nodos exteriores.
Si sólo se computan los valores de la variable de campo en los nodos del elemento, cómo se obtienen
los valores en otros puntos dentro de un elemento finito?. La respuesta a esta crucial pregunta contiene
la esencia del método de elemento finito: Los valores de la variable de campo calculados para los nodos
se usan para aproximar los valores en los puntos no–nodales (es decir, en todos los puntos al interior
del elemento) mediante un procedimiento de interpolación usando los valores nodales yá calculados en
forma previa. Para el ejemplo del elemento triangular con nodos en sus vértices, en el que los nodos
son todos exteriores, la variable de campo en cualquier otro punto dentro del elemento y en sus aristas
limitantes, se describe por la relación aproximada
donde φ1 , φ2 , y φ3 son los valores de la variable de campo en los puntos nodales, y N1 , N2 , y N3 son
las funciones de interpolación, también conocidas como funciones de forma o funciones de mezcla. En
el ámbito del método de elemento finito, los valores nodales de la variable de campo se tratan como
constantes desconocidas que serán determinadas. Las funciones de interpolación son más a menudo
formas polinómicas de las variables independientes, deducidas para satisfacer ciertas condiciones re-
queridas en los nodos. Estas condiciones son discutidas en detalle en los capı́tulos subsecuentes. El
punto relevante a ser enunciado aquı́ es que las funciones de interpolación son predeterminadas, con
formas de funciones conocidas de las variables independientes; y estas funciones describen la variación
de la variable de campo dentro del elemento finito.
Se dice que el elemento triangular descrito por la Ecuación (1.1) tiene 3 grados de libertad, en virtud
que se requieren tres valores nodales de la variable de campo para describir a la variable dependiente
en todos los puntos que conforman el elemento. Éste serı́a el caso si la variable de campo representa
una cantidad escalar, como la temperatura en un problema de trasmisión de calor , por ejemplo (véase
el Capı́tulo 7). Si el dominio de Figura 1.1 representa un cuerpo sólido delgado, sujeto a un estado
de tensión plana (véase el Capı́tulo 9); la variable de campo se adopta como el vector desplazamiento
interno y en éste caso deben computarse los valores de dos componentes en cada nodo. En este último
caso, el elemento triangular con nodos en sus tres vértices tiene 6 grados de libertad. En general, el
número de grados de libertad asociados con un elemento finito es igual al producto del número de
nodos y el número de valores del campo variable (y posiblemente sus derivadas) que deben calcularse
en cada nodo.
En esta instancia surge otra pregunta: Cómo este procedimiento basado en el elemento, puede
ser ahora aplicado sobre el dominio entero de interés ?. Como se bosqueja en la Figura 1.1(c), cada
elemento se conecta mediante sus nodos exteriores a los otros elementos circundantes a él. Las ecua-
ciones de elemento finito se formulan de modo que, en las conexiones nodales, el valor de la variable en
cualquier punto de conexión común es el mismo para cada elemento que concurre a dicha ubicación y
está conectado al nodo. Ası́, se asegura la continuidad de la variable de campo (variable dependiente
en el problema) en los puntos nodales. De hecho, las formulaciones de elemento finito son tales que
la continuidad del campo variable lo largo de los lı́mites de elementos aledaños que comparten dicho
lı́mite también se asegura. Este rasgo de imposición sobre el método en cuanto a la continuidad de la
4 CAPÍTULO 1. CONCEPTOS BÁSICOS
variable que está siendo modelada y aproximada mediante el método de elemento finito, evita la posi-
bilidad fı́sicamente inaceptable de existencia de huecos o discontinuidades en el dominio considerado.
En problemas estructurales, tales huecos fı́sicamente representarı́an la separación del material. En los
procesos de transmisión de calor, un “hueco” se manifestarı́a en sı́–mismo en la forma de posibilidad
de valores de magnitud de temperaturas diferentes correspondientes al mismo punto fı́sico, lo cual por
cierto es una completa aberración.
Aunque la continuidad del campo variable de elemento a elemento adyascente es inherente a la
formulación de elemento finito, la continuidad de los gradientes (i.e., derivadas) del campo de la variable
dependiente generalmente no existe. Ésta es una observación crı́tica hacia el método de elemento finito.
En la mayorı́a de los casos, dichas derivadas son de más interés que los valores del campo variable en
sı́–mismo. Por ejemplo, en los problemas estructurales la variable de campo es el desplazamiento,
pero el verdadero interés está más a menudo en la deformación y la tensión. Como la deformación
unitaria (al igual que la tensión interna) está definida en términos de las primeras derivadas espaciales
de las componentes de desplazamiento, por lo que se refiere a esta variable mecánica resultará que
no es contı́nua sobre la frontera lı́mite del elemento (exceptuando sus puntos nodales exteriores). Sin
embargo, las magnitudes de las discontinuidades de las derivadas de la función que se esté manipulando
pueden ser usadas para evaluar la exactitud de la solución y su convergencia hacia la solución real (o
exácta) a medida que el número de elementos se incrementa, como se ilustra a través del ejemplo que
a continuación mostraremos.
(a) Modelo discreto elaborado median- (b) Malla de elemento finito refinada
te uso de elementos cuadrados usando elementos más pequeños
Si las funciones de interpolación satisfacen ciertos requisitos matemáticos (véase el Capı́tulo 6),
una solución de elemento finito para un problema particular converge hacia la solución exacta del
problema a medida que el modelo utilizado sea asociado a una malla de elementos de menor tamaño.
Es decir, a manera que el número de elementos se incrementa y se disminuyen las dimensiones fı́sicas
1.2. CONCEPCIÓN DEL MÉTODO DE ELEMENTO FINITO ? 5
de los mismos, la solución de elemento finito cambia incrementalmente hacia la solución correcta. Los
cambios incrementales decrecen con el proceso de refinamiento de la malla y la aproximación hacia la
solución exacta se verifica que tiene forma asintótica.
Para ilustrar la convergencia del método, consideraremos un problema relativamente simple que
tiene una solución conocida. La Figura 1.3(a) muestra un cuerpo sólido de forma cilı́ndrica tronco–
cónica, empotrado en un extremo y sujeto a una carga de tracción en su otro extremo libre. Asumiendo
que el desplazamiento del punto de la carga aplicada sea de interés en su evaluación, se obtiene
una primera aproximación de solución considerando el cuerpo como cilı́ndrico de sección transversal
uniforme, teniendo ésta variable magnitud igual al valor medio del área transversal del cuerpo original
(véase la Figura 1.3b).
r0
rL A0+AL
A= 2
P
(a) Esquema del sistema fı́sico (b) Modelo de un elemento finito
(c) Modelo de dos elementos finitos (d) Modelo de cuatro elementos finitos
La varilla uniforme gruesa de la Figura 1.3(b) es un elemento finito barra (véase el Capı́tulo 2); ası́,
nuestra primera aproximación es un modelo de un solo elemento finito. La solución es obtenida usando
la teorı́a de la mecánica de materiales. Luego, modelamos el cilindro de sección transversal linealmente–
variable mediante dos barras uniformes en serie, como se muestra en la Figura 1.3(c). En el modelo de
dos elementos, cada uno de ellos es de longitud igual a la mitad de la longitud total del cilindro original
y tienen un área de sección transversal igual al área media de la correspondiente mitad de longitud del
cuerpo sólido original; las cuales puede obtenerse en base a la ecuación: r(x) = r0 − (x/L)(r0 − rL ), que
describe la variación del radio con respecto a la ubicación según dirección axial para el cuerpo tronco–
cónico original. El proceso de refinamiento de la malla se continúa, elaborando un modelo de cuatro
elementos como se ve en la Figura 1.3(d), y ası́ sucesivamente. En el modelo de cuatro elementos, la
longitud de los mismos es la cuarta parte de la longitud del cuerpo original y su sección transversal
6 CAPÍTULO 1. CONCEPTOS BÁSICOS
tiene área idéntica al área media del pedazo correspondiente del objeto cilı́ndrico tronco–cónico, como
se vé en la Figura recién mencionada. En todos los modelos presentados en la Figura 1.3, en trazo
punteado mostramos el objeto sólido cilı́ndrico tronco–cónico original.
Para este problema simple, el desplazamiento del extremo inferior del cilindro δ(x = L) para cada
uno de los modelos de elementos finitos es como se muestra en la Figura 1.4(a), dónde la lı́nea sólida
representa la solución exacta (conocida). La convergencia de las soluciones de elemento finito hacia la
solución exacta se indica claramente en la gráfica mencionada mediante una lı́nea curva segmentada.
Solución exacta
Solución aproximada
Solución exacta (cuatro elementos)
( x L )
(x/L)
Por otra parte, si trazamos el desplazamiento como una función de la posición a lo largo de la
longitud del cilindro, podemos observar la convergencia ası́ como también la naturaleza aproximada de
las soluciones de elementos finitos. La Figura 1.4(b) muestra la solución exacta y la solución del des-
plazamiento para el modelo de cuatro-elementos. En primer lugar, podemos notar que la variación del
desplazamiento en cada elemento es una aproximación lineal a la verdadera solución no–lineal. La va-
riación lineal en los sub–dominios es directamente atribuible al hecho que las funciones de interpolación
para un elemento–barra son lineales. Segundo, notamos que cuando la malla es refinada, la solución
del desplazamiento converge hacia la solución no–lineal en cada punto del dominio de definición del
problema.
El párrafo anterior discutió la convergencia del desplazamiento del extremo del cilindro tronco–cóni-
co. Como se verá en el Capı́tulo 2, el desplazamiento es la variable de campo primaria en los problemas
estructurales. En la mayorı́a de los problemas, sin embargo, estamos interesados principalmente en las
tensiones inducidas por las cargas especificadas. Las tensiones deben ser calculadas a través de las rela-
ciones tensión–deformación apropiadas (ley de Hocke generalizada), y las componentes de deformación
unitaria se obtienen de la solución del campo de desplazamientos. De aquı́, tensiones y deformaciones
se refieren, por guardar esta dependencia, como variables derivadas.
Por ejemplo, si trazamos la variación teórica de las tensiones internas para el ejemplo del cuerpo
ciı́ndrico tronco–cónico recién citado, obtenida desde la solución exácta, como también las soluciones
de elemento finito para los modelos de dos y cuatro elementos como se muestra en la Figura 1.5
(donde σ0 = P/A0 ), observamos que las tensiones son constantes en el interior de cada elemento y
representan una solución discontinua del problema por lo que se refiere a las tensiones y deformaciones.
También notamos que, a medida que aumenta el número de los elementos en el modelo, los saltos de
las discontinuidades en la tensión interna disminuyen en su magnitud. Este fenómeno es caracterı́stico
del método de elemento finito. La formulación del método de elemento finito para un problema dado
es tal que la variable de campo primaria es continua de elemento al elemento adyascente, pero las
variables derivadas (como la tensión interna, por ejemplo) no necesariamente poseen continuidad. En
1.2. CONCEPCIÓN DEL MÉTODO DE ELEMENTO FINITO ? 7
el proceso lı́mite de refinamiento de la malla, las variables derivadas se ponen cada vez más cercanas
al estado de continuidad si el número de elementos de la malla se incrementa.
4.5
4.0
Solución exacta
Dos elementos
3.5 Cuatro elementos
3.0
0
2.5
2.0
1.5
1.0
0 0.25 0.5 0.75 1.0
x/L
Nuestro ejemplo muestra cómo la solución de elemento finito converge a la solución exacta conocida
(la exactitud de la solución en este caso es aquella marcada por la teorı́a de la mecánica de los
materiales). Si nosotros conocemos la solución exacta, no estarı́amos aplicando el método de elemento
finito !. Ası́, se plantea la interrogante: Cómo aseguramos nosotros la exactitud de una solución obtenida
mediante el método de elemento finito para un problema con una solución teórica desconocida ?. La
respuesta a esta pregunta no es simple. Si nosotros no tuviéramos la lı́nea sólida en la Figura 1.5
que representa la solución exacta, todavı́a podrı́amos discernir la convergencia hacia la solución. La
convergencia de un método numérico (como el método de elemento finito) no está dada por ninguna
convicción de que los medios utilizados indiquen que la convergencia es precisamente hacia la solución
correcta (porque podrı́a suceder que la aproximación sea hacia un valor erróneo). Una persona que
usa la técnica de análisis de elemento finito, debe examinar la solución muy metódicamente por lo
que se refiere a: (1) la convergencia numérica, (2) la racionalidad (el resultado obtenido, tiene sentido
fı́sico ?), (3) si las leyes fı́sicas del problema están satisfechas (La estructura está en equilibrio ? Existe
balance entre el calor de entrada y el calor de salida ?), y (4) si los valores de las discontinuidades
de las variables derivadas en los lı́mites inter–elementos son razonables. Deben proponerse muchas de
tales preguntas y deben examinarse prioritariamente las respuestas antes de aceptar los resultados de
un análisis de elemento finito como representativos de una solución correcta y útil para los propósitos
del diseño.
donde debemos advertir que la igualdad debe ser tomada como “aproximadamente igual ”. De la teorı́a
de ecuaciones diferenciales, nosotros sabemos que la solución de una ecuación diferencial de primer–
orden contiene una constante de integración. La constante de integración debe ser determinada de
modo que una condición pre–establecida (una condición lı́mite o una condición inicial) sea satisfecha.
En el actual ejemplo, asumiremos que la condición especificada es x(0) = A (cte). Si escogemos un
paso de integración ∆x que sea un valor constante y pequeño (no se exige que el paso de integración
sea necesariamente constante), entonces podemos escribir
xi+1 = xi + ∆x i = 0, 1, . . . , N (1.5)
donde N es el número total de pasos requeridos para cubrir el dominio de definición de la variable
independiente. La Ecuación (1.4) entonces puede escribirse como
fi+1 = fi − xi ∆x f0 = A i = 0, 1, . . . , N (1.6)
La Ecuación (1.6) es conocida como una relación de recurrencia y provee una aproximación al valor de
la función desconocida f (x) en un número determinado de puntos discretos en el dominio del problema.
Para ilustrar este método alternativo, la Figura 1.6 muestra la solución exacta f (x) = 1 − x2 /2 , y
una solución de diferencias finitas obtenida con ∆x = 0,1 y considerando como valor de la constante
f0 = A = 1. La solución de diferencias finitas sólo es mostrada en los puntos discretos de evaluación
de la función. La manera de variación de la función entre los puntos calculados no es conocida en el
método de diferencias finitas; pero, mostramos una curva de ajuste en trazo segmentado que pasa a
través de los puntos calculados. Uno puede, sin embargo, interpolar linealmente los valores hallados para
producir una aproximación a la curva de la solución exacta; pero, la manera de efectuar este proceso de
interpolación no es algo relevante, ni de determinación “a priori ” en el método de diferencias finitas.
Para contrastar el método de diferencias finitas con el método de elemento finito, notamos que, en
el método de elemento finito, la variación de la variable de campo en el dominio fı́sico es una parte in-
tegral del procedimiento. Es decir, basados en las funciones de interpolación seleccionadas, la variación
de la variable de campo a lo largo de un elemento finito se especifica como una parte integral de la
formulación del problema. En el método de diferencias finitas, éste no es el caso: La variable de campo
se computa sólo en puntos especificados. La ramificación mayor de este contraste es que las derivadas
(hasta un cierto nivel) pueden evaluarse en el procedimiento de elemento finito, considerando que en
contraposición el método de diferencias finitas sólo proporciona datos en la propia variable (y nó en
cantidades que provengan de ella). En un problema estructural por ejemplo, ambos métodos propor-
cionan las soluciones del desplazamiento, pero la solución de elemento finito puede usarse para calcular
1.3. PROCEDIMIENTO GENERAL 9
1.0
0.8
f(x)
diferencia finitas
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1.0
x
las componentes de tensión directamente (en términos de la primera derivada). Para obtener datos
de tensión en el método de diferencias finitas se requieren consideraciones adicionales no inherentes al
modelo matemático mismo.
Hay también ciertas similitudes entre los dos métodos. Los puntos de integración en el método de
diferencias finitas son análogos a los nodos en un modelo de elemento finito. La variable de interés se
evalúa explı́citamente en dichos puntos. También, como el paso de integración (el tamaño del paso) en
el método de diferencias finitas es muy reducido si se desea un grado aceptable de exactitud, se espera
que la solución obtenida por éste método tenga convergencia hacia la solución exacta. Esto es similar a
la convergencia esperada de una solución de elemento finito del modo en el que la malla de elementos
es refinada. En ambos casos, el refinamiento representa la reducción del modelo matemático desde
términos finitos hasta infinitesimales. Y en ambos casos, las ecuaciones diferenciales son reducidas a
ecuaciones algebráicas.
Probablemente la manera más descriptiva de contrastar los dos métodos es notar que el método de
diferencias finitas modela la ecuacion diferencial del problema y usa la integración numérica para obte-
ner la solución en los puntos de discretización. El método de elemento finito modela el dominio entero
del problema y usa principios fı́sicos conocidos para desarrollar ecuaciones algebráicas que describen
las soluciones aproximadas. Ası́, el método de diferencias finitas modela las ecuaciones diferenciales
gobernantes del fenómeno en estudio, mientras que el método de elemento finito puede decirse que
modela el problema fı́sico más estrechamente desde el principio. Como se observará en el resto de este
texto, hay casos en que una combinación de los métodos de elemento finito y diferencias finitas es muy
útil y eficaz para obtener soluciones a los problemas que plantea la ingenierı́a, particularmente donde
los efectos dinámicos (transitorios) son importantes.
1.3.1. Pre–procesamiento
El paso de pre–procesamiento es, muy generalmente, descrito como el proceso de definir al modelo
e incluye las siguientes acciones a ejecutar:
Definir el dominio geométrico del problema.
Definir el tipo del elemento a ser utilizado (véase el Capı́tulo 6).
Definir las propiedades materiales de los elementos.
Definir las propiedades geométricas de los elementos (la longitud, el área, etc).
Definir las conectividades de los elementos (el modelo de la malla).
Definir las restricciones fı́sicas (condiciones de borde lı́mite).
Definir la solicitación o perturbación externa.
El paso de pre–procesamiento (la definición del modelo) es realmente crı́tico. En ningún caso es más
apropiado el ejemplo del axioma relacionado con los procesos relacionados con el uso de la computadora
“basura entrante, resulta en basura saliente”. Una solución de elemento finito perfectamente calculada
numéricamente no posee absolutamente ningún valor si corresponde a un pésimo proceso de modelado,
o a un problema erróneamente planteado.
1.3.3. Post–procesamiento
El análisis y la evaluación de los resultados de la solución es denominado como post–procesamiento.
El software asociado a este proceso importante del método de elemento finito contiene rutinas sofisti-
cadas usadas para ordenar, imprimir, y trazar los resultados seleccionados de una solución que haya
sido obtenida. Los ejemplos de operaciones que puede lograrse en esta fase incluyen
Clasificar las tensiones internas de elemento en orden de magnitud.
Verificar el equilibrio estático estructural.
Calcular factores de seguridad de diseño.
Bosquejar la deformación estructural producida.
Producir esquemas dinámicos de la respuesta del modelo.
Producir gráficas de color–codificado para el campo de temperaturas.
Puesto que los datos de la solución pueden manipularse de muchas maneras en la fase final de
post–procesamiento, el objetivo más importante es aplicar de modo legı́timo el juicio de ingenierı́a con
el objetivo de determinar las condiciones bajo las cuales los resultados de la solución obtenida son
fı́sicamente razonables.
1.4. BREVE HISTORIA DEL MÉTODO DE ELEMENTO FINITO 11
Figura 1.7: Malla de elementos finitos sobre una región rectangular que posee
una perforación circular central (elementos cuadriláteros).
La malla de elemento finito mostrada en la Figura 1.7 podrı́a representar el modelo de elemento
finito de variados problemas fı́sicos. Por ejemplo, para el análisis de tensiones internas, la geometrı́a
podrı́a representar una placa plana delgada con una perforación circular central sujeta a solicitación
aplicada en las aristas lı́mite mediante cargas contenidas en el propio plano de la placa. En este caso,
la solución de elemento finito se usarı́a para examinar los efectos de concentración de tensiones en la
vecindad del agujero. La malla de elementos mostrada también podrı́a representar el caso de flujo fluido
alrededor de un cilindro redondo sólido. Todavı́a podemos pensar otra aplicación, el modelo mostrado
podrı́a representar una aleta de transmisión de calor atravezada por una tuberı́a (el agujero) desde la
cual el calor se transfiere a la aleta para ser luego dispersado hacia el medio ambiente. En cada caso,
la formulación de las ecuaciones gobernantes del comportamiento fı́sico de los elementos en respuesta
a las influencias externas realmente es diferente para los problemas antes mencionados.
La Figura 1.8(a) muestra un módulo de una cercha tri–dimensional que alguna vez fué considerado
como un bloque–constructivo o módulo básico para la construcción de una estación espacial [52].
Diseñado para plegarse en forma de acordeón en un volumen pequeño para su transporte hacia la órbita
1.5. EJEMPLOS DE ANÁLISIS DE ELEMENTO FINITO 13
planetaria; el módulo cuando era desplegado, se extendı́a a las dimensiones globales de 1,4 m × 1,4 m
× 2,8 m. Conectando por sus extremos estos módulos, se podı́a armar una estructura esencialmente
de cualquier tamaño que se desee. La estructura se analizó mediante el método de elemento finito
para determinar las caracterı́sticas de vibración en función del número de módulos que participaban
para establecer la longitud global, cuando el mismo era variado. Como las conexiones entre los varios
miembros estructurales es del tipo alfiler o de juntas de rótula (articuladas), se usó en el modelo un
elemento axial simple solicitado en estado de tracción–compresión (véase el Capı́tulo 2). El modelo de
elemento finito de un módulo estaba compuesto de 33 elementos. Una muestra de la forma de vibración
(la deformación producida asociada al modo principal) de una cercha de cinco–módulos se muestra en
la Figura 1.8(b).
El ejemplo de la cercha plegable recientemente descrito involucra una estructura bastante grande
modelada por un número pequeño de elementos finitos relativamente largos. En contraste, la Figura 1.9
muestra que el modelo del elemento finito de un tubo muy delgado diseñado para uso en una aplicación
de un intercambiador de calor muy compácto. El tubo tiene diámetro interno de 0,976 pulg, espesor
de pared 0,00197 pulg, y longitud total de 36 pulg. El material considerado para la construcción del
tubo fué una aleación de cobre y titanio. Debido al espesor de pared, se encontraron que los tubos
del prototipo eran muy frágiles y difı́ciles de manejar sin causarles daño. Los objetivos del análisis
de elemento finito eran examinar la flexión, la torsión, y las cargas de pandeo permisibles. La figura
muestra la malla de elemento finito usada para modelar una sección del tubo de sólo 0,25 pulg de
longitud. Este modelo contiene 1920 elementos sólidos tridimensionales, y cada uno de ellos tiene ocho
nodos con 3 grados de libertad en cada nodo. Era requerido un número tan grande de elementos para
una estructura pequeña en la consideración de lograr exactitud computacional en los resultados. La
caracterı́stica involucrada aquı́ es la llamada relación de aspecto de los elementos, lo cual se define y
discute en los capı́tulos subsecuentes.
Como un último ejemplo, la Figura 1.10(a) representa al modelo de elemento finito del componente
de transporte de carga principal de un dispositivo de prótesis (para sustituir una mano). El dispositivo
tiene intención de acoplar esta mano a un brazo artificial. En el uso, la mano permitirı́a armar al
amputado de un brazo más corto para permitirle el levantamiento de pesas como parte de un programa
de salud. El modelo de elemento finito fue usado para determinar la distribución de tensiones en el
componente por lo que se refiere al rango de peso a cargarse por anticipado, para determinar un tamaño
apropiado de este dispositivo y también para seleccionar el material más óptimo. La Figura 1.10(b)
muestra un prototipo construı́do de la mano, una vez que el proceso de diseño ha sido completado.
14 CAPÍTULO 1. CONCEPTOS BÁSICOS
Los ejemplos mostrados en esta sección evidentemente no pueden ser resueltos manualmente debido
al elevado número de elementos contenidos en la malla que define el modelo de análisis; por ello para
obtener la solución de los problemas que plantean deberá hacerse uso de métodos computacionales
implementados en paquetes informáticos de cálculo numérico.
requiere para el uso exitoso de este texto. El estudio del requisito previo está basado en las asignaturas
más comunes de cualquier programa de estudios de ingenierı́a: el álgebra lineal, el cálculo a través de
las ecuaciones diferenciales, y la serie usual de estática, dinámica y mecánica de materiales. Aunque
no es requerido, el estudio previo de mecánica de fluidos y transmisión del calor serı́a de mucha ayuda.
Dado como conocido todo este cimiento temático académico, el método de elemento finito se desarrolla
en base a las leyes fı́sicas más utilizadas (el equilibrio, la conservación de masa, la conservación de la
energı́a, y otras). Se dá cierto énfasis al principio de energı́a potencial mı́nima (ver el Capı́tulo 2), y el
método de elemento finito de Galerkin (introducido y desarrollado en el Capı́tulo 5).
A medida que usted como lector progrese a través del texto, descubrirá que cubrimos una cantidad
significativa de la teorı́a de elemento finito además de la aplicación mediante ejemplos con graduación
de dificultad. Dada la disponibilidad de muchos poderosos y sofisticados paquetes de software de
elemento finito, por qué la necesidad de estudiar la teorı́a ?. El método de elemento finito es una
herramienta, y como cualquier otra herramienta, usarla sin la instrucción apropiada realmente puede
ser hasta peligroso. Mi premisa es que una instrucción apropiada en este contexto incluye comprender
la teorı́a básica de la formulación subyacente en los modelos de elementos finitos de problemas fı́sicos.
Como fué establecido previamente, el análisis crı́tico de los resultados de un modelo computacional
de elemento finito es esencial, desde que esos resultados pueden en el futuro convertirse en la base
para el diseño. El conocimiento de la teorı́a es necesaria para un apropiado proceso de modelado y la
evaluación de los resultados computacionales.
Análisis uni–dimensional
Capítulo 2
En éste capı́tulo primero estableceremos las propiedades de rigidez de los elementos más simples que
pueden elaborarse asociados con la metodologı́a que pretendemos desarrollar a lo largo de los capı́tulos
subsecuentes; dichos elementos son aquellos que se denominan uni–dimensionales, debido a que en
su descripción de comportamiento requieren solamente de una dirección espacial. Por simplicidad
en la descripción de las relaciones carga–desplazamiento, utilizaremos primero un sistema referencial
coordenado definido con origen en algún punto del elemento, marco referencial al que lo consideraremos
fijo y solidario respecto a este cuerpo sólido deformable.
También definiremos un segundo sistema coordenado referencial, respecto del cual se efectúa la
descripción de comportamiento de toda la estructura (y por supuesto, de todos los elementos que la
componen). Éste sistema coordenado lo escogeremos con origen en algún punto arbitrario del espacio
con orientación espacial de ejes coincidente a áquel que posee el sistema coordenado propio del elemento
tı́pico escogido para ser analizado. La coincidencia de orientación espacial de ambos sistemas de ejes
es una exigencia en el capı́tulo presente, pero no es una restricción de la aplicación general del método
de elemento finito.
En base a las ecuaciones que dan la descripción de las propiedades de rigidez del elemento tı́pico
respecto de su propio sistema coordenado, efectuaremos un proceso de transformación de éste sistema
hacia el marco global único y común estructural (aquı́ algunas veces se requiere efectuar una trans-
formación de rotación y traslación del sistema de ejes coordenados). Al efectuar este procedimiento,
habremos logrado establecer la descripción de las propiedades de rigidez del elemento tı́pico respecto
ahora del marco referencial de toda la estructura; de modo que yá sea posible aplicar el principio de
superposición para lograr establecer las relaciones carga–desplazamiento de todo el sistema o estruc-
tura. Esto último, tiene la analogı́a de armar un rompe–cabezas con las piezas individuales dispuestas
sobre un hipotético tablero.
En esencia, y de modo literal muy resumido, éste será el procedimiento general que aplicaremos
con el objetivo de establecer las propiedades de rigidez de los diversos elementos tı́picos con los cuales
se conforman los diversos sistemas que son objeto de análisis en éste documento.
17
18 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
2.1. Introducción
Las caracterı́sticas primarias de un elemento finito son incluidas en la matriz de rigidez de elemento.
Para un elemento finito estructural, la matriz de rigidez contiene información de su geometrı́a y su
conducta material, que indica la resistencia del elemento a la deformación cuando sobre el mismo
se aplica un estado de carga. La deformación producida puede incluir efectos axiales, flexionantes,
torsionales y cizallantes. Para elementos finitos usados en análisis no–estructural, como flujo fluido y
transmisión de calor, también se usa el término matriz de rigidez, ya que este arreglo de coeficientes
representa la resistencia del elemento al cambio de su estado cuando es expuesto a las influencias
externas.
Este capı́tulo desarrolla las caracterı́sticas de elemento finito de dos elementos estructurales uni–
dimensionales, relativamente simples: un resorte lineal elástico y un miembro en estado de tensión–
compresión. Se seleccionan estas entidades como elementos introductorios porque la conducta de cada
uno de ellos es relativamente muy conocida a partir de materia teórica que normalmente estudiamos en
ingenierı́a, como son los temas de estática y mecánica de materiales. Ası́, el “puente” hacia el método
de elemento finito no es obscurecido por teorı́as nuevas para el estudiante de ingenierı́a. Más bien,
construimos sobre los principios conocidos de la ingenierı́a, para introducir los fundamentos del método
de elemento finito. El resorte lineal y el miembro de tensión–compresión (de ahora en adelante referido
como un elemento–barra, también conocido en la literatura de elemento finito como: listón, baqueta,
o elemento de cercha) también se usan para introducir el concepto de funciones de interpolación.
Como mencionamos brevemente en el Capı́tulo 1, la premisa básica del método de elemento finito
es describir la variación contı́nua de la variable de campo (en este capı́tulo, el desplazamiento fı́sico)
en términos de los valores discretos adoptados en los puntos nodales del elemento para esta variable
de interés. En el interior de un elemento finito, ası́ como a lo largo de sus lı́mites (aplicable en los
problemas bi y tri– tridimensionales), la variable de campo se describe por funciones de interpolación
(véase el Capı́tulo 6) que deben satisfacer condiciones prescritas.
El análisis de elemento finito está basado, dependiendo del tipo de problema, en varios principios
fı́sico/matemáticos. En la presente introducción al método, mostramos varios de tales principios apli-
cables al análisis de elemento finito. Primero, y principalmente, para los sistemas de resortes y barras
interconectadas utilizamos el principio de equilibrio estático, pero — y ésto es esencial — incluı́mos la
deformación en el desarrollo del análisis; es decir, no estamos tratando con aspectos relacionados a la
mecánica de cuerpo rı́gido.
Para la extensión del método de elemento finito a los sistemas estructurales elásticos más complica-
dos, también establecemos y aplicamos el primer teorema de Castigliano [43] y el principio ampliamente
usado de la energı́a potencial mı́nima [16]. El primer teorema de Castigliano, en la forma aquı́ pre-
sentada, puede parecer ser nuevo para el lector. El primer teorema es la contraparte y complemento
del segundo teorema, el cual se encuentra más a menudo en el estudio elemental de la mecánica de
materiales [29]. Ambos teoremas relacionan desplazamientos y las fuerzas aplicadas, a las condiciones
de equilibrio de un sistema mecánico; en términos de la energı́a acumulada al interior del material por
el proceso de deformación. El uso aquı́ del primer teorema de Castigliano es para el propósito distinto
de introducir el concepto de energı́a potencial mı́nima sin tener que recurrir a principios matemáticos
del cálculo variacional, los cuales se ubican más allá del nivel matemático pensado para este texto.
La formulación del resorte lineal como un elemento finito es logrado con referencia a la Figura 2.1(a).
Como un resorte elástico soporta sólo carga axial, seleccionamos un sistema de coordenadas de elemento
(también conocido como sistema coordenado local ) como un eje x orientado a lo largo de la longitud del
resorte con origen en el extremo izquierdo, como se muestra. El sistema coordenado local es solidario
al elemento (está adherido a él) y es escogido ası́, por conveniencia geométrica y por simplicidad en la
descripción de la conducta del elemento. El elemento o el sistema de coordenadas local se contrasta con
el sistema de coordenadas global. El sistema de coordenadas global es aquel sistema en el cual se describe
la conducta de una estructura completa (de ahı́ su denominación de global). Por estructura completa se
significa al ensamblaje de muchos elementos finitos (en este punto, varios resortes), sistema para el cual
deseamos evaluar la respuesta a las condiciones de carga aplicada. En este capı́tulo, tratamos con casos
en los que los sistemas de coordenadas locales y globales son esencialmente el mismo, excepto para una
traslación del origen que podrı́a considerarse si fuese necesario. En casos bi y tri–dimensionales, sin
embargo, las distinciones de los sistemas coordenados es bastante diferente y requiere transformación
matemática de los sistemas coordenados de elemento hacia una base común. La base común es el
sistema de coordenadas global.
Fuerza, f
u1 u2 k
1
f1 k f2
2 x Deflexión, u2 u1
1
(a) Fuerzas y desplazamientos nodales de elemento (b) Curva carga – deformación
f1 = −k(u2 − u1 ) (2.3a)
f2 = k(u2 − u1 ) (2.3b)
notando que como alternativa escogimos tomar f = −f1 = f2 , como par de fuerzas que solicitan en
tracción al resorte. Estas dos últimas ecuaciones pueden ser expresadas en forma matricial (véase el
Apéndice A para una revisión de álgebra lineal) como
k −k u1 f1
= (2.4)
−k k u2 f2
donde
k −k 1 −1
[ke ] = =k (2.6)
−k k −1 1
es definida como la matriz de rigidez del elemento resorte en el sistema local coordenado (o sistema
propio del elemento), {u} es la matriz columna (vector) de desplazamientos nodales, y {f } es la matriz
columna (vector) de cargas nodales de elemento (en capı́tulos subsecuentes la notación matricial es
usada extensivamente. Una matriz general es designada entre paréntesis rectos [ ] y una matriz columna
(o vector) entre paréntesis rizados o llaves { }).
La Ecuación (2.6) muestra que la matriz de rigidez de elemento para el resorte elástico lineal
es una matriz de dimensión 2×2 (ésto es, 2 filas y 2 columnas). Esto corresponde al hecho que el
elemento exhibe dos desplazamientos nodales (o grados de libertad) y que los dos desplazamientos no
son independientes (es decir, el cuerpo es contı́nuo y elástico). Además, la matriz es simétrica. Una
matriz simétrica tiene términos fuera de la diagonal principal tales que kij = kji . La simetrı́a de la
matriz de rigidez es indicativa del hecho que el cuerpo es linealmente elástico y cada desplazamiento se
relaciona al otro por el mismo fenómeno fı́sico. Por ejemplo, si una fuerza F (positiva, o de tracción)
es aplicada al nodo 2 con el nodo 1 sostenido fijamente, el desplazamiento relativo de los dos nodos
es igual como si la fuerza fuese simétricamente aplicada (negativa, o de tracción) al nodo 1 con el
nodo 2 inmóvil o fijo. (Se ven contra–ejemplos a la simetrı́a en análisis de transferencia de calor y flujo
fluido en los Capı́tulos 7 y 8). Como se verá a medida que se desarrollen elementos estructurales más
complicados, éste es un resultado general: Un elemento que exhibe N grados de libertad, tiene una
matriz de rigidez correspondiente cuya dimensión es N ×N .
Luego, consideremos la solución del sistema de ecuaciones representada sintéticamente por la Ecua-
ción (2.5). En general, las fuerzas nodales se prescriben (son datos) y el objetivo es resolver para los
desplazamientos nodales desconocidos. Formalmente, la solución se representa por
donde [ke ]−1 es la inversa de la matriz de rigidez de elemento. Sin embargo, esta matriz inversa no existe,
desde que el determinante de la matriz de rigidez de elemento es idénticamente cero. Por consiguiente,
la matriz de rigidez de elemento es singular, y esto también demuestra ser un resultado general en
la mayorı́a de los casos. La importancia fı́sica de la naturaleza singular de la matriz de rigidez de
elemento es encontrada por re–examinación de la Figura 2.1(a), que muestra que ninguna restricción
del desplazamiento en cualquier lugar se ha impuesto en el movimiento del elemento resorte; es decir,
el resorte no está conectado a cualquier objeto fı́sico o apoyo que prevendrı́a o limitarı́a el movimiento
de cualquiera de sus nodos. Sin restricciones al movimiento, no es posible resolver individualmente
2.2. EL RESORTE LINEALMENTE ELÁSTICO 21
para los desplazamientos nodales. En cambio, sólo la diferencia en los desplazamientos nodales puede
determinarse, como esta diferencia representa el alargamiento o contracción del elemento–resorte debido
a efectos elásticos de su comportamiento.
Como será discutido en más detalle en la formulación general de las funciones de interpolación
(Capı́tulo 6) y la dinámica estructural (Capı́tulo 10), una formulación apropiada de elemento finito
debe permitir un valor constante de la variable de campo. En el ejemplo que aquı́ es de nuestro interés,
esto significa movimiento de cuerpo rı́gido para el resorte; es decir que se mueve estando deformado
simultáneamente. Nosotros podemos identificar la capacidad de movimiento de cuerpo rı́gido por lo que
se refiere a un solo resorte (el elemento) y en el contexto de varios elementos conectados (un sistema o
estructura completa) que tiene capacidad de movimiento común de conjunto. Para un solo, elemento no
restringido, si las fuerzas arbitrarias son aplicadas a cada nodo, el resorte no sólo se deforma axialmente,
además también surge una aceleración según la segunda ley de Newton. De esto concluimos que no
sólo existe deformación sino también movimiento total del elemento (el mismo se mueve como cuerpo
rı́gido, pero en configuración estable deformada). Si, en un sistema conectado de elementos resorte, la
respuesta del sistema completo es tal que los nodos 1 y 2 de un elemento particular se desplazan la
misma cantidad, no hay ninguna deformación elástica del resorte y por consiguiente ninguna fuerza
elástica en el mismo. Esta situación fı́sica debe ser incluida en la formulación del elemento. La capacidad
de movimiento de cuerpo rı́gido se indica matemáticamente por la singularidad de la matriz de rigidez
de elemento. Cuando la matriz de rigidez se formula en base a la deformación del elemento, no podemos
esperar calcular los desplazamientos nodales si no hay ninguna deformación del mismo.
La Ecuación (2.7) indica la operación matemática de invertir la matriz de rigidez para obtener
la solución. En el contexto de un elemento individual, la naturaleza singular de una matriz de rigi-
dez de elemento evita esta operación, ya que la inversa de una matriz singular no existe. Como se
ilustra profusamente en el resto del texto, la solución general de un problema de elemento finito, en
una descripción global, como algo opuesto al contexto de un solo elemento, involucra la solución de
ecuaciones del aspecto matemático mostrado en la Ecuación (2.5). Para modelos reales de elemento
finito, los cuales son de grandes dimensiones en términos de la matriz de rigidez de toda la estructura
involucrada, cuya dimensión sea N ×N, N >> 10 ; el calcular su inversa es un procedimiento verdade-
ramente muy ineficaz, y de pérdida de tiempo si se pretende usar una simple calculadora electrónica.
Para estas situaciones, otras técnicas de solución más eficaces están disponibles, y éstas se discuten
más adelante. (Muchos de los problemas propuestos de fin–de–capı́tulo incluidos en este texto son de
orden pequeño y pueden resolverse eficazmente mediante proceso de inversión matricial usando un
software con hoja de cálculo, o paquetes destinados especı́ficamente al cómputo numérico matricial,
como Mathlab, Mathematica, o MathCad).
U1 U2 U3
1 2
F1 F2 F3
1 2 3 x
k1 k2
desplazamientos nodales globales se identifican como U1 , U2 , y U3 ; dónde la letra mayúscula se usa para
indicar que las cantidades representadas son globales o desplazamientos del sistema (no confundirlos
con los desplazamientos de elemento individuales). Similarmente, las fuerzas nodales aplicadas son F1 ,
F2 , y F3 .
u(1)
1 u(1)
2 u(2)
1 u(2)
2
1 2
f1(1) f2(1) f2(2) f3(2)
1 2 2 3
k1 k2
(a) Elemento 1 (b) Elemento 2
Asumiendo que el sistema de dos elementos–resorte está en perfecto equilibrio estático, examinamos
los diagramas de cuerpo–libre de los resortes individualmente ahora referidos a sus propios sistemas
coordenados (sistemas coordenados locales), que se aprecian en las Figuras 2.3(a) y 2.3(b); las cuales
expresan las condiciones de equilibrio para cada resorte. Ahora, usando la Ecuación (2.4) para cada
uno de los esquemas, tenemos
" #( ) ( )
k1 −k1 u(1)
1 f1(1)
= (2.8a)
−k1 k1 u(1)
2 f2(1)
" #( ) ( )
k2 −k2 u(2)
1 f2(2)
= (2.8b)
−k2 k2 u(2)
2 f3(2)
u(1)
1 = U1 u(1)
2 = U2 u(2)
1 = U2 u(2)
2 = U3 (2.9)
Las condiciones de compatibilidad establecen el hecho fı́sico que los resortes están conectados en el
nodo 2, y permanecen conectados a éste nodo después de la deformación; y por tanto deben tener
el mismo desplazamiento nodal en éste punto de conexión. Ası́, la continuidad del desplazamiento
2.2. EL RESORTE LINEALMENTE ELÁSTICO 23
elemento–a–elemento se refuerza en las conexiones nodales. Sustituyendo las Ecuaciones (2.9) en las
Ecuaciones (2.8), obtenemos
" #( ) ( )
k1 −k1 U1 f1(1)
= (2.10a)
−k1 k1 U2 f2(1)
" #( ) ( )
k2 −k2 U2 f2(2)
= (2.10b)
−k2 k2 U3 f3(2)
Aquı́, usamos la notación fi(j) para representar la fuerza ejercida sobre el elemento j en el nodo i.
Las Ecuaciones (2.10) son las formulaciones de equilibrio estático para cada elemento–resorte ex-
presadas en términos de los desplazamientos globales especificados. En esta forma, las ecuaciones
claramente muestran que los elementos se conectan fı́sicamente en el nodo 2, y tienen el mismo despla-
zamiento U2 en ese nodo. Estas ecuaciones no son todavı́a adecuadas para una combinación dirécta,
porque los vectores desplazamiento no son los mismos. Nosotros podemos expandir ambas ecuaciones
matriciales hacia una dimensión 3×3 como sigue (expresando formalmente el hecho que el elemento 1
no se conecta al nodo 3 y el elemento 2 no está conectado al nodo 1):
k1 −k1 0 U1 f1(1)
−k1 k1 0 U2 = f2(1) (2.11a)
0 0 0 0 0
0 0 0 0 0
0 k2 −k2 U2 = f2(2) (2.11b)
0 −k2 k2
(2)
U3 f3
Aquı́, debemos notar que en la Ecuación (2.11a) podrı́amos poner U3 , y U1 en la Ecuación (2.11b) en
los vectores desplazamiento, en lugar de los valores nulos explı́citos que pusimos (verifique usted que
las ecuaciones desarrolladas son idénticas con este cambio propuesto).
La superposición de las ecuaciones anteriores, como se puede fácilmente comprobar, dá como re-
sultado
k1 −k1 0 U1 f1(1)
−k1 k1 + k2 −k2 U2 = f2(1) + f2(2) (2.12)
0 −k2 k2 U3 f3(2)
Luego, refiriéndonos a los diagramas de cuerpo–libre de cada uno de los nodos, diagramados en las
Figuras 2.3(c), 2.3(d), y 2.3(e); nos muestran que las condiciones de equilibrio estático imponen el
cumplimiento de las relaciones
f1(1) = F1 f2(1) + f2(2) = F2 f3(2) = F3 (2.13)
respectivamente. Sustituyendo en la Ecuación (2.12), obtenemos el resultado final
k1 −k1 0 U1 F1
−k1 k1 + k2 −k2 U2 = F2 (2.14)
0 −k2 k2 U3 F3
Esta ecuación podrı́amos escribirla en forma sintética: [K ]{U } = {F }, cuya forma es similar al de la
Ecuación (2.5). Sin embargo, la Ecuación (2.14) representa a las relaciones que gobiernan al sistema
compuesto de dos elementos resorte inter–conectados. Por consideración dirécta de las condiciones de
equilibrio estático, hemos obtenido la matriz de rigidez del sistema [K] (note el uso del sı́mbolo en
mayúscula), llamada también matriz de rigidez global, como
k1 −k1 0
[K ] = −k1 k1 + k2 −k2 (2.15)
0 −k2 k2
24 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
Notemos que la matriz de rigidez del sistema implı́citamente contiene las siguientes propiedades:
(1) es simétrica, como es el caso con todos los sistemas lineales referidos a sistemas coordenados
ortogonales (cuyos ejes están en cuadratura unos respecto a los otros); (2) es singular, desde que no se
aplican restricciones para prevenir el movimiento de cuerpo rı́gido de todo el sistema; y (3) la matriz
del sistema simplemente es una superposición de las matrices de rigidez de los elementos individuales
con la asignación apropiada de los desplazamientos nodales de elemento (desplazamientos locales) y
asociados con los coeficientes de rigidez a los desplazamientos nodales del sistema (desplazamientos
globales). El procedimiento del superposición se formaliza en el contexto de estructuras del tipo marco
en los párrafos siguientes.
Ejemplo 2.1.
Considere el sistema de dos elementos–resorte mostrado en la Figura 2.2, dado que el Nodo 1 se sujeta
a un apoyo fijo proporcionando la restricción del desplazamiento U1 = 0. Considerando los valores:
k1 = 50 lb/in, k2 = 75 lb/in, F1 = F2 = 75 lb; determinar los desplazamientos U2 y U3 .
> Solución
Sustituyendo los valores especificados en la Ecuación (2.14), tendremos
50 −50 0 0 F1
−50 125 −75 U2 = 75
0 −75 75 U3 75
y notamos que, debido al valor de desplazamiento nulo en el nodo 1, la fuerza nodal F1 se convierte en
una fuerza de reacción desconocida (ésta serı́a la reacción de apoyo). Formalmente, la primera ecuación
algebráica representada en esta ecuación matricial llega a ser
−50 U2 = F1
y esta es conocida como una ecuación de restricción, ya que representa la condición de equilibrio de
un nodo en el cual el desplazamiento está restringido (con un valor pre–establecido). Desechando la
primera ecuación del sistema anterior, la segunda y tercera ecuaciones se escriben ahora como
125 −75 U2 75
=
−75 75 U3 75
Este último cálculo resulta ser una verificación de la solución obtenida para los desplazamientos. >
2.2. EL RESORTE LINEALMENTE ELÁSTICO 25
El Ejemplo 2.1 ilustra el procedimiento general para la solución de modelos de elemento finito:
Formular las ecuaciones de equilibrio del sistema, aplicar las condiciones de restricción especificadas,
resolver la serie reducida de ecuaciones para los desplazamientos “activos”, y sustituir los desplazamien-
tos evaluados en las ecuaciones de restricción para obtener las reacciones desconocidas. Aunque no es
directamente aplicable para el elemento–resorte, en las formulaciones de elemento finito más generales,
los desplazamientos calculados son también sustituidos en las relaciones desplazamiento – deforma-
ción para obtener las deformaciones unitarias; y éstas son, a su vez, sustituidas en las ecuaciones de
tensión–deformación unitaria aplicables, para obtener los valores de tensión interna de elemento.
Ejemplo 2.2.
La Figura 2.4(a) muestra un sistema de tres resortes linealmente elásticos que soportan tres cuerpos
con pesos de igual magnitud, suspendidos en un plano vertical. Tratando los resortes como elementos
finitos, determine el desplazamiento vertical de cada cuerpo.
3k U1
3k 1
W 2
U2
2k 2k 2
3
W
U3
k 3
k
4
W U4
(a) (b)
> Solución
Para tratar esto como un problema de elemento finito, nosotros asignamos cierta codificación numérica
a los nodos y elementos como se muestra en la Figura 2.4(b) e ignoramos, por el momento, que el
desplazamiento U1 se conoce que tiene valor nulo por la restricción de movimiento impuesto por el
apoyo fijo. Por la Ecuación (2.6), la matriz de rigidez de cada elemento es (pre–procesamiento)
(1) 3k −3k
[k ] =
−3k 3k
2k −2k
[k (2) ] =
−2k 2k
(3) k −k
[k ] =
−k k
26 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
u(1)
1 = U1 u(1)
2 = u(2)
1 = U2 u(2)
2 = u(3)
1 = U3 u(3)
2 = U4
Procediendo como en el ejemplo anterior, escribimos las ecuaciones individuales de elemento como
(1)
3k −3k 0 0 U1
f1
−3k 3k 0 0 U2 f2(1)
0
=
0 0 0 U3 0
0 0 0 0 U4 0
0 0 0 0 U1 0
(2)
0 2 −2k 0 U2 = f1(2)
0 −2k 2k 0 U3 f2
0 0 0 0 U4 0
0 0 0 0 U1
0
0 0 0 0 U2 = 0
0 0 k −k U3 f1(3)
(3)
0 0 −k k U4 f2
Sumando estas ecuaciones, factorizando el parámetro k que es común a todas las matrices,
3 −3 0 0 U1 F1
−3 5 −2 0 U2 W
k
0 −2 3 −1 U3 W = (1)
0 0 −1 1 U4 W
donde hemos utilizado el hecho que la suma de las fuerzas de elemento en cada nodo debe igualar a la
fuerza aplicada a ese nodo y; en el nodo 1, la fuerza actuante es una reacción desconocida.
Aplicando la restricción de desplazamiento U1 = 0 (esto es todavı́a pre–procesamiento), obtenemos
−3 k U2 = F1 (2)
para los desplazamientos activos. Nuevamente, notemos que la Ecuación (3) se obtiene eliminando la
Ecuación de restricción (2) de la Ecuación (1) de comportamiento global estructural.
La solución simultánea (procesamiento) del sistema de ecuaciones algebráicas representadas por la
Ecuación (3) dá los desplazamientos incógnitas con los valores indicados a continuación
W 2W 3W
U2 = U3 = U4 =
k k k
y la Ecuación (2) proporciona como fuerza de reacción de apoyo (ésto es post–procesamiento)
F1 = −3 W
Debemos notar que la fuerza de reacción de apoyo equilibra la fuerza externa aplicada al sistema
(o peso total de los cuerpos que cuelgan). >
2.2. EL RESORTE LINEALMENTE ELÁSTICO 27
Notemos que la solución aquı́ hallada aplicando el método de elemento finito para el ejemplo ante-
rior, es exáctamente aquella que hubiésemos obtenido mediante aplicación de las conocidas ecuaciones
de equilibrio estático. También notemos que hemos seguido el procedimiento general de utilización del
método, consistente de los pasos secuenciales siguientes:
Formular las matrices de rigidez individuales de elemento.
Escribir las relaciones de transformación de desplazamientos locales hacia globales.
Ensamblar la ecuación de equilibrio global en forma matricial.
Reducir la ecuación matricial según las restricciones especificadas.
Resolver el sistema de ecuaciones para los desplazamientos nodales desconocidos (variables pri-
marias).
Resolver las ecuaciones de restricción para las fuerzas de reacción (variables secundarias) por
substitución–regresiva.
Ejemplo 2.3.
La Figura 2.5 muestra un sistema de tres elementos resorte lineales conectados en configuración en
serie. La numeración de nodos y elementos adoptada es como se indica. El nodo 1 se fija rı́gidamente
para prevenir el movimiento, y en el nodo 3 se dá un desplazamiento especificado δ como es mostrado.
Las fuerzas F2 = −F y F4 = 2F están aplicadas en los nodos 2 y 4. Determine el desplazamiento de
cada nodo y la fuerza requerida en el nodo 3 para las condiciones especificadas previas.
F2 = – F
1 3
3 4
1 k 2 3k 2k F4 = 2 F
2
> Solución
Este ejemplo incluye una condición de borde no–homogénea. En los ejemplos anteriores, las condicio-
nes lı́mite se representaron por desplazamientos de valor nulo. En este ejemplo, tenemos ambos tipos:
condiciones de desplazamiento nulo (homogéneo) y de magnitud especificada (no–homogéneo). El tra-
tamiento algebráico debe ser diferente como sigue. Las ecuaciones de equilibrio del sistema se expresan
en forma matricial (véase el Problema 2.6) como
k −k 0 0 U1
F1
F1
−k −3k
4k 0 U2 = F2 = −F
0 −3k 5k −2k
U3
F3 F3
0 0 −2k 2k U4 F4 2F
juntamente con −k U2 = F1 , como ecuación de restricción. Pero, la serie de ecuaciones obtenida cla-
ramente muestra que simplemente no podemos eliminar la fila y columna que corresponden al des-
plazamiento no–nulo especificado, porque éste parámetro aparece en las ecuaciones que gobiernan los
desplazamientos activos. Para ilustrar un procedimiento general, volvemos a escribir la última ecuación
reordenándola como
5k −3k −2k δ F3
−3k 4k 0 U2 = −F
−2k 0 2k U4 2F
la cual puede ser resuelta para los desplazamientos globales incógnita, dando como ecuación de solución
con tal que la matriz [KU U ]−1 exista. Desde que las condiciones de restricción han sido aplicadas
correctamente, la matriz inversa para nuestro ejemplo sı́ existe y está dada por
1
0
[KU U ]−1 = 4k
1
0
2k
Substituyendo en la Ecuación (5), la solución para los desplazamientos inicialmente desconocidos re-
sulta F
1
3δ
0 −F + 3kδ
− +
U2 4k 4
4k
{U } = = 1 = (6)
U4 0
2F + 2kδ F +δ
2k k
2.2. EL RESORTE LINEALMENTE ELÁSTICO 29
la cual, al ser substituida por la solución para los desplazamientos ahora yá conocidos, proporciona la
solución para las fuerzas asociadas con los desplazamientos no–homogéneos que han sido especificados.
Reemplazando la solución previamente obtenida, especificada por la Ecuación (6), tendremos
F 3δ
− +
3
4k 4 5
F3 = 5 k δ − 3k 2k = kδ − F (8)
F +δ
4 4
k
Elemento 1
( 3 F
) ( )
4k δ − 4 f1(1)
k −k 0 −k U2
= = =
−k k U2 k U2 −3k δ + F f2(1)
4 4
lo cual muestra que las fuerzas nodales sobre el elemento 1 son iguales y de sentido opuesto como es
requerido para el equilibrio estático.
Elemento 2
( 3F
− 4k − 34 k δ
( F
− 4k + 34 δ
) ) ( )
f2(2)
3k −3k U2 3k −3k
= = =
−3k 3k U3 −3k 3k δ 3F
+ 43 k δ f3(2)
4k
Elemento 3
( (3) )
2k −2k U3 2k −2k δ −2F f3
= F = =
−2k 2k U4 −2k 2k k +δ 2F f4(3)
u1 u2
1 2 x
x u(x)
L
La Figura 2.6 muestra una barra elástica de longitud L, la cual no está apoyada y es asociada a un
sistema coordenado uniaxial denotado como x con su origen coincidente con el extremo izquierdo del
elemento (el nodo 1). Este es el sistema coordenado de elemento o marco de referencia local.
Denotando el desplazamiento axial a lo largo de la longitud de la barra en cualquier posición como
u(x), definimos los nodos 1 y 2 a cada extremo de la barra como se muestra, e introducimos los
desplazamientos nodales u1 = u(x = 0) y u2 = u(x = L). Ası́, tenemos la variable contı́nua de campo
u(x), la cual será expresada (aproximadamente) en términos de las variables de desplazmiento nodales
u1 y u2 . Para lograr esta discretización, asumimos la existencia de funciones de interpolación N1 (x) y
N2 (x) (también conocidas como funciones de forma o funciones de mezcla) tales que
Nótese que usando los valores extremos de elemento u1 y u2 estamos planteando que podemos hallar
el desplazamiento u(x) de cualquier punto interno haciendo uso de las funciones N1 (x) y N2 (x); lo que
sin duda representa el planteamiento de un proceso de interpolación matemática. De ahı́ el nombre de
funciones de interpolación para las funciones de la variable independiente: N1 y N2 .
Debe darse énfasis a que, aunque se indica una igualdad en la Ecuación (2.16), la relación plan-
teada para los elementos finitos en general, es una aproximación. Para el elemento–barra la relación,
sin embargo, de hecho es exácta. Para determinar las funciones de interpolación, requerimos que los
valores de frontera lı́mite de u(x) (los desplazamientos nodales) sean idénticamente satisfechos por la
discretización planteada, de modo que
Las Ecuaciones (2.16) y (2.17) conducen a las siguientes condiciones lı́mite (nodales)
N1 (x) = a0 + a1 x (2.20)
N2 (x) = b0 + b1 x (2.21)
donde los coeficientes polinómicos serán determinados mediante el cumplimiento de las condiciones de
borde lı́mites (nodales). Notamos aquı́ que existe un número ilimitado de expresiones matemáticas que
podrı́an asumirse como funciones de interpolación, a sóla condición que satisfagan las condiciones de
borde requeridas. Las razones para escoger una forma lineal para las funciones de interpolación en este
caso particular se explican en detalle en el Capı́tulo 6.
La aplicación de las condiciones representadas por la Ecuación (2.18) nos proporciona: a0 = 1,
b0 = 0, mientras que la Ecuación (2.19) dá como resultados a1 = −(1/L) y b1 = 1/L. Por consiguiente,
las funciones de interpolación son
Como se encontrará subsecuentemente más conveniente, la Ecuación (2.23) puede expresarse en forma
matricial como
u1
u(x) = N1 (x) N2 (x) = [N ]{u} (2.24)
u2
donde [N ] es la matriz fila de funciones de interpolación y {u} es la matriz columna (vector) de
desplazamientos nodales.
Habiendo expresado el campo de desplazamientos en términos de las variables nodales, el resto de
la tarea obligada es determinar la relación existente entre los desplazamientos nodales y las fuerzas
aplicadas en dichos puntos extremos del elemento, para obtener la matriz de rigidez del elemento barra.
32 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
u2 − u1 δ
x = = (2.28)
L L
que muestra que el elemento barra es un elemento de deformación constante. Esto está de acuerdo
con la teorı́a de resistencia de materiales: El elemento tiene área de sección transversal constante, y
está sujeto a fuerzas constantes aplicadas en los puntos extremos (los nodos), ası́ que la tensión interna
no varı́a a lo largo de la longitud. La tensión axial interna, por la ley de Hocke, es entonces
u2 − u1
σx = Ex = E (2.29)
L
y la fuerza axial asociada es luego
EA
P = σx A = (u2 − u1 ) (2.30)
L
Tomando cuidado para observar la correcta convención de signo algebraico, la Ecuación (2.30) es
ahora usada para relacionar las fuerzas nodales aplicadas f1 y f2 a los desplazamientos nodales u1 y u2 .
Aquı́ podrı́amos apelar a la analogı́a de comportamiento mecánico del elemento–barra con el elemento
resorte (equivalente) el cual es mostrado en la Figura 2.1(a). Observando que, si la Ecuación (2.30)
tiene un signo positivo, el elemento está en solicitación de tracción y la fuerza nodal f2 debe estar en
la dirección de la coordenada positiva mientras que la fuerza nodal f1 debe ser igual en magnitud y de
sentido opuesto para preservar el equilibrio; por consiguiente,
AE
f1 = − (u2 − u1 ) (2.31a)
L
EA
f2 = (u2 − u1 ) (2.31b)
L
2.3. EL ELEMENTO BARRA 33
La comparación de la Ecuación (2.32) y la Ecuación (2.4) muestra que la matriz de rigidez para el
elemento–barra está dada por
E A 1 −1
[ke ] = (2.33)
L −1 1
Como es el caso con el resorte lineal, observamos que la matriz de rigidez para el elemento barra es
simétrica, singular, y de orden 2×2 en correspondencia con los dos desplazamientos nodales o grados de
libertad. Debe darse énfasis a que la matriz de rigidez dada por la Ecuación (2.33) está expresada en
el sistema coordenado de elemento (sistema local) que en este caso es uni–dimensional. La aplicación
de la formulación de éste elemento al análisis de estructuras bi y tri–dimensionales es considerada en
el próximo capı́tulo.
Ejemplo 2.4.
La Figura 2.7(a) muestra una barra elástica adelgazada sujeta a una carga de tracción P aplicada en
un extremo, estando fijo el otro extremo. El área de la seción–transversal varı́a linealmente desde A0
en el extremo fijo x = 0, hasta A0 /2 en x = L. Calcular el desplazamiento del extremo libre de la
barra (a) modelando este cuerpo como un solo elemento que tiene área de sección transversal igual al
área de la barra real en su punto medio a lo largo de su longitud, (b) usando dos elementos barra de
igual longitud y similarmente evaluando el área en el punto medio de cada uno de los segmentos, y (c)
usando la integración para obtener la solución exacta.
(x)
A0 1 A0 1 A0
x x x
7A0/8 L/2 A(x)
L 3A0/4 L 2
L x
5A0/8 L/2
A0/2 A20/2 3
P P P P
(a) Sistema original (b) Modelo de 1 elemento (c) Modelo de 2 elementos (d) Tensión normal
> Solución
(a) Para un solo elemento, el área de sección–transversal es 3A0 /4 y el “coeficiente de rigidez equiva-
lente” de resorte es
AE 3A0 E
k= =
L 4L
y la ecuación matricial gobernante del comportamiento de elemento (o estructura, en éste caso) es
3A0 E 1 −1 U1 F1
=
4L −1 1 U2 P
El elemento y los puntos nodales del mismo son como se muestran en la Figura 2.7(b). Los desplaza-
mientos de estos puntos (grados de libertad) se asumen estar definidos según la dirección x. Aplicando
34 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
4P L PL
U2 = = 1,333
3A0 E A0 E
como el desplazamiento en x = L.
(b) Dos elementos de igual longitud L/2 con los puntos nodales asociados a cada uno de ellos se
muestran en la Figura 2.7(c). Los desplazamientos de estos puntos (grados de libertad) se asumen
estar definidos según la dirección x. Para el elemento 1, A1 = 7A0 /8, por lo que el coeficiente de rigidez
de resorte equivalente resulta
A1 E 7A0 E 7A0 E
k1 = = =
L1 8(L/2) 4L
mientras que para el elemento 2 tenemos
Puesto que ninguna carga está aplicada en el centro de la barra, las ecuaciones de equilibrio para
el sistema de dos elementos es
k1 −k1 0 U1 F1
−k1 k1 + k2 −k2 U2 = 0
0 −k2 k2 U3 P
k2
U2 = U3
k1 + k2
Reemplazando en la segunda ecuación y evaluando para el desplazamiento incógnita, se comprueba
que resulta
k1 + k2 48P L PL
U3 = P = = 1,371
k1 k2 35A0 E A0 E
(c) Para obtener la solución exácta, nos referimos a la Figura 2.7(d), en la que mostramos un diagrama
de cuerpo–libre de una sección de la barra entre una posición arbitraria x y el extremo x = L de ella.
Para el equilibrio,
x
σx A = P y, puesto que A = A(x) = A0 1 −
2L
la variación de tensión axial interna a lo largo de la longitud de la barra se describe por
P
σx = x
A0 1 − 2L
Luego, aplicando la ley de Hocke, la deformación unitaria axial vendrá determinada por
σx P
x = = x
E E A0 1 − 2L
2.3. EL ELEMENTO BARRA 35
Pero, por definición sabemos que la deformación unitaria se relaciona con el campo de desplazamientos
interno u = u(x) mediante la conocida relación: x = du/dx. Ahora, como la barra está fijada en x = 0,
el desplazamiento en cualquier posición genérica x estará determinado por
Z x Z x x
P dx 2P L
u(x) = x dx = = − [ln(2L − x) − ln(2L)]
x
E A 1 − 2L E A
0 0 0 0 0
2P L 2P L 2L
=− [ln(2L − x) − ln(2L)] = ln
E A0 E A0 2L − x
Cuando evaluamos esta función en la posición correspondiente al extremo libre del cuerpo, donde
está aplicada la carga, obtenemos el valor del desplazamiento en dicho punto que en este caso corres-
ponde a la deformación neta del cuerpo tronco–cónico; entonces
2P L 2L
δ = u(x) = u(L) = ln
x=L E A0 2L − L
2P L PL
= ln 2 = 1,386
E A0 E A0
1.4 2.0
1.8 Soluc. Exácta
1.2 Soluc. Exácta Un elemento
Un elemento 1.6 Dos elementos
1.0 Dos elementos
1.4
0 0 A 0
P
1.2
0.8
u(x)
1.0
0.6 0.8
0.4 0.6
0.4
0.2
0.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x x
L L
(a) Campo de desplazamiento interno (b) Tensión normal axial interna
La comparación de los resultados de las partes (b) y (c) revela que la solución de dos elementos
exhibe un error porcentual relativo de sólo aproximadamente 1 por ciento comparado con la solución
exacta obtenida mediante la teorı́a de mecánica de los materiales. La Figura 2.8(a) muestra la variación
del desplazamiento a lo largo de la longitud para las tres soluciones obtenidas. Notemos que la solución
proporcionada por el modelo de dos elementos finitos es más cercana a la solución exacta que la
solución obtenida utilizando un solo elemento finito, ésto en el dominio de definición del problema;
pero, al comparar la solución proporcionada para la deformación neta (el desplazamiento del extremo
libre) con la solución exácta, se calcula que el error porcentual relativo es de 3,8 %, lo cual no es un
error demasiado grosero.
Es sumamente importante notar, sin embargo, que la tensión normal axial calculada para las so-
luciones de elemento finito varı́an significativamente respecto de aquellas que describen la solución
exacta. La tensión normal axial para la solución de dos–elementos se muestra en la Figura 2.8(b), jun-
to con la tensión calculada proveniente de la solución exacta. Note particularmente la discontinuidad
de valores de tensión calculados para los dos elementos en el punto de conexión nodal entre ellos. Esto
es tı́pico de las variables derivadas, o secundarias, como la deformación unitaria o la tensión interna,
36 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
cuando son calculadas en el método de elemento finito. A medida que más y más elementos pequeños
sean usados en el modelo, los valores de tales discontinuidades disminuyen, indicando la convergencia
hacia la solución exacta. En los análisis estructurales, el usuario del elemento finito está más a menudo
interesado en las tensiones que los desplazamientos, de aquı́ es esencial que se supervise la convergencia
de las variables secundarias. >
El ejemplo que resolvimos anteriormente claramente nos muestra que el método de elemento finito
es un método aproximado, y como tal posee falencias inherentes a su propia naturaleza. Sin embargo,
es posible corregir los errores que se presentan apelando a ciertas modificaciones en el proceso de
formulación de los elementos. Para el ejemplo anterior, se lograrı́a un resultado más óptimo si el
campo de desplazamientos interno en cada uno de los elementos se asume con variación cuadrática
respecto a la posición; con ello, las deformaciones unitarias y también las tensiones normales internas
tendrı́an variación lineal, porque las deformaciones se relacionan a los desplazamientos mediante la
primera derivada y las tensiones son diréctamente proporcionales a las deformaciones. Este modo de
aplicar el método de elemento finito conduce a los llamados elementos de orden superior, tema que
será tratado en capı́tulos posteriores.
F~ = Fx ı̂ + Fy ̂ + Fz k̂ y ~ = dxı̂ + dy ̂ + dz k̂
dr
el trabajo mecánico valorado en este sistema coordenado, viene determinado por la relación
Z x2 Z y2 Z z2
W = Fx dx + Fy dy + Fz dz (2.35)
x1 y1 z1
k
F
Fuerza, F
L0 F0
A=W=Ue
k
k F0
1
L0
0 Deformación, 0
(a) Carga sobre el resorte (b) Respuesta estática (c) Trabajo – Energı́a
Consideremos un cuerpo elástico sujeto a N fuerzas Fj para el cual la energı́a total de deformación
se expresa como
N Z δj
X
Ue = W = Fj dδj (2.40)
j=1 0
donde se asume que la fuerza original Fi es constante durante el cambio infinitesimal efectuado. El
término integral en la Ecuación (2.41) involucra el producto de cantidades infinitesimales, que producen
un infinitésimo de segundo orden que puede despreciarse por su orden de magnitud comparado con
el término que lo antecede en la misma ecuación. Despreciando el término identificado, la relación
anterior se reduce a
∆Ue
= Fi
∆δi
la cual en el lı́mite cuando ∆δi → 0, se convierte en
∂Ue
= Fi (2.42)
∂δi
Aplicando el teorema de Castigliano con respecto a cada uno de los desplazamientos, se obtiene
∂Ue AE
= (u1 − u2 ) = f1 (2.44a)
∂u1 L
∂Ue AE
= (u2 − u1 ) = f2 (2.44b)
∂u2 L
las cuales como puede observarse son idénticas a las Ecuaciones (2.31).
Nótese que en el ámbito del método de elemento finito, la formulación energética es una metodologı́a
alternativa al método de equilibrio de fuerzas que inicialmente aplicamos en los desarrollos iniciales
de éste capı́tulo. Está demostrado que los métodos de análisis basados en una concepción Newtoniana
y los métodos de análisis basados en una concepción de tipo energético, necesariamente conducen a
soluciones idénticas cuando son aplicados a un mismo problema, preservando obviamente las hipótesis
asumidas para ambos casos.
El primer teorema de Castigliano también es aplicable a los desplazamientos angulares. En el
caso de rotación, la derivada parcial de la energı́a de deformación con respecto a un desplazamiento
angular es igual al momento/torque aplicado al punto concerniente en el sentido de la rotación. El
ejemplo siguiente ilustra la aplicación de este concepto en términos de un miembro simple torsional,
más propiamente un eje de sección transversal circular.
2.4. ENERGÍA DE DEFORMACIÓN ELÁSTICA 39
Ejemplo 2.5.
Un eje circular sólido de radio R y longitud L, se somete a la acción de un torque constante T . El
eje está fijo en un extremo, como es mostrado en la Figura 2.10. Formule la energı́a de deformación
elástica en términos del ángulo de torsión en x = L, y mostrar que el primer teorema de Castigliano
proporciona la expresión correcta para el torque aplicado
L T
> Solución
De la teorı́a de mecánica de sólidos, la tensión cortante en cualquier sección transversal a lo largo de
la longitud del miembro está determinada por
Tr
τ=
J
donde r es la distancia radial variable medida desde el eje axial que atravieza el centro de la superficie
circular, y J es el momento polar de inercia de la sección transversal.
Considerando un comportamiento elástico del material, tendremos
τ
γ=
G
donde G es el módulo de elasticidad transversal (módulo de Young) del material, γ es la deformación
angular unitaria (desplazamiento angular por unidad de longitud); y entonces la energı́a de deformación
de estado torsional es por analogı́a con la energı́a de deformación en estado axial de solicitación
Z Z L Z
1 1 Tr Tr
Ue = τ γ dV = dA dx
2 V 2 0 A J GJ
LZ
T2 T2 L
Z
= r2 dA dx =
2 J2 G 0 A 2GJ
Nuevamente, invocando los resultados que proporciona la mecánica de sólidos, el ángulo de torsión
al final del miembro se conoce que es
TL
θ=
GJ
40 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
∂Ue GJ θ
=T =
∂θ L
que es exactamente la relación deducida por la teorı́a de la mecánica de sólidos.
El lector puede pensar que usamos un razonamiento circular en este ejemplo (volvimos al mismo
lugar desde el cual empezamos), ya que que utilizamos muchos resultados previamente conocidos. Sin
embargo, la formulación de la energı́a de deformación debe estar basada en relaciones conocidas de
tensión y deformación unitaria; y la aplicación del teorema de Castigliano es, de hecho, un concepto
diferente. >
Para sistemas linealmente elásticos, la formulación de la energı́a de deformación en términos de los
desplazamientos es relativamente dirécta. Como fué establecido previamente, la energı́a de deformación
para un sistema elástico es una función cuadrática de los desplazamientos. La naturaleza cuadrática
es explicada de modo simplista por los hechos que, en deformación elástica, la tensión interna es
proporcional a la fuerza (en este caso, al momento o torque), y la deformación unitaria es proporcional
a la tensión, y la deformación unitaria es proporcional al desplazamiento (o rotación). Y, desde que la
energı́a de deformación elástica es igual al trabajo mecánico expendido, resulta una función cuadrática.
Por consiguiente, la aplicación del primer teorema de Castigliano resulta en ecuaciones algebráicas
lineales que relacionan los desplazamientos a las fuerzas aplicadas. Esta declaración sigue del hecho que
la derivada de un término cuadrático es lineal. Los coeficientes de los desplazamientos en las ecuaciones
resultantes son los componentes de la matriz de rigidez del sistema para la cual la función de energı́a
de deformación se escribe. Tal concepción de análisis basado en conceptos energéticos es el método
más simple y dirécto para establecer la matriz de rigidez de muchos elementos finitos estructurales.
Ejemplo 2.6.
La Figura 2.11 muestra un sistema compuesto de cuatro resortes conectados en serie y paralelo.
(a) Aplicar el primer teorema de Castigliano para obtener la matriz de rigidez del sistema. Los miembros
verticales ubicados en los nodos 2 y 3 se consideran rı́gidos.
(b) Resolver el modelo matemático planteado para los desplazamientos activos y la fuerza transmitida
hacia el apoyo en el nodo 1 si
k2
k1 k3 F4
F2
1 2 3 4
k2
> Solución
(a) La energı́a total de deformación del sistema de cuatro resortes mostrado en la figura, se expresa en
2.4. ENERGÍA DE DEFORMACIÓN ELÁSTICA 41
términos de los desplazamientos nodales y los diversos coeficientes de rigidez de los mismos, mediante
la relación
Ue = 12 k1 (U2 − U1 )2 + 2 [ 21 k2 (U3 − U2 )2 ] + 12 k3 (U4 − U3 )2
Aplicando el teorema de Castigliano, usando cada uno de los desplazamientos nodales en turno sucesivo
∂Ue
= F1 = k1 (U2 − U1 )(−1) = k1 U1 − k1 U2
∂U1
∂Ue
= F2 = k1 (U2 − U1 ) + 2k2 (U3 − U2 )(−1) = −k1 U1 + (k1 + 2k2 )U2 − 2k2 U3
∂U2
∂Ue
= F3 = 2k2 (U3 − U2 ) + k3 (U4 − U3 )(−1) = −2k2 U2 + (2k2 + k3 )U3 − k3 U4
∂U3
∂Ue
= F4 = k3 (U4 − U3 ) = −k3 U3 + k3 U4
∂U4
Este sistema de ecuaciones lineales algebráico puede escribirse en forma de arreglo matricial según
k1 −k1 0 0 U1 F1
−k1 k1 + 2k2 −2k 2 0 U2 F2
=
0 −2k2 2k2 + k3 −k3 U3 F3
0 0 −k3 k3 U4 F4
y, la matriz de rigidez que es mostrada en la ecuación precedente fué hallada por aplicación del teorema
de Castigliano.
(b) Substituyendo los valores numéricos especificados en el planteamiento del problema, la ecuación
matricial gobernante del comportamiento estático del sistema llega a ser
4 −4 0 0 0
F1
−4 16 −12 0 U2 −30
0 −12 15 −3 U3 0 =
0 0 −3 3 U4 50
que resolveremos manipulando las ecuaciones para convertir la matriz de coeficientes (la matriz de rigi-
dez) en una forma triangular–superior; es decir, que todos los términos debajo de la diagonal principal
se vuelvan ceros.
Paso 1 Multiplicar la primera ecuación (fila) por 12, multiplicar la segunda ecuación (fila) por 16,
agregar las dos filas y reemplazar la segunda ecuación con la ecuación resultante de la suma,
para obtener
16 −12 0 U2 −30
0 96 −48 U3 = −360
0 −3 3 U4 50
Paso 2 Multiplicar la tercera ecuación por 32, añadir a la segunda ecuación, y reemplazar la tercera
ecuación con el resultado. Esto proporciona la forma triangularizada deseada
16 −12 0 U2 −30
0 96 −48 U3 = −360
0 0 48 U4 1240
42 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
En esta forma, las ecuaciones pueden resolverse ahora por “sustitución regresiva”, y será encontrado
que, a cada paso, hay sólo un valor desconocido. En este caso, la secuencia es
1240
U4 = = 25,83 mm
48
1
U3 = [ −360 + 48(25,83) ] = 9,17 mm
96
1
U2 = [ −30 + 12(9,17) ] = 5,0 mm
16
La fuerza de reacción en el nodo 1 se obtiene a partir de la ecuación de restricción que fué desechada
en la solución de los grados de libertad no–restringidos
F1 = −4 U2 = −4 (5,0) = −20 N
y, podemos observar que el sistema se encuentra en equilibrio estático, desde que la suma de fuerzas
externas actuantes sobre el sistema (incluı́da la reacción de apoyo) es idénticamente cero, como se
requiere. >
En este Ejemplo hemos mostrado una técnica de solución de la ecuación matricial gobernante del
comportamiento estático del sistema, la cual consiste en efectuar transformaciones con las filas de la
matriz de coeficientes (la matriz de rigidez global), de modo que al final de las mismas la apariencia
de dicho arreglo sea de forma triangular–superior; entonces luego dicho sistema ası́ transformado es
fácilmente resuelto evaluando el conjunto de ecuaciones de “atrás para adelante”, es decir en forma
inversa (o regresiva) sustituyendo en cada nueva instancia la solución hallada en la instancia anterior.
Π = Ue + UF (2.45)
donde debemos hacer notar que el término “fuerzas” externas es generalizado, porque también debe
incluir torques o momentos si estos existen.
En este texto, sólo trataremos con sistemas elásticos sujetos a sistemas de fuerzas conservativas. Una
fuerza conservativa se define como aquella que realiza trabajo mecánico de magnitud independiente
del camino de movimiento, de modo tal que el trabajo en este caso es reversible o recuperable. El
ejemplo más común de una fuerza no–conservativa, también llamada disipativa, es la fuerza de fricción
2.5. ENERGÍA POTENCIAL MÍNIMA 43
Π = Ue − W (2.47)
Como mostramos en los ejemplos siguientes y las aplicaciones a la mecánica del medio sólido de-
formable en el Capı́tulo 9, el término de energı́a de deformación Ue es una función cuadrática de los
desplazamientos del sistema y el término de trabajo mecánico W es una función lineal de los despla-
zamientos. Rigurosamente, la minimización de la energı́a potencial total es un problema de cálculo
variacional [30]. Nosotros no suponemos que el público intencional de este texto esté familiarizado con
el cálculo de variaciones. Más bien, aquı́ simplemente imponemos la minimización según el principio de
cálculo de derivación de funciones de variables múltiples. Si tenemos una expresión de energı́a potencial
total que es una función, digamos, de N desplazamientos Ui , i = 1, . . . , N ; ésto es
Π = Π(U1 , U1 , . . . , UN ) (2.48)
entonces la energı́a potencial total será minimizada si hacemos cumplir las relaciones
∂Π
=0 i = 1, 2, . . . , N (2.49)
∂Ui
La Ecuación (2.49) se mostrará como relación que representa a las N ecuaciones algebráicas que
conforman la aproximación de elemento finito a la solución de las ecuaciones diferenciales gobernantes
de la respuesta de un sistema estructural debida a la solicitación externa impuesta sobre dicho sistema.
Ejemplo 2.7.
Repetir el Ejemplo 2.6, hallando la solución aplicando el principio de energı́a potencial mı́nima.
> Solución
Del planteamiento de solución del Ejemplo mencionado anterior, la energı́a de deformación elástica es
Ue = 1
2 k1 (U2 − U1 )2 + 2 [ 21 k2 (U3 − U2 )2 ] + 1
2 k3 (U4 − U3 )2
UF = −W = −F1 U1 − F2 U2 − F3 U3 − F4 U4
Π = Ue − W = 1
2 k1 (U2 − U1 )2 + 2 [ 12 k2 (U3 − U2 )2 ]
+ 1
2 k3 (U4 − U3 )2 − F1 U1 − F2 U2 − F3 U3 − F4 U4
aplicando en secuencia la ecuación recursiva anterior incrementando el ı́ndice genérico desde su valor
inicial, tendremos
∂Π
= k1 (U2 − U1 )(−1) − F1 = k1 U1 − k1 U2 − F1 = 0
∂U1
∂Π
= k1 (U2 − U1 ) + 2 k2 (U3 − U2 )(−1) − F2
∂U2
= −k1 U1 + (k1 + k2 )U2 − 2 k2 U3 − F2 = 0
∂Π
= 2 k2 (U3 − U2 ) + k3 (U4 − U3 )(−1) − F3
∂U3
= −2 k2 U2 + (2 k2 + k3 )U3 − k3 U4 − F3 = 0
∂Π
= k3 (U4 − U3 ) − F4 = −k3 U3 + k3 U4 − F4 = 0
∂U4
El sistema de ecuaciones lineales algebráico, cuando se ordena en forma matricial dá como resultado
k1 −k1 0 0 U1 F1
−k1 k1 + 2k2
−2k 0 U F2
2 2
=
0 −2k2 2k2 + k3 −k3 U3 F3
0 0 −k3 k3 U4 F4
y podemos ver que resulta idéntico al resultado previamente obtenido. Consecuentemente ya no re-
solvemos el ejemplo numéricamente en virtud que este procedimiento yá fué efectuado en el ejemplo
previo a éste. >
En este punto, como complemento de carácter teórico, re–examinaremos la ecuación de energı́a
global del Ejemplo 2.7 para desarrollar una forma más general, la cual será de valor significativo en
los sistemas más complicados a ser discutidos en los capı́tulos siguientes. El vector de desplazamiento
global, o de sistema, es
U1
U2
{U } =
U3
U4
y, como fué deducido, la matriz de rigidez del sistema (o, matriz de rigidez global) es
k1 −k1 0 0
−k1 k1 + 2k2 −2k2 0
[K ] =
0 −2k2 2k2 + k3 −k3
0 0 −k3 k3
este procedimiento, podemos indicar que la matriz resultado del producto triple especificado en la
relación anterior, representa la energı́a de deformación de cualquier sistema elástico.
Ue = 12 {U }T [ K ]{U } (2.50)
2.6. Resumen
Dos elementos mecánicos lineales, el resorte elástico idealizado y un miembro elástico en solicitación
de tracción/compresión (el elemento–barra) han sido usados para introducir los conceptos básicos
involucrados en la formulación de las ecuaciones que gobiernan el comportamiento de un elemento
finito. Las ecuaciones del elemento se obtienen por una concepción de equilibrio estático dirécto y
también alternativamente por un método basado en la energı́a de deformación elástica, usando el
primer teorema de Castigliano. También se introdujo el principio de energı́a potencial mı́nima como una
alternativa equivalente de formulación de problemas estructurales a ser resueltos mediante un modelo
de elemento finito. El próximo capı́tulo muestra cómo el elemento–barra uni–dimensional puede usarse
para demostrar los procedimientos de ensamble en modelos compuestos de varios elementos finitos, en
el contexto de estructuras simples que se desarrollan en dos y tres dimensiones del espacio.
k1
Problemas propuestos
2.1. — 2.2. — 2.3. Para cada conjunto de resortes conectados del modo aquı́ mostrado, según
se aprecia en las Figuras P2.1 – P2.2 – P2.3; determine la matriz de rigidez global usando el
procedimiento de ensamblaje de sistemas desarrollado en la Sección 2.2.
k3
k1 k2 k3 k1 k2
4
1 2 3
1 2 3 4 k3
Figura P2.3
2.4. Para el conjunto de resortes conectados como se muestra, determinar la fuerza F3 requerida
para desplazar el nodo 2 un monto de δ = 0.75 pulg. hacia la derecha. También calcular el
desplazamiento del nodo 3. Considere que,
k1 = 50 lb/pulg y k2 = 25 lb/pulg
k1 k2 F3
1 2 3
Figura P2.4
46 CAPÍTULO 2. ANÁLISIS UNI–DIMENSIONAL
2.5. En el ensamble de resortes conectados como se muestra en la Figura, las fuerzas F2 y F4 van
a ser aplicadas de modo que la fuerza resultante sobre el elemento 2 sea cero y el nodo 4 se
desplace un monto δ = 1 pulg. Determinar (a) los valores requeridos de las fuerzas F2 y F4 , (b)
el desplazamiento del nodo 2, y (c) la fuerza de reacción en el nodo 1.
k1 F2 k2 k3
1 2 3 4 F4
k1 k3 30 lb/in k2 40 lb/in
Figura P2.5
2.6. Verificar la matriz de rigidez global hallada en el Ejemplo 2.3 usando (a) un procedimiento de
ensamble dirécto y (b) el primer teorema de Castigliano.
2.7. Dos vagones de un tren son conectados por el arreglo de resortes mostrados en la Figura. (a) De-
terminar la serie completa de ecuaciones de equilibrio para el sistema en la forma [ K ]{U } = {F }.
(b) Si k = 50 lb/pulg, F1 = 20 lb, y F2 = 15 lb, calcular el desplazamiento de cada vagón y la
fuerza actuante en cada resorte.
F1 2k
F2
k k
2k
Figura P2.7
2.8. Usar el primer teorema de Castigliano para obtener la representación matricial de las ecuaciones
de equilibrio estático para el sistema de resortes mostrado en la Figura.
k1 F2 k2 F3 k3 F4 k4
1 2 3 4 5
Figura P2.8
2.10. Una barra de acero sujeta a solicitación de compresión es modelada mediante dos elementos –
barra, como se muestra en la Figura. Determinar los desplazamientos nodales y la tensión nor-
mal en cada elemento. Que otras variables importantes se podrı́an calcular en este problema ?.
0.5 m 0.5 m
12 kN
1 1 2 2 3
8 kN/m
E = 207 GPa A = 500 mm2
Figura P2.10
Problemas propuestos 47
2.11. La Figura muestra el ensamble de dos elementos – barra, construidos de diferentes materiales.
Determinar los desplazamientos nodales, la tensión normal en cada elemento, y la fuerza de
reacción de apoyo.
E1 , L1 , A1
E2 , L2 , A2
2 X 104 lb
1 2
3
2 7
1 7 E2= 10 lb/in2
E1= 1.5 X 10 lb/in2 L2= 20 in
L1= 20 in A2= 2.25 in2
A1= 4 in2
Figura P2.11
2.12. Obtener la solución del modelo de cuatro elementos para la barra tronco–cónica analizada en
el Ejemplo 2.4. Bosquejar la gráfica de la tensión normal en los elementos versus la solución
exacta. Usar los siguientes valores numéricos:
2
E = 107 lb/pulg A0 = 4 pulg2 L = 20 pulg P = 4000 lb
2.13. Un peso W es suspendido de un plano vertical conectado a un resorte lineal que tiene coeficiente
de rigidez k conocido. Mostrar que la posición de equilibrio estático corresponde a un valor de
energı́a potencial total mı́nimo.
2.14. Para un elemento – barra se propone discretizar el campo de desplazamientos interno mediante
la relación funcional
u(x) = N1 (x) u1 + N2 (x) u2
con las siguientes funciones de interpolación
πx
N1 (x) = cos
2L
πx
N2 (x) = sin
2L
Son éstas funciones de interpolación válidas ?. (Sugerencia: Considere las variaciones de tensión
y deformación unitaria al interior del elemento).
2.15. El elemento torsional mostrado en la Figura tiene sección transversal circular y comportamiento
elástico. Los desplazamientos nodales son las rotaciones θ1 y θ2 , y las cargas nodales asociadas
son los torques aplicados T1 y T2 . Usar el principio de energı́a potencial mı́nima para deducir
las ecuaciones de comportamiento estático de elemento en forma matricial.
1 ,T
1 R
L
2 ,T
2
Figura P2.15
El método dirécto de rigidez
Capítulo 3
En el Capı́tulo anterior efectuamos una introducción a la formulación del método de elemento fi-
nito utilizando para ello los elementos uni–dimensionales: resorte, barra y eje, efectuando también un
análisis de sistemas pequeños conformados mediante interconexión de éstos elementos. Obtuvimos la
ecuación gobernante del sistema por aplicación de condiciones de equilibrio nodal. También efectua-
mos una introdución a los denominados métodos energéticos de formulación, y vimos que los mismos
producen soluciones idénticas a aquellas obtenidas mediante las leyes de equilibrio estático.
En éste capı́tulo mostraremos el modo en el que se pueden analizar sistemas estructurales de relativo
mayor tamaño compuestos particularmente de elementos barra interconectados entre sı́ en sus puntos
extremos. Efectuaremos una discusión mucho más detallada de las ecuaciones de equilibrio de los
puntos relevantes (nodos) de una estructura. Tambı́en mostraremos los procesos de transformación de
los sistemas de ejes coordenados utilizados en el análisis, que posteriormente nos permitirá efectuar el
procedimiento de ensamblaje estructural que se traduce en el planteamiento de la ecuación gobernante
del comportamiento del sistema. Mostraremos como se aplican las condiciones de borde o contorno
impuestas sobre el modelo matemático elaborado, y finalmente efectuaremos una descripción de la
obtención de variables secundarias de interés como son las deformaciones, reacciones de apoyo, y las
tensiones internas que aparecen en la estructura deformada.
3.1. Introducción
Los elementos que tienen geometrı́a y comportamiento uni–dimensional (elementos lineales simples)
discutidos en el Capı́tulo 2 introdujeron los conceptos de nodos, desplazamientos nodales, y matrices
de rigidez de elemento. En este capı́tulo, es considerada la creación de un modelo de elemento finito
de un sistema mecánico compuesto de cualquier número de elementos. La discusión se limita a las
estructuras tipo cercha o armadura que nosotros definimos como estructuras compuestas de miembros
elásticos rectos sujetos sólamente a fuerzas axiales. La satisfacción de esta restricción requiere que
todos los miembros de la estructura sean elementos barra y que los mismos se inter–conecten mediante
uniones articuladas tal que cada elemento sea libre de rotar alrededor de la junta de unión. Aunque
49
50 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
3 7
10 3 7
Y 2
4 6 8
2 6
X 4
1 5 9
(a) Interconexión de elementos (b) Detalle de conexión
1. El desplazamiento nodal del extremo cada uno de los elementos conectados en un punto común
debe ser el mismo que el desplazamiento del nodo de conexión al que concurren en el sistema de
coordenadas globales; la formulación matemática, como se verá, refuerza este requerimiento (al
que se denomina compatibilidad de desplazamiento).
2. Las caracterı́sticas fı́sicas (en este caso, la matriz de rigidez y la fuerza axial) de cada elemento debe
ser transformado, matemáticamente, al sistema coordenado global para representar las propiedades
estructurales en el sistema global en un marco de referencia matemáticamente consistente.
3. Los parámetros individuales de elemento de preocupación a ser evaluados (por ejemplo para el ele-
mento barra, la tensión axial) son determinados después de haber hallado la solución del problema
en el sistema de coordenadas global por la transformación de resultados retornando al marco de
referencia local de elemento (post–procesamiento).
Por qué estamos basando la formulación teórica sobre los desplazamientos?. Generalmente, un di-
seño ingenieril está más interesado en la tensión interna a la cual cada miembro de la cercha está sujeto,
para comparar este valor de tensión con una propiedad material conocida, como la tensión lı́mite de
fluencia del material (valor coincidente con su resistencia mecánica). La comparación de los valores de
tensión calculados con las propiedades materiales, puede llevar a cambios del diseño en cuanto al ma-
terial a utilizarse o las propiedades geométricas de los elementos individuales (en el caso del elemento
barra, el área de su sección transversal). La respuesta a las diversas preguntas radica en la naturaleza
de los problemas fı́sicos que se presentan. Es mucho más fácil estimar la solicitación (las fuerzas y mo-
mentos aplicados) a la que una estructura está sujeta, que las deformaciones ocurridas por esta causa
en tal estructura. Si las cargas externas se especifican, las relaciones entre la solicitación aplicada y
los desplazamientos producidos son formuladas en términos de la matriz de rigidez y resolvemos las
ecuaciones planteadas ası́ para los desplazamientos incógnita. La substitución–regresiva de los despla-
zamientos en las ecuaciones de elemento individual nos dan entonces las deformaciones y las tensiones
en cada elemento como es deseado. Éste es el método de rigidez y se usa exclusivamente en este texto.
3.2. ECUACIONES DE EQUILIBRIO NODAL 51
mación unitaria, la tensión, y las fuerzas de reacción en las ubicaciones restringidas (condiciones de
borde) donde la estructura se conecta hacia sus apoyos.
Se aconseja al lector que tome nota que sólo usamos el término variable secundaria en el sentido
matemático; la deformación y tensión sólo son secundarias en el sentido que son evaluadas después que
los valores de la solución general para los desplazamientos hayan sido calculados. En el ámbito de la
mecánica estructural, los valores de tensión y deformación son de importancia primaria en el diseño.
U6
F 3Y U5
2
2 3 F3 X U4 U6
x
2 U3 U5
1
Y 1
x
1 U2
1 X U1
(a) Interconexión de elementos (b) Especificación de desplazamientos
La conversión de ecuaciones de elemento desde las coordenadas propias del mismo (sistema local
coordenado) hacia las coordenadas globales y el ensamble de las ecuaciones de equilibrio globales se
describe primero en el caso bidimensional con referencia a la Figura 3.2(a). La figura muestra una
simple cercha bidimensional compuesta de dos miembros estructurales unidos por conexiones articu-
ladas y sujetos a las fuerzas externas aplicadas mostradas. Se toman las conexiones articuladas como
los nodos de los dos elementos considerados; se numeran los nodos y los elementos, ası́ como también
se muestra el sistema coordenado global seleccionado para efectuar el análisis. Los desplazamientos
globales correspondientes se muestran en la Figura 3.2(b). La convención usada aquı́ para los despla-
zamientos globales es que U2i−1 es el desplazamiento en la dirección global X del nodo i, y U2i es el
desplazamiento del nodo i en la dirección global Y . La convención adoptada no es de ningún motivo
restrictivo; la convención se selecciona tal que los desplazamientos en la dirección del eje X global sean
numerados de manera impar y los desplazamientos en la dirección del eje global Y sean numerados
de modo par. (Al usar software de elemento finito, usted hallará que los desplazamientos globales se
denotan como U V W ; o también como UX UY UZ ; etc). El ángulo de orientación espacial de elemento
θ se mide como positivo desde el eje X del sistema coordenado global, hacia el eje x del sistema coor-
denado local propio del elemento como se muestra también en la Figura 3.2(a). Los números asignados
a los nodos se muestran encerrados en cı́rculos
, i en cambio los números asignados a los elementos se
los muestra encerrados en rectángulos j . En las diversas variables que se asocian con los elementos,
se usa el número del mismo como superı́ndice encerrado entre paréntesis del modo (j ).
Para obtener las condiciones de equilibrio, se elaboran diagramas de cuerpo–libre de los nodos de
conexión y de los elementos, los cuales son mostrados en la Figura 3.3. Debe notarse que las fuerzas
externas se numeran con la misma convención adoptada para la codificación de los desplazamientos
globales. Para el nodo 1, según la Figura 3.3(a), tenemos las ecuaciones de equilibrio siguientes según
las direcciones globales X e Y , respectivamente:
F2 (2)
f2 F6
2
F1 F3 (2)
f3
1 F4 F5
(1)
f1
(1)
f3
(1)
f3
(2)
f3
(2)
f2 2
(1)
1
f1
Las Ecuaciones (3.1) a (3.3) simplemente representan las condiciones de equilibrio estático desde un
punto de vista de mecánica de cuerpo rı́gido. Asumiendo que las cargas externas F5 y F6 son conocidas,
estas seis ecuaciones de equilibrio nodales formalmente contienen a ocho incógnitas (las fuerzas). Desde
que la cercha del ejemplo es estáticamente indeterminada, podemos invocar condiciones de equilibrio
adicionales aplicables a la estructura en conjunto ası́ como también para aquellos elementos individua-
les que la componen [véase las Figuras 3.3(d) y 3.3(e)] y eventualmente resolver para todas las fuerzas.
Sin embargo, un procedimiento más sistemático se obtiene si la formulación es transformada para
que las incógnitas sean los desplazamientos nodales. Una vez que la transformación ha sido realizada,
encontramos que el número de incógnitas es exactamente el mismo que el número de ecuaciones de
equilibrio nodales. Además, la indeterminación estática se acomoda automáticamente a este procedi-
miento sin ningún tipo de análisis adicional. Como el lector puede recordar del estudio de mecánica
de materiales, la solución de sistemas estáticamente indeterminados requiere la especificación de una
o más relaciones de compatibilidad de desplazamiento; de aquı́, resulta evidente que la formulación de
desplazamiento del método de elemento finito incluye tales situaciones.
Para ilustrar la transformación a los desplazamientos, la Figura 3.2(a) muestra un elemento barra
conectado a los nodos i y j en una posición general en una estructura cercha bidimensional (2–D).
Como resultado de la carga externa actuante sobre la estructura, asumimos que los puntos nodales
extremos sufren desplazamientos en el propio plano, como es mostrado en Figura 3.2(b). Desde que el
54 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
(e)
(e)
u2(e) v2 U4
(e)
j u2 j
j
(e)
U3
deformado
X
x original (e)
U2
(e)
i v1
u1(e) u1(e) i
i (e)
U1
(a) Orientación espacial (b) Desplazamientos generales (c) Desplazamientos globales
elemento debe permanecer conectado a las juntas estructurales, los nodos del elemento conectados en
dichos puntos, deben sufrir los mismos desplazamientos. Esto significa que el elemento no sólo tiene
movimiento axial sino rotación también. Para describir la rotación, agregamos los desplazamientos v1 y
v2 a los nodos 1 y 2 del elemento, respectivamente, en la dirección perpendicular al eje x del elemento.
Debido a la asunción de conexiones articuladas lisas en los puntos de juntura, los desplazamientos
perpendiculares no son asociados con la rigidez del elemento; no obstante, estos desplazamientos deben
existir para que el elemento continúe conectado a la junta estructural para que los desplazamientos
del elemento sean compatibles con (es decir, iguales que) los desplazamientos de la junta. Aunque el
elemento sufre una rotación, para los propósitos de cálculo, en general se asume que el ángulo de la
orientación después de producirse la deformación es igual que en la estructura indeformada. Éste es
resultado de la asunción de deformaciones elásticas pequeñas, y se usa a lo largo del texto.
Para relacionar ahora los desplazamientos nodales de elemento referidos a las coordenadas locales
con los desplazamientos del mismo elemento referidos en las coordenadas globales, la Figura 3.2(c)
muestra los desplazamientos nodales de elemento en el sistema global usando la notación
De nuevo, note el uso de letras mayúscula para las cantidades globales y la notación usada como
exponente para referirse a un elemento individual. Como los desplazamientos nodales deben ser los
mismos en los dos sistemas coordenados, podemos igualar las componentes del vector de desplazamiento
global a los desplazamientos referidos al sistema coordenado de elemento (o sistema local) para obtener
las relaciones
u(e)
1 = U1(e) cos θ + U2(e) sin θ
(3.4a)
v1(e) = −U1(e) sin θ + U2(e) cos θ
u(e)
2 = U3(e) cos θ + U4(e) sin θ
(3.4b)
v2(e) = −U3(e) sin θ + U4(e) cos θ
Como yá lo dijimos, las componentes v de desplazamiento no son asociadas con la rigidez del
elemento, y por ende tampoco se asocian con las fuerzas actuantes sobre el mismo; de modo que
3.2. ECUACIONES DE EQUILIBRIO NODAL 55
δ (e) = u(e)
2 − u(e)
1 = (U3(e) − U1(e) ) cos θ + (U4(e) − U2(e) ) sin θ (3.5)
f (e) = k (e) δ (e) = k (e) {(U3(e) − U1(e) ) cos θ + (U4(e) − U2(e) ) sin θ} (3.6)
Utilizando la Ecuación (3.6) para el elemento 1, acorde con la Figura 3.3(d), y notando que para
éste elemento los desplazamientos propios se relacionan con los desplazamientos globales asignados a
la estructura mediante: U1(1) = U1 , U2(1) = U2 , U3(1) = U5 , y U4(1) = U6 , tenemos la fuerza actuante sobre
el elemento 1 como
f3(1) = −f1(1) = k (1) δ (1) = k (1) {(U5 − U1 ) cos θ1 + (U6 − U2 ) sin θ1 } (3.7)
f3(2) = −f2(2) = k (2) δ (2) = k (2) {(U5 − U3 ) cos θ2 + (U6 − U4 ) sin θ2 } (3.8)
Notemos que, al escribir las anteriores dos ecuaciones últimas, invocamos la condición que los despla-
zamientos del nodo 3 (U5 y U6 ) son los mismos para cada elemento. Para reiterar ésto, esta hipótesis
es actualmente un requerimiento; puesto que sobre una base fı́sica, la estructura debe permanecer
conectada mediante sus juntas después de la deformación. La compatibilidad de desplazamientos en
los nodos es un requerimiento fundamental del método de elemento finito.
Sustituyendo las Ecuaciones (3.7) y (3.8) en las condiciones nodales de equilibrio establecidas
previamente mediante las Ecuaciones (3.1) a (3.3), obtenemos
U4 F4
k (1) c2 θ1 + k (1) sθ1 cθ1 +
−k (1) c2 θ1 k (1) sθ1 cθ1 −k (2) sθ2 cθ2 −k (2) s2 θ2
U5 F5
(2) 2
k c θ2 k (2) sθ2 cθ2
(1) (1) 2
k sθ1 cθ1 + k s θ1 + U6 F6
(1) 2 (2) 2
−k sθ1 cθ1 −k s θ1 −k sθ2 cθ2 −k s θ2
(1) (2)
(2) (2) 2
k sθ2 cθ2 k s θ2
(3.15)
donde usamos la notación sθ ≡ sin θ, cθ ≡ cos θ, s2 θ ≡ sin2 θ, c2 θ ≡ cos2 θ; para abreviar espacio en
las columnas de la matriz de coeficientes del sistema.
56 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
Las seis ecuaciones algebraicas representadas por la Ecuación (3.15) expresan la serie completa de
condiciones de equilibrio para la estructura cercha de dos-elementos. Esta ecuación puede ser escrita
de manera compacta y sintética como
[ K ] {U } = {F } (3.16)
donde [ K ] es la matriz de rigidez estructural, {U } es el vector de desplazamientos nodales estructural,
y {F } es el vector de cargas nodales aplicadas estructural. Todos estos arreglos serán denominados
globales puesto que están referidos y descritos desde el sistema coordenado global utilizado en el
análisis del problema. Observamos que la matriz de rigidez global es un arreglo cuadrado de dimensión
6×6 correspondiente a los seis desplazamientos globales posibles que posee la estructura. La aplicación
de las condiciones de borde y la solución del sistema de ecuaciones algebráicas que están representadas
por la ecuación matricial obtenida como en el ejemplo se difieren por el momento hasta algo más
adelante, quedando estos temas pendientes de discusión.
o, en forma sintética
[ k (e) ]{u(e) } = {f (e) } (3.17a)
el objetivo presente es transformar estas ecuaciones de equilibrio hacia el sistema coordenado global
en la forma (e) (e)
U F
1(e) 1(e)
U F2
[ K (e) ] 2
(e) = (e) (3.18)
U F
3(e) 3(e)
U4 F4
u(e)
1 = U1(e) cos θ + U2(e) sin θ
u(e)
2 = U3(e) cos θ + U4(e) sin θ
3.3. TRANSFORMACIÓN DE ELEMENTO 57
o, en forma resumida
{u(e) } = [ R(e) ]{U (e) } (3.19a)
donde
cos θ sin θ 0 0
[ R(e) ] = (3.20)
0 0 cos θ sin θ
resumidamente
[ k (e) ][ R(e) ]{U (e) } = {f (e) } (3.21a)
A pesar de haber transformado las ecuaciones de equilibrio desde los desplazamientos locales hacia
los desplazamientos globales como incógnitas, las ecuaciones se expresan todavı́a en el sistema de
coordenadas de elemento. La primera relación de la Ecuación (3.21) desarrollada, es la condición
de equilibrio para el nodo 1 en el sistema coordenado local o de elemento. Si multiplicamos esta
ecuación por cos θ, obtendremos la ecuación de equilibrio para el nodo en la dirección X del sistema
coordenado global. Similarmente, multiplicando por sin θ, se obtiene la ecuación de equilibrio según
la dirección Y global. Aplicando exáctamente el mismo procedimiento con la segunda ecuación de
equilibrio en el nodo 2, se obtiene la ecuación de equilibrio nodal referido al sistema coordenado global.
Las mismas operaciones descritas anteriormente se obtienen si pre–multiplicamos ambos miembros de
la Ecuación 3.21 por [ R(e) ]T , la transpuesta de la matriz de transformación de desplazamientos; esto
es en forma desarrollada,
(e) (e)
U1
cos θ 0
(e) (e)f1 cos θ
ke −ke
(e)
U2 sin θ 0 f1 f1 sin θ
(e) T (e)
[R ] [R ] = = (3.22)
−ke ke U (e) 0 cos θ f2(e) f (e) cos θ
3(e) 2(e)
U4 0 sin θ f2 sin θ
Claramente el término del lado–derecho de la ecuación anterior representa al vector de fuerzas nodales
aplicadas al elemento referido al sistema coordenado global; y en esencia mediante la operación de
pre–multiplicación matricial efectuada en realidad hemos realizado un proceso de transformacion de
fuerzas nodales locales hacia fuerzas nodales globales del elemento. Esta transformación se describe
sintéticamente por la ecuación
{F (e) } = [ R(e) ]T {f (e) } (3.22a)
58 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
Ésta ecuación matricial representa a las condiciones de equilibrio estático nodal del elemento, expresada
en el sistema coordenado global. Comparando este resultado con la Ecuación (3.18), la matriz de rigidez
de elemento en el marco referencial global coordenado se ve que está dada por
y, el vector unitario que direcciona espacialmente al elemento desde el nodo i hasta el nodo j es por
tanto
1
λ̂ = [ (Xj − Xi )Î + (Yj − Yi )Ĵ ] = cos θX Î + cos θY Ĵ (3.27)
L
donde Î y Ĵ son vectores unitarios coincidentes con los ejes X e Y del sistema coordenado global, como
se aprecia en la Figura 3.5.
Recordando la definición del producto escalar de dos vectores y refiriéndonos de nuevo a las Figu-
ras 3.5 y 3.2(a), los valores trigonométricos requeridos para construir la matriz de transformación de
3.4. ENSAMBLE DE LA MATRIZ DE RIGIDEZ GLOBAL 59
Y
j
Yj
Y L
X
Yi
i
J
X
I Xi Xj
elemento también se determinan rápidamente a partir de las coordenadas nodales como los cosenos
directores que aparecen en la Ecuación (3.27).
Xj − Xi
cos θ = cos θX = λ̂◦Î = (3.28a)
L
Yj − Yi
sin θ = cos θY = λ̂◦Ĵ = (3.28b)
L
Ası́, la matriz de rigidez de un elemento barra en coordenadas globales puede ser completamente deter-
minada por la especificación de las coordenadas nodales (ubicación espacial de sus puntos extremos),
el área de su sección transversal, y el módulo de elasticidad del material del elemento.
Las ecuaciones anteriores son las relaciones de conectividad para la cercha, y explı́citamente indican
cómo cada elemento está conectado en la estructura. Por ejemplo la ecuación (3.30a) claramente
muestra que el elemento 1 no está asociado con los desplazamientos globales U3 y U4 (por tanto, no
está conectado al nodo global 2) y, de aquı́, no contribuye a los coeficientes de rigidez que afectan esos
desplazamientos. Esto significa que el elemento 1 no tiene efecto en la tercera y cuarta filas y columnas
de la matriz de rigidez global. Similarmente, el elemento 2 no contribuye en nada a la primera y
segunda filas y columnas del arreglo de coeficientes de rigidez globales.
En lugar de escribir las relaciones individuales de desplazamiento para cada uno de los elementos,
es preferible y conveniente poner todos los datos de desplazamientos globales de todos los elementos
en un arreglo tabular como el mostrado en la Tabla 3.1.
1 1 0
2 2 0
3 0 1
4 0 2
5 3 3
6 4 4
donde la simetrı́a conocida de la matriz de rigidez se ha usado implı́citamente para evitar la repetición.
Se demuestra prontamente que la matriz de rigidez global resultante es idéntica en cada término a
aquella obtenida en la Sección 3.2 mediante ecuaciones de equilibrio. Éste es el método dirécto de
rigidez; la matriz de rigidez global es “ensamblada” por la suma directa de los coeficientes de rigidez
de elemento individuales mediante la tabla de correspondencia de desplazamientos nodales que define
la conectividad del elemento.
Ejemplo 3.1.
Para la cercha plana mostrada en la Figura 3.2, θ1 = π/4, θ2 = 0, y las propiedades de elementos son
tales que k1 = A1 E1 /L1 , k2 = A2 E2 /L2 . Transformar la matriz de rigidez local de cada elemento hacia
el sistema global coordenado, y ensamblar la matriz de rigidez global para la estructura.
> Solución
√
Para el elemento 1, cos θ1 = sin θ1 = 2/2 y c2 θ1 = s2 θ1 = sθ1 cθ1 = 1/2; de modo que sustituyendo
estos valores en la Ecuación (3.25), resulta
1 1 −1 −1
k1
1 1 −1 −1
[ K (1) ] =
2 −1 −1 1 1
−1 −1 1 1
Para el elemento 2, cos θ2 = 1, sin θ2 = 0, lo cual dá la matriz de rigidez transformada
1 0 −1 0
0 0 0 0
[ K (2) ] = k2
−1 0 1
0
0 0 0 0
Ensamblando la matriz de rigidez diréctamente, usando las ecuaciones (3.30a) y (3.30b) se obtiene
La matriz de rigidez global de la estructura es entonces formada con los coeficientes anteriores, resul-
tando el arreglo de coeficientes siguiente
k1 /2 k1 /2 0 0 −k1 /2 −k1 /2
k1 /2
k1 /2 0 0 −k1 /2 −k1 /2
0 0 k2 0 −k2 0
[K ] =
0 0 0 0 0 0
−k1 /2 −k1 /2 −k2 0 k1 /2 + k2 k1 /2
−k1 /2 −k1 /2 0 0 k1 /2 k1 /2
>
62 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
La representación previamente descrita del método dirécto de rigidez es concreta pero embarazosa e
ineficaz en la práctica. El problema principal inherente al método presentado radica en el hecho que cada
término de la matriz de rigidez global se calcula secuencialmente y acompañada de esta construcción
secuencial se requiere que cada elemento sea considerado a cada paso. A continuación describiremos una
técnica que es mucho más eficaz y bien–adecuada a las operaciones que pueden efectuarse mediante una
computadora digital. En este segundo método, la matriz de rigidez de elemento para cada elemento
componente de la estructura es considerada en secuencia, y los coeficientes de la matriz de rigidez
de elemento se agregan a la matriz de rigidez global mediante las relaciones de conectividad nodal
establecidas para un elemento por vez. Ası́, todos los términos de una matriz de rigidez de elemento
individual se agregan a la matriz de rigidez global, después de lo cual no existe necesidad que dicho
elemento sea considerado posteriormente. Para ilustrar esta técnica alternativa, volvemos a escribir las
Ecuaciones (3.29) como
1 2 5 6
(1) (1) (1) (1)
k11 k12 k13 k14 1
(1) (1) (1) (1)
k21 k22 k23 k24 2
[ K (1) ] = (3.31a)
(1) (1) (1) (1)
k31 k32 k33 k34 5
(1) (1) (1) (1)
k41 k42 k43 k44 6
para el elemento 1, y
3 4 5 6
(2) (2) (2) (2)
k11 k12 k13 k14 3
(2) (2) (2) (2)
k21 k22 k23 k24 4
[ K (2) ] = (3.31b)
(2) (2) (2) (2)
k31 k32 k33 k34 5
(2) (2) (2) (2)
k41 k42 k43 k44 6
para el elemento 2. En este bosquejo de las matrices de rigidez para los dos elementos individuales, los
números a la derecha de cada fila y encima de cada columna indican el desplazamiento global asociado
con la fila y columna correspondiente de la matriz de rigidez del elemento. Ası́, combinamos la tabla
de correspondencia de desplazamientos nodales con las matrices de rigidez de elemento individuales.
Para las matrices de elemento, cada coeficiente individual se etiqueta ahora como asociado con una
posición en fila–columna especı́fica de la matriz de rigidez global y puede agregarse directamente a esa
ubicación.
(2)
Por ejemplo, la Ecuación (3.31b) muestra que el coeficiente k24 componente del elemento 2 debe ser
agregado al coeficiente de rigidez global K46 (y vı́a simetrı́a al K64 ). Ası́, podemos tomar cada elemento
en turno y agregar los coeficientes individuales de la matriz de rigidez del mismo a las posiciones
apropiadas en la matriz de rigidez global, aprovechando el proceso de etiquetado o codificación numérica
asignada (con los desplazamientos globales especificados por elemento) a filas y columnas de cada
matriz de rigidez de elemento.
Muchas veces este método es también llamado “método de ensamble dirécto por transporte de
coeficientes”. Y, en la aplicación de esta metodologı́a se procede secuencialmente en orden ascendente:
primero se ensamblan los coeficientes de la matriz de rigidez global del elemento 1 en la matriz de
rigidez global de estructura (que inicialmente se puede suponer que es una matriz nula, con todos sus
coeficientes iniciales igual a cero), luego se repite el proceso para la matriz de rigidez del elemento 2;
y ası́ siguiendo. Si en una instancia intermedia al transportar cierto coeficiente de la matriz de rigidez
de elemento hacia la matriz de rigidez global ocurre que la ubicación en esta última ya está ocupada
por determinado valor, se debe efectuar la suma algebráica con el número pre–existente y reemplazar
el resultado en dicha ubicación.
La forma de las Ecuaciones (3.31) son convenientes sólo para propósitos ilustrativos. Para cálculos
reales, la inclusión de los números de desplazamiento globales dentro la matriz de rigidez de elemen-
to es pesada computacionalmente. Una técnica de algoritmo dinámico conveniente para aplicaciones
3.4. ENSAMBLE DE LA MATRIZ DE RIGIDEZ GLOBAL 63
computacionales se describe posteriormente. Las convenciones siguientes se adoptan para una cercha
2–D modelada por elementos barra lineales:
1. Los nodos globales a los que cada elemento se conecta son denotados por i y j.
2. El origen del sistema coordenado de elemento (sistema coordenado local) se ubica con su origen
en coincidencia con el nodo i, y el eje x de este sistema tiene un sentido positivo en la dirección
del nodo i al nodo j.
3. Los desplazamientos globales en los nodos del elemento son: U2i−1 , U2i , U2j−1 , y U2j ; como fué de-
notado en la Sección 3.2.
Usando estas convenciones, toda la información requerida para definir la conectividad del elemento
y ensamblar la matriz de rigidez global es incluida en una tabla de conectividad nodal de elemento, que
establece un listado de los números de elemento en secuencia y muestra los números globales nodales
i y j a los que cada elemento se conecta. Para la cercha de dos elementos mostrado en la Figura 3.2,
los datos requeridos son como los mostrados en la Tabla 3.2.
1 1 3
2 2 3
Usando los datos nodales de la Tabla 3.2, nosotros definimos, para cada elemento, un vector fila
(de dimensión 1×4) denominado vector de ubicación de desplazamientos, como
{D(e) } = 2i − 1 2i 2j − 1 2j (3.32)
donde cada valor es el número del desplazamiento global que corresponde a las filas y columnas 1, 2,
3, 4 respectivamente, de la matriz de rigidez de elemento. Para la cercha plana de la Figura 3.2, los
vectores de ubicación de desplazamientos son
{D(1) } = 1 2 5 6 (3.33a)
(2)
{D } = 3 4 5 6 (3.33b)
Antes de proceder, notemos la cantidad de información que puede ser obtenida desde una simple
vista de la Tabla 3.2. Con la geometrı́a de la estructura definida, las coordenadas globales (X, Y ) de
cada nodo son especificadas. Usando éstos datos, la longitud de cada elemento y los cosenos directores
de orientación espacial relativa del elemento se calculan mediante las Ecuaciones (3.26) y (3.28),
respectivamente. La especificación del área de sección–transversal A y el módulo de elasticidad E de
cada elemento, permite el cómputo de la matriz de rigidez de elemento en el marco referencial global
usando la Ecuación (3.25). Finalmente, los coeficientes de la matriz de rigidez de elemento yá calculados
se agregan a la matriz de rigidez global mediante el procedimiento de ensamble a través del transporte
de coeficientes usando el vector de ubicación de desplazamientos especificado para el elemento que
está siendo manipulando.
En el contexto del ejemplo actual, el lector debe imaginar un arreglo cuadrado de una serie de
6×6 “buzones” que representan la matriz de rigidez global, cada uno de los cuales está originalmente
vacı́o (es decir, el coeficiente de rigidez inicialmente contenido es cero). Luego, entonces consideramos
64 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
la matriz de rigidez de un elemento individual en el marco referencial global. Usando luego el vector
de ubicación de desplazamientos establecido para el elemento (“las direcciones de correo”), los valores
individuales de la matriz de rigidez del elemento se ponen en el buzón apropiado. De este modo,
cada elemento es procesado en secuencia y sus caracterı́sticas de rigidez agregadas a la matriz global.
Después de que todos los elementos se procesan del modo anteriormente indicado, la serie de buzones
contiene la matriz de rigidez global estructural.
Como fué indicado, la matriz de rigidez global es una matriz singular; por consiguiente, una solu-
ción única a la Ecuación anterior no puede obtenerse directamente. Sin embargo, desarrollando estas
ecuaciones, nosotros no hemos tenido en cuenta todavı́a las restricciones impuestas sobre los despla-
zamientos del sistema, por las condiciones de apoyo que deben existir para evitar el movimiento de
cuerpo rı́gido. En este ejemplo, observamos que las condiciones de borde lı́mite de desplazamiento son
U1 = U2 = U3 = U4 = 0 (3.35)
dejando sólo U5 y U6 para ser determinados. Sustituyendo los valores de las condiciones de borde
anteriormente especificadas y expandiendo la Ecuación (3.34), tendremos formalmente
K15 U5 + K16 U6 = F1
K25 U5 + K26 U6 = F2
K35 U5 + K36 U6 = F3
(3.36)
K45 U5 + K46 U6 = F4
K55 U5 + K56 U6 = F5
K56 U5 + K66 U6 = F6
como las ecuaciones del sistema reducido (ésto es la serie particionada de ecuaciones matriciales,
escrita explı́citamente para los desplazamientos activos). En este ejemplo, F1 , F2 , F3 , y F4 son las
componentes de las fuerzas de reacción de apoyo en los nodos restringidos 1 y 2, mientras F5 y F6
son las componentes globales de fuerza externa aplicada en el nodo 3. Dadas las componentes de
fuerza externa aplicada a la estructura, las dos últimas relaciones del sistema de Ecuaciones (3.36)
pueden resolverse explı́citamente para los desplazamientos U5 y U6 . Los valores obtenidos para estos
dos desplazamientos se sustituyen entonces en las ecuaciones de restricción (las primeras cuatro de las
Ecuaciones (3.36)) y las componentes de fuerza de reacción en los apoyos pueden ser calculadas.
Un procedimiento más general de aplicación de condiciones de borde y evaluación de reacciones
de apoyo es como sigue. Suponiendo que el subı́ndice r denote a los desplazamientos restringidos y
el subı́ndice a denote a los desplazamientos no–restringidos (activos), el sistema de ecuaciones puede
3.6. DEFORMACIONES Y TENSIONES DE ELEMENTO 65
donde los valores de los desplazamientos restringidos {Ur } son conocidos (pero no necesariamente
nulos), como son las fuerzas externas aplicadas {Fa }. Ası́, las incógnitas o desplazamientos activos, se
obtienen mediante la ecuación inferior resultante de la partición como
[Kar ]{Ur } + [Kaa ]{Ua } = {Fa } (3.38a)
−1
{Ua } = [Kaa ] ( {Fa } − [Kar ]{Ur } ) (3.38b)
donde hemos asumido que los desplazamientos especificados {Ur } no son necesariamente cero, aunque
éste normalmente es el caso en una estructura tipo cercha. (De nuevo, notemos que, para la eficacia
numérica se aplican otros métodos distintos a la inversión matricial para obtener las soluciones formal-
mente representadas por la Ecuación (3.38b) para los desplazamientos incógnita). Dada la solución para
los desplazamientos activos, las reacciones se obtienen de la partición superior de la Ecuación (3.37)
como
{Fr } = [Krr ]{Ur } + [Kra ]{Ua } (3.39)
donde [Kra ] = [Kar ]T por la propiedad de simetrı́a de la matriz de rigidez global estructural.
Reemplazando esta última relación en la Ecuación 3.41, adoptando la notación resumida para las
funciones trigonométricas, se obtiene finalmente
(e)
U1
1 U2(e)
= (e) −cθ −sθ cθ sθ
(e)
(3.42)
L
U3(e)
(e)
U4
La tensión normal axial interna luego puede ser hallada apelando a la aplicación de la ley de Hooke
donde E (e) es el módulo de elasticidad de elemento. Ahora, considerando la última ecuación, la tensión
normal axial al interior del elemento en cualquier posición o sección de corte hipotética, en forma
desarrollada resulta tener el aspecto aquı́ mostrado
(e)
U1
E (e) U2(e)
σ = (e) −cθ −sθ cθ sθ
(e)
(3.44)
L
U3(e)
(e)
U4
Del modo en el cual el elemento barra fué formulado aquı́, un valor de tensión axial positivo
indica que el elemento está en tracción y un valor negativo indica compresión, por la convención usual
adoptada para las tensiones normales. Un aspecto interesante de tener en cuenta es que el valor de
tensión normal al interior del elemento es de valor constante, por los términos en la ecuación que
la definen. Nótese además que el cálculo de tensión normal interna indicado en la Ecuación (3.44)
debe ser realizado en una base elemento–por–elemento. Si se desea, las fuerzas del elemento pueden
ser obtenidas mediante la Ecuación (3.21), para poder efectuar una prueba de verificación del cálculo
numérico efectuado ya que para los coeficientes de este vector debe cumplirse la condición básica de
equilibrio estático f1(e) = −f2(e) .
Ejemplo 3.2.
La cercha de dos–elementos mostrada en la Figura 3.6 está sujeta a las cargas externas indicadas.
Usando la misma numeración para nodos y elementos como se aprecia en la Figura 3.2, determine las
componentes del desplazamiento del nodo 3, las componentes de las fuerzas de reacción en los nodos
1 y 2, y los desplazamientos locales de elemento, tensiones internas, y fuerzas nodales locales. Los ele-
mentos tienen módulo de elasticidad E1 = E2 = 107 lb/pulg2 y las áreas de sus secciones transversales
son A1 = A2 = 1,5 pulg2 .
> Solución
Las coordenadas nodales son tales que√ las orientaciones espaciales resultan θ1 = π/4 y θ2 = 0, y
las longitudes de elemento son L1 = 404 + 4062 ≈ 56,57 pulg, L2 = 40 pulg. La rigidez de resorte
equivalente para cada barra en este caso resulta ser
A1 E1 1,5×107
k1 = = = 2,65×105 lb/pulg
L1 56,57
A2 E2 1,5×107
k2 = = = 3,75×105 lb/pulg
L2 40
Como los ángulos de orientación de elemento y la numeración esquemática adoptada son iguales que
en el Ejemplo 3.1, podemos usar el resultado de ese ejemplo para escribir la matriz de rigidez global
3.6. DEFORMACIONES Y TENSIONES DE ELEMENTO 67
300 lb
( 0, 40 ) ( 40, 40 )
2 3 500 lb
2
1
( X, Y ) [ pulg ]
1
( 0, 0 )
como
1,325 1,325 0 0 −1,325 −1,325
1,325
1,325 0 0 −1,325 −1,325
0 0 3,750 0 −3,75 0 ×105 lb/pulg
[K ] =
0 0 0 0 0 0
−1,325 −1,325 −3,75 0 5,075 1,325
−1,325 −1,325 0 0 1,325 1,325
Incorporando las restricciones de desplazamiento U1 = U2 = U3 = U4 = 0 y la solicitación nodal
aplicada a la estructura, las ecuaciones globales de equilibrio resultan ser
1,325 1,325 0 0 −1,325 −1,325
0 F1
1,325 1,325 0 0 −1,325 −1,325
0
F2
5
0 0 3,75 0 −3,75 0
0
F 3
10 =
0 0 0 0 0 0 0 F4
−1,325 −1,325 −3,75 0 5,075 1,325 U 500
5
−1,325 −1,325 0 0 1,325 1,325 U 6 300
en las que las lı́neas divisorias al interior de los arreglos matriciales representan el proceso de particio-
namiento realizado, con el mismo esquema de la Ecuación (3.37). De aquı́, los desplazamientos activos
estan gobernados por la ecuación
5 5,075 1,325 U5 500
10 =
1,325 1,325 U6 300
La solución simultánea de ambas ecuaciones proporciona los siguientes valores para las componentes
de desplazamiento del nodo 3
U5 = 5,333×10−4 pulg U6 = 1,731×10−3 pulg
Como todos los desplazamientos restringidos tienen valor nulo, de la misma ecuación gobernante del
comportamiento del sistema, mostrada en forma particionada, se puede obtener la siguiente ecuación
para el cálculo de las fuerzas de reacción de apoyo
F1
−1,325 −1,325
−300
F2 5 −1,325 −1,325 0,5333 −300
−3
= 10 10 = lb
F3 −3,75 0 1,731 −200
F4 0 0 0
68 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
Podemos notar que una vez calculadas las reacciones de apoyo, es posible efectuar una verifica-
ción de las condiciones de equilibrio estático completo de la estructura. Si usted realiza el cálculo,
podrá comprobar que la suma de fuerzas según las direcciones globales X e Y son simultáneamente
nulas y la suma de momentos respecto al eje global Z es idénticamente cero.
Para el elemento 1, los desplazamientos nodales locales, o sea referidos al sistema coordenado propio
del mismo elemento, pueden calcularse aplicando la ecuación de transformación de desplazamientos,
{u(e) } = [ R(e) ]{U (e) }, lo que dá como resultado
U1 √ 0
(1) 2 1 1 0 0
u1 U2 0 0
−3
= [ R (1)
] = 10 = 10−3 pulg
u(1)
2
U5 2 0 0 1 1 0,5333 1,6
U6 1,731
La tensión normal axial interna de elemento puede calcularse ahora usando la Ecuación (3.44), que
aplicada al elemento 1 dá como resultado
U1
E (1) U2
σ = (1) −cθ1 −sθ1 cθ1 sθ1
(1)
L
U5
U6
√ 0
107
2 0 −3 2
= −1 −1 1 1 10 = 283,03 lb/pulg
56,57 2
0,5333
1,731
donde el signo positivo de esta tensión normal interna nos indica que el elemento está en estado de
solicitación de tracción.
Las fuerzas nodales locales de elemento, es decir aquellas fuerzas que actúan con dirección axial
en los puntos extremos del elemento, se calculan aplicando la Ecuación (3.17). Para el elemento 1
tendrı́amos
(1)
f1 1 −1 u(1)1 5 1 −1 0 −3 −424
= k 1 = 2,65×10 10 = lb
f2(1) −1 1 u(1)
2 −1 1 1,6 424
y los signos algebráicos de la solución obtenida para las fuerzas nodales axiales de elemento también
indican solicitación de tracción.
Para el elemento 2, la misma secuencia de procedimiento proporciona
(2)
U3
0
u1 U4 1 0 0 0 0 0
−3
(2)
= [R ] = 10 = 10−4 pulg
u(2)
2
U5
0 0 1 0 0,5333 0,5333
U6 1,731
U3
E (2) U4
σ (2) = −cθ 2 −sθ 2 cθ 2 sθ 2
L(2) U5
U6
0
107
0 −3 2
= −1 0 1 0 10 = 133,32 lb/pulg
40
0,5333
1,731
(2)
f1 1 −1 u(2) 1 5 1 −1 0 −3 −200
= k2 = 3,75×10 10 = lb
f2(2) −1 1 u(2)
2 −1 1 0,5333 200
que también indican que este elemento está sometido a tracción. >
3.7. UN EJEMPLO COMPLETO 69
6000 lb
U4
40 pulg 40 pulg U8 U12
2 3 4 7 6
4000 lb U3 U7 U11
40 pulg 2 4 6 8
Y
U2
U6 U10
1 5
2000 lb 1 U1 3 U5 5 U9 X
2000 lb
Se piensa que el método finito es un procedimiento de propósito general para analizar problemas
en los cuales la solución general no es conocida; sin embargo, es informativo en los ejemplos de este
capı́tulo (desde que el elemento barra posee una formulación exacta) verificar las soluciones en términos
de la tensión axial simplemente calculada como σ = F/A para un miembro axialmente cargado. Se
invita al lector a calcular la tensión axial por la fórmula simple de tensión normal para cada uno de los
ejemplos a fin de verificar que las soluciones basadas en el método de elemento finito son completamente
correctas.
Paso 1. Especificar el sistema global coordenado, asignar números a los nodos, especificar los elemen-
tos, establecer la conectividad de los diversos elementos, y definir las componentes de desplaza-
mientos nodales globales, como se muestra en la Figura 3.7(b).
Paso 2. Calcular los valores del coeficiente de rigidez de resorte equivalente k = E A/L para todos
los elementos
1,5(107 )
k (1) = k (3) = k (4) = k (5) = k (7) = k (8) = = 3,75×105 lb/pulg.
40
1,5(107 )
k (1) = k (3) = √ = 2,65×105 lb/pulg.
40 2
70 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
Paso 3. Transformar las matrices de rigidez locales de elemento hacia el sistema global coordenado.
Utilizando la Ecuación (3.25) con los valores de los siguientes ángulos de orientación espacial
tenemos
1 0 −1 0
5 0 0 0 0
(1) (3) (5) (7)
[ K ] = [ K ] = [ K ] = [ K ] = 3,75×10
−1 0 1 0
0 0 0 0
0 0 0 0 1 1 −1 −1
(4) (8) 5 0
1 0 −1 2,65×105 1 1 −1 −1
[ K ] = [ K ] = 3,75×10 [ K (2) ] =
−1
0 0 0 0 2 −1 1 1
0 −1 0 1 −1 −1 1 1
1 −1 −1 1
2,65×105
−1 1 1 −1
[ K (6 ] = −1 1
2 1 −1
1 −1 −1 1
Paso 4a. Construir la tabla de correspondencia de desplazamientos locales – hacia – globales. Con
referencia a la Figura 3.7(b), las relaciones de desplazamientos y conectividad se muestran en la
Tabla 3.3
1 1 1 0 0 0 0 0 0
2 2 2 0 0 0 0 0 0
3 0 0 1 0 0 0 0 0
4 0 0 2 0 0 0 0 0
5 3 0 0 1 1 0 0 0
6 4 0 0 2 2 0 0 0
7 0 3 3 3 0 3 1 0
8 0 4 4 4 0 4 2 0
9 0 0 0 0 3 1 0 1
10 0 0 0 0 4 2 0 2
11 0 0 0 0 0 0 3 3
12 0 0 0 0 0 0 4 4
Paso 4b. Alternativa y más eficazmente, formamos la tabla de conectividad de nodos de elemento
(Tabla 3.4), y los correspondientes vectores de ubicación de desplazamientos globales de elemento
para cada miembro componente de la estructura, lo cual es mostrado en la serie siguiente
{D(1) } = 1 2 5 6 {D(2) } = 1 2 7 8 {D(3) } = 3 4 7 8
{D(4) } = 5 6 7 8 {D(5) } = 5 6 9 10 {D(6) } = 9 10 7 8
{D(7) } = 7 8 11 12 {D(8) } = 9 10 11 12
3.7. UN EJEMPLO COMPLETO 71
1 1 3
2 1 4
3 2 4
4 3 4
5 3 5
6 5 4
7 4 6
8 5 6
Paso 5. Ensamblar la matriz de rigidez global ya sea por el Paso 4a o aplicando el Paso 4b. Los
coeficientes de la matriz de rigidez global estructural son:
(1) (2) (1) (2)
K11 =k11 +k11 =(3,75+2,65/2)105 K12 =k12 +k12 =(0+2,65/2)105
(1)
K13 =K14 =0 K15 =k13 =−3,75(105 )
(1) (2)
K16 =k14 =0 K17 =k13 =−(2,65/2)105
(2)
K18 =k14 =−(2,65/2)105 K19 =K1,10 =K1,11 =K1,12 =0
(1) (2)
K22 =k22 +k22 =0+(2,65/2)105 K23 =K24 =0
(1) (1)
K25 =k23 =0 K26 =k24 =0
(2) (2)
K27 =k23 =−(2,65/2)105 K28 =k24 =−(2,65/2)105
(1)
K29 =K2,10 =K2,11 =K2,12 =0 K33 =k33 =(3,75/2)105
(3)
K34 =k12 =0 K35 =K36 =0
(3) (3)
K37 =k13 =−3,75(105 ) K38 =k14 =0
(3)
K39 =K3,10 =K3,11 =K3,12 =0 K44 =k22 =0
(3)
K45 =K46 =0 K47 =k23 =0
(3)
K48 =k24 =0 K49 =K4,10 =K4,11 =K4,12 =0
(1) (4) (5) (1) (4) (5)
K55 =k33 +k11 +k11 =(3,75+0+3,75)105 K56 =k34 +k12 +k12 =0
(4) (4)
K57 =k13 =0 K58 =k14 =0
(5) (5)
K59 =k13 =−3,75(1−5 ) K5,10 =k14 =0
(2) (4) (5)
K5,11 =K5,12 =0 K66 =k44 +k22 +k22 =(0+3,75+0)105
(4) (4)
K67 =k23 =0 K68 =k24 =−3,75(105 )
(5) (5)
K69 =k23 =0 K6,10 =k24 =0
(2) (3) (4) (6) (7) (2) (3) (4) (6) (7)
K77 =k33 +k33 +k33 +k33 +k11 K78 =k34 +k34 +k34 +k34 +k12
=(2,65/2+3,75+0+2,65/2+3,75)105 =(2,65/2+0+0−2,65/2+0)105
(6)
K6,11 =K6,12 =0 K79 =k13 =−(2,65/2)105
(6) (7)
K7,10 =k23 =(2,65/2)105 K7,11 =k13 =(−3,75(105 )
(7) (2) (3) (4) (6) (7)
K7,12 =k14 =0 K88 =k44 +k44 +k44 +k44 +k22
(6)
K89 =k14 =(2,65/2)105 =(2,65/2+0+3,75+2,65/2+0)105
(6) (7)
K8,10 =k24 =−(2,65/2)105 K8,11 =k23 =0
(7) (5) (6) (8)
K8,12 =k24 =0 K99 =k33 +k11 +k11 =(3,75+2,65/2+0)105
72 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
Paso 6. Aplicar las restricciones como son dictadas por las condiciones de borde. En este ejemplo,
los nodos 1 y 2 están conectados a apoyos articulados fijos, de modo que las restricciones de
desplazamiento en estos puntos se traducen como
U1 = U2 = U3 = U4 = 0
Por tanto, las primeras cuatro ecuaciones en la ecuación matricial gobernante pueden desecharse.
Esto último es equivalente a eliminar las cuatro primeras filas y las cuatro primeras columnas de
la matrix 12×12 de coeficientes del sistema
[ K ]{U } = {F }
la eliminación de cierto número de ecuaciones es posible en este caso, por la consideración que
en este ejemplo particular los desplazamientos restringidos son todos nulos. Si es que no lo fue-
ran, las restricciones son consideradas como en la Ecuaciones (3.38b) y (3.39), en cuyo caso las
restriciones no–nulas imponen fuerzas adicionales sobre los desplazamientos no–restringidos. Las
fuerzas restringidas (reacciones de apoyo) no pueden ser obtenidas hasta que los desplazamien-
tos activos no hayan sido calculados. Pero, en el problema que estamos analizando si podemos
efectuar la reducción del problema por imposición de las condiciones de borde en los apoyos. Efec-
tuando la eliminación que fué mencionada anteriormente, las ecuaciones globales que gobiernan el
comportamiento resultan
U5
0
U6 −2000
U 0
7
U8 0
[ Kaa ] =
U9
2000
U10 0
U11 4000
U12 6000
como relaciones que describen el comportamiento del sistema estructural reducido solamente a los
desplazamientos activos.
Paso 7. Resolver las ecuaciones que corresponden a los desplazamientos no–restringidos. Para el ejem-
plo actual, las ecuaciones se resuelven usando un programa de hoja de cálculo, invirtiendo la
(relativamente pequeña) matriz global de rigidez para obtener
U5
0,02133
U6 0,04085
U −0,01600
7
U8 0,04619
= pulg
U9
0,04267
U10 0,15014
U11 −0,00533
U12 0,16614
Paso 8. Aplicar un proceso de substitución – regresiva de los datos solución de los desplazamientos
para poder calcular las fuerzas de reacción de apoyo. Utilizando la Ecuación (3.39), con la condición
3.7. UN EJEMPLO COMPLETO 73
{Ur } = {0}, usamos las cuatro ecuaciones previamente ignoradas para evaluar las componentes de
fuerza actuantes en los nodos 1 y 2 (conectados a los apoyos de la estructura). Las ecuaciones de
restricción son de la forma
y, substituyendo los valores de desplazamiento previamente hallados, se obtiene como solución para
las fuerzas de reacción de apoyo
F1
−12000
F2 −4000
= lb
F3 6000
F4 0
Se deja a la inquietud del lector la utilización de éstas componentes de fuerzas de reacción para
efectuar una verificación del equilibrio estático global de la estructura.
Paso 9. Calcular la deformación y la tensión en cada elemento. La tarea mayor de evaluación de
valores efectuada en el Paso 7 proporciona las componentes de desplazamiento de cada nodo en el
sistema coordenado global. Con esta información y los desplazamientos restringidos conocidos, los
desplazamientos de cada uno de los elementos en su propio sistema de coordenadas puede obtenerse
(aplicando las ecuaciones de transformación de desplazamientos); de aquı́, la deformación y la
tensión normal interna en cada elemento puede calcularse.
Por ejemplo, para el elemento 2 tenemos
u(2)
1 = U1 cos θ2 + U2 sin θ2 = 0
√
u(2)
2 = U7 cos θ2 + U8 sin θ2 = (−0,01600 + 0,04618) 2/2
= 0,02134 pulg
Sólo se presentan los resultados para elemento 2 como un ejemplo. En cualquier software de ele-
mento finito, los resultados para cada elemento están disponibles y pueden ser examinados por el
usuario del programa, como él desee (ésto es parte del post–procesamiento).
Se muestran resultados para cada uno de los ocho elementos en la Tabla 3.5; y por la convención
usual de signos para tensiones, los valores positivos indican tensión de tracción mientras los valores
negativos corresponden a tensión de compresión. Para obtener los resultados presentados para este
ejemplo, nosotros usamos el programa de hoja de cálculo MathCad para invertir la matriz de rigidez
y poder calcular los desplazamientos activos en un proceso semi–manual; también utilizamos un
paquete de software de elemento finito muy popular (Abaqus), y comprobamos que las soluciones
resultantes de cada procedimiento fueron idénticas.
>
El ejemplo que resolvimos muestra los pasos de procedimiento secuenciales a ser aplicados para
analizar una estructura cercha plana. Este procedimiento es general y puede aplicarse en el mismo
orden con el que aquı́ fué aplicado. Seguramente usted percibirá que existe ligera diferencia en el
procedimiento cuando se utilice un programa computacional para resolver éste problema. Sin embargo,
el tener conocimiento de las operaciones necesarias para hallar variables de inicio desconocidas es lo
que vale en la solución manual que presentamos en este ejemplo.
74 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
1 5,33 5333
2 3,77 3771
3 -4,0 -4000
4 1,33 1333
5 5,33 5333
6 -5,67 -5657
7 2,67 2667
8 4,00 4000
U3j-1
(e)
j
U3j U3j-2
U3i-1 y
Y U3i i x
z U3i-2
Z X
Las cerchas tridimensionales (3–D) también pueden modelarse usando el elemento barra, previnien-
do que las conexiones entre los elementos son tales que sólo se transmite carga axial. Estrictamente
hablando, esto requiere que todas las conexiones sean juntas de rótula esférica. Incluso cuando la
restricción de conexión no se satisface precisamente, el análisis de una cercha 3–D usando elementos
barra es a menudo de mucho valor porque permite obtener estimaciones preliminares de las tensiones
en los miembros componentes de la estructura, que en el contexto del diseño mecánico es valioso en la
determinación de las propiedades estructurales requeridas. Refiriéndonos a la Figura 3.8 que muestra
un elemento barra uni–dimensional conectado a los nodos i y j en un marco de referencia global tri–
dimensional, el vector unitario a lo largo del eje del elemento (es decir, el eje de referencia del mismo)
expresado en el sistema global es
(e) 1
λ̂ = [ (Xj − Xi )Î + (Yj − Yi )Ĵ + (Zi )K̂ ] (3.45)
L
o bien (e)
λ̂ = cos θX Î + cos θY Ĵ + cos θZ K̂ (3.45a)
3.8. CERCHAS TRI–DIMENSIONALES 75
Ası́, los desplazamientos del elemento se expresan en componentes en el sistema global coordenado
3–D como
u(e)
1 = U1(e) cos θx + U2(e) cos θy + U3(e) cos θz (3.46a)
(e) (e) (e) (e)
u 2 =U4 cos θx + U 5 cos θy + U 6 cos θz (3.46b)
Aquı́, usamos la anotación que los desplazamientos 1 y 4 de elemento están en la dirección X global,
los desplazamientos 2 y 5 se producen según la dirección Y global, y los desplazamientos de elemento
3 y 6 están en la dirección Z global.
En analogı́a con la Ecuación (3.19), las Ecuaciones (3.46) pueden ser expresadas matricialmente
como (e)
U1
(e)
U
(e) 2
(e)
u1 cos θx cos θy cos θz 0 0 0 U3
= (3.47)
u(e)
2 0 0 0 cos θx cos θy cos θz U4(e)
U (e)
5(e)
U6
o en forma resumida,
u(e) = [ R(e) ]{U (e) } (3.47a)
(e)
donde [ R ] es la matriz de transformación de desplazamientos, que mapea los desplazamientos de
dirección axial uni–dimensional de los nodos extremos del elemento hacia el sistema coordenado global
tri–dimensional.
Siguiendo un procedimiento idéntico al usado para el caso bi–dimensional en la Sección 3.3, la matriz
de rigidez de elemento en el sistema de coordenadadas locales (propio del elemento) se transforma hacia
las coordenadas globales 3–D mediante
El ensamble de la matriz de rigidez global (es decir, las ecuaciones de equilibrio estático), es idéntico
al procedimiento discutido para el caso bidimensional con la obvia excepción que serán considerados
tres desplazamientos para cada nodo.
Ejemplo 3.4.
La cercha de tres–miembos mostrada en la Figura 3.9(a) es conectada por juntas de rótula y abrazadera
esféricas, estando fijos los nodos 1, 2, y 3. Una fuerza de 5000 lb está aplicada en el nodo 4 en la dirección
Y negativa, como es mostrado. Cada uno de los tres miembros es idéntico y exhibe caracterı́stica de
una rigidez axial de 3×105 lb/pulg. Calcular las componentes del desplazamiento del nodo 4 usando
un modelo de elemento finito con elementos barra como fué descrito.
> Solución
Primero, debemos notar que la cercha 3–D con cuatro nodos tiene 12 posibles desplazamientos. Sin
76 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
(0, 0, -30) U5
Y 2 U6 U4
Y 2
Z X 2
4 Z X
(40, 0, 0) U11
U2
1
1 (0, 0, 30) U3 4 U10
5000 lb 1 U12
U1
(0, -30, 0) 3 U8 3
(x, y, z) [pulg] 3
U9 U7
(a) Disposición geométrica (b) Codificación estructural
embargo, desde que los nodos 1, 2, y 3 son fijos, se conocen nueve de los posibles desplazamientos
los que tienen valor nulo. Por consiguiente, necesitamos ensamblar sólo una porción de la matriz de
rigidez del sistema para resolver el problema y hallar los tres desplazamientos desconocidos. Utilizando
el esquema de numeración mostrado en Figura 3.9(b) y la tabla de correspondencia de desplazamientos
locales–hacia–globales (Tabla 3.6), necesitamos sólo considerar las ecuaciones
K10,10 K10,11 K10,12 U10 0
K11,10 K11,11 K11,12 U11 = −5000
K12,10 K12,11 K12,12 U12 0
1 1 0 0
2 2 0 0
3 3 0 0
4 0 1 0
5 0 2 0
6 0 3 0
7 0 0 1
8 0 0 2
9 0 0 3
10 4 4 4
11 5 5 5
12 6 6 6
Antes de ensamblar los coeficientes requeridos en la matriz de rigidez del sistema, debemos trans-
formar las matrices locales de rigidez individuales hacia el sistema de coordenadas globales como sigue.
Elemento 1.
(1) 1
λ̂ = [ (40 − 0)Î + (0 − 0)Ĵ + (0 − 30)K̂ ] = 0,8Î − 0,6K̂
50
3.8. CERCHAS TRI–DIMENSIONALES 77
aplicando la Ecuación (3.45); y por tanto cx = 0,8, cy = 0, y cz = −0,6. Ahora, aplicando la Ecua-
ción (3.49) resulta
0,64 0 −0,48 −0,64 0 0,48
0 0 0 0 0 0
5 −0,48 0 0,36 0,48 0 −0,36
(1)
[ K ] = 3×10 lb/pulg
−0,64 0 0,48 0,64 0 −0,48
0 0 0 0 0 0
0,48 0 −0,36 −0,48 0 0,36
Elemento 2.
(2) 1
λ̂ [ (40 − 0)Î + (0 − 0)Ĵ + (0 − (−30))K̂ ] = 0,8Î + 0,6K̂
=
50
0,64 0 0,48 −0,64 0 −0,48
0 0 0 0 0 0
(2) 5
0,48 0 0,36 −0,48 0 −0,36
[ K ] = 3×10 lb/pulg
−0,64 0 −0,48 0,64 0 0,48
0 0 0 0 0 0
−0,48 0 −0,36 0,48 0 0,36
Elemento 3.
1
(3)
λ̂ [ (40 − 0)Î + (0 − (−30))Ĵ + (0 − 0)K̂ ] = 0,8Î + 0,6Ĵ
=
50
0,64 0,48 0 −0,64 −0,48 0
0,48
0,36 0 −0,48 −0,36 0
(3) 5
0 0 0 0 0 0
[ K ] = 3×10 lb/pulg
−0,64 −0,48 0 0,64 0,48 0
−0,48 −0,36 0 0,48 0,36 0
0 0 0 0 0 0
Refiriéndonos a las últimas tres filas de la tabla de correspondencia de desplazamientos, los coefi-
cientes requeridos de la matriz de rigidez global se ensamblan como sigue
(1)
K10,10 = k44 (2)
+ k44 (3)
+ k44 = 3×105 (0,64 + 0,64 + 0,64) = 5,76×105 lb/pulg
(1)
K10,11 = K11,10 = k45 (2)
+ k45 (3)
+ k45 = 3×105 (0 + 0 + 0,48) = 1,44×105 lb/pulg
(1)
K10,12 = K12,10 = k46 (2)
+ k46 (3)
+ k46 = 3×105 (−0,48 + 0,48 + 0) = 0 lb/pulg
(1)
K11,11 = k55 (2)
+ k55 (3)
+ k55 = 3×105 (0 + 0 + 0,36) = 1,08×105 lb/pulg
(1)
K11,12 = K12,11 = k56 (2)
+ k56 (3)
+ k56 = 3×105 (0 + 0 + 0) = 0 lb/pulg
(1)
K12,12 = k66 (2)
+ k66 (3)
+ k66 = 3×105 (0,36 + 0,36 + 0) = 2,16×105 lb/pulg
El sistema de ecuaciones a ser resuelto para los desplazamientos del nodo 4 es entonces
5,76 1,44 0 U10 0
105 × 1,44 1,08 0 U11 = −5000
0 0 2,16 U12 0
requeridos para obtener los resultados de los elementos individuales son prontamente obtenidos por
las operaciones matriciales descritas previamente. Una vez que los desplazamientos globales activos
han sido calculados, las restantes variables secundarias (ası́ llamadas) como la deformación, la tensión
interna, y las fuerzas axiales de elemento, se pueden calcular fácilmente usando las matrices y funciones
de interpolación de desplazamiento desarrolladas en la formulación del problema original.
3.9. Resumen
Este capı́tulo desarrolla el procedimiento completo para realizar un análisis de elemento finito de
una estructura y lo ilustra mediante varios ejemplos. Aunque sólo el elemento axial simple se ha usado,
el procedimiento descrito es común al método de elemento finito para todos los tipos de análisis y
elementos utilizados en ellos, como se pondrá en claro en los capı́tulos subsecuentes. El método directo
de rigidez es de lejos la técnica más expedita para ensamblar las matrices de sistema requeridas para
el análisis de elemento finito, y los procedimientos de cálculo secuencial aquı́ presentados también son
muy adecuados para generar los algoritmos de programación que efectúen la valoración automática de
variables de interés durante la implementación del método en una computadora digital.
Problemas propuestos
3.1. En la cercha de dos elementos mostrada en la Figura 3.2; considere los valores siguientes de
orientación espacial: θ1 = 45◦ , θ2 = 15◦ , y de solicitación F5 = 5000 lb, F6 = 3000 lb.
a. Usando sólamente las ecuaciones de equilibrio estático de fuerzas, hallar las fuerzas que actúan
en cada elemento, como también las fuerzas de reacción de apoyo.
b. Asumiendo que cada miembro tiene coeficiente de rigidez axial de resorte equivalente igual a
ke = 52000 lb/pulg, hallar la deformación axial en cada miembro.
c. Utilizando la solución obtenida en el inciso b, calcular las componentes de desplazamiento del
nodo 3.
3.2. Calcular los desplazamientos X e Y del nodo 3 usando procedimientos de elemento finito y los
datos proporcionados en el Problema 3.1. También calcular la fuerza en cada elemento. Que
grado de exactitud se alcanza en la solución de éste problema comparado con los resultados
obtenidos en el Problema 3.1 anterior.
3.3. Verificar la exactitud de la Ecuación (3.25), efectuando la multiplicación directa de las matrices
involucradas en la relación que define este arreglo matricial.
3.4. Mostrar que la matriz de rigidez transformada hacia el sistema coordenado global del elemento
barra tı́pico perteneciente a una estructura cercha bi–dimensional, dada por la Ecuación (3.25),
es una matriz singular.
3.5. Cada uno de los elementos barra mostrados en la Figura P3.5 tiene sección transversal circular
sólida de diámetro d = 0,5 pulg. El material de composición de estos miembros es acero de bajo
contenido de carbón, con módulo de elasticidad lineal E = 3×107 psi. Las coordenadas nodales
indicadas están dadas referidas a un sistema global coordenado X–Y (en pulg). Determinar la
matriz de rigidez global de cada uno de estos elementos.
Problemas propuestos 79
(30, 30) 2
(5, 30)
2
Y (30, 15)
1 1
3.6. Repetir el Problema 3.5 para los elementos barra mostrados en la Figura P3.6. Para estos ele-
mentos d = 40 mm, E = 69 GPa, y las coordenadas nodales indicadas están en metros.
1 (-0.3, 3)
(0, 0) 2
2
(0.4, 0.2)
1
(0.1, 0.1)
1
2
(0.2, -0.2) (1, 2)
Figura P3.6
3.7. Para cada una de las estructuras mostradas en la Figura P3.7, construir una tabla de correspon-
dencia de desplazamientos locales–hacia–globales, como aquella mostrada en la Tabla (3.1).
3 6 5 10 7 14 9 18 11
2 3 5 7 9 11 13 15 17 19 21
1 12
1 2 4 4 8 6 12 8 16 10 20
(a)
4 7 6
6 10
8
4
3 5 5 9 7
3 6 6 12 8
3
2
5 8
4 11
2
2 5 13
1
3 10
2 7
1 1 4 9 7
(b)
3 (c)
80 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
6
9 10
5 8
6 12
3 9
7 11
2 5 17
3 13 15
1 11
1 2 4 4 8 7 14 10 16
(d)
Figura P3.7
2 2 3 6 5
1 3 5 7
1
4 4
(e)
Figura P3.7
3.8. Para cada una de las cerchas de la Figura P3.7, expresar los datos de conectividad para cada
elemento en la forma de la Ecuación (3.32)
3.9. Para cada uno de los elementos mostrados en la Figura P3.9, los desplazamientos globales han
sido calculados como U1 = 0,05 in, U2 = 0,02 in, U3 = 0,075 in, U4 = 0,09 in. El área transversal
de todos los elementos es A = 0,75 pulg2 y el módulo de elasticidad es E = 107 lb/pulg2 . Usando
ecuaciones de elemento finito calcule:
a. El desplazamiento axial en cada nodo ubicado en sus extremos.
b. La deformación neta producida en cada elemento.
c. La tensión normal axial interna.
d. Las fuerzas axiales nodales actuantes en los extremos.
Los valores de tensión normal interna calculados, coinciden con el valor simple σ (e) = f (e) /A(e) ?.
U4
U4 U3
U3 U4
U3 110 °
U2
45 ° U2 U2
30 °
U1 U1 U1
(a) (b) (c)
Figura P3.9
Problemas propuestos 81
3.10. La cercha plana mostrada en la Figura P3.10 está sujeta a una carga vertical descendente en
el nodo 2. Determine mediante el método dirécto de rigidez la deflexión del nodo 2 referido al
sistema coordenado global especificado, y la tensión axial interna en cada uno de los elementos.
Considere que para ambos elementos A = 0,5 pulg y E = 3×106 lb/pulg2 .
(0, 0) X (40, 0)
1 3
(30, -10)
2
1500 lb
Figura P3.10
3.11. La cercha plana mostrada en la Figura P3.11 está compuesta de miembros que tienen sección
transversal cuadrada de 15 mm × 15 mm, y módulo de elasticidad E = 69 GPa.
a. Ensamblar la matriz de rigidez global.
b. Calcular los desplazamientos nodales en el sistema coordenado global para las cargas aplicadas
mostradas.
c. Evaluar la tensión axial interna en cada uno de los elementos.
3 kN
3
5 kN
2 4
1.5 m
1.5 m
Figura P3.11
3.12. Repetir el Problema 3.11 asumiendo que los elementos 1 y 4 son removidos. En que porcentaje
varı́an las tensiones internas en los elementos que permanecen soportando la carga aplicada ?.
3.13. La cercha en voladizo mostrada en la Figura P3.13 fué construı́da para soportar un guinche
con sistema de cables (no mostrado) para levantar y descargar materiales de construcción en
un edificio en construcción. Los miembros componentes de esta cercha son listones de madera
pino, con dimensiones de sección transversal rectangular de 1.75 pulg × 3.5 pulg; y módulo de
elasticidad E = 2×106 lb/pulg2 . Usando el método dirécto de rigidez, calcular
a. Las componentes de desplazamiento global de todos los nodos no–restringidos.
b. Las tensiones axiales internas para todos los miembros componentes de la estructura.
82 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
2 5
3
5000 lb
Nodo X Y
Y 1 0 0
2 0 96
45°
3 96 96
X 4 96 151.4
1
5 192 96
(pulg)
Figura P3.13
3.14. La Figura P3.14 muestra una cercha plana de dos–miembros soportada por un resorte lineal
elástico. Los miembros componentes de la cercha son de sección transversal circular sólida te-
niendo diámetro d = 20 mm y módulo de elasticidad E = 80 GPa. El resorte lineal tiene
coeficiente de rigidez constante de k = 5 N/mm.
a. Ensamblar la matriz de rigidez global y calcular los desplazamientos globales de los nodos
no–restringidos.
b. Calcular las fuerzas de reacción de apoyo y verificar las condiciones de equilibrio estático para
la estructura.
c. Verificar el balance de energı́a del sistema. Está la energı́a de deformación interna acumulada
balanceada con el trabajo mecánico de las fuerzas aplicadas ?.
3m
15 kN
50°
4m
k
Figura P3.14
3.15. Repetir el Problema 3.14 si el resorte es removido, y evaluar la máxima variación porcentual
relativa en las tensiones internas de los elementos componentes de la estructura.
Problemas propuestos 83
3.16. Debido a una conexión de apoyo defectuosa, el nodo 1 en el Problema 3.13 se desplaza horizon-
talmente 0,5 mm a la izquierda cuando la carga es aplicada. Repita los cálculos efectuados para
esta nueva condición. La solución cambia, o no lo hace ?. Por qué razones?.
3.17. Dado el siguiente sistema de ecuaciones lineales simultáneas, representado por la ecuación ma-
tricial
10 −10 0 0 x1 F1
−10 20
−10 0 x F2
2
=
0 −10 20 −10
x3 F3
0 0 −10 10 x4 F4
x1 = 0 x3 = 1,5 F2 = 20 F4 = 35
y las condición especificada U2 = 0,5; determinar el valor de todas las incógnitas sobrantes apli-
cando el procedimiento sugerido en el problema anterior.
3.19. Para la cercha mostrada en la Figura P3.19, resolver para los desplazamientos globales en el
nodo 3 y las tensiones internas en cada elemento. Los miembros componentes de la estructura
tienen área de sección transversal A = 1,0 pulg2 y módulo de elasticidad E = 1,5×107 lb/pulg2 .
2400 lb
72 pulg 60°
3 1800 lb
1 60° 2
30°
Figura P3.19
3.20. Cada elemento barra mostrado en la Figura P3.20 es parte de una cercha tri–dimensional. Las
cooordenadas nodales (en pulg) están especificadas en un sistema global coordenado X–Y –Z.
Considerando un área de sección transversal A = 2 pulg2 y módulo de elasticidad E = 3×107
lb/pulg2 para los miembros componentes; calcular la matriz de rigidez global de cada uno de los
84 CAPÍTULO 3. EL MÉTODO DIRÉCTO DE RIGIDEZ
elementos.
2 (0, 80, 0) 1 (0, 0, 0)
2 (30, 30, -20)
(50, 40, 30) 2
1 (0, 0, 0)
1 (10, 10, 0) (0, 0, 0)
1 2 (0, -80, 0)
Figura P3.20
3.21. Verificar la Ecuación (3.49) mediante cálculo dirécto del producto matricial indicado.
3.22. Mostrar que la tensión axial interna en un elemento barra de una cercha tri–dimensional es dada
por
h i u(e) E
dN2 (x)
σ = E = E dN
1
dx) dx
1
(e) = −1 1 [ R(e) ]{U (e) }
u2 L
notar que esta relación es la misma que para un elemento barra de una cercha bi–dimensional.
3.23. Determinar las tensiones normales axiales y las fuerzas nodales locales para cada elemento
mostrado en la Figura P3.20, asumiendo que el nodo 1 está fijo y el nodo 2 tiene como despla-
zamientos globales U4 = U5 = U6 = 0,06 pulg.
3.24. Utilizar las Ecuaciones (3.46) para expresar la energı́a de deformación elástica acumulada en
el interior de un elemento barra en términos de los desplazamientos globales. Aplicar el pri-
mer teorema de Castigliano y mostrar que la matriz de rigidez global es idéntica a aquella que
está determinada por la Ecuación (3.49).
3.26. Ensamblar la matriz de rigidez global de la cercha tri–dimensional mostrada en la Figura P3.26
y evaluar las componentes de desplazamiento en el nodo 4. También, calcular la tensión axial
interna en cada elemento. Tomar como valores numéricos para el área de sección transversal
A = 1,5 pulg2 y el módulo de elasticidad E = 107 lb/pulg2 , para todos los elementos.
Y
1
3
X
2 4 Nodo X Y Z
Z FY=-1500 lb
1 0 0 0
2 0 0 30
3 40 0 0
4 30 -20 25
(pulg)
Figura P3.26
Capítulo
Elementos de f lexión
4
Un elemento estructural rectilı́neo esbelto puede ser cargado según la dirección axial propia del
mismo; pero también puede ser cargado transversalmente a ésta dirección. En tal caso decimos que
el elemento está sometido a una solicitación de tipo flexionante, cuyo efecto principal es producir una
deformación que lo convierte en un elemento que adquiere cierto monto de curvatura. Este comporta-
miento es caracterı́stico de un elemento estructural al que comúnmente se lo denomina viga.
En este Capı́tulo, desarrollaremos la formulación mediante el método de elemento finito de una
viga para establecer sus propiedades de rigidez, el tratamiento de la solicitación aplicada para conver-
tir cualquier carga distribuı́da en un sistema de cargas concentradas nodales, ya que ası́ lo requiere el
modelo matemático de análisis. Este procedimiento será efectuado en base a la utilización del concepto
de equivalencia de trabajo mecánico. Finalmente generalizaremos los conceptos elaborados para for-
mular ecuaciones que describan el comportamiento mecánico del elemento rectilı́neo general (en tres
dimensiones) el cual es miembro componente de las denominadas estructuras reticuladas de elementos
lineales.
4.1. Introducción
Los elementos uni–dimensionales, cargados axialmente discutidos en los Capı́tulos 2 y 3 son bas-
tante útiles en el análisis de la respuesta presentada por muchas estructuras simples. Sin embargo, la
restricción marcada que estos elementos no son capaces de transmitir efectos flexionantes (debido a
momentos actuantes) evitan su uso en modelar estructuras que son más comunes de encontrarse, las
que tienen juntas de unión soldadas o remachadas entre elementos. En este capı́tulo, la teorı́a elemental
de la viga es aplicada para desarrollar un elemento flexionante (viga) que sea capaz de exhibir pro-
piamente efectos de flexión transversal o curvatura durante su deformación. El elemento se presenta
primero como una lı́nea (elemento uni–dimensional) capaz de flexionarse en un plano espacial. En el
contexto de desarrollar las ecuaciones discretizadas para este elemento, presentamos un procedimiento
general para determinar las funciones de interpolación usando una forma polinómica supuesta para la
variable de campo. El desarrollo posteriormente se extiende a la flexión en dos–planos y los efectos de
85
86 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
carga axial y torsión se agregan finalmente para lograr plantear un elemento rectilı́neo de comporta-
miento de respuesta mecánica completamente general, con el que se puede modelar cualquier tipo de
estructura espacial conformada con elementos de geometrı́a rectilı́nea.
d
y +M
q = q(x) +M +V
y
ds
x
dx
+V
(a) Viga simplemente apoyada (b) Porción de viga deformada (c) Convención de signos
Las ramificaciones de hipótesis del punto 4 son ilustradas en la Figura 4.2 que muestra dos secciones
transversales que satisfacen la asunción y una sección adicional que no lo hace. Las secciones trans-
versales cuadradas y triangulares son simétricas sobre el plano x–y, y sólo se flexionan en ese mismo
plano. Por otro lado, la sección de forma en L no posee ninguna simetrı́a y su curvatura durante la
deformación queda fuera del plano x–y, incluso bajo carga aplicada sólamente en dicho mismo plano.
Con respecto a la figura, la hipótesis 2 puede cuantificarse aproximadamente para significar que la
deflexión máxima de la viga sea mucho menor que la dimensión de menor longitud lineal de su sección
transversal h. Una regla generalmente aplicable es que la deflexión máxima sea menor que 0,1 h.
Considerado una longitud diferencial dx de una viga después de flexionarse como en la Figura 4.1(b)
la superficie del fondo o inferior ha aumentado en longitud. De aquı́, debe necesariamente existir una
“capa” que debe permanecer indeformada durante la deflexión de la viga. Asumiendo que esta superficie
está localizada a una distancia ρ desde el centro de curvatura O y escogiendo esta capa (que recordemos,
es conocida como la superficie neutra) para corresponderse a la coordenada y = 0, la longitud después
de la deformación flexionante producida en cualquier posición y se expresa como
ds = (ρ − y)dθ (4.1)
4.2. TEORÍA ELEMENTAL DE VIGAS 87
ds − dx (ρ − y)dθ − ρdθ y
x = = =− (4.2)
dx ρdθ ρ
Desde los conceptos básicos de cálculo diferencial, el radio de curvatura de una linea curva plana
está dada por
3/2
dυ
1 + ( )2
dx
ρ= (4.3)
d2 υ
dx2
donde υ = υ(x) representa la curva de deflexión de la superficie neutra.
Siguiendo la teorı́a de la desviación pequeña, las pendientes de las rectas tangentes a la curva de
deformación son también pequeñas, de modo que la Ecuación (4.3) se aproxima por
1
ρ∼
= (4.4)
d2 υ
dx2
luego la deformación normal en la dirección del eje longitudinal como resultado de la flexión producida
en la viga es
d2 υ
x = −y 2 (4.5)
dx
y la tensión normal interna correspondiente es por la ley de Hooke
d2 υ
σx = E x = − E y (4.6)
dx2
donde E es el módulo de elasticidad del material de la viga. La Ecuación (4.6) muestra que, en una
sección transversal dada, la tensión normal varı́a linealmente con la distancia desde la superficie neutra.
Como ninguna fuerza axial neta está actuando en la sección transversal de la viga, la fuerza resultan-
te de la distribución de tensiones normales dada por la Ecuación (4.6) debe ser cero. Por consiguiente,
en cualquier posición axial x a lo largo de la longitud, tendremos
d2 υ
Z Z
Fx = σx dA = − E y 2 dA = 0 (4.7)
A A dx
Notando que en una sección transversal arbitraria la curvatura es constante, la Ecuación (4.7) implica
que necesariamente debe cumplirse Z
y dA = 0 (4.8)
A
88 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
lo cual se satisface si el plano x–z (y = 0) pasa a través del centroide del área plana. Ası́, obtenemos
el resultado muy conocido que la superficie neutra es perpendicular al plano de flexión y atraviesa el
centroide del área de la sección transversal de la viga.
Similarmente, el momento interno de flexión en una sección transversal debe ser equivalente al
momento resultante de la distribución de tensiones normales, entonces
d2 υ
Z Z
M (x) = − y σx dA = E 2 y 2 dA (4.9)
A dx A
El término integral en laR Ecuación (4.9) representa el momento de inercia de la sección transversal
alrededor del eje z, Iz = A y 2 dA, de modo que la expresión del momento flexionante se vuelve
d2 υ
M (x) = E Iz (4.10)
dx2
Combinando las Ecuaciones (4.6) y (4.10), obtenemos la ecuación de tensión normal interna para el
estado de solicitación flexionante de la viga:
M (x) y d2 υ
σx = − = −yE 2 (4.11)
Iz dx
Notemos que el signo negativo en la Ecuación (4.11) asegura que, cuando la viga se sujeta a un
momento flexionante positivo por la convención de signos mostrado en la Figura 4.1(c), se obtienen los
valores de tensión correctamente; es decir, compresivo (negativo) y traccionante (positivo), dependiendo
del signo del valor y de ubicación donde se calcula la tensión normal.
1 2 1 2
1 2
2
1 2 x
1 L
localizados en los extremos del elemento, y las variables nodales son los desplazamientos transversales
υ1 y υ2 en estos nodos, ası́ como también las pendientes (rotaciones) θ1 y θ2 . Las variables nodales como
son mostradas están definidas en la dirección positiva de los ejes utilizados como sistema coordenado,
y deberá ser notado que los ángulos de deformación (rotaciones) serán especificados en radianes. Por
conveniencia, el exponente (e) que indica las propiedades del elemento no se usa por el momento, como
debe ser entendido en el contexto que la discusión actual se aplica y circunscribe a un solo elemento.
Cuando elementos múltiples estén involucrados en los ejemplos a seguir, la notación del exponente se
restaura. La función de desplazamiento υ(x) va a ser discretizada tal que
υ(x = x1 ) = υ1
υ(x = x2 ) = υ2
(4.12a)
dυ
= θ1
dx x=x1
dυ
= θ2
dx x=x2
Antes de proceder, nosotros asumimos que el sistema de coordenadas de elemento (o sistema coor-
denado local) es escogido tal que x1 = 0 y x2 = L, para simplificar la presentación algebráicamente.
(Esto no es para nada restrictivo, desde que L = x2 − x1 en cualquier caso).
Considerado las cuatro condiciones de borde lı́mite y la naturaleza uni–dimensional del problema
por lo que se refiere a la variable independiente, asumimos la función de desplazamiento en la forma
υ(x) = a0 + a1 x + a2 x2 + a3 x3 (4.13)
El escoger una función cúbica para describir el campo de desplazamientos internos al elemento no
es arbitrario. Aunque los requisitos generales de las funciones de interpolación se discutirán en el
Capı́tulo 6, hacemos algunas observaciones pertinentes aquı́. Claramente, con la especificación de cuatro
condiciones de borde, podemos determinar no más de cuatro constantes en la función de desplazamiento
90 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
supuesta. Segundo, en vista de las Ecuaciones (4.10) y (4.13), la segunda derivada de la función de
desplazamiento asumida υ(x) es lineal; entonces, el momento flexionante varı́a linealmente, a lo sumo,
a lo largo de la longitud del elemento. Esto está de acuerdo con la hipótesis que las cargas sólo son
aplicadas en los nodos del elemento, como es indicado por el diagrama de momentos flexionantes de
un elemento cargado mostrado en la Figura 4.5. Si una carga distribuida se aplicara al elemento en
toda su longitud, el momento flexionante variarı́a por lo menos cuadráticamente.
MZ F1L + M2_ M1
F1L _ M1
F1 F2
M2
x
1 2
M1
_
M1
La aplicación de las condiciones de borde especificadas por las Ecuaciones (4.12a) en sucesión
dá como resultado el sistema
υ(x = 0) = υ1 = a0 (4.14a)
2 3
υ(x = L) = υ2 = a0 + a1 L + a2 L + a3 L (4.14b)
dυ
= θ1 = a1 (4.14c)
dx x=0
dυ
= θ2 = a1 + 2a2 L + 3a3 L2 (4.14d)
dx x=L
Las Ecuaciones (4.14) se resuelven simultáneamente para obtener los coeficientes en términos de las
variables de desplazamiento nodal como
a0 = υ1 (4.15a)
a1 = θ1 (4.15b)
3 1
a2 = (υ2 − υ1 ) − (θ2 − θ1 ) (4.15c)
L2 L
2 1
a3 = 3 (υ1 − υ2 ) + (θ1 + θ2 ) (4.15d)
L L2
Substituyendo las Ecuaciones (4.15) en la Ecuación (4.13) y colectando los coeficientes de las variables
nodales mediante procedimiento de factorización, resulta la siguiente expresión que aproxima el campo
de desplazamientos (deflexiones verticales) internos en el elemento
o, en notación matricial,
υ1
θ1
υ(x) = N1 N2 N3 N4 = [ N (x) ]{δ}] (4.17a)
υ2
θ2
Para el elemento de flexión (elemento viga), es conveniente introducir una nueva variable definida
como la coordenada de posición longitudinal adimensional, mediante
x
ξ= (4.19)
L
En función de esta nueva variable, la Ecuación (4.16) se transforma en
donde 0 ≤ ξ ≤ 1. Esta forma provee caracterı́stica de facilidad de evaluación de las integraciones re-
queridas para completar el desarrollo de las ecuaciones de elemento, que serán obtenidas en la próxima
sección.
Como fué discutido en el Capı́tulo 3, los desplazamientos son importantes, pero el ingeniero está más
a menudo interesado en examinar las tensiones asociadas con las condiciones de carga. Usando la
Ecuación (4.11) en conjunción con la Ecuación (4.17a), nos permite obtener la distribución de tensiones
normales internas para una sección transversal localizada en una posición axial genérica x, la cual se
presenta como
d2 [N ]
σx (x, y) = − y E {δ} (4.21)
dx2
Puesto que la tensión normal interna varı́a linealmente en una sección transversal (sección de corte),
los valores máximo y mı́nimo en cualquier posición axial ocurren sobre las superficies exteriores del
elemento, donde la distancia y medida desde la superficie neutra es de mayor magnitud, o la más grande.
Como es de costumbre, nosotros tomamos la tensión máxima como el valor más grande de tensión de
tracción (positiva) y el valor mı́nimo como el valor más grande (negativo) de tensión compresiva. De
aquı́, volvemos a escribir la Ecuación (4.21) como
d2 [N ]
σx = ymáx E {δ} (4.22)
dx2
y será entendido que la Ecuación (4.22) representa los valores máximo y mı́nimo de tensión normal
en cualquier sección transversal definida por la coordenada axial x. También ymax representa a las
distancias más grandes (una positiva, otra negativa) medidas desde la superficie neutra hasta las
92 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
superficies externas (superior e inferior) del elemento. Sustituyendo las funciones de interpolación y
llevando a cabo las diferenciaciones indicadas, obtenemos
12x 6 6x 4 6 12x
σ(x) = ymáx E − 2 υ1 + − θ1 + − 3 υ2
L3 L L2 L l2 L
(4.23)
6x 2
+ − θ2
L2 L
Observando que la Ecuación (4.23) indica una variación lineal de la tensión normal a lo largo de la
longitud del elemento y subsecuentemente, una vez que la solución de los desplazamientos se obtiene,
los valores nodales son constantes conocidas; ası́ necesitamos calcular sólamente los valores de tensión
en las secciones transversales correspondientes a los nodos; es decir, en x = 0 y x = L correspondientes
a los extremos del elemento. Los valores de tensión en las secciones nodales están dadas por
6 2
σx (x = 0) = ymáx E (υ 2 − υ 1 ) − (2θ 1 + θ 2 ) (4.24a)
L2 L
6 2
σx (x = L) = ymáx E (υ1 − υ2 ) + (2θ2 + θ1 ) (4.24b)
L2 L
Los cálculos de tensión normal interna se ilustran en los ejemplos desarrollados en páginas siguientes.
Nuevamente reconocemos la integral de área como el momento de inercia Iz alrededor del eje centroidal
perpendicular al plano de flexión, y por ello tendremos
2
E Iz L d 2 υ
Z
Ue = dx (4.28)
2 0 dx2
La Ecuación (4.28) representa la energı́a de deformación flexionante para cualquier sección trans-
versal constante de la viga que obedece las hipótesis de la teorı́a elemental de vigas rectas. Para la
energı́a de deformación del elemento finito que está desarrollándose, sustituimos la relación discretizada
de desplazamiento de la Ecuación (4.17) para obtener
L 2
d2 N1 d2 N2 d2 N3 d2 N4
Z
E Iz
Ue = υ 1 + θ 1 + υ 2 + θ2 dx (4.29)
2 0 dx2 dx2 dx2 dx2
4.4. MATRIZ DE RIGIDEZ DE ELEMENTO VIGA 93
como aproximación para la energı́a de deformación elástica. Enfatizamos que la Ecuación (4.29) es una
aproximación porque la función de desplazamiento discretizada no es en general una solución exacta
para el problema de flexión de una viga.
Aplicando el teorema de Castigliano a la función de energı́a de deformación con respecto al despla-
zamiento nodal υ1 , se obtiene la fuerza transversal en el nodo 1 como
L
d2 N1 d2 N2 d2 N3 d2 N4 d2 N1
Z
∂Ue
= F1 = E Iz 2
υ1 + 2
θ1 + 2
υ2 + θ2 dx (4.30)
∂υ1 0 dx dx dx dx2 dx2
mientras que la aplicación del teorema con respecto al desplazamiento rotacional en el mismo nodo
proporciona el momento flexionante como
L
d2 N1 d2 N2 d2 N3 d2 N4 d2 N2
Z
∂Ue
= M1 = E Iz 2
υ1 + 2
θ1 + 2
υ2 + θ2 dx (4.31)
∂θ1 0 dx dx dx dx2 dx2
L
d2 N1 d2 N2 d2 N3 d2 N4 d2 N4
Z
∂Ue
= M2 = E Iz 2
υ1 + 2
θ1 + 2
υ2 + θ2 dx (4.33)
∂θ2 0 dx dx dx dx2 dx2
Las Ecuaciones (4.30) a (4.33) algebráicamente relacionan los cuatro valores de desplazamientos nodales
a las cuatro fuerzas nodales aplicadas (aquı́ usamos el vocablo fuerza en un sentido general, el cual
también incluye a los momentos aplicados), y son de la forma
k11 k12 k13 k14 υ1 F1
k21
k22 k23 k24 θ M
1 1
k31
= (4.34)
k32 k33 k34
υ2 F2
k41 k42 k43 k44 θ2 M2
donde kmn (m, n = 1; 4) son los coeficientes de la matriz de rigidez de elemento. Por comparación
de las Ecuaciones (4.30) a (4.33) con las ecuaciones algebráicas representadas matricialmente por la
Ecuación (4.34), podemos apreciar que se cumple la ecuación recursiva
L
d 2 Nm d2 Nn
Z
kmn = knm = E Iz dx m, n = 1; 4 (4.35)
0 dx2 dx2
y la matriz de rigidez de elemento es matriz simétrica, como deberı́a esperarse por el comportamiento
mecánico lineal elástico del elemento.
Antes de calcular los coeficientes de la matriz de rigidez, es conveniente convertir la integración
hacia la variable de longitud adimensional ξ = x/L, notando que
Z L Z 1
d 1 d
f (x) dx = f (ξ) dξ y =
0 0 dx L dξ
L 1
d 2 Nm d2 Nn d2 Nm d2 Nn
Z Z
E Iz
kmn = knm = E Iz 2 2
dx = 3 dξ m, n = 1; 4 (4.36)
0 dx dx L 0 dξ 2 dξ 2
94 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
E Iz L 36 E Iz L 2
Z Z
2
k11 = 3 (12ξ − 6) dξ = (4ξ − 4ξ + 1)dξ
L 0 L3 0
36 E Iz 4 12 E Iz
= −2+1 =
L3 3 L3
Z L
E Iz 6 E Iz
k12 = k21 = 3 (12ξ − 6)(6ξ − 4)Ldξ =
L 0 L2
Z L
E Iz 12 E Iz
k13 = k31 = 3 (12ξ − 6)(6 − 12ξ)dξ = −
L 0 L3
Z L
E Iz 6 E Iz
k14 = k41 = 3 (12ξ − 6)(6ξ − 2)Ldξ =
L 0 L2
Continuando con las integraciones, se demuestra que los coeficientes de rigidez faltantes son
4 E Iz 6 E Iz 2 E Iz
k22 = k23 = k32 = − k24 = k42 =
L L2 L
12 E Iz 6 E Iz 4 E Iz
k33 = k34 = k43 =− k44 =
L3 L3 L
La matriz de rigidez completa para el elemento viga sometido a flexión es luego escrita como
12 6L −12 6L
E Iz 6L 4L2 −6L 2L2
[ ke ] = 3 (4.37)
L −12 −6L 12 −6L
2 2
6L 2L −6L 4L
La simetrı́a de la matriz de rigidez de elemento está clara, como previamente fué observado. De
nuevo, la matriz de rigidez de elemento se muestra que es singular desde que es posible el movimiento
de cuerpo rı́gido, a menos que el elemento esté restringido de alguna manera a poder desplazarse. La
matriz de rigidez como es dada por la Ecuación (4.37) es válida en cualquier sistema consistente de
unidades de medida, a condición que los grados de libertad rotacionales (pendientes de la deformación)
se expresen en radianes.
P0
F1 F2
M2 M0
1 2 x
M1
V
(a) Convención de cargas positivas
método de elemento finito
_ } P0
+
x
M1 V2 M2
1 2
M
M0 }
V1
+
(b) Convención de cargas positivas x
mecánica de sólidos
(c) Diagramas de esfuerzos internos
Ejemplo 4.1.
La Figura 4.7(a) muestra una viga estáticamentye indeterminada sujeta a una carga transversal apli-
cada en la mitad de su longitud. Usando dos elementos finitos de flexión, obtener una solución para la
deflexión vertical de la sección media longitudinal de la viga (donde la carga se aplica).
> Solución
Puesto que los elementos finitos de flexión requieren que las cargas aplicadas actúen exclusivamente
en sus puntos nodales (es decir, en sus extremos), los elementos a definirse necesariamente deberán ser
de longitud L/2 como se muestra en la Figura 4.7(b), para que la carga concentrada externa aplicada
96 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
y
P 1 2 3
3
2
x
1 1 2 2 3
L L
2 2 1
(a) Esquema de solicitación y apoyos (b) Codificación estructural adoptada
y
1 1 0 2 3 0 P
3PL/16 1
3
2 3 L/2 2 L/2 x
11P/16 5P/16
sea una carga nodal. Las matrices de rigidez de los elementos individuales son entonces
12 6L/2 −12 6L/2
E Iz 2 2
6L/2 4L /4 −6L/2 2L /4
[ k (1) ] = [ k (2) ] =
(L/2) 3 −12 −6L/2 12 −6L/2
6L/2 2L2 /4 −6L/2 4L2 /4
12 3L −12 3L
8 E Iz
3L L2 −3L L2 /2
=
L3 −12 −3L 12 −3L
3L L2 /2 −3L L2
Notemos que particularmente la longitud de cada elemento es L/2, y dicha longitud fué usada en la
valoración de las matrices de rigidez locales de ambos elementos. Las condiciones de borde apropiadas
en este caso, con referencia a la Figura 4.7(b) donde se definen los grados de libertad globales, son:
υ1 = θ1 = υ3 = 0; mientras que la correspondencia de desplazamientos locales–hacia–globales se
muestra en la Tabla 4.1.
1 1 0
2 2 0
3 3 1
4 4 2
5 0 3
6 0 4
(1)
96 E Iz (1)
4 E Iz (1) (2)
16 E Iz
K11 = k11 = K24 = k24 = K44 = k44 + k22 =
L3 L L
(1)
24 E Iz (2)
−24 E Iz
K12 = k12 = K25 = K26 = 0 K45 = k23 =
L2 L2
4.5. VECTOR DE CARGA DE ELEMENTO 97
(1)
−96 E Iz (1) (2)
192 E Iz (2)
4 E Iz
K13 = k13 = K33 = k33 + k11 = K46 = k24 =
L3 L3 L
(1)
24 E Iz (1) (2) (2)
96 E Iz
K14 = k14 = K34 = k34 + k12 =0 K55 = k33 =
L2 L3
(1)
8 E Iz (2)
−96 E Iz (2)
−24 E Iz
K22 = k22 = K35 = k13 = K56 = k34 =
L L3 L2
(1)
−24 E Iz (2)
24 E Iz (2)
8 E Iz
K23 = k23 = K36 = k14 = K66 = k44 =
L2 L2 L
Usando la forma general [ K ]{U } = {F } para el modelo matemático de análisis, obtenemos el sistema
de ecuaciones como
96 24L −96 24L 0 0 υ1 F1
24L 8L2 −24L 4L2 0 0 θ1
M1
E Iz
−96 −24L 192 0 −96 24L
υ2
F2
2
=
L3 24L 4L 0 16L2 −24L 4L2
θ2 M2
0 0 −96 −24L 96 24L υ F3
3
0 0 24L 4L2 24L 8L2 θ3 M3
Nótese que aquı́ hemos incorporado la fuerza vertical externa aplicada en el nodo 2, la cual actúa con
sentido opuesto al desplazamiento vertical asociado en dicho nodo estructural (que establece la dirección
positiva en dicho punto); por ello el signo negativo para esta acción actuante sobre la estructura en
nuestro modelo (F2 = −P ). La solución simultánea del sistema previamente planteado, como se puede
comprobar, dá como resultado
−7 P L3 −P L2 P L2
υ2 = θ2 = θ3 =
768 E Iz 128 E Iz 32 E Iz
Un bosquejo de la deformación de la viga se muestra en la Figura 4.7(c), donde la lı́nea horizontal
segmentada representa a la estructura indeformada, el que se obtiene considerando los valores de
solución para los desplazamientos nodales previamente encontrados.
Ahora, la substitución de los valores de los desplazamientos nodales hallados como solución previa
en las ecuaciones de restricción (aquellas ecuaciones que fueron eliminadas en el proceso de reducción
del modelo), permite evaluar las fuerzas de reacción de apoyo según el siguiente procedimiento
E Iz 11P
F1 = (−96υ2 + 24Lθ2 ) =
L3 16
E Iz 5P
F3 = 3 (−96υ2 − 24Lθ2 − 24Lθ3 ) =
L 16
E Iz 3P L
M1 = 3 (−24Lυ2 + 4L2 θ2 ) =
L 16
Estos valores obtenidos nos permiten efectuar una verificación del cumplimiento de las ecuaciones
de equilibrio estático global de toda la estructura considerada como sistema. Con referencia a la
Figura 4.7(d), tenemos sumando algebráicamente las fuerzas
P 11P 5P
Fy = −P + =0
16 16
98 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
y sumando los momentos con relación a un eje perpendicular a la Figura que pase por el nodo 1,
P 3P L L 5P
Mz = −P + L=0
16 2 16
Luego, la solución de elemento finito satisface exáctamente las ecuaciones de equilibrio para la estruc-
tura analizada, lo cual es indicativo que no se cometieron errores en el procedimiento de aplicación del
método; y también permite aceptar la solución obtenida como satisfactoria. >
El lector astuto puede desear comparar los resultados del Ejemplo 4.1 con aquellos obtenidos de
alguna tabla de deflexiones estandar de vigas en diferentes disposiciones de apoyo y solicitación; en
cuyo caso se encontrará que los resultados están en exacto acuerdo con la teorı́a elemental de vigas
rectas. En general, el método de elemento finito es un método aproximado; pero en el caso de flexión
simple, los resultados son exactos en ciertos casos. En este ejemplo, la ecuación de la deflexión de
la superficie neutra es una ecuación cúbica y, puesto que las funciones de interpolación utilizadas
son cúbicas también, los resultados son exactos. Cuando existen cargas distribuı́das, sin embargo, los
resultados no son necesariamente exactos, como se mostrará y discutirá luego.
L 12 12
Refiriéndonos a la Figura 4.8, el trabajo mecánico realizado por la carga distribuida original puede
expresarse como
Z L
W = q(x)υ(x)dx (4.39)
0
donde q(x) es la función de carga linealmente distribuı́da actuante al interior del elemento. Esta función
es una relación que describe de forma general la carga aplicada al interior del elemento. (Es posible
escribir esta función con la metodologı́a de uso de funciones singulares, que permite incorporar fuerzas
y momentos flectores puntuales que actúen a interior del elemento).
El objetivo aquı́ es determinar las cargas nodales equivalentes (cargas concentradas actuando en
los extremos del elemento), de modo que el trabajo mecánico efectuado por estas acciones sea el mismo
que el trabajo en el sistema original. El trabajo mecánico de las cargas equivalentes es
donde F1eq , F2eq son las fuerzas nodales equivalentes en los nodos 1 y 2, respectivamente; y M1eq , M2eq
son los momentos nodales equivalentes. Substituyendo la función de desplazamiento discretizada dada
4.6. CARGAS NODALES EQUIVALENTES 99
o, desarrollando términos
! !
Z L Z L
W = q(x)N1 (x)dx υ1 + q(x)N2 (x)dx θ1
0 0
! ! (4.41)
Z L Z L
+ q(x)N3 (x)dx υ2 + q(x)N4 (x)dx θ2
0 0
La comparación de las Ecuaciones (4.40) y (4.41), por la hipótesis planteada W = Weq , nos muestra
que se deberán cumplir las relaciones siguientes
Z L
F1eq = q(x)N1 (x)dx (4.42a)
0
Z L
M1eq = q(x)N2 (x)dx (4.42b)
0
Z L
F2eq = q(x)N3 (x)dx (4.42c)
0
Z L
M2eq = q(x)N4 (x)dx (4.42d)
0
De aquı́, el vector de fuerza nodal que representa una carga distribuı́da arbitraria en base a la
equivalencia de trabajo mecánico está dado por las Ecuaciones (4.42). Por ejemplo, para una carga
uniforme aplicada al elemento q(x) = q constante, incorporando las funciones de interpolación indicadas
por las Ecuaciones (4.18); y efectuando las integraciones como se indican en las ecuaciones precedentes
dá como resultado
qL
F1eq 2
2
qL
M1eq
12
= (4.43)
F2eq qL
2
M2eq
−qL2
12
Ejemplo 4.2.
La viga simplemente apoyada mostrada en la Figura 4.9(a) está solicitada por una fuerza transversal
linealmente distribuida uniforme. Usando dos elementos de igual longitud y las cargas nodales equiva-
lentes de trabajo mecánico, obtenga una solución de elemento finito para la deflexión a la mitad de la
longitud de la viga y comparar la solución obtenida con aquella que proporciona la teorı́a elemental
de vigas rectas.
> Solución
Primero numeramos los nodos y elementos como es mostrado por la Figura 4.9(b), y también notamos
que las condiciones de borde lı́mite son en este caso: υ1 = υ3 = 0. También podrı́amos notar la
condición de simetrı́a de deformación que impone θ2 = 0. Sin embargo, en este caso, permitimos que
100 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
y
q 1 2 3
3
2
x 1 1 2 2 3
1
L
(a) Esquema de solicitación y apoyos (b) Codificación estructural adoptada
2
qL
2 qL 2
q q 48
qL
48 48
1 2 2 3
1 2 3
L L qL qL qL
2 2 4 4
4
Figura 4.9: Viga simplemente apoyada cargada con fuerza distribuida uniforme
ese hecho ocurra como resultado del proceso de solución. Las matrices de rigidez de ambos elementos
son idénticas, dadas por
12 6L/2 −12 6L/2
E Iz 2 2
6L/2 4L /4 −6L/2 2L /4
[ k (1) ] = [ k (2) ] =
(L/2) 3 −12 −6L/2 12 −6L/2
6L/2 2L /4 −6L/2 4L2 /4
2
12 3L −12 3L
8 E Iz
3L L2 −3L L2 /2
=
L3 −12 −3L 12 −3L
3L L2 /2 −3L L2
(nuevamente notamos que hemos usado una longitud igual a L/2 para evaluar los términos de la
matriz de rigidez de ambos elementos). La correspondencia de desplazamientos de elemento y despla-
zamientos globales asignados a la estructura, es mostrada en la Tabla 4.2 que describe las relaciones
de conectividad.
1 1 0
2 2 0
3 3 1
4 4 2
5 0 3
6 0 4
Ahora, usando como referencia la Tabla 4.2, podemos proceder a efectuar el proceso de ensamblaje
de coeficientes de la matriz de rigidez global de la estructura. Se demuestra que la aplicación del pro-
cedimiento mencionado, produce la matriz de rigidez global de la viga que fácilmente usted comprueba
4.6. CARGAS NODALES EQUIVALENTES 101
que es
12 3L −12 3L 0 0
3L
L2 −3L L2 /2 0 0
8 E Iz −12 −3L 24 0 −12 3L
[K ] = 2 2 2
(♠)
L3 3L L /2 0 2L −3L L /2
0 0 −12 −3L 12 3L
0 0 3L L2 /2 3L L2
Las cargas de trabajo – equivalente para cada elemento se calculan en base al esquema presentado
en la Figura 4.9(c), y las cargas nodales equivalentes resultantes se muestran en la Figura 4.9(d). Estas
cargas se calcularon adaptando la Ecuación (4.43) tomando q(x) = −q (ya que la carga actúa en
dirección y negativa), y longitud igual a L/2. Observando que existen reacciones de apoyo en los nodos
1 y 3, en adición a las cargas nodales equivalentes calculadas como solicitación externa proveniente de
la carga distribuida aplicada, las ecuaciones globales de equilibrio llegan a ser
−qL
4 + F1
υ1
−qL2
θ
48
1
−qL
υ2
[K ] = 2
θ2
0
υ
−qL
3
4 + F
3
θ3
qL2
48
la cual, en solución simultánea de incógnitas (por ejemplo aplicando la técnica de inversión matricial
ayudado mediante un paquete de evaluación matricial como MathCad, que tiene capacidad de operar
algebráicamente arreglos matriciales), dá como resultado
qL3 5qL4 qL3
θ1 = − θ2 = 0 υ2 = − θ3 =
24EIz 384EIz 24EIz
Como es esperado, la pendiente de deformación de la viga a la mitad de su longitud es cero, y
puesto que la carga y las condiciones de apoyo son simétricas, la solución de la deformación también
es simétrica, como se aprecia por los valores de magnitud de las rotaciones en los extremos de la viga.
El desplazamiento nodal planteado como incógnita primaria del problema, obtenido como solución
mediante aplicación del análisis de elemento finito en este ejemplo, es exactamente el mismo resultado
que se obtiene por la teorı́a básica de la mecánica de sólidos deformables. Esto es debido a la aplicación
de las cargas nodales equivalentes obtenidas por el principio de equivalencia de trabajo mecánico. Sin
embargo, la forma de deflexión general como es dada por la solución de elemento finito no es la misma
que el resultado que proporciona la mecánica de sólidos deformables. La ecuación que describe la
desviación de la superficie neutra es una función cuadrática de la posición x y, puesto que las funciones
de interpolación usadas en el modelo de elemento finito son cúbicas, la curva de deflexión varı́a en algo
con respecto a la solución que podemos considerar exacta. >
102 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
Ejemplo 4.3.
En la Figura 4.10(a), la viga OC es soportada por una conexión articulada lisa en O y también en B
por una varilla elástica BD, que a su vez se inmoviliza en el punto D a través de conexiones articuladas.
Una carga concentrada F = 10 kN es aplicada en C. Determine la deflexión del punto C y la tensión
axial en el miembro BD. El módulo de elasticidad de la viga es 207 GPa (acero) y la dimensión de
la arista de su sección transversal cuadrada es de 40 mm. Para la varilla elástica BD, el módulo de
elasticidad es 69 GPa (aluminio) y el área de su sección transversal es de 78,54 mm2 .
D Y U7
200 mm 10 kN
U1 3 U5
U3
O B C
300 mm 300 mm U2 1 U4 2 U6 X
(a) Esquema de solicitación y apoyos (b) Sistema coordenado global y codificación estructural
u2(3)
2(3)
(1)
1 2(1) (1) u1(3) 1(2) 2(2) (2)
2 2
(1)
1(3) (2)
1 1
(c) Desplazamientos individuales de elemento
> Solución
Este es el primer ejemplo en el que nosotros usamos tipos múltiples de elementos, cuando la viga es
modelada con elementos de flexión y la varilla elástica como un elemento barra. Claramente, el miembro
horizontal está sujeto a cargas transversales flexionantes, de modo que las asunciones de elemento barra
no se aplican a este miembro. Por otro lado, el miembro de apoyo vertical está sujeto sólo a carga axial,
puesto que las conexiones articuladas no pueden transmitir momentos. Por consiguiente, usamos dos
tipos de elementos diferentes para simplificar el modelado y la solución.
El sistema coordenado global y las variables globales se muestran en la Figura 4.10(b), donde el
sistema es dividido en dos elementos flexionantes (1 y 2) y un elemento barra (3). Para los propósitos
de etiquetado numérico en la matriz de rigidez global, se usa el esquema de desplazamientos detallado
en la Tabla 4.3.
Aunque la notación mostrada en la Figura 4.10(b) puede parecer ser incoherente con la notación
previa, es más simple en términos de numerar sucesivamente los desplazamientos en las ecuaciones
globales. Por la asignación apropiada de desplazamientos de elemento con los desplazamientos globales
(relaciones de conectividad), la distinción entre los desplazamientos lineales y rotatorios está clara. Los
desplazamientos de elemento individual se muestran en la Figura 4.10(c), donde claramente indicamos
al elemento barra en su configuración general bi–dimensional; aunque, en este caso, nosotros sabemos
que υ1(3) = υ2(3) = 0 y esos desplazamientos se ignoran en el proceso de búsqueda de la solución.
La correspondencia de desplazamientos de elemento se muestra en la Tabla 4.4. Para los elementos
4.6. CARGAS NODALES EQUIVALENTES 103
1 U1 υ1(1) 0 0
2 U2 θ1(1) 0 0
3 U3 υ2(1) υ1(2) u(3)
1
4 U4 θ2(1) θ1(2) 0
5 U5 0 υ2(2) 0
6 U6 0 θ2(2) 0
7 U7 0 0 u(3)
2
1 1 0 0
2 2 0 0
3 3 1 1
4 4 2 0
5 0 3 0
6 0 4 0
7 0 0 3
viga, el momento de inercia alrededor del eje de z es calculado en base a la fórmula general conocida
para el momento de inercia centroidal de un rectángulo
bh3 a4 404
Iz = = = = 213333 mm4
12 12 12
Para los elementos 1 y 2,
EIz (207×103 )213333
k1 = k2 = = = 1635,6 N/mm
L3 3003
de modo que las matrices de rigidez locales de elemento son [por la Ecuación (4.37)]
12 1800 −12 1800
1800 360000 −1800 180000
[ k (1) ] = [ k (2) ] = 1635,6
−12 −1800
12 −1800
1800 180000 −1800 360000
1 2 3 4 5 6 7
6 6
19627,2 2,944×10 −19627,2 2,944×10 0 0 0 1
2,944×106 5,888×108 −2,944×106 2,944×108 0 0 0 2
−19627,2 −2,944×106 66350,4 0 −19627,2 2,944×106 −27096 3
[K ] =
2,944×106 2,944×108 0 11,78×108 −2,944×106 2,944×108 0
4
0 0 −19627,2 −2,944×106 19627,2 −2,944×106 0 5
0 0 2,944×106 2,944×108 −2,944×106 5,889×108 0 6
0 0 −27096 0 0 0 27096 7
(f)
Las condiciones de restricción de desplazamientos son U1 = U7 = 0, y el vector de fuerza externa
aplicada sobre la estructura es
F1
R1
M1 0
F2 0
M2 = 0
F3 −10000
M 0
3
F4 R4
donde usamos R para indicar una fuerza de reacción de apoyo. Si aplicamos las condiciones de res-
tricción de desplazamientos que se presentan en los apoyos de la estructura, obtenemos una ecuación
reducida pues debemos desechar dos filas y dos columnas del sistema original gobernante del compor-
tamiento mecánico. Reduciendo el sistema y hallando la solución para los desplazamientos activos, se
demuestra que obtenemos
θ1 = 9,3638×10−4 rad
υ2 = −0,73811 mm
θ2 = −0,0092538 rad
υ3 = −5,5523 mm
θ3 = −0,019444 rad
(Notemos que intencionalmente tomamos más dı́gitos decimales significativos de los necesarios, para
evitar inexactitudes de “redondeo de fracciones” en los cálculos secundarios). Para obtener la tensión
axial en el miembro BD, utilizamos la Ecuación (3.44) con θ(3) = π/2:
0
3 1
−0,7381
σBD = 69×10 0 −1 0 1 = 254,65 MPa
200
0
0
El resultado positivo nos indica tensión de tracción en este elemento barra (la varilla de aluminio).
Las fuerzas de reacción de apoyo se obtienen por substitución de los desplazamientos calculados en
la primera y séptima ecuaciones (las ecuaciones de restricción que fueron desechadas en el proceso de
reducción del sistema):
y dentro de la exactitud numérica usada en este ejemplo, el sistema está en equilibrio estático. Se insiste
al lector a verificar el equilibrio de momentos con relación al nodo del extremo izquierdo de la viga y
notar que, por estática simplemente, la fuerza en el elemento 3 debe ser 20,000 N y la tensión normal
axial calculada por σ = F/A resultarı́a ser 254,65 MPa (el mismo valor que se obtuvo anteriormente).
Las tensiones normales provenientes de la flexión en los nodos extremos de los elementos que
modelan la viga son evaluadas mediante las Ecuaciones (4.24), notando que para la sección transversal
cuadrada ymáx/mı́n = 20 mm. Para el elemento 1,
6 2
σx(1) (x = 0) = ±20(207)(103 ) (−0,738 − 0) − (−(2)0,00092 − 0,0092) ≈ 0 MPa
3002 300
en el nodo 1. Debemos notar que la tensión calculada en el nodo 1 deberı́a ser idénticamente cero,
porque este nodo coincide con un apoyo articulado fijo que no puede soportar momentos flexionantes.
Para el nodo 2 del elemento 1, hallamos
6 2
σx(1) (x = L) = ±20(207)(103 ) (0 + 0,738) + (−(2)0,00092 − 0,0092) ≈ ±281,3 MPa
3002 300
y los últimos resultados también son los esperados, porque el extremo derecho de la viga está libre
de momento flexionante externo actuante en dicho punto. Aquı́ necesitamos observar cuidadosamente
que las tensiones normales provenientes de la flexión en la viga tienen idéntico valor en el nodo de
interconexión 2 entre elementos (lo que indica cumplimiento continuidad de la tensión normal interna en
la frontera lı́mite de estos elementos, en este caso particular). Ésta no es la situación usual en el análisis
por elemento finito. La formulación efectuada requiere continuidad de desplazamientos (deformaciones
transversales) y rotaciones (pendientes de la deformación) pero, en general, no exige continuidad de
las derivadas de orden superior al primero para la variable de campo de interés que en este caso es la
deflexión transversal.
Puesto que el elemento de flexión desarrollado aquı́ está basado en una función de desplazamiento
cúbica, el elemento no exhibe a menudo continuidad en el momento flector interno que es diréctamente
proporcional a la segunda derivada de la deflexión (y por tanto tampoco en la tensión normal interna
— que depende del momento flector). La convergencia de las derivadas de funciones es relevante para
examinar la exactitud de una solución de elemento finito a un problema dado. Nosotros debemos
examinar el comportamiento numérico de las variables derivadas a medida que el modelo de elemento
finito esté asociado a una “malla” que es refinada (cuando se usa un mayor número de elementos de
tamaño más pequeño). >
El ejemplo anterior muestra que es posible combinar elementos de distinto comportamiento mecáni-
co en una sola estructura, y que el procedimiento de pasos secuenciales que maneja el método de
elemento finito se mantiene invariable; es decir que la secuencia de cálculo numérico al tratar cual-
quier elemento especı́fico es el mismo, obviamente respetando sus peculiaridades de comportamiento
mecánico. El proceso de ensamblaje estructural es exáctamente el mismo que fué aplicado en ejemplos
anteriores, y el procedimiento de solución del modelo matemático de análsis (ecuación gobernante) es
de hecho invariable.
106 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
i j
j
i ui j
uj
i
M
M M
(x)
x F
F F F
M(x) = M - F(x) M(x)
(a) Solicitación sobre elemento marco (b) Efecto mecánico del esfuerzo axial
Siendo este el caso, nosotros simplemente podemos agregar el comportamiento de rigidez axial del
elemento barra al comportamiento de rigidez flexionante del elemento viga para obtener la matriz
de rigidez de dimensión 6×6 (por el número de desplazamientos nodales – véase la Figura 4.11) del
4.7. EL ELEMENTO MARCO 107
que representa una superposición no–acoplada de rigideces axial y transveral del elemento lineal que
estamos desarrollando.
Agregando la capacidad axial al elemento viga se elimina la restricción que tales elementos se
alineen linealmente y habilitan el uso del elemento en el análisis de estructuras tipo marco plano
(llamados también pórticos) en el que las juntas de unión poseen resistencia a la flexión; es decir
transmiten momentos flexionantes. Para tales aplicaciones, la orientación del elemento en el sistema de
coordenadas global debe ser considerada, como era el caso con el elemento barra de una cercha plana.
2 U5
x U6
2 u2 Y U4
y
U2
1
u1
U3 U1 X
1
(a) Referidos al sistema coordenado local (b) Referidos al sistema coordenado global
La Figura 4.13(a) muestra un elemento marco orientado a un ángulo ψ arbitrario respecto del eje
X de un sistema referencial coordenado global, y muestra los desplazamientos nodales del elemento
referidos a su propio sistema coordenado. Aquı́, nosotros usamos el sı́mbolo ψ para indicar el ángulo de
orientación para evitar confusión con la deformación angular nodal, denotada como θ. La Figura 4.13(b)
muestra los desplazamientos globales asignados al elemento, donde de nuevo hemos adoptado un solo
sı́mbolo para los desplazamientos con un sub–ı́ndice que se incrementa numéricamente de nodo a nodo.
Antes de proceder, notemos que es conveniente aquı́ re–ordenar los coeficientes de la matriz de rigidez
de elemento dada por la Ecuación (4.44), para que el vector de desplazamientos de elemento en el
marco de referencia local sea establecido como
u1
υ1
θ1
{δ} = (4.45)
u2
υ2
θ2
108 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
10 lb/pulg Y
u2
B U5 U8
C U9 v2
U4
20 pulg 2
2 U7
U6
20 pulg
1
u1
U2
v1
1
O U3 U1 X
(a) Solicitación aplicada y apoyos (b) Codificación estructural adoptada (c) Elemento 1
y la rotación en B.
> Solución
Usando los datos especificados, el área de sección transversal es A = 1×1 = 1 pulg2 , y el momento de
inercia alrededor del eje z centroidal es Iz = a4 /12 = 14 /12 = 0,083 pulg4 . Denotando el miembro OB
como el elemento 1 y el miembro BC como el elemento 2, se vé que ambos tienen idéntica longitud
L = 20 pulg. Las matrices de rigidez locales, o sea referidas a los sistemas coordenados propios de
cada elemento son idénticas y están definidas por la Ecuación (4.46). Evaluando numéricamente con
los datos anteriores, tendrı́amos
5(105 ) −5(105 )
0 0 0 0
0 1250,4 12504 0 −1250,4 12504
0 12504 166720 0 −12504 83360
[ k (1) ] = [ k (2) ] =
−5(105 )
0 0 5(105 ) 0 0
0 −1250,4 −12504 0 1250,4 −12504
0 12504 83360 0 −12504 166720
Note particularmente cómo la matriz de rigidez del elemento 1 cambia como resultado de la rotación
de 90◦ . Los valores de los coeficientes individuales en la matriz de rigidez están inalterados, pero su
ubicación se altera. La posición de los términos en la matriz de rigidez se cambia para reflejar, bastante
simplemente, las direcciones de los desplazamientos transversales asociados con la flexión y aquellos
110 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
1 1 0
2 2 0
3 3 0
4 4 1
5 5 2
6 6 3
7 0 4
8 0 5
9 0 6
1 2 3 4 5 6 7 8 9
1250,4 0 −12504 1250,4 0 −12504 0 0 0 1
0 500000 0 0 −500000 0 0 0 0
2
−12504 0 166720 12504 0 833360 0 0 0 3
1250,4
0 1250 501250,4 0 12504 −500000 0 0
4
[ K ]=
0 −500000 0 0 501250,4 12504 0 −1250,4 12504
5
−12504 0 83360 12504 12504 333440 0 −12504 83360 6
0 0 0 −500000 0 0 500000 0 0
7
0 0 0 0 −1250,4 −12504 0 1250,4 −12504 8
0 0 0 0 12504 83360 0 −12504 166720 9
(?)
Usando la matriz de rigidez de sistema calculada previamente, y usando los resultados obtenidos
en la Ecuación (4.43) para las cargas nodales equivalentes a la carga distribuida actuante sobre el
elemento 2, las ecuaciones de equilibrio de la estructura congregadas por el ensamble efectuado son
RX1
U1
R
X1
U2
R
Y 1
U
3
M
R1
U
4
0
[ K ] U5 =
−100
U6
−333,3
U
7
R
X3
U
8
R − 100
Y 3
U9
MR3 + 333,3
donde denotamos las fuerzas en los nodos 1 y 3 como las componentes de reacción, debidas a las
restricciones de desplazamiento U1 = U2 = U3 = U7 = U8 = U9 = 0. Tomando éstas restricciones en
4.8. EL ELEMENTO RECTILÍNEO GENERAL 111
consideración, las ecuaciones a ser resueltas para los desplazamientos activos son entonces
501250,4 0 12504 U4 0
0 501250,4 12504 U5 = −100
12504 12504 333440 U6 333,3
U4 = 2,47974(10−5 ) pulg
U5 = −1,74704(10−4 ) pulg
U6 = −9,94058(10−4 ) rad
Como de costumbre, las componentes de reacción de apoyo pueden obtenerse sustituyendo los despla-
zamientos calculados en las seis ecuaciones de restricción que aún no fueron consideradas.
Para el elemento marco (viga con capacidad axial), el cálculo de tensión normal interna debe
tomar en consideración la superposición de tensión normal proveniente de la flexión y la tensión axial
directa. Para el elemento 1, por ejemplo, usamos la Ecuación (4.47) con ψ = π/2 para calcular los
desplazamientos locales del elemento como
u1 0 1 0 0 0 0 U1
0
υ1
−1 0 0 0 0 0
U2
0
θ1 0 0 1 0 0 0 U3 0
= =
U4 −1,74704(10 ) −4
u2 0 0 0 0 1 0
−2,47974(10−5 )
υ 2
0 0 0 −1 0 0
U5
−4
θ2 0 0 0 0 0 1 U6 −9,94058(10 )
Las tensiones normales provenientes de la flexión se calculan en los nodos 1 y 2 aplicando las Ecuacio-
nes (4.24) como
7 6 −5 2 −4 2
σx (x = 0) = ±0,5(10 ) (−2,47974)(10 ) − (−9,94058)(10 ) = ±495,2 lb/pulg
202 20
7 6 −5 2 −4 2
σx (x = L) = ±0,5(10 ) (2,47974)(10 ) + (2)(−9,94058)(10 ) = ±992,2 lb/pulg
202 20
y, la tensión axial calculada mediante σx = E = E(u2 − u1 )/L es
−1,74704(10−4 ) 2
σx = 107 = −87,35 lb/pulg
20
Luego, superponiendo los resultados, podemos indicar que la magnitud de tensión más grande
ocurre en el nodo 2, en el que la tensión normal axial compresiva se agrega a la porción compresiva de
la distribución de tensión normal proveniente de la flexión, para dar como resultado final
2
σxmáx = −1079,6 lb/pulg (compresión)
La evaluación de las tensiones internas en el elemento 2 es mucho más simple, puesto que no se
requiere efectuar una transformación de desplazamientos (los desplazamientos globales de solución ha-
llados, son al mismo tiempo desplazamientos de tipo local — referidos a su propio sistema coordenado);
y este breve cálculo se deja al lector como ejercicio complementario a éste ejemplo. >
y z
q z(x)
w1 w2
x
1 2
z x y1 y2
L
(a) Definición de sistema coordenado local (b) Desplazamientos nodales locales y solicitación
de rigidez de tal elemento y obtener la matriz de rigidez para el mismo, nosotros primero extendemos
el elemento marco de la sección anterior para incluir la flexión en dos planos mutuamente ortogonales,
y luego agregaremos la capacidad de respuesta torsional.
La Figura 4.15(a) muestra un elemento marco con un sistema coordenado local o propio del mismo,
en el que el eje x corresponde al eje longitudinal del miembro y se asume que pasa a travéz del
centroide de todas las secciones transversales a lo largo de la longitud. Los ejes y y z se asumen que
corresponden a los ejes principales de los momentos de inercia de área de la sección transversal [29]. Si
éste no es el caso, el tratamiento de la flexión simultánea en dos planos y la superposición de efectos de
comportamiento según los estados de solicitación al que está sometido este elemento no producirá los
resultados correctos que son esperados [43].
Para la flexión alrededor del eje z (i.e., el plano de flexión es el plano xy), la matriz de rigidez
de elemento está dada por la Ecuación (4.37). Para la flexión alrededor del eje y, el plano de flexión
es el plano xz, como se aprecia en la Figura 4.15(b), que muestra un elemento viga definido por los
nodos 1 y 2, sujeto a una carga distribuida qz (x) que suponemos actúa según la dirección z positiva.
Los desplazamientos nodales en la dirección z se denotan como w1 y w2 ; mientras que las rotaciones
nodales son θy1 y θy2 . Para este caso, es necesario agregar el subı́ndice del eje y a las rotaciones nodales
para identificar especı́ficamente el eje sobre el que las rotaciones son medidas. En este contexto, las
rotaciones correspondientes al plano de flexión xy se denotan θz1 y θz2 de aquı́ en adelante. También
es importante notar que, en la Figura 4.15(b), el eje y es perpendicular al plano de la página con el
sentido positivo ingresando hacia ella. Por consiguiente, las rotaciones mostradas son positivas sobre
el eje y por la regla de mano–derecha. Notando la diferencia en el sentido positivo de la rotación
relativa a los desplazamientos lineales, un desarrollo análogo a aquel usado para el elemento viga en
las Secciones 4.3 y 4.4 resulta en la matriz de rigidez de elemento para el plano xz, considerado plano
de flexión del elemento, como
12 −6L −12 −6L
E Iy −6L 4L2 6L 2L2
[ ke ]xz = 3 (4.49)
L −12 6L 12 6L
−6L 2L2 6L 4L2
Las únicas diferencias entre la matriz de rigidez en el plano de flexión xz y aquella correspondiente
para el plano xy, se ve que son los cambios de signo en los términos fuera de la diagonal–principal, y
el hecho que la rigidez caracterı́stica depende del momento de inercia de área Iy , que en este caso es
el pertinente por ser el eje y aquel alrededor del cual ocurre la curvatura del elemento.
Combinando en secuencia la matriz de rigidez de elemento en su comportamiento axial, la matriz
de rigidez flexionante en el plano xy, y la matriz de rigidez flexionante en el plano xz dada por la
Ecuación (4.49), las ecuaciones de equilibrio de elemento para un miembro estructural que soporta
4.8. EL ELEMENTO RECTILÍNEO GENERAL 113
flexión en dos planos, con rigidez axial adicionalmente, se escriben en forma matricial como
u1
fx1
u2 fx2
υ f
y
1
1
θz1 Mz1
[kbarra ] [ 0 ] [0]
υ2
fy2
[ ke ] = [ 0 ] [kviga ]xy [ 0 ]
= (4.50)
θz2
Mz2
[0] [ 0 ] [kviga ]xz w1 fz1
θ y1 My 1
w2 fz2
θ y2 My 2
qz L
f 2
qz
1 −qz L2
Mqz1
12
= qz L
(4.51)
fqz
2
2
Mqz2
qz L2
12
valores que se obtienen utilizando las relaciones integrales que involucran las funciones de interpola-
ción polinómicas yá especificadas anteriormente (polinomios de Hermite), con la única variante que
las funciones de interpolación asociadas a las rotaciones deben ser cambiadas de signo para tener con-
sistencia con la definición de estos desplazamientos angulares alrededor del eje y, tal como muestra la
Figura 4.15(b).
L
1
T1 G, J T2 Mx1 Mx2
2
x x
(a) Eje circular con solicitación torsional (b) Notación torsional de elemento finito
La adición de torsión al elemento rectilı́neo general es cumplida con referencia a la Figura 4.16(a)
que muestra un cilindro circular sujeto a solicitación torsional (al cual de forma común se lo denomina
eje) mediante momentos torsores o cuplas aplicadas en planos transversales al eje axial en sus secciones
ubicadas en los extremos de este miembro circular. Un elemento finito correspondiente se muestra en
la Figura 4.16(b), donde los nodos son los puntos 1 y 2 coincidentes con los centroides de las secciones
transversales extremas; el eje del cilindro es el eje x, y los momentos torsores son positivos acorde a la
regla de mano–derecha. Por la teorı́a elemental de mecánica de sólidos deformables, es bien conocido
114 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
que el ángulo de torsión por unidad de longitud de un cilindro circular elástico, de seccción transversal
uniforme, sujeto a un torque T está dado por
dθ T
φ= = (4.52)
dx GJ
donde J es el momento polar de inercia del área se sección transversal y G es el módulo de elasticidad
transversal del material, o módulo de Young. Como el ángulo de torsión por unidad de longitud es
constante, el ángulo total de rotación del elemento puede expresarse en términos de las rotaciones
nodales y el momento torsor aplicado como
TL
θx2 − θx1 = (4.53)
GJ
y, de aquı́ resulta
GJ
T = (θx2 − θx1 ) = kT (θx2 − θx1 ) (4.54)
L
donde kT = GJ/L es el coeficiente de rigidez de resorte torsional equivalente para el eje.
La comparación de la Ecuación (4.54) con la Ecuación (2.2) para un resorte linealmente elástico
y la consideración de la condición de equilibrio Mx1 + Mx2 = 0, permite directamente plantear las
ecuaciones de comportamiento mecánico de equilibrio de elemento:
G J 1 −1 θx1 Mx1
= (4.55)
L −1 1 θx2 Mx2
y la matriz de rigidez final para un elemento rectilı́neo general tri–dimensional se observa que es una
matriz simétrica de dimensión 12×12 compuesta de matrices individuales de rigidez que representan
comportamiento axial, comportamiento flexionante en dos planos, y comportamiento torsional.
El elemento rectilı́neo general puede utilizarse en los análisis de elemento finito de estructuras marco
tri–dimensionales. Como con la mayorı́a de los elementos finitos, es a menudo necesario transformar
4.9. COMENTARIOS FINALES 115
las matrices de elemento desde el sistema coordenado local hacia el sistema de coordenadas globales.
El procedimiento de transformación es bastante similar a aquel discutido para el elemento barra y los
elementos viga bi–dimensionales; excepto, por supuesto, que con un grado de complejidad algebráica
agregada, proveniente del tamaño de la matriz de rigidez y ciertos detalles requeridos para describir la
orientación espacial relativa del elemento. Estos detalles en el desarrollo de formulación completa del
comportamiento mecánico del elemento rectilı́neo general los dejamos a la curiosidad del lector. Un
libro que presenta el método dirécto de rigidez aplicado a elementos estructurales lineales es uno que
escribı́ hace tiempo atrás [20], el cual es una fuente de consulta documental muy recomendada.
Problemas propuestos
4.1. Se conectan dos elementos viga idénticos a un nodo común como es mostrado en la Figura P4.1.
Asumiendo que los desplazamientos nodales υi , θi son conocidos, usar la Ecuación (4.23) para
mostrar que la tensión normal interna σx es, en general, discontinua en el lı́mite común de ambos
elementos (i.e., en el nodo 2). Bajo que condicion(es) la tensión normal interna serı́a continua
entre elementos?.
1 2 3
Figura P4.1
4.2. Para el elemento viga cargado como muestra la Figura P4.2, construya los diagramas de esfuer-
zo cortante y momento flector internos. Cual es la importancia de estos diagramas respecto a
las Ecuaciones (4.10), (4.13), y la relación V = dM/dx de la teorı́a de la mecánica de sólidos
deformables básica ?.
116 CAPÍTULO 4. ELEMENTOS DE FLEXIÓN
F F
M2
1 2
M1 L
Figura P4.2
4.3. Para una viga uniformemente cargada como se muestra en la Figura P4.3, la teorı́a de mecánica
de sólidos dá la deflexión máxima como
5 q L4
υmáx = −
384 E Iz
en x = L/2. Tratar esta viga como un solo elemento finito y calcular la deformación máxima.
Cómo se comparan los valores hallados ?. Existe adecuada coincidencia según el criterio de com-
paración planteado ?.
y
q
x
L
Figura P4.3
4.4. El elemento viga mostrado en la Figura P4.4 se solicita mediante una carga linealmente variable
de intensidad máxima q0 . Usando el criterio del principio de equivalencia del trabajo mecánico,
determine las fuerzas y momentos nodales equivalentes.
y
q0
x
L
Figura P4.4
4.5. Use los resultados del Problema 4.4 para calcular la deflexión en el nodo 2 de la viga mostrada
en la Figura P4.5, si la estructura se trata como un solo elemento finito.
y
q0
1 E, Iz, L 2 x
Figura P4.5
Problemas propuestos 117
4.6. Para el elemento viga de la Figura P4.5, calcule la fuerza y el momento de reacción en el apoyo
empotrado (nodo 1). Calcule también la tensión normal máxima debida a la flexión, asumien-
do que la altura de la viga es 2h. Cómo se hace la comparación del valor de tensión obtenido
˙
con la tensión máxima proporcionada por la teorı́a básica de la mecánica de sólidos deformables?.
4.7. Repita el Problema 4.5 usando dos elementos de longitudes iguales. Para este problema, consi-
dere: E = 3×107 lb/pulg2 , Iz = 0,1 pulg4 , L = 10 pulg, q0 = 10 lb/pulg.
4.8. Considere la viga mostrada en la Figura P4.8. Cual es el número mı́nimo de elementos que pue-
den usarse para modelar este problema con adecuada exactitud? Construya el vector de carga
nodal global correspondiente a su respuesta.
P
q0
Figura P4.8
4.9. Cual es la justificación para escribir la Ecuación (4.26) en la forma de la Ecuación (4.27) ?.
4.10. — 4.11. — 4.12. — 4.13. — 4.14. — 4.15. Para cada viga mostrada en la figura asociada,
calcule la deflexión en los puntos nodales. El módulo de elasticidad es E = 107 lb/pulg2 y la
sección transversal de cada una de las vigas es como se muestra en cada figura. También calcule
la tensión normal flexionante máxima en cada caso. Use el método de elemento finito elaborando
para cada situación un modelo con el número mı́nimo de elementos.
10 lb/pulg
500 lb
10 pulg 30 pulg 2
p 0.5 pulg
u
l
g
18 pulg 6 pulg
40 mm
45 N/m 10 mm
0.5 pulg
30 mm
2 pulg
D2 = 1.0 pulg
D1 = 1.5 pulg 400 lb-pulg
0.25 pulg
2h h
x
Figura P4.16
4.17. Usar los resultados del Problema 4.16. para deducir el valor de la componente k11 de la matriz
de rigidez de elemento.
4.18. La matriz de rigidez completa en todos sus términos del elemento de sección transversal lineal-
mente variable del Problema 4.16. está dada por
243 156L −243 87L
Eth3
156L 56L
2
−156L 42L2
[ ke ] =
60L3 −243 −156L 243 −87L
87L 42L2 −87L 45L2
t =0.25 pulg 10 lb
10 lb
t =0.25 pulg
6 pulg 6 pulg
12 pulg
(a) (b)
Figura P4.18
d. Evaluar mediante ambos modelos la tensión normal máxima proveniente de la flexión y com-
parar los resultados.
4.19. Las seis ecuaciones de equilibrio estático de un elemento marco plano en el sistema coordenado
local (propio del mismo elemento) se expresan en forma matricial como [ ke ]{δe } = {fe }, con el
vector {δe } dado como en la Ecuación (4.45), [ ke ] por la Ecuación (4.46), y {fe } como el vector
de cargas nodales locales definido mediante
T
{fe } = f1x f1y M1 f2x f2y M2
Para un elemento orientado en un ángulo ψ arbitrario respecto al eje X del sistema global coor-
denado, convertir las ecuaciones de equilibrio hacia el sistema coordenado global y verificar el
cumplimiento de la Ecuación (4.48).
4.20. Use la Ecuación (4.47) para expresar la energı́a de deformación acumulada por un elemento
marco en términos de los desplazamientos globales. Aplicar el principio de la energı́a potencial
mı́nima para deducir la expresión matricial que describe el comportamiento mecánico de equi-
librio estático referida al sistema coordenado global. (Advertencia: Éste procedimiento puede
resultar muy tedioso, si es que usted no sistematiza su deducción).
4.21. La estructura marco bidimensional mostrada en la Figura P4.21 está compuesta de dos miembros
de acero estructural (E = 107 lb/pulg2 ) con dimensiones de sección rectangular de 2×4 pulg,
donde la dimensión de 2 pulg es perpendicular al plano de carga o plano que contiene la estruc-
tura. Todas las conexiones se tratan como uniones soldadas. Usando dos elementos marco–plano
y la numeración de nodos mostrada, determine
30 pulg
1200 lb
20 lb/pulg
3 2
1500 lb _ pulg
30 pulg
Figura P4.21
4.22. Repetir el Problema 4.21. para el caso en el que la conexión entre elementos en el nodo 2 es una
conexión articulada (que no transmite momentos flectores).
4.23. La estructura marco mostrada en la Figura P4.23 es la estructura soporte para un guinche lo-
calizado en el punto de aplicación de la carga W. Los apoyos en A y B son completamente fijos
(apoyos empotrados). Se sueldan otras conexiones a esta estructura. Asumiendo que la estruc-
tura va a ser modelada usando el número mı́nimo de elementos marco:
A B
Figura P4.23
4.24. Repetir el Problema 4.23. para la estructura marco plano mostrada en la Figura P4.24.
Figura P4.24
y y
1.5 m 500 N
z
6 cm
z x 3 cm
300 N
Figura P4.26
4.27. Repetir el Problema 4.26. para el caso en el cual las fuerzas concentradas aplicadas son reem-
plazadas por una carga uniformemente distribuida a lo largo de toda la longitud del miembro,
con componentes qy = 6 N/cm y qz = 4 N/cm actuando según las direcciones positivas y y z
del sistema coordenado, respectivamente.
El método de residuos ponderados
Capítulo 5
Hemos mencionado que el método de elemento finito es un procedimiento para obtener soluciones
aproximadas a problemas gobernados por ecuaciones diferenciales (tanto ordinarias como parciales).
En este Capı́tulo mostraremos como hallar la solución de una ecuación diferencial convirtiéndola en
una relación integral para la cual exigimos la minimización del error residual en el mismo dominio de
definición original. Con este objetivo presentaremos el método de residuos ponderados, y una variante
del mismo que da lugar al elemento finito de Galerkin.
La utilización de este tipo particular de elemento finito (formulado originalmente por Galerkin)
es relativamente sencilla en su formulación asociada a problemas que plantea la mecánica del medio
contı́nuo; por ello explotando esta singular ventaja mostramos aplicaciones de este método particular
a los elementos estructurales barra y viga, y también abordamos el problema de transmisón de calor
uni–dimensional, como una simple muestra de su potencialidad.
5.1. Introducción
Los capı́tulos 2, 3, y 4 introdujeron algunos de los conceptos básicos del método de elemento
finito en términos de los llamados elementos lineales. El resorte lineal elástico, el elemento barra, el
elemento viga, el elemento eje, y el elemento marco son elementos lineales porque sus propiedades
estructurales pueden describirse en función de una sola variable espacial que identifica la posición a lo
largo del eje longitudinal del elemento. Las relaciones fuerza–desplazamiento para los elementos lineales
son establecidas de manera dirécta, cuando estas relaciones se describen a partir de los conceptos
elementales de la mecánica de materiales sólidos deformables. Para extender el método de elemento
finito a situaciones más generales, particularmente las aplicaciones no–estructurales, se requieren de
técnicas matemáticas adicionales. En este capı́tulo, el método de errores residuales ponderados se
describe en general y se da énfasis al método de Galerkin de residuos ponderados [51] como una
herramienta para la formulación de modelos de elemento finito para esencialmente cualquier problema
de campo gobernado por una ecuación diferencial.
123
124 CAPÍTULO 5. EL MÉTODO DE RESIDUOS PONDERADOS
donde R(x) es el residuo (también llamado error residual). Notemos que el residuo también es una fun-
ción de los parámetros ci desconocidos. El método de residuos ponderados requiere que los parámetros
ci sean evaluados de modo que
Z b
wi (x)R(x)dx = 0 i = 1; n (5.4)
a
donde wi (x) representa n funciones de ponderación arbitrarias. Nosotros observamos que, en el proceso
de integración, la Ecuación (5.4) resulta en n ecuaciones algebráicas que pueden resolverse para los
n valores de los parámetros ci . En resumen, la Ecuación (5.4) expresa que la suma (integral) de los
errores residuales ponderados sobre el dominio del problema es cero (tratando que la aproximación
planteada coincida con la solución exacta). Debido a los requisitos impuestos sobre las funciones de
ensayo, la solución es al final exacta en los puntos extremos del dominio (las condiciones de borde
5.2. EL MÉTODO DE RESIDUOS PONDERADOS 125
lı́mite deben satisfacerse) pero, en general, en cualquier punto interior el error residual es distinto de
cero. Como se discute posteriormente en consecuencia, el MRP puede capturar la solución exácta bajo
ciertas condiciones, pero ésta ocurrencia es la excepción en lugar de la regla.
Existen muchas variaciones del MRP, y las técnicas varı́an principalmente en el modo en el que
los factores de ponderación son determinados o seleccionados. Las técnicas más comunes y usuales son
la colocación de punto, colocación del subdominio, mı́nimos cuadrados, y el método de Galerkin [51].
Como éste último es bastante simple en su uso y es prontamente adaptable al método de elemento
finito, nosotros discutiremos sólamente el método de Galerkin.
En el método de residuos ponderados de Galerkin, las funciones de ponderación son escogidas de
modo que sean idénticas a las funciones de ensayo o prueba; es decir,
produciendo de nuevo n ecuaciones algebráicas para la evaluación de los parámetros desconocidos. Los
ejemplos siguientes ilustran detalles del procedimiento.
Ejemplo 5.1.
Utilice el método de residuos ponderados de Galerkin para obtener una solución aproximada de la
ecuación diferencial
d2 y
− 10x2 = 5 0≤x≤1
dx2
con las condiciones de borde extremo y(0) = y(1) = 0.
> Solución
La presencia del término cuadrático en la ecuación diferencial sugiere que las funciones de prueba con
aspecto polinomial son apropiadas en este caso. Para condiciones de borde extremo lı́mite homogéneas
en x = a y x = b, la forma general
N (x) = (x − xa )p (x − xb )q
con p y q siendo números enteros positivos, automáticamente satisfacen las condiciones de borde y son
contı́nuas en el dominio xa ≤ x ≤ xb . Usando una función de prueba individual, la forma más simple
que satisface las condiciones de borde establecidas para el problema, resulta ser
N1 (x) = x(x − 1)
y ∗ (x) = c1 x(x − 1)
dy ∗
= c1 (2x − 1)
dx
d2 y ∗
= 2c1
dx2
126 CAPÍTULO 5. EL MÉTODO DE RESIDUOS PONDERADOS
respectivamente. (Vemos a estas alturas, que la solución de ensayo seleccionada no satisface la fı́sica del
problema, puesto que hemos obtenido una segunda derivada de valor constante. La ecuación diferencial
es tal que la segunda derivada debe ser una función cuadrática de x. No obstante estas observaciones,
continuaremos el ejemplo para ilustrar el procedimiento).
La substitución de la segunda derivada de y ∗ (x) en la ecuación diferencial proporciona el error
residual como
R(c1 , x) = 2c1 − 10x2 − 5
el cual claramente es no–nulo. La substitución en la Ecuación (5.6) nos dá
Z 1
x(x − 1)(2c1 − 10x2 − 5)dx = 0
0
la cual después de la integración arroja como resultado c1 = 4, de modo que la solución aproximada
se expresa como
y ∗ (x) = 4 x(x − 1)
Para este ejemplo relativamente simple, podemos comparar el resultado de solución aproximada
con la solución exacta, obtenida por doble integración de la ecuación diferencial planteada como sigue:
Z 2
10x3
Z
dy d y
= 2
dx = (10x2 + 5)dx = + 5x + C1
dx dx 3
10x3 5x4 5x2
Z Z
dy
y(x) = dx = + 5x + C1 dx = + + C1 x + C2
dx 3 6 2
Aplicando la condición de borde y(0) = 0 obtenemos C2 = 0; mientras que la condición y(1) = 0
proporciona
5 5
+ + C1 = 0
6 2
de donde C1 = −10/3. De aquı́, la solución exacta vendrá dada por
5 4 5 2 10
y(x) = x + x − x
6 2 3
x
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
- 0.2
- 0.4
y(x) - 0.6
- 0.8
Exacta
-1.0 Galerkin
- 1.2
Una comparación gráfica de las dos soluciones se bosqueja en la Figura 5.1, la cual muestra que
la solución aproximada está en acuerdo razonable con la solución exacta. Sin embargo, debemos notar
5.2. EL MÉTODO DE RESIDUOS PONDERADOS 127
que la solución aproximada de un solo término es simétrica sobre el intervalo de interés. Puede verse
que esto no es correcto, simplemente examinando la ecuación diferencial. La perturbación principal en
la ecuación diferencial (escrita en la forma d2 y/dx2 = 10x2 + 5) proviene del término cuadrático en x;
por consiguiente, es improbable que la solución sea simétrica. El ejemplo siguiente expande la solución
aproximada y muestra cómo el método converge hacia la solución exacta. >
Ejemplo 5.2.
Obtener una solución aproximada de dos–términos de clase Galerkin al problema planteado en el
Ejemplo 5.1, usando las funciones de prueba
> Solución
La solución aproximada con dos términos en este caso es
y ∗ (x) = c1 x(x − 1) + c2 x2 (x − 1)
y la segunda derivada es
d2 y ∗
= 2c1 + 2c2 (3x − 1)
dx2
Substituyendo en la ecuación diferencial, obtenemos el residuo siguiente:
Utilizando las funciones de prueba como funciones de ponderación, por el método de Galerkin, las
ecuaciones residuales resultan
Z 1
x(x − 1)[2c1 + 2c2 (3x − 1) − 10x2 − 5]dx = 0
0
Z 1
x2 (x − 1)[2c1 + 2c2 (3x − 1) − 10x2 − 5]dx = 0
0
19 5
c1 = c2 =
6 3
de modo que la solución aproximada de dos términos es
19 5 5 3 19
y ∗ (x) = x(x − 1) + x2 (x − 1) = x3 + x2 − x
6 3 3 2 6
Para comparación se dibujan: la solución exacta, y las soluciones de uno y dos términos, en una sola
gráfica que se muestra en la Figura 5.2. Las diferencias entre la solución exacta y la solución aproximada
Galerkin de dos términos, son escasamente discernibles porque la coincidencia es casi plena. >
128 CAPÍTULO 5. EL MÉTODO DE RESIDUOS PONDERADOS
x
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
- 0.2
- 0.4
y(x) - 0.6
- 0.8
Exacta
-1.0 1 Término
2 Términos
- 1.2
Ejemplo 5.3.
Use el método de residuos ponderados de Galerkin para obtener la aproximación de un–término a la
solución de la ecuación diferencial
d2 y
+ y = 4x 0≥x≥1
dx2
con las condiciones de borde: y(0) = 0, y(1) = 1.
> Solución
Aquı́ las condiciones de borde lı́mite no son homogéneas [y(1) 6= 0], de modo que se requiere una
modificación al método presentado. Al contrario del caso de condiciones de borde lı́mite homogéneas,
no es posible construir una solución de ensayo de la forma c1 N1 (x) que satisfaga las dos condiciones
lı́mite declaradas. En cambio, asumimos una solución de prueba como
y ∗ (x) = c1 N1 (x) + f (x)
donde N1 (x) satisface las condiciones de borde homogéneas, y f (x) se escoge de modo que satisfaga
las condiciones de borde no–homogéneas. (Notemos que, si ambas condiciones de borde fuesen no–
homogéneas, dos de tales funciones deberı́an ser incluidas). Una de las soluciones aproximadas serı́a
y ∗ (x) = c1 x(x − 1) + x
la cual satiaface y(0) = 0 y y(1) = 1, idénticamente como puede comprobarse.
La substitución en la ecuación diferencial proporciona el error residual
d2 y ∗
R(c1 , x) = + y ∗ − 4x = 2c1 + c1 x2 − c1 x + x − 4x = c1 x2 − c1 x + 2c1 − 3x
dx2
y, la integral residual ponderada llega a ser
Z 1 Z 1
N1 (x)R(c1 , x)dx = x(x − 1)(c1 x2 − c1 x + 2c1 − 3x)dx = 0
0 0
Aunque el manipuleo algebráico es tedioso, la integración es directa y proporciona c1 = 5/6. Por tanto,
la solución aproximada es
5 5 1
y ∗ (x) = x(x − 1) + x = x2 + x
6 6 6
5.2. EL MÉTODO DE RESIDUOS PONDERADOS 129
Como en el ejemplo anterior, tenemos el lujo de comparar la solución aproximada con la solución
exacta, la cual es
y(x) = 4 x − 3,565 sin x
1.2
Exacta
1.0
Galerkin
0.8
y(x) 0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x
La solución aproximada y la solución exacta se muestran en la Figura 5.3 para una comparación. De
nuevo, se observa coincidencia razonable de resultados; pero podrı́a mejorarse agregando una segunda
función de ensayo a la solución aproximada, como lo hicimos en el ejemplo anterior. >
Cómo uno sabe cuándo la solución del MRP es bastante exacta ?. Es decir, cómo nosotros determi-
namos si la solución está cerca de la solución exacta ?. Esta pregunta de convergencia debe plantearse
en todas las técnicas de solución aproximadas. Si nosotros no conocemos la solución exacta, y rara-
mente lo hacemos, debemos desarrollar algún criterio para determinar el grado de exactitud de nuestra
solución aproximada. En general, para el método de los residuos ponderados, el procedimiento es
continuar obteniendo soluciones mientras se incremente el número de funciones de ensayo y notar el
comportamiento de la solución al introducir estas variaciones. Si la solución cambia muy poco cuando
aumentamos el número de funciones de ensayo, podemos decir que la solución converge. Ahora, si la
solución converge a la solución correcta o verdadera es todavı́a otra pregunta a ser respondida. El deta-
lle de estas cuestiones está más allá del alcance de este libro, un gran cuerpo de la matemática teórica
está dirigida a las preguntas sobre convergencia y si la convergencia es hacia la solución correcta. En el
contexto de este trabajo, nosotros asumimos que una solución que converge, lo hace hacia la solución
correcta. Ciertas pruebas de carácter externo al procedimiento de la solución, pueden hacerse para
determinar la “racionalidad” de una solución numérica en el caso de problemas fı́sicos. Estas pruebas
incluyen el equilibrio estático y dinámico, el equilibrio de energı́a, el balance de calor y flujo fluido, y
muchos otros que se discutirán en los capı́tulos siguientes.
Seguramente usted tiene la impresión que en los ejemplos anteriores se aplicó una metodologı́a de
prueba de ensayo–error; sin embargo ello está muy alejado de la realidad. El encontrar funciones de
prueba para aplicación del método de residuos ponderados no es aleatorio como podrı́a suponerse,
porque si se usara ésta técnica se deberı́a tener una gran experiencia en el conocimiento del análisis
funcional. En los ejemplos anteriores, nosotros usamos funciones de prueba “preparadas” para satisfacer
automáticamente las condiciones de borde lı́mite, pero las mismas no se basaron en un procedimiento
sistemático. Mientras que absolutamente nada está equivocado con este procedimiento, presentamos
ahora un procedimiento basado en funciones de ensayo polinómicas, que dan un método sistemático
para aumentar el número de funciones de ensayo y, de ello, una ayuda en la examinación del problema
130 CAPÍTULO 5. EL MÉTODO DE RESIDUOS PONDERADOS
y ∗ (x) = c0 + c1 x + c2 x2 + c3 x3 + · · ·
> Solución
Para un primer ensayo, tomamos sólo la forma cuadrática
y ∗ (x) = c0 + c1 x + c2 x2
y ∗ (0) = 0 = c0
y ∗ (1) = 0 = c1 + c2
y ∗ (x) = c1 x + c2 x2 = c1 x − c1 x2 = c1 x(1 − x)
y ∗ (x) = c0 + c1 x + c2 x2 + c3 x3
y ∗ (0) = 0 = c0
y ∗ (1) = 0 = c1 + c2 + c3
y hemos obtenido dos funciones de ensayo, cada una que satisface idénticamente las condiciones de
borde lı́mite. La determinación de las constantes para la solución de dos–términos se deja como ejercicio
de fin de capı́tulo. En cambio, nosotros agregaremos el término con exponente 4, y examinaremos la
solución de ensayo.
y ∗ (x) = c0 + c1 x + c2 x2 + c3 x3 + c4 x4
y las condiciones de borde dan
c0 = 0
c1 + c2 + c3 + c4 = 0
Usamos la relación de restricción para eliminar (arbitrariamente) el coeficiente c4 [ es decir que
tomamos c4 = −(c1 + c2 + c3 ) ], para obtener
Si ponemos la expresión residual igual a cero, e igualamos los coeficientes de las potencias de x,
hallamos que el residuo es exactamente nulo si
10 5 5
c1 = − c2 = c3 = 0 c4 =
3 2 6
5 4 5 2 10
de modo que: y ∗ (x) = x + x − x ; y hemos obtenido la solución exacta. >
6 2 3
d2 y
+ f (x) = 0 a≤x≤b (5.7)
dx2
sujeta a las condiciones de borde
y(a) = ya y(b) = yb (5.7a)
El dominio del problema es dividido en M “elementos” [véase la Figura 5.3(a)] los que en conjunto
contienen M + 1 valores xi de la variable independiente, siendo que x1 = xa y xM +1 = xb para asegurar
la inclusión de los lı́mites globales. Una solución aproximada es asumida en la forma
M +1
X
y ∗ (x) = yi ni (x) (5.8)
i=1
n1(x)
1
x1 x2 x3 x4 x5
n2(x)
1
x1 x2 x3 x4 x5
n3(x)
1
x1 x2 x3 x4 x5
n4(x)
1
1 2 3 M
x1 x2 x3 xM xM+1
( xa ) ( xb ) x1 x2 x3 x4 x5
del dominio global del problema. Especı́ficamente, una función de ensayo ni (x) es no–nula en el in-
tervalo xi−1 < x < xi+1 , y para facilidad de ilustración, usamos funciones lineales definidas como
sigue:
x − xi−1
ni (x) = xi−1 ≤ x ≤ xi
xi − xi−1
xi+1 − x (5.9)
ni (x) = xi ≤ x ≤ xi+1
xi+1 − xi
ni (x) = 0 x < xi−1 x > xi+1
Claramente, en este caso, las funciones de ensayo son simplemente funciones de interpolación abso-
lutamente lineales, tales que el valor de la solución y ∗ (x) en xi < x < xi+1 es una combinación lineal
de los valores “nodales” adyascentes yi y yi+1 . Las primeras cuatro funciones de ensayo son como las
mostradas en la Figura 5.3(b); y observamos que, en el intervalo x2 < x < x3 por ejemplo, la solución
aproximada como es dada por la Ecuación (5.8) es
x3 − x x − x2
y ∗ (x) = y2 n2 (x) + y3 n3 (x) = y2 + y3 (5.10)
x3 − x2 x3 − x2
(Las funciones de ensayo usadas aquı́ son lineales, pero se pueden usar funciones de orden–superior
como mostraremos posteriormente cuando apliquemos la técnica al análisis del elemento viga).
La substitución de la solución asumida descrita por la Ecuación (5.8) en la Ecuación (5.7) gober-
nante del problema, proporciona el error residual como
M +1 2 ∗ M +1
d2
X d y X
R(x, yi ) = + f (x) = {yi ni (x)} + f (x) (5.11)
i=1
dx2 i=1
dx2
al cual aplicamos el método de residuos ponderados de Galerkin usando como funciones de ponderación
las funciones de interpolación aquı́ planteadas, para obtener
Zxb Zxb M +1
X d2
nj (x)R(x, yi )dx = nj (x) {yi ni (x)} + f (x) dx = 0 j = 1; M + 1 (5.12)
i=1
dx2
xa xa
5.3. EL MÉTODO DE ELEMENTO FINITO GALERKIN 133
En vista de la Ecuación (5.9) y la Figura 5.3(b), observamos que, en cualquier intervalo arbitrario
xj ≤ x ≤ xj+1 , sólo dos de las funciones de ensayo son no–nulas. Tomando esta observación en cuenta,
la ecuación anterior puede expresarse como
xZj+1
d2
nj (x) {y n
j j (x) + y n
j+1 j+1 (x)} + f (x) dx = 0 j = 1; M + 1 (5.13)
dx2
xj
[ K ]{y} = {F } (5.14)
d2 y
+ f (x) = 0 xj ≤ x ≤ xj+1 (5.15)
dx2
donde xj y xj+1 están contenidos en (a, b) y definen los nodos de un elemento finito. Las apropiadas
condiciones de borde aplicables a la Ecuación (5.15) son
y estos son los valores incógnita de la solución en los puntos extremos del subdominio. A continuación,
proponemos una solución aproximada de la forma
donde el super–ı́ndice (e) indica que la solución es para el elemento finito planteado anteriormente, y
las funciones de interpolación N1 y N2 son definidas ahora como
xj+1 − x
N1 (x) = xj ≤ x ≤ xj+1 (5.17a)
xj+1 − xj
x − xj
N2 (x) = xj ≤ x ≤ xj+1 (5.17b)
xj+1 − xj
Note usted la relación entre las funciones de interpolación definidas en las Ecuaciones (5.17) y las
funciones de ensayo en la Ecuación (5.9). Las funciones de interpolación corresponden a las porciones
superpuestas de las funciones de ensayo aplicables en el dominio de un solo elemento. También note
que las funciones de interpolación satisfacen las condiciones
N1 (x = xj ) = 1 N1 (x = xj+1 ) = 0
(5.18)
N2 (x = xj ) = 0 N2 (x = xj+1 ) = 1
134 CAPÍTULO 5. EL MÉTODO DE RESIDUOS PONDERADOS
de modo que las condiciones de borde (nodales) especificadas en la Ecuación (5.15a), se satisfacen
idénticamente. La substitución de la solución asumida en la Ecuación (5.15) dá el residuo como
d2 y (e) d2
R(e) (x, yj , yj+1 ) = + f (x) = [yj N1 (x) + yj+1 N2 (x)] + f (x) 6= 0 (5.19)
dx2 dx2
donde el super–ı́ndice es usado nuevamente para indicar que el residuo calculado es para el elemento.
Aplicando el criterio de los residuos ponderados de Galerkin, resulta en
xZj+1 xZj+1
d2 y (e)
(e)
Ni (x)R (x, yj , yj+2 )dx = Ni (x) + f (x) dx = 0 i = 1, 2 (5.20)
dx2
xj xj
o también
xZj+1 xZj+1
d2 y (e)
Ni (x) dx + Ni (x)f (x)dx = 0 i = 1, 2 (5.21)
dx2
xj xj
x xZj+1 xZj+1
dy (e) j+1 dNi dy (e)
Ni (x) − dx + Ni (x)f (x)dx = 0 i = 1, 2 (5.22)
dx xj dx dx
xj xj
Notemos que en el proceso de obtención de las Ecuaciones (5.23) hemos hecho uso de las propiedades
explı́citamente indicadas en la Ecuación (5.18) en la evaluación de las funciones de interpolación en
los nodos del elemento.
La integración por partes de uno de los términos de la Ecuación (5.21) trae tres beneficios desta-
cables
1. El orden más alto de las derivadas que aparecen en las ecuaciones de elemento ha sido reducido en
una unidad.
2. Como se observará explcitamente, la matriz de rigidez fue hecha simétrica. Si nosotros no integrára-
mos por partes, una de las funciones de ensayo en cada ecuación serı́a diferenciada dos veces y la
otra función de ensayo no se diferenciarı́a en absoluto.
3. La integración por partes introduce las condiciones de borde gradiente (de la pendiente de la
función) en los nodos del elemento. La importancia fı́sica de las condiciones lı́mite de gradiente se
pondrán claras en las aplicaciones fı́sicas subsecuentes.
5.3. EL MÉTODO DE ELEMENTO FINITO GALERKIN 135
Zx2 Zx2
dN1 dN1 dN2 dy (e)
y1 + y2 dx = N1 (x)f (x)dx + (5.24a)
dx dx dx dx x1
x1 x1
Zx2 Zx2
dN2 dN1 dN2 dy (e)
y1 + y2 dx = N2 (x)f (x)dx + (5.24b)
dx dx dx dx x2
x1 x1
Zx2
dNi dNj
kij = i, j = 1, 2 (5.26)
dx dx
x1
y las fuerzas nodales de elemento están dadas por los lados derechos de las Ecuaciones (5.24).
3 4
x3 x4 x5
y3(4) = y3(4)
dy (3) dy (4)
6=
dx x4 dx x4
Si la solución de elemento finito fuera la solución exacta, las primeras derivadas de cada elemento
indicadas en la Ecuación (5.27) deberı́an ser iguales y el valor de la expresión deberı́a ser cero. Sin
embargo, las soluciones de elemento finito raramente son exactas, por lo que éstos términos no son, en
general, cero. No obstante, en el procedimiento de ensamble, se asume que en todos los nodos interiores,
las condiciones de gradiente que posee la variable involucrada aparecen como iguales y opuestas en los
elementos adyacentes y ası́ la cancelación se efectúa a menos que una influencia externa actúe en el
nodo. Sin embargo, en los nodos de borde lı́mite globales, los términos gradiente pueden especificar
condiciones lı́mite o representar “reacciones” obtenidas mediante la fase de solución. De hecho, una
técnica muy poderosa para obtener exactitud de las soluciones de elemento finito es examinar la
magnitud de las discontinuidades del gradiente en los nodos o, más generalmente, en fronteras de
inter–elementos.
136 CAPÍTULO 5. EL MÉTODO DE RESIDUOS PONDERADOS
> Solución
Primero, notemos que la ecuación diferencial es equivalente a
d dy
x − 4x = 0
dx dx
la cual, luego de dos integraciones directas y la aplicación de las condiciones de borde, tiene la solución
exacta
3
y(x) = x2 − ln x − 1
ln 2
Para la solución de elemento finito, la aproximación más simple es usar un elemento con dos nodos
extremos para el cual la solución de elemento se asume como
x2 − x x − x1
y(x) = N1 (x) y1 + N2 (x) y2 = y1 + y2
x2 − x1 x2 − x1
donde y1 y y2 son los valores nodales. La ecuación residual para el elemento es
Zx2
d dy
Ni x − 4x dx = 0 i = 1, 2
dx dx
x1
que luego de ser integrada por partes en su primer término adopta la forma
x Zx2 Zx2
dy 2 dNi dy
Ni x − x dx − 4 x Ni dx = 0 i = 1, 2
dx x1 dx dx
x1 x1
Expandiendo las dos ecuaciones representadas por el último resultado después de la substitución
de las funciones de interpolación y las primeras derivadas de las mismas, resulta
Zx2 Zx2
1 dy x2 − x
2
x (y1 − y2 ) dx = −x1 −4 x dx
(x2 − x1 ) dx x1
x2 − x1
x1 x1
Zx2 Zx2
1 dy x − x1
x (y2 − y1 ) dx = x2 −4 x dx
(x2 − x1 )2 dx x2
x2 − x1
x1 x1
5.3. EL MÉTODO DE ELEMENTO FINITO GALERKIN 137
x1 = 1,5 x2 = 2 k = 3,5
Z 2
2−x
F1(2) = −4 x dx = −1,6666 . . .
1,5 2 − 1,5
Z 2
(2)
x − 1,5
F2 = −4 x dx = −1,8333 . . .
1,5 2 − 1,5
Las ecuaciones de elemento son entonces
dy
(1) −1,1667 + dx
2,5 −2,5 y1
x1
=
−2,5 2,5 y2(1)
dy
−1,3333 + 1,5
dx
x2
dy
(2) −1,1667 − 1,5 dx
3,5 −3,5 y1
x
= 2
−3,5 3,5 (2)
y2 −1,8333 + 2 dy
dx
x3
Denotando los valores nodales del sistema como Y1 , Y2 , Y3 en x = 1 – 1,5 – 2 respectivamente; las
ecuaciones ensambladas del sistema son
−1,1667 + dy
2,5 −2,5 0 Y1 dx
x1
−2,5 6 −3,5 Y2 = −3
0 −3,5 3,5 Y3 −1,8333 + 2 dy
dx
x3
Aunque los detalles de cálculo numérico se dejan como ejercicio para el lector al final del capı́tulo,
aquı́ presentamos una solución de cuatro–elementos para este ejemplo (nuevamente usando nodos
igualmente espaciados en xi = 1,0 – 1,25 – 1,5 – 1,75 – 2,0) cuyo desarrollo se traduce en el siguiente
sistema de ecuaciones globales
−0,5417 − dy
4,5 −4,5 0 0 0 Y1
dx
x1
−4,5 10 −5,5 0 0 Y
−1,25
2
0
−5,5 12 −6,5 0 Y3 =
−1,5
0 0 −6,5 14 −7,5 Y4 −1,75
0 0 0 −7,5 7,5 Y5 −0,9583 + 2 dy
dx
x5
Aplicando las condiciones de borde Y1 = Y5 = 0 y resolviendo el sistema 3×3, se obtienen los siguientes
resultados
Y2 = −0,4026 Y3 = −0,5047 Y4 = −0,3603
dy dy
= −2,350 = 1,831
dx x1 dx x5
x
0 1.25 1.5 1.75 2.0 2.5
0
- 0.1
- 0.2
y(x) - 0.3
- 0.4
Exacta
-0.5 2 Elementos
4 Elementos
- 0.6
Para efectuar la comparación pertinente, se muestran las gráficas obtenidas de las soluciones de dos
y cuatro–elementos juntamente con la solución exacta, en la Figura 5.6. La solución de dos–elementos
se ve que es una aproximación grosera excepto en los nodos del elemento y la discontinuidad de la
derivada se aprecia ser significante. La solución de cuatro–elementos tiene los valores calculados de
y(x) en los nodos, que son casi idénticos a la solución exacta. Con cuatro–elementos, las magnitudes de
las discontinuidades de las primeras derivadas en los nodos son reducidas, pero todavı́a aparentemente
muy aproximadas a los valores exactos. >
El ejemplo anterior claramente muestra la convergencia de la solución aproximada hacia la solución
exacta a medida que el número de elementos utilizado en el modelo de elemento finito de Galerkin se
incrementa. Nótese además que la solución aproximada se podrı́a optimizar si se utilizaran funciones
de prueba de orden–superior.
dσx d d2 u(x)
= (Ex ) = E =0 (5.28)
dx dx dx2
donde hemos asumido el módulo de elasticidad constante, y hemos expresado la deformación unitaria
interna en términos del campo de desplazamientos internos. Denotando con L la longitud del elemento,
el campo de desplazamientos es discretizado acorde con la Ecuación (2.16) como
x x
u(x) = N1 (x) u1 + N2 (x) u2 = 1 − u1 + u2 (5.29)
L L
Y, puesto que el dominio de interés es el volumen del elemento, la ecuación residual de Galerkin llega
a ser
ZZZ 2 ZL 2
d u) Ed u
Ni (x) E 2 dV = Ni Adx = 0 i = 1, 2 (5.30)
dx dx2
V 0
donde dV = A dx, siendo A el área de la sección transversal del elemento. Integrando por partes y
re–ordenando, obtenemos
ZL L
dNi du du
AE dx = Ni A E (5.31)
dx dx dx 0
0
ZL
dN1 d du
AE (N1 u1 + N2 u2 )dx = −A E = −A E |x=0 = −A σx |x=0 (5.32a)
dx dx dx x=0
0
ZL
dN2 d du
AE (N1 u1 + N2 u2 )dx = A E = A E |x=L = A σx |x=L (5.32b)
dx dx dx x=L
0
Interpretando los términos que aparecen en el lado derecho de las anteriores ecuaciones, notamos que,
para el elemento barra las condiciones de borde del gradiente simplemente representan las fuerzas
externas aplicadas en los puntos nodales extremos del elemento, ya que f = σ A.
Las Ecuaciones (5.32) fácilmente son combinadas para ser presentadas en forma matricial como
donde los coeficientes de la matriz son integrados independientemente (la integral de una matriz es
igual a la integral de los términos que la definen).
140 CAPÍTULO 5. EL MÉTODO DE RESIDUOS PONDERADOS
Efectuando primeramente las derivadas de las funciones de interpolación y luego integrando los
coeficientes de la matriz, como puede comprobarse se obtiene la siguiente relación
A E 1 −1 u1 f1
= (5.34)
L −1 1 u2 f2
la cual como puede verificarse es el mismo resultado que fué obtenido en el Capı́tulo 2 para el elemento
barra (véase la Ecuación (2.32) y compare). Esto simplemente ilustra la equivalencia del método de
Galerkin y los métodos de equilibrio y energı́a (teorema de Castigliano) usados anteriormente para la
fo