La selectividad y el rendimiento de un catalizador, y algunas ideas generales y conceptos importantes para la operación de reactores continuos tanque-agitados (CSTR)

El reactor CSTR

El reactor continuo tanque-agitado (CSTR) es el caballo de batalla de la Ingeniería Química, noble profesión que quedaría reducida a casi nada sin él.

El CSTR y otros reactores se tratan extensivamente durante el estudio de la Ingeniería Química, pero siempre quedan “Capítulos que se le olvidaron a Cervantes”. Por esto en éste artículo expondré algunos conceptos útiles para la operación de los CSTR, y algo acerca de selectividad, y rendimiento de los catalizadores, mediante dos ejemplos.

Algunos conceptos

Concepto de tiempo medio de residencia

El tiempo medio de residencia (TMR) es, estadísticamente, el Primer Momento de la Curva de Tiempos de Residencia (CDTR) de la mezcla reaccionante que pasa por el CSTR. Corresponde a la media aritmética de los tiempos de residencia (Jhonson, Norman and Leone, Fred), (Vogler), que es el indicador de tendencia del valor más frecuente, o más probable de ocurrir.

El tiempo de permanencia de cada molécula que pasa por un CSTR es diferente al de las demás. La CDTR es una muestra representativa de la distribución estadística de los tiempos individuales de residencia del universo de todas las moléculas que pasan por el CSTR. Ésta CDTR muestral tiene una media,  y un TMR, también muestral, que pertenece a un universo de medias.

Importancia del TMR 

El TMR relaciona el tamaño del CSTR con el valor del flujo volumétrico de alimentación. Esta relación se colige de la inspección de la ecuación que lo define, que se muestra a continuación (Vogler).

Producción del CSTR

La producción (throughput) P del CSTR, y el TMR, se relacionan a través del flujo volumétrico de alimentación al reactor. La ecuación para calcular la producción, que también es función de la concentración del producto en él, se muestra a continuación. 

Si se despeja el flujo volumétrico de alimentación FAo de la ecuación que define el TMR y se reemplaza el resultado en la ecuación que define la producción, se obtiene una ecuación que evidencia que la producción es directamente proporcional al volumen y a la concentración de producto, e inversamente proporcional al TMR, que se muestra a continuación.

En razón de lo importante que el TMR es, es necesario que los Ingenieros Químicos tengamos una clara idea de su significado y uso. Por eso he creado un ejemplo sencillo que muestra cómo éste parámetro influye en la producción y en la conversión del CSTR; y cómo su variación determina la velocidad de alimentación del reactor.

Para estos fines he propuesto una reacción de constante cinética k=1×10^-2 s^-1, energía de activación de 5×10^3 J-kmol^-1, que se realiza a 353.15 K en un CSTR. El factor pre-exponencial ko es 4.991493001 x 10-2 s^-1. 

La ecuación fundamental del modelo se muestra a continuación. 

De acuerdo a (Luyben), la ecuación de balance de masa de un CSTR isotérmico de volumen de reacción constante es la siguiente,

Como en la ecuación precedente se reconoce que el factor del primer término es igual a la inversa de Tau, puede escribirse la ecuación en función de éste parámetro, como sigue.

Si se ejecuta la multiplicación del paréntesis de la anterior ecuación puede escribirse de la siguiente forma,

Si se factoriza CA entre el segundo y tercer términos, la ecuación puede compactarse de la forma que a continuación se muestra.

Si en la ecuación anterior se define 

Y

Es posible escribir la separación de las variables en función de éstos parámetros, como sigue:

Integrando se obtiene la siguientes expresión.

Aplicando los límites de integración a la expresión anterior se obtiene la siguiente ecuación.

La anterior ecuación obtiene la expresión siguiente, mediante sencilla manipulación.

Si la última expresión se traslada a Excel™, se pueden obtener los gráficos que a continuación se presentan, donde tan sólo cambiando el valor de FAo se ha cambiado el valor de Tau, posibilitando mostrar que la conversión crece y la producción disminuye con el incremento Tau.

Primer experimento virtual: Investigación sobre el efecto de la variación del tiempo medio de residencia sobre la conversión del reaccionante A, y sobre la producción del reaccionante B

Para observar el efecto del valor del Tau sobre la conversión del reactivo A en un reactor de 1m^3 de volumen, se varió el valor de Tau desde 250s, hasta 1000s mediante la variación del flujo volumétrico de alimentación FAo, para el mismo intervalo de tiempo en todos los casos. 

Los resultados se muestran a continuación: para cada valor particular de Tau, en forma gráfica primero, y en forma tabular consolidada después.

Ilustración 1 Experimento virtual 1 Tau=1000s
Ilustración 2 Experimento virtual 1 Tau=500s
Ilustración 3 Experimento virtual 1, Tau=250

Tabla 1 Variación del Tau, la conversión XA, la producción P, y las concentraciones CA y CB con el flujo de alimentación al reactor FAo

Tabla 1 Variación de Tau, la conversión XA, la producción P, y las concentraciones CA y CB con el flujo de alimentación FAo
Ilustración 4 Variación de la conversión XA y la Producción P con el valor de Tau

Comentarios a los resultados del primer experimento virtual

En los gráficos de la Ilustración 1,Ilustración 2, e Ilustración 3 la concentración final del reaccionante A aumenta con la disminución de Tau, lo que permite concluir que, en general, (también se ilustrará esto para reacciones sucesivas que producen productos no deseados) para aumentar la conversión de A (disminuir la concentración final de A) se debe aumentar el valor de Tau.

El incremento de la conversión con el incremento del valor de Tau implica aumento de producto solamente cuando la reacción no produce productos indeseados. Lo que pasa con el aumento de Tau en el caso de dos reacciones sucesivas que producen producto, y producto indeseado se verá en un ejemplo propuesto por (Rusell T. W. Fraser, Denn, Morton M.) que he desarrollado más adelante. 

Con respecto a las reacciones sucesivas con producción de producto indeseado, se puede adelantar, sin embargo, que cuando éste es el caso el aumento de Tau indudablemente que aumenta la conversión, pero no aumenta la producción de producto porque en ese caso la conversión es igual a la suma de los rendimientos del producto deseado, y del producto indeseado.

Por esta razón, cuando un producto deseado compite con otro que no lo es, el parámetro de interés es el rendimiento, y no la conversión.

La selectividad y el rendimiento de los catalizadores

Voy a utilizar los problemas 12.16-12.20 propuestos por (Rusell T. W. Fraser, Denn, Morton M.) para construir esta parte del artículo, porque creo que constituyen herramientas valiosísimas que ilustran tanto la aplicación de la cinética química a un caso práctico, como el uso de la selectividad y del rendimiento de un catalizador y su aplicación a la operación de un CSTR.

Los autores citados usan la reacción de la Naftoquinona y su conversión a Anhídrido Ftálico, y a “otros” productos de combustión”, reacción que -a su vez- fue estudiada por (DeMaria, Longfield, Butler) en un lecho fluidizado catalítico. En el problema bajo consideración, se pregunta al estudiante cuál de dos catalizadores utilizar entre un catalizador A y un catalizador B, cada uno de diferentes características termodinámicas.

Indican los autores también, que “una sorprendente manera de modelar un lecho fluidizado” es considerarlo como un CSTR, lo que, evidentemente no es de sorprenderse, porque el supuesto fundamental del modelo de CSTR presupone “agitación perfecta”, condición que un lecho fluidizado cumple con largueza. 

Las reacciones sucesivas que se consideran son las siguientes, donde A representa la Naftoquinona, R representa al Anhídrido Ftálico, y S representa a los “productos no deseados de combustión” (CO2 y H2O).

Por los “productos no deseado” del problema, que son productos de combustión, no es difícil inferir que en el lecho fluidizado la Naftoquinona se somete a una pirólisis en presencia de un catalizador, y de O2, pirólisis que obtiene, en una primera reacción, Anhídrido Ftálico; y en una segunda reacción, consecutiva a la primera, CO2 y H2O como productos de combustión.

Por no disponer de software para dibujar estructuras químicas, a continuación muestro las reacciones utilizando fórmulas compactas, con el provisto de que es fácil ver las fórmulas estructurales de los compuestos en Wikipedia™.

Para referencia, la reacción es la siguiente:

El oxígeno está presente en exceso. Por eso la primera reacción es de primer orden sólo respecto de la Naftoquinona.

Los parámetros termodinámicos de los catalizadores que se ofertan al estudiante son los que me permito transcribir a continuación, en unidades SI.

Tabla2 Parámetros termodinámicos de los catalizadores

Como en éste caso los resultados que interesan son la selectividad y el rendimiento que resultarían del uso de cada uno de los catalizadores, he pensado que, antes que nada, es pertinente ahondar un poco en estos conceptos, y compararlos con el concepto conocido de conversión, lo que hago a renglón seguido. 

El paradigma “conversión” versus paradigma “rendimiento”

Los problemas citados son muy útiles para abordar, y aclarar, paradigmas que, como es el caso de la conversión, por haberse esgrimido continuamente, sin pensar ni ahondar mucho, a veces terminan por perder su significado. 

El paradigma que requiere aclaración es “conversión”, que la costumbre ha impuesto como único, pero que ni es único, ni es idéntico al paradigma “rendimiento”.

Veamos esto en mayor detalle.

Para una reacción en que puedan encontrarse presente tres (o más) compuestos, como en el caso del problema propuesto, la suma de las concentraciones de las especies será igual a la concentración original de A, a cualquier tiempo. Si se acepta lo anterior es posible escribir la siguiente igualdad.

Dividiendo la anterior ecuación para  CAo se obtiene la siguiente expresión.

Pasado el primer término del miembro derecho al miembro izquierdo obtiene la siguiente igualdad.

De ésta expresión puede apreciarse claramente que la conversión de A es igual a la suma de los rendimientos de R y de S, y que -si se quiere hilar un poco más fino- el rendimiento de R respecto de A, que es la magnitud de interés, es igual a la conversión de A menos el rendimiento de S respecto de A, como se puede apreciar de la ecuación siguiente.

La anterior igualdad se evidencia gráficamente en la divergencia entre conversión y rendimiento respecto de A en la Ilustración 7 y en la Ilustración 8.

Además, de las ecuaciones precedentes se puede apreciar que: (1) Conversión se refiere a [la fracción de] A transformada, y que: (2) Rendimiento se refiere a [la fracción de] R formada; y que, sin embargo de que las dos magnitudes puedan ser numéricamente iguales en el caso de producción de una sola especie química, no lo son cuando se produce más de una, y que aún siendo iguales tienen diferentes significados.

Los anteriores ecuaciones también muestran que cuando se produce más de una especie química “conversión” se refiere a toda la masa de A convertida, tanto en producto deseado, como en producto indeseado, y que fijarse en su valor para establecer la productividad de la reacción es engañoso; y que el indicador idóneo para esto es el rendimiento de cada especie que se produzca.   

El paradigma “selectividad”

La selectividad, que varia con la temperatura, es la capacidad del catalizador de privilegiar la producción del producto deseado versus la producción del producto indeseado. Una selectividad (teórica) del 100% significaría que sólo se produce el producto deseado. 

Para minimizar la producción de el (o los) productos indeseados se requiere un catalizador de alta selectividad respecto del producto deseado, que en este caso sería el catalizador B (véase la Ilustración 6) para el rango de temperaturas entre 500 K y 560 K, lo que -sin embargo- se contrapone con el bajo rendimiento del catalizador en el rango de temperaturas de alta selectividad. 

Ilustración 5 Rendimiento y selectividad del catalizador A
Ilustración 6 Rendimiento y selectividad del catalizador B

Nótese que, en la ilustración referida, la selectividad del catalizador B es muy cercana a 100% hasta una temperatura de 557 K, lo que quiere decir que en éste rango de temperaturas virtualmente no se produce S, que es lo deseable, pero que lo indeseable es que en este rango el rendimiento no llega ni a 40%, lo que no hace atractiva la posibilidad de usarlo.

¿Cuál es el significado de una selectividad cercana a 100%?

Veamos cuál es el valor de s cuando no se producen productos secundarios.

La definición de selectividad respecto de CR es la siguiente (Rusell T. W. Fraser, Denn, Morton M.)

Si recordamos que a cualquier tiempo mayor que cero,

Y si se reemplaza la última ecuación en la penúltima, es evidente que la selectividad no es otra cosa que la fracción de masa entre el producto CR y la suma de productos CR+CS.

Y que cuando CS=0, s=1 (100%)

Con esto es claro que se podría operar el CSTR con el catalizador B, a 557.15 K (ver Tabla 3), a una selectividad de 99.17% y a un rendimiento de sólo 34.36%.

Tabla 3 Selectividad y rendimiento del catalizador B a diferentes temperaturas, a Tau=7.07s

Una manera realista de evaluar ésta posibilidad, es reparar en que, sin embargo de que se obtendría una fracción CR/(CR+CS) igual a 0.9917, sólo se obtendría un rendimiento de producto R de 34.36%, en mezcla con un 65.64% de reaccionante A, lo que no es muy atractivo 

Una manera alternativa de operar el CSTR sería utilizar el catalizador A, a 617.15 K, para obtener una selectividad 88.86%, y un rendimiento de producto de 82.00% en mezcla con sólo 8% de A, lo que constituye una alternativa más atractiva que la anterior por lo que, a primera vista, usar el catalizador A sería más ventajoso que utilizar el catalizador B (Ver Tabla 4).

Tabla 4 Rendimiento y selectividad del catalizador A a Tau=7.07s

Si se resuelven las ecuaciones siguientes es posible comprobar lo afirmado respecto del rendimiento y selectividad del catalizador A cuando de alimenta el CSTR con una concentración CAo igual a 0.5 mol-m^-3.


Qué efecto produce aumentar el Tau en el caso de reacciones consecutivas que producen producto, y producto indeseado?

Como el tema es importante decidí abordarlo con un ejemplo , porque es necesario que éste efecto quede absolutamente claro.

Para estos fines se procedió a simular la operación de un CSTR a un Tau de 7.07 s, y a un Tau de 10.6 s.

Los resultados de las simulaciones se muestran a continuación.

Tabla 5 Catalizador A rendimiento y selectividad a diferentes temperaturas a Tau=7.07ss
Tabla 6 Catalizador A rendimiento y selectividad a diferentes temperaturas, Tau=10.6s
Ilustración 7 Operación d un CSTR a Tau=7.07s Obsérvese la divergencia entre la conversión y el rendimiento con el aumento de temperatura
Ilustración 8 Operación de un CSTR A Tau=10.6s . Obsérvese que la diferencia entre conversión y rendimiento es mayor cuanto mayor es Tau
Ilustración 9 Comparación entre conversiones en un CSTR a Tau0/.07 y a Tau=10.6s. Obsérvese que la conversión aumenta con el aumento de Tau
Ilustración 10 Rendimientos en un CSTR a diferentes Taus. Obsérvese que a mayor Tau mayor rendimiento hasta 610K, y que el gráfico guarda similitud con el gráfico de conversiones

Conclusiones generales a los resultados anteriores

Los resultados son específicos para el catalizador A del ejemplo. Para generalizarlos, en principio habría que hacer un plan experimental, y llevar a cabo el número de simulaciones (o experimentos) que sea del caso. Sin embargo, el que (Luyben) presente una curva similar a las de las de la Ilustración 7 e Ilustración 8 para el caso de reacciones simultáneas en la que la divergencia en la parte alta del rango de temperaturas es mucho más pronunciada que en los casos aquí referidos me lleva a pensar que el asunto es generalizable.

  1. Los resultados tabulares muestran que los valores de k1 y k2 son los mismos para los diferentes valores de Tau, lo que es correcto porque el valor de k es función exclusiva de la temperatura (véanse Tabla 5 y Tabla 6)
  2. Para rangos bajos de temperatura los valores de rendimiento y conversión son muy similares. La divergencia entre éstos valores se produce en el rango alto de temperaturas de las curvas. 
  3.  Las conversiones aumentan con el valor de Tau, al igual que en el caso de reacciones que sólo producen producto deseado, con las salvedades que quedan señaladas.
  4. El rendimiento es mayor cuanto mayor es el valor de Tau, exceptuando el cuartil superior de la curva. 

El algoritmo de simulación

El algoritmo de simulación está constituido por las tres ecuaciones de balance de masa de las especies A, R, y S. Las ecuaciones de energía no aparecen porque el proceso es isotérmico lo que sin embargo de constituir una simplificación notable, limita el ámbito de observación del analista del proceso, que queda ciego al proceso de arranque del CSTR

Para formular las ecuaciones se debe tener presente que (1) R y se forma a expensas de A, y (2) que S se forma a expensas de R porque se trata de reacciones consecutivas.

Referencias a las velocidades de desaparición y formación, en negritas.

Teniendo en cuenta los dos párrafos anteriores, es muy simple formular las ecuaciones de la manera que a continuación se muestra (Luyben).

Base del análisis: 1 m3

k1CA es la velocidad con que A desaparece y FAoCA/V es la velocidad con que A sale del CSTR.

La ecuación de balance de masa de producto R toma en cuenta la formación de R (igual a la velocidad de desaparición de A, pero con signo cambiado), y la velocidad a la que R desaparece para formar S.

Y, finalmente, la ecuación de balance de masa de la especie S (en realidad las especies CO2 y H2O) es la que se muestra a renglón seguido en donde la velocidad de formación de S es igual a la velocidad de desaparición de R, con signo cambiado, y FAoCR/V es la velocidad a la que S sale del reactor.

La tarea de integrar las ecuaciones, de escribir título y resultados, etc. se puede deducir fácilmente del programa fuente VBA que se incluye.

Programa de simulación

Option Explicit

Dim V, FAo, k1, k2, CRo, CAo, FA, FEV, T, Tmax As Variant

Dim Delltat, Printerval, DeltaCA, DeltaCR, Deltat, DeltaPrint, CR, CA As Variant

Dim Linea_Escritura_Titulos, Linea_Escritura_Resultados, CS, DeltaCS, TR, Tau As Variant

Dim E1, E2, ko1, ko2, R, Principio, u, D, Pi, Aspect_Ratio, A As Variant

Dim Rango, Catalizador  As String

Sub Rusell_Denn()

Principio = Timer

V = Worksheets(“Sheet1”).Cells(3, 3)

FAo = Worksheets(“Sheet1”).Cells(4, 3)

‘k1 = Worksheets(“Sheet1”).Cells(5, 3)

‘k2 = Worksheets(“Sheet1”).Cells(6, 3)

CRo = Worksheets(“Sheet1”).Cells(7, 3)

CAo = Worksheets(“Sheet1”).Cells(8, 3)

FEV = Worksheets(“Sheet1”).Cells(9, 3) ‘ Factor de expnasiónde volumen

Deltat = Worksheets(“Sheet1”).Cells(10, 3)

Tmax = Worksheets(“Sheet1”).Cells(11, 3)

Printerval = Worksheets(“Sheet1”).Cells(12, 3)

DeltaPrint = Worksheets(“Sheet1”).Cells(13, 3)

Catalizador = Worksheets(“Sheet1”).Cells(14, 3)

TR = Worksheets(“Sheet1”).Cells(15, 3)

‘E1 = Worksheets(“Sheet1”).Cells(16, 3)

‘E2 = Worksheets(“Sheet1”).Cells(17, 3)

‘ko1 = Worksheets(“Sheet1”).Cells(18, 3)

‘ko2 = Worksheets(“Sheet1”).Cells(19, 3)

R = Worksheets(“Sheet1”).Cells(20, 3)

Aspect_Ratio = Worksheets(“Sheet1”).Cells(21, 3)

Pi = 3.14159

Linea_Escritura_Titulos = 100

Linea_Escritura_Resultados = Linea_Escritura_Titulos + 1

Call Calculo_de_Constantes_Cineticas(ko1, ko2, k1, k2, E1, E2, TR, R, Catalizador)

Call Calculo_Velocidad_Alimentacion(V, Aspect_Ratio, u, Tau)

‘Definición de variables

‘FA = FEV * FAo

FA = FAo

CR = CRo: CA = CAo

‘Limpieza de pantalla

Rango = “A100:k5000”

Worksheets(“Sheet1”).Range(Rango).Clear

Call Escritura_Titulos(Linea_Escritura_Titulos)

CS = 0

DeltaCR = 0: DeltaCA = 0

For T = 0 To Tmax Step Deltat

        If T = 0 Then Call Escritura_Resultados(Linea_Escritura_Resultados, CA, CR, CS, T)

        If T >= Printerval Then

            Call Escritura_Resultados(Linea_Escritura_Resultados, CA, CR, CS, T)

            Printerval = Printerval + DeltaPrint

        End If

    DeltaCA = ((1 / V) * (FAo * CAo – FA * CA) – (k1 * CA)) * Deltat

    CA = CA + DeltaCA

    If CRo <> 0 Then

    DeltaCR = ((1 / V) * (FAo * CRo – FA * CR) – (k1 * CA) – k2 * CR) * Deltat

    Else

    DeltaCR = ((-FAo / V) * CR + k1 * CA – k2 * CR) * Deltat

    End If

    CR = CR + DeltaCR

    DeltaCS = ((-FAo / V) * CS + k2 * CR) * Deltat

    CS = CS + DeltaCS

Next T

Call Escritura_Resultados(Linea_Escritura_Resultados, CA, CR, CS, T)

Call Escritura_Parametros_Varios(Linea_Escritura_Resultados, CA, CR, CAo, Catalizador, V, FAo, Principio, k1, k2, Tau, u, TR)

End Sub

Sub Escritura_Titulos(Linea_Escritura_Titulos)

Worksheets(“Sheet1”).Cells(Linea_Escritura_Titulos, 1) = “Tiempo,s”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Titulos, 2) = “CA, kmol-m^-3”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Titulos, 3) = “CR, kmol-m^-3”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Titulos, 4) = “CS, kmol-m^-3”

End Sub

Sub Escritura_Resultados(Linea_Escritura_Resultados, CA, CR, CS, T)

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados, 1) = T

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados, 2) = CA

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados, 3) = CR

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados, 4) = CS

Linea_Escritura_Resultados = Linea_Escritura_Resultados + 1

End Sub

Sub Calculo_de_Constantes_Cineticas(ko1, ko2, k1, k2, E1, E2, TR, R, Catalizador)

If Catalizador = “A” Then

    E1 = 175680000#: E2 = 84490000#

    ko1 = 1.25E+15: ko2 = 251000#

ElseIf Catalizador = “B” Then

    E1 = 86240000#: E2 = 193400000#

    ko1 = 9125386#: ko2 = 1.5984470922E+15

End If

k1 = ko1 * Exp(-E1 / (R * TR))

k2 = ko2 * Exp(-E2 / (R * TR))

End Sub

Sub Calculo_Velocidad_Alimentacion(V, Aspect_Ratio, u, Tau)

D = ((4 * V) / (Pi * Aspect_Ratio)) ^ 1 / 3

A = (Pi * (D) ^ 2) / 4

u = FAo / A

Tau = V / FAo

End Sub

Sub Escritura_Parametros_Varios(Linea_Escritura_Resultados, CA, CR, CAo, Catalizador, V, FAo, Principio, k1, k2, Tau, u, TR)

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados, 1) = “Temperattura del Reactor,  TR= ” & Round((TR), 2) & “K”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados + 1, 1) = “Rendimiento (CR/Cao)*100 del catalizador ” & Catalizador & “= ” & Round((CR / CAo) * 100, 2) & “%”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados + 2, 1) = “Selectividad del catalizador, (CR / (CAo – CA))*100 ” & Catalizador & “= ” & Round((CR / (CAo – CA)) * 100, 2) & “%”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados + 3, 1) = “Tiempo medio de residencia, (V/FAo)=  ” & Round(Tau, 2) & “s”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados + 4, 1) = “Velocidad del flujo de alimentación,( FAo//A)=  ” & Round(u, 2) & ” m/s”

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados + 5, 1) = “Constantes cinéticas=  ” & “k1= ” & Format(k1, “scientific”) & “, k2= ” & Format(k2, “scientific”)

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados + 6, 1) = “XA=  ” & Round((1 – (CA / CAo)) * 100, 2)

Worksheets(“Sheet1”).Cells(Linea_Escritura_Resultados + 7, 1) = “Tiempo de ejecución=  ” & Round(((Timer – Principio) / 60), 2) & “min”

End Sub

Colofón

El artículo salió un poco largo. Me disculpo. Espero que les guste y aproveche. Saludos atentos.

Bibliografía

DeMaria, Longfield, Butler. Ind. Eng. Chem. (1961): 53, 259.

Jhonson, Norman and Leone, Fred. Statistics and Experimental Design. New York: John Wiley & Sons, Inc., 1964.

Luyben, William L. Chemical Reaction Engineering. New Jersey: John Wiley & Sons, 2007.

Rusell T. W. Fraser, Denn, Morton M. Introduction to Chemical Engineering Analysis. New York: John Wiley & Sons, Inc., 1972.

Vogler, Scott H. Elements of Chemical REaction Engineering. New Jersey: Prentice Hall, Inc., 1999.

Publicado en Sin categoría

Uso de Ecuaciones de Estado Estacionario y de Estado Transitorio Para el Diseño de un Reactor Tanque-Agitado Continuo: El Qué y el Para Qué

La verdad es que en mi trabajo de modelado matemático siempre les he hecho un poco de asco a las ecuaciones de estado estacionario, porque siempre he considerado que las únicas que contaban eran las conocidas ecuaciones de estado transitorio, que -al fin y al cabo- permiten modelar, y “fotografiar” el comportamiento transitorio de un reactor, y pronosticar una explosión, o una reacción descontrolada, lo que es mucho más descriptivo e inmensamente más útil que lo que se obtiene por métodos estáticos, que son simples valores de volumen de reactor y/o conversión de un reaccionante, que no dicen nada acerca del arranque, ni de la historia de los hechos entre éste y el establecimiento del estado estacionario del reactor.

Pero no es menos verdad que es siempre es posible reconsiderar los puntos de vista, y revaluar criterios. En eso consiste la sabiduría, dicen. Y en esta tarea encontré una especie de eslabón perdido entre las ecuaciones estacionarias y las ecuaciones transitorias, que a renglón seguido quiero compartir con ustedes, amables lectores.

El Problema

El problema, si así se le quiere llamar, es que en un programa se simulación de un reactor tanque-agitado continuo que considera una reacción exotérmica siempre existe la incógnita de cuál será la temperatura de estabilización de la reacción; cuánto, y cómo se habrá que enfriar el reactor, o si éste podrá manejarse de modo adiabático, sin enfriamiento.

Y para que el modelizador pueda contestarse estas preguntas existen dos caminos, a mi entender: El primero es establecer un mecanismo de control Proporcional-Integral (PI), o PID que manipule el flujo de fluido de enfriamiento a la chaqueta del reactor en una modalidad “cascada” (Harriott), que lo enfriará, y ver si es que con éste enfriamiento la reacción arranca, o no; y si esto no sucede, quitar el enfriamiento y ver si la reacción arranca, y se sostiene, con la evolución de calor que produce su exotermicidad.

Es claro que estos caminos no son muy ortodoxos que digamos, y que un poco son de fuerza bruta, pero deberían dar resultado, si es que el modelizador acierta con la zona termodinámica de viabilidad de la reacción, que por fuerza está delimitada por la temperatura de arranque, y por la temperatura de estabilización de la reacción, que puede ser superior, o inferior a la temperatura de arranque.

Lo anterior implica que se puede arrancar desde una temperatura superior a la de arranque, o desde una temperatura inferior a ésta(T. W. Fraser, Morton M. Denn).

A la vista de esta problemática, y como no tengo a quién consultar, el único camino que me quedó fue ponerme a leer (“Agarrá los libros, que no muerden”, decía Sarmiento, el ilustre argentino), y encontré el método de resolver el problema, e incidentalmente, la manera en la que los profesores especifican los parámetros cinéticos de los problemas, lo que paso a comentarles inmediatamente.

La Solución al Problema

Las ecuaciones transitorias de un reactor son las siguientes (W. L. Luyben).

En las ecuaciones anteriores el subíndice R se refiere al reactor, y el subíndice j se refiere a la chaqueta (jacket). Como el calor de reacción  es negativo en el caso de una  reacción exotérmica (convención norteamericana de signo) el término correspondiente se vuelve positivo, como debe ocurrir para el caso de una reacción exotérmica cuando el calor de reacción exotérmico se consigna como negativo.

En el caso de las ecuaciones anteriores el volumen del reactor, y el de la chaqueta, son constantes, lo que corresponde para el caso físico en que tanto el uno como la otra se hayan cargado hasta su nivel de operación antes del arranque.

Respecto de la anterior afirmación me permito agregar que esa es la manera en que se debe arrancar un reactor no adiabático, para evitar que el área disponible para el enfriamiento varíe con el tiempo, porque eso desequilibraría la transferencia de calor al arranque y podría causar un descontrol de la exotermicidad de la reacción (runaway reaction).

Si el reactor no se arrancase lleno, su volumen y el de la chaqueta no serían constantes, deberían permanecer dentro de los cocientes diferenciales, y deberían inicializarse como un producto, que debe tratarse, como si fuese una variable-producto según se indica en la bibliografía (W. L. Luyben), y que deberían actualizarse por división por la variable secundaria.

Lo anterior quiere decir que, si la variable producto es,

su valor actualizado en cada iteración sería 

y que su valor actualizado se calcularía mediante la siguiente división

Volviendo a nuestro problema. (T. W. Fraser, Morton M. Denn) toman las ecuaciones transitorias anteriores y las igualan a cero con lo que las convierten en ecuaciones de estado estacionario, de la siguiente forma.

La segunda cosa que hacen los profesores citados es trabajar sólo con las ecuaciones del reactor, porque la operación del mismo es adiabática, sin enfriamiento ni calentamiento.

Una vez realizada la anterior manipulación de las ecuaciones, (T. W. Fraser, Morton M. Denn)definen los siguientes parámetros (la nomenclatura con subíndices es mía, para homologarla con la de mi programa).

Substituyendo los parámetros anteriormente definidos, las ecuaciones pueden escribirse de la manera siguiente.

En estas ecuaciones, el subíndice “i-1” se refiere a las condiciones corriente arriba del reactor, y el subíndice “i” e refiere a las condiciones dentro del reactor. Insisto que la ecuación correspondiente a la chaqueta no se toma en cuenta por tratarse de una operación adiabática.

Despejando  en ambas ecuaciones se obtienen dos expresiones para la concentración como función de la temperatura del reactor.

Reemplazando el valor de los parámetros y demás variables estáticas en las funciones así definidas y calculando sus valores para diferentes temperaturas es posible obtener la solución gráfica y aritmética del sistema, para un conjunto de valores de diseño y de operación del reactor.

Yo utilicé los valores de (T. W. Fraser, Morton M. Denn) con la excepción del calor de reacción, ya que el suministrado por los estos autores obtenía una solución trivial (hay que enfatizar, y eso lo pueden leer de la publicación de la referencia, que en la mayoría de este tipo de problemas puede haber más de una solución), porque -al fin y al cabo- estamos hablando de raíces de ecuaciones.

A continuación se muestra la matriz Excel de datos para la solución de las ecuaciones.

Ilustración 1. Datos para la resolución de las ecuaciones estáticas
Ilustración 2. Tabla de la resolución de las ecuaciones estáticas

La solución de las ecuaciones se ha resaltado en color rojo.

La “solución gráfica” es la siguiente.

Ilustración 3. Solución “gráfica” de las ecuaciones estáticas

Como puede apreciarse el gráfico como tal no permite apreciar el valor numérico de la solución porque las ecuaciones no son muy bien condicionadas; sin embargo la Tabla de valores de Excel™ sí permite la valoración de la respuesta.

A la fecha de publicación del muy valioso libro de (T. W. Fraser, Morton M. Denn)no existía Excel™, ni computadores personales. Había “mainframes”, se programaba en Fortran IV, y se codificaba en tarjetas perforadas (eran los albores de la programación, pero aún con esas dificultades se modelaba!).

Comprobación de la validez de la solución gráfica

Para los fines de comprobar la validez de la solución obtenida, invalorable para la simulación porque de lo contrario no se sabe a qué temperatura arrancar, ni si la operación es o no adiabática, creé y ejecuté un programa que originalmente tenía por objeto resolver un problema de (W. L. Luyben)de tres reactores en serie, que es muy interesante, del que postergo cualquier comentario para cuando lo resuelva…

A continuación se presenta la Interface de Usuario del programa, que contiene datos que no pertenecen al problema en sí, la tabla de valores que se obtuvo, y la gráfica (dinámica) de la evolución de la temperatura del reactor con el tiempo.

Ilustración 4. Interfase de Usuario del programa dinámico
Ilustración 5. Tabla de resultados del programa dinámico
Ilustración 6. Arranque y evolución de temperaturas del reactor adiabático tanque-agitado continuo

Comentarios a los resultados

Me gustan mucho este tipo de disquisiciones donde se puede apreciar el progreso y sofisticación que las herramientas de simulación han experimentado a través del tiempo, así como la evolución, sindéresis y vigencia del pensamiento genial de maestros como Harriott, Luyben, Fraser y Denn, y cómo los diferentes tópicos que ellos abordaron, se complementan de manera natural, como siempre ocurre con aquello que es verdaderamente valedero.

Respecto de los resultados, pienso que son convincentes, porque son coherentes entre sí, válidos e interesantes en tanto demuestran la complementariedad entre lo estático y lo dinámico del análisis de procesos, aunque -si se piensa bien- lo estático no lo es tanto porque allí también varía el tiempo, porque varía la temperatura, que en un arranque es función de éste.

Respecto de las referencias citadas, debo enfatizar su interés y pertinencia. (T. W. Fraser, Morton M. Denn) ilustran su libro con problemas con datos reales de cinética, lo que no es muy común; (W. L. Luyben) enfoca las bases de la simulación dinámica en el primer libro citado; y el diseño y control de reactores, en la segunda obra abordando aspectos de estado estacionario, y aspectos de estado transitorio, asuntos todos que deberían ser el pan de cada día de todo Ingeniero Químico. Harriott, por su parte, expone los esquemas fundamentales de control de reactores, evaporadores y columnas de destilación. Qué más se puede pedir!

Felicidades. Espero que hayan encontrado interesante la disquisición ingenieril, y el poquito filosófica también.  

Publicado en Sin categoría

Dimensionamiento y Arranque de un Calentador Solar Doméstico

Como se muestra en la Ilustración 1, un calentador solar doméstico de agua está constituido por un colector solar, generalmente constituido por varios paneles conectados en paralelo (no más de seis de acuerdo a las recomendaciones constantes en https:/www.buildingsolar.com/References/ColFlowRate.htm); un tanque de almacenamiento del agua caliente, que sirve como repositorio para el almacenamiento del agua caliente que ha de abastecerse al hogar familiar donde se ha instalado; una bomba de circulación del agua, o de la mezcla anti congelante que circule en circuito del colector, en los países donde las temperaturas ambientes podrían congelar el agua dentro de los tubos del colector.

Hay calentadores solares directosen los que el agua caliente que sale de la batería de colectores simplemente se envía al usuario, y hay calentadores, como el de la siguiente ilustración en el que el calor recolectado por la batería de colectores se transfiere al agua contenida en el tanque de almacenamiento a través de un serpentín helicoidal, a veces denominado intercambiador de calor.

En el caso de los calentadores solares como el de la Ilustración 1el tanque de almacenamiento recibe agua fría del hogar que se calienta en el calentador solar por medio del intercambiador de calor, para enviarse hacia el hogar de acuerdo con la demanda que se vaya presentando.

Para los fines del modelo matemático que se presenta se ha propuesto un calentador eléctrico auxiliar a corriente máxima de 2 amperios resistencia de 2,000 Ohms y cos igual a 0.8 en el que la corriente se modula mediante un controlador PI que actúa solamente cuando existe demanda. 

Ilustración 1. Esquema del Calentador Solar Domestico Típico

Objetivo del modelo

En el caso que nos ocupa el modelo, y el programa, ilustran el arranque y dimensionamiento del área del colector y del tamaño del tanque de almacenamiento en función de la demanda que el Calentador Solar Doméstico ha de abastecer. En el caso presente, la demanda propuesta se muestra en la Ilustración 2. 

La correspondiente integración,  donde ml es la demanda instantánea, y t es tiempo, obtiene 323.1 litros demandados en 24 horas.

Ilustración 2. Esquema de la demanda de agua caliente sobre el calentador solar

El Modelo Matemático

            Balance de energía alrededor del tanque de almacenamiento

La ecuación fundamental del modelo es un balance integro-diferencial de potencia alrededor del tanque de almacenamiento, que se basa en el cálculo de la temperatura instantánea del agua contenida en el tanque, la relaciona con la potencias transferidas por el serpentín al agua del tanque, al agua circulada hacia el domicilio al que sirve el Calentador Solar, con la potencia que se pierde a través del aislamiento térmico del tanque, y con potencia que el calefactor eléctrico auxiliar aporta.

La ecuación es la siguiente.





En la ecuación ds es el diámetro nominal del serpentín; Tses la temperatura del agua en el serpentín; x es la dimensión a lo largo del eje del serpentín; L es el largo del serpentín; mles el flujo másico del agua que se transfiere al domicilio abastecido; Cp es el calor específico del agua; Ttes la temperatura del agua contenida en el tanque; TTies la temperatura inicial del agua en el tanque; Usy Uson los coeficientes globales de transferencia del serpentín y del tanque; ATes el área exterior del aislamiento térmico del tanque; Taes la temperatura ambiente; I es la corriente en amperios que circula por el circuito de calefacción auxiliar; R es la resistencia óhmica del elemento calefactor y  es el factor de potencia, y M es la masa de agua contenida en el tanque de almacenamiento.

Cálculo de la temperatura del agua que deja el colector

La temperatura del agua se obtiene del simple balance de energía entre el agua que sale de una placa colectora conectada en paralelo con otras placas. La siguiente expresión es válida para cualquier instante de tiempo.


En la ecuación anterior, mces el flujo que entra al cabezal (header) de NP placas conectadas en paralelo; Cp es el calor específico del agua; Ac es la suma de las superficies de todas las placas conectadas en paralelo; Ins(horas) es una función de usuario que permite calcular la energía que capta el banco de colectores a cada instante del día..

Si se despeja el valor de Ts de la ecuación anterior se obtiene la expresión para el cálculo de la temperatura de agua que circula en el serpentín, que varía con el tiempo, porque Ts es función de éste. La ecuación para el cálculo es la siguiente.


La eficiencia del colector  se puede calcular, para una placa colectora de reflectancia y la absorbancia conocidas, mediante la expresión siguiente (Salcines et.al.), que en este contexto tiene sólo valor referencial, porque el modelo que aquí se presenta no está referido a ninguna placa colectora en particular. La ecuación es la siguiente. Taes la temperatura ambiente. A y B son parámetros de la ecuación.

Balance de energía alrededor del serpentín

El balance de energía alrededor de un elemento infinitesimal del serpentín sirve para calcular la disminución de la temperatura del agua del serpentín, a lo que ésta circula por éste, que -a su vez- es consecuencia de la transferencia de calor del serpentín hacia el seno del agua contenida en el tanque de almacenamiento. La ecuación es la siguiente.

La expresión para el cálculo de la variación de la temperatura TSes la siguiente.

Irradiación solar y orientación de los paneles solares

De acuerdo los valores medidos por la Estación EMA de la Universidad de San Francisco de Quito, USFQ, a 0o11’47” S,78o26’6″S la radiación solar, en W/m2, en Cumbayá, Quito, D.M. (researchgate.net) puede representarse gráficamente como se muestra a continuación.

Ilustración 3. Radiación Solar en Cumbayá, Quito, D.M. Ecuador

Por otro lado del sitio globalsolaratlas.info el valor de Global Horizontal Irradiance (GHI) para Quito, Ecuador, es igual a 2,000 kW-h/año, valor que asumiendo un año de 365 días y días de brillo solar de 12 horas, obtiene 456.6 W/m2, lo que puede interpretarse como un valor medio anual. 

Por las razones antedichas, se decidió utilizar un valor equivalente al 50% (“50% de nubosidad”) del valor consignado por la USFQ en la ejecución del programa.

Respecto de la orientación azimutal de los paneles solares, de acuerdo a la información que se encuentra en la Internet, si el colector está ubicado en el hemisferio norte, debería orientarse hacia el sur; y si está el hemisferio sur, debería orientársele la norte, y el ángulo de elevación de los paneles solares debería ser igual a la latitud del sitio geográfico donde estos se hallen.

Análisis del efecto de la variación del volumen del tanque de almacenamiento

Como la demanda de agua del domicilio que sirve el calentador solar puede representarse por

y cuando se evalúa dentro del programa es igual a 323.1 litros/24 horas, cabe preguntarse cuál será el efecto de aumentar el tamaño del tanque de almacenamiento de un volumen de 100 L a un volumen de 300 L, lo que puede contestarse fácilmente cambiando el campo respectivo de Vt igual a 100 L a Vt igual a 300 L cuando se ejecuta el programa.

El resultado del análisis de la evolución de las temperaturas para los dos volúmenes de tanques se puede observar en las dos ilustraciones siguientes en donde se han superpuesto tres ejecuciones para Vt igual a 100 L, y otras tres para Vt igual a 300 L, y donde cada trazo representa 24 horas de operación, a Tti igual a la temperatura final del trazo anterior, con lo que la representación abarca 72 horas, a insolación constante para las horas de brillo del sol.

Ilustración 4. Evolución de los valores de temperaturas para el calentador solar equipado con un tanque de 100 L de volumen, período de arranque de tres días consecutivos
Ilustración 5. Evolución de los valores de temperatura para el calentador solar equipado con un tanque de 300 l de volumen, período de arranque de tres días consecutivos

Conclusión al análisis del efecto de la variación del volumen del tanque de almacenamiento

De la observación de las ilustraciones anteriores puede concluirse que aumentando el volumen del tanque de almacenamiento a 300 L el perfil de temperaturas que se obtiene al tercer día de operación es mucho más tendido y uniforme que en el caso del tanque de 100 L, por lo que, aunque la inversión fija del proyecto aumenta un poco, es más conveniente instalar un tanque de 300 L que uno de 100 L, para la demanda especificada en la Ilustración 2, lo que se debe a que la masa térmica del agua contenida en el tanque es mayor cuando la capacidad de almacenamiento del tanque es mayor, lo que hace que el orden de magnitud de las fluctuaciones de temperatura sea menor, confiriéndole mayor estabilidad térmica , como se puede inferir de la ecuación de balance de potencia alrededor del tanque.  

Debe anotarse también que el aumento del volumen del tanque no implica aumento del área del banco de colectores, lo que es ventajoso porque los colectores constituyen el equipo más importante y más caro de la instalación.

El Programa 

El programa se realizó en VBA para Excel™ de Office 635 para Mac.; se ejecuta en alrededor de 4.07 min, lo que no es malo, porque es un poco largo.

En Ingreso de datos y los resultados gráficos y tabulares se despliegan en la misma hoja. A continuación se muestra la matriz de ingreso de datos.

            Ingreso de datos

Ilustración 6. Matriz de Ingreso de Datos

            Resultados

Los resultados gráficos de una ejecución del programa pueden apreciarse a continuación.

Ilustración 7. Resultados gráficos de la ejecución del programa

Código VBA 

Como este blog es fundamentalmente para estudiantes, acostumbro a insertar el código VBA en todos los artículos, para que el lector que así lo desee pueda copiarlo y pegarlo en el editor de VBA de Excel™ para ejecutarlo, si así lo desea.

Este es el código.

Dim Ts, Tc, CpH2O, Ac, I, TiempoComienzo, TiempoTranscurrido, SumCorr, Ilim, Ti, Kc, DeltaTt, Td, Pi, ds, Deltax, Tti, At, Ta, R, fi, Us, L, ml, Theta, FN As Variant

Sub colector_solar()

TiempoComienzo = Timer

Pi = 3.14159: Mt = 0: At = 0

Ro = Worksheets(“Sheet2”).Cells(5, 5)

CpH2O = Worksheets(“Sheet2”).Cells(6, 5)

Ac = Worksheets(“Sheet2”).Cells(7, 5)

‘Dt = Worksheets(“Sheet2”).Cells(10, 5)

Ta = Worksheets(“Sheet2”).Cells(11, 5)

Vt = Worksheets(“Sheet2”).Cells(12, 5)

If Vt = 100 Then

    Ht = 1.1 ‘m

    Dt = 0.37

    Eat = 0.05 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

ElseIf Vt = 150 Then

    Ht = 1.5 ‘m

    Dt = 0.37

    Eat = 0.05 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

ElseIf Vt = 200 Then

    Ht = 1.39 ‘m

    Dt = 0.55

    Eat = 0.05 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

ElseIf Vt = 300 Then

    Ht = 1.75 ‘m

    Dt = 0.56

    Eat = 0.05 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

ElseIf Vt = 400 Then

    Ht = 1.6 ‘m

    Dt = 0.7

    Eat = 0.06 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

ElseIf Vt = 500 Then

    Ht = 1.95 ‘m

    Dt = 0.7

    Eat = 0.06 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

ElseIf Vt = 600 Then

    Ht = 1.98 ‘m

    Dt = 0.75

    Eat = 0.0575 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

 ElseIf Vt = 700 Then

    Ht = 2.24 ‘m

    Dt = 0.75

    Eat = 0.0575 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

ElseIf Vt = 800 Then

    Ht = 1.91 ‘m

    Dt = 1

    Eat = 0.1 ‘Espesor de aislamiento

    Call Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

End If

mc = Worksheets(“Sheet2”).Cells(13, 5)

L = Worksheets(“Sheet2”).Cells(14, 5)

Tti = Worksheets(“Sheet2”).Cells(15, 5)

Tc = Worksheets(“Sheet2”).Cells(16, 5)

ds = Worksheets(“Sheet2”).Cells(17, 5)

Us = Worksheets(“Sheet2”).Cells(18, 5)

‘Ins = Worksheets(“Sheet2”).Cells(19, 5)

‘Ut = Worksheets(“Sheet2”).Cells(20, 5)

ThetaMax = Worksheets(“Sheet2”).Cells(21, 5)

DeltaThetaPrint = Worksheets(“Sheet2”).Cells(22, 5)

Deltax = Worksheets(“Sheet2”).Cells(23, 5)

‘Parametros de eficiencia del colector

A = Worksheets(“Sheet2”).Cells(24, 5)

B = Worksheets(“Sheet2”).Cells(25, 5)

‘Parametros de impresion

Incr_Impr = Worksheets(“Sheet2”).Cells(26, 5)

R = Worksheets(“Sheet2”).Cells(27, 5)

Ilim = Worksheets(“Sheet2”).Cells(28, 5)

fi = Worksheets(“Sheet2”).Cells(29, 5)

BM = Worksheets(“Sheet2”).Cells(30, 5)

Ttset = Worksheets(“Sheet2”).Cells(31, 5)

Kc = Worksheets(“Sheet2”).Cells(32, 5)

Ti = Worksheets(“Sheet2”).Cells(33, 5)

Td = Worksheets(“Sheet2”).Cells(34, 5)

Solucion = Worksheets(“Sheet2”).Cells(35, 5)

FN = Worksheets(“Sheet2”).Cells(36, 5)

Fila_Comienzo_Titulos = 45: Fila_Comienzo_Resultados = Fila_Comienzo_Titulos + 1

‘Parametros generales de control del programa

SumTheta = 0 ‘Incializa totalizador de tiempo, segundos

Tt = Tti

Ts = 18 ‘Se fija Ts en Tti para que no arrastre vaores

Ats = (Pi * (ds) ^ 2 / 4) ‘Area transversal del serpentin

NP = Ac / 2.5

ml = 0: Sum_ml = 0

‘Limpieza de Sheet2

Worksheets(“Sheet2”).Range(“B38:n300”).Clear

Call Escritura_de_Titulos(Fila_Comienzo_Titulos)

If Solucion = “Serpentin” Then

Call Calculo_de_Ts(n, Ac, mc, CpH2O, Tc, Ts, Tt, SumTheta, x, Deltax, NP, Us, Pi, ds, TSR, horas)

End If

Call Calculo_Variacion_Temperatura_Serpentin(Us, Pi, ds, Tt, mc, CpH2O, Ts, Ats, Ro, Deltax, Mt, _ At, Ta, L, SumTheta, DeltaThetaPrint, ThetaMax, Fila_Comienzo_Resultados, n, Ut, _

Incr_Impr, NP, Dt, Eat, A, B, Tc, ml, Iset, R, Ttset, BM, fi, Tti, TSR, horas, Solucion, Ac)

End Sub

Sub Calculo_eficiencia(A, B, Tc, Ta, n, mc, Ts, Ac, CpH2O, horas)

If horas <> 0 Then

n = A + B * ((Tc – Ta) / Ins(horas, FN)) ‘n=Eficiencia

End If

End Sub

Sub Calculo_de_Ts(n, Ac, mc, CpH2O, Tc, Ts, Tt, SumTheta, x, Deltax, NP, Us, Pi, ds, TSR, horas)

    If horas <> 0 And horas <= 17.5 Then

   Ts = Tc + (Ac * (Ins(horas, FN) * n) / ((mc / NP) * CpH2O * 1000)) ‘Temperatura agua

   ‘                                                                            en el banco de colectores

   Else

    Ts = 18

    End If

End Sub

Sub Calculo_Variacion_Temperatura_Serpentin(Us, Pi, ds, Tt, mc, CpH2O, Ts, Ats, Ro, Deltax, Mt, _

At, Ta, L, SumTheta, DeltaThetaPrint, ThetaMax, Fila_Comienzo_Resultados, n, Ut, _

Incr_Impr, NP, Dt, Eat, A, B, Tc, ml, Iset, R, Ttset, BM, fi, Tti, TSR, horas, Solucion, Ac)

DeltaTheta = ((Ats * Ro) / mc) * Deltax  ‘DeltaTheta es intervalo de tiempo transcurrido, en segundos

   Do While SumTheta <= ThetaMax

        For x = 0 To L Step Deltax

                If SumTheta = 0 Then

                        Call Escritura_Resultados(Fila_Comienzo_Resultados, SumTheta, x, Ts, Tt, Incr_Impr, ThetaMax, x, mc, CpH2O, Tc, Ac, NP, ml, I, TSR, horas, Sum_ml)

                End If

                    If horas <= 17.5 Then

                        DeltaTs = -((Us * Pi * ds * Deltax * (Ts – Tt)) / (mc * CpH2O * 1000))

                        Ts = Ts + DeltaTs

                    End If

                horas = 6.4 + SumTheta / 3600

                SumTheta = SumTheta + DeltaTheta ‘SumTheta totaliza el tiempo transcurrido,segundos

                Tmean = (Tt + Ts) / 2

                If Tmean > 0 Then

                    Call Calculo_de_Ut(Dt, Eat, Ut, Tmean)

                Else

                    Stop

                End If

                Call Demanda_Domestica_de_Agua(SumTheta, ml)

                Sum_ml = Sum_ml + ml * DeltaTheta

                    If Solucion = “Serpentin” Then

                    Call Control_PID_de_Corriente(Ttset, Ilim, BM, I, Tt, ml, DeltaTheta, Kc, Ti, DeltaTt, Td)

                    End If

                If SumTheta < 17.56 Then

                Call Calculo_eficiencia(A, B, Tc, Ta, n, mc, Ts, Ac, CpH2O, horas)

                End If

                Call Actualizacion_Temperatura_Tanque(CpH2O, Ts, Us, Pi, Tti, At, I, R, fi, ds, _

                DeltaTheta, Mt, Ta, x, Tt, Ut, mc, ml, Solucion, horas, Ac, n)

                    If SumTheta >= DeltaThetaPrint Then

                            Call Escritura_Resultados(Fila_Comienzo_Resultados, SumTheta, Ut, Ts, Tt, Inc_Impr, ThetaMax, x, mc, CpH2O, Tc, Ac, NP, ml, I, TSR, horas, Sum_ml)

                            DeltaThetaPrint = Incr_Impr + DeltaThetaPrint

                    End If

                        If SumTheta < 43200 Then ‘43200 equivale a las 18.00 horas

                                    If Solucion = “Serpentin” Then

                                    Call Calculo_de_Ts(n, Ac, mc, CpH2O, Tc, Ts, Tt, SumTheta, x, Deltax, _

                                    NP, Us, Pi, ds, TSR, horas)

                                    End If

                        ElseIf Ts <= Tt And Solucion = “Serpentin” Then

                       Exit For

                    Else

                        Ts = 18

                    End If

        Next x

    Loop

End Sub

Sub Actualizacion_Temperatura_Tanque(CpH2O, Ts, Us, Pi, Tti, At, I, R, fi, ds, _

DeltaTheta, Mt, Ta, x, Tt, Ut, mc, ml, Solucion, horas, Ac, n)

If Solucion = “Serpentin” Then

    DeltaTt = (DeltaTheta / (1000 * Mt * CpH2O)) * ((Us * Pi * ds * x * (Ts – Tt)) – _

    (1000 * ml * CpH2O * (Tt – Tti)) – (Ut * At * (Tt – Ta)) + (I * R * Cos(fi)))

ElseIf Solucion = “Clasica” Then

    DeltaTt = (DeltaTheta / (1000 * Mt * CpH2O)) * ((Ins(horas, FN) * Ac * n) – _

    (1000 * ml * CpH2O * (Tt – Tti)) – (Ut * At * (Tt – Ta)) + (I * R * Cos(fi)))

End If

Tt = Tt + DeltaTt

End Sub

Sub Masa_y_Area_Exterior_Tanque(Pi, Dt, Ht, Ro, Eat, Vt, Mt, At)

At = Pi * (Dt + 2 * Eat) * Ht + 2 * (Pi * (Dt + 2 * Eat) ^ 2 / 4)

Mt = (Vt / 1000) * Ro

End Sub

Sub Escritura_de_Titulos(Fila_Comienzo_Titulos)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 2) = “Theta, h”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 3) = “Ut, w/m2-oC”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 4) = “Ts, oC”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 5) = “Tt, oC”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 6) = “x, m”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 7) = “ml,kg/s”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 8) = “I,Amp”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 9) = “Tranferencia de serpentin a tanque, W”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 10) = “Carga, W”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 11) = “Perdida, W”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 12) = “Calefaccion, W”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 13) = “Suma, W”

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Titulos, 14) = “Ins, W/m2”

End Sub

Sub Escritura_Resultados(Fila_Comienzo_Resultados, SumTheta, Ut, Ts, Tt, Incr_Imp, ThetaMax, x, mc, CpH2O, Tc, Ac, NP, ml, I, TSR, horas, Sum_ml)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 2) = 6.4 + (SumTheta / 3600)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 3) = Ut

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 4) = Ts

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 5) = Tt

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 6) = x

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 7) = ml

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 8) = I

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 9) = Round(Us * Pi * ds * x * (Ts – Tt), 1)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 10) = Round(ml * CpH2O * 1000 * (Tt – Tti), 1)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 11) = Round(Ut * At * (Tt – Ta), 1)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 12) = Round((R * I * Cos(fi)), 1)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 13) = Round((Us * Pi * ds * x * (Ts – Tt)) – (ml * CpH2O * 1000 * (Tt – Tti)) – (Ut * At * (Tt – Ta)) + (R * I * Cos(fi)), 1)

Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados, 14) = Ins(horas, FN)

Fila_Comienzo_Resultados = Fila_Comienzo_Resultados + 1

If SumTheta >= ThetaMax Then

       TiempoEjecucion = Round(Timer – TiempoComienzo)

       TiempoEjecucion = TiempoEjecucion / 60

       Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados + 1, 2) = “Tiempo de ejecucion: ” & Round(TiempoEjecucion, 2) & ”  minutos”

       Worksheets(“Sheet2”).Cells(Fila_Comienzo_Resultados + 2, 2) = “Litros Totales, demanda= ” & Round(Sum_ml, 1) & ” litros demandados en ” & Round((SumTheta / 3600), 1) & ” horas”

       End

End If

End Sub

Function k_term(Tmean)

k_term = 1E-08 * (Tmean) ^ 3 + 6E-07 * (Tmean) ^ 2 – 2E-05 * Tmean + 0.0224

End Function

Sub Calculo_de_Ut(Dt, Eat, Ut, Tmean)

r2 = (Dt / 2) + Eat

r1 = Dt / 2

Ut = (k_term(Tmean)) / (r2 * Application.WorksheetFunction.Ln(r2 / r1))

End Sub

Sub Demanda_Domestica_de_Agua(SumTheta, ml)

Tiempo_en_Horas = (SumTheta / 3600) + 6

If Tiempo_en_Horas >= 6.5 And Tiempo_en_Horas < 7.5 Then

ml = 0.053   ‘DUCHAS, kg/s

ElseIf Tiempo_en_Horas >= 7.5 And Tiempo_en_Horas < 8.5 Then

ml = 0.00525 ‘LAVAPLATOS, kg/s

ElseIf Tiempo_en_Horas > 8.5 And Tiempo_en_Horas < 13 Then

ml = 0

ElseIf Tiempo_en_Horas >= 13 And Tiempo_en_Horas < 14 Then

ml = 0.00525 ‘LAVAPLATOS, kg/s

ElseIf Tiempo_en_Horas >= 14 And Tiempo_en_Horas < 15 Then

ml = 0.021 ‘LAVADO ROPA, kg/s

ElseIf Tiempo_en_Horas >= 15 And Tiempo_en_Horas < 22 Then

ml = 0

ElseIf Tiempo_en_Horas >= 22 And Tiempo_en_Horas < 23 Then

ml = 0.00525 ‘LAVAPLATOS, kg/s

ElseIf Tiempo_en_Horas >= 23 And Tiempo_en_Horas <= 24 Then

ml = 0

End If

End Sub

Sub Control_PID_de_Corriente(Ttset, Ilim, BM, I, Tt, ml, DeltaTheta, Kc, Ti, DeltaTt, Td)

Corr = (Ttset – Tt) * DeltaTheta

SumCorr = SumCorr + Corr

‘Calefactor funciona solo cuando hay demanda

       If ml <> 0 Then

            I = Kc * (Ttset – Tt) + (Kc / Ti) * SumCorr – Td * (Tt / DeltaTheta)

        ElseIf ml = 0 Then

            I = 0

        End If

If I > Ilim Then

I = Ilim

ElseIf I <= 0 Then

I = 0

SumCorr = 0

End If

End Sub

Function Ins(horas, FN)

    If horas <= 17.56 And horas <> 0 Then

       Ins = (1 – FN) * (-0.0169 * (horas) ^ 6 + 1.2177 * (horas) ^ 5 – 35.286 * (horas) ^ 4 + 524.86 * (horas) ^ 3 – 4242.6 * (horas) ^ 2 + 17983 * horas – 31616)

   Else

        Ins = 0

    End If

End Function

Trabajos Consultados

Industrial Engineering Letters4.12 (2014): 72-79.

Delfin Silió Sacines, et. al. «Simulation of a Solar Domestic Water Heating System With Different Collector Efficiencies and Different Volume Storage Tanks.» RE&PQJ,Vol1.No.2(2004): 124-130.

Struckmann, Fabio. «Analysis of a Flat-plate Solar Collector.» (2008).

Publicado en Sin categoría

Diseño Básico de un Reactor Para el Caso de una Reacción Instantánea

Introducción

Se presenta el caso general de una reacción de velocidad finita, y luego el caso de una reacción infinitamente rápida, para realzar las diferencias entre el primer caso y el segundo.

El caso de reacción instantánea se inspira en un problema propuesto por (Smith & Corripio, 1995).

Caso más común: La reacción no es instantánea

Generalmente, lo primero que el diseñador hace es buscar, o establecer experimentalmente la expresión cinética que gobierna la reacción que se va a llevar a cabo en el reactor.

En el caso de una reacción suficientemente rápida para que se pueda clasificar como “instantánea”, la expresión cinética no existe, o es irrelevante, porque el valor de la constante de reacción es “infinito”.

Lo anterior significa que, como “infinito” es un concepto, y no un número, no es posible construir ninguna expresión cinética como la siguiente, por ejemplo, que es la que se aplica en la mayoría de los casos.

Ecuacion1

Como en la ecuación precedente Ci es la concentración del reaccionante i expresada en masa por unidad de volumen (M/L3), k es una constante cinética, y n es el orden de reacción respecto del reaccionante Ci.

Para el caso de la forma cinética de la ecuación anterior, las unidades de k se pueden obtener despejando k de la ecuación anterior y substituyendo las unidades de las variables que quedan a la derecha del signo igual.

Cuando se despeja k, se obtiene la siguiente ecuación:

Ecuacion2

Si se escribe la ecuación anterior en función de las unidades de cada facto, se obtiene la siguiente expresión para la unidades de k.

Ecuacion3

Si reemplazamos estas unidades en la ecuación cinética, obtendremos lo siguiente.

Ecuacion4

También se sabe que a la expresión anterior debe añadírsele los flujos másicos de reaccionantes que ingresan al reactor, y el flujo másico que sale de éste y el flujo másico de la especie que desaparece por reacción.

Un ejemplo

Supóngase que se trata de formular los balances para una reacción que se represente mediante la siguiente estequiometria:

Ecuacion5

y que, en la ecuación anterior, a y b representan el número de moles de las especies A, B, que se requiere para formar c moles de C.

Si se acepta que lo expresado anteriormente ocurre, y que la reacción se lleva a cabo en un reactor perfectamente agitado las velocidades de reacción estarán relacionadas a los coeficientes estequiométricos (Wikipedia, 2018) de la siguiente forma.

Ecuacion6

La expresión anterior es útil porque permite expresar la velocidad de desaparición/aparición de cualquiera de las especies que toman parte en la reacción si se conoce la expresión cinética de una de ellas, que generalmente es la del producto.

Utilizando estos conceptos es posible formular balances diferenciales de masa para cada una de las especies químicas que intervienen en la reacción.

Si se acepta que FA denota el flujo volumétrico de la especie A, que FB denota el flujo volumétrico de B, y que FC denota el flujo volumétrico que sale del reactor (que de acuerdo al tiempo medio de residencia V/F y a la velocidad de reacción, que son dos parámetros que determinan la conversión puede contener A, B y C), se puede formular la siguiente ecuación para las especies que intervienen en la reacción, en función de la velocidad de formación de la especie C, que -para fines del ejemplo- se asume que puede representarse por medio de la expresión siguiente.

Ecuacion7

Las ecuaciones para las especies A, B, y C pueden formularse de la manera siguiente.

Ecuacion8

Ecuacion9

Ecuacion10

Las condiciones iniciales para las ecuaciones precedentes pueden ser las siguientes: hay dos opciones para el valor inicial de V que dependen de si la simulación comienza con el reactor lleno, o con el reactor vacío, significando con ello que no abarque el tiempo de llenado y el tiempo de estabilización del nivel de sus contenidos, o que sí lo haga.

@ t=0 CA=CAo, CB=CBo, CC=0, V=V, ó V=0

La ecuación restante es la que relaciona el volumen de reacción con el tiempo, y es la siguiente.Ecuacion11

En la última ecuación dV/dt es diferente de cero sólo cuando el volumen inicial es cero, y/o cuando se produce una perturbación en FA, FB, o en ambos,  lo que el analista hace para estudiar la velocidad de recuperación del reactor respecto de diferentes valores de la perturbación a parámetros de control constantes, o vice-versa.

Segundo Caso: La Reacción es instantánea

