330 lines
11 KiB
C
330 lines
11 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / 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/>.
|
|
|
|
Application
|
|
setWaves
|
|
|
|
Description
|
|
Applies wave models to the entire domain for case initialisation using
|
|
level sets for second-order accuracy.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "fvCFD.H"
|
|
#include "levelSet.H"
|
|
#include "pointFields.H"
|
|
#include "timeSelector.H"
|
|
#include "wallDist.H"
|
|
#include "waveAlphaFvPatchScalarField.H"
|
|
#include "waveVelocityFvPatchVectorField.H"
|
|
#include "waveSuperposition.H"
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
timeSelector::addOptions(false, false);
|
|
|
|
Foam::argList::addOption
|
|
(
|
|
"U",
|
|
"name",
|
|
"name of the velocity field, default is \"U\""
|
|
);
|
|
|
|
Foam::argList::addOption
|
|
(
|
|
"alpha",
|
|
"name",
|
|
"name of the volume fraction field, default is \"alpha\""
|
|
);
|
|
|
|
#include "setRootCase.H"
|
|
#include "createTime.H"
|
|
|
|
instantList timeDirs = timeSelector::selectIfPresent(runTime, args);
|
|
|
|
#include "createMesh.H"
|
|
#include "readGravitationalAcceleration.H"
|
|
|
|
const pointMesh& pMesh = pointMesh::New(mesh);
|
|
|
|
forAll(timeDirs, timeI)
|
|
{
|
|
runTime.setTime(timeDirs[timeI], timeI);
|
|
|
|
Info<< "Time = " << runTime.timeName() << nl << endl;
|
|
|
|
mesh.readUpdate();
|
|
|
|
// Read the fields which are to be set
|
|
volScalarField alpha
|
|
(
|
|
IOobject
|
|
(
|
|
args.optionFound("alpha") ? args["alpha"] : "alpha",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::MUST_READ
|
|
),
|
|
mesh
|
|
);
|
|
volVectorField U
|
|
(
|
|
IOobject
|
|
(
|
|
args.optionFound("U") ? args["U"] : "U",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::MUST_READ
|
|
),
|
|
mesh
|
|
);
|
|
|
|
// Create modelled fields on both cells and points
|
|
volScalarField h
|
|
(
|
|
IOobject("h", runTime.timeName(), mesh),
|
|
mesh,
|
|
dimensionedScalar("0", dimLength, 0)
|
|
);
|
|
pointScalarField hp
|
|
(
|
|
IOobject("hp", 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 uLiq
|
|
(
|
|
IOobject("uLiq", runTime.timeName(), mesh),
|
|
mesh,
|
|
dimensionedVector("0", dimVelocity, vector::zero)
|
|
);
|
|
pointVectorField uLiqp
|
|
(
|
|
IOobject("uLiqp", runTime.timeName(), mesh),
|
|
pMesh,
|
|
dimensionedVector("0", dimLength, vector::zero)
|
|
);
|
|
|
|
// The number of wave patches
|
|
label nWaves = 0;
|
|
|
|
// Whether the alpha conditions refer to the liquid phase
|
|
bool liquid = false;
|
|
|
|
// Loop the patches, averaging and superimposing wave model data
|
|
forAll(mesh.boundary(), patchi)
|
|
{
|
|
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 scalar t = runTime.value();
|
|
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();
|
|
uGasp.primitiveFieldRef() += waves.UGas(t, pts) - waves.UMean();
|
|
uLiq.primitiveFieldRef() += waves.ULiquid(t, ccs) - waves.UMean();
|
|
uLiqp.primitiveFieldRef() += waves.ULiquid(t, pts) - waves.UMean();
|
|
|
|
// 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();
|
|
uLiq.boundaryFieldRef()[patchj] +=
|
|
waves.ULiquid(t, fcs) - waves.UMean();
|
|
}
|
|
|
|
++ 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());
|
|
}
|
|
}
|
|
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());
|
|
}
|
|
|
|
// Complete the average for the mean velocity
|
|
UMean = weightUMean/weight;
|
|
}
|
|
|
|
// Set the fields
|
|
alpha == levelSetFraction(h, hp, !liquid);
|
|
U == UMean + levelSetAverage(mesh, h, hp, uGas, uGasp, uLiq, uLiqp);
|
|
|
|
// Set the boundary fields
|
|
forAll(mesh.boundary(), patchi)
|
|
{
|
|
fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
|
|
fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
|
|
if (isA<waveAlphaFvPatchScalarField>(alphap))
|
|
{
|
|
alphap == refCast<waveAlphaFvPatchScalarField>(alphap).alpha();
|
|
Up == refCast<waveVelocityFvPatchVectorField>(Up).U();
|
|
}
|
|
}
|
|
|
|
// Output
|
|
Info<< "Writing " << alpha.name() << nl;
|
|
alpha.write();
|
|
Info<< "Writing " << U.name() << nl << endl;
|
|
U.write();
|
|
}
|
|
|
|
Info<< "End\n" << endl;
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|