MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HertzianViscoelasticNormalSpecies.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 
28 #include <cmath>
31 #include "Logger.h"
32 #include "Species/BaseSpecies.h"
33 
34 class BaseParticle;
35 
36 class BaseInteractable;
37 
39  : BaseNormalForce()
40 {
41  elasticModulus_ = 0;
42  dissipation_ = 0;
43 #ifdef DEBUG_CONSTRUCTOR
44  std::cout<<"HertzianViscoelasticNormalSpecies::HertzianViscoelasticNormalSpecies() finished"<<std::endl;
45 #endif
46 }
47 
52  : BaseNormalForce(p)
53 {
56 #ifdef DEBUG_CONSTRUCTOR
57  std::cout<<"HertzianViscoelasticNormalSpecies::HertzianViscoelasticNormalSpecies(const HertzianViscoelasticNormalSpecies &p) finished"<<std::endl;
58 #endif
59 }
60 
62 {
63 #ifdef DEBUG_DESTRUCTOR
64  std::cout<<"HertzianViscoelasticNormalSpecies::~HertzianViscoelasticNormalSpecies() finished"<<std::endl;
65 #endif
66 }
67 
71 void HertzianViscoelasticNormalSpecies::write(std::ostream& os) const
72 {
73  os << " elasticModulus " << elasticModulus_
74  << " dissipation " << dissipation_;
75 }
76 
81 {
82  std::string dummy;
83  is >> dummy >> elasticModulus_
84  >> dummy >> dissipation_;
85 }
86 
91 {
92  return "HertzianViscoelastic";
93 }
94 
97 {
98  if (elasticModulus >= 0)
99  elasticModulus_ = elasticModulus;
100  else
101  {
102  logger(ERROR, "setEffectiveElasticModulus(%) argument has to be non-negative!", elasticModulus);
103  }
104 }
105 
108 {
109  // Here is a very nice paper describing contact modelling
110  // http://people.ds.cam.ac.uk/jae1001/CUS/research/pfizer/Antypov_Elliott_EPL_2011.pdf
111  // see also: https://answers.launchpad.net/yade/+question/235934
112  if (elasticModulus < 0.0)
113  logger(ERROR,
114  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] elasticModulus % should be nonnegative",
115  elasticModulus);
116 
117  else if (rest < 0.0 || rest > 1.0)
118  logger(ERROR,
119  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] rest % should be between 0 and 1 (inclusive)",
120  rest);
121 
122  else
123  {
124  elasticModulus_ = elasticModulus;
125  if (rest > 0.0)
126  {
127  Mdouble logRestSquared = log(rest) * log(rest);
128  dissipation_ = sqrt(5.0 * logRestSquared / (logRestSquared + constants::sqr_pi));
129  }
130  else
131  dissipation_ = sqrt(5.0);
132 
133  //logger(INFO, "Effective elastic modulus %", elasticModulus_);
134  //logger(INFO, "Set dissipation % to match restitution coefficient %", dissipation_, rest);
135  }
136 }
137 
140 {
141  if (elasticModulus < 0.0)
142  logger(ERROR,
143  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] elasticModulus % should be nonnegative",
144  elasticModulus);
145 
146  else if (poissonRatio < 0.0 || poissonRatio > 1.0)
147  logger(ERROR,
148  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] poissonRatio % should be between 0 and 1 (inclusive)",
149  poissonRatio);
150  else
151  {
152  elasticModulus_ = elasticModulus;
153  auto mindlin = dynamic_cast<MindlinSpecies*>(getBaseSpecies());
154  logger.assert_debug(mindlin, "Please define HertzianViscoelasticMindlinSpecies to use this setter");
155  mindlin->setEffectiveShearModulus(elasticModulus / 2 * (1 + poissonRatio));
156  }
157 }
158 
161 {
162  if (elasticModulus < 0.0)
163  logger(ERROR,
164  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndEffectiveShearModulus] elasticModulus % should be nonnegative",
165  elasticModulus);
166 
167  else if (shearModulus < 0.0)
168  logger(ERROR,
169  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndEffectiveShearModulus] shearModulus % should be nonnegative",
170  shearModulus);
171  else
172  {
173  elasticModulus_ = elasticModulus;
174  auto mindlin = dynamic_cast<MindlinSpecies*>(getBaseSpecies()->getFrictionForce());
175  logger.assert_debug(mindlin, "Please define HertzianViscoelasticMindlinSpecies to use this setter");
176  mindlin->setEffectiveShearModulus(shearModulus);
177  }
178 }
179 
182 {
183  return elasticModulus_;
184 }
185 
188 {
189  if (dissipation >= 0)
190  {
191  dissipation_ = dissipation;
192  }
193  else
194  {
195  logger(ERROR, "setDissipation(%) argument has to be non-negative!", dissipation);
196  }
197 }
198 
201 {
202  return dissipation_;
203 }
204 
206 //Mdouble HertzianViscoelasticNormalSpecies::getCollisionTime(Mdouble mass)
207 //{
208 // if (mass <= 0)
209 // {
210 // std::cerr << "Error in getCollisionTime(Mdouble mass) mass is not set or has an unexpected value, (getCollisionTime(" << mass << "))" << std::endl;
211 // exit(-1);
212 // }
213 // if (elasticModulus_ <= 0)
214 // {
215 // std::cerr << "Error in getCollisionTime(Mdouble mass) stiffness is not set or has an unexpected value, (getCollisionTime(" << mass << "), with stiffness=" << elasticModulus_ << ")" << std::endl;
216 // exit(-1);
217 // }
218 // if (dissipation_ < 0)
219 // {
220 // std::cerr << "Error in getCollisionTime(Mdouble mass) dissipation is not set or has an unexpected value, (getCollisionTime(" << mass << "), with dissipation=" << dissipation_ << ")" << std::endl;
221 // exit(-1);
222 // }
223 // Mdouble tosqrt = elasticModulus_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass);
224 // if (tosqrt <= 0)
225 // {
226 // std::cerr << "Error in getCollisionTime(Mdouble mass) values for mass, stiffness and dissipation would lead to an overdamped system, (getCollisionTime(" << mass << "), with stiffness=" << elasticModulus_ << " and dissipation=" << dissipation_ << ")" << std::endl;
227 // exit(-1);
228 // }
229 // return constants::pi / std::sqrt(tosqrt);
230 //}
231 //
233 //Mdouble HertzianViscoelasticNormalSpecies::getRestitutionCoefficient(Mdouble mass)
234 //{
235 // return std::exp(-dissipation_ / mass * getCollisionTime(mass));
236 //}
237 //
239 //Mdouble HertzianViscoelasticNormalSpecies::getMaximumVelocity(Mdouble radius, Mdouble mass)
240 //{
241 // return radius * std::sqrt(elasticModulus_ / (.5 * mass));
242 //}
243 //
245 //void HertzianViscoelasticNormalSpecies::setStiffnessAndRestitutionCoefficient(Mdouble k_, Mdouble eps, Mdouble mass)
246 //{
247 // elasticModulus_ = k_;
248 // dissipation_ = -std::sqrt(2.0 * mass * elasticModulus_ / (constants::sqr_pi + mathsFunc::square(log(eps)))) * log(eps);
249 //}
250 //
252 //void HertzianViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(Mdouble tc, Mdouble eps, Mdouble mass)
253 //{
254 // dissipation_ = -mass / tc * std::log(eps);
255 // elasticModulus_ = .5 * mass * (mathsFunc::square(constants::pi / tc) + mathsFunc::square(dissipation_ / mass));
256 //}
257 //
261 //void HertzianViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(Mdouble collisionTime, Mdouble restitutionCoefficient, Mdouble mass1, Mdouble mass2)
262 //{
263 // Mdouble reduced_mass = mass1 * mass2 / (mass1 + mass2);
264 // setCollisionTimeAndRestitutionCoefficient(collisionTime, restitutionCoefficient, 2.0 * reduced_mass);
265 //}
266 
274 {
277 }
278 
285  Mdouble relativeVelocity) const
286 {
287  // Here is a very nice paper describing contact modelling
288  // http://people.ds.cam.ac.uk/jae1001/CUS/research/pfizer/Antypov_Elliott_EPL_2011.pdf
289  //caution: this function assumes the contact is elastic (no dissipation)
290  //Mdouble omega0 = 2.0/constants::sqrt_pi*sqrt(getEffectiveElasticModulus()/particleDensity)/relativeVelocity;
291  //Mdouble omega1 = sqrt(omega0*omega0-getDissipation()*getDissipation());
292  return 2.214 * pow(mathsFunc::square(particleDensity / getEffectiveElasticModulus()) / relativeVelocity, 0.2) *
293  particleDiameter;
294 }
295 
296 
BaseSpecies * getBaseSpecies() const
Definition: BaseForce.h:35
void setEffectiveElasticModulusAndPoissonRatio(Mdouble elasticModulus, Mdouble poissonRatio)
Allows the elastic modulus and the poisson ratio to be changed in order to compute the shear modulus...
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
void setEffectiveElasticModulusAndEffectiveShearModulus(Mdouble elasticModulus, Mdouble shearModulus)
Allows the elastic modulus and the shear modulus to be changed in order to compute the poisson ratio...
Mdouble elasticModulus_
The effective elastic modulus.
Mdouble getCollisionTime(Mdouble particleDiameter, Mdouble particleDensity, Mdouble relativeVelocity) const
Used in Species::getName to obtain a unique name for each Species.
void read(std::istream &is)
Reads the species properties from an input stream.
void setDissipation(Mdouble dissipation)
Allows the normal dissipation to be changed.
Mdouble log(Mdouble Power)
MindlinSpecies contains the parameters used to describe sliding friction.
BaseFrictionForce * getFrictionForce() const
Definition: BaseSpecies.h:150
HertzianViscoelasticNormalSpecies()
The default constructor.
Mdouble getDissipation() const
Allows the normal dissipation to be accessed.
static Mdouble average(Mdouble a, Mdouble b)
Returns the harmonic mean of two variables.
Definition: BaseSpecies.cc:110
void setEffectiveElasticModulus(Mdouble elasticModulus)
Allows the spring constant to be changed.
void mix(HertzianViscoelasticNormalSpecies *SBase, HertzianViscoelasticNormalSpecies *TBase)
Calculates collision time for two copies of a particle of given disp, k, mass.
void setEffectiveElasticModulusAndRestitutionCoefficient(Mdouble elasticModulus, Mdouble rest)
Allows the spring constant to be changed.
Mdouble getEffectiveElasticModulus() const
Allows the spring constant to be accessed.
Mdouble dissipation_
normal dissipation constant
std::string getBaseName() const
Used in Species::getName to obtain a unique name for each Species.
Defines the basic properties that a interactable object can have.
T square(const T val)
squares a number
Definition: ExtendedMath.h:106
void write(std::ostream &os) const
Writes the species properties to an output stream.
HertzianViscoelasticNormalSpecies contains the parameters used to describe a Hertzian normal force (T...
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
~HertzianViscoelasticNormalSpecies()
The default destructor.