waves: Split mean flow from wave perturbation modelling

In order to increase the flexibility of the wave library, the mean flow
handling has been removed from the waveSuperposition class. This makes
waveSuperposition work purely in terms of perturbations to a mean
background flow.

The input has also been split, with waves now defined as region-wide
settings in constant/waveProperties. The mean flow parameters are sill
defined by the boundary conditions.

The new format of the velocity boundary is much simpler. Only a mean
flow velocity is required.

    In 0/U:

        boundaryField
        {
            inlet
            {
                type            waveVelocity;
                UMean           (2 0 0);
            }
            // etc ...
        }

Other wave boundary conditions have not changed.

The constant/waveProperties file contains the wave model selections and
the settings to define the associated coordinate system and scaling
functions:

    In constant/waveProperties:

        origin          (0 0 0);
        direction       (1 0 0);
        waves
        (
            Airy
            {
                length      300;
                amplitude   2.5;
                phase       0;
                angle       0;
            }
        );
        scale           table ((1200 1) (1800 0));
        crossScale      constant 1;

setWaves has been changed to use a system/setWavesDict file rather than
relying on command-line arguments. It also now requires a mean velocity
to be specified in order to prevent ambiguities associated with multiple
inlet patches. An example is shown below:

    In system/setWavesDict:

        alpha   alpha.water;
        U       U;
        liquid  true;
        UMean   (1 0 0);
This commit is contained in:
Will Bainbridge
2018-12-06 12:05:07 +00:00
parent bd2f275e09
commit 967edc9425
25 changed files with 490 additions and 473 deletions

View File

