Introducción a los métodos

numéricos
en mecánica de los fluidos
computacional
 Antes de iniciar el estudio de los fundamentos de los métodos numéricos,
observemos que las ecuaciones de gobierno siempre describen la forma:

Término temporal Término convectivo Término difusivo Término fuente

A está ecuación se la denomina: “Ecuación de tranporte genérica” porque representa

el transporte de una propiedad cualquiera ɸ.
Para aplicar recuperar a partir de la ecuación anterior se debe reemplazar la propiedad
intensiva correspondiente a la extensiva transportada. En el caso de la masa,
se tiene:
Clasificación matemática de las ecuaciones de gobierno
Observemos por un momento a la ecuación de cantidad de movimiento

Esta es una ecuación en derivadas parciales (PDE, por sus siglas en inglés).

Se dice, de manera general, que es una PDE dado que hay más de una variable
Independiente involucrada.

Las derivadas involucradas ya no son respecto a una única variable.

Ejemplos EDP’s:

Ordinaria u=u(t)

Ordinaria u=u(t)

PDE u=u(x,y)

PDE u=u(x,t), c puede ser constante, o no.

Orden de una PDE

Es el orden de la mayor derivada involucrada en la PDE

PDE u=u(x,t) de segundo orden en el tiempo y el espacio

PDE u=u(x,y) de segundo orden en el espacio.

Linealidad

A cada lado del igual hay expresiones algebraicas que involucran a la función incógnita
(dependiente), sus derivadas, otras funciones de la variable independiente y/o constantes.

Se dice que una ecuación es lineal si los términos de dicha ecuación, que involucran a la
función dependiente, contienen solo productos de dicha función o sus derivadas con otras
funciones conocidas de la variable independiente o constantes.

En una PDE lineal la función dependiente no aparece elevada a ninguna potencia mayor
que 1 o es el argumento de cualquier función no lineal (seno, coseno, exp, log, etc).

ODE lineal
De manera general se tiene:

Para que esta ecuación sea de segundo orden ninguno de los términos que involucran
derivadas segundas puede ser nulo.

Ecuación homogénea: se dice que es homogénea si la función nula es solución.

u=0 es solución, por lo tanto es homogénea.

Una importante propiedad de las ecuaciones lineales es que dadas dos soluciones a la
ecuación u y v, y dos constantes ɑ y β la combinación lineal:

también será solución a la ecuación original.

Por ejemplo si consideramos la ecuación

Dadas dos soluciones, dado que

 Se tiene entonces que:

w también es entonces solución.

Es decir, la combinación lineal de soluciones también es solución.

La linealidad en una función implica que se satisfagan dos condiciones:

1) Aditividad

f(x+y) = f(x) + f(y)

2) Homogeneidad de grado 1

f(ax) = af(x) para todo a

Ejemplo. Indique si la función f(x) = 3x +2 es lineal

Probemos la propiedad 1: f(a+b) = f(a) + f(b)

f(a+b) = 3a+2 + 3b+2 (Satisface la aditividad)

C
2) probando la homogeneidad de orden 1 f(ax)=af(x)

f(ax) = 3ax + 2a

af(x) = a(3x+2)

Dado que: 3ax + 2a = a(3x+2), entonces se satisface hp1, por lo que se puede afirmar
que la función es lineal.

Considere ahora la función f(x) = sen(x) y verifique si es lineal.

Probemos entonces la propiedad de aditividad:

f(a+b) = f(a) + f(b)

f(a+b) = sin(a+b)
sin(a+b) ≠ sin(a) + sin(b)
f(a)+ f(b) = sin(a) + sin(b)

Claramente la igualdad anterior no satisface.

El operador derivada es lineal, lo que implica que:

Ejercicio 1. Indique si lo siguientes operadores son lineales o no

 Ecuación no homogénea: se dice que es no homogénea si la solución trivial nula no es
solución.

La solución u=0 no satisface a la ecuación

Ecuación no homogénea: se dice que es no homogénea si la solución trivial nula no es

Solución.

Ecuación no lineal:

La ecuación es claramente no lineal porque la variable dependiente

está elevada al cuadrado.
Otra manera de verificar la no linealidad de la anterior ecuación es considerando que
Si la ecuación fuera lineal entonces deberían ser lineal y homogénea en el sentido
de que la solución nula debe ser solución.

Suponga que u(t) y v(t) son dos soluciones no triviales:

La linealidad implica que w = u + v sea también solución.

w no satisface la ecuación por lo tanto no es lineal

 Ejercicio. Indique cuales de las siguientes ecuaciones son ordinarias, PDE, lineales y
no lineales.
g)
a)

b) h)

i)
c)

j)
d)

k)
e)

f) l)
EDPs cuasi-lineares
Considerando el sistema dado por:

Se considera:

A) lineal con coeficientes constantes si A y B están compuestas de constantes.

B) lineal con coeficientes variables si A y B son funciones de las variables

independientes.

C) SI B depende linealmente de U entonces el sistema es aún lineal

D) Quasi-lineales

Las derivadas de alto orden aparecen de manera lineal, esto quiere decir que no
aparecen multiplicadas entre sí, aparecen solo multiplicadas por coeficientes que son
función de las variables dependientes A=A(U).
Ec. de Euler Estas son las ecuaciones de Euler.

Gobiernan el flujo inviscido

Aparecen solo derivadas de orden

uno y todas de forma lineal
Clasificación de las EDP quasi-lineales

Para resolver una EDP (o un sistema de EDPS) se debe definir un dominio

sobre el cual se obtendrá su solución.

El dominio en este contexto no es otra cosa que la región en la cual se

buscará la solución.

SI la EDP tiene como variable independiente solo a x, el dominio será una línea.

Si tiene como variables independientes a x y t, el dominio seguirá siendo una

linea pero la distribución solución sobre esa linea cambiará con t.

SI la EDP tiene como variables independientes a x, y, entonces el dominio será

una superficie.

Si la EDP tiene como variables independientes a x, y, z, entonces el dominio

será una región en el espacio (volumen)
Con los fines de fijar los conceptos considere que está abordando la
solución a la ecuación:

En este caso la ecuación será resuelta en la región definida desde x=a

