class Stockmayer calculating transport coefs assuming Stockmayer potential function

.gitignore

@@ -0,0 +1,4 @@
+[._]*.s[a-v][a-z]
+[._]*.sw[a-p]
+[._]s[a-v][a-z]
+[._]sw[a-p]

Make/files

@@ -1,5 +1,8 @@
 diffusivityModel/Particle/Particle.C
 diffusivityModel/Neutral/Neutral.C
+
+diffusivityModel/Stockmayer/Stockmayer.C
+
 diffusivityModel/diffusivityModel/diffusivityModel.C
 
 eReactingFoam.C

Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -IdiffusivityModel/Particle \
     -IdiffusivityModel/Neutral \
+    -IdiffusivityModel/Stockmayer \
     -IdiffusivityModel/diffusivityModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

diffusivityModel/Neutral/Neutral.H

@@ -64,16 +64,16 @@ class Neutral
     // Private data
 
         //- Lenard-Jones parameter - potential well depth
-        //- epsilon / kB []
+        //- epsilon / kB [K]
         scalar wellDepth_;
 
-        //- Lenard-Jones parameter - collision diameter []
+        //- Lenard-Jones parameter - collision diameter [Angstrom]
         scalar diameter_;
 
-        //- Dipole moment []
+        //- Dipole moment [Debye]
         scalar dipoleMoment_;
 
-        //- Dipole polarizability []
+        //- Dipole polarizability [Angstrom^3]
         scalar alpha_;
 
         //- Quadrupole polarizability []
@@ -82,7 +82,7 @@ class Neutral
         //- Dispersion coefficient []
         scalar C6_;
 
-        //- Rotational relaxation collision number []
+        //- Rotational relaxation collision number [-]
         scalar Zrot_;
 
 

diffusivityModel/Neutral/NeutralI.H

@@ -46,14 +46,14 @@ inline Foam::scalar Foam::Neutral::wellDepth()
 inline Foam::scalar Foam::Neutral::diameter()
     const
 {
-    return diameter_;
+    return diameter_ * Angstrom;
 }
 
 
 inline Foam::scalar Foam::Neutral::dipoleMoment()
     const
 {
-    return dipoleMoment_;
+    return dipoleMoment_ * Debye;
 }
 
 

diffusivityModel/Particle/Particle.C

@@ -28,10 +28,33 @@ License
 #include "dictionary.H"
 #include "scalarList.H"
 
+#include "constants.H"
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 // const dataType Foam::Particle::staticData();
 
+//- Universal gas constant (default in [J/(kmol K)])
+const Foam::scalar Foam::Particle::RR = constant::physicoChemical::R.value()*1000;
+
+//- Elementary charge (default in [C])
+const Foam::scalar Foam::Particle::e = constant::electromagnetic::e.value();
+
+//- Avogadro number (default in [1/mol])
+const Foam::scalar Foam::Particle::NA = constant::physicoChemical::NA.value()*1000;
+
+//- Boltzmann constant (default in [J/K])
+const Foam::scalar Foam::Particle::k = constant::physicoChemical::k.value();
+
+//- 1 Angstom in meter
+const Foam::scalar Foam::Particle::Angstrom = 1e-10;
+
+//- 1 Debye in [ Coulomb x meter ]
+const Foam::scalar Foam::Particle::Debye = 1.0 / 2.99792458e29;
+
+//- Boltzmann constant (default in [J/K])
+const Foam::scalar Foam::Particle::pi = constant::mathematical::pi;
+
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 

diffusivityModel/Particle/Particle.H

@@ -103,6 +103,27 @@ public:
 
     // Static data members
 
+        //- Universal gas constant (default in [J/(kmol K)])
+        static const scalar RR;
+
+        //- Elementary charge (default in [C])
+        static const scalar e;
+
+        //- Avogadro number (default in [1/mol])
+        static const scalar NA;
+
+        //- Boltzmann constant (default in [J/K])
+        static const scalar k;
+
+        //- 1 Angstom in meter
+        static const scalar Angstrom;
+
+        //- 1 Angstom in meter
+        static const scalar Debye;
+
+        //- Circumference / Diameter
+        static const scalar pi;
+
 
     // Constructors
 

