C3M/src/CustomTransport.cpp at main · GalacticBobster/C3M · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
// Author: Ananyo Bhattacharya
// Affiliation: University of Michigan
// Email: ananyo@umich.edu

// C/C++ headers
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

// Cantera headers
#include <cantera/base/Solution.h>
#include <cantera/base/ct_defs.h>
#include <cantera/kinetics/Kinetics.h>
#include <cantera/kinetics/MultiRate.h>
#include <cantera/kinetics/ReactionData.h>
#include <cantera/kinetics/ReactionRate.h>
#include <cantera/thermo.h>

// C3M headers
#include "CustomTransport.hpp"

using namespace std;
using Eigen::MatrixXd;
using Eigen::VectorXd;

VectorXd handleCustomMolecularDiffusion(string PlanetName,
                                        Cantera::ThermoPhase* NetworkName,
                                        double Pres, double Temp,
                                        VectorXd mWt) {
  int nsp = NetworkName->nSpecies();
  VectorXd DiffConst = VectorXd::Zero(nsp);

  if (PlanetName == "JupiterAurora") {
    DiffConst = JupiterMolDiff(NetworkName, Pres, Temp, mWt);
  }
  return DiffConst;
}

VectorXd JupiterMolDiff(Cantera::ThermoPhase* NetworkName, double Pres,
                        double Temp, VectorXd mWt) {
  int nsp = NetworkName->nSpecies();
  VectorXd DiffConst = VectorXd::Zero(nsp);
  VectorXd col_freq = VectorXd::Zero(nsp);
  VectorXd mol_fr = VectorXd::Zero(nsp);
  VectorXd I = VectorXd::Ones(nsp);
  int elementIndex = 0;  // NetworkName->elementIndex("H2");
  NetworkName->getMoleFractions(&mol_fr[0]);
  double r = 2.7E-10;  // Collision radius of H2 (m)
  //Get total number density of H2
 double totalDensity = NetworkName->molarDensity(); //kmol/m^3
 double nH2 = totalDensity*1E3*6.022E23*mol_fr(elementIndex); //kmol/m^3 -> #/m^3
//Compute collision frequency
 col_freq = 6.022e23*2*3.14*1.38E-23*Temp*(mWt + (I*mWt(elementIndex)))/(mWt(elementIndex)*1E-3);
 col_freq = (col_freq.array()/mWt.array()).sqrt().matrix();
//Compute binary diffusion coefficient
 col_freq = col_freq*2*nH2*r*r;
 DiffConst = (1/col_freq.array()).matrix();
 DiffConst = (DiffConst.array()/mWt.array()).matrix()*1E3*1.38E-23*6.022e23*Temp;
// DiffConst(1) = 3.13E-5*7.339E11*pow(Temp, 0.765);
 DiffConst(0) = 0.0;
// std::cout << "Diffusion constant" << std::endl;
 std::cout << DiffConst.transpose() << std::endl;
 return DiffConst;
}