hasta x=b para tiempos mayores que cero.

Nótese que el dominio en este caso es una línea que se extiende de a

hasta b, por lo que en a y en b se deben definir además condiciones de
borde apropiadas.

Las condiciones de borde estableceran el comportamiento de la función

solución en los bordes.

En general las condiciones de borde pueden ser de diversos tipos, las más
simples son del tipo: Dririchlet (o esencial) y Neumann (o natural).
 Una condición del tipo Dirichlet implica que en el borde se conoce el valor
de la variable dependiente, es decir en este caso sería:

u(a) = k1 y u(b) = k2 (donde k1 y k2, son valores conocidos)

Por otro lado, la condición del tipo Neumann implica que en el borde se
conoce el valor de la derivada de la variable dependiente

du/dx (a) = k2 y du/dx(b) = k2

También se pueden imponer otros tipos de condiciones, como la de Robin

(o mixta) que impone en el borde una combinación entre la condición
del tipo Newmann la del tipo Dirichlet:

Aquí a y b, son factores de ponderación que se definen en función de la condición particular.

Adicional a las condiciones de borde, como en este caso la ecuación también involucra
al tiempo se debe proporcional una condición inicial.

La condición inicial, no es otra cosa que una solución de partida, es decir los valores
de la variable dependiente sobre el dominio de interés en t=0.

Cabe mencionar que la imposición de condiciones de borde e iniciales condicionan a la

evolución del problema, si se imponen condiciones de borde y/o iniciales no
consistentes se tendrá un problema “mal condicionado” lo que comprometerá que sea
factible arribar a una solución.

Teniendo el conjunto de ecuaciones a resolver, el dominio, las condiciones de borde e

iniciales ya se estaría en condición de abordar la solución.

Previo a esto haremos otra clasificación de las ecuaciones de acuerdo al tipo de

problemas que representan:

1) Problemas de avance

2) Problemas de equilibrio
Problemas de equilibrio (Ecuaciones elípticas)

Situaciones en estado estacionario

Flujo estacionario incompresible e irrotacional

Conducción de calor estacionario

Problemas de valores en el contorno

Cualquier perturbación en el interior

del dominio cambia a toda la solución.

Soluciones suaves

Solución en un punto influenciada por los alrededores.

Problemas de avance (marching problems)

Situaciones en estado no estacionario (transferencia de calor dependiente del

tiempo, cualquier flujo inestacionario, transporte de ondas, etc)

La ecuaciones pueden ser hiperbólicas o parabólicas.

No todos los problemas de avance son no estacionarios.

Parabólicas

Problemas no estacionarios que involucran difusión

Flujos viscosos no estacionarios y transferencia de calor no estacionaria.

Problemas de contorno y valor incial

Cualquier perturbación en el interior del dominio solo afectará eventos posteriores.

La solución avanza en el tiempo y se “difunde en el espacio”

Solución siempre suave para t>0

A medida que el tiempo avanza la solución tiende al estado estacionario.

Ecuaciones hiperbólicas

Vibraciones

Disipación despreciable

También acarrean problemas de contorno y valor inicial

Los problemas de flujo compresibles (inviscidos) con velocidades cercanas o
superiores a la del sonido son hiperbólicos.

La aparición de ondas (discontinuidades) es una manifestación de está naturaleza

matemática.

Por esto los métodos que se desarrollan para flujos compresibles con altas
velocidades deben ser capaces de “lidiar” con dichas discontinuidades.

Un perturbación solo afectará a cierta región de la solución.

En problemas hiperbólicos la información se propaga con velocidad finita.

En problemas parabólicos y elípticos la velocidad de propagación es infinita.

 Métodos de solución
●
Diferencias finitas
●
Elementos finitos
●
Volúmenes finitos
●
Flujo potencial
Métodos espectrales
●
Métodos sin mallas
●
Lattice boltzmann
●
Entre otros
Método de las diferencias finitas
●
Es el método más antiguo (Euler, 1768)

●
Requiere la discretización del dominio en puntos (nodos) en
los que se obtendrá una solución aproximada al problema
analizado.

●
Se utiliza la expansión en serie de Taylor para determinar
aproximaciones a los operadores diferenciales de las
ecuaciones.
●
Luego de la aproximación se obtienen ecuaciones para los
valores nodales de las incógnitas
●
Puede ser aplicado a cualquier tipo de malla, pero es más
usual en el contexto de mallas estructuradas.
●
Se pueden usar incluso mallas con distribución no regular, no
obstante hay un límite en la no uniformidad para mantener la
precisión.
●
El orden dependerá de la forma de la aproximación en
diferencias.
●
El método no es conservativo por construcción por lo que se
debe ser muy cuidadoso cuando se deriva o selecciona la
aproximación.
 Para iniciar el proceso se debe partir de la discretización geométrica del dominio

1D
Cada nodo se identifica con
Uno i (1D) o dos i,j (2D) índices.

2D

Los nodos vecinos se definen

Por medio del incremento
o reducción del índice respectivo
Se puede definir un paralelo entre la solución de la versión continua de una ecuación
y una solución aproximada discreta.

En la solución continua (si está existe se obtiene una expresión cerrada que determina
a toda la solución en cualquier punto sobre el dominio físico

Por ejemplo la ecuación:

Posee la solución exacta:

Solución que es válida para

Por el contrario una solución discreta ofrece únicamente valores de la solución en puntos
específicos que dependen de la discretización adoptada.
Aproximación de las derivadas
No obstante cuando se aproxima una solución se debe acotar el dominio, ya que no es
factible tomar puntos infinitos por tanto se puede buscar una manera de aproximar de
manera discreta al operador diferencial y obtener la solución a u en cada valor nodal.

De la definición de derivada sabemos que:

También sabemos que a partir de la expansión en serie de Taylor se puede aproximar el

valor de una función f en la vecidad de un punto a

Si reemplazamos a por un valor nodal xi:

Se ha aproximado a la derivada en el
punto xl utilizando la pendiente de la
curva secante que pasa por los puntos
x(l-1) y x(l+1).

Al hacer este tipo de aproximaciones

se introduce error.

El error se puede reducir reduciendo

la separación del intervalo, pero nunca
desaparecerá.
Ejemplo. Expansión en serie de Taylor

Obtenga la expansión en serie de Taylor de la función f(x) = ln(x) alrededor de x=1. Retenga
los tres primeros términos de la expansión.

Elabore una gráfica que compare la expansión al aumentar de termino en término para
f(x) en 0<=x<=1.5

Se tiene entonces:
La expansión a determinar implica evaluar hasta la tercera derivada de la función:

Reemplazando en la expansión en serie:

Ejercicio. Dada la función f(x)= xex encuentre la aproximación alrededor de x=0 reteniendo
5 términos.

Elabore una gráfica para las aproximaciones con 1, 2, 3, 4 y 5 y compare con el valor exacto
para -0.5 =<x <=0.5
De la expansión en serie es posible despejar la derivada del orden necesario en el punto
de interés, por ejemplo la de primer orden:
Término de orden superior

Término de orden n
Término de segundo orden

Origen del error de truncamiento

Al reescribir la anterior aproximación de la primera derivada se tiene:

Es evidente que el error de truncamiento en este caso está dado por:

Lo que en general se puede escribir como:

donde los factores de ponderación α dependen del tipo de esquema. Se debe notar en este
punto que en el esquema centrado m=1. Este factor m indica el orden de la aproximación
conseguida
Si se analiza la anterior expresión se observa que el término que domina al error
será en el que el tamaño de la discretización (Δx) este elevado a la m. x) este elevado a la m.

Se observa también que a medida que Δx→0 la importancia de los términos x→0 la importancia de los términos
asociados al error será menor, es decir que la aproximación hacia el valor exacto
mejorará.

El orden de la aproximación indica la rapidez con la que el error se reduce a medida

que el tamaño de la discretización se reduce.

Para los esquemas descentrados hacia adelante y atrás (que son de primer orden el
error cae con Δx→0 la importancia de los términos x) el error de truncamiento está dado por:
Si se ignoran los términos de mayor orden a 1 se puede entonces aproximar la derivada
en xi como:

Ahora, si introducimos la notación de Leveque para la aproximación

Aproximación hacia adelante

Aproximación hacia atrás

 Otra posibilidad es aproximar la derivada como un promedio de las dos aproximaciones
anteriores:

Aproximación centrada

Si se lleva a cabo el análisis del error en este caso se observa que:

De donde se concluye que en este caso m=2 y por lo tanto la aproximación es de segundo
orden.

Otra nomenclatura que se usa es decir que el esquema es de O(Δx→0 la importancia de los términos xm) (orden delta x a la m),
Lo que significa que error está decayendo a medida que se refina con esa tasa.
Es muy importante tener en cuenta que el orden no es una medida absoluta del error.
La figura muestra, mediante la pendiente de la recta tangente, como la aproximación de
centrada (segundo orden) es una mejora las hacia adelante y atrás (primer orden).
Por supuesto es posible obtener aproximaciones de mayor orden:

Aproximación de tercer orden, el error es proporcional a h³

Dada la función u=sin(x) y el punto xi=1 se desea aproximar du/dx en la vecindad de xi.

En este caso conocemos la derivada exacta du/dx = -cos(x), por lo tanto tenemos un
valor para comparar que tan buena es la aproximación

du/dx (xi) = 0.5403023

Utilizando las tres aproximaciones presentadas

Se seleccionan diferentes pasos (h), se realizan las aproximaciones y se calcula el
error:
Valor exacto

Valor aproximado
Ejercicio: realice la aproximación para la derivada anterior utilizando la fórmula de
3er orden:
Aproximación de la segunda derivada I

Segundas derivadas aparecen en los términos difusivos, para aproximarla se puede

usar la derivada numérica de la aproximación a la primera derivada:

Geométricamente la segunda derivada es la recta tangente a la curva que representa

a la primera derivada
Aproximación de la segunda derivada II

Por ejemplo partiendo de la aproximación de la primera derivada hacia adelante se

tiene:
Aproximación de la segunda derivada III

Usando la aproximación hacia adelante para las derivadas internas

Aproximación de la segunda derivada IV

Ejercicio: Derive una aproximación para la segunda derivada usando aproximación

hacia adelante para las derivada externa y hacia atrás para las derivadas internas

Otra aproximación para la segunda derivada puede ser centrada, en cuyo caso
tomando una molécula computacional:
Aproximación de la segunda derivada V

Con la molécula mostrada:

Aproximando la segunda derivada:

Aproximación de la segunda derivada VI

Se tiene entonces:

Ejercicio. Muestre que para una malla equidistante la anterior aproximación se

puede escribir
Aproximación de la segunda derivada VII

Otra forma de obtener aproximaciones para la segunda derivada es por medio de la

serie de Taylor
Cantidad de puntos requeridos

Existe una relación entre la cantidad de puntos necesarios para obtener una aproximación
De un dado orden para una derivada de un orden k.

Por ejemplo, para una derivada de primer orden al menos se requieren dos puntos, ya que
a través de dichos puntos para el polinomio de menor orden (recta) para el cual la primera
derivada no es nula.

Si se aplica el mismo razonamiento llegamos a que para la derivada de segundo orden

la mínima cantidad de puntos es de 3 puntos. Con lo que se tiene:

donde N es el número de puntos, p la precisión del método y k el orden de la derivada a

aproximar.
El enfoque semi-discreto

Como se debe intuir en este punto, si se utilizan operadores discretos para las derivadas
temporales y espaciales involucradas se puede reducir el problema de PDEs a un
sistema algebraico.

En ocasiones es ventajoso considerar las discrtizaciones de los operadores espaciales

y temporales de manera separada. En primer lugar se discretiza el espacio con lo cual
el sistema original de PDEs se reduce a uno de ODES.

Este formulación se conoce como formulación semi-discreta. A continuación cuando se

discretiza al término temporal el sistema queda completamente discretizado, para ello
se puede utilizar cualquier esquema apropiado para realizar el avance en el tiempo.
Cabe mencionar que algunos métodos no conducen a formas semi-discretas ya que
abordan a la parte espacial y temporal en simultaneo.

