10: Estrategia de Convergencia

Ciencia de la destilación (una mezcla de Química e Ingeniería Química)

Se trata de la Parte X, Estrategia de Convergencia de una serie de diez partes de artículos técnicos sobre la Ciencia de la Destilación, como se practica actualmente a nivel industrial. Véase también la Parte I, Panorama general para los comentarios introductorios, los alcances de la serie de artículos y la nomenclatura.

Este último artículo discute cómo abordar la resolución de ecuaciones no intrínsecas y sistemas de ecuaciones, ya que estos forman parte de la práctica de la ciencia de la destilación. Si bien se puede argumentar que este tema está técnicamente más relacionado con las matemáticas o la informática, las técnicas son herramientas indispensables para obtener soluciones operablemente precisas.

La mayoría de las relaciones en la Ciencia de la Destilación presentan ecuaciones que se escriben con la temperatura como variable inicial, y luego una resuelve para los otros parámetros. Por ejemplo, casi todas las relaciones de presión de vapor se escriben intrínsecamente con VP en función de la temperatura. Pero frecuentemente la declaración del problema comienza con una presión fija, por lo que necesita haber un procedimiento iterativo para resolver la temperatura.

Las relaciones más complejas que involucran dos variables iniciales se escriben en temperatura y presión, con la variable dependiente de la ecuación escrita en función de temperatura (absoluta) y presión (absoluta). Un ejemplo es una Ecuación de Estado, donde Z se escribe como f (T, P), siendo el volumen molar (V) retroresuelto a partir de la relación definitoria PV=ZRT. Algunas relaciones científicas se escriben en variables adimensionales “reducidas”, donde T, P y V se dividen por su valor de punto crítico (p. ej., T _r =T/T _c, P _r =P/P _c), lo que permite formas de ecuación más globales. La Ecuación de Estado de Peng-Robinson es un ejemplo de uso de propiedades reducidas, aunque la Ecuación 5-4 enmascara parcialmente eso al romper la ecuación en pedazos.

La técnica más básica para la convergencia de soluciones en una relación no intrínseca se basa en la “estrategia de juego”: (1) elegir el valor de una posible variable independiente y trabajarlo a través de una solución, comparándolo con la variable dependiente final deseada; luego (2) mover un poco la elección del valor y ver cómo eso ayuda. Esto es solo una forma metódica de “prueba y error”, y eventualmente conducirá a una solución. Sin embargo, tal estrategia es casi imposible de automatizar para facilitar la programación, como en una macro de hoja de cálculo o en un programa BASIC independiente. Además, puede ser difícil comenzar, puede conducir a una situación inestable que gira alrededor de una solución en lugar de converger de manera estable con la precisión requerida.

Método Newton-Raphson

Una mejor opción de convergencia es una modificación de Newton-Raphson, también conocido como Método de Newton. En las instrucciones siguientes, “X” es la variable independiente de la relación y “Y” es la variable dependiente. El problema es que tienes un valor objetivo deseado para “Y” y tienes que buscar iterativamente el valor correcto para “X”, y no hay forma de reescribir la relación para que “Y” sea la variable independiente. Los subíndices crecientes denotan la progresión en la mejora de la solución, de tal manera que “Y _n" es idéntico al valor objetivo deseado de “Y” = “Y _T “, y “X _n" es el valor para “X” que produce esa solución.

