/[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 1865 - (show annotations) (download) (as text)
Mon Sep 15 08:40:14 2008 UTC (11 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 4266 byte(s)
Still working on fixing helmholtz_a.
1 /* 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 */
19
20 #include <math.h>
21
22 #include "ideal.h"
23
24 #ifdef TEST
25 #include <assert.h>
26 #include <stdlib.h>
27 #include <stdio.h>
28 #endif
29
30 /*---------------------------------------------
31 IDEAL COMPONENT RELATIONS
32 */
33
34 /**
35 Zero-pressure specific heat function (ideal gas limit)
36
37 This is returned in the units of data->cp0star.
38 */
39 double helm_cp0(double T, const IdealData *data){
40 const IdealPowTerm *pt;
41 const IdealExpTerm *et;
42
43 unsigned i;
44 double sum = 0;
45 double term;
46
47 /* power terms */
48 pt = &(data->pt[0]);
49 #if 0
50 fprintf(stderr,"np = %d\n",data->np);
51 #endif
52 for(i = 0; i<data->np; ++i, ++pt){
53 term = pt->c * pow(T, pt->t);
54 #if 0
55 fprintf(stderr,"i = %d: ",i);
56 fprintf(stderr,"power term, c = %f, t = %f, val = %f\n",pt->c, pt->t, term);
57 #endif
58 sum += term;
59 }
60
61 /* 'exponential' terms */
62 et = &(data->et[0]);
63 for(i=0; i<data->ne; ++i, ++et){
64 #if 0
65 fprintf(stderr,"exp term\n");
66 #endif
67 double x = et->beta / T;
68 double e = exp(-x);
69 double d = (1-e)*(1-e);
70 term = et->b * x*x * e / d;
71 #if 0
72 fprintf(stderr,"exp term, b = %f, beta = %f, val = %f\n",et->b, et->beta, term);
73 #endif
74 sum += term;
75 }
76
77 #if 0
78 fprintf(stderr,"Mult by cp0* = %f\n",data->cp0star);
79 #endif
80 return data->cp0star * sum;
81 }
82
83 /*
84 Hypothesis:
85 in calculating the ideal component relations, we have some integration
86 constants that ultimately allow the arbitrary scales of h and s to
87 be specified. Hence we can ignore components of helm_ideal that
88 are either constant or linear in tau, and work them out later when we
89 stick in the values of data->c and data->m.
90 */
91
92 /**
93 Ideal component of helmholtz function
94 */
95 double helm_ideal(double tau, double delta, const IdealData *data){
96
97 const IdealPowTerm *pt;
98 const IdealExpTerm *et;
99
100 unsigned i;
101 double sum = log(delta) - log(tau) + data->c + data->m * tau;
102 double term;
103 double Tstar_on_tau = data->Tstar / tau;
104
105 /* power terms */
106 pt = &(data->pt[0]);
107 for(i = 0; i<data->np; ++i, ++pt){
108 double c = pt->c;
109 double t = pt->t;
110 if(t == 0){
111 term = c*log(tau);
112 }else{
113 #ifdef TEST
114 assert(t!=-1);
115 #endif
116 term = -c / (t*(t+1)) * pow(Tstar_on_tau,t);
117 fprintf(stderr,"i = %d, c = %f, t = %f, term = %f\n",i,c,t,term);
118 }
119 sum += term;
120 }
121
122 /* 'exponential' terms */
123 et = &(data->et[0]);
124 for(i=0; i<data->ne; ++i, ++et){
125 term = et->b * log(1 - exp(-et->beta / Tstar_on_tau));
126 fprintf(stderr,"exp i=%d, b=%f, beta=%f, term = %f\n",i,et->b, et->beta, term);
127 sum += term;
128 }
129
130 #ifdef TEST
131 fprintf(stderr,"phi0 = %f\n",sum);
132 #endif
133
134 return sum;
135 }
136
137 /**
138 Partial dervivative of ideal component of helmholtz residual function with
139 respect to tau.
140 */
141 double helm_ideal_tau(double tau, double delta, const IdealData *data){
142 const IdealPowTerm *pt;
143 const IdealExpTerm *et;
144
145 unsigned i;
146 double term;
147 double sum = -1./tau + data->m;
148 double Tstar_on_tau = data->Tstar / tau;
149
150 pt = &(data->pt[0]);
151 for(i = 0; i<data->np; ++i, ++pt){
152 double c = pt->c;
153 double t = pt->t;
154 if(t==0){
155 term = c / tau;
156 }else{
157 // term = -c / (t*(t+1)) * pow(Tstar_on_tau,t);
158 term = c/(t+1)*pow(Tstar_on_tau,t)/tau;
159 }
160 #ifdef TEST
161 if(isinf(term)){
162 fprintf(stderr,"Error with infinite-valued term with i = %d, c = %f, t = %f\n", i,c ,t);
163 abort();
164 }
165 #endif
166 sum += term;
167 }
168
169 /* 'exponential' terms */
170 et = &(data->et[0]);
171 for(i=0; i<data->ne; ++i, ++et){
172 term = et->b * et->beta / data->Tstar / (1 - exp(-et->beta/Tstar_on_tau));
173 #ifdef TEST
174 assert(!isinf(term));
175 #endif
176 sum += term;
177 }
178
179 #ifdef TEST
180 assert(!isinf(sum));
181 #endif
182 return sum;
183 }

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