622 lines
16 KiB
C
622 lines
16 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / 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 <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "diffusivityModel.H"
|
|
|
|
#include "IFstream.H"
|
|
#include "speciesTable.H"
|
|
#include "volFieldsFwd.H"
|
|
#include "scalarMatrices.H"
|
|
#include "Particle.H"
|
|
#include "Neutral.H"
|
|
#include "Ion.H"
|
|
#include "Stockmayer.H"
|
|
#include "Coulomb.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
// const dataType Foam::diffusivityModel::staticData();
|
|
|
|
|
|
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
Foam::scalar Foam::diffusivityModel::mixAvgD(const label k, const UList<scalar>& Di, const UList<scalar>& X, const UList<scalar>& Y)
|
|
{
|
|
scalarField XD((X)/Di);
|
|
scalarField YD((Y)/Di);
|
|
|
|
scalar XY = X[k] / (1. - Y[k]);
|
|
|
|
XD[k] = 0.0;
|
|
YD[k] = 0.0;
|
|
|
|
return 1. / (sum(XD) + XY*sum(YD));
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::diffusivityModel::mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W)
|
|
{
|
|
const label K = muI.size();
|
|
|
|
const scalar a = 1. / sqrt(8.);
|
|
|
|
scalarSquareMatrix Phi(K);
|
|
|
|
scalarField den(K);
|
|
|
|
for (label k = 0; k < K; k++)
|
|
{
|
|
for (label j = 0; j < K; j++)
|
|
{
|
|
const scalar Wj = W[j];
|
|
const scalar Wk = W[k];
|
|
|
|
const scalar muj = muI[j];
|
|
const scalar muk = muI[k];
|
|
|
|
Phi(k,j) = a * sqrt(1. + Wk/Wj) * sqr(1. + sqrt(muk/muj)*pow(Wj/Wk, 1./4.));
|
|
}
|
|
}
|
|
|
|
|
|
for (label k = 0; k < K; k++)
|
|
{
|
|
den[k] = sum(X * UList<scalar>(Phi[k], Phi.n()));
|
|
}
|
|
|
|
return sum(X * muI / den);
|
|
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::diffusivityModel::mixAvgK(const UList<scalar>& k, const UList<scalar>& X)
|
|
{
|
|
scalar sum1 = sum(k*X);
|
|
scalar sum2 = 1.0 / sum(X/k);
|
|
return (sum1 + sum2) / 2.0;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
|
|
:
|
|
thermo_(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()),
|
|
ccs_(thermo_.composition().species().size()*thermo_.composition().species().size())
|
|
{
|
|
|
|
const speciesTable &species_(thermo_.composition().species());
|
|
|
|
const volScalarField::Mesh &mesh = thermo_.composition().Y(0).mesh();
|
|
|
|
dictionary transports
|
|
(
|
|
IFstream
|
|
(
|
|
fileName(thermo_.lookup("foamTransportFile")).expand()
|
|
)()
|
|
);
|
|
|
|
Switch readIons = false;
|
|
Switch readNeutrals = false;
|
|
|
|
label nIon = 0;
|
|
label nNeutral = 0;
|
|
|
|
forAll(species_, i)
|
|
{
|
|
const word namei(species_[i]);
|
|
const label Wi = thermo_.composition().W(i);
|
|
const label zi = thermo_.composition().z(i);
|
|
|
|
dictionary tranDict(transports.subDict(species_[i]).subDict("transport"));
|
|
tranDict.add("name", namei);
|
|
tranDict.add("W", Wi);
|
|
tranDict.add("z", zi);
|
|
|
|
Particle p(tranDict);
|
|
|
|
// Info << tranDict << endl;
|
|
// Info << p << endl;
|
|
// Info << scalarList(tranDict["tranlib"]) << endl;
|
|
|
|
// Electron
|
|
if (species_[i] == "E-")
|
|
{
|
|
if (readIons || readNeutrals)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Electron should be the first species in the list"
|
|
<< abort(FatalError);
|
|
}
|
|
// new Electron (tranDict);
|
|
}
|
|
// Ions
|
|
else if (zi != 0)
|
|
{
|
|
if (readNeutrals)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Ions should be listed before any neutrals"
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
readIons = true;
|
|
|
|
ions_.set(nIon, new Ion(tranDict));
|
|
|
|
nIon++;
|
|
|
|
// new Ion (tranDict);
|
|
}
|
|
// Neutrals
|
|
else
|
|
{
|
|
readNeutrals = true;
|
|
|
|
neutrals_.set(nNeutral, new Neutral(tranDict));
|
|
|
|
nNeutral++;
|
|
}
|
|
}
|
|
|
|
ions_.resize(nIon);
|
|
|
|
neutrals_.resize(nNeutral);
|
|
nns_.resize(nNeutral*(nNeutral+1)/2);
|
|
|
|
// cns_.resize(nNeutral*nIon);
|
|
ccs_.resize((nIon+1)*(nIon+2)/2 - 1);
|
|
// ens_.resize(nNeutral);
|
|
// ee_;
|
|
|
|
forAll(species_, i)
|
|
{
|
|
const label zi = thermo_.composition().z(i);
|
|
|
|
// Electron
|
|
if (species_[i] == "E-")
|
|
{
|
|
// Create electron self diffusion interaction object
|
|
// new highOrderCoulomb (Electron, Electron);
|
|
|
|
for (label j = 1; j < nIon+1; j++)
|
|
{
|
|
// Create Coulomb potential interaction object
|
|
// new Coulomb (Electron, Ions[j]);
|
|
}
|
|
|
|
for (label j = nIon+1; j < species_.size(); j++)
|
|
{
|
|
// Create collision cross section interaction object
|
|
// new CrossSection (Electron, Neutrals[j]);
|
|
}
|
|
}
|
|
// Ions
|
|
else if (zi != 0)
|
|
{
|
|
for (label j = i; j < nIon+1; j++)
|
|
{
|
|
// Create Coulomb potential interaction object
|
|
// new Coulomb (Ions[i], Ions[j]);
|
|
}
|
|
|
|
for (label j = nIon+1; j < species_.size(); j++)
|
|
{
|
|
// Create (n,6,4) potential interaction object
|
|
// new n64 (Ions[i], Neutrals[j]);
|
|
}
|
|
}
|
|
// Neutrals
|
|
else
|
|
{
|
|
for (label j = i; j < species_.size(); j++)
|
|
{
|
|
// Create Stockmayer potential interaction object
|
|
// new Stockmayer (Neutrals[i], Neutrals[j]);
|
|
}
|
|
}
|
|
}
|
|
|
|
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++)
|
|
{
|
|
nns_.set(iPair, new Stockmayer(neutrals_[i], neutrals_[j]));
|
|
iPair++;
|
|
}
|
|
}
|
|
|
|
mu_.set
|
|
(
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
"mu",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("zero", dimArea/dimTime, 0.0)
|
|
)
|
|
);
|
|
|
|
k_.set
|
|
(
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
"lambda",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("zero", dimArea/dimTime, 0.0)
|
|
)
|
|
);
|
|
|
|
forAll(species_, i)
|
|
{
|
|
D_.set
|
|
(
|
|
i,
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
"D." + species_[i],
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("zero", dimArea/dimTime, 0.0)
|
|
)
|
|
);
|
|
}
|
|
|
|
}
|
|
|
|
Foam::diffusivityModel::diffusivityModel(const diffusivityModel& dm)
|
|
:
|
|
thermo_(dm.thermo_)
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
|
|
|
|
// Foam::autoPtr<Foam::diffusivityModel>
|
|
// Foam::diffusivityModel::New()
|
|
// {
|
|
// return autoPtr<diffusivityModel>(new diffusivityModel);
|
|
// }
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
|
|
|
Foam::diffusivityModel::~diffusivityModel()
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
|
|
|
void Foam::diffusivityModel::correct()
|
|
{
|
|
const volScalarField &T = thermo_.T();
|
|
const volScalarField &p = thermo_.p();
|
|
const volScalarField rho(thermo_.rho());
|
|
const volScalarField Cp(thermo_.Cp());
|
|
|
|
const speciesTable &species_(thermo_.composition().species());
|
|
|
|
const volScalarField Wbar(thermo_.composition().W());
|
|
|
|
const PtrList<volScalarField> &Y(thermo_.composition().Y());
|
|
|
|
const volScalarField rhoQc2(thermo_.composition().Qc2() * rho);
|
|
|
|
scalarField Wpure(species_.size());
|
|
forAll (Wpure, i)
|
|
{
|
|
Wpure[i] = thermo_.composition().W(i);
|
|
}
|
|
|
|
scalarSymmetricSquareMatrix Dij(species_.size());
|
|
|
|
scalarField localX(species_.size());
|
|
scalarField localY(species_.size());
|
|
scalarField localCv(species_.size());
|
|
|
|
scalarField localD(species_.size());
|
|
|
|
scalarField Dii(species_.size());
|
|
scalarField muI(species_.size());
|
|
scalarField kI(species_.size());
|
|
|
|
scalarField dof(species_.size());
|
|
scalarField::subField dofI(dof, ions_.size(), 1);
|
|
scalarField::subField dofN(dof, neutrals_.size(), ions_.size()+1);
|
|
|
|
dof[0] = 0;
|
|
forAll (ions_, i)
|
|
{
|
|
// dofI[i] = ions_[i].F();
|
|
dofI[i] = 0.0;
|
|
}
|
|
forAll (neutrals_, i)
|
|
{
|
|
dofN[i] = neutrals_[i].F();
|
|
}
|
|
|
|
scalarField Zrot(species_.size());
|
|
scalarField::subField ZrotI(Zrot, ions_.size(), 1);
|
|
scalarField::subField ZrotN(Zrot, neutrals_.size(), ions_.size()+1);
|
|
|
|
forAll (p, celli)
|
|
{
|
|
const scalar rhoi = rho[celli];
|
|
const scalar Cpi = Cp[celli];
|
|
const scalar pi = p[celli];
|
|
const scalar Ti = T[celli];
|
|
const scalar WbarI = Wbar[celli];
|
|
const scalar rhoQc2i = rhoQc2[celli];
|
|
|
|
forAll (species_, i)
|
|
{
|
|
localCv[i] = thermo_.composition().Cv(i,pi,Ti);
|
|
}
|
|
|
|
forAll (species_, i)
|
|
{
|
|
localY[i] = Y[i][celli];
|
|
localX[i] = Y[i][celli] * WbarI / Wpure[i];
|
|
}
|
|
|
|
Zrot[0] = 0;
|
|
forAll (ions_, i)
|
|
{
|
|
// ZrotI[i] = ions_[i].Zrot(Ti);
|
|
ZrotI[i] = 0.0;
|
|
}
|
|
forAll (neutrals_, i)
|
|
{
|
|
ZrotN[i] = neutrals_[i].Zrot(Ti);
|
|
}
|
|
|
|
// Electron - Electron
|
|
|
|
// muI[0] = nns_[idx].mu(pi,Ti);
|
|
muI[0] = SMALL;
|
|
|
|
// Dii[0] = nns_[idx].D(pi,Ti);
|
|
Dii[0] = SMALL;
|
|
Dij(0,0) = Dii[0];
|
|
|
|
label idx = 0;
|
|
// Electron - Ions
|
|
for (label J = 0; J < ions_.size(); J++)
|
|
{
|
|
label j = J + 1;
|
|
// Calculate Dij
|
|
// Dij(0,j) = nns_[idx].D(pi,Ti);
|
|
Dij(0,j) = SMALL;
|
|
Dij(j,0) = Dij(0,j);
|
|
idx++;
|
|
}
|
|
|
|
idx = 0;
|
|
// Electron - Neutrals
|
|
for (label J = 0; J < neutrals_.size(); J++)
|
|
{
|
|
label j = J + ions_.size() + 1;
|
|
// Calculate Dij
|
|
// Dij(0,j) = nns_[idx].D(pi,Ti);
|
|
Dij(0,j) = SMALL;
|
|
Dij(j,0) = Dij(0,j);
|
|
idx++;
|
|
}
|
|
|
|
// Ions
|
|
label idx1 = 0;
|
|
label idx2 = 0;
|
|
forAll (ions_, I)
|
|
{
|
|
label i = I + 1;
|
|
|
|
muI[i] = ccs_[idx1].mu(pi,Ti,rhoQc2i);
|
|
// muI[i] = SMALL;
|
|
|
|
Dii[i] = ccs_[idx1].D(pi,Ti,rhoQc2i);
|
|
// Dii[i] = SMALL;
|
|
Dij(i,i) = Dii[i];
|
|
idx1++;
|
|
|
|
// Ion - Ions
|
|
for (label J = I+1; J < ions_.size(); J++)
|
|
{
|
|
label j = J + 1;
|
|
// Calculate Dij
|
|
// Dij(i,j) = nns_[idx1].D(pi,Ti);
|
|
Dij(i,j) = SMALL;
|
|
Dij(j,i) = Dij(i,j);
|
|
idx1++;
|
|
}
|
|
|
|
// Ion - Neutrals
|
|
for (label j = ions_.size()+1; j < species_.size(); j++)
|
|
{
|
|
// Calculate Dij
|
|
// Dij(i,j) = nns_[idx2].D(pi,Ti);
|
|
Dij(i,j) = SMALL;
|
|
Dij(j,i) = Dij(i,j);
|
|
idx2++;
|
|
}
|
|
}
|
|
|
|
// Neutrals
|
|
idx = 0;
|
|
forAll (neutrals_, I)
|
|
{
|
|
label i = I + 1 + ions_.size();
|
|
|
|
muI[i] = nns_[idx].mu(pi,Ti);
|
|
|
|
Dii[i] = nns_[idx].D(pi,Ti);
|
|
Dij(i,i) = Dii[i];
|
|
idx++;
|
|
|
|
// Neutral - Neutrals
|
|
for (label j = i+1; j < species_.size(); j++)
|
|
{
|
|
// Calculate Dij
|
|
Dij(i,j) = nns_[idx].D(pi,Ti);
|
|
Dij(j,i) = Dij(i,j);
|
|
idx++;
|
|
}
|
|
}
|
|
|
|
|
|
// Pure Thermal conductivity
|
|
forAll (species_, i)
|
|
{
|
|
const scalar R = Neutral::RR / Wpure[i];
|
|
const scalar CvTrans = (3./2.)*R;
|
|
const scalar CvRot = (dof[i]/2.)*R;
|
|
const scalar CvVib = localCv[i] - CvTrans - CvRot;
|
|
|
|
const scalar rSc = rhoi * Dii[i] / muI[i];
|
|
|
|
const scalar A = 5./2. - rSc;
|
|
const scalar B = Zrot[i] + (2./Neutral::pi) * ((5./3.)*CvRot + rSc);
|
|
const scalar AB = (2./Neutral::pi)*(A/B);
|
|
|
|
const scalar fTrans = (5./2.) * (1.0 - AB * CvRot / CvTrans);
|
|
const scalar fRot = rSc*(1.0+AB);
|
|
const scalar fVib = rSc;
|
|
|
|
kI[i] = muI[i]*(fTrans*CvTrans + fRot*CvRot + fVib*CvVib);
|
|
|
|
}
|
|
|
|
// Pure Condition Test ( Xi = 1 )
|
|
Switch pure = false;
|
|
label pureSpecieI = -1;
|
|
forAll (species_, i)
|
|
{
|
|
if (1. - localX[i] < 1e-12)
|
|
{
|
|
pureSpecieI = i;
|
|
pure = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Mixture average D, mu and k
|
|
if (pure)
|
|
{
|
|
UList<scalar> Dpure(Dij[pureSpecieI], Dij.n());
|
|
|
|
forAll (species_, i)
|
|
{
|
|
D_[i][celli] = Dpure[i];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
forAll (species_, i)
|
|
{
|
|
UList<scalar> Di(Dij[i], Dij.n());
|
|
D_[i][celli] = mixAvgD(i, Di, localX, localY);
|
|
}
|
|
}
|
|
|
|
mu_()[celli] = mixAvgMu(muI, localX, Wpure);
|
|
k_()[celli] = mixAvgK(kI, localX) / rhoi / Cpi;
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
|
|
|
|
void Foam::diffusivityModel::operator=(const diffusivityModel& rhs)
|
|
{
|
|
// Check for assignment to self
|
|
if (this == &rhs)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Attempted assignment to self"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
|
|
// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
|
|
|
|
|
|
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
|
|
|
|
|
|
// ************************************************************************* //
|