El objetivo de este trabajo es definir un método semiexplícito para estimar el error en cantidades de interés derivado de la resolución por elementos finitos de un problema de elasticidad lineal. En este proceso, se identifican dos etapas: una estimación implícita del error del problema adjunto y otra, explícita, asociada al problema directo (primal). La parte implícita de la estimación requiere dos fases, cada una de las cuales consiste en hacer proyecciones sobre espacios de funciones burbuja. Estas proyecciones sobre espacios burbuja comportan operaciones locales y, por tanto, son de bajo coste. Las dos fases se encargan de estimar el error interior en los elementos y la contribución de las aristas (asociada a los saltos de tensiones normales). El carácter implícito de esta parte de la estimación se aplica a la resolución de los problemas locales en espacios de dimensión pequeña. Se analiza también la posibilidad de escoger espacios de funciones burbuja locales de una sola dimensión (fijando la forma de la burbuja y determinando únicamente el coeficiente escalar que la multiplica), lo cual, en la práctica, convierte esta fase en explícita. La parte explícita de la estimación inyecta en el residuo débil primal la aproximación del error adjunto estimada en la primera fase. Esta estrategia es similar a la empleada en el método DWR, pero utilizando una formulación débil del residuo e inyectando un error dual en vez de utilizar una solución posprocesada. La metodología propuesta se valida con su aplicación a un ejemplo numérico.
We aim at defining a semi-explicit approach to estimate the error in quantities of interest associated with the Finite Element solution of a linear elasticity problem. The advocated procedure is split in two parts, an implicit error estimate for the adjoint problem and an explicit estimate assessing the error in the direct (primal) problem. The implicit part of the estimate (on the adjoint problem) embraces two phases, each consisting in projecting the error on “bubble” functional spaces. The projections are low-cost operations due to their local nature. The two phases account the error in the interior of the elements and the contribution of the elementary edges (associated with the tractions jump). The implicit character of this part is provided by the solution of the local problems in low-dimensional functional spaces. We also analyze the particular case of selecting one-dimensional functional spaces (setting the shape of the bubble and computing a scalar coefficient), which, in practice, make this part of the process explicit. The explicit part of the estimate consists in injecting in the weak primal residual the approximation of the adjoint error obtained in the first phase. This approach is similar to the DWR method but using a weak formulation of the residual and injecting an implicitly estimated adjoint error rather than a post-processed solution. The proposed methodology is validated in a numerical example
Cuando se resuelve un problema con el método de los elementos finitos, la solución numérica es una aproximación de la solución exacta, que es desconocida. Existen diferentes métodos para estimar el error cometido (v. [1] y sus referencias). Este trabajo se centra en el estudio y el diseño de estimadores a posteriori (frente a los estimadores a priori, que determinan las propiedades de convergencia del método pero no dan un valor preciso del error que se ha cometido realmente) en cantidades de interés (por oposición a los estimadores en norma energética que son los más clásicos). Pensando en extender su dominio de aplicación a problemas de dinámica estructural, interesará utilizar estimadores cuyo coste computacional sea bajo. Por ello, este estudio renuncia, de entrada, a obtener estimadores de error que proporcionen cotas superiores o inferiores del error (v. [2], [3] y [4]) y se conforma con lograr aproximaciones del error razonables. Se busca obtener un estimador que sea rápido y económico en términos de coste computacional, y que pueda ser utilizado en problemas que requieran estimar el error en múltiples ocasiones sobre un mismo dominio, como ocurre en los procesos adaptativos [5] o en los problemas transitorios [6].
Para estimar el error en cantidades de interés definidas por el usuario, se introduce un problema auxiliar (adjunto o dual) que permite representar la cantidad de interés en términos de medidas energéticas del problema original (directo o primal) y el adjunto [7–9]. Las cantidades de interés se describen con funcionales lineales que se utilizan como término fuente (cargas aplicadas) en el problema adjunto. En este trabajo, se presentan representaciones residuales, las cuales utilizan la solución del problema original o primal y de un problema auxiliar o dual para estimar el error. Se muestran las representaciones exactas del error en la cantidad de interés y se proponen estimaciones utilizando varias funciones de tipo burbuja, las cuales se centran en el análisis del error en los desplazamientos y no en las tensiones.
Los estimadores de tipo recovery se basan en el posprocesamiento de la solución numérica. Se reconstruye una solución incrementando su regularidad y mejorando su precisión, lo cual sirve para estimar el error cometido en la solución numérica original. La super-convergent patch recovery, introducida en [10], es una de las opciones más populares y utilizadas para estimar el error en la norma energética y también en cantidades de interés. Aunque, en principio, no proporciona cotas del error, se han introducido variantes que sí las proporcionan (v. [11] y [12]). Se trata de una familia de estimadores eficiente y robusta. Sin embargo, su extensión a los problemas de dinámica es más compleja, puesto que, en sus versiones originales, no contempla reconstruir campos de desplazamientos o velocidades, únicamente de tensiones.
En este artículo, se estudia el posible uso de los estimadores residuales en este contexto, con vistas a su generalización al problema dinámico. Se opta por una estrategia residual semiexplícita, en la cual se espera que el coste sea moderado. En la perspectiva de la dinámica estructural, renunciar a las propiedades de acotación que proporcionan los estimadores implícitos no es grave (en todo caso, obtener cotas es muy costoso y estas no son muy precisas cuando el amortiguamiento es pequeño [13]). El carácter semiexplícito se asocia a una aproximación implícita del error en el problema adjunto y a una estimación explícita del error en la cantidad de interés a partir de una representación del error que consiste en evaluar el residuo débil primal (explícito) en la estimación del adjunto (o dual). Esta naturaleza semiexplícita hace que el método sea poco costoso, y se prevé extenderlo a problemas de dinámica estructural en que el coste computacional es un cuello de botella. En la resolución de la parte implícita (estimación del problema adjunto), y también para el cálculo del estimador explícito, es necesario introducir una malla refinada en que se calculan las aproximaciones del error y de los residuos. Se consideran dos estrategias distintas para construir esta malla refinada: una con una submalla estructurada y otra, no estructurada. Existen varias formulaciones de estimadores residuales implícitos, todas ellas basadas en resolver problemas locales con diferentes estrategias para descomponer el dominio e imponer las condiciones de contorno correspondientes (v. [14] y [15]). En este trabajo, se opta por construir la estimación con funciones positivas de soporte compacto, denominadas funciones de tipo burbuja (v. [14] y [16]), y analizar el error relativo y el índice de efectividad correspondiente.
La estructura del resto del artículo es la siguiente: en la sección 2, se presenta el problema a resolver, así como la forma de medir los errores de aproximación. En la sección 3, se introducen la cantidad de interés y el problema auxiliar utilizado en las representaciones del error que se dan en la sección siguiente. En la sección 4, se analiza un nuevo término que aparece en las representaciones del error en cantidades de interés cuando se utilizan submallas distorsionadas. En la sección 5, se presenta una estimación residual implícita del error [4], para la cual se requiere resolver un problema sobre cada elemento del dominio y un problema por cada arista. Para reducir el coste computacional de este procedimiento, en la sección 6 se presentan dos estimaciones del error utilizando funciones de tipo burbuja. Los resultados numéricos se presentan en la sección 7 y las conclusiones, en la sección 8.
2Planteamiento del problemaSea Ω⊂R2 un dominio abierto y acotado con frontera continua y función lipschitziana ∂Ω, compuesta por una parte del tipo Dirichlet y otra del tipo Neumann, es decir, ∂Ω=ΓD¯∪ΓN¯, con ΓD¯∩ΓN¯=∅. El dominio está ocupado por un material elástico lineal, que está sometido a la acción de una fuerza f en [L2(Ω)]2 y a tensiones de borde h en L2ΓN2. Los desplazamientos u y las tensiones σ=σu en el cuerpo elástico satisfacen la condición de equilibrio y las condiciones de contorno siguientes:
donde n denota la normal unitaria exterior a ΓN. El tensor de tensiones σ(u) está relacionado con el tensor de deformacionespor medio de la ley de Hookedonde C es el tensor de constantes elásticas. Con el método de los residuos ponderados de Galerkin, se obtiene la forma variacional del problema, que consiste en encontrar una función u en el espacio V tal quedonde a(u,v)=∫Ωσu:εvdΩ, lv=∫Ωf⋅vdΩ+∫ΓNh⋅vdΓ y V=w∈H1Ω2:w=0 en ΓD, con H1(Ω), que denota un espacio de Sobolev. La norma en energía se define por
Para aproximar la solución u, se construye un espacio de elementos finitos VH, inducido por una partición del dominio Ω en elementos Ωk,k=1,⋯,nelem, tales que Ω¯=∪kΩ¯k, con Ωk∩Ωj=∅, para k≠j. La solución de elementos finitos uH pertenece al espacio VH y satisface la ecuación
El error numérico asociado a uH, se denota por e y está definido por e=u−uH. Este error pertenece al espacio V cuando uH es interpolado en V y satisface la ecuación del error
donde RP se denomina residuo débil asociado a uH o residuo primal. Cuando no se dispone de una solución exacta del problema, se utiliza una malla más fina, denominada malla de referencia, para obtener una solución de referencia, uref. Con esta solución, el error se define por eref=uref−uH, donde uH es interpolada sobre la malla de referencia.3Cantidad de interés y problema dual
Sea J un funcional lineal. Se quiere calcular el error en la evaluación de J(uH). Para ello, se introduce un problema auxiliar, que consiste en encontrar una función u* en el espacio V, tal que
El error numérico asociado a una solución de elementos finitos uH* de este problema dual, calculada sobre VH, se denota por ε y está definido por ε=u*−uH*. Este error pertenece al espacio V cuando uH* es interpolado en V y satisface la ecuación del errordonde RD se denomina residuo débil asociado a uH* o residuo dual. Una de las cantidades de interés analizadas en este trabajo es un valor puntual de la solución. Se quiere utilizar la solución numérica del problema para obtener un valor suficientemente preciso de la solución en un punto x0∈Ω dado. Para lograrlo, se utiliza el funcionaldonde q=qxqyT representa un vector que indica la dirección y el sentido de una fuerza en el problema dual. La magnitud J(u)(≃qxuxx0+qyuy(x0)) representa el promedio de los desplazamientos u en la dirección de q, en un pequeño entorno de x0. La función de suavizado, Wr[17] se escoge del modo siguiente:con r>0, y la constante c=c(r,x0) se escoge de modo que ∫ΩWrx−x0dΩ=1.4Representaciones residuales del error
Si se dispone de la solución exacta del problema, el error en la cantidad de interés descrita a través del funcional J se puede representar de tres formas distintas: utilizando el producto escalar de la energía de deformación, a través del residuo primal y a través del residuo dual. Pero, si no se dispone de la solución exacta, se puede encontrar una solución de referencia, calculada sobre una malla suficientemente fina, que sea una buena aproximación de la solución y, en consecuencia, del error exacto. En este caso, si con la malla fina no se genera un subespacio funcional que contenga el subespacio asociado a la malla gruesa, aparecerá un nuevo término en las representaciones del error en cantidades de interés.
4.1Representaciones exactas del errorSi uH y uH* son aproximaciones en VH de los problemas primal y dual dados en las ecuaciones (2) y (5), el error en la cantidad de interés asociado a estas soluciones se puede representar de las tres formas siguientes [7]:
donde a(*,*) es el producto escalar energía de deformación; RP*, el residuo primal, y RD(*), el residuo dual, dados en las ecuaciones (2), (4) y (6). Para obtener estas representaciones, se parte de la validez de la ortogonalidad de Galerkin,Como e∈V, escogiendo v=e en la ecuación (5),
Teniendo en cuenta que u*=uH*+ε y la ortogonalidad de Galerkin, se sigue que
Por otro lado, como ε∈V, se puede utilizar en la ecuación (4) para obtener
Finalmente, utilizando la ecuación (6),
4.2Análisis del término RPuH*
Las representaciones del error que se muestran en el apartado 4.1 se basan en la propiedad de ortogonalidad de Galerkin. Cuando se calcula una solución aproximada sobre una malla gruesa y esta es interpolada sobre una malla refinada encajada a la anterior, y, además, los espacios funcionales asociados a ambas mallas están encajados, la solución interpolada satisface la propiedad. Sin embargo, si no se cumple la contención de los espacios, el término RPuH* puede ser no nulo, debido a que, al calcular las integrales involucradas en dicho término, se generan errores, ya que la solución interpolada no puede representarse de forma exacta sobre la malla fina. Este fenómeno se observa cuando se utilizan funciones bilineales típicas para elementos cuadriláteros. En esta sección, se analiza el término RPuH*, calculado sobre una malla distorsionada de elementos cuadriláteros encajada en una malla gruesa. Se quiere ver el efecto de la distorsión de una malla refinada en el cálculo de la integral de una función de forma proyectada sobre dicha malla y su consecuencia sobre el término RPuH*, dado en las representaciones del error en cantidades de interés.
Sea N1ξ,η=0,25 (1−ξ)(1−η) una función de forma típica en el cálculo de elementos finitos, definida sobre el cuadrado −1,1×[−1,1] en el sistema local ξη. Esta función se puede enviar a Ω=−1,1×−1,1 en un sistema global xy, teniendo en este la representación N1x,y=0,251−x1−y. Calculando, se obtiene que ∫ΩN1dΩ=0,25.
Si se escoge una partición de Ω que tenga cuatro elementos con sus aristas ortogonales, Ω=∪j=14Ωj, entonces, como consecuencia de la contención de los subespacios asociados a ambas mallas y de la bilinealidad de las Ni, se sigue que
Por el contrario, si se escoge una partición del dominio formada por elementos distorsionados, como la que se muestra en la figura 1, el valor de la integral será distinto de 14. En la figura 2, se pueden observar los contornos de las superficies correspondientes a la parte de N1 definida sobre los elementos de la malla distorsionada. La traza de la gráfica de la función de forma N1 definida sobre el contorno distorsionado está formada por arcos de curvas (cuadráticas), mientras que la traza correspondiente a la proyección está formada por segmentos de recta. Luego, al calcular la integral de N1 utilizando dicha aproximación, se generan errores de proyección. Observando la figura 2, se puede inferir que, al utilizar una malla refinada con elementos distorsionados, aunque sea encajada, el espacio asociado a la malla distorsionada no va a contener al espacio asociado a la malla grosera, pues no se puede representar N1 de forma exacta sobre la malla fina. Esto, a su vez, hace que, en la representación del error en la cantidad de interés, se genere un nuevo error, que es captado exactamente con el término RPuH*, como se evidencia en los resultados numéricos.
Cabe destacar que esta situación no se presenta cuando la malla refinada está formada por elementos no distorsionados. No obstante, en procesos de adaptatividad, utilizando elementos cuadriláteros con funciones de forma bilineales, se pueden generar submallas con elementos distorsionados, lo cual justifica el análisis anterior.
4.3Nuevas representaciones del errorAl proyectar una solución sobre una malla refinada que tenga elementos distorsionados, como la que se puede observar en la figura 1, se generan errores en el cálculo del término RPuH*. Estos errores inciden en la representación del error en cantidades de interés, lo cual se expresa en las representaciones siguientes:
Aquí, las soluciones aproximadas al problema, así como la de uH*, son interpoladas sobre la malla de referencia utilizada para calcular los errores.
Para verificar la primera representación, en la ecuación (5) se toma u*=uH*+ε, de donde av,uH*+ε=Jv. Luego, por la linealidad de a*,*, J(v)=a(v,uH*)+a(v,ε). Escogiendo v=e en la última ecuación, se tiene J(e)=a(e,uH*)+a(e,ε), pero
Por tanto,
Para verificar la segunda representación, se escoge v=εen la ecuación (4):
Por tanto,
Para verificar la tercera representación, se parte de la ecuación (6): RDv=Jv−av,uH*, donde se escoge v=e, para obtener
Por tanto,
Las representaciones residuales del error presentadas anteriormente permiten utilizar cualquier estimador del error en los desplazamientos para calcular aproximadamente el error en la cantidad de interés.
5Estimación implícita del errorUna forma de calcular el error de una solución de elementos finitos es resolver problemas locales para aproximar las dos contribuciones, la debida al error en el interior de los elementos y la debida al error sobre las aristas. La suma de las soluciones de estos problemas locales se denota por e¯=e¯int+e¯edg≈e.
Para calcular aproximadamente los errores interiores, se escribe e¯int≈∑k=1neleme¯kint, donde e¯kint representa la contribución del elemento Ωk al error interior. Luego, utilizando la ecuación del error (4) sobre cada uno de los elementos de la malla, se resuelve un problema con condiciones de contorno de Dirichlet homogéneas, tomando como término fuente la restricción del residuo débil al elemento Ωk. El problema local consiste en resolver, sobre una malla refinada como la que se muestra en la figura 3, la ecuación
Los operadores ak*,* y RkP* representan las restricciones de los operadores correspondientes al elemento Ωk. Para calcular aproximadamente el error debido a la contribución sobre las aristas, se escribe e¯edg≈∑l=1nedge¯ledg, donde e¯ledg representa la contribución de la arista l al error y nedg es el número de aristas. Para cada una de las aristas, se escoge un subdominio Λl⊂Ω que la contenga, como el representado en la figura 4. Este subdominio estará contenido en dos elementos si la arista es interior y en un solo elemento si es una arista de la frontera de Ω. El problema consiste en resolver, sobre una malla refinada como la que se puede observar en la figura 4, la ecuación:
tomando condiciones de Dirichlet homogéneas sobre el contorno. El error sobre arista encontrado como la solución a dicho problema se ortogonaliza de los errores interiores mediante la fórmula:donde los subíndices k y k' corresponden a los índices de los elementos Ωk y Ωk' que intersecan el subdominio Λl. Una vez calculadas las contribuciones sobre el interior y sobre las aristas, el error estimado se escribe como
Cuando se ha estimado el error implícitamente, se utiliza el residuo dual para obtener una estimación del error en la cantidad de interés, según la representación dada en la ecuación (10):
Para tener un mapa del error estimado sobre todo el dominio, la distribución se hace por elementos. El error interior RDe¯kint es asignado al elemento Ωk. El error sobre arista RD(e¯ledg) es asignado a uno o dos elementos; Si la arista es interior, la contribución se distribuye de forma equitativa entre los dos elementos que intersecan el subdominio Λl. Si la arista es frontera, la contribución es asignada completamente al elemento que contiene el subdominio.
Procediendo de forma análoga a la descrita anteriormente, se obtiene otra estimación del error en la cantidad de interés, pero utilizando el error del problema auxiliar:
Para escribir cada una de las dos representaciones anteriores, se recuperan los errores estimados elemento a elemento y se calculan los residuos correspondientes sobre el interior y sobre las aristas.
6Estimación con función burbujaSe extiende el uso de las funciones de tipo burbuja a problemas vectoriales, como los presentadas en [16], donde se estima el error en cantidades de interés utilizando funciones de tipo burbuja para problemas elípticos, siguiendo las ideas de [14]. Las estimaciones son de tipo residual dual ponderado, como las presentadas en [18] y [19]. En [16], se estima el error en cantidades de interés, y el error en desplazamientos se aproxima como la suma de múltiplos de funciones de tipo burbuja, definidas con soporte compacto sobre todos los elementos de una malla de elementos finitos. Esta suma consta de dos contribuciones al error, una sobre el interior de los elementos y otra sobre las aristas. Luego, el error en la cantidad de interés se estima utilizando el error aproximado y los residuos de los problemas primal y dual según las representaciones allí definidas. En este trabajo, se ensaya con varias funciones de tipo burbuja diferentes, para analizar tanto el comportamiento del error relativo como el de los índices de efectividad en cada caso. Se trabaja con dos tipos de funciones con soporte compacto. Las primeras se escogen con derivada no nula en la dirección normal sobre el contorno de cada elemento y las segundas, con derivada nula en la dirección normal sobre el contorno de cada elemento. El análisis de los resultados obtenidos permitirá determinar cuál de todas las anteriores es la función de tipo burbuja más adecuada para estimar el error.
6.1Estimación interior y sobre aristas con funciones del tipo burbujaPara estimar el error en cantidades de interés, se utilizan dos soluciones de elementos finitos. Una correspondiente al problema original o primal (ecuación 2) y otra correspondiente a un problema auxiliar o dual (ecuación 5), que se resuelven sobre una misma malla (aunque se pueden resolver sobre mallas diferentes). Se utilizan las representaciones del error definidas en las ecuaciones (9) y (10), utilizando el residuo de los problemas primal y dual. Los residuos de las representaciones operan sobre una malla refinada o de referencia, que se escoge encajada a la malla gruesa. La estimación de este error con funciones de tipo burbuja se hace en dos etapas. En la primera etapa, se estiman los errores interiores, eint, correspondientes a las contribuciones al error sobre el interior de los elementos. Esta aproximación se denota por e˜int. Para ello, se supone que, sobre cada elemento, el error interior se puede escribir como un múltiplo de una función burbuja, de soporte compacto sobre el elemento:
en que la función burbuja se da a priori, sin tener que resolver un problema de Dirichlet como en el caso de la estimación implícita; ν es el radio de Poisson, y las constantes ck se han de determinar. Esta hipótesis permite que el cálculo aproximado del error interior consista en realizar solo productos de matrices con vectores. Aunque se perderá un poco la precisión de la estimación del error, el cálculo será mucho más rápido, sobre todo cuando se tenga un gran número de elementos, o cuando haya que calcular el error en muchas ocasiones. Entonces, a partir de la ecuacióny escogiendo v=ψk, se sigue quede donde se puede obtener ck y el error interior queda expresado como
En la segunda etapa, se estima la contribución al error sobre las aristas, eedg. Esta aproximación se denota por e˜edg. Para ello, sobre cada subdominio Λl se escoge una función burbuja χl de soporte compacto y se supone que, sobre dicho subdominio, el error sobre arista se aproxima como un múltiplo de la burbuja. Esto es:
En la ecuación aeint+eedg,v=RPv, se utilizan las aproximaciones del error ae˜int+e˜edg,v=RPv. Por la linealidad del producto escalar, se sigue que ae˜edg,v=RPv−ae˜int,v. Luego, escogiendo v=χl y teniendo en cuenta la ortogonalidad de las funciones burbuja (por ser escogidas con soporte compacto), se tiene que
de donde se puede obtener dl y, denotando la restricción del residuo débil al subdominio Λl por RlP*, se sigue que
Una vez obtenida la estimación sobre aristas, se ortogonaliza a partir de la estimación interior con el proceso de Gram-Schmidt. La estimación del error primal, denotada por e˜, es definida como la suma de las dos contribuciones, es decir, e≈e˜=e˜int+e˜edg. Con esta estimación y la representación dada en la ecuación (10), se obtiene la estimación del error en cantidades de interés:
donde RDe˜int=∑k=1nelemRPψkRDψkaψk,ψk y RDe˜edg=∑l=1nedgRPχl−aχl,e˜intaχl,χlRDχl y las componentes de las funciones ψk=ψkx,1+νψky, para todo k=1,…,nelem y χl=χlx,1+νχly, para todo l=1,…,nedg, son del tipo burbuja, definidas sobre el interior y sobre las aristas de los elementos, respectivamente.
Para estimar el error del problema dual, se procede de manera análoga a como se ha hecho anteriormente. En la primera etapa, se estiman los errores interiores, ε˜int, y, en la segunda etapa, se estiman los errores sobre aristas, ε˜edg. La estimación del error dual, denotada por ε˜, se define como la suma de estas dos contribuciones, es decir, ε≈ε˜=ε˜int+ε˜edg. Con esta estimación y la representación dada en la ecuación (9), se obtiene otra estimación del error en cantidades de interés, definida por:
donde RPε˜int=∑k=1nelemRDψkRPψkaψk,ψk y RPε˜edg=∑l=1nedgRDχl−aχl,ε˜intaχl,χlRPχl.6.1.1Función burbuja
En este trabajo, se utilizan diferentes funciones de tipo burbuja para determinar cuál de ellas ofrece la mejor estimación del error, tanto en energía como en cantidades de interés. Estas funciones se definen en un sistema de coordenadas (ξ,η) local. Luego se utilizan transformaciones de las usuales en el método del elemento finito para enviar dichas funciones a cada uno de los elementos de la malla o elementos globales. En esta ocasión, sobre las aristas, se parte de estas mismas funciones, pero transformadas de manera adecuada con rotaciones, traslaciones y dilataciones en el sistema ξ,η, para que, al ser transformadas al elemento global, sean de soporte compacto sobre un subconjunto Λl que contenga la arista como el que se muestra en la figura 4. La burbuja sobre arista definida en el sistema local es enviada al elemento global Λl.
6.1.2Función burbuja escogida a posteriori (a partir de la solución de los problemas locales)Una de las opciones a la hora de escoger funciones de tipo burbuja es resolver problemas locales sobre cada elemento de la malla de elementos finitos, que es la estimación implícita descrita en el apartado 5. Con este proceso, se obtienen buenos resultados, aunque es costoso porque hay que resolver un problema por cada elemento y uno por cada arista.
6.1.3Función burbuja escogida a prioriEl primer tipo de función burbuja interior escogida a priori es como la presentada en [16]:
que tiene una derivada no nula en la dirección normal del contorno y un segundo tipo de función burbuja escogida a priori, de manera tal que tenga una derivada nula en el contorno:7Resultados numéricos
En este ejemplo, se tiene un dominio de material elástico lineal cuadrado, de dimensiones 1 × 1 m2, con un orificio cuadrado en el medio de 0,5 × 0,5 m2. El módulo de Young es E = 100MPa; el coeficiente de Poisson es ν = 0,1 y la densidad es ρ = 1.000kg/m3. Se imponen condiciones de Dirichlet homogéneas en el extremo izquierdo. Para el problema primal, se aplica una fuerza de 250N, repartida en el segmento derecho de la parte superior del dominio, que resulta de aplicar una presión de 1.000N/m sobre los 0,25 m del lado superior situados sobre la columna vertical. Las condiciones de contorno, la fuerza aplicada sobre el dominio y la deformada del problema primal (para la solución en los desplazamientos) se pueden observar en la figura 5. El problema se resolvió con el método de los elementos finitos, utilizando elementos cuadriláteros con funciones de forma bilineales, mediante Cast3M [20]. La cantidad de interés escogida para el análisis es el desplazamiento del punto x0=(0,75;0,75). Para ello, se resuelve un problema auxiliar, que consiste en aplicar sobre el punto x0 del dominio una fuerza j de módulo 1 N. Utilizando el funcional dado en la ecuación (7), se escoge q=12,1√2. Las condiciones de contorno son iguales a las del problema primal. Las condiciones de contorno, la fuerza aplicada y la deformada (para la solución en los desplazamientos) del problema dual se pueden observar en la figura 6.
El problema se resolvió sobre una malla muy fina de 202.800 elementos, sobre la cual se obtuvo un valor de referencia para la norma energética de los desplazamientos de uE=3,7827 y, para la cantidad de interés evaluada en dicha solución, de Ju=−5,45171×10−2. El problema se resolvió utilizando cuatro mallas más gruesas y, con el fin de no utilizar el espacio ocupado por la solución de referencia indicada anteriormente, en cada prueba se recalculó una nueva solución de referencia sobre la malla utilizada para la estimación. Los resultados pueden observarse en la tabla 1, donde el número de elementos de la malla de cálculo y la utilizada para la estimación son denotados por nelem y neles, respectivamente; el error en la cantidad de interés, por J(e), y las representaciones, por Rep1=RDe+RPuH*,Rep2=RPε+RPuH*, RP(ε) y RDe. El residuo primal evaluado en la solución del problema dual es denotado por RPuH*; la cantidad de interés evaluada en la solución, por JuH; la norma energética del error primal, de la solución, de la estimación implícita, de la estimación explícita con la burbuja ψ1 y de la estimación explícita con la burbuja ψ2, por eE, uHE, e¯E, e1E y e2E, respectivamente. Haciendo el análisis en energía en este problema, se observa que, con la estimación implícita y con la burbuja ψ1, se obtiene un índice de efectividad con una media del 81%, mientras que, con la burbuja ψ2, la media es del 61%. También se puede ver que, con una malla de 1.728 elementos, se obtiene un error relativo del 12%; es necesaria una malla de más de 120.000 elementos para disminuir este error a menos del 2%. Sin embargo, el propósito de este trabajo se centra en la cantidad de interés y, en este caso, el error relativo obtenido es del 2% con 1.728 elementos.
Resultados de referencia sobre diferentes mallas y norma energética de los errores estimados
| nelem | neles | Je=Rep1=Rep2 | JuH | eE | uHE | RDe=RPε | RPuH* | e¯E | e1E | e2E |
|---|---|---|---|---|---|---|---|---|---|---|
| 108 | 1296 | -5,08071E-3 | -4,67E-2 | 1,1089 | 3,55 | -5,49E-3 | 4,17E-4 | 0,8854 | 0,8765 | 0,6698 |
| 192 | 2304 | -3,66567E-3 | -4,88E-2 | 0,9279 | 3,62 | -3,95E-3 | 2,93E-4 | 0,7545 | 0,7466 | 0,5678 |
| 432 | 5184 | -2,31288E-3 | -5,09E-2 | 0,7177 | 3,69 | -2,48E-3 | 1,75E-4 | 0,5878 | 0,5873 | 0,4341 |
| 1728 | 20736 | -1,09475E-3 | -5,28E-2 | 0,4652 | 3,74 | -1,16E-3 | 0,74E-4 | 0,3770 | 0,3877 | 0,290 |
En la estimación del error, se ensayó con tres tipos diferentes de funciones de tipo burbuja para observar la influencia de estas sobre los resultados.
Ensayo 1. Sobre cada uno de los elementos de una malla gruesa, se resuelve un problema con condiciones de Dirichlet homogéneas, con el residuo como término fuente, y se escoge como burbuja la solución de este problema (una burbuja para cada elemento). Las estimaciones del error en la cantidad de interés dadas en las ecuaciones (14) y (15), Eα, el error relativo ηα=EαJu y los índices de efectividad Θα=EαJ(e)*100 %, para α=1,2, se pueden observar en la tabla 2. La media de los índices de efectividad en este ensayo es del 98,4% y se obtiene convergencia del orden de 1,2. Sobre la malla de 1.728 elementos, el error relativo estimado es inferior al 2%. En este ensayo, se resuelven nelem problemas sobre mallas de 12 elementos cada una y nedg problemas sobre mallas de 12 elementos cada una, para un total de (nelem+ nedg) problemas. Para la malla de 2.700 elementos y 2.880 nodos, esto representa un poco más de 5.580 problemas resueltos.
Estimación implícita sobre mallas uniformes
| nelem | E1 | E2 | η1 | η2 | Θ1 | Θ2 |
|---|---|---|---|---|---|---|
| 108 | -5,00915E-3 | -5,04821E-3 | 9,66E-02 | 9,74E-02 | 98,59 | 99,36 |
| 192 | -3,70759E-3 | -3,7284E-3 | 7,05E-02 | 7,09E-02 | 101,1 | 101,7 |
| 432 | -2,31759E-3 | -2,32338E-3 | 4,35E-02 | 4,36E-02 | 100,2 | 100,4 |
| 1728 | -1,01868E-3 | -1,01488E-3 | 1,88E-02 | 1,88E-02 | 93,05 | 92,70 |
En los dos ensayos siguientes, se resolvió el problema sobre cuatro mallas uniformes, con submallas no estructuradas para calcular las estimaciones e˜j y ε˜j del error primal y el dual, con las burbujas ψj, para j = 1, 2. Los resultados se pueden observar en las tablas 3 y 4, donde las estimaciones son E1j=RDe˜j y E2j=RPε˜j, los errores relativos son definidos por η1j=E1jJ(u) y η2j=E2jJ(u), y los índices de efectividad son definidos por Θ1j=E1jJ(e)*100% y Θ2j=E2jJ(e)*100%, para j=1,2.
Estimación sobre mallas uniformes con ψ1ξ,η
| nelem | E11 | E21 | η11 | η21 | Θ11 | Θ21 |
|---|---|---|---|---|---|---|
| 108 | -3,10714E-3 | -3,20829E-3 | 5,99E-2 | 6,19E-02 | 61,1 | 63,1 |
| 192 | -2,26008E-3 | -2,34881E-3 | 4,30E-02 | 4,46E-02 | 61,6 | 64,0 |
| 432 | -1,39022E-3 | -1,45637E-3 | 2,61E-02 | 2,73E-02 | 60,1 | 62,9 |
| 1728 | -5,97947E-4 | -6,4447E-4 | 1,10E-02 | 1,19E-02 | 54,6 | 58,8 |
Estimación sobre mallas uniformes con ψ2ξ,η
| nelem | E12 | E22 | η12 | η22 | Θ12 | Θ22 |
|---|---|---|---|---|---|---|
| 108 | -2,42236E-3 | -2,45699E-3 | 4,67E-02 | 4,74E-02 | 47,6 | 48,3 |
| 192 | -1,75857E-3 | -1,78471E-3 | 3,34E-02 | 3,39E-02 | 47,9 | 48,6 |
| 432 | -1,08154E-3 | -1,10251E-3 | 2,03E-02 | 2,07E-02 | 46,7 | 47,6 |
| 1.728 | -4,6204E-4 | -4,76812E-4 | 8,56E-03 | 8,83E-03 | 43,5 | 43,5 |
Ensayo 2. Se escogió la burbuja interior, definida por ψ1={ψx1,1+νψy1}, donde ψx1=ψy1=1−ξ21−η2. Con esta elección, se obtuvieron los resultados que se pueden observar en la tabla 3.
En este ensayo, la media de los índices de efectividad para las estimaciones en cantidad de interés es del 60,8% y se obtiene convergencia de orden 1. Sobre la malla de 1.728 elementos, el error relativo de las estimaciones es inferior al 2%.
Ensayo 3. Se escogió la burbuja interior definida por ψ2ξ,η={ψx2,1+νψy2}, donde ψx2=ψy2=14cosπξ+1cosπη+1. Con esta elección, se obtuvieron los resultados que se pueden observar en la tabla 4.
La media de los índices de efectividad en este ensayo es del 46,7% y se obtiene convergencia del orden de 1,2. Sobre la malla de 1.728 elementos, el error relativo de las estimaciones es inferior al 1%.
En la figura 7, se puede observar el mapa de los errores estimados con E2. Para dibujarlo, se procedió de la forma siguiente. Los errores interiores se asignaron por elementos tal como se construyeron. Los errores sobre aristas se distribuyeron de manera que, si el subdominio utilizado se interseca con dos elementos, la contribución se reparte a partes iguales para cada uno de los elementos involucrados, mientras que, si el subdominio está contenido en un solo elemento, como es el caso de las aristas frontera, la contribución se asigna completamente a dicho elemento. Los errores son mayores en las esquinas interiores, donde surgen singularidades, producto de las tensiones producidas por las fuerzas aplicadas y la geometría del dominio. Para analizar estos sectores, se podría utilizar un método de refinamiento adaptativo, utilizando las estimaciones presentadas aquí, para colocar un mayor número de elementos en estas zonas y tratar de disminuir el error de la solución aproximada de la cantidad de interés; sin embargo, no se ha procedido así en este trabajo, porque estamos pensando en utilizar estas estimaciones en problemas transitorios. En la figura 8 se puede ver en escala logarítmica, el error relativo de la solución, cuando el problema es resuelto utilizando mallas uniformes, donde J(e) representa el error en la cantidad de interés. El orden de convergencia de cada método tiene la misma magnitud que el error en la cantidad de interés de referencia, que en este problema es de 1.2.
El índice de efectividad de las estimaciones analizadas se puede observar en la figura 9, donde la notación es la indicada al comienzo del apartado. En dicha figura, puede apreciarse que la norma en energía del error estimado implícitamente tiene un índice de efectividad muy similar al obtenido con la estimación con la burbuja ψ1. También puede observarse que, para la cantidad de interés, es la estimación implícita la que presenta mejores índices de efectividad, seguida de la estimación con la burbuja ψ1 y, por último, de la estimación con la burbuja ψ2. Aunque los mejores resultados se obtienen con la estimación implícita, esta tiene el costo computacional más elevado, pues se resuelven nelem+nedg problemas de contorno mientras que, al estimar el error utilizando funciones burbuja, solo se realizan productos de matrices con vectores sobre cada elemento.
8Conclusiones
Se ha estimado el error ensayando con tres funciones de tipo burbuja, definidas con soporte compacto. La primera es implícita, escogida como la combinación ortogonal de dos funciones, en que una es la solución de un problema de Dirichlet sobre cada elemento de la malla y la otra es la solución de un problema sobre un subdominio que contiene cada arista de la malla. La segunda y la tercera fueron escogidas a priori, y una de las diferencias entre ellas es que las derivadas normales en el contorno del dominio son nulas para una y no nulas para la otra. Entre los tres tipos de funciones, las implícitas dieron mejores índices de efectividad, aunque el orden de convergencia fue similar con todas las estimaciones.
El problema analizado presenta singularidades debidas a las esquinas reentrantes presentes en el dominio, lo cual hace necesario un número grande de elementos para obtener una buena aproximación. Aunque, en energía, este hecho es muy evidente, pues se necesitan cerca de 20.000 elementos para alcanzar el régimen de convergencia asintótico, en el análisis de la cantidad de interés se requieren menos elementos.
A medida que el número de elementos de las mallas se incrementa, los errores relativos de las estimaciones decrecen, lo cual es deseable en un buen estimador. En el error estimado en la cantidad de interés, utilizando implícitamente y explícitamente diferentes tipos de funciones burbuja en el problema mecánico analizado, se obtuvo convergencia del orden de 1,2.
Para la cantidad de interés, el índice de efectividad medido con la estimación implícita tiene un media del 98,4%; con la burbuja ψ1, del 60,8% y, con la burbuja ψ2, del 46,7%.
Se ha analizado el efecto de proyectar una solución desde una malla gruesa de cuadriláteros y funciones de forma bilineales sobre una malla no estructurada y se ha observado que, en esta situación, aparece un nuevo término en las diferentes representaciones del error en cantidades de interés, el cual es fácil de calcular, una vez obtenida una solución del problema dual. Este término afecta el índice de efectividad, pero no el orden de convergencia, y es nulo cuando las mallas pertenecen al mismo espacio.
Una vez obtenidos los resultados del problema analizado, se puede sugerir el uso del método semiexplícito para estimar el error en cantidades de interés en los problemas mecánicos, en aquellos casos en que se requiera una buena precisión, o el método explícito utilizando funciones de tipo burbuja, si se quiere una mayor velocidad de cálculo, aunque sea a costa de perder un poco de precisión. Con el método explícito, no se requiere resolver problemas locales sobre todo el dominio, lo cual permite reducir el tiempo de cálculo, pues solo hay que realizar multiplicaciones de matrices y vectores. Puede ser de gran utilidad para los problemas de adaptatividad o para los problemas transitorios, donde es preciso estimar el error en múltiples ocasiones.
Este trabajo ha contado con el apoyo del Consejo de Desarrollo Científico, Humanístico, Tecnológico y de las Artes (CDCHTA) de la Universidad de Los Andes, a través del proyecto I-1204-09-05-C.