source: trunk/sample/miscellaneous/Discrete/propane_propene.mso @ 978

Last change on this file since 978 was 978, checked in by Argimiro Resende Secchi, 6 years ago

examples of discrete polynomial approximation
File size: 31.0 KB
Line 
1#*---------------------------------------------------------------------
2* EMSO Model Library (EML) Copyright (C) 2004 - 2016 ALSOC.
3*
4* This LIBRARY is free software; you can distribute it and/or modify
5* it under the therms of the ALSOC FREE LICENSE as available at
6* http://www.enq.ufrgs.br/alsoc.
7*
8* EMSO Copyright (C) 2004 - 2016 ALSOC, original code
9* from http://www.rps.eng.br Copyright (C) 2002-2004.
10* All rights reserved.
11*
12* EMSO is distributed under the therms of the ALSOC LICENSE as
13* available at http://www.enq.ufrgs.br/alsoc.
14*
15*----------------------------------------------------------------------
16* Author: Argimiro R. Secchi
17* COPPE/UFRJ - Group of Modeling, Simulation, Control,
18*                               and Optimization of Processes
19* $Id$
20*--------------------------------------------------------------------*#
21
22using "types";
23
24Model Coluna_completa_fix2
25       
26        PARAMETERS
27outer   PP      as Plugin(Brief="Physical Properties",Type="PP");
28outer   NComp   as Integer;
29        NTrays  as Integer(Brief="Number of trays");
30        N1              as Integer;
31        N2              as Integer;
32        Nf              as Integer;
33        Nfm1    as Integer;
34        Nfp1    as Integer;
35        M(N2)   as mol;
36
37        SET
38        N1 = NTrays + 1;
39        N2 = NTrays + 2;
40        Nfm1 = Nf - 1;
41        Nfp1 = Nf + 1;
42        M(1) = 20 * 'kmol';
43        M(2:N1) = 0.4 * 'kmol';
44        M(N2) = 25 * 'kmol';
45
46        VARIABLES
47        F as flow_mol;
48        D as flow_mol;
49        L(N2) as flow_mol;
50        V(N2) as flow_mol;
51        z(NComp) as fraction;
52        x(NComp,N2) as fraction;
53        y(NComp,N2) as fraction;
54        m(NComp,N2) as mol;
55        K(NComp,N2) as positive;
56#       H(N2) as enth_mol;
57#       h(N2) as enth_mol;
58#       hf        as enth_mol;
59        T(N2) as temperature;
60        Tn(N2) as temperature;
61        Tf        as temperature;
62        P(N2) as pressure;
63        Pf        as pressure;
64#       Qc        as heat_rate;
65#       Qr        as heat_rate;
66#       vf        as fraction;
67#       xf(NComp) as fraction;
68#       yf(NComp) as fraction;
69        RR        as positive;
70
71        EQUATIONS
72
73        Tn = (T-318*'K')*100;
74
75    "Condenser mass balance"
76    V(2) = L(1) + D + V(1);
77
78    "Condenser component mass balance"
79    diff(m(:,1)) = V(2)*y(:,2) - (L(1) + D)*x(:,1) - V(1)*y(:,1);
80
81#       "Reboiler energy balance"
82#   Qc = V(2)*H(2) - (L(1) + D)*h(1) - V(1)*H(1);
83
84        for j in [2:Nfm1,Nfp1:N1] do
85           "Tray mass balance"
86#           L(j) = L(j-1) + V(j+1) - V(j);
87            L(j) = L(j-1);
88
89                "Tray component mass balance"
90           diff(m(:,j)) = L(j-1)*x(:,j-1) + V(j+1)*y(:,j+1) - L(j)*x(:,j) - V(j)*y(:,j);
91       
92#               "Tray energy balance"
93#           V(j)*H(j) = L(j-1)*h(j-1) + V(j+1)*H(j+1) - L(j)*h(j);
94                V(j+1) = V(j);
95        end
96
97        "Feed tray mass balance"
98#    L(Nf) = F + L(Nfm1) + V(Nfp1) - V(Nf);
99    L(Nf) = F + L(Nfm1);
100
101        "Feed tray component mass balance"
102    diff(m(:,Nf)) = F*z + L(Nfm1)*x(:,Nfm1) + V(Nfp1)*y(:,Nfp1) - L(Nf)*x(:,Nf) - V(Nf)*y(:,Nf);
103
104#       "Feed tray energy balance"
105#       V(Nf)*H(Nf) = F*hf + L(Nfm1)*h(Nfm1) + V(Nfp1)*H(Nfp1) - L(Nf)*h(Nf);
106    V(Nfp1) = V(Nf);
107
108    "Reboiler mass balance"
109#   L(N2) = L(N1) - V(N2);
110    L(N2) = F - V(1) - D;
111
112    "Reboiler component mass balance"
113    diff(m(:,N2)) = L(N1)*x(:,N1) - L(N2)*x(:,N2) - V(N2)*y(:,N2);
114       
115#       "Reboiler energy balance"
116#   -Qr = L(N1)*h(N1) - L(N2)*h(N2) - V(N2)*H(N2);
117
118        for j in [1:N2] do
119                "Mol fraction"
120                m(:,j) = M(j) * x(:,j);
121               
122#               "Liquid Enthalpy"
123#               h(j) = PP.LiquidEnthalpy(T(j), P(j), x(:,j));
124
125#               "Vapour Enthalpy"
126#               H(j) = PP.VapourEnthalpy(T(j), P(j), y(:,j));
127               
128                "Chemical Equilibrium"
129                PP.LiquidFugacityCoefficient(T(j), P(j), x(:,j)) =
130                PP.VapourFugacityCoefficient(T(j), P(j), y(:,j))*K(:,j);
131               
132                y(:,j) = K(:,j)*x(:,j);
133
134        #       "Liquid mol fraction normalization"
135        #       sum(x(:,j)) = 1;
136
137        #       "Vapour mol fraction normalization"
138        #       sum(x(:,j)) = sum(y(:,j));
139                T(j) = PP.BubbleT(P(j),x(:,j));
140        end
141
142#       "Feed flash Calculation"
143#       [vf, xf, yf] = PP.FlashPH(Pf, hf, z);
144
145#       "Feed enthalpy"
146#       hf = (1-vf)*PP.LiquidEnthalpy(Tf, Pf, xf) + vf*PP.VapourEnthalpy(Tf, Pf, yf);
147
148        "Reflux ratio"
149        RR * (D + V(1)) = L(1);
150end
151
152Model Coluna_reduzida_fix2
153       
154        PARAMETERS
155outer   PP      as Plugin(Brief="Physical Properties",Type="PP");
156outer   Disc_top as Plugin(Type="Discrete");
157outer   Disc_bot as Plugin(Type="Discrete");
158outer   NComp   as Integer;
159        NTrays  as Integer(Brief="Number of trays");
160        N1              as Integer;
161        N2              as Integer;
162        Nf              as Integer;
163        Nfm1    as Integer;
164        Nfp1    as Integer;
165        np_top as Integer;
166        n_top as Integer;
167        n1_top as Integer;
168        n2_top as Integer;
169        Ns_top as Integer;
170        Ns1_top as Integer;
171        M_top(np_top)   as mol;
172        Am_top(Ns_top) as Real;
173        An_top(Ns_top) as Real;
174        Bm_top(Ns1_top) as Real;       
175        Bn_top(Ns1_top) as Real;
176        mV_top(n2_top) as Real(Brief="V matrix");
177        np_bot as Integer;
178        n_bot as Integer;
179        n1_bot as Integer;
180        n2_bot as Integer;
181        Ns_bot as Integer;
182        Ns1_bot as Integer;
183        M_bot(np_bot)   as mol;
184        Am_bot(Ns_bot) as Real;
185        An_bot(Ns_bot) as Real;
186        Bm_bot(Ns1_bot) as Real;       
187        Bn_bot(Ns1_bot) as Real;
188        mV_bot(n2_bot) as Real(Brief="V matrix");
189
190        r_top(np_top) as Real(Brief="raizes de hahn");
191        rs_top(np_top) as Real(Brief="raizes de hahn desnormal.");
192        w_top(np_top) as Real(Brief="quadrature weights");
193        Norm_top  as Real(Brief="scaling");
194        r_bot(np_bot) as Real(Brief="raizes de hahn");
195        rs_bot(np_bot) as Real(Brief="raizes de hahn desnormal.");
196        w_bot(np_bot) as Real(Brief="quadrature weights");
197        Norm_bot as Real(Brief="scaling");
198       
199        SET
200        NTrays = Disc_bot.UpperBound-Disc_top.LowerBound-1;
201        N2 = NTrays + 2;
202        Nf = Disc_top.UpperBound;
203        np_top = Disc_top.NodalPoints;
204        n_top = Disc_top.InternalPoints;
205        n1_top = n_top + 1;
206        n2_top = 2*n_top;
207        Ns_top= np_top * np_top;
208        Ns1_top= n_top * np_top;
209        Am_top = Disc_top.Aplus;       
210        An_top = Disc_top.Aminus;
211        Bm_top = Disc_top.Bplus;       
212        Bn_top = Disc_top.Bminus;
213        mV_top = Disc_top.Vmatrix;
214        np_bot = Disc_bot.NodalPoints;
215        n_bot = Disc_bot.InternalPoints;
216        n1_bot = n_bot + 1;
217        n2_bot = 2*n_bot;
218        Ns_bot= np_bot * np_bot;
219        Ns1_bot= n_bot * np_bot;
220        Am_bot = Disc_bot.Aplus;       
221        An_bot = Disc_bot.Aminus;
222        Bm_bot = Disc_bot.Bplus;       
223        Bn_bot = Disc_bot.Bminus;
224        mV_bot = Disc_bot.Vmatrix;
225
226        r_top = Disc_top.Roots;
227        w_top = Disc_top.Weights;
228        Norm_top = Disc_top.Scaling;
229        rs_top = r_top * Norm_top;
230        r_bot = Disc_bot.Roots;
231        w_bot = Disc_bot.Weights;
232        Norm_bot = Disc_bot.Scaling;
233        rs_bot = r_bot * Norm_bot;
234
235        M_top(1) = 20 * 'kmol';
236        M_top(2:np_top) = 0.4 * 'kmol';
237        M_bot(1:n1_bot) = 0.4 * 'kmol';
238        M_bot(np_bot) = 25 * 'kmol';
239       
240        ##############################
241        N1 = NTrays + 1;
242        N2 = NTrays + 2;
243        Nfm1 = Nf - 1;
244        Nfp1 = Nf + 1;
245    ###############################
246       
247        VARIABLES
248        F as flow_mol;
249        D as flow_mol;
250        B as flow_mol;
251        V0 as flow_mol;
252        z(NComp) as fraction;
253        Tf as temperature;
254        Pf as pressure;
255        L(N2) as flow_mol;
256        V(N2) as flow_mol;
257
258        L_top(np_top) as flow_mol;
259        V_top(np_top) as flow_mol;
260        T_top(np_top) as temperature;
261        P_top(np_top) as pressure;
262        x_top(NComp,np_top) as fraction(Lower=-1,Upper=2);
263        y_top(NComp,np_top) as fraction(Lower=-1,Upper=2);
264        K_top(NComp,np_top) as positive;
265        m_top(NComp,np_top) as mol(Lower=-1e3);
266        mb_top(NComp,np_top) as mol(Lower=-1e6);
267        RES(NComp,N2)    as flow_mol(Lower=-1e6);
268       
269        L_bot(np_bot) as flow_mol;
270        V_bot(np_bot) as flow_mol;
271        T_bot(np_bot) as temperature;
272        P_bot(np_bot) as pressure;
273        x_bot(NComp,np_bot) as fraction(Lower=-1,Upper=2);
274        y_bot(NComp,np_bot) as fraction(Lower=-1,Upper=2);
275        K_bot(NComp,np_bot) as positive;
276        m_bot(NComp,np_bot) as mol(Lower=-1e3);
277        mb_bot(NComp,np_bot) as mol(Lower=-1e6);
278        #RES_bot(NComp,np_bot)    as flow_mol;
279
280        x(NComp,N2) as fraction(Lower=-1,Upper=2);
281        y(NComp,N2) as fraction(Lower=-1,Upper=2);
282        T(N2) as temperature;
283        Tn(N2) as temperature;
284
285#       H(N2) as enth_mol;
286#       h(N2) as enth_mol;
287#       hf        as enth_mol;
288#       Qc        as heat_rate;
289#       Qr        as heat_rate;
290#       vf        as fraction;
291#       xf(NComp) as fraction;
292#       yf(NComp) as fraction;
293        RR        as positive;
294
295        EQUATIONS
296
297        Tn = (T-318*'K')*100;
298
299        "Reflux ratio"
300        RR * (D + V0) = L_top(1);
301
302    "Condenser mass balance"
303    V_top(1) = L_top(1) + D + V0;
304    V_top(2) = V_top(1);
305
306        for i in [1:NComp-1] do
307                "Condenser component mass balance"
308                diff(mb_top(i,1)) = sum(Am_top([1:np_top])*V_top*y_top(i,:)) - (L_top(1) + D)*x_top(i,1) - V0*y_top(i,1);
309        #       RES_top= sum(Am_top([1:np_top])*V_top*y_top(i,:)) - (L_top(1) + D)*x_top(i,1) - V0*y_top(i,1);#diff(mb_top(i,1));
310        end
311
312        mb_top(:,1) = m_top(:,1);
313       
314#       "Reboiler energy balance"
315#   Qc = V(2)*H(2) - (L(1) + D)*h(1) - V(1)*H(1);
316
317        for j in [2:n1_top] do
318           "Tray mass balance"
319#           L(j) = L(j-1) + V(j+1) - V(j);
320            L_top(j) = L_top(j-1);
321
322                for i in [1:NComp-1] do
323                        "Tray component mass balance"
324                        diff(mb_top(i,j)) = sum(Bn_top([1:np_top]+(j-2)*np_top)*L_top*x_top(i,:)) +
325                                                   sum(Bm_top([1:np_top]+(j-2)*np_top)*V_top*y_top(i,:)) -
326                                                   L_top(j)*x_top(i,j) - V_top(j)*y_top(i,j) -
327                                                   mV_top(j-1)*(L_top(j)*x_top(i,1) + V_top(j)*y_top(i,1)) -
328                                                   mV_top(j+n_top-1)*(L_top(j)*x_top(i,np_top) + V_top(j)*y_top(i,np_top));
329                    #RES_top(i,j)= diff(mb_top(i,j));
330                end
331               
332                mb_top(:,j) = m_top(:,j) + mV_top(j-1)*m_top(:,1) + mV_top(j+n_top-1)*m_top(:,np_top);
333       
334#               "Tray energy balance"
335#           V(j)*H(j) = L(j-1)*h(j-1) + V(j+1)*H(j+1) - L(j)*h(j);
336                V_top(j+1) = V_top(j);
337        end
338
339        for j in [2:n1_bot] do
340           "Tray mass balance"
341#           L(j) = L(j-1) + V(j+1) - V(j);
342            L_bot(j) = L_bot(j-1);
343
344                for i in [1:NComp-1] do
345                        "Tray component mass balance"
346                        diff(mb_bot(i,j)) = sum(Bn_bot([1:np_bot]+(j-2)*np_bot)*L_bot*x_bot(i,:)) +
347                                                   sum(Bm_bot([1:np_bot]+(j-2)*np_bot)*V_bot*y_bot(i,:)) -
348                                                   L_bot(j)*x_bot(i,j) - V_bot(j)*y_bot(i,j) -
349                                                   mV_bot(j-1)*(L_bot(j)*x_bot(i,1) + V_bot(j)*y_bot(i,1)) -
350                                                   mV_bot(j+n_bot-1)*(L_bot(j)*x_bot(i,np_bot) + V_bot(j)*y_bot(i,np_bot));
351                        #RES_bot(i,j)= diff(mb_bot(i,j));
352                end
353               
354                mb_bot(:,j) = m_bot(:,j) + mV_bot(j-1)*m_bot(:,1) + mV_bot(j+n_bot-1)*m_bot(:,np_bot);
355       
356#               "Tray energy balance"
357#           V(j)*H(j) = L(j-1)*h(j-1) + V(j+1)*H(j+1) - L(j)*h(j);
358                V_bot(j+1) = V_bot(j);
359        end
360
361        "Feed tray mass balance"
362#    L(Nf) = F + L(Nfm1) + V(Nfp1) - V(Nf);
363    L_top(np_top) = L_top(n1_top);
364
365        for i in [1:NComp-1] do
366                "Feed tray component mass balance"
367                diff(mb_top(i,np_top)) = F*z(i) + sum(An_top([1:np_top]+np_top+Ns1_top)*L_top*x_top(i,:)) +
368                                                    sum(Am_bot([1:np_bot])*V_bot*y_bot(i,:)) -
369                                                        L_bot(1)*x_bot(i,1) - V_top(np_top)*y_top(i,np_top);
370                #RES_top(i,np_top)= diff(mb_top(i,np_top));
371        end
372
373        mb_top(:,np_top) = m_top(:,np_top);
374
375        m_bot(:,1) = m_top(:,np_top);
376        mb_bot(:,1) = m_bot(:,1);
377
378        x_bot(:,1) = x_top(:,np_top);
379#       m_bot(:,1) = M_bot(1) * x_bot(:,1);
380        y_bot(:,1) = y_top(:,np_top);
381        K_bot(:,1) = K_top(:,np_top);
382        T_bot(1) = T_top(np_top);
383
384        L_bot(1) = F + L_top(np_top);
385
386#       "Feed tray energy balance"
387#       V(Nf)*H(Nf) = F*hf + L(Nfm1)*h(Nfm1) + V(Nfp1)*H(Nfp1) - L(Nf)*h(Nf);
388    V_bot(1) = V_top(np_top);
389    V_bot(2) = V_bot(1);
390
391    "Reboiler mass balance"
392#   L(N2) = L(N1) - V(N2);
393    L_bot(np_bot) = L_bot(n1_bot);
394    B = F - V0 - D;
395
396        for i in [1:NComp-1] do
397    "Reboiler component mass balance"
398                diff(mb_bot(i,np_bot)) = sum(An_bot([1:np_bot]+np_bot+Ns1_bot)*L_bot*x_bot(i,:)) -
399                                                        B*x_bot(i,np_bot) - V_bot(np_bot)*y_bot(i,np_bot);
400                #RES_bot(i,np_bot)= diff(mb_bot(i,np_bot));
401        end
402
403        mb_bot(:,np_bot) = m_bot(:,np_bot);
404
405#       "Reboiler energy balance"
406#   -Qr = L(N1)*h(N1) - L(N2)*h(N2) - V(N2)*H(N2);
407
408        for j in [1:np_top] do
409                sum(x_top(:,j)) = 1;
410               
411                "Mol fraction"
412                m_top(:,j) = M_top(j) * x_top(:,j);
413               
414#               "Liquid Enthalpy"
415#               h(j) = PP.LiquidEnthalpy(T(j), P(j), x(:,j));
416
417#               "Vapour Enthalpy"
418#               H(j) = PP.VapourEnthalpy(T(j), P(j), y(:,j));
419               
420                "Chemical Equilibrium"
421                PP.LiquidFugacityCoefficient(T_top(j), P_top(j), x_top(:,j)) =
422                PP.VapourFugacityCoefficient(T_top(j), P_top(j), y_top(:,j))*K_top(:,j);
423               
424                y_top(:,j) = K_top(:,j)*x_top(:,j);
425
426        #       "Liquid mol fraction normalization"
427        #       sum(x(:,j)) = 1;
428
429        #       "Vapour mol fraction normalization"
430        #       sum(x(:,j)) = sum(y(:,j));
431                T_top(j) = PP.BubbleT(P_top(j),x_top(:,j));
432        end
433
434        for j in [2:np_bot] do
435                sum(x_bot(:,j)) = 1;
436               
437                "Mol fraction"
438                m_bot(:,j) = M_bot(j) * x_bot(:,j);
439
440#               "Liquid Enthalpy"
441#               h(j) = PP.LiquidEnthalpy(T(j), P(j), x(:,j));
442
443#               "Vapour Enthalpy"
444#               H(j) = PP.VapourEnthalpy(T(j), P(j), y(:,j));
445               
446                "Chemical Equilibrium"
447                PP.LiquidFugacityCoefficient(T_bot(j), P_bot(j), x_bot(:,j)) =
448                PP.VapourFugacityCoefficient(T_bot(j), P_bot(j), y_bot(:,j))*K_bot(:,j);
449               
450                y_bot(:,j) = K_bot(:,j)*x_bot(:,j);
451
452        #       "Liquid mol fraction normalization"
453        #       sum(x(:,j)) = 1;
454
455        #       "Vapour mol fraction normalization"
456        #       sum(x(:,j)) = sum(y(:,j));
457                T_bot(j) = PP.BubbleT(P_bot(j),x_bot(:,j));
458        end
459
460#       "Feed flash Calculation"
461#       [vf, xf, yf] = PP.FlashPH(Pf, hf, z);
462
463#       "Feed enthalpy"
464#       hf = (1-vf)*PP.LiquidEnthalpy(Tf, Pf, xf) + vf*PP.VapourEnthalpy(Tf, Pf, yf);
465
466        x(:,1) = x_top(:,1);
467        y(:,1) = y_top(:,1);
468        T(1) = T_top(1);
469
470        x(:,Nf) = x_top(:,np_top);
471        y(:,Nf) = y_top(:,np_top);
472        T(Nf) = T_top(np_top);
473
474        x(:,N2) = x_bot(:,np_bot);
475        y(:,N2) = y_bot(:,np_bot);
476        T(N2) = T_bot(np_bot);
477       
478# Interpolations
479        for j in [2:Nfm1] do
480                for i in [1:NComp-1] do
481                        x(i,j) = Disc_top.Interpol(x_top(i,:),j);
482                        y(i,j) = Disc_top.Interpol(y_top(i,:),j);
483                end
484                sum(x(:,j)) = 1;
485                sum(y(:,j)) = 1;
486                T(j) = Disc_top.Interpol(T_top/'K',j)*'K';
487        end
488
489        for j in [Nfp1:N1] do
490                for i in [1:NComp-1] do
491                        x(i,j) = Disc_bot.Interpol(x_bot(i,:),j);
492                        y(i,j) = Disc_bot.Interpol(y_bot(i,:),j);
493                end
494                sum(x(:,j)) = 1;
495                sum(y(:,j)) = 1;
496                T(j) = Disc_bot.Interpol(T_bot/'K',j)*'K';
497        end
498       
499###########################Cálculo dos Residuos ##########################################
500       
501        "Reflux ratio"
502        RR * (D + V(1)) = L(1);
503
504    "Condenser mass balance"
505        V(1) = 0* 'kmol/d';
506        V(2) = L(1) + D + V(1);
507               
508        "Condenser component mass balance"
509        RES(:,1) = V(2)*y(:,2) - (L(1) + D)*x(:,1) - V(1)*y(:,1);
510       
511        for j in [2:Nfm1,Nfp1:N1] do
512       
513        "Tray mass balance"
514                L(j) = L(j-1);
515
516        "Tray component mass balance"
517                RES(:,j) = L(j-1)*x(:,j-1) + V(j+1)*y(:,j+1) - L(j)*x(:,j) - V(j)*y(:,j);
518       
519                V(j+1) = V(j);
520       
521        end
522
523        "Feed tray mass balance"
524    L(Nf) = F + L(Nfm1);
525
526        "Feed tray component mass balance"
527    RES(:,Nf) = F*z + L(Nfm1)*x(:,Nfm1) + V(Nfp1)*y(:,Nfp1) - L(Nf)*x(:,Nf) - V(Nf)*y(:,Nf);
528
529    V(Nfp1) = V(Nf);
530
531    "Reboiler mass balance"
532    L(N2) = F - V(1) - D;
533
534    "Reboiler component mass balance"
535    RES(:,N2) = L(N1)*x(:,N1) - L(N2)*x(:,N2) - V(N2)*y(:,N2);
536       
537       
538##################################################################################
539       
540end
541
542
543Model Coluna_redcol_fix2
544       
545        PARAMETERS
546outer   PP      as Plugin(Brief="Physical Properties",Type="PP");
547outer   Disc_top as Plugin(Type="Discrete");
548outer   Disc_bot as Plugin(Type="Discrete");
549outer   NComp   as Integer;
550        NTrays  as Integer(Brief="Number of trays");
551        N1              as Integer;
552        N2              as Integer;
553        Nf              as Integer;
554        Nfm1    as Integer;
555        Nfp1    as Integer;             
556        np_top as Integer;
557        n_top as Integer;
558        n1_top as Integer;
559        n2_top as Integer;
560        Ns_top as Integer;
561        Ns1_top as Integer;
562        M_top(np_top)   as mol;
563        Am_top(Ns_top) as Real;
564        An_top(Ns_top) as Real;
565        np_bot as Integer;
566        n_bot as Integer;
567        n1_bot as Integer;
568        n2_bot as Integer;
569        Ns_bot as Integer;
570        Ns1_bot as Integer;
571        M_bot(np_bot)   as mol;
572        Am_bot(Ns_bot) as Real;
573        An_bot(Ns_bot) as Real;
574
575        r_top(np_top) as Real(Brief="raizes de hahn");
576        rs_top(np_top) as Real(Brief="raizes de hahn desnormal.");
577        w_top(np_top) as Real(Brief="quadrature weights");
578        Norm_top  as Real(Brief="scaling");
579        r_bot(np_bot) as Real(Brief="raizes de hahn");
580        rs_bot(np_bot) as Real(Brief="raizes de hahn desnormal.");
581        w_bot(np_bot) as Real(Brief="quadrature weights");
582        Norm_bot as Real(Brief="scaling");
583       
584        SET
585        NTrays = Disc_bot.UpperBound-Disc_top.LowerBound-1;
586        N2 = NTrays + 2;
587        Nf = Disc_top.UpperBound;
588        np_top = Disc_top.NodalPoints;
589        n_top = Disc_top.InternalPoints;
590        n1_top = n_top + 1;
591        n2_top = 2*n_top;
592        Ns_top= np_top * np_top;
593        Ns1_top= n_top * np_top;
594        Am_top = Disc_top.Aplus;       
595        An_top = Disc_top.Aminus;
596        np_bot = Disc_bot.NodalPoints;
597        n_bot = Disc_bot.InternalPoints;
598        n1_bot = n_bot + 1;
599        n2_bot = 2*n_bot;
600        Ns_bot= np_bot * np_bot;
601        Ns1_bot= n_bot * np_bot;
602        Am_bot = Disc_bot.Aplus;       
603        An_bot = Disc_bot.Aminus;
604
605        r_top = Disc_top.Roots;
606        w_top = Disc_top.Weights;
607        Norm_top = Disc_top.Scaling;
608        rs_top = r_top * Norm_top;
609        r_bot = Disc_bot.Roots;
610        w_bot = Disc_bot.Weights;
611        Norm_bot = Disc_bot.Scaling;
612        rs_bot = r_bot * Norm_bot;
613
614        M_top(1) = 20 * 'kmol';
615        M_top(2:np_top) = 0.4 * 'kmol';
616        M_bot(1:n1_bot) = 0.4 * 'kmol';
617        M_bot(np_bot) = 25 * 'kmol';
618       
619        ##############################
620        N1 = NTrays + 1;
621        N2 = NTrays + 2;
622        Nfm1 = Nf - 1;
623        Nfp1 = Nf + 1;
624    ###############################
625
626        VARIABLES
627        F as flow_mol;
628        D as flow_mol;
629        V0 as flow_mol;
630        z(NComp) as fraction;
631        Tf as temperature;
632        Pf as pressure;
633        L(N2) as flow_mol;
634        V(N2) as flow_mol;
635
636        L_top(np_top) as flow_mol;
637        V_top(np_top) as flow_mol;
638        T_top(np_top) as temperature;
639        P_top(np_top) as pressure;
640        x_top(NComp,np_top) as fraction(Lower=-1,Upper=2);
641        y_top(NComp,np_top) as fraction(Lower=-1,Upper=2);
642        K_top(NComp,np_top) as positive;
643        m_top(NComp,np_top) as mol(Lower=-1e3);
644        RES(NComp,N2)    as flow_mol(Lower=-1e6);
645
646        L_bot(np_bot) as flow_mol;
647        V_bot(np_bot) as flow_mol;
648        T_bot(np_bot) as temperature;
649        P_bot(np_bot) as pressure;
650        x_bot(NComp,np_bot) as fraction(Lower=-1,Upper=2);
651        y_bot(NComp,np_bot) as fraction(Lower=-1,Upper=2);
652        K_bot(NComp,np_bot) as positive;
653        m_bot(NComp,np_bot) as mol(Lower=-1e3);
654
655        x(NComp,N2) as fraction(Lower=-1,Upper=2);
656        y(NComp,N2) as fraction(Lower=-1,Upper=2);
657        T(N2) as temperature;
658       
659
660#       H(N2) as enth_mol;
661#       h(N2) as enth_mol;
662#       hf        as enth_mol;
663#       Qc        as heat_rate;
664#       Qr        as heat_rate;
665#       vf        as fraction;
666#       xf(NComp) as fraction;
667#       yf(NComp) as fraction;
668        RR        as positive;
669
670        EQUATIONS
671
672        "Reflux ratio"
673        RR * (D + V0) = L_top(1);
674
675    "Condenser mass balance"
676    V_top(1) = V0;
677        V_top(2) = L_top(1) + D + V0;
678
679        for i in [1:NComp-1] do
680                "Condenser component mass balance"
681                diff(m_top(i,1)) = V_top(2)*sum(Am_top([1:np_top])*y_top(i,:)) - (L_top(1) + D)*x_top(i,1) - V0*y_top(i,1);
682        end
683
684#       "Reboiler energy balance"
685#   Qc = V(2)*H(2) - (L(1) + D)*h(1) - V(1)*H(1);
686
687        for j in [2:n1_top] do
688           "Tray mass balance"
689#           L(j) = L(j-1) + V(j+1) - V(j);
690            L_top(j) = L_top(j-1);
691
692                for i in [1:NComp-1] do
693                        "Tray component mass balance"
694                        diff(m_top(i,j)) = L_top(j-1)*sum(An_top([1:np_top]+(j-1)*np_top)*x_top(i,:)) +
695                                                   V_top(j+1)*sum(Am_top([1:np_top]+(j-1)*np_top)*y_top(i,:)) -
696                                                   L_top(j)*x_top(i,j) - V_top(j)*y_top(i,j);
697                end
698               
699#               "Tray energy balance"
700#           V(j)*H(j) = L(j-1)*h(j-1) + V(j+1)*H(j+1) - L(j)*h(j);
701                V_top(j+1) = V_top(j);
702        end
703
704        for j in [2:n1_bot] do
705           "Tray mass balance"
706#           L(j) = L(j-1) + V(j+1) - V(j);
707            L_bot(j) = L_bot(j-1);
708
709                for i in [1:NComp-1] do
710                        "Tray component mass balance"
711                        diff(m_bot(i,j)) = L_bot(j-1)*sum(An_bot([1:np_bot]+(j-1)*np_bot)*x_bot(i,:)) +
712                                                   V_bot(j+1)*sum(Am_bot([1:np_bot]+(j-1)*np_bot)*y_bot(i,:)) -
713                                                   L_bot(j)*x_bot(i,j) - V_bot(j)*y_bot(i,j);
714                end
715               
716#               "Tray energy balance"
717#           V(j)*H(j) = L(j-1)*h(j-1) + V(j+1)*H(j+1) - L(j)*h(j);
718                V_bot(j+1) = V_bot(j);
719        end
720
721        "Feed tray mass balance"
722#    L(Nf) = F + L(Nfm1) + V(Nfp1) - V(Nf);
723    L_top(np_top) = F + L_top(n1_top);
724
725        for i in [1:NComp-1] do
726                "Feed tray component mass balance"
727                diff(m_top(i,np_top)) = F*z(i) + L_top(n1_top)*sum(An_top([1:np_top]+np_top+Ns1_top)*x_top(i,:)) +
728                                                    V_bot(1)*sum(Am_bot([1:np_bot])*y_bot(i,:)) -
729                                                        L_top(np_top)*x_top(i,np_top) - V_bot(1)*y_bot(i,1);
730        end
731
732        m_bot(:,1) = m_top(:,np_top);
733
734        x_bot(:,1) = x_top(:,np_top);
735#       m_bot(:,1) = M_bot(1) * x_bot(:,1);
736        y_bot(:,1) = y_top(:,np_top);
737        K_bot(:,1) = K_top(:,np_top);
738        T_bot(1) = T_top(np_top);
739
740        L_bot(1) = L_top(np_top);
741
742#       "Feed tray energy balance"
743#       V(Nf)*H(Nf) = F*hf + L(Nfm1)*h(Nfm1) + V(Nfp1)*H(Nfp1) - L(Nf)*h(Nf);
744    V_bot(1) = V_top(np_top);
745    V_bot(2) = V_bot(1);
746
747    "Reboiler mass balance"
748#   L(N2) = L(N1) - V(N2);
749    L_bot(np_bot) = F - V0 - D;
750
751        for i in [1:NComp-1] do
752    "Reboiler component mass balance"
753                diff(m_bot(i,np_bot)) = L_bot(n1_bot)*sum(An_bot([1:np_bot]+np_bot+Ns1_bot)*x_bot(i,:)) -
754                                                        L_bot(np_bot)*x_bot(i,np_bot) - V_bot(np_bot)*y_bot(i,np_bot);
755        end
756
757#       "Reboiler energy balance"
758#   -Qr = L(N1)*h(N1) - L(N2)*h(N2) - V(N2)*H(N2);
759
760        for j in [1:np_top] do
761                sum(x_top(:,j)) = 1;
762               
763                "Mol fraction"
764                m_top(:,j) = M_top(j) * x_top(:,j);
765               
766#               "Liquid Enthalpy"
767#               h(j) = PP.LiquidEnthalpy(T(j), P(j), x(:,j));
768
769#               "Vapour Enthalpy"
770#               H(j) = PP.VapourEnthalpy(T(j), P(j), y(:,j));
771               
772                "Chemical Equilibrium"
773                PP.LiquidFugacityCoefficient(T_top(j), P_top(j), x_top(:,j)) =
774                PP.VapourFugacityCoefficient(T_top(j), P_top(j), y_top(:,j))*K_top(:,j);
775               
776                y_top(:,j) = K_top(:,j)*x_top(:,j);
777
778        #       "Liquid mol fraction normalization"
779        #       sum(x(:,j)) = 1;
780
781        #       "Vapour mol fraction normalization"
782        #       sum(x(:,j)) = sum(y(:,j));
783                T_top(j) = PP.BubbleT(P_top(j),x_top(:,j));
784        end
785
786        for j in [2:np_bot] do
787                sum(x_bot(:,j)) = 1;
788               
789                "Mol fraction"
790                m_bot(:,j) = M_bot(j) * x_bot(:,j);
791
792#               "Liquid Enthalpy"
793#               h(j) = PP.LiquidEnthalpy(T(j), P(j), x(:,j));
794
795#               "Vapour Enthalpy"
796#               H(j) = PP.VapourEnthalpy(T(j), P(j), y(:,j));
797               
798                "Chemical Equilibrium"
799                PP.LiquidFugacityCoefficient(T_bot(j), P_bot(j), x_bot(:,j)) =
800                PP.VapourFugacityCoefficient(T_bot(j), P_bot(j), y_bot(:,j))*K_bot(:,j);
801               
802                y_bot(:,j) = K_bot(:,j)*x_bot(:,j);
803
804        #       "Liquid mol fraction normalization"
805        #       sum(x(:,j)) = 1;
806
807        #       "Vapour mol fraction normalization"
808        #       sum(x(:,j)) = sum(y(:,j));
809                T_bot(j) = PP.BubbleT(P_bot(j),x_bot(:,j));
810        end
811
812#       "Feed flash Calculation"
813#       [vf, xf, yf] = PP.FlashPH(Pf, hf, z);
814
815#       "Feed enthalpy"
816#       hf = (1-vf)*PP.LiquidEnthalpy(Tf, Pf, xf) + vf*PP.VapourEnthalpy(Tf, Pf, yf);
817
818# Interpolations
819        x(:,1) = x_top(:,1);
820        y(:,1) = y_top(:,1);
821        T(1) = T_top(1);
822
823        x(:,Nf) = x_top(:,np_top);
824        y(:,Nf) = y_top(:,np_top);
825        T(Nf) = T_top(np_top);
826
827        x(:,N2) = x_bot(:,np_bot);
828        y(:,N2) = y_bot(:,np_bot);
829        T(N2) = T_bot(np_bot);
830
831        for j in [2:Nfm1] do
832                for i in [1:NComp-1] do
833                        x(i,j) = Disc_top.Interpol(x_top(i,:),j);
834                        y(i,j) = Disc_top.Interpol(y_top(i,:),j);
835                end
836                sum(x(:,j)) = 1;
837                sum(y(:,j)) = 1;
838                T(j) = Disc_top.Interpol(T_top/'K',j)*'K';
839        end
840
841        for j in [Nfp1:N1] do
842                for i in [1:NComp-1] do
843                        x(i,j) = Disc_bot.Interpol(x_bot(i,:),j);
844                        y(i,j) = Disc_bot.Interpol(y_bot(i,:),j);
845                end
846                sum(x(:,j)) = 1;
847                sum(y(:,j)) = 1;
848                T(j) = Disc_bot.Interpol(T_bot/'K',j)*'K';
849        end
850       
851        ########Cálculo dos Residuos para identificação se os Balanços de massa fecham######
852       
853        "Reflux ratio"
854        RR * (D + V(1)) = L(1);
855
856    "Condenser mass balance"
857     V(1) = 0* 'kmol/d';
858     V(2) = L(1) + D + V(1);
859               
860        "Condenser component mass balance"
861        RES(:,1) = V(2)*y(:,2) - (L(1) + D)*x(:,1) - V(1)*y(:,1);
862       
863        for j in [2:Nfm1,Nfp1:N1] do
864           
865        "Tray mass balance"
866                L(j) = L(j-1);
867
868        "Tray component mass balance"
869           RES(:,j) = L(j-1)*x(:,j-1) + V(j+1)*y(:,j+1) - L(j)*x(:,j) - V(j)*y(:,j);
870       
871                V(j+1) = V(j);
872        end
873
874        "Feed tray mass balance"
875    L(Nf) = F + L(Nfm1);
876
877        "Feed tray component mass balance"
878    RES(:,Nf) = F*z + L(Nfm1)*x(:,Nfm1) + V(Nfp1)*y(:,Nfp1) - L(Nf)*x(:,Nf) - V(Nf)*y(:,Nf);
879
880    V(Nfp1) = V(Nf);
881
882    "Reboiler mass balance"
883    L(N2) = F - V(1) - D;
884
885    "Reboiler component mass balance"
886    RES(:,N2) = L(N1)*x(:,N1) - L(N2)*x(:,N2) - V(N2)*y(:,N2);
887               
888##################################################################################
889
890end
891
892FlowSheet Test_Coluna_classica_peq
893       
894        PARAMETERS
895        PP      as Plugin(Brief="Physical Properties",
896                Type="PP",
897                Components = ["propylene", "propane"],
898                LiquidModel = "PR",
899                VapourModel = "PR"
900        );
901        Disc_top as Plugin(Type="Discrete", Boundary="BOTH",
902                                        InternalPoints=1, Lower=1, Upper=12,
903                                        Normal=23);
904        Disc_bot as Plugin(Type="Discrete", Boundary="BOTH",
905                                        InternalPoints=1, Lower=12, Upper=23,
906                                        Normal=23);
907        NComp   as Integer;
908
909        SET
910        NComp = PP.NumberOfComponents; 
911
912        DEVICES
913        col as Coluna_redcol_fix2;
914
915        SPECIFY
916        col.F = 1073.4 * 'kmol/d';
917        col.Tf = (46.11+273.15) * 'K';
918        col.Pf = 1860.60 * 'kPa';
919        col.z = [0.8973, 0.1027];
920        #col.P(1) = 1860.60*'kPa';
921        #col.Qc = 2800 * 'kW';
922        col.P_top = 1860.60*'kPa';
923        col.P_bot = 1860.60*'kPa';
924        col.D = 965 * 'kmol/d';
925        col.V0 = 0 * 'kmol/d';
926        #col.RR = 19.7;
927       
928        EQUATIONS
929       
930                #Disturbance
931        #if time < 2*'h' then
932                col.RR = 19.7;
933        #else if time < 50*'h' then
934        #       col.RR = 22;
935        #else if time < 75*'h' then
936        #       col.RR = 19.7;
937        #else if time < 110*'h' then
938        #       col.RR = 18.7;
939        #else
940        #       col.RR = 19.7;
941        #end
942    #end
943        #end
944        #end
945       
946        GUESS
947#       col.L = 19.7 * 965 * 'kmol/d';
948#       col.V = (19.7+1) * 965 * 'kmol/d';
949#       col.T = [320:(330-320)/col.N1:320] * 'K';
950#       col.T = 320 * 'K';
951#       col.y = 0.5;
952
953        INITIAL
954#       col.T = [320:(330-320)/col.N1:320] * 'K';
955        col.x_top(1,:) = [0.8:(0.5-0.8)/col.n1_top:0.5];
956        col.x_bot(1,2:col.np_bot) = [0.5:(0.2-0.5)/col.n_bot:0.2];
957#       col.x_top(2,:) = 1-col.x_top(1,:);
958#       col.x_bot(2,:) = 1-col.x_bot(1,:);
959
960        OPTIONS
961       
962        # Rodar caso Estacionario       
963        #Dynamic        = true;
964        #Dynamic        = false; 
965        TimeStep        = 0.2; #0.05;
966        TimeEnd         = 5; # colocar igual ao número de estágios internos
967        TimeUnit        = 'h';
968        #GuessFile = "RET_6_ESG_5";
969       
970        # Rodar caso Dinâmico   
971#       Dynamic         = true; 
972#       TimeStep        = 0.1;
973#       TimeEnd         = 12; # colocar igual ao número de estágios internos
974#       TimeUnit        = 'h';
975#       InitialFile = "C:\DOCUMENTOS\Dissertação\Simulações_EMSO\estacionario\Classica\RET_5_ESG_3";
976       
977#       DAESolver(RelativeAccuracy = 1e-12,AbsoluteAccuracy = 1e-14);
978#       NLASolver(RelativeAccuracy = 1e-12,AbsoluteAccuracy = 1e-14);
979        DAESolver(RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
980        NLASolver(RelativeAccuracy = 1e-7,AbsoluteAccuracy = 1e-9);
981#       DAESolver(RelativeAccuracy = 1e-10,AbsoluteAccuracy = 1e-12);
982#       NLASolver(RelativeAccuracy = 1e-10,AbsoluteAccuracy = 1e-12);
983
984       
985       
986#       Dynamic         = false; 
987#       TimeStep        = 0.2;
988#       TimeEnd         = 5; # colocar igual ao número de estágios internos
989#       TimeUnit        = 'h';
990#       InitialFile = "Test_Coluna_reduzida_fix10";
991#       GuessFile = "Test_Coluna_reduzida_fix20";
992#       GuessFile = "Test_Coluna_redcol_fix";
993#       InitialFile = "Test_Coluna_classica_fix";
994#       DAESolver(RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
995#       NLASolver(RelativeAccuracy = 1e-7,AbsoluteAccuracy = 1e-9);
996end
997
998
999FlowSheet Test_Coluna_momentos_peq
1000       
1001        PARAMETERS
1002        PP      as Plugin(Brief="Physical Properties",
1003                Type="PP",
1004                Components = ["propylene", "propane"],
1005                LiquidModel = "PR",
1006                VapourModel = "PR"
1007        );
1008        Disc_top as Plugin(Type="Discrete", Boundary="BOTH",
1009                                        InternalPoints=1, Lower=1, Upper=12,
1010                                        Normal=23);
1011        Disc_bot as Plugin(Type="Discrete", Boundary="BOTH",
1012                                        InternalPoints=1, Lower=12, Upper=23,
1013                                        Normal=23);
1014        NComp   as Integer;
1015
1016        SET
1017        NComp = PP.NumberOfComponents; 
1018
1019        DEVICES
1020        col as Coluna_reduzida_fix2;
1021
1022        SPECIFY
1023        col.F = 1073.4 * 'kmol/d';
1024        col.Tf = (46.11+273.15) * 'K';
1025        col.Pf = 1860.60 * 'kPa';
1026        col.z = [0.8973, 0.1027];
1027        #col.P(1) = 1860.60*'kPa';
1028        #col.Qc = 2800 * 'kW';
1029        col.P_top = 1860.60*'kPa';
1030        col.P_bot = 1860.60*'kPa';
1031        col.D = 965 * 'kmol/d';
1032        col.V0 = 0 * 'kmol/d';
1033#       col.RR = 19.7;
1034
1035        EQUATIONS
1036       
1037        #Disturbance
1038        #if time < 2*'h' then
1039                col.RR = 19.7;
1040        #else if time < 50*'h' then
1041        #       col.RR = 22;
1042        #else if time < 75*'h' then
1043        #       col.RR = 19.7;
1044        #else if time < 110*'h' then
1045        #       col.RR = 18.7;
1046        #else
1047        #       col.RR = 19.7;
1048        #end
1049    #end
1050        #end
1051        #end
1052       
1053       
1054
1055        GUESS
1056#       col.L = 19.7 * 965 * 'kmol/d';
1057#       col.V = (19.7+1) * 965 * 'kmol/d';
1058#       col.T = [320:(330-320)/col.N1:320] * 'K';
1059#       col.T = 320 * 'K';
1060#       col.y = 0.5;
1061
1062        INITIAL
1063#       col.T = [320:(330-320)/col.N2:350] * 'K';
1064        col.x_top(1,:) = [0.8:(0.5-0.8)/col.n1_top:0.5];
1065        col.x_bot(1,2:col.np_bot) = [0.5:(0.2-0.5)/col.n_bot:0.2];
1066#       col.x_top(2,:) = 1-col.x_top(1,:);
1067#       col.x_bot(2,:) = 1-col.x_bot(1,:);
1068
1069        OPTIONS
1070       
1071        # Rodar caso Estacionario       
1072        #Dynamic        = true;
1073        #Dynamic        = false; 
1074        TimeStep        = 0.2;
1075        TimeEnd         = 5; # colocar igual ao número de estágios internos
1076        TimeUnit        = 'h';
1077#       GuessFile = "Test_Coluna_momentos_peq";
1078#       GuessFile = "C:\DOCUMENTOS\Dissertação\Simulações_EMSO\estacionario\Momento\Momentos_8_8";
1079
1080        # Rodar caso Dinâmico   
1081#       Dynamic         = true; 
1082#       TimeStep        = 0.1;
1083#       TimeEnd         = 12; # colocar igual ao número de estágios internos
1084#       TimeUnit        = 'h';
1085#       InitialFile = "C:\DOCUMENTOS\Dissertação\Simulações_EMSO\estacionario\Momento\RET_5_ESG_3";
1086       
1087#       DAESolver(RelativeAccuracy = 1e-14,AbsoluteAccuracy = 1e-16);
1088#       NLASolver(RelativeAccuracy = 1e-12,AbsoluteAccuracy = 1e-14);
1089        DAESolver(RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
1090        NLASolver(RelativeAccuracy = 1e-7,AbsoluteAccuracy = 1e-9);
1091        #DAESolver(RelativeAccuracy = 1e-10,AbsoluteAccuracy = 1e-12);
1092        #NLASolver(RelativeAccuracy = 1e-10,AbsoluteAccuracy = 1e-12);
1093       
1094end
1095
1096FlowSheet Test_Coluna_completa_peq
1097       
1098        PARAMETERS
1099        PP      as Plugin(Brief="Physical Properties",
1100                Type="PP",
1101                Components = ["propylene", "propane"],
1102                LiquidModel = "PR",
1103                VapourModel = "PR"
1104        );
1105        NComp   as Integer;
1106
1107        DEVICES
1108        col as Coluna_completa_fix2;
1109       
1110        SET
1111        NComp = PP.NumberOfComponents; 
1112        col.NTrays = 21;
1113        col.Nf = 12;
1114
1115        SPECIFY
1116        col.F = 1073.4 * 'kmol/d';
1117        col.Tf = (46.11+273.15) * 'K';
1118        col.Pf = 1860.60 * 'kPa';
1119        col.z = [0.8973, 0.1027];
1120        #col.P(1) = 1860.60*'kPa';
1121        #col.Qc = 2800 * 'kW';
1122        col.P = 1860.60*'kPa';
1123        col.D = 965 * 'kmol/d';
1124        col.V(1) = 0 * 'kmol/d';
1125#       col.RR = 19.7;
1126       
1127        EQUATIONS
1128       
1129        #Disturbance
1130        #if time < 2*'h' then
1131                col.RR = 19.7;
1132        #else if time < 50*'h' then
1133        #       col.RR = 22;
1134        #else if time < 75*'h' then
1135        #       col.RR = 19.7;
1136        #else if time < 110*'h' then
1137        #       col.RR = 18.7;
1138        #else
1139        #       col.RR = 19.7;
1140        #end
1141    #end
1142        #end
1143        #end
1144       
1145        GUESS
1146#       col.L = 19.7 * 965 * 'kmol/d';
1147#       col.V = (19.7+1) * 965 * 'kmol/d';
1148#       col.T = [320:(330-320)/col.N1:320] * 'K';
1149#       col.T = 320 * 'K';
1150#       col.y = 0.5;
1151
1152        INITIAL
1153#       col.T = [320:(330-320)/col.N1:320] * 'K';
1154        col.x(1,:) = [0.8:(0.2-0.8)/col.N1:0.2];
1155        col.x(2,:) = 1-col.x(1,:);
1156
1157        OPTIONS
1158
1159                # Rodar caso Estacionario       
1160        Dynamic         = true;
1161        #Dynamic        = false; 
1162        TimeStep        = 0.2;
1163        TimeEnd         = 5; # colocar igual ao número de estágios internos
1164        TimeUnit        = 'h';
1165        #GuessFile = "Completo";
1166       
1167       
1168        # Rodar caso Dinâmico   
1169#       Dynamic         = true; 
1170#       TimeStep        = 0.1;
1171#       TimeEnd         = 12; # colocar igual ao número de estágios internos
1172#       TimeUnit        = 'h';
1173#       InitialFile = "C:\DOCUMENTOS\Dissertação\Simulações_EMSO\estacionario\Completo\Completo";
1174       
1175#       DAESolver(RelativeAccuracy = 1e-14,AbsoluteAccuracy = 1e-16);
1176#       NLASolver(RelativeAccuracy = 1e-12,AbsoluteAccuracy = 1e-14);
1177        DAESolver(RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
1178        NLASolver(RelativeAccuracy = 1e-7,AbsoluteAccuracy = 1e-9);
1179        #DAESolver(RelativeAccuracy = 1e-10,AbsoluteAccuracy = 1e-12);
1180        #NLASolver(RelativeAccuracy = 1e-10,AbsoluteAccuracy = 1e-12);
1181
1182
1183
1184
1185        #Dynamic        = false; 
1186        #TimeStep       = 0.1;
1187        #TimeEnd        = 12; # colocar igual ao número de estágios internos
1188        #TimeUnit       = 'h';
1189        #InitialFile = "Completo/Test_Coluna_completa_fix2";
1190        #DAESolver(RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
1191        #NLASolver(RelativeAccuracy = 1e-7,AbsoluteAccuracy = 1e-9);
1192end
Note: See TracBrowser for help on using the repository browser.