Enabling to pass surface composition for B' calculation

CMakeLists.txt

@@ -32,7 +32,7 @@ cmake_policy(SET CMP0048 NEW)
 set(CMAKE_CXX_STANDARD 14)
 
 project(mutation++
-    VERSION 1.1.3
+    VERSION 1.2.0
     LANGUAGES CXX
 )
 

examples/fortran/wrapper_test/wrapper_test.f90

@@ -74,7 +74,7 @@ program main
     do i = 1,295
         T = dble(i-1)*50.0 + 300.0
 
-        call mpp_set_state(T, P, var)
+        call mpp_set_state(P, T, var)
         call mpp_x(species_x)
 
         write(*,'(E12.4)',advance='no') T

interface/fortran/cwrapper.cpp

@@ -455,6 +455,15 @@ void NAME_MANGLE(surface_mass_balance)
     p_mix->surfaceMassBalance(p_Yke, p_Ykg, *T, *P, *Bg, *Bc, *hw, p_Xs);
 }
 
+//==============================================================================
+void NAME_MANGLE(surface_mass_balance_general)
+    (const double *const p_Yke, const double *const p_Ykg,  const double *const p_Ykc,
+     const double* const T, 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, p_Ykc);
+}
+
 //==============================================================================
 void NAME_MANGLE(get_composition)
     (F_STRING mixture, double* const p_Yk, F_STRLEN mixture_length)
@@ -481,6 +490,26 @@ void NAME_MANGLE(gasmixture_surface_mass_balance)
     p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), *T, *P, *Bg, *Bc, *hw, p_Xs);
 }
 
