diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..ebb0a8e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+[._]*.s[a-v][a-z]
+[._]*.sw[a-p]
+[._]s[a-v][a-z]
+[._]sw[a-p]
diff --git a/Make/files b/Make/files
index 43f1f5a..f2c3a82 100644
--- a/Make/files
+++ b/Make/files
@@ -1,5 +1,8 @@
diffusivityModel/Particle/Particle.C
diffusivityModel/Neutral/Neutral.C
+
+diffusivityModel/Stockmayer/Stockmayer.C
+
diffusivityModel/diffusivityModel/diffusivityModel.C
eReactingFoam.C
diff --git a/Make/options b/Make/options
index f223db6..71ab37f 100644
--- a/Make/options
+++ b/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 \
diff --git a/diffusivityModel/Neutral/Neutral.H b/diffusivityModel/Neutral/Neutral.H
index 4a438d8..f707a2d 100644
--- a/diffusivityModel/Neutral/Neutral.H
+++ b/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_;
diff --git a/diffusivityModel/Neutral/NeutralI.H b/diffusivityModel/Neutral/NeutralI.H
index e47be06..df32e98 100644
--- a/diffusivityModel/Neutral/NeutralI.H
+++ b/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;
}
diff --git a/diffusivityModel/Particle/Particle.C b/diffusivityModel/Particle/Particle.C
index 695fb45..840448d 100644
--- a/diffusivityModel/Particle/Particle.C
+++ b/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 * * * * * * * * * * * * //
diff --git a/diffusivityModel/Particle/Particle.H b/diffusivityModel/Particle/Particle.H
index 81492c7..7c06456 100644
--- a/diffusivityModel/Particle/Particle.H
+++ b/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
diff --git a/diffusivityModel/Stockmayer/Stockmayer.C b/diffusivityModel/Stockmayer/Stockmayer.C
new file mode 100644
index 0000000..7086a16
--- /dev/null
+++ b/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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "Stockmayer.H"
+
+#include "error.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+/*
+const Foam::scalarField Foam::Stockmayer::a(
+ autoPtr(
+ new IStringStream(
+ "6 (1.0548 0.15504 0.55909 2.1705 0.093193 1.5)"
+ )
+ )()
+ );
+
+const Foam::scalarField Foam::Stockmayer::b(
+ autoPtr(
+ 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 * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Stockmayer/Stockmayer.H b/diffusivityModel/Stockmayer/Stockmayer.H
new file mode 100644
index 0000000..4dacb50
--- /dev/null
+++ b/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 .
+
+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
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Stockmayer/StockmayerI.H b/diffusivityModel/Stockmayer/StockmayerI.H
new file mode 100644
index 0000000..b609061
--- /dev/null
+++ b/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 .
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * 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 * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index b276449..bd407ae 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/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 &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;
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index 594e940..f7fc2fa 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/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 D_;
+ //- Electron object
+ // autoPtr electron_;
+
+ //- Ion object list
+ // PtrList ions_;
+
+ //- Neutral object list
+ PtrList neutrals_;
+
+ //- Stockmayer potential list
+ PtrList nns_;
+
// Private Member Functions