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