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  std::cerr << "Error in setEffectiveElasticModulus" << std::endl;
103  exit(-1);
104  }
105 }
106 
109 {
110  // Here is a very nice paper describing contact modelling
111  // http://people.ds.cam.ac.uk/jae1001/CUS/research/pfizer/Antypov_Elliott_EPL_2011.pdf
112  // see also: https://answers.launchpad.net/yade/+question/235934
113  if (elasticModulus < 0.0)
114  logger(ERROR,
115  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] elasticModulus % should be nonnegative",
116  elasticModulus);
117 
118  else if (rest < 0.0 || rest > 1.0)
119  logger(ERROR,
120  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] rest % should be between 0 and 1 (inclusive)",
121  rest);
122 
123  else
124  {
125  elasticModulus_ = elasticModulus;
126  if (rest > 0.0)
127  {
128  Mdouble logRestSquared = log(rest) * log(rest);
129  dissipation_ = sqrt(5.0 * logRestSquared / (logRestSquared + constants::sqr_pi));
130  }
131  else
132  dissipation_ = sqrt(5.0);
133 
134  //logger(INFO, "Effective elastic modulus %", elasticModulus_);
135  //logger(INFO, "Set dissipation % to match restitution coefficient %", dissipation_, rest);
136  }
137 }
138 
141 {
142  if (elasticModulus < 0.0)
143  logger(ERROR,
144  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] elasticModulus % should be nonnegative",
145  elasticModulus);
146 
147  else if (poissonRatio < 0.0 || poissonRatio > 1.0)
148  logger(ERROR,
149  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndRestitutionCoefficient] poissonRatio % should be between 0 and 1 (inclusive)",
150  poissonRatio);
151  else
152  {
153  elasticModulus_ = elasticModulus;
154  auto mindlin = dynamic_cast<MindlinSpecies*>(getBaseSpecies());
155  logger.assert(mindlin, "Please define HertzianViscoelasticMindlinSpecies to use this setter");
156  mindlin->setEffectiveShearModulus(elasticModulus / 2 * (1 + poissonRatio));
157  }
158 }
159 
162 {
163  if (elasticModulus < 0.0)
164  logger(ERROR,
165  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndEffectiveShearModulus] elasticModulus % should be nonnegative",
166  elasticModulus);
167 
168  else if (shearModulus < 0.0)
169  logger(ERROR,
170  "[HertzianViscoelasticNormalSpecies::setEffectiveElasticModulusAndEffectiveShearModulus] shearModulus % should be nonnegative",
171  shearModulus);
172  else
173  {
174  elasticModulus_ = elasticModulus;
175  auto mindlin = dynamic_cast<MindlinSpecies*>(getBaseSpecies()->getFrictionForce());
176  logger.assert(mindlin, "Please define HertzianViscoelasticMindlinSpecies to use this setter");
177  mindlin->setEffectiveShearModulus(shearModulus);
178  }
179 }
180 
183 {
184  return elasticModulus_;
185 }
186 
189 {
190  if (dissipation >= 0)
191  {
192  dissipation_ = dissipation;
193  }
194  else
195  {
196  std::cerr << "Error in setDissipation(" << dissipation << ")" << std::endl;
197  exit(-1);
198  }
199 }
200 
203 {
204  return dissipation_;
205 }
206 
208 //Mdouble HertzianViscoelasticNormalSpecies::getCollisionTime(Mdouble mass)
209 //{
210 // if (mass <= 0)
211 // {
212 // std::cerr << "Error in getCollisionTime(Mdouble mass) mass is not set or has an unexpected value, (getCollisionTime(" << mass << "))" << std::endl;
213 // exit(-1);
214 // }
215 // if (elasticModulus_ <= 0)
216 // {
217 // std::cerr << "Error in getCollisionTime(Mdouble mass) stiffness is not set or has an unexpected value, (getCollisionTime(" << mass << "), with stiffness=" << elasticModulus_ << ")" << std::endl;
218 // exit(-1);
219 // }
220 // if (dissipation_ < 0)
221 // {
222 // std::cerr << "Error in getCollisionTime(Mdouble mass) dissipation is not set or has an unexpected value, (getCollisionTime(" << mass << "), with dissipation=" << dissipation_ << ")" << std::endl;
223 // exit(-1);
224 // }
225 // Mdouble tosqrt = elasticModulus_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass);
226 // if (tosqrt <= 0)
227 // {
228 // 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;
229 // exit(-1);
230 // }
231 // return constants::pi / std::sqrt(tosqrt);
232 //}
233 //
235 //Mdouble HertzianViscoelasticNormalSpecies::getRestitutionCoefficient(Mdouble mass)
236 //{
237 // return std::exp(-dissipation_ / mass * getCollisionTime(mass));
238 //}
239 //
241 //Mdouble HertzianViscoelasticNormalSpecies::getMaximumVelocity(Mdouble radius, Mdouble mass)
242 //{
243 // return radius * std::sqrt(elasticModulus_ / (.5 * mass));
244 //}
245 //
247 //void HertzianViscoelasticNormalSpecies::setStiffnessAndRestitutionCoefficient(Mdouble k_, Mdouble eps, Mdouble mass)
248 //{
249 // elasticModulus_ = k_;
250 // dissipation_ = -std::sqrt(2.0 * mass * elasticModulus_ / (constants::sqr_pi + mathsFunc::square(log(eps)))) * log(eps);
251 //}
252 //
254 //void HertzianViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(Mdouble tc, Mdouble eps, Mdouble mass)
255 //{
256 // dissipation_ = -mass / tc * std::log(eps);
257 // elasticModulus_ = .5 * mass * (mathsFunc::square(constants::pi / tc) + mathsFunc::square(dissipation_ / mass));
258 //}
259 //
263 //void HertzianViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(Mdouble collisionTime, Mdouble restitutionCoefficient, Mdouble mass1, Mdouble mass2)
264 //{
265 // Mdouble reduced_mass = mass1 * mass2 / (mass1 + mass2);
266 // setCollisionTimeAndRestitutionCoefficient(collisionTime, restitutionCoefficient, 2.0 * reduced_mass);
267 //}
268 
276 {
279 }
280 
287  Mdouble relativeVelocity) const
288 {
289  // Here is a very nice paper describing contact modelling
290  // http://people.ds.cam.ac.uk/jae1001/CUS/research/pfizer/Antypov_Elliott_EPL_2011.pdf
291  //caution: this function assumes the contact is elastic (no dissipation)
292  //Mdouble omega0 = 2.0/constants::sqrt_pi*sqrt(getEffectiveElasticModulus()/particleDensity)/relativeVelocity;
293  //Mdouble omega1 = sqrt(omega0*omega0-getDissipation()*getDissipation());
294  return 2.214 * pow(mathsFunc::square(particleDensity / getEffectiveElasticModulus()) / relativeVelocity, 0.2) *
295  particleDiameter;
296 }
297 
298 
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")
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 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:104
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.