Mathematical Model of a Methanol Synthesis Quench Reactor with Catalyst Deactivation

Molina-Cárdenas Gabriel Alberto¹, Rosales-Chinchilla Andrés Emilio², Papa-Annibalini José³

¹ Escuela de Ingeniería Química Facultad de Ingeniería Universidad Central de Venezuela, Caracas Correo: gabrielmolina87@gmail.com

² Escuela de Ingeniería Química Facultad de Ingeniería Universidad Central de Venezuela, Caracas Correo: rosales.andres@gmail.com

Información del artículo: recibido: octubre de 2012,
Reevaluado: mayo de 2013 y agosto de 2013,
Aceptado: septiembre de 2013

Resumen

En el presente trabajo se desarrolló un modelo heterogéneo basado en ecuaciones diferenciales con el fin de simular un reactor de síntesis de metanol con interenfriamiento en condiciones pseudo-estacionarias. Para ello se utilizó información abierta de la literatura con ajustes apropiados para el caso particular analizado. Se introdujo un factor de actividad que tiene en cuenta la desactivación del catalizador en función del tiempo y temperatura de operación y se usó una expresión para las velocidades de reacción, ya disponible, con parámetros ajustados para el catalizador en uso en el reactor modelado. El modelo es capaz de predecir satisfactoriamente el comportamiento de los perfiles de temperatura, así como el flujo de efluente y su composición a lo largo de la vida útil del catalizador, obteniendo una herramienta apropiada para el análisis y optimización del proceso, así como para programar los cambios operacionales necesarios para mantener el rendimiento. Esto se logra con una reducción significativa del tiempo de cálculo considerando todos los fenómenos importantes involucrados.

Descriptores: síntesis de metanol, reactor de lecho catalítico, simulación en estado estacionario, modelo matemático, desactivación de catalizador.

In this paper a heterogeneous model based on differential equations to simulate a methanol synthesis quench reactor in pseudo-steady state conditions was developed. To achieve the objective, open information available in the literature for our particular case was used, with the appropriate adjustments. In order to take into account the activity decay of the catalyst, an activity factor was introduced as a function of the time on run and the operation temperature. An available expression for reaction rates with parameters adjusted to our particular case was used. The model is able to successfully predict the behaviour of the temperature profiles, the effluent flow and composition for the whole catalyst life, yielding a useful tool for analysis, process optimization and scheduled periodic operational changes needed to hold the reactor performance. The objectives were achieved also with a significant reduction in the computing time, without losing the influence of any significant phenomena.

Keywords: methanol synthesis, catalytic bed reactor, steady-state simulation, mathematical model, catalyst deactivation.

Introducción

La simulación de reactores de lecho catalítico es de suma importancia tanto para el diseño como para la optimización de procesos químicos industriales en los que se ven involucrados y es una herramienta de gran valor para la operación segura de las plantas, mantener su capacidad productiva y programar las acciones que sean necesarias. En este trabajo se estudia el proceso de síntesis de metanol, donde cada planta industrial en operación cuenta con un reactor diseñado y operado en condiciones específicas que pueden variar dependiendo del catalizador utilizado, la energía disponible, la materia prima, etcétera. Debido a estas particularidades existe la necesidad de crear un modelo adecuado a cada caso específico, el cual considere todos los aspectos relevantes del proceso.

En este trabajo se presenta el desarrollo de un modelo para un reactor de síntesis de metanol instalado y funcionando, el cual utiliza un catalizador del tipo Cu-Zn-Al. Los modelos usados para este tipo de reactores son complejos, incorporando una serie de ecuaciones diferenciales no lineales y la predicción de propiedades físico-químicas necesarias durante el proceso de resolución, lo cual requiere simplificaciones razonables para reducir el tiempo de cálculo a valores convenientes sin introducir errores significativos. A estas dificultades hay que agregar las propias de un reactor constituido por lechos múltiples con enfriamiento entre etapas por agregado de reactantes fríos.

