From 603b3e65b25fd5b894d7ee7e0dda65b9927b456f Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 8 May 2014 11:44:40 +0100 Subject: [PATCH] fvOptions: Correct handling of density and add multiphase support --- src/fvOptions/Make/options | 1 + src/fvOptions/fvOptions/fvOption.C | 139 +++++++++++- src/fvOptions/fvOptions/fvOption.H | 91 ++++++++ src/fvOptions/fvOptions/fvOptionList.H | 31 ++- .../fvOptions/fvOptionListTemplates.C | 87 ++++++-- .../sources/derived/MRFSource/MRFSource.C | 31 ++- .../sources/derived/MRFSource/MRFSource.H | 16 +- .../actuationDiskSource/actuationDiskSource.C | 63 +++--- .../actuationDiskSource/actuationDiskSource.H | 20 +- .../effectivenessHeatExchangerSource.C | 3 +- .../effectivenessHeatExchangerSource.H | 30 ++- .../explicitPorositySource.C | 31 +-- .../explicitPorositySource.H | 12 +- .../pressureGradientExplicitSource.C | 13 +- .../pressureGradientExplicitSource.H | 18 +- .../radialActuationDiskSource.C | 63 +++--- .../radialActuationDiskSource.H | 18 +- .../derived/rotorDiskSource/rotorDiskSource.C | 87 +++++--- .../derived/rotorDiskSource/rotorDiskSource.H | 35 +-- .../rotorDiskSource/rotorDiskSourceI.H | 21 +- .../trimModel/fixed/fixedTrim.C | 18 +- .../trimModel/fixed/fixedTrim.H | 16 +- .../trimModel/targetCoeff/targetCoeffTrim.C | 207 ++++++++++-------- .../trimModel/targetCoeff/targetCoeffTrim.H | 27 ++- .../trimModel/trimModel/trimModel.H | 16 +- .../sources/general/codedSource/CodedSource.C | 19 ++ .../sources/general/codedSource/CodedSource.H | 11 +- .../semiImplicitSource/SemiImplicitSource.C | 20 +- .../semiImplicitSource/SemiImplicitSource.H | 19 +- .../interRegionExplicitPorositySource.C | 146 +++++++----- .../interRegionExplicitPorositySource.H | 16 +- .../interRegionHeatTransferModel.C | 25 ++- .../interRegionHeatTransferModel.H | 19 +- 33 files changed, 980 insertions(+), 389 deletions(-) diff --git a/src/fvOptions/Make/options b/src/fvOptions/Make/options index a98e673d25..d4edb59316 100644 --- a/src/fvOptions/Make/options +++ b/src/fvOptions/Make/options @@ -3,6 +3,7 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \ diff --git a/src/fvOptions/fvOptions/fvOption.C b/src/fvOptions/fvOptions/fvOption.C index 884e090567..011b8e6e9f 100644 --- a/src/fvOptions/fvOptions/fvOption.C +++ b/src/fvOptions/fvOptions/fvOption.C @@ -403,13 +403,21 @@ void Foam::fv::option::correct(volTensorField& fld) } -void Foam::fv::option::addSup(fvMatrix& eqn, const label fieldI) +void Foam::fv::option::addSup +( + fvMatrix& eqn, + const label fieldI +) { // do nothing } -void Foam::fv::option::addSup(fvMatrix& eqn, const label fieldI) +void Foam::fv::option::addSup +( + fvMatrix& eqn, + const label fieldI +) { // do nothing } @@ -425,18 +433,141 @@ void Foam::fv::option::addSup } -void Foam::fv::option::addSup(fvMatrix& eqn, const label fieldI) +void Foam::fv::option::addSup +( + fvMatrix& eqn, + const label fieldI +) { // do nothing } -void Foam::fv::option::addSup(fvMatrix& eqn, const label fieldI) +void Foam::fv::option::addSup +( + fvMatrix& eqn, + const label fieldI +) { // do nothing } +void Foam::fv::option::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + // do nothing +} + + +void Foam::fv::option::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + // do nothing +} + + +void Foam::fv::option::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + // do nothing +} + + +void Foam::fv::option::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + // do nothing +} + + +void Foam::fv::option::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + // do nothing +} + + +void Foam::fv::option::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + addSup(alpha*rho, eqn, fieldI); +} + + +void Foam::fv::option::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + addSup(alpha*rho, eqn, fieldI); +} + + +void Foam::fv::option::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + addSup(alpha*rho, eqn, fieldI); +} + + +void Foam::fv::option::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + addSup(alpha*rho, eqn, fieldI); +} + + +void Foam::fv::option::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + addSup(alpha*rho, eqn, fieldI); +} + + void Foam::fv::option::setValue(fvMatrix& eqn, const label fieldI) { // do nothing diff --git a/src/fvOptions/fvOptions/fvOption.H b/src/fvOptions/fvOptions/fvOption.H index e49462ed4a..771ef7e136 100644 --- a/src/fvOptions/fvOptions/fvOption.H +++ b/src/fvOptions/fvOptions/fvOption.H @@ -384,6 +384,97 @@ public: ); + // Add explicit and implicit contributions to compressible equations + + //- Scalar + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Vector + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Spherical tensor + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Symmetric tensor + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Tensor + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + + // Add explicit and implicit contributions to phase equations + + //- Scalar + virtual void addSup + ( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Vector + virtual void addSup + ( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Spherical tensor + virtual void addSup + ( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Symmetric tensor + virtual void addSup + ( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + //- Tensor + virtual void addSup + ( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + // Set values directly //- Scalar diff --git a/src/fvOptions/fvOptions/fvOptionList.H b/src/fvOptions/fvOptions/fvOptionList.H index d4d99df36e..cfeea7b381 100644 --- a/src/fvOptions/fvOptions/fvOptionList.H +++ b/src/fvOptions/fvOptions/fvOptionList.H @@ -132,18 +132,37 @@ public: ); //- Return source for equation - template + template tmp > operator() ( - const RhoType& rho, + const volScalarField& rho, GeometricField& fld ); //- Return source for equation with specified name - template + template tmp > operator() ( - const RhoType& rho, + const volScalarField& rho, + GeometricField& fld, + const word& fieldName + ); + + //- Return source for equation + template + tmp > operator() + ( + const volScalarField& alpha, + const volScalarField& rho, + GeometricField& fld + ); + + //- Return source for equation with specified name + template + tmp > operator() + ( + const volScalarField& alpha, + const volScalarField& rho, GeometricField& fld, const word& fieldName ); @@ -155,10 +174,6 @@ public: template void constrain(fvMatrix& eqn); - //- Apply constraints to equation with specified name - template - void constrain(fvMatrix& eqn, const word& fieldName); - // Flux manipulations diff --git a/src/fvOptions/fvOptions/fvOptionListTemplates.C b/src/fvOptions/fvOptions/fvOptionListTemplates.C index 24462ce59d..08406e4f88 100644 --- a/src/fvOptions/fvOptions/fvOptionListTemplates.C +++ b/src/fvOptions/fvOptions/fvOptionListTemplates.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,10 +80,8 @@ Foam::tmp > Foam::fv::optionList::operator() const dimensionSet ds = fld.dimensions()/dimTime*dimVolume; tmp > tmtx(new fvMatrix(fld, ds)); - fvMatrix& mtx = tmtx(); - forAll(*this, i) { option& source = this->operator[](i); @@ -111,10 +109,10 @@ Foam::tmp > Foam::fv::optionList::operator() } -template +template Foam::tmp > Foam::fv::optionList::operator() ( - const RhoType& rho, + const volScalarField& rho, GeometricField& fld ) { @@ -122,10 +120,10 @@ Foam::tmp > Foam::fv::optionList::operator() } -template +template Foam::tmp > Foam::fv::optionList::operator() ( - const RhoType& rho, + const volScalarField& rho, GeometricField& fld, const word& fieldName ) @@ -135,10 +133,8 @@ Foam::tmp > Foam::fv::optionList::operator() const dimensionSet ds = rho.dimensions()*fld.dimensions()/dimTime*dimVolume; tmp > tmtx(new fvMatrix(fld, ds)); - fvMatrix& mtx = tmtx(); - forAll(*this, i) { option& source = this->operator[](i); @@ -157,7 +153,63 @@ Foam::tmp > Foam::fv::optionList::operator() << fieldName << endl; } - source.addSup(mtx, fieldI); + source.addSup(rho, mtx, fieldI); + } + } + } + + return tmtx; +} + + +template +Foam::tmp > Foam::fv::optionList::operator() +( + const volScalarField& alpha, + const volScalarField& rho, + GeometricField& fld +) +{ + return this->operator()(alpha, rho, fld, fld.name()); +} + + +template +Foam::tmp > Foam::fv::optionList::operator() +( + const volScalarField& alpha, + const volScalarField& rho, + GeometricField& fld, + const word& fieldName +) +{ + checkApplied(); + + const dimensionSet ds = + alpha.dimensions()*rho.dimensions()*fld.dimensions()/dimTime*dimVolume; + + tmp > tmtx(new fvMatrix(fld, ds)); + fvMatrix& mtx = tmtx(); + + forAll(*this, i) + { + option& source = this->operator[](i); + + label fieldI = source.applyToField(fieldName); + + if (fieldI != -1) + { + source.setApplied(fieldI); + + if (source.isActive()) + { + if (debug) + { + Info<< "Applying source " << source.name() << " to field " + << fieldName << endl; + } + + source.addSup(alpha, rho, mtx, fieldI); } } } @@ -168,17 +220,6 @@ Foam::tmp > Foam::fv::optionList::operator() template void Foam::fv::optionList::constrain(fvMatrix& eqn) -{ - constrain(eqn, eqn.psi().name()); -} - - -template -void Foam::fv::optionList::constrain -( - fvMatrix& eqn, - const word& fieldName -) { checkApplied(); @@ -186,7 +227,7 @@ void Foam::fv::optionList::constrain { option& source = this->operator[](i); - label fieldI = source.applyToField(fieldName); + label fieldI = source.applyToField(eqn.psi().name()); if (fieldI != -1) { @@ -197,7 +238,7 @@ void Foam::fv::optionList::constrain if (debug) { Info<< "Applying constraint " << source.name() - << " to field " << fieldName << endl; + << " to field " << eqn.psi().name() << endl; } source.setValue(eqn, fieldI); diff --git a/src/fvOptions/sources/derived/MRFSource/MRFSource.C b/src/fvOptions/sources/derived/MRFSource/MRFSource.C index 5dc8ba9fbe..59fd4ff021 100644 --- a/src/fvOptions/sources/derived/MRFSource/MRFSource.C +++ b/src/fvOptions/sources/derived/MRFSource/MRFSource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -89,8 +89,7 @@ Foam::fv::MRFSource::MRFSource : option(name, modelType, dict, mesh), mrfPtr_(NULL), - UName_(coeffs_.lookupOrDefault("UName", "U")), - rhoName_(coeffs_.lookupOrDefault("rhoName", "rho")) + UName_(coeffs_.lookupOrDefault("UName", "U")) { initialise(); } @@ -104,19 +103,20 @@ void Foam::fv::MRFSource::addSup const label fieldI ) { - if (eqn.dimensions() == dimForce) - { - const volScalarField& rho = - mesh_.lookupObject(rhoName_); + // Add to rhs of equation + mrfPtr_->addCoriolis(eqn, true); +} - // use 'true' flag to add to rhs of equation - mrfPtr_->addCoriolis(rho, eqn, true); - } - else - { - // use 'true' flag to add to rhs of equation - mrfPtr_->addCoriolis(eqn, true); - } + +void Foam::fv::MRFSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + // Add to rhs of equation + mrfPtr_->addCoriolis(rho, eqn, true); } @@ -173,7 +173,6 @@ bool Foam::fv::MRFSource::read(const dictionary& dict) if (option::read(dict)) { coeffs_.readIfPresent("UName", UName_); - coeffs_.readIfPresent("rhoName", rhoName_); initialise(); diff --git a/src/fvOptions/sources/derived/MRFSource/MRFSource.H b/src/fvOptions/sources/derived/MRFSource/MRFSource.H index 8f99788bc9..7b5edb1895 100644 --- a/src/fvOptions/sources/derived/MRFSource/MRFSource.H +++ b/src/fvOptions/sources/derived/MRFSource/MRFSource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -79,9 +79,6 @@ protected: //- Velocity field name, default = U word UName_; - //- Density field name, default = rho - word rhoName_; - // Protected Member Functions @@ -135,6 +132,17 @@ public: ); + // Add explicit and implicit contributions to compressible equation + + //- Vector + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + // Flux manipulations //- Make the given absolute flux relative diff --git a/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C b/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C index be83564551..be1ec03c26 100644 --- a/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C +++ b/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -115,40 +115,45 @@ void Foam::fv::actuationDiskSource::addSup const label fieldI ) { - bool compressible = false; - if (eqn.dimensions() == dimForce) - { - compressible = true; - } - const scalarField& cellsV = mesh_.V(); vectorField& Usource = eqn.source(); const vectorField& U = eqn.psi(); if (V() > VSMALL) { - if (compressible) - { - addActuationDiskAxialInertialResistance - ( - Usource, - cells_, - cellsV, - mesh_.lookupObject("rho"), - U - ); - } - else - { - addActuationDiskAxialInertialResistance - ( - Usource, - cells_, - cellsV, - geometricOneField(), - U - ); - } + addActuationDiskAxialInertialResistance + ( + Usource, + cells_, + cellsV, + geometricOneField(), + U + ); + } +} + + +void Foam::fv::actuationDiskSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + const scalarField& cellsV = mesh_.V(); + vectorField& Usource = eqn.source(); + const vectorField& U = eqn.psi(); + + if (V() > VSMALL) + { + addActuationDiskAxialInertialResistance + ( + Usource, + cells_, + cellsV, + rho, + U + ); } } diff --git a/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H b/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H index 3d95d8f70f..23ce3d9b94 100644 --- a/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H +++ b/src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -187,10 +187,22 @@ public: } - // Public Functions + // Add explicit and implicit contributions - //- Source term to fvMatrix - virtual void addSup(fvMatrix& eqn, const label fieldI); + //- Source term to momentum equation + virtual void addSup + ( + fvMatrix& eqn, + const label fieldI + ); + + //- Source term to compressible momentum equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); // I-O diff --git a/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C b/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C index 8fc1053095..68e51f40f1 100644 --- a/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C +++ b/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -204,6 +204,7 @@ bool Foam::fv::effectivenessHeatExchangerSource::alwaysApply() const void Foam::fv::effectivenessHeatExchangerSource::addSup ( + const volScalarField& rho, fvMatrix& eqn, const label ) diff --git a/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H b/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H index b60e7f8cb7..9977d0fa1b 100644 --- a/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H +++ b/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -247,10 +247,32 @@ public: virtual bool alwaysApply() const; - // Public Functions + // Add explicit and implicit contributions - //- Source term to fvMatrix - virtual void addSup(fvMatrix& eqn, const label fieldI); + //- Scalar + virtual void addSup + ( + fvMatrix& eqn, + const label fieldI + ) + { + notImplemented + ( + "effectivenessHeatExchangerSource::addSup(eqn, fieldI): " + "only compressible solvers supported." + ); + } + + + // Add explicit and implicit contributions to compressible equation + + //- Scalar + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); // I-O diff --git a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C index c398085a0a..00cf50ba09 100644 --- a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C +++ b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -86,7 +86,6 @@ Foam::fv::explicitPorositySource::explicitPorositySource option(name, modelType, dict, mesh), porosityPtr_(NULL), UName_(coeffs_.lookupOrDefault("UName", "U")), - rhoName_(coeffs_.lookupOrDefault("rhoName", "rho")), muName_(coeffs_.lookupOrDefault("muName", "thermo:mu")) { initialise(); @@ -103,18 +102,23 @@ void Foam::fv::explicitPorositySource::addSup { fvMatrix porosityEqn(eqn.psi(), eqn.dimensions()); - if (eqn.dimensions() == dimForce) - { - const volScalarField& rho = - mesh_.lookupObject(rhoName_); - const volScalarField& mu = mesh_.lookupObject(muName_); + porosityPtr_->addResistance(porosityEqn); - porosityPtr_->addResistance(porosityEqn, rho, mu); - } - else - { - porosityPtr_->addResistance(porosityEqn); - } + eqn -= porosityEqn; +} + + +void Foam::fv::explicitPorositySource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + fvMatrix porosityEqn(eqn.psi(), eqn.dimensions()); + + const volScalarField& mu = mesh_.lookupObject(muName_); + porosityPtr_->addResistance(porosityEqn, rho, mu); eqn -= porosityEqn; } @@ -132,7 +136,6 @@ bool Foam::fv::explicitPorositySource::read(const dictionary& dict) if (option::read(dict)) { coeffs_.readIfPresent("UName", UName_); - coeffs_.readIfPresent("rhoName", rhoName_); coeffs_.readIfPresent("muName", muName_); return true; diff --git a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H index a89503ecf7..2692fd9759 100644 --- a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H +++ b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -143,13 +143,21 @@ public: // Add explicit and implicit contributions - //- Vector + //- Add implicit contribution to momentum equation virtual void addSup ( fvMatrix& eqn, const label fieldI ); + //- Add implicit contribution to compressible momentum equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + // I-O diff --git a/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C b/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C index 5fc86e9204..b760173749 100644 --- a/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C +++ b/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -203,6 +203,17 @@ void Foam::fv::pressureGradientExplicitSource::addSup } +void Foam::fv::pressureGradientExplicitSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + this->addSup(eqn, fieldI); +} + + void Foam::fv::pressureGradientExplicitSource::setValue ( fvMatrix& eqn, diff --git a/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.H b/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.H index 0b8390f9f9..c3d7904ef5 100644 --- a/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.H +++ b/src/fvOptions/sources/derived/pressureGradientExplicitSource/pressureGradientExplicitSource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -130,8 +130,20 @@ public: //- Correct the pressure gradient virtual void correct(volVectorField& U); - //- Add explicit contribution to equation - virtual void addSup(fvMatrix& eqn, const label fieldI); + //- Add explicit contribution to momentum equation + virtual void addSup + ( + fvMatrix& eqn, + const label fieldI + ); + + //- Add explicit contribution to compressible momentum equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); //- Set 1/A coefficient virtual void setValue diff --git a/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.C b/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.C index 791e55e5f1..4db0484d31 100644 --- a/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.C +++ b/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -68,40 +68,45 @@ void Foam::fv::radialActuationDiskSource::addSup const label fieldI ) { - bool compressible = false; - if (eqn.dimensions() == dimForce) - { - compressible = true; - } - const scalarField& cellsV = mesh_.V(); vectorField& Usource = eqn.source(); const vectorField& U = eqn.psi(); if (V_ > VSMALL) { - if (compressible) - { - addRadialActuationDiskAxialInertialResistance - ( - Usource, - cells_, - cellsV, - mesh_.lookupObject("rho"), - U - ); - } - else - { - addRadialActuationDiskAxialInertialResistance - ( - Usource, - cells_, - cellsV, - geometricOneField(), - U - ); - } + addRadialActuationDiskAxialInertialResistance + ( + Usource, + cells_, + cellsV, + geometricOneField(), + U + ); + } +} + + +void Foam::fv::radialActuationDiskSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + const scalarField& cellsV = mesh_.V(); + vectorField& Usource = eqn.source(); + const vectorField& U = eqn.psi(); + + if (V_ > VSMALL) + { + addRadialActuationDiskAxialInertialResistance + ( + Usource, + cells_, + cellsV, + rho, + U + ); } } diff --git a/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.H b/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.H index 3fb20e45d5..33cd215ce9 100644 --- a/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.H +++ b/src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -145,8 +145,20 @@ public: // Member Functions - //- Source term to fvMatrix - virtual void addSup(fvMatrix& eqn, const label fieldI); + //- Source term to momentum equation + virtual void addSup + ( + fvMatrix& eqn, + const label fieldI + ); + + //- Source term to compressible momentum equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); // I-O diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C index f37570ec7e..59c1517d69 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C +++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,10 +25,9 @@ License #include "rotorDiskSource.H" #include "addToRunTimeSelectionTable.H" -#include "mathematicalConstants.H" #include "trimModel.H" -#include "unitConversion.H" #include "fvMatrices.H" +#include "geometricOneField.H" #include "syncTools.H" using namespace Foam::constant; @@ -435,7 +434,6 @@ Foam::fv::rotorDiskSource::rotorDiskSource ) : option(name, modelType, dict, mesh), - rhoName_("none"), rhoRef_(1.0), omega_(0.0), nBlades_(0), @@ -465,8 +463,10 @@ Foam::fv::rotorDiskSource::~rotorDiskSource() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template void Foam::fv::rotorDiskSource::calculate ( + const RhoFieldType& rho, const vectorField& U, const scalarField& thetag, vectorField& force, @@ -475,8 +475,6 @@ void Foam::fv::rotorDiskSource::calculate ) const { const scalarField& V = mesh_.V(); - const bool compressible = this->compressible(); - tmp trho(rho()); // logging info scalar dragEff = 0.0; @@ -545,11 +543,7 @@ void Foam::fv::rotorDiskSource::calculate scalar tipFactor = neg(radius/rMax_ - tipEffect_); // calculate forces perpendicular to blade - scalar pDyn = 0.5*magSqr(Uc); - if (compressible) - { - pDyn *= trho()[cellI]; - } + scalar pDyn = 0.5*rho[cellI]*magSqr(Uc); scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi; vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl); @@ -571,7 +565,6 @@ void Foam::fv::rotorDiskSource::calculate } } - if (output) { reduce(AOAmin, minOp()); @@ -594,42 +587,69 @@ void Foam::fv::rotorDiskSource::addSup const label fieldI ) { - dimensionSet dims = dimless; - if (eqn.dimensions() == dimForce) - { - coeffs_.lookup("rhoName") >> rhoName_; - dims.reset(dimForce/dimVolume); - } - else - { - coeffs_.lookup("rhoRef") >> rhoRef_; - dims.reset(dimForce/dimVolume/dimDensity); - } - volVectorField force ( IOobject ( name_ + ":rotorForce", mesh_.time().timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE + mesh_ ), mesh_, - dimensionedVector("zero", dims, vector::zero) + dimensionedVector + ( + "zero", + eqn.dimensions()/dimVolume, + vector::zero + ) ); - const volVectorField& U = eqn.psi(); - - const vectorField Uin(inflowVelocity(U)); + // Read the reference density for incompressible flow + coeffs_.lookup("rhoRef") >> rhoRef_; + const vectorField Uin(inflowVelocity(eqn.psi())); trim_->correct(Uin, force); + calculate(geometricOneField(), Uin, trim_->thetag(), force); - calculate(Uin, trim_->thetag(), force); + // Add source to rhs of eqn + eqn -= force; + + if (mesh_.time().outputTime()) + { + force.write(); + } +} - // add source to rhs of eqn +void Foam::fv::rotorDiskSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + volVectorField force + ( + IOobject + ( + name_ + ":rotorForce", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedVector + ( + "zero", + eqn.dimensions()/dimVolume, + vector::zero + ) + ); + + const vectorField Uin(inflowVelocity(eqn.psi())); + trim_->correct(rho, Uin, force); + calculate(rho, Uin, trim_->thetag(), force); + + // Add source to rhs of eqn eqn -= force; if (mesh_.time().outputTime()) @@ -653,7 +673,6 @@ bool Foam::fv::rotorDiskSource::read(const dictionary& dict) coeffs_.lookup("fieldNames") >> fieldNames_; applied_.setSize(fieldNames_.size(), false); - // read co-ordinate system/geometry invariant properties scalar rpm(readScalar(coeffs_.lookup("rpm"))); omega_ = rpm/60.0*mathematical::twoPi; diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.H b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.H index 1ac30fcfb8..0b412ba3ec 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.H +++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -37,7 +37,6 @@ Description rotorDiskSourceCoeffs { fieldNames (U); // names of fields on which to apply source - rhoName rho; // density field if compressible case nBlades 3; // number of blades tipEffect 0.96; // normalised radius above which lift = 0 @@ -95,7 +94,6 @@ SourceFiles #include "bladeModel.H" #include "profileModelList.H" #include "volFieldsFwd.H" -#include "dimensionSet.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -149,10 +147,7 @@ protected: // Protected data - //- Name of density field - word rhoName_; - - //- Reference density for rhoName = 'none' + //- Reference density for incompressible case scalar rhoRef_; //- Rotational speed [rad/s] @@ -258,7 +253,7 @@ public: // Access - //- Return the reference density for rhoName = 'none' + //- Return the reference density for incompressible case inline scalar rhoRef() const; //- Return the rotational speed [rad/s] @@ -272,18 +267,14 @@ public: //- Return the rotor co-ordinate system (r, theta, z) inline const cylindricalCS& coordSys() const; - //- Return true if solving a compressible case - inline bool compressible() const; - - //- Return the density field [kg/m3] - inline tmp rho() const; - // Evaluation //- Calculate forces + template void calculate ( + const RhoFieldType& rho, const vectorField& U, const scalarField& thetag, vectorField& force, @@ -294,8 +285,20 @@ public: // Source term addition - //- Source term to fvMatrix - virtual void addSup(fvMatrix& eqn, const label fieldI); + //- Source term to momentum equation + virtual void addSup + ( + fvMatrix& eqn, + const label fieldI + ); + + //- Source term to compressible momentum equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); // I-O diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceI.H b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceI.H index 64f77e0ce3..e923790514 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceI.H +++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -51,23 +51,4 @@ const Foam::cylindricalCS& Foam::fv::rotorDiskSource::coordSys() const } -bool Foam::fv::rotorDiskSource::compressible() const -{ - return rhoName_ != "none"; -} - - -Foam::tmp Foam::fv::rotorDiskSource::rho() const -{ - if (compressible()) - { - return mesh_.lookupObject(rhoName_); - } - else - { - return volScalarField::null(); - } -} - - // ************************************************************************* // diff --git a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C index 4379c3399d..5522250772 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C +++ b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -86,7 +86,21 @@ Foam::tmp Foam::fixedTrim::thetag() const } -void Foam::fixedTrim::correct(const vectorField& U, vectorField& force) +void Foam::fixedTrim::correct +( + const vectorField& U, + vectorField& force +) +{ + // do nothing +} + + +void Foam::fixedTrim::correct +( + const volScalarField rho, + const vectorField& U, + vectorField& force) { // do nothing } diff --git a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H index 5b8a6c2cd7..cd41018074 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H +++ b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,7 +80,19 @@ public: virtual tmp thetag() const; //- Correct the model - virtual void correct(const vectorField& U, vectorField& force); + virtual void correct + ( + const vectorField& U, + vectorField& force + ); + + //- Correct the model for compressible flow + virtual void correct + ( + const volScalarField rho, + const vectorField& U, + vectorField& force + ); }; diff --git a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C index 4de719b1f6..7ac47f1b2b 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C +++ b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,9 +24,8 @@ License \*---------------------------------------------------------------------------*/ #include "targetCoeffTrim.H" +#include "geometricOneField.H" #include "addToRunTimeSelectionTable.H" -#include "unitConversion.H" -#include "mathematicalConstants.H" using namespace Foam::constant; @@ -42,17 +41,16 @@ namespace Foam // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // +template Foam::vector Foam::targetCoeffTrim::calcCoeffs ( + const RhoFieldType& rho, const vectorField& U, const scalarField& thetag, vectorField& force ) const { - rotor_.calculate(U, thetag, force, false, false); - - bool compressible = rotor_.compressible(); - tmp trho = rotor_.rho(); + rotor_.calculate(rho, U, thetag, force, false, false); const labelList& cells = rotor_.cells(); const vectorField& C = rotor_.mesh().C(); @@ -76,11 +74,7 @@ Foam::vector Foam::targetCoeffTrim::calcCoeffs if (useCoeffs_) { scalar radius = x[i].x(); - scalar coeff2 = coeff1*pow4(radius); - if (compressible) - { - coeff2 *= trho()[cellI]; - } + scalar coeff2 = rho[cellI]*coeff1*pow4(radius); // add to coefficient vector cf[0] += (fc & yawAxis)/(coeff2 + ROOTVSMALL); @@ -101,6 +95,99 @@ Foam::vector Foam::targetCoeffTrim::calcCoeffs } +template +void Foam::targetCoeffTrim::correctTrim +( + const RhoFieldType& rho, + const vectorField& U, + vectorField& force +) +{ + if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0) + { + word calcType = "forces"; + if (useCoeffs_) + { + calcType = "coefficients"; + } + + Info<< type() << ":" << nl + << " solving for target trim " << calcType << nl; + + const scalar rhoRef = rotor_.rhoRef(); + + // iterate to find new pitch angles to achieve target force + scalar err = GREAT; + label iter = 0; + tensor J(tensor::zero); + + vector old = vector::zero; + while ((err > tol_) && (iter < nIter_)) + { + // cache initial theta vector + vector theta0(theta_); + + // set initial values + old = calcCoeffs(rho, U, thetag(), force); + + // construct Jacobian by perturbing the pitch angles + // by +/-(dTheta_/2) + for (label pitchI = 0; pitchI < 3; pitchI++) + { + theta_[pitchI] -= dTheta_/2.0; + vector cf0 = calcCoeffs(rho, U, thetag(), force); + + theta_[pitchI] += dTheta_; + vector cf1 = calcCoeffs(rho, U, thetag(), force); + + vector ddTheta = (cf1 - cf0)/dTheta_; + + J[pitchI + 0] = ddTheta[0]; + J[pitchI + 3] = ddTheta[1]; + J[pitchI + 6] = ddTheta[2]; + + theta_ = theta0; + } + + // calculate the change in pitch angle vector + vector dt = inv(J) & (target_/rhoRef - old); + + // update pitch angles + vector thetaNew = theta_ + relax_*dt; + + // update error + err = mag(thetaNew - theta_); + + // update for next iteration + theta_ = thetaNew; + iter++; + } + + if (iter == nIter_) + { + Info<< " solution not converged in " << iter + << " iterations, final residual = " << err + << "(" << tol_ << ")" << endl; + } + else + { + Info<< " final residual = " << err << "(" << tol_ + << "), iterations = " << iter << endl; + } + + Info<< " current and target " << calcType << nl + << " thrust = " << old[0]*rhoRef << ", " << target_[0] << nl + << " pitch = " << old[1]*rhoRef << ", " << target_[1] << nl + << " roll = " << old[2]*rhoRef << ", " << target_[2] << nl + << " new pitch angles [deg]:" << nl + << " theta0 = " << radToDeg(theta_[0]) << nl + << " theta1c = " << radToDeg(theta_[1]) << nl + << " theta1s = " << radToDeg(theta_[2]) << nl + << endl; + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::targetCoeffTrim::targetCoeffTrim @@ -185,90 +272,24 @@ Foam::tmp Foam::targetCoeffTrim::thetag() const } -void Foam::targetCoeffTrim::correct(const vectorField& U, vectorField& force) +void Foam::targetCoeffTrim::correct +( + const vectorField& U, + vectorField& force +) { - if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0) - { - word calcType = "forces"; - if (useCoeffs_) - { - calcType = "coefficients"; - } + correctTrim(geometricOneField(), U, force); +} - Info<< type() << ":" << nl - << " solving for target trim " << calcType << nl; - const scalar rhoRef = rotor_.rhoRef(); - - // iterate to find new pitch angles to achieve target force - scalar err = GREAT; - label iter = 0; - tensor J(tensor::zero); - - vector old = vector::zero; - while ((err > tol_) && (iter < nIter_)) - { - // cache initial theta vector - vector theta0(theta_); - - // set initial values - old = calcCoeffs(U, thetag(), force); - - // construct Jacobian by perturbing the pitch angles - // by +/-(dTheta_/2) - for (label pitchI = 0; pitchI < 3; pitchI++) - { - theta_[pitchI] -= dTheta_/2.0; - vector cf0 = calcCoeffs(U, thetag(), force); - - theta_[pitchI] += dTheta_; - vector cf1 = calcCoeffs(U, thetag(), force); - - vector ddTheta = (cf1 - cf0)/dTheta_; - - J[pitchI + 0] = ddTheta[0]; - J[pitchI + 3] = ddTheta[1]; - J[pitchI + 6] = ddTheta[2]; - - theta_ = theta0; - } - - // calculate the change in pitch angle vector - vector dt = inv(J) & (target_/rhoRef - old); - - // update pitch angles - vector thetaNew = theta_ + relax_*dt; - - // update error - err = mag(thetaNew - theta_); - - // update for next iteration - theta_ = thetaNew; - iter++; - } - - if (iter == nIter_) - { - Info<< " solution not converged in " << iter - << " iterations, final residual = " << err - << "(" << tol_ << ")" << endl; - } - else - { - Info<< " final residual = " << err << "(" << tol_ - << "), iterations = " << iter << endl; - } - - Info<< " current and target " << calcType << nl - << " thrust = " << old[0]*rhoRef << ", " << target_[0] << nl - << " pitch = " << old[1]*rhoRef << ", " << target_[1] << nl - << " roll = " << old[2]*rhoRef << ", " << target_[2] << nl - << " new pitch angles [deg]:" << nl - << " theta0 = " << radToDeg(theta_[0]) << nl - << " theta1c = " << radToDeg(theta_[1]) << nl - << " theta1s = " << radToDeg(theta_[2]) << nl - << endl; - } +void Foam::targetCoeffTrim::correct +( + const volScalarField rho, + const vectorField& U, + vectorField& force +) +{ + correctTrim(rho, U, force); } diff --git a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H index 022d0c49a5..8925e03bea 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H +++ b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -125,13 +125,24 @@ protected: // Protected member functions //- Calculate the rotor force and moment coefficients vector + template vector calcCoeffs ( + const RhoFieldType& rho, const vectorField& U, const scalarField& alphag, vectorField& force ) const; + //- Correct the model + template + void correctTrim + ( + const RhoFieldType& rho, + const vectorField& U, + vectorField& force + ); + public: @@ -154,7 +165,19 @@ public: virtual tmp thetag() const; //- Correct the model - virtual void correct(const vectorField& U, vectorField& force); + virtual void correct + ( + const vectorField& U, + vectorField& force + ); + + //- Correct the model for compressible flow + virtual void correct + ( + const volScalarField rho, + const vectorField& U, + vectorField& force + ); }; diff --git a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H index ca9b6ed21e..6f078d7c98 100644 --- a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H +++ b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -120,7 +120,19 @@ public: virtual tmp thetag() const = 0; //- Correct the model - virtual void correct(const vectorField& U, vectorField& force) = 0; + virtual void correct + ( + const vectorField& U, + vectorField& force + ) = 0; + + //- Correct the model for compressible flow + virtual void correct + ( + const volScalarField rho, + const vectorField& U, + vectorField& force + ) = 0; }; diff --git a/src/fvOptions/sources/general/codedSource/CodedSource.C b/src/fvOptions/sources/general/codedSource/CodedSource.C index ca91212c6e..765326e5db 100644 --- a/src/fvOptions/sources/general/codedSource/CodedSource.C +++ b/src/fvOptions/sources/general/codedSource/CodedSource.C @@ -181,6 +181,25 @@ void Foam::fv::CodedSource::addSup } +template +void Foam::fv::CodedSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + if (debug) + { + Info<< "CodedSource<"<< pTraits::typeName + << ">::addSup for source " << name_ << endl; + } + + updateLibrary(redirectType_); + redirectFvOption().addSup(rho, eqn, fieldI); +} + + template void Foam::fv::CodedSource::setValue ( diff --git a/src/fvOptions/sources/general/codedSource/CodedSource.H b/src/fvOptions/sources/general/codedSource/CodedSource.H index 535c166063..ed22a6d2b0 100644 --- a/src/fvOptions/sources/general/codedSource/CodedSource.H +++ b/src/fvOptions/sources/general/codedSource/CodedSource.H @@ -42,7 +42,7 @@ Description setValue ( - fvMatrix& eqn, + fvMatrix& eqn, const label fieldI ) @@ -206,6 +206,15 @@ public: const label fieldI ); + //- Explicit and implicit matrix contributions + // to compressible equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + //- Set value virtual void setValue ( diff --git a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C index 6170b41046..bf5593e26b 100644 --- a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C +++ b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -193,4 +193,22 @@ void Foam::fv::SemiImplicitSource::addSup } +template +void Foam::fv::SemiImplicitSource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + if (debug) + { + Info<< "SemiImplicitSource<" << pTraits::typeName + << ">::addSup for source " << name_ << endl; + } + + return this->addSup(eqn, fieldI); +} + + // ************************************************************************* // diff --git a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H index f956b8fb3e..496151b099 100644 --- a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H +++ b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,10 +80,10 @@ namespace fv // Forward declaration of classes - template class SemiImplicitSource; + // Forward declaration of friend functions template @@ -93,6 +93,7 @@ Ostream& operator<< const SemiImplicitSource& ); + /*---------------------------------------------------------------------------*\ Class SemiImplicitSource Declaration \*---------------------------------------------------------------------------*/ @@ -184,7 +185,19 @@ public: // Evaluation //- Add explicit contribution to equation - virtual void addSup(fvMatrix& eqn, const label fieldI); + virtual void addSup + ( + fvMatrix& eqn, + const label fieldI + ); + + //- Add explicit contribution to compressible equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); // I-O diff --git a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C index 29d72dd112..b91bfab0bd 100644 --- a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C +++ b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -122,7 +122,6 @@ Foam::fv::interRegionExplicitPorositySource::interRegionExplicitPorositySource porosityPtr_(NULL), firstIter_(-1), UName_(coeffs_.lookupOrDefault("UName", "U")), - rhoName_(coeffs_.lookupOrDefault("rhoName", "rho")), muName_(coeffs_.lookupOrDefault("muName", "thermo:mu")) { if (active_) @@ -171,64 +170,108 @@ void Foam::fv::interRegionExplicitPorositySource::addSup fvMatrix nbrEqn(UNbr, eqn.dimensions()); - if (eqn.dimensions() == dimForce) - { - volScalarField rhoNbr + porosityPtr_->addResistance(nbrEqn); + + // convert source from neighbour to local region + fvMatrix porosityEqn(U, eqn.dimensions()); + scalarField& Udiag = porosityEqn.diag(); + vectorField& Usource = porosityEqn.source(); + + Udiag.setSize(eqn.diag().size(), 0.0); + Usource.setSize(eqn.source().size(), vector::zero); + + meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp(), Udiag); + meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp(), Usource); + + eqn -= porosityEqn; +} + + +void Foam::fv::interRegionExplicitPorositySource::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + initialise(); + + const fvMesh& nbrMesh = mesh_.time().lookupObject(nbrRegionName_); + + const volVectorField& U = eqn.psi(); + + volVectorField UNbr + ( + IOobject ( - IOobject - ( - "rho:UNbr", - nbrMesh.time().timeName(), - nbrMesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), + name_ + ":UNbr", + nbrMesh.time().timeName(), nbrMesh, - dimensionedScalar("zero", dimDensity, 0.0) - ); + IOobject::NO_READ, + IOobject::NO_WRITE + ), + nbrMesh, + dimensionedVector("zero", U.dimensions(), vector::zero) + ); - volScalarField muNbr + // map local velocity onto neighbour region + meshInterp().mapSrcToTgt + ( + U.internalField(), + plusEqOp(), + UNbr.internalField() + ); + + fvMatrix nbrEqn(UNbr, eqn.dimensions()); + + volScalarField rhoNbr + ( + IOobject ( - IOobject - ( - "mu:UNbr", - nbrMesh.time().timeName(), - nbrMesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), + "rho:UNbr", + nbrMesh.time().timeName(), nbrMesh, - dimensionedScalar("zero", dimViscosity, 0.0) - ); + IOobject::NO_READ, + IOobject::NO_WRITE + ), + nbrMesh, + dimensionedScalar("zero", dimDensity, 0.0) + ); - const volScalarField& rho = - mesh_.lookupObject(rhoName_); - - const volScalarField& mu = - mesh_.lookupObject(muName_); - - // map local rho onto neighbour region - meshInterp().mapSrcToTgt + volScalarField muNbr + ( + IOobject ( - rho.internalField(), - plusEqOp(), - rhoNbr.internalField() - ); + "mu:UNbr", + nbrMesh.time().timeName(), + nbrMesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + nbrMesh, + dimensionedScalar("zero", dimViscosity, 0.0) + ); - // map local mu onto neighbour region - meshInterp().mapSrcToTgt - ( - mu.internalField(), - plusEqOp(), - muNbr.internalField() - ); + const volScalarField& mu = + mesh_.lookupObject(muName_); - porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr); - } - else - { - porosityPtr_->addResistance(nbrEqn); - } + // map local rho onto neighbour region + meshInterp().mapSrcToTgt + ( + rho.internalField(), + plusEqOp(), + rhoNbr.internalField() + ); + + // map local mu onto neighbour region + meshInterp().mapSrcToTgt + ( + mu.internalField(), + plusEqOp(), + muNbr.internalField() + ); + + porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr); // convert source from neighbour to local region fvMatrix porosityEqn(U, eqn.dimensions()); @@ -257,7 +300,6 @@ bool Foam::fv::interRegionExplicitPorositySource::read(const dictionary& dict) if (option::read(dict)) { coeffs_.readIfPresent("UName", UName_); - coeffs_.readIfPresent("rhoName", rhoName_); coeffs_.readIfPresent("muName", muName_); // reset the porosity model? diff --git a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H index 84fac78b8f..c430e5725c 100644 --- a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H +++ b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -91,9 +91,6 @@ protected: //- Velocity field name, default = U word UName_; - //- Density field name (compressible case only), default = rho - word rhoName_; - //- Dynamic viscosity field name (compressible case only) // default = thermo:mu word muName_; @@ -154,6 +151,17 @@ public: ); + // Add explicit and implicit contributions to compressible equation + + //- Vector + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); + + // I-O //- Write data diff --git a/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C index 6dd8c3f073..f3c4517025 100644 --- a/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C +++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C @@ -24,7 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "interRegionHeatTransferModel.H" -#include "fluidThermo.H" +#include "basicThermo.H" #include "fvmSup.H" #include "zeroGradientFvPatchFields.H" #include "fvcVolumeIntegrate.H" @@ -216,7 +216,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup { if (he.dimensions() == dimEnergy/dimMass) { - if (mesh_.foundObject("thermophysicalProperties")) + if (mesh_.foundObject("thermophysicalProperties")) { const basicThermo& thermo = mesh_.lookupObject("thermophysicalProperties"); @@ -237,7 +237,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup } else { - FatalErrorIn + FatalErrorIn ( "void Foam::fv::interRegionHeatTransferModel::addSup" "(" @@ -245,11 +245,9 @@ void Foam::fv::interRegionHeatTransferModel::addSup " const label " ")" ) << " on mesh " << mesh_.name() - << " could not find object fluidThermo." - << " The available objects : " + << " could not find object basicThermo." + << " The available objects are: " << mesh_.names() - << " The semi implicit option can only be used for " - << "fluid-fluid inter region heat transfer models " << exit(FatalError); } } @@ -265,12 +263,23 @@ void Foam::fv::interRegionHeatTransferModel::addSup } +void Foam::fv::interRegionHeatTransferModel::addSup +( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI +) +{ + addSup(eqn, fieldI); +} + + void Foam::fv::interRegionHeatTransferModel::writeData(Ostream& os) const { os.writeKeyword("name") << this->name() << token::END_STATEMENT << nl; os.writeKeyword("nbrRegionName") << nbrRegionName_ << token::END_STATEMENT << nl; - os.writeKeyword("nbrModeleName") << nbrModelName_ + os.writeKeyword("nbrModelName") << nbrModelName_ << token::END_STATEMENT << nl; os.writeKeyword("master") << master_ << token::END_STATEMENT << nl; os.writeKeyword("semiImplicit") << semiImplicit_ << token::END_STATEMENT diff --git a/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H index 55aabd45e1..19f340b318 100644 --- a/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H +++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -157,7 +157,6 @@ public: // Member Functions - // Access //- Return the heat transfer coefficient @@ -169,8 +168,20 @@ public: //- Return access to the neighbour model inline interRegionHeatTransferModel& nbrModel(); - //-Source term to fvMatrix - virtual void addSup(fvMatrix& eqn, const label fieldI); + //- Source term to energy equation + virtual void addSup + ( + fvMatrix& eqn, + const label fieldI + ); + + //- Source term to compressible energy equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& eqn, + const label fieldI + ); //- Calculate heat transfer coefficient virtual void calculateHtc() = 0;