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 Contributing authors: 37 Rahul Mohanty (University of Edinburgh, P&G) 38 Tomaz M. Zorec (University of Ljubljana) 39 ------------------------------------------------------------------------- */ 40 #ifdef NORMAL_MODEL 41 NORMAL_MODEL(EDINBURGH,edinburgh,10) 42 #else 43 #ifndef NORMAL_MODEL_EDINBURGH_H_ 44 #define NORMAL_MODEL_EDINBURGH_H_ 45 #include "contact_models.h" 46 #include "normal_model_base.h" 47 #include <cmath> 48 #include "atom.h" 49 #include "force.h" 50 #include "update.h" 51 #include "global_properties.h" 52 53 namespace LIGGGHTS { 54 namespace ContactModels 55 { 56 template<> 57 class NormalModel<EDINBURGH> : public NormalModelBase 58 { 59 public: 60 NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup,class ContactModelBase *c) : 61 NormalModelBase(lmp, hsetup, c), 62 Yeff(NULL), 63 Geff(NULL), 64 CoeffRestLog(NULL), 65 betaeff(NULL), 66 kn2kc(NULL), 67 kn2k1(NULL), 68 cex(0.0), 69 dex(0.0), 70 f_adh(NULL), 71 gamma_surf(NULL), 72 history_offset(-1), 73 kc_offset(-1), 74 fo_offset(-1), 75 tangential_damping(false), 76 limitForce(false), 77 fixKc(false) 78 { 79 history_offset = hsetup->add_history_value("deltaMax", "1"); 80 hsetup->add_history_value("old_delta", "1"); 81 kc_offset = hsetup->add_history_value("kc", "1"); 82 fo_offset = hsetup->add_history_value("fo", "1"); 83 c->add_history_offset("kc_offset", kc_offset); 84 c->add_history_offset("fo_offset", fo_offset); 85 } 86 87 void registerSettings(Settings & settings) 88 { 89 settings.registerOnOff("tangential_damping", tangential_damping, true); 90 settings.registerOnOff("limitForce", limitForce,true); 91 settings.registerOnOff("fixKc", fixKc); 92 //TODO error->one(FLERR,"TODO here also check if right surface model used"); 93 } 94 95 inline void postSettings(IContactHistorySetup *, ContactModelBase *) {} 96 97 void connectToProperties(PropertyRegistry & registry) { 98 99 registry.registerProperty("Yeff", &MODEL_PARAMS::createYeff,"model edinburgh"); // only used for non-linear model 100 registry.registerProperty("Geff", &MODEL_PARAMS::createGeff,"model edinburgh"); // only used for non-linear model 101 registry.registerProperty("CoeffRestLog", &MODEL_PARAMS::createCoeffRestLog, "model edinburgh"); 102 registry.registerProperty("betaeff", &MODEL_PARAMS::createBetaEff,"model edinburgh"); 103 registry.registerProperty("kn2kc", &MODEL_PARAMS::createCoeffAdhesionStiffness,"model edinburgh"); 104 registry.registerProperty("kn2k1", &MODEL_PARAMS::createUnloadingStiffness, "model edinburgh"); 105 registry.registerProperty("cex", &MODEL_PARAMS::createAdhesionExponent, "model edinburgh"); 106 registry.registerProperty("dex", &MODEL_PARAMS::createOverlapExponent, "model edinburgh"); 107 registry.registerProperty("f_adh", &MODEL_PARAMS::createPullOffForce, "model edinburgh"); 108 registry.registerProperty("gamma_surf", &MODEL_PARAMS::createSurfaceEnergy, "model edinburgh"); 109 110 registry.connect("Yeff", Yeff,"model edinburgh"); 111 registry.connect("Geff", Geff,"model edinburgh"); 112 registry.connect("CoeffRestLog", CoeffRestLog, "model edinburgh"); 113 registry.connect("betaeff", betaeff,"model edinburgh"); 114 registry.connect("kn2kc", kn2kc,"model edinburgh"); 115 registry.connect("kn2k1", kn2k1,"model edinburgh"); 116 registry.connect("cex", cex,"model edinburgh"); 117 registry.connect("dex", dex,"model edinburgh"); 118 registry.connect("f_adh", f_adh,"model edinburgh"); 119 registry.connect("gamma_surf", gamma_surf,"model edinburgh"); 120 121 } 122 123 // effective exponent for stress-strain relationship 124 125 inline double stressStrainExponent() 126 { 127 return 1.5; 128 } 129 130 inline double calculate_deltan_p_max (double deltan_p, double * const history, int count_flag, const double f_0, double fTmp, const double k2, double dex, double dex_i, double deltan, double k_adh) 131 { 132 //calculating the maximum overlap 133 double deltan_p_max; 134 135 if (count_flag == 0 ) { 136 if (deltan_p > history[0]) { 137 deltan_p_max = deltan_p; 138 history[0] = deltan_p_max; 139 } else { 140 deltan_p_max = history[0]; 141 } 142 } else { 143 deltan_p_max = pow((pow(history[1], dex) + ((k_adh/k2)*pow(history[1], cex))), dex_i); 144 history[0] = deltan_p_max; 145 } 146 return history[0]; 147 } 148 149 inline double whichd(double fTmp, double k1, double deltan_e, double d1, double d2) 150 { 151 //calculation distance between the particles based on branch 152 if (fTmp >= k1 * deltan_e){ 153 return d1; 154 } else{ 155 return d2; 156 } 157 } 158 159 inline double calculate_k_adh(double d, double risq, double rjsq, const double g_surf, double f_min_lim, double deltan_pe_max, double deltan_p_max, const double k2, const double dex_i, const double cex) 160 { 161 double new_k_c, delta_min, a, f_min, dsq; 162 dsq = d*d; 163 a = (1./(2.*d)) * sqrt(4* dsq * risq - ((dsq - rjsq + risq)*(dsq - rjsq + risq))); 164 f_min = ( 1.5 * M_PI * g_surf * a); 165 166 if (f_min > f_min_lim) 167 { 168 f_min = f_min_lim; 169 delta_min = 0.5 * deltan_p_max; 170 } else { 171 delta_min = pow((-f_min + k2 * deltan_pe_max)/k2, dex_i); 172 } 173 new_k_c = f_min/pow(delta_min, cex); 174 return new_k_c; 175 } 176 177 inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) 178 { 179 const int itype = sidata.itype; 180 const int jtype = sidata.jtype; 181 const double deltan = sidata.deltan; 182 double ri = sidata.radi; 183 double rj = sidata.radj; 184 double reff=sidata.is_wall ? sidata.radi : ((ri*rj)/(ri+rj)); 185 #ifdef SUPERQUADRIC_ACTIVE_FLAG 186 if(sidata.is_non_spherical && atom->superquadric_flag) 187 reff = sidata.reff; 188 #endif 189 double meff=sidata.meff; 190 const double f_0 = f_adh[itype][jtype]; 191 double sqrtval = sqrt(reff*sidata.deltan); 192 double Sn=2.*Yeff[itype][jtype]*sqrtval; 193 double St=8.*Geff[itype][jtype]*sqrtval; 194 double kn, kt; 195 196 if(dex==1){ 197 kn = 2.0*Yeff[itype][jtype]*reff; 198 kt = Sn; 199 } else { 200 kn = (4./3.)*Yeff[itype][jtype]*sqrt(reff); 201 kt = St; 202 } 203 204 // convert Kn and Kt from pressure units to force/distance^2 205 kn /= force->nktv2p; 206 kt /= force->nktv2p; 207 208 const double dex_i = 1./dex; 209 const double k1 = kn; 210 const double k2 = kn*kn2k1[itype][jtype]; 211 212 double gamman, gammat; 213 gamman = sqrt(4.*meff*kn/(1.+(M_PI/CoeffRestLog[itype][jtype])*(M_PI/CoeffRestLog[itype][jtype]))); 214 gammat = sqrt(4.*meff*kn/(1.+(M_PI/CoeffRestLog[itype][jtype])*(M_PI/CoeffRestLog[itype][jtype]))); 215 216 if(sidata.contact_flags) *sidata.contact_flags |= CONTACT_NORMAL_MODEL; 217 double * const history = &sidata.contact_history[history_offset]; 218 double * const kc_history = &sidata.contact_history[kc_offset]; 219 double * const fo_history = &sidata.contact_history[fo_offset]; 220 221 if (!tangential_damping) gammat = 0.0; 222 223 double Fn_contact, deltan_e, deltan_ce, deltaMax_e; 224 double k_adh = 0.0; 225 226 if (fixKc == false) { 227 const double lambda = pow((1. - k1/k2), dex_i); 228 int count_flag = 0; 229 230 // get the history value -- maximal overlap 231 232 double deltan_p = lambda*deltan; 233 double deltan_p_max, deltan_pe_max; 234 235 deltan_e = pow(deltan, dex); 236 deltan_ce = pow(deltan, cex); 237 238 double fTmp = k2*(deltan_e - history[0]); 239 double d, d1, d2, r_sum; 240 241 deltan_p_max = calculate_deltan_p_max (deltan_p, &sidata.contact_history[history_offset], count_flag, f_0, fTmp, k2, dex, dex_i, deltan, k_adh); 242 243 // Normal force calculation for Edinburgh model 244 245 const double g_surf = gamma_surf[itype][jtype]; 246 247 r_sum = sidata.radsum; 248 d1 = sidata.is_wall ? ri : r_sum - deltan; 249 d2 = sidata.is_wall ? ri : r_sum - deltan_p_max; 250 double risq = ri*ri; 251 double rjsq = rj*rj; 252 253 temp_calc: 254 255 deltan_pe_max = pow(deltan_p_max, dex); 256 const double f_min_lim = (k2 * deltan_pe_max) * 0.5; 257 fTmp=k2*(deltan_e - deltan_pe_max); 258 259 if (fTmp >= k1 * deltan_e){ // loading 260 Fn_contact = k1 * deltan_e; 261 }else{ 262 d = whichd(fTmp, k1, deltan_e, d1, d2); 263 k_adh = calculate_k_adh(d, risq, rjsq, g_surf, f_min_lim, deltan_pe_max, deltan_p_max, k2, dex_i, cex); 264 kc_history[0] = k_adh; 265 if (fTmp > (-k_adh * deltan_ce)){ 266 Fn_contact = fTmp; 267 }else{ // cohesion part 268 269 if (deltan > history[1]){ 270 count_flag = 1; 271 deltan_p_max = calculate_deltan_p_max (deltan_p, &sidata.contact_history[history_offset], count_flag, f_0, fTmp, k2, dex, dex_i, deltan, k_adh); 272 goto temp_calc; 273 } 274 Fn_contact = -k_adh * deltan_ce; 275 } 276 } 277 history[1] = deltan; 278 }else{ 279 const double k_adh = kn2kc[itype][jtype] * kn ; 280 kc_history[0] = k_adh; 281 282 double deltaMax; // the 4th value of the history array is deltaMax 283 if (deltan > history[0]) { 284 history[0] = deltan; 285 deltaMax = deltan; 286 }else{ 287 deltaMax = history[0]; 288 } 289 290 deltan_e = pow(deltan, dex); 291 deltan_ce = pow(deltan, cex); 292 deltaMax_e = pow(deltaMax, dex); 293 294 const double fTmp = k2*(deltan_e-deltaMax_e)+k1*deltaMax_e; 295 296 if (fTmp >= k1 * deltan_e) // loading 297 { 298 Fn_contact = k1 * deltan_e; 299 } else { 300 if (fTmp > (-k_adh * deltan_ce)) 301 { 302 Fn_contact = fTmp; 303 } else { // cohesion part 304 Fn_contact = -k_adh * deltan_ce; 305 const double newDeltaMax = ((k2 + k_adh)/(k2-k1))*deltan; 306 history[0] = newDeltaMax; 307 } 308 } 309 } 310 311 const double Fn_damping = -gamman*sidata.vn; 312 double Fn = Fn_damping + Fn_contact + f_0; 313 fo_history[0]=f_0; 314 315 if(limitForce && (Fn<0.0) && k_adh == 0 && f_0 == 0.0) 316 { 317 Fn = 0.0; 318 } 319 sidata.Fn = Fn; 320 sidata.kn = kn; 321 sidata.kt = kt; 322 sidata.gamman = gamman; 323 sidata.gammat = gammat; 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 } else { 332 i_forces.delta_F[0] = sidata.Fn * sidata.en[0]; 333 i_forces.delta_F[1] = sidata.Fn * sidata.en[1]; 334 i_forces.delta_F[2] = sidata.Fn * sidata.en[2]; 335 336 j_forces.delta_F[0] = -i_forces.delta_F[0]; 337 j_forces.delta_F[1] = -i_forces.delta_F[1]; 338 j_forces.delta_F[2] = -i_forces.delta_F[2]; 339 } 340 } 341 342 void surfacesClose(SurfacesCloseData & scdata, ForceData&, ForceData&) 343 { 344 if(scdata.contact_flags) *scdata.contact_flags &= ~CONTACT_NORMAL_MODEL; 345 double * const history = &scdata.contact_history[history_offset]; 346 history[0] = 0.0; 347 history[1] = 0.0; 348 } 349 void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 350 void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 351 352 protected: 353 double **Yeff; 354 double **Geff; 355 double **CoeffRestLog; 356 double **betaeff; 357 double **kn2kc; 358 double **kn2k1; 359 double cex; 360 double dex; 361 double **f_adh; 362 double **gamma_surf; 363 int history_offset; 364 int kc_offset; 365 int fo_offset; 366 bool tangential_damping; 367 bool limitForce; 368 bool fixKc; 369 }; 370 371 } 372 373 } 374 #endif 375 #endif 376