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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2768 by jpye, Thu Mar 27 06:42:18 2014 UTC revision 2769 by jpye, Thu Mar 27 08:36:59 2014 UTC
# Line 166  double thcond1_chitilde(FluidState state Line 166  double thcond1_chitilde(FluidState state
166      double dpdrho_T = (*(state.fluid->dpdrho_T_fn))(state.T, state.rho, state.fluid->data, err);      double dpdrho_T = (*(state.fluid->dpdrho_T_fn))(state.T, state.rho, state.fluid->data, err);
167      //MSG("drhodp_T = %e",1/dpdrho_T);      //MSG("drhodp_T = %e",1/dpdrho_T);
168    
169      double chitilde = (p_c / SQ(rho_c) / T_c) * state.rho * state.T * (1. / dpdrho_T);      double chitilde = p_c * state.rho / SQ(rho_c) / dpdrho_T;
170      //MSG("chitilde = %f",chitilde);      //MSG("chitilde = %f",chitilde);
171      return chitilde;      return chitilde;
172  }  }
173    
174  double thcond1_xi(FluidState state, FpropsError *err){  
175    double thcond1_lamc(FluidState state, FpropsError *err){
176      if(state.fluid->thcond->type != FPROPS_THCOND_1){      if(state.fluid->thcond->type != FPROPS_THCOND_1){
177          *err = FPROPS_INVALID_REQUEST;          *err = FPROPS_INVALID_REQUEST;
178          return NAN;          return NAN;
179      }      }
180      // parameters specific to CO2...      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 */;      double xi0 = 1.5e-10 /* m */;
185      double Gamma = 0.052;      double Gamma = 0.052;
186      double T_r = 450; /* K */      double T_ref = 450; /* K */
187    
188      // 'universal' parameter'...      // '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;      double nu = 0.630;
191      double gamma = 1.2415;      double gamma = 1.2415;
192    
193      MSG("state: T=%f, rho=%f",state.T, state.rho);      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;      double T_orig = state.T;
202      if(T_orig >= 445.){      if(T_orig >= 445.){
203          state.T = 445.;          state.T = 445.;
204      }      }
205    #endif
206      FluidState state_r = state;      FluidState state_r = state;
207      state_r.T = 445.;      state_r.T = 445.;
208      MSG("state_r: T=%f, rho=%f",state_r.T, state_r.rho);      MSG("state_r: T=%f, rho=%f",state_r.T, state_r.rho);
209      MSG("chitilde(state) = %e", thcond1_chitilde(state,err));      MSG("chitilde(state) = %e", thcond1_chitilde(state,err));
210      MSG("chitilde(state_r) = %e", thcond1_chitilde(state_r,err));      MSG("chitilde(state_r) = %e", thcond1_chitilde(state_r,err));
211      MSG("chitilde(state_r)*T_r/T = %e", thcond1_chitilde(state_r,err)*T_r/state.T);      //MSG("chitilde(state_r)*T_ref/T = %e", thcond1_chitilde(state_r,err)*T_ref/state.T);
   
     double Delta_chitilde = thcond1_chitilde(state,err) - thcond1_chitilde(state_r,err) * T_r / state.T;  
     MSG("xi0 = %e, Delta_chitilde = %f", xi0, Delta_chitilde);  
     MSG("Delta_chitilde/Gamma = %e", Delta_chitilde / Gamma);  
     double xi = xi0 * pow(Delta_chitilde / Gamma, nu/gamma); /* m */  
     ASSERT(!isnan(xi));  
     if(T_orig >= 445.){  
         xi *= exp(-(T_orig - 445)/10.);  
     }  
     MSG("xi = %f",xi);  
     return xi;  
212            
213  }      double brackterm = (thcond1_chitilde(state,err) - thcond1_chitilde(state_r,err) * T_ref / state.T) / Gamma;
214    
215  double thcond1_lamc(FluidState state, FpropsError *err){      double lamc = 0;
216      if(state.fluid->thcond->type != FPROPS_THCOND_1){      if(brackterm<=0){
217          *err = FPROPS_INVALID_REQUEST;          /* according to Lemmon & Jacobsen, we should use lamc=0 whenever brackterm <= 0 */
218          return NAN;          MSG("brackterm<=0 -> lamc = 0");
219      }      }else{
220      //const ThermalConductivityData1 *k1 = &(state.fluid->thcond->data.k1);          double xi = xi0 * pow(brackterm, nu/gamma); /* m */
221            ASSERT(!isnan(xi));
222      /* use the cp/cv functions directly, to avoid bothering with saturation boundary checks (ie assume we're outside saturation region?) */  #if 0
223      double cp = (*(state.fluid->cp_fn))(state.T, state.rho, state.fluid->data, err);          if(T_orig >= 445.){
224      double cv = (*(state.fluid->cp_fn))(state.T, state.rho, state.fluid->data, err);              xi *= exp(-(T_orig - 445)/10.);
225      double xi = thcond1_xi(state, err);          }
226      double rho_c = state.fluid->data->rho_c;          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      /* this value should be stored in k1 */          double mu = visc1_mu(state, err); /* TODO ensure visc1_mu excludes crit enhancement! */
     double qtilde_D = 1./(4.0e-10); // [m^-1]  
237    
238      double Omegatilde = 2/PI * ((cp-cv)/cp * atan(qtilde_D * xi) + cv/cp*qtilde_D*xi);          lamc = k1->k_star * state.rho * cp * R_0 * K_BOLTZMANN *state.T/(6.*PI*xi*mu)* (Omegatilde - Omegatilde_0);
239      double Omegatilde_0 = 2/PI * (1 - exp(-1/(1./(qtilde_D*xi) + 1./3*SQ(qtilde_D*xi*rho_c/state.rho))));      }
   
     double R = 1.01; /* 'universal' amplitude parameter, see eq 35 of Vesovic et al, 1990 */  
     double mubar = visc1_mu(state, err); /* TODO make sure visc1_mu excludes critical enhancement */  
   
     double lamc = state.rho * cp * R * K_BOLTZMANN *state.T/(6*PI*mubar*xi)* (Omegatilde - Omegatilde_0);  
240      return lamc;      return lamc;
241  }  }
242    

Legend:
Removed from v.2768  
changed lines
  Added in v.2769

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