1. Si hay una manera de obtener una estimación muy aproximada de la solución, comience con eso como el primer ensayo = X ₀, y trabaje la relación científica a través de una primera solución posible, Y ₀. (Si no hay una estimación aproximada posible, comience con una inicial conocida “segura” X ₀ y calcule Y ₀.) Luego haga un paso ΔX extremadamente pequeño (digamos 0.1% de cambio en X ₀). En la siguiente iteración, X ₁ =X ₀ +ΔX sería evaluado para obtener Y ₁. Si el cambio en ΔY (i.e., ΔY= Y ₁ -Y ₀) es demasiado intrascendente para ser visto, aumenta el valor del paso ΔX, hasta que haya algún cambio notable en ΔY. (Esta es una consideración importante ya que algunas relaciones son complejas, por lo que es difícil elegir ese primer paso ΔX para obtener un cambio razonable en Y). No importa si el cambio a X (i.e., ΔX) es positivo o negativo: eso se resolverá más adelante.
2. Calcular la cantidad\(m=(\Delta Y/ \Delta X) = (Y_{1} -Y_{0})/(X_{1} -X_{0}) \nonumber\), que se aproxima a la primera derivada de la relación científica. “m” puede ser positiva o negativa (razón por la cual el signo en ΔX no importó). Un valor positivo para “m” indica que Y aumenta con el aumento de X; o un valor negativo para “m” indica el efecto contrario. Tomar nota del signo de “m”, para ver si su signo tiene sentido común.
3. Con Y _T como valor de solución objetivo, haga la siguiente conjetura para X como\[X_{2} =X_{1} + K \times (Y_{T} -Y_{1})/m \label{10-1} \]
4. La variable “K” es el factor de relajación, con valores de “K” más cercanos a la unidad dando una convergencia más rápida, pero a riesgo de sobreimpulso o inestabilidad de convergencia; y valores de “K” más cercanos a cero dando convergencia lenta pero estable. Hasta que se sepa algo sobre la relación científica, se recomienda utilizar un valor de “K” de aproximadamente 0.2-0.3 al principio, especialmente cuando se trata de relaciones logarítmicas. Siempre se puede cambiar el valor de la “K” una vez que vea lo sensible que es la relación científica.
5. Para la iteración “i th”, la ecuación de iteración se convierte en:\[ X_{i+1}=X_{i} +K \times (Y_{T} -Y_{i}) \times (X_{i} -X_{i-1}) / (Y_{i} -Y_{i-1}) \label{10-2} \]
6. Continuar iterando “n” veces hasta que la solución, Y _n, esté adecuadamente cerca de la solución objetivo, Y _T. Se sugiere que la precisión de iteración de bucle debe ser un decimal mayor que la precisión deseada por la solución final, para permitir el error de redondeo y la precisión del cálculo por computadora.

Al hacer cálculos automatizados por computadora, ya sea a través de una macro de hoja de cálculo o un programa BASIC independiente, es una buena idea verificar errores en cada iteración para ver si hay división por cero, y limitar las iteraciones de bucle numérico (digamos a un número alto como 100). De lo contrario, el programa de iteración puede “colgarse” en un bucle infinito. Si está obteniendo esencialmente valor cero para la cantidad\((Y_{n} -Y_{n-1}) \nonumber\) o requiriendo un número excesivo de bucles para converger, eso requiere ajustar el factor de relajación, “K”, arriba o abajo para ajustar la sensibilidad de convergencia.

Algunas ideas para hacer esa primera estimación aproximada, para comenzar la iteración, son:

