分享以及讨论 欧拉拉格朗日模型中 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是相同的话，我可以理解，因为这个是原生的部分，我并没有修改什么，但是现在就输出成这个样子，我就很费解了。。希望大神能帮忙看一下
