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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2308 - (show annotations) (download) (as text)
Thu Sep 2 03:56:14 2010 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 8062 byte(s)
Turn off FPE in solve_ph.c (otherwise causes problems on Ubuntu 10.04 with rankine_co2 model).
1 /*
2 ASCEND modelling environment
3 Copyright (C) 2004-2010 John Pye
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License
7 as published by the Free Software Foundation; either version 2
8 of the License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18 */
19
20 #include "solve_ph.h"
21 #include "sat.h"
22 #include "derivs.h"
23
24 #include <stdio.h>
25 #include <math.h>
26 #include <stdlib.h>
27
28 //#define FPE_DEBUG
29 #define ASSERT_DEBUG
30
31 #ifdef ASSERT_DEBUG
32 # include <assert.h>
33 #else
34 # define assert(ARGS...)
35 #endif
36
37 #include <setjmp.h>
38 #include <signal.h>
39
40 #ifdef FPE_DEBUG
41 #define _GNU_SOURCE
42 #include <fenv.h>
43 int feenableexcept (int excepts);
44 int fedisableexcept (int excepts);
45 int fegetexcept (void);
46 #endif
47
48 #define SQ(X) ((X)*(X))
49
50 //#define SOLVE_PH_DEBUG
51 #ifdef SOLVE_PH_DEBUG
52 # define MSG(STR,...) fprintf(stderr,"%s:%d: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
53 #else
54 # define MSG(ARGS...)
55 #endif
56
57 #define ERRMSG(STR,...) fprintf(stderr,"%s:%d: ERROR: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
58
59 int fprops_region_ph(double p, double h, const HelmholtzData *D){
60
61 double Tsat, rhof, rhog;
62 double p_c = fprops_pc(D);
63
64 if(p >= p_c)return FPROPS_NON;
65
66 int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);
67 if(res)return FPROPS_ERR;
68
69 double hf = helmholtz_h(Tsat, rhof, D);
70 if(h <= hf)return FPROPS_NON;
71
72 double hg = helmholtz_h(Tsat, rhog, D);
73 if(h >= hg)return FPROPS_NON;
74
75 return FPROPS_SAT;
76 }
77
78 typedef void SignalHandler(int);
79
80 #ifdef FPE_DEBUG
81 jmp_buf mark;
82 void fprops_fpe(int sig){
83 MSG("Catching signal %d!",sig);
84 feclearexcept(FE_DIVBYZERO|FE_INVALID|FE_OVERFLOW);
85 longjmp(mark, -1);
86 }
87 #endif
88
89 int fprops_solve_ph(double p, double h, double *T, double *rho, int use_guess, const HelmholtzData *D){
90 double Tsat, rhof, rhog, hf, hg;
91 double T1, rho1;
92 int subcrit_pressure = 0;
93 int liquid_iteration = 0;
94 double rhof_t;
95 double p_c = fprops_pc(D);
96
97 #ifdef FPE_DEBUG
98 feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
99 SignalHandler *old = signal(SIGFPE,&fprops_fpe);
100 (void)old;/* not used for anything at this stage */
101 int jmpret = setjmp(mark);
102 if(jmpret==0){
103 #endif
104
105 if(p < p_c){
106 MSG("Calculate saturation Tsat(p < p_c)");
107 int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);
108 if(res){
109 ERRMSG("Unable to solve saturation state");
110 return res;
111 }
112 hf = helmholtz_h(Tsat, rhof, D);
113 hg = helmholtz_h(Tsat, rhog, D);
114 MSG("hf = %f kJ/kg, hg = %f",hf/1e3,hg/1e3);
115
116 if(h > hf && h < hg){
117 MSG("SATURATION REGION");
118 /* saturation region... easy */
119 double x = (h - hf)/(hg - hf);
120 *rho = 1./(x/rhog + (1.-x)/rhof);
121 *T = Tsat;
122 return 0;
123 }
124
125 subcrit_pressure = 1;
126 if(h < hf){
127 liquid_iteration = 1;
128 if(!use_guess){
129 #if 1
130 double Tsat1, psat1, rhof1, rhog1;
131 MSG("SOLVING TSAT(HF)");
132 res = fprops_sat_hf(h, &Tsat1, &psat1, &rhof1, &rhog1, D);
133 if(res){
134 ERRMSG("Unable to solve Tsat(hf), returning T = %f, rho = %f.",Tsat,rhof);
135 *T = Tsat;
136 *rho = rhof;
137 return res;
138 }
139 *T = Tsat1;
140 *rho = rhof1;
141 #else
142 *T = Tsat;
143 *rho = rhof;
144 #endif
145 MSG("LIQUID GUESS: T = %f, rho = %f",*T, *rho);
146 }
147 }else if(!use_guess){
148 *T = 1.1 * Tsat;
149 *rho = rhog * 0.5;
150 MSG("GAS GUESS: T = %f, rho = %f",*T, *rho);
151 }
152 }else{
153 /* FIXME still some problems here at very high pressures */
154 if(!use_guess){
155 double hc = helmholtz_h_raw(D->T_c, D->rho_c, D);
156 assert(!__isnan(hc));
157 MSG("hc = %f",hc);
158 if(h < 0.9 * hc){
159 MSG("h < hc... using saturation Tsat(hf) for starting guess");
160 double Tsat1, psat1, rhof1, rhog1;
161 int res = fprops_sat_hf(h, &Tsat1, &psat1, &rhof1, &rhog1, D);
162 if(res){
163 MSG("Unable to solve Tsat(hf)");
164 /* accuracy of the estimate of Tsat(hf) doesn't matter
165 very much so we can ignore this error. */
166 }
167 *T = Tsat1;
168 *rho = rhof1;
169 }else{
170 *T = D->T_c * 1.01;
171 *rho = D->rho_c * 1.05;
172 }
173 MSG("SUPERCRITICAL GUESS: T = %f, rho = %f", *T, *rho);
174 }
175 }
176
177 MSG("STARTING NON-SAT ITERATION");
178 //*rho = 976.82687191126922;
179 //*T = 344.80371310850518;
180
181 T1 = *T;
182 rho1 = *rho;
183 assert(!__isnan(T1));
184 assert(!__isnan(rho1));
185
186 if(liquid_iteration){
187 double pt,rhogt;
188 int res = fprops_sat_T(D->T_t,&pt,&rhof_t, &rhogt, D);
189 if(res){
190 ERRMSG("Unable to solve triple point liquid density.");
191 return 1;
192 }
193 }
194
195 /* try our own home-baked newton iteration */
196 int i = 0;
197 double delta_T = 0;
198 double delta_rho = 0;
199 MSG("STARTING ITERATION");
200 while(i++ < 200){
201 double p1 = helmholtz_p_raw(T1,rho1,D);
202 assert(!__isnan(p1));
203 double h1 = helmholtz_h_raw(T1,rho1,D);
204 assert(!__isnan(h1));
205
206 MSG(" T = %f, rho = %f\tp = %f bar, h = %f kJ/kg", T1, rho1, p1/1e5, h1/1e3);
207 MSG(" p error = %f bar",(p1 - p)/1e5);
208 MSG(" h error = %f kJ/kg",(h1 - h)/1e3);
209
210 if(p1 < 0){
211 T1 -= (delta_T *= 0.5);
212 rho1 -= (delta_rho *= 0.5);
213 continue;
214 }
215
216 if(fabs(p1 - p) < 1e-4 && fabs(h1 - h) < 1e-8){
217 MSG("Converged to T = %f, rho = %f, in homebaked Newton solver", T1, rho1);
218 *T = T1;
219 *rho = rho1;
220 return 0;
221 }
222 /* calculate step, we're solving log(p1) in this code... */
223 double f = log(p1) - log(p);
224 double g = h1 - h;
225 assert(!__isnan(f));
226 assert(!__isnan(g));
227 assert(rho1 != 0);
228 assert(p1 != 0);
229 assert(!__isnan(fprops_non_dZdT_v('p', T1,rho1, D)));
230 double f_T = 1./p1 * fprops_non_dZdT_v('p', T1,rho1, D);
231 double f_rho = -1./p1/SQ(rho1) * fprops_non_dZdv_T('p', T1, rho1, D);
232 double g_T = fprops_non_dZdT_v('h', T1,rho1, D);
233 double g_rho = -1./SQ(rho1) * fprops_non_dZdv_T('h', T1, rho1, D);
234 assert(!__isnan(f_T));
235
236 if(__isnan(f_rho)){
237 ERRMSG(" rho1 = %f, T1 = %f",rho1, T1);
238 }
239 assert(!__isnan(f_rho));
240
241 assert(!__isnan(g_T));
242 assert(!__isnan(g_rho));
243
244 double det = g_rho * f_T - f_rho * g_T;
245 assert(det!=0);
246 assert(!__isnan(det));
247 MSG(" ���f/���T = %e\t\t���f/���rho = %e",f_T, f_rho);
248 MSG(" ���g/���T = %e\t\t���g/���rho = %e",g_T, g_rho);
249
250 delta_T = -1./det * (g_rho * f - f_rho * g);
251 delta_rho = -1./det * (f_T * g - g_T * f);
252 assert(!__isnan(delta_T));
253 assert(!__isnan(delta_rho));
254 MSG(" ��T = %f", delta_T);
255 MSG(" ��rho = %f", delta_rho);
256
257 if(subcrit_pressure){
258 if(h > hg){
259 /* vapour */
260 int i = 0;
261 while(rho1 + delta_rho > rhog && i++ < 20){
262 delta_rho *= 0.5;
263 }
264 }else{
265 /* liquid */
266 while(rho1 + delta_rho < rhof && i++ < 20){
267 delta_rho *= 0.5;
268 }
269 if(T1 + delta_T < D->T_t) delta_T = 0.5 * (T1 + D->T_t);
270 }
271 }else{
272 #if 0
273 /* supercritical... stay above critical temperature of density > rho_crit */
274 if(rho1 + delta_rho > D->rho_c && T1 + delta_T < D->T_c){
275 delta_T = 0.5 *(T1 + D->T_c);
276 }
277 #endif
278 }
279 /* don't go too dense */
280 if(liquid_iteration){
281 if(rho1 + delta_rho > rhof_t) delta_rho = rhof_t;
282 }
283
284 /* don't go too hot */
285 if(T1 + delta_T > 5000) delta_T = 5000 - T1;
286
287 /* avoid huge step */
288 while(fabs(delta_T / T1) > 0.6){
289 delta_T *= 0.5;
290 }
291 while(fabs(delta_rho / rho1) > 0.2){
292 delta_rho *= 0.5;
293 }
294
295 T1 = T1 + delta_T;
296 if(T1 < D->T_t)T1 = D->T_t;
297 rho1 = rho1 + delta_rho;
298 }
299 #ifdef FPE_DEBUG
300 }else{
301 /* an FPE occurred */
302 MSG("An FPE occurred");
303 abort();
304 }
305 #endif
306
307 *T = T1;
308 *rho = rho1;
309 ERRMSG("Iteration failed for p = %.12e, h = %.12e",p,h);
310 return 999;
311
312 #if 0
313 int res = fprops_nonsolver('p','h',p,h,T,rho,D);
314 ERRMSG("Iteration failed in nonsolver");
315 return res;
316 #endif
317 }
318

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