hydroscal_model.hpp
Go to the documentation of this file.
00001 /*
00002  *    This file is part of ACADO Toolkit.
00003  *
00004  *    ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
00005  *    Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
00006  *    Milan Vukov, Rien Quirynen, KU Leuven.
00007  *    Developed within the Optimization in Engineering Center (OPTEC)
00008  *    under supervision of Moritz Diehl. All rights reserved.
00009  *
00010  *    ACADO Toolkit is free software; you can redistribute it and/or
00011  *    modify it under the terms of the GNU Lesser General Public
00012  *    License as published by the Free Software Foundation; either
00013  *    version 3 of the License, or (at your option) any later version.
00014  *
00015  *    ACADO Toolkit is distributed in the hope that it will be useful,
00016  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018  *    Lesser General Public License for more details.
00019  *
00020  *    You should have received a copy of the GNU Lesser General Public
00021  *    License along with ACADO Toolkit; if not, write to the Free Software
00022  *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00023  *
00024  */
00025 
00026 
00027 
00056 //#define PSEUDO_STATES
00057 /* problem setup */
00058 
00059 /* #define  NMOS  2*/
00060 
00061 #define  NMOS   2
00062 #define  NP     17
00063 #define  NRC    0
00064 #define  NRCE   0
00065 
00066 #define  NPR    0
00067 #define  NLSQ   4 
00068 
00069 #ifdef PSEUDO_STATES
00070 #define  NXD    84
00071 #else
00072 #define  NXD    82
00073 #endif
00074 
00075 #define  NXA    122
00076 #define  NU     2
00077 #define  NW     1
00078 
00079 /* model parameters */
00080 
00081 #define  qF      1
00082 #define  Nfeed   20
00083 
00084 #define A1       23.48
00085 #define B1       3626.6
00086 #define C1       -34.29
00087 
00088 #define A2       22.437
00089 #define B2       3166.4
00090 #define C2       -80.15
00091 
00092 #define Rconst        8.3147
00093 
00094 #define T14s    (70)
00095 #define T28s    (88)
00096 
00097 #define QT14    1    
00098 #define QT28    1    
00099                      
00100 #define  R11  (2/2)
00101 #define  R22  (0.1/2)
00102 #define  RPENALTY  (0.1)
00103 
00104 
00105 typedef struct {
00106 
00107     double x[42]; 
00108     double n[42]; 
00109     double L[42]; 
00110     double V[42]; 
00111     double Temp[42];
00112 
00113     double refvol[42];
00114     double Pres[42];
00115     double P1[42];
00116     double P2[42];
00117     double y[42];
00118 
00119     double F; 
00120     double xF; 
00121 
00122 
00123     double  nvolB;
00124     double  nvolD;
00125     double  nvolM;
00126     double  nvolK;
00127     double  flowwidth;
00128 
00129     double  alphastrip;
00130     double  alpharect;
00131     double  Ptop;
00132     double  DPstrip;
00133     double  DPrect;
00134     double  Tfeed;
00135     double  Tcond;
00136 
00137 } Parameters;
00138 
00139 static Parameters par;
00140 
00141 
00142 
00143 
00144 static double T_x(double T, double x)
00145      /*
00146       Calculates derivative of T with respect to x, 
00147       according to implicit function theorem for
00148       0=F=x*P1(T)+(1-x)*P2(T)-P
00149       */
00150 {
00151      double P1,P2,P1_T,P2_T, F_T, F_x, T_x;
00152 
00153      P1   = exp(A1-B1/(C1+T));  
00154      P1_T = P1 * B1/(C1+T)/(C1+T); 
00155 
00156      P2   = exp(A2-B2/(C2+T));  
00157      P2_T = P2 * B2/(C2+T)/(C2+T);
00158 
00159      F_T=x*P1_T+(1-x)*P2_T;
00160      F_x=P1-P2;
00161 
00162      T_x=- 1/F_T*F_x;
00163      return T_x;
00164 }
00165 
00166 
00167 static void enthalpies_1(double Temp, double Pres, double *hL1, double *hV1, double *hL1_T)
00168      /* 
00169         gives enthalpies of pure liquid/gaseous Methanol at Temperature Temp and Pressure Pres 
00170      */
00171 
00172 {
00173   const double  h11    =   18.31    ;
00174   const double  h21    =   1.713E-2 ;
00175   const double  h31    =   6.399E-5 ;
00176   const double  Tc1    =   512.6    ;
00177   const double  Pc1    =   8.096E6  ;
00178   const double  OMEGA1 =   0.557    ;
00179 
00180   double TR1, PR1, Dhv1;
00181 
00182   TR1=Temp/Tc1;
00183   PR1=Pres/Pc1;
00184   Dhv1 = Rconst*Tc1*sqrt(1-PR1/pow(TR1,3))*(6.09648-1.28862*TR1+
00185       1.016*pow(TR1,7)+OMEGA1*(15.6875-13.4721*TR1+2.615*pow(TR1,7)));
00186   *hL1=(h11*(Temp-273.15)+h21*pow(Temp-273.15,2)+h31*pow(Temp-273.15,3))*4.186;
00187   *hV1=*hL1+Dhv1;
00188   *hL1_T=(h11+h21*2*(Temp-273.15)+h31*3*pow(Temp-273.15,2))*4.186;
00189 }
00190 
00191 static void enthalpies_2(double Temp, double Pres, double *hL2, double *hV2, double *hL2_T)
00192      /* 
00193         gives enthalpies of pure liquid/gaseous n-Propanol at Temperature Temp and Pressure Pres 
00194      */
00195 
00196 {
00197   const double  h12     =   31.92      ;
00198   const double  h22     =   4.49E-2    ;
00199   const double  h32     =   9.663E-5   ;
00200   const double  Tc2     =   536.7      ;
00201   const double  Pc2     =   5.166E6    ;
00202   const double  OMEGA2  =   0.612      ;
00203 
00204   double TR2, PR2, Dhv2;
00205 
00206   TR2=Temp/Tc2;
00207   PR2=Pres/Pc2;
00208   Dhv2 = Rconst*Tc2*sqrt(1-PR2/pow(TR2,3))*(6.09648-1.28862*TR2+
00209       1.016*pow(TR2,7)+OMEGA2*(15.6875-13.4721*TR2+2.615*pow(TR2,7)));
00210   *hL2=(h12*(Temp-273.15)+h22*pow(Temp-273.15,2)+h32*pow(Temp-273.15,3))*4.186;
00211   *hV2=*hL2+Dhv2;
00212   *hL2_T=(h12+h22*2*(Temp-273.15)+h32*3*pow(Temp-273.15,2))*4.186;
00213 }
00214 
00215 
00216 
00217 static void enthalpies(double Temp, double Pres, double x, double y, 
00218                        double *hL, double *hV, double *hL_x_total)
00219      /* 
00220         gives enthalpies  of mixed liquid/vapour
00221         - and (total) derivative of liquid -
00222         at Temperature Temp and Pressure Pres 
00223      */
00224 
00225 {
00226   double hL1, hV1, hL1_T, hL2, hV2, hL2_T, hL_x, hL_T;
00227   enthalpies_1(Temp, Pres, &hL1,&hV1, &hL1_T);
00228   enthalpies_2(Temp, Pres, &hL2,&hV2, &hL2_T);
00229   *hL=hL1*x+hL2*(1-x);
00230   *hV=hV1*y+hV2*(1-y);
00231   hL_x=hL1-hL2;
00232   hL_T=hL1_T*x+hL2_T*(1-x);
00233   *hL_x_total=hL_x
00234               +hL_T*T_x(Temp,x)
00235              ;
00236 }
00237 
00238 
00239 static double litre_per_kmol(double TK, double x_m)
00240      /* gives molar Volume in 
00241              litre / kmol
00242         from 
00243              Temp in Kelvin, 
00244              x_m the Methanol concentration
00245      */
00246 
00247 {
00248   double d_m, vm, d_p, vp, vol_p, vol_m, vol_sum, v_mix;//, F;
00249   d_m =  2.288/ pow(0.2685,(1+ pow((1 - TK/512.4),0.2453))); 
00250   vm = 1/d_m;
00251 
00252   d_p =  1.235/ pow(0.27136, (1+ pow((1- TK/536.4),0.24)));
00253   vp = 1/d_p;
00254 
00255   vol_m = x_m*vm;
00256   vol_p = (1-x_m)*vp;
00257   vol_sum = vol_m + vol_p;
00258   v_mix= vol_sum*1000;
00259 
00260   return v_mix;
00261 }
00262 
00263 
00264 static double ramp(double x)
00265 {
00266   /*
00267    smooth version of the "__/" function 1/2 ( x + |x| ) 
00268    */
00269   const double R=20; 
00270   double xo;
00271    if (x<-R)  xo=0;
00272   else if (x<R) xo=0.5*(x+log(2*cosh(x)));
00273   else xo=x;
00274 
00275   return xo;
00276 
00277 }
00278 static double kmol_per_sec_outflow(double n, double TK, double x_m, double refvol, 
00279                                    double flowwidth)
00280      /* gives Tray outflow in 
00281              kmol / sec
00282         from 
00283              n molar tray holdup
00284              Temp in Kelvin, 
00285              x_m the Methanol concentration
00286              refvol reference volume in litre
00287      */
00288 
00289 {
00290   double volumeholdup, volumeoutflow, outflow, l_per_kmol;
00291   const double alpha=500;
00292   l_per_kmol=litre_per_kmol(TK, x_m);
00293   volumeholdup=l_per_kmol * n;
00294   volumeoutflow = flowwidth*(
00295                              pow(1/alpha*ramp(alpha*(volumeholdup-refvol)),1.5)
00296                              /* -1/pow(alpha,1.5)*ramp(-alpha*(volumeholdup-refvol)) */
00297                              );
00298   outflow=volumeoutflow / l_per_kmol;
00299 
00300   /*
00301   acadoFPrintf(stderr,"%g  ",refvol/volumeholdup);
00302   */
00303 
00304   return outflow;
00305 }
00306 
00307 
00308 
00309 static void insert(double *xd, 
00310                    double *xa,  
00311                    double *u, 
00312                    double *p
00313 )
00314      /* 
00315         Extracts variables out of xd, xa, u, p, puts them into
00316         par.x, par.Temp, par.V, par.L, etc.
00317        
00318 
00319         Attention: par.V[41] and par.L[0] are only specified in prepare(), as they need y.
00320 
00321      */
00322 {
00323   double Lc, Vc;
00324   double Q;//, R;
00325   double L_lh,F_lh;
00326   //double xC;
00327   double holdup = 0.0;
00328   double dummy, hL, hV, deltaH, Qloss;
00329   long i;
00330   extern Parameters par;
00331   
00332   L_lh   = u[0];
00333   Q      = u[1];
00334 
00335   par.alpharect  = p[2];
00336   par.alphastrip = p[3];
00337 
00338   Qloss   = p[5];
00339 
00340   par.nvolB      = p[6];
00341   par.nvolD      = p[7];
00342   par.nvolM      = p[0];
00343   par.nvolK      = p[0];
00344   par.flowwidth  = p[4];
00345 
00346 
00347   par.Ptop  = p[8];
00348   par.DPstrip    = p[9];
00349   par.DPrect     = p[9];
00350 
00351 #ifdef PSEUDO_STATES
00352   F_lh      = xd[82];
00353   par.xF    = xd[83];
00354 #else
00355   F_lh      = p[10];
00356   par.xF    = p[11];
00357 #endif
00358 
00359 
00360   par.Tfeed = p[12] +273.15;
00361   par.Tcond = p[13] +273.15;
00362 
00363   /* Let's just put the variables in a nicer form */  
00364 
00365   for (i=0; i<= 39; i++) 
00366     {
00367     par.L[i+1]=xa[i];
00368     par.V[i+1]=xa[40+i]; 
00369     par.Temp[i]=xa[80+i];
00370     }
00371   par.Temp[40]=xa[120];
00372   par.Temp[41]=xa[121];
00373 
00374   for (i=0; i<= 41; i++) 
00375     par.x[i]=xd[i];
00376 
00377   for (i=1; i<= 40; i++) 
00378     par.n[i]=xd[41+i];
00379 
00380 
00381   /* Pressure  */
00382   par.Pres[41]=par.Ptop;
00383   for (i=41; i>=1; i--)
00384     {
00385       if (i>Nfeed)
00386         {par.Pres[i-1]=par.Pres[i]+par.DPrect;}
00387       else
00388         {par.Pres[i-1]=par.Pres[i]+par.DPstrip;}
00389     }
00390 
00391   /* vapor concentrations */
00392   for (i=0; i<= 41; i++)
00393     {
00394     par.P1[i] = exp(A1-B1/(C1+par.Temp[i]));  
00395     par.P2[i] = exp(A2-B2/(C2+par.Temp[i]));  
00396     par.y[i]=par.P1[i]*par.x[i]/par.Pres[i];
00397     }
00398  /* Tray efficiency */
00399   for (i=1; i<= Nfeed; i++)
00400     par.y[i]=par.alphastrip*par.y[i]+(1-par.alphastrip)*par.y[i-1]; 
00401 
00402   for (i=Nfeed+1; i<= 40; i++)
00403     par.y[i]=par.alpharect*par.y[i]+(1-par.alpharect)*par.y[i-1]; 
00404 
00405 
00406   par.y[41]=par.x[41]; /* Condenser outflow out of column is liquid */
00407 
00408   /* recalculate input fluxes from litres/hour into kmol/sec */
00409   par.F     = F_lh / (litre_per_kmol(par.Tfeed, par.xF)*3600); /* [kmol/s] */
00410   Lc     = L_lh / (litre_per_kmol(par.Tcond, par.x[41])*3600); /* [kmol/s] */
00411   /* Lc     = L_lh / (litre_per_kmol(par.Temp[41], par.x[41])*3600); */
00412   /* Lc     = R/(1+R) * par.V[40]; */
00413   par.L[41]=Lc;  
00414 
00415   /* calculate vapour flow */  
00416   i=0;
00417   enthalpies(par.Temp[i], par.Pres[i], par.x[i], par.y[i], &hL,&hV, &dummy);
00418   deltaH=hV-hL;
00419   Vc  = (Q-Qloss)/deltaH;
00420   par.V[0]=Vc;  
00421 
00422   for (i=0; i<= 41; i++) 
00423     { 
00424       if (i==0) holdup=par.nvolB;
00425       else if (i<Nfeed) holdup=par.nvolK;
00426       if (i==Nfeed) holdup=par.nvolM;
00427       else if (i<41) holdup=par.nvolK;
00428       else if (i==41) holdup=par.nvolD;
00429       
00430       par.refvol[i]=holdup;
00431     }
00432   par.n[0] =par.refvol[0]/litre_per_kmol(par.Temp[0], par.x[0]);
00433   par.n[41]=par.refvol[41]/litre_per_kmol(par.Temp[41], par.x[41]);
00434                        
00435     /* Volume of Reboiler is fixed */
00436   par.L[0]=1/litre_per_kmol(par.Temp[0], par.x[0])
00437            * (  litre_per_kmol(par.Temp[0], par.x[1]) * par.L[1]
00438               - litre_per_kmol(par.Temp[0], par.y[0]) * par.V[0]
00439            );
00440 
00441     /* Volume of Condenser is fixed */
00442   par.V[41]=1/litre_per_kmol(par.Temp[41], par.y[41])
00443             * (  litre_per_kmol(par.Temp[41], par.y[40]) * par.V[40]
00444                - litre_per_kmol(par.Temp[41], par.x[41]) * par.L[41]
00445             );
00446 
00447 }
00448 
00449 static void diffeq(
00450                    double *dx,
00451                    double *dn
00452 )
00453      /* 
00454         Calculates dx[42] and dn[40] from
00455         par.x, par.Temp, par.V, par.L, par.y
00456 
00457      */
00458 {
00459 
00460   extern Parameters par;
00461   long i;
00462   double Vmol, dVmol_x;
00463 
00464   /* Reboiler with fixed volume */
00465   i=0;
00466   Vmol=litre_per_kmol(par.Temp[i],  par.x[i]);
00467   dVmol_x=litre_per_kmol(par.Temp[i], 1.0)-litre_per_kmol(par.Temp[i], 0.0);
00468 
00469   dx[i] = 1 / (par.n[i]*(1 + dVmol_x * par.x[i] / Vmol) ) 
00470             * (par.L[i+1]*par.x[i+1] - par.V[i]*par.y[i] - par.L[i]*par.x[i]);
00471 
00472   /* Stripping + Rectifying + Feed */
00473 
00474   for (i=1; i<=40; i++) 
00475     {
00476       dx[i] = -(1/par.n[i])*par.V[i]*(par.y[i]-par.x[i]);
00477       dn[i-1] = -par.V[i]-par.L[i];
00478 
00479       dx[i] += (1/par.n[i])*par.L[i+1]*(par.x[i+1]-par.x[i]);
00480       dn[i-1] += par.L[i+1];
00481 
00482       dx[i] += (1/par.n[i])*par.V[i-1]*(par.y[i-1]-par.x[i]); 
00483       dn[i-1] += par.V[i-1]; 
00484       
00485       if (i==Nfeed) 
00486         {
00487           dx[i] += (1/par.n[i])*par.F*(par.xF-par.x[i]); 
00488           dn[i-1] += par.F;
00489         }
00490     }
00491 
00492   /* Condenser  with fixed volume */
00493   i=41;
00494   Vmol=litre_per_kmol(par.Temp[i],  par.x[i]);
00495   dVmol_x=litre_per_kmol(par.Temp[i], 1.0)-litre_per_kmol(par.Temp[i], 0.0);
00496 
00497   dx[i] = 1 / (par.n[i]*(1 + dVmol_x * par.x[i] / Vmol) ) 
00498             * (par.V[i-1]*par.y[i-1] - par.V[i]*par.y[i] - par.L[i]*par.x[i]);
00499 }
00500 
00501 
00502 static void ffcn( double *t, double *xd, double *xa, double *u, 
00503                   double *p, double *rhs ){
00504 
00505   long i;
00506 #ifdef PSEUDO_STATES
00507   double xdcopy[84];
00508 #else
00509   double xdcopy[82];
00510 #endif
00511   double xacopy[122];
00512   /* scaling */
00513   for(i=0;i<42;i++)
00514     xdcopy[i]=xd[i];
00515   for(i=42;i<82;i++)
00516     xdcopy[i]=xd[i]*1e-3;
00517 #ifdef PSEUDO_STATES
00518    xdcopy[82]=xd[82];
00519    xdcopy[83]=xd[83];
00520 #endif
00521   for(i=0;i<80;i++)
00522     xacopy[i]=xa[i]*1e-5;
00523   for(i=80;i<122;i++)
00524     xacopy[i]=xa[i]+273.15;
00525 
00526   insert(xdcopy, xacopy, u, p);   
00527   diffeq(rhs, &rhs[42]);
00528 
00529   /* scaling */
00530   for(i=42;i<82;i++)
00531     rhs[i]=rhs[i]/1e-3;
00532 
00533 #ifdef PSEUDO_STATES
00534   rhs[82]=0;
00535   rhs[83]=0;
00536 #endif
00537 
00538 }
00539 
00540 
00541 
00542 static void gfcn(double *t, double *xd, double *xa, double *u, 
00543                  double *p, double *rhs ){
00544 
00545   extern Parameters par;
00546 
00547   long i;
00548   double hL_x[42], hL[42], hV[42]; 
00549   double dx[42],dn[42];
00550   double Pfeed, yfeed, P1, P2, dummy;
00551   double hLfeed, hVfeed; 
00552 
00553 #ifdef PSEUDO_STATES
00554   double xdcopy[84];
00555 #else
00556   double xdcopy[82];
00557 #endif
00558   double xacopy[122];
00559 
00560   /* scaling */ 
00561   for(i=0;i<42;i++)
00562     xdcopy[i]=xd[i];
00563   for(i=42;i<82;i++)
00564     xdcopy[i]=xd[i]*1e-3;
00565 #ifdef PSEUDO_STATES
00566    xdcopy[82]=xd[82];
00567    xdcopy[83]=xd[83];
00568 #endif
00569 
00570   for(i=0;i<80;i++)
00571     xacopy[i]=xa[i]*1e-5;
00572   
00573   for(i=80;i<122;i++)
00574     xacopy[i]=xa[i]+273.15;
00575 
00576 
00577   
00578   insert(xdcopy, xacopy, u, p);   
00579   diffeq(dx, dn);
00580 
00581   /* enthalpies */
00582   for (i=0; i<= 41; i++)
00583     {
00584       enthalpies(par.Temp[i], par.Pres[i], par.x[i], par.y[i], &hL[i],&hV[i], &hL_x[i]);
00585     }
00586    
00587   /* entalphy of the feed */
00588  
00589   P1 = exp(A1-B1/(C1+par.Tfeed));
00590   P2 = exp(A2-B2/(C2+par.Tfeed));
00591   Pfeed=P1*par.xF+(1-par.xF)*P2;
00592   yfeed=P1*par.xF/Pfeed;         
00593   
00594   enthalpies(par.Tfeed, Pfeed, par.xF, yfeed, &hLfeed,&hVfeed, &dummy);
00595   
00596   /* enthalphy of the condenser reflux, when Temperature different from equilibrium */
00597       i=41;
00598       par.Temp[i]=par.Tcond;
00599       enthalpies(par.Temp[i], par.Pres[i], par.x[i], par.y[i], &hL[i],&hV[i], &hL_x[i]);
00600   
00601   
00602   /* ++++++++++++++++++++ Algebraic equations ++++++++++++++++++++++ */
00603   
00604   
00605   /* Stripping + Rectifying + Feed - algebraic */
00606   
00607   for (i=1; i<=40; i++)
00608     {
00609       rhs[i-1] = par.L[i]-kmol_per_sec_outflow(par.n[i], par.Temp[i], par.x[i], 
00610                                                par.refvol[i], par.flowwidth);
00611       /*
00612       rhs[40+i-1] =   par.V[i-1]*hV[i-1] 
00613                     - par.V[i]*hV[i] 
00614                     + par.L[i+1]*hL[i+1] 
00615                     - par.L[i]*hL[i] 
00616                     - par.n[i]*hL_x[i]*dx[i]
00617                     - dn[i-1]*hL[i];
00618       if (i==Nfeed) 
00619         {
00620           rhs[40+i-1] += par.F*hLfeed + (1-qF)*hVfeed; 
00621         }
00622       */
00623 
00624       rhs[40+i-1] =   par.V[i-1]*(hV[i-1] - hL[i])  
00625                     - par.V[i]*(hV[i] -hL[i]) 
00626                     + par.L[i+1]*(hL[i+1]-hL[i])  
00627         - par.n[i]*hL_x[i]*dx[i];
00628       if (i==Nfeed) 
00629         {
00630           rhs[40+i-1] += par.F*(hLfeed - hL[i]) + (1-qF)*hVfeed; 
00631         }
00632 
00633     }
00634   
00635   /* Temperatures - algebraic */
00636   for (i=0; i<= 41; i++)
00637     rhs[80+i]=1-par.P1[i]/par.Pres[i]*par.x[i]-(1-par.x[i])*par.P2[i]/par.Pres[i];
00638 
00639   /* scaling */
00640   for(i=0;i<40;i++)
00641     rhs[i]=rhs[i]/1e-5;
00642 }
00643 


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:24