MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HertzianSinterInteraction.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 
28 #include "InteractionHandler.h"
29 #include "Particles/BaseParticle.h"
30 #include "DPMBase.h"
32 #include "Math/ExtendedMath.h"
33 #include <iomanip>
34 #include <fstream>
35 #include <cmath> // std::max
36 
43  : BaseInteraction(P, I, timeStamp)
44 {
45  maxOverlap_ = 0;
46 #ifdef DEBUG_CONSTRUCTOR
47  std::cout<<"HertzianSinterInteraction::HertzianSinterInteraction() finished"<<std::endl;
48 #endif
49 }
50 
55  : BaseInteraction(p)
56 {
58 #ifdef DEBUG_CONSTRUCTOR
59  std::cout<<"HertzianSinterInteraction::HertzianSinterInteraction(const HertzianSinterInteraction &p finished"<<std::endl;
60 #endif
61 }
62 
64 
69 {
70 #ifdef DEBUG_DESTRUCTOR
71  std::cout<<"HertzianSinterInteraction::~HertzianSinterInteraction() finished"<<std::endl;
72 #endif
73 }
74 
79 void HertzianSinterInteraction::write(std::ostream& os) const
80 {
82  os << " maxOverlap " << maxOverlap_;
83 }
84 
89 void HertzianSinterInteraction::read(std::istream& is)
90 {
92  std::string dummy;
93  is >> dummy >> maxOverlap_;
94 }
95 
100 {
101  return "HertzianSinter";
102 }
103 
108 {
109  // Compute the relative velocity vector of particle P w.r.t. I
111  getP()->getVelocityAtContact(getContactPoint()) - getI()->getVelocityAtContact(getContactPoint()));
112  // Compute the projection of vrel onto the normal (can be negative)
114 
115  if (getOverlap() > 0) //if contact forces
116  {
117  const HertzianSinterNormalSpecies* species = getSpecies();
118  Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
119 
120  //calculate the overlap above which the max. unloading stiffness becomes active (the 'fluid branch')
121  static Mdouble maxFactor = 1
123  cbrt((species->getLoadingModulus() + species->getCohesionModulus()) /
124  species->getUnloadingModulusMax()));
125  Mdouble deltaStar = species->getPenetrationDepthMax() * effectiveDiameter / maxFactor;
126 
127  //increase max overlap if necessary
128  if (getOverlap() > getMaxOverlap())
129  {
131  std::cout << "," << getHandler()->getDPMBase()->getTime();
132  }
133  //limit max overlap if necessary
134  if (getMaxOverlap() > deltaStar)
135  setMaxOverlap(deltaStar);
136 
137  //calculate the unloading modulus
138  Mdouble loadingCohesionModulus = species->getLoadingModulus() + species->getCohesionModulus();
139  Mdouble unloadingModulus = loadingCohesionModulus
140  + (species->getUnloadingModulusMax() - loadingCohesionModulus) *
141  (getMaxOverlap() / deltaStar);
142 
143  //calculate the overlap where the force is minimal
144  Mdouble factor = 1 - mathsFunc::square(cbrt(loadingCohesionModulus / unloadingModulus));
145  Mdouble minOverlap = factor * maxOverlap_;
146 
147  //add dissipative force
148  Mdouble normalForce = -species->getDissipation() * getNormalRelativeVelocity();
149 
150  //compute elastic force
151  if (getOverlap() < minOverlap)
152  {
153  //decrease max overlap if in cohesive range
154  std::cout << "." << getHandler()->getDPMBase()->getTime();
155  setMaxOverlap(getOverlap() / factor);
156  }
157  else
158  {
159  Mdouble contactRadius = sqrt(2.0 * effectiveDiameter * (getOverlap() - minOverlap));
160  normalForce += 4. / 3. * unloadingModulus * contactRadius * (getOverlap() - minOverlap);
161  }
162 
163  setAbsoluteNormalForce(std::abs(normalForce)); //used for the friction force calculations;
164 
165  Mdouble contactRadius = sqrt(2.0 * effectiveDiameter * getOverlap());
166  normalForce -= 4. / 3. * species->getCohesionModulus() * contactRadius * getOverlap();
167 
168  setForce(getNormal() * normalForce);
169  setTorque(Vec3D(0.0, 0.0, 0.0));
170 
171  //now add the sintering model 'modified Frenkel' of the Pokula paper
172  //plasticOverlap_+=species->getSinterRate()*(deltaStar-plasticOverlap_)*getHandler()->getDPMBase()->getTimeStep();
173  //x/a=sqrt(2*a*del)/a
174  Mdouble x = 1e-10 + sqrt(2.0 * maxOverlap_ / effectiveDiameter);
175  //Mdouble x2 = x*x;
176  Mdouble dx = 0.5 /
177  x;//+ x*(-0.5 + x2* (0.15625 + x2*(-0.0208333 +x2*(-0.00325521 +x2*(0.000189887 +x2*0.0000542535)))));
178  Mdouble doverlap = x * dx * effectiveDiameter;
179  //doverlap = 0.5/(factor*factor*plasticOverlap_);
180  maxOverlap_ += species->getSinterRate() * doverlap * getHandler()->getDPMBase()->getTimeStep();
181  }
182  else
183  {
185  setForce(Vec3D(0.0, 0.0, 0.0));
186  setTorque(Vec3D(0.0, 0.0, 0.0));
187  }
188 }
189 
194 {
196 }
197 
202 {
203  if (getOverlap() > 0)
204  return 0.5 * (getSpecies()->getLoadingModulus() * mathsFunc::square(getOverlap()));
205  else
206  return 0.0;
208 }
209 
214 {
215  return static_cast<const HertzianSinterNormalSpecies*>(getBaseSpecies()->getNormalForce());
216 }
217 
222 {
223  return maxOverlap_;
224 }
225 
230 {
231  maxOverlap_ = maxOverlap;
232 }
233 
238 {
239  const HertzianSinterNormalSpecies* species = getSpecies();
240  Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
241  Mdouble deltaMaxFluid = species->getPenetrationDepthMax() * effectiveDiameter /
242  (1.0 - species->getLoadingModulus() / species->getUnloadingModulusMax());
243  if (getOverlap() > deltaMaxFluid)
244  return species->getUnloadingModulusMax();
245  else
246  return species->getLoadingModulus() +
247  (species->getUnloadingModulusMax() - species->getLoadingModulus()) * getMaxOverlap() / deltaMaxFluid;
248 }
249 
void write(std::ostream &os) const override
Interaction write function, which accepts an std::ostream as input.
const HertzianSinterNormalSpecies * getSpecies() const
Mdouble getEffectiveRadius() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
Mdouble getLoadingModulus() const
Returns the loading stiffness of the linear plastic-viscoelastic normal force.
InteractionHandler * getHandler() const
Gets a point to the interaction handlers to which this interaction belongs.
double Mdouble
Definition: GeneralDefine.h:34
virtual std::string getBaseName() const
Returns the name of the interaction.
Mdouble getDissipation() const
Allows the normal dissipation to be accessed.
const Vec3D & getRelativeVelocity() const
Returns a constant reference to a vector of relative velocity.
void setRelativeVelocity(Vec3D relativeVelocity)
set the relative velocity of the current of the interactions.
BaseNormalForce * getNormalForce() const
Definition: BaseSpecies.h:148
void setForce(Vec3D force)
set total force (this is used by the normal force, tangential forces are added use addForce) ...
const Vec3D & getContactPoint() const
Gets constant reference to contact point (vector).
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
void setNormalRelativeVelocity(Mdouble normalRelativeVelocit)
set the normal component of the relative velocity.
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.
Mdouble getSinterRate() const
Allows the normal dissipation to be accessed.
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
Mdouble getNormalRelativeVelocity() const
Returns a double which is the norm (length) of the relative velocity vector.
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.
Mdouble getCohesionModulus() const
Returns the cohesive stiffness of the linear plastic-viscoelastic normal force.
void computeNormalForce()
Calls computeSinterForce().
void setTorque(Vec3D torque)
set the total force (this is used by the normal force, tangential torques are added use addTorque) ...
HertzianSinterNormalSpecies contains the parameters used to describe a plastic-cohesive normal force ...
Computes normal forces in case of a linear plastic visco-elastic interaction.
void computeSinterForce()
Creates a copy of an object of this class. (Deep copy)
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
void setAbsoluteNormalForce(Mdouble absoluteNormalForce)
the absolute values of the norm (length) of the normal force
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
Defines the basic properties that a interactable object can have.
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 print function, which accepts an std::ostream as input.
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:797
Mdouble getUnloadingModulusMax() const
Returns the maximum unloading stiffness of the linear plastic-viscoelastic normal force...
Mdouble getPenetrationDepthMax() const
Returns the maximum penetration depth of the linear plastic-viscoelastic normal force.
~HertzianSinterInteraction() override
Destructor.
Mdouble getElasticEnergy() const override
Computes and returns the amount of elastic energy stored in the spring.