UNIVERSIDAD AUTONOMA METROPOLITANA ...148.206.53.84/tesiuami/UAM1134.pdfbrio de un sistema es un...

UNIVERSIDAD AUTONOMA METROPOLITANA

IZTAPALAPA

LICENCIATURA EN QUIMICA

QUIMICA ANALITICA

'gDeterminación de la composición de equilibrio: Algoritmo VCS"

21 J U N I O 1994

ALUMNA: María del Rocío Trejo Correa

MATRICULA: 84222264

ASESOR: Dr. Alberto Rojas H.

INDICE

INTRODUCCION

I. FUNDAMENTOS TEORICOS 1.1 Segunda Ley de la Termodinámica 7

1.. 2 Terminología 10

1.3 Formulación estequiométrica 11

1.3.1 Restricción de sistema cerrado 11

l . 4 Ilustración del método 17

1.5 Singularidad de las soluciones 20

11. PROBLEMAS NUMERICOS

2.1 Problema de minimizacibn 23

2.2 Método de minimización 23

111. ALGORITMOS

3.1 Algoritmo estequiométrico 24

3 .2 El algoritmo VCS 28

3.3 Formas de representar Is energía libre estándar 33

3.3.1 Datos de energía de formación a partir de constantes de

equilibrio. 33

3.3.2 Programa de computadora basado en el algoritmo VCS. 36

3.4 Variables 61

3.5 Aplicación del algoritmo 62

3.5.1 Salida de datos 65

IV. APLICACION A SISTEMAS DE DISOLUCIONES

4.1 Cálculo de concentraciones en equilibrio 67

4.2 Equilibrios sucesivos 68

4.3 Comparación de resultados 70

V. CONCLUSIONES

VI. BIBLIOGRAFIA

APENDICE A: Errores causados por una mala entrada de datos

1

7

21

24

67

72

72

73

EQUILIBRIO.

Se define un cuerpo en equilibrio como aquél en que todas las fuerzas en oposi-

ción se contrarrestan exactamente de modo que las propiedades macroscópicas del cue2

PO no cambian con el tiempo.

Todos los cuerpos aislados tienden a acercarse a una situación de equilibrio.

Una condición de equilibrio se llama estable porque una vez que un cuerpo se ha

ya alejado de su condición de equilibrio, volverá siempre a ésta. Además de la condi -

ción anterior existen las de equilibrio metaestable e inestable.

Un equilibrio metaestable es aquél en el que el sistema vuelve a su estado ini-

cial si se produce una perturbación, pero se sitúa en una condición de equilibrio di -

ferente, .si se produce una perturbacibn de magnitud suficiente. La condición de equi -

librio inestable es aquella en la que el sistema no vuelve a su condición original ,

si se le somete a una perturbación finita.

Visto en un contexto más amplio el equilibrio es un proceso dinámico a nivel mi -

croscópico, aunque se le trate como Condición estática a nivel macroscópico.

Hay varias clases de equilibrio, equilibrio térmico, mecánico, slectrostático,-

equilibrio de fase, equilibrio quimico, etc. Para cada clase de equilibrio hay una

propiedad termodinámica que dos sistemas deben tener en común.

El término equilibrio termodinámico es usado con respecto a todos l o s posibles

cambios macroscópicos en un sitema donde las moléculas son libres de interactuar en-

tre si.

Conjuntamente con la interacción quimica entre las sustancias iniciles trans

1

curre otra entre los productos de reacción, como resultado de lo cual vuelven a

formarse las sustancias iniciales que interactúan, forman productos, y vuelve a

repetirse el proceso de interacción. A medida que se desarrolla el proceso, la

velocidad de la reacción directa, es decir, la velocidad de formación de produc -

tos, disminuye y la velocidad de reacción inversa aumenta. Cuando ambas veloci-

dades se igualan, sobreviene un estada de equilibrio químico. De esta manera e-

quilibrio es dinámico y desplazable. Con una variación en las condiciones exte-

riores el equilibrio se desplaza en uno u otro sentido de la reacción. Una va--

riación infinitesimal en las condiciones exteriores origina también una va--

riación infinitesimal en el estado de equilibrio. Por consiguiente, las reac--

ciones químicas pueden ocurrir como procesos termodinámicoas de equilibrio, lo

que quiere decir que se les pueden aplicar las condiciones generales de equili-

brio termodinámico.

Como se ha visto, el concepto de equilibrio químico se basa en nociones,

tanto físicas como quimicas. En el sentido químico se ha desarrollado a partir

de ideas sobre e1 cambio químico y de afinidad, y en la época contemporánea, en

ideas sobre velocidad y mecanismos de reacción. También se desarrolló a partir

del concepto de balance de fuerzas en equilibrio estático en los sistemas mecá-

nicos.

Quizá la primera contribución importante en el desarrollo cuantitativo del

equilibrio, aparte de la afinidad, la constituyó el trabajo de Berthelot y St.-

Giles ( 1802) , quines demostraron que se obtenían las mismas composiciones de e

quilibrio en el sistema ácido acético-etanol-agua , a partir de cualquiera de

las dos direcciones en que se procediera; esterificación o hidrólisis. Esto con

dujo a Guldberg y Waage (1864) a una ecuación de equilibrio cuantitativa a tra-

vés de lo que se llamó "Ley de acción de masas", balanceando las "fuerzas qui-

micas" en el equilibrio.

-

2

La derivación de la expresión de la constante de equilibrio junto con l as ex-

presiones cinéticas que se usaban, supuestamente condujeron a la derivación ciné-ti -

ca de la constante de equilibrio tal como la considera, por ejemplo, Van't Hoff en

1898.

Los parámetros que modifican el equilibrio llamaron la atención de Van't Hoff

quien propuso el concepto de equilibrio móvil en base a la evidencia en Cuanto a -

que la elevación de temperatura favorece las reacciones endotérmicas. Le Chatelier

que consideraba las fuerzas que afectan el equilibrio postuló, sin pruebaSrel prin -

cipio de Le Chatelier, que consideraba la forma en que un sistema químico en equi-

librio estable se ajusta a factores externos que tienden a modificar su temperatu-

ra, presión o densidad molar o a un cambio de concentración de las especies dadas.

Este principio establece, que si el sistema sufre cambios que ocurrieran por sí mis

nos, provocarían una modificación en temperatura, presión y otros parámetros opues

tos causados por el efecto externo. En este principio Le Chatelier seguía una afir

nación de Van't Hoff, que consideraba solo efectos de temperatura.

"

-

-

Se han hecho formulaciones más precisas y en ocasiones se hace referencia al

teorema como "Principio de Le Chatelier- Braun" en reconocimiento a la última con

tribución de F. Braun (1887-1888).

-

La formulación termodinámica de 'la constante de equilibrio también la obtuvo

Van't Hoff usando ciclos reversibles. Gibbs (1876) proporcionó otro enfoque al ir

troducir el concepto de potencial químico.

Aunque l o s principios teóricos se establecieron a principios de siglo, el cál

Culo del equilibrio en forma práctica se dificultó a medida que se consideraban O

se estudiaban sistemas más complejos conectados con el desarrollo del procesamien-

to químico y con los métodos de análisis químico.

-

El estudio del equilibrio de una reacción química tiene dos intereses; Deter-

minar la composición, en equilibrio, de us sistema bajo ciertas condiciones de pre

SiÓn Y temperatura. El Otro propósito es el de establecer si el estado de equili" -

3

brio de un sistema es un modelo útil para un propósito en particular.

Cuando se habla de "equilibrio de un sistema cerrado" se hace referencia a un

estado con las características siguientes:

1) Independiente del tiempo

2) Independiente de la historia previa del sistema

3) Resistente a fluctuaciones en composición

4) Independiente de la posición dentro del sistema, excluyendo

a sistemas sujetos a eradientes de potencial que impliquen

flujo de masa o energía dentro de él.

El sistema químico de interés pasa por un cambio en su composición química al

aproximarse a un estado de equilibrio desde algún otros estado inicial, fuera de - e

quilibrio. Lo anterior se lleva a cabo por reacción o transferencia de masa o una

combinación de ambas. La posibilidad de la transferencia de masa implica que exis-

te más de una fase en el sistema. De esta forma el término "equilibrio químico" -

incluye tanto e1,de reacción, como el de fases. Sin embargo a nivel microscópico ,

se considera que las moleculas que se encuentran en una fase en un momento dado'no

son las mismas instantes después. Moléculas cercanas a la interfase, con velocida-

des suficientemente altas, vencen las fuerzas superficiales y pasan a la otra fase;

no obstante, el flujo promedio de las moléculas es le mismo en ambas direcciones y

no hay transferencia de materia neta entre las fases.

Una vez que se decide que vale la pena investigar el modelo de equilibrio, se

debe determinar el equilibrio enciertas condiciones y quizá la dependencia de éS--

tas. Esto requiere de la aplicación de la segunda ley de la termodinámica, ya que,

proporciona los criterios de requisito. Para la explotación de dichos criterios es

conveniente construir funciones de potencial que intervienen y dan información a-

propida i ecuaciones de estado,datos de energía libre, etcétera), para especificar los parámetros del problema. Lo anterior hace al problema escencialmente matemáti-

A

co, y como tal, puede considerarse como una optimización o como un problema que im -

plique resolución de un sistema de ecuaciones algebraicas simultáneas no lineales.

Los algoritmos para el cálculo de equilibrio químico anteriores a l o s prime--

ros anos de la década de l o s ! j o g s , estaban orientados principalmente a métodos de

cálculo manual, y giraban alrededor de expresiones de constantes de equilibrio es-

critas para un conjunto de ecuaciones estequiométricas en términos de concentracio -

nes, lo cual implica la solución de un conjunto de ecuaciones algebraicas no linea

les, casi siempre en un número igual a la diferencia entre el número de especies y

los elementos invclucrados.

Los cimientos para l o s algoritmos de propósitos generales para el cálculo de

equilibrio químico fueron colocados por Brinkley (1947) quien esbozó un procedi---

miento sistemático enfocado a la atención en un conjunto de especies "componentes"

cuyo número generalmente es igual al número de elementos M que intervienen. El al-

goritmo de Brinkley se derivó primero considerando las expresiones de las constan-

tes de equilibrio.

Considerando el problema del cálculo de equilibrio como un problema de optimi

zación no lineal, White ( 1 9 5 3 ) , desarrolló un algoritmo que resolvía el problema -

"minimizando" directamente la energía libre. Esto significaba que no se usaban e-

cuaciones o reacciones estequiométricas.

-

Desde los últimos anos de la década de los 50's existe un nuevo interés por

implementar con las computadoras algunas de las técnicas que utilizaban ecuaciones

estequiométricas y variables de grado de extensión de la reacción ( Villar, 1 9 5 9 ) .

Estso algoritmos estequiométricos también pueden ser motivados por considera-

ciones de energía libre. La única diferencia es que las variables independientes u

sadas Para efectuar la minimización de energía libre son el conjunto de variables

de grado de extensión de la reacción, mientras que en los métodos directos , las

restricciones sobre el balance de masa se manejan por cálculo utilizando los multi

plicadores de Lagrange.

-

-

5

Los métodos de cálculo a los que se ha hecho referencia se desarrollaron ori -

ginalmente para ser utilizados en sistemas ideales. La rnayoria de los métodos con -

sidera una sola fase gaseosa y cualquier número de fases condensadas Furas.

En el cálculo del equilibrio la descripción del sisterria puede hacerse de for -

ma estequiométrica, que es válida aun sin importar si el sistema está o no en e-

quilibrio y de forma termodinámica.

En este trabajo se presenta algoritmo basado en una descripción termodi.nij

mica que hace el cálculo del equilibrio minimizando la energia libre de Gibbs, cg

yas variables son las coordenadas de reacción de las posibles reacciones que se

lleven a cabo en el sistema. Este algoritmo es el que se debe a Villars-Cruise---

Smith (VCS), considerado como un algoritmo de minimización descendente. El propó-

sito principal de este trabajo es describir los fundamentos teóricos de una forma

general tanto en la parte teórica,como en la matemática, así como también es ha-

cer una comparación de los resultados obtenidos con el algoritmo y los obtenidos

usando un modelo de dos equilibrios para el cálculo de concentraciones en equili

brio de un sistema de mezclas, y principalmente para el cálculo del pH que se uti

liza en los cursos de ensenanza de quimica analítica de la UAM-I.

-

-

6

I. FUNDAMENTOS TEORICOS

1 . 1 Segunda Ley de la T e r m o d i n h i c a .

La segunda ley de la termodinknica proporciona variss funciones potenciales

que regulan la direccibn de los procesos naturales o espontáneos. La función de

potencial adecuada para cada situacih est$ definida por la selección de varia--

bles termodinámicas que se consideren independientes, éstas definen el estado de

un sistema y PGI^ tanto, las funciones se conocen como funciones de estado. Cual-

quier cambio en dichas funciones es independiente de la trayectoria del cambio.

Entre las funciones más importantes están la de Gibbs, la de Helmholtz y la

de entropía. Para cada una existe un enunciado de la segunda ley que incluye tan

to el criterio para que se produzca un proceso natural, como para su estado fi--

nal de equilibrio, debiendo tomarse en cuenta las restricciones.

Así pues para cada una de éstas se tiene:

1 ) dSad > o adiabático

2) dA 2 0 Func. Helmholtz

3) dG 2 0 Func. de Gibbs t,v

t 9.P

en donde d se refiere a un cambio infinitesimal, la desigualdad a un proceso es--

pontáneo y la igualdad al equilibrio. Dependiendo de las restricciones S se en---

cuentra en un máximo local, A y G en un mínimo. Debido a estas restricciones de

temperatura y presión, la más importante es la función de Gibbs.

Un sistema químico de un solo componente y además homogéneo, se describe ter

modinámicamente por alguno de los siguientes conjuntos naturales de funciones de

-

estado y variables independientes.

4 ) U = U(S,V,n)

5) H = H(S,P,n)

6) A = A(T,V,n)

7) G = G(T,P,n)

energía interna

entalpía

energía de Helmholtz

energía libre de Gibbs

Cada una de éstas también es homogénea de grado uno en cada número de moles,

7

Y da lugar a su correspondiente para la diferencial completa;

8 ) dU = TdS - PdV +qpi 7 dn i

9 ) dH = TdS + VdP + 1 pidni 10) dA = -SdT - PdV + xpidn

dl

L:i 00

G 1 i

11) dG = -SdT + VdP +>p.dn. i,zL 1 1

Debido a la propiedad de homogeneidad de éstas funciones, pi depende sólo del

estado intensivo del sistema, definido por T y P y la composición. Dado que la más

importante es G la descripción termodinámica se hará en términos de ésta función . De la ecuación 11 y de la definición de G , tomando derivadas de G y pi respec

to a T y P se tiene:

11 tenemos, con T y P constantes;

8

17) G(T,P,n) = "ipi i=l

diferenciando 17 y comparando con la ecuación 11

18) SdT -VdP ++ n. dp 1 i

que es la ecuación de Gibbs-Duhem. Todas estas acuaciones se aplican a un sistema

homogéneo, o a cada fase en el sistema heterogéneo. Dado que puede distinguirse u -

na especie de otra con la misma fórmula química, pero dependiendo de la fase don-

de se encuentre, el potencial químico y las cantidades molares parciales de una

fase determinada serán dadas por las variables que definen el estado de esa fase

y las implicaciones para G y p son las siguientes: i

1) Las ecuaciones 7,11,13,14 y 17 se aplican al sistema completo a cada fase, siempre que las cantidades extensivas G,H,S y V se relacionen al sistema completo o a la fase en consideración.

2) Las ecuaciones 12, 15,16 ,Gi. y se aplican a cada especie.

3) La ecuación 18 se aplica a cada fase.

i

Para que un sistema de fase o de fases múltiples esté en equilibrio, G debe

encontrarse en un mínimo global, el cual esta sujeto a las restricciones de sis-

tema cerrado y de no negatividad de las condiciones termodinámicas establecidas a

T y P fijos. Suponiendo n O se ignora la posibilidad en la restricción de no ne-

gatividad, de que n f O. En el equilibrio se trata con i

dGt,p = O (19)

condición necesaria pero no suficiente.

El problema fundamental es expresar G en función de n y buscar valores

ésta que hagan a G minima sujeta a las restricciones. En principio se supone i para

que

se conocen los valores los valores del vector de abundancia de elementos b,la tern

peratura y la presión, as1 como l o s datos de energía libre.

t -

9

Lo anterior lleva a cómo realizar la formulación del problema de minimiza-

ción. Existen varios tipos de formulaciones, en particular se trata con dos ti-

pos a los que se hace referencia como:

1. FORMULACION ESTEQUIOMETRICA. La restricción del sistema cerrado se trata p o r medio de ecuaciones estequiométricas, dando por resultado esencialmente un problema de minimización sin restricciones.

2. FORMULACION NO ESTEQUIOMETRICA. La restricción del sistema cerrado se trúta por medio de multiplicadores de Lagrange.

Antes de comenzar la explicación de las formulaciones, es necesario mencio-

nar algunas definiciones que se consideran pertinentes para el mejcr entendimien -

to de dichas formulaciones.

1.2 Terminología.

