Una nueva formulación para el problema del transporte por convección-difusión

Héctor Gómez Díaz ([email protected]) Escuela Técnica Superior de Ingenieros de Caminos, Canales y Puertos, Universidad de A Coruña Tutores: Ignasi Colominas Ezponda, Fermín L. Ambarina Martínez

´ PARA UNA NUEVA FORMULACION EL PROBLEMA DEL TRANSPORTE ´ ´ POR CONVECCION–DIFUSI ON

1. OBJETIVOS DE ESTE PROYECTO. Este trabajo se enmarca dentro de la l´ınea de investigaci´ on del Departamento de M´etodos Matem´aticos y Representaci´ on de la Universidad de A Coru˜ na en el ambito de la mec´anica computacional de fluidos. ´ La resoluci´on num´erica de los problemas de fluidos presenta grandes dificultades especialmente cuando se trata de situaciones en las que la velocidad con la que se mueve el fluido es elevada. El M´etodo de Elementos Finitos —el m´etodo num´erico m´as ampliamente utilizado en la resoluci´ on de ecuaciones diferenciales— se ha aplicado con ´exito en multitud de problemas pr´ acticos de ingenier´ıa, pero presenta grandes inconvenientes cuando se utiliza para problemas de fluidos con convecci´on elevada. Cuando se emplea el M´etodo de Elementos Finitos —MEF en lo que sigue— para resolver problemas de elasticidad se obtienen resultados muy precisos incluso con discretizaciones muy bastas y con la formulaci´ on m´ as sencilla posible que es la ponderaci´ on tipo Galerkin. Sin embargo, cuando intentamos abordar un problema de fluidos empleando un m´etodo num´erico, es necesario utilizar mallas fin´ısimas para obtener resultados aceptables. La primera impresi´ on que se puede obtener de esto es que los malos resultados son debidos a la no-linealidad y complejidad de las ecuaciones de Navier-Stokes, no obstante, al intentar resolver problemas simplificados y aparentemente m´ as sencillos los resultados no mejoran demasiado. Esto nos hace pensar que los problemas num´ericos que se producen se deben a la propia naturaleza de las ecuaciones de Navier-Stokes. En cierto modo, el paradigma que resume todo esto es la ecuaci´ on del transporte por convecci´on-difusi´ on, que se puede considerar formalmente como la versi´ on lineal y escalar de las ecuaciones de Navier-Stokes. Esta ecuaci´ on a pesar de su sencillez y linealidad presenta grandes dificultades a la hora de ser resuelta 1

—H´ector G´omez D´ıaz—

2

num´ericamente. Muchos autores —T. J. R. Hughes, A. Brooks, F. Brezzi, L. P. Franca, R. Codina, E. O˜ nate, etc— han desarrollado m´etodos de estabilizaci´ on para esta ecuaci´on del transporte, pero —quiz´ as exceptuando los m´etodos que propone Hughes en el ´ambito de las multiescalas— ninguno de ellos tiene una justificaci´on o sentido f´ısico evidentes. El objetivo b´asico de este trabajo es la estabilizaci´ on de la ecuaci´ on del transporte, pero de un modo justificable y claro. Para ello lo que se ha hecho es substituir la ecuaci´on constitutiva del problema —ecuaci´ on de Fick— por otra m´ as general. Se presentan varias ecuaciones constitutivas alternativas pero la que se estudiar´a en este trabajo con mayor profusi´ on es la ecuaci´ on de Cattaneo. Esta ecuaci´on introduce un t´ermino transitorio seg´ un el cual el flujo difusivo por unidad de masa se va aproximando asint´ oticamente a su valor m´ aximo que es el que precide la ecuaci´on de Fick. Sin embargo, el motivo principal para substituir la ecuaci´on de Fick por la ecuaci´on de Cattaneo es que la ecuaci´ on de Fick predice una velocidad no acotada del transporte de masa por difusi´ on, lo cual contradice, por ejemplo, la Teor´ıa de la Relatividad de Einstein. 2. CONTENIDO DEL PROYECTO. Este proyecto comienza con una breve descripci´ on del fen´ omeno de la difusi´ on y de los principales m´etodos que se han propuesto hasta el momento para estabilizar la ecuaci´on del transporte. Posteriormente, se desarrollan las ecuaciones que rigen el fen´ omeno del transporte empleando las ecuaciones constitutivas de Fick y Cattaneo. Adem´ as, se muestran otras ecuaciones constitutivas m´ as complejas y se analiza su significado. En un an´alisis posterior se demuestra que la ecuaci´ on de Fick nos conduce a la paradoja del transporte por difusi´ on a velocidad no acotada. Sequidamente, se presenta un ejemplo que nos confirma que la ecuaci´ on de Cattaneo resuelve este problema, y por u ´ltimo se analiza el problema unidimensional del transporte a partir de la ecuaci´on diferencial que se deriva de la ecuaci´ on de Cattaneo. Las conclusiones extra´ıdas en este apartado ser´ an fundamentales para el resto del trabajo. A continuaci´on, se desarrollan unas condiciones de estabilidad para el problema del transporte estacionario que se deriva de la ecuaci´ on de Cattaneo. Por medio de estas condiciones seremos capaces de conocer el tama˜ no de malla que debemos emplear para resolver un determinado problema mediante un esquema en diferencias no sesgado, es decir, que no emplee informaci´ on de la soluci´ on anal´ıtica para elegir los puntos de aproximaci´ on. Comenzaremos por desarrollar

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

3

