double_trap_kinetics.h

Go to the documentation of this file.

//---------------------------------------------------------------------------
//
//    FCST: Fuel Cell Simulation Toolbox
//
//    Copyright (C) 2010-13 by Energy Systems Design Laboratory, University of Alberta
//
//    This software is distributed under the MIT License.
//    For more information, see the README file in /doc/LICENSE
//
//    - Class: double_trap_kinetics.h
//    - Description: DoubleTrap Kinetics implements the double trap model proposed by
//      Wang et al. and re-developed by Moore et al.
//    - Developers: M. Moore, M. Bhaiya, M. Secanell, P. Wardlaw
//    - $Id: double_trap_kinetics.h 2605 2014-08-15 03:36:44Z secanell $
//
//---------------------------------------------------------------------------

#ifndef _FUELCELLSHOP__DOUBLETRAP_KINETICS_H
#define _FUELCELLSHOP__DOUBLETRAP_KINETICS_H

// Include openFCST routines:
#include <reactions/base_kinetics.h>
#include <utils/fcst_constants.h>

// Include deal.II classes
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
#include <deal.II/base/function.h>
#include <deal.II/lac/vector.h>
#include <deal.II/fe/fe_values.h>

//Include STL
#include <cmath>
#include <iostream>

using namespace dealii;

