Resumen

En este artículo presentamos formulaciones estabilizadas de elementos finitos para resolver la ecuación de convección-difusión-reacción escalar en el caso de difusión pequeña. Por un lado, resumimos la aplicación de las formulaciones ASGS y OSS basadas en el concepto de métodos variacionales multiescala. Por otro lado, discutimos aspectos de la utilización de elementos de alto orden, centrando nuestros experimentos numéricos en elementos cuadráticos, cúbicos y de cuarto orden. Asimismo, la aplicación del método OSS requiere de la introducción de una proyección, para lo cual introducimos una modificación del elemento simplicial de cuarto orden con una regla de integración numérica asociada.

Palabras clave: Convección-difusión-reacción, métodos de elementos finitos estabilizados, convección dominante, alto orden, cuadraturas nodales

Abstract

In this paper we present stabilized finite element formulations of to solve the scalar convection-diffusion-reaction equation in the case of small diffusion. On the one hand, we summarize the application of the ASGS and OSS formulations based on the concept of variational multiscale methods. On the other hand, we discuss aspects of the use of high order elements, focusing our numerical experiments on quadratic, cubic and fourth order elements. Likewise, the application of the OSS method requires the introduction of a projection, for which we introduce a modification of the fourth order simplicial element with an associated numerical integration rule.

Keywords: Convection-diffusion-reaction, stabilized finite element methods, dominant convection, high order, nodal quadratures

1 Introducción

El contar con equipos computacionales más robustos tanto en memoria RAM como en velocidad de los procesadores ha permitido entre otras cosas incentivar a numerosos investigadores al estudio de métodos numéricos complejos en fluidos que abordan modelos constitutivos cada vez más reales, buscando resolver las ecuaciones sin hacer hipótesis de simplificación por su complejidad. Dentro de este marco, la ecuación de convección-difusión-reacción (CDR) es un modelo matemático de simulación del transporte de una magnitud física que se utiliza entre otras cosas en la simulación del transporte de contaminantes.

Es bien conocido que el método estándar de elementos finitos de Galerkin aplicado a la ecuación de CDR presenta inestabilidades en la solución cuando el término convectivo es dominante frente al término difusivo, motivo por el cual distintos autores han desarrollado formulaciones estabilizadas. Dicha estabilización consiste en añadir términos que dependen de la malla de elementos finitos a los términos de Galerkin. Así, se pueden encontrar en la literatura diferentes métodos, como el método Streamline-upwind/Petrov-Galerkin (SUPG) [1,2,3], el método de Galerkin-mínimos cuadrados (GLS) [4,5,6], el método de las características-Galerkin (CG) [7] y el método de Taylor-Galerkin (TG) [8], entre otros. Una descripción y comparación de estos métodos se puede encontrar en [9], en donde se expone que todos los métodos esencialmente consisten en la adición de un término de estabilización a la formulación original de Galerkin, y que fundamentalmente excepto por menores modificaciones este término se puede escribir como un parámetro numérico, denominado tiempo intrínseco, multiplicado al producto entre el residuo de la ecuación diferencial a ser resuelta y un operador aplicado a la función de prueba.

En el presente trabajo nos centraremos en los métodos llamados variacionales multiescalas (VMS, por Variational Multi-Scale), introducidos por Hughes y otros [10,11] (véase [12] para una revisión de esta teoría). En particular, nos centraremos en la versión más común, a la que nos referiremos como método ASGS (por Algebraic Sub-Grid Scale), y al llamado método OSS (por Orthogonal Subscale Stabilization) [13,14,15,16,17].

Además de describir resumidamente los métodos ASGS y OSS, la contribución fundamental del presente trabajo es analizar su comportamiento con elementos finitos de alto orden, en particular cuadráticos, cúbicos y de cuarto orden, en todos los casos de clase . Para estudiar dicho comportamiento, presentaremos resultados de convergencia en malla usando soluciones manufacturadas y analizaremos cómo son de robustos estos métodos en presencia de capas límite de la solución. Asimismo, en el caso del método OSS es necesario llevar a cabo una proyección en el espacio de elementos finitos, para lo cual resulta especialmente conveniente usar reglas de integración numérica con los puntos de integración en los nodos. Para elementos simpliciales de cuarto orden esto nos lleva a proponer una modificación del elemento estándar, que describimos en detalle para el caso del elemento triangular.

El artículo está organizado como sigue. En la siguiente sección se describe el problema a resolver y las aproximaciones y formulaciones de los métodos variacionales estabilizados ASGS y OSS. En la Sección 3 presentamos el elemento finito triangular de cuarto orden modificado. Las pruebas de convergencia en malla con soluciones analíticas conocidas y las pruebas de estabilidad con capas límite, así como un ejemplo práctico, se presentan en la Sección 4. Finalmente, en base a los resultados obtenidos se extraen algunas conclusiones en la Sección 5.

2 Planteamiento del problema, aproximación y estabilización numérica

En esta sección vamos a describir la aproximación y estabilización numérica de la ecuación escalar transitoria de CDR, con condiciones de contorno y condiciones iniciales conocidas.

2.1 Problema de contorno

Sea un dominio acotado de y sea , con , el dominio temporal. El problema que nos planteamos consiste en encontrar tal que

(1)

con

(2)

sujeta a condiciones de contorno de Dirichlet

condiciones de contorno de Neumann

y condiciones iniciales

En estas ecuaciones, denota la derivada temporal y la derivada con respecto a la -ésima coordenada cartesiana . Hemos usado la notación de Eintein, de manera que los índices repetidos indican suma sobre todas las coordenadas espaciales. Los coeficientes en (2) son la difusión, la componente en la dirección del campo convectivo de velocidades (para flujo incompresible ) y el coeficiente de reacción. es el término fuente conocido y es la componente -ésima de la normal unitaria a la frontera, denotada por , la cual está dada por Las condiciones de contorno de Dirichlet y las condiciones de contorno de Neumann son conocidas en las fronteras y , respectivamente. Las condiciones iniciales y de contorno de Dirichlet verifican la siguiente condición de compatibilidad

2.2 Forma variacional o débil del problema

Sea el espacio de funciones de (funciones en con derivadas en ) que sea anulan en , y sea el espacio de funciones en tales que la norma es . Para simplificar la exposición, supongamos que .

Multiplicando (1) por una función e integrando el término de difusión por partes, se obtiene la forma débil o variacional del problema, la cual consiste en encontrar tal que

(3)

donde

y cumpliéndose la condición inicial en en . En lo que sigue, tomaremos para simplificar la notación.