las condiciones de estabilidad para la ecuaci´ on en diferencias que se obtiene empleando el MEF con elementos lineales de dos nodos y ponderaci´ on tipo Galerkin. Posteriormente, veremos que estas condiciones tambi´en son aplicables al m´etodo de diferencias finitas centradas, dado que conducen a la misma ecuaci´ on en diferencias que el m´etodo anteriormente mencionado. Adem´ as, se demuestra que estas condiciones tambi´en son aplicables cuando se emplean elementos de tres nodos y funciones de forma cuadr´ aticas. Posteriomente, se presentan una serie de ejemplos de la resoluci´ on num´erica de la ecuaci´on del transporte que se deriva de la ecuaci´ on de Cattaneo. Con estos ejemplos se pretende mostrar —adem´ as de la validez de las condiciones de estabilidad desarrolladas— el comportamiento num´erico de la nueva ecuaci´ on del transporte. En u ´ltimo lugar se recogen las concusiones extra´ıdas tras la realizaci´ on de este trabajo y se apuntan las futuras l´ıneas de investigaci´ on. ´ DEL FENOMENO ´ 3. DESCRIPCION F´ ISICO. Para comprender el fen´omeno f´ısico del transporte por convecci´ on-difusi´ on es necesario entender y distinguir dos procesos de movimiento de un soluto en un fluido: Difusi´on y convecci´on. Difusi´on es el proceso f´ısico debido al cual el soluto se mueve como resultado del movimiento intermolecular de las part´ıculas de ambas substancias. La convecci´on es el movimiento del soluto debido al movimiento del fluido, por lo tanto cuando el fluido permanezca en reposo no habr´ a convecci´ on. En este proceso las part´ıculas fluidas act´ uan como “portadores” de soluto. Este fen´ omeno es intuitivo y f´acil de comprender.

Fig.1- Transporte por difusi´on. Transporte por convecci´on.

—H´ector G´omez D´ıaz—

4

En ausencia de difusi´on, un soluto introducido en un fluido en movimiento ser´ a arrastrado por el flujo, pudiendo caracterizarse este movimiento por la ecuaci´ on de ondas unidireccional —problema hiperb´ olico—: ∂u ∂u +a =0 ∂t ∂x En esta ecuaci´on a es la velocidad del fluido y u es la concentraci´ on de soluto en el fluido. La soluci´on de esta ecuaci´ on es una onda no amortiguada desplaz´ andose con celeridad a. Ambos mecanismos son de gran importancia en cualquier proceso de transporte pero en muchas ocasiones uno de los dos tiene mayor peso que el otro. Existen numerosos aspectos de la ingenier´ıa en los que aparecen los problemas de transporte por convecci´on-difusi´on, entre ellos se encuentran: — vertidos de contaminantes en medios h´ıdricos y a la atm´ osfera. — flujo de fluidos en explotaci´ on de recursos hidr´ aulicos. — simulaci´on del comportamiento de reservas de petr´ oleo, gas natural, etc. — estudio de transferencias de calor y masa en problemas de ingenier´ıa qu´ımica y nuclear. — etc. La naturaleza de estos procesos tambi´en se ve claramente reflejada en la estructura de las ecuaciones de Navier-Stokes. Sin embargo, puesto que la ecuaci´ on del transporte puede considerarse formalmente la versi´ on lineal y escalar de las ecuaciones de Navier-Stokes, resulta m´ as sencillo estudiar la problem´ atica de estos procesos a trav´es del problema de transporte. ´ CLASICA ´ 4. FORMULACION DEL PROBLEMA La formulaci´on del problema del transporte se fundamenta —como la de todos los fen´omenos f´ısicos— en las ecuaciones de equilibrio y las ecuaciones constitutivas. Para el problema de convecci´ on–difusi´ on tendremos dos ecuaciones de equilibrio —conservaci´on de masa fluida y conservaci´ on de soluto— y una ecuaci´ on constitutiva. Esta ecuaci´ on es un aspecto fundamental de este trabajo. Empezaremos formulando el problema con la ecuaci´ on cl´ asica —ley de Fick— y posteriormente lo reformularemos empleando la ecuaci´ on de Cattaneo. Las ecuaciones del problema son: ∂ρ(x x, t) ∂ρ(x x, t) + a T (x x, t) + ρ(x x, t) · div(a a(x x, t)) = 0 ∂t ∂x x

(1)

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

ρ(x x, t)

∂u(x x, t) ∂u(x x, t) + ρ(x x, t)a aT (x x, t) + div(ρ(x x, t) · q (x x, t)) − f = 0 ∂t ∂x x ∂u (x x, t) ∂x x

q (x x, t) = −K K (x x, t)

5

(2)

(3)

La ecuaci´on (1) es la ecuaci´ on de conservaci´ on de masa fluida, la ecuaci´ on (2) es la ecuaci´on de conservaci´on de masa de soluto y la (3) es la ley de Fick o ecuaci´ on constitutiva del problema. 5. PROPUESTAS ANTERIORES PARA ´ DEL TRANSPORTE. ECUACION

ESTABILIZAR

LA

En este apartado nos referiremos en todo momento a la ecuaci´ on unidimensional del transporte que se deriva de la ley de Fick, ya que es la que se ha utilizado hasta ahora para estudiar el problema del transporte por convecci´ ondifusi´ on. Por lo tanto, el problema tipo que analizaremos ser´ a

a

du d2 u − k 2 = 0, dx dx

x ∈ (0, L)

u(0) = u0 u(L) = uL

(I.6)

Podemos hacernos una idea de lo importantes que pueden llegar a ser las oscilaciones cuando la convecci´ on es elevada sin m´ as que representar los valores obtenidos al aplicar a este problema el MEF con una ponderaci´ on tipo Galerkin. Vamos a emplear una aproximaci´ on mediante una malla de 50 elementos lineales de dos nodos, con lo que resulta un tama˜ no del elemento —para L = 1— de h = 0.02. Adem´as utilizaremos un valor para la velocidad de a = 50 y una difusividad k = 0.1. Los valores prescritos los consideraremos u0 = 0 y uL = 1. Con este ejemplo —Fig. 2—queda clara la inviabilidad de la resoluci´ on num´erica de un problema de transporte con convecci´ on grande. Los m´etodos que se han propuesto hasta ahora para estabilizar esta ecuaci´ on consisten esencialmente en a˜ nadir un t´ermino estabilizador a la ecuaci´ on diferencial

6

—H´ector G´omez D´ıaz—

1.2 Galerkin Analitica 1

0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.8

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Fig. 2- Ejemplo de soluci´on oscilatoria para un problema con convecci´on dominante.

