Modelado matemático de un reactor de síntesis de metanol con interenfriamiento y desactivación de catalizador

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.

Abstract

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.

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.

Reactor estudiado

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.

Figura 1.

Reactor de síntesis de metanol con interenfriamiento (Molina 2010)

Tabla 1.

Alimentación típica del reactor

Componente	Concentración (%mol)
CH4	15.11
CO	2.83
CO2	2.15
H2	79.86
N2	0.05

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 6×5 mm. La presión típica de operación ronda los 10 MPa y las temperaturas en los lechos están entre 503y 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:

Balance de energía:

donde:

Fi	= 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^3)
Θ	= fracción vacía del lecho una vez que contiene el catalizador (m^3/m^3)
D	= diámetro del lecho (m)
rmet	= velocidad de formación metanol según la reacción 8 (molmetanol/(kgcat.s))
rwgs	= velocidad de formación dióxido de carbono según la reacción 9 (molCO2/(kgcat.s))
ηmet	= factor de eficiencia para la reacción 8 (adimensional)
ηwgs	= factor de eficiencia para la reacción 9 (adimensional)
a	= factor de actividad del catalizador (adimensional)
T	= temperatura en el lecho (K)
Cpi	= capacidad calorífica de la especie “i” (Joule/(mol.K))
ΔHmet	= calor de reacción asociado con la reacción 8 (Joule/mol)
ΔHwgs	= 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

Para la síntesis de metanol se han propuesto varios modelos cinéticos. Para este trabajo se adoptó el propuesto por Vanden y Froment (1996), con base en el mecanismo que desarrollaron los autores se consideran las siguientes reacciones químicas y expresiones cinéticas.

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.

Tabla 2.

Constantes cinéticas y de equilibrio dados por Vanden Bussche y Froment (1996)

ki=A.e(ψBR.T)	A	B
k1 (mol/(Kg.s.bar^2))	1.07	36696
k2	3453.38	-
k3 (bar^-1/2)	0.499	17197
k4 (bar^-1)	6.62×10^−11	124119
k5 (mol/(Kg·s·bar))	1.22×10^−10	−94765
keqi=10(AT-B)	A	B
keq.met (bar^-2)	3066	10.592
keq.wgs	2073	2.029

Tabla 3.

Factores de ajuste cinético “ψ” por cada lecho del reactor

Lecho	ψ
1	0.8913
2	0.9187
3	0.9245
4	0.9384
5	0.9057

Modelo de difusión interna

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:

Ci	= concentración molar del compuesto “i” (mol/m^3)
kpseuMet	= constante pseudo-cinética para la reacción 8
kpseuWgs	= 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)
Defec_j	= difusividad efectiva del componente clave “j” en la matriz de catalizador (m^2/s)
Rp	= radio de la pastilla de catalizador (m)

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)
Kd	= factor de frecuencia de desactivación (4.39 E-3 días^–1)
Ed	= energía de activación de sinterización (91270 Joule/mol)
R	= constante universal de los gases (3.14 Joul/(mol.K))
ao	= factor de actividad de equilibrio (0.4 adimensional)
To	= 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.

Tabla 4.

Métodos usados para el cálculo de propiedades físico químicas y de transporte

Propiedad	Método
Fugacidad del component “i”	Ecuación de Berthelot(Valenzuela, 1995)
Capacidad calorífica del gas “i”	Gas ideal(Perry y Green, 2008)
Difusividad binaria entre la especie “i” y la “j”	Fuller-Schettler-Gidding(Perry y Green, 2008)
Difusividad de Knusend	Difusividad de Knusend(Missen et al., 1999)
Difusividad efectiva de “i” en el catalizador	Ecuación de Wilke(Froment y Bischoff, 1979)

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.

Para resolver el modelo se requirió una estructura donde cada submodelo se va solucionando en secuencia (figura 2). La simulación del reactor depende de la correcta solución e interacción entre los distintos submodelos. En el caso de las ecuaciones diferenciales ordinarias presentes en el modelo del lecho catalítico, estas presentaron problemas de rigidez, por lo que se usó un método adecuado para las mismas basado en una fórmula de diferenciación numérica (Shampine y Reichelt, 1997). Este método está incluido como una función en la herramienta de software matemático Matlab, donde se programó el modelo.

Figura 2.

Esquema de solución del modelo del reactor

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.

Tabla 5.

Comparación entre los resultados de la simulación y los datos de planta. Flujo y composición del gas efluente

		Año 1			Año 2			Año 3			Año 4			Año 5
Variable	Unidades	Planta	Modelo	Error rel.	Planta	Modelo	Error rel.	Planta	Modelo	Error rel.	Planta	Modelo	Error rel.	Planta	Modelo	Error rel.
Flujo	Kmol/h	80521.555	79155.120	1.697	82796.192	80464.507	2.816	81388.224	79684.357	2.094	80236.835	77937.860	2.865	81432.562	79069.297	2.902
CO	% conv.	82.237	87.951	6.948	81.448	87.483	7.409	80.390	86.460	7.551	79.727	84.443	5.915	78.762	83.787	6.380
CO2	% conv.	81.658	87.775	7.492	81.770	87.873	7.464	81.861	87.698	7.131	83.231	88.245	6.024	83.415	88.451	6.037
H2O	% conv.	93.765	94.938	1.251	93.703	94.928	1.307	93.703	94.717	1.082	92.141	94.076	2.100	92.407	93.844	1.555
CH4	%mol	16.543	16.687	0.867	16.239	16.378	0.856	16.280	16.397	0.723	15.966	16.038	0.452	15.987	16.043	0.355
N2	%mol	0.075	0.073	2.543	0.066	0.065	0.976	0.066	0.064	2.709	0.122	0.119	2.844	0.094	0.086	8.482
H2	%mol	76.017	75.789	0.300	76.358	76.127	0.302	76.274	76.087	0.245	76.839	76.739	0.130	76.922	76.927	0.006
Metanol	%mol	4.560	4.697	3.010	4.606	4.714	2.363	4.537	4.697	3.528	4.563	4.463	2.191	4.449	4.366	1.851

Figura 3.

Comparación del perfil de temperatura medido en planta y el predicho por el modelo para el primer año de operación

Figura 4.

Comparación del perfil de temperatura medido en planta y el predicho por el modelo para el tercer año de operación

Figura 5.

Comparación del perfil de temperatura medido en planta y el predicho por el modelo para el quinto año de operación

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/cm² (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.

Figura 6.

Comparación del desempeño del reactor según el modelo, manteniendo constantes las variables operacionales, y considerando la funcionalidad de la actividad con el tiempo en operación (□), con los datos reales de planta donde se observa el efecto de la manipulación de las variables operacionales (Δ)

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.

Conclusiones

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.