/[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 2739 - (show annotations) (download) (as text)
Thu Dec 12 12:33:26 2013 UTC (6 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 9405 byte(s)
cp0 problem appears to be when cp0red==1 in cp0_prepare with IDEAL_CP0 data type.
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 FPROPS_MSG
38 #else
39 # define MSG(ARGS...) ((void)0)
40 #endif
41
42 /*--------------------------------------------
43 PREPARATION OF IDEAL RUNDATA from FILEDATA
44 */
45
46 /**
47 This function prepares the ideal part of the helmholtz function, phi (\$f \phi = \frac{a}{R T}\f$).
48 If we have IdealData of type IDEAL_PHI0, then we just copy the data directly.
49 We assume that we don't need to normalise for Tstar or R in that case (it's
50 not even checked).
51
52 If we have IdealData of type IDEAL_CP0, then we have to do some conversions
53 in order to be able to obtain the required internal $\f\ \phi(\tau,\delta) \$f
54 function required here. The complication is that in general, cp0 can be
55 expressed with different normalising parameter cp0star and Tstar than those
56 used for the residual part function. So in that case, we normalise to R and
57 Tstar which are provided as parameters to this function. Except it doesn't \
58 work correctly for Peng-Robinson or Ideal EOS yet.
59
60 FIXME
61
62 we have changed the definition of the cp0 expression to cp/R = sum( c *T/T*)^t )
63
64 but now we need to fix the re-derive the conversion from the cp0 to alpha0
65 expressions
66
67 currently there is an error!
68
69 TODO check if ^^^ is still true
70 */
71 Phi0RunData *cp0_prepare(const IdealData *I, double R, double Tstar){
72 MSG("Tstar = %f",Tstar);
73 #ifdef CP0_DEBUG
74 if(I->type==IDEAL_CP0){
75 MSG("R=%f,Tstar=%f (cp0 data: cp0star=%f, Tstar=%f)",R,Tstar,I->data.cp0.cp0star, I->data.cp0.Tstar);
76 }else{
77 MSG("R=%f,Tstar=%f (cp0 data: Tstar=%f)",R,Tstar,I->data.phi0.Tstar);
78 }
79 #endif
80
81 Phi0RunData *N = FPROPS_NEW(Phi0RunData);
82 int i, add_const_term=1;
83 double Tred, cp0red, p;
84 N->c = 0;
85 N->m = 0;
86 switch(I->type){
87 case IDEAL_CP0:
88 N->np = I->data.cp0.np;
89 N->ne = I->data.cp0.ne;
90 Tred = I->data.cp0.Tstar;
91 cp0red = I->data.cp0.cp0star;
92 MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);
93 MSG("Tred = %f, Tstar = %f", Tred, Tstar);
94 MSG("cp0red = %f, R = %f", cp0red, R);
95
96 for(i=0; i < N->np; ++i)if(I->data.cp0.pt[i].t == 0)add_const_term = 0;
97 //MSG("add_const_term = %d",add_const_term);
98
99 N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np + add_const_term);
100 N->et = FPROPS_NEW_ARRAY(Phi0RunExpTerm, N->ne);
101 add_const_term = 1;
102
103 for(i=0; i < N->np; ++i){
104 p = - I->data.cp0.pt[i].t;
105 N->pt[i].p = p;
106 if(p == 0){
107 //if(add_const_term)fprintf(stderr,"Addding offset here\n");
108 N->pt[i].a = (I->data.cp0.pt[i].c - add_const_term);
109 add_const_term = 0;
110 }else{
111 N->pt[i].a = -I->data.cp0.pt[i].c / p / (p - 1) * pow(Tred / Tstar, p);
112 }
113 }
114
115 if(add_const_term){
116 //MSG("WARNING: adding constant term %d in cp0, is that what you really want?",i);
117 N->pt[i].a = -1;
118 N->pt[i].p = 0;
119 N->np++;
120 }
121
122 for(i=0; i < N->ne; ++i){
123 N->et[i].n = I->data.cp0.et[i].b;
124 N->et[i].gamma = I->data.cp0.et[i].beta * Tred / Tstar;
125 }
126
127 if(cp0red != R){
128 //MSG("WARNING: adjusting for R (=%f) != cp0red (=%f)...\n",R,cp0red);
129 double X = cp0red / R;
130 // scale for any differences in R and cpstar */
131 for(i=0; i < N->np; ++i)N->pt[i].a *= X;
132 for(i=0; i < N->ne; ++i)N->et[i].n *= X;
133 }
134
135 break;
136 case IDEAL_PHI0:
137 // TODO add checks for disallowed terms, eg p = 0 or p = 1?
138 N->np = I->data.phi0.np;
139 N->ne = I->data.phi0.ne;
140 //MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);
141
142 // power terms
143 N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np);
144 for(i=0; i< N->np; ++i){
145 double a = I->data.phi0.pt[i].a0;
146 double p = I->data.phi0.pt[i].p0;
147 N->pt[i].a = a;
148 N->pt[i].p = p;
149 }
150
151 // exponential terms
152 N->et = FPROPS_NEW_ARRAY(Phi0RunExpTerm,N->ne);
153 for(i=0; i < N->ne; ++i){
154 double n = I->data.phi0.et[i].n;
155 double gamma = I->data.phi0.et[i].gamma;
156 N->et[i].gamma = gamma;
157 N->et[i].n = n;
158 }
159 break;
160 }
161 return N;
162 }
163
164 void cp0_destroy(Phi0RunData *N){
165 if(N->pt)FPROPS_FREE(N->pt);
166 if(N->et)FPROPS_FREE(N->et);
167 FPROPS_FREE(N);
168 }
169
170 /*---------------------------------------------
171 IDEAL COMPONENT RELATIONS
172 */
173
174 /*
175 in calculating the ideal component relations, we have some integration
176 constants that ultimately allow the arbitrary scales of h and s to
177 be specified. Hence we can ignore components of helm_ideal that
178 are either constant or linear in tau, and work them out later when we
179 stick in the values of data->c and data->m.
180
181 See refstate.h, ReferenceState for the various ways we have for setting
182 c and m using reference state specifications.
183 */
184
185 //#define IDEAL_DEBUG
186 /**
187 Ideal component of helmholtz function
188 */
189 double ideal_phi(double tau, double delta, const Phi0RunData *data){
190
191 const Phi0RunPowTerm *pt;
192 const Phi0RunExpTerm *et;
193
194 unsigned i;
195 // FIXME what if rhostar != rhoc??
196 double sum = log(delta) + data->c + data->m * tau;
197 double term;
198
199 #ifdef IDEAL_DEBUG
200 fprintf(stderr,"\ttau = %f, delta = %f\n",tau,delta);
201 fprintf(stderr,"\tT ~ %f\n\n",Tstar_on_tau);
202 fprintf(stderr,"\tlog(delta) + c + m tau = %f (c=%f,m=%f)\n",sum,data->c, data->m);
203 fprintf(stderr,"sum = %f\n",sum);
204 #endif
205
206 /* power terms */
207 pt = &(data->pt[0]);
208 for(i = 0; i<data->np; ++i, ++pt){
209 double a = pt->a;
210 double p = pt->p;
211 if(p == 0){
212 term = a*log(tau);
213 #ifdef IDEAL_DEBUG
214 fprintf(stderr,"\ta log(tau) = %f (a=%f, t=%f)\n",term,a,p);
215 #endif
216 }else{
217 #ifdef IDEAL_DEBUG
218 assert(p!=-1);
219 #endif
220 term = a * pow(tau, p);
221 #ifdef IDEAL_DEBUG
222 fprintf(stderr,"\ta tau^p = %f (a=%f, p=%f)\n",term,a,p);
223 #endif
224 }
225 sum += term;
226 #ifdef IDEAL_DEBUG
227 fprintf(stderr,"sum = %f\n",sum);
228 #endif
229 }
230
231 /* Planck-Einstein terms */
232 et = &(data->et[0]);
233 for(i=0; i<data->ne; ++i, ++et){
234 term = et->n * log(1 - exp(-et->gamma * tau));
235 #ifdef IDEAL_DEBUG
236 fprintf(stderr,"\tn log(1-exp(-gamma*tau)) = %f, (n=%f, gamma=%f)\n",term,et->n,et->gamma);
237 #endif
238 sum += term;
239 }
240
241 #ifdef IDEAL_DEBUG
242 fprintf(stderr,"phi0 = %f\n",sum);
243 #endif
244
245 return sum;
246 }
247
248 /**
249 Partial dervivative of ideal component (phi0) of normalised helmholtz
250 residual function (phi), with respect to tau.
251 */
252 double ideal_phi_tau(double tau, double delta, const Phi0RunData *data){
253 const Phi0RunPowTerm *pt;
254 const Phi0RunExpTerm *et;
255
256 unsigned i;
257 double term;
258 double sum = data->m;
259
260 pt = &(data->pt[0]);
261 for(i = 0; i<data->np; ++i, ++pt){
262 double a = pt->a;
263 double p = pt->p;
264 if(p == 0){
265 term = a / tau;
266 //fprintf(stderr,"\tc/tau = %f\n",term);
267 }else{
268 term = a*p*pow(tau,p - 1);
269 //fprintf(stderr,"\tc / tau^p = %f (t=%f, c=%.3e, p=%f)\n",term,t,coeff,-t-1.);
270 }
271 #ifdef TEST
272 if(isinf(term)){
273 fprintf(stderr,"Error with infinite-valued term with i = %d, a = %f, p = %f\n", i,a ,p);
274 abort();
275 }
276 #endif
277 assert(!isnan(term));
278 sum += term;
279 }
280
281 /* Planck-Einstein terms */
282 et = &(data->et[0]);
283 for(i=0; i<data->ne; ++i, ++et){
284 double expo = exp(-et->gamma * tau);
285 term = et->n * et->gamma * expo / (1 - expo);
286 #ifdef TEST
287 assert(!isinf(term));
288 #endif
289 sum += term;
290 }
291
292 #ifdef TEST
293 assert(!isinf(sum));
294 #endif
295 return sum;
296 }
297
298
299
300 /**
301 Second partial dervivative of ideal component (phi0) of normalised helmholtz
302 residual function (phi), with respect to tau. This one is easy!
303 It's not a function of delta.
304
305 FIXME: although this one is easy, we want to pull Tstar out of the
306 ideal properties stuff, if that's possible.
307 */
308 double ideal_phi_tautau(double tau, const Phi0RunData *data){
309 const Phi0RunPowTerm *pt;
310 const Phi0RunExpTerm *et;
311
312 unsigned i;
313
314 double sum = 0;
315 double term;
316
317 #ifdef IDEAL_DEBUG
318 fprintf(stderr,"\ttau = %f\n",tau);
319 #endif
320
321 /* power terms */
322 pt = &(data->pt[0]);
323 for(i = 0; i<data->np; ++i, ++pt){
324 if(pt->p == 0){
325 term = pt->a;
326 }else{
327 term = -pt->a * pt->p * (pt->p - 1) * pow(tau, pt->p);
328 }
329 #ifdef IDEAL_DEBUG
330 fprintf(stderr,"\tpt[%d] = ap(p-1)*tau^p (a = %e, p = %e) = %f\n",i,pt->a, pt->p, term);
331 #endif
332 sum += term;
333 }
334
335 /* 'exponential' terms */
336 et = &(data->et[0]);
337 for(i=0; i<data->ne; ++i, ++et){
338 double x = et->gamma * tau;
339 double e = exp(-x);
340 double d = (1-e)*(1-e);
341 term = et->n * x*x * e / d;
342 #ifdef IDEAL_DEBUG
343 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);
344 #endif
345 sum += term;
346 }
347 /* note, at this point, sum == cp0/R - 1 */
348 #ifdef IDEAL_DEBUG
349 MSG("sum = %f, phi_tautau = %f",sum, -sum/SQ(tau));
350 #endif
351 return -sum/SQ(tau);
352 }

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