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_debug(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  logger(ERROR, "Arguments of setPlasticParameters do not make sense");
140  }
141  setLoadingStiffness(loadingStiffness);
142  setUnloadingStiffnessMax(unloadingStiffnessMax);
143  setCohesionStiffness(cohesionStiffness);
144  setPenetrationDepthMax(penetrationDepthMax);
145 }
146 
151 {
152  return loadingStiffness_;
153 }
154 
159 {
160  return unloadingStiffnessMax_;
161 }
162 
167 {
168  return cohesionStiffness_;
169 }
170 
175 {
176  return penetrationDepthMax_;
177 }
178 
183 {
184  loadingStiffness_ = loadingStiffness;
185 }
186 
191 {
192  unloadingStiffnessMax_ = unloadingStiffnessMax;
193 }
194 
199 {
200  cohesionStiffness_ = cohesionStiffness;
201 }
202 
207 {
208  penetrationDepthMax_ = penetrationDepthMax;
209 }
210 
217 {
218  if (getConstantRestitution()) mass = 1;
219  return 0.02 * constants::pi /
220  std::sqrt(unloadingStiffnessMax_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass));
221 }
222 
228 {
229  logger.assert_debug(dissipation >= 0, "setDissipation(%): dissipation should be non-negative",dissipation);
230  dissipation_ = dissipation;
231 }
232 
237 {
238  setLoadingStiffness(new_.k);
239  setDissipation(new_.disp);
240 }
241 
246 {
247  return dissipation_;
248 }
249 
257 void
260 {
261  if (getConstantRestitution()) mass = 1;
262  if (eps == 0.0)
263  {
265  dissipation_ = std::sqrt(2.0 * mass * loadingStiffness_);
266  }
267  else
268  {
269  dissipation_ = -mass / tc * std::log(eps);
271  .5 * mass * (mathsFunc::square(constants::pi / tc) + mathsFunc::square(dissipation_ / mass));
272  }
274 }
275 
284 {
285  if (getConstantRestitution()) mass = 1;
286  loadingStiffness_ = stiffness;
287  setRestitutionCoefficient(eps, mass);
288 }
289 
298 {
299  if (getConstantRestitution()) mass = 1;
300  if (eps == 0.0) {
301  dissipation_ = std::sqrt(2.0 * mass * getLoadingStiffness());
302  } else {
303  const Mdouble logEps = log(eps);
304  dissipation_ = -std::sqrt(2.0 * mass * getLoadingStiffness()
305  / (constants::sqr_pi + mathsFunc::square(logEps))) * logEps;
306  }
307 }
308 
309 
315 {
316  if (getConstantRestitution()) mass = 1;
317 
318  logger.assert_debug(mass > 0, "getCollisionTime(%): mass should be positive",mass);
319 
320  Mdouble elasticContribution = loadingStiffness_ / (.5 * mass);
321  Mdouble elastoDissipativeContribution = elasticContribution - mathsFunc::square(dissipation_ / mass);
322 
323  logger.assert_debug(elastoDissipativeContribution > -1e-8 * elasticContribution,
324  "getCollisionTime(%): values for mass, stiffness and dissipation lead to an overdamped system: reduce dissipation",mass);
325 
326  return constants::pi / std::sqrt(elastoDissipativeContribution);
327 }
328 
330 {
331  if (getConstantRestitution()) mass = 1;
332 
333  return std::exp(-dissipation_ / mass * getCollisionTime(mass));
334 }
335 
336 Mdouble LinearPlasticViscoelasticNormalSpecies::computeBondNumberMax(Mdouble harmonicMeanRadius, Mdouble gravitationalAcceleration) const {
337  if (getConstantRestitution()) {
338  //harmonicMeanRadius unused
339  const Mdouble plasticOverlapMax = penetrationDepthMax_;
340  const Mdouble overlapMaxCohesion = plasticOverlapMax / (1 + cohesionStiffness_ / unloadingStiffnessMax_);
341  const Mdouble cohesionAccelerationMax = cohesionStiffness_ * overlapMaxCohesion;
342  return cohesionAccelerationMax / gravitationalAcceleration;
343  } else {
344  const Mdouble plasticOverlapMax = penetrationDepthMax_ * harmonicMeanRadius;
345  const Mdouble overlapMaxCohesion = plasticOverlapMax / (1 + cohesionStiffness_ / unloadingStiffnessMax_);
346  const Mdouble cohesionForceMax = cohesionStiffness_ * overlapMaxCohesion;
347  auto species = dynamic_cast<const ParticleSpecies*>(getBaseSpecies());
348  logger.assert_debug(species,"computeBondNumberMax: species needs to be a ParticleSpecies");
349  const Mdouble gravitationalForce = gravitationalAcceleration * species->getMassFromRadius(harmonicMeanRadius);
350  return cohesionForceMax / gravitationalForce;
351  }
352 }
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")
Definition of different loggers with certain modules. A user can define its own custom logger here...
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:106
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.