johnpye/fprops/water.c

/*  ASCEND modelling environment
    Copyright (C) 2008 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*/

#include "water.h"

#define WATER_R 461.51805 /* J/kg路K */
#define WATER_TSTAR 647.096 /* K */

/**
Ideal gas data for Water/Steam, from IAPWS-95.
http://www.iapws.org/relguide/IAPWS95.pdf
*/
const IdealData ideal_data_water = {
    -8.32044648201 /* constant */
    , 6.6832105268 /* linear */
    , WATER_TSTAR /* Tstar */
    , WATER_R /* cpstar J/kgK */
    , 1 /* power terms */   
    , (const IdealPowTerm[]){
        {1. + 3.00632,  0}
    }
    , 5 /* exponential terms */
    , (const IdealExpTerm []){
        /* b, beta */
        {0.012436, 1.28728967*WATER_TSTAR}
        ,{0.97315,  3.53734222*WATER_TSTAR}
        ,{1.27950, 7.74073708*WATER_TSTAR}
        ,{0.96956, 9.24437796*WATER_TSTAR}
        ,{0.24873, 27.5075105*WATER_TSTAR}
    }
};


/**
Residual (non-ideal) property data for Water/Steam, from IAPWS-95.
http://www.iapws.org/relguide/IAPWS95.pdf
*/
const HelmholtzData helmholtz_data_water = {
    /* R */ WATER_R /* J/kg/K */
    , /* M */ 0.00000000000000 /* kg/kmol  -- need to look up value cited by IAPWS */
    , /* rho_star */322. /* kg/m鲁 */
    , /* T_star */ WATER_TSTAR /* K */
    , &ideal_data_water
    , 51 /* np */
    , (const HelmholtzPowTerm[]){
        /* a_i, t_i, d_i, l_i */
        {0.12533547935523E-1, -0.5, 1, 0}
        ,{0.78957634722828E1, 0.875, 1, 0}
        ,{-0.87803203303561E1, 1, 1, 0}
        ,{0.31802509345418, 0.5, 2, 0}
        ,{-0.26145533859358, 0.75, 2, 0}
        ,{-0.78199751687981E-2, 0.375, 3, 0}
        ,{0.88089493102134E-2, 1, 4, 0}
        ,{-0.66856572307965, 4, 1, 1}
        ,{0.20433810950965, 6, 1, 1}
        ,{-0.66212605039687E-4, 12, 1, 1}
        ,{-0.19232721156002, 1, 2, 1}
        ,{-0.25709043003438, 5, 2, 1}
        ,{0.16074868486251, 4, 3, 1}
        ,{-0.40092828925807E-1, 2, 4, 1}
        ,{0.39343422603254E-6, 13, 4, 1}
        ,{-0.75941377088144E-5, 9, 5, 1}
        ,{0.56250979351888E-3, 3, 7, 1}
        ,{-0.15608652257135E-4, 4, 9, 1}
        ,{0.11537996422951E-8, 11, 10, 1}
        ,{0.36582165144204E-6, 4, 11, 1}
        ,{-0.13251180074668E-11, 13, 13, 1}
        ,{-0.62639586912454E-9, 1, 15, 1}
        ,{-0.10793600908932, 7, 1, 2}
        ,{0.17611491008752E-1, 1, 2, 2}
        ,{0.22132295167546, 9, 2, 2}
        ,{-0.40247669763528, 10, 2, 2}
        ,{0.58083399985759, 10, 3, 2}
        ,{0.49969146990806E-2, 3, 4, 2}
        ,{-0.31358700712549E-1, 7, 4, 2}
        ,{-0.74315929710341, 10, 4, 2}
        ,{0.47807329915480, 10, 5, 2}
        ,{0.20527940895948E-1, 6, 6, 2}
        ,{-0.13636435110343, 10, 6, 2}
        ,{0.14180634400617E-1, 10, 7, 2}
        ,{0.83326504880713E-2, 1, 9, 2}
        ,{-0.29052336009585E-1, 2, 9, 2}
        ,{0.38615085574206E-1, 3, 9, 2}
        ,{-0.20393486513704E-1, 4, 9, 2}
        ,{-0.16554050063734E-2, 8, 9, 2}
        ,{0.19955571979541E-2, 6, 10, 2}
        ,{0.15870308324157E-3, 9, 10, 2}
        ,{-0.16388568342530E-4, 8, 12, 2}
        ,{0.43613615723811E-1, 16, 3, 3}
        ,{0.34994005463765E-1, 22, 4, 3}
        ,{-0.76788197844621E-1, 23, 4, 3}
        ,{0.22446277332006E-1, 23, 5, 3}
        ,{-0.62689710414685E-4, 10, 14, 4}
        ,{-0.55711118565645E-9, 50, 3, 6}
        ,{-0.19905718354408, 44, 6, 6}
        ,{0.31777497330738, 46, 6, 6}
        ,{-0.11841182425981, 50, 6, 6}
    }
    , 3 /* gaussian terms */
    , (const HelmholtzGausTerm[]){
        /* n, t, d, alpha, beta, gamma, epsilon */
        {-0.31306260323435E2, 0, 3, 20, 150, 1.21, 1}
        ,{0.31546140237781E2, 1, 3, 20, 150, 1.21, 1}
        ,{-0.25213154341695E4,4, 3, 20, 250, 1.25, 1}
    }
    , 2 /* critical terms */
    , (const HelmholtzCritTerm[]){
        /* n, a, b, beta, A, B, C, D */
        {-0.14874640856724, 3.5, 0.85, 0.3, 0.32, 0.2, 28, 700}
        ,{0.31806110878444, 3.5, 0.95, 0.3, 0.32, 0.2, 32, 800}
    }
};


