Optimización de un Reactor Catalítico de Lecho Fijo para la

Producción de Formaldehído, Basado en un Modelo

Bidimensional

José Santiago Rodríguez

Asesor: Dr. J.Gómez

Departamento de Ingeniería Química

Universidad de los Andes
Bogotá, Colombia
Junio, 2011

1
Agradecimientos

Primero, gracias a Dios por estar conmigo a lo largo de todo este proceso educativo y permitirme
concluir con un proyecto de grado tan exitoso como el presente. Gracias por cruzar en mi camino a
todas las personas importantes que me guiaron en este periodo educativo. En especial gracias a los
profesores Gustavo Orozco, por sus infinitas enseñanzas en ingeniería química y por su tan
apreciada amistad, Felipe Muñoz, por su apoyo en momentos difíciles y por brindarme las
herramientas necesarias para culminar el presente proyecto, Andrés Gonzales, por la
retroalimentación que me presto y finalmente gracias a Nicanor Quijano por su gran apoyo en este
último periodo de 6 meses.

Gracias a mi padre y a mi madre por darme la mejor educación y por hacer de mí el hombre que
soy. Sin ustedes nada de esto habría sido posible. Gracias a mi familia que siempre me ha apoyado,
a mis hermanos, Juliana y Sergio que han estado junto a mí y también gracias a Aura Tolosa y
Wilson Molina por su amistad incondicional.

“Es tan grande el placer que se experimenta al encontrar un hombre agradecido que vale la pena arriesgarse a no ser un
ingrato”. Lucio Anneo Séneca

2
Índice General

1) Notación 5
2) Introducción formaldehido en la industria 7
3) Estado del arte
a. Cinéticas reportadas para la oxidación de metanol a formaldehido. 8
b. Modelamiento de reactores Tubulares. 12
c. Clasificación de ecuaciones diferenciales parciales de segundo orden. 21
d. Métodos de resolución de ecuaciones diferenciales parciales. 22
4) Resultados y análisis de resultados
a. Análisis de estabilidad de la cinética de reacción. 28
b. Simulación en dos dimensiones. 36
c. Optimización económica 42
5) Conclusiones. 46
6) Bibliografía. 47

3
Tabla de Ilustraciones

Ilustración 1. Modelos de reactores catalíticos de lecho fijo ........................................................... 12

Ilustración 2. Función de suministro y consumo para un pellet catalítico isotérmico con cinética no
monotonica, extraído de Modelling, Simulations and optimizations of industrial fixed bed catalytic
reactors; p 145 .................................................................................................................................. 18
Ilustración 3. Perfiles cinética Tesser ............................................................................................... 30
Ilustración 4. Perfiles cinética Moustafa .......................................................................................... 31
Ilustración 5. Perfil de Temperatura obtenido con ode45 en Matlab ............................................... 32
Ilustración 6. Valores propios Moustafa........................................................................................... 33
Ilustración 7. Perfiles cinética Cozzolino. ........................................................................................ 34
Ilustración 8. Perfiles cinética Windes ............................................................................................. 35
Ilustración 9. Perfil de Presión de Metanol ...................................................................................... 36
Ilustración 10. Perfil de Presión de Formaldehído .......................................................................... 37
Ilustración 11. Perfil de Presión de Agua ......................................................................................... 37
Ilustración 12. Perfil de Presión de Monóxido de Carbono ............................................................. 38
Ilustración 13. Perfil de Temperatura............................................................................................... 38
Ilustración 14. Perfil de presión parcial de oxígeno......................................................................... 39
Ilustración 15. Perfil de presión parcial de metanol. ....................................................................... 39
Ilustración 16. Perfil de presión parcial de formaldehido. ............................................................... 40
Ilustración 17. Perfil de presión parcial de agua. ............................................................................ 40
Ilustración 18. Perfil de presión parcial de monóxido de carbono. ................................................. 41
Ilustración 19. Perfil de temperatura. ............................................................................................... 41
Ilustración 20. Perfil óptimo de Metanol .......................................................................................... 43
Ilustración 21. Perfil óptimo de Formaldehído................................................................................. 43
Ilustración 22. Perfil óptimo de Agua ............................................................................................... 43
Ilustración 23. Perfil óptimo de monóxido de carbono..................................................................... 43
Ilustración 24. Perfil óptimo de temperatura.................................................................................... 43
Ilustración 25. Óptimos locales. ....................................................................................................... 44
Ilustración 26. Perfil de formaldehído óptimo .................................................................................. 45
Ilustración 27. Perfil de Temperatura óptimo .................................................................................. 45

4
Nomenclatura

: Concentración.

Difusividad efectiva axial.

Difusividad efectiva radial.

Coordenada axial.

Coordenada radial.

Velocidad.

Tiempo.

Avance global de reacción.

Densidad del lecho catalítico.

Densidad de la mezcla catalítica.

Conductividad efectiva axial.

Conductividad efectiva radial.

Temperatura de la mezcla reactiva.

Temperatura de refrigerante.

Capacidad calorífica.

Número de componentes.

Número de reacciones.

Delta H de reacción.

: Velocidad de la iesima de reacción.

Fracción vacía.

5
Velocidad másica.

Viscosidad.

Diámetro de partícula catalítica.

Constante pre-exponencial.

Energía de activación.

Presión parcial de metanol.

Presión parcial de oxígeno.

Presión parcial de formaldehído.

Presión parcial de monóxido de carbono.

Presión parcial de agua.

Presión total.

Diámetro del reactor.

Radio del reactor.

Flujo molar.

Conversión.

Selectividad.

Rendimiento.

6
Introducción

La producción de formaldehído actualmente es de gran importancia para la industria mundial.

Dicho producto tiene principalmente 3 aplicaciones en la industria química; en la fabricación de
resinas para la elaboración de adhesivos de elementos de madera, como intermediario para la
producción de colorantes, papel, material fotográfico, vitaminas, EDTA, fertilizantes y como como
agente conservante en la industria cosmética [17]. Por lo tanto el formaldehído se constituye en
materia prima altamente demandada para la fabricación de productos de diferentes mercados y
puede ser clasificado como un producto intermedio en la industria petroquímica. Por esta razón a lo
largo del último siglo se han desarrollado una gran variedad de estudios del sistema reactivo para la
producción de formaldehído, con el propósito de obtener un procedimiento eficiente y seguro.
Actualmente el formaldehído se fabrica principalmente por oxidación parcial de metanol en
catalizadores de óxidos metálicos como el Fe-Mo o con catalizador de Ag [17]. Sin embargo la
oxidación parcial de este alcohol presenta serias complicaciones debido a la alta no linealidad de la
cinética de reacción, donde pueden presentarse fenómenos de multiplicidad de estados estacionarios
y puntos calientes en reactores multitubulares de lecho fijo usados en el proceso.

La no linealidad del mecanismo cinético de oxidación parcial de metanol a formaldehído hace de

este sistema muy complejo matemáticamente. La inestabilidad del sistema y posible multiplicidad
de estados estacionarios hacen atractivo desde el punto de vista ingenieril el estudio del mismo. Los
métodos numéricos deben ser correctamente empleados e inicializados debido a la inestabilidad ya
mencionada, adicionalmente los parámetros seleccionados para modelar un reactor de oxidación
parcial de metanol a formaldehído deben ser escogidos con mucho cuidado, debido a la sensibilidad
matemática de las ecuaciones a estos valores.

Con el objeto de determinar las condiciones de operación y diseño más favorables para la
producción de formaldehído a partir de la oxidación parcial de metanol en un reactor tubular, el
presente trabajo pretende realizar un estudio de la optimización de este sistema reactivo, analizando
los mecanismos cinéticos propuestos en la literatura con el propósito de posteriormente optimizar
un reactor, basado en un riguroso modelo bidimensional-2D, en donde sea posible determinar las
condiciones más favorables para la producción del formaldehído.

