compressible_multi_component_KG_equations_coupled.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
//
// FCST: Fuel Cell Simulation Toolbox
//
// Copyright (C) 2006-2015 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: compressible_multi_component_KG_equations_coupled.h
// - Description: This class describes steady-state compressible and isothermal Kerkhof-Geboers
//   fluid transport equations for a single-phase multi-component case coupled with fuel cell
//   physics
// - Developers: Valentin N. Zingan, Chad Balen and Marc Secanell, University of Alberta
//
// ----------------------------------------------------------------------------

#ifndef _FCST_FUELCELLSHOP_EQUATION_COMPRESSIBLE_MULTI_COMPONENT_KG_EQUATIONS_COUPLED_H_
#define _FCST_FUELCELLSHOP_EQUATION_COMPRESSIBLE_MULTI_COMPONENT_KG_EQUATIONS_COUPLED_H_

#include <equations/equation_base.h>
#include <layers/channel.h>
#include <layers/gas_diffusion_layer.h>
#include <layers/micro_porous_layer.h>
#include <layers/experimental_porous_layer.h>
#include <layers/catalyst_layer.h>
#include <utils/fcst_units.h>

namespace FuelCellShop
{
    namespace Equation
    {
        
        template<int dim>
        class CompressibleMultiComponentKGEquationsCoupled : public EquationBase<dim>
        {
        public:
            

            
            CompressibleMultiComponentKGEquationsCoupled(FuelCell::SystemManagement& system_management,boost::shared_ptr< FuelCell::ApplicationCore::ApplicationData > data = 
            boost::shared_ptr< FuelCell::ApplicationCore::ApplicationData >());
            
            virtual ~CompressibleMultiComponentKGEquationsCoupled();
            
            virtual void declare_parameters(ParameterHandler& param) const;
            
            virtual void initialize(ParameterHandler& param);
            
            

            
            virtual void assemble_cell_matrix(FuelCell::ApplicationCore::MatrixVector&                                 cell_matrices,
                                              const typename FuelCell::ApplicationCore::DoFApplication<dim>::CellInfo& cell_info,
                                              FuelCellShop::Layer::BaseLayer<dim>* const                               layer);  
            virtual void assemble_cell_residual(FuelCell::ApplicationCore::FEVector&                                     cell_residual,
                                                const typename FuelCell::ApplicationCore::DoFApplication<dim>::CellInfo& cell_info,
                                                FuelCellShop::Layer::BaseLayer<dim>* const                               layer);
            
            virtual void assemble_bdry_matrix(FuelCell::ApplicationCore::MatrixVector&                                 bdry_matrices,
                                              const typename FuelCell::ApplicationCore::DoFApplication<dim>::FaceInfo& bdry_info,
                                              FuelCellShop::Layer::BaseLayer<dim>* const                               layer);
            
            virtual void assemble_bdry_residual(FuelCell::ApplicationCore::FEVector&                                     bdry_residual,
                                                const typename FuelCell::ApplicationCore::DoFApplication<dim>::FaceInfo& bdry_info,
                                                FuelCellShop::Layer::BaseLayer<dim>* const                               layer);
            
            


            const std::vector<double>& get_inlet_outlet_velocity_max() const
            {
                return inlet_outlet_velocity_max;
            }
            
            const double get_inlet_outlet_velocity_mixture_max() const
            {
                return inlet_outlet_velocity_mixture_max;
            }
            
            const bool get_use_parabolic_profile() const
            {
                return use_parabolic_profile;
            }
            
            const double get_inlet_outlet_velocity_channel_base() const
            {
                return inlet_outlet_velocity_channel_base;
            }
            
            const double get_inlet_outlet_velocity_channel_roof() const
            {
                return inlet_outlet_velocity_channel_roof;
            }
            
            const int get_inlet_outlet_boundary_ID() const
            {
                return inlet_outlet_boundary_ID;
            }
            
            const std::string get_press_vel_comp_apply_to() const
            {
                return press_vel_comp_apply_to;
            }
            
            const std::string get_inlet_outlet_velocity_mixture_equation() const
            {
                return inlet_outlet_velocity_mixture_equation;
            }
            
            const std::map<unsigned int, std::string> get_gas_species_map() const
            {
                return gas_species_map;
            }
            
            const unsigned int get_num_of_species() const
            {
                return n_species;
            }
            
            virtual void print_equation_info() const;
            
            
        protected:
            

            
            template<typename INFO>
            void make_assemblers_generic_constant_data(const INFO&                                InFo,
                                                       FuelCellShop::Layer::BaseLayer<dim>* const layer);
            
