diff --git a/.gitignore b/.gitignore
index bad0f74f..b9725ab2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,3 +19,5 @@ docs/html/
# OS specifics
.DS_Store
+
+*/CMakeFiles/
diff --git a/data/gsi/smb_ACA_model.xml b/data/gsi/smb_ACA_model.xml
new file mode 100644
index 00000000..b490d2e0
--- /dev/null
+++ b/data/gsi/smb_ACA_model.xml
@@ -0,0 +1,111 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/data/thermo/elements.xml b/data/thermo/elements.xml
index 186801fb..562306b0 100644
--- a/data/thermo/elements.xml
+++ b/data/thermo/elements.xml
@@ -2,7 +2,7 @@
- 0.00055
+ 0.000548579903
@@ -62,7 +62,7 @@
- 12.011
+ 12.0107
diff --git a/src/fortran/cwrapper.cpp b/src/fortran/cwrapper.cpp
index 004d2383..ae48717a 100644
--- a/src/fortran/cwrapper.cpp
+++ b/src/fortran/cwrapper.cpp
@@ -141,6 +141,12 @@ void NAME_MANGLE(species_name)(int* index, F_STRING species, F_STRLEN species_le
string_to_char(p_mix->speciesName(*index-1), species, species_length);
}
+//==============================================================================
+bool NAME_MANGLE(has_electrons)()
+{
+ return p_mix->hasElectrons();
+}
+
//==============================================================================
double NAME_MANGLE(mixture_mw)()
{
@@ -366,6 +372,24 @@ double NAME_MANGLE(heavy_thermal_conductivity)()
return p_mix->heavyThermalConductivity();
}
+//==============================================================================
+double NAME_MANGLE(rotational_thermal_conductivity)()
+{
+ return p_mix->rotationalThermalConductivity();
+}
+
+//==============================================================================
+double NAME_MANGLE(vibrational_thermal_conductivity)()
+{
+ return p_mix->vibrationalThermalConductivity();
+}
+
+//==============================================================================
+double NAME_MANGLE(electronic_thermal_conductivity)()
+{
+ return p_mix->electronicThermalConductivity();
+}
+
//==============================================================================
double NAME_MANGLE(electron_thermal_conductivity)()
{
@@ -373,9 +397,9 @@ double NAME_MANGLE(electron_thermal_conductivity)()
}
//==============================================================================
-double NAME_MANGLE(internal_thermal_conductivity)(double T)
+double NAME_MANGLE(internal_thermal_conductivity)(double *T)
{
- return p_mix->internalThermalConductivity(T);
+ return p_mix->internalThermalConductivity(*T);
}
//==============================================================================
@@ -483,4 +507,17 @@ void NAME_MANGLE(convert_ys_to_ye)(
p_mix->convert(species_y, elements_y);
}
-//==============================================================================
\ No newline at end of file
+//==============================================================================
+double NAME_MANGLE(mixture_s_mass)()
+{
+ return p_mix->mixtureSMass();
+}
+
+//==============================================================================
+
+double NAME_MANGLE(mixture_equilibrium_sound_speed)()
+{
+ return p_mix->equilibriumSoundSpeed();
+}
+
+//==============================================================================
diff --git a/src/fortran/cwrapper.h b/src/fortran/cwrapper.h
index 22c4411e..5961c809 100644
--- a/src/fortran/cwrapper.h
+++ b/src/fortran/cwrapper.h
@@ -126,6 +126,11 @@ int NAME_MANGLE(species_index)(
void NAME_MANGLE(species_name)(
int* index, F_STRING species, F_STRLEN species_length);
+/**
+ * Returns true if the gas has electrons.
+ */
+bool NAME_MANGLE(has_electrons)();
+
/**
* Returns the mixture molecular weight in kg/mol.
*/
@@ -360,6 +365,24 @@ double NAME_MANGLE(equilibrium_thermal_conductivity)();
*/
double NAME_MANGLE(heavy_thermal_conductivity)();
+/**
+ * Returns the rotational thermal conductivity using the
+ * set algorithm.
+ */
+double NAME_MANGLE(rotational_thermal_conductivity)();
+
+/**
+ * Returns the vibrational thermal conductivity using the
+ * set algorithm.
+ */
+double NAME_MANGLE(vibrational_thermal_conductivity)();
+
+/**
+ * Returns the electronic thermal conductivity using the
+ * set algorithm.
+ */
+double NAME_MANGLE(electronic_thermal_conductivity)();
+
/**
* Returns the electron translational thermal conductivity.
*/
@@ -368,7 +391,7 @@ double NAME_MANGLE(electron_thermal_conductivity)();
/**
* Returns the internal energy thermal conductivity using Euken's formulas.
*/
-double NAME_MANGLE(internal_thermal_conductivity)(double T);
+double NAME_MANGLE(internal_thermal_conductivity)(double *T);
/**
* Returns the reactive thermal conductivity which accounts for reactions
@@ -483,6 +506,16 @@ void NAME_MANGLE(convert_ye_to_xe)(
void NAME_MANGLE(convert_ys_to_ye)(
const double* species_y, double* elements_y);
+/**
+* Returns the equilibrium sound speed
+*/
+double NAME_MANGLE(mixture_equilibrium_sound_speed)();
+
+/**
+* Returns the entropy per unit mass
+*/
+double NAME_MANGLE(mixture_s_mass)();
+
#ifdef __cplusplus
}
#endif
diff --git a/src/fortran/cwrapper_interface.f90 b/src/fortran/cwrapper_interface.f90
index 8bce48da..844fc307 100644
--- a/src/fortran/cwrapper_interface.f90
+++ b/src/fortran/cwrapper_interface.f90
@@ -105,7 +105,8 @@ real(kind=8) function mpp_heavy_thermal_conductivity()
real(kind=8) function mpp_electron_thermal_conductivity()
end function
- real(kind=8) function mpp_internal_thermal_conductivity()
+ real(kind=8) function mpp_internal_thermal_conductivity(Tv)
+ real(kind=8), intent(in) :: Tv
end function
real(kind=8) function mpp_reactive_thermal_conductivity()
@@ -113,6 +114,24 @@ real(kind=8) function mpp_reactive_thermal_conductivity()
real(kind=8) function mpp_sigma()
end function
+
+ logical function mpp_has_electrons()
+ end function
+
+ real(kind=8) function mpp_rotational_thermal_conductivity()
+ end function
+
+ real(kind=8) function mpp_vibrational_thermal_conductivity()
+ end function
+
+ real(kind=8) function mpp_electronic_thermal_conductivity()
+ end function
+
+ real(kind=8) function mpp_mixture_s_mass()
+ end function
+
+ real(kind=8) function mpp_mixture_equilibrium_sound_speed()
+ end function
end interface
diff --git a/src/gsi/CMakeLists.txt b/src/gsi/CMakeLists.txt
index b483766f..f2a864b8 100644
--- a/src/gsi/CMakeLists.txt
+++ b/src/gsi/CMakeLists.txt
@@ -33,6 +33,9 @@ add_sources(mutation++
GSIRateLawGammaConst.cpp
GSIRateLawGammaT.cpp
GSIRateLawERArrhenius.cpp
+ GSIRateLawER2Arrhenius.cpp
+ GSIRateLawHLArrhenius.cpp
+ GSIRateLawHLArrhenius_abla.cpp
GSIRateLawSublimation.cpp
GSIRateManagerDetailed.cpp
GSIRateManagerPhenomenological.cpp
diff --git a/src/gsi/GSIRateLawER2Arrhenius.cpp b/src/gsi/GSIRateLawER2Arrhenius.cpp
new file mode 100644
index 00000000..61cc9b56
--- /dev/null
+++ b/src/gsi/GSIRateLawER2Arrhenius.cpp
@@ -0,0 +1,117 @@
+/**
+ * @file GSIRateLawGammaT.cpp
+ *
+ * @brief Class which computes the reaction rate constant for an ER
+ * surface reaction based on an Arrhenius formula.
+ */
+
+/*
+ * Copyright 2018 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
+ * .
+ */
+
+
+#include "Thermodynamics.h"
+#include "Transport.h"
+
+#include "AutoRegistration.h"
+#include "Utilities.h"
+
+#include "GSIRateLaw.h"
+#include "SurfaceProperties.h"
+
+using namespace Mutation::Utilities::Config;
+
+namespace Mutation {
+ namespace GasSurfaceInteraction {
+
+class GSIRateLawER2Arrhenius : public GSIRateLaw
+{
+public:
+ GSIRateLawER2Arrhenius(ARGS args)
+ : GSIRateLaw(args),
+ m_surf_props(args.s_surf_props),
+ mv_react(args.s_reactants),
+ pos_T_trans(0),
+ pos_gas_r(0),
+ pos_site_r(1)
+ {
+ assert(args.s_node_rate_law.tag() == "ER_arrhenius");
+
+ args.s_node_rate_law.getAttribute( "pre_exp", m_pre_exp,
+ "The sticking coeffcient probability for the reaction "
+ "should be provided.");
+ args.s_node_rate_law.getAttribute( "T", m_T_act,
+ "The activation temperature for the reaction "
+ "should be provided for an adsorption reaction.");
+
+ // For the gas in the reactants
+ m_idx_gas = mv_react[pos_gas_r];
+ // Error if m_idx_gas > ns
+
+ // For the sites
+ int idx_site = mv_react[pos_site_r];
+ // Error if idx_site > ns
+
+ m_site_categ = m_surf_props.siteSpeciesToSiteCategoryIndex(idx_site);
+ m_n_sites = m_surf_props.nSiteDensityInCategory(m_site_categ);
+ }
+
+//==============================================================================
+
+ ~GSIRateLawER2Arrhenius( ){ }
+
+//==============================================================================
+
+ double forwardReactionRateCoefficient(
+ const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_Tsurf) const
+ {
+ const double Tsurf = v_Tsurf(pos_T_trans);
+
+ const int set_state_with_rhoi_T = 1;
+ m_thermo.setState(v_rhoi.data(), v_Tsurf.data(), set_state_with_rhoi_T);
+ const double thermal_speed =
+ m_transport.speciesThermalSpeed(m_idx_gas);
+
+ return m_pre_exp * thermal_speed /
+ (4. * m_n_sites * m_n_sites) * exp(-m_T_act / Tsurf);
+ }
+
+private:
+ const size_t pos_T_trans;
+ const size_t pos_gas_r;
+ const size_t pos_site_r;
+
+ int m_idx_gas;
+ int m_site_categ;
+ double m_n_sites;
+
+ double m_pre_exp;
+ double m_T_act;
+
+ const std::vector& mv_react;
+ const SurfaceProperties& m_surf_props;
+};
+
+ObjectProvider<
+ GSIRateLawER2Arrhenius, GSIRateLaw>
+ gsi_rate_law_er2_arrhenius("ER_2_arrhenius");
+
+ } // namespace GasSurfaceInteraction
+} // namespace Mutation
diff --git a/src/gsi/GSIRateLawERArrhenius.cpp b/src/gsi/GSIRateLawERArrhenius.cpp
index f4ff254a..4fbf1e78 100644
--- a/src/gsi/GSIRateLawERArrhenius.cpp
+++ b/src/gsi/GSIRateLawERArrhenius.cpp
@@ -1,7 +1,7 @@
/**
* @file GSIRateLawGammaT.cpp
*
- * @brief Class which computes the reaction rate constant for an adsorption
+ * @brief Class which computes the reaction rate constant for an ER
* surface reaction based on an Arrhenius formula.
*/
diff --git a/src/gsi/GSIRateLawHLArrhenius.cpp b/src/gsi/GSIRateLawHLArrhenius.cpp
new file mode 100644
index 00000000..35ce40cf
--- /dev/null
+++ b/src/gsi/GSIRateLawHLArrhenius.cpp
@@ -0,0 +1,110 @@
+/**
+ * @file GSIRateLawGammaT.cpp
+ *
+ * @brief Class which computes the reaction rate constant for an HL
+ * surface reaction based on an Arrhenius formula.
+ */
+
+/*
+ * Copyright 2018 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
+ * .
+ */
+
+
+#include "Thermodynamics.h"
+#include "Transport.h"
+
+#include "AutoRegistration.h"
+#include "Utilities.h"
+
+#include "GSIRateLaw.h"
+#include "SurfaceProperties.h"
+
+using namespace Mutation::Utilities::Config;
+
+namespace Mutation {
+ namespace GasSurfaceInteraction {
+
+class GSIRateLawLHArrhenius : public GSIRateLaw
+{
+public:
+ GSIRateLawLHArrhenius(ARGS args)
+ : GSIRateLaw(args),
+ m_surf_props(args.s_surf_props),
+ mv_react(args.s_reactants),
+ pos_T_trans(0),
+ pos_site_r(0)
+ {
+ assert(args.s_node_rate_law.tag() == "LH_arrhenius");
+
+ args.s_node_rate_law.getAttribute( "pre_exp", m_pre_exp,
+ "The sticking coeffcient probability for the reaction "
+ "should be provided.");
+ args.s_node_rate_law.getAttribute( "T", m_T_act,
+ "The activation temperature for the reaction "
+ "should be provided for an adsorption reaction.");
+
+ // For the sites
+ m_idx_site = mv_react[pos_site_r];
+
+ m_site_categ = m_surf_props.siteSpeciesToSiteCategoryIndex(m_idx_site);
+ m_n_sites = m_surf_props.nSiteDensityInCategory(m_site_categ);
+
+ int gas_index = m_surf_props.surfaceToGasIndex(m_idx_site);
+ m_mw_r = m_thermo.speciesMw()[gas_index] / NA;
+ }
+
+//==============================================================================
+
+ ~GSIRateLawLHArrhenius( ){ }
+
+//==============================================================================
+
+ double forwardReactionRateCoefficient(
+ const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_Tsurf) const
+ {
+ const double Tsurf = v_Tsurf(pos_T_trans);
+
+ return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*Tsurf/ (2.0*m_mw_r)) * exp(-m_T_act/ Tsurf);
+
+ }
+
+private:
+ const size_t pos_T_trans;
+ const size_t pos_site_r;
+
+ int m_idx_site;
+ int m_site_categ;
+ double m_n_sites;
+
+ double m_pre_exp;
+ double m_T_act;
+
+ double m_mw_r;
+
+ const std::vector& mv_react;
+ const SurfaceProperties& m_surf_props;
+};
+
+ObjectProvider<
+ GSIRateLawLHArrhenius, GSIRateLaw>
+ gsi_rate_law_lh_arrhenius("LH_arrhenius");
+
+ } // namespace GasSurfaceInteraction
+} // namespace Mutation
diff --git a/src/gsi/GSIRateLawHLArrhenius_abla.cpp b/src/gsi/GSIRateLawHLArrhenius_abla.cpp
new file mode 100644
index 00000000..cf5c3c47
--- /dev/null
+++ b/src/gsi/GSIRateLawHLArrhenius_abla.cpp
@@ -0,0 +1,87 @@
+/**
+ * @file GSIRateLawGammaT.cpp
+ *
+ * @brief Class which computes the reaction rate constant for a gas
+ * independent ablation surface reaction based on an Arrhenius formula.
+ */
+
+/*
+ * Copyright 2018 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
+ * .
+ */
+
+
+#include "Thermodynamics.h"
+#include "Transport.h"
+
+#include "AutoRegistration.h"
+#include "Utilities.h"
+
+#include "GSIRateLaw.h"
+#include "SurfaceProperties.h"
+
+using namespace Mutation::Utilities::Config;
+
+namespace Mutation {
+ namespace GasSurfaceInteraction {
+
+class GSIRateLawLHAblaArrhenius : public GSIRateLaw
+{
+public:
+ GSIRateLawLHAblaArrhenius(ARGS args)
+ : GSIRateLaw(args),
+ pos_T_trans(0)
+ {
+ assert(args.s_node_rate_law.tag() == "LH_abla_arrhenius");
+
+ args.s_node_rate_law.getAttribute( "pre_exp", m_pre_exp,
+ "The sticking coeffcient probability for the reaction "
+ "should be provided.");
+ args.s_node_rate_law.getAttribute( "T", m_T_act,
+ "The activation temperature for the reaction "
+ "should be provided for an adsorption reaction.");
+ }
+
+//==============================================================================
+
+ ~GSIRateLawLHAblaArrhenius( ){ }
+
+//==============================================================================
+
+ double forwardReactionRateCoefficient(
+ const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_Tsurf) const
+ {
+ const double Tsurf = v_Tsurf(pos_T_trans);
+
+ return m_pre_exp*exp(-m_T_act / Tsurf);
+ }
+
+private:
+ const size_t pos_T_trans;
+
+ double m_pre_exp;
+ double m_T_act;
+};
+
+ObjectProvider<
+ GSIRateLawLHAblaArrhenius, GSIRateLaw>
+ gsi_rate_law_lh_abla_arrhenius("LH_abla_arrhenius");
+
+ } // namespace GasSurfaceInteraction
+} // namespace Mutation
diff --git a/src/gsi/GSIRateManagerDetailed.cpp b/src/gsi/GSIRateManagerDetailed.cpp
index c2e4128d..a4ef532d 100644
--- a/src/gsi/GSIRateManagerDetailed.cpp
+++ b/src/gsi/GSIRateManagerDetailed.cpp
@@ -91,8 +91,8 @@ class GSIRateManagerDetailed :
}
// Setup NewtonSolver
- setMaxIterations(5);
- setWriteConvergenceHistory(true);
+ setMaxIterations(50);
+ setWriteConvergenceHistory(false);
setEpsilon(m_tol);
}
@@ -170,10 +170,6 @@ class GSIRateManagerDetailed :
if (is_surf_cov_steady_state)
computeSurfaceSteadyStateCoverage();
- // double B = mv_sigma(0);
- // cout << scientific << setprecision(100);
- // std::cout << "mpp kf1 = " << mv_kf(0)*B << " kf2 = " << mv_kf(1)*B << std::endl;
-
mv_rate = mv_kf;
m_reactants.multReactions(mv_nd, mv_rate);
@@ -217,7 +213,7 @@ class GSIRateManagerDetailed :
for (int i = 0; i < mn_site_sp; ++i) {
double m_X_unpert = v_X(i);
double pert = v_X(i) * m_pert;
- v_X(i) += pert;
+ v_X(i) += pert;
// Update Jacobian column
updateFunction(v_X);
@@ -256,21 +252,20 @@ class GSIRateManagerDetailed :
double norm()
{
- return mv_f_unpert.lpNorm();
+ //return mv_f_unpert.tail(mn_site_sp).lpNorm();
+ return mv_dX.lpNorm();
}
//=============================================================================
private:
void computeSurfaceSteadyStateCoverage(){
- // mv_nd.head(m_ns) = 1.;
- // mv_X = mv_nd.tail(mn_site_sp);
- mv_X.setConstant(5.e19);
+ mv_X.setZero();
+ for (int i = 0; i < mn_site_cat; ++i)
+ mv_X.setConstant(mv_sigma(i) / mv_sp_in_site[i]);
- // std::cout << "Before X = \n" << mv_X << std::endl;
mv_X = solve(mv_X);
applyTolerance(mv_X);
- // std::cout << "After X = \n" << mv_X << std::endl;
// Setting up the SurfaceSiteCoverage.
// This is not essential for efficiency.
@@ -281,7 +276,8 @@ class GSIRateManagerDetailed :
//=============================================================================
inline void applyTolerance(Eigen::VectorXd& v_x) const {
for (int i = 0; i < v_x.size(); i++)
- if (std::abs(v_x(i) / mv_sigma(0)) < m_tol) v_x(i) = 0.;
+ if (std::abs(v_x(i) / mv_sigma(0)) < 1e-30) v_x(i) = 0.;
+ //if (std::abs(v_x(i) / mv_sigma(0)) < m_tol) v_x(i) = 0.; //prevent system to converge
}
//=============================================================================
@@ -323,5 +319,8 @@ class GSIRateManagerDetailed :
ObjectProvider
gsi_rate_manager_detailed_mass("detailed_mass");
+ObjectProvider
+ gsi_rate_manager_detailed_mass_energy("detailed_mass_energy");
+
} // namespace GasSurfaceInteraction
} // namespace Mutation
diff --git a/src/gsi/SurfaceBalanceSolverMass.cpp b/src/gsi/SurfaceBalanceSolverMass.cpp
index e28c7191..36e49042 100644
--- a/src/gsi/SurfaceBalanceSolverMass.cpp
+++ b/src/gsi/SurfaceBalanceSolverMass.cpp
@@ -67,7 +67,7 @@ class SurfaceBalanceSolverMass :
mv_f_unpert(m_ns),
m_jac(m_ns, m_ns),
m_tol(1.e-13),
- m_pert(1.e-2),
+ m_pert(1.e-4),
mv_X_unpert(m_ns),
pos_T_trans(0),
set_state_with_rhoi_T(1),
diff --git a/src/gsi/SurfaceBalanceSolverMassEnergy.cpp b/src/gsi/SurfaceBalanceSolverMassEnergy.cpp
index b855ff59..eb82d880 100644
--- a/src/gsi/SurfaceBalanceSolverMassEnergy.cpp
+++ b/src/gsi/SurfaceBalanceSolverMassEnergy.cpp
@@ -75,8 +75,8 @@ class SurfaceBalanceSolverMassEnergy :
mv_f_unpert(m_neqns),
m_jac(m_neqns, m_neqns),
m_tol(1.e-12),
- m_pert_m(1.e-2),
- m_pert_T(1.e0),
+ m_pert_m(1.e-4),
+ m_pert_T(1.e-2),
pos_E(m_ns),
pos_T_trans(0),
m_phi(m_surf_state.getSolidProperties().getPhiRatio()),
@@ -443,5 +443,9 @@ ObjectProvider<
SurfaceBalanceSolverMassEnergy, Surface>
surface_balance_solver_phenomenological_mass_energy("phenomenological_mass_energy");
+ObjectProvider<
+ SurfaceBalanceSolverMassEnergy, Surface>
+ surface_balance_solver_detailed_mass_energy("detailed_mass_energy");
+
} // namespace GasSurfaceInteraction
} // namespace Mutation
diff --git a/src/gsi/SurfacePropertiesDetailed.cpp b/src/gsi/SurfacePropertiesDetailed.cpp
index 13cd46f2..4c1facec 100644
--- a/src/gsi/SurfacePropertiesDetailed.cpp
+++ b/src/gsi/SurfacePropertiesDetailed.cpp
@@ -301,7 +301,15 @@ class SurfacePropertiesDetailed : public SurfaceProperties
const int non_empty = v_species.size();
for (int i_sp = 0; i_sp < non_empty; ++i_sp) {
- int id_sp = m_thermo.speciesIndex(v_species[i_sp]);
+ string species = v_species[i_sp];
+ int id_sp = -1;
+ //identify species labelled with *, e.g. O*-s
+ if (species.back() == '*') {
+ id_sp = m_thermo.speciesIndex(species.substr(0, species.size()-1));
+ }
+ else {
+ id_sp = m_thermo.speciesIndex(v_species[i_sp]);
+ }
if (id_sp == -1) {
throw InvalidInputError("SurfaceProperties",
diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp
index f0df0e0d..a09fee35 100644
--- a/src/kinetics/Kinetics.cpp
+++ b/src/kinetics/Kinetics.cpp
@@ -31,6 +31,7 @@
#include
+
using namespace std;
using namespace Eigen;
using namespace Mutation::Thermodynamics;
diff --git a/src/thermo/ParticleRRHO.cpp b/src/thermo/ParticleRRHO.cpp
index 91d4b9ca..e1386358 100644
--- a/src/thermo/ParticleRRHO.cpp
+++ b/src/thermo/ParticleRRHO.cpp
@@ -80,15 +80,19 @@ ParticleRRHO::ParticleRRHO(const IO::XmlElement& xml_element)
int degeneracy;
double temperature;
-
+ //int dummy = 1;
+
for ( ; level_iter != iter->end(); ++level_iter) {
if (level_iter->tag() == "level") {
level_iter->getAttribute("degeneracy", degeneracy);
level_iter->getAttribute("energy", temperature);
// convert from 1/cm to K
+ // if (dummy == 1) {
m_electronic_energies.push_back(
std::make_pair(degeneracy, temperature * 1.4387));
+ //dummy +=1;
+ //}
}
}
}
diff --git a/src/transfer/OmegaCV.cpp b/src/transfer/OmegaCV.cpp
index 519834cd..9b6701ee 100644
--- a/src/transfer/OmegaCV.cpp
+++ b/src/transfer/OmegaCV.cpp
@@ -79,9 +79,12 @@ class OmegaCV : public TransferModel
{
static int i_transfer_model = 0;
switch (i_transfer_model){
- case 0:
+ case 0:
return compute_source_Candler();
- break;
+ break;
+ case 1:
+ return compute_source_preferential();
+ break;
default:
std::cerr << "The selected Chemistry-Vibration-Chemistry model is not implemented yet";
return 0.0;
@@ -92,8 +95,12 @@ class OmegaCV : public TransferModel
int m_ns;
double* mp_wrk1;
double* mp_wrk2;
+ const int m_mol = 6;
+ const char* m_molecules[6] = { "N2", "O2", "NO", "N2+", "O2+", "NO+"};
+ const double m_energy[6] = { 9.759, 5.115, 6.496, 8.712, 6.663, 10.85};
double const compute_source_Candler();
+ double const compute_source_preferential();
};
/**
@@ -109,7 +116,6 @@ class OmegaCV : public TransferModel
double const OmegaCV::compute_source_Candler()
{
-
// Getting Vibrational Energy
m_mixture.speciesHOverRT(NULL, NULL, NULL, mp_wrk1, NULL, NULL);
@@ -126,6 +132,40 @@ double const OmegaCV::compute_source_Candler()
return(c1*sum*m_mixture.T()*RU);
}
+ /**
+ * Preferential Model according to Gnoffo with
+ *
+ * For more information:
+ *
+ * Gnoffo, Gupta, Shinn, Conservation equations and Physics modeles for Hypersonic Air Flows un Thermal and Chemical Nonequilibrium,
+ * Nasa technical report, 1989
+ *
+ */
+
+double const OmegaCV::compute_source_preferential()
+{
+
+ // Getting Production Rate
+ m_mixture.netProductionRates(mp_wrk2);
+
+ double DissEnergy[m_ns] = { };
+ int index;
+ for (int i=0; i < m_mol; ++i) {
+ index = m_mixture.speciesIndex(m_molecules[i]);
+ if ( index >= 0)
+ DissEnergy[index] = m_energy[i]*1000.0*96.487;
+ }
+
+ // Inner Product
+ double c1 = 0.3E0;
+ double sum = 0.E0;
+
+ for(int i = 0 ; i < m_ns; ++i)
+ sum += c1 * mp_wrk2[i] * DissEnergy[i]/m_mixture.speciesMw(i);
+
+ return sum;
+ }
+
// Register the transfer model
Mutation::Utilities::Config::ObjectProvider<
OmegaCV, TransferModel> omegaCV("OmegaCV");
diff --git a/src/transfer/OmegaVT.cpp b/src/transfer/OmegaVT.cpp
index 9cd848a3..fcaa2a2c 100644
--- a/src/transfer/OmegaVT.cpp
+++ b/src/transfer/OmegaVT.cpp
@@ -30,6 +30,7 @@
#include "TransferModel.h"
#include
+
using namespace Mutation;
namespace Mutation {
@@ -97,6 +98,7 @@ class OmegaVT : public TransferModel
src += p_Y[iv]*rho*RU*T/mp_Mw[iv]*(mp_hveq[iv] - mp_hv[iv])/compute_tau_VT_m(iv-inv);
}
}
+
return src;
}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index f78e999c..0c3ccf6b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -34,6 +34,8 @@ add_executable(run_tests
test_energies.cpp
test_errors.cpp
test_gsi_detailed.cpp
+ test_gsi_detailed_ACA.cpp
+ test_gsi_detailed_mass_energy.cpp
test_gsi_mass.cpp
test_gsi_mass_energy.cpp
test_gsi_rates.cpp
diff --git a/tests/data/gsi/seb_FRC_nitridarion.xml b/tests/data/gsi/seb_FRC_nitridarion.xml
new file mode 100644
index 00000000..c773b3ab
--- /dev/null
+++ b/tests/data/gsi/seb_FRC_nitridarion.xml
@@ -0,0 +1,43 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/data/gsi/smb_FRC2_nitridarion.xml b/tests/data/gsi/smb_FRC2_nitridarion.xml
new file mode 100644
index 00000000..e966bac2
--- /dev/null
+++ b/tests/data/gsi/smb_FRC2_nitridarion.xml
@@ -0,0 +1,40 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/data/gsi/smb_full_ACA.xml b/tests/data/gsi/smb_full_ACA.xml
new file mode 100644
index 00000000..b490d2e0
--- /dev/null
+++ b/tests/data/gsi/smb_full_ACA.xml
@@ -0,0 +1,111 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tests/data/mixtures/seb_FRC_nitridation_NASA9_ChemNonEq1T.xml b/tests/data/mixtures/seb_FRC_nitridation_NASA9_ChemNonEq1T.xml
new file mode 100644
index 00000000..1a9e3f97
--- /dev/null
+++ b/tests/data/mixtures/seb_FRC_nitridation_NASA9_ChemNonEq1T.xml
@@ -0,0 +1,13 @@
+
+
+
+ N N2 C CN
+
+
+
+ C:0., N:1.
+
+
diff --git a/tests/data/mixtures/smb_FRC2_nitridation_NASA9_ChemNonEq1T.xml b/tests/data/mixtures/smb_FRC2_nitridation_NASA9_ChemNonEq1T.xml
new file mode 100644
index 00000000..9d24ca8f
--- /dev/null
+++ b/tests/data/mixtures/smb_FRC2_nitridation_NASA9_ChemNonEq1T.xml
@@ -0,0 +1,13 @@
+
+
+
+ N N2 C CN
+
+
+
+ C:0., N:1.
+
+
diff --git a/tests/data/mixtures/smb_FRC_nitridation_NASA9_ChemNonEq1T.xml b/tests/data/mixtures/smb_FRC_nitridation_NASA9_ChemNonEq1T.xml
new file mode 100644
index 00000000..eb048384
--- /dev/null
+++ b/tests/data/mixtures/smb_FRC_nitridation_NASA9_ChemNonEq1T.xml
@@ -0,0 +1,13 @@
+
+
+
+ N N2 C CN
+
+
+
+ C:0., N:1.
+
+
diff --git a/tests/data/mixtures/smb_full_ACA_NASA9_ChemNonEq1T.xml b/tests/data/mixtures/smb_full_ACA_NASA9_ChemNonEq1T.xml
new file mode 100644
index 00000000..c26eff87
--- /dev/null
+++ b/tests/data/mixtures/smb_full_ACA_NASA9_ChemNonEq1T.xml
@@ -0,0 +1,13 @@
+
+
+
+ O O2 N N2 C CO CO2 CN
+
+
+
+ C:0., N:0.79, O:0.21
+
+
diff --git a/tests/test_gsi_detailed.cpp b/tests/test_gsi_detailed.cpp
index b3b6c4ba..553d0701 100644
--- a/tests/test_gsi_detailed.cpp
+++ b/tests/test_gsi_detailed.cpp
@@ -245,8 +245,10 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]")
// MixtureOptions optspark("smb_oxidation_NASA9_ChemNonEq1T");
// Mixture mixpark(optspark);
+ const double tol = 10e-4;
+
size_t ns = 5;
- size_t nr = 2;
+ size_t nr = 5;
CHECK(mixpsmm.nSpecies() == ns);
CHECK(mixpsmm.nSurfaceReactions() == nr);
@@ -300,8 +302,16 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]")
v_surf_cov_mpp_frac = mixpsmm.getSurfaceProperties().getSurfaceSiteCoverageFrac();
v_surf_cov_mpp_frac *= B;
- std::cout << "Coverage MPP = " << v_surf_cov_mpp_frac(0) << " " << v_surf_cov_mpp_frac(1) << std::endl;
- std::cout << "Coverage HERE = " << v_surf_cov_frac(0) << " " << v_surf_cov_frac(1) << std::endl;
+ //std::cout << "Coverage MPP = " << v_surf_cov_mpp_frac(0) << " " << v_surf_cov_mpp_frac(1) << std::endl;
+ //std::cout << "Coverage HERE = " << v_surf_cov_frac(0) << " " << v_surf_cov_frac(1) << std::endl;
+
+ CHECK(v_surf_cov_mpp_frac(0) >= 0.0);
+ CHECK(v_surf_cov_mpp_frac(1) >= 0.0);
+ CHECK((B - v_surf_cov_mpp_frac.sum())/B == Approx(0.0).epsilon(tol));
+ if (v_surf_cov_mpp_frac(0)/B >= 1.0e-14)
+ CHECK((v_surf_cov_frac(0) - v_surf_cov_mpp_frac(0))/v_surf_cov_mpp_frac(0) == Approx(0.0).epsilon(tol));
+ else
+ CHECK(v_surf_cov_frac(0)/B <= 1.0e-14);
rates(0) = kfads * nO * v_surf_cov_frac(0);
rates(1) = kfdes * v_surf_cov_frac(1);
@@ -309,35 +319,49 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]")
rates(3) = kfer2 * nO *v_surf_cov_frac(1);
rates(4) = kfer3 * nO * v_surf_cov_frac(0);
- std::cout << "T = " << T << " P = " << P << std::endl;
- std::cout << "MPP Rates = " << ratesmpp(0) << " "
- << ratesmpp(1) << " "
- << ratesmpp(2) << " "
- << ratesmpp(3) << " "
- << ratesmpp(4) << std::endl;
- std::cout << "HERE Rates = " << rates(0) << " "
- << rates(1) << " "
- << rates(2) << " "
- << rates(3) << " "
- << rates(4) << std::endl;
-
- wdot(0) = mm(iO) / NA * (-rates(0) + rates(1) - rates(2) - rates(3) - rates(4));
+ //std::cout << "T = " << T << " P = " << P << std::endl;
+ //std::cout << "MPP Rates = " << ratesmpp(0) << " "
+ // << ratesmpp(1) << " "
+ // << ratesmpp(2) << " "
+ // << ratesmpp(3) << " "
+ // << ratesmpp(4) << std::endl;
+ //std::cout << "HERE Rates = " << rates(0) << " "
+ // << rates(1) << " "
+ // << rates(2) << " "
+ // << rates(3) << " "
+ // << rates(4) << std::endl;
+
+ for (int ii = 0; ii < 5; ++ii ) {
+ if (abs(ratesmpp(ii)) >= 1.0e-14)
+ CHECK((rates(ii) - ratesmpp(ii))/ratesmpp(ii) == Approx(0.0).epsilon(tol));
+ else
+ CHECK(abs(rates(ii)) <= 1.0e-14);
+ }
+
+ wdot(0) = - mm(iO) / NA * (-rates(0) + rates(1) - rates(2) - rates(3) - rates(4));
wdot(1) = 0.;
wdot(2) = 0.;
- wdot(3) = mm(iCO) / NA * (+ rates(2) + rates(4));
- wdot(4) = mm(iCO2) / NA * (+ rates(3));
-
- std::cout << "T = " << T << " P = " << P << std::endl;
- std::cout << "MPP Rates = " << wdotmpp(0) << " "
- << wdotmpp(1) << " "
- << wdotmpp(2) << " "
- << wdotmpp(3) << " "
- << wdotmpp(4) << std::endl;
- std::cout << "HERE Rates = " << wdot(0) << " "
- << wdot(1) << " "
- << wdot(2) << " "
- << wdot(3) << " "
- << wdot(4) << std::endl;
+ wdot(3) = - mm(iCO) / NA * (+ rates(2) + rates(4));
+ wdot(4) = - mm(iCO2) / NA * (+ rates(3));
+
+ //std::cout << "T = " << T << " P = " << P << std::endl;
+ //std::cout << "MPP Rates = " << wdotmpp(0) << " "
+ // << wdotmpp(1) << " "
+ // << wdotmpp(2) << " "
+ // << wdotmpp(3) << " "
+ // << wdotmpp(4) << std::endl;
+ //std::cout << "HERE Rates = " << wdot(0) << " "
+ // << wdot(1) << " "
+ // << wdot(2) << " "
+ // << wdot(3) << " "
+ // << wdot(4) << std::endl;
+
+ for (int ii = 0; ii < 4; ++ii ) {
+ if (abs(wdotmpp(ii)) >= 1e-14)
+ CHECK((wdot(ii) - wdotmpp(ii))/wdot(ii) == Approx(0.0).epsilon(tol));
+ else
+ CHECK(abs(wdot(ii)) <= 1.0e-14);
+ }
// Park
// mixpark.equilibrate(T, P);
@@ -346,8 +370,241 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]")
// mixpark.setSurfaceState(v_rhoi.data(), &T, set_state_rhoi_T);
// nO = mixpark.X()[iO] * mixpark.numberDensity();
- std::cout << "End" << std::endl;
- double in; std::cin >> in;
+ //std::cout << "End" << std::endl;
+ //double in; std::cin >> in;
+
+ }
+
+ P *= dP;
+ }
+
+ }
+
+/* SECTION("Nitridation Model.")
+ {
+ // Setting up M++
+ MixtureOptions optsFRC("smb_FRC_nitridation_NASA9_ChemNonEq1T");
+ Mixture mixFRC(optsFRC);
+
+ size_t ns = 4;
+ size_t nr = 4;
+
+ CHECK(mixFRC.nSpecies() == ns);
+ CHECK(mixFRC.nSurfaceReactions() == nr);
+
+ const size_t iN = 0;
+ const size_t iN2 = 1;
+ const size_t iCN = 3;
+
+ const int set_state_rhoi_T = 1;
+
+ const double tol = 10e-4;
+
+ ArrayXd v_rhoi(ns);
+ ArrayXd wdot(ns); ArrayXd wdotmpp(ns);
+ ArrayXd rates(nr); ArrayXd ratesmpp(nr);
+ wdot.setZero(); wdotmpp.setZero();
+
+ ArrayXd mm = mixFRC.speciesMw();
+
+ CHECK(mixFRC.getSurfaceProperties().isSurfaceCoverageSteady() == true);
+ ArrayXd v_surf_cov_mpp_frac(mixFRC.getSurfaceProperties().nSurfaceSpecies());
+ ArrayXd v_surf_cov_frac(mixFRC.getSurfaceProperties().nSurfaceSpecies());
+
+ // Equilibrium Surface
+ double P = 1.e-5;
+ double dP = 10.;
+ double T; // K
+ double dT = 200.; // K
+ for (int i = 0; i < 15; i++) {
+ for (int j = 0; j < 16; j++) {
+ T = (j+1) * dT;
+
+ mixFRC.equilibrate(T, P);
+ mixFRC.densities(v_rhoi.data());
+
+ mixFRC.setSurfaceState(v_rhoi.data(), &T, set_state_rhoi_T);
+ double nN = mixFRC.X()[iN] * mixFRC.numberDensity();
+
+ mixFRC.surfaceReactionRatesPerReaction(ratesmpp.data());
+ mixFRC.surfaceReactionRates(wdotmpp.data());
+
+ const double B = 6.022e18;
+ double F = 1./B * sqrt(RU * T / (2 * PI * mm(iN)));
+ double kfads = F*exp(-7500./T);
+ double kfdes = 2 * PI * mm(iN) / NA * KB * KB * T / (HP * HP * HP) / B;
+ kfdes *= exp(-73971.6/T);
+
+ double kfer1 = F * 9.0e+5 * exp(-20676./T);
+ double kfer2 = F * 1.1e+6 * exp(-18000./T);
+
+ v_surf_cov_frac(0) = (kfdes + nN*(kfer1+kfer2))/(nN*(kfads + kfer1 + kfer2) + kfdes)*B;
+ v_surf_cov_frac(1) = B - v_surf_cov_frac(0);
+
+ v_surf_cov_mpp_frac = mixFRC.getSurfaceProperties().getSurfaceSiteCoverageFrac();
+ v_surf_cov_mpp_frac *= B;
+
+ //Check total number of sites is respected and free spot is the same of analytical solution
+ CHECK(v_surf_cov_mpp_frac(0) >= 0.0);
+ CHECK(v_surf_cov_mpp_frac(1) >= 0.0);
+ CHECK((B - v_surf_cov_mpp_frac.sum())/B == Approx(0.0).epsilon(tol));
+ if (v_surf_cov_mpp_frac(0)/B >= 1.0e-14)
+ CHECK((v_surf_cov_frac(0) - v_surf_cov_mpp_frac(0))/v_surf_cov_mpp_frac(0) == Approx(0.0).epsilon(tol));
+ else
+ CHECK(v_surf_cov_frac(0)/B <= 1.0e-14);
+
+ //Check rates
+ rates(0) = kfads * nN * v_surf_cov_frac(0);
+ rates(1) = kfdes * v_surf_cov_frac(1);
+ rates(2) = kfer1 * nN * v_surf_cov_frac(1);
+ rates(3) = kfer2 * nN * v_surf_cov_frac(1);
+
+ for (int ii = 0; ii < 4; ++ii ) {
+ if (abs(ratesmpp(ii)) >= 1.0e-14)
+ CHECK((rates(ii) - ratesmpp(ii))/ratesmpp(ii) == Approx(0.0).epsilon(tol));
+ else
+ CHECK(abs(rates(ii)) <= 1.0e-14);
+ }
+
+ //Check chemical production
+ wdot(0) = - mm(iN) / NA * (-rates(0) + rates(1) - rates(3));
+ wdot(1) = - mm(iN2) / NA * (rates(3));
+ wdot(2) = 0.;
+ wdot(3) = -mm(iCN) / NA * rates(2);
+
+ for (int ii = 0; ii < 4; ++ii ) {
+ if (abs(wdotmpp(ii)) >= 1e-14)
+ CHECK((wdot(ii) - wdotmpp(ii))/wdot(ii) == Approx(0.0).epsilon(tol));
+ else
+ CHECK(abs(wdot(ii)) <= 1.0e-14);
+ }
+
+ }
+
+ P *= dP;
+ }
+
+ }*/
+
+
+ SECTION("Nitridation Model 2.")
+ {
+ // Setting up M++
+ MixtureOptions optsFRC("smb_FRC2_nitridation_NASA9_ChemNonEq1T");
+ Mixture mixFRC(optsFRC);
+
+ size_t ns = 4;
+ size_t nr = 6;
+
+ CHECK(mixFRC.nSpecies() == ns);
+ CHECK(mixFRC.nSurfaceReactions() == nr);
+
+ const size_t iN = 0;
+ const size_t iN2 = 1;
+ const size_t iCN = 3;
+
+ const int set_state_rhoi_T = 1;
+
+ const double tol = 10e-4;
+
+ ArrayXd v_rhoi(ns);
+ ArrayXd wdot(ns); ArrayXd wdotmpp(ns);
+ ArrayXd rates(nr); ArrayXd ratesmpp(nr);
+ wdot.setZero(); wdotmpp.setZero();
+
+ ArrayXd mm = mixFRC.speciesMw();
+ mm /= NA;
+
+ CHECK(mixFRC.getSurfaceProperties().isSurfaceCoverageSteady() == true);
+ ArrayXd v_surf_cov_mpp_frac(mixFRC.getSurfaceProperties().nSurfaceSpecies());
+ ArrayXd v_surf_cov_frac(mixFRC.getSurfaceProperties().nSurfaceSpecies());
+
+ // Equilibrium Surface, problem at T < 1000.0 @TODO
+ double P = 1.e-5;
+ double dP = 10.;
+ double T; // K
+ double dT = 200.; // K
+ for (int i = 0; i < 15; i++) {
+ for (int j = 0; j < 12; j++) {
+ T = (j+6) * dT;
+ //for (int i = 0; i < 1; i++) {
+ // for (int j = 0; j < 1; j++) {
+ // T = (j+4) * dT;
+
+ mixFRC.equilibrate(T, P);
+ mixFRC.densities(v_rhoi.data());
+
+ mixFRC.setSurfaceState(v_rhoi.data(), &T, set_state_rhoi_T);
+ double nN = mixFRC.X()[iN] * mixFRC.numberDensity();
+
+ mixFRC.surfaceReactionRatesPerReaction(ratesmpp.data());
+ mixFRC.surfaceReactionRates(wdotmpp.data());
+
+ const double B = 6.022e18; //[#/m2]
+ double F = 1./B * sqrt(KB * T / (2 * PI * mm(iN))); //[m3/s/#]
+ double F2D = sqrt(PI * KB * T / (2 * mm(iN))); //[m/s]
+ double kfads = F*exp(-2500./T); //[m3/s/#]
+ double kfdes = 2.0 * PI * mm(iN) * KB * KB * T * T / (HP * HP * HP) / B;
+ kfdes *= exp(-73971.6/T); //[1/s]
+ double kfer1 = F * 1.5 * exp(-7000./T); //[m3/s/#]
+ double kfer2 = F * 0.5 * exp(-2000./T); //[m3/s/#]
+ double kflh1 = sqrt(1.0/(B)) * F2D * 0.1 * exp(-21000.0/ T); //[m2/s]
+ double kflh2 = 1.0e8 * exp(-20676.0 / T); //[1/s]
+
+ double A = 2.0*kflh1;
+ double BB = (kfads + kfer1 + kfer2)*nN + kfdes + kflh2;
+ double C = kfads * B * nN;
+
+ v_surf_cov_frac(1) = (sqrt(BB*BB + 4.*A*C) - BB )/ (2.0 * A);
+ v_surf_cov_frac(0) = B - v_surf_cov_frac(1);
+
+ //v_surf_cov_frac(0) = (kfdes + nN*(kfer1+kfer2))/(nN*(kfads + kfer1 + kfer2) + kfdes)*B;
+ //v_surf_cov_frac(1) = B - v_surf_cov_frac(0);
+
+ v_surf_cov_mpp_frac = mixFRC.getSurfaceProperties().getSurfaceSiteCoverageFrac();
+ v_surf_cov_mpp_frac *= B;
+
+ //std::cout << "analytic:" << v_surf_cov_frac(0)/B << " " << v_surf_cov_frac(1)/B<< std::endl;
+ //std::cout << "MPP: " << v_surf_cov_mpp_frac(0)/B << " " << v_surf_cov_mpp_frac(1)/B<< std::endl;
+
+ //Check total number of sites is respected and free spot is the same of analytical solution
+ CHECK(v_surf_cov_mpp_frac(0) >= 0.0);
+ CHECK(v_surf_cov_mpp_frac(1) >= 0.0);
+ CHECK((B - v_surf_cov_mpp_frac.sum())/B == Approx(0.0).epsilon(tol));
+ CHECK((v_surf_cov_frac(0) - v_surf_cov_mpp_frac(0))/v_surf_cov_mpp_frac(0) == Approx(0.0).epsilon(tol));
+
+ //Check rates
+ rates(0) = kfads * nN * v_surf_cov_frac(0);
+ rates(1) = kfdes * v_surf_cov_frac(1);
+ rates(2) = kfer1 * nN * v_surf_cov_frac(1);
+ rates(3) = kfer2 * nN * v_surf_cov_frac(1);
+ rates(4) = kflh1 * v_surf_cov_frac(1)* v_surf_cov_frac(1);
+ rates(5) = kflh2 * v_surf_cov_frac(1);
+
+ //std::cout << "Rates are" << std::endl;
+ //std::cout << rates(0) << " " << ratesmpp(0) << std::endl;
+ //std::cout << rates(1) << " " << ratesmpp(1) << std::endl;
+ //std::cout << rates(2) << " " << ratesmpp(2) << std::endl;
+ //std::cout << rates(3) << " " << ratesmpp(3) << std::endl;
+ //std::cout << rates(4) << " " << ratesmpp(4) << std::endl;
+ //std::cout << rates(5) << " " << ratesmpp(5) << std::endl;
+
+ for (int ii = 0; ii < 6; ++ii ) {
+ CHECK((rates(ii) - ratesmpp(ii))/ratesmpp(ii) == Approx(0.0).epsilon(tol));
+ }
+
+ //Check chemical production
+ wdot(0) = - mm(iN) * (-rates(0) + rates(1) - rates(3));
+ wdot(1) = - mm(iN2) * (rates(3) + rates(4));
+ wdot(2) = 0.;
+ wdot(3) = -mm(iCN) * (rates(2) + rates(5));
+
+ for (int ii = 0; ii < 4; ++ii ) {
+ if (abs(wdotmpp(ii)) >= 1e-14)
+ CHECK((wdot(ii) - wdotmpp(ii))/wdot(ii) == Approx(0.0).epsilon(tol));
+ else
+ CHECK(abs(wdot(ii)) <= 1.0e-14);
+ }
}
@@ -356,4 +613,4 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]")
}
-}
\ No newline at end of file
+}
diff --git a/tests/test_gsi_detailed_ACA.cpp b/tests/test_gsi_detailed_ACA.cpp
new file mode 100644
index 00000000..6e00b93e
--- /dev/null
+++ b/tests/test_gsi_detailed_ACA.cpp
@@ -0,0 +1,156 @@
+/*
+ * Copyright 2014-2018 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
+ * .
+ */
+
+#include "mutation++.h"
+#include "Configuration.h"
+#include "TestMacros.h"
+#include
+#include
+
+#include "SurfaceProperties.h"
+
+using namespace Mutation;
+using namespace Catch;
+using namespace Eigen;
+
+void computeRates(double Rs) {
+ Rs = 3.0;
+}
+
+TEST_CASE("Detailed surface chemictry tests full ACA.","[gsi]")
+{
+ const double tol = 100. * std::numeric_limits::epsilon();
+ const double tol_det = 1.e2 * std::numeric_limits::epsilon();
+
+ Mutation::GlobalOptions::workingDirectory(TEST_DATA_FOLDER);
+
+ SECTION("Surface Species and Coverage.")
+ {
+ // Setting up M++
+ MixtureOptions opts("smb_full_ACA_NASA9_ChemNonEq1T");
+ Mixture mix(opts);
+
+ CHECK(mix.nSpecies() == 8);
+
+ // Check global options
+ CHECK(mix.nSurfaceReactions() == 20);
+ CHECK(mix.getSurfaceProperties().nSurfaceSpecies() == 5);
+ CHECK(mix.getSurfaceProperties().nSiteSpecies() == 4);
+
+ // Check Species
+ CHECK(mix.getSurfaceProperties().surfaceSpeciesIndex("s") == 8);
+ CHECK(mix.getSurfaceProperties().surfaceSpeciesIndex("N-s") == 9);
+ CHECK(mix.getSurfaceProperties().surfaceSpeciesIndex("O-s") == 10);
+ CHECK(mix.getSurfaceProperties().surfaceSpeciesIndex("O*-s") == 11);
+ CHECK(mix.getSurfaceProperties().surfaceSpeciesIndex("C-b") == 12);
+ CHECK(mix.getSurfaceProperties().surfaceSpeciesIndex("A") == -1);
+
+ // Check surface species association with gaseous species
+ CHECK( mix.getSurfaceProperties().surfaceToGasIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("N-s")) == 2);
+ CHECK( mix.getSurfaceProperties().surfaceToGasIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("O-s")) == 0);
+ CHECK( mix.getSurfaceProperties().surfaceToGasIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("O*-s")) == 0);
+ CHECK(mix.getSurfaceProperties().surfaceToGasIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("C-b")) == 4);
+
+ CHECK(mix.getSurfaceProperties().surfaceToGasIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("s")) == -2);
+ CHECK(mix.getSurfaceProperties().surfaceToGasIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("b")) == -1);
+ CHECK(mix.getSurfaceProperties().surfaceToGasIndex(100) == -1);
+
+ // Check site species map correctly to the site category
+ CHECK(mix.getSurfaceProperties().siteSpeciesToSiteCategoryIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("N-s")) == 0);
+ CHECK(mix.getSurfaceProperties().siteSpeciesToSiteCategoryIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("O-s")) == 0);
+ CHECK(mix.getSurfaceProperties().siteSpeciesToSiteCategoryIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("O*-s")) == 0);
+ CHECK(mix.getSurfaceProperties().siteSpeciesToSiteCategoryIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("C-b")) == -1);
+ CHECK(mix.getSurfaceProperties().siteSpeciesToSiteCategoryIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("s")) == 0);
+ CHECK(mix.getSurfaceProperties().siteSpeciesToSiteCategoryIndex(
+ mix.getSurfaceProperties().surfaceSpeciesIndex("b")) == -1);
+
+ CHECK(mix.getSurfaceProperties().nSiteDensityInCategory(0) == 6.022e18);
+ CHECK(mix.getSurfaceProperties().nSiteDensityInCategory(3) == -1);
+
+ // Check reaction rates
+ double P = 1500; //Pa
+ double T = 2000; // K
+
+ int iO = 0;
+ int iN = 2;
+ int ns = mix.nSpecies();
+
+ ArrayXd v_rhoi(ns);
+ ArrayXd v_nd(ns);
+ ArrayXd mm = mix.speciesMw();
+ const int set_state_rhoi_T = 1;
+
+ mix.equilibrate(T, P);
+ mix.densities(v_rhoi.data());
+ for (int i = 0; i < ns; ++i)
+ v_nd(i) = v_rhoi(i) / mm(i);
+
+ double B = 6.022e18;
+ double FO = 1./B * sqrt(RU * T / (2 * PI * mm(iO)));
+ double FO2 = 1./B * sqrt(RU * T / (2 * PI * 2 * mm(iO)));
+ double FN = 1./B * sqrt(RU * T / (2 * PI * mm(iN)));
+
+ //coefficients
+ double k1 = FO * 0.3;
+ double k2 = 2.0 * PI * mm(iO) / NA * KB * KB * T * T / (HP * HP * HP) / B * exp(-44277.6/T);
+ double k3 = FO * 100. * exp(-4000/T);
+ double k4 = FO * exp(-500/T);
+ double k5 = FO * 0.7;
+ double k6 = 2.0 * PI * mm(iO) / NA * KB * KB * T * T / (HP * HP * HP) / B * exp(-96500.6/T);
+ double k7 = FO * 1000. * exp(-4000/T);
+ double k8 = sqrt(1.0/B) * sqrt(NA * PI * KB * T / (2.0 * mm(iO)))* 1.0e-3 * exp(-15000.0/ T);
+ double k9 = sqrt(1.0/B) * sqrt(NA * PI * KB * T / (2.0 * mm(iO)))* 5.0e-5 * exp(-15000.0/ T);
+ double k10 = FN * exp(-2500./T);
+ double k11 = 2.0 * PI * mm(iN) / NA * KB * KB * T * T / (HP * HP * HP) / B * exp(-73971.6/T);
+ double k12 = FN * 1.5 * exp(-7000./T);
+ double k13 = FN * 0.5 * exp(-2000./T);
+ double k14 = sqrt(1.0/(B)) * sqrt(NA*PI*KB*T/ (2.0 * mm(iN)))* 0.1 * exp(-21000.0/ T);
+ double k15 = 1e8 * exp(-20676.0 / T);
+ double k16 = FO2 / B * exp(-8000./T);
+ double k17 = FO2 * 100. * exp(-4000./T);
+ double k18 = FO2 * exp(-500./T);
+ double k19 = FO2 / B * exp(-8000./T);
+ double k20 = FO2 * 1000. * exp(-4000./T);
+
+ double r;
+ computeRates(r);
+
+ //compute steady state
+ mix.setSurfaceState(v_rhoi.data(), &T, set_state_rhoi_T);
+ mix.getSurfaceProperties();
+
+
+
+ }
+
+
+}
diff --git a/tests/test_gsi_detailed_mass_energy.cpp b/tests/test_gsi_detailed_mass_energy.cpp
new file mode 100644
index 00000000..1c13c6c9
--- /dev/null
+++ b/tests/test_gsi_detailed_mass_energy.cpp
@@ -0,0 +1,135 @@
+/*
+ * Copyright 2014-2018 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
+ * .
+ */
+
+#include "mutation++.h"
+#include "Configuration.h"
+#include "TestMacros.h"
+#include
+#include
+
+using namespace Mutation;
+using namespace Catch;
+using namespace Eigen;
+
+TEST_CASE("Solution of the Detailed MassEnergyBalanceSolver is converged.", "[gsi]")
+{
+ const double tol = std::numeric_limits::epsilon();
+ Mutation::GlobalOptions::workingDirectory(TEST_DATA_FOLDER);
+
+ SECTION("Detailed Mass and Energy Balance Ablation.")
+ {
+ // Mixture
+ MixtureOptions opts("seb_FRC_nitridation_NASA9_ChemNonEq1T");
+ Mixture mix(opts);
+
+ // Setting up
+ const size_t set_state_with_rhoi_T = 1;
+ const size_t pos_T_trans = 0;
+ size_t ns = mix.nSpecies();
+ size_t nT = mix.nEnergyEqns();
+ size_t neq = ns + nT;
+
+ // Conditions T = 3000K and p = 100Pa
+ VectorXd Teq = VectorXd::Constant(nT, 3000.);
+ double Peq = 100.; // Pa
+ mix.equilibrate(Teq(pos_T_trans), Peq);
+
+ // Setting number of iterations for the solver
+ const int iter = 100;
+ mix.setIterationsSurfaceBalance(iter);
+
+ // Mass gradient
+ VectorXd xi_e(ns);
+ xi_e = Map(mix.X(), ns);
+ double dx = 1.e-3;
+ mix.setDiffusionModel(xi_e.data(), dx);
+
+ // Temperature gradient
+ VectorXd T_e = Teq;
+ mix.setGasFourierHeatFluxModel(T_e.data(), dx);
+
+ // Initial conditions of the surface are the ones in the first
+ // physical cell
+ VectorXd rhoi_s(ns);
+ mix.densities(rhoi_s.data());
+ VectorXd T_s = VectorXd::Constant(nT, 1800.);
+ mix.setSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T);
+
+ // Solve balance and request solution
+ mix.solveSurfaceBalance();
+ mix.getSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T);
+ double rho = rhoi_s.sum();
+
+ // Verifying the solution gives low residual in the balance equations
+ mix.setState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T);
+ VectorXd xi_s(ns);
+ xi_s = Map(mix.X(), ns);
+
+ // Compute diffusion velocities
+ VectorXd dxidx(ns);
+ dxidx = (xi_s - xi_e) / dx;
+ VectorXd vdi(ns);
+ double E = 0.;
+ mix.stefanMaxwell(dxidx.data(), vdi.data(), E);
+
+ // Conductive heat flux
+ VectorXd dTdx(nT);
+ dTdx = (T_s - T_e) / dx;
+ VectorXd lambda(nT);
+ mix.frozenThermalConductivityVector(lambda.data());
+
+ // Get surface production rates
+ VectorXd wdot(ns);
+ mix.setSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T);
+ mix.surfaceReactionRates(wdot.data());
+
+ // Blowing flux (should be zero for catalysis)
+ double mblow;
+ mix.getMassBlowingRate(mblow);
+
+ // Species and mixture enthalpies
+ VectorXd v_hi(ns*nT);
+ mix.getEnthalpiesMass(v_hi.data());
+ VectorXd v_h(nT);
+ v_h(pos_T_trans) = (rhoi_s/rho).dot(v_hi.head(ns));
+
+ // Surface radiation
+ const double sigma = 2.*pow(PI, 5)*pow(KB, 4)/(15*pow(C0, 2)*pow(HP, 3));
+ const double eps = .86;
+ VectorXd q_srad = VectorXd::Zero(nT);
+ q_srad(pos_T_trans) = sigma*eps*pow(T_s(pos_T_trans), 4);
+
+ // Chemical Energy Contribution
+ VectorXd v_hi_rhoi_vi= VectorXd::Zero(nT);
+ v_hi_rhoi_vi(pos_T_trans) = -v_hi.head(ns).dot(rhoi_s.cwiseProduct(vdi));
+
+ // Building balance functions
+ VectorXd F(neq);
+ F.head(ns) = (rhoi_s/rho)*mblow + rhoi_s.cwiseProduct(vdi) - wdot;
+ F.tail(nT) = -lambda.cwiseProduct(dTdx) - q_srad + mblow*v_h -
+ v_hi_rhoi_vi;
+
+ // Compute error
+ double err = F.lpNorm();
+ CHECK(err == Approx(0.0).margin(tol));
+ }
+
+}
diff --git a/thirdparty/catch/include/internal/catch_fatal_condition.hpp b/thirdparty/catch/include/internal/catch_fatal_condition.hpp
index 1d6674fc..b37cc85c 100644
--- a/thirdparty/catch/include/internal/catch_fatal_condition.hpp
+++ b/thirdparty/catch/include/internal/catch_fatal_condition.hpp
@@ -136,7 +136,7 @@ namespace Catch {
static bool isSet;
static struct sigaction oldSigActions [sizeof(signalDefs)/sizeof(SignalDefs)];
static stack_t oldSigStack;
- static char altStackMem[SIGSTKSZ];
+ static char altStackMem[32768];
static void handleSignal( int sig ) {
std::string name = "";
@@ -156,7 +156,7 @@ namespace Catch {
isSet = true;
stack_t sigStack;
sigStack.ss_sp = altStackMem;
- sigStack.ss_size = SIGSTKSZ;
+ sigStack.ss_size = 32768;
sigStack.ss_flags = 0;
sigaltstack(&sigStack, &oldSigStack);
struct sigaction sa = { 0 };
@@ -188,7 +188,7 @@ namespace Catch {
bool FatalConditionHandler::isSet = false;
struct sigaction FatalConditionHandler::oldSigActions[sizeof(signalDefs)/sizeof(SignalDefs)] = {};
stack_t FatalConditionHandler::oldSigStack = {};
- char FatalConditionHandler::altStackMem[SIGSTKSZ] = {};
+ char FatalConditionHandler::altStackMem[32768] = {};
} // namespace Catch