diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ff9ec2b..6c0242f7 100644 --- a/CMakeLists.txt +++ b/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 ) diff --git a/examples/fortran/wrapper_test/wrapper_test.f90 b/examples/fortran/wrapper_test/wrapper_test.f90 index 3a7361f8..458e66cf 100644 --- a/examples/fortran/wrapper_test/wrapper_test.f90 +++ b/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 diff --git a/interface/fortran/cwrapper.cpp b/interface/fortran/cwrapper.cpp index 777a78ab..4c19080e 100644 --- a/interface/fortran/cwrapper.cpp +++ b/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 Yke (ne,0); + std::vector Ykg (ne,0); + std::vector 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) diff --git a/interface/fortran/cwrapper.h b/interface/fortran/cwrapper.h index 0358caa9..287e894c 100644 --- a/interface/fortran/cwrapper.h +++ b/interface/fortran/cwrapper.h @@ -444,13 +444,22 @@ 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 */ @@ -458,7 +467,7 @@ 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, @@ -466,6 +475,15 @@ void NAME_MANGLE(gasmixture_surface_mass_balance) 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. diff --git a/src/apps/bprime.cpp b/src/apps/bprime.cpp index 3de7ec24..73832796 100644 --- a/src/apps/bprime.cpp +++ b/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 Yke(ne, 0); std::vector Ykg(ne, 0); + std::vector Ykc(ne, 0); std::vector 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; diff --git a/src/thermo/Thermodynamics.cpp b/src/thermo/Thermodynamics.cpp index ddad37ac..7febb2c8 100644 --- a/src/thermo/Thermodynamics.cpp +++ b/src/thermo/Thermodynamics.cpp @@ -1090,16 +1090,28 @@ 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 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) { @@ -1107,12 +1119,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 condensed phase to simulate infinite surface + const double LargeNumber = 100.0; + std::vector 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); diff --git a/src/thermo/Thermodynamics.h b/src/thermo/Thermodynamics.h index 9ba1b888..f66f8385 100644 --- a/src/thermo/Thermodynamics.h +++ b/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: