3.0 liquid pool

include/cantera/oneD/Boundary1D.h

@@ -14,6 +14,7 @@
 #include "cantera/thermo/SurfPhase.h"
 #include "cantera/kinetics/InterfaceKinetics.h"
 #include "StFlow.h"
+#include "Phase_liquid.h"
 
 namespace Cantera
 {
@@ -389,6 +390,43 @@ class ReactingSurf1D : public Boundary1D
     vector<double> m_fixed_cov;
 };
 
+/**
+ * This is a boundary condition from inlet 1D
+ * I want to change from gas boundary to liquid pool boundary:
+ * (Temperature, mass_fraction, lambda)
+ */
+
+class Inlet1D_new : public Inlet1D{
+public:
+
+    void addLiquidBcs(Phase_liquid &m){
+        pool=&m;
+    }
+
+    void setMdot(double mdot) {
+        m_mdot = mdot;
+    }
+    void showMdot() {
+        std::cout<<"Mass flow rate is "<<m_mdot<<std::endl;
+    }
+
+    void setMoleFractions(const std::string& xin);
+    void setMoleFractions(const double* xin);
+
+    double massFraction(size_t k) {
+        return m_yin[k];
+    }
+
+    void init();
+    void eval(size_t jg, double* xg, double* rg,
+                      integer* diagg, double rdt);
+
+    Phase_liquid* Pool(){return pool;}
+protected:
+    Phase_liquid* pool;
+};
+
+
 //! @} End of bdryGroup
 
 }

include/cantera/oneD/Data_post_analysis.h

@@ -0,0 +1,89 @@
+//
+// Created by Liang on 2023/4/21.
+//
+
+#ifndef BCS_DATA_POST_ANALYSIS_H
+#define BCS_DATA_POST_ANALYSIS_H
+#include "Sim1D.h"
+#include "StFlow.h"
+#include "Boundary1D.h"
+#include <stdio.h>
+#include <fstream>
+#include "cantera/base/Solution.h"
+#include <filesystem>
+
+namespace Cantera
+{/*this file is designed to do sensitivity analysis for liquid-pool counterflow flame results*/
+    class Data_post_analysis {
+    public:
+        Data_post_analysis(){
+            std::cout<<"\n..............................................................................\n";
+            std::cout<<"*************************Data post analysis is called*************************\n";
+            std::cout<<"..............................................................................\n";}
+        //For gas phase counterflow flame analysis
+        void Initialize(Sim1D &flame, StFlow &flow, shared_ptr<Solution> sol, std::string sol_id);
+        //For liquid pool conterflow flame analysis
+        void Initialize(Sim1D &flame, StFlow &flow,Inlet1D_new &inlet, Inlet1D_new &outlet, shared_ptr<Solution> sol, std::string sol_id);
+
+        //Restore flame with existed solution
+        void Restore(const std::string &filename, const std::string &file_id )
+        {m_flame->restore(filename, file_id);
+            std::cout<<m_flame->nDomains()<<std::endl;}
+
+        //print temperature for all points
+        void Temperature(){for(size_t n=0; n<points;n++){std::cout<<solution_vector[n*n_variable+c_offset_T]<<std::endl;}}
+
+        //return temeprature at point j
+        double Temperature(size_t j){return solution_vector[j*n_variable+c_offset_T];}
+
+        //i starts from the first specie, j is the loc_point
+        double MassFraction(size_t i, size_t j){return solution_vector[i+j*n_variable+c_offset_Y];}
+
+        //Output solution vector and molefraction to file
+        void Output_Solution(std::string Result_file_name);
+
+        //Output solution at liquid-gas phase
+        void Output_liquid_pool(std::string Liquid_pool_solution_name);
+
+        //Output reaction rate and constant to file
+        void Output_Reaction_rate(std::string Reaction_rate_file_name);
+
+        //Output heat_prod_rate to file
+        void Output_Heat_Prod_Rate(std::string Heat_Prod_Rate_filename);
+
+        //calculate sensitivity
+        void get_temperature_reaction_sensitivity(size_t point_loc);
+        void get_oxygen_sensitivity(size_t point_loc);
+        void get_sensitivity(size_t point_loc, std::string variable);
+        void get_sensitivity_all(std::string variable);
+
+        void Output_Sensitivity(std::vector<std::string> var_vec, std::string Sensitivity_file_name);
+
+        void perturb(size_t i, double dp ){m_flow->kinetics().setMultiplier(i,1+dp);}
+
+        vector<double> solve_adjoint( size_t n_params, vector<double> dgdx, double(*g_ptr)()=nullptr, double dp=1e-5 );
+
+        void solve_adjoint_all(size_t n_params, Array2D* dgdx_all, double(*g_ptr)()=nullptr, double dp=1e-5 );
+
+        void check_path(std::string path_name){
+            std::filesystem::path folderPath(path_name);
+            if (std::filesystem::exists(folderPath)) {
+            } else {
+                std::filesystem::create_directories(folderPath);
+            }
+        }
+    protected:
+        Sim1D *m_flame;
+        StFlow *m_flow;
+        Inlet1D_new *m_inlet, *m_outlet;
+        Phase_liquid *m_pool;
+        shared_ptr<Solution> m_sol;
+        vector<double> r_Mean_mw,solution_vector,Sens_reaction;
+        std::vector<std::string> reaction_name_vector;
+        size_t nsp, points, n_variable, n_reaction;
+        Array2D dfdp, Sens_reaction_all;
+        bool dfdp_flag;
+        std::string folder_name;
+    };
+}
+#endif //BCS_DATA_POST_ANALYSIS_H

include/cantera/oneD/Phase_liquid.h

@@ -0,0 +1,133 @@
+//
+// Created by Liang on 2023/1/8.
+//
+
+#ifndef BCS_PHASE_LIQUID_H
+#define BCS_PHASE_LIQUID_H
+#include "cantera/thermo/Phase.h"
+#include "StFlow.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <string>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+namespace Cantera
+{/**The class of boundary for liquid pool used in one-dimensional spatial domains.
+ * setup inlet molar fraction of liquid pool and get outlet conditions
+*/
+class Phase_liquid
+{
+public:
+    //!constructor
+    Phase_liquid();
+
+    //![K] set inlet temperature for liquid_phase, default temperature is 298.15K
+    void setTemperature_inlet(double t){temperature_inlet=t;};
+
+    //!return inlet temeprature for liquid_pool
+    double Temperature_inlet(){return temperature_inlet;};
+
+    //![K] set outlet temperature for liquid_phase, default temperature is 298.15K
+    void setTemperature_outlet(double t){temperature_outlet=t;
+       update();};
+
+    //![K] return inlet temperature for liquid_phase
+    double Temperature_outlet(){return temperature_outlet;};
+
+    //!set inlet mole fraction of liquid pool
+    void setMoleFraction(std::string name, const double mole);
+
+    //!set inlet mole fraction of liquid pool
+    void setMoleFraction(const std::string& fuel_composition);
+
+    //!calculate conditions (temperature, molar fraction) at liquid side of interface based on gas phase (support for mixture less than two components)
+    void Build_interface(double *xb, Cantera::StFlow *m_flow );
+
+    //!return inlet mole fraction by number
+    double MoleFraction_inlet(int m){return mole_fraction_inlet[m];}
+    //!return inlet mole fraction by name
+    double MoleFraction_inlet(std::string species_name){return MoleFraction_inlet(speciesIndex(species_name));}
+
+    //!return outlet mass fraction by number
+    double MassFraction_inlet(int m) {return mole_fraction_inlet[m]*molecular_weight[m]/sum_mole_weight_inlet;}
+    //!return outlet mass fraction by name
+    double MassFraction_inlet(std::string species_name){return MassFraction_inlet(speciesIndex(species_name));}
+
+    //!return outlet mole fraction by number
+    double MoleFraction_outlet(int m){return mole_fraction_outlet[m];}
+    //!return outlet mole fraction by name
+    double MoleFraction_outlet(std::string species_name){ return MoleFraction_outlet(speciesIndex(species_name));}
+
+    //!return outlet mass fraction by number
+    double MassFraction_outlet(int m){return mole_fraction_outlet[m]*molecular_weight[m]/sum_mole_weight_outlet;};
+    //!return outlet mass fraction by name
+    double MassFraction_outlet(std::string species_name){ return MassFraction_outlet(speciesIndex(species_name));};
+
+    //!update sum_molecular weight and sum_liquid enthalpy and temperature
+    void update();
+
+    //![J/kmol] liquid enthalpy of each single component at specified temperature
+    double Liquid_enthpy(std::string name, double temperature);
+
+    //![J/kmol] liquid enthalpy of mixture
+    double Liquid_enthpy(){return sum_liquid_ethalpy;}
+
+    //![J/kmol] enthalpy of component in gas phase
+    double enthpy_gas(std::string name, StFlow *m_flow);
+
+    //![J/kmol] enthalpy change between liquid and gas in molar basis
+    double Delta_H_mole(doublereal *xb, StFlow *m_flow);
+    double Delta_H_mole(double t);
+
+    //![J/kg] enthalpy change at outlet from liquid to gas in mass basis
+    double Delta_H_mass(doublereal *xb,StFlow *m_flow){return Delta_H_mole(xb,m_flow)/sum_mole_weight_outlet;}
+    double Delta_H_mass(double t){return Delta_H_mole(t)/sum_mole_weight_outlet;}
+
+    //![K] return temerpature at the liquid-gas interface
+    double Temperature_interface(double *xb, StFlow *m_flow);
+
+    //![J/kmol] vaporization enthalpy of each single component at specified temperature
+    double Vap_enthpy(std::string name, double temperature);
+
+    //![Pa] return of saturated pressure of specie (name) at temperature [K]
+    double SatPressure(std::string name, double temperature);
+
+    //![K] return of saturated temperature of specie (name) at pressure [Pa]
+    double SatTemperature(std::string name, double pressure);
+
+    //!return the number of species in liquid pool
+    size_t Nsp_liquid(){return nsp_liquid;}
+
+    //!return location of species at liquid phase
+    size_t speciesIndex(const std::string nameStr){
+        size_t loc=npos;
+        auto it=m_speciesIndices.find(nameStr);
+        if (it!=m_speciesIndices.end()){return it->second;}
+        return loc;
+    }
+    //!return name of species at liquid phase
+    std::string speciesName(const size_t n){return m_speciesName.find(n)->second;}
+
+protected:
+    vector<double> mole_fraction_inlet, mole_fraction_outlet, mole_fraction_gas;
+    vector<double> liquid_enthalpy,delta_H_sp;//[J/kmol]
+    std::vector<size_t> Loc_sp;//location of activated species in pool
+
+    std::vector<std::string> name{"C2H5OH", "NC7H16", "C7H16", "H2O"};
+    vector<double>    H_formation{-276*1e6, -224*1e6, -224*1e6, -285.83*1e6};//J/kmol
+    vector<double>             Cp{112*1e3, 224.64*1e3, 224.64*1e3, 75.37*1e3};//J/kmol/K}
+    vector<double> molecular_weight{46.07, 100.21, 100.21, 18.015};//kg/kmol
+
+    //!map of species names to indices
+    std::map<std::string, size_t> m_speciesIndices;
+    //!map of species indices to names
+    std::map<size_t, std::string> m_speciesName;
+
+    double sum_mole_weight_inlet, sum_mole_weight_outlet, sum_liquid_ethalpy,temperature_inlet, temperature_outlet, temperature_interface_gas;
+    double delta_H_sum;
+
+    size_t nsp_liquid ;//number of liquid species activated
+};
+}
+#endif //BCS_PHASE_LIQUID_H