namespace FuelCellShop
{
    namespace Kinetics
    {
        class DoubleTrapKinetics : 
        public  BaseKinetics
        {
        public:


            static const std::string concrete_name;         
            


            DoubleTrapKinetics(const bool);
            DoubleTrapKinetics();
            ~DoubleTrapKinetics();
            virtual void declare_parameters(ParameterHandler &param) const;
                                    
            virtual void initialize(ParameterHandler &param);
            
        


            virtual void current_density (std::vector<double>& );

            virtual void derivative_current (std::map< VariableNames, std::vector<double> >& );
        


            inline void OH_coverage(std::vector<double>& temp) const
            { 
                temp = theta_OH;
            }
            
            inline void O_coverage(std::vector<double>& temp) const
            { 
                temp = theta_O;
            }
            
            inline void G_DA_rate(std::vector<double>& temp) const
            {
                double C;
                if (reactants_map.at(oxygen_concentration)[0] < 0.0)
                    C = 1.0e-12;
                else
                    C = pow(((reactants_map.at(oxygen_concentration)[0])/c_ref_oxygen),0.5);

                std::vector<double> t;

                for(unsigned int i =0; i < intrinsic_rates[0].size(); i++)
                    t.push_back(C*intrinsic_rates[0][i]);

                temp = t;//free_energies[0];
            }
            
            inline void g_da_rate(std::vector<double>& temp) const
            {
                std::vector<double> t;

                for(unsigned int i =0; i < intrinsic_rates[1].size(); i++)
                    t.push_back(intrinsic_rates[1][i]*theta_O[i]);

                temp = t;//free_energies[1];
            }
            
            inline void G_RA_rate(std::vector<double>& temp) const
            {
                double C;
                if (reactants_map.at(oxygen_concentration)[0] < 0.0)
                    C = 1.0e-12;
                else
                    C = pow(((reactants_map.at(oxygen_concentration)[0])/c_ref_oxygen),0.5);

                double Ch;
                if(reactants_map.count(proton_concentration)){
                    Ch=(reactants_map.at(proton_concentration)[0])/c_ref_protons;
                }
                else
                    Ch = 1.0;


                std::vector<double> t;
                for(unsigned int i =0; i < intrinsic_rates[2].size(); i++)
                    t.push_back(C*Ch*intrinsic_rates[2][i]);

                temp = t;//free_energies[2];
            }
            
            inline void g_ra_rate(std::vector<double>& temp) const
            {
                std::vector<double> t;
                for(unsigned int i =0; i < intrinsic_rates[3].size(); i++)
                    t.push_back(theta_OH[i]*intrinsic_rates[3][i]);//free_energies[3];

                temp = t;
            }
            
            inline void G_RT_rate(std::vector<double>& temp) const
            {

                double Ch;
                if(reactants_map.count(proton_concentration)){
                    Ch=(reactants_map.at(proton_concentration)[0])/c_ref_protons;
                }
                else
                    Ch = 1.0;

                std::vector<double> t;
                for(unsigned int i =0; i < intrinsic_rates[4].size(); i++)
                    t.push_back(theta_O[i]*Ch*intrinsic_rates[4][i]);//free_energies[4];

                temp = t;
            }
            
            inline void g_rt_rate(std::vector<double>& temp) const
            {
                std::vector<double> t;
                for(unsigned int i =0; i < intrinsic_rates[5].size(); i++)
                    t.push_back(theta_OH[i]*intrinsic_rates[5][i]);

                temp = t;//free_energies[5];
            }
            
            inline void G_RD_rate(std::vector<double>& temp) const
            {
                double Ch;
                if(reactants_map.count(proton_concentration)){
                    Ch=(reactants_map.at(proton_concentration)[0])/c_ref_protons;
                }
                else
                    Ch = 1.0;

                std::vector<double> t;
                for(unsigned int i =0; i < intrinsic_rates[6].size(); i++)
                    t.push_back(Ch*theta_OH[i]*intrinsic_rates[6][i]);//free_energies[6];

                temp = t;
            }
            
            inline void g_rd_rate(std::vector<double>& temp) const
            {
                std::vector<double> t;
                for(unsigned int i =0; i < intrinsic_rates[5].size(); i++)
                    t.push_back(intrinsic_rates[7][i]*(1.0 -theta_OH[i] -theta_O[i]));

                temp = t;//free_energies[7];
            }
            
            void compute_energies();
        


            virtual void set_reaction_kinetics(const ReactionNames & name)
            {
                if (name == ORR)
                {
                    name_reaction_kinetics = name;
                }
                else
                {
                    const std::type_info& info = typeid(*this);
                    FcstUtilities::log << "Only ORR reaction is to be implemented in " << __FUNCTION__ << " called in Class " << info.name()  << std::endl;
                    exit(1);
                }
            }


            virtual bool has_coverage(const VariableNames& type){

                if((type == VariableNames::OH_coverage)or (type == VariableNames::O_coverage))
                    return true;

                return false;
            }

        protected: 
        


            virtual boost::shared_ptr<FuelCellShop::Kinetics::BaseKinetics > create_replica ()
            {
                return boost::shared_ptr<FuelCellShop::Kinetics::BaseKinetics > (new FuelCellShop::Kinetics::DoubleTrapKinetics () );
            }
            static DoubleTrapKinetics const* PROTOTYPE;
        
        private:


            inline double negative_concentration_correction(unsigned int index) const
            {
                double correction(1.0);
                
                if (reactants_map.at(oxygen_concentration)[index] <= 0.0)
                    correction = 0.0;
                
                return correction;
            }

            inline double oxygen_concentration_ratio(unsigned int index) const
            {
                double ratio;
                double exponential = 0.5;

                if (reactants_map.at(oxygen_concentration)[index] < 0.0)
                    ratio = 0.0;
                else
                    ratio = pow((reactants_map.at(oxygen_concentration)[index]/c_ref_oxygen), exponential);
                
                return ratio;
            }
            
            inline double proton_concentration_ratio(unsigned int index) const
            {
                double ratio = 1.0;
                
                if(reactants_map.count(proton_concentration))
                    ratio =(reactants_map.at(proton_concentration)[index])/c_ref_protons;
                
                return ratio;
            }
            virtual void init_kin_param()
            {
                Assert( !kin_param_initialized, ExcInternalError() );
                Assert( catalyst != NULL, ExcMessage("Catalyst object not initialized in the DoubleTrapKinetics object.") );
                Assert( catalyst->get_reaction_name() == ORR, ExcMessage("Catalyst object in the DoubleTrapKinetics not set to ORR reaction name.") );
                Assert( phi_m.is_initialized() && phi_s.is_initialized() && T.is_initialized(), ExcMessage("Either phi_m/phi_s/T is not set in the DoubleTrapKinetics object.") );
                Assert( reactants_map.find(oxygen_concentration) != reactants_map.end(), ExcMessage("Oxygen concentration is not set in the DoubleTrapKinetics object.") );

                std::vector<VariableNames> names(1, oxygen_concentration);
                names.push_back(proton_concentration);
                std::map< VariableNames, double > cref_map;
                catalyst->reference_concentration(names, cref_map);
                
                c_ref_oxygen = cref_map[oxygen_concentration];
                c_ref_protons = cref_map[proton_concentration];
                kin_param_initialized = true;
            }
        
            void d_OH_coverage_du (std::map< VariableNames, std::vector<double> >& ) const;
            
            void d_O_coverage_du (std::map< VariableNames, std::vector<double> >& ) const;
            
            void drate_du(std::vector<std::vector<double> >&, const VariableNames& ) const;
                       
            std::vector<std::vector<double> > free_energies;
            
            double G_DA_0, G_RA_0, G_RT_0, G_RD_0;
            
            double G_O, G_OH;
            
            std::vector<std::vector<double> > intrinsic_rates;
            
            std::vector<double> theta_OH, theta_O;
            
            double prefactor;
            
            double c_ref_oxygen;
            
            double c_ref_protons;

            double alpha;
        };        
    } //Kinetics 
} //FuelCellShop

#endif //_FUELCELLSHOP__DOUBLETRAP_KINETICS_H