diffusivityModel/Stockmayer/Stockmayer.C

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Stockmayer.H"
+
+#include "error.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+/*
+const Foam::scalarField Foam::Stockmayer::a(
+        autoPtr<IStringStream>(
+            new IStringStream(
+                "6 (1.0548 0.15504 0.55909 2.1705 0.093193 1.5)"
+                )
+            )()
+        );
+
+const Foam::scalarField Foam::Stockmayer::b(
+        autoPtr<IStringStream>(
+            new IStringStream(
+                "6 (1.0413 0.11930 0.43628 1.6041 0.095661 2.0)"
+                )
+            )()
+        );
+*/
+
+const Foam::scalar Foam::Stockmayer::a[6] = 
+{ 1.0548, 0.15504, 0.55909, 2.1705, 0.093193, 1.5 };
+
+const Foam::scalar Foam::Stockmayer::b[6] = 
+{ 1.0413, 0.11930, 0.43628, 1.6041, 0.095661, 2.0 };
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Stockmayer::Stockmayer(const Neutral& a, const Neutral& b)
+:
+    a_(a.dipoleMoment() > 0.0 ? b : a),
+    b_(a.dipoleMoment() > 0.0 ? a : b),
+    eij_(Foam::sqrt(a_.wellDepth()*b_.wellDepth())),
+    sigmaij_((a_.diameter()+b_.diameter())/2.0),
+    Wij_(1.0/(1.0/a_.W() + 1.0/b_.W())),
+    deltaij_(a_.dipoleMoment()*b_.dipoleMoment()/2.0/(eij_*a_.k)/pow3(sigmaij_))
+{
+    if (b_.dipoleMoment() > 0.0 && deltaij_ > 0.0)
+    {
+        scalar alphaStar = a_.alpha() / a_.diameter() / a_.diameter() / a_.diameter();
+
+        scalar muStar = b_.dipoleMoment() / Foam::sqrt(b_.wellDepth()*b_.diameter()*b_.diameter()*b_.diameter());
+
+        scalar ksi = 1.0 + alphaStar * muStar * Foam::sqrt(b_.wellDepth()/a_.wellDepth());
+
+        eij_ *= ksi * ksi;
+        sigmaij_ *= Foam::pow(ksi, -1.0/6.0);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::Stockmayer::~Stockmayer()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::Stockmayer::operator=(const Stockmayer& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+}
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

diffusivityModel/Stockmayer/Stockmayer.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::Stockmayer
+
+Description
+
+SourceFiles
+    StockmayerI.H
+    Stockmayer.C
+    StockmayerIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Stockmayer_H
+#define Stockmayer_H
+
+#include "scalar.H"
+#include "scalarField.H"
+#include "Neutral.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class Stockmayer;
+Istream& operator>>(Istream&, Stockmayer&);
+Ostream& operator<<(Ostream&, const Stockmayer&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class Stockmayer Declaration
+\*---------------------------------------------------------------------------*/
+
+class Stockmayer
+{
+    // Private data
+
+        //- Collision participant a
+        const Neutral &a_;
+
+        //- Collision participant b
+        const Neutral &b_;
+
+        //- Combined well depth
+        scalar eij_;
+
+        //- Combined diameter
+        scalar sigmaij_;
+
+        //- Combined molecular weight
+        scalar Wij_;
+
+        //- Stockmayer potential parameter
+        scalar deltaij_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        Stockmayer(const Stockmayer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const Stockmayer&);
+
+        //- Disallow default bitwise assignment
+        inline scalar power(const scalar a, const scalar b);
+
+        //- Disallow default bitwise assignment
+        inline scalar expon(const scalar a);
+
+
+public:
+
+    // Static data members
+
+        //- Coefficients for Omega(1,1) fitting
+        static const scalar a[6];
+
+        //- Coefficients for Omega(2,2) fitting
+        static const scalar b[6];
+
+
+    // Constructors
+
+        //- Construct null
+        Stockmayer();
+
+        //- Construct from components
+        Stockmayer(const Neutral& a, const Neutral& b);
+
+        //- Construct from Istream
+        Stockmayer(Istream&);
+
+        //- Construct as copy
+        // Stockmayer(const Stockmayer&);
+
+
+    //- Destructor
+    ~Stockmayer();
+
+
+    // Member Functions
+
+        // Access
+
+        // Check
+
+        // Edit
+
+        // Write
+
+        inline scalar Omega11(const scalar T);
+
+
+        inline scalar Omega22(const scalar T);
+
+
+        inline scalar f11(const scalar T);
+
+
+        inline scalar f22(const scalar T);
+
+
+        inline scalar D(const scalar p, const scalar T);
+
+
+        inline scalar mu(const scalar p, const scalar T);
+
+
+        inline scalar k(const scalar p, const scalar T);
+
+
+    // Member Operators
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, Stockmayer&);
+        friend Ostream& operator<<(Ostream&, const Stockmayer&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "StockmayerI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

diffusivityModel/Stockmayer/StockmayerI.H

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline Foam::scalar Foam::Stockmayer::power(const scalar a, const scalar b)
+{
+    return pow(a, b);
+}
+
+inline Foam::scalar Foam::Stockmayer::expon(const scalar a)
+{
+    return exp(a);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::Stockmayer::Omega11(const scalar Tstar)
+{
+    return a[0] * power(Tstar, -a[1]) + power(Tstar + a[2], -a[3]);
+}
+
+
+inline Foam::scalar Foam::Stockmayer::Omega22(const scalar Tstar)
+{
+    return b[0] * power(Tstar, -b[1]) + power(Tstar + b[2], -b[3]);
+}
+
+
+inline Foam::scalar Foam::Stockmayer::f11(const scalar Tstar)
+{
+    return 1.0 + sqr(deltaij_)*(expon(a[4]/Tstar)-expon(-a[5]/Tstar))/(2 + 2.5*deltaij_);
+}
+
+
+inline Foam::scalar Foam::Stockmayer::f22(const scalar Tstar)
+{
+    return 1.0 + sqr(deltaij_)*(expon(b[4]/Tstar)-expon(-b[5]/Tstar))/(2 + 2.5*deltaij_);
+}
+
+
+inline Foam::scalar Foam::Stockmayer::D(const scalar p, const scalar T)
+{
+    return (3.0/8.0) * sqrt(2.0*a_.pi*a_.NA*pow3(a_.k*T)/Wij_)
+        / (a_.pi * sqr(sigmaij_) * p * Omega11(T/eij_));
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

diffusivityModel/diffusivityModel/diffusivityModel.C

@@ -31,6 +31,7 @@ License
 #include "scalarMatrices.H"
 #include "Particle.H"
 #include "Neutral.H"
+#include "Stockmayer.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -51,7 +52,9 @@ License
 Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
 :
     thermo_(thermo),
-    D_(thermo_.composition().species().size())
+    D_(thermo_.composition().species().size()),
+    neutrals_(thermo_.composition().species().size()),
+    nns_(thermo_.composition().species().size()*thermo_.composition().species().size())
 {
 
     const speciesTable &species_(thermo_.composition().species());
@@ -69,6 +72,9 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
     Switch readIons = false;
     Switch readNeutrals = false;
 
+    label nIon = 0;
+    label nNeutral = 0;
+
     forAll(species_, i)
     {
         const word namei(species_[i]);
@@ -108,6 +114,7 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
             }
 
             readIons = true;
+            nIon++;
 
             // new Ion (tranDict);
         }
@@ -116,12 +123,72 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
         {
             readNeutrals = true;
 
-            Neutral n(tranDict);
-            Info << n << endl;
+            neutrals_.set(nNeutral, new Neutral(tranDict));
+
+            nNeutral++;
+        }
+    }
+
+    neutrals_.resize(nNeutral);
+    nns_.resize(nNeutral*(nNeutral+1)/2);
+
+    forAll(species_, i)
+    {
+        const label zi = thermo_.composition().z(i);
+
+        // Electron
+        if (species_[i] == "E-")
+        {
+            // Create electron self diffusion interaction object
+            // new highOrderCoulomb (Electron, Electron);
+
+            for (label j = 1; j < nIon+1; j++)
+            {
+                // Create Coulomb potential interaction object
+                // new Coulomb (Electron, Ions[j]);
+            }
+
+            for (label j = nIon+1; j < species_.size(); j++)
+            {
+                // Create collision cross section interaction object
+                // new CrossSection (Electron, Neutrals[j]);
+            }
+        }
+        // Ions
+        else if (zi != 0)
+        {
+            for (label j = i; j < nIon+1; j++)
+            {
+                // Create Coulomb potential interaction object
+                // new Coulomb (Ions[i], Ions[j]);
+            }
+
+            for (label j = nIon+1; j < species_.size(); j++)
+            {
+                // Create (n,6,4) potential interaction object
+                // new n64 (Ions[i], Neutrals[j]);
+            }
+        }
+        // Neutrals
+        else
+        {
+            for (label j = i; j < species_.size(); j++)
+            {
+                // Create Stockmayer potential interaction object
+                // new Stockmayer (Neutrals[i], Neutrals[j]);
+            }
+        }
+    }
+
+    label iPair = 0;
+    forAll(neutrals_, i)
+    {
+        for (label j = i; j < neutrals_.size(); j++)
+        {
+            nns_.set(iPair, new Stockmayer(neutrals_[i], neutrals_[j]));
+            iPair++;
         }
     }
-/*
-*/
 
     forAll(species_, i)
     {
@@ -178,15 +245,27 @@ void Foam::diffusivityModel::correct()
 
     scalarSymmetricSquareMatrix Dij(species_.size());
 
+    const volScalarField Cv(thermo_.Cv());
+    const volScalarField Wbar(thermo_.composition().W());
+
+    const PtrList<volScalarField> &Y(thermo_.composition().Y());
+
+    scalarField Wpure(species_.size());
+    forAll (Wpure, i)
+    {
+        Wpure[i] = thermo_.composition().W(i);
+    }
+
     forAll (p, celli)
     {
         const scalar pi = p[celli];
         const scalar Ti = T[celli];
-        scalarField Yi(species_.size());
+        const scalar Cvi = Cv[celli];
 
+        scalarField localY(species_.size());
         forAll (species_, i)
         {
-            Yi[i] = thermo_.composition().Y(i)[celli];
+            localY[i] = Y[i][celli];
         }
 
         label idx = 0;
@@ -197,7 +276,7 @@ void Foam::diffusivityModel::correct()
             {
                 // Calculate Dij
                 // Dij[idx] = interactions[idx].D(pi,Ti);
-                Dij(i,j) = idx;
+                Dij(i,j) = nns_[idx].D(pi,Ti);
                 Dij(j,i) = Dij(i,j);
                 idx++;
             }
@@ -206,7 +285,7 @@ void Foam::diffusivityModel::correct()
 
         }
 
-        Info << Yi << endl;
+        Info << localY << endl;
     }
 
     return;

diffusivityModel/diffusivityModel/diffusivityModel.H

@@ -46,6 +46,8 @@ namespace Foam
 {
 
 // Forward declaration of classes
+class Neutral;
+class Stockmayer;
 
 // Forward declaration of friend functions and operators
 class diffusivityModel;
@@ -62,9 +64,21 @@ class diffusivityModel
         //- Reference to thermo object
         const psiReactionThermo& thermo_;
 
-        //- Mixture averaged diffusivities field
+        //- Mixture averaged diffusivity fields
         PtrList<volScalarField> D_;
 
+        //- Electron object
+        // autoPtr<Electron> electron_;
+
+        //- Ion object list
+        // PtrList<Ion> ions_;
+
+        //- Neutral object list
+        PtrList<Neutral> neutrals_;
+
+        //- Stockmayer potential list
+        PtrList<Stockmayer> nns_;
+
 
     // Private Member Functions