From 378ae77e3d6af8ed0beeb05ddc8620bab12043dd Mon Sep 17 00:00:00 2001 From: grgbellasvki <44604564+grgbellasvki@users.noreply.github.com> Date: Thu, 28 Mar 2019 10:06:22 +0100 Subject: [PATCH 1/9] Updating Documentation. Work In Progress. --- docs/input-files.md | 555 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 555 insertions(+) diff --git a/docs/input-files.md b/docs/input-files.md index acb659bb..26c77cdb 100644 --- a/docs/input-files.md +++ b/docs/input-files.md @@ -20,6 +20,8 @@ - [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) ## Introduction @@ -33,6 +35,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 @@ -949,3 +952,555 @@ 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, \f$\partial \mathcal / \partial t\f$, +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: + +(**Fg** - **Fb**) . **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: +\f[ +\dot{\omega}_{A,\text{cat}} = \gamma_{A} m_{A} F_{A}^\downarrow \quad \text{and} \quad +\dot{\omega}_{A_2,\text{cat}} = - \gamma_{A} m_{A} F_{A}^\downarrow +\f] +with \f$ m_{A} \f$ being the mass of species \f$ A \f$. + +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$ F_{A}^\downarrow \f$, is equal to +\f[ +F_{A}^\downarrow = n_{A} \sqrt{\frac{k_B T_{\text{s}}}{2 \pi m_{A}}} +\f] +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 +\f[ +A + B \rightarrow AB, +\f] +with the chemical rate for the produced molecule \f$ AB \f$ equal to: +\f[ +\dot{\omega}_{AB,\text{cat}} = - \gamma_{AB} m_{A} F_{A}^\downarrow - \gamma_{BA} m_{B} F_{B}^\downarrow. +\f] +One atom of \f$ A \f$ recombines with one atom of \f$ B \f$ on the surface to produce a molecule \f$ AB \f$. +In other words the number of atoms of \f$ A \f$ that recombine into \f$ AB \f$ should be equal to the number of \f$ B \f$ +atoms that recombine into \f$ AB \f$. This restriction should be explicitly imposed in order to conserve mass: +\f[ +\gamma_{AB} F_{A}^\downarrow = \gamma_{BA} F_{B}^\downarrow. +\f] +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 \f$ A \f$ atoms and the opposite. In practice the two recombination number fluxes, +\f$ \gamma_{AB} F_{A}^\downarrow \f$ and \f$ \gamma_{BA} F_{B}^\downarrow \f$, 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 +\f[ +A + A \rightarrow A_2 \quad \text{and} \quad B + B \rightarrow B_2 +\f] +the gammas should be further constrained, as: +\f[ +0 \le \gamma_{A} + \gamma_{AB} \le 1 \quad \text{and} \quad 0 \le \gamma_{B} + \gamma_{BA} \le 1, +\f] +in order for mass to be conserved. This approach is consistent with similar approaches +considering the catalytic recombination occurring in Martian atmospheres, where \f$ O \f$ can recombine +into both \f$ O_2 \f$ and \f$ CO_2 \f$ due to catalytic reactions. + +An example input file of catalytic reactions in air including the formation of \f$ NO \f$ are below. +\code{.xml} + + + + + + + + + N:.001 + + + + O:.001 + + + + N:.002 O:.003 + + + + +\endcode + +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. +\code{.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. + + + + +\endcode + +An example of ablation model can be seen below. +\code{.xml} + + + + + + + + + + + + + + N:0.003 + + + + + + + + N:0.001 + + + + +\endcode + +#### Carbon Oxidation Model + + +The first ablation reaction presented above is the oxidation +of the solid carbon by atomic oxygen. The reaction reads as: + +\f[ +C_{\text{b}} + O \rightarrow CO +\f] + +and is exothermic, releasing \f$ 3.74 \f$ eV per molecule produced. Its reaction +rate coefficient is given by defining a recombination probability +\f$ \gamma_{CO} \f$. This probability is an Arrhenius type function of temperature +and is given by the formula: +\f[ +\gamma_{CO} = 0.63 \exp(-1160/T_{\text{s}}). +\f] +Carbon oxidation with molecular oxygen is also possible +(\f$ C_{\text{b}} + O_2 \rightarrow CO + O \f$), but it is considered a less significant +process. + +#### Carbon Nitridation Model + + +The second shown in the example is carbon nitridation, +\f[ +C_{\text{b}} + N \rightarrow CN, +\f] +an exothermic reaction with \f$ 0.35 \f$ eV of energy released per reacting atom. +The reaction rate is given using a constant recombination probability +\f$ \gamma_{CN} \f$ 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 \f$ C_3 \f$ is considered here. +It should be noted that this type of reactions are invertible, with formula: +\f[ +3C_\text{b} \rightleftharpoons C_3. +\f] +The chemical production rate of this reaction is equal to: +\f[ +\dot{\omega}_{\text{subl},C_3} = \beta_{C_3} ( \rho_{C_3, \text{equil}} - \rho_{C_3}) + \sqrt{ \frac{ k_B T_\text{s}}{2 \pi m_{C_3}}}. +\f] +The equilibrium partial density of the \f$ C_3 \f$ species is obtained from the +saturated vapor pressure of carbon, which is equal to +\f[ +p_{C_3, \text{sat}} = c \exp (-T_{\text{act}}/T_{\text{w}}) +\f] +with \f$ \beta_{C_3} \f$ being the evaporation coefficient, \f$ c \f$ the pre-exponential coefficient, +and \f$ T_{\text{act}} \f$ 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: +\f[ +\mathbf{j}_i \cdot \mathbf{n} = \dot{\omega}_{\text{cat},i}, +\f] +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: +\f[ +[ \rho_i (\mathbf{u}_{\text{g}} - \mathbf{u}_{\text{r}} ) + \mathbf{j}_i - \mathbf{F}_{\text{b},i} ] +\cdot \mathbf{n} = \dot{\omega}_{\text{s},i}, +\f] +with +\f$ \dot{\omega}_{\text{s},i} = \dot{\omega}_{\text{cat},i} + \dot{\omega}_{\text{abl},i} \f$ and +term \f$ \mathbf{F}_{\text{b},i} \f$ is the flux of species \f$ i \f$ 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 \f$ \mathbf{u}_{\text{r}} \f$. + +### Surface Energy Balance Solver + + +In order to determine the surface temperature a surface energy balance should be solved along +with the mass balances. It takes the form: +\f[ +[ \rho ( \mathbf{u}_{\text{g}} - \mathbf{u}_{\text{r}} ) H + \mathbf{q}_{\text{g}} + - \mathbf{F}_{\text{b},{e}} ] \cdot \mathbf{n} = + \dot{\Omega}_{\text{s},e} +\f] +where \f$ \mathbf{q}_{\text{g}} \f$ is the heat flux to the gas phase equal to: +\f[ + \mathbf{q}_{\text{g}} = -\lambda \nabla T + \sum_{i = 1}^{n_\text{s}} \mathbf{j}_i h_i. +\f] +The radiative heat flux can be also included and will be discussed along with the surface +radiation. + +Term \f$ \mathbf{F}_{\text{s},{e}} \f$ 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 \f$ \mathbf{q}_{\text{b,cond}} \f$, the second one is the +enthalpy entering the interface due to the movement of the surface with the recession velocity, +\f$ \mathbf{u}_{\text{r}} \rho_{\text{s}} h_{\text{s}} \f$ 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 \f$ \mathbf{u}_{\text{p}} \rho_{\text{p}} h_{\text{p}} \f$. The subscript p symbolizes +the pyrolysis gas properties, with the \f$ \rho_{\text{p}} h_{\text{p}} \f$ being actually the +sum \f$ \sum_{i=i}^{n_\text{s}} = \rho_i h_i \f$ for the pyrolysis gas densities. The surface +enthalpy \f$ h_s \f$ 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. +\code{.xml} + + + + + + + + + + + + + + + + + + N:0.003 + + + + + + + + N:0.001 + + + + + +\endcode + +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 \f$ T_{\text{s}} \f$ emitting energy following the Stefan-Boltzmann law, +\f[ +\dot{\Omega}_{\text{s},e} = \dot{\Omega}_{\text{rad}} = \epsilon \left( \sigma T_{\text{s}}^4 +- q_{\text{rad},g} \right) +\f] +with \f$ \sigma \f$ being the Stefan-Boltzmann constant and \f$ \epsilon \f$ the emissivity +of the surface. \f$ T_{\text{env}} \f$ is the surrounding environment temperature which can +be used to simulate the far field radiative heat flux based on the formula +\f$ \sigma T_{\text{env}}^4 \f$, which replaces the term \f$ q_{\text{rad},g} \f$. +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: +\f[ +q_{\text{cond}}^{\text{SS}} = -\mathbf{u}_{\text{r}}(\rho_\text{v} h_{\text{v}} - \rho_\text{s} h_{\text{s}}). +\f] +By replacing the formula above in the surface energy balance and after simplifications one gets: +\f[ +[ +\rho ( \mathbf{u}_{\text{g}} - \mathbf{u}_{\text{r}} ) H + \mathbf{q}_{\text{g}} ++ \mathbf{u}_{\text{r}}\rho_\text{v} h_{\text{v}} +] \cdot \mathbf{n} = +\dot{\Omega}_{\text{s},e}, +\f] +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 +\f[ \phi = ( \frac{\rho_\text{v}}{\rho_\text{s}} - 1) \f] 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. +@code{.xml} + + + + + ... + + +@endcode + +Below an example code of how to use the gas surface interaction features +of Mutation is presented. +@code{.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()); +@endcode + +Note that this example is not supposed to compile or run, but is there to indicate +the most important features of the library. + +*/ + + From 0eda4578707d5de812b221d533e882e2eff8e6e5 Mon Sep 17 00:00:00 2001 From: grgbellasvki <44604564+grgbellasvki@users.noreply.github.com> Date: Thu, 28 Mar 2019 12:10:22 +0100 Subject: [PATCH 2/9] Committing documentation before spellchecking --- docs/input-files.md | 255 ++++++++++++++++++++++---------------------- 1 file changed, 127 insertions(+), 128 deletions(-) diff --git a/docs/input-files.md b/docs/input-files.md index 26c77cdb..f3d0d589 100644 --- a/docs/input-files.md +++ b/docs/input-files.md @@ -22,7 +22,14 @@ - [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 @@ -979,7 +986,7 @@ three dimensional fluxes reduce to the normal fluxes at the interface while the source terms, such as chemical reactions, go to zero; only the surface sources remain. Generally, the set of balance equations can be written as: -(**Fg** - **Fb**) . **n** = Ωs, +(**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. @@ -1031,57 +1038,55 @@ means that all the particles of species A impinging the surface recombine at the 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: -\f[ -\dot{\omega}_{A,\text{cat}} = \gamma_{A} m_{A} F_{A}^\downarrow \quad \text{and} \quad -\dot{\omega}_{A_2,\text{cat}} = - \gamma_{A} m_{A} F_{A}^\downarrow -\f] -with \f$ m_{A} \f$ being the mass of species \f$ A \f$. + +Ω_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$ F_{A}^\downarrow \f$, is equal to -\f[ -F_{A}^\downarrow = n_{A} \sqrt{\frac{k_B T_{\text{s}}}{2 \pi m_{A}}} -\f] +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. +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 -\f[ -A + B \rightarrow AB, -\f] -with the chemical rate for the produced molecule \f$ AB \f$ equal to: -\f[ -\dot{\omega}_{AB,\text{cat}} = - \gamma_{AB} m_{A} F_{A}^\downarrow - \gamma_{BA} m_{B} F_{B}^\downarrow. -\f] -One atom of \f$ A \f$ recombines with one atom of \f$ B \f$ on the surface to produce a molecule \f$ AB \f$. -In other words the number of atoms of \f$ A \f$ that recombine into \f$ AB \f$ should be equal to the number of \f$ B \f$ -atoms that recombine into \f$ AB \f$. This restriction should be explicitly imposed in order to conserve mass: -\f[ -\gamma_{AB} F_{A}^\downarrow = \gamma_{BA} F_{B}^\downarrow. -\f] + +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 \f$ A \f$ atoms and the opposite. In practice the two recombination number fluxes, -\f$ \gamma_{AB} F_{A}^\downarrow \f$ and \f$ \gamma_{BA} F_{B}^\downarrow \f$, are compared and the limiting one +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 -\f[ -A + A \rightarrow A_2 \quad \text{and} \quad B + B \rightarrow B_2 -\f] -the gammas should be further constrained, as: -\f[ -0 \le \gamma_{A} + \gamma_{AB} \le 1 \quad \text{and} \quad 0 \le \gamma_{B} + \gamma_{BA} \le 1, -\f] +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 \f$ O \f$ can recombine -into both \f$ O_2 \f$ and \f$ CO_2 \f$ due to catalytic reactions. +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. -An example input file of catalytic reactions in air including the formation of \f$ NO \f$ are below. -\code{.xml} +```xml @@ -1103,7 +1108,7 @@ An example input file of catalytic reactions in air including the formation of \ -\endcode +``` 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 @@ -1111,7 +1116,8 @@ 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. -\code{.xml} + +```xml @@ -1153,10 +1159,9 @@ Mutation. -\endcode - +``` An example of ablation model can be seen below. -\code{.xml} +```xml @@ -1183,7 +1188,7 @@ An example of ablation model can be seen below. -\endcode +``` #### Carbon Oxidation Model @@ -1191,53 +1196,50 @@ An example of ablation model can be seen below. The first ablation reaction presented above is the oxidation of the solid carbon by atomic oxygen. The reaction reads as: -\f[ -C_{\text{b}} + O \rightarrow CO -\f] +C_b + O -> CO -and is exothermic, releasing \f$ 3.74 \f$ eV per molecule produced. Its reaction +and is exothermic, releasing 3.74 eV per molecule produced. Its reaction rate coefficient is given by defining a recombination probability -\f$ \gamma_{CO} \f$. This probability is an Arrhenius type function of temperature +γ_CO. This probability is an Arrhenius type function of temperature and is given by the formula: -\f[ -\gamma_{CO} = 0.63 \exp(-1160/T_{\text{s}}). -\f] + +γ_CO = 0.63 exp(-1160 / T). + Carbon oxidation with molecular oxygen is also possible -(\f$ C_{\text{b}} + O_2 \rightarrow CO + O \f$), but it is considered a less significant +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, -\f[ -C_{\text{b}} + N \rightarrow CN, -\f] -an exothermic reaction with \f$ 0.35 \f$ eV of energy released per reacting atom. + +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 -\f$ \gamma_{CN} \f$ like in the catalytic case. +γ_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 \f$ C_3 \f$ is considered here. +change processes like sublimation. The production of C3 is considered here. It should be noted that this type of reactions are invertible, with formula: -\f[ -3C_\text{b} \rightleftharpoons C_3. -\f] + +3C_b -> C3. + The chemical production rate of this reaction is equal to: -\f[ -\dot{\omega}_{\text{subl},C_3} = \beta_{C_3} ( \rho_{C_3, \text{equil}} - \rho_{C_3}) - \sqrt{ \frac{ k_B T_\text{s}}{2 \pi m_{C_3}}}. -\f] -The equilibrium partial density of the \f$ C_3 \f$ species is obtained from the + +Ω_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 -\f[ -p_{C_3, \text{sat}} = c \exp (-T_{\text{act}}/T_{\text{w}}) -\f] -with \f$ \beta_{C_3} \f$ being the evaporation coefficient, \f$ c \f$ the pre-exponential coefficient, -and \f$ T_{\text{act}} \f$ the activation temperature. + +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. @@ -1267,16 +1269,16 @@ 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: -\f[ -\mathbf{j}_i \cdot \mathbf{n} = \dot{\omega}_{\text{cat},i}, -\f] + +**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 +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 @@ -1298,48 +1300,46 @@ and occur only the surface of the material. In cases where the material is porou 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: -\f[ -[ \rho_i (\mathbf{u}_{\text{g}} - \mathbf{u}_{\text{r}} ) + \mathbf{j}_i - \mathbf{F}_{\text{b},i} ] -\cdot \mathbf{n} = \dot{\omega}_{\text{s},i}, -\f] -with -\f$ \dot{\omega}_{\text{s},i} = \dot{\omega}_{\text{cat},i} + \dot{\omega}_{\text{abl},i} \f$ and -term \f$ \mathbf{F}_{\text{b},i} \f$ is the flux of species \f$ i \f$ entering the interface due + +ρ_i (**u**_g - **u**_r) + **j**_i - **F**_b_i ] +\cdot **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 \f$ \mathbf{u}_{\text{r}} \f$. +orders of magnitude lower than the gas velocity **u**_r. -### Surface Energy Balance Solver +### 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: -\f[ -[ \rho ( \mathbf{u}_{\text{g}} - \mathbf{u}_{\text{r}} ) H + \mathbf{q}_{\text{g}} - - \mathbf{F}_{\text{b},{e}} ] \cdot \mathbf{n} = - \dot{\Omega}_{\text{s},e} -\f] -where \f$ \mathbf{q}_{\text{g}} \f$ is the heat flux to the gas phase equal to: -\f[ - \mathbf{q}_{\text{g}} = -\lambda \nabla T + \sum_{i = 1}^{n_\text{s}} \mathbf{j}_i h_i. -\f] + +( ρ ( **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$ \mathbf{F}_{\text{s},{e}} \f$ describes the energy exchanged between the interface +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 \f$ \mathbf{q}_{\text{b,cond}} \f$, the second 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, -\f$ \mathbf{u}_{\text{r}} \rho_{\text{s}} h_{\text{s}} \f$ and the third one appears only in cases of +**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 \f$ \mathbf{u}_{\text{p}} \rho_{\text{p}} h_{\text{p}} \f$. The subscript p symbolizes -the pyrolysis gas properties, with the \f$ \rho_{\text{p}} h_{\text{p}} \f$ being actually the -sum \f$ \sum_{i=i}^{n_\text{s}} = \rho_i h_i \f$ for the pyrolysis gas densities. The surface -enthalpy \f$ h_s \f$ is an input to the code, as an attribute to the surface_properties element. +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. -\code{.xml} +```xml @@ -1379,20 +1379,19 @@ An example input file for solving both mass and energy balance can be seen below -\endcode +``` 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 \f$ T_{\text{s}} \f$ emitting energy following the Stefan-Boltzmann law, -\f[ -\dot{\Omega}_{\text{s},e} = \dot{\Omega}_{\text{rad}} = \epsilon \left( \sigma T_{\text{s}}^4 -- q_{\text{rad},g} \right) -\f] -with \f$ \sigma \f$ being the Stefan-Boltzmann constant and \f$ \epsilon \f$ the emissivity -of the surface. \f$ T_{\text{env}} \f$ is the surrounding environment temperature which can +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 -\f$ \sigma T_{\text{env}}^4 \f$, which replaces the term \f$ q_{\text{rad},g} \f$. + σ 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, @@ -1406,23 +1405,23 @@ semi-infinite surface. By writing the steady-state energy equation for the solid 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: -\f[ -q_{\text{cond}}^{\text{SS}} = -\mathbf{u}_{\text{r}}(\rho_\text{v} h_{\text{v}} - \rho_\text{s} h_{\text{s}}). -\f] + +**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: -\f[ + [ -\rho ( \mathbf{u}_{\text{g}} - \mathbf{u}_{\text{r}} ) H + \mathbf{q}_{\text{g}} -+ \mathbf{u}_{\text{r}}\rho_\text{v} h_{\text{v}} -] \cdot \mathbf{n} = -\dot{\Omega}_{\text{s},e}, -\f] +ρ ( **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 -\f[ \phi = ( \frac{\rho_\text{v}}{\rho_\text{s}} - 1) \f] and is the attribute +φ = (ρ_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 @@ -1448,7 +1447,8 @@ adiabatic boundary condition for the internal energy modes. 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. -@code{.xml} + +```xml @@ -1456,11 +1456,12 @@ mechanism file. ... -@endcode +``` Below an example code of how to use the gas surface interaction features of Mutation is presented. -@code{.cpp} + +```cpp const int set_state_with_rhoi_T = 1; MixtureOptions opts("mixture_name"); @@ -1496,11 +1497,9 @@ of Mutation is presented. // Getting number of reactions and surface reaction rates per reaction. int m_nr = mix.nSurfaceReactions(); mix.getSurfaceReactionRatesPerReaction(wdot_reac.data()); -@endcode +``` Note that this example is not supposed to compile or run, but is there to indicate the most important features of the library. -*/ - From 2657e04c4e34c1c55cecf78995fcd272fd753304 Mon Sep 17 00:00:00 2001 From: grgbellasvki <44604564+grgbellasvki@users.noreply.github.com> Date: Thu, 28 Mar 2019 13:03:47 +0100 Subject: [PATCH 3/9] Fixed typos in GSI documentation --- docs/input-files.md | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/docs/input-files.md b/docs/input-files.md index f3d0d589..ccb1eab1 100644 --- a/docs/input-files.md +++ b/docs/input-files.md @@ -434,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 @@ -980,7 +980,7 @@ momentum and energy, which are imposed as boundary values for the full Navier-St 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, \f$\partial \mathcal / \partial t\f$, +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. @@ -1073,7 +1073,7 @@ atoms that recombine into AB. This restriction should be explicitly imposed in o 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 +γ_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, @@ -1301,11 +1301,10 @@ 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 ] -\cdot **n** = Ω_i, +ρ_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 +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 @@ -1334,7 +1333,7 @@ enthalpy entering the interface due to the movement of the surface with the rece **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 +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. @@ -1386,7 +1385,7 @@ 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) +Ω_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 @@ -1410,11 +1409,7 @@ is the virgin material (subscript v), the steady state heat flux is given by the 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, +(ρ ( **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. From c22dda78fb795147efb1cff4fda87bd5a5cfd00c Mon Sep 17 00:00:00 2001 From: MicheleCapriati Date: Mon, 8 Apr 2019 11:30:08 +0200 Subject: [PATCH 4/9] new energy terms --- src/gsi/CMakeLists.txt | 3 +++ src/gsi/GSIRateManagerGamma.cpp | 3 +++ src/gsi/GasFourierHeatFluxCalculator.cpp | 21 ++++++++++++++++++++- src/gsi/GasFourierHeatFluxCalculator.h | 13 +++++++++++++ src/gsi/GasSurfaceInteraction.cpp | 11 +++++++++++ src/gsi/GasSurfaceInteraction.h | 7 +++++++ src/gsi/Surface.h | 14 ++++++++++++++ tests/CMakeLists.txt | 1 + tests/test_gsi_mass_energy.cpp | 2 +- 9 files changed, 73 insertions(+), 2 deletions(-) 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/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/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2512e308..0bf73d89 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/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 +} From 2b9d897ce414515ec31778a3472ccd309685f1fb Mon Sep 17 00:00:00 2001 From: MicheleCapriati Date: Mon, 8 Apr 2019 11:36:51 +0200 Subject: [PATCH 5/9] new energy terms --- src/gsi/GSIRateLawGammaTTv.cpp | 101 ++++ src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp | 452 ++++++++++++++++++ src/gsi/SurfaceInelastic.cpp | 108 +++++ src/gsi/SurfaceInelastic.h | 86 ++++ tests/data/gsi/gsi_N2_C_ablation.xml | 28 ++ .../nitrogen-ions-carbon_9_olynick99_mech.xml | 52 ++ ...itridation_nitrogen-ions-carbon_9_RRHO.xml | 15 + tests/test_gsi_mass_energyTTv.cpp | 157 ++++++ 8 files changed, 999 insertions(+) create mode 100644 src/gsi/GSIRateLawGammaTTv.cpp create mode 100644 src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp create mode 100644 src/gsi/SurfaceInelastic.cpp create mode 100644 src/gsi/SurfaceInelastic.h create mode 100644 tests/data/gsi/gsi_N2_C_ablation.xml create mode 100644 tests/data/mechanisms/nitrogen-ions-carbon_9_olynick99_mech.xml create mode 100644 tests/data/mixtures/sebTTv_nitridation_nitrogen-ions-carbon_9_RRHO.xml create mode 100644 tests/test_gsi_mass_energyTTv.cpp 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/SurfaceBalanceSolverMassEnergyTTv.cpp b/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp new file mode 100644 index 00000000..ea954e2c --- /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(true); + 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(mv_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(); + + // Energy + m_thermo.getEnthalpiesMass(mv_hi.data()); + //double hmix = m_thermo.mixtureHMass(); + double hmix = mv_rhoi.dot(mv_hi.head(m_ns))/mv_rhoi.sum(); + + 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 entalphy of the mixture + 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); + } + +//============================================================================== + + void updateJacobian(Eigen::VectorXd& v_X) + { + // 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; + + // 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) = (mv_f-mv_f_unpert) / T_pert; + v_X(pos_E+1) = X_unpert; + for (int i=0; i(); + // 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..688ff217 --- /dev/null +++ b/src/gsi/SurfaceInelastic.cpp @@ -0,0 +1,108 @@ +/** + * @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) +{ + + //computing Vibrational traslational exchange + double T_tra = v_X(pos_E); + double T_vib = v_X(pos_E + 1); + Eigen::VectorXd h_tra(m_ns); + Eigen::VectorXd h_vib(m_ns); + Eigen::VectorXd h_el(m_ns); + 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); + double number_density = m_thermo.numberDensity(); + double inleastic_term = 0.; + double tol = 1.E-19; + if (m_eff_coll > tol) { + for(int i = 0; i < m_ns; ++i) { + //thermal speed of species i + double thermal_speed; + //if (i < m_thermo.hasElectrons()) thermal_speed = sqrt(RU*T_vib/(2.*PI*m_speciesMw(i))); + //else thermal_speed = sqrt(RU*T_tra/(2.*PI*m_speciesMw(i))); + 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 + double num_dens_i = v_X(i)*number_density; + + //characteristc time ^-1 + double one_over_tau = num_dens_i*thermal_speed; + //std::cout << "one_over_tau is" << one_over_tau << std::endl; + + //vibronic enthaply at Teq + double h_VE; + if (i < m_thermo.hasElectrons()) h_VE = (h_tra[i]*T_tra)*RU/m_speciesMw(i); //check it + else h_VE = (h_vib(i) + h_el(i))*T_tra*RU/m_speciesMw(i); + + inleastic_term += m_eff_coll*(h_VE-v_h(m_ns + i))*one_over_tau; + } + } + + //Compute vibrational chemical production + if (1.-m_beta > tol) 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.h b/src/gsi/SurfaceInelastic.h new file mode 100644 index 00000000..4c6a69a9 --- /dev/null +++ b/src/gsi/SurfaceInelastic.h @@ -0,0 +1,86 @@ +/** + * @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); + +//============================================================================== + +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/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_energyTTv.cpp b/tests/test_gsi_mass_energyTTv.cpp new file mode 100644 index 00000000..e6b6ca0c --- /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)); + } +} From d9f2887c84ad50e8adb3b1736accf5d110e5cc9e Mon Sep 17 00:00:00 2001 From: MicheleCapriati Date: Mon, 8 Apr 2019 20:35:20 +0200 Subject: [PATCH 6/9] fix minor stuff --- src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp | 9 +++++---- src/gsi/SurfaceInelastic.cpp | 2 +- tests/test_gsi_mass_energyTTv.cpp | 2 +- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp b/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp index ea954e2c..faa3c9f8 100644 --- a/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp +++ b/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp @@ -302,6 +302,7 @@ class SurfaceBalanceSolverMassEnergyTTv : 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++){ @@ -312,7 +313,8 @@ class SurfaceBalanceSolverMassEnergyTTv : updateFunction(v_X); // Update Jacobian column - m_jac.col(i_ns) = (mv_f-mv_f_unpert) / pert; + //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; @@ -330,10 +332,9 @@ class SurfaceBalanceSolverMassEnergyTTv : X_unpert = v_X(pos_E+1); v_X(pos_E+1) += T_pert; updateFunction(v_X); - m_jac.col(pos_E+1) = (mv_f-mv_f_unpert) / T_pert; + 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; - for (int i=0; i tol) { + double thermal_speed; for(int i = 0; i < m_ns; ++i) { //thermal speed of species i - double thermal_speed; //if (i < m_thermo.hasElectrons()) thermal_speed = sqrt(RU*T_vib/(2.*PI*m_speciesMw(i))); //else thermal_speed = sqrt(RU*T_tra/(2.*PI*m_speciesMw(i))); if (i < m_thermo.hasElectrons()) thermal_speed = sqrt(T_vib)*m_therm_vel_over_T(i); diff --git a/tests/test_gsi_mass_energyTTv.cpp b/tests/test_gsi_mass_energyTTv.cpp index e6b6ca0c..91b7a2c8 100644 --- a/tests/test_gsi_mass_energyTTv.cpp +++ b/tests/test_gsi_mass_energyTTv.cpp @@ -140,7 +140,7 @@ TEST_CASE 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); + //mix.getInelasticTerm(xi_s, v_hi, wdot, inelastic); // Building balance functions VectorXd F(neq); From 0679355580ef0fa9094a01202be3ede4edd59fc4 Mon Sep 17 00:00:00 2001 From: MicheleCapriati Date: Wed, 17 Jul 2019 11:50:29 +0200 Subject: [PATCH 7/9] inelastic fixed --- src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp | 15 ++--- src/gsi/SurfaceInelastic.cpp | 66 +++++++++++-------- src/gsi/SurfaceInelastic.h | 3 +- tests/CMakeLists.txt | 2 +- 4 files changed, 47 insertions(+), 39 deletions(-) diff --git a/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp b/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp index faa3c9f8..127f59f4 100644 --- a/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp +++ b/src/gsi/SurfaceBalanceSolverMassEnergyTTv.cpp @@ -122,7 +122,7 @@ class SurfaceBalanceSolverMassEnergyTTv : // Setup NewtonSolver setMaxIterations(10); - setWriteConvergenceHistory(true); + setWriteConvergenceHistory(false); setEpsilon(m_tol); } @@ -241,7 +241,7 @@ class SurfaceBalanceSolverMassEnergyTTv : void updateFunction(Eigen::VectorXd& v_X) { - applyTolerance(mv_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. @@ -271,10 +271,9 @@ class SurfaceBalanceSolverMassEnergyTTv : mv_surf_reac_rates); mv_f.head(m_ns) += mv_rhoi*mass_blow/mv_rhoi.sum(); - // Energy + // Total Energy m_thermo.getEnthalpiesMass(mv_hi.data()); - //double hmix = m_thermo.mixtureHMass(); - double hmix = mv_rhoi.dot(mv_hi.head(m_ns))/mv_rhoi.sum(); + double hmix = m_thermo.mixtureHMass(); mv_f(pos_E) += mv_hi.head(m_ns).dot(mv_Vdiff.cwiseProduct(mv_rhoi)); @@ -286,7 +285,7 @@ class SurfaceBalanceSolverMassEnergyTTv : if (mp_surf_rad != NULL) mv_f(pos_E) -= mp_surf_rad->surfaceNetRadiativeHeatFlux(); - //vibronic entalphy of the mixture + // Vibronic Energy double hVMix = mv_rhoi.dot(mv_hi.tail(m_ns))/mv_rhoi.sum(); mv_f(pos_E+1) += mass_blow*hVMix; @@ -294,8 +293,8 @@ class SurfaceBalanceSolverMassEnergyTTv : 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_f(pos_E+1) += + mp_surf_inelastic->surfaceInelasticTerm(v_X, mv_hi, mv_surf_reac_rates, mv_rhoi); } //============================================================================== diff --git a/src/gsi/SurfaceInelastic.cpp b/src/gsi/SurfaceInelastic.cpp index 97996c18..7527ab4d 100644 --- a/src/gsi/SurfaceInelastic.cpp +++ b/src/gsi/SurfaceInelastic.cpp @@ -61,46 +61,54 @@ SurfaceInelastic::~SurfaceInelastic(){} //============================================================================== -double SurfaceInelastic::surfaceInelasticTerm(const VectorXd& v_X, const VectorXd& v_h, const VectorXd& chem_souce) +double SurfaceInelastic::surfaceInelasticTerm(const VectorXd& v_X, const VectorXd& v_h, const VectorXd& chem_souce, const VectorXd& v_rhoi) { //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); - 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); + + 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(); - double inleastic_term = 0.; - double tol = 1.E-19; - if (m_eff_coll > tol) { - double thermal_speed; - for(int i = 0; i < m_ns; ++i) { - //thermal speed of species i - //if (i < m_thermo.hasElectrons()) thermal_speed = sqrt(RU*T_vib/(2.*PI*m_speciesMw(i))); - //else thermal_speed = sqrt(RU*T_tra/(2.*PI*m_speciesMw(i))); - 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 - double num_dens_i = v_X(i)*number_density; - - //characteristc time ^-1 - double one_over_tau = num_dens_i*thermal_speed; - //std::cout << "one_over_tau is" << one_over_tau << std::endl; - - //vibronic enthaply at Teq - double h_VE; - if (i < m_thermo.hasElectrons()) h_VE = (h_tra[i]*T_tra)*RU/m_speciesMw(i); //check it - else h_VE = (h_vib(i) + h_el(i))*T_tra*RU/m_speciesMw(i); - - inleastic_term += m_eff_coll*(h_VE-v_h(m_ns + i))*one_over_tau; - } - } + + 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 - if (1.-m_beta > tol) inleastic_term += (1-m_beta)*v_h(m_index)*chem_souce(m_index); + inleastic_term -= (1.-m_beta)*v_h(m_index)*chem_souce(m_index); return inleastic_term; } diff --git a/src/gsi/SurfaceInelastic.h b/src/gsi/SurfaceInelastic.h index 4c6a69a9..1a0830ea 100644 --- a/src/gsi/SurfaceInelastic.h +++ b/src/gsi/SurfaceInelastic.h @@ -62,7 +62,8 @@ class SurfaceInelastic /** * Function which returns the inelastic term. */ - double surfaceInelasticTerm(const Eigen::VectorXd& v_X, const Eigen::VectorXd& v_h, const Eigen::VectorXd& chem_souce); + double surfaceInelasticTerm(const Eigen::VectorXd& v_X, const Eigen::VectorXd& v_h, const Eigen::VectorXd& chem_souce, + const Eigen::VectorXd& v_rhoi); //============================================================================== diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0bf73d89..8b7e9436 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,7 +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_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 From 95a8dd9ef0aece78f14245dae4b533926abf62ab Mon Sep 17 00:00:00 2001 From: MicheleCapriati Date: Wed, 17 Jul 2019 11:59:02 +0200 Subject: [PATCH 8/9] new OmegaI --- src/transfer/OmegaI.cpp | 151 ++++++++++++++++++++++------------------ 1 file changed, 85 insertions(+), 66 deletions(-) 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 From cf244fe9bda71f4f98da22475e96dd950a92ac67 Mon Sep 17 00:00:00 2001 From: MicheleCapriati Date: Thu, 31 Oct 2019 10:08:06 +0100 Subject: [PATCH 9/9] imroved chemical term --- ...eBalanceSolverMassEnergyTTv_with_VDiff.cpp | 452 ++++++++++++++++++ src/gsi/SurfaceInelastic.cpp | 5 +- src/gsi/SurfaceInelastic_with_VDiff.cpp | 117 +++++ src/gsi/SurfaceInelastic_with_VDiff.h | 87 ++++ 4 files changed, 658 insertions(+), 3 deletions(-) create mode 100644 src/gsi/SurfaceBalanceSolverMassEnergyTTv_with_VDiff.cpp create mode 100644 src/gsi/SurfaceInelastic_with_VDiff.cpp create mode 100644 src/gsi/SurfaceInelastic_with_VDiff.h 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 index 7527ab4d..5d6fe040 100644 --- a/src/gsi/SurfaceInelastic.cpp +++ b/src/gsi/SurfaceInelastic.cpp @@ -48,7 +48,6 @@ SurfaceInelastic::SurfaceInelastic( 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.); @@ -65,9 +64,9 @@ double SurfaceInelastic::surfaceInelasticTerm(const VectorXd& v_X, const VectorX { //computing Vibrational traslational exchange + double inleastic_term = 0.; 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; @@ -108,7 +107,7 @@ double SurfaceInelastic::surfaceInelasticTerm(const VectorXd& v_X, const VectorX } //Compute vibrational chemical production - inleastic_term -= (1.-m_beta)*v_h(m_index)*chem_souce(m_index); + inleastic_term += (1.-m_beta)*chem_souce.dot( v_h.head(m_ns) - v_h.tail(m_ns)); return inleastic_term; } 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