Sin embargo, es muy importante notar que la gran mayoría de algoritmos y métodos
utilizan el enfoque semi-discreto, lo cual ofrece la posibilidad de utilizar operadores
discretos distintos en espacio y tiempo.

Al separar las discretizaciones espacial y temporal se obtiene una compresión más clara
del impacto en cuanto a precisión y estabilidad del esquema utilizado en el espacio
y del empleado en el tiempo.

Esta metodología también ofrece el beneficio de poder utilizar la amplia teoría existen sobre
ODEs.
Operadores en diferencias en forma matricial
Considérese el siguiente operador discreto que aproxima a la segunda derivada
involucrada en la siguiente ecuación en el punto j:

Considere ahora que dicho operador se utilizará en un dominio , el

cual se ha discretizado con 4 puntos interiores:

Supongamos también que en este dominio se han impuesto condiciones del tipo
Dirichlet en los bordes u(0)=ua y u(π)=ub, si se aplica el esquema introducido
previamente se tiene:
Note que hay una ecuación por nodo, y que en este caso se dice que el stencil utiliza tres
puntos dado que el operador discreto resultante en cada punto solo involucra a los dos
puntos vecinos.
Por ejemplo, si consideramos la ecuación resultante de aplicar el operador en el punto 2:

Se observa que solo involucra información del los puntos o nodos vecinos.

Stencil o molécula computacional es el término que se usa para designar los puntos
(información) utilizada en la construcción del operador discreto.

Dado que siempre se manejan gran cantidad de puntos es conveniente usar herramientas
del álgebra matricial, así si se introduce la siguiente terminología:
Al introducir dichos operadores el sistema resultante se puede escribir como:
Discretización temporal
Como se comentó el enfoque semi-discreto permite utilizar un esquema para
discretizar la parte espacial y otra para la temporal. Lo cual como se comento
introduce ciertas ventajas.

En problemas no estacionarios se deben tener en cuenta 4 coordenadas, aparte de

las 3 asociadas al espacio se agrega al tiempo.

Esta coordenada se puede discretizar en el contexto de las diferencias finitas tal como
se hacía en con el espacio.
Para discretizar a los términos temporales se pueden aplicar diferentes esquemas de
cálculo que como en el espacio tendrán diferente orden de aproximación dependiendo de
como se aproxime la derivada.

Dado el problema de valor inicial:

Condición inicial

Se busca obtener la solución u para tiempos mayores que t0.

Para ello la solución debe existir y ser única. Para que sea única, se requiere que dicha
solución sea continúa en el sentido de Lipschitz.
A diferencia de un problema de valores en el contorno (caso del espacio), donde se
debían dar condiciones de borde para resolver el problema. En el caso de los
problemas de valor inicial se debe proporcionar un valor inicial sobre todo el dominio
de cálculo.

A partir de está condición inicial se evalúan los estados posteriores del problema, es
decir la solución que se busca en tiempos posteriores discretos:

Donde el tiempo posterior depende del paso de tiempo Δx→0 la importancia de los términos t utilizado.

Para diferenciar de cuando se evalúa un operador diferencial discreto utilizado en el

espacio de uno empleado en el tiempo, para los temporales utilizaremos super-
índices:
SI se integra al problema de valor inicial dado por:

Se tiene:

El valor en n+1 dependerá de como se evalúe la integral. Lo cual dependerá de la cuadratura

que se utilice
Método de Euler explícito (fordward Euler)

Método de Euler implícito (backward Euler)

Punto medio (leapfrog)

Trapecio (Crank-Nicolson)
Es importante notar que cualquiera de los métodos presentados requieren de conocer el
valor del campo a evaluar en un punto de partida, es decir una condición inicial o valor
previo.

No obstante, algunos utilizan en el miembro del lado derecho información del pasado y
otros información del tiempo a calcular:
Evaluado en tiempo pasado

Evaluado en el tiempo en
que se va a calcular

Los esquemas que emplean meramente información de pasado para calcular en n+1
se denominan explícitos, los que el miembro del lado derecho debe ser evaluado en
n+1 se conocen como implícitos.
Los esquemas de discretización presentados producen buenos resultados con paso de
tiempo cortos.

No obstante, se debe analizar como se comportan a medida que el paso de tiempo se

incrementa. Ya que en problemas asociados a la fluidodinámica se encuentran
involucradas diversas escalas temporalres, lo que da origen al problema conocido
como la rigidez numérica (stiffness).

La rigidez impone restricciones sobre el paso de tiempo que se puede utilizar y da

origen a lo que se conoce como estabilidad numérica que será discutida
posteriormente.

Se debe destacar que los métodos explícitos son menos restrictivos en cuanto a uso de
memoria y son más fáciles de programar.
Ejemplo. Método de Euler explícito aplicado al problema de la pelota que cae

Como se vio en clase al aplicar la segunda ley de Newton a una pelota que cae por la
atmósfera se arriba a la ecuación diferencial dada por:
Sabemos que en tiempo t=0 la velocidad es V=0 (esta sería la condición inicial). Con lo
que el problema queda especificado como:

ODE

CI
Para aplicar el método de Euler hacia adelante a este problema reconoce que:
Si calculamos desde t=0 hasta t=5 s, usando un paso de tiempo de Δx→0 la importancia de los términos t=0.1.

Además considerando que para una esfera el arrastre se puede obtener a

partir de la siguiente figura:
Con fines prácticos supondremos que el Cd=0.41 (y en el rango considerado es constante).
La densidad del aire la supondremos estándar ρ=1.225 kg/m³) y la esfera tendrá diámetro
D=50 cm y m=1000 g de masa.

Con los datos anteriores se tiene que k = 0.054119 kg/m.

Para obtener la solución se aplica el esquema:

Con

Al discretizar el tiempo con los datos t0=0 s hasta t = 5.0 s con Δx→0 la importancia de los términos t=0.1. se tiene:

t = [0.0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60 0.70 0.80 0.90 1.00 1.10 1.20
1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00……,5]

t(n+1) = t(n) + dt
Se puede entonces evaluar la solución en n+1:

Evaluación en t=0.1

0.981 m/s

V(0.1)= 0.981 s Velocidad en t=0.1 s

Se puede entonces evaluar la solución en n+2 (t=0.1+0.1=0.2:

Es decir la velocidad en t = 0.2 s es de ~1,957 m/s m/s

Este esquema de cálculo se puede hacer recursivo utilizando cualquier lenguaje de
programación, por ejemplo en matlab se tiene (sphereEulerV1.m)
Se debe recordar que de acuerdo con la expansión en serie de Taylor, el método de Euler
hacia adelante o explícito es solo de primer orden.

Por otro lado el cálculo de la velocidad terminal arroja:

La gráfica muestra la solución
numérica en azul junto con la
velocidad terminal (roja)
14.07 m/s
Se observa como a medida que el
tiempo avanza la solución numérica
tiende al valor de la velocidad
terminal
Velocidad terminal (14.10 m/s)
Cabe mencionar que la solución
Analítica para este problema se
puede encontrar mediante integración
y está dada por :
A continuación se presenta una comparación de los resultados obtenidos
empleando el método de Euler junto con la solución analítica para el intervalo
t[0,5]
Es importante notar que aunque la opción de utilizar el esquema de Euler puede parecer la
más elemental, esta es suficiente para diversas aplicaciones.

Para realizar el avance en el tiempo existen diversos tipos de esquemas, dentro de los más
utilizados se encuentra la familia de métodos de Runge-Kutta. Estos métodos con la finalidad
de mejorar la aproximación emplean puntos intermedios entre n y n+1
Ecuación de convección lineal

La ecuación de convección lineal es una de las más simples que se pueden analizar:

Como sabemos para resolver cualquier EDP se deben imponer condiciones de borde
iniciales. De acuerdo a las BC que se apliquen se tienen dos fenómenos distintos:

1. En uno en un contorno se define una onda que ingresa al dominio, es decir se tiene
una condición de entrada. En el otro extremo no se impone ningún tipo de condición. Puesto que
la ecuación es de primer orden esto es consistente con su grado, sólo se requiere de una CB.
Este fenómeno se conoce como convección simple.

2. En el otro el flujo simulado es periódico. Para cualquier instante lo que entra por un extremo
debe salir por el otro. Este problema se conoce como de bi-convección. Este tipo de CB
periódicas resultan de utilidad en el análisis de métodos de solución.
Ecuación de convección lineal II

Para obtener su solución se puede aplicar el método separación de variables, de esta

forma se propone que la solución es la combinación lineal de dos funciones
de la forma:

Reemplazando en la ecuación original:

Ecuación de convección lineal III

El lado izquierdo es función únicamente de t, mientras el derecho lo es de x:

B es función de t y A de x, por lo tanto para mantener la igualdad la única opción es que

las derivadas sean constantes.
Ecuación de convección lineal III

Reemplazando en la solución propuesta u=A(x)B(T)

Esta ecuación representa la propagación de una onda con velocidad c. La forma de

la onda no cambia en el tiempo solo se desplaza dependiendo de c.
Ecuación de difusión lineal

La ecuación de difusión es otro modelo muy importante en el análisis de métodos de solución:

Donde ν es una constante positiva. Representa los flujos generados por el

movimiento molecular en un medio continuo.

Esta ecuación si u es la temperatura y ν la difusividad térmica, no es otra cosa que

la ecuación que gobierna la difusión de calor es una dimensión.

Admite condiciones de borde periódicas, de Dirichlet, de Neumman o mixtas.

La ecuación de difusión tiene solución no trivial estacionaria, es decir admite

solución sin la derivada temporal:
Ecuación de difusión lineal II

Ejercicio: Investigue la solución al problema estacionario sobre un dominio x(0,L) con

condiciones de borde u(0)=T1 y u(L)= T2

Otras soluciones estacionarias se pueden obtener al añadir un término fuente a la

ecuación:

En 2D la ecuación de difusión se escribe: (parabólica)

Añadiendo un término fuente se puede obtener de nuevo una solución estacionaria:

Ejercicio. Verifique la clasificación dada a las

(elíptica)
EDP mostradas.
Propiedades de los esquemas numéricos

Consistencia

La discretización debe tender a la solución continua a medida que el tamaño de la malla

tiende a cero.

La diferencia entre la ecuación discreta y la exacta define el error de truncamiento

Se puede evaluar usando expansión en serie de Taylor

Un esquema será consistente si el error de truncamiento tiende a cero cuando Δx→0 la importancia de los términos x→0 y
Δx→0 la importancia de los términos t→0.

El error de truncamiento es proporcional a una potencia del espaciamiento de la malla

Consistencia
La discretización debe tender a la versión continua a medida que el tamaño de la malla
tiende a cero.

La diferencia entre la ecuación discreta y la exacta define el error de truncamiento

Se puede evaluar usando expansión en serie de Taylor

Un esquema será consistente si el error de truncamiento tiende a cero cuando Δx→0 la importancia de los términos x→0 y
Δx→0 la importancia de los términos t→0.

El error de truncamiento es proporcional a una potencia del espaciamiento de la malla,

si el término más importante proporcional a (Δ x) (o en el tiempo) se dice que el método
n

es de orden n.

Idealmente los términos de orden n deben discretizarse con esquemas del mismo
orden.

Sin embargo en ocasiones resulta ventajoso discretizar ciertos términos con mejor
orden de aproximación que otros dependiendo de la aplicación.
Consistencia

Algunos métodos tienen errores de truncamiento que son función de la

relación Δx→0 la importancia de los términos x/Δx→0 la importancia de los términos t (o viceversa).

Que una aproximación sea consistente no implica que la solución al

sistema discreta se convierta en la exacta a medida que h→0, para que
esto suceda el método también debe ser estable
Si se considera a la ecuación de difusión y se le aplica el esquema FTCS se tiene:

Luego si se expande en serie de Taylor alrededor del punto (xi,tn) se tiene:

Lo términos resaltados corresponden a los errores de truncamiento de las aproximaciones

utilizadas.
 Lo anterior se puede escribir:

Estos errores se reducen a medida que el tamaño de la disscretización Δx→0 la importancia de los términos x y el paso de
tiempo Δx→0 la importancia de los términos t se hagan más chicos, siempre y cuando las derivadas de orden superior se mantengan
acotadas.

Lo que indica que el esquema FTCS es consistente.

Estabilidad

Se dice que un método es estable si este, a medida que el proceso

evoluciona no tiende a amplificar los errores presentes.

En problemas no estacionarios, si la solución exacta está acotada, un

método estable producirá una solución acotada.

En un proceso iterativo un método estable será aquel que no diverja.

En problemas con BC complejas y no linealidades la estabilidad se

analiza basándose en la experiencia y en la intuición.
Estabilidad

Es una propiedad de interés únicamente en problemas de avance

(hiperbólicos y parablólicos). El método más utilizado es el método de
von Neumann

La estabilidad es compleja de analizar, especialmente cuando están

involucradas condiciones de borde y no linealidades.

Algunos métodos tienen limitaciones en el paso de tiempo que se

puede emplear, a menos que se use sub-relajación.

Un método estable hará que los errores de redondeo introducidos a

medida que la simulación avanza se atenúen o al menos se mantengan
acotados.
Estabilidad

Los esquemas se pueden clasificar en:

1) incondicionalmente estables

