1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing authors:
17 Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL)
18 Leo Silbert (SNL), Gary Grest (SNL)
19 ----------------------------------------------------------------------- */
20
21 #include "pair_granular.h"
22
23 #include "atom.h"
24 #include "comm.h"
25 #include "error.h"
26 #include "fix.h"
27 #include "fix_dummy.h"
28 #include "fix_neigh_history.h"
29 #include "force.h"
30 #include "math_const.h"
31 #include "math_special.h"
32 #include "memory.h"
33 #include "modify.h"
34 #include "neigh_list.h"
35 #include "neigh_request.h"
36 #include "neighbor.h"
37 #include "update.h"
38
39 #include <cmath>
40 #include <cstring>
41
42 using namespace LAMMPS_NS;
43 using namespace MathConst;
44 using namespace MathSpecial;
45
46 #define PI27SQ 266.47931882941264802866 // 27*PI**2
47 #define THREEROOT3 5.19615242270663202362 // 3*sqrt(3)
48 #define SIXROOT6 14.69693845669906728801 // 6*sqrt(6)
49 #define INVROOT6 0.40824829046386307274 // 1/sqrt(6)
50 #define FOURTHIRDS 4.0/3.0 // 4/3
51 #define THREEQUARTERS 0.75 // 3/4
52
53 #define EPSILON 1e-10
54
55 enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR};
56 enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI};
57 enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY,
58 TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE,
59 TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE};
60 enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL};
61 enum {ROLL_NONE, ROLL_SDS};
62
63 /* ---------------------------------------------------------------------- */
64
PairGranular(LAMMPS * lmp)65 PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
66 {
67 single_enable = 1;
68 no_virial_fdotr_compute = 1;
69 centroidstressflag = CENTROID_NOTAVAIL;
70 finitecutflag = 1;
71
72 single_extra = 12;
73 svector = new double[single_extra];
74
75 neighprev = 0;
76
77 nmax = 0;
78 mass_rigid = nullptr;
79
80 onerad_dynamic = nullptr;
81 onerad_frozen = nullptr;
82 maxrad_dynamic = nullptr;
83 maxrad_frozen = nullptr;
84
85 limit_damping = nullptr;
86
87 history_transfer_factors = nullptr;
88
89 dt = update->dt;
90
91 // set comm size needed by this Pair if used with fix rigid
92
93 comm_forward = 1;
94
95 use_history = 0;
96 beyond_contact = 0;
97 nondefault_history_transfer = 0;
98 tangential_history_index = 0;
99 roll_history_index = twist_history_index = 0;
100
101 // create dummy fix as placeholder for FixNeighHistory
102 // this is so final order of Modify:fix will conform to input script
103
104 fix_history = nullptr;
105 fix_dummy = (FixDummy *) modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY");
106 }
107
108 /* ---------------------------------------------------------------------- */
109
~PairGranular()110 PairGranular::~PairGranular()
111 {
112 delete[] svector;
113 delete[] history_transfer_factors;
114
115 if (!fix_history) modify->delete_fix("NEIGH_HISTORY_GRANULAR_DUMMY");
116 else modify->delete_fix("NEIGH_HISTORY_GRANULAR");
117
118 if (allocated) {
119 memory->destroy(setflag);
120 memory->destroy(cutsq);
121 memory->destroy(cutoff_type);
122
123 memory->destroy(normal_coeffs);
124 memory->destroy(tangential_coeffs);
125 memory->destroy(roll_coeffs);
126 memory->destroy(twist_coeffs);
127
128 memory->destroy(Emod);
129 memory->destroy(poiss);
130
131 memory->destroy(normal_model);
132 memory->destroy(damping_model);
133 memory->destroy(tangential_model);
134 memory->destroy(roll_model);
135 memory->destroy(twist_model);
136 memory->destroy(limit_damping);
137
138 delete [] onerad_dynamic;
139 delete [] onerad_frozen;
140 delete [] maxrad_dynamic;
141 delete [] maxrad_frozen;
142 }
143
144 memory->destroy(mass_rigid);
145 }
146
147 /* ---------------------------------------------------------------------- */
148
compute(int eflag,int vflag)149 void PairGranular::compute(int eflag, int vflag)
150 {
151 int i,j,ii,jj,inum,jnum,itype,jtype;
152 double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz;
153 double radi,radj,radsum,rsq,r,rinv;
154 double Reff, delta, dR, dR2, dist_to_contact;
155
156 double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
157 double wr1,wr2,wr3;
158 double vtr1,vtr2,vtr3,vrel;
159
160 double knfac, damp_normal=0.0, damp_normal_prefactor;
161 double k_tangential, damp_tangential;
162 double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit;
163 double fs, fs1, fs2, fs3, tor1, tor2, tor3;
164
165 double mi,mj,meff;
166 double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;
167
168 // for JKR
169 double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E;
170 double t0, t1, t2, t3, t4, t5, t6;
171 double sqrt1, sqrt2, sqrt3;
172
173 // rolling
174 double k_roll, damp_roll;
175 double torroll1, torroll2, torroll3;
176 double rollmag, rolldotn, scalefac;
177 double fr, fr1, fr2, fr3;
178
179 // twisting
180 double k_twist, damp_twist, mu_twist;
181 double signtwist, magtwist, magtortwist, Mtcrit;
182 double tortwist1, tortwist2, tortwist3;
183
184 double shrmag,rsht,prjmag;
185 bool frameupdate;
186 int *ilist,*jlist,*numneigh,**firstneigh;
187 int *touch,**firsttouch;
188 double *history,*allhistory,**firsthistory;
189
190 bool touchflag = false;
191 const bool historyupdate = (update->setupflag) ? false : true;
192
193 ev_init(eflag,vflag);
194
195 // update rigid body info for owned & ghost atoms if using FixRigid masses
196 // body[i] = which body atom I is in, -1 if none
197 // mass_body = mass of each rigid body
198
199 if (fix_rigid && neighbor->ago == 0) {
200 int tmp;
201 int *body = (int *) fix_rigid->extract("body",tmp);
202 double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
203 if (atom->nmax > nmax) {
204 memory->destroy(mass_rigid);
205 nmax = atom->nmax;
206 memory->create(mass_rigid,nmax,"pair:mass_rigid");
207 }
208 int nlocal = atom->nlocal;
209 for (i = 0; i < nlocal; i++)
210 if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
211 else mass_rigid[i] = 0.0;
212 comm->forward_comm_pair(this);
213 }
214
215 double **x = atom->x;
216 double **v = atom->v;
217 double **f = atom->f;
218 int *type = atom->type;
219 double **omega = atom->omega;
220 double **torque = atom->torque;
221 double *radius = atom->radius;
222 double *rmass = atom->rmass;
223 int *mask = atom->mask;
224 int nlocal = atom->nlocal;
225
226 inum = list->inum;
227 ilist = list->ilist;
228 numneigh = list->numneigh;
229 firstneigh = list->firstneigh;
230 if (use_history) {
231 firsttouch = fix_history->firstflag;
232 firsthistory = fix_history->firstvalue;
233 }
234
235 for (ii = 0; ii < inum; ii++) {
236 i = ilist[ii];
237 itype = type[i];
238 xtmp = x[i][0];
239 ytmp = x[i][1];
240 ztmp = x[i][2];
241 itype = type[i];
242 radi = radius[i];
243 if (use_history) {
244 touch = firsttouch[i];
245 allhistory = firsthistory[i];
246 }
247 jlist = firstneigh[i];
248 jnum = numneigh[i];
249
250 for (jj = 0; jj < jnum; jj++) {
251 j = jlist[jj];
252 j &= NEIGHMASK;
253
254 delx = xtmp - x[j][0];
255 dely = ytmp - x[j][1];
256 delz = ztmp - x[j][2];
257 jtype = type[j];
258 rsq = delx*delx + dely*dely + delz*delz;
259 radj = radius[j];
260 radsum = radi + radj;
261
262 E = normal_coeffs[itype][jtype][0];
263 Reff = radi*radj/radsum;
264 touchflag = false;
265
266 if (normal_model[itype][jtype] == JKR) {
267 E *= THREEQUARTERS;
268 if (touch[jj]) {
269 R2 = Reff*Reff;
270 coh = normal_coeffs[itype][jtype][3];
271 a = cbrt(9.0*MY_PI*coh*R2/(4*E));
272 delta_pulloff = a*a/Reff - 2*sqrt(MY_PI*coh*a/E);
273 dist_pulloff = radsum-delta_pulloff;
274 touchflag = (rsq < dist_pulloff*dist_pulloff);
275 } else {
276 touchflag = (rsq < radsum*radsum);
277 }
278 } else {
279 touchflag = (rsq < radsum*radsum);
280 }
281
282 if (!touchflag) {
283 // unset non-touching neighbors
284 if (use_history) {
285 touch[jj] = 0;
286 history = &allhistory[size_history*jj];
287 for (int k = 0; k < size_history; k++) history[k] = 0.0;
288 }
289 } else {
290 r = sqrt(rsq);
291 rinv = 1.0/r;
292
293 nx = delx*rinv;
294 ny = dely*rinv;
295 nz = delz*rinv;
296
297 // relative translational velocity
298
299 vr1 = v[i][0] - v[j][0];
300 vr2 = v[i][1] - v[j][1];
301 vr3 = v[i][2] - v[j][2];
302
303 // normal component
304
305 vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n
306 vn1 = nx*vnnr;
307 vn2 = ny*vnnr;
308 vn3 = nz*vnnr;
309
310 // meff = effective mass of pair of particles
311 // if I or J part of rigid body, use body mass
312 // if I or J is frozen, meff is other particle
313
314 mi = rmass[i];
315 mj = rmass[j];
316 if (fix_rigid) {
317 if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
318 if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
319 }
320
321 meff = mi*mj / (mi+mj);
322 if (mask[i] & freeze_group_bit) meff = mj;
323 if (mask[j] & freeze_group_bit) meff = mi;
324
325 delta = radsum - r;
326 dR = delta*Reff;
327
328 if (normal_model[itype][jtype] == JKR) {
329 touch[jj] = 1;
330 R2=Reff*Reff;
331 coh = normal_coeffs[itype][jtype][3];
332 dR2 = dR*dR;
333 t0 = coh*coh*R2*R2*E;
334 t1 = PI27SQ*t0;
335 t2 = 8*dR*dR2*E*E*E;
336 t3 = 4*dR2*E;
337 // in case sqrt(0) < 0 due to precision issues
338 sqrt1 = MAX(0, t0*(t1+2*t2));
339 t4 = cbrt(t1+t2+THREEROOT3*MY_PI*sqrt(sqrt1));
340 t5 = t3/t4 + t4/E;
341 sqrt2 = MAX(0, 2*dR + t5);
342 t6 = sqrt(sqrt2);
343 sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*MY_PI*R2/(E*t6));
344 a = INVROOT6*(t6 + sqrt(sqrt3));
345 a2 = a*a;
346 knfac = normal_coeffs[itype][jtype][0]*a;
347 Fne = knfac*a2/Reff - MY_2PI*a2*sqrt(4*coh*E/(MY_PI*a));
348 } else {
349 knfac = E; // Hooke
350 Fne = knfac*delta;
351 a = sqrt(dR);
352 if (normal_model[itype][jtype] != HOOKE) {
353 Fne *= a;
354 knfac *= a;
355 }
356 if (normal_model[itype][jtype] == DMT)
357 Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff;
358 }
359
360 // NOTE: consider restricting Hooke to only have
361 // 'velocity' as an option for damping?
362
363 if (damping_model[itype][jtype] == VELOCITY) {
364 damp_normal = 1;
365 } else if (damping_model[itype][jtype] == MASS_VELOCITY) {
366 damp_normal = meff;
367 } else if (damping_model[itype][jtype] == VISCOELASTIC) {
368 damp_normal = a*meff;
369 } else if (damping_model[itype][jtype] == TSUJI) {
370 damp_normal = sqrt(meff*knfac);
371 } else damp_normal = 0.0;
372
373 damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal;
374 Fdamp = -damp_normal_prefactor*vnnr;
375
376 Fntot = Fne + Fdamp;
377 if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0;
378
379 //****************************************
380 // tangential force, including history effects
381 //****************************************
382
383 // For linear, mindlin, mindlin_rescale:
384 // history = cumulative tangential displacement
385 //
386 // For mindlin/force, mindlin_rescale/force:
387 // history = cumulative tangential elastic force
388
389 // tangential component
390 vt1 = vr1 - vn1;
391 vt2 = vr2 - vn2;
392 vt3 = vr3 - vn3;
393
394 // relative rotational velocity
395 wr1 = (radi*omega[i][0] + radj*omega[j][0]);
396 wr2 = (radi*omega[i][1] + radj*omega[j][1]);
397 wr3 = (radi*omega[i][2] + radj*omega[j][2]);
398
399 // relative tangential velocities
400 vtr1 = vt1 - (nz*wr2-ny*wr3);
401 vtr2 = vt2 - (nx*wr3-nz*wr1);
402 vtr3 = vt3 - (ny*wr1-nx*wr2);
403 vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
404 vrel = sqrt(vrel);
405
406 // if any history is needed
407 if (use_history) {
408 touch[jj] = 1;
409 history = &allhistory[size_history*jj];
410 }
411
412 if (normal_model[itype][jtype] == JKR) {
413 F_pulloff = 3*MY_PI*coh*Reff;
414 Fncrit = fabs(Fne + 2*F_pulloff);
415 } else if (normal_model[itype][jtype] == DMT) {
416 F_pulloff = 4*MY_PI*coh*Reff;
417 Fncrit = fabs(Fne + 2*F_pulloff);
418 } else {
419 Fncrit = fabs(Fntot);
420 }
421 Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;
422
423 //------------------------------
424 // tangential forces
425 //------------------------------
426 k_tangential = tangential_coeffs[itype][jtype][0];
427 damp_tangential = tangential_coeffs[itype][jtype][1] *
428 damp_normal_prefactor;
429
430 if (tangential_history) {
431 if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
432 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE) {
433 k_tangential *= a;
434 } else if (tangential_model[itype][jtype] ==
435 TANGENTIAL_MINDLIN_RESCALE ||
436 tangential_model[itype][jtype] ==
437 TANGENTIAL_MINDLIN_RESCALE_FORCE) {
438 k_tangential *= a;
439 // on unloading, rescale the shear displacements/force
440 if (a < history[3]) {
441 double factor = a/history[3];
442 history[0] *= factor;
443 history[1] *= factor;
444 history[2] *= factor;
445 }
446 }
447 // rotate and update displacements / force.
448 // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
449 if (historyupdate) {
450 rsht = history[0]*nx + history[1]*ny + history[2]*nz;
451 if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE ||
452 tangential_model[itype][jtype] ==
453 TANGENTIAL_MINDLIN_RESCALE_FORCE)
454 frameupdate = fabs(rsht) > EPSILON*Fscrit;
455 else
456 frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit;
457 if (frameupdate) {
458 shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
459 history[2]*history[2]);
460 // projection
461 history[0] -= rsht*nx;
462 history[1] -= rsht*ny;
463 history[2] -= rsht*nz;
464
465 // also rescale to preserve magnitude
466 prjmag = sqrt(history[0]*history[0] + history[1]*history[1] +
467 history[2]*history[2]);
468 if (prjmag > 0) scalefac = shrmag/prjmag;
469 else scalefac = 0;
470 history[0] *= scalefac;
471 history[1] *= scalefac;
472 history[2] *= scalefac;
473 }
474 // update history
475 if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
476 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
477 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
478 // tangential displacement
479 history[0] += vtr1*dt;
480 history[1] += vtr2*dt;
481 history[2] += vtr3*dt;
482 } else {
483 // tangential force
484 // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
485 history[0] -= k_tangential*vtr1*dt;
486 history[1] -= k_tangential*vtr2*dt;
487 history[2] -= k_tangential*vtr3*dt;
488 }
489 if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE ||
490 tangential_model[itype][jtype] ==
491 TANGENTIAL_MINDLIN_RESCALE_FORCE)
492 history[3] = a;
493 }
494
495 // tangential forces = history + tangential velocity damping
496 if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
497 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
498 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
499 fs1 = -k_tangential*history[0] - damp_tangential*vtr1;
500 fs2 = -k_tangential*history[1] - damp_tangential*vtr2;
501 fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
502 } else {
503 fs1 = history[0] - damp_tangential*vtr1;
504 fs2 = history[1] - damp_tangential*vtr2;
505 fs3 = history[2] - damp_tangential*vtr3;
506 }
507
508 // rescale frictional displacements and forces if needed
509 fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
510 if (fs > Fscrit) {
511 shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
512 history[2]*history[2]);
513 if (shrmag != 0.0) {
514 if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
515 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
516 tangential_model[itype][jtype] ==
517 TANGENTIAL_MINDLIN_RESCALE) {
518 history[0] = -1.0/k_tangential*(Fscrit*fs1/fs +
519 damp_tangential*vtr1);
520 history[1] = -1.0/k_tangential*(Fscrit*fs2/fs +
521 damp_tangential*vtr2);
522 history[2] = -1.0/k_tangential*(Fscrit*fs3/fs +
523 damp_tangential*vtr3);
524 } else {
525 history[0] = Fscrit*fs1/fs + damp_tangential*vtr1;
526 history[1] = Fscrit*fs2/fs + damp_tangential*vtr2;
527 history[2] = Fscrit*fs3/fs + damp_tangential*vtr3;
528 }
529 fs1 *= Fscrit/fs;
530 fs2 *= Fscrit/fs;
531 fs3 *= Fscrit/fs;
532 } else fs1 = fs2 = fs3 = 0.0;
533 }
534 } else { // classic pair gran/hooke (no history)
535 fs = damp_tangential*vrel;
536 if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel;
537 else Ft = 0.0;
538 fs1 = -Ft*vtr1;
539 fs2 = -Ft*vtr2;
540 fs3 = -Ft*vtr3;
541 }
542
543 if (roll_model[itype][jtype] != ROLL_NONE ||
544 twist_model[itype][jtype] != TWIST_NONE) {
545 relrot1 = omega[i][0] - omega[j][0];
546 relrot2 = omega[i][1] - omega[j][1];
547 relrot3 = omega[i][2] - omega[j][2];
548 }
549 //****************************************
550 // rolling resistance
551 //****************************************
552
553 if (roll_model[itype][jtype] != ROLL_NONE) {
554 // rolling velocity,
555 // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
556 // this is different from the Marshall papers,
557 // which use the Bagi/Kuhn formulation
558 // for rolling velocity (see Wang et al for why the latter is wrong)
559 vrl1 = Reff*(relrot2*nz - relrot3*ny);
560 vrl2 = Reff*(relrot3*nx - relrot1*nz);
561 vrl3 = Reff*(relrot1*ny - relrot2*nx);
562
563 int rhist0 = roll_history_index;
564 int rhist1 = rhist0 + 1;
565 int rhist2 = rhist1 + 1;
566
567 k_roll = roll_coeffs[itype][jtype][0];
568 damp_roll = roll_coeffs[itype][jtype][1];
569 Frcrit = roll_coeffs[itype][jtype][2] * Fncrit;
570
571 if (historyupdate) {
572 rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
573 frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit;
574 if (frameupdate) { // rotate into tangential plane
575 rollmag = sqrt(history[rhist0]*history[rhist0] +
576 history[rhist1]*history[rhist1] +
577 history[rhist2]*history[rhist2]);
578 // projection
579 history[rhist0] -= rolldotn*nx;
580 history[rhist1] -= rolldotn*ny;
581 history[rhist2] -= rolldotn*nz;
582 // also rescale to preserve magnitude
583 prjmag = sqrt(history[rhist0]*history[rhist0] +
584 history[rhist1]*history[rhist1] +
585 history[rhist2]*history[rhist2]);
586 if (prjmag > 0) scalefac = rollmag/prjmag;
587 else scalefac = 0;
588 history[rhist0] *= scalefac;
589 history[rhist1] *= scalefac;
590 history[rhist2] *= scalefac;
591 }
592 history[rhist0] += vrl1*dt;
593 history[rhist1] += vrl2*dt;
594 history[rhist2] += vrl3*dt;
595 }
596
597 fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
598 fr2 = -k_roll*history[rhist1] - damp_roll*vrl2;
599 fr3 = -k_roll*history[rhist2] - damp_roll*vrl3;
600
601 // rescale frictional displacements and forces if needed
602
603 fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
604 if (fr > Frcrit) {
605 rollmag = sqrt(history[rhist0]*history[rhist0] +
606 history[rhist1]*history[rhist1] +
607 history[rhist2]*history[rhist2]);
608 if (rollmag != 0.0) {
609 history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1);
610 history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2);
611 history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3);
612 fr1 *= Frcrit/fr;
613 fr2 *= Frcrit/fr;
614 fr3 *= Frcrit/fr;
615 } else fr1 = fr2 = fr3 = 0.0;
616 }
617 }
618
619 //****************************************
620 // twisting torque, including history effects
621 //****************************************
622
623 if (twist_model[itype][jtype] != TWIST_NONE) {
624 // omega_T (eq 29 of Marshall)
625 magtwist = relrot1*nx + relrot2*ny + relrot3*nz;
626 if (twist_model[itype][jtype] == TWIST_MARSHALL) {
627 k_twist = 0.5*k_tangential*a*a;; // eq 32 of Marshall paper
628 damp_twist = 0.5*damp_tangential*a*a;
629 mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2];
630 } else {
631 k_twist = twist_coeffs[itype][jtype][0];
632 damp_twist = twist_coeffs[itype][jtype][1];
633 mu_twist = twist_coeffs[itype][jtype][2];
634 }
635 if (historyupdate) {
636 history[twist_history_index] += magtwist*dt;
637 }
638 magtortwist = -k_twist*history[twist_history_index] -
639 damp_twist*magtwist; // M_t torque (eq 30)
640 signtwist = (magtwist > 0) - (magtwist < 0);
641 Mtcrit = mu_twist*Fncrit; // critical torque (eq 44)
642 if (fabs(magtortwist) > Mtcrit) {
643 history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist -
644 damp_twist*magtwist);
645 magtortwist = -Mtcrit * signtwist; // eq 34
646 }
647 }
648
649 // apply forces & torques
650
651 fx = nx*Fntot + fs1;
652 fy = ny*Fntot + fs2;
653 fz = nz*Fntot + fs3;
654
655 f[i][0] += fx;
656 f[i][1] += fy;
657 f[i][2] += fz;
658
659 tor1 = ny*fs3 - nz*fs2;
660 tor2 = nz*fs1 - nx*fs3;
661 tor3 = nx*fs2 - ny*fs1;
662
663 dist_to_contact = radi-0.5*delta;
664 torque[i][0] -= dist_to_contact*tor1;
665 torque[i][1] -= dist_to_contact*tor2;
666 torque[i][2] -= dist_to_contact*tor3;
667
668 if (twist_model[itype][jtype] != TWIST_NONE) {
669 tortwist1 = magtortwist * nx;
670 tortwist2 = magtortwist * ny;
671 tortwist3 = magtortwist * nz;
672
673 torque[i][0] += tortwist1;
674 torque[i][1] += tortwist2;
675 torque[i][2] += tortwist3;
676 }
677
678 if (roll_model[itype][jtype] != ROLL_NONE) {
679 torroll1 = Reff*(ny*fr3 - nz*fr2); // n cross fr
680 torroll2 = Reff*(nz*fr1 - nx*fr3);
681 torroll3 = Reff*(nx*fr2 - ny*fr1);
682
683 torque[i][0] += torroll1;
684 torque[i][1] += torroll2;
685 torque[i][2] += torroll3;
686 }
687
688 if (force->newton_pair || j < nlocal) {
689 f[j][0] -= fx;
690 f[j][1] -= fy;
691 f[j][2] -= fz;
692
693 dist_to_contact = radj-0.5*delta;
694 torque[j][0] -= dist_to_contact*tor1;
695 torque[j][1] -= dist_to_contact*tor2;
696 torque[j][2] -= dist_to_contact*tor3;
697
698 if (twist_model[itype][jtype] != TWIST_NONE) {
699 torque[j][0] -= tortwist1;
700 torque[j][1] -= tortwist2;
701 torque[j][2] -= tortwist3;
702 }
703 if (roll_model[itype][jtype] != ROLL_NONE) {
704 torque[j][0] -= torroll1;
705 torque[j][1] -= torroll2;
706 torque[j][2] -= torroll3;
707 }
708 }
709 if (evflag) ev_tally_xyz(i,j,nlocal,force->newton_pair,
710 0.0,0.0,fx,fy,fz,delx,dely,delz);
711 }
712 }
713 }
714 }
715
716 /* ----------------------------------------------------------------------
717 allocate all arrays
718 ------------------------------------------------------------------------- */
719
allocate()720 void PairGranular::allocate()
721 {
722 allocated = 1;
723 int n = atom->ntypes;
724
725 memory->create(setflag,n+1,n+1,"pair:setflag");
726 for (int i = 1; i <= n; i++)
727 for (int j = i; j <= n; j++)
728 setflag[i][j] = 0;
729
730 memory->create(cutsq,n+1,n+1,"pair:cutsq");
731 memory->create(cutoff_type,n+1,n+1,"pair:cutoff_type");
732 memory->create(normal_coeffs,n+1,n+1,4,"pair:normal_coeffs");
733 memory->create(tangential_coeffs,n+1,n+1,3,"pair:tangential_coeffs");
734 memory->create(roll_coeffs,n+1,n+1,3,"pair:roll_coeffs");
735 memory->create(twist_coeffs,n+1,n+1,3,"pair:twist_coeffs");
736
737 memory->create(Emod,n+1,n+1,"pair:Emod");
738 memory->create(poiss,n+1,n+1,"pair:poiss");
739
740 memory->create(normal_model,n+1,n+1,"pair:normal_model");
741 memory->create(damping_model,n+1,n+1,"pair:damping_model");
742 memory->create(tangential_model,n+1,n+1,"pair:tangential_model");
743 memory->create(roll_model,n+1,n+1,"pair:roll_model");
744 memory->create(twist_model,n+1,n+1,"pair:twist_model");
745 memory->create(limit_damping,n+1,n+1,"pair:limit_damping");
746
747 onerad_dynamic = new double[n+1];
748 onerad_frozen = new double[n+1];
749 maxrad_dynamic = new double[n+1];
750 maxrad_frozen = new double[n+1];
751 }
752
753 /* ----------------------------------------------------------------------
754 global settings
755 ------------------------------------------------------------------------- */
756
settings(int narg,char ** arg)757 void PairGranular::settings(int narg, char **arg)
758 {
759 if (narg == 1) {
760 cutoff_global = utils::numeric(FLERR,arg[0],false,lmp);
761 } else {
762 cutoff_global = -1; // will be set based on particle sizes, model choice
763 }
764
765 normal_history = tangential_history = 0;
766 roll_history = twist_history = 0;
767 }
768
769 /* ----------------------------------------------------------------------
770 set coeffs for one or more type pairs
771 ------------------------------------------------------------------------- */
772
coeff(int narg,char ** arg)773 void PairGranular::coeff(int narg, char **arg)
774 {
775 int normal_model_one, damping_model_one;
776 int tangential_model_one, roll_model_one, twist_model_one;
777
778 double normal_coeffs_one[4];
779 double tangential_coeffs_one[4];
780 double roll_coeffs_one[4];
781 double twist_coeffs_one[4];
782
783 double cutoff_one = -1;
784
785 if (narg < 2)
786 error->all(FLERR,"Incorrect args for pair coefficients");
787
788 if (!allocated) allocate();
789
790 int ilo,ihi,jlo,jhi;
791 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
792 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
793
794 //Defaults
795 normal_model_one = tangential_model_one = -1;
796 roll_model_one = ROLL_NONE;
797 twist_model_one = TWIST_NONE;
798 damping_model_one = VISCOELASTIC;
799 int ld_flag = 0;
800
801 int iarg = 2;
802 while (iarg < narg) {
803 if (strcmp(arg[iarg], "hooke") == 0) {
804 if (iarg + 2 >= narg)
805 error->all(FLERR,"Illegal pair_coeff command, "
806 "not enough parameters provided for Hooke option");
807 normal_model_one = HOOKE;
808 normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // kn
809 normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping
810 iarg += 3;
811 } else if (strcmp(arg[iarg], "hertz") == 0) {
812 if (iarg + 2 >= narg)
813 error->all(FLERR,"Illegal pair_coeff command, "
814 "not enough parameters provided for Hertz option");
815 normal_model_one = HERTZ;
816 normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // kn
817 normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping
818 iarg += 3;
819 } else if (strcmp(arg[iarg], "hertz/material") == 0) {
820 if (iarg + 3 >= narg)
821 error->all(FLERR,"Illegal pair_coeff command, "
822 "not enough parameters provided for Hertz/material option");
823 normal_model_one = HERTZ_MATERIAL;
824 normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E
825 normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping
826 normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio
827 iarg += 4;
828 } else if (strcmp(arg[iarg], "dmt") == 0) {
829 if (iarg + 4 >= narg)
830 error->all(FLERR,"Illegal pair_coeff command, "
831 "not enough parameters provided for Hertz option");
832 normal_model_one = DMT;
833 normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E
834 normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping
835 normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio
836 normal_coeffs_one[3] = utils::numeric(FLERR,arg[iarg+4],false,lmp); // cohesion
837 iarg += 5;
838 } else if (strcmp(arg[iarg], "jkr") == 0) {
839 if (iarg + 4 >= narg)
840 error->all(FLERR,"Illegal pair_coeff command, "
841 "not enough parameters provided for JKR option");
842 beyond_contact = 1;
843 normal_model_one = JKR;
844 normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E
845 normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping
846 normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio
847 normal_coeffs_one[3] = utils::numeric(FLERR,arg[iarg+4],false,lmp); // cohesion
848 iarg += 5;
849 } else if (strcmp(arg[iarg], "damping") == 0) {
850 if (iarg+1 >= narg)
851 error->all(FLERR, "Illegal pair_coeff command, "
852 "not enough parameters provided for damping model");
853 if (strcmp(arg[iarg+1], "velocity") == 0) {
854 damping_model_one = VELOCITY;
855 iarg += 1;
856 } else if (strcmp(arg[iarg+1], "mass_velocity") == 0) {
857 damping_model_one = MASS_VELOCITY;
858 iarg += 1;
859 } else if (strcmp(arg[iarg+1], "viscoelastic") == 0) {
860 damping_model_one = VISCOELASTIC;
861 iarg += 1;
862 } else if (strcmp(arg[iarg+1], "tsuji") == 0) {
863 damping_model_one = TSUJI;
864 iarg += 1;
865 } else error->all(FLERR, "Illegal pair_coeff command, "
866 "unrecognized damping model");
867 iarg += 1;
868 } else if (strcmp(arg[iarg], "tangential") == 0) {
869 if (iarg + 1 >= narg)
870 error->all(FLERR,"Illegal pair_coeff command, must specify "
871 "tangential model after tangential keyword");
872 if (strcmp(arg[iarg+1], "linear_nohistory") == 0) {
873 if (iarg + 3 >= narg)
874 error->all(FLERR,"Illegal pair_coeff command, "
875 "not enough parameters provided for tangential model");
876 tangential_model_one = TANGENTIAL_NOHISTORY;
877 tangential_coeffs_one[0] = 0;
878 // gammat and friction coeff
879 tangential_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
880 tangential_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
881 iarg += 4;
882 } else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||
883 (strcmp(arg[iarg+1], "mindlin") == 0) ||
884 (strcmp(arg[iarg+1], "mindlin_rescale") == 0) ||
885 (strcmp(arg[iarg+1], "mindlin/force") == 0) ||
886 (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) {
887 if (iarg + 4 >= narg)
888 error->all(FLERR,"Illegal pair_coeff command, "
889 "not enough parameters provided for tangential model");
890 if (strcmp(arg[iarg+1], "linear_history") == 0)
891 tangential_model_one = TANGENTIAL_HISTORY;
892 else if (strcmp(arg[iarg+1], "mindlin") == 0)
893 tangential_model_one = TANGENTIAL_MINDLIN;
894 else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
895 tangential_model_one = TANGENTIAL_MINDLIN_RESCALE;
896 else if (strcmp(arg[iarg+1], "mindlin/force") == 0)
897 tangential_model_one = TANGENTIAL_MINDLIN_FORCE;
898 else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)
899 tangential_model_one = TANGENTIAL_MINDLIN_RESCALE_FORCE;
900 tangential_history = 1;
901 if ((tangential_model_one == TANGENTIAL_MINDLIN ||
902 tangential_model_one == TANGENTIAL_MINDLIN_RESCALE ||
903 tangential_model_one == TANGENTIAL_MINDLIN_FORCE ||
904 tangential_model_one == TANGENTIAL_MINDLIN_RESCALE_FORCE) &&
905 (strcmp(arg[iarg+2], "NULL") == 0)) {
906 if (normal_model_one == HERTZ || normal_model_one == HOOKE) {
907 error->all(FLERR, "NULL setting for Mindlin tangential "
908 "stiffness requires a normal contact model that "
909 "specifies material properties");
910 }
911 tangential_coeffs_one[0] = -1;
912 } else {
913 tangential_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // kt
914 }
915 // gammat and friction coeff
916 tangential_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
917 tangential_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp);
918 iarg += 5;
919 } else {
920 error->all(FLERR, "Illegal pair_coeff command, "
921 "tangential model not recognized");
922 }
923 } else if (strcmp(arg[iarg], "rolling") == 0) {
924 if (iarg + 1 >= narg)
925 error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
926 if (strcmp(arg[iarg+1], "none") == 0) {
927 roll_model_one = ROLL_NONE;
928 iarg += 2;
929 } else if (strcmp(arg[iarg+1], "sds") == 0) {
930 if (iarg + 4 >= narg)
931 error->all(FLERR,"Illegal pair_coeff command, "
932 "not enough parameters provided for rolling model");
933 roll_model_one = ROLL_SDS;
934 roll_history = 1;
935 // kR and gammaR and rolling friction coeff
936 roll_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
937 roll_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
938 roll_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp);
939 iarg += 5;
940 } else {
941 error->all(FLERR, "Illegal pair_coeff command, "
942 "rolling friction model not recognized");
943 }
944 } else if (strcmp(arg[iarg], "twisting") == 0) {
945 if (iarg + 1 >= narg)
946 error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
947 if (strcmp(arg[iarg+1], "none") == 0) {
948 twist_model_one = TWIST_NONE;
949 iarg += 2;
950 } else if (strcmp(arg[iarg+1], "marshall") == 0) {
951 twist_model_one = TWIST_MARSHALL;
952 twist_history = 1;
953 iarg += 2;
954 } else if (strcmp(arg[iarg+1], "sds") == 0) {
955 if (iarg + 4 >= narg)
956 error->all(FLERR,"Illegal pair_coeff command, "
957 "not enough parameters provided for twist model");
958 twist_model_one = TWIST_SDS;
959 twist_history = 1;
960 // kt and gammat and friction coeff
961 twist_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
962 twist_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
963 twist_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp);
964 iarg += 5;
965 } else {
966 error->all(FLERR, "Illegal pair_coeff command, "
967 "twisting friction model not recognized");
968 }
969 } else if (strcmp(arg[iarg], "cutoff") == 0) {
970 if (iarg + 1 >= narg)
971 error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
972 cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp);
973 iarg += 2;
974 } else if (strcmp(arg[iarg], "limit_damping") == 0) {
975 ld_flag = 1;
976 iarg += 1;
977 } else error->all(FLERR, "Illegal pair_coeff command");
978 }
979
980 // error not to specify normal or tangential model
981 if ((normal_model_one < 0) || (tangential_model_one < 0))
982 error->all(FLERR, "Illegal pair_coeff command, "
983 "must specify normal or tangential contact model");
984
985 int count = 0;
986 double damp;
987 if (damping_model_one == TSUJI) {
988 double cor;
989 cor = normal_coeffs_one[1];
990 damp = 1.2728-4.2783*cor+11.087*square(cor)-22.348*cube(cor)+
991 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6);
992 } else damp = normal_coeffs_one[1];
993
994 if (ld_flag && normal_model_one == JKR)
995 error->all(FLERR,"Illegal pair_coeff command, "
996 "Cannot limit damping with JKR model");
997
998 if (ld_flag && normal_model_one == DMT)
999 error->all(FLERR,"Illegal pair_coeff command, "
1000 "Cannot limit damping with DMT model");
1001
1002 for (int i = ilo; i <= ihi; i++) {
1003 for (int j = MAX(jlo,i); j <= jhi; j++) {
1004 normal_model[i][j] = normal_model[j][i] = normal_model_one;
1005 normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = damp;
1006 if (normal_model_one != HERTZ && normal_model_one != HOOKE) {
1007 Emod[i][j] = Emod[j][i] = normal_coeffs_one[0];
1008 poiss[i][j] = poiss[j][i] = normal_coeffs_one[2];
1009 normal_coeffs[i][j][0] = normal_coeffs[j][i][0] =
1010 FOURTHIRDS*mix_stiffnessE(Emod[i][j],Emod[i][j],
1011 poiss[i][j],poiss[i][j]);
1012 } else {
1013 normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0];
1014 }
1015 if ((normal_model_one == JKR) || (normal_model_one == DMT))
1016 normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = normal_coeffs_one[3];
1017
1018 damping_model[i][j] = damping_model[j][i] = damping_model_one;
1019
1020 tangential_model[i][j] = tangential_model[j][i] = tangential_model_one;
1021 if (tangential_coeffs_one[0] == -1) {
1022 tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] =
1023 8*mix_stiffnessG(Emod[i][j],Emod[i][j],poiss[i][j],poiss[i][j]);
1024 } else {
1025 tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] =
1026 tangential_coeffs_one[0];
1027 }
1028 for (int k = 1; k < 3; k++)
1029 tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] =
1030 tangential_coeffs_one[k];
1031
1032 roll_model[i][j] = roll_model[j][i] = roll_model_one;
1033 if (roll_model_one != ROLL_NONE)
1034 for (int k = 0; k < 3; k++)
1035 roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = roll_coeffs_one[k];
1036
1037 twist_model[i][j] = twist_model[j][i] = twist_model_one;
1038 if (twist_model_one != TWIST_NONE && twist_model_one != TWIST_MARSHALL)
1039 for (int k = 0; k < 3; k++)
1040 twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_one[k];
1041
1042 cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one;
1043
1044 limit_damping[i][j] = limit_damping[j][i] = ld_flag;
1045
1046 setflag[i][j] = 1;
1047 count++;
1048 }
1049 }
1050
1051
1052 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
1053 }
1054
1055 /* ----------------------------------------------------------------------
1056 init specific to this pair style
1057 ------------------------------------------------------------------------- */
1058
init_style()1059 void PairGranular::init_style()
1060 {
1061 int i;
1062
1063 // error and warning checks
1064
1065 if (!atom->radius_flag || !atom->rmass_flag)
1066 error->all(FLERR,"Pair granular requires atom attributes radius, rmass");
1067 if (comm->ghost_velocity == 0)
1068 error->all(FLERR,"Pair granular requires ghost atoms store velocity");
1069
1070 // determine whether we need a granular neigh list, how large it needs to be
1071
1072 use_history = normal_history || tangential_history ||
1073 roll_history || twist_history;
1074
1075 // for JKR, will need fix/neigh/history to keep track of touch arrays
1076
1077 for (int i = 1; i <= atom->ntypes; i++)
1078 for (int j = i; j <= atom->ntypes; j++)
1079 if (normal_model[i][j] == JKR) use_history = 1;
1080
1081 size_history = 3*tangential_history + 3*roll_history + twist_history;
1082
1083 // determine location of tangential/roll/twist histories in array
1084
1085 if (roll_history) {
1086 if (tangential_history) roll_history_index = 3;
1087 else roll_history_index = 0;
1088 }
1089 if (twist_history) {
1090 if (tangential_history) {
1091 if (roll_history) twist_history_index = 6;
1092 else twist_history_index = 3;
1093 } else {
1094 if (roll_history) twist_history_index = 3;
1095 else twist_history_index = 0;
1096 }
1097 }
1098 for (int i = 1; i <= atom->ntypes; i++)
1099 for (int j = i; j <= atom->ntypes; j++)
1100 if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE ||
1101 tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE_FORCE) {
1102 size_history += 1;
1103 roll_history_index += 1;
1104 twist_history_index += 1;
1105 nondefault_history_transfer = 1;
1106 history_transfer_factors = new int[size_history];
1107 for (int ii = 0; ii < size_history; ++ii)
1108 history_transfer_factors[ii] = -1;
1109 history_transfer_factors[3] = 1;
1110 break;
1111 }
1112
1113 int irequest = neighbor->request(this,instance_me);
1114 neighbor->requests[irequest]->size = 1;
1115 if (use_history) neighbor->requests[irequest]->history = 1;
1116
1117 dt = update->dt;
1118
1119 // if history is stored and first init, create Fix to store history
1120 // it replaces FixDummy, created in the constructor
1121 // this is so its order in the fix list is preserved
1122
1123 if (use_history && fix_history == nullptr) {
1124 fix_history = (FixNeighHistory *) modify->replace_fix("NEIGH_HISTORY_GRANULAR_DUMMY",
1125 "NEIGH_HISTORY_GRANULAR"
1126 " all NEIGH_HISTORY "
1127 + std::to_string(size_history),1);
1128 fix_history->pair = this;
1129 }
1130
1131 // check for FixFreeze and set freeze_group_bit
1132
1133 int ifix = modify->find_fix_by_style("^freeze");
1134 if (ifix < 0) freeze_group_bit = 0;
1135 else freeze_group_bit = modify->fix[ifix]->groupbit;
1136
1137 // check for FixRigid so can extract rigid body masses
1138 // FIXME: this only catches the first rigid fix, there may be multiple.
1139
1140 fix_rigid = nullptr;
1141 for (i = 0; i < modify->nfix; i++)
1142 if (modify->fix[i]->rigid_flag) break;
1143 if (i < modify->nfix) fix_rigid = modify->fix[i];
1144
1145 // check for FixPour and FixDeposit so can extract particle radii
1146
1147 int ipour = modify->find_fix_by_style("^pour");
1148 int idep = modify->find_fix_by_style("^deposit");
1149
1150 // set maxrad_dynamic and maxrad_frozen for each type
1151 // include future FixPour and FixDeposit particles as dynamic
1152
1153 int itype;
1154 for (i = 1; i <= atom->ntypes; i++) {
1155 onerad_dynamic[i] = onerad_frozen[i] = 0.0;
1156 if (ipour >= 0) {
1157 itype = i;
1158 double radmax = *((double *) modify->fix[ipour]->extract("radius",itype));
1159 onerad_dynamic[i] = radmax;
1160 }
1161 if (idep >= 0) {
1162 itype = i;
1163 double radmax = *((double *) modify->fix[idep]->extract("radius",itype));
1164 onerad_dynamic[i] = radmax;
1165 }
1166 }
1167
1168 double *radius = atom->radius;
1169 int *mask = atom->mask;
1170 int *type = atom->type;
1171 int nlocal = atom->nlocal;
1172
1173 for (i = 0; i < nlocal; i++) {
1174 double radius_cut = radius[i];
1175 if (mask[i] & freeze_group_bit) {
1176 onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius_cut);
1177 } else {
1178 onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius_cut);
1179 }
1180 }
1181
1182 MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
1183 MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
1184
1185 // set fix which stores history info
1186
1187 if (size_history > 0) {
1188 int ifix = modify->find_fix("NEIGH_HISTORY_GRANULAR");
1189 if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID");
1190 fix_history = (FixNeighHistory *) modify->fix[ifix];
1191 }
1192 }
1193
1194 /* ----------------------------------------------------------------------
1195 init for one type pair i,j and corresponding j,i
1196 ------------------------------------------------------------------------- */
1197
init_one(int i,int j)1198 double PairGranular::init_one(int i, int j)
1199 {
1200 double cutoff=0.0;
1201
1202 if (setflag[i][j] == 0) {
1203 if ((normal_model[i][i] != normal_model[j][j]) ||
1204 (damping_model[i][i] != damping_model[j][j]) ||
1205 (tangential_model[i][i] != tangential_model[j][j]) ||
1206 (roll_model[i][i] != roll_model[j][j]) ||
1207 (twist_model[i][i] != twist_model[j][j])) {
1208 error->all(FLERR,"Granular pair style functional forms are different, "
1209 "cannot mix coefficients for types {} and {}. \n"
1210 "This combination must be set explicitly via a "
1211 "pair_coeff command",i,j);
1212 }
1213
1214 if (normal_model[i][j] == HERTZ || normal_model[i][j] == HOOKE)
1215 normal_coeffs[i][j][0] = normal_coeffs[j][i][0] =
1216 mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]);
1217 else
1218 normal_coeffs[i][j][0] = normal_coeffs[j][i][0] =
1219 mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]);
1220
1221 normal_coeffs[i][j][1] = normal_coeffs[j][i][1] =
1222 mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]);
1223 if ((normal_model[i][j] == JKR) || (normal_model[i][j] == DMT))
1224 normal_coeffs[i][j][3] = normal_coeffs[j][i][3] =
1225 mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]);
1226
1227 for (int k = 0; k < 3; k++)
1228 tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] =
1229 mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]);
1230
1231 if (roll_model[i][j] != ROLL_NONE) {
1232 for (int k = 0; k < 3; k++)
1233 roll_coeffs[i][j][k] = roll_coeffs[j][i][k] =
1234 mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]);
1235 }
1236
1237 if (twist_model[i][j] != TWIST_NONE && twist_model[i][j] != TWIST_MARSHALL) {
1238 for (int k = 0; k < 3; k++)
1239 twist_coeffs[i][j][k] = twist_coeffs[j][i][k] =
1240 mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]);
1241 }
1242 }
1243
1244 // It is possible that cut[i][j] at this point is still 0.0.
1245 // This can happen when
1246 // there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates
1247 // problems because neighbor.cpp uses min(cut[i][j]) to decide on the bin size
1248 // To avoid this issue, for cases involving cut[i][j] = 0.0 (possible only
1249 // if there is no current information about radius/cutoff of type i and j).
1250 // we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0.
1251
1252 double pulloff;
1253
1254 if (cutoff_type[i][j] < 0 && cutoff_global < 0) {
1255 if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) ||
1256 ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) ||
1257 // radius info about both i and j exist
1258 ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) {
1259 cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
1260 pulloff = 0.0;
1261 if (normal_model[i][j] == JKR) {
1262 pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j);
1263 cutoff += pulloff;
1264 }
1265
1266 if (normal_model[i][j] == JKR)
1267 pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j);
1268 cutoff = MAX(cutoff, maxrad_frozen[i]+maxrad_dynamic[j]+pulloff);
1269
1270 if (normal_model[i][j] == JKR)
1271 pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j], i, j);
1272 cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]+pulloff);
1273
1274 } else {
1275
1276 // radius info about either i or j does not exist
1277 // (i.e. not present and not about to get poured;
1278 // set to largest value to not interfere with neighbor list)
1279
1280 double cutmax = 0.0;
1281 for (int k = 1; k <= atom->ntypes; k++) {
1282 cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]);
1283 cutmax = MAX(cutmax,2.0*maxrad_frozen[k]);
1284 }
1285 cutoff = cutmax;
1286 }
1287 } else if (cutoff_type[i][j] > 0) {
1288 cutoff = cutoff_type[i][j];
1289 } else if (cutoff_global > 0) {
1290 cutoff = cutoff_global;
1291 }
1292
1293 return cutoff;
1294 }
1295
1296 /* ----------------------------------------------------------------------
1297 proc 0 writes to restart file
1298 ------------------------------------------------------------------------- */
1299
write_restart(FILE * fp)1300 void PairGranular::write_restart(FILE *fp)
1301 {
1302 int i,j;
1303 for (i = 1; i <= atom->ntypes; i++) {
1304 for (j = i; j <= atom->ntypes; j++) {
1305 fwrite(&setflag[i][j],sizeof(int),1,fp);
1306 if (setflag[i][j]) {
1307 fwrite(&normal_model[i][j],sizeof(int),1,fp);
1308 fwrite(&damping_model[i][j],sizeof(int),1,fp);
1309 fwrite(&tangential_model[i][j],sizeof(int),1,fp);
1310 fwrite(&roll_model[i][j],sizeof(int),1,fp);
1311 fwrite(&twist_model[i][j],sizeof(int),1,fp);
1312 fwrite(&limit_damping[i][j],sizeof(int),1,fp);
1313 fwrite(normal_coeffs[i][j],sizeof(double),4,fp);
1314 fwrite(tangential_coeffs[i][j],sizeof(double),3,fp);
1315 fwrite(roll_coeffs[i][j],sizeof(double),3,fp);
1316 fwrite(twist_coeffs[i][j],sizeof(double),3,fp);
1317 fwrite(&cutoff_type[i][j],sizeof(double),1,fp);
1318 }
1319 }
1320 }
1321 }
1322
1323 /* ----------------------------------------------------------------------
1324 proc 0 reads from restart file, bcasts
1325 ------------------------------------------------------------------------- */
1326
read_restart(FILE * fp)1327 void PairGranular::read_restart(FILE *fp)
1328 {
1329 allocate();
1330 int i,j;
1331 int me = comm->me;
1332 for (i = 1; i <= atom->ntypes; i++) {
1333 for (j = i; j <= atom->ntypes; j++) {
1334 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
1335 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
1336 if (setflag[i][j]) {
1337 if (me == 0) {
1338 utils::sfread(FLERR,&normal_model[i][j],sizeof(int),1,fp,nullptr,error);
1339 utils::sfread(FLERR,&damping_model[i][j],sizeof(int),1,fp,nullptr,error);
1340 utils::sfread(FLERR,&tangential_model[i][j],sizeof(int),1,fp,nullptr,error);
1341 utils::sfread(FLERR,&roll_model[i][j],sizeof(int),1,fp,nullptr,error);
1342 utils::sfread(FLERR,&twist_model[i][j],sizeof(int),1,fp,nullptr,error);
1343 utils::sfread(FLERR,&limit_damping[i][j],sizeof(int),1,fp,nullptr,error);
1344 utils::sfread(FLERR,normal_coeffs[i][j],sizeof(double),4,fp,nullptr,error);
1345 utils::sfread(FLERR,tangential_coeffs[i][j],sizeof(double),3,fp,nullptr,error);
1346 utils::sfread(FLERR,roll_coeffs[i][j],sizeof(double),3,fp,nullptr,error);
1347 utils::sfread(FLERR,twist_coeffs[i][j],sizeof(double),3,fp,nullptr,error);
1348 utils::sfread(FLERR,&cutoff_type[i][j],sizeof(double),1,fp,nullptr,error);
1349 }
1350 MPI_Bcast(&normal_model[i][j],1,MPI_INT,0,world);
1351 MPI_Bcast(&damping_model[i][j],1,MPI_INT,0,world);
1352 MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world);
1353 MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world);
1354 MPI_Bcast(&twist_model[i][j],1,MPI_INT,0,world);
1355 MPI_Bcast(&limit_damping[i][j],1,MPI_INT,0,world);
1356 MPI_Bcast(normal_coeffs[i][j],4,MPI_DOUBLE,0,world);
1357 MPI_Bcast(tangential_coeffs[i][j],3,MPI_DOUBLE,0,world);
1358 MPI_Bcast(roll_coeffs[i][j],3,MPI_DOUBLE,0,world);
1359 MPI_Bcast(twist_coeffs[i][j],3,MPI_DOUBLE,0,world);
1360 MPI_Bcast(&cutoff_type[i][j],1,MPI_DOUBLE,0,world);
1361 }
1362 }
1363 }
1364 }
1365
1366 /* ---------------------------------------------------------------------- */
1367
reset_dt()1368 void PairGranular::reset_dt()
1369 {
1370 dt = update->dt;
1371 }
1372
1373 /* ---------------------------------------------------------------------- */
1374
single(int i,int j,int itype,int jtype,double rsq,double,double,double & fforce)1375 double PairGranular::single(int i, int j, int itype, int jtype,
1376 double rsq, double /* factor_coul */,
1377 double /* factor_lj */, double &fforce)
1378 {
1379 double radi,radj,radsum;
1380 double r,rinv,delx,dely,delz, nx, ny, nz, Reff;
1381 double dR, dR2;
1382 double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
1383 double vtr1,vtr2,vtr3,vrel;
1384 double mi,mj,meff;
1385 double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3;
1386
1387 double knfac, damp_normal, damp_normal_prefactor;
1388 double k_tangential, damp_tangential;
1389 double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit;
1390 double fs, fs1, fs2, fs3;
1391
1392 // for JKR
1393 double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E;
1394 double delta, t0, t1, t2, t3, t4, t5, t6;
1395 double sqrt1, sqrt2, sqrt3;
1396
1397 // rolling
1398 double k_roll, damp_roll;
1399 double rollmag;
1400 double fr, fr1, fr2, fr3;
1401
1402 // twisting
1403 double k_twist, damp_twist, mu_twist;
1404 double signtwist, magtwist, magtortwist, Mtcrit;
1405
1406 double shrmag;
1407 int jnum;
1408 int *jlist;
1409 double *history,*allhistory;
1410
1411 int nall = atom->nlocal + atom->nghost;
1412 if ((i >= nall) || (j >= nall))
1413 error->all(FLERR,"Not enough atoms for pair granular single function");
1414
1415 double *radius = atom->radius;
1416 radi = radius[i];
1417 radj = radius[j];
1418 radsum = radi + radj;
1419 Reff = (radsum > 0.0) ? radi*radj/radsum : 0.0;
1420
1421 bool touchflag;
1422 E = normal_coeffs[itype][jtype][0];
1423 if (normal_model[itype][jtype] == JKR) {
1424 E *= THREEQUARTERS;
1425 R2 = Reff*Reff;
1426 coh = normal_coeffs[itype][jtype][3];
1427 a = cbrt(9.0*MY_PI*coh*R2/(4*E));
1428 delta_pulloff = a*a/Reff - 2*sqrt(MY_PI*coh*a/E);
1429 dist_pulloff = radsum+delta_pulloff;
1430 touchflag = (rsq <= dist_pulloff*dist_pulloff);
1431 } else touchflag = (rsq <= radsum*radsum);
1432
1433 if (!touchflag) {
1434 fforce = 0.0;
1435 for (int m = 0; m < single_extra; m++) svector[m] = 0.0;
1436 return 0.0;
1437 }
1438
1439 double **x = atom->x;
1440 delx = x[i][0] - x[j][0];
1441 dely = x[i][1] - x[j][1];
1442 delz = x[i][2] - x[j][2];
1443 r = sqrt(rsq);
1444 rinv = 1.0/r;
1445
1446 nx = delx*rinv;
1447 ny = dely*rinv;
1448 nz = delz*rinv;
1449
1450 // relative translational velocity
1451
1452 double **v = atom->v;
1453 vr1 = v[i][0] - v[j][0];
1454 vr2 = v[i][1] - v[j][1];
1455 vr3 = v[i][2] - v[j][2];
1456
1457 // normal component
1458
1459 vnnr = vr1*nx + vr2*ny + vr3*nz;
1460 vn1 = nx*vnnr;
1461 vn2 = ny*vnnr;
1462 vn3 = nz*vnnr;
1463
1464 // tangential component
1465
1466 vt1 = vr1 - vn1;
1467 vt2 = vr2 - vn2;
1468 vt3 = vr3 - vn3;
1469
1470 // relative rotational velocity
1471
1472 double **omega = atom->omega;
1473 wr1 = (radi*omega[i][0] + radj*omega[j][0]);
1474 wr2 = (radi*omega[i][1] + radj*omega[j][1]);
1475 wr3 = (radi*omega[i][2] + radj*omega[j][2]);
1476
1477 // meff = effective mass of pair of particles
1478 // if I or J part of rigid body, use body mass
1479 // if I or J is frozen, meff is other particle
1480
1481 double *rmass = atom->rmass;
1482 int *mask = atom->mask;
1483
1484 mi = rmass[i];
1485 mj = rmass[j];
1486 if (fix_rigid) {
1487 // NOTE: ensure mass_rigid is current for owned+ghost atoms?
1488 if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
1489 if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
1490 }
1491
1492 meff = mi*mj / (mi+mj);
1493 if (mask[i] & freeze_group_bit) meff = mj;
1494 if (mask[j] & freeze_group_bit) meff = mi;
1495
1496 delta = radsum - r;
1497 dR = delta*Reff;
1498 if (normal_model[itype][jtype] == JKR) {
1499 dR2 = dR*dR;
1500 t0 = coh*coh*R2*R2*E;
1501 t1 = PI27SQ*t0;
1502 t2 = 8*dR*dR2*E*E*E;
1503 t3 = 4*dR2*E;
1504 // in case sqrt(0) < 0 due to precision issues
1505 sqrt1 = MAX(0, t0*(t1+2*t2));
1506 t4 = cbrt(t1+t2+THREEROOT3*MY_PI*sqrt(sqrt1));
1507 t5 = t3/t4 + t4/E;
1508 sqrt2 = MAX(0, 2*dR + t5);
1509 t6 = sqrt(sqrt2);
1510 sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*MY_PI*R2/(E*t6));
1511 a = INVROOT6*(t6 + sqrt(sqrt3));
1512 a2 = a*a;
1513 knfac = normal_coeffs[itype][jtype][0]*a;
1514 Fne = knfac*a2/Reff - MY_2PI*a2*sqrt(4*coh*E/(MY_PI*a));
1515 } else {
1516 knfac = E;
1517 Fne = knfac*delta;
1518 a = sqrt(dR);
1519 if (normal_model[itype][jtype] != HOOKE) {
1520 Fne *= a;
1521 knfac *= a;
1522 }
1523 if (normal_model[itype][jtype] == DMT)
1524 Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff;
1525 }
1526
1527 if (damping_model[itype][jtype] == VELOCITY) {
1528 damp_normal = 1;
1529 } else if (damping_model[itype][jtype] == MASS_VELOCITY) {
1530 damp_normal = meff;
1531 } else if (damping_model[itype][jtype] == VISCOELASTIC) {
1532 damp_normal = a*meff;
1533 } else if (damping_model[itype][jtype] == TSUJI) {
1534 damp_normal = sqrt(meff*knfac);
1535 } else damp_normal = 0.0;
1536
1537 damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal;
1538 Fdamp = -damp_normal_prefactor*vnnr;
1539
1540 Fntot = Fne + Fdamp;
1541 if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0;
1542
1543 jnum = list->numneigh[i];
1544 jlist = list->firstneigh[i];
1545
1546 if (use_history) {
1547 if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr))
1548 error->one(FLERR,"Pair granular single computation needs history");
1549 allhistory = fix_history->firstvalue[i];
1550 for (int jj = 0; jj < jnum; jj++) {
1551 neighprev++;
1552 if (neighprev >= jnum) neighprev = 0;
1553 if (jlist[neighprev] == j) break;
1554 }
1555 history = &allhistory[size_history*neighprev];
1556 }
1557
1558 //****************************************
1559 // tangential force, including history effects
1560 //****************************************
1561
1562 // For linear, mindlin, mindlin_rescale:
1563 // history = cumulative tangential displacement
1564 //
1565 // For mindlin/force, mindlin_rescale/force:
1566 // history = cumulative tangential elastic force
1567
1568 // tangential component
1569 vt1 = vr1 - vn1;
1570 vt2 = vr2 - vn2;
1571 vt3 = vr3 - vn3;
1572
1573 // relative rotational velocity
1574 wr1 = (radi*omega[i][0] + radj*omega[j][0]);
1575 wr2 = (radi*omega[i][1] + radj*omega[j][1]);
1576 wr3 = (radi*omega[i][2] + radj*omega[j][2]);
1577
1578 // relative tangential velocities
1579 vtr1 = vt1 - (nz*wr2-ny*wr3);
1580 vtr2 = vt2 - (nx*wr3-nz*wr1);
1581 vtr3 = vt3 - (ny*wr1-nx*wr2);
1582 vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
1583 vrel = sqrt(vrel);
1584
1585 if (normal_model[itype][jtype] == JKR) {
1586 F_pulloff = 3*MY_PI*coh*Reff;
1587 Fncrit = fabs(Fne + 2*F_pulloff);
1588 } else if (normal_model[itype][jtype] == DMT) {
1589 F_pulloff = 4*MY_PI*coh*Reff;
1590 Fncrit = fabs(Fne + 2*F_pulloff);
1591 } else {
1592 Fncrit = fabs(Fntot);
1593 }
1594 Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit;
1595
1596 //------------------------------
1597 // tangential forces
1598 //------------------------------
1599 k_tangential = tangential_coeffs[itype][jtype][0];
1600 damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;
1601
1602 if (tangential_history) {
1603 if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) {
1604 k_tangential *= a;
1605 }
1606
1607 shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
1608 history[2]*history[2]);
1609
1610 // tangential forces = history + tangential velocity damping
1611 if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
1612 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
1613 tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
1614 fs1 = -k_tangential*history[0] - damp_tangential*vtr1;
1615 fs2 = -k_tangential*history[1] - damp_tangential*vtr2;
1616 fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
1617 } else {
1618 fs1 = history[0] - damp_tangential*vtr1;
1619 fs2 = history[1] - damp_tangential*vtr2;
1620 fs3 = history[2] - damp_tangential*vtr3;
1621 }
1622
1623 // rescale frictional forces if needed
1624 fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
1625 if (fs > Fscrit) {
1626 if (shrmag != 0.0) {
1627 fs1 *= Fscrit/fs;
1628 fs2 *= Fscrit/fs;
1629 fs3 *= Fscrit/fs;
1630 fs *= Fscrit/fs;
1631 } else fs1 = fs2 = fs3 = fs = 0.0;
1632 }
1633
1634 // classic pair gran/hooke (no history)
1635 } else {
1636 fs = damp_tangential*vrel;
1637 if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel;
1638 else Ft = 0.0;
1639 fs1 = -Ft*vtr1;
1640 fs2 = -Ft*vtr2;
1641 fs3 = -Ft*vtr3;
1642 fs = Ft*vrel;
1643 }
1644
1645 //****************************************
1646 // rolling resistance
1647 //****************************************
1648
1649 if ((roll_model[itype][jtype] != ROLL_NONE)
1650 || (twist_model[itype][jtype] != TWIST_NONE)) {
1651 relrot1 = omega[i][0] - omega[j][0];
1652 relrot2 = omega[i][1] - omega[j][1];
1653 relrot3 = omega[i][2] - omega[j][2];
1654
1655 // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
1656 // this is different from the Marshall papers,
1657 // which use the Bagi/Kuhn formulation
1658 // for rolling velocity (see Wang et al for why the latter is wrong)
1659
1660 vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1;
1661 vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2;
1662 vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3;
1663
1664 int rhist0 = roll_history_index;
1665 int rhist1 = rhist0 + 1;
1666 int rhist2 = rhist1 + 1;
1667
1668 // rolling displacement
1669 rollmag = sqrt(history[rhist0]*history[rhist0] +
1670 history[rhist1]*history[rhist1] +
1671 history[rhist2]*history[rhist2]);
1672
1673 k_roll = roll_coeffs[itype][jtype][0];
1674 damp_roll = roll_coeffs[itype][jtype][1];
1675 fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
1676 fr2 = -k_roll*history[rhist1] - damp_roll*vrl2;
1677 fr3 = -k_roll*history[rhist2] - damp_roll*vrl3;
1678
1679 // rescale frictional displacements and forces if needed
1680 Frcrit = roll_coeffs[itype][jtype][2] * Fncrit;
1681
1682 fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
1683 if (fr > Frcrit) {
1684 if (rollmag != 0.0) {
1685 fr1 *= Frcrit/fr;
1686 fr2 *= Frcrit/fr;
1687 fr3 *= Frcrit/fr;
1688 fr *= Frcrit/fr;
1689 } else fr1 = fr2 = fr3 = fr = 0.0;
1690 }
1691 } else fr1 = fr2 = fr3 = fr = 0.0;
1692
1693 //****************************************
1694 // twisting torque, including history effects
1695 //****************************************
1696
1697 if (twist_model[itype][jtype] != TWIST_NONE) {
1698 // omega_T (eq 29 of Marshall)
1699 magtwist = relrot1*nx + relrot2*ny + relrot3*nz;
1700 if (twist_model[itype][jtype] == TWIST_MARSHALL) {
1701 k_twist = 0.5*k_tangential*a*a;; //eq 32
1702 damp_twist = 0.5*damp_tangential*a*a;
1703 mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2];;
1704 } else {
1705 k_twist = twist_coeffs[itype][jtype][0];
1706 damp_twist = twist_coeffs[itype][jtype][1];
1707 mu_twist = twist_coeffs[itype][jtype][2];
1708 }
1709 // M_t torque (eq 30)
1710 magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;
1711 signtwist = (magtwist > 0) - (magtwist < 0);
1712 Mtcrit = mu_twist*Fncrit; // critical torque (eq 44)
1713 if (fabs(magtortwist) > Mtcrit) {
1714 magtortwist = -Mtcrit * signtwist; // eq 34
1715 } else magtortwist = 0.0;
1716 } else magtortwist = 0.0;
1717
1718 // set force and return no energy
1719
1720 fforce = Fntot*rinv;
1721
1722 // set single_extra quantities
1723
1724 svector[0] = fs1;
1725 svector[1] = fs2;
1726 svector[2] = fs3;
1727 svector[3] = fs;
1728 svector[4] = fr1;
1729 svector[5] = fr2;
1730 svector[6] = fr3;
1731 svector[7] = fr;
1732 svector[8] = magtortwist;
1733 svector[9] = delx;
1734 svector[10] = dely;
1735 svector[11] = delz;
1736 return 0.0;
1737 }
1738
1739 /* ---------------------------------------------------------------------- */
1740
pack_forward_comm(int n,int * list,double * buf,int,int *)1741 int PairGranular::pack_forward_comm(int n, int *list, double *buf,
1742 int /* pbc_flag */, int * /* pbc */)
1743 {
1744 int i,j,m;
1745
1746 m = 0;
1747 for (i = 0; i < n; i++) {
1748 j = list[i];
1749 buf[m++] = mass_rigid[j];
1750 }
1751 return m;
1752 }
1753
1754 /* ---------------------------------------------------------------------- */
1755
unpack_forward_comm(int n,int first,double * buf)1756 void PairGranular::unpack_forward_comm(int n, int first, double *buf)
1757 {
1758 int i,m,last;
1759
1760 m = 0;
1761 last = first + n;
1762 for (i = first; i < last; i++)
1763 mass_rigid[i] = buf[m++];
1764 }
1765
1766 /* ----------------------------------------------------------------------
1767 memory usage of local atom-based arrays
1768 ------------------------------------------------------------------------- */
1769
memory_usage()1770 double PairGranular::memory_usage()
1771 {
1772 double bytes = (double)nmax * sizeof(double);
1773 return bytes;
1774 }
1775
1776 /* ----------------------------------------------------------------------
1777 mixing of Young's modulus (E)
1778 ------------------------------------------------------------------------- */
1779
mix_stiffnessE(double Eii,double Ejj,double poisii,double poisjj)1780 double PairGranular::mix_stiffnessE(double Eii, double Ejj,
1781 double poisii, double poisjj)
1782 {
1783 return 1/((1-poisii*poisii)/Eii+(1-poisjj*poisjj)/Ejj);
1784 }
1785
1786 /* ----------------------------------------------------------------------
1787 mixing of shear modulus (G)
1788 ------------------------------------------------------------------------ */
1789
mix_stiffnessG(double Eii,double Ejj,double poisii,double poisjj)1790 double PairGranular::mix_stiffnessG(double Eii, double Ejj,
1791 double poisii, double poisjj)
1792 {
1793 return 1/((2*(2-poisii)*(1+poisii)/Eii) + (2*(2-poisjj)*(1+poisjj)/Ejj));
1794 }
1795
1796 /* ----------------------------------------------------------------------
1797 mixing of everything else
1798 ------------------------------------------------------------------------- */
1799
mix_geom(double valii,double valjj)1800 double PairGranular::mix_geom(double valii, double valjj)
1801 {
1802 return sqrt(valii*valjj);
1803 }
1804
1805
1806 /* ----------------------------------------------------------------------
1807 compute pull-off distance (beyond contact) for a given radius and atom type
1808 ------------------------------------------------------------------------- */
1809
pulloff_distance(double radi,double radj,int itype,int jtype)1810 double PairGranular::pulloff_distance(double radi, double radj,
1811 int itype, int jtype)
1812 {
1813 double E, coh, a, Reff;
1814 Reff = radi*radj/(radi+radj);
1815 if (Reff <= 0) return 0;
1816 coh = normal_coeffs[itype][jtype][3];
1817 E = normal_coeffs[itype][jtype][0]*THREEQUARTERS;
1818 a = cbrt(9*MY_PI*coh*Reff*Reff/(4*E));
1819 return a*a/Reff - 2*sqrt(MY_PI*coh*a/E);
1820 }
1821
1822 /* ----------------------------------------------------------------------
1823 transfer history during fix/neigh/history exchange
1824 only needed if any history entries i-j are not just negative of j-i entries
1825 ------------------------------------------------------------------------- */
1826
transfer_history(double * source,double * target)1827 void PairGranular::transfer_history(double* source, double* target)
1828 {
1829 for (int i = 0; i < size_history; i++)
1830 target[i] = history_transfer_factors[i]*source[i];
1831 }
1832
1833 /* ----------------------------------------------------------------------
1834 self-interaction range of particle
1835 ------------------------------------------------------------------------- */
1836
atom2cut(int i)1837 double PairGranular::atom2cut(int i)
1838 {
1839 double cut;
1840
1841 cut = atom->radius[i]*2;
1842 if(beyond_contact) {
1843 int itype = atom->type[i];
1844 if(normal_model[itype][itype] == JKR) {
1845 cut += pulloff_distance(cut, cut, itype, itype);
1846 }
1847 }
1848
1849 return cut;
1850 }
1851
1852 /* ----------------------------------------------------------------------
1853 maximum interaction range for two finite particles
1854 ------------------------------------------------------------------------- */
1855
radii2cut(double r1,double r2)1856 double PairGranular::radii2cut(double r1, double r2)
1857 {
1858 double cut = 0.0;
1859
1860 if(beyond_contact) {
1861 int n = atom->ntypes;
1862 double temp;
1863
1864 // Check all combinations of i and j to find theoretical maximum pull off distance
1865 for(int i = 0; i < n; i++){
1866 for(int j = 0; j < n; j++){
1867 if(normal_model[i][j] == JKR) {
1868 temp = pulloff_distance(r1, r2, i, j);
1869 if(temp > cut) cut = temp;
1870 }
1871 }
1872 }
1873 }
1874
1875 cut += r1 + r2;
1876
1877 return cut;
1878 }
1879