nueva página del texto (beta)
Español (pdf)
Articulo en XML
Traducción automática
Citado por SciELO 
Citado por Google 
Similares en
SciELO 
Similares en Google 
Permalink
INTRODUCCIÓN
Desde la creación de los primeros instrumentos para realizar las operaciones aritméticas más simples (v.g., el ábaco) hasta la invención de las súper computadoras, la humanidad ha buscado reducir las etapas de cálculo, el tiempo de los procesos de computación y por ende el consumo energético. Es por eso que el desarrollo de nuevos software, métodos y algoritmos cada vez más eficientes se ha convertido en un pilar fundamental en el progreso de cualquier tipo de organización. En el contexto de la matemática aplicada al estudio de la estabilidad de los reactores químicos, se han creado herramientas de cálculo versátiles para determinar condiciones seguras de operación de los mismos y prevenir incidentes en planta. Algunas de ellas básicas, como la metodología de sensibilidad paramétrica que, mediante un análisis geométrico de las variables de estado (a condiciones fijas y con variación de un parámetro), ha permitido definir criterios de estabilidad sobre la base de modelos matemáticos en estado estable. Otras se basan en respuestas operativas (por ensayo y error) que han permitido seleccionar las mejores condiciones de operación (Rangel y Bogoya, 1992; Rivera et. al, 2005) o en la aplicación de los principios matemáticos (v.g., de la primera y segunda derivada en el plano de fases (Varma et. al, 1999; Gaviria et. al, 2016) y de diferenciación directa para definir coeficientes de sensibilidad (Ojeda et. al, 2014)). Otros métodos más avanzados, como el análisis dinámico que, basado en una rigurosa y bien fundamentada teoría (v.g., singularidad y bifurcaciones (Kuznetsov, 2004)), ha permitido caracterizar en profundidad la topología de sistemas reactivos (los puntos críticos o bifurcaciones y las regiones de inestabilidad (Ojeda et. al, 2016)) a partir de modelos variantes en el tiempo (Elnashaie y Uhlig, 2007). Entre los fenómenos de inestabilidad inherentes a procesos reactivos exotérmicos, los de histéresis y/o multiplicidad de estados estacionarios son especialmente complejos y pueden surgir en sistemas de diferente naturaleza (Durán y Córdoba, 2007; Rincón et. al, 2009). De hecho, algunas de sus condiciones se consideran altamente riesgosas debido a que el más mínimo cambio en las condiciones de operación puede provocar un salto abrupto de las variables de estado y posiblemente un accidente.
Una ingeniosa metodología de cálculo fue propuesta por Hsu (1987), la cual nombró como mapeo de celdas (del inglés cell-mapping), cuyo fin es detectar los conjuntos o regiones de atracción de las soluciones del sistema (cuencas o dominios de atracción). Sobre esta metodología y sus principales características, Merillas Santos (2006) propuso, a grandes rasgos, un algoritmo para su implementación. Por otra parte, Taborda y Angulo (2015) propusieron un método de control inductivo de punto fijo para determinar y controlar dichas regiones de atracción. Más recientemente, Erazo et. al. (2017) describieron en detalle y gráficamente el método de mapeo de celdas y lo demostraron mediante su aplicación a un sistema físico. Convencionalmente, la metodología de análisis del comportamiento de determinado sistema dinámico no lineal implica su solución desde varios puntos iniciales y la definición de su “patrón de flujo” (v.g., la dirección y sentido de las trayectorias solución). Los puntos finales de las trayectorias definen el tiempo de integración necesario para alcanzarlos. Para verificar que se trata de un atractor se debe realizar el siguiente test: seleccionar una misma condición inicial e incrementar el tiempo de integración, si la solución final se conserva es un atractor. Ahora bien, para verificar a que dominio de atracción pertenece se debe realizar la siguiente prueba: tomar otra condición inicial (diferente a la anterior) e integrar con el mismo tiempo, si los puntos solución coinciden es porque corresponden al mismo dominio. Cerón et al. (2017) compararon gráficamente estas dos metodologías aplicadas a un caso matemático hipotético. De este trabajo se resalta la forma innovadora de presentar las cuencas de atracción mediante el mapeo de celdas. Hay que destacar que el algoritmo de mapeo de celdas permite analizar el espacio de estados definido de cualquier sistema en su totalidad y no por partes, reduciendo drásticamente el número de operaciones matemáticas con la configuración apropiada del algoritmo y la selección del paso de integración (Cao et. al, 2017; Erazo et. al, 2017).
En este trabajo se propone un algoritmo que permite implementar la metodología de mapeo de celdas. En él, se lleva a cabo la construcción de la estructura de celdas y la discriminación de los dominios o cuencas de atracción. Para ilustrar su funcionamiento, se seleccionaron dos casos de estudio de interés industrial: dos sistemas reactivos altamente exotérmicos, en los cuales se presenta el fenómeno de histéresis (a determinadas condiciones) y cuyos modelos se representan mediante ecuaciones diferenciales ordinarias. El primer caso corresponde a la hidrólisis de isocianato de metilo (Ojeda et al., 2016), una reacción con problemas de sensibilidad térmica (runaway) que ocasionó uno de los peores accidentes de la industria química: el denominado “desastre de Bhopal” (India) en 1984 (Kletz, 1999; Eckerman, 2005). El modelo dinámico en dos dimensiones de este sistema incluye los balances de materia y energía de la mezcla reactiva. El segundo caso corresponde a la hidrólisis de anhídrido acético (Gómez García et al., 2016), también considerada como una reacción altamente sensible a la temperatura y que ha generado múltiples incidentes, aunque en una menor escala (C&EN, 2017). Su modelo dinámico en tres dimensiones se resume a los balances de materia y energía de la mezcla reactiva y al balance de energía del fluido de servicio. El algoritmo propuesto fue implementado por los autores en el software MatLab®. Se comparan aquí los tiempos de integración requeridos tanto por la metodología convencional como por la de mapeo de celdas con el fin de evidenciar el potencial de esta última herramienta.
MARCO TEÓRICO
Esta sección está dividida en tres partes: en la primera se describe brevemente algunos conceptos básicos sobre sistemas dinámicos, en la segunda se presenta la teoría base para la construcción y manipulación de la estructura de celdas y en tercera se muestra la implementación del algoritmo por etapas.
Sistemas dinámicos clásicos
Un sistema dinámico es una función (o funciones) que describe el estado de las variables del sistema en el tiempo. En general, un sistema dinámico a tiempo continuo, con un espacio de estados de dimensión finita, puede modelarse mediante una o un conjunto de ecuaciones diferenciales como se muestra a continuación:
donde x es el vector de las variables de estado N-dimensional (v.g., concentración, temperatura, presión, etc.), t es el tiempo, ẋ es el campo vectorial (vector de las derivadas de x con respecto a t), μ es el campo paramétrico (vector de parámetros) k-dimensional (v.g., flujo de la mezcla reactiva, flujo del refrigerante, temperatura de entrada del refrigerante, etc.) y F es la función con valores vectoriales (funciones evaluadas que dependen de x, t y μ). En el diseño de reactores químicos, el sistema dinámico corresponde al modelo del reactor en estado transitorio (balances de materia y energía dependientes del tiempo).
A continuación se presentan algunos términos importantes relacionados con la temática y usados a lo largo de este trabajo: (1) Perfil: Evolución de una variable de estado con el tiempo; (2) Espacio de estados: Conjunto de posibles estados que puede alcanzar el sistema; cada estado corresponde a un único punto en el espacio de estados; (3) Estado estacionario: Punto, condición o solución en el cual las variables del sistema no cambian con tiempo; (4) Flujo: El movimiento continuo de los estados en el tiempo. Asi, un patrón de flujo es el conjunto de trayectorias que tienden a un punto o región invariante en el tiempo (atractor); (5) Dominio o cuenca de atracción: El conjunto de condiciones iniciales (o región en el espacio de estados) cuyas soluciones tienden a un atractor para un tiempo muy grande (existe uno por cada atractor); (6) Mapeo: Evolución del sistema dinámico a tiempo discreto.
Discretización del sistema, analogía del mapeo de puntos y celdas
La metodología del mapeo de puntos consiste en resolver el sistema dinámico en una serie de pasos a tiempo discreto (análogo a un mapa de Poincaré (Leine y Nijmeijer, 2004: Kuznetsov, 2004)). Esta es la principal diferencia con respecto a la metodología convencional en tiempo continuo. Dependiendo de la naturaleza y rigidez del modelo, existen diferentes alternativas para el mapeo de los puntos del flujo del sistema. Como los modelos bajo estudio corresponden a sistemas autónomos (no dependen explícitamente del tiempo) se deben tener en cuenta las siguientes etapas: (i) selección de un intervalo de tiempo (t); (ii) definición de un punto de partida dentro del plano de fases con límites definidos; (iii) solución simultánea del conjunto de ecuaciones diferenciales ordinarias (Ecuación 1). La solución después de la integración se define como la imagen de la condición inicial; y repetición del procedimiento hasta alcanzar a un atractor o una región atractora.
Para el entendimiento del algoritmo de mapeo de celdas, resulta esencial familiarizarse con la siguiente terminología: (1) Celda: Elemento definido por intervalos de los ejes del espacio de estados; (2) Imagen: Mapeo de una celda a otra. La condición inicial corresponde a las coordenadas del centro de la primera celda y se realiza la integración hasta alcanzar la segunda; (3) Celda de equilibrio: Celda cuya imagen es la misma celda. Un patrón de flujo tiende a esta celda; (4) Grupo o movimiento periódico: Conjunto de celdas que conforman una órbita cíclica. Un patrón de flujo tiende a él. A veces se lo denomina movimiento P-K, en donde K corresponde al periodo; (5) Celda periódica: Celda que hace parte de un grupo periódico; (6) Periodo: Número de celdas que conforman el grupo periódico; (7) Celda de transición: Celda o conjunto de celdas cuyos posteriores mapeos tienden a una celda de equilibrio o periódica. Se tiene en cuenta el número de pasos en alcanzarlas; (8) Celda sumidero (del inglés Sink Cell): Celda o conjunto de celdas cuya trayectoria solución se encuentra fuera de los límites del espacio de estados con celdas previamente acotado; (9) Grupo: Contador o número que se asigna a una celda de equilibrio o periódica. Al igual que a las celdas que hacen parte de su cuenca de atracción.
Construcción del espacio de estados con celdas
Exhibir las trayectorias o patrón de flujo de las soluciones del sistema desde cualquier punto de partida (plano de fases acotado) es una de las principales metas de esta metodología. Sin embargo, es evidente que físicamente existe una incertidumbre en las medidas experimentales y la precisión está limitada numéricamente por los redondeos. Conscientes de ello, se introduce un nuevo concepto: se consideran a las variables de estado como una colección de intervalos (celdas) y el plano de fases como una colección N-dimensional de celdas (espacio de estados con celdas). Por lo tanto, una celda se define como la entidad mínima e indivisible del plano de fases. La Figura 1 muestra esquemáticamente la diferencia que existe entre las metodologías del mapeo de puntos y el mapeo con celdas para un sistema bidimensional, en donde G es la función del mapeo de puntos, C es la función de mapeo de celdas, y n corresponde al número de pasos.
Fig. 1 Espacio de estados del (a) mapeo de puntos y (b) mapeo celda a celda (cell-mapping). Adaptado de Hsu (1987)
La estructura de celdas consiste de rectángulos o paralelepípedos rectangulares de tamaño uniforme (para sistemas bidimensionales o tridimensionales, respectivamente), que se pueden construir dividiendo las variables de estado xi (i=1,2,…,N) en intervalos uniformes hi, estos se reconocen por sus valores enteros sucesivos (correspondientes a las coordenadas de las celdas zi). En la Figura 2 se muestra un ejemplo de la conformación de la estructura de celdas para un sistema en dos dimensiones, en donde las coordenadas o posiciones de las celdas en el espacio de estados Zn se leen como lo presenta la Tabla 1.
Fig. 2 Construcción del espacio de estados con celdas y comportamiento del mapeo de celdas cuando se alcanza una celda de equilibrio y una celda periódica. Adaptado de Erazo et. al, (2017)
El número que aparece dentro de cada una de las celdas es simplemente un identificador (no se relaciona con la posición que ocupa en el plano de fases en Zn). Con esta información, es posible inferir que para cierto intervalo de celdas zi a lo largo del eje xi se debe cumplir con la siguiente condición:
donde zi (i=1,2,…,N) es un número entero que depende de la posición de la celda en cada eje (como se mostró en la Tabla 1) y z se define como el vector de celdas. Un punto x con componentes xi pertenece a una celda z con componentes zi si y solo si se cumple con la ecuación (2) para todo i. De esta forma se crea la estructura de celdas.
Otra opción consiste en: (i) definir el número de intervalos en cada eje, (ii) dividir los ejes del espacio de estados limitado con respecto a estos últimos y (iii) determinar el tamaño de los intervalos en cada eje. Para el desarrollo de este trabajo se seleccionó la segunda alternativa debido a su fácil implementación dentro de un algoritmo computacional. Por notación, C(z) se define como el vector de funciones de celda. El progreso que tiene el proceso de mapeo celda-a-celda puede describirse mediante la siguiente expresión:
Se dice que una celda de equilibrio (celda fija) se alcanza cuando una celda z satisface la igualdad:
Este es el caso de la celda 13, de color gris en la Figura 2.
A su vez, una celda periódica (de periodo K) o movimiento periódico se alcanza cuando una serie de K celdas diferentes z(j), j ϵ {1,…,K}, satisface las siguientes condiciones:
donde m es el número de veces que se aplica el mapeo de celdas. De forma análoga al caso anterior, las celdas 10, 16 y 17 de color negro (Figura 2) son celdas periódicas y hacen parte de un movimiento periódico (periodo K=3).
Siguiendo la teoría original del método de cell-mapping, de la Figura 2 es posible extraer la información que se muestra en la Tabla 2.
Sin embargo, con la metodología que se propone en este trabajo no es necesario guardar todos estos datos. Se demostrará que con el uso de la convención de colores (ver Tabla 4) se pueden determinar las celdas atractoras (grupo) y sus respetivas cuencas de atracción, el cual es el objetivo principal.
Método del punto central
El método del punto central se usa para aplicar el mapeo de celdas a partir de un mapeo de puntos. Asumiendo que el mapeo de puntos se realiza mediante la expresión x(n+1)=G(x(n)), el mapeo de las celdas sobre ZN se puede construir con los siguientes pasos. Inicialmente, se encuentra el punto central x(d)(n) de una determinada celda z(n) (como se muestra en la Figura 2 el mapeo por pasos siempre inicia desde el centro de la celda correspondiente). Los componentes de x(d)(n) vienen dados por:
donde zi(n) es el i-ésimo componente de z(n). Luego, se evalúa la imagen del mapeo del punto x(d)(n),
La celda en la que x(d)(n+1) se encuentra se toma como la celda z(n+1), tal como se muestra en la Figura 3. Este proceso define un mapeo de celdas C de la forma (4), cuyos componentes son:
donde x(d)(n) se relaciona con z(n) mediante la ecuación (8) e Int(y) denota el mayor entero de y.
Fig. 3 Método del punto central para crear un mapeo de celdas a partir de un mapeo de puntos. Adaptado de Hsu (1987)
Algoritmo de mapeo de celdas y convención de colores
A continuación, se describe en detalle el algoritmo propuesto (Tablas 1 y 3) para la implementación del método de cálculo del mapeo de celdas (Figuras 4 y 5). Se usará un sistema de dos ecuaciones diferenciales ordinarias y una convención de colores (Tabla 4) para un mejor entendimiento.
Tabla 4 Convención de colores de mapeo de puntos
| Valor | Color | Identificación |
|---|---|---|
| -2 | Blanco | Celda sumidero |
| -1 | Naranja | Celda que se está procesando |
| 0 | --- | Celda que no ha sido procesada |
| 1 | Gris | Grupo 1 |
| 2 | Negro | Grupo 2 |
| 3 | Azul | Grupo 3 |
| 4 | Rojo | Grupo 4 |
| … | --- | … |
| n | --- | Grupo n |
RESULTADOS
Se analizan dos casos: (1) La hidrólisis de isocianato de metilo (ICM, CH3NCO) para la formación de metilamina (CH3NH2) y dióxido de carbono (CO2); y (2) La hidrólisis de anhídrido acético ((CH3CO)2O) en presencia de ácido sulfúrico (H2SO4) como catalizador, para la formación de ácido acético (CH3COOH). Finalmente, es de importancia mencionar que todas las simulaciones se realizaron en un computador personal con una capacidad de 3.2 GHz de procesador y 8 Gb de RAM.
Caso 1
La hidrólisis de isocianato de metilo (ICM, CH3NCO) para la formación de metilamina (CH3NH2) y dióxido de carbono (CO2) fue el origen del conocido desastre de Bhopal (Union Carbide-India):
Esta reacción ocurrió por error en un tanque de almacenamiento de ICM. Los balances de materia y energía del reactor en estado transitorio se pueden expresar de la siguiente forma (Ojeda et al., 2016):
Con condiciones iniciales,
Donde m es masa de ICM (kg), c es la concentración de ICM (mol/kg), t es el tiempo (s), A es el factor de frecuencia de la reacción (s-1), Ea es la energía de activación (kJ/mol), R es la constante de los gases ideales (kJ/(mol∙K)), F es el flujo másico (kg/s), T es la temperatura (K), L es el coeficiente de transferencia de calor (W/K), CP es el calor específico de la mezcla reactiva J/(kg∙K) y ∆Hrxn es la entalpía de reacción (kJ/mol). Los subíndices 0 y a indican alimentación y ambiente (medio refrigerante), respectivamente. Las magnitudes de algunas de estas propiedades se pueden consultar en el trabajo original (Ojeda et al., 2016).
Con el fin de simplificar el modelo del sistema se introducen las siguientes variables y parámetros adimensionales (Tabla 6).
De esta manera los balances adimensionales adoptan la siguiente forma:
Las ecuaciones (14) y (15) estarán sujetas a determinadas condiciones iniciales según el desarrollo del mapeo de puntos (como se mostró en las Figuras 2 y 3). La Tabla 7 resume las condiciones de operación del sistema reactivo con multiplicidad de estados estacionarios, las cuales han sido comprobadas previamente en nuestro trabajo (Ojeda et al. 2016). Se seleccionó un espacio de estados acotado por los siguientes valores de las variables: u = [0 1] y θ = [0 0.1]. La variable u estrictamente se encuentra entre 0 y 1 debido a que es una relación de concentraciones del reactivo límite (ICM). Mientras que la variable θ tiene como valor mínimo 0 (un valor extremo, 0 K), el cual se cambia por un valor realista de 0.033 (Ta = 260 K) correspondiente a la temperatura del medio refrigerante y un valor máximo de 0.1 (T = 786.58 K), muy superior a las temperaturas de saturación de las especies involucradas (condición de runaway).
La Figura 6 (a-b) muestra los resultados de la simulación aplicando el método convencional (a) y el método de mapeo de celdas (b) para la hidrólisis de ICM. Mientras que en la Tabla 8 se presenta sus principales diferencias.
Fig. 6 Plano de fases de la hidrólisis de ICM con perfiles a diferentes condiciones iniciales (a) y mapeo de celdas (b)
Caso 2
El analiza la hidrólisis de anhídrido acético ((CH3CO)2O) en presencia de ácido sulfúrico (H2SO4) como catalizador, para la formación de ácido acético (CH3COOH) (Gómez García et al., 2016):
Esta reacción se llevó a cabo en un reactor de mezcla perfecta (CSTR de su acrónimo en inglés) equipado con un serpentín para enfriamiento. Los balances de materia y energía (de la mezcla reactiva y del fluido refrigerante) en estado transitorio adoptan la siguiente forma (Gómez García et al., 2016):
Los cuales están sujetos a las siguientes condiciones iniciales:
Donde CA es la concentración de anhídrido acético (mol/m3), Cs es la concentración de ácido sulfúrico (mol/m3), T es la temperatura de la mezcla reactiva (K), t es el tiempo (s), VR es el volumen del reactor (m3), es el flujo volumétrico (m3/s), A0 es el factor pre-exponencial (m3/(mol∙s)), E es la energía de activación (J/mol), R es la constante de gases ideales (J/(mol∙K)), ρ es la densidad de la mezcla reactiva (kg/m3), CP es el calor específico de la mezcla reactiva (J/(kg∙K)), ws es la capacitancia de la pared del reactor (J/K), ΔHrxn es la entalpía de la reacción (J/mol), and Ua es el coeficiente global de transferencia de calor (W/K). Los subíndices 0 y c indican condiciones de alimentación y refrigeración, respectivamente. Al igual que la hidrólisis de ICM, este sistema se puede simplificar incluyendo las variables y parámetros adimensionales que se muestran a continuación (Tabla 9):
Tabla 9 Definición de las variables y parámetros adimensionales. Caso hidrólisis de anhídrido acético
Por lo tanto, los balances de forma adimensional se expresan como sigue:
Las condiciones necesarias para la aparición de múltiples estados estacionarios se presentan en la Tabla 10. Toda esta información se puede verificar en el artículo de base (Gómez García et al., 2016), con sus respectivas propiedades de forma dimensional. Para este caso se analizó un espacio de estados definido por los intervalos: xA = [0 1], θ = [-4 6] y θc = [-4 6]. La conversión es una medida de la transformación que tiene el reactivo límite (anhídrido acético) y puede tomar valores entre 0 (transformación nula del reactivo) y 1 (transformación completa del reactivo). Por su parte θ y θc tienen como valor mínimo -4 (278.26 K), la cual es una temperatura inferior a la del medio refrigerante y un valor máximo de 6 (T = 365.5 K, condición de runaway).
La Figura 7 (a-b) muestra los patrones de flujo generados con el método convencional (a) y el mapeo de celdas (b). En la Tabla 11 se presenta las principales características de cada uno de estos. Al igual que la hidrólisis de ICM, el fenómeno de histéresis o de múltiples estados estacionarios se puede evitar, por ejemplo, reduciendo el número de Semenov (incrementando la remoción de calor) 
o aumentando la relación de los tiempos de residencia M > 33 (lo cual se puede obtener con la disminución del flujo volumétrico de la mezcla reactiva o el incremento del flujo del fluido refrigerante). Nuevamente se sugiere llevar a cabo un análisis dinámico para obtener más información sobre comportamiento del sistema con respecto a la fluctuación de parámetros (Gómez García et al., 2016).
Fig. 7 Plano de fases de la hidrólisis de anhídrido acético a diferentes condiciones iniciales (a) y mapeo de celdas (b)
CONCLUSIONES
En este trabajo se presentó un procedimiento detallado para la implementación del algoritmo de mapeo de celdas y la construcción del espacio de estados con celdas. Se seleccionaron dos sistemas reactivos químicos que presentan fenómenos de histéresis y/o multiplicidad de estados estacionarios: (i) la hidrólisis de isocianato de metilo, como un sistema de dos dimensiones (balances de materia y energía de la mezcla reactiva), y (ii) la hidrólisis de anhídrido acético, como un sistema de tres dimensiones (balances de materia y energía de la mezcla reactiva y balance de energía del fluido de servicio). En ambos casos se identificaron dos estados estacionarios estables y una separatríz originada por la variedad inestable de un punto silla (estado estacionario inestable), la cual divide el espacio de estados. Mediante la estructura de celdas y la convención de colores de la misma se logró reconocer los dominios de atracción de los sistemas reactivos. Se comparó el tiempo computacional o la velocidad de cálculo de las diferentes operaciones con la metodología convencional. Se demostró que el algoritmo de mapeo de celdas es una herramienta potente de cálculo significativamente superior, entre cuatro y nueve veces más rápida en comparación con la metodología convencional.