2) condicionalmente estables

3) incondicionalmente inestables
Número de Courant-Frederich-Lewy (CFL)

Se define como la relación entre el intervalo de tiempo y el tiempo de

residencia en una celda computacional.

Considere la ecuación de advección lineal discretizada con el esquema

FTCE:
Si se aplica el método de von Neumann para analizar la estabilidad (ver Krause, cap.
3)
Se demuestra que el esquema FTFS será estable sí y solo sí:

Considere ahora la discretización con el esquema FTBS

La condición es opuesta a la del esquema FTFS.

Si se aplica el esquema FTCS:

El análisis de von Newmann no arroja ninguna condición de estabilidad, por lo que

este esquema al ser aplicado a la ecuación de advección lineal es
incondicionalmente inestable.
En resumen, los esquemas FTFS y FTBS son condicionalmente estables mientras el
FTCS es incondicionalmente inestable al ser aplicado a la ecuación de advección
lineal.

Las condiciones de estabilidad de los esquemas FTFS y FTBS, como se observó son
opuestas. Es decir que el esquema será estable o inestable dependiendo de la
velocidad de advección.

De acuerdo con lo anterior, FTFS será estable si a<0, mientras que FTBS lo será
parra a>0.

La estabilidad, como se vio, exige que:

Esto implica que, los esquemas descentrados serán estables si la dirección de

propagación de la información, dada por el signo de la velocidad de la onda, es
opuesta a la dirección en la cual se realiza la aproximación.

Visto de otra forma, cuando el cambio dx desde el nodo con coordenada i se

haga en la dirección opuesta a la “corriente”, es decir corriente arriba (aguas
arriba) en la dirección upwind.
La condición de estabilidad está directamente relacionada con la naturaleza matemática
de las ecuaciones.

La condición de estabilidad exige entonces que en el diseño de un esquema de solución se

deba tener en cuenta la dirección de propagación de la onda.

Como se vio aparte de la condición upwind, en el caso de los esquema FTFS y FTBS
Requieren que se satisfaga:

Esta relación fue introducida por Courant, Fredrich y Lewy en 1928 y se conoce como
Condición CFL o número de Courant.

Como adelantó está es una de las condiciones más importantes en CFD. Esta condición
permite establecer el paso de tiempo máximo que se podrá utilizar con un esquema numérico
con una discretización espacial determinada.
El número de Courant, que no es otra cosa que el cociente entre el paso de tiempo
y el tiempo característico de la convección.

El tiempo característico para la convección, es el tiempo necesario para que una

perturbación producida en el campo de flujo se propague gracias a este mecanismo una
distancia Δx→0 la importancia de los términos x.
Convergencia

Un método es convergente si la solución discreta estable tiende a la exacta a

medida que el tamaño de la malla se reduce.

Para IVP lineales el teorema de equivamencia de Lax establece que: “

Dado un problema bien condicionado y una aproximación en diferencias
finitas que es consistente, la estabilidad es condición necesaria y suficiente
para que la solución converja”

En problemas no lineales, influenciados fuertemente por la condiciones de

borde, la estabilidad y la convergencia son difíciles de analizar. La
convergencia se analiza a partir de experimentos numéricos: por ejemplo, se
realizan cálculos con mallas cada vez más chicas.
Convergencia

Si el método es estable y las aproximaciones consistentes, la solución se

debe hacer independiente de la malla.

Para mallas lo suficientemente pequeñas la velocidad de convergencia es

dominada por el principal componente de error de truncamiento.
Conservación

Las leyes de gobierno son leyes de conservación.

Los métodos utilizados deben satisfacer la conservación a nivel discreto a

nivel local y global.

En estado estacionario y en ausencia de fuentes o sumideros, la cantidad de

la variable conservada que entra a un VC debe ser igual a la que sale.

En el método de los volúmenes finitos si se emplea la forma “fuerte” de las

leyes de conservación la conservación se garantiza tanto para cada volumen
finito como para todo el dominio de cálculo.

En otros métodos se debe ser cuidadoso con la selección de los esquemas

para garantizar la conservación.
Conservación

En presencia de fuentes o sumideros en el dominio, se debe asegurar que su

tratamiento numérico sea consistente de tal forma que se asegure que el
flujo neto a través de los contornos del dominio se equilibre con el efecto de
dichas fuentes.

Esta es una importante propiedad ya que impone restricciones sobre el error.

Al garantizar la conservación (masa, cantidad de movimiento, energía, etc) el

error sólo sufrirá una distribución diferente sobre el dominio.

Un esquema que no garantice la conservación producirá fuentes artificiales,

cambiando el balance local y global.

Esquemas no conservativos, pero consistentes y estables en el límite de

mallas pequeñas pueden conducir a soluciones correctas.
Conservación

Los errores debidos a la no conservación se aprecian en mallas “gruesas”

Se deben preferir en lo posible esquemas conservativos

Acotación (Boundedness)

En muchas ocasiones la solución debe estar contenida entre ciertos límites.

Propiedades que físicamente no pueden ser negativas deben mantenerse en

