MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MindlinRollingTorsionInteraction.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  unsigned timeStamp)
42  : BaseInteraction(P, I, timeStamp), MindlinInteraction(P, I, timeStamp)
43 {
46 #ifdef DEBUG_CONSTRUCTOR
47  std::cout<<"MindlinRollingTorsionInteraction::MindlinRollingTorsionInteraction() finished"<<std::endl;
48 #endif
49 }
50 
51 //used for mpi
54 {
57 #ifdef DEBUG_CONSTRUCTOR
58  std::cout<<"MindlinRollingTorsionInteraction::MindlinRollingTorsionInteraction() finished"<<std::endl;
59 #endif
60 }
61 
67 {
70 #ifdef DEBUG_CONSTRUCTOR
71  std::cout<<"MindlinRollingTorsionInteraction::MindlinRollingTorsionInteraction(const MindlinRollingTorsionInteraction& p) finished"<<std::endl;
72 #endif
73 }
74 
79 {
80 #ifdef DEBUG_DESTRUCTOR
81  std::cout<<"MindlinRollingTorsionInteraction::~MindlinRollingTorsionInteraction() finished"<<std::endl;
82 #endif
83 }
84 
88 void MindlinRollingTorsionInteraction::write(std::ostream& os) const
89 {
91  os << " rollingSpring " << rollingSpring_;
92  os << " torsionSpring " << torsionSpring_;
93 }
94 
99 {
101  std::string dummy;
102  is >> dummy >> rollingSpring_;
103  is >> dummy >> torsionSpring_;
104 }
105 
110 {
112 
113  const MindlinRollingTorsionSpecies* species = getSpecies();
114  //If tangential forces are present
115  if (getAbsoluteNormalForce() == 0.0) return;
116 
117  if (species->getRollingFrictionCoefficient() != 0.0)
118  {
119  Mdouble rollingStiffness = tangentialStiffnessZero_;
120  double effectiveRadius = getEffectiveRadius();
121 
122  //From Luding2008, objective rolling velocity (eq 15) w/o 2.0!
123  Vec3D rollingRelativeVelocity = - effectiveRadius *
125  getP()->getAngularVelocity() - getI()->getAngularVelocity());
126 
127  if (dynamic_cast<BaseParticle*>(getI()) == nullptr) //if particle-wall
128  rollingSpringVelocity_ = rollingRelativeVelocity;
129  else //if particle-particle
130  {
131  const Vec3D relativeVelocity = getP()->getVelocity() - getI()->getVelocity();
132  rollingSpringVelocity_ = rollingRelativeVelocity
133  - Vec3D::dot(rollingSpring_, relativeVelocity) / getDistance() * getNormal();
134  }
135 
136  //used to Integrate the spring
137  //rollingSpringVelocity_= rollingRelativeVelocity;
138  //integrate(getHandler()->timeStep_);
140 
141  //Calculate test force acting on P including viscous force
142  Vec3D rollingForce = -rollingStiffness * rollingSpring_ -
143  species->getRollingDissipation() * rollingRelativeVelocity;
144 
145  //tangential forces are modelled by a spring-damper of elasticity kt and viscosity dispt (sticking),
146  //but the force is limited by Coulomb friction (rolling):
147  Mdouble rollingForceSquared = rollingForce.getLengthSquared();
148  if (rollingForceSquared <=
150  {
151  //if sticking (|ft|<=mu*|fn|), apply the force
152  }
153  else
154  {
155  //if rolling, resize the tangential force such that |ft|=mu*|fn|
156  rollingForce *= species->getRollingFrictionCoefficient() * getAbsoluteNormalForce() /
157  std::sqrt(rollingForceSquared);
158  //resize the tangential spring accordingly such ft=-k*delta-nu*relVel
159  rollingSpring_ = -(rollingForce + species->getRollingDissipation() * rollingRelativeVelocity) /
160  rollingStiffness;
161  }
162  //Add (virtual) rolling force to torque
163  addTorque(effectiveRadius * Vec3D::cross(getNormal(), rollingForce));
164  } //end if rolling force
165 
166  if (species->getTorsionFrictionCoefficient() != 0.0)
167  {
168  Mdouble torsionStiffness = tangentialStiffnessZero_;
170  Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
171 
172  //From Luding2008, spin velocity (eq 16) w/o 2.0!
173  Vec3D torsionRelativeVelocity = effectiveDiameter * Vec3D::dot(getNormal(), getP()->getAngularVelocity() -
174  getI()->getAngularVelocity()) *
175  getNormal();
176 
177  //Integrate the spring
178  torsionSpringVelocity_ = torsionRelativeVelocity;
179  //integrate(getHandler()->timeStep_);
180  torsionSpring_ +=
181  Vec3D::dot(torsionSpring_ + torsionSpringVelocity_ * getHandler()->getDPMBase()->getTimeStep(),
182  getNormal()) * getNormal();
183 
184  //Calculate test force acting on P including viscous force
185  Vec3D torsionForce = -torsionStiffness * torsionSpring_ -
186  species->getTorsionDissipation() * torsionRelativeVelocity;
187 
188  //tangential forces are modelled by a spring-damper of elasticity kt and viscosity dispt (sticking),
189  //but the force is limited by Coulomb friction (torsion):
190  Mdouble torsionForceSquared = torsionForce.getLengthSquared();
191  if (torsionForceSquared <=
193  {
194  //if sticking (|ft|<=mu*|fn|), apply the force
195  }
196  else
197  {
198  //if torsion, resize the tangential force such that |ft|=mu*|fn|
199  torsionForce *= species->getTorsionFrictionCoefficient() * getAbsoluteNormalForce() /
200  std::sqrt(torsionForceSquared);
201  //resize the tangential spring accordingly such ft=-k*delta-nu*relVel
202  torsionSpring_ = -(torsionForce + species->getTorsionDissipation() * torsionRelativeVelocity) /
203  torsionStiffness;
204  }
205  //Add (virtual) rolling force to torque
206  addTorque(effectiveDiameter * torsionForce);
207  } //end if torsion force
208 }
209 
214 {
218 }
219 
224 {
228 }
229 
234 {
235  return static_cast<const MindlinRollingTorsionSpecies*>(getBaseSpecies()->getFrictionForce());
236 ;
237 }
238 
243 {
244  return "Friction";
245 }
246 
251 {
253  //rollingSpring_=-rollingSpring_;
254  //rollingSpringVelocity_=-rollingSpringVelocity_;
257 }
258 
263 {
264  MindlinInteraction::rotateHistory(rotationMatrix);
265  rollingSpring_ = rotationMatrix * rollingSpring_;
267  torsionSpring_ = rotationMatrix * torsionSpring_;
269 }
270 
272 {
273  return rollingSpring_;
274 }
275 
277 {
278  return torsionSpring_;
279 }
280 
282 {
283  rollingSpring_ = rollingSpring;
284 }
285 
287 {
288  torsionSpring_ = torsionSpring;
289 }
290 
291 
292 
293 
294 
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().
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_.
Mdouble getEffectiveRadius() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
void write(std::ostream &os) const override
Interaction print function, which accepts an std::ostream as input.
InteractionHandler * getHandler() const
Gets a point to the interaction handlers to which this interaction belongs.
double Mdouble
Definition: GeneralDefine.h:34
Computes the forces corresponding to sliding friction.
Vec3D torsionSpring_
Stores the amount of torsional spring compression. Set in integrate(), used in computing frictional f...
Mdouble getTorsionFrictionCoefficient() const
Allows the (dynamic) Coulomb torsion friction coefficient to be accessed.
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
void addTorque(Vec3D torque)
const MindlinRollingTorsionSpecies * getSpecies() const
Returns a const pointer of type FrictionSpecies*.
Mdouble getRollingFrictionCoefficient() const
Allows the (dynamic) Coulomb friction coefficient to be accessed.
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
Vec3D rollingSpring_
Stores the amount of rolling spring compression. Set in integrate(), used in computing frictional for...
Mdouble getRollingDissipation() const
Allows the tangential viscosity to be accessed.
void integrate(Mdouble timeStep) override
Increments the amount of compression in sliding spring.
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 computeFrictionForce()
Computes the forces arising due to all three types of friction, i.e., sliding, rolling and torsional...
Mdouble getElasticEnergy() const override
Returns the global amount of energy stored in all the springs (rolling, sliding and torsional)...
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
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_.
BaseFrictionForce * getFrictionForce() const
Definition: BaseSpecies.h:150
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
void rotateHistory(Matrix3D &rotationMatrix) override
When periodic particles are used, some interactions need certain history properties rotated (e...
void rotateHistory(Matrix3D &rotationMatrix) override
When periodic particles are used, some interactions need certain history properties rotated (e...
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().
void computeFrictionForce()
Computes the tangential force generated due to compression in the sliding spring. Does take into acco...
Mdouble getElasticEnergy() const override
Returns the amount of elastic energy stored in sliding spring.
Mdouble getTorsionFrictionCoefficientStatic() const
Allows the static Coulomb torsion friction coefficient to be accessed.
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
void integrate(Mdouble timeStep) override
Computes the amount of compression in all the springs, i.e., increments the rollingSpring_, slidingSpring_ (see MindlinInteraction.cc) and torsionSpring_.
std::string getBaseName() const
Returns interaction name/type.
This class allows one to take all three types of frictional interactions into account. The sliding, rolling and torsional frictional interaction. See.
Mdouble getDistance() const
Returns an Mdouble which is the norm (length) of distance vector.
Mdouble getTorsionDissipation() const
Allows the torsion viscosity to be accessed.
Mdouble getRollingFrictionCoefficientStatic() const
Allows the static Coulomb rolling friction coefficient to be accessed.
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.
Implementation of a 3D matrix.
Definition: Matrix.h:37
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
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
MindlinRollingTorsionSpecies contains the parameters used to describe sliding, rolling and torsional ...
void write(std::ostream &os) const override
Interaction write function, which accepts an std::ostream as input.
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.