o a la formulaci´on tipo Galerkin que se desarrollar´ a posteriormente —esto u ´ltimo resulta justificable desde un punto de vista matem´ atico ya que en principio tenemos total libertad para elegir las funciones de test, sin embargo, no tiene un significado f´ısico claro—. Algunos de estos m´etodos son: 1. Streamline Upwind Petrov Galerkin (SUPG) o M´ etodo de la Difusividad Artificial [A. N. Brooks, T. J. R. Hughes, 1982] 2. Galerkin Least Squares (GLS) [T. J. R. Hughes, L. P. Franca, G. M. Hulbert, 1989] 3. Subgrid Scale Method (SGS) [T. J. R. Hughes, 1995] [L. P. Franca, C. Farhat, 1994] 4. Characteristic Galerkin Method (CG) [O. Pironneau, 1982] [O. C. Zienkiewicz, R. Codina, 1995] 5. Taylor-Galerkin Method (TG) [J. Donea, 1984] 6. Equilibrio en subdominios finitos (FIC) [E. O˜ nate, 1998] Esencialmente los cinco primeros m´etodos consisten en la adici´ on de un t´ermino estabilizador a la formulaci´on original tipo Galerkin. Este t´ermino estabilizador es de la forma [R. Codina, 1996] Z Ωh

P(ωh )τ R(ph )dΩh

(4)

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

7

siendo — P(ωh ) un operador aplicado al espacio funcional de las funciones de test ωh . — τ un par´ametro con dimensiones temporales. — R(ph ) el residuo de la ecuaci´ on diferencial a resolver al aproximar la soluci´ on u mediante un elemento del espacio de las funciones de prueba ph . d2 uh duh −k 2 R(ph ) = a dx dx con uh (x) =

n X

αi pi (x)

i=1

Con este planteamiento com´ un arrancan todos los m´etodos antes enunciados. El m´etodo de la difusividad artificial (SUPG) fue el primero en emplearse para la resoluci´on de problemas de transporte, siguiendo al principio razonamientos meramente heur´ısticos, sin una base formal firme. Adem´ as, es el m´ as sencillo conceptualemte. 5.1. M´ etodo de la difusividad artificial (Streamline Upwind). Este m´etodo es muy sencillo conceptualmente pero presenta grandes inconvenientes, sobre todo cuando se intenta justificar desde un punto de vista te´ orico ya que modifica las ecuaciones del problema a˜ nadiendo difusividad. Adem´ as, s´olo es posible determinar la cantidad de difusividad que hay que a˜ nadir en el caso unidimensional estacionario y sin t´erminos fuente. 5.2. Formulaciones Petrov–Galerkin. La idea fundamental de este m´etodo, consiste en emplear unas funciones de test distintas de las de forma. Muchos autores [I. Christies, 1976], [J. C. Heinrich, O. C. Zienkiewicz, 1977], [J. C. Heinrich, O. C. Zienkiewicz, 1978], [J. C. Heinrich, O. C. Zienkiewicz, 1979], [J. C. Heinrich, 1980], [G. F. Carey, J. T. Oden, 1986] obtienen estas funciones de test mediante la introducci´ on de un “sesgo” en las funciones de forma en sentido contrario al avance del flujo. Esta es una de las formas de introducir el operador P de (4) sobre el espacio funcional de las funciones de test. De esta manera se consigue un efecto estabilizador similar al producido por el m´etodo de la difusividad artificial. A este tipo de formulaciones se le llama Petrov-Galerkin. En resumen: el m´etodo de Petrov-Galerkin, si bien es correcto desde un punto de vista matem´atico, presenta dos inconvenientes:

—H´ector G´omez D´ıaz—

8

1. Introduce difusi´on num´erica fuera de la l´ınea de corriente —crosswind diffusion— cuando se aplica a problemas bidimensionales o tridimensionales [T. R. J. Hughes, A. N. Brooks, 1979]. 2. Exige conocer la direcci´ on del flujo, con lo cual ya podr´ıamos resolver el problema unidimensional por diferencias finitas sin ning´ un problema sin m´ as que utilizar esquemas descentrados seg´ un la direcci´ on del flujo. 5.3. M´ etodo de Galerkin–Least Squares (GLS). Este m´etodo se obtiene particularizando la expresi´ on (4) para P = L, siendo L el operador que define la ecuaci´ on diferencial. En consecuencia, en el caso que nos ocupa,

L = −a

d d2 + k 2. dx dx

Otra forma de obtener este m´etodo de estabilizaci´ on [Natalia Camprub´ı, 1999] es promediando las matrices elementales que produce el m´etodo de Galerkin con las que produce el m´etodo de m´ınimos cuadrados mediante el par´ ametro temporal τ , lo cual resulta equivalente a emplear como funciones de test ω = p + τ L(p)

Adem´as, este m´etodo es equivalente al m´etodo de difusividad artificial — SUPG— para problemas unidimensionales formulados con elementos lineales de dos nodos. 5.4. Equilibrio en subdominios finitos (FIC). Este m´etodo fue propuesto por E. O˜ nate [E. O˜ nate, 1998] y con ´el pretende demostrar que los t´erminos estabilizadores obtenidos de forma heur´ıstica en los m´etodos anteriores surgen de forma natural en la ecuaci´ on diferencial que rige el problema si se realiza un balance de flujo sobre un subdominio finito del problema. Como inconveniente de este m´etodo destaca que modifica la ecuaci´ on de equilibrio del problema.

9

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

´ BASADA EN LA ECUACION ´ DE 6. UNA NUEVA FORMULACION CATTANEO. La nueva formulaci´on se obtiene a˜ nadiendo a las ecuaciones (1) y (2) la ecuaci´ on de Cattaneo, i.e.,   ∂qq (x x, t) ∂qq (x x, t) ∂u(x x, t) q (x x, t) + β (x x, t) + a (x x, t) = −K K (x x, t) (5) ∂t ∂x x ∂x x En la ecuaci´on anterior β (x x, t) es un tensor cuyas coordenadas en una base ortonormal tienen dimensiones de tiempo. Si hacemos las hip´otesis de medio homog´eneo e is´ otropo y suponemos que el campo de velocidades del fluido es constante el sistema de tres ecuaciones que constituyen la ecuaci´on de conservaci´ on de masa fluida, la ecuaci´ on de conservaci´ on de soluto y la ecuaci´on de Cattaneo se puede reescribir como una ecuaci´ on de segundo orden. Esta ecuaci´on es la siguiente: β

∂ 2u ∂u ∂f ∂∇u + 2βa a· − k∇2 u + βa a · ∇(a a · ∇u) + + a · ∇u = f + β βa a · ∇f (6) 2 ∂t ∂t ∂t ∂t En el caso de un problema unidimensional, la ecuaci´ on es