• Para encontrar la temperatura de saturación para una presión conocida (es decir, invertir una relación de presión de vapor), comience con el uso de una relación VP básica (inexacta, pero simple), como\(Ln(P) \alpha - \Delta H_{vap}/RT \nonumber\) usar la temperatura de ebullición ambiente conocida, T _b @ P = 1 atmósfera y una estimación razonable para latente de vaporización,\(\Delta H_{vap} \nonumber\). Entonces la estimación aproximada de inicio es\(T_{0}=T_{b}- \Delta H/[R \times Ln(P_{T})] \nonumber\), donde P _T es el objetivo de presión conocido, en atmósferas. Esto dará un T ₀ estimado alto para P _T > atmosférico, y un T ₀ estimado bajo para P _T < atmosférico. Es más importante comenzar con una estimación “segura” que cerrar la primera iteración, para evitar las trampas de las soluciones automatizadas mencionadas anteriormente.
• Los gráficos y tablas de presión de vapor son buenas formas de comenzar, aunque a veces difíciles de programar en una secuencia automatizada.
• Al iterar para la temperatura de saturación de una mezcla, comience con la suposición de una relación lineal entre las temperaturas de saturación de los dos fluidos a esa presión; luego interpolar entre esas dos temperaturas de componentes puros por fracción molar de la mezcla. Siempre hay curvatura entre las temperaturas de saturación de dos fluidos, por lo que asumir la linealidad no es correcta; pero inicia el bucle de iteración.
• El inicio de un bucle de iteración para los coeficientes de fugacidad solo debe hacerse después de una mirada superficial a la presión reducida, P _r, ya que estas relaciones implican calcular la compresibilidad, Z. Para los coeficientes de fugacidad líquida, evite comenzar o acercarse demasiado a P _r =0.9, debido a que el solución automatizada a ecuaciones cúbicas se vuelve demasiado sensible, y una de las raíces puede convertirse en un “número imaginario” (es decir, un número veces la raíz cuadrada de uno). Para los coeficientes de fugacidad de vapor, haga una comprobación inicial para ver qué ecuación Z = f (T _r, P _r) es la mejor para usar (es decir, Ecuación 6-6 o Ecuación 6-7): una es mejor para valores más bajos de P _r y otra mejor para valores más altos. Si intenta cruzar entre ecuaciones en medio de una solución, esa discontinuidad podría hacer que una iteración automatizada “explote”.
• Resista la tentación de usar Newton-Raphson para obtener las raíces de la relación cúbica Z: hay tres raíces, y solo dos son científicamente significativas. Además, es necesario conocer las tres raíces para poder juzgar cuál es Z _v (la más grande de las tres) y cuál es Z _L (la más pequeña de las tres). En su lugar, pasar por el procedimiento en la Parte V y resolver la ecuación cúbica para Z, utilizando el método trigonométrico exacto de la Ecuación 5-4.
• El inicio del bucle para los coeficientes de actividad nunca debe hacerse en ninguno de los extremos de las fracciones molares (cero o unidad). Las relaciones del Coeficiente de Actividad Líquida son de forma compleja y pueden tener\(\gamma/ \Delta X \nonumber\) valores que son casi cero (cerca de X =1) o muy altos (cerca de X = 0). Para separaciones masivas, comience en algún lugar entre 0.3< X <0.7, y use un valor moderado para el factor de relajación, “K” (digamos 1/3) en E quation (\ ref {10-2}). Para los sistemas de destilación de impurezas traza, ya sea comience más cerca de X =0.1 o en X =0.9, además de usar un parámetro de relajación conservador, digamos K=0.1.

Estrategia de Convergencia para Ecuaciones VP de Thek-Stiel modificadas

El “aguijón” de elaborar iteraciones de bucle anidado se ha eliminado dando las soluciones para la ecuación modificada de Thek-Stiel VP (Ecuación 4-3) en forma tabular (Tabla 4-1). La obtención de soluciones de presión de vapor para fluidos distintos de los clorosilanos y sus impurezas (ver Parte IV) requiere estrategias de convergencia más complejas que Newton-Raphson.

Estas soluciones VP involucran un sistema de bucles anidados: el valor del coeficiente de Watson “q” depende de evaluar la ecuación de presión de vapor candidata a T _r =0.7 para obtener la estimación del factor acéntrico del bucle; y obtener la estimación del bucle de “c” y “n” depende de encontrar un valor de α c que da la temperatura de saturación correcta a presión ambiente (también conocido como el punto de ebullición normal). A continuación se aconseja la siguiente estrategia de solución, para fluidos distintos a los del Cuadro 4-1 (o en caso de que se quiera verificar los valores de la tabla).

Primero, asegúrese de que todos los valores para T _b, T _c, P _c estén bien validados y que se conozca el grado de enlace de hidrógeno, según las definiciones de la Parte IV. Si el nuevo fluido a investigar es de un grupo de tabla periódica que no sea III, IV o V, la solución implicará una solución más intrincada que la de abajo, ya que “k” (es decir, el grado de enlace de hidrógeno) no está bien determinado para otros que no sean los Grupos III-V, y se necesitará una doble iteración. También asegúrese de que los datos de presión de vapor también estén bien validados. Una buena comprobación en eso es utilizar una extrapolación simplista de los datos VP conocidos, para ver si el valor aproximado del factor acéntrico coincide con los valores de T _c y P _c. Consulte la Parte III para más detalles sobre el factor acéntrico, y especialmente la Ecuación 3-3 Por ejemplo, si una extrapolación de datos VP usando la relación VP de Antoine (Ecuación 2-3) junto con valores para T _c y P _c, termina dando un valor de factor acéntrico superior a 0.3 para una molécula compacta, o un factor acéntrico por debajo de 0.1 para una molécula compleja, algo está mal con los datos (ya sea los datos VP o los valores de T _c y P _c).

Para la Ecuación 4-3, el valor del parámetro de enlace de hidrógeno, k, se establece mediante datos experimentales para hidruros, cloruros, oxicloruros, clorhidratos intermedios y varios hidrocarburos alifáticos de los Grupos III, IV y V. El valor de k oscila entre 0.12 y cero, y se tabula en la Parte IV para muchos fluidos. Las reglas para determinar “k” son un poco complejas y están disponibles a pedido del autor. En resumen, para los hidruros del Grupo III, IV y V, k=0.1197; 0.0914; y 0.1009 respectivamente. Para todos los compuestos completamente clorados o metilados de los Grupos III, IV y V, k = 0.0291. Para hidrocarburos alifáticos, k = 0. Para otros intermedios de hidrocloruro, se ha determinado una función cuadrática para interpolar entre todos los hidruros y todos los cloruros. Estas reglas para los valores de “k” se utilizan en las entradas del Cuadro 4-3.

Tener una estimación aproximada del calor latente de vaporización a presión ambiente es útil (pero no necesario) para la primera conjetura en la variable de calor T-S. Si no hay estimación disponible, consulte uno de los muchos manuales de compuestos químicos, como el Manual de Propiedad Química de Yaw o el Manual de Química de Lange y busque un valor de calor latente para un compuesto de complejidad similar. Siempre es mejor comenzar la solución de un nuevo fluido a la ecuación Thek-Stiel VP con una estimación baja para la variable térmica T-S, y luego aumentar.

La estrategia global de convergencia de la solución a la ecuación modificada de Thek-Stiel VP implica un bucle de iteración triple anidado, con el diagrama de flujo lógico de la estrategia que se muestra en la Figura 10-1. La solución se obtiene mejorando las soluciones candidatas para la variable térmica T-S, el factor acéntrico (ω) y el parámetro α _c, hasta que todos los bucles sean convergentes.

