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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1810 - (show annotations) (download) (as text)
Tue Jul 29 07:45:48 2008 UTC (12 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 2558 byte(s)
Added new code for calculating properties according to modified BWR correlation, also some initial
work on Helmholtz free energy correlation.
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 Modified Benedict-Webb-Rubin (MBWR) equation of state.
20
21 John Pye, 29 Jul 2008.
22 */
23
24 #include "mbwr.h"
25
26 #include <math.h>
27
28 /**
29 Function to calculate pressure from MBWR correlation, given temperature
30 and molar density.
31
32 @param T temperature in K
33 @param rhob molar density in mol/L
34 @return pressure in Pa
35 */
36 double mbwr_p(double T, double rhob, MbwrData *data){
37 int i;
38 double p = 0;
39 double Ti, Ti2, Ti3, Ti4;
40 double rhob2, rhobpow, sum10, rhob_r;
41 double alpha[15];
42
43 /* precalculate powers of T^-1 for faster evaluation */
44 Ti =1. / T;
45 Ti2 = Ti*Ti;
46 Ti3 = Ti2*Ti;
47 Ti4 = Ti3*Ti;
48
49 /* values of alpha in MBWR are functions of temperature */
50 double *B = data->beta;
51 #define R data->R
52 alpha[0] = R * T;
53 alpha[1] = B[0]*T + B[1]*sqrt(T) + B[2] + B[3]*Ti + B[4]*Ti2;
54 alpha[2] = B[5]*T + B[6] + B[7]*Ti + B[8]*Ti2;
55 alpha[3] = B[9]*T + B[10] + B[11]*Ti;
56 alpha[4] = B[12];
57 alpha[5] = B[13]*Ti + B[14]*Ti2;
58 alpha[6] = B[15]*Ti;
59 alpha[7] = B[16]*Ti + B[17]*Ti2;
60 alpha[8] = B[18]*Ti2;
61 alpha[9] = B[19]*Ti2 + B[20]*Ti3;
62 alpha[10] = B[21]*Ti2 + B[22]*Ti4;
63 alpha[11] = B[23]*Ti2 + B[24]*Ti3;
64 alpha[12] = B[25]*Ti2 + B[26]*Ti4;
65 alpha[13] = B[27]*Ti2 + B[28]*Ti3;
66 alpha[14] = B[29]*Ti2 + B[30]*Ti3 + B[31]*Ti4;
67 #undef R
68
69 /* add up the first sum in the MBWR correlation */
70 rhobpow = 1;
71 for(i=0;i<9;++i){
72 rhobpow *= rhob;
73 p += alpha[i]* rhobpow;
74 }
75
76 /* work out the second sum in the MBWR correlation */
77 sum10 = 0;
78 rhobpow = rhob;
79 rhob2 = rhob*rhob;
80 for(i=9; i<15; ++i){
81 rhobpow *= rhob2;
82 sum10 += alpha[i] * rhobpow;
83 }
84
85 /* calculate the exponential term and add it to 'p' */
86 #define RHO_C data->rhob_c
87 rhob_r = rhob/RHO_C;
88 p += exp(rhob_r*rhob_r) * sum10;
89 #undef RHO_C
90
91 return p * 1e-5; /* convert bar to Pa on return */
92 }
93
94

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