research/eReactingFoam-4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

n,6,4 resonant charge transfer, fit intercepts, powers of gamma_ij precalculation

Browse source
This commit is contained in:
ignis 2018-06-10 03:12:26 +09:00
parent 67fe50953d
commit 2b44e65696
5 changed files with 147 additions and 14 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
diffusivityModel/Ion/Ion.H
View file

@ -131,6 +131,8 @@ public:
inline scalar Zrot(scalar T) const;
inline bool hasRct(const word &target) const;
inline scalar Arct(const word &target) const;
inline scalar Brct(const word &target) const;

7
diffusivityModel/Ion/IonI.H
View file

@ -73,6 +73,13 @@ inline Foam::scalar Foam::Ion::Zrot(const scalar T)
}
inline bool Foam::Ion::hasRct(const word &target)
const
{
return rctTarget_.contains(target);
}
inline Foam::scalar Foam::Ion::Arct(const word &target)
const
{

93
diffusivityModel/N64/N64.C
View file

@ -26,18 +26,20 @@ License
#include "N64.H"
#include "error.H"
#include "scalarMatrices.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::scalar Foam::N64::a11[(M+1)*(M+2)/2] = {
0., -0.81508875, -0.06484549, 0.31926761, 0.35164219, -0.01484638,
1.21327441805, -0.81508875, -0.06484549, 0.31926761, 0.35164219, -0.01484638,
-0.02385399, -0.14337633, -0.18493828, 0.05466979, -0.0010004, 0.01936938,
-0.0130214, 0.12179535, -0.05237791, -0.00172218, -0.00447904, 0.02211965,
0.00964791, -0.07566612, 0.02527412
};
const Foam::scalar Foam::N64::a22[(M+1)*(M+2)/2] = {
0.00000000e+00, -8.59231154e-01, -7.02815662e-02, 2.79760234e-01,
1.3678302569, -8.59231154e-01, -7.02815662e-02, 2.79760234e-01,
3.59026813e-01, -3.25082644e-01, -8.27621218e-03, -1.12825149e-01,
-2.07675933e-01, 1.02433682e+00, -4.59772292e-04, 1.03264160e-02,
-1.35152522e-02, 1.69188798e-01, -1.21839345e+00, -1.88957374e-03,
@ -45,12 +47,19 @@ const Foam::scalar Foam::N64::a22[(M+1)*(M+2)/2] = {
5.07723424e-01
};
Foam::scalar Foam::N64::vec[(M+1)*(M+2)/2] = { 0.0 };
//Foam::scalar Foam::N64::vec[(M+1)*(M+2)/2] = { 0.0 };
Foam::scalarField Foam::N64::vec((M+1), 0.0);
const Foam::scalar Foam::N64::digamma3 = 0.922784335098467139393487909917597568957840664060;
const Foam::scalar Foam::N64::trigamma3 = sqr(Ion::pi) / 6. - 5./4.;
/*
const Foam::scalar Foam::N64::digamma3 = 0.922784335098;
const Foam::scalar Foam::N64::trigamma3 = 0.394934066848;
*/
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@ -71,15 +80,89 @@ Foam::N64::N64(const Ion& a, const Neutral& b)
eij_(5.2 * b_.alpha() * sqr(a_.z()) * (1.0 + ksiij_) / pow4(sigmaij_) / (a_.k/a_.e)),
C6ij_(0.0),
gammaij_(0.0),
Wij_(1.0/(1.0/a_.W() + 1.0/b_.W()))
Wij_(1.0/(1.0/a_.W() + 1.0/b_.W())),
b11_((M+1), 0.0),
b22_((M+1), 0.0)
{
scalar aij = a_.alpha() / b_.alpha();
scalar aji = 1.0 / aij;
C6ij_ = 2.0 * a_.C6() * b_.C6() / (aji * a_.C6() + aij * b_.C6());
gammaij_ = ((2.0/sqr(a_.z())) * C6ij_ + b_.alphaQ()) / (b_.alpha() * sqr(sigmaij_));
Info << "N64, gamma = " << gammaij_ << endl;
Info << a_.name() + " = " + b_.name() + " N64, gamma = " << gammaij_ << endl;
scalarField temp((M+1)*(M+2)/2, 0.0); // temporary feature vector with fixed gammaij value
transform(1.0, gammaij_, temp.data());
Info << "fixed gammaij feature" << endl << temp << endl;
scalarField t11 (temp * UList<scalar>(const_cast<scalar*>(a11), temp.size()));
scalarField t22 (temp * UList<scalar>(const_cast<scalar*>(a22), temp.size()));
scalarSquareMatrix mat11(M+1, 0.0);
scalarSquareMatrix mat22(M+1, 0.0);
label index = 0;
for (label s = 0; s < M+1; s++)
{
for (label i = 0; i <= s; i++)
{
mat11(i, s-i) = t11[index];
mat22(i, s-i) = t22[index];
index++;
}
}
for (label s = 0; s < M+1; s++)
{
b11_ += UList<scalar>(mat11[s], M+1);
b22_ += UList<scalar>(mat22[s], M+1);
}
if (false) // (a.hasRct(b.name()))
{
Info << "RCT" << endl;
scalar A = a.Arct(b.name());
scalar B = a.Brct(b.name());
scalar C = A/B - digamma3 - loge(eij_);
scalar D = sqr(B) / (a_.pi * sqr(sigmaij_)); // Formula in [A^2] (Angstrom)
scalar a0 = D * (trigamma3 + sqr(C));
scalar a1 = D * (-2.0 * C);
scalar a2 = D;
b11_[0] += a0;
b11_[1] += a1;
b11_[2] += a2;
}
/*
for (label s = 0; s < M+1; s++)
{
Info << UList<scalar>(mat11[s], M+1) << endl;
}
Info << b11_ << endl;
for (label s = 0; s < M+1; s++)
{
Info << UList<scalar>(mat22[s], M+1) << endl;
}
Info << b22_ << endl;
*/
}

51
diffusivityModel/N64/N64.H
View file

@ -93,11 +93,20 @@ class N64
//- Combined molecular weight
scalar Wij_;
//- Omega(1,1) Fit Coefficient with fixed gamma_ij
scalarField b11_;
//- Omega(2,2) Fit Coefficient with fixed gamma_ij
scalarField b22_;
//- Construct null
static const scalar a11[(M+1)*(M+2)/2];
static const scalar a22[(M+1)*(M+2)/2];
static scalar vec[(M+1)*(M+2)/2];
// static scalar vec[(M+1)*(M+2)/2];
static scalarField vec;
//- Construct null
static constexpr scalar K1 = 1.767;
@ -125,9 +134,28 @@ class N64
// primary template
template <label Order>
class BivarPolyFeat
struct UnivarPolyFeat
{
static scalar transform (scalar x1, scalar *vec)
{
scalar highest = UnivarPolyFeat<Order-1>::transform(x1,vec) * x1;
*(vec+Order) = highest;
return highest;
}
};
// convenience function
inline void transform (scalar x1, scalar *vec)
{
UnivarPolyFeat<M>::transform(x1,vec);
}
// primary template
template <label Order>
struct BivarPolyFeat
{
public:
static scalar *transform (scalar x1, scalar x2, scalar *vec)
{
scalar *src = BivarPolyFeat<Order-1>::transform(x1,x2,vec);
@ -149,6 +177,8 @@ class N64
BivarPolyFeat<M>::transform(x1,x2,vec);
}
// primary template
template <label N>
struct ScalarProduct
{
@ -229,6 +259,17 @@ public:
};
// partial specialization as end criteria
template <>
class N64::UnivarPolyFeat<0>
{
public:
static scalar transform (scalar x1, scalar *vec)
{
*vec = 1.0;
return 1.0;
}
};
template <>
class N64::BivarPolyFeat<0>
{
@ -241,10 +282,10 @@ class N64::BivarPolyFeat<0>
};
template <>
struct N64::ScalarProduct<0>
struct N64::ScalarProduct<1>
{
static scalar product (const scalar* first, const scalar* second) {
return 0.0;
return *first * *second;
};
};

8
diffusivityModel/N64/N64I.H
View file

@ -59,15 +59,15 @@ inline Foam::scalar Foam::N64::bivarPolyFeat(const scalar Tstar)
inline Foam::scalar Foam::N64::Omega11(const scalar Tstar)
{
transform (loge(Tstar), gammaij_, vec);
return scalarProduct (a11, vec);
transform (loge(Tstar), vec.data());
return ScalarProduct<M+1>::product (b11_.data(), vec.data());
}
inline Foam::scalar Foam::N64::Omega22(const scalar Tstar)
{
transform (loge(Tstar), gammaij_, vec);
return scalarProduct (a22, vec);
transform (loge(Tstar), vec.data());
return ScalarProduct<M+1>::product (b22_.data(), vec.data());
}