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