737 lines
17 KiB
C
737 lines
17 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 "multiComponentMixture.H"
|
|
|
|
#include "Particle.H"
|
|
#include "Electron.H"
|
|
#include "Neutral.H"
|
|
#include "Ion.H"
|
|
#include "Stockmayer.H"
|
|
#include "Coulomb.H"
|
|
#include "N64.H"
|
|
#include "CrossSection.H"
|
|
#include "GasState.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
const Foam::PtrList<Foam::gasHThermoPhysics> *Foam::diffusivityModel::thermoData_ = NULL;
|
|
|
|
|
|
// * * * * * * * * * * * * * 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));
|
|
}
|
|
|
|
|
|
void
|
|
Foam::diffusivityModel::mixAvgDi(UList<scalar>& Di, const scalarSymmetricSquareMatrix &Dij, const UList<scalar>& X, const UList<scalar>& Y)
|
|
{
|
|
|
|
// Pure Condition Test ( Xi = 1 )
|
|
Switch pure = false;
|
|
label pureSpecieI = -1;
|
|
forAll (X, i)
|
|
{
|
|
if (1. - X[i] < 1e-12)
|
|
{
|
|
pureSpecieI = i;
|
|
pure = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Mixture average D, mu and k
|
|
if (pure)
|
|
{
|
|
const scalar* Dpure = Dij[pureSpecieI];
|
|
|
|
forAll (X, i)
|
|
{
|
|
Di[i] = Dpure[i];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
forAll (X, k)
|
|
{
|
|
scalar sum1 = 0.0;
|
|
scalar sum2 = 0.0;
|
|
|
|
forAll (X, i)
|
|
{
|
|
if (i == k)
|
|
{
|
|
continue;
|
|
}
|
|
|
|
const scalar* Dk = Dij[k];
|
|
|
|
sum1 += X[i] / Dk[i];
|
|
sum2 += Y[i] / Dk[i];
|
|
}
|
|
|
|
sum2 *= X[k] / (1. - Y[k]);
|
|
|
|
Di[k] = inv(sum1 + sum2);
|
|
}
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::diffusivityModel::mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X)
|
|
{
|
|
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 = thermoData()[j].W();
|
|
const scalar Wk = thermoData()[k].W();
|
|
|
|
const scalar muj = muI[j];
|
|
const scalar muk = muI[k];
|
|
|
|
Phi(k,j) = a * sqr(1. + sqrt(muk/muj)*pow(Wj/Wk, 1./4.)) / sqrt(1. + Wk/Wj);
|
|
}
|
|
}
|
|
|
|
|
|
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;
|
|
}
|
|
|
|
|
|
void Foam::diffusivityModel::calculateMuD
|
|
(
|
|
UList<scalar> &muI,
|
|
scalarSymmetricSquareMatrix &Dij,
|
|
const GasState &state
|
|
)
|
|
{
|
|
const scalar pi = state.p();
|
|
const scalar Ti = state.T();
|
|
const scalar rhoQc2i = state.rhoQc2();
|
|
const scalarField &localY(state.Y());
|
|
const scalarField &localX(state.X());
|
|
|
|
label ionStart = 0;
|
|
|
|
label neutralStart = 0;
|
|
|
|
if (electron_.valid())
|
|
{
|
|
ionStart = 1;
|
|
|
|
neutralStart = ionStart + ions_.size();
|
|
|
|
// Electron - Electron
|
|
label idx = 0;
|
|
|
|
muI[0] = ecs_[idx].mu(pi,Ti,rhoQc2i);
|
|
|
|
Dij(0,0) = ecs_[idx].D(pi,Ti,rhoQc2i);
|
|
|
|
// Electron - Ions
|
|
idx = 1;
|
|
forAll (ions_, J)
|
|
{
|
|
label j = J + ionStart;
|
|
|
|
// Calculate Dij
|
|
Dij(0,j) = ecs_[idx].D(pi,Ti,rhoQc2i);
|
|
Dij(j,0) = Dij(0,j);
|
|
idx++;
|
|
}
|
|
|
|
// Electron - Neutrals
|
|
idx = 0;
|
|
forAll (neutrals_, J)
|
|
{
|
|
label j = J + neutralStart;
|
|
|
|
// Calculate Dij
|
|
Dij(0,j) = ens_[idx].D(pi,Ti);
|
|
Dij(j,0) = Dij(0,j);
|
|
idx++;
|
|
}
|
|
|
|
}
|
|
|
|
// Ions
|
|
label idx1 = 0;
|
|
forAll (ions_, I)
|
|
{
|
|
label i = I + ionStart;
|
|
|
|
muI[i] = ccs_[idx1].mu(pi,Ti,rhoQc2i);
|
|
|
|
Dij(i,i) = ccs_[idx1].D(pi,Ti,rhoQc2i);
|
|
|
|
idx1++;
|
|
|
|
// Ion - Ions
|
|
for (label J = I + 1; J < ions_.size(); J++)
|
|
{
|
|
label j = J + ionStart;
|
|
// Calculate Dij
|
|
Dij(i,j) = ccs_[idx1].D(pi,Ti,rhoQc2i);
|
|
Dij(j,i) = Dij(i,j);
|
|
idx1++;
|
|
}
|
|
|
|
}
|
|
|
|
label idx2 = 0;
|
|
forAll (ions_, I)
|
|
{
|
|
label i = I + ionStart;
|
|
|
|
// Ion - Neutrals
|
|
forAll (neutrals_, J)
|
|
{
|
|
label j = J + neutralStart;
|
|
|
|
// Calculate Dij
|
|
Dij(i,j) = cns_[idx2].D(pi,Ti);
|
|
Dij(j,i) = Dij(i,j);
|
|
idx2++;
|
|
}
|
|
}
|
|
|
|
// Neutrals
|
|
label idx = 0;
|
|
forAll (neutrals_, I)
|
|
{
|
|
label i = I + neutralStart;
|
|
|
|
muI[i] = nns_[idx].mu(pi,Ti);
|
|
Dij(i,i) = nns_[idx].D(pi,Ti);
|
|
idx++;
|
|
|
|
// Neutral - Neutrals
|
|
for (label J = I + 1; J < neutrals_.size(); J++)
|
|
{
|
|
label j = J + neutralStart;
|
|
|
|
// Calculate Dij
|
|
Dij(i,j) = nns_[idx].D(pi,Ti);
|
|
Dij(j,i) = Dij(i,j);
|
|
idx++;
|
|
}
|
|
}
|
|
|
|
|
|
}
|
|
|
|
|
|
void Foam::diffusivityModel::calculateK
|
|
(
|
|
UList<scalar> &kI,
|
|
const UList<scalar> &muI,
|
|
const UList<scalar> &Dii,
|
|
const GasState &state
|
|
)
|
|
{
|
|
label ionStart = 0;
|
|
|
|
if (electron_.valid())
|
|
{
|
|
ionStart = 1;
|
|
|
|
kI[0] = electron_->lambda(muI[0], Dii[0], state);
|
|
}
|
|
|
|
forAll (ions_, j)
|
|
{
|
|
label i = j + ionStart;
|
|
|
|
kI[i] = ions_[j].lambda(muI[i], Dii[i], state);
|
|
}
|
|
|
|
forAll (neutrals_, j)
|
|
{
|
|
label i = j + ionStart + ions_.size();
|
|
|
|
kI[i] = neutrals_[j].lambda(muI[i], Dii[i], state);
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
|
|
:
|
|
thermo_(thermo),
|
|
D_(thermo_.composition().species().size()),
|
|
mu_
|
|
(
|
|
IOobject
|
|
(
|
|
"mu",
|
|
thermo_.composition().Y(0).mesh().time().timeName(),
|
|
thermo_.composition().Y(0).mesh(),
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
thermo_.composition().Y(0).mesh(),
|
|
dimensionedScalar("zero", dimDynamicViscosity, 0.0)
|
|
),
|
|
k_
|
|
(
|
|
IOobject
|
|
(
|
|
"lambda",
|
|
thermo_.composition().Y(0).mesh().time().timeName(),
|
|
thermo_.composition().Y(0).mesh(),
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
thermo_.composition().Y(0).mesh(),
|
|
dimensionedScalar("zero", dimForce/dimTime/dimTemperature, 0.0)
|
|
),
|
|
ions_(thermo_.composition().species().size()),
|
|
neutrals_(thermo_.composition().species().size()),
|
|
ecs_(thermo_.composition().species().size()),
|
|
ens_(thermo_.composition().species().size()),
|
|
nns_(thermo_.composition().species().size()*thermo_.composition().species().size()),
|
|
ccs_(thermo_.composition().species().size()*thermo_.composition().species().size()),
|
|
cns_(thermo_.composition().species().size()*thermo_.composition().species().size())
|
|
{
|
|
|
|
typedef multiComponentMixture<gasHThermoPhysics> gasMix;
|
|
|
|
const gasMix &mcm (dynamicCast<const gasMix>(thermo_));
|
|
|
|
const speciesTable &species_(thermo_.composition().species());
|
|
|
|
const volScalarField::Mesh &mesh = thermo_.composition().Y(0).mesh();
|
|
|
|
if (thermoData_ == NULL)
|
|
{
|
|
thermoData_ = &(mcm.speciesData());
|
|
}
|
|
|
|
dictionary transports
|
|
(
|
|
IFstream
|
|
(
|
|
fileName(thermo_.lookup("foamTransportFile")).expand()
|
|
)()
|
|
);
|
|
|
|
dictionary thermos
|
|
(
|
|
IFstream
|
|
(
|
|
fileName(thermo_.lookup("foamChemistryThermoFile")).expand()
|
|
)()
|
|
);
|
|
|
|
Switch readIons = false;
|
|
Switch readNeutrals = false;
|
|
|
|
label nIon = 0;
|
|
label nNeutral = 0;
|
|
|
|
forAll(species_, i)
|
|
{
|
|
const word namei(species_[i]);
|
|
const scalar Wi = thermo_.composition().W(i);
|
|
const label zi = thermo_.composition().z(i);
|
|
|
|
dictionary thermoDict(thermos.subDict(species_[i]));
|
|
|
|
dictionary tranDict(transports.subDict(species_[i]).subDict("transport"));
|
|
tranDict.add("name", namei);
|
|
tranDict.add("W", Wi);
|
|
tranDict.add("z", zi);
|
|
|
|
const Particle::thermoType &thermoData(mcm.speciesData()[i]);
|
|
|
|
// Electron
|
|
if (species_[i] == "E-")
|
|
{
|
|
if (readIons || readNeutrals)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Electron should be the first species in the list"
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
electron_.set(new Electron (thermoDict, 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(thermoDict, tranDict));
|
|
|
|
nIon++;
|
|
|
|
// new Ion (tranDict);
|
|
}
|
|
// Neutrals
|
|
else
|
|
{
|
|
readNeutrals = true;
|
|
|
|
neutrals_.set(nNeutral, new Neutral(thermoDict, tranDict));
|
|
|
|
nNeutral++;
|
|
}
|
|
}
|
|
|
|
ions_.resize(nIon);
|
|
|
|
neutrals_.resize(nNeutral);
|
|
nns_.resize(nNeutral*(nNeutral+1)/2);
|
|
|
|
ccs_.resize((nIon+1)*(nIon)/2);
|
|
cns_.resize(nNeutral*nIon);
|
|
ecs_.resize(nIon+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;
|
|
|
|
if (electron_.valid())
|
|
{
|
|
ecs_.set(iPair, new Coulomb(electron_(), electron_()));
|
|
iPair++;
|
|
|
|
forAll(ions_, i)
|
|
{
|
|
ecs_.set(iPair, new Coulomb(electron_(), ions_[i]));
|
|
iPair++;
|
|
}
|
|
|
|
iPair = 0;
|
|
forAll(neutrals_, i)
|
|
{
|
|
ens_.set(iPair, new CrossSection(electron_(), neutrals_[i]));
|
|
iPair++;
|
|
}
|
|
}
|
|
|
|
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++;
|
|
}
|
|
}
|
|
|
|
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_),
|
|
mu_(dm.mu_),
|
|
k_(dm.k_)
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * 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 speciesTable &species_(thermo_.composition().species());
|
|
|
|
const PtrList<volScalarField> &Y(thermo_.composition().Y());
|
|
|
|
scalarSymmetricSquareMatrix Dij(species_.size());
|
|
|
|
scalarField localY(species_.size());
|
|
|
|
scalarField Di(species_.size());
|
|
scalarField Dii(species_.size());
|
|
scalarField muI(species_.size());
|
|
scalarField kI(species_.size());
|
|
|
|
forAll (p, celli)
|
|
{
|
|
const scalar pi = p[celli];
|
|
const scalar Ti = T[celli];
|
|
|
|
forAll (species_, i)
|
|
{
|
|
localY[i] = Y[i][celli];
|
|
}
|
|
|
|
GasState state (pi, Ti, localY);
|
|
|
|
const scalarField &localX = state.X();
|
|
|
|
calculateMuD (muI, Dij, state);
|
|
|
|
forAll (Dii, i)
|
|
{
|
|
Dii[i] = Dij(i,i);
|
|
}
|
|
|
|
// Pure Thermal conductivities
|
|
calculateK ( kI, muI, Dii, state);
|
|
|
|
mixAvgDi(Di, Dij, localX, localY);
|
|
forAll (Di, i)
|
|
{
|
|
D_[i][celli] = Di[i];
|
|
}
|
|
|
|
mu_[celli] = mixAvgMu(muI, localX);
|
|
k_[celli] = mixAvgK(kI, localX);
|
|
}
|
|
|
|
|
|
forAll(p.boundaryField(), patchi)
|
|
{
|
|
const volScalarField::Patch &pp = p.boundaryField()[patchi];
|
|
const volScalarField::Patch &Tp = T.boundaryField()[patchi];
|
|
|
|
forAll(pp, facei)
|
|
{
|
|
const scalar pi = pp[facei];
|
|
const scalar Ti = Tp[facei];
|
|
|
|
forAll (species_, i)
|
|
{
|
|
localY[i] = Y[i].boundaryField()[patchi][facei];
|
|
}
|
|
|
|
GasState state (pi, Ti, localY);
|
|
|
|
const scalarField &localX = state.X();
|
|
|
|
calculateMuD (muI, Dij, state);
|
|
|
|
forAll (Dii, i)
|
|
{
|
|
Dii[i] = Dij(i,i);
|
|
}
|
|
|
|
|
|
// Pure Thermal conductivities
|
|
calculateK ( kI, muI, Dii, state);
|
|
|
|
mixAvgDi(Di, Dij, localX, localY);
|
|
forAll (Di, i)
|
|
{
|
|
D_[i].boundaryFieldRef()[patchi][facei] = Di[i];
|
|
}
|
|
|
|
mu_.boundaryFieldRef()[patchi][facei] = mixAvgMu(muI, localX);
|
|
k_.boundaryFieldRef()[patchi][facei] = mixAvgK(kI, localX);
|
|
}
|
|
}
|
|
|
|
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 * * * * * * * * * * * * * * //
|
|
|
|
|
|
// ************************************************************************* //
|