7
Modelos Cinéticos
Los principales modelos cinéticos propuestos para estudiar el sistema reactivo de oxidación
del metanol a formaldehído son el modelo de Larry C. Windes, Marvin J. Schwedock
establecido en 1987 y el modelo de Cozzolino.M, Tesse.R, DiSerio.M en el año de 2007,
Deshmukh S.A.R.K, Van Sint Annaland M en el 2005, Andreasen.A, Lynggard.H,
Stegelman.C, Stoltze.P en el 2003 y Tesser.R, DiSerio.M, Santacesaria. En general las
cinéticas reconocen la oxidación selectiva de metanol a formaldehído con catalizadores de
óxidos de metalicos y siguen un mecanismo redox.

Modelo de L.C. Windes et al. [17]

Larry C. Windes, Marvin J. Schwedock realizan el estudio de la oxidación parcial de

methanol a formaldehido sobre un catalizador de FeO3-MoO3 en un reactor multitubular
que opera a 523 K y 1.5 atm de presión. En su artículo presentan diferentes modelos
matemáticos que tienen en cuenta difusión axial y radial, además de que también
consideran resistencias en la interfaz catalizador-fluido. Las reacciones postuladas en el
artículo son:

Con entalpias de reacción ∆H = -37.9 Kcal/mol y -55.7 Kcal/mol respectivamente. Las

constantes pre-exponenciales y energías de activación reportadas en el artículo son:

Rxn Ai mol/(s*m3*atm1/2) EA (cal/mol)

1 6250 19000
2 5.6 16000
Tabla 1. Parámetros Cinéticos

En su artículo proponen que el modelo cinético que siguen las reacciones es:

(1)

(2)

8
 Modelo de M.Cozzolino et al. [20]

Este modelo fue desarrollado a partir de un análisis cinético de un mecanismo redox

propuesto para la oxidación de metanol con el uso de un catalizador de vanadio.
Experimentalmente estudiaron el efecto de la temperatura, el tiempo de contacto y la
proporción de reactivos alimentados a un reactor multitubular de acero inoxidable operado
isotérmicamente. Realizaron 3 corridas diferentes, variando la temperatura, el tiempo de
residencia y la cantidad de catalizador presente en los tubos. Las reacciones postuladas en
su artículo son:

3CH3OH + O2 → CH2(OCH3)2 +2H2O

CH2(OCH3)2 + H2O ↔ 2CH3OH + HCHO
CH3OH + O2 → CO2 +2H2O
CH3OH + H2CO → HCOOCH3 + H2O
Las cinéticas reportadas en su artículo son:

(3)

(4)

(5)

(6)

Por último los parámetros cinéticos fueron determinados a partir de una regresión no lineal
ajustada a los datos experimentales.

Constante Ln(A) EA (kcal/mol)

K1 23.09 ± 2.6 20.4 ± 2.3

K2 20.6 ± 1.1 11.1 ±1.0

K3 15.6 ±4.8 14.8 ± 4.2

K4 39.0 ± 3.2 39.0 ± 3.2

Tabla 2. Parámetros Cinéticos

9
 Modelo de R.Tesser et al. [7]

Tesser R y sus colaboradores a partir del estudio del factor de efectividad de la oxidación
parcial de metanol a formaldehido desarrollan un modelo cinético que considera el
siguiente mecanismo en un catalizador de molibdeno de hierro (Fe-Mo).

Estas reacciones siguen un mecanismo redox y la experimentación para desarrollar las

expresiones cinéticas se realizó a 1.1 atm y a un temperatura de 563K

(7)

Los órdenes de reacción a, b reportados en su artículo son a=1 y b=1/2, adicionalmente se

considera que la segunda reacción sigue una cinética de pseudo primer orden con respecto
al formaldehido. A continuación se presentan los parámetros cinéticos reportados.

Tabla3. Parámetros Cinéticos

Modelo de Deshmukh S.A.R.K

Deshmukh et al estudian la oxidación parcial de metanol a formaldehido en un reactor

diferencial a temperatura de 230°C con un catalizador de Fe-Mo. En su estudio se realizan
las siguientes observaciones:

1. La oxidación parcial de metanol sigue una cinética de Langmuir-Hinshelwood y

depende de la concentración de metanol.
2. Se observa un comportamiento de inhibición de la formación de formaldehido en
presencia de concentraciones considerables de agua, debido a que se genera un
fenómeno de competitividad por los sitios activos del catalizador.
3. Se observa una cinética de primer orden para la formación de monóxido de carbono
con respecto a la concentración de formaldehido.
4. Se presentan las expresiones cinéticas desarrolladas a partir de un mecanismo de
Langmuir-Hinshelwood, basados en la oxidación de formaldehido, dimetil éter y
dimetoximetano.

10
 Modelo de Moustafa [22]

Moustafa presenta un modelo cinético que considera el fenómeno reactivo-difusivo para la

oxidación de metanol a formaldehido. A diferencia de otros modelos presentados en la
literatura, el modelo de Moustfa [22],[23] se caracteriza por considerar las resistencias a la
transferencia de masa y energía en la interface catalítica. Para ello calcula el factor de
efectividad y lo incluye en las ecuaciones de balance de materia y energía al simular el
reactor. Las reacciones que considera este modelo son:

Las respectivas velocidades de reacción son:

(8)

(9)

(10)

(11)

Las constantes se estiman a partir de la ecuación de Arrhenius y los parámetros cinéticos

se presentan en la referencia [22]

11
Modelamiento de Reactores Tubulares

En el último siglo una gran variedad de estudios acerca del modelamiento de reactores han sido
realizados, entre estos se destacan los efectuados por Deans and Lapidus en 1960; Deasch and
Froment en 1971; Varma en 1981; Fogler 1986 Entre muchos otros autores. Actualmente se
consideran dos clases de modelos para reactores catalíticos de lecho fijo, tal como se observa en la
Ilustración 1 se clasifican principalmente debido a el número de fases que se consideran presentes
en el reactor. De dicha clasificación se habla de modelos Pseudo- Homogéneos y Heterogéneos, en
el primero de ellos se considera que las condiciones del bulk son las mismas que las que las
condiciones en la superficie catalítica, mientras que en el modelo heterogéneo si se tienen en cuenta
la resistencias debidas a la transferencia de masa y energía presentes en la capa límite del pellet
catalítico. El concepto más práctico para determinar si es correcto aplicar un modelo Pseudo-
Homogéneo, pues presentara predicciones exactas de los datos experimentales, es el factor de
efectividad ya que se considera que las condiciones del bulk iguales a las de la superficie catalítica
cuando η = 1. Además de considerarse un modelo debido al número de fases, también se tiene en
cuenta las dimensiones del reactor, es decir, que tan válido es asumir un modelo matemático que
considere dos dimensiones (radio y longitud) o simplemente describa el comportamiento a lo largo
del reactor ignorando las consideraciones radiales.

Ilustración 1. Modelos de reactores catalíticos de lecho fijo

Modelo Principal.
La clasificación de modelos presentada en la Ilustración 1 de obtiene a partir de la realización de
diferentes suposiciones que simplifican el modelo matemático para la descripción del fenómeno
reactivo que se presenta en un reactor tubular de lecho catalítico fijo. Dicho modelo se obtiene a

12
partir de la aplicación de las ecuaciones de continuidad y energía sobre un elemento diferencial en
el reactor. Siendo estrictos también debería incluirse en el análisis la ecuación de Navier-Stockes
para considerar variaciones de momentum en el elemento diferencial. Sin embargo dicho acople le
agregaría aun una mayor complejidad matemática al modelo. Por lo tanto en la literatura consultada
se encontró principalmente el modelo descrito por la ecuación de masa y energía que se presentan
en las ecuaciones 12 y 13.

(12)

(13)

Con condiciones de iniciales y de frontera

para todo r y t cuando

para todo z y t cuando r=0

para todo z y t cuando r=Rt

A partir de las ecuaciones 12 y 13 es posible bajo ciertas suposiciones desarrollar los modelos
ilustrados en la Ilustración 1. Sin embargo antes de entrar en detalle sobre las suposiciones
realizadas en cada modelo, se explicara que significa físicamente cada término de las ecuaciones de
masa y energía, con el propósito de darle al lector mayores herramientas para entender el modelo
matemático y las simplificaciones que implican las suposiciones realizadas en los modelos pseudo-
homogeneo y heterogéneo. El primer término de las ecuaciones se asocia a la difusión axial de la
mezcla reactiva, en la literatura consultada este término se refiere a la dispersión axial que se
presenta en el reactor. Normalmente por simplificación matemática para la resolución del sistema
de ecuaciones diferenciales parciales este término se desprecia debido a que el coeficiente de
difusión axial efectivo tiene un orden de magnitud muy pequeño en reactores industriales. El
segundo término es el laplaciano de la propiedad que se analice (T o Ci) y representa la difusión
radial en el reactor. El tercer término representa el fenómeno de convección de la mezcla reactiva y
por cuestiones de simplicidad se impone un flujo pistón a lo largo del reactor (velocidad constante).
El cuarto término es de gran importancia para estudios dinámicos de reactores tubulares, este
término considera el estado transitorio del sistema reactivo. Finalmente el ultimo termino se refiere
la generación o consumo de la propiedad que se analice (calor o Ci) y aquí se incluye la presencia
de un catalizador solido o pseudo-homogeneo tal como se verá a continuación.

Modelo unidimensional Pseudo-Homogeneo

El más básico de los modelos de reactores catalíticos de lecho fijo, es el modelo unidimensional
Pseudo-Homogéneo. En este modelo se asume que solo existen gradientes de concentración y

13
temperatura en la dirección axial del reactor, por lo tanto se ignoran la dispersión radial y los
gradientes debidos a la resistencia en la interface solido-gas. Las ecuaciones de conservación de
masa y energía se obtienen a partir de un balance en estado estacionario sobre un componente A en
la sección transversal del reactor tubular.

(15)

(16)

Con condiciones iniciales

El sistema de ecuaciones resultante modela únicamente la fase fluida dentro del reactor e ignora la
presencia de partículas catalíticas solidas dentro del reactor. La principal suposición de este modelo
radica en no considerar la transferencia de masa y energía en la capa limite presente en la interface
solido gas, pues se asume que la resistencia a la transferencia es despreciable con el propósito de
simplificar la resolución del sistema.

Modelo unidimensional Heterogéneo

A diferencia del modelo pseudo-homogeneo el modelo heterogéneo considera gradientes de

concentración y temperatura en la interface solido-gas. Esta distinción permite que con el modelo
heterogéneo se obtengan datos mucho más exactos con respecto a los datos experimentales. Más
aun, el modelo heterogéneo es esencial para el estudio de estabilidad del reactor, pues el análisis de
los gradientes en la interface es fundamental para determinar posibles casos de bifurcaciones. Como
se vio anteriormente el cálculo del factor de efectividad considera diferencias en la superficie
catalítica e interior del pellet, y estas diferencias pueden ser cuantificadas a partir del modelo
heterogéneo. A partir de un balance de masa y energía sobre una sección transversal se obtienen las
siguientes expresiones.

(17)

(18)

(19)

(20)

Donde las ecuaciones 17 y 18 representan el fluido y las 19 y 20 el sólido, las primeras dos
expresiones tienen como condiciones iniciales cuando z=0, C=C0 y T=T0. Es de resaltar que en este
modelo los términos de generación se consideran únicamente en las ecuaciones de modelamiento
del pellet, las cuales se acoplan a las ecuaciones de fluido para resolver un sistema algebraico
diferencial. Por lo tanto la complejidad de este esquema es representativamente más alta en
comparación al modelo unidimensional pseudo-homogeneo. Adicionalmente si se desea ser
riguroso en la simulación del fenómeno físico, se debe incorporar el cálculo del factor de
efectividad en cada paso de integración del DAEs, con el propósito de considerar fenómenos de

14
transporte dentro de la partícula catalítica. Dicho cálculo del factor de efectividad de discutirá
posteriormente.

Modelo Bidimensional Pseudo-Homogéneo

Debido a que los modelos presentados previamente no consideran resistencias a la transferencia de

energía y masa en dirección radial, se desarrolló un modelo que permite predecir perfiles de
temperatura y conversión en secciones transversales. Las ecuaciones de continuidad y energía para
un componente reactivo A se presentan a continuación:

(21)

(22)

Con condiciones de frontera C=C0 y T=T0 cuando z=0, en r=0 y r=Rt, en r=0 y una
condición de robin en r=Rt, .

Modelo Bidimensional Heterogéneo

Con el propósito de considerar variaciones energéticas y de masa en dirección radial del catalizador,
en donde se tenga en cuenta la resistencia de transferencia debida a la capa limite presente en la
interface solido-vapor entre el catalizador y bulk, en 1971 se propuso el siguiente modelo:

o Para el fluido

(23)

(25)

o Para el catalizador

(27)

(28)

Con condiciones de frontera:

15
 Factor de efectividad

El cálculo del factor de efectividad (η) se desarrolló en un principio como un concepto practico que
permitiera determinar las velocidades de reacción en un pellet catalítico de los compuestos
involucrados en un sistema reactivo. En términos matemáticos el factor de efectividad es:

(29)

Si se escribe la ecuación 29 con el fin de definir la velocidad de reacción en función de las

condiciones de superficie se obtiene

(30)

Entonces es obvio que la velocidad de reacción sobre todo el pellet está relacionada con las
condiciones de superficie y también con el factor de efectividad. Pero debido a que las condiciones
de superficie dependen de las condiciones del bulk (resistencia convectiva), es válido concluir que
la velocidad de reacción es función de estas últimas y del factor de efectividad. Por lo tanto la
evaluación de η tiene como objetivo poder evaluar rp.

Para entender un poco más el concepto de factor de efectividad suponga que tiene una reacción con
cinética de primer orden (A B) y que dicha reacción sucede sobre un pellet catalítico de
geometría esférica de radio r, al realizar un balance de masa detallado obtenemos:

(31)

(32)

(33)

(34)

16
 (35)

Resolviendo la ecuación 35 con condiciones de frontera r=0; y r=R; C=Cs se obtiene el

siguiente perfil de concentración

Donde (36)

Manipulando la ecuación 36 y realizando de nuevo un balance es posible llegar a la siguiente

expresión [14]

(37)

Si se analiza detalladamente esta última expresión, es de notar que el factor de efectividad, para un
sistema reactivo isotérmico y de primer orden, solo es función del módulo de thiele ф y que para
thiele pequeños η tiende a 1 lo que indica que la velocidad de reacción predomina sobre la
trasferencia de masa y por ende en un sistema reactivo isotérmico donde se obtienen valores de η
iguales o cercanos a 1, la velocidad de reacción se considera como la etapa limitante. Para casos
contrarios en donde se obtienen ф>>1 el limitante es la transferencia de masa pues predomina sobre
el paso químico. Pero la pregunta que surge es cuando se consiguen módulos de Thiele pequeños, la
respuesta es que este comportamiento sucede cuando se utilizan pellets pequeños, con difusividad
efectiva de gran magnitud o cuando simplemente la reacción es muy lenta.

Para el caso de un sistema no isotérmico el cálculo del factor de efectividad es parecido al

procedimiento realizado anteriormente, sin embargo además de realizar un balance de materia se
debe hacer un balance de energía y resolver un sistema de ecuaciones diferenciales, debido a que
ambos balances están acoplados por la ecuación de Arrhenius. La solución a este sistema EDP
(Ecuaciones Diferenciales Parciales) se presenta en la siguiente ecuación y es notable que se asume
una cinética lineal. A partir de este desarrollo matemático es posible demostrar que η es función de
tres parámetros adimensionales.

o Modulo de Thiele
(38)
o Modulo de Prater
(39)
o Numero de Arrhenius
(40)

