Description
The chemical plant performed the hydro-dealkylation of toluene into benzene and methane. The flowsheet model was used to make decisions on choosing between alternative process units at various stages of the process. The resulting model is nonlinear and mixed integer.
Large Model of Type : MINLP
Category : GAMS Model library
Main file : hda.gms
$title Synthesis: Hydrodealkylation of Toluene (HDA,SEQ=123)
$onText
The chemical plant performed the hydro-dealkylation of toluene
into benzene and methane. The flowsheet model was used to make
decisions on choosing between alternative process units at
various stages of the process. The resulting model is nonlinear
and mixed integer.
Kocis, G R, and Grossmann, I E, Computational Experience with DICOPT
Solving Minlp Problems in Process Synthesis. Computers and Chemical
Engineering 13, 3 (1989), 307-315.
-----------------
- process units -
-----------------
absorber abs
compressors comp
distillation dist
flash flsh
furnace furn
utility-stream exchangers hec & heh
stream-stream exchanger exch
membrane separator memb
mixers (single inlet stream) mxr1
mixers mxr
pump pump
reactor rct
splitters (single outlet stream) spl1
splitters spl
expansion valve valve
-----------------------
- process unit models -
-----------------------
absorber kremser equation
(douglas, page 431)
neglect heat effects
assume: pure solvent
33 % tray efficiency
l / m g = 1.4 (benzene)
fixed recoveries for
h2 , ch4 , dip
compressors isentropic compression of ideal gas
(mccabe and smith, page 201)
distillation rmin from underwood equation
(douglas, page 441)
r = 1.2 * rmin (heuristic)
nmin from fenske equation
(douglas, page 441)
nt = 2 * nmin (gillilands approx.)
(douglas, page 443)
n = nt / efficiency
expansion valve isentropic expansion of ideal gas
(mccabe and smith, page 200)
flash ideal flash (raoult's law)
(mccabe and smith, page 484)
furnace 50 % efficiency
heat exchangers
membrane separator shortcut method (driving force
approximated as arithmetic mean)
(grossmann, design notes)
mixers (single inlet stream) linear model for heat & mass balances
(kocis and grossmann, 1988)
mixers nonlinear heat balance
(kocis and grossmann, 1988)
pump pout > pin
reactor plug flow reactor
1. isothermal
2. adiabatic - use average temp
splitters (single outlet stream) linear model for heat & mass balances
(kocis and grossmann, 1988)
splitters split fraction model
(kocis and grossmann, 1988)
Keywords: mixed integer nonlinear programming, chemical engineering, toluene
hydrodealkylation, process synthesis
$offText
Scalar
alpha 'compressor coefficient' / 0.3665 /
compeff 'compressor efficiency' / 0.750 /
gam 'ratio of cp to cv' / 1.300 /
abseff 'absorber tray efficiency' / 0.333 /
disteff 'column tray efficiency' / 0.500 /
uflow 'upper bound - flow logicals' / 50 /
upress 'upper bound - pressure logicals' / 4.0 /
utemp 'upper bound - temperature logicals' / 7.0 /
* want costs in units of $1000. per year
costelec 'electricity cost ($340. per kw-yr)' / 0.340 /
costqc 'cooling cost ($700 per 1e9 kj)' / 0.700 /
costqh 'heating cost ($8000 per 1e9 kj)' / 8.000 /
costfuel 'fuel cost furnace ($4 per 1e6 btu)' / 4.0 /;
$onText
user must fill in the data for the following sets to define
the number of process streams, components, and process units
$offText
Set
str 'process streams' / 1*35, 37*73 /
compon 'chemical components' / h2 'hydrogen', ch4 'methane', ben 'benzene', tol 'toluene', dip 'diphenyl' /
abs 'absorber' / 1 /
comp 'compressors' / 1*4 /
dist 'distillation' / 1*3 /
flsh 'flash separators' / 1*3 /
furn 'furnace' / 1 /
hec 'cooling units' / 1*2 /
heh 'heating units' / 1*4 /
exch 'exchanger' / 1 /
memb 'membrane separators' / 1*2 /
mxr1 'mixers (1 input stream)' / 1*5 /
mxr 'mixers' / 1*5 /
pump 'pump' / 1*2 /
rct 'reactors' / 1*2 /
spl1 'splitters (1 input stream)' / 1*6 /
spl 'splitters' / 1*3 /
valve 'expansion valves' / 1*6 /
$onText
user must fill in the data for the following sets to
define the input-output structure of the superstructure
by defining which streams are input and/or output streams
for the process units
$offText
* =======================
* = =
* = a b s o r b e r =
* = =
* =======================
ilabs(abs,str) 'abs-stream (inlet liquid) matches' / 1.67 /
olabs(abs,str) 'abs-stream (outlet liquid) matches' / 1.68 /
ivabs(abs,str) 'abs-stream (inlet vapor) matches' / 1.63 /
ovabs(abs,str) 'abs-stream (outlet vapor) matches' / 1.64 /
asolv(abs,compon) 'abs-solvent component matches' / 1.tol /
anorm(abs,compon) 'abs-comp matches (normal model)' / 1.ben /
asimp(abs,compon) 'abs-heavy component matches' / 1.(h2,ch4,dip) /
* ===========================
* = =
* = c o m p r e s s o r =
* = =
* ===========================
icomp(comp,str) 'compressor-stream (inlet) matches' / 1.5, 2.59, 3.64, 4.56 /
ocomp(comp,str) 'compressor-stream (outlet) matches' / 1.6, 2.60, 3.65, 4.57 /
* ===============================
* = =
* = d i s t i l l a t i o n =
* = =
* ===============================
idist(dist,str) 'dist-stream (inlet) matches' / 1.25, 2.30, 3.33 /
vdist(dist,str) 'dist-stream (vapor) matches' / 1.26, 2.31, 3.34 /
ldist(dist,str) 'dist-stream (liquid) matches' / 1.27, 2.32, 3.35 /
dl(dist,compon) 'dist-light components matches' / 1.h2 , 2.ch4, 3.ben /
dlkey(dist,compon) 'dist-light key component matches' / 1.ch4 , 2.ben, 3.tol /
dhkey(dist,compon) 'dist-heavy key component matches' / 1.ben , 2.tol, 3.dip /
dh(dist,compon) 'dist-heavy components matches' / 1.(tol,dip), 2.dip /
dkey(dist,compon) 'dist-key component matches';
dkey(dist,compon) = dlkey(dist,compon) + dhkey(dist,compon);
* =================
* = =
* = f l a s h =
* = =
* =================
Set
iflsh(flsh,str) 'flsh-stream (inlet) matches' / 1.17 , 2.46 , 3.39 /
vflsh(flsh,str) 'flsh-stream (vapor) matches' / 1.18 , 2.47 , 3.40 /
lflsh(flsh,str) 'flsh-stream (liquid) matches' / 1.19 , 2.48 , 3.41 /
fkey(flsh,compon) 'flash-key component matches' / 1.ch4, 2.ch4, 3.tol /
* =======================
* = =
* = f u r n a c e =
* = =
* =======================
ifurn(furn,str) 'furn-stream (inlet) matches' / 1.70 /
ofurn(furn,str) 'furn-stream (outlet) matches' / 1.9 /
* ===================
* = =
* = c o o l e r =
* = =
* ===================
ihec(hec,str) 'hec-stream (inlet) matches' / 1.71, 2.45 /
ohec(hec,str) 'hec-stream (outlet) matches' / 1.17, 2.46 /
* ===================
* = =
* = h e a t e r =
* = =
* ===================
iheh(heh,str) 'heh-stream (inlet) matches' / 1.24, 2.23, 3.37, 4.61 /
oheh(heh,str) 'heh-stream (outlet) matches' / 1.25, 2.44, 3.38, 4.73 /
* =========================
* = =
* = e x c h a n g e r =
* = =
* =========================
icexch(exch,str) 'exch-cold stream (inlet) matches' / 1.8 /
ocexch(exch,str) 'exch-cold stream (outlet) matches' / 1.70 /
ihexch(exch,str) 'exch-hot stream (inlet) matches' / 1.16 /
ohexch(exch,str) 'exch-hot stream (outlet) matches' / 1.71 /
* =======================
* = =
* = m e m b r a n e =
* = =
* =======================
imemb(memb,str) 'memb-stream (inlet) matches' / 1.3, 2.54 /
nmemb(memb,str) 'memb-stream (non-permeate) matches' / 1.4, 2.55 /
pmemb(memb,str) 'memb-stream (permeate) matches' / 1.5, 2.56 /
mnorm(memb,compon) 'normal components' / 1.(h2,ch4) , 2.(h2,ch4) /
msimp(memb,compon) 'simplified flux components' / 1.(ben,tol,dip), 2.(ben,tol,dip) /
* ==================
* = =
* = m i x e r =
* = (single input) =
* = =
* ==================
imxr1(mxr1,str) 'mixer-stream (inlet) matches' / 1.(2,6), 2.(11,13), 3.(27,48), 4.(34,40), 5.(49,50) /
omxr1(mxr1,str) 'mixer-stream (outlet) matches' / 1.7 , 2.14 , 3.30 , 4.42 , 5.51 /
mxr1spl1(mxr1,str,str) '1-mxr-inlet 1-spl-outlet matches' / 1.(2.2,6.3) , 2.(11.10,13.12)
3.(27.24,48.23), 4.(34.33,40.37)
5.(49.23,50.24) /
* =================
* = =
* = m i x e r =
* = =
* =================
imxr(mxr,str) '-stream (inlet) matches' / 1.(7,43,66,72), 2.(15,20)
3.(21,69) , 4.(51,62)
5.(57,60,65) /
omxr(mxr,str) 'mixer-stream (outlet) matches' / 1.8, 2.16, 3.22, 4.63, 5.72 /
* ===============
* = =
* = p u m p =
* = =
* ===============
ipump(pump,str) 'pump-stream (inlet) matches' / 1.42, 2.68 /
opump(pump,str) 'pump-stream (outlet) matches'/ 1.43, 2.69 /
* =====================
* = =
* = r e a c t o r =
* = =
* =====================
irct(rct,str) 'reactor-stream (inlet) matches' / 1.10 , 2.12 /
orct(rct,str) 'reactor-stream (outlet) matches' / 1.11 , 2.13 /
rkey(rct,compon) 'reactor-key component matches' / 1.tol, 2.tol /
* =======================
* = =
* = s p l i t t e r =
* = (single output) =
* = =
* =======================
ispl1(spl1,str) 'splitter-stream (inlet) matches' / 1.1 , 2.9 , 3.22
4.32 , 5.52 , 6.58 /
ospl1(spl1,str) 'splitter-stream (outlet) matches' / 1.(2,3) , 2.(10,12), 3.(23,24)
4.(33,37), 5.(53,54), 6.(59,61) /
* =======================
* = =
* = s p l i t t e r =
* = =
* =======================
ispl(spl,str) 'splitter-stream (inlet) matches' / 1.19, 2.18, 3.26 /
ospl(spl,str) 'splitter-stream (outlet) matches' / 1.(20,21), 2.(52,58), 3.(28,29) /
* =========================
* = =
* = e x p v a l v e =
* = =
* =========================
ival(valve,str) 'exp.valve-stream (inlet) matches' / 1.44, 2.38, 3.14
4.47, 5.29, 6.73 /
oval(valve,str) 'exp.valve-stream (outlet) matches' / 1.45, 2.39, 3.15
4.49, 5.50, 6.62 /;
Alias (str,str2), (compon,compon2);
$onText
heat capacities for pure components ( kj per kg-mol k ) -
vapor liquid
h2 30
ch4 40
ben & tol 220 230
dip 440 460
(h2, ch4, and tol vapor values from douglas
tol liquid value from mccabe and smith
ben liquid and vapor values assumed the same as tol
dip liquid and vapor values assumed 2 times ben
i.e. molecular weight of dip = 2 * ben )
approximate bulk stream cp by guessing a stream composition
and calculating a weighted average.
heat of vaporization of toluene taken from reid, prausnitz & sherwood
30890 kj per kg-mole
(same value assumed for ben if needed)
$offText
Parameter
cp(str) 'heat capacities ( kj per kgmole-k)'
heatvap(compon) 'heat of vaporization (kj per kg-mol)'
cppure(compon) 'pure component heat capacities' / h2 30, ch4 40, ben 225, tol 225, dip 450 /;
heatvap('tol') = 30890.;
Table gcomp(str,compon) 'guess composition values'
h2 ch4 ben tol dip
7 0.95 0.05
8 0.50 0.40 0.10
9 0.50 0.40 0.10
10 0.50 0.40 0.10
11 0.45 0.45 0.05 0.05
12 0.50 0.40 0.10
13 0.45 0.45 0.05 0.05
14 0.45 0.45 0.05 0.05
15 0.45 0.45 0.05 0.05
16 0.40 0.40 0.10 0.10
17 0.40 0.40 0.10 0.10
20 0.03 0.07 0.55 0.35
21 0.03 0.07 0.55 0.35
22 0.03 0.07 0.55 0.35
24 0.03 0.07 0.55 0.35
25 0.03 0.07 0.55 0.35
37 1.00
38 1.00
43 0.05 0.95
44 0.03 0.07 0.55 0.35
45 0.03 0.07 0.55 0.35
46 0.03 0.07 0.55 0.35
51 0.30 0.70
57 0.80 0.20
60 0.50 0.50
62 0.50 0.50
* Line for stream 63 estimated by Adrud based on Mixer data
63 0.47 0.40 0.01 0.12
65 0.50 0.50
66 1.00
69 1.00
70 0.50 0.40 0.10
71 0.40 0.40 0.10 0.10
72 0.50 0.50 ;
cp(str) = sum(compon, cppure(compon)*gcomp(str,compon));
display cp;
Parameter
anta(compon) 'antoine coefficient'
/ h2 13.6333, ch4 15.2243, ben 15.9008, tol 16.0137, dip 16.6832 /
antb(compon) 'antoine coefficient'
/ h2 164.9, ch4 897.84, ben 2788.51, tol 3096.52, dip 4602.23 /
antc(compon) 'antoine coefficient'
/ h2 3.19, ch4 -7.16, ben -52.36, tol -53.67, dip -70.42 /
perm(compon) 'permeability'
/ h2 55.0e-06, ch4 2.3e-06 /
cbeta(compon) 'constant values (exp(beta)) in absorber'
/ h2 1.0003, ch4 1.0008, dip 1.0e+04 /
aabs(compon) 'absorption factors'
/ ben 1.4, tol 4.0 /;
$onText
aabs(ben) fixed at 1.4 (hueristic) and
aabs(tol) calculated as: aabs(ben)*vp(tol)/vp(ben) at
t = 3.0 (temperature of absorber) ---> 4.0
units: cc (stp) / cm**2 - sec - cm hg
vol / area - time - pressure
1 g-mol = 22.4 l = 22.4 * 1000 cc
1000 g-mol = 1 kg-mole
1 m**2 = ( 100 cm )**2 = 1.0e+04 cm**2
76 cm hg = 1 atm = 0.101325 mpa
---> 750.062 cm hg = 1 mpa
60 sec = 1 min
convert units to: kg-mol / m**2 - min - mpa
1 psia = 6.894757e+03 pa = 6.894757e-03 mpa
$offText
perm(compon) = perm(compon)*(1./22400.)*1.0e+04*750.062*60./1000.;
Scalar eps1 'small number to avoid div. by zero' / 1.e-04 /;
Parameter
heatrxn(rct) 'heat of reaction (kj per kg-mol)'
/ 1 50100., 2 50100. /
f1comp(compon) 'feedstock compositions (h2 feed)'
/ h2 0.95, ch4 0.05, ben 0.00, tol 0.00, dip 0.00 /
f66comp(compon) 'feedstock compositions (tol feed)'
/ h2 0.00, ch4 0.00, ben 0.00, tol 1.00, dip 0.00 /
f67comp(compon) 'feedstock compositions (tol feed)'
/ h2 0.00, ch4 0.00, ben 0.00, tol 1.00, dip 0.00 /;
$onText
reaction: tol + h2 ---> ben + ch4 (desired)
2 ben ---> dip + h2 (undesired)
$offText
* =======================
* = =
* = a b s o r b e r =
* = =
* =======================
Positive Variable nabs(abs) 'number of absorber trays';
Variable gamma(abs,compon), beta(abs,compon);
* ===========================
* = =
* = c o m p r e s s o r =
* = =
* ===========================
Positive Variable
elec(comp) 'electricity requirement (kw)'
presrat(comp) 'ratio of outlet to inlet pressure'
* ===============================
* = =
* = d i s t i l l a t i o n =
* = =
* ===============================
nmin(dist) 'minimum number of trays in column'
ndist(dist) 'number of trays in column'
rmin(dist) 'minimum reflux ratio'
reflux(dist) 'reflux ratio'
distp(dist) 'column pressure'
distt(dist) 'column temperature'
avevlt(dist) 'average volatility'
* =================
* = =
* = f l a s h =
* = =
* =================
flsht(flsh) 'flash temperature (100 k)'
flshp(flsh) 'flash pressure (mega-pascal)'
eflsh(flsh,compon) 'vapor phase recovery in flash'
* =======================
* = =
* = f u r n a c e =
* = =
* =======================
qfuel(furn) 'heating requied (1.e+12 kj per yr)'
* ===================
* = =
* = c o o l e r =
* = =
* ===================
qc(hec) 'utility requirement (1.e+12 kj per yr)'
* ===================
* = =
* = h e a t e r =
* = =
* ===================
qh(heh) 'utility requirement (1.e+12 kj per yr)'
* =========================
* = =
* = e x c h a n g e r =
* = =
* =========================
qexch(exch) 'heat exchanged (1.e+12 kj per yr)'
* =======================
* = =
* = m e m b r a n e =
* = =
* =======================
a(memb) 'surface area for mass transfer (m**2)'
* =================
* = =
* = m i x e r =
* = (1 input) =
* = =
* =================
mxr1t(mxr1) 'mixer temperature (100 k)'
mxr1p(mxr1) 'mixer pressure (m-pa)'
* =================
* = =
* = m i x e r =
* = =
* =================
mxrt(mxr) 'mixer temperature (100 k)'
mxrp(mxr) 'mixer pressure (m-pa)'
* ===============
* = =
* = p u m p =
* = =
* ===============
* =====================
* = =
* = r e a c t o r =
* = =
* =====================
rctt(rct) 'reactor temperature (100 k)'
rctp(rct) 'reactor pressure (m-pa)'
rctvol(rct) 'reactor volume (cubic meter)'
krct(rct) 'rate constant'
conv(rct,compon) 'conversion of key component'
sel(rct) 'selectivity to Benzene'
consum(rct,compon) 'consumption rate of key'
q(rct) 'heat removed (1.e+9 kj per yr)'
* =======================
* = =
* = s p l i t t e r =
* = (1 output) =
* = =
* =======================
spl1t(spl1) 'splitter temperature (100 k)'
spl1p(spl1) 'splitter pressure (m-pa)'
* =======================
* = =
* = s p l i t t e r =
* = =
* =======================
splt(spl) 'splitter temperature (100 k)'
splp(spl) 'splitter pressure (m-pa)'
* =====================
* = =
* = s t r e a m s =
* = =
* =====================
vp(str,compon) 'vapor pressure (mega-pascal)'
e(str) 'split fraction'
f(str) 'stream flowrates (kg-mole per min)'
fc(str,compon) 'component flowrates (kg-mole per min)'
p(str) 'stream pressure (mega_pascal)'
t(str) 'stream temperature (100 k)';
Binary Variable y(str) 'binary variable';
Variable
const 'constant term in obj fcn'
profit 'overall profit ($1000 per year)';
* =======================
* = =
* = a b s o r b e r =
* = =
* =======================
Equation
absfact(abs,compon) 'absorbption factor equation'
gameqn(abs,compon) 'definition of gamma'
betaeqn(abs,compon) 'definition of beta'
abssvrec(abs,compon) 'recovery of solvent'
absrec(abs,compon) 'recovery of non-solvent'
abssimp(abs,compon) 'recovery of simplified components'
abscmb(abs,compon) 'overall component mass balance'
abslogc 'logical constraint'
abspl(abs) 'pressure relation for liquid'
abstl(abs) 'temperature relation for liquid'
abspv(abs) 'pressure relation for vapor'
absp(abs) 'pressure relation for inlets'
abst(abs) 'temperature relation at top';
absfact(abs,compon)$anorm(abs,compon)..
sum(str$ilabs(abs,str), f(str)*p(str)) =e=
sum(str$ivabs(abs,str), f(str))*aabs(compon)*
sum(str$ilabs(abs,str), vp.l(str,compon));
gameqn(abs,compon)$asolv(abs,compon)..
gamma(abs,compon) =e= log((1 - aabs(compon)**(nabs(abs)*abseff + eps1))/(1 - aabs(compon))); ;
betaeqn(abs,compon)$(not asimp(abs,compon))..
beta(abs,compon) =e= log((1 - aabs(compon)**((nabs(abs)*abseff) + 1))/(1 - aabs(compon))); ;
abssvrec(abs,compon)$asolv(abs,compon)..
sum(str$ovabs(abs,str), fc(str,compon))*exp(beta(abs,compon)) =e=
sum(str$ivabs(abs,str), fc(str,compon)) +
exp(gamma(abs,compon))*sum(str$ilabs(abs,str), fc(str,compon));
absrec(abs,compon)$anorm(abs,compon)..
sum(str$ovabs(abs,str), fc(str,compon))*exp(beta(abs,compon)) =e=
sum(str$ivabs(abs,str), fc(str,compon));
abssimp(abs,compon)$asimp(abs,compon)..
sum(str$ovabs(abs,str), fc(str,compon)) =e=
sum(str$ivabs(abs,str), fc(str,compon))/cbeta(compon);
abscmb(abs,compon)..
sum(str$ilabs(abs,str), fc(str,compon)) +
sum(str$ivabs(abs,str), fc(str,compon)) =e=
sum(str$olabs(abs,str), fc(str,compon)) +
sum(str$ovabs(abs,str), fc(str,compon));
abslogc..
f('63') + f('67') =l= 25.*y('63');
abspl(abs)..
sum(str$ilabs(abs,str), p(str)) =e=
sum(str$olabs(abs,str), p(str));
abstl(abs)..
sum(str$ilabs(abs,str), t(str)) =e= sum(str$olabs(abs,str), t(str));
abspv(abs)..
sum(str$ivabs(abs,str), p(str)) =e= sum(str$ovabs(abs,str), p(str));
absp(abs)..
sum(str$ilabs(abs,str), p(str)) =e= sum(str$ivabs(abs,str), p(str));
abst(abs)..
sum(str$ilabs(abs,str), t(str)) =e= sum(str$ovabs(abs,str), t(str));
* ===========================
* = =
* = c o m p r e s s o r =
* = =
* ===========================
Equation
compcmb(comp,compon) 'component balance for compressor'
comphb(comp) 'heat balance for compressor'
compelec(comp) 'energy balance for compressor'
ratio(comp) 'pressure ratio (out to in)';
compcmb(comp,compon)..
sum(str$ocomp(comp,str), fc(str,compon)) =e= sum(str$icomp(comp,str), fc(str,compon));
comphb(comp)..
sum(str$ocomp(comp,str), t(str)) =e= presrat(comp)*sum(str$icomp(comp,str), t(str));
compelec(comp)..
elec(comp) =e= alpha*(presrat(comp) - 1)
* sum(str$icomp(comp,str), 100.*t(str)*f(str)/60.)
* (1/compeff)*(gam/(gam - 1.));
ratio(comp)..
presrat(comp)**(gam/(gam - 1.)) =e= sum(str$ocomp(comp,str), p(str))
/ sum(str$icomp(comp,str), p(str));
* ===============================
* = =
* = d i s t i l l a t i o n =
* = =
* ===============================
Equation
distcmb(dist,compon) 'component mass balance'
antdistb(dist,str,compon) 'vapor pressure correlation (bot)'
antdistt(dist,str,compon) 'vapor pressure correlation (top)'
relvol(dist) 'average relative volatilty'
undwood(dist) 'minimum reflux ratio equation'
actreflux(dist) 'actual reflux ratio'
fenske(dist) 'minimum number of trays'
acttray(dist) 'actual number of trays'
distspec(dist,str,compon) 'recovery specification'
distheav(dist,compon) 'heavy components'
distlite(dist,compon) 'light components'
distpi(dist,str) 'inlet pressure relation'
distvpl(dist,str) 'bottom vapor pressure relation'
distvpv(dist,str) 'top vapor pressure relation'
distpl(dist,str) 'outlet pressure relation(liquid)'
distpv(dist,str) 'outlet pressure relation(vapor)';
distcmb(dist,compon)..
sum(str$idist(dist,str), fc(str,compon)) =e=
sum(str$vdist(dist,str), fc(str,compon)) +
sum(str$ldist(dist,str), fc(str,compon));
antdistb(dist,str,compon)$(ldist(dist,str)$dkey(dist,compon))..
log(vp(str,compon)*7500.6168) =e= anta(compon) - antb(compon)/(t(str)*100. + antc(compon));
antdistt(dist,str,compon)$(vdist(dist,str)$dkey(dist,compon))..
log(vp(str,compon)*7500.6168) =e= anta(compon) - antb(compon)/(t(str)*100. + antc(compon));
relvol(dist)..
avevlt(dist) =e= sqrt(sum(str$ ldist(dist,str) ,
sum(compon$dlkey(dist,compon), vp(str,compon)) /
sum(compon$dhkey(dist,compon), vp(str,compon))) *
sum(str$ vdist(dist,str) ,
sum(compon$dlkey(dist,compon), vp(str,compon)) /
sum(compon$dhkey(dist,compon), vp(str,compon))));
undwood(dist)..
sum(str$idist(dist,str), sum(compon$dlkey(dist,compon), fc(str,compon)))*
rmin(dist)*(avevlt(dist) - 1) =e= sum(str$idist(dist,str), sum(compon$dlkey(dist,compon), f(str)));
actreflux(dist).. reflux(dist) =e= 1.2*rmin(dist);
fenske(dist)..
nmin(dist)*log(avevlt(dist)) =e=
log(sum(str$vdist(dist,str),
sum(compon$dhkey(dist,compon), (f(str) + eps1)/(fc(str,compon) + eps1)))*
sum(str$ldist(dist,str),
sum(compon$dlkey(dist,compon), (f(str) + eps1)/(fc(str,compon) + eps1))));
acttray(dist)..
ndist(dist) =e= nmin(dist)*2./disteff;
distspec(dist,str,compon)$(vdist(dist,str)$dhkey(dist,compon))..
fc(str,compon) =l= 0.05*sum(str2$idist(dist,str2), fc(str2,compon));
distheav(dist,compon)$dh(dist,compon)..
sum(str$idist(dist,str), fc(str,compon)) =e= sum(str$ldist(dist,str), fc(str,compon));
distlite(dist,compon)$dl(dist,compon)..
sum(str$idist(dist,str), fc(str,compon)) =e= sum(str$vdist(dist,str), fc(str,compon));
distpi(dist,str)$idist(dist,str)..
distp(dist) =l= p(str);
distvpl(dist,str)$ldist(dist,str)..
distp(dist) =e= sum(compon$dhkey(dist,compon), vp(str,compon));
distvpv(dist,str)$((ord(dist) <> 1)$vdist(dist,str))..
distp(dist) =e= sum(compon$dlkey(dist,compon), vp(str,compon));
distpl(dist,str)$ldist(dist,str)..
distp(dist) =e= p(str);
distpv(dist,str) $ vdist(dist,str)..
distp(dist) =e= p(str);
* =================
* = =
* = f l a s h =
* = =
* =================
Equation
flshcmb(flsh,compon) 'component mass balance'
flshpr(flsh,str) 'flash pressure relation'
flshpi(flsh,str) 'inlet pressure relation'
flshpl(flsh,str) 'outlet pressure relation(liquid)'
flshpv(flsh,str) 'outlet pressure relation(vapor)'
flshti(flsh,str) 'inlet temp. relation'
flshtl(flsh,str) 'outlet temp. relation(liquid)'
flshtv(flsh,str) 'outlet temp. relation(vapor)'
flshrec(flsh,str,compon) 'vapor recovery relation'
flsheql(flsh,compon) 'equilibrium relation'
antflsh(flsh,str,compon) 'vapor pressure correlation';
flshcmb(flsh,compon)..
sum(str$iflsh(flsh,str), fc(str,compon)) =e=
sum(str$vflsh(flsh,str), fc(str,compon)) +
sum(str$lflsh(flsh,str), fc(str,compon));
antflsh(flsh,str,compon)$lflsh(flsh,str)..
log(vp(str,compon)*7500.6168) =e= anta(compon) - antb(compon)/(t(str)*100. + antc(compon));
flshrec(flsh,str,compon)$lflsh(flsh,str)..
sum(compon2$fkey(flsh,compon2), eflsh(flsh,compon2))*
(eflsh(flsh,compon)*sum(compon2$fkey(flsh,compon2), vp(str,compon2)) +
(1. - eflsh(flsh,compon))*vp(str,compon)) =e=
sum(compon2$fkey(flsh,compon2), vp(str,compon2))*eflsh(flsh,compon);
flsheql(flsh,compon)..
sum(str$vflsh(flsh,str), fc(str,compon)) =e=
sum(str$iflsh(flsh,str), fc(str,compon))*eflsh(flsh,compon);
flshpr(flsh,str)$lflsh(flsh,str)..
flshp(flsh)*f(str) =e= sum(compon, vp(str,compon)*fc(str,compon));
flshpi(flsh,str)$iflsh(flsh,str)..
flshp(flsh) =e= p(str);
flshpl(flsh,str)$lflsh(flsh,str)..
flshp(flsh) =e= p(str);
flshpv(flsh,str)$vflsh(flsh,str)..
flshp(flsh) =e= p(str);
flshti(flsh,str)$iflsh(flsh,str)..
flsht(flsh) =e= t(str);
flshtl(flsh,str)$lflsh(flsh,str)..
flsht(flsh) =e= t(str);
flshtv(flsh,str)$vflsh(flsh,str)..
flsht(flsh) =e= t(str);
* =======================
* = =
* = f u r n a c e =
* = =
* =======================
Equation
furnhb(furn) 'heat balance'
furncmb(furn,compon) 'component mass balance'
furnp(furn) 'pressure relation';
furnhb(furn)..
qfuel(furn) =e= (sum(str$ofurn(furn,str), cp(str)*f(str)*100.*t(str))
-sum(str$ifurn(furn,str), cp(str)*f(str)*100.*t(str)))
*3600.*8500.*1.0e-12/60.;
furncmb(furn,compon)..
sum(str$ofurn(furn,str), fc(str,compon)) =e= sum(str$ifurn(furn,str), fc(str,compon));
furnp(furn)..
sum(str$ofurn(furn,str), p(str)) =e= sum(str$ifurn(furn,str), p(str)) - 0.4826;
* ===================
* = =
* = c o o l e r =
* = =
* ===================
Equation
heccmb(hec,compon) 'component balance in cooler'
hechb(hec) 'heat balance for cooler'
hecp(hec) 'no pressure drop thru cooler';
heccmb(hec,compon)..
sum(str$ohec(hec,str), fc(str,compon)) =e= sum(str$ihec(hec,str), fc(str,compon));
hechb(hec)..
qc(hec) =e= (sum(str$ihec(hec,str), cp(str)*f(str)*100.*t(str))
-sum(str$ohec(hec,str), cp(str)*f(str)*100.*t(str)))
*3600.*8500.*1.0e-12/60.;
hecp(hec)..
sum(str$ihec(hec,str), p(str)) =e= sum(str$ohec(hec,str), p(str));
* ===================
* = =
* = h e a t e r =
* = =
* ===================
Equation
hehcmb(heh,compon) 'component balance in heater'
hehhb(heh) 'heat balance for heater'
hehp(heh) 'no pressure drop thru heater';
hehcmb(heh,compon)..
sum(str$oheh(heh,str), fc(str,compon)) =e= sum(str$iheh(heh,str), fc(str,compon));
hehhb(heh)..
qh(heh) =e= (sum(str$oheh(heh,str), cp(str)*f(str)*100.*t(str))
-sum(str$iheh(heh,str), cp(str)*f(str)*100.*t(str)))
*3600.*8500.*1.0e-12/60.;
hehp(heh).. sum(str$iheh(heh,str), p(str)) =e= sum(str$oheh(heh,str), p(str));
* =========================
* = =
* = e x c h a n g e r =
* = =
* =========================
Equation
exchcmbc(exch,compon) 'component balance (cold)'
exchcmbh(exch,compon) 'component balance (hot)'
exchhbc(exch) 'heat balance for cold stream'
exchhbh(exch) 'heat balance for hot stream'
exchdtm1(exch) 'delta t min condition'
exchdtm2(exch) 'delta t min condition'
exchpc(exch) 'pressure relation (cold)'
exchph(exch) 'pressure relation (hot)';
exchcmbc(exch,compon)..
sum(str$ocexch(exch,str), fc(str,compon)) =e= sum(str$icexch(exch,str), fc(str,compon));
exchcmbh(exch,compon)..
sum(str$ohexch(exch,str), fc(str,compon)) =e= sum(str$ihexch(exch,str), fc(str,compon));
exchhbc(exch)..
qexch(exch) =e= (sum(str$ocexch(exch,str), cp(str)*f(str)*100.*t(str))
-sum(str$icexch(exch,str), cp(str)*f(str)*100.*t(str)))
*3600.*8500.*1.0e-12/60.;
exchhbh(exch)..
qexch(exch) =e= (sum(str$ihexch(exch,str), cp(str)*f(str)*100.*t(str))
-sum(str$ohexch(exch,str), cp(str)*f(str)*100.*t(str)))
*3600.*8500.*1.0e-12/60.;
exchdtm1(exch)..
sum(str$ohexch(exch,str), t(str)) =g= sum(str$icexch(exch,str), t(str)) + 0.25;
exchdtm2(exch)..
sum(str$ocexch(exch,str), t(str)) =l= sum(str$ihexch(exch,str), t(str)) - 0.25;
exchpc(exch)..
sum(str$ocexch(exch,str), p(str)) =e= sum(str$icexch(exch,str), p(str));
exchph(exch)..
sum(str$ohexch(exch,str), p(str)) =e= sum(str$ihexch(exch,str), p(str));
* =======================
* = =
* = m e m b r a n e =
* = =
* =======================
Equation
memcmb(memb,str,compon) 'component mass balance'
flux(memb,str,compon) 'mass flux relation'
simp(memb,str,compon) 'mass flux relation (simplified)'
memtp(memb,str) 'temp relation for permeate'
mempp(memb,str) 'pressure relation for permeate'
memtn(memb,str) 'temp relation for non-permeate'
mempn(memb,str) 'pressure relation for non-permeate'
rec(memb,str) 'recovery spec'
pure(memb,str) 'purity spec';
memcmb(memb,str,compon)$imemb(memb,str)..
fc(str,compon) =e= sum(str2$pmemb(memb,str2), fc(str2,compon))
+ sum(str2$nmemb(memb,str2), fc(str2,compon));
flux(memb,str,compon)$(pmemb(memb,str)$mnorm(memb,compon))..
fc(str,compon) =e= a(memb)*perm(compon)/2.0
*(sum(str2$imemb(memb,str2), p(str2))
*(sum(str2$imemb(memb,str2), (fc(str2,compon) + eps1)/(f(str2) + eps1)) +
sum(str2$nmemb(memb,str2), (fc(str2,compon) + eps1)/(f(str2) + eps1))) -
2.*p(str)*(fc(str,compon) + eps1)/(f(str) + eps1));
simp(memb,str,compon)$(pmemb(memb,str)$msimp(memb,compon))..
fc(str,compon) =e= 0.0;
memtp(memb,str)$pmemb(memb,str)..
t(str) =e= sum(str2$imemb(memb,str2), t(str2));
mempp(memb,str)$pmemb(memb,str)..
p(str) =l= sum(str2$imemb(memb,str2), p(str2));
memtn(memb,str)$nmemb(memb,str)..
t(str) =e= sum(str2$imemb(memb,str2), t(str2));
mempn(memb,str)$nmemb(memb,str)..
p(str) =e= sum(str2$imemb(memb,str2), p(str2));
rec(memb,str)$pmemb(memb,str)..
fc(str,'h2') =g= 0.50*sum(str2$imemb(memb,str2), fc(str2,'h2'));
pure(memb,str)$pmemb(memb,str)..
fc(str,'h2') =g= 0.50*f(str);
* ==================
* = =
* = m i x e r =
* = (single input) =
* = =
* ==================
Equation
mxr1cmb(mxr1,str,compon) 'component balance in mixer'
mxr1hbl(mxr1,str) 'heat balance in mixer (ls than)'
mxr1hbg(mxr1,str) 'heat balance in mixer (grt than)'
mxr1pl(mxr1,str) 'pressure relation (less than)'
mxr1pg(mxr1,str) 'pressure relation (greater than)';
mxr1cmb(mxr1,str,compon)$omxr1(mxr1,str)..
fc(str,compon) =e= sum(str2$imxr1(mxr1,str2), fc(str2,compon));
mxr1hbl(mxr1,str)$imxr1(mxr1,str)..
sum(str2$omxr1(mxr1,str2), t(str2)) =l= t(str) + utemp*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));
mxr1hbg(mxr1,str)$imxr1(mxr1,str)..
sum(str2$omxr1(mxr1,str2), t(str2)) =g= t(str) - utemp*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));
mxr1pl(mxr1,str)$imxr1(mxr1,str)..
sum(str2$omxr1(mxr1,str2), p(str2)) =l= p(str) + upress*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));
mxr1pg(mxr1,str)$imxr1(mxr1,str)..
sum(str2$omxr1(mxr1,str2), p(str2)) =g= p(str) - upress*sum(str2$mxr1spl1(mxr1,str,str2), (1 - y(str2)));
* =================
* = =
* = m i x e r =
* = =
* =================
Equation
mxrcmb(mxr,compon) 'component balance in mixer'
mxrhb(mxr) 'heat balance in mixer'
mxrhbq(mxr) 'heat balance in quench'
mxrpi(mxr,str) 'inlet pressure relation'
mxrpo(mxr,str) 'outlet pressure relation';
mxrcmb(mxr,compon)..
sum(str$omxr(mxr,str), fc(str,compon)) =e= sum(str$imxr(mxr,str), fc(str,compon));
mxrhb(mxr)$(ord(mxr) <> 2)..
sum(str$imxr(mxr,str), f(str)*t(str)*cp(str)) =e= sum(str$omxr(mxr,str), f(str)*t(str)*cp(str));
mxrhbq(mxr)$(ord(mxr) = 2)..
f('16')*t('16') =e= f('15')*t('15') - (fc('20','ben') + fc('20','tol'))*heatvap('tol')/(100.* cp('15'));
mxrpi(mxr,str)$imxr(mxr,str)..
mxrp(mxr) =e= p(str);
mxrpo(mxr,str)$omxr(mxr,str)..
mxrp(mxr) =e= p(str);
* ===============
* = =
* = p u m p =
* = =
* ===============
Equation
pumpcmb(pump,compon) 'component balance'
pumphb(pump) 'heat balance'
pumppr(pump) 'pressure relation';
pumpcmb(pump,compon)..
sum(str$opump(pump,str), fc(str,compon)) =e= sum(str$ipump(pump,str), fc(str,compon));
pumphb(pump)..
sum(str$opump(pump,str), t(str)) =e= sum(str$ipump(pump,str), t(str));
pumppr(pump)..
sum(str$opump(pump,str), p(str)) =g= sum(str$ipump(pump,str), p(str));
* =====================
* = =
* = r e a c t o r =
* = =
* =====================
Equation
rctspec(rct,str) 'spec. on reactor feed stream'
rxnrate(rct) 'reaction rate constant'
rctconv(rct,str,compon) 'conversion of key component'
rctsel(rct) 'selectivity to benzene'
rctcns(rct,str,compon) 'consumption rate of key comp.'
rctmbtol(rct) 'mass balance in reactor (tol)'
rctmbben(rct) 'mass balance in reactor (ben)'
rctmbdip(rct) 'mass balance in reactor (dip)'
rctmbh2(rct) 'mass balance in reactor (h2)'
rctmbch4(rct) 'mass balance in reactor (ch4)'
rcthbadb(rct) 'heat balance (adiabatic)'
rcthbiso(rct) 'heat balance (isothermal)'
rctisot(rct) 'temp relation (isothermal)'
rctpi(rct,str) 'inlet pressure relation'
rctpo(rct,str) 'outlet pressure relation'
rcttave(rct) 'average temperature relation';
rctspec(rct,str)$irct(rct,str)..
fc(str,'h2') =g= 5.*(fc(str,'ben') + fc(str,'tol') + fc(str,'dip'));
rxnrate(rct)..
krct(rct) =e= 6.3e+10*exp(-26167./(rctt(rct)*100.));
rctconv(rct,str,compon)$(rkey(rct,compon)$irct(rct,str))..
1. - conv(rct,compon) =e= sqr(1./(1. + 0.372*krct(rct)*rctvol(rct)*sqrt(fc(str,compon)
/60 + eps1)*(f(str)/60. + eps1)**(-3./2.)));
rctsel(rct)..
(1. - sel(rct)) =e= 0.0036*(1. - conv(rct,'tol'))**(-1.544);
rctcns(rct,str,compon)$(rkey(rct,compon)$irct(rct,str))..
consum(rct,compon) =e= conv(rct,compon)*fc(str,compon);
rctmbtol(rct)..
sum(str$orct(rct,str), fc(str,'tol')) =e= sum(str$irct(rct,str), fc(str,'tol')) - consum(rct,'tol');
rctmbben(rct)..
sum(str$orct(rct,str), fc(str,'ben')) =e= sum(str$irct(rct,str), fc(str,'ben')) + consum(rct,'tol')*sel(rct);
rctmbdip(rct)..
sum(str$orct(rct,str), fc(str,'dip')) =e=
sum(str$irct(rct,str), fc(str,'dip')) +
consum( rct,'tol')*0.5 +
(sum(str$irct(rct,str), fc(str,'ben')) -
sum(str$orct(rct,str), fc(str,'ben')))*0.5;
rctmbh2(rct)..
sum(str$orct(rct,str), fc(str,'h2')) =e=
sum(str$irct(rct,str), fc(str,'h2')) -
consum(rct,'tol') +
sum(str$orct(rct,str), fc(str,'dip')) -
sum(str$irct(rct,str), fc(str,'dip'));
rctmbch4(rct)..
sum(str$orct(rct,str), fc(str,'ch4')) =e=
sum(str$irct(rct,str), fc(str,'ch4')) +
consum(rct,'tol');
rcthbadb(rct)$(ord(rct) = 1)..
heatrxn(rct)*consum(rct,'tol')/100. =e=
sum(str$orct(rct,str), cp(str)*f(str)*t(str)) -
sum(str$irct(rct,str), cp(str)*f(str)*t(str));
rcthbiso(rct)$(ord(rct) = 2)..
heatrxn(rct)*consum(rct,'tol')*60.*8500*1.0e-09 =e= q(rct);
rctisot(rct)$(ord(rct) = 2)..
sum(str$orct(rct,str), t(str)) =e= sum(str$irct(rct,str), t(str));
rctpi(rct,str)$irct(rct,str)..
rctp(rct) =e= p(str);
rctpo(rct,str)$orct(rct,str)..
p(str) =e= rctp(rct) - 0.20684;
rcttave(rct)..
rctt(rct) =e= (sum(str$irct(rct,str), t(str)) +
sum(str$orct(rct,str), t(str)))/2.;
* =======================
* = =
* = s p l i t t e r =
* = (single output) =
* = =
* =======================
Equation
spl1logc(spl1,str) 'splitter logical constraint'
spl1ysum(spl1) 'single choice'
spl1cmb(spl1,compon) 'component balance in splitter'
spl1pi(spl1,str) 'inlet pressure relation'
spl1po(spl1,str) 'outlet pressure relation'
spl1ti(spl1,str) 'inlet temperature relation'
spl1to(spl1,str) 'outlet temperature relation';
spl1logc(spl1,str)$ospl1(spl1,str)..
f(str) =l= uflow*y(str);
spl1ysum(spl1)..
sum(str$ospl1(spl1,str), y(str)) =e= 1;
spl1cmb(spl1,compon)..
sum(str$ospl1(spl1,str), fc(str,compon)) =e= sum(str$ispl1(spl1,str), fc(str,compon));
spl1pi(spl1,str)$ispl1(spl1,str)..
spl1p(spl1) =e= p(str);
spl1po(spl1,str)$ospl1(spl1,str)..
spl1p(spl1) =e= p(str);
spl1ti(spl1,str)$ispl1(spl1,str)..
spl1t(spl1) =e= t(str);
spl1to(spl1,str)$ospl1(spl1,str)..
spl1t(spl1) =e= t(str);
* =======================
* = =
* = s p l i t t e r =
* = =
* =======================
Equation
splcmb(spl,str,compon) 'component balance in splitter'
esum(spl) 'split fraction relation'
splpi(spl,str) 'inlet pressure relation'
splpo(spl,str) 'outlet pressure relation'
splti(spl,str) 'inlet temperature relation'
splto(spl,str) 'outlet temperature relation';
splcmb(spl,str,compon)$ospl(spl,str)..
fc(str,compon) =e= sum(str2$ispl(spl,str2), fc(str2,compon)*e(str));
esum(spl)..
sum(str$ospl(spl,str), e(str)) =e= 1.0;
splpi(spl,str)$ispl(spl,str)..
splp(spl) =e= p(str);
splpo(spl,str)$ospl(spl,str)..
splp(spl) =e= p(str);
splti(spl,str)$ispl(spl,str)..
splt(spl) =e= t(str);
splto(spl,str)$ospl(spl,str)..
splt(spl) =e= t(str);
* ==========================
* = =
* = e x p v a l v e =
* = =
* ==========================
Equation
valcmb(valve,compon) 'component mass balance'
valt(valve) 'temperature relation'
valp(valve) 'pressure relation';
valcmb(valve,compon)..
sum(str$oval(valve,str), fc(str,compon)) =e= sum(str$ival(valve,str), fc(str,compon));
valt(valve)..
sum(str$oval(valve,str), t(str)/(p(str)**((gam - 1.)/gam))) =e=
sum(str$ival(valve,str), t(str)/(p(str)**((gam - 1.)/gam)));
valp(valve)..
sum(str$oval(valve,str), p(str)) =l= sum(str$ival(valve,str), p(str));
* ===================
* = =
* = s t r e a m =
* = =
* ===================
Equation
h2feed(compon) 'define h2 feed component flows'
tolfeed(compon) 'define tol feed component flows'
tolabs(compon) 'define tol absorber compon flows'
specrec 'spec on h2 recycle'
specprod 'spec on ben product'
fbal(str) 'stream flow = sum of comp. flows';
specrec..
fc('72','h2') =g= 0.5*f('72');
specprod..
fc('31','ben') =g= 0.9997*f('31');
fbal(str)..
f(str) =e= sum(compon, fc(str,compon));
h2feed(compon)..
fc('1',compon) =e= f('1')*f1comp(compon);
tolfeed(compon)..
fc('66',compon) =e= f('66')*f66comp(compon);
tolabs(compon)..
fc('67',compon) =e= f('67')*f67comp(compon);
* =========================
* = =
* = o b j e c t i v e =
* = =
* =========================
Equation
abound 'fake bound on a(1)'
obj 'cost function definition';
abound..
a('1') =g= 0.0;
obj..
profit =e= 510.*(- 2.5 *f('1') - 14. *(f('66') + f('67'))
+ 19.9*f('31') + 11.84* f('35')
+ 1.08*(fc('4','h2') + fc('28','h2') +
fc('53','h2') + fc('55','h2'))
+ 3.37*(fc('4','ch4') + fc('28','ch4') +
fc('53','ch4') + fc('55','ch4')))
- 0.815*(elec('1') + elec('2') + elec('3'))
- 0.887* elec('4')
- 7.155*(y('3') + y('59') + y('63'))
- 4.866* y('54')
- sum(comp, costelec*elec(comp))
- (74.3*y('10') + 1.257*rctvol('1') )
- 1.25*(74.3*y('12') + 1.257*rctvol('2')) - 0.7*q('2')
- (1.126*y('24') + 0.375*ndist('1')) - (+ 1.55*ndist('2'))
- (3.9*y('33') + 1.12*ndist('3'))
- (43.24*y('3') + 49.0*f('3'))
- (43.24*y('54') + 49.0*f('54'))
- (3.0 *y('63') + 1.2 *nabs('1'))
- ( + 1171.7 * qfuel('1') + 4000.*qfuel('1'))
- sum(hec, 700 * qc(hec))
- sum(heh, 8000* qh(heh))
- const;
* + + + + + + + + + + + + + + + +
* + +
* + b o u n d s +
* + +
* + + + + + + + + + + + + + + + +
profit.up = 1.0e+08;
const.fx = 22.5;
* ===================
* = s t r e a m =
* ===================
y.lo(str) = 0.0;
y.up(str) = 1.0;
e.up(str) = 1.0;
e.up('20') = 0.50;
e.lo('21') = 0.50;
* e.up('52') = 0.50;
* e.lo('58') = 0.50;
f.lo(str) = 0.;
f.up(str) = 10.;
f.up('8') = 50.;
f.up('9') = 50.;
f.up('10') = 50.;
f.up('11') = 50.;
f.up('12') = 50.;
f.up('13') = 50.;
f.up('14') = 50.;
f.up('15') = 50.;
f.up('16') = 50.;
f.up('17') = 50.;
f.up('18') = 50.;
f.up('52') = 50.;
f.up('54') = 50.;
f.up('56') = 50.;
f.up('57') = 50.;
f.up('58') = 50.;
f.up('59') = 50.;
f.up('60') = 50.;
f.up('70') = 50.;
f.up('71') = 50.;
f.up('72') = 50.;
fc.lo(str,compon) = 0.;
fc.up(str,compon) = 10.;
fc.up('8',compon) = 30.;
fc.up('9',compon) = 30.;
fc.up('10',compon) = 30.;
fc.up('11',compon) = 30.;
fc.up('12',compon) = 30.;
fc.up('13',compon) = 30.;
fc.up('14',compon) = 30.;
fc.up('15',compon) = 30.;
fc.up('16',compon) = 30.;
fc.up('17',compon) = 30.;
fc.up('18',compon) = 30.;
fc.up('52',compon) = 30.;
fc.up('54',compon) = 30.;
fc.up('56',compon) = 30.;
fc.up('57',compon) = 30.;
fc.up('58',compon) = 30.;
fc.up('59',compon) = 30.;
fc.up('60',compon) = 30.;
fc.up('70',compon) = 30.;
fc.up('71',compon) = 30.;
fc.up('72',compon) = 30.;
p.lo(str) = 0.1;
p.up(str) = 4.0;
t.lo(str) = 3.00;
t.up(str) = 10.00;
t.lo('49') = 2.0;
t.lo('50') = 2.0;
t.lo('51') = 2.0;
vp.lo(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo(str)*100. + antc(compon)));
vp.up(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up(str)*100. + antc(compon)));
* =======================
* = a b s o r b e r =
* =======================
nabs.lo(abs) = 0.;
nabs.up(abs) = 40.;
gamma.lo(abs,'tol') = log((1 - aabs('tol')**( nabs.lo('1')*abseff + eps1))/(1 - aabs('tol')));
gamma.up(abs,'tol') = min(15, log((1 - aabs('tol')**( nabs.up('1')*abseff + eps1))/(1 - aabs('tol'))));
beta.lo(abs,compon) = log((1 - aabs(compon)**( nabs.lo('1')*abseff + eps1 + 1.))/(1 - aabs(compon)));
beta.up(abs,compon) = min(15, log((1 - aabs(compon)**( nabs.up('1')*abseff + eps1 + 1.))/(1 - aabs(compon))));
t.fx('67') = 3.0;
vp.fx('67',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.l('67')*100. + antc(compon)));
* ===========================
* = c o m p r e s s o r =
* ===========================
presrat.lo(comp) = 1.;
presrat.up(comp) = 8./3.;
elec.lo(comp)= 0.;
elec.up(comp)= 100.;
* ===============================
* = d i s t i l l a t i o n =
* ===============================
nmin.lo('1') = 0.;
nmin.up('1') = 4.;
nmin.lo('2') = 8.;
nmin.up('2') = 14.;
nmin.lo('3') = 0.;
nmin.up('3') = 4.;
ndist.lo(dist) = nmin.lo(dist)*2./disteff;
ndist.up(dist) = nmin.up(dist)*2./disteff;
rmin.lo('1') = 0.02;
rmin.up('1') = 0.10;
rmin.lo('2') = 0.50;
rmin.up('2') = 2.00;
rmin.lo('3') = 0.02;
rmin.up('3') = 0.10;
reflux.lo(dist) = rmin.lo(dist)*1.2;
reflux.up(dist) = rmin.up(dist)*1.2;
distp.lo(dist) = 0.1;
distp.up(dist) = 4.0;
* fix pressure at 1.02 mpa (150 psia)
distp.fx('1') = 1.020;
distp.up('2') = 0.400;
distp.up('3') = 0.250;
* fix top temperature at 320 k
t.fx('26') = 3.20;
t.lo('27') = (antb('ben')/(anta('ben') - log(distp.lo('1')*7500.6168)) - antc('ben'))/100.;
t.up('27') = (antb('ben')/(anta('ben') - log(distp.up('1')*7500.6168)) - antc('ben'))/100.;
t.lo('31') = (antb('ben')/(anta('ben') - log(distp.lo('2')*7500.6168)) - antc('ben'))/100.;
t.up('31') = (antb('ben')/(anta('ben') - log(distp.up('2')*7500.6168)) - antc('ben'))/100.;
t.lo('32') = (antb('tol')/(anta('tol') - log(distp.lo('2')*7500.6168)) - antc('tol'))/100.;
t.up('32') = (antb('tol')/(anta('tol') - log(distp.up('2')*7500.6168)) - antc('tol'))/100.;
t.lo('34') = (antb('tol')/(anta('tol') - log(distp.lo('3')*7500.6168)) - antc('tol'))/100.;
t.up('34') = (antb('tol')/(anta('tol') - log(distp.up('3')*7500.6168)) - antc('tol'))/100.;
t.lo('35') = (antb('dip')/(anta('dip') - log(distp.lo('3')*7500.6168)) - antc('dip'))/100.;
t.up('35') = (antb('dip')/(anta('dip') - log(distp.up('3')*7500.6168)) - antc('dip'))/100.;
vp.lo(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo(str)*100. + antc(compon)));
vp.up(str,compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up(str)*100. + antc(compon)));
avevlt.lo(dist) = sum(compon$dlkey(dist,compon),
sum(str $ldist(dist,str) ,vp.up(str,compon)))/
sum(compon$dhkey(dist,compon),
sum(str $ldist(dist,str) ,vp.up(str,compon)));
avevlt.up(dist) = sum(compon$dlkey(dist,compon),
sum(str $vdist(dist,str) , vp.lo(str,compon)))/
sum(compon$dhkey(dist,compon),
sum(str $vdist(dist,str) , vp.lo(str,compon)));
* ==================
* = f l a s h =
* ==================
eflsh.up(flsh,compon) = 1.;
eflsh.lo(flsh,compon) = 0.;
eflsh.lo('1','h2') = 0.975;
eflsh.lo('1','ch4') = 0.975;
eflsh.up('1','ben') = 0.075;
eflsh.up('1','tol') = 0.050;
eflsh.up('1','dip') = 0.001;
eflsh.lo('2','h2') = 0.975;
eflsh.lo('2','ch4') = 0.975;
eflsh.up('2','ben') = 0.20;
eflsh.up('2','tol') = 0.10;
eflsh.up('2','dip') = 0.01;
eflsh.lo('3','h2') = 0.99;
eflsh.lo('3','ch4') = 0.99;
eflsh.lo('3','ben') = 0.90;
eflsh.lo('3','tol') = 0.90;
eflsh.up('3','dip') = 0.40;
flshp.lo('1') = 2.5;
flshp.up('1') = 3.25;
flshp.lo('2') = 0.1;
flshp.up('2') = 3.25;
flshp.lo('3') = 0.1;
flshp.up('3') = 3.25;
flsht.lo('1') = 3.00;
flsht.up('1') = 4.50;
flsht.lo('2') = 3.00;
flsht.up('2') = 4.50;
flsht.lo('3') = 3.00;
flsht.up('3') = 6.50;
t.lo('19') = flsht.lo('1');
t.up('19') = flsht.up('1');
vp.lo('19',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo('19')*100. + antc(compon)));
vp.up('19',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up('19')*100. + antc(compon)));
t.lo('48') = flsht.lo('2');
t.up('48') = flsht.up('2');
vp.lo('48',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo('48')*100. + antc(compon)));
vp.up('48',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up('48')*100. + antc(compon)));
*t.lo('41') = flsht.lo('3');
t.lo('41') = t.lo('32');
t.up('41') = flsht.up('3');
vp.lo('41',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.lo('41')*100. + antc(compon)));
vp.up('41',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.up('41')*100. + antc(compon)));
* ======================
* = f u r n a c e =
* ======================
qfuel.up(furn) = 10.;
* ====================
* = c o o l e r =
* ====================
qc.up(hec) = 10.;
* ====================
* = h e a t e r =
* ====================
qh.up(heh) = 10.;
* ==========================
* = e x c h a n g e r =
* ==========================
qexch.up(exch) = 10.;
* =======================
* = m e m b r a n e =
* =======================
a.lo(memb) = 100;
a.up(memb) = 10000;
* ==================
* = m i x e r =
* = (1 output) =
* ==================
mxr1t.lo(mxr1) = 3.0;
mxr1t.up(mxr1) = 10.0;
mxr1p.lo(mxr1) = 0.1;
mxr1p.up(mxr1) = 4.0;
* ==================
* = m i x e r =
* ==================
mxrt.lo(mxr) = 3.0;
mxrt.up(mxr) = 10.0;
mxrp.lo(mxr) = 0.1;
mxrp.up(mxr) = 4.0;
* =====================
* = r e a c t o r =
* =====================
rctt.lo(rct) = 8.9427;
rctt.up(rct) = 9.7760;
rctp.fx(rct) = 3.4474;
rctvol.up(rct) = 200.;
krct.lo(rct) = 6.3e+10*exp(-26167./(rctt.lo(rct)*100.));
krct.up(rct) = 6.3e+10*exp(-26167./(rctt.up(rct)*100.));
conv.up(rct,'tol') = 0.973;
sel.up(rct) = 1.0 - 0.0036;
q.up('2') = 100.;
* =======================
* = s p l i t t e r =
* = (1 outlet) =
* =======================
spl1t.lo(spl1) = 3.00;
spl1t.up(spl1) = 10.00;
spl1p.lo(spl1) = 0.1;
spl1p.up(spl1) = 4.0;
* =======================
* = s p l i t t e r =
* =======================
splt.lo(spl) = 3.00;
splt.up(spl) = 10.00;
splp.lo(spl) = 0.1;
splp.up(spl) = 4.0;
* =========================
* = o t h e r s =
* =========================
f.fx('31') = 2.08;
t.lo('9') = 8.943;
t.up('16') = 8.943;
t.fx('1') = 3.0;
p.fx('1') = 3.930;
p.fx('66') = 3.930;
t.fx('66') = 3.00;
* + + + + + + + + + + + + + + + + + +
* + +
* + i n i t i a l +
* + +
* + p o i n t +
* + +
* + + + + + + + + + + + + + + + + + +
* ---------------------------------------
* _ b i n a r y v a r i a b l e s _
* _______________________________________
y.l(str) = 0.0;
y.l('2') = 1.0;
*y.l('10') = 1.0;
y.l('12') = 1.0;
y.l('24') = 1.0;
*y.l('33') = 1.0;
y.l('37') = 1.0;
*y.l('53') = 1.0;
y.l('54') = 1.0;
y.l('59') = 1.0;
*y.fx(str) = y.l(str);
* -----------------------------------------------
* _ c o n t i n u o u s v a r i a b l e s _
* _______________________________________________
t.l(str) = 3.0;
p.l(str) = 3.0;
f.l('1') = 4.5 - 2.0*y.l('54');
fc.l('1',compon) = f.l('1')*f1comp(compon);
* ---------------
* - b y p a s s -
* ---------------
fc.l('2',compon) = fc.l('1',compon)*y.l('2');
t.l('2') = t.l('1');
p.l('2') = p.l('1');
* -----------------------
* - m e m b r a n e 1 -
* -----------------------
fc.l('3',compon) = fc.l('1',compon)*y.l('3');
fc.l('4',compon) = fc.l('3',compon);
fc.l('5',compon) = fc.l('4',compon);
t.l('3') = t.l('1');
t.l('4') = t.l('1');
t.l('5') = t.l('1');
p.l('3') = p.l('3');
p.l('4') = p.l('3');
p.l('5') = p.l('3')*0.25;
* ---------------------------
* - c o m p r e s s o r 1 -
* ---------------------------
fc.l('6',compon) = fc.l('5',compon);
p.l('6') = p.l('1');
presrat.l('1') = (p.l('6')/p.l('5'))**0.23;
t.l('6') = t.l('5')*presrat.l('1');
* ---------------------
* - 1 - m i x e r 1 -
* ---------------------
fc.l('7',compon) = fc.l('1',compon) + fc.l('6',compon);
t.l('7') = t.l('2')*y.l('2') + t.l('6')*y.l('3');
p.l('7') = p.l('1');
* -----------------------
* - h 2 r e c y c l e -
* -----------------------
fc.l('57','h2') = 2.5*y.l('54');
fc.l('57','ch4') = 0.5*y.l('54');
fc.l('57','ben') = 0.0*y.l('54');
fc.l('57','tol') = 0.0*y.l('54');
fc.l('57','dip') = 0.0*y.l('54');
p.l('57') = p.l('1');
fc.l('60','h2') = 14.0*y.l('59') - 2.0*y.l('54');
fc.l('60','ch4') = 14.0*y.l('59') - 2.0*y.l('54');
fc.l('60','ben') = 0.0 *y.l('59') + 0.0*y.l('54');
fc.l('60','tol') = 0.0 *y.l('59') + 0.0*y.l('54');
fc.l('60','dip') = 0.0 *y.l('59') + 0.0*y.l('54');
p.l('60') = p.l('1');
fc.l('65','h2') = 1.0*y.l('63');
fc.l('65','ch4') = 2.5*y.l('63');
fc.l('65','ben') = 0.0*y.l('63');
fc.l('65','tol') = 0.5*y.l('63');
fc.l('65','dip') = 0.0*y.l('63');
p.l('66') = p.l('1');
fc.l('72',compon) = fc.l('60',compon) + fc.l('57',compon) + fc.l('65',compon);
p.l('72') = p.l('1');
fc.l('59',compon) = fc.l('60',compon);
fc.l('56',compon) = fc.l('57',compon);
fc.l('64',compon) = fc.l('65',compon);
* -------------------------
* - t o l r e c y c l e -
* -------------------------
fc.l('43',compon) = 0.0;
fc.l('43','ben') = 0.1;
fc.l('43','tol') = 2.0;
p.l('43') = p.l('1');
fc.l('42',compon) = fc.l('43',compon);
* -------------------
* - t o l f e e d -
* -------------------
f.l('66') = 2.1;
fc.l('66',compon) = f.l('66')*f66comp(compon);
p.l('66') = p.l('1');
* -----------------
* - m i x e r 1 -
* -----------------
fc.l('8',compon) = fc.l('7',compon) + fc.l('43',compon) + fc.l('66',compon) + fc.l('72',compon);
t.l('8') = 3.5;
p.l('8') = p.l('1');
* -----------------------
* - e x c h a n g e r 1
* -----------------------
fc.l('70',compon) = fc.l('8',compon);
t.l('70') = 8.693;
p.l('70') = p.l('1');
* ------------------
* - f u r n a c e -
* ------------------
fc.l('9',compon) = fc.l('8',compon);
t.l('9') = 8.943;
p.l('9') = p.l('1') - 0.4826;
* -----------------------------
* - 1 - s p l i t t e r 2 -
* -----------------------------
fc.l('10',compon) = fc.l('9',compon)*y.l('10');
fc.l('12',compon) = fc.l('9',compon)*y.l('12');
t.l('10') = t.l('9');
t.l('12') = t.l('9');
p.l('10') = p.l('9');
p.l('12') = p.l('9');
* ---------------------
* - r e a c t o r -
* ---------------------
f.l(str) = sum(compon, fc.l(str,compon));
rctt.l(rct) = 9.2;
krct.l(rct) = 6.3e+10*exp(-26167./(rctt.l(rct)*100.));
rctvol.l('1') = 100.*y.l('10');
rctvol.l('2') = 100.*y.l('12');
conv.l('1','tol') = 1. - sqr(1./(1. + 0.372*krct.l('1')*rctvol.l('1')*sqrt(fc.l('10','tol')/60. + eps1)*(f.l('10')/60. + eps1)**(-3./2.)));
conv.l('2','tol') = 1. - sqr(1./(1. + 0.372*krct.l('2')*rctvol.l('2')*sqrt(fc.l('12','tol')/60. + eps1)*(f.l('12')/60. + eps1)**(-3./2.)));
sel.l(rct) = 1. - 0.0036*(1. - conv.l(rct,'tol'))**(-1.544);
consum.l('1','tol') = conv.l('1','tol')*fc.l('10','tol');
consum.l('2','tol') = conv.l('2','tol')*fc.l('12','tol');
fc.l('11','tol') = fc.l('10','tol') - consum.l('1','tol');
fc.l('11','ben') = fc.l('10','ben') + consum.l('1','tol')* sel.l('1');
fc.l('11','dip') = fc.l('10','dip') + consum.l('1','tol')*(1 - sel.l('1'))/2.;
fc.l('11','h2') = fc.l('10','h2') - consum.l('1','tol')*(1 + sel.l('1'))/2.;
fc.l('11','ch4') = fc.l('10','ch4') + consum.l('1','tol');
fc.l('13','tol') = fc.l('12','tol') - consum.l('2','tol');
fc.l('13','ben') = fc.l('12','ben') + consum.l('2','tol')* sel.l('2');
fc.l('13','dip') = fc.l('12','dip') + consum.l('2','tol')*(1 - sel.l('2'))/2.;
fc.l('13','h2') = fc.l('12','h2') - consum.l('2','tol')*(1 + sel.l('2'))/2.;
fc.l('13','ch4') = fc.l('12','ch4') + consum.l('2','tol');
f.l(str) = sum(compon, fc.l(str,compon));
t.l('11') = (heatrxn('1')*consum.l('1','tol')/(100.*cp('10')) + f.l('10')*t.l('10'))/(f.l('11') + eps1);
t.l('13') = (heatrxn('2')*consum.l('2','tol')/(100.*cp('12')) + f.l('12')*t.l('12'))/(f.l('13') + eps1);
p.l('11') = p.l('9') - 0.2068;
p.l('13') = p.l('9') - 0.2068;
* -------------------------
* - 1 - m i x e r 2 -
* -------------------------
fc.l('14',compon) = fc.l('11',compon) + fc.l('13',compon);
t.l('14') = t.l('11')*y.l('10') + t.l('13')*y.l('12');
p.l('14') = p.l('11');
* ----------------
* - e x p 3 -
* ----------------
fc.l('15',compon) = fc.l('14',compon);
t.l('15') = t.l('14');
p.l('15') = p.l('14');
* --------------------
* - m i x e r 2 -
* --------------------
fc.l('16','h2') = fc.l('15','h2');
fc.l('16','ch4') = fc.l('15','ch4');
fc.l('16','ben') = fc.l('15','ben')*1.2;
fc.l('16','tol') = fc.l('15','tol')*1.2;
fc.l('16','dip') = fc.l('15','dip')*1.2;
t.l('16') = 8.943;
p.l('16') = p.l('15');
* -----------------------
* - e x c h a n g e r 1
* -----------------------
fc.l('71',compon) = fc.l('16',compon);
t.l('71') = 5.0;
p.l('71') = p.l('16');
* ------------------
* - c o o l e r 1 -
* ------------------
fc.l('17',compon) = fc.l('71',compon);
t.l('17') = 3.0;
p.l('17') = p.l('71');
* ---------------------
* - f l a s h 1 -
* ---------------------
eflsh.l('1','h2') = 0.995;
eflsh.l('1','ch4') = 0.99;
eflsh.l('1','ben') = 0.04;
eflsh.l('1','tol') = 0.01;
eflsh.l('1','dip') = 0.0001;
fc.l('18',compon) = fc.l('17',compon)*eflsh.l('1',compon);
fc.l('19',compon) = fc.l('17',compon) - fc.l('18',compon);
t.l('17') = 3.0;
t.l('18') = 3.0;
t.l('19') = 3.0;
p.l('18') = p.l('17');
p.l('19') = p.l('17');
vp.l('19',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.l('19')*100. + antc(compon)));
* -------------------------------------
* - q u e n c h s p l i t t e r -
* -------------------------------------
fc.l('20',compon) = fc.l('19',compon)*0.20;
fc.l('21',compon) = fc.l('19',compon)*0.80;
t.l('20') = t.l('19');
t.l('21') = t.l('19');
p.l('20') = p.l('19');
p.l('21') = p.l('19');
* --------------------
* - m i x e r 3 -
* --------------------
fc.l('69',compon) = 0.0*y.l('63');
fc.l('22',compon) = fc.l('21',compon) + fc.l('69',compon);
t.l('22') = t.l('21');
p.l('22') = p.l('21');
* ---------------------------
* - 1 - s p l i t t e r 3 -
* ---------------------------
fc.l('23',compon) = fc.l('22',compon)*y.l('23');
fc.l('24',compon) = fc.l('22',compon)*y.l('24');
t.l('23') = t.l('22');
t.l('24') = t.l('22');
p.l('23') = p.l('22');
p.l('24') = p.l('22');
* ------------------
* - h e a t e r 1 -
* ------------------
fc.l('25',compon) = fc.l('24',compon);
t.l('25') = t.l('24');
p.l('25') = p.l('24');
* ------------------------
* - c o l u m n - 1 -
* ------------------------
nmin.l('1') = 2.5*y.l('24');
ndist.l('1') = nmin.l('1')*4;
rmin.l('1') = 0.06;
reflux.l('1') = 1.2*rmin.l('1');
distp.l('1') = 1.02;
avevlt.l('1') = 250.;
fc.l('26','h2') = fc.l('25','h2') *1.00;
fc.l('26','ch4') = fc.l('25','ch4')*0.99;
fc.l('26','ben') = fc.l('25','ben')*0.01;
fc.l('26','tol') = fc.l('25','tol')*0.00;
fc.l('26','dip') = fc.l('25','dip')*0.00;
fc.l('27',compon) = fc.l('25',compon) - fc.l('26',compon);
p.l('26') = distp.l('1');
p.l('27') = distp.l('1');
t.l('26') = t.lo('26');
t.l('27') = t.lo('27') + 0.1 ;
vp.l('26','ch4') = (1./7500.6168)*exp(anta('ch4') - antb('ch4')/(t.l('26')*100. + antc('ch4')));
vp.l('26','ben') = (1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('26')*100. + antc('ben')));
vp.l('27','ch4') = (1./7500.6168)*exp(anta('ch4') - antb('ch4')/(t.l('27')*100. + antc('ch4')));
vp.l('27','ben') = (1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('27')*100. + antc('ben')));
* ------------------
* - h e a t e r 1 -
* ------------------
fc.l('44',compon) = fc.l('23',compon);
p.l('44') = p.l('23');
t.l('44') = 3.0*(3.24/0.1)**0.23;
* ----------------
* - e x p 3 -
* ----------------
fc.l('45',compon) = fc.l('44',compon);
p.l('45') = 0.1;
t.l('45') = 3.0;
* ------------------
* - c o o l e r 2 -
* ------------------
fc.l('46',compon) = fc.l('45',compon);
p.l('46') = p.l('45');
t.l('46') = t.l('45');
* ---------------------
* - f l a s h 2 -
* ---------------------
t.l('47') = t.l('46');
t.l('48') = t.l('46');
p.l('47') = p.l('46');
p.l('48') = p.l('46');
vp.l('48',compon) = (1./7500.6168)*exp(anta(compon) - antb(compon)/(t.l('48')*100. + antc(compon)));
eflsh.l('2','h2') = 0.99;
eflsh.l('2','ch4') = 0.99;
eflsh.l('2','ben') = 0.05;
eflsh.l('2','tol') = 0.02;
eflsh.l('2','dip') = 0.0001;
fc.l('47',compon) = fc.l('46',compon)*eflsh.l('2',compon);
fc.l('48',compon) = fc.l('46',compon) - fc.l('47',compon);
* ----------------
* - e x p 4 -
* ----------------
fc.l('49',compon) = fc.l('47',compon);
p.l('49') = p.l('47');
t.l('49') = t.l('47');
* -------------------------
* - 1 - m i x e r 3 -
* -------------------------
fc.l('30',compon) = fc.l('27',compon) + fc.l('48',compon);
t.l('30') = t.l('27')*y.l('24') + t.l('48')*y.l('23');
p.l('30') = p.l('27')*y.l('24') + p.l('48')*y.l('23');
* ------------------------
* - s p l i t t e r 3 -
* ------------------------
e.l('28') = 0.9;
e.l('29') = 1.0 - e.l('28');
fc.l('28',compon) = fc.l('26',compon)*e.l('28');
fc.l('29',compon) = fc.l('26',compon)*e.l('29');
p.l('28') = p.l('26');
p.l('29') = p.l('29');
t.l('28') = t.l('26');
t.l('29') = t.l('29');
* ----------------
* - e x p 5 -
* ----------------
fc.l('50',compon) = fc.l('29',compon);
p.l('50') = p.l('49');
t.l('50') = t.l('29')*(p.l('50')/p.l('29'))**0.23;
* ----------------------
* - m i x e r 5 -
* ----------------------
fc.l('51',compon) = fc.l('49',compon) + fc.l('50',compon);
p.l('51') = p.l('49');
t.l('51') = t.l('49')*y.l('23') + t.l('50')*y.l('24');
* ------------------------
* - c o l u m n - 2 -
* ------------------------
nmin.l('2') = 12.0;
ndist.l('2') = nmin.l('2')*4;
rmin.l('2') = 1.3;
reflux.l('2') = 1.2*rmin.l('2');
distp.l('2') = 0.1;
avevlt.l('2') = 2.5;
p.l('31') = distp.l('2');
p.l('32') = distp.l('2');
t.l('31') = t.lo('31') + 0.1;
t.l('32') = t.lo('32') + 0.1;
fc.l('31','h2') = fc.l('30','h2') *1.00;
fc.l('31','ch4') = fc.l('30','ch4')*1.00;
fc.l('31','ben') = fc.l('30','ben')*0.99;
fc.l('31','tol') = fc.l('30','tol')*0.01;
fc.l('31','dip') = fc.l('30','dip')*0.00;
fc.l('32',compon) = fc.l('30',compon) - fc.l('31',compon);
vp.l('31','ben') = ( 1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('31')*100. + antc('ben')));
vp.l('31','tol') = ( 1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('31')*100. + antc('tol')));
vp.l('32','ben') = ( 1./7500.6168)*exp(anta('ben') - antb('ben')/(t.l('32')*100. + antc('ben')));
vp.l('32','tol') = ( 1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('32')*100. + antc('tol')));
* ----------------------------
* - 1 - s p l i t t e r 4 -
* ----------------------------
fc.l('33',compon) = fc.l('32',compon)*y.l('33');
fc.l('37',compon) = fc.l('32',compon)*y.l('37');
t.l('33') = t.l('32');
t.l('37') = t.l('32');
p.l('33') = p.l('32');
p.l('37') = p.l('32');
* ------------------------
* - c o l u m n - 3 -
* ------------------------
nmin.l('3') = 2.5*y.l('33');
ndist.l('3') = nmin.l('3')*4;
rmin.l('3') = 0.03;
reflux.l('3') = 1.2*rmin.l('3');
distp.l('3') = 0.1;
avevlt.l('3') = 40.;
fc.l('34','h2') = fc.l('33','h2') *1.00;
fc.l('34','ch4') = fc.l('33','ch4')*1.00;
fc.l('34','ben') = fc.l('33','ben')*1.00;
fc.l('34','tol') = fc.l('33','tol')*0.99;
fc.l('34','dip') = fc.l('33','dip')*0.00;
fc.l('35',compon) = fc.l('33',compon) - fc.l('34',compon);
p.l('34') = distp.l('3');
p.l('35') = distp.l('3');
t.l('34') = t.lo('34') + 0.1;
t.l('35') = t.lo('35') + 0.1;
vp.l('34','tol') = (1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('34')*100. + antc('tol')));
vp.l('34','dip') = (1./7500.6168)*exp(anta('dip') - antb('dip')/(t.l('34')*100. + antc('dip')));
vp.l('35','tol') = (1./7500.6168)*exp(anta('tol') - antb('tol')/(t.l('35')*100. + antc('tol')));
vp.l('35','dip') = (1./7500.6168)*exp(anta('dip') - antb('dip')/(t.l('35')*100. + antc('dip')));
* --------------------
* - h e a t e r 3 -
* --------------------
fc.l('38',compon) = fc.l('37',compon);
p.l('38') = p.l('37');
t.l('38') = t.l('37');
* ----------------
* - e x p 2 -
* ----------------
fc.l('39',compon) = fc.l('38',compon);
p.l('39') = 0.1;
t.l('39') = t.l('38');
* ---------------------
* - f l a s h 3 -
* ---------------------
eflsh.l('2','h2') = 0.999;
eflsh.l('2','ch4') = 0.999;
eflsh.l('2','ben') = 0.999;
eflsh.l('2','tol') = 0.95;
eflsh.l('2','dip') = 0.25;
fc.l('40',compon) = fc.l('39',compon) * eflsh.l('2',compon);
fc.l('41',compon) = fc.l('39',compon) - fc.l('40',compon);
t.l('41') = 4.0;
vp.l('41',compon) = (1/7500.6168)*exp(anta(compon) - antb(compon)/(t.l('41')*100. + antc(compon)));
p.l('40') = p.l('39');
p.l('41') = p.l('39');
t.l('40') = t.l('41');
* ------------------------
* - s p l i t t e r 2 -
* ------------------------
e.l('52') = 0.2;
e.l('58') = 1.0 - e.l('52');
fc.l('52',compon) = fc.l('18',compon)*e.l('52');
fc.l('58',compon) = fc.l('18',compon)*e.l('58');
p.l('52') = p.l('18');
p.l('58') = p.l('18');
t.l('52') = t.l('18');
t.l('58') = t.l('18');
* ----------------------------
* - 1 - s p l i t t e r 5 -
* ----------------------------
fc.l('53',compon) = fc.l('52',compon)*y.l('53');
fc.l('54',compon) = fc.l('52',compon)*y.l('54');
t.l('53') = t.l('52');
t.l('54') = t.l('52');
p.l('53') = p.l('52');
p.l('54') = p.l('52');
* ------------------------
* - m e m b r a n e 2 -
* ------------------------
fc.l('55','h2') = fc.l('54','h2') *0.2;
fc.l('55','ch4') = fc.l('54','ch4')*0.9;
fc.l('55','ben') = fc.l('54','ben')*1.0;
fc.l('55','tol') = fc.l('54','tol')*1.0;
fc.l('55','dip') = fc.l('54','dip')*1.0;
fc.l('56',compon) = fc.l('54',compon) - fc.l('55',compon);
t.l('55') = t.l('54');
t.l('56') = t.l('54');
p.l('55') = p.l('54');
p.l('56') = p.l('54')*0.25;
* ----------------------------
* - 1 - s p l i t t e r 6 -
* ----------------------------
fc.l('61',compon) = fc.l('58',compon)*y.l('61');
t.l('61') = t.l('58');
p.l('61') = p.l('58');
* --------------------
* - h e a t e r 4 -
* --------------------
fc.l('73',compon) = fc.l('61',compon);
p.l('73') = p.l('61');
t.l('73') = t.l('61')*(3.24/0.1)**0.23;
* ----------------
* - e x p 6 -
* ----------------
fc.l('62',compon) = fc.l('73',compon);
p.l('62') = p.l('51');
t.l('62') = 3.0;
* -----------------
* - m i x e r 4 -
* -----------------
fc.l('63',compon) = fc.l('51',compon) - fc.l('62',compon);
t.l('63') = t.l('62');
p.l('63') = p.l('62');
* -----------------------
* - a b s o r b e r -
* -----------------------
f.l('63') = sum( compon, fc.l('63',compon));
f.l('67') = f.l('63')*0.20;
fc.l('67',compon) = f.l('67')*f67comp(compon);
nabs.l(abs) = 30.*y.l('63');
gamma.l(abs,'tol') = log((1 - aabs('tol')**(nabs.l(abs)*abseff + eps1))/(1 - aabs('tol')));
beta.l(abs,compon) = log((1 - aabs(compon)**(nabs.l(abs)*abseff + 1.0))/(1 - aabs(compon)));
fc.l('64',compon)$asimp('1',compon) = fc.l('63',compon)/cbeta(compon);
fc.l('64',compon)$anorm('1',compon) = fc.l('63',compon)/exp( beta.l('1',compon));
fc.l('64',compon)$asolv('1',compon) = (fc.l('63',compon) + fc.l('67',compon)*exp(gamma.l('1',compon)))
/exp(beta.l('1',compon));
fc.l('65','ben') = fc.l('64','ben');
fc.l('68',compon) = fc.l('63',compon) + fc.l('67',compon) - fc.l('64',compon);
f.l(str) = sum(compon, fc.l(str,compon));
display f.l, fc.l, t.l, p.l, rctvol.l, conv.l, consum.l, krct.l, sel.l;
display nmin.l, ndist.l, vp.l, gamma.l, beta.l, eflsh.l;
option limCol = 0, limRow = 100, domLim = 1000;
Model douglas / all /;
solve douglas maximizing profit using minlp;
$onText
option limRow = 0;
* y.fx(str) = 0.0;
* y.fx('2') = 1.0;
* y.fx('10') = 1.0;
* y.fx('24') = 1.0;
* y.fx('37') = 1.0;
* y.fx('53') = 1.0;
* y.fx('59') = 1.0;
* switch to reactor #2
y.l('10') = 0.0;
y.l('12') = 1.0;
* -----------------------------
* - 1 - s p l i t t e r 2 -
* -----------------------------
fc.l('10',compon) = fc.l('9',compon)*y.l('10');
fc.l('12',compon) = fc.l('9',compon)*y.l('12');
* ---------------------
* - r e a c t o r -
* ---------------------
f.l(str) = sum(compon, fc.l(str,compon));
krct.l(rct) = 6.3e+10*exp(-26167./(rctt.l(rct)*100.));
rctvol.l('1') = 100.*y.l('10');
rctvol.l('2') = 100.*y.l('12');
conv.l('1','tol') = 1. - sqr(1./(1. + 0.372*krct.l('1')*
rctvol.l('1')*sqrt(fc.l('10','tol')/60. + eps1)*
(f.l('10')/60. + eps1)**(-3./2.)));
conv.l('2','tol') = 1. - sqr(1./(1.+ 0.372*krct.l('2')*
rctvol.l('2')*sqrt(fc.l('12','tol')/60. + eps1)*
(f.l('12')/60. + eps1)**(-3./2.)));
sel.l(rct) = 1. - 0.0036*(1. - conv.l(rct,'tol'))**(-1.544);
consum.l('1','tol') = conv.l('1','tol')*fc.l('10','tol');
consum.l('2','tol') = conv.l('2','tol')*fc.l('12','tol');
fc.l('11','tol') = fc.l('10','tol') - consum.l('1','tol');
fc.l('11','ben') = fc.l('10','ben') + consum.l('1','tol')*sel.l('1');
fc.l('11','dip') = fc.l('10','dip') + consum.l('1','tol')*(1 - sel.l('1'))/2.;
fc.l('11','h2') = fc.l('10','h2') - consum.l('1','tol')*(1 + sel.l('1'))/2.;
fc.l('11','ch4') = fc.l('10','ch4') + consum.l('1','tol');
fc.l('13','tol') = fc.l('12','tol') - consum.l('2','tol');
fc.l('13','ben') = fc.l('12','ben') + consum.l('2','tol')*sel.l('2');
fc.l('13','dip') = fc.l('12','dip') + consum.l('2','tol')*(1 - sel.l('2'))/2.;
fc.l('13','h2') = fc.l('12','h2') - consum.l('2','tol')*(1 + sel.l('2'))/2.;
fc.l('13','ch4') = fc.l('12','ch4') + consum.l('2','tol');
f.l(str) = sum(compon, fc.l(str,compon));
t.l('11') = (heatrxn('1')*consum.l('1','tol')/(100.*cp('10')) + f.l('10')*t.l('10'))/(f.l('11') + eps1);
t.l('13') = (heatrxn('2')*consum.l('2','tol')/(100.*cp('12')) + f.l('12')*t.l('12'))/(f.l('13') + eps1);
solve douglas maximizing profit using nlp;
* add membrane #2
y.l('53') = 0.0;
y.l('54') = 1.0;
* ----------------------------
* - 1 - s p l i t t e r 5 -
* ----------------------------
fc.l('53',compon) = fc.l('52',compon)*y.l('53');
fc.l('54',compon) = fc.l('52',compon)*y.l('54');
* ------------------------
* - m e m b r a n e 2 -
* ------------------------
fc.l('55','h2') = fc.l('54','h2') *0.2;
fc.l('55','ch4') = fc.l('54','ch4')*0.9;
fc.l('55','ben') = fc.l('54','ben')*1.0;
fc.l('55','tol') = fc.l('54','tol')*1.0;
fc.l('55','dip') = fc.l('54','dip')*1.0;
fc.l('56',compon) = fc.l('54',compon) - fc.l('55',compon);
t.l('55') = t.l('54');
t.l('56') = t.l('54');
p.l('55') = p.l('54');
p.l('56') = p.l('54')*0.25;
solve douglas maximizing profit using nlp;
$offText