*****************************************************************
*** Chapter 8
*** Section: Phase and Chemical Equilibrium Problems - Equations of State
*** Test Problem 5
***
*** Peng-Robinson equation
*** Tangent Plane distance minimization
*** Carbon Dioxide -- Methane
*****************************************************************

SETS
      i   components        /1*2/


*****************************************************************
* P = pressure (bar)
* T = temperature (K)
* feedmf(i) = mole fraction of component i in candidate phase
* feedz = compressibility of candidate phase
* feedfc(i) = fugacity coefficient of component i in candidate phase
* b(i) = Peng-Robinson pure-component parameter
* a(i,j) = Peng-Robinson mixture parameter

PARAMETERS
P, T, feedmf(i), feedz, feedfc(i), b(i), a(i,j);

P = 60.8;
T = 220.0;

feedmf('1') = 0.20;
feedmf('2') = 0.80;

feedz = 0.445886;

feedfc('1') = -0.965195;
feedfc('2') = -0.323254;

b('1') = 0.0885973;
b('2') = 0.0890893;

a('1','1') = 0.884831; a('1','2') = 0.555442;
a('2','1') = 0.555442; a('2','2') = 0.427888;

*******************************************************************
* dist = tangent plane distance
* x(i) = mole fractio of component i in incipient phase
* z = compressibility of incipient phase
* amix = mixture A parameter (function of composition)
* bmix = mixture B parameter (function of composition)

VARIABLES dist, x(i), z, amix, bmix;


*******************************************************************
* obj = objective (tangent plane distance)
* eos = equation of state constraint
* defa = definition of Amix
* defb = definition of Bmix
* molesum = mole fractions sum to 1

EQUATIONS
      obj
      eos
      defa
      defb
      molesum;

obj.. dist =e= SUM(i, x(i)*LOG(x(i)))
              +z-1.0
	      -LOG(z-bmix)
	      -amix*LOG((z+(1+SQRT(2))*bmix)/(z+(1-SQRT(2))*bmix))/(2*SQRT(2)*bmix)
	      -SUM(i, x(i)*(LOG(feedmf(i))+feedfc(i)));

eos.. POWER(z,3) - (1-bmix)*POWER(z,2) + (amix - 3.0*POWER(bmix,2) - 2.0*bmix)*z
     - amix*bmix + POWER(bmix,3) + POWER(bmix,2) =e= 0;
defa.. amix - SUM(i,SUM(j, a(i,j)*x(i)*x(j))) =e= 0;
defb.. bmix - SUM(i, b(i)*x(i)) =e= 0;
molesum.. SUM(i, x(i)) =e= 1.0;

MODEL tpd /all/;