From a63c28d73780e79a8b380e637c192c38d1b36d13 Mon Sep 17 00:00:00 2001 From: Yeongdo Park Date: Mon, 5 Nov 2018 00:12:28 -0500 Subject: [PATCH] correct well depth value in Zrot method of class Ion and fixed errors in methods mu and D of class Coulomb --- diffusivityModel/Coulomb/Coulomb.C | 10 ++++-- diffusivityModel/Coulomb/Coulomb.H | 24 +++++++++---- diffusivityModel/Coulomb/CoulombI.H | 39 ++++++++++++++++++--- diffusivityModel/Interaction/Interaction.H | 10 +++--- diffusivityModel/Interaction/InteractionI.H | 5 +++ diffusivityModel/Ion/Ion.C | 3 +- diffusivityModel/Ion/Ion.H | 10 ++++++ diffusivityModel/Ion/IonI.H | 11 +++++- 8 files changed, 90 insertions(+), 22 deletions(-) diff --git a/diffusivityModel/Coulomb/Coulomb.C b/diffusivityModel/Coulomb/Coulomb.C index 4ccb5e2..4695f7c 100644 --- a/diffusivityModel/Coulomb/Coulomb.C +++ b/diffusivityModel/Coulomb/Coulomb.C @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "Coulomb.H" +#include "Ion.H" #include "Electron.H" #include "error.H" @@ -44,7 +45,8 @@ License Foam::Coulomb::Coulomb(const Ion& a, const Ion& b) : - eij_(abs(a.z()*b.z())*sqr(a.e)/(4.0*Foam::Ion::pi*Foam::Ion::eps0*Foam::Ion::k)), + eij_(abs(a.z()*b.z())*sqr(a.e) + /(4.0*Foam::Ion::pi*Foam::Ion::eps0*Foam::Ion::k)), Wij_(inv(inv(a.W()) + inv(b.W()))) { } @@ -56,7 +58,8 @@ Foam::Coulomb::Coulomb(const Electron& a, const Ion& b) // 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())) - eij_(abs(a.z()*b.z())*sqr(a.e)/(4.0*Foam::Ion::pi*Foam::Ion::eps0*Foam::Ion::k)), + eij_(abs(a.z()*b.z())*sqr(a.e) + /(4.0*Foam::Ion::pi*Foam::Ion::eps0*Foam::Ion::k)), Wij_(inv(inv(a.W()) + inv(b.W()))) { } @@ -68,7 +71,8 @@ Foam::Coulomb::Coulomb(const Electron& a, const Electron& b) // 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())) - eij_(abs(a.z()*b.z())*sqr(a.e)/(4.0*Foam::Ion::pi*Foam::Ion::eps0*Foam::Ion::k)), + eij_(abs(a.z()*b.z())*sqr(a.e) + /(4.0*Foam::Ion::pi*Foam::Ion::eps0*Foam::Ion::k)), Wij_(inv(inv(a.W()) + inv(b.W()))) { } diff --git a/diffusivityModel/Coulomb/Coulomb.H b/diffusivityModel/Coulomb/Coulomb.H index a9f3914..1a0e3f9 100644 --- a/diffusivityModel/Coulomb/Coulomb.H +++ b/diffusivityModel/Coulomb/Coulomb.H @@ -38,9 +38,10 @@ SourceFiles #include "scalar.H" #include "scalarField.H" -#include "Ion.H" -// #include "Electron.H" + #include "Interaction.H" +#include "GasState.H" +#include "constants.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -53,6 +54,7 @@ class Ostream; // Forward declaration of friend functions and operators class Coulomb; +class Ion; class Electron; Istream& operator>>(Istream&, Coulomb&); Ostream& operator<<(Ostream&, const Coulomb&); @@ -123,21 +125,29 @@ public: // Write - inline scalar Omega11(const scalar T); + inline scalar Omega11(const scalar T) const ; - inline scalar Omega22(const scalar T); + inline scalar Omega22(const scalar T) const ; - inline scalar D(const scalar p, const scalar T, const scalar rhoQc2); + inline scalar D(const scalar p, const scalar T, const scalar rhoQc2) const ; - inline scalar mu(const scalar p, const scalar T, const scalar rhoQc2); + inline scalar mu(const scalar p, const scalar T, const scalar rhoQc2) const ; - inline scalar k(const scalar p, const scalar T); + inline scalar k(const scalar p, const scalar T) const ; + inline scalar DebyeLength(const GasState &state) const ; + + + inline scalar DebyeLength(const scalar T, const scalar rhoQc2) const ; + + + inline scalar epsK(const scalar T, const scalar rhoQc2) const ; + // Member Operators diff --git a/diffusivityModel/Coulomb/CoulombI.H b/diffusivityModel/Coulomb/CoulombI.H index bdd7af0..ccc086b 100644 --- a/diffusivityModel/Coulomb/CoulombI.H +++ b/diffusivityModel/Coulomb/CoulombI.H @@ -38,6 +38,7 @@ License inline Foam::scalar Foam::Coulomb::Omega11(const scalar Tstar) + const { constexpr scalar A = 0.5; constexpr scalar B = -0.129; @@ -46,6 +47,7 @@ inline Foam::scalar Foam::Coulomb::Omega11(const scalar Tstar) inline Foam::scalar Foam::Coulomb::Omega22(const scalar Tstar) + const { constexpr scalar A = 0.5; constexpr scalar B = 0.1165; @@ -53,19 +55,46 @@ inline Foam::scalar Foam::Coulomb::Omega22(const scalar Tstar) } -inline Foam::scalar Foam::Coulomb::D(const scalar p, const scalar T, const scalar rhoQc2) +inline Foam::scalar Foam::Coulomb::epsK(const scalar T, const scalar rhoQc2) + const { - scalar sigmaij_ = sqrt(Foam::Ion::eps0 * Foam::Ion::k * T / (rhoQc2 + ROOTVSMALL)); // Debye length + return eij_ / DebyeLength(T, rhoQc2); +} - return Interaction::D(Wij_, p, T, (Foam::Particle::pi * sqr(sigmaij_) * Omega11(T/eij_))); + +inline Foam::scalar Foam::Coulomb::DebyeLength(const scalar T, const scalar rhoQc2) + const +{ + return sqrt(Foam::constant::electromagnetic::epsilon0.value() + * Foam::constant::physicoChemical::k.value() + * T / (rhoQc2 + ROOTVSMALL)); +} + + +inline Foam::scalar Foam::Coulomb::DebyeLength(const GasState &state) + const +{ + return sqrt(Foam::constant::electromagnetic::epsilon0.value() + * Foam::constant::physicoChemical::k.value() + * state.T() / (state.rhoQc2() + ROOTVSMALL)); +} + + +inline Foam::scalar Foam::Coulomb::D(const scalar p, const scalar T, const scalar rhoQc2) + const +{ + scalar sigmaij_ = DebyeLength(T, rhoQc2); + + return Interaction::D(Wij_, p, T, (Foam::Particle::pi * sqr(sigmaij_) * Omega11(T/eij_/sigmaij_))); } inline Foam::scalar Foam::Coulomb::mu(const scalar p, const scalar T, const scalar rhoQc2) + const { - scalar sigmaij_ = sqrt(Foam::Ion::eps0 * Foam::Ion::k * T / (rhoQc2 + ROOTVSMALL)); // Debye length + scalar sigmaij_ = DebyeLength(T, rhoQc2); - return Interaction::mu(2.0*Wij_, T, (Foam::Particle::pi * sqr(sigmaij_) * Omega22(T/eij_))); + return Interaction::mu(2.0*Wij_, T, (Foam::Particle::pi * sqr(sigmaij_) * Omega22(T/eij_/sigmaij_))); } diff --git a/diffusivityModel/Interaction/Interaction.H b/diffusivityModel/Interaction/Interaction.H index 9ed9b4c..a9f95ca 100644 --- a/diffusivityModel/Interaction/Interaction.H +++ b/diffusivityModel/Interaction/Interaction.H @@ -76,19 +76,19 @@ class Interaction protected: //- Disallow default bitwise assignment - inline scalar power(const scalar a, const scalar b); + inline scalar power(const scalar a, const scalar b) const; //- Disallow default bitwise assignment - inline scalar expon(const scalar a); + inline scalar expon(const scalar a) const; //- Disallow default bitwise assignment - inline scalar loge(const scalar a); + inline scalar loge(const scalar a) const; //- Disallow default bitwise assignment - inline scalar D(const scalar mij, const scalar p, const scalar T, const scalar OmegaD); + inline scalar D(const scalar mij, const scalar p, const scalar T, const scalar OmegaD) const; //- Disallow default bitwise assignment - inline scalar mu(const scalar mk, const scalar T, const scalar OmegaD); + inline scalar mu(const scalar mk, const scalar T, const scalar OmegaD) const; public: diff --git a/diffusivityModel/Interaction/InteractionI.H b/diffusivityModel/Interaction/InteractionI.H index 3726675..c2e4ed6 100644 --- a/diffusivityModel/Interaction/InteractionI.H +++ b/diffusivityModel/Interaction/InteractionI.H @@ -28,30 +28,35 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // inline Foam::scalar Foam::Interaction::power(const scalar a, const scalar b) + const { return pow(a, b); } inline Foam::scalar Foam::Interaction::expon(const scalar a) + const { return exp(a); } inline Foam::scalar Foam::Interaction::loge(const scalar a) + const { return log(a); } inline Foam::scalar Foam::Interaction::D(const scalar mij, const scalar p, const scalar T, const scalar OmegaD) + const { return (3.0/16.0) * sqrt(2.0*Foam::Particle::pi*Foam::Particle::NA*pow3(Foam::Particle::k*T)/mij) / (p * OmegaD); } inline Foam::scalar Foam::Interaction::mu(const scalar mk, const scalar T, const scalar OmegaD) + const { return (5.0/16.0) * sqrt(Foam::Particle::pi*mk*Foam::Particle::k*T/Foam::Particle::NA) / OmegaD; } diff --git a/diffusivityModel/Ion/Ion.C b/diffusivityModel/Ion/Ion.C index 8a41eab..08efd21 100644 --- a/diffusivityModel/Ion/Ion.C +++ b/diffusivityModel/Ion/Ion.C @@ -52,7 +52,8 @@ Foam::Ion::Ion(const dictionary& thermoDict, const dictionary& tranDict) n_(tranDict.lookupOrDefault("nStockmayer", 12)), rctTarget_(tranDict.lookupOrDefault("rctTargets", hashedWordList())), Arct_(tranDict.lookupOrDefault("Arct", scalarList())), - Brct_(tranDict.lookupOrDefault("Brct", scalarList())) + Brct_(tranDict.lookupOrDefault("Brct", scalarList())), + selfInteraction(*this, *this) { if (z() > 0) { diff --git a/diffusivityModel/Ion/Ion.H b/diffusivityModel/Ion/Ion.H index 2880e64..d3a16d5 100644 --- a/diffusivityModel/Ion/Ion.H +++ b/diffusivityModel/Ion/Ion.H @@ -40,6 +40,10 @@ SourceFiles #include "hashedWordList.H" #include "scalarList.H" +#include "Coulomb.H" + +#include "GasState.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -51,6 +55,7 @@ class Ostream; // Forward declaration of friend functions and operators class Ion; +class Coulomb; Istream& operator>>(Istream&, Ion&); Ostream& operator<<(Ostream&, const Ion&); @@ -83,6 +88,9 @@ class Ion //- Resonant charge transfer parameter B scalarList Brct_; + //- Resonant charge transfer parameter B + Coulomb selfInteraction; + // Private Member Functions @@ -126,6 +134,8 @@ public: inline scalar Zrot(scalar T) const; + inline scalar Zrot(const GasState &state) const; + inline bool hasRct(const word &target) const; inline scalar Arct(const word &target) const; diff --git a/diffusivityModel/Ion/IonI.H b/diffusivityModel/Ion/IonI.H index 2a6d23b..a372c66 100644 --- a/diffusivityModel/Ion/IonI.H +++ b/diffusivityModel/Ion/IonI.H @@ -53,12 +53,21 @@ inline Foam::scalar Foam::Ion::C6() inline Foam::scalar Foam::Ion::Zrot(const scalar T) const { - scalar wellDepth_ = 0.0; + scalar wellDepth_ = selfInteraction.epsK (T, 0.0); return Particle::Zrot() * zeta(298.0, wellDepth_) / zeta(T, wellDepth_); } +inline Foam::scalar Foam::Ion::Zrot(const GasState &state) + const +{ + scalar wellDepth_ = selfInteraction.epsK (state.T(), state.rhoQc2()); + + return Particle::Zrot() * zeta(298.0, wellDepth_) / zeta(state.T(), wellDepth_); +} + + inline bool Foam::Ion::hasRct(const word &target) const {