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