-
Notifications
You must be signed in to change notification settings - Fork 0
/
MINLP_solver.gms
1331 lines (1171 loc) · 53.6 KB
/
MINLP_solver.gms
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
$ontext
================================================================================
MINLP solver
David Esteban Bernal Neira-David Alejandro Liñan Romero
2019
================================================================================
$offtext
*-------------------------------------------------------------------------------
* Sección 1
* Conjuntos para definir operacion en estado estable
*-------------------------------------------------------------------------------
*Se usa solo el elemento inicial con el punto de colocacion inicial: estado estable
*Esto es equivalente a quitar j y N de la formulacion
sets j "1 punto de colocacion (0)" /1/
N "1 elemento finito" /1/;
*-------------------------------------------------------------------------------
* Sección 2
* Conjuntos, variables, parámetros y ecuaciones principales del sistema
*-------------------------------------------------------------------------------
*Conjuntos
set comp "lista de componentes que intervienen en el sistema" /iButene, Ethanol, nButene, ETBE/;
set Net "Todas las etapas de la columna reales y sobrantes incluyendo con y rev" /1*22/;
*Alias para Net, que se usa en las restricciones binarias:
alias(Net,Net1);
set F "1 representa etanol y 2 representa butenos" /1,2/;
*Variables principales
positive variables
L(N,j,Net) "Flujo de líquido [mol/min]"
V(N,j,Net) "Flujo de vapor [mol/min]"
x(N,j,comp,Net) "Porcentaje molar en el líqudio [%]"
y(N,j,comp,Net) "Porcentaje molar en el vapor [%]"
Temp(N,j,Net) "Temperatura de operación [K]"
P(N,j,Net) "Presión por etapa [bar]"
Z(N,j,Net) "Coeficiente de compresibilidad [-]"
RR(N,j) "Relación molar de reflujo [-]"
Qc(N,j) "Carga térmica del condensador [kJ/min]"
Qr(N,j) "Carga térmica del rehervidor [kJ/min]"
BR(N,j) "Boil up [-]"
;
*Parámetros hidráulicos
parameter
da "Diámetro de los agujeros [m]" /2E-3/
ep "Espesor del plato [m]" /0.002/
pitch "Distancia entre agujeros [m]" /0.009/
Sfactor "Factor de seguridad altura de la columna [-]" /0.15/
poro "Porosidad del plato [-]"
K0 "Coeficiente de orificio [-]"
;
poro=0.907*sqr(da/pitch);
K0=(880.6-(67.7*da/ep)+(7.32*((da/ep)**2))-(0.338*((da/ep)**3)))*1E-3;
*Variables hidraulicas
positive variable D "Diámetro de la columna [m]";
positive variable hw "Weir height [m]";
positive variable HS "Altura de cada plato [m]";
positive variable Htotal "Altura total de la columna [m]";
positive variable At "Area activa [m2]";
positive variable Ad "Area de derramadero [m2]";
positive variable Lw "Weir length [m]";
positive variable A0 "Area agujerada [m2]";
*Ecuaciones hidraulicas
Equation EqHwmin;
EqHwmin.. hw=g=0.05*HS;
Equation EqHwmax;
EqHwmax.. hw=l=HS/3;
equation EqAt;
EqAt.. At=e=sqr(D/2)*(pi-(1.854590-0.96));
equation EqAd;
EqAd.. Ad=e=sqr(D/2)*(0.5*(1.854590-0.96));
equation EqLw;
EqLw.. Lw=e=0.8*D;
equation EqA0;
EqA0.. A0=e=At*poro;
*Alimentacion butenos
parameter FB "Flujo de alimentación de butenos [mol/min]" /5.774/
parameter zb(N,j,comp) "Porcentaje molar en la alimentación de butenos";
zb(N,j,'iButene')=30;
zb(N,j,'nButene')=100-zb(N,j,'iButene');
zb(N,j,'Ethanol')=0;
zb(N,j,'ETBE')=0;
*Alimentacion etanol
parameter
FE "Flujo de alimentación de etanol [mol-h]" /1.7118/
ze(comp) "Porcentaje molar en la alimentación de etanol"
/
iButene 0
Ethanol 100
nButene 0
ETBE 0
/
;
*Parametros de operacion
parameter
Pop "Presión de operación condensador [bar]" /9.5/
TaliB "Temperatura de alimentación de butenos [K]" /323/
TaliE "Temperatura de alimentación etanol [K]" /342.38/
xBetbe "Composición molar de ETBE en fondos deseada" /83/
MCR "Retención constante en rehervidor y condensador [mol]" /1/
cR "Constante de los gases [m3*bar/K*mol]" /0.00008314/
;
*-------------------------------------------------------------------------------
* Sección 3
* Parametro de conversion de unidades
*-------------------------------------------------------------------------------
parameter
hora "Si estamos en análisis por minutos u hora [s]" /60/;
*-------------------------------------------------------------------------------
* Sección 4
* Restricciones de pureza
*-------------------------------------------------------------------------------
equations pureza0(Net);
pureza0(Net)$(ord(Net) eq card(Net)).. x('1','1','ETBE',Net)=g=xBetbe;
*-------------------------------------------------------------------------------
* Sección 5
* Cálculo de presiones de saturación por medio de la ecuación de Antoine
*-------------------------------------------------------------------------------
*Constantes de la ecuación de Antoine expandida
parameters
C1a(comp)
/
iButene 66.4970745
Ethanol 61.7910745
nButene 40.3230745
ETBE 52.67507454
/
C2a(comp)
/
iButene -4634.1
Ethanol -7122.3
nButene -4019.2
ETBE -5820.2
/
C3a(comp)
/
iButene 0
Ethanol 0
nButene 0
ETBE 0
/
C4a(comp)
/
iButene 0
Ethanol 0
nButene 0
ETBE 0
/
C5a(comp)
/
iButene -8.9575
Ethanol -7.1424
nButene -4.5229
ETBE -6.1343
/
C6a(comp)
/
iButene 1.3413E-5
Ethanol 2.8853E-6
nButene 4.8833E-17
ETBE 2.1405E-17
/
C7a(comp)
/
iButene 2
Ethanol 2
nButene 6
ETBE 6
/
;
positive variables Psat(N,j,comp,Net) presión de saturación (bar);
equations EqPsat(N,j,comp,Net);
EqPsat(N,j,comp,Net).. Psat(N,j,comp,Net)=e=exp( C1a(comp) + (C2a(comp)/(Temp(N,j,Net)+C3a(comp))) + (C4a(comp)*Temp(N,j,Net)) + (C5a(comp)*log(Temp(N,j,Net)) + (C6a(comp)*(Temp(N,j,Net)**C7a(comp)))) );
*-------------------------------------------------------------------------------
* Sección 6
* Cálculo de densidades de líquido por medio de la ecuación IK-CAPI
* Cálculo de densidades de líquido por medio de la ecuación DIPPR crítica
* Cálculo de densidades de gas por medio de ecuación de gas ideal corregida
*-------------------------------------------------------------------------------
*Constantes de la ecuación DIPPR
parameters
MW(comp) "Peso molecular [kg/kmol]"
/
iButene 56.10752
Ethanol 46.06904
nButene 56.10752
ETBE 102.17656
/
Tcrit(comp) "Temperatura crítica [K]"
/
iButene 417.9
Ethanol 516.2
nButene 419.6
ETBE 509.4
/
Pcrit(comp) "Presión crítica [bar]"
/
iButene 38.98675
Ethanol 60.35675
nButene 39.18675
ETBE 28.32675
/
C1rh(comp)
/
iButene 8.9711123119
Ethanol -2.932961888E-2
nButene 5.956235579
ETBE -1.323678817E-1
/
C2rh(comp)
/
iButene 0
Ethanol 6.9361857406E-4
nButene 0
ETBE 2.1486345729E-3
/
C3rh(comp)
/
iButene 0
Ethanol -1.962897037E-6
nButene 0
ETBE -6.092181735E-6
/
C4rh(comp)
/
iButene 0
Ethanol 2.089632106E-9
nButene 0
ETBE 6.4627035532E-9
/
C5rh(comp)
/
iButene 0
Ethanol 0
nButene 0
ETBE 0
/
C6rh(comp)
/
iButene -1.4666609E-10
Ethanol 0
nButene -9.3717935E-11
ETBE 0
/
C7rh(comp)
/
iButene 1.286186216E-12
Ethanol 0
nButene 8.150339357E-13
ETBE 0
/
C8rh(comp)
/
iButene -4.33826109E-15
Ethanol 0
nButene -2.72421122E-15
ETBE 0
/
C9rh(comp)
/
iButene 6.619652613E-18
Ethanol 0
nButene 4.115761136E-18
ETBE 0
/
C10rh(comp)
/
iButene -3.8362103001E-21
Ethanol 0
nButene -2.3593237507E-21
ETBE 0
/
C1r(comp)
/
iButene 1.1446
Ethanol 1.6288
nButene 1.0877
ETBE 0.66333
/
C2r(comp)
/
iButene 0.2724
Ethanol 0.27469
nButene 2.6454E-01
ETBE 2.6135E-01
/
C3r(comp)
/
iButene 0.28172
Ethanol 0.23178
nButene 0.2843
ETBE 0.28571
/
C4r(comp)
/
iButene 0
Ethanol 0
nButene 0
ETBE 0
/
;
positive variable Tcritm(N,j,Net);
equation EqTcritm(N,j,Net);
EqTcritm(N,j,Net).. Tcritm(N,j,Net) =e= (sqr(sum(comp,(x(N,j,comp,Net)/100)*Tcrit(comp)/(Pcrit(comp)**0.5))))/(sum(comp,(x(N,j,comp,Net)/100)*Tcrit(comp)/Pcrit(comp)));
positive variables rho(N,j,comp,Net) "Densidad molar por componente de líquido [mol/m^3]";
equation Eqrho(N,j,comp,Net);
Eqrho(N,j,comp,Net).. rho(N,j,comp,Net)=e=( C1r(comp)/(C2r(comp)**(1+((1-(Temp(N,j,Net)/Tcritm(N,j,Net)))**C4r(comp)))) )*1000;
positive variable rhoV(N,j,Net) "Densidad molar de vapor [mol/m^3]";
equation EqurhoV(N,j,Net);
EqurhoV(N,j,Net).. rhoV(N,j,Net)=e=P(N,j,Net)/(0.00008314*Temp(N,j,Net)*(Z(N,j,Net)));
*-------------------------------------------------------------------------------
* Sección 7
* Cálculo de tensión superficial por medio de la ecuación DIPPR crítica
*-------------------------------------------------------------------------------
*Constantes de la ecuación DIPPR
parameters
C1sig(comp)
/
iButene 0.05544
Ethanol 0.03764
nButene 0.055945
ETBE 0.071885
/
C2sig(comp)
/
iButene 1.2453
Ethanol -2.157E-5
nButene 1.2402
ETBE 2.1204
/
C3sig(comp)
/
iButene 0.0
Ethanol 1.025E-7
nButene 0
ETBE -1.5583
/
C4sig(comp)
/
iButene 0
Ethanol 0
nButene 0
ETBE 0.76657
/
;
positive variables sigma(N,j,Net) "Tensión superficial líquido vapor [N/m]";
equation Eqsigma(N,j,Net);
Eqsigma(N,j,Net).. sigma(N,j,Net)=e=sum(comp,(x(N,j,comp,Net)/100)*C1sig(comp)*(1-(Temp(N,j,Net)/Tcritm(N,j,Net)))**(C2sig(comp)+C3sig(comp)*(Temp(N,j,Net)/Tcritm(N,j,Net))+C4sig(comp)*((Temp(N,j,Net)/Tcritm(N,j,Net)))**2));
*-------------------------------------------------------------------------------
* Sección 8
* Cálculo de coeficientes de actividad por medio del modelo NRTL
*-------------------------------------------------------------------------------
table a_nrtl(comp,comp) Parámetro a de NRTL
iButene Ethanol nButene ETBE
iButene 0.0 0.0 0.0 0.0
Ethanol 0.0 0.0 0.0 0.0
nButene 0.0 0.0 0.0 0.0
ETBE 0.0 0.0 0.0 0.0
;
table b_nrtl(comp,comp) Parámetro b de NRTL
iButene Ethanol nButene ETBE
iButene 0.0 623.5810010 107.526499 219.73407
Ethanol 141.9632130 0.0 164.57256 187.104064
nButene -93.24546420 595.5299820 0.0 226.373398
ETBE -172.59152 344.481315 -177.88565 0.0
;
table c_nrtl(comp,comp) Parámetro c de NRTL
iButene Ethanol nButene ETBE
iButene 0.0 0.3 0.3 0.3
Ethanol 0.3 0.0 0.3 0.3
nButene 0.3 0.3 0.0 0.3
ETBE 0.3 0.3 0.3 0.0
;
alias (comp,comp1);
parameter alfa_nrtl(comp,comp);
alfa_nrtl(comp,comp1)$(ord(comp) ne ord(comp1))=c_nrtl(comp,comp1);
*Parámetros G y Tao
variables tao_nrtl(N,j,comp,comp1,Net);
equations Eq_tao_nrtl(N,j,comp,comp1,Net);
Eq_tao_nrtl(N,j,comp,comp1,Net).. tao_nrtl(N,j,comp,comp1,Net)=e=a_nrtl(comp,comp1) + (b_nrtl(comp,comp1)/Temp(N,j,Net));
variables g_nrtl(N,j,comp,comp1,Net);
equations Eq_g_nrtl(N,j,comp,comp1,Net);
Eq_g_nrtl(N,j,comp,comp1,Net).. g_nrtl(N,j,comp,comp1,Net)=e=exp( -alfa_nrtl(comp,comp1)*tao_nrtl(N,j,comp,comp1,Net));
*Coeficiente de actividad (gamma)
alias (comp,comp2,comp3);
variables gamma(N,j,comp,Net);
equations Eqgamma(N,j,comp,Net);
Eqgamma(N,j,comp,Net).. gamma(N,j,comp,Net)=e=
exp(sum(comp1,x(N,j,comp1,Net)*tao_nrtl(N,j,comp1,comp,Net)*
g_nrtl(N,j,comp1,comp,Net))/sum(comp1,x(N,j,comp1,Net)*
g_nrtl(N,j,comp1,comp,Net))+sum(comp1,x(N,j,comp1,Net)*
g_nrtl(N,j,comp,comp1,Net)/sum(comp2,x(N,j,comp2,Net)*
g_nrtl(N,j,comp2,comp1,Net))*(tao_nrtl(N,j,comp,comp1,Net)-
sum(comp2,x(N,j,comp2,Net)*tao_nrtl(N,j,comp2,comp1,Net)*
g_nrtl(N,j,comp2,comp1,Net))/sum(comp3,x(N,j,comp3,Net)*
g_nrtl(N,j,comp3,comp1,Net)))));
*-------------------------------------------------------------------------------
* Sección 9
* Cálculo de reacción química
*-------------------------------------------------------------------------------
Parameter
Nu(comp) "Coeficientes estequiométricos en la reacción"
/
iButene -1
Ethanol -1
nButene 0
ETBE 1
/
mcat "Masa del catalizador" /0.4/
;
variable Ketbe(N,j,Net) "Constante de equilibrio [-]";
$ontext
$offtext
equation EqKetbe(N,j,Net);
EqKetbe(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. Ketbe(N,j,Net) =e=
exp(10.387+4060.59/(Temp(N,j,Net))
-2.89055*log(Temp(N,j,Net))
-0.01915144*Temp(N,j,Net)
+0.0000528586*power(Temp(N,j,Net),2)
-0.0000000532977*power(Temp(N,j,Net),3));
positive variable Krate(N,j,Net) "Tasa de avance de reacción [mol/(kg_cat.min)]";
equation EqKrate(N,j,Net);
EqKrate(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. Krate(N,j,Net) =e= 7.41816E15*exp(-60400.0/(8.314*Temp(N,j,Net)))*hora/3600;
positive variable Ka(N,j,Net) "Tasa de adsorción";
equation EqKa(N,j,Net);
EqKa(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. Ka(N,j,Net) =e= exp(-1.0707+1323.1/Temp(N,j,Net));
variable Rx(N,j,Net) "Tasa de reacción [mol/(kg_cat.min)]";
equation EqRx(N,j,Net);
EqRx(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. Rx(N,j,Net)*(power(1+Ka(N,j,Net)*gamma(N,j,'Ethanol',Net)*x(N,j,'Ethanol',Net)/100,3))*Ketbe(N,j,Net) =e=
(Krate(N,j,Net)*(gamma(N,j,'Ethanol',Net)*x(N,j,'Ethanol',Net)/100))
*((Ketbe(N,j,Net)*gamma(N,j,'iButene',Net)*x(N,j,'iButene',Net)/100*gamma(N,j,'Ethanol',Net)*x(N,j,'Ethanol',Net)/100)
-(gamma(N,j,'ETBE',Net)*x(N,j,'ETBE',Net)/100));
*-------------------------------------------------------------------------------
* Sección 10
* Ecuación de estado (calculo de phi)
*-------------------------------------------------------------------------------
parameter
Omega(comp) "Factor acéntrico [-]"
/
iButene 0.19484
Ethanol 0.643558
nButene 0.184495
ETBE 0.316231
/
TcritSRK(comp) "Temperatura crítica de Soave-Redlich-Kwong [K]"
/
iButene 417.9
Ethanol 514
nButene 419.5
ETBE 509.4
/
mEOS(comp) "Parameter m in EOS"
biEOS(comp) "Parameter bi in EOS"
;
mEOS(comp)=0.48508+1.55171*Omega(comp)-0.15613*sqr(Omega(comp));
biEOS(comp)=0.08664*0.00008314*TcritSRK(comp)/Pcrit(comp);
positive variable alphaEOS(N,j,comp,Net);
equation EqAlphaEOS(N,j,comp,Net);
EqAlphaEOS(N,j,comp,Net).. alphaEOS(N,j,comp,Net) =e= sqr(1+mEOS(comp)*(1-(Temp(N,j,Net)/Tcritm(N,j,Net))**(1/2)));
positive variable aiEOS(N,j,comp,Net);
equation EqaiEOS(N,j,comp,Net);
EqaiEOS(N,j,comp,Net).. aiEOS(N,j,comp,Net) =e= alphaEOS(N,j,comp,Net)*0.42747*(sqr(0.00008314*TcritSRK(comp)))/Pcrit(comp);
positive variable bEOS(N,j,Net);
equation EqbEOS(N,j,Net);
EqbEOS(N,j,Net).. bEOS(N,j,Net) =e= sum(comp,(y(N,j,comp,Net)/100)*biEOS(comp));
positive variable aEOS(N,j,Net);
equation EqaEOS(N,j,Net);
EqaEOS(N,j,Net).. aEOS(N,j,Net) =e= sum(comp,sum(comp1, (y(N,j,comp,Net)/100)*(y(N,j,comp1,Net)/100)*(aiEOS(N,j,comp,Net)*aiEOS(N,j,comp1,Net))**0.5));
equation VaporZ(N,j,Net);
VaporZ(N,j,Net).. (Z(N,j,Net))**3-(Z(N,j,Net))**2+(Z(N,j,Net))
*((aEOS(N,j,Net)*P(N,j,Net)/((0.00008314*Temp(N,j,Net))**2))
-(bEOS(N,j,Net)*P(N,j,Net)/(0.00008314*Temp(N,j,Net)))
-(bEOS(N,j,Net)*P(N,j,Net)/(0.00008314*Temp(N,j,Net)))**2)
-((aEOS(N,j,Net)*P(N,j,Net)/((0.00008314*Temp(N,j,Net))**2)))
*(bEOS(N,j,Net)*P(N,j,Net)/(0.00008314*Temp(N,j,Net))) =e= 0;
positive variable phi(N,j,comp,Net);
equation EqPhi(N,j,comp,Net);
EqPhi(N,j,comp,Net).. phi(N,j,comp,Net) =e= exp(((Z(N,j,Net))-1)*biEOS(comp)/bEOS(N,j,Net)
-log((Z(N,j,Net))-bEOS(N,j,Net))
-(aEOS(N,j,Net)/bEOS(N,j,Net))
*(2*((aiEOS(N,j,comp,Net)/aEOS(N,j,Net))**(1/2))
-biEOS(comp)/bEOS(N,j,Net))*log(((Z(N,j,Net))
-bEOS(N,j,Net))/(Z(N,j,Net))));
*-------------------------------------------------------------------------------
* Sección 11
* Cálculo de entalpías
*-------------------------------------------------------------------------------
*Constantes de Cp (kJ/mol.K) gas ideal
parameters
C1c(comp)
/
iButene 0.016052191
Ethanol 0.00901418
nButene -0.00299356
ETBE -0.014651654
/
C2c(comp)
/
iButene 0.000280432
Ethanol 0.000214071
nButene 0.000353198
ETBE 0.000698631
/
C3c(comp)
/
iButene -0.00000010914988
Ethanol -0.000000083903472
nButene -0.00000019904047
ETBE -0.00000044791741
/
C4c(comp)
/
iButene 0.0000000000090979164
Ethanol 0.0000000000013732704
nButene 0.000000000044631288
ETBE 0.00000000011636811
/
C5c(comp)
/
iButene 0
Ethanol 0
nButene 0
ETBE 0
/
C6c(comp)
/
iButene 0
Ethanol 0
nButene 0
ETBE 0
/
;
parameter
Tref "Temperatura de referencia [K]" /298.15/
Hform(comp) "Entalpía de formación (kJ/mol)"
/
iButene -16.9147
Ethanol -234.963
nButene -0.125604
ETBE -313.9
/
Tb "Temperatura de ebullición de los componentes a P=9.5bar [K]"
/
iButene 341.7
Ethanol 421.9
nButene 342.6
ETBE 438.8
/
;
*Entalpía de la fase vapor (kJ/mol) -- Int(CpdT)
variable HVi(N,j,comp,Net),HV(N,j,Net);
equations EqHVi(N,j,comp,Net),EqHV(N,j,Net);
EqHVi(N,j,comp,Net).. HVi(N,j,comp,Net)=e=( (C1c(comp)*(Temp(N,j,Net)-Tref)) + ((C2c(comp)/2)*((Temp(N,j,Net)**2)-(Tref**2)))
+ ((C3c(comp)/3)*((Temp(N,j,Net)**3)-(Tref**3))) + ((C4c(comp)/4)*((Temp(N,j,Net)**4)-(Tref**4)))
+ ((C5c(comp)/5)*((Temp(N,j,Net)**5)-(Tref**5))) + ((C6c(comp)/6)*((Temp(N,j,Net)**6)-(Tref**6))) + Hform(comp)
+ (8.314/1000)*Temp(N,j,Net)*(Z(N,j,Net)-1)+(1+mEOS(comp))*((aEOS(N,j,Net)**0.5)/bEOS(N,j,Net))*log(Z(N,j,Net)/(Z(N,j,Net)+(bEOS(N,j,Net)*P(N,j,Net)/(0.00008314*Temp(N,j,Net))))));
EqHV(N,j,Net).. HV(N,j,Net)=e=sum(comp,HVi(N,j,comp,Net)*y(N,j,comp,Net)/100);
*Constantes de entalpía de vaporización (kJ/mol)
parameter
C1v(comp)
/iButene 32.614
Ethanol 55.789
nButene 33.774
ETBE 45.29
/
C2v(comp)
/iButene 0.38073
Ethanol 0.31245
nButene 0.5107
ETBE 0.27343
/
C3v(comp)
/iButene 0
Ethanol 0
nButene -0.17304
ETBE 0.21645
/
C4v(comp)
/iButene 0
Ethanol 0
nButene 0.05181
ETBE -0.11756
/
C5v(comp)
/iButene 0
Ethanol 0
nButene 0
ETBE 0
/
;
*Temperaturas reducidas
parameter Tred(comp);
Tred(comp)=Tb(comp)/Tcrit(comp);
parameter alphaEOSb(comp), aiEOSb(comp);
alphaEOSb(comp)=(1+mEOS(comp)*(1-(Tb(comp)/Tcrit(comp))**(1/2)))**2;
aiEOSb(comp)=alphaEOSb(comp)*0.42747*((0.00008314*TcritSRK(comp))**2)/Pcrit(comp);
positive variable Zboil(N,j,comp,Net);
equation VaporZb(N,j,comp,Net);
VaporZb(N,j,comp,Net).. (Zboil(N,j,comp,Net))**3-(Zboil(N,j,comp,Net))**2+(Zboil(N,j,comp,Net))
*((aiEOSb(comp)*P(N,j,Net)/((0.00008314*Tb(comp))**2))
-(biEOS(comp)*P(N,j,Net)/(0.00008314*Tb(comp)))
-(biEOS(comp)*P(N,j,Net)/(0.00008314*Tb(comp)))**2)
-((aiEOSb(comp)*P(N,j,Net)/((0.00008314*Tb(comp))**2)))
*(biEOS(comp)*P(N,j,Net)/(0.00008314*Tb(comp))) =e= 0;
*Entalpía de vaporización (kJ/mol)
parameter DHvap(comp), Hvib(comp);
DHVap(comp)=( C1v(comp)*( (1-Tred(comp))**( C2v(comp) + (C3v(comp)*Tred(comp)) + (C4v(comp)*(Tred(comp)**2)) + (C5v(comp)*(Tred(comp)**3)) ) ) );
HVib(comp)=( (C1c(comp)*(Tb(comp)-Tref)) + ((C2c(comp)/2)*((Tb(comp)**2)-(Tref**2))) + ((C3c(comp)/3)*((Tb(comp)**3)-(Tref**3))) + ((C4c(comp)/4)*((Tb(comp)**4)-(Tref**4))) + ((C5c(comp)/5)*((Tb(comp)**5)-(Tref**5))) + ((C6c(comp)/6)*((Tb(comp)**6)-(Tref**6))) + Hform(comp));
variable depHvib(N,j,comp,Net);
equation EqdepHvib(N,j,comp,Net);
EqdepHvib(N,j,comp,Net).. depHvib(N,j,comp,Net) =e= (8.314/1000)*Tb(comp)*(Zboil(N,j,comp,Net)-1)
+(1+mEOS(comp))*((aiEOSb(comp)**0.5)/biEOS(comp))
*log(Zboil(N,j,comp,Net)/(Zboil(N,j,comp,Net)+(biEOS(comp)*P(N,j,Net)/(0.00008314*Tb(comp)))));
*Constantes de Cp (kJ/mol.K) de líquido
parameter
C1l(comp)
/
iButene 0.08768
Ethanol 0.10264
nButene 0.18205
ETBE 0.11096
/
C2l(comp)
/iButene 0.0002171
Ethanol -0.00013963
nButene -0.001611
ETBE 0.00031422
/
C3l(comp)
/iButene -9.15300E-07
Ethanol -3.03410E-08
nButene 1.19630E-05
ETBE 1.74800E-07
/
C4l(comp)
/iButene 2.2660E-09
Ethanol 2.0386E-09
nButene -3.7454E-08
ETBE 0
/
C5l(comp)
/iButene 0
Ethanol 0
nButene 4.5027E-11
ETBE 0
/
;
*Entalpía de la fase liquida (kJ/mol)
variable HLi(N,j,comp,Net),HL(N,j,Net);
equation EqHLi(N,j,comp,Net),EqHL(N,j,Net);
EqHLi(N,j,comp,Net).. HLi(N,j,comp,Net)=e=HVib(comp)-DHVap(comp)
+((C1l(comp)*(Temp(N,j,Net)-Tb(comp))) + ((C2l(comp)/2)*((Temp(N,j,Net)**2)-(Tb(comp)**2)))
+((C3l(comp)/3)*((Temp(N,j,Net)**3)-(Tb(comp)**3))) + ((C4l(comp)/4)*((Temp(N,j,Net)**4)-(Tb(comp)**4)))
+((C5l(comp)/5)*((Temp(N,j,Net)**5)-(Tb(comp)**5))))+depHvib(N,j,comp,Net);
EqHL(N,j,Net).. HL(N,j,Net)=e=sum(comp,HLi(N,j,comp,Net)*x(N,j,comp,Net)/100);
*-------------------------------------------------------------------------------
* Sección 12
* Cálculo de entalpía de alimentación
*-------------------------------------------------------------------------------
*Entalpía de la alimentación de butenos
parameter HV_b(comp) "Entalpía de vapor de la alimentación [kJ/mol]"
Tred_b(comp) "Temperatura reducida alimentación [-]"
DHVap_b(comp) "Entalpía de vaporización alimentación [kJ/mol]"
HL_b(comp) "Entalpía de líquido de la alimentación [kJ/mol]";
HV_b(comp)=( (C1c(comp)*(TaliB-Tref)) + ((C2c(comp)/2)*((TaliB**2)-(Tref**2))) + ((C3c(comp)/3)*((TaliB**3)-(Tref**3))) + ((C4c(comp)/4)*((TaliB**4)-(Tref**4))) + ((C5c(comp)/5)*((TaliB**5)-(Tref**5))) + ((C6c(comp)/6)*((TaliB**6)-(Tref**6))) + Hform(comp));
Tred_b(comp)=TaliB/Tcrit(comp);
DHVap_b(comp)=( C1v(comp)*( (1-Tred_b(comp))**( C2v(comp) + (C3v(comp)*Tred_b(comp)) + (C4v(comp)*(Tred_b(comp)**2)) + (C5v(comp)*(Tred_b(comp)**3)) ) ) );
HL_b(comp)=HV_b(comp)-DHVap_b(comp);
parameter alphaEOSbut(comp), aiEOSbut(comp), aEOSbut(N,j), bEOSbut(N,j);
alphaEOSbut(comp)=(1+mEOS(comp)*(1-(TaliB/Tcrit(comp))**(1/2)))**2;
aiEOSbut(comp)=alphaEOSbut(comp)*0.42747*((0.00008314*TcritSRK(comp))**2)/Pcrit(comp);
bEOSbut(N,j)=sum(comp,(zb(N,j,comp)/100)*biEOS(comp));
aEOSbut(N,j)=sum(comp,sum(comp1, (zb(N,j,comp)/100)*(zb(N,j,comp1)/100)*(aiEOSbut(comp)*aiEOSbut(comp1))**0.5));
*Zbut se calcula para todas las etapas internas
positive variable Zbut(N,j,Net);
equation VaporZbut(N,j,Net);
VaporZbut(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. (Zbut(N,j,Net))**3-(Zbut(N,j,Net))**2+(Zbut(N,j,Net))
*((aEOSbut(N,j)*P(N,j,Net)/((0.00008314*TaliB)**2))
-(bEOSbut(N,j)*P(N,j,Net)/(0.00008314*TaliB))
-(bEOSbut(N,j)*P(N,j,Net)/(0.00008314*TaliB))**2)
-((aEOSbut(N,j)*P(N,j,Net)/((0.00008314*TaliB)**2)))
*(bEOSbut(N,j)*P(N,j,Net)/(0.00008314*TaliB)) =e= 0;
*HFB se calcula para todas las etapas internas
variable HFB(N,j,Net) "Entalpía de la alimentación de butenos";
equation EqHFB(N,j,Net);
EqHFB(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. HFB(N,j,Net) =e= sum(comp,(zb(N,j,comp)/100)*(HL_b(comp)+(8.314/1000)*TaliB*(Zbut(N,j,Net)-1)
+(1+mEOS(comp))*((aEOSbut(N,j)**0.5)/bEOSbut(N,j))
*log(Zbut(N,j,Net)/(Zbut(N,j,Net)+(bEOSbut(N,j)*P(N,j,Net)/(0.00008314*TaliB))))));
*Entalpía de la alimentación de etanol
parameter HV_e(comp) "Entalpía de vapor de la alimentación [kJ/mol]"
Tred_e(comp) "Temperatura reducida alimentación [K]"
DHVap_e(comp) "Entalpía de vaporización alimentación [kJ/mol]"
HL_e(comp) "Entalpía de líquido de la alimentación [kJ/mol]";
HV_e(comp)=( (C1c(comp)*(TaliE-Tref)) + ((C2c(comp)/2)*((TaliE**2)-(Tref**2))) + ((C3c(comp)/3)*((TaliE**3)-(Tref**3))) + ((C4c(comp)/4)*((TaliE**4)-(Tref**4))) + ((C5c(comp)/5)*((TaliE**5)-(Tref**5))) + ((C6c(comp)/6)*((TaliE**6)-(Tref**6))) + Hform(comp));
Tred_e(comp)=TaliE/Tcrit(comp);
DHVap_e(comp)=( C1v(comp)*( (1-Tred_e(comp))**( C2v(comp) + (C3v(comp)*Tred_e(comp)) + (C4v(comp)*(Tred_e(comp)**2)) + (C5v(comp)*(Tred_e(comp)**3)) ) ) );
HL_e(comp)=HV_e(comp)-DHVap_e(comp);
parameter alphaEOSeth(comp), aiEOSeth(comp), aEOSeth, bEOSeth;
alphaEOSeth(comp)=(1+mEOS(comp)*(1-(TaliE/Tcrit(comp))**(1/2)))**2;
aiEOSeth(comp)=alphaEOSeth(comp)*0.42747*((0.00008314*TcritSRK(comp))**2)/Pcrit(comp);
bEOSeth=sum(comp,(ze(comp)/100)*biEOS(comp));
aEOSeth=sum(comp,sum(comp1, (ze(comp)/100)*(ze(comp1)/100)*(aiEOSeth(comp)*aiEOSeth(comp1))**0.5));
*Zeth se calcula para todas las etapas internas
positive variable Zeth(N,j,Net);
equation VaporZeth(N,j,Net);
VaporZeth(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. (Zeth(N,j,Net))**3-(Zeth(N,j,Net))**2+(Zeth(N,j,Net))
*((aEOSeth*P(N,j,Net)/((0.00008314*TaliE)**2))
-(bEOSeth*P(N,j,Net)/(0.00008314*TaliE))
-(bEOSeth*P(N,j,Net)/(0.00008314*TaliE))**2)
-((aEOSeth*P(N,j,Net)/((0.00008314*TaliE)**2)))
*(bEOSeth*P(N,j,Net)/(0.00008314*TaliE)) =e= 0;
*HFE se alcula para todas las etapas internas
variable HFE(N,j,Net) "Entalpía de la alimentación de etanol [kJ/mol]";
equation EqHFE(N,j,Net);
EqHFE(N,j,Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1)).. HFE(N,j,Net) =e= sum(comp,(ze(comp)/100)*(HL_e(comp)+(8.314/1000)*TaliE*(Zeth(N,j,Net)-1)
+(1+mEOS(comp))*((aEOSeth**0.5)/bEOSeth)
*log(Zeth(N,j,Net)/(Zeth(N,j,Net)+(bEOSeth*P(N,j,Net)/(0.00008314*TaliE))))));
*-------------------------------------------------------------------------------
* Sección 13
* Definicion de parametros, restricciones y variables binarias
*-------------------------------------------------------------------------------
*Parametro que determina si las etapas de rxn deben considerarse como etapas de equilibrio
scalar CASE "0 indica que en las etapas de rxn si hay equilibrio"/0/;
*Existencia de catalizador
binary variable yc(Net) "1 indica que en la etapa si hay catalizador";
*Existencia de reflujo
parameter yr(Net) "1 indica que en la etapa si hay reflujo";
*Existencia de boil up
binary variable yb(Net) "1 indica que en la etapa si hay boil up";
*Permite saber si la etapa es real o sobrante
positive variable par(Net) "1 indica que la etapa es real fisicamente";
equation eqpar1,eqpar2(Net),eqpar3(Net);
eqpar1..par('1')=e=1;
eqpar2(Net)$(ord(Net) eq card(Net))..par(Net)=e=1;
eqpar3(Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1))..par(Net)=e=(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yr(Net1)))+(yb(Net))-(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yb(Net1)));
*Existencia de relaciones de equilibrio
positive variable ye(Net) "1 indica que la etapa es de equilibrio";
equation eqyeq(Net);
eqyeq(Net)$((ord(Net) ne card(Net)) and (ord(Net) ne 1))..ye(Net)=e=par(Net)*(1-yc(Net))*CASE+par(Net)*(1-CASE);
*Existencia de alimentacion
binary variable yf(Net,F) "1 indica que en la etapa hay alimentacion de F";
*-------------------------------------------------------------------------------
* Sección 14
* Restricciones logicas
*-------------------------------------------------------------------------------
scalar cmej /1/;
scalar NCmax "numero maximo de etapas reactivas" /3/;
equation logic1(Net) "The boil up stage is below the reflux stage";
logic1(Net)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yr(Net1)))=g=cmej*(yb(Net));
equation logic2 "There is one reflux stage";
logic2..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le (card(Net1)-1))),yr(Net1)))=e=cmej*1;
equation logic3 "There is one boil up stage";
logic3..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le (card(Net1)-1))),yb(Net1)))=e=cmej*1;
equation logic4(F)"There is one feed stage of EtOH and there is one feed stage of butenes";
logic4(F)..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le (card(Net1)-1))),yf(Net1,F)))=e=cmej*1;
equation logic6 "There is a maximum number of catalytic stages";
logic6.. cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le (card(Net1)-1))),yc(Net1))) =e=cmej*NCmax;
equation logic7(Net,F) "Both feed stages are below the reflux";
logic7(Net,F)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yr(Net1)))=g=cmej*yf(Net,F);
equation logic8(Net,F) "The boil up stage is below the feed stages";
logic8(Net,F)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yf(Net1,F)))=g=cmej*yb(Net);
equation logic9(Net) "The EtOH feed is above the butenes feed";
logic9(Net)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yf(Net1,'1')))=g=cmej*yf(Net,'2');
equation logic10(Net) "The catalytic stages are below the EtOH feed stage";
logic10(Net)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yf(Net1,'1')))=g=cmej*yc(Net);
equation logic11(Net) "The catalytic stages are above the butenes feed stage";
logic11(Net)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*((sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yf(Net1,'2')))-(yf(Net,'2')))=l=cmej*(1-yc(Net));
equation logic12(Net) "The catalytic stages are below the reflux stage";
logic12(Net)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*(sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yr(Net1)))=g=cmej*yc(Net);
equation logic13(Net) "The catalytic stages are above the boil up stage";
logic13(Net)$(ord(Net)>1 and ord(Net)<card(Net))..cmej*((sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yb(Net1)))-(yb(Net)))=l=cmej*(1-yc(Net));
*-------------------------------------------------------------------------------
* Sección 15
* Ecuaciones del condensador
*-------------------------------------------------------------------------------
*Condiciones iniciales (operación en estado estable)
equation BalMasaC0,BalMasaParcialC0(comp),SumaC0,EquilibrioC0(comp),BalEnergiaC0;
BalMasaC0.. 0=e=V('1','1','2')-V('1','1','1')*(1+RR('1','1'));
BalMasaParcialC0(comp).. 0=e=V('1','1','2')*y('1','1',comp,'2')-V('1','1','1')*x('1','1',comp,'1')*(1+RR('1','1'));
SumaC0.. sum(comp,y('1','1',comp,'1')-x('1','1',comp,'1'))=e=0;
EquilibrioC0(comp).. y('1','1',comp,'1')*P('1','1','1')*phi('1','1',comp,'1')=e=Psat('1','1',comp,'1')*gamma('1','1',comp,'1')*x('1','1',comp,'1');
BalEnergiaC0.. 0=e=V('1','1','2')*HV('1','1','2')-V('1','1','1')*(1+RR('1','1'))*HL('1','1','1')-QC('1','1');
*Flujo de liquido fijo
equation fixedL(N,j);
fixedL(N,j)..L(N,j,'1')=e=0;
*-------------------------------------------------------------------------------
* Sección 16
* Ecuaciones de la columna (Punto Inicial)
*-------------------------------------------------------------------------------
*Condiciones iniciales (operación en estado estable)
equations BalMasa0(Net,Net1),BalMasaParcial0(comp,Net,Net1),Suma0(Net),BalEnergia0(Net,Net1);
BalMasa0(Net,Net1)$((ord(Net)>1 and ord(Net)<card(Net)) and (ord(Net1) eq card(Net1)))..0=e=yf(Net,'1')*FE+yf(Net,'2')*FB+RR('1','1')*V('1','1','1')*yr(Net)+BR('1','1')*L('1','1',Net1)*yb(Net)+L('1','1',Net-1)+V('1','1',Net+1)-L('1','1',Net)-V('1','1',Net)+yc(Net)*(sum(comp,Nu(comp))*mcat*Rx('1','1',Net)) ;
BalMasaParcial0(comp,Net,Net1)$((ord(Net)>1 and ord(Net)<card(Net)) and (ord(Net1) eq card(Net1)))..0=e=yf(Net,'1')*FE*ze(comp)+yf(Net,'2')*FB*zb('1','1',comp)+RR('1','1')*V('1','1','1')*yr(Net)*x('1','1',comp,'1')+BR('1','1')*L('1','1',Net1)*yb(Net)*y('1','1',comp,Net1)+L('1','1',Net-1)*x('1','1',comp,Net-1)+V('1','1',Net+1)*y('1','1',comp,Net+1)-L('1','1',Net)*x('1','1',comp,Net)-V('1','1',Net)*y('1','1',comp,Net)+100*yc(Net)*(Nu(comp)*mcat*Rx('1','1',Net));
Suma0(Net)$(ord(Net)>1 and ord(Net)<card(Net)).. sum(comp,x('1','1',comp,Net)-y('1','1',comp,Net))=e=0;
BalEnergia0(Net,Net1)$(ord(Net)>1 and ord(Net)<card(Net) and ord(Net1) eq card(Net1))..0=e=yf(Net,'1')*FE*HFE('1','1',Net)+yf(Net,'2')*FB*HFB('1','1',Net)+RR('1','1')*V('1','1','1')*yr(Net)*HL('1','1','1')+BR('1','1')*L('1','1',Net1)*yb(Net)*HV('1','1',Net1)+L('1','1',Net-1)*HL('1','1',Net-1)+V('1','1',Net+1)*HV('1','1',Net+1)-L('1','1',Net)*HL('1','1',Net)-V('1','1',Net)*HV('1','1',Net);
*Relaciones de equilibrio
equations Equilibrio10(comp,Net);
Equilibrio10(comp,Net)$(ord(Net)>1 and ord(Net)<card(Net))..0=e=ye(Net)*((y('1','1',comp,Net)*P('1','1',Net)*phi('1','1',comp,Net))-(Psat('1','1',comp,Net)*gamma('1','1',comp,Net)*x('1','1',comp,Net)));
equations Equilibrio20(comp,Net);
Equilibrio20(comp,Net)$(ord(Net)>1 and ord(Net)<card(Net) and ord(comp) ne 1)..0=e=(sum(Net1$(ord(Net1) ge 2 and ord(Net1) le ord(Net)),yr(Net1)))*(1-ye(Net))*(y('1','1',comp,Net)-y('1','1',comp,Net+1));
equation Equilibrio30(comp,Net);
Equilibrio30(comp,Net)$(ord(Net)>1 and ord(Net)<card(Net) and ord(comp) ne 1)..0=e=(1-sum(Net1$((ord(Net1) ge 2) and (ord(Net1) le ord(Net))),yr(Net1)))*(1-ye(Net))*(x('1','1',comp,Net)-x('1','1',comp,Net-1));
Equation Equilibrio40(Net,Net1);
Equilibrio40(Net,Net1)$(ord(Net)>1 and ord(Net)<card(Net) and ord(Net1) eq card(Net1))..0=e=(1-ye(Net))*(V('1','1',Net)-V('1','1',Net+1)-BR('1','1')*L('1','1',Net1)*yb(Net));
*-------------------------------------------------------------------------------
* Sección 17
* Ecuaciones del rehervidor
*-------------------------------------------------------------------------------
*Condiciones iniciales (operación en estado estable)
equation BalMasaR0(Net),BalMasaParcialR0(comp,Net),SumaR0(Net),EquilibrioR0(comp,Net),BalEnergiaR0(Net);
BalMasaR0(Net)$(ord(Net) eq card(Net))..0=e=L('1','1',Net-1)-L('1','1',Net)*(1+BR('1','1'));
BalMasaParcialR0(comp,Net)$(ord(Net) eq card(Net))..0=e=L('1','1',Net-1)*x('1','1',comp,Net-1)-L('1','1',Net)*(x('1','1',comp,Net)+BR('1','1')*y('1','1',comp,Net));
SumaR0(Net)$(ord(Net) eq card(Net)).. sum(comp,y('1','1',comp,Net)-x('1','1',comp,Net))=e=0;
EquilibrioR0(comp,Net)$(ord(Net) eq card(Net))..y('1','1',comp,Net)*P('1','1',Net)*phi('1','1',comp,Net)=e=Psat('1','1',comp,Net)*gamma('1','1',comp,Net)*x('1','1',comp,Net);
BalEnergiaR0(Net)$(ord(Net) eq card(Net))..0=e=QR('1','1')+L('1','1',Net-1)*HL('1','1',Net-1)-L('1','1',Net)*HL('1','1',Net)-BR('1','1')*L('1','1',Net)*HV('1','1',Net);
*Variable fija del flujo de vapor en la ultima etapa
equation fixedV(N,j,Net);
fixedV(N,j,Net)$(ord(Net) eq card(Net))..V(N,j,Net)=e=0;
*-------------------------------------------------------------------------------
* Sección 18
* Relaciones hidráulicas para todas las etapas internas
*-------------------------------------------------------------------------------
*Caracteristicas del catalizador
scalar fracvol /0.3/;
scalar fracEnvelop /0.5/;
*Definición de velocidad de vapor
positive variables far(N,j,Net) "Factor de areación [-]";
equations Eqfa(N,j,Net);
Eqfa(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. par(Net)*(far(N,j,Net))=e=par(Net)*(0.981*exp(-0.411*((V(N,j,Net)/(rhoV(N,j,Net))/hora)*(rhoV(N,j,Net)*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000)**(0.5))/At));
positive variable hD(N,j,Net) "Altura del líquido por encima del divisor [m]";
equations EqhD(N,j,Net);
EqhD(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. (hD(N,j,Net))=e=(0.6*(((((L(N,j,Net)/sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100))/hora)/Lw))**(2/3)));
positive variable uhv(N,j,Net) "Velocidad del vapor por los agujeros [m/s]";
equations Equhv(N,j,Net);
Equhv(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. par(Net)*(uhv(N,j,Net))=e=par(Net)*((V(N,j,Net)/(rhoV(N,j,Net))/hora)/A0);
positive variable unv(N,j,Net) "Velocidad del vapor por el plato [m/s]";
equations Equnv(N,j,Net);
Equnv(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. par(Net)*unv(N,j,Net)=e=par(Net)*((V(N,j,Net)/(rhoV(N,j,Net))/hora)/At);
*Definicion de velocidad del liquido
positive variable ul(N,j,Net) "Velocidad del líquido en el derramadero [m/s]";
equations Equl(N,j,Net);
Equl(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. par(Net)*ul(N,j,Net)=e=par(Net)*((L(N,j,Net)/(sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100))/hora)/Ad);
*Carga de liquido
positive variable hcl(N,j,Net) "Altura del líquido libre en régimen de spray [m]"
equation Eqhcl(N,j,Net);
scalar consmach /1e-20/;
Eqhcl(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. par(Net)*hcl(N,j,Net)=e=par(Net)*((0.157*(poro**(-0.791))/(1+1.04E-4*(((((L(N,j,Net)+consmach)/sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100))/hora)/Lw)**(-0.59))
*(poro**(-1.791))))*(da**0.833)
*(996/(sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000))**(0.5*(1-0.91*da/poro)));
positive variable Csbf(N,j,Net);
equation EqCsbf(N,j,Net);
EqCsbf(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. par(Net)*(Csbf(N,j,Net))=e=par(Net)*(0.37*(((sqr(da)*sigma(N,j,Net)/(sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000)))**0.125)
*((((rhoV(N,j,Net))*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000)/(sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000))**0.1)
*((HS/hcl(N,j,Net))**0.5));
*Carga de liquido en etapas cataliticas
positive variable Lload(N,j,Net) "carga de liquido en etapas cataliticas [m_s]";
equation eqLload(N,j,Net);
eqLload(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net))..(1-fracvol)*((3.1415926/4)*(D**2))*Lload(N,j,Net)=e=ul(N,j,net)*Ad ;
*Factor de flujo de vapor en etapas cataliticas
positive variable Ffactor(N,j,Net) "Factor de flujo de vapor para etapas cataliticas [Pa**0.5]"
equation eqFfactor(N,j,Net);
eqFfactor(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. (1-fracvol)*(3.1415926/4)*(D**2)*(((rhov(N,j,net))*(sum(comp,(y(N,j,comp,Net)/100)*MW(comp)))*(1/1000))**(1/2))*(Ffactor(N,j,Net))=e=(V(N,j,net)*(1/60))*(sum(comp,(y(N,j,comp,Net)/100)*MW(comp)*(1/1000)));
*Caida de presion
positive variables DPL(N,j,Net) "Caída de presión por la presencia de líquido [bar]";
equations EqDPL(N,j,Net);
EqDPL(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. (DPL(N,j,Net))=e=((far(N,j,Net)*9.81*(sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000)*(hD(N,j,Net)+hw))/100000);
positive variables DPS(N,j,Net) "Caída de presión debido a la presencia de los agujeros - seco [bar]";
equations EqDPS(N,j,Net);
EqDPS(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. (DPS(N,j,Net))=e=((1/(2*sqr(K0)))*( (((sqr(V(N,j,Net)/(rhoV(N,j,Net))/hora)/A0)) )*((rhoV(N,j,Net))*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000)*(1-sqr(poro)))/100000);
positive variable DPq(N,j,Net) "Caída de presión en el derramadero [bar]";
equations EqDPq(N,j,Net);
EqDPq(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. DPq(N,j,Net)=e=(1/(100000))*1.62*((sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000))/(sqr(Lw*hw))*(sqr((L(N,j,Net)/sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100))/hora)+sqr((V(N,j,Net)/(rhoV(N,j,Net))/hora)));
positive variables DP(N,j,Net) "Caída de presión total [bar]";
positive variable dPcat(N,j,net) "caida de presion por catalizador en etapas cataliticas [bar]";
equations EqDP(N,j,Net),EqDPR(N,j,Net),EqdPcat(N,j,net),EqP(N,j,Net),EqPC(N,j,Net),EqPR(N,j,Net) "Definición de presión por etapa [bar]";
EqDPR(N,j,Net)$(ord(Net) eq card(Net)).. DP(N,j,Net)=e=DP(N,j,Net-1);
EqdPcat(N,j,net)$(ord(Net)>1 and ord(Net)<card(Net))..dPcat(N,j,net)=e=hs*fracEnvelop*(0.001)*( (5.69228924748553E-06)*((Lload(N,j,Net)*60*60)**3.05308055949085)*((Ffactor(N,j,Net))**7.851695947) + 1.367015225*((Ffactor(N,j,Net))**1.764157687) );
EqDP(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. (DP(N,j,Net))=e=(DPS(N,j,Net)+DPL(N,j,Net))+yc(Net)*dPcat(N,j,net);
EqP(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. P(N,j,Net)=e=P(N,j,Net-1)+par(Net)*DP(N,j,Net);
EqPC(N,j,Net)$(ord(Net) eq 1).. P(N,j,Net)=e=Pop;
EqPR(N,j,Net)$(ord(Net) eq card(Net)).. P(N,j,Net)=e=P(N,j,Net-1);
*Efectos indeseados en la columna
*Downflow flooding (inundación en los derramaderos)
equation DownFlood(N,j,Net);
DownFlood(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net))..0=g=((HD(N,j,Net)+((DP(N,j,Net)+DPq(N,j,Net))*100000)
/(9.81*(((sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000))
-(rhoV(N,j,Net)*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000))))-(HS))*par(Net);
*Entrainment flooding (inundación por arrastre de líquido)
equation EntrainFloodV(N,j,Net), EntrainFloodL(N,j,Net);
EntrainFloodV(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net))..par(Net)*((unv(N,j,Net))-(Csbf(N,j,Net)*(((((sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000))
-(rhoV(N,j,Net)*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000)))
/(rhoV(N,j,Net)*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000))**0.5))=l=0;
EntrainFloodL(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net))..par(Net)*((ul(N,j,Net))-((sigma(N,j,Net)*9.81*(((sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000))
-(rhoV(N,j,Net)*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000))
/((sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000)**2))**(1/4)))=l=0
*Weeping (lloriqueo)
equation Weep(N,j,Net);
Weep(N,j,Net)$(ord(Net)>1 and ord(Net)<card(Net)).. 0=g=(((0.68-0.12)/(((rhoV(N,j,Net)*sum(comp,MW(comp)*y(N,j,comp,Net)/100)/1000)
/((sum(comp,rho(N,j,comp,Net)*x(N,j,comp,Net)/100)*sum(comp,MW(comp)*x(N,j,comp,Net)/100)/1000)
*9.81*far(N,j,Net)*(hw+hd(N,j,Net))))**0.5))-(uhv(N,j,Net)))*par(Net);
*Catalyst flooding (inundación del empaque del catalizador)
equation catflood(N,j,net);
catflood(N,j,net)..yc(net)*(dPcat(N,j,net)-(12E-3)*hs*fracEnvelop)=l=0;
*Construcción de la columna
equation Size "Tamaño del equipo";
Size.. 1*Htotal =e= 1*((1+Sfactor)*sum(Net$(ord(Net)>1 and ord(Net)<card(Net)),HS*par(Net)));
equation ammountcat "Espacio disponible para el catalizador";
ammountcat..mcat=l=(fracvol)*((3.1415926/4)*(D**2))*(hs*fracEnvelop)*770;
equation DtoLratio "Relacion entre el diametro y la altura";
DtoLratio..htotal/D=l=20;
*-------------------------------------------------------------------------------