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