2.3 Aproximación numérica con el método estándar de Galerkin

La ecuación 3 es la que aproximamos numéricamente en el espacio usando una formulación de elementos finitos. La discretización temporal la llevaremos a cabo mediante un método de diferencias finitas. En principio, podríamos discretizar primero en espacio o en tiempo indistintamente, aunque en nuestro caso llevaremos a cabo primero la discretización temporal.

2.3.1 Discretización temporal en diferencias finitas

Consideremos una partición uniforme del intervalo de tiempo , con , siendo es el tamaño del paso de tiempo, considerado constante. Aquí y en adelante, usaremos el superíndice para indicar el nivel de tiempo .

Aunque la discretización en el tiempo puede realizarse mediante cualquier aproximación de la derivada temporal en diferencias finitas o elementos finitos, aquí utilizaremos la regla trapezoidal generalizada. Para una función genérica , sea

(4)

Aproximaremos

Esta aproximación es de segundo orden en en el caso en que (método de Crank-Nicolson) y de primer orden si . El caso corresponde al método de Euler implícito o el método de diferencias hacia atrás de primer orden, a menudo abreviado como BDF1. Puesto que nuestro interés en este trabajo se centra en la aproximación espacial, en los ejemplos numéricos hemos usado este último método.

En cada paso de tiempo, el problema 3 discretizado en el tiempo, consiste en encontrar tal que

(5)

2.3.2 Discretización de Galerkin en el espacio

Para aproximar numéricamente la ecuación 5 en el espacio por el método de elementos finitos, consideremos una partición del dominio , de manera que y para cualesquiera , , , donde es el número de elementos de la partición y es el dominio del elemento .

Si es un espacio conforme de elementos finitos para aproximar (), entonces el problema discreto que se conoce como el método estándar de Galerkin, que resuelve la ecuación escalar transitoria de CDR 1 con condiciones de Dirichlet y de Neumann homogéneas en la frontera, consiste en: dado encontrar tal que

(6)

para , siendo conocido proyectando la condición inicial en (con el producto escalar de ).

2.4 Métodos de estabilización numérica

Como se mencionó anteriormente, el método estándar de Galerkin presenta inestabilidades numéricas cuando el término convectivo es dominante con respecto al término difusivo. Los métodos que vamos a utilizar para la estabilización numérica se pueden enmarcar en el método variacional multiescalas (VMS) [10,11]. En particular, veremos dos opciones, correspondientes a dos elecciones del espacio de subescalas (véase también [13,15,16]).

2.4.1 El método variacional mutiescalas (VMS)

Aunque el método VMS puede considerar un número arbitrario de escalas en la solución, en la mayoría de las ocasiones basta con considerar dos escalas para diseñar un método de estabilización. La idea básica es descomponer para cada instante de tiempo la variable continua desconocida en una parte resoluble en el espacio de elementos finitos , con , y una parte en la escala de una submalla , con , la cual no puede ser capturada por la malla de elementos finitos, siendo , donde es cualquier espacio para completar en . Para evitar tecnicismos, podemos pensar en y como espacios dimensionalmente finitos, con una dimensión grande. Puesto que representa la componente de que no es reproducida por el espacio de elementos finitos, la llamamos espacio de las subescalas. De este modo, la variable continua desconocida podemos escribirla como , donde es la componente de en el espacio de elementos finitos y es su componente en el complemento (con respecto a un cierto producto interno) de . La misma descomposición es aplicable a las funciones de test.

Aplicando la descomposición mencionada a la ecuación continua 3 se obtiene:

(7)

Utilizando la linealidad de las formas involucradas, la ecuación 7 se puede descomponer en el siguiente sistema de ecuaciones:

(8)
(9)

Para simplificar el método resultante, vamos a asumir que , lo cual implica que la variación temporal de las subescalas en 8 y 9 es despreciable comparada con el resto de términos. En [17] las subescalas que satisfacen esta condición son llamadas cuasi-estáticas. La inclusión o no de la derivada temporal de las subescalas se analiza en detalle en [18], donde se muestra que es esencial para evitar las inestabilidades que pueden producirse para pasos de tiempo pequeños, que en este trabajo no consideraremos.

Con esta suposición el sistema de ecuaciones 8 y 9 se escribe de la siguiente manera:

(10)
(11)

La ecuación 10 corresponde a la escala resoluble en el espacio de elementos finitos y tiene tres términos a su lado izquierdo, donde el primero y segundo términos son la contribución temporal y espacial de del método estándar de Galerkin y el tercero toma en cuenta la influencia de la subescala en .

De la ecuación 11 se obtiene , que es la contribución de la subescala sobre la componente en el espacio de elementos finitos. Para evitar aproximar derivadas de , podemos considerar la siguiente integración por partes en cada elemento:

(12)

Sea el adjunto del operador , dado por

(13)

La ecuación 12 se puede escribir como

(14)

Reemplazando 14 en la ecuación 10, ésta se puede escribir como

(15)

para todo . Ahora bien integrando por partes los términos y de la ecuación 11, ésta se puede escribir como

(16)

Observemos que el primer término de la ecuación 16 se anula, ya que en la frontera de los elementos las componentes normales de los flujos de son iguales pero de signo contrario, debido a que los flujos difusivos deben ser continuos a través de los contornos ínter-elementales. Como resultado de esto, la ecuación 16 es equivalente a encontrar tal que

(17)
(18)

para y para una cierta función , que llamamos el esqueleto de . Es importante observar que la ecuación 17 se verifica para cualquier elemento ortogonal a La presencia del elemento en la ecuación 17 garantiza que pertenezca a , que es lo que implica la ecuación 16 sin el primer término. En otras palabras, garantiza que pertenezca al espacio de subescalas (y no a todo el espacio ).

Las ecuaciones 15 17 y 18 equivalen exactamente al sistema de ecuaciones 10 y 11 en donde hasta ahora no hemos hecho ninguna aproximación. El siguiente paso es hacer una selección de las funciones y una solución aproximada de la ecuación 17 La manera de escoger dichas funciones da lugar a las formulaciones ASGS y OSS que describimos más adelante. Los pasos comunes los describimos a continuación.

Para seleccionar vamos a usar el siguiente argumento. Para tener una aproximación discreta óptima de que dé valores nodales exactos, se podría pedir que las subescalas se anulen en los contornos de los elementos. En problemas unidimensionales esto da condiciones de contorno homogéneas para el problema 17 y 18 en los extremos de cada elemento. Para los casos de más de una dimensión espacial, esto no es posible, puesto que la solución del problema continuo tendría que ser polinómica en los contornos de los elementos para que pudiera coincidir con . Sin embargo, imponer que puede pensarse como una aproximación.

Suponer que equivale a considerar el espacio de las subescalas como un espacio de funciones burbuja, es decir, un espacio de funciones que se anulan en los contornos de los elementos (véase [19,6]). Con esta aproximación, los términos con las integrales en los contornos de los elementos de la ecuación 15 desaparecen, por lo cual el problema a ser resuelto consiste en encontrar para cada instante de tiempo una función tal que

(19)

en donde todavía hace falta determinar . Para ello, tenemos que resolver de forma aproximada el problema 17 y 18 en el espacio de las subescalas de funciones burbuja. Formalmente podemos escribir:

(20)

en cada elemento, donde es el operador inverso del operador con condiciones de Dirichlet homogéneas.

De la ecuación 19 vemos que solo se necesita la componente de sobre el espacio , donde es el espacio de funciones de la forma , con . Esto sugiere aproximar 20 como

(21)

donde es el parámetro de estabilización y que se calcula dentro del dominio de cada elemento. Es una aproximación algebraica al operador .

El diseño de es una de las piedras angulares en el desarrollo de los métodos de elementos finitos estabilizados. Algunas formulaciones algebraicas para el diseño del parámetro de estabilización para la ecuación escalar de CDR en 1D, basadas en el principio del máximo discreto, pueden encontrarse en [10,11,9]. Un posible argumento para obtener es el análisis de Fourier que se presenta en [17,20,21]. En este trabajo, para abordar la utilización de elementos finitos de alto orden, el diseño del que adoptaremos es:

(22)

donde son constantes algorítmicas, que para los ejemplos numéricos hemos adoptado y es coeficiente de difusión, es la norma de la velocidad, es el coeficiente de reacción, es el diámetro del elemento y es el orden polinomial de interpolación.

A falta de determinar la función , y por consiguiente el espacio de subescalas, el problema para la función se obtiene reemplazando la ecuación (21) en (19) y usando la expresión (22).

2.4.2 El método ASGS (Algebraic Sub-Grid Scales)

La idea del método ASGS es tomar como espacio de subescalas el espacio de funciones que son residuos de funciones de elementos finitos, es decir, el espacio de funciones de la forma , con . Esto nos lleva a que .

El problema que finalmente obtenemos a partir de (21) y (19) con es:

(23)

Esta ecuación variacional discreta, a la cual hay que añadirle las condiciones iniciales, la podemos ahora discretizar en el tiempo usando por ejemplo la regla trapezoidal descrita anteriormente. Para ello, habrá que sustituir por y evaluar el resto de funciones dependientes del tiempo que aparecen en .

El parámetro de estabilización para elementos finitos de alto orden está dado por la ecuación 22, el operador por 2 y el operador por la ecuación 13.

2.4.3 El método OSS (Orthogonal Subscales Stabilization)

En este caso, el espacio de subescalas no será el espacio de residuos de elementos finitos como en el caso anterior, sino que se toma como espacio el espacio ortogonal en el sentido de al espacio de elementos finitos, es decir,

(24)

que es una legítima elección que cumple con .

Tomando en cuenta (24), de (17) se sigue que

lo que significa que es una función del espacio de elementos finitos y por lo tanto numéricamente calculable, mientras que estará en el espacio ortogonal a . A este espacio escogido para las subescalas lo denominamos espacio de subescalas ortogonales.

Puesto que , de (21) se tiene que

(25)

Si llamamos la proyección sobre el espacio de elementos finitos pesada con dentro de cada elemento, de (25) se tiene que

(26)

Reemplazando (26) en (21) tenemos

(27)

Sea ahora

(28)

donde es la identidad en y es la proyección sobre el espacio ortogonal al espacio de elementos finitos. La ecuación (27) se puede escribir entonces como

(29)

Podemos reemplazar esta expresión en (15) y usar el hecho de que las integrales sobre los contornos de los elementos son cero, con lo cual obtenemos la formulación estabilizada de elementos finitos OSS que resuelve la ecuación escalar de CDR transitoria, la cual consiste en encontrar en cada instante de tiempo tal que

(30)

Podemos aplicar aquí los mismos comentarios que para 23. Sin embargo, hay una consideración importante a tener en cuenta, y es que puesto que es una función que pertenece al espacio de elementos finitos , su componente en el espacio ortogonal es nula, esto es, la ecuación variacional discreta anterior se reduce a

(31)

2.4.4 Propiedades de estabilidad y convergencia

Aunque los métodos ASGS y OSS difieren en la definición del espacio de subescalas, ambos tienen las mismas propiedades de estabilidad y convergencia. No es nuestro objetivo presentar el análisis numérico de estos métodos, el cual puede encontrarse en [22] (véase también [23,24] y [25] para la teoría general), sino mostrar en qué sentido dichos métodos son estables y convergen, justificar la expresión del parámetro de estabilización (y en particular de su dependencia en ) e indicar cómo debe ser la integración numérica para preservar la convergencia asociada a la integración exacta. Explotaremos este último punto en el apartado siguiente.

Para simplificar la exposición, consideremos el problema estacionario, con coeficientes constantes, y malla de elementos finitos uniforme, de tamaño . En este caso, el parámetro de estabilización dado por (22) será el mismo en todos los elementos.

Los métodos ASGS y OSS se puede demostrar que son estables y convergentes en la norma:

(32)

donde es la norma de . El punto clave para demostrar esta estabilidad es hacer uso del estimador inverso (ver por ejemplo [25])

válido para cualquier función de elementos finitos , siendo una constante que depende de la forma de la malla, pero no de su tamaño ni del orden de interpolación . La potencia 2 de es precisamente quien obliga a dividir por en la contribución del término difusivo en el parámetro de estabilización dado por (22). Con ello tenemos garantizada la estabilidad.

Para analizar la convergencia es necesario hacer uso de las propiedades de interpolación. Si es la norma del espacio de Sobolev , la solución del problema continuo es suficientemente regular y es su mejor aproximación en el espacio de elementos finitos, se tiene que

(33)

donde los casos de interés corresponden a y . Aquí y en lo que sigue, es una constante que no depende de la malla. Con este estimador de interpolación, es posible demostrar que la función de error de las formulaciones que consideramos, dada por

se comporta como

(34)

Que sea la función de error de los métodos ASGS y OSS en la norma (32) quiere decir que

En virtud de la expresión (32) de la norma de trabajo y de la función de error (34), se observa que el estimador de error presentado es óptimo tanto para advección dominante como para difusión dominante. Hay que notar sin embargo que si se considera la convergencia en , hay un rango de parámetros en los que se pierde la optimalidad (véase [22]).

El objetivo de haber presentado los resultados de estabilidad y convergencia en este trabajo es doble. Por un lado, justifica la validez de la expresión del parámetro de estabilización (22) para elementos de alto orden, que es el caso que nos ocupa. Por otro lado, nos permite observar que si la integración es numérica y no exacta (como obviamente siempre es el caso), hay un límite de exactitud mínimo para no perder orden de convergencia de la aproximación. En el apartado siguiente presentaremos una modificación del elemento triangular de cuarto orden con una cuadratura numérica asociada. En este caso, , con lo que de acuerdo con (34) la cuadratura debe ser como mínimo de cuarto orden en el caso de difusión dominante y de orden para convección dominante. De hecho, en presencia de término reactivo () y con reacción dominante se obtiene que la cuadratura mínima debe ser de quinto orden. Este hecho lo usamos a continuación.

3 Elemento finito triangular modificado de cuarto orden

3.1 Introducción

En el método OSS, para calcular la proyección ortogonal del residuo en la ecuación 31, lo que hacemos es calcular la proyección del residuo en el espacio de elementos finitos. Suponiendo para simplificar que los parámetros de estabilización son iguales para todos los elementos, debemos resolver el sistema

(35)

siendo . Cuando expresamos en función de los valores nodales y tomamos como las funciones de forma, (35) da lugar al sistema lineal

(36)

donde es conocido y es la matriz de masa de la partición de elementos finitos. Esta matriz de masa es en principio llena, por lo que es necesario resolver el sistema (36). Sin embargo, es bien sabido que se puede aproximar por una matriz diagonal usando una regla de integración numérica con los puntos de integración en los nodos (cuadratura cerrada) en vez de la integral exacta que aparece en (35). Con ello la resolución del sistema que aproxima a (36) es trivial. El punto clave es que la exactitud de la regla de integración numérica debe preservar la de la aproximación de elementos finitos que se tendría con la integración exacta.

En ocasiones las matrices de masa diagonales obtenidas a través de una cuadratura cerrada presentan en su diagonal valores nulos, lo cual hace que no sean invertibles e imposibilita su uso. En elementos lineales como el (simplicial, es decir, triangular para y tetraédrico para ), esta matriz se obtiene directamente utilizando integración cerrada, con todos sus pesos distintos de cero. Lo mismo sucede para elementos cúbicos, mientras que los elementos cuadráticos tienen que subdividirse en elementos lineales y usar cuadradura nodal para la división resultante; con este procedimiento puede obtenerse una matriz de masa diagonal que, si bien reduciría el orden de convergencia en caso de que se reemplazara por ella la matriz de masa original, puede usarse como precondicionador para resolver el sistema (36).

Atención especial merece el elemento de cuarto orden simplicial. Veremos a continuación que es posible modificar el elemento convencional, en el sentido de alterar la posición de los nodos centrales, y plantear una regla de integración numérica con los puntos de integración en los nodos y los pesos distintos de cero, lo cual nos permitirá aproximar la matriz de masa por una matriz diagonal invertible. Además, dicha regla de integración veremos que es capaz de integrar de manera exacta hasta un polinomio completo de quinto grado, con lo que de acuerdo con la discusión del apartado anterior no se disminuirá el orden de aproximación de la formulación final.

Aunque es posible aplicar la modificación que proponemos tanto para como para , haremos los desarrollos en el primer caso. Plantearemos el elemento modificado considerando móviles los tres nodos interiores del elemento original y parametrizando la posición de los mismos con una variable cuyo valor será incógnita, al igual que los pesos de la regla de integración nodal asociada; la posición de los 12 nodos restantes es la misma que la del elemento original.

3.2 Descripción del elemento P₄ modificado

En la Figura 1 se representa el elemento en el dominio de referencia, dado por el triángulo , , . Las posiciones de los nodos sobre los lados es la del elemento estándar, pero hemos introducido una variable para parametrizar la posición de los tres nodos interiores, obteniéndose así el elemento que llamaremos modificado, con los tres nodos interiores móviles en función del parámetro . Obsérvese que con tendríamos el elemento original.

Elemento P4-modificado con los 3 nodos interiores móviles
Figura 1: Elemento P4-modificado con los 3 nodos interiores móviles

3.3 Integral exacta

Consideremos un polinomio de grado de la forma:

(37)

donde los coeficientes son constantes. La integral exacta de sobre la región triangular de la Figura 1 está dada por:

(38)

Sea

(39)

Reemplazando (39) en (38) tenemos:

(40)

3.4 Integral numérica

La integral exacta vamos a aproximarla mediante una cuadratura cerrada sobre una región triangular con 15 puntos distribuidos de la siguiente manera: 3 puntos (1, 2, 3) en los vértices del triángulo, 9 puntos numerados del 4 al 12 distribuidos simétricamente en los lados del triángulo y los 3 puntos móviles interiores restantes (13, 14, 15) distribuidos simétricamente como se ve en la Figura 2.

Posición de los 15 puntos de integración numérica
Figura 2: Posición de los 15 puntos de integración numérica

Los pesos y están asociados a los puntos detallados en la Figura 3 para los 15 puntos de integración. Estos pesos, junto con el parámetro , son las incógnitas a determinar.

Pesos para los 15 puntos de integración numérica
Figura 3: Pesos para los 15 puntos de integración numérica

Con las posiciones y los pesos de los 15 puntos de integración que se detallan en las Figuras 2 y 3, planteamos la integral numérica, , como sigue:

(41)

De la ecuación (41), reemplazando dado en la ecuación (37) y agrupando términos, obtenemos:

(42)

Puesto queremos que la integral exacta y la integral numérica sean iguales, esto es , de las ecuaciones (40) y (42) se tiene que:

(43)

para todo , , donde es el grado del polinomio completo cuya integral exacta sobre una región triangular se quiere calcular numéricamente.

Las ecuaciones (39) y (43) nos permiten obtener un sistema de ecuaciones no lineal para determinar , , , y . Recordemos que los puntos , , son las coordenadas de los 12 puntos ubicados en los vértices y lados del triángulo respectivamente, y los puntos , , son las coordenadas de los 3 puntos interiores cuya ubicación está parametrizada con la variable .

3.5 Cuadratura cerrada para integración exacta de un polinomio de cuarto grado

Utilizando las ecuaciones (39) y (43) para fácilmente se puede plantear el sistema no lineal de 15 ecuaciones para integrar exactamente un polinomio completo de cuarto grado. De las 15 ecuaciones, obtenemos tres ecuaciones linealmente independientes en términos del parámetro , que son:

