MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FrictionInteraction.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 
27 #include "FrictionInteraction.h"
29 #include "Particles/BaseParticle.h"
30 #include "InteractionHandler.h"
31 #include <iomanip>
32 #include <fstream>
33 #include <DPMBase.h>
34 
41  : BaseInteraction(P, I, timeStamp), SlidingFrictionInteraction(P, I, timeStamp)
42 {
45 #ifdef DEBUG_CONSTRUCTOR
46  std::cout<<"FrictionInteraction::FrictionInteraction() finished"<<std::endl;
47 #endif
48 }
49 
50 //used for mpi
53 {
56 #ifdef DEBUG_CONSTRUCTOR
57  std::cout<<"FrictionInteraction::FrictionInteraction() finished"<<std::endl;
58 #endif
59 }
60 
66 {
69 #ifdef DEBUG_CONSTRUCTOR
70  std::cout<<"FrictionInteraction::FrictionInteraction(const FrictionInteraction& p) finished"<<std::endl;
71 #endif
72 }
73 
78 {
79 #ifdef DEBUG_DESTRUCTOR
80  std::cout<<"FrictionInteraction::~FrictionInteraction() finished"<<std::endl;
81 #endif
82 }
83 
87 void FrictionInteraction::write(std::ostream& os) const
88 {
90  os << " rollingSpring " << rollingSpring_;
91  os << " torsionSpring " << torsionSpring_;
92 }
93 
97 void FrictionInteraction::read(std::istream& is)
98 {
100  std::string dummy;
101  is >> dummy >> rollingSpring_;
102  is >> dummy >> torsionSpring_;
103 }
104 
109 {
111 
112  const FrictionSpecies* species = getSpecies();
113  //If tangential forces are present
114  if (getAbsoluteNormalForce() == 0.0) return;
115 
116 
117  if (species->getRollingFrictionCoefficient() != 0.0)
118  {
119  if (species->getRollingStiffness() != 0.0)
120  {
121  Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
122 
123  //From Luding2008, objective rolling velocity (eq 15) w/o 2.0!
124  Vec3D rollingRelativeVelocity = -effectiveDiameter *
126  getP()->getAngularVelocity() - getI()->getAngularVelocity());
127 
128  if (dynamic_cast<BaseParticle*>(getI()) == nullptr) //if particle-wall
129  rollingSpringVelocity_ = rollingRelativeVelocity;
130  else //if particle-particle
131  {
132  const Vec3D relativeVelocity = getP()->getVelocity() - getI()->getVelocity();
133  rollingSpringVelocity_ = rollingRelativeVelocity
134  - Vec3D::dot(rollingSpring_, relativeVelocity) / getDistance() * getNormal();
135  }
136 
137  //used to Integrate the spring
138  //rollingSpringVelocity_= rollingRelativeVelocity;
139  //integrate(getHandler()->timeStep_);
141 
142  //Calculate test force acting on P including viscous force
143  Vec3D rollingForce = -species->getRollingStiffness() * rollingSpring_ -
144  species->getRollingDissipation() * rollingRelativeVelocity;
145  if (getBaseSpecies()->getNormalForce()->getConstantRestitution()) rollingForce *=2.0*getEffectiveMass();
146 
147  //tangential forces are modelled by a spring-damper of elasticity kt and viscosity dispt (sticking),
148  //but the force is limited by Coulomb friction (rolling):
149  Mdouble rollingForceSquared = rollingForce.getLengthSquared();
150  if (rollingForceSquared <=
152  {
153  //if sticking (|ft|<=mu*|fn|), apply the force
154  }
155  else
156  {
157  //if rolling, resize the tangential force such that |ft|=mu*|fn|
158  rollingForce *= species->getRollingFrictionCoefficient() * getAbsoluteNormalForce() /
159  std::sqrt(rollingForceSquared);
160  //resize the tangential spring accordingly such ft=-k*delta-nu*relVel
161  rollingSpring_ = -(rollingForce + species->getRollingDissipation() * rollingRelativeVelocity) /
162  species->getRollingStiffness();
163  }
164  //Add (virtual) rolling force to torque
165  addTorque(effectiveDiameter * Vec3D::cross(getNormal(), rollingForce));
166  }
167  else //if no spring stiffness is set
168  {
169  std::cerr << "FrictionInteraction::computeFrictionForce(): Rolling stiffness is zero" << std::endl;
170  exit(-1);
171  }
172  } //end if rolling force
173 
174  if (species->getTorsionFrictionCoefficient() != 0.0)
175  {
176  if (species->getTorsionStiffness() != 0.0)
177  {
179  Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
180 
181  //From Luding2008, spin velocity (eq 16) w/o 2.0!
182  Vec3D torsionRelativeVelocity = effectiveDiameter * Vec3D::dot(getNormal(), getP()->getAngularVelocity() -
183  getI()->getAngularVelocity()) *
184  getNormal();
185 
186  //Integrate the spring
187  torsionSpringVelocity_ = torsionRelativeVelocity;
188  //integrate(getHandler()->timeStep_);
190  Vec3D::dot(torsionSpring_ + torsionSpringVelocity_ * getHandler()->getDPMBase()->getTimeStep(),
191  getNormal()) * getNormal();
192 
193  //Calculate test force acting on P including viscous force
194  Vec3D torsionForce = -species->getTorsionStiffness() * torsionSpring_ -
195  species->getTorsionDissipation() * torsionRelativeVelocity;
196  if (getBaseSpecies()->getNormalForce()->getConstantRestitution()) torsionForce *=2.0*getEffectiveMass();
197 
198  //tangential forces are modelled by a spring-damper of elasticity kt and viscosity dispt (sticking),
199  //but the force is limited by Coulomb friction (torsion):
200  Mdouble torsionForceSquared = torsionForce.getLengthSquared();
201  if (torsionForceSquared <=
203  {
204  //if sticking (|ft|<=mu*|fn|), apply the force
205  }
206  else
207  {
208  //if torsion, resize the tangential force such that |ft|=mu*|fn|
209  torsionForce *= species->getTorsionFrictionCoefficient() * getAbsoluteNormalForce() /
210  std::sqrt(torsionForceSquared);
211  //resize the tangential spring accordingly such ft=-k*delta-nu*relVel
212  torsionSpring_ = -(torsionForce + species->getTorsionDissipation() * torsionRelativeVelocity) /
213  species->getTorsionStiffness();
214  }
215  //Add (virtual) rolling force to torque
216  addTorque(effectiveDiameter * torsionForce);
217  }
218  else //if no spring stiffness is set
219  {
220  std::cerr << "FrictionInteraction::computeFrictionForce(): Torsion stiffness is zero" << std::endl;
221  exit(-1);
222  }
223  } //end if torsion force
224 }
225 
230 {
234 }
235 
240 {
241  if (getBaseSpecies()->getNormalForce()->getConstantRestitution()) {
245  } else {
249  }
250 }
251 
256 {
257  return static_cast<const FrictionSpecies*>(getBaseSpecies()->getFrictionForce());
258 ;
259 }
260 
265 {
266  return "Friction";
267 }
268 
273 {
275  //rollingSpring_=-rollingSpring_;
276  //rollingSpringVelocity_=-rollingSpringVelocity_;
279 }
280 
285 {
287  rollingSpring_ = rotationMatrix * rollingSpring_;
289  torsionSpring_ = rotationMatrix * torsionSpring_;
291 }
292 
294 {
295  return rollingSpring_;
296 }
297 
299 {
300  return torsionSpring_;
301 }
302 
304 {
305  rollingSpring_ = rollingSpring;
306 }
307 
309 {
310  torsionSpring_ = torsionSpring;
311 }
312 
313 
314 
315 
316 
void rotateHistory(Matrix3D &rotationMatrix) override
When periodic particles are used, some interactions need certain history properties rotated (e...
Mdouble getRollingFrictionCoefficient() const
Allows the (dynamic) Coulomb friction coefficient to be accessed.
const FrictionSpecies * getSpecies() const
Returns a const pointer of type FrictionSpecies*.
Mdouble getEffectiveRadius() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
Vec3D getRollingSpring() const
InteractionHandler * getHandler() const
Gets a point to the interaction handlers to which this interaction belongs.
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getTorsionFrictionCoefficient() const
Allows the (dynamic) Coulomb torsion friction coefficient to be accessed.
FrictionSpecies contains the parameters used to describe sliding, rolling and torsional friction...
void reverseHistory() override
A useful feature if one wants to return to the initial state of the springs. However, reverse history decrements the current state to the state corresponding to previous time step. Decrements the state or value of rollingSpring_, torsionSpring_ and slidingSpring_.
Computes the forces corresponding to sliding friction.
Mdouble getElasticEnergy() const override
Returns the global amount of energy stored in all the springs (rolling, sliding and torsional)...
This class allows one to take all three types of frictional interactions into account. The sliding, rolling and torsional frictional interaction. See.
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
void addTorque(Vec3D torque)
void reverseHistory() override
A useful feature if one wants to return to the initial state of the spring. However, reverse history decrements the current state to the state corresponding to previous time step. Decrements the value of slidingSpring_.
void setRollingSpring(Vec3D rollingSpring)
Mdouble getTorsionStiffness() const
Allows the torsion stiffness to be accessed.
Mdouble getElasticEnergy() const override
Returns the amount of elastic energy stored in sliding spring.
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
Mdouble getRollingDissipation() const
Allows the tangential viscosity to be accessed.
Mdouble getRollingFrictionCoefficientStatic() const
Allows the static Coulomb rolling friction coefficient to be accessed.
~FrictionInteraction() override
Destructor.
Stores information about interactions between two interactable objects; often particles but could be ...
void setTorsionSpring(Vec3D torsionSpring)
const Vec3D & getNormal() const
Gets the normal vector between the two interacting objects.
Mdouble getRollingStiffness() const
Allows the spring constant to be accessed.
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
std::string getBaseName() const
Returns interaction name/type.
void computeFrictionForce()
Computes the tangential force generated due to compression in the sliding spring. Does take into acco...
void computeFrictionForce()
Computes the forces arising due to all three types of friction, i.e., sliding, rolling and torsional...
BaseFrictionForce * getFrictionForce() const
Definition: BaseSpecies.h:150
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
Mdouble getDistance() const
Returns an Mdouble which is the norm (length) of distance vector.
Vec3D rollingSpringVelocity_
Stores the rate at which the rolling spring compresses or relaxes. Set in computeFrictionForce(), used in computing the amount of compression in rolling spring. Used in integrate().
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:316
Defines the basic properties that a interactable object can have.
void integrate(Mdouble timeStep) override
Increments the amount of compression in sliding spring.
Implementation of a 3D matrix.
Definition: Matrix.h:37
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
Mdouble getTorsionDissipation() const
Allows the torsion viscosity to be accessed.
Definition: Vector.h:49
void rotateHistory(Matrix3D &rotationMatrix) override
When periodic particles are used, some interactions need certain history properties rotated (e...
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1211
Vec3D torsionSpring_
Stores the amount of torsional spring compression. Set in integrate(), used in computing frictional f...
Vec3D torsionSpringVelocity_
Stores the rate at which the torsional spring compresses or relaxes. Set in computeFrictionForce(), used in computing the amount of compression in torsion spring. Used in integrate().
Mdouble getTorsionFrictionCoefficientStatic() const
Allows the static Coulomb torsion friction coefficient to be accessed.
Vec3D rollingSpring_
Stores the amount of rolling spring compression. Set in integrate(), used in computing frictional for...
void write(std::ostream &os) const override
Interaction print function, which accepts an std::ostream as input.
void write(std::ostream &os) const override
Interaction write function, which accepts an std::ostream as input.
void integrate(Mdouble timeStep) override
Computes the amount of compression in all the springs, i.e., increments the rollingSpring_, slidingSpring_ (see SlidingFrictionInteraction.cc) and torsionSpring_.
Mdouble getEffectiveMass() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
Vec3D getTorsionSpring() const
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
Mdouble getAbsoluteNormalForce() const
Returns the absolute value of the norm (length) of the Normal force vector.