1 /* ---------------------------------------------------------------------- 2 This is the 3 4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗ 5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝ 6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗ 7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║ 8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║ 9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝® 10 11 DEM simulation engine, released by 12 DCS Computing Gmbh, Linz, Austria 13 http://www.dcs-computing.com, office@dcs-computing.com 14 15 LIGGGHTS® is part of CFDEM®project: 16 http://www.liggghts.com | http://www.cfdem.com 17 18 Core developer and main author: 19 Christoph Kloss, christoph.kloss@dcs-computing.com 20 21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public 22 License, version 2 or later. It is distributed in the hope that it will 23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have 25 received a copy of the GNU General Public License along with LIGGGHTS®. 26 If not, see http://www.gnu.org/licenses . See also top-level README 27 and LICENSE files. 28 29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH, 30 the producer of the LIGGGHTS® software and the CFDEM®coupling software 31 See http://www.cfdem.com/terms-trademark-policy for details. 32 33 ------------------------------------------------------------------------- 34 Contributing author and copyright for this file: 35 Christoph Kloss (DCS Computing GmbH, Linz) 36 Christoph Kloss (JKU Linz) 37 Richard Berger (JKU Linz) 38 39 Copyright 2012- DCS Computing GmbH, Linz 40 Copyright 2009-2012 JKU Linz 41 ------------------------------------------------------------------------- */ 42 43 #ifdef NORMAL_MODEL 44 45 NORMAL_MODEL(HOOKE,hooke,0) 46 47 #else 48 49 #ifndef NORMAL_MODEL_HOOKE_H_ 50 #define NORMAL_MODEL_HOOKE_H_ 51 52 #include "contact_models.h" 53 #include <cmath> 54 #include "atom.h" 55 #include "force.h" 56 #include "update.h" 57 #include "normal_model_base.h" 58 59 namespace LIGGGHTS { 60 namespace ContactModels 61 { 62 template<> 63 class NormalModel<HOOKE> : public NormalModelBase 64 { 65 public: 66 NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup,class ContactModelBase *c) : 67 NormalModelBase(lmp, hsetup, c), 68 Yeff(NULL), 69 Geff(NULL), 70 coeffRestMax(NULL), 71 coeffRestLog(NULL), 72 coeffMu(NULL), 73 coeffStc(NULL), 74 charVel(0.0), 75 viscous(false), 76 tangential_damping(false), 77 limitForce(false), 78 ktToKn(false), 79 displayedSettings(false), 80 heating(false), 81 heating_track(false), 82 elastic_potential_offset_(-1), 83 elasticpotflag_(false), 84 fix_dissipated_(NULL), 85 dissipatedflag_(false), 86 overlap_offset_(0.0), 87 disable_when_bonded_(false), 88 bond_history_offset_(-1), 89 dissipation_history_offset_(-1), 90 cmb(c) 91 { 92 93 } 94 95 inline void registerSettings(Settings & settings) 96 { 97 settings.registerOnOff("viscous", viscous); 98 settings.registerOnOff("tangential_damping", tangential_damping, true); 99 settings.registerOnOff("limitForce", limitForce); 100 settings.registerOnOff("ktToKnUser", ktToKn); 101 settings.registerOnOff("heating_normal_hooke",heating,false); 102 settings.registerOnOff("heating_tracking",heating_track,false); 103 settings.registerOnOff("computeElasticPotential", elasticpotflag_, false); 104 settings.registerOnOff("computeDissipatedEnergy", dissipatedflag_, false); 105 settings.registerOnOff("disableNormalWhenBonded", disable_when_bonded_, false); 106 } 107 108 inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) 109 { 110 if (elasticpotflag_) 111 { 112 elastic_potential_offset_ = cmb->get_history_offset("elastic_potential_normal"); 113 if (elastic_potential_offset_ == -1) 114 { 115 elastic_potential_offset_ = hsetup->add_history_value("elastic_potential_normal", "0"); 116 hsetup->add_history_value("elastic_force_normal_0", "1"); 117 hsetup->add_history_value("elastic_force_normal_1", "1"); 118 hsetup->add_history_value("elastic_force_normal_2", "1"); 119 hsetup->add_history_value("elastic_torque_normal_i_0", "0"); 120 hsetup->add_history_value("elastic_torque_normal_i_1", "0"); 121 hsetup->add_history_value("elastic_torque_normal_i_2", "0"); 122 hsetup->add_history_value("elastic_torque_normal_j_0", "0"); 123 hsetup->add_history_value("elastic_torque_normal_j_1", "0"); 124 hsetup->add_history_value("elastic_torque_normal_j_2", "0"); 125 if (cmb->is_wall()) 126 hsetup->add_history_value("elastic_potential_wall", "0"); 127 cmb->add_history_offset("elastic_potential_normal", elastic_potential_offset_); 128 } 129 } 130 if (dissipatedflag_) 131 { 132 if (cmb->is_wall()) 133 { 134 fix_dissipated_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("dissipated_energy_wall", "property/atom", "vector", 0, 0, "dissipated energy")); 135 dissipation_history_offset_ = cmb->get_history_offset("dissipation_force"); 136 if (!dissipation_history_offset_) 137 error->one(FLERR, "Internal error: Could not find dissipation history offset"); 138 } 139 else 140 fix_dissipated_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("dissipated_energy", "property/atom", "vector", 0, 0, "dissipated energy")); 141 if (!fix_dissipated_) 142 error->one(FLERR, "Surface model has not registered dissipated_energy fix"); 143 } 144 if (disable_when_bonded_) 145 { 146 bond_history_offset_ = cmb->get_history_offset("bond_contactflag"); 147 if (bond_history_offset_ < 0) 148 error->one(FLERR, "Could not find bond history offset"); 149 overlap_offset_ = hsetup->add_history_value("overlap_offset", "0"); 150 } 151 } 152 153 inline void connectToProperties(PropertyRegistry & registry) { 154 registry.registerProperty("Yeff", &MODEL_PARAMS::createYeff); 155 registry.registerProperty("Geff", &MODEL_PARAMS::createGeff); 156 registry.registerProperty("charVel", &MODEL_PARAMS::createCharacteristicVelocity); 157 158 registry.connect("Yeff", Yeff,"model hooke"); 159 registry.connect("Geff", Geff,"model hooke"); 160 registry.connect("charVel", charVel,"model hooke"); 161 162 if(viscous) { 163 registry.registerProperty("coeffMu", &MODEL_PARAMS::createCoeffMu); 164 registry.registerProperty("coeffStc", &MODEL_PARAMS::createCoeffStc); 165 registry.registerProperty("coeffRestMax", &MODEL_PARAMS::createCoeffRestMax); 166 167 registry.connect("coeffMu", coeffMu,"model hooke viscous"); 168 registry.connect("coeffStc", coeffStc,"model hooke viscous"); 169 registry.connect("coeffRestMax", coeffRestMax,"model hooke viscous"); 170 //registry.connect("log(coeffRestMax)+coeffStc", logRestMaxPlusStc); 171 } else { 172 registry.registerProperty("coeffRestLog", &MODEL_PARAMS::createCoeffRestLog); 173 174 registry.connect("coeffRestLog", coeffRestLog,"model hooke viscous"); 175 } 176 177 // error checks on coarsegraining 178 if(force->cg_active()) 179 error->cg(FLERR,"model hooke"); 180 181 // enlarge contact distance flag in case of elastic energy computation 182 // to ensure that surfaceClose is called after a contact 183 if (elasticpotflag_) 184 //set neighbor contact_distance_factor here 185 neighbor->register_contact_dist_factor(1.01); 186 } 187 188 // effective exponent for stress-strain relationship 189 190 inline double stressStrainExponent() 191 { 192 return 1.; 193 } 194 195 void dissipateElasticPotential(SurfacesCloseData &scdata) 196 { 197 if (elasticpotflag_) 198 { 199 double * const elastic_energy = &scdata.contact_history[elastic_potential_offset_]; 200 if (scdata.is_wall) 201 { 202 // we need to calculate half an integration step which was left over to ensure no energy loss, but only for the elastic energy. The dissipation part is handled in fix_wall_gran_base.h. 203 double delta[3]; 204 scdata.fix_mesh->triMesh()->get_global_vel(delta); 205 vectorScalarMult3D(delta, update->dt); 206 // -= because force is in opposite direction 207 // no *dt as delta is v*dt of the contact position 208 elastic_energy[0] -= (delta[0]*(elastic_energy[1]) + 209 delta[1]*(elastic_energy[2]) + 210 delta[2]*(elastic_energy[3]))*0.5 211 // from previous half step 212 + elastic_energy[10]; 213 elastic_energy[10] = 0.0; 214 } 215 elastic_energy[1] = 0.0; 216 elastic_energy[2] = 0.0; 217 elastic_energy[3] = 0.0; 218 elastic_energy[4] = 0.0; 219 elastic_energy[5] = 0.0; 220 elastic_energy[6] = 0.0; 221 elastic_energy[7] = 0.0; 222 elastic_energy[8] = 0.0; 223 elastic_energy[9] = 0.0; 224 } 225 } 226 227 inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) 228 { 229 if (sidata.contact_flags) 230 *sidata.contact_flags |= CONTACT_NORMAL_MODEL; 231 const bool update_history = sidata.computeflag && sidata.shearupdate; 232 const int itype = sidata.itype; 233 const int jtype = sidata.jtype; 234 const double radi = sidata.radi; 235 const double radj = sidata.radj; 236 double reff=sidata.is_wall ? radi : (radi*radj/(radi+radj)); 237 #ifdef SUPERQUADRIC_ACTIVE_FLAG 238 if(sidata.is_non_spherical && atom->superquadric_flag) 239 reff = sidata.reff; 240 #endif 241 const double meff=sidata.meff; 242 double coeffRestLogChosen; 243 244 const double sqrtval = sqrt(reff); 245 246 if(!displayedSettings) 247 { 248 displayedSettings = true; 249 /* 250 if(ktToKn) 251 if(0 == comm->me) fprintf(screen," NormalModel<HOOKE>: will use user-modified ktToKn of 2/7.\n"); 252 if(tangential_damping) 253 if(0 == comm->me) fprintf(screen," NormalModel<HOOKE>: will apply tangential damping.\n"); 254 if(viscous) 255 if(0 == comm->me) fprintf(screen," NormalModel<HOOKE>: will apply damping based on Stokes number.\n"); 256 if(limitForce) 257 if(0 == comm->me) fprintf(screen," NormalModel<HOOKE>: will limit normal force.\n"); 258 */ 259 } 260 if (viscous) { 261 // Stokes Number from MW Schmeeckle (2001) 262 const double stokes=sidata.meff*sidata.vn/(6.0*M_PI*coeffMu[itype][jtype]*reff*reff); 263 // Empirical from Legendre (2006) 264 coeffRestLogChosen=log(coeffRestMax[itype][jtype])+coeffStc[itype][jtype]/stokes; 265 } else { 266 coeffRestLogChosen=coeffRestLog[itype][jtype]; 267 } 268 269 double kn = 16./15.*sqrtval*(Yeff[itype][jtype])*pow(15.*meff*charVel*charVel/(16.*sqrtval*Yeff[itype][jtype]),0.2); 270 double kt = kn; 271 if(ktToKn) kt *= 0.285714286; //2//7 272 const double coeffRestLogChosenSq = coeffRestLogChosen*coeffRestLogChosen; 273 //const double gamman=sqrt(4.*meff*kn/(1.+(M_PI/coeffRestLogChosen)*(M_PI/coeffRestLogChosen))); 274 const double gamman=sqrt(4.*meff*kn*coeffRestLogChosenSq/(coeffRestLogChosenSq+M_PI*M_PI)); 275 const double gammat = tangential_damping ? gamman : 0.0; 276 277 // convert Kn and Kt from pressure units to force/distance^2 278 kn /= force->nktv2p; 279 kt /= force->nktv2p; 280 281 const double Fn_damping = -gamman*sidata.vn; 282 if (disable_when_bonded_ && update_history && sidata.deltan < sidata.contact_history[overlap_offset_]) 283 sidata.contact_history[overlap_offset_] = sidata.deltan; 284 const double deltan = disable_when_bonded_ ? fmax(sidata.deltan-sidata.contact_history[overlap_offset_], 0.0) : sidata.deltan; 285 const double Fn_contact = kn*deltan; 286 double Fn = Fn_damping + Fn_contact; 287 288 //limit force to avoid the artefact of negative repulsion force 289 if(limitForce && (Fn<0.0) ) 290 { 291 Fn = 0.0; 292 } 293 sidata.Fn = Fn; 294 295 sidata.kn = kn; 296 sidata.kt = kt; 297 sidata.gamman = gamman; 298 sidata.gammat = gammat; 299 300 #ifdef NONSPHERICAL_ACTIVE_FLAG 301 double torque_i[3] = {0., 0., 0.}; 302 double Fn_i[3] = { Fn * sidata.en[0], Fn * sidata.en[1], Fn * sidata.en[2]}; 303 if(sidata.is_non_spherical) { 304 double xci[3]; 305 vectorSubtract3D(sidata.contact_point, atom->x[sidata.i], xci); 306 vectorCross3D(xci, Fn_i, torque_i); 307 } 308 #endif 309 310 // apply normal force 311 if (!disable_when_bonded_ || sidata.contact_history[bond_history_offset_] < 0.5) 312 { 313 314 if(heating) 315 { 316 sidata.P_diss += fabs(Fn_damping*sidata.vn); //fprintf(screen," contrib %f\n",fabs(Fn_damping*sidata.vn)); 317 if(heating_track && sidata.is_wall) cmb->tally_pw(fabs(Fn_damping*sidata.vn),sidata.i,jtype,0); 318 if(heating_track && !sidata.is_wall) cmb->tally_pp(fabs(Fn_damping*sidata.vn),sidata.i,sidata.j,0); 319 } 320 321 // energy balance terms 322 if (update_history) 323 { 324 // compute increment in elastic potential 325 if (elasticpotflag_) 326 { 327 double * const elastic_energy = &sidata.contact_history[elastic_potential_offset_]; 328 // correct for wall influence 329 if (sidata.is_wall) 330 { 331 double delta[3]; 332 sidata.fix_mesh->triMesh()->get_global_vel(delta); 333 vectorScalarMult3D(delta, update->dt); 334 // -= because force is in opposite direction 335 // no *dt as delta is v*dt of the contact position 336 //printf("pela %e %e %e %e\n", update->get_cur_time()-update->dt, deb, -sidata.radj, deb-sidata.radj); 337 elastic_energy[0] -= (delta[0]*elastic_energy[1] + 338 delta[1]*elastic_energy[2] + 339 delta[2]*elastic_energy[3])*0.5 340 // from previous half step 341 + elastic_energy[10]; 342 elastic_energy[10] = -(delta[0]*Fn_contact*sidata.en[0] + 343 delta[1]*Fn_contact*sidata.en[1] + 344 delta[2]*Fn_contact*sidata.en[2])*0.5; 345 } 346 elastic_energy[1] = -Fn_contact*sidata.en[0]; 347 elastic_energy[2] = -Fn_contact*sidata.en[1]; 348 elastic_energy[3] = -Fn_contact*sidata.en[2]; 349 elastic_energy[4] = 0.0; 350 elastic_energy[5] = 0.0; 351 elastic_energy[6] = 0.0; 352 elastic_energy[7] = 0.0; 353 elastic_energy[8] = 0.0; 354 elastic_energy[9] = 0.0; 355 } 356 // compute increment in dissipated energy 357 if (dissipatedflag_) 358 { 359 double * const * const dissipated = fix_dissipated_->array_atom; 360 double * const dissipated_i = dissipated[sidata.i]; 361 double * const dissipated_j = dissipated[sidata.j]; 362 const double F_diss = -Fn_damping; 363 dissipated_i[1] += sidata.en[0]*F_diss; 364 dissipated_i[2] += sidata.en[1]*F_diss; 365 dissipated_i[3] += sidata.en[2]*F_diss; 366 if (sidata.j < atom->nlocal && !sidata.is_wall) 367 { 368 dissipated_j[1] -= sidata.en[0]*F_diss; 369 dissipated_j[2] -= sidata.en[1]*F_diss; 370 dissipated_j[3] -= sidata.en[2]*F_diss; 371 } 372 else if (sidata.is_wall) 373 { 374 double * const diss_force = &sidata.contact_history[dissipation_history_offset_]; 375 diss_force[0] -= sidata.en[0]*F_diss; 376 diss_force[1] -= sidata.en[1]*F_diss; 377 diss_force[2] -= sidata.en[2]*F_diss; 378 } 379 } 380 #ifdef NONSPHERICAL_ACTIVE_FLAG 381 if ((dissipatedflag_ || elasticpotflag_) && sidata.is_non_spherical) 382 error->one(FLERR, "Dissipation and elastic potential do not compute torque influence for nonspherical particles"); 383 #endif 384 } 385 386 // apply normal force 387 if(sidata.is_wall) { 388 const double Fn_ = Fn * sidata.area_ratio; 389 i_forces.delta_F[0] += Fn_ * sidata.en[0]; 390 i_forces.delta_F[1] += Fn_ * sidata.en[1]; 391 i_forces.delta_F[2] += Fn_ * sidata.en[2]; 392 #ifdef NONSPHERICAL_ACTIVE_FLAG 393 if(sidata.is_non_spherical) { 394 //for non-spherical particles normal force can produce torque! 395 i_forces.delta_torque[0] += torque_i[0]; 396 i_forces.delta_torque[1] += torque_i[1]; 397 i_forces.delta_torque[2] += torque_i[2]; 398 } 399 #endif 400 } else { 401 const double fx = sidata.Fn * sidata.en[0]; 402 const double fy = sidata.Fn * sidata.en[1]; 403 const double fz = sidata.Fn * sidata.en[2]; 404 405 i_forces.delta_F[0] += fx; 406 i_forces.delta_F[1] += fy; 407 i_forces.delta_F[2] += fz; 408 409 j_forces.delta_F[0] += -fx; 410 j_forces.delta_F[1] += -fy; 411 j_forces.delta_F[2] += -fz; 412 #ifdef NONSPHERICAL_ACTIVE_FLAG 413 if(sidata.is_non_spherical) { 414 //for non-spherical particles normal force can produce torque! 415 double xcj[3], torque_j[3]; 416 double Fn_j[3] = { -Fn_i[0], -Fn_i[1], -Fn_i[2]}; 417 vectorSubtract3D(sidata.contact_point, atom->x[sidata.j], xcj); 418 vectorCross3D(xcj, Fn_j, torque_j); 419 420 i_forces.delta_torque[0] += torque_i[0]; 421 i_forces.delta_torque[1] += torque_i[1]; 422 i_forces.delta_torque[2] += torque_i[2]; 423 424 j_forces.delta_torque[0] += torque_j[0]; 425 j_forces.delta_torque[1] += torque_j[1]; 426 j_forces.delta_torque[2] += torque_j[2]; 427 } 428 #endif 429 } 430 } 431 else if (update_history) 432 { 433 sidata.contact_history[overlap_offset_] = sidata.deltan; 434 dissipateElasticPotential(sidata); 435 } 436 } 437 438 void surfacesClose(SurfacesCloseData &scdata, ForceData&, ForceData&) 439 { 440 if (scdata.contact_flags) 441 *scdata.contact_flags |= CONTACT_NORMAL_MODEL; 442 dissipateElasticPotential(scdata); 443 } 444 445 void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 446 void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 447 448 protected: 449 double ** Yeff; 450 double ** Geff; 451 double ** coeffRestMax; 452 double ** coeffRestLog; 453 double ** coeffMu; 454 double ** coeffStc; 455 double charVel; 456 457 bool viscous; 458 bool tangential_damping; 459 bool limitForce; 460 bool ktToKn; 461 bool displayedSettings; 462 bool heating; 463 bool heating_track; 464 int elastic_potential_offset_; 465 bool elasticpotflag_; 466 FixPropertyAtom *fix_dissipated_; 467 bool dissipatedflag_; 468 int overlap_offset_; 469 bool disable_when_bonded_; 470 int bond_history_offset_; 471 int dissipation_history_offset_; 472 class ContactModelBase *cmb; 473 }; 474 } 475 } 476 #endif // NORMAL_MODEL_HOOKE_H_ 477 #endif 478