Added switch to change mixing rule for MW relaxation time. Fixing issue #310

CHANGELOG.md

@@ -11,6 +11,46 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/) and this p
 - `Security` to invite users to upgrade in case of vulnerabilities
 
 ## [Unreleased]
+
+## [1.3.0] - 2026-07-23
+### Added  
+ * `Energy exchange mechanisms`
+    * VT transfer (Millikan and White (MW), Park correction)
+        - Enabling to select original mixing rule or Gnoffo at runtime     
+ * `GSI`
+    * Functions to get surface emissivity and radiative heat flux
+    * Function to set GSI solver convergence tolerance
+ * `Examples`
+    * 2T dissociation of O2
+
+## [1.2.0] - 2026-06-15
+### Changed
+ * `Included Tools`
+    * bprime
+        - Enabling to pass surface composition
+
+## [1.1.3] - 2023-08-23
+### Fixed
+ * `Thermodynamics`:
+    * Pure species thermo: 
+        * Update electron molecular weight
+
+## [1.1.2] - 2022-01-06
+### Fixed
+ * `Interfaces to other languages`
+    * Fix issue with Fortran compiler flags not being set
+
+## [1.1.1] - 2021-09-24
+### Fixed
+ * `Interfaces to other languages`
+    * Fix install path for fortran module and export name
+
+## [1.0.5] - 2021-07-19
+### Added
+ * `Interfaces to other languages`
+    * Python
+
+## [1.0.1] - 2020-04-26
 ### Added
  * `Thermodynamics`:
     * Databases: 
@@ -42,7 +82,7 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/) and this p
     * Surface Energy Balance
 
 * `Energy exchange mechanisms`
-    * VT transfer (Milikan and White (MW), Park correction)
+    * VT transfer (Millikan and White (MW), Park correction)
     * ET transfer
     * Chemical energy exchange:
         * Electron impact ionization and dissociation

docs/contributors.md

