diff --git a/Make/files b/Make/files
index 68e469c..bfc605f 100644
--- a/Make/files
+++ b/Make/files
@@ -3,6 +3,7 @@ diffusivityModel/Neutral/Neutral.C
diffusivityModel/Ion/Ion.C
diffusivityModel/Stockmayer/Stockmayer.C
+diffusivityModel/Coulomb/Coulomb.C
diffusivityModel/diffusivityModel/diffusivityModel.C
diff --git a/Make/options b/Make/options
index c4ab78a..e6da545 100644
--- a/Make/options
+++ b/Make/options
@@ -1,8 +1,10 @@
EXE_INC = \
+ -g \
-IdiffusivityModel/Particle \
-IdiffusivityModel/Neutral \
-IdiffusivityModel/Ion \
-IdiffusivityModel/Stockmayer \
+ -IdiffusivityModel/Coulomb \
-IdiffusivityModel/diffusivityModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
diff --git a/diffusivityModel/Coulomb/Coulomb.C b/diffusivityModel/Coulomb/Coulomb.C
new file mode 100644
index 0000000..138a6eb
--- /dev/null
+++ b/diffusivityModel/Coulomb/Coulomb.C
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 "Coulomb.H"
+
+#include "error.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::Coulomb::Coulomb(const Ion& a, const Ion& b)
+:
+ a_(a),
+ b_(b),
+ eij_(abs(a_.z()*b_.z())*sqr(a_.e)/(4.0*a_.pi*a_.eps0*a_.k)),
+ Wij_(1.0/(1.0/a_.W() + 1.0/b_.W()))
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::Coulomb::~Coulomb()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+
+
+// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
+
+void Foam::Coulomb::operator=(const Coulomb& rhs)
+{
+ // Check for assignment to self
+ if (this == &rhs)
+ {
+ FatalErrorInFunction
+ << "Attempted assignment to self"
+ << abort(FatalError);
+ }
+}
+
+// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Coulomb/Coulomb.H b/diffusivityModel/Coulomb/Coulomb.H
new file mode 100644
index 0000000..262b2c5
--- /dev/null
+++ b/diffusivityModel/Coulomb/Coulomb.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::Coulomb
+
+Description
+
+SourceFiles
+ CoulombI.H
+ Coulomb.C
+ CoulombIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Coulomb_H
+#define Coulomb_H
+
+#include "scalar.H"
+#include "scalarField.H"
+#include "Ion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class Coulomb;
+Istream& operator>>(Istream&, Coulomb&);
+Ostream& operator<<(Ostream&, const Coulomb&);
+
+
+/*---------------------------------------------------------------------------*\
+ Class Coulomb Declaration
+\*---------------------------------------------------------------------------*/
+
+class Coulomb
+{
+ // Private data
+
+ //- Collision participant a
+ const Ion &a_;
+
+ //- Collision participant b
+ const Ion &b_;
+
+ //- Combined well depth (depend on Length)
+ scalar eij_;
+
+ //- Combined molecular weight
+ scalar Wij_;
+
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ Coulomb(const Coulomb&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const Coulomb&);
+
+ //- Disallow default bitwise assignment
+ inline scalar power(const scalar a, const scalar b);
+
+ //- Disallow default bitwise assignment
+ inline scalar expon(const scalar a);
+
+ //- Disallow default bitwise assignment
+ inline scalar loge(const scalar a);
+
+
+public:
+
+ // Static data members
+
+
+ // Constructors
+
+ //- Construct null
+ Coulomb();
+
+ //- Construct from components
+ Coulomb(const Ion& a, const Ion& b);
+
+ //- Construct from Istream
+ Coulomb(Istream&);
+
+ //- Construct as copy
+ // Coulomb(const Coulomb&);
+
+
+ //- Destructor
+ ~Coulomb();
+
+
+ // Member Functions
+
+ // Access
+
+ // Check
+
+ // Edit
+
+ // Write
+
+ inline scalar Omega11(const scalar T);
+
+
+ inline scalar Omega22(const scalar T);
+
+
+ inline scalar D(const scalar p, const scalar T, const scalar rhoQc2);
+
+
+ inline scalar mu(const scalar p, const scalar T, const scalar rhoQc2);
+
+
+ inline scalar k(const scalar p, const scalar T);
+
+
+ // Member Operators
+
+
+ // Friend Functions
+
+ // Friend Operators
+
+ // IOstream Operators
+
+ friend Istream& operator>>(Istream&, Coulomb&);
+ friend Ostream& operator<<(Ostream&, const Coulomb&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CoulombI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/diffusivityModel/Coulomb/CoulombI.H b/diffusivityModel/Coulomb/CoulombI.H
new file mode 100644
index 0000000..f3dc1a2
--- /dev/null
+++ b/diffusivityModel/Coulomb/CoulombI.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::Coulomb::power(const scalar a, const scalar b)
+{
+ return pow(a, b);
+}
+
+
+inline Foam::scalar Foam::Coulomb::expon(const scalar a)
+{
+ return exp(a);
+}
+
+
+inline Foam::scalar Foam::Coulomb::loge(const scalar a)
+{
+ return log(a);
+}
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+
+inline Foam::scalar Foam::Coulomb::Omega11(const scalar Tstar)
+{
+ const scalar A = 0.5;
+ const scalar B = -0.129;
+ return (A * loge (Tstar) + B) / sqr(Tstar);
+}
+
+
+inline Foam::scalar Foam::Coulomb::Omega22(const scalar Tstar)
+{
+ const scalar A = 0.5;
+ const scalar B = 0.1165;
+ return (A * loge (Tstar) + B) / sqr(Tstar);
+}
+
+
+inline Foam::scalar Foam::Coulomb::D(const scalar p, const scalar T, const scalar rhoQc2)
+{
+ scalar sigmaij_ = sqrt(a_.eps0 * a_.k * T / (rhoQc2 + ROOTVSMALL)); // Debye length
+
+ return (3.0/8.0) * sqrt(2.0*a_.pi*a_.NA*pow3(a_.k*T)/Wij_)
+ / (a_.pi * sqr(sigmaij_) * p * Omega11(sigmaij_*T/eij_));
+}
+
+
+inline Foam::scalar Foam::Coulomb::mu(const scalar p, const scalar T, const scalar rhoQc2)
+{
+ scalar sigmaij_ = sqrt(a_.eps0 * a_.k * T / (rhoQc2 + ROOTVSMALL)); // Debye length
+
+ return (5.0/16.0) * sqrt(2.0*Wij_*a_.k*T/(a_.pi*a_.NA))
+ / (sqr(sigmaij_) * Omega22(sigmaij_*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 1235d2f..cb0b7b8 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -33,6 +33,7 @@ License
#include "Neutral.H"
#include "Ion.H"
#include "Stockmayer.H"
+#include "Coulomb.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -112,7 +113,8 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
D_(thermo_.composition().species().size()),
neutrals_(thermo_.composition().species().size()),
ions_(thermo_.composition().species().size()),
- nns_(thermo_.composition().species().size()*thermo_.composition().species().size())
+ nns_(thermo_.composition().species().size()*thermo_.composition().species().size()),
+ ccs_(thermo_.composition().species().size()*thermo_.composition().species().size())
{
const speciesTable &species_(thermo_.composition().species());
@@ -196,7 +198,7 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
nns_.resize(nNeutral*(nNeutral+1)/2);
// cns_.resize(nNeutral*nIon);
- // ccs_.resize((nIon+1)*(nIon+2)/2 - 1);
+ ccs_.resize((nIon+1)*(nIon+2)/2 - 1);
// ens_.resize(nNeutral);
// ee_;
@@ -249,6 +251,28 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
}
label iPair = 0;
+ forAll(ions_, i)
+ {
+ for (label j = i; j < ions_.size(); j++)
+ {
+ ccs_.set(iPair, new Coulomb(ions_[i], ions_[j]));
+ iPair++;
+ }
+ }
+
+ /*
+ iPair = 0;
+ forAll(ions_, i)
+ {
+ for (label j = 0; j < neutrals_.size(); j++)
+ {
+ cns_.set(iPair, new N64(ions_[i], neutrals_[j]));
+ iPair++;
+ }
+ }
+ */
+
+ iPair = 0;
forAll(neutrals_, i)
{
for (label j = i; j < neutrals_.size(); j++)
@@ -351,6 +375,8 @@ void Foam::diffusivityModel::correct()
const PtrList &Y(thermo_.composition().Y());
+ const volScalarField rhoQc2(thermo_.composition().Qc2() * rho);
+
scalarField Wpure(species_.size());
forAll (Wpure, i)
{
@@ -395,6 +421,7 @@ void Foam::diffusivityModel::correct()
const scalar pi = p[celli];
const scalar Ti = T[celli];
const scalar WbarI = Wbar[celli];
+ const scalar rhoQc2i = rhoQc2[celli];
forAll (species_, i)
{
@@ -458,11 +485,11 @@ void Foam::diffusivityModel::correct()
{
label i = I + 1;
- // muI[i] = nns_[idx1].mu(pi,Ti);
- muI[i] = SMALL;
+ muI[i] = ccs_[idx1].mu(pi,Ti,rhoQc2i);
+ // muI[i] = SMALL;
- // Dii[i] = nns_[idx1].D(pi,Ti);
- Dii[i] = SMALL;
+ Dii[i] = ccs_[idx1].D(pi,Ti,rhoQc2i);
+ // Dii[i] = SMALL;
Dij(i,i) = Dii[i];
idx1++;
@@ -546,6 +573,7 @@ void Foam::diffusivityModel::correct()
}
}
+ // Mixture average D, mu and k
if (pure)
{
UList Dpure(Dij[pureSpecieI], Dij.n());
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index ad67466..6337762 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -49,6 +49,7 @@ namespace Foam
class Neutral;
class Ion;
class Stockmayer;
+class Coulomb;
// Forward declaration of friend functions and operators
class diffusivityModel;
@@ -86,6 +87,9 @@ class diffusivityModel
//- Stockmayer potential list
PtrList nns_;
+ //- Stockmayer potential list
+ PtrList ccs_;
+
// Private Member Functions