分享以及讨论 欧拉拉格朗日模型中 ParcelCollector

[星星星星晴]

下面是我修改的ParcelCollector.C 及.H 基于OF4.1
目的是通过在domain中建立一个面（同心圆或多边形），对碰到该面的parcel进行一些统计。
原code中只计算了mass 以及massflowrate。
我在做喷雾，在validation中需要做一个sauter mean diameter（D32）VS axial diameter 的图。
在试过其他方法并没有特别好的效果的情况下，试着修改了一下parcelcollector这个cloudfunction。

#include "ParticleCollector.H"
#include "Pstream.H"
#include "surfaceWriter.H"
#include "unitConversion.H"
#include "Random.H"
#include "triangle.H"
#include "cloud.H"
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
template<class CloudType>
void Foam::ParticleCollector<CloudType>::makeLogFile
(
    const faceList& faces,
    const Field<point>& points,
    const Field<scalar>& area
)
{
    // Create the output file if not already created
    if (log_)
    {
        if (debug)
        {
            Info<< "Creating output file" << endl;
        }
 
        if (Pstream::master())
        {
            // Create directory if does not exist
            mkDir(this->writeTimeDir());
 
            // Open new file at start up
            outputFilePtr_.reset
            (
                new OFstream(this->writeTimeDir()/(type() + ".dat"))
            );
 
            outputFilePtr_()
                << "# Source     : " << type() << nl
                << "# Bins       : " << faces.size() << nl
                << "# Total area : " << sum(area) << nl;
 
            outputFilePtr_()
                << "# Geometry   :" << nl
                << '#'
                << tab << "Bin"
                << tab << "(Centre_x Centre_y Centre_z)"
                << tab << "Area"
                << nl;
 
            forAll(faces, i)
            {
                outputFilePtr_()
                    << '#'
                    << tab << i
                    << tab << faces[i].centre(points)
                    << tab << area[i]
                    << nl;
            }
 
            outputFilePtr_()
                << '#' << nl
                << "# Output format:" << nl;
 
            forAll(faces, i)
            {
                word id = Foam::name(i);
                word binId = "bin_" + id;
 
                outputFilePtr_()
                    << '#'
                    << tab << "Time"
                    << tab << binId
                    << tab << "mass[" << id << "]"
                    << tab << "massFlowRate[" << id << "]"
                    << tab << "D3[" << id << "]"
                    << tab << "D2[" << id << "]"
                    << endl;
            }
        }
    }
}
 
 
template<class CloudType>
void Foam::ParticleCollector<CloudType>::initPolygons
(
    const List<Field<point>>& polygons
)
{
    mode_ = mtPolygon;
 
    label nPoints = 0;
    forAll(polygons, polyI)
    {
        label np = polygons[polyI].size();
        if (np < 3)
        {
            FatalIOErrorInFunction(this->coeffDict())
                << "polygons must consist of at least 3 points"
                << exit(FatalIOError);
        }
 
        nPoints += np;
    }
 
    label pointOffset = 0;
    points_.setSize(nPoints);
    faces_.setSize(polygons.size());
    faceTris_.setSize(polygons.size());
    area_.setSize(polygons.size());
    forAll(faces_, facei)
    {
        const Field<point>& polyPoints = polygons[facei];
        face f(identity(polyPoints.size()) + pointOffset);
        UIndirectList<point>(points_, f) = polyPoints;
        area_[facei] = f.mag(points_);
 
        DynamicList<face> tris;
        f.triangles(points_, tris);
        faceTris_[facei].transfer(tris);
 
        faces_[facei].transfer(f);
 
        pointOffset += polyPoints.size();
    }
}
 
 
template<class CloudType>
void Foam::ParticleCollector<CloudType>::initConcentricCircles()
{
    mode_ = mtConcentricCircle;
 
    vector origin(this->coeffDict().lookup("origin"));
 
    radius_ = this->coeffDict().lookup("radius");
    nSector_ = readLabel(this->coeffDict().lookup("nSector"));
 
    label nS = nSector_;
 
    vector refDir;
    if (nSector_ > 1)
    {
        refDir = this->coeffDict().lookup("refDir");
        refDir -= normal_[0]*(normal_[0] & refDir);
        refDir /= mag(refDir);
    }
    else
    {
        // set 4 quadrants for single sector cases
        nS = 4;
 
        vector tangent = Zero;
        scalar magTangent = 0.0;
 
        Random rnd(1234);
        while (magTangent < SMALL)
        {
            vector v = rnd.vector01();
 
            tangent = v - (v & normal_[0])*normal_[0];
            magTangent = mag(tangent);
        }
 
        refDir = tangent/magTangent;
    }
 
    scalar dTheta = 5.0;
    scalar dThetaSector = 360.0/scalar(nS);
    label intervalPerSector = max(1, ceil(dThetaSector/dTheta));
    dTheta = dThetaSector/scalar(intervalPerSector);
 
    label nPointPerSector = intervalPerSector + 1;
 
    label nPointPerRadius = nS*(nPointPerSector - 1);
    label nPoint = radius_.size()*nPointPerRadius;
    label nFace = radius_.size()*nS;
 
 
 
    // add origin
    nPoint++;
 
    points_.setSize(nPoint);
    faces_.setSize(nFace);
    area_.setSize(nFace);
 
    coordSys_ = cylindricalCS("coordSys", origin, normal_[0], refDir, false);
 
    List<label> ptIDs(identity(nPointPerRadius));
 
    points_[0] = origin;
 
    // points
    forAll(radius_, radI)
    {
        label pointOffset = radI*nPointPerRadius + 1;
 
        for (label i = 0; i < nPointPerRadius; i++)
        {
            label pI = i + pointOffset;
            point pCyl(radius_[radI], degToRad(i*dTheta), 0.0);
            points_[pI] = coordSys_.globalPosition(pCyl);
        }
    }
 
    // faces
    DynamicList<label> facePts(2*nPointPerSector);
    forAll(radius_, radI)
    {
        if (radI == 0)
        {
            for (label secI = 0; secI < nS; secI++)
            {
                facePts.clear();
 
                // append origin point
                facePts.append(0);
 
                for (label ptI = 0; ptI < nPointPerSector; ptI++)
                {
                    label i = ptI + secI*(nPointPerSector - 1);
                    label id = ptIDs.fcIndex(i - 1) + 1;
                    facePts.append(id);
                }
 
                label facei = secI + radI*nS;
 
                faces_[facei] = face(facePts);
                area_[facei] = faces_[facei].mag(points_);
            }
        }
        else
        {
            for (label secI = 0; secI < nS; secI++)
            {
                facePts.clear();
 
                label offset = (radI - 1)*nPointPerRadius + 1;
 
                for (label ptI = 0; ptI < nPointPerSector; ptI++)
                {
                    label i = ptI + secI*(nPointPerSector - 1);
                    label id = offset + ptIDs.fcIndex(i - 1);
                    facePts.append(id);
                }
                for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
                {
                    label i = ptI + secI*(nPointPerSector - 1);
                    label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
                    facePts.append(id);
                }
 
                label facei = secI + radI*nS;
 
                faces_[facei] = face(facePts);
                area_[facei] = faces_[facei].mag(points_);
            }
        }
    }
}
 
 
template<class CloudType>
void Foam::ParticleCollector<CloudType>::collectParcelPolygon
(
    const point& p1,
    const point& p2
) const
{
    label dummyNearType = -1;
    label dummyNearLabel = -1;
 
    forAll(faces_, facei)
    {
        const label facePoint0 = faces_[facei][0];
 
        const point& pf = points_[facePoint0];
 
        const scalar d1 = normal_[facei] & (p1 - pf);
        const scalar d2 = normal_[facei] & (p2 - pf);
 
        if (sign(d1) == sign(d2))
        {
            // did not cross polygon plane
            continue;
        }
 
        // intersection point
        const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
 
        const List<face>& tris = faceTris_[facei];
 
        // identify if point is within poly bounds
        forAll(tris, triI)
        {
            const face& tri = tris[triI];
            triPointRef t
            (
                points_[tri[0]],
                points_[tri[1]],
                points_[tri[2]]
            );
 
            if (t.classify(pIntersect, dummyNearType, dummyNearLabel))
            {
                hitFaceIDs_.append(facei);
            }
        }
    }
}
 
 
template<class CloudType>
void Foam::ParticleCollector<CloudType>::collectParcelConcentricCircles
(
    const point& p1,
    const point& p2
) const
{
    label secI = -1;
 
    const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
    const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
 
    if (sign(d1) == sign(d2))
    {
        // did not cross plane
        return;
    }
 
    // intersection point in cylindrical co-ordinate system
    const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
 
    scalar r = pCyl[0];
 
    if (r < radius_.last())
    {
        label radI = 0;
        while (r > radius_[radI])
        {
            radI++;
        }
 
        if (nSector_ == 1)
        {
            secI = 4*radI;
        }
        else
        {
            scalar theta = pCyl[1] + constant::mathematical::pi;
 
            secI =
                nSector_*radI
              + floor
                (
                    scalar(nSector_)*theta/constant::mathematical::twoPi
                );
        }
    }
 
    hitFaceIDs_.append(secI);
}
 
 
template<class CloudType>
void Foam::ParticleCollector<CloudType>::write()
{
    const fvMesh& mesh = this->owner().mesh();
    const Time& time = mesh.time();
    scalar timeNew = time.value();
    scalar timeElapsed = timeNew - timeOld_;
 
    totalTime_ += timeElapsed;
 
    const scalar alpha = (totalTime_ - timeElapsed)/totalTime_;
    const scalar beta = timeElapsed/totalTime_;
 
    forAll(faces_, facei)
    {
        massFlowRate_[facei] =
            alpha*massFlowRate_[facei] + beta*mass_[facei]/timeElapsed;
        massTotal_[facei] += mass_[facei];
        d3Total_[facei] += d3_[facei];
        d2Total_[facei] += d2_[facei];
    }
 
    const label proci = Pstream::myProcNo();
 
    Field<scalar> faceMassTotal(mass_.size(), 0.0);
    this->getModelProperty("massTotal", faceMassTotal);
 
    Field<scalar> faceMassFlowRate(massFlowRate_.size(), 0.0);
    this->getModelProperty("massFlowRate", faceMassFlowRate);
 
    Field<scalar> faceD3Total(d3_.size(), 0.0);
    this->getModelProperty("d3Total", faceD3Total);
 
    Field<scalar> faceD2Total(d2_.size(), 0.0);
    this->getModelProperty("d2Total", faceD2Total);
 
 
    scalar sumTotalMass = 0.0;
    scalar sumAverageMFR = 0.0;
    scalar sumTotalD3 = 0.0;
    scalar sumTotalD2 = 0.0;
    forAll(faces_, facei)
    {
 
        scalarList allProcMass(Pstream::nProcs());
        allProcMass[proci] = massTotal_[facei];
        Pstream::gatherList(allProcMass);
        faceMassTotal[facei] += sum(allProcMass);
 
        scalarList allProcMassFlowRate(Pstream::nProcs());
        allProcMassFlowRate[proci] = massFlowRate_[facei];
        Pstream::gatherList(allProcMassFlowRate);
        faceMassFlowRate[facei] += sum(allProcMassFlowRate);
 
        scalarList allProcD3(Pstream::nProcs());
        allProcD3[proci] = d3Total_[facei];
        Pstream::gatherList(allProcD3);
        faceD3Total[facei] += sum(allProcD3);
 
        scalarList allProcD2(Pstream::nProcs());
        allProcD2[proci] = d2Total_[facei];
        Pstream::gatherList(allProcD2);
        faceD2Total[facei] += sum(allProcD2);
 
        sumTotalMass += faceMassTotal[facei];
        sumAverageMFR += faceMassFlowRate[facei];
        sumTotalD3 += faceD3Total[facei];
        sumTotalD2 += faceD2Total[facei];
 
        //- Display values for each face
        if (outputFilePtr_.valid())
        {
            outputFilePtr_()
                << time.timeName()
                << tab << facei
                << tab << faceMassTotal[facei]
                << tab << faceMassFlowRate[facei]
                << tab << faceD3Total[facei]
                << tab << faceD2Total[facei]
                << endl;
        }
    }
 
    Info<< "    sum(total mass) = "             << sumTotalMass     << nl
        << "    sum(average mass flow rate) = " << sumAverageMFR    << nl
        << "    sum(total d3) = "               << sumTotalD3       << nl
        << "    sum(total d2) = "               << sumTotalD2       << nl
        << "    radius.size = "                 << radius_.size()   << nl
 
        << endl;
 
    if (surfaceFormat_ != "none")
    {
        if (Pstream::master())
        {
            autoPtr<surfaceWriter> writer(surfaceWriter::New(surfaceFormat_));
 
            writer->write
            (
                this->writeTimeDir(),
                "collector",
                points_,
                faces_,
                "massTotal",
                faceMassTotal,
                false
            );
 
            writer->write
            (
                this->writeTimeDir(),
                "collector",
                points_,
                faces_,
                "massFlowRate",
                faceMassFlowRate,
                false
            );
 
             writer->write
            (
                this->writeTimeDir(),
                "collector",
                points_,
                faces_,
                "d3Total",
                faceD3Total,
                false
            );
 
             writer->write
            (
                this->writeTimeDir(),
                "collector",
                points_,
                faces_,
                "d2Total",
                faceD2Total,
                false
            );
        }
    }
 
    if (resetOnWrite_)
    {
        Field<scalar> dummy(faceMassTotal.size(), 0.0);
        this->setModelProperty("massTotal", dummy);
        this->setModelProperty("massFlowRate", dummy);
        this->setModelProperty("d3Total", dummy);
        this->setModelProperty("d2Total", dummy);
 
        timeOld_ = timeNew;
        totalTime_ = 0.0;
    }
    else
    {
        this->setModelProperty("massTotal", faceMassTotal);
        this->setModelProperty("massFlowRate", faceMassFlowRate);
        this->setModelProperty("d3Total", faceD3Total);
        this->setModelProperty("d2Total", faceD2Total);
    }
 
    forAll(faces_, facei)
    {
        mass_[facei] = 0.0;
        massTotal_[facei] = 0.0;
        massFlowRate_[facei] = 0.0;
 
        d3_[facei] = 0.0;
        d3Total_[facei] = 0.0;
        d2_[facei] = 0.0;
        d2Total_[facei] = 0.0;
    }
 
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::ParticleCollector<CloudType>::ParticleCollector
(
    const dictionary& dict,
    CloudType& owner,
    const word& modelName
)
:
    CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
    mode_(mtUnknown),
    parcelType_(this->coeffDict().lookupOrDefault("parcelType", -1)),
    removeCollected_(this->coeffDict().lookup("removeCollected")),
    points_(),
    faces_(),
    faceTris_(),
    nSector_(0),
    radius_(),
    coordSys_(false),
    normal_(),
    negateParcelsOppositeNormal_
    (
        readBool(this->coeffDict().lookup("negateParcelsOppositeNormal"))
    ),
    surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
    resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
    totalTime_(0.0),
    mass_(),
    massTotal_(),
    massFlowRate_(),
    d3_(),
    d2_(),
    d3Total_(),
    d2Total_(),
    log_(this->coeffDict().lookup("log")),
    outputFilePtr_(),
    timeOld_(owner.mesh().time().value()),
    hitFaceIDs_()
{
    normal_ /= mag(normal_);
 
    word mode(this->coeffDict().lookup("mode"));
    if (mode == "polygon")
    {
        List<Field<point>> polygons(this->coeffDict().lookup("polygons"));
 
        initPolygons(polygons);
 
        vector n0(this->coeffDict().lookup("normal"));
        normal_ = vectorField(faces_.size(), n0);
    }
    else if (mode == "polygonWithNormal")
    {
        List<Tuple2<Field<point>, vector>> polygonAndNormal
        (
            this->coeffDict().lookup("polygons")
        );
 
        List<Field<point>> polygons(polygonAndNormal.size());
        normal_.setSize(polygonAndNormal.size());
 
        forAll(polygons, polyI)
        {
            polygons[polyI] = polygonAndNormal[polyI].first();
            normal_[polyI] = polygonAndNormal[polyI].second();
            normal_[polyI] /= mag(normal_[polyI]) + ROOTVSMALL;
        }
 
        initPolygons(polygons);
    }
    else if (mode == "concentricCircle")
    {
        vector n0(this->coeffDict().lookup("normal"));
        normal_ = vectorField(1, n0);
 
        initConcentricCircles();
    }
    else
    {
        FatalIOErrorInFunction(this->coeffDict())
            << "Unknown mode " << mode << ".  Available options are "
            << "polygon, polygonWithNormal and concentricCircle"
            << exit(FatalIOError);
    }
 
    mass_.setSize(faces_.size(), 0.0);
    massTotal_.setSize(faces_.size(), 0.0);
    massFlowRate_.setSize(faces_.size(), 0.0);
 
    d3_.setSize(faces_.size(), 0.0);
    d2_.setSize(faces_.size(), 0.0);
    d3Total_.setSize(faces_.size(), 0.0);
    d2Total_.setSize(faces_.size(), 0.0);
 
    makeLogFile(faces_, points_, area_);
}
 
 
template<class CloudType>
Foam::ParticleCollector<CloudType>::ParticleCollector
(
    const ParticleCollector<CloudType>& pc
)
:
    CloudFunctionObject<CloudType>(pc),
    mode_(pc.mode_),
    parcelType_(pc.parcelType_),
    removeCollected_(pc.removeCollected_),
    points_(pc.points_),
    faces_(pc.faces_),
    faceTris_(pc.faceTris_),
    nSector_(pc.nSector_),
    radius_(pc.radius_),
    coordSys_(pc.coordSys_),
    normal_(pc.normal_),
    negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
    surfaceFormat_(pc.surfaceFormat_),
    resetOnWrite_(pc.resetOnWrite_),
    totalTime_(pc.totalTime_),
    mass_(pc.mass_),
    massTotal_(pc.massTotal_),
    massFlowRate_(pc.massFlowRate_),
    d3_(pc.d3_),
    d2_(pc.d2_),
    d3Total_(pc.d3Total_),
    d2Total_(pc.d2Total_),
    log_(pc.log_),
    outputFilePtr_(),
    timeOld_(0.0),
    hitFaceIDs_()
{}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::ParticleCollector<CloudType>::~ParticleCollector()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
template<class CloudType>
void Foam::ParticleCollector<CloudType>::postMove
(
    parcelType& p,
    const label celli,
    const scalar dt,
    const point& position0,
    bool& keepParticle
)
{
    if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
    {
        return;
    }
 
    // slightly extend end position to avoid falling within tracking tolerances
    const point position1 = position0 + 1.0001*(p.position() - position0);
 
    hitFaceIDs_.clear();
 
    switch (mode_)
    {
        case mtPolygon:
        {
            collectParcelPolygon(position0, position1);
            break;
        }
        case mtConcentricCircle:
        {
            collectParcelConcentricCircles(position0, position1);
            break;
        }
        default:
        {
        }
    }
 
 
    forAll(hitFaceIDs_, i)
    {
        label facei = hitFaceIDs_[i];
        scalar m = p.nParticle() * p.mass();
        scalar d3p = p.nParticle() * p.d() * p.d()* p.d();
        scalar d2p = p.nParticle() * p.d() * p.d();
        if (negateParcelsOppositeNormal_)
        {
            vector Uhat = p.U();
            Uhat /= mag(Uhat) + ROOTVSMALL;
            if ((Uhat & normal_[facei]) < 0)
            {
                m *= -1.0;
            }
        }
 
        // add mass contribution
        mass_[facei] += m;
        d2_[facei] += d2p;
        d3_[facei] += d3p;
 
        if (nSector_ == 1)
        {
            mass_[facei + 1] += m;
            mass_[facei + 2] += m;
            mass_[facei + 3] += m;
 
            d2_[facei + 1] += d2p;
            d2_[facei + 2] += d2p;
            d2_[facei + 3] += d2p;
 
            d3_[facei + 1] += d3p;
            d3_[facei + 2] += d3p;
            d3_[facei + 3] += d3p;
 
        }
 
        if (removeCollected_)
        {
            keepParticle = false;
        }
 
    }
 
}
 

#ifndef ParticleCollector_H
#define ParticleCollector_H
 
#include "CloudFunctionObject.H"
#include "cylindricalCS.H"
#include "face.H"
#include "Switch.H"
#include "OFstream.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace Foam
{
 
/*---------------------------------------------------------------------------*\
                      Class ParticleCollector Declaration
\*---------------------------------------------------------------------------*/
 
template<class CloudType>
class ParticleCollector
:
    public CloudFunctionObject<CloudType>
{
public:
 
    enum modeType
    {
        mtPolygon,
        mtConcentricCircle,
        mtUnknown
    };
 
 
private:
 
    // Private Data
 
        // Typedefs
 
            //- Convenience typedef for parcel type
            typedef typename CloudType::parcelType parcelType;
 
        //- Collector mode type
        modeType mode_;
 
        //- Index of parcel types to collect (-1 by default = all particles)
        const label parcelType_;
 
        //- Flag to remove collected particles
        Switch removeCollected_;
 
        //- List of points
        Field<point> points_;
 
        //- List of faces
        List<face> faces_;
 
 
        // Polygon collector
 
            //- Triangulation of faces
            List<List<face>> faceTris_;
 
        // Concentric circles collector
 
            //- Number of sectors per circle
            label nSector_;
 
            //- List of radii
            List<scalar> radius_;
 
            //- Cylindrical co-ordinate system
            cylindricalCS coordSys_;
 
 
        //- Face areas
        Field<scalar> area_;
 
        //- Polygon normal vector per face
        Field<vector> normal_;
 
        //- Remove mass of parcel travelling in opposite direction to normal_
        bool negateParcelsOppositeNormal_;
 
        //- Surface output format
        const word surfaceFormat_;
 
        //- Flag to indicate whether data should be reset/cleared on writing
        Switch resetOnWrite_;
 
        //- Total time
        scalar totalTime_;
 
        //- Mass storage
        List<scalar> mass_;
 
        //- Mass total storage
        List<scalar> massTotal_;
 
        //- Mass flow rate storage
        List<scalar> massFlowRate_;
 
        List<scalar> d3_;
        List<scalar> d2_;
        List<scalar> d3Total_;
        List<scalar> d2Total_;
 
        //- Flag to indicate whether data should be written to file
        Switch log_;
 
        //- Output file pointer
        autoPtr<OFstream> outputFilePtr_;
 
        //- Last calculation time
        scalar timeOld_;
 
        //- Work list to store which faces are hit
        mutable DynamicList<label> hitFaceIDs_;
 
 
    // Private Member Functions
 
        //- Helper function to create log files
        void makeLogFile
        (
            const faceList& faces,
            const Field<point>& points,
            const Field<scalar>& area
        );
 
        //- Initialise polygon collectors
        void initPolygons(const List<Field<point>>& polygons);
 
        //- Initialise concentric circle collectors
        void initConcentricCircles();
 
        //- Collect parcels in polygon collectors
        void collectParcelPolygon
        (
            const point& p1,
            const point& p2
        ) const;
 
        //- Collect parcels in concentric circle collectors
        void collectParcelConcentricCircles
        (
            const point& p1,
            const point& p2
        ) const;
 
 
protected:
 
    // Protected Member Functions
 
        //- Write post-processing info
        void write();
 
 
public:
 
    //- Runtime type information
    TypeName("particleCollector");
 
 
    // Constructors
 
        //- Construct from dictionary
        ParticleCollector
        (
            const dictionary& dict,
            CloudType& owner,
            const word& modelName
        );
 
        //- Construct copy
        ParticleCollector(const ParticleCollector<CloudType>& pc);
 
        //- Construct and return a clone
        virtual autoPtr<CloudFunctionObject<CloudType>> clone() const
        {
            return autoPtr<CloudFunctionObject<CloudType>>
            (
                new ParticleCollector<CloudType>(*this)
            );
        }
 
 
    //- Destructor
    virtual ~ParticleCollector();
 
 
    // Member Functions
 
        // Access
 
            //- Return const access to the reset on write flag
            inline const Switch& resetOnWrite() const;
 
 
        // Evaluation
 
            //- Post-move hook
            virtual void postMove
            (
                parcelType& p,
                const label celli,
                const scalar dt,
                const point& position0,
                bool& keepParticle
            );
};
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace Foam
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#include "ParticleCollectorI.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#ifdef NoRepository
    #include "ParticleCollector.C"
#endif
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#endif
 
// ************************************************************************* //

[星星星星晴]

现在的情况是这样的，在相同位置有不同的particlecollector，但是得到的数据就完全不同。

 particleCollector1
    {
        type            particleCollector;
        mode            concentricCircle;
        origin          (0.0 0.0 0.975);
        radius          (0.19999);
        nSector         1;
        refDir          (1 0 0);
        normal          (0 0 1);
        negateParcelsOppositeNormal no;
        removeCollected no;
        surfaceFormat   vtk;
        resetOnWrite    no;
        log             yes;
   }

 particleCollector2
    {
        type            particleCollector;
        mode            concentricCircle;
        origin          (0.0 0.0 0.975);
        radius          (0.0005 0.0015 0.0025 0.0035 
0.0045 0.0055 0.0065 0.0075 
0.0085 0.0095 0.0105 0.0115 
0.0125 0.0135 0.0145 0.0155 
0.0165 0.0175 0.0185 0.0195 
0.0205 0.0215 0.0225 0.0235 
0.0245 0.0255 0.0265 0.0275 
0.0285 0.0295 0.0305 0.05 
0.075 0.1 0.15 0.19999);
        nSector         1;
        refDir          (1 0 0);
        normal          (0 0 1);
        negateParcelsOppositeNormal no;
        removeCollected no;
        surfaceFormat   vtk;
        resetOnWrite    no;
        log             yes;
    }

 particleCollector3
    {
        type            particleCollector;
        mode            concentricCircle;
        origin          (0.0 0.0 0.975);
        radius          (0.0005 0.0015 0.0025 0.0035 0.05 0.075 0.1 0.15 0.19999);
        nSector         1;
        refDir          (1 0 0);
        normal          (0 0 1);
        negateParcelsOppositeNormal no;
        removeCollected no;
        surfaceFormat   vtk;
        resetOnWrite    no;
        log             yes;
    }

输出结果如下。

    sum(total mass) = 1.75085e-05
    sum(average mass flow rate) = 0.0875421
    sum(total d3) = 3.34387e-08
    sum(total d2) = 0.000615982
    radius.size = 1
    nSector_ = 1

    sum(total mass) = 3.51273e-05
    sum(average mass flow rate) = 0.175637
    sum(total d3) = 6.70882e-08
    sum(total d2) = 0.00122567
    radius.size = 36
    nSector_ = 1

    sum(total mass) = 1.51653e-05
    sum(average mass flow rate) = 0.0758266
    sum(total d3) = 2.89636e-08
    sum(total d2) = 0.000605451
    radius.size = 9
    nSector_ = 1

如果说 mass和 massflow rate是相同的话，我可以理解，因为这个是原生的部分，我并没有修改什么，但是现在就输出成这个样子，我就很费解了。。希望大神能帮忙看一下

[chengan.wang]

@星星星星晴 相同的位置统计数据不一样，这个您找到原因了么？有点晕，不知道该怎么统计数据了。还有您的验证validation，sauter mean diameter（D32）VS axial diameter 的图，方便发出来做为一个标准验证么？

[星星星星晴]

@chengan-wang 我在别的帖子和你说过啦，我后来不用particlecollector这个原生的，就是输出文件，对文件的数据分析，python matlab，excel都行，剩下就是统计了

[chengan.wang]

@星星星星晴 您好，我有几个参数想请教您，
1.您的算例中radius分别设定1个，36个和9个数据，代表采集不同同心圆的数据？
2.nSector，自带算例设成10， 您设的是1。源代码里面Number of sectors per circle，这个怎么考虑
3.refDir代表啥意思呢？好像跟nSector数量还有关系
打扰了

[星星星星晴]

@chengan-wang 不好意思，我研究了我能用的部分我就再没看了，如果你研究明白了，可以分享一下。自己查code，做几个实验code

[chengan.wang]

@星星星星晴 感觉nSector不能设为1，最后输出的sum(total mass)，sum(average mass flow rate)是nSector>1的4倍左右。

    if (nSector_ > 1)
    {
        refDir = this->coeffDict().lookup("refDir");
        refDir -= normal_[0]*(normal_[0] & refDir);
        refDir /= mag(refDir);
    }
    else
    {
        // set 4 quadrants for single sector cases
        nS = 4;
 
        vector tangent = vector::zero;
        scalar magTangent = 0.0;
 
        Random rnd(1234);
        while (magTangent < SMALL)
        {
            vector v = rnd.vector01();
 
            tangent = v - (v & normal_[0])*normal_[0];
            magTangent = mag(tangent);
        }
 
        refDir = tangent/magTangent;
    }

可能跟源代码中nS = 4;有关系