+//==============================================================================
+void NAME_MANGLE(gasmixture_surface_mass_balance_general)
+    (F_STRING edge, F_STRING pyro, F_STRING surf,
+     const double* const T, const double* const P, const double* const Bg,
+     double* const Bc, double* const hw, double *const p_Xs,
+     F_STRLEN edge_length, F_STRLEN pyro_length, F_STRLEN surf_length)
+{
+    const int ne = p_mix->nElements();
+
+    std::vector<double> Yke (ne,0);
+    std::vector<double> Ykg (ne,0);
+    std::vector<double> Ykc (ne,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->getComposition(char_to_string(surf, surf_length), Ykc.data(), Composition::MASS);
+
+    p_mix->surfaceMassBalance(Yke.data(), Ykg.data(), *T, *P, *Bg, *Bc, *hw, p_Xs, Ykc.data());
+}
+
 //==============================================================================
 void NAME_MANGLE(source_energy_transfer)
      (double *const p_source_transfer)

interface/fortran/cwrapper.h

@@ -444,28 +444,46 @@ void NAME_MANGLE(diffusion_matrix)(double* const p_Dij);
 double NAME_MANGLE(sigma)();
 
 /**
- * Solves the surface mass balance at an ablating surface.
+ * Solves the surface mass balance at a carbon ablating surface.
  */
 void NAME_MANGLE(surface_mass_balance)
     (const double *const p_Yke, const double *const p_Ykg, const double* const T,
      const double* const P, const double* const Bg, double* const Bc,
      double* const hw, double *const p_Xs);
 
+/**
+ * Solves the surface mass balance at a general ablating surface.
+ */
+void NAME_MANGLE(surface_mass_balance_general)
+    (const double *const p_Yke, const double *const p_Ykg, const double *const p_Ykc,
+     const double* const T, const double* const P, const double* const Bg,
+     double* const Bc, double* const hw, double *const p_Xs);
+
+
 /**
  * Gets the composition for components of the mixture
  */
 void NAME_MANGLE(get_composition)
     (F_STRING mixture, double* const p_Yk, F_STRLEN mixture_length);
 
 /**
- * Solves the surface mass balance at an ablating surface provided gas mixture names.
+ * Solves the surface mass balance at a carbon ablating surface provided gas mixture names.
  */
 void NAME_MANGLE(gasmixture_surface_mass_balance)
     (F_STRING edge, F_STRING pyro,
      const double* const T, const double* const P, const double* const Bg,
      double* const Bc, double* const hw, double *const p_Xs,
      F_STRLEN edge_length, F_STRLEN pyro_length);
 
+/**
+ * Solves the surface mass balance at a general ablating surface provided gas mixture names.
+ */
+void NAME_MANGLE(gasmixture_surface_mass_balance_general)
+    (F_STRING edge, F_STRING pyro, F_STRING surf,
+     const double* const T, const double* const P, const double* const Bg,
+     double* const Bc, double* const hw, double *const p_Xs,
+     F_STRLEN edge_length, F_STRLEN pyro_length, F_STRLEN surf_length);
+
 /**
  * Returns the pointer to the energy transfer between the internal
  * energy modes.

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 CondensedPhase
  * 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 `CondensedPhase` 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 condensed 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 Options {
@@ -70,8 +69,10 @@ typedef struct Options {
     std::string mixture;
     std::string boundary_layer_comp;
     std::string pyrolysis_composition;
+    std::string condensed_phase_composition;
 
     bool pyrolysis_exist = false;
+    bool condensed_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                 condensed 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 CondensedPhase"
         << 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
+         << "CondensedPhase - corresponds to the condensed 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.condensed_phase_composition = getOption(argc, argv, "-cp");
+        opts.condensed_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> Ykc(ne, 0);
     std::vector<double> Xw(ns, 0);
 
     // Run conditions
@@ -279,6 +293,13 @@ int main(int argc, char* argv[]) {
         mix.getComposition(opts.pyrolysis_composition, Ykg.data(),
                            Composition::MASS);
 
+    if (opts.condensed_phase_exist)
+        mix.getComposition(opts.condensed_phase_composition, Ykc.data(), Composition::MASS);
+    else {
+        int ic = mix.elementIndex("C");
+        Ykc[ic] = 1.0;
+    }
+
     cout << setw(10) << "\"Tw[K]\"" << setw(15) << "\"B'c\"" << setw(15)
          << "\"hw[MJ/kg]\"";
     for (int i = 0; i < ns; ++i)
@@ -287,7 +308,7 @@ int main(int argc, char* argv[]) {
 
     for (double T = T1; T < T2 + 1.0e-6; T += dt) {
         mix.surfaceMassBalance(Yke.data(), Ykg.data(), T, P, Bg, Bc, hw,
-                               Xw.data());
+                               Xw.data(), Ykc.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];
         cout << endl;

src/thermo/Thermodynamics.cpp

@@ -1090,29 +1090,45 @@ 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 T, const double P, const double Bg, double &Bc, double &hw,
+    double *const p_Xs, const double *const p_Ykc)
 {
     const int ne = nElements();
     const int ng = nGas();
 
     double p_Xw [ne];
-    double* p_X  = (p_Xs != NULL ? p_Xs : mp_work1);
+    double* p_X  = (p_Xs != nullptr ? p_Xs : mp_work1);
     double* p_h  = mp_work2;
 
+    const double *p_Ykc_to_use = p_Ykc;
+    std::vector<double> default_Ykc(ne);
+    if (p_Ykc == nullptr) {
+        for (int i = 0; i < ne; ++i) {
+            default_Ykc[i] = 0.0;
+        }    
+        int ic = elementIndex("C");
+        default_Ykc[ic] = 1.0;
+        p_Ykc_to_use = default_Ykc.data();
+    }
+
     // Initialize the wall element fractions to be the pyrolysis gas fractions
     double sum = 0.0;
     for (int i = 0; i < ne; ++i) {
         p_Xw[i] = p_Yke[i] + Bg*p_Ykg[i];
         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 condensed phase to simulate infinite surface
+    const double LargeNumber = 100.0; 
+    std::vector<int> condensedPhaseElements;
+    const double tol = 1.0e-16;
+    for (int i = 0; i < ne; ++i) {
+        p_Xw[i] += LargeNumber*p_Ykc_to_use[i];
+        sum += LargeNumber*p_Ykc_to_use[i];
+        if (abs(p_Ykc_to_use[i]) > tol)
+            condensedPhaseElements.push_back(i);
+    }
 
     for (int i = 0; i < ne; ++i)
         p_Xw[i] /= sum;
@@ -1123,22 +1139,35 @@ 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 i = 0; i < ne; ++i) {
+        p_Yw[i] = 0.0;
+    }
 
     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 = condensedPhaseElements.size();
+    for (int i=0; i < ncp; ++i ) {
+        int idx_cp = condensedPhaseElements[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-15);
 
     // Compute the gas enthalpy
     speciesHOverRT(T, p_h);

src/thermo/Thermodynamics.h

@@ -880,11 +880,13 @@ class Thermodynamics //: public StateModelUpdateHandler
      * @param hw On return, enthalpy of the mixed gas at the surface in J/kg.
      * @param p_Xs (optional) On return, vector of species mole fractions at
      * the surface.
+     * @param p_Ykc (optional) Element mass fractions of the condensed phase.
+     * If not value is passed, a carbon surface is assumed.
      */
     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,
-        double *const p_Xs = NULL);
+        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 = nullptr, const double *const p_Ykc = nullptr);
 
 private: