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.

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

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:

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.

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

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:

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.

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 F_A denota el flujo volumétrico de la especie A, que F_B denota el flujo volumétrico de B, y que F_C 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.

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

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 C_A=C_Ao, C_B=C_Bo, C_C=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.

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 F_A, F_B, 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.

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 m³/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 h_set, 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.

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.

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.

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

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

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

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).

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.

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.

É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

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 x_1,i es la siguiente:

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.

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

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.

y

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

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

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.

Los 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.

Objetivo del artículo

El objetivo de este artículo es mostrar cómo se debe calcular la composición de un flash adiabático, y exponer el papel que juegan el punto de burbuja y el punto de rocío en el cálculo del Flash.

El artículo está inspirado en unas hojas del libro del Profesor Skogestad, de la Politécnica de Noruega a juzgar por el sitio web (http://folk.ntnu.no/skoge/bok/mer/flash_english_edition_2009), que está escrito en un lenguaje sencillo y muy entendible.

Restricciones al modelo matemático

Se ha tratado a los componentes de la mezcla mediante la Ley de Raoult; se ha estimado sus presiones de vapor mediante el modelo de Antoine, .

Donde

Mezcla que se alimenta al tambor Flash

A continuación se indica la mezcla que se utilizó, así como los respectivos coeficientes de Antoine, que se encuentran en la fuente citada.

Figura 1. Datos del Problema

Para estos datos se escribió un pequeño programa, que se transcribe más adelante, que maneja tres subrutinas, quq calculan el punto de burbuja, el punto de rocío , y la composición del Flash.

Resultados

Los resultados del programa se muestran a continuación.

Figura 2. Resultados del programa

Cálculo del punto de burbuja

En este caso se calculó el punto de burbuja a una presión dada, para lo que se trabajó con la ecuación siguiente (Skogestad);

La solución del problema se reduce a encontrar la temperatura que satisfaga la ecuación, lo que se realiza mediante un For-Next  y una condición de convergencia. También se muestra en el programa, y en la Figura 1, la composición de la primera burbuja que se forma.

Cálculo del punto de rocío

Para este propósito se considera que la mezcla de alimentación está constituida de vapor, y se busca la temperatura a la que se condensa la primera gota de líquido, lo que significa que en la composición de la mezcla que está expresada en fracciones molares Z_i, éstas deben expresarse como y_i,

Lo anterior significa que, para fines de cálculo se utiliza la expresión (Skogestad) siguiente:

El cálculo, dada la presión, implica encontrar una temperatura que satisfaga la ecuación anterior de acuerdo a un criterio de convergencia que se muestra en el programa. En el programa también se muestran las fracciones molares de la primera gota que se condensa.

Cálculo del Flash

Para los fines consiguientes se utilizó la ecuación siguiente con un criterio de convergencia adecuado que se puede observar en el programa, y que obtiene el V/F. Los cálculos de la composición pueden observarse en el programa.

El programa

Option Base 1

Option Explicit

Dim ArrayZ(1 To 5), ArrayA(1 To 5), ArrayB(1 To 5), ArrayC(1 To 5), ArrayY(1 To 5), Arrayp(1 To 5), p As Variant

Dim ArrayK(1 To 5), ArrayX(1 To 5) As Variant

Dim Summx As Variant

Dim kk, Nc, Fila, N, PP As Integer

Dim SumFPB, SumFPPB, T, Tn, Tv, M, MM, MMM, MMMM, Sump, j, Psiv, Psin, FsumFPsi, F2Psi, Sumx, Psi, Sumy, SumPsi, SumFPsi, SumFLV, SumPFLV As Variant

Sub destilacion_flash()

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

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

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

ArrayZ(4) = Worksheets(“Hoja1”).Cells(5, 3)

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

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

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

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

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

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

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

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

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

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

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

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

Log (2) / Log(10)

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

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

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

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

T = T + 273.2

Tv = T

N = 0: Sump = 0

‘Limpieza

Worksheets(“Hoja1”).Range(“A25:g33000”).Clear

Fila = 25

Call Punto_de_Burbuja(p, ArrayA(), ArrayC(), ArrayZ(), ArrayB(), Tv, Tn, Fila)

Call Dew_Point(31, ArrayA(), ArrayB(), ArrayC(), ArrayZ(), ArrayY(), 5)

Call Tanteo(ArrayZ(), ArrayK(), ArrayA(), ArrayB(), ArrayC(), ArrayX(), 4, 375, 37)

End Sub

Sub Punto_de_Burbuja(p, ArrayA(), ArrayC(), ArrayZ(), ArrayB(), Tv, Tn, Fila)

SumFPPB = 0

For T = 354 To 355 Step 0.1

SumFPB = 0:

For kk = 1 To Nc

MM = ArrayZ(kk) * (10 ^ (ArrayA(kk) – (ArrayB(kk) / (T + ArrayC(kk)))))

SumFPB = SumFPB + MM

Next kk

If Abs(SumFPB – p) <= 0.05 Then

N = N + 1

Worksheets(“Hoja1”).Cells(Fila, 1) = “El Punto de Burbuja es  ” & T & ” oK” & ” ó ” & “T= ” & (T – 273.2) & ” oC ”

Exit For

End If

Next T

Sump = 0

For kk = 1 To Nc

Arrayp(kk) = ArrayZ(kk) * (10 ^ (ArrayA(kk) – (ArrayB(kk) / (T + ArrayC(kk)))))

Worksheets(“Hoja1″).Cells(Fila + kk, 1) = ” y(” & kk & “)= “: Worksheets(“Hoja1”).Cells(Fila + kk, 2) = Arrayp(kk) / p

Sump = Sump + (Arrayp(kk) / p)

‘ArrayY(kk) = Arrayp(kk) / p

Next kk

Worksheets(“Hoja1”).Cells(Fila + kk, 2) = Sump

Fila = Fila + 6

Exit Sub

Fila = Fila + 1

End Sub

Sub Tanteo(ArrayZ(), ArrayK(), ArrayA(), ArrayB(), ArrayC(), ArraX(), p, T, Fila)

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

For Psi = 0.02 To 1 Step 0.02

Sumx = 0: Sumy = 0

For kk = 1 To Nc

ArrayK(kk) = (10 ^ (ArrayA(kk) – ((ArrayB(kk) / (T + ArrayC(kk)))))) / p

ArrayX(kk) = ArrayZ(kk) / (1 + Psi * (ArrayK(kk) – 1))

ArrayY(kk) = ArrayX(kk) * ArrayK(kk)

Sumy = Sumy + ArrayY(kk)

Sumx = Sumx + ArrayX(kk)

Next kk

If Abs(Sumx – 1) < 0.001 Then Exit For

Next Psi

Worksheets(“Hoja1”).Cells(Fila + 1, 1) = “V/F=” & Psi: Worksheets(“Hoja1”).Cells(Fila + 1, 2) = “T= ” & T & ” oK o ” & T – 273.2 & “oC”

Worksheets(“Hoja1”).Cells(Fila + 1, 3) = “p=” & p

For kk = 1 To Nc

Worksheets(“Hoja1”).Range(“A43:B43”).HorizontalAlignment = xlLeft

Worksheets(“Hoja1”).Cells(Fila + 1 + kk, 1) = “x(” & kk & “)= ” & ArrayX(kk): Worksheets(“Hoja1”).Cells(Fila + 1 + kk, 2) = “y(” & kk & “)= ” & ArrayY(kk)

Next kk

Worksheets(“Hoja1”).Cells(Fila + 1 + kk, 1) = Sumx: Worksheets(“Hoja1”).Cells(Fila + 1 + kk, 2) = Sumy

End Sub

Sub Dew_Point(Fila, ArrayA(), ArrayB(), ArrayC(), ArrayZ(), Arrayp(), p)

For T = 354.3 To 410

Sump = 0: Summx = 0

For kk = 1 To Nc

Sump = Sump + ArrayZ(kk) / 10 ^ (ArrayA(kk) – (ArrayB(kk) / (T + ArrayC(kk))))

Next kk

If Abs(Sump – (1 / p)) < 0.005 Then

Worksheets(“Hoja1”).Cells(Fila, 1) = “El Punto de Rocío es T= ” & T & “oK” & ” ó ” & T – 273.2 & “oC”

Exit For

End If

Next T

For kk = 1 To Nc

Worksheets(“Hoja1”).Cells(Fila + kk, 2) = “x(” & kk & “)=” & (ArrayZ(kk) * p) / 10 ^ (ArrayA(kk) – (ArrayB(kk) / (T + ArrayC(kk))))

Summx = Summx + (ArrayZ(kk) * p) / 10 ^ (ArrayA(kk) – (ArrayB(kk) / (T + ArrayC(kk))))

Next kk

Worksheets(“Hoja1”).Range(“b36”).HorizontalAlignment = xlLeft

Worksheets(“Hoja1”).Cells(Fila + kk, 2) = Summx

End Sub

Un poco de filosofía

Cuando uno adquiere el compromiso de escribir un blog, también adquiere el de mantenerlo, y a veces éste último es el más difícil de cumplir. Por las estadísticas de visitas al aitio Web, el tema de diseño de plantas es lo que más gusta, pero es el tema más difícil de abordar porque toma un tiempo muy grande.

Con respecto al tema del diseño de plantas, el amable lector Ingeniero Químico debe reconocer que para diseñar una planta (en realidad, un instalación industrial constituida de varias plantas) hay que primero saber diseñar los equipos, y que para diseñarlos hay que tener muy claros los conceptos de los fenómenos de transferencia, y de la termodinámica.

Por eso se me ocurrió el tema que queda señalado, para lo que hube de ponerme a mirar los conceptos que les menciono con mucha claridad, que en mi formación (del siglo pasado), no los recibí ni en pregrado, ni en postgrado.

El tema y su tratamiento

La primera cosa que me ocurrió fue que había que elegir un problema. Y para eso busqué primero un algoritmo, en el libro Smith, que usé en mi post grado, pero la verdad es que me desilusioné, porque no es un libro algorítmico, porque es de la época en que el IBM 1130, que tenía menos memoria que mi iPhone.

Entonces me acoré de un libro que, en su tiempo me pareció complicado, el famoso Multicomponent Distillation, de (Holland), que sí es algorítmico, y por lo tanto hecho para aplicar los métodos que allí se exponen en la computadora, aunque coincidentemente hayan sido publicados, los dos, en 1962.

Leyendo a Holland me enteré que él se había basado en el famoso libro de Brown, para tratar el tema, y entonces me dije que quizá podría ser de interés para mis lectores, y por eso lo abordé.

El concepto

Para Holland –y para todos, por lo menos hasta que nos damos cuenta (y hablo por mi), la destilación flash es una destilación multicomponente de una sola etapa. ¿Verdad de Perogrullo?: Si y no, porque una vez que se lee la afirmación de Holland, es obvia, y antes no lo es.

¿Genialidad de Holland de presentar el tema de los multicomponentes en forma algorítmica?, o sabía que venía el IBM 360, que haría aplicable todo lo que él propone en su libro? Nunca lo sabremos. La verdad es que en mi postgrado trabajábamos esos problemas con unos calculadores carísimos a tubos de vacío, de cuatro operaciones.

Más conceptos

Uno debe aproximarse A las mezclas multicomponentes en la dirección que va desde el líquido hacia el vapor si de una expansión flash se trata.

Si así se procede entonces lo que primero se encuentra es el punto de burbuja, que es la temperatura a la que se forma la primera burbuja (realmente es el comienzo de la ebullición). Pero como se trata de una mezcla no se empieza a evaporar hasta que no se llega al punto de rocío, que es la temperatura a la que se comienza a formar vapor. Entre estas dos temperaturas hay una zona de transición.

El punto de burbuja

Aquí yo voy a utilizar los subíndices que usé en mi programa, para evitar confusiones.

La primera ecuación es una verdad termodinámica más grande que una casa, que lo único que dice es que la suma de las fracciones molares en fase líquida multiplicada por el coeficiente de reparto K(zz), que es igual a y(kk)/x(kk).

En el caso que nos ocupa, y habiendo sido Holland profesor del Departamento de Ingeniería Química del Agricultural and Mechanical College of Texas, el ejemplo, que él tomó de (Brown) corresponde a una mezcla de hidrocarburos n-C₂, n-C₃, i-C₃, n-C₄, i-C₄, n-C₅, i-C₅, y n-C₆.

En esos tiempos se trabajaba en psia, y las tablas de Holland, que él toma del Natural Gasoline Asssociation of América, se presentan en forma de expansiones polinomiales de K(zz)=f(T), a 50 psia, donde la temperatura está en grados Rankine.

En este punto me di cuenta que las tablas de Holland que yo estaba mirando no abarcaban los compuestos iso y entonces procedí a modificar la mezcla sumando las fracciones molares de los hidrocarburos normales  a las de los iso y tomando esa suma como si se tratase de sólo de hidrocarburos normales.

Después me di cuenta que algunas combinaciones de coeficientes estaban bien, porque obtenían los valores de Brown y otras no, pero ya estaba metido en el baile y debí hacerles algunas “correcciones” (educated guesses), hasta que obtuve números de K(kk) que hacían lógica, ya que mientras más pesado es el hidrocarburo menor es su volatilidad, y menor es su K(kk).

De todo esto resultó que de siete hidrocarburos que aparecen en el problema de Brown yo me reduje a cinco, como les indico a continuación.

Figura 1. Valores de K para la estimación de Punto de burbuja de la alimentación

La condición para poder calcular el punto de burbuja, que se convierte en ecuación es la siguiente.

El problema se reduce, entonces a encontrar la temperatura que satisfaga la ecuación anterior, que –escrita como función- queda como se muestra a continuación, y se conoce como la curva de temperatura de punto de burbuja de la mezcla de alimentación al tambor flash.

Esto significa, que la función puede graficarse (de hecho Hollan lo muestra en su página 16, Cáp.2, donde manifiesta que en vez del método numérico de Newton-Raphson se debe usar el método de Regula Falsi, que funciona con dos aproximaciones.

El gráfico obtenido fue el siguiente.

Figura 2. Trazo de la Función de Punto de Burbuja

Lo interesante de todo esto, por lo menos para mí, es que la discusión acerca de qué método numérico usar para determinar el Punto de Burbuja de la mezcla no es importante porque cualquiera puede leer el  gráfico anterior, y saber ipso-facto, que el punto de burbuja es 570^oR “y piquillo”.

Sin embargo de lo anteriormente expuesto, Holland anota en la página 15, que es posible obtener el punto de burbuja usando el método de Newton-Raphson.

Y, claro, después de todas las peripecias que he narrado, yo también apliqué el método anterior, usando como primera aproximación 400^oR, obteniendo, en ¡98 iteraciones!, en la primera versión de este artículo tenía yo una equivocación en la derivada y en la subrutina de cáculo de punto de burbuja, por eso el cambio de texto.

Pto.Burbuja570.497733705646oRL=98

Las concentraciones que usé para la alimentación fueron las siguientes:

Figura 3. Concentraciones de alimentación

Los coeficientes polinomiales que mencioné antes fueron los siguientes:

Figura 4. Coeficientes polinomiales para el vaso de K de cada hidrocarburo

El programa

Option Base 1

Option Explicit

Dim Arrayx(1 To 10) As Variant

Dim ArrayA(1 To 10, 1 To 10) As Variant

Dim ArrayK(1 To 20) As Variant

Dim ArrayKp(1 To 20) As Variant

Dim i, j, NC, kk, L, Fila As Integer

Dim T As Double

Dim Sum, Sumr, Tv, Tn, MM, MMM, MMMM, Sump As Variant

Dim De_Nuevo As String

Dim HH As Double

Sub DESTILACION()

Fila = 36: Sum = 0: Sumr = 0

‘Limpieza

Worksheets(“Hoja1”).Range(“A36:f300”).Clear

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

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

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

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

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

‘x(6) = Worksheets(“Hoja1”).Cells(7, 3)

‘x(7) = Worksheets(“Hoja1”).Cells(8, 3)

NC = 5

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ArrayA(4, 5) = Worksheets(“Hoja1”).Cells(33, 3)

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

T = T + 459.67

For kk = 1 To NC

ArrayK(kk) = 0

Next kk

Call Calculo_de_K(T, ArrayA(), ArrayK(), NC)

Call Imprimir_Resultados(Fila, ArrayK(), kk, NC)

Call Calculo_Funcion_Temperatura_Burbuja_Rocio(ArrayK(), Arrayx(), T, Fila, Sum, Sumr)

Tv = 700

Call Calculo_del_Punto_Burbuja(Tv, Tn, ArrayA(), ArrayK(), ArrayKp(), Fila, Arrayx())

End Sub

Sub Calculo_de_K(T, ArrayA(), ArrayK(), NC)

For kk = 1 To NC

ArrayK(kk) = T * (ArrayA(1, kk) + ArrayA(2, kk) * T + ArrayA(3, kk) * T ^ 2 + ArrayA(4, kk) * T ^ 3) ^ 3

Next kk

End Sub

Sub Imprimir_Resultados(Fila, ArrayK(), kk, NC)

For kk = 1 To NC

Worksheets(“Hoja1”).Cells(Fila, 1) = “K(” & kk & “)=” & ArrayK(kk)

Fila = Fila + 1

Next kk

End Sub

Sub Calculo_Funcion_Temperatura_Burbuja_Rocio(ArrayK(), Arrayx(), T, Fila, Sum, Sumr)

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

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

‘Worksheets(“Hoja1”).Cells(Fila, 3) = “Fr(T)”

Fila = Fila + 1

For T = 100 To 600 Step 10

For kk = 1 To NC

ArrayK(kk) = T * (ArrayA(1, kk) + ArrayA(2, kk) * T + ArrayA(3, kk) * T ^ 2 + ArrayA(4, kk) * T ^ 3) ^ 3

Sum = Sum + ArrayK(kk) * Arrayx(kk)

‘Sumr = Sumr + (1 / ArrayK(kk)) * Arrayx(kk)

Next kk

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

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

‘Worksheets(“Hoja1”).Cells(Fila, 3) = Sumr – 1

Sum = 0: Sumr = 0: Fila = Fila + 1

Next T

End Sub

Sub Calculo_del_Punto_Burbuja(Tv, Tn, ArrayA(), ArrayK(), ArrayKp(), Fila, Arrayx())

Sum = 0: Sump = 0: L = 0

De_Nuevo:

For kk = 1 To NC

ArrayK(kk) = Tv * (ArrayA(1, kk) + ArrayA(2, kk) * Tv + ArrayA(3, kk) * Tv ^ 2 + ArrayA(4, kk) * Tv ^ 3) ^ 3

Sum = Sum + ArrayK(kk) * Arrayx(kk)

Next kk

For kk = 1 To NC

MM = (ArrayA(1, kk) + ArrayA(2, kk) * Tv + ArrayA(3, kk) * Tv ^ 2 + ArrayA(4, kk) * Tv ^ 3) * 2

MMM = ArrayA(2, kk) + 2 * ArrayA(3, kk) * Tv + 3 * ArrayA(4, kk) * Tv ^ 2

MMMM = (ArrayA(1, kk) + ArrayA(2, kk) * Tv + ArrayA(3, kk) * Tv ^ 2 + ArrayA(4, kk) * Tv ^ 3) ^ 3

ArrayKp(kk) = 3 * Tv * MM * MMM + MMMM

Sump = Sump + ArrayKp(kk) * Arrayx(kk)

Next kk

Tn = Tv – (((Sum – 1) / Sump))

If Abs(Tn – Tv) < 0.01 Then

Worksheets(“Hoja1”).Cells(Fila, 1) = “Pto.Burbuja” & Tn & ” oR” & “L=” & L

Exit Sub

Else

Tv = Tn: Sum = 0: Sump = 0

L = L + 1

End If

If L > 100 Then

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

Exit Sub

Else

GoTo De_Nuevo

End If

End Sub

Trabajos Citados

Buford D. Smith et al. (1963). Design of Equilibrium Stage Processes. New York, New York, U.S.A.: McGraw-Hill Book Company.

Charles Holland (1963). Multicomponent Distillation. Englewood Cliffs, New Jersey: Prentice Hall, Inc.

George Granger Brown et al. (1960, 7th printing). Unit Operations. (J. W. Sons, Ed.) New York, New York, U.S.A.: John Wiley & Sons.