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