MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SinterNormalSpecies.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 "SinterNormalSpecies.h"
29 #include "Logger.h"
30 #include "Species/BaseSpecies.h"
31 
32 class BaseParticle;
33 
34 class BaseInteractable;
35 
37  : BaseNormalForce()
38 {
40  loadingStiffness_ = 0.0;
42  cohesionStiffness_ = 0.0;
44  dissipation_ = 0.0;
45  sinterAdhesion_ = 0.0;
46  sinterRate_ = 0.0;
48  complianceZero_ = 0.0;
49  surfTension_ = 0.0;
50  constantC1_ = 0.0;
51  separationDis_ = 0.0;
52 #ifdef DEBUG_CONSTRUCTOR
53  std::cout<<"SinterNormalSpecies::SinterNormalSpecies() finished"<<std::endl;
54 #endif
55 }
56 
61  : BaseNormalForce(p)
62 {
76 #ifdef DEBUG_CONSTRUCTOR
77  std::cout<<"SinterNormalSpecies::SinterNormalSpecies(const SinterNormalSpecies &p) finished"<<std::endl;
78 #endif
79 }
80 
82 {
83 #ifdef DEBUG_DESTRUCTOR
84  std::cout<<"SinterNormalSpecies::~SinterNormalSpecies() finished"<<std::endl;
85 #endif
86 }
87 
91 void SinterNormalSpecies::write(std::ostream& os) const
92 {
93  os << " loadingStiffness " << loadingStiffness_;
94  os << " maxUnloadingStiffness " << unloadingStiffnessMax_;
95  os << " cohesionStiffness " << cohesionStiffness_;
96  os << " maxPenetration " << penetrationDepthMax_;
97  os << " dissipation " << dissipation_;
98  os << " sinterAdhesion " << sinterAdhesion_;
99  os << " inverseSinterViscosity " << inverseSinterViscosity_;
100  os << " sinterRate " << sinterRate_;
101  os << " complianceZero " << complianceZero_;
102  os << " surfTension " << surfTension_;
103  os << " constantC1 " << constantC1_;
104  os << " separationDis " << separationDis_;
105  os << " sinterType " << (unsigned) sinterType_;
106 }
107 
111 void SinterNormalSpecies::read(std::istream& is)
112 {
113  std::string dummy;
114  is >> dummy >> loadingStiffness_;
115  is >> dummy >> unloadingStiffnessMax_;
116  is >> dummy >> cohesionStiffness_;
117  is >> dummy >> penetrationDepthMax_;
118  is >> dummy >> dissipation_;
119  is >> dummy >> sinterAdhesion_;
120  is >> dummy >> inverseSinterViscosity_;
121  is >> dummy >> sinterRate_;
122  is >> dummy >> complianceZero_;
123  is >> dummy >> surfTension_;
124  is >> dummy >> constantC1_;
125  is >> dummy >> separationDis_;
126  unsigned type;
127  is >> dummy >> type;
128  sinterType_ = (SINTERTYPE) type;
129 }
130 
135 {
136  return "Sinter";
137 }
138 
145 {
157 
158 }
159 
166 void SinterNormalSpecies::setPlasticParameters(Mdouble loadingStiffness, Mdouble unloadingStiffnessMax,
167  Mdouble cohesionStiffness, Mdouble penetrationDepthMax)
168 {
169  if (loadingStiffness <= 0 || unloadingStiffnessMax < loadingStiffness || cohesionStiffness < 0 ||
170  penetrationDepthMax < 0)
171  {
172  std::cerr << "Error: arguments of setPlasticParameters do not make sense" << std::endl;
173  exit(-1);
174  }
175  setLoadingStiffness(loadingStiffness);
176  setUnloadingStiffnessMax(unloadingStiffnessMax);
177  setCohesionStiffness(cohesionStiffness);
178  setPenetrationDepthMax(penetrationDepthMax);
179 }
180 
185 {
186  return loadingStiffness_;
187 }
188 
193 {
194  return unloadingStiffnessMax_;
195 }
196 
201 {
202  return cohesionStiffness_;
203 }
204 
209 {
210  return penetrationDepthMax_;
211 }
212 
217 {
218  loadingStiffness_ = loadingStiffness;
219 }
220 
225 {
226  unloadingStiffnessMax_ = unloadingStiffnessMax;
227 }
228 
233 {
234  cohesionStiffness_ = cohesionStiffness;
235 }
236 
241 {
242  penetrationDepthMax_ = penetrationDepthMax;
243 }
244 
250 {
251  if (unloadingStiffnessMax_ / (.5 * mass) < mathsFunc::square(dissipation_ / mass))
252  {
253  std::cerr << "Dissipation too high; max. allowed " << sqrt(2.0 * unloadingStiffnessMax_ * mass) << std::endl;
254  exit(-1);
255  }
256  return 0.02 * constants::pi /
257  std::sqrt(unloadingStiffnessMax_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass));
258 }
259 
265 {
266  if (dissipation < 0)
267  {
268  logger(ERROR, "setDissipation(%)", dissipation);
269  exit(-1);
270  }
271  dissipation_ = dissipation;
272 }
273 
279 {
280  if (sinterAdhesion < 0)
281  {
282  logger(ERROR, "setSinterAdhesion(%)", sinterAdhesion);
283  exit(-1);
284  }
285  sinterAdhesion_ = sinterAdhesion;
286 }
287 
293 {
294  if (sinterRate < 0)
295  {
296  logger(ERROR, "setSinterRate(%)", sinterRate);
297  }
298  sinterRate_ = sinterRate;
299 }
300 
302 {
303  if (complianceZero < 0)
304  {
305  logger(ERROR, "complianceZero(%)", complianceZero);
306  }
307  complianceZero_ = complianceZero;
308 }
309 
311 {
312  if (surfTension < 0)
313  {
314  logger(ERROR, "SurfTension(%)", surfTension);
315  }
316  surfTension_ = surfTension;
317 }
318 
320 {
321  if (complianceOne < 0)
322  {
323  logger(ERROR, "ComplianceOne(%)", complianceOne);
324  }
325  constantC1_ = complianceOne;
326 }
327 
329 {
330  if (separationDis < 0)
331  {
332  logger(ERROR, "SeparationDis(%)", separationDis);
333  }
334  separationDis_ = separationDis;
335 }
336 
337 
338 
339 
345 {
346  sinterType_ = sinterType;
348  logger(INFO, "Sintertype set to CONSTANT_RATE");
350  logger(INFO, "Sintertype set to PARHAMI_MCKEEPING");
352  logger(INFO, "Sintertype set to TEMPERATURE_DEPENDENT_FRENKEL");
354  logger(INFO, "Sintertype set to REGIME_SINTERING");
355  else
356  logger(ERROR, "Sintertype not understood");
357 }
358 
364 {
365  //assertOrDie(inverseSinterViscosity < 0, "inverseSinterViscosity should be non-negative!");
366  if (inverseSinterViscosity < 0)
367  {
368  logger(ERROR, "setInverseSinterViscosity(%)", inverseSinterViscosity);
369  }
371  inverseSinterViscosity_ = inverseSinterViscosity;
372 }
373 
374 void SinterNormalSpecies::setSinterForceAndTime(Mdouble adhesionForce, Mdouble sinterTime, Mdouble radius)
375 {
378  setSinterAdhesion(adhesionForce / radius);
379  setInverseSinterViscosity(1.0 / (93.75 * adhesionForce * sinterTime / std::pow(radius, 5)));
380  logger(INFO, "Set sintering parameters: adhesion force f_a=%*r, sinter viscosity is nu = contactRadius^4/%",
382 }
383 
389  (Mdouble alpha, Mdouble beta, Mdouble atomicVolume /*Omega*/, Mdouble surfaceEnergy /*gamma_s*/,
390  Mdouble thicknessDiffusion /*deltaB*D0B*/, Mdouble activationEnergy /*QB*/, Mdouble temperature /*T*/)
391 {
392  setSinterType(SINTERTYPE::PARHAMI_MCKEEPING);
393  const Mdouble boltzmannConstant /*k_B*/ = 1.38064852e-23;
394  const Mdouble gasConstant /*R_g*/ = 8.314459848;
395  const Mdouble thicknessDiffusionVacancy /*DB*/ =
396  thicknessDiffusion * exp(-activationEnergy / gasConstant / temperature);
397  const Mdouble diffusionParameter /*DeltaB*/ =
398  atomicVolume / boltzmannConstant / temperature * thicknessDiffusionVacancy;
399  setInverseSinterViscosity(constants::pi / (2.0 * beta * diffusionParameter));
400  setSinterAdhesion(alpha / beta * constants::pi * surfaceEnergy);
401  logger(INFO,
402  "Set sintering parameters: the adhesion force f_a=%*r, the rate of change of the plastic overlap is d(delta0)/dt=f_ep/(%*a^4)",
403  getSinterAdhesion(),
404  getInverseSinterViscosity());
405 }
406 
411 {
412  setLoadingStiffness(new_.k);
413  setDissipation(new_.disp);
414 }
415 
420 {
421  return dissipation_;
422 }
423 
428 {
429  return sinterAdhesion_;
430 }
431 
436 {
438 }
439 
444 {
445  return sinterRate_;
446 }
447 
449 {
450  return complianceZero_;
451 }
452 
454 {
455  return surfTension_;
456 }
457 
459 {
460  return constantC1_;
461 }
462 
464 {
465  return separationDis_;
466 }
467 
469 {
470  return sinterType_;
471 }
472 
474 {
475  return temperatureDependentSinterRate_(temperature);
476 }
477 
478 std::function<double(double temperature)> SinterNormalSpecies::getTemperatureDependentSinterRate() const
479 {
481 }
482 
484  std::function<double(double temperature)> temperatureDependentSinterRate)
485 {
486  temperatureDependentSinterRate_ = temperatureDependentSinterRate;
487 }
488 
489 
498 {
499  if (eps == 0.0)
500  {
502  dissipation_ = std::sqrt(2.0 * mass * loadingStiffness_);
503  }
504  else
505  {
506  dissipation_ = -mass / tc * std::log(eps);
508  .5 * mass * (mathsFunc::square(constants::pi / tc) + mathsFunc::square(dissipation_ / mass));
509  }
511 }
512 
520 {
521  loadingStiffness_ = stiffness;
522  if (eps == 0.0)
523  {
524  dissipation_ = std::sqrt(2.0 * mass * stiffness);
525  }
526  else
527  {
528  dissipation_ =
529  -std::sqrt(2.0 * mass * stiffness / (constants::sqr_pi + mathsFunc::square(log(eps)))) * log(eps);
530  }
532 }
533 
536 {
537  if (mass <= 0)
538  {
539  std::cerr << "Error in getCollisionTime(" << mass
540  << ") mass is not set or has an unexpected value, (getCollisionTime(" << mass << "))" << std::endl;
541  exit(-1);
542  }
543  if (loadingStiffness_ <= 0)
544  {
545  std::cerr << "Error in getCollisionTime(" << mass << ") stiffness=" << loadingStiffness_
546  << " is not set or has an unexpected value, (getCollisionTime(" << mass << "), with stiffness="
547  << loadingStiffness_ << ")" << std::endl;
548  exit(-1);
549  }
550  if (dissipation_ < 0)
551  {
552  std::cerr << "Error in getCollisionTime(" << mass << ") dissipation=" << dissipation_
553  << " is not set or has an unexpected value, (getCollisionTime(" << mass << "), with dissipation="
554  << dissipation_ << ")" << std::endl;
555  exit(-1);
556  }
557  Mdouble tosqrt = loadingStiffness_ / (.5 * mass) - mathsFunc::square(dissipation_ / mass);
558  if (tosqrt <= -1e-8 * loadingStiffness_ / (.5 * mass))
559  {
560  std::cerr << "Error in getCollisionTime(" << mass
561  << ") values for mass, stiffness and dissipation would lead to an overdamped system, (getCollisionTime("
562  << mass << "), with stiffness=" << loadingStiffness_ << " and dissipation=" << dissipation_ << ")"
563  << std::endl;
564  exit(-1);
565  }
566  return constants::pi / std::sqrt(tosqrt);
567 }
void setParhamiMcKeeping(Mdouble alpha, Mdouble beta, Mdouble atomicVolume, Mdouble surfaceEnergy, Mdouble thicknessDiffusion, Mdouble activationEnergy, Mdouble temperature)
Sets the sinterAdhesion_ and inverseSinterViscosity_ based on the Parhami-McKeeping parameters...
Mdouble unloadingStiffnessMax_
the maximum elastic constant (k_2^max) for plastic deformations
SINTERTYPE sinterType_
sinterType options determin how the rate of sintering, d(equilibriumOverlap)/dt is computed PARHAMI_M...
SinterNormalSpecies()
The default constructor.
void setSeparationDis(Mdouble separationDis)
void mix(SinterNormalSpecies *S, SinterNormalSpecies *T)
creates default values for mixed species
void setDissipation(Mdouble dissipation)
Sets the linear dissipation coefficient of the linear plastic-viscoelastic normal force...
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
return type specifically for fuctions returning k and disp at once
Definition: Helpers.h:45
LL< Log::INFO > INFO
Info log level.
Definition: Logger.cc:55
void setCohesionStiffness(Mdouble cohesionStiffness)
Sets the cohesive stiffness of the linear plastic-viscoelastic normal force.
void setTemperatureDependentSinterRate(std::function< double(double temperature)> temperatureDependentSinterRate)
Mdouble exp(Mdouble Exponent)
Definition: ExtendedMath.cc:84
Mdouble loadingStiffness_
(normal) spring constant (k_1)
Mdouble cohesionStiffness_
the adhesive spring constant (k^c) for plastic deformations
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
void read(std::istream &is)
Reads the species properties from an input stream.
void setPenetrationDepthMax(Mdouble penetrationDepthMax)
Sets the maximum penetration depth of the linear plastic-viscoelastic normal force.
void setPlasticParameters(Mdouble loadingStiffness, Mdouble unloadingStiffnessMax, Mdouble cohesionStiffness, Mdouble penetrationDepthMax)
Sets all parameters of the linear plastic-viscoelastic normal force at once.
Mdouble getSinterAdhesion() const
Accesses sinterAdhesion_.
Mdouble computeTimeStep(Mdouble mass)
Returns the optimal time step to resolve a collision of two particles of a given mass.
void setStiffnessAndRestitutionCoefficient(Mdouble k_, Mdouble eps, Mdouble mass)
Sets k, disp such that it matches a given tc and eps for a collision of two copies of P...
void setSinterAdhesion(Mdouble sinterAdhesion)
Sets sinterAdhesion_.
Mdouble penetrationDepthMax_
the depth (relative to the normalized radius) at which k_2^max is used (phi_f)
Mdouble getCollisionTime(Mdouble mass)
Calculates collision time for two copies of a particle of given disp, k, mass.
Mdouble log(Mdouble Power)
void setCollisionTimeAndRestitutionCoefficient(Mdouble tc, Mdouble eps, Mdouble mass)
Set k, disp such that is matches a given tc and eps for a collision of two different masses...
SINTERTYPE getSinterType() const
void setInverseSinterViscosity(Mdouble inverseSinterViscosity)
Sets inverseSinterViscosity_.
void write(std::ostream &os) const
Writes the species properties to an output stream.
void setSinterRate(Mdouble sinterRate)
Sets sinterRate_.
void setSinterForceAndTime(Mdouble adhesionForce, Mdouble sinterTime, Mdouble radius)
void setSurfTension(Mdouble complianceZero)
MERCURY_DEPRECATED void setLoadingStiffnessAndDissipation(helpers::KAndDisp new_)
Allows the spring and dissipation constants to be changed simultaneously.
Mdouble disp
Definition: Helpers.h:49
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble dissipation_
linear dissipation coefficient
std::function< double(double temperature)> temperatureDependentSinterRate_
static Mdouble average(Mdouble a, Mdouble b)
Returns the harmonic mean of two variables.
Definition: BaseSpecies.cc:110
Mdouble getComplianceZero() const
Mdouble getLoadingStiffness() const
Returns the loading stiffness of the linear plastic-viscoelastic normal force.
Mdouble getCohesionStiffness() const
Returns the cohesive stiffness of the linear plastic-viscoelastic normal force.
void setConstantC1(Mdouble constantC1_)
Mdouble getSeparationDis() const
~SinterNormalSpecies()
The default destructor.
void setLoadingStiffness(Mdouble loadingStiffness)
Sets the loading stiffness of the linear plastic-viscoelastic normal force.
std::string getBaseName() const
Used in Species::getName to obtain a unique name for each Species.
void setComplianceZero(Mdouble complianceZero)
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
void setUnloadingStiffnessMax(Mdouble unloadingStiffnessMax)
Sets the maximum unloading stiffness of the linear plastic-viscoelastic normal force.
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.
Mdouble getSurfTension() const
void setSinterType(SINTERTYPE sinterType)
Sets sinterRate_.
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
Mdouble getInverseSinterViscosity() const
Accesses inverseSinterViscosity_.
Mdouble getConstantC1() const
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
Mdouble getSinterRate() const
Accesses sinterRate_.