Commit 17b2e50e authored by henry's avatar henry
Browse files

Added shallowWaterFoam solver and tutorial case.

parent ef537003
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Global
CourantNo
Description
Calculates and outputs the maximum Courant Number.
\*---------------------------------------------------------------------------*/
scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
scalar waveCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField SfUfbyDelta =
mesh.surfaceInterpolation::deltaCoeffs()
*mag(phi)/fvc::interpolate(h);
CoNum = max(SfUfbyDelta/mesh.magSf())
.value()*runTime.deltaT().value();
meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
.value()*runTime.deltaT().value();
// Gravity wave Courant number
waveCoNum =
0.5*max
(
mesh.surfaceInterpolation::deltaCoeffs()
*sqrt(fvc::interpolate(h))
).value()*sqrt(magg).value()*runTime.deltaT().value();
}
Info<< "Courant number mean: " << meanCoNum
<< " max: " << CoNum << endl;
Info<< "Gravity wave Courant number max: " << waveCoNum
<< endl;
// ************************************************************************* //
shallowWaterFoam.C
EXE = $(FOAM_APPBIN)/shallowWaterFoam
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume
Info<< "Reading field h\n" << endl;
volScalarField h
(
IOobject
(
"h",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field h0 if present\n" << endl;
volScalarField h0
(
IOobject
(
"h0",
runTime.findInstance("polyMesh", "points"),
mesh,
IOobject::READ_IF_PRESENT
),
mesh,
dimensionedScalar("h0", dimLength, 0.0)
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Creating field hU\n" << endl;
volVectorField hU
(
IOobject
(
"hU",
runTime.timeName(),
mesh
),
h*U,
U.boundaryField().types()
);
Info<< "Creating field hTotal for post processing\n" << endl;
volScalarField hTotal
(
IOobject
(
"hTotal",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
h+h0
);
hTotal.write();
# include "createPhi.H"
Info<< "Creating Coriolis Force" << endl;
const dimensionedVector F("F", ((2.0*Omega) & gHat)*gHat);
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Global
createPhi
Description
Creates and initialises the face-flux field phi.
\*---------------------------------------------------------------------------*/
#ifndef createPhi_H
#define createPhi_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(hU) & mesh.Sf()
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
Info<< "\nReading environmentalProperties" << endl;
IOdictionary environmentalProperties
(
IOobject
(
"environmentalProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
const dimensionedVector g(environmentalProperties.lookup("g"));
const Switch rotating(environmentalProperties.lookup("rotating"));
const dimensionedVector Omega =
rotating ? environmentalProperties.lookup("Omega")
: dimensionedVector("Omega", -dimTime, vector(0,0,0));
const dimensionedScalar magg = mag(g);
const dimensionedVector gHat = g/magg;
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
shallowWaterFoam
Description
Transient solver for inviscid shallow-water equations with rotation.
If the geometry is 3D then it is assumed to be one layers of cells and the
component of the velocity normal to gravity is removed.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readEnvironmentalProperties.H"
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
runTime++;
Info<< "\n Time = " << runTime.timeName() << nl << endl;
#include "readPISOControls.H"
#include "CourantNo.H"
for (int ucorr=0; ucorr<nOuterCorr; ucorr++)
{
surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
fvVectorMatrix hUEqn
(
fvm::ddt(hU)
+ fvm::div(phiv, hU)
);
hUEqn.relax();
if (momentumPredictor)
{
if (rotating)
{
solve(hUEqn + (F ^ hU) == -magg*h*fvc::grad(h + h0));
}
else
{
solve(hUEqn == -magg*h*fvc::grad(h + h0));
}
// Constrain the momentum to be in the geometry if 3D geometry
if (mesh.nGeometricD() == 3)
{
hU -= (gHat & hU)*gHat;
hU.correctBoundaryConditions();
}
}
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
surfaceScalarField hf = fvc::interpolate(h);
volScalarField rUA = 1.0/hUEqn.A();
surfaceScalarField ghrUAf = magg*fvc::interpolate(h*rUA);
surfaceScalarField phih0 = ghrUAf*mesh.magSf()*fvc::snGrad(h0);
if (rotating)
{
hU = rUA*(hUEqn.H() - (F ^ hU));
}
else
{
hU = rUA*hUEqn.H();
}
phi = (fvc::interpolate(hU) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, h, hU, phi)
- phih0;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix hEqn
(
fvm::ddt(h)
+ fvc::div(phi)
- fvm::laplacian(ghrUAf, h)
);
if (ucorr < nOuterCorr-1 || corr < nCorr-1)
{
hEqn.solve();
}
else
{
hEqn.solve(mesh.solver(h.name() + "Final"));
}
if (nonOrth == nNonOrthCorr)
{
phi += hEqn.flux();
}
}
hU -= rUA*h*magg*fvc::grad(h + h0);
// Constrain the momentum to be in the geometry if 3D geometry
if (mesh.nGeometricD() == 3)
{
hU -= (gHat & hU)*gHat;
}
hU.correctBoundaryConditions();
}
}
U == hU/h;
hTotal == h + h0;
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //
# GMT common arguments shelf
-B0/0
-JX18c/18c
-R0/1/0/1
-jX18c/18c
EOF
F
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.1 0 0);
boundaryField
{
sides
{
type slip;
}
inlet
{
type fixedValue;
value uniform (0.1 0 0);
}
outlet
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object h;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 0 0 0 0 0];
internalField nonuniform List<scalar>
400
(
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01