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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2734 - (show annotations) (download) (as text)
Tue Dec 10 12:11:46 2013 UTC (6 years, 7 months ago) by jpye
File MIME type: text/x-csrc
File size: 7463 byte(s)
testing ideal calcs, cf example in Bejan Tsatsaronis & Moran.
1 /* ASCEND modelling environment
2 Copyright (C) 2008-2013 John Pye
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 Ideal-gas components of helmholtz fundamental functions, calculated using
18 terms in cp0 in a standard power series form. For details see the
19 publications cited in the various fluid *.c files.
20
21 John Pye, Jul 2008.
22 */
23
24 #include <math.h>
25
26 #include "ideal.h"
27 #include "cp0.h"
28 #include "refstate.h"
29
30 #include <assert.h>
31 #include <stdlib.h>
32 #include <stdio.h>
33
34 #define SQ(X) ((X)*(X))
35
36 #ifndef FPROPS_NEW
37 # error "Where is FPROPS_NEW??"
38 #endif
39
40 #ifndef FPROPS_ARRAY_COPY
41 # error "where is FPROPS_ARRAY_COPY??"
42 #endif
43
44 //static void ideal_set_reference_std(FluidData *D, const ReferenceStateStd *R);
45
46 #define IDEAL_DEBUG
47 #ifdef IDEAL_DEBUG
48 # include "color.h"
49 # define MSG FPROPS_MSG
50 # define ERRMSG FPROPS_ERRMSG
51 #else
52 # define MSG(ARGS...) ((void)0)
53 # define ERRMSG(ARGS...) ((void)0)
54 #endif
55
56 /*--------------------------------------------
57 PREPARATION OF IDEAL RUNDATA from FILEDATA
58 */
59
60 PureFluid *ideal_prepare(const EosData *E, const ReferenceState *ref){
61 PureFluid *P = FPROPS_NEW(PureFluid);
62 P->data = FPROPS_NEW(FluidData);
63 #define D P->data
64
65 //MSG("...");
66
67 /* metadata */
68 P->name = E->name;
69 P->source = E->source;
70 P->type = FPROPS_IDEAL;
71
72 switch(E->type){
73 case FPROPS_CUBIC:
74 //MSG("Cubic");
75 D->M = E->data.cubic->M;
76 D->R = R_UNIVERSAL / D->M;
77 D->T_t = 0; /* TODO how will we flag this object so that sat.c doesn't try to solve? */
78 D->T_c = 0; /* TODO we need a temperature for scaling against, what should it be if critical point is not specific, and how will it be provided? */
79 D->p_c = 0;
80 D->rho_c = 0;
81 D->omega = 0;
82 D->Tstar = 1;
83 D->rhostar = 1;
84 D->cp0 = cp0_prepare(E->data.cubic->ideal, D->R, D->Tstar);
85 D->corr.helm = NULL;
86
87 //MSG("ref0 type = %d", E->data.cubic->ref0.type);
88 D->ref0 = E->data.cubic->ref0;
89 if(ref == NULL){
90 ref = &(E->data.cubic->ref);
91 }
92 break;
93 case FPROPS_HELMHOLTZ:
94 MSG("Helmholtz");
95 D->M = E->data.helm->M;
96 if(E->data.helm->R == 0){
97 P->data->R = R_UNIVERSAL / D->M;
98 }else{
99 P->data->R = E->data.helm->R;
100 }
101 D->T_t = 0;
102 D->T_c = E->data.helm->T_c;
103 D->p_c = 0;
104 D->rho_c = E->data.helm->rho_c;
105 D->omega = 0;
106 D->Tstar = 1;
107 D->rhostar = 1;
108 D->cp0 = cp0_prepare(E->data.helm->ideal, D->R, D->Tstar);
109 D->corr.helm = NULL;
110
111 if(ref == NULL){
112 ref = &(E->data.helm->ref);
113 }
114 break;
115 default:
116 ERRMSG("Unsupported source data type in ideal_prepare");
117 FPROPS_FREE(P->data);
118 FPROPS_FREE(P);
119 return NULL;
120 }
121
122 /* function pointers... more to come still? */
123 #define FN(VAR) P->VAR##_fn = &ideal_##VAR
124 FN(p); FN(u); FN(h); FN(s); FN(a); FN(g); FN(cp); FN(cv); FN(w);
125 FN(dpdrho_T);
126 FN(sat);
127 #undef FN
128 #undef D
129
130 //MSG("Setting reference state...");
131 // set the reference point
132 switch(ref->type){
133 case FPROPS_REF_PHI0:
134 MSG("Applying PHI0 reference data");
135 P->data->cp0->c = ref->data.phi0.c;
136 P->data->cp0->m = ref->data.phi0.m;
137 break;
138 case FPROPS_REF_REF0:
139 //MSG("Applying ref0 reference state");
140 switch(P->data->ref0.type){
141 case FPROPS_REF_TPHG:
142 {
143 //MSG("TPHG");
144 ReferenceState *ref0 = &(P->data->ref0);
145 //MSG("T0 = %f, p0 = %f, h0 = %f, g0 = %f",ref0->data.tphg.T0,ref0->data.tphg.p0,ref0->data.tphg.h0,ref0->data.tphg.g0);
146 FpropsError res = FPROPS_NO_ERROR;
147 double rho0 = ref0->data.tphg.p0 / P->data->R / ref0->data.tphg.T0;
148 double T0 = ref0->data.tphg.T0;
149 double s0 = (ref0->data.tphg.h0 - ref0->data.tphg.g0) / T0;
150 double h0 = ref0->data.tphg.h0;
151
152 P->data->cp0->c = 0;
153 P->data->cp0->m = 0;
154 //MSG("T0 = %f, rho0 = %f",T0,rho0);
155 //MSG("btw, p = %f", P->data->R * T0 *rho0); // is OK
156 res = FPROPS_NO_ERROR;
157 double h1 = ideal_h(T0, rho0, P->data, &res);
158 double s1 = ideal_s(T0, rho0, P->data, &res);
159 if(res)ERRMSG("error %d",res);
160 //MSG("h1 = %f",h1);
161 P->data->cp0->c = -(s0 - s1)/P->data->R;
162 P->data->cp0->m = (h0 - h1)/P->data->R/P->data->Tstar;
163
164 h0 = ideal_h(T0,rho0, P->data, &res);
165 if(res)ERRMSG("error %d",res);
166 //MSG("new h0(T0,rho0) = %f", h0);
167 double g0 = ideal_g(T0,rho0, P->data, &res);
168 if(res)ERRMSG("error %d",res);
169 //MSG("new g0(T0,rho0) = %f", g0);
170 //MSG("DONE");
171 }
172 break;
173 default:
174 ERRMSG("Unsupported type of reference state (ref0) in ideal_prepare");
175 FPROPS_FREE(P->data); FPROPS_FREE(P);
176 return NULL;
177 }
178 break;
179 default:
180 ERRMSG("Unsupported type of reference state requested in ideal_prepare.\n");
181 FPROPS_FREE(P->data);
182 FPROPS_FREE(P);
183 return NULL;
184 }
185
186 return P;
187 }
188
189 #define DEFINE_TAU double tau = data->Tstar / T
190 #define DEFINE_TAUDELTA DEFINE_TAU; double delta = rho / data->rhostar
191
192 double ideal_p(double T, double rho, const FluidData *data, FpropsError *err){
193 return data->R * T * rho;
194 }
195
196 double ideal_h(double T, double rho, const FluidData *data, FpropsError *err){
197 DEFINE_TAUDELTA;
198 return data->R * T * (1 + tau * ideal_phi_tau(tau,delta,data->cp0));
199 }
200
201 double ideal_s(double T, double rho, const FluidData *data, FpropsError *err){
202 DEFINE_TAUDELTA;
203 double pht = ideal_phi_tau(tau,delta,data->cp0);
204 double ph = ideal_phi(tau,delta,data->cp0);
205 return data->R * (tau * ideal_phi_tau(tau,delta,data->cp0) - ideal_phi(tau,delta,data->cp0));
206 }
207
208 double ideal_u(double T, double rho, const FluidData *data, FpropsError *err){
209 return ideal_h(T,rho,data,err) - data->R * T;
210 }
211
212 double ideal_a(double T, double rho, const FluidData *data, FpropsError *err){
213 return ideal_h(T,rho,data,err) - T * (data->R + ideal_s(T,rho,data,err));
214 }
215
216 double ideal_g(double T, double rho, const FluidData *data, FpropsError *err){
217 //MSG("g(T=%f,rho=%f)...",T,rho);
218 double h = ideal_h(T,rho,data,err);
219 double s = ideal_s(T,rho,data,err);
220 //MSG("h = %f, T = %f, s = %f, h-T*s = %f",h,T,s,h-T*s);
221 return h - T * s;
222 }
223
224 /**
225 Note that this function is called by ALL fluid types via 'fprops_cp0' which
226 means that it needs to include the scaling temperature within the structure;
227 we can't just define Tstar as a constant for ideal fluids.
228 */
229 double ideal_cp(double T, double rho, const FluidData *data, FpropsError *err){
230 DEFINE_TAU;
231 double res = data->R * (1. - SQ(tau) * ideal_phi_tautau(tau,data->cp0));
232 return res;
233 }
234
235 double ideal_cv(double T, double rho, const FluidData *data, FpropsError *err){
236 DEFINE_TAU;
237 return - data->R * SQ(tau) * ideal_phi_tautau(tau,data->cp0);
238 }
239
240 double ideal_w(double T, double rho, const FluidData *data, FpropsError *err){
241 DEFINE_TAU;
242 double w2onRT = 1. - 1. / (SQ(tau) * ideal_phi_tautau(tau,data->cp0));
243 return sqrt(data->R * T * w2onRT);
244 }
245
246 double ideal_dpdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
247 return data->R * T;
248 }
249
250 double ideal_sat(double T,double *rhof_ret, double *rhog_ret, const FluidData *data, FpropsError *err){
251 MSG("Ideal gas: saturation calculation is not possible");
252 *err = FPROPS_RANGE_ERROR;
253 return 0;
254 }
255
256
257
258

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