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

SETS
      i   components        /1*3/


*****************************************************************
* 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 = 76.0;
T = 270.0;

feedmf('1') = 0.15;
feedmf('2') = 0.30;
feedmf('3') = 0.55;

feedz = 0.448135;

feedfc('1') = 0.468354;
feedfc('2') = -0.067012;
feedfc('3') = -1.0288;

b('1') = 0.0815247;
b('2') = 0.0907391;
b('3') = 0.13705;

a('1','1') = 0.142724; a('1','2') = 0.206577; a('1','3') = 0.342119;
a('2','1') = 0.206577; a('2','2') = 0.323084; a('2','3') = 0.547748;
a('3','1') = 0.342119; a('3','2') = 0.547748; a('3','3') = 0.968906;


*******************************************************************
* 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/;