MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MindlinInteraction.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 "MindlinInteraction.h"
29 #include "Particles/BaseParticle.h"
30 #include "InteractionHandler.h"
31 #include <iomanip>
32 #include <DPMBase.h>
33 
40  : BaseInteraction(P, I, timeStamp)
41 {
43  //k_edit
44  //setting the "previousSlidingSpring_" initially to zero, as before a contact
45  //this parameter should (obviously) not carry a value.
47  //similarly, setting the turning point tangential forces to zero, as a "new"
48  //interaction will clearly have no such history
51  //k_edit
52  //setting the K_t0 parameter by default to zero...
54  //...and its previous value...
56 
58  //k_edit
59  //...and, similarly, the "normal" tangential stiffness, K_t
61  //k_edit
62  //setting the tangential force direction by default to -1 (i.e. acting against the motion!)
64  //k_edit
65  //setting to zero the temporary 'holders' which store intermediate values of tangential force
66  //and displacement during multiple-step force calculations
71  //Similarly for temporary turning point force values...
74  //...and the turning point displacement values
77  //k_edit
78  //setting the variable which stores the minimum displacement for simple loading to zero
80  //k_edit
81  //setting the priorLoadingFlag to zero as, when an interaction is created for the first time,
82  //there will, by definition, have been no prior loading
83  priorLoadingFlag_ = false;
84  //ensuring that the initial tangential velocity of a new interaction is zero until updated
86 
87 #ifdef DEBUG_CONSTRUCTOR
88  std::cout<<"MindlinInteraction::MindlinInteraction() finished"<<std::endl;
89 #endif
90 }
91 
96  : BaseInteraction(p)
97 {
99  //k_edit
100  //setting the K_t0 parameter to that possessed by the interaction
101  //to be copied....
103  //...and its previous value...
105  //...and the same for the "normal" tangential stiffness parameter
107  //k_edit
108  //allowing turning point forces to be copied when an interaction is copied
111  //k_edit
112  //allowing the direction of the tangential force to be "remembered" also
113  //by the copied interaction
115  //k_edit
116  //copying the temporary 'holders' which store intermediate values of tangential force
117  //and displacement during multiple-step force calculations
122  //and similarly for the temporary turning point variables
125  //...and the displacement equivalents
128  //k_edit
129  //copying the variable which stores the minimum displacement for simple loading
131  //k_edit
133  //copying the initial tangential velocity
135 
137 #ifdef DEBUG_CONSTRUCTOR
138  std::cout<<"MindlinInteraction::MindlinInteraction(const MindlinInteraction& p) finished"<<std::endl;
139 #endif
140 }
141 
146 {
147  //I don't see why we are restricted in the use of interactions
148 //#ifdef MERCURY_USE_MPI
149 // logger(FATAL,"MindlinInteractions are currently not implemented in parallel MercuryDPM");
150 //#endif
151 }
152 
157 {
158 #ifdef DEBUG_DESTRUCTOR
159  std::cout<<"MindlinInteraction::~MindlinInteraction() finished"<<std::endl;
160 #endif
161 }
162 
166 void MindlinInteraction::write(std::ostream& os) const
167 {
168  //BaseInteraction::write(os);
169  os << " slidingSpring " << slidingSpring_;
170 }
171 
175 void MindlinInteraction::read(std::istream& is)
176 {
177  //BaseInteraction::read(is);
178  std::string dummy;
179  is >> dummy >> slidingSpring_;
180 }
181 
182 //k_new
183 
184 
185 //k_edit
186 //Allows the K_t0 parameter (for Mindlin model) to be set
188 {
189  tangentialStiffnessZero_ = newKt0;
190 }
191 
192 //k_edit
193 //Allows the K_t0 parameter (for Mindlin model) to be accessed
195 {
197 }
198 
199 //k_edit
200 //Allows the "normal" K_t parameter (for Mindlin model) to be accessed
202 {
203  return tangentialStiffness_;
204 }
205 
206 //k_edit
207 //A quick way to update the K_t0 parameter when necessary
208 //**************k_TODO at present, only valid for particle-wall interactions and/or interactions*********************
209 //between similar particles; need to include the reduced parameters for dissimilar particles!
210 //takes as arguments the relevant particle radius and the relevant shear modulus
211 //2.0/3.0 as a factor from Di Renzo and Di Maio (2005) eqn (15)
213 {
214  //K_t0 = 8G_eq * sqrt(R_eq * delta_n) (Di Renzo and Di Maio)
215  tangentialStiffnessZero_ = 8 * shearMod * sqrt(rad * getOverlap());
216 }
217 
218 //k_new
219 //A general function to update the tangential stiffness corresponding to an interaction for *any* loading path
220 //Takes as arguments:
221 // i) the relevant getSlidingFrictionCoefficientStatic() for the current species,
222 // ii) The vector direction corresponding to the particles current tangential motion relative to its initial loading
223 // iii) A boolean to determine whether the turning point force should be considered (false if a 'virgin loading' OR
224 // if the current tangential force exceeds the magnitude of the relevant turning point force
225 // iv) a boolean to indicate whether loading or unloading
226 void MindlinInteraction::updateK_t(const Mdouble fric, const Vec3D direction, const bool useTurningPoint,
227  const bool isLoading)
228 {
229  //creating a simple positive or negative integer 'multiplier', m, whose values are either +1 or -1
230  //that can be used to alter the relevant equations depending on if we are loading or unloading
231  Mdouble m;
232  //A parameter to temporarily hold the relevant turning point so that this function can work both for
233  //loading and loading.
234  Vec3D turningPoint;
235  //Changing the directions and turning points used in calculation depending on if we are loading or unloading
236  if (isLoading)
237  {
238  m = 1;
239  turningPoint = tangentialForceTurningPointUL_;
240  }
241  else
242  {
243  m = -1;
244  turningPoint = tangentialForceTurningPointLU_;
245  }
246  //if the current part of the interaction can be treated as a virgin loading
247  //(i.e. the turning point can be ignored), we simply update K_t in the 'standard' manner
248  if (!useTurningPoint)
249  {
250  //determining the direction of the current tangential force relative to the
251  //current direction of motion
252  Mdouble forceDirection = Vec3D::dot(Vec3D::getUnitVector(tangentialForceTemp_), direction);
253  //finding the scalar ratio f_t1 / (mu f_n1)
254  //(note for a first loading f_t1 is simply equal to f_t, and similarly f_n1 f_n = const
255  Mdouble forceRatio = tangentialForceTemp_.getLength() / (fric * getAbsoluteNormalForce());
256  //changing the direction of the force ratio, if necessary, based on its relative orientation
257  //with respect to the current tangential motion.
258  if (forceDirection < 0)
259  forceRatio = -forceRatio;
260  //using the scalar ratio to update tangential stiffness
261  tangentialStiffness_ = tangentialStiffnessZero_ * cbrt(1 - forceRatio);
262  }
263  //alternatively, if we are reloading and within the limits where the turning point force must be considered
264  else
265  {
266  Vec3D resultantForce = m * (tangentialForceTemp_ - turningPoint);
267  //QUESTION: Does it make sense to use the total force here? Or should I just use the tangentialForceTemp?
268  Mdouble forceDirection = Vec3D::dot(Vec3D::getUnitVector(resultantForce), direction);
269  Mdouble forceRatio = resultantForce.getLength() / (2 * fric * getAbsoluteNormalForce());
270  //changing the direction of the force ratio, if necessary, based on its relative orientation
271  //with respect to the current tangential motion.
272  if (forceDirection < 0)
273  forceRatio = -forceRatio;
274  //using the scalar ratio to update tangential stiffness
275  tangentialStiffness_ = tangentialStiffnessZero_ * cbrt(1 - m * forceRatio);
276  }
277 }
278 
280 {
281  //If tangential forces are absent
282  if (getAbsoluteNormalForce() == 0.0) return;
283 
284  const MindlinSpecies* species = getSpecies();//dynamic_cast
285 
286  //Checking if the relevant particle species has a non-zero friction coefficient
287  if (species->getSlidingFrictionCoefficient() != 0.0)
288  {
289  //Compute the tangential component of relativeVelocity_
290  Vec3D tangentialRelativeVelocity = getRelativeVelocity() - getNormal() * getNormalRelativeVelocity();
291  {
292  // *************************************************************************************************************************
293  // **********************************************GENERAL CALCULATIONS*******************************************************
294  // *************************************************************************************************************************
295  //used to Integrate the spring
296  //if particle-wall
297  if (dynamic_cast<BaseParticle*>(getI()) == nullptr)
298  {
299  slidingSpringVelocity_ = tangentialRelativeVelocity;
300  }
301  //if particle-particle
302  else
303  {
304  slidingSpringVelocity_ = (tangentialRelativeVelocity -
305  Vec3D::dot(slidingSpring_, getP()->getVelocity() - getI()->getVelocity()) *
306  getNormal() / getDistance());
307  }
308  //integrate(getHandler()->timeStep_);
310 
311  //1) Calculating the current value of K_t0
312  //This is identical for both unloading and loading (under constant normal force) and hence can
313  //potentially later be moved for neatness/efficiency
314  Mdouble r1 = 1.0 * getEffectiveRadius();
315  Mdouble shearModulus = species->getEffectiveShearModulus();
316  updateTangentialStiffnessZero(r1, shearModulus);
317 
318  //2) Calculating the relevant tangential dissipation constant based on the current stiffness (and other relevant parameters)
319  Mdouble slidingDissipationCoefficient =
321 
322  //3) Using the parameters determined above, calculating the tangential force as per Di Maio and Di Renzo eqn. (37)
324 
325  //applying the Coulomb criterion (Di Maio and Di Renzo eqn. (36)) to ensure falsely large tangential forces are avoided
326  //if all is OK, i.e. if tangential force is smaller than mu*f_n, we simply add the tangential force calculated.
328  {
329  tangentialForce_ -= slidingDissipationCoefficient * tangentialRelativeVelocity;
331  }
332  //otherwise, we take the force as mu*f_n and recalculate the spring length accordingly
333  else
334  {
337  //adding the force as mu*f_n
339  //updating the spring length accordingly
340  slidingSpring_ = -(tangentialForce_ / (2.0/3.0 * tangentialStiffnessZero_));
341  }
342  }
343  }
345 }
346 
351 {
353 }
354 
359 {
360 
361  //new_edit replacing outdated 'getSlidingStiffness()'
362  //return 0.5 * getSpecies()->getSlidingStiffness() * slidingSpring_.getLengthSquared();
364 }
365 
370 {
372  //k_edit Indeed it should - todo done :)
373  //k_edit
374  //Updating sliding spring to give positive / negative values dependent on the direction of the extension
376 }
377 
382 {
383  return tangentialForce_;
384 }
385 
386 //k_edit
387 //Overwrites the function of the same name in the BaseInteraction class
388 //so that the tangential force direction can be known in the base interaction
389 //class and hence used to correctly output a **signed** tangential force value.
390 //(the value was previously always positive!)
392 {
394 }
395 
400 {
401  return static_cast<const MindlinSpecies*>(getBaseSpecies()->getFrictionForce());
402 ;
403 }
404 
409 {
410  return "Mindlin";
411 }
412 
417 {
421 }
422 
424 {
425  slidingSpring_ = rotationMatrix * slidingSpring_;
427  tangentialForce_ = rotationMatrix * tangentialForce_;
428 }
429 
430 //k_edit
439 {
441 }
442 
443 //k_edit
450 {
451  absoluteNormalForcePrevious_ = absoluteNormalForcePrevious;
452 }
Mdouble getSlidingDissipation() const
Allows the tangential viscosity to be accessed.
Mdouble getEffectiveRadius() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
void updateTangentialStiffnessZero(Mdouble rad, double shearMod)
Mdouble getSlidingFrictionCoefficient() const
Allows the (dynamic) Coulomb friction coefficient to be accessed.
std::string getBaseName() const
Returns the type/name of interaction (sliding friction interaction)
Mdouble getEffectiveShearModulus() const
Allows the shear modulus to be accessed.
static Vec3D getUnitVector(const Vec3D &a)
Returns a unit Vec3D based on a.
Definition: Vector.cc:345
InteractionHandler * getHandler() const
Gets a point to the interaction handlers to which this interaction belongs.
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getTangentialOverlap() const override
Returns the amount of tangential overlap which is needed by BaseInteraction::writeToFstat().
Vec3D slidingSpringVelocity_
Stores the rate at which the sliding spring compressed or relaxed. Set in the member function compute...
Vec3D tangentialForceTurningPointULTemp_
Vec3D tangentialDisplacementTurningPointLU_
Mdouble absoluteNormalForcePrevious_
const Vec3D & getRelativeVelocity() const
Returns a constant reference to a vector of relative velocity.
Computes the forces corresponding to sliding friction.
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
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 slidingSpring_
Stores the amount of sliding spring ( ) compression from the expression . Set in the member function ...
MindlinSpecies contains the parameters used to describe sliding friction.
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 ...
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.
const Mdouble getTangentialForceDirection() const
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_.
Mdouble getNormalRelativeVelocity() const
Returns a double which is the norm (length) of the relative velocity vector.
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 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.
const Vec3D getTangentialForce() const override
Returns the sliding friction force vector.
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
void addForce(Vec3D force)
add an force increment to the total force.
Vec3D tangentialForce_
Computes the tangential force such that . Set and computed in computeFrictionForce().
void setTangentialStiffnessZero(Mdouble newKt0)
Mdouble getAbsoluteNormalForcePrevious() const
Returns the absolute value of the norm (length) of the previous Normal force vector.
Mdouble getDistance() const
Returns an Mdouble which is the norm (length) of distance vector.
const MindlinSpecies * getSpecies() const
Returns a const pointer of type MindlinSpecies*.
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.
Mdouble getTangentialStiffnessZero()
Mdouble tangentialStiffnessZeroPrevious_
Implementation of a 3D matrix.
Definition: Matrix.h:37
Definition: Vector.h:49
void updateK_t(Mdouble fric, Vec3D direction, bool useTurningPoint, bool isLoading)
Vec3D tangentialDisplacementTurningPointUL_
MindlinInteraction()
Empty constructor.
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 tangentialForceTurningPointLUTemp_
void setAbsoluteNormalForcePrevious(Mdouble absoluteNormalForcePrevious)
allows the previous normal force to be (re)set from external classes
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 getEffectiveMass() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
~MindlinInteraction() override
Destructor.
Mdouble getAbsoluteNormalForce() const
Returns the absolute value of the norm (length) of the Normal force vector.