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

calculating D, mu, k of Neutral mixture

Browse source
This commit is contained in:
ignis 2018-04-30 04:02:16 +09:00
parent 7baaa4f389
commit 80c5f0d65a
8 changed files with 212 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -2,3 +2,6 @@
[._]*.sw[a-p]
[._]s[a-v][a-z]
[._]sw[a-p]
Make
!Make/file
!Make/option

2
diffusivityModel/Neutral/Neutral.H
View file

@ -136,6 +136,8 @@ public:
inline scalar Zrot() const;
inline scalar Zrot(scalar T) const;
// Check
// Edit

14
diffusivityModel/Neutral/NeutralI.H
View file

@ -84,6 +84,20 @@ inline Foam::scalar Foam::Neutral::Zrot()
return Zrot_;
}
inline Foam::scalar Foam::Neutral::Zrot(const scalar T)
const
{
scalar a3 = sqrt(pow3(pi));
scalar a1 = a3/2.0;
scalar a2 = sqr(pi) / 4.0 + 2.0;
scalar t1 = sqrt(wellDepth_/T);
scalar t2 = wellDepth_/T;
scalar t3 = t1 * t2;
return 1.0 + a1 * t1 + a2 * t2 + a3 * t3;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //

6
diffusivityModel/Particle/Particle.H
View file

@ -68,8 +68,8 @@ public:
enum Geometry
{
ATOM = 0,
LINEAR = 1,
NONLINEAR = 2
LINEAR = 2,
NONLINEAR = 3
};
@ -150,6 +150,8 @@ public:
inline scalar z() const;
inline scalar F() const;
inline Geometry geometry() const;
// Check

7
diffusivityModel/Particle/ParticleI.H
View file

@ -57,6 +57,13 @@ inline Foam::scalar Foam::Particle::z()
}
inline Foam::scalar Foam::Particle::F()
const
{
return (geometry_);
}
inline Foam::Particle::Geometry Foam::Particle::geometry()
const
{

180
diffusivityModel/diffusivityModel/diffusivityModel.C
View file

@ -43,6 +43,62 @@ License
// * * * * * * * * * * * * * 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 * * * * * * * * * * * //
@ -190,6 +246,40 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
}
}
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
@ -240,12 +330,10 @@ void Foam::diffusivityModel::correct()
{
const volScalarField &T = thermo_.T();
const volScalarField &p = thermo_.p();
const volScalarField rho(thermo_.rho());
const speciesTable &species_(thermo_.composition().species());
scalarSymmetricSquareMatrix Dij(species_.size());
const volScalarField Cv(thermo_.Cv());
const volScalarField Wbar(thermo_.composition().W());
const PtrList<volScalarField> &Y(thermo_.composition().Y());
@ -256,17 +344,29 @@ void Foam::diffusivityModel::correct()
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());
forAll (p, celli)
{
const scalar rhoi = rho[celli];
const scalar pi = p[celli];
const scalar Ti = T[celli];
const scalar Cvi = Cv[celli];
const scalar WbarI = Wbar[celli];
scalarField localX(species_.size());
scalarField localY(species_.size());
scalarField localD(species_.size());
forAll (species_, i)
{
localCv[i] = thermo_.composition().Cv(i,pi,Ti);
}
forAll (species_, i)
{
@ -278,25 +378,75 @@ void Foam::diffusivityModel::correct()
forAll (species_, i)
{
for (label j = i; j < species_.size(); j++)
muI[i] = nns_[idx].mu(pi,Ti);
Dii[i] = nns_[idx].D(pi,Ti);
Dij(i,i) = Dii[i];
idx++;
for (label j = i+1; j < species_.size(); j++)
{
// Calculate Dij
// Dij[idx] = interactions[idx].D(pi,Ti);
Dij(i,j) = nns_[idx].D(pi,Ti);
Dij(j,i) = Dij(i,j);
idx++;
}
// DMat = 0;
//
}
forAll (species_, i)
{
// localD[i] = mixAvgD(Dij.row(i));
const scalar R = Neutral::RR / Wpure[i];
const scalar CvTrans = (3./2.)*R;
const scalar CvRot = (neutrals_[i].F()/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 Zrot = neutrals_[i].Zrot(Ti);
const scalar B = Zrot + (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);
}
Info << localY << endl;
Switch pure = false;
label pureSpecieI = -1;
forAll (species_, i)
{
if (1. - localX[i] < 1e-12)
{
pureSpecieI = i;
pure = true;
break;
}
}
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);
}
return;

15
diffusivityModel/diffusivityModel/diffusivityModel.H
View file

@ -67,6 +67,12 @@ class diffusivityModel
//- Mixture averaged diffusivity fields
PtrList<volScalarField> D_;
//- Mixture averaged viscosity field
autoPtr<volScalarField> mu_;
//- Mixture averaged thermal conductivity field
autoPtr<volScalarField> k_;
//- Electron object
// autoPtr<Electron> electron_;
@ -88,6 +94,15 @@ class diffusivityModel
//- Disallow default bitwise assignment
void operator=(const diffusivityModel&);
//- Disallow default bitwise assignment
scalar mixAvgD(const label k, const UList<scalar>& Di, const UList<scalar>& X, const UList<scalar>& Y);
//- Disallow default bitwise assignment
scalar mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W);
//- Disallow default bitwise assignment
scalar mixAvgK(const UList<scalar>& k, const UList<scalar>& X);
public:

2
eReactingFoam.C
View file

@ -91,6 +91,8 @@ int main(int argc, char *argv[])
while (pimple.loop())
{
diff.correct();
#include "UEqn.H"
#include "YEqn.H"
#include "EEqn.H"