MercuryDPM  Alpha
 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-2014, 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 }
54 {
57 #ifdef DEBUG_CONSTRUCTOR
58  std::cout<<"FrictionInteraction::FrictionInteraction(const FrictionInteraction& p) finished"<<std::endl;
59 #endif
60 }
65 {
66 #ifdef DEBUG_DESTRUCTOR
67  std::cout<<"FrictionInteraction::~FrictionInteraction() finished"<<std::endl;
68 #endif
69 }
73 void FrictionInteraction::write(std::ostream& os) const
74 {
76  os << " rollingSpring " << rollingSpring_;
77  os << " torsionSpring " << torsionSpring_;
78 }
82 void FrictionInteraction::read(std::istream& is)
83 {
85  std::string dummy;
86  is >> dummy >> rollingSpring_;
87  is >> dummy >> torsionSpring_;
88 }
93 {
95 
96  const FrictionSpecies* species = getSpecies();
97  //If tangential forces are present
98  if (getAbsoluteNormalForce() == 0.0) return;
99 
100 
101  if (species->getRollingFrictionCoefficient() != 0.0)
102  {
103  if (species->getRollingStiffness() != 0.0)
104  {
105  Mdouble effectiveDiameter = 2.0*getEffectiveRadius();
106 
107  //From Luding2008, objective rolling velocity (eq 15) w/o 2.0!
108  Vec3D rollingRelativeVelocity = -effectiveDiameter *
109  Vec3D::cross(getNormal(), getP()->getAngularVelocity() - getI()->getAngularVelocity());
110 
111  if (dynamic_cast<BaseParticle*>(getI())==0) //if particle-wall
112  rollingSpringVelocity_= rollingRelativeVelocity;
113  else //if particle-particle
114  {
115  const Vec3D relativeVelocity = getP()->getVelocity() - getI()->getVelocity();
116  rollingSpringVelocity_ = rollingRelativeVelocity
117  - Vec3D::dot(rollingSpring_, relativeVelocity) / getDistance() * getNormal();
118  }
119 
120  //used to Integrate the spring
121  //rollingSpringVelocity_= rollingRelativeVelocity;
122  //integrate(getHandler()->timeStep_);
124 
125  //Calculate test force acting on P including viscous force
126  Vec3D rollingForce = - species->getRollingStiffness() * rollingSpring_ - species->getRollingDissipation() * rollingRelativeVelocity;
127 
128  //tangential forces are modelled by a spring-damper of elasticity kt and viscosity dispt (sticking),
129  //but the force is limited by Coulomb friction (rolling):
130  Mdouble rollingForceSquared = rollingForce.getLengthSquared();
131  if (rollingForceSquared <= mathsFunc::square(species->getRollingFrictionCoefficientStatic() * getAbsoluteNormalForce()))
132  {
133  //if sticking (|ft|<=mu*|fn|), apply the force
134  }
135  else
136  {
137  //if rolling, resize the tangential force such that |ft|=mu*|fn|
138  rollingForce *= species->getRollingFrictionCoefficient() * getAbsoluteNormalForce() / std::sqrt(rollingForceSquared);
139  //resize the tangential spring accordingly such ft=-k*delta-nu*relVel
140  rollingSpring_ = -(rollingForce + species->getRollingDissipation() * rollingRelativeVelocity) / species->getRollingStiffness();
141  }
142  //Add (virtual) rolling force to torque
143  addTorque(effectiveDiameter * Vec3D::cross(getNormal(), rollingForce));
144  }
145  else //if no spring stiffness is set
146  {
147  std::cerr << "FrictionInteraction::computeFrictionForce(): Rolling stiffness is zero" << std::endl;
148  exit(-1);
149  }
150  } //end if rolling force
151 
152  if (species->getTorsionFrictionCoefficient() != 0.0)
153  {
154  if (species->getTorsionStiffness() != 0.0)
155  {
157  Mdouble effectiveDiameter = 2.0*getEffectiveRadius();
158 
159  //From Luding2008, spin velocity (eq 16) w/o 2.0!
160  Vec3D torsionRelativeVelocity = effectiveDiameter * Vec3D::dot(getNormal(), getP()->getAngularVelocity() - getI()->getAngularVelocity()) * getNormal();
161 
162  //Integrate the spring
163  torsionSpringVelocity_= torsionRelativeVelocity;
164  //integrate(getHandler()->timeStep_);
165  torsionSpring_ += Vec3D::dot(torsionSpring_ + torsionSpringVelocity_ * getHandler()->getDPMBase()->getTimeStep(), getNormal()) * getNormal();
166 
167  //Calculate test force acting on P including viscous force
168  Vec3D torsionForce = - species->getTorsionStiffness() * torsionSpring_ - species->getTorsionDissipation() * torsionRelativeVelocity;
169 
170  //tangential forces are modelled by a spring-damper of elasticity kt and viscosity dispt (sticking),
171  //but the force is limited by Coulomb friction (torsion):
172  Mdouble torsionForceSquared = torsionForce.getLengthSquared();
173  if (torsionForceSquared <= mathsFunc::square(species->getTorsionFrictionCoefficientStatic() * getAbsoluteNormalForce()))
174  {
175  //if sticking (|ft|<=mu*|fn|), apply the force
176  }
177  else
178  {
179  //if torsion, resize the tangential force such that |ft|=mu*|fn|
180  torsionForce *= species->getTorsionFrictionCoefficient() * getAbsoluteNormalForce() / std::sqrt(torsionForceSquared);
181  //resize the tangential spring accordingly such ft=-k*delta-nu*relVel
182  torsionSpring_ = -(torsionForce + species->getTorsionDissipation() * torsionRelativeVelocity) / species->getTorsionStiffness();
183  }
184  //Add (virtual) rolling force to torque
185  addTorque(effectiveDiameter * torsionForce);
186  }
187  else //if no spring stiffness is set
188  {
189  std::cerr << "FrictionInteraction::computeFrictionForce(): Torsion stiffness is zero" << std::endl;
190  exit(-1);
191  }
192  } //end if torsion force
193 }
198 {
202 }
207 {
211 }
216 {
217  return dynamic_cast<const FrictionSpecies*>(getBaseSpecies());
218 }
223 {
224  return "Friction";
225 }
230 {
232  //rollingSpring_=-rollingSpring_;
233  //rollingSpringVelocity_=-rollingSpringVelocity_;
236 }
237 
242 {
244  rollingSpring_=rotationMatrix*rollingSpring_;
246  torsionSpring_=rotationMatrix*torsionSpring_;
248 }
249 
251 {
252  return rollingSpring_;
253 }
void read(std::istream &is)
Interaction read function, which accepts an std::istream as input.
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.cc:291
void rotateHistory(Matrix3D &rotationMatrix)
When periodic particles are used, some interactions need certain history properties rotated (e...
void reverseHistory()
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_.
Mdouble getRollingFrictionCoefficient() const
Allows the (dynamic) Coulomb friction coefficient to be accessed.
void write(std::ostream &os) const
Interaction write function, which accepts an std::ostream as input.
const FrictionSpecies * getSpecies() const
Returns a const pointer of type FrictionSpecies*.
void integrate(Mdouble timeStep)
Computes the amount of compression in all the springs, i.e., increments the rollingSpring_, slidingSpring_ (see SlidingFrictionInteraction.cc) and torsionSpring_.
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.
Mdouble getTorsionFrictionCoefficient() const
Allows the (dynamic) Coulomb torsion friction coefficient to be accessed.
void reverseHistory()
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_.
FrictionSpecies contains the parameters used to describe sliding, rolling and torsional friction...
Mdouble getElasticEnergy() const
Returns the global amount of energy stored in all the springs (rolling, sliding and torsional)...
Computes the forces corresponding to sliding friction.
double Mdouble
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:52
void addTorque(Vec3D torque)
Mdouble getTorsionStiffness() const
Allows the torsion stiffness to be accessed.
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:167
Mdouble getRollingDissipation() const
Allows the tangential viscosity to be accessed.
Mdouble getRollingFrictionCoefficientStatic() const
Allows the static Coulomb rolling friction coefficient to be accessed.
T square(T val)
squares a number
Definition: ExtendedMath.h:91
void read(std::istream &is)
Interaction read function, which accepts an std::istream as input.
Stores information about interactions between two interactable objects; often particles but could be ...
Mdouble getRollingStiffness() const
Allows the spring constant to be accessed.
std::string getBaseName() const
Returns interaction name/type.
virtual ~FrictionInteraction()
Destructor.
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...
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
const Vec3D & getNormal() const
Gets the normal vector between the two interacting objects.
void integrate(Mdouble timeStep)
Increments the amount of compression in sliding spring.
void rotateHistory(Matrix3D &rotationMatrix)
When periodic particles are used, some interactions need certain history properties rotated (e...
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:255
FrictionInteraction(BaseInteractable *P, BaseInteractable *I, Mdouble timeStamp)
Constructor.
BaseInteractable * getI()
void write(std::ostream &os) const
Interaction print function, which accepts an std::ostream as input.
Mdouble getDistance() const
Returns an Mdouble which is the norm (length) of distance vector.
Mdouble getElasticEnergy() const
Returns the amount of elastic energy stored in sliding spring.
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().
Defines the basic properties that a interactable object can have.
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
Implementation of a 3D matrix.
Definition: Matrix.h:36
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
Mdouble getTorsionDissipation() const
Allows the torsion viscosity to be accessed.
Implementation of a 3D vector (by Vitaliy).
Definition: Vector.h:45
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:543
Mdouble getTimeStep() const
Allows the time step dt to be accessed.
Definition: DPMBase.cc:465
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...
Mdouble getAbsoluteNormalForce() const
Returns the absolute value of the norm (length) of the Normal force vector.