Formatting

Conflicts:
	applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C
This commit is contained in:
Hrvoje Jasak 2015-01-27 12:12:53 +00:00
parent 2154e0dee1
commit 15a24604c7
15 changed files with 233 additions and 344 deletions

View file

@ -0,0 +1,203 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
potentialDyMFoam
Description
Transient solver for potential flow with dynamic mesh.
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::validOptions.insert("resetU", "");
argList::validOptions.insert("writep", "");
# include "setRootCase.H"
# include "createTime.H"
# include "createDynamicFvMesh.H"
# include "createFields.H"
# include "initTotalVolume.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
# include "readPISOControls.H"
# include "checkTotalVolume.H"
Info<< "Time = " << runTime.timeName() << nl << endl;
bool meshChanged = mesh.update();
reduce(meshChanged, orOp<bool>());
p.internalField() = 0;
if (args.optionFound("resetU"))
{
U.internalField() = vector::zero;
}
# include "volContinuity.H"
# include "meshCourantNo.H"
// Solve potential flow equations
adjustPhi(phi, U, p);
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
p.storePrevIter();
fvScalarMatrix pEqn
(
fvm::laplacian
(
dimensionedScalar
(
"1",
dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
1
),
p
)
==
fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
else
{
p.relax();
}
}
Info<< "continuity error = "
<< mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
<< endl;
U = fvc::reconstruct(phi);
U.correctBoundaryConditions();
Info<< "Interpolated U error = "
<< (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
/sum(mesh.magSf())).value()
<< endl;
// Calculate velocity magnitude
{
volScalarField magU = mag(U);
Info<< "mag(U): max: " << gMax(magU.internalField())
<< " min: " << gMin(magU.internalField()) << endl;
}
if (args.optionFound("writep"))
{
// Find reference patch
label refPatch = -1;
scalar maxMagU = 0;
// Go through all velocity patches and find the one that fixes
// velocity to the largest value
forAll (U.boundaryField(), patchI)
{
const fvPatchVectorField& Upatch = U.boundaryField()[patchI];
if (Upatch.fixesValue())
{
// Calculate mean velocity
scalar u = sum(mag(Upatch));
label patchSize = Upatch.size();
reduce(u, sumOp<scalar>());
reduce(patchSize, sumOp<label>());
if (patchSize > 0)
{
scalar curMag = u/patchSize;
if (curMag > maxMagU)
{
refPatch = patchI;
maxMagU = curMag;
}
}
}
}
if (refPatch > -1)
{
// Calculate reference pressure
const fvPatchVectorField& Upatch = U.boundaryField()[refPatch];
const fvPatchScalarField& pPatch = p.boundaryField()[refPatch];
scalar patchE = sum(mag(pPatch + 0.5*magSqr(Upatch)));
label patchSize = Upatch.size();
reduce(patchE, sumOp<scalar>());
reduce(patchSize, sumOp<label>());
scalar e = patchE/patchSize;
Info<< "Using reference patch " << refPatch
<< " with mag(U) = " << maxMagU
<< " p + 0.5*U^2 = " << e << endl;
p.internalField() = e - 0.5*magSqr(U.internalField());
p.correctBoundaryConditions();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

View file

@ -64,7 +64,7 @@ int main(int argc, char *argv[])
# include "setDeltaT.H"
runTime++;
Info<< "deltaT = " << runTime.deltaT().value() << nl << endl;
Info<< "Time = " << runTime.timeName() << nl << endl;
bool meshChanged = mesh.update();
@ -93,7 +93,7 @@ int main(int argc, char *argv[])
if (meshChanged)
{
# include "CourantNo.H"
# include "compressibleCourantNo.H"
}
// --- PIMPLE loop

View file

@ -1,8 +1,8 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
);
solve(UEqn == -fvc::grad(p));

View file

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

View file

@ -1,10 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-ldynamicMesh \
-lmeshTools \
-llduSolvers

View file

@ -1,20 +0,0 @@
{
# include "rhoEqn.H"
}
{
scalar sumLocalContErr =
sum
(
mag(rho.internalField() - (psi*p)().internalField())
)/sum(rho.internalField());
scalar globalContErr =
sum(rho.internalField() - (psi*p)().internalField())
/sum(rho.internalField());
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr << endl;
}

View file

@ -1,103 +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 T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Calculating field e from T\n" << endl;
volScalarField e
(
IOobject
(
"e",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
Cv*T,
T.boundaryField().types()
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField psi
(
IOobject
(
"psi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
1.0/(R*T)
);
psi.oldTime();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh
),
psi*p
);
# include "compressibleCreatePhi.H"
Info<< "Creating field phid\n" << endl;
surfaceScalarField phid
(
IOobject
(
"phid",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
phi/fvc::interpolate(p),
phi.boundaryField().types()
);

View file

@ -1,23 +0,0 @@
Info<< "Reading thermodynamicProperties\n" << endl;
IOdictionary thermodynamicProperties
(
IOobject
(
"thermodynamicProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar R
(
thermodynamicProperties.lookup("R")
);
dimensionedScalar Cv
(
thermodynamicProperties.lookup("Cv")
);

View file

@ -1,18 +0,0 @@
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar mu
(
transportProperties.lookup("mu")
);

View file

@ -1,138 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
sonicFoamAutoMotion
Description
Transient solver for trans-sonic/supersonic, laminar flow of a
compressible gas with mesh motion..
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "motionSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readThermodynamicProperties.H"
# include "readTransportProperties.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
autoPtr<Foam::motionSolver> motionPtr = motionSolver::New(mesh);
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readPISOControls.H"
# include "compressibleCourantNo.H"
mesh.movePoints(motionPtr->newPoints());
# include "rhoEqn.H"
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
- fvm::laplacian(mu, U)
);
solve(UEqn == -fvc::grad(p));
solve
(
fvm::ddt(rho, e)
+ fvm::div(phi, e)
- fvm::laplacian(mu, e)
==
- p*fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(rho, U))
+ mu*magSqr(symm(fvc::grad(U)))
);
T = e/Cv;
psi = 1.0/(R*T);
// --- PISO loop
for (int corr = 0; corr < nCorr; corr++)
{
U = UEqn.H()/UEqn.A();
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)*
(
(fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)
)
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho/UEqn.A(), p)
);
pEqn.solve();
phi = pEqn.flux();
}
# include "compressibleContinuityErrs.H"
U -= fvc::grad(p)/UEqn.A();
U.correctBoundaryConditions();
}
rho = psi*p;
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

