Skip to content
Snippets Groups Projects
Commit b70ad9e5 authored by Andrew Heather's avatar Andrew Heather
Browse files

Merge branch 'master' into thermo

parents 3e63edf5 53b0d4e7
Branches
Tags
No related merge requests found
Showing
with 166 additions and 155 deletions
......@@ -15,12 +15,10 @@
(
UEqn
==
-fvc::reconstruct
fvc::reconstruct
(
(
fvc::snGrad(pd)
+ ghf*fvc::snGrad(rho)
) * mesh.magSf()
fvc::interpolate(rho)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
);
}
......@@ -59,27 +59,6 @@
fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
Info<< "Creating field pd\n" << endl;
volScalarField pd
(
IOobject
(
"pd",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
p = pd + rho*gh + pRef;
thermo->correct();
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
{
bool closedVolume = pd.needReference();
bool closedVolume = p.needReference();
rho = thermo->rho();
......@@ -17,38 +17,35 @@
)
);
phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
fvScalarMatrix pEqn
(
fvm::ddt(psi, pd)
+ fvc::ddt(psi)*pRef
+ fvc::ddt(psi, rho)*gh
fvm::ddt(psi, p)
+ fvc::div(phi)
- fvm::laplacian(rhorUAf, pd)
- fvm::laplacian(rhorUAf, p)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pdEqn.solve(mesh.solver(pd.name() + "Final"));
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pdEqn.solve(mesh.solver(pd.name()));
pEqn.solve(mesh.solver(p.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi += pdEqn.flux();
phi += pEqn.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
U.correctBoundaryConditions();
p == pd + rho*gh + pRef;
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
#include "rhoEqn.H"
......@@ -62,6 +59,4 @@
/fvc::domainIntegrate(thermo->psi());
rho = thermo->rho();
}
pd == p - (rho*gh + pRef);
}
......@@ -13,14 +13,11 @@
(
UEqn()
==
-fvc::reconstruct
fvc::reconstruct
(
(
fvc::snGrad(pd)
+ ghf*fvc::snGrad(rho)
) * mesh.magSf()
fvc::interpolate(rho)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
).initialResidual();
maxResidual = max(eqnResidual, maxResidual);
......@@ -39,15 +39,14 @@ Description
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readEnvironmentalProperties.H"
#include "createFields.H"
#include "initContinuityErrs.H"
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readEnvironmentalProperties.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
......@@ -55,17 +54,17 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readSIMPLEControls.H"
# include "initConvergenceCheck.H"
#include "readSIMPLEControls.H"
#include "initConvergenceCheck.H"
pd.storePrevIter();
p.storePrevIter();
rho.storePrevIter();
// Pressure-velocity SIMPLE corrector
{
# include "UEqn.H"
# include "hEqn.H"
# include "pEqn.H"
#include "UEqn.H"
#include "hEqn.H"
#include "pEqn.H"
}
turbulence->correct();
......@@ -76,7 +75,7 @@ int main(int argc, char *argv[])
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
# include "convergenceCheck.H"
#include "convergenceCheck.H"
}
Info<< "End\n" << endl;
......
......@@ -51,38 +51,16 @@
)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
Info<< "Creating field pd\n" << endl;
volScalarField pd
(
IOobject
(
"pd",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
p = pd + rho*gh + pRef;
thermo->correct();
label pdRefCell = 0;
scalar pdRefValue = 0.0;
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
pd,
p,
mesh.solutionDict().subDict("SIMPLE"),
pdRefCell,
pdRefValue
pRefCell,
pRefValue
);
......
{
rho = thermo->rho();
volScalarField rUA = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
......@@ -7,60 +9,57 @@
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
phi -= buoyancyPhi;
surfaceScalarField buoyancyPhi =
rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
phi += buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
fvScalarMatrix pEqn
(
fvm::laplacian(rhorUAf, pd) == fvc::div(phi)
fvm::laplacian(rhorUAf, p) == fvc::div(phi)
);
pdEqn.setReference(pdRefCell, pdRefValue);
pEqn.setReference(pRefCell, p[pRefCell]);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pdEqn.solve().initialResidual();
eqnResidual = pEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pdEqn.solve();
pEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
}
// Calculate the conservative fluxes
phi -= pdEqn.flux();
phi -= pEqn.flux();
// Explicitly relax pressure for momentum corrector
pd.relax();
p.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf);
U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
U.correctBoundaryConditions();
}
}
#include "continuityErrs.H"
p == pd + rho*gh + pRef;
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
}
rho = thermo->rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
<< endl;
pd == p - (rho*gh + pRef);
}
......@@ -52,8 +52,24 @@ void Foam::timeActivatedExplicitMulticomponentPointSource::updateAddressing()
forAll(pointSources_, sourceI)
{
const pointSourceProperties& psp = pointSources_[sourceI];
bool foundCell = false;
label cid = mesh_.findCell(psp.location());
if (cid < 0)
if (cid >= 0)
{
foundCell = mesh_.pointInCell(psp.location(), cid);
}
reduce(foundCell, orOp<bool>());
if (!foundCell)
{
label cid = mesh_.findNearestCell(psp.location());
if (cid >= 0)
{
foundCell = mesh_.pointInCell(psp.location(), cid);
}
}
reduce(foundCell, orOp<bool>());
if (!foundCell)
{
FatalErrorIn
(
......@@ -180,7 +196,11 @@ Foam::timeActivatedExplicitMulticomponentPointSource::Su
)
{
const label cid = cellOwners_[sourceI];
sourceField[cid] += dt*psp.fieldData()[i].second()/V[cid];
if (cid >= 0)
{
sourceField[cid] +=
dt*psp.fieldData()[i].second()/V[cid];
}
}
}
}
......@@ -235,7 +255,11 @@ Foam::timeActivatedExplicitMulticomponentPointSource::Su()
)
{
const label cid = cellOwners_[sourceI];
sourceField[cid] += dt*psp.fieldData()[i].second()/V[cid];
if (cid >= 0)
{
sourceField[cid] +=
dt*psp.fieldData()[i].second()/V[cid];
}
}
}
}
......
/*--------------------------------*- 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 alphat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
floor
{
type compressible::alphatWallFunction;
value uniform 0;
}
ceiling
{
type compressible::alphatWallFunction;
value uniform 0;
}
fixedWalls
{
type compressible::alphatWallFunction;
value uniform 0;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......@@ -10,11 +10,12 @@ FoamFile
version 2.0;
format ascii;
class volScalarField;
object pd;
location "0";
object mut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0];
dimensions [1 -1 -1 0 0 0 0];
internalField uniform 0;
......@@ -22,21 +23,20 @@ boundaryField
{
floor
{
type fixedFluxBuoyantPressure;
type mutWallFunction;
value uniform 0;
}
ceiling
{
type fixedFluxBuoyantPressure;
type mutWallFunction;
value uniform 0;
}
fixedWalls
{
type fixedFluxBuoyantPressure;
type mutWallFunction;
value uniform 0;
}
}
// ************************************************************************* //
......@@ -16,26 +16,26 @@ FoamFile
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 100000;
internalField uniform 1e5;
boundaryField
{
floor
{
type calculated;
value uniform 100000;
type fixedFluxBuoyantPressure;
value uniform 1e5;
}
ceiling
{
type calculated;
value uniform 100000;
type fixedFluxBuoyantPressure;
value uniform 1e5;
}
fixedWalls
{
type calculated;
value uniform 100000;
type fixedFluxBuoyantPressure;
value uniform 1e5;
}
}
......
......@@ -42,7 +42,7 @@ laplacianSchemes
{
default none;
laplacian(muEff,U) Gauss linear corrected;
laplacian((rho*(1|A(U))),pd) Gauss linear corrected;
laplacian((rho*(1|A(U))),p) Gauss linear corrected;
laplacian(alphaEff,h) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
......@@ -62,7 +62,7 @@ snGradSchemes
fluxRequired
{
default no;
pd ;
p ;
}
......
......@@ -25,7 +25,7 @@ solvers
relTol 0;
}
pd
p
{
solver PCG;
preconditioner DIC;
......@@ -33,7 +33,7 @@ solvers
relTol 0.1;
}
pdFinal
pFinal
{
solver PCG;
preconditioner DIC;
......
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
......@@ -10,11 +10,12 @@ FoamFile
version 2.0;
format ascii;
class volScalarField;
object pd;
location "0";
object mut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0];
dimensions [1 -1 -1 0 0 0 0];
internalField uniform 0;
......@@ -22,21 +23,20 @@ boundaryField
{
floor
{
type fixedFluxBuoyantPressure;
type mutWallFunction;
value uniform 0;
}
ceiling
{
type fixedFluxBuoyantPressure;
type mutWallFunction;
value uniform 0;
}
fixedWalls
{
type fixedFluxBuoyantPressure;
type mutWallFunction;
value uniform 0;
}
}
// ************************************************************************* //
......@@ -16,26 +16,26 @@ FoamFile
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 100000;
internalField uniform 1e5;
boundaryField
{
floor
{
type calculated;
value uniform 100000;
type fixedFluxBuoyantPressure;
value uniform 1e5;
}
ceiling
{
type calculated;
value uniform 100000;
type fixedFluxBuoyantPressure;
value uniform 1e5;
}
fixedWalls
{
type calculated;
value uniform 100000;
type fixedFluxBuoyantPressure;
value uniform 1e5;
}
}
......
......@@ -41,7 +41,7 @@ laplacianSchemes
{
default none;
laplacian(muEff,U) Gauss linear corrected;
laplacian((rho*(1|A(U))),pd) Gauss linear corrected;
laplacian((rho*(1|A(U))),p) Gauss linear corrected;
laplacian(alphaEff,h) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
......@@ -61,7 +61,7 @@ snGradSchemes
fluxRequired
{
default no;
pd ;
p ;
}
......
......@@ -17,7 +17,7 @@ FoamFile
solvers
{
pd
p
{
solver PCG;
preconditioner DIC;
......@@ -69,16 +69,16 @@ solvers
SIMPLE
{
nNonOrthogonalCorrectors 0;
pdRefCell 0;
pdRefValue 0;
pRefCell 0;
pRefValue 1e5;
}
relaxationFactors
{
rho 1;
pd 0.3;
p 0.3;
U 0.7;
h 0.7;
h 0.1;
k 0.7;
epsilon 0.7;
R 0.7;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment