FEAT: pull calculation of lambda

AetherModel · Nov 22, 2024 · c3d04d0 · c3d04d0
1 parent a454d25
commit c3d04d0
Show file tree

Hide file tree

Showing 2 changed files with 70 additions and 2 deletions.
diff --git a/include/ions.h b/include/ions.h
@@ -118,6 +118,9 @@ class Ions {
     /// Specific heat (constant volume):
     arma_cube Cv_scgc;
 
+    /// Heat Conduction:
+    arma_cube lambda;
+
   };
 
   // bulk quantities (states):
@@ -169,6 +172,9 @@ class Ions {
   /// Specific heat (constant volume):
   arma_cube Cv_scgc;
 
+  /// Head Conduction (bulk):
+  arma_cube lambda;
+
   // Electrodynamics:
   /// Electric potential:
   arma_cube potential_scgc;
@@ -260,6 +266,11 @@ class Ions {
    **/
   void calc_specific_heat();
 
+  /**********************************************************************
+     \brief Calculate the individual and bulk thermal conductivities
+   **/
+  void calc_lambda();
+
   /**********************************************************************
      \brief Sets the boundary conditions of the ions
      \param grid The grid that the ions are defined on

diff --git a/src/ions.cpp b/src/ions.cpp
@@ -46,6 +46,8 @@ Ions::species_chars Ions::create_species(Grid grid) {
   tmp.heating_sources_total.zeros();
   tmp.Cv_scgc.set_size(nLons, nLats, nAlts);
   tmp.Cv_scgc.zeros();
+  tmp.lambda.set_size(nLons, nLats, nAlts);
+  tmp.lambda.zeros();
 
   tmp.par_velocity_vcgc = make_cube_vector(nLons, nLats, nAlts, 3);
   tmp.perp_velocity_vcgc = make_cube_vector(nLons, nLats, nAlts, 3);
@@ -107,6 +109,8 @@ Ions::Ions(Grid grid, Planets planet) {
 
   Cv_scgc.set_size(nLons, nLats, nAlts);
   Cv_scgc.zeros();
+  lambda.set_size(nLons, nLats, nAlts);
+  lambda.zeros();
   gamma_scgc.set_size(nLons, nLats, nAlts);
   gamma_scgc.ones();
   sound_scgc.set_size(nLons, nLats, nAlts);
@@ -345,8 +349,61 @@ void Ions::calc_cMax() {
 }
 
 // ----------------------------------------------------------------------
-// Calculate a bunch of derived products:
-//   - Specific Heat at Constant Volume (Cv)
+// Calculate thermal conduction (lambda)
+// ----------------------------------------------------------------------
+
+void Ions::calc_lambda() {
+
+  int64_t iIon, jIon;
+
+  std::string function = "Ions::calc_specific_heat";
+  static int iFunction = -1;
+  report.enter(function, iFunction);
+
+  lambda.zeros();
+  precision_t Mi, Mj;
+  arma_cube density_ratio, ratios;
+
+  // This sets the size:
+  ratios = density_scgc;
+
+  for (iIon = 0; iIon < nSpecies; iIon++) {
+
+    Mi = species[iIon].mass / cAMU;
+    ratios.zeros();
+
+    for (jIon = 0; jIon < nSpecies; jIon++) {
+      if (jIon != iIon) {
+        Mj = species[jIon].mass / cAMU;
+        density_ratio = species[jIon].density_scgc / 
+          species[iIon].density_scgc;
+        density_ratio.clamp(0.001, 1000.0);
+        ratios = ratios + density_ratio * 
+          (species[jIon].charge * species[jIon].charge / 
+           species[iIon].charge / species[iIon].charge) *
+           sqrt(Mj / (Mi + Mj)) *
+          (3 * Mi * Mi + 1.6 * Mi * Mj + 1.3 * Mj * Mj)/
+          ((Mi + Mj) * (Mi + Mj));
+      }
+    species[iIon].lambda = 
+      3.1e6 / sqrt(Mi) / pow(species[iIon].charge, 4) *
+      pow(species[iIon].temperature_scgc, 2.5) % (1 + 1.75 * ratios) * cE;
+    }
+    lambda = lambda + species[iIon].lambda % species[iIon].density_scgc;
+  }
+  lambda = lambda / density_scgc;
+
+  //lambda1d = 25.0 * cKB * pow(temp1d, 2.5) * (cKB / species[iIon].mass)
+  //           / species[iIon].nu_ion_ion[iIon] / 8.0;
+
+
+  report.exit(function);
+  return;
+}
+
+
+// ----------------------------------------------------------------------
+// Calculate Specific Heat at Constant Volume (Cv)
 // ----------------------------------------------------------------------
 
 void Ions::calc_specific_heat() {