Graaf et al. (1990) modelaron la síntesis de metanol usando un catalizador comercial de Cu-Zn-Al y demostraron que, para los tamaños de pastillas usados, la reacción está limitada por difusión interna dentro de la matriz del catalizador. Lommerts et al. (2000) viendo la complejidad y el tiempo de simulación de un modelo que considere la difusión interna, comparó varias alternativas de cálculo y determinó que la simplificación de las expresiones cinéticas con la introducción de un factor de efectividad en función del módulo de Thiele representaba correctamente el fenómeno. Lovik (2001) desarrolló un modelo con el que realizó estimaciones y optimización del proceso de síntesis de metanol considerando la desactivación del catalizador, generando así un modelo pseudo-estacionario riguroso. Velardi y Barresi (2002) propusieron un modelo para la síntesis de metanol en un proceso con tres lechos catalíticos analizando el efecto del cambio en los puntos de alimentación a determinados tiempos. Shahrokhi y Baghmisheh (2005) realizaron la simulación de un reactor de lecho fijo para síntesis de metanol, tanto en estado estacionario como dinámico, usando distintas formas de cálculo para el modelo de difusión interna, determinando que todas conducían a una simulación adecuada del comportamiento del reactor. De los modelos analizados y de los resultados de diversos investigadores se demuestra que los fenómenos más importantes a tener en cuenta para el desarrollo de un modelo adecuado son: la difusión interna dentro de la matriz del catalizador, la actividad del mismo y el uso de una cinética apropiada.

El presente trabajo se enfoca en modelar el comportamiento de un reactor de síntesis de metanol de lecho empacado, adiabático, con interenfriamiento y operado en estado pseudo-estacionario, teniendo en cuenta  todos los fenómenos antes descritos,  que pueda utilizarse en el análisis, optimización y seguimiento de un caso real en estado de producción. En el modelo se incorporaran las simplificaciones apropiadas sugeridas por la literatura con el propósito de reducir el tiempo de cálculo a un valor razonable sin perder la capacidad de reproducir el comportamiento físico del reactor. Se describe el funcionamiento del reactor, sus complejidades internas y se muestra la estrategia usada para resolver las ecuaciones implicadas. Finalmente se analizan los resultados y su importancia para programar las actividades de operación, mantenimiento y cambio del catalizador.

Como se puede ver en la figura 1, el reactor consta de 5 lechos catalíticos con interenfriamiento. El flujo de alimentación es dividido en 6 corrientes, la primera de ellas es precalentada con el efluente del reactor  y luego alimentada por el tope del mismo, su temperatura se ajusta aproximadamente a 503 K por mezclado con la cantidad necesaria de alimentación fresca, cuya composición típica es la especificada en la tabla 1.

Al entrar esta corriente a la primera etapa del reactor, reacciona produciendo calor con el consiguiente aumento de temperatura, luego el gas que sale del primer lecho es enfriado mezclando la cantidad adecuada de alimentación fresca a menor temperatura antes de entrar al segundo lecho catalítico, esquema que se repite hasta el último lecho. El catalizador es de óxido de cobre y zinc soportados sobre alúmina con forma cilíndrica de 6x5 mm. La presión típica de operación ronda los 10 MPa y las temperaturas en los lechos están entre 503 y 533 K. Aunque estas condiciones son típicas en el esquema de operación se hacen ajustes progresivos a la presión del reactor y flujos de interenfriamiento para mantener la productividad del reactor y mitigar perturbaciones.

Modelo del reactor

El modelo consta de una serie de sub modelos que, acoplados entre sí, logran la simulación del equipo. Para hacerlo se hicieron las siguientes suposiciones:

1) se consideró un modelo heterogéneo con variaciones solo en la dirección axial del reactor dada la elevada relación entre el diámetro del reactor y el de las pastillas de catalizador,

2) la caída de presión es despreciable de modo que esta puede considerarse uniforme en todo el reactor,

3) los efectos de difusión interna y desactivación del catalizador se consideran mediante factores de eficiencia y

4) los lechos se consideran adiabáticos.

