From 8182a0cff62cb363f3f10cd6e9ca02849a6d0c35 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Tue, 12 Jan 2021 17:57:17 +0000 Subject: [PATCH] ThermophysicalTransportModels::MaxwellStefanFourier: New multi-component thermophysical transport model for laminar flow 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. 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 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 --- .../matrices/LUscalarMatrix/LUscalarMatrix.C | 17 +- .../matrices/LUscalarMatrix/LUscalarMatrix.H | 8 +- ...uidReactionThermophysicalTransportModels.C | 3 + .../laminar/Fickian/Fickian.C | 31 +- .../laminar/Fickian/Fickian.H | 4 +- .../laminar/FickianFourier/FickianFourier.H | 21 +- .../laminar/MaxwellStefan/MaxwellStefan.C | 697 ++++++++++++++++++ .../laminar/MaxwellStefan/MaxwellStefan.H | 230 ++++++ .../MaxwellStefanFourier.C | 75 ++ .../MaxwellStefanFourier.H | 163 ++++ 10 files changed, 1224 insertions(+), 25 deletions(-) create mode 100644 src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.C create mode 100644 src/ThermophysicalTransportModels/laminar/MaxwellStefan/MaxwellStefan.H create mode 100644 src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.C create mode 100644 src/ThermophysicalTransportModels/laminar/MaxwellStefanFourier/MaxwellStefanFourier.H 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 + +// ************************************************************************* //