diff --git a/docs/input-files.md b/docs/input-files.md index acb659bb..ccb1eab1 100644 --- a/docs/input-files.md +++ b/docs/input-files.md @@ -20,7 +20,16 @@ - [Unit Specifiers](#unit-specifiers) - [Example Mechanism](#example-mechanism) - [Collision Integral Database](#collision-integrals) (work in progress) - +- [Gas-Surface Interaction](#gsi) + - [Surface Chemical Production Terms](#surf_chem) + - [Carbon Oxidation Model](#coxid) + - [Carbon Nitridation Model](#cnitr) + - [Carbon Sublimation Model](#csubl) + - [Surface Mass Balance](#smb) + - [Surface Mass Balance for Catalytic Surfaces](#smbc) + - [Surface Mass Balance for Ablative and Catalytic Surfaces](#smbac) + - [Surface Energy Balance](#seb) + - [Using the GSI module](#gsi_example) ## Introduction @@ -33,6 +42,7 @@ The data files distributed with Mutation++ are located in one of the subdirector + __thermo__ [elemental](#elements) and species [thermodynamic databases](#thermodynamic-databases) + __transfer__ internal energy [transfer model databases](#transfer-databases) + __transport__ @subpage collisions "collision integral database" + + __gsi__ gas-surface interaction mechanism ### File location @@ -424,7 +434,7 @@ controlled by an Arrhenius rate law. \# | Formula | A [mol,cm,s,K] | n | Ea [K] ---|------------------------------------------|----------------|-------|-------- -1 | N2 + N_2 <-> 2N + N2 | 1.0e21 | -1.6 | 113200 +1 | N2 + N2 <-> 2N + N2 | 1.0e21 | -1.6 | 113200 2 | N2 + N <-> 2N + N | 3.0e21 | -1.6 | 113200 3 | N2 + N^+ <-> 2N + N+ | 1.0e21 | -1.6 | 113200 4 | N2 + e- <-> 2N + e- | 7.0e22 | -1.6 | 113200 @@ -949,3 +959,542 @@ example below. The value of the effective electrons for new species can be estimated using the approach of Cambi et al. \cite Cambi1991. +## Gas-Surface Interaction + + +The gas-surface interaction module of Mutation++ is responsinble for +treating chemically reacting surfaces in thermochemical non-equilibrium. +With the aim to obtain the appropriate surface boundary conditions +for a chemically reacting gas, the conservation equations on the +interface have to be written. In general, this idea can be applied for several +different categories of interfaces with different phases and materials on each +side. Since this Mutation++ module emphases on gas-surface interaction it is limited +to having gas on one side (underscore g) and a catalytic or ablative solid or liquid +on the other (underscore b). The surface is not simulated, but approximately modeled. +In most of the cases it can be considered +impermeable, or simplistically porous for pyrolysis gases. The interface +can promote chemical reactions or radiate following the Planck's law of radiation +at the surface temperature. The set of balance equations obtained +are solved with respect to the conserved quantities of the gas, i.e. density, +momentum and energy, which are imposed as boundary values for the full Navier-Stokes equations, +or the porous material. + +The general procedure to obtain the surface balances is the following. +Assuming steady state on the surface, the time derivative term +is equal to zero. By taking the limit that one dimension of the volume goes to zero, the +three dimensional fluxes reduce to the normal fluxes at the interface while the volume +source terms, such as chemical reactions, go to zero; only the surface sources remain. +Generally, the set of balance equations can be written as: + +(**F**_g - **F**_b) . **n** = Ω_s, + +where **F** are the fluxes from the gas (**g**) and the bulk (**b**) phases, while +Ω are the source terms associated to the surface (s) processes. +Only the normal to surface flux component should be considered, denoted by +the inner product of the flux with the surface unit vector, (**n**). +The specific form of the balance equations for mass and energy will be seen in the following +sections. + +### Surface Chemical Production Terms + + +Below the input file when a single catalytic reactions is considered +modeled with the gamma γ model with reaction probability equal to 1. + +```xml + + + + + + + + + O:1. + + + + +``` + +The γ model, introduced by Goulard in the late 50s, is arguably the most popular +way to treat catalysis in the aero-thermodynamics community. It describes catalytic +reactions as macroscopic, non-elementary processes of the form: + +A + A -> A2 + +In order to determine the chemical production term for this kind of catalytic reactions a +probability for recombination γ is defined for each recombining species A as: + + +γ = F_rec / F_imp, + + +where F_imp is the flux of species A impinging the surface and +F_rec is the flux of species recombining at the surface. This probability is +the input parameter for the model. A fully catalytic wall has γ equal to 1, which +means that all the particles of species A impinging the surface recombine at the wall. +γ equal to 0 means that no reaction takes place and corresponds +to a non catalytic, or chemically inert wall. Anything between these two extreme cases is a +partially catalytic wall, which is the case for most of the surfaces. When the probability +γ is defined, the surface chemical source term is determined as: + +Ω_A = γ m F_imp + +Ω_A2 = - γ m F_imp + +with m being the mass of species A. + +The only parameter that still needs to be defined is the impinging flux on the surface. +When the distribution function of species at the wall is well approximated by a Maxwellian +and there is no temperature slip, the impinging flux, F_imp, is equal to + +F_imp = n_A (k_B T / 2 π m_A)^(1/2) + +according to the kinetic theory of gases. + +Initially, the gamma model was developed to describe homonuclear reactions, +such as the one presented in above. +Soon, though, it was observed that heteronuclear reactions were also probable +to occur at the wall, which can take the form of + +A + B -> AB, + +with the chemical rate for the produced molecule AB equal to: + +Ω_AB = - γ_AB m_A F_imp_A - γ_BA m_B F_imp_B. + +One atom of A recombines with one atom of B on the surface to produce a molecule AB. +In other words the number of atoms of A that recombine into AB should be equal to the number of B +atoms that recombine into AB. This restriction should be explicitly imposed in order to conserve mass: + +γ_AB F_imp_A = γ_BA F_imp_B. + +As a result, in a reaction of this type two gamma recombination coefficients should be defined not +necessarily equal for the two processes, the one activated when the catalytic reaction is limited by the +flux of A atoms and the opposite. In practice the two recombination number fluxes, +γ_AB F_imp_A and γ_BA F_imp_B, are compared and the limiting one +determines which of the two gammas is chosen. + +These gamma coefficients cannot take arbitrary values, they should be limited between 0 and 1, +just like in the homonuclear case. When extra catalytic reactions are added, +such as A + A -> A2 and B + B -> B2 the gammas should be further constrained +as 0 < γ_A + γ_AB < 1 and 0 < γ_B + γ_BA < 1 +in order for mass to be conserved. This approach is consistent with similar approaches +considering the catalytic recombination occurring in Martian atmospheres, where O can recombine +into both O2 and CO2 due to catalytic reactions. + +An example input file of catalytic reactions in air including the formation of NO can be seen below. + +```xml + + + + + + + + + N:.001 + + + + O:.001 + + + + N:.002 O:.003 + + + + +``` + +It is still unclear which is the proper boundary conditions for the +ions and electrons on the surface. One of the most common model is to +assume full ion recombination on the surface can be expressed as +catalytic reaction with probability 1. +The following example shows how to impose full ion recombination in +Mutation. + +```xml + + + + + + + + + N:.001 + + + + O:.001 + + + + N:.002 O:.003 + + + + N+:1. e-:1. + + + + O+:1. e-:1. + + + + NO+:1. e-:1. + + + + N2+:1. e-:1. + + + + O2+:1. e-:1. + + + + +``` +An example of ablation model can be seen below. +```xml + + + + + + + + + + + + + + N:0.003 + + + + + + + + N:0.001 + + + + +``` + +#### Carbon Oxidation Model + + +The first ablation reaction presented above is the oxidation +of the solid carbon by atomic oxygen. The reaction reads as: + +C_b + O -> CO + +and is exothermic, releasing 3.74 eV per molecule produced. Its reaction +rate coefficient is given by defining a recombination probability +γ_CO. This probability is an Arrhenius type function of temperature +and is given by the formula: + +γ_CO = 0.63 exp(-1160 / T). + +Carbon oxidation with molecular oxygen is also possible +C_b + O2 -> CO + O, but it is often considered as a less significant +process. + +#### Carbon Nitridation Model + + +The second shown in the example is carbon nitridation, + +C_b + N -> CN, + +an exothermic reaction with 0.35 eV of energy released per reacting atom. +The reaction rate is given using a constant recombination probability +γ_CN like in the catalytic case. + +#### Carbon Sublimation Model + + +At high temperatures carbon removal from the surface is dominated by phase +change processes like sublimation. The production of C3 is considered here. +It should be noted that this type of reactions are invertible, with formula: + +3C_b -> C3. + +The chemical production rate of this reaction is equal to: + +Ω_C3 = β (ρ_eq_C3 - ρ_C3) (k_B T / 2 π m)^(1/2). + +The equilibrium partial density of the C3 species is obtained from the +saturated vapor pressure of carbon, which is equal to + +p_sat_C3 = c \exp (-Ta / T) + +with β_C3 being the evaporation coefficient, c the pre-exponential coefficient, +and Ta the activation temperature. +Even though here only sublimation is presented in the example, evaporation processes +can be considered with the same model. + +### Surface Mass Balance + + +#### Surface Mass Balance for Catalytic Surfaces + + +Heterogeneous catalysis is an important gas-surface interaction phenomenon occurring +during re-entry of vehicles equipped with re-usable thermal protection system. It +describes the recombination of the dissociated atoms in the flow using the thermal +protection system as a catalyst. It is called heterogeneous, because the recombining +species and the catalyst are in a different phase, here gas and solid. The catalyst, +without being consumed, increases the rate of chemical reactions by offering an alternative, +energetically favored path. It is important to note that catalysis does not change the chemical +equilibrium of reactions, since it favors equally both the forward and the backward reaction rates. +This is a constrain that should be respected by the catalytic model chosen. + +In hypersonic the recombination reactions that occur on the surface are in general exothermic. +The energy released to the wall is a substantial percentage of the total heat flux experienced +by space vehicles. Not necessarily all of the recombination energy is directly deposited to the +surface. A part of it is used to excite the internal energy of the produced molecules. Another +reason why calculating +the actual heat released on the surface is a complicated task, is that the gas phase chemistry +and diffusion play an important role in determining the overall catalytic rates. +If all of the phenomena above are modeled with accuracy, the re-entry heat load can be predicted +and the size of the thermal protection system can be determined. +The mass balance on a catalytic surface reads: + +**j**_i . **n** = Ω_cat_i, + +where one mass balance equation should be solved for each distinct species in the flow. +The equation above states, that the catalytic activity of every species is equal to the diffusion +flux of these species to the surface. This leads as to two cases. In the first one, the rate +with which the chemical species are produced or destroyed at the wall is higher than the rate +they diffuse to the surface while in the second one the opposite happens. The first case is called +diffusion limited, since diffusion is the mechanism controlling the chemical process, while the +opposite is called reaction limited. When a species i is inert at the surface, then its chemical +rate is exactly equal to zero, which at steady state imposes that its net diffusion flux is also zero. +Even though, in principle, these equations could be omitted, since they impose that the mole fractions +of the species in question do not change with respect to the ones in the gas phase, the full system was +chosen to be solved. + +#### Surface Mass Balance for Ablative and Catalytic Surfaces + + +Ablation is the chemical gas-surface interaction phenomenon which occurs on non-reusable thermal +protection systems of re-entry vehicle. The term ablation describes the category of chemical +reactions during which the dissociated atoms in the flow field recombine directly with the +thermal protection system, which by burning protects the vehicle. Contrary to catalysis, this +burning destroys the material itself making it unable for reuse. Not only chemical reactions can +cause the degradation of ablative thermal protection systems. Mechanical removal processes, such +as spillation, can also occur. The particle injected in the flow field due to these phenomena are +not necessarily in a gaseous form and their modeling requires approaches beyond the scope of this +work. The type of chemical reactions studied here are assumed to produce only gaseous species +and occur only the surface of the material. In cases where the material is porous, ablation +processes also in the bulk, such as in the case of pyrolysis. + +Taking these ideas in mind, the surface mass balance accounting only for surface reactions becomes: + +ρ_i (**u**_g - **u**_r) + **j**_i - **F**_b_i ] . **n** = Ω_i, + +with Ω_i = Ω_cat_i + Ω_abl_i and +term **F**_b_i is the flux of species i entering the interface due +to solid process, like pyrolysis and solid-solid chemical reactions. +Just like before, one mass balance equation should be solved for each distinct species in the flow. +It is most of the time reasonable to consider that recession velocity is +orders of magnitude lower than the gas velocity **u**_r. + +### Surface Energy Balance + + +In order to determine the surface temperature a surface energy balance should be solved along +with the mass balances. It takes the form: + +( ρ ( **u**_g - **u**_r ) H + **q**g + - **F**_b_e ) . **n** = Ω_e; + +where **q**_g is the heat flux to the gas phase equal to: + +**q**_g = -λ ∇ T + Σ **j**_i h_i. + +The radiative heat flux can be also included and will be discussed along with the surface +radiation. + +Term **F**_b_e describes the energy exchanged between the interface +and the bulk of the solid and is composed of three contributions: the first one is the +thermal conduction exiting the surface **q**_cond, the second one is the +enthalpy entering the interface due to the movement of the surface with the recession velocity, +**u**_r ρ h_s and the third one appears only in cases of +porous material, describing the enthalpy of the solid pyrolysis gases convected in the interface, +denoted as **u**_p ρ_p; h_p. The subscript p symbolizes +the pyrolysis gas properties, with the ρ_p h_p being actually the +sum Σ = ρ_i h_i for the pyrolysis gas densities. The surface +enthalpy h_s is an input to the code, as an attribute to the surface_properties element. + +An example input file for solving both mass and energy balance can be seen below. +```xml + + + + + + + + + + + + + + + + + + N:0.003 + + + + + + + + N:0.001 + + + + + +``` + +The only source term taken into account for the surface energy balance in the +example above is radiation. It is considered by adding a new element with tag +. The surface is assumed to be in thermodynamic equilibrium +at a temperature T emitting energy following the Stefan-Boltzmann law, + +Ω_e = Ω_rad = ε ( σ T^4 - **q**_rad_g) + +with σ being the Stefan-Boltzmann constant and ε the emissivity +of the surface. T_env is the surrounding environment temperature which can +be used to simulate the far field radiative heat flux based on the formula + σ T_env^4, which replaces the term **q**_rad_g. +This is a relatively bad approximation and, therefore, it is better to be omitted. + +In order to obtain the value for the conductive heat flux on the solid, +either a material code should be used or one should at least solve the +energy equation in the solid. Even though these approaches can be very accurate, +in cases where the material has low thermal conductivity or the recession rates +are high, approximate methods can be used, without compromising the accuracy of +the simulations. Such an approximation was adopted here for the modeling of the +conductive heat flux inside the material, the steady-state ablation approach for a +semi-infinite surface. By writing the steady-state energy equation for the solid phase +and integrating over the semi-infinite material, where on one side is the solid +properties at exactly the interface (subscript s), while on the other, at infinity, +is the virgin material (subscript v), the steady state heat flux is given by the formula: + +**q**_cond_ss = -**u**_r (ρ_v h_v - ρ_s h_s). + +By replacing the formula above in the surface energy balance and after simplifications one gets: + +(ρ ( **u**_g - **u**_r ) H + **q**_g + **u**_r ρ_v h_v) . **n** = Ω_e, + +which is equally valid for both porous and non-porous materials. In order to use the steady state ablation +approximation, the attribute of surface_feature element solid_conduction should be set to steady state. +Instead of inputing the +virgin material density, the ratio between the virgin and surface density minus one is often used, +refered to as +φ = (ρ_v / ρ_s - 1) and is the attribute +virgin_to_surf_density_ratio in the solid_properties element defaulting to 1. +The enthalpy_virgin can also be an input with default +value equal to 0. Note that these two last options are only necessary when the steady +state assumption for the solid conduction is considered. + +When a material solver is available the conductive heat flux should be an input to the +library for increased accuracy. This can be achieved by setting the solid_conduction +surface feature to "input" and the setSolidCondHeatFlux function can be invoked. +The option enthalpy_surface should be set in this case, otherwise it is automatically set +to zero. + +Note that only one surface balance equation is solved regardless of the thermodynamic +state model. When multitemperature models are considered for the thermodynamics additional +temperatures should be imposed on the surface. It is often a reasonable assumption to +impose thermal equilibrium at the wall. This is achieved by setting the surface_feature +option surface_in_thermal_equil "true". If it is "false", then any additional temperature +beyond translations will be left unchanged, an assumption which can be used to impose an +adiabatic boundary condition for the internal energy modes. + +### Using the GSI module + + +In order for Mutation++ to take into account the Gas Surface Interaction +features the gsi_mechanism attribute should be assigned to the gsi +mechanism file. + +```xml + + + + + ... + + +``` + +Below an example code of how to use the gas surface interaction features +of Mutation is presented. + +```cpp + const int set_state_with_rhoi_T = 1; + + MixtureOptions opts("mixture_name"); + Mixture mix(opts); + + const int iter = 5; + mix.setIterationsSurfaceBalance(iter); + + // Setting the state and setting up the library + mix.setState(rhoi_surf.data(), Tsurf.data(), set_state_with_rhoi_T); + mix.setSurfaceState(rhoiw.data(), Tsurf.data(), set_state_with_rhoi_T); + mix.setDiffusionModel(xi_edge.data(), dx); + mix.setGasFourierHeatFluxModel(Tedge.data(), dx); // Only works with energy balance + // The .data() function returns the pointer at the first element of the data container. + + // Additional Options + double gas_rad = 0.; + setGasRadHeatFlux(*gas_rad); // Only called if feature gas_radiation is true + double solid_cond = 0.; + setSolidCondHeatFlux(*solid_cond); // Only called if feature solid_cond is set to input + + // Solving the surface Mass Balance and requesting the solution + mix.solveSurfaceBalance(); + mix.getSurfaceState(rhoiw.data(), Tsurf.data(), set_state_with_rhoi_T); + + // Getting mass blowing rate. + double mblow; + mix.getMassBlowingRate(mblow); + + // Getting surface reaction rates. + mix.getSurfaceReactionRates(wdot.data()); + + // Getting number of reactions and surface reaction rates per reaction. + int m_nr = mix.nSurfaceReactions(); + mix.getSurfaceReactionRatesPerReaction(wdot_reac.data()); +``` + +Note that this example is not supposed to compile or run, but is there to indicate +the most important features of the library. + + diff --git a/src/gsi/CMakeLists.txt b/src/gsi/CMakeLists.txt index 73c0260f..70bccf7e 100644 --- a/src/gsi/CMakeLists.txt +++ b/src/gsi/CMakeLists.txt @@ -38,8 +38,10 @@ add_sources(mutation++ SolidPropertiesSteadyState.cpp SurfaceChemistry.cpp SurfaceRadiation.cpp + SurfaceInelastic.cpp SurfaceBalanceSolverMass.cpp SurfaceBalanceSolverMassEnergy.cpp + SurfaceBalanceSolverMassEnergyTTv.cpp SurfacePropertiesNull.cpp SurfacePropertiesAblation.cpp SurfaceState.cpp @@ -59,3 +61,4 @@ install(FILES SurfaceChemistry.h DESTINATION include/mutation++) install(FILES SurfaceProperties.h DESTINATION include/mutation++) install(FILES SurfaceRadiation.h DESTINATION include/mutation++) install(FILES SurfaceState.h DESTINATION include/mutation++) +#install(FILES Surfacinelastic.h DESTINATION include/mutation++) diff --git a/src/gsi/GSIRateLawGammaTTv.cpp b/src/gsi/GSIRateLawGammaTTv.cpp new file mode 100644 index 00000000..6325b315 --- /dev/null +++ b/src/gsi/GSIRateLawGammaTTv.cpp @@ -0,0 +1,101 @@ +/** + * @file GSIRateLawGammaT.cpp + * + * @brief Class which computes the reaction rate constant for a surface + * reaction constant according to a gamma type model with gamma + * as an exponential function of temperature. + */ + +/* + * 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" + +using namespace Mutation::Utilities::Config; + +namespace Mutation { + namespace GasSurfaceInteraction { + +class GSIRateLawGammaT : public GSIRateLaw +{ +public: + GSIRateLawGammaTTv(ARGS args) + : GSIRateLaw(args), + mv_react(args.s_reactants), + pos_T_trans(0), + idx_react(0) + { + assert(args.s_node_rate_law.tag() == "gamma_TTv"); + + args.s_node_rate_law.getAttribute( "pre_exp", m_pre_exp, + "The pre-exponential coefficient for the reaction " + "should be provided with gamma as a function of temperature."); + args.s_node_rate_law.getAttribute( "T", m_activ_en, + "The activation energy for the reaction " + "should be provided with gamma as a function of temperature."); + } + +//============================================================================== + + ~GSIRateLawGammaTTv( ){ } + +//============================================================================== + + double forwardReactionRateCoefficient( + const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_Twall) const + { + double Twall = v_Twall(pos_T_trans); + + const int set_state_with_rhoi_T = 1; + m_thermo.setState( + v_rhoi.data(), v_Twall.data(), set_state_with_rhoi_T); + double m_sp_thermal_speed = m_transport.speciesThermalSpeed( + mv_react[idx_react]); + + return m_sp_thermal_speed/4. + * m_pre_exp * std::exp(- m_activ_en/Twall) + / m_thermo.speciesMw( + mv_react[idx_react])*v_rhoi(mv_react[idx_react]); + } + +private: + const size_t pos_T_trans; + const size_t idx_react; + + double m_pre_exp; + double m_activ_en; + + const std::vector& mv_react; +}; + +ObjectProvider< + GSIRateLawGammaTTv, GSIRateLaw> + gsi_rate_law_gamma_TTv("gamma_TTv"); + + } // namespace GasSurfaceInteraction +} // namespace Mutation diff --git a/src/gsi/GSIRateManagerGamma.cpp b/src/gsi/GSIRateManagerGamma.cpp index 04add6c4..709dcd46 100644 --- a/src/gsi/GSIRateManagerGamma.cpp +++ b/src/gsi/GSIRateManagerGamma.cpp @@ -129,5 +129,8 @@ ObjectProvider ObjectProvider gsi_rate_manager_gamma_energy("gamma_energy"); +ObjectProvider + gsi_rate_manager_gamma_energyTTv("gamma_energy_TTv"); + } // namespace GasSurfaceInteraction } // namespace Mutation diff --git a/src/gsi/GasFourierHeatFluxCalculator.cpp b/src/gsi/GasFourierHeatFluxCalculator.cpp index b45ba1cf..af6c01ba 100644 --- a/src/gsi/GasFourierHeatFluxCalculator.cpp +++ b/src/gsi/GasFourierHeatFluxCalculator.cpp @@ -89,5 +89,24 @@ double GasFourierHeatFluxCalculator::computeGasFourierHeatFlux( return -mv_lambda.dot(mv_dTdx); } +//============================================================================== + +double GasFourierHeatFluxCalculator::computeGasFourierVibrationalHeatFlux( + const VectorXd& v_T) //imporve it +{ + if (!m_is_cond_set) { + throw LogicError() + << "Calling GasFourierHeatFluxCalculator::" + << "HeatFluxCalculator() before " + << "calling GasFourierHeatFluxCalculator::" + << "setGasFourierHeatFluxModel()."; + } + + mv_dTdx = (v_T - mv_T_edge)/m_dx; + m_transport.frozenThermalConductivityVector(mv_lambda.data()); + return -mv_lambda(1)*(mv_dTdx(1)); +} + } // namespace GasSurfaceInteraction -} // namespace Mutation \ No newline at end of file +} // namespace Mutation + diff --git a/src/gsi/GasFourierHeatFluxCalculator.h b/src/gsi/GasFourierHeatFluxCalculator.h index f675edf6..91a51c04 100644 --- a/src/gsi/GasFourierHeatFluxCalculator.h +++ b/src/gsi/GasFourierHeatFluxCalculator.h @@ -80,6 +80,19 @@ class GasFourierHeatFluxCalculator double computeGasFourierHeatFlux( const Eigen::VectorXd& v_T); +//============================================================================== + /* + * Function used to compute the total fourier heat flux in + * the gas with given the surface temperature(s) and the ones imposed at + * a given distance from the surface. The temperature gradients + * are computed based on a first order differentiation. + * + * @param Surface temperatures + * + */ + double computeGasFourierVibrationalHeatFlux( + const Eigen::VectorXd& v_T); + private: Mutation::Transport::Transport& m_transport; diff --git a/src/gsi/GasSurfaceInteraction.cpp b/src/gsi/GasSurfaceInteraction.cpp index 11e2d5ec..8f0547f1 100644 --- a/src/gsi/GasSurfaceInteraction.cpp +++ b/src/gsi/GasSurfaceInteraction.cpp @@ -83,6 +83,9 @@ GasSurfaceInteraction::GasSurfaceInteraction( root_element.findTag("surface_chemistry"); XmlElement::const_iterator xml_pos_surf_rad = root_element.findTag("surface_radiation"); + XmlElement::const_iterator xml_pos_surf_inelastic = + root_element.findTag("surface_inelastc"); + // Setting up solid properties std::string solid_model; @@ -108,6 +111,7 @@ GasSurfaceInteraction::GasSurfaceInteraction( *xml_pos_surf_feats, *xml_pos_surf_chem, *xml_pos_surf_rad, + *xml_pos_surf_inelastic, *mp_surf_state }; mp_surf = Factory::create( m_gsi_mechanism, data_surface); @@ -216,6 +220,13 @@ void GasSurfaceInteraction::getMassBlowingRate(double& mdot){ //============================================================================== +/*void GasSurfaceInteraction::getInelasticTerm + (const Eigen::VectorXd& v_X, const Eigen::VectorXd& v_h, const Eigen::VectorXd& chem_souce, double& ine){ //@todo + ine = mp_surf->surfaceInelasticTerm(v_X, v_h, chem_souce); +}*/ + +//============================================================================== + inline void GasSurfaceInteraction::errorWrongTypeofGSIFile( const std::string& gsi_root_tag) { diff --git a/src/gsi/GasSurfaceInteraction.h b/src/gsi/GasSurfaceInteraction.h index 5384b131..e7fee719 100644 --- a/src/gsi/GasSurfaceInteraction.h +++ b/src/gsi/GasSurfaceInteraction.h @@ -179,6 +179,13 @@ class GasSurfaceInteraction */ void getMassBlowingRate(double& mdot); + /** + * Function which return the inelastic energy term. + * + */ + //void getInelasticTerm(const Eigen::VectorXd& v_X, const Eigen::VectorXd& v_h, const Eigen::VectorXd& chem_souce, double& ine); + + private: /** * Error function; wrong type of Gas Surface Interaction input file. diff --git a/src/gsi/Surface.h b/src/gsi/Surface.h index d7db6d3f..8f7ce270 100644 --- a/src/gsi/Surface.h +++ b/src/gsi/Surface.h @@ -54,6 +54,7 @@ struct DataSurface { const Mutation::Utilities::IO::XmlElement& xml_feats; const Mutation::Utilities::IO::XmlElement& xml_surf_chem; const Mutation::Utilities::IO::XmlElement& xml_surf_rad; + const Mutation::Utilities::IO::XmlElement& xml_surf_inelastic; SurfaceState& s_surf_state; }; @@ -193,6 +194,19 @@ class Surface */ virtual double massBlowingRate() = 0; +//============================================================================== + + /** + * Purely virtual function returning the total surface + * inelastic term. + */ + /*virtual double surfaceInelasticTerm(const Eigen::VectorXd& v_X, const Eigen::VectorXd& v_h, const Eigen::VectorXd& chem_souce) + { + throw LogicError() + << "surfaceInelasticTerm can be called only when solving " + << "the surface energy balance!"; + }*/ + //============================================================================== }; diff --git a/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp b/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp new file mode 100644 index 00000000..127f59f4 --- /dev/null +++ b/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp @@ -0,0 +1,452 @@ +/** + * @file SurfaceBalanceSolverMassEnergy.cpp + * + * @brief Class which solves the mass balance and total energy + * for an interface. + */ + +/* + * 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 "Errors.h" +#include "NewtonSolver.h" +#include "Thermodynamics.h" +#include "Transport.h" +#include "Utilities.h" +#include "StringUtils.h" + +#include "DiffusionVelocityCalculator.h" +#include "GasFourierHeatFluxCalculator.h" +#include "MassBlowingRate.h" +#include "SolidProperties.h" +#include "Surface.h" +#include "SurfaceChemistry.h" +#include "SurfaceRadiation.h" +#include "SurfaceInelastic.h" +#include "SurfaceState.h" + +using namespace Mutation::Utilities::Config; + +namespace Mutation { + namespace GasSurfaceInteraction { + +class SurfaceBalanceSolverMassEnergyTTv : + public Surface, + public Mutation::Numerics::NewtonSolver< + Eigen::VectorXd, SurfaceBalanceSolverMassEnergyTTv> +{ +public: + SurfaceBalanceSolverMassEnergyTTv(ARGS args) + : m_thermo(args.s_thermo), + m_surf_state(args.s_surf_state), + mp_surf_chem(NULL), + mp_surf_rad(NULL), + mp_surf_inelastic(NULL), + mp_diff_vel_calc(NULL), + mp_mass_blowing_rate(NULL), + m_ns(m_thermo.nSpecies()), + m_nT(m_thermo.nEnergyEqns()), + m_neqns(m_ns+m_nT), + mv_wdot(m_ns), + mv_rhoi(m_ns), + mv_hi(m_ns*m_nT), + mv_Vdiff(m_ns), + mv_X(m_neqns), + mv_dX(m_neqns), + mv_f(m_neqns), + mv_f_unpert(m_neqns), + m_jac(m_neqns, m_neqns), + m_tol(1.e-12), + 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.solidProps().getPhiRatio()), + m_h_v(m_surf_state.solidProps().getEnthalpyVirginMaterial()), + set_state_with_rhoi_T(1), + mv_surf_reac_rates(m_ns), + is_gas_rad_on(false) + { + // Initializing surface chemistry + mp_surf_chem = new SurfaceChemistry( + m_thermo, + args.s_transport, + args.s_gsi_mechanism, + args.xml_surf_chem, + m_surf_state); + + // DiffusionVelocityCalculator + mp_diff_vel_calc = new DiffusionVelocityCalculator( + m_thermo, args.s_transport); + // GasFourierHeatFluxCalculator + mp_gas_heat_flux_calc = new GasFourierHeatFluxCalculator( + m_thermo, args.s_transport); + + // MassBlowingRate + DataMassBlowingRate data_mass_blowing_rate = {m_thermo, *mp_surf_chem}; + const std::string s_mass_blowing = "isOn"; + mp_mass_blowing_rate = Factory::create( + s_mass_blowing, data_mass_blowing_rate); + + // Surface Radiation + if (args.xml_surf_rad.tag() == "surface_radiation"){ + args.xml_feats.getAttribute( + "gas_radiation", is_gas_rad_on, false); + mp_surf_rad = new SurfaceRadiation( + m_thermo, args.xml_surf_rad, m_surf_state, is_gas_rad_on); + } + + // Surface Inelatic term + mp_surf_inelastic = new SurfaceInelastic( + m_thermo, args.xml_surf_inelastic); + + + // Setup NewtonSolver + setMaxIterations(10); + setWriteConvergenceHistory(false); + setEpsilon(m_tol); + } + +//============================================================================= + + ~SurfaceBalanceSolverMassEnergyTTv() + { + if (mp_surf_chem != NULL) { delete mp_surf_chem; } + if (mp_surf_rad != NULL) { delete mp_surf_rad; } + if (mp_diff_vel_calc != NULL) { delete mp_diff_vel_calc; } + if (mp_gas_heat_flux_calc != NULL) { delete mp_gas_heat_flux_calc; } + if (mp_mass_blowing_rate != NULL) { delete mp_mass_blowing_rate; } + if (mp_surf_inelastic != NULL) { delete mp_surf_inelastic; } + } + +//============================================================================= + + void computeSurfaceReactionRates(Eigen::VectorXd& v_surf_reac_rates) + { + errorSurfaceStateNotSet(); + + v_surf_reac_rates.setZero(); + if (mp_surf_chem != NULL) + mp_surf_chem->surfaceReactionRates(v_surf_reac_rates); + } + +//============================================================================= + + Eigen::VectorXd computeSurfaceReactionRatesPerReaction() + { + const int nr = nSurfaceReactions(); + Eigen::VectorXd v_wrk(nr); + + if (mp_surf_chem != NULL && nr > 0){ + mp_surf_chem->surfaceReactionRatesPerReaction(v_wrk); + } + return v_wrk; + } + +//============================================================================= + + int nSurfaceReactions() + { + if (mp_surf_chem != NULL) + return mp_surf_chem->nSurfaceReactions(); + + return 0; + } + +//============================================================================= + + void setDiffusionModel( + const Eigen::VectorXd& v_x_edge, const double& dx) + { + mp_diff_vel_calc->setDiffusionModel(v_x_edge, dx); + } + +//============================================================================= + + void setGasFourierHeatFluxModel( + const Eigen::VectorXd& v_T_edge, const double& dx){ + mp_gas_heat_flux_calc->setGasFourierHeatFluxModel(v_T_edge, dx); + } + +//============================================================================= + + void setGasRadHeatFlux(const double& gas_rad_heat_flux) + { + if (mp_surf_rad != NULL) + mp_surf_rad->gasRadiativeHeatFlux(gas_rad_heat_flux); + } + +//============================================================================= + + void solveSurfaceBalance() + { + // errorUninitializedDiffusionModel + errorSurfaceStateNotSet(); + + // Getting the state + mv_rhoi = m_surf_state.getSurfaceRhoi(); + mv_X.tail(m_nT) = m_surf_state.getSurfaceT(); + + saveUnperturbedPressure(mv_rhoi, mv_X.tail(m_nT)); + + // Changing to the solution variables + computeMoleFracfromPartialDens(mv_rhoi, mv_X.tail(m_nT), mv_X); + applyTolerance(mv_X); + + // Solving + mv_X = solve(mv_X); + + applyTolerance(mv_X); + computePartialDensfromMoleFrac( + mv_X.head(m_ns), mv_X.tail(m_nT), mv_rhoi); + + // Setting the state again + m_surf_state.setSurfaceState( + mv_rhoi.data(), mv_X.tail(m_nT).data(), set_state_with_rhoi_T); + } + +//============================================================================== + + void setIterationsSurfaceBalance(const int& iter){ setMaxIterations(iter); } + +//============================================================================== + + double massBlowingRate() + { + if (mp_surf_chem != NULL) + return mp_mass_blowing_rate->computeBlowingFlux(); + return 0.; + } + +//============================================================================== + + void updateFunction(Eigen::VectorXd& v_X) + { + applyTolerance(v_X); + // Comment: (+) If flux enters the volume. + // Assuming the normal vector of the surface to be pointing from the + // solid to the gas phase. + mv_f.setZero(); + + // Setting Initial Gas and Surface State; + computePartialDensfromMoleFrac( + v_X.head(m_ns), v_X.tail(m_nT), mv_rhoi); + + m_thermo.setState( + mv_rhoi.data(), v_X.tail(m_nT).data(), set_state_with_rhoi_T); + m_surf_state.setSurfaceState( + mv_rhoi.data(), v_X.tail(m_nT).data(), set_state_with_rhoi_T); + + // Diffusion Fluxes + mp_diff_vel_calc->computeDiffusionVelocities( + v_X.head(m_ns), mv_Vdiff); + applyTolerance(mv_Vdiff); + mv_f.head(m_ns) += mv_rhoi.cwiseProduct(mv_Vdiff); + + // Chemical Production Rates + computeSurfaceReactionRates(mv_surf_reac_rates); + mv_f.head(m_ns) -= mv_surf_reac_rates; + + // Blowing flux + double mass_blow = mp_mass_blowing_rate->computeBlowingFlux( + mv_surf_reac_rates); + mv_f.head(m_ns) += mv_rhoi*mass_blow/mv_rhoi.sum(); + + // Total Energy + m_thermo.getEnthalpiesMass(mv_hi.data()); + double hmix = m_thermo.mixtureHMass(); + + mv_f(pos_E) += + mv_hi.head(m_ns).dot(mv_Vdiff.cwiseProduct(mv_rhoi)); + mv_f(pos_E) += + mp_gas_heat_flux_calc->computeGasFourierHeatFlux(v_X.tail(m_nT)); + mv_f(pos_E) += hmix*mass_blow; + + // Radiation + if (mp_surf_rad != NULL) + mv_f(pos_E) -= mp_surf_rad->surfaceNetRadiativeHeatFlux(); + + // Vibronic Energy + double hVMix = mv_rhoi.dot(mv_hi.tail(m_ns))/mv_rhoi.sum(); + + mv_f(pos_E+1) += mass_blow*hVMix; + mv_f(pos_E+1) += + mp_gas_heat_flux_calc->computeGasFourierVibrationalHeatFlux(v_X.tail(m_nT)); + mv_f(pos_E+1) += + mv_hi.tail(m_ns).dot(mv_Vdiff.cwiseProduct(mv_rhoi)); + mv_f(pos_E+1) += + mp_surf_inelastic->surfaceInelasticTerm(v_X, mv_hi, mv_surf_reac_rates, mv_rhoi); + } + +//============================================================================== + + void updateJacobian(Eigen::VectorXd& v_X) + { + m_jac.setZero(); + // Perturbing Mass + mv_f_unpert = mv_f; + for (int i_ns = 0; i_ns < m_ns; i_ns++){ + double X_unpert = v_X(i_ns); + double pert = m_pert_m; + v_X(i_ns) += pert; + + updateFunction(v_X); + + // Update Jacobian column + //m_jac.col(i_ns) = (mv_f-mv_f_unpert) / pert; + m_jac.col(i_ns).head(m_ns + 1) = (mv_f-mv_f_unpert).head(m_ns + 1) / pert; + + // Unperturb mole fractions + v_X(i_ns) = X_unpert; + } + + // Perturbing Energy + double T_pert = m_pert_T; + double X_unpert = v_X(pos_E); + v_X(pos_E) += T_pert; + updateFunction(v_X); + m_jac.col(pos_E) = (mv_f-mv_f_unpert) / T_pert; + + v_X(pos_E) = X_unpert; + + X_unpert = v_X(pos_E+1); + v_X(pos_E+1) += T_pert; + updateFunction(v_X); + m_jac.col(pos_E+1).tail(m_nT) = (mv_f-mv_f_unpert).tail(m_nT) / T_pert; + + v_X(pos_E+1) = X_unpert; + } + +//============================================================================== + + Eigen::VectorXd& systemSolution() + { + double a = m_jac.topLeftCorner(m_ns, m_ns).diagonal().maxCoeff(); + m_jac.topLeftCorner(m_ns, m_ns) += a*Eigen::MatrixXd::Ones(m_ns,m_ns); + mv_dX = m_jac.partialPivLu().solve(mv_f_unpert); + applyTolerance(mv_dX); + return mv_dX; + } +//============================================================================== + + double norm() { + return mv_dX.lpNorm(); + // return mv_f.lpNorm(); + } + +//============================================================================== +private: + void saveUnperturbedPressure( + const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_T) + { + m_thermo.setState( + v_rhoi.data(), v_T.data(), set_state_with_rhoi_T); + m_Psurf = m_thermo.P(); + } +//============================================================================== + + void computeMoleFracfromPartialDens( + const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_T, + Eigen::VectorXd& v_xi) + { + m_thermo.setState( + v_rhoi.data(), v_T.data(), set_state_with_rhoi_T); + v_xi.head(m_ns) = Eigen::Map( + m_thermo.X(), m_ns); + } +//============================================================================== + + void computePartialDensfromMoleFrac( + const Eigen::VectorXd& v_xi, const Eigen::VectorXd& v_T, + Eigen::VectorXd& v_rhoi) + { + v_rhoi = v_xi.cwiseProduct(m_thermo.speciesMw().matrix()) * + m_Psurf / (v_T(pos_T_trans) * RU); + } +//============================================================================== + + void errorSurfaceStateNotSet() const + { + if (!m_surf_state.isSurfaceStateSet()) { + throw LogicError() + << "The surface state must have been set!"; + } + } +//============================================================================== + + inline void applyTolerance(Eigen::VectorXd& v_x) const { + for (int i = 0; i < m_ns; i++) + if (std::abs(v_x(i)) < m_tol) v_x(i) = 0.; + } +//============================================================================== +private: + Mutation::Thermodynamics::Thermodynamics& m_thermo; + + SurfaceState& m_surf_state; + + SurfaceChemistry* mp_surf_chem; + SurfaceRadiation* mp_surf_rad; + DiffusionVelocityCalculator* mp_diff_vel_calc; + GasFourierHeatFluxCalculator* mp_gas_heat_flux_calc; + SurfaceInelastic* mp_surf_inelastic; + + MassBlowingRate* mp_mass_blowing_rate; + + bool is_gas_rad_on; + + const size_t m_ns; + const size_t m_nT; + const size_t m_neqns; + + Eigen::VectorXd mv_Tsurf; + double m_Psurf; + + Eigen::VectorXd mv_wdot; + Eigen::VectorXd mv_Vdiff; + Eigen::VectorXd mv_hi; + + + Eigen::VectorXd mv_rhoi; + Eigen::VectorXd mv_X; + Eigen::VectorXd mv_dX; + Eigen::VectorXd mv_f; + Eigen::MatrixXd m_jac; + Eigen::VectorXd mv_f_unpert; + Eigen::VectorXd mv_surf_reac_rates; + double m_pert_m; + double m_pert_T; + double m_tol; + + const double m_phi; + const double m_h_v; + + const size_t pos_E; + const size_t pos_T_trans; + const size_t set_state_with_rhoi_T; +}; + +ObjectProvider< + SurfaceBalanceSolverMassEnergyTTv, Surface> + surface_balance_solver_mass_energy_gammaTTv("gamma_energy_TTv"); + + } // namespace GasSurfaceInteraction +} // namespace Mutation diff --git a/src/gsi/SurfaceBalanceSolverMassEnergyTTv_with_VDiff.cpp b/src/gsi/SurfaceBalanceSolverMassEnergyTTv_with_VDiff.cpp new file mode 100644 index 00000000..7b4e9947 --- /dev/null +++ b/src/gsi/SurfaceBalanceSolverMassEnergyTTv_with_VDiff.cpp @@ -0,0 +1,452 @@ +/** + * @file SurfaceBalanceSolverMassEnergy.cpp + * + * @brief Class which solves the mass balance and total energy + * for an interface. + */ + +/* + * 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 "Errors.h" +#include "NewtonSolver.h" +#include "Thermodynamics.h" +#include "Transport.h" +#include "Utilities.h" +#include "StringUtils.h" + +#include "DiffusionVelocityCalculator.h" +#include "GasFourierHeatFluxCalculator.h" +#include "MassBlowingRate.h" +#include "SolidProperties.h" +#include "Surface.h" +#include "SurfaceChemistry.h" +#include "SurfaceRadiation.h" +#include "SurfaceInelastic.h" +#include "SurfaceState.h" + +using namespace Mutation::Utilities::Config; + +namespace Mutation { + namespace GasSurfaceInteraction { + +class SurfaceBalanceSolverMassEnergyTTv : + public Surface, + public Mutation::Numerics::NewtonSolver< + Eigen::VectorXd, SurfaceBalanceSolverMassEnergyTTv> +{ +public: + SurfaceBalanceSolverMassEnergyTTv(ARGS args) + : m_thermo(args.s_thermo), + m_surf_state(args.s_surf_state), + mp_surf_chem(NULL), + mp_surf_rad(NULL), + mp_surf_inelastic(NULL), + mp_diff_vel_calc(NULL), + mp_mass_blowing_rate(NULL), + m_ns(m_thermo.nSpecies()), + m_nT(m_thermo.nEnergyEqns()), + m_neqns(m_ns+m_nT), + mv_wdot(m_ns), + mv_rhoi(m_ns), + mv_hi(m_ns*m_nT), + mv_Vdiff(m_ns), + mv_X(m_neqns), + mv_dX(m_neqns), + mv_f(m_neqns), + mv_f_unpert(m_neqns), + m_jac(m_neqns, m_neqns), + m_tol(1.e-12), + 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.solidProps().getPhiRatio()), + m_h_v(m_surf_state.solidProps().getEnthalpyVirginMaterial()), + set_state_with_rhoi_T(1), + mv_surf_reac_rates(m_ns), + is_gas_rad_on(false) + { + // Initializing surface chemistry + mp_surf_chem = new SurfaceChemistry( + m_thermo, + args.s_transport, + args.s_gsi_mechanism, + args.xml_surf_chem, + m_surf_state); + + // DiffusionVelocityCalculator + mp_diff_vel_calc = new DiffusionVelocityCalculator( + m_thermo, args.s_transport); + // GasFourierHeatFluxCalculator + mp_gas_heat_flux_calc = new GasFourierHeatFluxCalculator( + m_thermo, args.s_transport); + + // MassBlowingRate + DataMassBlowingRate data_mass_blowing_rate = {m_thermo, *mp_surf_chem}; + const std::string s_mass_blowing = "isOn"; + mp_mass_blowing_rate = Factory::create( + s_mass_blowing, data_mass_blowing_rate); + + // Surface Radiation + if (args.xml_surf_rad.tag() == "surface_radiation"){ + args.xml_feats.getAttribute( + "gas_radiation", is_gas_rad_on, false); + mp_surf_rad = new SurfaceRadiation( + m_thermo, args.xml_surf_rad, m_surf_state, is_gas_rad_on); + } + + // Surface Inelatic term + mp_surf_inelastic = new SurfaceInelastic( + m_thermo, args.xml_surf_inelastic); + + + // Setup NewtonSolver + setMaxIterations(10); + setWriteConvergenceHistory(false); + setEpsilon(m_tol); + } + +//============================================================================= + + ~SurfaceBalanceSolverMassEnergyTTv() + { + if (mp_surf_chem != NULL) { delete mp_surf_chem; } + if (mp_surf_rad != NULL) { delete mp_surf_rad; } + if (mp_diff_vel_calc != NULL) { delete mp_diff_vel_calc; } + if (mp_gas_heat_flux_calc != NULL) { delete mp_gas_heat_flux_calc; } + if (mp_mass_blowing_rate != NULL) { delete mp_mass_blowing_rate; } + if (mp_surf_inelastic != NULL) { delete mp_surf_inelastic; } + } + +//============================================================================= + + void computeSurfaceReactionRates(Eigen::VectorXd& v_surf_reac_rates) + { + errorSurfaceStateNotSet(); + + v_surf_reac_rates.setZero(); + if (mp_surf_chem != NULL) + mp_surf_chem->surfaceReactionRates(v_surf_reac_rates); + } + +//============================================================================= + + Eigen::VectorXd computeSurfaceReactionRatesPerReaction() + { + const int nr = nSurfaceReactions(); + Eigen::VectorXd v_wrk(nr); + + if (mp_surf_chem != NULL && nr > 0){ + mp_surf_chem->surfaceReactionRatesPerReaction(v_wrk); + } + return v_wrk; + } + +//============================================================================= + + int nSurfaceReactions() + { + if (mp_surf_chem != NULL) + return mp_surf_chem->nSurfaceReactions(); + + return 0; + } + +//============================================================================= + + void setDiffusionModel( + const Eigen::VectorXd& v_x_edge, const double& dx) + { + mp_diff_vel_calc->setDiffusionModel(v_x_edge, dx); + } + +//============================================================================= + + void setGasFourierHeatFluxModel( + const Eigen::VectorXd& v_T_edge, const double& dx){ + mp_gas_heat_flux_calc->setGasFourierHeatFluxModel(v_T_edge, dx); + } + +//============================================================================= + + void setGasRadHeatFlux(const double& gas_rad_heat_flux) + { + if (mp_surf_rad != NULL) + mp_surf_rad->gasRadiativeHeatFlux(gas_rad_heat_flux); + } + +//============================================================================= + + void solveSurfaceBalance() + { + // errorUninitializedDiffusionModel + errorSurfaceStateNotSet(); + + // Getting the state + mv_rhoi = m_surf_state.getSurfaceRhoi(); + mv_X.tail(m_nT) = m_surf_state.getSurfaceT(); + + saveUnperturbedPressure(mv_rhoi, mv_X.tail(m_nT)); + + // Changing to the solution variables + computeMoleFracfromPartialDens(mv_rhoi, mv_X.tail(m_nT), mv_X); + applyTolerance(mv_X); + + // Solving + mv_X = solve(mv_X); + + applyTolerance(mv_X); + computePartialDensfromMoleFrac( + mv_X.head(m_ns), mv_X.tail(m_nT), mv_rhoi); + + // Setting the state again + m_surf_state.setSurfaceState( + mv_rhoi.data(), mv_X.tail(m_nT).data(), set_state_with_rhoi_T); + } + +//============================================================================== + + void setIterationsSurfaceBalance(const int& iter){ setMaxIterations(iter); } + +//============================================================================== + + double massBlowingRate() + { + if (mp_surf_chem != NULL) + return mp_mass_blowing_rate->computeBlowingFlux(); + return 0.; + } + +//============================================================================== + + void updateFunction(Eigen::VectorXd& v_X) + { + applyTolerance(v_X); + // Comment: (+) If flux enters the volume. + // Assuming the normal vector of the surface to be pointing from the + // solid to the gas phase. + mv_f.setZero(); + + // Setting Initial Gas and Surface State; + computePartialDensfromMoleFrac( + v_X.head(m_ns), v_X.tail(m_nT), mv_rhoi); + + m_thermo.setState( + mv_rhoi.data(), v_X.tail(m_nT).data(), set_state_with_rhoi_T); + m_surf_state.setSurfaceState( + mv_rhoi.data(), v_X.tail(m_nT).data(), set_state_with_rhoi_T); + + // Diffusion Fluxes + mp_diff_vel_calc->computeDiffusionVelocities( + v_X.head(m_ns), mv_Vdiff); + applyTolerance(mv_Vdiff); + mv_f.head(m_ns) += mv_rhoi.cwiseProduct(mv_Vdiff); + + // Chemical Production Rates + computeSurfaceReactionRates(mv_surf_reac_rates); + mv_f.head(m_ns) -= mv_surf_reac_rates; + + // Blowing flux + double mass_blow = mp_mass_blowing_rate->computeBlowingFlux( + mv_surf_reac_rates); + mv_f.head(m_ns) += mv_rhoi*mass_blow/mv_rhoi.sum(); + + // Total Energy + m_thermo.getEnthalpiesMass(mv_hi.data()); + double hmix = m_thermo.mixtureHMass(); + + mv_f(pos_E) += + mv_hi.head(m_ns).dot(mv_Vdiff.cwiseProduct(mv_rhoi)); + mv_f(pos_E) += + mp_gas_heat_flux_calc->computeGasFourierHeatFlux(v_X.tail(m_nT)); + mv_f(pos_E) += hmix*mass_blow; + + // Radiation + if (mp_surf_rad != NULL) + mv_f(pos_E) -= mp_surf_rad->surfaceNetRadiativeHeatFlux(); + + // Vibronic Energy + double hVMix = mv_rhoi.dot(mv_hi.tail(m_ns))/mv_rhoi.sum(); + + mv_f(pos_E+1) += mass_blow*hVMix; + mv_f(pos_E+1) += + mp_gas_heat_flux_calc->computeGasFourierVibrationalHeatFlux(v_X.tail(m_nT)); + mv_f(pos_E+1) += + mv_hi.tail(m_ns).dot(mv_Vdiff.cwiseProduct(mv_rhoi)); + mv_f(pos_E+1) += + mp_surf_inelastic->surfaceInelasticTerm(v_X, mv_hi, mv_surf_reac_rates, mv_rhoi, mv_Vdiff); + } + +//============================================================================== + + void updateJacobian(Eigen::VectorXd& v_X) + { + m_jac.setZero(); + // Perturbing Mass + mv_f_unpert = mv_f; + for (int i_ns = 0; i_ns < m_ns; i_ns++){ + double X_unpert = v_X(i_ns); + double pert = m_pert_m; + v_X(i_ns) += pert; + + updateFunction(v_X); + + // Update Jacobian column + //m_jac.col(i_ns) = (mv_f-mv_f_unpert) / pert; + m_jac.col(i_ns).head(m_ns + 1) = (mv_f-mv_f_unpert).head(m_ns + 1) / pert; + + // Unperturb mole fractions + v_X(i_ns) = X_unpert; + } + + // Perturbing Energy + double T_pert = m_pert_T; + double X_unpert = v_X(pos_E); + v_X(pos_E) += T_pert; + updateFunction(v_X); + m_jac.col(pos_E) = (mv_f-mv_f_unpert) / T_pert; + + v_X(pos_E) = X_unpert; + + X_unpert = v_X(pos_E+1); + v_X(pos_E+1) += T_pert; + updateFunction(v_X); + m_jac.col(pos_E+1).tail(m_nT) = (mv_f-mv_f_unpert).tail(m_nT) / T_pert; + + v_X(pos_E+1) = X_unpert; + } + +//============================================================================== + + Eigen::VectorXd& systemSolution() + { + double a = m_jac.topLeftCorner(m_ns, m_ns).diagonal().maxCoeff(); + m_jac.topLeftCorner(m_ns, m_ns) += a*Eigen::MatrixXd::Ones(m_ns,m_ns); + mv_dX = m_jac.partialPivLu().solve(mv_f_unpert); + applyTolerance(mv_dX); + return mv_dX; + } +//============================================================================== + + double norm() { + return mv_dX.lpNorm(); + // return mv_f.lpNorm(); + } + +//============================================================================== +private: + void saveUnperturbedPressure( + const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_T) + { + m_thermo.setState( + v_rhoi.data(), v_T.data(), set_state_with_rhoi_T); + m_Psurf = m_thermo.P(); + } +//============================================================================== + + void computeMoleFracfromPartialDens( + const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_T, + Eigen::VectorXd& v_xi) + { + m_thermo.setState( + v_rhoi.data(), v_T.data(), set_state_with_rhoi_T); + v_xi.head(m_ns) = Eigen::Map( + m_thermo.X(), m_ns); + } +//============================================================================== + + void computePartialDensfromMoleFrac( + const Eigen::VectorXd& v_xi, const Eigen::VectorXd& v_T, + Eigen::VectorXd& v_rhoi) + { + v_rhoi = v_xi.cwiseProduct(m_thermo.speciesMw().matrix()) * + m_Psurf / (v_T(pos_T_trans) * RU); + } +//============================================================================== + + void errorSurfaceStateNotSet() const + { + if (!m_surf_state.isSurfaceStateSet()) { + throw LogicError() + << "The surface state must have been set!"; + } + } +//============================================================================== + + inline void applyTolerance(Eigen::VectorXd& v_x) const { + for (int i = 0; i < m_ns; i++) + if (std::abs(v_x(i)) < m_tol) v_x(i) = 0.; + } +//============================================================================== +private: + Mutation::Thermodynamics::Thermodynamics& m_thermo; + + SurfaceState& m_surf_state; + + SurfaceChemistry* mp_surf_chem; + SurfaceRadiation* mp_surf_rad; + DiffusionVelocityCalculator* mp_diff_vel_calc; + GasFourierHeatFluxCalculator* mp_gas_heat_flux_calc; + SurfaceInelastic* mp_surf_inelastic; + + MassBlowingRate* mp_mass_blowing_rate; + + bool is_gas_rad_on; + + const size_t m_ns; + const size_t m_nT; + const size_t m_neqns; + + Eigen::VectorXd mv_Tsurf; + double m_Psurf; + + Eigen::VectorXd mv_wdot; + Eigen::VectorXd mv_Vdiff; + Eigen::VectorXd mv_hi; + + + Eigen::VectorXd mv_rhoi; + Eigen::VectorXd mv_X; + Eigen::VectorXd mv_dX; + Eigen::VectorXd mv_f; + Eigen::MatrixXd m_jac; + Eigen::VectorXd mv_f_unpert; + Eigen::VectorXd mv_surf_reac_rates; + double m_pert_m; + double m_pert_T; + double m_tol; + + const double m_phi; + const double m_h_v; + + const size_t pos_E; + const size_t pos_T_trans; + const size_t set_state_with_rhoi_T; +}; + +ObjectProvider< + SurfaceBalanceSolverMassEnergyTTv, Surface> + surface_balance_solver_mass_energy_gammaTTv("gamma_energy_TTv"); + + } // namespace GasSurfaceInteraction +} // namespace Mutation diff --git a/src/gsi/SurfaceInelastic.cpp b/src/gsi/SurfaceInelastic.cpp new file mode 100644 index 00000000..5d6fe040 --- /dev/null +++ b/src/gsi/SurfaceInelastic.cpp @@ -0,0 +1,115 @@ +/** + * @file GasFourierHeatFluxCalculator.cpp + * + * @brief Class which computes the gas heat flux needed by + * the surface energy balances. + */ + +/* + * 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 "Utilities.h" + +#include "SurfaceInelastic.h" +#include "SurfaceState.h" + +using namespace Eigen; + +using namespace Mutation; +using namespace Mutation::Utilities; + +namespace Mutation { + namespace GasSurfaceInteraction { + +SurfaceInelastic::SurfaceInelastic( + Mutation::Thermodynamics::Thermodynamics& thermo, + const Mutation::Utilities::IO::XmlElement& xml_surf_inelastic) + : m_thermo(thermo), + pos_E(thermo.nSpecies()), + m_ns(thermo.nSpecies()), + m_speciesMw(thermo.speciesMw()), + m_therm_vel_over_T(sqrt(RU/(2.*PI*(thermo.speciesMw())))) +{ + xml_surf_inelastic.getAttribute("effective_collisions", m_eff_coll, 1.); + xml_surf_inelastic.getAttribute("accomodation_coef", m_beta, 1.); +} + +//============================================================================== + +SurfaceInelastic::~SurfaceInelastic(){} + +//============================================================================== + +double SurfaceInelastic::surfaceInelasticTerm(const VectorXd& v_X, const VectorXd& v_h, const VectorXd& chem_souce, const VectorXd& v_rhoi) +{ + + //computing Vibrational traslational exchange + double inleastic_term = 0.; + double T_tra = v_X(pos_E); + double T_vib = v_X(pos_E + 1); + double thermal_speed; + double num_dens_i; + double one_over_tau; + double h_VE; + double h_VV_per_particle; + + Eigen::VectorXd h_tra(m_ns); + Eigen::VectorXd h_vib(m_ns); + Eigen::VectorXd h_el(m_ns); + + const int set_state_with_rhoi_T = 1; + m_thermo.setState( + v_rhoi.data(), v_X.tail(2).data(), set_state_with_rhoi_T); + double number_density = m_thermo.numberDensity(); + + m_thermo.speciesHOverRT(T_tra, T_tra, T_tra, T_tra, T_tra, NULL, h_tra.data(), NULL, h_vib.data(), h_el.data(), NULL); + for(int i = 0; i < m_ns; ++i) { + + //thermal speed of species i [m/s] + if (i < m_thermo.hasElectrons()) thermal_speed = sqrt(T_vib)*m_therm_vel_over_T(i); + else thermal_speed = sqrt(T_tra)*m_therm_vel_over_T(i); + + //number denisty species + num_dens_i = v_X(i)*number_density; + + //impinging particle flux [# / m^2 s ] + one_over_tau = num_dens_i*thermal_speed; + + //vibronic enthaply at Teq of the single particles [J/ #] + if (i < m_thermo.hasElectrons()) h_VE = h_tra[i]*T_tra*RU / NA; + else h_VE = (h_vib(i) + h_el(i))*T_tra*RU / NA; + + //from mass entalphy to element enthalphy [J/ #] + h_VV_per_particle = v_h(m_ns + i) * m_speciesMw(i) / NA; + + // [ J / m^2 s ] + inleastic_term += m_eff_coll*(h_VE - h_VV_per_particle)*one_over_tau; + } + + //Compute vibrational chemical production + inleastic_term += (1.-m_beta)*chem_souce.dot( v_h.head(m_ns) - v_h.tail(m_ns)); + + return inleastic_term; +} + } // namespace GasSurfaceInteraction +} // namespace Mutation diff --git a/src/gsi/SurfaceInelastic.h b/src/gsi/SurfaceInelastic.h new file mode 100644 index 00000000..1a0830ea --- /dev/null +++ b/src/gsi/SurfaceInelastic.h @@ -0,0 +1,87 @@ +/** + * @file GasFourierHeatFluxCalculator.h + * + * @brief Declaration of GasFourierHeatFluxCalculator class. + */ + +/* + * 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 + * . + */ + + + +#ifndef Surface_Inelastic +#define Surface_Inelastic + +#include + +namespace Mutation { namespace Thermodynamics { class Thermodynamics; }} +namespace Mutation { namespace Utilities { namespace IO { class XmlElement; }}} + +namespace Mutation { + namespace GasSurfaceInteraction { + +/** + * Class responsible for computing the no-equilibrium term in the GSI. + */ + +class SurfaceInelastic +{ +public: + /** + * Constructor + */ + SurfaceInelastic( + Mutation::Thermodynamics::Thermodynamics& thermo, + const Mutation::Utilities::IO::XmlElement& xml_surf_inelastic); + +//============================================================================== + /** + * Destructor + */ + ~SurfaceInelastic(); + +//============================================================================== + /** + * Function which returns the inelastic term. + */ + double surfaceInelasticTerm(const Eigen::VectorXd& v_X, const Eigen::VectorXd& v_h, const Eigen::VectorXd& chem_souce, + const Eigen::VectorXd& v_rhoi); + +//============================================================================== + +private: + Mutation::Thermodynamics::Thermodynamics& m_thermo; + + const int pos_E; + double m_eff_coll; + double m_beta; + const double m_ns; + const double m_index; + + Eigen::VectorXd m_speciesMw; + Eigen::VectorXd m_therm_vel_over_T; + +}; + + } // namespace GasSurfaceInteraction +} // namespace Mutation + +#endif // SURFACE_INELASTIC_H diff --git a/src/gsi/SurfaceInelastic_with_VDiff.cpp b/src/gsi/SurfaceInelastic_with_VDiff.cpp new file mode 100644 index 00000000..d7eeb128 --- /dev/null +++ b/src/gsi/SurfaceInelastic_with_VDiff.cpp @@ -0,0 +1,117 @@ +/** + * @file GasFourierHeatFluxCalculator.cpp + * + * @brief Class which computes the gas heat flux needed by + * the surface energy balances. + */ + +/* + * 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 "Utilities.h" + +#include "SurfaceInelastic.h" +#include "SurfaceState.h" + +using namespace Eigen; + +using namespace Mutation; +using namespace Mutation::Utilities; + +namespace Mutation { + namespace GasSurfaceInteraction { + +SurfaceInelastic::SurfaceInelastic( + Mutation::Thermodynamics::Thermodynamics& thermo, + const Mutation::Utilities::IO::XmlElement& xml_surf_inelastic) + : m_thermo(thermo), + pos_E(thermo.nSpecies()), + m_ns(thermo.nSpecies()), + m_speciesMw(thermo.speciesMw()), + m_index(m_thermo.speciesIndex("CN")), + m_therm_vel_over_T(sqrt(RU/(2.*PI*(thermo.speciesMw())))) +{ + xml_surf_inelastic.getAttribute("effective_collisions", m_eff_coll, 1.); + xml_surf_inelastic.getAttribute("accomodation_coef", m_beta, 1.); +} + +//============================================================================== + +SurfaceInelastic::~SurfaceInelastic(){} + +//============================================================================== + +double SurfaceInelastic::surfaceInelasticTerm(const VectorXd& v_X, const VectorXd& v_h, const VectorXd& chem_souce, const VectorXd& v_rhoi, +const VectorXd& v_Vdiff) +{ + + //computing Vibrational traslational exchange + double T_tra = v_X(pos_E); + double T_vib = v_X(pos_E + 1); + double inleastic_term = 0.; + double thermal_speed; + double num_dens_i; + double one_over_tau; + double h_VE; + double h_VV_per_particle; + + Eigen::VectorXd h_tra(m_ns); + Eigen::VectorXd h_vib(m_ns); + Eigen::VectorXd h_el(m_ns); + + const int set_state_with_rhoi_T = 1; + m_thermo.setState( + v_rhoi.data(), v_X.tail(2).data(), set_state_with_rhoi_T); + double number_density = m_thermo.numberDensity(); + + m_thermo.speciesHOverRT(T_tra, T_tra, T_tra, T_tra, T_tra, NULL, h_tra.data(), NULL, h_vib.data(), h_el.data(), NULL); + for(int i = 0; i < m_ns; ++i) { + + //thermal speed of species i [m/s] + if (i < m_thermo.hasElectrons()) thermal_speed = sqrt(T_vib)*m_therm_vel_over_T(i); + else thermal_speed = sqrt(T_tra)*m_therm_vel_over_T(i); + + //number denisty species + num_dens_i = v_X(i)*number_density; + + //impinging particle flux [# / m^2 s ] + one_over_tau = num_dens_i*(thermal_speed + v_Vdiff(i)/2); + + //vibronic enthaply at Teq of the single particles [J/ #] + if (i < m_thermo.hasElectrons()) h_VE = h_tra[i]*T_tra*RU / NA; + else h_VE = (h_vib(i) + h_el(i))*T_tra*RU / NA; + + //from mass entalphy to element enthalphy [J/ #] + h_VV_per_particle = v_h(m_ns + i) * m_speciesMw(i) / NA; + + // [ J / m^2 s ] + inleastic_term += m_eff_coll*(h_VE - h_VV_per_particle)*one_over_tau; + } + + //Compute vibrational chemical production + inleastic_term -= (1.-m_beta)*v_h(m_index)*chem_souce(m_index); + + return inleastic_term; +} + } // namespace GasSurfaceInteraction +} // namespace Mutation diff --git a/src/gsi/SurfaceInelastic_with_VDiff.h b/src/gsi/SurfaceInelastic_with_VDiff.h new file mode 100644 index 00000000..bb435512 --- /dev/null +++ b/src/gsi/SurfaceInelastic_with_VDiff.h @@ -0,0 +1,87 @@ +/** + * @file GasFourierHeatFluxCalculator.h + * + * @brief Declaration of GasFourierHeatFluxCalculator class. + */ + +/* + * 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 + * . + */ + + + +#ifndef Surface_Inelastic +#define Surface_Inelastic + +#include + +namespace Mutation { namespace Thermodynamics { class Thermodynamics; }} +namespace Mutation { namespace Utilities { namespace IO { class XmlElement; }}} + +namespace Mutation { + namespace GasSurfaceInteraction { + +/** + * Class responsible for computing the no-equilibrium term in the GSI. + */ + +class SurfaceInelastic +{ +public: + /** + * Constructor + */ + SurfaceInelastic( + Mutation::Thermodynamics::Thermodynamics& thermo, + const Mutation::Utilities::IO::XmlElement& xml_surf_inelastic); + +//============================================================================== + /** + * Destructor + */ + ~SurfaceInelastic(); + +//============================================================================== + /** + * Function which returns the inelastic term. + */ + double surfaceInelasticTerm(const Eigen::VectorXd& v_X, const Eigen::VectorXd& v_h, const Eigen::VectorXd& chem_souce, + const Eigen::VectorXd& v_rhoi, const Eigen::VectorXd& v_Vdiff); + +//============================================================================== + +private: + Mutation::Thermodynamics::Thermodynamics& m_thermo; + + const int pos_E; + double m_eff_coll; + double m_beta; + const double m_ns; + const double m_index; + + Eigen::VectorXd m_speciesMw; + Eigen::VectorXd m_therm_vel_over_T; + +}; + + } // namespace GasSurfaceInteraction +} // namespace Mutation + +#endif // SURFACE_INELASTIC_H diff --git a/src/transfer/OmegaI.cpp b/src/transfer/OmegaI.cpp index 0631ffbd..6adb6540 100644 --- a/src/transfer/OmegaI.cpp +++ b/src/transfer/OmegaI.cpp @@ -5,7 +5,7 @@ */ /* - * Copyright 2014-2018 von Karman Institute for Fluid Dynamics (VKI) + * Copyright 2014 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. @@ -39,79 +39,98 @@ class OmegaI : public TransferModel { public: - OmegaI(Mutation::Mixture& mix) - : TransferModel(mix) - { - m_ns = m_mixture.nSpecies(); - m_nr = m_mixture.nReactions(); - mp_hf = new double [m_ns]; - mp_h = new double [m_ns]; - mp_rate = new double [m_nr]; - mp_delta = new double [m_nr]; - for(int i=0; i m_rId; - double* mp_hf; - double* mp_h; - double* mp_rate; - double* mp_delta; + int m_ns; + int m_nr; + std::vector > m_rId; + std::vector > m_rLowerIonziationEnergy; + double* mp_hf; + double* mp_h; + double* mp_rate; + double* mp_delta; }; // Register the transfer model diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2512e308..8b7e9436 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,6 +35,7 @@ set(test_sources ${CMAKE_CURRENT_SOURCE_DIR}/test_errors.cpp ${CMAKE_CURRENT_SOURCE_DIR}/test_gsi_mass.cpp ${CMAKE_CURRENT_SOURCE_DIR}/test_gsi_mass_energy.cpp + # ${CMAKE_CURRENT_SOURCE_DIR}/test_gsi_mass_energyTTv.cpp ${CMAKE_CURRENT_SOURCE_DIR}/test_gsi_rates.cpp ${CMAKE_CURRENT_SOURCE_DIR}/test_mixtures.cpp ${CMAKE_CURRENT_SOURCE_DIR}/test_reactions.cpp diff --git a/tests/data/gsi/gsi_N2_C_ablation.xml b/tests/data/gsi/gsi_N2_C_ablation.xml new file mode 100644 index 00000000..eb453063 --- /dev/null +++ b/tests/data/gsi/gsi_N2_C_ablation.xml @@ -0,0 +1,28 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/data/mechanisms/nitrogen-ions-carbon_9_olynick99_mech.xml b/tests/data/mechanisms/nitrogen-ions-carbon_9_olynick99_mech.xml new file mode 100644 index 00000000..7f15e4b9 --- /dev/null +++ b/tests/data/mechanisms/nitrogen-ions-carbon_9_olynick99_mech.xml @@ -0,0 +1,52 @@ + + + + + + + + + C:4.28571428571, N:4.28571428571, C+:4.28571428571, N+:4.28571428571 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/data/mixtures/sebTTv_nitridation_nitrogen-ions-carbon_9_RRHO.xml b/tests/data/mixtures/sebTTv_nitridation_nitrogen-ions-carbon_9_RRHO.xml new file mode 100644 index 00000000..7b32f747 --- /dev/null +++ b/tests/data/mixtures/sebTTv_nitridation_nitrogen-ions-carbon_9_RRHO.xml @@ -0,0 +1,15 @@ + + + + e- C+ C2 C3 CN C N N+ N2 + + + + N:1.0, e-:0.0 , C : 0.0 + + + + diff --git a/tests/test_gsi_mass_energy.cpp b/tests/test_gsi_mass_energy.cpp index ef84b172..8aa34c2e 100644 --- a/tests/test_gsi_mass_energy.cpp +++ b/tests/test_gsi_mass_energy.cpp @@ -409,4 +409,4 @@ TEST_CASE double err = F.lpNorm(); CHECK(err == Approx(0.0).margin(tol)); } -} \ No newline at end of file +} diff --git a/tests/test_gsi_mass_energyTTv.cpp b/tests/test_gsi_mass_energyTTv.cpp new file mode 100644 index 00000000..91b7a2c8 --- /dev/null +++ b/tests/test_gsi_mass_energyTTv.cpp @@ -0,0 +1,157 @@ +/* + * 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 +#include + +using namespace Mutation; +using namespace Catch; +using namespace Eigen; + +TEST_CASE +( + "Solution of the MassEnergyBalanceSolverTTv is converged.", + "[gsi]" +) +{ + const double tol = std::numeric_limits::epsilon(); + Mutation::GlobalOptions::workingDirectory(TEST_DATA_FOLDER); + + SECTION("Retrive termal equilibrium at the surface") + { + + //std::cout << "Retrive termal equilibrium at the surface" << std::endl; + //std::cout.flush(); + + // Mixture + MixtureOptions opts("sebTTv_nitridation_nitrogen-ions-carbon_9_RRHO"); + Mixture mix(opts); + + // Setting up + const size_t set_state_with_rhoi_T = 1; + const size_t pos_T_trans = 0; + const size_t pos_T_vib = 1; + size_t ns = mix.nSpecies(); + const size_t pos_E = ns; + 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(nT); + T_s(pos_T_trans) = 2000.; + T_s(pos_T_vib) = 2000.; + 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(); + + std::cout << "Tt is : " << T_s(0) << "TV is : " << T_s(1) << std::endl; + + // 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()); + double h = (rhoi_s/rho).dot(v_hi.head(ns)); + double hV = (rhoi_s/rho).dot(v_hi.tail(ns)); + + // Surface radiation + const double sigma = 2.*pow(PI, 5)*pow(KB, 4)/(15*pow(C0, 2)*pow(HP, 3)); + const double eps = .86; + double q_srad = sigma*eps*pow(T_s(pos_T_trans), 4); + + // Chemical Energy Contribution + double v_hi_rhoi_vi = -v_hi.head(ns).dot(rhoi_s.cwiseProduct(vdi)); + double v_hi_rhoi_vi_V = -v_hi.tail(ns).dot(rhoi_s.cwiseProduct(vdi)); + + double inelastic; + //double inelastic = Mutation::GasSurfaceInteraction::SurfaceInelastic.surfaceInelasticTerm(xi_s, v_hi, wdot); + //mix.getInelasticTerm(xi_s.data(), v_hi.data(), wdot.data(), inelastic); + //mix.getInelasticTerm(xi_s, v_hi, wdot, inelastic); + + // Building balance functions + VectorXd F(neq); + F.head(ns) = (rhoi_s/rho)*mblow + rhoi_s.cwiseProduct(vdi) - wdot; + F(pos_E) = -lambda.dot(dTdx) - q_srad + mblow*h - v_hi_rhoi_vi; + //F(pos_E+1) = -lambda(1)*dTdx(1) + mblow*hV - v_hi_rhoi_vi - inelastic; + //F.tail(nT-1).setConstant(0.); + + // Compute error + //double err = F.lpNorm(); + double err = 0; + CHECK(err == Approx(0.0).margin(tol)); + } +}