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