D, mu and k boundary fields update in correct()
This commit is contained in:
ignis
parent
6033bedc3a
commit
c5027f53ba
2 changed files with 286 additions and 133 deletions
|
|
@ -57,11 +57,68 @@ Foam::scalar Foam::diffusivityModel::mixAvgD(const label k, const UList<scalar>&
|
|||
|
||||
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 UList<scalar>& W)
|
||||
{
|
||||
const label K = muI.size();
|
||||
|
|
@ -105,6 +162,122 @@ Foam::scalar Foam::diffusivityModel::mixAvgK(const UList<scalar>& k, const UList
|
|||
}
|
||||
|
||||
|
||||
void Foam::diffusivityModel::calculateMuD
|
||||
(
|
||||
UList<scalar> &muI,
|
||||
scalarSymmetricSquareMatrix &Dij,
|
||||
const scalar rhoi,
|
||||
const scalar Cpi,
|
||||
const scalar pi,
|
||||
const scalar Ti,
|
||||
const scalar WbarI,
|
||||
const scalar rhoQc2i,
|
||||
const scalarField &localY,
|
||||
const scalarField &localX
|
||||
)
|
||||
{
|
||||
|
||||
const label ionStart = 1;
|
||||
|
||||
const label neutralStart = 1 + 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
|
||||
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++;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
||||
|
||||
|
||||
|
|
@ -388,6 +561,7 @@ Foam::diffusivityModel::~diffusivityModel()
|
|||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
void Foam::diffusivityModel::correct()
|
||||
{
|
||||
const volScalarField &T = thermo_.T();
|
||||
|
|
@ -417,6 +591,7 @@ void Foam::diffusivityModel::correct()
|
|||
|
||||
scalarField localD(species_.size());
|
||||
|
||||
scalarField Di(species_.size());
|
||||
scalarField Dii(species_.size());
|
||||
scalarField muI(species_.size());
|
||||
scalarField kI(species_.size());
|
||||
|
|
@ -449,10 +624,7 @@ void Foam::diffusivityModel::correct()
|
|||
const scalar WbarI = Wbar[celli];
|
||||
const scalar rhoQc2i = rhoQc2[celli];
|
||||
|
||||
forAll (species_, i)
|
||||
{
|
||||
localCv[i] = thermo_.composition().Cv(i,pi,Ti);
|
||||
}
|
||||
thermo_.composition().Cv(localCv,pi,Ti);
|
||||
|
||||
forAll (species_, i)
|
||||
{
|
||||
|
|
@ -460,111 +632,24 @@ void Foam::diffusivityModel::correct()
|
|||
localX[i] = Y[i][celli] * WbarI / Wpure[i];
|
||||
}
|
||||
|
||||
calculateMuD ( muI, Dij, rhoi, Cpi, pi, Ti, WbarI, rhoQc2i, localY, localX);
|
||||
|
||||
forAll (Dii, i)
|
||||
{
|
||||
Dii[i] = Dij(i,i);
|
||||
}
|
||||
|
||||
|
||||
Zrot[0] = 0;
|
||||
forAll (ions_, i)
|
||||
{
|
||||
// ZrotI[i] = ions_[i].Zrot(Ti);
|
||||
ZrotI[i] = 0.0;
|
||||
ZrotI[i] = ions_[i].Zrot(Ti);
|
||||
}
|
||||
forAll (neutrals_, i)
|
||||
{
|
||||
ZrotN[i] = neutrals_[i].Zrot(Ti);
|
||||
}
|
||||
|
||||
// Electron - Electron
|
||||
label idx = 0;
|
||||
|
||||
muI[0] = ecs_[idx].mu(pi,Ti,rhoQc2i);
|
||||
// muI[0] = SMALL;
|
||||
|
||||
Dii[0] = ecs_[idx].D(pi,Ti,rhoQc2i);
|
||||
// Dii[0] = SMALL;
|
||||
Dij(0,0) = Dii[0];
|
||||
|
||||
// Electron - Ions
|
||||
idx = 1;
|
||||
for (label J = 0; J < ions_.size(); J++)
|
||||
{
|
||||
label j = J + 1;
|
||||
// Calculate Dij
|
||||
Dij(0,j) = ecs_[idx].D(pi,Ti,rhoQc2i);
|
||||
// Dij(0,j) = SMALL;
|
||||
Dij(j,0) = Dij(0,j);
|
||||
idx++;
|
||||
}
|
||||
|
||||
// Electron - Neutrals
|
||||
idx = 0;
|
||||
for (label J = 0; J < neutrals_.size(); J++)
|
||||
{
|
||||
label j = J + ions_.size() + 1;
|
||||
// Calculate Dij
|
||||
Dij(0,j) = ens_[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) = ccs_[idx1].D(pi,Ti,rhoQc2i);
|
||||
// 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) = cns_[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)
|
||||
{
|
||||
|
|
@ -587,42 +672,94 @@ void Foam::diffusivityModel::correct()
|
|||
|
||||
}
|
||||
|
||||
// Pure Condition Test ( Xi = 1 )
|
||||
Switch pure = false;
|
||||
label pureSpecieI = -1;
|
||||
forAll (species_, i)
|
||||
mixAvgDi(Di, Dij, localX, localY);
|
||||
forAll (Di, 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);
|
||||
}
|
||||
D_[i][celli] = Di[i];
|
||||
}
|
||||
|
||||
mu_()[celli] = mixAvgMu(muI, localX, Wpure);
|
||||
k_()[celli] = mixAvgK(kI, localX) / rhoi / Cpi;
|
||||
}
|
||||
|
||||
|
||||
forAll(p.boundaryField(), patchi)
|
||||
{
|
||||
const volScalarField::Patch &rhop = rho.boundaryField()[patchi];
|
||||
const volScalarField::Patch &Cpp = Cp.boundaryField()[patchi];
|
||||
const volScalarField::Patch &pp = p.boundaryField()[patchi];
|
||||
const volScalarField::Patch &Tp = T.boundaryField()[patchi];
|
||||
const volScalarField::Patch &Wbarp = Wbar.boundaryField()[patchi];
|
||||
const volScalarField::Patch &rhoQc2p = rhoQc2.boundaryField()[patchi];
|
||||
|
||||
forAll(rhop, facei)
|
||||
{
|
||||
const scalar rhoi = rhop[facei];
|
||||
const scalar Cpi = Cpp[facei];
|
||||
const scalar pi = pp[facei];
|
||||
const scalar Ti = Tp[facei];
|
||||
const scalar WbarI = Wbarp[facei];
|
||||
const scalar rhoQc2i = rhoQc2p[facei];
|
||||
|
||||
thermo_.composition().Cv(localCv,pi,Ti);
|
||||
|
||||
forAll (species_, i)
|
||||
{
|
||||
localY[i] = Y[i].boundaryField()[patchi][facei];
|
||||
localX[i] = localY[i] * WbarI / Wpure[i];
|
||||
}
|
||||
|
||||
calculateMuD ( muI, Dij, rhoi, Cpi, pi, Ti, WbarI, rhoQc2i, localY, localX);
|
||||
|
||||
forAll (Dii, i)
|
||||
{
|
||||
Dii[i] = Dij(i,i);
|
||||
}
|
||||
|
||||
|
||||
Zrot[0] = 0;
|
||||
forAll (ions_, i)
|
||||
{
|
||||
ZrotI[i] = ions_[i].Zrot(Ti);
|
||||
}
|
||||
forAll (neutrals_, i)
|
||||
{
|
||||
ZrotN[i] = neutrals_[i].Zrot(Ti);
|
||||
}
|
||||
|
||||
// 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);
|
||||
|
||||
}
|
||||
|
||||
mixAvgDi(Di, Dij, localX, localY);
|
||||
forAll (Di, i)
|
||||
{
|
||||
D_[i].boundaryFieldRef()[patchi][facei] = Di[i];
|
||||
}
|
||||
|
||||
mu_().boundaryFieldRef()[patchi][facei] = mixAvgMu(muI, localX, Wpure);
|
||||
k_().boundaryFieldRef()[patchi][facei] = mixAvgK(kI, localX) / rhoi / Cpi;
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -112,14 +112,30 @@ class diffusivityModel
|
|||
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);
|
||||
inline 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);
|
||||
inline void mixAvgDi(UList<scalar>& Di, const scalarSymmetricSquareMatrix &Dij, const UList<scalar>& X, const UList<scalar>& Y);
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
scalar mixAvgK(const UList<scalar>& k, const UList<scalar>& X);
|
||||
inline scalar mixAvgMu(const UList<scalar>& muI, const UList<scalar>& X, const UList<scalar>& W);
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
inline scalar mixAvgK(const UList<scalar>& k, const UList<scalar>& X);
|
||||
|
||||
inline void calculateMuD
|
||||
(
|
||||
UList<scalar> &muI,
|
||||
scalarSymmetricSquareMatrix &Dij,
|
||||
const scalar rhoi,
|
||||
const scalar Cpi,
|
||||
const scalar pi,
|
||||
const scalar Ti,
|
||||
const scalar WbarI,
|
||||
const scalar rhoQc2i,
|
||||
const scalarField &localY,
|
||||
const scalarField &localX
|
||||
);
|
||||
|
||||
public:
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue