HorizontalMixer.h
Go to the documentation of this file.
1 #include "Mercury3D.h"
3 #include "Walls/RestrictedWall.h"
4 #include "Walls/InfiniteWall.h"
9 
10 #ifndef HorizontalMixer_h
11 #define HorizontalMixer_h
12 
13 class HorizontalMixer : public Mercury3D
14 {
15 private:
23  HorizontalScrew* screw = nullptr;
35  std::vector<unsigned> bladeMounts_;
40 
41 public:
45  bool haveOuterWalls = true;
46 
47 
48 public:
49 
55  {
59  }
60 
62  void setScrewWalls(const ParticleSpecies* s, Mdouble screwCenter,
63  Mdouble screwBaseHeight, Mdouble screwBaseRadius, Mdouble screwTopHeight,
64  Mdouble windingLength, Mdouble minR, Mdouble lowerR, Mdouble diffR, Mdouble thickness)
65  {
66  //Screws are dx=2200 apart, outerTopHeight height 1236, bottom height 110 (plus cap)
67  //winding width 411
68  Vec3D start = Vec3D(screwCenter, 0, screwBaseHeight);
69  Mdouble length = screwTopHeight-screwBaseHeight;
70  Mdouble numTurns = length/windingLength;
71  //create a screw that increases its length clockwise, just like (sin(t),cos(t), t).
72  // screw starts at height 0.11 until 1.197, distance 1.1 from the center in negative x-direction, making one turn each 0.411.
73  // rotation speed 1.0, thickness as small as possible, radius function is piecewise linear, decreasing from 0.6 to 0.4.
74  // screw starts at height 0.11 until 1.197, making one turn each 0.411.
76  (HorizontalScrew(start, length, minR, lowerR, diffR, numTurns, rotationSpeed/2.0/constants::pi, thickness, s));
77  //todo: change orientation of one screw
79  if (screw) screw->move_time(40);
80 
81  Mdouble sinA2Max = 0.25;
82  auto baseScrew = wallHandler.copyAndAddObject
83  (HorizontalBaseScrew(Vec3D(screwCenter,0,0),Vec3D(0,0,1),
84  {{Vec3D(0,0,-1),Vec3D(0,0,screwBaseHeight)},
85  {Vec3D(-1,0,0),Vec3D(screwBaseRadius,0,0)}},s,sinA2Max,timeMin));
86  baseScrew->setAngularVelocity(Vec3D(0,0,rotationSpeed));
87  //baseScrew->rotate(Vec3D(0,0,100));
88 
89  }
90 
92  virtual void setScrewCore(const ParticleSpecies* s, Mdouble screwCenter, Mdouble screwBaseHeight,
93  Mdouble coreTopHeight, Mdouble coreTopRadius,
94  Mdouble coreBottomHeight, Mdouble coreBottomRadius)
95  {
96  //the inner, thinner core cylinder in the screw
97  auto screwCore = wallHandler.copyAndAddObject
98  (AxisymmetricIntersectionOfWalls(Vec3D(screwCenter,0,0),Vec3D(0,0,1),
99  {{Vec3D(-.3,0,-1),Vec3D(coreTopRadius,0,coreTopHeight)},//slightly slant top, so particles roll off
100  {Vec3D(-1,0,0),Vec3D(coreTopRadius,0,0)}},s));
101  screwCore->setAngularVelocity(Vec3D(0,0,rotationSpeed));
102 
103  //the bottom, thicker core cylinder in the screw
104  auto screwBottom = wallHandler.copyAndAddObject
105  (AxisymmetricIntersectionOfWalls(Vec3D(screwCenter,0,0),Vec3D(0,0,1),
106  {{Vec3D(0,0,-1),Vec3D(0,0,coreBottomHeight)},
107  {Vec3D(-1,0,0),Vec3D(coreBottomRadius,0,0)}},s));
108  screwBottom->setAngularVelocity(Vec3D(0,0,rotationSpeed));
109  }
110 
112  virtual void setOuterWalls(const ParticleSpecies* s, Mdouble outerBaseRadius, Mdouble screwCenter, Mdouble outerTopCenter,
113  Mdouble outerTopRadius, Mdouble outerTopHeight)
114  {
115  //cylindrical wall
116  Vec3D normal = (Vec3D(outerTopCenter,0,outerTopHeight)-Vec3D(screwCenter,0,0));
117  normal.normalise();
118  Vec3D position = Vec3D(screwCenter,0,0);
119  Vec3D normalWall = Vec3D(outerTopHeight,0,outerBaseRadius-outerTopRadius);
120  normalWall.normalise();
121  Vec3D positionWall = Vec3D(outerBaseRadius,0,0);
122  AxisymmetricIntersectionOfWalls outerWall(position, normal, {{normalWall,positionWall}},s);
123  auto rightSide = wallHandler.copyAndAddObject(outerWall);
124 
125  //bottom plate
126  auto bottomPlate = wallHandler.copyAndAddObject
127  (InfiniteWall(Vec3D(0,0,-1),Vec3D(0,0,getZMin()),s));
128  }
129 
130  void setupInitialConditions() override
131  {
134 
135  // First, we define the screw
136  // bottom radius 0.974, outerTopHeight 1210, center-bottom x=1.21, center-outerTopHeight x=1.3385
137  Mdouble screwCenter = 1.21; //changed from 1.1
138  Mdouble screwBaseHeight = 0.11; //same
139  Mdouble screwTopHeight = 1.5; //changed from 1.155 //max. Height of the screw/somehow, screw is still felt on top of the core
140  Mdouble windingLength = 0.5; //changed from 0.411
141  Mdouble minR = 0.54;//change from 0.4 //radius in upper part
142  Mdouble lowerR = 1.195; //change from 0.6 ??//radius in upper part
143  Mdouble diffR = -0.6; //changed from -0.3 //radius in upper part
144  Mdouble thickness = 0.5*particleRadius; //needs change?
145 
146  //bottom screw core radius 270, outerTopHeight height 336
147  //outerTopHeight screw core radius ~170, outerTopHeight height 1236
148  Mdouble coreTopHeight = 1.5; //changed from 1236
149  Mdouble coreTopRadius = 0.17; //same
150  Mdouble coreBottomHeight = .491; //changed from 0.356
151  Mdouble coreBottomRadius = 0.27; //same
152  Mdouble screwBaseRadius = 1.0; //changed from 0.9
153 
154  //the cylindrical container:
155  Mdouble outerBaseRadius =1.208; //changed from 0.974
156  Mdouble outerTopCenter = screwCenter;// should be 1.3385, but for that I need to correct the implementation of orientation
157  Mdouble outerTopRadius =1.329; //changed from 1.21
158  Mdouble outerTopHeight = 1.995; // changed from 2.252
159 
160  setXMax(screwCenter+outerTopRadius);
161  setYMax(outerTopRadius);
162  setZMax(outerTopHeight);
163  setXMin(screwCenter-outerTopRadius);
164  setYMin(-getYMax());
165  setZMin(0.0);
166 
167  setScrewWalls(s, screwCenter, screwBaseHeight, screwBaseRadius, screwTopHeight, windingLength, minR, lowerR, diffR, thickness);
168  setScrewCore(s, screwCenter, screwBaseHeight, coreTopHeight, coreTopRadius, coreBottomHeight, coreBottomRadius);
169  if (haveOuterWalls) setOuterWalls(s, outerBaseRadius, screwCenter, outerTopCenter, outerTopRadius, outerTopHeight);
170 
171  //introduceSingleParticle(Vec3D(screwCenter,0.5*coreTopRadius,coreTopHeight));
172  //introduceParticlesAtWall();
174 
175  setGravity(Vec3D(0,0,-9.8));
176  }
177 
179  /* Simple run settings
180  * Nx*Ny*Nz particles are created evenly spaced between [xmin,xmax]*[ymin,ymax]*[zmin,zmax] and checked for contact with the screw
181  */
185  p0.setSpecies(s);
186  p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
188  p0.setPosition(p+Vec3D(0,0,particleRadius));
190  logger(INFO, "Inserted single particle ");
191  }
192 
194  /* Simple run settings
195  * Nx*Ny*Nz particles are created evenly spaced between [xmin,xmax]*[ymin,ymax]*[zmin,zmax] and checked for contact with the screw
196  */
200  p0.setSpecies(s);
201  p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
202  p0.setRadius(2.0*particleRadius);
204  p1.setSpecies(s);
205  p1.setVelocity(Vec3D(0.0, 0.0, 0.0));
207 
208  //number of particles that fit in domain
209  Mdouble distance;
210  Vec3D normal;
211  Vec3D p;
212  Mdouble minDistance;
213  unsigned counter = 0;
214  for (p.X = getXMin() + particleRadius; p.X < getXMax(); p.X += 2.0 * particleRadius)
215  for (p.Y = getYMin() + particleRadius; p.Y < getYMax(); p.Y += 2.0 * particleRadius)
216  for (p.Z = getZMin() + particleRadius; p.Z < fillHeight_; p.Z += 2.0 * particleRadius)
217  {
218  minDistance = p0.getRadius();
219  p0.setPosition(p);
220  for (auto w : wallHandler) {
221  //if touching the wall
222  if (w->getDistanceAndNormal(p0, distance, normal) && distance<minDistance)
223  {
224  minDistance=distance;
225  if (distance<0) break;
226  p1.setPosition(p0.getPosition()+(distance-1.0001*particleRadius)*normal);
227  }
228  }
229  if (minDistance<p0.getRadius() && minDistance>0)
230  {
232  counter++;
233  }
234  }
235  logger(INFO, "Inserted particles: %", counter);
236  }
237 
238  void introduceParticlesInDomain(Mdouble polydispersity=1) {
239  /* Simple run settings
240  * Nx*Ny*Nz particles are created evenly spaced between [xmin,xmax]*[ymin,ymax]*[zmin,zmax] and checked for contact with the screw
241  */
245  p0.setSpecies(s);
246  p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
248 
249  //CHANGED BY BERT need ~ 1.5 times cubic packing fraction to get HCP packing, assume particles are more HCP packed than cubic
250  Mdouble distance;
251  Vec3D normal;
252  Vec3D p;
253  Mdouble minDistance;
254  unsigned counter = 0;
255  for (p.X = getXMin() + particleRadius; p.X < getXMax(); p.X += 2.0 * particleRadius)
256  for (p.Y = getYMin() + particleRadius; p.Y < getYMax(); p.Y += 2.0 * particleRadius)
257  for (p.Z = getZMin() + particleRadius; p.Z < fillHeight_; p.Z += 2.0 * particleRadius) //Changed Cubic to HCP here
258  {
259  bool touch = false;
260  p0.setPosition(p);
261  for (auto w : wallHandler) {
262  //if touching the wall
263  if (w->getDistanceAndNormal(p0, distance, normal))
264  {
265  touch = true;
266  break;
267  }
268  }
269  if (!touch) {
272  counter++;
273  }
274  }
275  logger(INFO, "Inserted particles: %", counter);
276  }
277 
278  void printTime () const override
279  {
280  logger(INFO, "t=%3.6, tmax=%3.6, EneRatio=%3.6", getTime(), getTimeMax(),
282  }
283 
284  void actionsAfterTimeStep () override {
285  if (screw) screw->move_time(getTimeStep());
286  }
287 
288  void writeScript() {
289  logger(INFO,"Writing matlab script % to color the particles",getName() + ".m");
290  helpers::writeToFile(getName() + ".m", "cd " + helpers::getPath() + "\n"
291  "%% read in first file, to get the initial positions\n"
292  "f = fopen('" + getName() + "_7.vtu');\n"
293  "% header\n"
294  "line = textscan(f,'%s %s %s %s %s %s %s',1,'Delimiter','\\n');\n"
295  "% number of particles\n"
296  "N = textscan(line{5}{1}(24:end),'%d',1); N=N{1};\n"
297  "% positions\n"
298  "P = textscan(f,'%f %f %f',N);\n"
299  "%scatter(P{1},P{2})\n"
300  "fclose(f);\n"
301  "%% define a new speciesIndex, based on position, to color particles\n"
302  "index = 1000*P{1};\n"
303  "%% read in second file, a write out again with modified index\n"
304  "f = fopen('" + getName() + "_250.vtu');\n"
305  "g = fopen('Particle.vtu','w');\n"
306  "% header\n"
307  "line = textscan(f,'%s',3*N+15,'Delimiter','\\n');\n"
308  "for i=1:length(line{1}), fprintf(g,'%s\\n',line{1}{i}); end\n"
309  "% i/o indSpecies\n"
310  "textscan(f,'%f',N,'Delimiter','\\n');\n"
311  "fprintf(g,'%f\\n',index);\n"
312  "% footer\n"
313  "line = textscan(f,'%s','Delimiter','\\n');\n"
314  "for i=1:length(line{1}), fprintf(g,'%s\\n',line{1}{i}); end\n"
315  "fclose(f);\n"
316  "fclose(g);");
317  }
318 };
319 
320 #endif
@ MULTIPLE_FILES
each time-step will be written into/read from separate files numbered consecutively: name_....
@ NO_FILE
file will not be created/read
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ INFO
Use AxisymmetricIntersectionOfWalls to Screw Screw::read Screw::read Screw::read define axisymmetric ...
Definition: AxisymmetricIntersectionOfWalls.h:126
std::enable_if<!std::is_pointer< U >::value, U * >::type copyAndAddObject(const U &object)
Creates a copy of a Object and adds it to the BaseHandler.
Definition: BaseHandler.h:379
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
void setAngularVelocity(const Vec3D &angularVelocity)
set the angular velocity of the BaseInteractble.
Definition: BaseInteractable.cc:360
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
Definition: BaseInteractable.cc:350
virtual void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Definition: BaseInteractable.h:239
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Definition: BaseInteractable.h:218
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:348
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species)
Definition: BaseParticle.cc:553
void setSpecies(const ParticleSpecies *species)
Definition: BaseParticle.cc:818
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin.
Definition: DPMBase.h:619
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax.
Definition: DPMBase.h:626
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1427
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
Definition: DPMBase.h:1483
void setYMin(Mdouble newYMin)
Sets the value of YMin, the lower bound of the problem domain in the y-direction.
Definition: DPMBase.cc:1034
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin.
Definition: DPMBase.h:632
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1250
const std::string & getName() const
Returns the name of the file. Does not allow to change it though.
Definition: DPMBase.cc:399
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:808
Mdouble getKineticEnergy() const
Returns the global kinetic energy stored in the system.
Definition: DPMBase.cc:1544
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1447
void setYMax(Mdouble newYMax)
Sets the value of YMax, the upper bound of the problem domain in the y-direction.
Definition: DPMBase.cc:1191
void setZMin(Mdouble newZMin)
Sets the value of ZMin, the lower bound of the problem domain in the z-direction.
Definition: DPMBase.cc:1058
void setXMax(Mdouble newXMax)
Sets the value of XMax, the upper bound of the problem domain in the x-direction.
Definition: DPMBase.cc:1165
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created.
Definition: DPMBase.h:1437
void setZMax(Mdouble newZMax)
Sets the value of ZMax, the upper bound of the problem domain in the z-direction.
Definition: DPMBase.cc:1217
void setParticlesWriteVTK(bool writeParticlesVTK)
Sets whether particles are written in a VTK file.
Definition: DPMBase.cc:942
RNG random
This is a random generator, often used for setting up the initial conditions etc.....
Definition: DPMBase.h:1432
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax.
Definition: DPMBase.h:638
void setXMin(Mdouble newXMin)
Sets the value of XMin, the lower bound of the problem domain in the x-direction.
Definition: DPMBase.cc:1010
Mdouble getTimeMax() const
Returns the maximum simulation duration.
Definition: DPMBase.cc:888
void setGravity(Vec3D newGravity)
Sets a new value for the gravitational acceleration.
Definition: DPMBase.cc:1383
Mdouble getElasticEnergy() const
Returns the global elastic energy within the system.
Definition: DPMBase.cc:1530
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin.
Definition: DPMBase.h:644
void setFileType(FileType fileType)
Sets the type of file needed to write into or read from. File::fileType_.
Definition: File.cc:215
A HorizontalBaseScrew is a copy of AxisymmetricIntersectionOfWalls, with an additional,...
Definition: HorizontalBaseScrew.h:39
Definition: HorizontalMixer.h:14
HorizontalMixer(Mdouble particleRadius, Mdouble rotationSpeed, Mdouble timeMin, Mdouble fillHeight)
Constructor, turns off fstat output by default.
Definition: HorizontalMixer.h:53
void introduceParticlesAtWall()
Definition: HorizontalMixer.h:193
void printTime() const override
Displays the current simulation time and the maximum simulation duration.
Definition: HorizontalMixer.h:278
void actionsAfterTimeStep() override
A virtual function which allows to define operations to be executed after time step.
Definition: HorizontalMixer.h:284
void introduceParticlesInDomain(Mdouble polydispersity=1)
Definition: HorizontalMixer.h:238
virtual void setOuterWalls(const ParticleSpecies *s, Mdouble outerBaseRadius, Mdouble screwCenter, Mdouble outerTopCenter, Mdouble outerTopRadius, Mdouble outerTopHeight)
sets other walls that define the outer boundary
Definition: HorizontalMixer.h:112
virtual void setScrewCore(const ParticleSpecies *s, Mdouble screwCenter, Mdouble screwBaseHeight, Mdouble coreTopHeight, Mdouble coreTopRadius, Mdouble coreBottomHeight, Mdouble coreBottomRadius)
sets four walls, leftScrewCore, rightScrewCore, leftScrewBottom, rightScrewBottom
Definition: HorizontalMixer.h:92
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: HorizontalMixer.h:130
void introduceSingleParticle(Vec3D p)
Definition: HorizontalMixer.h:178
std::vector< unsigned > bladeMounts_
The index number of all mounted blades (the blades mounts are numbered 0-11, with the i-th blade moun...
Definition: HorizontalMixer.h:35
Mdouble fillHeight_
Definition: HorizontalMixer.h:39
void setScrewWalls(const ParticleSpecies *s, Mdouble screwCenter, Mdouble screwBaseHeight, Mdouble screwBaseRadius, Mdouble screwTopHeight, Mdouble windingLength, Mdouble minR, Mdouble lowerR, Mdouble diffR, Mdouble thickness)
sets four walls, leftScrew, rightScrew, leftBaseScrew, rightBaseScrew
Definition: HorizontalMixer.h:62
Mdouble timeMin
The time where the screw begins to spin.
Definition: HorizontalMixer.h:31
Mdouble rotationSpeed
The rotation speed of the screw.
Definition: HorizontalMixer.h:27
void writeScript()
Definition: HorizontalMixer.h:288
Mdouble particleRadius
The mean radius of the particles in the feeder.
Definition: HorizontalMixer.h:19
bool haveOuterWalls
Definition: HorizontalMixer.h:45
HorizontalScrew * screw
Pointer to the right screw.
Definition: HorizontalMixer.h:23
This function defines an Archimedes' screw in the z-direction from a (constant) starting point,...
Definition: HorizontalScrew.h:39
void move_time(Mdouble dt)
Rotate the HorizontalScrew for a period dt, so that the offset_ changes with omega_*dt.
Definition: HorizontalScrew.cc:274
A infinite wall fills the half-space {point: (position_-point)*normal_<=0}.
Definition: InfiniteWall.h:48
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:37
void clear() override
Empties the whole ParticleHandler by removing all BaseParticle.
Definition: ParticleHandler.cc:977
Definition: ParticleSpecies.h:37
Mdouble getRandomNumber()
This is a random generating routine can be used for initial positions.
Definition: RNG.cc:143
A spherical particle is the most simple particle used in MercuryDPM.
Definition: SphericalParticle.h:37
Definition: Vector.h:51
Mdouble Y
Definition: Vector.h:66
Mdouble Z
Definition: Vector.h:66
Mdouble X
the vector components
Definition: Vector.h:66
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
void setWriteVTK(FileType)
Sets whether walls are written into a VTK file.
Definition: WallHandler.cc:467
const Mdouble pi
Definition: ExtendedMath.h:45
std::string getPath()
Definition: FileIOHelpers.cc:224
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: FileIOHelpers.cc:58