17
El conocimiento de estos 3 parámetros ha permitido estudiar más a fondo el concepto de factor de
efectividad, se ha observado que para reacciones exotérmicas (β+) el incremento en la velocidad de
reacción acompañado del aumento en la temperatura puede balancear la baja velocidad de remoción
debido a la caída en la concentración del reactivo, obteniendo η>1 lo cual puede significar la
desactivación del sitio activo del catalizador debido a las altas temperaturas o también puede
evidenciar bajas selectividades en consecuencia a la razón de las constantes cinéticas y al aumento
de la temperatura.

Cuando el sistema reactivo presenta una cinética no lineal, el cálculo del factor de efectividad
resulta más complicado, Para casos como los desarrollados anteriormente, es válido emplear un
modelo de resistencias en serie para describir los fenómenos de transferencia dentro del pellet,
mientras que para cinéticas no lineales no es correcto aplicar dicho modelo. Como ejemplo
considere la siguiente velocidad de reacción

(41)

Ahora si se realiza un balance de masa sobre un pellet no poroso en estado estacionario se obtendría

(42)

La solución analítica de la ecuación (42) no proporcionará mucha información acerca del problema,
sin embargo resolver esta expresión gráficamente, para obtener CA para diferentes valores de
CAB, puede ser más conveniente en análisis del cálculo de factores de efectividad en sistemas no
lineales. En la Ilustración (2) se puede observar dos funciones, la primera de ellas es denominada la
función de suministro de reactivos S(CA) y es equivalente al lado izquierdo de la ecuación (42), de
forma similar la función de consumo de reactivos C(CA) es igual al lado derecho de la misma
ecuación. Al graficar estas dos funciones versus CA, se obtiene la Ilustración (2) en donde es claro
que las intersecciones son la solución al balance de materia.

Ilustración 2. Función de suministro y consumo para un pellet catalítico isotérmico con cinética
no monotonica, extraído de Modelling, Simulations and optimizations of industrial fixed
bed catalytic reactors; p 145

18
En el caso en que CAB=CAB1, se obtengan valores de CA=CA1, el cual es muy cercano a cero,
indica que la etapa controlante es la transferencia de masa ya que el catalizador es no poroso. Por
otro lado cuando CAB=CAB2, se obtienen valores de CA=CA2 cercanos a CAB2 lo que hace
referencia a que la cinética controla el proceso. Sin embargo al analizar el caso en que CAB=CAB3
se encuentra que la solución puede ser cualquiera de las tres raíces de la ecuación cubica y por ende
no es fácil determinar que fenómeno es el limitante.

En cuanto al cálculo de los factores de efectividad, se deben analizar los tres casos que se
describieron en el párrafo anterior:

1. Para CAB=CAB1

(15)
2. Para CAB=CAB2

(16)
3. Para CAB=CAB3

(17)

(18)

(19)

Es importante notar que las conclusiones referentes a las etapas controlantes de cada caso se
evidencian también en el valor del factor de efectividad, cuando el factor de efectividad es mayor
que la unidad, el proceso se ve limitado por el paso cinético, mientras que cuando es menor a uno,
la etapa controlante es la transferencia de masa. En la siguiente sección se analizara el fenómeno de
multiplicidad de estados estacionarios, no obstante cabe resaltar que para el caso 3, cuando
CAB=CAB3, existe bifurcación en el sistema.

Calculo de coeficientes de transporte

Adicional al desarrollo del sistema de ecuaciones diferenciales parciales que propone el modelo
bidimensional pseudo-homogeneo, se deben estimar los coeficientes de transporte de la mezcla
reactiva. La difusividad y conductividad efectivas radiales son parámetros muy importantes para el
análisis del sistema reactivo. A pesar de ser simplemente una constante, detrás de ese coeficiente se
resume un análisis de los fenómenos de transporte que suceden dentro del reactor. Por tal razón es
de vital importancia ser riguroso en el cálculo de dichos parámetros para simular correctamente el
reactor. Para tal propósito se recomienda la ecuación de Wilke

19
 (43)

En esta ecuación se obtiene un promedio de las diferentes combinaciones binarias. El caculo de

estas difusividades binarias se realiza a partir de la ecuación de Fuller-Schettler-Giddings

(44)

Donde

(45)

Para la ecuación de Fuller debe ingresarse la temperatura en K y la presión en bar.

Dado que la ecuación de Wilke es función de la composición de cada componente, a lo largo del
reactor se tendrán diferentes difusividades de mezcla según sea la posición radial y axial. Por tal
razón debe calcularse la difusividad efectiva para cada componente en la mezcla en cada nodo
discretizado del sistema de ecuaciones diferenciales.

20
Clasificación de Ecuaciones Diferenciales [1][2][5]
Las ecuaciones diferenciales parciales se clasifican según tres criterios. EL primer criterio se
refiere a la linealidad de la ecuación diferencial, en segundo lugar se encuentra el orden de la
ecuación y finalmente las ecuaciones diferenciales se clasifican según las condiciones de frontera a
las que están sujetas. Considere los siguientes 3 ejemplos

Orden Ecuación diferencial Linealidad

1° Lineal

2° Cuasi-Lineal

3° No-Lineal

Para el caso concreto de una ecuación diferencial de segundo orden se tiene la forma canónica

(46)

Donde los coeficientes pueden ser constantes o funciones de x,y y dependiendo del valor que
tomen estas funciones las ecuaciones pueden clasificarse en:

Elíptica si ,ej: ecuación de Laplace

Parabólica si , ej: ecuación de Fourier
Hiperbólica si , ej: ecuación de onda

Finalmente las condiciones de frontera se clasifican en

1. Condiciones de Dirichlet: los valores de la variable dependiente son fijos o son funciones
de las variables independientes.

2. Condiciones de Newman: La derivada de la variable dependiente evaluada en un punto es

constante.

3. Condiciones de Robin: la derivada de la variable dependiente en un punto es igual a una

función de la variable dependiente.

21
Resolución de Ecuaciones Diferenciales Parciales

Existen diferentes métodos para la resolución de sistemas de ecuaciones diferenciales (PDEs), en la

presente sección se explicaran tan solo unos cuantos fundamentos con el propósito de motivar al
lector a indagar y profundizar en el campo de métodos numéricos empleados para resolver PDEs.
Dichos sistemas de ecuaciones son frecuentes en problemas científicos e ingenieriles, sin embargo
el enfoque de este documento será únicamente hacia un problema típico de ingeniería química como
lo es el caso del modelamiento de un reactor tubular. Como se discuto previamente los modelos
matemáticos para la simulación de reactores tubulares pueden consistir en sistemas de ecuaciones
diferenciales parciales según las suposiciones que se realicen, para resolver dicho sistema se
proponen los siguientes métodos.

Métodos de colocación

Uno de los métodos más recomendados para resolver sistemas de ecuaciones como los descritos en
la sección del modelamiento de reactores es el método de colocación [10]. Con el objetivo presentar
una explicación sencilla y concreta del método de resolución, suponga que se desea resolver la
siguiente ecuación diferencial ordinaria de segundo orden.

(47)

Con condiciones de frontera

en x=0

en x=1

Como es de notar esta ecuación es un balance de masa adimensional de un una especie reactiva en
un pellet catalítico esférico. Si se supone una solución polinomial de la forma

(48)

Las condiciones de frontera deben cumplirse para el polinomio que se asuma. Posteriormente se
determina la función R1(x) al sustituir el polinomio en la ecuación diferencial

(49)

Para determinar el valor de a1 se resuelve

(50)

Donde la función de residuos ponderados W(x) se define según el método de resolución que se
desee emplear. Existen 4 métodos de resolución propuestos para esta ecuación.

22
 1. Método de colocación: se busca que sea igual a cero, donde xi es cualquier
punto interior.
2. Método de subdominio: se toma la función de residuales ponderados W(x)=1.
3. Método de Garlekin: se toma la función de residuales ponderados .

4. Método de mínimos cuadrados: se toma la función de residuales ponderados .

El método más sencillo es el método de colocación debido a que se evita el cálculo de la integral
pero precisamente es el más inexacto. Sin embargo el algoritmo para solucionar la ecuación
diferencial por Garlekin o mínimos cuadrados adquiere mayor exactitud pero más peso
computacional. Al resolver la ecuación se obtiene [15]

