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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2776 - (hide annotations) (download) (as text)
Thu Mar 27 12:48:14 2014 UTC (6 years, 6 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 jpye 2654 /* 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 jpye 2661 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 jpye 2654 *//** @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 jpye 2776 #include "thcond.h"
47     #include "visc.h"
48 jpye 2654
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 jpye 2776 ExtBBoxFunc fprops_mu_calc;
73     ExtBBoxFunc fprops_lam_calc;
74 jpye 2654 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 jpye 2663 static symchar *fprops_symbols[3];
88 jpye 2654 #define COMPONENT_SYM fprops_symbols[0]
89     #define TYPE_SYM fprops_symbols[1]
90 jpye 2663 #define SOURCE_SYM fprops_symbols[2]
91 jpye 2654
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 jpye 2776 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 jpye 2654
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 jpye 2776 CALCFN(fprops_mu,2,1);
157     CALCFN(fprops_lam,2,1);
158 jpye 2654 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 jpye 2663 struct Instance *compinst, *typeinst, *srcinst;
178     const char *comp, *type = NULL, *src = NULL;
179 jpye 2654
180     fprops_symbols[0] = AddSymbol("component");
181     fprops_symbols[1] = AddSymbol("type");
182 jpye 2663 fprops_symbols[2] = AddSymbol("source");
183 jpye 2654
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 jpye 2689 //CONSOLE_DEBUG("TYPE: %s",type?type:"(null)");
211 jpye 2657 if(type && strlen(type)==0)type = NULL;
212 jpye 2654 }
213    
214 jpye 2663 /* 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 jpye 2654 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 jpye 2776 /**
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 jpye 2654
482 jpye 2776 /* 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 jpye 2654 /**
490 jpye 2776 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 jpye 2654 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 jpye 2660 return 7;
608 jpye 2654 }
609    
610     if(p < PCRIT(FLUID)){
611     double T_sat, rho_f, rho_g;
612 jpye 2660
613 jpye 2654 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 jpye 2660 return 8;
624 jpye 2654 }
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 jpye 2660 if(err){
656     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve for (p,h): %s",fprops_error(err));
657     return 9;
658     }
659 jpye 2654 /* 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 jpye 2660 return 0;
673 jpye 2654 }
674    
675    
676    

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