Agradecimientos

AGRADECIMIENTOS

Sirvan estas breves palabras como muestra de mi sentido y sincero

agradecimiento hacia todas aquellas personas que han colaborado en la realización de
esta tesina.

Gracias a todos y cada uno de los miembros de LaCàN que han contribuido de
una u otra manera a que la tesina salga adelante. En particular, al director Antonio
Huerta, por ofrecerme la posibilidad de trabajar bajo su dirección dentro de este
dinámico grupo de investigación.

Muchísimas gracias en especial a Rubén Sevilla, por la dedicación de su tiempo

en unas múltiples, largas y provechosas sesiones de tutoría, introduciéndome en el
mundo de la investigación científica. Sin su ayuda, esta tesina hoy no existiría.

Gracias también a toda esa otra gente, mi gente, que desde un ámbito más
externo me ha brindado sus ánimos, apoyo y comprensión durante el largo proceso de
confección de este trabajo.

A todos ellos, una vez más, gracias.

 1. Introducción y objetivos

1 INTRODUCCIÓN Y OBJETIVOS

1.1 Motivación

La modelización numérica de problemas de mecánica de fluidos es actualmente

uno de los ámbitos de mayor actividad investigadora en el campo de la ingeniería
computacional. Sus aplicaciones son múltiples y variadas: las industrias aeroespacial y
de la automoción, la meteorología, la hidrología, la oceanografía, la protección del
medio ambiente, el asesoramiento sobre la seguridad de centrales nucleares o incluso el
flujo arterial [1,2]. En cada uno de estos casos, las ecuaciones que gobiernan el
comportamiento del fluido dependen de manera particular del problema investigado,
puesto que es éste el que determina si ciertas propiedades del fluido, como la viscosidad
o la compresibilidad, pueden ser o no despreciadas para el estudio.

Un problema particular de los habitualmente tratados en mecánica de fluidos son

las ecuaciones de Euler. Se trata de un sistema de ecuaciones hiperbólicas de primer
orden no lineales que expresan las leyes de conservación para fluidos no viscosos. Las
ecuaciones de Euler constituyen una descripción muy ajustada de los flujos
compresibles, y, en consecuencia, son habitualmente adoptadas para la modelización de
este tipo de problemas.

La resolución numérica de las ecuaciones de Euler, así como de todos los

problemas hiperbólicos no lineales en general, comporta diferentes dificultades. Por una
parte, el carácter hiperbólico de las ecuaciones dificulta la prescripción de condiciones
de contorno, ya que la información se propaga en estos casos de forma direccional. Por
otro lado, la no linealidad del problema conduce a la posible aparición de
discontinuidades en la solución, incluso si las condiciones iniciales son continuas. Para
obtener unos buenos resultados, el método numérico utilizado ha de ser capaz de
reproducir correctamente estas discontinuidades.

El método de Galerkin discontinuo es en la actualidad uno de los métodos

numéricos más investigados, especialmente en problemas de mecánica de fluidos y de
electromagnetismo. Desde su aparición en 1973 [3], el método de Galerkin discontinuo
ha demostrado ciertas ventajas en la resolución numérica de problemas frente a otros
métodos de elementos finitos convencionales [4]. En particular, la formulación del
método de Galerkin discontinuo, si se combina con algoritmos explícitos de integración
temporal, conduce a la resolución de un sistema resultante cuya matriz de masa global
es diagonal por bloques, cada uno de ellos correspondiente a un elemento de la malla.
Este hecho permite una resolución computacionalmente mucho más eficiente del
problema, facilitando además en gran medida el uso de interpolaciones de alto orden en
la discretización de la solución.

La principal idea en que se basa el método de Galerkin discontinuo es el

planteamiento de la forma débil del problema elemento a elemento. A cambio, la
continuidad de la solución se impone en la forma débil mediante flujos a través de las

1
 1. Introducción y objetivos

caras de los elementos. Desde el punto de vista computacional, y análogamente a la

formulación de otros métodos de elementos finitos convencionales, la interpolación de
la solución no se define en el propio elemento físico, sino que se plantea en un elemento
de referencia. El vínculo entre cada uno de los elementos físicos de la malla y el
elemento de referencia se establece a través de la transformación isoparamétrica. De
esta manera -mediante el elemento de referencia y la transformación isoparamétrica- se
consigue una metodología de resolución de problemas de formulación sencilla y
económicamente muy eficiente.

No obstante, algunos autores indican que el uso de la transformación

isoparamétrica puede ser una fuente de error en la resolución de problemas en dominios
de geometrías curvadas [5]. Según sostienen, este error aumenta con la no linealidad del
problema y con el grado de interpolación adoptado para la solución. En cualquier caso,
conviene remarcar que el error no aparece exclusivamente con el método de Galerkin
discontinuo, sino que también aplica a los métodos de elementos finitos convencionales.

Para evitar la aparición de este tipo de error, en esta tesina se propone el uso de
una formulación alternativa, particularizada para el caso de Galerkin discontinuo, en
que, para cada uno de los elementos curvos, la interpolación de la solución se plantea en
el propio elemento físico. Se evita de esta manera una incorrecta utilización de la
transformación isoparamétrica en la resolución del problema y, en consecuencia, se
garantizan unos menores errores en los resultados obtenidos. Además, el aumento de
coste computacional que implica esta metodología alternativa no es excesivamente
remarcable, puesto que solo es necesario modificar el cálculo de los elementos curvos.

1.2 Objetivos

El objetivo principal de esta tesina es la implementación de un código de

Galerkin discontinuo para la resolución numérica de las ecuaciones de Euler a través de
una formulación en que, para cada uno de los elementos curvos, las funciones de forma
se planteen en el propio elemento físico. Se pretende validar el código implementado
mediante algún test clásico del ámbito de los problemas de flujo compresible y
comparar también sus resultados con los proporcionados por una resolución análoga
realizada con un código de Galerkin discontinuo convencional, con elementos
isoparamétricos.

Con el fin de abordar el problema de una manera gradual y ordenada,

previamente a la resolución de las propias ecuaciones de Euler se resolverán otros
problemas hiperbólicos más sencillos. El grado de complejidad de las ecuaciones
resueltas se incrementará progresivamente, pudiéndose además observar de esta manera
las distintas dificultades que van apareciendo para la resolución numérica de cada tipo
de problema. En concreto, se plantearán, por orden sucesivo, los problemas de
convección lineal pura (problema hiperbólico escalar y lineal), de electrodinámica o
Maxwell (problema hiperbólico vectorial y lineal) y de flujo compresible o Euler
(problema hiperbólico vectorial y no lineal).

2
 1. Introducción y objetivos

1.3 Estructura de la tesina

En los capítulos 2 y 3 se presentan, con una visión bastante general, los distintos
conceptos teóricos en que se basa la presente tesina. En concreto, el capítulo 2 consiste
en una introducción a las ecuaciones de Euler y a otros conceptos teórico-matemáticos
relacionados, como las leyes básicas de conservación y otros problemas hiperbólicos
también tratados en esta tesina. Por otra parte, en el capítulo 3 se describen los
principales fundamentos del método de Galerkin discontinuo.

El capítulo 4 describe y justifica con detalle la nueva formulación de Galerkin

discontinuo implementada en esta tesina y presenta ya algunos sencillos resultados que
permiten validar el código. Pese a emplearse una explicación particularizada para el
método de Galerkin discontinuo, las justificaciones y planteamientos descritos son
también extensibles a otros métodos de elementos finitos convencionales.

En el capítulo 5 se muestran los resultados obtenidos a partir de la aplicación de

la nueva formulación implementada a los distintos problemas hiperbólicos considerados
en esta tesina (convección lineal pura, electrodinámica o Maxwell y flujo compresible o
Euler).

Finalmente, el capítulo 6 recoge las principales conclusiones del estudio

realizado e indica además algunas líneas futuras de investigación sugeridas por los
resultados derivados de esta tesina.

3
 2. Las ecuaciones de Euler

2 LAS ECUACIONES DE EULER

2.1 Leyes básicas de conservación

La base física de la mayoría de sistemas hiperbólicos son las denominadas leyes

de conservación. Se trata de unas ecuaciones dependientes del tiempo que expresan la
conservación de diferentes cantidades (o variables de estado), como por ejemplo la
masa, el momento o la energía, en dinámica de fluidos.

Más concretamente, las variables que aparecen en las ecuaciones de

conservación son las densidades de las variables de estado. Estas funciones representan
la distribución espacial del conjunto de variables para un instante t, y generalmente
cambian con el tiempo. El hecho de que se trate de unas variables conservadas significa
que la integral en todo el dominio de su función de densidad ha de ser constante con el
tiempo. Además de las funciones de densidad, las leyes de conservación también
contienen funciones referentes a los flujos de cada una de las variables de estado, en las
que intervienen además otras variables secundarias.

A continuación se presentan las ecuaciones de conservación correspondientes a

la masa, momento y energía para un fluido compresible no viscoso, expresadas en su
forma diferencial.

• Ley de conservación de la masa o ecuación de continuidad:

∂ρ
+ ∇·( ρ v ) = 0 (1)
∂t

• Ley de conservación del momento o ecuación de movimiento:

∂ ( ρv)
= ρ b − ∇·( pId + ρ v ⊗ v ) (2)
∂t

• Ley de conservación de la energía:

∂(ρE)
= −∇·( ( p + ρ E ) v ) + v·ρ b (3)
∂t

En las ecuaciones anteriores, ρ representa la función de densidad de masa, v es la

función de velocidad, E es la energía total, p es la presión en el fluido y b son las
fuerzas externas por unidad de masa actuantes sobre cada punto.

5
 2. Las ecuaciones de Euler

2.2 Algunos problemas hiperbólicos

Los problemas hiperbólicos formados por leyes conservativas son sistemas de

ecuaciones en derivadas parciales dependientes del tiempo y con una estructura
particularmente sencilla. En general, todo problema de este tipo se puede expresar de
forma compacta mediante la ecuación general

∂U
+ ∇·F = Q , (4)
∂t

donde U es el vector de variables de estado, F es la matriz de flujos, en que cada

columna es el vector de flujos de las variables conservadas en cada dirección cartesiana,
y Q es el vector de término fuente o acciones exteriores.

Típicamente las funciones de flujo son no lineales respecto a U, de manera que

el sistema asociado es también no lineal. En general no es posible obtener soluciones
analíticas para estas ecuaciones, y es por lo tanto necesario recurrir a los métodos
numéricos para obtener soluciones aproximadas.

Existen varias razones para estudiar este tipo particular de ecuaciones [2]. Por
una parte, muchos problemas propios de la ciencia y la ingeniería están relacionados
con la conservación de ciertas cantidades y conducen a ecuaciones de este tipo. Por otra
parte, desde el punto de vista computacional existen ciertas dificultades en la resolución
de este tipo de sistemas (discontinuidades, condiciones de contorno) que no se dan en
ningún otro problema y que son motivo de intensa investigación para el desarrollo de
nuevos métodos numéricos que las superen de una manera eficiente.

Los distintos problemas hiperbólicos analizados en esta tesina difieren en la

expresión que adoptan los vectores U, F y Q de la ecuación 4, en función de las
distintas variables de estado consideradas en el problema. A continuación se presentan
las ecuaciones particulares y se describen las principales peculiaridades de cada caso.

2.2.1 Convección lineal pura

La ecuación que, junto con las correspondientes condiciones iniciales y de

contorno, gobierna el problema transitorio de convección lineal pura es:

ut + ( a ( x, t )·∇ ) u = s ( x, t ) , (5)

donde u es la cantidad conservada (escalar) y a es la velocidad de propagación,

