MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SinterInteraction.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 "SinterInteraction.h"
28 #include "InteractionHandler.h"
29 #include "Particles/BaseParticle.h"
30 #include "DPMBase.h"
33 
40  : BaseInteraction(P, I, timeStamp)
41 {
42  plasticOverlap_ = 0;
43 #ifdef DEBUG_CONSTRUCTOR
44  std::cout<<"SinterInteraction::SinterInteraction() finished"<<std::endl;
45 #endif
46 }
47 
52  : BaseInteraction(p)
53 {
55 #ifdef DEBUG_CONSTRUCTOR
56  std::cout<<"SinterInteraction::SinterInteraction(const SinterInteraction &p finished"<<std::endl;
57 #endif
58 }
59 
61 = default;
62 
67 {
68 #ifdef DEBUG_DESTRUCTOR
69  std::cout<<"SinterInteraction::~SinterInteraction() finished"<<std::endl;
70 #endif
71 }
72 
77 void SinterInteraction::write(std::ostream& os) const
78 {
80  os << " plasticOverlap " << plasticOverlap_;
81 }
82 
87 void SinterInteraction::read(std::istream& is)
88 {
90  std::string dummy;
91  is >> dummy >> plasticOverlap_;
92 }
93 
97 std::string SinterInteraction::getBaseName() const
98 {
99  return "Sinter";
100 }
101 
106 {
107  // Compute the relative velocity vector of particle P w.r.t. I
109  getP()->getVelocityAtContact(getContactPoint()) - getI()->getVelocityAtContact(getContactPoint()));
110  // Compute the projection of vrel onto the normal (can be negative)
112 
113  if (getOverlap() > 0) //if contact forces
114  {
115  const SinterNormalSpecies* species = getSpecies();
116 
117  // calculate the effective diameter, equal to the radius for two equal-sized particles
118  const Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
119 
120  // [1] Calculate the overlap above which the max. unloading stiffness
121  //becomes active (the 'fluid branch'). In this interaction, it is the same than delta*
122 
123  const Mdouble d_fluid_0 = (species->getUnloadingStiffnessMax()
124  / (species->getUnloadingStiffnessMax() - species->getLoadingStiffness()))
125  * species->getPenetrationDepthMax() * effectiveDiameter;
126 
127  //const Mdouble deltaStar = species->getPenetrationDepthMax() * effectiveDiameter;
128 
129  //[2] Compute the rate d(k2/k1)/d(delta0). To obtain this parameter, the linear relationship
130  //between unloading stiffness and maximum overlap is used.
131  const Mdouble dk = (species->getUnloadingStiffnessMax() / species->getLoadingStiffness() - 1.0) / d_fluid_0;
132 
133  //increase max overlap if necessary
134  //Here, two relationships are used. The linear relationship between unloadingstiffness max and
135  //the equilibrium overlap based on the unloading stiffness.
136  //k2 = k1*(1+dk*d), k1*d = k2*(d-d0)
137  //k1*d = k1*(1 + dk*d)*(d-d0)
138  //d = d - d0 +dk*d^2 -dk*d*d0
139  //dk*d^2 = dk *d * d0 + d0
140  //dk*d^2 = (dk*d + 1)*d0
141  //dkd^2/(dk*d + 1) = d0
142  //d0 = d/(1 + 1/(dk*d))
143 
144  //[3] Compute the equilibirum overlap is:
145  const Mdouble d0 = getOverlap()/(1.0 + 1.0/(dk*getOverlap()));
146 
147  const Mdouble minPlasticOverlap = std::min(d0,d_fluid_0);
148 
149  //[4] Determine the plastic overlap
150  plasticOverlap_ = std::max(minPlasticOverlap, plasticOverlap_);
151 
152  //[5] Compute the unloading Stiffness \hat{k2}.
153  const Mdouble unloadingStiffness = species->getLoadingStiffness() * (1.0 + dk * plasticOverlap_);
154 
155  //[6] Compute the elastic force
156  Mdouble normalForce = unloadingStiffness * (getOverlap() - plasticOverlap_);
157 
158  //[7] Add the adhesive force (distinct from sintering)
159  Mdouble nonSinterAdhesiveForce = -species->getCohesionStiffness() * getOverlap();
160 
161  if (normalForce < nonSinterAdhesiveForce)
162  {
163  plasticOverlap_ = (1.0 + species->getCohesionStiffness() / unloadingStiffness) * getOverlap();
164  normalForce = nonSinterAdhesiveForce;
165  }
166 
167  //[[8] Add dissipative force (distinct from sintering)
168  normalForce -= species->getDissipation() * getNormalRelativeVelocity();
169 
170  //[9] Sintering effect:
171  Mdouble adhesiveForce = species->getSinterAdhesion() * effectiveDiameter;
172 
173  //[10] Now set the interaction force equal to this normal force (friction and adhesive forces will be added later)
174  setForce(getNormal() * ((normalForce - adhesiveForce)));
175  setTorque(Vec3D(0.0, 0.0, 0.0));
176  //used for tangential force calculations; don't add adhesive force components
177  setAbsoluteNormalForce(std::abs(normalForce));
178 
179  //[11] Approaches - Increase plastic overlap due to sintering
180  const Mdouble baseNum = 9*constants::pi*species->getComplianceZero()*species->getSurfTension()/(effectiveDiameter);
181  const Mdouble a0_R = std::pow(baseNum,1.0/3.0);
182  const Mdouble realOverlap = std::sqrt(getOverlap()/effectiveDiameter);
183 
184  DPMBase* dpmBase = getHandler()->getDPMBase();
185  Mdouble rateOverlap;
186  // sinter adhesion force fa=sinterAdhesion_*radius in sinter rate:
188  {
189  rateOverlap = normalForce * species->getSinterRate() /
190  (0.375 * species->getSinterAdhesion() * mathsFunc::square(
191  getOverlap() / effectiveDiameter));
192  if (species->getSinterRate() == 0) rateOverlap = 0;
193  }
194  else if (species->getSinterType() == SINTERTYPE::CONSTANT_RATE)
195  {
196  rateOverlap = 2.0 * normalForce * species->getSinterRate() / species->getSinterAdhesion();
197  if (species->getSinterRate() == 0) rateOverlap = 0;
198  }
200  {
201  ThermalParticle* tp = dynamic_cast<ThermalParticle*>(getP());
202  ThermalParticle* ti = dynamic_cast<ThermalParticle*>(getI());
203  logger.assert(tp && ti,
204  "warning contact partners have to be ThermalParticle's if this sinter species is used");
205  double temperature =
206  2.0 * tp->getTemperature() * ti->getTemperature() / (tp->getTemperature() + ti->getTemperature());
207  rateOverlap = 2.0 * normalForce * species->getTemperatureDependentSinterRate(temperature) /
208  species->getSinterAdhesion();
209  }
210  else if (species->getSinterType() == SINTERTYPE::REGIME_SINTERING)
211  {
212  if (realOverlap < a0_R) {
213  //Here, Sintering rate has unit of [1/s]
214  rateOverlap = normalForce / species->getSinterAdhesion();
215  //rateOverlap = effectiveDiameter* species->getSinterRate();
216 
217  }else {
218  const Mdouble val0 = std::pow((63 * pow(constants::pi, 3.0) / 16.0), 1.0 / 7.0);
219  const Mdouble val1 = pow((species->getSeparationDis()*10 / effectiveDiameter), 2.0 / 7.0);
220  const Mdouble val2 = (2.0 / 7.0) * (2.0 * species->getConstantC1() * species->getSurfTension() /
221  (effectiveDiameter));
222  const Mdouble val3 = 1.0 / (std::pow(getOverlap(), 5.0 / 2.0));
223  const Mdouble val4 = std::pow(effectiveDiameter, 7.0 / 2.0);
224 
225  rateOverlap = (std::pow(val0 * val1, 7.0) * val2 * val3 * val4) * normalForce / species->getSinterAdhesion();
226 
227  const Mdouble aVisc_R = std::pow(63.0 * std::pow(constants::pi, 3.0), 1.0 / 5.0) *
228  std::pow((1.0 / 8.0) * (species->getSeparationDis()/33.0) / (effectiveDiameter), 2.0 / 5.0);
229 
230  if (realOverlap > aVisc_R) {
231  rateOverlap = 2.0*normalForce * species->getSinterRate() / species->getSinterAdhesion();
232 // logger(INFO," RealOver % aVisc_R %",realOverlap,aVisc_R);
233  }
234  }
235  }
236  else
237  {
238  rateOverlap = 0;
239  //missing: add the sintering model 'modified Frenkel' of the Pokula paper
240  }
241  plasticOverlap_ = std::max(0.0, std::min(d_fluid_0, plasticOverlap_ + rateOverlap * dpmBase->getTimeStep()));
242 
243  /*//change particle radius by dr
244  Mdouble dr;
245  BaseParticle* PParticle = dynamic_cast<BaseParticle*>(getP());
246  BaseParticle* IParticle = dynamic_cast<BaseParticle*>(getI());
247  if (dpmBase->getSystemDimensions()==2) {
248  //2D: increase the radius of each particle such that the particle area
249  //increases by the same amount that the contact area decreases
250  //Particle circumference C = 2 pi r increased by dr => dA = 2 pi r dr
251  //Contact line L = 2*sqrt(2*r*o) indented by do/2 => dA = sqrt(2*r*o) do
252  //Thus, dr = sqrt(0.5*o/r)/pi do.
253  dr = sqrt(0.5*plasticOverlap_/effectiveDiameter)/3.14 *doverlap;
254  } else {
255  //3D: increase the radius of each sphere such that the particle volume
256  //increases by the same amount that the contact volume decreases
257  //Particle surface area S = 4 pi r^2 increased by dr => dA = 4 pi r^2 dr
258  //Contact area L = pi 2*r*o indented by do/2 => dA = pi r o do
259  //Thus, dr = 0.25*o/r do
260  dr = 0.25*plasticOverlap_/effectiveDiameter *doverlap;
261  }
262  if (PParticle==nullptr) { //if P is a wall
263  IParticle->setRadius(IParticle->getRadius()+dr);//should be twice that amount
264  } else if (IParticle==nullptr) { //if I is a wall
265  PParticle->setRadius(PParticle->getRadius()+dr);
266  } else { //if both P and I are particles
267  PParticle->setRadius(PParticle->getRadius()+dr);
268  IParticle->setRadius(IParticle->getRadius()+dr);
269  }*/
270  }
271  else
272  {
274  setForce(Vec3D(0.0, 0.0, 0.0));
275  setTorque(Vec3D(0.0, 0.0, 0.0));
276  }
277 }
278 
283 {
284  if (getOverlap() > 0)
286  else
287  return 0.0;
289 }
290 
295 {
296  return static_cast<const SinterNormalSpecies*>(getBaseSpecies()->getNormalForce());
297 ;
298 }
299 
304 {
305  return plasticOverlap_;
306 }
307 
312 {
313  plasticOverlap_ = plasticOverlap;
314 }
315 
320 {
321  const SinterNormalSpecies* species = getSpecies();
322  Mdouble effectiveDiameter = 2.0 * getEffectiveRadius();
323  Mdouble deltaMaxFluid = species->getPenetrationDepthMax() * effectiveDiameter /
324  (1.0 - species->getLoadingStiffness() / species->getUnloadingStiffnessMax());
325  if (getOverlap() > deltaMaxFluid)
326  return species->getUnloadingStiffnessMax();
327  else
328  return species->getLoadingStiffness() + (species->getUnloadingStiffnessMax() - species->getLoadingStiffness()) *
329  getPlasticOverlap() / deltaMaxFluid;
330 }
331 
333 {
334  return 2;
335 }
336 
337 std::string SinterInteraction::getTypeVTK(unsigned i) const
338 {
339  return "Float32";
340 }
341 
342 std::string SinterInteraction::getNameVTK(unsigned i) const
343 {
344  if (i == 0)
345  return "plasticOverlap";
346  else
347  return "neckRadius";
348 }
349 
350 std::vector<Mdouble> SinterInteraction::getFieldVTK(unsigned i) const
351 {
352  if (i == 0)
353  return std::vector<Mdouble>(1, plasticOverlap_);
354  else
355  return std::vector<Mdouble>(1, sqrt(2.0 * getEffectiveRadius() * plasticOverlap_));
356 }
357 
358 
359 //set xlabel 't'; set ylabel 'f_n = f_ep-f_a'; p 'SingleParticleCT.fstat' u 1:9 title 'f_g=f_a', 'SingleParticleCT_noGravity.fstat' u 1:9 title 'f_g=0'/
360 
361 //set xlabel 't'; set ylabel 'delta/r'; p 'SingleParticleCT.fstat' u 1:($7/1e-6) title 'f_g=f_a', 'SingleParticleCT_noGravity.fstat' u 1:($7/1e-6) title 'f_g=0'
Mdouble getEffectiveRadius() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:72
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
const Vec3D & getRelativeVelocity() const
Returns a constant reference to a vector of relative velocity.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Mdouble getSinterAdhesion() const
Accesses sinterAdhesion_.
void setRelativeVelocity(Vec3D relativeVelocity)
set the relative velocity of the current of the interactions.
std::vector< Mdouble > getFieldVTK(unsigned i) const override
Mdouble getTemperature() const
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.
Mdouble getUnloadingStiffness() const
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.
SINTERTYPE getSinterType() const
void setPlasticOverlap(Mdouble plasticOverlap)
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.
~SinterInteraction() override
Destructor.
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble getComplianceZero() const
std::string getTypeVTK(unsigned i) const override
Mdouble getLoadingStiffness() const
Returns the loading stiffness of the linear plastic-viscoelastic normal force.
void computeNormalForce()
Creates a copy of an object of this class. (Deep copy)
Mdouble getCohesionStiffness() const
Returns the cohesive stiffness of the linear plastic-viscoelastic normal force.
void setTorque(Vec3D torque)
set the total force (this is used by the normal force, tangential torques are added use addTorque) ...
Mdouble getSeparationDis() const
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
std::string getNameVTK(unsigned i) const override
void setAbsoluteNormalForce(Mdouble absoluteNormalForce)
the absolute values of the norm (length) of the normal force
Mdouble getPenetrationDepthMax() const
Returns the maximum penetration depth of the linear plastic-viscoelastic normal force.
SinterNormalSpecies contains the parameters used to describe a plastic-cohesive normal force (Stefan ...
std::function< double(double temperature)> getTemperatureDependentSinterRate() const
unsigned getNumberOfFieldsVTK() const override
Mdouble getPlasticOverlap() const
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
virtual std::string getBaseName() const
Returns the name of the interaction.
Defines the basic properties that a interactable object can have.
Mdouble getUnloadingStiffnessMax() const
Returns the maximum unloading stiffness of the linear plastic-viscoelastic normal force...
Mdouble getDissipation() const
Allows the normal dissipation to be accessed.
void write(std::ostream &os) const override
Interaction write function, which accepts an std::ostream as input.
Definition: Vector.h:49
Mdouble getSurfTension() const
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.
Computes normal forces in case of a linear plastic visco-elastic interaction.
Mdouble getConstantC1() const
const SinterNormalSpecies * getSpecies() const
Mdouble getElasticEnergy() const override
Computes and returns the amount of elastic energy stored in the spring.
Mdouble getSinterRate() const
Accesses sinterRate_.