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  logger(ERROR, "setStiffness(%) argument has to be non-negative!", new_k);
100  }
101 }
102 
105 {
106  return stiffness_;
107 }
108 
111 {
112  setStiffness(new_.k);
113  setDissipation(new_.disp);
114 }
115 
118 {
119  if (dissipation >= 0)
120  {
121  dissipation_ = dissipation;
122  }
123  else
124  {
125  logger(ERROR, "Error in setDissipation(%)", dissipation);
126  }
127 }
128 
131 {
132  return dissipation_;
133 }
134 
138 {
139  if (getConstantRestitution()) mass = 1;
140  if (mass <= 0)
141  {
142  logger(ERROR, "Warning in getCollisionTime(%) mass is not set or has an unexpected value, "
143  "(getCollisionTime(%)", mass, mass);
144  }
145  if (stiffness_ <= 0)
146  {
147  logger(ERROR, "Error in getCollisionTime(%) stiffness=% is not set or has an unexpected value "
148  "(getCollisionTime(%), with stiffness=%)",
149  mass, stiffness_, mass, stiffness_);
150  }
151  if (dissipation_ < 0)
152  {
153  logger(ERROR, "Error in getCollisionTime(%) dissipation=% is not set or has an unexpected value "
154  "(getCollisionTime(%), with dissipation=%)",
155  mass, dissipation_, mass, dissipation_);
156  }
157  Mdouble tosqrt = stiffness_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass);
158  if (tosqrt <= 0)
159  {
160  logger(ERROR, "Warning in getCollisionTime(%) values for mass, stiffness and dissipation would lead to an "
161  "overdamped system, (getCollisionTime(%), with stiffness=% and dissipation=%",
162  mass, mass, stiffness_, dissipation_);
163  }
164  return constants::pi / std::sqrt(tosqrt);
165 }
166 
169 {
170  if (getConstantRestitution()) mass = 1;
171  return std::exp(-dissipation_ / mass * getCollisionTime(mass));
172 }
173 
176 {
177  return radius * std::sqrt(stiffness_ / (.5 * mass));
178 }
179 
187 {
188  if (getConstantRestitution()) mass = 1;
189  stiffness_ = stiffness;
190  setRestitutionCoefficient(eps, mass);
191 }
192 
200 {
201  if (getConstantRestitution()) mass = 1;
202  if (eps == 0.0) {
203  dissipation_ = std::sqrt(2.0 * mass * getStiffness());
204  } else {
205  const Mdouble logEps = log(eps);
206  dissipation_ = -std::sqrt(2.0 * mass * getStiffness()
207  / (constants::sqr_pi + mathsFunc::square(logEps))) * logEps;
208  }
209 }
210 
211 void
213 {
214  Mdouble mass = p->getSpecies()->getDensity() * p->getVolume();
216 }
217 
218 
226 {
227  if (getConstantRestitution()) mass = 1;
228  if (eps == 0.0)
229  {
230  stiffness_ = .5 * mass * mathsFunc::square(constants::pi / tc);
231  dissipation_ = std::sqrt(2.0 * mass * stiffness_);
232  }
233  else
234  {
235  dissipation_ = -mass / tc * mathsFunc::log(eps);
236  stiffness_ = .5 * mass * (mathsFunc::square(constants::pi / tc)
237  + mathsFunc::square(dissipation_ / mass));
238  }
239  logger(INFO,
240  "setCollisionTimeAndRestitutionCoefficient: set stiffness to % and dissipation to % to obtain a collision time of % and a restitution coefficient of % for two particles of mass %",
241  getStiffness(), getDissipation(), tc, eps, mass);
242 }
243 
251  Mdouble restitutionCoefficient,
252  Mdouble mass1, Mdouble mass2)
253 {
254  Mdouble reduced_mass = mass1 * mass2 / (mass1 + mass2);
255  setCollisionTimeAndRestitutionCoefficient(collisionTime, restitutionCoefficient, 2.0 * reduced_mass);
256 }
257 
263 void
265 {
268 }
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")
Definition of different loggers with certain modules. A user can define its own custom logger here...
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:106
void write(std::ostream &os) const
Writes the species properties to an output stream.
const Mdouble sqr_pi
Definition: ExtendedMath.h:47