Condiciones de contorno

Hola

Voy a escribir un post acerca de una parte fundamental en la resolución de ecuaciones diferenciales. Se trata de las condiciones de contorno.

Como ya sabréis, las condiciones de contorno permiten definir una ecuación diferencial a un caso concreto, esto es, la ecuación pasa de tener infinitas familias de soluciones a tener una sola. Las condiciones de contorno son definidas de esta forma porque típicamente describen el comportamiento de un sistema en la zona donde se acaba el sistema que se quiere analizar: por ejemplo, si quiero describir la transferencia de calor a través de una placa metálica, el volumen en el cual me interesa analizar como fluye el calor utilizando mi ecuación diferencial será una zona determinada que yo me trace, digamos por simplicidad, un rectángulo o cuadrado que abarque todo el espesor de la placa. Se introduce un dibujo para ilustrar el ejemplo:

En la imagen se muestra en linea sólida el contorno físico de nuestra placa metálica, y en linea discontinua el límite de nuestro "volumen de control" que vamos a considerar. Típicamente en la literarura se denomina al dominio que se va a analizar con una ecuación Ω, y a cada uno de las superficies que delimitan dicho volumen se las denomina con la letra Γ . En este caso, se ha denotado como Γ1, Γ2, Γ3 y Γ4 a las superficies que delimitan el contorno.

 En este punto uno se plantea la pregunta de cuantas condiciones de contorno debe imponer y que forma tienen que tener. Esta pregunta no tiene una respuesta inmediata: depende de que fenómeno estemos describiendo.

Las ecuaciones diferenciales de la forma $u_t = K \nabla^2 u$, con $\nabla^2 = \frac{\partial}{\partial x_i} $ son denominadas ecuaciones parabólicas. Este tipo de ecuaciones se caracteriza porque debemos definir una condición de contorno por cada derivada que aparezca para cada una de las variables que este derivada.

Por ejemplo, si tenemos $u_t = D(u_{xx} + u_{zz})$, entonces tendremos que definir 4 condiciones de contorno, puesto que tenemos 2 variables espaciales independientes ("x" y "z") derivadas 2 veces. Estas condiciones deben estar dadas de forma que cubran TODAS las fronteras que posee nuestro volumen. Además, la variable tiempo también es una variable independiente del problema y esta derivada una vez. Por lo tanto necesitamos una condición más para esta variable. Típicamente se denomina a las condiciones de contorno para la variable tiempo "condiciones iniciales", puesto que marcan el estado de un sistema en un punto del tiempo, normalmente en el estado inicial de la simulación (t=0). Ejemplos clásicos de este tipo de ecuaciones son la ley de Fourier del calor o la ley de Fick de transferencia de masa.

Las ecuaciones del tipo $u_{tt}=u_{xx}$ se denominan ecuaciones hiperbólicas. Notad que la diferencia con respecto al caso parabólico es que la variable tiempo esta derivada dos veces. Para la resolución de estas ecuaciones sólo son necesarias dos condiciones iniciales.

Las ecuaciones tipo $\nabla^2 u = f(u_i) $ son denominadas ecuaciones elípticas. Estas ecuaciones necesitan definir i condiciones de contorno, esto es, una condición de contorno para cada variable dimensional independiente. Notad que en este tipo de problemas no hay variable tiempo.

¿Que forma deben tener las condiciones de contorno?. En la literatura estan clasificadas en tres tipos:
- Condición de contorno tipo "Dirichlet": Son aquellas que dan el valor de una variable en una de las fronteras. Ej: $T(y=0)=T_0$

- Condición de contorno tipo "Newman": Son aquellas que dan el valor de LA DERIVADA de una variable en una de las fronteras. Ej: $\frac{dT}{dx}(x=0)=0$

- Condición de contorno tipo "Robin" o "mixta": Son aquellas que dan el el valor de LA SUMA del valor de una función y de su derivada. Ej: $T(x=0)+\frac{dT}{dx}(x=0)=0$

Aplicaremos todo lo explicado al ejemplo que se trataba anteriormente. En primer lugar, el problema es un problema bidimensional, puesto que el dibujo es 2D y no se especifica nada de 3 dimensiones. Al ser una ecuación de Fourier, la ecuación es parabólica, por lo que exige definir condiciones de contorno en Γ1, Γ2, Γ3 y Γ4 . Estableceremos el origen de coordenadas en el extremo inferior izquierdo para determinar las condiciones apropiadamente.

 $$T_t= K\left( T_{xx}+T{zz} \right)$$

Recordad que nuestro objetivo es hallar T(x,z,t)

Conocemos las temperaturas de la placa superior e inferior. Esto nos da dos condiciones de contorno:

$$T(x,0,t)=T_{cold}$$
$$T(x,H,t)=T_{hot}$$

siendo H el espesor de la placa.

Falta definir 2 condiciones de contorno para los bordes Γ3 y Γ4. No conocemos el valor de la temperatura en dichos contornos, pero si podemos suponer que si nuestra placa se extiende uniformemente en la dirección x, es de esperar que el calor fluya para cualquier x en dirección z. Esto se traduce en que para una misma altura z, el calor va a fluir de la misma forma sea cual sea el punto en el eje x. Eso implica que no existirá variación de la temperatura al movernos en el eje x:


$$\frac{\partial T}{\partial x}(0,z,t)=0$$
$$\frac{\partial T}{\partial x}(L,z,t)=0$$

siendo L la longitud del trozo de placa que hemos tomado.

Puesto que la ley de Fourier depende del tiempo, debemos establecer una condición inicial. Podemos suponer que la temperatura en todo el dominio $\Omega$ sea una temperatura $T_0$ que tiene que ser dada, la cual puede ser una constante o una función.

$$T(x,z,0)=T_0$$

Con esto el problema quedaría cerrado y podría ser resuelto.

Un detalle más. Se trata de la diferencia entre condición de contorno homogenea y condición de contorno inhomogenea. Una condición de contorno es homogenea cuando esta igualada a cero (ej: T(0,y)=0). En caso contrario se denomina inhomogenea (Ej: T(0,y)=f(x,y)).
Un punto imporante que debe ser entendido es que las condiciones de contorno inhomogeneas introducen UNA FUENTE DE EXCITACIÓN EN EL SISTEMA. Esto es, cuando se dice que la temperatura en una frontera toma un valor $$T_0$$, estoy asumiendo que estoy introduciendo calor al sistema en ese punto. Esto de un modo u otro va a suponer que el sistema RESPONDERÁ ante dicha entrada con un comportamiento dado, que en este ejemplo será que el perfil de temperaturas en el volumen de control será uno específico. Una condición de contorno homogenea NO introduce excitación en el problema.

Sobre la divergencia y el rotacional (II)

Hola a todos

Voy a realizar alguna aclaración más acerca de estos dos operadores, de gran importancia en el mundo de la ciencia y la ingeniería como la divergencia y el rotacional. La idea me ha surgido al realizar unos ejercicios de comprensión sobre el significado de estos, y pienso que será util para aquellos que tenga que tratar con ellos o bien durante su formación o bien en su trabajo.

Se comentó en el anterior post que la divergencia tiene que ver con un balance de flujo a través de una superficie. Quizá lo más riguroso sería acudir a su definición matemática estricta:

$$\nabla \cdot F = \frac{1}{\bigtriangleup V}\int_S \vec{F} \cdot \vec{n} \cdot dS $$

Como se puede observar, la divergencia se define como el flujo de campo neto que atraviesa una superficie asociada a un volumen determinado. Vamos a poner un ejemplo de cual sería la divergencia para un campo concreto. En concreto, utilizaremos un campo en coordenadas esféricas:

 $\vec{F}=r^2 \vec{r}$

Aplicando la definición de divergencia para un sistema de referencia de coordenadas esféricas, se puede obtener el valor exacto:

$\nabla \cdot \vec{F}=\frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \vec{F} \right) = 4r $

El resultado nos dice que la divergencia es siempre positiva en todo el campo, salvo en el centro que debe ser nula. Efectivamente, cuanto más nos acercamos las lineas de campo son cada vez menores.

Se puede analizar las zonas del campo en las cuales la divergencia es positiva, negativa o nula. Recordad que la divergencia tiene que ver con un balance en un volumen determinado.  Como se puede observar, la divergencia es siempre positiva en todo el campo. Conforme nos alejamos del origen de coordenadas, el flujo a través de las distintas superficies que englobarían volúmenes de tipo esférico (escogemos estos por ser un sistema de coordenadas esféricas) dentro del campo (en el plano 2D serían circunferencias concéntricas) cada vez es mayor.

Al representar gráficamente un detalle de este campo y pintar las lineas de campo a una distancia r y otra r +dr, se observa que efectivamente las lineas de campo son mayores al alejarnos del centro y por lo tanto mayores en la superficie S+dS que en la superficie S. La divergencia aplicada a este diferencial de volumen que abarca ambas superficies debe ser mayor que cero.

Pero este balance visual puede llevar a error algunas veces, por lo que hay que ser precavido a la hora de evaluar. Lo ilustraremos ahora realizando el mismo análisis con otro campo distinto:

$\vec{F}=\frac{1}{r^2} \vec{r}$

La representación del campo es:

Como se aprecia, las lineas de campo son decrecientes conforme nos alejamos del centro. Si hacemos un análisis visual basándonos en el ejemplo anterior, a priori podría pensarse que la divergencia de este campo es siempre menor que cero para cualquier r, haciéndose 0 en el infinito. Si ampliamos en imagen un detalle en una zona del campo:

Parece que las lines de campo son menores en S+dS que en S en todo el plano, por lo que el balance sería menor que cero. Sin embargo, esto no es así. Calculando la divergencia analíticamente:

$\nabla \cdot \vec{F}=\frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \vec{F} \right) = 0 $

El resultado del cálculo analítico arroja un resultado 0 para cualquier punto del campo. ¿Que es lo que ha fallado en nuestra deducción?

La respuesta a esta pregunta no es trivial. Si observamos de nuevo la definición de divergencia, se puede apreciar que es la integral del flujo superficial a través de una superficie. Al escoger un volumen como el del detalle anterior, se observa que el valor del campo es mayor en un punto de la superficie S que en una posición de la superficie S+dS (porque al tener mayor radio, el valor del campo es menor). En concreto la variación del campo es inversamente proporcional a $s^2$. Pero hay que destacar que LA SUPERFICIE A TRAVÉS DE LA QUE CIRCULA EL FLUJO DEL CAMPO TAMBIÉN ES MAYOR EN S+dS QUE EN S. Concretamente, es mayor proporcionalmente con el cuadrado de la distancia.

Se tiene en conjunto que, LA DISMINUCIÓN DEL VALOR DEL CAMPO AL AUMENTAR EL RADIO SE VE COMPENSADA CON EL AUMENTO DE SUPERFICIE AL AUMENTAR DICHO RADIO. Digamos que, aunque F es más intenso en la superficie S que en S+dS, la superficie que atravesa es menor que la de S+dS, compensando así los valores en ambos puntos. En este caso, esta compensación es exactamente igual, y por tanto el flujo a través de ambas superficies es exactamente el mismo. El balance neto entonces es exactamente cero.

Por lo tanto, para resumir:

- La divergencia de un campo vectorial es un escalar.
- La divergencia implica un BALANCE NETO de flujo en un volumen determinado a través de su superficie.
- Al calcular la divergencia en un volumen dado (por ejemplo un rectángulo, o un sector circular) hay que tener en cuenta tanto el valor de campo en cada superficie como el tamaño de las superficies que atraviesa. Recordar que es un flujo superficial.

El próximo post introduciré un ejemplo similar para el operador rotacional

Introducción a los elementos finitos

Hola a todos

Hoy voy a realizar una brevísima introducción sobre uno de los métodos numéricos de mayor uso dentro del mundo de la ciencia e ingeniería: el método de elementos finitos.

Para los que no estais familiarizados, el método de elementos finitos consiste en realizar una aproximación discreta de una función continua que sea solución de una ecuación o sistema de ecuaciones que detallan el comportamiento de un sistema físico o matemático.

Voy a ilustrar un ejemplo unidimensional muy intuitivo para ver cual es la mecánica que se sigue:

Supongamos que se desea resolver una ecuación como esta:

$ \partial_x (p \partial_x y)+ (q+\lambda \sigma) \partial_x y= 0$

Esta ecuación es una expresión general para el problema de Sturm-Liouville, pero en este artículo esto es intrascendental. Se puede ver como una ecuación diferencial de segundo orden con coeficientes variables.

El método de elementos finitos propone la división de la función "y" en n elementos geométricos. Los extremos de cada elemento donde se juntan 2 o más elementos se denominan nodos. Se va a considerar que la función "y" es conocida EXACTAMENTE en todos los nodos. El valor de la misma en los elementos debe ser interpolado. Para realizar esto, se propone para cada elemento una "función nodal" determinada: esta función puede tener diversas formas, pero siempre debe cumplir que sea 1 en el nodo al que corresponde dicha función y 0 en el resto. Puesto que se ha discretizado la función en n+1 nodos donde el valor de la función es conocida, se puede definir la función y como la suma de cada uno de las funciones nodales multiplicadas por el valor de la función en cada una de ellas:

$y=\Sigma_{n=0}^{n+1} \phi_n y_n$

Puesto que esta solución (aunque aproximada) es solución de la ecuación diferencial, podemos sustituirla en la ecuación dada:

$ \partial_x (p \partial_x (\Sigma_{n=0}^{n+1} \phi_n y_n))+ (q+\lambda \sigma) \partial_x (\Sigma_{n=0}^{n+1} \phi_n y_n)= E $

Notad que ahora la ecuación NO esta igualada a 0, lo cual es lógico puesto que hemos "aproximado" la solución real por una solución aproximada y discreta. Este error de aproximación viene dado por la función error E.

Se puede observar que, al hacer esta aproximación se ha podido definir EXACTAMENTE cual es la expresión de la función error de nuestra propuesta numérica. Ahora, sólo hay que variar la función error con un criterio adecuado para hacer tender la función error a 0. Si se consigue, la solución aproximada se acercará a la real. Para hacer esto se recurrirá a lo que se denomina "método de los pesos ponderados".

Este método consiste en hayar una serie de funciones peso tal que la integral de su producto escalar con respecto a la función error E sea cero.

$\int w(x) E dx = 0$

Dicho de otro modo, se quiere encontrar una serie de funciones ortogonales a la función error. La propiedad de ortogonalidad y las funciones peso son importantes: si se elige como función peso las mismas funciones de interpolación nodales que se han introducido en la solución aproximada tal que sean ortogonales a la función peso, entonces es inmediato observar que las funciones de interpolación hayadas harán que la función error E tenga un valor muy próximo a cero y por lo tanto nuestra aproximación sea válida como solución final.

El concepto es muy simple en su fundamento, pero su implementación computacional se muestra mucho más tediosa. Sabemos que la solución aproximada será tanto más exacta cuanto mayor número de nodos se consideren en la solución. Sin embargo, se observa que habrá que resolver un sistema de ecuaciones algebráico con un número de ecuaciones igual al número de nodos que se consideren (¡¡ habrá que cumplir la condición de ortogonalidad para cada función peso !!), aumentando por lo tanto el esfuerzo de computación.

En posts siguientes explicaré más en profundidad algún ejemplo y otras observaciones más detalladas acerca de este método.