GRINS-0.8.0
absorption_coeff.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // GRINS - General Reacting Incompressible Navier-Stokes
5 //
6 // Copyright (C) 2014-2017 Paul T. Bauman, Roy H. Stogner
7 // Copyright (C) 2010-2013 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 
26 // This class
27 #include "grins/absorption_coeff.h"
28 
29 // GRINS
30 #include "grins/variable_warehouse.h"
31 #include "grins/physical_constants.h"
32 #include "grins/math_constants.h"
33 
34 #if GRINS_HAVE_ANTIOCH
35 #include "grins/antioch_chemistry.h"
36 #endif
37 
38 #if GRINS_HAVE_CANTERA
39 #include "grins/cantera_mixture.h"
40 #endif
41 
42 // libMesh
43 #include "libmesh/fe.h"
44 #include "libmesh/elem.h"
45 #include "libmesh/fe_interface.h"
46 
47 namespace GRINS
48 {
49  template<typename Chemistry>
50  AbsorptionCoeff<Chemistry>::AbsorptionCoeff(SharedPtr<Chemistry> & chem, SharedPtr<HITRAN> & hitran,
51  libMesh::Real nu_min, libMesh::Real nu_max,
52  libMesh::Real desired_nu, const std::string & species,
53  libMesh::Real thermo_pressure)
54  : _chemistry(chem),
55  _hitran(hitran),
56  _nu(desired_nu),
57  _T_var(GRINSPrivate::VariableWarehouse::get_variable_subclass<PrimitiveTempFEVariables>("Temperature")),
58  _P_var(GRINSPrivate::VariableWarehouse::get_variable_subclass<PressureFEVariable>("Pressure")),
59  _Y_var(GRINSPrivate::VariableWarehouse::get_variable_subclass<SpeciesMassFractionsVariable>("SpeciesMassFractions")),
60  _T0(296), // [K]
61  _Pref(1), // [atm]
62  _rad_coeff(Constants::second_rad_const * 100) // [cm K]
63  {
64  // sanity checks
65  if ( (nu_min>nu_max) || (desired_nu>nu_max) || (desired_nu<nu_min) )
66  {
67  std::stringstream ss;
68  ss <<"Invalid specification of wavenumber range:" <<std::endl
69  <<"nu_min: " <<nu_min <<std::endl
70  <<"nu_max: " <<nu_max <<std::endl
71  <<"desired_nu: " <<desired_nu <<std::endl;
72  libmesh_error_msg(ss.str());
73  }
74 
75  _species_idx = _chemistry->species_index(species);
76  unsigned int data_size = _hitran->get_data_size();
77 
78  bool min_flag = false;
79 
80  for (unsigned int i=0; i<data_size; i++)
81  {
82  if (_hitran->nu0(i) > nu_min)
83  {
84  _min_index = i;
85  min_flag = true;
86  break;
87  }
88  }
89 
90  if (!min_flag)
91  {
92  std::stringstream ss;
93  ss <<"Minimum wavenumber " <<nu_min <<" is greater than the maximum wavenumber in provided HITRAN data";
94  libmesh_error_msg(ss.str());
95  }
96 
97  bool max_flag = false;
98 
99  for (unsigned int i=data_size-1; i>=0; i--)
100  {
101  if (_hitran->nu0(i) < nu_max)
102  {
103  _max_index = i;
104  max_flag = true;
105  break;
106  }
107  }
108 
109  if (!max_flag)
110  _max_index = data_size-1;
111 
112  if (thermo_pressure == -1.0) {
113  _calc_thermo_pressure = true;
114  _thermo_pressure = -1.0; // not needed in this case
115  } else {
116  _calc_thermo_pressure = false;
117  _thermo_pressure = thermo_pressure;
118  }
119 
120  this->init_voigt();
121  }
122 
123  template<typename Chemistry>
124  libMesh::Real AbsorptionCoeff<Chemistry>::operator()(const libMesh::FEMContext& context, const libMesh::Point& qp_xyz, const libMesh::Real /*t*/)
125  {
126  START_LOG("operator()","AbsorptionCoeff");
127  libMesh::Real T,p,thermo_p; // temperature, hydrostatic pressure, thermodynamic pressure
128  std::vector<libMesh::Real> Y(_chemistry->n_species()); // mass fractions
129 
130  for (unsigned int s=0; s<_chemistry->n_species(); s++)
131  context.point_value(_Y_var.species(s), qp_xyz, Y[s]);
132 
133  context.point_value(_T_var.T(), qp_xyz, T); // [K]
134 
135  if (_calc_thermo_pressure) {
136  libmesh_not_implemented();
137  } else {
138  thermo_p = _thermo_pressure;
139  }
140 
141  context.point_value(_P_var.p(), qp_xyz, p); // [Pa]
142 
143  libMesh::Real P = p + thermo_p; // total pressure [Pa]
144  libmesh_assert_greater(P,0.0);
145 
146  libMesh::Real kv = 0.0;
147 
148  for (unsigned int i=_min_index; i<=_max_index; i++)
149  kv += this->kv(T,P,Y,i);
150 
151  STOP_LOG("operator()","AbsorptionCoeff");
152  return kv;
153  }
154 
155  template<typename Chemistry>
156  void AbsorptionCoeff<Chemistry>::operator()( const libMesh::FEMContext & /*context*/,
157  const libMesh::Point & /*p*/,
158  const libMesh::Real /*time*/,
159  libMesh::DenseVector<libMesh::Real> & /*output*/)
160  {
161  libmesh_not_implemented();
162  }
163 
164  template<typename Chemistry>
165  void AbsorptionCoeff<Chemistry>::derivatives( libMesh::FEMContext & context,
166  const libMesh::Point & qp_xyz,
167  const libMesh::Real & JxW,
168  const unsigned int qoi_index,
169  const libMesh::Real /*time*/)
170  {
171  START_LOG("derivatives()","AbsorptionCoeff");
172 
173  libMesh::DenseSubVector<libMesh::Number> & dQdT = context.get_qoi_derivatives(qoi_index, _T_var.T());
174  libMesh::DenseSubVector<libMesh::Number> & dQdP = context.get_qoi_derivatives(qoi_index, _P_var.p());
175  libMesh::DenseSubVector<libMesh::Number> & dQdYs = context.get_qoi_derivatives(qoi_index, _Y_var.species(_species_idx));
176 
177  // need to map the physical coordinates of QP to reference coordinates
178  libMesh::Elem & main_elem = context.get_elem();
179  libMesh::Point qp_ref = libMesh::FEInterface::inverse_map(main_elem.dim(),main_elem.type(),&main_elem,qp_xyz);
180 
181  std::vector<libMesh::Point> qp(1);
182  qp[0] = qp_ref;
183 
184  libMesh::UniquePtr< libMesh::FEBase > T_fe = libMesh::FEGenericBase<libMesh::Real>::build(context.get_elem().dim(), context.get_element_fe(_T_var.T())->get_fe_type() );
185  const std::vector<std::vector<libMesh::Real> > & T_phi =T_fe->get_phi();
186  T_fe->reinit(&(context.get_elem()),&qp);
187 
188  libMesh::UniquePtr< libMesh::FEBase > P_fe = libMesh::FEGenericBase<libMesh::Real>::build(context.get_elem().dim(), context.get_element_fe(_P_var.p())->get_fe_type() );
189  const std::vector<std::vector<libMesh::Real> > & P_phi = P_fe->get_phi();
190  P_fe->reinit(&(context.get_elem()),&qp);
191 
192  libMesh::UniquePtr< libMesh::FEBase > Ys_fe = libMesh::FEGenericBase<libMesh::Real>::build(context.get_elem().dim(), context.get_element_fe(_Y_var.species(_species_idx))->get_fe_type() );
193  const std::vector<std::vector<libMesh::Real> > & Ys_phi = Ys_fe->get_phi();
194  Ys_fe->reinit(&(context.get_elem()),&qp);
195 
196  libMesh::Real T,p,thermo_p; // temperature, hydrostatic pressure, thermodynamic pressure
197  std::vector<libMesh::Real> Y(_chemistry->n_species()); // mass fractions
198 
199  context.point_value(_T_var.T(), qp_xyz, T); // [K]
200 
201  if (_calc_thermo_pressure) {
202  libmesh_not_implemented();
203  } else {
204  thermo_p = _thermo_pressure;
205  }
206 
207  context.point_value(_P_var.p(), qp_xyz, p); // [Pa]
208 
209  libMesh::Real P = p + thermo_p; // total pressure [Pa]
210  libmesh_assert_greater(P,0.0);
211 
212  // all mass fractions needed to get M_mix
213  for (unsigned int s=0; s<_chemistry->n_species(); s++)
214  context.point_value(_Y_var.species(s), qp_xyz, Y[s]);
215 
216  for (unsigned int i=_min_index; i<=_max_index; i++)
217  {
218  // no velocity dependence
219 
220  // temperature deriv
221  for (unsigned int j=0; j<dQdT.size(); j++)
222  dQdT(j) += d_kv_dT(T,P,Y,i)*JxW * T_phi[j][0];
223 
224  // pressure deriv
225  for (unsigned int j=0; j<dQdP.size(); j++)
226  dQdP(j) += d_kv_dP(T,P,Y,i)*JxW * P_phi[j][0];
227 
228  // mass fraction deriv
229  for (unsigned int j=0; j<dQdYs.size(); j++)
230  dQdYs(j) += d_kv_dY(T,P,Y,i)*JxW * Ys_phi[j][0];
231  }
232 
233  STOP_LOG("derivatives()","AbsorptionCoeff");
234  }
235 
236  template<typename Chemistry>
237  libMesh::UniquePtr<libMesh::FEMFunctionBase<libMesh::Real> > AbsorptionCoeff<Chemistry>::clone() const
238  {
239  libmesh_not_implemented();
240  }
241 
242 
243  template<typename Chemistry>
244  libMesh::Real AbsorptionCoeff<Chemistry>::kv(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
245  {
246  // linestrength
247  libMesh::Real S = this->Sw(T,P,i);
248 
249  // Voigt profile [cm^-1]
250  libMesh::Real phi_V = this->voigt(T,P,Y,i);
251 
252  libMesh::Real M_mix = _chemistry->M_mix(Y);
253  libMesh::Real X = _chemistry->X(_species_idx,M_mix,Y[_species_idx]);
254 
255  // absorption coefficient [cm^-1]
256  return S*(P/Constants::atmosphere_Pa)*X*phi_V;
257  }
258 
259  template<typename Chemistry>
260  libMesh::Real AbsorptionCoeff<Chemistry>::d_kv_dT(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
261  {
262  libMesh::Real dS = dS_dT(T,P,i);
263  libMesh::Real dV = d_voigt_dT(T,P,Y,i);
264 
265  libMesh::Real S = this->Sw(T,P,i);
266  libMesh::Real V = this->voigt(T,P,Y,i);
267  libMesh::Real M_mix = _chemistry->M_mix(Y);
268  libMesh::Real X = _chemistry->X(_species_idx,M_mix,Y[_species_idx]);
269 
270  return X*(P/Constants::atmosphere_Pa) * ( S*dV + dS*V );
271  }
272 
273  template<typename Chemistry>
274  libMesh::Real AbsorptionCoeff<Chemistry>::d_kv_dP(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
275  {
276  libMesh::Real dS = dS_dP(T,P,i);
277  libMesh::Real dV = d_voigt_dP(T,P,Y,i);
278 
279  libMesh::Real S = this->Sw(T,P,i);
280  libMesh::Real V = this->voigt(T,P,Y,i);
281  libMesh::Real M_mix = _chemistry->M_mix(Y);
282  libMesh::Real X = _chemistry->X(_species_idx,M_mix,Y[_species_idx]);
283 
284  return X * ( S*(P/Constants::atmosphere_Pa)*dV + S*(1.0/Constants::atmosphere_Pa)*V + dS*P*V );
285  }
286 
287  template<typename Chemistry>
288  libMesh::Real AbsorptionCoeff<Chemistry>::d_kv_dY(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
289  {
290  libMesh::Real dX = dX_dY(Y);
291  libMesh::Real dV = d_voigt_dY(T,P,Y,i);
292 
293  libMesh::Real S = this->Sw(T,P,i);
294  libMesh::Real V = this->voigt(T,P,Y,i);
295  libMesh::Real M_mix = _chemistry->M_mix(Y);
296  libMesh::Real X = _chemistry->X(_species_idx,M_mix,Y[_species_idx]);
297 
298  return S*(P/Constants::atmosphere_Pa) * ( X*dV + dX*V );
299  }
300 
301 
302  template<typename Chemistry>
303  libMesh::Real AbsorptionCoeff<Chemistry>::Sw(libMesh::Real T, libMesh::Real P, unsigned int i)
304  {
305  // isotopologue
306  unsigned int iso = _hitran->isotopologue(i);
307 
308  // linecenter wavenumber
309  libMesh::Real nu = this->get_nu(P,i);
310 
311  // linestrength
312  libMesh::Real sw0 = _hitran->sw(i);
313 
314  // lower state energy of transition
315  libMesh::Real E = _hitran->elower(i);
316 
317  // partition function at reference temp
318  libMesh::Real QT0 = _hitran->partition_function(_T0,iso);
319 
320  // partition function at current temp
321  libMesh::Real QT = _hitran->partition_function(T,iso);
322 
323  libMesh::Real a = sw0;
324  libMesh::Real b = (QT0/QT);
325  libMesh::Real c = std::exp(-E*_rad_coeff*( (1.0/T) - (1.0/_T0) ));
326  libMesh::Real d = ( 1.0-std::exp(-_rad_coeff*nu/T) );
327  libMesh::Real e = std::pow(1.0-std::exp(-_rad_coeff*nu/_T0),-1.0);
328 
329  libMesh::Real S = a*b*c*d*e;
330 
331  // convert linestrength units to [cm^-2 atm^-1]
332  libMesh::Real loschmidt = (Constants::atmosphere_Pa)/(T*Constants::Boltzmann*1.0e6);
333  S = S*loschmidt;
334 
335  return S;
336  }
337 
338  template<typename Chemistry>
339  libMesh::Real AbsorptionCoeff<Chemistry>::dS_dT(libMesh::Real T, libMesh::Real P, unsigned int i)
340  {
341  libMesh::Real iso = _hitran->isotopologue(i);
342  libMesh::Real sw = _hitran->sw(i);
343  libMesh::Real E = _hitran->elower(i);
344  libMesh::Real QT0 = _hitran->partition_function(_T0,iso);
345  libMesh::Real QT = _hitran->partition_function(T,iso);
346  libMesh::Real dQT = this->dQ_dT(T,iso);
347  libMesh::Real nu = this->get_nu(P,i);
348 
349  libMesh::Real constants = sw * QT0 * std::pow( (1.0-std::exp(-_rad_coeff*nu/_T0)), -1.0) * Constants::atmosphere_Pa/(Constants::Boltzmann*1.0e6);
350  libMesh::Real A = 1.0/T;
351  libMesh::Real dA = -1.0/(T*T);
352 
353  libMesh::Real B = 1.0/QT;
354  libMesh::Real dB = -1.0/(QT*QT) * dQT;
355 
356  libMesh::Real C = std::exp(-_rad_coeff*E* (1.0/T - 1.0/_T0) );
357  libMesh::Real dC = std::exp(-_rad_coeff*E* (1.0/T - 1.0/_T0) ) * (_rad_coeff*E/(T*T));
358 
359  libMesh::Real D = 1.0 - std::exp(-_rad_coeff*nu/T);
360  libMesh::Real dD = -std::exp(-_rad_coeff*nu/T) * (_rad_coeff*nu/(T*T));
361 
362  return constants * ( A*B*C*dD + A*B*dC*D + A*dB*C*D + dA*B*C*D );
363  }
364 
365  template<typename Chemistry>
366  libMesh::Real AbsorptionCoeff<Chemistry>::dS_dP(libMesh::Real T, libMesh::Real P, unsigned int i)
367  {
368  libMesh::Real iso = _hitran->isotopologue(i);
369  libMesh::Real sw = _hitran->sw(i);
370  libMesh::Real E = _hitran->elower(i);
371  libMesh::Real QT0 = _hitran->partition_function(_T0,iso);
372  libMesh::Real QT = _hitran->partition_function(T,iso);
373  libMesh::Real nu = this->get_nu(P,i);
374  libMesh::Real dnu = this->d_nu_dP(i);
375 
376  libMesh::Real constants = sw * QT0/QT * std::exp(-_rad_coeff*E*(1.0/T - 1.0/_T0)) * Constants::atmosphere_Pa/(Constants::Boltzmann*1.0e6 * T);
377 
378  libMesh::Real A = 1.0 - std::exp(-_rad_coeff*nu/T);
379  libMesh::Real dA = -std::exp(-_rad_coeff*nu/T) * (-_rad_coeff/T) * dnu;
380 
381  libMesh::Real B = std::pow( (1.0-std::exp(-_rad_coeff*nu/T)), -1.0);
382  libMesh::Real dB = -std::pow( (1.0-std::exp(-_rad_coeff*nu/T)), -2.0) * -std::exp(-_rad_coeff*nu/_T0) * (-_rad_coeff/_T0) * dnu;
383 
384  return constants * ( A*dB + dA*B );
385  }
386 
387 
388  template<typename Chemistry>
389  libMesh::Real AbsorptionCoeff<Chemistry>::nu_D(libMesh::Real T, libMesh::Real P, unsigned int i)
390  {
391  libMesh::Real k = Constants::Boltzmann;
392  libMesh::Real c = Constants::c_vacuum;
393  libMesh::Real NA = Constants::Avogadro;
394  libMesh::Real M = _chemistry->M(_species_idx);
395  libMesh::Real nu = this->get_nu(P,i);
396 
397  return (nu/c)*std::sqrt( ( 8.0*k*T*std::log(2.0) )/( M/NA ) );
398  }
399 
400  template<typename Chemistry>
401  libMesh::Real AbsorptionCoeff<Chemistry>::d_nuD_dT(libMesh::Real T, libMesh::Real P, unsigned int i)
402  {
403  libMesh::Real k = Constants::Boltzmann;
404  libMesh::Real c = Constants::c_vacuum;
405  libMesh::Real NA = Constants::Avogadro;
406  libMesh::Real M = _chemistry->M(_species_idx);
407  libMesh::Real nu = this->get_nu(P,i);
408 
409  return (nu/c) * std::sqrt( (8.0*k*std::log(2.0))/(M/NA) ) * 0.5 * 1.0/(std::sqrt(T));
410  }
411 
412  template<typename Chemistry>
413  libMesh::Real AbsorptionCoeff<Chemistry>::d_nuD_dP(libMesh::Real T, unsigned int i)
414  {
415  libMesh::Real k = Constants::Boltzmann;
416  libMesh::Real c = Constants::c_vacuum;
417  libMesh::Real NA = Constants::Avogadro;
418  libMesh::Real M = _chemistry->M(_species_idx);
419  libMesh::Real dnu = this->d_nu_dP(i);
420 
421  return (1.0/c)*std::sqrt( ( 8.0*k*T*std::log(2.0) )/( M/NA ) ) * dnu;
422  }
423 
424 
425  template<typename Chemistry>
426  libMesh::Real AbsorptionCoeff<Chemistry>::nu_C(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
427  {
428  libMesh::Real g_self = _hitran->gamma_self(i);
429  libMesh::Real g_air = _hitran->gamma_air(i);
430  libMesh::Real n = _hitran->n_air(i);
431 
432  libMesh::Real M_mix = _chemistry->M_mix(Y);
433  libMesh::Real X = _chemistry->X(_species_idx,M_mix,Y[_species_idx]);
434 
435  return 2.0*(P/Constants::atmosphere_Pa)*std::pow(_T0/T,n)*( X*g_self + (1.0-X)*g_air );
436  }
437 
438  template<typename Chemistry>
439  libMesh::Real AbsorptionCoeff<Chemistry>::d_nuC_dT(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
440  {
441  libMesh::Real g_self = _hitran->gamma_self(i);
442  libMesh::Real g_air = _hitran->gamma_air(i);
443  libMesh::Real n = _hitran->n_air(i);
444 
445  libMesh::Real M_mix = _chemistry->M_mix(Y);
446  libMesh::Real X = _chemistry->X(_species_idx,M_mix,Y[_species_idx]);
447 
448  return 2.0*(P/Constants::atmosphere_Pa) * n*std::pow(_T0/T,n-1.0)*(-_T0/(T*T)) * ( X*g_self + (1.0-X)*g_air );
449  }
450 
451  template<typename Chemistry>
452  libMesh::Real AbsorptionCoeff<Chemistry>::d_nuC_dP(libMesh::Real T, std::vector<libMesh::Real> Y, unsigned int i)
453  {
454  libMesh::Real g_self = _hitran->gamma_self(i);
455  libMesh::Real g_air = _hitran->gamma_air(i);
456  libMesh::Real n = _hitran->n_air(i);
457 
458  libMesh::Real M_mix = _chemistry->M_mix(Y);
459  libMesh::Real X = _chemistry->X(_species_idx,M_mix,Y[_species_idx]);
460 
461  return 2.0*(1.0/Constants::atmosphere_Pa)*std::pow(_T0/T,n)*( X*g_self + (1.0-X)*g_air );
462  }
463 
464  template<typename Chemistry>
465  libMesh::Real AbsorptionCoeff<Chemistry>::d_nuC_dY(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
466  {
467  libMesh::Real g_self = _hitran->gamma_self(i);
468  libMesh::Real g_air = _hitran->gamma_air(i);
469  libMesh::Real n = _hitran->n_air(i);
470  libMesh::Real dX = dX_dY(Y);
471 
472  return 2.0*(P/Constants::atmosphere_Pa)*std::pow(_T0/T,n)*( dX*g_self - dX*g_air );
473  }
474 
475 
476  template<typename Chemistry>
477  libMesh::Real AbsorptionCoeff<Chemistry>::voigt(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
478  {
479  libMesh::Real nu_D = this->nu_D(T,P,i);
480 
481  libMesh::Real a = this->voigt_a(T,P,Y,i);
482  libMesh::Real w = this->voigt_w(T,P,i);
483 
484  // Voigt coefficient
485  libMesh::Real V = 0.0;
486 
487  for(int i=0; i<4; i++)
488  {
489  libMesh::Real Ai = _voigt_coeffs[0][i];
490  libMesh::Real Bi = _voigt_coeffs[1][i];
491  libMesh::Real Ci = _voigt_coeffs[2][i];
492  libMesh::Real Di = _voigt_coeffs[3][i];
493 
494  libMesh::Real aAi = a-Ai;
495  libMesh::Real wBi = w-Bi;
496  libMesh::Real aAi2 = std::pow(a-Ai,2.0);
497  libMesh::Real wBi2 = std::pow(w-Bi,2.0);
498 
499  V += ( Ci*aAi + Di*wBi )/( aAi2 + wBi2 );
500  }
501 
502  libMesh::Real phi_V = (2.0*std::sqrt(std::log(2.0)))/(std::sqrt(Constants::pi)*nu_D)*V;
503 
504  return phi_V;
505  }
506 
507  template<typename Chemistry>
508  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_dT(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
509  {
510  libMesh::Real nu_D = this->nu_D(T,P,i);
511  libMesh::Real dnu_D = d_nuD_dT(T,P,i);
512 
513  libMesh::Real a = this->voigt_a(T,P,Y,i);
514  libMesh::Real da = this->d_voigt_a_dT(T,P,Y,i);
515 
516  libMesh::Real w = this->voigt_w(T,P,i);
517  libMesh::Real dw = this->d_voigt_w_dT(T,P,i);
518 
519  // Voigt coefficient
520  libMesh::Real V = 0.0;
521  libMesh::Real dV = 0.0;
522 
523  for(int i=0; i<4; i++)
524  {
525  libMesh::Real Ai = _voigt_coeffs[0][i];
526  libMesh::Real Bi = _voigt_coeffs[1][i];
527  libMesh::Real Ci = _voigt_coeffs[2][i];
528  libMesh::Real Di = _voigt_coeffs[3][i];
529 
530  libMesh::Real aAi = a-Ai;
531  libMesh::Real wBi = w-Bi;
532  libMesh::Real aAi2 = std::pow(a-Ai,2.0);
533  libMesh::Real wBi2 = std::pow(w-Bi,2.0);
534 
535  V += ( Ci*aAi + Di*wBi )/( aAi2 + wBi2 );
536 
537  dV += ( (aAi2+wBi2)*(Ci*da + Di*dw) - (Ci*aAi + Di*wBi)*(2.0*aAi*da + 2.0*wBi*dw) )/( std::pow(aAi2+wBi2,2.0) );
538  }
539 
540  libMesh::Real constants = 2.0*std::sqrt(std::log(2.0))/std::sqrt(Constants::pi);
541  libMesh::Real C = 1.0/nu_D;
542  libMesh::Real dC = -1.0/(nu_D*nu_D) * dnu_D;
543 
544  return constants * ( C*dV + dC*V );
545  }
546 
547  template<typename Chemistry>
548  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_dP(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
549  {
550  libMesh::Real nu_D = this->nu_D(T,P,i);
551  libMesh::Real dnu_D = d_nuD_dP(T,i);
552 
553  libMesh::Real a = this->voigt_a(T,P,Y,i);
554  libMesh::Real da = this->d_voigt_a_dP(T,P,Y,i);
555 
556  libMesh::Real w = this->voigt_w(T,P,i);
557  libMesh::Real dw = this->d_voigt_w_dP(T,P,i);
558 
559  // Voigt coefficient
560  libMesh::Real V = 0.0;
561  libMesh::Real dV = 0.0;
562 
563  for(int i=0; i<4; i++)
564  {
565  libMesh::Real Ai = _voigt_coeffs[0][i];
566  libMesh::Real Bi = _voigt_coeffs[1][i];
567  libMesh::Real Ci = _voigt_coeffs[2][i];
568  libMesh::Real Di = _voigt_coeffs[3][i];
569 
570  libMesh::Real aAi = a-Ai;
571  libMesh::Real wBi = w-Bi;
572  libMesh::Real aAi2 = std::pow(a-Ai,2.0);
573  libMesh::Real wBi2 = std::pow(w-Bi,2.0);
574 
575  V += ( Ci*aAi + Di*wBi )/( aAi2 + wBi2 );
576 
577  dV += ( (aAi2+wBi2)*(Ci*da + Di*dw) - (Ci*aAi + Di*wBi)*(2.0*aAi*da + 2.0*wBi*dw) )/( std::pow(aAi2+wBi2,2.0) );
578  }
579 
580  libMesh::Real constants = 2.0*std::sqrt(std::log(2.0))/std::sqrt(Constants::pi);
581  libMesh::Real C = 1.0/nu_D;
582  libMesh::Real dC = -1.0/(nu_D*nu_D) * dnu_D;
583 
584  return constants * ( C*dV + dC*V );
585  }
586 
587  template<typename Chemistry>
588  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_dY(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
589  {
590  libMesh::Real nu_D = this->nu_D(T,P,i);
591 
592  libMesh::Real a = this->voigt_a(T,P,Y,i);
593  libMesh::Real da = this->d_voigt_a_dY(T,P,Y,i);
594 
595  libMesh::Real w = this->voigt_w(T,P,i);
596 
597  // Voigt coefficient
598  libMesh::Real dV = 0.0;
599 
600  for(int i=0; i<4; i++)
601  {
602  libMesh::Real Ai = _voigt_coeffs[0][i];
603  libMesh::Real Bi = _voigt_coeffs[1][i];
604  libMesh::Real Ci = _voigt_coeffs[2][i];
605  libMesh::Real Di = _voigt_coeffs[3][i];
606 
607  libMesh::Real aAi = a-Ai;
608  libMesh::Real wBi = w-Bi;
609  libMesh::Real aAi2 = std::pow(a-Ai,2.0);
610  libMesh::Real wBi2 = std::pow(w-Bi,2.0);
611 
612  dV += ( (aAi2+wBi2)*(Ci*da) - (Ci*aAi + Di*wBi)*(2.0*aAi*da) )/( std::pow(aAi2+wBi2,2.0) );
613  }
614 
615  libMesh::Real constants = 2.0*std::sqrt(std::log(2.0))/(std::sqrt(Constants::pi)*nu_D);
616 
617  return constants * dV;
618  }
619 
620  template<typename Chemistry>
621  void AbsorptionCoeff<Chemistry>::init_voigt() {
622  _voigt_coeffs.resize(4);
623  for (int i=0; i<4; i++)
624  _voigt_coeffs[i].resize(4);
625 
626  _voigt_coeffs[0][0] = -1.2150;
627  _voigt_coeffs[0][1] = -1.3509;
628  _voigt_coeffs[0][2] = -1.2150;
629  _voigt_coeffs[0][3] = -1.3509;
630  _voigt_coeffs[1][0] = 1.2359;
631  _voigt_coeffs[1][1] = 0.3786;
632  _voigt_coeffs[1][2] = -1.2359;
633  _voigt_coeffs[1][3] = -0.3786;
634  _voigt_coeffs[2][0] = -0.3085;
635  _voigt_coeffs[2][1] = 0.5906;
636  _voigt_coeffs[2][2] = -0.3085;
637  _voigt_coeffs[2][3] = 0.5906;
638  _voigt_coeffs[3][0] = 0.0210;
639  _voigt_coeffs[3][1] = -1.1858;
640  _voigt_coeffs[3][2] = -0.0210;
641  _voigt_coeffs[3][3] = 1.1858;
642  }
643 
644  template<typename Chemistry>
645  libMesh::Real AbsorptionCoeff<Chemistry>::voigt_a(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
646  {
647  libMesh::Real nu_c = this->nu_C(T,P,Y,i);
648  libMesh::Real nu_D = this->nu_D(T,P,i);
649 
650  return std::sqrt(std::log(2.0))*nu_c/nu_D;
651  }
652 
653  template<typename Chemistry>
654  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_a_dT(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
655  {
656  libMesh::Real nu_c = this->nu_C(T,P,Y,i);
657  libMesh::Real nu_D = this->nu_D(T,P,i);
658  libMesh::Real dnu_c = d_nuC_dT(T,P,Y,i);
659  libMesh::Real dnu_D = d_nuD_dT(T,P,i);
660 
661  return std::sqrt(std::log(2.0)) * (nu_D*dnu_c - nu_c*dnu_D)/(nu_D*nu_D);
662  }
663 
664  template<typename Chemistry>
665  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_a_dP(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
666  {
667  libMesh::Real nu_c = this->nu_C(T,P,Y,i);
668  libMesh::Real nu_D = this->nu_D(T,P,i);
669  libMesh::Real dnu_c = d_nuC_dP(T,Y,i);
670  libMesh::Real dnu_D = d_nuD_dP(T,i);
671 
672  return std::sqrt(std::log(2.0)) * (nu_D*dnu_c - nu_c*dnu_D)/(nu_D*nu_D);
673  }
674 
675  template<typename Chemistry>
676  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_a_dY(libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> Y, unsigned int i)
677  {
678  libMesh::Real nu_D = this->nu_D(T,P,i);
679  libMesh::Real dnu_c = d_nuC_dY(T,P,Y,i);
680 
681  return std::sqrt(std::log(2.0))/nu_D * dnu_c;
682  }
683 
684  template<typename Chemistry>
685  libMesh::Real AbsorptionCoeff<Chemistry>::voigt_w(libMesh::Real T, libMesh::Real P, unsigned int i)
686  {
687  libMesh::Real nu = this->get_nu(P,i);
688  libMesh::Real nu_D = this->nu_D(T,P,i);
689 
690  return 2.0*std::sqrt(std::log(2.0))*(_nu-nu)/nu_D;
691  }
692 
693  template<typename Chemistry>
694  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_w_dT(libMesh::Real T, libMesh::Real P, unsigned int i)
695  {
696  libMesh::Real nu = this->get_nu(P,i);
697  libMesh::Real nu_D = this->nu_D(T,P,i);
698  libMesh::Real dnu_D = d_nuD_dT(T,P,i);
699 
700  return 2.0*std::sqrt(std::log(2.0))*(_nu-nu)/(nu_D*nu_D) * -dnu_D;
701  }
702 
703  template<typename Chemistry>
704  libMesh::Real AbsorptionCoeff<Chemistry>::d_voigt_w_dP(libMesh::Real T, libMesh::Real P, unsigned int i)
705  {
706  libMesh::Real nu = this->get_nu(P,i);
707  libMesh::Real dnu = d_nu_dP(i);
708  libMesh::Real nu_D = this->nu_D(T,P,i);
709  libMesh::Real dnu_D = d_nuD_dP(T,i);
710 
711  return 2.0*std::sqrt(std::log(2.0)) * ( (nu_D)*(-dnu) - (_nu-nu)*(dnu_D) )/(nu_D*nu_D);
712  }
713 
714  template<typename Chemistry>
715  libMesh::Real AbsorptionCoeff<Chemistry>::get_nu(libMesh::Real P, unsigned int i)
716  {
717  libMesh::Real nu = _hitran->nu0(i);
718  libMesh::Real d_air = _hitran->delta_air(i);
719 
720  return nu + d_air*((P/Constants::atmosphere_Pa)/_Pref);
721  }
722 
723  template<typename Chemistry>
724  libMesh::Real AbsorptionCoeff<Chemistry>::d_nu_dP(unsigned int i)
725  {
726  libMesh::Real d_air = _hitran->delta_air(i);
727 
728  return d_air/_Pref * 1.0/Constants::atmosphere_Pa;
729  }
730 
731  template<typename Chemistry>
732  libMesh::Real AbsorptionCoeff<Chemistry>::dX_dY(std::vector<libMesh::Real> Y)
733  {
734  libMesh::Real MW = _chemistry->M(_species_idx);
735  libMesh::Real MW_mix = _chemistry->M_mix(Y);
736 
737  libMesh::Real Ys = Y[_species_idx];
738 
739  return 1.0/MW * ( MW_mix - Ys*(MW_mix*MW_mix)/MW );
740  }
741 
742 
743  template<typename Chemistry>
744  libMesh::Real AbsorptionCoeff<Chemistry>::dQ_dT(libMesh::Real T, unsigned int iso)
745  {
746  libMesh::Real deriv = _hitran->partition_function_derivative(T,iso);
747  return deriv;
748  }
749 
750 #if GRINS_HAVE_ANTIOCH
751  template class AbsorptionCoeff<AntiochChemistry>;
752 #endif
753 
754 #if GRINS_HAVE_CANTERA
755  template class AbsorptionCoeff<CanteraMixture>;
756 #endif
757 } //namespace GRINS
GRINS::AbsorptionCoeff::d_voigt_a_dP
libMesh::Real d_voigt_a_dP(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Voigt a parameter pressure derivative.
Definition: absorption_coeff.C:665
GRINS::AbsorptionCoeff::d_voigt_dY
libMesh::Real d_voigt_dY(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Voigt profile mass fraction derivative.
Definition: absorption_coeff.C:588
GRINS::AbsorptionCoeff::d_nuC_dY
libMesh::Real d_nuC_dY(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Collisional broadening mass fraction derivative.
Definition: absorption_coeff.C:465
cantera_mixture.h
math_constants.h
GRINS::AbsorptionCoeff::operator()
virtual libMesh::Real operator()(const libMesh::FEMContext &context, const libMesh::Point &qp_xyz, const libMesh::Real t)
Calculate the absorption coefficient at a quadratue point.
Definition: absorption_coeff.C:124
GRINS::Constants::pi
const libMesh::Real pi
Definition: math_constants.h:36
GRINS::AbsorptionCoeff::derivatives
virtual void derivatives(libMesh::FEMContext &context, const libMesh::Point &qp_xyz, const libMesh::Real &JxW, const unsigned int qoi_index, const libMesh::Real time)
Calculate the derivative of the absorption coefficient at a QP with respect to all state variables...
Definition: absorption_coeff.C:165
GRINS::AbsorptionCoeff::d_nuD_dT
libMesh::Real d_nuD_dT(libMesh::Real T, libMesh::Real P, unsigned int i)
Doppler broadening temperature derivative.
Definition: absorption_coeff.C:401
GRINS::AbsorptionCoeff::_species_idx
unsigned int _species_idx
Index for the species of interest.
Definition: absorption_coeff.h:136
GRINS::Constants::second_rad_const
const libMesh::Real second_rad_const
Second Radiation Constant hc/k, [m K].
Definition: physical_constants.h:67
GRINS::AbsorptionCoeff::d_voigt_w_dP
libMesh::Real d_voigt_w_dP(libMesh::Real T, libMesh::Real P, unsigned int i)
Voigt w parameter pressure derivative.
Definition: absorption_coeff.C:704
GRINS::AbsorptionCoeff::Sw
libMesh::Real Sw(libMesh::Real T, libMesh::Real P, unsigned int i)
Linestrength [cm^-2 atm^-1].
Definition: absorption_coeff.C:303
GRINS::AbsorptionCoeff::voigt
libMesh::Real voigt(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Calculate the Voigt profile [cm^-1].
Definition: absorption_coeff.C:477
GRINS::AbsorptionCoeff::voigt_a
libMesh::Real voigt_a(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Voigt a parameter.
Definition: absorption_coeff.C:645
GRINS::AbsorptionCoeff::_thermo_pressure
libMesh::Real _thermo_pressure
Thermodynamic Pressure [atm].
Definition: absorption_coeff.h:130
GRINS::AbsorptionCoeff::get_nu
libMesh::Real get_nu(libMesh::Real P, unsigned int i)
Pressure shift of linecenter wavenumber.
Definition: absorption_coeff.C:715
variable_warehouse.h
GRINS::AbsorptionCoeff::d_kv_dP
libMesh::Real d_kv_dP(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Absorption coefficient pressure derivative.
Definition: absorption_coeff.C:274
GRINS::AbsorptionCoeff::d_kv_dY
libMesh::Real d_kv_dY(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Absorption coefficient mass fraction derivative.
Definition: absorption_coeff.C:288
GRINS::AbsorptionCoeff::dS_dT
libMesh::Real dS_dT(libMesh::Real T, libMesh::Real P, unsigned int i)
Linestrength temperature derivative.
Definition: absorption_coeff.C:339
GRINS
GRINS namespace.
Definition: arrhenius_catalycity.h:31
GRINS::AbsorptionCoeff::kv
libMesh::Real kv(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Absorption coefficient [cm^-1].
Definition: absorption_coeff.C:244
GRINS::AbsorptionCoeff::d_voigt_w_dT
libMesh::Real d_voigt_w_dT(libMesh::Real T, libMesh::Real P, unsigned int i)
Voigt w parameter temperature derivative.
Definition: absorption_coeff.C:694
physical_constants.h
GRINS::AbsorptionCoeff::d_nuC_dP
libMesh::Real d_nuC_dP(libMesh::Real T, std::vector< libMesh::Real > Y, unsigned int i)
Collisional broadening pressure derivative.
Definition: absorption_coeff.C:452
GRINS::AbsorptionCoeff::d_nuD_dP
libMesh::Real d_nuD_dP(libMesh::Real T, unsigned int i)
Doppler broadening pressure derivative.
Definition: absorption_coeff.C:413
GRINS::AbsorptionCoeff::d_nu_dP
libMesh::Real d_nu_dP(unsigned int i)
Derivative of pressure-shifted linecenter wavenumber.
Definition: absorption_coeff.C:724
GRINS::AbsorptionCoeff::_max_index
unsigned int _max_index
Index of maximum wavenumber.
Definition: absorption_coeff.h:127
GRINS::AbsorptionCoeff::_calc_thermo_pressure
bool _calc_thermo_pressure
Flag for whether Thermodynamic Pressure is calculated or constant.
Definition: absorption_coeff.h:133
GRINS::AbsorptionCoeff::dX_dY
libMesh::Real dX_dY(std::vector< libMesh::Real > Y)
Mole fraction derivative with respect to mass fraction.
Definition: absorption_coeff.C:732
GRINS::AbsorptionCoeff
Evaluates the Beer-Lambert Law at a given point in space.
Definition: absorption_coeff.h:63
GRINS::Constants::Avogadro
const libMesh::Real Avogadro
Avogadro's number, [particles/mol].
Definition: physical_constants.h:43
GRINS::AbsorptionCoeff::d_nuC_dT
libMesh::Real d_nuC_dT(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Collisional broadening temperature derivative.
Definition: absorption_coeff.C:439
GRINS::AbsorptionCoeff::_chemistry
SharedPtr< Chemistry > _chemistry
Antioch/Cantera object.
Definition: absorption_coeff.h:102
GRINS::AbsorptionCoeff::nu_C
libMesh::Real nu_C(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Collisional broadening [cm^-1].
Definition: absorption_coeff.C:426
GRINS::Constants::c_vacuum
const libMesh::Real c_vacuum
Speed of light in a vacuum, [m/s].
Definition: physical_constants.h:53
GRINS::AbsorptionCoeff::d_voigt_dP
libMesh::Real d_voigt_dP(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Voigt profile pressure derivative.
Definition: absorption_coeff.C:548
GRINS::AbsorptionCoeff::d_voigt_dT
libMesh::Real d_voigt_dT(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Voigt profile temperature derivative.
Definition: absorption_coeff.C:508
GRINS::AbsorptionCoeff::dQ_dT
libMesh::Real dQ_dT(libMesh::Real T, unsigned int iso)
Partition Function derivative (finite difference)
Definition: absorption_coeff.C:744
antioch_chemistry.h
GRINS::Constants::atmosphere_Pa
const libMesh::Real atmosphere_Pa
1 atmosphere in Pascals, [Pa/atm]
Definition: physical_constants.h:58
GRINS::AbsorptionCoeff::_hitran
SharedPtr< HITRAN > _hitran
HITRAN.
Definition: absorption_coeff.h:105
GRINS::AbsorptionCoeff::voigt_w
libMesh::Real voigt_w(libMesh::Real T, libMesh::Real P, unsigned int i)
Voigt w parameter.
Definition: absorption_coeff.C:685
absorption_coeff.h
GRINS::AbsorptionCoeff::d_voigt_a_dY
libMesh::Real d_voigt_a_dY(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Voigt a parameter mass fraction derivative.
Definition: absorption_coeff.C:676
GRINS::AbsorptionCoeff::init_voigt
void init_voigt()
Initialize the coeff matrix for calculating the Voigt profile.
Definition: absorption_coeff.C:621
GRINS::AbsorptionCoeff::AbsorptionCoeff
AbsorptionCoeff()
User should not call empty constructor.
GRINS::AbsorptionCoeff::d_voigt_a_dT
libMesh::Real d_voigt_a_dT(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Voigt a parameter temperature derivative.
Definition: absorption_coeff.C:654
GRINS::AbsorptionCoeff::dS_dP
libMesh::Real dS_dP(libMesh::Real T, libMesh::Real P, unsigned int i)
Linestrength pressure derivative.
Definition: absorption_coeff.C:366
GRINS::AbsorptionCoeff::d_kv_dT
libMesh::Real d_kv_dT(libMesh::Real T, libMesh::Real P, std::vector< libMesh::Real > Y, unsigned int i)
Absorption coefficient temperature derivative.
Definition: absorption_coeff.C:260
GRINS::AbsorptionCoeff::clone
virtual libMesh::UniquePtr< libMesh::FEMFunctionBase< libMesh::Real > > clone() const
Not used.
Definition: absorption_coeff.C:237
GRINS::Constants::Boltzmann
const libMesh::Real Boltzmann
Boltzmann Constant, [J/K].
Definition: physical_constants.h:48
GRINS::AbsorptionCoeff::_min_index
unsigned int _min_index
Index of minimum wavenumber.
Definition: absorption_coeff.h:124
GRINS::AbsorptionCoeff::nu_D
libMesh::Real nu_D(libMesh::Real T, libMesh::Real P, unsigned int i)
Doppler broadening [cm^-1].
Definition: absorption_coeff.C:389
Generated on Tue Dec 19 2017 12:47:28 for GRINS-0.8.0 by  doxygen 1.8.9.1