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

Diff of /branches/sid/models/johnpye/fprops/sat.c

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

revision 2287 by jpye, Sun Aug 15 01:19:21 2010 UTC revision 2288 by jpye, Sun Aug 15 06:16:13 2010 UTC
# Line 48  Line 48 
48  #include <fenv.h>  #include <fenv.h>
49  #endif  #endif
50    
51    #ifdef SAT_DEBUG
52    # define MSG(STR,...) fprintf(stderr,"%s:%d: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
53    #else
54    # define MSG(ARGS...)
55    #endif
56    
57  /**  /**
58      Estimate of saturation pressure using H W Xiang ''The new simple extended      Estimate of saturation pressure using H W Xiang ''The new simple extended
59      corresponding-states principle: vapor pressure and second virial      corresponding-states principle: vapor pressure and second virial
# Line 273  double fprops_pc(const HelmholtzData *d) Line 279  double fprops_pc(const HelmholtzData *d)
279      Calculate the critical pressure using the T_c and rho_c values in the HelmholtzData.      Calculate the critical pressure using the T_c and rho_c values in the HelmholtzData.
280  */  */
281  int fprops_triple_point(double *p_t_out, double *rhof_t_out, double *rhog_t_out, const HelmholtzData *d){  int fprops_triple_point(double *p_t_out, double *rhof_t_out, double *rhog_t_out, const HelmholtzData *d){
282        MSG("Startin triple point routine");
283      static const HelmholtzData *d_last = NULL;      static const HelmholtzData *d_last = NULL;
284      static double p_t, rhof_t, rhog_t;      static double p_t, rhof_t, rhog_t;
285      if(d == d_last){      if(d == d_last){
# Line 281  int fprops_triple_point(double *p_t_out, Line 288  int fprops_triple_point(double *p_t_out,
288          *rhog_t_out = rhog_t;          *rhog_t_out = rhog_t;
289          return 0;          return 0;
290      }      }
291        MSG("Calculating saturation for T = %f",d->T_t);
292      int res = fprops_sat_T(d->T_t, &p_t, &rhof_t, &rhog_t,d);      int res = fprops_sat_T(d->T_t, &p_t, &rhof_t, &rhog_t,d);
293      if(res)return res;      if(res)return res;
294      else{      else{
# Line 337  int fprops_sat_p(double p, double *Tsat_ Line 345  int fprops_sat_p(double p, double *Tsat_
345          double dpdT_sat = (hg - hf) / T1 / (1./rhog - 1./rhof);          double dpdT_sat = (hg - hf) / T1 / (1./rhog - 1./rhof);
346          //fprintf(stderr,"\t\tdpdT_sat = %f bar/K\n",dpdT_sat/1e5);          //fprintf(stderr,"\t\tdpdT_sat = %f bar/K\n",dpdT_sat/1e5);
347          double delta_T = -(p1 - p)/dpdT_sat;          double delta_T = -(p1 - p)/dpdT_sat;
348          if(T1 + delta_T < d->T_t)T1 = 0.5 * (d->T_t + T1);          if(T1 + delta_T < d->T_t - 1e-4)T1 = 0.5 * (d->T_t + T1);
349          else T1 += delta_T;          else T1 += delta_T;
350      }      }
351      *Tsat_out = T1;      *Tsat_out = T1;

Legend:
Removed from v.2287  
changed lines
  Added in v.2288

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22