MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SlidingFrictionInteraction.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 
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)
42 {
45 #ifdef DEBUG_CONSTRUCTOR
46  std::cout<<"SlidingFrictionInteraction::SlidingFrictionInteraction() finished"<<std::endl;
47 #endif
48 }
49 
50 //used for mpi
52  : BaseInteraction()
53 {
56 #ifdef DEBUG_CONSTRUCTOR
57  std::cout<<"SlidingFrictionInteraction::SlidingFrictionInteraction() finished"<<std::endl;
58 #endif
59 }
60 
65  : BaseInteraction(p)
66 {
69 #ifdef DEBUG_CONSTRUCTOR
70  std::cout<<"SlidingFrictionInteraction::SlidingFrictionInteraction(const SlidingFrictionInteraction& p) finished"<<std::endl;
71 #endif
72 }
73 
78 {
79 #ifdef DEBUG_DESTRUCTOR
80  std::cout<<"SlidingFrictionInteraction::~SlidingFrictionInteraction() finished"<<std::endl;
81 #endif
82 }
83 
87 void SlidingFrictionInteraction::write(std::ostream& os) const
88 {
89  //BaseInteraction::write(os);
90  os << " slidingSpring " << slidingSpring_;
91 }
92 
96 void SlidingFrictionInteraction::read(std::istream& is)
97 {
98  //BaseInteraction::read(is);
99  std::string dummy;
100  is >> dummy >> slidingSpring_;
101 }
102 
107 {
108  //If tangential forces are absent
109  if (getAbsoluteNormalForce() == 0.0) return;
110 
111  const SlidingFrictionSpecies* species = getSpecies();//dynamic_cast
113 
114  if (species->getSlidingFrictionCoefficient() != 0.0)
115  {
116  //Compute the tangential component of relativeVelocity_
117  // relativeVelocity = [v_p + (r_p-c) x w_p] - [v_i + (r_i-c) x w_i]
118  // tangentialRelativeVelocity = relativeVelocity - (relativeVelocity . n) n
119  Vec3D tangentialRelativeVelocity = getRelativeVelocity() - getNormal() * getNormalRelativeVelocity();
120 
121  if (species->getSlidingStiffness() != 0.0)
122  {
124  {
125  computeSlidingSpring(tangentialRelativeVelocity);
126  }
127  else
128  {
129  computeSlidingSpringSuperQuadric(tangentialRelativeVelocity);
130  }
131 
132  //Calculate test force acting on P including viscous force
134  species->getSlidingDissipation() * tangentialRelativeVelocity;
135  if (getBaseSpecies()->getNormalForce()->getConstantRestitution()) tangentialForce_ *=2.0*getEffectiveMass();
136 
137  //tangential forces are modelled by a spring-damper of elasticity kt and viscosity dispt (sticking),
138  //but the force is limited by Coulomb friction (sliding):
139  Mdouble tangentialForceSquared = tangentialForce_.getLengthSquared();
140  if (tangentialForceSquared <=
142  {
143  //if sticking (|ft|<=mu*|fn|), apply the force
145  }
146  else
147  {
148  //if sliding, resize the tangential force such that |ft|=mu*|fn|
150  std::sqrt(tangentialForceSquared);
152  //resize the tangential spring accordingly such ft=-k*delta-nu*relVel
153  slidingSpring_ = -(tangentialForce_ + species->getSlidingDissipation() * tangentialRelativeVelocity) /
154  species->getSlidingStiffness();
155  }
156  }
157  else //if no spring stiffness is set
158  {
159 // if (species->getSlidingDissipation()==0.0)
160 // {
161 // std::cerr << "SlidingFrictionInteraction::getForce(): warning: both sliding stiffness and dissipation are zero" << std::endl;
162 // }
163  Mdouble tangentialRelativeVelocitySquared = tangentialRelativeVelocity.getLengthSquared();
164  if (tangentialRelativeVelocitySquared * mathsFunc::square(species->getSlidingDissipation()) <=
166  tangentialForce_ = -species->getSlidingDissipation() * tangentialRelativeVelocity;
167  else //if sliding, set force to Coulomb limit
169  std::sqrt(tangentialRelativeVelocitySquared)) * tangentialRelativeVelocity;
170 
172  }
173  }
174 
175  /* And if species->getSlidingFrictionCoefficient() == 0.0, then there is no
176  * friction force at all. */
177 }
178 
179 void SlidingFrictionInteraction::computeSlidingSpring(const Vec3D& tangentialRelativeVelocity)
180 {
181  //used to Integrate the spring
182  if (dynamic_cast<BaseParticle*>(getI()) == nullptr) //if particle-wall
183  {
184  slidingSpringVelocity_ = tangentialRelativeVelocity;
185  }
186  else //if particle-particle
187  {
188  slidingSpringVelocity_ = (tangentialRelativeVelocity -
189  Vec3D::dot(slidingSpring_, getP()->getVelocity() - getI()->getVelocity()) *
190  getNormal() / getDistance());
191  }
192  // v_s = v_t - [xi . (v_p-v_i)/|r_pi|] n
193 
194 
195  //integrate(getHandler()->timeStep_);
196  // xi = xi' + dt v_s
198  // Stefan does [EJECE-12/2008] sth. like xi = xi' - (xi . n) n + dt*v_t, see SlidingFrictionSuperQuadricInteraction
199 }
200 
202 {
203  //used to Integrate the spring
204  if (dynamic_cast<BaseParticle*>(getI()) == nullptr) //if particle-wall
205  {
206  slidingSpringVelocity_ = tangentialRelativeVelocity;
207  }
208  else //if particle-particle
209  {
210  slidingSpringVelocity_ = tangentialRelativeVelocity;
211  }
212 
213  // From Stefan [EJECE-12/2008]:
214  // xi = xi' - (xi . n) n and renormalising, then add v_t * dt
215  logger(DEBUG, "sliding spring, normal before rotation-step: %, %", slidingSpring_, getNormal());
216  if (slidingSpring_.getLength() > 1e-10)
217  {
218  const Mdouble springLength = slidingSpring_.getLength();
221  slidingSpring_ *= springLength;
222  logger.assert(std::abs(slidingSpring_.getLength() - springLength) < 1e-10, "Spring length not the same after rotation");
223  }
224  logger(DEBUG, "sliding spring after rotation-step: %\n", slidingSpring_);
226  logger.assert(std::abs(Vec3D::dot(slidingSpring_, getNormal())) < 1e-10, "sliding spring not perpendicular to normal");
227 }
228 
229 
234 {
236 }
237 
242 {
243  if (getBaseSpecies()->getNormalForce()->getConstantRestitution()) {
245  } else {
247  }
248 }
249 
254 {
256  return -slidingSpring_.getLength();
257 }
258 
263 {
264  return tangentialForce_;
265 }
266 
271 {
272  return static_cast<const SlidingFrictionSpecies*>(getBaseSpecies()->getFrictionForce());
273 ;
274 }
275 
280 {
281  return "SlidingFriction";
282 }
283 
286 {
287  slidingSpring_ = slidingSpring;
288 }
289 
294 {
298 }
299 
301 {
302  slidingSpring_ = rotationMatrix * slidingSpring_;
304  tangentialForce_ = rotationMatrix * tangentialForce_;
305 }
306 
308 {
309  return slidingSpring_;
310 }
311 
313 {
314  slidingSpring_ += displacement;
315 }
316 
318 {
319  isSuperQuadricInteraction_ = isSuperQuadricInteraction;
320 }
void rotateHistory(Matrix3D &rotationMatrix) override
When periodic particles are used, some interactions need certain history properties rotated (e...
Mdouble getSlidingFrictionCoefficientStatic() const
Allows the static Coulomb friction coefficient to be accessed.
Vec3D slidingSpringVelocity_
Stores the rate at which the sliding spring compressed or relaxed. Set in the member function compute...
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
InteractionHandler * getHandler() const
Gets a point to the interaction handlers to which this interaction belongs.
double Mdouble
Definition: GeneralDefine.h:34
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
Mdouble getTangentialOverlap() const override
Returns the amount of tangential overlap which is needed by BaseInteraction::writeToFstat().
const Vec3D & getRelativeVelocity() const
Returns a constant reference to a vector of relative velocity.
Computes the forces corresponding to sliding friction.
void setSlidingSpring(Vec3D slidingSpring)
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
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_.
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
Stores information about interactions between two interactable objects; often particles but could be ...
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
Definition: Vector.cc:331
const Vec3D & getNormal() const
Gets the normal vector between the two interacting objects.
Vec3D tangentialForce_
Computes the tangential force such that . Set and computed in computeFrictionForce().
~SlidingFrictionInteraction() override
Destructor.
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
void computeSlidingSpringSuperQuadric(const Vec3D &tangentialRelativeVelocity)
void setIsSuperQuadricInteraction(bool isSuperQuadricInteraction)
void computeFrictionForce()
Computes the tangential force generated due to compression in the sliding spring. Does take into acco...
Mdouble getNormalRelativeVelocity() const
Returns a double which is the norm (length) of the relative velocity vector.
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.
const SlidingFrictionSpecies * getSpecies() const
Returns a const pointer of type SlidingFrictionSpecies*.
Vec3D slidingSpring_
Stores the amount of sliding spring ( ) compression from the expression . Set in the member function ...
const Vec3D getTangentialForce() const override
Returns the sliding friction force vector.
std::string getBaseName() const
Returns the type/name of interaction (sliding friction interaction)
void addForce(Vec3D force)
add an force increment to the total force.
SlidingFrictionSpecies contains the parameters used to describe sliding friction. ...
Mdouble getDistance() const
Returns an Mdouble which is the norm (length) of distance vector.
Mdouble getSlidingStiffness() const
Allows the spring constant to be accessed.
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:316
Mdouble getSlidingDissipation() const
Allows the tangential viscosity to be accessed.
void moveSlidingSpring(Vec3D displacement)
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
Definition: Vector.h:49
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
void write(std::ostream &os) const override
Interaction write function, which accepts an std::ostream as input.
void computeSlidingSpring(const Vec3D &tangentialRelativeVelocity)
Mdouble getEffectiveMass() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
Mdouble getSlidingFrictionCoefficient() const
Allows the (dynamic) Coulomb friction coefficient to be accessed.
Mdouble getAbsoluteNormalForce() const
Returns the absolute value of the norm (length) of the Normal force vector.