Apuntes_CFD_Jose_Miguel_Perez (2parcial)

Vista previa del material en texto

DR
AF
T
Notas sobre mecánica de fluidos computacional
Borrador: Versión 0.2.1
Notas sobre mecánica de fluidos computacional
JOSÉ MIGUEL PÉREZ
November 27, 2013
DR
AF
T
DR
AF
T
C O N T E N T S
1 introducción 7
i parte 1 9
2 volúmenes finitos 11
2.1 Formulación integral de la ecuación del transporte 11
2.1.1 Carácter conservativo 13
2.2 Objetivo 13
2.3 Definición del método de los volúmenes finitos 14
2.3.1 Discretización del dominio del problema 14
2.3.2 Discretización de la ecuación del transporte 15
2.3.3 Discretización temporal 28
2.4 Sistema algebraico de ecuaciones 29
2.5 Condiciones de contorno 31
2.6 Errores 31
3 problemas volúmenes finitos 33
3.1 Problemas unidimensionales 33
3.1.1 Ecuación de difusión reacción estacionaria 33
3.1.2 Ecuación de convección difusión reacción estacionaria 39
3.1.3 Ecuación de difusión reacción dependiente del tiempo 45
3.2 Problemas bidimensionales 46
3.3 Mallas no estructuradas 47
3.4 Carácter conservativo del esquema numérico 49
ii parte 2 51
4 integradores de riemann 53
4.1 Sistemas hiperbólicos 53
4.1.1 EDPs de primer orden 53
4.1.2 Definición de sistema hiperbólico 55
4.2 Sistemas Lineales 55
4.2.1 Ecuación de convección/advección lineal 55
4.2.2 Sistema lineal de ecuaciones de tipo hiperbólico 57
4.3 Sistemas no lineales 60
4.4 Aplicación a los Volúmenes Finitos 64
4.5 Esquema numérico para las ecuaciones de Euler 69
4.6 Condición de CFL y condiciones de contorno 71
4.7 Alto orden 72
iii appendix 77
a ecuaciones prototipo 79
3
DR
AF
T
DR
AF
T
P R Ó L O G O
El objetivo del documento consiste en describir la fase de diseño de herramientas soft-
ware en el marco de la dinámica de fluidos computacional (computational fluid dy-
namics, CFD).
Este documento está en fase de escritura y corrección. Diríjase a la dirección de
correo jomipp@torroja.dmt.upm.es para realizar cualquier consulta, sugerencias y/o
corrección de erratas.
5
DR
AF
T
DR
AF
T
1
I N T R O D U C C I Ó N
Las fases de desarrollo de un código depende de muchos factores; alcance del código,
existencia de códigos previos, etc. A grandes rasgos el desarrollo de códigos CFD
suelen implicar las siguientes frases:
1. Análisis del problema
a) Estudio y descripción completa del problema físico a resolver
b) Especificación funcional de la herramienta, estimación de los requisitos com-
putacionales, definición de los criterios de validación, etc.
c) Formulación del problema matemático: Definición de las ecuaciones en de-
rivadas parciales EDP (o EDO), tipo (hiperbólica/parabólica/elíptica), de-
scripción de las condiciones iniciales IC y condiciones de contorno CC.
2. Diseño
a) Discretización del dominio
b) Discretización de las ecuaciones
i. Discretización espacial
ii. Discretización temporal
c) Algoritmos numéricos: resolución de sistemas lineales, interpoladores, etc.
3. Implementación del código
4. Verificación y validación del código
5. Utilización y mantenimiento
Estas notas se centran en la fase de diseño, por lo que se supone que el lector es
capaz de realizar la fase de análisis del problema. Con el objetivo de ilustrar el proceso
a seguir, se utilizan problemas bien conocidos de la mecánica de fluidos. Las fases de
implementación y validación del código se verán en un documento posterior.
7
DR
AF
T
DR
AF
TPart IPA RT E 1
DR
AF
T
DR
AF
T
2
V O L Ú M E N E S F I N I T O S
En este capítulo se realiza una introducción al método de los volúmenes finitos. Este
método se utilizan en muchos programas de dinámica de fluidos computacional (por
ejemplo, Fluent, OpenFoam, Flow3D, etc), dada su fácil aplicabilidad en mallas no
estructuradas, su carácter conservativo y su robustez numérica.
En cuanto al contenido de este capítulo, la parte correspondiente al Teorema del
Transporte de Reynolds se describe y demuestra en el capítulo 2 de Liñán et al. (pági-
nas 7 y 9). La definición del método de los volúmenes finitos y la discusión al respecto
de su carácter conservativo se encuentran en la Sección 4.1 de Ferziger and Peric (2002).
Por otro lado, la reconstrucción de los flujos numéricos se puede ver en más detalle en
el Capítulo 3 de Hrvoje (1996). Esta reconstrucción se encuentra también en la Sección
4.4 de Ferziger and Peric (2002).
2.1 formulación integral de la ecuación del transporte
El teorema del transporte de Reynolds permite calcular la variación con el tiempo de
una magnitud fluida intensiva ligada a un volumen fluido VF, a partir de la variación
de dicha magnitud en un volumen de control VC que se mueve con velocidad uVC,
más el flujo convectivo a través de dicho volumen de control. Más concretamente.
Teorema 1 (Teorema del transporte de Reynolds). Dado un volumen de control fluido VF
que se mueve con velocidad u y un volumen de control VC con velocidad uVC, se cumple la
siguiente relación,
d
dt
Z
VF(t)
F(x, t)dV =
d
dt
Z
VC(t)
F(x, t)dV +
Z
∂VC(t)
F(x, t) (u � uVC) dS , (1)
donde F es una propiedad intensiva del fluido.
⇤
Una demostración del teorema anterior se encuentra en el Capítulo 2 de (Liñán et al.).
Toda ley de conservación de la fluido-dinámica se puede formular a partir de la ley de
conservación de las propiedades definidas en los volúmenes fluidos: r (conservación
de la masa), ru (conservación del momento), E (conservación de la energía). Por lo
tanto, dichas leyes se pueden escribir como,
d
dt
Z
VC
FdV +
Z
∂VC
F · dS =
Z
VC
SdV , (2)
donde F puede ser r, ru y/o E y F es igual a F (u � uVC). En la expresión anterior se
ha añadido la posibilidad de que existan términos fuentes en el volumen de control.
11
DR
AF
T
volúmenes finitos
Definiendo F = rf y asumiendo por simplicidad que uVC = 0, se llega a la siguiente
expresión,
Z
VC
∂rf
∂t
dV +
Z
VC
r · (rfu) dV =
Z
VC
r · (Grf) dV +
Z
VC
SfdV , (3)
donde r es la densidad, G el coeficiente de difusión y Sf es el término fuente que en
general depende de f. En la expresión anterior se ha asumido que VC es constante
en el tiempo. Las integrales de los términos convectivo y difusivo pueden reescribirse
como integrales sobre la superficie exterior del volumen de control mediante la uti-
lización del teorema de Gauss o teorema de la divergencia.
Teorema 2 (Teorema de Gauss). El teorema de la divergencia afirma que la divergencia del
campo fu en un volumen de control dado VC es igual al flujo del campo a través de la superficie
que encierra dicho volumen de control, S,
Z
VC
div (fu) dV =
Z
S
fn · udS . (4)
n es el vector normal a la superficie y n · u es la componente de u en la dirección normal al
elemento de superficie dS (ver Figura 1).
⇤
Figure 1.: Volumen de control V definido por la superficie ∂V. El vector normal a la
superficie apunta hacia el exterior del volumen de control.
Por lo tanto, las integrales de los términos convectivo y difusivo definidos en (3)
pueden reescribirse como integrales sobre la superficie,
Z
VC
∂rf
∂t
dV +
Z
S
n · (rfu) dS =
Z
S
n · (Grf) dS +
Z
VC
SfdV . (5)
Notar que n · (rfu) y n · (Grf) son respectivamente los flujos convectivo y difusivo.
En el primer caso el flujo se produce por el efecto de la componente de velocidad n · u.
En el segundo el flujo se produce por la difusión a nivel molecular.
12
DR
AF
T
2.2 objetivo
Es posible obtener la formulación diferencial de las leyes de conservación a partir de
la formulación integral siempre que la solución sea lo suficientemente diferenciable.
El objeto de estudio de este capítulo es la ecuación del transporte (5). Dada la
complejidad de dicha ecuación, ésta se desglosara en distintos problemas:
1. Difusión - Reacción estacionaria:
Z
S
n · (Grf) dS +
Z
VC
SfdV = 0.
2. Convección - Difusión:
Z
S
n · (rfu) dS =
Z
S
n · (Grf) dS.
3. Difusión no estacionaria:
Z
VC
∂rf
∂t
dV =
Z
S
n · (Grf) dS.
4. Convección no estacionaria:
Z
VC
∂rf
∂t
dV +
Z
S
n · (rfu) dS = 0.
El primer y segundo caso representan un problema de tipo elíptico. Como ejem-
plo de la primera se tiene la ecuación de Poisson que sueleaparecer en los métodos
de proyección en presión utilizados en la resolución de las Navier-Stokes incompresi-
ble. La tercera ecuación es del tipo parabólico. Esta ecuación aparece en problemas
de propagación del calor. La cuarta ecuación es de tipo hiperbólico y aparece en el
movimiento de ondas o transporte de sustancias por convección.
La teoría necesaria para estudiar estos problemas se desarrollará en este capítulo
para el caso en el que los flujos convectivos son lineales. Un ejemplo unidimensional
se verá en el Capítulo 3 para los tres primeros tipos de ecuaciones. El caso no lineal se
verá en el Capítulo 4.
Véase el Apéndice A para más información sobre la casuística de ecuaciones asoci-
adas a la ecuación del transporte, (5).
2.1.1 Carácter conservativo
La Ecuación (5) define un ley de conservación sobre cada volumen de control. Despre-
ciando el término fuente esta ecuación afirma que,
2
4
Ritmo de cambio de
la variable en el
VC
3
5 =
[Flujo neto por convención en VC] +
[Flujo neto por difusión en VC] +
[Ritmo de creación de la variable en VC]
Como se verá, esta propiedad conservativa es una de las propiedades más intere-
santes que presenta el método de los volúmenes finitos. En general, los métodos
numéricos que poseen la propiedad de conservación a nivel local también garantizan
la conservación global con tal de que estos sean consistentes.
2.2 objetivo
El objetivo de este capítulo consiste en resolver la ecuación general del transporte en
un dominio dado para unas condiciones iniciales y de contorno dadas. Para ello se
13
(Problema tipo elíptico)
(Problema hiperbólico)
DR
AF
T
volúmenes finitos
utilizará la formulación integral de la ecuación evaluada en un conjunto de volúmenes
de control.
En cuanto al orden del método, como regla general el orden de la discretización
debe ser igual o superior al grado de la ecuación que se pretende resolver. Dado que
las ecuaciones de la fluido-dinámica son de grado dos, se propondrá un esquema a
segundo orden en espacio y tiempo.
2.3 definición del método de los volúmenes finitos
El objetivo de cualquier técnica de discretización, y por tanto del método de los
volúmenes finitos, consiste en transforma un sistema de EDPs (EDOs) en un sistema
algebraico de ecuaciones en el cual las incógnitas representan la solución del problema
original en una posiciones predeterminadas para un tiempo dado.
Para ello es necesario discretizar el dominio, para lo que se deberá especificar los
puntos en los que se resolverá el problema así como los puntos que permiten definir
el contorno.
Definición 1 (Volúmenes finitos). El método de los volúmenes finitos se basa en la for-
mulación integral de las ecuaciones a resolver. En dicho método el dominio se divide en un
recubrimiento de volúmenes de control en los cuales se deben evaluar las ecuaciones integrales.
En cuanto a las variables computacionales, éstas se definen como los valores medios de las
variables problema en cada uno de los volúmenes de control. ⇤
El procedimiento seguido en la resolución de los problemas es el siguiente:
• Aproximación: Integración formal de las ecuaciones sobre los volúmenes de con-
trol y utilización de las variables computacionales.
• Discretización: Sustitución de los integrandos en función de las variables com-
putacionales y resolución de las integrales por cuadraturas.
• Solución del sistema algebraico lineal que se obtiene del paso anterior.
2.3.1 Discretización del dominio del problema
El dominio se divide en un conjunto arbitrario de volúmenes de control VP que se cons-
truyen cada uno alrededor de un punto P (por ejemplo su centroide). Dos volúmenes
de control adyacentes, VP y VN comparten una única cara, véase Figura 2. En la figura,
el vector Sf se aplica en el centroide f de la cara, tiene la dirección normal a ésta,
apunta hacia el exterior del volumen de control y su módulo es igual al área de la
superficie. El vector d une los centroides P y N. El recubrimiento del dominio del
problema con volúmenes de control cumple las siguientes propiedades:
• El dominio del problema, V, es igual a la unión de todos los volúmenes de
control, V = [VC.
• Los volúmenes de control no intersectan: VP \ VN = 0 , 8P, N, tal que P 6= N.
Esto quiere decir que la medida de la intersección entre dos volúmenes adya-
centes es nula en el domino del problema, permitiéndose por lo tanto que dichos
volúmenes adyacentes tengan una cara o puntos en común.
14
DR
AF
T
2.3 definición del método de los volúmenes finitos
Figure 2.: Volumen de control tridimensional VP que tiene al nodo P como centroide.
• En este documento se considera sólo el caso en el que las caras de los volúmenes
de control son planas.
En cuanto a los centroides de los que se ha hablado con anterioridad, estos cumplen
las siguientes propiedades:
• Centroide del volumen de control P:
Z
VP
(x � xP) dV = 0
• Centroide de de la superficie f :
Z
S f
�
x � x f
�
dS = 0
2.3.2 Discretización de la ecuación del transporte
Aproximación de las integrales de volumen
Se definen como variables computacionales o incógnitas del método de los volúmenes
finitos los valores medios de la función en cada uno de los volúmenes de control,
f =
1
VP
Z
VP
f(x)dV . (6)
Con esta definición la ley de conservación (2) queda como sigue:
d
dt
f =
1
VP
✓
�
Z
∂VP
F · dS +
Z
VP
SdV
◆
, (7)
donde F representa los flujos y tomando r = 1.
El objetivo ahora consiste en expresar el lado derecho de la expresión anterior en fun-
ción de las variables problema, f. A esto se le llama cierre de las ecuaciones. Para ello
se tienen que dar los siguientes pasos: (a) reconstrucción de las variables problema
sobre las superficies de los volúmenes de control a partir de las variables computa-
cionales (por ejemplo, mediante procesos de interpolación), (b) evaluación de los flujos
sobre las superficies, (c) discretización de las integrales por cuadraturas (punto medio,
trapezoidal, Simpson, cuadratura de Gauss, etc).
15
P=centroide del volumen de control Vp
f=centroide de la cara
Sf=vector de normal de módulo igual al área de la cara
DR
AF
T
volúmenes finitos
Para el punto (a) se suelen utilizar polinomios interpolantes (caso unidimensional)
o generalizaciones de estos (casos bi- y tridimensional)1 en los que se impone que el
valor medio de éstos sobre los volúmenes de control VC coinciden con el valor de la
variable computacional f en dicho volumen de control. Recordar que la utilización de
un polinomio interpolante de grado n � 1 tiene un error de orden n.
Estos polinomios interpolantes puede usarse para definir el valor f sobre el centroide
del volumen de control, es decir fP. Este valor es igual a f en el caso de considerar
un polinomio interpolante de primer orden (variación lineal). En ese caso, el error
cometido es orden dos. Dicho error se puede obtener aplicando un desarrollo de
Taylor entorno al centoride,
f =
1
VP
Z
VP
f(x)dV =
1
VP
Z
VP
(fP + (x � xP)(rf)P + . . .) dV = fP + O(Dx2) .
ya que
R
VP
(x � xP) dV = 0. Por lo tanto, en un esquema de orden dos, en el que
se asume una variación lineal de f, los valores de las variables se pueden almacenar
en el centroide de los volúmenes de control, lo que define un problema con un perfil
constante a trozos (ver Figura 3).
Figure 3.: Representación de la variable computacional f constante a trozos en un
problema unidimensional. VP representa uno de los volúmenes de control
para el cual el valor promediado de f es f.
Aproximación de las integrales de superficie
Tal y como se ha dicho, en las integrales de superficie aparecen los términos convec-
tivos y difusivos,
Z
VP
r · a dV =
Z
SP
a · dS = Â
f
Z
Sf
a · dS ,
donde a representa el flujo convectivo o difusivo, SP la superficie que encierra el vol-
umen de control, dS el elemento infinitesimal de dicha superficie que apunta hacia el
exterior y Sf es una de las caras planas que define SP.
1 En el caso bidimensional un polinomio interpolante de cuarto orden viene dado por la siguienteexpresión:
a0 + a1x + a2y + a3xy + a4x2 + a5y2 + a6xy2 + a7x2y + a8a2y2
16
DR
AF
T
2.3 definición del método de los volúmenes finitos
Una aproximación de segundo orden viene dada por,
Z
Sf
a · dS = a f
Z
Sf
dS +
 
Z
Sf
(x � x f )dS
!
: (ra) f + . . . ⇡ Sf · af .
En esta expresión se asume que el punto f se encuentra en el centro de la cara. La
aproximación anterior se puede obtener con tal de considerar la regla del punto medio.
Por lo tanto, las integrales de volumen que contienen los términos convectivo y
difusivo se pueden escribir como,
Z
VP
r · a dV =
Z
SP
a · dS = Â
f
Sf · af . (8)
Aproximación de los flujos difusivos
Aplicando la Ecuación (8), los flujos difusivos a = rGrf se aproximan como,
Z
VP
r · (rGrf) ⇡ Â
f
Sf · (rGrf) f ⇡ Â
f
S f (rGf) f · (rf) f .
Aproximación de los flujos convectivos
En cuanto a los flujos convectivos, tomando a = ruf, se tiene que,
Z
VP
r · (ruf) ⇡ Â
f
Sf · (ruf) f .
Gradientes
La expresión del gradiente a segundo orden se obtiene de la Ecuación (8) con tal de
identificar la integral de la izquierda con el valor medio del gradiente evaluado en P
multiplicado por el volumen VP,
(rf)P =
1
VP Âf
�
Sf f f
�
. (9)
Existe otra alternativa a la hora de calcular el rfP válida a orden dos para el caso
unidimensional. Antes de ello se debe definir la notación seguida a la hora de definir
los volúmenes de control en el caso unidimensional, la cual se corresponde con la
notación más extendida en la literatura. Dado un nodo P asociado a una celda (ver
Figura 4) se definen los nodos asociados a las celdas adyacentes como W (oeste) y E
(este). La frontera de la celda VP con las celdas adyacentes se representan con las letras
minúsculas w y e respectivamente. La primera se usa para la frontera izquierda y la
segunda para la derecha. La distancia entre los nodos se representan por dxWP y dxPE,
mientras que la distancia del nodo P a sus contornos izquierdo y derecho se representa
respectivamente por dxwP y dxPe. Finalmente la longitud de la celda P es Dx.
Pues bien, realizando dos desarrollos de Taylor entorno al centroide P para una
malla equi-espaciada,
17
DR
AF
T
volúmenes finitos
Figure 4.: Definición del volumen de control VP. El nodo P es el centro del volumen
de control. Los nodos E y W son los nodos de las celdas adyacentes. La
longitud de la celda es dx y sus fronteras son w y e.
fE = fP + Dx
✓
∂f
∂x
◆
P
+
1
2
Dx2
✓
∂2f
∂x2
◆
P
+ O(Dx3) ,
fW = fP � Dx
✓
∂f
∂x
◆
P
+
1
2
Dx2
✓
∂2f
∂x2
◆
P
+ O(Dx3) ,
y restando ambas expresiones se obtiene que,
✓
∂f
∂x
◆
P
=
fE � fW
2Dx
+ O(Dx2) .
Esta expresión es una alternativa a la Ecuación (9).
Términos fuente
Con este concepto se engloban todos los términos que no pueden tratarse como tér-
minos de convección, difusión o temporales. Tal y como se verá en los ejemplos, este
término puede tener dos contribuciones, una debida a la presencia de un término
fuente sobre el dominio computacional y otra debida al efecto de los contornos. Esto
implicará una dependencia del término fuente con f.
SVP :=
Z
VP
Sf (f) dV = Su + SPfP .
Esta aproximación es exacta si Sf varía linealmente dentro del volumen de control, en
caso contrario es de segundo orden.
Sistema semi-discretizado
Usando las expresiones anteriores en la Ecuación (6) para r 6= 1, se llega a,
d
dt
(rf)P + Â
f
FCf � Â
f
FDf = Su + SPfP ,
donde FCf = S f · (ruf) f y FDf = S f · (rGrf) f . Tal y como se dijo, todas las variables
se calculan y almacenan en el centroides de los volúmenes de control, por lo tanto
los flujos convectivos/difusivos (los cuales se definen en las caras del VC) deben ser
estimados de alguna manera.
18
DR
AF
T
2.3 definición del método de los volúmenes finitos
Reconstrucción de las variables en los contornos de los volúmenes de control
Supongamos que deseamos construir el valor de fe mediante la utilización de poli-
nomios de grado cero. En ese caso imponemos que:
fP =
1
VP
Z
VP
PP(x)dx .
donde PP(x) = aP es el polinomio interpolante utilizado en el volumen de control VP.
Igualmente, uno puede hacer lo propio en el volumen de control E,
fE =
1
VE
Z
VE
PE(x)dx .
donde PE(x) = aE es el polinomio interpolante utilizado en el volumen de control VE.
Como se puede ver, se tienen dos estimaciones de fe, una dada por PP(xe) y otra
dada por PE(xe). En principio dichas estimaciones (de primer orden) son distintas. Ob-
viamente, la discrepancia observada disminuye al refinar la malla. Esta reconstrucción
se representa esquemáticamente en la Figura 5.
Figure 5.: Reconstrucción de fe mediante la utilización de polinomios interpolantes
de grado cero PP y PE. Como se puede ver dicha reconstrucción es distinta a
cada lado de xe, PP(xe) 6= PE(xe).
Para obtener una reconstrucción de segundo orden y de paso evitar la discontinuidad
observada en el caso anterior se pueden utilizar polinomios interpolantes de primer
grado. Dado que en este caso se necesitan evaluar dichos polinomios sobre dos puntos
de control, será necesario expandir el dominio de validez de éstos a dos celdas com-
putacionales adyacentes, ver Figura 6. Tomando como referencia el nodo P, existen
dos posibilidades que nos permiten estimar el valor de fe. Una de ellas utilizaría la
información correspondiente a los volúmenes de control W y P (polinomio P1 de la
figura). La otra utilizaría la información contenida en los volúmenes de control P y
E (polinomio P2 de la figura). Aunque ambas reconstrucciones son de segundo orden
y convergentes la una a la otra cuando se refine la malla, parece razonable elegir P2
como polinomio interpolante a la hora de evaluar fe.
Imponiendo que P2 debe reproducir fP y fE en los volúmenes de control VP y VE, se
llega al siguiente sistema de dos ecuaciones con dos incógnitas,
19
Habrá otros volúmenes al norte (N) y al oeste (W)
∆x
∆x
DR
AF
T
volúmenes finitos
Figure 6.: Reconstrucción de fe mediante la utilización de polinomios interpolantes
de grado uno. Dominio de P1 2 VW [ VP y dominio de P2 2 VP [ VE.
fP =
1
VP
Z
VP
(a2 + b2x) dV = a2 + b2
1
2VP
x2
�
�
e
w = a2 + b2xP ,
fE =
1
VE
Z
VE
(a2 + b2x) dV = a2 + b2
1
2VE
x2
�
�
e
w = a2 + b2xE ,
donde las incógnitas son a2 y b2 y VP = xe � xw. La solución del sistema anterior es,
a2 =
xEfP � xPfE
xE � xP
, b2 =
fE � fP
xE � xP
.
lo que conduce a al siguiente expresión para el polinomio interpolante P2,
P2(x) =
xEfP � xPfE
xE � xP
+
fE � fP
xE � xP
x . (10)
Finalmente, la evaluación de P2 en xe proporciona el valor buscado,
fe := P2(xe) = fefE + (1 � fe) fP , con fe =
xe � xP
xE � xP
. (11)
Esta expresión será utilizada a la largo de este capítulo a la hora de realizar las
reconstrucciones de orden dos de las variables computacionales en las caras de los
volúmenes de control. Dicho orden se debe a que se ha considerado un polinomio
interpolante de grado uno. Una expresión analítica del error se puede obtener con tal
de realizar dos desarrollos de Taylor entorno a P.
Demostración 1. Realizamos los siguientes desarrollos de Taylor,
fe = fP + (xe � xP) f0P +
1
2
(xe � xP)2 f00P + HOT
fE = fP + (xE � xP) f0P +
1
2
(xE � xP)2 f00P + HOT
Nota: HOT = high order terms. Despejando f0P de la segunda expresión,
f0P =
fE � fP
xE � xP
� 1
2
(xE � xP) f00P + HOT
20
Si utilizásemos grado 0 tendríamos polinomios de la forma p(x)=a
DR
AF
T
2.3 definición del método de los volúmenes finitos
y sustituyendo dicho resultado en la primera se obtiene,
fe ⇡ fP + (xe � xP)

fE � fP
xE � xP
� 1
2
(xE � xP) f00P
�
+
1
2
(xe � xP)2 f00P =
fefE + (1 � fe) fP +
1
2
⇥
(xe � xP)2 � (xe � xP)(xE � xP)
⇤
f00P =
fefE + (1 � fe) fP �
1
2
(xe � xP) (xE � xe) f00P .
Como se puede ver el error es de orden dos. En mallas equi-espaciadas, o en mallas en las
que el ratio entre dos celdas adyacentes es orden 1, Dxi/Dxi+1 ⇡ 1 el valor de dicho error es
⇡ (DxP/8)f00P.
⇤
Aunque esta reconstrucción suele ser apropiada en la mayoría de los casos, existen
variables en las que dicha reconstrucciónno es física. Supongamos por ejemplo que
deseamos reconstruir la difusividad térmica en la frontera de dos celdas (por ejemplo
en e) a partir de la expresión anterior. Dicha reconstrucción asume que la difusividad
térmica en e es igual a (1 � fe)kP cuando kE = 0. Esto no tiene sentido físico puesto
que la condición kE = 0 debería implicar que no se produce difusividad térmica en
el volumen de control VE. Es más, la difusividad en e debe estar acotada por la di-
fusividad más baja de entre las dos celdas adyacentes a la cara e. Una expresión que
reproduce la física viene dada por,
ke =
kEkP
fekE + (1 � fe) kP
.
Así por ejemplo, si kE ⌧ kP, entonces ke ⇡ KE/(1 � fe), como era de esperar. El
procedimiento de reconstrucción de las variables sobre la superficie de separación de
dos volúmenes de control se puede generalizar a alto orden. Suponiendo por ejemplo
un polinomio interpolante de grado 3, f(x) = a0 + a1x+ a2x2 + a3x3, se obtienen recon-
strucciones de cuarto orden. Para fijar los parámetros ai se necesitan cuatro volúmenes
de control, por ejemplo; W, P, E y EE. El último volumen de control se define como el
este (derecha) del volumen de control definido por E. Imponiendo que,
fW =
1
VW
Z
VW
f(x)dx ,
fP =
1
VP
Z
VP
f(x)dx ,
fE =
1
VE
Z
VE
f(x)dx ,
fEE =
1
VEE
Z
VEE
f(x)dx ,
se obtiene un sistema de cuatro ecuaciones con cuatro incógnitas (los parámetros ai).
Para una malla equi-espaciada estas incógnitas son,
21
DR
AF
T
volúmenes finitos
a0 = fW ,
a1 =
18fP � 11fW � 9fE + 2fEE
6Dx
,
a2 =
�5fP + 2fW + 4fE � fEE
2Dx2
,
a3 =
3fP � fW � 3fE + fEE
6Dx3
.
Sustituyendo los coeficientes anteriores en el polinomio interpolante y evaluando
éste sobre xe se obtiene la siguiente expresión,
fe =
27fP + 27fE � 3fW � 3fEE
48
.
El polinomio interpolante se puede utilizar para obtener una estimación de la derivada
en xe, para ello basta con derivar el polinomio con respecto a x y evaluarlo en xe,
✓
∂f
∂x
◆
e
=
�27fP + 27fE + fW � fEE
24Dx
.
Reconstrucción de los flujos difusivos
La reconstrucción del término S · (rf) f que aparece en los flujos difusivos depende
del tipo de malla. El caso más sencillo se tiene al considerar mallas ortogonales. En este
tipo de mallas el vector que une los centroides de dos celdas adyacentes es ortogonal a
la superficie de separación y además corta por el centroide de ésta, ver Figura 7-Derecha.
En este caso la siguiente expresión proporciona una aproximación de segundo orden,
S · (rf) f = |S|
fF � fP
|d| , (12)
donde VP y VF son dos volúmenes de control adyacentes. La expresión anterior se
puede obtener en el caso unidimensional a partir de la reconstrucción realizada en
(10). Derivando (10) con respecto a x y evaluando en e se obtiene (12) para el caso
particular en el que f = e y F = E. El error de truncación de la expresión (12) se puede
obtener a partir de dos desarrollos de Taylor evaluados en e.
fP = fe � (xe � xP) f0e +
1
2
(xe � xP)2 f00e �
1
6
(xe � xP)3 f000e + HOT
fE = fe + (xE � xe) f0e +
1
2
(xE � xe)2 f00e +
1
6
(xE � xe)3 f000e + HOT
Restando amabas expresiones y despejando f0e se obtiene,
f0e =
fE � fP
xE � xP
+
✓
(xe � xP)2 � (xE � xe)2
2 (xE � xP)
◆
f00e +
✓
(xe � xP)3 + (xE � xe)3
6 (xE � xP)
◆
f00e + HOT
Aunque formalmente el error no es orden dos, en mallas equi-espaciadas o en mallas
en las que el ratio de crecimiento entre dos celdas adyacentes es orden uno, se obtiene
orden dos conforme se refina la malla.
22
DR
AF
T
2.3 definición del método de los volúmenes finitos
Otra posibilidad a la hora de calcular el gradiente en la frontera entre dos volúmenes
de control consiste en calcular el gradiente en cada nodo y después interpolar sus
valores en la cara que éstos comparten. Es decir,
(rf) f = fx (rf)F + (1 � fx)rfP , (13)
donde fx =
f P
PF
=
|x f � xP|
|d| .
Figure 7.: Diferencia entre malla ortogonal (derecha) y no ortogonal (izquierda). Rep-
resentación de la cara f que comparten los nodos P y F. En la figura se
representan dos vectores. El vector que une los nodos P y F y el vector
normal a la superficie S.
Aunque los dos métodos presentados son de segundo orden, el segundo de ellos
utiliza dominios mayores, por lo que el error cometido es mayor.
En el caso de considerar mallas no ortogonales se debe proceder de otra manera a la
hora de calcular el termino S · (rf) f . Para ello se descompone S en una componente
ortogonal D y otra no ortogonal, k. Sustituyendo esta descomposición en S · (rf) f se
obtienen dos términos: el término ortogonal y la corrección no ortogonal,
S · (rf) f = |D|
fF � fP
|d| + k (rf) f , (14)
donde S = D + k. Se pueden considerar distintas descomposiciones de los vectores
anteriores: (a) D =
d · S
d · d d, (b) D =
d
|d| |S|, (c) D =
d
d · S |S|
2, siendo la más popular
la tercera de ellas dado que minimiza la contribución no ortogonal.
Una vez definido el vector D (o k), se realiza la corrección no-ortogonal, para lo que
se utiliza la Ecuación (13). Los términos de la derecha de la ecuación se pueden obtener
a partir de (9). La reconstrucción anterior puede generar problemas de estabilidad
dado que se consideran celdas computacionales grandes. La componente no ortogonal
puede generar problemas de acotación de las soluciones. En casos extremos dicha
contribución se debe acotar o incluso eliminar.
Otro problema que surge a la hora de evaluar los flujos difusivos se produce cuando
la malla tiene Skewness. Este problema surge cuando el punto de intersección del
vector que conecta dos nodos adyacentes con la superficie que separa los volúmenes
de control asociados a éstos no coincide con el centro de la cara. En las Figuras 8
se representa con Di la distancia entre el centroide de la cara y el punto de corte de
23
DR
AF
T
volúmenes finitos
ésta con el vector d, punto fi. Una medida del skewness viene dada por la siguiente
expresión,
y =
|Di|
|di|
. (15)
El skewness añade difusión numérica a la solución y reduce la precisión del método
llegado a generar soluciones no acotadas. El error debido al skewness puede reducirse
mediante la introducción de la desviación entre el punto fi y el punto f en el modelo.
Figure 8.: Izquierda: Malla skewness ortogonal Derecha: Malla skewness no-ortogonal.
El vector Di va de fi a f .
A priori, el procedimiento para calcular (rf) f sigue los siguientes tres pasos que se
deberán evaluar en orden inverso: 3, 2 y 1:
1. (rf) f = f̃xrfF +
�
1 � f̃x
�
rfP
2. (rf)P =
1
VP
 f
�
Sf f f
�
3. f f = f fi + Di · (rf) fi donde,
• f fi = fxfF + (1 � fx) fP
• (rf) fi = fxrfF + (1 � fx)rfP
siendo fx = fiP/(PN).
El problema que aparece se debe a que para evaluar el punto 2 es necesario saber el
valor de f f el cual necesitar conocer de antemano 2. Este problema se puede resolver
de forma iterativa. A continuación se exponen dos modelos.
Algoritmo 1 (Método iterativo 1). : Ver Figura 9-Izquierda
1. Aproximación inicial: (rf) fi =
fF � fP
d
2. Se calcula (rf)P =
1
VP
 f S f
⇣
f fi + Di · (rf) fi
⌘
3. Se calcula (rf) fi = fxrfF + (1 � fx)rfP ir a 2 (hasta convergencia)
24
DR
AF
T
2.3 definición del método de los volúmenes finitos
⇤
Algoritmo 2 (Método iterativo 2). : Ver Figura 9-Derecha
1. Aproximación inicial: f f =
fP + fF
2
2. Se calcula (rf)P =
1
VP
 f S f f f
3. Se calcula f f =
�
fP + DP f ·rfP
�
+
�
fF + DF f ·rfF
�
2
ir a 2 (hasta convergencia)
⇤
Una vez conocido (rf)P se puede evaluar (rf) f y por lo tanto evaluar la expresión
(14). En cuanto al coeficiente f̃x, éste se puede obtener a partir de la relación entre f f ,
fF y fP dada por: f f = f̃xfF +
�
1 � f̃x
�
fP.
Figure 9.: Figuras utilizadas en los modelos 1 (izquierda) y 2 (derecha).
Para evitar este tipo de problemas se han propuesto distintas mallas que simplifican
la reconstrucción de los flujos. Así por ejemplo en la Figura 10 se muestra la malla
de Voronoi la cual es un ejemplo de malla ortogonal no-estructurada. Dicha malla se
construye a partir de una triangularización de Delaunay, siendo los vértices de cadavolumen de control Voronoi los circuncentros de los triángulos de Delaunay que tienen
uno de sus vértices en el nodo v del Voronoi.
25
DR
AF
T
volúmenes finitos
Figure 10.: Izquierda: Delaunay, Derecha: Voronoi
Reconstrucción de los flujos convectivos
En la Figura 11 se representa el volumen de control utilizado en la evaluación de los
flujos convectivos S f · (ruf) f para el caso unidimensional.
Figure 11.: P es el centroide de la celda VP. Sus nodos adyacentes son W y E. La
longitud de la celda es dx y sus fronteras son w y e. Finalmente u es la
velocidad, la cual se considera positiva en la figura.
A priori el valor de f que aparece en el flujo se puede reconstruir utilizando la
reconstrucción de segundo orden vista con anterioridad: f f = fxfF + (1 � fx) fP, ver
Figura 12.
Sin embargo, esta reconstrucción puede generar soluciones no acotadas cuando u es
lo suficientemente grande como para romper la simetría del problema. En esos casos
se debe considerar un esquema upwinding del tipo,
f f =
⇢
fP si u > 0
fE si u < 0
Este esquema es tan sólo de primer orden, por lo que es un esquema difusivo. Se
puede obtener el valor de la difusividad numérica con tal de reemplaza el siguiente
desarrollo de Taylor,
fe = fP + (xe � xP) f0P + HOT
26
DR
AF
T
2.3 definición del método de los volúmenes finitos
Figure 12.: Interpolación sobre la cara f a partir de un esquema de segundo orden.
Figure 13.: Interpolación sobre la cara f a partir de un esquema upwinding de primer
orden.
27
DR
AF
T
volúmenes finitos
en el flujo convectivo, f c = rfu · n,
f c = rffP + rf (xe � xP) f0P + . . .
Como se puede ver el segundo término de la derecha representa un flujo de difusión
cuya difusividad numérica es rfDx/2. La extensión del esquema upwinding a más
dimensiones es directa con tal de repetir el esquema en cada una de las direcciones
espaciales por separado. Este procedimiento sin embargo genera resultados erróneos
cuando el campo de velocidad no es perpendicular a las caras que definen las celdas.
En este caso aparece una difusión numérica que difumina la solución. Este problema
se puede atenuar con tal de aumentar la resolución de la malla, aunque eso sí a costa
de aumentar el esfuerzo computacional.
El número de Peclet es número adimensional utilizado en procesos de transporte.
Este número se define como la relación entre el coeficiente de convección y el de di-
fusión y puede utilizarse para determinar cuándo es necesario aplicar un esquema
upwinding o cuando basta con considerar un esquema central. Dicho número se de-
fine como,
Pe = ru/(G/Dx) =
DxruG
k
,
donde G es el coeficiente de difusión.
2.3.3 Discretización temporal
Después de la discretización espacial se procede a la discretización temporal. Entre los
esquemas clásicos que se pueden considerar están: Euler explícito/implícito y Runge-
Kutta.
Integrando en el tiempo la Ecuación (3), se llega a,
Z t+Dt
t
dt
Z
VC
∂rf
∂t
dV +
Z t+Dt
t
dt
Z
VC
r · (rfu � Grf) dV =
Z t+Dt
t
dt
Z
VC
SfdV , (16)
Aplicando la regla de Leibniz2 en la primera integral y la discretización espacial
vista con anterioridad en los flujos y el término fuente se llega a,
Z t+Dt
t
dt
d
dt
(rf)P VP +
Z t+Dt
t
"
Â
f
S f · (ruf) f � Â
f
S f ·
�
rGfrf
�
f
#
dt =
Z t+Dt
t
SfVP .
(17)
La primera integral es directa,
VP (rf)P
�
�
t+Dt
t +
Z t+Dt
t
"
Â
f
S f · (ruf) f � Â
f
S f ·
�
rGfrf
�
f
#
dt =
Z t+Dt
t
SfVP (18)
En cuanto al resto de integrales se pueden considerar distintas cuadraturas:
2 Regla de Leibniz: Sea una función f (x, y) tal que fy(x, y) existe y es continua, entonces:
d
dy
✓
Z b(y)
a(y)
f (x, y)dx
◆
=
Z b(y)
a(y)
fy(x, y)dx + f (b(y), y)b0(y)� f (a(y), y)a0(y)
.
28
DR
AF
T
2.4 sistema algebraico de ecuaciones
• Euler explícita:
R t+Dt
t ydt = y
nDt (atrasada)
• Euler implícita:
R t+Dt
t ydt = y
n+1Dt (adelantada)
• Crank-Nicolson :
R t+Dt
t ydt =
1
2
Dt
�
yn+1 + yn
�
(trapezoidal)
Donde se ha asumido que yn = y(tn). Una forma general de escribir las expresiones
anteriores consiste en definir un parámetro q 2 [0, 1] tal que,
Z t+Dt
t
ydt = qDtyn+1 + (1 � q)Dtyn .
Reemplazando las integrales temporal por la expresión anterior se obtiene,
(rf)n+1 = (rf)n +
Dtq
VP
(ruf)n+1f · S f +
Dt(1 � q)
VP
(ruf)nf · S f
�Dtq
VP
(Grf)n+1f · S f �
Dt(1 � q)
VP
(Grf)nf · S f = DtqS
n+1
f + Dt(1 � q)S
n
f .
(19)
Se recupera el esquema Euler implícito para q = 1, el Euler explícito para q = 0 y el
Crank-Nicolson para q = 1/2.
2.4 sistema algebraico de ecuaciones
El sistema de ecuaciones (19) puede escribirse de forma simplificada como,
aiifn+1i = Â
i
fijfnj � Â
i,i 6=j
aijfn+1j + bi . (20)
Los coeficientes aij se deben a la derivada temporal, fij representan los coeficientes
que se derivan de los términos convectivos y difusivos y finalmente bi representa los
términos fuentes. El sistema es implícito si existe algún aij distinto de cero con i 6= j.
En caso contrario es explícito. El sistema es estacionario en el caso que de todos los
coeficientes aij = 0, 8i, j.
Así pues, en el caso no-estacionario e implícito el objetivo consiste en invertir la
matriz A con coeficientes aij para así obtener fn+1i . Para ello existen diversas técnicas
de inversión de matrices: (a) Métodos directos válidos para problemas pequeños (LU,
algoritmo de Thomas, etc.) (b) Métodos iterativos validos para problemas grandes (por
ejemplo, método de Gauss-Seidel). La convergencia de los métodos iterativos depen-
den de la estructura de la matriz. Dado que éstas se derivan del esquema numérico,
la convergencia de dichos algoritmos estará relacionada con la estabilidad de los es-
quemas numéricos considerados, para los parámetros de discretización elegidos: Dx y
Dt.
Entre las propiedades que debe cumplir la matriz se mencionaran las dos más impor-
tantes: (1) las matrices deben ser diagonal-dominantes, (2) los coeficientes del esquema
deben tener el signo apropiado. Respecto al primer punto, la definición de matrices
diagonal-dominantes viene dada por,
29
DR
AF
T
volúmenes finitos
Â
j
|aij|  |aii| , 8 i, j ,
9! i tal que Â
j
|aij| < |aii| , 8j .
Los términos de difusión crean matrices diagonal dominantes en mallas ortogonales.
Esto no es así en el caso general lo que hace que los métodos iterativos fallen. Para
evitar este problema se pueden plantear esquemas numéricos que hagan que las matri-
ces sean diagonal dominantes. Los esquemas que se muestran en este documento son
diagonal-dominantes, dándose la segunda condición en los nodos que limitan con el
contorno del problema.
En cuanto al signo de los coeficientes, es de esperar que un incremento de fni o f
n
j
produzcan un incremento de fn+1i . También es de esperar que si aumenta el valor de
fn+1j entonces f
n+1
i se incremente. Como se puede ver en (20) lo anterior implica que
los coeficientes aii y fij deben tener el mismo signo mientras que los coeficientes aii y
aij deben tener signo contrario. Como se verá en el Capítulo 3 esta condición permitirá
definir una condición CFL que garantice la convergencia del sistema de estudio.
En el caso de considerar problemas estacionarios el sistema (20) se reduce a,
� fiifi = Â
i,i 6=j
fijfnj + bi .
Definiendo f̃ij = � fij si i = j y f̃ij = fij para i 6= j, se tiene,
f̃iifi = Â
i,i 6=j
f̃ijfnj + bi .
El cambio de signo permite utilizar la notación seguida en el Capítulo 3 en los prob-
lemas estacionarios. En forma matricial el sistema queda como sigue para el caso
unidimensional,
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
f̃11 � f̃12
. . .
� f̃21 f̃22 � f̃23
. . .
. . . . . . . . . . . .
. . . . . . . . . . . . . . .
� f̃W f̃P � f̃E
. . . . . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . � f̃n,n�1 f̃nn
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
f1
fW
fP
fE
fn
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
=
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
bSbW
bP
bE
bn
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
Nuevamente se debe imponer que el sistema matricial definido por f̃ij sea diagonal-
dominante y que los coeficientes f̃ij tengan todos el mismo signo. Como se verá en
el Capítulo 3 esta condición permitirá definir un límite para el número de Peclet en
problemas convectivos-difusivos en los que se consideran esquemas centrales.
30
DR
AF
T
2.5 condiciones de contorno
2.5 condiciones de contorno
Las caras de los volúmenes de control que limitan con el contorno del dominio com-
putacional se deben tratar de forma especial. Sobre estas caras se pueden imponer
distintas condiciones de contorno matemáticas: Dirichlet, Neumann y Robin.
Recordando las expresiones de los flujos convectivos y difusivos, las condiciones de
contorno anteriores deberán definirse de tal manera que permitan reconstruir dichos
flujos sobre las caras de los volúmenes de control que limitan con el contorno físico del
problema:  f Ff f y  f
�
rGf
�
f S · (rf) f . Estas condiciones de contorno se especifican
en al Tabla 1
CC/Problema Convección Difusión
Dirichlet Fbfb S (rf)b = |S|
fb � fP
|dn|
Neumann fb = fP + dn · (rf)b S · (rf)b
Table 1.: Condiciones de contorno utilizadas en la reconstrucción de los flujos convec-
tivos y difusivos.
En cuanto a las condiciones de contorno de tipo Robin, estas son una combinación
de las condiciones de tipo Dirichlet y Neumann.
Por otro lado, cuando se computa numéricamente la solución de las Navier-Stokes
en un dominio no acotado es necesario utilizar condiciones de contorno artificiales
con el objetivo de limitar el dominio computacional. Estas condiciones de contorno
se denominan transparentes cuando la solución en el dominio acotado es igual a la
solución que se obtendría si se considerar el espacio completo. Por desgracia en mu-
chos casos de interés las condiciones de contorno existentes no son transparentes. Así
por ejemplo, las condiciones de contorno basadas en la linealización de la solución en
las proximidades de las fronteras (reales o artificiales) son apropiadas para problemas
de flujo de entrada y radiación, pero no para problemas con flujo saliente. Este prob-
lema se agrava cuando se consideran esquemas numéricos de alto orden en los que
la molécula computacional aumenta, dado que en estas circunstancias el tratamiento
incorrecto del contorno afecta una mayor extensión del dominio computacional. En
general, las discretizaciones que se suelen utilizar a la hora de resolver las condiciones
de contorno suele generar reflexiones que se propagan a lo largo del dominio com-
putacional destruyendo la estabilidad del método.
2.6 errores
Los errores que se comete en la resolución numérica pueden deberse al modelo o al
método usado en la resolución de éste. Estos últimos se definen como errores de
discretización y pueden ser de tres tipos.
1. Errores causados por la discretización del dominio de definición: insuficiente
resolución de la malla, no ortogonalidad de la malla y skewness.
31
DR
AF
T
volúmenes finitos
2. Errores causados por la discretización de las ecuaciones. En el caso de un método
de volúmenes finitos a segundo orden el error se puede representar como una
difusividad numérica.
3. Errores de resolución del sistema lineal de ecuaciones; procesos iterativos.
El error de las soluciones numéricas se obtiene a través de la comparación entre
la solución exacta y la solución numérica, siempre y cuando la primera sea conocida.
Existen dos formas de reducir el error del resultado. Aumentando la resolución o
incrementando el orden del método.
32
DR
AF
T
3
P R O B L E M A S V O L Ú M E N E S F I N I T O S
3.1 problemas unidimensionales
Con el objetivo de sentar las bases, en este capítulo se procede a resolver una serie
de problemas unidimensionales. Dichos problemas se han sacado del libro de Ver-
steeg and Malalasekera (1995): La ecuación de difusión-reacción estacionaria unidi-
mensional se describe en las Secciones 4.1, 4.2 y 4.3, los problemas de convección-
difusión estacionarios se describen en el Capítulo 5 mientras que la descripción del
problema de difusión-reacción no estacionario se encuentra en el Capítulo 8, (páginas
168-172). En cuanto al caso bidimensional, este se recoge en la Sección 4.4 (páginas
99-100). Estos problemas han sido extendidos en este capítulo para más condiciones
de contorno.
3.1.1 Ecuación de difusión reacción estacionaria
Se considera la ecuación de difusión reacción estacionaria unidimensional definida
sobre [A, B] 2 R,
Z
S
n · (krf) dS +
Z
VC
SfdV = 0 . (21)
Por integración directa se tiene que,
k
df
dx
�
�
�
�
b
a
+
Z b
a
Sfdx = 0 . (22)
Los valores a y b hacen referencia a los límite izquierdo y derecho del volumen de
control de integración.
El primer paso consiste en dividir el dominio en un conjunto discreto de volúmenes
de control o celdas. Por ejemplo, se puede posicionar una serie de nodos en el interior
del dominio (círculos negros de la Figura 14) y situar las fronteras de las celdas en
el punto medio entre los nodos (lineas rectas). Más concretamente, dado un dominio
[A, B] éste se discretiza con una malla {xi; i = 1, . . . N}, siendo N el número de nodos.
Entorno a cada nodo se define una celda Ci de la siguiente manera:
Ci =
✓
xi + xi�1
2
,
xi + xi+1
2
◆
.
Como puede verse, la longitud en cada celda Ci es Dxi = (xi+1 � xi�1)/2. Normal-
mente se usan mallas equi-espaciadas tal que Dxi = Dx, 8i = 1, . . . , N.
33
DR
AF
T
problemas volúmenes finitos
Figure 14.: Discretización del dominio AB en celdas o volúmenes de control en el caso
unidimensional. El nodo P es el centro del volumen de control de estudio
VP. Los nodos E y W son los nodos de las celdas adyacentes. La longitud
de la celda es dx y sus fronteras son w y e.
En lugar de utilizar la notación anterior con índices i, i � 1 e i + 1, se utilizará la
notación descrita en la Figura 14. El siguiente paso consiste en la discretización de las
integrales de volumen: SfVP = Su + SpfP. Sustituyendo la expresión anterior en (22)
se llega a,
✓
k
df
dx
◆
e
�
✓
k
df
dx
◆
w
+ SPdxwe = 0 , (23)
donde los dos primeros términos están evaluados en los contornos de la celda. Por
lo tanto, se necesita el valor del coeficiente de difusión y la derivada de f sobre el
contorno de VP. Estos valores no se conocen a priori dado que los valores de f sólo se
conocen sobre los nodos. Tal y como se vio, a segundo orden de aproximación se tiene
que,
✓
df
dx
◆
e
=
fE � fP
dxPE
+ O(Dx2) ,
✓
df
dx
◆
w
=
fP � fW
dxWP
+ O(Dx2) .
y
kw = fwkW + (1 � fw)kP , fw =
wP
WP
,
ke = fekE + (1 � fe)kP , fe =
eP
EP
,
Sustituyendo las expresiones anteriores en la Ecuación (23),
ke
fE � fP
dxPE
� kw
fP � fW
dxWP
+ Sdxwe = 0 . (24)
y reagrupando, se obtiene que
✓
ke
dxPE
+
kw
dxWP
� Sp
◆
fP =
✓
kw
dxWP
◆
fW +
✓
ke
dxPE
◆
fE + Sdxwe . (25)
34
DR
AF
T
3.1 problemas unidimensionales
Se puede utilizar la Ecuación (24) para obtener las expresiones para los volúmenes de
control definidos sobre el contorno con tal de sustituir W por A en el caso de considerar
el contorno izquierdo y E por B en el caso de considerar el contorno derecho:
Figure 15.: Condiciones de contorno Dirichlet en los contornos izquierdo (A) y dere-
cho (B) en problemas unidimensionales.
Figure 16.: Condiciones de contorno Neumann en los contornos izquierdo (A) y dere-
cho (B) en problemas unidimensionales.
• Dirichlet en A (Ver Figura 15-Izquierda)
ke
fE � fP
dxPE
� kw
fP � fA
dxAP
+ SdxAe = 0 . (26)
• Neumann en A (Ver Figura 16-Izquierda)
ke
fE � fP
dxPE
� qA + SdxAe = 0 . (27)
• Dirichlet en B (Ver Figura 15-Derecha)
kB
fB � fP
dxPB
� kw
fP � fW
dxWP
+ SdxwB = 0 . (28)
• Neumann en B (Ver Figura 16-Derecha)
qB � kw
fP � fW
dxWP
+ SdxwB = 0 . (29)
Las expresiones anteriores se puede escribir de forma más compacta como:
aPfP = aWfW + aEfE + Su .
Los resultados se recogen en la Tabla 2. Como se puede comprobar en dicha tabla
el término fuente tiene la forma,SVP = Su + SPfb donde b hace referencia al contorno
físico.
35
DR
AF
T
problemas volúmenes finitos
Volumen control aW aE SP Su
Interior
kw
dxWP
ke
dxPE
0 Sdxwe
Dirichlet A 0
ke
dxPE
� kA
dxPA
SdxAe +
kA
dxPA
fA
Neumann A 0
ke
dxPE
0 SdxAe � qA
Dirichlet B
kw
dxWP
0 � kB
dxPB
SdxwB +
kB
dxPB
fB
Neumann B
kw
dxWP
0 0 SdxwB + qB
Table 2.: Coeficientes ecuación difusión - reacción estacionario unidimensional, donde
aP = aE + aW � SP.
Independientemente de las condiciones de contorno consideras la Tabla 2 define un
sistema matricial de la forma descrita en (30).
0
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
@
a1 �ae1 0 . . .
�aw2 a2 �ae2 0 . . .
0 �aw3 a3 �ae3 0 . . .
· · ·
. . . 0 �awi ai �aei 0 . . .
· · ·
. . . 0 �awN�1 aN�1 �aeN�1
. . . 0 �awN aN
1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A
·
0
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
@
f1
f2
f3
fi
fN�1
fN
1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A
=
0
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
@
s1
s2
s3
si
sN�1
sN
1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A
(30)
Durante la construcción del sistema se debe tener especial cuidado en los nodos que
limitan con la frontera en cuyo caso se debe considerar las condiciones de contorno a
imponer.
La solución del sistema anterior implica invertir la matriz para obtener fi. Tal y
como se dijo, existen diversas técnicas de resolución, aunque dado que la matriz es
tridiagonal se podría utilizar el algoritmo de Thomas.
Problem 1 (Difusión del calor sobre una barra unidimensional). En este primer ejemplo
se estudia el problema de propagación del calor por conducción sobre una barra unidimensional.
d
dx
✓
k
dT
dx
◆
= 0 .
donde T es el campo de temperaturas y k es la conductividad térmica del material. Se supone
que, dx = 0.1 m, A = 0.001 m2, kA/dx = 100 W/K, TA = 100 Ky TB = 500 K.
⇤
36
DR
AF
T
3.1 problemas unidimensionales
Nota: En las expresiones anteriores el área A se simplifico. Este área se recupera
en el problema con el objetivo de que los números que aparecen en la matriz sean
razonables.
Figure 17.: Discretización problema 1. Discretización del dominio en 5 volúmenes de
control.
En primer lugar se divide el dominio computacional en 5 volúmenes de control, ver
Figura 17. Los nodos 2, 3 y 4 tiene nodos a su izquierda y derecha, y en ellos se
puede aplicar directamente la ecuación discretizada (volúmenes de control internos de
la Tabla 2). Por otro lado, en los contornos se consideran condiciones de tipo Dirichlet
(ver Tabla 2). Asumiendo que el medio es homogéneo, ke = kw = k, para todos los
nodos y que los nodos están distribuidos equi-espaciadamente con separación dx, se
llega a:
Volumen control aW aE SP Su
Interior
k
dx
k
dx
0 0
Dirichlet A 0
k
dx
� 2k
dx
2k
dx
TA
Dirichlet B
k
dx
0 � 2k
dx
2k
dx
TB
Table 3.: Esquema numérico a utilizar en el problema 1, donde aW = aE =
k
dx
y
aP = aW + aE
La presencia de las condiciones de contorno implican la modificación del coeficiente
que multiplica TP, así como añaden un término en el lado derecho del sistema matricial,
Su. Sustituyendo el valor de los parámetros se llega a:
2
6
6
6
6
6
6
6
6
6
6
6
6
4
300 �100 0 0 0
�100 200 �100 0 0
0 �100 200 �100 0
0 0 �100 200 �100
0 0 0 �100 300
3
7
7
7
7
7
7
7
7
7
7
7
7
5
·
2
6
6
6
6
6
6
6
6
6
6
6
6
4
T1
T2
T3
T4
T5
3
7
7
7
7
7
7
7
7
7
7
7
7
5
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
200TA
0
0
0
0
200TB
3
7
7
7
7
7
7
7
7
7
7
7
7
5
37
DR
AF
T
problemas volúmenes finitos
Como se puede ver, la matriz anterior es diagonal dominante. Para TA = 100 y
TB = 500 la solución se puede obtener por eliminación Gausiana,
T =
2
6
6
6
6
6
4
140
220
300
380
460
3
7
7
7
7
7
5
mientras que la solución analítica viene dada por T = 800x + 100. Como se puede ver
en la Figura 18, existe una buena correspondencia entre ambos resultados.
0 0.1 0.2 0.3 0.4 0.5100
150
200
250
300
350
400
450
500
X
T
Figure 18.: Problema 1. Analítico (líneas) vs. Numérico (círculos).
Problem 2 (Difusión del calor sobre barra con Q uniforme). En este caso se considera un
problema análogo al anterior pero en el que existe un término fuente de calor constante sobre el
dominio.
d
dx
✓
k
dT
dx
◆
= q .
T es el campo de temperaturas, k es la conductividad térmica y q la fuente de calor. Se supone
que, dx = 4 ⇥ 10�3 m, A = 1 m2, k = 0.5 W/(mK), q = 1 ⇥ 106 W/m3, TA = 100 K y
TB = 200 K. ⇤
Utilizando de nuevo los resultados presentados en la Tabla 2 para el caso en el que k
y dx son constantes se obtiene el resultado que se muestra en la Tabla 4. Sustituyendo
el valor de los parámetros en los coeficientes se llega al siguiente sistema matricial,
38
DR
AF
T
3.1 problemas unidimensionales
Volúmenes control aW aE SP Su
Interior
k
dx
k
dx
0 qdx
Dirichlet A 0
k
dx
� 2k
dx
qdx +
2k
dx
TA
Dirichlet B
k
dx
0 � 2k
dx
qdx +
2k
dx
TB
Table 4.: Esquema numérico a utilizar en el problema 2, donde aW = aE =
k
dx
y
aP = aW + aE
2
6
6
6
6
6
6
6
6
6
6
6
6
4
375 �125 0 0 0
�125 250 �125 0 0
0 �125 250 �125 0
0 0 �125 250 �125
0 0 0 �125 375
3
7
7
7
7
7
7
7
7
7
7
7
7
5
·
2
6
6
6
6
6
6
6
6
6
6
6
6
4
T1
T2
T3
T4
T5
3
7
7
7
7
7
7
7
7
7
7
7
7
5
=
2
6
6
6
6
6
6
6
6
6
6
6
6
4
29000
4000
4000
4000
4000
54000
3
7
7
7
7
7
7
7
7
7
7
7
7
5
Nuevamente, la solución se puede obtener por eliminación Gausiana,
T =
2
6
6
6
6
6
4
150
218
254
258
230
3
7
7
7
7
7
5
Mientras que la solución analítica viene dada por,
T =
✓
TB � TA
L
+
q
2k
(L � x)
◆
x + TA .
Un estudio de convergencia muestra que la aproximación numérica es de segundo
orden, ver Figura 19.
3.1.2 Ecuación de convección difusión reacción estacionaria
Tal y como se verá en esta sección, los esquemas centrales generan soluciones no aco-
tadas en problemas con términos convectivos cuando el número de Peclet es suficien-
temente grande. Esto se debe a que la solución numérica trata por igual ambas direc-
ciones, cuando en verdad el término convectivo marca una dirección de avance de la
información. En esta sección se estudiaran este tipo de problemas utilizando para ello
una discretización central y otra upwinding en la que se tiene en cuenta la dirección
del flujo.
39
DR
AF
T
problemas volúmenes finitos
-3
-2.5
-2
-1.5
-1
-0.5
-1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
Lo
g 1
0(
Er
ro
r)
Log10(∆ z)
fitting
y = m x + n, m=2.000000, n=0.795880
Figure 19.: Problema 2. Convergencia, Upwinding
Partiendo de la ecuación de convección-difusión estacionaria unidimensional definida
sobre el dominio [A, B]
Z
S
n · (rfu) dS =
Z
S
n · (krf) dS +
Z
VC
SfdV ,
y evaluándola sobre los volúmenes de control definidos en la Figura 20, se llega a la
siguiente ecuación semi-discretizada,
(ruf)e � (ruf)w =
✓
k
df
dx
◆
e
�
✓
k
df
dx
◆
w
+ Sfdxwe .
Figure 20.: P es el centroide de la celda VP. Sus nodos adyacentes son W y E. La
longitud de la celda es dx y sus fronteras son w y e. u es la velocidad, la
cual se considera positiva.
Los términos difusivos se puede discretizar tal y como se hizo en la Sección 2.3.2,
✓
df
dx
◆
e
=
fE � fP
dxPE
+ O(Dx2) ,
✓
df
dx
◆
w
=
fP � fW
dxWP
+ O(Dx2) .
En cuanto a los términos convectivos, se van a estudiar dos alternativas: (1) central
y (2) upwinding.
1. Esquema central:
40
DR
AF
T
3.1 problemas unidimensionales
Partiendo de,
Fefe � Fwfw =
ke
dxPE
(fE � fP)�
kw
dxWP
(fP � fW) + Sfdxwe ,
donde Fe = Fw := (ru)w es el flujo de masa. Se reconstruyen los valores de fe y
fw a partir de,
fw = fwfW + (1 � fw)fP , fw =
wP
WP
,
fe = fefE + (1 � fe)fP , fe =
eP
EP
.
Sustituyendo estas expresiones en la ecuación anterior,
Fe fefE + Fe (1 � fe) fP � FwfwfW � Fw (1 � fw) fP =
ke
dxPE
(fE � fP)�
kw
dxWP
(fP � fW) + Sdxwe ,
y reagrupando, se obtiene,
✓✓
kw
dxWP
+ fwFw
◆
+
✓
ke
dxPE
� feFe
◆
+ (Fe � Fw)
◆
fP =
✓
kw
dxWP
+ fwFw
◆
fW +
✓
ke
dxPE
� feFe
◆
fE + Sdxwe .
(31)
2. Esquema upwinding:
En el esquema anterior el valor de f sobre w se obtiene pesando los valores de
fP y fW de igual manera. Sin embargo, tal y como se puede ver en la Figura
20, el hecho de que el flujo vaya de W a P implica que el valor de fW debería
tener mayor peso que el de fP a la hora de calcular el valor de f sobre w. A
continuación se presenta un modelo sencillo que tiene en cuenta la dirección del
flujo a la hora de estimar el flujo convectivo sobre las fronteras de la celda P; e y
w.
si u > 0: Se considera que fw = fW y fe = fP
si u < 0: Se considera que fw = fP y fe = fE
El caso de u = 0, se corresponde con la discretización central. Por lo tanto, en
función de la dirección de u se tienen dos esquemas.
a) Si u > 0:
✓✓
kw
dxW P
+ Fw
◆
+
ke
dxPE
+ (Fe � Fw )
◆
fP =
✓
kw
dxW P
+ Fw
◆
fW +
✓
ke
dxPE
◆
fE + Sdxwe
b) Si u < 0:
✓
kw
dxW P
+
✓
ke
dxPE
� Fe
◆
+ (Fe � Fw )
◆
fP =
✓
kw
dxW P
◆
fW +
✓
ke
dxPE
� Fe
◆
fE + Sdxwe
41
DR
AF
T
problemas volúmenes finitos
En cuanto al tratamiento del contorno sólo se va a considerar condiciones de tipo
Dirichlet para el caso u > 0.
• Dirichlet en A:
Fe fP � FA fA = ke
fE � fP
dxPE
� k A
fP � fA
dx AP
+ Sdx Ae = 0 .
Reagrupando,
✓
ke
dxPE
+
k A
dx AP
+ Fe
◆
fP =
ke
dxPE
fE + 0fW +
✓
FA +
k A
dx AP
◆
fA + Sdx Ae
• Dirichlet en B:
FB fP � Fw fW = k B
fB � fP
dxPB
� kw
fP � fW
dxW P
+ SdxwB = 0 .
Reagrupando,
✓
kw
dxW P
+ Fw +
k B
dxPB
+ FB � Fw
◆
fP =
✓
kw
dxW P
+ Fw
◆
fW +
k B
dxPB
fB + SdxwB
La Tabla 5 recoge los resultados obtenidos.
aW aE S p Su
Interior
u > 0
kw
dxW P
+ Fw
ke
dxPE
Fw � Fe Sdxwe
u < 0
kw
dxW P
ke
dxPE
� Fe Fw � Fe Sdxwe
Dirichlet 0
ke
dxPE
�
✓
k A
dx AP
+ Fe
◆
Sdx Ae+
u > 0, A
✓
FA +
k A
dx AP
◆
fA
Dirichlet
kw
dxW P
+ Fw 0 �
k B
dxPB
SdxwB+
u > 0, B �(FB � Fw )
k B
dxPB
fB
Table 5.: Coeficientes del esquema upwinding, donde aP = aW + aE � S p
El error que se comete en los esquemas centrales es de segundo orden mientras que
los esquemas upwinding es de primer orden.
42
DR
AF
T
3.1 problemas unidimensionales
VC aW aE S p Su
Interior D + F/2 D � F/2 0 0
Dirichlet A 0 D � F/2 �(2D + F) (2D + F)fA
Dirichlet B D + F/2 0 �(2D + F) (2D � F)fB
Table 6.: Esquema central, donde D = k/dx y S = 0.
VC aW aE S p Su
Interior D + F D 0 0
Dirichlet A 0 D �(2D + F) (2D + F)fA
Dirichlet B D + F 0 �2D 2DfB
Table 7.: Esquema upwinding, donde D = k/dx, S = 0 y u > 0.
Con el objetivo de poder comparar más fácilmente los esquemas central y upwinding
se asume que ke = kw = k, dx es constante, Fe = Fw = F, D = k/dx y F/ D es el
número de Peclet (ratio entre convección y difusión), véase Tablas 6 y 7.
Tal y como puede verse en dichas tablas, los coeficientes en el esquema upwinding
son siempre positivos mientras que el esquema central el coeficiente aE puede ser
negativo cuando D < F/2. La condición que se debe cumplir para que esto no
suceda se traduce en que F/D < 2, es decir, el número de Peclet debe ser menor que
dos.
Problem 3 (Convección-difusión estacionaria). Se considera la variación de f sobre un
dominio L debido a un proceso de difusión y otro de convección forzada debida al arrastre que
produce un campo de velocidad uniforme u.
Asumiendo que fA = 1 y fB = 0, se tiene que la solución viene dada por,
f � fA
fB � fA
=
exp (rux/k) � 1
exp (ru L/k) � 1 .
⇤
Al igual que en el caso anterior el dominio computacional se divide en 5 celdas.
En cuanto al valor de u, se consideraran 3 posibles configuraciones:
• Caso 1: u = 0.1m/s, r = 1, k = 0.1 y dx = 0.2, Ncel das = 10
• Caso 2: u = 2.5m/s, r = 1, k = 0.1 y dx = 0.2, Ncel das = 10
• Caso 3: u = 2.5m/s, r = 1, k = 0.1 y dx = 0.05, Ncel das = 100
Los números de Peclet para cada uno de ellas son respectivamente (ru/(k/dx )):
0.2, 5, 1.25. Los resultados se muestran y comentan en las Figuras 21, 22 y 23.
43
DR
AF
T
problemas volúmenes finitos
0 0.2 0.4 0.6 0.8 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
X
T
0 0.1 0.2 0.3 0.4 0.5 0.6 0.70.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
1.05
X
T
Figure 21.: Caso 1: Pe = 0.2. Izquierda: Central, Derecha Upwinding. Como se pueden
ver ambos esquemas reproducen la solución analítica aunque el esquema
upwinding es ligeramente difusivo.
0 0.2 0.4 0.6 0.8 1
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
X
T
0 0.2 0.4 0.6 0.8 10.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
1.05
X
T
Figure 22.: Caso 2: Pe = 5.0. Izquierda: Central, Derecha Upwinding. El esquema
central es inestable mientras que el esquema upwinding es difusivo.
0 0.2 0.4 0.6 0.8 10.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
X
T
0 0.2 0.4 0.6 0.8 10.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
1.05
X
T
Figure 23.: Caso 3: Pe = 1.5. Izquierda: Central, Derecha Upwinding. Aumentado el
número de nodos se vuelve a recuperar la solución analítica, viéndose el
comportamiento difusivo en el esquema upwinding.
44
DR
AF
T
3.1 problemas unidimensionales
-2
-1.8
-1.6
-1.4
-1.2
-1
-0.8
-0.6
-2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8
Lo
g 1
0(
Er
ro
r)
Log10(∆ z)
fitting
y = m x + n, m=0.700664, n=-0.047444
Figure 24.: Convergencia caso upwinding
3.1.3 Ecuación de difusión reacción dependiente del tiempo
Problem 4. Convección difusión no estacionaria: Problema no-estacionario de difusión de calor
sobre el dominio [A , B ] donde TA y qB están fijos. ⇤
Se parte de la ecuación,
rc
∂T
∂ t
=
∂
∂x
✓
k
∂T
∂x
◆
+ S .
Utilizando la Ecuación (19) se llega a,
rcVP T n+1p = rcVP T np + q
✓
k
∂T
∂x
�
�
e � k
∂T
∂x
�
�
w
◆n+1
AD t+
(1 + q )
✓
k
∂T
∂x
�
�
e � k
∂T
∂x
�
�
w
◆n
AD t + SVP D t .
Reemplazando los flujos en w y e, dividiendo por AD t y reagrupando se obtiene,
✓
rcDx
D t
+ q (aE + aW )
◆
T n+1P =
✓
rcDx
D t
� (1 � q ) (aE + aW )
◆
T nP +
aE
⇣
q T n+1E + (1 � q )T nE
⌘
+ aW
⇣
q T n+1W + (1 � q )T nW
⌘
+ SDx ,
donde aE =
ke
dxPE
y aW =
kw
dxW P
.
Tal y como se dijo en la Sección 2.4, los coeficientes deben tener el mismo signo por
razones de convergencia. Asumiendo que dicho signo es positivo, se puede ver que el
coeficiente que multiplica T nP puede ser negativo. Imponiendo que dicho coeficiente
es positivo se tiene que,
rcDx
D t
> (1 � q ) (aE + aW )
En el caso de considera Euler explícito se llega a que D t < rcDx2 /(2k), mientras
que en el caso de considerar Crank-Nicolson D t < rcDx2 /k. Como se puede ver,
existe un cota para el paso de tiempo que cumple la condición de positividad de los
coeficientes. Esta condición puede verse como una condición de CFL.
45
DR
AF
T
problemas volúmenes finitos
Si se considera un esquema Euler implícito, los coeficientes siempre son positivos
por lo que el esquema numérico es incondicionalmente estable. Esto sin embargo no
significa que cualquier paso de tiempo sea válido puesto que el esquema es tan sólo
de primer orden y por lo tanto los errores que se comenten cuando D t es grande son
grandes.
3.2 problemas bidimensionales
El proceso seguido en la discretización de las ecuaciones unidimensionales puede gen-
eralizarse fácilmente al caso bidimensional. Para mostrar el procedimiento a seguir se
utilizará la ecuación de difusión bidimensional,
∂
∂x
✓
k
∂f
∂x
◆
+
∂
∂y
✓
k
∂f
∂y
◆
+ S = 0 . (32)
En la Figura 25 se muestra una porción de la malla bidimensional. Ahora además
de los nodos W y E, definidos a lo largo del eje x, se deben considerar los nodos N
(norte) y S (sur), a lo largo del eje y. La zona gris representa la celda computacional
asociada al nodo P. Las fronteras de esta celda con las celdas adyacentes N , S, E y W
reciben respectivamente los nombres n, s, e y w. Aunque la formulación del problema
se hará de forma general, al final los resultados se particularizarán para el caso en el
queDx y Dy son constantes para todas las celdas.
Figure 25.: Discretización utilizada en el problema bidimensional en el caso de consid-
erar una malla cartesiana.
La integración formal de la Ecuación (32) viene dada por,
ke Le
✓
∂f
∂x
◆
e
� kw Lw
✓
∂f
∂x
◆
w
+ kn Ln
✓
∂f
∂y
◆
n
� ks Ls
✓
∂f
∂y
◆
s
+ S A = 0 .
Como se puede ver los dos primeros términos tienen la misma estructura que la que
se obtuvo en el caso unidimensional para el problema de propagación de calor a lo
largo del eje x. Lo mismo se puede argumentar de los dos términos siguientes para
46
DR
AF
T
3.3 mallas no estructuradas
el caso de considerar un problema de propagación del calor a lo largo del eje y. La
diferencia con respecto al problema unidimensional es que ahora aparece unas nuevos
parámetros Le , Lw , Ln y Ls que representan la longitud de las caras de las celdas.
Finalmente, S representa el área de la celda P.
Es necesario proporcionar una expresión de los flujos de la variable f a través de
las caras de la celda de control. Para ello se utilizara la misma formulación (válida a
orden 2) que la considerada en el caso unidimensional:
flujo a través de la cara oeste
kw Lw
∂f
∂x
�
�
�
�
w
= kw Lw
fP � fW
dxW P
,
flujo a través de la cara este
ke Le
∂f
∂x
�
�
�
�
e
= ke Le
fP � fE
dxPE
,
flujo a través de la cara sur
ks Ls
∂f
∂x
�
�
�
�
s
= ks Ls
fP � fS
dxSP
,
flujo a través de la cara norte
kn Ln
∂f
∂x
�
�
�
�
n
= kn Ln
fP � fN
dxP N
.
Asumiendo que la distribución de puntos es equi-espaciada en ambas direcciones se
tiene que Le = Lw = Dx y Ln = Ls = Dy. En ese caso los valores de k sobre las caras
se pueden obtener a partir del valor medio. Por ejemplo, kw = kW +k P2 .
Sustituyendo las expresiones anteriores en la ecuación semi-discretizada y reagru-
pando se obtiene que,
✓
ke Le
dxPE
+
kw Lw
dxW P
+
kn Ln
dxP N
+
ks Ls
dxSP
� S p
◆
fP =
✓
ke Le
dxPE
◆
fE +
✓
kw Lw
dxW P
◆
fW +
✓
kn Ln
dxP N
◆
fN +
✓
ks Ls
dxSP
◆
fS + Su .
Esta expresión se puede escribir de forma más compacta como,
aP fP = aE fE + aW fW + aS fS + a N fN + Su .
Este sistema de ecuaciones nos conduce a un sistema penta-diagonal, el cual puede
resolverse de manera eficiente mediante algoritmos específicos para este tipo de matri-
ces.
3.3 mallas no estructuradas
Una de las principales ventajas del método de volúmenes finitos en su fácil imple-
mentación sobre mallas no estructuradas, las cuales suelen ser más convenientes a la
47
DR
AF
T
problemas volúmenes finitos
Figure 26.: Volumen de control en el tridimensional cartesiano.
hora de resolver problemas definidos sobre dominios complejos. En la Figura 27 se
muestran dos ejemplos, en el primero se usa una triangularización de Delaunay en un
caso bidimensional mientras que en el segundo se muestra la proyección de la malla
no estructurada tridimensional sobre el ala y fuselaje de un avión. Dos de las comple-
jidades que implica esta generalización tiene que ver con la identificación de las celdas
vecinas a cada celda a la hora de realizar el cálculo de las derivadas así como proyectar
dichas derivadas a lo largo de la normal a las caras de las celdas.
Debido a que las mallas no estructuradas no presentan un patrón de conectividad
predeterminado, es necesario identificar los vecinos que tiene cada celda o volumen
de control antes de comenzar a realizar el cálculo. Así el volumen de control esta
delimitado con un número de caras planas arbitrario, tal que cada una de ellas o está
compartida con un único vecino o pertenece al contorno exterior. El vector que define
cada cara, S, se construye de tal manera que apunta hacia el exterior del volumen de
control, es normal a la cara y su longitud es igual al área de la superficie. De la misma
manera los vectores que definen el contorno del problema apuntan hacia fuera del
dominio computacional.
48
DR
AF
T
3.4 carácter conservativo del esquema numérico
X
Y
0 0.5 1
0
0.5
1
1.5
Figure 27.: Izquierda: Malla Delaunay no estructurada generada para un pentágono.
Derecha: Proyección de una malla tridimensional no estructurada definida
entorno a un avión sobre su ala y fuselaje.
3.4 carácter conservativo del esquema numérico
Se considera un problema estacionario de difusión unidimensional en el que se aplican
dos fuentes de calor en ambos extremos del dominio, ver Figura 28.
Dado que el esquema es conservativo, cabría esperar que la suma de los balances
de energía en todos y cada uno de los volúmenes de control sea igual al balance de
energía total, Dq = qB � q A . Es decir, Â5i=1 k
df
dx
�
�
ei
wi
= qB � q A . Descomponiendo el
sumatorio y discretizando los flujos se llega a,
Figure 28.: Problema utilizado en la demostración del carácter conservativo de los
esquemas basados en volúmenes finitos.
5
Â
i=1
k
df
dx
�
�
ei
wi
= k
df
dx
�
�
e1
A + k
df
dx
�
�
e2
w2
+ k
df
dx
�
�
e3
w3
+ k
df
dx
�
�
e4
w4
+ k
df
dx
�
�
B
w5
=

ke1
fe2 � fe1
Dx
� q A
�
+

ke2
fe3 � fe2
Dx
� kw2
f2 � f1
Dx
�
+

ke3
fe4 � fe3
Dx
� kw3
f3 � f2
Dx
�
+

ke4
fe5 � fe4
Dx
� kw4
f4 � f3
Dx
�
+

qB � kw5
f5 � f4
Dx
�
.
49
DR
AF
T
problemas volúmenes finitos
Teniendo en cuenta que ke1 = kw2 , ke2 = kw3 ,... se llega a que los flujos en las caras
internas de los volúmenes de control cancelan a pares. Es decir, el flujo que entra en un
volumen de control a través de una cara es igual al flujo que sale a través de esa misma
cara en el volumen de control adyacente. Esto demuestra que el esquema numérico es
conservativo, es decir, Â5i=1 k
df
dx
�
�
ei
wi
= qB � q A .
Nota: Cuando se escribe un código se aprovecha esta circunstancia y tan sólo se
calcula uno de los flujos pues el otro es igual pero con el signo cambiado.
50
DR
AF
TPart IIPA RT E 2
DR
AF
T
DR
AF
T
4
I N T E G R A D O R E S D E R I E M A N N
4.1 sistemas hiperbólicos
Los sistemas hiperbólicos aparecen en el movimiento de ondas o transporte de sustan-
cias por convección, y suelen generar discontinuidades como por ejemplo ondas de
choque. El cálculo del flujo a través de dichas ondas resulta difícil debido a los fuertes
gradientes existentes. Para ello existen dos técnicas de resolución: (a) shock capturing
y (b) shock fitting. En la primera de ellas las ondas de choque surgen de la apli-
cación del esquema numérico, mientras que en el segundo caso las ondas de choque
se resuelven explícitamente utilizando para ello las relaciones que modelan el choque
(por ejemplo, relaciones de Rankine-Hugoniot). Aunque los métodos shock capturing
suelen ser difusivos, tienden a generar oscilaciones numéricas (fenómeno de Gibbs) y
pueden llegar a predecir velocidades erróneas de propagación de los choques, debido
a su fácil implementación éstos suelen utilizarse más. En éste capítulo se estudiaran
estos métodos.
La razón por la que se tienen que estudiar estos métodos numéricos se debe a que los
modelos upwinding descritos en los capítulos anteriores dejan de ser válidos cuando
los flujos convectivos son no lineales. Por otro lado, a lo largo de este capítulo se
estudiaran problemas unidimensionales. La generalización al caso bi- y tridimensional
se verá más adelante. En estos últimos casos se podrá utilizar la teoría desarrollada
para problemas unidimensionales.
El contenido de esta sección se encuentra en páginas 41-46 de Toro (2009). La parte
correspondiente a la resolución de sistema lineales, Secciones 4.2 y 4.3, se encuentran
en las páginas 47-86 de Toro (2009). La reconstrucción de los flujos de Godunov se
describe en más detalle en las Secciones 6.1 y 6.2 de Toro (2009) mientras que la aproxi-
mación del flujo de Roe se encuentra en la Sección 11.1. En cuanto a la implementación
del modelo numérico para las ecuaciones de Euler, éste se puede encontrar en el Capí-
tulo 11 (Sección 11.2.2) de Toro (2009). Finalmente la información correspondiente a las
reconstrucciones de alto orden se han sacado de distintos artículos de investigación.
4.1.1 EDPs de primerorden
Se considera un sistema de m ecuaciones con m incógnitas (ui) que dependen del
espacio x y el tiempo t,
U t + A Ux + B = 0 , (33)
donde,
53
DR
AF
T
integradores de riemann
U =
2
6
6
6
4
u1
u2
...
um
3
7
7
7
5
, B =
2
6
6
6
4
b1
b2
...
bm
3
7
7
7
5
, A =
2
6
6
6
4
a11 a12 . . . a1m
a21 a22 . . . a2m
...
...
...
...
am1 am2 . . . amm
3
7
7
7
5
,
El sistema anterior se puede clasificar en:
• Sistema lineal con coeficientes constantes: si ai j y bi son constantes
• Sistema lineal con coeficientes variables: si ai j y bi dependen de x y t
• Sistema lineal: si ai j depende de x y t, y bi es función lineal de u j
• Sistema cuasi-lineal: si A = A(U)
• Sistema homogéneo: si B = 0
A la hora de resolver el sistema anterior de debe especificar el dominio espacial de
la EDP, [xl , xr ], las condiciones de contorno definidas en xl y xr y la condición inicial.
Esto nos define un problema de valor inicial (PVI) o problema de Cauchy.
Como ejemplos de este tipo de sistemas tenemos la ecuación de convección escalar,
∂u
∂ t
+ a
∂u
∂x
= 0, la ecuación de Burgers no viscosa,
∂u
∂ t
+ u
∂u
∂x
= 0 y las ecuaciones
de Euler unidimensionales para un gas,
∂r
∂ t
+
∂ru
∂x
= 0 ,
∂ru
∂ t
+
∂ (ruu + p)
∂x
= 0 ,
∂rE
∂ t
+
∂ru H
∂x
= 0 ,
donde H = E + p/r la entalpía total del sistema.
Por otro lado, las leyes de conservación son sistemas de ecuaciones en derivadas
parciales que se pueden escribir de la siguiente forma:
U t + F (U)x = 0 , (34)
donde,
U =
2
6
6
6
4
u1
u2
...
um
3
7
7
7
5
, F (U) =
2
6
6
6
4
f1
f2
...
f m
3
7
7
7
5
.
U el vector de variables conservadas, F el vector de flujos y en general f i son funciones
de u j . Es decir, las leyes de conservación pertenecen al conjunto de sistemas homogé-
neos. Estas leyes surgen de la forma integral de las leyes físicas. En cuanto a los flujos,
estos pueden ser lineales o no lineales.
El Jacobiano del sistema se define como,
A (U) =
∂F
∂U
=
2
6
6
6
6
6
6
6
6
4
∂ f1 /∂u1 ∂ f1 /∂u2 . . . ∂ f1 /∂um
∂ f2 /∂u1 ∂ f2 /∂u2 . . . ∂ f2 /∂um
...
...
...
∂ f m /∂u1 ∂ f m /∂u2 . . . ∂ f m /∂um
3
7
7
7
7
7
7
7
7
5
54
DR
AF
T
4.2 sistemas lineales
Dicho Jacobiano permite escribir la ley de conservación (34) como un sistema de
Ecuaciones (33). Para ello, dado,
∂F (U)
∂x
=
∂F
∂U
∂U
∂x
,
entonces,
U t + A (U) Ux = 0 . (35)
Así por ejemplo, la ecuación de convección escalar se puede escribir de forma con-
servativa con tal de tomar f (u) = au, mientras que para la ecuación de Burgers
f (u) = 1/2u2.
4.1.2 Definición de sistema hiperbólico
Definición 2. Sistema hiperbólico: Un sistema de ecuaciones del tipo (35) se denomina hiper-
bólico en un punto (x , t) si la matriz A tiene m valores propios reales siendo el correspondi-
ente conjunto de vectores propios {K i} linealmente independiente. El sistema es estrictamente
hiperbólico si todos los valores propios son distintos. ⇤
Obviamente, un sistema estrictamente hiperbólico es también hiperbólico. Como
contra ejemplo se tienen las ecuaciones de Cauchy-Riemann, ∂u/∂x ± ∂v/∂y = 0.
En este capítulo se estudiara el problema en orden creciente de complejidad: (a)
Ecuación de convección/advección lineal, (b) Sistema hiperbólico lineal de ecuaciones,
(c) Ecuación de convección/advección no lineal.
4.2 sistemas lineales
4.2.1 Ecuación de convección/advección lineal
La ecuación de convección es una EDP escalar, linear y con coeficientes constantes.
Para su estudio se considera el siguiente problema de Cauchy,
ED P ut + a ux = 0, �• < x < • , t > 0 ,
C I u(x , 0) = u0 (x) siendo a la velocidad de propagación
Definición 3. Características: Las características son curvas en el plano x � t a lo largo de
las cuales la EDP se convierte en una EDO. ⇤
Estas curvas tienen la propiedad de propagar la información. El valor de u es cons-
tante a lo largo de las características. Al evaluar u sobre las curvas características
x = x( t), se obtiene una función que sólo depende de t, u = u(x( t) , t), por lo que
se verifica que la variación de u a lo largo de x( t) viene dada por,
du
dt
=
∂u
∂ t
+
dx
dt
∂u
∂x
.
Suponiendo que dx/dt = a, entonces u es constante a lo largo de x( t). Esto conduce
a soluciones de la forma (véase Figura 29),
u(x , t) = u0 (x0 ) = u0 (x � at) . (36)
55
DR
AF
T
integradores de riemann
Figure 29.: Característica de un ecuación lineal con coeficientes constantes.
Esto significa que la solución del problema se corresponde con una traslación de la
condición inicial con velocidad a.
Figure 30.: Condición inicial del problema de Cauchy.
Como caso particular se puede considerar el caso en el que la condición inicial es
discontinua, tal y como se muestra en la Figura 31. Esto define un problema de Rie-
mann,
ED P ut + a ux = 0, �• < x < • ,
t > 0 ,
C I u(x , 0) = u0 (x) =
⇢
u L si x < 0
uR si x > 0
Tal y como sucedía antes, la solución se corresponde con la traslación con velocidad
a de la condición inicial,
u(x , t) = u0 (x0 ) = u0 (x � at) =
⇢
u L si x � at < 0 o a > x/ t
uR si x � at > 0 o a < x/ t
56
DR
AF
T
4.2 sistemas lineales
Figure 31.: Condición inicial para el problema de Riemann para la ecuación de convec-
ción lineal.
En la Figura 32 se muestra la solución del problema. Como se puede ver, la carac-
terística que pasa por el origen define dos regiones en los que u vale u L y uR respecti-
vamente.
Figure 32.: Solución del problema de Riemann de la ecuación de convección lineal.
4.2.2 Sistema lineal de ecuaciones de tipo hiperbólico
Consideramos un sistema de ecuaciones con coeficientes constantes del tipo (33),
U t + A Ux = 0 .
Si el sistema es hiperbólico, entonces la matriz A se puede diagonalizar A = KLK�1,
donde AKi = l i K i ,
57
DR
AF
T
integradores de riemann
L =
0
B
B
B
@
l1 . . . 0
0 . . . 0
...
...
...
0 . . . lm
1
C
C
C
A
Además, el sistema de vectores propios K =
⇥
K1 , K2 , . . . , Km
⇤
define una base de
vectores. Sin pérdida de generalidad, se asumirá que l1 < . . . < l i < . . . < lm1
El sistema de ecuaciones original se puede desacoplar con tal de multiplicarlo por
K�1 y utilizar la diagonalización de A,
Wt + LWx = 0 .
A las componentes del vector W = K�1 U se les denominan variables características.
Cada una de estas variables definen un EDP de la forma,
∂wi
∂t
+ li
∂wi
∂x
= 0, i = 1, . . . , m ,
las cuales tienen asociada una curva característica, x(t), a lo largo de la cual wi es
constante. La forma de definir la curva característica es idéntica a la utilizada en el
caso anterior, es decir, dx/dt = li. La solución formal de dichas ecuaciones viene
dada por,
wi(x, t) = w0i (x0) = w
0
i (x � lit)
donde w0i se obtiene de la definición de W evaluada en t0.
La solución final se obtiene invirtiendo la transformación inicial, U = KW, es decir,
U(x, t) =
m
Â
i=1
wi(x, t)Ki . (37)
Como se puede ver, la condición de hiperbolicidad permite descomponer la solu-
ción en la base de vectores propios. Físicamente el procedimiento anterior permite
descomponer la solución como una combinación de ondas que viajan con la velocidad
determinada por su correspondiente valor propio.
De nuevo se puede plantear el problema Riemann esta vez para cada una de las
ecuaciones características. Sea,
EDP Ut + A Ux = 0, �• < x < •,
t > 0 ,
CI U(x, 0) = U0(x) =
⇢
uL si x < 0
uR si x > 0
donde la condición inicial se representa en la Figura 31. Como se ha visto, la solución
viene dada por,
U(x, t) =
m
Â
i=1
wi(x, t)Ki =
I
Â
i=1
biKi +
m
Â
i=I+1
aiKi , lI < x/t < lI+1 (38)
1 En este capítulo no se considerar el caso degenerado en el que hay al menos dos li iguales. Las dificul-
tades que surgen cuando se estudian sistemas degenerados no se estudian en este capítulo.
58
DR
AF
T
4.2 sistemas lineales
donde los pesos de la expansión o variables características se obtienen de,
wi(x, t) = w0i (x0) = w
0
i (x � lit) =
(
w0,Li = ai si x � lit < 0 o li > x/t
w0,Ri = bi si x � lit > 0 o li < x/t
Típicamente la solución