kill(all)$

Bp calculation of mixture (with relativ volatilities,ideal),version 1.0
peter.vlasschaert@gmail.com, 09/08/2020

youtube :Upstream LNG Technology by Prof. Pavitra Sandilya
============================================
=>  Lecture 61: Tutorial on equillibrium in multicomponent systems

ref :youtube
========
https://www.youtube.com/watch?v=....
NPrAjRaGhNQ&list=PLbRMhDVUMngd8Lw6T6edxXqeXLQ2DR7db&index=65

load packages, option.

load(lrats)$
ratprint:false$

conversion : bar to mmHg pressure ( Antoine equation ,'vapor pressure')
==========================

MKS :system
------------------

1.5bar = 150000 Pa  (N/m^2)  → mm Hg
info: ρ Hg = 13600 kg/m^3, g =9.81 m/s^2, h = ? mHg  ( *1000 ) = ? mmHg
         

pa:[ρ =13600,g = 9.81];

p1: p = ρ*h*g;

p1a: solve(p1,h)[1],numer;

FORMULA : 1
==========

p1b: lhs(lratsubst(pa,p1a)) = 1000*rhs(lratsubst(pa,p1a));

p1c:lratsubst([p=150000],p1b),numer  /*mm Hg */;

DEFINITION : Antione equation ( T = [°C] , psat = [mmHg])
==========================================
→  log(psat) = A- B/(T+C), math : log = ln

FORMULA : 2
==========

p2a: psat(i,m1,m2,m3,TT):=exp((m1-m2/(TT+m3)));

FORMULA : 3
==========

p2b: Tsat(i,m1,m2,m3,pp):= m2/(m1-log(pp))-m3;

→ relativ volatility : symbol 'r'

FORMULA : 4a
===========

p3a : αk[ki,kj]:=ki/kj;

p3a1:αkr[ki,kr]:=ki/kr;

FORMULA : 4b
===========

p3b : αp[pi,pj]:=pi/pj;

p3b1:αpr[pi,pr]:=pi/pr;

given : problem :determine buble point temperature ( = ideal mixture)

----------------------------------------
data : three components ( =nc)
----------------------------------------

=======================================
component  ,bp( °C), A,B,C, MW ,mole faction.
=======================================
1e) n-pentane ;36, 15.8365,2477.07,233.21,72,0.2
2e) n-hexane ,69,15.9155,2738.42,226.20, 86,0.4
3e) n-heptane ,98.5,15.8770.2911.32,226.65, 100,0.4
=======================================
given : total pressure of the mixture :pp=1.5bar
=======================================
nc = number of components

1e) step : build a data matrix
=====================

MC1:[36, 15.8365,2477.07,233.21,72,0.2];
MC2:[69,15.9155,2738.42,226.20, 86,0.4];
MC3:[98.5,15.877,2911.32,226.65, 100,0.4];

MC:matrix(MC1,MC2,MC3);

2e) step :how to find ref component (=r) , from the data.
====================================

nc :3 /*number of components*/;

FORMULA : 5
===========
need column : 5 & 6 , 'MC'
MW[avg] : molecular average weight

p4:MW[avg]= sum(MC[i,5]*MC[i,6],i,1,nc);

?  Which row I can find ref component

a) create list  MW for all components in the mixture
------------------------------------------------------------------

MWL:makelist(MC[i,5],i,1,nc);

b) absolute value of difference between MWL & MW[avg]
--------------------------------------------------------------------------

MWU:abs(MWL-rhs(p4));

c) find minimum of MWU 
----------------------------------

lmin(MWU);

rem : lmin,lmax also be found in the list you want, numbers.

fref[mr]:=block(
for i:1 thru length(mr)do(
  if mr[i]= lmin(mr) then
       ref:i
    ));

d) ref = r ,see 'above'
---------------------------

fref[MWU];

ref;

3e) step : how to find temperature of the mixture
===================================
use : FORMULA 3 (p2b) , @ 1.5 bar,see begin 'above'

see :conversion : bar to mmHg pressure

need use :column 2,3,4 for data Antione equation

lTsat1: makelist(Tsat(i,MC[i,2],MC[i,3],MC[i,4],pp),i,1,nc);

lTsat:lratsubst([pp=rhs(p1c)],lTsat1),numer
;

result :3e) step :temperature to use : ref (temperature)
----------------------------------------------------------------------
update temperature::'ref'

Tref:lTsat[ref];

4e) step  use  FORMULA 2 (p2a) , @ 1.5 bar,see begin 'above'
=============================================

lpsat:makelist(psat(i,MC[i,2],MC[i,3],MC[i,4],Tref),i,1,nc);

5e) step use relativ volatilities of components 
=================================
use :FORMULA : 4b  & use 4e) step

lrvol:makelist(αpr[lpsat[i],rhs(p1c)],i,1,nc);

6e) step use  kr (formula), r=ref
=======================
use :column 6 of MC use 5e)

kr=1/sum(lrvol[i]*MC[i,6],i,1,nc);

7e) step calculation new temperature (psat(ref))
==================================
psat(ref) = pp*kr

prefsat:rhs(p1c)*rhs(%);

8e) step calculation update temperature with respect to (psat(ref))
============================================
use :'p2b',use column :2,3,4

trefsat:Tsat(ref,MC[ref,2],MC[ref,3],MC[ref,4],prefsat);

9e) step calculation update psat with the temperature 8e)
============================================
use :p2a

luppsat:makelist(psat(i,MC[i,2],MC[i,3],MC[i,4],trefsat),i,1,nc);

10e) step calculation : equilibrium values , K[j] , j =1...nc
============================================

lkval: luppsat/rhs(p1c);

11e) step  check : 'bp equation'
=======================
keqsum:sum (K[j]*x[j],j,1,nc) ?  = 1 , or cvg to one

keqsum:sum(lkval[i]*MC[i,6],i,1,nc);

result : bubble point temperature , trefsat=77.37099058331643 °C

update : cvg for temperature ( keqsum = 1)
===============================

when not satisfied go back " keqsum ≠ 1 " : make update calculation :
start  from step  5e)