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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2671 - (show annotations) (download) (as text)
Thu Jan 24 00:35:44 2013 UTC (6 years, 7 months ago) by jpye
File MIME type: text/x-csrc
File size: 8272 byte(s)
Added testing of saturation curve for all fluids in their 'best' correlation. Six fluids are showing occasional failures, all of them PR EOS.
1 /* ASCEND modelling environment
2 Copyright (C) 2008-2009 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 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 "cp0.h"
27
28 #include <assert.h>
29 #include <stdlib.h>
30 #include <stdio.h>
31
32 #define SQ(X) ((X)*(X))
33
34 //#define CP0_DEBUG
35 #ifdef CP0_DEBUG
36 # include "color.h"
37 # define MSG(FMT, ...) \
38 color_on(stderr,ASC_FG_BRIGHTRED);\
39 fprintf(stderr,"%s:%d: ",__FILE__,__LINE__);\
40 color_on(stderr,ASC_FG_BRIGHTBLUE);\
41 fprintf(stderr,"%s: ",__func__);\
42 color_off(stderr);\
43 fprintf(stderr,FMT "\n",##__VA_ARGS__)
44 #else
45 # define MSG(ARGS...) ((void)0)
46 #endif
47
48 /*--------------------------------------------
49 PREPARATION OF IDEAL RUNDATA from FILEDATA
50 */
51
52 /*
53 FIXME
54
55 we have changed the definition of the cp0 expression to cp/R = sum( c *T/T*)^t )
56
57 but now we need to fix the re-derive the conversion from the cp0 to alpha0
58 expressions
59
60 currently there is an error!
61 */
62 Phi0RunData *cp0_prepare(const IdealData *I, double R, double Tstar){
63 Phi0RunData *N = FPROPS_NEW(Phi0RunData);
64 int i, add_const_term=1;
65 double Tred, cp0red, p;
66 N->c = 0;
67 N->m = 0;
68 switch(I->type){
69 case IDEAL_CP0:
70 N->np = I->data.cp0.np;
71 N->ne = I->data.cp0.ne;
72 Tred = I->data.cp0.Tstar;
73 cp0red = I->data.cp0.cp0star;
74 MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);
75 MSG("Tred = %f, Tstar = %f", Tred, Tstar);
76 MSG("cp0red = %f, R = %f", cp0red, R);
77
78 for(i=0; i < N->np; ++i)if(I->data.cp0.pt[i].t == 0)add_const_term = 0;
79 MSG("add_const_term = %d",add_const_term);
80
81 N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np + add_const_term);
82 N->et = FPROPS_NEW_ARRAY(Phi0RunExpTerm, N->ne);
83 add_const_term = 1;
84
85 for(i=0; i < N->np; ++i){
86 p = - I->data.cp0.pt[i].t;
87 N->pt[i].p = p;
88 if(p == 0){
89 //if(add_const_term)fprintf(stderr,"Addding offset here\n");
90 N->pt[i].a = (I->data.cp0.pt[i].c - add_const_term);
91 add_const_term = 0;
92 }else{
93 N->pt[i].a = -I->data.cp0.pt[i].c / p / (p - 1) * pow(Tred / Tstar, p);
94 }
95 }
96
97 if(add_const_term){
98 MSG("WARNING: adding constant term %d in cp0, is that what you really want?",i);
99 N->pt[i].a = -1;
100 N->pt[i].p = 0;
101 N->np++;
102 }
103
104 for(i=0; i < N->ne; ++i){
105 N->et[i].n = I->data.cp0.et[i].b;
106 N->et[i].gamma = I->data.cp0.et[i].beta * Tred / Tstar;
107 }
108
109 if(cp0red != R){
110 MSG("WARNING: adjusting for R (=%f) != cp0red (=%f)...\n",R,cp0red);
111 double X = cp0red / R;
112 // scale for any differences in R and cpstar */
113 for(i=0; i < N->np; ++i)N->pt[i].a *= X;
114 for(i=0; i < N->ne; ++i)N->et[i].n *= X;
115 }
116
117 break;
118 case IDEAL_PHI0:
119 // TODO add checks for disallowed terms, eg p = 0 or p = 1?
120 N->np = I->data.phi0.np;
121 N->ne = I->data.phi0.ne;
122 MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);
123
124 // power terms
125 N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np);
126 for(i=0; i< N->np; ++i){
127 double a = I->data.phi0.pt[i].a0;
128 double p = I->data.phi0.pt[i].p0;
129 N->pt[i].a = a;
130 N->pt[i].p = p;
131 }
132
133 // exponential terms
134 N->et = FPROPS_NEW_ARRAY(Phi0RunExpTerm,N->ne);
135 for(i=0; i < N->ne; ++i){
136 double n = I->data.phi0.et[i].n;
137 double gamma = I->data.phi0.et[i].gamma;
138 N->et[i].gamma = gamma;
139 N->et[i].n = n;
140 }
141 break;
142 }
143 return N;
144 }
145
146 /*---------------------------------------------
147 IDEAL COMPONENT RELATIONS
148 */
149
150 /*
151 Hypothesis:
152 in calculating the ideal component relations, we have some integration
153 constants that ultimately allow the arbitrary scales of h and s to
154 be specified. Hence we can ignore components of helm_ideal that
155 are either constant or linear in tau, and work them out later when we
156 stick in the values of data->c and data->m.
157 */
158
159 //#define IDEAL_DEBUG
160 /**
161 Ideal component of helmholtz function
162 */
163 double ideal_phi(double tau, double delta, const Phi0RunData *data){
164
165 const Phi0RunPowTerm *pt;
166 const Phi0RunExpTerm *et;
167
168 unsigned i;
169 // FIXME what if rhostar != rhoc??
170 double sum = log(delta) + data->c + data->m * tau;
171 double term;
172
173 #ifdef IDEAL_DEBUG
174 fprintf(stderr,"\ttau = %f, delta = %f\n",tau,delta);
175 fprintf(stderr,"\tT ~ %f\n\n",Tstar_on_tau);
176 fprintf(stderr,"\tlog(delta) + c + m tau = %f (c=%f,m=%f)\n",sum,data->c, data->m);
177 fprintf(stderr,"sum = %f\n",sum);
178 #endif
179
180 /* power terms */
181 pt = &(data->pt[0]);
182 for(i = 0; i<data->np; ++i, ++pt){
183 double a = pt->a;
184 double p = pt->p;
185 if(p == 0){
186 term = a*log(tau);
187 #ifdef IDEAL_DEBUG
188 fprintf(stderr,"\ta log(tau) = %f (a=%f, t=%f)\n",term,a,p);
189 #endif
190 }else{
191 #ifdef IDEAL_DEBUG
192 assert(p!=-1);
193 #endif
194 term = a * pow(tau, p);
195 #ifdef IDEAL_DEBUG
196 fprintf(stderr,"\ta tau^p = %f (a=%f, p=%f)\n",term,a,p);
197 #endif
198 }
199 sum += term;
200 #ifdef IDEAL_DEBUG
201 fprintf(stderr,"sum = %f\n",sum);
202 #endif
203 }
204
205 /* Planck-Einstein terms */
206 et = &(data->et[0]);
207 for(i=0; i<data->ne; ++i, ++et){
208 term = et->n * log(1 - exp(-et->gamma * tau));
209 #ifdef IDEAL_DEBUG
210 fprintf(stderr,"\tn log(1-exp(-gamma*tau)) = %f, (n=%f, gamma=%f)\n",term,et->n,et->gamma);
211 #endif
212 sum += term;
213 }
214
215 #ifdef IDEAL_DEBUG
216 fprintf(stderr,"phi0 = %f\n",sum);
217 #endif
218
219 return sum;
220 }
221
222 /**
223 Partial dervivative of ideal component (phi0) of normalised helmholtz
224 residual function (phi), with respect to tau.
225 */
226 double ideal_phi_tau(double tau, double delta, const Phi0RunData *data){
227 const Phi0RunPowTerm *pt;
228 const Phi0RunExpTerm *et;
229
230 unsigned i;
231 double term;
232 double sum = data->m;
233
234 pt = &(data->pt[0]);
235 for(i = 0; i<data->np; ++i, ++pt){
236 double a = pt->a;
237 double p = pt->p;
238 if(p == 0){
239 term = a / tau;
240 //fprintf(stderr,"\tc/tau = %f\n",term);
241 }else{
242 term = a*p*pow(tau,p - 1);
243 //fprintf(stderr,"\tc / tau^p = %f (t=%f, c=%.3e, p=%f)\n",term,t,coeff,-t-1.);
244 }
245 #ifdef TEST
246 if(isinf(term)){
247 fprintf(stderr,"Error with infinite-valued term with i = %d, a = %f, p = %f\n", i,a ,p);
248 abort();
249 }
250 #endif
251 assert(!isnan(term));
252 sum += term;
253 }
254
255 /* Planck-Einstein terms */
256 et = &(data->et[0]);
257 for(i=0; i<data->ne; ++i, ++et){
258 double expo = exp(-et->gamma * tau);
259 term = et->n * et->gamma * expo / (1 - expo);
260 #ifdef TEST
261 assert(!isinf(term));
262 #endif
263 sum += term;
264 }
265
266 #ifdef TEST
267 assert(!isinf(sum));
268 #endif
269 return sum;
270 }
271
272
273
274 /**
275 Second partial dervivative of ideal component (phi0) of normalised helmholtz
276 residual function (phi), with respect to tau. This one is easy!
277 It's not a function of delta.
278
279 FIXME: although this one is easy, we want to pull Tstar out of the
280 ideal properties stuff, if that's possible.
281 */
282 double ideal_phi_tautau(double tau, const Phi0RunData *data){
283 const Phi0RunPowTerm *pt;
284 const Phi0RunExpTerm *et;
285
286 unsigned i;
287
288 double sum = 0;
289 double term;
290
291 #ifdef CP0_DEBUG
292 fprintf(stderr,"\ttau = %f\n",tau);
293 #endif
294
295 /* power terms */
296 pt = &(data->pt[0]);
297 for(i = 0; i<data->np; ++i, ++pt){
298 if(pt->p == 0){
299 term = pt->a;
300 }else{
301 term = -pt->a * pt->p * (pt->p - 1) * pow(tau, pt->p);
302 }
303 #ifdef CP0_DEBUG
304 fprintf(stderr,"\tpt[%d] = ap(p-1)*tau^p (a = %e, p = %e) = %f\n",i,pt->a, pt->p, term);
305 #endif
306 sum += term;
307 }
308
309 /* 'exponential' terms */
310 et = &(data->et[0]);
311 for(i=0; i<data->ne; ++i, ++et){
312 double x = et->gamma * tau;
313 double e = exp(-x);
314 double d = (1-e)*(1-e);
315 term = et->n * x*x * e / d;
316 #ifdef CP0_DEBUG
317 fprintf(stderr,"\tet[%d] = n x^2 exp(-x)/(1 - exp(-x))^2 (n = %e, x=gamma*tau, gamma = %e) = %f\n",i,et->n, et->gamma, term);
318 #endif
319 sum += term;
320 }
321 /* note, at this point, sum == cp0/R - 1 */
322 MSG("sum = %f, phi_tautau = %f",sum, -sum/SQ(tau));
323 return -sum/SQ(tau);
324 }

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