Departamento de Mecnica Estructural y Construcciones Industriales - ETS Ingenieros Industriales Madrid a

Teor General del a Mtodo de los Elementos Finitos e

Francisco Beltrn a

Notas de Clase / Curso de Doctorado 1998-99

Presentacin o

Estas notas se han concebido como material de apoyo didctico dentro del curso de a doctorado \Teor General del Mtodo de los Elementos Finitos", que imparte el Depara e tamento de Mecnica Estructural y Construcciones Industriales de la ETS de Ingenieros a Industriales de Madrid. Se pretende dar al alumno la posibilidad de contrastar con ellas sus apuntes de clase y, de esta manera, ayudarle a comprender mejor las ideas transmitidas por el profesor. De acuerdo con los objetivos del curso de doctorado, se proporciona una panormica a general de los aspectos del mtodo de los elementos nitos (MEF) necesarios para iniciar al e alumno en su aplicacin industrial prctica. En este sentido, el documento puede resultar o a util tambin para aquellos que, al margen del curso, busquen una formacin bsica que e o a les permita utilizar programas basados en el MEF conociendo las l neas generales de la tecnolog numrica y sus limitaciones. a e No se trata de remplazar los muchos libros de texto que, desde diferentes opticas, abordan el MEF. Por el contrario, la idea ha sido componer un resumen introductorio, escrito en un lenguaje asequible, que sirva de punto de partida para la consulta de esos libros. As para facilitar esta labor, en las pginas nales se incluye una lista de referencias , a bibliogrcas donde el alumno interesado puede ampliar los conceptos expuestos. a

Collado-Villalba, septiembre de 1999

Para Alba y Hctor e

II

 Indice General
1 Introduccin o 1.1 Perspectiva histrica . . . . . . . . . . . . . . . . o 1.1.1 Or genes . . . . . . . . . . . . . . . . . . . 1.1.2 Evolucin . . . . . . . . . . . . . . . . . . o 1.1.3 Presente . . . . . . . . . . . . . . . . . . . 1.2 Panormica de aplicaciones industriales actuales a 1 1 1 3 4 5 8 8 9 11 12 14 16 16 17 19 19 19 20 22 24 25 27 27 27 30 32

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

2 Fundamentos matemticos (I) a 2.1 Introduccin . . . . . . . . . . . . . . . . . . . . . . . . . o 2.2 Mtodos de residuos ponderados . . . . . . . . . . . . . e 2.3 Mtodo de Galerkin . . . . . . . . . . . . . . . . . . . . e 2.4 Formulacin dbil . . . . . . . . . . . . . . . . . . . . . . o e 2.5 Analog mecnica: el principio de los trabajos virtuales a a 2.6 Soluciones dbiles aproximadas . . . . . . . . . . . . . . e 2.7 Teorema de Lax-Milgram . . . . . . . . . . . . . . . . . 2.8 Propiedad de aproximacin optima . . . . . . . . . . . . o 3 Fundamentos matemticos (II) a 3.1 Introduccin . . . . . . . . . . . . . . . . . . . o 3.2 Principios variacionales . . . . . . . . . . . . 3.3 Mtodo de Rayleigh-Ritz . . . . . . . . . . . . e 3.4 El problema elstico: notacin . . . . . . . . a o 3.5 Principios variacionales en Elasticidad . . . . 3.6 Ecuaciones de Euler y tipos de condiciones de

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . contorno

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

4 Programacin del MEF (I) o 4.1 Introduccin . . . . . . . . . . . . . . . . . . . . . o 4.2 La \receta" del MEF . . . . . . . . . . . . . . . . 4.3 Clculos por el MEF: datos y resultados . . . . . a 4.4 Flujo general en un programa de EF para clculo a

. . . . . . . . . . . . lineal

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

5 Tecnolog de elementos (I) a 35 5.1 Introduccin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 o 5.2 Formulacin convencional en desplazamientos . . . . . . . . . . . . . . . . . 36 o III

5.3 Ejemplo: elemento cuadriltero . . . . . . . . . . a 5.4 La transformacin isoparamtrica . . . . . . . . . o e 5.5 Integracin numrica . . . . . . . . . . . . . . . . o e 5.6 Algunas familias corrientes de funciones de forma

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

37 39 42 43 46 46 47 50 53 54 56 60 62 62 64 64 65 66 67 72 72 73 74 77 77 78 79 82 83 83 83 85 87 87 88 89 92 94 96 99

6 Tecnolog de elementos (II) a 6.1 Introduccin . . . . . . . . . . . . . . . . . . . . . . o 6.2 Formulacin C 1 en elementos viga . . . . . . . . . o 6.3 Formulacin C 1 en elementos placa . . . . . . . . . o 6.4 Dicultad de la formulacin C 1 en elementos placa o 6.5 Ejemplo: elemento triangular de 3 nodos . . . . . . 6.6 Formulacin C 0 en elementos viga . . . . . . . . . o 6.7 Formulacin C 0 en elementos placa . . . . . . . . . o 7 Tecnolog de elementos (III) a 7.1 Introduccin . . . . . . . . . . . . . . . . . . . . . o 7.2 Elementos no conformes . . . . . . . . . . . . . . 7.2.1 Elasticidad bidimensional . . . . . . . . . 7.2.2 Flexin de placas . . . . . . . . . . . . . . o 7.3 Bloqueo por deformacin isocrica . . . . . . . . o o 7.4 Bloqueo del elemento cuadriltero convencional . a 7.5 Soluciones heur sticas a los problemas de bloqueo 7.5.1 Integracin reducida . . . . . . . . . . . . o 7.5.2 Formulacin B . . . . . . . . . . . . . . . o 7.6 La prueba de la parcela . . . . . . . . . . . . . . 8 Tecnolog de elementos (IV) a 8.1 Introduccin . . . . . . . . . . . . . o 8.2 Principio variacional multicampo . 8.3 Discretizacin por el MEF . . . . . o 8.4 Las condiciones de Babuka-Brezzi s

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

9 Procedimientos de clculo (I) a 9.1 Introduccin . . . . . . . . . . . . . . o 9.2 Resolucin de sistemas de ecuaciones o 9.3 Eliminacin de Gauss . . . . . . . . o 9.4 Factorizaciones de Crout y Cholesky 9.4.1 Factorizacin de Crout . . . . o 9.4.2 Factorizacin de Cholesky . . o 9.5 Mtodo de resolucin frontal . . . . e o 9.6 Mtodos iterativos . . . . . . . . . . e 9.7 Iteracin de Jacobi y sus variantes . o 9.8 Mtodo del gradiente conjugado . . . e 9.9 Relajacin dinmica . . . . . . . . . o a

. . . . . lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

IV

10 Procedimientos de clculo (II) a 10.1 Introduccin . . . . . . . . . . . . . . . o 10.2 El problema de la elasticidad dinmica a 10.3 Discretizacin por el MEF . . . . . . . o 10.4 Procedimientos tipo Newmark . . . . . 10.5 Operador de Hilber-Hughes-Taylor .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

101 101 101 103 106 109 110 110 111 112 112 113 113 114 115 117 117 118 118 120 123 126 127 128 129 131 133 133 133 135 137 140 140 141 142 142 143 143 145

11 Estimacin a posteriori del error o 11.1 Introduccin . . . . . . . . . . . . . . . . . . . . . . . . . . . o 11.2 Conceptos bsicos . . . . . . . . . . . . . . . . . . . . . . . a 11.3 Indicadores y estimadores de error . . . . . . . . . . . . . . 11.3.1 Densidad de energ de deformacin (SED) . . . . . a o 11.3.2 Normas de residuos . . . . . . . . . . . . . . . . . . . 11.3.3 Estimadores del tipo Z 2 . . . . . . . . . . . . . . . . 11.3.4 Estimadores basados en diferencias entre funcionales 11.4 Procesos adaptables . . . . . . . . . . . . . . . . . . . . . . 12 Conceptos bsicos de la Mecnica de Slidos a a o 12.1 Introduccin . . . . . . . . . . . . . . . . . . . . . o 12.2 Fuentes de no linealidad en Mecnica de Slidos . a o 12.3 Tensor gradiente de deformacin . . . . . . . . . o 12.4 Teorema de descomposicin polar . . . . . . . . . o 12.5 Medidas de la deformacin . . . . . . . . . . . . . o 12.6 Tasa o velocidad de deformacin . . . . . . . . . o 12.7 Equilibrio y trabajo virtual . . . . . . . . . . . . 12.8 Medidas de tensin . . . . . . . . . . . . . . . . . o 12.9 Tasa o velocidad de cambio de la tensin . . . . . o 12.10Particin aditiva de la velocidad de deformacin o o 13 Procedimientos de clculo (III) a 13.1 Introduccin . . . . . . . . . . . . . . . . . . . . o 13.2 La \receta" del MEF en problemas no lineales . 13.3 El mtodo de Newton y sus variantes . . . . . . e 13.4 Problemas no estacionarios . . . . . . . . . . . 13.5 Procedimientos de integracin expl o cita . . . . . 13.5.1 Ideas generales . . . . . . . . . . . . . . 13.5.2 Flujo general . . . . . . . . . . . . . . . 13.5.3 Ventajas e inconvenientes . . . . . . . . 13.5.4 Medidas de calidad de la solucin . . . . o 13.5.5 Campo de aplicaciones industriales . . . 13.5.6 Tcnicas de aceleracin del clculo . . . e o a Bibliograf a

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

Cap tulo 1

Introduccin o
1.1 Perspectiva histrica o

Los mtodos de elementos nitos constituyen hoy en d el procedimiento habitual de e a clculo en Mecnica Estructural y Mecnica de Slidos en general. Su uso est tambin a a a o a e muy extendido en la resolucin de problemas de Transferencia de Calor, y empieza a cobrar o importancia en otras areas, como la Mecnica de Fluidos o el Electromagnetismo. a El conocimiento de estas tcnicas numricas resulta actualmente casi imprescindible e e para aquellos que se desenvuelven en el ambito de la Ingenier Civil y la Ingenier a a Mecnica, ya que la mayor parte de los anlisis de tensiones que se llevan a cabo en la a a industria estn basados en ellas. a A pesar de su gran difusin actual, los procedimientos de elementos nitos tal y como o los entendemos hoy en d son relativamente modernos. Su nacimiento y desarrollo es a una consecuencia de la disponibilidad de herramientas electrnicas de clculo cada vez o a ms potentes. Puede decirse, por tanto, que estas tcnicas son un resultado ms de la a e a revolucin informtica de nales del siglo XX. o a

1.1.1

Or genes

La rese~a histrica del mtodo de los elementos nitos (MEF) hay que iniciarla en la n o e dcada de los cincuenta, cuando el recin nacido ordenador digital hac por n posible e e a el clculo automtico de estructuras de barras sin recurrir a tediosos procedimientos de a a relajacin, como el de Cross o el de Kani. Se concibi entonces una nueva tcnica de o o e clculo, inabordable sin la ayuda del ordenador, que fue bautizada con el nombre de a \clculo matricial de estructuras", en reconocimiento del papel que desempe~a el algebra a n matricial en su formalismo matemtico. a Recordemos que el clculo matricial de estructuras1 se basa en la idea de dividir la a estructura en barras, dentro de las cuales se conoce la solucin exacta en funcin de ciertos o o coecientes que se hacen coincidir con los movimientos de los nodos extremos. Dichos coecientes se obtienen planteando el equilibrio de todos los nodos de la estructura y
1

Ver, por ejemplo, Alarcn, Alvarez y Gmez-Lera. Clculo Matricial de Estructuras, Revert, 1986. o o a e

resolviendo el sistema de ecuaciones que resulta. De esta manera, conocidos los coecientes o movimientos nodales, se desciende de nuevo al nivel local de cada barra y se obtiene la solucin de esfuerzos y movimientos en el conjunto de la estructura por agregacin de o o soluciones locales. El MEF naci como una generalizacin de esta idea bsica del clculo matricial. Alo o a a guien que trabajaba con sistemas estructurales complejos, que no se idealizaban bien mediante entramados de barras, pens que pod dividir su estructura en zonas o \eleo a mentos" ms complejos que una simple barra. Estos elementos estar conectados entre a an s tambin en nodos pero, a diferencia con el clculo matricial, dentro de ellos slo conoc e a o a la solucin de manera aproximada en funcin de los movimientos nodales. Al igual que o o en el clculo matricial, a partir de las soluciones locales se pod plantear el equilibrio de a a los nodos y obtener los movimientos nodales resolviendo un sistema de ecuaciones. Estos movimientos nodales den la solucin dentro de cada uno de los \elementos" en que se an o hab dividido la estructura y, por agregacin, la solucin en toda ella. Lo que ocurr es a o o a que, ahora, esta solucin no era la exacta, sino una aproximacin. o o La partida de nacimiento del MEF, en la que se publica por primera vez la idea anterior, est fechada en 1956. Se trata de un art a culo histrico aparecido en una revista relacionada o con la industria aeronutica2 . a As pues, el MEF naci en el ambito del clculo de estructuras y esto ha impregnado o a toda la terminolog asociada al mismo. En un principio se present como un procedia o miento de clculo ms, entre los muchos desarrollados por ingenieros ocupados en resolver a a problemas prcticos. Sin embargo, durante los a~os sesenta los investigadores descubriea n ron que la esencia de lo que hab sido una mera generalizacin del clculo matricial pod a o a a utilizarse, no slo para resolver problemas de clculo de estructuras, sino tambin probleo a e mas de campo en general, tales como problemas de elasticidad o de conduccin de calor. o La idea bsica segu siendo la misma: la divisin del dominio de clculo en peque~os a a o a n subdominios y la aproximacin en ellos de la variable de campo en funcin de su valor en o o puntos privilegiados llamados nodos. Aparec as el MEF moderno. a Por otro lado, tras el xito en las primeras aplicaciones, se comprob que a pesar de e o haber sido desarrollado con mentalidad prctica (ingenieril), el mtodo ten hondas ra a e a ces matemticas, en la l a nea del procedimiento de Ritz para obtener soluciones aproximadas de ecuaciones diferenciales3 o dentro de los llamados mtodos de residuos ponderados. e En su aplicacin a la elasticidad, el mtodo pod interpretarse tambin como una forma o e a e aproximada de resolver las condiciones de equilibrio derivadas del clsico principio de los a trabajos virtuales. Esta generalidad empez a atraer el inters de los matemticos, los cuales contribuyeron o e a decisivamente a explicar con rigor las bases del MEF. Sin embargo, debe hacerse notar que la contribucin de los matemticos al MEF ha ido siempre muy por detrs de las o a a aplicaciones prcticas. El MEF naci como una herramienta ingenieril y sus l a o neas bsicas a de desarrollo han estado siempre muy vinculadas a la presin de la industria por resolver o problemas. En muchas etapas de su evolucin se ha concebido y aplicado con xito una o e
M.J. Turner, R.W. Clough, H.C. Martin y L.J. Topp. Stiness and Deection Analysis of Complex Structures. Journal of Aeronautical Sciences, vol.23, 9. 1956. 3 El mtodo de Ritz data de 1909. e
2

determinada tcnica numrica antes de encontrar su justicacin matemtica rigurosa. De e e o a hecho, es sintomtico que el primer libro importante en que se analiza el MEF desde el a punto de vista matemtico se publicara en 19734 , cuando el mtodo llevaba al menos quince a e a~os emplendose en la industria y hab alcanzado una gran madurez en su aplicacin a n a a o problemas lineales.

1.1.2

Evolucin o

El MEF alcanza su mayor de edad hacia nales de los sesenta, con la aparicin de los a o primeros programas comerciales5 . En ese momento entra en franca competencia con el unico mtodo de clculo numrico disponible hasta entonces para problemas de campo: el e a e mtodo de diferencias nitas. En el ambito del anlisis de tensiones en slidos, el MEF se e a o impuso rpidamente, ya que est libre de las restricciones de tipo geomtrico que dicultan a a e el uso de los procedimientos clsicos de diferencias nitas en este campo. a Al nal de la dcada de los sesenta el MEF hab demostrado ya su potencia y su vere a satilidad, pero su empleo estaba todav muy restringido dentro de la industria aerospacial a y de defensa, debido al alt simo precio de los ordenadores de entonces. Empiezan a aparecer en aquel momento los llamados \centros de clculo", compa~ a nas que vend tiempo de ordenador a usuarios que carec de los \grandes" ordenadores an an necesarios para resolver problemas industriales. Los centros de clculo se organizaban a alrededor de un ordenador en el que se encontraban instalados, entre otros, los programas de elementos nitos. Los ingenieros del centro proporcionaban al usuario la documentacin o necesaria para preparar la entrada de datos a los programas e interpretar los resultados que se produc an. El usuario preparaba sus datos y los remit al centro de clculo, inicialmente a a mediante paquetes de tarjetas perforadas y, ms tarde, mediante cheros que se enviaban a a travs de una l e nea telefnica. Los datos se procesaban en el ordenador del centro de o clculo y los resultados le llegaban al usuario al cabo de unos d normalmente en forma a as, de tremendos listados de nmeros que tardaban tambin varios d en ser comprobados u e as e interpretados. Los centros de clculo tuvieron su auge en la dcada de los setenta. Contribuyeron de a e manera muy importante a la popularizacin del MEF en industrias como la del automvil, o o la nuclear y la de grandes obras civiles. Por otro lado, los centros de clculo universitarios a pusieron la infraestructura necesaria para el enorme esfuerzo investigador que se llev a o cabo en esta dcada. Si los a~os sesenta fueron la poca de los pioneros, los a~os setenta e n e n son los de los grandes desarrollos del MEF, tanto en tecnologa de elementos como en procedimientos de clculo y aumento de prestaciones. El nmero de publicaciones sobre a u el mtodo creci exponencialmente y el MEF se aplic progresivamente a problemas cada e o o vez ms complejos, como el clculo de transitorios o el estudio de respuestas no lineales. a a Puede decirse que al nal de la dcada el desarrollo de las tcnicas numricas casi se pone e e e por delante de la potencia de clculo que son capaces de proporcionar los ordenadores. a Los centros de clculo inician su declive con la aparicin de los llamados \mini" ora o denadores, a principios de los ochenta. Los avances tecnolgicos permitieron poner en el o
4 5

Strang y Fix. An Analysis of the Finite Element Method, Prentice-Hall, 1973. El primer programa comercial de elementos nitos aparece en 1968.

mercado mquinas comparables a aquellas de que disponan los centros de clculo, pero a a a precios mucho ms bajos y con unos costes de mantenimiento y explotacin muy inferiores. a o El avance se hizo vertiginoso hacia el nal de la dcada, con la aparicin de las primeras e o \estaciones de trabajo", ordenadores pensados para un solo usuario, con una potencia de clculo nada despreciable, dotadas de capacidades grcas y con un precio peque~o. Como a a n consecuencia, los ordenadores se trasladan desde los centros de clculo a las ocinas de los a ingenieros y stos ganan autonom para usar el MEF y experimentar con l. e a e Durante la dcada de los ochenta el desarrollo de las tcnicas de elementos nitos no e e fue tan espectacular como en los setenta, se empez a alcanzar un cierto grado de madurez. o El esfuerzo investigador puntero se concentr ms en estos a~os en aplicaciones dentro del o a n a mbito no lineal, las cuales pod empezar a ser utilizadas de manera rutinaria gracias a an los avances en la potencia de clculo. a Donde s hubo un avance importante fue en la popularizacin del MEF y en su facilidad o de uso, tanto por el abaratamiento espectacular de los ordenadores, como por las capacidades grcas que proporcionaban. En la dcada de los ochenta empiezan a comercializarse a e pre y post-procesadores grcos para los clculos de elementos nitos, siendo ste un paso a a e muy importante de cara a poder abordar de manera rutinaria y con un m nimo de garant a clculos tridimensionales con geometr complejas, como las que aparecen en el dise~o a as n mecnico. a

1.1.3

Presente

