diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C index b0857fa78e..f46c77e3e2 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C +++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C @@ -52,66 +52,65 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer ) : porosityModel(name, modelType, mesh, dict, cellZoneName), - coordSys_(coeffs_, mesh), - D_("D", dimless/sqr(dimLength), tensor::zero), - F_("F", dimless/dimLength, tensor::zero), + D_(cellZoneIds_.size()), + F_(cellZoneIds_.size()), rhoName_(coeffs_.lookupOrDefault("rho", "rho")), muName_(coeffs_.lookupOrDefault("mu", "thermo:mu")), nuName_(coeffs_.lookupOrDefault("nu", "nu")) { - // local-to-global transformation tensor - const tensor& E = coordSys_.R(); dimensionedVector d(coeffs_.lookup("d")); - if (D_.dimensions() != d.dimensions()) - { - FatalIOErrorIn - ( - "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer" - "(" - "const word&, " - "const word&, " - "const fvMesh&, " - "const dictionary&" - ")", - coeffs_ - ) << "incorrect dimensions for d: " << d.dimensions() - << " should be " << D_.dimensions() - << exit(FatalIOError); - } + dimensionedVector f(coeffs_.lookup("f")); adjustNegativeResistance(d); - - D_.value().xx() = d.value().x(); - D_.value().yy() = d.value().y(); - D_.value().zz() = d.value().z(); - D_.value() = (E & D_ & E.T()).value(); - - dimensionedVector f(coeffs_.lookup("f")); - if (F_.dimensions() != f.dimensions()) - { - FatalIOErrorIn - ( - "Foam::porosityModels::DarcyForchheimer::DarcyForchheimer" - "(" - "const word&, " - "const word&, " - "const fvMesh&, " - "const dictionary&" - ")", - coeffs_ - ) << "incorrect dimensions for f: " << f.dimensions() - << " should be " << F_.dimensions() - << exit(FatalIOError); - } - adjustNegativeResistance(f); - // leading 0.5 is from 1/2*rho - F_.value().xx() = 0.5*f.value().x(); - F_.value().yy() = 0.5*f.value().y(); - F_.value().zz() = 0.5*f.value().z(); - F_.value() = (E & F_ & E.T()).value(); + if (coordSys_.R().uniform()) + { + forAll (cellZoneIds_, zoneI) + { + D_[zoneI].setSize(1, tensor::zero); + F_[zoneI].setSize(1, tensor::zero); + + D_[zoneI][0].xx() = d.value().x(); + D_[zoneI][0].yy() = d.value().y(); + D_[zoneI][0].zz() = d.value().z(); + + D_[zoneI][0] = coordSys_.R().transformTensor(D_[zoneI][0]); + + // leading 0.5 is from 1/2*rho + F_[zoneI][0].xx() = 0.5*f.value().x(); + F_[zoneI][0].yy() = 0.5*f.value().y(); + F_[zoneI][0].zz() = 0.5*f.value().z(); + + F_[zoneI][0] = coordSys_.R().transformTensor(F_[zoneI][0]); + } + + } + else + { + forAll(cellZoneIds_, zoneI) + { + const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]]; + + D_[zoneI].setSize(cells.size(), tensor::zero); + F_[zoneI].setSize(cells.size(), tensor::zero); + + forAll(cells, i) + { + D_[zoneI][i].xx() = d.value().x(); + D_[zoneI][i].yy() = d.value().y(); + D_[zoneI][i].zz() = d.value().z(); + + F_[zoneI][i].xx() = f.value().x(); + F_[zoneI][i].yy() = f.value().y(); + F_[zoneI][i].zz() = f.value().z(); + } + + D_[zoneI] = coordSys_.R().transformTensor(D_[zoneI], cells); + F_[zoneI] = coordSys_.R().transformTensor(F_[zoneI], cells); + } + } } @@ -132,7 +131,7 @@ void Foam::porosityModels::DarcyForchheimer::correct const scalarField& V = mesh_.V(); scalarField& Udiag = UEqn.diag(); vectorField& Usource = UEqn.source(); - + if (UEqn.dimensions() == dimForce) { const volScalarField& rho = @@ -163,7 +162,7 @@ void Foam::porosityModels::DarcyForchheimer::correct const scalarField& V = mesh_.V(); scalarField& Udiag = UEqn.diag(); vectorField& Usource = UEqn.source(); - + apply(Udiag, Usource, V, rho, mu, U); } diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H index ab3007d204..856089cee6 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H +++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.H @@ -53,7 +53,6 @@ SourceFiles #define DarcyForchheimer_H #include "porosityModel.H" -#include "coordinateSystem.H" #include "dimensionedTensor.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -75,14 +74,12 @@ private: // Private data - //- Local co-ordinate system - coordinateSystem coordSys_; //- Darcy coefficient [1/m2] - dimensionedTensor D_; + List D_; //- Forchheimer coefficient [1/m] - dimensionedTensor F_; + List F_; //- Name of density field word rhoName_; @@ -161,7 +158,7 @@ public: virtual void correct ( const fvVectorMatrix& UEqn, - volTensorField& AU + volTensorField& AU ) const; diff --git a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C index 184206a010..27b10a79ef 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C +++ b/src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C @@ -36,18 +36,19 @@ void Foam::porosityModels::DarcyForchheimer::apply const vectorField& U ) const { - const tensor& D = D_.value(); - const tensor& F = F_.value(); - forAll(cellZoneIds_, zoneI) { + const tensorField& dZones = D_[zoneI]; + const tensorField& fZones = F_[zoneI]; + const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]]; forAll(cells, i) { const label cellI = cells[i]; - - const tensor Cd = mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F; + const label j = this->fieldIndex(i); + const tensor Cd = + mu[cellI]*dZones[j] + (rho[cellI]*mag(U[cellI]))*fZones[j]; const scalar isoCd = tr(Cd); @@ -67,16 +68,20 @@ void Foam::porosityModels::DarcyForchheimer::apply const vectorField& U ) const { - const tensor& D = D_.value(); - const tensor& F = F_.value(); - forAll(cellZoneIds_, zoneI) { + const tensorField& dZones = D_[zoneI]; + const tensorField& fZones = F_[zoneI]; + const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]]; forAll(cells, i) { const label cellI = cells[i]; + const label j = this->fieldIndex(i); + const tensor D = dZones[j]; + const tensor F = fZones[j]; + AU[cellI] += mu[cellI]*D + (rho[cellI]*mag(U[cellI]))*F; } } diff --git a/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C b/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C index 59cef74368..b8267f7152 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C +++ b/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C @@ -50,19 +50,18 @@ void Foam::porosityModels::fixedCoeff::apply const scalar rho ) const { - const tensor& alpha = alpha_.value(); - const tensor& beta = beta_.value(); - forAll(cellZoneIds_, zoneI) { + const tensorField& alphaZones = alpha_[zoneI]; + const tensorField& betaZones = beta_[zoneI]; + const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]]; forAll(cells, i) { const label cellI = cells[i]; - - const tensor Cd = rho*(alpha + beta*mag(U[cellI])); - + const label j = fieldIndex(i); + const tensor Cd = rho*(alphaZones[j] + betaZones[j]*mag(U[cellI])); const scalar isoCd = tr(Cd); Udiag[cellI] += V[cellI]*isoCd; @@ -79,16 +78,21 @@ void Foam::porosityModels::fixedCoeff::apply const scalar rho ) const { - const tensor& alpha = alpha_.value(); - const tensor& beta = beta_.value(); forAll(cellZoneIds_, zoneI) { + const tensorField& alphaZones = alpha_[zoneI]; + const tensorField& betaZones = beta_[zoneI]; + const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]]; forAll(cells, i) { const label cellI = cells[i]; + const label j = fieldIndex(i); + const tensor alpha = alphaZones[j]; + const tensor beta = betaZones[j]; + AU[cellI] += rho*(alpha + beta*mag(U[cellI])); } } @@ -107,62 +111,59 @@ Foam::porosityModels::fixedCoeff::fixedCoeff ) : porosityModel(name, modelType, mesh, dict, cellZoneName), - coordSys_(coeffs_, mesh), - alpha_("alpha", dimless/dimTime, tensor::zero), - beta_("beta", dimless/dimLength, tensor::zero) + alpha_(cellZoneIds_.size()), + beta_(cellZoneIds_.size()) { - // local-to-global transformation tensor - const tensor& E = coordSys_.R(); - dimensionedVector alpha(coeffs_.lookup("alpha")); - if (alpha_.dimensions() != alpha.dimensions()) - { - FatalIOErrorIn - ( - "Foam::porosityModels::fixedCoeff::fixedCoeff" - "(" - "const word&, " - "const word&, " - "const fvMesh&, " - "const dictionary&" - ")", - coeffs_ - ) << "incorrect dimensions for alpha: " << alpha.dimensions() - << " should be " << alpha_.dimensions() - << exit(FatalIOError); - } + dimensionedVector beta(coeffs_.lookup("beta")); adjustNegativeResistance(alpha); - - alpha_.value().xx() = alpha.value().x(); - alpha_.value().yy() = alpha.value().y(); - alpha_.value().zz() = alpha.value().z(); - alpha_.value() = (E & alpha_ & E.T()).value(); - - dimensionedVector beta(coeffs_.lookup("beta")); - if (beta_.dimensions() != beta.dimensions()) - { - FatalIOErrorIn - ( - "Foam::porosityModels::fixedCoeff::fixedCoeff" - "(" - "const word&, " - "const word&, " - "const fvMesh&, " - "const dictionary&" - ")", - coeffs_ - ) << "incorrect dimensions for beta: " << beta.dimensions() - << " should be " << beta_.dimensions() - << exit(FatalIOError); - } - adjustNegativeResistance(beta); - beta_.value().xx() = beta.value().x(); - beta_.value().yy() = beta.value().y(); - beta_.value().zz() = beta.value().z(); - beta_.value() = (E & beta_ & E.T()).value(); + if (coordSys_.R().uniform()) + { + forAll (cellZoneIds_, zoneI) + { + alpha_[zoneI].setSize(1, tensor::zero); + beta_[zoneI].setSize(1, tensor::zero); + + alpha_[zoneI][0].xx() = alpha.value().x(); + alpha_[zoneI][0].yy() = alpha.value().y(); + alpha_[zoneI][0].zz() = alpha.value().z(); + alpha_[zoneI][0] = coordSys_.R().transformTensor(alpha_[zoneI][0]); + + beta_[zoneI][0].xx() = beta.value().x(); + beta_[zoneI][0].yy() = beta.value().y(); + beta_[zoneI][0].zz() = beta.value().z(); + beta_[zoneI][0] = coordSys_.R().transformTensor(beta_[zoneI][0]); + } + } + else + { + forAll(cellZoneIds_, zoneI) + { + const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]]; + + alpha_[zoneI].setSize(cells.size(), tensor::zero); + beta_[zoneI].setSize(cells.size(), tensor::zero); + + forAll(cells, i) + { + alpha_[zoneI][i].xx() = alpha.value().x(); + alpha_[zoneI][i].yy() = alpha.value().y(); + alpha_[zoneI][i].zz() = alpha.value().z(); + + beta_[zoneI][i].xx() = beta.value().x(); + beta_[zoneI][i].yy() = beta.value().y(); + beta_[zoneI][i].zz() = beta.value().z(); + } + + alpha_[zoneI] = + coordSys_.R().transformTensor(alpha_[zoneI], cells); + + beta_[zoneI] = coordSys_.R().transformTensor(beta_[zoneI], cells); + } + } } @@ -183,7 +184,7 @@ void Foam::porosityModels::fixedCoeff::correct const scalarField& V = mesh_.V(); scalarField& Udiag = UEqn.diag(); vectorField& Usource = UEqn.source(); - + scalar rho = 1.0; if (UEqn.dimensions() == dimForce) { @@ -211,7 +212,7 @@ void Foam::porosityModels::fixedCoeff::correct { coeffs_.lookup("rhoRef") >> rho; } - + apply(Udiag, Usource, V, U, rho); } diff --git a/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H b/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H index 61f4c36224..09cf5cfa79 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H +++ b/src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.H @@ -43,7 +43,6 @@ SourceFiles #define fixedCoeff_H #include "porosityModel.H" -#include "coordinateSystem.H" #include "dimensionedTensor.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -65,14 +64,11 @@ private: // Private data - //- Local co-ordinate system - coordinateSystem coordSys_; - //- Model alpha coefficient [1/s] - dimensionedTensor alpha_; + List alpha_; //- Model beta coefficient [1/m] - dimensionedTensor beta_; + List beta_; // Private Member Functions @@ -138,7 +134,7 @@ public: virtual void correct ( const fvVectorMatrix& UEqn, - volTensorField& AU + volTensorField& AU ) const; diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C index 77de5b1ffe..3b9a9c6375 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C +++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C @@ -66,6 +66,17 @@ void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist) } +Foam::label Foam::porosityModel::fieldIndex(const label i) const +{ + label index = 0; + if (!coordSys_.R().uniform()) + { + index = i; + } + return index; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::porosityModel::porosityModel @@ -83,7 +94,8 @@ Foam::porosityModel::porosityModel coeffs_(dict.subDict(modelType + "Coeffs")), active_(true), zoneName_(cellZoneName), - cellZoneIds_() + cellZoneIds_(), + coordSys_(coordinateSystem::New(mesh, coeffs_)) { if (zoneName_ == word::null) { @@ -158,7 +170,7 @@ void Foam::porosityModel::porosityModel::addResistance ( const fvVectorMatrix& UEqn, volTensorField& AU, - bool correctAUprocBC + bool correctAUprocBC ) const { if (cellZoneIds_.empty()) diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H index 1c61e29f18..019c063f87 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H +++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H @@ -40,6 +40,7 @@ SourceFiles #include "dictionary.H" #include "fvMatricesFwd.H" #include "runTimeSelectionTables.H" +#include "coordinateSystem.H" #include "dimensionedVector.H" #include "keyType.H" @@ -90,6 +91,9 @@ protected: //- Cell zone Ids labelList cellZoneIds_; + //- Local co-ordinate system + coordinateSystem coordSys_; + // Protected Member Functions @@ -111,6 +115,9 @@ protected: volTensorField& AU ) const = 0; + //- Return label index + label fieldIndex(const label index) const; + public: