diff --git a/src/waves/Make/files b/src/waves/Make/files index 4d4f0361a..7ec670464 100644 --- a/src/waves/Make/files +++ b/src/waves/Make/files @@ -9,5 +9,6 @@ waveSuperposition/waveSuperposition.C derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C +derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C LIB = $(FOAM_LIBBIN)/libwaves diff --git a/src/waves/Make/options b/src/waves/Make/options index a3ae8da83..a3458eaa3 100644 --- a/src/waves/Make/options +++ b/src/waves/Make/options @@ -1,7 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude LIB_LIBS = \ -lfiniteVolume \ - -lmeshTools + -lmeshTools \ + -ldynamicMesh diff --git a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C new file mode 100644 index 000000000..951491cdb --- /dev/null +++ b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C @@ -0,0 +1,251 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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 "wavePressureFvPatchScalarField.H" +#include "waveVelocityFvPatchVectorField.H" +#include "addToRunTimeSelectionTable.H" +#include "levelSet.H" +#include "volFields.H" +#include "fvMeshSubset.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::wavePressureFvPatchScalarField::wavePressureFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(p, iF), + UName_("U"), + rhoName_("rho") +{ + refValue() = Zero; + refGrad() = Zero; + valueFraction() = Zero; +} + + +Foam::wavePressureFvPatchScalarField::wavePressureFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + mixedFvPatchScalarField(p, iF), + UName_(dict.lookupOrDefault("U", "U")), + rhoName_(dict.lookupOrDefault("rho", "rho")) +{ + if (dict.found("value")) + { + fvPatchScalarField::operator=(scalarField("value", dict, p.size())); + } + else + { + fvPatchScalarField::operator=(patchInternalField()); + } + + refValue() = *this; + refGrad() = Zero; + valueFraction() = Zero; +} + + +Foam::wavePressureFvPatchScalarField::wavePressureFvPatchScalarField +( + const wavePressureFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + mixedFvPatchScalarField(ptf, p, iF, mapper), + UName_(ptf.UName_), + rhoName_(ptf.rhoName_) +{} + + +Foam::wavePressureFvPatchScalarField::wavePressureFvPatchScalarField +( + const wavePressureFvPatchScalarField& ptf +) +: + mixedFvPatchScalarField(ptf), + UName_(ptf.UName_), + rhoName_(ptf.rhoName_) +{} + + +Foam::wavePressureFvPatchScalarField::wavePressureFvPatchScalarField +( + const wavePressureFvPatchScalarField& ptf, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(ptf, iF), + UName_(ptf.UName_), + rhoName_(ptf.rhoName_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp Foam::wavePressureFvPatchScalarField::p() const +{ + const waveVelocityFvPatchVectorField& Up = + refCast + ( + patch().lookupPatchField(UName_) + ); + + const scalar t = db().time().timeOutputValue(); + + return + levelSetAverage + ( + patch(), + Up.waves().height(t, patch().Cf()), + Up.waves().height(t, patch().patch().localPoints()), + Up.waves().pGas(t, patch().Cf())(), + Up.waves().pGas(t, patch().patch().localPoints())(), + Up.waves().pLiquid(t, patch().Cf())(), + Up.waves().pLiquid(t, patch().patch().localPoints())() + ); +} + + +Foam::tmp Foam::wavePressureFvPatchScalarField::pn() const +{ + const waveVelocityFvPatchVectorField& Up = + refCast + ( + patch().lookupPatchField(UName_) + ); + + const scalar t = db().time().timeOutputValue(); + + const fvMeshSubset& subset = Up.faceCellSubset(); + const fvMesh& meshs = subset.subMesh(); + const label patchis = findIndex(subset.patchMap(), patch().index()); + + const scalarField ps = + levelSetAverage + ( + meshs, + Up.waves().height(t, meshs.cellCentres())(), + Up.waves().height(t, meshs.points())(), + Up.waves().pGas(t, meshs.cellCentres())(), + Up.waves().pGas(t, meshs.points())(), + Up.waves().pLiquid(t, meshs.cellCentres())(), + Up.waves().pLiquid(t, meshs.points())() + ); + + tmp tResult(new scalarField(patch().size())); + scalarField& result = tResult.ref(); + + if (patchis != -1) + { + forAll(meshs.boundary()[patchis], is) + { + const label fs = is + meshs.boundary()[patchis].patch().start(); + const label cs = meshs.boundary()[patchis].faceCells()[is]; + const label f = subset.faceMap()[fs]; + const label i = patch().patch().whichFace(f); + result[i] = ps[cs]; + } + } + + return tResult; +} + + +void Foam::wavePressureFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + + const fvPatchVectorField& Up = + patch().lookupPatchField(UName_); + + if (!isA(Up)) + { + FatalErrorInFunction + << "The corresponding condition for the velocity " + << "field " << UName_ << " on patch " << patch().name() + << " is not of type " << waveVelocityFvPatchVectorField::typeName + << exit(FatalError); + } + + const waveVelocityFvPatchVectorField& Uwp = + refCast(Up); + + if (Uwp.pName() != internalField().name()) + { + FatalErrorInFunction + << "The corresponding condition for the velocity " + << "field " << UName_ << " on patch " << patch().name() + << " does not have the pressure set to " << internalField().name() + << exit(FatalError); + } + + const scalarField p(this->p()), pn(this->pn()); + + const scalarField out(pos0(Uwp.U() & patch().Sf())); + + valueFraction() = out; + refValue() = p; + refGrad() = (p - pn)*patch().deltaCoeffs(); + + if (internalField().dimensions() == dimPressure) + { + const fvPatchField& rhop = + patch().lookupPatchField(rhoName_); + refValue() *= rhop; + refGrad() *= rhop; + } + + mixedFvPatchScalarField::updateCoeffs(); +} + + +void Foam::wavePressureFvPatchScalarField::write(Ostream& os) const +{ + mixedFvPatchScalarField::write(os); + writeEntryIfDifferent(os, "U", "U", UName_); + writeEntryIfDifferent(os, "rho", "rho", rhoName_); +} + + +// * * * * * * * * * * * * * * Build Macro Function * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField(fvPatchScalarField, wavePressureFvPatchScalarField); +} + +// ************************************************************************* // diff --git a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H new file mode 100644 index 000000000..5b9f44fa9 --- /dev/null +++ b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H @@ -0,0 +1,185 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 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::wavePressureFvPatchScalarField + +Group + grpGenericBoundaryConditions + +Description + This boundary condition provides a wavePressure condition. This sets the + pressure to a value specified by a superposition of wave models. All the + parameters are looked up from the corresponding velocity condition. + +Usage + \table + Property | Description | Req'd? | Default + U | name of the velocity field | no | U + rho | name of the density field | no | rho + \endtable + + Example of the boundary condition specification: + \verbatim + + { + type wavePressure; + U U; + rho rho; + } + \endverbatim + + +SourceFiles + wavePressureFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef wavePressureFvPatchScalarField_H +#define wavePressureFvPatchScalarField_H + +#include "mixedFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class wavePressureFvPatchScalarField Declaration +\*---------------------------------------------------------------------------*/ + +class wavePressureFvPatchScalarField +: + public mixedFvPatchScalarField +{ + // Private data + + //- Name of the velocity field + const word UName_; + + //- Name of the density field + const word rhoName_; + + +public: + + //- Runtime type information + TypeName("wavePressure"); + + + // Constructors + + //- Construct from patch and internal field + wavePressureFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + wavePressureFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given mixedTypeFvPatchField + // onto a new patch + wavePressureFvPatchScalarField + ( + const wavePressureFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + wavePressureFvPatchScalarField + ( + const wavePressureFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new wavePressureFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + wavePressureFvPatchScalarField + ( + const wavePressureFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new wavePressureFvPatchScalarField + ( + *this, + iF + ) + ); + } + + + // Member functions + + // Evaluation functions + + //- Return the current modelled pressure field on the patch faces + tmp p() const; + + //- Return the current modelled pressure field in the neighbour cell + tmp pn() const; + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +#endif + +// ************************************************************************* // diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C index f9f0d7f30..7a6a278f9 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C +++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C @@ -24,10 +24,11 @@ License \*---------------------------------------------------------------------------*/ #include "waveVelocityFvPatchVectorField.H" +#include "wavePressureFvPatchScalarField.H" #include "addToRunTimeSelectionTable.H" #include "levelSet.H" -#include "surfaceFields.H" #include "volFields.H" +#include "fvMeshSubset.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -39,7 +40,10 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField : directionMixedFvPatchVectorField(p, iF), phiName_("phi"), - waves_(db()) + pName_(word::null), + waves_(db()), + faceCellSubset_(nullptr), + faceCellSubsetTimeIndex_(-1) { refValue() = Zero; refGrad() = Zero; @@ -56,7 +60,10 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField : directionMixedFvPatchVectorField(p, iF), phiName_(dict.lookupOrDefault("phi", "phi")), - waves_(db(), dict) + pName_(dict.lookupOrDefault("p", word::null)), + waves_(db(), dict), + faceCellSubset_(nullptr), + faceCellSubsetTimeIndex_(-1) { if (dict.found("value")) { @@ -83,7 +90,10 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField : directionMixedFvPatchVectorField(ptf, p, iF, mapper), phiName_(ptf.phiName_), - waves_(ptf.waves_) + pName_(ptf.pName_), + waves_(ptf.waves_), + faceCellSubset_(nullptr), + faceCellSubsetTimeIndex_(-1) {} @@ -94,7 +104,10 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField : directionMixedFvPatchVectorField(ptf), phiName_(ptf.phiName_), - waves_(ptf.waves_) + pName_(ptf.pName_), + waves_(ptf.waves_), + faceCellSubset_(nullptr), + faceCellSubsetTimeIndex_(-1) {} @@ -106,12 +119,41 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField : directionMixedFvPatchVectorField(ptf, iF), phiName_(ptf.phiName_), - waves_(ptf.waves_) + pName_(ptf.pName_), + waves_(ptf.waves_), + faceCellSubset_(nullptr), + faceCellSubsetTimeIndex_(-1) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +const Foam::fvMeshSubset& +Foam::waveVelocityFvPatchVectorField::faceCellSubset() const +{ + const fvMesh& mesh = patch().boundaryMesh().mesh(); + const label timeIndex = mesh.time().timeIndex(); + + if + ( + !faceCellSubset_.valid() + || (mesh.changing() && faceCellSubsetTimeIndex_ != timeIndex) + ) + { + faceCellSubset_.reset(new fvMeshSubset(mesh)); + faceCellSubset_->setCellSubset(patch().faceCells()); + faceCellSubsetTimeIndex_ = timeIndex; + + // Ask for the tetBasePtIs to trigger all processors to build them. + // Without this, processors that do not contain this patch will + // generate a comms mismatch. + faceCellSubset_->subMesh().tetBasePtIs(); + } + + return faceCellSubset_(); +} + + Foam::tmp Foam::waveVelocityFvPatchVectorField::U() const { const scalar t = db().time().timeOutputValue(); @@ -130,6 +172,45 @@ Foam::tmp Foam::waveVelocityFvPatchVectorField::U() const } +Foam::tmp Foam::waveVelocityFvPatchVectorField::Un() const +{ + const scalar t = db().time().timeOutputValue(); + + const fvMeshSubset& subset = faceCellSubset(); + const fvMesh& meshs = subset.subMesh(); + const label patchis = findIndex(subset.patchMap(), patch().index()); + + const vectorField Us = + levelSetAverage + ( + meshs, + waves_.height(t, meshs.cellCentres())(), + waves_.height(t, meshs.points())(), + waves_.UGas(t, meshs.cellCentres())(), + waves_.UGas(t, meshs.points())(), + waves_.ULiquid(t, meshs.cellCentres())(), + waves_.ULiquid(t, meshs.points())() + ); + + tmp tResult(new vectorField(patch().size())); + vectorField& result = tResult.ref(); + + if (patchis != -1) + { + forAll(meshs.boundary()[patchis], is) + { + const label fs = is + meshs.boundary()[patchis].patch().start(); + const label cs = meshs.boundary()[patchis].faceCells()[is]; + const label f = subset.faceMap()[fs]; + const label i = patch().patch().whichFace(f); + result[i] = Us[cs]; + } + } + + return tResult; +} + + void Foam::waveVelocityFvPatchVectorField::updateCoeffs() { if (updated()) @@ -137,33 +218,60 @@ void Foam::waveVelocityFvPatchVectorField::updateCoeffs() return; } - const vectorField UWave(U()); - - const scalarField& phip = - patch().lookupPatchField(phiName_); - - const scalarField out(pos0(phip)); - - // Where inflow, fix all velocity components to values specified by the - // wave model. - refValue() = (1 - out)*UWave; - valueFraction() = (1 - out)*symmTensor::I; - - // Where outflow, set the normal component of the velocity to a value - // consistent with phi, but scale it to get the volumentic flow rate - // specified by the wave model. Tangential components are extrapolated. - const scalar QPhip = gSum(out*phip); - const scalar QWave = gSum(out*(UWave & patch().Sf())); - const vectorField nBySf(patch().Sf()/sqr(patch().magSf())); - if (QPhip > VSMALL) + if (pName_ != word::null) { - refValue() += out*(QWave/QPhip)*phip*nBySf; + const fvPatchScalarField& pp = + patch().lookupPatchField(pName_); + + if (!isA(pp)) + { + FatalErrorInFunction + << "The corresponding pressure condition for the pressure " + << "field " << pName_ << " on patch " << patch().name() + << " is not of type " + << wavePressureFvPatchScalarField::typeName + << exit(FatalError); + } + + const vectorField U(this->U()), Un(this->Un()); + const scalarField out(pos0(U & patch().Sf())); + + // Where inflow, set all velocity components to values specified by the + // wave model. Where outflow, set the tangential values and the normal + // gradient. + valueFraction() = symmTensor::I - out*sqr(patch().nf()); + refValue() = U; + refGrad() = (U - Un)*patch().deltaCoeffs(); } else { - refValue() += out*QWave*nBySf; + const scalarField& phip = + patch().lookupPatchField(phiName_); + + const vectorField U(this->U()); + const scalarField out(pos0(phip)); + + // Where inflow, fix all velocity components to values specified by the + // wave model. + refValue() = (1 - out)*U; + valueFraction() = (1 - out)*symmTensor::I; + + // Where outflow, set the normal component of the velocity to a value + // consistent with phi, but scale it to get the volumentic flow rate + // specified by the wave model. Tangential components are extrapolated. + const scalar QPhip = gSum(out*phip); + const scalar QWave = gSum(out*(U & patch().Sf())); + const vectorField nBySf(patch().Sf()/sqr(patch().magSf())); + if (QPhip > VSMALL) + { + refValue() += out*(QWave/QPhip)*phip*nBySf; + } + else + { + refValue() += out*QWave*nBySf; + } + valueFraction() += out*sqr(patch().nf()); } - valueFraction() += out*sqr(patch().nf()); directionMixedFvPatchVectorField::updateCoeffs(); directionMixedFvPatchVectorField::evaluate(); @@ -177,6 +285,7 @@ void Foam::waveVelocityFvPatchVectorField::write { directionMixedFvPatchVectorField::write(os); writeEntryIfDifferent(os, "phi", "phi", phiName_); + writeEntryIfDifferent(os, "p", word::null, pName_); waves_.write(os); } diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H index 1e3bec546..9089d902d 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H +++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H @@ -33,6 +33,26 @@ Description corresponding phase fraction condition looks this condition up and re-uses the wave modelling. + Flow reversal will occur in the event that the amplitude of the velocity + oscillation is greater than the mean flow. This triggers special handling, + the form of which depends on whether a pressure field has been specified or + not. + + If a pressure field is not specified, then the proportion of the patch over + which the flow is reversed functions in a manner similar to the + flowRateOutletVelocity condition; i.e., the velocity is extrapolated and + then scaled to match the required outlet flow rate. Numerically, this is + still a fixedValue constraint on the normal velocity, just one which tends + to avoid instability. The corresponding pressure condition should be + fixedFluxPressure. + + If a pressure field is specified, then the normal velocity condition becomes + fixedGradient on outlet faces. This gradient is calculated numerically by + evaluating the wave model on both the patch face and the adjacent cell. The + pressure boundary in this case should be a wavePressure condition. This will + do the opposite; it will fix the pressure value on outlet faces, and the + gradient otherwise. + Usage \table Property | Description | Req'd? | Default @@ -43,6 +63,7 @@ Usage scale | scale factor along the mean flow direction | no | None crossScale | scale factor across the mean flow direction | no | None phi | Name of the flux field | no | phi + p | Name of the pressure field | no | \endtable Example of the boundary condition specification: @@ -91,6 +112,8 @@ SourceFiles namespace Foam { +class fvMeshSubset; + /*---------------------------------------------------------------------------*\ Class waveVelocityFvPatchVectorField Declaration \*---------------------------------------------------------------------------*/ @@ -104,9 +127,18 @@ class waveVelocityFvPatchVectorField //- Name of the flux field const word phiName_; + //- Name of the pressure field + const word pName_; + //- Wave superposition const waveSuperposition waves_; + //- Mesh subset corresponding to the patch adjacent cells + mutable autoPtr faceCellSubset_; + + //- Time index for keeping the subset up to date + mutable label faceCellSubsetTimeIndex_; + public: @@ -179,18 +211,30 @@ public: // Access + //- Access the name of the pressure field + const word& pName() const + { + return pName_; + } + //- Access the wave models const waveSuperposition& waves() const { return waves_; } + //- Access the face-cell subset + const fvMeshSubset& faceCellSubset() const; + // Evaluation functions - //- Return the current modelled velocity field + //- Return the current modelled velocity field on the patch faces tmp U() const; + //- Return the current modelled velocity field in the neighbour cell + tmp Un() const; + //- Update the coefficients associated with the patch field virtual void updateCoeffs();