/[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 2261 by jpye, Tue Jul 27 01:56:55 2010 UTC revision 2262 by jpye, Thu Aug 5 08:34:43 2010 UTC
# Line 26  Line 26 
26  */  */
27    
28  #include "sat.h"  #include "sat.h"
   
29  #include "helmholtz_impl.h"  #include "helmholtz_impl.h"
 # include <assert.h>  
30    
31    # include <assert.h>
32  #include <math.h>  #include <math.h>
33  #ifdef TEST  #include <stdio.h>
 # include <stdio.h>  
 # include <stdlib.h>  
 #endif  
   
 #include <gsl/gsl_multiroots.h>  
34    
35  #define SQ(X) ((X)*(X))  #define SQ(X) ((X)*(X))
36    
37    #if 0
38    #define _GNU_SOURCE
39    #include <fenv.h>
40    #endif
41    
42  /**  /**
43      Estimate of saturation pressure using H W Xiang ''The new simple extended      Estimate of saturation pressure using H W Xiang ''The new simple extended
44      corresponding-states principle: vapor pressure and second virial      corresponding-states principle: vapor pressure and second virial
# Line 145  double fprops_rhog_T_chouaieb(double T, Line 144  double fprops_rhog_T_chouaieb(double T,
144      return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));      return D->rho_c * exp(PPP * (pow(alpha,NNN) - exp(1-alpha)));
145  }  }
146    
147    int fprops_sat_T(double T, double *psat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
148  /**      double tau = d->T_c / T;
149      Maxwell phase criterion, first principles      double delf = 1.1 * fprops_rhof_T_rackett(T,d) / d->rho_c;
150  */      double delg = 0.9 * fprops_rhog_T_chouaieb(T,d) / d->rho_c;
151  void phase_criterion(double T, double rho_f, double rho_g, double p_sat, double *eq1, double *eq2, double *eq3, const HelmholtzData *D){      //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
152  #ifdef TEST  
153      fprintf(stderr,"PHASE CRITERION: T = %f, rho_f = %f, rho_g = %f, p_sat = %f\n", T, rho_f, rho_g, p_sat);      int i = 0;
154  #endif      while(i++ < 50){
155      double delta_f, delta_g, tau;          //fprintf(stderr,"%s: iter %d: rhof = %f, rhog = %f\n",__func__,i,delf*d->rho_c, delg*d->rho_c);
156      tau = D->T_c / T;          double phirf = helm_resid(tau,delf,d);
157            double phirf_d = helm_resid_del(tau,delf,d);
158      delta_f = rho_f / D->rho_c;          double phirf_dd = helm_resid_deldel(tau,delf,d);
159      delta_g = rho_g/ D->rho_c;          double phirg = helm_resid(tau,delg,d);
160            double phirg_d = helm_resid_del(tau,delg,d);
161  #ifdef TEST          double phirg_dd = helm_resid_deldel(tau,delg,d);
162      assert(!isnan(delta_f));  
163      assert(!isnan(delta_g));  #define J(FG) (del##FG * (1. + del##FG * phir##FG##_d))
164      assert(!isnan(p_sat));        #define K(FG) (del##FG * phir##FG##_d + phir##FG + log(del##FG))
165      assert(!isinf(delta_f));  #define J_del(FG) (1 + 2 * del##FG * phir##FG##_d + SQ(del##FG) * phir##FG##_dd)
166      assert(!isinf(delta_g));  #define K_del(FG) (2 * phir##FG##_d + del##FG * phir##FG##_dd + 1./del##FG)
167      assert(!isinf(p_sat));                double Jf = J(f);
168  #endif          double Jg = J(g);
169            double Kf = K(f);
170      *eq1 = (p_sat - helmholtz_p(T, rho_f, D));          double Kg = K(g);
171      *eq2 = (p_sat - helmholtz_p(T, rho_g, D));          double Jf_del = J_del(f);
172      *eq3 = helmholtz_g(T, rho_f, D) - helmholtz_g(T,rho_g, D);          double Jg_del = J_del(g);
173            double Kf_del = K_del(f);
174  #ifdef TEST          double Kg_del = K_del(g);
175      fprintf(stderr,"eq1 = %e\t\teq2 = %e\t\teq3 = %e\n", *eq1, *eq2, *eq3);  
176  #endif          double DELTA = Jg_del * Kf_del - Jf_del * Kg_del;
177  }  
178    #define gamma 1
179            delf += gamma/DELTA * ((Kg - Kf) * Jg_del - (Jg - Jf) * Kg_del);
180  typedef struct PhaseSolve_struct{          delg += gamma/DELTA * ((Kg - Kf) * Jf_del - (Jg - Jf) * Kf_del);
181      double tau;  
182      const HelmholtzData *D;          if(fabs(Kg - Kf) + fabs(Jg - Jf) < 1e-8){
183      const double *scaling;              //fprintf(stderr,"%s: CONVERGED\n",__func__);
184  } PhaseSolve;              *rhof_out = delf * d->rho_c;
185                *rhog_out = delg * d->rho_c;
186  /**              *psat_out = helmholtz_p_raw(T, *rhog_out, d);
187      Maxwell phase criterion, from equations from IAPWS-95.              return 0;
188            }
     We define here pi = p_sat / (R T rho_c)  
 */  
 static int phase_resid(const gsl_vector *x, void *params, gsl_vector *f){  
     const double tau = ((PhaseSolve *)params)->tau;  
     const double pi = gsl_vector_get(x,0);  
     double delta_f = gsl_vector_get(x,1);  
     double delta_g = exp(gsl_vector_get(x,2));  
     const HelmholtzData *D = ((PhaseSolve *)params)->D;  
     const double *scaling = ((PhaseSolve *)params)->scaling;  
   
     assert(!isinf(pi));  
     assert(!isinf(delta_f));  
     assert(!isinf(delta_g));  
   
     gsl_vector_set(f, 0, scaling[0] * (delta_f + delta_f * delta_f * helm_resid_del(tau, delta_f, D) - pi));  
     gsl_vector_set(f, 1, scaling[1] * (1 + delta_g * helm_resid_del(tau, delta_g, D) - pi/delta_g));  
     gsl_vector_set(f, 2, scaling[2] * (pi * (delta_f - delta_g) - delta_f * delta_g * (log(delta_f / delta_g)  + helm_resid(delta_g,tau,D) - helm_resid(delta_f,tau,D))));  
   
   
     if(isnan(gsl_vector_get(f,2)) || isnan(gsl_vector_get(f,2)) || isnan(gsl_vector_get(f,2))){  
         fprintf(stderr,"NaN encountered... ");  
189      }      }
190        *rhof_out = delf * d->rho_c;
191        *rhog_out = delg * d->rho_c;
192        *psat_out = helmholtz_p_raw(T, *rhog_out, d);
193        fprintf(stderr,"%s: NOT CONVERGED for '%s' with T = %e (rhof=%f, rhog=%f)\n",__func__,d->name,T,*rhof_out,*rhog_out);
194        return 1;
195    
     //gsl_vector_set(f, 2, pi * (delta_f - delta_g) * delta_f*delta_g* (log(delta_f/delta_g) + helm_resid(delta_g,tau,D) - helm_resid(delta_f,tau,D)));  
   
     return GSL_SUCCESS;  
196  }  }
197    
198  static int print_state(size_t iter, gsl_multiroot_fsolver * s){  int fprops_sat_p(double p, double *Tsat_out, double *rhof_out, double * rhog_out, const HelmholtzData *d){
199      fprintf(stderr,"iter = %3u: pi = %.3f delf = %.3f delg = %.3f  "      fprintf(stderr,"%s: NOT YET IMPLEMENTED\n",__func__);
200          "E = % .3e % .3e %.3e\n",      return 1;
         iter,  
         gsl_vector_get(s->x, 0),  
         gsl_vector_get(s->x, 1),  
         gsl_vector_get(s->x, 2),  
         gsl_vector_get (s->f, 0),  
         gsl_vector_get (s->f, 1),  
         gsl_vector_get (s->f, 2));  
201  }  }
202    
 #if 0  
 int phase_solve(double T, double *p_sat, double *rho_f, double *rho_g, const HelmholtzData *D){  
     gsl_multiroot_fsolver *s;  
     int status;  
     const size_t n = 3;  
     double tau = D->T_c / T;  
     size_t i, iter = 0;  
     const gsl_multiroot_fsolver_type *ftype;  
   
     const double scaling[3] = {1., 1., 0.1};  
   
     PhaseSolve p = {  
         tau  
         , D  
         , scaling  
     };  
   
     gsl_multiroot_function f = {  
         &phase_resid,  
         n, &p  
     };  
   
     /* TODO use the first guess 'provided' if psat, rho_f, rho_g less than zero? */  
   
     double x_init[3] = {  
         /* first guesses, such as they are... */  
         fprops_psat_T_xiang(T, D) / D->R / T / D->rho_c  
         ,fprops_rhof_T_rackett(T, D) / D->rho_c  
         ,log(fprops_rhog_T_chouaieb(T, D) / D->rho_c)  
     };  
   
     gsl_vector *x = gsl_vector_alloc (n);  
     for(i=0; i<3; ++i)gsl_vector_set(x, i, x_init[i]);  
     ftype = gsl_multiroot_fsolver_broyden;  
     s = gsl_multiroot_fsolver_alloc(ftype, n);  
     gsl_multiroot_fsolver_set(s, &f, x);  
   
     //print_state(iter, s);  
   
     double pi_xiang = fprops_psat_T_xiang(T, D) / D->R / T / D->rho_c;  
     double delg_chouaieb = fprops_rhog_T_chouaieb(T, D) / D->rho_c;  
     double delf_rackett = fprops_rhof_T_rackett(T, D) / D->rho_c;  
   
     assert(!isnan(pi_xiang));  
     assert(!isnan(delg_chouaieb));  
     assert(!isnan(delf_rackett));  
   
     do{  
         iter++;  
   
         status = gsl_multiroot_fsolver_iterate(s);  
   
         print_state(iter, s);  
         if(status)break;  
   
         status = gsl_multiroot_test_residual(s->f, 1e-7);  
   
         double pi = gsl_vector_get(s->x,0);  
         double delf = gsl_vector_get(s->x,1);  
         double delg = exp(gsl_vector_get(s->x,2));  
   
   
 #define VAR_PI 0  
 #define VAR_DELF 1  
 #define VAR_DELG 2  
 #define CHECK_RESET(COND, VAR, NEWVAL)\  
     if(COND){\  
         gsl_vector_set(s->x,VAR,NEWVAL);\  
         fprintf(stderr,"RESET %s to %s = %f (FAILED %s)\n", #VAR, #NEWVAL, NEWVAL, #COND);\  
     }  
   
 #define CHECK_RESET_CONTINUE(COND, VAR, NEWVAL)\  
     if(COND){\  
         gsl_vector_set(s->x,VAR,NEWVAL);\  
         fprintf(stderr,"RESET %s to %s = %f (FAILED %s)\n", #VAR, #NEWVAL, NEWVAL, #COND);\  
         continue;\  
     }  
   
         //CHECK_RESET_CONTINUE(pi < 0, VAR_PI, pi_xiang);  
         //CHECK_RESET_CONTINUE(delg < 0, VAR_DELG, delg_chouaieb);  
         //CHECK_RESET_CONTINUE(delf < 0, VAR_DELF, delf_rackett);  
   
 #if 0  
         CHECK_RESET_CONTINUE(delf < 1, VAR_DELF, delf_rackett);  
         CHECK_RESET_CONTINUE(delg > 1, VAR_DELG, delg_chouaieb);  
   
         CHECK_RESET_CONTINUE(pi > 2. * pi_xiang, VAR_PI, pi_xiang)  
         else CHECK_RESET_CONTINUE(pi < 0.5 * pi_xiang, VAR_PI, pi_xiang)  
   
         CHECK_RESET_CONTINUE(delg < 0.8 * delg_chouaieb, VAR_DELG, delg_chouaieb)  
         else CHECK_RESET_CONTINUE(delg > 1.5 * delg_chouaieb, VAR_DELG, delg_chouaieb);  
   
         //if(gsl_vector_get(s->x,1) < D->rho_c)gsl_vector_set(s->x,1,2 * D->rho_c);  
         //if(gsl_vector_get(s->x,2) > D->rho_c)gsl_vector_set(s->x,2,0.5 * D->rho_c);  
   
 #endif  
     }while(status == GSL_CONTINUE && iter < 40);  
   
     if(status!=GSL_SUCCESS)fprintf(stderr,"%s:%d: warning: (%d) status = %s\n", __FILE__,__LINE__, status, gsl_strerror (status));  
   
     //fprintf(stderr,"SOLUTION: pi = %f\n", gsl_vector_get(s->x, 0));  
   
   
     //fprintf(stderr,"          p  = %f\n", p_sat);  
   
     *p_sat = gsl_vector_get(s->x, 0) * T * D->R * D->rho_c;  
     *rho_f = gsl_vector_get(s->x, 1) * D->rho_c;  
     *rho_g = gsl_vector_get(s->x, 2) * D->rho_c;  
203    
     gsl_multiroot_fsolver_free(s);  
     gsl_vector_free(x);  
   
     return status;  
 }  
 #endif  
   
 #if 0  
       
   
 #ifdef TEST  
     fprintf(stderr,"PHASE CRITERION: T = %f, rho_f = %f, rho_g = %f, p_sat = %f\n", T, rho_f, rho_g, p_sat);  
 #endif  
     double delta_f, delta_g, tau;  
     tau = D->T_c / T;  
   
     delta_f = rho_f / D->rho_c;  
     delta_g = rho_g/ D->rho_c;  
   
 #ifdef TEST  
     assert(!isnan(delta_f));  
     assert(!isnan(delta_g));  
     assert(!isnan(p_sat));        
     assert(!isinf(delta_f));  
     assert(!isinf(delta_g));  
     assert(!isinf(p_sat));        
 #endif  
   
     *eq1 = (p_sat - helmholtz_p(T, rho_f, D));  
     *eq2 = (p_sat - helmholtz_p(T, rho_g, D));  
     *eq3 = helmholtz_g(T, rho_f, D) - helmholtz_g(T,rho_g, D);  
   
 #ifdef TEST  
     fprintf(stderr,"eq1 = %e\t\teq2 = %e\t\teq3 = %e\n", *eq1, *eq2, *eq3);  
 #endif  
 }  
   
   
   
 #if 0  
 static int phase_deriv(const gsl_vector * x, void *params, gsl_matrix * J){  
     const double tau = ((PhaseSolve *)params)->tau;  
     const double pi = gsl_vector_get(x,0);  
     const double delta_f = gsl_vector_get(x,1);  
     const double delta_g = gsl_vector_get(x,2);  
     const HelmholtzData *D = ((PhaseSolve *)params)->D;  
   
     double dE1ddelf = helm_resid_del(tau, delta_f, D) + delta_f * helm_resid_deldel(tau, delta_f, D) - pi ;  
     double dE1ddelg = 0;  
     double dE1dpi = - 1./ delta_f;  
   
     double dE2ddelf = 0;  
     double dE2ddelg = helm_resid_del(tau, delta_g, D) + delta_g * helm_resid_deldel(tau, delta_g, D) - pi ;  
     double dE2dpi = - 1./ delta_g;  
   
     double dE3ddelf = pi - delta_g * (log(delta_f/delta_g) + helm_resid(delta_g, tau, D) - helm_resid(delta_f, tau, D) + 1 - delta_f * helm_resid_del(delta_f, tau, D));  
     double dE3ddelg = -pi - delta_f * (log(delta_f/delta_g) + helm_resid(delta_g, tau, D) - helm_resid(delta_f, tau, D) - 1 + delta_g * helm_resid_del(delta_g, tau, D));  
     double dE3dpi = delta_f - delta_g;  
   
     gsl_matrix_set (J, 0, 0, dE1ddelf);  
     gsl_matrix_set (J, 0, 1, dE1ddelg);  
     gsl_matrix_set (J, 0, 2, dE1dpi);  
     gsl_matrix_set (J, 1, 0, dE2ddelf);  
     gsl_matrix_set (J, 1, 1, dE2ddelg);  
     gsl_matrix_set (J, 1, 2, dE2dpi);  
     gsl_matrix_set (J, 2, 0, dE3ddelf);  
     gsl_matrix_set (J, 2, 1, dE3ddelg);  
     gsl_matrix_set (J, 2, 2, dE3dpi);  
   
     return GSL_SUCCESS;  
 }  
   
 static int phase_residderiv(  
     const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J  
 ){  
     phase_resid(x, params, f);  
     phase_deriv(x, params, J);  
     return GSL_SUCCESS;  
 }  
   
 static int print_state_fdf(size_t iter, gsl_multiroot_fdfsolver * s){  
     fprintf(stderr,"iter = %3u: delf = %.3f delg = %.3f pi = %.3f "  
         "E = % .3e % .3e %.3e\n",  
         iter,  
         gsl_vector_get (s->x, 0),  
         gsl_vector_get (s->x, 1),  
         gsl_vector_get (s->x, 2),  
         gsl_vector_get (s->f, 0),  
         gsl_vector_get (s->f, 1),  
         gsl_vector_get (s->f, 2));  
 }  
   
   
 double phase_solve_fdf(double T, const HelmholtzData *D){  
     gsl_multiroot_fdfsolver *s;  
     int status;  
     const size_t n = 3;  
     double tau = D->T_c / T;  
     size_t i, iter = 0;  
     const gsl_multiroot_fdfsolver_type *fdftype;  
   
     PhaseSolve p = {tau, D};  
     gsl_multiroot_function_fdf f = {  
         &phase_resid,  
         &phase_deriv,  
         &phase_residderiv,  
         n, &p  
     };  
   
     double x_init[3] = {  
         /* first guesses, such as they are... */  
         fprops_psat_T_xiang(T, D) / D->R / T / D->rho_c  
         ,fprops_rhof_T_rackett(T, D) / D->rho_c  
         ,fprops_rhog_T_chouaieb(T, D) / D->rho_c  
     };  
   
     gsl_vector *x = gsl_vector_alloc (n);  
     for(i=0; i<3; ++i)gsl_vector_set(x, i, x_init[i]);  
     fdftype = gsl_multiroot_fdfsolver_hybridsj;  
     s = gsl_multiroot_fdfsolver_alloc (fdftype, n);  
     gsl_multiroot_fdfsolver_set (s, &f, x);  
   
     print_state_fdf(iter, s);  
   
     do{  
         iter++;  
         status = gsl_multiroot_fdfsolver_iterate (s);  
         print_state_fdf(iter, s);  
         if(status)break;  
         if(gsl_vector_get(s->x,2) > D->rho_c)gsl_vector_set(s->x,2,D->rho_c);  
         status = gsl_multiroot_test_residual(s->f, 1e-7);  
     }while(status == GSL_CONTINUE && iter < 1000);  
   
     fprintf(stderr,"status = %s\n", gsl_strerror (status));  
   
     fprintf(stderr,"SOLUTION: pi = %f\n", gsl_vector_get(s->f, 0));  
   
     double p_sat = gsl_vector_get(s->f, 0) * T * D->R * D->rho_c;  
   
     fprintf(stderr,"          p  = %f\n", p_sat);  
     gsl_multiroot_fdfsolver_free (s);  
     gsl_vector_free (x);  
   
     return 0;  
 }  
 #endif  
   
   
 void solve_saturation(double T, const HelmholtzData *D){  
     double rho_f, rho_g, p_sat;  
   
     /* first guesses, such as they are... */  
     p_sat = fprops_psat_T_xiang(T, D);  
     rho_f = fprops_rhof_T_rackett(T, D);  
     rho_g = fprops_rhog_T_chouaieb(T, D);  
   
       
   
 #ifdef TEST  
     fprintf(stderr,"Rackett liquid density: %f\n", rho_f);    
     fprintf(stderr,"Chouaieb vapour density: %f\n", rho_g);  
 #endif  
   
     double delta_f, delta_g, eq1, eq2, eq3, tau;  
           
     int i = 40;  
     while(--i > 0){  
         delta_f = rho_f / D->rho_c;  
         delta_g = rho_g/ D->rho_c;  
   
 #ifdef TEST  
         assert(!isnan(delta_f));  
         assert(!isnan(delta_g));  
         assert(!isnan(p_sat));        
         assert(!isinf(delta_f));  
         assert(!isinf(delta_g));  
         assert(!isinf(p_sat));        
 #endif  
           
         phase_criterion(T, rho_f, rho_g, p_sat, &eq1, &eq2, &eq3, D);  
               
 #ifdef TEST      
         fprintf(stderr,"p_sat = %f\trho_f = %f\trho_g = %f\teq1 = %e\t\teq2 = %e\t\teq3 = %e\n",p_sat, rho_f, rho_g, eq1, eq2, eq3);  
         assert(!isnan(eq1));  
         assert(!isnan(eq2));  
         assert(!isnan(eq3));  
         assert(!isinf(eq1));          
         assert(!isinf(eq2));          
         assert(!isinf(eq3));          
 #endif  
   
         double p1 = D->R * T * D->rho_c * delta_f * delta_g / (delta_f - delta_g) * (helm_resid(delta_f,tau,D) - helm_resid(delta_g,tau,D) + log(delta_f/delta_g));  
         if(p1/p_sat > 1.5){  
             p1 = p_sat * 1.5;  
         }else if(p1/p_sat < 0.7){  
             p1 = p_sat * 0.7;  
         }  
         p_sat = p1;  
         if(p_sat < D->p_t) p_sat = D->p_t;  
         if(p_sat > D->p_c) p_sat = D->p_c;  
   
         double rho_f1 = rho_f - (helmholtz_p(T, rho_f, D) - p_sat) / helmholtz_dpdrho_T(T, rho_f, D);  
         double rho_g1 = rho_g - (helmholtz_p(T, rho_g, D) - p_sat) / helmholtz_dpdrho_T(T, rho_g, D);  
   
         if(rho_g1 > 0){  
             rho_g = rho_g1;  
         }  
         if(rho_f1 > 0){  
             rho_f = rho_f1;  
         }  
         if(fabs(rho_f - rho_g) < 1e-5){  
 #ifdef TEST  
             fprintf(stderr,"%s:%d: rho_f = rho_g\n", __FILE__, __LINE__);  
             exit(1);  
 #else  
             break;  
 #endif  
         }  
   
   
   
     }  
   
     //return eq3;  
 }  
   
   
 #endif  
204    
205    
206    

Legend:
Removed from v.2261  
changed lines
  Added in v.2262

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