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 36 Christoph Kloss (DCS Computing GmbH, Linz) 37 Christoph Kloss (JKU Linz) 38 Richard Berger (JKU Linz) 39 40 Copyright 2012- DCS Computing GmbH, Linz 41 Copyright 2009-2012 JKU Linz 42 ------------------------------------------------------------------------- */ 43 44 #ifdef NORMAL_MODEL 45 NORMAL_MODEL(HOOKE_STIFFNESS,hooke/stiffness,1) 46 #else 47 #ifndef NORMAL_MODEL_HOOKE_STIFFNESS_H_ 48 #define NORMAL_MODEL_HOOKE_STIFFNESS_H_ 49 #include "contact_models.h" 50 #include "normal_model_base.h" 51 52 namespace LIGGGHTS { 53 namespace ContactModels 54 { 55 template<> 56 class NormalModel<HOOKE_STIFFNESS> : public NormalModelBase 57 { 58 public: 59 NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup, class ContactModelBase * c) : 60 NormalModelBase(lmp, hsetup, c), 61 k_n(NULL), 62 k_t(NULL), 63 gamma_n(NULL), 64 gamma_t(NULL), 65 tangential_damping(false), 66 absolute_damping(false), 67 limitForce(false), 68 displayedSettings(false), 69 elastic_potential_offset_(0), 70 elasticpotflag_(false), 71 fix_dissipated_(NULL), 72 dissipatedflag_(false), 73 dissipation_history_offset_(0) 74 { 75 76 } 77 78 void registerSettings(Settings & settings) 79 { 80 settings.registerOnOff("tangential_damping", tangential_damping, true); 81 settings.registerOnOff("absolute_damping", absolute_damping); 82 settings.registerOnOff("limitForce", limitForce); 83 settings.registerOnOff("computeElasticPotential", elasticpotflag_, false); 84 settings.registerOnOff("computeDissipatedEnergy", dissipatedflag_, false); 85 } 86 87 inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) 88 { 89 if (elasticpotflag_) 90 { 91 elastic_potential_offset_ = cmb->get_history_offset("elastic_potential_normal"); 92 if (elastic_potential_offset_ == -1) 93 { 94 elastic_potential_offset_ = hsetup->add_history_value("elastic_potential_normal", "0"); 95 hsetup->add_history_value("elastic_force_normal_0", "1"); 96 hsetup->add_history_value("elastic_force_normal_1", "1"); 97 hsetup->add_history_value("elastic_force_normal_2", "1"); 98 hsetup->add_history_value("elastic_torque_normal_i_0", "0"); 99 hsetup->add_history_value("elastic_torque_normal_i_1", "0"); 100 hsetup->add_history_value("elastic_torque_normal_i_2", "0"); 101 hsetup->add_history_value("elastic_torque_normal_j_0", "0"); 102 hsetup->add_history_value("elastic_torque_normal_j_1", "0"); 103 hsetup->add_history_value("elastic_torque_normal_j_2", "0"); 104 if (cmb->is_wall()) 105 hsetup->add_history_value("elastic_potential_wall", "0"); 106 cmb->add_history_offset("elastic_potential_normal", elastic_potential_offset_); 107 } 108 } 109 if (dissipatedflag_) 110 { 111 if (cmb->is_wall()) 112 { 113 fix_dissipated_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("dissipated_energy_wall", "property/atom", "vector", 0, 0, "dissipated energy")); 114 dissipation_history_offset_ = cmb->get_history_offset("dissipation_force"); 115 if (!dissipation_history_offset_) 116 error->one(FLERR, "Internal error: Could not find dissipation history offset"); 117 } 118 else 119 fix_dissipated_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("dissipated_energy", "property/atom", "vector", 0, 0, "dissipated energy")); 120 if (!fix_dissipated_) 121 error->one(FLERR, "Surface model has not registered dissipated_energy fix"); 122 } 123 } 124 125 void connectToProperties(PropertyRegistry & registry) { 126 registry.registerProperty("k_n", &MODEL_PARAMS::createKn); 127 registry.registerProperty("k_t", &MODEL_PARAMS::createKt); 128 129 registry.connect("k_n", k_n,"model hooke/stiffness"); 130 registry.connect("k_t", k_t,"model hooke/stiffness"); 131 132 if(absolute_damping) { 133 registry.registerProperty("gamman_abs", &MODEL_PARAMS::createGammanAbs); 134 registry.registerProperty("gammat_abs", &MODEL_PARAMS::createGammatAbs); 135 registry.connect("gamman_abs", gamma_n,"model hooke/stiffness"); 136 registry.connect("gammat_abs", gamma_t,"model hooke/stiffness"); 137 } else { 138 registry.registerProperty("gamman", &MODEL_PARAMS::createGamman); 139 registry.registerProperty("gammat", &MODEL_PARAMS::createGammat); 140 registry.connect("gamman", gamma_n,"model hooke/stiffness"); 141 registry.connect("gammat", gamma_t,"model hooke/stiffness"); 142 } 143 144 // error checks on coarsegraining 145 if(force->cg_active()) 146 error->cg(FLERR,"model hooke/stiffness"); 147 148 // enlarge contact distance flag in case of elastic energy computation 149 // to ensure that surfaceClose is called after a contact 150 if (elasticpotflag_) 151 //set neighbor contact_distance_factor here 152 neighbor->register_contact_dist_factor(1.01); 153 } 154 155 // effective exponent for stress-strain relationship 156 157 inline double stressStrainExponent() 158 { 159 return 1.; 160 } 161 162 void dissipateElasticPotential(SurfacesCloseData &scdata) 163 { 164 if (elasticpotflag_) 165 { 166 double * const elastic_energy = &scdata.contact_history[elastic_potential_offset_]; 167 if (scdata.is_wall) 168 { 169 // 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. 170 double delta[3]; 171 scdata.fix_mesh->triMesh()->get_global_vel(delta); 172 vectorScalarMult3D(delta, update->dt); 173 // -= because force is in opposite direction 174 // no *dt as delta is v*dt of the contact position 175 elastic_energy[0] -= (delta[0]*(elastic_energy[1]) + 176 delta[1]*(elastic_energy[2]) + 177 delta[2]*(elastic_energy[3]))*0.5 178 // from previous half step 179 + elastic_energy[10]; 180 elastic_energy[10] = 0.0; 181 } 182 elastic_energy[1] = 0.0; 183 elastic_energy[2] = 0.0; 184 elastic_energy[3] = 0.0; 185 elastic_energy[4] = 0.0; 186 elastic_energy[5] = 0.0; 187 elastic_energy[6] = 0.0; 188 elastic_energy[7] = 0.0; 189 elastic_energy[8] = 0.0; 190 elastic_energy[9] = 0.0; 191 } 192 } 193 194 inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) 195 { 196 if (sidata.contact_flags) 197 *sidata.contact_flags |= CONTACT_NORMAL_MODEL; 198 const bool update_history = sidata.computeflag && sidata.shearupdate; 199 const int itype = sidata.itype; 200 const int jtype = sidata.jtype; 201 double meff=sidata.meff; 202 203 double kn = k_n[itype][jtype]; 204 double kt = k_t[itype][jtype]; 205 double gamman, gammat; 206 207 if(!displayedSettings) 208 { 209 displayedSettings = true; 210 211 /* 212 if(limitForce) 213 if(0 == comm->me) fprintf(screen," NormalModel<HOOKE_STIFFNESS>: will limit normal force.\n"); 214 */ 215 } 216 if(absolute_damping) 217 { 218 gamman = gamma_n[itype][jtype]; 219 gammat = gamma_t[itype][jtype]; 220 } 221 else 222 { 223 gamman = meff*gamma_n[itype][jtype]; 224 gammat = meff*gamma_t[itype][jtype]; 225 } 226 227 if (!tangential_damping) gammat = 0.0; 228 229 // convert Kn and Kt from pressure units to force/distance^2 230 kn /= force->nktv2p; 231 kt /= force->nktv2p; 232 233 const double Fn_damping = -gamman*sidata.vn; 234 const double Fn_contact = kn*sidata.deltan; 235 double Fn = Fn_damping + Fn_contact; 236 237 //limit force to avoid the artefact of negative repulsion force 238 if(limitForce && (Fn<0.0) ) 239 { 240 Fn = 0.0; 241 } 242 243 sidata.Fn = Fn; 244 245 sidata.kn = kn; 246 sidata.kt = kt; 247 sidata.gamman = gamman; 248 sidata.gammat = gammat; 249 250 #ifdef NONSPHERICAL_ACTIVE_FLAG 251 double Fn_i[3] = { Fn * sidata.en[0], Fn * sidata.en[1], Fn * sidata.en[2]}; 252 double torque_i[3] = {0.0, 0.0, 0.0}; //initialized here with zeros to avoid compiler warnings 253 if(sidata.is_non_spherical) { 254 double xci[3]; 255 vectorSubtract3D(sidata.contact_point, atom->x[sidata.i], xci); 256 vectorCross3D(xci, Fn_i, torque_i); 257 } 258 #endif 259 260 // energy balance terms 261 if (update_history) 262 { 263 // compute increment in elastic potential 264 if (elasticpotflag_) 265 { 266 double * const elastic_energy = &sidata.contact_history[elastic_potential_offset_]; 267 // correct for wall influence 268 if (sidata.is_wall) 269 { 270 double delta[3]; 271 sidata.fix_mesh->triMesh()->get_global_vel(delta); 272 vectorScalarMult3D(delta, update->dt); 273 // -= because force is in opposite direction 274 // no *dt as delta is v*dt of the contact position 275 //printf("pela %e %e %e %e\n", update->get_cur_time()-update->dt, deb, -sidata.radj, deb-sidata.radj); 276 elastic_energy[0] -= (delta[0]*elastic_energy[1] + 277 delta[1]*elastic_energy[2] + 278 delta[2]*elastic_energy[3])*0.5 279 // from previous half step 280 + elastic_energy[10]; 281 elastic_energy[10] = -(delta[0]*Fn_contact*sidata.en[0] + 282 delta[1]*Fn_contact*sidata.en[1] + 283 delta[2]*Fn_contact*sidata.en[2])*0.5; 284 } 285 elastic_energy[1] = -Fn_contact*sidata.en[0]; 286 elastic_energy[2] = -Fn_contact*sidata.en[1]; 287 elastic_energy[3] = -Fn_contact*sidata.en[2]; 288 elastic_energy[4] = 0.0; 289 elastic_energy[5] = 0.0; 290 elastic_energy[6] = 0.0; 291 elastic_energy[7] = 0.0; 292 elastic_energy[8] = 0.0; 293 elastic_energy[9] = 0.0; 294 } 295 // compute increment in dissipated energy 296 if (dissipatedflag_) 297 { 298 double * const * const dissipated = fix_dissipated_->array_atom; 299 double * const dissipated_i = dissipated[sidata.i]; 300 double * const dissipated_j = dissipated[sidata.j]; 301 const double F_diss = -Fn_damping; 302 dissipated_i[1] += sidata.en[0]*F_diss; 303 dissipated_i[2] += sidata.en[1]*F_diss; 304 dissipated_i[3] += sidata.en[2]*F_diss; 305 if (sidata.j < atom->nlocal && !sidata.is_wall) 306 { 307 dissipated_j[1] -= sidata.en[0]*F_diss; 308 dissipated_j[2] -= sidata.en[1]*F_diss; 309 dissipated_j[3] -= sidata.en[2]*F_diss; 310 } 311 else if (sidata.is_wall) 312 { 313 double * const diss_force = &sidata.contact_history[dissipation_history_offset_]; 314 diss_force[0] -= sidata.en[0]*F_diss; 315 diss_force[1] -= sidata.en[1]*F_diss; 316 diss_force[2] -= sidata.en[2]*F_diss; 317 } 318 } 319 #ifdef NONSPHERICAL_ACTIVE_FLAG 320 if ((dissipatedflag_ || elasticpotflag_) && sidata.is_non_spherical) 321 error->one(FLERR, "Dissipation and elastic potential do not compute torque influence for nonspherical particles"); 322 #endif 323 } 324 325 // apply normal force 326 if(sidata.is_wall) { 327 const double Fn_ = Fn * sidata.area_ratio; 328 i_forces.delta_F[0] += Fn_ * sidata.en[0]; 329 i_forces.delta_F[1] += Fn_ * sidata.en[1]; 330 i_forces.delta_F[2] += Fn_ * sidata.en[2]; 331 #ifdef NONSPHERICAL_ACTIVE_FLAG 332 if(sidata.is_non_spherical) { 333 //for non-spherical particles normal force can produce torque! 334 i_forces.delta_torque[0] += torque_i[0]; 335 i_forces.delta_torque[1] += torque_i[1]; 336 i_forces.delta_torque[2] += torque_i[2]; 337 } 338 #endif 339 } else { 340 i_forces.delta_F[0] += sidata.Fn * sidata.en[0]; 341 i_forces.delta_F[1] += sidata.Fn * sidata.en[1]; 342 i_forces.delta_F[2] += sidata.Fn * sidata.en[2]; 343 344 j_forces.delta_F[0] += -i_forces.delta_F[0]; 345 j_forces.delta_F[1] += -i_forces.delta_F[1]; 346 j_forces.delta_F[2] += -i_forces.delta_F[2]; 347 #ifdef NONSPHERICAL_ACTIVE_FLAG 348 if(sidata.is_non_spherical) { 349 //for non-spherical particles normal force can produce torque! 350 double xcj[3], torque_j[3]; 351 double Fn_j[3] = { -Fn_i[0], -Fn_i[1], -Fn_i[2]}; 352 vectorSubtract3D(sidata.contact_point, atom->x[sidata.j], xcj); 353 vectorCross3D(xcj, Fn_j, torque_j); 354 355 i_forces.delta_torque[0] += torque_i[0]; 356 i_forces.delta_torque[1] += torque_i[1]; 357 i_forces.delta_torque[2] += torque_i[2]; 358 359 j_forces.delta_torque[0] += torque_j[0]; 360 j_forces.delta_torque[1] += torque_j[1]; 361 j_forces.delta_torque[2] += torque_j[2]; 362 } 363 #endif 364 } 365 } 366 367 void surfacesClose(SurfacesCloseData &scdata, ForceData&, ForceData&) 368 { 369 if (scdata.contact_flags) 370 *scdata.contact_flags |= CONTACT_NORMAL_MODEL; 371 dissipateElasticPotential(scdata); 372 } 373 374 void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 375 void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 376 377 protected: 378 double ** k_n; 379 double ** k_t; 380 double ** gamma_n; 381 double ** gamma_t; 382 383 bool tangential_damping; 384 bool absolute_damping; 385 bool limitForce; 386 bool displayedSettings; 387 int elastic_potential_offset_; 388 bool elasticpotflag_; 389 FixPropertyAtom *fix_dissipated_; 390 bool dissipatedflag_; 391 int dissipation_history_offset_; 392 }; 393 } 394 } 395 #endif // NORMAL_MODEL_HOOKE_STIFFNESS_H_ 396 #endif 397