Definición de términos, en su mayoría rel.acionados con un sistema cerrado.

Especie química: Entidad química diferenciable de otras por; su fórmula molecu--

lar 3 por su estructura (isómeros), si esto falta por la fase en que se encuentra.

Sustancia química: Entidad química diferenciable por su fórmula molecular o es-- tructura molecular.

Componente : Uno de un conjunto de C especies del sistema químico, cuyo conjun- to de vectores fórmula (a,, a2, a3,..., a ) satisface C = rango(A)

Sistema químico: Conjunto de especies químicas y elementos denotados por un conjunto ordenado de elemenetos contenidos en ‘el de la siguiente - forma

C

A es la fórmula molecular y designación estructural y de fase y E es el elemento . La lista de elementos incluye isótopos y car ga protónica si se incluyen especies iónicas, y la designación (Xl, X2, ... , X es para cada especie inerte.

-

n

Vector fórmula: Vector de subindices de l o s elementos en la fórmula molecular de una especie, como en

C6H5N02 a = ( 6 , 5 , 1 , 2 )

Matriz de fórmula A: Matriz de MxN, donde la columna i es a ; A es la matriz de coeficientes de las ecuaciones de ahundancia de elemen-

10

tos.

Vector de abundancia de elementos b: Vector de números reales (casi siempre no ne gativos) que representan el número de moles de elementos en una

-

cantidad base del sistema químico.

\lector de abundancia de especies n : Vector de números reales no negativos que re presentan el número de moles de las especies en una cantidad base

-

del sistema químico, y n = (nl,n 2 , . . . ,nN); ni>/O, n también indi ca la composición.

-

Sistema químico cerrado: Aquel para el cual todos los n posibles satisfacen la e-- cuación 22, para un valor dado de b.

Vector cambio de abundancia de especies: dn = n - n , los cambios en el número de moles entre los estados composicionales 1 y 2 del sistema quimi ~

L I

co cerrado. Este vector debe satisfacer la acuación 23. -

l. 3 Formulación estequiométrica.

La restriccioh que se impone sobre el cálculo de equilibrio químico debido al

requerimiento de conservación de los elementos en un sistema cerrado que sufre una

transformación química, está intimamente ligada a lo que se conoce generalmente cc

mo estequiometría química, ya sea expresada en términos de ecuaciones de conserva-

ción o de ecuaciones químicas.

El propósito específico de ‘esta formulación es desarrollar la estequiometria

química para un sistema cerrado en forma adecuada para incorporarla a un algoritmo

de cálculo de equilibrio.

Ya que las ecuaciones de conservación son en sí ecuaciones algebraicas linea-

les, el uso del algebra lineal proporciona las herramientas necesarias para el cál -

culo.

1.3.1 Restricción de sistema cerrado.

Un sistema cerrado tiene una masa fija, es decir, no intercambia materia, pe-

ro si puede intercambiar energía. La importancia en los cálculos de equilibrio de

este tipo de sistemas, estriba en las condiciones de equilibrio de la termodinámi-

ca que se aplican principalmente a este sistema. Cualquier descripción de sistema

cerrado es una Ley de Conservación de la Masa. Un sistema cerrado se define por un

conjunto de ecuaciones de abundancia de elementos que expresan la conservación de

11

los mismos, y que constituyen las especies del sistema. Para cada elemento existe

una ecuación con la siguiente forma:

K" N

)-tini = b k i=l

; k= 1,2,3, ..., M (20)

a es el subindice del k-ésimo elemento en la fórmula molecular de la especie i:

n es el número de moles de i y b es el número de moles fijo del k-ésimo elemen-

to; M es el número de elementos y N es el número de especies. La ecuación ante--

ki

i k

rior también puede escribirse de la forma que pueda expresar cambio de composi---

ción;

N

x a ki dn i - O - ; k= 1,2,3, ..., M (21)

i=l

dn. es el cambio del número de moles de la i-ésima especie, entre dos estados de

composición. En forma matricial las ecuaciones respectivas son:

1

A es la matriz de fórmula, n es el vector de abundancia de la especie (vector co-

lumna) y b es el vector de abundancia de elementos.

Ejemplo 1 .

Suponga un estado inicial con NH3 y 0 , , en relación 4:7. Intervienen las

especies; NH3, O,, NO, NO2 y H 2 0 . De la relación molar 4:7, se establece que la

base del sistema es

b = 4 , b H = 1 2 , b 0 = 1 4 N donde las ecuaciones 20's correspondientes para cada elemento en el orden N,H y 0

serán

1 n + O no, + 1 nNO

3 n + O no, + O n + O n + 2 n

NH 3 + nNOz + O n H = b = 4 2 0

NH 3 NO NO 2 Hz0 = b = 12

12

O n + 2 n + I n NH 3 O 2 NO

la ecuación (2'7) correspondiente es

~ 1 0 1 1 0 I

A

+ 2 n + NO z

4

12

14

b

El número máximo de ecuaciones linealmente independientes, que es el mismo que

el de renglones (o columr,as) linealmente independientes en la matriz A , esta dado

por el rango de A . Precisamente uno de los propósitos de la estequiometría química

es determinar este número máximo de ecuaciones linealmente independientes.

Las ecuaciones de conservación generalmente no proporcionan toda la informa--

ción requerida para la determinación de la composición n . Las relaciones adiciona--

les entre las variables que se requieren para determianr cualquier estado composi--

cional es el número de grados de libertad estequiométricos Fs, que es la diferen--

cia entre el número de variables N que se utilizan para describir la composición y

el número máximo de ecuaciones linealmente independientes. La diferencia entre N y

el númro máximo de ecuaciones de abundancia de elementos linealmente independientes

se representa por R. Fs y R en general no son iguales. La estequiometría química si

permite determinar valores de Fs y R para un sistema dado y establecer un conjunto

permisible de ecuaciones químicas.

La generación de la estequiometría química y de las ecuaciones químicas a par-

tir de las ecuaciones de conservación se da a través de un tratamiento general, co-

menzando con la solución de la ecuación (20) que es

c R

13

O

n es cualquier solución particualr ( una composisción inicial), (vl,v 2,...,v ) es

cualquier conjunto de R soluciones linealmente independientes de la ecuación homo-

génea correspondiente a (22), las § . son un conjunto de parámetros reales, cada

v es denominado vector estequiométrico cuya definición general es:

Vector estequiométrico: Cualquier vector diferente de cero, de N números reales

que satisfaga Av = O por tanto

R

J

.i

A V = O (vj f O ) j = 1,2,3,. . . ,R (25) j

o también M

k = 1 , 2 , 3 , ..., M j = 1,2,3, ..., R (26)

y Vij f O cuando menso en un i en cada j . R es el número máximo de soluciones linealmente independientes de (25) y esta dado por

R = N - C (27)

donde

C = rango ( A ) (28)

Generalmente, aunque no siempre, C = M. Una forma alterna de escribir la e-

cuación (24) es O

ni = ni +. 7 R vij Sj i = 1,2, ..., N (24b)

O

para n fijo

i = 1 , 2 , ..., N j = 1,2, ..., R ( 2 9 )

por tanto -t es la rapidez de cambio de número de moles de especie i , n con reg

pecto a 4 La ecuación (24) se puede considerar como una transformación lineal de

las N variables independientes n a las R variables independientes s.

ij i'

j'

14

El significado químico de la ecuación ( 2 4 ) , es que cualquier estado composi-

cional del sistema n se puede representar en términos de cualquier estado particu-

lar n y de una combinación lineal de un conjunto de R vectores linealmente inde - O

pendientes V que satisfacen a la ecuación (23). La ecuación (25) conduce al con-

cepto de ecuaciones químicas. TJna forma abreviada de escribir esta ecuación es j

= o i=l

j = 1,2,...,R (30)

v tiene doble subíndice para identificar su asociación tanto con los componentes

como con la reacción, así v designa el número estequiométrico del componente i

en la reacción j, y a cada reacción por separado se le asocia una coordenada

de reacción.

ij

ij

'j

Para la utilización de estos conceptos en casos reales es conveniente determi -

nar R y el conjunto de R vectores estequiométricos linealmente independientes Iv.1

a este conjunto se le llama conjunto completo de vectores estequiométricos para un J

sistema con matriz de fórmula A . Una forma concisa de escribir cualquier conjunto

de vectores estequiométricos es definiendo una matriz N cuyas columnas son los v .

vectores. Ij

N = ( VI, V2,...7V 9

(31)

cuando q = R y todos los Y son linealmente independientes, N es "completa" y se

definen los siguientes conceptos: j

N , MATRIZ ESTEQUIOMETRICA COMPLETA. N x R cuyas R columnas son vectores estequio métricos linealmente independientes, con R = N- rango( A), lo cual permite escri bir la ecuación (25) como

-

A N = O (32)

Análogamente, un conjunto completo de ecuaciones químicas, es un conjunto de

ecuaciones (30) donde las yij forman una matriz estequiométrica completa N.

Este conjunto se genera a partir de la lista tle especies que se supone o se de

15

muestra que están presentes sin requerir del conocimiento de las reacciones que se

llevan a cabo.

El significado de C = rango (A) es, que al darse R valores de n es posible

obtener C valores de n resolviendo las ecuaciones (22), siempre que sean lineal.--

mente independientes los vectores fórmula de estos C valores, lo cual es equivalen -

te a dividir las especies en 30s grupos; componentes ( en número de C ) y no compo

nentes ( en númro de R ) , considerándose a l o s componentes como partes de la cons-

trucción de los no componentes y requiriendose una ecuación para cada uno de éstos

últimos.

i

i

Los algoritmos de procedimiento esteqaiométricos determinan simultáneamente ,

el rango (A) y un conjunto completo de ecuaciones químicas. Este procedimiento es

similar al que se utiliza en la solución de ecuaciones algebraicas lineales por la

reducción de Gauss-Jordan. Implica, par tanto, la reducción de la matriz fórmula A

a una unitaria por medio de operaciones elementales de renglones e intercambio de

columnas.

La matriz unitaria esta representada por;

IC es una matriz de identidad CxC y Z una matriz CxR, donde al menos uno de los

elementos es distinto de cero. C es el rango de A ' y a la vez también lo es de

A . La matriz estequiométrica completa se forma a partir de A ' anexándole la ma--

triz identidad RxR a -Z

N = -Z

IR

( 3 3 )

aquí N es una matriz estequiométrica para A , ya que el procedimiento Gauss- Jordan

esencialmente construye las columnas de Z en A ' para que se satisfaga

A Z = A C R ( 3 4 )

16

por tanto,

Z = A c' AR (35)

además de operaciones con renglones se requieren intercambios de col-umnas para obte-

ner la forma de matriz unitaria , que depende de como se hayan ordenado las especies

en el principio, como columnas de A .

Las etapas del procedimiento son las siguientes:

1.- Escribir la matriz de fórmula A para el sistema dado identificando cada columna por la especie que representa.

2.- Formar una matriz unitaria tan grande como sea posible en la por-- ción superior izquierda de A por medio de operaciones elementales con renglones e irtercambio de columnas (si es necesario). Si las columnas se intercambian, el orden de las especies también se debe intercambiar y el resultado es A ' .

3.- Posterior a estas etapas se establece que: a) Rango ( A ) = C, número de componentes que es igual al número de

cifras 1 en la diagonal principal de A ' . b) El conjunto de componentes está dado por las especies indicadas

sobre las columnas de la matriz unitaria. c) El número máximo de ecuaciones estequiométricas linealmente in-

dependientes esta dado por R = N - C . d) Se obtiene N en forma canónica de la submatriz Z de (32) de a-

cuerdo con (33), cada ecuación en un conjunto permisible de e- cuaciones químicas se obtiene de una columna N escribiendo

f 1 1 J i=l

A . v . . = O (26 b.)

y reordenando.

1.4 Ilustración del método.

Para el caso de la hidrólisis de ácido fosfórico, se sabe que el sistema con--

tiene las siguientes especies ( H Po H2Poi, HPO;, H+, OH-, H 0 ) y l o s ele--

mentos participantes son (H,O,P,k); donde k representa la carga iónica.

3 4 ' 2

Siguiendo el procedimiento con N = 7 y M = 3

17

1.- La matriz A es entonces;

3 2 1 0 1 1 2

A = 4 4 4 4 0 1 1

1 1 1 1 0 0 0

o -1 -2 -3 1 -1 o (1) (2) (3) (4) ( 5 ) ( 6 ) ( 7 )

donde cada columna corresponde a

el primer renglón corresponde a los hidrógenos presentes en cada molécula de acuer

do a como esth numeradas, el segundo renglón a los oxígenos, el tercero al fbsfo-

-

ro y el cuarto a la carga.

2.- Formando la matriz unitaria tan grande como es posible, tenemos A ' por medio

de operaciones entre renglones e intercambio de columnas;

A '

18

3.- a) C= rango (A) = 3; que es el número de elementos 1 en la diagonal princi-

pal.

b) Los componentes para el caso particular son H20 (7) , Poi3 (4) Y H+ ( 5 ) .

c) El númro de especies no componentes es R = 7-3 = 4 . Reordenando las es-

pecies según 10 especificado en A ' , la matriz estequiométrica corilpleta N

es ;

N =

2

-3

O

1

O

O

O

-1 1 -1

1 -2 1

o o -1

O 0 0

1 0 0

0 1 0

O 0 1

-Z

IR

donde I = RxR = 4x4 ; es una matriz identidad de 4x4. Observese que en

N hay 7 renglones después de aplicar la fórmula, cada uno corresponde a

una especie según su designación en A ' , y cada columna corresponde a una

reacción. Hay 4 reacciones o ecuaciones estequiométricas 'representadas

por las especies participantes, estas son;

3H2P04 """_ """_

H3P04 """_ """_

2H2P04 """_ -"""

OH- + H3P04 """_ """_

PO; + 2 H3POq

H+ + H2POi

HPO; + H3P04

H PO- + H20 2 4

19

1.5 Singularidad de las soluciones.

Antes de describir el cálculo y como lleva est0 a la construcción del algorit-

mo, es importante hacer notar algo referente a las soluciones.

Smith (1980) revisó las condiciones para la existencia de una solución para un

problema de equilibrio químico. Se supone que las restricciones de abundancia de

elementos y de no negatividad se satisfacen al menos para un vector de composición

n, y que todos los bk son finitos, además bk f O al menos para un elemento, también

se hace necesario que la función G sea continua en n.

No sólo es importante la existencia de la solución, sino también el número de

soluciones, es decir, cuántos vectores n satisfacen las restricciones de elementos

y las condiciones de equilibrio; 10 cual surge porque en algunos casos puede no ha-

ber singularidad. Esto se conecta típicamente a la formación incipiente de alguna

fase.

Para el caso del sistema formado por una sola fase de solución ideal, el pro--

blema de equilibrio tiene solución única. Una condición suficiente de singul.aridad

es G continua y estrictamente convexa de n . Para una sola fase, la convexidad de -

pende de la forma cuadrática de

donde 3 ’ G son los elementos de una matriz llamada Hessiana. ”

3 nian j

La singularidad se establece si‘ Q(&) > O para todas las composiciones permi-

tidas de n y las variaciones de dn que no desaparecen. Los elementos de la matriz

están dados por:

donde ij es la delta de Kroeneker

Rearreglando la ecuación tenemos

20

puesto que n. > O , Q es positivo a menos que la cantidad entre paréntesis sea cero

para cada i, pero como b es diferente de cero para al menos un valor de k, Q de k

be ser positivo y entonces se cumple la singularidad.

1

-

Para un sistema ideal de fases múltiples no es necesario que se cumpla la Sin-

gularidad corno emcontraron Hancock y Motzkin ( 1 9 6 0 ) .

11. PROBLEMAS NUMERICOS

Los procedimientos que se tratan para resolver el problema de equilibrio quími -

CO implican el considerar dos problemas numéricos; 1) La minimización de la función

f(x) sujeta a ciertas restricciones y 2) La solución de sistemas de ecuaciones alge -

braicas no lineales o trascendentales. Entonces se desea resolver

mín f(x)

X€ A.

donde a- es el conjunto de restricciones, g y x son N-vectores. Se supone que f y

g son diferenciables dos veces en forma continua. Cuando -h = R , las condicio- N

nes necesarias de primer orden para que f(x) presente un mínimo local en x* son

que este valor de x satisfaga las ecuaciones lineales ;

l%l X* = . O

si x* es una solución de (41), tambtifi es solución de la minimización.

Para que x* sea mínimo local de f, la ecuación (40) debe cumplirse y la ma--

triz Hessiana a'f/ axi ax sea positiva y deferenciable en x* . j 8

La mayoría de los algoritmos numéricos para resolver problemas lineales proce-

21

den por resolución de una secuencia de problemas cuyo grado de dificultad no es ma-

yor que el de resolver conjuntos de ecuaciones lineales. A cada etapa de esta se-

cuencia se le llama iteración. La secuencia termina cuando los componentes de x cam

bian de una iteración a otra en menos de cierta cantidad especificada.

El éxito depende de la complejidad de cada iteración y las propiedades de con-

vergencia. Así pues l o s algoritmos deben partir de una aproximación inicial x de

la solución y calcular la secuencia por medio de

O

