From 8144cfb5160c1386fd457d097f8ed0d945455bed Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 19 Apr 2018 10:45:12 +0100 Subject: [PATCH] waves: Control additions An "inletOutlet" switch has been added to the wave velocity boundary condition to allow the boundary to be fixed, as is possible for the corresponding alpha condition. A "heightAboveWave" option has been added to the wave superposition class to calculate velocity based on the height above the wave, rather than above the origin. This may improve initialisation but it may also generate divergence in the initial velocity field. The alpha condition has also been completed so that it applies a modelled gradient when the flow points out and a wave pressure condition is in use. --- .../waveAlpha/waveAlphaFvPatchScalarField.C | 91 +++++++++++++++++-- .../waveAlpha/waveAlphaFvPatchScalarField.H | 29 ++++-- .../wavePressureFvPatchScalarField.C | 3 +- .../wavePressureFvPatchScalarField.H | 11 ++- .../waveVelocityFvPatchVectorField.C | 77 ++++++++-------- .../waveVelocityFvPatchVectorField.H | 47 ++++++---- .../waveSuperposition/waveSuperposition.C | 29 +++++- .../waveSuperposition/waveSuperposition.H | 4 + 8 files changed, 217 insertions(+), 74 deletions(-) diff --git a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C index 1c240eb6a..519d46c11 100644 --- a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C +++ b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,11 +24,13 @@ License \*---------------------------------------------------------------------------*/ #include "waveAlphaFvPatchScalarField.H" +#include "wavePressureFvPatchScalarField.H" #include "waveVelocityFvPatchVectorField.H" #include "addToRunTimeSelectionTable.H" #include "levelSet.H" #include "surfaceFields.H" #include "volFields.H" +#include "fvMeshSubset.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -139,6 +141,50 @@ Foam::tmp Foam::waveAlphaFvPatchScalarField::alpha() const } +Foam::tmp Foam::waveAlphaFvPatchScalarField::alphan() 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 alphas + ( + levelSetFraction + ( + meshs, + Up.waves().height(t, meshs.cellCentres())(), + Up.waves().height(t, meshs.points())(), + !liquid_ + ) + ); + + 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] = alphas[cs]; + } + } + + return tResult; +} + + void Foam::waveAlphaFvPatchScalarField::updateCoeffs() { if (updated()) @@ -146,18 +192,49 @@ void Foam::waveAlphaFvPatchScalarField::updateCoeffs() return; } - refValue() = alpha(); + const fvPatchVectorField& Up = + patch().lookupPatchField(UName_); - if (inletOutlet_) + if (!isA(Up)) { - const scalarField& phip = - patch().lookupPatchField("phi"); + FatalErrorInFunction + << "The corresponding condition for the velocity " + << "field " << UName_ << " on patch " << patch().name() + << " is not of type " << waveVelocityFvPatchVectorField::typeName + << exit(FatalError); + } - valueFraction() = 1 - pos0(phip); + const waveVelocityFvPatchVectorField& Uwp = + refCast(Up); + + const fvPatchScalarField& pp = + patch().lookupPatchField(Uwp.pName()); + + if (isA(pp)) + { + const scalarField alpha(this->alpha()), alphan(this->alphan()); + const scalarField out(pos0(Uwp.U() & patch().Sf())); + + valueFraction() = out; + refValue() = alpha; + refGrad() = (alpha - alphan)*patch().deltaCoeffs(); } else { - valueFraction() = 1; + refValue() = alpha(); + + if (inletOutlet_) + { + const scalarField& phip = + patch().lookupPatchField("phi"); + const scalarField out(pos0(phip)); + + valueFraction() = 1 - out; + } + else + { + valueFraction() = 1; + } } mixedFvPatchScalarField::updateCoeffs(); diff --git a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H index 049804d64..58d7c519c 100644 --- a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H +++ b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,6 +32,21 @@ Description fraction to that specified by a superposition of wave models. All the parameters are looked up from the corresponding velocity condition. + 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 the inletOutlet flag and whether a wave + pressure condition is being used. + + If a wave pressure condition is not being used, the inletOutlet switches + between a fixedValue and an inletOutlet condition, with the value given by + the wave model. If fixedValue, the result may be more accurate, but it + might also be unstable. + + If a wave pressure condition is being used, then the normal phase fraction + 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. + Usage \table Property | Description | Req'd? | Default @@ -76,16 +91,12 @@ class waveAlphaFvPatchScalarField // Private data //- Name of the velocity field - word UName_; + const word UName_; //- Is this alpha field that of the liquid under the wave? const Switch liquid_; - //- Act as an inlet/outlet patch? If false, the alpha field will be set - // by the wave model even during outflow. This may apply the wave model - // more accurately, but it might also be unstable. If true, the alpha - // boundary condition will be zero gradient when the flow reverses, as - // is usual. + //- Act as an inlet/outlet patch? const Switch inletOutlet_; @@ -172,6 +183,10 @@ public: //- Return the current modelled phase fraction field tmp alpha() const; + //- Return the current modelled phase fraction field in the + // neighbour cell + tmp alphan() const; + //- Update the coefficients associated with the patch field virtual void updateCoeffs(); diff --git a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C index a75eb70c4..a1dbd4c2c 100644 --- a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C +++ b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -216,7 +216,6 @@ void Foam::wavePressureFvPatchScalarField::updateCoeffs() } const scalarField p(this->p()), pn(this->pn()); - const scalarField out(pos0(Uwp.U() & patch().Sf())); valueFraction() = out; diff --git a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H index 5b9f44fa9..1055fe033 100644 --- a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H +++ b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,6 +32,15 @@ Description pressure to a value specified by a superposition of wave models. All the parameters are looked up from the corresponding velocity condition. + This functions like an outletInlet condition. Faces on which the flow is + leaving the domain have a value set by the wave model. Faces on which the + flow is entering the domain have the gradient set. This gradient is + calculated numerically by evaluating the wave model on both the patch face + and the adjacent cell. + + Use of this boundary condition triggers a consistent behaviour in the + corresponding velocity and phase-fraction conditions. + Usage \table Property | Description | Req'd? | Default diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C index 44ef9f725..608ab8cdc 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C +++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C @@ -40,7 +40,8 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField : directionMixedFvPatchVectorField(p, iF), phiName_("phi"), - pName_(word::null), + pName_("p"), + inletOutlet_(true), waves_(db()), faceCellSubset_(nullptr), faceCellSubsetTimeIndex_(-1) @@ -60,7 +61,8 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField : directionMixedFvPatchVectorField(p, iF), phiName_(dict.lookupOrDefault("phi", "phi")), - pName_(dict.lookupOrDefault("p", word::null)), + pName_(dict.lookupOrDefault("p", "p")), + inletOutlet_(dict.lookupOrDefault("inletOutlet", true)), waves_(db(), dict), faceCellSubset_(nullptr), faceCellSubsetTimeIndex_(-1) @@ -91,6 +93,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField directionMixedFvPatchVectorField(ptf, p, iF, mapper), phiName_(ptf.phiName_), pName_(ptf.pName_), + inletOutlet_(ptf.inletOutlet_), waves_(ptf.waves_), faceCellSubset_(nullptr), faceCellSubsetTimeIndex_(-1) @@ -105,6 +108,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField directionMixedFvPatchVectorField(ptf), phiName_(ptf.phiName_), pName_(ptf.pName_), + inletOutlet_(ptf.inletOutlet_), waves_(ptf.waves_), faceCellSubset_(nullptr), faceCellSubsetTimeIndex_(-1) @@ -120,6 +124,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField directionMixedFvPatchVectorField(ptf, iF), phiName_(ptf.phiName_), pName_(ptf.pName_), + inletOutlet_(ptf.inletOutlet_), waves_(ptf.waves_), faceCellSubset_(nullptr), faceCellSubsetTimeIndex_(-1) @@ -220,21 +225,11 @@ void Foam::waveVelocityFvPatchVectorField::updateCoeffs() return; } - if (pName_ != word::null) + const fvPatchScalarField& pp = + patch().lookupPatchField(pName_); + + if (isA(pp)) { - 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())); @@ -247,32 +242,41 @@ void Foam::waveVelocityFvPatchVectorField::updateCoeffs() } else { - 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) + if (inletOutlet_) { - refValue() += out*(QWave/QPhip)*phip*nBySf; + 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)*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 volumetric + // 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()); } else { - refValue() += out*QWave*nBySf; + refValue() = U; + valueFraction() = symmTensor::I; } - valueFraction() += out*sqr(patch().nf()); } directionMixedFvPatchVectorField::updateCoeffs(); @@ -287,7 +291,8 @@ void Foam::waveVelocityFvPatchVectorField::write { directionMixedFvPatchVectorField::write(os); writeEntryIfDifferent(os, "phi", "phi", phiName_); - writeEntryIfDifferent(os, "p", word::null, pName_); + writeEntryIfDifferent(os, "p", "p", pName_); + writeEntryIfDifferent(os, "inletOutlet", true, inletOutlet_); waves_.write(os); } diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H index eff095204..aba710d79 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H +++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H @@ -30,32 +30,41 @@ Group Description This boundary condition provides a waveVelocity condition. This sets the velocity to that specified by a superposition of wave models. The - corresponding phase fraction condition looks this condition up and re-uses - the wave modelling. + corresponding phase fraction and pressure conditions look this condition up + and re-use 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. + the form of which depends on the inletOutlet flag and whether a wave + pressure condition is being used. - 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 + If a wave pressure condition is not being used, and inletOutlet is false, + then this is a standard fixed value condition, with the value supplied by + the wave model. If flow reversal occurs this state may be unstable. The + corresponding pressure condition should be fixedFluxPressure. + + If a wave pressure condition is not being used, and inletOutlet is true or + 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. Again, 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. + If a wave pressure condition is being used, 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 + phi | Name of the flux field | no | phi + p | Name of the pressure field | no | p + inletOutlet | does the condition behave like inletOutlet | no | true origin | origin of the wave coordinate system | yes | direction | direction of the mean flow | yes | speed | speed of the mean flow | yes | @@ -63,8 +72,7 @@ Usage ramp | ramping function for the mean flow speed | no | None 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 | + heightAboveWave | use with the height above the wave | no | false \endtable Example of the boundary condition specification: @@ -132,6 +140,9 @@ class waveVelocityFvPatchVectorField //- Name of the pressure field const word pName_; + //- Act as an inlet/outlet patch? + const Switch inletOutlet_; + //- Wave superposition const waveSuperposition waves_; diff --git a/src/waves/waveSuperposition/waveSuperposition.C b/src/waves/waveSuperposition/waveSuperposition.C index 62b9a4731..2e3fe8e3b 100644 --- a/src/waves/waveSuperposition/waveSuperposition.C +++ b/src/waves/waveSuperposition/waveSuperposition.C @@ -182,7 +182,8 @@ Foam::waveSuperposition::waveSuperposition(const objectRegistry& db) waveAngles_(), ramp_(), scale_(), - crossScale_() + crossScale_(), + heightAboveWave_(false) {} @@ -196,7 +197,8 @@ Foam::waveSuperposition::waveSuperposition(const waveSuperposition& waves) waveAngles_(waves.waveAngles_), ramp_(waves.ramp_, false), scale_(waves.scale_, false), - crossScale_(waves.crossScale_, false) + crossScale_(waves.crossScale_, false), + heightAboveWave_(false) {} @@ -229,7 +231,8 @@ Foam::waveSuperposition::waveSuperposition dict.found("crossScale") ? Function1::New("crossScale", dict) : autoPtr>() - ) + ), + heightAboveWave_(dict.lookupOrDefault("heightAboveWave", false)) { const PtrList waveEntries(dict.lookup("waves")); @@ -285,6 +288,11 @@ Foam::tmp Foam::waveSuperposition::ULiquid vectorField xyz(p.size()); transformation(p, axes, u, xyz); + if (heightAboveWave_) + { + xyz.replace(2, height(t, p)); + } + return UMean(t) + (velocity(t, xyz) & axes); } @@ -302,6 +310,11 @@ Foam::tmp Foam::waveSuperposition::UGas axes = tensor(- axes.x(), - axes.y(), axes.z()); + if (heightAboveWave_) + { + xyz.replace(2, height(t, p)); + } + return UMean(t) + (velocity(t, xyz) & axes); } @@ -317,6 +330,11 @@ Foam::tmp Foam::waveSuperposition::pLiquid vectorField xyz(p.size()); transformation(p, axes, u, xyz); + if (heightAboveWave_) + { + xyz.replace(2, height(t, p)); + } + return pressure(t, xyz); } @@ -358,6 +376,11 @@ void Foam::waveSuperposition::write(Ostream& os) const { crossScale_->writeData(os); } + if (heightAboveWave_) + { + os.writeKeyword("heightAboveWave") << heightAboveWave_ + << token::END_STATEMENT << nl; + } } diff --git a/src/waves/waveSuperposition/waveSuperposition.H b/src/waves/waveSuperposition/waveSuperposition.H index 4bb0e638c..6c6b154dc 100644 --- a/src/waves/waveSuperposition/waveSuperposition.H +++ b/src/waves/waveSuperposition/waveSuperposition.H @@ -78,6 +78,10 @@ class waveSuperposition //- Scaling perpendicular to the flow direction const autoPtr> crossScale_; + //- Calculate wave properties using the height above the wave (true) or + // the height above the origin (false)? + const Switch heightAboveWave_; + // Private Member Functions