canteraTransport working

src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModels.C

@@ -54,6 +54,13 @@ namespace Foam
         gasHThermoPhysics
     );
 
+    makeChemistryModel
+    (
+        chemistryModel,
+        psiChemistryModel,
+        canteraGasHThermoPhysics
+    );
+
     makeChemistryModel
     (
         chemistryModel,

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModels.C

@@ -54,6 +54,13 @@ namespace Foam
         gasHThermoPhysics
     );
 
+    makeChemistryModel
+    (
+        chemistryModel,
+        rhoChemistryModel,
+        canteraGasHThermoPhysics
+    );
+
     makeChemistryModel
     (
         chemistryModel,

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C

@@ -36,6 +36,7 @@ namespace Foam
     // Chemistry solvers based on sensibleEnthalpy
     makeChemistrySolverTypes(psiChemistryModel, constGasHThermoPhysics);
     makeChemistrySolverTypes(psiChemistryModel, gasHThermoPhysics);
+    makeChemistrySolverTypes(psiChemistryModel, canteraGasHThermoPhysics);
     makeChemistrySolverTypes
     (
         psiChemistryModel,
@@ -49,6 +50,7 @@ namespace Foam
     makeChemistrySolverTypes(psiChemistryModel, icoPoly8HThermoPhysics);
     makeChemistrySolverTypes(rhoChemistryModel, constGasHThermoPhysics);
     makeChemistrySolverTypes(rhoChemistryModel, gasHThermoPhysics);
+    makeChemistrySolverTypes(rhoChemistryModel, canteraGasHThermoPhysics);
     makeChemistrySolverTypes
     (
         rhoChemistryModel,

src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C

@@ -41,12 +41,14 @@ namespace Foam
 
 makeChemistryReader(constGasHThermoPhysics);
 makeChemistryReader(gasHThermoPhysics);
+makeChemistryReader(canteraGasHThermoPhysics);
 makeChemistryReader(constIncompressibleGasHThermoPhysics);
 makeChemistryReader(incompressibleGasHThermoPhysics);
 makeChemistryReader(icoPoly8HThermoPhysics);
 
 makeChemistryReaderType(foamChemistryReader, constGasHThermoPhysics);
 makeChemistryReaderType(foamChemistryReader, gasHThermoPhysics);
+makeChemistryReaderType(foamChemistryReader, canteraGasHThermoPhysics);
 makeChemistryReaderType
 (
     foamChemistryReader,

src/thermophysicalModels/reactionThermo/psiReactionThermo/psiReactionThermos.C

@@ -241,6 +241,15 @@ makeReactionMixtureThermo
     gasHThermoPhysics
 );
 
+makeReactionMixtureThermo
+(
+    psiThermo,
+    psiReactionThermo,
+    hePsiThermo,
+    reactingMixture,
+    canteraGasHThermoPhysics
+);
+
 makeReactionMixtureThermo
 (
     psiThermo,

src/thermophysicalModels/specie/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I/home/ignis/OpenFOAM/cantera/include \
+    -I/home/ignis/OpenFOAM/cantera/include/cantera \
+    -I$(LIB_SRC)/OpenFOAM/lnInclude
+
+LIB_LIBS = \
+    -L/home/ignis/OpenFOAM/cantera/lib \
+    -lcantera_shared \
+    -lOpenFOAM \
+    -lboost_system

src/thermophysicalModels/specie/include/reactionTypes.H

@@ -48,6 +48,8 @@ namespace Foam
 
     typedef Reaction<gasHThermoPhysics> gasHReaction;
 
+    typedef Reaction<canteraGasHThermoPhysics> canteraGasHReaction;
+
     typedef Reaction<constIncompressibleGasHThermoPhysics>
         constIncompressibleGasHReaction;
 

src/thermophysicalModels/specie/reaction/reactions/makeReactions.C

@@ -86,6 +86,7 @@ namespace Foam
     // sensible enthalpy based reactions
     makeReactions(constGasHThermoPhysics, constGasHReaction)
     makeReactions(gasHThermoPhysics, gasHReaction)
+    makeReactions(canteraGasHThermoPhysics, canteraGasHReaction)
     makeReactions
     (
         constIncompressibleGasHThermoPhysics,

src/thermophysicalModels/specie/transport/cantera/canteraTransport.C

@@ -26,14 +26,48 @@ License
 #include "canteraTransport.H"
 #include "IOstreams.H"
 
+#include "cantera/transport.h"
+#include "cantera/IdealGasMix.h"
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Thermo>
+typename Foam::canteraTransport<Thermo>::canteraManager Foam::canteraTransport<Thermo>::cm_;
+
+
+    template<class Thermo>
+Foam::canteraTransport<Thermo>::canteraManager::canteraManager()
+    :
+        gas_ (NULL),
+        tr_ (NULL),
+        nCanteraSp_ (0)
+{
+    try
+    {
+        gas_ = new Cantera::IdealGasMix("./cantera-transport.cti", "gas");
+        tr_ = Cantera::newTransportMgr("Mix", gas_);
+        nCanteraSp_ = gas_->nSpecies();
+    }
+    catch (Cantera::CanteraError e)
+    {
+        gas_ = NULL;
+        tr_ = NULL;
+        nCanteraSp_ = 0;
+    }
+}
+
+    template<class Thermo>
+Foam::canteraTransport<Thermo>::canteraManager::~canteraManager()
+{
+    delete tr_;
+    delete gas_;
+}
+
+    template<class Thermo>
 Foam::canteraTransport<Thermo>::canteraTransport(Istream& is)
-:
-    Thermo(is),
-    As_(readScalar(is)),
-    Ts_(readScalar(is))
+    :
+        Thermo(is),
+        X_(cm_.gas()->nSpecies())
 {
     is.check("canteraTransport<Thermo>::canteraTransport(Istream&)");
 }
@@ -43,9 +77,11 @@ template<class Thermo>
 Foam::canteraTransport<Thermo>::canteraTransport(const dictionary& dict)
 :
     Thermo(dict),
-    As_(readScalar(dict.subDict("transport").lookup("As"))),
-    Ts_(readScalar(dict.subDict("transport").lookup("Ts")))
-{}
+    X_(cm_.gas()->nSpecies(), 0.0)
+{
+    label si = cm_.gas()->speciesIndex(this->specie::name());
+    if (si >= 0) X_[si] = this->nMoles();
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -59,8 +95,6 @@ void Foam::canteraTransport<Thermo>::write(Ostream& os) const
     Thermo::write(os);
 
     dictionary dict("transport");
-    dict.add("As", As_);
-    dict.add("Ts", Ts_);
     os  << indent << dict.dictName() << dict;
 
     os  << decrIndent << token::END_BLOCK << nl;
@@ -75,7 +109,7 @@ Foam::Ostream& Foam::operator<<
     const canteraTransport<Thermo>& st
 )
 {
-    os << static_cast<const Thermo&>(st) << tab << st.As_ << tab << st.Ts_;
+    os << static_cast<const Thermo&>(st);
 
     os.check
     (
@@ -86,4 +120,29 @@ Foam::Ostream& Foam::operator<<
 }
 
 
+template<class Thermo>
+Foam::scalar Foam::canteraTransport<Thermo>::mu
+(
+ const scalar p,
+ const scalar T
+ ) const
+{
+    cm_.gas()->setState_TPX(T, p, X_.cdata());
+    scalar mu_ = cm_.tr()->viscosity();
+    return mu_;
+}
+
+
+template<class Thermo>
+inline Foam::scalar Foam::canteraTransport<Thermo>::kappa
+(
+ const scalar p, const scalar T
+ ) const
+{
+    cm_.gas()->setState_TPX(T, p, X_.cdata());
+    scalar kappa_ = cm_.tr()->thermalConductivity();
+    return kappa_;
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/specie/transport/cantera/canteraTransport.H

@@ -46,6 +46,14 @@ SourceFiles
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "scalarField.H"
+
+namespace Cantera
+{
+    class IdealGasMix;
+    class Transport;
+}
+
 namespace Foam
 {
 
@@ -98,21 +106,33 @@ class canteraTransport
 :
     public Thermo
 {
+    // Private class
+
+    class canteraManager
+    {
+        Cantera::IdealGasMix *gas_;
+        Cantera::Transport *tr_;
+        label nCanteraSp_;
+
+        public:
+
+        canteraManager();
+        ~canteraManager();
+        inline Cantera::IdealGasMix *gas() {return gas_;}
+        inline Cantera::Transport *tr() {return  tr_;}
+        inline label nCanteraSp() {return  nCanteraSp_;}
+    };
+
+
     // Private data
 
-        // Sutherland's coefficients
-        scalar As_, Ts_;
+        scalarField X_;
+
+        static canteraManager cm_;
 
 
     // Private Member Functions
 
-        //- Calculate the Sutherland coefficients
-        //  given two viscosities and temperatures
-        inline void calcCoeffs
-        (
-            const scalar mu1, const scalar T1,
-            const scalar mu2, const scalar T2
-        );
 
 
 public:
@@ -123,8 +143,7 @@ public:
         inline canteraTransport
         (
             const Thermo& t,
-            const scalar As,
-            const scalar Ts
+            const scalarField &X
         );
 
         //- Construct from two viscosities
@@ -163,10 +182,10 @@ public:
         }
 
         //- Dynamic viscosity [kg/ms]
-        inline scalar mu(const scalar p, const scalar T) const;
+        scalar mu(const scalar p, const scalar T) const;
 
         //- Thermal conductivity [W/mK]
-        inline scalar kappa(const scalar p, const scalar T) const;
+        scalar kappa(const scalar p, const scalar T) const;
 
         //- Thermal diffusivity of enthalpy [kg/ms]
         inline scalar alphah(const scalar p, const scalar T) const;

src/thermophysicalModels/specie/transport/cantera/canteraTransportI.H

@@ -27,21 +27,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class Thermo>
-inline void Foam::canteraTransport<Thermo>::calcCoeffs
-(
-    const scalar mu1, const scalar T1,
-    const scalar mu2, const scalar T2
-)
-{
-    scalar rootT1 = sqrt(T1);
-    scalar mu1rootT2 = mu1*sqrt(T2);
-    scalar mu2rootT1 = mu2*rootT1;
-
-    Ts_ = (mu2rootT1 - mu1rootT2)/(mu1rootT2/T1 - mu2rootT1/T2);
-
-    As_ = mu1*(1.0 + Ts_/T1)/rootT1;
-}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -50,13 +35,11 @@ template<class Thermo>
 inline Foam::canteraTransport<Thermo>::canteraTransport
 (
     const Thermo& t,
-    const scalar As,
-    const scalar Ts
+    const scalarField &X
 )
 :
     Thermo(t),
-    As_(As),
-    Ts_(Ts)
+    X_(X)
 {}
 
 
@@ -68,9 +51,9 @@ inline Foam::canteraTransport<Thermo>::canteraTransport
     const scalar mu2, const scalar T2
 )
 :
-    Thermo(t)
+    Thermo(t),
+    X_(cm_.nCanteraSp())
 {
-    calcCoeffs(mu1, T1, mu2, T2);
 }
 
 
@@ -82,8 +65,7 @@ inline Foam::canteraTransport<Thermo>::canteraTransport
 )
 :
     Thermo(name, st),
-    As_(st.As_),
-    Ts_(st.Ts_)
+    X_(st.X_)
 {}
 
 
@@ -128,28 +110,6 @@ Foam::canteraTransport<Thermo>::New
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Thermo>
-inline Foam::scalar Foam::canteraTransport<Thermo>::mu
-(
-    const scalar p,
-    const scalar T
-) const
-{
-    return As_*::sqrt(T)/(1.0 + Ts_/T);
-}
-
-
-template<class Thermo>
-inline Foam::scalar Foam::canteraTransport<Thermo>::kappa
-(
-    const scalar p, const scalar T
-) const
-{
-    scalar Cv_ = this->Cv(p, T);
-    return mu(p, T)*Cv_*(1.32 + 1.77*this->R()/Cv_);
-}
-
-
 template<class Thermo>
 inline Foam::scalar Foam::canteraTransport<Thermo>::alphah
 (
@@ -173,8 +133,7 @@ Foam::canteraTransport<Thermo>::operator=
 {
     Thermo::operator=(st);
 
-    As_ = st.As_;
-    Ts_ = st.Ts_;
+    X_ = st.X_;
 
     return *this;
 }
@@ -193,8 +152,10 @@ inline void Foam::canteraTransport<Thermo>::operator+=
     molr1 /= this->nMoles();
     scalar molr2 = st.nMoles()/this->nMoles();
 
-    As_ = molr1*As_ + molr2*st.As_;
-    Ts_ = molr1*Ts_ + molr2*st.Ts_;
+    // X_ = molr1*X_ + molr2*st.X_;
+    // X_ *= molr1;
+    // X_ += molr2*st.X_;
+    X_ += st.X_;
 }
 
 
@@ -211,8 +172,10 @@ inline void Foam::canteraTransport<Thermo>::operator-=
     molr1 /= this->nMoles();
     scalar molr2 = st.nMoles()/this->nMoles();
 
-    As_ = molr1*As_ - molr2*st.As_;
-    Ts_ = molr1*Ts_ - molr2*st.Ts_;
+    // X_ = molr1*X_ - molr2*st.X_;
+    // X_ *= molr1;
+    // X_ -= molr2*st.X_;
+    X_ -= st.X_;
 }
 
 
@@ -246,8 +209,8 @@ inline Foam::canteraTransport<Thermo> Foam::operator+
     return canteraTransport<Thermo>
     (
         t,
-        molr1*st1.As_ + molr2*st2.As_,
-        molr1*st1.Ts_ + molr2*st2.Ts_
+        // molr1*st1.X_ + molr2*st2.X_
+        st1.X_ + st2.X_
     );
 }
 
@@ -270,8 +233,8 @@ inline Foam::canteraTransport<Thermo> Foam::operator-
     return canteraTransport<Thermo>
     (
         t,
-        molr1*st1.As_ - molr2*st2.As_,
-        molr1*st1.Ts_ - molr2*st2.Ts_
+        // molr1*st1.X_ - molr2*st2.X_
+        st1.X_ - st2.X_
     );
 }
 
@@ -286,8 +249,7 @@ inline Foam::canteraTransport<Thermo> Foam::operator*
     return canteraTransport<Thermo>
     (
         s*static_cast<const Thermo&>(st),
-        st.As_,
-        st.Ts_
+        s*st.X_
     );
 }