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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2290 - (hide annotations) (download) (as text)
Sun Aug 15 10:11:54 2010 UTC (13 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 18639 byte(s)
Trying to catch errors associated with pressure below triple point for CO2.
1 jpye 1822 /* 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, write to the Free Software
16     Foundation, Inc., 59 Temple Place - Suite 330,
17     Boston, MA 02111-1307, USA.
18     *//** @file
19     Wrapper for helmholtz.c to allow access to the routine from ASCEND.
20     */
21    
22     /* include the external function API from libascend... */
23 jpye 2018 #include <ascend/compiler/extfunc.h>
24 jpye 1822
25     /* include error reporting API as well, so we can send messages to user */
26 jpye 2018 #include <ascend/utilities/error.h>
27 jpye 1822
28     /* for accessing the DATA instance */
29 jpye 2018 #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 jpye 1822
42     /* the code that we're wrapping... */
43     #include "helmholtz.h"
44 jpye 2262 #include "sat.h"
45 jpye 2272 #include "solve_ph.h"
46 jpye 1822
47 jpye 1832 /* for the moment, species data are defined in C code, we'll implement something
48     better later on, hopefully. */
49     #include "ammonia.h"
50 jpye 1849 #include "nitrogen.h"
51 jpye 1995 #include "hydrogen.h"
52     #include "water.h"
53 jpye 2110 #include "carbondioxide.h"
54 hongke 2280 #include "methane.h"
55     #include "carbonmonoxide.h"
56     #include "ethanol.h"
57     #include "acetone.h"
58     #include "carbonylsulfide.h"
59     #include "decane.h"
60     #include "hydrogensulfide.h"
61     #include "isohexane.h"
62     #include "isopentane.h"
63     #include "krypton.h"
64     #include "neopentane.h"
65     #include "nitrousoxide.h"
66     #include "nonane.h"
67     #include "sulfurdioxide.h"
68     #include "toluene.h"
69     #include "xenon.h"
70     #include "butane.h"
71     #include "butene.h"
72     #include "cisbutene.h"
73     #include "isobutene.h"
74     #include "transbutene.h"
75     #include "dimethylether.h"
76     #include "ethane.h"
77     #include "parahydrogen.h"
78     #include "isobutane.h"
79     #include "r41.h"
80     #include "r116.h"
81     #include "r141b.h"
82     #include "r142b.h"
83     #include "r218.h"
84     #include "r245fa.h"
85 jpye 1832
86 jpye 1822 #ifndef ASC_EXPORT
87     # error "Where is ASC_EXPORT?"
88     #endif
89    
90    
91     /*------------------------------------------------------------------------------
92     FORWARD DECLARATIONS
93     */
94    
95 jpye 1825 ExtBBoxInitFunc helmholtz_prepare;
96 jpye 1822 ExtBBoxFunc helmholtz_p_calc;
97 jpye 1825 ExtBBoxFunc helmholtz_u_calc;
98 jpye 1903 ExtBBoxFunc helmholtz_s_calc;
99 jpye 1825 ExtBBoxFunc helmholtz_h_calc;
100 jpye 1903 ExtBBoxFunc helmholtz_a_calc;
101 jpye 2124 ExtBBoxFunc helmholtz_g_calc;
102 jpye 2236 ExtBBoxFunc helmholtz_phsx_vT_calc;
103 jpye 2272 ExtBBoxFunc helmholtz_Tvsx_ph_calc;
104 jpye 1822
105     /*------------------------------------------------------------------------------
106     GLOBALS
107     */
108    
109     /* place to store symbols needed for accessing ASCEND's instance tree */
110     static symchar *helmholtz_symbols[1];
111     #define COMPONENT_SYM helmholtz_symbols[0]
112    
113 jpye 1825 static const char *helmholtz_p_help = "Calculate pressure from temperature and density, using Helmholtz fundamental correlation";
114     static const char *helmholtz_u_help = "Calculate specific internal energy from temperature and density, using Helmholtz fundamental correlation";
115 jpye 1903 static const char *helmholtz_s_help = "Calculate specific entropy from temperature and density, using Helmholtz fundamental correlation";
116 jpye 1825 static const char *helmholtz_h_help = "Calculate specific enthalpy from temperature and density, using Helmholtz fundamental correlation";
117 jpye 1903 static const char *helmholtz_a_help = "Calculate specific Helmholtz energy from temperature and density, using Helmholtz fundamental correlation";
118 jpye 2124 static const char *helmholtz_g_help = "Calculate specific Gibbs energy from temperature and density, using Helmholtz fundamental correlation";
119 jpye 1825
120 jpye 2236 static const char *helmholtz_phsx_vT_help = "Calculate p, h, s, x from temperature and density, using FPROPS/Helmholtz eqn";
121    
122 jpye 2272 static const char *helmholtz_Tvsx_ph_help = "Calculate T, v, s, x from pressure and enthalpy, using FPROPS/Helmholtz eqn";
123 jpye 2124
124 jpye 2272
125 jpye 1822 /*------------------------------------------------------------------------------
126     REGISTRATION FUNCTION
127     */
128    
129     /**
130     This is the function called from "IMPORT helmholtz"
131    
132     It sets up the functions contained in this external library
133     */
134     extern
135     ASC_EXPORT int helmholtz_register(){
136     int result = 0;
137    
138 jpye 1920 ERROR_REPORTER_HERE(ASC_USER_WARNING,"FPROPS is still EXPERIMENTAL. Use with caution.\n");
139 jpye 1822
140 jpye 1825 #define CALCFN(NAME,INPUTS,OUTPUTS) \
141     result += CreateUserFunctionBlackBox(#NAME \
142     , helmholtz_prepare \
143     , NAME##_calc /* value */ \
144     , (ExtBBoxFunc*)NULL /* derivatives not provided yet*/ \
145     , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
146     , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
147     , INPUTS,OUTPUTS /* inputs, outputs */ \
148     , NAME##_help /* help text */ \
149     , 0.0 \
150     ) /* returns 0 on success */
151 jpye 1822
152 jpye 1909 #define CALCFNDERIV(NAME,INPUTS,OUTPUTS) \
153     result += CreateUserFunctionBlackBox(#NAME \
154     , helmholtz_prepare \
155     , NAME##_calc /* value */ \
156     , NAME##_calc /* derivatives */ \
157     , (ExtBBoxFunc*)NULL /* hessian not provided yet */ \
158     , (ExtBBoxFinalFunc*)NULL /* finalisation not implemented */ \
159     , INPUTS,OUTPUTS /* inputs, outputs */ \
160     , NAME##_help /* help text */ \
161     , 0.0 \
162     ) /* returns 0 on success */
163    
164     CALCFNDERIV(helmholtz_p,2,1);
165 jpye 1825 CALCFN(helmholtz_u,2,1);
166 jpye 1903 CALCFN(helmholtz_s,2,1);
167 jpye 1825 CALCFN(helmholtz_h,2,1);
168 jpye 1903 CALCFN(helmholtz_a,2,1);
169 jpye 2124 CALCFN(helmholtz_g,2,1);
170 jpye 2236 CALCFN(helmholtz_phsx_vT,2,4);
171 jpye 2272 CALCFN(helmholtz_Tvsx_ph,2,4);
172 jpye 1825
173     #undef CALCFN
174    
175 jpye 1822 if(result){
176     ERROR_REPORTER_HERE(ASC_PROG_NOTE,"CreateUserFunction result = %d\n",result);
177     }
178     return result;
179     }
180    
181 jpye 1825 /**
182     'helmholtz_prepare' just gets the data member and checks that it's
183     valid, and stores it in the blackbox data field.
184 jpye 1822 */
185 jpye 1825 int helmholtz_prepare(struct BBoxInterp *bbox,
186 jpye 1822 struct Instance *data,
187     struct gl_list_t *arglist
188     ){
189     struct Instance *compinst;
190     const char *comp;
191    
192     helmholtz_symbols[0] = AddSymbol("component");
193    
194     /* get the data file name (we will look for this file in the ASCENDLIBRARY path) */
195     compinst = ChildByChar(data,COMPONENT_SYM);
196     if(!compinst){
197     ERROR_REPORTER_HERE(ASC_USER_ERROR
198     ,"Couldn't locate 'component' in DATA, please check usage of HELMHOLTZ."
199     );
200     return 1;
201     }
202     if(InstanceKind(compinst)!=SYMBOL_CONSTANT_INST){
203     ERROR_REPORTER_HERE(ASC_USER_ERROR,"DATA member 'component' must be a symbol_constant");
204     return 1;
205     }
206     comp = SCP(SYMC_INST(compinst)->value);
207     CONSOLE_DEBUG("COMPONENT: %s",comp);
208     if(comp==NULL || strlen(comp)==0){
209     ERROR_REPORTER_HERE(ASC_USER_ERROR,"'component' is NULL or empty");
210     return 1;
211     }
212    
213 jpye 1849 if(strcmp(comp,"ammonia")==0){
214     bbox->user_data = (void*)&helmholtz_data_ammonia;
215     }else if(strcmp(comp,"nitrogen")==0){
216     bbox->user_data = (void*)&helmholtz_data_nitrogen;
217 jpye 1995 }else if(strcmp(comp,"hydrogen")==0){
218     bbox->user_data = (void*)&helmholtz_data_hydrogen;
219     }else if(strcmp(comp,"water")==0){
220     bbox->user_data = (void*)&helmholtz_data_water;
221 jpye 2110 }else if(strcmp(comp,"carbondioxide")==0){
222     bbox->user_data = (void*)&helmholtz_data_carbondioxide;
223 hongke 2280 }else if(strcmp(comp,"methane")==0){
224     bbox->user_data = (void*)&helmholtz_data_methane;
225     }else if(strcmp(comp,"carbonmonoxide")==0){
226     bbox->user_data = (void*)&helmholtz_data_carbonmonoxide;
227     }else if(strcmp(comp,"ethanol")==0){
228     bbox->user_data = (void*)&helmholtz_data_ethanol;
229     }else if(strcmp(comp,"acetone")==0){
230     bbox->user_data = (void*)&helmholtz_data_acetone;
231     }else if(strcmp(comp,"carbonylsulfide")==0){
232     bbox->user_data = (void*)&helmholtz_data_carbonylsulfide;
233     }else if(strcmp(comp,"decane")==0){
234     bbox->user_data = (void*)&helmholtz_data_decane;
235     }else if(strcmp(comp,"hydrogensulfide")==0){
236     bbox->user_data = (void*)&helmholtz_data_hydrogensulfide;
237     }else if(strcmp(comp,"isohexane")==0){
238     bbox->user_data = (void*)&helmholtz_data_isohexane;
239     }else if(strcmp(comp,"isopentane")==0){
240     bbox->user_data = (void*)&helmholtz_data_isopentane;
241     }else if(strcmp(comp,"krypton")==0){
242     bbox->user_data = (void*)&helmholtz_data_krypton;
243     }else if(strcmp(comp,"neopentane")==0){
244     bbox->user_data = (void*)&helmholtz_data_neopentane;
245     }else if(strcmp(comp,"nitrousoxide")==0){
246     bbox->user_data = (void*)&helmholtz_data_nitrousoxide;
247     }else if(strcmp(comp,"nonane")==0){
248     bbox->user_data = (void*)&helmholtz_data_nonane;
249     }else if(strcmp(comp,"sulfurdioxide")==0){
250     bbox->user_data = (void*)&helmholtz_data_sulfurdioxide;
251     }else if(strcmp(comp,"toluene")==0){
252     bbox->user_data = (void*)&helmholtz_data_toluene;
253     }else if(strcmp(comp,"xenon")==0){
254     bbox->user_data = (void*)&helmholtz_data_xenon;
255     }else if(strcmp(comp,"butane")==0){
256     bbox->user_data = (void*)&helmholtz_data_butane;
257     }else if(strcmp(comp,"butene")==0){
258     bbox->user_data = (void*)&helmholtz_data_butene;
259     }else if(strcmp(comp,"cisbutene")==0){
260     bbox->user_data = (void*)&helmholtz_data_cisbutene;
261     }else if(strcmp(comp,"isobutene")==0){
262     bbox->user_data = (void*)&helmholtz_data_isobutene;
263     }else if(strcmp(comp,"transbutene")==0){
264     bbox->user_data = (void*)&helmholtz_data_transbutene;
265     }else if(strcmp(comp,"dimethylether")==0){
266     bbox->user_data = (void*)&helmholtz_data_dimethylether;
267     }else if(strcmp(comp,"ethane")==0){
268     bbox->user_data = (void*)&helmholtz_data_ethane;
269     }else if(strcmp(comp,"parahydrogen")==0){
270     bbox->user_data = (void*)&helmholtz_data_parahydrogen;
271     }else if(strcmp(comp,"isobutane")==0){
272     bbox->user_data = (void*)&helmholtz_data_isobutane;
273     }else if(strcmp(comp,"r41")==0){
274     bbox->user_data = (void*)&helmholtz_data_r41;
275     }else if(strcmp(comp,"r116")==0){
276     bbox->user_data = (void*)&helmholtz_data_r116;
277     }else if(strcmp(comp,"r141b")==0){
278     bbox->user_data = (void*)&helmholtz_data_r141b;
279     }else if(strcmp(comp,"r142b")==0){
280     bbox->user_data = (void*)&helmholtz_data_r142b;
281     }else if(strcmp(comp,"r218")==0){
282     bbox->user_data = (void*)&helmholtz_data_r218;
283     }else if(strcmp(comp,"r245fa")==0){
284     bbox->user_data = (void*)&helmholtz_data_r245fa;
285     }else{
286 jpye 1995 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Component name was not recognised. Check the source-code for for the supported species.");
287 jpye 1822 }
288    
289 jpye 1849 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Prepared component '%s' OK.\n",comp);
290 jpye 1822 return 0;
291     }
292    
293 jpye 1825 /*------------------------------------------------------------------------------
294     EVALULATION ROUTINES
295     */
296    
297 jpye 2124 #define CALCPREPARE(NIN,NOUT) \
298 jpye 1825 /* a few checks about the input requirements */ \
299 jpye 2124 if(ninputs != NIN)return -1; \
300     if(noutputs != NOUT)return -2; \
301 jpye 1825 if(inputs==NULL)return -3; \
302     if(outputs==NULL)return -4; \
303     if(bbox==NULL)return -5; \
304     \
305     /* the 'user_data' in the black box object will contain the */\
306     /* coefficients required for this fluid; cast it to the required form: */\
307 jpye 2272 const HelmholtzData *helmholtz_data = (const HelmholtzData *)bbox->user_data
308 jpye 1825
309 jpye 1822 /**
310 jpye 1825 Evaluation function for 'helmholtz_p'.
311 jpye 1909 @param inputs array with values of inputs T and rho.
312     @param outputs array with just value of p
313     @param jacobian, the partial derivative df/dx, where
314     each row is df[i]/dx[j] over each j for the y_out[i] of
315     matching index. The jacobian array is 1-D, row major, i.e.
316     df[i]/dx[j] -> jacobian[i*ninputs+j].
317 jpye 1822 @return 0 on success
318     */
319     int helmholtz_p_calc(struct BBoxInterp *bbox,
320     int ninputs, int noutputs,
321     double *inputs, double *outputs,
322     double *jacobian
323     ){
324 jpye 2124 CALCPREPARE(2,1);
325 jpye 1822
326 jpye 2236 /* first input is temperature, second is density */
327 jpye 1909 if(bbox->task == bb_func_eval){
328     outputs[0] = helmholtz_p(inputs[0], inputs[1], helmholtz_data);
329     }else{
330     //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
331     jacobian[0*1+0] = helmholtz_dpdT_rho(inputs[0], inputs[1], helmholtz_data);
332     jacobian[0*1+1] = helmholtz_dpdrho_T(inputs[0], inputs[1], helmholtz_data);
333     }
334 jpye 1822
335     /* no need to worry about error states etc. */
336 jpye 1825 return 0;
337     }
338 jpye 1822
339 jpye 1825
340     /**
341     Evaluation function for 'helmholtz_u'
342     @param jacobian ignored
343     @return 0 on success
344     */
345     int helmholtz_u_calc(struct BBoxInterp *bbox,
346     int ninputs, int noutputs,
347     double *inputs, double *outputs,
348     double *jacobian
349     ){
350 jpye 2124 CALCPREPARE(2,1);
351 jpye 1825
352 jpye 2236 /* first input is temperature, second is density */
353 jpye 1923 if(bbox->task == bb_func_eval){
354     outputs[0] = helmholtz_u(inputs[0], inputs[1], helmholtz_data);
355     }else{
356     jacobian[0*1+0] = helmholtz_dudT_rho(inputs[0], inputs[1], helmholtz_data);
357     jacobian[0*1+1] = helmholtz_dudrho_T(inputs[0], inputs[1], helmholtz_data);
358     }
359 jpye 1825
360     /* no need to worry about error states etc. */
361 jpye 1822 return 0;
362     }
363 jpye 1825
364    
365     /**
366 jpye 2124 Evaluation function for 'helmholtz_s'
367 jpye 1825 @param jacobian ignored
368     @return 0 on success
369     */
370 jpye 1903 int helmholtz_s_calc(struct BBoxInterp *bbox,
371     int ninputs, int noutputs,
372     double *inputs, double *outputs,
373     double *jacobian
374     ){
375 jpye 2124 CALCPREPARE(2,1);
376 jpye 1903
377 jpye 2236 /* first input is temperature, second is density */
378 jpye 1903 outputs[0] = helmholtz_s(inputs[0], inputs[1], helmholtz_data);
379    
380     /* no need to worry about error states etc. */
381     return 0;
382     }
383    
384    
385     /**
386     Evaluation function for 'helmholtz_h'
387     @param jacobian ignored
388     @return 0 on success
389     */
390 jpye 1825 int helmholtz_h_calc(struct BBoxInterp *bbox,
391     int ninputs, int noutputs,
392     double *inputs, double *outputs,
393     double *jacobian
394     ){
395 jpye 2124 CALCPREPARE(2,1);
396 jpye 1825
397 jpye 2236 /* first input is temperature, second is density */
398 jpye 1920 if(bbox->task == bb_func_eval){
399     outputs[0] = helmholtz_h(inputs[0], inputs[1], helmholtz_data);
400     }else{
401     //ERROR_REPORTER_HERE(ASC_USER_NOTE,"JACOBIAN CALCULATION FOR P!\n");
402     jacobian[0*1+0] = helmholtz_dhdT_rho(inputs[0], inputs[1], helmholtz_data);
403     jacobian[0*1+1] = helmholtz_dhdrho_T(inputs[0], inputs[1], helmholtz_data);
404     }
405 jpye 1825
406     /* no need to worry about error states etc. */
407     return 0;
408     }
409    
410    
411 jpye 1903 /**
412 jpye 2124 Evaluation function for 'helmholtz_a'
413 jpye 1903 @param jacobian ignored
414     @return 0 on success
415     */
416     int helmholtz_a_calc(struct BBoxInterp *bbox,
417     int ninputs, int noutputs,
418     double *inputs, double *outputs,
419     double *jacobian
420     ){
421 jpye 2124 CALCPREPARE(2,1);
422 jpye 1825
423 jpye 2236 /* first input is temperature, second is density */
424 jpye 1903 outputs[0] = helmholtz_a(inputs[0], inputs[1], helmholtz_data);
425 jpye 1825
426 jpye 1903 /* no need to worry about error states etc. */
427     return 0;
428     }
429    
430    
431 jpye 2124 /**
432     Evaluation function for 'helmholtz_g'
433     @param jacobian ignored
434     @return 0 on success
435     */
436     int helmholtz_g_calc(struct BBoxInterp *bbox,
437     int ninputs, int noutputs,
438     double *inputs, double *outputs,
439     double *jacobian
440     ){
441     CALCPREPARE(2,1);
442 jpye 1903
443 jpye 2236 /* first input is temperature, second is density */
444 jpye 2124 outputs[0] = helmholtz_g(inputs[0], inputs[1], helmholtz_data);
445    
446     /* no need to worry about error states etc. */
447     return 0;
448     }
449    
450    
451 jpye 2236
452 jpye 2124 /**
453 jpye 2236 Evaluation function for 'helmholtz_phsx_vT'
454     @return 0 on success
455     */
456     int helmholtz_phsx_vT_calc(struct BBoxInterp *bbox,
457     int ninputs, int noutputs,
458     double *inputs, double *outputs,
459     double *jacobian
460     ){
461     CALCPREPARE(2,4);
462    
463     double rho = 1./inputs[0];
464     double T = inputs[1];
465     double p_sat, rho_f, rho_g;
466    
467     if(T < helmholtz_data->T_c){
468     int res = fprops_sat_T(T, &p_sat, &rho_f, &rho_g, helmholtz_data);
469    
470     if(rho < rho_f && rho > rho_g){
471     /* saturated */
472     double vf = 1./rho_f;
473     double vg = 1./rho_g;
474     double x = (inputs[0] - vf) /(vg - vf);
475     double sf = helmholtz_s(T,rho_f, helmholtz_data);
476     double hf = helmholtz_h(T,rho_f, helmholtz_data);
477     double sg = helmholtz_s(T,rho_g, helmholtz_data);
478     double hg = helmholtz_h(T,rho_g, helmholtz_data);
479     outputs[0] = p_sat;
480     outputs[1] = hf + x * (hg-hf);
481     outputs[2] = sf + x * (sg-sf);
482     outputs[3] = x;
483     /* maybe there was an error solving the saturation state? */
484     return res;
485     }
486     }
487    
488     /* non-saturated */
489     outputs[0] = helmholtz_p(T,rho, helmholtz_data);
490     outputs[1] = helmholtz_h(T,rho, helmholtz_data);
491     outputs[2] = helmholtz_s(T,rho, helmholtz_data);
492 jpye 2274 outputs[3] = rho < helmholtz_data->rho_c ? 1 : 0;
493 jpye 2236 return 0;
494     }
495    
496    
497     /**
498 jpye 2272 Evaluation function for 'helmholtz_Tvsx_ph'
499 jpye 2124 @return 0 on success
500     */
501 jpye 2272 int helmholtz_Tvsx_ph_calc(struct BBoxInterp *bbox,
502 jpye 2124 int ninputs, int noutputs,
503     double *inputs, double *outputs,
504     double *jacobian
505     ){
506 jpye 2272 CALCPREPARE(2,4);
507 jpye 2124
508 jpye 2272 static const HelmholtzData *last = NULL;
509     static double p,h,T,v,s,x;
510     int res;
511     if(last == helmholtz_data && p == inputs[0] && h == inputs[1]){
512     outputs[0] = T;
513     outputs[1] = v;
514     outputs[2] = s;
515     outputs[3] = x;
516     return 0;
517     }
518 jpye 2290
519 jpye 2272 p = inputs[0];
520     h = inputs[1];
521 jpye 2124
522 jpye 2290 double hft, pt, rhoft,rhogt;
523     res = fprops_triple_point(&pt,&rhoft,&rhogt,helmholtz_data);
524     if(res){
525     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to solve triple point for %s.",helmholtz_data->name);
526     return 5;
527     }
528     hft = helmholtz_h(helmholtz_data->T_t, rhoft, helmholtz_data);
529     if(h < hft){
530     ERROR_REPORTER_HERE(ASC_PROG_ERR
531     ,"Input enthalpy %f kJ/kg is below triple point liquid enthalpy %f kJ/kg"
532     ,h/1e3,hft/1e3
533     );
534     return 6;
535     }
536    
537     if(p < pt){
538     ERROR_REPORTER_HERE(ASC_PROG_ERR
539     ,"Input pressure %f bar is below triple point pressure %f bar"
540     ,p/1e5,pt/1e5
541     );
542     outputs[0] = helmholtz_data->T_t;
543     outputs[1] = 1./ rhoft;
544     outputs[2] = helmholtz_s_raw(helmholtz_data->T_t, rhoft, helmholtz_data);
545     outputs[3] = 0;
546     return 6;
547     }
548    
549 jpye 2272 if(p < fprops_pc(helmholtz_data)){
550     double T_sat, rho_f, rho_g;
551     res = fprops_sat_p(p, &T_sat, &rho_f, &rho_g, helmholtz_data);
552 jpye 2290 if(res){
553     ERROR_REPORTER_HERE(ASC_PROG_ERR
554     , "Failed to solve saturation state of %s for p = %f bar < pc (= %f bar)"
555     , helmholtz_data->name, p/1e5,fprops_pc(helmholtz_data)/1e5
556     );
557     outputs[0] = helmholtz_data->T_t;
558     outputs[1] = 1./rhoft;
559     outputs[2] = helmholtz_s_raw(helmholtz_data->T_t, rhoft, helmholtz_data);
560     outputs[3] = 0;
561     return 1;
562     }
563 jpye 2272 double hf = helmholtz_h(T_sat,rho_f, helmholtz_data);
564     double hg = helmholtz_h(T_sat,rho_g, helmholtz_data);
565 jpye 2124
566 jpye 2272 if(hf < h && h < hg){
567     /* saturated */
568     double vf = 1./rho_f;
569     double vg = 1./rho_g;
570     double sf = helmholtz_s(T_sat,rho_f, helmholtz_data);
571     double sg = helmholtz_s(T_sat,rho_g, helmholtz_data);
572     T = T_sat;
573     x = (h - hf) /(hg - hf);
574     v = vf + x * (vg-vf);
575     s = sf + x * (sg-sf);
576     last = helmholtz_data;
577     outputs[0] = T;
578     outputs[1] = v;
579     outputs[2] = s;
580     outputs[3] = x;
581 jpye 2290 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Saturated state, p=%f bar, h = %f kJ/kg",p/1e5,h/1e3);
582 jpye 2272 return 0;
583     }
584     }
585    
586     double rho;
587     res = fprops_solve_ph(p,h, &T, &rho, 0, helmholtz_data);
588     /* non-saturated */
589     v = 1./rho;
590     s = helmholtz_s(T,rho, helmholtz_data);
591     x = (v > 1./helmholtz_data->rho_c) ? 1 : 0;
592     last = helmholtz_data;
593     outputs[0] = T;
594     outputs[1] = v;
595     outputs[2] = s;
596     outputs[3] = x;
597 jpye 2290 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Non-saturated state, p = %f bar, h = %f kJ/kg",p/1e5,h/1e3);
598 jpye 2275 return res;
599 jpye 2124 }
600    
601 jpye 2272
602    

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