/*
    Test suite. These tests attempt to validate the current code using 
    a few sample figures output by REFPROP 7.0.

    To run the test, compile and run as follows:

    ./test.py water
*/
#ifdef TEST

/* 
    some code from freesteam, http://freesteam.sf.net/, which has been thoroughly
    validated already.
*/

const double n0[] = {
    0.0 /* placeholder */,
    -8.32044648201, 6.6832105268, 3.00632 /* const, linear, ln(tau) coeffs */
    , 0.012436, 0.97315, 1.27950, 0.96956, 0.24873 /* exponential coeffs */
};

const double gamma0[] = {
    0.0, 0.0, 0.0, 0.0,
    1.28728967,
    3.53734222,
    7.74073708,
    9.24437796,
    27.5075105
};

enum Limits{
    eGamma1 = 4,
    eGamma2 = 9,
};

#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ideal_impl.h"
#include "helmholtz_impl.h"

double phi0(const double delta, const double tau){
    int i;
    double sum = 0;
    for (i = eGamma1; i < eGamma2; i++)
    {
        sum += n0[i]*log(1-exp(-tau*gamma0[i]));
    }
    sum += log(delta) + n0[1] + n0[2]*tau + n0[3]*log(tau);
    return sum;
}

int main(void){
    double rho, T;
    const HelmholtzData *d;

    d = &helmholtz_data_water;
    double maxerr = 0;

/* a simple macro to actually do the testing */
#define ASSERT_TOL(FN,PARAM1,PARAM2,PARAM3,VAL,TOL) {\
        double cval; cval = FN(PARAM1,PARAM2,PARAM3);\
        double err; err = cval - (double)(VAL);\
        double relerrpc = (cval-(VAL))/(VAL)*100;\
        if(fabs(relerrpc)>maxerr)maxerr=fabs(relerrpc);\
        if(fabs(err)>fabs(TOL)){\
            fprintf(stderr,"ERROR in line %d: value of '%s(%f,%f,%s)' = %0.8f,"\
                " should be %f, error is %.10e (%.2f%%)!\n"\
                , __LINE__, #FN,PARAM1,PARAM2,#PARAM3, cval, VAL,cval-(VAL)\
                ,relerrpc\
            );\
            exit(1);\
        }else{\
            fprintf(stderr,"    OK, %s(%f,%f,%s) = %8.2e with %.8f%% err.\n"\
                ,#FN,PARAM1,PARAM2,#PARAM3,VAL,relerrpc\
            );\
        }\
    }


    fprintf(stderr,"phi0 TESTS\n");
    for(T = 300; T <= 900; T+= 100){
        for(rho = 900; rho >= 0.9; rho*=0.5){
            double delta = rho / d->rho_star;
            double tau = d->T_star / T;
            double p0 = phi0(delta,tau);
            
            ASSERT_TOL(helm_ideal, tau, delta, d->ideal, p0, p0*1e-5);
        }
    }

    fprintf(stderr,"IAPWS95 TABLE 6 TESTS\n");
    T = 500.; /* K */
    rho = 838.025; /* kg/m鲁 */
    double tau = d->T_star / T;
    double delta = rho / d->rho_star;

    ASSERT_TOL(helm_ideal, tau, delta, d->ideal, 0.204797733E1, 1e-6);
    ASSERT_TOL(helm_ideal_tau, tau, delta, d->ideal, 0.904611106E1, 1e-6);
    ASSERT_TOL(HELM_IDEAL_DELTAU, tau, delta, d->ideal, 0., 1e-6);
    
    double phitt = helm_ideal_tautau(tau, d->ideal);
    double val = (-0.193249185E1);
    double err = phitt - val;
    if(fabs(err) > 1e-6){
        fprintf(stderr,"ERROR in helm_ideal_tautau near line %d\n",__LINE__);
        exit(1);
    }else{
        fprintf(stderr,"    OK, helm_ideal_tautau(%f) = %8.2e with %.6f%% err.\n"
            ,tau,val,err/val*100.
        );
    }

    /* FIXME still need to implement helm_ideal_del, helm_ideal_deldel, helm_ideal_deltau */

    ASSERT_TOL(helm_resid, tau, delta, d, -0.342693206E1, 1e-8);
    ASSERT_TOL(helm_resid_del, tau, delta, d, -0.364366650, 1e-8);
    ASSERT_TOL(helm_resid_deldel, tau, delta, d, 0.856063701, 1e-8);
    ASSERT_TOL(helm_resid_tau, tau, delta, d, -0.581403435E1, 1e-8);
    ASSERT_TOL(helm_resid_tautau, tau, delta, d, -0.223440737E1, 1e-8);
    ASSERT_TOL(helm_resid_deltau, tau, delta, d, -0.112176915e1, 1e-8);

    fprintf(stderr,"Tests completed OK (maximum error = %0.2f%%)\n",maxerr);
    exit(0);
}