m es un índice de iteración, w es el parámetro de etapa, lo que determina la dis--

tancia entre las iteraciones sicesLvas en la dirección definida por En gene-

ral los algoritmos utilizan (44) y difieren por la forma de calcular dx . (m)

La forma conveniente de determinar es encontrando el valor de w que mini-

miza aproximadamente a f(~(~)+ w&(~)) en cada iteración.

Todos los problemas de minimización son de la forma

mín G(x) (45)

y se ha visto que

es una función cuyo mínimo determina la solución de ecuaciones no lineales

g(x) = 0 (47)

Un procedimiento simple es cuando se sabe que w = 1 para proporcionar una esti-

mación del valor Óptimo, primero se calcula el valor

N

si esta cantidad es menor o igual a cero, todavía no se ha pasado el valor de minimi

zación de w y se produce la siguiente iteración, con w (m) = 1 en la ecuación ( 4 4 ) .

-

si la cantidad en (48) es mayor que cero se establece que

22

Y se emplea en los programas de aplicacibn general en computadora como en el algorit

mo VCS, del cual se hablará más tarde.

-

Si todo x de la solución de la ecuación (45) es positivo, se escoge w de mane- i

ra que todo valor de xi permanezca positivo.

2.1 Problema de minimización.

El problema de minimización que se considera, está restringido a la igualdad li -

neal y a la no negatividad

pin f(x)

X€ .t-L

donde es el conjunto de restricciones y esta definido por:

a menos que ( 3f/ a x . ) = O. El algoritmo que satisface la ecuacijn (48) se co-

noce como método descendente. 1 (m) X

2.2 Métodos de minimización.

Entre los métodos de minimización de funciones existen los métodos no restringi

dos de primero y segundo orden. Los métodos restringidos cuentan con el método de

los multiplicadores de Lagrange, además de otros, como el de la primera variación de

la función o de primer orden, el de segundo orden, análogo al de la segunda varia--

ción de la función, el método de proyección y el de conversión a problemas no res--

rativo de Newton - Raphson y el método de variación de parámetros.

111. ALGORITMOS

Una forma de clasificar los algoritmos es según como se utilizan las constan--

tes de restricción de elementos en los cálculos. Los algoritmos que eliminan estas

restricciones son los algoritmos estequiométricos. Estos tratan principalmente, el

