/[ascend]/branches/sid/models/johnpye/fprops/ttse.c
ViewVC logotype

Contents of /branches/sid/models/johnpye/fprops/ttse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3025 - (show annotations) (download) (as text)
Thu Jul 23 18:43:12 2015 UTC (3 years, 2 months ago) by sid
File MIME type: text/x-csrc
File size: 13671 byte(s)


1
2 #include <math.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5
6
7 #include "rundata.h"
8 #include "ttse.h"
9
10
11
12 #define TTSE_DEBUG //sid change
13 #ifdef TTSE_DEBUG
14 # include "color.h"
15 # define MSG FPROPS_MSG
16 # define ERRMSG FPROPS_ERRMSG
17 #else
18 # define MSG(ARGS...) ((void)0)
19 # define ERRMSG(ARGS...) ((void)0)
20 #endif
21
22
23 #include <ascend/general/ascMalloc.h>
24 #include <ascend/utilities/error.h>
25
26 #define INIT -23456789.123
27
28 //#define FO
29 #define SO
30 inline TtseMatrix alloc_matrix(int tp, int rhop) {
31
32 TtseMatrix matrix = ASC_NEW_ARRAY (TtseMatrix, tp*rhop); //tp rows and rhop columns
33 /*
34 matrix[10][24] = 99;
35 printf("%f\n", *((double*)matrix+(tp*10)+24));
36 matrix[20][44] = -99;
37 printf("%f\n", *((double*)matrix+(tp*20)+44));
38 matrix[30][64] = 199;
39 printf("%f\n", *((double*)matrix+(tp*30)+64));
40 */
41 return matrix;
42 }
43
44
45 inline void remove_matrix(TtseMatrix mat , int tp){
46
47 ASC_FREE(mat);
48 }
49
50
51 void alloc_tables(Ttse * table)
52 {
53
54 table->satFRho = ASC_NEW_ARRAY ( double , NSAT);
55 table->satFdRhodt = ASC_NEW_ARRAY (double, NSAT);
56 table->satFd2RhodT2 = ASC_NEW_ARRAY (double, NSAT);
57 table->satGRho = ASC_NEW_ARRAY (double, NSAT);
58 table->satGdRhodt = ASC_NEW_ARRAY (double, NSAT);
59 table->satGd2RhodT2 = ASC_NEW_ARRAY (double, NSAT);
60
61
62
63 table->s = alloc_matrix(NTP,NRHOP);
64 table->dsdt = alloc_matrix(NTP,NRHOP);
65 table->d2sdt2 = alloc_matrix(NTP,NRHOP);
66 table->dsdrho = alloc_matrix(NTP,NRHOP);
67 table->d2sdrho2 = alloc_matrix(NTP,NRHOP);
68 table->d2sdtdrho = alloc_matrix(NTP,NRHOP);
69
70
71 table->p = alloc_matrix(NTP,NRHOP);
72 table->dpdt = alloc_matrix(NTP,NRHOP);
73 table->d2pdt2 = alloc_matrix(NTP,NRHOP);
74 table->dpdrho = alloc_matrix(NTP,NRHOP);
75 table->d2pdrho2 = alloc_matrix(NTP,NRHOP);
76 table->d2pdtdrho = alloc_matrix(NTP,NRHOP);
77
78
79 table->u = alloc_matrix(NTP,NRHOP);
80 table->dudt = alloc_matrix(NTP,NRHOP);
81 table->d2udt2 = alloc_matrix(NTP,NRHOP);
82 table->dudrho = alloc_matrix(NTP,NRHOP);
83 table->d2udrho2 = alloc_matrix(NTP,NRHOP);
84 table->d2udtdrho = alloc_matrix(NTP,NRHOP);
85
86
87 table->g = alloc_matrix(NTP,NRHOP);
88 table->dgdt = alloc_matrix(NTP,NRHOP);
89 table->d2gdt2 = alloc_matrix(NTP,NRHOP);
90 table->dgdrho = alloc_matrix(NTP,NRHOP);
91 table->d2gdrho2 = alloc_matrix(NTP,NRHOP);
92 table->d2gdtdrho = alloc_matrix(NTP,NRHOP);
93
94
95
96 table->h = alloc_matrix(NTP,NRHOP);
97 table->dhdt = alloc_matrix(NTP,NRHOP);
98 table->d2hdt2 = alloc_matrix(NTP,NRHOP);
99 table->dhdrho = alloc_matrix(NTP,NRHOP);
100 table->d2hdrho2 = alloc_matrix(NTP,NRHOP);
101 table->d2hdtdrho = alloc_matrix(NTP,NRHOP);
102 }
103
104
105 void remove_tables(Ttse *table)
106 {
107
108 ASC_FREE(table->satFRho );
109 ASC_FREE(table->satFdRhodt );
110 ASC_FREE(table->satFd2RhodT2 );
111
112 ASC_FREE(table->satGRho );
113 ASC_FREE(table->satGdRhodt );
114 ASC_FREE(table->satGd2RhodT2 );
115
116 remove_matrix(table->dsdt,NTP);
117 remove_matrix(table->d2sdt2,NTP);
118 remove_matrix(table->dsdrho,NTP);
119 remove_matrix(table->d2sdrho2,NTP);
120 remove_matrix(table->d2sdtdrho,NTP);
121
122 remove_matrix(table->dpdt,NTP);
123 remove_matrix(table->d2pdt2,NTP);
124 remove_matrix(table->dpdrho,NTP);
125 remove_matrix(table->d2pdrho2,NTP);
126 remove_matrix(table->d2pdtdrho,NTP);
127
128
129 remove_matrix(table->dudt,NTP);
130 remove_matrix(table->d2udt2,NTP);
131 remove_matrix(table->dudrho,NTP);
132 remove_matrix(table->d2udrho2,NTP);
133 remove_matrix(table->d2udtdrho,NTP);
134
135 remove_matrix(table->dgdt,NTP);
136 remove_matrix(table->d2gdt2,NTP);
137 remove_matrix(table->dgdrho,NTP);
138 remove_matrix(table->d2gdrho2,NTP);
139 remove_matrix(table->d2gdtdrho,NTP);
140
141 remove_matrix(table->dhdt,NTP);
142 remove_matrix(table->d2hdt2,NTP);
143 remove_matrix(table->dhdrho,NTP);
144 remove_matrix(table->d2hdrho2,NTP);
145 remove_matrix(table->d2hdtdrho,NTP);
146
147 }
148
149
150
151 /*
152 This will load the binary file from tables/ for the liquid of interest and the EOS and populate the matrices.
153 If the files are not present in tables/ then build_tables() should be used.
154 */
155
156 void load_tables(PureFluid *P){
157
158
159 }
160
161 /*
162 After building the tables once this should be called to save the files in binary inside tables/
163 */
164 void save_tables(PureFluid *P){
165
166
167 }
168
169 /*
170 Actual building of tables is done here.
171 */
172 void build_tables(PureFluid *P){
173
174 #ifndef PT
175 #define PT P->table
176
177 int i,j;
178 FpropsError err = FPROPS_NO_ERROR;
179
180 double Tt = P->data->T_t;
181 double Tc = P->data->T_c;
182 double dt = (Tc - Tt)/(NSAT);
183
184 MSG("triple point and critical temperature --> %f %f",Tt,Tc);
185
186 for(i=0; i<NSAT; ++i)
187 {
188
189 double T = Tt + i*dt;
190
191 double rhof,rhog;
192 P->sat_fn(T,&rhof,&rhog,P->data,&err);
193
194
195 //The fluid saturation line rho's and 1st & 2nd derivatives of rho with respect to T
196
197 double dpdT_rho = P->dpdT_rho_fn(T,rhof,P->data,&err);
198 double dpdrho_T = P->dpdrho_T_fn(T,rhof,P->data,&err);
199
200 double d2pdrho2_T = P->d2pdrho2_T_fn(T,rhof,P->data,&err);
201 double d2pdrhodT = P->d2pdTdrho_fn(T,rhof,P->data,&err);
202 double d2pdT2_rho = P->d2pdT2_rho_fn(T,rhof,P->data,&err);
203
204
205 double ddrho_drhodT_p_constT = ( dpdT_rho*d2pdrho2_T - dpdrho_T*d2pdrhodT ) / pow(dpdrho_T,2);
206 double ddT_drhodT_p_constrho = ( dpdT_rho*d2pdrhodT - dpdrho_T*d2pdT2_rho ) / pow(dpdrho_T,2);
207
208 double drhodT_p = (-dpdT_rho )/(dpdrho_T);
209
210
211 PT->satFRho[i] = rhof;
212 PT->satFdRhodt[i] = drhodT_p;
213 PT->satFd2RhodT2[i] = ddT_drhodT_p_constrho + ddrho_drhodT_p_constT * drhodT_p;
214
215
216 //The Vapour saturation line rho's and 1st & 2nd derivatives of rho with respect to T
217
218 dpdT_rho = P->dpdT_rho_fn(T,rhog,P->data,&err);
219 dpdrho_T = P->dpdrho_T_fn(T,rhog,P->data,&err);
220
221 d2pdrho2_T = P->d2pdrho2_T_fn(T,rhog,P->data,&err);
222 d2pdrhodT = P->d2pdTdrho_fn(T,rhog,P->data,&err);
223 d2pdT2_rho = P->d2pdT2_rho_fn(T,rhog,P->data,&err);
224
225 ddrho_drhodT_p_constT = ( dpdT_rho*d2pdrho2_T - dpdrho_T*d2pdrhodT ) / pow(dpdrho_T,2);
226 ddT_drhodT_p_constrho = ( dpdT_rho*d2pdrhodT - dpdrho_T*d2pdT2_rho ) / pow(dpdrho_T,2);
227
228 drhodT_p = (-dpdT_rho )/(dpdrho_T);
229
230
231 PT->satGRho[i] = rhog;
232 PT->satGdRhodt[i] = drhodT_p;
233 PT->satGd2RhodT2[i] = ddT_drhodT_p_constrho + ddrho_drhodT_p_constT * drhodT_p;
234
235
236 }
237
238
239 double tmin,tmax,rhomin,rhomax;
240
241 //Pseudo values for water
242 //Should be implemented elsewhere per fluid
243 PT->tmin = 200;
244 PT->tmax = 4200;
245 PT->rhomin = 400;
246 PT->rhomax = 4400;
247
248
249
250 tmin = PT->tmin;
251 tmax = PT->tmax;
252 rhomin = PT->rhomin;
253 rhomax = PT->rhomax;
254
255 dt = (tmax-tmin)/NTP;
256 double drho = (rhomax-rhomin)/NRHOP;
257
258
259 MSG("DTemp is %f and DRho is %f",dt,drho);
260 MSG("BUILDING TABLES");
261
262 clock_t start = clock();
263
264 for( i = 0; i < NTP; i++)
265 for( j = 0; j < NRHOP; j++){
266
267 double t = tmin+i*dt;
268 double rho = rhomin+j*drho;
269
270 PT->p[i][j] = P->p_fn( t, rho , P->data, &err);
271 PT->dpdt[i][j] = P->dpdT_rho_fn( t, rho , P->data, &err);
272 PT->dpdrho[i][j] = P->dpdrho_T_fn( t, rho , P->data, &err);
273 #ifdef SO
274 PT->d2pdt2[i][j] = P->d2pdT2_rho_fn( t, rho , P->data, &err);
275 PT->d2pdrho2[i][j] = P->d2pdrho2_T_fn( t, rho , P->data, &err);
276 PT->d2pdtdrho[i][j] = P->d2pdTdrho_fn( t, rho , P->data, &err);
277 #endif
278
279 PT->h[i][j] = P->h_fn( t, rho , P->data, &err);
280 PT->dhdt[i][j] = P->dhdT_rho_fn( t, rho , P->data, &err);
281 PT->dhdrho[i][j] = P->dhdrho_T_fn( t, rho , P->data, &err);
282 #ifdef SO
283 PT->d2hdt2[i][j] = P->d2hdT2_rho_fn( t, rho , P->data, &err);
284 PT->d2hdrho2[i][j] = P->d2hdrho2_T_fn( t, rho , P->data, &err);
285 PT->d2hdtdrho[i][j] = P->d2hdTdrho_fn( t, rho , P->data, &err);
286 #endif
287
288 PT->s[i][j] = P->s_fn( t, rho , P->data, &err);
289 PT->dsdt[i][j] = P->dsdT_rho_fn( t, rho , P->data, &err);
290 PT->dsdrho[i][j] = P->dsdrho_T_fn( t, rho , P->data, &err);
291 #ifdef SO
292 PT->d2sdt2[i][j] = P->d2sdT2_rho_fn( t, rho , P->data, &err);
293 PT->d2sdrho2[i][j] = P->d2sdrho2_T_fn( t, rho , P->data, &err);
294 PT->d2sdtdrho[i][j] = P->d2sdTdrho_fn( t, rho , P->data, &err);
295 #endif
296 PT->u[i][j] = P->u_fn( t, rho , P->data, &err);
297 PT->dudt[i][j] = P->dudT_rho_fn( t, rho , P->data, &err);
298 PT->dudrho[i][j] = P->dudrho_T_fn( t, rho , P->data, &err);
299 #ifdef SO
300 PT->d2udt2[i][j] = P->d2udT2_rho_fn( t, rho , P->data, &err);
301 PT->d2udrho2[i][j] = P->d2udrho2_T_fn( t, rho , P->data, &err);
302 PT->d2udtdrho[i][j] = P->d2udTdrho_fn( t, rho , P->data, &err);
303 #endif
304
305 PT->g[i][j] = P->g_fn( t, rho , P->data, &err);
306 PT->dgdt[i][j] = P->dgdT_rho_fn( t, rho , P->data, &err);
307 PT->dgdrho[i][j] = P->dgdrho_T_fn( t, rho , P->data, &err);
308 #ifdef SO
309 PT->d2gdt2[i][j] = P->d2gdT2_rho_fn( t, rho , P->data, &err);
310 PT->d2gdrho2[i][j] = P->d2gdrho2_T_fn( t, rho , P->data, &err);
311 PT->d2gdtdrho[i][j] = P->d2gdTdrho_fn( t, rho , P->data, &err);
312 #endif
313 }
314
315
316 clock_t end = clock();
317 double msec = (double)(end - start) / (CLOCKS_PER_SEC/1000);
318 MSG("Tables built in %f seconds", msec/1000);
319
320
321
322 #undef PT
323 #endif
324 }
325
326
327 /*
328 Second Order Taylor series expansion
329 */
330 #ifdef SO
331 #define EVALTTSEFN(VAR) \
332 double evaluate_ttse_##VAR(PureFluid *P , double t, double rho){\
333 double rho_f, rho_g;\
334 FpropsError err;\
335 int i,j;\
336 double tmin = P->data->T_t;\
337 double tmax = P->data->T_c;\
338 if(t >= tmin && t< tmax) {\
339 double dt = (tmax-tmin)/NSAT;\
340 i = (int)round(((t - tmin)/(tmax - tmin)*(NSAT)));\
341 double delt = t - ( tmin + i*dt);\
342 rho_f = P->table->satFRho[i] + delt*P->table->satFdRhodt[i] + 0.5*delt*delt*P->table->satFd2RhodT2[i];\
343 rho_g = P->table->satGRho[i] + delt*P->table->satGdRhodt[i] + 0.5*delt*delt*P->table->satGd2RhodT2[i];\
344 if(rho_g < rho && rho < rho_f){\
345 double x = rho_g*(rho_f/rho - 1)/(rho_f - rho_g);\
346 double Qf = P->VAR##_fn( t,rho_f,P->data,&err);\
347 double Qg = P->VAR##_fn( t,rho_g,P->data,&err);\
348 return x*Qg + (1-x)*Qf;\
349 }\
350 }\
351 tmin = P->table->tmin;\
352 tmax = P->table->tmax;\
353 double rhomin = P->table->rhomin;\
354 double rhomax = P->table->rhomax;\
355 double dt = (tmax-tmin)/NTP;\
356 double drho = (rhomax-rhomin)/NRHOP;\
357 i = (int)round(((t-tmin)/(tmax-tmin)*(NTP)));\
358 j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP)));\
359 double delt = t - ( tmin + i*dt);\
360 double delrho = rho - ( rhomin + j*drho);\
361 double ttse##VAR = P->table->VAR[i][j]\
362 + delt*P->table->d##VAR##dt[i][j] + 0.5*delt*delt*P->table->d2##VAR##dt2[i][j]\
363 + delrho*P->table->d##VAR##drho[i][j] + 0.5*delrho*delrho*P->table->d2##VAR##drho2[i][j]\
364 + delrho*delt*P->table->d2##VAR##dtdrho[i][j];\
365 return ttse##VAR;\
366 }
367 #endif
368
369 /*
370 First Order Taylor series expansion
371 */
372 #ifdef FO
373 #define EVALTTSEFNFO(VAR) \
374 double evaluate_ttse_##VAR(PureFluid *P , double t, double rho){\
375 int i,j;\
376 double tmin = P->table->tmin; double tmax = P->table->tmax;\
377 double rhomin = P->table->rhomin; double rhomax= P->table->rhomax;\
378 double dt = (tmax-tmin)/NTP;\
379 double drho = (rhomax-rhomin)/NRHOP;\
380 i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));\
381 j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));\
382 double delt = t - ( tmin + i*dt);\
383 double delrho = rho - ( rhomin + j*drho);\
384 double ttse##VAR = P->table->VAR[i][j]\
385 + delt*P->table->d##VAR##dt[i][j] \
386 + delrho*P->table->d##VAR##drho[i][j] ;\
387 return ttse##VAR;\
388 }
389 #endif
390 //Second order accurate evaluation
391
392 #ifdef SO
393 EVALTTSEFN(p);
394 EVALTTSEFN(h);
395 EVALTTSEFN(s);
396 EVALTTSEFN(g);
397 EVALTTSEFN(u);
398 #endif
399
400
401 //First order accurate evaluation
402 #ifdef FO
403 EVALTTSEFNFO(p);
404 EVALTTSEFNFO(h);
405 EVALTTSEFNFO(s);
406 EVALTTSEFNFO(g);
407 EVALTTSEFNFO(u);
408 #endif
409
410 /*
411 double evaluate_ttse_p(PureFluid *P , double t, double rho)
412 {
413
414 #ifndef PT
415 #define PT P->table
416
417 int i,j;
418 FpropsError err = FPROPS_NO_ERROR;
419
420 double tmin = PT->tmin;
421 double tmax = PT->tmax;
422 double rhomin = PT->rhomin;
423 double rhomax= PT->rhomax;
424
425 double dt = (tmax-tmin)/NTP;
426 double drho = (rhomax-rhomin)/NRHOP;
427
428 i = (int)round(((t-tmin)/(tmax-tmin)*(NTP-1)));
429 j = (int)round(((rho-rhomin)/(rhomax-rhomin)*(NRHOP-1)));
430
431 double delt = t - ( tmin + i*dt);
432 double delrho = rho - ( rhomin + j*drho);
433
434 MSG("%f %f",i ,j ,delt ,delrho );
435 MSG("%f %f %f %f %f %f", PT->p[i][j], PT->dpdt[i][j], PT->d2pdt2[i][j], PT->dpdrho[i][j], PT->d2pdrho2[i][j], PT->d2pdtdrho[i][j]);
436
437
438 double ttse_p = PT->p[i][j] + delt*PT->dpdt[i][j] + 0.5*delt*delt*PT->d2pdt2[i][j]
439 + delrho*PT->dpdrho[i][j] + 0.5*delrho*delrho*PT->d2pdrho2[i][j]
440 + delrho*delt*PT->d2pdtdrho[i][j];
441
442
443 MSG("%e %e", ttse_p , P->p_fn(t, rho, P->data,&err) );
444
445
446 return(ttse_p)
447 #undef PT
448 #endif
449 }
450
451 */
452
453
454 void ttse_prepare(PureFluid *P){
455
456 #ifdef TTSE_DEBUG
457 //FILE *F1 = fopen("ttse.txt","w");
458 //fprintf(F1,"%f %f\n",t, P->p_fn( t, rho , P->data,&err) );
459 #endif
460
461
462 if(!P->table->usettse)
463 return;
464
465
466 MSG("Inside TTSE");
467
468 alloc_tables(P->table);
469
470 build_tables(P);
471
472 #ifdef TTSE_DEBUG
473 //fclose(F1);
474 #endif
475
476
477 }
478
479
480 void ttse_clean(PureFluid *P){
481 remove_tables(P->table);
482
483
484 }
485

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