La decda de los noventa corresponde al momento que vivimos. Se caracteriza por un abaa ratamiento de los ordenadores impensable hace slo unos a~os. Desde el punto de vista de o n lo que hace falta para el clculo por elementos nitos, puede decirse que hoy en d resulta a a normalmente ms caro el programa de clculo que el ordenador que se necesita para ejea a cutarlo. Todo lo contrario de lo que suced hace apenas una dcada, cuando el vendedor a e de ordenadores (el \hardware") prcticamente regalaba los programas (el \software") al a hacer una venta. Adems, con mucho, los mayores gastos asociados a un anlisis por a a elementos nitos no son ya los correspondientes al anlisis mismo (amortizacin del ordea o nador y licencia de uso del programa) sino los de preparacin del modelo e interpretacin o o de resultados. El abaratamiento de ordenadores y programas ha contribuido a que la difusin de las o herramientas de elementos nitos sea tremenda. Cualquier ocina tcnica, por peque~a e n que sea, las tiene a su alcance. Hay que decir a este respecto que la difusin de las o herramientas no siempre se corresponde con la adecuada formacin para su uso. Hoy en o d resulta relativamente frecuente que se lleven a cabo clculos por personal que desconoce a a casi absolutamente los fundamentos del MEF y sus limitaciones y que, por tanto, es incapaz de evaluar la bondad de los resultados que est obteniendo. a Otro aspecto importante del momento actual es la integracin del clculo por elementos o a nitos con otras ramas de lo que se ha dado en llamar Ingenier Asistida por Ordenador a (\Computer Aided Engineering" - CAE). En la actualidad es normal la integracin del o clculo por elementos nitos (\Finite Element Analysis" - FEA) y el dibujo asistido por a ordenador (\Computer Aided Design" - CAD), con el objetivo, siempre, de reducir los 4

tiempos de proyecto o de puesta de producto en el mercado. Las tcnicas de clculo no lineal han alcanzado una madurez suciente como para e a poder ser empleadas por la industria de forma rutinaria. No tienen an la difusin alu o canzada por los mtodos de clculo lineal y requieren de ordenadores ms potentes, pero e a a se emplean ya ampliamente en campos tales como el estudio de la resistencia a impacto de veh culos (\crashworthiness"), el dise~o de procesos de conformado de piezas metlicas n a (forja, estampacin, extrusin, laminacin) y el proyecto de componentes elastomricos. o o o e El objetivo es tambin el mismo, reducir al mximo el nmero de pruebas con prototipos e a u reales para acortar los plazos de dise~o o de puesta en el mercado. n No se est viviendo una poca de grandes avances en cuanto a las tcnicas de clculo a e e a por el mtodo de los elementos nitos. Se sigue investigando, pero el MEF ha alcanzado e ya un grado de madurez que no se presta a progresos espectaculares como los vividos en las dcadas anteriores. Desde el punto de vista del que redacta estas notas son cuatro e las l neas de investigacin a lo largo de las cuales se est desarrollando el MEF en la o a actualidad: Adaptacin de algoritmos de clculo a las nuevas arquitecturas de ordenadores, con o a objeto de aumentar la velocidad de clculo y, por tanto, el tama~o mximo de los a n a problemas abordables. Desarrollo de medidas error, mallados autoadaptativos y elementos de altas prestaciones, con objeto de aumentar la precisin y abilidad de los resultados obtenidos o por usuarios inexpertos en entornos de clculo integrados con el CAD. a Desarrollo de nuevos elementos y tcnicas de solucin encaminados a aumentar la e o eciencia, robustez y abilidad de los clculos en el ambito no lineal. a Modelos numricos de leyes de comportamiento de materiales, sobre todo para la e prediccin del fallo y para la representacin del comportamiento de nuevos materiao o les.

1.2

Panormica de aplicaciones industriales actuales a

Hoy en d la aplicacin industrial mayoritaria del MEF es el clculo de tensiones en slidos a o a o y estructuras. En esta parcela prcticamente no se usa otro procedimiento numrico. Para a e problemas muy concretos, tales como los relacionados con dominios innitos (acstica, u suelos) o el estudio de fracturas, es posible que en un futuro el Mtodo de los Elementos e de Contorno (MEC) pueda desplazar al MEF, por ser intr nsecamente ms adecuado. Sin a embargo, el conocimiento y el uso del MEC, no ya en la industria, sino incluso dentro de los ambientes docentes, son m nimos. No parece, ni siquiera a medio plazo, que el MEC pueda jugar un papel signicativo en la prctica industrial6 . a
6 Frente al gran nmero de programas basados en el MEF que existen hoy en d en el mercado, u a varias decenas, el autor slo conoce dos programas comerciales basados en el MEC. Esto da idea de o la desproporcin actual entre el uso que hace la industria de una y otra tcnica numrica. o e e

Dentro del clculo de tensiones hay que distinguir entre dos tipos generales de aplia caciones: el clculo lineal y el no lineal. La gran mayor de los usuarios del MEF en a a la actualidad, en torno al 80%, realiza clculos lineales. Las tcnicas de clculo lineal a e a estn lo sucientemente maduras y probadas como para que puedan emplearse de modo a generalizado sin apenas incertidumbres en cuanto a los recursos necesarios para llegar al resultado7 . El clculo lineal de tensiones, tanto esttico como dinmico, se utiliza sobre todo en a a a la fase de dise~o o de proyecto, donde se busca hacer un uso eciente del material y, en n ocasiones, justicar el cumplimiento de una normativa o cdigo de buena prctica. Su uso o a est muy difundido en el proyecto de elementos mecnicos y estructuras complejas. Se a a utiliza mucho tambin en el estudio de vibraciones (p.ej. acstica o ingenier s e u a smica). Por otro lado, los clculos lineales por elementos nitos juegan un papel destacado en los a procesos de licenciamiento o certicacin de componentes en la industria nuclear, \oo shore" o aeronutica. a El clculo y la visualizacin de los resultados permite al ingeniero entender mejor el a o funcionamiento de sus dise~os y, en consecuencia, optimizarlos. En este sentido, el clculo n a lineal ha sustituido casi completamente a los ensayos y pruebas de prototipos en que se basaba buena parte del dise~o mecnico hace slo unas dcadas. No porque el clculo sea n a o e a ms barato, que muchas veces no lo es, sino porque es mucho ms rpido e interactivo. a a a Permite realizar muchas pruebas del tipo \>qu pasar si. . . ?" en poco tiempo, lo que e a facilita enormemente la compenetracin entre el proyectista y su dise~o. o n El clculo no lineal de tensiones comienza a tener un peso especco grande dentro de las a aplicaciones prcticas del MEF. La industria ha impulsado mucho la investigacin en esta a o l nea con el objetivo de que, a medio plazo, se puedan llegar a eliminar las incertidumbres que afectan hoy en d a los clculos no lineales. Aunque se ha avanzado bastante en la a a ultima dcada, todav existen areas en las que abordar un clculo no lineal tiene una e a a cierta componente de investigacin, ya que no se conocen a priori los recursos que sern o a necesarios para alcanzar el resultado. Esto, junto con los mayores requisitos de formacin o y de infraestructura informtica que se imponen al usuario, ha retrasado la difusin de los a o clculos no lineales. a Sin embargo, en determinados sectores industriales la no linealidad de los clculos a no puede evitarse, ya que es parte intr nseca del comportamiento que intenta simularse. Es el caso normalmente de la industria de defensa (bal stica terminal), la ingenier de a determinados procesos de fabricacin (conformado de metales y vidrio), la industria de o componentes elastomricos (juntas de goma, soportes de caucho-metal), las aplicaciones e geotcnicas o el estudio de la seguridad a impacto de veh e culos (\crashworthiness"). Es en estas areas donde se encuentra ms difundido el clculo no lineal de tensiones utilizando a a el MEF. En c rculos ms minoritarios, el clculo no lineal de tensiones se utiliza tambin en la a a e investigacin de causas de accidentes (ingenier forense) y en la obtencin de las cargas o a o ultimas o l mites resistentes de las estructuras. Esta clase de estudios contrasta con los de proyecto en el sentido de que se busca una descripcin lo ms ajustada posible del o a
7

El unico recurso signicativo hoy en d es el tiempo del calculista. a

comportamiento real, mientras que en los clculos de dimensionamiento de estructuras lo a que se busca es, simplemente, garantizar la seguridad. Fuera del ambito del clculo de tensiones, el MEF est muy difundido asimismo en el a a estudio de problemas de transferencia de calor, sobre todo en ingeniera mecnica (motores a y sistemas de refrigeracin). En este campo es tambin el MEF prcticamente la unica o e a herramienta numrica que se utiliza. Habr que distinguir igualmente entre clculos e a a lineales y no lineales y, en general, aplican los mismos comentarios hechos para el caso del clculo de tensiones. Quiz en los problemas de transferencia de calor son menos las a a incertidumbres cuando se aborda un clculo no lineal. a La herramienta de clculo tradicional dentro de la Mecnica de Fluidos ha sido el a a Mtodo de Diferencias Finitas (MDF). El MEF se encuentra menos difundido aqu pore que la representacin de la geometr no tiene tanta importancia como en Mecnica de o a a Slidos y porque en muchas de las aplicaciones de inters industrial los problemas tienen o e 8 . En aplicaciones que requieren unicamente clculos lineales, tales como carcter no lineal a a el estudio del ujo en medios porosos (aguas subterrneas), la difusin de contaminantes, a o o la propagacin de ondas de gravedad (oleaje), el MEF se encuentra bastante extendido; o mientras que los avances en las tcnicas de clculo no lineal han hecho que el MEF sea e a cada vez ms competitivo con el MDF en otras aplicaciones. Parece que las tcnicas de a e elementos nitos ganan progresivamente terreno, aunque sea disfrazadas de otros nombres, como \dominios nitos" o \volmenes nitos". u La utilizacin del MEF en problemas de Electromagnetismo a escala industrial es relao tivamente reciente, aunque existen ya numerosos programas comerciales disponibles. Las aplicaciones incluyen el proyecto de mquinas elctricas (motores, generadores, transfora e madores) y el estudio de componentes (aisladores, interruptores).

En los problemas clsicos de Mecnica de Fluidos la variable de campo es el vector de velocidades. La a a no linealidad aparece en las ecuaciones del movimiento por el trmino convectivo de la derivada total de e la velocidad.

Cap tulo 2

Fundamentos matemticos (I) a

2.1 Introduccin o

Desde el punto de vista matemtico, el Mtodo de los Elementos nitos (MEF) puede a e entenderse como un procedimiento para resolver numricamente problemas planteados e mediante ecuaciones diferenciales. En esto es similar a otros procedimientos, como el Mtodo de Diferencias Finitas (MDF) o el Mtodo de los Elementos de Contorno (MEC). e e La forma ms elegante de explicar los fundamentos matemticos del MEF parte de la a a teor de espacios normados y utiliza los conceptos del anlisis funcional1 . Este es el marco a a en el que hay que situarse si se quieren estudiar con rigor las bases del MEF e investigar sobre sus propiedades matemticas. a Sin embargo, desde el punto de vista pedaggico, iniciar el estudio del MEF situndose o a en este marco puramente matemtico tiene serios inconvenientes para los tcnicos. Entre a e ellos se pueden mencionar los siguientes: Para la mayor de los ingenieros la teor de espacios normados y el anlisis funcional a a a resultan demasiado generales, abstractas y alejadas de las aplicaciones prcticas en a las que estn interesados. a Se requiere un tiempo y esfuerzo considerable para manejar con soltura los conceptos del anlisis funcional. Se corre el riesgo de desanimar a los estudiantes que se acercan a por primera vez al MEF y de fomentar entre ellos la idea de que el mtodo es slo e o una gran teor matemtica, dif de entender, y sin relacin aparente con la forma a a cil o en que luego se resuelven los problemas reales. El iniciar el estudio del MEF desde un marco puramente matemtico no se corresa ponde con la evolucin histrica del mtodo, el cual fue concebido por ingenieros con o o e la idea de resolver problemas concretos, y en cuyo desarrollo las aplicaciones han ido siempre por delante de las justicaciones matemticas generales y elegantes. a Por otro lado, la mayor de los estudiantes no tienen necesidad de conocer con toda a generalidad las bases matemticas del MEF, ya que en su vida profesional se limitarn a a
1

Ver las referencias 6 y 7 de la bibliograf complementaria. a

a ser usuarios de programas comerciales de clculo. Estos alumnos slo necesitan tener a o claros los conceptos matemticos indispensables para hacer un uso prctico correcto de a a las tcnicas numricas, ya probadas, que incorporan los programas comerciales. e e Por las razones anteriores, y por limitaciones de espacio, se ha decidido buscar una solucin de compromiso para explicar los fundamentos matemticos del MEF, sin cargar o a el peso en la generalidad y la elegancia matemtica. Se ha elegido una aproximacin a o que muestre la base matemtica del MEF en un lenguaje lo menos oscuro posible para el a estudiante medio y poniendo nfasis en la l e nea ingenieril de desarrollo del mtodo. e El objetivo es transmitir ideas y conceptos, ms que desarrollos y formulaciones. Las a ideas permitirn luego al estudioso penetrar en aparatos matemticos ms complicados, a a a que lo unico que hacen es generalizar estas ideas y presentarlas de manera ms elegante. a En este cap tulo va a introducirse el MEF desde el punto de vista matemtico como a un caso particular del mtodo de residuos ponderados de Galerkin y, tambin, como un e e procedimiento de obtener una solucin aproximada a un problema planteado de forma o \dbil". e

2.2

Mtodos de residuos ponderados e

Sea el siguiente problema modelo2 . Encontrar una funcin u : ! <, con = [0; 1], tal que: o u;xx + f = 0 en =]0; 1[ (2.1)

siendo f : ! < una funcin lisa o \suave", y de manera que se cumplan las siguientes o condiciones de contorno: u;x (0) = h u(1) = g h2< g2< (2.2) (2.3)

La aproximacin clsica por diferencias nitas a la resolucin numrica de este probleo a o e ma es muy directa, no requiere apenas elaboracin. La \receta" podr ser la siguiente: o a 1. Tomar n puntos de =]0; 1[. 2. Aproximar el valor de u;xx en los n puntos en funcin del valor de u en esos puntos o y de las condiciones de contorno. 3. Particularizar la ecuacin de campo u;xx + f = 0 en cada uno de los n puntos o seleccionados. Se obtiene as un sistema de n ecuaciones lineales cuyas n incgnitas o son los valores de u en los puntos seleccionados.
2 Este problema corresponde, desde el punto vista mecnico, al equilibrio de una barra sometida a a solicitacin axil. o

4. La aproximacin a la solucin u en se obtiene interpolando entre los valores de u o o calculados para los puntos seleccionados. La solucin puede renarse aumentando o el nmero n de puntos. u Puede verse que el procedimiento de solucin anterior se basa en una aproximacin o o por puntos a la funcin incgnita. La aproximacin se extiende luego a todo el dominio o o o de clculo por interpolacin entre los valores obtenidos para esos puntos. a o Podr pensarse en utilizar directamente una funcin extendida a todo el dominio de a o clculo y denida a partir de ciertos coecientes incgnita. La receta ser entonces: a o a 1. Aproximar la solucin u en todo el dominio de clculo mediante, por ejemplo: o a uh = N0 +
n X

c j Nj

(2.4)

j=1

o donde cj 2 < son coecientes incgnita, N0 : ! < es una funcin conocida que o e cumple las condiciones de contorno 2.2 y 2.3 y Nj : ! < son funciones tambin 3 conocidas pero que cumplen las condiciones de contorno homogneas . e o 2. Determinar los coecientes incgnita cj de modo que la ecuacin diferencial 2.1 se o cumpla \en promedio" o de forma ponderada en todo el dominio . 3. Una forma de realizar lo anterior es tomar n funciones wk ; k = 1; : : : ; n y plantear las n ecuaciones:
Z
1

wk (uh + f ) dx = 0 ;xx

k = 1; : : : ; n

(2.5)

Este sistema proporciona un medio para determinar los coecientes incgnita cj ; j = o 1; : : : ; n. El procedimiento anterior es una muestra de lo que se conoce con el nombre de mtodo e de los residuos ponderados, ya que el segundo factor de la cantidad subintegral de 2.5 recibe el nombre de funcin residuo4 . o El mtodo de los residuos ponderados es muy general y muy potente. Se trata de e una tcnica que engloba procedimientos tan dispares a primera vista como el MEF, el e MEC o el ajuste por m nimos cuadrados. La generalidad estriba en que hay muchas posibilidades para elegir las funciones de aproximacin Nj y las funciones de ponderacin o o wk . El criterio de que la funcin N0 cumpla las dos condiciones de contorno 2.2 y 2.3 se o ha jado arbitrariamente y esta funcin podr seleccionarse tambin con otro criterio. o a e

Las Nj pueden interpretarse como una base de un cierto espacio de funciones. El signicado matemtico de la ecuacin 2.5 es el de hacer ortogonal la funcin residuo a un espacio a o o de funciones cuya base son las funciones wk .
4

10

2.3

Mtodo de Galerkin e

El mtodo de (Bubnov)Galerkin es el mtodo de residuos ponderados que corresponde a e e la formulacin ms clsica del MEF. Segn este mtodo, siguiendo el razonamiento de la o a a u e h seccin anterior, la aproximacin u se construye como: o o uh = gh + donde: o cj 2 < son coecientes incgnita,
n X

c j Nj

(2.6)

j=1

Nj : ! < son funciones que cumplen Nj (1) = 0.

gh : ! < es una funcin que cumple la condicin de contorno 2.3, gh (1) = g, o o

Adems, las funciones de ponderacin wk : ! <, se eligen iguales a las Nj : wj = a o Nj ; j = 1; : : : ; n. Ntese que 2.6 puede interpretarse como que uh pertenece a un espacio de funciones o h , el espacio de las funciones de aproximacin, construido por traslacin de otro espacio o o S V h cuya base son las funciones Nj . Entonces, el mtodo de Galerkin se caracteriza porque e utiliza como espacio de funciones de ponderacin el mismo espacio V h . o Las ecuaciones 2.5 de residuos ponderados quedan entonces:
Z |
1

Nk (uh + f) dx + ;xx
{z }

ponderacin en el dominio o

ponderacin en el contorno o

Nk (uh + h)jx=0 ;x
{z

=0

k = 1; : : : ; n

(2.7)

las cuales, integrando por partes y teniendo en cuenta que Nk (1) = 0, dan lugar a:
Z
1

Nk;x uh dx = ;x

1 0

Nk f dx + Nk hjx=0

k = 1; : : : ; n

(2.8)

y, sustituyendo 2.6:
Z
1 n X

Nk;x (

cj Nj;x ) dx =

j=1

Nk f dx + Nk hjx=0

h Nk;x g;x dx

k = 1; : : : ; n (2.9)

donde los segundos miembros son conocidos. Es decir, se obtiene un sistema de n ecuaciones lineales con n incgnitas que permite o determinar los coecientes cj de 2.6. El mtodo de Galerkin data de 1915 y constituye una primera manera de justicar el e MEF desde el punto de vista matemtico. La aportacin del MEF moderno al mtodo a o e de Galerkin consiste en una forma sistemtica, fcilmente automatizable, de construir las a a funciones Nj y gh a partir de funciones denidas localmente. A modo de ejemplo, la sistemtica anterior se aplica a continuacin a la resolucin a o o aproximada del problema modelo. El proceso ser el siguiente: a

11

1. Se divide el dominio ]0,1[ en subintervalos o \elementos": 0

@

i+1

2. Dentro de cada elemento se denen dos funciones locales a y b, por ejemplo:

@ @

a@

b 1
@ @ @

i+1

3. Las funciones Nj se denen combinando estas dos funciones elementales:

gh

@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @

Ni i

g 1

i+1

Ntese que Nj (1) = 0 8j y que los coecientes cj cobran el sentido del valor de uh en o los puntos o \nodos" que se han utilizado para denir los subdominios o \elementos". Una vez denidas las funciones elementales, las funciones de aproximacin quedan o denidas automticamente. a 4. Las integrales que aparecen en 2.9 se hacen elemento a elemento, lo cual facilita tambin la automatizacin. Los dominios de denicin de las funciones Nj son e o o locales. En este caso se extienden, como mucho, a dos elementos. o 5. Los coecientes cj , que denen la solucin aproximada, se obtienen resolviendo el sistema 2.9.

2.4

Formulacin dbil o e

El problema modelo planteado en forma de ecuacin diferencial 2.1 con condiciones de o contorno 2.2 y 2.3 recibe el nombre de forma \fuerte" del problema, porque impone las condiciones ms exigentes a la funcin que se trata de obtener, en cuanto al orden de a o derivabilidad y en cuanto al cumplimiento de la ecuacin y de las condiciones de contorno o punto a punto dentro del dominio de clculo. a Puede pensarse en plantear el problema de modo que se exija a la funcin incgnita un o o orden menor de derivabilidad y, en vez del cumplimiento punto a punto, un cumplimiento 12

en promedio de la ecuacin diferencial y de las condiciones de contorno. Por contraste con o lo anterior, esta ser una formulacin \dbil" del mismo problema. a o e Una formulacin dbil de nuestro problema modelo ser la siguiente. Se denen dos o e a espacios o conjuntos de funciones, S = fu j u 2 H 1 ; u(1) = gg y V = fw j w 2 H 1 ; w(1) = 0g

donde H 1 es el conjunto de funciones con derivada de cuadrado integrable en = [0; 1], esto es: H = fu : ! < j
1

(u;x )2 dx < 1g

A partir de estos dos espacios, el problema se formula como, dados f : ! <, g 2 < y h 2 < como antes, encontrar una funcin u 2 S tal que 8w 2 V se cumpla: o
Z
1 0

w;x u;x dx =

wf dx + w(0)h

(2.10)

A primera vista puede parecer que el problema planteado as tiene poco que ver con el planteado en 2.1 a 2.3. Sin embargo, cuando las funciones que intervienen son sucientemente lisas o \suaves", en el sentido de que estn denidas las derivadas necesarias, a resulta que la formulacin dbil es equivalente a la formulacin fuerte. En efecto: o e o 1. Si u es solucin de la forma fuerte, tambin lo es de la dbil, ya que se cumplir: o e e a
Z
1

w(uh + f) dx = 0 ;xx

8w 2 V

e integrando por partes: 0=

Z
1

w;x u;x dx

wf dx [wu;x ]1 0

8w 2 V

y como u;x (0) = h, por ser u solucin de la forma fuerte, y w(1) = 0, por ser o w 2 V, resulta:
Z
1 0

w;x u;x dx =

wf dx + w(0)h

8w 2 V

y como adems u(1) = g, por ser u solucin de la forma fuerte, se tiene que u 2 S, a o luego u es solucin de la forma dbil. o e 2. Si u es solucin de la forma dbil, tambin lo es de la fuerte, ya que se cumplir: o e e a
Z
1

w;x u;x dx =

wf dx + w(0)h 13

8w 2 V

Integrando por partes y teniendo en cuenta que w(1) = 0 8w 2 V, resulta: 0=

Z
1 0

w(u;xx + f) dx + w(0)[u;x (0) + h]

8w 2 V

(2.11)

Como u(1) = g, por ser u 2 S, bastar con probar que la relacin anterior implica a o necesariamente que: u;xx + f = 0 y u;x (0) + h = 0 En efecto, se puede tomar w = (u;xx + f ) x (x 1) 2 V y entonces, por 2.11: 0=
Z
1 0

en

x(1 x) (u;xx + f)2 dx

|
0

{z

}|

{z

u;xx + f = 0 en

(2.12)

Como la pertenencia de w a V no impone ningn requisito sobre el valor w(0), 2.12 u o implica tambin, por 2.11, que u;x (0) = h. Luego u es solucin de la forma fuerte. e Entonces, si la solucin de la forma fuerte es unica y la solucin de la forma dbil es o o e tambin unica, los dos resultados anteriores implican que ambas soluciones sern iguales e a y que, por tanto es indiferente plantear el problema de una u otra manera.

2.5

Analog mecnica: el principio de los trabajos virtuales a a

Para concretar las ideas anteriores sobre formas fuertes y dbiles, estudiemos el equilibrio e de una barra sometida a solicitacin axil. o p-

N + dN

dx La forma fuerte del problema se obtiene planteando directamente el equilibrio de la rebanada. Si N es el esfuerzo axil; u, el desplazamiento a lo largo de la directriz; E, el mdulo de elasticidad; A la seccin de la barra y p las fuerzas distribuidas por unidad de o o longitud, se tiene: 14

N = A = AE = AEu;x

dN = AEu;xx dx

dN + pdx = 0

(equilibrio de la rebanada)

u;xx +

p =0 AE

(ecuacin diferencial) o

Suponiendo una barra de longitud unidad, las condiciones de contorno podran ser: u;x (0) = h (axil impuesto)

u(1) = g (desplazamiento impuesto) La forma dbil del problema se obtiene planteando el equilibrio mediante el principio de e los trabajos virtuales. La aplicacin de dicho principio a este caso dice que en un campo o virtual de desplazamientos w denido a lo largo de la barra, el trabajo de las fuerzas exteriores (las fuerzas distribuidas p y el axil impuesto en uno de los extremos) es igual al trabajo de las fuerzas interiores (tensiones ). Campo \virtual" de desplazamientos signica que el campo debe de ser continuo a lo largo de la barra, para respetar la compatibilidad interna (no \romper" la barra), y debe anularse en los contornos (extremos) donde se han impuesto desplazamientos, para respetar la compatibilidad externa. En este caso, debe ser w(1) = 0. Ntese que precisamente estas o dos condiciones son las que se necesitan para que w 2 V, ya que la continuidad implica la pertenencia a H 1 . El trabajo de las fuerzas exteriores para los desplazamientos virtuales w es:
Z
1

wf dx + hAEw(0)
0

Y el trabajo de las fuerzas interiores:

Z
1

virtual real

|{z} | {z }

w;x u;x E A dx

Igualando ambos trabajos para cualquier campo virtual w de desplazamientos se obtiene la forma dbil que se ha dado en la seccin anterior. e o

15

2.6

Soluciones dbiles aproximadas e

El MEF puede entenderse tambin desde el punto de vista matemtico como un procedie a miento de obtener soluciones aproximadas al problema planteado de forma dbil. e La idea es aproximar los espacios S y V, que son de dimensin innita, por espacios o S h y V h de dimensin nita: o Sh S y Vh V

Esta es la forma matemtica de interpretar el proceso de discretizacin. a o Una vez denidos estos dos espacios S h y V h mediante sus correspondientes bases de funciones, se resuelve la forma dbil del problema dentro ellos. Es decir, se busca una e h h funcin u 2 S que haga que se cumpla: o
Z
1 0 h w;x uh dx = ;x

wh f dx + wh (0) h

8wh 2 V h

(2.13)

Como en el caso del mtodo de Galerkin, la aportacin del MEF moderno dentro de e o esta interpretacin matemtica del mismo consiste en una forma sistemtica de construir o a a los espacios S h y V h . La relacin con el mtodo de Galerkin es clara, ya que 2.6 puede interpretarse como o e que las funciones Nj ; j = 1; : : : ; n forman la base del espacio V h y que el espacio S h se construye por traslacin de V h mediante la funcin g h . Las incgnitas son los coecientes o o o h en trminos de la base del espacio de funciones. cj que denen la funcin u o e Desde el punto de vista mecnico que proporciona el principio de los trabajos virtuales, a 2.13 se interpreta como que la solucin aproximada uh es la que satisface las condiciones o de equilibrio correspondientes a un conjunto restringido, V h , de desplazamientos virtuales.

2.7

Teorema de Lax-Milgram

El teorema de Lax-Milgram es un resultado bsico de la teor de problemas de contorno a a formulados variacionalmente. Desde el punto de vista que nos interesa, es un teorema de existencia y unicidad de la solucin de problemas planteados de forma dbil. o e Se menciona aqu el teorema simplemente como muestra de otros resultados similares. Estos resultados prueban que tiene sentido desde el punto de vista matemtico buscar la a solucin de problemas planteados de forma dbil, tal y como hace el MEF, ya que si se o e dan determinadas condiciones, esta solucin existe y es unica. o Con objeto de dar ms generalidad al enunciado del teorema, se cambiar la notacin a a o empleada en el enunciado del problema modelo 2.1 de la manera siguiente: a(w; u)
Z
1

w;x u;x dx

l(w)

w f dx + w(0) h 16

Con esta notacin, la forma dbil del problema es: encontrar una funcin u 2 S tal o e o que 8w 2 V se cumpla: a(w; u) = l(w) El teorema de Lax-Milgram aplica a problemas con condiciones de contorno esenciales homogneas, es decir, cuando S V. e El teorema establece que si la aplicacin a : S S ! <: o 1. Es bilineal: a(c1 u + c2 v; w) = c1 a(u; w) + c2 a(v; w) 2. Es simtrica: a(u; w) = a(w; u) e 3. Es continua: j a(u; w) j m k u k k w k 0<m<1 >0 c1 ; c2 2 <

4. Es denida positiva (el ptica): a(u; u) k u k2

y la aplicacin l : S ! < es lineal, entonces existe solucin para la forma dbil del o o e problema y esta solucin es unica. o p En el enunciado del teorema k k representa una norma, por ejemplo: k k= a(; ).

2.8

Propiedad de aproximacin optima o

La solucin uh del problema aproximado 2.13 y, por tanto, las soluciones que proporciona o el MEF, tienen una propiedad importante. Resulta que la solucin uh es optima en el o sentido de que es la mejor solucin posible dentro del nivel de aproximacin utilizado, es o o decir, dentro del espacio S h . Con la notacin introducida en la seccin anterior, el problema 2.13 puede redenirse o o h 2 S h tal que 8w h 2 V h se cumpla: como: encontrar una funcin u o a(wh ; uh ) = l(wh ) donde a y l cumplen con las hiptesis del teorema de Lax-Milgram. o Entonces, si se dene la funcin error en la solucin aproximada como: o o e = uh u se cumple que: 1. El error e es ortogonal al espacio V h V, es decir: a(wh ; e) = 0 8wh 2 V h

2. La aproximacin es optima en trminos de la norma denida por a (norma energtica): o e e a(e; e)

| {z }

a(~h u; uh u) u ~

8~h 2 S h u

(norma del error)2

17

3. Como corolario del punto 1 anterior, cuando S h V h (condiciones de contorno esenciales homogneas), resulta que la norma (energ del error e es igual al error e a) en la norma a(u; u) a(uh ; uh ):
| {z }

a(e; e)

= a(u; u) a(uh ; uh )
|
error en la (norma)2

(norma del error)2

{z

En efecto, por ser V h V se cumple: a(wh ; u) = l(wh ) y como tambin se cumple: e a(wh ; uh ) = l(wh ) 8wh 2 V h 8wh 2 V h

resulta que restando y aprovechando la bilinealidad de a: a(wh ; u uh ) = 0 8wh 2 V h

lo que constituye el primer resultado buscado. Adems, 8wh 2 V h , se tiene que (por la a elipticidad de a):
|

a(e + wh ; e + wh ) = a(e; e) +2 a(wh ; e) + a(wh ; wh )

0

{z

| {z }
0

| {z }
=0

{z

luego se cumple: a(e; e) = a(uh u; uh u) a(e + wh ; e + wh ) = a(~h u; uh u) ~ u 8~h 2 S h u

que es el segundo resultado buscado. El corolario se obtiene porque cuando S h V h , resulta que a(uh ; e) = 0, por ser uh 2 V h , y entonces: a(u; u) = a(uh e; uh e) = a(uh ; uh ) 2a(uh ; e) + a(e; e) = a(uh ; uh ) + a(e; e) Ntese que esto implica que, o a(uh ; uh ) a(u; u) esto es, que la solucin aproximada tiene una norma (energ inferior a la de la solucin o a) o exacta.

18

Cap tulo 3

Fundamentos matemticos (II) a

3.1 Introduccin o

En el cap tulo anterior se han introducido dos formas de entender el MEF desde el punto de vista matemtico. Por un lado, el MEF puede considerarse como un procedimiento de a residuos ponderados de tipo Galerkin y, por otro lado, puede ser entendido como un medio de obtener una solucin aproximada a un problema de campo formulado de manera dbil. o e En este cap tulo se proporciona una tercera v para entender el MEF, la v variaa a cional. El MEF puede tambin interpretarse como un caso de aplicacin del mtodo de e o e Ritz, es decir, como una forma aproximada de obtener la solucin de problemas de campo o que responden a un principio variacional. Este enfoque es muy importante dentro de las tcnicas actuales de desarrollo de elementos nitos. e No se debe perder de vista que la sistemtica del MEF, similar a la del clculo matricial a a de estructuras, fue previa a sus interpretaciones o justicaciones matemticas. El objetivo a de estos primeros cap tulos es proporcionar varias formas diferentes de interpretar una misma cosa, el formalismo del MEF, desde el punto de vista matemtico. a

3.2

Principios variacionales

Muchos problemas de la F sica con aplicacin prctica en Ingenier responden a un prino a a cipio variacional. Pueden encontrarse ejemplos en la Mecnica Clsica (principio de Haa a milton, principio de la m nima accin), en la Elasticidad (principio de la m o nima energ a potencial) e incluso en la Optica (principio de Fermat). Dentro de nuestro ambito de inters, la formulacin variacional, cuando existe, corres e o ponde a un nivel de abstraccin superior a la formulacin dbil introducida en el cap o o e tulo anterior. Es equivalente a la formulacin dbil, pero se presta de forma ms directa a la o e a obtencin de soluciones aproximadas. o Para el tipo de problemas en que estamos interesados, la formulacin variacional se o plantea normalmente del modo siguiente: 1. Dado un problema concreto se dene un espacio de posibles funciones solucin: S. o 19

2. El principio variacional establece que la solucin es la funcin que, dentro de ese o o espacio, hace m nimo, mximo o, simplemente, estacionario un funcional. a Un funcional es una \funcin de funciones", con valores reales: : S ! <. o Los funcionales que intervienen en los principios variacionales aplicables a problemas de Elasticidad tienen normalmente un signicado energtico. En este sentido, el principio e variacional introduce un propsito para la solucin del problema: hacer m o o nima, mxima a o estacionaria una magnitud energtica. e

3.3

Mtodo de Rayleigh-Ritz e

La utilidad de los principios variacionales desde el punto de vista del clculo numrico a e es que proporcionan una forma sencilla de obtener soluciones aproximadas. Una vez establecido el principio variacional aplicable, se procede de la forma siguiente: 1. Se aproxima el espacio S de posibles soluciones por uno de dimensin nita: S h S. o 2. Se escribe la solucin uh como combinacin lineal de las funciones base del espacio o o h. S 3. Se determinan los coecientes de esta combinacin lineal con la condicin de que o o h minimice, maximice o haga estacionario el funcional correspondiente, segn estau u blezca el principio variacional aplicable. Lo anterior constituye en esencia el mtodo de Rayleigh-Ritz. Este mtodo fue ine e troducido hacia 1909, con el objetivo de obtener soluciones numricas de problemas de e contorno utilizando los principios variacionales asociados a los mismos. La mayor dicultad del procedimiento est en la construccin del espacio S h de manera adecuada. El MEF a o no es ms que el mtodo de Ritz, pero con una forma sistemtica y general de construir a e a h. el espacio de funciones S Veamos un ejemplo de la forma clsica de utilizar el mtodo de Ritz. Sea una viga a e biapoyada de luz L, con seccin transversal de inercia I y mdulo de elasticidad E. La o o viga est sometida a una carga q por unidad de longitud perpendicular a su directriz. Se a busca la variacin de la echa u a lo largo de la directriz de la viga. o q
? u ? ? ? ? ? ? ? ? ? ? u

x L

El espacio S de soluciones posibles o admisibles est formado por todas las funciones a u : [0; L] ! < tales que cumplan las condiciones de contorno en los dos extremos de la ~ viga: u(0) = 0 y u(L) = 0 y que, adems, no quiebren su directriz (giro u;x continuo en ~ ~ a ~ ]0; L[). 20

Uno de los principios variacionales aplicables es el principio de la mnima energ a potencial, el cual establece que, dentro de S, la solucin del problema es la que hace o m nimo el funcional potencial total:
L 1 L E I (~;xx )2 dx u u q dx ~ (3.1) 2 0 0 Entonces, una solucin aproximada del problema puede obtenerse deniendo el espacio o S h como aquel que tiene por base las funciones:

(~) = u

sin(

ix ) L

i = 1; 2; : : : ; n

siendo n un nmero nito. u As las funciones uh 2 S h se pueden poner como: , ~ u = ~ y uh = ~;xx

n X i=1 h n X i=1

ci sin(

ix ) L

ci 2 <

ci (

i 2 ix ) sin( ) L L

ci 2 <

Sustituyendo estas expresiones en 3.1 se obtiene una expresin discreta del funcional o en la que ste slo depende de los coecientes ci . De este modo, el problema de la e o minimizacin del funcional se convierte en un problema algebraico. Para obtener los o coecientes ci que denen la solucin aproximada se plantea el sistema de ecuaciones: o @ =0 @ci i = 1; 2; : : : ; n

En este problema concreto, tomando n = 3, se obtiene: uh = x 1 4 q L4 3x )] [sin( ) + sin( 5 EI L L 243

Lo que da una echa mxima de: a uh = 0:0130172 max Frente al valor exacto de: umax = 0:0130208 q L4 EI q L4 EI

21

3.4

El problema elstico: notacin a o

El problema de la Elasticidad es uno de los primeros a los que se aplicaron las tcnicas e de elementos nitos. En lo que sigue, este problema ser utilizado con frecuencia como a veh culo concreto para explicar conceptos asociados al MEF. En esta seccin se plantea o el problema elstico, tanto en su forma fuerte como en una forma dbil, y se establece la a e notacin que ser utilizada. o a Se tiene un cuerpo elstico que ocupa un volumen 2 <3 y que est limitado por a a una supercie S = Sd [ St, de modo que Sd \ St = ;.1 La supercie S tiene una normal n ni .2 Sobre la supercie Sd se imponen desplazamientos d di ; mientras que sobre la t t a u supercie St se imponen tensiones i . Adems, sobre el volumen , acta un campo b bi de fuerzas por unidad de volumen. Las incgnitas del problema son los campos de desplazamientos, u ui ; de deformao ciones, e eij ; y de tensiones, s sij . Todos ellos denidos en el volumen . La formulacin fuerte del problema se hace en base al siguiente sistema de ecuaciones o diferenciales: 1. Ecuaciones de compatibilidad en . 1 e = (ru + rtu) 2 2. Ecuaciones de compatibilidad en Sd . u=d o ui = di (3.3) o 1 eij = (ui;j + uj;i ) 2 (3.2)

3. Ecuaciones constitutivas (o de respuesta del material) en . s = De 4. Ecuaciones de equilibrio en . div s + b = 0 5. Ecuaciones de equilibrio en St . t sn s n = o t sij nj = i (3.6) o sij;j + bi = 0 (3.5) o sij = Dijkl ekl (3.4)

Ntese que las ecuaciones diferenciales requieren que la funcin vectorial u sea derivable o o de orden 2 en , y que las funciones tensoriales s y e sean derivables de orden 1 en . t Entonces, puede ocurrir que la solucin fuerte no exista, en funcin de los datos b, d y . o o
1 Una barra sobre el s mbolo de un conjunto denota su clausura, esto es, la unin del conjunto y su o frontera: = [ S. 2 Salvo indicacin en contra, la notacin indicial supone que los o o ndices var de 1 a 3. Se utiliza an tambin el convenio de suma por e ndices repetidos.

22

Una de las posibles formulaciones dbiles del problema elstico corresponde a la aplicae a cin del principio de los trabajos virtuales para imponer las condiciones de equilibrio. De o esta manera las ecuaciones 3.5 y 3.6 son sustituidas por una ecuacin integral que impone o menos requisitos de derivabilidad. Para claricar el signicado de los distintos campos que intervienen en esta formulacin o dbil se introduce la notacin siguiente: e o 1. Los campos independientes o primarios, aquellos que pueden variarse libremente, se representarn con una tilde sobre su s a mbolo: ~, u, . . . s ~ 2. Los campos dependientes o derivados, aquellos que se obtienen a partir de campos primarios, se representarn utilizando como super a ndice el s mbolo del campo del 1 u u u t que derivan: s = D e , e = 2 (r + r )~ , . . . u 3. Los campos solucin del problema elstico se representarn sin tilde ni super o a a ndice: s, e, u. Adems, para simplicar la escritura de las integrales, se introduce tambin la siguiente a e notacin: o 1. Integrales de volumen: (f) (f ; g) 2. Integrales de supercie: [f ]St [f; g]St
Z Z

f d f : g d

Z Z

f d f : g d

St

St

Con esta notacin, una posible formulacin dbil del problema elstico queda como sio o e a gue. Dados los espacios S (desplazamientos admisibles) y V (variaciones o desplazamientos virtuales): ~ S = f~ 2 H 1 j u = d en Sd g u V = f~ 2 H 1 j ~ = 0 en Sd g u u encontrar u 2 S, de manera que: u t; u (su ; eu ) = (b; ~ ) + [ ~ ]St El conjunto H 1 es el espacio de las funciones que tienen derivadas primeras de cuadrado integrable en . 23 8~ 2 V u (3.7)

3.5

Principios variacionales en Elasticidad

A partir de la forma dbil 3.7 puede deducirse el principio de la m e nima energ potencial. a Se dene en el espacio S el funcional p : S ! <, o funcional \potencial total", como: 1 ~ t ~ p(~ ) = (su ; eu ) (b; u) [; u]St u 2 Y entonces, recordando la denicin de primera variacin de un funcional: o o p (~ ) u resulta que: u t; u u p (~ ) = (su ; eu ) (b; ~ ) [ ~ ]St Por lo tanto, la solucin u de la forma dbil 3.7 o principio de los trabajos virtuales o e hace que p (u) = 0. Es decir, hace estacionario el funcional p dentro del espacio S. As el principio de los trabajos virtuales, un planteamiento dbil del problema elstico, , e a es equivalente a buscar la estacionaridad del funcional p sobre los campos de desplaza miento u que cumplen la restriccin u = d en Sd . Este es precisamente el principio ~ o ~ variacional de la m nima energ potencial: De entre todos los campos de desplazamiento a que satisfacen la restriccin u = d en Sd (campos admisibles), la solucin u del problema o ~ o elstico es la que hace estacionario (m a nimo) el funcional potencial total. El principio de la m nima energ potencial proporciona una nueva manera de intera pretar la forma dbil del problema elstico a que da lugar el principio de los trabajos e a virtuales. A travs del mtodo de Ritz, el principio de la m e e nima energ potencial puede a utilizarse en la obtencin de soluciones aproximadas. Se trata de minimizar el funcional o potencial total dentro de un espacio S h S de dimensin nita. La solucin aproximada o o uh se dene a partir de la combinacin lineal de ciertas funciones base. Los coecientes o de esta combinacin se obtienen haciendo estacionario el funcional potencial total, lo que o da lugar a un sistema de ecuaciones idntico al que proporciona el mtodo de Galerkin. e e Como ya se ha apuntado en numerosas ocasiones, esta tcnica de aproximacin numrica e o e es muy anterior al desarrollo del MEF. La aportacin del MEF no es ms que una foro a ma sistemtica de construir los espacios de aproximacin S h y las funciones base para la a o denicin de uh . o El principio variacional de la m nima energ potencial, basado en un funcional con un a solo campo independiente u, es el que se utiliza para justicar la formulacin clsica del ~ o a MEF, tambin llamada formulacin \en desplazamientos". Sin embargo, existen muchos e o otros principios variacionales que pueden utilizarse en la resolucin del problema elstico. o a Empleando estos principios variacionales de la misma manera esbozada ms arriba se a obtienen formulaciones del MEF distintas de la clsica, algunas de las cuales tienen mua cha importancia hoy en d en el desarrollo de los llamados elementos nitos de \altas a prestaciones". Conviene mencionar los siguientes principios variacionales: 24 d (~ + ~ )j=0 u u d con ~ 2 V y 2 < u

1. Principio de la energ complementaria m a nima. Es el principio variacional conjugado del principio de la m nima energ potencial. El funcional asociado tiene tambin a e un solo campo independiente, un campo de tensiones ~. Este principio variacional s da lugar a la formulacin del MEF \en tensiones". Esta formulacin ha tenido un o o desarrollo muy peque~o porque se requiere que los campos de tensiones ~ cumplan n s las ecuaciones de equilibrio 3.5 y 3.6, lo cual plantea serias dicultades en la prctica. a 2. Principios variacionales mixtos o hbridos. Corresponden a puntos intermedios en tre el principio de la m nima energ potencial y el de la energ complementaria a a m nima. Se caracterizan porque los funcionales asociados tienen varios campos independientes. Dan lugar a formulaciones del MEF llamadas \mixtas" o \hbridas" por su carcter intermedio entre la formulacin en desplazamientos y la formulacin en a o o tensiones. Merecen citarse el principio de Hellinger-Reissner, con dos campos independientes u y ~; y el principio de Hu-Washizu, con cuatro campos independientes, ~ s ~ s ~ t, u, ~, e y ~ siendo este ultimo un campo de tensiones en el contorno. Los principios variacionales mixtos o multicampo juegan un papel muy importante en la moderna tecnolog de desarrollo de elementos nitos. Constituyen la manera de a justicar matemticamente el empleo de determinados articios numricos muy ecientes3 . a e

3.6

Ecuaciones de Euler y tipos de condiciones de contorno

En esta seccin se muestra cmo el principio de la m o o nima energ potencial conduce a a la forma \fuerte" de las ecuaciones de equilibrio 3.5 y 3.6 si se cumplen los requisitos de derivabilidad necesarios y que, por tanto, el planteamiento variacional y el planteamiento clsico en base a ecuaciones diferenciales son equivalentes. a El principio de la m nima energ potencial establece que la solucin u del problema a o elstico es la que hace estacionario el funcional potencial total: a 1 ~ t; ~ u p(~ ) = (su ; eu ) (b; u) [ u]St 2 ~ cuando ste se dene sobre los campos de desplazamientos u que cumplen con la condicin: e o ~ u = d en Sd (3.8)

Las restricciones 3.8 se conocen como condiciones esenciales de contorno asociadas al principio variacional. Forman parte de la denicin del conjunto de funciones admisibles. o As cuando se utilizan principios variacionales para la obtencin de soluciones aproxi, o madas, stas cumplen exactamente las condiciones esenciales de contorno del principio e variacional utilizado. Las ecuaciones diferenciales de Euler-Lagrange asociadas a un principio variacional se obtienen igualando a cero la primera variacin del funcional asociado. En el caso del o principio de la m nima energ potencial la primera variacin del funcional es: a o
3

Ver cap tulos 7 y 8.

25

p (~ ) = (su ; eu ) (b; ~ ) [ ~ ]St u t; u u con ~ = 0 en Sd . u Utilizando el teorema de la divergencia, la relacin anterior puede ponerse como: o u u t u p (~ ) = [su n; ~ ]St (div su ; ~ ) (b; ~ ) [; ~ ]St = [su n; ~ ]St (div su +b; ~ ) u u t u u Y como resulta que ~ es arbitrario, salvo que ~ = 0 en Sd , la condicin p = 0 u u o implica que: 1. Ecuacin de Euler (equilibrio en ): o div s + b = 0 en

2. Condiciones naturales de contorno (equilibrio en St ): sn = t en St

Las dos relaciones anteriores son precisamente las ecuaciones de equilibrio 3.5 y 3.6 del planteamiento fuerte del problema.

26

Cap tulo 4

Programacin del MEF (I) o

4.1 Introduccin o

En los cap tulos anteriores se ha visto que el MEF moderno no es un procedimiento novedoso desde el punto de vista matemtico. Puede asimilarse al procedimiento de residuos a ponderados de Galerkin o al mtodo de Ritz, mtodos numricos publicados medio siglo e e e antes de que aparecieran los elementos nitos. Lo innovador del MEF es la forma sistemtica y general que se da a estos procedimiena tos numricos, de manera que resultan fcilmente automatizables empleando un ordenador e a digital. El uso prctico del MEF es impensable sin la ayuda de los ordenadores. Todos los a desarrollos tericos y sus justicaciones matemticas acaban convertidos en l o a neas de un programa de ordenador. Una de las mejores formas de estudiar el MEF es analizar y tratar de entender la estructura interna de un programa de clculo basado en el mismo. a O, mejor todav programar de principio a n un programa que resuelva un problema a, sencillo. Esto proporciona una visin prctica del mtodo, de cmo se hacen en realidad o a e o las cosas, que es casi imposible conseguir de otro modo. El objetivo del presente cap tulo es ver cules son los componentes bsicos de un a a programa de clculo basado en el MEF y, a la vez, analizar una de las posibles arquitecturas a internas. Como primer paso, se proporciona la sistemtica bsica o \receta" del MEF, a a que fue previa a cualquiera de sus justicaciones matemticas. a

4.2

La \receta" del MEF

La sistemtica bsica o receta del MEF se explica aqu utilizando como veh a a culo el problema elstico denido en el cap a tulo anterior. El problema elstico se plantea de forma dbil como encontrar una funcin u 2 S de a e o manera que: (su ; eu ) = (b; ~ ) + [; ~ ]St u t u siendo: 27 8~ 2 V u (4.1)

~ S = f~ 2 H 1 j u = d en Sd g u y V = f~ 2 H 1 j ~ = 0 en Sd g u u La receta es la siguiente: 1. Dividir el dominio de clculo en subdominios o elementos: a =

[
i

\
i

i = ;

o 2. Dentro de cada elemento i , aproximar los campos incgnita que intervienen en la ecuacin integral 4.1: o ~ u = N ae ~ = N ae u eu = L N ae = B ae eu = L N ae = B ae su = D B ae donde N es una matriz de funciones conocidas (funciones de forma); ae y ae son vectores de coecientes incgnita, cuyo valor se hace coincidir normalmente con el o valor de los desplazamientos de ciertos puntos del elemento llamados nodos; L es una matriz de operadores (derivadas) y D es una matriz de coecientes elsticos. a Cada coeciente incgnita o movimiento nodal en ae es un grado de libertad para el o u campo aproximado u. El nmero de componentes de ae es el nmero de grados de ~ u libertad del elemento. 3. Las integrales que aparecen en 4.1 se hacen sumando las contribuciones de los distintos elementos, es decir, 8~ 2 V: u u t; u 0 = (su ; eu ) (b; ~ ) [ ~ ]St =
X
i

(4.2)

f(su ; eu )i (b; ~ )i [ ~ ]Sti g u t; u

4. La contribucin de cada elemento a la suma anterior es, sustituyendo 4.2: o aei t (Bt D B)i aei aei t (Nt b)i aei t [Nt]Sti = aei t (Kei aei fei ) t donde Kei (Bt D B)i es la matriz de rigidez del elemento i; y fei (Nt b)i + [Nt]Sti es el vector de cargas del elemento i. t

28

5. La suma de las contribuciones de los distintos elementos se hace a travs del llamado e proceso de ensamblaje. Ntese que algunos de los coecientes incgnita ae denidos a nivel de cada elemento o o sern normalmente compartidos por los elementos adyacentes, ya que existirn nodos a a comunes. Podr entonces establecerse una numeracin de carcter global para estos a o a coecientes. Llamaremos a al vector de los coecientes incgnita denidos para todos o los elementos. As un mismo coeciente tendr una unica numeracin global, como , a o componente del vector a, y una o varias numeraciones locales, como componente de uno o varios vectores locales ae . El proceso de ensamblaje consiste en ir acumulando las contribuciones de cada elemento en matrices y vectores de carcter global, esto es, ordenados segn la numea u racin global de grados de libertad. La forma discreta de 4.1 queda: o
X

elementos

aei t (Kei aei fei ) at (K a f ) = 0

8 a

(4.3)

donde K es la matriz de rigidez global y f es el vector global de cargas. 6. Como la relacin anterior ha de cumplirse 8 a, se llega nalmente al sistema de o ecuaciones lineales: Ka f = 0 el cual permite determinar los coecientes incgnita a. El nmero de componentes o u de a es el nmero total de grados de libertad del problema. u 7. De acuerdo con el paso 2, los coecientes a denen la solucin dentro de todos o los dominios elementales y, por agregacin, denen la solucin buscada en todo el o o dominio de clculo . a De todos los puntos anteriores, quiz el ms dif de ver para el novicio es el del a a cil ensamblaje (paso 5). El concepto se visualiza fcilmente con un ejemplo. Sea un problema a unidimensional en el que se utilizan los tres elementos nitos de dos nodos representados en la gura, y supongamos que slo se tiene un grado de libertad por nodo. o

1
x

2
x

3
x

4
x

La numeracin global de nodos y elementos es la que se muestra en la gura. En el o elemento 1 la matriz de rigidez es: 29

K1 = el vector de cargas:

"

K111 K112 K121 K122

( )

f1 =

f11 f12

y el vector de grados de libertad o incgnitas nodales: o a1 =

(

a11 a12

Y anlogamente para los elementos 2 y 3. La relacin entre grados de libertad locales a o y globales es: a11 a1 a21 a2 a31 a3 a12 a2 a22 a3 a32 a4 La suma de contribuciones de los elementos que se escribe en 4.3 es: (a11 a12 )
(" #( ) ( ))

K111 K112 K121 K122

(a21

a11 f11 + a12 f (" # ( 12 ) ( )) K211 K212 a21 f21 a22 ) + K221 K222 a22 f (" # ( 22 ) ( )) K311 K312 a31 f31 =0 (a31 a32 ) K321 K322 a32 f32
9 8 38 > a1 > > f11 > > > < = < 7> a > > f +f 7 2 12 21 7 5 > a3 > > f22 + f31 > > > > > > : ; : 99 >> >> >> == >> >> >> ;;

la cual puede ponerse como:

82 > K111 > >6 < K 6 (a1 a2 a3 a4 ) 6 121 >4 0 > > :

K112 0 0 K122 + K211 K212 0 K222 + K311 K312 K221 K322 0 K321

=0

a4

f32

El proceso mediante el cual se construyen directamente la matriz de rigidez global K y el vector global de cargas f es el proceso de ensamblaje. La numeracin global de grados o de libertad sirve para saber en qu posicin de la matriz K y el vector f hay que acumular e o los elementos de las matrices y vectores correspondientes a los elementos.

4.3

Clculos por el MEF: datos y resultados a

Los datos bsicos que se requieren para llevar a cabo un clculo por el MEF son los a a siguientes: 30

1. Denicin de la geometr del dominio de clculo y discretizacin del mismo. o a a o Esto se hace dando una lista de nodos y de elementos. Cada nodo es un punto dentro del dominio de clculo y se dene mediante un nmero de orden o etiqueta a u identicativa (nmero de nodo) y sus coordenadas en el sistema de referencia elegido. u Cada elemento corresponde a uno de los subdominios en que se divide el dominio de clculo. Se dene mediante un nmero de orden o etiqueta identicativa (nmero de a u u elemento) y una lista de nmeros de nodo, la cual se conoce tambin con el nombre u e de conectividad del elemento. Como se ver ms adelante, la geometr del elemento a a a queda completamente denida a partir de la formulacin interna del elemento y de o las coordenadas de sus nodos. El conjunto de nodos y elementos constituye lo que se conoce como mallado o malla de elementos nitos. 2. Atributos o propiedades de los elementos. Estas propiedades dependen de la clase de elemento nito que se est utilizando. As e como en todos los casos hay que identicar el material que constituye el subdominio o elemento, existen otros atributos que necesitan o no ser especicados en funcin o del problema y el tipo de elemento. Ejemplos t picos son: el espesor, la seccin o transversal y su orientacin en el espacio, las direcciones de anisotrop . . . o a, 3. Propiedades de los materiales. Cada material se identica mediante un nmero o etiqueta. A cada material se u le asocia un modelo matemtico para representar su comportamiento (elasticidad, a plasticidad, hiperelasticidad, . . . ) y se denen en cada caso los parmetros numricos a e del modelo matemtico elegido (p.ej. mdulo de elasticidad, coeciente de Poisson, a o tensin de uencia, . . . ). El modelo matemtico de comportamiento del material se o a conoce muchas veces con el nombre de modelo constitutivo o ley de comportamiento. 4. Condiciones de contorno. En problemas mecnicos se distinguen dos clases de condiciones de contorno: las a condiciones de contorno en desplazamientos y las condiciones de contorno en fuerzas. Las primeras son restricciones de tipo cinemtico y corresponden normalmente a a las condiciones de contorno que hemos llamado esenciales en el planteamiento variacional. En los manuales de usuario de los programas de elementos nitos estas restricciones se conocen como condiciones de contorno propiamente dichas. Se caracterizan porque afectan directamente a la variable de campo del problema. Ejemplos t picos son los desplazamientos restringidos o los desplazamientos impuestos. Estas condiciones se aplican directamente sobre los nodos, limitando o anulando sus desplazamientos, y se dan mediante una lista de nodos a los que se asocia un cdigo o que corresponde a la condicin de contorno que se desea aplicar. o Las condiciones de contorno en fuerzas son las que normalmente se conoce como acciones en los manuales de usuario de los programas de elementos nitos. Se trata

31

de fuerzas aplicadas sobre nodos (cargas puntuales), presiones sobre la supercie de los elementos o fuerzas distribuidas por unidad de volumen en los elementos. En otro tipo de problemas, no mecnicos, tambin puede hacerse la distincin entre a e o estas dos clases de condiciones de contorno. Por ejemplo, en problemas de transferencia de calor, las condiciones de contorno en desplazamientos corresponden a temperaturas impuestas; y las condiciones de contorno en fuerzas o acciones corresponden a ujos o fuentes de calor. 5. Otros datos. Los datos que se mencionan en los puntos anteriores son imprescindibles en cualquier clculo por elementos nitos. Existen otras clases de datos que pueden no ser a necesarios en funcin del tipo de problema que se trate de resolver. Se puede peno sar, por ejemplo, en datos de condiciones iniciales del dominio (tensin, velocidad, o temperatura, . . . ), en relaciones impuestas entre el movimiento de diferentes nodos (v nculos cinemticos) . . . a A partir de los datos anteriores, los resultados bsicos que proporciona un clculo por a a elementos nitos corresponden a dos grandes categor as: 1. Variables nodales. Son los resultados que denen la variable de campo incgnita bsica o sus derivadas o a con respecto al tiempo. Dependiendo del problema, son los desplazamientos, giros, velocidades, temperaturas, etc. . . . Se trata de valores que se obtienen directamente en los nodos del mallado y su orden de aproximacin suele ser bastante bueno, aun o con discretizaciones muy gruesas. En esta categora se incluyen tambin las reaccio e nes en los puntos a los que se aplican condiciones de contorno en desplazamientos. 2. Variables elementales. Son los resultados que corresponden a campos derivados del campo incgnita bsica o a a travs de derivadas espaciales. Son, por ejemplo, los campos de deformaciones, e tensiones, ujo de calor, etc. . . . Se trata de resultados calculados en puntos internos de los elementos, aunque a veces se extrapolen luego a los nodos. Su orden de aproximacin es peor que el de las variables nodales: si se utilizan discretizaciones o gruesas pueden cometerse errores importantes.

4.4

Flujo general en un programa de EF para clculo lineal a

Para llegar desde los datos bsicos a los resultados bsicos descritos en la seccin anterior, a a o se pasa por el proceso principal del programa de clculo. En un clculo de naturaleza a a lineal, como el utilizado para explicar la \receta" del MEF que se da en la seccin 4.2, o el problema de obtener la solucin se reduce a resolver un sistema de ecuaciones lineales. o Entonces, lo que debe hacer el programa de clculo es, a partir de los datos, montar a la matriz de coecientes del sistema y su vector de trminos independientes (esto es, la e 32

matriz de rigidez global y el vector global de cargas), resolver el sistema de ecuaciones y, nalmente, obtener los resultados necesarios a partir de la solucin del sistema de o ecuaciones. El ujo general de un programa de elementos nitos para clculo lineal podr ser el a a siguiente: 1. Entrada de datos En esta etapa se realiza la lectura y/o generacin de los datos y su acoplamiento deno tro de la estructura de datos del programa. Se realizan tambin las comprobaciones e o chequeos bsicos de consistencia de los datos introducidos por el usuario. a 2. Tareas preliminares En esta fase realizan las labores previas al ensamblaje del sistema global de ecuaciones. Por ejemplo, la numeracin global de grados de libertad, el clculo del semiano a cho de banda, las comprobaciones de que se dispone de espacio suciente (memoria, disco), etc. . . . 3. Construccin de la matriz de rigidez global K y del vector global de cargas f o El proceso puede estructurarse en dos etapas: un bucle que recorre los elementos y, luego, la imposicin de las acciones nodales directas: o Para cada elemento:

(a) Calcular matriz de rigidez elemental y vector elemental de cargas. (b) Ensamblar los elementos de la matriz y del vector en los lugares correspondientes de la matriz de rigidez global y del vector global de cargas.

Inclusin en el vector global de cargas de las acciones (cargas) nodales introduo cidas directamente. 4. Resolver el sistema K a = f 5. Salida de resultados nodales: a 6. Elaboracin y salida de resultados elementales o Esta parte se organiza tambin mediante un bucle que recorre los elementos: e Para cada elemento:

(a) Recuperar resultados nodales correspondientes a sus nodos ae . (b) Calcular resultados elementales a partir de los resultados nodales: su = D B ae . (c) Salida de resultados elementales.

Ntese que el ncleo del programa es el clculo de la matriz de rigidez elemental y del o u a vector elemental de cargas. En ese punto intervienen la llamada tecnolog de elementos, o a

33

tcnicas de construccin del espacio de aproximacin, y los modelos matemticos de leyes e o o a de comportamiento de los materiales. El otro punto importante es la forma de resolver el sistema de ecuaciones, que es la etapa donde ms tiempo de ordenador se gasta en los problemas lineales. La tecnolog a a de elementos, los modelos matemticos de leyes de comportamiento y los procedimientos a de solucin forman los tres pilares de las tcnicas de elementos nitos. o e El resto de las etapas del programa, las de organizacin y transferencia de informacin, o o constituyen lo que se conoce como interfase con el usuario. Son etapas muy importantes, sobre todo desde el punto de vista de la facilidad de uso del programa, y se abordan hoy en d con apoyo grco. Sin embargo se han dejado fuera del alcance de estas notas, por a a no corresponder estrictamente al clculo con elementos nitos. a

34

Cap tulo 5

Tecnolog de elementos (I) a

5.1 Introduccin o

En los cap tulos anteriores hemos hablado de los fundamentos matemticos del MEF. a Hemos visto cmo, dado un problema de campo denido en un determinado dominio, se o puede encontrar una solucin aproximada buscndola dentro de un grupo o familia de o a funciones denidas en el dominio y que cumplan una serie de condiciones. Hemos visto tambin que estos procedimientos de aproximacin son anteriores al dee o sarrollo del MEF y que lo realmente innovador del MEF es su \receta" para construir de forma sucientemente general y sistemtica esas familias de funciones dentro de las cuales a se busca luego la solucin aproximada. o De acuerdo con la \receta" del MEF, una funcin dentro de una de esas familias, funo cin denida en el dominio completo, se obtiene agregando funciones denidas en dominios o locales, elementos nitos o subdominios del dominio completo. La denicin de dichas funciones locales constituye el ingrediente fundamental del o MEF, la calidad de la aproximacin depende en gran medida de ellas. Las funciones o locales se conocen por distintos nombres, el ms extendido es quiz el de funciones de a a forma, porque proporcionan las formas que puede adoptar localmente el campo incgnita, o o porque esas mismas funciones se utilizan muchas veces tambin para denir la geometr e a del dominio local. Un nombre ms adecuado desde el punto de vista matemtico ser el a a a de funciones de aproximacin local. o Dentro de la teor general del MEF, la parte que trata sobre cmo denir estas a o funciones de forma, sobre las condiciones que deben cumplir y sobre cmo obtener a partir o de ellas la aproximacin extendida al dominio completo, se conoce globalmente como o tecnolog de elementos. Este ser el objeto de los prximos cap a a o tulos. Debe decirse que la eleccin las funciones de forma est condicionada no slo por la o a o geometr de los dominios locales o elementos nitos, sino tambin por el tipo de problema a e de campo que se intenta resolver (p.ej. elasticidad, plasticidad, transferencia de calor, electromagnetismo) y por la manera de abordarlo (nmero de campos independientes). u Como veremos, el papel fundamental que juegan dentro del esquema general del MEF

35

estas funciones locales o, abusando del lenguaje, los elementos1 , es un papel de interpolacin desde los nodos de la geometr de las acciones y, sobre todo, de las variables de o a, campo. Por eso a veces estas funciones se conocen tambin como funciones de interpolae cin. o

5.2

Formulacin convencional en desplazamientos o

La formulacin ms clsica del MEF recibe el nombre de formulacin en desplazamientos. o a a o La formulacin toma el nombre de su aplicacin a problemas de elasticidad, que fue una o o de las primeras areas en las que se utiliz el MEF y de donde se ha extra mucho de su o do lenguaje. En otros campos de aplicacin, habr que llamarla formulacin en temperaturas (proo a o blemas de transferencia de calor), en velocidades (mecnica de uidos), etc. a La formulacin se basa en la aproximacin de un solo campo independiente: los deso o plazamientos en sentido generalizado. Dentro de cada dominio elemental o elemento nito, el campo independiente se aproxima o interpola mediante funciones de forma que cumplen las condiciones siguientes: 1. Cada funcin de forma est asociada a un punto dentro del elemento, de manera o a que vale 1 en ese punto y se anula en los puntos asociados a las otras funciones de forma. Estos puntos especiales dentro de cada elemento reciben el nombre de nodos. 2. Las funciones de forma correspondientes a cada elemento estn denidas de manera a que hay continuidad del campo independiente al pasar de un elemento a otro, es decir, al agregar las aproximaciones locales para construir la aproximacin extendida o al dominio completo se obtiene una funcin continua. o 3. Las funciones de forma son capaces de representar de manera exacta campos con un m nimo orden de variacin, el cual depende del tipo de problema. Por ejemplo, en o elasticidad las funciones de forma son capaces de representar exactamente campos de desplazamientos con variacin lineal en el espacio; en problemas de exin de o o placas, campos de echas con variacin cuadrtica, etc. o a 4. Las funciones de forma dentro de cada elemento son sucientemente lisas o \suaves" como para permitir el clculo de las derivadas que sean necesarias para el planteaa miento de la forma dbil del problema. e La primera condicin no es estrictamente necesaria, pero facilita la interpretacin de o o los resultados y se ha convertido en algo que muchas veces se sobrentiende. Las otras tres condiciones derivan del fundamento matemtico del mtodo y, en conjuna e to, son una condicin suciente para la convergencia hacia la solucin exacta del problema o o al aumentar el nmero de elementos o grado de discretizacin. u o
1 Dentro del MEF las funciones locales son parte tan intr nseca de los elementos que ambos conceptos resultan dif ciles de disociar en el lenguaje habitual. Por ejemplo, cuando se habla de un elemento cuadriltero de ocho nodos, se entiende que se est utilizando un conjunto muy concreto de ocho funciones de a a forma.

36

5.3

Ejemplo: elemento cuadriltero a

Cumpliendo las cuatro condiciones que se dan en la seccin anterior, hay muchas posibio lidades para las funciones de forma. Veamos como ejemplo una que se emplea mucho en la prctica para elementos rectangulares. a Sea el elemento rectangular e <2 representado en la gura:

4
z

3
z

Y
6 z z

1
-

2 X

Parece razonable situar un nodo en cada esquina del rectngulo. Esto dar lugar a a a cuatro funciones de forma. Sea (x; y) el campo independiente que se quiere aproximar dentro del elemento, por ejemplo, el campo de temperaturas. Parece tambin razonable utilizar una aproximacin e o para basada en polinomios, que son las funciones ms sencillas que pueden encontrarse: a ' 1 + 2 x + 3 y + 4 xy + 5 x2 : : : i 2 <

Si unicamente se dispone de cuatro funciones de forma, slo se tienen cuatro parmetros o a independientes para aproximar dentro del elemento: los coecientes que multipliquen a cada una de las funciones. De este modo, nuestra aproximacin a ser, como mucho: o a ' 1 + 2 x + 3 y +
| {z }

4 xy

necesarios arbitrario Los tres primeros trminos son necesarios para tener una aproximacin lineal completa e o del campo dentro del elemento, pero el ultimo trmino puede elegirse arbitrariamente. e Se ha seleccionado un trmino simtrico en x e y, el ms simple posible, para no dar ms e e a a peso en la aproximacin a ninguna direccin del plano <2 . o o De este modo, se puede evaluar la aproximacin 5.1 a en los cuatro nodos que hemos o elegido para el elemento, ya que se conocen sus coordenadas (xi ; yi ) para i = 1; : : : ; 4. Resulta entonces:

| {z }

(5.1)

37

8 9 > 1 > > > > < > =

2

> 3 > > > > > : ;

=6

6 6 4

1 1 1 1

x1 x2 x3 x4

y1 y2 y3 y4

x1 y1 x2 y2 x3 y3 x4 y4

9 38 > 1 > > > = 7> > < 7 2 7 5 > 3 > > > > > ; :

donde i = (xi ; yi ), con i = 1; : : : ; 4. La relacin anterior puede escribirse en notacin o o matricial: ae = C q Salvo conguraciones patolgicas del elemento, la matriz C ser regular, con lo que se o a pueden determinar los coecientes q de la interpolacin a partir del valor ae en los nodos o del elemento del campo que se quiere aproximar: q = C1 ae Entonces, resulta nalmente: (x; y) '
h | i

1 x y xy
{z

modos de la aprox. Donde N N(x; y) es la matriz de funciones de forma. Cada columna de esta matriz es una funcin de forma. o Obsrvese que, de este modo, se cumplen las condiciones que se hab exigido en la e an seccin anterior: o 1. Cada funcin de forma vale 1 en su nodo asociado y se anula en los tres restantes. o 2. La aproximacin global a (x; y) que se construya por agregacin de aproximaciones o o locales de esta clase ser continua a travs de las fronteras entre los elementos, ya a e que la aproximacin local a lo largo de uno de los lados del elemento slo depende o o de los valores i en los nodos extremos de ese lado. 3. La aproximacin local es \suave" dentro del elemento, en el sentido de que estn o a denidas sus derivadas de cualquier orden. 4. La aproximacin es capaz de interpolar exactamente campos (x; y) con variacin o o lineal en el dominio de clculo. a En este ejemplo, la matriz N de funciones de forma es:
6 6 6 t N =6 6 6 4 2
1 4 (1 1 4 (1 + 1 4 (1 + 1 4 (1 xxc a )(1 xxc a )(1 xxc a )(1 xxc a )(1

q = P q = PC1 ae = N ae

+ +

yyc b ) yyc b ) yyc b ) yyc b )

3 7 7 7 7 7 7 5

38

donde: a = 1 (x2 x1 ) = 1 (x3 x4 ) b = 1 (y4 y1 ) = 1 (y3 y2 ) 2 2 2 2 xc = 1 (x1 + x2 + x3 + x4 ) 4 yc = 1 (y1 + y2 + y3 + y4 ) 4

Ntese que se podr haber determinado directamente, ms o menos a estima, las o an a funciones de forma que hemos obtenido. De hecho, sta ha sido la forma ms habitual e a de hacerlo: escribir directamente funciones que cumplan las condiciones que se dan en la seccin anterior. o

5.4

La transformacin isoparamtrica o e

Como se puede ver en el ejemplo, las funciones de forma dependen de la geometra del elemento. Recordemos adems que, dentro de la sistemtica del MEF, las funciones de a a forma y sus derivadas aparecen dentro de integrales extendidas a todo el elemento o a su contorno. Desde el punto de vista prctico, lo anterior presenta dicultades de cara a la automaa tizacin del clculo ya que, en principio, la geometr de un elemento puede ser diferente o a a de la de los dems, lo que da lugar a funciones de forma distintas de un elemento a otro a y a un tratamiento diferente las integrales que hay que efectuar. Ntese, por ejemplo, o lo engorroso de invertir la matriz C en un caso general para determinar las funciones de forma de cada elemento de la malla; o de las dicultades algor tmicas derivadas de que el dominio de integracin cambie de unos elementos a otros. o Desde el punto de vista de la automatizacin ser muy ventajoso trabajar siempre, o a independientemente del elemento concreto, con el mismo dominio de interpolacin y con o las mismas funciones de forma. De esta manera se simplicar la evaluacin de estas a o funciones y de sus derivadas y se conseguir realizar las integrales siempre sobre el mismo a dominio de integracin, operacin programable mucho ms fcilmente. o o a a >Cmo se compagina entonces la necesaria heterogeneidad geomtrica de los elementos o e que forman el mallado con la automatizacin sencilla del clculo? La respuesta est en la o a a llamada transformacin isoparamtrica de la geometr de los elementos. Veamos la forma o e a en que se trabaja en la prctica con elementos bidimensionales (<2 ) de cuatro lados. El a procedimiento es totalmente general, aplicable igualmente a elementos tridimensionales (<3 ) o unidimensionales (<). La secuencia de operaciones es la siguiente: 1. Se proyecta el dominio elemental o elemento nito en un dominio estndar, siempre a el mismo:

39

 PROYECCION x = x(; ) , y = y(; ) Y

6 x B B x B B B B ( ( Bx ( ( ( ( x( ( ( @

(-1,1)
x 6

(1,1)
x -

@ x x

(-1,-1)

(1,-1)

Coordenadas locales (; )

Coordenadas reales (x; y)

2. La proyeccin puede interpretarse tambin como un cambio de variables dentro de o e cada dominio local: Cambio de variables
(

x = x(; ) y = y(; )

(5.2)

3. Las funciones de interpolacin o funciones de forma se denen en coordenadas locales, o por ejemplo: 2 1 3 2 3 (1 )(1 ) 6 4 7 N 6 1 (1 + )(1 ) 7 6 1 7 6 4 7 6 N2 7 7 Nt = 6 1 (5.3) 6 (1 + )(1 + ) 7 = 6 N 7 6 4 7 4 3 5 4 1 5 N4 (1 )(1 + )
4

4. Las derivadas con respecto a las coordenadas reales (x; y) se calculan a partir de la matriz jacobiana del mapeo 5.2. Si es el campo o funcin incgnita, se tiene: o o (x; y) = (x(; ); y(; )) ' con lo que: @ @ @ @ @ X @Ni @ X @Ni @ = + ' + i i @x @ @x @ @x @x i @ @x i @ @ @ @ @ @ @ X @Ni @ X @Ni = + ' + i i @y @ @y @ @y @y i @ @y i @ o, en forma matricial:
X
i

i Ni = N ae

40

8 < :

@ @x @ @y

9 = ;

=4

@ @x @ @y

@ @x @ @y

8 P 9 2 @Ni < i i @ = = J1 4 ' J1 P : i @Ni ; i @ 2 3 5

38 < 5 :

@ @ @ @

9 = ;

= J1
@N @ @N @

8 < : 3

@ @ @ @

9 = ;

'

5 ae

donde J es la matriz jacobiana del mapeo 5.2, es decir:

@x @ @x @ @y @ @y @

J=4

5. Las integrales necesarias dentro de la sistemtica o \receta" del MEF se realizan a en el dominio o cuadrado estndar utilizando el cambio de variable denido por el a mapeo 5.2, es decir, si f(x; y) es la funcin que debe integrarse en el elemento e : o
Z

f(x; y) dx dy =

f(; ) det J d d

La proyeccin 5.2 se conoce como transformacin isoparamtrica cuando se utilizan las o o e propias funciones de forma para interpolar la geometr es decir, cuando las coordenadas a, reales (x; y) se obtienen a partir de las coordenadas locales (; ) como: x = x1 N1 (; ) + x2 N2 (; ) + x3 N3 (; ) + x4 N4 (; ) y y = y1 N1 (; ) + y2 N2 (; ) + y3 N3 (; ) + y4 N4 (; ) Esta es prcticamente la unica proyeccin que se utiliza en las aplicaciones. Desde a o el punto de vista terico podr emplearse tambin transformaciones subparamtricas, o an e e cuando el orden de la proyeccin es inferior al de las funciones de forma, o transformaciones o superparamtricas, cuando el orden es superior. e Es importante darse cuenta de que, con la forma de operar que se describe ms arriba, a la geometr real del elemento slo aparece en la formulacin a travs de la matriz jacobiana a o o e J. Un punto asimismo importante es que, al estar las funciones de forma denidas en el cuadrado estndar, si el elemento est muy distorsionado (tiene angulos muy distintos a a del recto), la interpolacin de la variable de campo en el espacio real (x; y), puede ser o muy distinta de la que tiene lugar en coordenadas locales (; ). Esto es, las funciones de forma reales pueden llegar a ser incluso cualitativamente muy diferentes de los polinomios denidos en coordenadas locales. Una interpolacin polinmica lineal en el dominio o o 41

estndar no tiene por qu ser tambin lineal en el espacio real. Es posible que ni siquiera a e e sea completamente polinmica. o La distorsin de los elementos da lugar as en la prctica a una degradacin de la o a o calidad de la solucin aproximada. En general, esto se debe a que las funciones de forma o en el dominio f sico dejan de ser exclusivamente polinomios y aparecen funciones racionales. Dichas funciones racionales dan lugar a un espacio de aproximacin en el que algunas de las o derivadas espaciales, necesarias por ejemplo para interpolar las tensiones en un problema de elasticidad, ya no pertenecen al propio espacio de aproximacin, lo que da lugar a una o prdida en el orden de convergencia. e

5.5

Integracin numrica o e

Las integrales que aparecen en la sistemtica o \receta" del MEF, que son integrales a extendidas a los dominios elementales o elementos nitos, se hacen numricamente. Lo e ms corriente es emplear la cuadratura de Gauss, que es muy eciente para la integracin a o de funciones polinmicas. o Tras el cambio de variables que supone la transformacin isoparamtrica, las integrales o e se extienden a un dominio estndar: una recta, un cuadrado o un hexaedro, segn la a u dimensionalidad del problema.
u u u(1,1,1)

(-1,1)
u u u

(1,1)
u

-1

1
u u u(1,-1,1)

(-1,-1)

(1,-1)

En el dominio estndar unidimensional, la cuadratura de Gauss se escribe como: a

Z
1

f() d =

n X 1

wi f(i )

donde: n = nmero de puntos de integracin u o wi = peso asociado al punto de integracin i o o i = coordenada del punto de integracin i

42

Con n puntos se integra exactamente un polinomio de grado 2n 1. Las coordenadas de los puntos de integracin y los pesos asociados se tabulan en los textos de clculo o a numrico2 . e Anlogamente, en el dominio estndar bidimensional: a a
Z
1

f (; ) dd =

j=1 i=1

n n XX

wj wi f (i ; j )

y en el dominio tridimensional:
Z
1

1 1 1

f(; ; ) ddd =

m=1 j=1 i=1

n n n X XX

wm wj wi f (i ; j ; m )

Ntese que en esta clase de integracin numrica las integrales se calculan evaluando o o e la funcin subintegral en una serie de puntos y sumando estos valores tras multiplicarlos o por unos coecientes de peso. En principio, se debe escoger un nmero n de puntos en cada direccin estndar suu o a ciente para obtener un valor exacto de la integral. En la prctica, en la mayor parte a de los casos, se subintegra, esto es, se utilizan menos puntos de los necesarios para el clculo exacto de la integral. Las razones para proceder as se mostrarn en los cap a a tulos siguientes.

5.6

Algunas familias corrientes de funciones de forma

Se presentan a continuacin las dos familias de funciones de forma ms corrientes en la o a prctica. Para su denicin se utiliza el dominio bidimensional estndar, en el plano a o a (; ). Su generalizacin al dominio estndar tridimensional, o su particularizacin al o a o dominio unidimensional, son relativamente inmediatas. La familia que ms se usa en las aplicaciones es la familia serend a pita, deducida originalmente por mera observacin y llamada as por referencia al Pr o ncipe de Serendip3 . Los elementos de orden primero, segundo y tercero dentro de esta familia se representan en la gura. Su principal ventaja es que, salvo excepciones, slo tienen nodos en el contorno o del elemento, que es donde hacen realmente falta para la conexin con otros elementos. o
Ver por ejemplo la referencia 3 de la bibliograf recomendada. a Personaje de novela clebre por sus descubrimientos fortuitos, creado por Horacio Walpole, novelista e ingls del siglo XVIII. e
3 2

43

Familia Serend pita 4

u

3
u

4
u

7
u

3
u

4
u

10
u

9
u

3
u u 8 u 7

8 u
u u u u

6 u
u

11 u 12 u
u u u

1 Orden 1

5 Orden 2

5 6 Orden 3

Las funciones de forma serend pitas son las siguientes: 1. Primer orden. Funciones bilineales. 1 Ni = (1 + i )(1 + i ) 4 donde (i ; i ) son las coordenadas nodales. 2. Segundo orden. Funciones bicuadrticas. a 1 Ni = (1 + i )(1 + i )(i + i 1) 4 1 Ni = (1 2 )(1 + i ) 2 1 Ni = (1 2 )(1 + i ) 2 3. Tercer orden. Funciones bicbicas. u Ni = 1 (1 + i )(1 + i )[9( 2 + 2 ) 10)] 32 9 (1 2 )(1 + i )(1 + 9i ) 32 9 Ni = (1 2 )(1 + i )(1 + 9i ) 32 Ni = i = 1; : : : ; 4 (nodos de vrtice) e i = 5; 6; 9; 10 i = 7; 8; 11; 12 i = 1; : : : ; 4 (nodos de vrtice) e i=5y7 i=6y8 i = 1; : : : ; 4

La otra familia que se utiliza en las aplicaciones es la llamada familia lagrangiana. En esta familia, las funciones de forma se obtienen por producto de polinomios de Lagrange en una variable. Esta tcnica de produccin de las funciones de forma conduce, cuando se e o emplea en elementos de alto orden, a la introduccin de muchos trminos polinmicos de o e o orden superior. Estos trminos no mejoran la calidad de la aproximacin, porque no dan e o lugar a polinomios completos. 44

Es raro en la prctica encontrar elementos lagrangianos de orden superior a 2. Se a utilizan cuando los elementos serend pitos presentan algn tipo de dicultad de uso, por u ejemplo, cuando se requiere una representacin precisa de la geometr en el interior del o a elemento (problemas de inestabilidad de estructuras laminares) o cuando se trabaja en problemas de contacto. Los elementos de orden primero y segundo dentro de esta familia se representan en la gura. El elemento de primer orden es idntico al correspondiente a la familia serend e pita. Familia Lagrangiana 4
u

3
u

4
u

7
u u9

3
u u 6

8 u

1 Orden 1

5 Orden 2

Las funciones de forma lagrangianas son las siguientes: 1. Primer orden. Funciones bilineales. 1 Ni = (1 + i )(1 + i ) 4 donde (i ; i ) son las coordenadas nodales. 2. Segundo orden. Funciones bicuadrticas. a 1 Ni = (i + )(i + ) 4 i = 1; : : : ; 4 (nodos de vrtice) e i=5y7 i=6y8 i=9 i = 1; : : : ; 4

1 Ni = (1 2 )(i + ) 2 1 Ni = (1 2 )(i + ) 2 Ni = (1 2 )(1 2 )

45

Cap tulo 6

Tecnolog de elementos (II) a

6.1 Introduccin o

Dentro de las aplicaciones prcticas del MEF en ingenier mecnica y en ingenier de a a a a estructuras, tienen mucha importancia los problemas de exin de placas y los relacionados o con estructuras laminares. Se pueden encontrar multitud de ejemplos en el proyecto de vasijas de presin, buques, aeronaves, puentes, cubiertas y otras grandes estructuras. o Los problemas de clculo de tensiones en placas y lminas son problemas complejos, a a donde slo existe solucin anal o o tica para geometr y distribuciones de carga sencillas. as Antes de que el MEF estuviera disponible, los proyectistas trabajaban con soluciones aproximadas o recurr a los ensayos a escala. El MEF, con su potencia para representar an geometr y distribuciones complejas de cargas, es hoy en da la herramienta de clculo as a ms avanzada con que se cuenta para abordar clculos de este tipo. a a La importancia tecnolgica que tiene el clculo de placas y lminas se pone de manio a a esto por el gran esfuerzo de desarrollo de elementos nitos para esta clase de problemas llevado a cabo por los investigadores del MEF durante ms de dos dcadas. Hoy en d a e a siguen apareciendo peridicamente nuevas formulaciones de elementos lmina, sobre todo o a para aplicaciones no lineales, lo que es s ntoma de la dicultad del problema y de que queda todav camino por recorrer. a En este cap tulo se pretende introducir los conceptos bsicos relacionados con una a clase particular de elementos nitos: los elementos placa y los elementos lmina. Se haa blar fundamentalmente de exin de placas, introduciendo el tema mediante un problema a o anlogo, pero ms sencillo, como es la exin de vigas. a a o No se entrar en problemas de clculo de lminas, debido a la dicultad de las teoras de a a a deformacin de estos elementos estructurales. Sin embargo, el tratamiento del problema o con el MEF sigue la pauta general. Baste con decir que en el clculo de lminas se a a resuelve simultneamente un problema de exin y otro de elasticidad bidimensional en a o las coordenadas locales de la supercie media de la lmina. En este sentido, muchos de a los conceptos que se van a introducir en relacin con la exin de placas son directamente o o transplantables a las lminas, en particular, la diferencia entre las formulaciones C 0 y C 1 . a El lector no interesado en problemas de clculo de esfuerzos en placas y lminas debe a a

46

saber que el tratamiento que se da a los problemas de exin de placas en el MEF es aplicao ble a otros problemas f sicos en cuya formulacin fuerte aparezcan ecuaciones diferenciales o de cuarto orden.

6.2

Formulacin C 1 en elementos viga o

La llamada formulacin C 1 para elementos nitos tipo viga corresponde a la teor clsica o a a a de exin de vigas segn Euler-Bernoulli. Esta teor es aplicable a vigas esbeltas1 y est o u a basada en la hiptesis de que durante la deformacin de la viga las secciones transversales o o permanecen planas y normales a la directriz deformada (hiptesis de Navier-Bernoulli)2 . o A partir de esta hiptesis, la deformacin de la viga queda denida por el campo w de o o echas o desplazamientos transversales de su directriz. La hiptesis de Navier-Bernoulli es una suposicin sobre la forma en que se produce la o o deformacin de la viga, que permite reducir la dimensionalidad del problema real. Para o conocer qu tensiones y deformaciones aparecen en todos los puntos de la viga basta con e estudiar los desplazamientos transversales de una lnea (directriz) trazada a lo largo de la misma. Se pasa de un problema elstico tridimensional, denido en todo el volumen de la a viga, a un problema unidimensional mucho ms sencillo, denido slo en su directriz. a o f h h h h h h h h h h h h h h h h h hh

Q0
?

? ?

M0 x= 0

w 1

Dada la viga de longitud unidad representada en la gura, con M0 ; Q0 2 < y f : ]0; 1[! <, la forma fuerte correspondiente a la teor clsica de exin se plantea como a a o sigue. Encontrar la funcin de echas o desplazamientos transversales w : [0; 1] ! <, o de modo que se cumplan las ecuaciones siguientes: 1. Ecuaciones de compatibilidad. En el contorno: w(1) = 0
1 2

En la prctica, vigas con una relacin canto/luz inferior a 1/5, aproximadamente. a o Esta hiptesis equivale a despreciar la deformacin por cortante. o o

47

w;x (1) = 0 y en el dominio, curvatura = w;xx en

2. Ecuaciones constitutivas o de comportamiento del material, M momento = EI con E; I 2 < 3. Ecuaciones de equilibrio. En el contorno: M(0) = M0 Q(0) = Q0 y en el dominio, Q cortante = d2 M =f dx2 dM dx en en en

Por otro lado, la forma dbil del problema puede enunciarse de la manera siguiente. Dados e los espacios de funciones: S = fw 2 H 2 () j w(1) = 0; w;x (1) = 0g ~ ~ ~ y V = fw 2 H 2 () j w(1) = 0; w;x (1) = 0g S encontrar una funcin w 2 S tal que: o (M w ; w ) = (w; f) + Q0 w(0) M0 w;x (0) 8 w 2 V (6.1)

e El desarrollo de elementos nitos de clase C 1 se hace en base a esta forma dbil del problema, siguiendo los pasos de la \receta" o sistemtica del MEF (seccin 4.2). Los dos a o primeros pasos ser an: 1. Se divide el dominio en subdominios o elementos, =
S
i

i .

Dentro de cada subdominio i se dene un sistema local de coordenadas, con coordenada local : 1 1 x = x1 (1 ) + x2 (1 + ) 2 2 o 2 = (x x1 ) 1 x2 x1 48

nodo 1
x

nodo 2
x -

-1

+1

La transformacin de coordenadas anterior resulta en que: o J= dx L x1 x2 = + = d 2 2 2

donde L es la longitud del elemento. 2. Se aproxima la echa o funcin w dentro de cada elemento como: o wh = A + B + C 2 + D 3 (6.2)

Ntese que se utiliza una aproximacin cbica ya que, para que sea wh 2 H 2 () o o u (derivada segunda de cuadrado integrable en ), se requiere que la primera derivada e de wh (giro) sea continua a travs de las fronteras de los elementos3 . Esto se consigue tomando como grados de libertad no slo las echas de los nodos, sino tambin los o e giros. Como los nodos son compartidos por los elementos adyacentes, se obtiene as la continuidad de giros al pasar de elemento a elemento. Se tienen entonces cuatro grados de libertad por elemento, una echa y un giro por nodo, lo que da lugar a una funcin cbica (cuatro coecientes) para la aproximacin de la echa dentro de o u o cada elemento. La aproximacin anterior para la echa, corresponde a la siguiente aproximacin o o para el giro: 4C 6D 2 2 h w;x = B + + L L L La particularizacin en los nodos de las expresiones anteriores para la echa y o el giro proporciona las ecuaciones necesarias para determinar los coecientes A, B, C y D en funcin de los grados de libertad o movimientos nodales: o AB +C D 2B 4C 6D + L L L A+B +C +D 2B 4C 6D + + L L L = w1 (echa en nodo 1) = 1 (giro en nodo 1) = w2 (echa en nodo 2) = 2 (giro en nodo 2)

3 La continuidad de la primera derivada de la variable de campo es lo que se conoce como continuidad C 1 en el lenguaje del MEF.

49

Lo que da: L 1 L 1 w1 + w2 + 1 2 2 2 8 8 3 L L 3 B = w1 + w2 1 2 4 4 8 8 L L C = 1 + 2 8 8 1 L L 1 D = w1 w2 + 1 + 2 4 4 8 8 A = La sustitucin de los coecientes A, B, C y D en 6.2 proporciona las funciones o de forma del elemento, que resultan ser polinomios de Hermite:
h 9 8 > w1 > > > > = i< >
1

wh =

1 2

3 + 1 3 4 4

L 8 (1

2 + 3 )

1 2

+ 3 1 3 4 4

L 8 (1

+ 2 + 3 )

= N ae Para entrar en las integrales de la forma dbil 6.1 y obtener la matriz de rigidez e del elemento hay que calcular tambin las curvaturas y los momentos: e h = L N ae = [
h d2 ] N ae = dx2
6 L2 6 1 3 L2 L + L 1 L 3 + L

> > w2 > > > > ; :

ae = B ae

M h = EI = EI B ae = D B ae De acuerdo con los puntos siguientes de la sistemtica del MEF4 , la sustitucin de las a o aproximaciones locales a w, y M en 6.1 permite construir la matriz de rigidez y el vector de cargas de cada elemento. Ntese que las funciones wh 2 S h , obtenidas por agregacin o o de las aproximaciones locales, son funciones cbicas a trozos. u

6.3

Formulacin C 1 en elementos placa o

La formulacin C 1 para elementos nitos tipo placa puede entenderse como una geneo ralizacin de la formulacin C 1 para elementos viga. Corresponde a la teor clsica de o o a a a exin de placas de Love-Kirchho. Esta teor es aplicable a placas delgadas5 y est o a basada en la hiptesis de que durante la deformacin las secciones transversales al plano o o medio de la placa permanecen planas y normales al plano medio deformado (hiptesis de o Love-Kirchho)6 . Con esta hiptesis, la deformacin de la placa queda denida por el o o campo w de echas o desplazamientos transversales del plano medio.
4 5

Ver seccin 4.2. o En la prctica, placas con una relacin canto/luz inferior a 1/15 1/20, aproximadamente. a o o 6 Esta hiptesis equivale a despreciar la deformacin por cortante. o o

50

La hiptesis de Love-Kirchho es un supuesto sobre la cinemtica de la deformacin o a o de la placa, que permite reducir la dimensionalidad del problema real. Para conocer qu e tensiones y deformaciones aparecen en todos los puntos de la placa basta con estudiar los desplazamientos transversales de un plano. Se pasa de un problema elstico tridimensional, a denido en todo el volumen de la placa, a un problema bidimensional ms sencillo, denido a slo en el plano medio de la placa. La variable de campo es pues el campo de echas o w : <2 ! <. La notacin ms corriente en las aplicaciones del MEF a problemas de exin de placas o a o es la que se da en la gura:

33 - -

Mxy My

Mxy y
3 y -

-- 3 3

Mx My
++

+ +

Mxy

Mx

33 - ? ?

Mxy

Las deformaciones generalizadas son las curvaturas:

8 2 > @ w > @x2 > > < 2 9 > > > > = > > > > ;

c=

Las tensiones generalizadas son los momentos por unidad de longitud: M=

8 9 > Mx > < =
y > : M > ; xy

@y2 > > > > @ 2w : 2 @x@y

Los desplazamientos generalizados son la echa y los giros. Los giros son derivadas primeras de la echa, por la hiptesis de Love-Kirchho: o
8 9 > w > < = > x > : ; y 8 > > > < > > > :

w @w @y
@w @x

w=

9 > > > = > > > ;

51

La forma dbil del problema, punto de partida de la sistemtica del MEF, es similar a la e a correspondiente a la de las vigas de Euler-Bernoulli. En este caso, el dominio de clculo a es bidimensional, <2 , y corresponde al plano medio de la placa. La forma dbil puede e enunciarse de la manera siguiente: Dados los espacios de funciones: S = fw 2 H 2 () j w = d en Sd g ~ ~ donde d son los movimientos impuestos en la parte Sd del contorno de la placa, y ~ V = f w 2 H 2 () j w = 0 en Sd g ~ encontrar una funcin w 2 S tal que: o (Mw ; cw ) = (q; w) + [t; w]St 8 w 2 V (6.3)

donde q son fuerzas generalizadas por unidad de supercie y t son esfuerzos impuestos sobre la parte St del contorno de la placa. Las aproximaciones locales a los campos que intervienen en 6.3, necesarias para la aplicacin de la \receta" del MEF, son las siguientes: o wh = N ae 8 > N >
h

> <

ch =

@y > > > : @N 8 @x 2 > @2 > @x > > < 2

@N

9 > > > =

Mh = donde, para placas istropas: o

@y 2 > > > > > > > : 2 @2 > ; @x@y 8 9 > Mx > < = y > : M > ; xy

> > > ; 9 > > > > =

ae

N ae = LN ae = B ae

= D c = DB ae

1 6 D= 4 1 12(1 2 ) 0 0 Et3

0 0
1 2

3 7 5

siendo E 2 < el mdulo de elasticidad; 2 <, el coeciente de Poisson; y t 2 <, el espesor. o

52

6.4

Dicultad de la formulacin C 1 en elementos placa o

Como en el caso de los elementos viga, la formulacin convencional del MEF para el o problema clsico de la exin de placas exige la continuidad de giros (continuidad C 1 ) a a o travs de las fronteras de los elementos. Si no existe esta continuidad, entonces se tiene e 2H que la echa wh 6 2 (). Desde el punto de vista de la sistemtica del MEF, si no existe a la continuidad de giros la integral que aparece en el primer miembro de 6.3 deja de ser la suma de integrales extendidas a los dominios elementales, ya que el integrando tiene un valor singular a lo largo de las fronteras entre elementos7 . La dicultad de la formulacin C 1 en elementos placa consiste en que, en el caso de o un dominio de clculo bidimensional, como es el plano medio de la placa, es imposible a conseguir la continuidad de giros a travs de las fronteras de los elementos si se utilizan los e tres grados de libertad naturales por nodo (la echa y los dos giros) y funciones polinmicas o para construir la aproximacin local al campo de echas wh . o u Ntese que, para asegurar la continuidad tanto de wh como de su derivada segn la o h h normal a una l nea de separacin entre dos elementos, w;n , los valores de wh y w;n tendr o an que estar denidos un vocamente a lo largo de esa l nea en funcin de los movimientos o nodales, ya que stos son compartidos por ambos elementos. e Vamos a ver que si se dene wh como un polinomio dentro de un elemento y se supone h h que tanto la echa wh como los giros w;x y w;y en la frontera dependen exclusivamente de los movimientos nodales (echa y dos giros), entonces se llega a una contradiccin. o En efecto, sea un elemento rectangular como el de la gura y denamos wh dentro del elemento como un polinomio en x e y. 4 3
u u

y
6 -

1 u x

u 2

Sea el lado 4 3 del elemento rectangular. La normal a este lado es la direccin y. o h es un polinomio en x e y dentro del elemento, a lo largo de este lado Entonces, como w se tendr, al ser y constante: a wh = A1 + A2 x + A3 x2 + : : : h w;y = B1 + B2 x + B3 x2 + : : :
La curvatura adopta un valor innito en las fronteras. Este valor innito, aunque se produzca en una zona de rea nula, da lugar a una contribucin innita a la energ de deformacin de la placa representada a o a o por el primer miembro de 6.3.
7

53

h con el nmero preciso de coecientes en cada expresin para deteminar wh y w;y en funcin u o o de los parmetros o variables nodales asociadas al lado 4 3. As si slo hay dos nodos, a , o est permitido que wh sea un polinomio de tercer grado a lo largo de 4 3, ya que se a h conocen cuatro parmetros: wh y w;x en los nodos 3 y 4. Del mismo modo, slo est a o a h tenga variacin lineal, ya que unicamente se conocen dos parmetros: o a permitido que w;y h w;y en los nodos 3 y 4. Se puede realizar un ejercicio similar a lo largo del lado 4 1:

wh = C1 + C2 y + C3 y2 + C4 y 3 h w;x = D1 + D2 y y se tiene entonces que, de acuerdo con nuestra suposicin: o

h 1. A lo largo de 4 3, w;y slo depende de los parmetros nodales del lado 4 3; con o a h unicamente depende de los parmetros nodales del lado 4 3. a lo que w;yx h 2. A lo largo de 4 1, w;x slo depende de los parmetros nodales del lado 4 1; con o a h unicamente depende de los parmetros nodales del lado 4 1. a lo que w;xy

3. En el nodo 4, para valores arbitrarios de los parmetros nodales, no se cumplir que a a h h w;yx = w;xy , lo cual es una contradiccin si wh se dene como un polinomio en x e o y dentro del elemento. Es decir, si wh es un polinomio dentro del elemento, los valores de las derivadas segn la u direccin normal a los lados (giros segn la normal) no pueden depender exclusivamente o u de los parmetros nodales del lado correspondiente y, por tanto, no est asegurada la a a continuidad de giros a travs de las fronteras entre elementos8 . e En consecuencia, en problemas de exin de placas formulados segn la teor clsica o u a a es muy dif construir funciones de forma N que cumplan los requisitos de continuidad cil C 1 que les impone la formulacin convencional del MEF. o

6.5

Ejemplo: elemento triangular de 3 nodos

Para ilustrar los conceptos anteriores, se presenta en esta seccin un elemento sencillo utio lizado en problemas de exin de placas delgadas9 . Se trata de un elemento no conforme, o esto es, que no cumple los requisitos de continuidad C 1 . Sin embargo, esto no se consider o obstculo para que el elemento se emplease en la prctica. La utilizacin de elementos a a o no conformes en problemas de exin de placas ha sido bastante habitual, a pesar de que o violan una de las condiciones de la formulacin convencional del MEF. Una justicacin o o de esta prctica se dar en los cap a a tulos siguientes, al hablar de la prueba de la parcela. El elemento que se presenta como ejemplo tiene 9 grados de libertad, tres por nodo, que se organizan como:
Lo unico que comparten dos elementos vecinos es el valor de los parmetros nodales de los nodos a comunes. 9 O.C. Zienkiewicz. El Mtodo de los Elementos Finitos, Revert, 1982. e e
8

54

ae =

8 > > > > > > > > > > > > > > < > > > > > > > > > > > > > > :

w1 x1 y1 w2 x2 y2 w3 x3 y3

9 > > > > > > > > > > > > > > = > > > > > > > > > > > > > > ;

Las funciones de forma se denen, como es habitual en los elementos nitos triangulares, en coordenadas de area L1 , L2 y L3 . Las coordenadas de area o coordenadas triangulares de un punto interior P son las siguientes: a rea 4P23 a rea 4123 a rea 4P31 a rea 4123 a rea 4P12 a rea 4123 3 S E
u

L1 = L2 = L3 =

y
6

ES E S E S EP S Z S Z Z S u X Z S 1 X X X X X X X Z S X XZ S XZ u S X 2

o tambin se pueden denir como: e x = L1 x1 + L2 x2 + L3 x3 y = L1 y1 + L2 y2 + L3 y3 1 = L1 + L2 + L3 donde (xi ; yi ); i = 1; 2; 3, son las coordenadas cartesianas de los vrtices del tringulo. e a En este elemento hay nueve funciones de forma para la interpolacin de la echa wh , o correspondientes a los nueve grados de libertad. Estas funciones de forma son: 55

Nt =

8 > > L1 + L2 L2 + L2 L3 L1 L2 L1 L2 > 1 1 2 3 > > > > b (L2 L + 1 L L L ) b (L2 L + 1 L L L ) > 3 1 2 2 > 1 3 2 1 2 3 2 1 2 3 > > > > c (L2 L + 1 L L L ) c (L2 L + 1 L L L ) > 3 1 2 > 2 1 3 2 1 2 3 2 1 2 3 > > > > > L2 + L2 L3 + L2 L1 L2 L2 L2 L2 > 2 2 3 1 <
1 1 2 2 > b1 (L2 L3 + 2 L1 L2 L3 ) b3 (L2 L1 + 2 L1 L2 L3 ) > > > > c1 (L2 L3 + 1 L1 L2 L3 ) c3 (L2 L1 + 1 L1 L2 L3 ) > 2 2 > 2 2 > > > 2 L + L2 L L L2 L L2 > L3 + L3 1 > 3 1 3 2 3 2 > > > > > b2 (L2 L1 + 1 L1 L2 L3 ) b1 (L2 L2 + 1 L1 L2 L3 ) > 3 3 > 2 2 > > : 1 1 2 2

9 > > > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > > > ;

c2 (L3 L1 + 2 L1 L2 L3 ) c1 (L3 L2 + 2 L1 L2 L3 )

con b1 = y2 y3 , c1 = x3 x2 , etc. Derivando estas nueve funciones de forma con respecto a x e y se obtienen las funciones de interpolacin de los giros. Derivando nuevamente se obtienen las funciones de o interpolacin de las curvaturas (matriz B). El clculo de todas estas derivadas se hace o a teniendo en cuenta que:

@ @x @ @y

@ @ @ 1 + b2 + b3 ) (b1 @L2 @L3 24 @L1 @ @ @ 1 + c2 + c3 ) (c1 @L2 @L3 24 @L1

donde 4 es el area del tringulo. a

6.6

Formulacin C 0 en elementos viga o

Utilizando la teor clsica de exin de placas existen dicultades serias para obtener a a o funciones de forma N que cumplan los requisitos de la formulacin convencional del MEF, o esto es, que tengan continuidad C 1 . En los a~os 70 y principios de los 80 se concibieron soluciones muy ingeniosas para, sin n abandonar la teor clsica, obtener funciones de forma con continuidad C 1 . Estas solua a ciones dieron lugar, en general, a elementos de formulacin complicada, con prestaciones o no muy brillantes y cuyo uso, en general, ha sido abandonado hoy en da. Al mismo tiempo se empez a pensar que, si la dicultad nac de las hiptesis de la o a o teor clsica de exin de placas, la dicultad pod sortearse formulando el problema de a a o a acuerdo con otra teor En consecuencia, se empezaron a utilizar teor de la exin que a. as o hasta entonces hab tenido slo una importancia marginal. Se trata de teor en las que an o as se tiene en cuenta la deformacin por corte: la teor de exin de vigas de Timoshenko o a o (1921) y la teor de exin de placas de Reissner-Mindlin (1944). Cuando se emplea la a o 56

formulacin convencional del MEF para abordar el problema de la exin planteado de o o acuerdo con una de estas teor slo se requiere continuidad C 0 para el campo de echas. as, o Como contrapartida aparecen los giros como campos independientes y para stos tambin e e 0 se requiere continuidad C . En conjunto, la formulacin de elementos nitos se simplica o a a. mucho, lo que ha hecho que sean estos elementos de tipo C 0 los ms utilizados hoy en d 0 Para introducir los conceptos bsicos, en esta seccin se plantea la formulacin C para a o o elementos viga. Estos mismos conceptos se extienden a placas en la seccin siguiente. o La formulacin C 0 para elementos nitos tipo viga corresponde a la teor de exin de o a o vigas de Timoshenko. Esta teor fue introducida en 1921 para el estudio de vibraciones a transversales de barras prismticas, con objeto de mejorar las predicciones de la teora a clsica en problemas de impacto o de vibraciones de alta frecuencia. a La teor est basada en la hiptesis de que durante la deformacin de la viga las a a o o secciones transversales permanecen planas, pero no normales a la directriz deformada. Con esta hiptesis, se necesitan dos campos independientes para denir la deformacin de o o la viga: el campo w de echas o desplazamientos transversales de su directriz y el campo de rotaciones de la normal a la directriz. La rotacin de la directriz se escribe como: o dw = rotacin de la directriz = o dx donde es la distorsin por cortante. o

f h h h h h h h h h h h h h h h h h hh

Q0
?

? ?

M0 x= 0

w 1

Dada la viga de longitud unidad representada en la gura, con M0 ; Q0 2 < y f : ]0; 1[! <, la forma fuerte del problema correspondiente a la teor de Timoshenko se a plantea como sigue. Encontrar las funciones w : [0; 1] ! < y : [0; 1] ! <, de modo que se cumplan las ecuaciones siguientes: 1. Ecuaciones de compatibilidad; en el contorno w(1) = 0 (1) = 0 57

y en el dominio, deform. de corte = w;x curvatura = ;x en en

2. Ecuaciones constitutivas o de comportamiento del material, Q M con E; I; G; k; A 2 < 3. Ecuaciones de equilibrio; en el contorno M(0) = M0 Q(0) = Q0 y en el dominio, Q = dQ dx dM dx en en cortante = GkA momento = EI en en

= f

Por otro lado, la forma dbil del problema puede enunciarse ahora de la manera siguiente. e Dados los espacios de funciones: Sw S Vw V = fw 2 H 1 () j w(1) = 0g ~ ~ ~ 2 H 1 () j (1) = 0g ~ = f = fw 2 H 1 () j w(1) = 0g Sw = f 2 H 1 () j (1) = 0g S

encontrar las funciones (w; ) 2 Sw S tales que: (M ; ) + (Qw; ; w; ) = (w; f) + Q0 w(0) M0 (0) 8 w 2 Vw V (6.4) e El desarrollo de elementos nitos de clase C 0 se hace en base a esta forma dbil del problema, siguiendo los pasos de la \receta" del MEF (seccin 4.2). Los dos primeros o pasos ser an: 1. Se divide el dominio en subdominios o elementos, =
S
i

i .

Dentro de cada subdominio i se dene un sistema local de coordenadas, con coordenada local : 1 1 x = x1 (1 ) + x2 (1 + ) 2 2 58

o =

2 (x x1 ) 1 x2 x1 nodo 2
x -

nodo 1
x

-1

+1

La transformacin de coordenadas anterior resulta en que: o J= dx x1 x2 L = + = d 2 2 2

donde L es la longitud del elemento. 2. En el caso de un elemento de dos nodos, se aproxima la echa y el giro dentro de cada elemento como:
(

w=

"

1 2 (1

) 0

1 2 (1

0 )

1 2 (1 +

)
1 2 (1

0 + )

8 9 > > # > w1 > > < > =

1

> w2 > > > > > : ;

= N ae

La eleccin anterior asegura la continuidad de echas y giros a travs de las fronteras o e de los elementos. Para entrar en las integrales de la forma dbil 6.4 y obtener la matriz de rigidez e del elemento hay que calcular tambin las deformaciones y tensiones generalie zadas:
(

e=

"

0 @ @x
(

@ @x

N ae =
"

"

0
1 L

1 L 0 1 1 2 (1 ) L

1 2 (1

1 L

+ )

8 9 > > # > w1 > > < > =

1

s=

M Q

EI 0 0 G kA

> w2 > > > > > : ;

= B ae

B ae = D B ae

De acuerdo con los puntos siguientes de la sistemtica del MEF10 , la sustitucin de las a o aproximaciones locales en 6.4 permite construir la matriz de rigidez y el vector de cargas de cada elemento.
10

Ver seccin 4.2. o

59

6.7

Formulacin C 0 en elementos placa o

La formulacin C 0 para elementos nitos tipo placa corresponde a la teor de exin de o a o placas de Reissner-Mindlin. Esta teora es aplicable a placas gruesas o poco esbeltas11 y est basada en la hiptesis de que en el proceso de deformacin las secciones transversales a o o al plano medio de la placa permanecen planas, aunque no necesariamente normales al plano medio deformado. Con esta hiptesis, la deformacin de la placa queda denida por o o el campo w de echas o desplazamientos transversales del plano medio y el campo de giros (x ; y ) de las bras inicialmente normales al plano medio. La variable de campo es pues el campo de desplazamientos generalizado w : ! <3 : w=
8 9 > w > < = > x > : ; y

Los campos de deformaciones generalizadas y de tensiones generalizadas son, respectivamente:

8 > > > > > <

e =

s =

> > > > > : 8 > Mx > > > M > y <

y;x x;y y;y x;x w;x y w;y + x

9 > > > > > = > > > > > ;

9 > > > > > = > > > > > ;

xy > > Q > > x > :

Qy

siendo s = D e, donde D es una matriz de constantes elsticas. a La forma dbil del problema, punto de partida de la \receta" del MEF, es similar a e la de las vigas de Timoshenko. Ahora el dominio de clculo <2 es bidimensional y a corresponde al plano medio de la placa. La forma dbil es la siguiente: e Dados los espacios de funciones: S = fw 2 H 1 () j w = d en Sd g ~ ~ y V = fw 2 H 1 () j w = 0; en Sd g encontrar w 2 S tal que: (sw ; ew ) = (w; q) + [w; t]St
11

8 w 2 V

En la prctica, placas con una relacin canto/luz inferior a 1/8 1/10, aproximadamente. a o o

60

con la notacin habitual. o El desarrollo de elementos nitos de clase C 0 se hace en base a esta forma dbil del e problema, siguiendo la sistemtica del MEF. La eleccin de funciones de forma es ahora a o ms sencilla, ya que no tienen que cumplir requisitos de continuidad C 1 . Sin embargo a la formulacin del elemento debe hacerse con cuidado, para evitar el problema conocido o como bloqueo de cortante. El mecanismo de este bloqueo se traduce en una respuesta excesivamente r gida de los elementos cuando la relacin canto/luz de la placa es peque~a o n (placas esbeltas). Este y otros problemas que presenta la formulacin convencional del o MEF, se describen en el cap tulo siguiente.

61

Cap tulo 7

Tecnolog de elementos (III) a

7.1 Introduccin o

En los cap tulos anteriores se han presentado las l neas generales de la formulacin cono vencional del MEF y sus fundamentos matemticos. La formulacin convencional se caa o racteriza por basarse en la aproximacin de una sola variable de campo (escalar, vectorial o o tensorial) a partir de la cual se obtienen de forma fuerte el resto de campos incgnita o del problema. La formulacin convencional impone tres requisitos a las funciones de forma o funciones o de aproximacin local de la variable de campo incgnita bsica: o o a Suavidad dentro de los elementos o dominios locales. Continuidad a travs de las fronteras entre elementos. e Completitud de la aproximacin hasta un orden que depende del problema de campo o que se resuelve. El cumplimiento de estas tres condiciones y el seguimiento punto por punto de la \receta" o sistemtica del MEF introducida en el Cap a tulo 5, garantiza la convergencia de las soluciones aproximadas hacia la solucin exacta del problema de campo. Desde el punto o de vista matemtico, el mtodo se convierte en un caso particular de los procedimientos a e de Rayleigh-Ritz-Galerkin y se cumplen las condiciones sucientes para la convergencia. Sin embargo, en los programas actuales de clculo basados en el MEF resulta muy a corriente encontrar elementos cuya formulacin no sigue estrictamente las l o neas de la formulacin convencional: se viola alguna de las condiciones sucientes de convergencia o o se ejecuta de modo irregular alguno de los pasos de la sistemtica del MEF. a Los elementos anteriores comenten lo que hacia mediados de los a~os 70 se baun tiz con el nombre de cr o menes variacionales: la formulacin del elemento se aparta o heur sticamente de la formulacin convencional, pero el elemento funciona muy bien, ino cluso mejor que los elementos leg timos equivalentes. Desde el punto de vista matemtico, a el procedimiento deja de ser un mtodo de Rayleigh-Ritz-Galerkin y, en principio, su e convergencia no est garantizada. a 62

Hoy en d se conoce en la mayor de los casos la justicacin matemtica de aquellos a a o a cr menes variacionales. Sin embargo, los elementos criminales empezaron a utilizarse sin que hubiera una garant matemtica clara de su aceptabilidad. Simplemente, proporcioa a naban a juicio de los analistas mejores soluciones con un menor coste de clculo. a Para concretar las observaciones anteriores, se pueden mencionar dos ejemplos clsicos a de este modo de proceder: Elementos no conformes. En este tipo de elementos el campo incgnita bsica no es o a continuo a travs de las fronteras entre los elementos, esto es, se viola la segunda de e las condiciones sucientes de convergencia mencionadas ms arriba. Como resultado a de esta falta de continuidad, alguna de las integrales que aparecen en la forma dbil e 1 del problema se hace singular en las fronteras entre los elementos . A pesar de ello, se opera como si dicha integral, extendida a todo el dominio de clculo, pudiera a todav descomponerse en suma de integrales extendidas a los distintos elementos. a Integracin reducida. Se sigue puntualmente la sistemtica del MEF, pero para o a evaluar numricamente las integrales necesarias, se emplean menos puntos de los e precisos para integrar exactamente los polinomios que aparecen como funciones subintegrales. Dentro del desarrollo del MEF, desde mediados de los a~os 60 los investigadores eran n conscientes de que en muchos casos se obten mejores soluciones, con menos esfuerzo de an clculo, si se comet alguno de estos cr a a menes variacionales. Al comienzo de los 70 se vio que en determinados problemas era prcticamente ima prescindible cometer un \crimen" si se quer llegar a soluciones m a nimamente aceptables (p.ej. en problemas de plasticidad o de elasticidad con materiales incompresibles). En esos problemas la formulacin convencional proporcionaba malos resultados en muchos casos. o Desde el punto de vista matemtico, no se tuvo una forma de explicar por qu los a e elementos \criminales" funcionaban bien hasta el nal de los a~os 70, cuando empez a n o verse que la mayor de estos elementos respond a formulaciones variacionales con varios a an campos independientes. Todav hoy la convergencia de los elementos que responden a un a principio variacional mixto o h brido de esta clase es motivo de investigacin. o Desde el punto de vista de los usuarios del MEF, hasta principios de los a~os 80, y n aun hoy en d en muchos casos, casi la unica herramienta disponible para justicar el a empleo de los elementos que se apartan de la formulacin convencional ha sido la prueba o de la parcela. Dicha prueba fue introducida por Irons en 1965 y, a lo largo del tiempo, su signicado matemtico real ha sido objeto de polmica. Sin embargo, durante los a~os 70 a e n su satisfaccin era aceptada como condicin necesaria y suciente para que un elemento o o diera lugar a soluciones convergentes, independientemente de cmo hubiese sido formulado. o Dentro de este cap tulo se va a presentar primero una muestra de los problemas en los que la formulacin convencional del MEF proporciona malos resultados. Luego se hablar o a de los remedios heur sticos o cr menes variacionales que se desarrollaron para sortear estos
1 En un problema de elasticidad, el signicado f sico de la singularidad es que la densidad de energ de a deformacin aproximada se hace innita a lo largo de las fronteras entre los elementos. o

63

problemas y, tambin, de la prueba de la parcela, como justicacin para el empleo de e o los elementos resultantes. En el cap tulo siguiente se introducir la forma matemtica a a de justicar los remedios heur sticos en base a principios variacionales con varios campos independientes.

7.2

Elementos no conformes

Se muestran en esta seccin ejemplos tomados de dos areas de aplicacin del MEF en las o o que el uso de elementos no conformes resulta corriente.

7.2.1

Elasticidad bidimensional

En problemas de elasticidad plana los elementos cuadrilteros convencionales de cuatro a nodos proporcionan una aproximacin muy pobre si la exin es el esfuerzo dominante. o o Cuando se utilizan estos elementos en problemas de exin se obtiene una respuesta muy o r gida y es necesaria una discretizacin extraordinariamente na para que los resultados o sean m nimamente aceptables. Este problema era conocido desde los comienzos de la aplicacin del MEF a problemas o de elasticidad. En 1973 Wilson mostr que el comportamiento en exin del elemento o o convencional cuadriltero de cuatro nodos pod mejorarse mucho si se introduc dos a a an nuevas funciones de forma asociadas a grados de libertad internos. Es decir, si se obtena la aproximacin local al campo de desplazamientos como: o uh =
4 X i=1

Ni ae + i

2 X i=1

Ni i

y anlogamente para la otra componente v h del desplazamiento. En la relacin anterior a o Ni son las cuatro funciones de forma bilineales convencionales, ae son los desplazamientos i nodales, Ni son las dos funciones de forma a~adidas y i son los dos grados de libertad n internos asociados. Las dos nuevas funciones estn pensadas para que, junto con las funciones convencioa nales, den lugar a modos de deformacin aproximados que sean capaces de reproducir un o estado tensional de exin pura. En las coordenadas (; ) del dominio cuadrado estndar o a las funciones a~adidas son: n
N1 = (1 2 ) N2 = (1 2 )

El inconveniente es que las dos nuevas funciones de forma se asocian a grados de libertad internos del elemento, esto es, a grados de libertad no compartidos con los elementos vecinos, con lo cual no existe necesariamente continuidad del campo de desplazamientos aproximado al pasar de un elemento a otro. Esto da lugar a la violacin de una de las o condiciones que la formulacin convencional impone a las funciones de forma. o 64

A pesar de lo anterior, puede formularse el elemento siguiendo la sistemtica del MEF a de construccin de la matriz de rigidez elemental y el vector elemental de cargas. Se genera o as un elemento no conforme, pero los resultados que se obtienen con l son espectacu e larmente buenos, mucho mejores que los obtenidos con el elemento leg timo, incluso en mallas relativamente gruesas. En consecuencia, este elemento basado en la idea de Wilson y elaborado un poco ms adelante por Taylor, se incorpor rpidamente a la mayor de a o a a los programas comerciales de clculo2 . a

7.2.2

Flexin de placas o

En la teor clsica de exin de placas delgadas, el estado de deformacin de la placa se a a o o representa completamente a partir del campo de echas o desplazamientos laterales w del plano medio de la placa. La formulacin convencional del MEF aplicada a esta clase de problemas requiere la o continuidad a travs de las fronteras entre elementos de los giros o primeras derivadas de la e aproximacin wh al campo de echas. Es lo que en el Cap o tulo 6 se denomin continuidad o 1 del campo aproximado wh . C Determinar funciones de forma que cumplan las condiciones de continuidad C 1 es relativamente complejo. De hecho, como se vio en el Cap tulo 6, esto no es posible si se utilizan funciones polinmicas sencillas y slo tres grados de libertad por nodo, la echa o o y los dos giros. Debido a la importancia tecnolgica de los problemas de clculo de placas, se ha dedio a cado un esfuerzo investigador muy importante a resolver esta dicultad. Mantenindose e dentro de la teor clsica de exin de placas, se han propuesto distintas soluciones: a a o Hiperelementos. Se trata de elementos que cumplen los requisitos de continuidad C 1 en base a la introduccin de grados de libertad no f o sicos en los nodos (p.ej. derivadas de la echa de orden superior al primero). Funciones de forma denidas a trozos. Se hace cumplir la continuidad C 1 a travs e de las fronteras de los elementos en base a disminuir la suavidad de las funciones de aproximacin de la echa en el interior de los mismos. Para ello se divide cada eleo mento en subelementos y se impone la continuidad C 1 tanto a travs de las fronteras e interiores como exteriores de estos subelementos. Elementos no conformes. En realidad no es una solucin, se trata simplemente de o utilizar elementos en los que no se cumple la continuidad de giros a travs de las e fronteras de los elementos. Las dos primeras soluciones cumplen los requisitos de la formulacin convencional, pero o dan lugar a elementos de formulacin muy compleja y que no siempre funcionan mejor que o
Este elemento se conoce con el nombre de elemento de Wilson-Taylor o elemento de modos incompatibles, y constituye uno de los primeros elementos de altas prestaciones. Elementos de altas prestaciones son aquellos capaces de suministrar resultados con un nivel de aproximacin suciente desde el punto de o vista ingenieril utilizando mallados relativamente gruesos.
2

65

los elementos no conformes. De hecho, los hiperelementos y los elementos con funciones de forma denidas a trozos han sido desterrados de la prctica actual. a

7.3

Bloqueo por deformacin isocrica o o

A principios de la dcada de los 70 se comenzaron a detectar dicultades en la aplicae cin de la formulacin convencional del MEF a problemas de plasticidad o de elasticidad o o con materiales poco compresibles. Estos problemas tienen en comn el que la deformau cin es isocrica, esto es, tiene lugar manteniendo el volumen del material prcticamente o o a constante3 . Las dicultades eran muy serias. En problemas de plasticidad no se calculaban correctamente las cargas l mite, cuya existencia se demostraba tericamente. Bien se obten o an cargas l mite que exced el nivel superior calculado mediante el teorema cinemtico del an a clculo plstico, o bien no se obten ninguna carga l a a a mite, esto es, se pod cargar indea nidamente la estructura o medio continuo sin llegar al colapso. Adems, las soluciones a apenas mejoraban con el renamiento de las mallas. En problemas de elasticidad con materiales poco compresibles se obtenan soluciones muy r gidas, sin apenas deformacin o con deformaciones muy inferiores a las tericas o o o experimentales. Las soluciones tampoco mejoraban perceptiblemente con el renamiento de los mallados, esto es, no se produc convergencia a efectos prcticos. a a Para designar estas situaciones se acu~o el trmino de bloqueo o atenazamiento de man e lla4 . Pareci en aquel momento que el MEF no pod proporcionar resultados aceptables o a en dos areas con aplicaciones industriales muy importantes. En 1974, en un art culo histrico que abri las puertas a muchas aplicaciones no lineales o o del MEF, Nagtegaal5 dio una explicacin a estos fenmenos de bloqueo y propuso una o o solucin que pasaba por el abandono de la formulacin convencional de los elementos. o o Desde entonces hasta hoy se han propuesto otras soluciones al problema, que se han utilizado luego con ms o menos xito. Como resultado de este esfuerzo rara vez se a e emplean hoy en d elementos de formulacin convencional en las aplicaciones del MEF a o dentro del ambito no lineal. Con posterioridad se han descubierto nuevos tipos de bloqueo que aparecen en otras clases de problemas. Hay que mencionar los dos siguientes: Bloqueo de cortante6 . Aparece en los elementos placa de tipo C 0 cuando se utilizan para estudiar la exin de placas muy esbeltas. o Bloqueo de membrana7 . Aparece en los elementos lmina cuando se utilizan para a
No todos los modelos de plasticidad dan lugar a deformaciones plsticas isocricas, pero s los modelos a o aplicados habitualmente al estudio de plasticidad de metales o al ujo plstico de suelos arcillosos. a 4 El fenmeno se conoce como \volumetric locking" en la literatura anglosajona. o 5 J.C. Nagtegaal, D.M. Parks y J.R. Rice. On Numerically Accurate Finite Element Solutions in the Fully Plastic Range. Computer Methods in Applied Mechanics and Engineering, vol.4, pgs. 153-178, 1974. a Este art culo es uno de los ms citados de la historia del MEF. a 6 Conocido como \shear locking" en la literatura anglosajona. 7 Conocido como \membrane locking" en la literatura anglosajona.
3

66

estudiar la deformacin de lminas de simple o doble curvatura. o a Desde el punto de vista matemtico, estos dos bloqueos responden a la misma razn a o que el bloqueo volumtrico o por deformacin isocrica. e o o Lo caracter stico de los bloqueos es que no se produce un cambio perceptible de la solucin con el renamiento de la malla, o el cambio es tan peque~o que es como si no o n existiera a efectos prcticos. a La razn de que se produzcan los bloqueos est en la eleccin de las funciones de forma. o a o Una vez elegidas las funciones de forma, estn determinadas la formas o modos que puede a adoptar el campo de deformaciones aproximado dentro del elemento. Imagnese que dentro de esos modos predeterminados no hay ninguno que permita la deformacin del elemento o sin cambio de volumen, esto es, que para cualquier combinacin de movimientos nodales el o campo aproximado de deformacin que se obtiene dentro del elemento implica un cambio o de volumen. Entonces, si el material es de naturaleza incompresible (resistencia innita al cambio de volumen), el campo aproximado de deformacin dar lugar a la aparicin de o a o una resistencia considerable al movimiento de los nodos, lo que origina el bloqueo. Un razonamiento similar puede hacerse para el caso en que los modos de deformacin o 0 de un elemento placa de tipo C no incluyan modos sin deformacin de cortante (bloqueo o de cortante), o para el caso en que los modos de deformacin de un elemento lmina o a no incluyan modos de exin sin deformacin de la supercie media o deformacin de o o o membrana (bloqueo de membrana). En todos esos casos, la eleccin de las funciones de forma impone, a travs del como e portamiento del material, restricciones cinemticas al campo incgnita bsica que son a o a superiores a las f sicamente razonables. En los prrafos siguientes se trata de explicar el mecanismo por el que se produce a el bloqueo volumtrico en el elemento cuadriltero convencional de cuatro nodos. El e a mecanismo es idntico desde el punto de vista matemtico en los casos de bloqueo de e a cortante y bloqueo de membrana.

7.4

Bloqueo del elemento cuadriltero convencional a

Sea el elemento cuadriltero isoparamtrico convencional de cuatro nodos. Para simplicar a e el algebra, se supondr que se trata de un elemento cuadrado de dos unidades de lado, a esto es, que el elemento ocupa el dominio estndar isoparamtrico. a e

67

y (-1,1) 4 x
6

(1,1)
x3 -

1 x (-1,-1)

x2

(1,-1)

Se toma el origen de los ejes locales en el centro de areas del elemento. Entonces, cuando el elemento se utiliza para el estudio de problemas de elasticidad bidimensional, las dos componentes del campo de desplazamientos se aproximan dentro del elemento como: u = donde:
6 6 6 6 6 6 6 6 6 6 6 6 6 t N =6 6 6 6 6 6 6 6 6 6 6 6 6 4 2
1 4 (1 h

u v

= N ae

x)(1 y) 0

1 4 (1

+ x)(1 y) 0

1 4 (1

+ x)(1 + y) 0

1 4 (1

x)(1 + y) 0

7 7 7 7 1 (1 x)(1 y) 7 4 7 7 7 7 0 7 7 7 1 (1 + x)(1 y) 7 4 7 7 7 7 0 7 7 7 1 (1 + x)(1 + y) 7 4 7 7 7 7 0 7 5

ae =

8 9 > > > > u > > 1 > > > > > > > > > > > > > > v > > 1 > > > > > > > > > > > > > > u > > 2 > > > > > > > > > > > > > > < =

v2

1 4 (1

x)(1 + y)

> > > > > > > > > > > > > > > > > > > > > > > > > > :

> > u3 > > > > > > > > v3 > > > > > > > > u4 > > > > > > > v ;
4

>

Por otro lado, en un problema de deformacin plana, las componentes del campo de o deformaciones se aproximan como:
8 > "x > > < 6 6 6 6 6 =6 > 6 > 6 > 6 ; 6 4 9 > > > = 2
@ @x

0
@ @y

"y e = > "z > > : "xy

h

0 0
1 @ 2 @y

0
1 @ 2 @x

7 7 7 7 7 7 N ae = L N ae = B ae 7 7 7 7 5

68

donde la matriz B es:

2 6 6 6 6 6 6 6 6 6 6 4

1 (1 4 0 0

y)

0 1 (1 x) 4 0

1 4 (1

y)

0 1 (1 + x) 4 0
1 8 (1

1 4 (1

+ y) 0 0
1 4 (1

0 + x) 0
1 8 (1

1 (1 4 0 0
1 8 (1

+ y)

0 0

1 (1 x) 1 (1 y) 1 (1 + x) 8 8 8

y)

1 8 (1

+ x)

+ y)

x)

1 (1 + y) 8

7 7 7 7 1 (1 x) 7 4 7 7 7 7 0 7 5

Cada columna de la matriz anterior representa un posible modo de deformacin del o elemento. Sin embargo, las ocho columnas de la matriz no son linealmente independientes entre s slo hay cinco columnas linealmente independientes, lo que da lugar unicamente , o a cinco modos de deformacin independientes. Estos modos son: o
8 9 > 1 > > > > > < =

m1 =

Es importante darse cuenta de que, una vez seleccionadas las funciones de forma N, el elemento ya no puede deformarse arbitrariamente. La eleccin de las funciones de forma o que se ha hecho implica que en el interior del elemento el campo de deformaciones aproximado slo puede obtenerse como combinacin lineal de esos cinco modos de deformacin. o o o La matriz B puede descomponerse en suma de dos matrices Bv y Bd , la primera correspondiente a la parte volumtrica de la deformacin y la segunda a la parte desviadora. e o Estas dos matrices son:
8 9 > 1 > > > > > 1< 1 = Bv0 Bv = 3> 1 > > > > > : ;

0 > 0 > > > > > : 0 ;

m2 =

8 9 > 0 > > > > > < =

1 > 0 > > > > > : 0 ;

m3 =

8 9 > 0 > > > > > < =

0 > 0 > > > > > : 1 ;

m4 =

8 9 > y > > > > > < =

0 > 0 > > > 1 > > : ; 2x

m5 =

8 9 > 0 > > > > > < =

x > 0 > > > 1 > > : ; 2y

donde Bv0 es el vector la:

h
1+y 4 1+x 4

1y 4

1x 4

1+y 4

1+x 4

1y 4

1x 4

y Bd = B Bv A partir de esta descomposicin, resulta que en un problema de deformacin plana se o o cumple que: D Bv = D Bd = E Bv (1 2) E Bd (1 + )

69

siendo D la matriz de coecientes elsticos: a

6 6 6 6 6 E(1 ) 6 D= 6 (1 + )(1 2) 6 6 6 6 4 2

1
(1) (1)

(1)

(1) (1)

0 0 0
(12) (1)

3 7 7 7 7 7 7 7 7 7 7 7 5

1
(1)

1 0

y E es el mdulo de elasticidad y el coeciente de Poisson. o De esta manera la matriz de rigidez del elemento queda: Ke = (Bt DB)e = (fBt + Bt gDfBv + Bd g)e v d = (Bt DBv )e + (Bt DBv )e + (Bt DBd )e + (Bt DBd )e v d v d E E = (Bt Bv )e + (Bt Bd )e (1 2) v (1 + ) d = Kev + Ked al ser Bt Bv = Bt Bd = 0, ya que la suma de las tres primeras las de Bd es idnticamente e v d nula. Resulta entonces que la matriz de rigidez elemental puede descomponerse en dos partes, una correspondiente a la parte volumtrica de la deformacin y otra a la parte desviadoe o ra. Cada una de ellas es el producto de una matriz por un coeciente derivado de las propiedades elsticas del material. Ntese que el coeciente correspondiente a la parte a o volumtrica tiende a innito cuando el coeciente de Poisson tiende a 1 . e 2 Como consecuencia, si para simplicar imaginamos una malla con un solo elemento, el sistema global de ecuaciones del MEF: Ke ae = (Kev + Ked ) ae = fe tender a ser: a Kev ae = 0 cuando el coeciente de Poisson tiende a 1 , ya que fe es constante y Ked tiende a un valor 2 nito. En consecuencia, las unicas soluciones con desplazamientos nodales ae distintos de cero que admitir el sistema sern los autovectores de valor propio nulo de la matriz Kev . a a Para ver cules son estos autovectores se debe profundizar en la estructura de la matriz a Kev . Esta matriz es:

70

Kev = = =

E (Bt Bv )e (1 2) v E (Bt Bv0 )e 3(1 2) v0 X E Bt Bv0 wi v0 3(1 2) i

(7.1) (7.2) (7.3)

donde i var desde 1 al nmero de puntos de integracin necesarios, y wi es el peso a u o asignado a cada uno de esos puntos. Como se tiene que el vector la Bv0 puede descomponerse en:
i 1h 1 1 1 1 1 1 1 1 4 i yh 1 0 1 0 1 0 1 0 4 i xh 0 1 0 1 0 1 0 1 4

Bv0 = + +

resulta que la matriz Kev , que es de 8 8, tiene, como mucho, rango 3, ya que el valor de Bv0 en cualquier punto de integracin dentro del elemento puede obtenerse siempre como o combinacin lineal de slo tres vectores. o o As resulta que la matriz Kev tiene al menos cinco autovectores linealmente indepen, dientes con autovalor nulo. Una base del espacio de autovectores con autovalor nulo puede ser la siguiente8 :
8 > > > > > > > > > > > > < > > > > > > > > > > > > :

a1 =

1 0 1 0 1 0 1 0

9 > > > > > > > > > > > > = > > > > > > > > > > > > ;

a2 =

8 > > > > > > > > > > > > < > > > > > > > > > > > > :

0 1 0 1 0 1 0 1

9 > > > > > > > > > > > > = > > > > > > > > > > > > ;

a3 =

8 > > > > > > > > > > > > < > > > > > > > > > > > > :

0 0 0 1 1 1 1 0

9 > > > > > > > > > > > > = > > > > > > > > > > > > ;

a4 =

8 > > > > > > > > > > > > < > > > > > > > > > > > > :

0 0 1 0 1 1 0 1

9 > > > > > > > > > > > > = > > > > > > > > > > > > ;

a5 =

8 > > > > > > > > > > > > < > > > > > > > > > > > > :

0 0 0 0 1 0 1 0

9 > > > > > > > > > > > > = > > > > > > > > > > > > ;

Los tres primeros vectores corresponden a movimientos de slido r o gido del elemento, es decir, a movimientos de conjunto en los que no hay deformaciones. El cuarto, a un movimiento en el que la deformacin "xx es igual y de signo contrario a la deformacin o o "yy . Por ultimo, el quinto vector corresponde a una deformacin de corte puro. o
8

Puede comprobarse que todos estos vectores son ortogonales a Bv0 .

71

Modo a4

Modo a5

Cuando el material tienda a ser incompresible, el elemento slo podr deformarse segn o a u uno de estos dos ultimos modos. Cuando las condiciones de contorno y las acciones no se correspondan con una deformacin segn uno de estos modos, tendr lugar el bloqueo de o u a la malla.

7.5
7.5.1

Soluciones heur sticas a los problemas de bloqueo

Integracin reducida o

Ntese que si en 7.3 se emplea un solo punto de integracin para formar la matriz Kev , o o entonces su rango pasa de 3 a 1, es decir, aparecen dos nuevos modos (autovectores) segn u los cuales el elemento puede deformarse sin cambio de volumen. Si se utiliza el centro de areas del elemento como unico punto de integracin, los dos o nuevos modos son:
8 9 > 1 > > > > > 0 > > > > > > > > > > > > 1 > > > > > > < = 8 > > > > > > > > > > > > < > > > > > > > > > > > > :

a6 =

> 1 > > > > > > > 0 > > > > > > > > 1 > > > > > > > > : ;

a7 =

0 1 0 1 0 1 0 1

9 > > > > > > > > > > > > = > > > > > > > > > > > > ;

A A A A HH

Modo a6

Modo a7 H H

H H

Esta es una buena solucin para el bloqueo de este elemento: subintegrar la parte o volumtrica de la matriz de rigidez y mantener la integracin correcta o completa de la e o parte desviadora. Esta tcnica se conoce como integracin reducida selectiva, y da lugar a e o 72

un elemento relativamente robusto, potente y de uso general el problemas de plasticidad y elasticidad incompresible. Desde el punto de vista de la programacin y del ahorro de tiempo de clculo, interesa o a subintegrar tambin la parte desviadora de la matriz de rigidez del elemento. Esto da e lugar a un elemento completamente subintegrado o de integracin reducida propiamente o dicha. Se trata de un elemento muy barato, que proporciona tiempos de ejecucin signio cativamente menores que los del elemento convencional, pero que presenta problemas en cierta medida opuestos a los del bloqueo. Ahora existen en el elemento dos modos de deformacin, a6 y a7 a los que no se opone ninguna fuerza interna de respuesta del material. o Son los llamados modos de energa nula o modos de \hourglassing". El elemento puede deformarse segn esos modos sin producir ninguna tensin, como si se tratase de modos u o de slido r o gido. Los modos de energ nula pueden originar dicultades si son excitados por las acciones a impuestas y las condiciones de contorno. En esos casos se calculan ms o menos bien a las tensiones, pero los desplazamientos estn completamente mal y las deformadas son a impresentables. A pesar de este problema, los elementos de integracin reducida se utilizan mucho, o 9 sobre todo en programas de integracin expl o cita , donde es vital reducir el tiempo de ejecucin dedicado a los clculos en los elementos. o a Se han desarrollado tcnicas de control de \hourglassing", con un coste computacional e intermedio entre el de la integracin reducida y el de la integracin reducida selectiva. Se o o trata de tcnicas heur e sticas que se basan en introducir una rigidez articial que se opone a la deformacin de energ nula. A primera vista esto parece una solucin sumamente o a o articiosa, pero anando con cuidado la rigidez articial los elementos resultantes funcionan bien para un rango muy amplio de problemas. De hecho, esta clase de elementos se encuentran ampliamente representados en las bibliotecas de los programas comerciales de clculo y su uso est muy generalizado, aunque no les faltan detractores. a a En los programas comerciales, todos los elementos de integracin reducida incorpoo ran algn tipo de control de \hourglassing", sintonizado para el tipo de clculos que se u a supone que resuelve el usuario medio del programa. Fuera de este rango, el control de \hourglassing" puede dar lugar a soluciones articialmente r gidas o exibles.

7.5.2

Formulacin B o

El bloqueo se debe a que no es adecuada la matriz de interpolacin de deformaciones B o que se deriva de las funciones de forma seleccionadas N, por dar lugar a demasiadas restricciones sobre el movimiento. Puede entonces intentar resolverse la dicultad cambiando la matriz B por otra matriz B (lase b-barra) que no d lugar a tantas restricciones. e e La solucin consistir en utilizar en el clculo de la matriz de rigidez elemental una o a a matriz B, en lugar de la matriz \correcta" B. La idea es que si se tiene:
9

Vese el Cap a tulo 10.

73

B = Bv + Bd entonces, como es la parte Bv la que produce el bloqueo, el bloqueo podr evitarse a trabajando con una matriz B tal que: B = Bv + Bd El inconveniente es que en este caso se dejar de cumplir estrictamente la sistemtica a a del MEF, ya que, evidentemente, ser: a B 6L N = B = La clave de este procedimiento est en elegir acertadamente las funciones contenidas a culo propuso una solucin de este estilo. En la referencia 1 de o en Bv . Nagtegaal en su art la bibliograf recomendada pueden verse otras propuestas. Ntese que si en el elemento a o cuadriltero de la seccin anterior se toma para Bv el valor en el centro de areas de la a o matriz Bv , es decir: Bv = Bv (0; 0) = constante se obtiene el elemento de integracin reducida selectiva. De este modo, la formulacin B o o permite generalizar los \trucos" de la integracin reducida. o Las formulaciones B se conocen tambin en la literatura con el nombre de mtodos de e e proyeccin de deformaciones. o

7.6

La prueba de la parcela

La prueba de la parcela o \patch test" fue propuesta por Irons a mediados de los a~os n 60 como comprobacin de la formulacin de los elementos y como test de convergencia o o aplicable a elementos no conformes. Durante los a~os 70 el cumplimiento de la prueba se tom como condicin necesaria n o o y suciente de convergencia y permiti justicar el uso de elementos que se apartaban de o la formulacin convencional del MEF. Sin embargo, hacia el nal de los 70, una serie de o experimentos numricos sugirieron que la prueba podr no ser una condicin necesaria. e a o Adems, al poco tiempo, mediante un contraejemplo, alguien vio que la prueba podra no a proporcionar tampoco una condicin suciente. o A la vista de esos resulados, Irons y otros notables defensores de la utilidad de la prueba, la repensaron y llegaron a la conclusin de que, salvo en ciertas situaciones patolgicas, o o su cumplimiento s que era una condicin necesaria y suciente de convergencia para las o soluciones aproximadas obtenidas mediante el MEF. La polmica continu y el signicado matemtico real de la prueba sigue siendo objeto e o a de cierta discusin. El lector interesado puede encontrar ms detalles en las referencias 1 o a y 3 de la bibliograf recomendada. a 74

Lo que s parece claro es que la prueba de la parcela es una forma muy sencilla de comprobar la condicin de consistencia o completitud de la formulacin10 . Y que si la o o o formulacin es consistente y estable11 , entonces la formulacin es convergente. o En cualquier caso, el cumplimiento de la prueba de la parcela ha sido durante muchos a~os la unica justicacin al empleo de elementos que comenten crmenes variacionales n o en su formulacin. o La idea f sica de la prueba es muy sencilla. Pinsese por ejemplo en un problema de e elasticidad. Si se considera el dominio de clculo dividido cada vez en un nmero mayor a u de subdominios e , parece lgico pensar que el campo de desplazamientos u solucin o o exacta del problema podr aproximarse progresivamente mejor mediante una funcin que a o var linealmente dentro de cada subdominio. e En una situacin que tienda al l o mite de disminucin del tama~o de los subdominios, o n llegar un momento en que la solucin exacta no podr distinguirse de su aproximacin a o a o mediante una funcin lineal. Entonces, si una parcela o grupo de elementos nitos de o forma arbitraria es capaz de representar de forma exacta una variacin lineal del campo o incgnita bsica u dentro de la parcela, resulta razonable pensar que la aproximacin o a ~ o obtenida por el MEF tender a la solucin exacta al aumentar o renar la discretizacin. a o o La realizacin de la prueba de la parcela es muy sencilla en el contexto de un programa o de elementos nitos. Una forma clsica de la misma es la siguiente: a
x h h h h h h S h h h h h h h S h h h h x S A S A S A S A S A S A S A S A S x i!S ! H A ! C H C H A ! ! H H C ! ! A ! H H C A ! ! H H ! C A ! ! H C A H H ! ! x P P C H A P P H H C P P P P H AAx C H P P P P C P P P P C PC x

El cumplimiento de la condicin de consistencia signica que la formulacin discreta se aproxima tanto o o como se quiera a la formulacin continua de partida al aumentar el nivel de discretizacin. o o 11 Una formulacin es estable cuando peque~ as perturbaciones en los datos no producen grandes perturo n baciones en la solucin aproximada. o

10

75

1. Se toma una parcela arbitraria de elementos, con forma no regular y con, al menos, un punto interior i. 2. Se toma una funcin polinmica que represente la variable de campo incgnita bsica o o o a u en la parcela. El orden del polinomio depende del tipo de problema, y debe ser tal que proporcione un campo de tensiones generalizadas constante dentro de la parcela. En Elasticidad, se toma una funcin lineal. o 3. Se determina el valor de la funcin anterior en todos los nodos de la parcela. o 4. Utilizando el programa de elementos nitos en el que se encuentre implementado el elemento que se prueba, se imponen en los nodos del contorno de la parcela los valores la variable incgnita bsica que les correspondan segn la funcin que se haya o a u o tomado (esto es, se imponen exclusivamente condiciones esenciales de contorno). 5. Se considera que el elemento cumple la prueba de la parcela si el programa de elementos nitos calcula el valor correcto de la funcin u en todos los nodos interiores o de la parcela y la tensin generalizada (constante) correcta en todos los elementos o de la parcela. La forma ms elaborada de la prueba es similar a la anterior, pero en el contorno de la a parcela se imponen no slo condiciones de contorno esenciales, sino tambin las condiciones o e de contorno naturales (tensiones generalizadas) compatibles con la funcin de campo u o seleccionada.

76

Cap tulo 8

Tecnolog de elementos (IV) a

8.1 Introduccin o

Para cerrar los temas dedicados a la tecnolog de elementos, se trata en este cap a tulo de acercarse a las tcnicas modernas de desarrollo de elementos nitos. e Ya se ha comentado que los elementos elaborados siguiendo las directrices de la formulacin convencional se emplean pocas veces en la prctica industrial diaria. Y ello o a a pesar de que las soluciones que proporcionan estos elementos tienen garantizada su convergencia, al hacer que la sistemtica del MEF se convierta en un procedimiento de a Rayleigh-Ritz-Galerkin. Los elementos que se utilizan en la prctica incorporan casi siempre algn \truco" a u que los aparta de la formulacin convencional. En el cap o tulo anterior se han introducido algunos de ellos: Integracin reducida con control de hourglassing. o Integracin reducida selectiva. o Proyeccin de deformaciones o deformaciones supuestas (B). o Modos incompatibles. En el cap tulo anterior se han mostrado algunas razones por las que estos \trucos" empezaron a utilizarse. La razn ultima es que en la prctica industrial hacen falta o a elementos nitos de altas prestaciones: elementos simples, en general con funciones de forma de orden bajo, que con mallados relativamente gruesos sean capaces de proporcionar soluciones con una precisin aceptable desde el punto de vista ingenieril. Y el hecho es que o la formulacin convencional del MEF no proporciona esas prestaciones en muchos casos o de alto inters industrial. e Los \trucos" aparecen as para dar respuesta a esa necesidad industrial imperiosa de llevar a cabo cada vez clculos ms complejos en el menor tiempo posible. Pero desa a de el punto de vista matemtico plantean numerosos interrogantes. >Por qu funcionan a e (convergen) los elementos con truco? >Por qu parece que funcionen incluso mejor que e 77

los elementos convencionales? >Convergen de verdad en todos los casos? >Podran fallar estos elementos en algn caso particular? u El fundamento matemtico de los elementos que se apartan de la formulacin convena o cional empez a vislumbrarse hacia nales de los 70 y principios de los 80. Se descrubri o o entonces que alguno de los elementos de integracin reducida utilizado para problemas o de anlisis de tensiones pod obtenerse sin ningn articio si se part de un principio a a u a variacional diferente del principio de la m nima energ potencial que utiliza la formulacin a o convencional. En problemas de elasticidad la formulacin clsica del MEF se basa en el principio de o a la m nima energ potencial. Habiendo muchos otros principios variacionales para plantear a el mismo problema, >por qu ha de esperarse que los elementos derivados del principio de e la m nima energ potencial sean los mejores? a Los elementos de altas prestaciones pueden derivarse, en general, de principios variacionales con varios campos primarios o independientes. Desde el punto de vista matemtico, a la tecnolog moderna de desarrollo de elementos est basada en esta idea: el empleo de a a principios variacionales multicampo distintos del principio de la m nima energ potencial. a La convergencia de las soluciones que proporcionan los elementos nitos formulados de este modo es mucho ms dif de estudiar que la de los elementos convencionales. a cil Tanto ms dif cuantos ms campos independientes haya en el principio variacional de a cil a partida. En general, la solucin del problema ya no hace m o nimo o mximo un potencial, a sino simplemente estacionario. El estudio de las propiedades de convergencia de esta clase de elementos es un campo de investigacin activo dentro de lo que es hoy en d la faceta matemtica del MEF. o a a El objetivo de este cap tulo es presentar una muestra de la tecnolog moderna de a desarrollo y/o justicacin de elementos de altas prestaciones. Para ello se va a justicar, o a partir de un principio variacional multicampo, el elemento cuadriltero no conforme de a Wilson y Taylor. Recurdese que este elemento fue introducido en la Seccin 7.2.1 de una e o manera heur stica. Luego se mencionarn las condiciones de Babuka-Brezzi, cuyo cumplimiento permite a s asegurar la convergencia de las soluciones obtenidas con cierta clase de elementos. Las condiciones de Babuka-Brezzi son un resultado clsico de la investigacin sobre la convers a o gencia de elementos nitos con formulacin mixta, esto es, con dos campos independientes. o

8.2

Principio variacional multicampo

Uno de los principios variacionales que puede utilizarse para plantear el problema elstico a es el principio de Hu-Washizu. Segn este principio, los campos solucin del problema u o elstico, u (desplazamientos), e (deformaciones), s (tensiones) y t (presiones en el contora no), son los que hacen estacionario el funcional1 : 1 ~ ~ ~ e e u ~ s t) s t; w (~ ; e; ~; ~ = (se ; ~) + (~; eu ~) (b; u) [t; u]St + [~ d u]Sd 2
1

(8.1)

Se emplea la notacin introducida en la seccin 3.4 para el problema elstico. o o a

78

cuando ste se dene sobre campos u, e, ~ y ~ sin ninguna restriccin, ms que las necee ~ ~ s t o a sarias para que las integrales tengan sentido. Puede pensarse en imponer a los campos u y e de 8.1 las restricciones: ~ ~

~ u = d en ~ e = eu + " ~ Entonces se obtiene un nuevo funcional:

Sd en

(8.2) (8.3)

u ~s s (~ ; "; ~) = =

1 e ~ ~ ~ s ~ (s ; e) + (~; eu e) (b; u) [t; u]St 2 1 u u 1 ~ ~ ~ ~ s ~ (s ; e ) + (su ; ") + (s" ; ") (~; ") (b; u) [t; u]St 2 2

y la solucin del problema elstico har estacionario este funcional entre los campos que o a a cumplan con las restricciones 8.2 y 8.3. La primera variacin del funcional s da: o s = (su + s" ~; ") (div(su + s" ) + b; ~ ) (~; ~) + [(su + s" )n t; ~ ]St s ~ " s u u La condicin de estacionaridad, s = 0, proporciona las ecuaciones de Euler correso pondientes a este principio variacional: div(su + s" ) + b = (su + s" )n t = su + s" = " = ~

0 0 ~ s 0

en en en en

(equilibrio en ) St (equilibrio en St)

8.3

Discretizacin por el MEF o

La discretizacin por el MEF del principio variacional anterior da lugar a elementos de o deformaciones supuestas2 . La sistemtica del MEF es ahora la siguiente: a 1. El dominio de clculo se divide en parcelas o elementos: = a
S
i i .

2. Dentro de cada elemento i , se aproximan los campos que intervienen en el funcional s : ~ u = N ae

2

(8.4)

En la literatura anglosajona estos elementos se conocen como elementos ANS (\Assumed Natural Strain").

79

eu =

1 u (r + rt)~ = B ae 2 u s = D B ae " = G ce ~ ~ = S pe s

(8.5) (8.6) (8.7) (8.8)

donde N son las funciones de forma serend pitas habituales, G es la llamada matriz de modos de deformacin mejorada y S es la matriz de modos de tensin. Para no o o 3 y evitar la aparicin de matrices tropezar con los llamados principios de limitacin o o singulares, en la eleccin de las funciones G y S se deben respetar las condiciones o siguientes: Las columnas de G deben ser linealmente independientes. Las columnas de G deben ser tambin linealmente independientes de las coe lumnas de B. Las columnas de S deben incluir los modos de tensin constante en el elemento. o Las matrices G y S han de ser ortogonales en el sentido de que:
Z

St G d = 0

3. Las aproximaciones locales 8.4 a 8.8 dan lugar una aproximacin global al campo de o ~ desplazamientos u que es continua en el dominio de clculo ; y a aproximaciones a globales a los campos " y ~ continuas a trozos en . As el funcional s puede ~ s , ponerse como suma de funcionales similares a l pero denidos en cada elemento: e s (~ ; "; ~) = u ~s
X
i

si (~ ; "; ~) u ~s

4. El funcional s se har estacionario cuando todos los funcionales si se hagan esa tacionarios. De este modo, basta con estudiar las condiciones de estacionaridad de uno de estos funcionales extendido a un elemento. Para simplicar la notacin, a o partir de este punto se suprime el sub ndice i indicativo del dominio elemental. 5. La sustitucin de las expresiones 8.4 a 8.8 en el funcional elemental proporciona su o expresin discreta: o 1 1 t s = at (Bt DB) ae +at (Bt DG) ce + ct (Gt DG) ce pt (St G) ce (bt N) ae [t N]St ae e e 2 e 2 e Ntese que el cuarto sumando es nulo por la condicin de ortogonalidad impuesta a o o las matrices G y S. La matriz D es la matriz de coecientes elsticos que permite a obtener tensiones a partir de las deformaciones.
Los principios de limitacin hacen referencia a que si las funciones de aproximacin a los campos o o independientes de tensiones y de deformaciones son del mismo orden, entonces con estos elementos de formulacin mixta se obtienen los mismos resultados que con los elementos de formulacin convencional y, o o por tanto, no merece la pena utilizarlos. Los principios de limitacin detuvieron el desarrollo de elementos o nitos con varios campos independientes durante las primeras etapas de la evolucin del MEF. o
3

80

6. Las condiciones de estacionaridad del funcional son: @s @ae @s @ce = (Bt DB) ae + (Bt DG) ce (bt N) [t N]St = 0 = (Gt DB) ae + (GtDG) ce = 0
t

Si se denen las matrices: Ku (Bt DB) Ke (Gt DG) Q (Bt DG) fb (bt N) t ft [t N]St las dos ecuaciones anteriores pueden escribirse:
" #( ) ( )

Ku Q Qt Ke

ae ce

fb + ft 0

(8.9)

7. Por las condiciones para la eleccin de G, la matriz Ke es regular. Entonces el o sistema de ecuaciones 8.9 se puede condensar a: [Ku QK1 Qt ]ae = fb + ft e donde la matriz en el primer miembro es la matriz de rigidez elemental en el sentido habitual y el segundo miembro corresponde al vector de cargas. Obsrvese que la e matriz de rigidez elemental se construye sumando a la matriz correspondiente a la formulacin convencional en desplazamientos un trmino corrector que depende de o e los modos de deformaciones supuestas. El vector de cargas es el mismo que en la formulacin convencional. o 8. Conocidos la matriz de rigidez elemental y el vector elemental de cargas, se procede al ensamblaje del sistema global de ecuaciones como en la formulacin convencional. o El desarrollo anterior puede utilizarse para justicar el elemento de Wilson y Taylor en elasticidad bidimensional, ya que la matriz de rigidez correspondiente a este elemento se obtiene si se elige la matriz G como:
8 9 > "x > < = > y > : " ; xy

"= ~

"

= G ce

0 0 0 6 7 G=4 0 0 0 5 0 0

donde y son las coordenadas en el dominio cuadrado estndar. a

81

8.4

Las condiciones de Babuka-Brezzi s

La formulacin clsica del MEF, basada en la aproximacin de un solo campo indepeno a o diente, exige una serie de condiciones a las funciones de forma o funciones de aproximacin o 4 local de ese campo independiente. Si se cumplen esas condiciones , entonces la convergencia de las soluciones aproximadas hacia la solucin exacta est garantizada. o a Cuando se abandona la formulacin convencional para trabajar con elementos que o responden a principios variacionales multicampo, resulta mucho ms dif determinar las a cil ~ condiciones que deben cumplir las aproximaciones locales a los campos primarios (p.ej. u, " ~) si se quiere tener garantizada la convergencia. Se deben analizar ahora, adems de ~o s a la aproximacin en s a cada campo, las relaciones entre las aproximaciones a uno y a otro o campo. El estudio de las propiedades de convergencia de los elementos formulados a partir de principios variacionales multicampo es un area abierta de investigacin que utiliza un o lenguaje dif para los no iniciados en anlisis funcional. cil a En este campo hay un resultado ya clsico que se desarroll en el estudio del problema a o de la elasticidad incompresible con dos campos independientes: el campo de desplazamien~ tos u y el campo de presiones p (primer invariante del tensor de tensiones). Ese resultado ~ son las condiciones de Babuka-Brezzi. Condiciones que, si son cumplidas, garantizan la s existencia y unicidad de la solucin aproximada y una tasa optima de convergencia. Se o presentan a continuacin como muestra de otros resultados similares o de otros resultados o que ser deseable tener disponibles. a En el estudio de problemas de elasticidad incompresible, partiendo de un principio variacional con los dos campos independientes citados y siguiendo un proceso anlogo al a descrito en la seccin anterior al estudiar el elemento de Wilson y Taylor, se llega una o ecuacin similar a la 8.9: o
"

Ku Q Qt 0

#(

ae pe

f h

(8.10)

donde pe son ahora las n variables nodales asociadas al campo de presiones5 . Un elemento nito cumple las condiciones de Babuka-Brezzi cuando se verica que: s = 1. Para cualquier conjunto de desplazamientos nodales ae 60 que cumpla Qt ae = 0, se tiene que: at Ku ae > 0 e 2. La matriz Qt tiene n las linealmente independientes. Si se quiere tener garantizada la convergencia, las funciones de aproximacin local o del campo de desplazamientos y del campo de presiones deben elegirse de manera que se cumplan estas dos condiciones.
4 Recordemos que esas condiciones eran tres: suavidad en el interior del elemento, continuidad a travs de e las fronteras entre elementos y completitud hasta un orden determinado que depende del tipo de problema. 5 Los nodos asociados al campo de presiones pueden ser diferentes de los asociados al campo de desplazamientos.

82

Cap tulo 9

Procedimientos de clculo (I) a

9.1 Introduccin o

Dentro del esquema general de la sistemtica del MEF para clculos lineales, una vez a a construida la matriz de rigidez global y el vector global de cargas por ensamblaje de las contribuciones elementales, el paso siguiente es resolver un sistema de ecuaciones lineales. Ello permite obtener las variables nodales. Los procedimientos numricos para la resolucin de este sistema son, en principio, los e o mismos que se utilizan en otras ramas de la tcnica. Lo que ocurre es que, al ser el tiempo e de ordenador dedicado a esta fase del clculo una fraccin muy importante del tiempo a o total utilizado para resolver el problema, los investigadores han ido creando una serie de tcnicas especiales adaptadas a las caracter e sticas de los sistemas de ecuaciones a que da lugar el MEF. As se ha desarrollado una rama de especializacin dentro del estudio de los mtodos , o e de elementos nitos: la que trata de la resolucin del sistema de ecuaciones y, en problemas o con variacin de la solucin en el tiempo, del acoplamiento con la integracin en el tiempo. o o o En este cap tulo y en el siguiente se trata de proporcionar una panormica de los a procedimientos de solucin que ms se emplean en la prctica. En este cap o a a tulo se tratar a de los mtodos de resolucin de sistemas de ecuaciones lineales, que tienen su aplicacin e o o en problemas estacionarios o sin variacin de la solucin en el tiempo. En el cap o o tulo siguiente se estudiar cmo se acoplan estos mtodos con la integracin en el tiempo, para a o e o resolver problemas no estacionarios o transitorios.

9.2

Resolucin de sistemas de ecuaciones lineales o

Dentro de la prctica del MEF se emplean dos grandes familias de procedimientos para a resolver los sistemas de ecuaciones lineales a que da lugar el mtodo. Esquemticamente, e a estos procedimientos son: 1. Mtodos de solucin directa, por ejemplo: e o Eliminacin de Gauss o 83

 Factorizaciones (Cholesky, Crout) Mtodo frontal e 2. Mtodos iterativos, por ejemplo: e Mtodo de Jacobi e

Mtodo del gradiente conjugado e Relajacin dinmica o a

En los mtodos de solucin directa, dado el sistema de ecuaciones que debe resolverse, e o se conoce a priori el nmero de operaciones necesarias para obtener la solucin. Son los u o procedimientos con ms tradicin dentro de la tecnolog del MEF y su rango de aplicacin a o a o es general, tanto en programas de clculo lineal como en los de clculo no lineal. Todos los a a procedimientos que se utilizan dentro de esta categor son elaboraciones de la eliminacin a o de Gauss. En los mtodos iterativos, por el contrario, no se conoce a priori el nmero de opee u raciones necesarias para llegar a la solucin. Son procedimientos que estn adquiriendo o a actualmente mucha difusin porque dan lugar a menos necesidades de almacenamiento o en el ordenador (memoria, disco) y, por tanto, permiten abordar problemas ms grandes a manteniendo los mismos recursos. La dicultad con que tropiezan los mtodos iterativos es que la velocidad con la que e llegan a la solucin depende de cada problema concreto, en particular, de las caractersticas o de la matriz de coecientes del sistema. Entonces, resulta que a veces son mucho ms a ecientes que los mtodos de solucin directa pero por contra, en otros casos, pueden llegar e o requerir un nmero bastante mayor de operaciones. Esto es especialmente perturbador u cuando se abordan clculos no lineales, en los que las matrices de coecientes de los a diferentes sistemas de ecuaciones que se resuelven en un mismo problema pueden tener caracter sticas muy diferentes de una fase a otra del clculo. a De este modo, parece que los mtodos iterativos no son procedimientos de propsito e o tan general como los de solucin directa. Es normal que los programas comerciales de o clculo los incorporen como opcin. a o Sea el sistema de n ecuaciones lineales: Ka = f (9.1)

donde K es una matriz de coecientes de n n, a es el vector de incgnitas y f es el vector o de trminos independientes. Dentro de las aplicaciones del MEF, K representa la matriz e de rigidez global, f es el vector global de cargas y a es el vector de movimientos nodales. A primera vista parece que la resolucin del sistema 9.1 es un asunto balad Se trata o . simplemente de invertir la matriz de coecientes, para dar: a = K1 f Sin embargo, cuando el orden de K es de varias decenas de miles, se comprende que hay que prestar mucha atencin a esta parte del clculo para obtener la solucin del sistema en o a o 84

un tiempo de ordenador razonable y con unos recursos de memoria y disco no demasiado grandes. En particular, la inversin de la matriz de coecientes ser una forma muy o a ineciente de resolver este sistema, ya que no se sacara ningn partido de la estructura u interna de la matriz K derivada de la sistemtica del MEF. La sistemtica del MEF da a a lugar a una matriz de rigidez global, en general simtrica, en la que la mayor de los e a coecientes son nulos (matriz dispersa), concentrndose los coecientes no nulos alrededor a de la diagonal principal (estructura en banda). Este es el tipo de matriz de coecientes en el que estamos interesados.

9.3

Eliminacin de Gauss o

El algoritmo de la eliminacin de Gauss para resolver el sistema 9.1 puede esquematizarse o de la forma siguiente: 1. Eliminacin y sustitucin hacia adelante1 : o o De i = 1 a i = n 1 (recorre columnas) (recorre las) (dentro de cada la recorre parte de las columnas)

De l = i + 1 a l = n De j = i a j = n klj klj Otro j fl fl

kli kii fi kli kii kij

Otro l Otro i 2. Sustitucin hacia detrs2 : o a De i = n a i = 1 fi Otro i El algoritmo anterior se visualiza fcilmente por medio del ejemplo siguiente. Sea el a sistema de ecuaciones:
2 6 6 6 4
1 2

(recorre las)
kij fj kii

fi

Pn

j=i+1; i6 =n

4 1 0 2 1 6 2 0 0 2 6 1 2 0 1 4

9 8 9 38 > a1 > > 1 > > > > > > = < 7> a > > < 2 = 7 2 = 7 5 > a3 > > 3 > > > > > > > > > ; : ; :

a4

En esta fase se reduce la matriz de coecientes a una matriz triangular superior En esta fase se remplaza el vector de trminos independientes por el vector solucin del sistema. e o

85

La primera fase del algoritmo, la de eliminacin y sustitucin hacia adelante, transforo o ma progresivamente el sistema del siguiente modo: i=1 9 8 9 2 38 4 1 0 2 > a1 > > 1 > > > > > < = < = 6 0 5:75 2 0:5 7 > a > > 2:25 > 6 7 2 = 6 7 4 0 2 6 1 5 > a3 > > 3 > > > > > > > > > 0 0:5 1 3 : a4 ; : 3:5 ; i=2
2 6 6 6 4

4 1 0 2 0 5:75 2 0:5 0 0 5:304 0:826 0 0 0:826 2:957 4 1 0 2 0 5:75 2 0:5 0 0 5:304 0:826 0 0 0 2:828

i=3

2 6 6 6 4

Ntese que el algoritmo mantiene la simetr de las ecuaciones activas, con lo que, o a desde el punto de vista de la programacin, no se incrementan durante el proceso las o necesidades de almacenamiento. Es decir, el algoritmo se puede programar unicamente a partir del espacio reservado para los datos, que son los coecientes no repetidos de la matriz K y el vector f . La segunda fase del algoritmo, la sustitucin hacia detrs, permite obtener las incgnitas o a o del modo siguiente: i=4 3:893 a4 = 2:828 i=3 3:783 + 0:826 a4 a3 = 5:304 i=2 2:25 0:5 a4 + 2:0 a3 a2 = 5:75 i=1 1:0 2 a4 + 1:0 a2 a1 = 4:0 El procedimiento de eliminacin de Gauss es numricamente muy robusto y forma la o e base de todos los mtodos de solucin directa que se utilizan en la actualidad en el ambito e o del MEF. Lo que cambia de un caso a otro es la estructura de datos (forma de almacenar K y f : en banda, en \sky-line", por bloques, . . . ) y la organizacin de las distintas o operaciones. As en la literatura aparecen multitud de procedimientos, aparentemente , distintos (p.ej. resolucin por bloques, resolucin en \sky-line", etc. . . . ), pero cuya o o esencia matemtica es la misma. a

9 8 9 38 > a1 > > 1 > > > > > > > > > 7< a = < 2:25 = 7 2 = 7 5 > a3 > > 3:783 > > > > > > > > > ; : ; :

9 8 9 38 > a1 > > 1 > > > > > > < = < 7> a > > 2:25 = 7 2 = 7 5 > a3 > > 3:783 > > > > > > > > > : ; : ;

a4

3:304

a4

3:893

86

9.4

Factorizaciones de Crout y Cholesky

Las factorizaciones de Crout y Cholesky son variantes de la eliminacin de Gauss. Su o objetivo es separar las operaciones sobre la matriz de coecientes del sistema de las operaciones sobre el vector de trminos independientes. Esto puede resultar interesante cuando, e por ejemplo, se quiere calcular una misma estructura (misma matriz de coecientes) para diferentes hiptesis de carga (distintos vectores de trminos independientes). Estas faco e torizaciones son una manera de tener almacenada la primera parte de la eliminacin de o Gauss, la ms costosa, para usarla repetidas veces. a

9.4.1

Factorizacin de Crout o

Se conoce tambin con el nombre de factorizacin de Cholesky modicada. Consiste en e o poner la matriz K como producto de tres matrices: K = L D Lt donde L es una matriz triangular inferior y D es una matriz diagonal. La matriz D contiene los trminos de la diagonal principal de la matriz de coecientes reducida tras e la primera fase de la eliminacin de Gauss o sustitucin hacia adelante. Por otro lado, la o o matriz Lt contiene las las de esa matriz de coecientes reducida tras dividirlas por los trminos de la diagonal contenidos en D. e En el ejemplo utilizado en la seccin anterior, la factorizacin de Crout ser o o a:
2 6 6 4

K=6

1 0 0 0:25 1 0 0 0:348 1 0:5 0:087 0:156

0 0 0 1

32 76 76 76 54

4 0 0 0 0 5:75 0 0 0 0 5:304 0 0 0 0 2:828

32 76 76 76 54

1 0:25 0 0:5 0 1 0:348 0:087 0 0 1 0:156 0 0 0 1

3 7 7 7 5

Es decir, las ecuaciones reducidas tras la primera etapa de la eliminacin de Gauss se o pueden escribir como: D Lt a = f e , donde f es el vector reducido de trminos independientes. As el algoritmo de la primera fase de la eliminacin de Gauss o sustitucin hacia adelante sirve para determinar las o o matrices D y L. Ntese que, una vez factorizada K, se tiene que: o K a = L (DLt ) a = L f = f de manera que el proceso de resolucin del sistema de ecuaciones puede descomponerse o en: 1. Obtener el vector reducido de trminos independientes f : e f = L1 f o La matriz L1 no se calcula. El vector f se obtiene por sustitucin hacia adelante. 87

2. Obtener el vector solucin a: o a = LtD1 f Las matrices Lt y D1 tampoco se calculan. El vector a se obtiene por sustitucin o hacia detrs. a De esta manera, cuando se tiene un mismo modelo sometido a varios conjuntos de acciones diferentes, una vez factorizada la matriz K, basta con repetir los pasos 1 y 2 anteriores para cada una de los distintos vectores f que se tengan. El nmero de opeu raciones necesario para ejecutar estos dos pasos, sustitucin hacia adelante y sustitucin o o hacia detrs, es muy peque~o si se compara con el nmero de operaciones necesarias para a n u factorizar la matriz.

9.4.2

Factorizacin de Cholesky o

La factorizacin de Cholesky se desarroll antes que la factorizacin de Crout, pero se o o o utiliza menos hoy en d porque el algoritmo de descomposicin de la matriz de coeciena, o tes utiliza ra ces cuadradas. Adems, no puede utilizarse en problemas con matrices de a coecientes que tengan autovalores negativos, ya que se obtienen entonces ra cuadraces das de nmeros negativos. Sin embargo, se ha empleado mucho esta descomposicin en u o programas de clculo matricial de estructuras, donde la matriz de rigidez global es denida a positiva. La factorizacin de Cholesky pone la matriz K como producto de dos matrices: o K = Ut U donde U es una matriz triangular superior que resulta ser: U = D1=2 Lt El algoritmo de obtencin de U puede esquematizarse como: o
p

u11 = u1j = uii

k11 j = 2; : : : ; n i = 2; : : : ; n

v u i1 u X u2 = tkii ki
k=1

k1j u11

uij =

i1 X 1 uki ukj ) (kij uii k=1

i = 2; : : : ; n

j = i + 1; : : : ; n

uij = 0

i>j

88

En el ejemplo utilizado en la seccin anterior, el resultado de la factorizacin de Choo o lesky ser a:

2 6 6 4

K=6

2 0 0 0 0:5 2:398 0 0 0 0:834 2:303 0 1:0 0:208 0:359 1:682

32 76 76 76 54

2 0:5 0 1:0 0 2:398 0:834 0:208 0 0 2:303 0:359 0 0 0 1:682

3 7 7 7 5

Al igual que en el caso de la factorizacin de Crout, el proceso de resolucin del sistema o o de ecuaciones puede descomponerse en dos pasos: 1. Obtener el vector reducido de trminos independientes f : e f = Ut f o La matriz Ut no se calcula. El vector f se obtiene por sustitucin hacia adelante. 2. Obtener el vector solucin a: o a = U1 f La matriz U1 tampoco se calculan. El vector a se obtiene por sustitucin hacia o detrs. a

9.5

Mtodo de resolucin frontal e o

Es un mtodo desarrollado por Irons para resolver las ecuaciones resultantes del MEF e (1970) y es el procedimiento ms utilizado actualmente por los programas comerciales a modernos. Fue pensado para optimizar los requerimientos de almacenamiento en memoria en unas circunstancias en las que la memoria del ordenador era mucho ms cara que en la a actualidad. En el fondo, el mtodo de resolucin frontal es tambin un procedimiento de eliminacin e o e o de Gauss en el que se hace la sustitucin hacia adelante al mismo tiempo que se ensambla o la matriz de rigidez global. Es decir, no se espera a tener ensamblado el sistema de ecuaciones para empezar a resolverlo, sino que se empieza a resolver el sistema mientras se est todav ensamblando. a a Se gana en eciencia con respecto a otras implementaciones de la sustitucin de Gauss o pero, como contrapartida, se pierde modularidad en la programacin y se trabaja con una o estructura de datos bastante complicada. La mayor ventaja del mtodo es que se minimiza e el uso de la memoria. La mecnica detallada y general del mtodo es muy elaborada. El lector interesado a e puede consultar la literatura3 . En los prrafos que siguen se intenta, mediante un ejemplo, a dar una idea intuitiva de cmo funciona el procedimiento. o Sea la malla de cuatro elementos representada en la gura y, para simplicar, supongamos que slo se tiene un grado de libertad por nodo. o
3

E. Hinton y D.R. Owen. Finite Element Programming. Academic Press. 1977.

89

El proceso del mtodo de resolucin frontal puede esquematizarse como sigue: e o 7 x 4 8

x x9

4 x

x6

2 frente

1. Se toma el elemento 1 y se hace una particin de los grados de libertad asociados al o mismo, a, en: ( ) ac a= ab donde ac rene los grados de libertad no compartidos con otros elementos (nodo 1) y u ab , los grados de libertad compartidos (nodos 2,4 y 5). Esto da lugar a una particin o de la matriz de rigidez elemental y del vector elemental de cargas del elemento 1:
"

Kcc Kcb Kbc Kbb

#(

ac ab

fc fb

(9.2)

2. Las ecuaciones globales que involucran a los grados de libertad ac ya estn completas, a puesto que el resto de los elementos no les a~aden ms trminos. De este modo, los n a e grados de libertad internos ac pueden ser eliminados a este nivel: ac = K1 fc K1 Kcb ab = K1 [fc Kcb ab ] cc cc cc lo que da lugar, sustituyendo en 9.2, al siguiente sistema para los grados de libertad externos: (Kbb Kbc K1 Kcb )ab = fb Kbc K1 fc cc cc o, escrito de otra manera: Kbb ab = f b 90

3. Las matrices y vectores necesarias para obtener ac a partir de ab (esto es, K1 , fc cc y Kbb ) se almacenan en disco para ser utilizados posteriormente en el proceso de sustitucin hacia detrs. o a 7 x 4 8
x x9

4 x

x6

2 frente

4. Se calcula la matriz de rigidez K0 y el vector elemental de cargas f 0 del siguiente elemento, el elemento 2. 5. La matriz K0 y el vector f 0 se ensamblan en Kbb y f b para obtener una matriz de rigidez y un vector de cargas conjunto K y f . 6. Se determinan los nuevos grados de libertad internos ac y externos ab del sistema conjunto y se repite el proceso de los puntos 1 a 3 anteriores. 7. Se opera de la misma manera con el elemento 3. 8. Al repetir el proceso para el ultimo elemento, el elemento 4, resulta que se puede resolver el sistema de ecuaciones 9.2 sin pasar adelante ninguna incgnita, esto es, no o o hay grados de libertad externos ab . Entonces se inicia la vuelta o sustitucin hacia detrs, recuperando de disco la informacin correspondiente a los distintos elementos a o en orden inverso, esto es, 3-2-1 y obteniendo en cada caso los grados de libertad que se hab ido eliminando, los grados de libertad ac . an

91

7 x 4

8
x

x9

4 x

x6

2 frente

En cada paso o incorporacin de un nuevo elemento, se llama frente activo o frente o de onda a los grados de libertad externos ab que resultan de la particin de la matriz del o sistema en ese momento. La cantidad de memoria necesaria para desarrollar el proceso depende de la longitud mxima o nmero mximo de grados de libertad en el frente activo. Esta longitud mxima, a u a a al contrario de lo que sucede en otros procedimientos, no depende de cmo se numeren los o nodos, sino del orden o nmero asignado a los elementos. u En las aplicaciones, resulta util saber que el tiempo de ordenador que emplea el mtodo e frontal es proporcional al nmero de ecuaciones (nmero de grados de libertad) y al cuau u drado de la media cuadrtica de la longitud del frente activo durante el proceso. Es muy a importante entonces trabajar con numeraciones de los elementos que minimicen dicha media cuadrtica. Los programas comerciales de clculo, con objeto de dejar libertad al a a usuario para numerar los elementos, incorporan normalmente un optimizador que renumera internamente los elementos a los efectos de la resolucin del sistema mediante este o procedimiento.

9.6

Mtodos iterativos e

El sistema de ecuaciones que se trata de resolver puede escribirse: r = Ka f = 0 (9.3)

y puede interpretarse desde el punto de vista mecnico como la bsqueda del equilibrio a u 92

entre las fuerzas exteriores, f, aplicadas sobre los nodos y las fuerzas interiores, K a, generadas por los movimientos nodales a. En este sentido, el vector r o vector de residuos representa el desequilibrio entre fuerzas interiores y fuerzas exteriores para un conjunto de movimientos nodales diferente de la solucin del sistema de ecuaciones. o Un procedimiento iterativo general para resolver el sistema 9.3 puede ponerse como: ai+1 = ai + i i i = i K1 ri + i i1 a

(9.4) (9.5)

donde el ndice i corresponde al nmero de iteracin. Los distintos procedimientos iterau o tivos se obtienen particularizando las expresiones 9.4 y 9.5. El vector i tiene el sentido de un vector de mejora de la solucin en la iteracin i. El o o i escalar ser un optimizador de la correccin a lo largo de la direccin dada por i . a o o o La matriz Ka es una aproximacin a la matriz de coecientes del sistema K. Se conoce con el nombre de matriz precondicionadora y puede ser cualquier matriz, normalmente denida positiva, comprendida entre la matriz identidad I y la matriz de coecientes del sistema K. Para que el procedimiento sea en realidad un procedimiento iterativo, la matriz Ka ha de ser ms fcilmente invertible (factorizable) que K. a a El escalar i se introduce a veces para intentar acelerar la convergencia y el escalar i sirve para introducir en el vector de correccin de la solucin i una correccin relacionada o o o con la correccin en la iteracin anterior. o o El procedimiento iterativo recogido en las relaciones 9.4 y 9.5 se detiene en cuanto se alcanza el criterio de convergencia elegido. Existen dos categor de criterios de esta as clase: 1. Criterios en desplazamientos. Son criterios de la forma: k k =
q
i

q
t

i i < "1 kai+1 k

q

i Ka i < "2

ai+1t Ka ai+1

2. Criterios en fuerzas, de la forma: kri k =

p

rit ri < "3 max(kRk; kf k)

donde R es el vector de reacciones. Los escalares "1 , "2 y "3 se conocen con el nombre de tolerancias. Lo que hace atractivo un procedimiento iterativo como el descrito en los prrafos a anteriores es que para calcular el vector de residuos ri no es necesario llegar a ensamblar la matriz de coecientes K. El vector de residuos puede obternerse ensamblando las contribuciones de los elementos al equilibrio de cada nodo: 93

ri =

^
e

(Ke ai fe ) e

Esto hace que un procedimiento iterativo requiera menos recursos de memoria y disco que un procedimiento de solucin directa cuando se programa en el ordenador. o

9.7

Iteracin de Jacobi y sus variantes o

El procedimiento iterativo ms simple es la iteracin de Jacobi. Se obtiene de 9.4 y 9.5 a o poniendo: i = 1 i = 1 i = 0 Ka = I El procedimiento queda entonces: ai+1 = ai ri = ai (K ai f ) = f + (I K) ai Este procedimiento tan simple es divergente siempre que los autovalores de K sean mayores o iguales que 2. En efecto, si a es la solucin exacta, la aproximacin ai se puede o o poner como: ai = a + "i donde "i es el error en la iteracin i. o entonces se tiene que: ai+1 = f + (I K) (a + "i ) = f + a + "i K "i K a = a + (I K) "i esto es, "i+1 = (I K) "i (9.6)

En los sistemas de ecuaciones que se generan en el MEF aplicado a problemas lineales, la matriz de coecientes K es simtrica y, generalmente, denida positiva. En este caso, e la matriz K, de orden n, tiene n autovalores reales y positivos, j ; j = 1; : : : ; n. Sea zj ; j = 1; : : : ; n una base de los autovectores de la matriz K. El vector error "i se podr poner como combinacin lineal de los vectores de esta base: a o "i =
n X

cj zj

(9.7)

j=1

94

con:

K zj = j zj zt zj = 0 k k 6j =
n X

De este modo, si se sustituye 9.7 en 9.6, se llega a: "i+1 =

n X j=1

cj zj

j=1

n X

cj j zj =

j=1

cj (1 j ) zj

Entonces, si se quiere que la componente j del vector "i+1 sea menor que la correspondiente al vector "i , debe ser: j1 j j < 1 es decir, 0 < j < 2 j = 1; : : : ; n

En otro caso, la iteracin de Jacobi diverge. o Ntese que la convergencia ser tanto ms rpida cuanto ms agrupados estn los o a a a a e autovalores j en torno a 1. Puede conseguirse agrupar los autovalores en torno a 1 con el escalado o precondicionamiento de la matriz de coecientes K. Los autovalores mximo y m a nimo de K estn acotados por: a
X
j6 =i

max max(Kii +
i

jKij j)

min min(Kii
i

X
j6 =i

jKij j)

Aprovechando esta propiedad, los autovalores pueden ser agrupados en torno a 1 si se escalan las ecuaciones de modo que los trminos de la diagonal principal de la matriz K e sean todos iguales a 1. Esto es equivalente a utilizar como matriz precondicionadora, Ka , una matriz diagonal con la diagonal principal de K. Adems, para mejorar la convergencia se puede introducir el acelerador en 9.5, con a un valor: = 2 max + min

Tambin se puede mejorar la convergencia encontrando en cada iteracin la longitud e o i i o ptima de la correccin segn , esto es, introduciendo el escalar en 9.4. El parmetro o u a 95

se calcula en cada iteracin mediante una bsqueda de la correccin optima en la direccin o u o o denida por i . En la literatura anglosajona esta bsqueda se conoce como \line search". u Para ver cmo podr calcularse , resulta util imaginarse que el sistema de ecuaciones o a K a = f deriva de hacer estacionario el potencial: = ya que: @ = Ka f = r @a Ntese que el sentido que tiene el residuo es el de gradiente de . Resolver el sistema de o ecuaciones es pues encontrar el vector a para el cual el potencial se hace estacionario (residuo nulo). Entonces, si dentro del proceso iterativo se utiliza 9.4 para obtener la solucin mejorada o i como i son jos, resulta que se puede determinar i con la y se considera que tanto a condicin de que el potencial se haga estacionario con respecto a i . Esto es, sustituyendo o 9.4 en 9.8: 1 i (a + i i )t K (ai + i i ) (ai + i i )t f 2 1 it t t t 1 2 t a K ai + i ai K i + i i K i ai f i i f 2 2 1 t a K a at f 2 (9.8)

= =

@ t t t = ai K i + i i K i i f @ de donde se obtiene, igualando a cero el gradiente de : i (f K ai ) i ri = = it t i K i K i

i
t t

Ntese que el denominador del cociente anterior puede obtenerse teniendo en cuenta que: o K i = ri+1 ri para i = 1.

9.8

Mtodo del gradiente conjugado e

El mtodo del gradiente conjugado encaja en el esquema iterativo denido por 9.4 y 9.5 e tomando: i = 1 Ka = I el escalar i se obtiene mediante \line search" y el coeciente i se escoge de modo que:

96

 i K j = 0

si i 6j =

(9.9)

Es decir, cada nuevo vector de correccin ha de ser ortogonal en K a todos los anteriores. o En el mtodo del gradiente conjugado precondicionado se utiliza una matriz Ka diferente e de la matriz identidad. Ntese que, de acuerdo con el esquema iterativo, la solucin en la iteracin i + 1 puede o o o ponerse como: ai+1 = a0 +
i X j=0

j j

Entonces, si i + 1 = n, el procedimiento anterior conduce a una correccin nula para o la solucin aproximada, ya que al llegar a este punto se habrn tomado n vectores de o a correccin j que, por ser K una matriz regular, dan lugar a una base del espacio de o dimensin n en el que se busca la solucin a. Al tener que ser la correccin siguiente o o o ortogonal a todos los vectores de esta base, la correccin no puede ser ms que el vector o a nulo. De este modo, la idea del mtodo del gradiente conjugado es garantizar que se e obtiene la solucin del sistema en n ciclos como mximo. o a La determinacin de los coecientes i se hace a partir de los resultados siguientes: o 1. La forma en que se determinan los coecientes i mediante \line search" implica la ortogonalidad del residuo a la direccin de correccin, esto es: o o i ri+1 = 0 En efecto, los coecientes i se determinan haciendo estacionario el potencial , luego: @ t t t t t = ai K i + i i K i i f = i ri + i i K i = 0 @ y se cumple adems que: a ri+1 ri = K(ai + i i ) f K ai + f = i K i Por tanto, la forma en que se determinan los coecientes i implica que: 0 = i ri + i i
t t t 1 i+1 ri ) = i ri+1 (r i t

2. Cuando se cumple la condicin 9.9, entonces los sucesivos vectores de residuos (grao dientes) son ortogonales entre s esto es4 : , ri+1 rk = 0
4 t

k = 0; 1; : : : ; i

De aqu el nombre de gradiente conjugado.

97

Esta propiedad es la que garantiza la convergencia en, como mximo, n iteraciones. a En efecto, se tiene que: ai+1 = ak+1 + y entonces, ri+1 k = (K(ak+1 +
t

j=k+1

i X

j j

k = 0; 1; : : : ; i 1
i X

j=k+1

j j ) f )t k

y por la condicin de ortogonalidad entre los distintos vectores de correccin: o o ri+1 k = rk+1 k = 0 Finalmente, por ser k = rk + k k1 : 0 = ri+1 k = ri+1 rk + k ri+1 k1 = ri+1 rk
t t t t t t

k = 0; 1; : : : ; i 1

3. Los coecientes se determinan imponiendo el cumplimiento de la condicin 9.9, es o decir, se tiene que: i = ri + i i1 luego i K i1 = ri K i1 + i i1 K i1 y se busca que el primer miembro de la ecuacin anterior sea cero, entonces debe o ser: t ri K i1 i = i1t K i1 La expresin de i se puede simplicar teniendo en cuenta que, segn el ep o u grafe 1 anterior, se tiene: ri ri1 = i1 K i1 t i1 ri = 0 Entonces, por la ortogonalidad de los vectores de residuos entre s y de los vectores de residuos con los vectores de correccin, se llega a: o i = ri ri i1t ri1 r
t t t t

Para llegar a la expresin nal anterior se ha hecho uso de que: o i1 = ri1 + i1 i2

t t t

De esta forma, en el clculo de los coecientes slo intervienen los vectores de a o residuos correspondientes a la iteracin actual y a la iteracin anterior. o o 98

9.9

Relajacin dinmica o a

La solucin del sistema de ecuaciones 9.1 puede considerarse que es el lmite al que tiende o la respuesta del sistema dinmico amortiguado cticio: a Ma + Ca + Ka = f _ (9.10)

partiendo de unas condiciones iniciales arbitrarias. De este modo, para resolver el sistema 9.1 pueden emplearse las mismas tcnicas e numricas de integracin que se utilizan para obtener la respuesta en problemas no ese o tacionarios. Dichas tcnicas se interpretan entonces como procedimientos iterativos de e resolver 9.1. Desde este punto de vista, cada paso de integracin o \salto" en el tiempo es o un ciclo de mejora de la solucin. o Las matrices de masas, M, y de amortiguamiento, C, pueden elegirse de modo que la integracin de 9.10 resulte lo ms cmoda posible. En particular, si se elige una matriz de o a o masas diagonal, es posible utilizar tcnicas de integracin expl e o cita, las cuales no requieren el ensamblaje de la matriz K de coecientes del sistema. Un ejemplo de ciclo de clculo a ser el siguiente: a 1. Inicializacin o i = 0 a a _
0 1 2

= 0 = 0

2. Bucle en elementos: clculo de fuerzas internas pi a Contribucin del elemento e: o pi = Ke ai Ce ae _ e e Ensamblaje de pi en pi e 3. Bucle en nodos: clculo de aceleraciones e integracin en el tiempo (mejora de la a o solucin) o Clculo de aceleraciones (residuos): a ai = M1 (f + pi ) Integracin en el tiempo (mejora de la solucin) por diferencias centrales o o _ ai+ 2 ai+1 4. i i + 1 99
1

i 1 2

_ = ai 2 + t ai i i+ 1 _ = a + t a 2

5. Volver a 2 Como se ver en el cap a tulo siguiente, este esquema de integracin expl o cita se emplea mucho para obtener la respuesta en problemas no estacionarios. El ciclo de clculo descrito a es bsicamente el mismo en problemas lineales y no lineales. Se trata de un procedimiento a muy robusto, eciente y que requiere muy pocos recursos de ordenador, al no ensamblar ninguna matriz. Su unico inconveniente es que el ciclo de clculo es slo condicionalmente estable. Para a o que la integracin sea estable se requiere que: o t tcrit y, en problemas lineales, resulta que: 2 q tcrit = ( 1 + 2 ) !max donde, !max es la mxima frecuencia natural de vibracin del modelo de elementos nitos a o y es la fraccin del amortiguamiento cr o tico a esa frecuencia. En problemas dinmicos reales esta es una restriccin muy fuerte, ya que implica a o que el paso de integracin t debe ser inferior al tiempo que tarda la onda de tensin o o ms rpida en atravesar el elemento nito ms peque~o5 . Sin embargo, cuando se usan a a a n procedimientos de relajacin dinmica, las matrices de masas y de amortiguamiento, M y o a C, pueden ajustarse convenientemente para optimizar el tiempo de clculo. a

En un material como el acero la velocidad de la onda de tensin ms rpida es del orden de 5200 o a a m/s. Esto implica que para un tama~o de elemento de 1 cm, el paso cr n tico de integracin es de unos 2 o microsegundos.

100

Cap tulo 10

Procedimientos de clculo (II) a

10.1 Introduccin o

En los cap tulos anteriores se ha tratado exclusivamente de problemas en los que las diferentes variables de campo dependen de la posicin en el espacio pero no del tiempo. o Sin embargo, las ideas que se han introducido y, en particular, la sistemtica del MEF, es a directamente aplicable tambin a problemas no estacionarios, en los que el tiempo aparece e como variable independiente. La idea bsica es la misma. Dado un dominio de clculo en el que se busca la solucin a a o a lo largo de un intervalo de tiempo, el MEF proporciona la discretizacin espacial del o dominio, esto es, dene los campos solucin en funcin de un nmero nito de parmetros, o o u a por ejemplo, los movimientos nodales. La variacin en el tiempo de la solucin se obtiene o o entonces estudiando la variacin en el tiempo de dichos parmetros. De este modo, el MEF o a permite en general pasar de un sistema de ecuaciones diferenciales en derivadas parciales, en el que aparecen como variables independientes tanto la posicin como el tiempo, a un o sistema de ecuaciones diferenciales ordinarias, en el que la unica variable independiente es el tiempo. Este ultimo sistema se resuelve utilizando un procedimiento de integracin de o diferencias nitas. En las secciones siguientes se pone primero el contexto matemtico y luego se comentan a los procedimientos de integracin en el tiempo ms usuales. Al igual que en cap o a tulos anteriores, se utiliza el problema elstico como veh a culo para expresar ideas que el lector puede fcilmente generalizar a otros problemas. a

10.2

El problema de la elasticidad dinmica a

El problema de la elasticidad dinmica o de la dinmica estructural se plantea del modo a a siguiente. Se tiene un cuerpo elstico que ocupa un volumen 2 <3 y que est limitado por una a a supercie S = Sd [ St, de modo que Sd \ St = ;. La supercie S tiene una normal n ni . El cuerpo tiene una densidad , denida como funcin de punto en , y se encuentra o _ _ en unas condiciones iniciales de posicin u0 u0i y velocidad u0 u0i , denidas tambin o e 101

como campos vectoriales en . Sobre la supercie Sd se imponen desplazamientos d di ; mientras que sobre la t t a u supercie St se imponen tensiones i . Adems, sobre el volumen , acta un campo b bi de fuerzas por unidad de volumen. Estos tres campos son funcin del tiempo en el o intervalo [0; T ] dentro del cual se busca la solucin: o d : Sd ]0; T [ ! <3 : St ]0; T [ ! <3 t b : ]0; T [ ! <3 Las incgnitas del problema son los campos de desplazamientos, u ui ; de deformao ciones, e eij ; y de tensiones, s sij . Todos ellos denidos en ]0; T [. La formulacin fuerte corresponde a un problema de valores iniciales y de contorno, y o se hace en base al siguiente sistema de ecuaciones diferenciales. Un punto sobre el s mbolo de una variable indica derivada con respecto al tiempo. 1. Ecuaciones de compatibilidad en ]0; T [ 1 e = (ru + rt u) 2 o 1 eij = (ui;j + uj;i ) 2

2. Ecuaciones de compatibilidad en Sd ]0; T [ u=d o ui = di

3. Ecuaciones constitutivas (o de respuesta del material) en ]0; T [ s = De o sij = Dijkl ekl

4. Ecuaciones del movimiento en ]0; T [ div s + b = u 5. Ecuaciones de equilibrio en St ]0; T [ sn s n = t 6. Condiciones iniciales en u(0) = u0 u(0) = u0 _ _ o o ui (0) = u0i ui (0) = u0i _ _ o t sij nj = i (10.2) o sij;j + bi = ui (10.1)

102

El punto de partida de los mtodos de elementos nitos es la formulacin dbil del e o e problema anterior. Una de las posibles formulaciones dbiles corresponde a la aplicacin e o del principio de D'Alembert y del principio de los trabajos virtuales para imponer las condiciones de equilibrio dinmico. De esta manera las ecuaciones 10.1 y 10.2 son sustituidas a por una ecuacin integral que impone menos requisitos de derivabilidad. o As una posible formulacin dbil del problema queda como sigue. Dados los espacios , o e St (desplazamientos admisibles) y V (variaciones o desplazamientos virtuales): ~ St = f~ (t) 2 H 1 (); t 2]0; T [ j u(t) = d(t) en Sd g u V = f~ 2 H 1 () j ~ = 0 en Sd g u u encontrar u(t) 2 St de manera que: u u u t u (su ; eu ) = (b; ~ ) ( ; ~ ) + [; ~ ]St u (u(0) u0 ; ~ ) = 0 _ _ u (u(0) u0 ; ~ ) = 0 8~ 2 V u t 2 ]0; T [

Ntese que las variaciones ~ 2 V no dependen del tiempo. o u

10.3

Discretizacin por el MEF o

La discretizacin por el MEF de la forma dbil que se da en la seccin anterior sigue la o e o misma sistemtica utilizada para los problemas estacionarios. a 1. Dividir el dominio de clculo en subdominios o elementos: a =
[
i

\
i

i = ;

o 2. Dentro de cada elemento i , aproximar los campos incgnita que intervienen en las ecuaciones integrales de la forma dbil: e u = N ae (t) _ _ u = N ae (t) u = N ae (t) ~ = N ae u eu = L N ae = B ae eu = L N ae = B ae (t) su = D B ae (t) 103

donde N es una matriz de funciones conocidas, que dependen del espacio y no del tiempo (funciones de forma); ae y ae son vectores de coecientes incgnita, o cuyo valor se hace coincidir normalmente con el valor de los desplazamientos en los nodos; L es una matriz de operadores (derivadas) y D es una matriz de coecientes elsticos. Ntese que, a diferencia de lo que ocurr en los problemas estacionarios, a o a los coecientes ae son ahora funcin del tiempo. o 3. Las integrales que aparecen en las ecuaciones de la forma dbil se hacen sumando e las contribuciones de los distintos elementos, es decir, 8~ 2 V: u u t; u u u 0 = (su ; eu ) (b; ~ ) + ( ; ~ ) [ ~ ]St =
X
i

f(su ; eu )i (b; ~ )i + ( ; ~ )i [; u]Sti g u t ~ u u

y, para las condiciones iniciales: u 0 = (u(0) u0 ; ~ ) = _ _ u 0 = (u(0) u0 ; ~ ) =

X X
i i

u (u(0) u0 ; ~ )i _ _ u (u(0) u0 ; ~ )i

4. La contribucin de cada elemento las sumas anteriores es, sustituyendo las aproxio maciones del paso 2: aei t (Bt D B)i aei aei t (Nt b)i + aei t (Nt N)i aei aei t [Nt]Sti t aei t (Kei aei + Mei aei fei ) donde Kei (Bt D B)i es la matriz de rigidez del elemento i; fei (Nt b)i + [Nt]Sti es el vector de cargas del elemento i; y Mei (Nt N)i es la matriz de t masas del elemento i. En cuanto a las condiciones iniciales, la contribucin de cada elemento es: o aei t (Nt N)i aei0 aei t (Nt u0 )i = aei t(Uei aei0 nuei ) =

aei t (Nt N)i aei0 aei t (Nt u0 )i = aei t(Uei aei0 nvei ) _ _ _ _ donde Uei (Nt N)i , nuei (Nt u0 )i y nvei (Nt u0 )i

5. La suma de las contribuciones de los distintos elementos se hace a travs del proceso e de ensamblaje. Como en el caso de los problemas estacionarios, el proceso de ensamblaje consiste en ir acumulando las contribuciones de cada elemento en matrices y vectores de 104

carcter global, esto es, ordenados segn la numeracin global de grados de libertad. a u o Las ecuaciones discretas de la forma dbil del problema quedan: e
X

elementos

aei t (Kei aei + Mei aei fei ) at (K a + M a f ) = 0

X

elementos elementos
X

aei t (Uei aei0 nuei ) at (Ua0 nu ) = 0 _ aei t (Uei aei0 nvei ) at (Ua0 nv ) = 0 8 a

donde K es la matriz de rigidez global, f es el vector global de cargas y M es la matriz global de masas. 6. Como las relaciones anteriores han de cumplirse 8 a, se llega nalmente al sistema de ecuaciones diferenciales ordinarias: Ma + Ka f = 0 con las condiciones iniciales: a(0) a0 = U1 nu _ _ a(0) a0 = U1 nv el cual permite determinar las funciones incgnita a(t). o 7. De acuerdo con el paso 2, las funciones a(t) denen la solucin dentro de todos o los dominios elementales y, por agregacin, denen la solucin buscada en todo el o o dominio de clculo . a As pues, la aplicacin de la \receta" o sistemtica del MEF conduce al sistema global o a de ecuaciones diferenciales ordinarias 10.3, que hay que integrar en el tiempo a partir de las condiciones iniciales 10.4 y 10.51 . En problemas de dinmica de estructuras aparece normalmente un tercer sumando en a el miembro de la izquierda de 10.3, el correspondiente al amortiguamiento viscoso. El sistema de ecuaciones y sus condiciones iniciales quedan: _ Ma + Ca + Ka f = 0 a(0) = U1 nu a(0) = U _
1

(10.3)

(10.4) (10.5)

(10.6) (10.7) (10.8)

nv

1 Ntese que el sistema 10.3 tiene una estructura similar a la del que gobierna el comportamiento de un o sistema de masas puntuales conectadas por muelles. El papel de la sistemtica del MEF ha sido entonces a el de permitir pasar de un sistema continuo (ecuaciones diferenciales en derivadas parciales denidas en ) a un sistema discreto (masas y muelles) energticamente equivalente. e

105

donde C es la matriz de amortiguamiento. Una hiptesis corriente es suponer que la matriz de amortiguamiento puede obtenerse o como2 : C = M+ K ; 2 <

Esta forma de la matriz C hubiese aparecido en el desarrollo anterior si se sustituyen _ las fuerzas de inercia u por la suma ( u + u), y si la relacin tensin-deformacin, o o o _ _ s = D e, se sustituye por s = D e + D e, donde e 1 (ru + rt u). _ _ 2

10.4

Procedimientos tipo Newmark

El sistema de ecuaciones diferenciales ordinarias 10.6 se conoce con el nombre de sistema de ecuaciones semidiscretas, ya que para llegar a l se ha discretizado el espacio pero no e el tiempo. Para integrar este sistema se utilizan normalmente procedimientos de diferencias nitas, que dan lugar a la discretizacin en el tiempo. Se obtienen as los valores de las o funciones a(t) en instantes de tiempo separados por lo que se conoce como paso de integracin t. o
t -

t0

t1

t2

t3

:::

En general, se utilizan los algoritmos de integracin de tipo Newmark, los cuales se o pueden poner como:

_ Mn+1 + Can+1 + Kan+1 = fn+1 a _ an+1 = an + t an + _ an+1 t2

(10.9)

[(1 2) an + 2 an+1 ] (10.10) 2 _ = an + t [(1 ) an + an+1 ] (10.11)

donde los sub ndices indican el instante de tiempo en que se particularizan las diferentes e o funciones vectoriales (p.ej. an+1 a(tn+1 )). Recurdese que la funcin f (t) es dato. A partir de las relaciones 10.9 a 10.11 se plantea un procedimiento en el que se progresa paso a paso en el tiempo a partir de las condiciones iniciales. Esquemticamente, este a procedimiento de integracin paso a paso podr ser: o a 1. Calcular las aceleraciones iniciales a0 a partir de la relacin 10.9 y de los desplazao mientos y velocidades iniciales: _ a0 = M1 (f0 Ca0 Ka0 )
2

Esta forma de la matriz de amortiguamiento se conoce como amortiguamiento de Rayleigh.

106

2. Inicializar n 0 _ o 3. Poner an+1 y an+1 en funcin de an+1 , utilizando 10.10 y 10.11 _ 4. Sustituir las expresiones an+1 y an+1 de la etapa anterior en 10.9. Se obtiene as un sistema algebraico de ecuaciones lineales que permite obtener los coecientes an+1 5. A partir de an+1 y de los valores ya conocidos del instante de tiempo anterior, an , _ _ an , an , calcular an+1 y an+1 utilizando 10.10 y 10.11 6. Hacer n n + 1 7. Volver a 3 Un procedimiento de integracin tipo Newmark, como el esquematizado en los prrafos o a anteriores, es estable para cualquier tama~o del paso de integracin t si se cumple que: n o 2 1 2

Sin embargo, el procedimiento es slo condicionalmente estable cuando: o 1 2 y > 2

Que el procedimiento es condicionalmente estable signica que el paso de integracin o ha de ser inferior a un valor cr tico: t tcrit para que la integracin numrica de las ecuaciones se mantenga estable. El paso cr o e tico de integracin es funcin tanto de los parmetros y del procedimiento de Newmark como o o a del mallado de elementos nitos y de las propiedades f sicas de los materiales del modelo. De la familia de procedimientos tipo Newmark, los dos ms populares son el proa cedimiento de la aceleracin promedio y el procedimiento de diferencias centrales. Sus o caracter sticas se dan en la tabla. Mtodo e Aceleracin promedio o Diferencias centrales Tipo Impl cito Expl cito
1 4

1 2 1 2

Estabilidad Incondicional Condicional

Orden precisin o 2 2

En la tabla anterior el orden de precisin indica la potencia de t a la cual es proporo cional el error del operador de diferencias nitas. El procedimiento de diferencias centrales es un procedimiento expl cito en el sentido de que si se escoge una matriz de masas M y una matriz de amortiguamiento C que sean diagonales, entonces resulta que el sistema de ecuaciones que se obtiene para progresar con 107

la solucin desde un instante tn al instante siguiente tn+1 , tiene su matriz de coecientes o diagonal. Esto permite obtener la solucin de aceleraciones en tn+1 de manera explcita a o partir de las resultantes de fuerzas internas y externas que actan sobre los nodos en ese u instante. En la prctica, la aceleracin se obtiene en cada nodo simplemente dividiendo a o la resultante de las fuerzas que actan sobre el mismo por la masa nodal. u El algoritmo de integracin expl o cita resulta pues mucho ms sencillo de implementar a que el algoritmo general o impl cito, al no tenerse que ensamblar matrices ni resolver sistemas de ecuaciones con matriz de coecientes llena. Su unico, pero gran inconveniente es que da lugar a un procedimiento de integracin condicionalmente estable, en el que el o paso cr tico de integracin suele ser muy peque~o. En problemas de dinmica estructural o n a este paso cr tico es del orden del tiempo que tarda la onda de tensin ms rpida en o a a atravesar el elemento ms peque~o de la malla3 . En consecuencia, se requiere un nmero a n u enormemente grande de pasos de integracin para completar el estudio de un transitorio o largo. Los procedimientos expl citos tienen su campo natural de aplicacin en los problemas o de impacto, donde los transitorios de inters suelen ser relativamente cortos y donde la e respuesta requiere de por s utilizar pasos de integracin peque~os para denir correcta o n mente su contenido de altas frecuencias. Sin embargo, debido a su robustez para abordar problemas muy no lineales, los algoritmos expl citos se estn empleando cada vez ms en a a otras areas, por ejemplo, en problemas de conformado de metales. Existen as en el mercado dos grandes de familias de programas comerciales de ele mentos nitos para problemas no estacionarios: los programas impl citos y los programas expl citos. En la tabla se resumen las caracter sticas generales de ambos grupos de programas. Debe recordarse siempre que las ecuaciones que resuelven unos y otros son las mismas y que, en principio, cualquier problema puede resolverse empleando ambas clases de programas. La diferencia est en que, segn el caso, una tcnica puede resultar a u e signicativamente ms eciente que la otra. a Impl citos Estables 8t t limitado por la precisin requerida o Ensamblan sistema de ecuaciones Ecientes en dinmica \lenta" a (sismos, viento, oleaje . . . ) Expl citos Debe ser t tcrit No ensamblan (necesitan menos recursos) Ecientes en dinmica \rpida" a a (impactos, explosiones . . . )

3 En un material como el acero, la onda de tensin ms rpida se propaga a una velocidad de unos o a a 5200 m/s. Esto signica que en una malla con elementos de, por ejemplo, 1 cm de lado, el paso cr tico de integracin es del orden de 2 microsegundos. o

108

10.5

Operador de Hilber-Hughes-Taylor

En algunos problemas resulta conveniente introducir en el algoritmo de integracin en el o tiempo un amortiguamiento de naturaleza numrica. Esto ayuda, por ejemplo, a eliminar e el \ruido" en la respuesta asociado a las frecuencias ms altas del modelo de clculo, que a a normalmente son parsitas. a Para ello se utiliza muchas veces una variante del procedimiento de la aceleracin o promedio que se conoce con el nombre de operador de Hilber-Hughes-Taylor. En este operador, la ecuacin 10.9 del algoritmo general de tipo Newmark se sustituye por: o _ _ Mn+1 + (1 + )Can+1 Can + (1 + )Kan+1 Kan = (1 + )fn+1 fn a donde 2 < es un nuevo parmetro tal que: a 1 3 0 1 2 2 (1 )2 4

= =

Este operador da lugar a un algoritmo impl cito incondicionalmente estable y permite introducir una cierta disipacin o amortiguamiento numrico sin prdida de precisin. El o e e o parmetro es el que determina la magnitud de la disipacin. La disipacin es mxima a o o a 1 para = 3 y se anula para = 0. Para = 0 el operador de Hilber-Hughes-Taylor se transforma en el operador de la aceleracin promedio. o

109

Cap tulo 11

Estimacin a posteriori del error o

11.1 Introduccin o

El clculo por el MEF aparece, cada vez ms, integrado en aplicaciones llamadas de a a ingenier asistida por ordenador (CAE), en las que se trata de reunir en un mismo entorno a informtico todas las herramientas necesarias para el proyecto. De este modo, existen a en el mercado paquetes de programas con un mdulo de dibujo asistido por ordenador o (CAD), que sirve para introducir la geometr otro mdulo de clculo (FEA), que evala a, o a u el funcionamiento del prototipo y, a veces, hasta un tercer mdulo que prepara la salida o para las mquinas de control numrico (CAM). La idea es que no tiene mucho sentido a e generar varios conjuntos de modelos, uno para el dibujo o representacin, otro para el o clculo, un tercero para la fabricacin, etc. El objetivo es disponer de una base de datos a o comn a la que acceden las distintas aplicaciones. Es lo que se conoce a veces como u ingenier concurrente. a En un entorno de trabajo \integrado" de este estilo, resulta posible que proyectistas o ingenieros novicios, con muy poco conocimiento del MEF, construyan modelos de clculo y a produzcan resultados. Resultados enmarcados de tal modo por la impresionante potencia grca de los ordenadores actuales, que puede que el usuario no sea consciente siquiera de a que sus comprobaciones tensionales, o de otro tipo, se basan en una aproximacin y de o que la calidad de esta aproximacin es, en principio, desconocida para l. o e As se detecta una cierta tendencia a convertir el MEF en una especie de \caja ne, gra" para la produccin de resultados y esto no est exento de peligro, sobre todo en o a empresas poco tecnicadas u ocinas tcnicas peque~as, donde lo ms probable es que no e n a exista ningn especialista capaz de corregir los abusos a los que se prestan los programas u comerciales. Por otro lado, aun utilizando el MEF de forma consciente y conociendo sus limitaciones, muchas veces resulta dif dar con el mallado ms adecuado o estimar la precisin de los cil a o resultados obtenidos. Esto es especialmente cierto cuando se aborda por primera vez un problema determinado. Las situaciones anteriores tienen un denominador comn: resulta deseable disponer de u medidas cuantitativas de la calidad de las soluciones aproximadas calculadas mediante el

110

MEF. La introduccin de esta clase de medidas permitir establecer un control de calidad o a sistemtico, basado en tolerancias prestablecidas. a El atractivo de esta idea ha hecho que, desde principios de la dcada de los ochenta, se e dedicara un esfuerzo investigador considerable a establecer procedimientos de estimacin o del error en las soluciones que da el MEF. El objetivo es obtener lo que se vino en llamar estimadores a posteriori del error o ndices que utilizan los propios resultados del clculo a para proporcionar una medida de la calidad de los mismos. En la actualidad existen algunos mtodos que permiten obtener aproximaciones muy e razonables al error cometido. En general, dichos mtodos se concentran en el error de e discretizacin que lleva consigo el MEF, ya que se supone que la geometr y las acciones o a exteriores han sido correctamente representadas y no se tienen en cuenta los errores debidos a la aritmtica nita del ordenador. e Los primeros ndices de esta clase aparecieron en los programas comerciales de clculo a hacia el comienzo de la dcada de los noventa. No tienen todav un uso muy extendido e a dentro del conjunto de los usuarios, probablemente porque el usuario medio no sabe interpretarlos bien y porque, salvo en muy pocos casos, no van acompa~ados de procedimientos n automticos de mejora de la solucin para cuando la calidad de la misma no sea aceptable. a o Es decir, el usuario medio no sabe decidir si el ndice de calidad calculado corresponde o no a una calidad suciente y, cuando piensa que la calidad no es suciente, tampoco sabe cmo utilizar los o ndices obtenidos para mejorar su solucin. o La aplicacin ms natural de la estimacin del error son los llamados procesos de o a o clculo adaptables, en los que, partindose de una solucin de calidad muy baja, sta se a e o e va enriqueciendo por sucesivos renamientos hasta que las medidas de calidad cumplen unos requisitos jados de antemano. Son este tipo de procesos los ms apropiados para ser a incluidos casi como cajas negras dentro de un paquete integrado de CAE. Su \adaptividad" hace que las soluciones que proporcionan sean relativamente insensibles aun a graves errores en la preparacin del modelo inicial. Errores tales que podr hacer inservibles o an los resultados de un clculo convencional. a Se trata en este cap tulo de introducir primero ciertos conceptos bsicos, para luego a presentar alguno de los estimadores de error disponibles y el esquema bsico de un proceso a de clculo adaptable. a

11.2

Conceptos bsicos a

El primer concepto que debe introducirse es el de error de discretizacin. Este error se o dene como: ~ er u u en

donde u es la solucin exacta para el campo incgnita bsica, u es la aproximacin obtenida o o a ~ o por el MEF y es el dominio de denicin del problema. El error de discretizacin es pues o o una funcin de punto y no es calculable, ya que su conocimiento implicara el conocimiento o de la solucin exacta. Las medidas de la calidad se basan entonces en estimaciones de la o solucin exacta u o en las manifestaciones del campo de errores er . o 111

La magnitud del error se dene habitualmente como una norma del campo er . Ejemplos de normas son la norma L2 : ker k2 2 L
Z

et er d r

y la norma energtica asociada a la forma dbil punto de partida de las ecuaciones del e e MEF. En un problema de elasticidad la norma energtica es: e ker k2
Z

er : "er d

donde er y "er son los campos de tensiones y deformaciones derivados del campo de desplazamientos er . ~ Normalmente, ante una solucin aproximada u, se tiene inters en dos aspectos del o e error de discretizacin cometido. Por un lado, interesa tener una medida global de la o calidad de la solucin aproximada. Esto permite decidir sobre si la solucin es aceptable o o o no. Por otro lado, tambin es importante conocer la distribucin del error a lo largo del e o dominio, ya que esto se~ala las zonas en que debe incidirse para mejorar la calidad global n de la solucin. o As se suele llamar estimador de error a cualquier medida que aproxime la magnitud , del campo de errores en el dominio del problema y que, por tanto, proporcione un ndice global de la calidad de la misma. Por contraste, se llama indicador de error a cualquier medida sobre la distribucin del error a lo largo del dominio de clculo y que, por tanto, o a proporcione informacin para renar la solucin de modo eciente. o o Los estimadores de error son expresin directa de la calidad de una solucin aproo o ximada, mientras que los indicadores proporcionan un mapa de la distribucin relativa o del error. Ambas categor no son excluyentes. Por ejemplo, medidas locales absolutas as del error pueden emplearse como indicadores de cara al renamiento de una malla de elementos nitos.

11.3

Indicadores y estimadores de error

En esta seccin se pasa revista a algunos de los indicadores y estimadores de error dispoo nibles actualmente.

11.3.1

Densidad de energ de deformacin (SED) a o

Los gradientes de la densidad de energ de deformacin (\Strain Energy Density" - SED), a o entendiendo este concepto en sentido generalizado, se han utilizado mucho para identicar las zonas del mallado en donde es ms conveniente el renamiento, esto es, como a indicadores de error. La idea es la siguiente. El MEF es un procedimiento optimo en sentido energtico. e Tiende a minimizar la diferencia entre la norma energtica exacta y la aproximada: e

112

kuk2 k~ k2 = u

Z |

: " d

u : "u d

en elasticidad Por lo tanto, se obtendrn mejores soluciones cuanto mejor se represente la norma energtica a e de la solucin exacta. As la solucin aproximada ser tanto mejor cuanto mejor se aproo , o a xime la funcin subintegral en el clculo su norma energtica. En problemas de elasticidad o a e esta cantidad subintegral es una densidad de energa de deformacin u : "u . En general, o har falta mayor discretizacin en una zona cuanto mayor sea el gradiente de densidad a o de energ de deformacin en ella. Por ejemplo, con elementos de deformacin constante, a o o esto signica incrementar el nmero de los mismos en las zonas con mayores saltos en la u densidad de energ al pasar de un elemento a otro. a Mediante esta tcnica, relativamente sencilla de implementar, se obtiene informacin e o sobre dnde resulta ms conveniente renar la discretizacin, pero no un o a o ndice de calidad global de la solucin. Es decir, se obtiene un indicador de error, pero no un estimador, en o el sentido explicado en la seccin anterior. o

{z

11.3.2

Normas de residuos

~ La solucin aproximada u se obtiene por el MEF a partir de la forma dbil de las ecuaciones o e ~ que rigen el problema. Como u no es la solucin exacta, si se introduce en las ecuaciones de o la forma fuerte, no las cumplir y dar lugar a una funcin de punto o campo de residuos a a o ru denido en el dominio . Por ejemplo, en el caso del problema elstico, el campo de residuos se obtiene al a sustituir u en la forma fuerte de las ecuaciones de equilibrio: ~ = div u + b = ru 60 u r : ! <3 El campo de residuos es una manifestacin o consecuencia de que exista error en o la solucin aproximada. Se han propuesto entonces diversos estimadores de error que o intentan acotar superiormente la norma energtica del campo de errores ker k utilizando e normas del campo de residuos. Dichos estimadores se encuentran entre los mejores de que se dispone actualmente. Su principal inconveniente es el enorme costo en tiempo de clculo que requiere el clculo de los residuos, que puede ser del orden del necesario para a a la construccin de la matriz de rigidez global del problema. o Por otro lado, estos estimadores se construyen en base a la suma de contribuciones locales, que pueden ser empleadas a su vez como indicadores de la distribucin del error, o para gu un proceso de renamiento de la solucin aproximada. ar o

11.3.3

Estimadores del tipo Z 2

Los estimadores del tipo Z 2 aparecieron como un intento de reducir los tiempos de clculo a asociados a la estimacin del error basada en normas de residuos y, sobre todo, con la idea o 113

de una fcil implementacin dentro de la estructura habitual de los programas de clculo. a o a Su fundamento es la evaluacin de la diferencia entre las tensiones que resultan del o clculo con elementos nitos y una proyeccin continua o \alisada" de dichas tensiones, la a o cual se postula que tiene un orden de precisin superior. Entindase aqu el concepto de o e \tensin" en el sentido generalizado de un \gradiente del campo incgnita bsica". o o a Se han propuesto diversos procedimientos para la obtencin de la proyeccin continua o o de las tensiones calculadas por el MEF en una solucin alisada. El ms sencillo consiste o a en extrapolar el valor de la tensin a los nodos, promediar en cada nodo el valor de la o tensin suministrado por los elementos que se conectan al mismo y utilizar estas medias o nodales y las funciones de forma N para denir el campo alisado de tensiones. De este modo, dentro de cada elemento el campo alisado es: = N e donde e representa el vector de promedios nodales de tensin. o denido de esta manera puede postularse que es una mejor aproximacin o El campo de la solucin exacta en tensiones que el campo u que se deriva directamente de la o solucin calculada en desplazamientos u. Resulta razonable entonces aproximar el error o ~ en tensiones como la diferencia: u u y entonces estimar la norma energtica del error como la energ de deformacin asociada e a o a esta diferencia. Esto es: ker k2
Z

( u )t D1 ( u ) d

donde D es una matriz de coecientes elsticos. La integral anterior se calcula como suma a de integrales extendidas a los dominios elementales y entonces cada sumando puede ser utilizado como indicador local del error en el elemento correspondiente.

11.3.4

Estimadores basados en diferencias entre funcionales

Cuando un problema admite una formulacin variacional, normalmente pueden utilizarse o diferentes principios variacionales para resolverlo. Los funcionales asociados a estos principios se derivan unos de otros utilizando la tcnica de los multiplicadores de Lagrange. e De un funcional a otro cambia el nmero de campos independientes y cambian tambin u e las restricciones que delimitan los espacios de funciones de los que se toman estos campos. Sin embargo, los funcionales gozan de la propiedad de que, si se sustituyen en ellos los campos solucin exacta, entonces todos adoptan el mismo valor. Esto es, la solucin o o exacta del problema hace que todos los funcionales valgan lo mismo. As las diferencias entre los valores de los funcionales para una solucin aproximada , o constituyen una manifestacin del error de discretizacin y puede ser aprovechada para o o estimarlo. Frente a estimadores de otras clases, los estimadores basados en diferencias entre funcionales tienen la ventaja de que se evalan elemento a elemento de forma independiente, u 114

sin ser necesario el clculo de promedios nodales ni de saltos interelementales. Desde el a punto de vista terico, no requieren para su aplicacin ninguna condicin de suavidad o o o o lisura en la solucin exacta, la cual puede exhibir discontinuidades. o Por contra, la evaluacin de esta clase de estimadores requiere disponer de varios o campos independientes, para poder evaluar al menos un funcional multicampo. Esto hace que no puedan utilizarse con facilidad en elementos basados en la formulacin convencional, o donde slo existe un campo independiente. o El punto clave de un estimador de esta clase es la eleccin de los funcionales por cuya o diferencia se estima el error. Si estos funcionales se escogen de manera que converjan hacia su valor exacto uno por la derecha y otro por la izquierda, entonces su diferencia ser una a cota superior del error en la energ y, por tanto, un estimador excelente. a Como los estimadores se calculan elemento a elemento, las contribuciones elementales pueden utilizarse como indicadores de error para gu el renamiento de la solucin. ar o

11.4

Procesos adaptables

Aun cuando la estimacin del error cometido en un anlisis determinado y el conocimiento o a de su distribucin a lo largo de la malla empleada son aspectos esenciales, stos no pueden o e considerarse ms que un primer paso. Conocida la calidad de la solucin aproximada, el a o ingeniero debe poseer las herramientas necesarias para poder renar esta solucin de modo o eciente en el caso de que no se considere aceptable la calidad obtenida. Este es precisamente el papel de los procesos adaptables, los cuales tienen como objetivo proporcionar de modo ms o menos automtico soluciones acordes con la calidad exigida. a a El diagrama de ujo de un proceso adaptable es el siguiente:

Solucin o inicial
@ Calidad @ - @ @ suciente? @ @@ ? @

Fin

No
?

Enriquecer discretizacin o

Solucin mejorada o

115

El nivel de calidad de la solucin se evala mediante un estimador del error global. o u Si el nivel de calidad no se considera suciente, la discretizacin se enriquece tratando de o uniformizar la distribucin del error a lo largo de la malla que dan los indicadores de error. o Existen dos formas bsicas de enriquecer una malla de elementos nitos, esto es, de a introducir ms grados de libertad en aquellas zonas que se identican como generadoras a de ms error. En un clculo convencional, en el que se utiliza un tipo de elemento detera a minado, el unico medio de enriquecer la malla es subdividiendo los elementos en elementos ms peque~os. Esta clase de renamiento se conoce con el nombre de renamiento h. a n Hay otra alternativa, que es hacer ms rica la discretizacin incrementando el orden de a o los polinomios que constituyen las funciones de aproximacin o funciones de forma dentro o de cada elemento. Es lo que se conoce como renamiento p. De esta manera, la malla permanece siempre con el mismo nmero de elementos. u Se puede pensar tambin en formas mixtas de renamiento h-p, las cuales combinan e las dos clases de renamiento anteriores.

116

Cap tulo 12

Conceptos bsicos de la Mecnica a a de Slidos o

12.1 Introduccin o

Se inician con este cap tulo los temas dedicados a la aplicacin del MEF a problemas no o lineales. Se dice que un problema es lineal cuando existe proporcionalidad entre las acciones y las respuestas, es decir, si la respuesta es una funcin lineal de las acciones. Un problema o es no lineal en cualquier otro caso. La aplicacin del MEF a problemas no lineales es una generalizacin de las ideas deo o sarrolladas para problemas lineales. El principio bsico de discretizacin del dominio de a o clculo es el mismo, aunque la naturaleza de los problemas da lugar ahora a sistemas a discretos no lineales, ms dif a ciles de tratar numricamente y para los cuales se han ido e desarrollando algoritmos y procedimientos de clculo mucho ms complejos que los utilia a zados en problemas lineales. Sin embargo, a la hora de abordar la aplicacin del MEF a esta clase de problemas o el primer escollo con que se encuentra el estudioso no es la complejidad de los algoritmos de clculo, sino que la teor bsica de su problema, aquella que es previa a la aplicacin a a a o del MEF, es tambin mucho ms compleja que la que hay detrs de los problemas lineales e a a convencionales y, adems, no le resulta familiar. a El aumento de la potencia de los ordenadores y el desarrollo paralelo del MEF ha permitido utilizar intensivamente conceptos y teor que, hasta hace pocos a~os, hab as n an tenido poca aplicacin prctica y que, por tanto, hab recibido en general una atencin o a an o relativamente escasa. De este modo, se da la paradoja de que el MEF pone al alcance del ingeniero la aplicacin prctica de herramientas tericas muy poderosas para describir la o a o realidad f sica, pero el ingeniero se encuentra con que no ha recibido suciente formacin o bsica acerca de esas teor y se siente incmodo utilizndolas. a as o a Intentando sortear esta dicultad, se introducen en este captulo algunos conceptos bsicos de la Mecnica de Slidos, los cuales resultan necesarios para el estudio de proa a o blemas no lineales dentro de esta disciplina. Estos conceptos son previos a la aplicacin o

117

del MEF pero, por razones obvias, deben conocerse antes de emplear sta o cualquier otra e tcnica numrica, ya que forman parte de la propia denicin del problema. e e o Una vez ms se ha elegido la Mecnica de Slidos como veh a a o culo para explicar la aplicacin del MEF, por ser ste el campo de utilizacin ms corriente. Los interesados o e o a en otras disciplinas debern hacer un esfuerzo de abstraccin an mayor que el exigido en a o u los cap tulos anteriores.

12.2

Fuentes de no linealidad en Mecnica de Slidos a o

En Mecnica de Slidos y Estructuras existen cuatro fuentes de no linealidad: a o 1. Comportamiento del material El comportamiento del material deja de ser lineal cuando las tensiones ya no son proporcionales a las deformaciones. Esto ocurre, por ejemplo, cuando se utilizan modelos de plasticidad, hiperelasticidad o da~o para representar la respuesta del n material. 2. Magnitud de la deformacin o La cinemtica del movimiento del slido o de la estructura puede hacer que las a o deformaciones dejen de ser una funcin lineal de los desplazamientos. Esto ocurre o cuando las acciones dan lugar a movimientos que ya no son \peque~os" con respecto n a las dimensiones del dominio de clculo. Entonces, adems, el propio dominio de a a clculo cambia a lo largo del problema. a 3. Condiciones de contorno Las restricciones al movimiento pueden alterarse con la aplicacin de las acciones, o dando lugar tambin a una falta de proporcionalidad en la respuesta. Es el caso de e los problemas en los que el contacto entre diferentes componentes se establece o se pierde a lo largo de la historia de cargas. 4. Fuerzas aplicadas Las fuerzas aplicadas pueden ser tambin una fuente de no linealidad cuando stas e e se vinculan a la deformacin del dominio. Es por ejemplo el caso de una presin o o constante aplicada sobre un area creciente con la deformacin o el de una fuerza o puntual que sigue una direccin que var con la deformacin de la estructura. o a o Un problema deja de ser lineal cuando en l se maniesta una o varias de las causas e anteriores.

12.3

Tensor gradiente de deformacin o

En un problema estructural general, el analista describe la conguracin inicial de la o estructura y busca su movimiento o deformacin a lo largo de la historia de cargas. Si o 118

X 2 <3 representa la conguracin inicial y t 2 < es el tiempo o un parmetro de carga, o a se busca la conguracin x 2 <3 en el instante denido por t: o x = x(X; t) (12.1)

Es decir, una part cula A de material inicialmente en la posicin X ocupar una nueva o a posicin x cuando el parmetro de carga valga t. El problema general consiste en encontrar o a la denicin de x en funcin de t. o o t = t1

t = t0
3 dx r A x + dx * x

A0r

0
dX X + dX r A X 6

A0r

x = x(X; t)

Si se supone que el material no se crea ni se destruye durante el movimiento y que el movimiento no da lugar tampoco discontinuidades (grietas) dentro del dominio de clculo, a entonces resulta que la correspondencia denida por 12.1 es biunvoca y que la relacin o 12.1 puede ser invertida. Esto es, que 12.1 dene tambin la posicin inicial X a partir de e o la posicin actual x y del parmetro de carga t. o a En el entorno de la part cula A el movimiento del medio continuo tiene dos componentes, una componente de slido r o gido, que no da lugar a una respuesta tensional del material, y una componente de deformacin, a la cual se oponen fuerzas internas derivadas o de la estructura del material. En el estudio del movimiento conviene pues separar estas dos componentes, ya que movilizan respuestas muy diferentes. Esto se consigue analizando la evolucin de dos part o culas prximas A y A0 . o 119

En la conguracin inicial la diferencia de posicin entre ambas part o o culas vendr a dada por AA0 = dX. Mientras que en la conguracin actual la diferencia de posicin o o o ser: AA0 = dx. La relacin entre dX y dx puede obtenerse a partir de 12.1: a @x dX F dX (12.2) @X donde F es un tensor de segundo orden que se conoce con el nombre de tensor gradiente de deformacin. El tensor F es una funcin de punto. o o dx =

12.4

Teorema de descomposicin polar o

El estudio del tensor gradiente de deformacin F es el que permite separar el movimiento o del entorno de un punto en la parte que da lugar a deformacin y en la que slo es un o o movimiento de slido r o gido. La longitud inicial dL del segmento AA0 viene dada por: dL2 = dXt dX y la longitud actual dl: dl2 = dxt dx Se dene entonces el cociente de alargamiento como: dl = dL
s

dxt dx dXt dX

El sentido f sico del cociente de alargamiento es claro. Si = 1, entonces no ha habido deformacin en torno a A en la direccin marcada por dX en la conguracin inicial, slo o o o o movimiento de slido r o gido. Como a partir de 12.2 resulta que: dl2 = dxtdx = dXt Ft FdX el cociente de alargamiento se puede poner como: dXt dXt Ft FdX dXt =p Ft F p = Nt Ft FN (12.3) dXt dX dXt dX dXt dX donde N es un vector unitario en la direccin de dX. o La relacin anterior permite, dada una part o cula A situada inicialmente en la posicin o X, medir el cociente de alargamiento segn cualquier direccin N que salga de ella. En u o cada punto A el cociente de alargamiento es pues una funcin de la direccin N. Se o o puede entonces variar N buscando puntos estacionarios (mximos o m a nimos) de . Esto es: 2 = Hacer estacionario: 2 (N) = Nt Ft FN 120 con la condicin o NtN 1 = 0

Utilizando la tcnica de los multiplicadores de Lagrange, los puntos estacionarios para e se obtienen de: (N; p) = Nt Ft FN p(Nt N 1) @ = Ft FN pN = 0 @N @ = Nt N 1 = 0 @p (12.4) (12.5) (12.6)

donde p es un multiplicador de Lagrange que, premultiplicando 12.5 por Nt y teniendo en cuenta 12.3, resulta ser igual a 2 . De este modo, segn 12.5 los puntos estacionarios de cumplen: u (FtF 2 I)N = 0 (12.7)

La relacin 12.7 dene un problema de autovalores que, por ser Ft F una matriz o simtrica de orden 3, proporciona tres direcciones NI , NII y NIII , ortogonales entre e s llamadas direcciones principales de deformacin. Y proporciona tambin tres cocientes , o e de alargamiento, I , II y III , segn esas direcciones, llamados alargamientos principales. u En la conguracin actual, denida por x, a los tres vectores unitarios NI , NII y NIII o les corresponden otros tres versores nI , nII y nIII , de modo que: I nI II nII III nIII = FNI = FNII = FNIII

y los tres versores nI , nII y nIII forman tambin un triedro ortogonal, ya que: e nt nII = I II t 1 1 t t NI F FNII = N NII = 0 I II I I

e igualmente para los productos nt nIII y nt nIII . I II e a De este modo, al formar (nI , nII , nIII ) un triedro ortogonal, resulta que ste podr obtenerse mediante una rotacin de slido r o o gido a partir del triedro (NI , NII , NIII ). Sea R la transformacin que representa esta rotacin, entonces: o o nI nII nIII = RNI = RNII = RNIII

y se cumplir que Rt = R1 . a Un vector dXI segn NI en la conguracin inicial se transforma en el vector dxI en u o la conguracin actual, tal que: o 121

dxI = I R dXI y anlogamente: a

dxII dxIII

= II R dXII = III R dXIII

Cualquier vector dX en la conguracin inicial se puede poner como: o dX = dXI + dXII + dXIII = (Nt dX)NI + (Nt dX)NII + (Nt dX)NIII I II III por ser NI , NII , NIII un triedro ortonormal, con lo que se tiene: dx = dxI + dxII + dxIII = I R dXI + II R dXII + III R dXIII = I R(Nt dX)NI + II R(Nt dX)NII + III R(Nt dX)NIII I II III = I RNI (Nt dX) + II RNII (Nt dX) + III RNIII (Nt dX) I II III = R(I NI Nt + II NII Nt + III NIII Nt ) dX I II III

La relacin anterior da lo que se conoce como descomposicin polar del tensor gradiente o o de desplazamiento F. El tensor de segundo orden: U I NI Nt + II NII Nt + III NIII Nt I II III es el tensor de alargamiento por la derecha y, con l, la descomposicin polar de F se e o escribe: dx = F dX = RU dX (12.8)

Mediante un desarrollo similar al anterior puede llegarse tambin a una descomposicin e o polar de la forma: dx = F dX = VR dX donde el tensor V es el tensor de alargamiento por la izquierda: V I nI nt + II nII nt + III nIII nt I II III El sentido f sico que tiene la descomposicin polar 12.8 es que el movimiento en el o entorno de un punto A dado por F puede descomponerse en dos transformaciones sucesivas, un alargamiento puro en tres direcciones ortogonales, seguido de una rotacin pura de o slido r o gido. En la descomposicin 12.9 sucede al revs, se produce primero la rotacin o e o pura y luego el alargamiento segn las tres direcciones ya giradas. u 122 (12.9)

12.5

Medidas de la deformacin o

Las medidas de deformacin utilizadas en problemas generales se entienden mejor si se o considera primero el concepto de deformacin en una dimensin y luego se generaliza ste o o e a tres dimensiones utilizando el teorema de descomposicin polar. o En una dimensin se ha introducido ya una medida de la deformacin, el cociente de o o alargamiento . En un problema de estiramiento unidimensional con deformacin uniforme o resulta que: = L L0

donde L es la longitud actual y L0 la longitud inicial. Sin embargo, esta medida no es muy util, por ejemplo, para estudiar la deformacin o elstica de un componente de acero. En un acero estructural tpico el mdulo de elasticidad a o vale aproximadamente 210 GPa y la tensin de uencia, 260 MPa. Entonces, el valor del o cociente de alargamiento al alcanzarse el l mite de uencia ser: a L0 (1 + 0:260 ) L0 + L 210 = = = 1:0012 L0 L0 Es decir, el primer d gito de inters se encuentra en la cuarta cifra signicativa y eso no e es deseable. Para evitar este inconveniente, se denen en general medidas de deformacin como o funciones f del cociente de alargamiento que cumplen las tres condiciones siguientes: " = f() 1. La deformacin es nula cuando el cociente de alargamiento es igual a la unidad: o f (1) = 0 2. La derivada de la deformacin es igual a uno cuando el cociente de alargamiento es o igual a la unidad: df (1) = 1 d 3. La deformacin es creciente con el alargamiento y a cada valor de le corresponde o estrictamente un valor de la deformacin: o df >0 d 8 > 0

123

Las dos primeras condiciones tienen por objeto garantizar que todas las medidas de deformacin sean iguales para \peque~as" deformaciones ( prximo a 1). En efecto, o n o desarrollando " en serie de Taylor alrededor de = 1, se tiene: " = f(1) + ( 1) df 1 d2 f (1) + ( 1)2 2 (1) + : : : d d 2

que para prximo a 1 da, teniendo en cuenta las condiciones impuestas a f: o " ( 1) lo cual constituye la denicin de deformacin a la que estar ms habituado el lector. o o a a En cuanto a la tercera condicin, sta es un poco arbitraria, ya que podr elegirse o e a tambin la condicin contraria: e o df <0 d 8 > 0

lo que dar lugar a deformaciones positivas al acortarse el material. Esto es lo ms cmodo a a o en problemas de Mecnica de Suelos. a Cumpliendo las tres condiciones anteriores existen muchas posibilidades para elegir la funcin f medida de deformaciones. Las tres ms corrientes son las que siguen: o a 1. Deformacin nominal o deformacin de Biot. o o fN () = ( 1) Llamada tambin deformacin ingenieril. Es la ms intuitiva. En un problema de e o a estiramiento unidimensional con deformacin uniforme vale: o fN = 2. Deformacin logar o tmica. fL () = ln() Esta medida es de uso corriente en problemas de plasticidad de metales ya que, cuando se representan las curvas de tensin-deformacin equivalente correspondientes a o o los ensayos de traccin simple, de compresin simple y de torsin pura, se obtieo o o nen curvas casi coincidentes si se emplean medidas de deformacin logar o tmica para calcular la deformacin equivalente. En un problema de estiramiento unidimensional o con deformacin uniforme la deformacin logar o o tmica vale: fL = ln( L ) L0 L L0 L0

124

3. Deformacin de Green. o 1 fG () = (2 1) 2 La generalizacin de esta medida a problemas tridimensionales es muy fcil de obo a tener a partir del tensor gradiente de deformacin. Aprovechando la condicin de o o que todas las medidas de deformacin adoptan los mismos valores para \peque~as" o n deformaciones, esta medida se ha utilizado mucho en problemas de \grandes" desplazamientos pero \peque~as" deformaciones. En un problema de estiramiento unin dimensional con deformacin uniforme la deformacin de Green vale: o o 1 L2 fG = ( 2 1) 2 L0 La eleccin de una medida de deformacin u otra es un asunto estrictamente de comoo o didad de empleo, a partir de dos consideraciones: 1. La facilidad con que la deformacin puede ser calculada a partir de los desplazao mientos, ya que stos ultimos suelen ser las variables bsicas dentro de los modelos e a de clculo. a 2. Lo apropiado de la medida de deformacin para servir de base a un modelo de como portamiento del material, ya que la deformacin es el v o nculo entre la cinemtica a o movimiento del continuo y la respuesta del material. En este sentido, por ejemplo, parece que la deformacin logar o tmica es la ms adecuada para modelos de a plasticidad de metales. Habiendo denido la idea bsica de medida de deformacin en una dimensin, sta se a o o e puede generalizar a tres dimensiones haciendo uso del teorema de descomposicin polar. o Alrededor del punto material A la parte del movimiento del continuo que da lugar a deformacin se caracteriza completamente mediante seis variables: los tres cocientes de o alargamiento principal, I , II , III , y la orientacin de las tres direcciones principales de o alargamiento en la conguracin de referencia, NI , NII , NIII . o Se pueden denir entonces las deformaciones principales como: "I = f(I ) "II = f(II ) "III = f(III ) y el tensor de deformaciones como: " "I NI Nt + "II NII Nt + "III NIII Nt I II III Este tensor caracteriza completamente el estado de deformacin en el entorno del punto o material. 125

12.6

Tasa o velocidad de deformacin o

Muchos materiales con inters prctico presentan una respuesta tensional que, en cada inse a tante, depende del camino de deformacin que se haya recorrido. En estos casos el modelo o matemtico de la ley de comportamiento del material suele formularse de forma incremena tal, es decir, se relaciona un incremento de tensin con el incremento de deformacin. Es o o el caso, en general, de los modelos de plasticidad. La denicin incremental de la deformacin se hace a partir del tensor velocidad de o o deformacin, tambin llamado tasa de deformacin1 . o e o La velocidad o tasa de cambio de posicin v de un punto material A es: o v= @x @t

o Entonces, la diferencia de velocidad entre dos puntos vecinos A y A0 en la conguracin actual es: dv = @v dx = L dx @x

en donde L recibe el nombre de tensor gradiente de velocidades. De este modo se tiene, por un lado: dv = L dx = L F dX y, por otro lado: dv = @(dx) @(F dX) _ = = F dX @t @t

luego el tensor gradiente de velocidades L cumple que: _ LF = F (12.10)

El tensor gradiente de velocidades L se descompone en la suma del tensor velocidad de deformacin D y el tensor velocidad de giro : o @v 1 1 @v + ( )t ) (L + Lt ) = ( @x 2 2 @x 1 (L Lt ) 2

D = =

Es el tensor D el que se emplea para denir el comportamiento de los materiales segn u modelos matemticos de naturaleza incremental. a
1

En la terminolog anglosajona este tensor se conoce con el nombre de \strain rate tensor". a

126

12.7

Equilibrio y trabajo virtual

El problema del clculo de tensiones en slidos y estructuras se plantea como encontrar la a o solucin de movimientos, deformaciones y tensiones en un dominio de conguracin inicial o o 3 0 2 < que es sometido a una cierta historia de carga. En este contexto se entiende como \carga" una serie de sucesos para los cuales se busca la respuesta del slido o estructura. o Tanto en problemas lineales como no lineales, el punto de partida clsico para la aplia cacin de las tcnicas de elementos nitos es la forma dbil de las ecuaciones de equilibrio: o e e el principio de los trabajos virtuales. El principio de los trabajos virtuales tiene una naturaleza muy general y resulta vlido a como expresin de las condiciones de equilibrio independientemente de la linealidad o o no del problema. Simplemente, en problemas geomtricamente no lineales, se debe ser e cuidadoso con los dominios de integracin. Con la notacin empleada en los cap o o tulos anteriores, el principio de los trabajos virtuales puede enunciarse como: Dados los espacios S (desplazamientos admisibles) y V (variaciones o desplazamientos virtuales): ~ S = f~ 2 H 1 () j u = d en Sd g u V = f~ 2 H 1 () j ~ = 0 en Sd g u u Encontrar u 2 S, de manera que: u t u ( u ; eu ) = (b; ~ ) + [; ~ ]St 8~ 2 V u (12.11)

En el caso de problemas geomtricamente no lineales hay que ser consciente de que los e dominios , Sd y St corresponden a la conguracin \actual" del slido o estructura y que o o corresponde a la tensin verdadera o de Cauchy. Tambin est impl o e a cito en la ecuacin o del trabajo virtual 12.11 que el desplazamiento u tiene una magnitud \peque~a", para n no alterar los dominios de integracin. o Precisamente, para no dejar lugar a dudas respecto a esto ultimo, se suele trabajar con velocidades virtuales v en vez de con desplazamientos virtuales u. Si el problema es esttico, el tiempo en el que se basa la denicin de velocidad virtual es un tiempo cticio, a o asociado al avance por el proceso de carga. Matemticamente, se trata simplemente de un cambio de notacin, desde u a v. La a o ecuacin 12.11 se escribe ahora: o ( u ; Dv ) = (b; ~ ) + [ ~ ]St v t; v 8~ 2 V v (12.12)

Desde el punto de vista f sico, la ecuacin 12.12 se lee como que la potencia desarrollada o por las fuerzas interiores con la velocidad de deformacin virtual es igual a la potencia o desarrollada por las fuerzas exteriores en el campo de velocidad virtual. 127

12.8

Medidas de tensin o

La ecuacin de los trabajos virtuales constituye la expresin \dbil" de las condiciones de o o e equilibrio y es el punto de partida de la aplicacin del MEF. o En la ecuacin de los trabajos virtuales 12.12, el trmino correspondiente a la potencia o e de las fuerzas interiores se expresa en funcin del tensor de tensiones verdaderas o tensiones o de Cauchy, , y del tensor de velocidad de deformacin virtual, D, es decir: o ( u ; Dv )
Z

u : Dv d

donde u representa las tensiones actuales en el slido y Dv la velocidad de deformacin o o virtual. Ambos campos denidos en el volumen actual . Resulta a veces conveniente, o ms cmodo, realizar la integral anterior sobre el dominio a o 0 que corresponde al volumen original o de referencia del slido. Para ello se realiza el o cambio de variables denido por la transformacin 12.1. Es decir: o
Z

@x (x) : D (x) d = (X) : D (X) det( ) d0 = @X 0

u v u v

J u (X) : Dv (X) d0

donde J es el determinante del jacobiano de la transformacin 12.1: o J det( d @x ) = det(F) = @X d0

El tensor J u es una medida de tensin \cticia" que se conoce con el nombre de o tensor de tensiones de Kirchho. Por otro lado, una medida de deformacin que se utiliza mucho, por lo fcil que resulta o a su clculo, es la deformacin de Green. El tensor de deformaciones de Green es: a o 1 "g (Ft F I) 2 Puede entonces resultar cmodo calcular la integral correspondiente a la potencia de o las fuerzas interiores utilizando la tasa o velocidad de deformacin de Green: o 1 _ _ _ "g = (Ft F + Ft F) 2 A partir de 12.10, se tiene que: 1 "g = Ft(Lt + L)F = Ft D F _ 2 o _ D = Ft "g F1 Entonces, en la ecuacin de los trabajos virtuales, se tendr: o a

128

: D d =

J : D d0 =

_ J : (Ft "g F1 ) d0 =

_ J(F1 Ft ) : "g d0

(12.13) En donde se ha hecho uso de que si A, B y C son tres matrices cuadradas, se cumple: (A B) : Ct = (C A) : Bt = (B C) : At El tensor de tensiones que aparece en el ultimo trmino de 12.13 es otra medida cticia e de tensin que se conoce con el nombre de segundo tensor de tensiones de Piola-Kirchho: o S J(F1 Ft ) Ntese que en la ecuacin de los trabajos virtuales deben emplearse medidas conjugao o das de tensin y deformacin para el clculo de la potencia desarrollada por las fuerzas o o a interiores en el campo de velocidades virtual. En este sentido, dependiendo de la medida de velocidad de deformacin que se utilice, debe tomarse una u otra medida de tensin. o o Los casos analizados en esta seccin son los ms corrientes en la prctica y se resumen en o a a la tabla siguiente. Deformacin o Velocidad de deformacin D o Velocidad de deformacin D o _ Tasa de deformacin de Green "g o Tensin conjugada o Verdadera Kirchho J 2o Piola-Kirchho S Dominio de integracin o Actual De referencia 0 De referencia 0

12.9

Tasa o velocidad de cambio de la tensin o

Como ya se ha se~alado, muchos modelos constitutivos o de comportamiento del material n se formulan de manera incremental. En esta clase de modelos se obtiene un incremento de tensin en funcin de un incremento de deformacin. o o o Algunas medidas de tensin, como el tensor de tensiones verdaderas o el tensor de o Kirchho, estn asociadas a direcciones espaciales dentro de la conguracin actual . a o Entonces resulta necesario separar el incremento de tensin en dos partes, una debida o a la variacin de las direcciones espaciales (por el movimiento del slido o cambio en la o o conguracin actual) y otra ocasionada por la respuesta del material al incremento de o deformacin. La primera se deriva de la cinemtica del movimiento y la segunda es la que o a proporciona el modelo de comportamiento del material. La diferencia entre ambas componentes del incremento de tensin se visualiza mediante o el siguiente ejemplo. Sea una barra de seccin A sometida a una fuerza axil P constante o de traccin en cada uno de sus extremos. La barra gira alrededor de su punto medio, de o manera que en el instante t = t1 se encuentra paralela al eje X y en el instante t = t2 es paralela al eje Y . En el instante t = t1 la tensin verdadera vale xx = P y yy = 0. o A Mientras que en el instante t = t2 vale xx = 0 y yy = P . Resulta entonces que en el A 129

= intervalo [t1 ; t2 ] ha habido un cambio de las componentes del tensor de tensiones: dxx 60 y dyy 60. Sin embargo, esta variacin de las tensiones tiene que ver con el movimiento = o de slido r o gido de la barra y no con la respuesta del material a la deformacin. o t = t2 P
6

t = t1

Y
6 ? -

Sea e , = 1; 2; 3 una base de tres vectores r gidamente unidos a un punto material A del slido. Los tres vectores forman un triedro r o gido que gira igual que lo hace el material en A. En A se tiene un tensor gradiente de velocidades L: L=D+ y se cumplir que: a e_ = e Entonces, si se considera un tensor de segundo orden T denido en A con arreglo a la conguracin actual del slido (p.ej. el tensor de tensiones verdaderas), ste se podr o o e a poner como: T = T e et y su derivada con respecto al tiempo: 130

_ _ _ _ T = T e et + T e et + T e et t t _ = T e e + T e e + T e et t = T + T + T t En la relacin anterior, el primer sumando del segundo miembro se conoce con el o nombre de tasa de variacin corrotacional del tensor T, o tasa de Jaumann T. o En el caso de que T sea el tensor de tensiones verdaderas, la tasa de Jaumann es la parte de la variacin de las tensiones que corresponde a la respuesta del material. Los o otros dos sumandos proporcionan la variacin de T debida al giro de slido r o o gido del entorno del punto material A.
r r

12.10

Particin aditiva de la velocidad de deformacin o o

En muchos problemas de inters industrial las deformaciones elsticas o recuperables de e a los materiales son muy peque~as si se comparan con las deformaciones plsticas o no n a recuperables. Este hecho permite simplicar la descripcin de la deformacin del medio o o continuo. Se puede suponer que en cualquier instante del proceso de deformacin el material o tiene un estado natural de referencia. Dicho estado de referencia sera aquel al que tiende el material al descargarlo de todas las acciones exteriores que actan sobre l. As en un u e , instante del proceso de deformacin se puede aislar un elemento de volumen y descargarlo o de todas las tensiones que se ejercen sobre el mismo. La transformacin asociada a este o 1 proceso de descarga ser Fel , la inversa de la deformacin recuperable. a o La deformacin no recuperable ser la deformacin restante: o a o Fpl = F1 F el De la relacin anterior se sigue una descomposicin de la transformacin representada por o o o el tensor gradiente de deformacin en una parte recuperable y otra no recuperable: o F = Fel Fpl Por otro lado, el tensor gradiente de velocidades con respecto a la conguracin actual o es: _ L = F F1 _ _ = (Fel Fpl + Fel Fpl )(F1 F1 )
pl el

_ _ = Fel F1 + Fel Fpl F1 F1 el pl el = Lel + Fel Lpl F1 el 131

por analog con la denicin de L. a o El teorema de descomposicin polar puede tambin escribirse, separando la deformao e cin recuperable de la no recuperable, como: o

F = Vel Vpl R = Fel Fpl Y si la deformacin recuperable es peque~a, Vel se podr poner como: o n a Vel = I + "el a con deformaciones principales j"I j; j"II j; j"III j << 1. Entonces se tendr: Fel Vel = (1 + "elI ) nI nt + (1 + "elII ) nII nt + (1 + "elIII ) nIII nt I II III y L = Lel + (I + "el ) Lpl (I + "el )1 Lel + Lpl Tomando entonces la parte simtrica de L, se llega a: e _ _ _ " "el + "pl que es la descomposicin aditiva clsica del incremento de deformacin en una parte o a o elstica o recuperable y en otra parte no recuperable. Como se ha visto, esta descoma posicin slo es estrictamente vlida en un problema general cuando la deformacin recuo o a o perable es peque~a. n

132

Cap tulo 13

Procedimientos de clculo (III) a

13.1 Introduccin o

En este cap tulo se presenta la sistemtica general del MEF en problemas no lineales. a Al igual que en problemas lineales, la sistemtica del MEF permite pasar del problema a planteado sobre un medio continuo a un problema discreto energticamente equivalente. e Las ecuaciones que denen el problema discreto no son lineales. Para resolverlas existen en la prctica dos grandes familias de procedimientos, los procedimientos implcitos, a basados en el mtodo de Newton-Raphson y sus variantes, y los procedimientos explcitos, e que son una generalizacin de los algoritmos de integracin expl o o cita introducidos en el cap tulo 10 para problemas no estacionarios. En este cap tulo se presentan tambin las e l neas generales de ambas categor de procedimientos. as Estas dos familias de procedimientos denen las dos clases de programas de clculo a disponibles hoy en d en el mercado para abordar problemas no lineales. a

13.2

La \receta" del MEF en problemas no lineales

La sistemtica del MEF en problemas no lineales sigue las mismas etapas que la presentada a en el cap tulo 4 para problemas lineales. Al igual que entonces, el punto de partida es una forma dbil del problema de campo que se trata de resolver. e En el caso de problemas de Mecnica de Slidos, el punto de partida puede ser el a o principio de los trabajos virtuales, enunciado como: Dados los espacios S (desplazamientos admisibles) y V (variaciones o velocidades virtuales): S = f~ 2 H 1 () j u = d en Sd g u ~ V = f~ 2 H 1 () j ~ = 0 en Sd g v v 133

Encontrar u 2 S, de manera que: v t v ( u ; Dv ) = (b; ~ ) + [; ~ ]St 8~ 2 V v (13.1)

con la notacin utilizada en los cap o tulos anteriores. A partir de esta forma dbil de las condiciones de equilibrio, la receta del MEF es la e siguiente: 1. Dividir el dominio de clculo en subdominios o elementos: a =
[
i

\
i

i = ;

o 2. Dentro de cada elemento i , aproximar los campos incgnita que intervienen en la ecuacin integral 13.1: o ~ u = N ae ~ = N ae v Dv = ae ~ u = funcin no lineal de u = f (ae ) o donde N es la matriz de funciones de forma; ae y ae son vectores de coecientes incgnita, cuyo valor se hace coincidir normalmente con el valor de los desplazamieno tos o velocidades de los nodos; D es un vector que almacena las componentes del tensor velocidad de deformacin y es una matriz de funciones. o 3. Las integrales que aparecen en 13.1 se hacen sumando las contribuciones de los distintos elementos, es decir, 8~ 2 V: v 0 = ( u ; Dv ) (b; ~ ) [; ~ ]St = v t v
X
i

f( u ; Dv )i (b; ~ )i [; ~ ]Sti g u t v

4. La contribucin de cada elemento a la suma anterior es, sustituyendo las aproximao ciones del paso 2: t (u )i aei (bt N)i aei [tN]Sti aei Ntese que las tensiones no pueden ser linealizadas en funcin de los desplazao o mientos nodales aei , como ocurr en los problemas de elasticidad lineal. El primer a sumando corresponde a las fuerzas interiores del elemento y los dos ultimos suman dos, a las fuerzas exteriores aplicadas sobre el elemento. 5. La suma de las contribuciones de los distintos elementos se hace a travs del proceso e de ensamblaje.
t

134

El proceso de ensamblaje consiste en este caso en ir acumulando las contribuciones de cada elemento al vector global de fuerzas interiores Fint(a) y al vector global de fuerzas exteriores Fext (a). Como resultado del ensamblaje, se llega a la forma discreta de 13.1: at fFint (a) Fext(a)g at G(a) = 0 8 a

6. Como la relacin anterior ha de cumplirse 8 a, se llega nalmente al sistema global o de ecuaciones no lineales: G(a) Fint (a) Fext (a) = 0 el cual debe permitir determinar los coecientes incgnita a. o 7. De acuerdo con el paso 2, los coecientes a denen la solucin dentro de todos o los dominios elementales y, por agregacin, denen la solucin buscada en todo el o o dominio de clculo . a (13.2)

13.3

El mtodo de Newton y sus variantes e

Una de las formas ms populares de resolver el sistema de ecuaciones no lineales 13.2 a es el mtodo de Newton-Raphson o alguna de sus variantes, como el mtodo de Newton e e modicado o los procedimientos cuasi-Newton. El procedimiento de Newton es un mtodo iterativo relativamente sencillo de visualizar e en una dimensin, esto es, cuando el vector de movimientos nodales, a, tiene una sola o componente. Sea Fext una accin exterior constante. En este caso el problema representado por o 13.2 se reduce a encontrar el movimiento nodal a tal que: Fint (a) = Fext Es decir, el movimiento nodal a que hace que la fuerza de respuesta del sistema, o fuerza interna, sea igual a la accin exterior aplicada o fuerza externa. o Normalmente se opera de forma incremental, dividiendo la accin total Fext en suma o de incrementos: Fext =
X
n

(F n+1 F n ) =

X
n

F n

y buscando sucesivos puntos de equilibrio (an ; F n ) entre acciones externas y fuerzas internas.

135

F
6
n K1

Kn Fext F
n+1

Fn

: n r1

Fint (a)

an+1 1

n+1

Para pasar de un punto de equilibrio (an ; F n ) al siguiente, se procede iterativamente. Si K n representa la pendiente de la tangente a la curva F = Fint (a) en el punto (an ; F n ), se obtiene una primera aproximacin al punto (an+1 ; F n+1 ) mediante: o an+1 = an + 1 F n+1 F n Kn

La diferencia F n+1 Fint (an+1 ) se conoce con el nombre de residuo: 1

n r1 = F n+1 Fint(an+1 ) 1

La aproximacin al punto (an+1 ; F n+1 ) se rena entonces mediante: o an+1 = an+1 + 2 1

n r1 n K1

n donde K1 es la pendiente de la tangente a la curva F = Fint(a) en el punto dado por n+1 (a1 ; Fint (an+1 )). 1 n Se calcula ahora un nuevo residuo r2 = F n+1 Fint (an+1 ) y se sigue renando la 2 aproximacin hasta que el residuo se haga inferior a un umbral de tolerancia prestablecido. o Este esquema iterativo corresponde al procedimiento clsico de Newton o procedimiena to de Newton completo. En el procedimiento de Newton modicado se utiliza en todas las iteraciones la pendiente o rigidez inicial K n correspondiente al incremento. En un procedimiento cuasi Newton se emplea en cada iteracin una rigidez intermedia entre las o dos anteriores. Por ejemplo, la pendiente correspondiente a la secante que resulta de la iteracin anterior, es decir, el cociente: o

136

Fint(an+1 ) F n i n+1 ai an El procedimiento de Newton completo es, dentro de esta familia, el que tiene un orden de precisin ms elevado y una convergencia ms rpida, esto es, el que requiere de menos o a a a iteraciones para llegar a la solucin. Por contra, es tambin el que requiere ms trabajo o e a por iteracin ya que el clculo de las pendientes Kin precisa de un esfuerzo de clculo o a a signicativo. El formalismo bsico del procedimiento de Newton, esta vez ya en un problema mula tidimensional, es el siguiente. 1. En el incremento n, tras la iteracin i se tiene una aproximacin an al vector solucin o o i o de movimientos nodales. Sea cn la diferencia entre esta aproximacin y la solucin o o i exacta. G(an ) = rn i i n n G(ai + ci ) = 0 2. Desarrollando en serie de Taylor el trmino izquierdo de la ultima igualdad alrededor e de la solucin aproximada: o G(an ) + i @G(an ) n i ci + : : : = 0 @a

3. Si se desprecian todos los trminos del desarrollo menos los dos primeros se llega a: e @G(an ) n i ci ' G(an ) rn i i @a donde
@G(an ) i @a

(13.3)

es la matriz de rigidez tangente o matriz jacobiana de G.

4. La ecuacin 13.3 representa un sistema de ecuaciones lineales que proporciona la o correccin cn buscada. o i an = an + cn i+1 i i 5. El proceso dado por los puntos anteriores se repite hasta que en una iteracin j una o n norma del residuo rj es inferior a una tolerancia prestablecida.

13.4

Problemas no estacionarios

La sistemtica del MEF descrita en la seccin 13.2 es directamente aplicable a problemas no a o estacionarios o cuyos campos solucin son funcin del tiempo. En el caso de la Mecnica de o o a Slidos el principio de los trabajos virtuales puede utilizarse como forma dbil de partida o e si se incluye el trabajo de las fuerzas de inercia. Empleando la notacin de cap o tulos anteriores, dicha forma dbil se puede enunciar como: e 137

Dados los espacios S (desplazamientos admisibles) y V (variaciones o velocidades virtuales): ~ St = f~ (t) 2 H 1 (); t 2]0; T [ j u(t) = d(t) en Sd g u V = f~ 2 H 1 () j ~ = 0 en Sd g v v Encontrar u(t) 2 St , de manera que: u v ( u ; Dv ) = (b; ~ ) ( ; ~ ) + [; ~ ]St v t v y (u(0) u0 ; ~ ) = 0 v _ _ v (u(0) u0 ; ~ ) = 0 8~ 2 V v t 2 ]0; T [ (13.5) (13.6) 8~ 2 V v (13.4)

_ donde D es el tensor velocidad de deformacin y u0 y u0 representan desplazamientos o y velocidades iniciales, respectivamente. Ntese que las variaciones ~ 2 V no dependen del tiempo. o v La aplicacin de la \receta" del MEF descrita en la seccin 13.2 a la forma dbil 13.4 da o o e lugar a la semidiscretizacin espacial del problema y conduce, al igual que en los problemas o lineales1 a un sistema de ecuaciones diferenciales ordinarias cuya variable independiente es el tiempo. Este sistema de ecuaciones es ahora no lineal y puede escribirse como: _ _ Fext (a; a) Fint (a; a) = M a (13.7)

donde las funciones a(t) representan los movimientos nodales y M es la matriz de masas del sistema. El sistema 13.7 se integra, partiendo de las condiciones iniciales que derivan de 13.5 y 13.6, utilizando un procedimiento de diferencias nitas. Esta es la semidiscretizacin o temporal, por la que se obtienen las funciones a(t) en puntos discretos t1 , t2 : : : del eje de tiempos separados por un paso de integracin t: a(t1 ) a1 , a(t2 ) a2 , : : : o Como ya se coment en el cap o tulo 10 para problemas lineales, la integracin de 13.7 o partiendo de las condiciones iniciales se hace normalmente mediante procedimientos de la familia de Newmark. Estos procedimientos relacionan las aceleraciones, velocidades y movimientos nodales en el instante n con los correspondientes al instante n + 1 mediante las ecuaciones:
1

Ver cap tulo 10.

138

_ an+1 = an + t an +

t2 [(1 2) an + 2 an+1 ] 2

(13.8) (13.9)

_ an+1 = an + t [(1 ) an + an+1 ] _

donde y son dos parmetros que dependen del procedimiento concreto y los subndices a indican el instante de tiempo en que se particularizan las diferentes funciones vectoriales (p.ej. an+1 a(tn+1 )). A partir de las relaciones 13.8, 13.9 y de la particularizacin de 13.7 para el instante o tn+1 , se plantea un procedimiento en el que se progresa paso a paso en el tiempo a partir de las condiciones iniciales. Esquemticamente, este procedimiento de integracin paso a a o paso podr ser: a o 1. Calcular las aceleraciones iniciales a0 a partir de la relacin 13.7 y de los desplaza mientos y velocidades iniciales: _ _ a0 = M1 (Fext (a0 ; a0 ) Fint (a0 ; a0 ) 2. Inicializar n 0 _ 3. Poner an+1 y an+1 en funcin de an+1 , utilizando 13.8 y 13.9 o 4. Sustituir las expresiones an+1 y an+1 de la etapa anterior en la particularizacin de _ o 13.7 para el instante tn+1 : _ _ Fext(an+1 ; an+1 ) Fint (an+1 ; an+1 ) = M an+1 Se obtiene as un sistema algebraico de ecuaciones no lineales: _ _ a Fext (n+1 ; an ; an ; an ) Fint (n+1 ; an ; an ; an ) M an+1 = 0 a (13.11) (13.10)

que permite obtener los coecientes an+1 , ya que los valores de desplazamientos, velocidades y aceleraciones en el instante tn son conocidos 5. A partir de an+1 y de los valores ya conocidos del instante de tiempo anterior, an , _ an , an , calcular an+1 y an+1 utilizando 13.8 y 13.9 _ 6. Hacer n n + 1 7. Volver a 3 En la prctica se utilizan fundamentalmente dos procedimientos de la familia de Newa mark, el procedimiento de la aceleracin promedio ( = 0:5 y 2 = 0:5) y el procedimiento o de diferencias centrales ( = 0:5 y = 0)). El procedimiento de la aceleracin promedio o es un procedimiento impl cito, en el que cada ciclo de clculo, para pasar del instante de a 139

tiempo tn al instante tn+1 se resuelve un sistema no lineal de ecuaciones 13.10 del mismo tipo que el 13.2 utilizando alguna variante del mtodo de Newton: e G (n+1 ) Fext(n+1 ; an ; an ; an ) Fint (n+1 ; an ; an ; an ) M an+1 = 0 _ _ a a a

(13.12)

Por el contrario, el procedimiento de diferencias centrales, cuando se utiliza una matriz de masas M diagonal, da lugar a un procedimiento expl cito ya que el sistema de ecuaciones 13.10, tras sustituir 13.8 y 13.9, queda en la forma: _ _ an+1 = M1 fFext (an ; an ; an ) Fint (an ; an ; an )g (13.13)

cita" Es decir, puede obtenerse la aceleracin an+1 en en instante tn+1 de manera \expl o a partir de los valores conocidos de desplazamientos, velocidades y aceleraciones en el instante tn . Los procedimientos de integracin expl o cita en problemas no lineales se discuten con ms detalle en la seccin siguiente. Como ya se apunt en el cap a o o tulo 10, su principal inconveniente es que son slo condicionalmente estables, esto es, requieren que el paso o de integracin t sea inferior a un paso cr o tico tcrit. Sin embargo, en problemas no lineales, los esquemas de integracin expl o cita son mucho ms robustos que los esquemas a impl citos, ya que no sufren de los potenciales problemas de convergencia asociados a los procedimientos de la familia de Newton.

13.5
13.5.1

Procedimientos de integracin expl o cita

Ideas generales

En problemas mecnicos, la idea f a sica que hay detrs de los procedimientos de integracin a o expl cita es la de obtener la aceleracin de cada nodo dividiendo la resultante de fuerzas o sobre el nodo por la masa del mismo. Es decir, la idea f sica corresponde a la aplicacin o de la segunda ley de Newton a cada uno de los nodos del mallado. La matriz de masas diagonal se interpreta como que la masa del continuo est concentrada en los nodos. a El procedimiento consigue as el desacoplamiento de las ecuaciones globales del mo vimiento en ecuaciones locales. Este desacoplamiento slo resulta vlido si el paso de o a integracin t es lo sucientemente peque~o como para que dos nodos contiguos no teno n gan tiempo material de \comunicarse" durante el tiempo t. De aqu deriva la existencia de un l mite tcrit para el paso de integracin si se quiere que la integracin se mantenga o o estable. Este es el signicado f sico de la llamada condicin de Courant: t tcrit . o El tiempo que necesitan dos nodos para comunicarse est relacionado con el tiempo que a tarda la onda de tensin ms rpida en recorrer el espacio que los separa. En sistemas o a a no amortiguados, el paso cr tico de integracin resulta ser el m o nimo periodo propio de vibracin del sistema discreto2 dividido por . o
El conjunto de masas puntuales conectadas por muelles a que conduce la semidiscretizacin espacial o mediante el MEF.
2

140

La resultante de fuerzas sobre cada nodo es funcin, sobre todo, de la respuesta teno sional del material a la deformacin. En un procedimiento expl o cito las ecuaciones constitutivas, o de respuesta del material, deben estar formuladas de forma que al nal de cada paso de integracin t el tensor de tensiones se pueda actualizar de manera explcita, o esto es, no iterativa, en funcin de los incrementos de deformacin, de las velocidades de o o deformacin y de otras variables de estado. o Por otro lado, una de las principales ventajas de los procedimientos explcitos en el tratamiento de problemas no lineales es que, al desacoplar las ecuaciones globales del movimiento en ecuaciones locales correspondientes al movimiento de cada nodo, se pueden tratar los contactos tambin a nivel local. Es decir, no hay que buscar iterativamente en e cada paso de integracin una conguracin global estable de los contactos, como sucede en o o los procedimientos impl citos, sino que la interaccin a travs de las potenciales supercies o e de contacto se obtiene nodo a nodo, con independencia del estado del contacto en otras zonas de la malla. Las fuerzas de interaccin se evalan as tambin a nivel local y se o u e hacen intervenir en la ecuacin del movimiento de cada nodo de la misma manera que las o fuerzas internas. Este tratamiento local o desacoplado de los contactos proporciona una gran robustez a los procedimientos expl citos.

13.5.2

Flujo general

El ujo general en el ciclo de clculo que permite avanzar la solucin un paso de integracin a o o t puede esquematizarse de la manera siguiente: 1. Clculos nodales (bucle en nodos) a (a) Ecuacin del movimiento de los nodos: o an = M1 fFext Fintg (b) Integracin por diferencias centrales: o an+ 1 _
2

= an 1 + t an _
2 2

_ an+1 = an + t an+ 1 2. Clculos en los elementos (bucle en elementos) a (a) Calcular incremento de deformacin " en el elemento durante t o (b) Calcular las tensiones y variables de estado en el elemento en t + t a partir del incremento de deformacin " y del modelo de comportamiento del material o (c) Calcular y ensamblar la contribucin del elemento al vector de fuerzas internas o Fint en t + t, por integracin alrededor de los nodos a partir de las tensiones o 3. t t + t y volver a 1

141

13.5.3

Ventajas e inconvenientes

Las principales ventajas del ciclo de clculo de integracin expl a o cita son su sencillez, que lo hace fcilmente paralelizable, y su gran robustez en problemas no lineales, ya que a desaparecen los problemas asociados a la falta de convergencia que tanto penalizan a veces los procedimientos impl citos de la familia de Newton. El inconveniente ms importante, y casi unico, es que, como se ha se~alado, el ciclo a n de clculo es slo condicionalmente estable, requirindose un paso de integracin t rea o e o lativamente peque~o para progresar en el tiempo. Para que la integracin sea estable se n o requiere que: t tcrit y resulta que: 2 q ( 1 + 2 ) tcrit = !max donde, !max es la mxima frecuencia natural de vibracin del modelo de elementos nitos a o o en el instante t en el que se calcula el paso cr tico de integracin3 y es la fraccin del o amortiguamiento cr tico a esa frecuencia. Esta es una restriccin muy fuerte, ya que implica que el paso de integracin t debe ser o o inferior en cada instante al tiempo que tarda la onda de tensin ms rpida en atravesar el o a a elemento nito ms peque~o. Ntese que un aumento de la rigidez o del amortiguamiento a n o disminuyen el paso cr tico, mientras que un aumento de la densidad o del tama~o de los n elementos, lo incrementan.

13.5.4

Medidas de calidad de la solucin o

La restriccin que supone el que el paso de integracin deba ser inferior al cr o o tico hace que la calidad numrica de la solucin est asegurada siempre que se calcule bien ese e o e paso cr tico tcrit . Un paso de integracin inferior al cr o tico supone que el esquema de integracin es capaz de capturar todas las frecuencias posibles en la respuesta dinmica o a del sistema discreto. Si el paso cr tico se calcula mal y en algn momento se emplea un paso de integracin u o t superior a tcrit, el procedimiento se hace automticamente inestable. La inestabilia dad tiene el efecto de generar energ no f a sica dentro del mallado. Cuando esto ocurre, normalmente la solucin \explota": se obtienen grandes distorsiones en la malla y resultao dos claramente inaceptables. Sin embargo, otras veces el exceso de energ que se genera a se convierte localmente en energ de deformacin plstica, la deformacin plstica reblana o a o a dece la respuesta del material y el paso cr tico vuelve a ser superior al paso de integracin o que se est utilizando. En estos casos a simple vista el analista no distingue nada y podra a dar por buena su solucin. o
3 En un problema no lineal las frecuencias propias del modelo pueden cambiar a lo largo del tiempo del problema.

142

El procedimiento de integracin expl o cita no conserva la energ La calidad de la a. solucin obtenida puede medirse por lo bien o mal que se haya conservado la energa a o lo largo del problema. El analista debe vigilar siempre la contabilidad de la energa y vericar que se cumple: Edeformacin + Ecintica + Edisipacin : : : Tfuerzas exteriores = constante o e o donde E representa energ y T , trabajo. a Un aumento en el primer miembro de la ecuacin anterior indica siempre que se ha o producido inestabilidad, por lo que no ser admisible a menos que el analista determine a que la inestabilidad tiene carcter local y que no inuye en los parmetros que le interesan a a de la solucin del problema. Sin embargo, es normal, sobre todo en problemas con muchos o ciclos de integracin, que se produzca una ligera disminucin del primer miembro de la o o ecuacin, a causa de la disipacin numrica asociada al procedimiento de integracin. o o e o

13.5.5

Campo de aplicaciones industriales

El campo natural de aplicacin de los procedimientos de integracin expl o o cita se encuentra en aquellos problemas en los que necesariamente el paso de integracin debe ser peque~o o n para capturar correctamente la respuesta. Es el caso de los problemas dominados por las altas frecuencias, como los problemas de impacto. En esta clase de problemas la restriccin o de que el paso de integracin deba ser inferior al paso cr o tico no es un inconveniente demasiado grande porque, en cualquier caso, la naturaleza del problema exigira ir a pasos de integracin peque~os para no ltrar las frecuencias en donde se concentra la respuesta. o n Los procedimientos de integracin expl o cita tambin tienen aplicacin, por su robuse o tez, en problemas con una no linealidad muy acusada. Es el caso de los problemas con condiciones de contacto muy complicadas, problemas con inestabilidades locales acusadas, reblandecimiento del material o fuerte dependencia de sus propiedades con respecto a la temperatura. Estas caracter sticas conuyen muchas veces en los problemas de conformado de metales.

13.5.6

Tcnicas de aceleracin del clculo e o a

Las tcnicas de aceleracin del clculo se hacen necesarias cuando se trata de resolver un e o a problema cuasiesttico o dinmico \lento" con procedimientos de integracin expl a a o cita. Si no se emplean esta clase de tcnicas de aceleracin, se obtienen normalmente tiempos de e o ejecucin demasiado largos para ser viables. o Ntese que un problema puede considerarse cuasiesttico cuando los tiempos de aplio a cacin de la carga son superiores a unas dos veces el periodo fundamental Tf del sistema o discreto. Esto siginica que, aproximadamente, se obtiene la misma solucin aplicando la o carga en un tiempo innito (problema esttico) que en un tiempo 2 Tf . a De este modo, una primera tcnica para reducir el tiempo de clculo en el ordenador e a correspondiente a un problema cuasiesttico consiste en aplicar la carga ms rpidamente, a a a hasta aplicarla en un tiempo de entre Tf y 2 Tf . El analista debe comprobar luego que

143

la energ cintica del problema se mantiene baja (inferior al 5% de la total) cuando las a e deformaciones empiezan a ser signicativas. Sin embargo, la tcnica anterior no es aplicable cuando las propiedades del material e dependen de la velocidad de deformacin (viscoelasticidad, viscoplasticidad). En este caso, o acelerando la aplicacin de la carga se falsear la respuesta del material. Cuando existe o a dependencia de la respuesta del material de la velocidad de deformacin debe recurrirse, o para reducir los tiempos de clculo en el ordenador, a tcnicas de escalado de masa. a e Ntese que el paso cr o tico de integracin es proporcional a la ra cuadrada de la densidad o z y que en problemas cuasiestticos la solucin es independiente de la densidad. Entonces, a o para reducir el nmero de ciclos de clculo necesarios para llegar al tiempo nal del u a problema, puede aumentarse articialmente la densidad o la masa del sistema de manera que el periodo fundamental Tf del sistema discreto llegue a ser 0.5 o 1 veces el tiempo de aplicacin de la carga o el tiempo caracter o stico del problema.

144

Bibliograf recomendada a

1. T.J.R. Hughes. The Finite Element Method. Prentice Hall. 1987. 2. E. O~ate. Clculo de Estructuras por el Mtodo de los Elementos Finitos. CIMNE. n a e 1992. 3. O.C. Zienkiewicz y R.L.Taylor. El Mtodo de los Elementos Finitos. McGraw-Hill. e 1994.

Bibliograf complemetaria a

1. M.A. Criseld. Finite Elements and Solution Procedures for Structural Analysis. Pineridge Press. 1986. 2. M.A. Criseld. Non-linear Finite Element Analysis of Solids and Structures. John Wiley. 1991. 3. Hibbitt, Karlsson & Sorensen, Inc. ABAQUS Theory Manual v.5.6. 1996. 4. T. Mura y T. Koya. Variational Methods in Mechanics. Oxford University Press. 1992. 5. NAFEMS. A Finite Element Primer. National Agency for Finite Elements and Standards, Reino Unido. 1987. 6. J.T. Oden y G.F. Carey. Finite Elements - Mathematical Aspects. Prentice-Hall. 1983. 7. J.N. Reddy. Applied Functional Analysis and Variational Methods in Engineering. McGraw-Hill. 1986. 8. I.M. Smith. Programming the Finite Element Method. John Wiley. 1982. 9. G. Strang y G.J. Fix. An Analysis of the Finite Element Method. Prentice-Hall. 1973. 10. K. Washizu. Variational Methods in Elasticity and Plasticity. Segunda edicin. o Pergamon Press. 1974. 11. O.C. Zienkiewicz y K. Morgan. Finite Elements and Approximation. John Wiley. 1983. 145