/[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 2776 - (show annotations) (download) (as text)
Thu Mar 27 12:48:14 2014 UTC (6 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 18839 byte(s)
add thcond,visc to ascend wrapper (only CO2 implemented/working at this stage)
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.\n");
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 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s'%s%s%s OK.\n"
233 ,comp
234 ,type?" type '":""
235 ,type?type:""
236 ,type?"'":""
237 );
238 return 0;
239 }
240
241 /*------------------------------------------------------------------------------
242 EVALULATION ROUTINES
243 */
244
245 #define CALCPREPARE(NIN,NOUT) \
246 /* a few checks about the input requirements */ \
247 if(ninputs != NIN)return -1; \
248 if(noutputs != NOUT)return -2; \
249 if(inputs==NULL)return -3; \
250 if(outputs==NULL)return -4; \
251 if(bbox==NULL)return -5; \
252 \
253 /* the 'user_data' in the black box object will contain the */\
254 /* coefficients required for this fluid; cast it to the required form: */\
255 const PureFluid *FLUID = (const PureFluid *)bbox->user_data;\
256 FpropsError err=FPROPS_NO_ERROR;
257
258 /**
259 Evaluation function for 'fprops_p'.
260 @param inputs array with values of inputs T and rho.
261 @param outputs array with just value of p
262 @param jacobian, the partial derivative df/dx, where
263 each row is df[i]/dx[j] over each j for the y_out[i] of
264 matching index. The jacobian array is 1-D, row major, i.e.
265 df[i]/dx[j] -> jacobian[i*ninputs+j].
266 @return 0 on success
267 */
268 int fprops_p_calc(struct BBoxInterp *bbox,
269 int ninputs, int noutputs,
270 double *inputs, double *outputs,
271 double *jacobian
272 ){
273 CALCPREPARE(2,1);
274
275 /* first input is temperature, second is density */
276 if(bbox->task == bb_func_eval){
277 FluidState S = fprops_set_Trho(inputs[0],inputs[1], FLUID, &err);
278 outputs[0] = fprops_p(S, &err);
279 }else{
280 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
281 FluidState S = fprops_set_Trho(inputs[0],inputs[1], FLUID, &err);
282 jacobian[0*1+0] = fprops_dpdT_rho(S, &err);
283 jacobian[0*1+1] = fprops_dpdrho_T(S, &err);
284 }
285
286 /* no need to worry about error states etc. */
287 return 0;
288 }
289
290
291 /**
292 Evaluation function for 'fprops_u'
293 @param jacobian ignored
294 @return 0 on success
295 */
296 int fprops_u_calc(struct BBoxInterp *bbox,
297 int ninputs, int noutputs,
298 double *inputs, double *outputs,
299 double *jacobian
300 ){
301 CALCPREPARE(2,1);
302 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
303
304 /* first input is temperature, second is density */
305 if(bbox->task == bb_func_eval){
306 outputs[0] = fprops_u(S, &err);
307 }else{
308 jacobian[0*1+0] = fprops_dudT_rho(S, &err);
309 jacobian[0*1+1] = fprops_dudrho_T(S, &err);
310 }
311
312 /* no need to worry about error states etc. */
313 return 0;
314 }
315
316
317 /**
318 Evaluation function for 'fprops_s'
319 @param jacobian ignored
320 @return 0 on success
321 */
322 int fprops_s_calc(struct BBoxInterp *bbox,
323 int ninputs, int noutputs,
324 double *inputs, double *outputs,
325 double *jacobian
326 ){
327 CALCPREPARE(2,1);
328 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
329
330 /* first input is temperature, second is density */
331 outputs[0] = fprops_s(S, &err);
332
333 /* no need to worry about error states etc. */
334 return 0;
335 }
336
337
338 /**
339 Evaluation function for 'fprops_h'
340 @param jacobian ignored
341 @return 0 on success
342 */
343 int fprops_h_calc(struct BBoxInterp *bbox,
344 int ninputs, int noutputs,
345 double *inputs, double *outputs,
346 double *jacobian
347 ){
348 CALCPREPARE(2,1);
349 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
350
351 /* first input is temperature, second is density */
352 if(bbox->task == bb_func_eval){
353 outputs[0] = fprops_h(S, &err);
354 }else{
355 //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
356 jacobian[0*1+0] = fprops_dhdT_rho(S, &err);
357 jacobian[0*1+1] = fprops_dhdrho_T(S, &err);
358 }
359
360 /* no need to worry about error states etc. */
361 return 0;
362 }
363
364
365 /**
366 Evaluation function for 'fprops_a'
367 @param jacobian ignored
368 @return 0 on success
369 */
370 int fprops_a_calc(struct BBoxInterp *bbox,
371 int ninputs, int noutputs,
372 double *inputs, double *outputs,
373 double *jacobian
374 ){
375 CALCPREPARE(2,1);
376 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
377
378 /* first input is temperature, second is density */
379 outputs[0] = fprops_a(S, &err);
380
381 /* no need to worry about error states etc. */
382 return 0;
383 }
384
385
386 /**
387 Evaluation function for 'fprops_g'
388 @param jacobian ignored
389 @return 0 on success
390 */
391 int fprops_g_calc(struct BBoxInterp *bbox,
392 int ninputs, int noutputs,
393 double *inputs, double *outputs,
394 double *jacobian
395 ){
396 CALCPREPARE(2,1);
397 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
398
399 /* first input is temperature, second is density */
400 outputs[0] = fprops_g(S, &err);
401
402 /* no need to worry about error states etc. */
403 return 0;
404 }
405
406
407 /**
408 Evaluation function for 'fprops_cp'
409 @param jacobian ignored
410 @return 0 on success
411 */
412 int fprops_cp_calc(struct BBoxInterp *bbox,
413 int ninputs, int noutputs,
414 double *inputs, double *outputs,
415 double *jacobian
416 ){
417 CALCPREPARE(2,1);
418 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
419
420 /* first input is temperature, second is density */
421 outputs[0] = fprops_cp(S, &err);
422
423 /* no need to worry about error states etc. */
424 return 0;
425 }
426
427
428 /**
429 Evaluation function for 'fprops_cv'
430 @param jacobian ignored
431 @return 0 on success
432 */
433 int fprops_cv_calc(struct BBoxInterp *bbox,
434 int ninputs, int noutputs,
435 double *inputs, double *outputs,
436 double *jacobian
437 ){
438 CALCPREPARE(2,1);
439 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
440
441 /* first input is temperature, second is density */
442 outputs[0] = fprops_cv(S, &err);
443
444 /* no need to worry about error states etc. */
445 return 0;
446 }
447
448
449 /**
450 Evaluation function for 'fprops_w'
451 @param jacobian ignored
452 @return 0 on success
453 */
454 int fprops_w_calc(struct BBoxInterp *bbox,
455 int ninputs, int noutputs,
456 double *inputs, double *outputs,
457 double *jacobian
458 ){
459 CALCPREPARE(2,1);
460 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
461
462 /* first input is temperature, second is density */
463 outputs[0] = fprops_w(S, &err);
464
465 /* no need to worry about error states etc. */
466 return 0;
467 }
468
469 /**
470 Evaluation function for 'fprops_mu'
471 @param jacobian ignored
472 @return 0 on success
473 */
474 int fprops_mu_calc(struct BBoxInterp *bbox,
475 int ninputs, int noutputs,
476 double *inputs, double *outputs,
477 double *jacobian
478 ){
479 CALCPREPARE(2,1);
480 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
481
482 /* first input is temperature, second is density */
483 outputs[0] = fprops_mu(S, &err);
484
485 /* no need to worry about error states etc. */
486 return 0;
487 }
488
489 /**
490 Evaluation function for 'fprops_lam'
491 @param jacobian ignored
492 @return 0 on success
493 */
494 int fprops_lam_calc(struct BBoxInterp *bbox,
495 int ninputs, int noutputs,
496 double *inputs, double *outputs,
497 double *jacobian
498 ){
499 CALCPREPARE(2,1);
500 FluidState S = fprops_set_Trho(inputs[0], inputs[1], FLUID, &err);
501
502 /* first input is temperature, second is density */
503 outputs[0] = fprops_lam(S, &err);
504
505 /* no need to worry about error states etc. */
506 return 0;
507 }
508
509
510 /**
511 Evaluation function for 'fprops_phsx_vT'
512 @return 0 on success
513 */
514 int fprops_phsx_vT_calc(struct BBoxInterp *bbox,
515 int ninputs, int noutputs,
516 double *inputs, double *outputs,
517 double *jacobian
518 ){
519 CALCPREPARE(2,4);
520
521 double rho = 1./inputs[0];
522 double T = inputs[1];
523 double p_sat, rho_f, rho_g;
524
525 if(T < TCRIT(FLUID)){
526 fprops_sat_T(T, &p_sat, &rho_f, &rho_g, FLUID, &err);
527
528 if(rho < rho_f && rho > rho_g){
529 /* saturated */
530 double vf = 1./rho_f;
531 double vg = 1./rho_g;
532 double x = (inputs[0] - vf) /(vg - vf);
533 FluidState Sf = fprops_set_Trho(T, rho_f, FLUID, &err);
534 double sf = fprops_s(Sf, &err);
535 double hf = fprops_h(Sf, &err);
536 FluidState Sg = fprops_set_Trho(T, rho_g, FLUID, &err);
537 double sg = fprops_s(Sg, &err);
538 double hg = fprops_h(Sg, &err);
539 outputs[0] = p_sat;
540 outputs[1] = hf + x * (hg-hf);
541 outputs[2] = sf + x * (sg-sf);
542 outputs[3] = x;
543 /* maybe there was an error solving the saturation state? */
544 return err;
545 }
546 }
547
548 /* non-saturated */
549 FluidState S = fprops_set_Trho(T, rho, FLUID, &err);
550 outputs[0] = fprops_p(S, &err);
551 outputs[1] = fprops_h(S, &err);
552 outputs[2] = fprops_s(S, &err);
553 outputs[3] = rho < RHOCRIT(FLUID) ? 1 : 0;
554 return 0;
555 }
556
557
558 /**
559 Evaluation function for 'fprops_Tvsx_ph'
560 @return 0 on success
561 */
562 int fprops_Tvsx_ph_calc(struct BBoxInterp *bbox,
563 int ninputs, int noutputs,
564 double *inputs, double *outputs,
565 double *jacobian
566 ){
567 CALCPREPARE(2,4);
568
569 static const PureFluid *last = NULL;
570 static double p,h,T,v,s,x;
571 if(last == FLUID && p == inputs[0] && h == inputs[1]){
572 outputs[0] = T;
573 outputs[1] = v;
574 outputs[2] = s;
575 outputs[3] = x;
576 return 0;
577 }
578
579 p = inputs[0];
580 h = inputs[1];
581
582 double hft, pt, rhoft,rhogt;
583 fprops_triple_point(&pt,&rhoft,&rhogt,FLUID,&err);
584 if(err){
585 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve triple point for %s.",FLUID->name);
586 return 5;
587 }
588 FluidState Sft = fprops_set_Trho(TTRIP(FLUID),rhoft,FLUID,&err);
589 hft = fprops_h(Sft, &err);
590 if(h < hft){
591 ERROR_REPORTER_HERE(ASC_PROG_ERR
592 ,"Input enthalpy %f kJ/kg is below triple point liquid enthalpy %f kJ/kg"
593 ,h/1e3,hft/1e3
594 );
595 return 6;
596 }
597
598 if(p < pt){
599 ERROR_REPORTER_HERE(ASC_PROG_ERR
600 ,"Input pressure %f bar is below triple point pressure %f bar"
601 ,p/1e5,pt/1e5
602 );
603 outputs[0] = TTRIP(FLUID);
604 outputs[1] = 1./ rhoft;
605 outputs[2] = FLUID->s_fn(TTRIP(FLUID), rhoft, FLUID->data, &err);
606 outputs[3] = 0;
607 return 7;
608 }
609
610 if(p < PCRIT(FLUID)){
611 double T_sat, rho_f, rho_g;
612
613 fprops_sat_p(p, &T_sat, &rho_f, &rho_g, FLUID, &err);
614 if(err){
615 ERROR_REPORTER_HERE(ASC_PROG_ERR
616 , "Failed to solve saturation state of %s for p = %f bar < pc (= %f bar)"
617 , FLUID->name, p/1e5,PCRIT(FLUID)/1e5
618 );
619 outputs[0] = TTRIP(FLUID);
620 outputs[1] = 1./rhoft;
621 outputs[2] = FLUID->s_fn(TTRIP(FLUID), rhoft, FLUID->data, &err);
622 outputs[3] = 0;
623 return 8;
624 }
625
626 FluidState Sf = fprops_set_Trho(T_sat,rho_f,FLUID,&err);
627 FluidState Sg = fprops_set_Trho(T_sat,rho_g,FLUID,&err);
628 double hf = fprops_h(Sf, &err);
629 double hg = fprops_h(Sg, &err);
630
631 if(hf < h && h < hg){
632 /* saturated */
633 double vf = 1./rho_f;
634 double vg = 1./rho_g;
635 double sf = fprops_s(Sf, &err);
636 double sg = fprops_s(Sg, &err);
637 T = T_sat;
638 x = (h - hf) /(hg - hf);
639 v = vf + x * (vg-vf);
640 s = sf + x * (sg-sf);
641 last = FLUID;
642 outputs[0] = T;
643 outputs[1] = v;
644 outputs[2] = s;
645 outputs[3] = x;
646 #ifdef ASC_FPROPS_DEBUG
647 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Saturated state, p=%f bar, h = %f kJ/kg",p/1e5,h/1e3);
648 #endif
649 return 0;
650 }
651 }
652
653 double rho;
654 fprops_solve_ph(p,h, &T, &rho, 0, FLUID, &err);
655 if(err){
656 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for (p,h): %s",fprops_error(err));
657 return 9;
658 }
659 /* non-saturated */
660 v = 1./rho;
661 FluidState S = fprops_set_Trho(T,rho,FLUID,&err);
662 s = fprops_s(S, &err);
663 x = (v > 1./RHOCRIT(FLUID)) ? 1 : 0;
664 last = FLUID;
665 outputs[0] = T;
666 outputs[1] = v;
667 outputs[2] = s;
668 outputs[3] = x;
669 #ifdef ASC_FPROPS_DEBUG
670 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Non-saturated state, p = %f bar, h = %f kJ/kg",p/1e5,h/1e3);
671 #endif
672 return 0;
673 }
674
675
676

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