View file

@ -2,19 +2,19 @@
dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
// Pressure bounds
dimensionedScalar pMin("pMin", p.dimensions(), 0);
dimensionedScalar pMax("pMax", p.dimensions(), 0);
dimensionedScalar pMin("pMin", dimPressure, 0);
dimensionedScalar pMax("pMax", dimPressure, GREAT);
fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
// Temperature bounds
dimensionedScalar TMin("TMin", T.dimensions(), 0);
dimensionedScalar TMax("TMax", T.dimensions(), 0);
dimensionedScalar TMin("TMin", dimTemperature, 0);
dimensionedScalar TMax("TMax", dimTemperature, GREAT);
fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
// Velocity bound
dimensionedScalar UMax("UMax", U.dimensions(), 0);
dimensionedScalar UMax("UMax", dimVelocity, GREAT);
fieldBounds.lookup(U.name()) >> UMax.value();
dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

View file

@ -2,19 +2,19 @@
dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
// Pressure bounds
dimensionedScalar pMin("pMin", p.dimensions(), 0);
dimensionedScalar pMax("pMax", p.dimensions(), 0);
dimensionedScalar pMin("pMin", dimPressure, 0);
dimensionedScalar pMax("pMax", dimPressure, GREAT);
fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
// Temperature bounds
dimensionedScalar TMin("TMin", T.dimensions(), 0);
dimensionedScalar TMax("TMax", T.dimensions(), 0);
dimensionedScalar TMin("TMin", dimTemperature, 0);
dimensionedScalar TMax("TMax", dimTemperature, GREAT);
fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
// Velocity bound
dimensionedScalar UMax("UMax", U.dimensions(), 0);
dimensionedScalar UMax("UMax", dimVelocity, GREAT);
fieldBounds.lookup(U.name()) >> UMax.value();
dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

View file

@ -1,19 +1,20 @@
// Read field bounds
dictionary fieldBounds = mesh.solutionDict().subDict("fieldBounds");
// Pressure bounds
dimensionedScalar pMin("pMin", p.dimensions(), 0);
dimensionedScalar pMax("pMax", p.dimensions(), 0);
dimensionedScalar pMin("pMin", dimPressure, 0);
dimensionedScalar pMax("pMax", dimPressure, GREAT);
fieldBounds.lookup(p.name()) >> pMin.value() >> pMax.value();
fieldBounds.lookup("p") >> pMin.value() >> pMax.value();
// Temperature bounds
dimensionedScalar TMin("TMin", T.dimensions(), 0);
dimensionedScalar TMax("TMax", T.dimensions(), 0);
dimensionedScalar TMin("TMin", dimTemperature, 0);
dimensionedScalar TMax("TMax", dimTemperature, GREAT);
fieldBounds.lookup(T.name()) >> TMin.value() >> TMax.value();
fieldBounds.lookup("T") >> TMin.value() >> TMax.value();
// Velocity bound
dimensionedScalar UrelMax("UrelMax", Urel.dimensions(), 0);
dimensionedScalar UMax("UMax", dimVelocity, GREAT);
fieldBounds.lookup(Urel.name()) >> UrelMax.value();
dimensionedScalar smallUrel("smallUrel", dimVelocity, 1e-10);
fieldBounds.lookup(U.name()) >> UMax.value();
dimensionedScalar smallU("smallU", dimVelocity, 1e-10);

View file

@ -1616,7 +1616,7 @@ int main(int argc, char *argv[])
// interior boundaries are handled via faceSets
// cell zones will only be written if there is more than one
if (writeZones && cellGroupZoneID.size()>1)
if (writeZones && cellGroupZoneID.size() > 1)
{
Info<< "Adding Zones" << endl;
List<pointZone*> pz(0);

View file

@ -38,7 +38,7 @@ FoamFile
// Tolerance used in matching faces. Absolute tolerance is span of
// face times this factor. To load incorrectly matches meshes set this
// to a higher value.
matchTolerance 1E-3;
matchTolerance 1e-3;
// Do a synchronisation of coupled points after creation of any patches.
// Note: this does not work with points that are on multiple coupled patches