Transferencia de calor por elementos finitos con R (I). Simulador

El pasado mes de enero la borrasca Filomena nos regaló en Madrid la mayor nevada que yo haya visto, y vivo aquí desde 1999. Este hecho unido a que por otros motivos llevaba un tiempo dándole al coco sobre formas de aislamiento en climas fríos, me inspiró definitivamente para tratar de hacer algunas simulaciones de transmisión del calor (clic sobre la imagen para verla en todo su esplendor).

En este artículo vamos a modelar la Ecuación del calor por elementos finitos valiéndonos de las facultades matriciales de R. Es un tipo de cálculos intensivos que en cualquier lenguaje no vectorizado requeriría compilar bucles anidados para lograr una velocidad de ejecución aceptable.

Esta ecuación en derivadas parciales rige la transferencia de energía inducida cuando cuerpos a diferente temperatura entran en contacto. En esa situación se produce un flujo de calor desde el cuerpo más caliente al más frío a un ritmo proporcional al gradiente térmico, flujo que solo cesa cuando se igualan las temperaturas. En coordenadas cartesianas:

La anterior ecuación solo tiene solución cerrada para geometrías muy básicas, pero podemos resolver estructuras arbitrarias complejas con aproximaciones de las derivadas parciales. Discretizando espacio y tiempo en elementos finitos, el cálculo de la temperatura en cualquier punto en un instante se puede obtener de forma sencilla a partir de la temperatura en el instante discreto inmediatamente anterior y de las condiciones térmicas de su entorno.

Existen métodos implícitos para formular estas aproximaciones que son incondicionalmente estables, resolviéndose con un sistema de ecuaciones que implica llegar a una matriz tridiagonal que debe invertirse para resolver el sistema. Sin embargo vamos a utilizar una formulación explícita mucho más simple, con el único coste de requerir una condición de estabilidad. Más sobre estas diferencias de métodos en 'Numerical Modeling Of Earth Systems'.

La siguiente ecuación explícita define, para el caso unidimensional, la temperatura en el instante discreto (n+1) para la localización también discreta i del espacio, en función de su temperatura y de las de celdas adyacentes en el instante anterior (n). El primer término modela la transferencia de calor por conducción y el segundo permite introducir fuentes internas de calor (hacer clic para verla con más detalle):

La fórmula se adapta a cambios de conductividad al pasar de un medio a otro, simplificándose en la forma estándar para material uniforme cuando se asume una conductividad (k) constante en el espacio. El término T_i+1 - 2T_i + T_i-1 podemos verlo como un filtro detector de gradientes de temperatura. Al sumarse a la temperatura local (ponderado por la difusividad térmica del material), tiende a reducir las diferencias respecto a celdas vecinas, garantizando siempre que los flujos de calor intercambiados entre celdas contiguas sean iguales a ambos lados de la frontera.

Extender el cálculo al caso 2D o 3D resulta trivial, ya que al ser los términos en derivadas parciales independientes para cada eje solo hay que añadir los sumandos correspondientes a las nuevas dimensiones. En nuestro simulador nos vamos a plantar en un desarrollo 2D ya que permite una buena flexibilidad para modelar estructuras, pero sin entrar en tiempos de cálculo excesivamente largos ni definiciones laboriosas del espacio simulado.

~~~

La condición de estabilidad a la que hacíamos referencia antes relaciona la velocidad de propagación del algoritmo (en cada iteración la propagación del calor no puede ser más rápida que una unidad espacial discreta por unidad de tiempo discreto, es decir Δx/Δt) y la difusividad (α) de los materiales participantes en la simulación. Concretando, como se explica en 'Stability of Finite Difference Schemes on the Diffusion Equation with Discontinuous Coefficients' es condición suficiente de estabilidad que se cumpla:

Otra forma de interpretar este umbral es como el valor máximo de corrección para que el algoritmo no llegue a oscilar. Si el término T_i+1 - 2T_i + T_i-1 se multiplica por una cantidad menor a 1/2 (caso 1D), la corrección de las temperaturas T_i > T_i+1 = T_i-1 nunca arrojará unas nuevas temperaturas donde T_i < T_i+1 = T_i-1. Es decir T_i se acercará a T_i+1 = T_i-1 en cada iteración pero sin llegar nunca a producirse un cruce que se traduzca en oscilaciones.

