From 69f0ba41e6661687cceeacb4fc288508156e00d9 Mon Sep 17 00:00:00 2001 From: michele capriati Date: Wed, 29 Jan 2020 14:05:50 +0100 Subject: [PATCH 01/14] Added nitrogen desorption k --- src/gsi/CMakeLists.txt | 1 + src/gsi/GSIRateLawDesorptionTST2.cpp | 105 +++++++++++++++++++++++++++ src/gsi/GSIRateManagerDetailed.cpp | 13 +++- 3 files changed, 117 insertions(+), 2 deletions(-) create mode 100644 src/gsi/GSIRateLawDesorptionTST2.cpp diff --git a/src/gsi/CMakeLists.txt b/src/gsi/CMakeLists.txt index b483766f..4328a9f4 100644 --- a/src/gsi/CMakeLists.txt +++ b/src/gsi/CMakeLists.txt @@ -30,6 +30,7 @@ add_sources(mutation++ GSIReactionDetailed.cpp GSIRateLawAdsorption.cpp GSIRateLawDesorptionTST.cpp + GSIRateLawDesorptionTST2.cpp GSIRateLawGammaConst.cpp GSIRateLawGammaT.cpp GSIRateLawERArrhenius.cpp diff --git a/src/gsi/GSIRateLawDesorptionTST2.cpp b/src/gsi/GSIRateLawDesorptionTST2.cpp new file mode 100644 index 00000000..f2b514bf --- /dev/null +++ b/src/gsi/GSIRateLawDesorptionTST2.cpp @@ -0,0 +1,105 @@ +/** + * @file GSIRateLawGammaT.cpp + * + * @brief Class which computes the reaction rate constant for a desorption + * surface reaction based on simple transition state theory. + */ + +/* + * 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 GSIRateLawDesorptionTST2 : public GSIRateLaw +{ +public: + GSIRateLawDesorptionTST2(ARGS args) + : GSIRateLaw(args), + m_surf_props(args.s_surf_props), + mv_react(args.s_reactants), + pos_T_trans(0), + idx_react(0) + { + assert(args.s_node_rate_law.tag() == "desorption_tst2"); + + args.s_node_rate_law.getAttribute("T", m_T_des, + "Activation temperature should be provided for every desorption " + "reaction."); + + m_mass_des = m_thermo.speciesMw( + args.s_surf_props.surfaceToGasIndex(mv_react[idx_react])) / NA; + m_site_categ = m_surf_props.siteSpeciesToSiteCategoryIndex( + mv_react[idx_react]); + m_n_sites = m_surf_props.nSiteDensityInCategory(m_site_categ); + // m_n_sites = args.s_surf_props.nSiteDensity(); + } + +//============================================================================== + + ~GSIRateLawDesorptionTST2( ){ } + +//============================================================================== + + double forwardReactionRateCoefficient( + const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_Tsurf) const + { + const double Tsurf = v_Tsurf(pos_T_trans); + + const double pre_exp = 2 * PI * m_mass_des * KB * KB * Tsurf + / (HP * HP * HP * m_n_sites); + + return pre_exp * exp(-m_T_des / Tsurf); + } + +private: + const size_t idx_react; + int m_site_categ; + + double m_T_des; + + double m_n_sites; + double m_mass_des; + + const double pos_T_trans; + + const std::vector& mv_react; + const SurfaceProperties& m_surf_props; +}; + +ObjectProvider< + GSIRateLawDesorptionTST2, GSIRateLaw> + gsi_rate_law_desorption_tst2("desorption_tst2"); + + } // namespace GasSurfaceInteraction +} // namespace Mutation diff --git a/src/gsi/GSIRateManagerDetailed.cpp b/src/gsi/GSIRateManagerDetailed.cpp index c2e4128d..fbe5bb67 100644 --- a/src/gsi/GSIRateManagerDetailed.cpp +++ b/src/gsi/GSIRateManagerDetailed.cpp @@ -92,7 +92,7 @@ class GSIRateManagerDetailed : // Setup NewtonSolver setMaxIterations(5); - setWriteConvergenceHistory(true); + setWriteConvergenceHistory(false); setEpsilon(m_tol); } @@ -264,7 +264,16 @@ class GSIRateManagerDetailed : // mv_nd.head(m_ns) = 1.; // mv_X = mv_nd.tail(mn_site_sp); - mv_X.setConstant(5.e19); + //mv_X.setConstant(5.e17); + //mv_X.setConstant(3.011e18); + + //mv_X.setZero(); @TODO + //for (int i = 0; i < mv_sigma.size(); ++i) { + // mv_X(0) += mv_sigma(i) / mv_sigma.size(); + // mv_X(element_of_site(i)) += mv_sigma(i) / mv_sigma.size(); + // } + + mv_X.setConstant(mv_sigma(0) / (mv_sigma.size() + 1)); // std::cout << "Before X = \n" << mv_X << std::endl; mv_X = solve(mv_X); From 720f801fa430117c2bca9f00d8ccd6b632ebde00 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 29 Jan 2020 14:19:07 +0100 Subject: [PATCH 02/14] added test case nitridation --- src/gsi/SurfaceBalanceSolverMass.cpp | 2 +- tests/data/gsi/smb_FRC_nitridarion.xml | 31 +++ .../smb_FRC_nitridation_NASA9_ChemNonEq1T.xml | 13 ++ tests/test_gsi_detailed.cpp | 196 +++++++++++++++--- 4 files changed, 208 insertions(+), 34 deletions(-) create mode 100644 tests/data/gsi/smb_FRC_nitridarion.xml create mode 100644 tests/data/mixtures/smb_FRC_nitridation_NASA9_ChemNonEq1T.xml 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/tests/data/gsi/smb_FRC_nitridarion.xml b/tests/data/gsi/smb_FRC_nitridarion.xml new file mode 100644 index 00000000..68b6e7ce --- /dev/null +++ b/tests/data/gsi/smb_FRC_nitridarion.xml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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/test_gsi_detailed.cpp b/tests/test_gsi_detailed.cpp index b3b6c4ba..7a48bc62 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,114 @@ 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); + } } @@ -356,4 +486,4 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") } -} \ No newline at end of file +} From e294a0990a915eb0bc2b024b3bdea2bc2a2e0982 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 29 Jan 2020 16:57:07 +0100 Subject: [PATCH 03/14] SEB Detailed Provider --- src/gsi/GSIRateManagerDetailed.cpp | 3 +++ src/gsi/SurfaceBalanceSolverMassEnergy.cpp | 8 ++++++-- tests/CMakeLists.txt | 1 + 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/gsi/GSIRateManagerDetailed.cpp b/src/gsi/GSIRateManagerDetailed.cpp index fbe5bb67..4a1dd6b7 100644 --- a/src/gsi/GSIRateManagerDetailed.cpp +++ b/src/gsi/GSIRateManagerDetailed.cpp @@ -332,5 +332,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/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/tests/CMakeLists.txt b/tests/CMakeLists.txt index f78e999c..0bdcb77f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -34,6 +34,7 @@ add_executable(run_tests test_energies.cpp test_errors.cpp test_gsi_detailed.cpp + test_gsi_detailed_mass_energy.cpp test_gsi_mass.cpp test_gsi_mass_energy.cpp test_gsi_rates.cpp From eddb415b12113eb984aea6f20fbb971ae947c310 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 29 Jan 2020 16:58:31 +0100 Subject: [PATCH 04/14] SEB detailed covergency test --- tests/data/gsi/seb_FRC_nitridarion.xml | 43 ++++++ .../seb_FRC_nitridation_NASA9_ChemNonEq1T.xml | 13 ++ tests/test_gsi_detailed_mass_energy.cpp | 135 ++++++++++++++++++ 3 files changed, 191 insertions(+) create mode 100644 tests/data/gsi/seb_FRC_nitridarion.xml create mode 100644 tests/data/mixtures/seb_FRC_nitridation_NASA9_ChemNonEq1T.xml create mode 100644 tests/test_gsi_detailed_mass_energy.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/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/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)); + } + +} From 2fad1e34217e8b8ad3c44e47ca42071807a883a3 Mon Sep 17 00:00:00 2001 From: capriatim Date: Thu, 1 Oct 2020 15:39:18 +0200 Subject: [PATCH 05/14] fix electron mass for matching the one of US3D --- data/thermo/elements.xml | 4 ++-- src/fortran/cwrapper.cpp | 30 +++++++++++++++++++++++++++--- src/fortran/cwrapper.h | 25 ++++++++++++++++++++++++- src/fortran/cwrapper_interface.f90 | 15 ++++++++++++++- src/thermo/ParticleRRHO.cpp | 5 +++++ 5 files changed, 72 insertions(+), 7 deletions(-) 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..7525efdb 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,4 @@ void NAME_MANGLE(convert_ys_to_ye)( p_mix->convert(species_y, elements_y); } -//============================================================================== \ No newline at end of file +//============================================================================== diff --git a/src/fortran/cwrapper.h b/src/fortran/cwrapper.h index 22c4411e..d21d506f 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 diff --git a/src/fortran/cwrapper_interface.f90 b/src/fortran/cwrapper_interface.f90 index 8bce48da..e79db172 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,18 @@ 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 end interface diff --git a/src/thermo/ParticleRRHO.cpp b/src/thermo/ParticleRRHO.cpp index 91d4b9ca..881b4fd0 100644 --- a/src/thermo/ParticleRRHO.cpp +++ b/src/thermo/ParticleRRHO.cpp @@ -80,6 +80,8 @@ ParticleRRHO::ParticleRRHO(const IO::XmlElement& xml_element) int degeneracy; double temperature; + + // int dummy = 0; for ( ; level_iter != iter->end(); ++level_iter) { if (level_iter->tag() == "level") { @@ -87,8 +89,11 @@ ParticleRRHO::ParticleRRHO(const IO::XmlElement& xml_element) level_iter->getAttribute("energy", temperature); // convert from 1/cm to K + //if (dummy == 0) { m_electronic_energies.push_back( std::make_pair(degeneracy, temperature * 1.4387)); + //dummy = 1; + //} } } } From bec63c93eac312edec696172ddf6934d65285e33 Mon Sep 17 00:00:00 2001 From: capriatim Date: Tue, 27 Oct 2020 13:16:28 +0100 Subject: [PATCH 06/14] part --- src/thermo/ParticleRRHO.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/thermo/ParticleRRHO.cpp b/src/thermo/ParticleRRHO.cpp index 881b4fd0..4c16f872 100644 --- a/src/thermo/ParticleRRHO.cpp +++ b/src/thermo/ParticleRRHO.cpp @@ -80,20 +80,15 @@ ParticleRRHO::ParticleRRHO(const IO::XmlElement& xml_element) int degeneracy; double temperature; - - // int dummy = 0; - + 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 == 0) { m_electronic_energies.push_back( std::make_pair(degeneracy, temperature * 1.4387)); - //dummy = 1; - //} } } } From 4045a28d938bd33cba51efdb295c33e31e410b30 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 3 Feb 2021 14:52:51 +0100 Subject: [PATCH 07/14] wrapper for cabaret --- src/fortran/cwrapper.cpp | 13 +++++++++++++ src/fortran/cwrapper.h | 10 ++++++++++ 2 files changed, 23 insertions(+) diff --git a/src/fortran/cwrapper.cpp b/src/fortran/cwrapper.cpp index 7525efdb..ae48717a 100644 --- a/src/fortran/cwrapper.cpp +++ b/src/fortran/cwrapper.cpp @@ -508,3 +508,16 @@ void NAME_MANGLE(convert_ys_to_ye)( } //============================================================================== +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 d21d506f..5961c809 100644 --- a/src/fortran/cwrapper.h +++ b/src/fortran/cwrapper.h @@ -506,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 From 26c6edf15878876bb28974814ba6f334eee94891 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 3 Feb 2021 14:53:18 +0100 Subject: [PATCH 08/14] wrapper for cabaret --- src/fortran/cwrapper_interface.f90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/fortran/cwrapper_interface.f90 b/src/fortran/cwrapper_interface.f90 index e79db172..844fc307 100644 --- a/src/fortran/cwrapper_interface.f90 +++ b/src/fortran/cwrapper_interface.f90 @@ -126,6 +126,12 @@ real(kind=8) function mpp_vibrational_thermal_conductivity() 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 From 45cfb169cd5f187db5eac2bc03519c633087ea5c Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 3 Feb 2021 14:54:17 +0100 Subject: [PATCH 09/14] preferential CV --- src/transfer/OmegaCV.cpp | 46 +++++++++++++++++++++++++++++++++++++--- src/transfer/OmegaVT.cpp | 2 ++ 2 files changed, 45 insertions(+), 3 deletions(-) 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; } From dcdb4ca3782378011699a9d14aab01029b2bc800 Mon Sep 17 00:00:00 2001 From: capriatim Date: Wed, 3 Feb 2021 14:55:53 +0100 Subject: [PATCH 10/14] PRATA FRC --- src/gsi/CMakeLists.txt | 2 + src/gsi/GSIRateLawHLArrhenius.cpp | 124 ++++++++++++++++ src/gsi/GSIRateLawHLArrhenius_abla.cpp | 119 +++++++++++++++ src/gsi/GSIRateManagerDetailed.cpp | 2 +- src/kinetics/Kinetics.cpp | 1 + src/thermo/ParticleRRHO.cpp | 4 + tests/data/gsi/smb_FRC2_nitridarion.xml | 40 ++++++ ...smb_FRC2_nitridation_NASA9_ChemNonEq1T.xml | 13 ++ tests/test_gsi_detailed.cpp | 135 ++++++++++++++++++ 9 files changed, 439 insertions(+), 1 deletion(-) create mode 100644 src/gsi/GSIRateLawHLArrhenius.cpp create mode 100644 src/gsi/GSIRateLawHLArrhenius_abla.cpp create mode 100644 tests/data/gsi/smb_FRC2_nitridarion.xml create mode 100644 tests/data/mixtures/smb_FRC2_nitridation_NASA9_ChemNonEq1T.xml diff --git a/src/gsi/CMakeLists.txt b/src/gsi/CMakeLists.txt index 4328a9f4..78a83042 100644 --- a/src/gsi/CMakeLists.txt +++ b/src/gsi/CMakeLists.txt @@ -34,6 +34,8 @@ add_sources(mutation++ GSIRateLawGammaConst.cpp GSIRateLawGammaT.cpp GSIRateLawERArrhenius.cpp + GSIRateLawHLArrhenius.cpp + GSIRateLawHLArrhenius_abla.cpp GSIRateLawSublimation.cpp GSIRateManagerDetailed.cpp GSIRateManagerPhenomenological.cpp diff --git a/src/gsi/GSIRateLawHLArrhenius.cpp b/src/gsi/GSIRateLawHLArrhenius.cpp new file mode 100644 index 00000000..1c379edf --- /dev/null +++ b/src/gsi/GSIRateLawHLArrhenius.cpp @@ -0,0 +1,124 @@ +/** + * @file GSIRateLawGammaT.cpp + * + * @brief Class which computes the reaction rate constant for an adsorption + * 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_gas_r(0), + pos_site_r(1) + { + 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 gas in the reactants + m_idx_gas = mv_react[pos_gas_r]; + //std::cout << "m_idx_gas is " << m_idx_gas << std::endl; //it is number 10 + // 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); + } + +//============================================================================== + + ~GSIRateLawLHArrhenius( ){ } + +//============================================================================== + + 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); + + //std::cout << "th speed is" << thermal_speed << std::endl; + double m_n = 0.0140067 / NA; //hardcoded for now + //return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*v_Tsurf[0]/ (8.0*m_n)) * exp(-21000.0/ v_Tsurf[0]); + //return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*v_Tsurf[0]/ (8.0*m_n)) * exp(-m_T_act/ v_Tsurf[0]); + return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*v_Tsurf[0]/ (2.0*m_n)) * exp(-m_T_act/ v_Tsurf[0]); + + //return m_pre_exp * thermal_speed * sqrt(1.0/m_n_sites)/ + // (4.) * 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< + 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..ae751a59 --- /dev/null +++ b/src/gsi/GSIRateLawHLArrhenius_abla.cpp @@ -0,0 +1,119 @@ +/** + * @file GSIRateLawGammaT.cpp + * + * @brief Class which computes the reaction rate constant for an adsorption + * 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), + 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() == "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."); + + // 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); + } + +//============================================================================== + + ~GSIRateLawLHAblaArrhenius( ){ } + +//============================================================================== + + 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 * sqrt(NA/m_n_sites)/ + // (4.) * exp(-m_T_act / Tsurf); + + return m_pre_exp*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< + 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 4a1dd6b7..6f79f9fa 100644 --- a/src/gsi/GSIRateManagerDetailed.cpp +++ b/src/gsi/GSIRateManagerDetailed.cpp @@ -91,7 +91,7 @@ class GSIRateManagerDetailed : } // Setup NewtonSolver - setMaxIterations(5); + setMaxIterations(50); setWriteConvergenceHistory(false); setEpsilon(m_tol); } 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 4c16f872..e1386358 100644 --- a/src/thermo/ParticleRRHO.cpp +++ b/src/thermo/ParticleRRHO.cpp @@ -80,6 +80,7 @@ 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") { @@ -87,8 +88,11 @@ ParticleRRHO::ParticleRRHO(const IO::XmlElement& xml_element) 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/tests/data/gsi/smb_FRC2_nitridarion.xml b/tests/data/gsi/smb_FRC2_nitridarion.xml new file mode 100644 index 00000000..aeb8cd14 --- /dev/null +++ b/tests/data/gsi/smb_FRC2_nitridarion.xml @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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/test_gsi_detailed.cpp b/tests/test_gsi_detailed.cpp index 7a48bc62..3539cb1b 100644 --- a/tests/test_gsi_detailed.cpp +++ b/tests/test_gsi_detailed.cpp @@ -486,4 +486,139 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") } + + 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(); + + CHECK(mixFRC.getSurfaceProperties().isSurfaceCoverageSteady() == true); + ArrayXd v_surf_cov_mpp_frac(mixFRC.getSurfaceProperties().nSurfaceSpecies()); + ArrayXd v_surf_cov_frac(mixFRC.getSurfaceProperties().nSurfaceSpecies()); + + // Equilibrium Surface, problema a P = 1.e-5 e T = 1000.0 @TODO + double P = 1.e-4; + double dP = 10.; + double T; // K + double dT = 200.; // K + for (int i = 0; i < 14; i++) { + for (int j = 0; j < 12; j++) { + T = (j+6) * dT; + //std::cout << T << " " << P << std::endl; + //T = 2407.; + //P = 1500.; + + 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(-2500./T); + double kfdes = 2.0 * PI * mm(iN) / NA * KB * KB * T * T / (HP * HP * HP) / B; + kfdes *= exp(-73971.6/T); + double kfer1 = F * 1.5 * exp(-7000./T); + double kfer2 = F * 0.5 * exp(-2000./T); + //double kflh1 = 4.96155568504e-12; //sqrt(1.0/(B)) * sqrt(PI*KB*T/ (8.0 * mm(iN)))* 0.1 * exp(-21000.0/ T); + //double kflh1 = sqrt(NA/(B)) * sqrt(PI*KB*T/ (8.0 * mm(iN)))* 0.1 * exp(-21000.0/ T); + double kflh1 = sqrt(NA/(B)) * sqrt(PI*KB*T/ (2.0 * mm(iN)))* 0.1 * exp(-21000.0/ T); + double kflh2 = 1e8 * exp(-20676.0 / T); + + //std::cout << "kflh1 is " << kflh1 << std::endl; + + 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 << v_surf_cov_frac(0) << " " << v_surf_cov_frac(1)<< std::endl; + //std::cout << v_surf_cov_mpp_frac(0) << " " << v_surf_cov_mpp_frac(1)<< 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)); + 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); + rates(4) = kflh1 * v_surf_cov_frac(1)* v_surf_cov_frac(1); + rates(5) = kflh2 * v_surf_cov_frac(1); + + /*std::cout << "HERE" << 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 ) { + 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) + rates(4)); + wdot(2) = 0.; + wdot(3) = -mm(iCN) / NA * (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); + } + + } + + P *= dP; + } + + } + } From e4c1553fae079d99bbdf56298d4370f3b8d031d5 Mon Sep 17 00:00:00 2001 From: capriatim Date: Fri, 16 Aug 2024 13:35:30 +0200 Subject: [PATCH 11/14] fix in catch_fatal_condition.hpp --- thirdparty/catch/include/internal/catch_fatal_condition.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 From b33d2308174bdafd914de7087b7e1fa0514edd8e Mon Sep 17 00:00:00 2001 From: capriatim Date: Fri, 16 Aug 2024 14:16:17 +0200 Subject: [PATCH 12/14] full ACA model implementation --- .gitignore | 2 + data/gsi/smb_ACA_model.xml | 111 +++++++++++++ src/gsi/GSIRateLawER2Arrhenius.cpp | 117 +++++++++++++ src/gsi/GSIRateLawERArrhenius.cpp | 2 +- src/gsi/GSIRateLawHLArrhenius.cpp | 36 ++-- src/gsi/GSIRateLawHLArrhenius_abla.cpp | 38 +---- src/gsi/GSIRateManagerDetailed.cpp | 29 +--- src/gsi/SurfacePropertiesDetailed.cpp | 10 +- tests/CMakeLists.txt | 1 + tests/data/gsi/smb_FRC2_nitridarion.xml | 2 +- tests/data/gsi/smb_FRC_nitridarion.xml | 31 ---- tests/data/gsi/smb_full_ACA.xml | 111 +++++++++++++ .../smb_full_ACA_NASA9_ChemNonEq1T.xml | 13 ++ tests/test_gsi_detailed.cpp | 72 ++++---- tests/test_gsi_detailed_ACA.cpp | 156 ++++++++++++++++++ 15 files changed, 576 insertions(+), 155 deletions(-) create mode 100644 data/gsi/smb_ACA_model.xml create mode 100644 src/gsi/GSIRateLawER2Arrhenius.cpp delete mode 100644 tests/data/gsi/smb_FRC_nitridarion.xml create mode 100644 tests/data/gsi/smb_full_ACA.xml create mode 100644 tests/data/mixtures/smb_full_ACA_NASA9_ChemNonEq1T.xml create mode 100644 tests/test_gsi_detailed_ACA.cpp 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/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 index 1c379edf..35ce40cf 100644 --- a/src/gsi/GSIRateLawHLArrhenius.cpp +++ b/src/gsi/GSIRateLawHLArrhenius.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 HL * surface reaction based on an Arrhenius formula. */ @@ -49,8 +49,7 @@ class GSIRateLawLHArrhenius : public GSIRateLaw m_surf_props(args.s_surf_props), mv_react(args.s_reactants), pos_T_trans(0), - pos_gas_r(0), - pos_site_r(1) + pos_site_r(0) { assert(args.s_node_rate_law.tag() == "LH_arrhenius"); @@ -61,17 +60,14 @@ class GSIRateLawLHArrhenius : public GSIRateLaw "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]; - //std::cout << "m_idx_gas is " << m_idx_gas << std::endl; //it is number 10 - // Error if m_idx_gas > ns - // For the sites - int idx_site = mv_react[pos_site_r]; - // Error if idx_site > ns + m_idx_site = mv_react[pos_site_r]; - m_site_categ = m_surf_props.siteSpeciesToSiteCategoryIndex(idx_site); + 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; } //============================================================================== @@ -85,32 +81,22 @@ class GSIRateLawLHArrhenius : public GSIRateLaw { 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); - - //std::cout << "th speed is" << thermal_speed << std::endl; - double m_n = 0.0140067 / NA; //hardcoded for now - //return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*v_Tsurf[0]/ (8.0*m_n)) * exp(-21000.0/ v_Tsurf[0]); - //return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*v_Tsurf[0]/ (8.0*m_n)) * exp(-m_T_act/ v_Tsurf[0]); - return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*v_Tsurf[0]/ (2.0*m_n)) * exp(-m_T_act/ v_Tsurf[0]); + return m_pre_exp *sqrt(1.0/m_n_sites) * sqrt(PI*KB*Tsurf/ (2.0*m_mw_r)) * exp(-m_T_act/ Tsurf); - //return m_pre_exp * thermal_speed * sqrt(1.0/m_n_sites)/ - // (4.) * 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_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; diff --git a/src/gsi/GSIRateLawHLArrhenius_abla.cpp b/src/gsi/GSIRateLawHLArrhenius_abla.cpp index ae751a59..cf5c3c47 100644 --- a/src/gsi/GSIRateLawHLArrhenius_abla.cpp +++ b/src/gsi/GSIRateLawHLArrhenius_abla.cpp @@ -1,8 +1,8 @@ /** * @file GSIRateLawGammaT.cpp * - * @brief Class which computes the reaction rate constant for an adsorption - * surface reaction based on an Arrhenius formula. + * @brief Class which computes the reaction rate constant for a gas + * independent ablation surface reaction based on an Arrhenius formula. */ /* @@ -46,11 +46,7 @@ class GSIRateLawLHAblaArrhenius : public GSIRateLaw public: GSIRateLawLHAblaArrhenius(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) + pos_T_trans(0) { assert(args.s_node_rate_law.tag() == "LH_abla_arrhenius"); @@ -60,17 +56,6 @@ class GSIRateLawLHAblaArrhenius : public GSIRateLaw 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); } //============================================================================== @@ -84,31 +69,14 @@ class GSIRateLawLHAblaArrhenius : public GSIRateLaw { 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 * sqrt(NA/m_n_sites)/ - // (4.) * exp(-m_T_act / Tsurf); - return m_pre_exp*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< diff --git a/src/gsi/GSIRateManagerDetailed.cpp b/src/gsi/GSIRateManagerDetailed.cpp index 6f79f9fa..a4ef532d 100644 --- a/src/gsi/GSIRateManagerDetailed.cpp +++ b/src/gsi/GSIRateManagerDetailed.cpp @@ -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,30 +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.e17); - //mv_X.setConstant(3.011e18); - - //mv_X.setZero(); @TODO - //for (int i = 0; i < mv_sigma.size(); ++i) { - // mv_X(0) += mv_sigma(i) / mv_sigma.size(); - // mv_X(element_of_site(i)) += mv_sigma(i) / mv_sigma.size(); - // } - - mv_X.setConstant(mv_sigma(0) / (mv_sigma.size() + 1)); + 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. @@ -290,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 } //============================================================================= 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/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0bdcb77f..0c3ccf6b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -34,6 +34,7 @@ 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 diff --git a/tests/data/gsi/smb_FRC2_nitridarion.xml b/tests/data/gsi/smb_FRC2_nitridarion.xml index aeb8cd14..e966bac2 100644 --- a/tests/data/gsi/smb_FRC2_nitridarion.xml +++ b/tests/data/gsi/smb_FRC2_nitridarion.xml @@ -33,7 +33,7 @@ - + diff --git a/tests/data/gsi/smb_FRC_nitridarion.xml b/tests/data/gsi/smb_FRC_nitridarion.xml deleted file mode 100644 index 68b6e7ce..00000000 --- a/tests/data/gsi/smb_FRC_nitridarion.xml +++ /dev/null @@ -1,31 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 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/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 3539cb1b..553d0701 100644 --- a/tests/test_gsi_detailed.cpp +++ b/tests/test_gsi_detailed.cpp @@ -380,7 +380,7 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") } - SECTION("Nitridation Model.") +/* SECTION("Nitridation Model.") { // Setting up M++ MixtureOptions optsFRC("smb_FRC_nitridation_NASA9_ChemNonEq1T"); @@ -484,7 +484,7 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") P *= dP; } - } + }*/ SECTION("Nitridation Model 2.") @@ -513,22 +513,23 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") 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, problema a P = 1.e-5 e T = 1000.0 @TODO - double P = 1.e-4; + // 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 < 14; i++) { + for (int i = 0; i < 15; i++) { for (int j = 0; j < 12; j++) { T = (j+6) * dT; - //std::cout << T << " " << P << std::endl; - //T = 2407.; - //P = 1500.; + //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()); @@ -539,19 +540,16 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") 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(-2500./T); - double kfdes = 2.0 * PI * mm(iN) / NA * KB * KB * T * T / (HP * HP * HP) / B; - kfdes *= exp(-73971.6/T); - double kfer1 = F * 1.5 * exp(-7000./T); - double kfer2 = F * 0.5 * exp(-2000./T); - //double kflh1 = 4.96155568504e-12; //sqrt(1.0/(B)) * sqrt(PI*KB*T/ (8.0 * mm(iN)))* 0.1 * exp(-21000.0/ T); - //double kflh1 = sqrt(NA/(B)) * sqrt(PI*KB*T/ (8.0 * mm(iN)))* 0.1 * exp(-21000.0/ T); - double kflh1 = sqrt(NA/(B)) * sqrt(PI*KB*T/ (2.0 * mm(iN)))* 0.1 * exp(-21000.0/ T); - double kflh2 = 1e8 * exp(-20676.0 / T); - - //std::cout << "kflh1 is " << kflh1 << std::endl; + 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; @@ -566,17 +564,14 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") v_surf_cov_mpp_frac = mixFRC.getSurfaceProperties().getSurfaceSiteCoverageFrac(); v_surf_cov_mpp_frac *= B; - //std::cout << v_surf_cov_frac(0) << " " << v_surf_cov_frac(1)<< std::endl; - //std::cout << v_surf_cov_mpp_frac(0) << " " << v_surf_cov_mpp_frac(1)<< std::endl; + //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)); - 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((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); @@ -586,26 +581,23 @@ TEST_CASE("Detailed surface chemictry tests.","[gsi]") rates(4) = kflh1 * v_surf_cov_frac(1)* v_surf_cov_frac(1); rates(5) = kflh2 * v_surf_cov_frac(1); - /*std::cout << "HERE" << 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;*/ + //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 ) { - 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) + rates(4)); + 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) / NA * (rates(2) + rates(5)); + wdot(3) = -mm(iCN) * (rates(2) + rates(5)); for (int ii = 0; ii < 4; ++ii ) { if (abs(wdotmpp(ii)) >= 1e-14) 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(); + + + + } + + +} From 3b0270d3ba801981cbfa710e6d9168c9cb9932e5 Mon Sep 17 00:00:00 2001 From: capriatim Date: Fri, 16 Aug 2024 14:17:13 +0200 Subject: [PATCH 13/14] full ACA model implementation --- src/gsi/GSIRateLawDesorptionTST2.cpp | 105 --------------------------- 1 file changed, 105 deletions(-) delete mode 100644 src/gsi/GSIRateLawDesorptionTST2.cpp diff --git a/src/gsi/GSIRateLawDesorptionTST2.cpp b/src/gsi/GSIRateLawDesorptionTST2.cpp deleted file mode 100644 index f2b514bf..00000000 --- a/src/gsi/GSIRateLawDesorptionTST2.cpp +++ /dev/null @@ -1,105 +0,0 @@ -/** - * @file GSIRateLawGammaT.cpp - * - * @brief Class which computes the reaction rate constant for a desorption - * surface reaction based on simple transition state theory. - */ - -/* - * 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 GSIRateLawDesorptionTST2 : public GSIRateLaw -{ -public: - GSIRateLawDesorptionTST2(ARGS args) - : GSIRateLaw(args), - m_surf_props(args.s_surf_props), - mv_react(args.s_reactants), - pos_T_trans(0), - idx_react(0) - { - assert(args.s_node_rate_law.tag() == "desorption_tst2"); - - args.s_node_rate_law.getAttribute("T", m_T_des, - "Activation temperature should be provided for every desorption " - "reaction."); - - m_mass_des = m_thermo.speciesMw( - args.s_surf_props.surfaceToGasIndex(mv_react[idx_react])) / NA; - m_site_categ = m_surf_props.siteSpeciesToSiteCategoryIndex( - mv_react[idx_react]); - m_n_sites = m_surf_props.nSiteDensityInCategory(m_site_categ); - // m_n_sites = args.s_surf_props.nSiteDensity(); - } - -//============================================================================== - - ~GSIRateLawDesorptionTST2( ){ } - -//============================================================================== - - double forwardReactionRateCoefficient( - const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_Tsurf) const - { - const double Tsurf = v_Tsurf(pos_T_trans); - - const double pre_exp = 2 * PI * m_mass_des * KB * KB * Tsurf - / (HP * HP * HP * m_n_sites); - - return pre_exp * exp(-m_T_des / Tsurf); - } - -private: - const size_t idx_react; - int m_site_categ; - - double m_T_des; - - double m_n_sites; - double m_mass_des; - - const double pos_T_trans; - - const std::vector& mv_react; - const SurfaceProperties& m_surf_props; -}; - -ObjectProvider< - GSIRateLawDesorptionTST2, GSIRateLaw> - gsi_rate_law_desorption_tst2("desorption_tst2"); - - } // namespace GasSurfaceInteraction -} // namespace Mutation From 58fe18279c71ec6b4d255dd311667ea5a0a07d35 Mon Sep 17 00:00:00 2001 From: capriatim Date: Fri, 16 Aug 2024 14:18:09 +0200 Subject: [PATCH 14/14] full ACA model implementation --- src/gsi/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/CMakeLists.txt b/src/gsi/CMakeLists.txt index 78a83042..f2a864b8 100644 --- a/src/gsi/CMakeLists.txt +++ b/src/gsi/CMakeLists.txt @@ -30,10 +30,10 @@ add_sources(mutation++ GSIReactionDetailed.cpp GSIRateLawAdsorption.cpp GSIRateLawDesorptionTST.cpp - GSIRateLawDesorptionTST2.cpp GSIRateLawGammaConst.cpp GSIRateLawGammaT.cpp GSIRateLawERArrhenius.cpp + GSIRateLawER2Arrhenius.cpp GSIRateLawHLArrhenius.cpp GSIRateLawHLArrhenius_abla.cpp GSIRateLawSublimation.cpp