/[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 2225 - (show annotations) (download) (as text)
Mon Jul 26 03:40:18 2010 UTC (9 years, 11 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 <assert.h>
5 #include <math.h>
6 #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
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<FPROPS_MAX_SUCCSUBS;++i){
41 if(fprops_rho_pT(p,T,FPROPS_PHASE_LIQUID,use_guess, d, &rhof)){
42 MSG(" increasing p, error with rho_f\n");
43 p *= 1.005;
44 continue;
45 }
46
47 if(fprops_rho_pT(p,T,FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
48 MSG(" decreasing p, error with rho_g\n");
49 p *= 0.95;
50 continue;
51 }
52
53 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 */
56
57 if(fabs(rhof - rhog) < 1e-8){
58 MSG("FPROPS: densities converged to same value\n");
59 *p_out = p;
60 *rhof_out = rhof;
61 *rhog_out = rhog;
62 return 1;
63 }
64
65 p_new = (helmholtz_a(T, rhof, d) - helmholtz_a(T, rhog, d)) / (1./rhog - 1./rhof);
66 delta_p = p_new - p;
67
68 //MSG(" delta_p = %f bar\n",delta_p/1e5);
69
70 /* convergence test */
71 if(abs(delta_p/p) < 1e-6){
72 MSG("CONVERGED...\n");
73 /* note possible need for delp non-change test for certain fluids */
74
75 p = p_new;
76
77 /* find vapour density, using guess */
78 if(fprops_rho_pT(p, T, FPROPS_PHASE_VAPOUR, use_guess, d, &rhog)){
79 MSG("FPROPS: fails to estimate vapour density\n");
80 *rhof_out = rhof;
81 *rhog_out = rhog;
82 *p_out = p;
83 return 2;
84 }
85 /* find liquid phase, using guess */
86 fprops_rho_pT(p, T, FPROPS_PHASE_LIQUID, use_guess, d, &rhof);
87
88 /* re-evaluate exact p for the value of rhof calculated */
89 p = helmholtz_p(T,rhof,d);
90
91 if(p > 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