# Diff of /trunk/models/johnpye/fprops/sat2.c

revision 2224 by jpye, Mon Jul 26 01:14:48 2010 UTC revision 2225 by jpye, Mon Jul 26 03:40:18 2010 UTC
# Line 5  Line 5
5  #include <math.h>  #include <math.h>
6  #include <stdio.h>  #include <stdio.h>
7
8    //#define PHASE_DEBUG
9
10    #ifndef PHASE_DEBUG
11    # define MSG(...)
12    #else
13    # define MSG(ARGS...) fprintf(stderr,"FPROPS: " ARGS)
14    #endif
15
16
17  /**  /**
18      solve Maxwell phase criterion using successive substitution      solve Maxwell phase criterion using successive substitution
19
20      @return p      @return 0 on success, else error code.
21        @param p output, pressure
22      @param rhof output, liquid density      @param rhof output, liquid density
23      @param rhof output, vapour density      @param rhof output, vapour density
@param err error code, returns non-zero if unable to solve.
24  */  */
25  double fprops_sat_succsubs(double T, double *rhof_out, double *rhog_out, const HelmholtzData *d, int *err){  int fprops_sat_T(double T, double *p_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
26      double p, p_new, delta_p, delta_p_old;      double p, p_new, delta_p, delta_p_old;
27
28      /* first guess using acentric factor */      /* first guess using acentric factor */
29      p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));      p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
30
31      fprintf(stderr,"Initial estimate: psat(%f) = %f bar\n",T,p/1e5);      MSG("Initial estimate: psat(%f) = %f bar\n",T,p/1e5);
32
33      int i,j;      int i,j;
34      char use_guess = 0;      char use_guess = 0;
# Line 30  double fprops_sat_succsubs(double T, dou Line 39  double fprops_sat_succsubs(double T, dou
39
40      for(i=0;i<FPROPS_MAX_SUCCSUBS;++i){      for(i=0;i<FPROPS_MAX_SUCCSUBS;++i){
41          if(fprops_rho_pT(p,T,FPROPS_PHASE_LIQUID,use_guess, d, &rhof)){          if(fprops_rho_pT(p,T,FPROPS_PHASE_LIQUID,use_guess, d, &rhof)){
42              fprintf(stderr,"  increasing p, error with rho_f\n");              MSG("  increasing p, error with rho_f\n");
43              p *= 1.005;              p *= 1.005;
44              continue;              continue;
45          }          }
46
47          if(fprops_rho_pT(p,T,FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){          if(fprops_rho_pT(p,T,FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
48              fprintf(stderr,"  decreasing p, error with rho_g\n");              MSG("  decreasing p, error with rho_g\n");
49              p *= 0.95;              p *= 0.95;
50              continue;              continue;
51          }          }
52
53          fprintf(stderr,"Iter %d: p = %f bar, rhof = %f, rhog = %f\n",i, p/1e5, rhof, rhog);          MSG("Iter %d: p = %f bar, rhof = %f, rhog = %f\n",i, p/1e5, rhof, rhog);
54
55          use_guess = 1; /* after first run, start re-using current guess */          use_guess = 1; /* after first run, start re-using current guess */
56
57          if(fabs(rhof - rhog) < 1e-8){          if(fabs(rhof - rhog) < 1e-8){
58              fprintf(stderr,"FPROPS: densities converged to same value\n");              MSG("FPROPS: densities converged to same value\n");
59              *err = 1;              *p_out = p;
60              *rhof_out = rhof;              *rhof_out = rhof;
61              *rhog_out = rhog;              *rhog_out = rhog;
62              return p;              return 1;
63          }          }
64
65          p_new = (helmholtz_a(T, rhof, d) - helmholtz_a(T, rhog, d)) / (1./rhog - 1./rhof);          p_new = (helmholtz_a(T, rhof, d) - helmholtz_a(T, rhog, d)) / (1./rhog - 1./rhof);
66          delta_p = p_new - p;          delta_p = p_new - p;
67
68          //fprintf(stderr," delta_p = %f bar\n",delta_p/1e5);          //MSG(" delta_p = %f bar\n",delta_p/1e5);
69
70          /* convergence test */          /* convergence test */
71          if(abs(delta_p/p) < 1e-6){          if(abs(delta_p/p) < 1e-6){
72              fprintf(stderr,"CONVERGED...\n");              MSG("CONVERGED...\n");
73              /* note possible need for delp non-change test for certain fluids */              /* note possible need for delp non-change test for certain fluids */
74
75              p = p_new;              p = p_new;
76
77              /* find vapour density, using guess */              /* find vapour density, using guess */
78              if(fprops_rho_pT(p, T, FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){              if(fprops_rho_pT(p, T, FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
79                  fprintf(stderr,"FPROPS: fails to estimate vapour density\n");                  MSG("FPROPS: fails to estimate vapour density\n");
*err = 2;
80                  *rhof_out = rhof;                  *rhof_out = rhof;
81                  *rhog_out = rhog;                  *rhog_out = rhog;
82                  return p;                  *p_out = p;
83                    return 2;
84              }              }
85              /* find liquid phase, using guess */              /* find liquid phase, using guess */
86              fprops_rho_pT(p, T, FPROPS_PHASE_LIQUID, use_guess, d, &rhof);              fprops_rho_pT(p, T, FPROPS_PHASE_LIQUID, use_guess, d, &rhof);
# Line 80  double fprops_sat_succsubs(double T, dou Line 89  double fprops_sat_succsubs(double T, dou
89              p = helmholtz_p(T,rhof,d);              p = helmholtz_p(T,rhof,d);
90
91              if(p > d->p_c || rhof < d->rho_c){              if(p > d->p_c || rhof < d->rho_c){
92                  fprintf(stderr,"FPROPS: invalid converged value of p, beyond p_crit\n");                  MSG("FPROPS: invalid converged value of p, beyond p_crit\n");
93                  *err = 3;                  *p_out = p;
94                  *rhof_out = rhof;                  *rhof_out = rhof;
95                  *rhog_out = rhog;                  *rhog_out = rhog;
96                  return p;                  return 3;
97
98              }              }
99
100              fprintf(stderr,"  recalc p = %f bar\n",p/1e5);              MSG("  recalc p = %f bar\n",p/1e5);
101                *p_out = p;
102              *rhof_out = rhof;              *rhof_out = rhof;
103              *rhog_out = rhog;              *rhog_out = rhog;
104              return p;              return 0;
105          }          }
106
107          delta_p_old = delta_p;          delta_p_old = delta_p;
# Line 102  double fprops_sat_succsubs(double T, dou Line 112  double fprops_sat_succsubs(double T, dou
112                  delta_p *= 0.5;                  delta_p *= 0.5;
113                  if(fabs(delta_p) < 0.4 * p)break;                  if(fabs(delta_p) < 0.4 * p)break;
114              }              }
115              fprintf(stderr,"FPROPS: delta_p was too large, has been shortened\n");              MSG("FPROPS: delta_p was too large, has been shortened\n");
116          }          }
117
118          p = p + delta_p;          p = p + delta_p;
119      }      }
120
121      if(i = FPROPS_MAX_SUCCSUBS){      MSG("FPROPS: too many iterations in %s, or need to try alternative method\n",__func__);
122          fprintf(stderr,"FPROPS: too many iterations in %s, need to try alternative method\n",__func__);      *p_out = p;
123          *err = 1;      *rhof_out = rhof;
124          *rhof_out = rhof;      *rhog_out = rhog;
125          *rhog_out = rhog;      return 1;
return p;
}
126  }  }
127
128
# Line 234  int fprops_rho_pT(double p, double T, ch Line 242  int fprops_rho_pT(double p, double T, ch
242                      }                      }
243                      if (vlog < -5.0 && T > d->T_c){                      if (vlog < -5.0 && T > d->T_c){
244                          use_liquid_calc = 0;                          use_liquid_calc = 0;
245                          fprintf(stderr,"FPROPS: switching to vapour calc\n");                          MSG("FPROPS: switching to vapour calc\n");
246                          break; /* switch to vapour calc */                          break; /* switch to vapour calc */
247                      }                      }
248                  }                  }
# Line 249  int fprops_rho_pT(double p, double T, ch Line 257  int fprops_rho_pT(double p, double T, ch
257          }          }
258
259          if(use_liquid_calc){          if(use_liquid_calc){
260              fprintf(stderr,"FPROPS: liquid iteration did not converge\n");              MSG("FPROPS: liquid iteration did not converge\n");
261              return 1;              return 1;
262          }          }
263      }      }
# Line 288  int fprops_rho_pT(double p, double T, ch Line 296  int fprops_rho_pT(double p, double T, ch
296
297      // not converged      // not converged
298      *rho=1./exp(vlog);      *rho=1./exp(vlog);
299      fprintf(stderr,"FPROPS: liquid iteration did not converge\n");      MSG("FPROPS: liquid iteration did not converge\n");
300      return 1;      return 1;
301  }  }
302

Legend:
 Removed from v.2224 changed lines Added in v.2225