β

2 ∂ 2u ∂ 2u 2 ) ∂ u + a ∂u + ∂u = 0 ∀x ∈ 0

∂u (x, 0) = g(x) ∀x ∈ < ∂t

(7)

Adimensionalizando la ecuaci´ on (7) se puede obtener un n´ umero adimensional similar al de P´eclet que jugar´a un papel muy importante en el resto del trabajo. Este n´ umero es: He =

a0 l k − βa2

(II.73)

siendo a0 y l una velocidad y un espacio caracter´ıstico respectivamente. En el caso unidimensional estacionario sin t´erminos fuente —que estudiaremos ampliamente— la ecuaci´on es la siguiente: a

du d2 u − (k − βa2 ) 2 = 0 dx dx

(8)

—H´ector G´omez D´ıaz—

10

´ PARA INTRODUCIR LA LEY DE CATTANEO. 7. MOTIVACION Para justificar la introducci´on de la ley de Cattaneo resolveremos un problema de difusi´on pura utilizando la ecuaci´ on cl´ asica de convecci´ on-difusi´ on y la ecuaci´ on que se deriva de la ley de Cattaneo. Veremos que al emplear la ecuaci´ on cl´ asica se producen propagaciones de masa a velocidades no acotadas —lo cual contradice el sentido com´ un— mientras que cuando se emplea la ecuaci´ on que se deriva de la ley de Cattaneo la masa se prpaga a una determinada velocidad finita y calculable. El problema que estudiaremos es el de un vertido puntual en un medio fluido sin convecci´on y un dominio infinito. Este problema est´ a gobernado —cuando se utiliza la ecuaci´on de Fick— por la siguiente ecuaci´ on y condici´ on inicial: ∂ 2u ∂u =k 2 ∂t ∂x

∀x ∈ 0

u(x, t = 0) = δ(x) ∀x ∈ < l´ım u(x, t) = 0, t > 0

x→±∞

y su soluci´on que se puede obtener aplicando una transformada de Fourier es

u(x, t) = √

x2 1 e− 4kt 4πkt

∀x ∈ 0

Si fijamos un tiempo t = τ > 0, podemos definir

u ˜(x) = u(x, τ ) = √

x2 1 e− 4kτ 4πkτ

que es la funci´on de Gauss y por lo tanto cumple u ˜(x) > 0 ∀x ∈ 0. Adem´as, como ya hab´ıamos visto, en el instante inicial u(x, 0) = 0 ∀x 6= 0, es decir, no hay masa de soluto en ning´ un punto que no sea el origen de coordenadas. Con todo esto podemos ver que si fijamos un punto cualquiera x0 se cumple: u(x0 , τ ) > 0 ∀τ > 0

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

11

Fig. 3- Condici´on inicial del problema.

y por lo tanto la velocidad media de las part´ıculas que en el instante t = τ se encuentran en el punto x = x0 ser´ıa v ¯ = xτ0 que no est´ a acotada porque el razonamiento anterior es v´alido ∀τ > 0 y ∀x0 ∈ 0

La soluci´on de este problema se ha obtenido realizando uyna doble transformaci´on integral (Laplace y Fourier) y resulta ser:    √  c c2 t2 −x2      I1 2k  √  c2 t 2  − c c c 1 2 2 2 +  e 2k δ(|x| − ct) + I0 √ , |x| ≤ ct 2 2k t 2k 2k c t − x c2 t2 −x2 u(x, t) =       0, |x| > ct Para ilustrar estos resultados presentaremos un par de figuras en las que se puede apreciar que efectivamente al emplear la ecuaci´ qon de Cattaneo la velocidad difusiva media de las part´ıculas de contaminante es βk .

—H´ector G´omez D´ıaz—

12

Fig. 4- Soluci´on para un tiempo t empleando la ecuaci´on de Cattaneo.

Este ejemplo es suficiente para justificar la introducci´ on de la ecuaci´ on de Cattaneo desde un punto de vista te´ orico; la importancia pr´ actica de esta nueva ecuaci´ on se ver´a cuando se aplique al problema de convecci´ on–difusi´ on. A partir de ahora nos centraremos ya en este problema considerando exclusivamente el caso unidimesional estacionario, para el que se desarrollan unas condiciones de estabilidad. ´ 8. CONSECUENCIAS DE LA REFORMULACION. Las consecuencias m´as importantes de la reformulaci´ on son las siguientes: 1. La propagaci´on de masa se produce a velocidad acotada 2. Aparece una velocidad cr´ıtica del fluido a partir de la cual: 2.1. El problema estacionario no tiene sentido. 2.2. El problema transitorio en un dominio finito se debe formular como un problema de valores iniciales ya que como problema de contorno est´ a mal planteado. Estas dos proposiciones se pueden demostrar de modo te´ orico pero adem´ as se pueden justificar desde un punto de vista f´ısico del siguinete modo: - En situaciones de velocidad supracr´ıtica el soluto se ve arrastrado por el flujo y no es posible alcanzar una situaci´ on de equilibrio por lo que el problema estacionario est´ a mal planteado en este caso. - El problema transitorio en un dominio a mal planteado como p finito est´ problema de contorno cuando |a| > k/β puesto que el soluto no puede propagarse hacia aguas arriba, lo que anula desde un punto de vista conceptual la validez de una condici´ on de contorno aguas abajo.

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

13

Se pueden encontrar situaciones similares a esta en problemas de propagaci´ on de ondas gravitatorias en flujos en l´ amina libre o de propagaci´ on de ondas sonoras en un fluido. 3. Se pueden establecer condiciones de estabilidad para el problema estacionario ya que se pueden definir espacios caracter´ısticos a partir de las celeridades de las dos ondas soluci´ on, que han de ser mayores que el tama˜ no de malla. 4. La soluci´on num´erica es estable pr´ acticamente para cualquier valor posible de la velocidad del fluido. Por tanto, al introducir la ley de Cattaneo no s´olo se pueden explicar las causas del mal funcionamiento de un m´etodo num´erico sino que las soluciones son estables en la pr´ actica totalidad del rango de definici´on de la velocidad del fluido a 9. CONDICIONES DE ESTABILIDAD. Como ya hemos comentado, al introducir la ecuaci´ on de Cattaneo aparece una velocidad del transporte de masa por difusi´ on acotada que da un nuevo significado a los par´ ametros del flujo. Esto tambi´en suceder´ a con las condiciones de estabilidad de los esquemas num´ericos que a partir de ahora ser´ an m´ as claras e intuitivas. Pensemos en un problema de transporte gen´erico con unas condiciones de contorno q que limiten la velocidad del fluido a, de modo que se ha de cumplir |a| < βk para que el problema tenga sentido. Por lo tanto el problema a resolver ser´ıa:

a

du d2 u − (k − βa2 ) 2 = 0; dx dx

x ∈ (0, L)

+ Condiciones de Contorno donde k > 0, β > 0. En cualquier punto del dominio, es decir, ∀x ∈ (0, L) existen dos ondas que se propagan en direcciones opuestas: q  k – Una onda que se propaga hacia aguas abajo con una celeridad β +a q  k – Una onda que se propaga hacia aguas arriba con una celeridad β −a A partir de este punto x ∈ (0, L) podemos definir dos espacios caracter´ısticos, a saber:

—H´ector G´omez D´ıaz—

–β –β

q q

14



k β −a

es el espacio caracter´ıstico aguas arriba.



k β +a

es el espacio caracter´ıstico aguas abajo.

Una vez hechas estas consideraciones vayamos a las condiciones de estabilidad en s´ı. Lo que tienen que asegurar estas condiciones es que el tama˜ no de malla sea menor que los espacios caracter´ısticos correspondientes para que el esquema num´erico pueda “captar” el fen´ omeno. En un esquema no sesgado involucramos puntos tanto aguas arriba como aguas abajo del punto de aproximaci´ on de la ecuaci´ on diferencial (xi ), por lo que se ha de cumplir:  q k a) h < λβ β +a . b) h < λβ

q



k β −a .

o lo que es equivalente  h < min λβ

s

 s  k k + a , λβ −a β β

El m´ınimo de estos dos valores ser´ a el primero o el segundo de ellos seg´ un sea a un n´ umero negativo o positivo. Haciendo una discusi´ on seg´ un el signo de a y refiriendo la condici´on al n´ umero elemental He que denotaremos por γ, se obtiene q q k que la condici´on de estabilidad ∀a ∈ (− β , βk ) es: |a|  |γ| < λ q k 2 β + |a| En los casos en que existen condiciones de estabilidad anal´ıticas se ha comprobado que las condiciones de estabilidad desarrolladas tienden asint´ oticamente a las anal´ıticas cuando la velocidad tiende a la velocidad cr´ıtica (la velocidad a partir de la cual no se propagan part´ıculas de soluto hacia aguas arriba que no exist´ıa empleando la ley de Fick). En los casos en que no se pueden establecer condiciones de estabilidad anal´ıticamente se ha comprobado su buen funcionamiente mediante ensayos num´ericos.

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

15

´ 10. ENSAYOS NUMERICOS. En este apartado realizaremos algunos ensayos num´ericos con la ecuaci´ on unidimensional del transporte que se deriva de la ley de Cattaneo bajo dos condiciones de contorno tipo Dirichlet. Por tanto, el problema que trataremos de resolver ser´a: d2 u du a − (k − βa2 ) 2 = 0; x ∈ (0, L) dx dx u(0) = u0 u(L) = uL

(9)

En todos los ejemplos que se presentan en este cap´ıtulo se utiliza una discretizaci´ on en 20 elementos —con lo que h = 0.05— y una difusividad k = 1. Se emplear´ a una formulaci´on tipo Galerkin con elementos de dos nodos y funciones de forma lineales. Los ejemplos se presentar´an en grupos de tres figuras. En cada uno de estos grupos el tiempo de relajacion β siempre permanecer´ a constante —al igual que el tama˜ no de malla y la difusividad—, pero iremos aumentando el valor de la velocidad del fluido a hasta que se produzcan oscilaciones. 10.1. GRUPO 1 Este grupo viene definido por β = 0.01. Con estos valores de k y β la velocidad cr´ıtica —que en este cap´ ıtulo denotaremos por ac — a partir de la cual el problema q no tiene sentido es ac = βk = 10. Seg´ un las condiciones de estabilidad desarrolladas en el cap´ıtulo anterior se producir´an oscilaciones cuando s h > λβ

 k −a β

ya que vamos a considerar valores de la velocidad a positivos. Como hemos visto en el cap´ıtulo anterior podemos suponer λ = 4, por lo tanto la condici´ on anterior queda del siguiente modo: s h > 4β

s  k k h −a ⇒a> − = 8.75 β β 4β

16

—H´ector G´omez D´ıaz—

Por lo tanto el m´etodo dejar´ a de funcionar “correctamente” cuando a > 8.75, es decir, cuando a > 0.875ac . Con esto podemos ver que este esquema num´erico permite resolver el problema del transporte para el 87.5 % de los valores del rango de definici´on de a. 1.4 Galerkin Analitica 1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Fig. 5- Soluci´on anal´ıtica frente a soluci´on MEF–Galerkin con elementos de dos nodos empleando la ecuaci´on de Cattaneo.

En la primera figura que se muestra —Fig. 5— se representa la soluci´ on num´erica y anal´ıtica para una velocidad a = 7. D´emonos cuenta de que aunque este es un problema con un n´ umero He elemental peque˜ no lo estamos resolviendo perfectamente para un valor de la velocidad relativamente grande ya que la velocidad cr´ıtica es ac = 10. A continuaci´on se presenta la figura 6 que nos muestra la soluci´ on num´erica y anal´ıtica para la velocidad l´ımite a partir de la cual el esquema num´erico empezar´a a fallar, es decir a = 8.75. Como veremos en la siguiente figura (Fig. 6) la soluci´ on num´erica de los problemas m´as convectivos que el anterior es oscilante.

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

17

1.4 Galerkin Analitica 1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Fig. 6- Soluci´on anal´ıtica frente a soluci´on MEF–Galerkin con elementos de dos nodos empleando la ecuaci´on de Cattaneo. 1.2 Galerkin Analitica 1

0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.8

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Fig. 6- Soluci´on anal´ıtica frente a soluci´on MEF–Galerkin con elementos de dos nodos empleando la ecuaci´on de Cattaneo.

