MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LinearPlasticViscoelasticNormalSpecies.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, 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 
30 
31 class BaseParticle;
32 
33 class BaseInteractable;
34 
36  : BaseNormalForce()
37 {
38  loadingStiffness_ = 0.0;
40  cohesionStiffness_ = 0.0;
42  dissipation_ = 0.0;
44 #ifdef DEBUG_CONSTRUCTOR
45  std::cout<<"LinearPlasticViscoelasticNormalSpecies::LinearPlasticViscoelasticNormalSpecies() finished"<<std::endl;
46 #endif
47 }
48 
54  : BaseNormalForce(p)
55 {
62 #ifdef DEBUG_CONSTRUCTOR
63  std::cout<<"LinearPlasticViscoelasticNormalSpecies::LinearPlasticViscoelasticNormalSpecies(const LinearPlasticViscoelasticNormalSpecies &p) finished"<<std::endl;
64 #endif
65 }
66 
68 {
69 #ifdef DEBUG_DESTRUCTOR
70  std::cout<<"LinearPlasticViscoelasticNormalSpecies::~LinearPlasticViscoelasticNormalSpecies() finished"<<std::endl;
71 #endif
72 }
73 
78 {
79  os << " loadingStiffness " << loadingStiffness_;
80  os << " maxUnloadingStiffness " << unloadingStiffnessMax_;
81  if (doConstantUnloadingStiffness_) os << " doConstantUnloadingStiffness " << doConstantUnloadingStiffness_;
82  os << " cohesionStiffness " << cohesionStiffness_;
83  os << " maxPenetration " << penetrationDepthMax_;
84  os << " dissipation " << dissipation_;
85 }
86 
91 {
92  std::string dummy;
93  is >> dummy >> loadingStiffness_;
94  is >> dummy >> unloadingStiffnessMax_;
95  helpers::readOptionalVariable<bool>(is,"doConstantUnloadingStiffness",doConstantUnloadingStiffness_);
96  is >> dummy >> cohesionStiffness_;
97  is >> dummy >> penetrationDepthMax_;
98  is >> dummy >> dissipation_;
99 }
100 
105 {
106  return "LinearPlasticViscoelastic";
107 }
108 
116 {
122  logger.assert(S->doConstantUnloadingStiffness_==T->doConstantUnloadingStiffness_,"you cannot mix species where doConstantUnloadingStiffness is not the same");
124 }
125 
132 void
134  Mdouble cohesionStiffness, Mdouble penetrationDepthMax)
135 {
136  if (loadingStiffness <= 0 || unloadingStiffnessMax < loadingStiffness || cohesionStiffness < 0 ||
137  penetrationDepthMax < 0 || penetrationDepthMax > 1)
138  {
139  std::cerr << "Error: arguments of setPlasticParameters do not make sense" << std::endl;
140  exit(-1);
141  }
142  setLoadingStiffness(loadingStiffness);
143  setUnloadingStiffnessMax(unloadingStiffnessMax);
144  setCohesionStiffness(cohesionStiffness);
145  setPenetrationDepthMax(penetrationDepthMax);
146 }
147 
152 {
153  return loadingStiffness_;
154 }
155 
160 {
161  return unloadingStiffnessMax_;
162 }
163 
168 {
169  return cohesionStiffness_;
170 }
171 
176 {
177  return penetrationDepthMax_;
178 }
179 
184 {
185  loadingStiffness_ = loadingStiffness;
186 }
187 
192 {
193  unloadingStiffnessMax_ = unloadingStiffnessMax;
194 }
195 
200 {
201  cohesionStiffness_ = cohesionStiffness;
202 }
203 
208 {
209  penetrationDepthMax_ = penetrationDepthMax;
210 }
211 
218 {
219  if (getConstantRestitution()) mass = 1;
220  return 0.02 * constants::pi /
221  std::sqrt(unloadingStiffnessMax_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass));
222 }
223 
229 {
230  logger.assert(dissipation >= 0, "setDissipation(%): dissipation should be non-negative",dissipation);
231  dissipation_ = dissipation;
232 }
233 
238 {
239  setLoadingStiffness(new_.k);
240  setDissipation(new_.disp);
241 }
242 
247 {
248  return dissipation_;
249 }
250 
258 void
261 {
262  if (getConstantRestitution()) mass = 1;
263  if (eps == 0.0)
264  {
266  dissipation_ = std::sqrt(2.0 * mass * loadingStiffness_);
267  }
268  else
269  {
270  dissipation_ = -mass / tc * std::log(eps);
272  .5 * mass * (mathsFunc::square(constants::pi / tc) + mathsFunc::square(dissipation_ / mass));
273  }
275 }
276 
285 {
286  if (getConstantRestitution()) mass = 1;
287  loadingStiffness_ = stiffness;
288  setRestitutionCoefficient(eps, mass);
289 }
290 
299 {
300  if (getConstantRestitution()) mass = 1;
301  if (eps == 0.0) {
302  dissipation_ = std::sqrt(2.0 * mass * getLoadingStiffness());
303  } else {
304  const Mdouble logEps = log(eps);
305  dissipation_ = -std::sqrt(2.0 * mass * getLoadingStiffness()
306  / (constants::sqr_pi + mathsFunc::square(logEps))) * logEps;
307  }
308 }
309 
310 
316 {
317  if (getConstantRestitution()) mass = 1;
318 
319  logger.assert(mass > 0, "getCollisionTime(%): mass should be positive",mass);
320 
321  Mdouble elasticContribution = loadingStiffness_ / (.5 * mass);
322  Mdouble elastoDissipativeContribution = elasticContribution - mathsFunc::square(dissipation_ / mass);
323 
324  logger.assert(elastoDissipativeContribution > -1e-8 * elasticContribution,
325  "getCollisionTime(%): values for mass, stiffness and dissipation lead to an overdamped system: reduce dissipation",mass);
326 
327  return constants::pi / std::sqrt(elastoDissipativeContribution);
328 }
329 
331 {
332  if (getConstantRestitution()) mass = 1;
333 
334  return std::exp(-dissipation_ / mass * getCollisionTime(mass));
335 }
336 
337 Mdouble LinearPlasticViscoelasticNormalSpecies::computeBondNumberMax(Mdouble harmonicMeanRadius, Mdouble gravitationalAcceleration) const {
338  if (getConstantRestitution()) {
339  //harmonicMeanRadius unused
340  const Mdouble plasticOverlapMax = penetrationDepthMax_;
341  const Mdouble overlapMaxCohesion = plasticOverlapMax / (1 + cohesionStiffness_ / unloadingStiffnessMax_);
342  const Mdouble cohesionAccelerationMax = cohesionStiffness_ * overlapMaxCohesion;
343  return cohesionAccelerationMax / gravitationalAcceleration;
344  } else {
345  const Mdouble plasticOverlapMax = penetrationDepthMax_ * harmonicMeanRadius;
346  const Mdouble overlapMaxCohesion = plasticOverlapMax / (1 + cohesionStiffness_ / unloadingStiffnessMax_);
347  const Mdouble cohesionForceMax = cohesionStiffness_ * overlapMaxCohesion;
348  auto species = dynamic_cast<const ParticleSpecies*>(getBaseSpecies());
349  logger.assert(species,"computeBondNumberMax: species needs to be a ParticleSpecies");
350  const Mdouble gravitationalForce = gravitationalAcceleration * species->getMassFromRadius(harmonicMeanRadius);
351  return cohesionForceMax / gravitationalForce;
352  }
353 }
void setCohesionStiffness(Mdouble cohesionStiffness)
Sets the cohesive stiffness of the linear plastic-viscoelastic normal force.
bool getConstantRestitution() const
BaseSpecies * getBaseSpecies() const
Definition: BaseForce.h:35
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
return type specifically for fuctions returning k and disp at once
Definition: Helpers.h:45
Mdouble exp(Mdouble Exponent)
Definition: ExtendedMath.cc:84
Mdouble getUnloadingStiffnessMax() const
Returns the maximum unloading stiffness of the linear plastic-viscoelastic normal force...
MERCURY_DEPRECATED void setLoadingStiffnessAndDissipation(helpers::KAndDisp new_)
Allows the spring and dissipation constants to be changed simultaneously.
Mdouble getCollisionTime(Mdouble mass) const
Calculates collision time for two copies of a particle of given disp, k, mass.
Mdouble log(Mdouble Power)
Mdouble cohesionStiffness_
the adhesive spring constant (k^c) for plastic deformations
Mdouble penetrationDepthMax_
the depth (relative to the normalized radius) at which k_2^max is used (phi_f)
void setPenetrationDepthMax(Mdouble penetrationDepthMax)
Sets the maximum penetration depth of the linear plastic-viscoelastic normal force.
Mdouble getCohesionStiffness() const
Returns the cohesive stiffness of the linear plastic-viscoelastic normal force.
Mdouble computeTimeStep(Mdouble mass)
Returns the optimal time step to resolve a collision of two particles of a given mass.
Mdouble computeBondNumberMax(Mdouble harmonicMeanRadius, Mdouble gravitationalAcceleration) const
1) Computes the maximum plastic overlap delta_p* = phi*r 2) Computes the overlap at which the maximum...
Mdouble disp
Definition: Helpers.h:49
Mdouble unloadingStiffnessMax_
the maximum elastic constant (k_2^max) for plastic deformations
Mdouble getRestitutionCoefficient(Mdouble mass) const
Calculates restitution coefficient for two copies of given disp, k, mass.
const Mdouble pi
Definition: ExtendedMath.h:45
static Mdouble average(Mdouble a, Mdouble b)
Returns the harmonic mean of two variables.
Definition: BaseSpecies.cc:110
void setCollisionTimeAndRestitutionCoefficient(Mdouble tc, Mdouble eps, Mdouble mass)
Set k, disp such that is matches a given tc and eps for a collision of two different masses...
Mdouble getDissipation() const
Allows the normal dissipation to be accessed.
void setLoadingStiffness(Mdouble loadingStiffness)
Sets the loading stiffness of the linear plastic-viscoelastic normal force.
void setRestitutionCoefficient(double eps, Mdouble mass)
Sets disp to obtain a restitution coefficient eps for a collision of two particles of mass m...
void setDissipation(Mdouble dissipation)
Sets the linear dissipation coefficient of the linear plastic-viscoelastic normal force...
Mdouble loadingStiffness_
(normal) spring constant (k_1)
void setPlasticParameters(Mdouble loadingStiffness, Mdouble unloadingStiffnessMax, Mdouble cohesionStiffness, Mdouble penetrationDepthMax)
Sets all parameters of the linear plastic-viscoelastic normal force at once.
void setStiffnessAndRestitutionCoefficient(Mdouble k_, Mdouble eps, Mdouble mass)
Sets k, disp such that it matches a given tc and eps for a collision of two copies of P...
Defines the basic properties that a interactable object can have.
LinearPlasticViscoelasticNormalSpecies contains the parameters used to describe a plastic-cohesive no...
std::string getBaseName() const
Used in Species::getName to obtain a unique name for each Species.
Mdouble getPenetrationDepthMax() const
Returns the maximum penetration depth of the linear plastic-viscoelastic normal force.
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
void mix(LinearPlasticViscoelasticNormalSpecies *S, LinearPlasticViscoelasticNormalSpecies *T)
creates default values for mixed species
void setUnloadingStiffnessMax(Mdouble unloadingStiffnessMax)
Sets the maximum unloading stiffness of the linear plastic-viscoelastic normal force.
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
void write(std::ostream &os) const
Writes the species properties to an output stream.
void read(std::istream &is)
Reads the species properties from an input stream.
Mdouble getLoadingStiffness() const
Returns the loading stiffness of the linear plastic-viscoelastic normal force.