/[ascend]/trunk/models/johnpye/fprops/pengrob.c
ViewVC logotype

Diff of /trunk/models/johnpye/fprops/pengrob.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2663 by jpye, Fri Jan 18 00:34:51 2013 UTC revision 2664 by jpye, Fri Jan 18 06:02:15 2013 UTC
# Line 48  PropEvalFn pengrob_dpdrho_T; Line 48  PropEvalFn pengrob_dpdrho_T;
48  PropEvalFn pengrob_alphap;  PropEvalFn pengrob_alphap;
49  PropEvalFn pengrob_betap;  PropEvalFn pengrob_betap;
50  SatEvalFn pengrob_sat;  SatEvalFn pengrob_sat;
51    //SatEvalFn pengrob_sat_akasaka;
52    
53  static double MidpointPressureCubic(double T, const FluidData *data, FpropsError *err);  static double MidpointPressureCubic(double T, const FluidData *data, FpropsError *err);
54    
55  //#define PR_DEBUG  #define PR_DEBUG
56  #define PR_ERRORS  #define PR_ERRORS
57    
58  #ifdef PR_DEBUG  #ifdef PR_DEBUG
# Line 185  PureFluid *pengrob_prepare(const EosData Line 186  PureFluid *pengrob_prepare(const EosData
186  #undef I  #undef I
187  #undef D  #undef D
188  #undef C  #undef C
189        //P->sat_fn = &pengrob_sat_akasaka;
190    
191      return P;      return P;
192  }  }
193    
# Line 245  double pengrob_h(double T, double rho, c Line 248  double pengrob_h(double T, double rho, c
248      DEFINE_SQRTALPHA;      DEFINE_SQRTALPHA;
249      DEFINE_A;      DEFINE_A;
250      DEFINE_V;      DEFINE_V;
251        if(rho > 1./PD->b){
252            MSG("Density exceeds limit value 1/b = %f",1./PD->b);
253            *err = FPROPS_RANGE_ERROR;
254            return 0;
255        }
256      double h0 = ideal_h(T,rho,data,err);      double h0 = ideal_h(T,rho,data,err);
257      double p = pengrob_p(T, rho, data, err);      double p = pengrob_p(T, rho, data, err);
258      double Z = p * v / (data->R * T);      double Z = p * v / (data->R * T);
# Line 278  double pengrob_s(double T, double rho, c Line 286  double pengrob_s(double T, double rho, c
286    
287  double pengrob_a(double T, double rho, const FluidData *data, FpropsError *err){  double pengrob_a(double T, double rho, const FluidData *data, FpropsError *err){
288      // FIXME maybe we can improve this with more direct maths      // FIXME maybe we can improve this with more direct maths
     double b = PD->b;  
     if(rho > 1./b){  
         MSG("Density exceeds limit value 1/b = %f",1./b);  
         *err = FPROPS_RANGE_ERROR;  
     }  
289      double h = pengrob_h(T,rho,data,err);      double h = pengrob_h(T,rho,data,err);
290      double s = pengrob_s(T,rho,data,err); // duplicated calculation of p!      double s = pengrob_s(T,rho,data,err); // duplicated calculation of p!
291      double p = pengrob_p(T,rho,data,err); // duplicated calculation of p!      double p = pengrob_p(T,rho,data,err); // duplicated calculation of p!
# Line 303  double pengrob_u(double T, double rho, c Line 306  double pengrob_u(double T, double rho, c
306  }  }
307    
308  double pengrob_g(double T, double rho, const FluidData *data, FpropsError *err){  double pengrob_g(double T, double rho, const FluidData *data, FpropsError *err){
309        if(rho > 1./PD->b){
310            MSG("Density exceeds limit value 1/b = %f",1./PD->b);
311            *err = FPROPS_RANGE_ERROR;
312        }
313    #if 0
314      double h = pengrob_h(T,rho,data,err);      double h = pengrob_h(T,rho,data,err);
315      double s = pengrob_s(T,rho,data,err); // duplicated calculation of p!      double s = pengrob_s(T,rho,data,err); // duplicated calculation of p!
316        if(isnan(h))MSG("h is nan");
317        if(isnan(s))MSG("s is nan");
318      return h - T*s;      return h - T*s;
319  //  previous code from Richard, probably fine but need to check  #else
320  //  DEFINE_A;      //  previous code from Richard, probably fine but need to check
321  //  DEFINE_V;      DEFINE_SQRTALPHA;
322  //  double p = pengrob_p(T, rho, data, err);      DEFINE_A;
323  //  double Z = p*v/(data->R * T);      DEFINE_V;
324  //  double B = p*PD->b/(data->R * T);      double p = pengrob_p(T, rho, data, err);
325  //  double A = p * a / SQ(data->R * T);      double Z = p*v/(data->R * T);
326  //  return -log(fabs(Z-B))-(A/(sqrt(8)*B))*log(fabs((Z+(1+sqrt(2))*B)/(Z+(1-sqrt(2))*B)))+Z-1;      double B = p*PD->b/(data->R * T);
327        double A = p * a / SQ(data->R * T);
328        return log(fabs(Z-B))-(A/(sqrt(8)*B))*log(fabs((Z+(1+sqrt(2))*B)/(Z+(1-sqrt(2))*B)))+Z-1;
329    #endif
330  }  }
331    
332  /**  /**
# Line 492  double pengrob_sat(double T,double *rhof Line 505  double pengrob_sat(double T,double *rhof
505      double sqrt2 = sqrt(2);      double sqrt2 = sqrt(2);
506    
507      double p = fprops_psat_T_acentric(T, data);      double p = fprops_psat_T_acentric(T, data);
508        MSG("Initial guess: p = %f from acentric factor",p);
509    
510      int i = 0;      int i = 0;
511      double Zg, Z1, Zf, vg, vf;      double Zg, Z1, Zf, vg, vf;
512    
513        FILE *F1 = fopen("pf.txt","w");
514    
515      // FIXME test upper iteration limit required      // FIXME test upper iteration limit required
516      while(++i < 100){      while(++i < 300){
517          MSG("iter %d: p = %f, rhof = %f, rhog = %f", i, p, 1/vf, 1/vg);          MSG("iter %d: p = %f, rhof = %f, rhog = %f", i, p, 1/vf, 1/vg);
518          // Peng & Robinson eq 17          // Peng & Robinson eq 17
519          double sqrtalpha = 1 + PD->kappa * (1 - sqrt(T / PD_TCRIT));          double sqrtalpha = 1 + PD->kappa * (1 - sqrt(T / PD_TCRIT));
# Line 529  double pengrob_sat(double T,double *rhof Line 546  double pengrob_sat(double T,double *rhof
546              double ff = FUG(Zf,A,B);              double ff = FUG(Zf,A,B);
547              double fg = FUG(Zg,A,B);              double fg = FUG(Zg,A,B);
548              double fratio = ff/fg;              double fratio = ff/fg;
549              //MSG("    ff = %f, fg = %f, fratio = %f", ff, fg, fratio);              MSG("    ff = %f, fg = %f, fratio = %f", ff, fg, fratio);
550    
551              //double hf = pengrob_h(T, 1/vf, data, err);              //double hf = pengrob_h(T, 1/vf, data, err);
552              //double hg = pengrob_h(T, 1/vg, data, err);              //double hg = pengrob_h(T, 1/vg, data, err);
# Line 540  double pengrob_sat(double T,double *rhof Line 557  double pengrob_sat(double T,double *rhof
557                  *rhog_ret = 1 / vg;                  *rhog_ret = 1 / vg;
558                  p = pengrob_p(T, *rhog_ret, data, err);                  p = pengrob_p(T, *rhog_ret, data, err);
559                  MSG("Solved for T = %f: p = %f, rhof = %f, rhog = %f", T, p, *rhof_ret, *rhog_ret);                  MSG("Solved for T = %f: p = %f, rhof = %f, rhog = %f", T, p, *rhof_ret, *rhog_ret);
560                    fclose(F1);
561                  return p;                  return p;
562              }              }
563                fprintf(F1,"%f\t%f\n",p,fratio);
564              p *= fratio;              p *= fratio;
565          }else{          }else{
566              /* In this case we need to adjust our guess p(T) such that we get              /* In this case we need to adjust our guess p(T) such that we get
# Line 549  double pengrob_sat(double T,double *rhof Line 568  double pengrob_sat(double T,double *rhof
568              p = MidpointPressureCubic(T, data, err);              p = MidpointPressureCubic(T, data, err);
569              if(*err){              if(*err){
570                  ERRMSG("Failed to solve for a midpoint pressure");                  ERRMSG("Failed to solve for a midpoint pressure");
571                    fclose(F1);
572                  return p;                  return p;
573              }              }
574              MSG("    single root: Z = %f. new pressure guess: %f", Zf, p);              MSG("    single root: Z = %f. new pressure guess: %f", Zf, p);
# Line 556  double pengrob_sat(double T,double *rhof Line 576  double pengrob_sat(double T,double *rhof
576      }      }
577      MSG("Did not converge");      MSG("Did not converge");
578      *err = FPROPS_SAT_CVGC_ERROR;      *err = FPROPS_SAT_CVGC_ERROR;
579        fclose(F1);
580      return 0;      return 0;
581  }  }
582    
583  /*  /*
584      FIXME can we generalise this to work with other cubic EOS as well?      FIXME can we generalise this to work with other *cubic* EOS as well?
585      Currently not easily done since pointers are kept at a higher level in our      Currently not easily done since pointers are kept at a higher level in our
586      data structures than FluidData.      data structures than FluidData.
587  */  */
# Line 621  double MidpointPressureCubic(double T, c Line 642  double MidpointPressureCubic(double T, c
642      return 0.5*(p1 + p2);      return 0.5*(p1 + p2);
643  }  }
644    
   
645  static double resid_dpdrho_T(double rho, void *user_data){  static double resid_dpdrho_T(double rho, void *user_data){
646  #define D ((MidpointSolveData *)user_data)  #define D ((MidpointSolveData *)user_data)
647      return pengrob_dpdrho_T(D->T,rho,D->data,D->err);      return pengrob_dpdrho_T(D->T,rho,D->data,D->err);
648  #undef D  #undef D
649  }  }
650    
651    
652    #if 0
653    
654    /* we would like to try the Akasaka approach and use it with cubic EOS.
655    but at this stage it's not working correctly. Not sure why... */
656    
657    /**
658        Solve saturation condition for a specified temperature using approach of
659        Akasaka, but adapted for general use to non-helmholtz property correlations.
660        @param T temperature [K]
661        @param psat_out output, saturation pressure [Pa]
662        @param rhof_out output, saturated liquid density [kg/m^3]
663        @param rhog_out output, saturated vapour density [kg/m^3]
664        @param d helmholtz data object for the fluid in question.
665        @return 0 on success, non-zero on error (eg algorithm failed to converge, T out of range, etc.)
666    */
667    double pengrob_sat_akasaka(double T, double *rhof_out, double * rhog_out, const FluidData *data, FpropsError *err){
668        if(T < data->T_t - 1e-8){
669            ERRMSG("Input temperature %f K is below triple-point temperature %f K",T,data->T_t);
670            return FPROPS_RANGE_ERROR;
671        }
672    
673        if(T > data->T_c){
674            ERRMSG("Input temperature is above critical point temperature");
675            return FPROPS_RANGE_ERROR;
676        }
677    
678        // we're at the critical point
679        if(fabs(T - data->T_c) < 1e-9){
680            *rhof_out = data->rho_c;
681            *rhog_out = data->rho_c;
682            return data->p_c;
683        }
684    
685        // FIXME at present step-length multiplier is set to 0.4 just because of
686        // ONE FLUID, ethanol. Probably our initial guess data isn't good enough,
687        // or maybe there's a problem with the acentric factor or something like
688        // that. This factor 0.4 will be slowing down the whole system, so it's not
689        // good. TODO XXX.
690    
691        // initial guesses for liquid and vapour density
692        double rhof = fprops_rhof_T_rackett(T,data);
693        double rhog= fprops_rhog_T_chouaieb(T,data);
694        double R = data->R;
695        double pc = data->p_c;
696    
697        MSG("initial guess rho_f = %f, rho_g = %f\n",rhof,rhog);
698        MSG("calculating for T = %.12e",T);
699    
700        int i = 0;
701        while(i++ < 70){
702            if(rhof > 1/PD->b){
703                MSG("limit rhof to 1/b");
704                rhof = 1/PD->b;
705            }
706    
707            MSG("iter %d: T = %f, rhof = %f, rhog = %f",i,T, rhof, rhog);
708    
709            double pf = pengrob_p(T,rhof,data,err);
710            double pg = pengrob_p(T,rhog,data,err);
711            double gf = pengrob_g(T,rhof,data,err);
712            double gg = pengrob_g(T,rhog,data,err);
713            double dpdrf = pengrob_dpdrho_T(T,rhof,data,err);
714            double dpdrg = pengrob_dpdrho_T(T,rhog,data,err);
715            if(*err){
716                ERRMSG("error returned");
717            }
718    
719            MSG("gf = %f, gg = %f", gf, gg);
720    
721            // jacobian for [F;G](rhof, rhog) --- derivatives wrt rhof and rhog
722            double F = (pf - pg)/pc;
723            double G = (gf - gg)/R/T;
724    
725            if(fabs(F) + fabs(G) < 1e-12){
726                //fprintf(stderr,"%s: CONVERGED\n",__func__);
727                *rhof_out = rhof;
728                *rhog_out = rhog;
729                return pengrob_p(T, *rhog_out, data, err);
730                /* SUCCESS */
731            }
732    
733            double Ff = dpdrf/pc;
734            double Fg = -dpdrg/pc;
735            MSG("Ff = %e, Fg = %e",Ff,Fg);
736    
737            double Gf = dpdrf/rhof/R/T;
738            double Gg = -dpdrg/rhog/R/T;
739            MSG("Gf = %e, Gg = %e",Gf,Gg);
740    
741            double DET = Ff*Gg - Fg*Gf;
742            MSG("DET = %f",DET);
743    
744            MSG("F = %f, G = %f", F, G);
745    
746            // 'gamma' needs to be increased to 0.5 for water to solve correctly (see 'test/sat.c')
747    #define gamma 1
748            rhof += gamma/DET * (Fg*G - Gg*F);
749            rhog += gamma/DET * ( Gf*F - Ff*G);
750    #undef gamma
751    
752            if(rhog < 0)rhog = -0.5*rhog;
753            if(rhof < 0)rhof = -0.5*rhof;
754        }
755        *rhof_out = rhof;
756        *rhog_out = rhog;
757        *err = FPROPS_SAT_CVGC_ERROR;
758        ERRMSG("Not converged: with T = %e (rhof=%f, rhog=%f).",T,*rhof_out,*rhog_out);
759        return pengrob_p(T, rhog, data, err);
760    }
761    
762    #endif
763    
764    
Legend:
Removed from v.2663  
changed lines
  Added in v.2664
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22