Commit f97160ba authored by henry's avatar henry
Browse files

Choosing Andy's

parents d29e7288 372cd1b0
label Nparcels = dieselSpray.size();
reduce(Nparcels, sumOp<label>());
Info<< "\nNumber of parcels in system | "
Info<< "\nNumber of parcels in system.... | "
<< Nparcels << endl
<< "Injected liquid mass....... | "
<< "Injected liquid mass........... | "
<< 1e6*dieselSpray.injectedMass(runTime.value()) << " mg" << endl
<< "Liquid Mass in system...... | "
<< "Liquid Mass in system.......... | "
<< 1e6*dieselSpray.liquidMass() << " mg" << endl
<< "SMD, Dmax.................. | "
<< "SMD, Dmax...................... | "
<< dieselSpray.smd()*1e6 << " mu, "
<< dieselSpray.maxD()*1e6 << " mu"
<< endl;
scalar evapMass =
dieselSpray.injectedMass(runTime.value()) - dieselSpray.liquidMass();
scalar evapMass =
dieselSpray.injectedMass(runTime.value())
- dieselSpray.liquidMass();
scalar gasMass = fvc::domainIntegrate(rho).value();
......@@ -26,6 +27,6 @@
scalar addedMass = gasMass - gasMass0;
Info<< "Added gas mass = " << 1e6*addedMass << " mg" << nl
<< "Evaporation Continuity Error| "
Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"
<< nl << "Evaporation Continuity Error... | "
<< 1e6*(addedMass - evapMass) << " mg" << endl;
......@@ -31,6 +31,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "cyclicPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
......@@ -75,26 +76,54 @@ int main(int argc, char *argv[])
}
// Give patch area
Info<< " Patch area = " << gSum(mesh.Sf().boundaryField()[patchi]) << endl;
if (isType<cyclicPolyPatch>(mesh.boundaryMesh()[patchi]))
{
Info<< " Cyclic patch area: " << nl;
label nFaces = mesh.boundaryMesh()[patchi].size();
vector sum1 = vector::zero;
vector sum2 = vector::zero;
for (label i=0; i<nFaces/2; i++)
{
sum1 += mesh.Sf().boundaryField()[patchi][i];
sum2 += mesh.Sf().boundaryField()[patchi][i+nFaces/2];
}
reduce(sum1, sumOp<vector>());
reduce(sum2, sumOp<vector>());
Info<< " - half 1 = " << sum1 << ", " << mag(sum1) << nl
<< " - half 2 = " << sum2 << ", " << mag(sum2) << nl
<< " - total = " << (sum1 + sum2) << ", "
<< mag(sum1 + sum2) << endl;;
}
else
{
Info<< " Patch area = "
<< gSum(mesh.Sf().boundaryField()[patchi]) << endl;
}
if (fieldHeader.headerClassName() == "volScalarField")
// Read field and calc integral
if (fieldHeader.headerClassName() == volScalarField::typeName)
{
Info<< " Reading volScalarField " << fieldName << endl;
volScalarField field(fieldHeader, mesh);
Info<< " Reading " << volScalarField::typeName << " "
<< fieldName << endl;
volScalarField field(fieldHeader, mesh);
vector sumField = gSum
(
mesh.Sf().boundaryField()[patchi]
* field.boundaryField()[patchi]
*field.boundaryField()[patchi]
);
Info<< " Integral of " << fieldName << " over patch "
<< patchName << '[' << patchi << ']' << " = "
<< sumField << nl;
}
else if (fieldHeader.headerClassName() == "surfaceScalarField")
else if
(
fieldHeader.headerClassName() == surfaceScalarField::typeName
)
{
Info<< " Reading surfaceScalarField " << fieldName << endl;
Info<< " Reading " << surfaceScalarField::typeName << " "
<< fieldName << endl;
surfaceScalarField field(fieldHeader, mesh);
scalar sumField = gSum(field.boundaryField()[patchi]);
......@@ -106,8 +135,10 @@ int main(int argc, char *argv[])
else
{
FatalError
<< "Only possible to integrate volScalarFields "
<< "and surfaceScalarFields" << nl << exit(FatalError);
<< "Only possible to integrate "
<< volScalarField::typeName << "s "
<< "and " << surfaceScalarField::typeName << "s"
<< nl << exit(FatalError);
}
}
else
......
......@@ -53,7 +53,6 @@ void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
vsf0.value()
);
vsf.correctBoundaryConditions();
vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
}
}
......
......@@ -27,7 +27,7 @@ Global
Description
Reset the timestep to maintain a constant maximum courant Number.
Reduction of time-step is imediate but increase is damped to avoid
Reduction of time-step is immediate, but increase is damped to avoid
unstable oscillations.
\*---------------------------------------------------------------------------*/
......@@ -45,7 +45,7 @@ if (adjustTimeStep)
maxDeltaT
)
);
Info<< "deltaT = " << runTime.deltaT().value() << endl;
}
......
......@@ -33,7 +33,7 @@ Description
if (adjustTimeStep)
{
if (CoNum > SMALL)
if ((runTime.timeIndex() == 0) && (CoNum > SMALL))
{
runTime.setDeltaT
(
......
......@@ -30,19 +30,19 @@ License
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(unitInjector, 0);
addToRunTimeSelectionTable
(
injectorType,
unitInjector,
dictionary
);
defineTypeNameAndDebug(unitInjector, 0);
addToRunTimeSelectionTable
(
injectorType,
unitInjector,
dictionary
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -74,26 +74,36 @@ Foam::unitInjector::unitInjector
// check if time entries for soi and eoi match
if (mag(massFlowRateProfile_[0][0]-TProfile_[0][0]) > SMALL)
{
FatalError << "unitInjector::unitInjector(const time& t, const dictionary dict) " << endl
<< " start-times do not match for TemperatureProfile and massFlowRateProfile."
<< abort(FatalError);
FatalErrorIn
(
"unitInjector::unitInjector(const time& t, const dictionary dict)"
)<< "start-times do not match for TemperatureProfile and "
<< " massFlowRateProfile." << nl << exit (FatalError);
}
if (mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]-TProfile_[TProfile_.size()-1][0]) > SMALL)
if
(
mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]
- TProfile_[TProfile_.size()-1][0])
> SMALL
)
{
FatalError << "unitInjector::unitInjector(const time& t, const dictionary dict) " << endl
<< " end-times do not match for TemperatureProfile and massFlowRateProfile."
<< abort(FatalError);
FatalErrorIn
(
"unitInjector::unitInjector(const time& t, const dictionary dict)"
)<< "end-times do not match for TemperatureProfile and "
<< "massFlowRateProfile." << nl << exit(FatalError);
}
// convert CA to real time
forAll(massFlowRateProfile_, i)
{
massFlowRateProfile_[i][0] = t.userTimeToTime(massFlowRateProfile_[i][0]);
massFlowRateProfile_[i][0] =
t.userTimeToTime(massFlowRateProfile_[i][0]);
velocityProfile_[i][0] = massFlowRateProfile_[i][0];
injectionPressureProfile_[i][0] = massFlowRateProfile_[i][0];
}
forAll(TProfile_, i)
{
TProfile_[i][0] = t.userTimeToTime(TProfile_[i][0]);
......@@ -105,14 +115,14 @@ Foam::unitInjector::unitInjector
{
// correct the massFlowRateProfile to match the injected mass
massFlowRateProfile_[i][1] *= mass_/integratedMFR;
CdProfile_[i][0] = massFlowRateProfile_[i][0];
CdProfile_[i][1] = Cd_;
}
// Normalize the direction vector
direction_ /= mag(direction_);
setTangentialVectors();
// check molar fractions
......@@ -124,9 +134,9 @@ Foam::unitInjector::unitInjector
if (mag(Xsum - 1.0) > SMALL)
{
Info << "Warning!!!\n unitInjector::unitInjector(const time& t, Istream& is)"
<< "X does not add up to 1.0, correcting molar fractions."
<< endl;
WarningIn("unitInjector::unitInjector(const time& t, Istream& is)")
<< "X does not sum to 1.0, correcting molar fractions."
<< nl << endl;
forAll(X_, i)
{
X_[i] /= Xsum;
......@@ -169,18 +179,18 @@ Foam::label Foam::unitInjector::nParcelsToInject
const scalar time1
) const
{
scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
label nParcels = label(mInj/averageParcelMass_ + 0.49);
return nParcels;
}
const Foam::vector Foam::unitInjector::position(const label n) const
{
return position_;
}
Foam::vector Foam::unitInjector::position
(
const label n,
......@@ -212,7 +222,7 @@ Foam::vector Foam::unitInjector::position
scalar iAngle = 2.0*mathematicalConstant::pi*rndGen.scalar01();
return
(
(
position_
+ iRadius
* (
......@@ -220,22 +230,25 @@ Foam::vector Foam::unitInjector::position
+ tangentialInjectionVector2_*sin(iAngle)
)
);
}
return position_;
}
Foam::label Foam::unitInjector::nHoles() const
{
return 1;
}
Foam::scalar Foam::unitInjector::d() const
{
return d_;
}
const Foam::vector& Foam::unitInjector::direction
(
const label i,
......@@ -245,6 +258,7 @@ const Foam::vector& Foam::unitInjector::direction
return direction_;
}
Foam::scalar Foam::unitInjector::mass
(
const scalar time0,
......@@ -255,7 +269,7 @@ Foam::scalar Foam::unitInjector::mass
{
scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
// correct mass if calculation is 2D
// correct mass if calculation is 2D
if (twoD)
{
mInj *= 0.5*angleOfWedge/mathematicalConstant::pi;
......@@ -264,82 +278,80 @@ Foam::scalar Foam::unitInjector::mass
return mInj;
}
Foam::scalar Foam::unitInjector::mass() const
{
return mass_;
}
const Foam::scalarField& Foam::unitInjector::X() const
{
return X_;
}
Foam::List<Foam::unitInjector::pair> Foam::unitInjector::T() const
{
return TProfile_;
}
Foam::scalar Foam::unitInjector::T(const scalar time) const
{
return getTableValue(TProfile_, time);
}
Foam::scalar Foam::unitInjector::tsoi() const
{
return massFlowRateProfile_[0][0];
}
Foam::scalar Foam::unitInjector::teoi() const
{
return massFlowRateProfile_[massFlowRateProfile_.size()-1][0];
}
Foam::scalar Foam::unitInjector::massFlowRate
(
const scalar time
) const
Foam::scalar Foam::unitInjector::massFlowRate(const scalar time) const
{
return getTableValue(massFlowRateProfile_, time);
}
Foam::scalar Foam::unitInjector::injectionPressure
(
const scalar time
) const
Foam::scalar Foam::unitInjector::injectionPressure(const scalar time) const
{
return getTableValue(injectionPressureProfile_, time);
}
Foam::scalar Foam::unitInjector::velocity
(
const scalar time
) const
Foam::scalar Foam::unitInjector::velocity(const scalar time) const
{
return getTableValue(velocityProfile_, time);
}
Foam::List<Foam::unitInjector::pair> Foam::unitInjector::CdProfile() const
{
return CdProfile_;
}
Foam::scalar Foam::unitInjector::Cd
(
const scalar time
) const
Foam::scalar Foam::unitInjector::Cd(const scalar time) const
{
return Cd_;
}
Foam::scalar Foam::unitInjector::fractionOfInjection(const scalar time) const
{
return integrateTable(massFlowRateProfile_, time)/mass_;
}
Foam::scalar Foam::unitInjector::injectedMass
(
const scalar t
) const
Foam::scalar Foam::unitInjector::injectedMass(const scalar t) const
{
return mass_*fractionOfInjection(t);
}
......@@ -351,7 +363,6 @@ void Foam::unitInjector::correctProfiles
const scalar referencePressure
)
{
scalar A = 0.25*mathematicalConstant::pi*pow(d_, 2.0);
scalar pDummy = 1.0e+5;
......@@ -365,14 +376,17 @@ void Foam::unitInjector::correctProfiles
}
}
Foam::vector Foam::unitInjector::tan1(const label n) const
Foam::vector Foam::unitInjector::tan1(const label) const
{
return tangentialInjectionVector1_;
}
Foam::vector Foam::unitInjector::tan2(const label n) const
Foam::vector Foam::unitInjector::tan2(const label) const
{
return tangentialInjectionVector2_;
}
// ************************************************************************* //
......@@ -63,14 +63,10 @@ void spray::inject()
{
Np = max(1, Np);
scalar mp = mass/Np/nHoles;
// constT is only larger than zero for the first
// constT is only larger than zero for the first
// part of the injection
scalar constT = max
(
0.0,
it->tsoi() - time0
);
scalar constT = max(0.0, it->tsoi() - time0);
// deltaT is the duration of injection during this timestep
scalar deltaT = min
......@@ -103,9 +99,10 @@ void spray::inject()
axisOfWedgeNormal_,
rndGen_
);
scalar diameter = injection().d0(i, toi);
vector direction = injection().direction(i, n, toi, diameter);
vector direction =
injection().direction(i, n, toi, diameter);
vector U = injection().velocity(i, toi)*direction;
scalar symComponent = direction & axisOfSymmetry_;
......@@ -117,13 +114,13 @@ void spray::inject()
scalar ddev = breakup().yDot0();
label injectorCell = mesh_.findCell(injectionPosition);
# include "findInjectorCell.H"
if (injectorCell >= 0)
{
scalar liquidCore = 1.0;
// construct the parcel that is to be injected
parcel* pPtr = new parcel
......@@ -148,19 +145,16 @@ void spray::inject()
fuels_->components()
);
injectedLiquidKE_ += 0.5*pPtr->m()*pow(mag(U), 2.0);
injectedLiquidKE_ += 0.5*pPtr->m()*magSqr(U);
scalar dt = time - toi;
pPtr->stepFraction() =
(runTime_.deltaT().value() - dt)
/runTime_.deltaT().value();
bool keepParcel = pPtr->move
(
*this
);
bool keepParcel = pPtr->move(*this);
if (keepParcel)
{
addParticle(pPtr);
......
......@@ -65,7 +65,7 @@ void spray::evolve()
inject();
atomizationLoop();
breakupLoop();
UInterpolator_.clear();
rhoInterpolator_.clear();
pInterpolator_.clear();
......@@ -89,12 +89,7 @@ void spray::move()
void spray::breakupLoop()
{
for
(
spray::iterator elmnt = begin();
elmnt != end();
++elmnt
)
forAllIter(spray::iterator, *this, elmnt)
{
// interpolate...
vector velocity = UInterpolator().interpolate
......@@ -128,12 +123,7 @@ void spray::breakupLoop()
void spray::atomizationLoop()
{
for
(
spray::iterator elmnt = begin();
elmnt != end();
++elmnt
)
forAllIter(spray::iterator, *this, elmnt)
{
// interpolate...
vector velocity = UInterpolator().interpolate
......
......@@ -108,31 +108,19 @@ void reitzKHRT::breakupParcel
// frequency of the fastest growing KH-wave
scalar omegaKH =
(
0.34 + 0.38*pow(weGas, 1.5)
)/
(
(1 + ohnesorge)*(1 + 1.4*pow(taylor, 0.6))
)*sqrt(sigma/(rhoLiquid*pow(r, 3)) );
(0.34 + 0.38*pow(weGas, 1.5))
/((1 + ohnesorge)*(1 + 1.4*pow(taylor, 0.6)))
*sqrt(sigma/(rhoLiquid*pow(r, 3)));
// ... and the corresponding KH wave-length.
// corresponding KH wave-length.
scalar lambdaKH =
9.02*r*
(
1.0 + 0.45*sqrt(ohnesorge)
)*
(
1.0 + 0.4*pow(taylor, 0.7)
)/
pow
(
(
1.0 + 0.865*pow(weGas, 1.67)
),
0.6