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

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

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

revision 2289 by jpye, Sun Aug 15 06:34:28 2010 UTC revision 2290 by jpye, Sun Aug 15 10:11:54 2010 UTC
# Line 46  Line 46 
46  #ifdef THROW_FPE  #ifdef THROW_FPE
47  #define _GNU_SOURCE  #define _GNU_SOURCE
48  #include <fenv.h>  #include <fenv.h>
49    int feenableexcept (int excepts);
50    int fedisableexcept (int excepts);
51    int fegetexcept (void);
52  #endif  #endif
53    
54  #ifdef SAT_DEBUG  #ifdef SAT_DEBUG
# Line 53  Line 56 
56  #else  #else
57  # define MSG(ARGS...)  # define MSG(ARGS...)
58  #endif  #endif
59    #define ERRMSG(STR,...) fprintf(stderr,"%s:%d: ERROR: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
60    
61  /**  /**
62      Estimate of saturation pressure using H W Xiang ''The new simple extended      Estimate of saturation pressure using H W Xiang ''The new simple extended
# Line 373  int fprops_sat_p(double p, double *Tsat_ Line 377  int fprops_sat_p(double p, double *Tsat_
377  /**  /**
378      Calculate Tsat based on a value of hf. This value is useful in setting      Calculate Tsat based on a value of hf. This value is useful in setting
379      first guess Temperatures when solving for the coordinates (p,h).      first guess Temperatures when solving for the coordinates (p,h).
380      Secant method.      This function uses the secant method for the iterative solution.
381  */  */
382  int fprops_sat_hf(double hf, double *Tsat_out, double *psat_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){  int fprops_sat_hf(double hf, double *Tsat_out, double *psat_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
383      double T1 = 0.4 * d->T_t + 0.6 * d->T_c;      double T1 = 0.4 * d->T_t + 0.6 * d->T_c;
# Line 381  int fprops_sat_hf(double hf, double *Tsa Line 385  int fprops_sat_hf(double hf, double *Tsa
385      double h1, h2, p, rhof, rhog;      double h1, h2, p, rhof, rhog;
386      int res = fprops_sat_T(T2, &p, &rhof, &rhog, d);      int res = fprops_sat_T(T2, &p, &rhof, &rhog, d);
387      if(res){      if(res){
388          fprintf(stderr,"%s:%d: Failed to solve psat(T_t)\n",__func__,__LINE__);          ERRMSG("Failed to solve psat(T_t = %.12e) for %s",T2,d->name);
389          return 1;          return 1;
390      }      }
391      double tol = 1e-6;      double tol = 1e-6;
392      h2 = helmholtz_h(T2,rhof,d);      h2 = helmholtz_h(T2,rhof,d);
393        if(hf < h2){
394            ERRMSG("Value given for hf = %.12e is below that calculated for triple point liquid hf_t = %.12e",hf,h2);
395            return 2;
396        }
397    
398      int i = 0;      int i = 0;
399      while(i++ < 40){      while(i++ < 60){
400          assert(T1 >= d->T_t);          assert(T1 >= d->T_t - 1e-4);
401          assert(T1 <= d->T_c);          assert(T1 <= d->T_c);
402          res = fprops_sat_T(T1, &p, &rhof, &rhog, d);          res = fprops_sat_T(T1, &p, &rhof, &rhog, d);
403          if(res){          if(res){
404              fprintf(stderr,"%s:%d: Failed to solve psat(T = %.12e)\n",__func__,__LINE__,T1);              ERRMSG("Failed to solve psat(T = %.12e) for %s",T1,d->name);
405              return 1;              return 1;
406          }          }
407          h1 = helmholtz_h(T1,rhof, d);          h1 = helmholtz_h(T1,rhof, d);
# Line 403  int fprops_sat_hf(double hf, double *Tsa Line 412  int fprops_sat_hf(double hf, double *Tsa
412              *rhog_out = rhog;              *rhog_out = rhog;
413              return 0;              return 0;
414          }          }
415            if(h1 == h2){
416                MSG("With %s, got h1 = h2 = %.12e, but hf = %.12e!",d->name,h1,hf);
417                return 2;
418            }
419    
420          double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);          double delta_T = -(h1 - hf) * (T1 - T2) / (h1 - h2);
421          T2 = T1;          T2 = T1;
422          h2 = h1;          h2 = h1;
423          while(T1 + delta_T > d->T_c)delta_T *= 0.5;          while(T1 + delta_T > d->T_c)delta_T *= 0.5;
         while(T1 + delta_T < d->T_t)delta_T *= 0.5;  
424          T1 += delta_T;          T1 += delta_T;
425            if(T1 < d->T_t)T1 = d->T_t;
426          if(i==20 || i==30)tol*=100;          if(i==20 || i==30)tol*=100;
427      }      }
428      fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);      fprintf(stderr,"Failed to solve Tsat for hf = %f (got to T = %f)\n",hf,T1);
Legend:
Removed from v.2289  
changed lines
  Added in v.2290
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22