MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LinearViscoelasticNormalSpecies.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>
30 #include "Logger.h"
31 #include "Particles/BaseParticle.h"
32 
33 class BaseInteractable;
34 
36  : BaseNormalForce()
37 {
38  stiffness_ = 0.0;
39  dissipation_ = 0.0;
40 #ifdef DEBUG_CONSTRUCTOR
41  std::cout<<"LinearViscoelasticNormalSpecies::LinearViscoelasticNormalSpecies() finished"<<std::endl;
42 #endif
43 }
44 
49  : BaseNormalForce(p)
50 {
53 #ifdef DEBUG_CONSTRUCTOR
54  std::cout<<"LinearViscoelasticNormalSpecies::LinearViscoelasticNormalSpecies(const LinearViscoelasticNormalSpecies &p) finished"<<std::endl;
55 #endif
56 }
57 
59 {
60 #ifdef DEBUG_DESTRUCTOR
61  std::cout<<"LinearViscoelasticNormalSpecies::~LinearViscoelasticNormalSpecies() finished"<<std::endl;
62 #endif
63 }
64 
68 void LinearViscoelasticNormalSpecies::write(std::ostream& os) const
69 {
70  os << " stiffness " << stiffness_
71  << " dissipation " << dissipation_;
72 }
73 
78 {
79  std::string dummy;
80  is >> dummy >> stiffness_
81  >> dummy >> dissipation_;
82 }
83 
88 {
89  return "LinearViscoelastic";
90 }
91 
94 {
95  if (new_k >= 0)
96  stiffness_ = new_k;
97  else
98  {
99  std::cerr << "Error in set_k" << std::endl;
100  exit(-1);
101  }
102 }
103 
106 {
107  return stiffness_;
108 }
109 
112 {
113  setStiffness(new_.k);
114  setDissipation(new_.disp);
115 }
116 
119 {
120  if (dissipation >= 0)
121  {
122  dissipation_ = dissipation;
123  }
124  else
125  {
126  std::cerr << "Error in setDissipation(" << dissipation << ")" << std::endl;
127  exit(-1);
128  }
129 }
130 
133 {
134  return dissipation_;
135 }
136 
140 {
141  if (getConstantRestitution()) mass = 1;
142  if (mass <= 0)
143  {
144  std::cerr << "Warning in getCollisionTime(" << mass
145  << ") mass is not set or has an unexpected value, (getCollisionTime(" << mass << "))" << std::endl;
146  }
147  if (stiffness_ <= 0)
148  {
149  std::cerr << "Warning in getCollisionTime(" << mass << ") stiffness=" << stiffness_
150  << " is not set or has an unexpected value, (getCollisionTime(" << mass << "), with stiffness="
151  << stiffness_ << ")" << std::endl;
152  }
153  if (dissipation_ < 0)
154  {
155  std::cerr << "Warning in getCollisionTime(" << mass << ") dissipation=" << dissipation_
156  << " is not set or has an unexpected value, (getCollisionTime(" << mass << "), with dissipation="
157  << dissipation_ << ")" << std::endl;
158  }
159  Mdouble tosqrt = stiffness_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass);
160  if (tosqrt <= 0)
161  {
162  std::cerr << "Warning in getCollisionTime(" << mass
163  << ") values for mass, stiffness and dissipation would lead to an overdamped system, (getCollisionTime("
164  << mass << "), with stiffness=" << stiffness_ << " and dissipation=" << dissipation_ << ")"
165  << std::endl;
166  }
167  return constants::pi / std::sqrt(tosqrt);
168 }
169 
172 {
173  if (getConstantRestitution()) mass = 1;
174  return std::exp(-dissipation_ / mass * getCollisionTime(mass));
175 }
176 
179 {
180  return radius * std::sqrt(stiffness_ / (.5 * mass));
181 }
182 
190 {
191  if (getConstantRestitution()) mass = 1;
192  stiffness_ = stiffness;
193  setRestitutionCoefficient(eps, mass);
194 }
195 
203 {
204  if (getConstantRestitution()) mass = 1;
205  if (eps == 0.0) {
206  dissipation_ = std::sqrt(2.0 * mass * getStiffness());
207  } else {
208  const Mdouble logEps = log(eps);
209  dissipation_ = -std::sqrt(2.0 * mass * getStiffness()
210  / (constants::sqr_pi + mathsFunc::square(logEps))) * logEps;
211  }
212 }
213 
214 void
216 {
217  Mdouble mass = p->getSpecies()->getDensity() * p->getVolume();
219 }
220 
221 
229 {
230  if (getConstantRestitution()) mass = 1;
231  if (eps == 0.0)
232  {
233  stiffness_ = .5 * mass * mathsFunc::square(constants::pi / tc);
234  dissipation_ = std::sqrt(2.0 * mass * stiffness_);
235  }
236  else
237  {
238  dissipation_ = -mass / tc * mathsFunc::log(eps);
239  stiffness_ = .5 * mass * (mathsFunc::square(constants::pi / tc)
240  + mathsFunc::square(dissipation_ / mass));
241  }
242  logger(INFO,
243  "setCollisionTimeAndRestitutionCoefficient: set stiffness to % and dissipation to % to obtain a collision time of % and a restitution coefficient of % for two particles of mass %",
244  getStiffness(), getDissipation(), tc, eps, mass);
245 }
246 
254  Mdouble restitutionCoefficient,
255  Mdouble mass1, Mdouble mass2)
256 {
257  Mdouble reduced_mass = mass1 * mass2 / (mass1 + mass2);
258  setCollisionTimeAndRestitutionCoefficient(collisionTime, restitutionCoefficient, 2.0 * reduced_mass);
259 }
260 
266 void
268 {
271 }
bool getConstantRestitution() const
Mdouble stiffness_
(normal) spring constant
virtual Mdouble getVolume() const
Get Particle volume function, which required a reference to the Species vector. It returns the volume...
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
void setRestitutionCoefficient(double eps, Mdouble mass)
Sets disp to obtain a restitution coefficient eps for a collision of two particles of mass m...
Mdouble getStiffness() const
Allows the spring constant to be accessed.
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
return type specifically for fuctions returning k and disp at once
Definition: Helpers.h:45
void setDissipation(Mdouble dissipation)
Allows the normal dissipation to be changed.
Mdouble exp(Mdouble Exponent)
Definition: ExtendedMath.cc:84
std::string getBaseName() const
Used in Species::getName to obtain a unique name for each Species.
Mdouble getDissipation() const
Allows the normal dissipation to be accessed.
Mdouble getRestitutionCoefficient(Mdouble mass) const
Calculates restitution coefficient for two copies of given disp, k, mass.
void mix(LinearViscoelasticNormalSpecies *SBase, LinearViscoelasticNormalSpecies *TBase)
creates default values for mixed species
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...
~LinearViscoelasticNormalSpecies()
The default destructor.
Mdouble log(Mdouble Power)
void setCollisionTimeAndRestitutionCoefficient(Mdouble tc, Mdouble eps, BaseParticle *p)
Sets k, disp such that it matches a given tc and eps for a collision of two copies of particle p...
LinearViscoelasticNormalSpecies contains the parameters used to describe a linear elastic-dissipative...
void setStiffness(Mdouble new_k)
Allows the spring constant to be changed.
Mdouble disp
Definition: Helpers.h:49
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
Mdouble getCollisionTime(Mdouble mass) const
Calculates collision time for two copies of a particle of given disp, k, mass.
Mdouble getMaximumVelocity(Mdouble radius, Mdouble mass) const
Calculates the maximum velocity allowed for a collision of two copies of P (for higher velocities par...
void read(std::istream &is)
Reads the species properties from an input stream.
LinearViscoelasticNormalSpecies()
The default constructor.
Defines the basic properties that a interactable object can have.
Mdouble getDensity() const
Allows density_ to be accessed.
MERCURY_DEPRECATED void setStiffnessAndDissipation(helpers::KAndDisp new_)
Allows the spring and dissipation constants to be changed simultaneously.
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.
const Mdouble sqr_pi
Definition: ExtendedMath.h:47