/[ascend]/trunk/models/johnpye/fprops/sat2.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2235 - (show annotations) (download) (as text)
Thu Jul 29 08:45:54 2010 UTC (9 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 7546 byte(s)
we'll satisfy ourselves with a finite difference approximation of helm_resid_deldeldel for the moment, and
move on with the saturation curve problem.
1 #include "sat2.h"
2
3 #include "helmholtz_impl.h"
4 #include <assert.h>
5 #include <math.h>
6 #include <stdio.h>
7
8 //#define PHASE_DEBUG
9 #define PHASE_ERRORS
10
11 #ifndef PHASE_DEBUG
12 # define MSG(...)
13 #else
14 # define MSG(ARGS...) fprintf(stderr,"FPROPS: " ARGS)
15 #endif
16
17 #ifndef PHASE_DEBUG
18 # define ERR(...)
19 #else
20 # define ERR(ARGS...) fprintf(stderr,"FPROPS: " ARGS)
21 #endif
22
23
24 /**
25 solve Maxwell phase criterion using successive substitution
26
27 @return 0 on success, else error code.
28 @param p output, pressure
29 @param rhof output, liquid density
30 @param rhof output, vapour density
31 */
32 int fprops_sat_T(double T, double *p_out, double *rhof_out, double *rhog_out, const HelmholtzData *d){
33 double p, p_new, delta_p, delta_p_old;
34
35 /*
36 Estimate of saturation pressure assuming algebraic form
37 log10(p)=A + B/T
38 See Reid, Prausnitz and Poling, 4th Ed., section 2.3.
39 */
40 p = d->p_c * pow(10, -7./3 * (1.+d->omega) * (d->T_c / T - 1.));
41
42 MSG("Initial estimate: psat(%f) = %f bar\n",T,p/1e5);
43
44 int i,j;
45 char use_guess = 0;
46
47 #define FPROPS_MAX_SUCCSUBS 12
48
49 double rhof = -1, rhog = -1;
50
51 /* this approach follows approximately the approach of Soave, page 204:
52 http://dx.doi.org/10.1016/0378-3812(86)90013-0
53 */
54 for(i=0;i<FPROPS_MAX_SUCCSUBS;++i){
55 if(fprops_rho_pT(p,T,FPROPS_PHASE_LIQUID,use_guess, d, &rhof)){
56 MSG(" increasing p, error with rho_f\n");
57 p *= 1.005;
58 continue;
59 }
60
61 if(fprops_rho_pT(p,T,FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
62 MSG(" decreasing p, error with rho_g\n");
63 p *= 0.95;
64 continue;
65 }
66
67 MSG("Iter %d: p = %f bar, rhof = %f, rhog = %f\n",i, p/1e5, rhof, rhog);
68
69 use_guess = 1; /* after first run, start re-using current guess */
70
71 if(fabs(rhof - rhog) < 1e-8){
72 ERR("densities converged to same value\n");
73 *p_out = p;
74 *rhof_out = rhof;
75 *rhog_out = rhog;
76 return 1;
77 }
78
79 double delta_a = helmholtz_a(T, rhof, d) - helmholtz_a(T, rhog, d);
80 MSG(" delta_a = %f\n",delta_a);
81
82 p_new = delta_a / (1./rhog - 1./rhof);
83 delta_p = p_new - p;
84
85 MSG(" delta_p = %f bar\n",delta_p/1e5);
86
87 /* convergence test */
88 if(fabs(delta_p/p) < 1e-6){
89 MSG("CONVERGED...\n");
90 /* note possible need for delp non-change test for certain fluids */
91
92 p = p_new;
93
94 /* find vapour density, using guess */
95 if(fprops_rho_pT(p, T, FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
96 ERR("failed to estimate vapour density\n");
97 *rhof_out = rhof;
98 *rhog_out = rhog;
99 *p_out = p;
100 return 2;
101 }
102 /* find liquid phase, using guess */
103 fprops_rho_pT(p, T, FPROPS_PHASE_LIQUID, use_guess, d, &rhof);
104
105 /* re-evaluate exact p for the value of rhof calculated */
106 p = helmholtz_p(T,rhof,d);
107
108 if(p > d->p_c || rhof < d->rho_c){
109 ERR("invalid converged value of p, beyond p_crit\n");
110 *p_out = p;
111 *rhof_out = rhof;
112 *rhog_out = rhog;
113 return 3;
114
115 }
116
117 MSG(" recalc p = %f bar\n",p/1e5);
118 *p_out = p;
119 *rhof_out = rhof;
120 *rhog_out = rhog;
121 return 0;
122 }
123
124 delta_p_old = delta_p;
125
126 if(fabs(delta_p) > 0.4 * p){
127 /* reduce the size of delta_p steps */
128 for(j=0;j<10;++j){
129 delta_p *= 0.5;
130 if(fabs(delta_p) < 0.4 * p)break;
131 }
132 MSG("FPROPS: delta_p was too large, has been shortened\n");
133 }
134
135 p = p + delta_p;
136 }
137
138 if(i==FPROPS_MAX_SUCCSUBS){
139 ERR("too many iterations (%d) in %s(T = %f), or need to try alternative method\n",i,__func__,T);
140 *p_out = p;
141 *rhof_out = rhof;
142 *rhog_out = rhog;
143 return 4;
144 }
145
146 ERR("unknown error with T = %f in %s\n",T,__func__);
147 *p_out = p;
148 *rhof_out = rhof;
149 *rhog_out = rhog;
150 return 9;
151 }
152
153
154 /*
155 Iterate to find density as a function of pressure and temperature
156
157 @param use_guess whether (1) or not (0) to use first guess value of *rho.
158 @param rho initial guess (in), and return value (out)
159 @return non-zero on error
160 */
161 int fprops_rho_pT(double p, double T, char phase, char use_guess, const HelmholtzData *d, double *rho){
162
163 double tol = 1e-8;
164 double vlog, Dguess, dpdrho;
165 double vlog_prev;
166 int i;
167
168 #if 0
169 if(p < 1e-11 /*Pa*/){
170 if(T > 0){
171 return /* estimate based on ideal gas equation */
172 }
173 }
174 #endif
175
176 double vc = 1./d->rho_c;
177 if(vc <= 0)vc = 1;
178 double vclog = log(vc);
179
180 char use_liquid_calc = 0;
181
182 if(p < d->p_c * (1. + 7.*(T/d->T_c - 1.)) && T > d->T_c){
183 use_liquid_calc = 0;
184 }else if(p > d->p_c || T > d->T_c){
185 use_liquid_calc = 1;
186 }else if(phase == 'L'){
187 use_liquid_calc = 1;
188 }else{
189 use_liquid_calc = 0;
190 }
191
192 char bad_guess = 0;
193 if(use_guess){
194 if(*rho > 0)vlog = log(1. / *rho);
195 else bad_guess = 1;
196 }
197
198 if(!use_guess || bad_guess){
199 if(p > d->p_c || T > d->T_c){
200 /* supercritical density estimate */
201 vlog = log(0.5 * vc);
202 }else if(use_liquid_calc){
203 /* first guess liquid density */
204 Dguess = fprops_rhof_T_rackett(T, d);
205
206 /* TODO sure first guess is within permitted upper bounds */
207
208 /* if Rackett estimate gives dp/drho < 1, increase guess */
209 int k = 0;
210 while(k++ < 8){
211 double dpdrho = helmholtz_dpdrho_T(T,Dguess,d);
212 if(dpdrho > 0)break;
213 Dguess *= 1.1;
214 };
215 vlog = log(1./Dguess);
216
217 /* upper bound our first guess by log(0.5*v_crit) */
218 double vlog_max = log(0.5 * vc);
219 if(vlog > vlog_max)vlog = vlog_max;
220 }else{
221 /* for gas, first guess is to assume ideal gas */
222 vlog=log(d->R * T / p);
223 }
224
225 if(!use_liquid_calc)vlog=log(d->R*T/p);
226 }
227
228
229 if (use_liquid_calc){
230
231 /* ==== liquid phase iteration ==== */
232
233 for(i=0; i<100; ++i){
234
235 vlog_prev=vlog;
236
237 *rho=1./exp(vlog);
238 double p2 = helmholtz_p(T,*rho,d);
239 double dpdrho = helmholtz_dpdrho_T(T,*rho,d);
240
241 if(dpdrho < 0){
242 // reduce spec volume if we have dp/drho < 0
243 vlog=vlog-0.1;
244 }else{
245 double newt;
246 double dpdlv = -*rho * dpdrho;
247 newt = (p2-p) / dpdlv; // first-order Newton method
248
249 /* loosen the tolerance if we're not getting there */
250 if(i == 10)tol = tol*10;
251 if(i == 15)tol = tol*10;
252
253 if(fabs(newt/p) < 0.0001*tol){
254 // iteration has converged
255 *rho = 1./exp(vlog - newt);
256 return 0;
257 }else{
258 // next guess
259 vlog = vlog_prev - newt;
260
261 // keep step-sizes under control...
262 if(fabs(vlog-vlog_prev) > 0.1 && T < 1.5*d->T_c){
263 vlog = vlog_prev + (vlog-vlog_prev > 0 ? 0.1 : -0.1);
264 }
265 if(vlog > vclog && T < d->T_c){
266 vlog=0.5*(vlog_prev+vclog);
267 }
268 if (vlog < -5.0 && T > d->T_c){
269 use_liquid_calc = 0;
270 MSG("FPROPS: switching to vapour calc\n");
271 break; /* switch to vapour calc */
272 }
273 }
274 }
275
276 }
277
278 // iteration has not converged
279 *rho = 1./exp(vlog);
280 if (T > d->T_c){
281 use_liquid_calc = 0;
282 }
283
284 if(use_liquid_calc){
285 MSG("FPROPS: liquid iteration did not converge\n");
286 return 1;
287 }
288 }
289
290 /* ==== vapour phase iteration ==== */
291
292 double plog=log(p);
293 for(i=0; i<30; ++i){
294 vlog_prev = vlog;
295 *rho = 1./exp(vlog);
296 double p2 = helmholtz_p(T,*rho,d);
297 double dpdrho = helmholtz_dpdrho_T(T,*rho,d);
298 if(dpdrho < 0. || p2 <= 0.){
299 // increase spec vol if dpdrho < 0 or guess pressure < 0
300 vlog += 0.1;
301 }else{
302 double dpdlv = - *rho * dpdrho;
303 double fvdpl = (log(p2) - plog)*p2/dpdlv; // first-order Newton method
304
305 if(fabs(fvdpl) > 25.)fvdpl=0.1e-5;
306
307 // loosen tol if we're not getting there
308 if(i == 10)tol *= 10;
309 if(i == 15)tol *= 10;
310
311 if(fabs(fvdpl) < 0.0001*tol){
312 // converged, success
313 *rho = 1./exp(vlog);
314 return 0;
315 }else{
316 // next guess
317 vlog=vlog-fvdpl;
318 }
319 }
320 }
321
322 // not converged
323 *rho=1./exp(vlog);
324 MSG("FPROPS: liquid iteration did not converge\n");
325 return 1;
326 }
327
328

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