            virtual void make_assemblers_cell_constant_data(const typename FuelCell::ApplicationCore::DoFApplication<dim>::CellInfo& cell_info);
            
            virtual void make_assemblers_bdry_constant_data(const typename FuelCell::ApplicationCore::DoFApplication<dim>::FaceInfo& bdry_info);
            
            virtual void make_assemblers_cell_variable_data(const typename FuelCell::ApplicationCore::DoFApplication<dim>::CellInfo& cell_info,
                                                            FuelCellShop::Layer::BaseLayer<dim>* const                               layer);
            
            virtual void make_assemblers_bdry_variable_data(const typename FuelCell::ApplicationCore::DoFApplication<dim>::FaceInfo& bdry_info,
                                                            FuelCellShop::Layer::BaseLayer<dim>* const                               layer);
            
            

            
            virtual void make_internal_cell_couplings();
            
            virtual void make_matrix_block_indices();
            
            virtual void make_residual_indices();
            
            
            // DATA //
            

            
            std::vector<bool> inertia_in_channels;
            
            std::vector<bool> shear_stress_in_channels;
            
            std::vector<bool> gravity_in_channels;
            
            std::vector<bool> inertia_in_porous_media;
            
            std::vector<bool> shear_stress_in_porous_media;
            
            std::vector<bool> gravity_in_porous_media;
            
            std::vector<std::string> drag_in_porous_media;
            
            std::vector<bool> normal_velocity_is_suppressed_weakly;
            
            

            
            unsigned int n_species;
            
            bool coupled_with_fuel_cell_physics;
            
            std::map<unsigned int, std::string> gas_species_map;
            
            double R_universal;
            
            double T_mixture;
            
            double P_in;
            
            Tensor<1,dim> gravity_acceleration;
            
            SymmetricTensor<2,dim> unit;
            
            std::vector<double> molar_mass;
            
            std::vector<double> dynamic_viscosity;
            
            std::vector<double> bulk_viscosity;
            
            std::vector<double> collision_diameter;
            
            std::vector<double> theta;
            
            std::vector<double> eta;
            
            Table< 2, double > maxwell_stefan_isobaric_diffusion_coefficient;
            
            std::vector<double> inlet_outlet_velocity_max;
            
            double      inlet_outlet_velocity_mixture_max;
            bool        use_parabolic_profile;
            double      inlet_outlet_velocity_channel_base;
            double      inlet_outlet_velocity_channel_roof;
            int         inlet_outlet_boundary_ID;
            std::string press_vel_comp_apply_to;
            std::string inlet_outlet_velocity_mixture_equation;
            
            std::vector<double> maxwell_constant;
            
            

            
            std::vector< FEValuesExtractors::Scalar > density;
            
            std::vector< FEValuesExtractors::Vector > velocity;
            
            

            
            

            
            

            
            

            
            std::vector< std::vector<double> > density_cell_old;
            
            std::vector< std::vector< Tensor<1,dim> > > grad_density_cell_old;
            
            std::vector< std::vector< Tensor<1,dim> > > velocity_cell_old;
            
            std::vector< std::vector<double> > div_velocity_cell_old;
            
            std::vector< std::vector< SymmetricTensor<2,dim> > > grads_velocity_cell_old;
            
            std::vector< std::vector< Tensor<1,dim> > > mass_flux_cell_old;
            
            std::vector< std::vector< SymmetricTensor<2,dim> > > momentum_flux_cell_old;
            
            std::vector< std::vector<double> > pressure_cell_old;
            
            std::vector< std::vector< SymmetricTensor<2,dim> > > shear_stress_cell_old;
            
            std::vector< std::vector<double> > partial_viscosity_old;
            
            std::vector< std::vector<double> > bulk_viscosity_old;
            
            std::vector< std::vector< std::vector<double> > > paramMatrix_old;
            
            std::vector< FullMatrix<double> > PInv_old;
            
            std::vector< std::vector< Tensor<1,dim> > > drag_cell_old;
            
            std::vector< std::vector< Tensor<1,dim> > > diffusion_cell_old;
            
            std::vector< std::vector< Tensor<1,dim> > > gravity_cell_old;
            
            

            
            std::vector< std::vector< std::vector<double> > > phi_density_cell;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > grad_phi_density_cell;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > phi_velocity_cell;
            
            std::vector< std::vector< std::vector<double> > > div_phi_velocity_cell;
            
