Commit c103331a authored by Andrew Heather's avatar Andrew Heather
Browse files

INT: integration updates

parent 3603cf28
......@@ -74,4 +74,5 @@ Info<< "Maximum ddtAlpha : " << ddtAlphaNum << endl;
Info<< "Maximum DiffNum : " << DiNum << endl;
// ************************************************************************* //
// Volumatric flux
surfaceScalarField& phi = fluid.phi();
......@@ -114,9 +114,6 @@
p_rgh = p - rho*gh;
}
// Total volumetric flux
surfaceScalarField& phi = fluid.phi();
// Mass flux
surfaceScalarField& rhoPhi = fluid.rhoPhi();
......
......@@ -29,7 +29,8 @@ Group
Description
Solver for n incompressible, non-isothermal immiscible fluids with
phase-change. Uses a VOF (volume of fluid) phase-fraction based interface capturing approach.
phase-change. Uses a VOF (volume of fluid) phase-fraction based interface
capturing approach.
The momentum, energy and other fluid properties are of the "mixture" and a
single momentum equation is solved.
......@@ -54,6 +55,8 @@ Description
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
......@@ -61,12 +64,15 @@ int main(int argc, char *argv[])
pimpleControl pimple(mesh);
#include "createFields.H"
#include "createFieldRefs.H"
#include "createFvOptions.H"
#include "createTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
//#include "alphaCourantNo.H"
#include "setInitialDeltaT.H"
turbulence->validate();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
......@@ -92,6 +98,7 @@ int main(int argc, char *argv[])
#include "UEqn.H"
#include "YEqns.H"
#include "TEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
......@@ -108,9 +115,7 @@ int main(int argc, char *argv[])
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
......
......@@ -62,9 +62,6 @@ class DTRMParticle
:
public particle
{
private:
// Private data
//- Size in bytes of the fields
......@@ -98,8 +95,6 @@ public:
:
public particle::trackingData
{
// Interpolators for continuous phase fields
const interpolationCell<scalar>& aInterp_;
......@@ -135,10 +130,6 @@ public:
volScalarField& Q
);
// Public data
// Member functions
......@@ -151,7 +142,6 @@ public:
inline const UPtrList<reflectionModel>& reflection() const;
inline scalar& Q(label celli);
};
// Static data members
......@@ -268,7 +258,6 @@ public:
// Edit
//- Return access to the target position
inline point& p1();
......@@ -293,7 +282,6 @@ public:
// Member Operators
//- Overridable function to handle the particle hitting a processorPatch
void hitProcessorPatch
(
......
......@@ -60,7 +60,6 @@ Foam::DTRMParticle::DTRMParticle
if (is.format() == IOstream::ASCII)
{
is >> p0_ >> p1_ >> I0_ >> I_ >> dA_ >> transmissiveId_;
DebugVar(transmissiveId_);
}
else
{
......@@ -95,7 +94,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const DTRMParticle& p)
}
// Check state of Ostream
os.check("Ostream& operator<<(Ostream&, const DTRMParticle&)");
os.check(FUNCTION_NAME);
return os;
}
......
......@@ -54,20 +54,17 @@ namespace Foam
0
);
namespace radiation
template<>
const char* Foam::NamedEnum
<
Foam::radiation::laserDTRM::powerDistributionMode,
3
>::names[] =
{
template<>
const char* Foam::NamedEnum
<
Foam::radiation::laserDTRM::powerDistributionMode,
3
>::names[] =
{
"Gaussian",
"manual",
"uniform"
};
}
"Gaussian",
"manual",
"uniform"
};
}
const Foam::NamedEnum
......@@ -82,35 +79,32 @@ const Foam::NamedEnum
Foam::scalar Foam::radiation::laserDTRM::calculateIp(scalar r, scalar theta)
{
const scalar t = mesh_.time().value();
const scalar power = laserPower_->value(t);
switch(mode_)
{
case pdGaussian:
{
scalar I0 =
laserPower_->value(t)/(mathematical::twoPi*sqr(sigma_));
return(I0*exp(-sqr(r)/2.0/sqr(sigma_)));
scalar I0 = power/(mathematical::twoPi*sqr(sigma_));
return I0*exp(-sqr(r)/2.0/sqr(sigma_));
break;
}
case pdManual:
{
return(laserPower_->value(t)*powerDistribution_()(theta, r));
return power*powerDistribution_()(theta, r);
break;
}
case pdUniform:
{
return
(
laserPower_->value(t)/(mathematical::pi*sqr(focalLaserRadius_))
);
return power/(mathematical::pi*sqr(focalLaserRadius_));
break;
}
default:
{
FatalErrorInFunction
<< "Unhandled type " << powerDistypeNames_
<< abort(FatalError);
return(0);
<< "Unhandled type " << powerDistypeNames_
<< abort(FatalError);
return (0);
}
}
}
......@@ -125,7 +119,7 @@ Foam::tmp<Foam::volVectorField> Foam::radiation::laserDTRM::nHatfv
const dimensionedScalar deltaN
(
"deltaN",
1e-7/pow(average(mesh_.V()), 1.0/3.0)
1e-7/cbrt(average(mesh_.V()))
);
const volVectorField gradAlphaf
......@@ -134,8 +128,8 @@ Foam::tmp<Foam::volVectorField> Foam::radiation::laserDTRM::nHatfv
- alpha1*fvc::grad(alpha2)
);
// Face unit interface normal
return gradAlphaf/(mag(gradAlphaf)+ deltaN);
// Face unit interface normal
return gradAlphaf/(mag(gradAlphaf)+ deltaN);
}
......@@ -203,7 +197,7 @@ void Foam::radiation::laserDTRM::initialise()
nParticles_ = ndr_*ndTheta_;
switch(mode_)
switch (mode_)
{
case pdGaussian:
{
......@@ -285,8 +279,8 @@ void Foam::radiation::laserDTRM::initialise()
if (cellI != -1)
{
// Create a new particle
DTRMParticle* pPtr = new DTRMParticle
(mesh_, p0, p1, Ip, cellI, dAi, -1);
DTRMParticle* pPtr =
new DTRMParticle(mesh_, p0, p1, Ip, cellI, dAi, -1);
// Add to cloud
DTRMCloud_.addParticle(pPtr);
......@@ -294,7 +288,7 @@ void Foam::radiation::laserDTRM::initialise()
if (returnReduce(cellI, maxOp<label>()) == -1)
{
WarningIn("void Foam::radiation::laserDTRM::initialise()")
WarningInFunction
<< "Cannot find owner cell for particle at position " << p0
<< endl;
}
......@@ -303,18 +297,16 @@ void Foam::radiation::laserDTRM::initialise()
}
else
{
FatalErrorIn("void Foam::radiation::laserDTRM::initialise()")
FatalErrorInFunction
<< "Current functionality limited to 3-D cases"
<< exit(FatalError);
}
if (debug)
{
Info << "Total Power in the laser : " << power << endl;
Info << "Total Area in the laser : " << area << endl;
Info << "Number of particles in the laser : "
<< returnReduce(DTRMCloud_.size(), sumOp<label>()) << endl;
}
DebugInfo
<< "Total Power in the laser : " << power << endl
<< "Total Area in the laser : " << area << endl
<< "Number of particles in the laser : "
<< returnReduce(DTRMCloud_.size(), sumOp<label>()) << endl;
}
......@@ -534,12 +526,6 @@ Foam::radiation::laserDTRM::laserDTRM
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::radiation::laserDTRM::~laserDTRM()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::radiation::laserDTRM::read()
......@@ -661,8 +647,6 @@ void Foam::radiation::laserDTRM::calculate()
{
nHat[cellI] += nHatPhase[cellI];
}
}
}
reflectionModelId++;
......@@ -687,8 +671,8 @@ void Foam::radiation::laserDTRM::calculate()
Q_
);
Info << "Move particles..."
<< returnReduce(DTRMCloud_.size(), sumOp<label>()) << endl;
Info<< "Move particles..."
<< returnReduce(DTRMCloud_.size(), sumOp<label>()) << endl;
DTRMCloud_.move(DTRMCloud_, td, mesh_.time().deltaTValue());
......@@ -697,13 +681,10 @@ void Foam::radiation::laserDTRM::calculate()
if (debug)
{
Info<< "Final number of particles..."
<< returnReduce(DTRMCloud_.size(), sumOp<label>()) << endl;
Info<< "Final number of particles..."
<< returnReduce(DTRMCloud_.size(), sumOp<label>()) << endl;
OFstream osRef
(
type() + ":particlePath.obj"
);
OFstream osRef(type() + ":particlePath.obj");
label vertI = 0;
List<pointField> positions(Pstream::nProcs());
......@@ -727,7 +708,7 @@ void Foam::radiation::laserDTRM::calculate()
Pstream::gatherList(p0);
Pstream::scatterList(p0);
for (label proci = 0; proci < Pstream::nProcs(); proci++)
for (label proci = 0; proci < Pstream::nProcs(); ++proci)
{
const pointField& pos = positions[proci];
const pointField& pfinal = p0[proci];
......@@ -748,8 +729,8 @@ void Foam::radiation::laserDTRM::calculate()
if (mesh_.time().outputTime())
{
reflectingCellsVol.write();
nHat.write();
reflectingCellsVol.write();
nHat.write();
}
}
......@@ -762,32 +743,29 @@ void Foam::radiation::laserDTRM::calculate()
Foam::tmp<Foam::volScalarField> Foam::radiation::laserDTRM::Rp() const
{
return tmp<volScalarField>
return tmp<volScalarField>::New
(
new volScalarField
IOobject
(
IOobject
(
"zero",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
"zero",
mesh_.time().timeName(),
mesh_,
dimensionedScalar
(
"zero",
dimPower/dimVolume/pow4(dimTemperature),
0.0
)
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedScalar
(
"zero",
dimPower/dimVolume/pow4(dimTemperature),
0.0
)
);
}
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
Foam::radiation::laserDTRM::Ru() const
{
return Q_.internalField();
......
......@@ -237,7 +237,7 @@ public:
//- Destructor
virtual ~laserDTRM();
virtual ~laserDTRM() = default;
// Member functions
......
......@@ -52,11 +52,7 @@ Foam::radiation::localDensityAbsorptionEmission::alpha(word alphaName) const
{
if (!mesh_.foundObject<volScalarField>(alphaName))
{
FatalErrorIn
(
"const Foam::volScalarField& "
"Foam::radiation::localDensityAbsorptionEmission::alpha() const"
)
FatalErrorInFunction
<< "Unable to retrieve density field " << alphaName << " from "
<< "database. Available objects:" << mesh_.sortedNames()
<< exit(FatalError);
......@@ -83,13 +79,6 @@ Foam::radiation::localDensityAbsorptionEmission::localDensityAbsorptionEmission
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::radiation::localDensityAbsorptionEmission::
~localDensityAbsorptionEmission()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
......
......@@ -33,8 +33,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef radiationLocalDensityAbsorptionEmission_H
#define radiationLocalDensityAbsorptionEmission_H
#ifndef radiation_localDensityAbsorptionEmission_H
#define radiation_localDensityAbsorptionEmission_H
#include "absorptionEmissionModel.H"
......@@ -94,7 +94,7 @@ public:
//- Destructor
virtual ~localDensityAbsorptionEmission();
virtual ~localDensityAbsorptionEmission() = default;
// Member Functions
......
......@@ -119,7 +119,7 @@ Foam::scalar Foam::radiation::Fresnel::rho
(sqr(sqrt(sqrP) + sqrt(n1)*sinTheta1*tanTheta1) + sqrQ)
)*rhoP;
return (rhoP + rhoN)/2;
return 0.5*(rhoP + rhoN);
}
......@@ -133,5 +133,4 @@ Foam::vector Foam::radiation::Fresnel::R
}
// ************************************************************************* //
......@@ -35,8 +35,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef Fresnel_H
#define Fresnel_H
#ifndef radiation_Fresnel_H
#define radiation_Fresnel_H
#include "reflectionModel.H"
......@@ -48,14 +48,13 @@ namespace radiation
{
/*---------------------------------------------------------------------------*\
Class Fresnel Declaration
Class Fresnel Declaration
\*---------------------------------------------------------------------------*/
class Fresnel
:
public reflectionModel
{
// Private data
//- Coefficients dictionary
......@@ -73,12 +72,8 @@ public:
//- Runtime type information
TypeName("Fresnel");
// Constructors
//- Construct from components
Fresnel(const dictionary& dict, const fvMesh& mesh);
//- Construct from components
Fresnel(const dictionary& dict, const fvMesh& mesh);
//- Destructor
virtual ~Fresnel();
......@@ -90,11 +85,8 @@ public:
virtual vector R(const vector& incident, const vector& n) const;
//- Return reflectivity from medium1 to medium2 and a incident angle.
// nk1 = (n1 - i k1) from medium 1.
virtual scalar rho
(
const scalar incidentAngle
) const;
// nk1 = (n1 - i k1) from medium 1.
virtual scalar rho(const scalar incidentAngle) const;
};
......
......@@ -57,12 +57,6 @@ Foam::radiation::FresnelLaser::FresnelLaser
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::radiation::FresnelLaser::~FresnelLaser()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::radiation::FresnelLaser::rho
......@@ -75,7 +69,7 @@ Foam::scalar Foam::radiation::FresnelLaser::rho
scalar rho =
0.5
* (
(1 + sqr(1 - epsilon_*cosTheta))/(1 + sqr(1 + epsilon_*cosTheta))
(1 + sqr(1 - epsilon_*cosTheta))/(1 + sqr(1 + epsilon_*cosTheta))
+
(sqr(epsilon_) - 2*epsilon_*cosTheta + 2*sqr(cosTheta))
/
......
......@@ -39,8 +39,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef FresnelLaser_H
#define FresnelLaser_H
#ifndef radiation_FresnelLaser_H
#define radiation_FresnelLaser_H
#include "reflectionModel.H"
......@@ -71,15 +71,11 @@ public:
//- Runtime type information
TypeName("FresnelLaser");
// Constructors