10.2. GRUPO 2 Este grupo de ejemplos viene q caracterizado por β = 1. En consecuencia la velocidad cr´ıtica ser´a ac = βk = 1 y la velocidad l´ımite para la cual este esquema num´erico proporciona una soluci´ on estable seg´ un los requerimientos de

18

—H´ector G´omez D´ıaz—

q h = 0.9875. Por lo tanto las condiciones de estabilidad desarrolladas es a < βk − 4β debe ser a < 0.9875ac y en consecuencia la soluci´ on que proporciona este esquema num´erico ser´a admisible para el 98.75 % de los valores posibles de a. Al igual que en los grupos anteriores mostraremos tres figuras que consideramos representativas del comportamiento de la soluci´ on de problema. La primera de ellas (7) ser´a la soluci´on del problema para una velocidad de a = 0.97 1.4 Galerkin Analitica 1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

7- Soluci´on anal´ıtica frente a soluci´on MEF–Galerkin con elementos de dos nodos empleando la ecuaci´on de Cattaneo.

La figura (8) nos muestra el comportamiento de la soluci´ on para una velocidad a = 0.9875 que es el l´ımite impuesto por las condiciones de establilidad. Por u ´ltimo presentamos la figura (9) donde la soluci´ on ya es oscilante. El valor de la velocidad elegida en este caso es a = 0.995.

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

19

1.4 Galerkin Analitica 1.2

1

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Fig. 8- Soluci´on anal´ıtica frente a soluci´on MEF–Galerkin con elementos de dos nodos empleando la ecuaci´on de Cattaneo. 1.2 Galerkin Analitica 1

0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Fig. 9- Soluci´on anal´ıtica frente a soluci´on MEF–Galerkin con elementos de dos nodos empleando la ecuaci´on de Cattaneo.

10.3. CONCLUSIONES ´ NUMERICOS

EXTRAIDAS

DE

LOS

ENSAYOS

— La primera conclusi´on que podemos extraer a la vista de los resultados num´ericos obtenidos es que, en contra de lo que parec´ıa en un principio, alejarnos de la ley de Fick —es decir, aumentar β— estabiliza el problema por el siguiente

—H´ector G´omez D´ıaz—

20

motivo: Aunque al aumentar β estemos haciendo el problema m´ as convectivo, a cambio surge una cota superior para la velocidad decreciente con β. Debemos entender bien esta idea. L´ ogicamente si aumentamos β y mantenemos a constante estamos desestabilizando el problema, pero a la vez estamos cambiando la naturaleza del problema, ya que el t´ermino difusivo est´ a perdiendo peso. Por lo tanto si aumentamos β manteniendo a constante estamos transformando el problema en otro m´ as convectivo y pueden ocurrir dos cosas: q 1) El problema deja de tener sentido ya que a > βk q 2) El problema sigue teniendo sentido ya que a < βk pero es un problema m´ as convectivo —en el sentido de que el n´ umero He elemental γ ha aumentado—. Por lo tanto vemos que el aumentar β manteniendo a constante no nos da una informaci´on clara de como cambia el problema al aumentar el tiempo de relajaci´ on del mismo. Una forma razonable de ver qu´e sucede cuando crece el tiempo de relajaci´on del medio es aumentar β pero no manteniendo a constante sino que a q k ser´ a una funci´on de β. En concreto a = µac = µ β donde µ ∈ (0, 1). Si hacemos esto podemos ver que las condiciones de estabilidad son menos restrictivas contra mayor es β. As´ı, si en la condici´ on s h < λβ

 k −a β

q tomamos a = µ βk , la condici´on de estabilidad queda de siguiente modo: s h < λβ

Si denominamos

k (1 − µ) β

s f (β) = λβ

k (1 − µ) β

la condici´on de estabilidad se transforma en h < f (β) Adem´as, podemos ver que la funci´ on f (β) es una funci´ on creciente para cualquier problema que cumpla µ < 1 —requisito que satisface cualquier problema

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

21

q bien planteado ya que es equivalente a la condici´ on a < βk —. Por lo tanto como hab´ıamos anunciado, los requerimientos de estabilidad son menos restrictivos cuanto mayor es β. Resumiendo: Al utilizar la ecuaci´ on de Cattaneo para el problema de convecci´on–difusi´on podemos obtener la soluci´ on del problema de modo estable —incluso utilizando la ponderaci´ on m´ as sencilla posible para el MEF, que es la formulaci´on tipo Galerkin— pr´ acticamente para todos los valores posibles de la velocidad. En consecuencia, desde un punto de vista pr´ actico, hemos estabilizado el problema, pues el rango de valores de la velocidad para los cuales el problema es inestable es totalmente despreciable, a´ un para valores muy peque˜ nos del tiempo de relajaci´on β. — La segunda conclusi´on es una reflexi´ on general. Se trata de hacer notar al lector que lo realmente importante no es que el m´etodo num´erico nos d´e o no la soluci´ on correcta para cualquier valor de los par´ ametros. Lo realmente importante es que una vez introducida la ley de Cattaneo el problema tiene sentido porque la velocidad del transporte de masa por difusi´ on est´ a acotada y que cuando un esquema num´erico falla conocemos las causas de su mal funcionamiento. Estas causas provienen directamente de la naturaleza ondulatoria del problema y est´ an relacionadas con la velocidad a la que se propagan las ondas de difusi´ on. Para comprender esto podr´ıamos decir —aunque abusando del lenguaje— que tenemos que emplear un esquema num´erico que pueda “captar” las ondas de difusi´ on. 11. CONCLUSIONES. En este trabajo se ha comenzado por situar al lector en el punto de partida del mismo, ya que es, en cierto modo, la continuaci´ on de otras Tesinas realizadas en Escuela T´ecnica Superior de Ingenieros de Caminos, Canales y Puertos de A Coru˜ na. Sin embargo, aunque hemos dicho que se puede considerar como la continuaci´on de otros trabajos, el enfoque es radicalmente diferente. En 1998, se desarroll´o un trabajo [Figueroa, A.] en el que se intentaba estabilizar la ecuaci´on del transporte mediante los m´etodos propuestos por diversos autores en los u ´ltimos a˜ nos —M´etodo de Difusividad Artificial, Formulaciones Petrov– Galerkin, Galerkin–Least Squares,...— para posteriormente realizar un c´ odigo que resolviese casos reales. Un a˜ no m´ as tarde, en 1999, se realiza otro Proyecto T´ecnico [Camprub´ı, N.] en el cual se aplicaba el M´etodo de Elementos Finitos, con diferentes ponderaciones, al problema del transporte. La idea de este trabajo es que las oscilacones que se producen al resolver num´ericamente la ecuaci´ on del transporte son consecuencia de la arbitrariedad conque se eligen las funciones de test en la ponderaci´on tipo Galerkin.