~~~

Con las herramientas matemáticas construidas vamos a enumerar los materiales que nos planteamos simular:

• 'solid': cualquier tipo de material sólido definido por su difusividad (α) y su conductividad (k) térmicas. Serán la base principal de nuestras simulaciones.
• 'source': se trata de sólidos a los que añadimos una fuente de calor interna (q). El parámetro admite valores negativos en cuyo caso estaríamos modelando un sumidero de calor.
• 'fluid': los fluidos se definirán igual que los sólidos por sus parámetros conductivos de difusividad (α) y conductividad (k). Sin embargo para respetar su comportamiento convectivo se acudirá a una aproximación que se explica más adelante.
• 'boundary': materiales sobre los que imponemos una condición de contorno de temperatura constante. Es decir, aunque intercambien calor con objetos vecinos nunca varían su temperatura (p.ej. aire exterior a una construcción).
• 'insulate': material aislante ideal. Se modela forzando a replicar en su periferia la temperatura de los objetos colindantes, bloqueando así todo intercambio de calor al anularse el gradiente térmico.
• 'isoflux': material que "deja pasar" el flujo de calor, ya sea entrante o saliente, existente en su periferia. Se modela forzando a replicar el gradiente térmico local existente en los objetos colindantes. NOTA: material no suficientemente probado.

Para definir la simulación de forma intuitiva y flexible, la escena se dibuja sobre un archivo de imagen llamado 'heatobjects.png' donde el color identifica el tipo de material. No hay limitación en el número de materiales pero algunos por su implementación tienen restricciones:

• 'solid', 'source', 'boundary': pueden tener cualquier geometría desde píxeles sueltos a formas complejas. Objetos inconexos pueden compartir el mismo material ahorrándonos descriptivo.
• 'fluid', 'insulate', 'isoflux': deben obligatoriamente definirse como materiales independientes contiguos (es decir cada objeto ha de ser único, con un color y descripción propios), y en el caso de 'insulate' y 'isoflux' por sencillez del código deben ser piezas rectangulares, teniendo un tamaño mínimo de 2x2 celdas. Para formas más complejas deberán combinarse varios de estos rectángulos (NOTA: cada uno debe definirse como un material independiente) o en el caso de buscar un aislante acudir a sólidos de muy baja difusividad (α).

El algoritmo iterativo no se ejecuta en las celdas perimetrales del espacio simulado, por lo que no conviene que los materiales que actualizan su temperatura con él ('solid', 'source', 'fluid') solapen con los bordes, ya que introducirían errores. Lo ideal es que los cuatro lados del perímetro consistan en materiales 'boundary', 'insulate' o 'isoflux' (p.ej. el aire que rodea un edificio o aislamientos para evitar efecto de bordes).

También es recomendable para evitar pequeños efectos de bordes, que los materiales que se adaptan a los valores de temperatura en sus vecinos ('insulate' y ' isoflux'), se especifiquen en último lugar en los ficheros descriptivos, después de todos los demás materiales. Así se garantiza que tras cada iteración los valores de temperatura en ellos queden como los requeridos para modelar la función buscada (flujo nulo o continuidad de flujo respectivamente) en el siguiente paso.

~~~

Aquí se muestra un ejemplo de 'heatobjects.png':

La primera fila de píxeles queda reservada ya que se requiere para establecer el orden de los materiales presentes en la simulación mediante píxeles de diferentes colores. Los colores que identifican los objetos presentes en la escena deben indicarse en la esquina superior izquierda encerrados por dos píxeles delimitadores que serán ignorados (pueden ser de cualquier color, pero ambos deben ser iguales):

Los parámetros térmicos de cada material se definen en el fichero de configuración 'heatobjects.csv', siendo especificados estrictamente en el mismo orden en que aparecían ordenados de izquierda a derecha los píxeles de color en la imagen anterior, ya que dicho orden es lo que asocia cada material simulado con su definición térmica. Los parámetros de cada material son:

• Temp: temperatura inicial a la que se encuentra el material (K o ºC)
• alpha: difusividad térmica (m2/s)
• k: conductividad térmica (W/(m·K))
• q: fuente de calor interna (W/m3) (se ignora si el objeto no es 'source')
• type: tipo de material a elegir entre 'boundary', 'solid', 'insulate', 'fluid' o 'source'
• desc: descriptivo alfanumérico opcional del material (p.ej. "suelo")

