Enabling to pass surface composition for B' calculation

interface/fortran/cwrapper.cpp

@@ -452,7 +452,11 @@ void NAME_MANGLE(surface_mass_balance)
      const double* const P, const double* const Bg, double* const Bc,
      double* const hw, double *const p_Xs)
 {
-    p_mix->surfaceMassBalance(p_Yke, p_Ykg, *T, *P, *Bg, *Bc, *hw, p_Xs);
+    const int ne = p_mix->nElements();
+    std::vector<double> Ykc (ne,0);
+    int ic = p_mix->elementIndex("C");
+    Ykc[ic] = 1.0; 
+    p_mix->surfaceMassBalance(p_Yke, p_Ykg, Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs);
 }
 
 //==============================================================================
@@ -474,11 +478,14 @@ void NAME_MANGLE(gasmixture_surface_mass_balance)
 
     std::vector<double> Yke (ne,0);
     std::vector<double> Ykg (ne,0);
+    std::vector<double> Ykc (ne,0);
+    int ic = p_mix->elementIndex("C");
+    Ykc[ic] = 1.0; 
 
     p_mix->getComposition(char_to_string(edge, edge_length), Yke.data(), Composition::MASS);
     p_mix->getComposition(char_to_string(pyro, pyro_length), Ykg.data(), Composition::MASS);
 
-    p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), *T, *P, *Bg, *Bc, *hw, p_Xs);
+    p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), Ykc.data(), *T, *P, *Bg, *Bc, *hw, p_Xs);
 }
 
 //==============================================================================

src/apps/bprime.cpp

@@ -46,17 +46,16 @@ using namespace Mutation::Utilities;
  *
  *    bprime \f$-T\f$ \f$T_1\f$:\f$\Delta T\f$:\f$T_2\f$ \f$-p\f
  *    \f$p\f$ \f$-b\f \f$B'_g\f$ \f$-m\f mixture \f$-bl\f BL \f$-py \f Pyrolysis