@@ -15,7 +15,8 @@ Pierre Schrooyen <br>
 Laurent Soucasse <br>
 Alessandro Turchi <br>
 [Ruben Di Battista](https://rdb.is) <br>
-Matthew Goodson - [Corvid Technologies](https://www.corvidtec.com)
+Matthew Goodson - [Corvid Technologies](https://www.corvidtec.com) <br>
+Michele Capriati
 
 Have you contributed but don't see your name?  Open a PR on Github and add
 yourself in this file, we will be happy to add you.

examples/c++/O2_dissociation2T/CMakeLists.txt

@@ -0,0 +1 @@
+../CMakeLists.txt.shared

examples/c++/O2_dissociation2T/O2_dissociation2T.cpp

@@ -0,0 +1,174 @@
+/**
+ * @file O2_dissociation2T.cpp
+ *
+ * @brief O<sub>2</sub> Dissociation with 2T model example program.
+ */
+
+/*
+ * Copyright 2014-2020 von Karman Institute for Fluid Dynamics (VKI)
+ *
+ * This file is part of MUlticomponent Thermodynamic And Transport
+ * properties for IONized gases in C++ (Mutation++) software package.
+ *
+ * Mutation++ is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * Mutation++ is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Mutation++.  If not, see
+ * <http://www.gnu.org/licenses/>.
+ */
+
+/// [example_code]
+// Must include this header file to use the Mutation++ library
+#include "mutation++.h"
+
+using namespace Mutation;
+using namespace Mutation::Thermodynamics;
+using namespace std;
+
+// Prints the header of the output table
+void printHeader(const Mixture& mix)
+{
+    cout << setw(8) << "Time[s]";
+    cout << setw(12) << "T[K]";
+    cout << setw(12) << "Tv[K]";
+    cout << setw(12) << "P[Pa]";
+    cout << setw(12) << "rho[kg/m^3]";
+    cout << setw(12) << "e[J/kg]";
+    for (int i = 0; i < mix.nSpecies(); ++i)
+        cout << setw(12) << "Y_" + mix.speciesName(i);
+    cout << endl;
+}
+
+// Prints the mixture properties at the given time
+void printResults(double time, const Mixture& mix)
+{
+    cout << setw(8) << time;
+    cout << setw(12) << mix.T();
+    cout << setw(12) << mix.Tv();
+    cout << setw(12) << mix.P();
+    cout << setw(12) << mix.density();
+    cout << setw(12) << mix.mixtureEnergyMass();
+    for (int i = 0; i < mix.nSpecies(); ++i)
+        cout << setw(12) << mix.Y()[i];
+    cout << endl;
+}
+
+// Main entry point
+int main()
+{
+    // Initial conditions are defined here
+    const double T_init   = 10000.0;  // K
+    const double Tv_init  = 300.0;  // K
+    const double rho_init = 0.3899740209756568; // kg/m^3
+
+    // First create the mixture object
+    MixtureOptions opts;
+    opts.setSpeciesDescriptor("O O2");        // 2 species O2 mixture
+    opts.setThermodynamicDatabase("RRHO");    // Thermo database is RRHO
+    opts.setMechanism("O2");                  // O2 + M = 2O + M
+    opts.setStateModel("ChemNonEqTTv");       // chemical noneq. w/ 2T
+    Mixture mix(opts);                        // Init. the mixture with opts
+
+    // Setup arrays
+    double rhoi [mix.nSpecies()];
+    double wdot [mix.nSpecies()];
+    double TToPass [2];
+    double evib;
+    double source;
+    double EnergiesMass[2];
+
+    // Set state of mixture to the initial conditions
+    rhoi[mix.speciesIndex("O")]  = 0.0;
+    rhoi[mix.speciesIndex("O2")] = 1.0*rho_init;
+
+    TToPass[0] = T_init;
+    TToPass[1] = Tv_init;
+    mix.setState(rhoi, TToPass, 1); 
+    mix.mixtureEnergies(EnergiesMass);
+
+    EnergiesMass[0] = EnergiesMass[0]*mix.density();
+    EnergiesMass[1] = EnergiesMass[1]*mix.density();
+
+    // Write the results header and initial conditions
+    printHeader(mix);
+    printResults(0.0, mix);
+
+    //time integration RK4
+    double dt   = 1.0e-10;
+    double tol  = 1.0e-14;
+    double time = 0.0;
+    double k1[mix.nSpecies()],k2[mix.nSpecies()],k3[mix.nSpecies()],k4[mix.nSpecies()], W[mix.nSpecies()], drhoi[mix.nSpecies()];
+    double s1, s2, s3, s4, s[2], dEnergiesMass;
+    double conv = 1.0;
+    s[0] = EnergiesMass[0];
+    s[1] = EnergiesMass[1];
+    while (conv > tol) {
+
+        dt = std::min(dt*1.0002, 1e-5);
+        // Get the species production rates
+        mix.setState(rhoi, EnergiesMass);
+        mix.netProductionRates(k1);
+        mix.energyTransferSource(&s1);
+
+       for (int i = 0; i < mix.nSpecies(); ++i)
+            W[i] = rhoi[i] + 0.5*dt*k1[i];
+	s[1] = EnergiesMass[1] + 0.5*s1*dt;	
+
+        mix.setState(W, s);
+        mix.netProductionRates(k2);
+        mix.energyTransferSource(&s2);
+
+       for (int i = 0; i < mix.nSpecies(); ++i)
+            W[i] = rhoi[i] + 0.5*dt*k2[i];
+	s[1] = EnergiesMass[1] + 0.5*s2*dt;	
+
+        mix.setState(W, s);
+        mix.netProductionRates(k3);
+        mix.energyTransferSource(&s3);
+
+       for (int i = 0; i < mix.nSpecies(); ++i)
+            W[i] = rhoi[i] + dt*k3[i];
+	s[1] = EnergiesMass[1] + s3*dt;	
+
+        mix.setState(W, s);
+        mix.netProductionRates(k4);
+        mix.energyTransferSource(&s4);
+
+        conv = 0.0;
+        for (int i = 0; i < mix.nSpecies(); ++i) {
+            drhoi[i] = (k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]) / 6.0 * dt;
+            rhoi[i] += drhoi[i];
+            if (std::abs(drhoi[i]) > conv) conv = std::abs(drhoi[i]);
+        }
+        dEnergiesMass =  (s1 + 2.0*s2 + 2.0*s3 + s4) / 6.0 * dt;
+	EnergiesMass[1] +=  dEnergiesMass;
+        time += dt;
+
+        // Set the new state of the mixture (using conserved variables)
+        mix.setState(rhoi, EnergiesMass);
+        printResults(time, mix);
+    }
+
+    // Now get the equilibrium values
+    mix.equilibriumComposition(mix.T(), mix.P(), rhoi); // puts X_i in rhoi
+    mix.convert<X_TO_Y>(rhoi, rhoi);                    // converts X_i to Y_i
+
+    cout << "Equilibrium mass fractions at " << mix.T() << " K and "
+         << mix.P() << " Pa:" << endl;
+
+    for (int i = 0; i < mix.nSpecies(); ++i)
+        cout << setw(15) << mix.speciesName(i) << ":"
+             << setw(15) << rhoi[i] << endl;
+
+    return 0;
+}
+/// [example_code]
+

src/transfer/CMakeLists.txt

@@ -21,6 +21,8 @@
 
 add_sources(mutation++
     MillikanWhite.cpp
+    MillikanWhiteRelaxationTimeGnoffo.cpp
+    MillikanWhiteRelaxationTimeOriginal.cpp
     OmegaCE.cpp
     OmegaCElec.cpp
     OmegaCV.cpp
@@ -32,4 +34,5 @@ add_sources(mutation++
 add_headers(mutation++
     MillikanWhite.h
     TransferModel.h
+    MillikanWhiteRelaxationTime.h
 )

src/transfer/MillikanWhite.cpp

@@ -27,8 +27,10 @@
 
 
 #include "MillikanWhite.h"
+#include "MillikanWhiteRelaxationTime.h"
 #include "Thermodynamics.h"
 #include "Utilities.h"
+#include "AutoRegistration.h"
 
 #include <iostream>
 #include <fstream>
@@ -42,6 +44,18 @@ using namespace Mutation::Thermodynamics;
 namespace Mutation {
     namespace Transfer {
 
+static std::string name_model = "Original";
+
+void setMillikanWhiteModel(const std::string& model)
+{
+    name_model = model;
+}
+
+const std::string& getMillikanWhiteModel()
+{
+    return name_model;
+}
+
 
 struct MillikanWhiteModelData::Impl
 {
@@ -50,8 +64,9 @@ struct MillikanWhiteModelData::Impl
     double omegav = 1.0E-20; // m^2
     Eigen::ArrayXd a;
     Eigen::ArrayXd b;
+    Eigen::ArrayXd mu;
 
-    Impl(size_t size) : a(size), b(size) { }
+    Impl(size_t size) : a(size), b(size), mu(size) { }
 };
 
 
@@ -73,6 +88,7 @@ MillikanWhiteModelData::MillikanWhiteModelData(
         double mu = 1000.0*m_impl->mw*mwi/(m_impl->mw + mwi); // reduced mass in g/mol
         m_impl->a[i] = 1.16E-3*std::sqrt(mu)*theta_power;
         m_impl->b[i] = 0.015*std::pow(mu, 0.25);
+        m_impl->mu[i] = mu/1000.;
     }
 }
 
@@ -113,6 +129,10 @@ Eigen::ArrayXd& MillikanWhiteModelData::b() { return m_impl->b; }
 
 const Eigen::ArrayXd& MillikanWhiteModelData::b() const { return m_impl->b; }
 
+Eigen::ArrayXd& MillikanWhiteModelData::mu() { return m_impl->mu; }
+
+const Eigen::ArrayXd& MillikanWhiteModelData::mu() const { return m_impl->mu; }
+
 MillikanWhiteModelData& MillikanWhiteModelData::setReferenceCrossSection(double omegav)
 {
     assert(omegav >= 0.0);
@@ -130,29 +150,19 @@ double MillikanWhiteModelData::limitingCrossSection(const double& T) const
 
 MillikanWhiteModel::MillikanWhiteModel(const MillikanWhiteModelData& data) :
     m_data(data) 
-{ }
-
+{ 
+    m_relaxationTime = std::shared_ptr<MillikanWhiteRelaxationTime>(
+        Config::Factory<MillikanWhiteRelaxationTime>::create(name_model, m_data));
+}
 
 double MillikanWhiteModel::relaxationTime(
     const Mutation::Thermodynamics::Thermodynamics& thermo) const
 {
-    // Millikan-White model for average relaxation time
-    const Eigen::Map<const Eigen::ArrayXd> Xh(
-        thermo.X()+(thermo.hasElectrons() ? 1 : 0), thermo.nHeavy());
-    const double T_fac = std::pow(thermo.T(), -1.0/3.0);
-    const double p_atm = thermo.P() / ONEATM;
-    const double tau_mw = 
-        (Xh*(m_data.a()*(T_fac - m_data.b()) - 18.42).exp()).sum() / 
-        (Xh.sum()*p_atm);
-
-    // Park correction
-    const double ni = thermo.numberDensity() * thermo.X()[m_data.speciesIndex()];
-    const double ci = std::sqrt(8*RU*thermo.T()/(PI*m_data.molecularWeight()));
-    const double tau_park = 1.0/(ni * ci * m_data.limitingCrossSection(thermo.T()));
-
-    return tau_mw + tau_park;
+    return m_relaxationTime->relaxationTime(thermo, m_data);
 }
 
+MillikanWhiteModel::~MillikanWhiteModel() = default;
+
 
 struct MillikanWhiteModelDB::Data
 {