johnpye/fprops/helmholtz.c

/*  ASCEND modelling environment
    Copyright (C) 2006 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.
*//** @file
    Implementation of the reduced molar Helmholtz free energy equation of state.

    For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
    Fundamentalgleichung f端r Ammoniak.

    John Pye, 29 Jul 2008.
*/

#include <math.h>

#include "helmholtz.h"

/* forward decls */

static double helm_ideal(double tau, double delta, HelmholtzData *data);
static double helm_ideal_tau(double tau, double delta, HelmholtzData *data);
static double helm_resid(double tau, double delta, HelmholtzData *data);
static double helm_resid_del(double tau, double delta, HelmholtzData *data);
static double helm_resid_tau(double tau, double delta, HelmholtzData *data);

/**
    Function to calculate pressure from Helmholtz free energy EOS, given temperature
    and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return pressure in Pa???
*/
double helmholtz_p(double T, double rho, HelmholtzData *data){
    
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

    return data->R * T * rho * (1. + delta * helm_resid_del(tau,delta,data));
}

/**
    Function to calculate internal energy from Helmholtz free energy EOS, given
    temperature and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return internal energy in ???
*/
double helmholtz_u(double T, double rho, HelmholtzData *data){
    
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

    return data->R * T * tau * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data));
}

/**
    Function to calculate enthalpy from Helmholtz free energy EOS, given
    temperature and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return enthalpy in ????
*/
double helmholtz_h(double T, double rho, HelmholtzData *data){
    
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

    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));
}

/*---------------------------------------------
  IDEAL COMPONENT RELATIONS
*/

/**
    Ideal component of helmholtz function
*/  
double helm_ideal(double tau, double delta, HelmholtzData *data){
    double *a0;

    double tau13 = pow(tau,1./3.);
    double taum32 = pow(tau,-3./2.);
    double taum74 = pow(tau,-7./4.);

    a0 = &(data->a0[0]);
    return log(delta) + a0[0] + a0[1] * tau - log(tau) + a0[2] * tau13 + a0[3] * taum32 + a0[4] * taum74;
}

/**
    Partial dervivative of ideal component of helmholtz residual function with 
    respect to tau.
*/  
double helm_ideal_tau(double tau, double delta, HelmholtzData *data){
    double *a0;

    double taum114 = pow(tau,-11./4.);
    double taum74 = tau * taum114;
    double taum23 = pow(tau,-2./3.);

    a0 = &(data->a0[0]);

    return a0[1] - 1./tau + a0[2]/3.* taum23 - a0[3] * 3./4. * taum74 - a0[4] * 7./4. * taum114;
}   

/**
    Residual part of helmholtz function. Note: we have NOT prematurely
    optimised here ;-)
*/
double helm_resid(double tau, double delta, HelmholtzData *data){
    
    double sum;
    double phir = 0;
    unsigned i;

    HelmholtzATD *atd = &(data->atd[0]);
    
    for(i=0; i<5; ++i){
        phir += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
        ++atd;
    }

    sum = 0;
    for(i=5; i<10; ++i){
        sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
        ++atd;
    }
    phir += exp(-delta) * sum;

    sum = 0; 
    for(i=10; i<17; ++i){
        sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
        ++atd;
    }
    phir += exp(-delta*delta) * sum;

    sum = 0;
    for(i=17; i<21; ++i){
        sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d);
        ++atd;
    }
    phir += exp(-delta*delta*delta) * sum;

    return phir;
}

/**
    Derivative of the helmholtz residual function with respect to
    delta.
*/  
double helm_resid_del(double tau,double delta,HelmholtzData *data){
    
    double sum;
    double phir = 0;
    unsigned i;
    double XdelX;

    HelmholtzATD *atd = &(data->atd[0]);
    
    for(i=0; i<5; ++i){
        phir += atd->a * pow(tau, atd->t) * pow(delta, atd->d - 1) * atd->d;
        ++atd;
    }

    sum = 0;
    XdelX = delta;
    for(i=5; i<10; ++i){
        sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d) * (atd->d - XdelX);
    }
    phir += exp(-delta) * sum;

    sum = 0; 
    XdelX = 2*delta*delta;
    for(i=10; i<17; ++i){
        sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d) * (atd->d - XdelX);
    }
    phir += exp(-delta*delta) * sum;

    sum = 0;
    XdelX = 3*delta*delta*delta;
    for(i=17; i<21; ++i){
        sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d) * (atd->d - XdelX);
    }
    phir += exp(-delta*delta*delta) * sum;

    return phir;
}

/**
    Derivative of the helmholtz residual function with respect to
    tau.
*/          
double helm_resid_tau(double tau,double delta,HelmholtzData *data){
    
    double sum;
    double phir = 0;
    unsigned i;

    HelmholtzATD *atd = &(data->atd[0]);
    
    for(i=0; i<5; ++i){
        if(atd->t != 0){
            phir += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
        }
        ++atd;
    }

    sum = 0;
    for(i=5; i<10; ++i){
        if(atd->t != 0){
            sum += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
        }
        ++atd;
    }
    phir += exp(-delta) * sum;

    sum = 0; 
    for(i=10; i<17; ++i){
        if(atd->t != 0){
            sum += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
        }
        ++atd;
    }
    phir += exp(-delta*delta) * sum;

    sum = 0;
    for(i=17; i<21; ++i){
        if(atd->t != 0){
            sum += atd->a * pow(tau, atd->t - 1) * pow(delta, atd->d) * atd->t;
        }
        ++atd;
    }
    phir += exp(-delta*delta*delta) * sum;

    return phir;
}   


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(-deltadelta) 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(-deltadeltadelta) * 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 = 2deltadelta;
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(-deltadelta) sum;
193
194	sum = 0;
195	XdelX = 3deltadelta*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(-deltadeltadelta) * 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(-deltadelta) 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(-deltadeltadelta) * sum;
249
250	return phir;
251	}
252
253
254
255