LCOV - code coverage report
Current view: top level - src/Model/Propellant - FullFeedback.H (source / functions) Coverage Total Hit
Test: coverage_merged.info Lines: 98.2 % 56 55
Test Date: 2025-08-12 17:45:17 Functions: 100.0 % 3 3

            Line data    Source code
       1              : //
       2              : // This SCP model implements the method described in the following paper
       3              : //
       4              : // .. bibliography::
       5              : //    :list: none
       6              : //    :filter: False
       7              : // 
       8              : //    meier2024diffuse
       9              : //
      10              : // This method imitates inheritance :code:`Propellant::Propellant` using static polymorphism
      11              : // 
      12              : 
      13              : #ifndef MODEL_PROPELLANT_FULLFEEDBACK_H
      14              : #define MODEL_PROPELLANT_FULLFEEDBACK_H
      15              : 
      16              : #include "Propellant.H"
      17              : 
      18              : namespace Model
      19              : {
      20              : namespace Propellant
      21              : {
      22              : class FullFeedback
      23              : {
      24              : public:
      25              :     static constexpr const char* name = "fullfeedback";
      26              : 
      27            3 :     FullFeedback() = default;
      28              :     FullFeedback(IO::ParmParse &pp, std::string name) 
      29              :     {
      30              :         pp_queryclass(name,*this);
      31              :     }
      32              : 
      33         2550 :     void set_pressure(const Set::Scalar a_P) 
      34              :     {
      35         2550 :         this->Pbar = a_P/Pref;
      36              :         
      37         2550 :         zeta_2 = zeta_fit_a + Pbar * zeta_fit_b;
      38              : 
      39         2550 :         if (arrhenius_dependency == 1) zeta_1 = zeta_2;
      40            0 :         else zeta_1 = zeta_0;
      41              : 
      42              :         // Pressure-dependent post calculation variables
      43         2550 :         k1 = a1 * Pbar + b1 - zeta_1 / zeta;
      44         2550 :         k2 = a2 * Pbar + b2 - zeta_1 / zeta;
      45         2550 :         k3 = 4.0 * log((c1*Pbar*Pbar + a3*Pbar + b3) - k1 / 2.0 - k2 / 2.0);
      46              : 
      47         2550 :         return;
      48              :     }
      49              :     
      50              :     AMREX_FORCE_INLINE
      51              :     Set::Scalar get_qdot(const Set::Scalar mdot, const Set::Scalar phi) 
      52              :     {
      53      1747840 :         Set::Scalar qflux = k1 * phi +
      54      1747840 :             k2 * (1.0 - phi) +
      55      1747840 :             (zeta_1 / zeta) * exp(k3 * phi * (1.0 - phi));
      56      1747840 :         Set::Scalar mlocal = (mlocal_ap) * phi + (mlocal_htpb) * (1.0 - phi) + mlocal_comb * phi * (1.0 - phi);
      57      1747840 :         Set::Scalar mdota = fabs(mdot);
      58      1747840 :         Set::Scalar mbase = tanh(4.0 * mdota / (mlocal));
      59              : 
      60      1747840 :         return mbase * qflux;
      61              :     }
      62              : 
      63              : 
      64              :     AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      65              :     Set::Scalar get_K (const Set::Scalar phi)
      66              :     {
      67      1747840 :         return k_ap * phi + k_htpb * (1.0 - phi);
      68              :     }
      69              : 
      70              :     AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      71              :     Set::Scalar get_rho (const Set::Scalar phi)
      72              :     {
      73      1747840 :         return rho_ap * phi + rho_htpb * (1.0 - phi); 
      74              :     }
      75              : 
      76              :     AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      77              :     Set::Scalar get_cp (const Set::Scalar phi)
      78              :     {
      79      1747840 :         return cp_ap * phi + cp_htpb * (1.0 - phi);
      80              :     }
      81              : 
      82              :     AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      83              :     Set::Scalar get_L (Set::Scalar phi, Set::Scalar T) 
      84              :     {
      85      1747840 :         Set::Scalar L = NAN;
      86      1747840 :         if (mob_ap == 1) L = m_ap * Pbar * exp(-E_ap / T) * phi;
      87      1747840 :         else L = m_ap * exp(-E_ap / T) * phi;
      88      1747840 :         L += m_htpb * exp(-E_htpb / T) * (1.0 - phi);
      89      1747840 :         if (T <= bound) L = 0.0;
      90      1747840 :         return L;
      91              :     }
      92              : 
      93              : 
      94            1 :     static void Parse(FullFeedback& value, IO::ParmParse& pp)
      95              :     {
      96              :         // AP/HTPB interface length
      97            4 :         pp.query_default("phi.zeta", value.zeta, "10.0_um", Unit::Length()); 
      98              :         // Reference interface length for heat integration
      99            4 :         pp.query_default("phi.zeta_0", value.zeta_0, "10.0_um", Unit::Length()); 
     100              :         
     101              :         // Constant offset for $\zeta$ adjustment with pressure
     102            4 :         pp.query_default("phi.zeta_fit_a",value.zeta_fit_a, "45.0_um", Unit::Length());
     103              :         // Scaling factors for $\zeta$ adjustment with pressure
     104            4 :         pp.query_default("phi.zeta_fit_b",value.zeta_fit_b, "-6.42_um", Unit::Length());
     105              :         
     106              :         // Surrogate heat flux model paramater - AP
     107            2 :         pp_query_required("a1", value.a1, Unit::Less()); 
     108              :         // Surrogate heat flux model paramater - HTPB
     109            2 :         pp_query_required("a2", value.a2,Unit::Less()); 
     110              :         // Surrogate heat flux model paramater - Total
     111            2 :         pp_query_required("a3", value.a3,Unit::Less()); 
     112              :         // Surrogate heat flux model paramater - AP
     113            2 :         pp_query_required("b1", value.b1, Unit::Less()); 
     114              :         // Surrogate heat flux model paramater - HTPB
     115            2 :         pp_query_required("b2", value.b2, Unit::Less()); 
     116              :         // Surrogate heat flux model paramater - Total
     117            2 :         pp_query_required("b3", value.b3,Unit::Less()); 
     118              :         // Surrogate heat flux model paramater - Total
     119            2 :         pp_query_required("c1", value.c1, Unit::Less()); 
     120              : 
     121              :         // Whether to use pressure to determined the reference Zeta 
     122            2 :         pp_query_default("pressure.dependency", value.arrhenius_dependency, true); 
     123              : 
     124              :         // Reference pressure for pressure nondimensionalization
     125            4 :         pp_query_default("Pref",value.Pref,"1_MPa",Unit::Pressure());
     126              :         
     127              :         // AP volumetric mass flux reference value 
     128            4 :         pp.query_default("mlocal_ap", value.mlocal_ap, "0.0", Unit::Density()/Unit::Time()); 
     129              :         // combined volumetric mass flux reference value 
     130            4 :         pp_query_default("mlocal_comb", value.mlocal_comb, "0.0", Unit::Density()/Unit::Time()); 
     131              :         // HTPB volumetric mass flux reference value
     132            4 :         pp.query_default("mlocal_htpb",value.mlocal_htpb, "5000_kg/m^3/s", Unit::Density()/Unit::Time());
     133              : 
     134              :         // Set::Scalar massfraction = 0.8;
     135              :         // value.mlocal_htpb = 685000.0 - 850e3 * massfraction; // Replicating but needs to be revised!
     136              :         // Util::Message(INFO,value.mlocal_htpb);
     137              : 
     138              :         // AP Density
     139            2 :         pp_query_required("rho_ap", value.rho_ap, Unit::Density());
     140              :         // HTPB Density
     141            2 :         pp_query_required("rho_htpb", value.rho_htpb, Unit::Density());
     142              :         // AP Thermal Conductivity
     143            2 :         pp_query_required("k_ap", value.k_ap, Unit::ThermalConductivity());
     144              :         // HTPB Thermal Conductivity
     145            2 :         pp_query_required("k_htpb", value.k_htpb, Unit::ThermalConductivity());
     146              :         // AP Specific Heat
     147            2 :         pp_query_required("cp_ap", value.cp_ap,Unit::SpecificHeatCapacity());
     148              :         //HTPB Specific Heat
     149            2 :         pp_query_required("cp_htpb", value.cp_htpb,Unit::SpecificHeatCapacity()); 
     150              : 
     151              :         // AP Pre-exponential factor for Arrhenius Law
     152            2 :         pp_query_required("m_ap", value.m_ap,Unit::Velocity()); 
     153              :         // HTPB Pre-exponential factor for Arrhenius Law
     154            2 :         pp_query_required("m_htpb", value.m_htpb,Unit::Velocity()); 
     155              :         // COMBINED AP Activation Energy, divided by Rg, for Arrhenius Law (has units of temperature)
     156            2 :         pp_query_required("E_ap", value.E_ap,Unit::Temperature()); 
     157              :         // COMBINED HTPB Activation Energy, divided by Rg, for Arrhenius Law(has units of temperature)
     158            2 :         pp_query_required("E_htpb", value.E_htpb,Unit::Temperature()); 
     159              :         // Whether to include pressure to the arrhenius law [??]n
     160            2 :         pp_query_default("mob_ap", value.mob_ap,false); 
     161              :         // HTPB Activation Energy for Arrhenius Law
     162            4 :         pp_query_default("bound", value.bound,"0.0", Unit::Temperature()); 
     163            1 :     }
     164              : 
     165              :     // IO Variables
     166              :     Set::Scalar zeta = NAN, zeta_0 = NAN; 
     167              :     Set::Scalar a1 = NAN, a2 = NAN, a3 = NAN;
     168              :     Set::Scalar b1 = NAN, b2 = NAN, b3 = NAN;
     169              :     Set::Scalar c1 = NAN;
     170              :     //Set::Scalar h1 = NAN, h2 = NAN;
     171              :     bool arrhenius_dependency = false;
     172              :     Set::Scalar mlocal_ap = NAN, mlocal_comb = NAN, mlocal_htpb = NAN;
     173              : 
     174              :     //Thermal properties
     175              :     Set::Scalar k_ap = NAN, k_htpb = NAN;
     176              :     Set::Scalar rho_ap = NAN, rho_htpb = NAN;
     177              :     Set::Scalar cp_ap = NAN, cp_htpb = NAN;
     178              : 
     179              :     // Regression variables
     180              :     int mob_ap = -1;
     181              :     Set::Scalar m_ap = NAN, m_htpb = NAN, m_comb = NAN;
     182              :     Set::Scalar E_ap = NAN, E_htpb = NAN;
     183              :     Set::Scalar bound = NAN;
     184              :     Set::Scalar Pref = NAN;
     185              : 
     186              :     // Calculated variables
     187              :     Set::Scalar Pbar = NAN;
     188              :     Set::Scalar zeta_1 = NAN, zeta_2 = NAN;
     189              :     Set::Scalar zeta_fit_a = NAN, zeta_fit_b = NAN;
     190              :     Set::Scalar k1=NAN, k2=NAN, k3=NAN, k4=NAN;
     191              : };
     192              : 
     193              : }
     194              : }
     195              : 
     196              : 
     197              : 
     198              : #endif
     199              : 
     200              : 
        

Generated by: LCOV version 2.0-1