Las temperaturas pueden especificarse y simularse indistintamente en K o °C, ya que en la Ecuación del calor T sólo aparece en forma diferencial. Eso sí, debe elegirse la misma unidad para todos los materiales y será en la que se expresen los resultados.

Como el calor no entiende de fronteras, tener los perfiles del espacio simulado nos será muy útil para delimitar los distintos materiales superponiéndolos a los mapas de temperatura resultado de las simulaciones. Los calculamos a partir del archivo de imagen anterior:

Adicionalmente en el fichero 'heatparams.csv' se especifican los parámetros de simulación:

• cell size (m): tamaño de las celdas discretas (el mismo para las dos dimensiones)
• time step (s): salto temporal entre iteraciones
• number of iterations: número total de iteraciones a simular
• number of snapshots: número de archivos PNG de salida con la representación instantánea de la temperatura en el espacio simulado
• gamma: si distinto de 1 deslinealización de la salida para resaltar diferencias de temperatura (OUT = IN^1/gamma)

~~~

Hay que apuntar que estrictamente nuestras ecuaciones modelan la transmisión de calor por conducción, mecanismo propio de los materiales sólidos. Sin embargo en los fluidos el calor se propaga por convección, bastante más complicada de modelar.

En fluidos paralelos a una superficie de material sólido a menudo se establece una condición de contorno caracterizada por una temperatura en el infinito del fluido (T_∞), calculando de forma empírica un coeficiente de convección (h) que modele localmente el intercambio de calor entre la pared y el fluido:

Para modelar fluidos caracterizados por una temperatura en el infinito (T_∞) de acuerdo a este modelo, los especificaremos como 'boundary', forzando una temperatura constante en todo el fluido. Serían equivalentes a fluidos de coeficiente de convección (h) infinito, lo que vendría a ser una convección forzada muy rápida.

Pero el anterior esquema no sirve para aproximar la temperatura en fluidos finitos con formas arbitrarias y en contacto con múltiples otras superficies igualmente irregulares. Como no entra en el alcance de este ejercicio abordar mecánica de fluidos, pero no queremos renunciar a poder simular ese tipo de fluidos, me he sacado de la manga una aproximación de "convección instantánea" consistente en lo siguiente:

• Tratar cada fluido como un sólido más (definido por su difusividad y conductividad), calculando el intercambio de calor con todos sus vecinos de acuerdo a las ecuaciones de conducción estándar.
• Tras cada iteración se calcula la temperatura media del fluido reasignándola en todas las celdas que lo constituyen, lo cual es consistente con el Principio de conservación de la energía.

Esta aproximación implica asumir que el intercambio de calor en las fronteras del fluido tiene lugar a un ritmo mucho más lento que la redistribución del calor dentro del propio fluido por convección. Esto equivale de nuevo a una convección forzada muy rápida también en todos los fluidos especificados como 'fluid'. Como bien me ha apuntado Salvador López Alfonso, la aproximación se hace más consistente con el mundo real cuanto más reducidas sean las dimensiones de los objetos fluidos.

NOTA: si vamos a modelar más de un objeto 'fluid', aunque compartan sus parámetros termodinámicos cada uno debe definirse como un material independiente (color y parámetros), o de lo contrario al promediar su temperatura estaríamos introduciendo trasvases de calor entre zonas de la simulación no comunicadas.

~~~

Por último añadiremos al código la posibilidad de calcular la densidad del flujo de calor o gradiente del campo escalar simulado de temperaturas:

Y también puede ser interesante obtener mapas de energía acumulada, que no siempre será máxima en los cuerpos que se encuentren a mayor temperatura ya que intervienen variables adicionales (densidad y calor específico). Asumiendo una energía de referencia nula para una temperatura de 0 K:

~~~

En la segunda parte del artículo, 'Transferencia de calor por elementos finitos con R (II). Simulaciones', se usa el código de simulación para mostrar varios ejemplos como paredes con diferentes conductividades, un disipador de calor, un logo corporativo metálico con fuentes de calor,...

Repositorio con el código R y simulaciones: GitHub.