- *
+ *    \f$-cp \f CondencedPhase
  * This program generates a so-called "B-prime" table for a given temperature
  * range and stepsize in K, a fixed pressure in Pa, a value of \f$B'_g\f$
  * (pyrolysis mass flux, nondimensionalized by the boundary layer edge mass
- * flux), a given mixture name.  Currently, this program assumes the char
- * composition to be solid graphite (which must be included in the mixture
- * file). The `BL` and `Pyrolysis` arguments are the names of the [elemental
- * compositions](@ref compositions) in the mixture file which represent the
- * boundary layer edge and pyrolysis gases respectively. The produced table
- * provides values of \f$B'_c\f$, the wall enthalpy in MJ/kg, and the species
- * mole fractions at the wall versus temperature.
+ * flux), a given mixture name.  
+ * The `BL`, `Pyrolysis`, and `CondencedPhase` arguments are the names of the 
+ * [elemental compositions](@ref compositions) in the mixture file which represent 
+ * the boundary layer edge and pyrolysis gases, and the condenced phase. 
+ * The produced table provides values of \f$B'_c\f$, the wall enthalpy in MJ/kg, 
+ * and the species mole fractions at the wall versus temperature.
  */
 // Simply stores the command line options
 typedef struct {
@@ -70,8 +69,10 @@ typedef struct {
     std::string mixture;
     std::string boundary_layer_comp;
     std::string pyrolysis_composition;
+    std::string condenced_phase_composition;
 
     bool pyrolysis_exist = false;
+    bool condenced_phase_exist = false;
 } Options;
 
 // Checks if an option is present
@@ -122,12 +123,15 @@ void printHelpMessage(const char* const name) {
     cout << tab
          << "-py                 pyrolysis composition name (default = null)"
          << endl;
+    cout << tab
+         << "-cp                 condenced phase composition name (default = carbon)"
+         << endl;
 
     cout << endl;
     cout << "Example:" << endl;
     cout
         << tab << name
-        << " -T 300:100:5000 -P 101325 -b 10 -m carbonPhenol -bl BLedge -py Gas"
+        << " -T 300:100:5000 -P 101325 -b 10 -m carbonPhenol -bl BLedge -py Gas -cp CondencedPhase"
         << endl;
     cout << endl;
     cout << "Mixture file:" << endl;
@@ -140,6 +144,9 @@ void printHelpMessage(const char* const name) {
     cout << tab
          << "Gas - corresponds to the pyrolysis elemental gas composition"
          << endl;
+    cout << tab
+         << "CondencedPhase - corresponds to the condenced phase composition"
+         << endl;
     cout << endl;
 
     exit(0);
@@ -246,6 +253,12 @@ Options parseOptions(int argc, char** argv) {
         opts.pyrolysis_exist = true;
     }
 
+    if (optionExists(argc, argv, "-cp")) {
+        opts.condenced_phase_composition = getOption(argc, argv, "-cp");
+        opts.condenced_phase_exist = true;
+    }
+
+
     return opts;
 }
 
@@ -263,6 +276,7 @@ int main(int argc, char* argv[]) {
 
     std::vector<double> Yke(ne, 0);
     std::vector<double> Ykg(ne, 0);
+    std::vector<double> Ycp(ne, 0);
     std::vector<double> Xw(ns, 0);
 
     // Run conditions
@@ -279,14 +293,21 @@ int main(int argc, char* argv[]) {
         mix.getComposition(opts.pyrolysis_composition, Ykg.data(),
                            Composition::MASS);
 
+    if (opts.condenced_phase_exist)
+        mix.getComposition(opts.condenced_phase_composition, Ycp.data(), Composition::MASS);
+    else {
+        int ic = mix.elementIndex("C");
+        Ycp[ic] = 1.0; 
+    }
+
     cout << setw(10) << "\"Tw[K]\"" << setw(15) << "\"B'c\"" << setw(15)
          << "\"hw[MJ/kg]\"";
     for (int i = 0; i < ns; ++i)
         cout << setw(25) << "\"" + mix.speciesName(i) + "\"";
     cout << endl;
 
     for (double T = T1; T < T2 + 1.0e-6; T += dt) {
-        mix.surfaceMassBalance(Yke.data(), Ykg.data(), T, P, Bg, Bc, hw,
+        mix.surfaceMassBalance(Yke.data(), Ykg.data(), Ycp.data(), T, P, Bg, Bc, hw,
                                Xw.data());
         cout << setw(10) << T << setw(15) << Bc << setw(15) << hw / 1.0e6;
         for (int i = 0; i < ns; ++i) cout << setw(25) << Xw[i];

src/thermo/Thermodynamics.cpp

@@ -1090,8 +1090,9 @@ void Thermodynamics::elementFractions(
 //==============================================================================
 
 void Thermodynamics::surfaceMassBalance(
-    const double *const p_Yke, const double *const p_Ykg, const double T, 
-    const double P, const double Bg, double &Bc, double &hw, double *const p_Xs)
+    const double *const p_Yke, const double *const p_Ykg, const double *const p_Ycp, 
+    const double T, const double P, const double Bg, double &Bc, double &hw,
+    double *const p_Xs)
 {
     const int ne = nElements();
     const int ng = nGas();
@@ -1107,12 +1108,16 @@ void Thermodynamics::surfaceMassBalance(
         sum += p_Xw[i];
     }
 
-    // Use "large" amount of carbon to simulate infinite char
-    int ic = elementIndex("C");
-    //double carbon = std::min(1000.0, std::max(100.0,1000.0*Bg));
-    double carbon = std::max(100.0*Bg, 200.0);
-    p_Xw[ic] += carbon;
-    sum += carbon;
+    // Use "large" amount of condences phase to simulate infinite surface
+    double LargeNumber = 100.0; 
+    std::vector<int> condencedPhaseElements;
+    double tol = 1.0e-16;
+    for (int i = 0; i < ne; ++i) {
+        p_Xw[i] += LargeNumber*p_Ycp[i];
+        sum += LargeNumber*p_Ycp[i];
+        if (abs(p_Ycp[i]) > tol) 
+            condencedPhaseElements.push_back(i);   
+    }
 
     for (int i = 0; i < ne; ++i)
         p_Xw[i] /= sum;
@@ -1123,22 +1128,32 @@ void Thermodynamics::surfaceMassBalance(
 
     // Compute the gas mass fractions at the wall
     double mwg = 0.0;
-    double ywc = 0.0;
+    double p_Yw[ne];
 
     for (int j = 0; j < ng; ++j) {
         mwg += speciesMw(j) * p_X[j];
-        ywc += elementMatrix()(j,ic) * p_X[j];
-        //for (int i = 0; i < ne; ++i)
-        //    p_Yw[i] += elementMatrix()(j,i) * p_X[j];
+        for (int i = 0; i < ne; ++i)
+            p_Yw[i] += elementMatrix()(j,i) * p_X[j];
     }
 
-    //for (int i = 0; i < ne; ++i)
-    //    p_Yw[i] *= atomicMass(i) / mwg;
-    ywc *= atomicMass(ic) / mwg;
+    for (int i = 0; i < ne; ++i)
+        p_Yw[i] *= atomicMass(i) / mwg;
+
+    double sum_Ye = 0.0;
+    double sum_Yg = 0.0;
+    double sum_Yw = 0.0;
+    double sum_YCp = 1.0;  //Assumed equal to one
+    int ncp = condencedPhaseElements.size();
+    for (int i=0; i < ncp; ++i ) {
+        int idx_cp = condencedPhaseElements[i];
+        sum_Ye += p_Yke[idx_cp];
+        sum_Yg += p_Ykg[idx_cp];
+        sum_Yw += p_Yw[idx_cp];
+    }    
 
     // Compute char mass blowing rate
-    Bc = (p_Yke[ic] + Bg*p_Ykg[ic] - ywc*(1.0 + Bg)) / (ywc - 1.0);
-    Bc = std::max(Bc, 0.0);
+    Bc = (Bg*(sum_Yg - sum_Yw) + sum_Ye - sum_Yw)/(sum_Yw - sum_YCp);
+    Bc = std::max(Bc, 1.0e-8);
 
     // Compute the gas enthalpy
     speciesHOverRT(T, p_h);

src/thermo/Thermodynamics.h

@@ -873,6 +873,7 @@ class Thermodynamics //: public StateModelUpdateHandler
      * 
      * @param p_Yke Element mass fractions of boundary layer edge.
      * @param p_Ykg Element mass fractions of the pyrolysis gas.
+     * @param p_Ycp Element mass fractions of the condenced phase.
      * @param T Temperature at the surface in K.
      * @param P Pressure at the surface in Pa.
      * @param Bg Non-dimensional pyrolysis gas mass blowing rate.
@@ -882,8 +883,8 @@ class Thermodynamics //: public StateModelUpdateHandler
      * the surface.
      */
     void surfaceMassBalance(
-        const double *const p_Yke, const double *const p_Ykg, const double T, 
-        const double P, const double Bg, double &Bc, double &hw,
+        const double *const p_Yke, const double *const p_Ykg, const double *const p_Ycp, 
+        const double T, const double P, const double Bg, double &Bc, double &hw,
         double *const p_Xs = NULL);
 
 private:

thirdparty/pybind11/CMakeLists.txt

@@ -5,7 +5,7 @@
 # All rights reserved. Use of this source code is governed by a
 # BSD-style license that can be found in the LICENSE file.
 
-cmake_minimum_required(VERSION 3.4)
+cmake_minimum_required(VERSION 3.5)
 
 # The `cmake_minimum_required(VERSION 3.4...3.22)` syntax does not work with
 # some versions of VS that have a patched CMake 3.11. This forces us to emulate