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 Christoph Kloss (DCS Computing GmbH, Linz) 37 Christoph Kloss (JKU Linz) 38 Richard Berger (JKU Linz) 39 40 Copyright 2012- DCS Computing GmbH, Linz 41 Copyright 2009-2012 JKU Linz 42 ------------------------------------------------------------------------- */ 43 44 #ifdef NORMAL_MODEL 45 NORMAL_MODEL(HERTZ_STIFFNESS,hertz/stiffness,4) 46 #else 47 #ifndef NORMAL_MODEL_HERTZ_STIFFNESS_H_ 48 #define NORMAL_MODEL_HERTZ_STIFFNESS_H_ 49 #include "contact_models.h" 50 #include "normal_model_base.h" 51 #include "global_properties.h" 52 #include <cmath> 53 54 namespace LIGGGHTS { 55 namespace ContactModels 56 { 57 template<> 58 class NormalModel<HERTZ_STIFFNESS> : public NormalModelBase 59 { 60 public: 61 NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup, class ContactModelBase * c) : 62 NormalModelBase(lmp, hsetup, c), 63 k_n(NULL), 64 k_t(NULL), 65 gamma_n(NULL), 66 gamma_t(NULL), 67 tangential_damping(false), 68 limitForce(false), 69 displayedSettings(false) 70 { 71 72 } 73 74 void registerSettings(Settings & settings) 75 { 76 settings.registerOnOff("tangential_damping", tangential_damping, true); 77 settings.registerOnOff("limitForce", limitForce); 78 } 79 80 inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) {} 81 82 void connectToProperties(PropertyRegistry & registry) 83 { 84 registry.registerProperty("k_n", &MODEL_PARAMS::createKn); 85 registry.registerProperty("k_t", &MODEL_PARAMS::createKt); 86 registry.registerProperty("gamma_n", &MODEL_PARAMS::createGamman); 87 registry.registerProperty("gamma_t", &MODEL_PARAMS::createGammat); 88 89 registry.connect("k_n", k_n,"model hertz/stiffness"); 90 registry.connect("k_t", k_t,"model hertz/stiffness"); 91 registry.connect("gamma_n", gamma_n,"model hertz/stiffness"); 92 registry.connect("gamma_t", gamma_t,"model hertz/stiffness"); 93 94 // error checks on coarsegraining 95 if(force->cg_active()) 96 error->cg(FLERR,"model hertz/stiffness"); 97 } 98 99 // effective exponent for stress-strain relationship 100 101 inline double stressStrainExponent() 102 { 103 return 1.5; 104 } 105 106 inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) 107 { 108 const int itype = sidata.itype; 109 const int jtype = sidata.jtype; 110 const double meff = sidata.meff; 111 double reff = sidata.is_wall ? sidata.radi : (sidata.radi*sidata.radj/(sidata.radi+sidata.radj)); 112 #ifdef SUPERQUADRIC_ACTIVE_FLAG 113 if(sidata.is_non_spherical && atom->superquadric_flag) { 114 reff = sidata.reff; 115 } 116 #endif 117 118 const double polyhertz = sqrt(reff*sidata.deltan); 119 double kn = polyhertz*k_n[itype][jtype]; 120 double kt = polyhertz*k_t[itype][jtype]; 121 const double gamman = polyhertz*meff*gamma_n[itype][jtype]; 122 const double gammat = tangential_damping ? polyhertz*meff*gamma_t[itype][jtype] : 0.0; 123 124 if(!displayedSettings) 125 { 126 displayedSettings = true; 127 128 /* 129 if(limitForce) 130 if(0 == comm->me) fprintf(screen," NormalModel<HERTZ_STIFFNESS>: will limit normal force.\n"); 131 */ 132 } 133 // convert Kn and Kt from pressure units to force/distance^2 134 kn /= force->nktv2p; 135 kt /= force->nktv2p; 136 137 const double Fn_damping = -gamman*sidata.vn; 138 const double Fn_contact = kn*sidata.deltan; 139 double Fn = Fn_damping + Fn_contact; 140 141 //limit force to avoid the artefact of negative repulsion force 142 if(limitForce && (Fn<0.0) ) 143 { 144 Fn = 0.0; 145 } 146 147 sidata.Fn = Fn; 148 sidata.kn = kn; 149 sidata.kt = kt; 150 sidata.gamman = gamman; 151 sidata.gammat = gammat; 152 153 #ifdef NONSPHERICAL_ACTIVE_FLAG 154 double torque_i[3] = {0.0, 0.0, 0.0}; //initialized here with zeros to avoid compiler warnings 155 double Fn_i[3] = { Fn * sidata.en[0], Fn * sidata.en[1], Fn * sidata.en[2]}; 156 if(sidata.is_non_spherical) { 157 double xci[3]; 158 vectorSubtract3D(sidata.contact_point, atom->x[sidata.i], xci); 159 vectorCross3D(xci, Fn_i, torque_i); 160 } 161 #endif 162 // apply normal force 163 if(sidata.is_wall) { 164 const double Fn_ = Fn * sidata.area_ratio; 165 i_forces.delta_F[0] += Fn_ * sidata.en[0]; 166 i_forces.delta_F[1] += Fn_ * sidata.en[1]; 167 i_forces.delta_F[2] += Fn_ * sidata.en[2]; 168 #ifdef NONSPHERICAL_ACTIVE_FLAG 169 if(sidata.is_non_spherical) { 170 //for non-spherical particles normal force can produce torque! 171 i_forces.delta_torque[0] += torque_i[0]; 172 i_forces.delta_torque[1] += torque_i[1]; 173 i_forces.delta_torque[2] += torque_i[2]; 174 } 175 #endif 176 } else { 177 i_forces.delta_F[0] += sidata.Fn * sidata.en[0]; 178 i_forces.delta_F[1] += sidata.Fn * sidata.en[1]; 179 i_forces.delta_F[2] += sidata.Fn * sidata.en[2]; 180 181 j_forces.delta_F[0] += -i_forces.delta_F[0]; 182 j_forces.delta_F[1] += -i_forces.delta_F[1]; 183 j_forces.delta_F[2] += -i_forces.delta_F[2]; 184 #ifdef NONSPHERICAL_ACTIVE_FLAG 185 if(sidata.is_non_spherical) { 186 //for non-spherical particles normal force can produce torque! 187 double xcj[3], torque_j[3]; 188 double Fn_j[3] = { -Fn_i[0], -Fn_i[1], -Fn_i[2]}; 189 vectorSubtract3D(sidata.contact_point, atom->x[sidata.j], xcj); 190 vectorCross3D(xcj, Fn_j, torque_j); 191 192 i_forces.delta_torque[0] += torque_i[0]; 193 i_forces.delta_torque[1] += torque_i[1]; 194 i_forces.delta_torque[2] += torque_i[2]; 195 196 j_forces.delta_torque[0] += torque_j[0]; 197 j_forces.delta_torque[1] += torque_j[1]; 198 j_forces.delta_torque[2] += torque_j[2]; 199 } 200 #endif 201 } 202 } 203 204 void surfacesClose(SurfacesCloseData&, ForceData&, ForceData&){} 205 void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 206 void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 207 208 protected: 209 double ** k_n; 210 double ** k_t; 211 double ** gamma_n; 212 double ** gamma_t; 213 214 bool tangential_damping; 215 bool limitForce; 216 bool displayedSettings; 217 }; 218 } 219 } 220 #endif // NORMAL_MODEL_HERTZ_STIFFNESS_H_ 221 #endif 222