—H´ector G´omez D´ıaz—

22

En este trabajo el planteamiento es esencialmente diferente. La idea es substituir la ecuaci´on constitutiva del problema, puesto que nos conduce a la paradoja del transporte por difusi´ on a velocidad no acotada. La nueva ecuaci´ on que emplearemos ser´a la ecuaci´on de Cattaneo. Seg´ un esta ecuaci´ on, el proceso de difusi´ on es de naturaleza ondulatoria y, por tanto, la masa —energ´ıa en el caso de transmisi´on de calor— se transporta a una velocidad acotada. La idea de aplicar la ecuaci´on de Cattaneo al problema de convecci´ on–difusi´ on es muy innovadora, ya que —hasta donde hemos podido saber— esta ecuaci´ on s´ olo ha sido aplicada a problemas termodin´amicos o de difusi´ on pura. A lo largo del trabajo se analizan con detalle las consecuencias de emplear la ecuaci´ on de Cattaneo en lugar de la de Fick. En primer lugar, observamos que, para cualquier problema, existe una velocidad a partir de la cual el problema de convecci´on–difusi´on no se puede tratar como un problema de contorno y hay que formularlo como un problema de valores iniciales. Este resultado es tambi´en muy satisfactorio, ya que intuitivamente se comprende que el problema de convecci´ on– difusi´ on es, en esencia, un problema de valores iniciales. Adem´ as, este cambio de naturaleza en el problema nos recuerda al que se produce en la propagaci´ on de ondas de gravedad en un fluido en l´ amina libre —gobernado por el n´ umero de Froude— o en la transmisi´on de ondas sonoras en un fluido —gobernado por el n´ umero de Mach—. Tambi´en para el problema de convecci´ on–difusi´ on encontraremos un n´ umero adimensional que determina la naturaleza del problema. Otro aspecto que ha resultado esencial es el mencionado cambio de naturaleza del problema. Hemos pasado de un problema parab´ olico sin una interpretaci´ on clara, a un problema ondulatorio mucho m´ as intuitivo. Esta capacidad para reflexionar sobre el problema en t´erminos de ondas nos ha permitido desarrollar unas condiciones de estabilidad para el problema. Finalmente, se ha dedicado un cap´ıtulo al an´ alisis de los resultados num´ericos que se obtienen al introducir la ecuaci´ on de Cattaneo en el problema. Estos ejemplos han servido para constatar la validez de las condiciones de estabilidad desarrolladas y para demostrar que —en contra de la impresi´ on inicial— el tiempo de relajaci´on es un par´ametro estabilizador del problema. En consecuencia, como conclusi´ on final se puede decir que al substituir la ecuaci´ on de Fick por la ecuaci´on de Cattaneo no s´ olo llegamos a un problema con un mayor sentido f´ısico, sino que, adem´ as, es un problema estable pr´ acticamente para cualquier valor de la velocidad. Adem´ as, tiene la ventaja de tener unas condiciones de estabilidad muy sencillas que nos permiten calcular el tama˜ no de malla que debemos emplear.

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

23

´ 12. L´ INEAS FUTURAS DE INVESTIGACION. Como primera l´ınea de investigaci´ on se puede considerar el estudio m´ as profundo de problemas multidimensionales. Este aspecto se ha dejado un poco de lado en este trabajo para intentar comprender perfectamente el caso unidimensional y tener as´ı una base conceptual s´ olida. A partir de este estudio se podr´ıa desarrollar un c´odigo de elementos finitos que resolviese el problema multidimensional de convecci´on–difusi´ on con la ecuaci´ on de Cattaneo y aplicarlo a casos reales. Esto servir´ıa para comparar los resultados con los de experiencias anteriores, por ejemplo la de 1998 [Figueroa, A.] en la que se realiz´ o un c´ alculo en el puerto de A Coru˜ na. Otra posible l´ınea de investigaci´ on —bastante m´ as compleja que la anterior— ser´ıa el estudio de las relaciones entre las condiciones de estabilidad desarrolladas y los esquemas de aproximaci´on ´ optima para casos unidimensionales. La idea es que con el conocimiento de estas condiciones de estabilidad estamos en una situaci´ on privilegiada para comprender las formulaciones de Petrov–Galerkin y el M´etodo de Galerkin Least Squares. El estudio de este u ´ltimo m´etodo se antoja fundamental para conseguir estabilizar los problemas de fluidos con convecci´ on importante. La causa del buen funcionamiento de este m´etodo para el problema de convecci´on–difusi´on es un tema totalmente abierto, pues reproduce la soluci´ on exacta del problema unidimensional, habi´endose cometido un crimen variacional en su formulaci´on. Finalmente, la l´ınea de investigaci´ on m´ as ambiciosa posible ser´ıa intentar introducir el concepto de relajaci´ on en las ecuaciones de Navier–Stokes. Esto consistir´ıa en modificar las ecuaciones constitutivas del problema para obtener unas nuevas ecuaciones din´amicas estables desde un punto de vista num´erico. Este estudio ser´ıa muy complejo y cabr´ıa preguntarse si no aparecer´ıa, tambi´en en este caso, una cota superior de la velocidad del fluido a partir de la cual la naturaleza del problema cambia. 13. BIBLIOGRAF´ IA A. BROOKS, T.J.R. HUGHES (1982), “Streamline upwind Petrov–Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations”. Computer Methods in Applied Mechanics and Engineering, 32, pp: 199-259 G.F. CAREY, J.T. ODEN (1981), “Finite Elements. Fluid Mechanics. Volume VI”. Prentice Hall.

—H´ector G´omez D´ıaz—

24

