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    Tomaz M. Zorec (University of Ljubljana)
39 ------------------------------------------------------------------------- */
40 #ifdef NORMAL_MODEL
41 NORMAL_MODEL(EDINBURGH,edinburgh,10)
42 #else
43 #ifndef NORMAL_MODEL_EDINBURGH_H_
44 #define NORMAL_MODEL_EDINBURGH_H_
45 #include "contact_models.h"
46 #include "normal_model_base.h"
47 #include <cmath>
48 #include "atom.h"
49 #include "force.h"
50 #include "update.h"
51 #include "global_properties.h"
52 
53 namespace LIGGGHTS {
54 namespace ContactModels
55 {
56   template<>
57   class NormalModel<EDINBURGH> : public NormalModelBase
58   {
59   public:
60     NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup,class ContactModelBase *c) :
61       NormalModelBase(lmp, hsetup, c),
62       Yeff(NULL),
63       Geff(NULL),
64       CoeffRestLog(NULL),
65       betaeff(NULL),
66       kn2kc(NULL),
67       kn2k1(NULL),
68       cex(0.0),
69       dex(0.0),
70       f_adh(NULL),
71       gamma_surf(NULL),
72       history_offset(-1),
73       kc_offset(-1),
74       fo_offset(-1),
75       tangential_damping(false),
76       limitForce(false),
77       fixKc(false)
78     {
79       history_offset = hsetup->add_history_value("deltaMax", "1");
80       hsetup->add_history_value("old_delta", "1");
81       kc_offset = hsetup->add_history_value("kc", "1");
82       fo_offset = hsetup->add_history_value("fo", "1");
83       c->add_history_offset("kc_offset", kc_offset);
84       c->add_history_offset("fo_offset", fo_offset);
85     }
86 
87     void registerSettings(Settings & settings)
88     {
89       settings.registerOnOff("tangential_damping", tangential_damping, true);
90       settings.registerOnOff("limitForce", limitForce,true);
91       settings.registerOnOff("fixKc", fixKc);
92       //TODO error->one(FLERR,"TODO here also check if right surface model used");
93     }
94 
95 	inline void postSettings(IContactHistorySetup *, ContactModelBase *) {}
96 
97     void connectToProperties(PropertyRegistry & registry) {
98 
99       registry.registerProperty("Yeff", &MODEL_PARAMS::createYeff,"model edinburgh");                  // only used for non-linear model
100       registry.registerProperty("Geff", &MODEL_PARAMS::createGeff,"model edinburgh");                  // only used for non-linear model
101       registry.registerProperty("CoeffRestLog", &MODEL_PARAMS::createCoeffRestLog, "model edinburgh");
102       registry.registerProperty("betaeff", &MODEL_PARAMS::createBetaEff,"model edinburgh");
103       registry.registerProperty("kn2kc", &MODEL_PARAMS::createCoeffAdhesionStiffness,"model edinburgh");
104       registry.registerProperty("kn2k1", &MODEL_PARAMS::createUnloadingStiffness, "model edinburgh");
105       registry.registerProperty("cex", &MODEL_PARAMS::createAdhesionExponent, "model edinburgh");
106       registry.registerProperty("dex", &MODEL_PARAMS::createOverlapExponent, "model edinburgh");
107       registry.registerProperty("f_adh", &MODEL_PARAMS::createPullOffForce, "model edinburgh");
108       registry.registerProperty("gamma_surf", &MODEL_PARAMS::createSurfaceEnergy, "model edinburgh");
109 
110       registry.connect("Yeff", Yeff,"model edinburgh");
111       registry.connect("Geff", Geff,"model edinburgh");
112       registry.connect("CoeffRestLog", CoeffRestLog, "model edinburgh");
113       registry.connect("betaeff", betaeff,"model edinburgh");
114       registry.connect("kn2kc", kn2kc,"model edinburgh");
115       registry.connect("kn2k1", kn2k1,"model edinburgh");
116       registry.connect("cex", cex,"model edinburgh");
117       registry.connect("dex", dex,"model edinburgh");
118       registry.connect("f_adh", f_adh,"model edinburgh");
119       registry.connect("gamma_surf", gamma_surf,"model edinburgh");
120 
121     }
122 
123     // effective exponent for stress-strain relationship
124 
125     inline double stressStrainExponent()
126     {
127       return 1.5;
128     }
129 
130     inline double calculate_deltan_p_max (double deltan_p, double * const history, int count_flag, const double f_0, double fTmp, const double k2, double dex, double dex_i, double deltan, double k_adh)
131     {
132       //calculating the maximum overlap
133       double deltan_p_max;
134 
135       if (count_flag == 0 ) {
136            if (deltan_p > history[0]) {
137             deltan_p_max = deltan_p;
138             history[0] = deltan_p_max;
139            } else {
140             deltan_p_max = history[0];
141           }
142       } else {
143         deltan_p_max = pow((pow(history[1], dex) + ((k_adh/k2)*pow(history[1], cex))), dex_i);
144         history[0] = deltan_p_max;
145       }
146       return history[0];
147     }
148 
149     inline double whichd(double fTmp, double k1, double deltan_e, double d1, double d2)
150     {
151       //calculation distance between the particles based on branch
152       if (fTmp >= k1 * deltan_e){
153         return d1;
154       } else{
155         return d2;
156       }
157     }
158 
159     inline double calculate_k_adh(double d, double risq, double rjsq, const double g_surf, double f_min_lim, double deltan_pe_max, double deltan_p_max, const double k2, const double dex_i, const double cex)
160     {
161       double new_k_c, delta_min, a, f_min, dsq;
162        dsq = d*d;
163        a = (1./(2.*d)) * sqrt(4* dsq * risq - ((dsq - rjsq + risq)*(dsq - rjsq + risq)));
164        f_min = ( 1.5 * M_PI * g_surf * a);
165 
166       if (f_min > f_min_lim)
167       {
168         f_min = f_min_lim;
169         delta_min = 0.5 * deltan_p_max;
170       } else {
171         delta_min = pow((-f_min + k2 * deltan_pe_max)/k2, dex_i);
172       }
173       new_k_c = f_min/pow(delta_min, cex);
174       return new_k_c;
175     }
176 
177     inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces)
178     {
179       const int itype = sidata.itype;
180       const int jtype = sidata.jtype;
181       const double deltan = sidata.deltan;
182       double ri = sidata.radi;
183       double rj = sidata.radj;
184       double reff=sidata.is_wall ? sidata.radi : ((ri*rj)/(ri+rj));
185 #ifdef SUPERQUADRIC_ACTIVE_FLAG
186         if(sidata.is_non_spherical && atom->superquadric_flag)
187             reff = sidata.reff;
188 #endif
189       double meff=sidata.meff;
190       const double f_0 = f_adh[itype][jtype];
191       double sqrtval = sqrt(reff*sidata.deltan);
192       double Sn=2.*Yeff[itype][jtype]*sqrtval;
193       double St=8.*Geff[itype][jtype]*sqrtval;
194       double kn, kt;
195 
196       if(dex==1){
197         kn = 2.0*Yeff[itype][jtype]*reff;
198         kt = Sn;
199       } else {
200         kn = (4./3.)*Yeff[itype][jtype]*sqrt(reff);
201         kt = St;
202         }
203 
204       // convert Kn and Kt from pressure units to force/distance^2
205       kn /= force->nktv2p;
206       kt /= force->nktv2p;
207 
208       const double dex_i = 1./dex;
209       const double k1 = kn;
210       const double k2 = kn*kn2k1[itype][jtype];
211 
212       double gamman, gammat;
213       gamman = sqrt(4.*meff*kn/(1.+(M_PI/CoeffRestLog[itype][jtype])*(M_PI/CoeffRestLog[itype][jtype])));
214       gammat = sqrt(4.*meff*kn/(1.+(M_PI/CoeffRestLog[itype][jtype])*(M_PI/CoeffRestLog[itype][jtype])));
215 
216       if(sidata.contact_flags) *sidata.contact_flags |= CONTACT_NORMAL_MODEL;
217       double * const history = &sidata.contact_history[history_offset];
218       double * const kc_history = &sidata.contact_history[kc_offset];
219       double * const fo_history = &sidata.contact_history[fo_offset];
220 
221       if (!tangential_damping) gammat = 0.0;
222 
223       double Fn_contact, deltan_e, deltan_ce, deltaMax_e;
224       double k_adh = 0.0;
225 
226       if (fixKc == false) {
227         const double lambda = pow((1. - k1/k2), dex_i);
228         int count_flag = 0;
229 
230         // get the history value -- maximal overlap
231 
232         double deltan_p = lambda*deltan;
233         double deltan_p_max, deltan_pe_max;
234 
235         deltan_e = pow(deltan, dex);
236         deltan_ce = pow(deltan, cex);
237 
238         double fTmp = k2*(deltan_e - history[0]);
239         double  d, d1, d2, r_sum;
240 
241         deltan_p_max = calculate_deltan_p_max (deltan_p, &sidata.contact_history[history_offset], count_flag, f_0, fTmp, k2, dex, dex_i, deltan, k_adh);
242 
243         // Normal force calculation for Edinburgh model
244 
245         const double g_surf = gamma_surf[itype][jtype];
246 
247         r_sum = sidata.radsum;
248         d1 = sidata.is_wall ? ri : r_sum - deltan;
249         d2 = sidata.is_wall ? ri  : r_sum - deltan_p_max;
250         double risq = ri*ri;
251         double rjsq = rj*rj;
252 
253         temp_calc:
254 
255         deltan_pe_max = pow(deltan_p_max, dex);
256         const double f_min_lim = (k2 * deltan_pe_max) * 0.5;
257         fTmp=k2*(deltan_e - deltan_pe_max);
258 
259         if (fTmp >= k1 * deltan_e){ // loading
260           Fn_contact = k1 * deltan_e;
261         }else{
262           d = whichd(fTmp, k1, deltan_e, d1, d2);
263           k_adh = calculate_k_adh(d, risq, rjsq, g_surf, f_min_lim, deltan_pe_max, deltan_p_max, k2, dex_i, cex);
264           kc_history[0] = k_adh;
265          if (fTmp > (-k_adh * deltan_ce)){
266             Fn_contact = fTmp;
267           }else{  // cohesion part
268 
269             if (deltan > history[1]){
270               count_flag = 1;
271               deltan_p_max = calculate_deltan_p_max (deltan_p, &sidata.contact_history[history_offset], count_flag, f_0, fTmp, k2, dex, dex_i, deltan, k_adh);
272               goto temp_calc;
273             }
274             Fn_contact = -k_adh * deltan_ce;
275           }
276         }
277         history[1] = deltan;
278       }else{
279         const double k_adh = kn2kc[itype][jtype] * kn ;
280         kc_history[0] = k_adh;
281 
282         double deltaMax; // the 4th value of the history array is deltaMax
283         if (deltan > history[0]) {
284             history[0] = deltan;
285             deltaMax = deltan;
286         }else{
287           deltaMax = history[0];
288         }
289 
290         deltan_e = pow(deltan, dex);
291         deltan_ce = pow(deltan, cex);
292         deltaMax_e = pow(deltaMax, dex);
293 
294         const double fTmp = k2*(deltan_e-deltaMax_e)+k1*deltaMax_e;
295 
296         if (fTmp >= k1 * deltan_e) // loading
297         {
298           Fn_contact = k1 * deltan_e;
299         } else {
300          if (fTmp > (-k_adh * deltan_ce))
301           {
302             Fn_contact = fTmp;
303           } else { // cohesion part
304             Fn_contact = -k_adh * deltan_ce;
305             const double newDeltaMax = ((k2 + k_adh)/(k2-k1))*deltan;
306             history[0] = newDeltaMax;
307           }
308         }
309       }
310 
311       const double Fn_damping = -gamman*sidata.vn;
312       double Fn = Fn_damping + Fn_contact + f_0;
313       fo_history[0]=f_0;
314 
315       if(limitForce && (Fn<0.0) && k_adh == 0 && f_0 == 0.0)
316       {
317         Fn = 0.0;
318       }
319       sidata.Fn = Fn;
320       sidata.kn = kn;
321       sidata.kt = kt;
322       sidata.gamman = gamman;
323       sidata.gammat = gammat;
324 
325       // apply normal force
326       if(sidata.is_wall) {
327         const double Fn_ = Fn * sidata.area_ratio;
328         i_forces.delta_F[0] = Fn_ * sidata.en[0];
329         i_forces.delta_F[1] = Fn_ * sidata.en[1];
330         i_forces.delta_F[2] = Fn_ * sidata.en[2];
331       } else {
332         i_forces.delta_F[0] = sidata.Fn * sidata.en[0];
333         i_forces.delta_F[1] = sidata.Fn * sidata.en[1];
334         i_forces.delta_F[2] = sidata.Fn * sidata.en[2];
335 
336         j_forces.delta_F[0] = -i_forces.delta_F[0];
337         j_forces.delta_F[1] = -i_forces.delta_F[1];
338         j_forces.delta_F[2] = -i_forces.delta_F[2];
339       }
340     }
341 
342     void surfacesClose(SurfacesCloseData & scdata, ForceData&, ForceData&)
343     {
344       if(scdata.contact_flags) *scdata.contact_flags &= ~CONTACT_NORMAL_MODEL;
345       double * const history = &scdata.contact_history[history_offset];
346       history[0] = 0.0;
347       history[1] = 0.0;
348     }
349     void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){}
350     void endPass(SurfacesIntersectData&, ForceData&, ForceData&){}
351 
352   protected:
353     double **Yeff;
354     double **Geff;
355     double **CoeffRestLog;
356     double **betaeff;
357     double **kn2kc;
358     double **kn2k1;
359     double cex;
360     double dex;
361     double **f_adh;
362     double **gamma_surf;
363     int history_offset;
364     int kc_offset;
365     int fo_offset;
366     bool tangential_damping;
367     bool limitForce;
368     bool fixKc;
369   };
370 
371 }
372 
373 }
374 #endif
375 #endif
376