El bucle de iteración más interno es el\(\alpha_{c} \nonumber\) parámetro (un parámetro de consistencia termodinámica de la Parte IV), que se ajusta hasta que la ecuación VP da la temperatura de saturación correcta a presión atmosférica, también conocido como el punto de ebullición normal T _b. Ese bucle se inicia usando el método de estimación de Reidel: calculando el valor aproximado de\(\psi_{b}= (-35 +36/T_{br} + 42LnT_{br} -T_{br}^6 \nonumber\), y luego enchufando ese valor aproximado de\(\psi_{b} \nonumber\) en la siguiente ecuación:

\[ \alpha_{c} \approx (0.315 \psi_{b} + LnP_{c})/(0.838 \psi_{b} + LnT_{rb}) \label{10-3} \]

La mitad de los bucles anidados es el factor acéntrico, que se puede iniciar convenientemente usando la aproximación del factor acéntrico de Lee-Kister (ω) (que no es muy precisa, pero una buena estimación de inicio de bucle):

La aproximación de Lee-Kister para el factor acéntrico (ω) utiliza convenientemente solo P _c y T _br =T _b /T _c:

\[\omega \cong \frac{-lnP_{c}+A+B/T_{br} +C*LnT_{br} +D*T_{br}^6}{E+F/T_{br} +GLnT_{br}+H*T_{br}^6} \label{10-4} \]

A = -5.92714 E = 15.2518

B = 6.09648 F = -15.6875

C = 1.28862 G = -13.4721

D = -0.169347 H=0.43577

Con base en este valor aproximado estimado de ω, se puede determinar un valor del coeficiente de Watson “q” (de la Ecuación 4-3) por la relación

\[ q = 0.37028 + 0.065404*ω \label{10-5} \]

El bucle más externo es la variable térmica T-S, “T-S ΔH _v”, que en realidad no tiene mucho significado científico, sino que forma parte de la relación Thek-Stiel. Este parámetro se itera en función de la función de error comparando la salida de la ecuación VP con los diversos puntos de datos VP conocidos (o estimados). Para obtener el mejor ajuste general, sugiero incluir datos cercanos a T _r = 0.7 (donde se evalúa el factor acéntrico). Al obtener las diversas soluciones T-S VP Equation mostradas en la Tabla 4-1, elegí no automatizar este bucle externo, sino confiar en la “estrategia de juego” más básica ⇒ siempre comenzando bajo, para monitorear mejor la convergencia. Para la automatización informática, utilicé un programa BASIC (el software real utilizado fue Liberty Basic versión 4.03, seleccionado por familiaridad previa). La suposición inicial para la iteración del bucle externo se basó en el _vap ΔH atmosférico dado en los manuales de química y de fuentes de internet, o en base a compuestos similares donde no se disponía de una estimación publicada de ΔH _vap atmosférico. Esta elección de conjetura inicial se basó en el criterio de “conjetura segura”, y se opuso a una estimación científica.

La Figura 10-1 a continuación muestra esquemáticamente la estrategia de convergencia previamente utilizada con éxito. Después de asegurar que todos los parámetros de entrada y datos son razonables, comience con la suposición inicial para la variable de calor T-S, “T-S ΔH _v”. Después estime una ω inicial por E quation (\ ref {10-4}); a partir de esa estimación un inicio\(\alpha_{c} \nonumber\) por Ecuación 10-3. Ahora calcule los valores para las constantes “A”, “B ₀ -B ₃”, “c” y “n” en la Ecuación 4-3, y vea si se satisface la convergencia del bucle más interno. Si no, iterar el valor de\(\alpha_{c} \nonumber\) hasta que el VP correcto se muestre en T _b de 1.0 atmósfera (es decir, el punto de ebullición normal). Este bucle más interno es mejor converger acercándose al VP @ T _b desde la atmósfera P<1, observándose la convergencia a ± 1.0E-4 atmósferas dentro de aproximadamente diez iteraciones vía Newton-Raphson, usando un factor de relajación de 0.4 en E quation (\ ref {10-2} ). Tenga en cuenta que si la primera conjetura de\(\alpha_{c} \nonumber\) muestra un VP @ T _b de >1.0 atmósferas, simplemente siga reduciendo el valor de\(\alpha_{c} \nonumber\) en un 10% hasta que el enfoque de convergencia sea a partir de valores más bajos de T _b. Tratar de converger desde el “lado alto” puede ser inestable.

Una vez que converja el bucle más interno, comience a cambiar el valor del factor acéntrico (ω) para converger el bucle medio. Este bucle converge con bastante facilidad con una simple sustitución continua. No se requiere Newton-Raphson. Simplemente ponga el último valor de ω, calcule las constantes VP candidatas (“A”, “B ₀ ‑B ₃”, “c” y “n” en la Ecuación 4-3) y determine un nuevo valor para ω basado en el VP a T _r =0.7 (ver Ecuación 3-3 para la definición de ω).

El bucle externo se converge escogiendo ese valor de la variable térmica T-S, “T-S ΔH _v”, que mejor satisface la función de error de los datos VP (es decir, VP predichos en comparación con los VP de punto de datos conocidos). Encontré el mejor ajuste general usando un enfoque de mínimos cuadrados que enfatizó los datos cerca de T _r =0.7, y que comparó el error de la ecuación VP en un nivel relativo (a diferencia de un nivel absoluto). Si se utilizan errores absolutos de la ecuación VP (a diferencia de los relativos), eso sobreenfatiza los datos a temperaturas más altas.