From d95e1b3c17b1c0616d711de610369aacf125bf27 Mon Sep 17 00:00:00 2001 From: ignis Date: Mon, 14 May 2018 14:31:40 +0900 Subject: [PATCH] added coulomb interaction --- Make/files | 1 + Make/options | 2 + diffusivityModel/Coulomb/Coulomb.C | 83 +++++++++ diffusivityModel/Coulomb/Coulomb.H | 172 ++++++++++++++++++ diffusivityModel/Coulomb/CoulombI.H | 108 +++++++++++ .../diffusivityModel/diffusivityModel.C | 40 +++- .../diffusivityModel/diffusivityModel.H | 4 + 7 files changed, 404 insertions(+), 6 deletions(-) create mode 100644 diffusivityModel/Coulomb/Coulomb.C create mode 100644 diffusivityModel/Coulomb/Coulomb.H create mode 100644 diffusivityModel/Coulomb/CoulombI.H 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