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 ------------------------------------------------------------------------- */ 39 40 #ifdef TANGENTIAL_MODEL 41 TANGENTIAL_MODEL(TANGENTIAL_LUDING,tan_luding,4) 42 #else 43 #ifndef TANGENTIAL_MODEL_LUDING_H_ 44 #define TANGENTIAL_MODEL_LUDING_H_ 45 #include "contact_models.h" 46 #include "tangential_model_base.h" 47 #include <cmath> 48 #include "update.h" 49 #include "global_properties.h" 50 #include "atom.h" 51 52 namespace LIGGGHTS { 53 namespace ContactModels 54 { 55 template<> 56 class TangentialModel<TANGENTIAL_LUDING> : public TangentialModelBase 57 { 58 double ** coeffFrict; 59 double ** coeffMu; 60 double ** coeffFricVisc; 61 double ** kT2kcMax; 62 int history_offset; 63 int kc_offset; 64 int fo_offset; 65 66 public: 67 TangentialModel(LAMMPS * lmp, IContactHistorySetup * hsetup,class ContactModelBase *c) : TangentialModelBase(lmp, hsetup, c), 68 coeffFrict(NULL), 69 coeffMu(NULL), 70 coeffFricVisc(NULL), 71 kT2kcMax(NULL), 72 heating(false), 73 heating_track(false), 74 // kc_Stiffness(NULL), 75 cmb(c) 76 { 77 history_offset = hsetup->add_history_value("shearx", "1"); 78 hsetup->add_history_value("sheary", "1"); 79 hsetup->add_history_value("shearz", "1"); 80 kc_offset = cmb->get_history_offset("kc_offset"); 81 fo_offset = cmb->get_history_offset("fo_offset"); 82 } 83 84 inline void registerSettings(Settings& settings){ 85 settings.registerOnOff("heating_tangential_history",heating,false); 86 settings.registerOnOff("heating_tracking",heating_track,false); 87 //TODO error->one(FLERR,"TODO here also check if right surface model used"); 88 } 89 90 inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) {} 91 92 inline void connectToProperties(PropertyRegistry & registry) 93 { 94 registry.registerProperty("coeffFrict", &MODEL_PARAMS::createCoeffFrict); 95 registry.registerProperty("coeffMu", &MODEL_PARAMS::createCoeffMu); 96 registry.registerProperty("coeffFricVisc", &MODEL_PARAMS::createCoeffFricVisc); 97 registry.registerProperty("kT2kcMax", &MODEL_PARAMS::createCoeffFrictionStiffness); 98 99 registry.connect("coeffFrict", coeffFrict,"tangential_model tan_luding"); 100 registry.connect("coeffMu", coeffMu,"tangential_model tan_luding"); 101 registry.connect("coeffFricVisc", coeffFricVisc,"tangential_model tan_luding"); 102 registry.connect("kT2kcMax", kT2kcMax,"tangential_model tan_luding"); 103 } 104 105 inline void surfacesIntersect(const SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) 106 { 107 const double enx = sidata.en[0]; 108 const double eny = sidata.en[1]; 109 const double enz = sidata.en[2]; 110 const double deltan = sidata.deltan; 111 112 double k_t = sidata.kt; 113 double kt = k_t * kT2kcMax[sidata.itype][sidata.jtype]; // tangential stiffness based on the ratio input 114 // shear history effects 115 if(sidata.contact_flags) *sidata.contact_flags |= CONTACT_TANGENTIAL_MODEL; 116 double * const shear = &sidata.contact_history[history_offset]; 117 // double * const K_adh = &sidata.contact_history[kc_Stiffness]; 118 if (sidata.shearupdate && sidata.computeflag) { 119 const double dt = update->dt; 120 shear[0] += sidata.vtr1 * dt; 121 shear[1] += sidata.vtr2 * dt; 122 shear[2] += sidata.vtr3 * dt; 123 // rotate shear displacements 124 double rsht = shear[0]*enx + shear[1]*eny + shear[2]*enz; 125 shear[0] -= rsht * enx; 126 shear[1] -= rsht * eny; 127 shear[2] -= rsht * enz; 128 } 129 130 const double shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); 131 132 const double xmuS = coeffFrict[sidata.itype][sidata.jtype]; //coefficient of sliding friction 133 const double xmuD = xmuS * coeffMu[sidata.itype][sidata.jtype]; //coefficient of dynamic friction 134 const double gammat = sidata.gammat * coeffFricVisc[sidata.itype][sidata.jtype]; 135 136 // error->one(FLERR,"model-specific data not allows in contact_interface.h, do use contact history instead"); 137 const double kc = sidata.contact_history[kc_offset];//K_adh;//0.;//sidata.kc; 138 const double f_adh = sidata.contact_history[fo_offset]; 139 double Ft1 = -(kt * shear[0]); 140 double Ft2 = -(kt * shear[1]); 141 double Ft3 = -(kt * shear[2]); 142 143 // rescale frictional displacements and forces if needed 144 const double Ft_shear = kt * shrmag; 145 const double Ft_frictionS = xmuS * (fabs(sidata.Fn + kc * deltan - f_adh)); 146 147 // energy loss from sliding or damping 148 if (Ft_shear > Ft_frictionS) { 149 if (shrmag != 0.0) { 150 const double Ft_frictionD = xmuD * (fabs(sidata.Fn + kc * deltan - f_adh)); 151 const double ratio = Ft_frictionD/Ft_shear; 152 153 Ft1 *= ratio; 154 Ft2 *= ratio; 155 Ft3 *= ratio; 156 157 shear[0] = -(Ft1) /kt; 158 shear[1] = -(Ft2) /kt; 159 shear[2] = -(Ft3) /kt; 160 } else Ft1 = Ft2 = Ft3 = 0.0; 161 } else { 162 Ft1 -= (gammat*sidata.vtr1); 163 Ft2 -= (gammat*sidata.vtr2); 164 Ft3 -= (gammat*sidata.vtr3); 165 } 166 167 // forces & torques 168 const double tor1 = eny * Ft3 - enz * Ft2; 169 const double tor2 = enz * Ft1 - enx * Ft3; 170 const double tor3 = enx * Ft2 - eny * Ft1; 171 172 #ifdef NONSPHERICAL_ACTIVE_FLAG 173 double torque_i[3]; 174 if(sidata.is_non_spherical) { 175 double xci[3]; 176 double Ft_i[3] = { Ft1, Ft2, Ft3 }; 177 vectorSubtract3D(sidata.contact_point, atom->x[sidata.i], xci); 178 vectorCross3D(xci, Ft_i, torque_i); 179 } else { 180 torque_i[0] = -sidata.cri * tor1; 181 torque_i[1] = -sidata.cri * tor2; 182 torque_i[2] = -sidata.cri * tor3; 183 } 184 #endif 185 186 // return resulting forces 187 if(sidata.is_wall) { 188 const double area_ratio = sidata.area_ratio; 189 i_forces.delta_F[0] += Ft1 * area_ratio; 190 i_forces.delta_F[1] += Ft2 * area_ratio; 191 i_forces.delta_F[2] += Ft3 * area_ratio; 192 #ifdef NONSPHERICAL_ACTIVE_FLAG 193 i_forces.delta_torque[0] += torque_i[0] * area_ratio; 194 i_forces.delta_torque[1] += torque_i[1] * area_ratio; 195 i_forces.delta_torque[2] += torque_i[2] * area_ratio; 196 #else 197 i_forces.delta_torque[0] += -sidata.cri * tor1 * area_ratio; 198 i_forces.delta_torque[1] += -sidata.cri * tor2 * area_ratio; 199 i_forces.delta_torque[2] += -sidata.cri * tor3 * area_ratio; 200 #endif 201 } else { 202 i_forces.delta_F[0] += Ft1; 203 i_forces.delta_F[1] += Ft2; 204 i_forces.delta_F[2] += Ft3; 205 206 j_forces.delta_F[0] += -Ft1; 207 j_forces.delta_F[1] += -Ft2; 208 j_forces.delta_F[2] += -Ft3; 209 #ifdef NONSPHERICAL_ACTIVE_FLAG 210 double torque_j[3]; 211 if(sidata.is_non_spherical) { 212 double xcj[3]; 213 vectorSubtract3D(sidata.contact_point, atom->x[sidata.j], xcj); 214 double Ft_j[3] = { -Ft1, -Ft2, -Ft3 }; 215 vectorCross3D(xcj, Ft_j, torque_j); 216 } else { 217 torque_j[0] = -sidata.crj * tor1; 218 torque_j[1] = -sidata.crj * tor2; 219 torque_j[2] = -sidata.crj * tor3; 220 } 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 j_forces.delta_torque[0] += torque_j[0]; 226 j_forces.delta_torque[1] += torque_j[1]; 227 j_forces.delta_torque[2] += torque_j[2]; 228 #else 229 230 i_forces.delta_torque[0] += -sidata.cri * tor1; 231 i_forces.delta_torque[1] += -sidata.cri * tor2; 232 i_forces.delta_torque[2] += -sidata.cri * tor3; 233 234 j_forces.delta_torque[0] += -sidata.crj * tor1; 235 j_forces.delta_torque[1] += -sidata.crj * tor2; 236 j_forces.delta_torque[2] += -sidata.crj * tor3; 237 #endif 238 } 239 } 240 241 inline void surfacesClose(SurfacesCloseData & scdata, ForceData&, ForceData&) 242 { 243 // unset non-touching neighbors 244 // TODO even if shearupdate == false? 245 if(scdata.contact_flags) *scdata.contact_flags &= ~CONTACT_TANGENTIAL_MODEL; 246 double * const shear = &scdata.contact_history[history_offset]; 247 shear[0] = 0.0; 248 shear[1] = 0.0; 249 shear[2] = 0.0; 250 } 251 252 inline void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 253 inline void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 254 protected: 255 bool heating; 256 bool heating_track; 257 class ContactModelBase *cmb; 258 }; 259 } 260 } 261 #endif // TANGENTIAL_MODEL_LUDING_H_ 262 #endif 263