MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LinearPlasticViscoelasticInteraction.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 
28 #include "InteractionHandler.h"
29 #include "Particles/BaseParticle.h"
30 #include <iomanip>
31 #include <fstream>
33 #include <cmath> // std::max
34 
41  unsigned timeStamp)
42  : BaseInteraction(P, I, timeStamp)
43 {
44  maxOverlap_ = 0;
45 #ifdef DEBUG_CONSTRUCTOR
46  std::cout<<"LinearPlasticViscoelasticInteraction::LinearPlasticViscoelasticInteraction() finished"<<std::endl;
47 #endif
48 }
49 
50 //used for mpi
52  : BaseInteraction()
53 {
54  maxOverlap_ = 0;
55 #ifdef DEBUG_CONSTRUCTOR
56  std::cout<<"LinearPlasticViscoelasticInteraction::LinearPlasticViscoelasticInteraction() finished"<<std::endl;
57 #endif
58 }
59 
65  : BaseInteraction(p)
66 {
68 #ifdef DEBUG_CONSTRUCTOR
69  std::cout<<"LinearPlasticViscoelasticInteraction::LinearPlasticViscoelasticInteraction(const LinearPlasticViscoelasticInteraction &p finished"<<std::endl;
70 #endif
71 }
72 
77 {
78 #ifdef DEBUG_DESTRUCTOR
79  std::cout<<"LinearPlasticViscoelasticInteraction::~LinearPlasticViscoelasticInteraction() finished"<<std::endl;
80 #endif
81 }
82 
87 void LinearPlasticViscoelasticInteraction::write(std::ostream& os) const
88 {
90  os << " maxOverlap " << maxOverlap_;
91 }
92 
98 {
100  std::string dummy;
101  is >> dummy >> maxOverlap_;
102  //helpers::readOptionalVariable<Mdouble>(is, "maxOverlap", maxOverlap_);
103 }
104 
109 {
110  return "LinearPlasticViscoelastic";
111 }
112 
117 {
118  // Compute the relative velocity vector of particle P w.r.t. I
120  getP()->getVelocityAtContact(getContactPoint()) - getI()->getVelocityAtContact(getContactPoint()));
121  // Compute the projection of vrel onto the normal (can be negative)
123 
124  if (getOverlap() > 0) //if contact forces
125  {
127 
128  //calculate the overlap above which the max. unloading stiffness becomes active (the 'fluid branch')
129  // Modified by Paolo
130  // Previously:
131  // const Mdouble effectiveDiameter = species->getConstantRestitution()?1.0:(2.0 * getEffectiveRadius());
132  // This has been modified so that, independently on whether or not the users sets constant restitution,
133  // this parameter is set to its relative value, i.e. the penetration depth max divided by the radius of a single particle.
134  const Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
135  const Mdouble deltaStar = (species->getUnloadingStiffnessMax()
136  / (species->getUnloadingStiffnessMax() - species->getLoadingStiffness()))
137  * species->getPenetrationDepthMax() * effectiveDiameter;
138 
139  //increase max overlap if necessary
140  if (getOverlap() > getMaxOverlap())
141  setMaxOverlap(std::min(deltaStar, getOverlap()));
142 
143  //calculate the unloading stiffness
144  const Mdouble unloadingStiffness =
146  (species->getLoadingStiffness() +
147  (species->getUnloadingStiffnessMax() - species->getLoadingStiffness()) * (getMaxOverlap() / deltaStar));
148 
149  //calculate the overlap where the force is zero
150  const Mdouble equilibriumOverlap =
151  (unloadingStiffness - species->getLoadingStiffness()) / unloadingStiffness * maxOverlap_;
152 
153  //compute elastic force
154  Mdouble normalForce = unloadingStiffness * (getOverlap() - equilibriumOverlap);
155 
156  //decrease max overlap if necessary
157  const Mdouble cohesiveForce = -species->getCohesionStiffness() * getOverlap();
158  if (normalForce < cohesiveForce)
159  {
160  setMaxOverlap((unloadingStiffness + species->getCohesionStiffness())
161  / (unloadingStiffness - species->getLoadingStiffness()) * getOverlap());
162  //only necessary because the timeStep is finite:
163  normalForce = cohesiveForce;
164  }
165 
166  //add dissipative force
167  normalForce -= species->getDissipation() * getNormalRelativeVelocity();
168 
169  if (species->getConstantRestitution()) normalForce *= 2.0*getEffectiveMass();
170  setAbsoluteNormalForce(std::abs(normalForce)); //used for further corce calculations;
171  setForce(getNormal() * normalForce);
172  setTorque(Vec3D(0.0, 0.0, 0.0));
173  }
174  else
175  {
177  setForce(Vec3D(0.0, 0.0, 0.0));
178  setTorque(Vec3D(0.0, 0.0, 0.0));
179  }
180 }
181 
186 {
187  Mdouble energy = getOverlap() > 0 ? 0.5 * (getSpecies()->getLoadingStiffness() * mathsFunc::square(getOverlap())) : 0.0;
189  return energy;
191 }
192 
197 {
199 ;
200 }
201 
206 {
207  return maxOverlap_;
208 }
209 
214 {
215  maxOverlap_ = maxOverlap;
216 }
217 
222 {
224  // Modified by Paolo: read comment in computeNormalForce().
225  // previously: const Mdouble effectiveDiameter = species->getConstantRestitution()?1.0:(2.0 * getEffectiveRadius());
226  const Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
227  Mdouble deltaMaxFluid = species->getPenetrationDepthMax() * effectiveDiameter /
228  (1.0 - species->getLoadingStiffness() / species->getUnloadingStiffnessMax());
229  if (getOverlap() > deltaMaxFluid)
230  return species->getUnloadingStiffnessMax();
231  else
232  return species->getLoadingStiffness() +
233  (species->getUnloadingStiffnessMax() - species->getLoadingStiffness()) * getMaxOverlap() / deltaMaxFluid;
234 }
235 
bool getConstantRestitution() const
virtual std::string getBaseName() const
Returns the name of the interaction.
Mdouble getEffectiveRadius() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getUnloadingStiffnessMax() const
Returns the maximum unloading stiffness of the linear plastic-viscoelastic normal force...
const Vec3D & getRelativeVelocity() const
Returns a constant reference to a vector of relative velocity.
Mdouble getElasticEnergy() const override
Computes and returns the amount of elastic energy stored in the spring.
void setRelativeVelocity(Vec3D relativeVelocity)
set the relative velocity of the current of the interactions.
BaseNormalForce * getNormalForce() const
Definition: BaseSpecies.h:148
void setForce(Vec3D force)
set total force (this is used by the normal force, tangential forces are added use addForce) ...
const Vec3D & getContactPoint() const
Gets constant reference to contact point (vector).
BaseInteractable * getI()
Returns a pointer to the second object involved in the interaction (often a wall or a particle)...
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
void setNormalRelativeVelocity(Mdouble normalRelativeVelocit)
set the normal component of the relative velocity.
Stores information about interactions between two interactable objects; often particles but could be ...
const Vec3D & getNormal() const
Gets the normal vector between the two interacting objects.
void computeNormalForce()
Creates a copy of an object of this class. (Deep copy)
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
Mdouble getCohesionStiffness() const
Returns the cohesive stiffness of the linear plastic-viscoelastic normal force.
Mdouble getNormalRelativeVelocity() const
Returns a double which is the norm (length) of the relative velocity vector.
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
void setTorque(Vec3D torque)
set the total force (this is used by the normal force, tangential torques are added use addTorque) ...
Mdouble getDissipation() const
Allows the normal dissipation to be accessed.
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
void setAbsoluteNormalForce(Mdouble absoluteNormalForce)
the absolute values of the norm (length) of the normal force
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
void write(std::ostream &os) const override
Interaction write function, which accepts an std::ostream as input.
const LinearPlasticViscoelasticNormalSpecies * getSpecies() const
Defines the basic properties that a interactable object can have.
LinearPlasticViscoelasticNormalSpecies contains the parameters used to describe a plastic-cohesive no...
Definition: Vector.h:49
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:104
void write(std::ostream &os) const override
Interaction print function, which accepts an std::ostream as input.
Computes normal forces in case of a linear plastic visco-elastic interaction.
Mdouble getEffectiveMass() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
Mdouble getLoadingStiffness() const
Returns the loading stiffness of the linear plastic-viscoelastic normal force.