            std::vector< std::vector< std::vector< SymmetricTensor<2,dim> > > > grads_phi_velocity_cell;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > delta_mass_flux_cell;
            
            std::vector< std::vector< std::vector< SymmetricTensor<2,dim> > > > delta_momentum_flux_cell;
            
            std::vector< std::vector< std::vector<double> > > delta_pressure_cell;
            
            std::vector< std::vector< std::vector<double> > > delta_partial_viscosity_cell;
            
            std::vector< std::vector< std::vector<double> > > delta_bulk_viscosity_cell;
            
            std::vector< std::vector< std::vector< SymmetricTensor<2,dim> > > > delta_shear_stress_cell;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > delta_drag_cell;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > delta_diffusion_cell;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > delta_gravity_cell;
            
            

            
            

            
            std::vector< std::vector<double> > density_bdry_old;
            
            std::vector< std::vector< Tensor<1,dim> > > grad_density_bdry_old;
            
            std::vector< std::vector< Tensor<1,dim> > > velocity_bdry_old;
            
            std::vector< std::vector<double> > div_velocity_bdry_old;
            
            std::vector< std::vector< SymmetricTensor<2,dim> > > grads_velocity_bdry_old;
            
            std::vector< std::vector< Tensor<1,dim> > > mass_flux_bdry_old;
            
            std::vector< std::vector< SymmetricTensor<2,dim> > > momentum_flux_bdry_old;
            
            std::vector< std::vector<double> > pressure_bdry_old;
            
            std::vector< std::vector< SymmetricTensor<2,dim> > > shear_stress_bdry_old;
            
            std::vector< std::vector<double> > partial_viscosity_bdry_old;
            
            std::vector< std::vector<double> > bulk_viscosity_bdry_old;
            
            std::vector< std::vector< std::vector<double> > > paramMatrix_bdry_old;
            
            std::vector< FullMatrix<double> > PInv_bdry_old;
            
            

            
            std::vector< std::vector< std::vector<double> > > phi_density_bdry;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > grad_phi_density_bdry;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > phi_velocity_bdry;
            
            std::vector< std::vector< std::vector<double> > > div_phi_velocity_bdry;
            
            std::vector< std::vector< std::vector< SymmetricTensor<2,dim> > > > grads_phi_velocity_bdry;
            
            std::vector< std::vector< std::vector< Tensor<1,dim> > > > delta_mass_flux_bdry;
            
            std::vector< std::vector< std::vector< SymmetricTensor<2,dim> > > > delta_momentum_flux_bdry;
            
            std::vector< std::vector< std::vector<double> > > delta_pressure_bdry;
            
            std::vector< std::vector< std::vector< SymmetricTensor<2,dim> > > > delta_shear_stress_bdry;
            
            std::vector< std::vector< std::vector<double> > > delta_partial_viscosity_bdry;
            
            std::vector< std::vector< std::vector<double> > > delta_bulk_viscosity_bdry;
            
            

            
            std::vector< std::vector< std::vector< SymmetricTensor<2,dim> > > > n_BY_phi_velocity_bdry_S;
            


            template <typename LAYER, typename INFO>
            inline void update_porosity(const LAYER* ptr, 
                                        const INFO info, 
                                        std::vector<double>& porosity)
            {
                if( ptr->get_porosity_is_constant() ) 
                    ptr->get_porosity(porosity);
                else
                    ptr->get_porosity(porosity, info.get_fe_val_unsplit().get_quadrature_points());
            }
            
            template <typename LAYER, typename INFO>
            inline void update_permeability(const LAYER* ptr, 
                                            const INFO info, 
                                            std::vector< SymmetricTensor<2,dim> >& permeability_INV,
                                            std::vector< SymmetricTensor<2,dim> >& Forchheimer_permeability)
            {
                if( ptr->get_permeability_is_constant() )
                {
                    ptr->get_permeability_INV(permeability_INV);
                    ptr->get_Forchheimer_permeability(Forchheimer_permeability);
                }
                else
                {
                    ptr->get_permeability_INV(permeability_INV, info.get_fe_val_unsplit().get_quadrature_points());
                    ptr->get_Forchheimer_permeability(Forchheimer_permeability, info.get_fe_val_unsplit().get_quadrature_points());
                }
            }
            template <typename LAYER, typename INFO>
            inline void update_tortuosity(const LAYER* ptr, 
                                          const INFO info, 
                                          std::vector< SymmetricTensor<2,dim> >& tortuosity )
            {
                if( ptr->get_tortuosity_is_constant() )
                    ptr->get_tortuosity(tortuosity);
                else
                    ptr->get_tortuosity(tortuosity, info.get_fe_val_unsplit().get_quadrature_points());
            }
            
            template <typename LAYER>
            inline void update_gas_properties(const LAYER* ptr)
            {
                //TODO: these parameters should be coming from the same location not different classes
                T_mixture = ptr->get_gas_mixture()->get_temperature();
                P_in = ptr->get_gas_mixture()->get_total_pressure();
                
                for(unsigned int s = 0; s < n_species; ++s)
                {
                    molar_mass[s] = ptr->get_gas_mixture()->get_gases()[s]->get_molar_mass() * 1.0e3; // Convert kg/mol to 1000 g/mol
                    dynamic_viscosity[s] = ptr->get_gas_mixture()->get_gases()[s]->get_dynamic_viscosity(T_mixture) * 1.0e1; // Convert kg/m*s to 10 g/cm*s
                    collision_diameter[s] = ptr->get_gas_mixture()->get_gases()[s]->get_collision_diameter() * 1e-10; // Convert Angstroms to 1e-10 meters
                    maxwell_stefan_isobaric_diffusion_coefficient = ptr->get_gas_mixture()->get_ChapmanEnskog_isobaric_diffusion_coefficients(T_mixture); // [m^2/s]
                }
            }
            
            template <typename LAYER>
            inline void update_inv_diffusion_coefficients(const LAYER* ptr, 
                                                          std::vector< Table< 2, SymmetricTensor<2,dim> > >& inv_pDeff) const
            {
                //************************************
                // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                // !!!!!!!!!!!!!! THIS NEEDS TO BE IMPROVED !!!!!!!!!!!!!!
                // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                //************************************   
                Table< 2, Tensor<2,dim> > Deff;
                
                //I would like to replace this...
                ptr->effective_gas_diffusivity(Deff); // m^2/s
                // with
                /*
                ptr->get_T_and_p(T,p);
                SolutionVariable temperature = SolutionVariable(T, this->n_q_points_cell , temperature_of_REV);
                ptr->set_temperature( temperature );
                SolutionVariable pressure = SolutionVariable(1.0/101325.0, this->n_q_points_cell, temperature_of_REV); // Value of 1.0/101325.0 is equivalent to 1 Pa, i.e. we return D(T)
                const std::vector< PureGas *> gases = ptr->get_gases ();
                */
                /* If saturation:
                if (p_liquid_water.indices_exist)
                {
                    ptr->set_capillary_pressure( FuelCellShop::SolutionVariable(&cell_info.values[last_iter_cell][p_liquid_water.solution_index], capillary_pressure) );
                    deriv_flags.push_back(capillary_pressure);
                }
                */              
                
                unsigned int n_q_points_cell = inv_pDeff.size();
                
                for(unsigned int s = 0; s < n_species; ++s)
                    for(unsigned int s1 = 0; s1 < n_species; ++s1)
                    {
                        // std::vector< double > Deff_noniso;
                        // ptr->compute_gas_diffusion(gases[s], gases[s1]);
                        // ptr->effective_gas_diffusivity(Deff_noniso);            // m^2/s
                        
                        if( s1 != s )  
                        {
                            for(unsigned int q = 0; q < this->n_q_points_cell; ++q)   
                            {                                                               
                                Tensor<2, dim> pDeff_CGS = P_in*Deff(s,s1)*Units::convert(1.,Units::C_UNIT2, Units::UNIT2);
                                // Units are Pa*cm^2/s, however Pa are NOT CGS, ( 1 Pa = 10 g/(cm s^2):
                                pDeff_CGS *= 10;
                                // Now, compute 1/(p*Deff):                                    
                                Tensor<2, dim> invTest;
                                invTest[0][0] = 1 / ( pDeff_CGS[0][0] );
                                invTest[1][1] = 1 / ( pDeff_CGS[1][1] );
                                #if deal_II_dimension == 3
                                    invTest[2][2] = 1 / ( pDeff_CGS[2][2] );
                                #endif
                                inv_pDeff[q](s,s1) = invTest;
                            }
                        }
                    }                    

                    //Print out Deff* for debugging
//                     bool printDeff = false;
//                     
//                     if( printDeff )
//                     {
//                         FcstUtilities::log << "-----------------------------------" <<std::endl;
//                         for( unsigned int s1 = 0; s1 < n_species; ++s1 )
//                         {
//                             for( unsigned int s2 = 0; s2 < n_species; ++s2 )
//                             {
//                                 if( s1 != s2 )
//                                 {
//                                     FcstUtilities::log << "Name layer is: "<<ptr->name_layer()<<std::endl;
//                                     FcstUtilities::log << "Species " << s1 << " and " << s2 << ": " <<std::endl;
//                                     FcstUtilities::log << "Deff[0][0]: " << Deff(s1,s2)[0][0]*Units::convert(1.,Units::C_UNIT2, Units::UNIT2) << std::endl;
//                                     FcstUtilities::log << "Deff[1][1]: " << Deff(s1,s2)[1][1]*Units::convert(1.,Units::C_UNIT2, Units::UNIT2) << std::endl;
//                                     
//                                     FcstUtilities::log << "inv_pDeff[0][0]: " << inv_pDeff[0](s1,s2)[0][0] << std::endl;
//                                     FcstUtilities::log << "inv_pDeff[1][1]: " << inv_pDeff[0](s1,s2)[1][1] << std::endl;
//                                 }
//                             }
//                         }
//                         FcstUtilities::log << "-----------------------------------" <<std::endl;
//                     }
            }
            
