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 ------------------------------------------------------------------------- */ 38 39 #ifdef NORMAL_MODEL 40 NORMAL_MODEL(LUDING,luding,12) 41 #else 42 #ifndef NORMAL_MODEL_LUDING_H_ 43 #define NORMAL_MODEL_LUDING_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 52 namespace LIGGGHTS { 53 namespace ContactModels 54 { 55 template<> 56 class NormalModel<LUDING> : public NormalModelBase 57 { 58 public: 59 NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup,class ContactModelBase *c) : 60 NormalModelBase(lmp, hsetup, c), 61 K_elastic(NULL), 62 CoeffRestLog(NULL), 63 kn2k1(NULL), 64 kn2kc(NULL), 65 phiF(NULL), 66 f_adh(NULL), 67 limitForce(false) 68 { 69 history_offset = hsetup->add_history_value("deltaMax", "0"); 70 kc_offset = hsetup->add_history_value("kc", "1"); 71 fo_offset = hsetup->add_history_value("fo", "1"); 72 c->add_history_offset("kc_offset", kc_offset); 73 c->add_history_offset("fo_offset", fo_offset); 74 } 75 76 inline void registerSettings(Settings & settings){ 77 settings.registerOnOff("tangential_damping", tangential_damping, true); 78 settings.registerOnOff("limitForce", limitForce, true); 79 } 80 inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) {} 81 82 inline void connectToProperties(PropertyRegistry & registry) { 83 registry.registerProperty("K_elastic", &MODEL_PARAMS::createLoadingStiffness,"model luding"); 84 registry.registerProperty("CoeffRestLog", &MODEL_PARAMS::createCoeffRestLog,"model luding"); 85 registry.registerProperty("kn2k1", &MODEL_PARAMS::createUnloadingStiffness,"model luding"); 86 registry.registerProperty("kn2kc", &MODEL_PARAMS::createCoeffAdhesionStiffness,"model luding"); 87 registry.registerProperty("phiF", &MODEL_PARAMS::createCoeffPlasticityDepth,"model luding"); 88 registry.registerProperty("f_adh", &MODEL_PARAMS::createPullOffForce,"model luding"); 89 90 registry.connect("K_elastic", K_elastic,"model luding"); 91 registry.connect("CoeffRestLog", CoeffRestLog,"model luding"); 92 registry.connect("kn2kc", kn2kc,"model luding"); 93 registry.connect("kn2k1", kn2k1, "model luding"); 94 registry.connect("phiF", phiF,"model luding"); 95 registry.connect("f_adh", f_adh,"model luding"); 96 97 // error checks on coarsegraining 98 if(force->cg_active()) 99 error->cg(FLERR,"model luding"); 100 } 101 102 // effective exponent for stress-strain relationship 103 104 inline double stressStrainExponent() 105 { 106 return 1.; 107 } 108 109 inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) 110 { 111 const int itype = sidata.itype; 112 const int jtype = sidata.jtype; 113 const double deltan = sidata.deltan; 114 double ri = sidata.radi; 115 double rj = sidata.radj; 116 double reff=sidata.is_wall ? sidata.radi : (ri*rj/(ri+rj)); 117 #ifdef SUPERQUADRIC_ACTIVE_FLAG 118 if(sidata.is_non_spherical && atom->superquadric_flag) 119 reff = sidata.reff; 120 #endif 121 double meff=sidata.meff; 122 double kn = K_elastic[itype][jtype]; 123 double kt = kn; 124 125 // convert Kn and Kt from pressure units to force/distance^2 126 kn /= force->nktv2p; 127 kt /= force->nktv2p; 128 129 const double k1 = kn; 130 const double k2Max = kn * kn2k1[itype][jtype]; 131 132 const double kc = kn2kc[itype][jtype] * kn; 133 const double f_0 = f_adh[itype][jtype]; 134 135 double gamman, gammat; 136 137 gamman = sqrt(4.*meff*kn/(1.+(M_PI/CoeffRestLog[itype][jtype])*(M_PI/CoeffRestLog[itype][jtype]))); 138 gammat = sqrt(4.*meff*kn/(1.+(M_PI/CoeffRestLog[itype][jtype])*(M_PI/CoeffRestLog[itype][jtype]))); 139 140 if (!tangential_damping) gammat = 0.0; 141 142 // get the history value -- maximal overlap 143 if(sidata.contact_flags) *sidata.contact_flags |= CONTACT_NORMAL_MODEL; 144 double * const history = &sidata.contact_history[history_offset]; 145 double * const kc_history = &sidata.contact_history[kc_offset]; 146 double * const fo_history = &sidata.contact_history[fo_offset]; 147 double deltaMax; // the 4th value of the history array is deltaMax 148 if (deltan > history[0]) { 149 history[0] = deltan; 150 deltaMax = deltan; 151 } else{ 152 deltaMax = history[0]; 153 } 154 155 // k2 dependent on the maximum overlap 156 // this accounts for an increasing stiffness with deformation - to capture nonlinearity 157 const double deltaMaxLim =(k2Max/(k2Max-k1))*phiF[itype][jtype]*2*reff; 158 159 double k2, fHys; 160 161 if (deltaMax >= deltaMaxLim) // big overlap ... no kn at all 162 { 163 k2 = k2Max; 164 const double fTmp = k2*(deltan-deltaMaxLim)+k1*deltaMaxLim;//k2*(deltan-delta0); 165 if (fTmp >= -kc*deltan) { // un-/reloading part (k2) 166 fHys = fTmp; 167 } else { // cohesion part 168 fHys = -kc*deltan; 169 const double newDeltaMax = ((k2 + kc)/(k2-k1))*deltan; 170 history[0] = newDeltaMax; 171 } 172 } else { 173 k2 = k1 + (k2Max - k1) * deltaMax/deltaMaxLim; 174 const double fTmp = k2*(deltan-deltaMax)+k1*deltaMax;//k2*(deltan-delta0); 175 if (fTmp >= k1*deltan) { // loading part (k1) 176 fHys = k1*deltan; 177 } else { // un-/reloading part (k2) 178 if (fTmp > -kc*deltan) { 179 fHys = fTmp; 180 } else { // cohesion part 181 fHys = -kc*deltan; 182 const double newDeltaMax = ((k2 + kc)/(k2-k1))*deltan; 183 history[0] = newDeltaMax; 184 } 185 } 186 } 187 188 const double Fn_damping = -gamman*sidata.vn; 189 double Fn = fHys + Fn_damping + f_0; 190 191 if(limitForce && (Fn<0.0) && kc == 0 && f_0 == 0.0){ 192 Fn = 0.0; 193 } 194 sidata.Fn = Fn; 195 sidata.kn = kn; 196 sidata.kt = kt; 197 kc_history[0] = kc; 198 fo_history[0] = f_0; 199 sidata.gamman = gamman; 200 sidata.gammat = gammat; 201 202 #ifdef NONSPHERICAL_ACTIVE_FLAG 203 double Fn_i[3] = { Fn * sidata.en[0], Fn * sidata.en[1], Fn * sidata.en[2]}; 204 double torque_i[3] = {0.0, 0.0, 0.0}; //initialized here with zeros to avoid compiler warnings 205 if(sidata.is_non_spherical) { 206 double xci[3]; 207 vectorSubtract3D(sidata.contact_point, atom->x[sidata.i], xci); 208 vectorCross3D(xci, Fn_i, torque_i); 209 } 210 #endif 211 212 // apply normal force 213 if(sidata.is_wall) { 214 const double Fn_ = Fn * sidata.area_ratio; 215 i_forces.delta_F[0] = Fn_ * sidata.en[0]; 216 i_forces.delta_F[1] = Fn_ * sidata.en[1]; 217 i_forces.delta_F[2] = Fn_ * sidata.en[2]; 218 #ifdef NONSPHERICAL_ACTIVE_FLAG 219 if(sidata.is_non_spherical) { 220 //for non-spherical particles normal force can produce torque! 221 i_forces.delta_torque[0] += torque_i[0]; 222 i_forces.delta_torque[1] += torque_i[1]; 223 i_forces.delta_torque[2] += torque_i[2]; 224 } 225 #endif 226 } else { 227 i_forces.delta_F[0] = sidata.Fn * sidata.en[0]; 228 i_forces.delta_F[1] = sidata.Fn * sidata.en[1]; 229 i_forces.delta_F[2] = sidata.Fn * sidata.en[2]; 230 231 j_forces.delta_F[0] = -i_forces.delta_F[0]; 232 j_forces.delta_F[1] = -i_forces.delta_F[1]; 233 j_forces.delta_F[2] = -i_forces.delta_F[2]; 234 #ifdef NONSPHERICAL_ACTIVE_FLAG 235 if(sidata.is_non_spherical) { 236 //for non-spherical particles normal force can produce torque! 237 double xcj[3], torque_j[3]; 238 double Fn_j[3] = { -Fn_i[0], -Fn_i[1], -Fn_i[2]}; 239 vectorSubtract3D(sidata.contact_point, atom->x[sidata.j], xcj); 240 vectorCross3D(xcj, Fn_j, torque_j); 241 242 i_forces.delta_torque[0] += torque_i[0]; 243 i_forces.delta_torque[1] += torque_i[1]; 244 i_forces.delta_torque[2] += torque_i[2]; 245 246 j_forces.delta_torque[0] += torque_j[0]; 247 j_forces.delta_torque[1] += torque_j[1]; 248 j_forces.delta_torque[2] += torque_j[2]; 249 } 250 #endif 251 } 252 } 253 254 inline void surfacesClose(SurfacesCloseData & scdata, ForceData&, ForceData&) 255 { 256 if(scdata.contact_flags) *scdata.contact_flags &= ~CONTACT_NORMAL_MODEL; 257 double * const history = &scdata.contact_history[history_offset]; 258 history[0] = 0.0; 259 } 260 261 void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 262 void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 263 264 protected: 265 double **K_elastic; 266 double **CoeffRestLog; 267 double **kn2k1; 268 double **kn2kc; 269 double **phiF; 270 double **f_adh; 271 272 int history_offset; 273 int kc_offset; 274 int fo_offset; 275 276 bool tangential_damping; 277 bool limitForce; 278 }; 279 } 280 } 281 #endif // NORMAL_MODEL_LUDING_H_ 282 #endif 283