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 TANGENTIAL_MODEL 45 TANGENTIAL_MODEL(TANGENTIAL_NO_HISTORY,no_history,1) 46 #else 47 #ifndef TANGENTIAL_MODEL_NO_HISTORY_H_ 48 #define TANGENTIAL_MODEL_NO_HISTORY_H_ 49 #include "contact_models.h" 50 #include "tangential_model_base.h" 51 #include <algorithm> 52 #include <cmath> 53 #include "global_properties.h" 54 #include "fix_property_atom.h" 55 56 namespace LIGGGHTS { 57 namespace ContactModels 58 { 59 template<> 60 class TangentialModel<TANGENTIAL_NO_HISTORY> : public TangentialModelBase 61 { 62 public: 63 static const int MASK = CM_CONNECT_TO_PROPERTIES | CM_SURFACES_INTERSECT; 64 65 TangentialModel(LAMMPS * lmp, IContactHistorySetup *hsetup, class ContactModelBase *cmb) : 66 TangentialModelBase(lmp, hsetup, cmb), 67 coeffFrict(NULL), 68 disable_when_bonded_(false), 69 bond_history_offset_(-1), 70 dissipation_history_offset_(-1), 71 dissipatedflag_(false), 72 fix_dissipated_(NULL) 73 { 74 75 } 76 77 inline void registerSettings(Settings &settings) 78 { 79 settings.registerOnOff("disableTangentialWhenBonded", disable_when_bonded_, false); 80 settings.registerOnOff("computeDissipatedEnergy", dissipatedflag_, false); 81 } 82 83 inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) 84 { 85 if (disable_when_bonded_) 86 { 87 bond_history_offset_ = cmb->get_history_offset("bond_contactflag"); 88 if (bond_history_offset_ < 0) 89 error->one(FLERR, "Could not find bond history offset"); 90 } 91 if (dissipatedflag_) 92 { 93 if (cmb->is_wall()) 94 { 95 fix_dissipated_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("dissipated_energy_wall", "property/atom", "vector", 0, 0, "dissipated energy")); 96 dissipation_history_offset_ = cmb->get_history_offset("dissipation_force"); 97 if (!dissipation_history_offset_) 98 error->one(FLERR, "Internal error: Could not find dissipation history offset"); 99 } 100 else 101 fix_dissipated_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("dissipated_energy", "property/atom", "vector", 0, 0, "dissipated energy")); 102 if (!fix_dissipated_) 103 error->one(FLERR, "Surface model has not registered dissipated_energy fix"); 104 } 105 } 106 107 inline void connectToProperties(PropertyRegistry & registry){ 108 registry.registerProperty("coeffFrict", &MODEL_PARAMS::createCoeffFrict); 109 registry.connect("coeffFrict", coeffFrict,"tangential_model history"); 110 } 111 112 inline void surfacesIntersect(const SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces) { 113 if(sidata.contact_flags) *sidata.contact_flags |= CONTACT_TANGENTIAL_MODEL; 114 const double xmu = coeffFrict[sidata.itype][sidata.jtype]; 115 const double enx = sidata.en[0]; 116 const double eny = sidata.en[1]; 117 const double enz = sidata.en[2]; 118 const double vrel = sqrt(sidata.vtr1*sidata.vtr1 + sidata.vtr2*sidata.vtr2 + sidata.vtr3*sidata.vtr3); 119 120 // force normalization 121 const double Ft_friction = xmu * fabs(sidata.Fn); 122 double gamma = 0.0; 123 124 if (Ft_friction < sidata.gammat*vrel) 125 gamma = Ft_friction/vrel; 126 else 127 gamma = sidata.gammat; 128 129 // tangential force due to tangential velocity damping 130 131 const double Ft1 = -gamma*sidata.vtr1; 132 const double Ft2 = -gamma*sidata.vtr2; 133 const double Ft3 = -gamma*sidata.vtr3; 134 135 // forces & torques 136 137 const double tor1 = (eny*Ft3 - enz*Ft2); 138 const double tor2 = (enz*Ft1 - enx*Ft3); 139 const double tor3 = (enx*Ft2 - eny*Ft1); 140 141 #ifdef NONSPHERICAL_ACTIVE_FLAG 142 double torque_i[3]; 143 if(sidata.is_non_spherical) { 144 double xci[3]; 145 double Ft_i[3] = { Ft1, Ft2, Ft3 }; 146 vectorSubtract3D(sidata.contact_point, atom->x[sidata.i], xci); 147 vectorCross3D(xci, Ft_i, torque_i); 148 } else { 149 torque_i[0] = -sidata.cri * tor1; 150 torque_i[1] = -sidata.cri * tor2; 151 torque_i[2] = -sidata.cri * tor3; 152 } 153 #endif 154 // return resulting forces 155 if (!disable_when_bonded_ || sidata.contact_history[bond_history_offset_] < 0.5) 156 { 157 // energy balance terms 158 if (dissipatedflag_ && sidata.computeflag && sidata.shearupdate) 159 { 160 const double crj = sidata.is_wall ? 0.0 : sidata.crj; 161 // compute increment in dissipated energy 162 double * const * const dissipated = fix_dissipated_->array_atom; 163 double * const dissipated_i = dissipated[sidata.i]; 164 double * const dissipated_j = dissipated[sidata.j]; 165 dissipated_i[1] += -Ft1; 166 dissipated_i[2] += -Ft2; 167 dissipated_i[3] += -Ft3; 168 dissipated_i[4] += sidata.cri*tor1; 169 dissipated_i[5] += sidata.cri*tor2; 170 dissipated_i[6] += sidata.cri*tor3; 171 if (sidata.j < atom->nlocal && !sidata.is_wall) 172 { 173 dissipated_j[1] -= -Ft1; 174 dissipated_j[2] -= -Ft2; 175 dissipated_j[3] -= -Ft3; 176 dissipated_j[4] += crj*tor1; 177 dissipated_j[5] += crj*tor2; 178 dissipated_j[6] += crj*tor3; 179 } 180 else if (sidata.is_wall) 181 { 182 double * const diss_force = &sidata.contact_history[dissipation_history_offset_]; 183 diss_force[0] -= -Ft1; 184 diss_force[1] -= -Ft2; 185 diss_force[2] -= -Ft3; 186 } 187 } 188 if(sidata.is_wall) { 189 const double area_ratio = sidata.area_ratio; 190 i_forces.delta_F[0] += Ft1 * area_ratio; 191 i_forces.delta_F[1] += Ft2 * area_ratio; 192 i_forces.delta_F[2] += Ft3 * area_ratio; 193 #ifdef NONSPHERICAL_ACTIVE_FLAG 194 i_forces.delta_torque[0] += torque_i[0] * area_ratio; 195 i_forces.delta_torque[1] += torque_i[1] * area_ratio; 196 i_forces.delta_torque[2] += torque_i[2] * area_ratio; 197 #else 198 i_forces.delta_torque[0] += -sidata.cri * tor1 * area_ratio; 199 i_forces.delta_torque[1] += -sidata.cri * tor2 * area_ratio; 200 i_forces.delta_torque[2] += -sidata.cri * tor3 * area_ratio; 201 #endif 202 } else { 203 i_forces.delta_F[0] += Ft1; 204 i_forces.delta_F[1] += Ft2; 205 i_forces.delta_F[2] += Ft3; 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 i_forces.delta_torque[0] += -sidata.cri * tor1; 230 i_forces.delta_torque[1] += -sidata.cri * tor2; 231 i_forces.delta_torque[2] += -sidata.cri * tor3; 232 233 j_forces.delta_torque[0] += -sidata.crj * tor1; 234 j_forces.delta_torque[1] += -sidata.crj * tor2; 235 j_forces.delta_torque[2] += -sidata.crj * tor3; 236 #endif 237 } 238 } 239 } 240 241 inline void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){} 242 inline void endPass(SurfacesIntersectData&, ForceData&, ForceData&){} 243 244 inline void surfacesClose(SurfacesCloseData &scdata, ForceData&, ForceData&) 245 { 246 if(scdata.contact_flags) *scdata.contact_flags &= ~CONTACT_TANGENTIAL_MODEL; 247 } 248 249 private: 250 double ** coeffFrict; 251 bool disable_when_bonded_; 252 int bond_history_offset_; 253 int dissipation_history_offset_; 254 bool dissipatedflag_; 255 FixPropertyAtom *fix_dissipated_; 256 }; 257 } 258 } 259 #endif // TANGENTIAL_MODEL_NO_HISTORY_H_ 260 #endif 261