1 //-----------------------------------------------------------------------bl- 2 //-------------------------------------------------------------------------- 3 // 4 // Antioch - A Gas Dynamics Thermochemistry Library 5 // 6 // Copyright (C) 2014-2016 Paul T. Bauman, Benjamin S. Kirk, 7 // Sylvain Plessis, Roy H. Stonger 8 // 9 // Copyright (C) 2013 The PECOS Development Team 10 // 11 // This library is free software; you can redistribute it and/or 12 // modify it under the terms of the Version 2.1 GNU Lesser General 13 // Public License as published by the Free Software Foundation. 14 // 15 // This library is distributed in the hope that it will be useful, 16 // but WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 // Lesser General Public License for more details. 19 // 20 // You should have received a copy of the GNU Lesser General Public 21 // License along with this library; if not, write to the Free Software 22 // Foundation, Inc. 51 Franklin Street, Fifth Floor, 23 // Boston, MA 02110-1301 USA 24 // 25 //-----------------------------------------------------------------------el- 26 27 28 #ifndef ANTIOCH_STAT_MECH_THERMO_H 29 #define ANTIOCH_STAT_MECH_THERMO_H 30 31 // Antioch 32 #include "antioch/macro_micro_thermo_base.h" 33 #include "antioch/antioch_exceptions.h" 34 35 // C++ 36 #include <vector> 37 #include <cmath> 38 #include <limits> 39 40 namespace Antioch 41 { 42 43 template<typename CoeffType=double> 44 class StatMechThermodynamics : public MacroMicroThermoBase<CoeffType,StatMechThermodynamics<CoeffType> > 45 { 46 public: 47 StatMechThermodynamics(const ChemicalMixture<CoeffType> & chem_mixture)48 StatMechThermodynamics( const ChemicalMixture<CoeffType>& chem_mixture ) 49 : MacroMicroThermoBase<CoeffType,StatMechThermodynamics<CoeffType> >(chem_mixture) 50 {} 51 ~StatMechThermodynamics()52 virtual ~StatMechThermodynamics(){} 53 54 /** 55 * @returns species vibrational/electronic specific heat 56 * constant volume. 57 */ 58 template<typename StateType> 59 StateType cv_ve (const unsigned int species, const StateType& Tv) const; 60 61 /** 62 * @returns mixture vibrational/electronic specific heat 63 * constant volume. 64 */ 65 template<typename VectorStateType> 66 typename enable_if_c< 67 has_size<VectorStateType>::value, 68 typename Antioch::value_type<VectorStateType>::type 69 >::type 70 cv_ve (const typename Antioch::value_type<VectorStateType>::type& Tv, 71 const VectorStateType& mass_fractions) const; 72 73 /** 74 * @returns species total specific heat at constant volume, J/kg. Note that the input 75 * translational/rotational temperature is currently not used, as these 76 * modes are assumed to be fully excited. However, the API is here 77 * in case this assumption is removed later. 78 */ 79 template<typename StateType> 80 StateType cv (const unsigned int species, const StateType& T, const StateType& Tv) const; 81 82 /** 83 * @returns mixture total specific heat at constant volume, J/kg. Note that the input 84 * translational/rotational temperature is currently not used, as these 85 * modes are assumed to be fully excited. However, the API is here 86 * in case this assumption is removed later. 87 */ 88 template<typename VectorStateType> 89 typename enable_if_c< 90 has_size<VectorStateType>::value, 91 typename Antioch::value_type<VectorStateType>::type 92 >::type 93 cv (const typename Antioch::value_type<VectorStateType>::type& T, 94 const typename Antioch::value_type<VectorStateType>::type& Tv, 95 const VectorStateType& mass_fractions) const; 96 97 /** 98 * @returns mixture total specific heat at constant pressure, J/kg. 99 */ 100 template<typename VectorStateType> 101 typename enable_if_c< 102 has_size<VectorStateType>::value, 103 typename Antioch::value_type<VectorStateType>::type 104 >::type 105 cp (const typename Antioch::value_type<VectorStateType>::type& T, 106 const typename Antioch::value_type<VectorStateType>::type& Tv, 107 const VectorStateType& mass_fractions) const; 108 109 /** 110 * @returns species total enthalpy, J/kg. 111 */ 112 template<typename StateType> 113 StateType h_tot (const unsigned int species, const StateType& T, const StateType& Tv) const; 114 115 /** 116 * @returns species total enthalpy, J/kg. Assumes thermal equilibrium 117 * of translational/rotational and vibrational/electronic temperature. 118 */ 119 template<typename StateType> 120 StateType h_tot (const unsigned int species, const StateType& T) const; 121 122 /** 123 * @returns mixture specific enthalpy, J/kg. 124 */ 125 template<typename VectorStateType> 126 typename enable_if_c< 127 has_size<VectorStateType>::value, 128 typename Antioch::value_type<VectorStateType>::type 129 >::type 130 h_tot (const typename Antioch::value_type<VectorStateType>::type& T, 131 const typename Antioch::value_type<VectorStateType>::type& Tv, 132 const VectorStateType& mass_fractions) const; 133 134 /** 135 * @returns mixture specific enthalpy, J/kg. Assumes thermal equilibrium 136 * of translational/rotational and vivbrational/electronic temperature. 137 */ 138 template<typename VectorStateType> 139 typename enable_if_c< 140 has_size<VectorStateType>::value, 141 typename Antioch::value_type<VectorStateType>::type 142 >::type 143 h_tot (const typename Antioch::value_type<VectorStateType>::type& T, 144 const VectorStateType& mass_fractions) const; 145 146 /** 147 * @returns species total internal energy, J/kg. 148 */ 149 template<typename StateType> 150 StateType e_tot (const unsigned int species, const StateType& T, const StateType& Tv) const; 151 152 /** 153 * @returns mixture total internal energy, J/kg. 154 */ 155 template<typename VectorStateType> 156 typename enable_if_c< 157 has_size<VectorStateType>::value, 158 typename Antioch::value_type<VectorStateType>::type 159 >::type 160 e_tot (const typename Antioch::value_type<VectorStateType>::type& T, 161 const typename Antioch::value_type<VectorStateType>::type& Tv, 162 const VectorStateType& mass_fractions) const; 163 164 /** 165 * @returns species total internal energy, J/kg. Assumes thermal equilibrium 166 * of translational/rotational and vibrational/electronic temperature. 167 */ 168 template<typename StateType> 169 StateType e_tot (const unsigned int species, const StateType& T) const; 170 171 /** 172 * @returns mixture total internal energy, J/kg. Assumes thermal equilibrium 173 * of translational/rotational and vivbrational/electronic temperature. 174 */ 175 template<typename VectorStateType> 176 typename enable_if_c< 177 has_size<VectorStateType>::value, 178 typename Antioch::value_type<VectorStateType>::type 179 >::type 180 e_tot (const typename Antioch::value_type<VectorStateType>::type& T, 181 const VectorStateType& mass_fractions) const; 182 183 /** 184 * @returns species translational/rotational energy, J/kg. 185 */ 186 template<typename StateType> 187 StateType e_tr (const unsigned int species, const StateType& T) const; 188 189 /** 190 * @returns mixture translational/rotational energy, J/kg. 191 */ 192 template<typename VectorStateType> 193 typename enable_if_c< 194 has_size<VectorStateType>::value, 195 typename Antioch::value_type<VectorStateType>::type 196 >::type 197 e_tr (const typename Antioch::value_type<VectorStateType>::type& T, 198 const VectorStateType& mass_fractions) const; 199 200 /** 201 * @returns species vibrational energy, J/kg. 202 */ 203 template<typename StateType> 204 StateType e_vib (const unsigned int species, const StateType& Tv) const; 205 206 /** 207 * @returns mixture vibrational energy, J/kg. 208 */ 209 template<typename VectorStateType> 210 typename enable_if_c< 211 has_size<VectorStateType>::value, 212 typename Antioch::value_type<VectorStateType>::type 213 >::type 214 e_vib (const typename Antioch::value_type<VectorStateType>::type& Tv, 215 const VectorStateType& mass_fractions) const; 216 217 /** 218 * @returns species electronic energy, J/kg. 219 */ 220 template<typename StateType> 221 StateType e_el (const unsigned int species, const StateType& Te) const; 222 223 /** 224 * @returns mixture electronic energy, J/kg. 225 */ 226 template<typename VectorStateType> 227 typename enable_if_c< 228 has_size<VectorStateType>::value, 229 typename Antioch::value_type<VectorStateType>::type 230 >::type 231 e_el (const typename Antioch::value_type<VectorStateType>::type& Te, 232 const VectorStateType& mass_fractions) const; 233 234 /** 235 * @returns species vibrational/electronic energy, J/kg. 236 */ 237 template<typename StateType> 238 StateType e_ve (const unsigned int species, const StateType& Tv) const; 239 240 /** 241 * @returns mixture vibrational/electronic energy, J/kg. 242 */ 243 template<typename VectorStateType> 244 typename enable_if_c< 245 has_size<VectorStateType>::value, 246 typename Antioch::value_type<VectorStateType>::type 247 >::type 248 e_ve (const typename Antioch::value_type<VectorStateType>::type& Te, 249 const VectorStateType& mass_fractions) const; 250 251 /** 252 * @returns the minimum valid mixture vibrational/electronic energy, J/kg. 253 * We allow for a user-specified minimum vibrational temperature, Tv. 254 * This in turn sets an a priori minimum valid species 255 * vibrational/electronic temperature. However, the minimum valid 256 * value for a mixture depends on the species mass fractions and must 257 * be computed. That is what this method does. 258 */ 259 CoeffType e_ve_min () const; 260 261 /** 262 * @returns mixture vibrational/electronic energy (same as 263 * e_ve()) and specific heat (same as cv_ve()), calculated 264 * together for efficiency. 265 */ 266 template<typename VectorStateType> 267 void e_and_cv_ve (const typename Antioch::value_type<VectorStateType>::type& Tv, 268 const VectorStateType& mass_fractions, 269 typename Antioch::value_type<VectorStateType>::type &e_ve, 270 typename Antioch::value_type<VectorStateType>::type &cv_ve) const; 271 272 /** 273 * @returns species formation energy, J/kg. 274 */ 275 CoeffType e_0 (const unsigned int species) const; 276 277 /** 278 * @returns mixture formation energy, J/kg. 279 */ 280 template<typename VectorStateType> 281 typename enable_if_c< 282 has_size<VectorStateType>::value, 283 typename Antioch::value_type<VectorStateType>::type 284 >::type 285 e_0 (const VectorStateType& mass_fractions) const; 286 287 /** 288 * Computes the mixture vibrational temperature (K) from input 289 * vibrational-electronic energy per unit mass (J/kg). 290 */ 291 template<typename VectorStateType> 292 typename enable_if_c< 293 has_size<VectorStateType>::value, 294 typename Antioch::value_type<VectorStateType>::type 295 >::type 296 Tv_from_e_ve (const typename Antioch::value_type<VectorStateType>::type& e_ve, 297 const VectorStateType& mass_fractions, 298 typename Antioch::value_type<VectorStateType>::type Tv = -1) const; 299 300 /** 301 * Computes the mixture temperature (K) from input 302 * total energy per unit mass (J/kg). 303 */ 304 template<typename VectorStateType> 305 typename enable_if_c< 306 has_size<VectorStateType>::value, 307 typename Antioch::value_type<VectorStateType>::type 308 >::type 309 T_from_e_tot (const typename Antioch::value_type<VectorStateType>::type& e_tot, 310 const VectorStateType& mass_fractions, 311 typename Antioch::value_type<VectorStateType>::type T = -1) const; 312 313 /** 314 * Computes the mixture temperature (K) from input 315 * translational/rotational energy per unit mass (J/kg). 316 */ 317 template<typename VectorStateType> 318 typename enable_if_c< 319 has_size<VectorStateType>::value, 320 typename Antioch::value_type<VectorStateType>::type 321 >::type 322 T_from_e_tr (const typename Antioch::value_type<VectorStateType>::type& e_tr, 323 const VectorStateType& mass_fractions, 324 typename Antioch::value_type<VectorStateType>::type T = -1) const; 325 326 /** 327 * Computes the mixture temperature (K) from input 328 * enthalpy per unit mass (J/kg). 329 */ 330 template<typename VectorStateType> 331 typename enable_if_c< 332 has_size<VectorStateType>::value, 333 typename Antioch::value_type<VectorStateType>::type 334 >::type 335 T_from_h_tot (const typename Antioch::value_type<VectorStateType>::type& h_tot, 336 const VectorStateType& mass_fractions, 337 typename Antioch::value_type<VectorStateType>::type T = -1) const; 338 339 /** 340 * Computes the mixture temperature (K) from input 341 * vibrational/electronic temperature Tv (K) and 342 * enthalpy per unit mass (J/kg). 343 */ 344 template<typename VectorStateType> 345 typename enable_if_c< 346 has_size<VectorStateType>::value, 347 typename Antioch::value_type<VectorStateType>::type 348 >::type 349 T_from_h_tot_Tv (const typename Antioch::value_type<VectorStateType>::type& h_tot, 350 const typename Antioch::value_type<VectorStateType>::type& Tv, 351 const VectorStateType& mass_fractions, 352 typename Antioch::value_type<VectorStateType>::type T = -1) const; 353 354 /** 355 * Computes the mixture specific entropy (J/kg-K) from input 356 * temperature (K) and pressure (Pa). 357 */ 358 template<typename VectorStateType> 359 typename enable_if_c< 360 has_size<VectorStateType>::value, 361 typename Antioch::value_type<VectorStateType>::type 362 >::type 363 s (const typename Antioch::value_type<VectorStateType>::type& T, 364 const typename Antioch::value_type<VectorStateType>::type& p, 365 const VectorStateType& mass_fractions) const; 366 367 // Friend the base class so we can make the CRTP implementation functions 368 // private. 369 friend class MacroMicroThermoBase<CoeffType,StatMechThermodynamics<CoeffType> >; 370 371 private: 372 373 //! Implemenation of species vibrational specific heat, [J/kg-K] 374 template<typename StateType> 375 StateType cv_vib_impl (const unsigned int species, const StateType & T) const; 376 377 //! Implemenation of normalized species vibrational specific heat. 378 template<typename StateType> 379 StateType cv_vib_over_R_impl (const unsigned int species, const StateType & T) const; 380 381 //! Implemenation of species electronic specific heat, [J/kg-K] 382 template<typename StateType> 383 StateType cv_el_impl (const unsigned int species, const StateType & T) const; 384 385 //! Default constructor 386 /*! Private to force to user to supply a ChemicalMixture object.*/ 387 StatMechThermodynamics(); 388 }; 389 390 391 /* ------------------------- Inline Functions -------------------------*/ 392 template<typename CoeffType> 393 template<typename StateType> 394 inline cv_vib_impl(const unsigned int species,const StateType & T)395 StateType StatMechThermodynamics<CoeffType>::cv_vib_impl (const unsigned int species, 396 const StateType& T) const 397 { 398 return this->cv_vib_over_R(species,T) * (this->_chem_mixture.chemical_species()[species])->gas_constant(); 399 } 400 401 template<typename CoeffType> 402 template<typename StateType> 403 inline cv_vib_over_R_impl(const unsigned int species,const StateType & T)404 StateType StatMechThermodynamics<CoeffType>::cv_vib_over_R_impl (const unsigned int species, 405 const StateType& T) const 406 { 407 using std::exp; 408 409 // convenience 410 const ChemicalSpecies<CoeffType>& chem_species = *(this->_chem_mixture.chemical_species()[species]); 411 const std::vector<CoeffType>& theta_v = chem_species.theta_v(); 412 const std::vector<unsigned int>& ndg_v = chem_species.ndg_v(); 413 414 antioch_assert_equal_to(ndg_v.size(), theta_v.size()); 415 416 // Use an input datum to make sure we get the size right 417 StateType cv_vib_over_R = Antioch::zero_clone(T); 418 419 if (theta_v.empty()) 420 return cv_vib_over_R; 421 422 for (unsigned int level=0; level<ndg_v.size(); level++) 423 { 424 typedef typename Antioch::raw_value_type<StateType>::type raw_type; 425 426 const StateType expval = exp(theta_v[level]/T); 427 428 const StateType expvalminusone = expval - raw_type(1); 429 430 cv_vib_over_R += (static_cast<CoeffType>(ndg_v[level])* 431 theta_v[level]*theta_v[level]*expval/(expvalminusone*expvalminusone)/(T*T)); 432 } 433 434 return cv_vib_over_R; 435 } 436 437 438 439 template<typename CoeffType> 440 template<typename StateType> 441 inline cv_el_impl(const unsigned int species,const StateType & T)442 StateType StatMechThermodynamics<CoeffType>::cv_el_impl (const unsigned int species, 443 const StateType& T) const 444 { 445 using std::exp; 446 447 // convenience 448 const ChemicalSpecies<CoeffType>& chem_species = *(this->_chem_mixture.chemical_species()[species]); 449 const std::vector<CoeffType>& theta_e = chem_species.theta_e(); 450 const std::vector<unsigned int>& ndg_e = chem_species.ndg_e(); 451 452 antioch_assert_equal_to(ndg_e.size(), theta_e.size()); 453 454 StateType cv_el = Antioch::zero_clone(T); 455 456 // Really < 2? Yes, b/c theta_e[0] = 0.0 always. See 457 // antioch_default_electronic_data.dat 458 if (theta_e.size() < 2) return cv_el; 459 460 typedef typename Antioch::raw_value_type<StateType>::type raw_type; 461 const raw_type one = static_cast<raw_type>(1); 462 463 const StateType Teinv = one/T; 464 const StateType Te2inv = Teinv*Teinv; 465 466 StateType 467 num = Antioch::zero_clone(T), dnum = Antioch::zero_clone(T), 468 den = Antioch::zero_clone(T), dden = Antioch::zero_clone(T); 469 470 for (unsigned int level=0; level<theta_e.size(); level++) 471 { 472 const StateType 473 expval = exp (-theta_e[level] * Teinv), 474 den_l = static_cast<raw_type>(ndg_e[level])*expval, 475 num_l = den_l*theta_e[level], 476 dden_l = num_l*Te2inv, 477 dnum_l = dden_l*theta_e[level]; 478 479 num += num_l; 480 den += den_l; 481 482 dden += dden_l; 483 dnum += dnum_l; 484 } 485 486 const StateType invden = one/den; 487 488 cv_el = chem_species.gas_constant()*(dnum - num*dden*invden) * invden; 489 490 return cv_el; 491 } 492 493 template<typename CoeffType> 494 template<typename StateType> 495 inline cv_ve(const unsigned int species,const StateType & Tv)496 StateType StatMechThermodynamics<CoeffType>::cv_ve (const unsigned int species, 497 const StateType& Tv) const 498 { 499 return (this->cv_vib(species, Tv) + this->cv_el(species, Tv)); 500 } 501 502 template<typename CoeffType> 503 template<typename VectorStateType> 504 inline 505 typename enable_if_c< 506 has_size<VectorStateType>::value, 507 typename Antioch::value_type<VectorStateType>::type 508 >::type cv_ve(const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions)509 StatMechThermodynamics<CoeffType>::cv_ve (const typename Antioch::value_type<VectorStateType>::type& Tv, 510 const VectorStateType& mass_fractions) const 511 512 { 513 return (this->cv_vib(Tv, mass_fractions) + this->cv_el(Tv, mass_fractions)); 514 } 515 516 template<typename CoeffType> 517 template<typename StateType> 518 inline cv(const unsigned int species,const StateType &,const StateType & Tv)519 StateType StatMechThermodynamics<CoeffType>::cv (const unsigned int species, 520 const StateType& /* T */, 521 const StateType& Tv) const 522 { 523 return (this->cv_tr(species) + this->cv_ve(species, Tv)); 524 } 525 526 template<typename CoeffType> 527 template<typename VectorStateType> 528 inline 529 typename enable_if_c< 530 has_size<VectorStateType>::value, 531 typename Antioch::value_type<VectorStateType>::type 532 >::type cv(const typename Antioch::value_type<VectorStateType>::type &,const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions)533 StatMechThermodynamics<CoeffType>::cv (const typename Antioch::value_type<VectorStateType>::type& /* T */, 534 const typename Antioch::value_type<VectorStateType>::type& Tv, 535 const VectorStateType& mass_fractions) const 536 { 537 return (this->cv_tr(mass_fractions) + this->cv_ve(Tv, mass_fractions)); 538 } 539 540 template<typename CoeffType> 541 template<typename VectorStateType> 542 inline 543 typename enable_if_c< 544 has_size<VectorStateType>::value, 545 typename Antioch::value_type<VectorStateType>::type 546 >::type cp(const typename Antioch::value_type<VectorStateType>::type & T,const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions)547 StatMechThermodynamics<CoeffType>::cp (const typename Antioch::value_type<VectorStateType>::type& T, 548 const typename Antioch::value_type<VectorStateType>::type& Tv, 549 const VectorStateType& mass_fractions) const 550 { 551 return (this->cv(T,Tv,mass_fractions) + this->_chem_mixture.R(mass_fractions)); 552 } 553 554 template<typename CoeffType> 555 template<typename StateType> 556 inline h_tot(const unsigned int species,const StateType & T,const StateType & Tv)557 StateType StatMechThermodynamics<CoeffType>::h_tot (const unsigned int species, 558 const StateType& T, 559 const StateType& Tv) const 560 { 561 const ChemicalSpecies<CoeffType>& chem_species = *(this->_chem_mixture.chemical_species()[species]); 562 return (this->e_tot(species, T, Tv) + chem_species.gas_constant()*T); 563 } 564 565 566 template<typename CoeffType> 567 template<typename StateType> 568 inline h_tot(const unsigned int species,const StateType & T)569 StateType StatMechThermodynamics<CoeffType>::h_tot (const unsigned int species, 570 const StateType& T) const 571 { 572 return this->h_tot(species, T, T); 573 } 574 575 template<typename CoeffType> 576 template<typename VectorStateType> 577 inline 578 typename enable_if_c< 579 has_size<VectorStateType>::value, 580 typename Antioch::value_type<VectorStateType>::type 581 >::type h_tot(const typename Antioch::value_type<VectorStateType>::type & T,const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions)582 StatMechThermodynamics<CoeffType>::h_tot (const typename Antioch::value_type<VectorStateType>::type& T, 583 const typename Antioch::value_type<VectorStateType>::type& Tv, 584 const VectorStateType& mass_fractions) const 585 { 586 typename Antioch::value_type<VectorStateType>::type 587 h_tot = mass_fractions[0]*this->h_tot(0, T, Tv); 588 589 for( unsigned int s = 1; s < this->_chem_mixture.n_species(); s++ ) 590 { 591 h_tot += mass_fractions[s]*this->h_tot(s, T, Tv); 592 } 593 594 return h_tot; 595 } 596 597 template<typename CoeffType> 598 template<typename VectorStateType> 599 inline 600 typename enable_if_c< 601 has_size<VectorStateType>::value, 602 typename Antioch::value_type<VectorStateType>::type 603 >::type h_tot(const typename Antioch::value_type<VectorStateType>::type & T,const VectorStateType & mass_fractions)604 StatMechThermodynamics<CoeffType>::h_tot (const typename Antioch::value_type<VectorStateType>::type& T, 605 const VectorStateType& mass_fractions) const 606 { 607 return this->h_tot(T, T, mass_fractions); 608 } 609 610 template<typename CoeffType> 611 template<typename StateType> 612 inline e_tot(const unsigned int species,const StateType & T,const StateType & Tv)613 StateType StatMechThermodynamics<CoeffType>::e_tot (const unsigned int species, 614 const StateType& T, 615 const StateType& Tv) const 616 { 617 return (this->e_tr(species, T) + this->e_ve(species, Tv) + this->e_0(species)); 618 } 619 620 621 template<typename CoeffType> 622 template<typename VectorStateType> 623 inline 624 typename enable_if_c< 625 has_size<VectorStateType>::value, 626 typename Antioch::value_type<VectorStateType>::type 627 >::type e_tot(const typename Antioch::value_type<VectorStateType>::type & T,const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions)628 StatMechThermodynamics<CoeffType>::e_tot (const typename Antioch::value_type<VectorStateType>::type& T, 629 const typename Antioch::value_type<VectorStateType>::type& Tv, 630 const VectorStateType& mass_fractions) const 631 { 632 return (this->e_tr(T, mass_fractions) + this->e_ve(Tv, mass_fractions) + this->e_0(mass_fractions)); 633 } 634 635 template<typename CoeffType> 636 template<typename StateType> 637 inline e_tot(const unsigned int species,const StateType & T)638 StateType StatMechThermodynamics<CoeffType>::e_tot (const unsigned int species, 639 const StateType& T) const 640 { 641 return this->e_tot(species, T, T); 642 } 643 644 template<typename CoeffType> 645 template<typename VectorStateType> 646 inline 647 typename enable_if_c< 648 has_size<VectorStateType>::value, 649 typename Antioch::value_type<VectorStateType>::type 650 >::type e_tot(const typename Antioch::value_type<VectorStateType>::type & T,const VectorStateType & mass_fractions)651 StatMechThermodynamics<CoeffType>::e_tot (const typename Antioch::value_type<VectorStateType>::type& T, 652 const VectorStateType& mass_fractions) const 653 { 654 return this->e_tot(T, T, mass_fractions); 655 } 656 657 template<typename CoeffType> 658 template<typename StateType> 659 inline e_tr(const unsigned int species,const StateType & T)660 StateType StatMechThermodynamics<CoeffType>::e_tr (const unsigned int species, 661 const StateType& T) const 662 { 663 return this->cv_tr(species)*T; 664 } 665 666 template<typename CoeffType> 667 template<typename VectorStateType> 668 inline 669 typename enable_if_c< 670 has_size<VectorStateType>::value, 671 typename Antioch::value_type<VectorStateType>::type 672 >::type e_tr(const typename Antioch::value_type<VectorStateType>::type & T,const VectorStateType & mass_fractions)673 StatMechThermodynamics<CoeffType>::e_tr (const typename Antioch::value_type<VectorStateType>::type& T, 674 const VectorStateType& mass_fractions) const 675 { 676 typename Antioch::value_type<VectorStateType>::type 677 e_tr = mass_fractions[0]*this->e_tr(0, T); 678 679 for( unsigned int s = 1; s < this->_chem_mixture.n_species(); s++ ) 680 { 681 e_tr += mass_fractions[s]*this->e_tr(s, T); 682 } 683 684 return e_tr; 685 } 686 687 template<typename CoeffType> 688 template<typename StateType> 689 inline e_vib(const unsigned int species,const StateType & Tv)690 StateType StatMechThermodynamics<CoeffType>::e_vib (const unsigned int species, 691 const StateType& Tv) const 692 { 693 using std::exp; 694 695 // convenience 696 const ChemicalSpecies<CoeffType>& chem_species = *(this->_chem_mixture.chemical_species()[species]); 697 const std::vector<CoeffType>& theta_v = chem_species.theta_v(); 698 const std::vector<unsigned int>& ndg_v = chem_species.ndg_v(); 699 700 antioch_assert_equal_to(ndg_v.size(), theta_v.size()); 701 702 StateType e_vib = 0.0; 703 704 if (theta_v.empty()) return e_vib; 705 706 for (unsigned int level=0; level<ndg_v.size(); level++) 707 e_vib += ndg_v[level]*chem_species.gas_constant()*theta_v[level]/(exp(theta_v[level]/Tv) - 1.); 708 709 return e_vib; 710 } 711 712 template<typename CoeffType> 713 template<typename VectorStateType> 714 inline 715 typename enable_if_c< 716 has_size<VectorStateType>::value, 717 typename Antioch::value_type<VectorStateType>::type 718 >::type e_vib(const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions)719 StatMechThermodynamics<CoeffType>::e_vib (const typename Antioch::value_type<VectorStateType>::type& Tv, 720 const VectorStateType& mass_fractions) const 721 { 722 typename Antioch::value_type<VectorStateType>::type 723 e_vib = mass_fractions[0]*this->e_vib(0, Tv); 724 725 for( unsigned int s = 1; s < this->_chem_mixture.n_species(); s++ ) 726 { 727 e_vib += mass_fractions[s]*this->e_vib(s, Tv); 728 } 729 730 return e_vib; 731 } 732 733 template<typename CoeffType> 734 template<typename StateType> 735 inline e_el(const unsigned int species,const StateType & Te)736 StateType StatMechThermodynamics<CoeffType>::e_el (const unsigned int species, 737 const StateType& Te) const 738 { 739 using std::exp; 740 741 // convenience 742 const ChemicalSpecies<CoeffType>& chem_species = *(this->_chem_mixture.chemical_species()[species]); 743 const std::vector<CoeffType>& theta_e = chem_species.theta_e(); 744 const std::vector<unsigned int>& ndg_e = chem_species.ndg_e(); 745 746 antioch_assert_equal_to(ndg_e.size(), theta_e.size()); 747 748 StateType e_el = Antioch::zero_clone(Te); 749 750 if (theta_e.size() < 2) return e_el; 751 752 StateType num = Antioch::zero_clone(Te), den = Antioch::zero_clone(Te); 753 754 for (unsigned int level=0; level<theta_e.size(); level++) 755 { 756 const StateType expval = exp (-theta_e[level] / Te); 757 num += static_cast<StateType>(ndg_e[level])*theta_e[level]*expval; 758 den += static_cast<StateType>(ndg_e[level])*expval; 759 } 760 761 return chem_species.gas_constant() * num / den; 762 } 763 764 template<typename CoeffType> 765 template<typename VectorStateType> 766 inline 767 typename enable_if_c< 768 has_size<VectorStateType>::value, 769 typename Antioch::value_type<VectorStateType>::type 770 >::type e_el(const typename Antioch::value_type<VectorStateType>::type & Te,const VectorStateType & mass_fractions)771 StatMechThermodynamics<CoeffType>::e_el (const typename Antioch::value_type<VectorStateType>::type& Te, 772 const VectorStateType& mass_fractions) const 773 { 774 typename Antioch::value_type<VectorStateType>::type 775 e_el = mass_fractions[0]*this->e_el(0, Te); 776 777 for( unsigned int s = 1; s < this->_chem_mixture.n_species(); s++ ) 778 { 779 e_el += mass_fractions[s]*this->e_el(s, Te); 780 } 781 782 return e_el; 783 } 784 785 template<typename CoeffType> 786 template<typename StateType> 787 inline e_ve(const unsigned int species,const StateType & Tv)788 StateType StatMechThermodynamics<CoeffType>::e_ve (const unsigned int species, 789 const StateType& Tv) const 790 { 791 return (this->e_vib(species,Tv) + this->e_el(species,Tv)); 792 } 793 794 template<typename CoeffType> 795 template<typename VectorStateType> 796 inline 797 typename enable_if_c< 798 has_size<VectorStateType>::value, 799 typename Antioch::value_type<VectorStateType>::type 800 >::type e_ve(const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions)801 StatMechThermodynamics<CoeffType>::e_ve (const typename Antioch::value_type<VectorStateType>::type& Tv, 802 const VectorStateType& mass_fractions) const 803 { 804 return (this->e_vib(Tv,mass_fractions) + this->e_el(Tv,mass_fractions)); 805 } 806 807 template<typename CoeffType> 808 inline e_ve_min()809 CoeffType StatMechThermodynamics<CoeffType>::e_ve_min () const 810 { 811 antioch_not_implemented(); 812 } 813 814 template<typename CoeffType> 815 template<typename VectorStateType> 816 inline e_and_cv_ve(const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions,typename Antioch::value_type<VectorStateType>::type & e_ve,typename Antioch::value_type<VectorStateType>::type & cv_ve)817 void StatMechThermodynamics<CoeffType>::e_and_cv_ve (const typename Antioch::value_type<VectorStateType>::type& Tv, 818 const VectorStateType& mass_fractions, 819 typename Antioch::value_type<VectorStateType>::type &e_ve, 820 typename Antioch::value_type<VectorStateType>::type &cv_ve) const 821 { 822 e_ve = this->e_ve (Tv, mass_fractions); 823 cv_ve = this->cv_ve(Tv, mass_fractions); 824 } 825 826 template<typename CoeffType> 827 inline e_0(const unsigned int species)828 CoeffType StatMechThermodynamics<CoeffType>::e_0 (const unsigned int species) const 829 { 830 const ChemicalSpecies<CoeffType>& chem_species = *(this->_chem_mixture.chemical_species()[species]); 831 return chem_species.formation_enthalpy(); 832 } 833 834 template<typename CoeffType> 835 template<typename VectorStateType> 836 inline 837 typename enable_if_c< 838 has_size<VectorStateType>::value, 839 typename Antioch::value_type<VectorStateType>::type 840 >::type e_0(const VectorStateType & mass_fractions)841 StatMechThermodynamics<CoeffType>::e_0 (const VectorStateType& mass_fractions) const 842 { 843 typename Antioch::value_type<VectorStateType>::type 844 e_0 = mass_fractions[0]*this->e_0(0); 845 846 for( unsigned int s = 1; s < this->_chem_mixture.n_species(); s++ ) 847 { 848 e_0 += mass_fractions[s]*this->e_0(s); 849 } 850 851 return e_0; 852 } 853 854 template<typename CoeffType> 855 template<typename VectorStateType> 856 inline 857 typename enable_if_c< 858 has_size<VectorStateType>::value, 859 typename Antioch::value_type<VectorStateType>::type 860 >::type Tv_from_e_ve(const typename Antioch::value_type<VectorStateType>::type & e_ve,const VectorStateType & mass_fractions,typename Antioch::value_type<VectorStateType>::type Tv)861 StatMechThermodynamics<CoeffType>::Tv_from_e_ve 862 (const typename Antioch::value_type<VectorStateType>::type& e_ve, 863 const VectorStateType& mass_fractions, 864 typename Antioch::value_type<VectorStateType>::type Tv) const 865 { 866 antioch_not_implemented(); 867 } 868 869 template<typename CoeffType> 870 template<typename VectorStateType> 871 inline 872 typename enable_if_c< 873 has_size<VectorStateType>::value, 874 typename Antioch::value_type<VectorStateType>::type 875 >::type T_from_e_tot(const typename Antioch::value_type<VectorStateType>::type & e_tot,const VectorStateType & mass_fractions,typename Antioch::value_type<VectorStateType>::type T)876 StatMechThermodynamics<CoeffType>::T_from_e_tot (const typename Antioch::value_type<VectorStateType>::type& e_tot, 877 const VectorStateType& mass_fractions, 878 typename Antioch::value_type<VectorStateType>::type T) const 879 { 880 using std::abs; 881 using std::max; 882 using std::min; 883 884 typedef typename Antioch::value_type<VectorStateType>::type StateType; 885 886 // Cache the translational/rotational specific heat - this will be used repeatedly 887 // and involves (2*NS-1) flops to compute, and since this has no functional 888 // dependence on temperature it will not change throughout the Newton iteration. 889 const typename Antioch::value_type<VectorStateType>::type 890 Cv_tr = this->cv_tr(mass_fractions); 891 892 // Similarly for mixture formation energy 893 const typename Antioch::value_type<VectorStateType>::type 894 E_0 = this->e_0(mass_fractions); 895 896 // if the user does not provide an initial guess for the temperature 897 // assume it is all in translation/rotation to compute a starting value. 898 if (T < 0) 899 { 900 T = (e_tot - E_0) / Cv_tr; 901 T = min(max(T,StateType(10.)),StateType(20000.)); 902 903 // FIXME: Use Antioch::Limits or similar? (i.e., don't 904 // hardcode min and max T) 905 906 // make sure the initial guess is valid 907 //T = max(T, Limits::CompNSLimits::T_min()); 908 T = max(T, StateType(10.)); 909 T = min(T, StateType(2.e4)); 910 } 911 912 // compute the translational/rotational temperature of the mixture using Newton-Rhapson iteration 913 CoeffType delta_T = std::numeric_limits<StateType>::max(); 914 const unsigned int max_iterations = 100; 915 916 const CoeffType dT_reltol= std::numeric_limits<CoeffType>::epsilon() * 100; 917 918 // NOTE: FIN-S uses a hardcoded, absolute tolerance on delta_T of 919 // 1e-8. Using a relative tolerance here of 100*epsilon. 920 for (unsigned int iter = 0; 921 abs(delta_T/T) > dT_reltol && 922 T >= 0.; 923 ++iter) 924 { 925 if (iter == max_iterations) 926 throw FailedNewtonTTvInversion ("ERROR: failed to converge T_from_e_tot!"); 927 928 // compute the residual, defined as the mismatch between the input e_tot and 929 // the value corresponding to the current T iterate 930 typename Antioch::value_type<VectorStateType>::type 931 Re = 0., dRevdT = 0.; 932 933 this->e_and_cv_ve(T, mass_fractions, Re, dRevdT); 934 Re += this->e_tr(T, mass_fractions) + E_0 - e_tot; 935 dRevdT += Cv_tr; 936 937 delta_T = -Re / dRevdT; 938 T += delta_T; 939 } 940 941 return T; 942 } 943 944 template<typename CoeffType> 945 template<typename VectorStateType> 946 inline 947 typename enable_if_c< 948 has_size<VectorStateType>::value, 949 typename Antioch::value_type<VectorStateType>::type 950 >::type T_from_e_tr(const typename Antioch::value_type<VectorStateType>::type & e_tr,const VectorStateType & mass_fractions,typename Antioch::value_type<VectorStateType>::type T)951 StatMechThermodynamics<CoeffType>::T_from_e_tr 952 (const typename Antioch::value_type<VectorStateType>::type& e_tr, 953 const VectorStateType& mass_fractions, 954 typename Antioch::value_type<VectorStateType>::type T ) const 955 { 956 antioch_not_implemented(); 957 } 958 959 template<typename CoeffType> 960 template<typename VectorStateType> 961 inline 962 typename enable_if_c< 963 has_size<VectorStateType>::value, 964 typename Antioch::value_type<VectorStateType>::type 965 >::type T_from_h_tot(const typename Antioch::value_type<VectorStateType>::type & h_tot,const VectorStateType & mass_fractions,typename Antioch::value_type<VectorStateType>::type T)966 StatMechThermodynamics<CoeffType>::T_from_h_tot 967 (const typename Antioch::value_type<VectorStateType>::type& h_tot, 968 const VectorStateType& mass_fractions, 969 typename Antioch::value_type<VectorStateType>::type T) const 970 { 971 antioch_not_implemented(); 972 } 973 974 template<typename CoeffType> 975 template<typename VectorStateType> 976 inline 977 typename enable_if_c< 978 has_size<VectorStateType>::value, 979 typename Antioch::value_type<VectorStateType>::type 980 >::type T_from_h_tot_Tv(const typename Antioch::value_type<VectorStateType>::type & h_tot,const typename Antioch::value_type<VectorStateType>::type & Tv,const VectorStateType & mass_fractions,typename Antioch::value_type<VectorStateType>::type T)981 StatMechThermodynamics<CoeffType>::T_from_h_tot_Tv 982 (const typename Antioch::value_type<VectorStateType>::type& h_tot, 983 const typename Antioch::value_type<VectorStateType>::type& Tv, 984 const VectorStateType& mass_fractions, 985 typename Antioch::value_type<VectorStateType>::type T) const 986 { 987 antioch_not_implemented(); 988 } 989 990 991 template<typename CoeffType> 992 template<typename VectorStateType> 993 inline 994 typename enable_if_c< 995 has_size<VectorStateType>::value, 996 typename Antioch::value_type<VectorStateType>::type 997 >::type s(const typename Antioch::value_type<VectorStateType>::type & T,const typename Antioch::value_type<VectorStateType>::type & p,const VectorStateType & mass_fractions)998 StatMechThermodynamics<CoeffType>::s 999 (const typename Antioch::value_type<VectorStateType>::type& T, 1000 const typename Antioch::value_type<VectorStateType>::type& p, 1001 const VectorStateType& mass_fractions) const 1002 { 1003 antioch_not_implemented(); 1004 } 1005 1006 } // end namespace Antioch 1007 1008 #endif // ANTIOCH_STAT_MECH_THERMO_H 1009