waves: Pressure-velocity boundary formulation

A wavePressure boundary condition has been added, and the Airy-type wave
models have been extended to generate the unsteady pressure field. This
provides another option for specifying wave motion at a boundary.

If a waveVelocity condition is used in isolation, then any outlet flow
will be extrapolated and scaled to match the required flow rate. This is
similar to how a flowRateOutletVelocity condition works.

0/U:

    <patchName>
    {
        type        waveVelocity;
        // wave parameters ...
    }

0/p_rgh:

    <patchName>
    {
        type        fixedFluxPressure;
    }

If a waveVelocity is used in conjunction with the new wavePressure
condition, then one will set the value and the other the gradient, as
appropriate for the direction of the flow.

0/U:

    <patchName>
    {
        type        waveVelocity;
        // wave parameters ...
        p           p_rgh;
    }

0/p_rgh:

    <patchName>
    {
        type        wavePressure;
    }

This new pressure-velocity formulation is less stable, but generates
more accurate waveforms on patches where the velocity reverses. It is
also necessary for sub-surface cases where fixing the velocity around
the entire domain generates a continuity error.

This work was supported by Alice Gillespie, on behalf of M3 Wave
This commit is contained in:
Will Bainbridge
2017-09-28 08:50:35 +01:00
parent 63c18662dd
commit d9fc7eb887
6 changed files with 623 additions and 31 deletions

View File

@ -9,5 +9,6 @@ waveSuperposition/waveSuperposition.C
derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C
derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C
derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C
LIB = $(FOAM_LIBBIN)/libwaves

View File

@ -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

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(p, iF),
UName_("U"),
rhoName_("rho")
{
refValue() = Zero;
refGrad() = Zero;
valueFraction() = Zero;
}
Foam::wavePressureFvPatchScalarField::wavePressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
mixedFvPatchScalarField(p, iF),
UName_(dict.lookupOrDefault<word>("U", "U")),
rhoName_(dict.lookupOrDefault<word>("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<scalar, volMesh>& 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<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(ptf, iF),
UName_(ptf.UName_),
rhoName_(ptf.rhoName_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::wavePressureFvPatchScalarField::p() const
{
const waveVelocityFvPatchVectorField& Up =
refCast<const waveVelocityFvPatchVectorField>
(
patch().lookupPatchField<volVectorField, scalar>(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::scalarField> Foam::wavePressureFvPatchScalarField::pn() const
{
const waveVelocityFvPatchVectorField& Up =
refCast<const waveVelocityFvPatchVectorField>
(
patch().lookupPatchField<volVectorField, scalar>(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<scalarField> 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<volVectorField, scalar>(UName_);
if (!isA<waveVelocityFvPatchVectorField>(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<const waveVelocityFvPatchVectorField>(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<scalar>& rhop =
patch().lookupPatchField<volScalarField, scalar>(rhoName_);
refValue() *= rhop;
refGrad() *= rhop;
}
mixedFvPatchScalarField::updateCoeffs();
}
void Foam::wavePressureFvPatchScalarField::write(Ostream& os) const
{
mixedFvPatchScalarField::write(os);
writeEntryIfDifferent<word>(os, "U", "U", UName_);
writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
}
// * * * * * * * * * * * * * * Build Macro Function * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField(fvPatchScalarField, wavePressureFvPatchScalarField);
}
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
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
<patchName>
{
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<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
wavePressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given mixedTypeFvPatchField
// onto a new patch
wavePressureFvPatchScalarField
(
const wavePressureFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
wavePressureFvPatchScalarField
(
const wavePressureFvPatchScalarField&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField>
(
new wavePressureFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
wavePressureFvPatchScalarField
(
const wavePressureFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchScalarField> clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchScalarField>
(
new wavePressureFvPatchScalarField
(
*this,
iF
)
);
}
// Member functions
// Evaluation functions
//- Return the current modelled pressure field on the patch faces
tmp<scalarField> p() const;
//- Return the current modelled pressure field in the neighbour cell
tmp<scalarField> pn() const;
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -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<word>("phi", "phi")),
waves_(db(), dict)
pName_(dict.lookupOrDefault<word>("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::vectorField> Foam::waveVelocityFvPatchVectorField::U() const
{
const scalar t = db().time().timeOutputValue();
@ -130,6 +172,45 @@ Foam::tmp<Foam::vectorField> Foam::waveVelocityFvPatchVectorField::U() const
}
Foam::tmp<Foam::vectorField> 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<vectorField> 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<surfaceScalarField, scalar>(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<volScalarField, scalar>(pName_);
if (!isA<wavePressureFvPatchScalarField>(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<surfaceScalarField, scalar>(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<word>(os, "phi", "phi", phiName_);
writeEntryIfDifferent<word>(os, "p", word::null, pName_);
waves_.write(os);
}

View File

@ -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<fvMeshSubset> 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<vectorField> U() const;
//- Return the current modelled velocity field in the neighbour cell
tmp<vectorField> Un() const;
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();