CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

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

    OpenFOAM
    2
    7
    529
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • 星
      星星星星晴 最后由 编辑

      下面是我修改的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
      
      // ************************************************************************* //
      

      m.sui20@foxmail.com

      1 条回复 最后回复 回复 引用
      • 星
        星星星星晴 最后由 编辑

        现在的情况是这样的,在相同位置有不同的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是相同的话,我可以理解,因为这个是原生的部分,我并没有修改什么,但是现在就输出成这个样子,我就很费解了。。希望大神能帮忙看一下

        m.sui20@foxmail.com

        chengan.wang 1 条回复 最后回复 回复 引用
        • chengan.wang
          chengan.wang @星星星星晴 最后由 编辑

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

          星 1 条回复 最后回复 回复 引用
          • 星
            星星星星晴 @chengan.wang 最后由 编辑

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

            m.sui20@foxmail.com

            chengan.wang 1 条回复 最后回复 回复 引用
            • chengan.wang
              chengan.wang @星星星星晴 最后由 编辑

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

              星 1 条回复 最后回复 回复 引用
              • 星
                星星星星晴 @chengan.wang 最后由 编辑

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

                m.sui20@foxmail.com

                chengan.wang 1 条回复 最后回复 回复 引用
                • chengan.wang
                  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;有关系

                  1 条回复 最后回复 回复 引用
                  • First post
                    Last post