#endif

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 "water.h"
21
22	#define WATER_R 461.51805 /* J/kg路K */
23	#define WATER_TSTAR 647.096 /* K */
24
25	/**
26	Ideal gas data for Water/Steam, from IAPWS-95.
27	http://www.iapws.org/relguide/IAPWS95.pdf
28	*/
29	const IdealData ideal_data_water = {
30	-8.32044648201 /* constant */
31	, 6.6832105268 /* linear */
32	, WATER_TSTAR /* Tstar */
33	, WATER_R /* cpstar J/kgK */
34	, 1 /* power terms */
35	, (const IdealPowTerm[]){
36	{1. + 3.00632, 0}
37	}
38	, 5 /* exponential terms */
39	, (const IdealExpTerm []){
40	/* b, beta */
41	{0.012436, 1.28728967*WATER_TSTAR}
42	,{0.97315, 3.53734222*WATER_TSTAR}
43	,{1.27950, 7.74073708*WATER_TSTAR}
44	,{0.96956, 9.24437796*WATER_TSTAR}
45	,{0.24873, 27.5075105*WATER_TSTAR}
46	}
47	};
48
49
50	/**
51	Residual (non-ideal) property data for Water/Steam, from IAPWS-95.
52	http://www.iapws.org/relguide/IAPWS95.pdf
53	*/
54	const HelmholtzData helmholtz_data_water = {
55	/* R / WATER_R / J/kg/K */
56	, /* M / 0.00000000000000 / kg/kmol -- need to look up value cited by IAPWS */
57	, /* rho_star /322. / kg/m鲁 */
58	, /* T_star / WATER_TSTAR / K */
59	, &ideal_data_water
60	, 51 /* np */
61	, (const HelmholtzPowTerm[]){
62	/* a_i, t_i, d_i, l_i */
63	{0.12533547935523E-1, -0.5, 1, 0}
64	,{0.78957634722828E1, 0.875, 1, 0}
65	,{-0.87803203303561E1, 1, 1, 0}
66	,{0.31802509345418, 0.5, 2, 0}
67	,{-0.26145533859358, 0.75, 2, 0}
68	,{-0.78199751687981E-2, 0.375, 3, 0}
69	,{0.88089493102134E-2, 1, 4, 0}
70	,{-0.66856572307965, 4, 1, 1}
71	,{0.20433810950965, 6, 1, 1}
72	,{-0.66212605039687E-4, 12, 1, 1}
73	,{-0.19232721156002, 1, 2, 1}
74	,{-0.25709043003438, 5, 2, 1}
75	,{0.16074868486251, 4, 3, 1}
76	,{-0.40092828925807E-1, 2, 4, 1}
77	,{0.39343422603254E-6, 13, 4, 1}
78	,{-0.75941377088144E-5, 9, 5, 1}
79	,{0.56250979351888E-3, 3, 7, 1}
80	,{-0.15608652257135E-4, 4, 9, 1}
81	,{0.11537996422951E-8, 11, 10, 1}
82	,{0.36582165144204E-6, 4, 11, 1}
83	,{-0.13251180074668E-11, 13, 13, 1}
84	,{-0.62639586912454E-9, 1, 15, 1}
85	,{-0.10793600908932, 7, 1, 2}
86	,{0.17611491008752E-1, 1, 2, 2}
87	,{0.22132295167546, 9, 2, 2}
88	,{-0.40247669763528, 10, 2, 2}
89	,{0.58083399985759, 10, 3, 2}
90	,{0.49969146990806E-2, 3, 4, 2}
91	,{-0.31358700712549E-1, 7, 4, 2}
92	,{-0.74315929710341, 10, 4, 2}
93	,{0.47807329915480, 10, 5, 2}
94	,{0.20527940895948E-1, 6, 6, 2}
95	,{-0.13636435110343, 10, 6, 2}
96	,{0.14180634400617E-1, 10, 7, 2}
97	,{0.83326504880713E-2, 1, 9, 2}
98	,{-0.29052336009585E-1, 2, 9, 2}
99	,{0.38615085574206E-1, 3, 9, 2}
100	,{-0.20393486513704E-1, 4, 9, 2}
101	,{-0.16554050063734E-2, 8, 9, 2}
102	,{0.19955571979541E-2, 6, 10, 2}
103	,{0.15870308324157E-3, 9, 10, 2}
104	,{-0.16388568342530E-4, 8, 12, 2}
105	,{0.43613615723811E-1, 16, 3, 3}
106	,{0.34994005463765E-1, 22, 4, 3}
107	,{-0.76788197844621E-1, 23, 4, 3}
108	,{0.22446277332006E-1, 23, 5, 3}
109	,{-0.62689710414685E-4, 10, 14, 4}
110	,{-0.55711118565645E-9, 50, 3, 6}
111	,{-0.19905718354408, 44, 6, 6}
112	,{0.31777497330738, 46, 6, 6}
113	,{-0.11841182425981, 50, 6, 6}
114	}
115	, 3 /* gaussian terms */
116	, (const HelmholtzGausTerm[]){
117	/* n, t, d, alpha, beta, gamma, epsilon */
118	{-0.31306260323435E2, 0, 3, 20, 150, 1.21, 1}
119	,{0.31546140237781E2, 1, 3, 20, 150, 1.21, 1}
120	,{-0.25213154341695E4,4, 3, 20, 250, 1.25, 1}
121	}
122	, 2 /* critical terms */
123	, (const HelmholtzCritTerm[]){
124	/* n, a, b, beta, A, B, C, D */
125	{-0.14874640856724, 3.5, 0.85, 0.3, 0.32, 0.2, 28, 700}
126	,{0.31806110878444, 3.5, 0.95, 0.3, 0.32, 0.2, 32, 800}
127	}
128	};
129
130
131	/*
132	Test suite. These tests attempt to validate the current code using
133	a few sample figures output by REFPROP 7.0.
134
135	To run the test, compile and run as follows:
136
137	./test.py water
138	*/
139	#ifdef TEST
140
141	/*
142	some code from freesteam, http://freesteam.sf.net/, which has been thoroughly
143	validated already.
144	*/
145
146	const double n0[] = {
147	0.0 /* placeholder */,
148	-8.32044648201, 6.6832105268, 3.00632 /* const, linear, ln(tau) coeffs */
149	, 0.012436, 0.97315, 1.27950, 0.96956, 0.24873 /* exponential coeffs */
150	};
151
152	const double gamma0[] = {
153	0.0, 0.0, 0.0, 0.0,
154	1.28728967,
155	3.53734222,
156	7.74073708,
157	9.24437796,
158	27.5075105
159	};
160
161	enum Limits{
162	eGamma1 = 4,
163	eGamma2 = 9,
164	};
165
166	#include <assert.h>
167	#include <stdlib.h>
168	#include <stdio.h>
169	#include <math.h>
170	#include "ideal_impl.h"
171	#include "helmholtz_impl.h"
172
173	double phi0(const double delta, const double tau){
174	int i;
175	double sum = 0;
176	for (i = eGamma1; i < eGamma2; i++)
177	{
178	sum += n0[i]log(1-exp(-taugamma0[i]));
179	}
180	sum += log(delta) + n0[1] + n0[2]tau + n0[3]log(tau);
181	return sum;
182	}
183
184	int main(void){
185	double rho, T;
186	const HelmholtzData *d;
187
188	d = &helmholtz_data_water;
189	double maxerr = 0;
190
191	/* a simple macro to actually do the testing */
192	#define ASSERT_TOL(FN,PARAM1,PARAM2,PARAM3,VAL,TOL) {\
193	double cval; cval = FN(PARAM1,PARAM2,PARAM3);\
194	double err; err = cval - (double)(VAL);\
195	double relerrpc = (cval-(VAL))/(VAL)*100;\
196	if(fabs(relerrpc)>maxerr)maxerr=fabs(relerrpc);\
197	if(fabs(err)>fabs(TOL)){\
198	fprintf(stderr,"ERROR in line %d: value of '%s(%f,%f,%s)' = %0.8f,"\
199	" should be %f, error is %.10e (%.2f%%)!\n"\
200	, __LINE__, #FN,PARAM1,PARAM2,#PARAM3, cval, VAL,cval-(VAL)\
201	,relerrpc\
202	);\
203	exit(1);\
204	}else{\
205	fprintf(stderr," OK, %s(%f,%f,%s) = %8.2e with %.8f%% err.\n"\
206	,#FN,PARAM1,PARAM2,#PARAM3,VAL,relerrpc\
207	);\
208	}\
209	}
210
211
212	fprintf(stderr,"phi0 TESTS\n");
213	for(T = 300; T <= 900; T+= 100){
214	for(rho = 900; rho >= 0.9; rho*=0.5){
215	double delta = rho / d->rho_star;
216	double tau = d->T_star / T;
217	double p0 = phi0(delta,tau);
218
219	ASSERT_TOL(helm_ideal, tau, delta, d->ideal, p0, p0*1e-5);
220	}
221	}
222
223	fprintf(stderr,"IAPWS95 TABLE 6 TESTS\n");
224	T = 500.; /* K */
225	rho = 838.025; /* kg/m鲁 */
226	double tau = d->T_star / T;
227	double delta = rho / d->rho_star;
228
229	ASSERT_TOL(helm_ideal, tau, delta, d->ideal, 0.204797733E1, 1e-6);
230	ASSERT_TOL(helm_ideal_tau, tau, delta, d->ideal, 0.904611106E1, 1e-6);
231	ASSERT_TOL(HELM_IDEAL_DELTAU, tau, delta, d->ideal, 0., 1e-6);
232
233	double phitt = helm_ideal_tautau(tau, d->ideal);
234	double val = (-0.193249185E1);
235	double err = phitt - val;
236	if(fabs(err) > 1e-6){
237	fprintf(stderr,"ERROR in helm_ideal_tautau near line %d\n",__LINE__);
238	exit(1);
239	}else{
240	fprintf(stderr," OK, helm_ideal_tautau(%f) = %8.2e with %.6f%% err.\n"
241	,tau,val,err/val*100.
242	);
243	}
244
245	/* FIXME still need to implement helm_ideal_del, helm_ideal_deldel, helm_ideal_deltau */
246
247	ASSERT_TOL(helm_resid, tau, delta, d, -0.342693206E1, 1e-8);
248	ASSERT_TOL(helm_resid_del, tau, delta, d, -0.364366650, 1e-8);
249	ASSERT_TOL(helm_resid_deldel, tau, delta, d, 0.856063701, 1e-8);
250	ASSERT_TOL(helm_resid_tau, tau, delta, d, -0.581403435E1, 1e-8);
251	ASSERT_TOL(helm_resid_tautau, tau, delta, d, -0.223440737E1, 1e-8);
252	ASSERT_TOL(helm_resid_deltau, tau, delta, d, -0.112176915e1, 1e-8);
253
254	fprintf(stderr,"Tests completed OK (maximum error = %0.2f%%)\n",maxerr);
255	exit(0);
256	}
257
258	#endif
259