Resolviendo este sistema de ecuaciones en términos del parámetro , se tiene:

(44)
(45)
(46)

Para obtener la ecuación parametrizada para , podemos usar la primera ecuación del sistema de 15 ecuaciones mencionado anteriormente. Reemplazando en las ecuaciones (39) y (43), tenemos respectivamente:

(47)
(48)

Igualando las ecuaciones (47) y (48) y reemplazando las ecuaciones (44), (45) y (46), tenemos:

(49)

Las ecuaciones parametrizadas (44), (45), (46) y (49), verifican las 11 ecuaciones restantes del sistema no lineal de 15 ecuaciones planteadas para Por lo tanto, con las ecuaciones parametrizadas para y hemos obtenido un conjunto infinito de soluciones de parámetro . Por ejemplo para tenemos: , donde observamos que todos los pesos son distintos de cero. Como observación, notemos que para , nos da, que es la cuadratura cerrada con los puntos de integración numérica coincidiendo con los nodos del elemento original, donde no todos los pesos son distintos de cero.

3.6 Cuadratura cerrada para integración exacta de un polinomio de quinto grado

En el apartado anterior hemos obtenido un conjunto infinito de cuadraturas cerradas que integran en forma exacta hasta un polinomio completo de cuarto grado sobre una región triangular con puntos de integración coincidentes con los nodos del elemento modificado, con sus tres nodos interiores móviles dependientes del parámetro , y con valores de para los cuales todos los pesos son distintos de cero. Ahora nos proponemos obtener de este conjunto infinito de cuadraturas, al menos una cuadratura que integre exactamente un polinomio completo de quinto grado, también con todos sus pesos distintos de cero.

Utilizando las ecuaciones (39) y (43) para obtenemos 6 ecuaciones no lineales adicionales al sistema no lineal anterior de 15 ecuaciones. De estas 6 ecuaciones adicionales, solo 1 es linealmente independiente con el sistema anterior, la misma que se obtiene reemplazando y en las ecuaciones (39) y (43):

(50)

Reemplazando en la ecuación (50) las soluciones parametrizadas de y dadas por las ecuaciones (49), (44), (45) y (46) respectivamente, se obtiene la ecuación:

(51)

cuya solución es:

(52)

Reemplazando las soluciones parametrizadas de y en las 5 ecuaciones restantes, es decir, para y en la ecuaciones (39) y (43), se obtiene la misma ecuación (51), con lo cual la solución es única para todas las ecuaciones del sistema.

De las dos soluciones de , vamos a encontrar la cuadratura cerrada para:

(53)

Reemplazando el valor de la ecuación (53) en las ecuaciones parametrizadas (49), (44), (45) y (46) se obtienen los pesos y , todos distintos de cero:

Con todo esto, hemos obtenido un elemento modificado, con las posiciones de los nodos interiores dadas por , y cuya cuadratura cerrada sobre una región triangular con puntos de integración coincidentes con los nodos del elemento modificado integra de forma exacta hasta un polinomio completo de quinto grado, siendo todos los pesos distintos de cero.

La matriz de masa aproximada que se obtiene con la cuadratura cerrada sobre el elemento modificado es una matriz diagonal con todos sus elementos distintos de cero. Este elemento modificado es el que utilizamos en todos nuestros ejemplos numéricos para el cálculo de la proyección ortogonal del residuo en el método OSS.

En la Tabla 1 presentamos la cuadratura cerrada del elemento modificado que acabamos de obtener.


Tabla 1. Cuadratura cerrada del elemento modificado
Puntos Pesos
1 0 0
2 1 0
3 0 1
4 1/4 0
5 1/2 0
6 3/4 0
7 3/4 1/4
8 1/2 1/2
9 1/4 3/4
10 0 3/4
11 0 1/2
12 0 1/4
13
14
15



4 Resultados de los experimentos numéricos

4.1 Pruebas de convergencia en malla

Las pruebas de convergencia en malla que hemos llevado a cabo consisten en calcular el error en norma , es decir, la norma de la diferencia entre entre la solución exacta y su aproximación numérica. Considerando solamente refinamiento en , si el error del método se comporta como el error de interpolación, situación que puede considerarse óptima, este error debe comportarse como

(54)

que corresponde al mismo comportamiento que el error de interpolación (33) con (y absorbiendo en la constante la dependencia con ). Consideraremos los casos .

Si graficamos la ecuación (54) (en el caso de la igualdad) en el plano, con en las ordenadas y en las abscisas, tenemos una línea recta en la que es la pendiente teórica óptima de convergencia. Las pruebas de convergencia en malla consisten entonces en verificar que la pendiente calculada sea precisamente . Para ello, seleccionamos como solución exacta a la función polinómica espacio-temporal

e imponemos el término fuente en la ecuación escalar de CDR transitoria (1) para que esta sea precisamente la solución exacta del problema:

donde las constantes de difusión , las componentes de la velocidad de convección y la reacción están dadas en cada uno de los ejemplos, tomando para todos los casos. Elegimos también las condiciones de contorno (todas de Dirichlet) e iniciales para que la solución sea la que hemos elegido.

En las Figuras 4, 5 y 6 presentamos los resultados de los experimentos numéricos de las pruebas de convergencia en malla con solución analítica conocida, tanto para el método ASGS como para el OSS. En cada gráfica se puede apreciar con línea continua la pendiente teórica de convergencia y con línea con apéndices sobre ella la pendiente calculada. Presentamos los resultados para los distintos grados polinómicos de las funciones de forma estudiados, es decir, elementos lineales, cuadráticos, cúbicos y de cuarto orden triangulares y elementos cuadrangulares , respectivamente.

El dominio computacional es el cuadrado y el intervalo de tiempo es . El tamaño del paso de tiempo lo hemos tomado como . Como se verá a continuación, en este caso el error está dominado por la aproximación espacial, no se ve afectado por la aproximación de la discretización temporal. La malla de elementos finitos consiste en triángulos o cuadrados formando una malla regular. Para cada grado polinómico de las funciones de forma o hemos calculado los errores en norma para ocho mallas de tamaño , siendo el número de partes en que se ha dividido cada lado del cuadrado. En la Tabla 2 se muestra el número de elementos y el número de nodos para cada refinamiento, tanto para elementos triangulares como para elementos cuadrangulares. El integrador temporal utilizado es BDF1.

Tabla 2. Refinamiento para elementos triangulares y cuadrangulares
Triángulos Número de nodos
elem.
15 450 256 961 2116 3721
20 800 441 1681 3721 6561
25 1250 676 2601 5776 10201
30 1800 961 3721 8281 14641
35 2450 1296 5041 11236 19881
40 3200 1681 6561 14641 25921
45 4050 2116 8281 18496 32761
50 5000 2601 10201 22801 40401
Cuadrados Número de nodos
elem.
15 225 256 961 2116 3721
20 400 441 1681 3721 6561
25 625 676 2601 5776 10201
30 900 961 3721 8281 14641
35 1250 1296 5041 11236 19881
40 1600 1681 6561 14641 25921
45 2025 2116 8281 18496 32761
50 2500 2601 10201 22801 40401


En el encabezado de las gráficas de las pruebas de convergencia en malla se adjuntan dos filas con los valores de las pendientes calculadas; la primera fila corresponde a los valores de las pendientes de las rectas que pasan por los primeros 5 puntos y la segunda fila son las pendientes de las rectas que pasan por los últimos 5 puntos de los 8 puntos del refinamiento de malla evaluados. De esta manera podemos visualizar numéricamente la tendencia de la convergencia de las pendientes calculadas hacia el valor teórico. Para el cálculo de dichas pendientes se ha utilizado el método de los mínimos cuadrados. La simbología corresponde a las pendientes para elementos triangulares lineales, cuadráticos, cúbicos y de cuarto orden, respectivamente, mientras que para elementos cuadrangulares lineales, cuadráticos, cúbicos y de cuarto orden la simbología es . Los valores de las pendientes teóricas se corresponden con , , , .

Los experimentos numéricos que presentamos son los siguientes: en la Figura 4 convección dominante y una pequeña reacción , en la Figura 5 convección y reacción del mismo orden, con y , y en la Figura 6 una pequeña convección y reacción dominante . Para todos los ejemplos tomamos , como se ha indicado anteriormente.


Review 640781029593-Ej1 triangulos ASGS.png Review 640781029593-Ej1 triangulos OSS.png
Review 640781029593-Ej1 cuadrados ASGS.png Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. k=10⁻³, a=\left[\cos60∘,sin60∘\right], s=10⁻³.
Figura 4: Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. ,
Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. k=10⁻³, a=\left[\cos60∘,sin60∘\right], s=1. Review 640781029593-Ej2 triangulos ASGS.png
Review 640781029593-Ej2 triangulos OSS.png Review 640781029593-Ej2 cuadrados ASGS.png
Figura 5: Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. , ,
Review 640781029593-Ej3 triangulos ASGS.png Review 640781029593-Ej3 triangulos OSS.png
Review 640781029593-Ej3 cuadrados ASGS.png Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. k=10⁻³, a=\left[\cos60∘,sin60∘\right]×10⁻³, s=1.
Figura 6: Convergencia ASGS-OSS, elementos triangulares y cuadrangulares. , ,

Para los dos métodos de estabilización ASGS y OSS, los valores de las constantes algorítmicas se han calibrado a . También es necesario indicar que para el método OSS con elementos triangulares de cuarto orden , hemos usado en el mallado el elemento modificado mencionado en la Sección 3.

En las gráficas de las figuras de convergencia en malla se observa tanto gráfica como numéricamente que las pendientes calculadas son mayores o tienden al valor de la pendiente teórica . Los resultados muestran claramente que las formulaciones ASGS y OSS son capaces de aproximar correctamente el problema de CDR escalar transitorio con todas las combinaciones de convección y reacción dominantes.

4.2 Pruebas con capas límite

A continuación presentamos un grupo de pruebas cualitativas para observar la estabilidad de los métodos de elementos finitos ASGS y OSS en la ecuación escalar transitoria de CDR. Se analizan cuatro casos en el dominio , combinando los valores de término fuente , y con el coeficiente de reacción , y además variando la dirección de la velocidad, considerando un valor de para todos los casos. El valor del coeficiente de difusión es para todos los ejemplos, con el fin de tener convección dominante y poder probar la estabilidad de los métodos. El tamaño del paso de tiempo considerado constante es y el tamaño del elemento para todas las figuras es . En todos los casos excepto el caso de la Figura 8, se imponen condiciones de Dirichlet homogéneas, observando inestabilidades locales en el contorno, como era de esperar, pero siendo las soluciones globalmente estables. Para el caso de la Figura 8, las condiciones de Dirichlet en el contorno son: , , . La condición inicial para todos los casos es en todo .

En cada figura se presentan ocho gráficos que corresponden a los cuatro grados polinómicos de las funciones de forma analizadas para cada método estabilizado de elementos finitos ASGS y OSS. La solución se presenta en , tiempo en el cual se ha alcanzado el estado estacionario usando el integrador temporal BDF1.

En la Figura 7 presentamos el caso con término fuente , velocidad de convección en la dirección del eje , coeficiente de reacción y condiciones de Dirichlet homogéneas. En la Figura 8 se encuentra el caso con término fuente , velocidad de convección inclinada coeficiente de reacción y las condiciones de Dirichlet en el contorno dadas anteriormente. En la Figura 9 encontramos el caso con término fuente velocidad de convección inclinada coeficiente de reacción y condiciones de Dirichlet homogéneas. Y, finalmente, en la Figura 10 tenemos el caso con término fuente velocidad de convección inclinada coeficiente de reacción y condiciones de Dirichlet homogéneas.

Capas_Ej1/Ej1_capas_ASGS_q1 Capas_Ej1/Ej1_capas_OSS_q1
Capas_Ej1/Ej1_capas_ASGS_q2 Capas_Ej1/Ej1_capas_OSS_q2
Capas_Ej1/Ej1_capas_ASGS_q3 Capas_Ej1/Ej1_capas_OSS_q3
Capas_Ej1/Ej1_capas_ASGS_q4 Capas_Ej1/Ej1_capas_OSS_q4
Figura 7: Capas límite con ASGS y OSS, elementos cuadrangulares. , , , y condiciones de Dirichlet homogéneas en todo el contorno.


Capas_Ej2/Ej2_capas_ASGS_q1 Capas_Ej2/Ej2_capas_OSS_q1
Capas_Ej2/Ej2_capas_ASGS_q2 Capas_Ej2/Ej2_capas_OSS_q2
Capas_Ej2/Ej2_capas_ASGS_q3 Capas_Ej2/Ej2_capas_OSS_q3
Capas_Ej2/Ej2_capas_ASGS_q4 Capas_Ej2/Ej2_capas_OSS_q4
Figura 8: Capas límite con ASGS y OSS, elementos cuadrangulares. , , , y condiciones de Dirichlet homogéneas en todo el contorno.


