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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2266 - (show annotations) (download) (as text)
Thu Aug 5 13:53:01 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 8013 byte(s)
The (p,h) function is converging almost everywhere. The home-made Newton iteration wasn't starting well for supercritical states
when the starting guess had rho = rho_c, so changing slightly away from fixed almost all the problems.
1 /*
2 ASCEND modelling environment
3 Copyright (C) 2004-2010 John Pye
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License
7 as published by the Free Software Foundation; either version 2
8 of the License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18 *//** @file
19 The following are the functions
20
21 ⎰ ∂z ⎱ ___\ VTn
22 ⎱ ∂v ⎰T /
23
24 etc., for saturation and non-saturation regions.
25 */
26
27 #include "derivs.h"
28 #include "helmholtz.h"
29 #include "sat.h"
30
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <math.h>
34
35 //#define FPE_DEBUG
36 #ifdef FPE_DEBUG
37 # include <assert.h>
38 #else
39 # define assert(ARGS...)
40 #endif
41
42
43 #define SQ(X) ((X)*(X))
44
45 /*------------------------------------------------------------------------------
46 EXPORTED FUNCTION(S)
47 */
48
49 /**
50 Calculates the derivative
51
52 ⎰ ∂z ⎱
53 ⎱ ∂x ⎰y
54
55 @param S steam state, already calculated using steam_ph, etc.
56 @param x in above equation, character one of 'pTvuhsgaf'.
57 @param y in above equation, character one of pTvuhsgaf.
58 @param z in above equation, character one of pTvuhsgaf.
59 Note that Helmholtz free energy can be signified by either 'a' or 'f'.
60
61 @NOTE that the variables ARE NOT IN ALPHABETICAL ORDER.
62
63 @return the numerical value of the derivative (∂z/∂x)y.
64 */
65 double fprops_deriv(FPROPS_CHAR z, FPROPS_CHAR x, FPROPS_CHAR y, double T, double rho, const HelmholtzData *D){
66 StateData S;
67 S.rho = rho; S.T = T;
68 char sat = 0;
69
70 if(T < D->T_c){
71 int res = fprops_sat_T(T,&S.psat, &S.rhof, &S.rhog, D);\
72 if(res){
73 fprintf(stderr,"Failed to calculation saturation props\n");
74 exit(1);
75 }
76 if(rho < S.rhof && rho > S.rhog){
77 /* we're in the saturation region */
78 sat = 1;
79 }
80 }
81
82 double ZTV, ZVT, XTV, XVT, YTV, YVT;
83 if(sat){
84 /* saturated. first use clapeyron equation for (∂p/∂T)v */
85 double hf, hg;
86 hf = helmholtz_h(S.T, S.rhof,D);
87 hg = helmholtz_h(S.T, S.rhog,D);
88 S.dpdT_sat = (hg - hf) / T / (1./S.rhog - 1./S.rhof);
89
90 fprintf(stderr,"Saturation region derivatives not yet implemented.\n");
91 exit(1);
92
93 /* ...then call the partial deriv routines, which use dpdT_sat */
94 #define TVSAT(X) fprops_sat_dZdT_v(X,&S)
95 #define VTSAT(X) fprops_sat_dZdv_T(X,&S)
96 ZTV = TVSAT(z); ZVT = VTSAT(z);
97 XTV = TVSAT(x); XVT = VTSAT(x);
98 YTV = TVSAT(y); YVT = VTSAT(y);
99 }else{
100 /* non-saturated */
101 #define TVNON(X) fprops_non_dZdT_v(X,T,rho,D)
102 #define VTNON(X) fprops_non_dZdv_T(X,T,rho,D)
103 ZTV = TVNON(z); ZVT = VTNON(z);
104 XTV = TVNON(x); XVT = VTNON(x);
105 YTV = TVNON(y); YVT = VTNON(y);
106 }
107 #undef TVNON
108 #undef VTNON
109
110 double deriv = ((ZTV*YVT-ZVT*YTV)/(XTV*YVT-XVT*YTV));
111 //fprintf(stderr,"Calculated (∂%c/∂%c)%c = %g\n",z,x,y,deriv);
112 return deriv;
113 }
114
115 /*------------------------------------------------------------------------------
116 Non-saturation derivatives... easy.
117 */
118
119 /*
120 FIXME the following macros avoid calculating unneeded results eg within VT3
121 but at the level of freesteam_deriv, there is wasted effort, because eg 'p'
122 will be calculated several times in different calls to VT3.
123 */
124
125 #define p helmholtz_p_raw(T,rho,D)
126 #define cv helmholtz_cv(T,rho,D)
127 #define v (1./rho)
128 #define s helmholtz_s_raw(T,rho,D)
129 #define alphap helmholtz_alphap(T,rho,D)
130 #define betap helmholtz_betap(T,rho,D)
131
132 double fprops_non_dZdv_T(FPROPS_CHAR x, double T, double rho, const HelmholtzData *D){
133 double res;
134 switch(x){
135 case 'p': res = -p*betap; break;
136 case 'T': res = 0; break;
137 case 'v': res = 1; break;
138 case 'u': res = p*(T*alphap-1.); break;
139 case 'h': res = p*(T*alphap-v*betap); break;
140 case 's': res = p*alphap; break;
141 case 'g': res = -p*v*betap; break;
142 case 'a':
143 case 'f': res = -p; break;
144 default:
145 fprintf(stderr,"%s (%s:%d): Invalid variable '%c'\n", __func__,__FILE__,__LINE__,x);
146 exit(1);
147 }
148 #if 1
149 if(__isnan(res)){
150 fprintf(stderr,"calculating '%c'\n",x);
151 }
152 #endif
153 assert(!__isnan(res));
154 //fprintf(stderr,"(∂%c/∂v)T = %f\n",x,res);
155 return res;
156 }
157
158 double fprops_non_dZdT_v(FPROPS_CHAR x, double T, double rho, const HelmholtzData *D){
159 double res;
160 switch(x){
161 case 'p': res = p*alphap; break;
162 case 'T': res = 1; break;
163 case 'v': res = 0; break;
164 case 'u': res = cv; break;
165 case 'h': res = cv + p*v*alphap; break;
166 case 's': res = cv/T; break;
167 case 'g': res = p*v*alphap - s; break;
168 case 'a':
169 case 'f': res = -s; break;
170 default:
171 fprintf(stderr,"%s (%s:%d): Invalid variable '%c'\n", __func__,__FILE__,__LINE__,x);
172 exit(1);
173 }
174 #if 0
175 if(__isnan(res)){
176 fprintf(stderr,"calculating '%c'\n",x);
177 }
178 #endif
179 assert(!__isnan(res));
180 //fprintf(stderr,"(∂%c/∂T)v = %f\n",x,res);
181 return res;
182 }
183 #undef p
184 #undef cv
185 #undef v
186 #undef s
187 #undef alphap
188 #undef betap
189
190
191 /*------------------------------------------------------------------------------
192 Saturation region derivatives... a bit harder.
193
194 FPROPS uses T,rho as coordinates, so we need to work out what these
195 derivatives are going to be within the saturation region.
196 */
197
198 /*
199 ⎰ ∂z ⎱ = ⎰∂z_f⎱ (1 - x) + ⎰∂z_g⎱ x
200 ⎱ ∂T ⎰v ⎱ ∂T ⎰ ⎱ ∂T ⎰
201
202 Need to review the theory of this section, might not be correct.
203 */
204 double fprops_sat_dZdT_v(FPROPS_CHAR z, const StateData *S){
205 double res;
206 switch(z){
207 /* a couple of easy cases, first: */
208 case 'p': return S->dpdT_sat;
209 case 'T': return 1.;
210 }
211
212 double drhofdT = fprops_drhofdT(S);
213 double drhogdT = fprops_drhogdT(S);
214
215 double dzfdT, dzgdT;
216
217 assert(S->rhof!=0);
218 assert(S->rhog!=0);
219 double dvfdT = -1./SQ(S->rhof) * drhofdT;
220 assert(!isnan(dvfdT));
221 double dvgdT = -1./SQ(S->rhog) * drhogdT;
222 assert(!isnan(dvgdT));
223
224 #define TVNON(X,RHO) fprops_non_dZdT_v(X,S->T,RHO,S->D)
225 #define VTNON(X,RHO) fprops_non_dZdv_T(X,S->T,RHO,S->D)
226 dzfdT = VTNON(z,S->rhof)*dvfdT + TVNON(z,S->rhof);
227 dzgdT = VTNON(z,S->rhog)*dvgdT + TVNON(z,S->rhog);
228
229 assert(!isnan(dzfdT));
230 assert(!isnan(dzgdT));
231 double x;
232
233 /* FIXME this is not the correct solution, this is for x const, not rho const. */
234
235 x = 1./(x/S->rhog + (1-x)/S->rhof);
236 res = dzfdT*(1-x) + dzgdT*x;
237 //fprintf(stderr,"(∂%c/∂T)x = %g\n",z,res);
238 return res;
239 }
240
241 /*
242 These derivatives are simply the gradient within the two-phase region,
243 and is very simply calculated as
244
245 ⎰ ∂z ⎱ = (z_g - z_f) / (rho_g - rho_f) where z in {v,u,h,z,s,g,a}.
246 ⎱ ∂v ⎰T
247
248 or, otherwise,
249
250 ⎰ ∂T ⎱ , ⎰ ∂p ⎱ = 0
251 ⎱ ∂x ⎰T ⎱ ∂x ⎰T
252
253
254 Need to double-check theory in this section.
255 */
256 double fprops_sat_dZdv_T(FPROPS_CHAR z, const StateData *S){
257 switch(z){
258 case 'p': return 0;
259 case 'T': return 0;
260 }
261 double zf, zg;
262 #define ZFG(Z,P,T) \
263 zf = helmholtz_##Z(S->T,S->rhof,S->D);\
264 zg = helmholtz_##Z(S->T,S->rhog,S->D)
265 switch(z){
266 case 'v': zf = 1./S->rhof; zg = 1./S->rhog; break;
267 case 'u': ZFG(u,p,T); break;
268 case 'h': ZFG(h,p,T); break;
269 case 's': ZFG(s,p,T); break;
270 case 'g': ZFG(g,p,T); break;
271 case 'a': case 'f': ZFG(a,p,T); break;
272 default:
273 fprintf(stderr,"%s (%s:%d): Invalid character x = '%c'\n", __func__,__FILE__,__LINE__,z);
274 exit(1);
275 }
276 #undef ZFG
277 //fprintf(stderr,"(∂%c/∂x)T = %g\n",z,zg-zf);
278 return zg - zf;
279 }
280
281 /*------------------------------------------------------------------------------
282 DERIVATIVES OF rhof and rhog with temperature
283 */
284
285 double fprops_drhofdT(const StateData *S){
286 double dpdT = TVNON('p',S->rhof);
287 double dpdrho = -1./SQ(S->rhof) * VTNON('p',S->rhof);
288 return (S->dpdT_sat - dpdT)/dpdrho;
289 }
290
291 double fprops_drhogdT(const StateData *S){
292 double dpdT = TVNON('p',S->rhog);
293 double dpdrho = -1./SQ(S->rhog) * VTNON('p',S->rhog);
294 return (S->dpdT_sat - dpdT)/dpdrho;
295 }
296
297
298

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