diff --git a/applications/utilities/preProcessing/setWaves/setWaves.C b/applications/utilities/preProcessing/setWaves/setWaves.C index 895b3a990..c0885d845 100644 --- a/applications/utilities/preProcessing/setWaves/setWaves.C +++ b/applications/utilities/preProcessing/setWaves/setWaves.C @@ -30,111 +30,15 @@ Description #include "fvCFD.H" #include "levelSet.H" +#include "pointFields.H" #include "timeSelector.H" +#include "wallDist.H" #include "waveAlphaFvPatchScalarField.H" #include "waveVelocityFvPatchVectorField.H" #include "waveSuperposition.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -void addWaves -( - const waveSuperposition& waves, - const bool liquid, - volScalarField& alpha, - volVectorField& U -) -{ - const scalar t = alpha.time().value();; - const fvMesh& mesh = alpha.mesh(); - const pointMesh& pMesh = pointMesh::New(mesh); - - // Height fields - const scalarField heightC(waves.height(t, mesh.cellCentres())); - const scalarField heightP(waves.height(t, mesh.points())); - - // Velocity fields - const DimensionedField - UGasC - ( - IOobject("UGasC", mesh.time().timeName(), mesh), - mesh, - dimVelocity, - waves.UGas(t, mesh.cellCentres()) - ); - const DimensionedField - UGasP - ( - IOobject("UGasP", mesh.time().timeName(), mesh), - pMesh, - dimVelocity, - waves.UGas(t, mesh.points()) - ); - const DimensionedField - ULiquidC - ( - IOobject("ULiquidC", mesh.time().timeName(), mesh), - mesh, - dimVelocity, - waves.ULiquid(t, mesh.cellCentres()) - ); - const DimensionedField - ULiquidP - ( - IOobject("ULiquidP", mesh.time().timeName(), mesh), - pMesh, - dimVelocity, - waves.ULiquid(t, mesh.points()) - ); - - // Convert from the level set to volume-averaged fields and sum up - alpha.ref() += levelSetFraction(mesh, heightC, heightP, !liquid); - U.ref() += - levelSetAverage - ( - mesh, - heightC, - heightP, - UGasC, - UGasP, - ULiquidC, - ULiquidP - ); - - // Now set the boundary fields - forAll(alpha.boundaryField(), patchi) - { - fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi]; - fvPatchVectorField& Up = U.boundaryFieldRef()[patchi]; - - const fvPatch& patch = alphap.patch(); - - // Height fields - const scalarField heightF(waves.height(t, patch.Cf())); - const scalarField heightP(waves.height(t, patch.patch().localPoints())); - - // Velocity fields - const vectorField UGasC(waves.UGas(t, mesh.cellCentres())); - const vectorField UGasP(waves.UGas(t, mesh.points())); - const vectorField ULiquidC(waves.ULiquid(t, mesh.cellCentres())); - const vectorField ULiquidP(waves.ULiquid(t, mesh.points())); - - alphap == alphap + levelSetFraction(patch, heightF, heightP, !liquid); - Up == Up - + levelSetAverage - ( - patch, - heightC, - heightP, - UGasC, - UGasP, - ULiquidC, - ULiquidP - ); - } -} - - int main(int argc, char *argv[]) { timeSelector::addOptions(false, false); @@ -161,6 +65,8 @@ int main(int argc, char *argv[]) #include "createMesh.H" #include "readGravitationalAcceleration.H" + const pointMesh& pMesh = pointMesh::New(mesh); + forAll(timeDirs, timeI) { runTime.setTime(timeDirs[timeI], timeI); @@ -169,7 +75,7 @@ int main(int argc, char *argv[]) mesh.readUpdate(); - // Read the phase fraction and velocity fields + // Read the fields which are to be set volScalarField alpha ( IOobject @@ -193,20 +99,58 @@ int main(int argc, char *argv[]) mesh ); - // Zero the fields - alpha = dimensionedScalar("0", alpha.dimensions(), 0); - U = dimensionedVector("0", U.dimensions(), vector::zero); - forAll(alpha.boundaryField(), patchi) - { - alpha.boundaryFieldRef()[patchi] == 0; - U.boundaryFieldRef()[patchi] == vector::zero; - } + // Create modelled fields on both cells and points + volScalarField height + ( + IOobject("height", runTime.timeName(), mesh), + mesh, + dimensionedScalar("0", dimLength, 0) + ); + pointScalarField heightp + ( + IOobject("heightp", runTime.timeName(), mesh), + pMesh, + dimensionedScalar("0", dimLength, 0) + ); + volVectorField uGas + ( + IOobject("uGas", runTime.timeName(), mesh), + mesh, + dimensionedVector("0", dimVelocity, vector::zero) + ); + pointVectorField uGasp + ( + IOobject("uGasp", runTime.timeName(), mesh), + pMesh, + dimensionedVector("0", dimLength, vector::zero) + ); + volVectorField uLiquid + ( + IOobject("uLiquid", runTime.timeName(), mesh), + mesh, + dimensionedVector("0", dimVelocity, vector::zero) + ); + pointVectorField uLiquidp + ( + IOobject("uLiquidp", runTime.timeName(), mesh), + pMesh, + dimensionedVector("0", dimLength, vector::zero) + ); - // Loop the patches, looking for wave conditions - forAll(alpha.boundaryField(), patchi) + // The number of wave patches + label nWaves = 0; + + // Whether the alpha conditions refer to the liquid phase + bool liquid = false; + + // The mean velocity of one of the wave patches + vector UMeanp = vector::zero; + + // Loop the patches, averaging and superimposing wave model data + forAll(mesh.boundary(), patchi) { - const fvPatchScalarField& alphap = alpha.boundaryField()[patchi]; - const fvPatchVectorField& Up = U.boundaryField()[patchi]; + fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi]; + fvPatchVectorField& Up = U.boundaryFieldRef()[patchi]; const bool isWave = isA(alphap); @@ -223,16 +167,175 @@ int main(int argc, char *argv[]) << " and vice-versa." << exit(FatalError); } - if (isWave) + if (!isWave) { - Info<< "Adding waves from patch " << Up.patch().name() << endl; - addWaves + continue; + } + + Info<< "Adding waves from patch " << Up.patch().name() << endl; + + const waveSuperposition& waves = + refCast(Up).waves(); + + UMeanp = waves.UMean(); + + const bool liquidp = + refCast(alphap).liquid(); + if (nWaves > 0 && liquidp != liquid) + { + FatalErrorInFunction + << "All " << waveAlphaFvPatchScalarField::typeName + << "patch fields must be configured for the same phase," + << " i.e., the liquid switch must have the same value." + << exit(FatalError); + } + liquid = liquidp; + + const scalar t = runTime.value(); + const pointField& ccs = mesh.cellCentres(); + const pointField& pts = mesh.points(); + + // Internal field superposition + height.primitiveFieldRef() += waves.height(t, ccs); + heightp.primitiveFieldRef() += waves.height(t, pts); + uGas.primitiveFieldRef() += waves.UGas(t, ccs) - UMeanp; + uGasp.primitiveFieldRef() += waves.UGas(t, pts) - UMeanp; + uLiquid.primitiveFieldRef() += waves.ULiquid(t, ccs) - UMeanp; + uLiquidp.primitiveFieldRef() += waves.ULiquid(t, pts) - UMeanp; + + // Boundary field superposition + forAll(mesh.boundary(), patchj) + { + const pointField& fcs = mesh.boundary()[patchj].Cf(); + height.boundaryFieldRef()[patchj] += waves.height(t, fcs); + uGas.boundaryFieldRef()[patchj] += waves.UGas(t, fcs) - UMeanp; + uLiquid.boundaryFieldRef()[patchj] += + waves.ULiquid(t, fcs) - UMeanp; + } + + ++ nWaves; + } + + // Warn and skip to the next time if no wave patches were found + if (nWaves == 0) + { + WarningInFunction + << "No " << waveAlphaFvPatchScalarField::typeName << " or " + << waveVelocityFvPatchVectorField::typeName << " patch fields " + << "were found. No waves have been set." << endl; + + continue; + } + + // Create the mean velocity field + volVectorField UMean + ( + IOobject("UMean", runTime.timeName(), mesh), + mesh, + dimensionedVector("UMean", dimVelocity, UMeanp) + ); + + if (nWaves > 1) + { + // Create weighted average fields for the mean velocity + volScalarField weight + ( + IOobject("weight", runTime.timeName(), mesh), + mesh, + dimensionedScalar("0", dimless/dimLength, 0) + ); + volVectorField weightUMean + ( + IOobject("weightUMean", runTime.timeName(), mesh), + mesh, + dimensionedVector("0", dimVelocity/dimLength, vector::zero) + ); + + // Loop the patches, inverse-distance weighting the mean velocities + forAll(mesh.boundary(), patchi) + { + fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi]; + fvPatchVectorField& Up = U.boundaryFieldRef()[patchi]; + + const bool isWave = isA(alphap); + + if (!isWave) + { + continue; + } + + const waveSuperposition& waves = + refCast(Up).waves(); + + UMeanp = waves.UMean(); + + const volScalarField w ( - refCast(Up).waves(), - refCast(alphap).liquid(), - alpha, - U + 1 + /( + wallDist(mesh, labelList(1, patchi)).y() + + dimensionedScalar("ySmall", dimLength, SMALL) + ) ); + weight += w; + weightUMean += + w*dimensionedVector("UMeanp", dimVelocity, UMeanp); + } + + // Complete the average for the mean velocity + UMean = weightUMean/weight; + } + + // Set the internal fields + alpha.ref() = levelSetFraction(mesh, height, heightp, !liquid); + U.ref() = + UMean + + levelSetAverage + ( + mesh, + height, + heightp, + uGas, + uGasp, + uLiquid, + uLiquidp + ); + + // Set the boundary fields + forAll(mesh.boundary(), patchi) + { + fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi]; + fvPatchVectorField& Up = U.boundaryFieldRef()[patchi]; + + const bool isWave = isA(alphap); + + if (!isWave) + { + alphap == + levelSetFraction + ( + mesh.boundary()[patchi], + height.boundaryField()[patchi], + heightp.boundaryField()[patchi].patchInternalField(), + !liquid + ); + Up == + UMean.boundaryField()[patchi] + + levelSetAverage + ( + mesh.boundary()[patchi], + height.boundaryField()[patchi], + heightp.boundaryField()[patchi].patchInternalField()(), + uGas.boundaryField()[patchi], + uGasp.boundaryField()[patchi].patchInternalField()(), + uLiquid.boundaryField()[patchi], + uLiquidp.boundaryField()[patchi].patchInternalField()() + ); + } + else + { + alphap == refCast(alphap).alpha(); + Up == refCast(Up).U(); } } diff --git a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C index de3ea9fc7..5466911c6 100644 --- a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C +++ b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C @@ -30,30 +30,6 @@ License #include "surfaceFields.H" #include "volFields.H" -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -Foam::tmp Foam::waveAlphaFvPatchScalarField::alpha() const -{ - const scalar t = db().time().timeOutputValue(); - - const waveVelocityFvPatchVectorField& Up = - refCast - ( - patch().lookupPatchField(UName_) - ); - const waveSuperposition& waves = Up.waves(); - - return - levelSetFraction - ( - patch(), - waves.height(t, patch().Cf()), - waves.height(t, patch().patch().localPoints()), - !liquid_ - ); -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::waveAlphaFvPatchScalarField::waveAlphaFvPatchScalarField @@ -141,6 +117,28 @@ Foam::waveAlphaFvPatchScalarField::waveAlphaFvPatchScalarField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +Foam::tmp Foam::waveAlphaFvPatchScalarField::alpha() const +{ + const scalar t = db().time().timeOutputValue(); + + const waveVelocityFvPatchVectorField& Up = + refCast + ( + patch().lookupPatchField(UName_) + ); + const waveSuperposition& waves = Up.waves(); + + return + levelSetFraction + ( + patch(), + waves.height(t, patch().Cf()), + waves.height(t, patch().patch().localPoints()), + !liquid_ + ); +} + + void Foam::waveAlphaFvPatchScalarField::updateCoeffs() { if (updated()) diff --git a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H index a22499270..c97b49547 100644 --- a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H +++ b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H @@ -89,12 +89,6 @@ class waveAlphaFvPatchScalarField const bool inletOutlet_; - // Private Member Functions - - //- Return the current modelled phase fraction field - tmp alpha() const; - - public: //- Runtime type information @@ -175,6 +169,9 @@ public: // Evaluation functions + //- Return the current modelled phase fraction field + tmp alpha() const; + //- Update the coefficients associated with the patch field virtual void updateCoeffs(); diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C index 22da4097b..c6877d3df 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C +++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C @@ -29,26 +29,6 @@ License #include "surfaceFields.H" #include "volFields.H" -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -Foam::tmp Foam::waveVelocityFvPatchVectorField::U() const -{ - const scalar t = db().time().timeOutputValue(); - - return - levelSetAverage - ( - patch(), - waves_.height(t, patch().Cf()), - waves_.height(t, patch().patch().localPoints()), - waves_.UGas(t, patch().Cf())(), - waves_.UGas(t, patch().patch().localPoints())(), - waves_.ULiquid(t, patch().Cf())(), - waves_.ULiquid(t, patch().patch().localPoints())() - ); -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField @@ -127,6 +107,24 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +Foam::tmp Foam::waveVelocityFvPatchVectorField::U() const +{ + const scalar t = db().time().timeOutputValue(); + + return + levelSetAverage + ( + patch(), + waves_.height(t, patch().Cf()), + waves_.height(t, patch().patch().localPoints()), + waves_.UGas(t, patch().Cf())(), + waves_.UGas(t, patch().patch().localPoints())(), + waves_.ULiquid(t, patch().Cf())(), + waves_.ULiquid(t, patch().patch().localPoints())() + ); +} + + void Foam::waveVelocityFvPatchVectorField::updateCoeffs() { if (updated()) diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H index 2d259703e..d0ce0e4c5 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H +++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H @@ -104,12 +104,6 @@ class waveVelocityFvPatchVectorField const waveSuperposition waves_; - // Private Member Functions - - //- Return the current modelled velocity field - tmp U() const; - - public: //- Runtime type information @@ -190,6 +184,9 @@ public: // Evaluation functions + //- Return the current modelled velocity field + tmp U() const; + //- Update the coefficients associated with the patch field virtual void updateCoeffs(); diff --git a/src/waves/waveSuperposition/waveSuperposition.C b/src/waves/waveSuperposition/waveSuperposition.C index ca762e68c..67c068acd 100644 --- a/src/waves/waveSuperposition/waveSuperposition.C +++ b/src/waves/waveSuperposition/waveSuperposition.C @@ -245,7 +245,7 @@ Foam::tmp Foam::waveSuperposition::ULiquid vectorField xyz(p.size()); transformation(p, axes, u, xyz); - return direction_*speed_ + (velocity(t, xyz) & axes); + return UMean() + (velocity(t, xyz) & axes); } @@ -262,7 +262,7 @@ Foam::tmp Foam::waveSuperposition::UGas axes = tensor(- axes.x(), - axes.y(), axes.z()); - return direction_*speed_ + (velocity(t, xyz) & axes); + return UMean() + (velocity(t, xyz) & axes); } diff --git a/src/waves/waveSuperposition/waveSuperposition.H b/src/waves/waveSuperposition/waveSuperposition.H index 099a657df..b2fefe788 100644 --- a/src/waves/waveSuperposition/waveSuperposition.H +++ b/src/waves/waveSuperposition/waveSuperposition.H @@ -139,6 +139,12 @@ public: //- Get the gas velocity at a given time and global positions tmp UGas(const scalar t, const vectorField& p) const; + //- Get the mean flow velocity + inline vector UMean() const + { + return direction_*speed_; + } + //- Write void write(Ostream&) const; };