/[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 2679 by jpye, Mon Jan 21 06:23:14 2013 UTC revision 2680 by jpye, Mon Jan 28 06:30:25 2013 UTC
# Line 262  void fprops_triple_point(double *p_t_out Line 262  void fprops_triple_point(double *p_t_out
262    
263  typedef struct{  typedef struct{
264      const PureFluid *P;      const PureFluid *P;
265      double p;      double logp;
266      FpropsError *err;      FpropsError *err;
267      double Terr;      double Terr;
268    #ifdef SAT_DEBUG
269        int neval;
270    #endif
271  } SatPResidData;  } SatPResidData;
272    
273  static ZeroInSubjectFunction sat_p_resid;  static ZeroInSubjectFunction sat_p_resid;
274  static double sat_p_resid(double T, void *user_data){  static double sat_p_resid(double rT, void *user_data){
275  #define D ((SatPResidData *)user_data)  #define D ((SatPResidData *)user_data)
276      double p, rhof, rhog;      double p, rhof, rhog, resid;
277      fprops_sat_T(T, &p, &rhof, &rhog, D->P, D->err);      fprops_sat_T(1/rT, &p, &rhof, &rhog, D->P, D->err);
278      if(*(D->err))D->Terr = T;      if(*(D->err))D->Terr = 1/rT;
279      MSG("T = %f --> p = %f, rhof = %f, rhog = %f, RESID %f", T, p, rhof, rhog, (p - D->p));      resid = log(p) - D->logp;
280        MSG("T = %f --> p = %f, rhof = %f, rhog = %f, RESID %f", 1/rT, p, rhof, rhog, resid);
281      //if(*(D->err))MSG("Error: %s",fprops_error(*(D->err)));      //if(*(D->err))MSG("Error: %s",fprops_error(*(D->err)));
282      //if(*(D->err))return -1;      //if(*(D->err))return -1;
283      return p - D->p;  #ifdef SAT_DEBUG
284        D->neval++;
285    #endif
286        return resid;
287  #undef D  #undef D
288  }  }
289    
# Line 284  static double sat_p_resid(double T, void Line 291  static double sat_p_resid(double T, void
291  /**  /**
292      Solve saturation conditions as a function of pressure.      Solve saturation conditions as a function of pressure.
293    
294      TODO Currently, we will just attempt a Brent solver (zeroin) but hopefully      Currently this is just a Brent solver. We've tried to improve it slightly
295      we can do better later. In particular with cubic EOS this approach seems      by solving for the residual of log(p)-log(p1) as a function of 1/T, which
296      very inefficient. At the very least we should be able to manage a Newton      should make the function a bit more linear.
297      solver...  
298  */        TODO Shouldn't(?) be hard at all to improve this to use a Newton solver via
299        the Clapeyron equation?
300    
301        dp/dT = h_fg/(T*v_fg)
302        
303        where
304    
305        h_fg = h(T, rhog) - h(T, rhof)
306        v_fg = 1/rhog - 1/rhof
307        rhof, rhog are evaluated at T.
308    
309        We guess T, calculate saturation conditions at T, then evaluate dp/dT,
310        use Newton solver to converge, while checking that we remain within
311        Tt < T < Tc. It may be faster to iterate using 1/T as the free variable,
312        and log(p) as the residual variable.
313    
314        TODO Better still would be to reexamine the EOS and find a strategy similar to
315        the Akasaka algorithm that works with with rhof, rhog as the free variables,
316        see helmholtz.c...? Method above uses nested iteration on T inside p, so is
317        going to be slooooooow, even if it's fairly reliable.
318    */
319  void fprops_sat_p(double p, double *T_sat, double *rho_f, double *rho_g, const PureFluid *P, FpropsError *err){  void fprops_sat_p(double p, double *T_sat, double *rho_f, double *rho_g, const PureFluid *P, FpropsError *err){
320      if(*err){      if(*err){
321          MSG("ERROR FLAG ALREADY SET");          MSG("ERROR FLAG ALREADY SET");
# Line 303  void fprops_sat_p(double p, double *T_sa Line 330  void fprops_sat_p(double p, double *T_sa
330      /* FIXME what about checking triple point pressure? */      /* FIXME what about checking triple point pressure? */
331            
332    
333      SatPResidData D = {P, p, err, 0};      SatPResidData D = {
334            P, log(p), err, 0
335    #ifdef SAT_DEBUG
336            ,0
337    #endif
338        };
339      MSG("Solving saturation conditions at p = %f", p);      MSG("Solving saturation conditions at p = %f", p);
340      double p1, T, resid;      double p1, rT, T, resid;
341      int errn;      int errn;
342      double Tt = P->data->T_t;      double Tt = P->data->T_t;
343      if(Tt == 0)Tt = 0.2* P->data->T_c;      if(Tt == 0)Tt = 0.2* P->data->T_c;
344      errn = zeroin_solve(&sat_p_resid, &D, Tt, P->data->T_c, 1e-5, &T, &resid);      errn = zeroin_solve(&sat_p_resid, &D, 1./P->data->T_c, 1./Tt, 1e-10, &rT, &resid);
345      if(*err){      if(*err){
346          MSG("FPROPS error within zeroin_solve iteration ('%s', p = %f, p_c = %f): %s"          MSG("FPROPS error within zeroin_solve iteration ('%s', p = %f, p_c = %f): %s"
347              , P->name, p, P->data->p_c, fprops_error(*err)              , P->name, p, P->data->p_c, fprops_error(*err)
# Line 325  void fprops_sat_p(double p, double *T_sa Line 357  void fprops_sat_p(double p, double *T_sa
357          }          }
358          *err = FPROPS_NO_ERROR;          *err = FPROPS_NO_ERROR;
359      }      }
360        T = 1./rT;
361      fprops_sat_T(T, &p1, rho_f, rho_g, P, err);      fprops_sat_T(T, &p1, rho_f, rho_g, P, err);
362      if(!*err)*T_sat = T;      if(!*err)*T_sat = T;
363      MSG("Got p1 = %f, p = %f", p1, p);      MSG("Got p1 = %f, p = %f in %d iterations", p1, p, D.neval);
364  }  }
365    
366    
Legend:
Removed from v.2679  
changed lines
  Added in v.2680
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22