independiente de la solución u y con divergencia nula. De hecho, son estas dos
propiedades de a las que permiten pasar de una fórmula como la de la ecuación de
continuidad (ecuación 1) a esta expresión más simplificada. Las variables U, F y Q que
permiten expresar el problema en los términos de la ecuación 4 quedan definidas como:

U=u; F = au ; Q = s ( x, t )

6
 2. Las ecuaciones de Euler

La principal característica de esta ecuación es que espacio y tiempo están

relacionados por las denominadas líneas características: para una solución u(x,t) y un
punto en el espacio y tiempo (x,t) determinados, se cumple dx dt = a . En consecuencia,
la información u se propaga con una velocidad a como una función del punto espacio-
tiempo considerado. La cantidad u se desplaza pues siguiendo una trayectoria definida
como curva característica. En general, toda cantidad conservada y propagada a través de
una línea característica se denomina variable Riemann o variable característica.

u (x,t)

cte
+
at
x=

u (x,0)

x
Figura 1. Representación esquemática del concepto de línea característica

Al tratarse de una ecuación lineal, las curvas características del problema de

convección pura son curvas fijas en el plano, independientes de la solución. Pese a todo,
la resolución numérica del problema requiere una doble discretización, en espacio y
tiempo. El número de Courant C, que relaciona ambas discretizaciones, y se define
como C = a∆t/h, condiciona la estabilidad de los métodos numéricos explícitos en la
resolución de este problema.

Otra propiedad interesante del problema lineal de convección pura es que,

gracias al concepto de las líneas características, se puede obtener una solución analítica
a través de un planteamiento Lagrangiano del movimiento, en que cada nodo se mueve
asociado a una partícula material, de manera que se elimina todo efecto convectivo.

2.2.2 Electrodinámica (ecuaciones de Maxwell)

Pese a no ser un problema propio de la mecánica de fluidos, las ecuaciones de

Maxwell, formuladas para la descripción de los fenómenos electromagnéticos,
constituyen un problema hiperbólico cuyo estudio resulta de gran interés para
aplicaciones tan diversas como la medicina, la levitación de trenes, o la aviación militar
[6].

7
 2. Las ecuaciones de Euler

Las ecuaciones de Maxwell, formuladas para el vacío, y aplicables también de

manera aproximada para el aire, son:

∂E
ε = ∇×H (ley de Maxwell-Ampère) (6)
∂t
∂H
µ = −∇ × E (ley de Faraday) (7)
∂t
∇ ⋅ (ε E ) = 0 (ley de Gauss) (8)

∇ ⋅ ( µH ) = 0 (ley de Gauss para el campo magnético) (9)

donde E i H son las intensidades del campo eléctrico y magnético

respectivamente, ε es la permitividad eléctrica del medio y µ es su permeabilidad
magnética.

El sistema formado por las ecuaciones 6 a 9 puede ser reformulado de una

manera compacta con la ecuación 4 mediante la siguiente definición de las variables U,
F y Q:

⎛ 0 − H3 H 2 ⎞
⎜ ⎟
⎜ H3 0 − H1 ⎟
⎛ εE ⎞ ⎜ − H 2 H1 0 ⎟
U=⎜ ⎟; F = ( F1 F2 F3 ) = ⎜ ⎟; Q=0
⎝ µH ⎠ ⎜ 0 E3 − E2 ⎟
⎜ − E3 0 E1 ⎟
⎜⎜ ⎟
⎝ E2 − E1 0 ⎟⎠

Para el caso particular de 2 dimensiones, debido a las relaciones de

independencia entre algunas de las ecuaciones, el problema puede ser desacoplado en
dos subproblemas, denominados transverso eléctrico y transverso magnético. Las
variables U y F que definen los respectivos problemas según la expresión de la ecuación
4 son:

• Transverso eléctrico

⎛ ε E1 ⎞ ⎛ 0 − H3 ⎞
⎜ ⎟
U = ⎜ ε E2 ⎟ ; F = ( F1 F2 ) = ⎜ H3 0 ⎟ (10)
⎜ ⎟ ⎜ E −E ⎟
⎝ µ H3 ⎠ ⎝ 2 1 ⎠

• Transverso magnético

⎛ µ H1 ⎞ ⎛ 0 E3 ⎞
⎜ ⎟
U = ⎜ µ H2 ⎟ ; F = ( F1 F2 ) = ⎜ − E3 0 ⎟ (11)
⎜ ⎟ ⎜ −H H ⎟
⎝ ε E3 ⎠ ⎝ 2 1⎠

8
 2. Las ecuaciones de Euler

De la misma manera que en el problema escalar de la convección pura, en el

caso de las ecuaciones de Maxwell la información también se propaga a lo largo de las
líneas características. No obstante, en este caso aparece una dificultad añadida, ya que la
incógnita U es ahora una variable vectorial. Este hecho conduce a unas ecuaciones del
problema acopladas entre sí, lo que provoca que, para un punto en el espacio-tiempo
determinado, existan distintas líneas características (tantas como componentes tiene el
vector de incógnitas) a través de las cuales se propagan distintas variables Riemann
conservadas. Estas variables de Riemann, además, no coinciden exactamente con los
términos que aparecen en el vector de incógnitas del problema, sino que son variables
secundarias calculables a partir de estos. El carácter vectorial de problema repercute
pues en un considerable aumento de la complejidad de su resolución numérica.

t
C-

C0

P (x,t) C+

Figura 2. Líneas características para un problema de tipo vectorial

2.3 Flujo compresible (ecuaciones de Euler)

Las ecuaciones de Euler, que describen los problemas de flujo compresible,

corresponden a las leyes de conservación de masa, momento y energía para fluidos
compresibles no viscosos, descritas por las ecuaciones 1 a 3. El sistema se completa con
una ecuación adicional, la ecuación de estado, que define las propiedades
termodinámicas del fluido, estableciendo una relación entre su presión, densidad,
velocidad y energía. Para el caso de los gases ideales, la ecuación de estado es

⎛ 1 2⎞
p = ( γ − 1) ρ ⎜ E − v ⎟ ,
⎝ 2 ⎠

donde γ es una constante característica del fluido que relaciona sus calores específicos a
temperatura y volumen constante.

El sistema de Euler, formado por las ecuaciones mencionadas, puede reescribirse

en la forma compacta de la ecuación 4 mediante la definición de las siguientes
variables:

9
 2. Las ecuaciones de Euler

⎛ ρ ⎞ ⎛ ρ vi ⎞ ⎛ 0 ⎞
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
U = ⎜ ρv ⎟ ; F = ( F1 F2 F3 ) , con Fi = ⎜ ρ vvi + pei ⎟ ; Q = ⎜ ρb ⎟
⎜ ρE ⎟ ⎜( ρE + p) v ⎟ ⎜ v ⋅ ρb ⎟
⎝ ⎠ ⎝ i⎠ ⎝ ⎠

En los problemas de flujo compresible, se define una nueva variable, el número

de Mach M, que relaciona la velocidad del fluido con la velocidad del sonido c:

v γp
M= , donde c =
c ρ

El número de Mach permite distinguir entre dos tipos de problemas de flujo

compresible: los problemas de flujo supersónico y los de flujo subsónico. En los
primeros, el número de Mach es superior a 1, y por tanto la velocidad del fluido es
superior a la del sonido. En los problemas de flujo subsónico, correspondientes al caso
contrario, el número de Mach es inferior a 1.

La distinción entre estos dos tipos de problemas tiene importantes repercusiones

en su resolución numérica [1,2,7]. Puesto que las ecuaciones de Euler definen un
sistema hiperbólico, la información se propaga, como en los otros problemas descritos
anteriormente, a través de las líneas características. En este caso, pero, al tratarse de un
sistema no lineal, las curvas características no son fijas, sino que varían con el tiempo.
En la formulación completa del problema a resolver, las condiciones de contorno sólo
pueden ser prescritas en aquellos contornos en que la información entra en el dominio.
Para problemas supersónicos, el caso es relativamente sencillo, puesto que se deben
prescribir todos las componentes de U en los contornos de entrada del flujo. En los
problemas de flujo subsónico, en cambio, la situación es más compleja, ya que los
contornos de entrada de información no siempre coinciden con los de entrada de flujo.
La prescripción de condiciones de contorno se convierte pues en una dificultad
numérica añadida en la resolución del problema.

P (x,t) P (x,t)

C+
C+ C0 C-
C0 C-

x
Flujo supersónico Flujo subsónico

Figura 3. Líneas características para problemas de flujo subsónico y supersónico

10
 2. Las ecuaciones de Euler

Otra dificultad que aparece en la resolución de un sistema hiperbólico no lineal,

como son las ecuaciones de Euler, es la posible aparición de discontinuidades en la
solución, incluso partiendo de condiciones iniciales suaves y continuas [1,2,7]. Los
métodos de elementos finitos convencionales no son capaces de reproducir
correctamente estas soluciones discontinuas, por lo que es necesario recurrir a
estrategias alternativas que permitan captarlas con precisión. En el caso concreto de las
ecuaciones de Euler, las discontinuidades sólo aparecen en condiciones de flujo
supersónico, en forma de choques. No obstante, en esta tesina los problemas de flujo
supersónico no serán tratados.

11
 3. Galerkin discontinuo

3 GALERKIN DISCONTINUO

El método de Galerkin discontinuo es un método numérico, situado entre los

métodos de elementos finitos y los métodos de volumen finito, que combina aspectos
positivos de ambos. Se trata de una técnica que, debido a sus características particulares,
resulta especialmente útil para la utilización de altos grados de interpolación.

La principal idea en que se basa el método de Galerkin discontinuo es el

planteamiento de la forma débil del problema elemento a elemento. Equivalentemente,
también puede decirse que se utilizan unas funciones de forma para la interpolación de
la solución discontinuas elemento a elemento. A cambio, la continuidad de la solución
se impone en la forma débil del problema mediante flujos a través de las caras de los
elementos.

El método de Galerkin discontinuo tiene varias ventajas. En primer lugar, el

hecho de que la comunicación entre elementos se produzca únicamente a través de
flujos entre caras permite, si se combina con algoritmos explícitos de integración
temporal, la paralelización del cálculo, ya que la matriz de masa global del sistema
resultante es diagonal por bloques, de manera que la solución en un elemento no se
reconstruye en función de sus elementos vecinos. En todo caso, conviene remarcar que
este hecho no significa que cada elemento pueda resolverse de manera independiente.
La estructura particular de la matriz del sistema no sólo posibilita la paralelización del
cálculo, sino que además también simplifica la adpaptatividad h-p del método,
facilitando la utilización de interpolaciones de alto orden en la discretización de la
solución.

Gracias a estas interesantes propiedades, el método de Galerkin discontinuo

permite una gran flexibilidad en el diseño de mallas y en la elección de los espacios de
aproximación. Incluso puede plantearse el caso de una descomposición mixta del
dominio, en la que para cada subdominio se utilice distinto tipo de interpolación.

El inconveniente más destacable del método de Galerkin discontinuo es la mayor

complejidad y coste computacional que requiere, puesto que deben calcularse integrales
en las caras de cada elemento y, además, existen nodos duplicados en la malla,
correspondientes a los contornos de cada elemento. Este segundo hecho, de todas
maneras, es despreciable en mallas de alto orden, puesto que, en estos casos, el número
de nodos en el contorno de cada elemento es prácticamente irrelevante frente al número
de nodos del interior.

3.1 Formulación débil

A continuación se deriva la formulación débil para a un problema hiperbólico

tipo según el método de Galerkin discontinuo. Por cuestiones de simplicidad, el

13
 3. Galerkin discontinuo

desarrollo se efectúa en base a una ecuación escalar, siendo el mismo planteamiento

