/[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 2266 - (show annotations) (download) (as text)
Thu Aug 5 13:53:01 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 6182 byte(s)
The (p,h) function is converging almost everywhere. The home-made Newton iteration wasn't starting well for supercritical states
when the starting guess had rho = rho_c, so changing slightly away from fixed almost all the problems.
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
27 #define FPE_DEBUG
28 #ifdef FPE_DEBUG
29 # include <assert.h>
30 #else
31 # define assert(ARGS...)
32 #endif
33
34 #include <setjmp.h>
35 #include <signal.h>
36 #define _GNU_SOURCE
37 #include <fenv.h>
38
39 #define SQ(X) ((X)*(X))
40
41 #ifdef SOLVE_PH_DEBUG
42 # define MSG(STR,...) fprintf(stderr,"%s:%d: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
43 #else
44 # define MSG(ARGS...)
45 #endif
46
47 #define ERRMSG(STR,...) fprintf(stderr,"%s:%d: ERROR: " STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
48
49
50 int fprops_region_ph(double p, double h, const HelmholtzData *D){
51
52 double Tsat, rhof, rhog;
53
54 if(p >= D->p_c)return FPROPS_NON;
55
56 int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);
57
58 double hf = helmholtz_h(Tsat, rhof, D);
59 if(h <= hf)return FPROPS_NON;
60
61 double hg = helmholtz_h(Tsat, rhog, D);
62 if(h >= hg)return FPROPS_NON;
63
64 return FPROPS_SAT;
65 }
66
67 typedef void SignalHandler(int);
68
69 jmp_buf mark;
70 void fprops_fpe(int sig){
71 feclearexcept(FE_DIVBYZERO|FE_INVALID|FE_OVERFLOW);
72 longjmp(mark, -1);
73 }
74
75 int fprops_solve_ph(double p, double h, double *T, double *rho, int use_guess, const HelmholtzData *D){
76 double Tsat, rhof, rhog, hf, hg;
77 double T1, rho1;
78 int subcrit_pressure = 0;
79
80 //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
81 SignalHandler *old = signal(SIGFPE,&fprops_fpe);
82 int jmpret = setjmp(mark);
83 if(jmpret==0){
84
85 if(p < D->p_c){
86 int res = fprops_sat_p(p, &Tsat, &rhof, &rhog, D);
87 if(res){
88 ERRMSG("Unable to solve saturation state");
89 return res;
90 }
91 hf = helmholtz_h(Tsat, rhof, D);
92 hg = helmholtz_h(Tsat, rhog, D);
93 MSG("hf = %f kJ/kg, hg = %f",hf/1e3,hg/1e3);
94
95 if(h > hf && h < hg){
96 MSG("SATURATION REGION");
97 /* saturation region... easy */
98 double x = (h - hf)/(hg - hf);
99 *rho = 1./(x/rhog + (1.-x)/rhof);
100 *T = Tsat;
101 return 0;
102 }
103
104 subcrit_pressure = 1;
105 if(!use_guess){
106 *T = 1.1 * Tsat;
107 if(h <= hf){
108 *rho = rhof;
109 MSG("LIQUID GUESS: T = %f, rho = %f",*T, *rho);
110 }else{
111 *rho = rhog * 0.5;
112 MSG("GAS GUESS: T = %f, rho = %f",*T, *rho);
113 }
114 }
115 }else{
116 if(!use_guess){
117 *T = D->T_c * 1.03;
118 *rho = D->rho_c * 0.99;
119 }
120 }
121
122 MSG("STARTING NON-SAT ITERATION");
123 //*rho = 976.82687191126922;
124 //*T = 344.80371310850518;
125
126 T1 = *T;
127 rho1 = *rho;
128 assert(!__isnan(T1));
129 assert(!__isnan(rho1));
130
131 /* try our own home-baked newton iteration */
132 int i = 0;
133 double delta_T = 0;
134 double delta_rho = 0;
135 MSG("STARTING ITERATION");
136 while(i++ < 60){
137 double p1 = helmholtz_p_raw(T1,rho1,D);
138 assert(!__isnan(p1));
139 double h1 = helmholtz_h_raw(T1,rho1,D);
140 assert(!__isnan(h1));
141
142 MSG(" T = %f, rho = %f\tp = %f bar, h = %f kJ/kg", T1, rho1, p1/1e5, h1/1e3);
143 MSG(" p error = %f bar",(p1 - p)/1e5);
144 MSG(" h error = %f kJ/kg",(h1 - h)/1e3);
145
146 if(p1 < 0){
147 T1 -= (delta_T *= 0.5);
148 rho1 -= (delta_rho *= 0.5);
149 continue;
150 }
151
152 if(fabs(p1 - p) < 1e-6 && fabs(h1 - h) < 1e-8){
153 MSG("Converged to T = %f, rho = %f, in homebaked Newton solver", T1, rho1);
154 *T = T1;
155 *rho = rho1;
156 return 0;
157 }
158 /* calculate step */
159 double f = log(p1) - log(p);
160 double g = h1 - h;
161 assert(!__isnan(f));
162 assert(!__isnan(g));
163 assert(rho1 != 0);
164 assert(p1 != 0);
165 assert(!__isnan(fprops_non_dZdT_v('p', T1,rho1, D)));
166 double f_T = 1./p1 * fprops_non_dZdT_v('p', T1,rho1, D);
167 double f_rho = -1./p1/SQ(rho1) * fprops_non_dZdv_T('p', T1, rho1, D);
168 double g_T = fprops_non_dZdT_v('h', T1,rho1, D);
169 double g_rho = -1./SQ(rho1) * fprops_non_dZdv_T('h', T1, rho1, D);
170 assert(!__isnan(f_T));
171
172 if(__isnan(f_rho)){
173 ERRMSG(" rho1 = %f, T1 = %f",rho1, T1);
174 }
175 assert(!__isnan(f_rho));
176
177 assert(!__isnan(g_T));
178 assert(!__isnan(g_rho));
179
180 double det = g_rho * f_T - f_rho * g_T;
181 assert(det!=0);
182 assert(!__isnan(det));
183 MSG(" ���f/���T = %e\t\t���f/���rho = %e",f_T, f_rho);
184 MSG(" ���g/���T = %e\t\t���g/���rho = %e",g_T, g_rho);
185
186 delta_T = -1./det * (g_rho * f - f_rho * g);
187 delta_rho = -1./det * (f_T * g - g_T * f);
188 assert(!__isnan(delta_T));
189 assert(!__isnan(delta_rho));
190 MSG(" ��T = %f", delta_T);
191 MSG(" ��rho = %f", delta_rho);
192
193 if(subcrit_pressure){
194 if(h > hg){
195 /* vapour */
196 int i = 0;
197 while(rho1 + delta_rho > rhog && i++ < 20){
198 delta_rho *= 0.5;
199 }
200 }else{
201 /* liquid */
202 while(rho1 + delta_rho < rhof && i++ < 20){
203 delta_rho *= 0.5;
204 }
205 if(T1 + delta_T < D->T_t) delta_T = 0.5 * (T1 + D->T_t);
206 if(T1 < D->T_t) delta_T = +1;
207 }
208 }else{
209 #if 0
210 /* supercritical... stay above critical temperature of density > rho_crit */
211 if(rho1 + delta_rho > D->rho_c && T1 + delta_T < D->T_c){
212 delta_T = 0.5 *(T1 + D->T_c);
213 }
214 #endif
215 }
216 /* don't go too dense */
217 if(rho1 + delta_rho > 2000) delta_rho = 2000;
218
219 /* don't go too hot */
220 if(T1 + delta_T > 5000) delta_T = 5000 - T1;
221
222 /* avoid huge step */
223 while(fabs(delta_T / T1) > 0.6){
224 delta_T *= 0.5;
225 }
226 while(fabs(delta_rho / rho1) > 0.2){
227 delta_rho *= 0.5;
228 }
229
230 T1 = T1 + delta_T;
231 rho1 = rho1 + delta_rho;
232 }
233 }else{
234 /* an FPE occurred */
235 MSG("An FPE occurred");
236 }
237
238 *T = T1;
239 *rho = rho1;
240 ERRMSG("Iteration failed");
241 return 999;
242
243 int res = fprops_nonsolver('p','h',p,h,T,rho,D);
244 ERRMSG("Iteration failed in nonsolver");
245 return res;
246 }
247

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