El método de los elementos finitos adaptable. Estimaciones de error

Revista internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 9, 4, 417-437( 1993) EL MÉTODO DE LOS ELEMENTOS FINITOS ADAPTABLE ESTIMACIONES DE ERROR CASTOR GONZALEZ G.* RAUL MANZANILLA** Y RAUL FEIJOO*** *Universidad Centrooccidental Lisandro Alvarado, Departamento de Procesos Agroindustriales, Apdo. 400, Barquisimeto, Edo. Lara, Barquisimeto, Venezuela **Intevep, S.A. Análisis Matemático y Estadístico, Apdo. 76343, Caracas 1070-A, Los Teques, Venezuela ***L.N. C. C. Rua Lauro Muller, 455-Caiza Postal 56018-22290, Río de Janeiro, Brasil RESUMEN Actualmente el desarrollo de software para el análisis de sistemas de ingeniería a través del método de los elementos finitos está notablemente influenciado por la solución al problema "EL METODO DE LOS ELEMENTOS FINITOS ADAPTABLE". Una estimación confiable del error es el punto central en el planteamiento de la solución a este problema. A partir de 1987 con la publicación de Zienkiewicz y Zhul, numerosos trabajos han sido publicados tratando de proponer estimadores de error y estrategias de refinamiento que sean computacionalmente aceptables para resolver este problema. En esta publicación se presenta una solución al método de los elementos finitos adaptable, basada en un nuevo tipo de estimador de error, obtenido de los trabajos de Johnson y Hansbo Este nuevo estimador constituye una generalización al propuesto por Zienkiewicz y Zhu y la idea central es comparar ambos estimadores. SUMMARY At the present time, the development of software for analysis of engineering systems, through the finite elements method, is notably influenced by the solution of the problem of "ADAPTABLE FINITE ELEMENT METHOD". A trustly estimation of the error is the central pointin getting a solution to this problem. After the publication of Zienkiewicz and Recibido: Abril 1993 OUniversitat Politicnica de Catalunya (España) ISSN 0213-1315 C. GONZALEZ G;, R. MANZANILLA Y R. FEIJOO Zhu', in 1987, many publication have been made trying to propose estimators of the error and computational strategies of refinement that are acceptable to resolve this prob1em.-In this work is presented a solution to the problem of adaptable finite element method, which is based on a new type of error estimator, obtained from the work of Johnson and Hansbo' and that is a generalization of the one propose by Zienkiewicz and Zhu. Aditionally a comparison among both is presented. INTRODUCCION El método de los elementos finitos adaptable, además de dar respuesta numérica a un problema específico pretende resolver uno de los aspectos más dificultosos del modelaje numérico, como lo es la validación de los resultados y en consecuencia la confiabilidad de la solución. Dentro de este panorama, la evaluación del error cometido en la solución aproximada juega un papel fundamental. Las estimaciones de error clásicas, tanto por diferencias finitas como para elementos finitos, se basan en cotas superiores de error que tienen convergencia asintótica cuando el tamaño de malla tiende a cero. Desafortunadamente, estas estimaciones dan poca información sobre la magnitud del error para un tamaño de malla fijo. Los primeros trabajos para establecer un error a-posteriori (obtenido de la solución aproximada), fueron publicados por Babuska y coa~tores~"*~~~~"~~*~~~. LOS métodos presentados en estos trabajos conducen a buenos indicadores de error pero son computacionalmente costosos y no son generales. El estimador de error propuesto por Zienkiewicz y Zhu (Z2), puede ser considerado como una técnica de interpolación, un tanto más general por se independiente de la formulación del problema, aunque menos precisa. Este estimador fué obtenido en forma intuitiva y funciona solamente cuando la derivada de la solución aproximada, presenta saltos en la interface de los elementos. El objetivo de este trabajo es implantar un nuevo tipo de estimador basado en los trabajos de Johnson y Hansbo2(~&H), el cual constituye una generalización del método Z2 y es aplicable a una gama más extensa de problemas. Para obtener resultados numéricos, que permitan la comparación de los estimadores, se resuelve el método de los elementos finitos adaptable aplicado al problema de Poisson para campos escalares. Un enfoque para encontrar soluciones aproximadas es el siguiente: dado un problema P cuya solución exacta u pertenece a un espacio de funciones H con norma 1 1, construir un subespacio Hh de H, la solución aproximada uh será la función de Hh que satisface: A grandes rasgos el método de los elementos finitos es una técnica sistemática para encontrar el espacio de aproximación Hh. Para construir este espacio, el método define una malla Th, con tamaño de malla h, la cual constituye una partición de la región de definición del problema, sobre esta malla se define un conjunto de funciones de interpolación que constituyen una base del espacio Hh. La Figura la muestra una representación gráfica de este problema. a)SoluciÓn de elementos finitos b) Método de los elementos finitos adaptable Figura 1. Representación gráfica del método de los elementos finitos y el método de los elementos finitos adaptable. Por otro lado, el método de los elementos finitos adaptable puede ser planteado en los siguientes términos: Sea P un problema que se pretende resolver por elementos finitos y que tiene solución exacta u perteneciente al espacio H, con norma 1 1. Denotemos por {Pi} a la familia de problemas de elementos finitos asociada al problema P, que consisten en aproximaciones polinomiales a trozos sobre una familia de mallas {T~,) con tamaño de malla local hg. Cada malla junto con los polinomios asociados a ella, define un espacio Hh, C H, teniendo así una familia de espacios de elementos finitos {Hhi). Sea {uh,} la familia de soluciones aproximadas correspondiente a la familia de mallas {rhi ). Finalmente sea 77 una tolerancia prefijada. Basados en estos datos, el objetivo consiste en implantar un algoritmo que conduzca a una malla 7hopt que define Hhopt, de tal manera que la solución correspondiente de elementos finitos uhopt cumpla con: y el trabajo para calcular rhopt sea mínimo, esto último define la eficiencia del método. Una medida del trabajo computacional estaría representada por el número de grados de libertad. La Figura lb muestra una representación gráfica de este problema. En otras palabras, el objetivo del método de los elementos finitos adaptable es obtener un algoritmo para resolver en forma aproximada el siguiente problema de optimización lineal: encontrar una malla óptima rhopt con la restricción )u - uhOpt) < 7. C. GONZALEZ G., R. MANZANLLA Y R. FEIJOO Este problema es difícil de resolver ya que la solución exacta u no está disponible y no existe una relación directa entra h y Iu - uhl. Teniendo presente que no se conoce el error 1% - uhl, el punto central del método de los elementos finitos adaptable es disponer de estimaciones de error confiables para poder implantar el algoritmo de refinamiento. Más adelante enfocaremos nuestra atención hacia la búsqueda de estimaciones de error. ESQUEMA DE REFINAMIENTO Los procedimientos para el refinamiento de malla parten de una discretización muy somera, que se enriquece sucesivamente hasta alcanzar una solución con una tolerancia predefinida. Cuando se dispone de una estimación del error, los métodos de refinamiento de malla adaptable se basan fundamentalmente en lograr una distribución uniforme del error. Es decir, evaluar la estimación del error a nivel de cada uno de los elementos de la malla y refinar según su magnitud. Se pueden plantear otros esquemas donde no es necesario conocer el error23, por ejemplo, para el caso donde el problema pueda ser modelado como la minimización de un funcional, se puede implantar un proceso de refinamiento en la dirección del gradiente. Esto consiste en dividir cada uno de los elementos y evaluar la variación del funcional en esa dirección y posteriormente refinar según la magnitud de estas varia~iones~~. Esta técnica es computacionalmente costosa. El algoritmo de refinamiento propuesto sigue un esquema de distribución uniforme del error. Para encontrar una malla óptima rhopt, en la cual se satisface que le1 = Iu - uhl 5 q, y donde el trabajo computacional para conseguir ~h~~t sea cercano al mínimo, se siguen los siguientes pasos: P1) Dada la región R de definición del problema, construir una malla rhUiejo lo más "gruesa" posible, pero que de una aproximación aceptable de R, esta malla define un hiUiej, y asocia polinomios de grado menor o igual que p a cada elemento. Sea N, y N, el número de nodos y el número de elementos, respectivamente. Calcular la solución de elementos finitos uh en esta malla. P2) Dada la tolerancia porcentual q (en el sentido le1 / 1~1): P2.1) Calcular el error permisible a nivel de cada elemento: P2.2) En cada elemento calcular: Si t; es menor que uno se desrefina y si (; es mayor que uno se refina la malla. P3) Refinar o desrefinar en cada elemento según: p representa el orden de convergencia de la solución aproximada. Así tendremos construida una nueva malla T~,,,~,. P4) Calcular el error porcentual actual qactua, en la nueva malla: Si qactual 5 q parar, sino repetir el proceso. Se pretende que este algoritmo sea confiable y eficiente en el sentido de que el control del error quede garantizado y que el trabajo computacional sea lo más pequeño posible. En la práctica le1 no está disponible y por esta razón debemos utilizar estimadores del error. EL PROBLEMA DE POISSON: TRANSFERENCIA DE CALOR El problema de Poisson para campos escalares se modela de la siguiente manera. Dada una región R del plano con contorno I' = ru U I', como se muestra en la Figura 2. Figura 2. Región Q. Sea L2(R) el espacio de las funciones cuadrado integrables dotado con la norma y producto interno C. GONZALEZ G., R. MANZANILLA Y R. FEIJOO y sea H1(52) el espacio de Sobolev de orden uno con la norma El problema de conducción de calor para régimen estable en dos dimensiones se modela matemáticamente de la siguiente manera: encontrar una función u que satisface donde u0 E L2(ru) es la temperatura prescrita en I', kx,k, son las conductividades térmicas. n es la normal a I' f son las fuentes de calor internas en 52. El flujo de calor asociado a este problema está dado por Al problema anterior se le puede asociar una formulación variacional. Para ello es necesario considerar los siguientes espacios afines de funciones y los operadores funcionales: La formulación variacional, basada en estas definiciones, se expresa en los siguientes términos: encontrar u E H tal que: EL MÉTODO DE LOS ELEMENTOS FINITOS ADAPTABLE Adicionalmente, la norma energía asociada a este problema está definida por: Una vez formulado el problema variacional (M), para implantar la solución aproximada se necesita construir un espacio de dimensión finita Hh al cual pertenece la solución aproximada uh. Para ello, en el método de los elementos finitos se introduce una malla rh que consiste en una partición de 52 en elementos finitos K de diámetro hK. Adicionalmente sea Pi(K) el espacio de los polinomios de grado uno sobre K E rh, y sea Vh el espacio definido por: Los correspondientes espacios de elementos finitos Hh y HhO quedan definidos por: De esta forma se suele formular el problema de elementos finitos siguiente: Encontrar uh E Hh tal que Aplicaremos el método de los elementos finitos adaptable, descrito anteriormente, a este problema particular de Poisson, utilizando la norma energía l.lB para la evaluación del error. Definido el problema que nos ocupa y teniendo presente-que la so1ució.n exacta u no es conocida, y por lo tanto el error le1 = lu - uhl no está disponible, a continuación se plantean estimadores de error que permitan implantar el método de los elementos finitos adaptable. De aquí en adelante se considera a 0 como un dominio con contorno poligonal. ESTIMADORES DE ERROR Los primeros trabajos para establecer un error a posteriori; es decir, obtenido de la solución aproximada uh, se formularon en base a criterios de superconvergencia. Babuska y coautores se basan en esta orientación. Posteriormente se origina el método Z2, considerado como una técnica de proyección. Más recientemente Johnson y Hansbo proponen una estimación de error que da origen al método que denominaremos J & H, este Último método pretende englobar los de Babuska y Z2. Otro tipo de estimador de error es el propuesto por Beltran y Alarcon*' basado en la diferencia de ciertos funcionales que definen formulaciones diferentes de un mismo problema. A continuación daremos una reseña de los trabajos de Babuska y Kelly en lo que respecta a estimadores de error, y posteriormente plantearemos los estimadores de Z2 y J&H. C. GONZ~LEZ G., R. MANZANILLA Y R. FEIJOO Para la deducción de los estimadores denotaremos Sh = {S} la colección de lados S de los elementos K E rh que pertenecen a I', para cada S E Sh, el cual es común a dos elementos finitos, definimos ns como la normal unitaria a s. En los trabajos de Bab~ska~~'-~~~, kellyl1J2 y coautores se pone de manifiesto que el error en la solución de elementos finitos está relacionado por una parte con el residuo: y por otro lado con el salto en el flujo normal asociado a la frontera inter-elemento l)nz al), de finido por: T ((n$.all= lim n$.alim ns .a x4x+ x4xPara ver esto, tomemos la formulación variacional (23) que establece: Si las funciones que definen el espacio de elementos finitos Hh son los polinomios Pi(K), tendremos discontinuidad en la derivada de uh sobre la interface de los elementos. Integrando por partes (26) a nivel de cada elemento, obtenemos: donde: O si s E I', salto del flujo sis E R 6, + nzah si s E r, Esta última expresión pone de manifiesto la aparición del residuo de la ecuación diferencial y el salto del flujo. Adicionalmente se puede establecer, Kelly"*12, que el error satisface el problema auxiliar: Donde R incluye tanto el residuo como el salto del flujo. Ésta es una forma de evaluar el error, pero evidentemente es poco práctica y computacionalmente costosa. Estimador de error de Zienkiewicz-Zhu (Z2) En general los estimadores de error de Babuska requieren de la solución de un problema auxiliar y su implantación computacional es costosa. Si se trabaja con elementos de Pl(K), y en el caso particular de que las fuentes de calor f sean nulas, de acuerdo con los operadores que modelan el problema de Poisson, el residuo STDSuh - f es nulo y el salto pasa a jugar un rol preponderante en el cálculo del error, esto puede generalizarse a otros problemas. Zienkiewicz y Zhu introdujeron un estimador que toma en cuenta este hecho, adicionalmente es simple y fácil de incorporar en los códigos existentes de elementos finitos. Sin embargo las bases sobre las cuales el estimador fué derivado son eurísticas. A continuación daremos los lineamientos generalesdel estimador Z2. La definición de la norma energía dada por (19) puede ser escrita como: Siendo uh y Ch la solución y el flujo aproximados por elementos finitos, la norma del error e = u - uh será: La solución de elementos finitos uh puede ser escrita como T Uh = NhUh (33) donde Nh es el vector de las funciones de interpolación relativas a rh y Uh es el vector que contiene los valores nodales de uh. En consecuencia la aproximación para el flujo es : El estimador de error Z2 se basa en sustituir el flujo exacto a por una función pos-procesada a* obtenida de ah. Una alternativa es el flujo suavizado: donde los valores nodales E; se obtienen a través de la proyección de ah es el espacio de elementos finitos Hh es decir, ~r* debe satisfacer: donde P es un vector de funciones peso que constituyen una bas de Hh; en particular P puede ser Nh. En la práctica los valores nodales puede obtenerse como el promedio de los valores de ah de los elementos que concurren en los nodos. La estimación del error toma la forma: 432 C. GONZ~LEZ G., R. MANZANILLA Y R. FEIJOO Los resultados numéricos para el estimador z2, = lelZz, en (37) con tomado como el promedio de los ah que concurren en un nodo, y el estimador (él = lelJkH dado por la ecuación (57) pueden observarse en la Tabla 1. Esta tabla fué obtenida para un refinamiento uniforme. Tabla I. Resultados numéricos del ejemplo 1. ~~~~~~,,,,~,=1.255656819. La ecuación (58) para este ejemplo particular, sin singularidad toma la forma: -112 lelB 5 CNn (63) En la Figura 5 puede observarse que este orden de convergencia es seguido tanto por el estimador Z2 como por el de J & H. EJEMPLO 2. - + - ll e ll Z&Z --*- -llellJ&H -1 40.5 0.7 0,9 1.1 1.3 1.5 1.7 1,9 2.1 2.3 2.5 LOG N Figura 5. Tasa de convergencia de los estimadores Z2 y J & H. EJEMPLO 2. Sea la región SZ = [(T, 8), O 5 T I. R, O I. 0 5 3/2n] mostrada en la Figura 6 y el siguiente problema de valor de contorno: Figura 6. Región circular con singularidad. En este caso la matriz S es la misma dada en (10) y D, = TD. La solución exacta con Icx = Ic, es: La norma de u toma el valor lulB = @R'/~. Los resultados numéricos para un refinamiento uniforme con R = 10 se muestran en la Tabla 11. Tabla 11. Resultados numéricos del ejemplo 2. IIuIIB~~~~~~=~ .909318031. En este caso tenemos una singularidad X = 213, la ecuación (58) será: En la Figura 5 se observa nuevamente que el orden de convergencia establecido por la teoría es seguido por los estimadores en estudio, para un refinamiento uniforme. Seguidamente aplicaremos el método de los elementos finitos adaptable al problema del ejemplo 2. En la Figura 7 se muestra la secuencia de mallas obtenida con el estimador Z2 y el estimador J & H. Podemos observar que el orden de convergencia con el refinamiento adaptable tiende al máximo de 112. El valor de la tolerancia es 7 = 0.2. EJEMPLO 3. Sea una chimenea cuadrada como muestra la Figura 8a. Se quiere conocer la distribución de temperatura en la pared. Por simetría se estudiará solamente un cuarto de la chimenea como se muestra en la Figura 8b. C. GONZALEZ G., R. MANZANILLA Y R. FEIJOO LOG N Figura 7. Secuencia de mallas adaptables con los estimadores de Z2 y J & H. Orden de convergencia. Región circular. Figura 8. Distribiición de temperatura en una chimenea. La ecuación diferencial es la misma del ejemplo 1 y las condiciones de contorno se dan en la Figura 8. En la Figura 9 se muestran los resultados; esta figura corresponde al método de los elementos finitos adaptable con 7 = 0.1. Se puede apreciar que el estimador de J & H concuerda mejor con lo establecido por la teoría de aproximación. Figura 9. Secuencia de mallas adaptables con los estimadores de Z2 y J & H. Orden de convergencia. Chimenea. C. GONZÁLEZ G.. R. MANZANILLA Y R. FEIJOO CONCLUSIONES Tanto las estimaciones de Zienkiewicz-Zhu como la de Johnson-Hansbo son adecuadas para implantar el Método de los Elementos Finitos Adaptable (M.E.F.A.) en el caso de problemas elípticos. Como puede observarse en las tablas 1, 11 y en la Figura 5, las estimaciones de error de Z2 y J & H proporcionan valores cercanos a los exactos, el estimador de J & H constituye una cota superior del error, mientras que el estimador de Z2 no lo es. En cuanto al orden de convergencia en problemas con singularidad, que constituye otra forma de comparar los estimadores, el estimador de J & H presenta un comportamiento similar al establecido por la teoría de aproximación; es decir, en el refinamiento adaptable el orden de convergencia tiende al correspondiente a un problema sin singularidad. El estimador de Z2 presenta este mismo comportamiento en el ejemplo 2 pero difiere en el ejemplo 3. En lineas generales se puede decir que el estimador de Z2 es más eficiente (más fácil de implantar y obtiene el error deseado con menos trabajo computacional), pero menos confiable que el basado en la propuesta de J & H. Adicionalmente este último tiene mayor generalidad, cubriendo una gama de problemas más extensa. REFERENCIAS 1. O.C. Zienkiewicz y J.Z. Zhu, "A simpleerror estimator and adaptive procedure for practica1 engineering analysis", Int. J. Num. Meth. Eng. Vol. 24, pp. 337-357, (1987). 2. C. Johnson y P. Hansbo, "Adaptive finite element methods in computational mechanics", Preprint No 1992-041ISSN347-2809. The University of Goteborg. (1992). 3. 1. Babuska y A. Miller, "The post-processing approach in the finite element method. Part 1: Calculation of displacements, stresses and other higher derivatives of the displacements", Int. J. Num. Meth. Eng., Vol. 20, pp. 1085-1109, (1984). 4. 1. Babuska y A. Miller, "The post-processing approach in the finite element method. Part 11: The calculation of stress intcnsity factor", Int. J. Num. Meth. Eng., Vol. 20, pp. 1111-1129, (1984). 5. 1. Babuska y A. Miller, "The post-processing approach in the finite element method. Part 111: A posteriori error estimates and adaptive mesh selection", Int. J. Num. Meth. Eng., Vol. 20, pp. 2311-2324, (1984). 6. 1. Babuska y W.C. Rheinboldt, "A posteriori error estimates for the finite element method", Int. J. Num. Meth. Eng., Vol. 12, pp. 1597-1615, (1978). 7. Y. Liu y H. Elmaraghy, "Assessment of discretized erros and adaptive refinement with quadrilateral finite elements" , Int. J. Num. Meth. Eng., Vol. 33, pp. 781-798, (1992). 8. M. Ainsworth, J.Z. Zhu, A. Craig y O.C. Zienkiewicz, "Analysis of the Zienkiewicz-Zhu a posteriori error estimator in the finite element method", Int. J. Num. Meth. Eng., Vol. 28, pp. 2161-2174, (1989). 9. E. Rank y O.C. Zienkiewicz, "A simple error estimator in the finite element method", Common. Appl. Num. Methods, Vol. 3, pp. 243-249, (1987). 10. O.C. Zienkiewicz y J.Z. Zhu, "The three R's of engineering analysis and error estimation and adaptivity", Comp. Methods Appl. Mech. and Eng., Vol. 82, pp. 95-113, (1990). 11. D.W. Kelly y J.P. de S.R. Gago, O.C. Zienkiewicz y 1. Babuska, "A posteriori error and adaptive procedure in the finite element method. Part 1: Error analysis. Part 11: Adaptive mesh refinement", Int. J. Num. Meth. Eng., Vol. 19, pp. 1593-1619, (1983). 12. D. Kelly, "The self-equilibration of residuals and complementary a posteriori error estimates in the finite element method", Int. J. Num. Meth. Eng., Vol. 20, pp. 1491-1506, (1984). 13. 1. Babuska, O.C. Zienkiewicz, J.P. de S.R. Gago y E.R. de A. Oliveira (eds.), "Accuracy Estimates and Adaptive Refinement in Finite Element Computations", Wiley, New York, (1986). 14. 1. Babuska y W.C. Rheinboldt, "Error estimates for adaptive finite element computations", SIAM J. Numer. Anal., Vol. 15, pp. 736-754, (1978). 15. 1. Babuska y A. Miller, "A posteriori error estimates and adaptive techniques for the finite element method", Tech. Note BN-968, Inst. for Phys. Sci. and Techn., Univ. of Maeyland, (1981). 16. F.J. Beltran y E. Alarcon, "Estimación del error de discretización con principios variacionales multicampo: 1. Elasticidad", Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingenieráa, Vol. 8, pp. 25-43, (1992). 17. J.Z. Zhu y O.C. Zienkiewicz, "Superconvergence recovery technique and a posteriori error estimators'", Int. J. Num. Meth. Eng., Vol. 30, pp. 1321-1339, (1990). 18. R. Arcangeli y J.L. Gout, "Sur 17evaluation de l'erreur d'interpolation de Lagrange dans un ouvert de Rfl", R.A.I.R.O. Analyse Numerique, Vol. 10, pp. 5-27, (1976). 19. G.G. Castor, "Introducción matemática al método de los elementos finitos7', Tesis de pos-grado, U.C.L.A. (1986). 20. E.A. Fancello, A.C. Salgado y R.A. Feijóo, "Aranha: generador de malhas 2D para elementos finitos triangulares de 3 e 6 nos", Publicación interna LNCC/CNPq, Vol. 21, (1990). 21. R.E.d Ewing, "A posteriori error estimation", Comp. Methods Appl. Mech. and Eng., Vol. 82, pp. 59-72, (1990). 22. B.A. Szabó, "The use of a priori estimates in engineering computations", Comp. Methods Appl. Mech. and Eng., Vol. 82, pp. 139-154, (1990). 23. G.G. Castor, "El método de los elementos finitos adaptable", Tesis de doctorado, U.C.V, Por presentar. 24. G. Carey y T. Oden, "Finite elements: A second course. Volumen 11", Prentice Hall, Inc. Englewood Cliffs, New Jersey, (1983).