Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
mattijs
2012-01-31 12:34:04 +00:00
66 changed files with 744 additions and 505 deletions

View File

@ -1,3 +0,0 @@
channelFoam.C
EXE = $(FOAM_APPBIN)/channelFoam

View File

@ -1,14 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/LES/LESModel \
-I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lincompressibleLESModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools

View File

@ -1,154 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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
channelFoam
Description
Incompressible LES solver for flow in a channel.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "LESModel.H"
#include "IFstream.H"
#include "OFstream.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readTransportProperties.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "createGradP.H"
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "readPISOControls.H"
#include "CourantNo.H"
sgsModel->correct();
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ sgsModel->divDevBeff(U)
==
flowDirection*gradP
);
if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
volScalarField rAU(1.0/UEqn.A());
for (int corr=0; corr<nCorr; corr++)
{
U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
#include "continuityErrs.H"
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
// Correct driving force for a constant mass flow rate
// Extract the velocity in the flow direction
dimensionedScalar magUbarStar =
(flowDirection & U)().weightedAverage(mesh.V());
// Calculate the pressure gradient increment needed to
// adjust the average flow-rate to the correct value
dimensionedScalar gragPplus =
(magUbar - magUbarStar)/rAU.weightedAverage(mesh.V());
U += flowDirection*rAU*gragPplus;
gradP += gragPplus;
Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
<< "pressure gradient = " << gradP.value() << endl;
runTime.write();
#include "writeGradP.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,43 +0,0 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::LESModel> sgsModel
(
incompressible::LESModel::New(U, phi, laminarTransport)
);

View File

@ -1,24 +0,0 @@
dimensionedScalar gradP
(
"gradP",
dimensionSet(0, 1, -2, 0, 0),
0.0
);
IFstream gradPFile
(
runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
);
if (gradPFile.good())
{
gradPFile >> gradP;
Info<< "Reading average pressure gradient" <<endl
<< endl;
}
else
{
Info<< "Initializing with 0 pressure gradient" <<endl
<< endl;
};

View File

@ -1,29 +0,0 @@
Info<< "\nReading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
);
dimensionedScalar nu
(
transportProperties.lookup("nu")
);
// Read centerline velocity for channel simulations
dimensionedVector Ubar
(
transportProperties.lookup("Ubar")
);
dimensionedScalar magUbar = mag(Ubar);
vector flowDirection = (Ubar/magUbar).value();

View File

@ -1,19 +0,0 @@
if (runTime.outputTime())
{
OFstream gradPFile
(
runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
);
if (gradPFile.good())
{
gradPFile << gradP << endl;
}
else
{
FatalErrorIn(args.executable())
<< "Cannot open file "
<< runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
<< exit(FatalError);
};
};

View File

@ -9,4 +9,6 @@
UrelEqn().relax();
solve(UrelEqn() == -fvc::grad(p));
sources.constrain(UrelEqn());
solve(UrelEqn() == -fvc::grad(p) + sources(Urel));

View File

@ -1,5 +1,5 @@
volScalarField rAUrel(1.0/UrelEqn().A());
Urel = rAUrel*UrelEqn().H();
Urel = rAUrel*(UrelEqn() == sources(Urel))().H();
if (pimple.nCorrPISO() <= 1)
{
@ -37,3 +37,4 @@ p.relax();
// Momentum corrector
Urel -= rAUrel*fvc::grad(p);
Urel.correctBoundaryConditions();
sources.correct(Urel);

View File

@ -5,8 +5,6 @@ tmp<fvVectorMatrix> UEqn
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
sources(U)
);
UEqn().relax();
@ -17,5 +15,5 @@ volScalarField rAU(1.0/UEqn().A());
if (pimple.momentumPredictor())
{
solve(UEqn() == -fvc::grad(p));
solve(UEqn() == -fvc::grad(p) + sources(U));
}

View File

@ -1,4 +1,4 @@
U = rAU*UEqn().H();
U = rAU*(UEqn() == sources(U))().H();
if (pimple.nCorrPISO() <= 1)
{
@ -36,3 +36,4 @@ p.relax();
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);

View File