@ -30,14 +30,17 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "argList.H"
#include "levelSet.H"
#include "pointFields.H"
#include "timeSelector.H"
#include "uniformDimensionedFields.H"
#include "volFields.H"
#include "wallDist.H"
#include "waveAlphaFvPatchScalarField.H"
#include "waveVelocityFvPatchVectorField.H"
#include "waveSuperposition.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -45,18 +48,27 @@ int main(int argc, char *argv[])
{
timeSelector::addOptions(false, false);
Foam::argList::addOption
#include "addDictOption.H"
#include "addRegionOption.H"
argList::addOption
(
"alpha",
"name",
"name of the volume fraction field, default is \"alpha\""
);
argList::addOption
(
"U",
"name",
"name of the velocity field, default is \"U\""
);
Foam::argList::addOption
argList::addBoolOption
(
"alpha",
"name",
"name of the volume fraction field, default is \"alpha\""
"gas",
"the volume fraction field is that that of the gas above the wave"
);
#include "setRootCase.H"
@ -64,11 +76,31 @@ int main(int argc, char *argv[])
instantList timeDirs = timeSelector::selectIfPresent(runTime, args);
#include "createMesh.H"
#include "createNamedMesh.H"
const word dictName("setWavesDict");
#include "setSystemMeshDictionaryIO.H"
Info<< "Reading " << dictName << "\n" << endl;
IOdictionary setWavesDict(dictIO);
#include "readGravitationalAcceleration.H"
const pointMesh& pMesh = pointMesh::New(mesh);
// Parse options
const word alphaName = setWavesDict.lookupOrDefault<word>("alpha", "alpha");
const word UName = setWavesDict.lookupOrDefault<word>("U", "U");
const bool liquid = setWavesDict.lookupOrDefault<bool>("liquid", true);
const dimensionedVector UMean
(
"UMean",
dimVelocity,
setWavesDict.lookup("UMean")
);
// Get the wave models
const waveSuperposition& waves = waveSuperposition::New(mesh);
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
@ -83,7 +115,7 @@ int main(int argc, char *argv[])
(
IOobject
(
args.optionFound("alpha") ? args["alpha"] : "alpha",
alphaName,
runTime.timeName(),
mesh,
IOobject::MUST_READ
@ -94,7 +126,7 @@ int main(int argc, char *argv[])
(
IOobject
(
args.optionFound("U") ? args["U"] : "U",
UName,
runTime.timeName(),
mesh,
IOobject::MUST_READ
@ -125,7 +157,7 @@ int main(int argc, char *argv[])
(
IOobject("uGasp", runTime.timeName(), mesh),
pMesh,
dimensionedVector("0", dimLength, vector::zero)
dimensionedVector("0", dimVelocity, vector::zero)
);
volVectorField uLiq
(
@ -137,165 +169,31 @@ int main(int argc, char *argv[])
(
IOobject("uLiqp", runTime.timeName(), mesh),
pMesh,
dimensionedVector("0", dimLength, vector::zero)
dimensionedVector("0", dimVelocity, vector::zero)
);
// The number of wave patches
label nWaves = 0;
// Offset
const vector offset = UMean.value()*t;
// Whether the alpha conditions refer to the liquid phase
bool liquid = false;
// Cell centres and points
const pointField& ccs = mesh.cellCentres();
const pointField& pts = mesh.points();
// Loop the patches, averaging and superimposing wave model data
forAll(mesh.boundary(), patchi)
// Internal field
h.primitiveFieldRef() = waves.height(t, ccs + offset);
hp.primitiveFieldRef() = waves.height(t, pts + offset);
uGas.primitiveFieldRef() = waves.UGas(t, ccs + offset);
uGasp.primitiveFieldRef() = waves.UGas(t, pts + offset);
uLiq.primitiveFieldRef() = waves.ULiquid(t, ccs + offset);
uLiqp.primitiveFieldRef() = waves.ULiquid(t, pts + offset);
// Boundary fields
forAll(mesh.boundary(), patchj)
{
fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
const bool isWave = isA<waveAlphaFvPatchScalarField>(alphap);
if (isA<waveVelocityFvPatchVectorField>(Up) != isWave)
{
FatalErrorInFunction
<< "The alpha condition on patch " << Up.patch().name()
<< " is " << alphap.type() << " and the velocity condition"
<< " is " << Up.type() << ". Wave boundary conditions must"
<< " be set in pairs. If the alpha condition is "
<< waveAlphaFvPatchScalarField::typeName
<< " then the velocity condition must be "
<< waveVelocityFvPatchVectorField::typeName
<< " and vice-versa." << exit(FatalError);
}
if (!isWave)
{
continue;
}
Info<< "Adding waves from patch " << Up.patch().name() << endl;
const waveSuperposition& waves =
refCast<waveVelocityFvPatchVectorField>(Up).waves();
const bool liquidp =
refCast<waveAlphaFvPatchScalarField>(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 pointField& ccs = mesh.cellCentres();
const pointField& pts = mesh.points();
// Internal field superposition
h.primitiveFieldRef() += waves.height(t, ccs);
hp.primitiveFieldRef() += waves.height(t, pts);
uGas.primitiveFieldRef() += waves.UGas(t, ccs) - waves.UMean(t);
uGasp.primitiveFieldRef() += waves.UGas(t, pts) - waves.UMean(t);
uLiq.primitiveFieldRef() += waves.ULiquid(t, ccs) - waves.UMean(t);
uLiqp.primitiveFieldRef() += waves.ULiquid(t, pts) - waves.UMean(t);
// Boundary field superposition
forAll(mesh.boundary(), patchj)
{
const pointField& fcs = mesh.boundary()[patchj].Cf();
h.boundaryFieldRef()[patchj] += waves.height(t, fcs);
uGas.boundaryFieldRef()[patchj] +=
waves.UGas(t, fcs) - waves.UMean(t);
uLiq.boundaryFieldRef()[patchj] +=
waves.ULiquid(t, fcs) - waves.UMean(t);
}
++ nWaves;
}
// Create the mean velocity field
volVectorField UMean
(
IOobject("UMean", runTime.timeName(), mesh),
mesh,
dimensionedVector("UMean", dimVelocity, Zero)
);
if (nWaves == 0)
{
// Warn and skip to the next time if there are no wave patches
WarningInFunction
<< "No " << waveAlphaFvPatchScalarField::typeName << " or "
<< waveVelocityFvPatchVectorField::typeName << " patch fields "
<< "were found. No waves have been set." << endl;
continue;
}
else if (nWaves == 1)
{
// Set a mean velocity equal to that on the only wave patch
forAll(mesh.boundary(), patchi)
{
const fvPatchVectorField& Up = U.boundaryField()[patchi];
if (!isA<waveVelocityFvPatchVectorField>(Up))
{
continue;
}
const waveSuperposition& waves =
refCast<const waveVelocityFvPatchVectorField>(Up).waves();
UMean ==
dimensionedVector("UMean", dimVelocity, waves.UMean(t));
}
}
else if (nWaves > 1)
{
// Set the mean velocity by distance weighting from the wave patches
// 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)
{
const fvPatchVectorField& Up = U.boundaryField()[patchi];
if (!isA<waveVelocityFvPatchVectorField>(Up))
{
continue;
}
const waveSuperposition& waves =
refCast<const waveVelocityFvPatchVectorField>(Up).waves();
const volScalarField w
(
1
/(
wallDist(mesh, labelList(1, patchi)).y()
+ dimensionedScalar("ySmall", dimLength, small)
)
);
weight += w;
weightUMean +=
w*dimensionedVector("wUMean", dimVelocity, waves.UMean(t));
}
// Complete the average for the mean velocity
UMean = weightUMean/weight;
const pointField& fcs = mesh.boundary()[patchj].Cf();
h.boundaryFieldRef()[patchj] = waves.height(t, fcs + offset);
uGas.boundaryFieldRef()[patchj] = waves.UGas(t, fcs + offset);
uLiq.boundaryFieldRef()[patchj] = waves.ULiquid(t, fcs + offset);
}
// Set the fields
@ -306,10 +204,13 @@ int main(int argc, char *argv[])
forAll(mesh.boundary(), patchi)
{
fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
if (isA<waveAlphaFvPatchScalarField>(alphap))
{
alphap == refCast<waveAlphaFvPatchScalarField>(alphap).alpha();
}
fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
if (isA<waveVelocityFvPatchVectorField>(Up))
{
Up == refCast<waveVelocityFvPatchVectorField>(Up).U();
}
}