M.C. CATTANEO (1958), “Sur Une Forme de L’´equation de la Chaleur ´ Eliminant le Paradoxe D’une Propagation Instantane´e”. Comptes Rendus De L’Acad´emie des Sciences–Series I–Mathematics, 247, pp: 431-433 J.S. CHEN, W.K. LIU, T. BELYTSCHKO (1988), “Arbitrary Lagrangian– Eulerian methods for materials with memory and friction”. T.E. Tezduyar, T.J.R. Hughes, (Eds). Recent Developments in Computational Fluid Dynamics. AMD, 95 I. CHRISTIES, D.F. GRIFFITHS, A.R. MITCHELL, O.C. ZIENKIEWICZ (1976), “Finite Element Method for Second Order Differential Equations with Significant First Derivates”. International Journal for Numerical Methods in Engineering, 10, pp: 1389-1396 R. CODINA (1993), “A Finite Element Formulation for the Numerical Solution of the Convection–Diffussion Equation”. International Center for numerical methods in engineering. Monograph CIMNE No 14, Barcelona. R. CODINA, J. BLASCO, G.C. BUSCAGLIA, A. HUERTA (2001), “Implementation of a stabilized finite element formulation for the incompressible Navier-Stokes equations based on a pressure gradient projection”. International Journal for Numerical Methods in Fluid, 37, pp: 419-444 J. DONEA (1984), “A Taylor–Galerkin Method for Convection Transport Problems”. International Journal for Numerical Methods in Engineering, 20, pp: 20-101 A. FICK (1855), “Uber diffusion”. Chemie, 94, pp: 59-86

Poggendorff’s Annalen der Physik und

L.P. FRANCA, S.L. FREY (1992), “Stabilized finite element methods: II. The incompressible Navier–Stokes equations”. Computer Methods in Applied Mechanics and Engineering, 89, pp: 141-219 M.E. GURTIN, A.C. PIPKIN (1968), “A General Theory of Heat Conduction with Finite Wave Speed”. Arch. Ration. Mech. Anal., 31, pp: 113-126 J.C. HEINRICH (1980), “On quadratic elements in finite element solutions of steady–state convection–diffusion equations”. International Journal of Numerical Methods in Engineering, 15, pp: 1041-1052 J.C. HEINRICH, O.C. ZIENKIEWICZ (1977), “Quadratic Finite Element Schemes for Two–Dimensional Convective Transport Problems”. International Journal for Numerical Methods in Engineering, 11, pp: 1831-1844 J.C. HEINRICH, O.C. ZIENKIEWICZ (1979), “The Finite Element Method and upwinding techniques in the numerical solution of convective dominated flow

Una nueva formulaci´ on para el problema del transporte por convecci´on–difusi´on

25

problems”. Finite Element Methods for Convection Dominated Flows. Ed. T.J.R. Hughes. AMD, 34. American Society of Numerical Engineering, New York. A. HUERTA (1988), “Una Formulaci´ on Arbitrariamente Lagrangiana–Euleriana para Grandes Desplazamientos del Contorno”. Centro Internacional de M´etodos Num´ericos en Ingenier´ıa, Barcelona. A. HUERTA, J. DONEA (2002), “Time-accurate solution of stabilized convection-diffusion-reaction equations: I. Time and space discretization”. Comunications in Numerical Methods in Engineering, 18, pp: 565-573 T.J.R. HUGHES (1995), “Multiscale phenomena: Green’s functions, the Dirichlet–to–Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods”. Computer Methods in Applied Mechanics and Engineering, 127, pp: 387-401 T.J.R. HUGHES (1987), “The Finite Element Method: Linear Static and Dynamic Finite Element Analysis”. Prentice Hall. T.J.R. HUGHES, A.N. BROOKS (1979), “A multidimensional upwind scheme with no crosswind diffusion”. Finite Element Methods for Convection Dominated Flows. Ed. T.J.R. Hughes. AMD, 34. American Society of Numerical Engineering, New York. T.J.R. HUGHES, L.P. FRANCA, G.M. HULBERT (1989), “A new finite element formulation for computational fluid dynamics: VIII. The Galerkin–Least Squares method for advective–diffusive equations”. Computer Methods in Applied Mechanics and Engineering, 73, pp: 173-189 T.J.R. HUGHES, W.K. LIU, T.K. ZIMMERMAN (1981), “Lagrangian– Eulerian finite element formulation for incompressible viscous flows”. Computer Methods in Applied Mechanics and Engineering, 54, pp: 341-355 ˜ E. ONATE (1998), “Derivation of Stabilized equations for numerical solution of advective–difussive transport and fluid flow problems”. Computer Methods in Applied Mechanics and Engineering, 151, pp: 233-265 ˜ E. ONATE (1997), “Formulaci´ on de Elementos Finitos estabilizada para problemas de transporte convectivo y flujo incompresible mediante t´ecnicas de c´ alculo infinitesimal”. M´etodos Num´ericos en Ingenier´ıa. R. Abascal, J. Dominguez, G. Bugueda (Eds.) SEMNI, Barcelona. ˜ E. ONATE (1996), “On the Stabilization of the Numerical Solution for Advective– Diffusive Transport on Fluid Flow Problem”. International Center for Numerical Methods in Engineering. Publication CIMNE, No 81, Barcelona. ˜ E. ONATE, J. GARC´ IA, S. IDELSOHN (1997), “Computation of the Stabilization Parameter for the Finite Element Solution of Advective–Difusive

—H´ector G´omez D´ıaz—

26

Problems”. International Center for Numerical Methods in Engineering. Publication CIMNE, No 100, Barcelona. O. PIRONNEAU (1989), “Finite element method for fluids”. Universit´e Pierre et Marie Curie & INRIA. John Wiley and Sons. O. PIRONNEAU (1982), “On the transport-diffusion algorithm and its applications to the Navier-Stokes equations”. Numerische Mathematik, 38, pp: 309-332 O.C. ZIENKIEWICZ, R. CODINA (1995), “A general algorithm for compressible and incompressible flow. Part I: the split, characteristic based scheme”. International Journal for Numerical Methods in Fluids, 20, pp: 887913 O.C. ZIENKIEWICZ, R.L. TAYLOR (1981), “The Finite Element Method: The Basis”. Butterworth–Heinemann. O.C. ZIENKIEWICZ, R.L. TAYLOR (1991), “The Finite Element Method: Fluid Dynamics”. Butterworth–Heineman.