            template <typename LAYER>
            inline void update_partial_viscosities(const LAYER*                                       ptr,
                                                   const unsigned int&                                quadraturePoints,
                                                   const std::vector<double>&                         porosity,
                                                   const std::vector< std::vector<double> >&          density,
                                                   std::vector< std::vector< std::vector<double> > >& paramMatrix,
                                                   std::vector< FullMatrix<double> >&                 PInv,
                                                   std::vector< std::vector<double> >&                partialViscosity,
                                                   std::vector< std::vector<double> >&                bulkViscosity) const
            {

                partialViscosity = ptr->get_gas_mixture()->get_isothermal_nonisobaric_partial_viscosity(T_mixture,          // [K]
                                                                                                        density,            // [g/cm^3]
                                                                                                        molar_mass,         // [g/mol]
                                                                                                        dynamic_viscosity,  // [g/cm*s]
                                                                                                        collision_diameter, // [m]
                                                                                                        porosity, 
                                                                                                        paramMatrix,
                                                                                                        PInv); // returns in [g/cm*s]
                
                bulkViscosity.clear(); //reset bulk viscosity vector
                bulkViscosity.resize(quadraturePoints, std::vector<double>(this->n_species, 0.0)); //size appropriately
                for(unsigned int q = 0; q < quadraturePoints; ++q)
                {
                    double porosityInv = 1 / porosity[q];
                    for(unsigned int s = 0; s < this->n_species; ++s)
                        bulkViscosity[q][s] = ptr->get_gas_mixture()->get_gases()[s]->get_bulk_viscosity(partialViscosity[q][s]); // [g/cm*s]
                }
            }
            
            template <typename LAYER>
            inline void update_delta_partial_viscosities(const LAYER*                                             ptr,
                                                         const std::vector< std::vector<double> >&                partialViscosity,
                                                         const std::vector< std::vector< std::vector<double> > >& paramMatrix,
                                                               std::vector< FullMatrix<double> >&                 PInv,
                                                         const std::vector<double>&                               porosity,
                                                         const std::vector< std::vector< std::vector<double> > >& deltaDensity,
                                                         const std::vector< std::vector<double> >&                density,
                                                               std::vector< std::vector< std::vector<double> > >& deltaPartialViscosity) const
            {
                ptr->get_gas_mixture()->get_isothermal_nonisobaric_delta_partial_viscosity(T_mixture,          // [K]
                                                                                           paramMatrix,
                                                                                           PInv,
                                                                                           porosity,
                                                                                           molar_mass,         // [g/mol]
                                                                                           dynamic_viscosity,  // [g/cm*s]
                                                                                           collision_diameter, // [m]
                                                                                           deltaDensity,
                                                                                           density,            // [g/cm^3]
                                                                                           deltaPartialViscosity); // returns in cgs units
            }
            
            template <typename LAYER>
            inline void update_delta_bulk_viscosities(const LAYER*                                             ptr,
                                                      const std::vector< std::vector< std::vector<double> > >& deltaPartialViscosity,
                                                            std::vector< std::vector< std::vector<double> > >& deltaBulkViscosity) const
            {
                ptr->get_gas_mixture()->get_delta_bulk_viscosity(ptr->get_gas_mixture()->get_gases(), 
                                                                 deltaPartialViscosity, 
                                                                 deltaBulkViscosity);
            }
            
        };
                 
    } // Equation
                 
} // FuelCellShop
             
#endif