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

Contents of /branches/pallav/models/johnpye/fprops/cp0.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2823 - (show annotations) (download) (as text)
Tue Feb 17 06:38:36 2015 UTC (4 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 8386 byte(s)
restored files from non-svn backup

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 void cp0_destroy(Phi0RunData *N){
147 if(N->pt)FPROPS_FREE(N->pt);
148 if(N->et)FPROPS_FREE(N->et);
149 FPROPS_FREE(N);
150 }
151
152 /*---------------------------------------------
153 IDEAL COMPONENT RELATIONS
154 */
155
156 /*
157 Hypothesis:
158 in calculating the ideal component relations, we have some integration
159 constants that ultimately allow the arbitrary scales of h and s to
160 be specified. Hence we can ignore components of helm_ideal that
161 are either constant or linear in tau, and work them out later when we
162 stick in the values of data->c and data->m.
163 */
164
165 //#define IDEAL_DEBUG
166 /**
167 Ideal component of helmholtz function
168 */
169 double ideal_phi(double tau, double delta, const Phi0RunData *data){
170
171 const Phi0RunPowTerm *pt;
172 const Phi0RunExpTerm *et;
173
174 unsigned i;
175 // FIXME what if rhostar != rhoc??
176 double sum = log(delta) + data->c + data->m * tau;
177 double term;
178
179 #ifdef IDEAL_DEBUG
180 fprintf(stderr,"\ttau = %f, delta = %f\n",tau,delta);
181 fprintf(stderr,"\tT ~ %f\n\n",Tstar_on_tau);
182 fprintf(stderr,"\tlog(delta) + c + m tau = %f (c=%f,m=%f)\n",sum,data->c, data->m);
183 fprintf(stderr,"sum = %f\n",sum);
184 #endif
185
186 /* power terms */
187 pt = &(data->pt[0]);
188 for(i = 0; i<data->np; ++i, ++pt){
189 double a = pt->a;
190 double p = pt->p;
191 if(p == 0){
192 term = a*log(tau);
193 #ifdef IDEAL_DEBUG
194 fprintf(stderr,"\ta log(tau) = %f (a=%f, t=%f)\n",term,a,p);
195 #endif
196 }else{
197 #ifdef IDEAL_DEBUG
198 assert(p!=-1);
199 #endif
200 term = a * pow(tau, p);
201 #ifdef IDEAL_DEBUG
202 fprintf(stderr,"\ta tau^p = %f (a=%f, p=%f)\n",term,a,p);
203 #endif
204 }
205 sum += term;
206 #ifdef IDEAL_DEBUG
207 fprintf(stderr,"sum = %f\n",sum);
208 #endif
209 }
210
211 /* Planck-Einstein terms */
212 et = &(data->et[0]);
213 for(i=0; i<data->ne; ++i, ++et){
214 term = et->n * log(1 - exp(-et->gamma * tau));
215 #ifdef IDEAL_DEBUG
216 fprintf(stderr,"\tn log(1-exp(-gamma*tau)) = %f, (n=%f, gamma=%f)\n",term,et->n,et->gamma);
217 #endif
218 sum += term;
219 }
220
221 #ifdef IDEAL_DEBUG
222 fprintf(stderr,"phi0 = %f\n",sum);
223 #endif
224
225 return sum;
226 }
227
228 /**
229 Partial dervivative of ideal component (phi0) of normalised helmholtz
230 residual function (phi), with respect to tau.
231 */
232 double ideal_phi_tau(double tau, double delta, const Phi0RunData *data){
233 const Phi0RunPowTerm *pt;
234 const Phi0RunExpTerm *et;
235
236 unsigned i;
237 double term;
238 double sum = data->m;
239
240 pt = &(data->pt[0]);
241 for(i = 0; i<data->np; ++i, ++pt){
242 double a = pt->a;
243 double p = pt->p;
244 if(p == 0){
245 term = a / tau;
246 //fprintf(stderr,"\tc/tau = %f\n",term);
247 }else{
248 term = a*p*pow(tau,p - 1);
249 //fprintf(stderr,"\tc / tau^p = %f (t=%f, c=%.3e, p=%f)\n",term,t,coeff,-t-1.);
250 }
251 #ifdef TEST
252 if(isinf(term)){
253 fprintf(stderr,"Error with infinite-valued term with i = %d, a = %f, p = %f\n", i,a ,p);
254 abort();
255 }
256 #endif
257 assert(!isnan(term));
258 sum += term;
259 }
260
261 /* Planck-Einstein terms */
262 et = &(data->et[0]);
263 for(i=0; i<data->ne; ++i, ++et){
264 double expo = exp(-et->gamma * tau);
265 term = et->n * et->gamma * expo / (1 - expo);
266 #ifdef TEST
267 assert(!isinf(term));
268 #endif
269 sum += term;
270 }
271
272 #ifdef TEST
273 assert(!isinf(sum));
274 #endif
275 return sum;
276 }
277
278
279
280 /**
281 Second partial dervivative of ideal component (phi0) of normalised helmholtz
282 residual function (phi), with respect to tau. This one is easy!
283 It's not a function of delta.
284
285 FIXME: although this one is easy, we want to pull Tstar out of the
286 ideal properties stuff, if that's possible.
287 */
288 double ideal_phi_tautau(double tau, const Phi0RunData *data){
289 const Phi0RunPowTerm *pt;
290 const Phi0RunExpTerm *et;
291
292 unsigned i;
293
294 double sum = 0;
295 double term;
296
297 #ifdef CP0_DEBUG
298 fprintf(stderr,"\ttau = %f\n",tau);
299 #endif
300
301 /* power terms */
302 pt = &(data->pt[0]);
303 for(i = 0; i<data->np; ++i, ++pt){
304 if(pt->p == 0){
305 term = pt->a;
306 }else{
307 term = -pt->a * pt->p * (pt->p - 1) * pow(tau, pt->p);
308 }
309 #ifdef CP0_DEBUG
310 fprintf(stderr,"\tpt[%d] = ap(p-1)*tau^p (a = %e, p = %e) = %f\n",i,pt->a, pt->p, term);
311 #endif
312 sum += term;
313 }
314
315 /* 'exponential' terms */
316 et = &(data->et[0]);
317 for(i=0; i<data->ne; ++i, ++et){
318 double x = et->gamma * tau;
319 double e = exp(-x);
320 double d = (1-e)*(1-e);
321 term = et->n * x*x * e / d;
322 #ifdef CP0_DEBUG
323 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);
324 #endif
325 sum += term;
326 }
327 /* note, at this point, sum == cp0/R - 1 */
328 MSG("sum = %f, phi_tautau = %f",sum, -sum/SQ(tau));
329 return -sum/SQ(tau);
330 }

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