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

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