Si la reacción es instantánea entonces, como ya se ha dicho, no es posible formular ecuaciones cinéticas como las que se muestran en la sección anterior, y se debe asumir (por lo menos eso es lo que se ha hecho en este artículo) que las especies que intervienen en la reacción forman el producto instantáneamente, de acuerdo a la estequiometria de que se trate, o de acuerdo a proporciones empíricas preestablecidas.

La estequiometria

La estequiometria es muy sencilla y se puede encontrar en la literatura muy fácilmente. Se basa en que, para obtener nitrato de amonio se debe hacer reaccionar amoníaco en estado gaseoso con ácido nítrico, lo que -como es de imaginar- desprende una cantidad de energía considerable.

Otra manera de producir nitrato de amonio consiste en producir previamente una solución de amoníaco al 30% y hacer que ésta reaccione con una solución de ácido nítrico de 55%.

Este último camino es el que se ha utilizado en este artículo.

La estequiometria se resume en la siguiente tabla.

 

 

 

EstequiometriaTablaEcel

Ilustración 1. Cálculos estequiométricos 

Breve explicación de la Ilustración 1

  1. Se parte de una estadística acerca del uso de fertilizantes químicos que estimaba que, en el año 2014, se requerían 105,000 toneladas métricas (TM) anuales (www.fao.org/faostat/en/#data/RF/visualize, s.f.) de fertilizantes químicos.
  2. Se estimó, para propósitos de construir el caso que, de la cantidad citada, se requerirían 96,000 TM/año de nitrato de amonio de 100%, que la instalación industrial trabajaría 8,000 horas/año, y que cada reactor produciría un décimo del total requerido, lo que implicaría 1,200 kg/h, que equivalen a 3.33 x 10-1 kg/s de nitrato de amonio.
  3. Como el reactor produce una solución de nitrato de amonio de 66.5% la cantidad producida por año es igual a (1,200+603.8)kg/h=14,430 TM/año.
  4. De esta última cantidad, como la tabla se construyó en Excel, se calcula mediante estequiometria simple las masas de 100% de amoniaco y ácido nítrico.
  5. Se calcula luego la masa de agua necesaria para reducir los reaccionantes a soluciones de 30%, y 55%, respectivamente.
  6. Se suman las cantidades de reactivo y las cantidades de agua obteniéndose la masa de nitrato de amonio de la que se partió, y la concentración porcentual de la solución en peso.

Conceptualización del programa

El programa simula la alimentación de las masas por segundo de las soluciones, a las concentraciones citadas, asumiéndose que la alimentación se convierte instantáneamente en la solución de nitrato de amonio ya mencionada en cuanto llega a la masa reaccionante contenida en el reactor perfectamente agitado, pero -al mismo tiempo- la solución acuosa del producto tiene la densidad de una solución de nitrato de amonio al 66.5%.

Los datos se presentan en kg/s, y los resultados del cálculo se presentan en m3/s. Estos flujos se obtienen de la división de las masas alimentadas para sus respectivas densidades, que -en cada caso- se obtuvieron de fuentes bibliográficas.

El flujo volumétrico de producto se obtiene (reacción instantánea) dividiendo las masas alimentadas para la densidad de la solución de nitrato de amonio, que no es igual a las densidades de las soluciones que se alimentan.

El programa tiene un valor de hset, que corresponde al nivel de operación del reactor, fijándose de esta manera el volumen del reactor, que puede cambiarse, y que -en la práctica- se maneja mediante una alarma de nivel alto (LAH), cuyo valor deseado, que puede ser cambiado por el operador, determina el nivel de la mezcla reaccionante  en el reactor.

Para fines de correr el programa el nivel se fijó en 0.50 m, pero podría fijarse cualquier otro valor lo que tan solo implicaría cambiar el dato en la interface de usuario, que se presenta a continuación. En realidad el nivel de la mezcla reaccionante debe especificarse para estimar apropiadamente la potencia de agitación, así como el tipo y número de paletas de agitación, y de agitadores.

InterfaceUsuario

Ilustración 2. Interface de usuario del programa

La densidad de la solución de nitrato de amonio se obtuvo por extrapolación lineal y ajuste de la recta resultante, debido a que en la bibliografía sólo se presentan densidades hasta el 50% en peso, porque el nitrato de amonio en solución se comercializa en soluciones de 50%, como máximo.

El gráfico de extrapolación y ajuste se muestra a continuación.

ExtrapolacionLineal

Ilustración 3. Extrapolación y ajuste de los valores extrapolados de densidad de nitrato de amonio a diferentes concentraciones porcentuales

Resultados

Los resultados de la ejecución del programa se muestran a continuación, en forma de tabla, y en forma gráfica.

ResultadosTabulares

Ilustración 4. Resultados de la ejecución del programa

ResultadosGraficos

Ilustración 5. Gráfica de los resultados de la ejecución del programa

ResultadosGraficos2

Ilustración 6. Resultados de la ejecución del programa: Se puede observar el nivel de operación

Control del Reactor

En el reactor se controlan los flujos mediante un controlador PID, y dos controladores PI. La subrutina que los maneja es general para un controlador PID.

Los controladores de alimentación son dos controladores que deberían haberse simulado como de proporción, pero que terminaron operando independientemente (cuestión de una línea extra de código, que no se añadió).

El programa de simulación

A continuación se muestra el código del programa.

Option Explicit

Dim F31, F32, Tiempo, Deltah As Variant

Dim DeltaFHNO3, TdFHNO3, DeltaFNH3, TdFNH3, DeltaFNH4NO3, TdFNH4NO3, DeltaF1, TdF1, DeltaF2, TdF2, DeltaF3, TdF3 As Variant

Dim Kc2, Ti2, F2set As Variant

Dim L, vH2O, Ds, CpR, Tempwo, RoH2O, U, DeltaHr, CpH2O, Taster As Variant

Dim Tiempo_Marcador As Boolean

Dim ffF, KcF, TiF, SumCorrF, Ffset, Ff, pF, FNH3set As Variant

Dim FHNO3set, FHNO3, KcHNO3, TiHNO3 As Variant

Dim siga, VR, hf, Constaprint, SumF3, SumFNH4NO3, FlujoMedioF3, FlujoMedioNH4NO3, RoNH4NO3 As Variant

Dim DeltaV, V, Fs, F, p, Corr, hRset, BD As Variant

Dim Pi, F1, F2, F3, Kv1, Kv2, Kv3, Deltap1, Deltap2, Deltap3, ff1, ff2, ff3, Ro1, Ro2, Ro3, T As Variant

Dim ffHNO3, ffNH3, ffNH4NO3, FNH3, FNH4NO3, FsetHNO3, FsetNH3, FNH4NO3set As Variant

Dim Fset, Kv, Kc, Ti, SumCorr, Deltap As Variant

Dim TiNH3, KcNH4NO3, TiNH4NO3, DeltaPHNO3, DeltaPNH3, DeltaPNH4NO3, KcNH3, RoHNO3, RoNO3, SumCorrHNO3, SumCorrNH3, SumCorrNH4NO3, SumCorr1, SumCorr2, SumCorr3 As Variant

Dim Deltat, Tmax, DR, hR, hR_Fisica As Variant

Dim KvHNO3, KvNH3, KvNH4NO3, Deltaprint, RoNH3, Fila As Variant

Sub Nirato_de_Amonio2()

‘Ingreso de variables y parámetros

KvHNO3 = Worksheets(“Hoja2”).Cells(2, 3)

KvNH3 = Worksheets(“Hoja2”).Cells(3, 3)

KvNH4NO3 = Worksheets(“Hoja2”).Cells(4, 3)

L = Worksheets(“Hoja2”).Cells(8, 3)

Ds = Worksheets(“Hoja2”).Cells(9, 3)

CpR = Worksheets(“Hoja2”).Cells(10, 3)

Tempwo = Worksheets(“Hoja2”).Cells(11, 3)

vH2O = Worksheets(“Hoja2”).Cells(12, 3)

RoH2O = Worksheets(“Hoja2”).Cells(13, 3)

DeltaPHNO3 = Worksheets(“Hoja2”).Cells(14, 3)

DeltaPNH3 = Worksheets(“Hoja2”).Cells(15, 3)

CpH2O = Worksheets(“Hoja2”).Cells(17, 3)

U = Worksheets(“Hoja2”).Cells(18, 3)

DeltaHr = Worksheets(“Hoja2”).Cells(19, 3)

RoHNO3 = Worksheets(“Hoja2”).Cells(20, 3)

RoNH3 = Worksheets(“Hoja2”).Cells(21, 3)

RoNH4NO3 = Worksheets(“Hoja2”).Cells(22, 3)

DR = Worksheets(“Hoja2”).Cells(23, 3)

hR_Fisica = Worksheets(“Hoja2”).Cells(24, 3)

DeltaPNH4NO3 = Worksheets(“Hoja2”).Cells(27, 3)

hR = Worksheets(“Hoja2”).Cells(28, 3)

DR = Worksheets(“Hoja2”).Cells(29, 3)

ffHNO3 = Worksheets(“Hoja2”).Cells(33, 3)

ffNH3 = Worksheets(“Hoja2”).Cells(34, 3)

hRset = Worksheets(“Hoja2”).Cells(36, 3)

BD = Worksheets(“Hoja2”).Cells(37, 3)

Deltat = Worksheets(“Hoja2”).Cells(40, 3)

Tmax = Worksheets(“Hoja2”).Cells(41, 3)

Deltaprint = Worksheets(“Hoja2”).Cells(42, 3)

Constaprint = Worksheets(“Hoja2”).Cells(43, 3)

FHNO3set = Worksheets(“Hoja2”).Cells(44, 3)

KcHNO3 = Worksheets(“Hoja2”).Cells(45, 3)

TiHNO3 = Worksheets(“Hoja2”).Cells(46, 3)

KcNH3 = Worksheets(“Hoja2”).Cells(47, 3)

TiNH3 = Worksheets(“Hoja2”).Cells(48, 3)

FNH3set = Worksheets(“Hoja2”).Cells(49, 3)

KcNH4NO3 = Worksheets(“Hoja2”).Cells(50, 3)

TiNH4NO3 = Worksheets(“Hoja2”).Cells(51, 3)

TdFNH4NO3 = Worksheets(“Hoja2”).Cells(52, 3)

FNH4NO3set = Worksheets(“Hoja2”).Cells(53, 3)

Pi = 3.14159

SumCorrHNO3 = 0: SumCorrNH3 = 0: SumCorrNH4NO3 = 0

FlujoMedioNH4NO3 = 0: DeltaFHNO3 = 0: TdFHNO3 = 0: DeltaFNH3 = 0: TdFNH3 = 0: DeltaFNH4NO3 = 0

‘========================================

‘Limpieza

Worksheets(“Hoja3”).Range(“a1:j17”).Clear

Worksheets(“Hoja3”).Cells(1, 2) = “KcNH4NO3= ” & KcNH4NO3 & ” TiNH4NO3= ” & TiNH4NO3 & ”  TdFNH4NO3= ” & TdFNH4NO3 & ”  KvNH4NO3 = ” & KvNH4NO3 & “DeltaPNH4NO3= ” & DeltaPNH4NO3

Worksheets(“Hoja3”).Range(“A3:j15”).HorizontalAlignment = xlCenter

Call Recuadro

‘Limpieza2

Worksheets(“Hoja3”).Range(“a15:j17”).Clear

Fila = 3

Call Titulos(Fila)

ffHNO3 = 0#: ffNH3 = 0: ffNH4NO3 = 0: hR = 0: siga = 1: _

VR = 0: SumCorr = 0: Pi = 3.14159: SumFNH4NO3 = 0: Tiempo_Marcador = True

‘************************************************

For T = 0 To Tmax Step Deltat

If T = 0 Then

FHNO3 = 0: FNH3 = 0: FNH4NO3 = 0

Call Resultados(T, hR, FHNO3, FNH3, _

FNH4NO3, ffHNO3, ffNH3, ffNH4NO3, RoHNO3, RoNH3, RoNH4NO3)

GoTo Fin

End If

Call Balance_de_Masa2(ffHNO3, ffNH3, ffNH4NO3, RoHNO3, RoNH3, RoNH4NO3, Deltat _

, DeltaPHNO3, DeltaPNH3, DeltaPNH4NO3, KvHNO3, KvNH3, KvNH4NO3, hR, _

DR, FHNO3, FNH3, FNH4NO3, hRset, BD, siga, VR, hR_Fisica, Kc, Ti, SumCorr, _

Pi, FlujoMedioNH4NO3, SumFNH4NO3, T, Taster, FHNO3set, KcHNO3, TiHNO3, SumCorrHNO3, KcNH3, TiNH3, SumCorrNH3, FNH3set _

, KcNH4NO3, TiNH4NO3, SumCorrNH4NO3, FNH4NO3set, DeltaFHNO3, TdFHNO3, DeltaFNH3, TdFNH3, DeltaFNH4NO3, TdFNH4NO3)

If (hR >= hRset) And (Tiempo_Marcador = True) Then

Taster = T

Tiempo_Marcador = False

ElseIf hR >= hR_Fisica Then

End

End If

If T >= Deltaprint Then

Call Resultados(T, hR, FHNO3, FNH3, _

FNH4NO3, ffHNO3, ffNH3, ffNH4NO3, RoHNO3, RoNH3, RoNH4NO3)

Deltaprint = Deltaprint + Constaprint

End If

Fin:

Next T

‘========================================

Call Resultados(T, hR, FHNO3, FNH3, _

FNH4NO3, ffHNO3, ffNH3, ffNH4NO3, RoHNO3, RoNH3, RoNH4NO3)

FlujoMedioNH4NO3 = (SumFNH4NO3 / RoNH4NO3) / (Tmax – Taster)

Worksheets(“Hoja3”).Cells(Fila, 4) = “Flujo medio NH4NO3=  ” & Round(FlujoMedioNH4NO3, 7) & ”  m^3/s ”

Worksheets(“Hoja3”).Cells(Fila + 1, 3) = “Producción es : ” & Round(SumFNH4NO3 / (Tmax – Taster), 2) & ” kg/s de NH4NO3 de 66.5%” & ”  Produccion Anual Estimada es ” & Round(SumFNH4NO3 / (Tmax – Taster), 2) * 3600 * 8000 / 1000 & ” TM/año”

Worksheets(“Hoja3”).Cells(Fila + 2, 3) = “Taster= ” & Round(Taster, 0) & ” s”

End Sub

Sub Titulos(Fila)

Worksheets(“Hoja3”).Cells(Fila, 2) = “T,s”

Worksheets(“Hoja3”).Cells(Fila, 3) = “hR,m”

Worksheets(“Hoja3”).Cells(Fila, 4) = “FHNO3 alimentado al reactor, m3/s”

Worksheets(“Hoja3”).Cells(Fila, 5) = “ffHNO3”

Worksheets(“Hoja3”).Cells(Fila, 6) = “FNH3 alimentado al reactor, m3/s”

Worksheets(“Hoja3”).Cells(Fila, 7) = “ffNH3”

Worksheets(“Hoja3”).Cells(Fila, 8) = “NH4NO3 formado en reactor, m3/s”

Worksheets(“Hoja3”).Cells(Fila, 9) = “FNH4NO3 sale de reactor, m3/s”

Worksheets(“Hoja3”).Cells(Fila, 10) = “ffNH4NO3”

Fila = Fila + 1

End Sub

Sub Resultados(T, hR, FHNO3, FNH3, FNH4NO3, ffHNO3, ffNH3, ffNH4NO3, RoHNO3, RoNH3, RoNH4NO3)

Worksheets(“Hoja3”).Cells(Fila, 2) = T

Worksheets(“Hoja3”).Cells(Fila, 3) = hR

Worksheets(“Hoja3”).Cells(Fila, 4) = (FHNO3 / RoHNO3)

Worksheets(“Hoja3”).Cells(Fila, 5) = ffHNO3

Worksheets(“Hoja3”).Cells(Fila, 6) = (FNH3 / RoNH3)

Worksheets(“Hoja3”).Cells(Fila, 7) = ffNH3

Worksheets(“Hoja3”).Cells(Fila, 8) = (FHNO3 + FNH3) / RoNH4NO3

Worksheets(“Hoja3”).Cells(Fila, 9) = (FNH4NO3 / RoNH4NO3)

Worksheets(“Hoja3”).Cells(Fila, 10) = ffNH4NO3

Fila = Fila + 1

End Sub

Sub Balance_de_Masa2(ff1, ff2, ff3, Ro1, Ro2, _

Ro3, Deltat, Deltap1, Deltap2, Deltap3, _

Kv1, Kv2, Kv3, h, D, F1, F2, F3, hset, BD, siga, V, _

hf, Kc, Ti, SumCorr, Pi, FlujoMedioF·, SumF3, T, Taster, F1set, Kc1, Ti1, SumCorr1, Kc2, Ti2, SumCorr2, F2set, Kc3, Ti3, _

SumCorr3, F3set, DeltaF1, TdF1, DeltaF2, TdF2, DeltaF3, TdF3)

Call ControlPID_Flujo(ff1, Kc1, Ti1, SumCorr1, F1set, F1, Deltat, TdF1, DeltaF1, T, 0)

F1 = Kv1 * ff1 * Sqr(Deltap1) ‘en kg/s

Call ControlPID_Flujo(ff2, Kc2, Ti2, SumCorr2, F2set, F2, Deltat, TdF2, DeltaF2, T, 0)

F2 = Kv2 * ff2 * Sqr(Deltap2) ‘ en kg/s

If (h >= hset) Then   ‘

F31 = F3

Call ControlPID_Flujo(ff3, Kc3, Ti3, SumCorr3, F3set, F3, Deltat, TdF3, DeltaF3, T, 0)

F3 = Kv3 * ff3 * Sqr(Deltap3) ‘ En kg/s Deltap en kPa^0.5_

F32 = F3

If (F32 = 0) Or (F31 = 0) Then

DeltaF3 = 0

Else

DeltaF3 = F31 – F32

End If

SumF3 = SumF3 + F3 * Deltat   ‘Suma la producción en kg

End If

‘Actualización de h

DeltaV = (((F1 + F2) / Ro3) – (F3 / Ro3)) * Deltat ‘ Fs En m3/s

V = V + DeltaV

h = V / (Pi * D ^ 2 / 4)

End Sub

Sub ControlPID_Flujo(ffF, KcF, TiF, SumCorrF, Ffset, Ff, Deltat, TdF, DeltaF, Tiempo, Bi)

SumCorrF = SumCorrF + (Ffset – Ff) * Deltat

pF = ((KcF * (Ffset – Ff)) + ((KcF / TiF) * SumCorrF) – (KcF * TdF) * (DeltaF / Deltat)) + Bi

If pF > 103.42 Then

pF = 103.42

SumCorrF = 0

ElseIf pF < 20.68 Then

pF = 20.68

SumCorrF = 0

End If

ffF = (pF – 20.68) / (103.42 – 20.68)

End Sub

Trabajos Citados

(s.f.). Recuperado el February de 2018, de http://www.fao.org/faostat/en/#data/RF/visualize.

Smith, C., & Corripio, A. (1995). Control Automático de Procesos: Teoría y Práctica. Mexico: Limusa-Noriega Editores.

Wikipedia. (24th de February de 2018). Obtenido de https://en.wikipedia.org/wiki/Reaction_rate

 

Publicado en Sin categoría

Un Ejemplo de un Método Acortado Para Estimar las Concentraciones en Cada Plato

Introducción

Los métodos acortados pueden ser útiles para hacer estimaciones rápidas sin tener que recurrir a los métodos rigurosos tales como, por ejemplo, Lewis-Matheson, acerca del cuál escribí un artículo en este blog hace algún tiempo.

 

En este caso pensé que podría ser interesante realizar la implementación computacional de un caso que se presenta como aplicación a la solución de un método acortado en  (John H. Perry et. al, 1963) , en las páginas 13-34 y 13-35.

El problema es típico de lo que se hace en la Universidad, y se puede resumir en la siguiente tabla, que es una adaptación de la original del artículo de (John H. Perry et. al, 1963).

Primera Tabla

En este caso es obvio que A es el componente más volátil y que estará presente en su mayoría en el destilado, pero que una fracción muy pequeña podría estar presente en el destilado. También puede suponerse que una muy pequeña fracción del componente D, que es el menos volátil, estará presente en los fondos.

Esta suposición es la que origina el problema, porque -en la realidad- la magnitud de la traza debe estimarse por tanteo, y el valor correcto (que yo lo tomé de la solución del problema que consta en la referencia) debe validar la composición del destilado, que es asumida.

Para resolver el problema se hace un cálculo desde la caldereta que es la etapa cero (j=0) hasta la etapa número seis, información que en principio no se conoce, pero que luego de calcular con el valor de convergencia del ejemplo,

se verifica (¡imagínense la magnitud del trabajo en la época de la publicación, en que por declaración explícita del autor él no disponía de computador para realizar los cálculos!)

El autor del problema también indica que la etapa de alimentación es la cuarta, contando de abajo hacia arriba porque en la cuarta etapa la composición de B y C será muy similar, como se puede comprobar examinando los resultados de la solución del problema.

Supongo que el o los autores tal vez hayan poseído la columna del problema, en la que podían variar la alimentación, la etapa de alimentación, y la tasa de reflujo, y que hayan realizado el trabajo experimental que les permitió desarrollar el método que aquí se presenta informatizado.

En el planteo del problema se proponen la siguiente configuración de la operación de la columna que se alimenta con líquido saturado, y que se opera con una tasa de reflujo igual a cuatro, la misma que es adiabática, y en la que prevalece “constant molal overflow” (no encontré la traducción en la Internet) y el único curso que hice en destilación lo realicé en USA, en el siglo pasado.

Configuración Columna2

De la configuración operacional, y de la siguiente ecuación permiten calcular la composición de vapor de los fondos, y también la composición del vapor de las bandejas.

Ecuacion1

Ecuacion2Ésta última ecuación, y sus coeficientes son una consecuencia de balance de masa sobre la sección de empobrecimiento de la columna, que se obtiene, de acuerdo al siguiente razonamiento, que comienza por establecer la retención o hold-up de los fondos mediante un sencillo balance de masa a estado estacionario, que es el siguiente

Ecuacion3

De acuerdo a la anterior ecuación, la retención, sin el puntito, es igual a 50 moles, lo que quiere decir que la proporción entre el líquido que baja en la zona de empobrecimiento es 300/50=6, y que la proporción entre el vapor que asciende y la retención en fondos es igual a 250/50=5, lo que quiere decir que la ecuación de cálculo para la x1,i es la siguiente:

Ecuacion4

Siempre de acuerdo al artículo citado, los valores de y para la sección de empobrecimiento se calculan utilizando la siguientes ecuaciones.

Para el caso de los componentes del vapor en equilibrio en cada plato de la sección de empobrecimiento es la ecuación es la siguiente.

Ecuacion5Y Para el caso de las fracciones molares en fase líquida la ecuación es.

Ecuacion6

Para el caso de la sección de enriquecimiento las ecuaciones y son las siguientes obteniéndose loa coeficientes mediante un balance de masa.

Cabe enfatizar que la ecuación general de la sección de empobrecimiento está anclada a la composición de los fondos, y la de enriquecimiento está anclada a la composición asumida del destilado, como se puede observar.

Ecuacion7

y

Ecuacion8

Para el cálculo de los valores de YDcalc es imprescindible haber ejecutado el programa, para contar con los valores y6,i como se indica en la ecuación siguiente, ya que éste valor  es necesario calcular el valor de YDcalci, lo que se realiza por medio de la ecuación siguiente.

Ecuacion9

El valor de x7,i se calcula utilizando la siguiente ecuación.

Ecuacion10

El Programa

A continuación muestro el programa, para que el uso y rol de cada una de las ecuaciones se aclare pueda observar.

Option Explicit

Dim ArrayAlfa(1 To 4), Arrayx(0 To 7, 1 To 4), Arrayy(0 To 6, 1 To 4), ArrayXw(1 To 4) As Variant

Dim SumAlfax, Fila, Sumx, Sumy, Sumaalfax, ArrayYD(1 To 4), ArrayYDcal(1 To 4) As Variant

Dim i, j, k, No_vale As Integer

Sub metodo_acortado()

ArrayAlfa(1) = Worksheets(“Hoja1”).Cells(7, 3)

ArrayAlfa(2) = Worksheets(“Hoja1”).Cells(8, 3)

ArrayAlfa(3) = Worksheets(“Hoja1”).Cells(9, 3)

ArrayAlfa(4) = Worksheets(“Hoja1”).Cells(10, 3)

Arrayx(0, 1) = Worksheets(“Hoja1”).Cells(11, 3) ‘Valor estimado de la traza

Arrayx(0, 2) = Worksheets(“Hoja1”).Cells(12, 3) ‘El subíndice 1 se refiere a fondos

Arrayx(0, 3) = Worksheets(“Hoja1”).Cells(13, 3)

Arrayx(0, 4) = Worksheets(“Hoja1”).Cells(14, 3)

For i = 1 To 4

ArrayXw(i) = Arrayx(0, i)

Next i

 

ArrayYD(1) = 0.5

ArrayYD(2) = 0.48

ArrayYD(3) = 0.02

ArrayYD(4) = 0

 

‘Calculo de producto de alfa por fracción molar de fondos

Fila = 3

‘Limpieza

Worksheets(“Hoja2”).Range(“A3:L12”).Clear

Call Escritura_Titulos(Fila)

Call Cálculos_de_Caldereta(0, ArrayAlfa(), Arrayx(), Arrayy(), ArrayXw()) ‘j=0

Call Escritura_Resultados_Caldereta(0, Fila, Arrayx(), Arrayy()) ‘j=0

Call Calculo_X_y_Y(Arrayx(), Arrayy(), ArrayXw(), ArrayYD())

Call Calculo_de_Ycalc(Arrayy(), ArrayYD(), Arrayx(), ArrayYDcal())

 

End Sub

Sub Escritura_Titulos(Fila)

Worksheets(“Hoja2”).Cells(Fila, 1) = “j”

For i = 1 To 5

If i <= 4 Then

Worksheets(“Hoja2”).Cells(Fila, 1 + i) = “x(” & “j” & “,” & ” i ” & “)”

ElseIf i = 5 Then

Worksheets(“Hoja2”).Cells(Fila, 1 + i) = “Sumax”

End If

Next i

For k = 1 To 5

If k <= 4 Then

Worksheets(“Hoja2”).Cells(Fila, 6 + k) = “y( ” & “j” & “,” & “i” & “)”

 

ElseIf k = 5 Then

Worksheets(“Hoja2”).Cells(Fila, 6 + k) = “Sumay”

End If

Next k

Fila = Fila + 1

End Sub

Sub Escritura_Resultados_Caldereta(j, Fila, Arrayx(), Arrayy())

Sumx = 0

Worksheets(“Hoja2”).Cells(Fila, 1) = j

For i = 1 To 4

Worksheets(“Hoja2”).Cells(Fila, i + 1) = Arrayx(j, i)

Sumx = Sumx + Arrayx(j, i)

Next i

Worksheets(“Hoja2”).Cells(Fila, 6) = Sumx

Sumy = 0

For k = 7 To 10

Worksheets(“Hoja2”).Cells(Fila, k) = Arrayy(j, k – 6)

Sumy = Sumy + Arrayy(j, k – 6)

Next k

Worksheets(“Hoja2”).Cells(Fila, 11) = Sumy

Fila = Fila + 1

End Sub

Sub Escritura_Resultados_Etapas(j, Fila, Arrayx(), Arrayy())

Worksheets(“Hoja2”).Cells(Fila, 1) = j

Sumx = 0

For i = 1 To 4

Worksheets(“Hoja2”).Cells(Fila, i + 1) = Arrayx(j, i) ‘ / 5

Sumx = Sumx + Arrayx(j, i) ‘ / 5

Next i

Worksheets(“Hoja2”).Cells(Fila, 6) = Sumx

Sumy = 0

k = 7

For i = 1 To 4

Worksheets(“Hoja2”).Cells(Fila, k) = Arrayy(j, i)

Sumy = Sumy + Arrayy(j, i)

k = k + 1

Next i

Worksheets(“Hoja2”).Cells(Fila, 11) = Sumy

Fila = Fila + 1

End Sub

Sub Cálculos_de_Caldereta(j, ArrayAlfa(), Arrayx(), Arrayy(), ArrayXw())

SumAlfax = 0

For i = 1 To 4

SumAlfax = SumAlfax + ArrayAlfa(i) * Arrayx(j, i)

ArrayXw(i) = Arrayx(j, i)

Next i

For i = 1 To 4

Arrayy(j, i) = ((ArrayAlfa(i) * Arrayx(j, i)) / SumAlfax)

Next i

End Sub

Sub Calculo_X_y_Y(Arrayx(), Arrayy(), ArrayXw(), ArrayYD())

For j = 1 To 6

If j <= 4 Then

For i = 1 To 4

Arrayx(j, i) = ((5 / 6) * Arrayy(j – 1, i) + (1 / 6) * ArrayXw(i))

Next i

Sumaalfax = 0

For i = 1 To 4

Sumaalfax = Sumaalfax + 6 * ArrayAlfa(i) * Arrayx(j, i)

Next i

For i = 1 To 4

Arrayy(j, i) = (6 * ArrayAlfa(i) * Arrayx(j, i)) / (Sumaalfax)

Next i

ElseIf j > 4 Then

For i = 1 To 4

Arrayx(j, i) = (5 / 4) * Arrayy(j – 1, i) – ArrayYD(i) / 4

Next i

Sumaalfax = 0

For i = 1 To 4

Sumaalfax = Sumaalfax + 4 * ArrayAlfa(i) * Arrayx(j, i)

Next i

For i = 1 To 4

Arrayy(j, i) = (4 * ArrayAlfa(i) * Arrayx(j, i)) / Sumaalfax

Next i

End If

Call Escritura_Resultados_Etapas(j, Fila, Arrayx(), Arrayy())

Next j

End Sub

Sub Calculo_de_Ycalc(Arrayy(), ArrayYD(), Arrayx(), ArrayYDcal())

For i = 1 To 4

Arrayx(7, i) = (5 / 4) * Arrayy(6, i) – ArrayYD(i) / 4

Next i

For i = 1 To 4

ArrayYDcal(i) = 5 * Arrayy(6, i) – 4 * Arrayx(7, i)

Worksheets(“Hoja2”).Cells(Fila, 1) = “YDcalc(” & i & “)=” & ArrayYDcal(i)

Fila = Fila + 1

Next i

End Sub

La interface de usuario

A continuación se muestra la interface de usuario, que es muy modesta. Los otros valores los definí dentro del programa, lo que – a decir verdad- no es muy buena práctica.

Interface de usuario

Los resultados

Resultados

Como se es posible apreciar las fracciones molares de los compuestos B y D son muy similares en la cuarta etapa, y es por eso que los autores la escogieron como etapa de alimentación, como ya hemos mencionado antes.

Trabajo citado

John H. Perry et. al. (1963). Chemical Engineers’ Handbook. Newy York: McGraw-Hill Book Company.

 

 

Publicado en Sin categoría

El Método de Lewis y Matheson para el Cálculo de Columnas de Destilación

La problemática del cálculo de las columnas de destilación de multicomponentes consiste en que para definir la columna que a uno le interesa, es necesario conocer la composición y la temperatura de cada plato o etapa, y en que cuando se trata de multicomponentes las sumas de las fracciones molares tanto en la fase líquida como en la de vapor deben ser iguales a la unidad.

En el caso de que los componentes sean sólo dos el problema se reduce porque si una de las fracciones es x1 la otra necesariamente deberá ser x2-1, lo que no ocurre en el caso de los sistemas multicomponentes en donde la suma de tres o más fracciones consiste de un número infinito de conjuntos de términos cuya suma puede ser igual a uno, sin embargo de que, de ellos, habrá uno solo que sea termodinámicamente correcto.

Por las razones anteriores Lewis y Matheson (http://www.spq.pt/magazines/RPQ/282/article/753/pdf) crearon un método que consiste en especificar la composición del destilado, y la alimentación a la columna, que se basa en el cálculo de la composición de cada plato partiendo desde el condensador, en el caso de la sección de enriquecimiento; y desde la caldereta en el caso de la sección de empobrecimiento, hasta el palto de alimentación.

El método se basa en un cálculo iterativo de una composición del destilado sobre la que se sigue iterando hasta que la composición de la etapa de alimentación sea exactamente igual para la sección de enriquecimiento y para la sección de empobrecimiento.

El método se basa también, para el caso de un condensador total, en que para el plato número uno, la composición del destilado en fase de vapor es igual a las fracciones molares del destilado en fase líquida, por el simple hecho de que lo que sale hacia afuera de la columna vía la corriente D es de igual composición a lo que sale del plato superior en fase de vapor.

A partir de este valor de las composiciones en fase de vapor se calcula, mediante las volatilidades relativas de los componentes, que se consideran constantes, las composiciones de la fase líquida, hasta llegar al plato de alimentación desde ambos extremos de la columna.

Después de haber implementado los cálculos desde el plato superior, número uno hasta el plato de alimentación el algoritmoo especifica que se debe comenzar un cálculo independiente desde los fondos de la columna, que permite obtener la composición de la corriente B, y las fracciones molares de los fondos, lo que a su vez, permite calcular la composición de la etapa o plato de alimentación independientemente desde abajo.

Una vez que se ha efectuado lo que queda explicado existe un método (Holland, 1963) que permite, mediante el uso de valores de las XDi del fin del cálculo, que muestran el el programa que acompaña este artículo, y –en último término- en último resultado que se obtiene, que es el de convergencia, que más adelante se consigna en la Figura 1.

El método se reduce entonces a calcular un conjunto de valores que constituyen nuevas aproximaciones a los valores de las fracciones molares de destilación, que deben insertarse como datos para repetir el cálculo, hasta obtener los resultados finales de convergencia.

En general no he encontrado una regla general que asegure la convergencia. En este caso particular se usaron los valores de un perfil de temperaturas totalmente supuesto, y los valores de alfa correspondientes al quinto plato.

El método es de convergencia lenta.

Figura 1

Figura 1. Composiciones de convergencia de los platos de la columna

La figura anterior constituye la última iteración del programa. Esto se sabe porque los “nuevos valores” de XND(1), XND(2), y XND(3) son iguales a las fracciones de vapor correspondientes al plato número 1, o a las fracciones molares del plato número uno en fase líquida (véase la Figura 2).

Figura 2

Figura 2. Interfase de usuario del programa

Para el despliegue de los contenidos de la Figura 1 se utilizó un macro que muestra las composiciones de la etapa cinco calculadas desde los fondos hasta la etapa de alimentación, y desde el condensador hasta la etapa mencionada, en color rojo, donde se puede apreciar que los dos cálculos coinciden en todas sus cifras.

Como comentario general debo indicar que en la Internet se pueden encontrar infinidad de tablas de valores para los coeficientes de Antoine, que difieren y que, como en el cálculo como forman parte de un exponente ejercen gran influencia en los resultados del programa. Yo utilicé los del National Institute of Standarization and Chemistry (http://webbook.nist.gov/cgi/cbook.cgi?ID=C112403&Mask=FFF) que me permito recomendar altamente.

El Programa Fuente

Option Explicit

Dim ArrayTc(1 To 11), ArrayTcn(1 To 11)

Dim ArrayxS(5 To 11, 1 To 3), ArrayyS(5 To 11, 1 To 3), ArrayxSB(1 To 11, 1 To 3), SumAlfa, Z, SumBXB As Variant

Dim ArrayXDN(1 To 3) As Variant

Dim Arrayv_over_bi(1 To 5, 1 To 3), Arrayv_over_di(1 To 5, 1 To 3), Array_bi_over_d(1 To 1, 1 To 3), Arraybi_overdi(1 To 3) As Variant

Dim Nombre As String

Dim ArrayLR(0 To 4), ArrayLF(5 To 5), ArrayLS(6 To 11), ArrayXf(1 To 3) As Variant

Dim F, D, B, p, R, L, Lf, V, Num, Denom, Fila, SumyAlfa, ArrayNombre(1 To 3), SumB, SumBi, Sumd, M, HH As Variant

Dim i, j As Integer

Dim Arrayy(1 To 11, 1 To 3), Arrayx(1 To 11, 1 To 3), ArrayXD(1 To 3) As Variant

Dim ArrayT(0 To 10), ArrayAA(0 To 10, 1 To 3), ArrayZ(1 To 10, 1 To 3), ArrayAlfa(1 To 3) As Variant

Dim Arrayd(1 To 3), ArrayA(1 To 3), ArrayB(1 To 3), ArrayC(1 To 3), ArrayL(0 To 10), ArrayV(1 To 10), ArrayXB(1 To 3), ArrayXBN(1 To 3) As Variant

Sub Thiele_Geddes_Holland()

F = Worksheets(“Hoja1”).Cells(2, 3)

‘ArrayXB(1) = Worksheets(“Hoja1”).Cells(3, 3)

‘ArrayXB(2) = Worksheets(“Hoja1”).Cells(4, 3)

‘ArrayXB(3) = Worksheets(“Hoja1”).Cells(5, 3)

D = Worksheets(“Hoja1”).Cells(6, 3)

B = Worksheets(“Hoja1”).Cells(7, 3)

ArrayA(1) = Worksheets(“Hoja1”).Cells(8, 3)

ArrayB(1) = Worksheets(“Hoja1”).Cells(9, 3)

ArrayC(1) = Worksheets(“Hoja1”).Cells(10, 3)

ArrayA(2) = Worksheets(“Hoja1”).Cells(11, 3)

ArrayB(2) = Worksheets(“Hoja1”).Cells(12, 3)

ArrayC(2) = Worksheets(“Hoja1”).Cells(13, 3)

ArrayA(3) = Worksheets(“Hoja1”).Cells(14, 3)

ArrayB(3) = Worksheets(“Hoja1”).Cells(15, 3)

ArrayC(3) = Worksheets(“Hoja1”).Cells(16, 3)

p = Worksheets(“Hoja1”).Cells(17, 3)

R = Worksheets(“Hoja1”).Cells(18, 3)

L = Worksheets(“Hoja1”).Cells(19, 3)

Lf = Worksheets(“Hoja1”).Cells(20, 3)

V = Worksheets(“Hoja1”).Cells(21, 3)

ArrayXD(1) = Worksheets(“Hoja1”).Cells(22, 3)

ArrayXD(2) = Worksheets(“Hoja1”).Cells(23, 3)

ArrayXD(3) = Worksheets(“Hoja1”).Cells(24, 3)

‘F = Worksheets(“Hoja1”).Cells(25, 3)

ArrayXf(1) = Worksheets(“Hoja1”).Cells(26, 3)

ArrayXf(2) = Worksheets(“Hoja1”).Cells(27, 3)

ArrayXf(3) = Worksheets(“Hoja1”).Cells(28, 3)

‘Perfil de temperaturas, corrientes líquidas, y corrientes de vapor de la columna

‘Etapa Rectificación

ArrayT(0) = 373: ArrayLR(0) = 50

ArrayT(1) = 383: ArrayLR(1) = 50: ArrayV(1) = 100 ‘Etapa Rectificación

ArrayT(2) = 393: ArrayLR(2) = 50: ArrayV(2) = 100 ‘Etapa Rectificación

ArrayT(3) = 403.2: ArrayLR(3) = 50: ArrayV(3) = 100 ‘Etapa Rectificación

ArrayT(4) = 413.2: ArrayLR(4) = 50: ArrayV(4) = 100 ‘Etapa Rectificación

ArrayT(5) = 423.2: ArrayLF(5) = 150: ArrayV(5) = 100 ‘Etapa de alimentación

ArrayT(6) = 433.2: ArrayLS(6) = 150: ArrayV(6) = 100 ‘Etapa de Agotamiento

ArrayT(7) = 443.2: ArrayLS(7) = 150: ArrayV(7) = 100 ‘Etapa de Agotamiento

ArrayT(8) = 453.2: ArrayLS(8) = 150: ArrayV(8) = 100 ‘Etapa de Agotamiento

ArrayT(9) = 463.2: ArrayLS(9) = 150: ArrayV(9) = 100 ‘Etapa de Agotamiento

ArrayT(10) = 473.2: ArrayLS(10) = 150: ArrayV(10) = 100 ‘Etapa de Agotamiento

ArrayLS(11) = 150

‘ Encerado vectores

For j = 1 To 11

For i = 1 To 3

Arrayy(j, i) = 0

Arrayx(j, i) = 0

ArrayxSB(j, i) = 0

Next i

Next j

For j = 11 To 5 Step -1

For i = 1 To 3

ArrayyS(j, i) = 0

ArrayxS(j, i) = 0

Next i

Next j

For i = 1 To 3

ArrayxS(11, i) = 0

Arrayd(i) = 0

Next i

For j = 1 To 11

ArrayTc(j) = 0

ArrayTcn(j) = 0

Next j

‘Limpieza

Worksheets(“Hoja1”).Range(“a50:g100”).Clear

‘Cálculos de K’s y alfas

Fila = 50

Range(“a50:a53”).HorizontalAlignment = xlRight

Range(“b50:b53”).HorizontalAlignment = xlLeft

For j = 5 To 5

Worksheets(“Hoja1”).Cells(Fila, 1) = “Valores de alfa calculados a T(” & j & “)=”: Worksheets(“Hoja1”).Cells(Fila, 2) = ArrayT(j) & “oK”

Fila = Fila + 1

For i = 1 To 3

Denom = K(j, 3, p, ArrayT(), ArrayA(), ArrayB(), ArrayC())

ArrayAlfa(i) = K(j, i, p, ArrayT(), ArrayA(), ArrayB(), ArrayC()) / Denom

If i = 1 Then ArrayNombre(1) = “n-pentano”

If i = 2 Then ArrayNombre(2) = “n-exano”

If i = 3 Then ArrayNombre(3) = “n-decano”

Worksheets(“Hoja1”).Cells(Fila, 1) = “Alfa(” & i & “)=”: Worksheets(“Hoja1”).Cells(Fila, 2) = ArrayAlfa(i)

Worksheets(“Hoja1”).Cells(Fila, 3) = ArrayNombre(i)

Fila = Fila + 1

Next i

Next j

Fila = 54

Call Escritura_Titulos(Fila, ArrayNombre())

‘Calculos para zona de rectificación

‘Para la etapa 1 de la zona de rectificación,de tratamiento único

‘================================================================

SumyAlfa = 0: Fila = 56

For i = 1 To 3

Arrayy(1, i) = ArrayXD(i)

SumyAlfa = SumyAlfa + Arrayy(1, i) / ArrayAlfa(i)

Next i

For i = 1 To 3

Arrayx(1, i) = ((Arrayy(1, i) / ArrayAlfa(i))) / SumyAlfa

Next i

ArrayTc(1) = 373.2 ‘oK

Call Calculo_Temperaturas(ArrayTc(), ArrayTcn(), ArrayA(), ArrayB(), ArrayC(), Arrayx(), 1, p)

Call Escritura_Resultados(Arrayx(), Arrayy(), ArrayTcn(), Fila, 1)

Fila = Fila + 1

‘Cálculo para las etapas 2 a 5 de la zona de rectificación

‘=========================================================

For j = 2 To 5

SumAlfa = 0

 

For i = 1 To 3

 

Arrayy(j, i) = ((ArrayLR(j – 1) * Arrayx(j – 1, i)) / ArrayV(j)) + (D * ArrayXD(i) / ArrayV(j))

SumAlfa = SumAlfa + (Arrayy(j, i) / ArrayAlfa(i))

Next i

For i = 1 To 3

Arrayx(j, i) = (Arrayy(j, i) / ArrayAlfa(i)) / SumAlfa

Next i

ArrayTc(j) = 473.2 ‘oK

Call Calculo_Temperaturas(ArrayTc(), ArrayTcn(), ArrayA(), ArrayB(), ArrayC(), Arrayx(), j, p)

Call Escritura_Resultados(Arrayx(), Arrayy(), ArrayTcn(), Fila, j)

Fila = Fila + 1

Next j

‘Cálculos de los fondos (Etapa 11)

‘=================================

SumB = 0

 

For j = 11 To 11

For i = 1 To 3

ArrayxSB(j, i) = (F * ArrayXf(i) – D * ArrayXD(i)) ‘ArrayxSB=B*XBi

SumB = SumB + ArrayxSB(j, i)

Next i

For i = 1 To 3

ArrayxS(j, i) = ArrayxSB(j, i) / SumB

Next i

SumB = 0

For i = 1 To 3

SumB = SumB + ArrayAlfa(i) * ArrayxS(j, i)

Next i

For i = 1 To 3

ArrayyS(j, i) = (ArrayAlfa(i) * ArrayxS(j, i)) / SumB

Next i

ArrayTc(j) = 473.2 ‘oK

Call Calculo_Temperaturas(ArrayTc(), ArrayTcn(), ArrayA(), ArrayB(), ArrayC(), ArrayxS(), j, p)

Call Escritura_ResultadosS(ArrayxS(), ArrayyS(), 67, j)

Next j

‘Cálculos de las concentraciones en la sección de agotamiento

‘============================================================

Fila = 66

For j = 10 To 5 Step -1

For i = 1 To 3

If j <> 5 Then

ArrayxS(j, i) = ((ArrayV(j) * ArrayyS(j + 1, i)) / ArrayLS(j)) + ((ArrayxSB(11, i) / ArrayLS(j)))

ElseIf j = 5 Then

ArrayxS(j, i) = (((ArrayV(j) * ArrayyS(j + 1, i)) / ArrayLF(j))) + ((ArrayxSB(11, i) / ArrayLF(j)))

End If

SumBi = SumBi + ArrayAlfa(i) * ArrayxS(j, i)

Next i

For i = 1 To 3

ArrayyS(j, i) = ((ArrayAlfa(i) * ArrayxS(j, i)) / SumBi)

Next i

ArrayTc(j) = 473.2 ‘oK

Call Calculo_Temperaturas(ArrayTc(), ArrayTcn(), ArrayA(), ArrayB(), ArrayC(), ArrayxS(), j, p)

Call Escritura_ResultadosS(ArrayxS(), ArrayyS(), Fila, j)

Fila = Fila – 1

Next j

 

‘Cálculo de los valores de XD(i) para la siguiente iteración

‘==================================================

Sumd = 0: SumB = 0

For i = 1 To 3

Arrayd(i) = (F * ArrayXf(i)) / (1 + (((Arrayy(5, i) * (ArrayxSB(11, i)) / (D * ArrayXD(i) * (ArrayyS(5, i)))))))

Sumd = Sumd + Arrayd(i)

Next i

For i = 1 To 3

ArrayXDN(i) = Arrayd(i) / Sumd

Next i

Call Escritura_de_Nuevas_XDN(68, ArrayXDN())

For j = 1 To 11

ArrayTc(j) = 100

Next j

For j = 1 To 11

Next j

End Sub

Function K(j, i, p, ArrayT(), ArrayA(), ArrayB(), ArrayC())

K = 10 ^ ((ArrayA(i) – ((ArrayB(i) / (ArrayT(j) + ArrayC(i))))) / p)

End Function

Sub Escritura_Titulos(Fila, ArrayNombre())

Range(“a53:g67”).HorizontalAlignment = xlCenter

Worksheets(“Hoja1”).Cells(Fila, 1) = “Etapa”

Worksheets(“Hoja1”).Cells(Fila, 2) = ArrayNombre(1)

Worksheets(“Hoja1”).Cells(Fila, 3) = ArrayNombre(2)

Worksheets(“Hoja1”).Cells(Fila, 4) = ArrayNombre(3)

Worksheets(“Hoja1”).Cells(Fila, 5) = ArrayNombre(1)

Worksheets(“Hoja1”).Cells(Fila, 6) = ArrayNombre(2)

Worksheets(“Hoja1”).Cells(Fila, 7) = ArrayNombre(3)

Fila = Fila + 1

Worksheets(“Hoja1”).Cells(Fila, 1) = “j”

Worksheets(“Hoja1”).Cells(Fila, 2) = “y(j,1)”

Worksheets(“Hoja1”).Cells(Fila, 3) = “y(j,2)”

Worksheets(“Hoja1”).Cells(Fila, 4) = “y(j,3)”

Worksheets(“Hoja1”).Cells(Fila, 5) = “x(j,1)”

Worksheets(“Hoja1”).Cells(Fila, 6) = “x(j,2)”

Worksheets(“Hoja1”).Cells(Fila, 7) = “x(j,3)”

Worksheets(“Hoja1”).Cells(Fila, 8) = “Temp,oK”

End Sub

Sub Escritura_Resultados(Arrayx(), Arrayy(), ArrayTcn(), Fila, j)

Worksheets(“Hoja1”).Cells(Fila, 1) = j

Worksheets(“Hoja1”).Cells(Fila, 2) = Arrayy(j, 1)

Worksheets(“Hoja1”).Cells(Fila, 3) = Arrayy(j, 2)

Worksheets(“Hoja1”).Cells(Fila, 4) = Arrayy(j, 3)

Worksheets(“Hoja1”).Cells(Fila, 5) = Arrayx(j, 1)

Worksheets(“Hoja1”).Cells(Fila, 6) = Arrayx(j, 2)

Worksheets(“Hoja1”).Cells(Fila, 7) = Arrayx(j, 3)

Worksheets(“Hoja1”).Cells(Fila, 8) = ArrayTcn(j)

End Sub

Sub Escritura_de_Nuevas_XDN(Fila, ArrayXDN())

Worksheets(“Hoja1”).Cells(Fila, 1) = “NUEVOS VALORES DE XD(I)”

Fila = Fila + 1

For i = 1 To 3

Worksheets(“Hoja1”).Cells(Fila, 1) = “XDN(” & i & “)=”: Worksheets(“Hoja1”).Cells(Fila, 2) = ArrayXDN(i)

Fila = Fila + 1

Next i

End Sub

Sub Escritura_ResultadosS(ArrayxS(), ArrayyS(), Fila, j)

Worksheets(“Hoja1”).Cells(Fila, 1) = j

Worksheets(“Hoja1”).Cells(Fila, 2) = ArrayyS(j, 1)

Worksheets(“Hoja1”).Cells(Fila, 3) = ArrayyS(j, 2)

Worksheets(“Hoja1”).Cells(Fila, 4) = ArrayyS(j, 3)

Worksheets(“Hoja1”).Cells(Fila, 5) = ArrayxS(j, 1)

Worksheets(“Hoja1”).Cells(Fila, 6) = ArrayxS(j, 2)

Worksheets(“Hoja1”).Cells(Fila, 7) = ArrayxS(j, 3)

Worksheets(“Hoja1”).Cells(Fila, 8) = ArrayTcn(j)

End Sub

Function Fu(ArrayTc(), ArrrayA(), ArrayB(), ArrayC(), Arrayx(), j, p)

Fu = 0

For i = 1 To 3

Fu = Fu + Arrayx(j, i) * 10 ^ (ArrayA(i) – ((ArrayB(i) / (ArrayTc(j) + ArrayC(i)))))

Next i

Fu = Fu – p

End Function

Function Fp(ArrayTc(), ArrrayA(), ArrayB(), ArrayC(), Arrayx(), j)

Fp = 0

For i = 1 To 3

HH = 10 ^ (ArrayA(i) – (ArrayB(i) / (ArrayTc(j) + ArrayC(i))))

Fp = Fp + Arrayx(j, i) * HH * Log(10) * (ArrayB(i) / (ArrayTc(j) – ArrayC(i)) ^ 2)

Next i

End Function

Sub Calculo_Temperaturas(ArrayTc(), ArrayTcn(), ArrayA(), ArrayB(), ArrayC(), Arrayx(), j, p)

‘ArrayTcn(j) = 0

For M = 1 To 30

ArrayTcn(j) = ArrayTc(j) – (Fu(ArrayTc(), ArrayA(), ArrayB(), ArrayC(), Arrayx(), j, p) / Fp(ArrayTc(), ArrayA(), ArrayB(), ArrayC(), Arrayx(), j))

If Abs(ArrayTcn(j) – ArrayTc(j)) < 0.001 Then

Exit For

Else

ArrayTc(j) = ArrayTcn(j)

End If

Next M

End Sub

Breve explicación del programa

La primera cosa que se debe decir es que el macro se intitula Thiele Gedddes, porque ese fue el algoritmo que se intentó presentar, pero luego se decidió implementar el algoritmo de Lewis y Matheson, que es el que en realidad se presenta.

El segundo aspecto que se debe aclarar es que el programa que aquí se muestra tiene dos funciones y una subrutina que permiten el cálculo de la temperatura de punto de burbuja de cada plato; y que, por la cercanía de algunos valores, de temperatura de la sección de agotamiento, pienso que tal vez esta sección debió haber tenido menos de cinco platos (véase la Figura 1).

Por lo demás, es pertinente agregar que las composiciones son “Arrays”, que es lo que (Holland, 1963) implica en su libro.

Por la bibliografía que se revisó entiendo que el Profesor Holland publicó una edición posterior, de la cual no dispongo, por lo que hube de construir el algoritmo a partir del ejemplo tabular que se presenta en el libro de la referencia.

Las composiciones de la sección de enriquecimiento se manejan por medio de dos matrices x(j,i), y y(j,i). Para la sección de agotamiento (stripping) se manejan  las matrices xS(i,j) e yS(i,j).

El programa calcula los coeficientes de volatilidad alfa a partir de los valores de Ai, B y Ci de la Ecuación de Antoine, de la manera que ya se indicó.

Los lazos anidados de cálculo no son más que lo que (Holland, 1963) presenta en su libro en forma tabular, porque en esos tiempos recién nacía el IBM 1130, que después fue reemplazado por el modelo 360, y -aunque no lo crean- tampoco había ni Lotus, el Excel tardaría todavía veinte o más años en hacer su aparición y ni se pensaba en los computadores personales.

Referencia

Holland, C. D. (1963). Multicomponent Distillation. Englewoods Cliffs, New Jersey: Prentice-Hall, Inc.

 

 

 

Publicado en DISEÑO DE PLANTAS, MODELADO Y SIMULACIÓN DE PROCESOS

El Método de Thiele Geddes para Cálculo de Columnas de Destilación de Multicomponentes

En este método se cuentan las etapas de equilibrio desde el tambor de reflujo hacia abajo, en la columna, tanto en la sección de enriquecimiento, como en la sección (Holland, 1963) de agotamiento En una primera aproximación el método es tabular, y esa parte puede llevarse a cabo en una hoja Excel™, aunque es laboriosa y hay que poner mucha atención para no cometer errores.

Esta parte “tabular” se encuentra en las páginas 56-58 de la referencia bibliográfica. Voy a tratar de reproducirla en este artículo, aunque por su extensión, tendré que hacerlo por fragmentos.

Figura1

Tabla 1

Figura2

Tabla 2

Figura3

Tabla 4

Figura5

Tabla 5

 

 

A partir de las tablas mostradas (Holland, 1963) recomienda embarcarse en un refinamiento tabular de los coeficientes de reparto K, cosa que yo no hice, porque me pareció demasiado laborioso, y demasiado “manual.

En vez de eso me adelanté al Capítulo 4 de la referencia, a la sección intitulada Development of the Methos of Convergence for Use with the Calculational Procedure of Thiele and Geddes (Holland, 1963) en la que el autor recomienda el método de Newton-Raphson para determinar un multiplicador que él denomina Θ, que sirve para corregir las aproximaciones de la tabulación.

Para los fines consiguientes (Holland, 1963) formula la siguiente expresión,

Ecuacion1

Ecuación 1

Y la derivada de la expresión anterior, que se muestra a continuación

Ecuacion2

Ecuación 2

En ambas expresiones el subíndice “ca” significa el valor calculado en la tabla, y el Θ es el factor de corrección, que se utiliza de la forma que se ilustra a continuación, para obtener los valores corregidos de .

La expresión del valor corregido de los valores de di es la que se muestra a continuación.

Picture3

Ecuación 3

La expresión para corregir los valores de bi se muestra a continuación.

Ecuacion4

Ecuación 4

Las fracciones molares corregidas pueden calcularse, de acuerdo a lo que consta en la página 85 de (Holland, 1963), mediante los siguientes algoritmos.

Ecuacion5

Ecuación 5

 

y

 

Ecuacion6

Ecuación 6

Lo anterior lo implementé, en forma de un algoritmo de Newton Raphson, que el lo que recomienda (Holland, 1963) por medio de un programa que muestro a continuación, conjuntamente con los resultados para las xji. El cálculo de las yji se lo dejo a algún amable lector que tenga el espíritu de trabajo de hacerlo.

El programa

Option Explicit

Option Base 1

Dim D, Theta, Thetan, Sumdc, Sumbc, Sumx As Variant

Dim k, N, Fila, j, i As Integer

Dim Arrayl(1 To 3, 1 To 3), Arrayv(1 To 3, 1 To 3), Arrayx(1 To 3, 1 To 3) As Variant

Dim ArrayFX(1 To 4), Arrayb(1 To 3), Arrayd(1 To 3), ArrayCorr(1 To 3), Arraybc(1 To 3), Arraydn(1 To 3), Arraydc(1 To 3) As Variant

Sub Direct_Iteration()

ArrayFX(1) = Worksheets(“Hoja2”).Cells(2, 3)

ArrayFX(2) = Worksheets(“Hoja2”).Cells(3, 3)

ArrayFX(3) = Worksheets(“Hoja2”).Cells(4, 3)

Arrayb(1) = Worksheets(“Hoja2”).Cells(5, 3)

Arrayb(2) = Worksheets(“Hoja2”).Cells(6, 3)

Arrayb(3) = Worksheets(“Hoja2”).Cells(7, 3)

Arrayd(1) = Worksheets(“Hoja2”).Cells(8, 3)

Arrayd(2) = Worksheets(“Hoja2”).Cells(9, 3)

Arrayd(3) = Worksheets(“Hoja2”).Cells(10, 3)

D = Worksheets(“Hoja2”).Cells(11, 3)

Theta = Worksheets(“Hoja2”).Cells(12, 3)

Arrayl(1, 1) = Worksheets(“Hoja2”).Cells(13, 3)

Arrayl(1, 2) = Worksheets(“Hoja2”).Cells(14, 3)

Arrayl(1, 3) = Worksheets(“Hoja2”).Cells(15, 3)

Arrayl(2, 1) = Worksheets(“Hoja2”).Cells(16, 3)

Arrayl(2, 2) = Worksheets(“Hoja2”).Cells(17, 3)

Arrayl(2, 3) = Worksheets(“Hoja2”).Cells(18, 3)

Arrayl(3, 1) = Worksheets(“Hoja2”).Cells(19, 3)

Arrayl(3, 2) = Worksheets(“Hoja2”).Cells(20, 3)

Arrayl(3, 3) = Worksheets(“Hoja2”).Cells(21, 3)

 

Arrayv(1, 1) = Worksheets(“Hoja2”).Cells(22, 3)

Arrayv(1, 2) = Worksheets(“Hoja2”).Cells(23, 3)

Arrayv(1, 3) = Worksheets(“Hoja2”).Cells(24, 3)

Arrayv(2, 1) = Worksheets(“Hoja2”).Cells(25, 3)

Arrayv(2, 2) = Worksheets(“Hoja2”).Cells(26, 3)

Arrayv(2, 3) = Worksheets(“Hoja2”).Cells(27, 3)

Arrayv(3, 1) = Worksheets(“Hoja2”).Cells(28, 3)

Arrayv(3, 2) = Worksheets(“Hoja2”).Cells(29, 3)

Arrayv(3, 3) = Worksheets(“Hoja2”).Cells(30, 3)

 

‘Limpieza

Worksheets(“Hoja2”).Range(“a40:d53”).Clear

For N = 1 To 10

Thetan = Theta – (g(Theta, k, D, ArrayFX(), Arrayb(), Arrayd()) / (gprima(Theta, k, ArrayFX(), Arrayb(), Arrayd())))

If Abs(Thetan – Theta) < 0.001 Then

Sumdc = 0: Sumbc = 0

For k = 1 To 3

Arraydc(k) = (ArrayFX(k) / (1 + Thetan * (Arrayb(k) / Arrayd(k))))

Sumdc = Sumdc + Arraydc(k)

Arraybc(k) = Thetan * ((Arrayb(k) / Arrayd(k)) * Arraydc(k))

Sumbc = Sumbc + Arraybc(k)

Worksheets(“Hoja2”).Cells(40 + k, 1) = “dc(” & k & “)=” & Round(Arraydc(k), 4)

Worksheets(“Hoja2”).Cells(40 + k, 2) = “bc(” & k & “)=” & Round(Arraybc(k), 4)

Next k

Worksheets(“Hoja2”).Cells(40 + k, 1) = Round(Sumdc, 4): Worksheets(“Hoja2”).Cells(40 + k, 2) = Round(Sumbc, 4)

Exit For

Else

Theta = Thetan

End If

Next N

‘Escritura de las fracciones molares

Fila = 47

Columns(“A:D”).HorizontalAlignment = xlCenter

Worksheets(“Hoja2″).Cells(Fila – 2, 1) = ” j”

For i = 1 To 3

Worksheets(“Hoja2”).Cells(Fila – 2, i + 1) = i

Next i

For i = 1 To 3

Worksheets(“Hoja2”).Cells(Fila – 1, i + 1) = Arraydc(i) / Sumdc

Next i

 

 

Worksheets(“Hoja2”).Cells(Fila – 1, 1) = “0”

For j = 1 To 3

 

For i = 1 To 3

Sumx = Sumx + ((Arrayl(j, i) / Arrayd(i))) * Arraydc(i)

Next i

Worksheets(“Hoja2”).Cells(Fila, 1) = j

For i = 1 To 3

 

Arrayx(j, i) = (((Arrayl(j, i) / Arrayd(i))) * Arraydc(i)) / Sumx

 

Worksheets(“Hoja2”).Cells(Fila, i + 1) = Arrayx(j, i)

Next i

 

 

Sumx = 0: Fila = Fila + 1

Next j

 

Worksheets(“Hoja2”).Cells(Fila, 1) = 4

For i = 1 To 3

Worksheets(“Hoja2”).Cells(Fila, i + 1) = Arraybc(i) / Sumbc

Next i

 

 

End Sub

Function g(Theta, k, D, ArrayFX(), Arrayb(), Arrayd())

g = 0

For k = 1 To 3

g = g + ArrayFX(k) / (1 + Theta * (Arrayb(k) / Arrayd(k)))

Next k

g = g – D

End Function

Function gprima(Theta, k, ArrayFX(), Arrayb(), Arrayd())

gprima = 0

For k = 1 To 3

gprima = gprima – ((Arrayb(k) / Arrayd(k)) * ArrayFX(k)) / (1 + Theta * (Arrayb(k) / Arrayd(k))) ^ 2

Next k

End Function

Los resultados

Figura 6

Tabla 6. Resultados

Comentarios a los resultados

Los resultados son coherentes ya que el compuesto 1 es el menos volátil y el compuesto 3 es el más volátil.

Trabajo Citado

Holland, C. D. (1963). Multicomponent Distillation (Vol. 1). (I. Prentice Hall, Ed.) Englewood Cliffs, New Jersey, U.S.A.: Prentice Hall, Inc.

 

Publicado en DISEÑO DE PLANTAS, MODELADO Y SIMULACIÓN DE PROCESOS

Enter your email address to follow this blog and receive notifications of new posts by email.

Únete a otros 1.656 seguidores

Categorías
Artículos y comentarios sobre modelado y simulación de plantas, y equipos de la industria química, con ejemplos
PREGUNTAS O INQUIETUDES
PUEDEN ENVIAR SUS PREGUNTAS/INQUIETUDES RESPECTO DE ARTÍCULOS DEL BLOG, O TEMAS AFINES A gasteaux@hotmail.com ASUNTO: BLOG DE INGENIERÍA QUÍMICA