extensible a problemas de tipo vectorial.

La ecuación del problema estudiado, en su forma continua y diferencial, es:

∂u ∂fk (u )
+ = q, en Ω × ]0,T [ (12)
∂t ∂xk

donde se adopta la notación de Einstein, es decir, donde un índice repetido

indica una sumatoria sobre todos los valores que este índice puede adoptar.

Premultiplicando la ecuación 12 por una función de test w e integrando en un

elemento Ωe, se obtiene:

∂fk (u )
∫Ωe
wut d Ω + ∫ w
Ωe ∂xk
d Ω = ∫ wqd Ω
Ωe

Integrando por partes y aplicando el teorema de la divergencia de Gauss resulta:

∂w
∫Ωe
wut d Ω − ∫
Ωe ∂x
k
fk (u )d Ω + ∫ wfk (u )nk d Γ = ∫ wqd Ω
∂Ωe Ωe
(13)

El flujo normal entre caras de elementos se aproxima por un flujo numérico fn :

fi (u )ni = fn (u ) ≈ fn (u , u + )

siendo u+ el valor de la solución en el elemento adyacente en la dirección de

propagación de la información (figura 4).

Figura 4. Notación propia del método de Galerkin discontinuo

Existen distintas posibilidades para la definición del flujo numérico. Todas ellas
reproducen, ya sea de manera exacta o aproximada, la conservación de las variables de
Riemann a lo largo de las líneas características, mediante lo que se denominan Riemann

14
 3. Galerkin discontinuo

solvers. Una de las opciones más utilizadas es el Riemann solver propuesto por Lax-
Friedrichs [8]:

1 α
fn (u, u+ ) = ⎡⎣ fn (u ) + fn (u+ ) ⎤⎦ − u ,
2 2

donde α es una viscosidad numérica y u es el salto a través del contorno:

u = u+ - u

Introduciendo el concepto del flujo numérico en la ecuación 13, se obtiene la

expresión de la forma débil del problema para el método de Galerkin discontinuo:

∂w
∫Ωe
wut d Ω − ∫
Ωe ∂xk
fk (u )d Ω + ∫ wfn (u, u+ )d Γ = ∫ wqd Ω
∂Ωe Ωe
(14)

3.2 Discretización espacial

Para hallar la forma semi-discreta del problema, se considera una interpolación

de la solución a partir de sus valores nodales:

nen
u (x, t ) = ∑ N j (x) u j (t ) ,
j

donde Nj es la función de forma asociada al nodo j (vale 1 en este punto y 0 en el

resto de nodos), uj es el valor de la solución en el nodo j, y nen es el número de nodos
del elemento

Sustituyendo la discretización adoptada en la forma débil del problema

(ecuación 14), e imponiendo su cumplimiento para un conjunto de funciones de test
{wj} iguales a las funciones de forma utilizadas en la interpolación de la solución
(wj=Nj):

⎛ ∂u ⎞
(∑ N u ) dΩ +
∂Ni
∫ Ni ⎜ ∑ N j j ⎟ d Ω − ∫ fk j j
Ωe
⎝ j ∂t ⎠ Ωe ∂x
k j
(15)
+ ∫ ∂Ωe
Ni fn
(∑ N u , ∑ N u ) dΓ = ∫
j
j j
j
j
+
j
Ωe
Ni qd Ω i, j = 1...n_nod

se obtiene, tras las pertinentes reagrupaciones, el sistema de ecuaciones diferenciales

ordinarias que permite encontrar la solución del problema en cada elemento:

15
 3. Galerkin discontinuo

∂u div
M + f ( u ) + f flux ( u,u+ ) - s = 0 (16)
∂t

donde ( M )ij = ∫Ω e
Ni N j d Ω
ui = ui

(∑ N u ) dΩ
∂Ni
fidiv = ∫ fk j j (17)
Ωe ∂x
k j

fiflux = ∫
∂Ωe (i n
j
)
N f ∑ N u , ∑ N u dΓ j j
j
j
+
j (18)

si = ∫ Ni qd Ω
Ωe

En el planteamiento descrito, la interpolación de la solución se aplica de manera

estricta en el cálculo de todos los términos de la forma débil, incluidos los
correspondientes al flujo. Esta metodología se conoce con el nombre de full quadrature
[1].

Existe una formulación alternativa para el cálculo de los flujos, tal vez menos
rigurosa desde el punto de vista matemático, pero que ha demostrado ser válida en la
resolución de distintos problemas. En este caso, es el flujo el que, en lugar de calcularse
a partir de la solución interpolada (ecuación 15), se interpola a partir de sus valores
nodales:
nen
fk ( x ) = ∑ N j ( x ) fk , j con fk , j =fk ( u j )
j
nen
fn , j ( x ) = ∑ N j ( x ) fn , j con fn , j =fn ( u j , u+j )
j

En este caso, la expresión de los vectores fdiv y fflux que aparecen en la ecuación
16 se modifica, y en lugar de calcularse como indican las ecuaciones 17 y 18, su cálculo
se realiza según:

( ∑ N f (u )) d Ω = C f (u )
∂Ni
f div ( u ) = ∫ j k, j j
∂xk
xk k
Ωe
j

f flux ( u,u ) = ∫ N ( ∑ N f ( u, u ) ) d Γ = M
+
∂Ωe
i j n, j
+
cara fn ( u,u+ ) ,
j

∂Ni
donde (C ) = ∫
xk
ij Ωe ∂xk
N jdΩ

( Mcara )ij = ∫∂Ω e

Ni N j d Γ

fk,i = fk ( ui )

fn,i = fn ( ui , ui+ )

16
 3. Galerkin discontinuo

Como se observa, el cálculo de los flujos responde a un procedimiento mucho

más sistemático en este caso, lo que conduce a unos algoritmos con un coste
computacional ligeramente inferior al correspondiente a la full quadrature.

Para problemas en los que la incógnita sea una variable vectorial, la forma semi-
discreta se deriva análogamente extendiendo los procedimientos descritos a cada una de
las componentes del vector de incógnitas, obteniéndose una forma matricial resultante
de la misma estructura, pero con matrices en lugar de vectores:

∂U
M + Fdiv (U) + Fflux (U, U+ ) - S = 0 (19)
∂t

3.3 Integración temporal

Para la discretización temporal del problema, es habitual combinar el método de

Galerkin discontinuo con algoritmos explícitos de integración en el tiempo [1]. No
obstante, actualmente muchos autores recomiendan el uso de esquemas implícitos de
integración temporal, que evitan el principal inconveniente de los algoritmos explícitos,
que es la necesidad de usar pasos de tiempo muy reducidos para poder cumplir las
condiciones de estabilidad [9].

En particular, en esta tesina se ha optado por la utilización de esquemas

explícitos tradicionales, los cuales se describen a continuación, particularizados para
una ecuación diferencial ordinaria de expresión genérica:

du
= F (t, u )
dt

• Esquema RK4:

El algoritmo de cálculo del clásico método explícito de Runge-Kutta de cuarto

orden es:

∆t ∆t ⎛ ∆t ⎞ ∆t ⎛ ∆t ⎞
u n +1 = u n + F ( t n , un + β1 ) + F ⎜ t n + , un + β2 ⎟ + F ⎜ t n + , un + β3 ⎟ +
6 3 ⎝ 2 ⎠ 3 ⎝ 2 ⎠
∆t
(
+ F t n + ∆t , un + β4
6
)
donde un + β1 = un

17
 3. Galerkin discontinuo

∆t ⎛ n ∆t n + β1 ⎞
u n + β2 = u n + F ⎜t + ,u ⎟
2 ⎝ 2 ⎠
∆t ⎛ ∆t ⎞
un + β3 = u n + F ⎜ t n + , u n + β2 ⎟
2 ⎝ 2 ⎠
(
un + β4 = un + ∆tF t n + ∆t , un + β3 )
Alternativamente, el algoritmo puede también definirse de manera esquemática
mediante la siguiente tabla de Butcher propia de los esquemas Runge-Kutta:

0
1/2 1/2
1/2 0 1/2
1 0 0 1
1/6 1/3 1/3 1/6

Este método se ha empleado en la resolución de todos los problemas de

electromagnetismo y en los casos de convección en que las interpolaciones espaciales
adoptadas para la solución eran de bajo orden.

• Esquemas explícitos de Padé:

El algoritmo genérico de cálculo de un esquema de Padé de ntg pasos se escribe

como:

un + β1 = un
un + βi = un + βi ∆t un + βi−1 i = 2,..., ntg + 1
n + βntg +1
un +1 = u ,

siendo β1 = 0
βi = 1/(ntg+2-i) i = 2,…,ntg +1

Pese a no tratarse de un método de Runge-Kutta, los esquemas explícitos de

Padé también pueden describirse mediante una tabla de Butcher:

0 0
1/ntg 1/ntg 0
1/(ntg-1) 0 1/(ntg-1) 0
… … 0
1/2 0 ... 0 1/2 0
0 … … 0 1

18
 3. Galerkin discontinuo

En esta tesina, se han utilizado esquemas de Padé de distinto orden para la

resolución de ciertos ejemplos del problema de convección lineal pura, coincidiendo
con el uso de interpolaciones espaciales de alto orden para la solución.

• Esquema TVD-RK3:

Los esquemas TVD (Variación Total Decreciente) son algoritmos de integración

temporal que evitan la aparición de oscilaciones en la solución impidiendo que se den
nuevos extremos en la evolución de las distintas incógnitas.

En esta tesina se ha empleado el esquema TVD-RK3, que aplica el concepto de

los métodos TVD al clásico algoritmo de Runge-Kutta de tercer orden. El algoritmo de
cálculo de este método de integración temporal es:

un + β1 = un
un + β2 = u n + β1 − ∆t F ( t n , u n + β1 )
3 1 ∆t
un + β3 = un + un + β2 − F ( t n , un + β2 )
4 4 4
1 2 2∆t
u n + 1 = u n + u n + β2 − F ( t n , u n + β2 )
3 3 3

El método TVD-RK3 se ha utilizado en la resolución de los problemas de flujo

compresible de esta tesina. Existen estudios detallados sobre las condiciones de
estabilidad de este método [10], cuyas conclusiones han sido aplicadas para determinar
el paso de tiempo óptimo para cada resolución.

19
 4. Implementación de una nueva formulación

4 IMPLEMENTACIÓN DE UNA NUEVA FORMULACIÓN

4.1 Justificación

Para la obtención de las diferentes matrices que aparecen en la forma matricial

del problema estudiado (ecuaciones 16 o 19) deben calcularse, para cada elemento, las
integrales de ciertas funciones, ya sea en todo su dominio o en sus distintas caras. Como
resultado de la discretización del problema, por la que las distintas variables se
interpolan en todo el dominio de un elemento a partir de sus valores nodales, en la
mayoría de ocasiones las funciones a integrar contienen, de alguna u otra manera, las
funciones de forma utilizadas para la interpolación. El cálculo de las integrales se
realiza mediante cuadraturas numéricas.

Con el fin de simplificar el cálculo de estas integrales, las funciones de forma

suelen plantearse habitualmente en un elemento de referencia, de geometría sencilla y
con propiedades muy ventajosas. La cuadratura numérica también se calcula en este
mismo elemento de referencia, y el valor de la integral se particulariza para cada
elemento a través de la denominada transformación isoparamétrica, que establece la
relación entre la geometría de cada elemento en las coordenadas cartesianas (x-y) de la
malla de discretización y la correspondiente a las coordenadas locales (ξ-η) del
elemento de referencia:

nen
⎫
x (ξ ,η ) = ∑ N j (ξ ,η ) xnod,j ⎪
j =1 ⎪
nen ⎬ transformación isoparamétrica
y (ξ ,η ) = ∑ N j (ξ ,η ) ynod,j ⎪
j =1
⎪⎭

siendo nen el número de nodos del elemento, que depende del grado de
interpolación p adoptado para la solución (en el caso de elementos bidimensionales
triangulares, nen=(p+2)(p+1)/2).

Figura 5. Representación conceptual de la transformación isoparamétrica

21
 4. Implementación de una nueva formulación

Así pues, el cálculo de las distintas integrales se plantea a través de una

transformación de la integral sobre el dominio x-y del elemento físico al dominio ξ-η del
elemento de referencia, mediante la aplicación de la transformación isoparamétrica:

∫ f ( x (ξ ,η ) , y (ξ ,η ) ) dxdy = ∫ f (ξ ,η ) J (ξ ,η ) dξ dη (20)
Ωe ΩRef

donde f es la función a integrar y J(ξ,η) es la matriz jacobiana de la

transformación isoparamétrica:

⎛ ∂x ∂x ⎞ ⎛ ∂N j (ξ ,η ) x ∂N j (ξ ,η ) ⎞
nen nen

⎜ ∂ξ ∂η ⎟ ⎜ ∑ ∂ξ
nod,j ∑ ∂η
xnod,j ⎟
J (ξ ,η ) = ⎜ ⎟ = ⎜ j =1 j =1 ⎟ (21)
⎜ ∂y ∂y ⎟ ⎜ nen ∂N j (ξ ,η ) j (
nen ∂N ξ ,η
) ⎟
⎜ ∂ξ ∂η ⎟ ⎜⎜ ∑ ynod,j ∑ ynod,j ⎟⎟
⎝ ⎠ ⎝ j =1 ∂ξ j =1 ∂η ⎠

Para poder computar la integral según la expresión de la ecuación 20, es

necesario seguir una serie de pasos:

i) obtención de los puntos y pesos de Gauss de una cuadratura bidimensional de