Capas_Ej3/Ej3_capas_ASGS_q1 Capas_Ej3/Ej3_capas_OSS_q1
Capas_Ej3/Ej3_capas_ASGS_q2 Capas_Ej3/Ej3_capas_OSS_q2
Capas_Ej3/Ej3_capas_ASGS_q3 Capas_Ej3/Ej3_capas_OSS_q3
Capas_Ej3/Ej3_capas_ASGS_q4 Capas_Ej3/Ej3_capas_OSS_q4
Figura 9: Capas límite con ASGS y OSS, elementos cuadrangulares. , , y condiciones de Dirichlet homogéneas en todo el contorno.


Capas_Ej4/Ej4_capas_ASGS_q1 Capas_Ej4/Ej4_capas_OSS_q1
Capas_Ej4/Ej4_capas_ASGS_q2 Capas_Ej4/Ej4_capas_OSS_q2
Capas_Ej4/Ej4_capas_ASGS_q3 Capas_Ej4/Ej4_capas_OSS_q3
Capas_Ej4/Ej4_capas_ASGS_q4 Capas_Ej4/Ej4_capas_OSS_q4
Figura 10: Capas límite con ASGS y OSS, elementos cuadrangulares. , , , y condiciones de Dirichlet homogéneas en todo el contorno.


En todos las figuras se observa que el comportamiento es el esperado en concordancia con los datos de cada figura. Además observamos que en todas la figuras existen inestabilidades en los bordes correspondientes a las capas límite, las cuales es importante notar que se atenúan al incrementar el refinamiento polinómico.

En la Figura 11 encontramos una comparación de los métodos ASGS y OSS, realizando cortes sobre las gráficas de capas límite para los dos métodos ASGS y OSS del ejemplo de la Figura 7, y graficando los dos cortes sobre un mismo plano cartesiano. Los cortes que encontramos en la Figura 11 son cortes longitudinales y cortes trasversales. Como era de esperar, aparecen oscilaciones localizadas en las capas límites que no se propagan al interior del dominio. Para eliminarlas habría que introducir algún método de captura de discontinuidades. Se observa también que los picos espúreos del método OSS son mayores que los del método ASGS. Este comportamiento es conocido, como también lo es la mayor exactitud del método OSS, y los resultados presentados aquí corroboran este hecho para elementos de alto orden.

Cortes longitudinales paralelos al eje Cortes transversales paralelos al eje
Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q1 Capas_Ej1/Cortes_transversales/ASGS_OSS_Q5
Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q2 Capas_Ej1/Cortes_transversales/ASGS_OSS_Q2
Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q3 Capas_Ej1/Cortes_transversales/ASGS_OSS_Q3
Capas_Ej1/Cortes_longitudinales/Long_ASGS_OSS_Q4 Capas_Ej1/Cortes_transversales/ASGS_OSS_Q4
Figura 11: Comparación ASGS y OSS, elementos cuadrangulares. Cortes longitudinales y transversales del ejemplo de la Figura 7.


4.3 Transporte de una concentración

Presentamos a continuación los resultados numéricos para un ejemplo menos académico, consistente en el transporte y difusión de una concentración en parte de la superficie de un ducto de directriz . Suponemos que tenemos un flujo constante con velocidad unidireccional en un ducto, como se ilustra en la Figura 12. También asumimos que sobre una región de la pared superior del ducto tenemos una concentración constante de una variable , y sobre el resto de la periferia del ducto, incluyendo la sección transversal a la entrada del flujo, asumimos el valor . Queremos estudiar como se propaga esta concentración dentro del flujo considerado.

Villota Codina 2018a 1462 Fig12.png
Figura 12: Ducto en 3D


Considerando el ancho mucho mayor que la altura y debido a que el flujo es constante en toda la sección transversal del ducto, podemos simplificar el planteamiento del problema a dos dimensiones como vemos en la Figura 13.

Villota Codina 2018a 2163 Fig13.png
Figura 13: Ducto en 2D


El modelo matemático que resuelve el problema planteado es:

Para nuestro problema, el dominio es el rectángulo de dimensiones y y el intervalo de tiempo es .

Las condiciones de Dirichlet en el contorno donde se imponen son:

(55)

En el extremo derecho () no se impone ninguna condición de Dirichlet, sino que se considera libre.

El intervalo de tiempo se discretiza mediante una partición uniforme de tamaño . Como en los ejemplos anteriores, el integrador temporal utilizado es BDF1. La malla de elementos finitos utilizada consiste en 850 elementos cuadrangulares formando una malla regular, con 918 nodos para el elemento 3535 nodos para el 7852 nodos para el y 13869 nodos para el

En la Figura 14 se muestran los resultados de la simulación numérica de este problema en para los dos métodos ASGS y OSS y para los cuatro grados polinómicos de las funciones de forma, y observando en todos los casos que los resultados son los esperados en concordancia con los datos del problema. Se visualiza el efecto de la convección dominante, ya que se observa que la concentración que se prescribe como condición de contorno se transporta al interior de dominio con poca difusión. En ninguno de los casos considerados se aprecian oscilaciones, que es el objetivo de los métodos de estabilización.

Ejemplo_practico/ASGS_Q1 Ejemplo_practico/OSS_Q1
Ejemplo_practico/ASGS_Q2 Ejemplo_practico/OSS_Q2
Ejemplo_practico/ASGS_Q3 Ejemplo_practico/OSS_Q3
Ejemplo_practico/ASGS_Q4 Ejemplo_practico/OSS_Q4
Figura 14: Transporte de una concentración

5 Conclusiones

En este artículo hemos presentado una revisión de los métodos estabilizados ASGS y OSS y su estudio con elementos finitos de alto orden, cuadráticos, cúbicos y de cuarto orden, para resolver la ecuación escalar de CDR en problemas con convección y reacción dominantes. En el cálculo del parámetro de estabilización hemos incluido el grado del polinomio, para de esta manera poder utilizar elementos finitos de alto orden y mantener el orden de aproximación óptimo tanto en como en . Para poder usar de forma eficiente el método OSS con el elemento simplicial de cuarto orden, hemos presentado el procedimiento para construir un elemento modificado, el cual permite aproximar la matriz de masa por una matriz diagonal con todos sus elementos distintos de cero; hemos usado este elemento modificado en todos los ejemplos con elementos de cuarto orden con el método OSS.