Método Expresión para a1 a1( =2) y1(x)

Colocación

Subdominio

Garlekin

Mínimos cuadrados

Tabla 4. Métodos de residuales ponderados.

En general para polinomios de orden mayor se propone [15]

(51)

Donde Ti es una función de ensayo y se ajusta T0 tal que se cumplan las condiciones de frontera.
Para este caso se generaliza la ecuación de residuales ponderados a un sistema de ecuaciones no
lineal que al resolverse con cualquiera de los 4 métodos descritos previamente es posible obtener el
valor de los coeficientes ai. Para el caso de ecuaciones diferenciales parciales es posible aplicar los
métodos 1-4 tras discretizar la ecuación diferencial en una coordenada, para más detalles se
recomienda revisar [15]

Diferencias finitas

Otro método de resolución de ecuaciones diferenciales muy utilizado es la discretización de las

ecuaciones diferenciales en diferencias finitas. Dichas diferencias finitas consisten en expresar
algebraicamente las derivadas parciales a partir de las series de Taylor. En la tabla X se presentan
las principales aproximaciones encontradas en la literatura [1 para discretizar derivadas.

23
 Primera Segunda
Aproximación
Derivada Derivada

Por diferencias
hacia
Adelante

Por diferencias
hacia
atrás

Por diferencias
centradas

Tabla 5. Aproximaciones por diferencias finitas.

Con el propósito de ejemplificar como se discretiza una ecuación diferencial con las
aproximaciones presentadas en la tabla anterior, considere la ecuación de Fourier que describe la
transferencia de calor en una placa plana [1]:

(52)
Con condición inicial

para todo x cuando t=0

Con condiciones de frontera

para todo t cuando x=0

para todo t cuando x=L

Para construir la ecuación en diferencias finitas aproximamos la derivada temporal con diferencias
hacia adelante y la derivada espacial con diferencias centradas

(53)

24
Donde el subíndice i se refiere al avance espacial, j al avance temporal, k el paso en el tiempo y h el
paso en la dirección x. Al realizar esta discretización se pasa de un problema diferencial a un
sistema de ecuaciones algebraico donde la ecuación anterior describe el comportamiento de la
temperatura en los nodos internos de la placa. Para las fronteras se realiza el mismo procedimiento
y se discretizan las derivadas en caso de tener condiciones de Neuman o Robin.

Sin embargo este procedimiento suele tener problemas de estabilidad numérica y se debe garantizar
un principio de estabilidad que relaciona el paso en x y el paso en t . Por lo tanto se han
desarrollado aproximaciones en diferencias finitas mucho más confiables y que no se encuentran
limitadas por el principio de estabilidad mencionado. La aproximación de Crank Nicholson presenta
la ventaja de que es incondicionalmente estable y hay una extensa investigación de esta en la
literatura [10][1]. Las expresiones para aproximar con Crank Nicholson son:

(54)

(55)

(56)

(57)

Siguiendo el mismo procedimiento realizado anteriormente se encuentra una expresión homologa a

la ecuación 53 pero incondicionalmente estable para describir la temperatura de los nodos internos
de la malla. Para la resolución del sistema de ecuaciones algebraico resultante se implementan
diferentes métodos numéricos que serán descritos en la próxima sección.

Métodos de resolución de sistemas de ecuaciones algebraicas no lineales

La resolución de sistemas de ecuaciones diferenciales parciales como según el método

seleccionado, resulta en la solución de un sistema de ecuaciones algebraico, generalmente no lineal,
que puede ser tratado con métodos numéricos fundamentados en teoremas de convergencia. Por tal
motivo se considera de gran importancia discutir algunos de estos métodos, para familiarizar al
lector con los procedimientos y algoritmos utilizados en el caso de estudio. Específicamente en
problemas de ingeniería química los sistemas de ecuaciones algebraicos que resultan de discretizar
ecuaciones diferenciales parciales se caracterizan por ser no lineales [3]. Por tal razón esta sección
se centrara en explicar dos métodos clásicos para resolver estos sistemas.

o Método de Newton-Raphson

El algoritmo de Newton Raphson que se basa en la continuidad de las funciones a resolver, es uno
de los más conocidos y útiles para resolver sistemas de ecuaciones no lineales [1], [8]. Siguiendo la

25
línea explicativa que se ha desarrollado en este documento, se realizara una breve descripción a
partir de un ejemplo concreto, con el propósito de aclarar el algoritmo.

Suponga que se desea resolver el sistema de ecuaciones

Al expresar estas funciones en serie de Taylor alrededor de un punto se obtiene

Suponiendo que es una solución al sistema y si se consideran pequeños cambios

próximos al punto de solución, se obtiene la siguiente formula de recurrencia

Por tanto si la matriz jacobiana es invertible es posible despejar para proporcionar la

siguiente aproximación a la solución del sistema. En un esquema el algoritmo seria [1]:

Paso 1: Evaluar la función

Paso 2: Evaluar la matriz Jacobiana

Paso 3: Calcular el paso de iteración, resolviendo el sistema lineal

Paso 4: Calcular el siguiente punto

(58)

26
Es de notar que esta fórmula de recurrencia es la generalización de la formula iterativa para
funciones univariables. El método al igual que en el caso univariable, presenta inconsistencias
cuando la matriz jacobiana no es invertible y debe aplicarse algún ajuste al algoritmo para su
correcto funcionamiento. Sin embargo es un método confiable y software ya desarrollado como
Matlab y Mathematica de Wolfram permiten su fácil implementación.

o Método de relajación

Otro método muy común para la resolución de sistemas de ecuaciones algebraicas no lineales, es el
método de relajación. Este método iterativo consiste en encontrar una función tal que

Por cuestiones de estabilidad se introduce un parámetro de relajación, escalar que oscila entre cero y
la unidad, que disminuye el paso de iteración. Al incluir dicho parámetro la fórmula de recurrencia
para el caso de una función univariable

La anterior formula de recurrencia es fácil de generalizar a un sistema de ecuaciones multivariables.

Sin embargo este método es muy instable y depende significativamente del punto inicial de
iteración para converger, por tal razón debe seleccionarse con cuidado el punto de partida y
adicionalmente se recomienda utilizar parámetros de relajación cercanos a la unidad.

27
Análisis Cinético
Con el propósito de analizar las diferentes cinetecas reportadas en la literatura para la
deshidrogenación oxidativa del metanol a formaldehido se simulo un reactor unidimensional
siguiendo el modelo psudo-homogéneo descrito en secciones anteriores. El estudio consistió
principalmente en analizar los perfiles de concentración y temperatura para cada una de las cinéticas
trabajadas y adicionalmente para algunos casos concretos se analizaron problemas de rigidez e
inestabilidad en las ecuaciones de balance de materia y energía. Las ecuaciones de balance para el
modelo matemático fueron:

(59)

(60)

(61)

El cálculo de propiedades de mezcla se realizó a partir de ecuaciones ideales.

Por lo tanto estas tres propiedades se recalculan en cada paso de integración para considerar
variaciones en la temperatura y composición de la mezcla reactiva a lo largo del reactor.

Cinética de R.Tesser, M. Diserio y E. Santacesaria

Empleando la cinética propuesta por el Doctor Tesser y sus colegas se resolvió un sistema de 7
ecuaciones diferenciales ordinarias acopladas. Para solucionar dicho sistema se utilizó un
método explícito de integración como Runge Kutta cuarto orden implementado en el lenguaje
de programación c++. En la tabla 5 se presentan los parámetros cinéticos utilizados en la
simulación.

28
 Parámetros Cinéticos
K1 = 1.85e9 exp(-24100/RT) kmol/(h*kg*atm)
K2 = 2.98e5 exp(-15800/RT) kmol/(h*kg*atm)
K3 = 2.42e5 exp(-16050/RT) kmol/(h*kg*atm)
ΔH1= -37.9 kcal/gmol
ΔH2= -55.7 kcal/gmol
Tabla 5. Parámetros Cinéticos Tesser

Donde la energía de activación se encuentra en kcal/kmol. En este punto es importante aclarar que
la cinética utilizada para simular el sistema reactivo no considera la adsorción del agua en el
catalizador de molibdato de hierro, pues al incluir este término en la velocidad de reacción 1 se
presentan problemas con los parámetros de adsorción dado que no predicen datos reales. Por lo
tanto se simulo el sistema con la cinética propuesta inicialmente en [20].

Parámetro Valor
2
U (kcal/(m s K)) 0.017197
dt (m) 0.015
Vs (m/s) 2.1794
L (m) 1.2
dp (m) 0.003
ε 0.45
T0 (k) 523
P0 (atm) 1.25
Tr (k) 523
Tabla 6. Parámetros de diseño

En la tabla 6 se presentan los parámetros de diseño utilizados en la simulación de la cinética del

profesor Tesser. Los resultados obtenidos para el modelo unidimensional descrito anteriormente se
pueden observar en la Ilustración 3. Es de notar que el reactor presenta un punto caliente
aproximadamente a los 20 cm de longitud y que para ese punto se agota casi en su totalidad el
reactivo limite. Además es de destacar que tras realizar diferentes corridas variando ciertos
parámetros de diseño el sistema siempre presento resultados confiables y estables.

29
 Ilustración 3. Perfiles cinética Tesser

Al comparar los resultados con los obtenidos por el profesor Tesser y sus colegas en su artículo, se
observan diferencias en cuanto a la localización del punto caliente. Los resultados reportados
muestran que el punto caliente se encuentra alrededor de los 7.5 cm, esta diferencia se atribuye
principalmente al cálculo del factor de efectividad que es incluido en la simulación de Tesser y
adicionalmente debido a que en su investigación se empleó una cinética que considera la absorción
del agua en el catalizador.

Cinética de T. Moustafa, M. Abou-Elreesh y S.Fateen

Para el modelo del profesor Moustafa se resolvió un sistema de 9 ecuaciones diferenciales

ordinarias no lineales de primer orden. Al igual que para la cinética anterior se empleó inicialmente
el método de Runge Kutta cuarto orden pero al observar ciertos problemas de estabilidad de las
ecuaciones diferenciales cuando alguna de las variables dependiente se aproximaba a cero se realizó
un análisis de rigidez y estabilidad del sistema. Los parámetros cinéticos y de diseño para este
modelo se presentan en las tablas 7,8.

Reacción Ei (KJ/mol) K0j

(kmol/(h*kg*atm))
1 79.5 3.96 x 106
2 66.9 7.85 x 103
3 50.2 1.32 x 102
4 62.8 7.74 x 104
Tabla 7. Parámetros cinéticos Moustafa

Parámetro Valor
U (kcal/(m2 s K)) 0.017197
dt (m) 0.021
Vs (m/s) 2.1794
L (m) 1.2
T0 (k) 500
P0 (atm) 1.25
Tr (k) 500

30
Para este modelo los perfiles obtenidos con el método explicito se presentan en la Ilustración 4. Sin
embargo al simular el reactor con una temperatura de alimentación superior a 500 K se presentan
problemas de estabilidad a los 80 cm de longitud, pues algunas variables se van a infinito. Por tal
razón se programó en leguaje c++ un método explicito más robusto como es el caso del método del
trapecio pero se siguieron presentando los mismos problemas. Por lo tanto se decidió utilizar un
método aún más robusto y para ello se simulo el reactor en Matlab con un ode15s pero nuevamente
se presentaron problemas.

Ilustración 4. Perfiles cinética Moustafa

Al comparar los resultados con los presentados por el profesor Moustafa en su artículo [22], se
observan tendencias similares en cuanto a los perfiles reportados. Aunque existe una diferencia a
destacar y es el que el perfil de temperatura contiene un punto de inflexión que en nuestro caso no
se observa, lo cual puede atribuirse al cálculo del factor de efectividad y adicionalmente a que el
modelo desarrollado por el profesor Moustafa es un modelo 2D que considera dispersión radial y
por ende puede verse afectada la tendencia de la temperatura a lo largo del reactor. Sin embargo al
correr la simulación en Matlab, (un software que permite manipular resultados imaginarios) con una

31
longitud equivalente a la del reactor simulada en el artículo del profesor Moustafa, con un ode45 se
obtiene la Ilustración 5 en donde se observa el punto de inflexión mencionado anteriormente. No
obstante al analizar los resultados numéricos obtenidos se encontró que a partir de los 57 cm de
longitud se obtienen números complejos (omega del orden de 10-4) de los cuales solo se grafica la
parte real con la herramienta computacional. Adicionalmente se observa que el punto caliente
presentado en el artículo equivale a 590 K mientras que el nuestro equivale a 542 K lo cual se debe
a que no se simula el reactor con la misma composición de entrada ni el mismo coeficiente global
de transferencia de calor, dado que no se encuentran reportados estos valores en la bibliografía.

Perfil de Temperatura Moustafa

545

540

535

530
Temperatura (K)

525

520

515

510

505

500
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Longitud (m)

Ilustración 5. Perfil de Temperatura obtenido con ode45 en Matlab

Para analizar en detalle los resultados imaginarios mencionados previamente, se realizó un análisis
de los valores propios del sistema de ecuaciones diferenciales ordinarias y los resultados se
presentan en la Ilustración 6. Se destaca que solo en las primeras iteraciones los valores propios del
balance de dióxido de carbono toman valores positivos. Desde este punto de vista el sistema se
consideraría estable y no rígido debido a que ningún valor propio es mucho mayor que los demás,
no obstante algunos valores propios tienen parte real positiva cercana a cero, de lo cual se concluye
que el sistema presenta instabilidad numérica. En conclusión se encuentra en esta cinética varias
complicaciones a la hora de pensar en simular un reactor en 2 dimensiones y más aún si la idea es
realizar una optimización económica del reactor seria probable encontrar problemas en los procesos
iterativos.

32
 Estabilidad y Rigidez
0.1
Lambda1
0.05 Lambda2
Lambda3
0 Lambda4
Lambda5
-0.05
Lambda6
Lambda7
-0.1

Valor Propio
-0.15

-0.2

-0.25

-0.3

-0.35

-0.4
0 5 10 15 20 25 30 35 40 45 50
Iteración

Ilustración 6. Valores propios Moustafa

Cinética de M. Cozzolino, R. Tesser, M. Di Serio, P. D’Onofrio y E. Santacesaria.

En el 2007 el profesor Cozzolino y sus colegas publican una cinética para la oxidación de metanol a
formaldehido en un catalizador de soporte de vanadio. Para la simulación de este modelo se
resolvió un sistema de 10 ecuaciones diferenciales ordinaria con Runge Kutta cuarto orden
programado en lenguaje c++. En las siguientes tablas se presentan los parámetros cinéticos
repostados para este modelo y las especificaciones diseño bajo las cuales se simulo el reactor.

Constante Ln(A) o Ln(b0) EA o ΔHrxn (Kcal/mol)

k1 23.9 20.4
k2 20.6 11.1
k3 15.6 14.8
k4 39.2 27.8
kox 24.3 19.0
bM -28.8 -28.9
bW -23.2 -25.3
Tabla 9. Parámetros Cinéticos Cozzolino

Parámetro Valor
U (kcal/(m2 s K)) 0.017197
dt (m) 0.021
Vs (m/s) 2.1794
L (m) 0.9
dp (m) 0.003
ε 0.45
T0 (k) 450
P0 (atm) 1.25
Tr (k) 450
Tabla 10. Parámetros de diseño del reactor.

33
Los perfiles de concentración, presiones parciales y presión total en el reactor se presentan en la
Ilustración 7.

Ilustración 7. Perfiles cinética Cozzolino.

Los resultados obtenidos son coherentes con los obtenidos en cinéticas anteriores pero al igual que
el modelo de Moustafa este sistema de ecuaciones diferenciales presenta inestabilidad y
adicionalmente se considera un problema que la cinética modele 7 sustancias, dado que aumentaría
la complejidad para resolver un problema de optimización debido al número de variables a manejar,
ya que al discretizar el sistema de ecuaciones y resolver el problema con programación no lineal
tendrían que manipularse una magnitud considerable de variables.

34
 Cinética de L. Windes, M. Schwedock y W. Ray

En el caso de la cinética propuesta por el profesor Windes en 1989 se resolvió un sistema de 7

ecuaciones diferenciales en c++ con el método explícito de Runge Kutta y los perfiles obtenidos se
presentan a continuación.

Ilustración 8. Perfiles cinética Windes

Para concluir, tras analizar las diferentes cinéticas encontradas en la literatura se decidió trabajar
con la cinética del profesor Tesser debido a la estabilidad numérica que presenta el sistema, el
número de sustancias que involucra, disponibilidad de datos en la literatura [17] y adicionalmente
debido a la fecha de publicación del artículo de Tesser.

35
Simulación en Dos Dimensiones
Una vez seleccionada la cinética de estudio, fue posible realizar la simulación del reactor en dos
dimensiones. Para resolver el sistema de ecuaciones diferenciales parciales se utilizaron 2 de los
métodos descritos en la sección de resolución PDEs. El primero de ellos consistió en resolver
sistema de ecuaciones algebraico, resultante del proceso de discretización con aproximaciones de
Crank Nicholson, con el método de Newton Raphson en MATLAB. La simulación se realizó con
los parámetros presentados en la Tabla 11.

Parámetro Magnitud
0.017197
0.015
2.1794
0.7
0.003
0.45
523
1.25
523
Tabla 11. Parámetros simulación 2D

Se discretizó el sistema con una malla de 9000 nodos en MATLAB y los perfiles de presión parcial
y temperatura se presentan a continuación

Ilustración 9. Perfil de Presión de Metanol

36
Ilustración 10. Perfil de Presión de Formaldehído

Ilustración 11. Perfil de Presión de Agua

37
Ilustración 12. Perfil de Presión de Monóxido de Carbono

Ilustración 13. Perfil de Temperatura

38
Los resultados presentados en las Ilustraciones 9-13 se obtienen de una simulación del reactor con
transferencia de calor en la pared, con un fluido a 523 K. Las condiciones de alimentación son
equivalentes a las utilizadas en el modelo unidimensional de Tesser, presentado previamente en la
sección de análisis cinético. Se observa en los 6 perfiles que en la frontera se obtiene un
perfil equivalente al perfil unidimensional obtenido al resolver el sistema de 7 ecuaciones
ordinarias. También se evidencia la presencia del punto caliente en los primeros cm del reactor pero
a diferencia del caso unidimensional la temperatura alcanza los 798.4 K.

Con el objetivo de validar los resultados obtenidos con el método simultáneo, se realizó un
programa en lenguaje c++ para resolver el sistema de ecuaciones algebraicas con un método
explícito de relajación y verificar la validez de los resultados obtenidos en el programa realizado en
MATLAB. Para ello, se discretizó el problema en una malla de 36000 nodos y se implementó en
Dev c++/Open_GL una animación que permitió caracterizar el proceso iterativo.

Ilustración 14. Perfil de presión parcial de oxígeno.

Ilustración 15. Perfil de presión parcial de metanol.

39
Ilustración 16. Perfil de presión parcial de formaldehido.

Ilustración 17. Perfil de presión parcial de agua.

40
 Ilustración 18. Perfil de presión parcial de monóxido de carbono.

Ilustración 19. Perfil de temperatura.

En resumen, una vez seleccionada la expresión que presento mayor estabilidad, se utilizó para
modelar la oxidación parcial de metanol a formaldehido en un reactor tubular simulado en 2D. Para
ello se emplearon diferencias finitas de Crank Nicholson, un algoritmo de relajación y un método
simultaneo de resolución de sistemas de ecuaciones no lineales. En los dos programas
implementados, fue necesario inicializar los nodos cercanos a una solución factible, para garantizar
la convergencia de los métodos, puesto que ambos métodos numéricos, en especial el método de
relajación, son fuertemente dependientes del punto inicial [16]. Por lo tanto se programó una
subrutina que inicializara todos los nodos con el valor obtenido a partir de una solución
unidimensional, y de esta forma asegurar la convergencia de las simulaciones.

41
Optimización económica
La optimización en el diseño y operación de reactores químicos se enfoca en la formulación de una
función objetivo y una descripción matemática del reactor [8]. Para el caso de un reactor tubular
puede describirse el sistema siguiendo alguno de los modelos explicados en secciones previas. Sin
embargo la selección de dicho modelo puede complicar o no el tratamiento numérico del problema.
Siguiendo un modelo unidimensional, la función objetivo está sujeta a un sistema de ecuaciones
algebraico/diferencial. Este tipo de problemas puede optimizarse siguiendo métodos secuenciales
que resuelven iterativamente el ODEs embebido. Para este tipo de problemas también aplican
estrategias simultáneas, que radican en discretizar el ODEs y resolver con programación no lineal el
sistema algebraico resultante [18]. En cuanto a modelos bidimensionales la estrategia simultánea es
frecuentemente implementada [24], con el inconveniente de que involucra un número elevado de
variables en el conjunto de restricciones.

En la optimización de reactores químicos existen alrededor de 13 funciones objetivo típicas [8]. En

este trabajo se estudiaron 4 posibles funciones objetivo, la conversión, la selectividad, el
rendimiento y una función económica. Inicialmente se realizaron barridos, cambiando algunas de
las variables de interés (temperatura de alimentación, presión de metanol y temperatura de
refrigerante) y se observó el comportamiento de cada una de las funciones al cambiar dichas
variables. Finalmente se decidió estudiar el problema con una función económica presentada en 65

(62)

(63)

(64)

(65)

Una vez definida la función objetivo se procedió a seleccionar una estrategia de solución. Se
decidió solucionar el problema a partir de un método simultáneo, y entre los algoritmos de
optimización para la resolución de problemas no lineales se contempló la implementación de
métodos de punto interior (barrera) y el método de programación cuadrática sucesiva (SQP). Dicha
selección se basó principalmente en la disponibilidad de software que permitiera abarcar problemas
con una cantidad considerable de restricciones. Por lo tanto, en lo referente a los métodos del punto
interior se estudió la posibilidad de acoplar el código realizado en el programa Dev c++ con IPOPT,
un software libre recientemente desarrollado en la universidad de Carnegie Mellon [11], que
presenta la ventaja de que permite trabajar problemas de optimización no convexos con un gran
número de variables y adicionalmente ya se ha usado en problemas semejantes a este [24]. En
cuanto al método SQP, Matlab, a partir de la versión 2010a, tiene la opción de resolver problemas

42
no lineales con este algoritmo, y también presenta la posibilidad de paralelizar el código, lo cual
sería conveniente debido a la cantidad de variables involucradas en el problema.

Ilustración 20. Perfil óptimo de Metanol Ilustración 21. Perfil óptimo de Formaldehído.

Ilustración 22. Perfil óptimo de Agua Ilustración 23. Perfil óptimo de monóxido de
carbono

Ilustración 24. Perfil óptimo de temperatura

43
A partir de la herramienta fmincon, disponible en MATLAB, se realizó una primera optimización
paralelizada con 3604 restricciones (6 nodos radiales, 100 nodos axiales, 4 restricciones de
desigualdad, 5 balances de masa y 1 balance de energía), en donde las variables independientes del
problema fueron la longitud y diámetro del reactor, mientras que las condiciones de alimentación al
reactor se mantuvieron fijas. Para este caso se obtuvo que el optimizador se detenía en el momento
que el paso de iteración era del orden de 1x10-6 y las restricciones se cumplían. Sin embargo la
condición de primer orden de optimalidad alcanzaba un valor del orden de 1x10-3. Razón por la
cual se inicializó el problema desde 80 dimensiones aleatorias distintas con el propósito de
encontrar un máximo global en un rango de diámetros y longitudes fijo. Para este estudio la
función objetivo (19) paso de un valor de -0.0367 a un valor de 1.0974 y las variables
independientes tomaron un valor de 8.835 cm de diámetro y 14.08 cm de longitud. Adicionalmente
se calculó el valor de la conversión y rendimiento en el óptimo encontrado y los valores pasaron
respectivamente de 1 y 0.015 a 0.997 y 0.4846.

El perfil de temperatura obtenido al optimizar la función económica con las dimensiones del reactor
se presenta en las Ilustraciones 20-24. En este perfil se observa que la optimización resultó en
disminuir la longitud del reactor de un valor de 70 hasta 14.08 cm y aumentó el diámetro de 1.5 a
8.835 cm, lo cual se atribuye a que las condiciones de alimentación están fijas y para evitar el
consumo de formaldehido, debido a la segunda reacción, se requiere disminuir la longitud del
reactor. Adicionalmente las variables independientes no toman ninguno de los valores extremos
fijos de la simulación por lo que el máximo obtenido no se encuentra en las fronteras del problema.

Ilustración 25. Óptimos locales.

Siendo la temperatura una variable relevante en el proceso, se decidió realizar una optimización con
la temperatura de alimentación y la temperatura de refrigerante como variables independientes, ya
que la temperatura de operación del proceso afecta fuertemente la selectividad, favoreciendo la
reacción principal. Nuevamente se inicializó el problema en 150 temperaturas aleatorias y se
optimizó. Los resultados obtenidos de las 150 iteraciones se presentan en la Ilustración 25.

44
Se observa que para las 150 inicializaciones realizadas, el problema de optimización converge en 18
posibles máximos locales, probablemente debido a la no convexidad del problema, y que para una
temperatura de alimentación de 503.5 K y una temperatura de refrigerante de 515.7 K se presenta
un máximo global en el intervalo de temperaturas permitido. Con el objetivo de verificar los
resultados obtenidos, se simuló el sistema con condiciones por encima y por debajo de las obtenidas
y la función objetivo en ambos casos tomo valores inferiores en comparación con el optimizado de
0.001709. Los perfiles de temperatura y presión parcial obtenidos se presentan en las Ilustraciones
26,27.

Ilustración 26. Perfil de formaldehído óptimo

Ilustración 27. Perfil de Temperatura óptimo

45
Conclusiones
Se realizó una optimización económica para la producción de formaldehído en un reactor tubular,
simulado con un modelo pseudo-homogéneo bidimensional y la cinética de Tesser et al. Para ello se
siguieron dos procedimientos diferentes; inicialmente se optimizo una función económica variando
las dimensiones del reactor y posteriormente la optimización consistió en variar las temperaturas de
alimentación y refrigerante para maximizar la misma función económica. En el primer
procedimiento se encuentra un aumento significativo en la función objetivo al diseñar el reactor con
un diámetro de 8.8 cm y 14.07 cm de longitud. Dichas dimensiones favorecen la producción de
formaldehido sobre agua y monóxido de carbono, por lo cual se logra aumentar la selectividad
hacia la primera reacción del mecanismo reactivo, y minimizan la formación de productos alternos.
En una segunda optimización aumento en menor proporción la función objetivo, mas sin embargo
nuevamente se consiguió incrementar la selectividad del proceso. En este caso se ha encontrado que
manteniendo fijos los parámetros de diseño en 70 cm de longitud y 1.5 cm de diámetro se debe
alimentar la mezcla reactiva a 503.5 K y la temperatura de refrigerante debe mantenerse en 515.7 K.
Es importante resaltar que en ambos procedimientos la optimización numérica se basó en el método
programación cuadrática sucesiva, debido a la no linealidad del sistema algebraico que involucra las
restricciones del problema. Sin embargo fue necesario trabajar el problema en paralelo, puesto que
éste abarca un número considerable de restricciones. Alrededor de 9000 restricciones de igualdad
fueron establecidas y entre 3 y 4 restricciones de desigualdad. Finalmente las simulaciones
llevadas a cabo arrojaron resultados específicos con respecto a las condiciones de operación y
diseño bajo las cuales se optimiza un reactor tubular para la producción de metanol a formaldehído.

46
Bibliografía

[1]A.Constantinides. Numerical methods for chemical engineers with MATLAB applications.

Prentice Hall, 3rd edition, 1998.

[2]A.Quarteroni and A.Valli. Numerical Aproximation of Partial Di erential Equations. Springer,

1997.

[3]B.A.Finlanson. Non linear Analysis in Chemical Engineering. McGraw-Hill, 1980.

[4]D.A.Robb and P.Harriott. The kinetics of methanol oxidation on a supported silver catalyst.
Journal of catalysis, 35, pag 176-183, 1974.

[5]D.Rosenberg. Methods for the Numerical Solution of Partial Di erential Equations. American
Elsevier Publishing Company, 1969.

[6]E.Santacesaria and M.Morbidelli. Kinetics of the catalytic oxidation of methanol to

formaldehyde. Chemical Engineering Science, 36, pag 909-918, 1981.

[7]M.Di Serio E.Tesser and E.Santacesaria. Catalytic oxidation of methanol to formaldehyde: an

example of kinetics with transport phenomena in a packed-bed reactor. Catalysis Today, 77, pag
325-333, 2003.

[8]E.T.Himmelblau. Optimization of chemical processes. McGraw-Hill, 2nd edition, 2001.

[9]F.Frusteri F.Arena and A.Parmaliana. Kinetics of the partial oxidation of methane to

formaldehyde on silica catalyst. AIChE Journal, 46,pag 2285-2294, 2000.

[10]G.F.Froment and K.B.Bischo . Chemical Reactor Analysis. John Wiley and Sons, 2 edition,
1990.

[11]Short Tutorial: Getting Started With Ipopt in 90 Minutes. Andreas w achter. Technical report,
Carnegie Mellon, 2009.

[12]J.H.Mathews. Numerical methods:using matlab. Prentice Hall, 3rd edition, 1999.

[13]J.Schotborgh. An alisis del reactor multitubular para la producci on de formaldeh do por medio
de modelos unidimensionales. Proyecto de grado. Universidad de los Andes.

[14]J.Smith. Chemical Engineering Kinetics. McGrawHill, 3rd edition, 1999.

[15]J.Villadsen and M.L.Michelsen. Solution of Di erential Equation Models by Polynomial

Aproximations. Prentice Hall, 1978.

[16]K.J.Beers. Numerical Methods for Chemical Engineering. Cambrige University Press, 2007.

47
[17]L.C.Windes. Steady state and dynamic modeling of a packed bed reactor for the partial
oxidation of methanol to formaldehyde. Chem.Eng.Comm, 78, pag 143, 1989.

[18]L.T.Biegler. Optimization strategies for complex process models. Advances in chemical

engineering, 18, 1992.

[19]M.B.Cutlip and M.Shacham. Problem solving chemical engineering with numerical methods.
Prentice Hall, 1st edition, 1998.

[20]M.Di Serio P.DOnofrio M.Cozzolino, R.Tesser and E.Santaseraria. Kinetics of the oxidative
dehydrogenation of metanol to formaldehyde by supported vanadium based nanocatalyst. Catalyst
Today, 128: 191-200, 2007.

[21]S.Elnashaie and S.Elshishini. Modelling Simulations and optimizations of industrial xed bed
catalytic reactors. Gordon and Breach Sci, 1994.

[22]T.M.Moustafa. Simulation of the industrial packed bed catalytic reactor for the partial oxidation
of methanol to formaldehyde. Dev.Chem.Eng.Mineral, 11 (3/4).

[23]M.Abou-Elreesh T.M.Moustafa and S.K.Fateen. Modeling, simulation, and optimization of the

catalytic reactor for methanol oxidative dehydrogenation. Excerpt from the Proceedings of the
COMSOL Conference, 2007.

[24]V.M.Zavala and L.T.Biegler. Large scale non-linear programing for the operation of ldpe
tubular
reactors. journal 18th European Symposium on Computer Aided Process Engineering, 2008.

48