@ -9,9 +9,11 @@ tmp<fvVectorMatrix> UEqn
UEqn().relax();
sources.constrain(UEqn());
rAU = 1.0/UEqn().A();
if (pimple.momentumPredictor())
{
solve(UEqn() == -fvc::grad(p));
solve(UEqn() == -fvc::grad(p) + sources(U));
}

View File

@ -1,4 +1,4 @@
U = rAU*UEqn().H();
U = rAU*(UEqn() == sources(U))().H();
if (pimple.nCorrPISO() <= 1)
{
@ -46,3 +46,4 @@ fvc::makeRelative(phi, U);
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);

View File

@ -3,8 +3,6 @@ tmp<fvVectorMatrix> UEqn
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
sources(U)
);
@ -14,5 +12,5 @@ sources.constrain(UEqn());
if (pimple.momentumPredictor())
{
solve(UEqn() == -fvc::grad(p_gh));
solve(UEqn() == -fvc::grad(p_gh) + sources(U));
}

View File

@ -1,7 +1,7 @@
volScalarField rAU(1.0/UEqn().A());
surfaceScalarField rAUf(rAU.name() + 'f', fvc::interpolate(rAU));
U = rAU*UEqn().H();
U = rAU*(UEqn() == sources(U))().H();
if (pimple.nCorrPISO() <= 1)
{
@ -41,3 +41,4 @@ p = p_gh + (g & (mesh.C() + zeta - refLevel));
U -= rAU*fvc::grad(p_gh);
U.correctBoundaryConditions();
sources.correct(U);

View File

@ -34,4 +34,5 @@
// Momentum corrector
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
}

View File

@ -34,4 +34,5 @@
// Momentum corrector
Urel -= rAUrel*fvc::grad(p);
Urel.correctBoundaryConditions();
sources.correct(Urel);
}

View File

@ -34,4 +34,5 @@
// Momentum corrector
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
}

View File

@ -7,7 +7,6 @@
==
rho.dimensionedInternalField()*g
+ parcels.SU(U)
+ sources(rho, U)
);
sources.constrain(UEqn);
@ -16,5 +15,5 @@
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
solve(UEqn == -fvc::grad(p) + sources(rho, U));
}

View File

@ -6,7 +6,7 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
U = rAU*UEqn.H();
U = rAU*(UEqn == sources(rho, U))().H();
if (pZones.size() > 0)
{
@ -60,6 +60,7 @@
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
rho = thermo.rho();
rho = max(rho, rhoMin);

View File

@ -7,7 +7,6 @@
rho.dimensionedInternalField()*g
+ coalParcels.SU(U)
+ limestoneParcels.SU(U)
+ sources(rho, U)
);
UEqn.relax();
@ -16,6 +15,6 @@
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
solve(UEqn == -fvc::grad(p) + sources(rho, U));
K = 0.5*magSqr(U);
}

View File

@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
U = rAU*UEqn.H();
U = rAU*(UEqn == sources(rho, U))().H();
if (pimple.transonic())
{
@ -74,6 +74,8 @@ else
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
K = 0.5*magSqr(U);
dpdt = fvc::ddt(p);

View File

@ -7,7 +7,6 @@
==
rho.dimensionedInternalField()*g
+ parcels.SU(U)
+ sources(rho, U)
);
UEqn.relax();
@ -18,7 +17,7 @@
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
solve(UEqn == -fvc::grad(p) + sources(rho, U));
K = 0.5*magSqr(U);
}

View File

@ -6,7 +6,7 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
U = rAU*UEqn.H();
U = rAU*(UEqn == sources(rho, U))().H();
if (pZones.size() > 0)
{
@ -59,6 +59,8 @@
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
K = 0.5*magSqr(U);
dpdt = fvc::ddt(p);

View File

@ -30,7 +30,7 @@
}
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;

View File

