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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2774 - (show annotations) (download) (as text)
Thu Mar 27 11:42:30 2014 UTC (3 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 8915 byte(s)
thcond and mu accessible from Python
1 /*
2 * FPROPS fluid property calculation library
3 * Copyright (C) 2011 - John Pye
4 *
5 * ASCEND is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * ASCEND 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 ASCEND; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor,
18 * Boston, MA 02110-1301 USA
19 */
20
21 /*
22 * This file should contain general functions which are common to all equations of state.
23 */
24
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <string.h>
28 #include "rundata.h"
29 #include "filedata.h"
30 #include "fprops.h"
31 #include "helmholtz.h"
32 #include "ideal.h"
33 #include "sat.h"
34 //#include "redkw.h"
35 #include "pengrob.h"
36 #include "visc.h"
37 #include "thcond.h"
38 //#include "mbwr.h"
39
40 #define FPR_DEBUG
41 #ifdef FPR_DEBUG
42 # include "color.h"
43 # define MSG FPROPS_MSG
44 # define ERRMSG FPROPS_ERRMSG
45 #else
46 # define MSG(ARGS...) ((void)0)
47 # define ERRMSG(ARGS...) ((void)0)
48 #endif
49
50 #include <stdio.h>
51
52 // FIXME convert to use of ideal + departure function here in this file (not at the correlation level)
53 // ... or doesn't that work? get further with PR EOS first?
54
55
56 int fprops_corr_avail(const EosData *E, const char *corrtype){
57 if(corrtype==NULL){
58 /* if none specify, return the 'best' we can manage (not yet very civilised) */
59 switch(E->type){
60 case FPROPS_HELMHOLTZ:
61 case FPROPS_IDEAL:
62 return E->type;
63 case FPROPS_CUBIC:
64 return FPROPS_PENGROB;
65 default:
66 return 0;
67 }
68 }else if(strcmp(corrtype,"helmholtz")==0){
69 switch(E->type){
70 case FPROPS_HELMHOLTZ:
71 return FPROPS_HELMHOLTZ;
72 default:
73 return 0;
74 }
75 }else if(strcmp(corrtype,"pengrob")==0){
76 switch(E->type){
77 case FPROPS_HELMHOLTZ:
78 case FPROPS_CUBIC:
79 return FPROPS_PENGROB;
80 default:
81 return 0;
82 }
83 }else if(strcmp(corrtype,"ideal")==0){
84 switch(E->type){
85 case FPROPS_IDEAL:
86 case FPROPS_HELMHOLTZ:
87 case FPROPS_CUBIC:
88 return FPROPS_IDEAL;
89 default:
90 return 0;
91 }
92 }
93 return 0;
94 }
95
96
97 PureFluid *fprops_prepare(const EosData *E,const char *corrtype){
98 PureFluid *P = NULL;
99 FpropsError err = FPROPS_NO_ERROR;
100 MSG("Working with EosData name '%s', source '%s", E->name, E->source);
101 MSG("Chosen correlation: %d (requested %s)", fprops_corr_avail(E,corrtype),corrtype);
102 switch(fprops_corr_avail(E,corrtype)){
103 case FPROPS_HELMHOLTZ:
104 P = helmholtz_prepare(E,NULL);
105 break;
106 case FPROPS_PENGROB:
107 P = pengrob_prepare(E,NULL);
108 break;
109 case FPROPS_IDEAL:
110 P = ideal_prepare(E,NULL);
111 break;
112 default:
113 ERRMSG("Invalid EOS data, unimplemented correlation type requested");
114 return NULL;
115 }
116 /* next: add preparation of viscosity, thermal conductivity, surface tension, ... */
117
118 MSG("Preparing viscosity data...");
119 P->visc = visc_prepare(E,P,&err);
120 if(err){
121 ERRMSG("Invalid viscosity data for '%s",P->name);
122 /* visc_prepare should return NULL if there was an error, so result is
123 same as when there is no viscosity data at all */
124 }
125
126 MSG("Preparing thermal conductivity data...");
127 err = FPROPS_NO_ERROR;
128 thcond_prepare(P,E->thcond,&err);
129 if(err){
130 ERRMSG("Invalid viscosity data for '%s",P->name);
131 /* visc_prepare should return NULL if there was an error, so result is
132 same as when there is no viscosity data at all */
133 }
134
135
136
137 return P;
138 }
139
140 FluidState fprops_set_Trho(double T, double rho, const PureFluid *fluid, FpropsError *err){
141 FluidState state = {T,rho,fluid};
142 return state;
143 }
144
145 // FIXME XXX not all properties are mass-weighted in the saturation region...
146
147 /*
148 TODO need some kind of support for freezing line, especially for water
149 since for temperatures 0.degC up to 0.01 degC liquid phase is possible
150 for ambient pressure, but this is BELOW triple line. Hmmm.
151
152 Also sublimation curve needs to be added.
153 */
154
155 #define EVALFN(VAR) \
156 double fprops_##VAR(FluidState state, FpropsError *err){\
157 double p, rho_f, rho_g;\
158 if(state.T >= state.fluid->data->T_t && state.T < state.fluid->data->T_c){\
159 fprops_sat_T(state.T, &p, &rho_f, &rho_g, state.fluid, err);\
160 if(*err){\
161 MSG("Got error %d from saturation calc in %s\n",*err,__func__);\
162 /*return state.fluid->VAR##_fn(state.fluid->data->T_c,state.fluid->data->rho_c,state.fluid->data,err);*/\
163 return 0;\
164 }\
165 if(rho_g < state.rho && state.rho < rho_f){\
166 double x = rho_g*(rho_f/state.rho - 1)/(rho_f - rho_g);\
167 double Qf = state.fluid->VAR##_fn(state.T,rho_f,state.fluid->data,err);\
168 double Qg = state.fluid->VAR##_fn(state.T,rho_g,state.fluid->data,err);\
169 return x*Qg + (1-x)*Qf;\
170 }\
171 }\
172 return state.fluid->VAR##_fn(state.T,state.rho,state.fluid->data,err);\
173 }
174
175 #define EVALFN_SATUNDEFINED(VAR) \
176 double fprops_##VAR(FluidState state, FpropsError *err){\
177 double p, rho_f, rho_g;\
178 if(state.T >= state.fluid->data->T_t && state.T < state.fluid->data->T_c){\
179 fprops_sat_T(state.T, &p, &rho_f, &rho_g, state.fluid, err);\
180 if(*err){\
181 MSG("Got error %d from sat calc in %s\n",*err,__func__);\
182 return state.fluid->data->rho_c;\
183 }\
184 if(rho_g < state.rho && state.rho < rho_f){\
185 *err = FPROPS_VALUE_UNDEFINED;\
186 }\
187 }\
188 return state.fluid->VAR##_fn(state.T,state.rho,state.fluid->data,err);\
189 }
190
191 EVALFN(p); EVALFN(u); EVALFN(h); EVALFN(s); EVALFN(a); EVALFN(g);
192 EVALFN_SATUNDEFINED(cp); EVALFN_SATUNDEFINED(cv);
193 EVALFN_SATUNDEFINED(w);
194 EVALFN(dpdrho_T);
195
196 EVALFN(alphap); EVALFN(betap);
197 // EVALFN(dpdT_rho);
198 //EVALFN(dpdrho_T); EVALFN(d2pdrho2_T); EVALFN(dhdT_rho); EVALFN(dhdrho_T);
199 //EVALFN(dudT_rho); EVALFN(dudrho)T);
200
201 #if 0
202 double fprops_alphap(FluidState state, FpropsError *err){
203 *err = FPROPS_NOT_IMPLEMENTED;
204 return 0;
205 }
206 double fprops_betap(FluidState state, FpropsError *err){
207 *err = FPROPS_NOT_IMPLEMENTED;
208 return 0;
209 }
210 #endif
211
212 /// TODO reimplement with function pointer?
213 double fprops_mu(FluidState state, FpropsError *err){
214 if(NULL!=state.fluid->visc){
215 switch(state.fluid->visc->type){
216 case FPROPS_VISC_1:
217 return visc1_mu(state,err);
218 default:
219 break;
220 }
221 }
222 *err = FPROPS_NOT_IMPLEMENTED;
223 return NAN;
224 }
225
226 /// TODO reimplement with function pointer?
227 double fprops_lam(FluidState state, FpropsError *err){
228 if(NULL!=state.fluid->thcond){
229 switch(state.fluid->thcond->type){
230 case FPROPS_THCOND_1:
231 return thcond1_lam(state,err);
232 default:
233 break;
234 }
235 }
236 *err = FPROPS_NOT_IMPLEMENTED;
237 return NAN;
238 }
239
240
241 double fprops_cp0(FluidState state, FpropsError *err){
242 return ideal_cp(state.T,0,state.fluid->data,err);
243 }
244
245 double fprops_x(FluidState state, FpropsError *err){
246 double p, rho_f, rho_g;
247 if(state.T >= state.fluid->data->T_t && state.T < state.fluid->data->T_c){
248 fprops_sat_T(state.T, &p, &rho_f, &rho_g, state.fluid, err);
249 if(*err)return 0;
250 if(state.rho > rho_f)return 0;
251 if(state.rho < rho_g)return 1;
252 return rho_g*(rho_f/state.rho - 1)/(rho_f - rho_g);
253 }
254 fprintf(stderr,"Temperature is <T_t or >T_c\n");
255 *err = FPROPS_VALUE_UNDEFINED;
256 return 0;
257 }
258
259 double fprops_dpdT_rho(FluidState state, FpropsError *err){
260 *err = FPROPS_NOT_IMPLEMENTED;
261 return 0;
262 }
263
264 double fprops_d2pdrho2_T(FluidState state, FpropsError *err){
265 *err = FPROPS_NOT_IMPLEMENTED;
266 return 0;
267 }
268
269 double fprops_dhdT_rho(FluidState state, FpropsError *err){
270 *err = FPROPS_NOT_IMPLEMENTED;
271 return 0;
272 }
273
274 double fprops_dhdrho_T(FluidState state, FpropsError *err){
275 *err = FPROPS_NOT_IMPLEMENTED;
276 return 0;
277 }
278
279 double fprops_dudT_rho(FluidState state, FpropsError *err){
280 *err = FPROPS_NOT_IMPLEMENTED;
281 return 0;
282 }
283 double fprops_dudrho_T(FluidState state, FpropsError *err){
284 *err = FPROPS_NOT_IMPLEMENTED;
285 return 0;
286 }
287
288 char *fprops_error(FpropsError err){
289 switch (err) {
290 case FPROPS_NO_ERROR:return NULL;
291 case FPROPS_NUMERIC_ERROR:return "FPROPS encountered a numerical error.";
292 case FPROPS_SAT_CVGC_ERROR:return "FPROPS unable to converge solution in saturation region.";
293 case FPROPS_RANGE_ERROR:return "FPROPS had a range error on one of its inputs.";
294 case FPROPS_DATA_ERROR:return "FPROPS encountered a data error.";
295 case FPROPS_NOT_IMPLEMENTED:return "FPROPS feature not yet implemented.";
296 case FPROPS_INVALID_REQUEST:return "FPROPS encountered an invalid request.";
297 case FPROPS_VALUE_UNDEFINED:return "FPROPS reports the request value is locally undefined.";
298 default:return "Unrecognised error";
299 }
300 }
301
302 char *fprops_corr_type(EosType type){
303 switch(type){
304 case FPROPS_IDEAL:
305 return "ideal";
306 case FPROPS_CUBIC:
307 return "cubic";
308 case FPROPS_PENGROB:
309 return "pengrob";
310 case FPROPS_REDKW:
311 return "redkw";
312 case FPROPS_SOAVE:
313 return "soave";
314 case FPROPS_HELMHOLTZ:
315 return "helmholtz";
316 case FPROPS_MBWR:
317 return "mbwr";
318 }
319 return NULL;
320 }
321

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