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

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