@ -155,7 +155,6 @@
phasei++;
}
Dp = mag(Dp);
adjustPhi(phi0, U, p);
while (pimple.correctNonOrthogonal())
{

View File

@ -203,10 +203,15 @@ void Foam::Time::setControls()
)
);
if (timeDict.readIfPresent("deltaT", deltaT_))
// Read and set the deltaT only if time-step adjustment is active
// otherwise use the deltaT from the controlDict
if (controlDict_.lookupOrDefault<Switch>("adjustTimeStep", false))
{
deltaTSave_ = deltaT_;
deltaT0_ = deltaT_;
if (timeDict.readIfPresent("deltaT", deltaT_))
{
deltaTSave_ = deltaT_;
deltaT0_ = deltaT_;
}
}
timeDict.readIfPresent("deltaT0", deltaT0_);

View File

@ -36,7 +36,6 @@ void Foam::interpolation2DTable<Type>::readTable()
// Read data from file
reader_()(fName, *this);
//IFstream(fName)() >> *this;
if (this->empty())
{
@ -361,7 +360,6 @@ void Foam::interpolation2DTable<Type>::write(Ostream& os) const
{
reader_->write(os);
}
//*this >> os;
}
@ -654,100 +652,6 @@ Type Foam::interpolation2DTable<Type>::operator()
}
}
label loY2 = 0;
label hiY2 = 0;
nY = matrix::operator[](hiX).second().size();
minLimit = matrix::operator[](loX).second()[0].first();
maxLimit = matrix::operator[](loX).second()[nY-1].first();
if (lookupValue < minLimit)
{
switch (boundsHandling_)
{
case interpolation2DTable::ERROR:
{
FatalErrorIn
(
"Foam::interpolation2DTable<Type>::operator[]"
"(const scalar, const scalar) const"
) << "value (" << lookupValue << ") underflow" << nl
<< exit(FatalError);
break;
}
case interpolation2DTable::WARN:
{
WarningIn
(
"Foam::interpolation2DTable<Type>::operator[]"
"(const scalar, const scalar) const"
) << "value (" << lookupValue << ") underflow" << nl
<< " Continuing with the first entry"
<< endl;
// fall-through to 'CLAMP'
}
case interpolation2DTable::CLAMP:
{
loY2 = 0;
loY2 = 1;
break;
}
}
}
else if (lookupValue >= maxLimit)
{
switch (boundsHandling_)
{
case interpolation2DTable::ERROR:
{
FatalErrorIn
(
"Foam::interpolation2DTable<Type>::operator[]"
"(const scalar, const scalar) const"
) << "value (" << lookupValue << ") overflow" << nl
<< exit(FatalError);
break;
}
case interpolation2DTable::WARN:
{
WarningIn
(
"Foam::interpolation2DTable<Type>::operator[]"
"(const scalar, const scalar) const"
) << "value (" << lookupValue << ") overflow" << nl
<< " Continuing with the last entry"
<< endl;
// fall-through to 'CLAMP'
}
case interpolation2DTable::CLAMP:
{
loY2 = nY-1;
loY2 = nY;
break;
}
}
}
else
{
// Finds the lo and hi of Y on the high x
for (label i = 0; i < nY; ++i)
{
if
(
lookupValue >= matrix::operator[](hiX).second()[i].first()
)
{
loY2 = hiY2 = i;
}
else
{
hiY2 = i;
break;
}
}
}
if (loX == hiX)
{
// we are at the end of the table - or there is only a single entry

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -45,7 +45,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class IObasicSourceList Declaration
Class IObasicSourceList Declaration
\*---------------------------------------------------------------------------*/
class IObasicSourceList
@ -80,6 +80,8 @@ public:
{}
// Member Functions
//- Read dictionary
virtual bool read();
};

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,6 +26,7 @@ License
#include "basicSource.H"
#include "fvMesh.H"
#include "fvMatrices.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -161,7 +162,7 @@ void Foam::basicSource::setCellSet()
}
case smMapRegion:
{
if(active_)
if (active_)
{
Info<< indent << "- selecting inter region mapping" << endl;
const fvMesh& secondaryMesh =
@ -170,7 +171,6 @@ void Foam::basicSource::setCellSet()
const boundBox secondaryBB = secondaryMesh.bounds();
if (secondaryBB.overlaps(primaryBB))
{
// Dummy patches
wordList cuttingPatches;
HashTable<word> patchMap;
@ -218,7 +218,7 @@ void Foam::basicSource::setCellSet()
}
// Set volume information
if(selectionMode_ != smMapRegion)
if (selectionMode_ != smMapRegion)
{
V_ = 0.0;
forAll(cells_, i)
@ -303,6 +303,7 @@ Foam::autoPtr<Foam::basicSource> Foam::basicSource::New
return autoPtr<basicSource>(cstrIter()(name, modelType, coeffs, mesh));
}
Foam::basicSource::~basicSource()
{
if (!secondaryToPrimaryInterpPtr_.empty())
@ -364,6 +365,36 @@ void Foam::basicSource::checkApplied() const
}
void Foam::basicSource::correct(volScalarField& fld)
{
// do nothing
}
void Foam::basicSource::correct(volVectorField& fld)
{
// do nothing
}
void Foam::basicSource::correct(volSphericalTensorField& fld)
{
// do nothing
}
void Foam::basicSource::correct(volSymmTensorField& fld)
{
// do nothing
}
void Foam::basicSource::correct(volTensorField& fld)
{
// do nothing
}
void Foam::basicSource::addSup(fvMatrix<scalar>& eqn, const label fieldI)
{
// do nothing

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -47,6 +47,7 @@ SourceFiles
#define basicSource_H
#include "fvMatricesFwd.H"
#include "volFieldsFwd.H"
#include "cellSet.H"
#include "autoPtr.H"
#include "meshToMesh.H"
@ -60,7 +61,6 @@ namespace Foam
class fvMesh;
/*---------------------------------------------------------------------------*\
Class basicSource Declaration
\*---------------------------------------------------------------------------*/
@ -330,6 +330,24 @@ public:
// Evaluation
// Correct
//- Scalar
virtual void correct(volScalarField& fld);
//- Vector
virtual void correct(volVectorField& fld);
//- Spherical tensor
virtual void correct(volSphericalTensorField& fld);
//- Symmetric tensor
virtual void correct(volSymmTensorField& fld);
//- Tensor
virtual void correct(volTensorField& fld);
// Add explicit and implicit contributions
//- Scalar

View File

@ -135,4 +135,5 @@ inline const Foam::word Foam::basicSource::mapRegionName() const
return mapRegionName_;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -95,6 +95,11 @@ public:
// Member Functions
//- Correct
template<class Type>
void correct(GeometricField<Type, fvPatchField, volMesh>& fld);
// Sources
//- Return source for equation

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,6 +25,39 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::basicSourceList::correct
(
GeometricField<Type, fvPatchField, volMesh>& fld
)
{
const word& fieldName = fld.name();
forAll(*this, i)
{
basicSource& source = this->operator[](i);
label fieldI = source.applyToField(fieldName);
if (fieldI != -1)
{
source.setApplied(fieldI);
if (source.isActive())
{
if (debug)
{
Info<< "Correcting source " << source.name()
<< " for field " << fieldName << endl;
}
source.correct(fld);
}
}
}
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type> > Foam::basicSourceList::operator()
(

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -69,53 +69,6 @@ void Foam::pressureGradientExplicitSource::writeGradP() const
}
void Foam::pressureGradientExplicitSource::update(fvMatrix<vector>& eqn)
{
volVectorField& U = const_cast<volVectorField&>(eqn.psi());
const volScalarField& rAU =
mesh_.lookupObject<volScalarField>("(1|A(" + U.name() + "))");
// Integrate flow variables over cell set
scalar magUbarAve = 0.0;
scalar rAUave = 0.0;
const scalarField& cv = mesh_.V();
forAll(cells_, i)
{
label cellI = cells_[i];
scalar volCell = cv[cellI];
magUbarAve += (flowDir_ & U[cellI])*volCell;
rAUave += rAU[cellI]*volCell;
}
// Collect across all processors
reduce(magUbarAve, sumOp<scalar>());
// Volume averages
magUbarAve /= V_;
rAUave /= V_;
// Calculate the pressure gradient increment needed to adjust the average
// flow-rate to the desired value
scalar gradPplus = (mag(Ubar_) - magUbarAve)/rAUave;
// Apply correction to velocity field
forAll(cells_, i)
{
label cellI = cells_[i];
U[cellI] += flowDir_*rAU[cellI]*gradPplus;
}
// Update pressure gradient
gradP_.value() += gradPplus;
Info<< "Uncorrected Ubar = " << magUbarAve << tab
<< "Pressure gradient = " << gradP_.value() << endl;
writeGradP();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
@ -123,7 +76,7 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
const word& sourceName,
const word& modelType,
const dictionary& dict,
const fvMesh& mesh
const fvMesh& mesh
)
:
basicSource(sourceName, modelType, dict, mesh),
@ -133,6 +86,23 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
flowDir_(Ubar_/mag(Ubar_))
{
coeffs_.lookup("fieldNames") >> fieldNames_;
if (fieldNames_.size() != 1)
{
FatalErrorIn
(
"Foam::pressureGradientExplicitSource::"
"pressureGradientExplicitSource"
"("
"onst word&, "
"const word&, "
"const dictionary&, "
"const fvMesh&"
")"
) << "Source can only be applied to a single field. Current "
<< "settings are:" << fieldNames_ << exit(FatalError);
}
applied_.setSize(fieldNames_.size(), false);
// Read the initial pressure gradient from file if it exists
@ -160,8 +130,6 @@ void Foam::pressureGradientExplicitSource::addSup
const label fieldI
)
{
update(eqn);
DimensionedField<vector, volMesh> Su
(
IOobject
@ -178,7 +146,53 @@ void Foam::pressureGradientExplicitSource::addSup
UIndirectList<vector>(Su, cells_) = flowDir_*gradP_.value();
eqn -= Su;
eqn += Su;
}
void Foam::pressureGradientExplicitSource::correct(volVectorField& U)
{
const volScalarField& rAU =
mesh_.lookupObject<volScalarField>("(1|A(" + U.name() + "))");
// Integrate flow variables over cell set
scalar magUbarAve = 0.0;
scalar rAUave = 0.0;
const scalarField& cv = mesh_.V();
forAll(cells_, i)
{
label cellI = cells_[i];
scalar volCell = cv[cellI];
magUbarAve += (flowDir_ & U[cellI])*volCell;
rAUave += rAU[cellI]*volCell;
}
// Collect across all processors
reduce(magUbarAve, sumOp<scalar>());
reduce(rAUave, sumOp<scalar>());
// Volume averages
magUbarAve /= V_;
rAUave /= V_;
// Calculate the pressure gradient increment needed to adjust the average
// flow-rate to the desired value
scalar gradPplus = (mag(Ubar_) - magUbarAve)/rAUave;
// Apply correction to velocity field
forAll(cells_, i)
{
label cellI = cells_[i];
U[cellI] += flowDir_*rAU[cellI]*gradPplus;
}
// Update pressure gradient
gradP_.value() += gradPplus;
Info<< "Pressure gradient source: uncorrected Ubar = " << magUbarAve
<< ", pressure gradient = " << gradP_.value() << endl;
writeGradP();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -118,7 +118,10 @@ public:
// Member Functions
// Access
// Evaluate
//- Correct the pressure gradient
virtual void correct(volVectorField& U);
//- Add explicit contribution to equation
virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);

View File

@ -170,7 +170,7 @@ $(derivedFvPatchFields)/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvP
$(derivedFvPatchFields)/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/outletMappedUniformInlet/outletMappedUniformInletFvPatchFields.C
$(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
$(derivedFvPatchFields)/phaseHydrostaticPressure/phaseHydrostaticPressureFvPatchScalarField.C
fvsPatchFields = fields/fvsPatchFields
$(fvsPatchFields)/fvsPatchField/fvsPatchFields.C

View File

@ -0,0 +1,202 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
\*---------------------------------------------------------------------------*/
#include "phaseHydrostaticPressureFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "uniformDimensionedFields.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::phaseHydrostaticPressureFvPatchScalarField::
phaseHydrostaticPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(p, iF),
phaseName_("alpha"),
rho_(0.0),
pRefValue_(0.0),
pRefPoint_(vector::zero)
{
this->refValue() = 0.0;
this->refGrad() = 0.0;
this->valueFraction() = 0.0;
}
Foam::phaseHydrostaticPressureFvPatchScalarField::
phaseHydrostaticPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
mixedFvPatchScalarField(p, iF),
phaseName_(dict.lookupOrDefault<word>("phaseName", "alpha")),
rho_(readScalar(dict.lookup("rho"))),
pRefValue_(readScalar(dict.lookup("pRefValue"))),
pRefPoint_(dict.lookup("pRefPoint"))
{
this->refValue() = pRefValue_;
if (dict.found("value"))
{
fvPatchScalarField::operator=
(
scalarField("value", dict, p.size())
);
}
else
{
fvPatchScalarField::operator=(this->refValue());
}
this->refGrad() = 0.0;
this->valueFraction() = 0.0;
}
Foam::phaseHydrostaticPressureFvPatchScalarField::
phaseHydrostaticPressureFvPatchScalarField
(
const phaseHydrostaticPressureFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
mixedFvPatchScalarField(ptf, p, iF, mapper),
phaseName_(ptf.phaseName_),
rho_(ptf.rho_),
pRefValue_(ptf.pRefValue_),
pRefPoint_(ptf.pRefPoint_)
{}
Foam::phaseHydrostaticPressureFvPatchScalarField::
phaseHydrostaticPressureFvPatchScalarField
(
const phaseHydrostaticPressureFvPatchScalarField& ptf
)
:
mixedFvPatchScalarField(ptf),
phaseName_(ptf.phaseName_)
{}
Foam::phaseHydrostaticPressureFvPatchScalarField::
phaseHydrostaticPressureFvPatchScalarField
(
const phaseHydrostaticPressureFvPatchScalarField& ptf,
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchScalarField(ptf, iF),
phaseName_(ptf.phaseName_),
rho_(ptf.rho_),
pRefValue_(ptf.pRefValue_),
pRefPoint_(ptf.pRefPoint_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::phaseHydrostaticPressureFvPatchScalarField::updateCoeffs()
{
if (this->updated())
{
return;
}
const scalarField& alphap =
patch().lookupPatchField<volScalarField, scalar>
(
phaseName_
);
const uniformDimensionedVectorField& g =
db().lookupObject<uniformDimensionedVectorField>("g");
// scalar rhor = 1000;
// scalarField alphap1 = max(min(alphap, 1.0), 0.0);
// valueFraction() = alphap1/(alphap1 + rhor*(1.0 - alphap1));
valueFraction() = max(min(alphap, 1.0), 0.0);
refValue() =
pRefValue_
+ rho_*((g.value() & patch().Cf()) - (g.value() & pRefPoint_));
mixedFvPatchScalarField::updateCoeffs();
}
void Foam::phaseHydrostaticPressureFvPatchScalarField::write(Ostream& os) const
{
fvPatchScalarField::write(os);
if (phaseName_ != "alpha")
{
os.writeKeyword("phaseName")
<< phaseName_ << token::END_STATEMENT << nl;
}
os.writeKeyword("rho") << rho_ << token::END_STATEMENT << nl;
os.writeKeyword("pRefValue") << pRefValue_ << token::END_STATEMENT << nl;
os.writeKeyword("pRefPoint") << pRefPoint_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::phaseHydrostaticPressureFvPatchScalarField::operator=
(
const fvPatchScalarField& ptf
)
{
fvPatchScalarField::operator=
(
valueFraction()*refValue()
+ (1 - valueFraction())*ptf
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchScalarField,
phaseHydrostaticPressureFvPatchScalarField
);
}
// ************************************************************************* //

View File

@ -0,0 +1,223 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
Class
Foam::phaseHydrostaticPressureFvPatchScalarField
Description
Phase hydrostatic pressure boundary condition calculated as
pRefValue + rho*g.(x - pRefPoint)
where rho is provided and assumed uniform
applied according to the phase-fraction field provided:
1 -> fix value to that provided
0 -> zero-gradient
SourceFiles
phaseHydrostaticPressureFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef phaseHydrostaticPressureFvPatchScalarField_H
#define phaseHydrostaticPressureFvPatchScalarField_H
#include "mixedFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class phaseHydrostaticPressureFvPatch Declaration
\*---------------------------------------------------------------------------*/
class phaseHydrostaticPressureFvPatchScalarField
:
public mixedFvPatchScalarField
{
protected:
// Protected data
//- Name of phase-fraction field
word phaseName_;
//- Constant density in the far-field
scalar rho_;
//- Reference pressure
scalar pRefValue_;
//- Reference pressure location
vector pRefPoint_;
public:
//- Runtime type information
TypeName("phaseHydrostaticPressure");
// Constructors
//- Construct from patch and internal field
phaseHydrostaticPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
phaseHydrostaticPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given
// phaseHydrostaticPressureFvPatchScalarField onto a new patch
phaseHydrostaticPressureFvPatchScalarField
(
const phaseHydrostaticPressureFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
phaseHydrostaticPressureFvPatchScalarField
(
const phaseHydrostaticPressureFvPatchScalarField&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField >
(
new phaseHydrostaticPressureFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
phaseHydrostaticPressureFvPatchScalarField
(
const phaseHydrostaticPressureFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchScalarField> clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchScalarField>
(
new phaseHydrostaticPressureFvPatchScalarField(*this, iF)
);
}
// Member functions
// Access
//- Return the phaseName
const word& phaseName() const
{
return phaseName_;
}
//- Return reference to the phaseName to allow adjustment
word& phaseName()
{
return phaseName_;
}
//- Return the constant density in the far-field
scalar rho() const
{
return rho_;
}
//- Return reference to the constant density in the far-field
// to allow adjustment
scalar& rho()
{
return rho_;
}
//- Return the reference pressure
scalar pRefValue() const
{
return pRefValue_;
}
//- Return reference to the reference pressure to allow adjustment
scalar& pRefValue()
{
return pRefValue_;
}
//- Return the pressure reference location
const vector& pRefPoint() const
{
return pRefPoint_;
}
//- Return reference to the pressure reference location
// to allow adjustment
vector& pRefPoint()
{
return pRefPoint_;
}
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
// Member operators
virtual void operator=(const fvPatchScalarField& pvf);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -26,6 +26,7 @@ License
#include "PatchPostProcessing.H"
#include "Pstream.H"
#include "stringListOps.H"
#include "ListOps.H"
#include "ListListOps.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -55,9 +56,12 @@ void Foam::PatchPostProcessing<CloudType>::write()
{
forAll(patchData_, i)
{
List<List<scalar> > procTimes(Pstream::nProcs());
procTimes[Pstream::myProcNo()] = times_[i];
Pstream::gatherList(procTimes);
List<List<string> > procData(Pstream::nProcs());
procData[Pstream::myProcNo()] = patchData_[i];
Pstream::gatherList(procData);
if (Pstream::master())
@ -100,18 +104,33 @@ void Foam::PatchPostProcessing<CloudType>::write()
procData,
accessOp<List<string> >()
);
sort(globalData);
List<scalar> globalTimes;
globalTimes = ListListOps::combine<List<scalar> >
(
procTimes,
accessOp<List<scalar> >()
);
labelList indices;
sortedOrder(globalTimes, indices);
string header("# Time currentProc " + parcelType::propHeader);
patchOutFile<< header.c_str() << nl;
forAll(globalData, dataI)
forAll(globalTimes, i)
{
patchOutFile<< globalData[dataI].c_str() << nl;
label dataI = indices[i];
patchOutFile
<< globalTimes[dataI] << ' '
<< globalData[dataI].c_str()
<< nl;
}
}
patchData_[i].clearStorage();
times_[i].clearStorage();
}
}
@ -128,6 +147,7 @@ Foam::PatchPostProcessing<CloudType>::PatchPostProcessing
CloudFunctionObject<CloudType>(dict, owner, typeName),
maxStoredParcels_(readScalar(this->coeffDict().lookup("maxStoredParcels"))),
patchIDs_(),
times_(),
patchData_()
{
const wordList allPatchNames = owner.mesh().boundaryMesh().names();
@ -167,6 +187,7 @@ Foam::PatchPostProcessing<CloudType>::PatchPostProcessing
}
patchData_.setSize(patchIDs_.size());
times_.setSize(patchIDs_.size());
}
@ -179,6 +200,7 @@ Foam::PatchPostProcessing<CloudType>::PatchPostProcessing
CloudFunctionObject<CloudType>(ppm),
maxStoredParcels_(ppm.maxStoredParcels_),
patchIDs_(ppm.patchIDs_),
times_(ppm.times_),
patchData_(ppm.patchData_)
{}
@ -203,9 +225,11 @@ void Foam::PatchPostProcessing<CloudType>::postPatch
const label localPatchI = applyToPatch(patchI);
if (localPatchI != -1 && patchData_[localPatchI].size() < maxStoredParcels_)
{
times_[localPatchI].append(this->owner().time().value());
OStringStream data;
data<< this->owner().time().timeName() << ' ' << Pstream::myProcNo()
<< ' ' << p;
data<< Pstream::myProcNo() << ' ' << p;
patchData_[localPatchI].append(data.str());
}
}

View File

@ -61,6 +61,9 @@ class PatchPostProcessing
//- List of patch indices to post-process
labelList patchIDs_;
//- List of time for each data record
List<DynamicList<scalar> > times_;
//- List of output data per patch
List<DynamicList<string> > patchData_;

View File

@ -39,6 +39,7 @@ usage()
usage: ${0##*/} [OPTION]
options:
-root <dir> specify root folder to run tests from
-default sets up a default scheme on all schemes
-help print the usage
@ -52,16 +53,22 @@ USAGE
unset DEFAULT_SCHEMES
ROOT="./"
# parse options
while [ "$#" -gt 0 ]
do
case "$1" in
-r | -root)
[ "$#" -ge 2 ] || usage "'$1' option requires an argument"
ROOT="$2"
shift
;;
-h | -help)
usage
;;
-d | -default)
DEFAULT_SCHEMES=true
shift
;;
-*)
usage "unknown option: '$*'"
@ -70,6 +77,7 @@ do
break
;;
esac
shift
done
@ -123,7 +131,7 @@ done
[ -f "$MAIN_CONTROL_DICT" ] || usage "main controlDict not found"
TUTORIALS_DIR=.
TUTORIALS_DIR=$ROOT
TEST_RUN_DIR=../tutorialsTest
FV_SCHEMES=\
" \
@ -198,8 +206,10 @@ then
done
fi
cp -f $FOAM_TUTORIALS/Allrun .
./Allrun
sed -e :a -e '/\\$/N; s/\\\n//; ta' Allrun > temp
APPLICATIONS=\
`grep "applications=" temp | sed 's/applications=\"\([A-Za-z \t]*\)\"/\1/g'`

View File

@ -0,0 +1,35 @@
/*--------------------------------*- 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 dictionary;
location "constant";
object sourcesProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
momentumSource
{
type pressureGradientExplicitSource;
active on; //on/off switch
timeStart 0.0; //start time
duration 1000000.0; //duration
selectionMode all; //cellSet // points //cellZone
pressureGradientExplicitSourceCoeffs
{
fieldNames (U);
Ubar ( 0.1335 0 0 );
gradPini gradPIni [0 1 -2 0 0 ] 0;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,20 @@
/*--------------------------------*- 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 dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType LESModel;
// ************************************************************************* //

View File

@ -15,7 +15,7 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application channelFoam;
application pimpleFoam;
startFrom startTime;

View File

@ -33,41 +33,25 @@ solvers
relTol 0;
}
U
"(U|k)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
relTol 0.1;
}
k
"(U|k)Final"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
B
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
nuTilda
{
solver PBiCG;
preconditioner DILU;
$U;
tolerance 1e-05;
relTol 0;
}
}
PISO
PIMPLE
{
nOuterCorrectors 1;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 1001;

View File

@ -128,7 +128,7 @@ relaxationFactors
}
equations
{
"(U|k|epsilon|omega|R|nuTilda)" 1;
"(U|k|epsilon|omega|R|nuTilda).*" 1;
}
}