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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2806 - (show annotations) (download) (as text)
Wed Jul 23 01:24:18 2014 UTC (4 years, 2 months ago) by jpye
File MIME type: text/x-csrc
File size: 18868 byte(s)
reduce debug output
1 /* ASCEND modelling environment
2 Copyright (C) 2008 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 *//** @file
17 Wrapper for FPROPS to allow access from ASCEND.
18 */
19
20 //#define ASC_FPROPS_DEBUG
21
22 /* include the external function API from libascend... */
23 #include <ascend/compiler/extfunc.h>
24
25 /* include error reporting API as well, so we can send messages to user */
26 #include <ascend/utilities/error.h>
27
28 /* for accessing the DATA instance */
29 #include <ascend/compiler/child.h>
30 #include <ascend/general/list.h>
31 #include <ascend/compiler/module.h>
32 #include <ascend/compiler/childinfo.h>
33 #include <ascend/compiler/parentchild.h>
34 #include <ascend/compiler/slist.h>
35 #include <ascend/compiler/type_desc.h>
36 #include <ascend/compiler/packages.h>
37 #include <ascend/compiler/symtab.h>
38 #include <ascend/compiler/instquery.h>
39 #include <ascend/compiler/instmacro.h>
40 #include <ascend/compiler/instance_types.h>
41
42 /* the code that we're wrapping... */
43 #include "fprops.h"
44 #include "sat.h"
45 #include "solve_ph.h"
46 #include "thcond.h"
47 #include "visc.h"
48
49 /* for the moment, species data are defined in C code, we'll implement something
50 better later on, hopefully. */
51 #include "fluids.h"
52
53 #ifndef ASC_EXPORT
54 # error "Where is ASC_EXPORT?"
55 #endif
56
57
58 /*------------------------------------------------------------------------------
59 FORWARD DECLARATIONS
60 */
61
62 ExtBBoxInitFunc asc_fprops_prepare;
63 ExtBBoxFunc fprops_p_calc;
64 ExtBBoxFunc fprops_u_calc;
65 ExtBBoxFunc fprops_s_calc;
66 ExtBBoxFunc fprops_h_calc;
67 ExtBBoxFunc fprops_a_calc;
68 ExtBBoxFunc fprops_g_calc;
69 ExtBBoxFunc fprops_cp_calc;
70 ExtBBoxFunc fprops_cv_calc;
71 ExtBBoxFunc fprops_w_calc;
72 ExtBBoxFunc fprops_mu_calc;
73 ExtBBoxFunc fprops_lam_calc;
74 ExtBBoxFunc fprops_phsx_vT_calc;
75 ExtBBoxFunc fprops_Tvsx_ph_calc;
76
77 #define TCRIT(FLUID) (FLUID->data->T_c)
78 #define TTRIP(FLUID) (FLUID->data->T_t)
79 #define RHOCRIT(FLUID) (FLUID->data->rho_c)
80 #define PCRIT(FLUID) (FLUID->data->p_c)
81
82 /*------------------------------------------------------------------------------
83 GLOBALS
84 */
85
86 /* place to store symbols needed for accessing ASCEND's instance tree */
87 static symchar *fprops_symbols[3];
88 #define COMPONENT_SYM fprops_symbols[0]
89 #define TYPE_SYM fprops_symbols[1]
90 #define SOURCE_SYM fprops_symbols[2]
91
92 static const char *fprops_p_help = "Calculate pressure from temperature and density, using FPROPS";
93 static const char *fprops_u_help = "Calculate specific internal energy from temperature and density, using FPROPS";
94 static const char *fprops_s_help = "Calculate specific entropy from temperature and density, using FPROPS";
95 static const char *fprops_h_help = "Calculate specific enthalpy from temperature and density, using FPROPS";
96 static const char *fprops_a_help = "Calculate specific Helmholtz energy from temperature and density, using FPROPS";
97 static const char *fprops_g_help = "Calculate specific Gibbs energy from temperature and density, using FPROPS";
98 static const char *fprops_cp_help = "Calculate isobaric specific heat from temperature and density, using FPROPS";
99 static const char *fprops_cv_help = "Calculate isochoric specific heat from temperature and density, using FPROPS";
100 static const char *fprops_w_help = "Calculate speed of sound from temperature and density, using FPROPS";
101 static const char *fprops_mu_help = "Calculate viscosity from temperature and density, using FPROPS";
102 static const char *fprops_lam_help = "Calculate thermal conductivity sound from temperature and density, using FPROPS";
103
104 static const char *fprops_phsx_vT_help = "Calculate p, h, s, x from temperature and density, using FPROPS/Helmholtz eqn";
105
106 static const char *fprops_Tvsx_ph_help = "Calculate T, v, s, x from pressure and enthalpy, using FPROPS/Helmholtz eqn";
107
108 /*------------------------------------------------------------------------------
109 REGISTRATION FUNCTION
110 */
111
112 /**
113 This is the function called from "IMPORT fprops"
114
115 It sets up the functions contained in this external library
116 */
117 extern
118 ASC_EXPORT int fprops_register(){
119 int result = 0;
120
121 ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.");
122
123 #define CALCFN(NAME,INPUTS,OUTPUTS) \
124 result += CreateUserFunctionBlackBox(#NAME \
125 , asc_fprops_prepare \
126 , NAME##_calc /* value */ \
127 , (ExtBBoxFunc*)NULL /* derivatives not provided yet*/ \
128 , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
129 , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
130 , INPUTS,OUTPUTS /* inputs, outputs */ \
131 , NAME##_help /* help text */ \
132 , 0.0 \
133 ) /* returns 0 on success */
134
135 #define CALCFNDERIV(NAME,INPUTS,OUTPUTS) \
136 result += CreateUserFunctionBlackBox(#NAME \
137 , asc_fprops_prepare \
138 , NAME##_calc /* value */ \
139 , NAME##_calc /* derivatives */ \
140 , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
141 , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
142 , INPUTS,OUTPUTS /* inputs, outputs */ \
143 , NAME##_help /* help text */ \
144 , 0.0 \
145 ) /* returns 0 on success */
146
147 CALCFNDERIV(fprops_p,2,1);
148 CALCFN(fprops_u,2,1);
149 CALCFN(fprops_s,2,1);
150 CALCFN(fprops_h,2,1);
151 CALCFN(fprops_a,2,1);
152 CALCFN(fprops_g,2,1);
153 CALCFN(fprops_cp,2,1);
154 CALCFN(fprops_cv,2,1);
155 CALCFN(fprops_w,2,1);
156 CALCFN(fprops_mu,2,1);
157 CALCFN(fprops_lam,2,1);
158 CALCFN(fprops_phsx_vT,2,4);
159 CALCFN(fprops_Tvsx_ph,2,4);
160
161 #undef CALCFN
162
163 if(result){
164 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
165 }
166 return result;
167 }
168
169 /**
170 'fprops_prepare' just gets the data member and checks that it's
171 valid, and stores it in the blackbox data field.
172 */
173 int asc_fprops_prepare(struct BBoxInterp *bbox,
174 struct Instance *data,
175 struct gl_list_t *arglist
176 ){
177 struct Instance *compinst, *typeinst, *srcinst;
178 const char *comp, *type = NULL, *src = NULL;
179
180 fprops_symbols[0] = AddSymbol("component");
181 fprops_symbols[1] = AddSymbol("type");
182 fprops_symbols[2] = AddSymbol("source");
183
184 /* get the component name */
185 compinst = ChildByChar(data,COMPONENT_SYM);
186 if(!compinst){
187 ERROR_REPORTER_HERE(ASC_USER_ERROR
188 ,"Couldn't locate 'component' in DATA, please check usage of FPROPS."
189 );
190 return 1;
191 }
192 if(InstanceKind(compinst)!=SYMBOL_CONSTANT_INST){
193 ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'component' must be a symbol_constant");
194 return 1;
195 }
196 comp = SCP(SYMC_INST(compinst)->value);
197 if(comp==NULL || strlen(comp)==0){
198 ERROR_REPORTER_HERE(ASC_USER_ERROR,"'component' is NULL or empty");
199 return 1;
200 }
201
202 /* get the component correlation type (FPROPS doesn't mind if none given) */
203 typeinst = ChildByChar(data,TYPE_SYM);
204 if(typeinst){
205 if(InstanceKind(typeinst)!=SYMBOL_CONSTANT_INST){
206 ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'type' must be a symbol_constant");
207 return 1;
208 }
209 type = SCP(SYMC_INST(typeinst)->value);
210 //CONSOLE_DEBUG("TYPE: %s",type?type:"(null)");
211 if(type && strlen(type)==0)type = NULL;
212 }
213
214 /* get the source data string (FPROPS doesn't mind if none given) */
215 srcinst = ChildByChar(data,SOURCE_SYM);
216 if(srcinst){
217 if(InstanceKind(srcinst)!=SYMBOL_CONSTANT_INST){
218 ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'source' must be a symbol_constant");
219 return 1;
220 }
221 src = SCP(SYMC_INST(srcinst)->value);
222 CONSOLE_DEBUG("SOURCE: %s",src?src:"(null)");
223 if(src && strlen(src)==0)src = NULL;
224 }
225
226 bbox->user_data = (void *)fprops_fluid(comp,type,src);
227 if(bbox->user_data == NULL){
228 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name/type was not recognised. Check the source-code for for the supported species.");
229 return 1;
230 }
231
232 #ifdef ASC_FPROPS_DEBUG
233 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s'%s%s%s OK.\n"
234 ,comp
235 ,type?" type '":""
236 ,type?type:""
237 ,type?"'":""
238 );
239 #endif
240 return 0;
241 }
242
243 /*------------------------------------------------------------------------------
244 EVALULATION ROUTINES
245 */
246
247 #define CALCPREPARE(NIN,NOUT) \
248 /* a few checks about the input requirements */ \
249 if(ninputs != NIN)return -1; \
250 if(noutputs != NOUT)return -2; \
251 if(inputs==NULL)return -3; \
252 if(outputs==NULL)return -4; \
253 if(bbox==NULL)return -5; \
254 \
255 /* the 'user_data' in the black box object will contain the */\
256 /* coefficients required for this fluid; cast it to the required form: */\
257 const PureFluid *FLUID = (const PureFluid *)bbox->user_data;\
258 FpropsError err=FPROPS_NO_ERROR;
259
260 /**
261 Evaluation function for 'fprops_p'.
262 @param inputs array with values of inputs T and rho.
263 @param outputs array with just value of p
264 @param jacobian, the partial derivative df/dx, where
265 each row is df[i]/dx[j] over each j for the y_out[i] of
266 matching index. The jacobian array is 1-D, row major, i.e.
267 df[i]/dx[j] -> jacobian[i*ninputs+j].
268 @return 0 on success
269 */
270 int fprops_p_calc(struct BBoxInterp *bbox,
271 int ninputs, int noutputs,
272 double *inputs, double *outputs,
273 double *jacobian
274 ){
275 CALCPREPARE(2,1);
276
277 /* first input is temperature, second is density */
278 if(bbox->task == bb_func_eval){
279 FluidState S = fprops_set_Trho(inputs[0],inputs[1], FLUID, &err);
280 outputs[0] = fprops_p(S, &err);
281 }else{
282 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
283 FluidState S = fprops_set_Trho(inputs[0],inputs[1], FLUID, &err);
284 jacobian[0*1+0] = fprops_dpdT_rho(S, &err);
285 jacobian[0*1+1] = fprops_dpdrho_T(S, &err);
286 }
287
288 /* no need to worry about error states etc. */
289 return 0;
290 }
291
292
293 /**
294 Evaluation function for 'fprops_u'
295 @param jacobian ignored
296 @return 0 on success
297 */
298 int fprops_u_calc(struct BBoxInterp *bbox,
299 int ninputs, int noutputs,
300 double *inputs, double *outputs,
301 double *jacobian
302 ){
303 CALCPREPARE(2,1);
304 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
305
306 /* first input is temperature, second is density */
307 if(bbox->task == bb_func_eval){
308 outputs[0] = fprops_u(S, &err);
309 }else{
310 jacobian[0*1+0] = fprops_dudT_rho(S, &err);
311 jacobian[0*1+1] = fprops_dudrho_T(S, &err);
312 }
313
314 /* no need to worry about error states etc. */
315 return 0;
316 }
317
318
319 /**
320 Evaluation function for 'fprops_s'
321 @param jacobian ignored
322 @return 0 on success
323 */
324 int fprops_s_calc(struct BBoxInterp *bbox,
325 int ninputs, int noutputs,
326 double *inputs, double *outputs,
327 double *jacobian
328 ){
329 CALCPREPARE(2,1);
330 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
331
332 /* first input is temperature, second is density */
333 outputs[0] = fprops_s(S, &err);
334
335 /* no need to worry about error states etc. */
336 return 0;
337 }
338
339
340 /**
341 Evaluation function for 'fprops_h'
342 @param jacobian ignored
343 @return 0 on success
344 */
345 int fprops_h_calc(struct BBoxInterp *bbox,
346 int ninputs, int noutputs,
347 double *inputs, double *outputs,
348 double *jacobian
349 ){
350 CALCPREPARE(2,1);
351 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
352
353 /* first input is temperature, second is density */
354 if(bbox->task == bb_func_eval){
355 outputs[0] = fprops_h(S, &err);
356 }else{
357 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
358 jacobian[0*1+0] = fprops_dhdT_rho(S, &err);
359 jacobian[0*1+1] = fprops_dhdrho_T(S, &err);
360 }
361
362 /* no need to worry about error states etc. */
363 return 0;
364 }
365
366
367 /**
368 Evaluation function for 'fprops_a'
369 @param jacobian ignored
370 @return 0 on success
371 */
372 int fprops_a_calc(struct BBoxInterp *bbox,
373 int ninputs, int noutputs,
374 double *inputs, double *outputs,
375 double *jacobian
376 ){
377 CALCPREPARE(2,1);
378 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
379
380 /* first input is temperature, second is density */
381 outputs[0] = fprops_a(S, &err);
382
383 /* no need to worry about error states etc. */
384 return 0;
385 }
386
387
388 /**
389 Evaluation function for 'fprops_g'
390 @param jacobian ignored
391 @return 0 on success
392 */
393 int fprops_g_calc(struct BBoxInterp *bbox,
394 int ninputs, int noutputs,
395 double *inputs, double *outputs,
396 double *jacobian
397 ){
398 CALCPREPARE(2,1);
399 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
400
401 /* first input is temperature, second is density */
402 outputs[0] = fprops_g(S, &err);
403
404 /* no need to worry about error states etc. */
405 return 0;
406 }
407
408
409 /**
410 Evaluation function for 'fprops_cp'
411 @param jacobian ignored
412 @return 0 on success
413 */
414 int fprops_cp_calc(struct BBoxInterp *bbox,
415 int ninputs, int noutputs,
416 double *inputs, double *outputs,
417 double *jacobian
418 ){
419 CALCPREPARE(2,1);
420 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
421
422 /* first input is temperature, second is density */
423 outputs[0] = fprops_cp(S, &err);
424
425 /* no need to worry about error states etc. */
426 return 0;
427 }
428
429
430 /**
431 Evaluation function for 'fprops_cv'
432 @param jacobian ignored
433 @return 0 on success
434 */
435 int fprops_cv_calc(struct BBoxInterp *bbox,
436 int ninputs, int noutputs,
437 double *inputs, double *outputs,
438 double *jacobian
439 ){
440 CALCPREPARE(2,1);
441 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
442
443 /* first input is temperature, second is density */
444 outputs[0] = fprops_cv(S, &err);
445
446 /* no need to worry about error states etc. */
447 return 0;
448 }
449
450
451 /**
452 Evaluation function for 'fprops_w'
453 @param jacobian ignored
454 @return 0 on success
455 */
456 int fprops_w_calc(struct BBoxInterp *bbox,
457 int ninputs, int noutputs,
458 double *inputs, double *outputs,
459 double *jacobian
460 ){
461 CALCPREPARE(2,1);
462 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
463
464 /* first input is temperature, second is density */
465 outputs[0] = fprops_w(S, &err);
466
467 /* no need to worry about error states etc. */
468 return 0;
469 }
470
471 /**
472 Evaluation function for 'fprops_mu'
473 @param jacobian ignored
474 @return 0 on success
475 */
476 int fprops_mu_calc(struct BBoxInterp *bbox,
477 int ninputs, int noutputs,
478 double *inputs, double *outputs,
479 double *jacobian
480 ){
481 CALCPREPARE(2,1);
482 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
483
484 /* first input is temperature, second is density */
485 outputs[0] = fprops_mu(S, &err);
486
487 /* no need to worry about error states etc. */
488 return 0;
489 }
490
491 /**
492 Evaluation function for 'fprops_lam'
493 @param jacobian ignored
494 @return 0 on success
495 */
496 int fprops_lam_calc(struct BBoxInterp *bbox,
497 int ninputs, int noutputs,
498 double *inputs, double *outputs,
499 double *jacobian
500 ){
501 CALCPREPARE(2,1);
502 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
503
504 /* first input is temperature, second is density */
505 outputs[0] = fprops_lam(S, &err);
506
507 /* no need to worry about error states etc. */
508 return 0;
509 }
510
511
512 /**
513 Evaluation function for 'fprops_phsx_vT'
514 @return 0 on success
515 */
516 int fprops_phsx_vT_calc(struct BBoxInterp *bbox,
517 int ninputs, int noutputs,
518 double *inputs, double *outputs,
519 double *jacobian
520 ){
521 CALCPREPARE(2,4);
522
523 double rho = 1./inputs[0];
524 double T = inputs[1];
525 double p_sat, rho_f, rho_g;
526
527 if(T < TCRIT(FLUID)){
528 fprops_sat_T(T, &p_sat, &rho_f, &rho_g, FLUID, &err);
529
530 if(rho < rho_f && rho > rho_g){
531 /* saturated */
532 double vf = 1./rho_f;
533 double vg = 1./rho_g;
534 double x = (inputs[0] - vf) /(vg - vf);
535 FluidState Sf = fprops_set_Trho(T, rho_f, FLUID, &err);
536 double sf = fprops_s(Sf, &err);
537 double hf = fprops_h(Sf, &err);
538 FluidState Sg = fprops_set_Trho(T, rho_g, FLUID, &err);
539 double sg = fprops_s(Sg, &err);
540 double hg = fprops_h(Sg, &err);
541 outputs[0] = p_sat;
542 outputs[1] = hf + x * (hg-hf);
543 outputs[2] = sf + x * (sg-sf);
544 outputs[3] = x;
545 /* maybe there was an error solving the saturation state? */
546 return err;
547 }
548 }
549
550 /* non-saturated */
551 FluidState S = fprops_set_Trho(T, rho, FLUID, &err);
552 outputs[0] = fprops_p(S, &err);
553 outputs[1] = fprops_h(S, &err);
554 outputs[2] = fprops_s(S, &err);
555 outputs[3] = rho < RHOCRIT(FLUID) ? 1 : 0;
556 return 0;
557 }
558
559
560 /**
561 Evaluation function for 'fprops_Tvsx_ph'
562 @return 0 on success
563 */
564 int fprops_Tvsx_ph_calc(struct BBoxInterp *bbox,
565 int ninputs, int noutputs,
566 double *inputs, double *outputs,
567 double *jacobian
568 ){
569 CALCPREPARE(2,4);
570
571 static const PureFluid *last = NULL;
572 static double p,h,T,v,s,x;
573 if(last == FLUID && p == inputs[0] && h == inputs[1]){
574 outputs[0] = T;
575 outputs[1] = v;
576 outputs[2] = s;
577 outputs[3] = x;
578 return 0;
579 }
580
581 p = inputs[0];
582 h = inputs[1];
583
584 double hft, pt, rhoft,rhogt;
585 fprops_triple_point(&pt,&rhoft,&rhogt,FLUID,&err);
586 if(err){
587 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve triple point for %s.",FLUID->name);
588 return 5;
589 }
590 FluidState Sft = fprops_set_Trho(TTRIP(FLUID),rhoft,FLUID,&err);
591 hft = fprops_h(Sft, &err);
592 if(h < hft){
593 ERROR_REPORTER_HERE(ASC_PROG_ERR
594 ,"Input enthalpy %f kJ/kg is below triple point liquid enthalpy %f kJ/kg"
595 ,h/1e3,hft/1e3
596 );
597 return 6;
598 }
599
600 if(p < pt){
601 ERROR_REPORTER_HERE(ASC_PROG_ERR
602 ,"Input pressure %f bar is below triple point pressure %f bar"
603 ,p/1e5,pt/1e5
604 );
605 outputs[0] = TTRIP(FLUID);
606 outputs[1] = 1./ rhoft;
607 outputs[2] = FLUID->s_fn(TTRIP(FLUID), rhoft, FLUID->data, &err);
608 outputs[3] = 0;
609 return 7;
610 }
611
612 if(p < PCRIT(FLUID)){
613 double T_sat, rho_f, rho_g;
614
615 fprops_sat_p(p, &T_sat, &rho_f, &rho_g, FLUID, &err);
616 if(err){
617 ERROR_REPORTER_HERE(ASC_PROG_ERR
618 , "Failed to solve saturation state of %s for p = %f bar < pc (= %f bar)"
619 , FLUID->name, p/1e5,PCRIT(FLUID)/1e5
620 );
621 outputs[0] = TTRIP(FLUID);
622 outputs[1] = 1./rhoft;
623 outputs[2] = FLUID->s_fn(TTRIP(FLUID), rhoft, FLUID->data, &err);
624 outputs[3] = 0;
625 return 8;
626 }
627
628 FluidState Sf = fprops_set_Trho(T_sat,rho_f,FLUID,&err);
629 FluidState Sg = fprops_set_Trho(T_sat,rho_g,FLUID,&err);
630 double hf = fprops_h(Sf, &err);
631 double hg = fprops_h(Sg, &err);
632
633 if(hf < h && h < hg){
634 /* saturated */
635 double vf = 1./rho_f;
636 double vg = 1./rho_g;
637 double sf = fprops_s(Sf, &err);
638 double sg = fprops_s(Sg, &err);
639 T = T_sat;
640 x = (h - hf) /(hg - hf);
641 v = vf + x * (vg-vf);
642 s = sf + x * (sg-sf);
643 last = FLUID;
644 outputs[0] = T;
645 outputs[1] = v;
646 outputs[2] = s;
647 outputs[3] = x;
648 #ifdef ASC_FPROPS_DEBUG
649 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Saturated state, p=%f bar, h = %f kJ/kg",p/1e5,h/1e3);
650 #endif
651 return 0;
652 }
653 }
654
655 double rho;
656 fprops_solve_ph(p,h, &T, &rho, 0, FLUID, &err);
657 if(err){
658 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for (p,h): %s",fprops_error(err));
659 return 9;
660 }
661 /* non-saturated */
662 v = 1./rho;
663 FluidState S = fprops_set_Trho(T,rho,FLUID,&err);
664 s = fprops_s(S, &err);
665 x = (v > 1./RHOCRIT(FLUID)) ? 1 : 0;
666 last = FLUID;
667 outputs[0] = T;
668 outputs[1] = v;
669 outputs[2] = s;
670 outputs[3] = x;
671 #ifdef ASC_FPROPS_DEBUG
672 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Non-saturated state, p = %f bar, h = %f kJ/kg",p/1e5,h/1e3);
673 #endif
674 return 0;
675 }
676
677
678

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