A continuación se describen los sub modelos implicados:

Modelo de un lecho catalítico

El modelo consta de un balance de masa integro-diferencial por cada una de las especies reactantes y un balance integro-diferencial de energía, los cuales fueron deducidos usando las leyes de conservación correspondientes (Froment y Bischoff, 1979; Missen et al., 1999; Rosner, 2000). Esto generó la siguiente serie de ecuaciones diferenciales.

Balances de masa:

donde:

F_i = flujo molar de la especie "i" en el lecho (mol/s)

z = longitud medida desde la entrada del gas al lecho (m)

ρ_cat = densidad del catalizador (kg/m³)

Θ = fracción vacía del lecho una vez que contiene el catalizador (m³/m³)

D = diámetro del lecho (m)

r_met = velocidad de formación metanol según la reacción 8 (molmetanol/(kgcat.s))

r_wgs = velocidad de formación dióxido de carbono según la reacción 9 (molCO₂/(kgcat.s))

η_wgs = factor de eficiencia para la reacción 9 (adimensional)

a = factor de actividad del catalizador (adimensional)

T = temperatura en el lecho (K)

Cp_i = capacidad calorífica de la especie "i" (Joule/ (mol.K))

ΔH_met = calor de reacción asociado con la reacción 8 (Joule/mol)

ΔH_wgs = calor de reacción asociado con la reacción 9 (Joule/mol)

Para cada ecuación la condición inicial corresponde al valor de la variable a la entrada del lecho catalítico que corresponda. Como se aprecia, se requiere el cálculo de los factores de eficiencia "η_met," y "η_wgs" para la difusión interna, las velocidades de reacción "r_met,", "r_wgs" y el factor de actividad "a" para considerar la desactivación por sinterización con el tiempo en operación. La escala de tiempos involucrada en este último fenómeno es suficientemente grande, comparada con el tiempo de contacto, como para permitir un enfoque en estado pseudo-estacionario, esto se logra incorporando el factor de actividad del catalizador indicado por el modelo de sinterización (sección del modelo de desactivación).

Modelo cinético

donde:
P_i = presión parcial de la especie "i" en el lecho (bar).

En la tabla 2 se especifican los valores base adoptados para los parámetros cinéticos extraídos del trabajo de Vanden y Froment (1996). Si bien la cinética se describe por estas ecuaciones, para lograr resultados satisfactorios se incorpora un factor de corrección para las energías de activación "ψ", el cual no está presente en el trabajo de Vanden y Froment (1996). Este parámetro disminuye levemente la energía de activación necesaria para la reacción y se incorpora en cada uno de los lechos para mejorar las predicciones del modelo. La razón principal de esta incorporación es que el catalizador del reactor analizado no es del mismo fabricante que el de la investigación de Vanden y Froment. En la tabla 3 se muestran los factores utilizados que se determinaron con el método de Nelder-Mead (Mathews y Fink, 2004) minimizando la diferencia entre el perfil de temperatura predicho por el modelo y reflejado en planta.

Este modelo fue tomado de Lommerts et al. (2000). En su trabajo se plantea una simplificación de las ecuaciones cinéticas para las reacciones (8) y (9), posteriormente se resuelven analíticamente las ecuaciones de balance diferencial de masa y energía correspondientes a la pastilla del catalizador considerando una geometría esférica, en donde se obtienen los factores de eficiencia mediante una serie de ecuaciones algebraicas. Para lograr el diámetro esférico equivalente para un pellet cilíndrico se debe usar la relación algebraica apropiada (Rosner, 2000). A continuación se muestra el modelo.

donde:

C_i = concentración molar del compuesto "i" (mol/m³)

k_pseuMet = constante pseudo-cinética para la reacción 8

k_pseuWgs = constante pseudo-cinética para la reacción 9

φ_Met = módulo de Thiele modificado para la pseudo-reacción 8 (adimensional)

φ_wgs= módulo de Thiele modificado para la pseudo-reacción 9 (adimensional)

D_{efec_j} = difusividad efectiva del componente clave "j" en la matriz de catalizador (m²/s)

Modelo de desactivación

Este modelo fue tomado de Lovik (2001), en donde se aprecia que la desactivación depende de dos factores, la temperatura y el tiempo de operación del reactor, lo que quiere decir que se toma en cuenta el efecto de sinterización térmica que sufre el catalizador durante su vida útil. Lovik usó el modelo como una ecuación diferencial ordinaria, la cual es función de la temperatura y el tiempo de operación.
Dado que en la operación del reactor los perfiles de temperatura en el mismo están controlados, se supone que no cambian de forma importante a lo largo del tiempo, por tanto, la ecuación diferencial de Lovik se integró analíticamente y solo en función del tiempo de operación, generando las ecuaciones algebraicas 18 y 19, de donde se obtiene el factor de actividad "a", lo cual hace al modelo pseudo-estacionario.

donde:

A = factor de actividad relativa (adimensional)

K_d = factor de frecuencia de desactivación (4.39 E-3 días^-1)

E_d = energía de activación de sinterización (91270 Joule/mol)

a_o = factor de actividad de equilibrio (0.4 adimensional)

T_o = temperatura de referencia para el cálculo del factor de actividad (513 K)

t = tiempo de operación del reactor (días)

Modelo del punto de mezcla (quench)

Para el punto de mezcla se supuso que no hay efecto fluido dinámico apreciable, por tanto solo se realizó un balance de masa y energía que considera la corriente de alimentación fresca y el efluente del lecho anterior. Estos se muestran a continuación:

Balance de masa:

Balance de energía:

donde los sub índices:

q = referente al flujo de interenfriamiento

a = referente al flujo del lecho anterior

s = referente al flujo al lecho siguiente

Adicionalmente se requirió un método para el cálculo de propiedades tanto termodinámicas como de transporte. En la tabla 4 se listan los métodos utilizados para dichos cálculos.

Solución realizada para la simulación del reactor

Se supuso estado estacionario y en forma superpuesta el fenómeno de sinterización, lo cual es posible porque las escalas de tiempos para el segundo fenómeno son varios órdenes de magnitud superior. Se consideró que bajo las condiciones de operación a un tiempo en uso dado del catalizador, el reactor puede simularse suponiendo que ha alcanzado el estado estacionario, pero utilizando el factor de actividad calculado para el tiempo en operación del catalizador.

Resultados obtenidos de las simulaciones

Los resultados de la simulación se compararon con un reactor industrial. Se cotejó el perfil de temperatura, el flujo de efluente, la composición de productos y conversión de reactivos a lo largo de la vida útil del catalizador, la cual ronda los 5 años. En la tabla 5 se muestran las comparaciones de estos valores para un día dentro de los años de operación considerados y se puede apreciar la buena semejanza entre los resultados de la simulación y los datos de planta, con errores relativos menores que 3% para el flujo de efluente, menores que 8% en las conversiones de reactantes y menores que 4% en las concentraciones de metanol. En las figuras 3, 4 y 5 se pueden observar los perfiles de temperatura del reactor para tres días de operación a lo largo de la vida útil del reactor, nuevamente se ve la buena semejanza entre los datos de planta y la simulación realizada.

Finalmente se realizó un análisis con el modelo, que  consistió en dejar fijas todas las variables de operación inciales en el proceso y progresivamente ir aumentando el tiempo de operación en el modelo. Después, se compararon los resultados del modelo con el esquema de operación utilizado en planta, donde se modifican distintas variables, en este caso la presión del reactor aumentó 2.75 kg/cm2 (269.7 kPa) para el segundo año, 0.47 (46.1 kPa), 0.33 (32.4 kPa) y 1.16 (113.8 kPa) para el tercero, cuarto y quinto año, respectivamente; los flujos de interenfriamiento no presentan un patrón específico a lo largo de la operación, ya que estos cambian debido a la aparición de perturbaciones con el fin de mitigarlas, la composición del gas de alimentación es la misma y la temperatura de alimentación del interenfriamiento aumenta aproximadamente 2.5 K cada año. Como puede observarse en la figura 6, si no se recurriera al ajuste de condiciones operativas de tiempo en tiempo, el rendimiento del reactor disminuiría constantemente debido fundamentalmente al fenómeno de sinterización. Los datos reales del reactor están superpuestos sobre estas gráficas unidos entre sí por líneas de tendencia, las cuales  muestran el efecto de la manipulación de las  variables operativas. Esta manipulación se hace con la finalidad de mantener el nivel de rendimiento de la planta.

En la figura 6 se muestra además, que a medida que pasa el tiempo, el modelo predice un aumento en el flujo de salida, lo cual puede atribuirse al cambio de actividad que redunda en una reducción en la conversión por paso, lo cual es congruente con la disminución sostenida de rendimiento que se muestra en la figura 6 (parte derecha). Por lo tanto, si no se hacen modificaciones en las variables operacionales la cantidad de gas de síntesis que se transforma en metanol disminuiría en forma sostenida. Para evitarlo se recurre fundamentalmente a un aumento de la presión. Las velocidades de reacción son sensibles a la presión y, como la reacción ocurre con una disminución de moles, la posición de equilibrio se desplaza hacia los productos, lo cual también contribuye a aumentar la velocidad neta de producción. La variable que en toda manipulación se trata de mantener a niveles estables es la temperatura de operación, ya que si se permitiera la variación hacia valores más elevados se incentivaría en forma importante el fenómeno de sinterización, lo que disminuiría en forma sensible el tiempo máximo en operación para el catalizador.

Por otra parte, la manipulación de las variables operacionales hace que los flujos totales de salida y la producción de metanol presenten saltos importantes, tal como se muestra en la figura 6. En consecuencia, el modelo desarrollado permite predecir los momentos aconsejables para introducir cambios operacionales. Como se mencionó, el parámetro de control que más afecta el rendimiento es la presión, lo cual coincide con lo que enseña la termodinámica. Es mucho menos aconsejable modificar el perfil de temperatura, ya que su aumento acerca la restricción de equilibrio y aumenta la velocidad de sinterización. Si analizamos la evolución anual de los rendimientos con el modelo desarrollado, se podría en principio predecir el tiempo adicional en operación que le queda al catalizador antes de llegar a condiciones de operación insostenibles para la tecnología desplegada.

Se logró simular la operación de un reactor instalado y en funcionamiento para la obtención de metanol usando un catalizador del tipo Cu-Zn-Al. Se demostró que el uso de información de acceso libre en la literatura es una base de información sólida para ajustar el comportamiento real si se usa con criterios sanos y bien fundamentados que también son accesibles. Se logró simular el comportamiento de dicho reactor a lo largo de 5 años con el conocimiento de las condiciones operacionales a la fecha, utilizando la misma información cinética a lo largo de todo el periodo y aplicando solo un factor de actividad, cuyo cálculo se propone en función del tiempo en corrida. El modelo desarrollado tiene el potencial de permitir programar los cambios necesarios de variables operacionales para mantener el nivel de producción en línea con la sala de control, evitando los saltos observados con el protocolo en uso. Incluso permitiría programar con tiempo las operaciones de recambio del catalizador.

Agradecimientos

Los autores expresan su agradecimiento al Centro de Desarrollo Científico y Humanístico (CDCH) de la UCV por haber financiado el proyecto 08.00.6126.2005.

Referencias

Froment G. y Bischoff K. Chemical Reactor Analysis and Desig, New York, Wiley, 1979, pp. 347-359.         [ Links ]

Graaf G.H, Scholtens H., Stamhuis E.J., Beenackers A.A.C.M. Intra-Particle Diffusion Limitation in Low-Pressure Methanol Synthesis. Chemical Enginering Science, volumen 45, 1990: 773-783.         [ Links ]

Lommerts B., Graaf G., Beenackers A. Mathematical Modeling of Internal Mass Transport Limitations in Methanol Synthesis. Chemical Engineering Science, volumen 55, 2000: 5589-5598.         [ Links ]

Lovik Ingvild. Modeling, Estimation and Optimization of the Methanol Synthesis with Catalyst Deactivation, tesis (doctorado en ingeniería), Noruega, Norwegian University of Science and Technology, marzo 2001, pp. 14-26.         [ Links ]

Mathews J.y Fink K. Numerical Methods Using Matlab, 4a ed., Nueva Jersey, USA, Prentice Hall, 2004, pp. 430-436.         [ Links ]

Missen R., Mins C., Saville B. Chemical Reaction Engineering and Kinetics, New York, Wiley, 1999, pp. 64-80.         [ Links ]

Molina C. y Gabriel A. Desarrollo de un modelo matemático para la simulación de reactor de síntesis de metanol (grado de ingeniero químico), Caracas, Universidad Central de Venezuela, Facultad de ingeniería, 2010, pp. 10-17.         [ Links ]

Perry R.H y Green D.W. Perry Chemical Engineering Handbook, 8a ed., Nueva York, McGraw-Hill, 2008, pp. 174-181.         [ Links ]

Rosner D.E. Transport Processes in Chemically Reacting Flow Systems, Mineola, Nueva York, Dover Publications, 2000, pp. 27-44.         [ Links ]

Shahrokhi M. y Baghmisheh G. Modeling, Simulation and Control of a Methanol Synthesis Fixed-Bed Reactor. Chemical Engineering Science, volumen 60, 2005: 4275-4286.         [ Links ]

Shampine L.F. y Reichelt M.W. The MATLAB ODE Suite. SIAM Journal on Scientific Computing, volumen 18, 1997: 1-22.         [ Links ]

Valenzuela C.C. Química general. Introducción a la química teórica, Salamanca España, Ediciones Universidad de Salamanca, 1995, pp. 271.         [ Links ]

Vanden-Bussche y Froment G. A Steady-State Kinetic Model for Methanol Synthesis and the Water Gas Shift Reaction on a Commercial Cu/Zn/Al2O3 Catalyst. Journal of Catalysis, volumen 161, 1996: 1-10.         [ Links ]

Velardi S. y Barresi A. Methanol Synthesis in a Forced Unsteady-State Reactor Network. Chemical Engineering Science, volumen 57, 2002: 2995-3004.         [ Links ]

Este artículo se cita:

Citación estilo Chicago
Molina-Cárdenas, Gabriel Alberto, Andrés Emilio Rosales-Chinchilla, José Papa-Annibalini. Modelado matemático de un reactor de síntesis de metanol con interenfriamiento y desactivación de catalizador. Ingeniería Investigación y Tecnología, XV, 04 (2014): 505-515.

Citación estilo ISO 690 ]]> Molina-Cárdenas G.A., Rosales-Chinchilla A.E., Papa-Annibalini J. Modelado matemático de un reactor de síntesis de metanol con interenfriamiento y desactivación de catalizador. Ingeniería Investigación y Tecnología, volumen XV (número 4), octubre-diciembre 2014: 505-515.

Semblanza de los autores

Gabriel Alberto Molina-Cárdenas. Es ingeniero químico egresado de la Universidad Central de Venezuela,  ejerce su profesión en una compañía reconocida de Venezuela en el área de procesos. Es investigador del desarrollo de modelos y simulación de procesos industriales.

Andrés Emilio Rosales-Chinchilla. Es profesor asistente dedicado en exclusiva y jefe del Departamento de Termodinámica y Fenómenos de Transporte de la Escuela de Ingeniería Química de la Universidad Central de Venezuela. Cuenta con varias ponencias en congresos nacionales e internacionales.

José Papa-Annibalini. Es profesor titular jubilado por la Universidad Central de Venezuela, es investigador activo desde hace más de treinta años en el área de ingeniería. Cuenta con varias publicaciones en revistas nacionales e internacionales, así como ponencias en congresos locales e internacionales.