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.

 

 

 

Acerca de

Professor of modeling and simulation, and process design at Escuela Politécnica Nacional, in Quito, Ecuador. . In the past I was a P4, P5, and D1 at the Organization for the Prohibition of Chemical Weapons, located in the Kingdom pf the Netherlands

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

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

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

Únete a otros 1.771 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
A %d blogueros les gusta esto: