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 Rahul Mohanty (University of Edinburgh, P&G) 37 Tomaz M. Zorec (University of Ljubljana) 38 ------------------------------------------------------------------------- */ 39 #ifdef NORMAL_MODEL 40 NORMAL_MODEL(THORNTON_NING,thornton_ning,8) 41 #else 42 #ifndef NORMAL_THORNTON_NING_H_ 43 #define NORMAL_THORNTON_NING_H_ 44 #include "contact_models.h" 45 #include "normal_model_base.h" 46 #include <cmath> 47 #include "atom.h" 48 #include "force.h" 49 #include "update.h" 50 #include "global_properties.h" 51 #include "math_extra_liggghts.h" 52 namespace LIGGGHTS { 53 54 namespace ContactModels 55 { 56 template<> 57 class NormalModel<THORNTON_NING> : public NormalModelBase 58 { 59 public: 60 61 NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup,class ContactModelBase *c) : 62 NormalModelBase(lmp, hsetup, c), 63 Yeff(NULL), 64 Geff(NULL), 65 betaeff(NULL), 66 limitForce(false), 67 displayedSettings(false) 68 { 69 history_offset = hsetup->add_history_value("tn_virgin_flag", "1"); 70 hsetup->add_history_value("delta_old", "1"); 71 hsetup->add_history_value("delta_max", "1"); 72 hsetup->add_history_value("force_old", "1"); 73 hsetup->add_history_value("force_max", "1"); 74 hsetup->add_history_value("adhesion_flag", "1"); 75 hsetup->add_history_value("detaching_delta", "1"); 76 hsetup->add_history_value("detaching_flag", "1"); 77 hsetup->add_history_value("detaching_force", "1"); 78 hsetup->add_history_value("yielding_flag", "1"); 79 kc_offset = hsetup->add_history_value("kc", "1"); 80 fo_offset = hsetup->add_history_value("fo", "1"); 81 c->add_history_offset("kc_offset", kc_offset); 82 c->add_history_offset("fo_offset", fo_offset); 83 84 } 85 86 void registerSettings(Settings & settings) 87 { 88 settings.registerOnOff("tangential_damping", tangential_damping, true); 89 settings.registerOnOff("limitForce", limitForce); 90 } 91 92 inline void postSettings(IContactHistorySetup *, ContactModelBase *) {} 93 94 void connectToProperties(PropertyRegistry & registry) { 95 registry.registerProperty("Yeff", &MODEL_PARAMS::createYeff,"model thornton_ning"); 96 registry.registerProperty("Geff", &MODEL_PARAMS::createGeff,"model thornton_ning"); 97 registry.registerProperty("betaeff", &MODEL_PARAMS::createBetaEff,"model thornton_ning"); 98 registry.registerProperty("gamma_surf", &MODEL_PARAMS::createSurfaceEnergy, "model thornton_ning"); 99 registry.registerProperty("yield_ratio", &MODEL_PARAMS::createYieldRatio, "model thornton_ning"); 100 // registry.registerProperty("coeffRestLog", &MODEL_PARAMS::createCoeffRestLog); 101 102 registry.connect("Yeff", Yeff,"model thornton_ning"); 103 registry.connect("Geff", Geff,"model thornton_ning"); 104 registry.connect("betaeff", betaeff,"model thornton_ning"); 105 registry.connect("gamma_surf", gamma_surf,"model thornton_ning"); 106 registry.connect("yield_ratio", yield_ratio, "model thornton_ning"); 107 // registry.connect("coeffRestLog", coeffRestLog,"model thornton_ning_oblique"); 108 } 109 110 // effective exponent for stress-strain relationship 111 112 inline double stressStrainExponent() 113 { 114 return 1.5; 115 } 116 117 /* ------------------------ CALL FUNCTIONS --------------------------------*/ 118 119 inline double calculate_fl(double force_old, double fc, int adhesion_flag) 120 { 121 122 double fl; 123 124 if (adhesion_flag == 1) 125 { 126 fl = force_old + 2*(fc-sqrt(fc * (force_old + fc))); 127 } else { 128 fl = force_old + 2*(fc+sqrt(fc * (force_old + fc))); 129 } 130 131 return fl; 132 } 133 134 inline double calculate_elastic_force_differential(double a, double E, double fl, double fc) 135 { 136 137 double df; 138 if (fc != 0) 139 { 140 double up = 3*sqrt(fl) - 3*sqrt(fc); 141 double down = 3*sqrt(fl) - sqrt(fc); 142 df = 2*E*a*up/down; 143 } else { 144 df = 2*E*a; 145 } 146 147 return df; 148 } 149 150 inline double calculate_plastic_force_differential(double a_yield, double E, double fl, double fc, double yield_stress, double R) 151 { 152 153 double df; 154 if (fc != 0) 155 { 156 double up = 3 * M_PI * R * yield_stress *sqrt(fl) - 2 * a_yield * E * sqrt(fc); 157 double down = 3*sqrt(fl) - sqrt(fc); 158 df = up/down; 159 } else { 160 df = 2*E*a_yield; 161 } 162 163 return df; 164 } 165 166 inline double calculate_a(double fl, double E, double R, double delta, double gamma_s) 167 { 168 double a; 169 170 if (gamma_s == 0) 171 { 172 a = sqrt(R * delta); 173 } else { 174 a = pow((3. * R * fl)/(4. * E), (1./3.)); 175 } 176 177 return a; 178 } 179 180 inline double calculate_delta_f_ratio(double a, double ac) 181 { 182 double a_over_ac = a/ac; 183 double daf = pow(3, (1./3.))*pow(a_over_ac, 2)*(1-(4./3.)*pow(a_over_ac, -(3./2.))); 184 185 return daf; 186 } 187 188 inline void branch_2_force_calc(double force_old, double fc, int adhesion_flag, double E, double rp, double delta, double dn, double gamma_s, double *f_df) 189 { 190 double fl_r = calculate_fl(force_old, fc, adhesion_flag); 191 double a = calculate_a(fl_r, E, rp, delta, gamma_s); 192 193 /* --------------- For checking the model --------------------*/ 194 // double ac = calculate_a(fc, E, rp, delta, gamma_s); 195 // double daf = calculate_delta_f_ratio(a, ac); 196 /* --------------- End model checking ------------------------*/ 197 198 double df = calculate_elastic_force_differential(a, E, fl_r, fc) * dn; 199 200 f_df[0] = force_old + df; 201 f_df[1] = df; 202 203 } 204 205 /* ------------------------- END OF CALL FUNCTIONS ------------------------*/ 206 207 inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) 208 { 209 const int itype = sidata.itype; 210 const int jtype = sidata.jtype; 211 double ri = sidata.radi; 212 double rj = sidata.radj; 213 double reff=sidata.is_wall ? sidata.radi : (ri*rj/(ri+rj)); 214 #ifdef SUPERQUADRIC_ACTIVE_FLAG 215 if(sidata.is_non_spherical && atom->superquadric_flag) 216 reff = sidata.reff; 217 #endif 218 double meff=sidata.meff; 219 double Eeff = Yeff[itype][jtype]; 220 const double gamma_s = gamma_surf[itype][jtype]; 221 const double sqrtFiveOverSix = 0.91287092917527685576161630466800355658790782499663875; 222 223 // load history 224 if(sidata.contact_flags) *sidata.contact_flags |= CONTACT_NORMAL_MODEL; 225 double * const history = &sidata.contact_history[history_offset]; 226 double * const kc_history = &sidata.contact_history[kc_offset]; 227 double * const fo_history = &sidata.contact_history[fo_offset]; 228 229 int virgin_loading_flag = history[0]; 230 double delta_old = history[1]; 231 double delta_max = history[2]; 232 double force_old = history[3]; 233 double force_max = history[4]; 234 int adhesion_flag = history[5]; 235 double detaching_delta = history[6]; 236 int detaching_flag = history[7]; 237 double detaching_force = history[8]; 238 int yielding_flag = history[9]; 239 240 double fc = gamma_s * reff * M_PI * (3./2.); 241 fo_history[0] = fc; 242 const double delta = sidata.deltan; 243 const double a_yield = ri * yield_ratio[itype]; 244 const double yield_stress = ((2.*Eeff*a_yield)/(M_PI*reff)) - sqrt(2.*gamma_s*Eeff/(M_PI*a_yield)); 245 const double limit_yield_stress = pow(2.*Eeff*Eeff*gamma_s/(M_PI*M_PI*reff), (1./3.)); 246 247 int plastic_from_start = 0; 248 249 if (yield_stress < limit_yield_stress) plastic_from_start = 1; 250 251 if (yield_stress <= 0) error->all(FLERR,"Invalid yield stress, please check surface energy and yield ratio!"); 252 253 double f, fl, fl_max; //forces 254 double df; // differentials 255 double a; // areas, plastic radius 256 257 // Initialization at contact 258 259 double k_t, gammat; 260 const double sqrtval = sqrt(reff*sidata.deltan); 261 // double coeffRestLogChosen = coeffRestLog[itype][jtype]; 262 263 k_t = 8.*Geff[itype][jtype]*sqrtval; 264 gammat = 2.*sqrtFiveOverSix*betaeff[itype][jtype]*sqrt(k_t*meff); 265 266 if ((virgin_loading_flag != 1)) 267 { 268 force_old = -(8./9.)*fc; 269 delta_old = 0; 270 virgin_loading_flag = 1; // update virgin loading flag 271 force_max = force_old; 272 delta_max = 0; 273 adhesion_flag = 0; 274 detaching_delta = 0; 275 detaching_flag = 0; 276 detaching_force = -(5./9.)*fc; 277 yielding_flag = 0; 278 } 279 280 if (detaching_flag == 1) 281 { 282 if (delta >= detaching_delta) 283 { 284 detaching_flag = 0; 285 delta_old = detaching_delta; 286 force_old = detaching_force; 287 goto UNLOADING_RELOADING; 288 289 } else { 290 f = 0; 291 } 292 } else { 293 UNLOADING_RELOADING: 294 const double dn = delta - delta_old; 295 296 if (delta >= delta_max) // loading or unloading/reloading 297 { 298 delta_max = delta; 299 300 fl = calculate_fl(force_old, fc, 0); 301 a = calculate_a(fl, Eeff, reff, delta, gamma_s); 302 if (yielding_flag == 0) if (a >= a_yield) yielding_flag = 1; 303 304 if (plastic_from_start == 1 || yielding_flag == 1) // if particle already yielded go to plastic, else load elastically 305 { 306 307 df = calculate_plastic_force_differential(a_yield, Eeff, fl, fc, yield_stress, reff) * dn; 308 f = force_old + df; 309 310 } else { 311 312 double f_df[2]; 313 branch_2_force_calc(force_old, fc, adhesion_flag, Eeff, reff, delta, dn, gamma_s, f_df); 314 f = f_df[0]; 315 df = f_df[1]; 316 317 } 318 } else { 319 if (yielding_flag == 0) // if particle already yielded, the Reff must be modified, if not we use the original Reff 320 { 321 322 double f_df[2]; 323 branch_2_force_calc(force_old, fc, adhesion_flag, Eeff, reff, delta, dn, gamma_s, f_df); 324 f = f_df[0]; 325 df = f_df[1]; 326 327 } else { 328 329 // Calculate modified R and Fc based on that 330 331 fl_max = calculate_fl(force_max, fc, 0); 332 reff = reff * fl_max / (force_max + sqrt(4*fc*fl_max)); 333 fc = (3./2.) * M_PI * gamma_s * reff; 334 fo_history[0] = -fc; 335 double f_df[2]; 336 branch_2_force_calc(force_old, fc, adhesion_flag, Eeff, reff, delta, dn, gamma_s, f_df); 337 f = f_df[0]; 338 df = f_df[1]; 339 340 } 341 } 342 343 if (gamma_s == 0 && f < 0) // the case of no cohesion 344 { 345 f = 0; 346 detaching_flag = 1; 347 detaching_delta = delta; 348 detaching_force = 0; 349 } 350 351 if (f < -fc) //to prevent stepping into nothingness 352 { 353 double f_df[2]; 354 adhesion_flag = 1 - adhesion_flag; 355 branch_2_force_calc(force_old, fc, adhesion_flag, Eeff, reff, delta, dn, gamma_s, f_df); 356 f = f_df[0]; 357 df = f_df[1]; 358 } 359 360 if ((adhesion_flag == 1) && (f + 3*df >= -(5./9.)*fc)) // to make sure particles detach and later reload from same point 361 { 362 f = 0; 363 364 detaching_flag = 1; 365 detaching_delta = delta_old; 366 detaching_force = force_old; 367 } 368 } 369 // put force where it belongs 370 sidata.Fn = f; 371 sidata.kt = k_t; 372 kc_history[0] = 0.0; 373 sidata.gammat = gammat; 374 //sidata.detaching_flag = detaching_flag; 375 376 force_old = f; 377 force_max = (f > force_max) ? f : force_max; 378 379 // save history values 380 history[0] = virgin_loading_flag; 381 history[1] = delta; 382 history[2] = delta_max; 383 history[3] = f; 384 history[4] = force_max; 385 history[5] = adhesion_flag; 386 history[6] = detaching_delta; 387 history[7] = detaching_flag; 388 history[8] = detaching_force; 389 history[9] = yielding_flag; 390 391 // apply normal force 392 393 if(sidata.is_wall) { 394 const double Fn_ = f * sidata.area_ratio; 395 i_forces.delta_F[0] = Fn_ * sidata.en[0]; 396 i_forces.delta_F[1] = Fn_ * sidata.en[1]; 397 i_forces.delta_F[2] = Fn_ * sidata.en[2]; 398 } else { 399 i_forces.delta_F[0] = sidata.Fn * sidata.en[0]; 400 i_forces.delta_F[1] = sidata.Fn * sidata.en[1]; 401 i_forces.delta_F[2] = sidata.Fn * sidata.en[2]; 402 403 j_forces.delta_F[0] = -i_forces.delta_F[0]; 404 j_forces.delta_F[1] = -i_forces.delta_F[1]; 405 j_forces.delta_F[2] = -i_forces.delta_F[2]; 406 } 407 } 408 409 void surfacesClose(SurfacesCloseData & scdata, ForceData&, ForceData&){ 410 if(scdata.contact_flags) *scdata.contact_flags &= ~CONTACT_NORMAL_MODEL; 411 double * const history = &scdata.contact_history[history_offset]; 412 history[0] = 0.0; 413 history[1] = 0.0; 414 history[2] = 0.0; 415 history[3] = 0.0; 416 history[4] = 0.0; 417 history[5] = 0.0; 418 history[6] = 0.0; 419 history[7] = 0.0; 420 history[8] = 0.0; 421 history[9] = 0.0; 422 } 423 void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 424 void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 425 426 protected: 427 double ** Yeff; 428 double ** Geff; 429 double ** betaeff; 430 double ** gamma_surf; 431 double * yield_ratio; 432 double ** coeffRestLog; 433 int history_offset; 434 int kc_offset; 435 int fo_offset; 436 437 bool tangential_damping; 438 bool limitForce; 439 bool displayedSettings; 440 class ContactModelBase *cmb; 441 }; 442 443 } 444 445 } 446 #endif 447 #endif 448