Los experimentos numéricos de las pruebas de convergencia en malla para elementos lineales y de alto orden revelan que la convergencia de los dos métodos estabilizados es óptima, para todos los casos de convección y reacción dominantes, tanto para elementos triangulares como para elementos cuadrangulares, y para todos los grados polinómicos de las funciones de forma que hemos considerado. Hemos observado que las pendientes de convergencia del método OSS son ligeramente mayores que las del método ASGS, lo cual indica mayor precisión del método OSS. Además, también se concluye que para un número de nodos dado, los errores son menores cuando mayor es el grado del polinomio de forma utilizado, lo cual indica que para obtener un menor coste computacional para una precisión dada podría ser conveniente utilizar elementos finitos de alto orden. Sin embargo, no ha sido nuestro objetivo estudiar con detalle el coste computacional de los métodos presentados.

En términos generales, los resultados de las pruebas de capas límite han sido los esperados y están en concordancia con las condiciones impuestas en cada caso. Asimismo, hemos observado que el método OSS produce más oscilaciones localizadas en las capas límites, como ya era conocido para el caso de elementos lineales y cuadráticos.

Finalmente, los resultados del ejemplo de la simulación numérica del transporte de una concentración han sido satisfactorios, de acuerdo con el fenómeno físico planteado. En este caso, las diferencias de los resultados entre los métodos ASGS y OSS son insignificantes.

Agradecimientos

R. Codina agradece la ayuda recibida a través del programa ICREA Academia, del Gobierno de Catalunya.

Referencias

[1] A.N. Brooks and T.J.R. Hughes. (1982) "Streamline Upwind / Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equation", Volume 32. Computer Methods in Applied Mechanics and Engineering 199–259

[2] D.W. Kelly and S. Nakazawa and O.C. Zienkiewicz and J.C. Heinrich. (1980) "A note on upwinding and anisotropic balancing dissipation in finite element approximations to convective diffusion problems", Volume 15. International Journal for Numerical Methods in Engineering 1705–1711

[3] C. Johnson and U. Nävert and J. Pitkäranta. (1984) "Finite element methods for linear hyperbolic equations", Volume 45. Computer Methods in Applied Mechanics and Engineering 285-312

[4] T.J.R. Hughes and L.P. Franca and G.M. Hulbert. (1989) "A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations", Volume 73. Computer Methods in Applied Mechanics and Engineering 173–189

[5] L. Franca and S.L. Frey and T.J.R. Hughes. (1992) "Stabilized finite element methods: I. Application to the advective-diffusive model", Volume 95. Computer Methods in Applied Mechanics and Engineering 253–276

[6] F. Brezzi and L.P. Franca and T.J.R. Hughes and A. Russo. (1997) "", Volume 145. Computer Methods in Applied Mechanics and Engineering 329–339

[7] J. Douglas and T. Russel. (1982) "Numerical methods for convection dominated problems based on combining the method of characteristics with finite elements or finite difference procedures", Volume 9. SIAM Journal on Numerical Analysis 871–885

[8] J. Donea. (1984) "A Taylor-Galerkin method for convection transport problems", Volume 20. International Journal for Numerical Methods in Engineering 101–119

[9] R. Codina. (1998) "Comparison of some finite element methods for solving the diffusion-convection-reaction equation", Volume 156. Computer Methods in Applied Mechanics and Engineering 185–210

[10] T. J. R. Huges. (1995) "Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods", Volume 127. Computer Methods in Applied Mechanics and Engineering 387–401

[11] T. J. R Hughes and G. O. Feijóo and L. Mazzei and J. B. Quincy. (1998) "The variational multiscale method-a paradign for computational mechanics", Volume 166. Computer Methods in Applied Mechanics and Engineering 3–24

[12] R. Codina and S. Badia and J. Baiges and J. Principe. (Por aparecer) "Variational Multiscale Methods in Computational Fluid Dynamics, in Encyclopedia of Computational Mechanics". John Wiley & Sons Ltd

[13] R. Codina. (2011) "Finite Element Approximation of the Convection-Diffusion Equation: Subgrid-Scale Spaces, Local Instabilities and Anisotropic Space-Time Discretizations, in Bail 2010 - Boundary and Interior Layers, Computational and Asymptotic Methods, Lecture Notes in Computational Science and Engineering, Eds. C. Clavero, J.L. Gracia and F.J. Lisbona", Volume 81. Springer 85–97

[14] R. Codina. (2000) "Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods", Volume 190. Computer Methods in Applied Mechanics and Engineering 1579–1599

[15] R. Codina. (2000) "On stabilized finite element methods for linear systems of convection-diffusion-reaction equations", Volume 188. Computer Methods in Applied Mechanics and Engineering 61–82

[16] R. Codina. (2001) "A stabilized finite element method for generalized stationary incompressible flows", Volume 190. Computer Methods in Applied Mechanics and Engineering 2681–2706

[17] R. Codina. (2002) "Stabilized finite element approximation of transient incompressible flows using orthogonal subscales", Volume 191. Computer Methods in Applied Mechanics and Engineering 4295–4321

[18] R. Codina and J. Principe and O. Guasch and S. Badia. (2007) "Time dependent subscales in the stabilized finite element approximation of incompressible flow problems", Volume 196. Computer Methods in Applied Mechanics and Engineering 2413–2430

[19] C. Baiocchi and F. Brezzi and L.P. Franca. (1993) "Virtual bubbles and Galerkin/least-squares type methods (Ga.L.S)", Volume 105. Computer Methods in Applied Mechanics and Engineering 125–141

[20] R. Codina and J. Blasco. (2002) "Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales", Volume 4. Computing and Visualization in Science 167–174

[21] J. Principe and R. Codina. (2010) "On the stabilization parameter in the subgrid scale approximation of scalar convection-diffusion-reaction equations on distorted meshes", Volume 199. Computer Methods in Applied Mechanics and Engineering 1386–1402

[22] R. Codina. (2018) "On convergence of stabilized finite element methods for the convection-diffusion equation", Volume 75. SeMA Journal 591–606

[23] P. Houston and E. Süli. (2001) "Stabilised -finite element approximation of partial differential equations with nonnegative characteristic form", Volume 66. Computing 99–119

[24] G. Lube and G. Rapin. (2006) "Residual-based stabilized higher-order FEM for advection-dominated problems", Volume 195. Computer Methods in Applied Mechanics and Engineering 4124–4138

[25] Ch. Schwab. (1998) "- and -Finite element methods. Theory and application to solid and fluid mechanics". Oxford University Press

Back to Top

Document information

Published on 08/02/19
Accepted on 24/01/19
Submitted on 19/06/18

Volume 35, Issue 1, 2019
DOI: 10.23967/j.rimni.2019.01.003
Licence: CC BY-NC-SA license

Document Score

0

Views 114
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?