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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1825 - (show annotations) (download) (as text)
Thu Aug 21 05:03:19 2008 UTC (15 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 6161 byte(s)
Added functions for enthalpy and internal energy.
Test suites still needed!
1 /* ASCEND modelling environment
2 Copyright (C) 2006 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 *//** @file
19 Implementation of the reduced molar Helmholtz free energy equation of state.
20
21 For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
22 Fundamentalgleichung f端r Ammoniak.
23
24 John Pye, 29 Jul 2008.
25 */
26
27 #include <math.h>
28
29 #include "helmholtz.h"
30
31 /* forward decls */
32
33 static double helm_ideal(double tau, double delta, HelmholtzData *data);
34 static double helm_ideal_tau(double tau, double delta, HelmholtzData *data);
35 static double helm_resid(double tau, double delta, HelmholtzData *data);
36 static double helm_resid_del(double tau, double delta, HelmholtzData *data);
37 static double helm_resid_tau(double tau, double delta, HelmholtzData *data);
38
39 /**
40 Function to calculate pressure from Helmholtz free energy EOS, given temperature
41 and mass density.
42
43 @param T temperature in K
44 @param rho mass density in kg/m続
45 @return pressure in Pa???
46 */
47 double helmholtz_p(double T, double rho, HelmholtzData *data){
48
49 double tau = data->T_star / T;
50 double delta = rho / data->rho_star;
51
52 return data->R * T * rho * (1. + delta * helm_resid_del(tau,delta,data));
53 }
54
55 /**
56 Function to calculate internal energy from Helmholtz free energy EOS, given
57 temperature and mass density.
58
59 @param T temperature in K
60 @param rho mass density in kg/m続
61 @return internal energy in ???
62 */
63 double helmholtz_u(double T, double rho, HelmholtzData *data){
64
65 double tau = data->T_star / T;
66 double delta = rho / data->rho_star;
67
68 return data->R * T * tau * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data));
69 }
70
71 /**
72 Function to calculate enthalpy from Helmholtz free energy EOS, given
73 temperature and mass density.
74
75 @param T temperature in K
76 @param rho mass density in kg/m続
77 @return enthalpy in ????
78 */
79 double helmholtz_h(double T, double rho, HelmholtzData *data){
80
81 double tau = data->T_star / T;
82 double delta = rho / data->rho_star;
83
84 return data->R * T * (1 + tau * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data)) + delta*helm_resid_del(tau,delta,data));
85 }
86
87 /*---------------------------------------------
88 IDEAL COMPONENT RELATIONS
89 */
90
91 /**
92 Ideal component of helmholtz function
93 */
94 double helm_ideal(double tau, double delta, HelmholtzData *data){
95 double *a0;
96
97 double tau13 = pow(tau,1./3.);
98 double taum32 = pow(tau,-3./2.);
99 double taum74 = pow(tau,-7./4.);
100
101 a0 = &(data->a0[0]);
102 return log(delta) + a0[0] + a0[1] * tau - log(tau) + a0[2] * tau13 + a0[3] * taum32 + a0[4] * taum74;
103 }
104
105 /**
106 Partial dervivative of ideal component of helmholtz residual function with
107 respect to tau.
108 */
109 double helm_ideal_tau(double tau, double delta, HelmholtzData *data){
110 double *a0;
111
112 double taum114 = pow(tau,-11./4.);
113 double taum74 = tau * taum114;
114 double taum23 = pow(tau,-2./3.);
115
116 a0 = &(data->a0[0]);
117
118 return a0[1] - 1./tau + a0[2]/3.* taum23 - a0[3] * 3./4. * taum74 - a0[4] * 7./4. * taum114;
119 }
120
121 /**
122 Residual part of helmholtz function. Note: we have NOT prematurely
123 optimised here ;-)
124 */
125 double helm_resid(double tau, double delta, HelmholtzData *data){
126
127 double sum;
128 double phir = 0;
129 unsigned i;
130
131 HelmholtzATD *atd = &(data->atd[0]);
132
133 for(i=0; i<5; ++i){
134 phir += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
135 ++atd;
136 }
137
138 sum = 0;
139 for(i=5; i<10; ++i){
140 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
141 ++atd;
142 }
143 phir += exp(-delta) * sum;
144
145 sum = 0;
146 for(i=10; i<17; ++i){
147 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
148 ++atd;
149 }
150 phir += exp(-delta*delta) * sum;
151
152 sum = 0;
153 for(i=17; i<21; ++i){
154 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
155 ++atd;
156 }
157 phir += exp(-delta*delta*delta) * sum;
158
159 return phir;
160 }
161
162 /**
163 Derivative of the helmholtz residual function with respect to
164 delta.
165 */
166 double helm_resid_del(double tau,double delta,HelmholtzData *data){
167
168 double sum;
169 double phir = 0;
170 unsigned i;
171 double XdelX;
172
173 HelmholtzATD *atd = &(data->atd[0]);
174
175 for(i=0; i<5; ++i){
176 phir += atd->a * pow(tau, atd->t) * pow(delta, atd->d - 1) * atd->d;
177 ++atd;
178 }
179
180 sum = 0;
181 XdelX = delta;
182 for(i=5; i<10; ++i){
183 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d) * (atd->d - XdelX);
184 }
185 phir += exp(-delta) * sum;
186
187 sum = 0;
188 XdelX = 2*delta*delta;
189 for(i=10; i<17; ++i){
190 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d) * (atd->d - XdelX);
191 }
192 phir += exp(-delta*delta) * sum;
193
194 sum = 0;
195 XdelX = 3*delta*delta*delta;
196 for(i=17; i<21; ++i){
197 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d) * (atd->d - XdelX);
198 }
199 phir += exp(-delta*delta*delta) * sum;
200
201 return phir;
202 }
203
204 /**
205 Derivative of the helmholtz residual function with respect to
206 tau.
207 */
208 double helm_resid_tau(double tau,double delta,HelmholtzData *data){
209
210 double sum;
211 double phir = 0;
212 unsigned i;
213
214 HelmholtzATD *atd = &(data->atd[0]);
215
216 for(i=0; i<5; ++i){
217 if(atd->t != 0){
218 phir += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
219 }
220 ++atd;
221 }
222
223 sum = 0;
224 for(i=5; i<10; ++i){
225 if(atd->t != 0){
226 sum += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
227 }
228 ++atd;
229 }
230 phir += exp(-delta) * sum;
231
232 sum = 0;
233 for(i=10; i<17; ++i){
234 if(atd->t != 0){
235 sum += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
236 }
237 ++atd;
238 }
239 phir += exp(-delta*delta) * sum;
240
241 sum = 0;
242 for(i=17; i<21; ++i){
243 if(atd->t != 0){
244 sum += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
245 }
246 ++atd;
247 }
248 phir += exp(-delta*delta*delta) * sum;
249
250 return phir;
251 }
252
253
254
255

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