diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C index 22e47a654e..6474c26346 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -45,6 +45,14 @@ Foam::LUscalarMatrix::LUscalarMatrix() {} +Foam::LUscalarMatrix::LUscalarMatrix(const label n) +: + scalarSquareMatrix(n), + comm_(Pstream::worldComm), + pivotIndices_(n) +{} + + Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix) : scalarSquareMatrix(matrix), @@ -404,6 +412,13 @@ void Foam::LUscalarMatrix::printDiagonalDominance() const } +void Foam::LUscalarMatrix::decompose() +{ + pivotIndices_.setSize(m()); + LUDecompose(*this, pivotIndices_); +} + + void Foam::LUscalarMatrix::decompose(const scalarSquareMatrix& M) { scalarSquareMatrix::operator=(M); diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H index 7ded91b5bf..9226e849fe 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -99,6 +99,9 @@ public: //- Construct null LUscalarMatrix(); + //- Construct given number of rows/columns + LUscalarMatrix(const label n); + //- Construct from and perform LU decomposition of the matrix M LUscalarMatrix(const scalarSquareMatrix& M); @@ -113,6 +116,9 @@ public: // Member Functions + //- Perform the LU decomposition of the matrix + void decompose(); + //- Perform the LU decomposition of the matrix M void decompose(const scalarSquareMatrix& M); diff --git a/src/ThermophysicalTransportModels/fluidReactionThermo/fluidReactionThermophysicalTransportModels.C b/src/ThermophysicalTransportModels/fluidReactionThermo/fluidReactionThermophysicalTransportModels.C index d5b3416ada..c82f62bace 100644 --- a/src/ThermophysicalTransportModels/fluidReactionThermo/fluidReactionThermophysicalTransportModels.C +++ b/src/ThermophysicalTransportModels/fluidReactionThermo/fluidReactionThermophysicalTransportModels.C @@ -48,6 +48,9 @@ makeLaminarThermophysicalTransportModel(unityLewisFourier); #include "FickianFourier.H" makeLaminarThermophysicalTransportModel(FickianFourier); +#include "MaxwellStefanFourier.H" +makeLaminarThermophysicalTransportModel(MaxwellStefanFourier); + // -------------------------------------------------------------------------- // // RAS models diff --git a/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.C b/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.C index 026ce20ea7..268e5b49fd 100644 --- a/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.C +++ b/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.C @@ -39,8 +39,7 @@ namespace Foam // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Fickian:: -Fickian +Fickian::Fickian ( const word& type, const momentumTransportModel& momentumTransport, @@ -55,6 +54,7 @@ Fickian ), D_(this->thermo().composition().species().size()), + DT_ ( this->coeffDict_.found("DT") @@ -67,8 +67,7 @@ Fickian // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -bool -Fickian::read() +bool Fickian::read() { if ( @@ -103,8 +102,7 @@ Fickian::read() template -tmp -Fickian::DEff +tmp Fickian::DEff ( const volScalarField& Yi ) const @@ -128,8 +126,7 @@ Fickian::DEff template -tmp -Fickian::DEff +tmp Fickian::DEff ( const volScalarField& Yi, const label patchi @@ -149,8 +146,7 @@ Fickian::DEff template -tmp -Fickian::q() const +tmp Fickian::q() const { tmp tmpq ( @@ -226,8 +222,10 @@ Fickian::q() const template -tmp -Fickian::divq(volScalarField& he) const +tmp Fickian::divq +( + volScalarField& he +) const { tmp tmpDivq ( @@ -322,8 +320,7 @@ Fickian::divq(volScalarField& he) const template -tmp -Fickian::j +tmp Fickian::j ( const volScalarField& Yi ) const @@ -350,8 +347,10 @@ Fickian::j template -tmp -Fickian::divj(volScalarField& Yi) const +tmp Fickian::divj +( + volScalarField& Yi +) const { if (DT_.size()) { diff --git a/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.H b/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.H index bbd6e1cf71..326aa080db 100644 --- a/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.H +++ b/src/ThermophysicalTransportModels/laminar/Fickian/Fickian.H @@ -28,8 +28,8 @@ Description Base class for multi-component Fickian based temperature gradient heat flux models with optional Soret thermal diffusion of species. - Currently the diffusion coefficients are constant but temperature and - pressure dependency will be added. + The mixture diffusion coefficients are specified as Function2s of + pressure and temperature but independent of composition. The heat flux source is implemented as an implicit energy correction to the temperature gradient based flux source. At convergence the energy diff --git a/src/ThermophysicalTransportModels/laminar/FickianFourier/FickianFourier.H b/src/ThermophysicalTransportModels/laminar/FickianFourier/FickianFourier.H index ed0655bcc7..2d329ebaf8 100644 --- a/src/ThermophysicalTransportModels/laminar/FickianFourier/FickianFourier.H +++ b/src/ThermophysicalTransportModels/laminar/FickianFourier/FickianFourier.H @@ -26,10 +26,10 @@ Class Description Multi-component Fickian and Fourier based temperature gradient heat flux - model laminar flow with optional Soret thermal diffusion of species. + models with optional Soret thermal diffusion of species for laminar flow. - Currently the diffusion coefficients are constant but temperature and - pressure dependency will be added. + The mixture diffusion coefficients are specified as Function2s of + pressure and temperature but independent of composition. The heat flux source is implemented as an implicit energy correction to the temperature gradient based flux source. At convergence the energy @@ -41,8 +41,19 @@ Usage { model FickianFourier; - D (1e-5 2e-5 1e-5); // [m^2/s] - DT (1e-5 2e-5 1e-5); // [kg/m/s] Optional + D // [m^2/s] + { + O2 1e-2; + O3 5e-2; + N2 1e-2; + } + + DT // [kg/m/s] Optional + { + O2 1e-2; + O3 5e-2; + N2 1e-2; + } } \endverbatim diff --git a/src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.C b/src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.C new file mode 100644 index 0000000000..4358752d22 --- /dev/null +++ b/src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.C @@ -0,0 +1,697 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2021 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 . + +\*---------------------------------------------------------------------------*/ + +#include "MaxwellStefan.H" +#include "fvcDiv.H" +#include "fvcLaplacian.H" +#include "fvcSnGrad.H" +#include "fvmSup.H" +#include "surfaceInterpolate.H" +#include "Function2Evaluate.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +MaxwellStefan::MaxwellStefan +( + const word& type, + const momentumTransportModel& momentumTransport, + const thermoModel& thermo +) +: + BasicThermophysicalTransportModel + ( + type, + momentumTransport, + thermo + ), + + D_(this->thermo().composition().species().size()), + + DT_ + ( + this->coeffDict_.found("DT") + ? this->thermo().composition().species().size() + : 0 + ), + + Dii_(this->thermo().composition().species().size()), + jexp_(this->thermo().composition().species().size()), + + W(this->thermo().composition().species().size()), + + YPtrs(W.size()), + DijPtrs(W.size()), + + Y(W.size()), + X(W.size()), + DD(W.size()), + A(W.size() - 1), + B(A.m()), + invA(A.m()), + D(W.size()) +{ + const basicSpecieMixture& composition = this->thermo().composition(); + + // Set the molecular weights of the species + forAll(W, i) + { + W[i] = composition.Wi(i); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool MaxwellStefan::read() +{ + if + ( + BasicThermophysicalTransportModel::read() + ) + { + const basicSpecieMixture& composition = this->thermo().composition(); + const speciesTable& species = composition.species(); + + const dictionary& Ddict = this->coeffDict_.subDict("D"); + + // Read the array of specie binary mass diffusion coefficient functions + forAll(species, i) + { + D_[i].setSize(species.size()); + + forAll(species, j) + { + if (j >= i) + { + const word nameij(species[i] + '_' + species[j]); + const word nameji(species[j] + '_' + species[i]); + + word Dname; + + if (Ddict.found(nameij) && Ddict.found(nameji)) + { + if (i != j) + { + WarningInFunction + << "Binary diffusivities for Both " + << nameij << " and " << nameji + << " provided, using " << nameij << endl; + } + + Dname = nameij; + } + else if (Ddict.found(nameij)) + { + Dname = nameij; + } + else if (Ddict.found(nameji)) + { + Dname = nameji; + } + else + { + FatalIOErrorInFunction(Ddict) + << "Binary diffusivity for pair " << nameij + << " or " << nameji << " not provided" + << exit(FatalIOError); + } + + D_[i].set + ( + j, + Function2::New(Dname, Ddict).ptr() + ); + } + } + } + + // Optionally read the List of specie Soret thermal diffusion + // coefficient functions + if (this->coeffDict_.found("DT")) + { + const dictionary& DTdict = this->coeffDict_.subDict("DT"); + + forAll(species, i) + { + DT_.set(i, Function2::New(species[i], DTdict).ptr()); + } + } + + return true; + } + else + { + return false; + } +} + + +template +tmp MaxwellStefan::DEff +( + const volScalarField& Yi +) const +{ + const basicSpecieMixture& composition = this->thermo().composition(); + + return volScalarField::New + ( + "DEff", + this->momentumTransport().rho()*Dii_[composition.index(Yi)] + ); +} + + +template +tmp MaxwellStefan::DEff +( + const volScalarField& Yi, + const label patchi +) const +{ + const basicSpecieMixture& composition = this->thermo().composition(); + + return + this->momentumTransport().rho().boundaryField()[patchi] + *Dii_[composition.index(Yi)].boundaryField()[patchi]; +} + + +template +tmp +MaxwellStefan::q() const +{ + tmp tmpq + ( + surfaceScalarField::New + ( + IOobject::groupName + ( + "q", + this->momentumTransport().alphaRhoPhi().group() + ), + -fvc::interpolate(this->alpha()*this->kappaEff()) + *fvc::snGrad(this->thermo().T()) + ) + ); + + const basicSpecieMixture& composition = this->thermo().composition(); + const label d = composition.defaultSpecie(); + + const PtrList& Y = composition.Y(); + + if (Y.size()) + { + surfaceScalarField sumJ + ( + surfaceScalarField::New + ( + "sumJ", + Y[0].mesh(), + dimensionedScalar(dimMass/dimArea/dimTime, 0) + ) + ); + + surfaceScalarField sumJh + ( + surfaceScalarField::New + ( + "sumJh", + Y[0].mesh(), + dimensionedScalar(sumJ.dimensions()*dimEnergy/dimMass, 0) + ) + ); + + forAll(Y, i) + { + if (i != d) + { + const volScalarField hi + ( + composition.HE(i, this->thermo().p(), this->thermo().T()) + ); + + const surfaceScalarField ji(this->j(Y[i])); + sumJ += ji; + + sumJh += ji*fvc::interpolate(hi); + } + } + + { + const label i = d; + + const volScalarField hi + ( + composition.HE(i, this->thermo().p(), this->thermo().T()) + ); + + sumJh -= sumJ*fvc::interpolate(hi); + } + + tmpq.ref() += sumJh; + } + + return tmpq; +} + + +template +tmp MaxwellStefan::divq +( + volScalarField& he +) const +{ + tmp tmpDivq + ( + fvm::Su + ( + -fvc::laplacian(this->alpha()*this->kappaEff(), this->thermo().T()), + he + ) + ); + + const basicSpecieMixture& composition = this->thermo().composition(); + const label d = composition.defaultSpecie(); + + const PtrList& Y = composition.Y(); + + if (!Y.size()) + { + tmpDivq.ref() -= + correction(fvm::laplacian(this->alpha()*this->alphaEff(), he)); + } + else + { + tmpDivq.ref() -= fvm::laplacian(this->alpha()*this->alphaEff(), he); + + volScalarField heNew + ( + volScalarField::New + ( + "he", + he.mesh(), + dimensionedScalar(he.dimensions(), 0) + ) + ); + + surfaceScalarField sumJ + ( + surfaceScalarField::New + ( + "sumJ", + he.mesh(), + dimensionedScalar(dimMass/dimArea/dimTime, 0) + ) + ); + + surfaceScalarField sumJh + ( + surfaceScalarField::New + ( + "sumJh", + he.mesh(), + dimensionedScalar(sumJ.dimensions()*he.dimensions(), 0) + ) + ); + + forAll(Y, i) + { + if (i != d) + { + const volScalarField hi + ( + composition.HE(i, this->thermo().p(), this->thermo().T()) + ); + + heNew += Y[i]*hi; + + const surfaceScalarField ji(this->j(Y[i])); + sumJ += ji; + + sumJh += ji*fvc::interpolate(hi); + } + } + + { + const label i = d; + + const volScalarField hi + ( + composition.HE(i, this->thermo().p(), this->thermo().T()) + ); + + heNew += Y[i]*hi; + + sumJh -= sumJ*fvc::interpolate(hi); + } + + tmpDivq.ref() += + fvc::laplacian(this->alpha()*this->alphaEff(), heNew); + + tmpDivq.ref() += fvc::div(sumJh*he.mesh().magSf()); + } + + return tmpDivq; +} + + +template +tmp MaxwellStefan::j +( + const volScalarField& Yi +) const +{ + const basicSpecieMixture& composition = this->thermo().composition(); + return + BasicThermophysicalTransportModel::j(Yi) + + jexp_[composition.index(Yi)]; +} + + +template +tmp MaxwellStefan::divj +( + volScalarField& Yi +) const +{ + const basicSpecieMixture& composition = this->thermo().composition(); + return + BasicThermophysicalTransportModel::divj(Yi) + + fvc::div(jexp_[composition.index(Yi)]*Yi.mesh().magSf()); +} + + +template +void MaxwellStefan:: +transformDiffusionCoefficient() +{ + const basicSpecieMixture& composition = this->thermo().composition(); + const label d = composition.defaultSpecie(); + + // Calculate the molecular weight of the mixture and the mole fractions + scalar Wm = 0; + + forAll(W, i) + { + X[i] = Y[i]/W[i]; + Wm += X[i]; + } + + Wm = 1/Wm; + X *= Wm; + + // i counter for the specie sub-system without the default specie + label is = 0; + + // Calculate the A and B matrices from the binary diffusion coefficients + // and specie mole fractions + forAll(X, i) + { + if (i != d) + { + A(is, is) = -X[i]*Wm/(DD(i, d)*W[d]); + B(is, is) = -(X[i]*Wm/W[d] + (1 - X[i])*Wm/W[i]); + + // j counter for the specie sub-system without the default specie + label js = 0; + + forAll(X, j) + { + if (j != i) + { + A(is, is) -= X[j]*Wm/(DD(i, j)*W[i]); + + if (j != d) + { + A(is, js) = + X[i]*(Wm/(DD(i, j)*W[j]) - Wm/(DD(i, d)*W[d])); + + B(is, js) = X[i]*(Wm/W[j] - Wm/W[d]); + } + } + + if (j != d) + { + js++; + } + } + + is++; + } + } + + // LU decompose A and invert + A.decompose(); + A.inv(invA); + + // Calculate the generalized Fick's law diffusion coefficients + multiply(D, invA, B); +} + + +template +void MaxwellStefan:: +transformDiffusionCoefficientFields() +{ + const basicSpecieMixture& composition = this->thermo().composition(); + const label d = composition.defaultSpecie(); + + // For each cell or patch face + forAll(*(YPtrs[0]), pi) + { + forAll(W, i) + { + // Map YPtrs -> Y + Y[i] = (*YPtrs[i])[pi]; + + // Map DijPtrs -> DD + forAll(W, j) + { + DD(i, j) = (*DijPtrs[i][j])[pi]; + } + } + + // Transform DD -> D + transformDiffusionCoefficient(); + + // i counter for the specie sub-system without the default specie + label is = 0; + + forAll(W, i) + { + if (i != d) + { + // j counter for the specie sub-system + // without the default specie + label js = 0; + + // Map D -> DijPtrs + forAll(W, j) + { + if (j != d) + { + (*DijPtrs[i][j])[pi] = D(is, js); + + js++; + } + } + + is++; + } + } + } +} + + +template +void MaxwellStefan::transform +( + List>& Dij +) +{ + const basicSpecieMixture& composition = this->thermo().composition(); + const PtrList& Y = composition.Y(); + const volScalarField& Y0 = Y[0]; + + forAll(W, i) + { + // Map composition.Y() internal fields -> YPtrs + YPtrs[i] = &Y[i].primitiveField(); + + // Map Dii_ internal fields -> DijPtrs + DijPtrs[i][i] = &Dii_[i].primitiveFieldRef(); + + // Map Dij internal fields -> DijPtrs + forAll(W, j) + { + if (j != i) + { + DijPtrs[i][j] = &Dij[i][j].primitiveFieldRef(); + } + } + } + + // Transform binary diffusion coefficients internal field DijPtrs -> + // generalized Fick's law diffusion coefficients DijPtrs + transformDiffusionCoefficientFields(); + + forAll(Y0.boundaryField(), patchi) + { + forAll(W, i) + { + // Map composition.Y() patch fields -> YPtrs + YPtrs[i] = &Y[i].boundaryField()[patchi]; + + // Map Dii_ patch fields -> DijPtrs + DijPtrs[i][i] = &Dii_[i].boundaryFieldRef()[patchi]; + + // Map Dij patch fields -> DijPtrs + forAll(W, j) + { + if (j != i) + { + DijPtrs[i][j] = &Dij[i][j].boundaryFieldRef()[patchi]; + } + } + } + + // Transform binary diffusion coefficients patch field DijPtrs -> + // generalized Fick's law diffusion coefficients DijPtrs + transformDiffusionCoefficientFields(); + } +} + + +template +void MaxwellStefan::correct() +{ + BasicThermophysicalTransportModel::correct(); + + const basicSpecieMixture& composition = this->thermo().composition(); + const label d = composition.defaultSpecie(); + + const PtrList& Y = composition.Y(); + const volScalarField& p = this->thermo().T(); + const volScalarField& T = this->thermo().T(); + const volScalarField& rho = this->momentumTransport().rho(); + + List> Dij(Y.size()); + + // Evaluate the specie binary mass diffusion coefficient functions + // and initialise the explicit part of the specie mass flux fields + forAll(Y, i) + { + if (i != d) + { + if (jexp_.set(i)) + { + jexp_[i] = Zero; + } + else + { + jexp_.set + ( + i, + surfaceScalarField::New + ( + "jexp" + Y[i].name(), + Y[i].mesh(), + dimensionedScalar(dimensionSet(1, -2, -1, 0, 0), 0) + ) + ); + } + } + + Dii_.set(i, evaluate(D_[i][i], dimViscosity, p, T)); + + Dij[i].setSize(Y.size()); + + forAll(Y, j) + { + if (j > i) + { + Dij[i].set(j, evaluate(D_[i][j], dimViscosity, p, T)); + } + else if (j < i) + { + Dij[i].set(j, Dij[j][i].clone()); + } + } + } + + //- Transform the binary mass diffusion coefficients into the + // the generalized Fick's law diffusion coefficients + transform(Dij); + + // Accumulate the explicit part of the specie mass flux fields + forAll(Y, j) + { + if (j != d) + { + const surfaceScalarField snGradYj(fvc::snGrad(Y[j])); + + forAll(Y, i) + { + if (i != d && i != j) + { + jexp_[i] -= fvc::interpolate(rho*Dij[i][j])*snGradYj; + } + } + } + } + + // Optionally add the Soret thermal diffusion contribution to the + // explicit part of the specie mass flux fields + if (DT_.size()) + { + const surfaceScalarField gradTbyT(fvc::snGrad(T)/fvc::interpolate(T)); + + forAll(Y, i) + { + if (i != d) + { + jexp_[i] -= fvc::interpolate + ( + evaluate(DT_[i], dimDynamicViscosity, p, T) + )*gradTbyT; + } + } + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.H b/src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.H new file mode 100644 index 0000000000..78e7e34db2 --- /dev/null +++ b/src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.H @@ -0,0 +1,230 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2021 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 . + +Class + Foam::turbulenceThermophysicalTransportModels::MaxwellStefan + +Description + Base class for multi-component Maxwell Stefan generalized Fick's law + diffusion coefficients based temperature gradient heat flux model with + optional Soret thermal diffusion of species. + + The binary diffusion coefficients are specified as Function2s of + pressure and temperature but independent of composition. + + The heat flux source is implemented as an implicit energy correction to the + temperature gradient based flux source. At convergence the energy + correction is 0. + + References: + \verbatim + Taylor, R., & Krishna, R. (1993). + Multicomponent mass transfer (Vol. 2). + John Wiley & Sons. + + Merk, H. J. (1959). + The macroscopic equations for simultaneous heat and mass transfer + in isotropic, continuous and closed systems. + Applied Scientific Research, + Section A, 8(1), 73-99. + \endverbatim + +SourceFiles + MaxwellStefan.C + +\*---------------------------------------------------------------------------*/ + +#include "Function2.H" +#include "LUscalarMatrix.H" + +#ifndef MaxwellStefan_H +#define MaxwellStefan_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class MaxwellStefan Declaration +\*---------------------------------------------------------------------------*/ + +template +class MaxwellStefan +: + public BasicThermophysicalTransportModel +{ + // Private data + + // Model coefficients + + //- Array of specie binary mass diffusion coefficient functions + // [m^2/s] + List>> D_; + + //- List of specie Soret thermal diffusion coefficient + // functions [kg/m/s] + PtrList> DT_; + + //- Generalized Fick's law diffusion coefficients field list. + // This is the diagonal of the mass diffusion coefficient matrix + PtrList Dii_; + + //- List of fields of the explicit part of the mass diffusion flux + // of the species + PtrList jexp_; + + + // Workspace for diffusion coefficient transformation + + //- Molecular weights of the species + scalarField W; + + //- List of mass-fraction field pointers + // for the internal mesh or a patch + List YPtrs; + + //- Matrix of mass diffusion coefficient field pointers + // for the internal mesh or a patch + SquareMatrix DijPtrs; + + //- Mass-fractions at a cell or face + scalarField Y; + + //- Mole-fractions at a cell or face + scalarField X; + + //- Binary mass diffusion coefficient matrix at a cell or face + scalarSquareMatrix DD; + + //- Matrix form of the coefficients in the Maxwell-Stefan equation + // at a cell or face + LUscalarMatrix A; + + //- Matrix of composition coefficients at a cell or face + // used to transform the binary mass diffusion coefficients into + // the generalized Fick's law diffusion coefficients + scalarSquareMatrix B; + + //- Inverse of A + scalarSquareMatrix invA; + + //- Matrix of the generalized Fick's law diffusion coefficients + // at a cell or face + scalarSquareMatrix D; + + + // Private member functions + + //- Transform the Binary mass diffusion coefficient matrix DD into the + // matrix of the generalized Fick's law diffusion coefficients D + void transformDiffusionCoefficient(); + + //- Calls transformDiffusionCoefficient() for each cell and patch face + // and maps the resulting D into DijPtrs + void transformDiffusionCoefficientFields(); + + //- Calls transformDiffusionCoefficientFields() for the internal mesh + // and each patch and maps the resulting DijPtrs into Dii_ and Dij + void transform(List>& Dij); + + +public: + + typedef typename BasicThermophysicalTransportModel::alphaField + alphaField; + + typedef typename + BasicThermophysicalTransportModel::momentumTransportModel + momentumTransportModel; + + typedef typename BasicThermophysicalTransportModel::thermoModel + thermoModel; + + + // Constructors + + //- Construct from a momentum transport model and a thermo model + MaxwellStefan + ( + const word& type, + const momentumTransportModel& momentumTransport, + const thermoModel& thermo + ); + + + //- Destructor + virtual ~MaxwellStefan() + {} + + + // Member Functions + + //- Read thermophysicalTransport dictionary + virtual bool read(); + + //- Effective mass diffusion coefficient + // for a given specie mass-fraction [kg/m/s] + // This is the self-diffusion coefficient + virtual tmp DEff(const volScalarField& Yi) const; + + //- Effective mass diffusion coefficient + // for a given specie mass-fraction for patch [kg/m/s] + virtual tmp DEff + ( + const volScalarField& Yi, + const label patchi + ) const; + + //- Return the heat flux [W/m^2] + virtual tmp q() const; + + //- Return the source term for the energy equation + virtual tmp divq(volScalarField& he) const; + + //- Return the specie flux for the given specie mass-fraction [kg/m^2/s] + virtual tmp j(const volScalarField& Yi) const; + + //- Return the source term for the given specie mass-fraction equation + virtual tmp divj(volScalarField& Yi) const; + + //- Update the multi-component diffusivities and flux corrections + virtual void correct(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "MaxwellStefan.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.C b/src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.C new file mode 100644 index 0000000000..168fb1dfcf --- /dev/null +++ b/src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.C @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2020-2021 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 . + +\*---------------------------------------------------------------------------*/ + +#include "MaxwellStefanFourier.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace laminarThermophysicalTransportModels +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +MaxwellStefanFourier:: +MaxwellStefanFourier +( + const momentumTransportModel& momentumTransport, + const thermoModel& thermo +) +: + MaxwellStefan> + ( + typeName, + momentumTransport, + thermo + ) +{ + read(); + this->correct(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool +MaxwellStefanFourier::read() +{ + return MaxwellStefan + < + unityLewisFourier + >::read(); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace laminarThermophysicalTransportModels +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.H b/src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.H new file mode 100644 index 0000000000..bc81273f92 --- /dev/null +++ b/src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.H @@ -0,0 +1,163 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2020-2021 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 . + +Class + Foam::laminarThermophysicalTransportModels::MaxwellStefanFourier + +Description + Multi-component Maxwell Stefan generalized Fick's law diffusion coefficients + and Fourier based temperature gradient heat flux model with optional Soret + thermal diffusion of species for laminar flow. + + The binary diffusion coefficients are specified as Function2s of + pressure and temperature but independent of composition. + + The heat flux source is implemented as an implicit energy correction to the + temperature gradient based flux source. At convergence the energy + correction is 0. + +Usage + \verbatim + laminar + { + model MaxwellStefanFourier; + + D // [m^2/s] + { + O2_O2 1e-2; + O3_O3 5e-2; + N2_N2 1e-2; + O3_O2 5e-2; + O3_N2 5e-2; + O2_N2 1e-2; + + O2 1e-2; + O3 5e-2; + N2 1e-2; + } + + DT // [kg/m/s] Optional + { + O2 1e-2; + O3 5e-2; + N2 1e-2; + } + } + \endverbatim + +SourceFiles + MaxwellStefanFourier.C + +\*---------------------------------------------------------------------------*/ + +#include "MaxwellStefan.H" +#include "unityLewisFourier.H" + +#ifndef MaxwellStefanFourier_H +#define MaxwellStefanFourier_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace laminarThermophysicalTransportModels +{ + +/*---------------------------------------------------------------------------*\ + Class MaxwellStefanFourier Declaration +\*---------------------------------------------------------------------------*/ + +template +class MaxwellStefanFourier +: + public MaxwellStefan + < + unityLewisFourier + > +{ + +protected: + + // Protected data + + // Model coefficients + + //- Turbulent Schmidt number [] + dimensionedScalar Sct_; + + +public: + + typedef typename laminarThermophysicalTransportModel::alphaField + alphaField; + + typedef typename + laminarThermophysicalTransportModel::momentumTransportModel + momentumTransportModel; + + typedef typename laminarThermophysicalTransportModel::thermoModel + thermoModel; + + + //- Runtime type information + TypeName("MaxwellStefanFourier"); + + + // Constructors + + //- Construct from a momentum transport model and a thermo model + MaxwellStefanFourier + ( + const momentumTransportModel& momentumTransport, + const thermoModel& thermo + ); + + + //- Destructor + virtual ~MaxwellStefanFourier() + {} + + + // Member Functions + + //- Read thermophysicalTransport dictionary + virtual bool read(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace laminarThermophysicalTransportModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "MaxwellStefanFourier.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //