AxisymmetricHopper.h
Go to the documentation of this file.
1 //Copyright (c) 2013-2023, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include <string.h>
27 
28 #include "Chute.h"
30 #include "Math/ExtendedMath.h"
32 using namespace std;
33 
35 class AxisymmetricHopper : public Chute
36 {
37 public:
38 
40  {
41  setXMin(0);
42  setYMin(0);
43  setZMin(0);
44  setXMax(20);
45  setYMax(20);
46  setZMax(10);
47  setTimeMax(0);
48  FunnelMinRadius = -1; //radius at funnel outlet
49  FunnelMaxRadius = -1; //radius at upper funnel surface
50  FunnelHeight = -1; //height/thickness of funnel
51  FunnelInflowHeight = -1; //height/thickness of region where particles are included
52  num_created = 0;
53  max_failed = 100;
54  FunnelPointOnAxis.setZero();
55  }
56 
57  void setupInitialConditions() override
58  {
59  //do not write first time step, as there are zero particles in it
60  setLastSavedTimeStep(1);
61 
62  //Rudi's chute: rmin=H=12.5, rmax=34, dndt=95 500/s= 747 /sqrt(d/g)
63  //Rudi's large chute: rmin=14, H=25, rmax=55.5, dndt=215 000/s = 1682 /sqrt(d/g)
64  //tmax=3s=380 sqrt(d/g)
65  //real chute: dndt=700 000/s = 5477 /sqrt(d/g)
66 
67  if (Vec3D::getLength(FunnelPointOnAxis) == 0.0)
68  FunnelPointOnAxis = Vec3D(0.5 * (getXMax() + getXMin()), 0.5 * (getYMax() + getYMin()), 0);
69  Mdouble PossibleRadius = 0.5 * (getXMax() - getXMin());
70  if (FunnelMaxRadius == -1)
71  FunnelMaxRadius = PossibleRadius;
72  if (FunnelMinRadius == -1)
73  FunnelMinRadius = FunnelMaxRadius / 3.;
74  if (FunnelHeight == -1)
75  FunnelHeight = (FunnelMaxRadius - FunnelMinRadius) / mathsFunc::sin(45. * constants::pi / 180.); //60 deg (2./cos(60.*constants::pi/180.))
76  setZMax(getZMax() + FunnelHeight);
77  if (FunnelInflowHeight == -1)
78  FunnelInflowHeight = min(FunnelHeight / 4., 25. * getInflowParticleRadius());
79 
80  //create particle properties
81  Mdouble tc = .1;
82 
83  auto species = speciesHandler.copyAndAddObject(LinearViscoelasticFrictionSpecies());
84  species->setDensity(6.0 / constants::pi);
85  //species->setCollisionTimeAndRestitutionCoefficient(tc, 2.0 / 3.0, 0.5 * species->getMassFromRadius(0.5 * (getMinInflowParticleRadius() + getMaxInflowParticleRadius()))));
86  species->setCollisionTimeAndRestitutionCoefficient(tc, 2.0 / 3.0, species->getMassFromRadius(0.5 * (getMinInflowParticleRadius() + getMaxInflowParticleRadius())));
87  setTimeStep(tc / 10.);
88  setSaveCount(helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimeStep(50, getTimeMax(), getTimeStep()));
89  species->setSlidingStiffness(2. / 7. * species->getStiffness());
90  species->setSlidingDissipation(2. / 7. * species->getDissipation());
91  species->setSlidingFrictionCoefficient(0.4);
92  species->setRollingStiffness(2. / 7. * species->getStiffness());
93  species->setRollingDissipation(2. / 7. * species->getDissipation());
94  species->setRollingFrictionCoefficient(0.1);
95 
96  //set walls
98  w0.setSpecies(speciesHandler.getObject(0));
99  //for a prism wall, define points ABC (or more points)
100  //as the corners of the prism base (ordered in clockwise
101  //direction) and D as the 'periodic direction' of the prism
102  //(upwards: Vec3D::Cross(A-B,D) has to point into the wall)
103  Vec3D A(FunnelMinRadius, 0, getZMax() - FunnelHeight);
104  Vec3D B(FunnelMaxRadius, 0, getZMax());
105  Vec3D C(FunnelMaxRadius, 0, getZMax() - FunnelHeight);
106  Vec3D D(0, 1, 0); //Periodic direction of the prism
107  w0.addObject(Vec3D::cross(A - B, D), A);
108  w0.addObject(Vec3D::cross(B - C, D), B);
109  w0.addObject(Vec3D::cross(C - A, D), C);
110  w0.setPosition(FunnelPointOnAxis);
111  w0.setAxis(Vec3D(0, 0, 1));
112  wallHandler.copyAndAddObject(w0);
113 
114  write(std::cout, false);
115  }
116 
117  void actionsBeforeTimeStep() override
118  {
119  add_particles();
120  cleanChute();
121  }
122 
124  {
125  //the following formula yields polydispersed particle radii:
126  P0.setRadius(random.getRandomNumber(getMinInflowParticleRadius(), getMaxInflowParticleRadius()));
127  //P0.computeMass(speciesHandler);
128  Mdouble R = random.getRandomNumber(0, FunnelMaxRadius - P0.getRadius());
129  Mdouble A = random.getRandomNumber(0, 2 * constants::pi);
130  double xhat = R * mathsFunc::cos(A);
131  double zhat = random.getRandomNumber(getZMax() - FunnelInflowHeight + P0.getRadius(), getZMax() - P0.getRadius());
132  double c = mathsFunc::cos(getChuteAngle());
133  double s = mathsFunc::sin(getChuteAngle());
134  P0.setPosition(FunnelPointOnAxis + Vec3D(xhat * c - zhat * s, R * mathsFunc::sin(A), zhat * c + xhat * s));
135  P0.setVelocity(Vec3D(0, 0, 0));
136  P0.setSpecies(speciesHandler.getObject(0));
137  }
138 
140  {
141  unsigned int failed = 0;
142 
143  //try max_failed times to find new insertable particle
144  while (failed <= max_failed)
145  {
146  create_inflow_particle();
147  if (checkParticleForInteraction (P0))
148  {
149  particleHandler.copyAndAddObject(P0);
150  failed = 0;
151  num_created++;
152  }
153  else
154  {
155  failed++;
156  }
157  };
158  }
159 
160  virtual void cleanChute()
161  {
162  //clean outflow every 100 time steps
163  static int count = 0, maxcount = 100;
164  if (count > maxcount)
165  {
166  count = 0;
167  // delete all outflowing particles
168  for (unsigned int i = 0; i < particleHandler.getNumberOfObjects();)
169  {
170  if (particleHandler.getObject(i)->getPosition().Z < getZMin())
171  {
172  particleHandler.removeObject(i);
173  }
174  else
175  i++;
176  }
177  }
178  else
179  count++;
180  }
181 
182  bool readNextArgument(int& i, int argc, char *argv[]) override
183  {
184  if (!strcmp(argv[i], "-inflowHeight"))
185  {
186  setInflowHeight(atof(argv[i + 1]));
187  setZMax(atof(argv[i + 1]));
188  }
189  else if (!strcmp(argv[i], "-FunnelMinRadius"))
190  {
191  set_FunnelMinRadius(atof(argv[i + 1]));
192  }
193  else if (!strcmp(argv[i], "-FunnelMaxRadius"))
194  {
195  set_FunnelMaxRadius(atof(argv[i + 1]));
196  }
197  else if (!strcmp(argv[i], "-FunnelHeight"))
198  {
199  set_FunnelHeight(atof(argv[i + 1]));
200  }
201  else if (!strcmp(argv[i], "-FunnelInflowHeight"))
202  {
203  set_FunnelInflowHeight(atof(argv[i + 1]));
204  }
205  else
206  return Chute::readNextArgument(i, argc, argv); //if argv[i] is not found, check the commands in Mercury3D
207  return true; //returns true if argv[i] is found
208  }
209 
211  {
212  FunnelMinRadius = new_;
213  }
214 
216  {
217  FunnelMaxRadius = new_;
218  }
219 
221  {
222  FunnelHeight = new_;
223  }
224 
226  {
227  if (FunnelHeight < new_)
228  {
229  logger(ERROR, "FunnelInflowHeight=% < FunnelHeight=%", new_);
230  }
231  FunnelInflowHeight = new_;
232  }
233 
234 
235 protected:
242  unsigned int max_failed;
243  unsigned int num_created;
244 };
dominoes D
Definition: Domino.cpp:76
Species< LinearViscoelasticNormalSpecies, FrictionSpecies > LinearViscoelasticFrictionSpecies
Definition: LinearViscoelasticFrictionSpecies.h:34
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ ERROR
@ R
Definition: StatisticsVector.h:42
@ A
Definition: StatisticsVector.h:42
Definition: AxisymmetricHopper.h:36
void set_FunnelInflowHeight(Mdouble new_)
Definition: AxisymmetricHopper.h:225
void add_particles()
Definition: AxisymmetricHopper.h:139
void set_FunnelMinRadius(Mdouble new_)
Definition: AxisymmetricHopper.h:210
Mdouble FunnelMinRadius
Definition: AxisymmetricHopper.h:237
void set_FunnelHeight(Mdouble new_)
Definition: AxisymmetricHopper.h:220
void set_FunnelMaxRadius(Mdouble new_)
Definition: AxisymmetricHopper.h:215
Vec3D FunnelPointOnAxis
Definition: AxisymmetricHopper.h:236
bool readNextArgument(int &i, int argc, char *argv[]) override
Interprets the i^th command-line argument.
Definition: AxisymmetricHopper.h:182
void create_inflow_particle()
Definition: AxisymmetricHopper.h:123
Mdouble FunnelHeight
Definition: AxisymmetricHopper.h:239
Mdouble FunnelMaxRadius
Definition: AxisymmetricHopper.h:238
SphericalParticle P0
Definition: AxisymmetricHopper.h:241
unsigned int num_created
Definition: AxisymmetricHopper.h:243
Mdouble FunnelInflowHeight
Definition: AxisymmetricHopper.h:240
virtual void cleanChute()
Definition: AxisymmetricHopper.h:160
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: AxisymmetricHopper.h:57
AxisymmetricHopper()
Definition: AxisymmetricHopper.h:39
unsigned int max_failed
Definition: AxisymmetricHopper.h:242
void actionsBeforeTimeStep() override
A virtual function which allows to define operations to be executed before the new time step.
Definition: AxisymmetricHopper.h:117
Use AxisymmetricIntersectionOfWalls to Screw Screw::read Screw::read Screw::read define axisymmetric ...
Definition: AxisymmetricIntersectionOfWalls.h:126
void setAxis(Vec3D a)
Definition: AxisymmetricIntersectionOfWalls.cc:152
virtual void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Definition: BaseInteractable.h:239
Creates chutes with different bottoms. Inherits from Mercury3D (-> MercuryBase -> DPMBase).
Definition: Chute.h:65
bool readNextArgument(int &i, int argc, char *argv[]) override
This method can be used for reading object properties from a string.
Definition: Chute.cc:555
void addObject(Vec3D normal, Vec3D point)
Adds a wall to the set of infinite walls, given a normal vector pointing into the wall (i....
Definition: IntersectionOfWalls.cc:138
void setSpecies(const ParticleSpecies *species)
sets species of subwalls as well
Definition: IntersectionOfWalls.cc:72
A spherical particle is the most simple particle used in MercuryDPM.
Definition: SphericalParticle.h:37
Definition: Vector.h:51
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
Mdouble getLength() const
Calculates the length of this Vec3D: .
Definition: Vector.cc:320
const Mdouble pi
Definition: ExtendedMath.h:45
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
unsigned int getSaveCountFromNumberOfSavesAndTimeMaxAndTimeStep(unsigned int numberOfSaves, Mdouble timeMax, Mdouble timeStep)
Returns the correct saveCount if the total number of saves, the final time and the time step is known...
Definition: FormulaHelpers.cc:96
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44