número de variables independientes desconocidas como ( N' - M ) . As í mismo, se hace

referencia a algoritmos que tratan con las restricciones de abundancia de elementos

como lo hacen los algoritmos no estequiométricos, para los cuales el número de va--

riables es ( N' + M ) , para sistemas ideales se reduce a ( M + ) , donde es el

número de fases del sistena.

En este trabajo solo se hsrá referenci,a, principalmente, a l o s algoritmos este -

quiométricos.

3.1 Algoritmo estequiornétrico.

Como ya se vió, los números de moles n se relacionan a los grados de avance

O de conversión'de la reacción $ de las ecuaciones estequiométricas R, que son las

variables independientes por la ecuación (24), por tanto se puede escribir

G = G(T,P,§) (51)

Y el problema consiste en minimizar G , para valores fijos de T yP,en términos de

las RSi. Puesto que éstas ííltimas cantidades son independientes, las condiciones de

primer orden necesarias para un mínimo de G son

I a § /T,P,§ifj

con j = 1,2,...,R y hay R = N - C ecuaciones en el conjunto (52) puesto que N

donde

24

- I S 1 1 T,P,n - Pi kfi

considerándose también que;

donde vijpi = AG. y su valaor negativo es llamado afinidad de DeDonder. Las L J

las condiciones de equili-brio "clásicas". Cuando se introducen expresiones apropia--

das para p . en las ecuaciones en términos de energía libre y de número de moles, la

solución de éstas ecuaciones proporciona la composición del sistema en equilibrio.

1

A partir de una estimación de n de la solución de la ecuación (55), se obtienen

los números de moles en la siguiente iteración por medio de:

j =1

La expansión de la ecuación (55) en series de Taylor en n(m), despreciando los

términos de segundo orden y superiores e igualando a cero, da el Método de Newton-

Raphson, dando como resultado la siguiente ecuación

$ f Yij

donde j= 1 , 2 , ..., R . El: indice (m) indica la evaluación de n en la etapa m, para solucines ideales

es necesario introducir la expresión para el potencial químico

25

pi = p z + RT In i t

n

n -

de donde se deduce

sustituyendo(59) y (29

f d§im) 1=1

j= 1,2, ..., R además

- 1 "I n t en ( 5 7 ) se tiene

v = Yij - j i=l

(59)

La ecuación (44) se resulve para despejar y el resultado es utilizado en 1.a

ecuación (56) para determinar n. Este procedimiento es repetido hasta lograr la con -

vergencia.

Específicamente los algoritmos estequiométricos eliminan las restricciones de

abundancia de elementos, esto se lleva a cabo transformando las N incógnitas corres -

pondientes a los números de moles n que están restringidos por las M ecuaciones de

abundancia de elementos a un nuevo conjunto de variables de conversión de reacción,

iguales en número a R = (N' - M ) .

Los cambios en los números de moles (dm)) están relacionados a partir de cual -

quier estimación de que satisfaga la restricción de abundancia de elementos, a

las nuevas variables por

i=l

La matriz N tiene R = (N' - M) columnas linealmente independientes y está relacio--

nada a la matriz A por

26

\- N ’

k= 1,2, ..., M L ki IJ 1=1

a v . . = O j= 1,2, ..., R (63)

Considerando la función de Gibbs como una función de variables de reacción 5,

el. problema químico de equilibrio es minimizar G ( 5 ) . Las condiciones necesa;.ias pa-

ra esto son

lo cual equivale a

AG = N ~ ( $ 1 = O t (65)

En los algoritmos estequiométricos el número de variables es (NI- M) sin importar si

las fases son ideales, y para sistemas no ideales siempre se tienen menos variables.

Existen varios algoritmos que utilizan en la formulación la estequiometría de

las reacciones. Uno de los primeros algoritmos estequiométricos se debió a Naphtali-

(1959,1960,1961), quien sugirió utilizar mhtodos de primer orden; Villars diseno( 19-

59-1960) un algoritmo para resolver conjuntos de ecuaciones no lineales; Cruise(1964)

mejoró dicho algoritmo; Smith y Missen (1966 y 1968, respectivamente) mejoraron el

algoritmo Villars-Cruise como un método de minimización que mejora la convergencia.

Dentro de los algoritmos estequiométricos que utilizan la minimización de ener-

gía libre, existen l o s algori.tmos de primero y segundo orden. Naphtali (1959-1961) ,

sugirió el método de priemr orden para minimizar G ( $ ) . Las variables en cada itera--

ción son ajustadas por las cantidades d s , donde

los números de moles se ajustan por medio de la ecuación (62). Este algoritmo conver -

ge con bastante lentitud, especialmente cerca de la solución.

Los algoritmos de segundo orden como el de Hutchinson (1962) y otros, han suge-

27

rid0 la aplicación del método de Newton-Raphson a las ecuaciones (65), dando como re

sultado la siguiecte expresión:

-

lo cual requiere de la solución de un conjunto de R = (N' - M) ecuaciones lineales -

en cada iteración. Puesto que N' por lo general es grande comparada con M, la solu--

ción nemérica de estas ecuacior?es puede ser una parte del algoritmo que se lleva mu-

cho tiempo. Debido a la cuestiór de rapidez del cálculo, estos algoritmos no tienen

un uso muy difundido.

3.2 El algoritmo Villars-Cruise-Smith.

El algoritmo Villars-Cruise-Smith (VCS) está entre un algoritmo de primer orden

y uno de segundo orden, este algoritmo también se considera de estequiometria optimi -

zada. El algoritmo principia con la ecuación (67). Para una sola fase la matriz He--

ssiana está dada por

donde d es una delta de Kronecker. La ecuación (67) también se puede escribir como kj

7 N'

donde

v . = f Yki -

1

k=l ( 6 9 )

Haciendo uso de la libertad para seleccionar N de manera que la matriz Hessiana

28

si se logra que los primeros términos en la

tonces la Hessiana se invierte con facilidad

Es posible seleccionar N escogiendo [ Y . J

pecto al producto interno y a las normas de

de la ecuación (67) se invierta con facilidad, el uso de (67) resulta muy atractivo

ecuación (68) desaparezcan para ifj, en -

I de tal manera que sea ortonormal res-

los vectores partiendo de,

es decir,

Y

N' V kiVkj = d . .

1 J k=l n

r N'

vi J = 2 vkivkj . . n k=1

respectivamente.

El algoritmo VCS constituye básicamente un compromiso entre el cálculo de N en

cada iteración y el calcularla sólo una vez al principio del procedimiento. Observe

se que si la matriz es de la forma canónica, v ki~kj para ifj es cero cuando k se

refiere a una especie no componente (k M), ya que cada una de estas especies tie-

-

ne un coeficiente distinto de cero sólo en un vector estequiométrico. Cuando i = j,

V ki' kj v =1 para estos valores de k.

Los elementos de la Hessiana, numerando los componentes desde 1 hasta M y los

no componentes desde (M+1) hasta N, son expresados como

1 a 2 G = 1J d. . + 2 Vkivkj "

. - vivj RT 3Si 3Sj nj+M n i,j= 1,2, ..., R (73)

k=l k n t

escogiendo como especies componentes las que tengan los números de moles más altos,

el segundo miembro del lado derecho de la ecuación (72) tiende a ser pequeno y el -

primero grande, el último término desaparece si cualquiera de los ? o v desapare- i j

cen y en cualquier caso es pequeño debido a n en el denominador. t

29

Formando la matriz N de ésta manera, es posible suponer razonablemente que la

matriz Uessiana sea diagonal y se invierta directamente con facilidad obteniendose

así el algoritmo VCS consiste en utilizar la ecuación (67) con la ecuación (74) e

iterativamente se ajusta cada ecuación por una cantidad

/-I a c y

I RT

( 7 5 )

con j= 1 , 2 , ..., R.

En cada iteración los números de moles de cada especie son examinados para ve-

rificar que las especies componentes sean aquellas con los números de moles más al-

En el caso de un sistema ideal de fases múltiples la ecuación (75) se transfor -

ma en

(Siempre que al menos una especie para la que Y f O se encuentre en una fase de multicomponentes)

ij

RT

en caso contrario

aquí B representa la fase, dzR es la unidad si la especie k está en cualquier fase

de multicomponentes R y es cero en caso contrario, y es la unidad si la especie k

esta en la fase particular R de especies múltiples y es cero en caso contrario.

30

Se hace énfasis de que el algoritmo VCS es adecuado para problemas, espe-

cialmente para aquellos en que intervienen fases con una sola especie, como los que

surgen en aplicaciones metalúrgicas, y se debe al hecho de que las restricciones de

no negatividad en el número de moles de la especie son fácilmente manejadas en este

algoritmo.

En sus orígenes este algoritmo no se consideró como un método de minimización

de energia libre, sino más bien, como uno para resolución de ecuaciones no linea--

les representadas por las condiziones clásicas de equilibrio.

Villars propuso (1959-1960) ei uso de la ecuación (75) implicando una matriz N

seleccionada al azar. Ajusta cada ecuación estequiométrica individual por turno, y

recalcula la composición del sistema antes de ajustar la siguiente ecuación. Cruise

(1964) incorporó la selección optirnizada de N antes descrita y también se avocó al

ajuste simultáneo de todas las ecuaciones estequiométricas por medio de la ecua --

ción (75) en cada iteración antes de recalcular la composición del sistema, así mis

mo, encontró mejoras sustanciales en la rapidez del cálculo y convergencia al méto-

-

do de Villars, con estas dos modificaciones Smith (1966) y Missen (1968) reformula-

ron el método como un algoritmo de minimización e incorporaron w, el parámetro de e

tapa, y cuyo cálculo ya se analizó. Las modificaciones anteriores resultan en un al

goritmo rápido y libre de problemas de convergencia.

-

-

El algoritmo VCS es también llamado un algoritmo descendente de minimización ,

ya que para una sola fase ideal

y la ecuación (75) en cada iteración da como resultado

31

una ventaja importante del cálculo con este algoritmo es el hecho de que no hay e-

cuaciones lineales que resolver en cada iteración.

A continuación se presenta el diagrama de flujo del algoritmo VCS en la fig.1

ILeer e imprimir datos I

\1 Im- m+l [

-l# Verificar bases óptimas;calcular N

si es necesario L

Calcula a partir de la ec., ( 7 5 )

Calcular tamano del pará

metro de etapa w (m) - I L.

no

1 Imprimir resultados 1

fig.1

32

3.3 Formas de representar datos de energía libre estándar.

Cualquiera que sea el método que se utilice para calcular la composición de e-

quilibrio de un sistema, es necesario asignar un valor numérico al potencial quími-

co estándar p . . 1

Principalmente hay cuatro formas en que se dispone de información;

1) Como energías libres estándar de formación a partir de los elemen-

tos constituyentes (G ) , basadas en la convención de Raoult O de

Henry.

O

f

O O

2) Como valores de la función de energía libre, (G - H )/T O de la O O

expresión, ( G - H )/T. 298

O

3) Como entropías absolutas convencionales (S ) junto con entalpías -

O

4) Como potenciales estándar de electrodo (E ) .

3.3.1 Datos de energía libre a partir de información de,constantes de equilibrio.

Para la aplicaión particular del algoritmo a soluciones ideales es necesario

contar con datos de energía libre de formación confiables, casi siempre estos datos

se encuentran disponibles en forma de constantes de equilibrio por lo que es necesa

rio hacer la transformación a datos de potenciales químicos estándar p de las espe

ties individuales. Puesto que no es tan obvia la utilización de esta información en

i -

los algoritmos se tratará de describir a continuación, como obtener un conjunto con -

sistente de potenciales químicos estándar individuales para las especies participan -

tes.

Se inicia primero con el caso de una sola .ecuación estequiométrica (R = 1 ) ;

Yipi = o ( 7 9 )

i=l O

la constante de equilibrio esta relacionada a A G por medio de

O

AG = -RT 1nK ( 8 0 )

33

y G a su vez satisface la siguiente relación

N

i= 1

entonces lo que se desea determinar es el valor adecuado de p o a partir de K y A G O .

Se observa que los diferentes conjuntos de p . O dan lugar al mismo valor de @Go para

cualquier ecuación estequiométrica en que intervienen subconjuntos de las especies,-

i a

1

dando origen a soluciones idénticas a los problemas de equilibrio en que intervienen

dichas especies. Lo anterior significa que cualquier valor de p0 que satisfaga la e-

cuación (81) puede ser usado para el cálculo de la composición de equilibrio y ade--

l

más se dice que es consistente con el valor

En general se desea encontrar un valor

$- 1=1

y . .po 1J 1

donde se supone que R satisface R= rango(N)

de A G O .

p o que satisfaga

j= 1 , 2 , ..., R (82)

= N - rango ( A ) cada una de las ecuacio-

nes que se forman de la expresión anterior constituyen un conjunto de R ecuaciones

con N incógnitas 'con N > R y tienen un número infinito de soluciones.CUa1quiera de

éstas soluciones es un valor consistente de ,uo .

Pueden obtenerse soluciones particulares usando cualquier algoritmo adecuado pa -

ra resolver conjuntos indeterminados de ecuaciones lineales.

Ejemplo 3.

Para el caso de ácido fosfórico suponemos que se cuenta con datos en forma de

constantes de disociación;

H3P04 ====== -

H2P04 + H+ 9 K1

K2

34

aquí las K. son las constantes de equilibrio de las reacciones tal y corno están es--

critas, K es el producto ibnico del agua. Para encontrar el valor 1.1' consistente pa

ra las siete especies del sistema, primero se calcula AG..= -RT l a , y se forma la

1

W -

1 1

matriz ( N , - @,Go) t

K = 7.5~10-~ 1

K = 6.2~10 -8 2

K = 4.7~10 -13 3

K = 10 -4

4

-1 1 o

(Nt , A G O ) =

o o -1

o -1 1

O 0 0

la matriz unitaria correspondiente es

1 0 0 0

0 1 0 0

0 0 1 0

1 0 o o 1

-3

-2

-1

1

O 0

O 0

O 0

1 -1

-1 o

-1 o

-1 o

o -1

( 4 ) (7)

AGO = - 12.07364 KJjmol

Go = - 40.95275 ' I

Go = - 75. 72752 I '

A Go = - 79. 54630 ' I

-12.07364

-40.95275

-75.72752

-79.54630

35

Escribiendo las ecuaciones siguiendo el orden de las especies estipulado en la

matriz unitaria ( el cual esta indicado por los números entre paréntesis), se obtie -

nen las siguientes ecuaciones químicas:

3 H+ + poi3 """ """ GY. + a GZ +

z H+ + poi3 _____- - - - - - " H2POi n G; + n G; H3P04 GG

H+ 1- poi3 """ _""_ H POq2 GG

h20 ______ H+ + OH- - A G i """

se observa, por tanto, que el método establece potenciales químicos est.?idar o ener--

gias libres de formación para F< = 4 especies no componentes ( H PO HzPOi, HP04 y

OH-) , y les asigna el valor de cero a las especies componentes con C = 3 ( HzO, H.' y

-2 3 4'

3.3.2. Programa de computadora basado en el algoritmo VCS.

A continuación se presenta el programa compilado de calculo de concentraciones

de las especies en equilibrio basado en el algoritmo de minimización de energía li-

bre de Villars-Cruise-Smith.

El programa se denominó MINELI y esta escrito en lenguaje Fortran, la versión

que se presenta fue modificada para trabajar con solucionesideales en fase líquida,

las instrucciones de cambio también están dadas por los autores del programa. Este

programa es muy largo por lo que en secciones posteriores ,se darán algunas defini-

ciones de variables importantes , tanto para el programa, como para la formación del

archivo de datos.

A lo largo del programa se ven una serie de comentarios marcados con la letra

C que van dando una indicación del proceso del cálculo.

36

$LARGE PROGRAM MINELI IMPLICIT REAL*8(A-H,O-Z) COMMON/ST~lCH/SC(l5,149),FE(l~~),FF(lSO),W(lSO).DNG(l5O) 1,DNL(l5O),Wr(l5O),DG(l5O),DS(~~~),FEL(l5O),GA(l5),G~(l5),TGl 2,TLl,TG,TL,DTG,DTL,BM(150,1~),T,P,TING,TINL,IF,IEST,DA(lSO) 3 , l N D ( l 5 O ) , S P ( 3 , 1 5 O ) , I C ( l ~ ~ ) , I R ( 1 5 0 ) , ~ ~ ( ~ ~ ~ ) , M , N , N C , N R , ~ , N E INTEGER SI,SP MAXIT=500 CALL VCS(l,I,O,MAXm STOP END SUBROUTINE VCS(IRP,IPR.IPl,MAXIT)

C ALGORITMO DE EQUILIBRIO ESTEQUIOMETRICO EN SISTEMA IDEAL C USANDO EL METODO VCS C EL PROGRAMA PUEDE MANHARE HASTA 150 ESPECIES Y 15 ELE- C "ENTOS, PERO LAS ESPECIFICACIONES DIMENSIONALES PUEDEN C VARIARSE. C MAXIT ES EL NUMERO MAXIM0 DE ITERACIONES PERMISIBLES SI NO SE C ALCANZA LA CONVERGENCIA,CASI SIEMPRE EL NUMERO DE ITERACIONES ES C MENOR DE 100,PERO MUY CERCA DE LOS LIMITES DE UNA FASE DE ESPECIES C MULTIPLES PUEDEN REQUERIRSE MILES DE ITERACIONES. C SE PUEDE MANEJAR CUALQUIER NUMERO DE FASES DE UNA ESPECIE C Y DOS FASES DE ESPECIES MULTIPLES LA FASE 1 ES C NOMINALMENTE UN GAS, YA QUE ALOG(P)SE SUMA A LOS DATOS DE C POTENCIAL QUIMICO ESTANDAR, ESTO PUEDE EVITARSE HACIENDO c p = l . c LA FASE DOS ES NOMINALMENTE IJN LIQUIDO, o CUALQUIER C FASE PARA LA CUAL LOS DATOS DE POTENCIALFS QUIMICOS C ESTANDAR SON INDEPENDIENTES DE P. C LA FASE DE ESPECIES MULTIPLES ESTARA AUSENTE SI C NT.LT.1.E-10. SI LA FASE DE ESPECIES MULTIPLES ESTA C AUSENTE EN EL EQUILIBRIO, EL VALOR DG/RT SE REFIERE A

C EN EL EQUILIBRIO ESTABLECIDO. C DEBE PROPORCIONARSE UNA RUTINA DE PROGRAMACION LINEAL C PARA LA OPCION DE LA ESTKMACION INICIAL DE LA COMPOSICION C DE EQUUBRIO (VER SUBRUTINA INESQ

C I-SIGMA(X(I)),DONDE X@ SON LAS FRACCIONES MOL VIRTUALES

IMPLICIT REAL*8(A-H,O-Z) REAL*8 LOUT,LIST COMMON/STOICH/SC(l5,149),FE(150),FF(150),W(l5O),DNG(l50) l,DNL(l5O),WT(l5O),DG(l5O),DS(l5O),FEL(1SO),GA(l5),GAI(l5),TGl 2,TLl,TG,TL,DTG,DTL,BM(I50,I5),T,P,TING,TINL,IF,IEST,DA(150) 3,lND(150),SP(3,150),IC(l5O),IR(150),SI(l5O),M,N,NC,NR,MR,NE INTEGER SI, SP,XM,STATE,STATl ,STAM INTEGER OUT(2),VF(13) LOGICAL LINDEP,SOLDEL,IFIRST,CONV,ICONV LOGICAL DAT,LIQ,LEC,ITL,ILT,LBO,IM,MACH,FORCED COMMON/FORCER/LIQ DIMENSION PS(l5O),WX(lSO),E(l5),MJ(15),EA(15),TI(8O),FFO(lSO), lLIST(30),LIST1(14) DIMENSION GA"(25),DB(lO),V(lO),PHY(lO) DIMENSION LOUT (10) DIMENSION SA(l5),SS(l5),SM(l5,15) DIMENSION YY(15,l) DIMENSION YY(15)

DIMENSION XY(150) DIMENSION X(16), Z(16), B(16) DIMENSION AW(lSO), AA(15.15) DIMENSION BB(15), PSOL(15O),DSOL(I5),CC(150),RW(623),I(303) DIMENSION AX(17,152) EQUIVALENCE (X,YY)

EQUIVALENCE (FFO,FF),(PS,WX,MJ,DS) EXP(X)=DEXPO() ALOG(X) = DLOG(X) A B S 0 = D A B S 0

EQUIVALENCE @(I) , YY(1))

?

37

DATA OUT/'SI(','I) ' / DATA OUT2/'SI(I) ' 1 DATA ~[~~~/'2','3','4','S','6','7','8~,'9','10','11','12', *'13','14','15'/ DATA VF/'(lX,'.'3A4,','1X,',' ','F6.3',',15',',T40',',D12',

*'.5,6','X,D1','2.5.','D12.'.'5) ' 1 DATA LIST/'5H213','SH42X','SH313','5H39X','5H413','5H36X','SH513', 1'5H33X','5H613','5H3OX','5H7~','~H2~','SH8I3','SH24X','5H9l3', 2'5H21X','5H1013','5Hl8X','5HllI3','~Hl~X','~Hl2l3','~Hl2X', 3'SH13I3','5H9X','5Hl4I3','5H6X','5HlS~','SH3X','SHl6~','5HlX'/ DATA LOUT/'SH(','SH IX','SH','3A4','SH','5H 1PD1','5H 2.5', 1'5H 2X','SH D12.','5HS)'/ DO 9191 J= 1,150 FE(J)=O. DA(J)=O. IF(IRP.NE.O)W(J)- 1 .DO FEL(J) =O.

9191 WT(J)=O.DO C C INSTRUCCIONES DE FORMATO DE ENTRADA C 101 FORMAT(713) 102 FORMAT(3A4,1X,lSF2.0,11,F9.3) 1022 FORMAT(3A4,IX,Il,F9.3) 2022 FORMAT(lSF5.3) 103 FORMAT(8E10.4) 901 FORMAT(8OAI)

902 FORMAT(//IS,BH SPECIES,I8,9H ELEMENTS,E16,11H COMPONENTS/ C FORMATOS DE SALDA

115,15H PHASE1 SPECIES,IlO,lSH PHASE2 SPECIES,I8, 222H SINGLE SPECIES PHASES//9H FRESSURE,F22.3,4H ATM/ 312H TEMPERATURE,F19,3,2H K)

903 FORMAT(14H PHASE2 INERTS.Fl7.3) 904 FORMAT(14H PHASE1 INERTS,F17.3) 905 FORMAT(/IX,40H MODIFIED LINEAR PROGRAMMING ESTIMATE OF,

906 FORMAT(/lX,29H USER ESTIMATE OF EQUILIBRIUM) 907 FORMAT(/8H SPECIES,SX,lSH FORMULA VECTOR,lOX 1,15H STAN.CHEM.POT.,3X,17H EQUILIBRIUM EST./)

9070 FORMAT(/BH SPECIES,SX,ISH FORMULA VECTOR,T40 1,15H STAN.CHEM.POT.,3X,17H EQUILIBRIUM EST./)

91 11 FORMAT(I2IH ELEMENTAL ABUNDANCES,lJX,8H CORRECT,lOX, 114H FROM ESTIMATE //(26X,A2,1PD20.12,D20.12))

911 FORMAT(I2lH ELEMENTAL ABUNDANCES,SX,A2,1PD20.12/(26X,A2,D20.12)) 9 12 FORMAT( 15A2) 913 FORMAT(23H VCS CALCULATION METHOD/) 914 FORMAT(/lX,80Al)

112H EQUILIBRIUM)

C FUAR VALORES INICIALES, LEER DATOS,IMPRIMIR TFULOS IF(IRP.EQ.O)NRUNS= I NID = O IF(IRP.EQ.0) GO TO 620 READ(1,lOl)NRUNS

620 DO 9999 U= 1 ,NRUNS ICONV= .FALSE. IF(IRP.EQ.O)GO TO 623 READ(l,lOl)M,NE,NSl,NL1,IF,IEST,ICE NC=NE LF(ICE.LT.O)GO TO 621 READ(l,lO2)((SP(J,I),J=l,3),(BM(I,J),J=l,1S),SI(I),FFO(I),I=1,M) GO TO 622

621 READ(l,1022)((SP(J,I),J= 1,3),SI(I),FFO(I),I= 1,M) READ(I,2022)((BM(1,J),J=l,lS),I=l,M)

622 IF(IEST.LT.O)GO TO 293 READ(1,103)(W(I),I=I,M)

293 READ(l,103)(GAI(T),I=l,NE) READ(l,lO3)T,P,TING,TINL READ(I,912)(E(I),I=I,NE)

38

READ(I,901) TI

NOPT=O IM=.FALSE.

623 l T = O

C C FORMATOS DE LAS VARIABLES INICIALES C

LOUT(4)=LIST(2+NE-I) LOUT(S)=LIST(2+NE) VF(4)=LISTI(NE-I)

C C M P R M I R TITULOS

WRITE(4,913) WRlTE(4,914)TI

C C CALCULAR NUMERO Y TIPO DE FASES C

NS=O NL=O DO 855 I= l ,M INDO = I IF(SI(I).EQ.O)NS=NS+ 1

855 IF(SI(I).EQ.2)NL=NL+ 1 IF(NS.NE.NSl)WRITE(4,105) IF(NL.NE.NLl)WRITE(4,106)

105 FORMAT(49H NUMBER OF SINGLE-SPECIES PHASES DOES NOT COMPUTE) 106 FORMAT(42H NUMBER OF PHASE2 SPECIES DOES NOT COMPUTE)

NG=M-NL-NS IF(NG.NE.I)GO TO 852 NS=NS+I NG=O DO 856 I= 1 ,M IND(l)=I IF(SI(I).NE.I)GO TO 856 IF(IF)8561,8562,8563

GO TO 856

GO TO 856

SI(I) = o

8561 FFO(I)=FFOO+ 1.9872DOfT*ALOG(P)

8562 FFO(I)= FFO(I) + ALOG(P)

8563 FFO(I)=FFO(T)+.0083143DO+T+ALOG(P)

WRITE(4,107)(SP(J,I),J=1,3) 856 CONTINUE 852 IF(NL.NE.I)GO TO 854

NS=NS+ I NL=O DO 858 I= l ,M IF(SI(I).NE.2) GO TO 858

WRlTE(4,108)(SP(J,I),J=1,3) SI(I)=O

858 CONTINUE 107 FORMAT(14H THIS SPECIES:,3A4,19H IS THE ONLY GAS.IT

108 FORMAT(14H THIS SPECIES:,3A4,22H IS THE ONLY LIQUIDAT

854 N=M-NC

137H WILL THEREFORE BE TREATED AS A SOLID)

137H WILL THEREFORE BE TREATED AS A SOLID)

LIQ=NL.GT.O NR=N MR=M NCI =NC+ 1 DO 889 I= l ,M IR(I)=NC+I

889 CONTINUE TEST=-I.D-lO IF(IEST.GE.O)GO TO 30001 DO 30000 I= 1 , h 4 R

3oooO W(I)=-FFO(I)

39

TEST=-1 .E20 30001 CONTINUE C C DETERMINAR EL NUMERO DE COMPONENTES C

CALL BASOPT(.TRUE..NOFT,AW,SA,SM,SS,AA,TEST,IT,CONV) IF(IPR.NE.O)WRlTE(4,902)M,NE,NC,NG,NL,NS,P,T NC1 =NC+ 1 IF(IPR.NE.O.AND.NG.GT.O)WRlTE(4,904)TING IF(IPR.NE.O.AND.LIQWRlTE(4,903)TINL

C C AJUSTAR PARA EL TIPO DE DATOS DE POT.QUIM.ESTANDAR C 910 RT= 1.9872DO+T

IF(1F) 165,163,164 165 TF= 1000.DORT 169 FORMAT(28H STAN.CHEM.POT.IN KCAL./MOLE)

164 TF=l.DO/(O.O083143DO+T) 168 FORMAT(28H STAN. CHEM. POT. IN N/MOLE)

163 TF=O. 167 FORMAT(23H STAN.CHEM.POT.IS MU/RT) 166 DAT=TF.NE.O

DO 888 I= l ,M IR(I)=NC+I DA(I)= FFO(I) IF(DAT)FFO(I)= FFO(I)+TF IF(SI(Q.NE.1)GO TO 8866 FFO(I) = FFOO + ALOG(P) IF (NG.EQ.1) SI(I)=O

GO TO 166

GO TO 166

8866 IF(SI(I).EQ.O)FEO=FF(I) 888 IF(SI(I).EQ.Z)SI(I)=-l

C C EVALUAR ESTIMACION DE EQUILIBRIO SI SE REQUIERE C

IF (1EST.GE.O)GO TO 159 CALL INEST(L,.FALSE.,NOPT,AW,SA,SM,SS,AA,TEST,IT,NCl)

494 CALL ELAB 11591 CALL ELCORR C C IMPRIMIR DATOS DE ESTIMACION DE EQUILIBRIO

159 CONTINUE IF(IEST.GE.O)CALL ELAB IF(IPR.EQ.0) GO TO 1955 WRJTE(4,911 I)(E(I),GAI(I),GA(I),I= 1 ,NE) lF(IEST.LT.O)WRRE(4,905) IF(IEST.GE.O)WRlTE(4,906) IF(IF.LT.O)WRlTE(4,169) IF(IF.EQ.O)WRlTE(4,167) IF(lF.GT.O)WRlTE(4,168) IF(ICE.LT.O)GO TO 78 WRlTE(4,907) WRlTE(4,890)(E(I),I= 1 ,NE),OUT

DO 865 I= 1 ,M IXA=SI(D

DO 860 J = 1,NE 860 MJ(J)=BM(I,J) 865 wRITE(4,LOVn(SP(J,r),J=1,3),(MI(n.J= l,NE),IXA,DA(I),W(T)

GO TO 1955 78 CONTINUE

890 FORMAT(' ',13X,17A3)

IF(SI(I).EQ.-l) R A = 2

WRlTE(4,9070) WRRE(4,891)(E(I),l= l,NE),OUT2

891 FORMAT(' ',13X,17(1X,A5))

40

DO 77 I = l , M IXA=SI(I) IF(SI(I).EQ.-l)EA=2

77 WRlTE(4,VF)(SP(J,I),J= 1,3),(BM(I,J),J= l,NE),IXA,DA(I),W(I) C WRlTE(4,VF)W(I) C C EVALUAR TOTAL DE MOLES DE GAS Y DE LIQUIDO C 1955 TG=TING

TL=TINL DO 1 I=l,MR IF (SI(I))3,1,2

03 TL=TL+ W(I)

02 TG=TG+W(I) O1 CONTINUE

G O T 0 1

EVALUAR TODOS LOS POTENCIALES QUIh4ICOS

CALL DFE(W,O,O) IF(IT.GT.0) GO TO 430 IF(IEST.LT.O)GO TO 441

DETERMINAR ESPECIE BASE,EVALUAR EsTEQUIOMETRIA

429 CALL BASOPT(.FALSE.,NOFT,AW,SA,SM,SS,AA,TEST,IT,CONV)

6714 WRITE(4,6713)

6713 FORMAT(44H CONVERGENCE TO NUMBER OF POSITIVE N(I) LESS

24192 El= 1

C C EVALUAR VECTOR MAYOR-MENOR INICIAL C 323 JM=O

DO 140 I=I ,NR L= R(I) IF(W(L).GT.O) GO TO 400 lM=JM+ 1 IC(I)=-l GO TO 140

400 IC(I)= 1 140 CONTINUE

IF(.NOT.CONV)GO TO 24192

GO TO 857

17H THAN C/35H CHECK RESULTSTO FOLLOW CAREFULLYII)

ILT= .FALSE.

LEC=.FALSE. LBO= .TRUE. IF(lT.GT.O)GO TO 96

C C EVALUAR TODOS LOS CAMBIOS DE ENERGIA LIBRE DE REACCION C 430 CALL DELTAG(0,l ,NR)

IF(.NOT.LBO)GO TO 212 1430 LBO=.FALSE.

IF(.NOT.IM)GO TO 93

GO TO 90 9090 I T I = O

C C ESTABLECER VALORES INICIALES PARA LA ITERACION C EVALUAR AJUSTES DE LAREACCION C

93 ITI=4*(IT1/4)-ITl IF(lTI.NE.O)GO TO 90 CALL DFE(W,O, 1) CALL DELTAG( 1,l ,NR)

90 DO i707 I= 1 , M R 1707 FEL(I) = FE(I)

41

IF(ILT.OR.IM)GO TO 11608 CALL ST2(SOLDEL,*429)

LEC= .FALSE. DTG=O. IF(LIQ)DTL=O. DO 1092 I= 1,NC

11608 ITL= .FALSE.

1092 DS(T)=O.

C C CICLO PRINCIPAL DEL CALCULO A PARTIR DE AQUI A LA INSTRUC.NO.33 C

9058 L=IR(I) 9059 IF(IC(1))134,34,22

I= 1

134 IF(DG(I).GE.O.)GO TO 3331 JM = JM- 1 rc(9 = 1 IM= .FALSE. ILT= .FALSE. GO TO 3331

C C ESPECIES MENORES C

34 IF(ITI.NE.0) GO TO 333 DGG = DG(9

IF(DGG.LT.82,)GO TO 36 WW=W(L)

1373 IF(WW.LT.l.E-26)GO TO 1374 W(L) = WW* 1 . E-6 GO TO 1135

36 IF(ABS(DGG).LE.l .E-5)GO TO 3331 C = A L O G O - D G G

W(L) = EXP(C) IF(WT(L).LE.W(l))GO TO 1135 WT(L)=WW*l.D+6

DS(L) = DX

GO TO 18

IF(C.LE.(-73.6))GO TO 1373

1135 DX=WT(L)-WW

IF(WW.LE.(-DX))ITL=.TRUE.

22 IF(ABS(DG(I)).LE.l.E-6)GO TO 3331 DX = PS(L) WT(L) = W(L) + DX

C C VERIFICAR MOLES POSITIVAS EN ESPECIES MAYORES ,

c IF(WT(L).LE.O.)GO TO 85

C C CALCULAR CAMBIOS EN EL NUMERO DE MOLES PARA GAS LIQUIDO c

18 DO 14 K=I,NC 14 DS(K)=DS(K)+SC(K,I)*DX

DTG = DTG + DNG(I)*DX IF(LIQ)DTL=DTL+DNL(I)*DX GO TO 33

333 1 W(L) = W(L) 333 DS(L)=O 33 I=I+l

1232 IF(I.LE.NR)GO TO 9058

C c L I M ~ A R REDUCCION DE ESPECIES’BASE A 99 POR CIENTO C 1587 P+=S

DO 1086 I= 1,NC XX=-DS(I)/W(I)

P

42

1086 IF(PAR.LT.XX) PAR=XX PAR= 1 ./PAR IF(PAR.LE.1.01 .AND.PAR.GT.O.)GO TO 29 PAR= l . LK= NC GO TO 1589

C C REDUCCCION EN BASE DEMASIADO GRANDE C REDUCIR TAMANO DE ETAPA GENERAL C

29 PAR= .99*PAR ITL= .FALSE. DO 1088 I= 1 ,MK

DTG = DTG*PAR IF(LIQ)DTL=DTL*PAR ITL= .FALSE. LK= MR

1088 IF(DS(I).NE.O.)DS(I)=DS(I)*PAR

1589 DO 1591 I = 1,LK 1591 WT(I)=W(I)+DS(I)

TG 1 =TG + DTG IF(LIQ)TLl=TL+DTL GO TO 95

C C MOLES NO POSITIVAS DE LAS ESPECIES MAYORES C

C C FIJAR ESPECIES DE FASES PURAS EN CERO C

85 IF(SI(L).NE.O)GO TO 86

DX=-W(L) 1123 DO 1119J=l ,NC

WT(J)=W(J)+SC(J,I)*DX IF(WT(J).LE.O.) GO TO 1120

DO 1124J=I,NC 1 124 W(J) = hT(J)

GO TO 1125 1120 DX=DX/2. .

1 1 19 CONTINUE

GO TO 1123 1 125 W(L) = W(L) + DX

TG=TG+DNG(I)*DX IT‘ ‘-IQ)TL=TL+DNL(I)*DX

IC(r)=-l I r N(L).GT.O.)GO TO 333 I

CALL DFE(W,O,ITI) CALL DELTAG(IT1,l ,NR) DO5111 LL=l,MR

5 1 1 1 FEL(LL)= FE(LL) GO TO 1218

C C CAMBIAR MAYOR AMENOR C 1118 IC(I)=O 1218 JM=JM+I

IM=JM.EQ.NR IF(IM.AND.ITI.NE.O.)GO TO 227 IF(SI(L).EQ.O)GO TO 3331 GO TO 18

c C DESCARTAR ESPECIES MENORES INFERIORES A I .E-32 MOLES C 1374 MM=O

CALL DELETE(L,0,*1366) 11375 JM.=JM-I

GO TO 1232 C

43

C CORTAR AJUSTE DE REACCION A 113 CPARA MOLES POSITNAS DE LA ESPECIE MAYOR C

86 DX=-.9*W(L) DS(L) = DX W(L) = W(L) + DX IF(W(L).LT..Ol*W(NC))GO TO 1 1 18 GO TO 18

C C POTENCIALES QUIMICOS TENTATNOS C

95 CALL DFE(WT, 1 ,IT0 c C IMPRIMIR RESULTADOS INTERMEDIOS C

IF(IP1 .EQ.O.)GO TO 28616 LI=ITI-1 WRITE(4,7167)IT,Ll

128H EVALUATION OF STOICHIOMETRY.13) 7167 FORMAT(/lOH ITERATION ,I3,22H ITERATIONS SINCE LAST

WRITE(4,4731) 4731 FORMAT(8H SPECIES,19X,13HWITIALMOLES,6X,llHFINAL MOLES,

18X,6H MU/RT,16X,lOH DELTAWRT) WRrrE(4,802)((SPcJ,I),J=1,3),W(I),WT(I),FE(I),I= 1,NC)

8055 FORMAT(lX,3A4,13X,IPD14.7.3D19.7) DO 4732 I=NCl,MR LI=I-NC

4732 WRITE(4,8055)(SP(J,I),J= 1,3),W(l),WT(I),FE(r),DG(LI) 286 16 CONTINUE

FORCED= .FALSE. IF(IM.OR.ILT.OR.lTL)GO TO 18617

C C FORZAR CONVERGENCIA C

CALL FORCE(ITI.FORCED) C C RESTABLECER VALORES AL FINAL DE LA ITERACION C CALCULAR LOS NUEVOS CAMBIOS EN LA ENERGIA LIBRE DE LA REACCION C 18617 lT=IT+l

IT l=ITI+I 142 IF(ITI.EQ.O.)CALL DELTAG(O,l,NR)

IF(lTI.NE.O.)CALL DELTAG(-I,l,NR) IF(F0RCED)GO TO 44 14 TG=TGI IF(LIQ)TL=TLl DO 10 I= 1,MR

10 IF(DS(I).NE.O.)W(I)=WT(I) C C FIJAR LA FASE DE ESPECIES MULTIPLES MICROSCOPICA(NT.LT.1 . € - l o ) C EN CERO C 4414 CONTINUE

XM=O IF(NG . EQ.O.OR.TG.EQ.O.0R.TG.GT. 1 .E- 1O)GO TO 995 1 DO 9952 I= 1 , M R IF(SI(I).NE. l)GO TO 9952 W(I)=O. W(I) =o. DS(I)=O.

9952 CONTINUE TG=O. TGl =O. DTG=O. XM= 1

9951 1FrNL.EQ.O.OR.TL.EQ.O.OR.TL.GT.l.E-IO) GO TO 9554 DO 3953 I= 1,MR

44

lF(SI(I).NE.(-1))GO TO 9953 w(l)=o. W ( I ) = O . DStI)=O.

9953 CONTINUE TL=O. TLl =O. DTL = O. XM= 1

9554 IF(XM.EQ.O)GO TO 441 DO 9956 I= 1 ,NR L=IR(I)

9956 IF(W(L).EQ.O.)lC(I)=-l CALL BASOPT(.FALSE.,NOPT,AW,SA,SM,SS,AA,TEST,IT,CONV) CALL DFE(W,O,O) CALL DELTAG(0,l ,NR) IF(.NOT.CONV) GO TO 441 WRlTE(4,9993) GO TO 6714

9993 FORMAT(I32H DELETION OF MULTISPECIES PHASE.) C C VERIFICAR SI LA BASE ES OPTIMA C 441 CONTINUE

IF(NC.EQ. 1 .)GO TO 4442 DO 4441 I=2,NC IF(W(1-i).LT.W(I))GO TO 4442

DO 426 I= 1 ,NR

J=NC

4441 CONTINUE

L = IR0

427 IF(W(L).LE.W(J))GO TO 426 IF(SC(J,I).NE.O.)GO TO 429 J=J-1 GO TO 427

426 CONTINUE

4442 DO 4443 I = 1, NR

DO 4443 J = 1 ,NC IF(W(L).LE.W(J))GO TO 4443 IF(SC(J,I).NE.O.) GO TO 429

GO TO 1241

L= IR(0

4443 CONTINUE 1241 IF(lT.EQ.O.)GO TO 24192

C C REEVALUAR EL VECTOR MAAYOR-MENOR,SI ES NECESARIO C

IF(ITI.NE.O.)GO TO 212

DO 3731 I= l ,NC 3731 WX(I)=.Ol*W()

1920 L=IR(I) I=1

IF(W(L).GT.O.)GO TO 21 1 IF(DG(I)) 20,20,1191

21 1 IF(SI(L).EQ.O.)GO TO 13 DO 1 lK= 1,NC I F ( S I g < ) . N E . O . A N D . S C ( K , I ) . N E . O . . A N D . W ( L ) ) G O TO 13

IF(SI(L).EQ.l.AND.W(L).GT..IO*TG)GO TO 13

221 JM=O

11 CONTINUE

IF(Sl(L).EQ.(-l).AND.W(L).GT..lO*TL)GO TO 13 IC(I)=O. GO TO 119

1191 IC(I)=-l 119 JM=JM+ 1

GO TO 19 13 IF(IC(I).GT.O.)GO TO 19

?

45

20 IC(I) = 1 IF(ITI.EQ.O.)GO TO 19

C C POT.QUIM.,DELTAG PARA EL NUEVO MAYOR C

IF(W(L).GT.O.)GO TO 1216 FE(L) = FF&) GO TO 1919

1216 IF(SI(L))1213,1919,1215 1213 FE(L)=FF(L)+ALOG(W(L)/TL)

1215 FE(L)=FF(L)+ALOG(W(L)/TG) 1919 CALL DELTAG(O,I,I)

GO TO 1919

19 I = l + l IF(I.LE.NR)GO TO 1920 IM=JM.EQ.NR

212 IF(IM) GO TO 227 C C VERIFICAR EQUILIBRIO PARA LAS ESPECIES MAYORES C 222 DO 15 I=I,NR ff(IC(I).LE.O.OR.ABS(DG(I)).LE.I.E-O6)GO TO 15 IF(IT.LT.MAXIT)GO TO 3857 lCONV= .TRUE. GO TO 857

3857 ILT= .FALSE. GO TO 93

15 CONTINUE C C VERIFICAR EQUILIBRIO PARA LAS ESPECIES MENORES C

IF(JM.EQ.O.)GO TO 96 227 IF(lTI.EQ.O.)GO TO 228

CALL DFE(W,O,I) 228 CALL DELTAG(I,l,NR)

DO 315 I=I,NR

IF(lT.LT.MAXIT)GO TO 2857 ICONV = .TRUE. GO TO 857

lF(IC(T).NE.O.OR.ABS@G(I)).LE.I.E-O5)GO TO 315 '

C ' C CONVERGENCIA ENTRE LAS MAYORES PERO NO EN LOS MENORES C 2857 ILT= .TRUE.

IF(ITI.EQ.0.) GO TO 90 GO TO 9090

3 15 CONTINUE C C EVALUAR VALORES FINALES O INTERMEDIOS C DE LAS ABUNDANCIAS DE LOS ELEMENTOS C 96 CALL ELAB

IF(.NOT.LEC)GO TO 1368 C C VERIFICAR EXACTITUD DE LA ABUNDANCIA DE LOS ELEMENTOS C

DO 1362 I=I ,NC IF(GAI(T).EQ.O.ODO)GO TO 1367

GO TO 1362 IF(ABS(GA(I)-GAIO).GT.SD-O8*GAI(I)) GO TO 1363

1367 IF(ABS(GA(I)).GT.I.D-1O)GO TO 1363 1362 CONTINUE

GO TO 1365 1368 LEC= .TRUE.

C C CORREGIR LAS ABUNDANCIAS DE LOS ELEMENTOS C

4

46

1363 CALL ELCORR GO TO 1955

1365 IF(MR.EQ.M)GO TO 857 C C VOLVER A VERIFICAR LAS ESPECIES DESCARTDAS C 1366 MRl =MR+l

NR1 =NR+ 1 DO 849 I=MRl,M

CALL DELTAG(O,NRl,N) I = N R I NPB=O. X T I =o. XT2=0. IF(TG.GT.O.)XTl=ALOG(l.E+32*TG) IF(LIQ.AND.TL.GT.O.)XT2=ALOG(l.E+32*TL)

IF(SI(L))850,850l,8501

GO TO 18531

849 FEO = FF(I)

1968 L=IR(I)

8501 IF~G.GT.O.AND.DG(I).LT.XTI.OR.TG.EQ.O.AND.DG(I).LT.O.) GOT0 853

850 IF(TL.GT.O.AND.DG(I).LT.XT2.OR.TL.EQ.O.AND.DG(l).LT.O.) GO TO 853 18531 I=1+1

1FO.LE.N)GO TO 1968 IF(NPB)857,857,1969

C C REINSERTAR ESPECIES DESCARTADAS C 853 W(L)=l.E-20 IF(Sl(L))8531,8533,8532

GO TO 8533 8531 TL=TL+ W(MR)

8532 TG=TG+W(MR) 8533 IC(I)=O. JM=JM+l NR=NR+ 1 MR=MR+l NPB = 1 CALL DELETE(L, 1 ,*1366)

C CALL DELETE(L,I,MRl) GO TO 18531

1969 ET= .TRUE. CALL DFE(W,O, 1) CALL DELTAG(O,NRl,NR) GO TO 9090

C C VALORES FINALES DELTAGIRT C 857 J=MR+l K=NR+ 1 DO 847 I = 1 ,NR J=J-1 K=K-1

847 DG(J) = DG(K) C C EVALUAR FRACCIONES MOL FINALES C

DO 1810 I=l,MR IF(SI(I).NE.O)GO TO 7153 IF(W(I).EQ.O.)WT(I)=O. IF(W(I).GT.O.)WT(I)= 1.

7153 LF(SI~.EQ.l.AND.W(I).NE.O.)WT(I)=W(I)iTG 1810 IF(SI(I).LT.O.AND.W(r).NE.O.)WT(I)=W(I)iTL C C ORDENAR ESPECIES DEPENDIENTES EN ORDEN DECRECIENTE C

DO 8475 I = NC 1 , M R

47

R(I) = I

DO 8476 L= NC I .MR CALL AMAX(XY .K,L,MR) IF(K.EQ.L)GO TO 8476 CALL DSWO(Y,K,L) CALL SWITCH(IR,K,L)

IF(IPR.EQ.0.) RETURN

8475 XY(I) = W(I)

8476 CONTINUE

C C ESPECIFICACIONES DE FORMATO DE SALIDA C 801 FORMAT(] 1H ITERATION=,I5/15H EVALUATIONS OF

115H STOICHIOMETRY=,I4//1X,7HSPECIES,In<,l7HEQUILIBRIUM MOLES 2,4X,13HMOLE FRACTION,4X,14HDG/RT REACTION0

802 FORMAT(IX,3A4,13X,IPD14 7,5X,D14.7,5X,D11.4) 803 FORMAT(I16H G/RT=,IPD15.7/22H TOTAL PHASE1 MOLES= ,ID11.4) 804 FORMAT(IX,3A4) 805 FORMAT(IX,3A4,13X,1PD14.7,5X,D14.7) 806 FORMAT(22H TOTAL PHASE2 MOLE=, 1 PDI 1.4)

L

C IMPRIMIR RESULTADOS c

WRlTE(4,5 138) 5138 FORMAT(///)

IF(IC0NV) WRITE(45139) WRlTE(4,801)IT,NOPT

5139 FORMAT(36H CONVERGENCE CRITERION NOT SATISFIED) 1809 WRITE(4.805)((SP(J,I),J= 1,3),W(I),WT(T),I= 1,NC)

DO 18990 I=NCI,MR L=IR(I)

18990 WRlTE(4,802)(SP(J,L),J=I,3),W(L),WT(L),DG(L) IF(MR.EQ.M)GO TO 1873 WRITE(4,9074)

MRI =MR+ 1 WRlTE(4,804)((SP(J,I),J=1,3),1=MRl,M)

IF(TG.GT.O.O.AND.TING.GT.O.O)G=G+TING*ALOG(TING/TG) IF(LIQ.AND.TINL.GT.O.O)G=G+TINL*ALOG(TINL/TL) DO 1847 I= 1 ,h4R

1847 G=G+W(I)*FE(I) 374 WRITE(4,803)G,TG

9074 FORMAT(/23H LESS THAN 1 . E 3 2 MOLES0

1873 G=O.

IF(NL.GT.O.)WRITE(4,806)TL WRlTE(4,91I)(E(I),GA(I),I= 1,NE) IF(NRUNS.NE. 1)RETURN

C C REORDENAR DATOS DE ENTRADA C C DO 1375 I= l ,M

C DO 1378 J = I , M DO 1378 J=I,M

L=IND(J) K1 = J IF(L.EQ.I)GO TO 1379

1378 CONTINUE 1379 CALL DSW(W,I,KI)

CALL DSW@A,I,KI) CALL DSW(WT,I,Kl) CALL DSWZ(SP,I,Kl) CALL SWITCH(SI,I,Kl) CALL SWITCH(IND,I,Kl) DO 1376 J = l , N E

1376 CALL DSW(BM(I,J),I,Kl) 1375 CONTINUE DO 1377 I= l , M lF(SIg.EQ.(-l))SI(r)=2

40

1377 FF(I)=DA(I) 9999 CONTINUE

RETURN END

C C ELIMINA O ADICIONA UNA ESPECIE A LOS CALCULOS

SUBROUTINE DELETE(L,MM,*)

COMMON/STOICH/SC(l5,149),FE(150),FF(150),W(l5O),DNG(l50) l,DNL(l5O),WT(l5O),DG(l5O),DS(l5O),FEL(15O),GA(l5),GAI(l5),TGl 2,TLl,TG,TL,DTG,DTL.BM(150,15),T.P,TING,TINL.IF,lEST,DA(150) 3 ,1ND(lSO) ,SP(3 ,15O) , IC(15O) , iR(15O) ,SI ( ,NE

IMPLICIT REAL*8(A-H,O-Z)

INTEGER SP,SI

COMMONlFORCEIULIQ LOGICAL LIQ

SE INTRODUJO AQUI COMMONlFORCER AUNQUE NO ESTA EN LA FUENTE

16700 NNL=L-NC IF(L.EQ.MR)GO TO 16796

C C REORDENAR DATOS CUANDO SE ADICIONA O ELIMINA UNA ESPECIE C

CALL DSW(WT,MR,L) CALL DSW(W,MR,L) CALL DSWZ(SP,MR,L) CALL DSW(FF,MR,L) IF(SI(MR).EQ.O)CALL DSW(FE,MR,L) CALL SWlTCH(SI,MR,L) CALL DSW(FEL,MR,L) CALL SWlTCH(IC,NR,NNL) CALL DSW(DA,MR,L) CALL SWITCH(IND,MR,L) DO 16701 K= l,NC T1 =SC(K.NR) SC(K,NR)=SC(K,NNL)

16701 SC(K,NNL)=Tl CALL DSW(DNG,NR,NNL) IF(LIQ.CALL DSW(DNL,NR,NNL) DO 3030 J=l ,NE CALL DSW(BM(l,J),MR,L)

3030 CONTINUE 16796 IF(”.NE.O.)RETURN L

C PROCEDIMIENMTOS ADICIONALES AL ELIMINAR UNA ESPECIE C

CALL DSW(DG,NR,NNL) CALL DSW(DS,MR,L)

L

C PROCEDIMIENTOS ESPECIALES PARA DESCARTAR UNA FASE L

IF(SI(MR)) 16798,16797,16799

GO TO 16797 16798 TL=TL-W(MR)

16799 TG=TG-W(MR) 16797 NR=NR-1

MR=MR-1 IF(NR.EQ.O.)RETURNI RETURN END

SUBROUTINE FORCE(IT1,FORCED) C C FORZADOR DE CONVERGENCIA C

MPLICIT REAL*8(A-H,O-Z) COMMONISTOICHISC(l5,149),FE(l5O),FF(l5O),W(l5O),DNG(l5O) 1,DNL(l5O),WT(l5O),DG(l5O),DS(l5O),FEL(l5O),GA(l5),G~(l5),TGl

L

49

2.TL1,TG,TL,DTG,DTL,BM(150,1~).T,P,TING,TlNL,IF.lEST,DA(150) 3.IND(150).SP(3,150),1C(150),IR(1~~),Sl(150),M,N,NC,NR,MR,NE INTEGER SI.SP COMMON/FORCER/LlQ LOGICAL LIQ.FORCED

C C CALCULAR PENDIENTE AL FINAL DE LA ETAPA c 13700 S2=0.

13702 IF@S(I).NE.O.)S2=S2+FE(I)*DS(I)

C C CALCULAR PENDIENTE ORIGINAL

DO 13702 I= I .MR

IF(S2.LE.O.)RETURN

S1 =o. DO 13704 I= 1,MR

13704 IF(DS(I).NE.O.)Sl = S 1 +FEL(I)+DS(I) IF(S1 .GE.O.) RETURN

C C AJUSTAR PARABOLA C

AL=SI/(SI-SZ) IF(AL.GE..95)RETURN

C C AJUSTAR NUMERO DE MOLES,F'OT.QUIM C

13705 IF@S(I).NE.O.)W(I)=W(I)+AL+DS(I) DO 13705 I= 1 , M R

TG=TG+AL*DTG IF(LIQ) TL=TL+AL+DTL CALL DFE(W,O,ITI) FORCED= .TRUE. RETURN END

SUBROUTINE ST2(SOLDEL,*) C C CALCULA AJUSTES DE LA REACCION C

IMPLICF REAL+8(A-H.0-Z) COMMON/STOICH/SC(l5,149),FE(l5O),FF(l5O),W(l5O),DNG(l5O) l,DNL(I5O),WT(150),DG(150),DS(l5O),FEL(ISO),GA(15),GAI(I5),TG1 2,TL1,TG,TL,DTG,DTL,BM(I50,I5),T,P.TING,TlNL,~F,lEST.DA(lSO) 3,lND(l5O),SP(3,15O),IC(l5O),IR(l5O),Sl(l5O),M,N,NC,NR,~,NE INTEGER S I S P COMMON/FORCER/LIQ LOGICAL LIQJOLDEL A B S 0 = D A B S 0

14700 SOLDEL=.FALSE. DO 14701 I= l ,NR

IF(W(L).EQ.O.AND.Sl(L).NE.O.)GO TO 47101

IF(IC(I).LE.O.AND.DG(I).GE.O.)GO TO 14701 IF(SI(L).NE.O)GO TO 14703 s=o.o GO TO 1471 1

14703 S = 1 ./W(L) 14711 DO 14704 J= l ,NC 14704 IF(SI(J).NE.O.) S=S+SC(J,I)**2/W(J)

L=IR(I)

IF(ABS@G(I)).LE. 1 .D-6)GO TO 14701

IF(TG.GT.0.) S=S-DNG(I)*+2/TG IF(TL.GT.0.) S=S-DNL(I)++2/TL IF(S.EQ.O.)GO TO 47103

GO TO 14701 DS(L)=-DG(I)/S

C

50

c REACCION EXCLUSIVAMENTE ENTRE FASES CONDENSADAS C DESCARTAR UN SOLIDO Y RECALCULAR LA BASE C 47103 CONTINUE

IF(DG(T).LE.O.)GO TO 47104 DSS = W(L) K = L DO 47106 J= l ,NC lF(SC(J,I).LE.O.)GO TO 47106 XX = W(J)/SC(J,I) IF(XX.GE.DSS) GO TG 47106 DSS = XX K=J

47106 CONTINUE DSS=-DSS GO TO 47109

47104 DSS=l.E+10 DO 47105 J= l ,NC IF(SC(J,I).GE.O.)GO TO 47105 XX=-W(r)/SC(J,r) IF(XX.GE.DSS)GO TO 47105 DSS = XX K = J

47105 CONTINIJE 47109 CONTINUE

IF(DSS.EQ.O.)GO TO 14701 W(L)=W(L)+DSS DO 47108 J = l , N C

47108 W(J)=W(J)+DSS*SC(J,I) w(K)=o. 1FN.NE.L) SOLDEL= .TRUE. GO TO 14701

C C FASE DE ESPECIES MULTIPLES CON N T - O C 47101 IF(DG(I).GE.-l.E-4) GO TO 47102

DS(L)=l.E-lO IC(0 = 1 GO TO 14701 .

47102 DS(L)=O. 14701 CONTINUE

IF(S0LDEL) RETURN 1 RETURN END

SUBRUTINA PARA CORREGIR ABUNDANCIA DE ELEMENTOS

SUBROUTINE ELCORR IMPLICIT REAL*8(A-H,O-Z) COMMON/STOICH/SC(1S,149),FE(150),FF(150),W(l5O),DNG(l50) I,DNL(l5O),WT(l50),DG(150),DS(150),FEL(150).GA(15),GAI(l5),TG1 2,TL1,TG,TL,DTG,DTL,BM(150,15),T,P,TING,TI[NL,I,l~T,DA(150) 3,lND(150),SP(3,15O),IC(l5O),IR(l5O),Sl(l5O),M,N,NC,NR,~,NE INTEGER SI,SP DIMENSION X(l6),XY(l5O),AA(l5,15),YY(15,1) DIMENSION X(16),XY(l5O),AA(15,15),YY(15)

EQUIVALENCEO<,YY) EQUIVALENCE(X(I).YY(l))

18700 DO 18701 I=I ,NC X(l)=GA(I)-GAI(T) CALL UNPACK(BM,I,NC,XY) DO 18701 J= l ,NC

CALL MLEQU(AA,lS,NC,YY,l)

DO 18702 I= I ,NC

18701 AA(J.r)=XY(J)

PAR,= .5

XX=-X(I)/W@)

51

18702 IF(PAR.LT.XX) PAR=XX PAR = 1 ./PAR IF(PAR.LE. 1 .O1 .AND.PAR.GT.O.)GO TO 18729 PAR= I . GO TO 18759

18729 PAR = .99*PAR 18759 DO 18704 I = 1 .NC 18704 W(T) = W(r) + PAR"X(0

RETURN END SUBROUTINE DELTAG(L,J,KP)

C C CALCULA LOS CAMBIOS DE ENERGIA DE REACCION C

MPLICIT REAL*8(A-H,O-Z) COMMON/STOlCH/SC(l~,149),FE(l5O),FF(150),W(lSO),DNG(lSO) 1.DNL(l5O),Wr(l5O),DG(l5O),DS(l5O),FEL(l5C),GA(l5),GAI(l5),TGl 2,TLl,TG.TL,DTG,T)~L.BM(I50,15),T,P,TING,'I'!NL,IF,IEST.DA(150) 3,1ND(150),SP(3,150),1C(150),IR(150),S1(150~,M,N,NC,NR,MR.NE INTEGER SI,SP EXP(X) = D E X P O

15700 IF(L) 15710,15711,15712 C C SOLO MAYORES C 15710 DO 15701 K= 1,NR

IF(IC(Kj.EQ.0.) GO TO 15701 LL= IR0 DG(K) = FE(LL) DO 15713 LL=l,NC

15713 DG(K)=DG(K)+SC(LL,K)*FE(LL) 15701 CONTINUE

GOT0 15715 C C TODAS LAS REACCIONES C 1571 1 DO 15702 K=J,KP

LL= IR&) DG(K) = FE(LL) . DO 15702 LL= 1,NC

GO TO 15715 15702 DG(K)=DG(K)+SC(LL,K)*FE(LL)

C C SOLO MENORES C 15712 DO 15703 K=l,NR

IF(IC(K).GT.O.)GO TO 15703 LL = IR(K) DG(K) = FE(LL) DO 15714 LL=I,NC

15714 DG(K)=DG(K)+SC(LL,K)*FE(LL) 15703 CONTINUE c C ESPECIES DE FASES MULTIPLES CON NT=O c 15715 CONTINUE

SDELl =O. SDEL2=0. DO 15716 I=l ,NR LL= R(r) IF(W(LL).NE.O.)GO TO 15716 IF(DCi(r).GTSO.)DG(r)=50.

IF(SI(LL).EQ.I) SDELl =SDELl +EXP(-DG(I)) IF(SI(LL).EQ.-I) SDEL:!=SDEL:!+EXP(-DG(r))

DO 15717 I = 1,NR

IF(DG(r).LT.-50.)DG(I)=-50.

15716 CONTINUE

52

LL= R(I) IF(W(LL).NE.O.)GO TO 15717 IF(SI(LL).EQ. l)DG(I)= 1 .-SDELI IF(SI(LL).EQ.-l)DG(I)= 1 .-SDEL2

15717 CONTINUE RETURN END

C C C

SUBROUTINE BASOPT(IFIRST,NOPT,AW,SA,SM,SS,AA,TEST,I,CONV)

SELECCIONA BSE OFTIMA,CALCULA FSTEQUIOPMETRIA

IMPLICIT REAL*8(A-H,O-Z) INTEGER SP,SI,XM COMMONISTOICHISC(l5,149),FE(150),FF(150),W(l5O),DNG(l50) l,DNL(l5O),WT(l5O),DG(l5O),DS(l5O),FEL(15O),GA(l5),GAI(l5),TGl 2,TLl,TG,TL,DTG,DTL,BM(150,15),T,P,TING,TINL,IF,IEST,DA(150) 3,IND(l5O),SP(3 ,15O),IC(15O),IR(15O),SI( ,MR,NE LOGICAL IFIRST,LINDEP LOGICAL C O W DIMENSION AW(l5O),SA(l5),SM(l5,15),SS(15),AA(l5.15) A B S 0 = DABSO() CONV=.FALSE.

12700 NOPT=NOPT+ 1 DO 12701 I = l , M R

12701 AW(Ij=W(I) JR=O

12750 JR=JR+l C C DETERMINA NUMERO DE MOLES RESTANTE MAS GRANDE C 12796 CALL AMAX(AW.K,JR,MR)

IF (AW(K).EQ.O.)CONV=.TRUE. IF(C0NV)RETURN IF(AW(K).NE.TEST)GO TO 12759 NC=JR-I N=M-NC NR=N DO 12760 I= 1 ,M

GO TO 12758 12760 IR(l)=NC+I

12759 CONTINUE AW(K)=TEST

C C VERIFICAR INDEPENDENCIA LINEAL CON LAS SPECIES ANTERIORES C C LOGICAL FUNCTION LINDEP(BM,JR.K,NC) C

JL=JR-I SA(JR)=O. CALL UNPACK(BM,K,NC,DG) DO 13748 J= l ,NC

13748 SM(J,JR)=DG(J) IF(JL.EQ.O.)GO TO 13731 DO 13758 J = l , J L SS(J)=O. DO 13749 I= 1,NC

13749 SS(J)=SS(J)+SM(I,JR)*SM(I,J) 13758 SS(J)=SS(J)/SA(J)

DO 13746 J = l , J L DO 13746 L= 1,NC

13746 SM(L,JR)= SM(L,JR)-SS(J)+SM(L,J) 13731 DO 13747 ML=l,NC 13747 IF(ABS(SM(ML,JR)).GT.l.E-I‘T)SA(JR)=SA(JR)+SM(ML,JR)**2 C

53

C SI NORMA DEL NUEVO RENGLON.LT.lO*+-3,RECHAZAR

IF(SA(JR).LT.l .E-06) GO TO 13790 LINDEP= .FALSE. C O T O 12791

13790 LINDEP=.TRUE. C C END OF LINDEP C C 12791 IF(L1NDEP) GO TO 12796 C C REORDENAR DATOS

CALL DSW(W,JR,K) CALL DSW(WT,JR,K) CALL DSWZ(SP,JR,K) CALL DSW(FF,JR,K) CALL DSW(FE.JR,K) CALL DSW(AW,JR,K) CALL SWlTCH(SI,JR,K) CALL DSW(DA,JR,K) CALL DSW(FEL,JR,K) CALL SWITCH(IND,JR,K) DO 2929 J = 1 ,NE CALL DSW(BM(1,J)JR.K)

2929 CONTINUE 12758 IF(JR.LT.NC) GO TO 12750

C C EVALUAR EsTEQUIOMETRIA C

IF(1FIRST) RETURN

DO 12754 J= l ,NC CALL UNPACK(BM,J,NC,DG) DO 12754 I= 1,NC

12754 AA(I,J)=DG(I) DO 12753 I= l ,N CALL UNPACK(BM,IRO,NC,DC) DO 12753 J= 1.NC

12753 SC(J,I)=DG(J) CALL MLEQU(AA,lS,NC,SC,N) DO 12755 I= l,N K= IR0

C C EVALUAR VALORES DELTA N C

IF(SI(K)) 12763,12762,12761 12762 DNG(I)=O.

GO TO 12765 12763 DNL(I)= I .

DNG(I)=O. GO TO 12764

12761 DNGO= 1 . 12765 DNL(I)=O. 12764 DO 12755 J = 1 ,NC

IF(ABS(SC(J,I)).LE.l.E-06) SC(J,I)=O. rF(sr(n) 12756,12755,12757

12756 DNL(I)=DNL(I)+SC(J,I)

12757 DNC(I)=DNG(I)+SC(J,I) 12755 CONTINUE

RETURN END

GO TO 12755

SUBROUTINE ELAB C C CALCULA ABUNDANCIA DE ELEMENTOS

54

C IMPLICIT REAL*8(A-H,O-Z) COMMON/STO~CH/SC(l5.149),FE(l5O),FF(l5O),W(l5O),DNG(l5O) 1,DNL(I5O),WT(I50),DG(150),DS(150).FEL(l5O),GA(l5),G~(l5),TGl 2.TL1,TG,TL,DTG.DTL,BM(150,15),T,P,TING,TINL,IF.IEST,DA(150) 3,IND(150),SP(3,150),IC(150),IR(150),SI(l5O),M,N,NC,NR,MR,NE INTEGER SI,XM,SP

17700 DO 17702 I=l ,NE 17702 GA(I)=O.DO

DO 17703 I= 1,NE DO 17703 I= l,M

17703 GA(J)=GA(J)+BM(I,J)*W(I) RETURN END

C C SUBRUTINA PARA EVALUAR LOS POTENCIA C

SUBROUTINE DFE(Z.KK.LL) IMPLICIT REAL *8(A-H,O-Z)

,LES QUD VllCOS

COMMON/STOICH/SC(l5,149),FE(l5O),FF(l5O),W(l5O),DNG(l5O) 1,DNL(150),~(150),DG(150),DS(150),FEL(150),GA(15)~GA1(15),TG1 2,TLl,TG,TL,DTG,DTL,BM(I50,15),T,P,TING,TINL,IF,IEST,DA(150) 3,IND(l5O),SP(3,150),IC(l5O),IR(ISO),SI(l5O),M,N,NC,NR,MR,NE INTEGER S1,SP COMMONIFORCER/LIQ LOGICAL LIQ DIMENSION Z(1)

C SE CAMBIO ESTA LINEA, PERMITIENDO QUE Z TENGA NC ELEMENTOS C DIMENSION Z(150) C DIMENSION Z(NC)

A L O G O = D L O G O IF(KK.GT.O)GO TO 402 X =TG IF(LIQ) Y =TL GO TO 403

402 X=TGl IF(L1Q) Y =TL1

403 IF(X.GT.0.) X=ALOG(X) IF(.NOT.LIQ) GO TO 404 IF(Y .GT.O.) Y = ALOG(Y)

Ll =NC GO TO 56

404 IF(LL.EQ.0) GO TO 55

55 L1 =MR C C TODAS LAS ESPECIES O COMPONENTES C 56 DO 1 I=I ,L1

IF(Z(I).NE.O.)GO TO 1 1 11 FE(I) = FF(I) G O T 0 1

1111 1F(S1(I))11,1,14 C I 1 FE(I)=FF(I)-Y +ALOG(Z(I)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C INSTRUCCION ADECUADA PARA SOLUCIONES

11 FE(I)=FF(I)+ALOG(Z(I))-ALOG(Z(l)*O.O18016DO)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IF (l.EQ.I)FE(I)=FF(I)+ALOG(Z(T))-Y

GO TO 1 14 FE(I) = FF(I) + ALOG(Z(I))-X O1 CONTINUE

IFCL) 10,99,12 C C SOLO MAYORES C

10 DO 2 I= 1,NR

55

IF(IC(I).EQ.O.)GO TO 2

IF(Z(L).NE.O.)GO TO 11 12 FE(L) = FF(L) GO TO 2

11 12 IF(SI(L))21,2.23

L = R(T)

C 21 FE(L)=FF(L)-Y +ALOG(Z(L)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C INSTRUCCION ADECUADA PARA SOLUCIONES

21 FE(L)=FF(L)+ALOG(Z(L))-ALOG(Z(l)*O.O18016DO) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

GO TO 2 23 FE(L) = FF(L) + ALOG(Z(L))-X 02 CONTINUE

RETURN L

C SOLO MENORES

12 DO 3 I = 1,NR IF(IC(I).NE.O.)GO TO 3

IF(Z(L).NE.O.)GO TO I I13 FE(L) = FF(L) GO TO 3

1 113 IF(SI(L))3 I ,3,32

L=IR(I)

C 31 FE(L)=FF(L)-Y +ALOG(Z(L)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C INSTRUCCION ADECUADA PARA SOLUCIONES

31 FE(L)=FF(L)+ALOG(Z(L))-ALOG(Z(l)*O.O18016DO) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

GO TO 3 32 FE(L) = FF(L)+ ALOG(Z(L))-X 03 CONTINUE 99 RETURN

END

SUBROUTINE UNPACK(BM,K,NC,X) DIMENSION X(15),BM(150,15) DOUBLE PRECISION X,BM DO 3 L=l ,NC

03 X(L) = BM(K,L) RETURN END

SUBROUTINE MLEQU(C,IDEM,N,B,M)

DIMENSION C(IDEM,IDEM),B(IDEM,M) DO 58 I= l ,N IF(C(I,I).NE.O.)GO TO 51 I P l = I + l DO 53 K=IPI ,N IF(C(K,I).NE.O.)GO TO 54

53 CONTINUE WRITE(435)

55 FORMAT(19H NO UNIQUE SOLUTION) C CALLEXIT

54 DO 57 J=I,N 57 C(I,J)=C(I,J)+CO<,J)

DO 67 J=l ,M 67 B(I,J) = B(1,J) + B(K,J) 51 DO 58 L=I ,N

R=C(L,I)/C(I,I)

IMPLICIT REAL*8(A-H,O-Z)

IF(L.EQ.l.OR.C(L,I).EQ.O.) GO TO 58

DO 62 J = l , N

DO 64 J=l ,M 62 C(L,J)=C(L,J)-C(I,J)*R

64 B(L,J)=B(L,J)-B(I,J)*R

56

58 CONTINUE DO 63 I=I.N DO63 J=I .M

63 B(I,J)=-B(I,J)/C(I,I) RETURN END

C C C

SUBRUTINA PARA REORDENAR LOS ELEMENTOS VECTORES

SUBROUTINE SWITCH(X,Jl,J2) INTEGER X(l),T T=X(JI) X(Jl)=X(J2) X(J2)=T RETURN END

SUBROUTINE AMAX(X,K,J,N) DIMENSION X(1) DOUBLE PRECISION X,T K=J BIG = X(J) DO 1 I=J,N IFO<(l).LE.BIG)GO TO 1 K=I BIG=X(I)

O1 CONTINUE

C

RErURN END

SUBROUTINE DSW(X,JI ,J2) DOUBLE PRECISION X,T DIMENSION X(1) T=X(JI) X(J 1) =X(J2) X(J2)=T RETURN END

SUBROUTINE DSWZ(X,Jl,J2) INTEGER X,T DIMENSION X(3,l) T=X(I ,JI) X(l,JI)=X(l,J2) X( 1 ,J2) =T T=X(Z,JI) X(2,JI)=X(2,12) X(2.J2)=T T=X(3.JI) X(3,JI)=X(3,12) X(3,J2)=T RETURN END

SUBROUTINE INEST(L.IFIRST,NOFT,AW,SA,SM,SS,AA,TEST,IT,NCI)

C ESTIMA LA COMPOSICION DE EQUILIBRIO C

IMPLICIT REAL *8(A-H,O-2) COMMON~STOICH/SC(l5,149),FE(l5O),FF(l5O),W(l5O),DNG(l5O) ~,~~~(~~~),Wr(l~~),DG(l5O),DS(l5O)~FEL(15O),GA(l5),GAI(l5),TGl 2,TL~,TG,TL,DTG,DTL,BM(150,1S),T,P,TING,T~L,~,l~T,DA(lSO) 3,IND(~50),SP(3,150),IC(1~0),n\(lSO),SI(1SO),M,N,NC,NR,MR,NE INTEGER SP,SI,XM COMMON/FORCER/LIQ

57

LOGICAL LINDEP,IFIRST,LIQ,CONV DIMENSION BB(lS) ,PS~L( lSO) .DSOL(lSO~,CC(lSO) ,RW(6~3) ,~(3~~) 1,AX(17,1S2),AW(150),SA(lS),SM~~~,~~~~SS~l~~.AA~l~.~~~ 2.FFO(lS0) A L O G O = DLOGP) EXP(X) = DEXP(X) ABS@)= DABS@)

11700 K L = O . LT=O. T L l = O .

C C FOJAR DATOS PARA RUTINA DE PROG. LIN. C

MI = O . M2= NE IA=17 DO 11703 I= 1.NE

DO 11701 J = l , M CALL UNPACK(BM,J,NE,DG) DO 11701 I= 1 ,NE

11701 AX(I,J)=DG(T) DO 11702 J = 1,M

11703 BB(I)=GAI(I)

11702 CC(J)=-FFQ C C UTILIZAR RUTINA DE PROG. LIN C C CALL ROUTINE TO SOLVE MAX(CC*PSOL)SUCH THAT AX*PSÓL=BB C

DO 11706 I= l ,M DS(I)=O. W(l) = PSOL(I) IF(PSOL(Q.EQ.0.) W(I)=l.D-lO

11706 CONTINUE

1 1707 WT(I) = W(T)

C C CALCULAR TOTAL DE MOLES GASEOSAS Y LIQUIDAS C POTENCIALES QUIMICOS DE LA BASE C 11717TGI=TING

DO 11707 I=I ,M

CALL BASOPT(.FALSE.,NO?'T,AW,SA,SM,SS,AA,TEST,I",CON~

IF(LIQ TLI =TINL DO 11711 I= l ,NC IF(SI(I).GT.O)TGI =TG1 +WT(I)

DO 11712 I= 1,NC IF(SI(I).EQ. 1 ) FE(I)=FF(l)+ ALOG(WT(I)/TGI)

DO 11725 I=NCl,M

CALL DELTAG(0, I ,N)

1171 1 IF(SI(D.LT.O)TLl =TLI +WT(I)

11712 IF(SI(I).LT.O)FE(I)=FF(I)+ALOG(WT(I)/TLl)

1 1725 FE(I) = FF(I)

C C ESTIMAR AJUSTES DE REACCION C 11726 DTG=O.

IF(L1Q) DTL=O. IF(TG1 .EQ.O.) GO TO 11753 XTI=ALOG(I.E+32+TGl) XT2=ALOG(I.E-32*TGI)

11753 I F v L l .EQ.O.)GO TO 11754 XT3=ALOG(I.E+32*TLl) XT4=ALOG(I .€-321TTLI)

11754 DO 11718 I=I ,NR L= lR(9 IF(SI(L)) 11721,11718,11720

11720 IF(DG(I).GE.XTI .OR.DG(I).LE.XT2)GO TO 11718

58

DS(L)=TGl*EXP(-DG(I)) GO TO 11733

11721 IF(DG(I).GE.XT3.OR.DG(I).LE.XT4) GO TO 11718

11733 DO 11724 K=l ,NC 11724 DSO<)=DS(K)+SC(~,n*Ds(L)

DTG=DTG+DNG(I)*DS(L) IF(LIQ)DTL=DTL+DNL(I)*DS(L)

DS(L)=TLl*EXP(-DG(T))

11718 CONTINUE C C MANTENER LA ESPECIE DE BASE POSITNA C

PAR= .S DO 11786 I=I ,NC

PAR= 1 . IPAR IF(PAR.LE.1 .AND.PAR.GT.O.)GO TO 11729 PAR= l . GO TO 11788

11729 PAR=O.B*PAR C C CALCULAR NUEVOS NUMEKOS DE MOLES C 11788 DO 11794 I=l ,NC 11794 W(I)=W(I)+PAR*DS(I)

DO 11791 I=NCI.M 11791 IF(DS(I).NE.O.)W(I)=OSO+PAR

TG=TGl+DTG*PAR IF(LIQ) TL=TLI +DTL*PAR IF(LT.GT.0.) RETURN

11786 IF(PAR.LT.(-DS(l)/WT(I))) PAR= -DS(l)/W(I)

C C SECCION PARA FORZAR CONVERGENCIA C

CALL DFE(W,O,O) s=o.o DO 11731 I=l .MR

IF(S) 11773,494,11735 1 173 1 IF(DS(I).NE.O.)S=S + DS(I)*FE(I)

11773 IF(KL)494.494,11738 11735 IF(KL.GT.O)GO TO 11738 C C PROBAR LA MITAD DEL TAMA;O DE ETAPA C

S1 =S PAR=.S*PAR KL= 1 GO TO 1 1788

C C MUSTAR PARABOLA POR LA MITAD DE LA ETAPA C 11738 X L = .5*(1 .-S/(Sl-S))

C C MALA DIRECCION,REDUCR TAMA;O DE ETAPA A.2

IF(XL.GE.O.)GO TO 17138

c PAR = PAR* .2 GO TO 11740

17138 IF(XL.LE.1 .)GO TO 11742 L

C ETAPA DEMASIADO GRANDE,CONSIDERAR ETAPA ORIGINAL COMPLETA C

PAR = PAR*2. GO TO 1 1740

C C ACEPTAR RESULTADOS DE FORZAR CONVERGENCIA C 11742 PAR=PAR*2.*XL

59

11740 LT= 1 GO TO1 1788

494 RETURN END

60

para la formación del archivo de datos.

NRUNS.- número de problemas que han de correrse.

M .- es el símbolo para el número de especies

NE.- número de elementos

NL1.- número de especies en la fase 2 ( fase líquida)

SP.- nombre de la especie

S1 .- tipo de fase, es decir, si es líquida o gaseosa siguiendo la especificación:

O una sola fase

1 multiespecies fase gas

2 multiespecies fase líquida

FF0.- potencial químico para las especies

IF .- tipo de datos del potencial químico; -1 en Kcal/mol

O en MU/RT

1 en KJ/mol

IEST .- estimación inicial; O si la da el usuario -1 si la hacé la máquina

ICE.- tipo de vector fórmula; O para elementos enteros

-1 para los no enteros

BM .- Vector fórmula, se forma con el número de átomos del elemento que hay en la fór -

mula molecular de la especie presente.

W(1) .- Vector de estimación inicial de las especies, formado por el número de moles iniciales de las especies participantes.

I

GAI(1) .- Vector de abundancia de elementos.

T .- temperatura en grados Kelvin

P .- presión dada en atmósferas.

TING .- total de moles de gas inerte presente.

61

T1NL.- total de moles de liquido inerte.

E(1) .- nombre del elemento representado por su símbolo o con cierto número de carac -

teres especificados en el formato de entrada.

TI .- título arbitrario

Cuando se corren problemas que involucran iones, las cargas se leen como uno de

los elementos y deben ser colocadas al final del vectro fórmula.

Si el usuario hace su propia estimación inicial, debe satisfacer la ecuación

BM (1,J) * W(1) = GAI (J) ( 8 3 )

lo que significa que se debe de cumplir el balance de materia (vector fórmula por la

estimación inicial debe ser igual al vector de abundancia de elementos).

En la llamada a la subrutina VCS(IRP,IPR,IPl,MAXIT) , los primeros tres argumen -

tos controlan la entrada y salida. Si IRP = O no se leen datos de otra forma debe

ser igual a 1. Si IPR = O no se imprimen datos ni resultados, de otro modo IPR = 1 ,

si se imprimen resultados. Si IP1 = 1 se imprimen resultados intermedios en cada ite

ración, MAXIT es el número de máximo de iteraciones permitidas si no se alcanza la

-

convergencia.

3.5 Aplicación del algoritmo. Ejemplo para la formación del archivo de datos.

Aunque el algoritmo en general es utilizado para mezclas de gases, también es u

sad0 para el balanceo de ecuaciones redox, en el caso particualr se trata con solu--

ciones ideales, el tratamiento de transformación de datos de constantes de equili--

brio a potenciales estándar para este tipo de sistemas ya se trató en una sección an

terior. Siendo consistentes este ejemplo también tratará sobre ácido fosfórico.

-

-

Para utilizar el algoritmo en este tipo de sistemas, en la subrutina DFE del al -

goritmo se deben realizar los siguientees cambios;

La instrucción 11 FE(I)= FF(I)-Y+ ALOG(Z(I)) se sustituye por

11 FE(1) = FF(I)+ ALOG(Z(1))- ALOG(Z(l)*O.Ol8OlSDO)

IF(I.EQ.1) FE(1) = FF(I)+ ALOG(Z(1))- Y

62

la instrucción 21 FE(L) = FF(L)- Y+ ALOG(Z(L)) es sustituida p o r

21 FE(L) = FF(L)+ ALOG(Z(L))- ALOG(Z(l)*O.O180l6DO)

por Último la instrucción 31 FE(L) = FF(L)- Y+ ALOG(Z(L)) se reemplaza por

31 FE(L) = FF(L)+ ALOG(Z(L))- ALOG(Z(l)*O.Ol8016DO)

La entrada de datos debe tener el siguiente formato

1

7,4,0,7,1,0,-1

HZ0 2 -79.54630

H 2 0.00000

PO4 2 0.00000

OH 2 0.00000

H3P04 2 -128.75394

H2P04 2 -116.68027

HP04 2 -75.72752

2.,0.,1.,0.

l.,O.,O.,l.

0.,1.,4.,-3.

1.,0.,1.,-1.

3.,1.,4.,0.

2.,1.,4.,-1.

1.,1.,4.,-2.

55.506,0.000022,0.000001,1.E-5,1.E-7,0.000001

111.012,1.E-7,55.506,0.0

298.15,1.0,0.0,0.0

H F O P

HIDROLISIS DE ACID0 FOSFORICO

63

en el formato de entrada el primer renglón es NRUNS, que indica que sólo se correrá

un problema. En el segundo renglón, el primer número indica las especies participan -

tes ( 7 ) , el siguiente, el número de elementos (41, uno de ellos es la carga y l o s - o

t r o s corresponden a 0,H jr P. El siguiente número indica las especies presentes en

la fase gas (O), el otro número inidica las especies en fase líquida ( 7 ) . El número

1 que sigue representa el tipo de datos, es decir, las unidades en que se trabaja y

corresponden, en este caso, a KJ/mol. El cero indica que la estimación inicial la

dará el usuario y el -1, que los elementos del vector fórmula son no enteros.

El siguiente bloque de datos corresponde a los nombres de las especies repre--

sentadas por sus fórmulas y el número 2 que les sigue, indica que se encuentran en

la fase líquida. Por último, la cifra que sigue corresponde al valor del potencial,

calculado de la manera que se explicó anteriormente, a partir de los valores repor

tados de las constantes de equilibrio. El valor del potencial para las especies que

se consideran componentes es cero.

-

El bloque de datos siguiente es el que indica la abundancia de elementos pre--

sentes en cada especie; la primera columan corresponde a los Hidrógenos, la segun-

da a Fósforos, la tercera a Oxígenos y la última a las cargas de cada especie.

En el renglón que sigue se anotan las estimaciones iniciales según el usuario,

la siguiente línea corresponde a las abundancias de los elementos que deben cumplir

con el balance de materia. El renglón posterior especifica datos de presión y tempe

ratura del sistema, así como también número de moles de gas y líquido inertes res--

pectivamente. El penúltimo renglón da los nombres de los elementos, donde F es para

designar al Fósforo y P indica la designación de la carga y por último se asigna un

título al problema para efectos de identificación.

-

64

3.5.1 Salida de datos.

RESULTADOS DEL ARCHIVO DE ENTRADA DE DATOS

VCS CALCULATION METHOD


7 SPECIES 4ELEMENTS 3COMPONENTS O PHASEl SPECIES 7 PHASE2 SPECIES O SINGLE SPECIES PHASES

PRESSURE 1.000 ATM TEMPERATURE 298.150 K PHASE2 INERTS . O00

ELEMENTAL ABUNDANCES CORRECT FROM ESTIMATE

H 1.113120000000D+02 1.113120470000D+02

O 5.590600000000D+01 5.590616200000D+01 F 1.000000000000D-01 1.000380000000D-01

P 0.000000000000D+00 -8.700000000000D-05

USER ESTIMATE OF EQUILIBRIUM STAN.CHEM. POT. IN KJ/MOLE

SPECIES FORMULA VECTOR STAN.CHEM.POT. EQUILIBRIUM EST.

H F O P SI(1) H20 2.000 .O00 1.000 .O00 -.79540D+02 .55506D+02 H3P04 3.000 1.000 4.000 ,000 -.12875D+03 .78000D-01 H 1.000 .O00 .O00 1.000 .00000D+00 .22000D-01 OH 1.000 .O00 1.000 -1.000 .00000D+00 .10000D-04 PO4 .O00 1.000 4.000 -3.000 .00000D+00 .10000D-05 H2P04 2.000 1.000 4.000 -1.000 -.11668D+03 .22000D-01 HP04 1.000 1.000 4.000 -2.000 -.75720D+02 .37000D-04

ITERATION= 6 EVALUATIONS OF STOICHIOMETRY= 4

SPECIES EQUILIBRIUM MOLES MOLE FRACTION DG/RT REACTION

H20 5.5506000D+01 9.9776851D-01 H3P04 7.5861852D-02 1.3636826D-03 H 2.4138214D-02 4.3390534D-04 H2P04 2.4138081D-02 4.3390295D-04 -4.3244D-10 HP04 6.6676815D-08 1.1985736D-09 1.9131D-06 OH 4.8000306D-13 8.6284714D-15 9.2578D-07 PO4 1.4977485D-19 2.6923329D-21 2.8388D-06

G/RT= -1.7865802D+03 TOTAL PHASEl MOLES=,O.OOOOD+OO TOTAL PHASE2 MOLE=,lP .5563D+02

ELEMENTAL ABUNDANCES H 1.113120000000D+02 F 1~000000000000D-01 O 5.590600000000D+01 P -3.007389341026D-15

65

Los datos de salida muestran que la convergencia es alcanzada después de seis

iteraciones y cuatro cálculos de la matriz estequiométrica, así como también puede

observarse el resultado del cálculo del número de moles de las especies en equili-

brio, de sus fracciones mol correspondientes y finalmente el cambio / G / R T para la

ecuación estequiométrica en la cual un mol de la especie indicada , formada de las

primeras tres especies consideradas componentes (H20 ,H3P04, H) , aparece debajo de

DG/RT REACTION, los valores se obtienen después de recalcular la estequimetria.

También se obtiene el valor de G/RT mínimo para el sistema, el númer3 total

de moles en cada fase y la abundancia de elementos al final del cálculo, la cual

coincide con el valor inicial dado.

66

IV. APLICACION A SISTEMAS DE DISOLUCIONES

4.1 Cálculo de concentraciones en equilibrio.

El cálculo de las concetraciones en equilibrio de especies presentes en una diso

lución, se torna cada vez más dificil a medida que aumenta el número de éstas. El au-

mento puede deberse no solo a que se trabaja con mezclas, sino también a que las espe

ties pueden reaccionar entre ellas y las especies producidas también reaccionen entre

ellas mismas y las especies originales, si las concentraciones lo permiten.

-

-

En un sistema de disoluciones de distintas especies donde es importante determi-

nar las concentraciones de equilibrio, se hace necesario contar con un método O mode-

lo que facilite ésta tarea. Herramientas importantes son los programas de computadora

como el que se describió anteriormente, ya que manejan varias especies e inclusive va

rias fases al mismo tiempo. Otros métodos son los métodos gráficos, de entre éstos se

encuentran los más utilizados como el de diagramas de distribución (Hogfeldt,1979 :--

Vicente-Perez, 1979), los diagramas logarítmicos, los de zonas de predominio lineales

propuestos y aplicados por Charlot (1967). y los bidimensionales propuestos por Pour--

baix (1966) y Rossotti (1961).

-

Debido a que la presentación de dichos diagramas se hace de manera independiente

y excluyente y a que pcr tanto no son muy difundidos, los profesores de la Universi"

dad Autónoma Metropolitana- Iztapalapa, Dr. Ignacio González y Dr. Alberto Rojas del

departamento de Química , realizaran una serie de trabajos sobre &tos diagramas como

un auxiliar en la ensenanza de la química de la3 disoluciones para los alumnos de Qui

mica Analítica.

Basándose en una parte fundamental para la elaboración de estos diagramas, como

es el establecimiento correcto de los equilibrios representativos, se calcularon las

concentraciones de equilibrio que fueron alimentadas como aproximaciones iniciales al

cálculo en computadora de las concentraciones que hacían que el sistema tuviera la e-

nergía libre mínima.

El propósito de esto era comprobar que por medio del modelo de equilibrios repre

67

Equilibrios sucesivos.


DIAGRAMAS DE ZONAS DE PI?F,DOMINIO

H3P04 H2POi

I I -I"- pH HPO;'

2.2 7.2 12.4

H3P0; HPO;' Pi53

Concentración 10- M

H PO + H20 ====== H2POi + 3 4 + H30

K = 10 M' = 10 -1 -2.2 d

resolviendo para dl queda la siguiente concentración

1 = 0.2216

68

Kd 2 = oc' (0.022) = 10 -7.2

(1-W)

resolviendo para alfa dos

&, = 10 -2.77

IHPOi*I = IH O + I = 3.715 "10 -5 3

HP0; ' + H ~ O ======== poi3 + H30 +

i 3.715*10-~ - -

eq 3. 715*10-5( 1-OL) 0C3.715~10- 5 3. 715*10-5

69

de aquí se obtiene el valor de alfa tres

-3.98 w 3 = 10

IPOq31 = IH O + I = 10 -8.45

4 3

Así considerando la dilución de la solución tenemos para esta ap

siguientes resultados en el equilibrio

IH PO I = 0.078 3 4

IH PO-/ = 0.022 2 4

IH PO;'\= 3.715*10-5

I Poi3 I = 3.5481*10-'

IH30 I = 0.022

I OH 1 = 4.57*10

+

-13

IH20 I = 55

roximació In los

La consideración de la concentración del agua tan alta es arbitraria. Todas las

concentraciones están dadas en molaridad.

4.3 Comparación de resultados.

En la siguiente tabla se presentan los resultados obtenidos con el modelo de e-

quilibrios sucesivos y los obtenidos con el algoritmo.

Cabe reiterar que los cálculos fueron hechos considerando la concentración ini-

cial de 0.1 M y que se partió del ácido fosfórico .

70

T A B L A 1

ESPECIE

H20

H3P04

H

OH

PO4

H2P04

HP04

M E S

55.506

O. 078

o. 022

4.57*10-'~

3.5481*10-'

o. 022

3. 715*10-5

A. MINELI

55.506000

7. 5861852*10-2

2. 4138214*10-2

4. 8000360*10-13

1.4~7485*10-~'

2. 4138081*10-2

6. 6676815*10-8

Todas las concentraciones están expresadas en molaridad. MES : Modelo de Equilibrios Sucesivos. A. MINELI: Algoritmo de Minimización Energía Libre

En la tabla anterior se puede observar que los valores obtenidos con uno y otro

método son muy cercanos, sin embargo en la concentraciones de las especies HP04 y la

de PO4 hay uan variación considerable.

Cabe mencionar dos aspectos importantes; 1) en la hoja de salida de los resulta

dos del algoritmo se observa que la abundancia de elementos checa la inicio y al fi-

-

nal del cálculo, la carga, que debe ser neutra, aparece como P= -3.0078...*10 -15 la

cual se considera como buena ya que el programa tiene como tolerancia 10 . 2) Las

cantidades dadas en la aproximación inicial para OH-,POi3 y HPOi' fueron dadas consi

derando cantidades apropiadas al formato de entrada del algoritmo y no los datos cal-

-14

-

culados con el modelo de equilibrios sucesivos. lo que pudo haber originado la varia

ción en los resultados. Puesto que la aproximación inicial dada al programa es preci

samente el valor calculado de concentraciones de equilibrio obtenidas con el modelo

de equilibrios sucesivos, puede decirse que tal modelo es una herramienta tan buena

como el programa de computadora, al menos para las condiciones en que se realizó el

cálculo.

-

71

V. CONCLUSIONES

Tanto en uno, como en otro método de los manejados aquí, es importante el cono -

cimiento de la naturaleza de las sustancias presentes. En el algoritmo es más difi-

cil analizar los resultados, por que el problema se convierte en un problema numéri -

co, en cambio por medio del modelo de diagramas de zonas de predominio y de los e--

quilibrios representativos es más fácil visualizar los equilibrios que influyen de

manera determinante para el cálculo de las concentraciones de equilibrio.

De cualquier manera los dos métodos son importantes y Gtiles y su validez radi -

ca, no solo en el método mismo, sino en la destreza con que uno los maneje e inter-

prete.

VI. BIBLIOGRAFIA

Kreysig, Erwing, (1988), Advanced Engineering Mathematics.

Smith,J.M., y H.C, Van Ness, (1975), Introducción a la Termodinámica en Ingeniería Química, 2a. Ed. en español, McGraw-Hill,

Guerasimiv, Y.A, (1977), Curso de Química-Física, 2a. Ed., Tomo I, MIR, Moscú.

Smith,W.R y R.V. Missen, (1987), Análisis del equilibrio en reacciones químicas, la Ed., en espaflol, Limusa, México.

Smith.W.R y R.V. Missen, (1980), Chemical Reaction Equilibrium Analysis: Theory and Algorithms, Ed. John Wiley & Sons. Inc.

Laitenen,H.A, (1960), Chemical Analysis, McGraw-Hil1,Nueva York.

Gonzále2,M.I y Rojas A,(1988), Contactos, Vol. 111, No 2

González,M.I y Rojas A, (1988),Contactos, Vol. 111, No 3

González, M.1 y Rojas A, (1989), Contactos, Vol. IV, No 1

Huheey, James E, (1981), Química Inorgánica, 2a. Ed., en español, Harla, México.

72

APENDICE A : Errores causados por una mala entrada de datos

A continuación se presenta una salida de datos que presenta varios errores de-

bido a que la entrada de los datos no respetó el formato.

VCS CALCULATION METHOD

298.15,1,0.0,0.0

7 ESPECIES 3 ELEMENTS 3COMPONENTS O PHASEl SPECIES 7 PHASE2 SPECIES O SINGLE SPECIES PHASES

PRESSURE .O00 ATM TEMPERATURE .005K PHASE2 INERTS . O00

ELEMENTAL ABUNDANCES CORRECT FROM ESTIMATE

10 1.000000000000D-08 1.000000000000D-09 O. 0.000000000000D+00 4.000000000000D-10 O9 1.000000000000D-08 1.800000000000D-09

MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM STAN. CHEM. POT. IN KJ/MOLE

SPECIES FORMULA VECTOR STAN.CHEM.POT. EQUILIBRIUM EST.

10 o. o9 SI(1) H3P04 3.000 1.000 4.000 2 -.12875D+03 .10000D-09 H2P04 2.000 1.000 4.000 2 -.11668D+03 .10000D-ll OH 1.000 .O00 1.000 2 .00000D+00 .25255D-08

PO4 .O00 1.000 4.000 2 .00000D+00 .10000D-09 HP04 1.000 1.000 4.000 2 -.75720D+02 a10000D-09 H20 2.000 .O00 1.000 2 .00000D+00 .10000D-09 ITERATION = 18 EVALUATIONS OF STOICHIOMETRY= 10

H 1.000 .O00 ,000 2 .00000D+00 .10000D-09

SPECIES EQUILIBRIUM MOLES MOLE FRACTION DG/RT REACTION

OH H2P04 HP04

9.9995198D-09 9.999980D-01 1.9600589D-14 1.9601492D-06 1.1556387D-24 1.1556920D-16

LESS THAN 1.E-32 MOLES

H PO4 H3P04 H20

G/RT= -5.5013850D-08 TOTAL PHASEl MOLES=,O.OOOOD+OO TOTAL PHASE2 MOLE=,lP .1000D-07

ELEMENTAL ABUNDANCES 10 9.999377702735D-09 O. 2.765765622337D-14 O9 9.999433018047D-09

73

El archivo de entrada correspondiente a l o s resultados anteriores es el siguien -

te :

1 7,3,0,7,1,-1,-1,0 OH 2 0.00000 PO4 2 0.00000 H 2 0.00000 H3P04 2 -128.75394 H2P04 2 -116.68027 HP04 2 -75.72752 H20 2 0.00000 1.,0.,1.,0.,0.,0.,0. 0.,1.,4.,0.,0.,0.,0. 1.,0.,0.,0.,0.,0.,0. 3.,1.,4.,0.,0.,0.,0. 2.,1.,4.,0.,0.,0.,0. 1.,1.,4.,0.,0.,0.,0. 2.,0.,1.,0.,0.,0.,0. 0.00000001,0.00000001,0.022,0.078,0.~22,0.00003715,50 50 100.0999,0.1,50.0999 298.15,1,0.0,0.0 H P O HIDROLISIS DE ACID0 FOSFORICO

En esta entrada no se está tomando en cuenta la carga, e-l O que aparece al íiltimo de

la segunda línea de datos está demás. El agua se esta considerando como un componen-

te, el vectro fórmula tiene cuatro columnas demás considerando que solo se están to-

mando en cuenta tres elementos. También.sobra el renglón de la estimación inicial ya

que se esta indicando con -1 (el primer -1 localizado en el segundo renglón del ar-

chivo) que la estimación inicial la da la máquina.

El no respetar el formato de entrada confunde al compilador y hace que realice,

como se ve en la salida de resultados, cálculos erróneos. El error más evidente co--

rresponde al cálculo de abundancia de elementos que no se conserva al final del cál -

culo. Otros errores que ocasiona son; el que no aparezca el título arbitrario,los da

tos de presión y temperatura y los nombres de los elementos. Cabe hacer notar que el

-

orden de las especies es diferente al del ejemplo dado en la sección 3.5, lo cual, a

su vez, hace evidente la importancia que tiene el darle cierto orden a las especies,

el asignarle el potencial adecuado a cada especie y el análisis de los valores G / RT

que se obtienen.

74

UNIVERSIDAD AUTONOMA METROPOLITANA ...148.206.53.84/tesiuami/UAM1134.pdfbrio de un sistema es un...

Documents

Transcript of UNIVERSIDAD AUTONOMA METROPOLITANA ...148.206.53.84/tesiuami/UAM1134.pdfbrio de un sistema es un...