correct well depth value in Zrot method of class Ion and fixed errors in methods mu and D of class Coulomb

This commit is contained in:
Yeongdo Park 2018-11-05 00:12:28 -05:00
parent 9f0fc16367
commit a63c28d737
8 changed files with 90 additions and 22 deletions

View file

@ -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())))
{
}

View file

@ -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

View file

@ -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_)));
}

View file

@ -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:

View file

@ -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;
}

View file

@ -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)
{

View file

@ -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;

View file

@ -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
{