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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2774 - (show annotations) (download) (as text)
Thu Mar 27 11:42:30 2014 UTC (3 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 9285 byte(s)
thcond and mu accessible from Python
1 /* ASCEND modelling environment
2 Copyright (C) 2014 John Pye
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, see <http://www.gnu.org/licenses/>.
16 */
17
18
19 #include "thcond.h"
20 #include "visc.h"
21 #include <math.h>
22
23 #ifndef PI
24 # define PI M_PI
25 #endif
26
27 #define K_BOLTZMANN 1.3806488e-23
28
29 //#define THCOND_DEBUG
30 #ifdef THCOND_DEBUG
31 # include "color.h"
32 # include "test.h"
33 # define MSG FPROPS_MSG
34 # define ERRMSG FPROPS_ERRMSG
35 #else
36 //# define ASSERT(ARGS...) (void)0)
37 # define MSG(ARGS...) ((void)0)
38 # define ERRMSG(ARGS...) ((void)0)
39 #endif
40
41 void thcond_prepare(PureFluid *P, const ThermalConductivityData *K, FpropsError *err){
42 MSG("Preparing thermal conductivity: currently we are just reusing the FileData pointer; no changes");
43 MSG("K.type: %d",K->type);
44 MSG("K.source: %s",K->source);
45 switch(K->type){
46 case FPROPS_THCOND_1:
47 MSG("K.data.k1.nc: %d",K->data.k1.nc);
48 P->thcond = K;
49 MSG("P.thcond.type = %d",P->thcond->type);
50 return;
51 case FPROPS_THCOND_NONE:
52 *err = FPROPS_NOT_IMPLEMENTED;
53 return;
54 }
55 }
56
57 /*-----------------------------IDEAL PART ------------------------------------*/
58
59 static double thcond1_cs(const ThermalConductivityData1 *K, double Tstar){
60 double res = 0;
61 int i;
62 MSG("K: %p",K);
63 MSG("Tstar = %f (1/Tstar = %f)",Tstar,1/Tstar);
64 MSG("nc = %d",K->nc);
65 for(i=0; i < K->nc; ++i){
66 //MSG("b[%d] = %e, i = %d, term = %f",i,K->ct[i].b, K->ct[i].i, K->ct[i].b * pow(Tstar, K->ct[i].i));
67 res += K->ct[i].b * pow(Tstar, K->ct[i].i);
68 }
69 MSG("res = %f",res);
70 return res;
71 }
72
73 double thcond1_lam0(FluidState state, FpropsError *err){
74 if(state.fluid->thcond->type != FPROPS_THCOND_1){*err = FPROPS_INVALID_REQUEST; return NAN;}
75 const ThermalConductivityData1 *k1 = &(state.fluid->thcond->data.k1);
76 double lam0 = 0;
77
78 MSG("k1: %p",k1);
79
80 // TODO FIXME need to re-factor this to be standardised and only use data from filedata.h structures.
81
82 if(0==strcmp(state.fluid->name,"carbondioxide")){
83 MSG("lam0 for carbondioxide");
84
85 //#define USE_CP0_FOR_LAM0
86 #ifdef USE_CP0_FOR_LAM0
87 double cp0 = fprops_cp0(state,err);
88 double R = state.fluid->data->R;
89 //MSG("cp0 = %f, R = %f", cp0, R);
90 double M = state.fluid->data->M;
91 double sigma = 0.3751; // nm!
92 double opr2 = 0.177568/0.475598 * cp0/R * 1 / SQ(sigma) / sqrt(M);
93 //MSG("1 + r^2 = %f (by cp0)", opr2_2);
94 //MSG("5/2 *(cp0(T)/R - 1) = %f\n", 5./2*(fprops_cp0(state,err)/state.fluid->data->R - 1));
95 #else
96 int i;
97 double sum1 = 0;
98 double c[] = {2.387869e-2, 4.350794, -10.33404, 7.981590, -1.940558};
99 for(i=0; i<5; ++i){
100 sum1 += c[i] * pow(state.T/100, 2-(i+1));
101 }
102 double cint_over_k = 1.0 + exp(-183.5/state.T)*sum1;
103
104 //MSG("cint/k = %f",cint_over_k);
105 //MSG("1 + r^2 = %f (by cint/k)",1+0.4*cint_over_k);
106 double opr2 = 1+0.4*cint_over_k;
107 //double r = sqrt(0.4*cint_over_k);
108 //MSG("r = %f",r);
109 #endif
110 // FIXME convert the other way, convert the cint/k stuff to simply cp0/R, should be more generalised that way.
111
112 double CS_star = thcond1_cs(k1, k1->eps_over_k/state.T);
113 //MSG("CS_star = %f", CS_star);
114
115 lam0 = 0.475598 * sqrt(state.T) * opr2 / CS_star;
116
117 // 0.177568 (mW/m/K)*(nm^2)
118 //lam0 = 0.177568 * sqrt(state.T) / sqrt(state.fluid->data->M) / SQ(sigma) * cp0 / R / CS_star;
119
120 }else if(0==strcmp(state.fluid->name,"nitrogen")){
121 // this uses a rather different approach/formulation; see http://ascend4.org/FPROPS/Thermal_conductivity
122 MSG("lam0 for nitrogen");
123 double N1 = 1.511;
124 double N2 = 2.117, t2 = -1.;
125 double N3 = -3.332, t3 = -0.7;
126 double tau = k1->T_star / state.T;
127 lam0 = N1 * (visc1_mu0(state,err)/1e-6) + N2 * pow(tau,t2) + N3 * pow(tau,t3);
128 }else{
129 ERRMSG("lam0 not implemented");
130 *err = FPROPS_NOT_IMPLEMENTED;
131 return 0;
132 }
133 MSG("lam0(T=%f) = %e",state.T, lam0);
134 return lam0 * k1->k_star;
135 }
136
137
138 /*---------------------------RESIDUAL PART -----------------------------------*/
139
140 double thcond1_lamr(FluidState state, FpropsError *err){
141 if(state.fluid->thcond->type != FPROPS_THCOND_1){
142 *err = FPROPS_INVALID_REQUEST;
143 return NAN;
144 }
145 const ThermalConductivityData1 *k1 = &(state.fluid->thcond->data.k1);
146
147 // value for the residual thermal conductivity
148 double lamr = 0;
149 double tau = k1->T_star / state.T;
150 double del = state.rho / k1->rho_star;
151 int i;
152 for(i=0; i < k1->nr; ++i){
153 double lamri = k1->rt[i].N * pow(tau, k1->rt[i].t) * pow(del, k1->rt[i].d);
154 if(0 == k1->rt[i].l){
155 lamr += lamri;
156 }else{
157 lamr += lamri * exp(-pow(del, k1->rt[i].l));
158 }
159 }
160 MSG("lamr(rho=%f) = %e",state.rho, lamr);
161 return lamr * k1->k_star;
162 }
163
164 /*---------------------------CRITICAL ENHANCEMENT-----------------------------*/
165 /**
166 Reduced symmetrised compressibility, as described/defined in Vesovic et al.,
167 1990, J Phys Chem Ref Data 19(3). This function is used in the calculation
168 of the critical enhancement of thermal conductivity, in particular for
169 Carbon Dioxide.
170 */
171 double thcond1_chitilde(FluidState state, FpropsError *err){
172 if(state.fluid->thcond->type != FPROPS_THCOND_1){
173 *err = FPROPS_INVALID_REQUEST;
174 return NAN;
175 }
176 double p_c = state.fluid->data->p_c;
177 double rho_c = state.fluid->data->rho_c;
178 //double T_c = state.fluid->data->T_c;
179 MSG("p_c = %e, rho_c = %f",p_c, rho_c);
180 /* FIXME we use dpdrho_T directly; assume that we have checked if we're in two-phase region or not */
181 double dpdrho_T = (*(state.fluid->dpdrho_T_fn))(state.T, state.rho, state.fluid->data, err);
182 MSG("dpdrho_T = %f",dpdrho_T);
183
184 double chitilde = p_c * state.rho / SQ(rho_c) / dpdrho_T;
185 MSG("chitilde(T=%f,rho=%f) = %f",state.T,state.rho,chitilde);
186 return chitilde;
187 }
188
189
190 double thcond1_lamc(FluidState state, FpropsError *err){
191 if(state.fluid->thcond->type != FPROPS_THCOND_1){
192 *err = FPROPS_INVALID_REQUEST;
193 return NAN;
194 }
195 //const ThermalConductivityData1 *k1 = &(state.fluid->thcond->data.k1);
196
197 /* parameters specific to CO2 */
198 double qt_D = 4.0e-10; // [m]
199 double xi0 = 1.5e-10 /* m */;
200 double Gamma = 0.052;
201 double T_ref = 450; /* K */
202
203 // 'universal' parameters, 'theoretically based parameters'
204 double R_0 = 1.01; /* see eq 35 of Vesovic et al, 1990 or eq 7 of Lemmon & Jacobsen 2004. */
205 double nu = 0.630;
206 double gamma = 1.2415;
207
208 MSG("state: T=%f, rho=%f",state.T, state.rho);
209 //const ThermalConductivityData1 *k1 = &(state.fluid->thcond->data.k1);
210
211 /* use the cp/cv functions directly, to avoid bothering with saturation boundary checks (ie assume we're outside saturation region?) */
212 double cp = (*(state.fluid->cp_fn))(state.T, state.rho, state.fluid->data, err);
213 double cv = (*(state.fluid->cv_fn))(state.T, state.rho, state.fluid->data, err);
214
215 MSG("cp = %f",cp);
216 MSG("cv = %f",cv);
217
218 #if 0
219 double T_orig = state.T;
220 if(T_orig >= 445.){
221 state.T = 445.;
222 }
223 #endif
224 FluidState state_r = state;
225 state_r.T = T_ref;
226 MSG("state_r: T=%f, rho=%f",state_r.T, state_r.rho);
227 //MSG("chitilde(state) = %e", thcond1_chitilde(state,err));
228 //MSG("chitilde(state_r) = %e", thcond1_chitilde(state_r,err));
229 //MSG("chitilde(state_r)*T_ref/T = %e", thcond1_chitilde(state_r,err)*T_ref/state.T);
230
231 double brackterm = (thcond1_chitilde(state,err) - thcond1_chitilde(state_r,err) * T_ref / state.T) / Gamma;
232
233 double lamc = 0;
234 if(brackterm<=0){
235 /* according to Lemmon & Jacobsen, we should use lamc=0 whenever brackterm <= 0 */
236 MSG("brackterm<=0 -> lamc = 0");
237 }else{
238 double xi = xi0 * pow(brackterm, nu/gamma); /* m */
239 MSG("xi = %e",xi);
240 //ASSERT(!isnan(xi));
241 #if 0
242 if(T_orig >= 445.){
243 xi *= exp(-(T_orig - 445)/10.);
244 }
245 MSG("xi = %f",xi);
246 #endif
247
248 double xioq = xi / qt_D;
249 MSG("xioq = %f",xioq);
250
251 //double xi = thcond1_xi(state, err);
252 double rho_c = state.fluid->data->rho_c;
253 double Omegatilde = 2/PI *( ((cp-cv)/cp) * atan(xioq) + (cv/cp)*xioq);
254 double Omegatilde_0 = 2/PI * (1 - exp(-1/(1./xioq + 1./3*SQ(xioq)*SQ(rho_c/state.rho))));
255 MSG("Omegatilde = %e",Omegatilde);
256 MSG("Omegatilde_0 = %e",Omegatilde_0);
257
258 double mu = visc1_mu(state, err); /* TODO ensure visc1_mu excludes crit enhancement! */
259
260 lamc = state.rho * cp * R_0 * K_BOLTZMANN *state.T/(6.*PI*xi*mu)* (Omegatilde - Omegatilde_0);
261 }
262 return lamc;
263 }
264
265 /*----------------------------------OVERALL RESULT----------------------------*/
266
267 double thcond1_lam(FluidState state, FpropsError *err){
268 // if we are here, we should be able to assume that state, should be able to remove following test (convert to assert)
269 if(state.fluid->thcond->type != FPROPS_THCOND_1){
270 *err = FPROPS_INVALID_REQUEST;
271 return NAN;
272 }
273
274 //const ThermalConductivityData1 *k1 = &(state.fluid->thcond->data.k1);
275
276 // value for the conductivity at the zero-density limit
277 double lam0 = thcond1_lam0(state,err);
278
279 double lamr = thcond1_lamr(state,err);
280
281 // FIXME need to check if k1->crit is used, etc etc
282 double lamc = thcond1_lamc(state,err);
283 MSG("lamc = %e",lamc);
284
285 return lam0 + lamr + lamc;
286 }
287
288
289
290
291
292
293
294
295
296

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