Los modelos dinámicos de sistemas con parámetros distribuidos en una dimensión dan origen a ecuaciones diferenciales parciales que tienen como variables independientes al tiempo y a la dimensión espacial considerada. Generalmente, la complejidad de estas ecuaciones hace imposible encontrar una solución analítica; para esos casos se desarrollaron varios métodos numéricos. Estos métodos pueden ser clasificados en dos grandes grupos de acuerdo a la forma en que describen el movimiento del fluido involucrado. El primer grupo está compuesto por los algoritmos eulerianos (Bank, 1994; Lang, 1998): los nodos se mantienen en coordenadas fijas mientras que el fluido los atraviesa. El segundo grupo está compuesto por los algoritmos lagrangianos (Baines, 1994), en los cuales los nodos se mueven con el fluido. Una clasificación complementaria es desde el punto de vista de la existencia o no de una malla de cálculo. Entre los métodos clásicos que emplean malla se encuentran los métodos de diferencias finitas, elementos finitos y volúmenes finitos, lo cuales tienen versiones para el enfoque euleriano y el lagrangiano (Neunzert y Struckmeier, 1997). En años recientes se desarrollaron métodos sin mallas o con mallas adaptativa debido a su capacidad para manejar cambios geométricos en el dominio de interés. Una razón adicional para desarrollar esos métodos es que la generación de la malla es la parte que más tiempo consume en los métodos que la utilizan, típicamente la generación de la malla consumen más del 70 % del tiempo total de cómputo.
En este trabajo, se utiliza una variante del método de volúmenes finitos para resolver sistemas de ecuaciones diferenciales que componen el modelo dinámico de un intercambiador de calor de doble tubo. Esta técnica utiliza el enfoque lagrangiano –los volúmenes se mueven a la misma velocidad que el fluido–. Para ello, el dominio se divide en volúmenes finitos. En algunos casos, no es necesaria la existencia de una malla para distribuir esos volúmenes finitos. Luego, se crea el modelo correspondiente a cada volumen finito considerando que se mueve junto con el fluido. El método propuesto se compara a través de dos casos de estudio con el método de diferencias finitas. De la comparación surge que el método propuesto es más estable y más fácil de implementar. Si bien el presente trabajo se restringe al estudio de sistemas con parámetros distribuidos en una única dirección, el método propuesto tiene un amplio campo de aplicación en la simulación de equipos de la industria química (Narataruksa et al., 2008; Luciu et al., 2009).
El sistema a estudiar es el mostrado en la Figura 1. Se trata de un intercambiador de calor de doble tubo donde se calienta agua fría, que entra en el tubo interno a 25 °C, utilizando agua caliente, que ingresa en el tubo externo a 60 °C, en contracorriente.
El modelo dinámico da lugar a un sistema de ecuaciones diferenciales parciales, donde las variables independientes son el tiempo t y la coordenada espacial z. La coordenada radial no se considera porque se supone flujo pistón. Se supone también que ambas áreas transversales –la interna y la externa– son iguales. Con esas consideraciones, el modelo correspondiente es:
donde P es la densidad del fluido, Cp es la capacidad calorífica del fluido, T es la temperatura del fluido en el tubo interno, Ts es la temperatura del fluido en el tubo externo, U es el coeficiente global de transferencia de energía, d es el diámetro del tubo interno, v es la velocidad del fluido en el tubo interno, vs es la velocidad del fluido en el tubo externo (se la toma como positiva), t es el tiempo y z es la coordenada longitudinal.
La forma tradicional de resolver el sistema de ecuaciones presentado en la sección anterior es utilizando el método de diferencias finitas. En ese método, el intercambiador se descompone en varias secciones de longitud , cilíndricas para el tubo interno, anulares para el tubo externo. Esas secciones están fijas como lo muestra la Figura 2. A su vez, el tiempo se discretiza utilizando un intervalo . Como resultado de la discretización, se tienen las matrices de temperatura T y Ts. La temperatura Tk,i es la temperatura de la sección diferencial i, con centro en (i+0.5) , del tubo interno en el tiempo k . La temperatura Tsk,i es la temperatura de la sección diferencial i, con centro en (i+0.5) , del tubo externo en el tiempo k . Tanto el índice k como el índice i comienzan en 0.
Como resultado de ambas discretizaciones (la espacial y la temporal), el sistema original de ecuaciones diferenciales parciales se transforma en un sistema de ecuaciones algebraicas.
Dicho sistema está compuesto por las siguientes ecuaciones
La resolución del sistema algebraico así obtenido requiere que se especifiquen los perfiles iniciales en el interior de ambos tubos y las condiciones de frontera. Como parte de la condición inicial se supone que ambos tubos están llenos. Si bien puede relajarse esta última condición, ello requiere de una programación adicional que no es contemplada por el método original. Con los perfiles iniciales como datos, se despejan de las ecuaciones anteriores Tk+1,i y Tsk+1,i para todo i, tomando k = 0. De este modo se estiman los perfiles correspondientes a k = 1, los cuales se toman como datos para calcular los perfiles para k = 2, y así sucesivamente.
El programa para obtener la solución fue implementado en la aplicación Mathcad®. La Tabla 1 presenta los valores adoptados para los parámetros del modelo. La Figura 3 presenta la evolución del perfil de temperatura en el tubo interno, mientras que la Figura 4 presenta la evolución del perfil de temperatura en el tubo externo. En ambas figuras se exhibe la evolución para los cuatro primeros minutos –cada perfil está espaciado por un minuto–. Como puede apreciarse en dichas figuras, el estado inicial que se asumió es considerar ambos tubos llenos con agua a 25 °C. Como condiciones de frontera se toman las temperaturas de entrada de la corriente fría y de la caliente. La Figura 5 presenta el perfil que se logra a los 400 s, ese perfil es el correspondiente al estado estacionario que finalmente alcanza el equipo. La diferencia de temperatura entre ambas corrientes en el estado estacionario estimado varía punto a punto, estando los valores en el intervalo [6.30 °C, 6.73 °C]. Por otra parte, resolviendo el modelo estacionario, que surge del modelo dinámico anulando las derivadas temporales, se obtiene que la diferencia de temperatura en ese estado es constante en cada punto, y es igual a 6.6 °C. La diferencia entre las dos simulaciones se debe a que la simulación dinámica involucra más cálculos, amplificando el efecto de los errores en cada uno de ellos. La magnitud del impacto de los errores de los sucesivos cálculos depende de los valores de y . El valor de adoptado en la Tabla 1 está cerca del valor crítico, un valor mayor causa inestabilidad en el sistema –aparecen oscilaciones de gran amplitud–; y no puede obtenerse una solución apropiada.
En el método que se propone en este trabajo también se consideran volúmenes de longitud z –cilíndricos para el tubo interno, anulares para el tubo externo–; pero en este caso, dichos volúmenes se mueven a lo largo de los tubos a la misma velocidad que el fluido. En otras palabras, mientras que el método de las diferencias finitas utiliza la representación euleriana, el método propuesto utiliza la representación lagrangiana. La primera consecuencia de este enfoque es que el modelo correspondiente a cada volumen finito móvil es más simple que el correspondiente a cada sección del método analizado anteriormente. En efecto, el modelo para el volumen finito móvil j es el siguiente:
El modelo es ahora más simple porque esta vez sólo se considera el calor que intercambian las dos corrientes (el segundo término de las ecuaciones que modelan las secciones fijas, ecuaciones (1) y (2)). Al moverse los volúmenes finitos a la misma velocidad que el fluido, nada ingresa ni sale de ellos; por lo tanto, desaparece el término convectivo (el primer término de las ecuaciones (1) y (2)). Esta forma de plantear el modelo para cada volumen finito es la principal diferencia con otros métodos que utilizan la representación lagrangiana. En el método propuesto, el modelo no está constituido por funciones de aproximación como es el caso en otros métodos, sino que es el modelo que surge de tomar el volumen finito móvil como volumen de control para plantear los balances requeridos.
Una forma alternativa de obtener este modelo es procediendo de la misma manera que en el método de volúmenes finitos fijos (integrando las ecuaciones (1) y (2) en el volumen finito) pero considerando que esta vez los volúmenes se mueven con el fluido.
Una importante consecuencia de considerar que los volúmenes finitos se mueven con el fluido es que el modelo tiene ahora una única variable independiente: el tiempo. Esto implica una reducción importante en el número de cálculos necesarios y, en consecuencia, un menor efecto de los errores numéricos. De lo analizado, se debe esperar que el método propuesto sea más estable que el método de las diferencias finitas estudiado en la sección anterior. Por los mismos motivos, el método propuesto es más fácil de implementar que el método presentado en la sección anterior.
El modelo discretizado aplicando diferencias finitas en la dimensión temporal es:
Nuevamente, la resolución del sistema algebraico así obtenido requiere que se especifiquen los perfiles iniciales en el interior de ambos tubos. La simulación se inicia suponiendo que ambos tubos están llenos y a una temperatura de 25 °C, como en el caso anterior; sin embargo, esta vez el método permite, sin mayor dificultad, que se pueda tomar como estado inicial al intercambiador con los dos tubos vacíos. Como condiciones de frontera se mantienen las temperaturas de entrada de la corriente fría y de la caliente. Con los perfiles iniciales como datos, se despejan de las ecuaciones anteriores Tk+1,j y Tsk+1,j para todo volumen j, tomando k = 0. De ese modo se estiman los perfiles correspondientes a k = 1, los cuales se toman como datos para calcular los perfiles para k = 2, y así sucesivamente.
Si bien el modelo se puede implementar en cualquier aplicación matemática o con cualquier lenguaje de programación, se eligió trabajar con el entorno de desarrollo NetLogo (Wilensky, 2010). Dicho entorno está orientado a la simulación de fenómenos naturales y sociales. NetLogo fue creado en 1999 por Uri Wilensky, y desde entonces es desarrollado por el Centro para el Aprendizaje Conectado y el Modelado Basado en Computadoras (Center for Connected Learning and Computer-Based Modeling). Este entorno es particularmente apropiado para modelar sistemas complejos que evolucionan en el tiempo. El modelado se realiza mediante instrucciones a cientos o miles de agentes que operan en forma independiente. Este enfoque hace posible analizar la conexión que existe entre la conducta de los agentes individuales y la conducta de la comunidad de agentes que emerge de las muchas interacciones individuales. Para mantener la simplicidad de la programación, se fijó el paso de discretización temporal de la siguiente manera:
donde p es un entero positivo. La incorporación de esta restricción asegura que en cada t siempre un volumen finito del tubo interno está en contacto con un único volumen finito del tubo externo (estos es, ambos volúmenes quedan alineados). Cuando p = 1, cada volumen finito del tubo interno hace contacto con todos los volúmenes finitos del tubo externo, es el valor de p que produce el menor error. Cuando p = 2, cada volumen finito con índice j par del tubo interno hace contacto sólo con los volúmenes finitos del tubo externo con índice j par; a su vez, cada volumen finito con índice j impar del tubo interno hace contacto sólo con los volúmenes finitos del tubo externo con índice j impares; por lo tanto, en este caso se tiene un error mayor que el obtenido para p = 1.
Los resultados obtenidos para el ejemplo que se está analizando fueron similares a los obtenidos por el método de diferencias finitas. En cuanto a la estabilidad, no puede realizarse un estudio completo debido a que la restricción planteada por la ecuación (9) no permite variar y en forma independiente. Por ese motivo, para analizar la estabilidad se presenta otro ejemplo en la sección siguiente. Lo que sí puede afirmarse para este ejemplo es que el método propuesto es al menos tan estable como el de las diferencias finitas porque produjo resultados similares para valores de = 0.4 m y = 0.8 m/s (comparables con los valores de reportados en la Tabla 1). Esos resultados se obtuvieron tomando un valor de p = 2, que no es el caso más favorable para el método propuesto
La estabilidad de un método iterativo depende de la forma en que evoluciona el error a lo largo de las sucesivas iteraciones que involucra el método. Si el error decae, el método es estable; si no, es inestable. Para problemas dependientes del tiempo, la estabilidad garantiza que el método produzca una solución acotada siempre que la solución exacta de la ecuación diferencial sea también acotada. En esta sección, se analiza el caso en que puede determinarse en forma analítica la estabilidad del método de diferencias finitas, y se muestra que el método propuesto es más estable.
El análisis de estabilidad de von Neumann (también conocido como análisis de estabilidad de Fourier) permite determinar la estabilidad del método de diferencias finitas cuando se aplica a la resolución de ecuaciones diferenciales parciales lineales. Para que este estudio sea necesario y suficiente para garantizar la estabilidad, se deben cumplir una serie de condiciones. Un problema que cumple con todas las condiciones requeridas por el estudio es el modelo de una onda unidireccional (one way wave):
resuelto por el siguiente esquema de diferencias finitas:
con c > 0. Para el citado problema, el análisis de von Neumann establece que el método será estable si se cumple:
Para el caso del intercambiador de calor en estudio, un modelo equivalente al de una onda unidireccional se puede conseguir si se considera que en el tubo externo condensa vapor.
En ese caso, el modelo para el tubo interno es:
resuelto por el siguiente esquema de diferencias finitas:
con v > 0. Para el citado problema, el análisis de von Neumann establece que el método será estable si se cumple:
La Tabla 2 presenta los valores adoptados para los parámetros del modelo. Se supone como perfil inicial para el tubo interno que el agua está en su interior a 25 °C; y además se supone que el tubo externo está lleno con vapor suturado a 100 ºC. La Figura 6 muestra la evolución para 20 s que en este caso tiene la temperatura en el tubo interno, los perfiles se tomaron cada 5 s. La solución fue obtenida con el máximo t que permite que el método sea estable de acuerdo con el análisis de von Neumann.
En el método que se propone en este trabajo, el modelo para cada volumen finito móvil en el tubo interior es el mismo que antes:
El modelo discretizado aplicando diferencias finitas en la dimensión temporal es:
Al no existir en este caso los correspondientes volúmenes finitos para el tubo exterior, no fue necesario plantar la restricción de la ecuación (9); por lo tanto, y fueron fijados en forma independiente.
Nuevamente la solución alcanzada fue similar a la obtenida por medio del método de diferencias finitas. Sin embargo esta vez, el método propuesto es más estable ya que se pudo incrementar hasta diez veces el valor original, produciendo un error en la temperatura de salida del tubo interno del 2 % con respecto a la solución inicial.
El método planteado tiene ventajas adicionales. Si se modela un único volumen finito móvil, y se grafica la evolución que tiene su temperatura en función de la posición que ocupa desde que el volumen finito ingresa hasta que sale del tubo interno, se obtiene el perfil correspondiente al estado estacionario que se muestra en la Figura 6. Esto se debe a que todos los volúmenes finitos que le siguen se ven sometidos a las mismas condiciones.
Por otra parte, dada la naturaleza del modelo empleado, también es de esperar que el método propuesto reproduzca mejor la evolución de las perturbaciones que tengan lugar en las fronteras, sobre todo las perturbaciones de tipo escalón; lo cual fue confirmado por experiencias que no se incluyen en el presente.
En el presente trabajo se presentó un método para resolver sistemas de ecuaciones diferenciales parciales que modelan un intercambiador de calor de doble tubo. El método fue comparado con el método de diferencias finitas, mostrando claras ventajas desde el punto de vista de la estabilidad de la solución, como así también desde el punto de vista de la simplicidad del modelo y su implementación.
Si bien es ampliamente aceptado que los métodos con malla móvil son más complejos de implementar que los de malla fija, en los problemas estudiados en este trabajo se observó lo contrario. Ello puede deberse a que la mayoría de los trabajos realizados en el área se refieren a problemas en dos o tres dimensiones, mientras que el trabajo presentado se restringió a problemas de una única dimensión. Además, se consideró que la densidad del fluido era constante, lo que permitió trabajar con volúmenes finitos de tamaño constante. A pesar de las condiciones asumidas, son muchos los modelos de sistemas de interés en la industria química que cumplen con esa condición.
Baines M.J., (1994), Moving finite elements, Monographs on Numerical Analysis, The Clarendon Press Oxford University Press, New York.
Bank R.E., (1994) PLTMG: a software package for solving elliptic partial differential equations, vol. 15 of Frontiers in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA.
Lang J., (1998) “Adaptive FEM for reaction-diffusion equations”, Appl. Numer. Math., 26, 105–116.
Luciu R., T. Mateescu, V. Cotorobai y T. Mare, (2009) “Number and Convection Heat Transfer Coefficient for a Coaxial Heat Exchanger Using Al2O3-Water pH=5 Nanofluid”, Bulletin of
the Polytechnic Institute of Jassy, Constructions, Architecture Section, 55 (59), 2, 71-80.
Narataruksa P., P. Triratana, K. Pana-Suppamassadu, P.J. Heggs y S. Tia, (2008) “Dynamic simulation of plate and frame heat exchanger undergoing food fouling: coconut milk fouling case study”, Science Asia, 34, 229–237.
Neunzert H. y J. Struckmeier, (1997) “Boltzmann Simulation by Particle Methods, in: Boltzmann's Legacy 150 years after his birth”, Accademia Nazionale dei Lincei, 131, 59-78. Wilensky U., NetLogo, (2010) http://ccl.northwestern.edu/netlogo/index.shtml.