antioch-0.4.0
molecular_binary_diffusion.h
Go to the documentation of this file.
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 #include "antioch_config.h"
29 #ifdef ANTIOCH_HAVE_GSL // if we do not have it, we don't even define the stuff
30 
31 #ifndef ANTIOCH_MOLECULAR_BINARY_DIFFUSION_H
32 #define ANTIOCH_MOLECULAR_BINARY_DIFFUSION_H
33 
34 
35 //Antioch
37 #include "antioch/math_constants.h"
38 #include "antioch/cmath_shims.h"
41 #include "antioch/gsl_spliner.h"
43 
44 //C++
45 #include <vector>
46 
47 namespace Antioch{
48 
69  template <typename CoeffType, typename Interpolator = GSLSpliner>
70  class MolecularBinaryDiffusion : public BinaryDiffusionBase<MolecularBinaryDiffusion<CoeffType,Interpolator>,CoeffType>
71  {
72  typedef unsigned int Species;
73 
74  public:
75  MolecularBinaryDiffusion(const TransportSpecies<CoeffType> & si, const TransportSpecies<CoeffType> & sj);
76  MolecularBinaryDiffusion(const std::vector<TransportSpecies<CoeffType> >& species);
77  ~MolecularBinaryDiffusion();
78 
107  void reset_coeffs(const std::vector<TransportSpecies<CoeffType> > & spec);
108 
109  template <typename StateType>
110  const
111  ANTIOCH_AUTO(StateType)
112  Stockmayer(const StateType & T)
113  ANTIOCH_AUTOFUNC(StateType, ant_sqrt(T) / _interp.interpolated_value(T) )
114 
116  template <typename StateType>
117  const
118  ANTIOCH_AUTO(StateType)
119  binary_diffusion(const StateType & T, const StateType & molar_density) const
120  ANTIOCH_AUTOFUNC(StateType,_coefficient * _interp.interpolated_value(T) / molar_density )
121 
123  const CoeffType & reduced_dipole_moment() const;
124 
126  const CoeffType & reduced_mass() const;
127 
129  const CoeffType & reduced_LJ_diameter() const;
130 
132  const CoeffType & reduced_LJ_depth() const;
133 
135  const CoeffType & xi() const;
136 
137  void print(std::ostream & out = std::cout) const;
138 
139  friend std::ostream & operator<< (std::ostream & out, const MolecularBinaryDiffusion<CoeffType,Interpolator> & bimol_diff)
140  {
141  bimol_diff.print(out);
142  return out;
143  }
144 
146  friend class BinaryDiffusionBase<MolecularBinaryDiffusion<CoeffType,Interpolator>,CoeffType>;
147 
148  protected:
149 
150  void reset_coeffs_impl(const TransportSpecies<CoeffType> & si, const TransportSpecies<CoeffType> & sj);
151 
153 
159  template <typename StateType>
160  void extrapolate_max_temp_impl(const StateType & T);
161 
163  template <typename StateType>
164  const
165  ANTIOCH_AUTO(StateType)
166  op_impl(const StateType & T, const StateType & molar_density) const
167  ANTIOCH_AUTOFUNC(StateType, this->binary_diffusion(T, molar_density))
168 
169  private:
170 
171  MolecularBinaryDiffusion();
172 
173  // convenient method to avoid code duplication
174  void set_coeffs(const std::vector<TransportSpecies<CoeffType> > & spec);
175 
176  template <typename SpeciesType>
177  const CoeffType composed_xi(const TransportSpecies<SpeciesType> & si, const TransportSpecies<SpeciesType> & sj);
178 
179  void build_interpolation();
180 
181  void build_spline(const StockmayerPotential<CoeffType> & surface);
182 
183 
184  Interpolator _interp;
185 
186 // parameter for the couple
187  Species _i,_j;
188  CoeffType _reduced_mass;
189  CoeffType _xi; // simplify calculations
190  CoeffType _reduced_LJ_diameter;
191  CoeffType _reduced_LJ_depth;
192  CoeffType _reduced_dipole_moment;
193 
194  CoeffType _coefficient;
195 
196  };
197 
198  template <typename CoeffType, typename Interpolator>
199  inline
200  MolecularBinaryDiffusion<CoeffType,Interpolator>::~MolecularBinaryDiffusion()
201  {
202  return;
203  }
204 
205  template <typename CoeffType, typename Interpolator>
206  inline
207  MolecularBinaryDiffusion<CoeffType,Interpolator>::MolecularBinaryDiffusion(const TransportSpecies<CoeffType> & si, const TransportSpecies<CoeffType> & sj)
208  : BinaryDiffusionBase<MolecularBinaryDiffusion<CoeffType,Interpolator>,CoeffType>(),
209  _i(si.species()),
210  _j(sj.species()),
211  _reduced_mass((si.species() == sj.species())?CoeffType(0.5) * si.M():(si.M() * sj.M()) / (si.M() + sj.M())), // kg/mol
212  _xi((si.polar() == sj.polar())?1:this->composed_xi(si,sj)),
213  _reduced_LJ_diameter(CoeffType(0.5) * (si.LJ_diameter() + sj.LJ_diameter()) * Units<CoeffType>("ang").get_SI_factor() * ant_pow(_xi,-CoeffType(1)/CoeffType(6))), // 1/2 * (sigma_1 + sigma_2) * xi^(-1/6)
214  _reduced_LJ_depth(ant_sqrt(si.LJ_depth() * sj.LJ_depth()) * _xi * _xi), // sqrt(eps_1 * eps_2) * xi^2
215  _reduced_dipole_moment((CoeffType(1e-7L) * ant_pow(Constants::light_celerity<CoeffType>(),2)) *
216  (si.dipole_moment() * sj.dipole_moment() * ant_pow(Units<CoeffType>("D").get_SI_factor(),2 )) /
217  (2 * _reduced_LJ_depth * Constants::Boltzmann_constant<CoeffType>() * ant_pow(_reduced_LJ_diameter,3) )), // mu^2 / (2*eps*sigma^3)
218 /* = ~ 7.16 10^-25,
219  float can't take it,
220  we cheat on the kb / nAvo division that makes float cry
221 
222  3/16 * sqrt ( 2 * kb / (Navo * pi) )
223 */
224  _coefficient(CoeffType(0.1875e-25L) * ant_sqrt(2 * Constants::Boltzmann_constant<CoeffType>() * CoeffType(1e25) /
225  (_reduced_mass * Constants::Avogadro<CoeffType>() * CoeffType(1e-25L) * Constants::pi<CoeffType>())
226  )
227  / ( _reduced_LJ_diameter * _reduced_LJ_diameter )
228  )
229  {
230 
231  this->build_interpolation();
232 
233  return;
234  }
235 
236  template <typename CoeffType, typename Interpolator>
237  inline
238  MolecularBinaryDiffusion<CoeffType,Interpolator>::MolecularBinaryDiffusion(const std::vector<TransportSpecies<CoeffType> >& species)
239  : BinaryDiffusionBase<MolecularBinaryDiffusion<CoeffType,Interpolator>,CoeffType>(),
240 #ifndef NDEBUG
241  _i(0),
242  _j(0),
243  _reduced_mass(-1),
244  _reduced_dipole_moment(-1),
245  _xi(-1),
246  _reduced_LJ_diameter(-1),
247  _reduced_LJ_depth(-1),
248  _coefficient(0.)
249 #else
250  _i(species[0].species()),
251  _j(species[1].species()),
252  _reduced_mass((species[0].species() == species[1].species())?CoeffType(0.5) * species[0].M():(species[0].M() * species[1].M()) / (species[0].M() + species[1].M())), // kg/mol
253  _xi((species[0].polar() == species[1].polar())?1:this->composed_xi(species[0],species[1])),
254  _reduced_LJ_diameter(CoeffType(0.5) * (species[0].LJ_diameter() + species[1].LJ_diameter()) * Units<CoeffType>("ang").get_SI_factor() * ant_pow(_xi,-CoeffType(1)/CoeffType(6))), // 1/2 * (sigma_1 + sigma_2) * xi^(-1/6)
255  _reduced_LJ_depth(ant_sqrt(species[0].LJ_depth() * species[1].LJ_depth()) *_xi * _xi), // sqrt(eps_1 * eps_2) * xi^2
256  _reduced_dipole_moment((CoeffType(1e-7L) * ant_pow(Constants::light_celerity<CoeffType>(),2)) *
257  (species[0].dipole_moment() * species[1].dipole_moment() * ant_pow(Units<CoeffType>("D").get_SI_factor(),2 )) /
258  (2 * _reduced_LJ_depth * Constants::Boltzmann_constant<CoeffType>() * ant_pow(_reduced_LJ_diameter,3) )), // mu^2 / (2*eps*sigma^3)
259  _coefficient(CoeffType(0.1875e-25L) * ant_sqrt(2 * Constants::Boltzmann_constant<CoeffType>() * CoeffType(1e25) /
260  (_reduced_mass * Constants::Avogadro<CoeffType>() * CoeffType(1e-25L) * Constants::pi<CoeffType>())
261  )
262  / ( _reduced_LJ_diameter * _reduced_LJ_diameter )
263  )
264 #endif
265  {
266 #ifndef NDEBUG
267  antioch_assert_equal_to(species.size(),2);
268 
269  this->set_coeffs(species);
270 #endif
271 
272  this->build_interpolation();
273  }
274 
275  template <typename CoeffType, typename Interpolator>
276  inline
277  void MolecularBinaryDiffusion<CoeffType,Interpolator>::reset_coeffs(const std::vector<TransportSpecies<CoeffType> >& species)
278  {
279  antioch_assert_equal_to(species.size(),2);
280 
281  this->reset_coeffs(species[0],species[1]);
282  }
283 
284  template <typename CoeffType, typename Interpolator>
285  inline
286  void MolecularBinaryDiffusion<CoeffType,Interpolator>::reset_coeffs_impl(const TransportSpecies<CoeffType>& si,
287  const TransportSpecies<CoeffType>& sj)
288  {
289  _i = si.species();
290  _j = sj.species();
291  _reduced_mass = (si.species() == sj.species())?CoeffType(0.5) * si.M():(si.M() * sj.M()) / (si.M() + sj.M());
292  _xi = (si.polar() == sj.polar())?1:this->composed_xi(si,sj);
293  _reduced_LJ_diameter = CoeffType(0.5) * (si.LJ_diameter() + sj.LJ_diameter()) * Units<CoeffType>("ang").get_SI_factor() * ant_pow(_xi,-CoeffType(1)/CoeffType(6));
294  _reduced_LJ_depth = ant_sqrt(si.LJ_depth() * sj.LJ_depth()) *_xi * _xi;
295  _reduced_dipole_moment = CoeffType(1e-7L) * ant_pow(Constants::light_celerity<CoeffType>(),2) * // * 1/(4*pi * eps_0) = 10^-7 * c^2
296  (si.dipole_moment() * sj.dipole_moment()) * ant_pow(Units<CoeffType>("D").get_SI_factor(),2)
297  / ((2 * _reduced_LJ_depth * Constants::Boltzmann_constant<CoeffType>() * ant_pow(_reduced_LJ_diameter,3)));
298 /* = ~ 7.16 10^-25,
299  float can't take it,
300  we cheat on the kb / nAvo division that makes float cry
301 
302  3/16 * sqrt ( 2 * kb / (Navo * pi) )
303 */
304  _coefficient = CoeffType(0.1875e-25L) * ant_sqrt(2 * Constants::Boltzmann_constant<CoeffType>() * CoeffType(1e25) /
305  (_reduced_mass * Constants::Avogadro<CoeffType>() * CoeffType(1e-25) * Constants::pi<CoeffType>())
306  )
307  / ( _reduced_LJ_diameter * _reduced_LJ_diameter );
308  }
309 
310  template <typename CoeffType, typename Interpolator>
311  inline
312  void MolecularBinaryDiffusion<CoeffType,Interpolator>::set_coeffs(const std::vector<TransportSpecies<CoeffType> >& species)
313  {
314  antioch_assert_equal_to(species.size(),2);
315 
316  this->reset_coeffs(species[0],species[1]);
317  }
318 
319  template <typename CoeffType, typename Interpolator>
320  template <typename SpeciesType>
321  inline
322  const CoeffType MolecularBinaryDiffusion<CoeffType,Interpolator>::composed_xi(const TransportSpecies<SpeciesType> & si, const TransportSpecies<SpeciesType> & sj)
323  {
324  const TransportSpecies<SpeciesType> &n = (si.polar())?sj:si;
325  const TransportSpecies<SpeciesType> &p = (si.polar())?si:sj;
326 
327  CoeffType pol = n.polarizability() / ant_pow(n.LJ_diameter(),3); //ang^3 / ang^3 -> cancel out
328  CoeffType dipole = p.dipole_moment() * Units<CoeffType>("D").get_SI_factor()
329  / ant_sqrt(4 * Constants::pi<CoeffType>() * Constants::vacuum_permittivity<CoeffType>()
330  * p.LJ_depth() * ant_pow(p.LJ_diameter(),3) );
331 
332  return 1 + CoeffType(0.25L) * pol * dipole * ant_sqrt(p.LJ_depth()/n.LJ_depth());
333  }
334 
335  template <typename CoeffType, typename Interpolator>
336  inline
337  void MolecularBinaryDiffusion<CoeffType,Interpolator>::build_spline(const StockmayerPotential<CoeffType> & surface)
338  {
339  std::vector<CoeffType> interp_surf(surface.log_temperature().size(),0);
340  std::vector<CoeffType> rescaled_temp(surface.log_temperature().size(),0);
341  for(unsigned int iT = 0; iT < surface.temperature().size(); iT++)
342  {
343  Interpolator spline(surface.delta(),surface.omega_1_1()[iT]);
344  interp_surf[iT] = ant_sqrt(surface.temperature()[iT] * _reduced_LJ_depth) /
345  spline.interpolated_value(_reduced_dipole_moment); // splining sqrt(T) / Omega<(1,1)>(log(T*))
346  rescaled_temp[iT] = surface.temperature()[iT] * _reduced_LJ_depth;
347  }
348 
349  _interp.spline_delete();
350  _interp.spline_init(rescaled_temp,interp_surf);
351  }
352 
353  template <typename CoeffType, typename Interpolator>
354  template <typename StateType>
355  inline
356  void MolecularBinaryDiffusion<CoeffType,Interpolator>::extrapolate_max_temp_impl(const StateType& Tmax)
357  {
358  StockmayerPotential<CoeffType> surface;
359  // Stockmayer is where the test is performed
360  surface.extrapolate_to(Tmax/_reduced_LJ_depth);
361  build_spline(surface);
362  }
363 
364  template <typename CoeffType, typename Interpolator>
365  inline
366  void MolecularBinaryDiffusion<CoeffType,Interpolator>::build_interpolation()
367  {
368  build_spline(StockmayerPotential<CoeffType>());
369  }
370 
371  template <typename CoeffType, typename Interpolator>
372  inline
373  void MolecularBinaryDiffusion<CoeffType,Interpolator>::print(std::ostream & out) const
374  {
375  out << "Molecular binary diffusion of couple " << _i << "," << _j << ":\n"
376  << " reduced mass = " << _reduced_mass << "\n"
377  << " reduced dipole moment = " << _reduced_dipole_moment << "\n"
378  << " xi = " << _xi << "\n"
379  << " reduced LJ diameter = " << _reduced_LJ_diameter << "\n"
380  << " reduced LJ depth = " << _reduced_LJ_depth << "\n"
381  << " [coefficient = " << _coefficient << "]"
382  << std::endl;
383  }
384 
385  template <typename CoeffType, typename Interpolator>
386  inline
387  const CoeffType & MolecularBinaryDiffusion<CoeffType,Interpolator>::reduced_dipole_moment() const
388  {
389  return _reduced_dipole_moment;
390  }
391 
392  template <typename CoeffType, typename Interpolator>
393  inline
394  const CoeffType & MolecularBinaryDiffusion<CoeffType,Interpolator>::reduced_mass() const
395  {
396  return _reduced_mass;
397  }
398 
399  template <typename CoeffType, typename Interpolator>
400  inline
401  const CoeffType & MolecularBinaryDiffusion<CoeffType,Interpolator>::reduced_LJ_diameter() const
402  {
403  return _reduced_LJ_diameter;
404  }
405 
406  template <typename CoeffType, typename Interpolator>
407  inline
408  const CoeffType & MolecularBinaryDiffusion<CoeffType,Interpolator>::reduced_LJ_depth() const
409  {
410  return _reduced_LJ_depth;
411  }
412 
413  template <typename CoeffType, typename Interpolator>
414  inline
415  const CoeffType & MolecularBinaryDiffusion<CoeffType,Interpolator>::xi() const
416  {
417  return _xi;
418  }
419 
420 
421 }
422 
423 #endif // ANTIOCH_BIMOL_DIFF
424 
425 #endif // ifdef ANTIOCH_HAVE_GSL
unsigned int Species
#define antioch_assert_equal_to(expr1, expr2)
CoeffType Avogadro()
Avogadro's number, particles per mole.
CoeffType Boltzmann_constant()
Boltzmann constant 1.380 6488 x 10-23 J/K (http://physics.nist.gov/cgi-bin/cuu/Value?k)
CoeffType light_celerity()
light celerity, m/s
#define ANTIOCH_AUTOFUNC(Type, Expr)
const ANTIOCH_AUTO(StateType) KineticsTheoryThermalConductivity< ThermoEvaluator
The parameters are reduced parameters.
CoeffType pi()
pi

Generated on Thu Jul 7 2016 11:09:45 for antioch-0.4.0 by  doxygen 1.8.8