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

Revision 2225 - (show annotations) (download) (as text)
Mon Jul 26 03:40:18 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 6995 byte(s)
```Added python bindings for new function.
We can see that it's not converging for about 10% of cases.
```
 1 #include "sat2.h" 2 3 #include "helmholtz_impl.h" 4 #include 5 #include 6 #include 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 19 20 @return 0 on success, else error code. 21 @param p output, pressure 22 @param rhof output, liquid density 23 @param rhof output, vapour density 24 */ 25 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; 27 28 /* first guess using acentric factor */ 29 p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.)); 30 31 MSG("Initial estimate: psat(%f) = %f bar\n",T,p/1e5); 32 33 int i,j; 34 char use_guess = 0; 35 36 #define FPROPS_MAX_SUCCSUBS 200 37 38 double rhof = -1, rhog = -1; 39 40 for(i=0;i d->p_c || rhof < d->rho_c){ 92 MSG("FPROPS: invalid converged value of p, beyond p_crit\n"); 93 *p_out = p; 94 *rhof_out = rhof; 95 *rhog_out = rhog; 96 return 3; 97 98 } 99 100 MSG(" recalc p = %f bar\n",p/1e5); 101 *p_out = p; 102 *rhof_out = rhof; 103 *rhog_out = rhog; 104 return 0; 105 } 106 107 delta_p_old = delta_p; 108 109 if(fabs(delta_p) > 0.4 * p){ 110 /* reduce the size of delta_p steps */ 111 for(j=0;j<10;++j){ 112 delta_p *= 0.5; 113 if(fabs(delta_p) < 0.4 * p)break; 114 } 115 MSG("FPROPS: delta_p was too large, has been shortened\n"); 116 } 117 118 p = p + delta_p; 119 } 120 121 MSG("FPROPS: too many iterations in %s, or need to try alternative method\n",__func__); 122 *p_out = p; 123 *rhof_out = rhof; 124 *rhog_out = rhog; 125 return 1; 126 } 127 128 129 /* 130 Iterate to find density as a function of pressure and temperature 131 132 @param use_guess whether (1) or not (0) to use first guess value of *rho. 133 @param rho initial guess (in), and return value (out) 134 @return non-zero on error 135 */ 136 int fprops_rho_pT(double p, double T, char phase, char use_guess, const HelmholtzData *d, double *rho){ 137 138 double tol = 1e-8; 139 double vlog, Dguess, dpdrho; 140 double vlog_prev; 141 int i; 142 143 #if 0 144 if(p < 1e-11 /*Pa*/){ 145 if(T > 0){ 146 return /* estimate based on ideal gas equation */ 147 } 148 } 149 #endif 150 151 double vc = 1./d->rho_c; 152 if(vc <= 0)vc = 1; 153 double vclog = log(vc); 154 155 char use_liquid_calc = 0; 156 157 if(p < d->p_c * (1. + 7.*(T/d->T_c - 1.)) && T > d->T_c){ 158 use_liquid_calc = 0; 159 }else if(p > d->p_c || T > d->T_c){ 160 use_liquid_calc = 1; 161 }else if(phase == 'L'){ 162 use_liquid_calc = 1; 163 }else{ 164 use_liquid_calc = 0; 165 } 166 167 char bad_guess = 0; 168 if(use_guess){ 169 if(*rho > 0)vlog = log(1. / *rho); 170 else bad_guess = 1; 171 } 172 173 if(!use_guess || bad_guess){ 174 if(p > d->p_c || T > d->T_c){ 175 /* supercritical density estimate */ 176 vlog = log(0.5 * vc); 177 }else if(use_liquid_calc){ 178 /* first guess liquid density */ 179 Dguess = fprops_rhof_T_rackett(T, d); 180 181 /* TODO sure first guess is within permitted upper bounds */ 182 183 /* if Rackett estimate gives dp/drho < 1, increase guess */ 184 int k = 0; 185 while(k++ < 8){ 186 double dpdrho = helmholtz_dpdrho_T(T,Dguess,d); 187 if(dpdrho > 0)break; 188 Dguess *= 1.1; 189 }; 190 vlog = log(1./Dguess); 191 192 /* upper bound our first guess by log(0.5*v_crit) */ 193 double vlog_max = log(0.5 * vc); 194 if(vlog > vlog_max)vlog = vlog_max; 195 }else{ 196 /* for gas, first guess is to assume ideal gas */ 197 vlog=log(d->R * T / p); 198 } 199 200 if(!use_liquid_calc)vlog=log(d->R*T/p); 201 } 202 203 204 if (use_liquid_calc){ 205 206 /* ==== liquid phase iteration ==== */ 207 208 for(i=0; i<100; ++i){ 209 210 vlog_prev=vlog; 211 212 *rho=1./exp(vlog); 213 double p2 = helmholtz_p(T,*rho,d); 214 double dpdrho = helmholtz_dpdrho_T(T,*rho,d); 215 216 if(dpdrho < 0){ 217 // reduce spec volume if we have dp/drho < 0 218 vlog=vlog-0.1; 219 }else{ 220 double newt; 221 double dpdlv = -*rho * dpdrho; 222 newt = (p2-p) / dpdlv; // first-order Newton method 223 224 /* loosen the tolerance if we're not getting there */ 225 if(i == 10)tol = tol*10; 226 if(i == 15)tol = tol*10; 227 228 if(fabs(newt/p) < 0.0001*tol){ 229 // iteration has converged 230 *rho = 1./exp(vlog - newt); 231 return 0; 232 }else{ 233 // next guess 234 vlog = vlog_prev - newt; 235 236 // keep step-sizes under control... 237 if(fabs(vlog-vlog_prev) > 0.1 && T < 1.5*d->T_c){ 238 vlog = vlog_prev + (vlog-vlog_prev > 0 ? 0.1 : -0.1); 239 } 240 if(vlog > vclog && T < d->T_c){ 241 vlog=0.5*(vlog_prev+vclog); 242 } 243 if (vlog < -5.0 && T > d->T_c){ 244 use_liquid_calc = 0; 245 MSG("FPROPS: switching to vapour calc\n"); 246 break; /* switch to vapour calc */ 247 } 248 } 249 } 250 251 } 252 253 // iteration has not converged 254 *rho = 1./exp(vlog); 255 if (T > d->T_c){ 256 use_liquid_calc = 0; 257 } 258 259 if(use_liquid_calc){ 260 MSG("FPROPS: liquid iteration did not converge\n"); 261 return 1; 262 } 263 } 264 265 /* ==== vapour phase iteration ==== */ 266 267 double plog=log(p); 268 for(i=0; i<30; ++i){ 269 vlog_prev = vlog; 270 *rho = 1./exp(vlog); 271 double p2 = helmholtz_p(T,*rho,d); 272 double dpdrho = helmholtz_dpdrho_T(T,*rho,d); 273 if(dpdrho < 0. || p2 <= 0.){ 274 // increase spec vol if dpdrho < 0 or guess pressure < 0 275 vlog += 0.1; 276 }else{ 277 double dpdlv = - *rho * dpdrho; 278 double fvdpl = (log(p2) - plog)*p2/dpdlv; // first-order Newton method 279 280 if(fabs(fvdpl) > 25.)fvdpl=0.1e-5; 281 282 // loosen tol if we're not getting there 283 if(i == 10)tol *= 10; 284 if(i == 15)tol *= 10; 285 286 if(fabs(fvdpl) < 0.0001*tol){ 287 // converged, success 288 *rho = 1./exp(vlog); 289 return 0; 290 }else{ 291 // next guess 292 vlog=vlog-fvdpl; 293 } 294 } 295 } 296 297 // not converged 298 *rho=1./exp(vlog); 299 MSG("FPROPS: liquid iteration did not converge\n"); 300 return 1; 301 } 302 303

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