/[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 1857 - (show annotations) (download) (as text)
Thu Sep 11 14:46:51 2008 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 4073 byte(s)
p, u, h working for nitrogen, still problems with s.
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->a0 * pow(T, pt->t0);
54 #if 0
55 fprintf(stderr,"i = %d: ",i);
56 fprintf(stderr,"power term, a = %f, t = %f, val = %f\n",pt->a0, pt->t0, 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 a = pt->a0;
109 double t = pt->t0;
110 if(t == 0){
111 term = a*log(tau);
112 }else{
113 #ifdef TEST
114 assert(t!=-1);
115 #endif
116 term = -a / (t*(t+1)) * pow(Tstar_on_tau,t);
117 }
118 sum += pt->a0 * pow(tau, pt->t0);
119 }
120
121 /* 'exponential' terms */
122 et = &(data->et[0]);
123 for(i=0; i<data->ne; ++i, ++et){
124 sum += et->b * log(1 - exp(-et->beta / Tstar_on_tau));
125 }
126
127 return sum;
128 }
129
130 /**
131 Partial dervivative of ideal component of helmholtz residual function with
132 respect to tau.
133 */
134 double helm_ideal_tau(double tau, double delta, const IdealData *data){
135 const IdealPowTerm *pt;
136 const IdealExpTerm *et;
137
138 unsigned i;
139 double term;
140 double sum = -1./tau + data->m;
141 double Tstar_on_tau = data->Tstar / tau;
142
143 pt = &(data->pt[0]);
144 for(i = 0; i<data->np; ++i, ++pt){
145 double a = pt->a0;
146 double t = pt->t0;
147 if(t==0){
148 term = a / tau;
149 }else{
150 // term = -a / (t*(t+1)) * pow(Tstar_on_tau,t);
151 term = a*pow(Tstar_on_tau,t)*tau/(t+1);
152 }
153 #ifdef TEST
154 if(isinf(term)){
155 fprintf(stderr,"Error with infinite-valued term with i = %d, a = %f, t = %f\n", i,a ,t);
156 abort();
157 }
158 #endif
159 sum += term;
160 }
161
162 /* 'exponential' terms */
163 et = &(data->et[0]);
164 for(i=0; i<data->ne; ++i, ++et){
165 term = et->b * et->beta / data->Tstar / (1 - exp(-et->beta/Tstar_on_tau));
166 #ifdef TEST
167 assert(!isinf(term));
168 #endif
169 sum += term;
170 }
171
172 #ifdef TEST
173 assert(!isinf(sum));
174 #endif
175 return sum;
176 }

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