include/cantera/oneD/Sim1D.h

@@ -9,6 +9,7 @@
 #define CT_SIM1D_H
 
 #include "OneDim.h"
+#include "cantera/base/Array.h"
 
 namespace Cantera
 {
@@ -332,6 +333,8 @@ class Sim1D : public OneDim
      * @f]
      */
     void solveAdjoint(const double* b, double* lambda);
+
+    void solveAdjoint_all(const Array2D* b_all, Array2D* lambda_all);
 
     void resize() override;
 

include/cantera/oneD/StFlow.h

@@ -345,6 +345,70 @@ class StFlow : public Domain1D
     size_t rightExcessSpecies() const {
         return m_kExcessRight;
     }
+    //! return heat conductivity vector (added for boundary calculation) Unit W/m-K
+    vector<double> heat_conductivity()
+    {return m_tcon;}
+
+    //! return grid vector (added for boundary calculation)
+    vector<double> grid_dz()
+    {return m_dz;}
+
+    //! return enthalpy vector (added for boundary calculation) Unit J/kmol
+    vector<double> Enthalpy()
+    {vector<double> Enthalpy_inlet(m_enthp.nRows());
+        double* itr= &Enthalpy_inlet[0];
+        m_enthp.getColumn(0, itr );
+    return Enthalpy_inlet;}
+
+    //!return heat capacity (added for boundary calculation)
+    vector<double> Cp()
+    {return m_cp;}
+
+    //!return molecular weight  (added for boundary calculation)
+    doublereal Mw(size_t k)
+    {return m_wt[k];}
+
+    //!return mean molecular weight  (added for boundary calculation)
+    doublereal Mean_Mw(size_t j)
+    {return m_wtm[j];}
+
+    //!return dT/dz (added for boundary calculation)
+    doublereal dT_dz(const doublereal* x, size_t j){
+        return dTdz(x, j);
+    }
+
+    //!return dY/dz (added for boundary calculation)
+    doublereal dY_dz(const doublereal* x, size_t k , size_t j){
+        return dYdz(x, k, j);
+    }
+
+    //!diffusive_flux (added for boundary calculation)
+    doublereal diffusive_flux(size_t k, size_t j)const{
+        return m_flux(k,j);//number of species: k, number of points: j;
+    }
+    //!rho_u (added for boundary calculation)
+    doublereal rho_u(doublereal* x, size_t j)const{
+        return m_rho[j]*x[index(c_offset_U,j)];//number of species: k, number of points: j;
+    }
+    //!molar fraction of individual component (added for boundary calculation)
+    doublereal molar_fraction(const doublereal* x, size_t k, size_t j) const {
+        return m_wtm[j]*Y(x,k,j)/m_wt[k];
+    }
+
+    /*doublereal& Y(doublereal* x, size_t k, size_t j) {
+        return x[index(c_offset_Y + k, j)];
+    }*/
+
+    double Eneg_terms(size_t k, size_t j){
+    std::vector<vector<double>> Eneg_terms={Eneg_RHS_term_1, Eneg_RHS_term_2, Eneg_RHS_term_3, Eneg_RHS_term_4};
+    return Eneg_terms[k][j];
+    }
+
+    // k: term index; j: point index, n: species index
+    double Species_terms(size_t k, size_t j, size_t n){
+    std::vector<Array2D> Species_terms={Sp_RHS_term_1, Sp_RHS_term_2, Sp_RHS_term_3};
+    return Species_terms[k].value(j, n);
+    }
 
 protected:
     AnyMap getMeta() const override;
@@ -381,6 +445,7 @@ class StFlow : public Domain1D
             m_wtm[j] = m_thermo->meanMolecularWeight();
             m_cp[j] = m_thermo->cp_mass();
             m_thermo->getPartialMolarEnthalpies(&m_hk(0, j));
+            m_thermo->getPartialMolarEnthalpies(&m_enthp(0,j));//enthalpy update was added for boundary condition
         }
     }
 
@@ -495,6 +560,7 @@ class StFlow : public Domain1D
     // species thermo properties
     vector<double> m_wt;
     vector<double> m_cp;
+    Array2D m_enthp;//added for boundary condition calculation
 
     // transport properties
     vector<double> m_visc;
@@ -558,6 +624,12 @@ class StFlow : public Domain1D
     //! to `j1`, based on solution `x`.
     virtual void updateTransport(double* x, size_t j0, size_t j1);
 
+    //output of all terms in energy equation
+    vector<double> Eneg_RHS_term_1, Eneg_RHS_term_2, Eneg_RHS_term_3, Eneg_RHS_term_4;
+
+    //output of all iterms in each specie equation
+    Array2D Sp_RHS_term_1, Sp_RHS_term_2, Sp_RHS_term_3;
+
 public:
     //! Location of the point where temperature is fixed
     double m_zfixed = Undef;