lagrangian: Corrected patch data
The KinematicCloud::patchData method has been made consistent on moving meshes and/or when the time-step is being sub-cycled. It has also been altered to calculate the normal component of a moving patch's velocity directly from the point motions. This prevents an infinite loop occuring due to inconsistency between the velocity used to calculate a rebound and that used when tracking. Some minor style improvements to the particle class have also been made.
This commit is contained in:
@ -63,7 +63,7 @@ bool Foam::DSMCParcel<ParcelType>::move(TrackData& td, const scalar trackTime)
|
||||
|
||||
if (p.onBoundaryFace() && td.keepParticle)
|
||||
{
|
||||
if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
|
||||
if (isA<processorPolyPatch>(pbMesh[p.patch()]))
|
||||
{
|
||||
td.switchProcessor = true;
|
||||
}
|
||||
|
||||
@ -40,7 +40,7 @@ namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::particle::tetGeometry
|
||||
void Foam::particle::stationaryTetGeometry
|
||||
(
|
||||
vector& centre,
|
||||
vector& base,
|
||||
@ -59,23 +59,23 @@ void Foam::particle::tetGeometry
|
||||
}
|
||||
|
||||
|
||||
Foam::barycentricTensor Foam::particle::tetTransform() const
|
||||
Foam::barycentricTensor Foam::particle::stationaryTetTransform() const
|
||||
{
|
||||
vector centre, base, vertex1, vertex2;
|
||||
tetGeometry(centre, base, vertex1, vertex2);
|
||||
stationaryTetGeometry(centre, base, vertex1, vertex2);
|
||||
|
||||
return barycentricTensor(centre, base, vertex1, vertex2);
|
||||
}
|
||||
|
||||
|
||||
void Foam::particle::tetReverseTransform
|
||||
void Foam::particle::stationaryTetReverseTransform
|
||||
(
|
||||
vector& centre,
|
||||
scalar& detA,
|
||||
barycentricTensor& T
|
||||
) const
|
||||
{
|
||||
barycentricTensor A = tetTransform();
|
||||
barycentricTensor A = stationaryTetTransform();
|
||||
|
||||
const vector ab = A.b() - A.a();
|
||||
const vector ac = A.c() - A.a();
|
||||
@ -117,27 +117,11 @@ void Foam::particle::movingTetGeometry
|
||||
const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
|
||||
const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
|
||||
|
||||
scalar f0 = stepFraction_, f1 = fraction;
|
||||
|
||||
// Old and new points and cell centres are not sub-cycled. If we are sub-
|
||||
// cycling, then we have to account for the timestep change here by
|
||||
// modifying the fractions that we take of the old and new geometry.
|
||||
if (mesh_.time().subCycling())
|
||||
{
|
||||
const TimeState& tsNew = mesh_.time();
|
||||
const TimeState& tsOld = mesh_.time().prevTimeState();
|
||||
|
||||
const scalar tFrac =
|
||||
(
|
||||
(tsNew.value() - tsNew.deltaTValue())
|
||||
- (tsOld.value() - tsOld.deltaTValue())
|
||||
)/tsOld.deltaTValue();
|
||||
|
||||
const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue();
|
||||
|
||||
f0 = tFrac + dtFrac*f0;
|
||||
f1 = dtFrac*f1;
|
||||
}
|
||||
const Pair<scalar> s = stepFractionSpan();
|
||||
const scalar f0 = s[0] + stepFraction_*s[1], f1 = fraction*s[1];
|
||||
|
||||
centre[0] = ccOld + f0*(ccNew - ccOld);
|
||||
base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
|
||||
@ -497,8 +481,8 @@ void Foam::particle::locate
|
||||
|
||||
if (direction == nullptr)
|
||||
{
|
||||
const polyPatch& p = mesh_.boundaryMesh()[patch(facei_)];
|
||||
direction = &p.faceNormals()[patchFace(patch(facei_), facei_)];
|
||||
const polyPatch& p = mesh_.boundaryMesh()[patch()];
|
||||
direction = &p.faceNormals()[p.whichFace(facei_)];
|
||||
}
|
||||
|
||||
const vector n = *direction/mag(*direction);
|
||||
@ -646,14 +630,7 @@ Foam::scalar Foam::particle::trackToFace
|
||||
|
||||
while (true)
|
||||
{
|
||||
if (mesh_.moving())
|
||||
{
|
||||
f *= trackToMovingTri(f*displacement, f*fraction, tetTriI);
|
||||
}
|
||||
else
|
||||
{
|
||||
f *= trackToStationaryTri(f*displacement, f*fraction, tetTriI);
|
||||
}
|
||||
f *= trackToTri(f*displacement, f*fraction, tetTriI);
|
||||
|
||||
if (tetTriI == -1)
|
||||
{
|
||||
@ -688,19 +665,20 @@ Foam::scalar Foam::particle::trackToStationaryTri
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
|
||||
Info<< "Particle " << origId() << endl
|
||||
<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
|
||||
}
|
||||
|
||||
// Get the tet geometry
|
||||
vector centre;
|
||||
scalar detA;
|
||||
barycentricTensor T;
|
||||
tetReverseTransform(centre, detA, T);
|
||||
stationaryTetReverseTransform(centre, detA, T);
|
||||
|
||||
if (debug)
|
||||
{
|
||||
vector o, b, v1, v2;
|
||||
tetGeometry(o, b, v1, v2);
|
||||
stationaryTetGeometry(o, b, v1, v2);
|
||||
Info<< "Tet points o=" << o << ", b=" << b
|
||||
<< ", v1=" << v1 << ", v2=" << v2 << endl
|
||||
<< "Tet determinant = " << detA << endl
|
||||
@ -717,7 +695,7 @@ Foam::scalar Foam::particle::trackToStationaryTri
|
||||
|
||||
// Calculate the hit fraction
|
||||
label iH = -1;
|
||||
scalar muH = detA == 0 ? VGREAT : mag(1/detA);
|
||||
scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
|
||||
for (label i = 0; i < 4; ++ i)
|
||||
{
|
||||
if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
|
||||
@ -778,7 +756,7 @@ Foam::scalar Foam::particle::trackToStationaryTri
|
||||
// Set the proportion of the track that has been completed
|
||||
stepFraction_ += fraction*muH*detA;
|
||||
|
||||
return 1 - muH*detA;
|
||||
return iH != -1 ? 1 - muH*detA : 0;
|
||||
}
|
||||
|
||||
|
||||
@ -795,7 +773,8 @@ Foam::scalar Foam::particle::trackToMovingTri
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Info<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
|
||||
Info<< "Particle " << origId() << endl
|
||||
<< "Tracking from " << x0 << " to " << x0 + x1 << endl;
|
||||
}
|
||||
|
||||
// Get the tet geometry
|
||||
@ -823,7 +802,7 @@ Foam::scalar Foam::particle::trackToMovingTri
|
||||
|
||||
// Calculate the hit fraction
|
||||
label iH = -1;
|
||||
scalar muH = detA[0] == 0 ? VGREAT : mag(1/detA[0]);
|
||||
scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
|
||||
for (label i = 0; i < 4; ++ i)
|
||||
{
|
||||
const Roots<3> mu = hitEqn[i].roots();
|
||||
@ -899,7 +878,25 @@ Foam::scalar Foam::particle::trackToMovingTri
|
||||
<< "End global coordinates = " << position() << endl;
|
||||
}
|
||||
|
||||
return 1 - muH*detA[0];
|
||||
return iH != -1 ? 1 - muH*detA[0] : 0;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::particle::trackToTri
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
label& tetTriI
|
||||
)
|
||||
{
|
||||
if (mesh_.moving())
|
||||
{
|
||||
return trackToMovingTri(displacement, fraction, tetTriI);
|
||||
}
|
||||
else
|
||||
{
|
||||
return trackToStationaryTri(displacement, fraction, tetTriI);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -926,6 +923,47 @@ void Foam::particle::constrainToMeshCentre()
|
||||
}
|
||||
|
||||
|
||||
void Foam::particle::patchData(vector& n, vector& U) const
|
||||
{
|
||||
if (!onBoundaryFace())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Patch data was requested for a particle that isn't on a patch"
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if (mesh_.moving())
|
||||
{
|
||||
Pair<vector> centre, base, vertex1, vertex2;
|
||||
movingTetGeometry(1, centre, base, vertex1, vertex2);
|
||||
|
||||
n = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
|
||||
n /= mag(n);
|
||||
|
||||
// Interpolate the motion of the three face vertices to the current
|
||||
// coordinates
|
||||
U =
|
||||
coordinates_.b()*base[1]
|
||||
+ coordinates_.c()*vertex1[1]
|
||||
+ coordinates_.d()*vertex2[1];
|
||||
|
||||
// The movingTetGeometry method gives the motion as a displacement
|
||||
// across the time-step, so we divide by the time-step to get velocity
|
||||
U /= mesh_.time().deltaTValue();
|
||||
}
|
||||
else
|
||||
{
|
||||
vector centre, base, vertex1, vertex2;
|
||||
stationaryTetGeometry(centre, base, vertex1, vertex2);
|
||||
|
||||
n = triPointRef(base, vertex1, vertex2).normal();
|
||||
n /= mag(n);
|
||||
|
||||
U = Zero;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::particle::transformProperties(const tensor&)
|
||||
{}
|
||||
|
||||
@ -983,17 +1021,22 @@ void Foam::particle::correctAfterInteractionListReferral(const label celli)
|
||||
// position method is called. A referred particle should never be tracked,
|
||||
// so this approximate topology is good enough. By using the nearby cell we
|
||||
// minimise the error associated with the incorrect topology.
|
||||
coordinates_ = barycentric(1, 0, 0, 0);
|
||||
if (mesh_.moving())
|
||||
{
|
||||
NotImplemented;
|
||||
Pair<vector> centre;
|
||||
FixedList<scalar, 4> detA;
|
||||
FixedList<barycentricTensor, 3> T;
|
||||
movingTetReverseTransform(0, centre, detA, T);
|
||||
coordinates_ += (pos - centre[0]) & T[0]/detA[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
vector centre;
|
||||
scalar detA;
|
||||
barycentricTensor T;
|
||||
tetReverseTransform(centre, detA, T);
|
||||
coordinates_ = barycentric(1, 0, 0, 0) + ((pos - centre) & T/detA);
|
||||
stationaryTetReverseTransform(centre, detA, T);
|
||||
coordinates_ += (pos - centre) & T/detA;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -175,7 +175,7 @@ private:
|
||||
// Tetrahedra functions
|
||||
|
||||
//- Get the vertices of the current tet
|
||||
void tetGeometry
|
||||
void stationaryTetGeometry
|
||||
(
|
||||
vector& centre,
|
||||
vector& base,
|
||||
@ -188,7 +188,7 @@ private:
|
||||
// cartesian position in the global coordinate system. The
|
||||
// conversion is x = A & y, where x is the cartesian position, y is
|
||||
// the barycentric position and A is the transformation tensor.
|
||||
barycentricTensor tetTransform() const;
|
||||
barycentricTensor stationaryTetTransform() const;
|
||||
|
||||
//- Get the reverse transform associated with the current tet. The
|
||||
// conversion is detA*y = (x - centre) & T. The variables x, y and
|
||||
@ -197,7 +197,7 @@ private:
|
||||
// multiplied by its determinant, detA. This separation allows
|
||||
// the barycentric tracking algorithm to function on inverted or
|
||||
// degenerate tetrahedra.
|
||||
void tetReverseTransform
|
||||
void stationaryTetReverseTransform
|
||||
(
|
||||
vector& centre,
|
||||
scalar& detA,
|
||||
@ -454,10 +454,48 @@ public:
|
||||
//- Return current tet face particle is in
|
||||
inline label tetPt() const;
|
||||
|
||||
//- Return current face particle is on otherwise -1
|
||||
inline label face() const;
|
||||
|
||||
//- Return the fraction of time-step completed
|
||||
inline scalar stepFraction() const;
|
||||
|
||||
//- Return the fraction of time-step completed
|
||||
inline scalar& stepFraction();
|
||||
|
||||
//- Return the originating processor ID
|
||||
inline label origProc() const;
|
||||
|
||||
//- Return the originating processor ID
|
||||
inline label& origProc();
|
||||
|
||||
//- Return the particle ID on the originating processor
|
||||
inline label origId() const;
|
||||
|
||||
//- Return the particle ID on the originating processor
|
||||
inline label& origId();
|
||||
|
||||
|
||||
// Check
|
||||
|
||||
//- Return the step fraction change within the overall time-step.
|
||||
// Returns the start value and the change as a scalar pair. Always
|
||||
// return Pair<scalar>(0, 1), unless sub-cycling is in effect, in
|
||||
// which case the values will reflect the span of the sub-cycle
|
||||
// within the time-step.
|
||||
inline Pair<scalar> stepFractionSpan() const;
|
||||
|
||||
//- Return the current fraction within the timestep. This differs
|
||||
// from the stored step fraction due to sub-cycling.
|
||||
inline scalar currentTimeFraction() const;
|
||||
|
||||
//- Return the indices of the current tet that the
|
||||
// particle occupies.
|
||||
inline tetIndices currentTetIndices() const;
|
||||
|
||||
//- Return the current tet transformation tensor
|
||||
inline barycentricTensor currentTetTransform() const;
|
||||
|
||||
//- Return the normal of the tri on tetFacei_ for the
|
||||
// current tet.
|
||||
inline vector normal() const;
|
||||
@ -467,15 +505,6 @@ public:
|
||||
// on oldPoints
|
||||
inline vector oldNormal() const;
|
||||
|
||||
//- Return current face particle is on otherwise -1
|
||||
inline label face() const;
|
||||
|
||||
//- Return the particle current time
|
||||
inline scalar currentTime() const;
|
||||
|
||||
|
||||
// Check
|
||||
|
||||
//- Is the particle on a face?
|
||||
inline bool onFace() const;
|
||||
|
||||
@ -485,87 +514,80 @@ public:
|
||||
//- Is the particle on a boundary face?
|
||||
inline bool onBoundaryFace() const;
|
||||
|
||||
//- Which patch is particle on // <-- !!!
|
||||
inline label patch(const label facei) const;
|
||||
|
||||
//- Which face of this patch is this particle on // <-- !!!
|
||||
inline label patchFace(const label patchi, const label facei) const;
|
||||
//- Return the index of patch that the particle is on
|
||||
inline label patch() const;
|
||||
|
||||
//- Return current particle position
|
||||
inline vector position() const;
|
||||
|
||||
//- Return the fraction of time-step completed
|
||||
inline scalar& stepFraction();
|
||||
|
||||
//- Return the fraction of time-step completed
|
||||
inline scalar stepFraction() const;
|
||||
// Track
|
||||
|
||||
//- Return const access to the originating processor id
|
||||
inline label origProc() const;
|
||||
//- Track along the displacement for a given fraction of the overall
|
||||
// step. End when the track is complete, or when a boundary is hit.
|
||||
// On exit, stepFraction_ will have been incremented to the current
|
||||
// position, and facei_ will be set to the index of the boundary
|
||||
// face that was hit, or -1 if the track completed within a cell.
|
||||
// The proportion of the displacement still to be completed is
|
||||
// returned.
|
||||
scalar track
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction
|
||||
);
|
||||
|
||||
//- Return the originating processor id for manipulation
|
||||
inline label& origProc();
|
||||
//- As particle::track, but also stops on internal faces.
|
||||
scalar trackToFace
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction
|
||||
);
|
||||
|
||||
//- Return const access to the particle id on originating processor
|
||||
inline label origId() const;
|
||||
//- As particle::trackToFace, but also stops on tet triangles. On
|
||||
// exit, tetTriI is set to the index of the tet triangle that was
|
||||
// hit, or -1 if the end position was reached.
|
||||
scalar trackToTri
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
label& tetTriI
|
||||
);
|
||||
|
||||
//- Return the particle id on originating processor for manipulation
|
||||
inline label& origId();
|
||||
//- As particle::trackToTri, but for stationary meshes
|
||||
scalar trackToStationaryTri
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
label& tetTriI
|
||||
);
|
||||
|
||||
//- As particle::trackToTri, but for moving meshes
|
||||
scalar trackToMovingTri
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
label& tetTriI
|
||||
);
|
||||
|
||||
//- As non-templated particle::trackToFace, but with additional
|
||||
// boundary handling.
|
||||
template<class TrackData>
|
||||
void trackToFace
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
TrackData& td
|
||||
);
|
||||
|
||||
//- Set the constrained components of the particle position to the
|
||||
// mesh centre.
|
||||
void constrainToMeshCentre();
|
||||
|
||||
|
||||
// Track
|
||||
// Patch data
|
||||
|
||||
//- Track along the displacement for a given fraction of the overall
|
||||
// step. End when the track is complete, or when a boundary is hit.
|
||||
// On exit, stepFraction_ will have been incremented to the current
|
||||
// position, and facei_ will be set to the index of the boundary
|
||||
// face that was hit, or -1 if the track completed within a cell.
|
||||
// The proportion of the displacement still to be completed is
|
||||
// returned.
|
||||
scalar track
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction
|
||||
);
|
||||
|
||||
//- As particle::track, but also stops on internal faces.
|
||||
scalar trackToFace
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction
|
||||
);
|
||||
|
||||
//- As particle::trackToFace, but also stops on tet triangles. On
|
||||
// exit, tetTriI is set to the index of the tet triangle that was
|
||||
// hit, or -1 if the end position was reached.
|
||||
scalar trackToStationaryTri
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
label& tetTriI
|
||||
);
|
||||
|
||||
//- As particle::trackToTri, but for moving meshes
|
||||
scalar trackToMovingTri
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
label& tetTriI
|
||||
);
|
||||
|
||||
//- As non-templated particle::trackToFace, but with additional
|
||||
// boundary handling.
|
||||
template<class TrackData>
|
||||
void trackToFace
|
||||
(
|
||||
const vector& displacement,
|
||||
const scalar fraction,
|
||||
TrackData& td
|
||||
);
|
||||
|
||||
//- Set the constrained components of the particle position to the
|
||||
// mesh centre.
|
||||
void constrainToMeshCentre();
|
||||
//- Get the normal and velocity of the current patch location
|
||||
void patchData(vector& n, vector& U) const;
|
||||
|
||||
|
||||
// Transformations
|
||||
|
||||
@ -72,68 +72,19 @@ inline Foam::label Foam::particle::tetPt() const
|
||||
}
|
||||
|
||||
|
||||
inline Foam::tetIndices Foam::particle::currentTetIndices() const
|
||||
{
|
||||
return tetIndices(celli_, tetFacei_, tetPti_);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::vector Foam::particle::normal() const
|
||||
{
|
||||
return currentTetIndices().faceTri(mesh_).normal();
|
||||
}
|
||||
|
||||
|
||||
inline Foam::vector Foam::particle::oldNormal() const
|
||||
{
|
||||
return currentTetIndices().oldFaceTri(mesh_).normal();
|
||||
}
|
||||
|
||||
|
||||
inline Foam::label Foam::particle::face() const
|
||||
{
|
||||
return facei_;
|
||||
}
|
||||
|
||||
|
||||
inline bool Foam::particle::onFace() const
|
||||
{
|
||||
return facei_ >= 0;
|
||||
}
|
||||
|
||||
|
||||
inline bool Foam::particle::onInternalFace() const
|
||||
{
|
||||
return onFace() && mesh_.isInternalFace(facei_);
|
||||
}
|
||||
|
||||
|
||||
inline bool Foam::particle::onBoundaryFace() const
|
||||
{
|
||||
return onFace() && !mesh_.isInternalFace(facei_);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::vector Foam::particle::position() const
|
||||
{
|
||||
if (mesh_.moving())
|
||||
{
|
||||
return movingTetTransform(0)[0] & coordinates_;
|
||||
}
|
||||
else
|
||||
{
|
||||
return tetTransform() & coordinates_;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
inline Foam::scalar& Foam::particle::stepFraction()
|
||||
inline Foam::scalar Foam::particle::stepFraction() const
|
||||
{
|
||||
return stepFraction_;
|
||||
}
|
||||
|
||||
|
||||
inline Foam::scalar Foam::particle::stepFraction() const
|
||||
inline Foam::scalar& Foam::particle::stepFraction()
|
||||
{
|
||||
return stepFraction_;
|
||||
}
|
||||
@ -163,25 +114,96 @@ inline Foam::label& Foam::particle::origId()
|
||||
}
|
||||
|
||||
|
||||
inline Foam::scalar Foam::particle::currentTime() const
|
||||
inline Foam::Pair<Foam::scalar> Foam::particle::stepFractionSpan() const
|
||||
{
|
||||
return mesh_.time().value() + stepFraction_*mesh_.time().deltaTValue();
|
||||
if (mesh_.time().subCycling())
|
||||
{
|
||||
const TimeState& tsNew = mesh_.time();
|
||||
const TimeState& tsOld = mesh_.time().prevTimeState();
|
||||
|
||||
const scalar tFrac =
|
||||
(
|
||||
(tsNew.value() - tsNew.deltaTValue())
|
||||
- (tsOld.value() - tsOld.deltaTValue())
|
||||
)/tsOld.deltaTValue();
|
||||
|
||||
const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue();
|
||||
|
||||
return Pair<scalar>(tFrac, dtFrac);
|
||||
}
|
||||
else
|
||||
{
|
||||
return Pair<scalar>(0, 1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
inline Foam::label Foam::particle::patch(const label facei) const
|
||||
inline Foam::scalar Foam::particle::currentTimeFraction() const
|
||||
{
|
||||
return mesh_.boundaryMesh().whichPatch(facei);
|
||||
const Pair<scalar> s = stepFractionSpan();
|
||||
|
||||
return s[0] + stepFraction_*s[1];
|
||||
}
|
||||
|
||||
|
||||
inline Foam::label Foam::particle::patchFace
|
||||
(
|
||||
const label patchi,
|
||||
const label facei
|
||||
) const
|
||||
inline Foam::tetIndices Foam::particle::currentTetIndices() const
|
||||
{
|
||||
return mesh_.boundaryMesh()[patchi].whichFace(facei);
|
||||
return tetIndices(celli_, tetFacei_, tetPti_);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::barycentricTensor Foam::particle::currentTetTransform() const
|
||||
{
|
||||
if (mesh_.moving())
|
||||
{
|
||||
return movingTetTransform(0)[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
return stationaryTetTransform();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
inline Foam::vector Foam::particle::normal() const
|
||||
{
|
||||
return currentTetIndices().faceTri(mesh_).normal();
|
||||
}
|
||||
|
||||
|
||||
inline Foam::vector Foam::particle::oldNormal() const
|
||||
{
|
||||
return currentTetIndices().oldFaceTri(mesh_).normal();
|
||||
}
|
||||
|
||||
|
||||
inline bool Foam::particle::onFace() const
|
||||
{
|
||||
return facei_ >= 0;
|
||||
}
|
||||
|
||||
|
||||
inline bool Foam::particle::onInternalFace() const
|
||||
{
|
||||
return onFace() && mesh_.isInternalFace(facei_);
|
||||
}
|
||||
|
||||
|
||||
inline bool Foam::particle::onBoundaryFace() const
|
||||
{
|
||||
return onFace() && !mesh_.isInternalFace(facei_);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::label Foam::particle::patch() const
|
||||
{
|
||||
return mesh_.boundaryMesh().whichPatch(facei_);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::vector Foam::particle::position() const
|
||||
{
|
||||
return currentTetTransform() & coordinates_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -703,143 +703,39 @@ void Foam::KinematicCloud<CloudType>::patchData
|
||||
(
|
||||
const parcelType& p,
|
||||
const polyPatch& pp,
|
||||
const scalar trackFraction,
|
||||
const tetIndices& tetIs,
|
||||
vector& nw,
|
||||
vector& Up
|
||||
) const
|
||||
{
|
||||
label patchi = pp.index();
|
||||
label patchFacei = pp.whichFace(p.face());
|
||||
p.patchData(nw, Up);
|
||||
|
||||
vector n = tetIs.faceTri(mesh_).normal();
|
||||
n /= mag(n);
|
||||
|
||||
vector U = U_.boundaryField()[patchi][patchFacei];
|
||||
|
||||
// Unless the face is rotating, the required normal is n;
|
||||
nw = n;
|
||||
|
||||
if (!mesh_.moving())
|
||||
// If this is a wall patch, then there may be a non-zero tangential velocity
|
||||
// component; the lid velocity in a lid-driven cavity case, for example. We
|
||||
// want the particle to interact with this velocity, so we look it up in the
|
||||
// velocity field and use it to set the wall-tangential component.
|
||||
if (isA<wallPolyPatch>(pp))
|
||||
{
|
||||
// Only wall patches may have a non-zero wall velocity from
|
||||
// the velocity field when the mesh is not moving.
|
||||
const label patchi = pp.index();
|
||||
const label patchFacei = pp.whichFace(p.face());
|
||||
|
||||
if (isA<wallPolyPatch>(pp))
|
||||
// We only want to use the boundary condition value onlyif it is set
|
||||
// by the boundary condition. If the boundary values are extrapolated
|
||||
// (e.g., slip conditions) then they represent the motion of the fluid
|
||||
// just inside the domain rather than that of the wall itself.
|
||||
if (U_.boundaryField()[patchi].fixesValue())
|
||||
{
|
||||
Up = U;
|
||||
const vector Uw1 = U_.boundaryField()[patchi][patchFacei];
|
||||
const vector& Uw0 =
|
||||
U_.oldTime().boundaryField()[patchi][patchFacei];
|
||||
|
||||
const scalar f = p.currentTimeFraction();
|
||||
|
||||
const vector Uw = Uw0 + f*(Uw1 - Uw0);
|
||||
|
||||
const tensor nnw = nw*nw;
|
||||
|
||||
Up = (nnw & Up) + Uw - (nnw & Uw);
|
||||
}
|
||||
else
|
||||
{
|
||||
Up = Zero;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
vector U00 = U_.oldTime().boundaryField()[patchi][patchFacei];
|
||||
|
||||
vector n00 = tetIs.oldFaceTri(mesh_).normal();
|
||||
|
||||
// Difference in normal over timestep
|
||||
vector dn = Zero;
|
||||
|
||||
if (mag(n00) > SMALL)
|
||||
{
|
||||
// If the old normal is zero (for example in layer
|
||||
// addition) then use the current normal, meaning that the
|
||||
// motion can only be translational, and dn remains zero,
|
||||
// otherwise, calculate dn:
|
||||
|
||||
n00 /= mag(n00);
|
||||
|
||||
dn = n - n00;
|
||||
}
|
||||
|
||||
// Total fraction through the timestep of the motion,
|
||||
// including stepFraction before the current tracking step
|
||||
// and the current trackFraction
|
||||
// i.e.
|
||||
// let s = stepFraction, t = trackFraction
|
||||
// Motion of x in time:
|
||||
// |-----------------|---------|---------|
|
||||
// x00 x0 xi x
|
||||
//
|
||||
// where xi is the correct value of x at the required
|
||||
// tracking instant.
|
||||
//
|
||||
// x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
|
||||
//
|
||||
// i.e. the motion covered by previous tracking portions
|
||||
// within this timestep, and
|
||||
//
|
||||
// xi = x0 + t*(x - x0)
|
||||
// = t*x + (1 - t)*x0
|
||||
// = t*x + (1 - t)*(s*x + (1 - s)*x00)
|
||||
// = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
|
||||
//
|
||||
// let m = (s + t - s*t)
|
||||
//
|
||||
// xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
|
||||
//
|
||||
// In the same form as before.
|
||||
|
||||
scalar m =
|
||||
p.stepFraction()
|
||||
+ trackFraction
|
||||
- (p.stepFraction()*trackFraction);
|
||||
|
||||
// When the mesh is moving, the velocity field on wall patches
|
||||
// will contain the velocity associated with the motion of the
|
||||
// mesh, in which case it is interpolated in time using m.
|
||||
// For other patches the face velocity will need to be
|
||||
// reconstructed from the face centre motion.
|
||||
|
||||
const vector& Cf = mesh_.faceCentres()[p.face()];
|
||||
|
||||
vector Cf00 = mesh_.faces()[p.face()].centre(mesh_.oldPoints());
|
||||
|
||||
if (isA<wallPolyPatch>(pp))
|
||||
{
|
||||
Up = U00 + m*(U - U00);
|
||||
}
|
||||
else
|
||||
{
|
||||
Up = (Cf - Cf00)/mesh_.time().deltaTValue();
|
||||
}
|
||||
|
||||
if (mag(dn) > SMALL)
|
||||
{
|
||||
// Rotational motion, nw requires interpolation and a
|
||||
// rotational velocity around face centre correction to Up
|
||||
// is required.
|
||||
|
||||
nw = n00 + m*dn;
|
||||
|
||||
// Cf at tracking instant
|
||||
vector Cfi = Cf00 + m*(Cf - Cf00);
|
||||
|
||||
// Normal vector cross product
|
||||
vector omega = (n00 ^ n);
|
||||
|
||||
scalar magOmega = mag(omega);
|
||||
|
||||
// magOmega = sin(angle between unit normals)
|
||||
// Normalise omega vector by magOmega, then multiply by
|
||||
// angle/dt to give the correct angular velocity vector.
|
||||
omega *= Foam::asin(magOmega)/(magOmega*mesh_.time().deltaTValue());
|
||||
|
||||
// Project position onto face and calculate this position
|
||||
// relative to the face centre.
|
||||
vector facePos =
|
||||
p.position()
|
||||
- ((p.position() - Cfi) & nw)*nw
|
||||
- Cfi;
|
||||
|
||||
Up += (omega ^ facePos);
|
||||
}
|
||||
|
||||
// No further action is required if the motion is
|
||||
// translational only, nw and Up have already been set.
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -588,8 +588,6 @@ public:
|
||||
(
|
||||
const parcelType& p,
|
||||
const polyPatch& pp,
|
||||
const scalar trackFraction,
|
||||
const tetIndices& tetIs,
|
||||
vector& normal,
|
||||
vector& Up
|
||||
) const;
|
||||
|
||||
@ -296,17 +296,15 @@ bool Foam::KinematicParcel<ParcelType>::move
|
||||
}
|
||||
else
|
||||
{
|
||||
// Abandon the track, and move to the end of the sub-step. If the
|
||||
// the mesh is moving, this will implicitly move the parcel.
|
||||
if (mesh.moving())
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Tracking was abandoned on a moving mesh. Parcels may "
|
||||
<< "move unphysically as a result." << endl;
|
||||
}
|
||||
// At present the only thing that sets active_ to false is a stick
|
||||
// wall interaction. We want the position of the particle to remain
|
||||
// the same relative to the face that it is on. The local
|
||||
// coordinates therefore do not change. We still advance in time and
|
||||
// perform the relevant interactions with the fixed particle.
|
||||
p.stepFraction() += f;
|
||||
}
|
||||
|
||||
|
||||
const scalar dt = (p.stepFraction() - sfrac)*trackTime;
|
||||
|
||||
// Avoid problems with extremely small timesteps
|
||||
@ -325,7 +323,7 @@ bool Foam::KinematicParcel<ParcelType>::move
|
||||
|
||||
if (p.onBoundaryFace() && td.keepParticle)
|
||||
{
|
||||
if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
|
||||
if (isA<processorPolyPatch>(pbMesh[p.patch()]))
|
||||
{
|
||||
td.switchProcessor = true;
|
||||
}
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -177,7 +177,7 @@ void Foam::ParticleErosion<CloudType>::postPatch
|
||||
vector Up;
|
||||
|
||||
// patch-normal direction
|
||||
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
|
||||
this->owner().patchData(p, pp, nw, Up);
|
||||
|
||||
// particle velocity reletive to patch
|
||||
const vector& U = p.U() - Up;
|
||||
|
||||
@ -237,7 +237,7 @@ bool Foam::LocalInteraction<CloudType>::correct
|
||||
vector nw;
|
||||
vector Up;
|
||||
|
||||
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
|
||||
this->owner().patchData(p, pp, nw, Up);
|
||||
|
||||
// Calculate motion relative to patch velocity
|
||||
U -= Up;
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -171,7 +171,7 @@ bool Foam::MultiInteraction<CloudType>::correct
|
||||
if (p.face() != origFacei)
|
||||
{
|
||||
origFacei = p.face();
|
||||
patchi = p.patch(p.face());
|
||||
patchi = p.patch();
|
||||
|
||||
// Interaction model has moved particle off wall?
|
||||
if (patchi == -1)
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -74,7 +74,7 @@ bool Foam::Rebound<CloudType>::correct
|
||||
vector nw;
|
||||
vector Up;
|
||||
|
||||
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
|
||||
this->owner().patchData(p, pp, nw, Up);
|
||||
|
||||
// Calculate motion relative to patch velocity
|
||||
U -= Up;
|
||||
|
||||
@ -145,7 +145,7 @@ bool Foam::StandardWallInteraction<CloudType>::correct
|
||||
vector nw;
|
||||
vector Up;
|
||||
|
||||
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
|
||||
this->owner().patchData(p, pp, nw, Up);
|
||||
|
||||
// Calculate motion relative to patch velocity
|
||||
U -= Up;
|
||||
|
||||
@ -84,7 +84,7 @@ bool Foam::solidParticle::move
|
||||
|
||||
if (onBoundaryFace() && td.keepParticle)
|
||||
{
|
||||
if (isA<processorPolyPatch>(pbMesh[patch(face())]))
|
||||
if (isA<processorPolyPatch>(pbMesh[patch()]))
|
||||
{
|
||||
td.switchProcessor = true;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user