Adding leapfrog integration parts to the move method.

This commit is contained in:
graham
2009-09-08 21:07:22 +01:00
parent 2d017d6cc4
commit def1e7b59f
6 changed files with 117 additions and 48 deletions

View File

@ -1620,7 +1620,7 @@ void Foam::ReferredCellList<ParticleType>::referParticles
const List<DynamicList<ParticleType*> >& cellOccupancy
)
{
Info<< "Refer particles" << endl;
Info<< " Refer particles" << endl;
// Clear all existing referred particles

View File

@ -254,23 +254,25 @@ void Foam::InteractingKinematicCloud<ParcelType>::evolve()
resetSourceTerms();
}
const scalar deltaT = mesh_.time().deltaT().value();
forAllIter(typename InteractingKinematicCloud<ParcelType>, *this, iter)
{
ParcelType& p = iter();
p.U() += 0.5*deltaT*p.f()/p.mass();
}
// Sympletic leapfrog integration of particle forces:
// + apply half deltaV with stored force
// + move positions with new velocity
// + calculate forces in new position
// + apply half deltaV with new force
td.part() = ParcelType::trackData::LEAPFROG_VELOCITY_STEP;
Cloud<ParcelType>::move(td);
td.part() = ParcelType::trackData::LINEAR_TRACK;
Cloud<ParcelType>::move(td);
// td.part() = ParcelType::trackData::ROTATIONAL_TRACK;
// Cloud<ParcelType>::move(td);
this->collision().collide();
forAllIter(typename InteractingKinematicCloud<ParcelType>, *this, iter)
{
ParcelType& p = iter();
p.U() += 0.5*deltaT*p.f()/p.mass();
}
td.part() = ParcelType::trackData::LEAPFROG_VELOCITY_STEP;
Cloud<ParcelType>::move(td);
postEvolve();
}

View File

@ -229,48 +229,73 @@ bool Foam::InteractingKinematicParcel<ParcelType>::move(TrackData& td)
const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
const scalar deltaT = mesh.time().deltaT().value();
scalar tEnd = (1.0 - p.stepFraction())*deltaT;
const scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
if (td.part() == TrackData::LEAPFROG_VELOCITY_STEP)
{
// Apply correction to position for reduced-D cases
meshTools::constrainToMeshCentre(mesh, p.position());
// First and last leapfrog velocity adjust part, required
// before and after tracking and force calculation
// Set the Lagrangian time-step
scalar dt = min(dtMax, tEnd);
p.U() += 0.5*deltaT*p.f()/p.mass();
}
else if (td.part() == TrackData::LINEAR_TRACK)
{
scalar tEnd = (1.0 - p.stepFraction())*deltaT;
const scalar dtMax = tEnd;
// Remember which cell the Parcel is in since this will change if a
// face is hit
label cellI = p.cell();
dt *= p.trackToFace(p.position() + dt*U_, td);
tEnd -= dt;
p.stepFraction() = 1.0 - tEnd/deltaT;
// Avoid problems with extremely small timesteps
if (dt > ROOTVSMALL)
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{
// Update cell based properties
p.setCellValues(td, dt, cellI);
// Apply correction to position for reduced-D cases
meshTools::constrainToMeshCentre(mesh, p.position());
if (td.cloud().cellValueSourceCorrection())
// Set the Lagrangian time-step
scalar dt = min(dtMax, tEnd);
// Remember which cell the Parcel is in since this will change if a
// face is hit
label cellI = p.cell();
dt *= p.trackToFace(p.position() + dt*U_, td);
tEnd -= dt;
p.stepFraction() = 1.0 - tEnd/deltaT;
// Avoid problems with extremely small timesteps
if (dt > ROOTVSMALL)
{
p.cellValueSourceCorrection(td, dt, cellI);
// Update cell based properties
p.setCellValues(td, dt, cellI);
if (td.cloud().cellValueSourceCorrection())
{
p.cellValueSourceCorrection(td, dt, cellI);
}
p.calc(td, dt, cellI);
}
p.calc(td, dt, cellI);
}
if (p.onBoundary() && td.keepParticle)
{
if (isType<processorPolyPatch>(pbMesh[p.patch(p.face())]))
if (p.onBoundary() && td.keepParticle)
{
td.switchProcessor = true;
if (isType<processorPolyPatch>(pbMesh[p.patch(p.face())]))
{
td.switchProcessor = true;
}
}
}
}
else if (td.part() == TrackData::ROTATIONAL_TRACK)
{
Info<< "No rotational tracking implementation" << endl;
}
else
{
FatalErrorIn
(
"InteractingKinematicParcel<ParcelType>::move(TrackData& td)"
)
<< td.part()
<< " is an invalid part of the tracking method."
<< abort(FatalError);
}
return td.keepParticle;
}

View File

@ -125,6 +125,18 @@ public:
:
public Particle<ParcelType>::trackData
{
public:
enum trackPart
{
LEAPFROG_VELOCITY_STEP,
LINEAR_TRACK,
ROTATIONAL_TRACK
};
private:
// Private data
//- Reference to the cloud containing this particle
@ -133,7 +145,6 @@ public:
//- Particle constant properties
const constantProperties& constProps_;
// Interpolators for continuous phase fields
//- Density interpolator
@ -148,6 +159,10 @@ public:
//- Local gravitational or other body-force acceleration
const vector& g_;
// label specifying which part of the integration
// algorithm is taking place
trackPart part_;
public:
@ -161,7 +176,8 @@ public:
const interpolation<scalar>& rhoInterp,
const interpolation<vector>& UInterp,
const interpolation<scalar>& muInterp,
const vector& g
const vector& g,
trackPart part = LINEAR_TRACK
);
@ -185,8 +201,14 @@ public:
// phase dynamic viscosity field
inline const interpolation<scalar>& muInterp() const;
// Return const access to the gravitational acceleration vector
//- Return const access to the gravitational acceleration vector
inline const vector& g() const;
//- Return the part of the tracking operation taking place
inline trackPart part() const;
//- Return access to the part of the tracking operation taking place
inline trackPart& part();
};

View File

@ -55,7 +55,8 @@ inline Foam::InteractingKinematicParcel<ParcelType>::trackData::trackData
const interpolation<scalar>& rhoInterp,
const interpolation<vector>& UInterp,
const interpolation<scalar>& muInterp,
const vector& g
const vector& g,
trackPart part
)
:
Particle<ParcelType>::trackData(cloud),
@ -64,7 +65,8 @@ inline Foam::InteractingKinematicParcel<ParcelType>::trackData::trackData
rhoInterp_(rhoInterp),
UInterp_(UInterp),
muInterp_(muInterp),
g_(g)
g_(g),
part_(part)
{}
@ -205,6 +207,24 @@ Foam::InteractingKinematicParcel<ParcelType>::trackData::g() const
}
template<class ParcelType>
inline
typename Foam::InteractingKinematicParcel<ParcelType>::trackData::trackPart
Foam::InteractingKinematicParcel<ParcelType>::trackData::part() const
{
return part_;
}
template<class ParcelType>
inline
typename Foam::InteractingKinematicParcel<ParcelType>::trackData::trackPart&
Foam::InteractingKinematicParcel<ParcelType>::trackData::part()
{
return part_;
}
// * * * * * * * InteractingKinematicParcel Member Functions * * * * * * * //
template <class ParcelType>

View File

@ -32,7 +32,7 @@ License
template<class CloudType>
void Foam::PairCollision<CloudType>::buildCellOccupancy()
{
Info<< "Build cell occupancy" << endl;
Info<< " Build cell occupancy" << endl;
forAll(cellOccupancy_, cO)
{