el rango posible, por ejemplo, la densidad, la energía cinética de la
turbulencia. Otras variables deben estar contenidas, entre el rango 0 y 100,%
por ejemplo la concentración de una especie química no puede superar el
100%. En la ecuación de calor, en ausencia de fuentes, los valores mínimo y
máximo de la temperatura deben estar contenidos dentro de los valores del
contorno.

Esta es una propiedad difícil de garantizar, sólo algunos de los esquemas de

primer orden la aseguran.

Los esquemas de alto orden dan origen a soluciones no acotadas. Esto sólo
sucede en mallas gruesas.
Acotación (Boundedness)

La presencia de “overshoots” y “undershoots” en la solución, son indicios de

que el error en la solución es grande y de que se requiere refinamiento.

Los esquemas que producen soluciones no acotadas, pueden ser estables y

convergentes, pero deben ser evitados en lo posible.
.
Factibilidad (Fiabilidad)

No es propiamente una propiedad de los esquemas, no obstante es

importante.

Relacionada con el conjunto de hipótesis y simplificaciones que se efectúan

para modelar un problema.

Si se parte de hipótesis incorrectas es probable que no se consiga

convergencia incluso cuando el método la garantice.
Precisión

No hay que perder de vista que las soluciones numéricas son solo
aproximaciones.

Además de los errores que se introducen en las condiciones de borde, en la

programación, etc. Los métodos de solución siempre sufren de tres tipos de
error sistemáticos:

Error de modelado: Dado por la diferencia entre el fenómeno fı́sico real que
se desea resol ver y el que se resuelve aproximadamente. Es introducido por
las diferentes hipótesis que se formulan.

Errores de discretización: Error introducido por la diferencia entre la

solución exacta a las ecuaciones en su forma continua y la aproximada de
las ecuaciones discretas.
Precisión

Error de convergencia iterativa: Se define como la diferencia que existe

entre la solución iterativa dada al sistema de ecuaciones y la solución exacta
del mismo.
Análisis del error
Si se conoce una solución de referencia el error absoluto se puede definir como:

E(u)=ua-ue

No obstante los métodos de aproximación no producen soluciones continuas, sino

discretas definidas en puntos nodales i, de aquí que:

Para una malla con N elementos se podría medir el error utilizando la norma-1:

No obstante esta medida puede dar una mala apreciación del error, dado que medirá
el error como la simple suma aritmética de los errores en cada punto, por lo que que
si la cantidad de puntos aumenta esta magnitud podría aumentar aunque el error real
en algunos lugares caiga.
Análisis del error

Para remediar este comportamiento se puede ponderar cada error nodal por
el tamaño de la celda, lo que en el caso de una malla uniforme conduce a:

La métrica anterior mide al error en términos de la norma L1, en general el

error se puede medir en términos de la p-norma, definida por:

En general existen diversas maneras de estimar el error.

Análisis del error
Una pregunta recurrente es, cómo debo medir el error?

Esto dependerá de si conoce o no una solución de referencia para comparar. Si la

conoce podrá aplicar una norma apropiada y computar su valor, pero esto por sí solo
no es suficiente también deberá comparar con una malla refinada y observar si el error
describe un comportamiento asintótico

Malla 50x50 Malla 100 x 100

Análisis del error

[Link]

NASA_2

Error2
Análisis del error
En ausencia de la solución de referencia se pueden comparar varias soluciones
obtenidas con mallas cada vez más finas.

No obstante hay que tener cuidado, ya que obtener un comportamiento convergente no

garantiza que la solución es correcta, solo indica que se está arribando a una función
solución que satisface a las ecuaciones y a las condiciones de borde.

Es perfectamente posible que un código este mostrando una buena convergencia a una
solución incorrecta.

Por lo tanto, es de trascendente importancia que el analista esté familiarizado con el

problema con el que trata para discernir si su solución tiene sentido o no.
 Proceso de solución
Se desea resolver la ecuación de difusión lineal en un dominio 1D:

Sujeta a las CI y CB:

Para resolver a la ecuación lineal numéricamente por el método de las diferencias

finitas primero se debe discretizar al dominio de cómputo, en el espacio y en el
tiempo.

Para una malla uniforme en el espacio se tiene:

En el tiempo también se debe discretizar:

El objetivo es obtener la evolución de la variable dependiente en el tiempo y su

distrubución en el espacio.
Para esta primera aproximación, para la derivada de orden 2 se usa la aproximación centrada:

Mientras que para la de primer orden se usa una aproximación hacia adelante

Luego de reemplazar en la ecuación:

Este esquema se conoce como FTCS (forward in time, centered in space), despejando se tiene:

Hay que recordar que en cualquier caso el cómputo parte de valores de la variable dependiente
conocidos en el tiempo n.

El esquema obtenido permite obtener la solución para cada punto interno de la malla
con la que se ha discretizado el dominio.

Dependiendo del tipo de CB puede ser necesario resolver ecuaciones adicionales para los
bordes.

Si en un punto correspondiente a un borde se fija el valor de la variable dependiente (dirichlet)

No hace falta resolver una ecuación adicional.

Por otro lado, si en un contorno se impone una condición del tipo Neumann (o mixta), donde
Interviene una derivada, sí será necesaria la solución de una ecuación adicional.
Ejercicio. Obtenga las aproximaciones FTFS (forward in time, forward in space) y FTBS (forward
in time, backward in space) para el problema de la ecuación de difusión lineal en 1D.

Nótese que la aproximación FTCS para obtener la solución en el tiempo n+1, utiliza únicamente
Datos del tiempo n, este tipo de formulación se conoce como explícita.

Si en lugar de utilizar en la aproximación de la derivada espacial datos en n, se emplean en

El tiempo posterior n+1, se tiene:

Se obtiene la formulación BTCS (backward in time, centered in space).

La formulación ahora involucra los valores nodales i-1, i e i+1 en el instante futuro n+1.
Valores que son incógnitas.
Por lo que ya no es posible obtener la solución en el punto i en el tiempo n+1 a partir de lo
que se conoce en el punto i en n. La formulación entonces ya no es explícita, y se
denomina entonces implícita.

En este escenario la solución implica que se debe resolver un sistema de ecuaciones

algebraícas acoplado de la forma:

La dimensión del vector solución será igual a la cantidad de nodos interiores

utilizados en la malla.

La matriz de coeficientes [A] es una matriz dispersa, ya que no estará completamente

llena, dado que el esquema de discretización utiliza para ensamblar cada ecuación nodal
Una pequeña cantidad de información en comparación a la cantidad total de nodos.

Si para numerar los nodos se es cuidadoso la matriz A será en banda, ya que en cada
ecuación nodal solo intervendrán los valores adyacentes.

El ancho de banda depende del orden del esquema utilizado.

 Para el esquema centrado de 2do orden la matriz es tridiagonal:

La estructura diagonal de la matriz de coeficientes es una gran ventaja, tanto desde el

punto de vista de las operaciones requeridas para resolver el sistema como
de la memoria necesaria para su almacenamiento.

Esto permite la utilización de eficiente métodos iterativos desarrollados especialmente

para resolver este tipo de problemas.

Debe notarse que la resolución de un sistema de ecuaciones algebraicas acopladas

siempre conlleva un costo computacional sustancialmente mayor que el
correspondiente a un conjunto de ecuaciones desacopladas.

El tipo de formulación utilizada (explícita o implícita) no tiene implicancias en cuanto al

orden de aproximación del método numérico.
Los esquemas de diferencias finitas implementados en este ejemplo son en ambos casos de
primer orden, ya que el mínimo orden de aproximación es O(∆t).

La ventaja de utilizar una u otra formulación está relacionada con otros aspectos de la
solución numérica. La resolución de esquemas explícitos es mucho más simple que la de
los esquemas implícitos, pero estos últimos en general presentan mucho mejor desempeño
en cuanto a la estabilidad, lo cual muchas veces permite “avanzar” más rápido en el
tiempo, tal como veremos más adelante.
Solucioń a la ecuación de convección lineal

Reemplazando estas aproximaciones de las derivadas en la ecuación:

Este esquema se conoce como FTCS (forward in time, centered in space)

De esta forma se obtuvo un esquema que permite determinar la solución en cada punto
De la malla computacional a la ecuación de convección lineal

No obstante, antes de utilizar este esquema se debe ver si es consistente y estable.

El esquema FTCS es consistente pero es incondicionalmente inestable.

Ejercicio. Investigue por qué el FTCS es incondicionalmente inestable

Debido a que el esquema solo utiliza información del paso de tiempo anterior j para
determinar la solución en j+1, se dice que es explícito.
Sabemos que la solución a la ecuación de convección está dada por:

Es decir, corresponde a una onda con forma determinada que se propaga con
velocidad c.
Ejemplo de aplicación: método diferencias finitas en 1D

La ecuación mostrada no es otra que la ecuación que gobierna la conducción de calor

en un sólido, es decir un proceso de difusión, donde ɑ es el coeficiente de difusión. Que
en el campo de la transferencia se conoce como conductividad térmica por unidad de
longitud (ɑ= 5 W / (m² K)).

Se desea determinar la distribución de temperatura en una barra de material

homogéneo que tiene una longitud total de 2 cm, se sabe que:

T(0) = 100 °C = 317.15 K

T(L) = 400 °C = 673.15 K

 Este problema de carácter estacionario cuenta con solución
analítica:
Sabemos que en x=0, T=373.15K
Y que en x=0.02 T=673.15K

Por lo que de:

Se deduce que C2=373.15

Mientras que de la 2da BC:

Evaluando en los
extremos del dominio
Con lo que la solución final resulta:

Nota: es buena práctica revisar que la solución respete las BC.

En este caso la distribución de

temperatura el lineal y va aumentando
de la zona de menor temperatura
hacia de la mayor
Deseamos ahora aproximar la solución mediante la técnica de las diferencias finitas.

Se debe entonces aproximar la segunda derivada

Usando la aproximación
centrada de segundo orden

Es decir la ecuación discreta para i-ésimo punto estará dada por:

Ecuación que involucra información únicamente
De los nodos inmediatamente vecinos a i
Ya tenemos entonces la ecuación en su versión discreta, ahora tenemos que discretizar
el dominio donde obtendremos la solución. Para este ejemplo utilizaremos 6 puntos,
es decir:

Note que hemos identificado de forma distinta a los nodos en el contorno y a los
del interior del dominio.

Note también que los valores de la variable dependiente en el contorno en este caso
Están dados por las condiciones de borde
Para cada nodo se debe satisfacer:

Para 6 puntos

Para determinar las coordenadas de lo nodos se puede escribir:

Iniciando en i=1

X2= 0 + 0.02/5

X3 = 0.02/5 + 0.06/5

X4 = (0.02/6) + 0.06/6 + 0.06/5

Y así sucesivamente, en matlab se puede hacer x=0:0.02/5:6 y se obtiene

Es decir, en este caso el espaciamiento de la malla o es igual a 0.004 m,
De esta forma el cociente

Evaluemos en cada nodo

Hay una ecuación por

cada nodo interior!!!!
 Los valores de T en el contorno del dominio
son conocidos.

El resto de valores son desconocidos en el

Interior del dominio.

Pasando al lado derecho los valores nodales que se conocen, es decir los del contorno:

Las incógnitas son los valores en el

Interior del dominio (en este caso)
Recordando que un sistema de ecuaciones se puede escribir en forma
matricial:

Entonces el sistema de ecuaciones nodales del problema se puede escribir

Matriz de coeficientes Vector solución Vector del lado derecho

Es evidente que el tamaño de la matriz de coeficientes, de los vectores solución y del lado
Derecho son función del tamaño de la malla.

Además los coeficientes de la matriz son función de la ecuación particular.

En este caso, los elementos de la diagonal son: Aii = -2, mientras que los que están por
Fuera de la diagonal Aij = 1

Para resolver un sistema de ecuaciones escrito en forma matricial existen diversos

métodos de básicamente dos familias: directos e iterativos

En matlab si se tiene la matriz A que este caso está dada por:

Y el vector b dado por:

Para determinar la solución se hace:

[A]b = c por lo tanto b = [A]⁻¹b ¹b

Lo que en matlab (o octave) resulta:

También se puede hacer inv(A)B en matlab o

A\B. Ejercicio. Pruebe estos métodos por su cuenta para
solucionar el sistema de ecuaciones de este problema.