orden p en el elemento de referencia:

( (ξg,1 ) ( (ξ
,ηg,1 ) ,w g,1 , g,2 )
,ηg,2 ) ,w g,2 …… ( (ξ
g,n_pg ,ηg,n_pg ) ,w g,n_pg )
donde n_pg es el número de puntos que conforman la cuadratura.

ii) obtención de las funciones de forma en el elemento de referencia y evaluación

de éstas y sus derivadas en los puntos de Gauss correspondientes:

Ni (ξg,1 ,ηg,1 ) , Ni (ξg,2 ,ηg,2 ) … Ni (ξg,n_pg ,ηg,n_pg ) i = 1… nen

∂Ni ∂N ∂Ni
∂ξ
( ξg,1 ,ηg,1 ) , i (ξg,2 ,ηg,2 ) …
∂ξ ∂ξ
(ξg,n_pg ,ηg,n_pg ) i = 1… nen
∂Ni ∂N ∂Ni
∂η
( ξg,1 ,ηg,1 ) , i (ξg,2 ,ηg,2 ) …
∂η ∂η
(ξg,n_pg ,ηg,n_pg ) i = 1… nen
A partir de los datos obtenidos en i) y ii), la integral en el elemento de referencia
a través de la cuadratura numérica:
n_pg

∫ f (ξ ,η ) J (ξ ,η ) d ξ dη = ∑ f (ξg ,i ,ηg ,i ) J (ξg ,i ,ηg ,i ) wi (22)

ΩRef
i =1

donde el valor de f en los puntos de Gauss se obtiene como interpolación de los

valores nodales a través de las funciones de forma previamente calculadas.

22
 4. Implementación de una nueva formulación

- en el caso en que el elemento sea de lados rectos, el jacobiano de la

transformación es constante en todo el elemento y, en consecuencia, puede salir
fuera de la integral, de manera que se cumple:
n_pg

∫ Ωe
f ( x, y ) dxdy = J ∑ f (ξ
i =1
g ,i ,ηg ,i ) wi = J ∫
ΩRef
f (ξ ,η ) dξ dη

- si por el contrario el elemento es de lados curvos, tal simplificación no es

aplicable y la integral debe calcularse según indica la ecuación 22.

Desde un punto de vista computacional, esta metodología de cálculo resulta

realmente eficiente, puesto que únicamente es necesario calcular las funciones de forma
y sus derivadas, así como los puntos y pesos de la cuadratura numérica, en un único
elemento. Además, la utilización del elemento de referencia simplifica tanto la
obtención y evaluación de las propias funciones de forma y derivadas, que pueden
calcularse fácilmente a partir de los polinomios ortogonales de Legendre, como el
cálculo de integrales, pues se pueden aplicar de manera casi directa los puntos y pesos
de integración de las cuadraturas de Gauss.

No obstante, esta formulación (de ahora en adelante, formulación

isoparamétrica) entraña alguna carencia conceptual. El problema aparece cuando se
pretenden integrar derivadas cartesianas de las funciones a integrar. A continuación se
deducen las operaciones a realizar para el cálculo de la integral siguiendo el mismo
planteamiento:

∂f ( x (ξ ,η ) , y (ξ ,η ) )
∫Ωe ∂xk
dxdy =

⎛ ∂f (ξ ,η ) ∂ξ ∂f (ξ ,η ) ∂η ⎞
=∫ ⎜ + ⎟ J (ξ ,η ) d ξ dη =
⎝ ∂ξ ∂xk ∂η ∂xk ⎠
ΩRef

= ∑⎜
( g ,i g ,i ) ∂ξ + ∂f (ξg ,i ,ηg ,i ) ∂η
n_pg ⎛ ∂f ξ ,η ⎞
⎟ J ( ξ g , i ,η g , i ) w i
i =1 ⎜ ∂ξ ∂xk ∂η ∂xk ⎟
⎝ ⎠

donde ∂ξ ∂xk y ∂η ∂xk se obtienen invirtiendo la matriz jacobiana de la

transformación isoparamétrica (ecuación 21).

Es importante destacar que, según el procedimiento anterior, en el cómputo de

los términos ∂ξ ∂xk y ∂η ∂xk se realiza la inversión de una matriz cuyas distintas
componentes son sumatorias de polinomios. Los términos ∂ξ ∂xk y ∂η ∂xk son pues
dos funciones racionales. Por otra parte, el cálculo de la integral se realiza a través de
una cuadratura de Gauss, que sólo garantiza una integración exacta de términos
polinómicos. Existe pues una fuente de error en el cálculo anterior.

En todo caso, conviene remarcar que el error sólo aparece cuando el elemento
considerado tiene algún lado curvo, puesto que, de lo contrario, las funciones a invertir

23
 4. Implementación de una nueva formulación

son derivadas de polinomios de primer orden y, en consecuencia, los términos

considerados sólo intervienen como una constante.

Algunos autores indican que el error cometido por la formulación isoparamétrica

aumenta con la no linealidad del problema y con el grado de interpolación adoptado de
la solución [5]. Según afirman, una formulación alternativa, que evite el uso incorrecto
de la transformación isoparamétrica en el cálculo de las integrales, proporcionará unos
resultados ostensiblemente mejores.

4.2 Descripción

La formulación implementada en la presente tesina pretende evitar la aparición

del error descrito en el apartado anterior, mediante un planteamiento en que, para los
elementos de lados curvos, tanto las funciones de forma como sus derivadas se calculan
directamente en el elemento físico, y en que la integración de las distintas funciones se
realiza siempre en el propio dominio x-y. A modo de abreviación, a esta nueva
formulación se la denominará formulación x-y. Conviene remarcar que el hecho de que
las modificaciones se introduzcan únicamente en el cálculo de los elementos curvos
permite que el aumento de coste computacional que la nueva formulación requiere no
sea muy significativo.

A continuación se describen las dos principales e importantes diferencias que

introduce la formulación x-y en el cálculo de los elementos curvos respecto a la
metodología propia de la formulación isoparamétrica.

i) Funciones de forma en el dominio x-y

Debido a la necesidad de evitar un cálculo indirecto (mediante un cambio de

variable como es la transformación isoparamétrica) de las derivadas cartesianas de las
funciones a integrar, las funciones de forma adoptadas para la interpolación de la
solución, así como también sus derivadas, se plantean, para cada uno de los elementos
curvos de la malla, en el propio elemento físico (dominio x-y).

El cálculo de las funciones de forma se realiza en base a su propia definición: se

considera un conjunto de polinomios genéricos de grado p (tantos como nodos tenga el
elemento) y se impone que cada uno de ellos adopte un valor unidad en las coordenadas
correspondientes a su nodo asociado y un valor nulo en las coordenadas del resto de
nodos. Los coeficientes de cada polinomio se obtienen tras la resolución del sistema
lineal de ecuaciones resultante. Una vez obtenidas las funciones de forma, el cálculo de
sus derivadas se obtiene de manera directa, mediante simple aplicación de las normas
del cálculo diferencial.

La figura 6 compara, a través de su representación gráfica, las funciones de

forma asociadas al mismo nodo de un mismo elemento obtenidas, respectivamente, a
través del elemento de referencia y tras la posterior aplicación de la transformación

24
 4. Implementación de una nueva formulación

isoparamétrica o bien directamente en el dominio x-y. Se observa que, pese a que ambas
cumplen las propiedades exigidas, existen notables diferencias entre ellas.

a)

b) c)

Figura 6. Representación de una misma función de forma de un elemento curvo (a)

obtenida con una formulación isoparamétrica (b) o con una formulación x-y (c)

ii) Integrales en el dominio x-y

Puesto que las funciones de forma se definen en cada uno de los elementos
físicos curvos de la malla, las distintas integrales que intervienen en la resolución del
problema han de ser calculadas también en este dominio. Este hecho implica que la
cuadratura numérica debe plantearse en las coordenadas globales x-y, lo que tiene unas
importantes repercusiones numéricas. A continuación se describen con detalle estas
implicaciones, particularizadas por cuestión de simplicidad para dominios en dos
dimensiones.

En primer lugar, los puntos y pesos de la cuadratura deben calcularse de manera

particular para cada elemento físico en las propias coordenadas globales x-y. Su
obtención se realiza aplicando una transformación isoparamétrica a los puntos y pesos
que conforman la cuadratura del orden requerido en el elemento de referencia:

nen
xg ,i = ∑ NRef, j (ξg ,i ,ηg ,i ) xnod,j i = 1,..., n_pg
j =1

25
 4. Implementación de una nueva formulación

nen
yg ,i = ∑ NRef, j (ξg ,i ,ηg ,i ) ynod,j i = 1,..., n_pg
j =1

w xy ,i = w Ref,i J (ξg ,i ,ηg ,i ) i = 1,..., n_pg

Puesto que esta transformación es directa y no aparecen matrices inversas en

ningún momento, su uso no introduce error alguno en la nueva formulación.

En segundo lugar, y como consecuencia de esta transformación de puntos y

pesos previa, debe considerarse que el orden de la cuadratura a utilizar para calcular de
manera exacta la integral deseada no es el mismo que el de una formulación
isoparamétrica convencional.

Las funciones de mayor orden que han de ser integradas para calcular los
distintos términos que aparecen en la forma semi-discreta de un problema con Galerkin
discontinuo corresponden a las matrices de masa, y son polinomios de grado igual o
inferior a 2p.

Cuando la integral se calcula en el elemento de referencia (formulación

isoparamétrica), las funciones a integrar se evalúan sobre puntos que corresponden a su
propio dominio de definición (coordenadas ξ-η). En este caso, que no reviste mayor
dificultad, la integración exacta de un polinomio de grado 2p se conseguirá con una
cuadratura de Gauss de orden 2p en cada dirección, esto es, con p+1 puntos de Gauss en
cada dirección.

Si, por el contrario, la cuadratura se plantea en las coordenadas globales x-y del
elemento físico, la función a integrar se evalúa en unos puntos que se han obtenido a
través de una transformación isoparamétrica. Los pesos que definirán la cuadratura,
además, también resultan de una transformación. Así, la función que de hecho pretende
ser integrada es una composición de funciones:

f ( x, y ) = f ( x (ξ ,η ) , y (ξ ,η ) ) J (ξ ,η )

Esta expresión resulta ser un polinomio de grado n·p+2(p-1), donde n es el grado

de la función f definida en las coordenadas x-y. Particularizando el valor de n para el
caso de las matrices de masa (n=2p), la formulación x-y deberá emplear una cuadratura
de Gauss que integre exactamente un polinomio de grado 2p2+2(p-1). Considerando las
propiedades de las cuadraturas gaussianas, esto significa que deberán emplearse p2+ p
puntos de integración en cada dirección. Como puede apreciarse, se trata de una
cantidad claramente superior a los p+1 puntos que requiere la formulación
isoparamétrica, especialmente para interpolaciones de alto orden.

Las modificaciones que introduce la formulación x-y repercuten en un código de

calculo con un coste computacional ligeramente superior al de una formulación
isoparamétrica. No obstante, si bien es cierto que las funciones deben evaluarse un
mayor número de puntos de Gauss en cada uno de los elementos curvos, conviene

26
 4. Implementación de una nueva formulación

recordar que los cambios en el código solo son necesarios para este tipo de elementos,
ya que el error a evitar no puede aparecer en elementos de lados rectos. Dado que el
número de elementos curvos en las mallas de discretización habituales suele ser muy
inferior al de elementos de lados rectos, se puede considerar que el coste computacional
no es un factor que marque grandes diferencias entre las dos metodologías.

La formulación x-y ha sido implementada en distintos códigos de Matlab para la

resolución de los tres problemas tratados en esta tesina. En todos los casos, los códigos
están diseñados con un carácter generalista, que puede ser aplicado de manera indistinta
a la resolución de diferentes ejemplos prácticos, y que permite además la elección de los
valores que adoptan ciertas variables características de cada problema y también de
diversas opciones de cálculo, como el grado de interpolación de la solución, el flujo
numérico considerado o el esquema de integración temporal a utilizar.

4.3 Validación

De cara a comprobar la validez y precisión de sus resultados, la formulación x-y

implementada en la presente tesina ha sido sometida a algunos tests de validación.
Además, con el objetivo de evaluar el grado de mejora que representa respecto a la
habitual formulación isoparamétrica, en algunos de los ejemplos resueltos se han
comparado las soluciones proporcionadas por ambas metodologías.

Test de cuadratura

Un primer proceso de validación, previo a la resolución de cualquier problema,

ha consistido en comprobar que el cálculo de integrales se realiza correctamente a través
de las cuadraturas planteadas. Para ello se han utilizado unos elementos un tanto
especiales, como los que aparecen en la figura 7. Se trata de elementos bidimensionales
de tres lados, dos de los cuales rectos y ortogonales entre sí, y el tercero siguiendo la
curva de un polinomio de expresión conocida. De esta manera es posible calcular por
métodos analíticos la integral exacta de cualquier polinomio en todo su dominio.

Figura 7. Algunos elementos usados en el test de cuadratura

27
 4. Implementación de una nueva formulación

Por otra parte, la integral también se ha calculado numéricamente a través del

planteamiento propio de la formulación x-y. De hecho, este cálculo se ha realizado
utilizando diferentes grados de interpolación en el elemento, siempre inferiores o
iguales al grado del polinomio que de fine su tercera cara. Posteriormente, además, se
han aplicado sucesivos refinamientos de malla en el dominio del elemento, repitiéndose
en cada caso el cálculo de la integral (figura 8).

Malla 0 Malla 1 Malla 2

Malla 3 Malla 4

Figura 8. Sucesivos refinamientos de malla en un elemento del test de cuadratura interpolado con p =1

En analogía al grado del polinomio que debería integrarse correctamente en el

cálculo de una matriz de masa, las funciones integradas en este test de comprobación
han comprendido toda una base de polinomios de grado igual al doble del grado de
interpolación utilizado (2p).

La comparación de los resultados proporcionados por las distintas integraciones

numéricas y analíticas conduce a las siguientes conclusiones:

- Una función únicamente puede integrarse correctamente en el dominio

considerado si el grado de interpolación utilizado es igual o superior al grado
del polinomio que define su cara curvada. En estas circunstancias, los errores
cometidos por la cuadratura implementada son nulos, y por lo tanto se puede
afirmar que la integración es exacta.

- Cuando la interpolación utilizada es de grado inferior al del polinomio que

define la cara curvada del elemento, el error de la integración disminuye al
aumentar el grado de interpolación. En este caso, además, a medida que la

28
 4. Implementación de una nueva formulación

malla del elemento se refina, el error disminuye de una manera aproximada

de acuerdo con las fórmulas correspondientes a las cuadraturas compuestas
de Gauss. Como muestra de algunos de estos resultados, en la figura 9 se
representan dos de los gráficos obtenidos, en los que se comparan las
variaciones teóricas y realmente computadas del error a medida que se refina
la malla. Los dos gráficos corresponden a distintas funciones integradas y a
un distinto grado de interpolación.

Variación del error en el cálculo de la Variación del error en el cálculo de la

integral de una función polinóm ica de integral de una función polinóm ica de
grado 3 con interpolación p =2 grado 7 con interpolación p =4

1.E+00 1.E+00
1.E-01
1.E-01
1.E-02
error / error malla 0

error / error malla 0

1.E-02 1.E-03
1.E-04
1.E-03
1.E-05
1.E-04 1.E-06
1.E-07
1.E-05
1.E-08
1.E-06 1.E-09
0 1 2 3 4 5 0 1 2 3 4 5
nº de malla nº de malla

Variación computada del error Variación computada del error

Variación teórica del error Variación teórica del error
.
Figura 9. Gráficos de variación teórica y computada del error
de integración en función del refinamiento de malla

A partir de estos resultados, se considera probada la validez de la cuadratura

implementada en la formulación x-y.

Test de convergencia

Un segundo test que ha servido para validar la formulación implementada en

esta tesina ha sido la resolución de un sencillo problema de convección lineal pura,
gobernado por la ecuación 5, en el que una onda periódica se propaga a través de un
canal semicircular.

En este ejemplo se ha adoptado una velocidad de transmisión de la onda paralela

a los lados curvos del canal, esto es, tangente a la semicircunferencia definida en cada
punto, y de magnitud variable, igual al radio. En uno de los lados rectos del canal se
impone la condición de contorno de entrada de la onda (condición Dirichlet), mientras

29
 4. Implementación de una nueva formulación

que en el lado opuesto, por donde sale la onda, se considera una condición absorbente,
que no impone restricción alguna. Por otra parte, el tratamiento del término del flujo de
las ecuaciones se ha realizado mediante el planteamiento de interpolación de sus valores
nodales.

Figura 10. Representación esquemática del problema

de convección lineal pura analizado

Como se ha comentado en la descripción general de los problemas de

convección lineal pura (aparatado 2.2.1), al tratarse de un problema lineal y escalar, es
posible encontrar una solución analítica para este ejemplo. Así, para cada una de las
soluciones numéricas obtenidas, se ha calculado su error, computándose en términos de
la norma L2 y de la norma del máximo. En la presentación de los resultados, sin
embargo, y siguiendo el mismo criterio que se adoptará a lo largo de toda la tesina, se
mostrarán sólo los valores correspondientes a la norma L2.

a)

b)

Figura 11. Representación gráfica de dos soluciones numéricas con la formulación x-y
para el problema de propagación de una onda periódica en un canal semicircular:
a) malla fina, p=3; b) malla grosera , p=9

30
 4. Implementación de una nueva formulación

Puesto que el problema analizado es de carácter transitorio, la solución es

dependiente del tiempo. La figura 11 muestra, para un instante t determinado, el aspecto
que ofrece la solución proporcionada por dos resoluciones numéricas con la
formulación x-y, con distintas mallas y grado de interpolación. Como se observa, la
propagación de las ondas se captura en ambos casos perfectamente, sin que aparezcan
problemas de difusión ni oscilaciones en la solución.

El problema se ha resuelto numéricamente de manera alternativa mediante la

formulación isoparamétrica y la formulación x-y Conviene remarcar que, en cada una de
las resoluciones, se ha tenido que garantizar que el paso de tiempo adoptado en la
discretización temporal del problema era lo suficientemente pequeño para que su error
pudiera ser despreciado frente al correspondiente a la discretización espacial.

Con el objetivo de analizar la convergencia del nuevo código, se han adoptado

dos planteamientos distintos. Por una parte, para un mismo grado de interpolación, el
ejemplo se ha resuelto utilizando sucesivos refinamientos de malla. Se obtiene de esta
manera información sobre la convergencia en h del método. La figura 12 muestra los
gráficos correspondientes para unos grados de interpolación iguales a 1, 2 y 3

Convergencia en h para p =1 Convergencia en h para p =2

0.0 -1.0

-0.5
-1.5
)

)
L2

L2

-1.0
log (error

log (error

-2.0
-1.5

-2.5
-2.0
1.94 2.90
-2.5 -3.0
-1.5 -1 -0.5 0 0.5 -1.5 -1 log (h) -0.5 0
log (h)
Isoparamétrica xy Isoparamétrica xy

Convergencia en h para p =3
-1.5

-2.0
) L2

-2.5
log (error

-3.0

-3.5 3.14

-4.0
-1.3 -1.1 -0.9 -0.7 -0.5
log (h)
Isoparamétrica xy

Figura 12. Gráficos de convergencia en h para el problema de propagación

de una onda periódica en un canal semicircular con p = 1, 2 y 3

31
 4. Implementación de una nueva formulación

A partir de los resultados anteriores, se puede considerar que la formulación

implementada obtiene unas convergencias muy aproximadas a las que teóricamente
deberían cumplirse. En teoría, si se desprecia toda influencia del error de integración
temporal, para un grado de interpolación p dado, el logaritmo del error cometido por el
método debería disminuir en función del logaritmo del tamaño h de la malla a razón de
una pendiente igual a p+1 [11]. En los casos de interpolación con grado 1 y 2, los
resultados obtenidos encajan perfectamente con la teoría. En el caso de p=3, en cambio,
la convergencia obtenida dista ligeramente del valor teórico a alcanzar. No obstante, la
diferencia podría atribuirse a que fueran necesarias discretizaciones espaciales y/o
temporales aún más finas que las utilizadas para llegar a la convergencia asintótica.
Debido al elevado tiempo que estos cómputos requieren, estos datos no han podido ser
incluidos en esta tesina.

Otro resultado que puede observarse en las gráficas de convergencia en h es que,

en los casos analizados, las diferencias entre la formulación isoparamétrica y la
formulación x-y son muy pequeñas. El hecho se debe a que se ha trabajado con unos
grados de interpolación especialmente bajos, donde la diferencia entre formulaciones es
muy poco significativa.

Por otro lado, y dado que para la resolución de problemas con el método de
Galerkin discontinuo resulta especialmente interesante el trabajo con interpolaciones de
alto orden, se ha contrastado también la convergencia en p de ambas formulaciones.
Para ello, el problema se ha resuelto con diferentes grados de interpolación, utilizando
siempre una misma malla de discretización del dominio. Puesto que, como debería
esperarse de la teoría, las diferencias entre las dos formulaciones son más evidentes
cuanto más grandes son los elementos, se ha utilizado en este caso la malla más grosera
de las utilizadas en el análisis anterior, formada toda ella por elementos curvos (figura
11b).

Un aspecto referente a este análisis que merece ser remarcado es que, puesto que
se pretendía trabajar con altos grados de interpolación, se ha empleado una distribución
particular de nodos en los elementos. En lugar de trabajar con los habituales nodos
equiespaciados, se han adoptado en este caso los denominados puntos de Fekete. Según
indican algunos autores, la utilización de estos puntos proporciona, para interpolaciones
de alto orden, mejoras considerables en los resultados [12]. Como prueba de ello se
presenta la figura 13, en que se comparan los errores cometidos por la formulación
isoparamétrica convencional en la resolución del problema de convección lineal pura
descrito, utilizando las dos distribuciones de nodos comentadas.

Una simple observación del gráfico permite constatar, cómo, efectivamente, para
interpolaciones de alto orden, los nodos de Fekete muestran mejores resultados que la
distribución equiespaciada habitual. No obstante, conviene remarcar que, aún
proporcionando mejores resultados, los nodos de Fekete, a partir de cierto grado de
interpolación, no siguen la tendencia de disminución del error evidenciada hasta el
momento sino que el este parece estancarse. Probablemente el fenómeno pueda
atribuirse al error introducido por el uso de la transformación isoparamétrica.

32
 4. Implementación de una nueva formulación

Convergencia en p
1.0

0.5

0.0
Isoparamétrica
)

-0.5 equiespaciada
L2
log (error

-1.0

-1.5 Isoparamétrica
Fekete
-2.0

-2.5

-3.0
0 2 4 6 8 10 12 14 16
p

Figura 13. Contraste de resultados entre el uso de nodos equiespaciados o de

Fekete en a resolución del problema de propagación de una onda periódica
en un canal semicircular mediante la formulación isoparamétrica

Los resultados de la figura 13 avalan plenamente la conveniencia de trabajar con

una distribución de nodos de Fekete para altos grados de interpolación. En
consecuencia, para contrastar los datos referentes a la convergencia en p de las dos
formulaciones comparadas en esta tesina se ha considerado apropiada y necesaria la
utilización de estos puntos.

La figura 14 muestra y compara los gráficos de convergencia en p obtenidos con

las formulaciones isoparamétrica y x-y. El gráfico se presenta en dos formas distintas:
además de representar la evolución del error con del grado de interpolación p utilizado,
se ha considerado interesante graficar también esta variación en función del número de
grados de libertad (ndof) utilizados en cada caso, puesto que es ésta una variable mucho
más indicativa de la diferencia de coste computacional entre las distintas resoluciones.

Con el fin de limitar el número de ejecuciones del código, debido a la limitación

de tiempo disponible para la realización de este estudio, se ha adoptado a lo largo de
toda la tesina el criterio de realizar el número de resoluciones mínimas necesarias hasta
alcanzar por primera vez un error igual o inferior a una cierta tolerancia dada. En este
caso, las resoluciones con la formulación x-y se han detenido al encontrar por primera
vez un error menor al mínimo que se había logrado conseguir, para el mismo ejemplo,
con la formulación isoparamétrica (en la figura 13, p=11).

Los resultados de la figura 14 permiten observar como, confirmándose la

hipótesis principal de esta tesina, la formulación x-y proporciona unos mejores
resultados que la formulación isoparamétrica. La tolerancia mínima impuesta para la
solución se ha alcanzado con la formulación x-y con un grado de interpolación inferior
en 2 unidades al requerido por la formulación isoparamétrica o, equivalentemente, con
138 grados de libertad de menos, lo que representa aproximadamente un 30% menos
respecto al total de nodos empleados por la formulación isoparamétrica.

33
 4. Implementación de una nueva formulación

Convergencia en p Convergencia en p
0.0 0.0

-0.5 -0.5

-1.0 -1.0
)

)
L2

L2
-1.5 -1.5
log (error

log (error
-2.0 -2.0

-2.5 -2.5

-3.0 -3.0

-3.5 -3.5
0 2 4 6 8 10 12 0 100 200 300 400 500
p ndof

Isoparamétrica Fekete Isoparamétrica Fekete

xy Fekete xy Fekete

Figura 14. Gráficos de convergencia en p para el problema de

propagación de una onda periódica en un canal semicircular

Se observa también en la figura que las diferencias entre formulaciones son más
evidentes conforme el grado de interpolación de la solución aumenta, llegándose a
alcanzar, para interpolaciones de alto orden, unos errores hasta tres veces menores con
la formulación x-y que con la isoparamétrica.

Tras la validación del código mediante los distintos tests implementados, y tras
observar también unos primeros y satisfactorios resultados en la comparación de las dos
formulaciones contrastadas, se considera suficientemente probada la justificación y
validez de la nueva formulación implementada y se plantea a continuación su uso en
otros problemas o aplicaciones de mayor complejidad.

34
 5. Aplicaciones y resultados

5 APLICACIONES Y RESULTADOS

5.1 Maxwell

Las ecuaciones de Maxwell (ecuaciones 6 a 9), formuladas para la descripción

de los fenómenos electromagnéticos, constituyen un problema hiperbólico lineal con
incógnita vectorial. Información más detallada sobre las distintas formulaciones de estas
ecuaciones y sus propiedades particulares puede consultarse en el apartado 2.2.2.
Conviene en todo caso remarcar la principal diferencia entre este problema y el de
convección lineal pura previamente estudiado, que es el paso de una incógnita de tipo
escalar en la convección a una de tipo vectorial en electromagnetismo. El hecho tiene
importantes repercusiones, puesto que el hecho de trabajar con ecuaciones acopladas
entre sí aumenta en gran medida la complejidad de la resolución numérica.

Una de las aplicaciones de las ecuaciones de Maxwell, destinada a un uso

fundamentalmente militar, es el cálculo de la denominada Radar Cross Section (RCS)
de los aviones [6]. La RCS es un parámetro que describe la extensión a la que un objeto
refleja una onda electromagnética incidente. Se trata de una medida de la intensidad de
la señal detectada por un radar tras la emisión sobre un objeto de una onda de potencia
determinada. En sentido estricto, la RCS se define como el área de la sección transversal
de una esfera perfectamente conductora que reproduce la misma magnitud de reflexión
que la observada con objeto real. La RCS es una propiedad física inherente a cada
objeto y depende de la forma y materiales del objeto y de la longitud de onda y el
ángulo de incidencia de la onda interceptada.

La potencia transmitida por un radar al emitir una onda disminuye a medida que
ésta se propaga, ya que la onda se va expandiendo en una área transversal cada vez
mayor, de manera proporcional al cuadrado de la distancia que atraviesa. Cuando la
onda alcanza el objeto, sólo una parte de ésta lo intercepta. Además, en general, no toda
la onda interceptada se releja, sino que, en función de los materiales, parte de la onda
puede ser absorbida por el obstáculo. Finalmente, una vez reflejada, la onda se expande
de nuevo en el trayecto de vuelta hasta el radar, de manera que éste recibe sólo una
porción de la energía que inicialmente había emitido.

Figura 15. Ondas emitidas por el radar y reflejadas por el objeto

35
 5. Aplicaciones y resultados

Técnicamente, la RCS es un área, y se expresa en metros cuadrados (m2). Sin

embargo, existe una manera habitual y alternativa de expresar este parámetro, que es en
decibelios por metro cuadrado (dBm2). La relación entre ambas unidades responde a:

(
RSC ( dBm2 ) = 10 log10 RSC ( m2 ) )
Una manera habitual de representar la RCS de un objeto es mediante un gráfico
que muestra el valor de la RCS, expresado en dBm2, para cada punto del contorno del
objeto, que se parametriza a través de un ángulo Φ.

Figura 16. Representación gráfica habitual de la RCS

El cálculo analítico de la RCS únicamente es posible para objetos de geometría

muy sencilla. Para casos más complejos, como vehículos, naves o aviones, es necesario
recurrir a simulaciones numéricas. En estos casos, el valor de la RCS (en m2) se calcula
para cada punto del contorno del objeto a partir de la siguiente expresión,
particularizada para el problema del transverso eléctrico en dos dimensiones:

2
H3I
RCS ( m2 ) = lim 2π r 2
r →∞
H3S

donde H3I y H3S son, respectivamente, las intensidades de la tercera componente

del campo magnético incidente y reflejado en ese punto. Para el problema del transverso
magnético, el valor de la RCS se calcula análogamente sustituyendo la componente
analizada del campo magnético (H3) por la que se corresponde en el campo eléctrico
(E3).

En esta tesina se presenta la resolución de dos ejemplos distintos del problema

de las ecuaciones de Maxwell en dos dimensiones. En el primero, en el que el objeto

36
 5. Aplicaciones y resultados

interceptado por las ondas es un sencillo círculo, y que por tanto dispone de solución
analítica, se evalúa el error cometido por el código y se contrasta con el que se obtiene
con una formulación isoparamétrica. El segundo ejemplo corresponde a una forma más
compleja como es una ala de avión, y, lejos de todo cálculo de errores, simplemente se
muestra como prueba de la aplicabilidad del código a geometrías más complejas.

En el caso correspondiente al círculo, se ha simulado la emisión de una onda

longitud de onda igual a 1 propagándose en dirección horizontal de izquierda a derecha.
La llegada de la onda incidente y su posterior reflexión se simulan con la imposición de
una condición de conductor eléctrico perfecto (PEC) sobre el contorno del obstáculo,
que prescribe unos campos eléctrico y magnético totales paralelo y perpendicular,
respectivamente, a la normal exterior al objeto (n):

n × ET = 0 , donde ET = EI + ES
n ⋅ HT = 0 , donde HT = HI + HS

Las ecuaciones de Maxwell se plantean sobre un dominio circular de gran

extensión alrededor del objeto, en el contorno exterior del cual se considera una
condición de contorno absorbente, que permite una salida libre de la onda. Se trata éste
de un contorno artificial, que no existe físicamente en la realidad, pero que debe
emplearse numéricamente para limitar la extensión del dominio analizado.

ΓPEC

Γabs

Figura 17. Representación esquemática del primer

ejemplo de electromagnetismo analizado

En este problema, y también en el otro ejemplo de Maxwell que se explicará más

adelante, el tratamiento de los flujos se hace mediante el planteamiento de interpolación
de flujos nodales (véase apartado 3.2).

Las figuras 18 y 19 muestran las soluciones obtenidas con la formulación x-y

para las distintas componentes del vector de incógnitas de los dos subproblemas
independientes (transverso eléctrico y transverso magnético) en que se descompone el
problema global de electromagnetismo planteado en dos dimensiones (ecuaciones 10 y
11, respectivamente).

37
 5. Aplicaciones y resultados

Comp. 1 (εE1) Comp. 2 (εE2)

Comp. 3 (µH3)

Figura 18. Solución numérica con la formulación x-y para las distintas
componentes del vector de incógnitas del problema del transverso
eléctrico para el ejemplo de un obstáculo circular (p = 10)

38
 5. Aplicaciones y resultados

Comp. 1 (µH1) Comp. 2 (µH2)

Comp. 3 (εE3)

Figura 19. Solución numérica con la formulación x-y para las distintas
componentes del vector de incógnitas del problema del transverso
magnético para el ejemplo de un obstáculo circular (p = 12)

39
 5. Aplicaciones y resultados

Por otra parte, la figura 20 muestra, para cada uno de los dos problemas
anteriores, las representaciones gráficas de la RCS obtenidas numérica y analíticamente.
Se aprecia como, en ambos casos, la calidad de las soluciones obtenidas es muy alta,
proporcionando unos valores para la RCS que coinciden prácticamente de manera
exacta con la solución analítica.

Analítica
10 GD(XY)

0
RCS

-5

-10

-15

-20

-pi -3pi/4 -pi/2 -pi/4 0 pi/4 pi/2 3pi/4 pi

φ
a)

16 Analítica
GD(XY)

14

12

10
RCS

-pi -3pi/4 -pi/2 -pi/4 0 pi/4 pi/2 3pi/4 pi

φ
b)

Figura 20. Representación gráfica de las soluciones numérica y analítica

para el valor de la RCS en el ejemplo de un obstáculo circular:
a) transverso eléctrico (p=10); b) transverso magnético (p=12)

40
 5. Aplicaciones y resultados

Finalmente, con el fin de comparar el comportamiento de las dos formulaciones

contrastadas en esta tesina, en la figura 21 se muestran, para uno de los dos problemas
resueltos (transverso magnético), los dos gráficos de convergencia en p (en términos del
grado de interpolación y del número de grados de libertad) correspondientes al error
cometido en el cálculo de la RCS. Nuevamente, se ha adoptado un criterio de tolerancia
para detener el número de cálculos en cada formulación, exigiendo un logaritmo del
error inferior a -2.5.
Convergencia en p Convergencia en p
0.5 0.5

0.0 0.0

-0.5 -0.5

log (RCS error L2)

log (RCS error L2)

-1.0 -1.0

-1.5 -1.5

-2.0 -2.0

-2.5 -2.5

-3.0 -3.0

-3.5 -3.5
0 2 4 6 8 10 12 14 16 0 2000 4000 6000 8000 10000
p ndof

xy Isoparamétrica xy Isoparamétrica

Figura 21. Gráficos de convergencia en p en la resolución del problema transverso para

el ejemplo de un obstáculo circular con formulaciones isoparamétrica y x-y

Se puede apreciar que, de manera análoga al problema de convección estudiado

previamente, la formulación x-y presenta mejores resultados que la formulación
isoparamétrica, especialmente para altos grados de interpolación. La tolerancia
requerida se ha alcanzado sobradamente con la formulación x-y con una interpolación
de grado 12, mientras que la formulación isoparamétrica ha necesitado llegar a una
interpolación de grado 15 o, equivalentemente, ha necesitado utilizar aproximadamente
unos 3000 nodos de más. Para p=12, incluso, la diferencia entre el error de ambas
formulaciones alcanza aproximadamente un orden de magnitud.

Por otra parte, destaca también en la figura 21 el comportamiento un tanto

irregular del error en las interpolaciones de más bajo orden. El fenómeno puede
atribuirse claramente a la imposibilidad de captar, con una combinación de malla
grosera como la utilizada y interpolaciones tan bajas, las importantes ondas que presenta
la solución. Los resultados obtenidos en estos casos no gozan de un mínimo de calidad
que permita tenerlos en consideración. No obstante, es interesante su inclusión en la
gráfica para mostrar, una vez más, que las diferencias entre las dos formulaciones con
interpolaciones de bajo orden no son muy remarcables.

41
 5. Aplicaciones y resultados

Tras la resolución del sencillo y típico caso del obstáculo circular, se ha

considerado a continuación la resolución de un segundo ejemplo de electromagnetismo,
también bidimensional, pero con geometrías algo más complejas. En este caso, el
obstáculo interceptado por las ondas representa, una sección de un ala de avión.
Nuevamente, la onda se propaga en dirección horizontal de izquierda a derecha, y su
longitud característica es ahora de 0.5. Las condiciones de contorno se prescriben de
manera análoga al caso anterior (figura 17).

Las figuras 22 y 23 muestran las soluciones obtenidas para las distintas

componentes de los vectores de incógnitas de los problemas del transverso eléctrico y
magnético, respectivamente, para este nuevo caso.

Comp. 1 (εE1) Comp. 2 (εE2)

Comp. 3 (µH3)

Figura 22. Solución numérica con la formulación x-y para las distintas
componentes del vector de incógnitas del problema del transverso
eléctrico para el ejemplo de la ala de avión (p = 7)

42
 5. Aplicaciones y resultados

Comp. 1 (µH1) Comp. 2 (µH2)

Comp. 3 (εE3)

Figura 23. Solución numérica con la formulación x-y para las distintas
componentes del vector de incógnitas del problema del transverso
magnético para el ejemplo de la ala de avión (p= 7)

Por otra parte, la figura 24 muestra las representaciones gráficas de la RCS

correspondientes cada uno de los dos problemas anteriores. En este caso, en ausencia de
una solución analítica exacta, se representan únicamente los valores obtenidos
numéricamente.

43
 5. Aplicaciones y resultados

-10 GD(XY)

-12

-14

-16
RCS

-18

-20

-22

-24

-26

-pi -3pi/4 -pi/2 -pi/4 0 pi/4 pi/2 3pi/4 pi

φ
a)

GD(XY)
4

0
RCS

-2

-4

-6

-8

-pi -3pi/4 -pi/2 -pi/4 0 pi/4 pi/2 3pi/4 pi

φ
b)

Figura 24. Representación gráfica de la solucion numérica

para el valor de la RCS en el ejemplo de una ala de avión:
a) transverso eléctrico; b) transverso magnético

44
 5. Aplicaciones y resultados

Los resultados obtenidos son coherentes con los publicados por P.D. Ledger
[13]. Se confirma por lo tanto la aplicabilidad del método implementado para la
resolución de este tipo de problemas, posibilitando el análisis de formas todavía más
complejas e interesantes, como por ejemplo aviones enteros.

5.2 Euler

Las ecuaciones de Euler, formuladas para la descripción y análisis de los flujos

compresibles, tienen aplicaciones múltiples y variadas, y conciernen áreas tan diversas
como la industria, las ciencias ambientales o incluso la medicina [1]. De hecho, los
sectores de la industria aeroespacial y de la automoción son hoy en día los principales
promotores de la investigación en este campo. El interés en un mejor conocimiento
sobre el comportamiento de los fluidos compresibles es muy grande, puesto que es esta
la base de la mayoría de los avances del mundo de la aerodinámica.

A partir de la resolución de las ecuaciones de Euler es posible conocer todas las

propiedades características del fluido en cualquier punto del dominio estudiado. Es
precisamente en base estos datos, como por ejemplo los referentes a la presión,
velocidad o temperatura alrededor de un obstáculo, como se diseñan la mayoría de
innovaciones para la mejora del comportamiento aerodinámico de aviones, naves o
vehículos.

En el estudio de las ecuaciones de Euler aparece nuevamente toda la

problemática inherente a la resolución numérica de cualquier problema hiperbólico de
tipo vectorial, como el estudiado en el apartado anterior. En este caso, además, para
condiciones de flujo supersónicas, existe una dificultad añadida: la posible existencia de
discontinuidades en la solución debido al carácter no lineal del sistema hiperbólico que
las ecuaciones constituyen (para mayor detalle, véase apartado 2.3). No obstante, en esta
tesina los problemas de flujo supersónico no se han considerado.

En esta tesina se plantea la resolución de un sencillo ejemplo de las ecuaciones

de Euler, que se utiliza con frecuencia como test de validación de códigos para flujo
compresible. Se trata de la simulación del comportamiento de un fluido ideal que
atraviesa un obstáculo circular en condiciones subsónicas. Pese a ser un caso no natural,
el ejemplo considerado es de gran utilidad desde el punto de vista computacional,
puesto que, como que las ecuaciones de Euler son para fluidos no viscosos, la
resolución de este caso permite evaluar la difusión numérica introducida por el método.

El planteamiento del problema considerado es, en ciertos aspectos, similar al

primer caso de electromagnetismo estudiado. Las ecuaciones de Euler se imponen en un
dominio circular de gran extensión alrededor de un obstáculo circular. En el contorno
del obstáculo se aplica una condición de contorno de tipo muro, que en este caso
significa una restricción de no penetración. Por otra parte, el contorno exterior del
dominio estudiado se divide en dos partes: en la mitad del lado izquierdo se imponen las
condiciones de entrada de flujo, mientras que en el lado derecho la condición es de
salida del flujo.

45
 5. Aplicaciones y resultados

Γmuro

Γentrada Γsalida

Figura 25. Representación esquemática

ejemplo de flujo compresible analizado

En la resolución de este problema, y a diferencia de los otros casos estudiados

hasta el momento, el tratamiento de los términos de los vectores de flujos se ha hecho
mediante un planteamiento de full-quadrature. De hecho, algunos autores afirman que,
por el contrario, la adopción de un planteamiento de interpolación de flujos nodales es
en este caso inviable, puesto que los resultados que con él se obtienen son del todo
insatisfactorios [14]. Un riguroso tratamiento –mediante la full-quadrature- de los
términos de los flujos resulta pues determinante en este problema.

El problema considerado es de tipo transitorio. Sin embargo, tras un cierto

intervalo de tiempo dado, la solución adquiere un carácter estacionario. Generalmente
es esta la solución que más interesa, puesto que es la que corresponde a unas
condiciones ya estabilizadas del flujo. En la resolución del problema, la detección de
este tiempo preciso en que la solución puede considerarse estacionaria se ha efectuado a
través de un criterio de parada.

El problema previamente descrito se ha resuelto, con distintas mallas, para

ciertos grados de interpolación determinados. Debido al limitado tiempo disponible para
la ejecución de los cálculos, se ha optado por unos grados de interpolación bajos, que
requieren un menor coste computacional. Sin embargo, esta elección tiene una
desventaja asociada: puesto que las diferencias entre las dos formulaciones contrastadas
en la tesina sólo son importantes a partir de interpolaciones de orden relativamente alto,
los resultados correspondientes a una formulación u otra en este caso serán muy
similares, perdiendo gran interés su comparación. Es por este motivo que la exposición
resultados se centrará primordialmente en la presentación y descripción de las
soluciones, incluyendo únicamente de manera puntual algún contraste de resultados
proporcione información realmente interesante.

La figura 26 muestra el aspecto de las soluciones correspondientes a las distintas

variables principales que aparecen en el vector de incógnitas del problema (densidad,
componentes de la velocidad y energía), obtenidas con una de las resoluciones de
problema mediante la formulación x-y. La solución corresponde a una interpolación con
p=2.

46
 5. Aplicaciones y resultados

Densidad Energía

velocidad x velocidad y

Figura 26. Representación de la solución numérica obtenida para las distintas

variables principales del vector de incógnitas para el problema de flujo
compresible subsónico con un obstáculo circular (p=2)

presión número de Mach

Figura 27. Representación de la solución numérica obtenida para las

variables de presión y número de Mach en el problema de flujo
compresible subsónico con un obstáculo circular (p=2)

47
 5. Aplicaciones y resultados

Por otra parte, la figura 27 muestra la solución correspondiente a otras variables

características del problema, que, si bien no son componentes del vector de incógnitas,
se pueden calcular a partir de éste y tienen un significado tanto o más importante que las
variables anteriormente mostradas. En particular, se muestran las soluciones que
corresponden a la presión y el número de Mach.

Pese a que las ecuaciones de Maxwell no disponen de solución analítica, existen

varias maneras de comprobar la validez de los resultados obtenidos. A partir de
derivaciones teóricas, basadas precisamente en la conservación de las distintas
variables, se conoce el valor que deben adoptar ciertas cantidades de interés, calculables
a partir de las variables principales del problema. Es el caso, por ejemplo, de la entropía
o de la pérdida total de presión, que ha de tener un valor nulo y unidad,
respectivamente, en todo el contorno del objeto, o bien del coeficiente de presión en
superficie, del cual es posible hallar una solución analítica.

Error en la entropía Error en la pérdida de presión

0.012 0.012
x-y x-y
0.01 Isoparamétrica 0.01 Isoparamétrica

0.008 0.008
error
error

0.006 0.006

0.004 0.004

0.002 0.002

0 0
-0.5 -0.3 -0.1 0.1 0.3 0.5 -0.5 -0.3 -0.1 0.1 0.3 0.5
x x

Error en el coeficiente de presión

0.25
x-y
Isoparamétrica Normas L2 del error
0.2
x-y Isoparamétrica

0.15 Entropía 0.0011 0.0012

error

Pérdida
0.0034 0.0036
0.1 de presión
Coeficiente
0.1045 0.1091
0.05 de presión

0
-0.5 -0.3 -0.1 0.1 0.3 0.5
x

Figura 28. Comparativa de errores para distintas cantidades de interés

cometidos con formulación x-y e isoparamétrica en la resolución del
problema de flujo compresible subsónico con obstáculo circular (p=2).

48
 5. Aplicaciones y resultados

La figura 28 muestra una representación gráfica del valor del error obtenido para
las tres variables anteriormente descritas a lo largo de todo el contorno del obstáculo
circular. En este caso, con el fin de mostrar una comparación entre las formulaciones
contrastadas en esta tesina, se han añadido también los datos correspondientes a las
mismas variables obtenidas con la formulación isoparamétrica. En la tabla que
acompaña los gráficos se muestra, de cara a una valoración más global, el valor de la
norma L2 en todo el contorno del objeto para cada uno de los errores graficados.

La observación de estos resultados constata las pocas diferencias que existe en

los distintos errores calculados para la resolución de este caso con un grado de
interpolación tan bajo (p=2). Si bien en algunos puntos concretos del dominio sí existen
diferencias algo apreciables entre las dos formulaciones, la media de los errores para
todo el contorno del objeto, computada tanto en la norma L2 como en la norma del
máximo, es muy similar en los dos casos, aunque siempre ligeramente menor en la
solución correspondiente a la formulación x-y.

Finalmente, un último aspecto que merece ser remarcado en referencia a las

distintas resoluciones del problema efectuadas es que, aún utilizando mallas de gran
refinamiento, en ninguno de los casos en que se ha utilizado una interpolación lineal
(p=1) ha sido posible obtener una solución que reprodujera de manera correcta el
comportamiento del fluido. Como prueba de ello se muestra, en la figura 29, solución
correspondiente al número de Mach obtenido en uno de estos casos.

Figura 29. Representación de la solución numérica obtenida

para el número de Mach en el problema de flujo compresible
subsónico con un obstáculo circular con p=1

Como se observa, la solución parece reproducir perfectamente el

comportamiento del fluido por delante y en los contornos superior e inferior del
obstáculo. Sin embargo, la solución que se obtiene para la parte trasera del obstáculo,
que es la zona de cálculo más problemática, no coincide en absoluto con la que debería
aparecer. Este fenómeno, que impide la correcta captación de la solución detrás del
objeto con interpolaciones lineales, no es un hecho exclusivo de la formulación x-y sino
que es algo característico del método de Galerkin discontinuo [15].

49
 5. Aplicaciones y resultados

A partir de los resultados mostrados previamente, se puede considerar que la

formulación x-y ha presentado un comportamiento satisfactorio en la resolución del
típico y clásico test de validación de códigos de flujo compresible. A continuación, y
tras resolver a modo de constatación algún otro ejemplo subsónico de geometría más
flujo subsónico, el siguiente paso lógico y natural que debería plantearse es la
resolución de problemas de flujo supersónico, con las nuevas dificultades que implicaría
la correcta captación de choques. Sin embargo, la resolución de este otro tipo de
problemas, con sus nuevas problemáticas asociadas, queda ya fuera del alcance de esta
tesina.

50
 6. Conclusiones

6 CONCLUSIONES

A lo largo de esta tesina se ha pretendido probar que la formulación

isoparamétrica habitualmente utilizada en la mayoría de métodos numéricos presenta
ciertas carencias conceptuales, que pueden repercutir en la calidad de los resultados que
proporciona. Para ello, se ha implementado una formulación alternativa, la formulación
x-y, que, desde el punto de vista de la rigurosidad matemática, plantea la resolución del
problema evitando las incorrecciones de la formulación isoparamétrica.

El planteamiento descrito se ha particularizado para el estudio concreto del

método de Galerkin discontinuo aplicado a la resolución de distintos problemas
hiperbólicos de dificultad creciente, hasta llegar a las ecuaciones de flujo compresible
de Euler. No obstante, los resultados y conclusiones que de él se derivan son totalmente
extensibles a otros métodos y problemas, puesto que las limitaciones observadas en la
formulación isoparamétrica no están asociadas a ninguno de estos dos aspectos.

Tras justificar en un primer momento la necesidad de una nueva formulación

que evite los problemas de la transformación isoparamétrica, y tras describir también las
modificaciones necesarias para que esto ocurra, la nueva formulación x-y implementada
en esta tesina ha sido testeada en la resolución de distintos tests clásicos con el objetivo
de validar su aplicabilidad.

Los resultados obtenidos en la resolución de los distintos problemas planteados

han demostrado un comportamiento satisfactorio de la nueva formulación
implementada. El código ha proporcionado soluciones de calidad para los distintos
problemas resueltos, demostrando tratarse de un método válido para la resolución de
todo tipo de problemas en general.

En comparación la formulación isoparamétrica, la formulación x-y ha

demostrado presentar, en general, un mejor comportamiento. No obstante, las
diferencias entre una y otra metodología varían en función de los parámetros de cálculo
utilizados. Así, se ha comprobado que es especialmente para interpolaciones de alto
orden que la calidad de los resultados proporcionados por una y otra metodología
presenta una mayor diferencia. La formulación x-y obtiene en estos casos unos errores
considerablemente menores que la isoparamétrica, habiéndose detectado diferencias de
hasta un orden de magnitud. Sin embargo, para grados de interpolación bajos, la mejora
de resultados proporcionada por la nueva formulación es mucho menos significativa.

El hecho de que el estudio se haya planteado con el método de Galerkin

discontinuo ha resultado especialmente interesante para la obtención de estos resultados
sobre la convergencia en p. Puesto que el método de Galerkin discontinuo facilita en
gran medida el trabajo con interpolaciones de alto orden, su utilización en esta tesina ha
permitido emplear en algunas resoluciones grados de interpolación especialmente altos.

51
 6. Conclusiones

Es precisamente en estos casos donde se han observado las mayores diferencias

previamente descritas.

Por otra parte, se ha detectado también que el tamaño de los elementos de la

malla de discretización es otro factor que condiciona las diferencias de comportamiento
entre las dos metodologías. En particular, cuando los elementos de la malla son grandes,
y por lo tanto sus lados contienen grandes curvaturas, los resultados de la formulación
x-y son remarcablemente mejores que los obtenidos con la formulación isoparamétrica.
Sin embargo, a medida que la malla se va refinando y el tamaño de los elementos
disminuye, las diferencias son cada vez menores. El hecho puede atribuirse a que, al
tratarse de elementos más pequeños, la curvatura de los lados es cada vez menor, y por
lo tanto, el error de la formulación isoparamétrica asociado a las geometrías curvas
disminuye.

Desde el punto de vista computacional, las diferencias que existen entre una y
otra metodología no son especialmente significativas. Si bien es cierto que las
modificaciones que implica la formulación x-y para un correcto cálculo de los elementos
curvos tienen un coste considerable, hay que remarcar que sólo es necesario introducir
los cambios este tipo particular de elementos. En general, el número de elementos con
lados curvos en una malla es muy inferior al de elementos rectos, de manera que, en
términos globales, la diferencia de coste computacional entre ambas metodologías no es
apenas remarcable.

La formulación x-y implementada en esta tesina ha sido aplicada, siguiendo un

orden de dificultad creciente, a los problemas de convección lineal pura,
electromagnetismo y flujo compresible. No obstante, en este último caso únicamente se
han considerado unas condiciones de flujo subsónico. Una continuación de este estudio
debería extender la aplicabilidad del método a la resolución también de problemas de
flujo supersónico, superando las nuevas dificultades que ello entraña, como es la
correcta captación de choques mediante el uso de viscosidad artificial.

Finalmente, y en línea con la estrategia de disminución de errores que se ha

perseguido en esta tesina, conviene mencionar otra línea de investigación que puede
reportar resultados realmente interesantes. Se trata del estudio y desarrollo de nuevos
métodos numéricos de elementos finitos basados en el uso de modelos geométricos
exactos, independientes de las mallas de discretización. Según sugieren algunos autores,
un planteamiento como el descrito evita no solo el error debido a la transformación
isoparamétrica estudiado en esta tesina, sino también otras carencias asociadas a la mala
definición geométrica del dominio en estudio propias de toda discretización física en
elementos [16].

52
 7. Referencias bibliográficas

7 REFERENCIAS BIBLIOGRÁFICAS

[1] J. Donea and A. Huerta: Finite Element Methods for Flow Problems, Wiley,
2003
.
[2] R. J. LeVeque: Numerical Methods for Conservation Laws, Birkhäuser, 1990.

[3] W. H. Reed and T. R. Hill: Triangular mesh methods for the neutron transport
equation, Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory,
1973.

[4] B. Cockburn: Discontinuous Galerkin methods for Computational Fluid

Dynamics, in: E. Stein, R. de Borst, T.J.R. Hughes (Eds.), Encyclopedia of
Computational Mechanics, Vol. 3, Chapter 4, Wiley, 2004.

[5] J. Bonet and R.D. Wood: Nonlinear Continuum Mechanics for Finite Element
Analysis, Cambridge, 1997.

[6] A. Taflove: Computational Electrodynamics: TheFinite-Difference Time-

Domain Method, Artech House, Inc., 1995.

[7] C. Hirsch: Numerical computation of internal and external flows, Vol. 2, Wiley,
1990.

[8] E.F. Toro: Riemann solvers and numerical methods for fluid dynamics, Springer,
1997.

[9] P.-O. Persson and J. Peraire: An Efficient Low Memory Implicit DG Algortihm
for Time Dependent Problems, Proc. of the 44th AIAA Aerospace Sciences
Meeting and Exhibit, 2006.

[10] B. Cockburn and C.-W. Shu: Runge-Kutta discontinuous Galerkin methods for
convection-dominated problems, Journal of Scientific Computing, Vol. 16, No.
3, 2001.

[11] B. Szabó and I. Babuška: Finite Element Analysis, Wiley, 1991

[12] Q. Chen and I. Babuška: Approximate optimal points for polynomial

interpolation of real functions in an interval and in a triangle, Computer
Methods in Applied Mechanics and Engineering, Vol 128, No. 3-4, 1995, 405-
417.

[13] P. D. Ledger, K. Morgan, O. Hassan: Frequency and time domain

electromagnetic scattering simulations employing higher order edge elements,
Computer Methods in Applied Mechanics and Engineering, Vol. 194, No. 2-5,
2005, 105-125.

53
 7. Referencias bibliográficas

[14] P.G. Koen Hillewaert, N. Chevaugeon, J.-F. Remacle: Hierarchic multigrid

iteration strategy for the discontinuous Galerkin solution of the steady Euler
equations, International Journal for Numerical Methods in Fluids, Vol. 51, 2006,
1157-1176

[15] F. Bassi and S. Rebay: High-order accurate discontinuous finite element solution
of the 2D Euler equations, Journal of Computationa Physics, Vol. 138, 1997,
251-285

[16] A. Huerta, R. Sevilla and S. Fernández-Méndez: NURBS-Enhanced Finite

Element Method (NEFEM), Proc. of the 14th International Conference on Finite
Elements (FEF), 2007.

54