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 This file is from LAMMPS, but has been modified. Copyright for
36 modification:
37
38 Copyright 2013- TU Graz
39
40 Copyright of original file:
41 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
42 http://lammps.sandia.gov, Sandia National Laboratories
43 Steve Plimpton, sjplimp@sandia.gov
44
45 Copyright (2003) Sandia Corporation. Under the terms of Contract
46 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
47 certain rights in this software. This software is distributed under
48 the GNU General Public License.
49 ------------------------------------------------------------------------- */
50
51 #include <cmath>
52 #include "force.h"
53 #include <stdlib.h>
54 #include <string.h>
55 #include "neighbor.h"
56 #include "neigh_list.h"
57 #include "neigh_request.h"
58 #include "pair_gran.h"
59 #include "fix_fiber_spring_simple.h"
60 #include "atom.h"
61 #include "error.h"
62 #include "modify.h"
63 #include "update.h"
64 #include "domain.h"
65 #include "error.h"
66 #include "group.h"
67 #include "memory.h"
68 #include "pair_gran.h"
69 #include "math_extra.h"
70 #include "fix_property_global.h"
71 #include "fix_property_atom.h"
72 #include "fix_scalar_transport_equation.h"
73 #include "properties.h"
74 #include "respa.h"
75
76 using namespace LAMMPS_NS;
77 using namespace FixConst;
78
79 #define SMALL 1.0e-10
80
81 /* ---------------------------------------------------------------------- */
82
FixFiberSpringSimple(LAMMPS * lmp,int narg,char ** arg)83 FixFiberSpringSimple::FixFiberSpringSimple(LAMMPS *lmp, int narg, char **arg) :
84 Fix(lmp, narg, arg)
85 {
86 if (narg < 9) error->all(FLERR,"Illegal fix fiber/simpleSpring command");
87
88 scalar_flag = 1;
89 vector_flag = 1;
90 size_vector = 4;
91 global_freq = 1;
92 extscalar = 1;
93 extvector = 1;
94
95 group2 = NULL;
96
97 Yeff = NULL;
98
99 if (strcmp(arg[3],"tether") == 0) {
100 if (narg != 9) error->all(FLERR,"Illegal fix spring command");
101 k_spring = atof(arg[4]);
102 xflag = yflag = zflag = 1;
103 if (strcmp(arg[5],"NULL") == 0) xflag = 0;
104 else xc = atof(arg[5]);
105 if (strcmp(arg[6],"NULL") == 0) yflag = 0;
106 else yc = atof(arg[6]);
107 if (strcmp(arg[7],"NULL") == 0) zflag = 0;
108 else zc = atof(arg[7]);
109 r0 = atof(arg[8]);
110 if (r0 < 0) error->all(FLERR,"R0 < 0 for fix spring command");
111
112 } else if (strcmp(arg[3],"couple") == 0) {
113 if (narg != 10) error->all(FLERR,"Illegal fix spring command");
114
115 int n = strlen(arg[4]) + 1;
116 group2 = new char[n];
117 strcpy(group2,arg[4]);
118 igroup2 = group->find(arg[4]);
119 if (igroup2 == -1)
120 error->all(FLERR,"Fix spring couple group ID does not exist");
121 if (igroup2 == igroup)
122 error->all(FLERR,"Two groups cannot be the same in fix spring couple");
123 group2bit = group->bitmask[igroup2];
124
125 k_spring = atof(arg[5]);
126 xflag = yflag = zflag = 1;
127 if (strcmp(arg[6],"NULL") == 0) xflag = 0;
128 else xc = atof(arg[6]);
129 if (strcmp(arg[7],"NULL") == 0) yflag = 0;
130 else yc = atof(arg[7]);
131 if (strcmp(arg[8],"NULL") == 0) zflag = 0;
132 else zc = atof(arg[8]);
133 r0 = atof(arg[9]);
134 if (r0 < 0) error->all(FLERR,"R0 < 0 for fix spring command");
135
136 } else error->all(FLERR,"Illegal fix spring command");
137
138 ftotal[0] = ftotal[1] = ftotal[2] = ftotal[3] = 0.0;
139 }
140
141 /* ---------------------------------------------------------------------- */
142
~FixFiberSpringSimple()143 FixFiberSpringSimple::~FixFiberSpringSimple()
144 {
145 delete [] group2;
146 memory->destroy(Yeff);
147 }
148
149 /* ---------------------------------------------------------------------- */
150
setmask()151 int FixFiberSpringSimple::setmask()
152 {
153 int mask = 0;
154 mask |= POST_FORCE;
155 mask |= THERMO_ENERGY;
156 mask |= POST_FORCE_RESPA;
157 mask |= MIN_POST_FORCE;
158 return mask;
159 }
160
161 /* ---------------------------------------------------------------------- */
162
init()163 void FixFiberSpringSimple::init()
164 {
165 // recheck that group 2 has not been deleted
166
167 if (group2) {
168 igroup2 = group->find(group2);
169 if (igroup2 == -1)
170 error->all(FLERR,"Fix spring couple group ID does not exist");
171 group2bit = group->bitmask[igroup2];
172 }
173
174 //Get pointer to the fixes that have the material properties
175 int max_type = 1;
176 memory->destroy(Yeff);
177 memory->create(Yeff,max_type+1,max_type+1,"Yeff");
178
179 Y1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulus","property/global","peratomtype",max_type,0,force->pair_style));
180 v1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("poissonsRatio","property/global","peratomtype",max_type,0,force->pair_style));
181 charVel1=static_cast<FixPropertyGlobal*>(modify->find_fix_property("characteristicVelocity","property/global","scalar",0,0,force->pair_style));
182
183 //pre-calculate parameters for possible contact material combinations
184 for(int i=1;i< max_type+1; i++)
185 {
186 for(int j=1;j<max_type+1;j++)
187 {
188 double Yi=Y1->compute_vector(i-1);
189 double Yj=Y1->compute_vector(j-1);
190 double vi=v1->compute_vector(i-1);
191 double vj=v1->compute_vector(j-1);
192 Yeff[i][j] = 1./((1.-pow(vi,2.))/Yi+(1.-pow(vj,2.))/Yj);
193
194 }
195 }
196 charVel = charVel1->compute_scalar();
197
198 if(!force->pair_match("gran", 0)) error->all(FLERR,"Please use a granular pair style for fix fiber/simpleSpring");
199 pair_gran = static_cast<PairGran*>(force->pair_match("gran", 0));
200
201 masstotal = group->mass(igroup);
202 // if (styleflag == COUPLE) masstotal2 = group->mass(igroup2);
203
204 if (strstr(update->integrate_style,"respa"))
205 nlevels_respa = ((Respa *) update->integrate)->nlevels;
206 }
207
208 /* ---------------------------------------------------------------------- */
209
setup(int vflag)210 void FixFiberSpringSimple::setup(int vflag)
211 {
212 if (strstr(update->integrate_style,"verlet"))
213 post_force(vflag);
214 else {
215 ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
216 post_force_respa(vflag,nlevels_respa-1,0);
217 ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
218 }
219 }
220
221 /* ---------------------------------------------------------------------- */
222
min_setup(int vflag)223 void FixFiberSpringSimple::min_setup(int vflag)
224 {
225 post_force(vflag);
226 }
227
228 /* ---------------------------------------------------------------------- */
229
post_force(int vflag)230 void FixFiberSpringSimple::post_force(int vflag)
231 {
232
233 //calculated from the material properties
234 double deltaEquilibrium= 0.75;
235 double snapDistance = 0.90;
236 double kn = 1000;
237
238 int i,j,ii,jj,inum,jnum,itype,jtype;
239 double xtmp,ytmp,ztmp,delx,dely,delz;
240 double radi,radj,radsum,rsq,r,rinv;
241 double meff,ccel;
242 int *ilist,*jlist,*numneigh,**firstneigh;
243 //int *contact_flag,**first_contact_flag;
244 //double *all_contact_hist,**first_contact_hist;
245
246 double **x = atom->x;
247 double **f = atom->f;
248 double *radius = atom->radius;
249 double *rmass = atom->rmass;
250 double *mass = atom->mass;
251 int *type = atom->type;
252 int nlocal = atom->nlocal;
253
254 inum = pair_gran->list->inum;
255 ilist = pair_gran->list->ilist;
256 numneigh = pair_gran->list->numneigh;
257 firstneigh = pair_gran->list->firstneigh;
258
259 //************** LOOP OVER PARTICLES *********
260 // loop over neighbors of my atoms
261
262 for (ii = 0; ii < inum; ii++)
263 {
264 i = ilist[ii];
265 xtmp = x[i][0];
266 ytmp = x[i][1];
267 ztmp = x[i][2];
268 radi = radius[i];
269 //contact_flag = first_contact_flag[i];
270 //all_contact_hist = first_contact_hist[i];
271 jlist = firstneigh[i];
272 jnum = numneigh[i];
273
274 for (jj = 0; jj < jnum; jj++)
275 {
276 j = jlist[jj];
277 j &= NEIGHMASK;
278
279 delx = xtmp - x[j][0]; //delta vector pointing from j to i
280 dely = ytmp - x[j][1];
281 delz = ztmp - x[j][2];
282 rsq = delx*delx + dely*dely + delz*delz;
283 r = sqrt(rsq);
284
285 radj = radius[j];
286 radsum = radi + radj;
287
288 if (r <= radsum*snapDistance) //check if there is sufficient overlap
289 {
290 // fprintf(screen, "SNAPPED!! \n");
291
292 //double deltan=radsum-r;
293 rinv = 1.0/r;
294
295 double mi,mj;
296 if (rmass)
297 {
298 mi = rmass[i];
299 mj = rmass[j];
300 }
301 else
302 {
303 itype = type[i];
304 jtype = type[j];
305 mi = mass[itype];
306 mj = mass[jtype];
307 }
308
309 meff = mi*mj/(mi+mj);
310 //deriveContactModelParams(i,j,meff,deltan,kn,kt,gamman,gammat,xmu,rmu);
311 deriveNormalStiffness(i,j,meff,kn);
312
313 //Calculate the force and subtract to model bonding
314 ccel = kn * radsum*(1.0-deltaEquilibrium) * rinv ;
315 // ccel = kn*(radsum-r)*rinv;
316 f[i][0] -= delx*ccel;
317 f[i][1] -= dely*ccel;
318 f[i][2] -= delz*ccel;
319
320 if (j < nlocal)
321 {
322 f[j][0] += delx*ccel;
323 f[j][1] += dely*ccel;
324 f[j][2] += delz*ccel;
325 }
326
327 // fprintf(screen, "r: %3.3g, critDistance: %3.3g \n", r, radsum*snapDistance);
328 // fprintf(screen, "NormalStiffness: %3.3g \n", kn);
329 // fprintf(screen, "i: %2.1d f: %3.3g %3.3g %3.3g \n", i, f[i][0],f[i][1],f[i][2]);
330
331 }
332
333 }
334 }
335
336 }
337
338 /* ---------------------------------------------------------------------- */
339
post_force_respa(int vflag,int ilevel,int iloop)340 void FixFiberSpringSimple::post_force_respa(int vflag, int ilevel, int iloop)
341 {
342 if (ilevel == nlevels_respa-1) post_force(vflag);
343 }
344
345 /* ---------------------------------------------------------------------- */
346
min_post_force(int vflag)347 void FixFiberSpringSimple::min_post_force(int vflag)
348 {
349 post_force(vflag);
350 }
351
352 /* ----------------------------------------------------------------------
353 energy of stretched spring
354 ------------------------------------------------------------------------- */
355
compute_scalar()356 double FixFiberSpringSimple::compute_scalar()
357 {
358 return espring;
359 }
360
361 /* ----------------------------------------------------------------------
362 return components of total spring force on fix group
363 ------------------------------------------------------------------------- */
364
compute_vector(int n)365 double FixFiberSpringSimple::compute_vector(int n)
366 {
367 return ftotal[n];
368 }
369
deriveNormalStiffness(int & ip,int & jp,double & meff,double & kn)370 inline void FixFiberSpringSimple::deriveNormalStiffness(int &ip, int &jp,double &meff, double &kn)
371 {
372
373 int itype = atom->type[ip];
374 int jtype = atom->type[jp];
375 double rj = atom->radius[jp];
376 double ri = atom->radius[ip];
377 double reff=ri*rj/(ri+rj);
378
379 kn = 16./15.*sqrt(reff)*(Yeff[itype][jtype])*pow(15.*meff*charVel*charVel/(16.*sqrt(reff)*Yeff[itype][jtype]),0.2);
380
381 // convert Kn and Kt from pressure units to force/distance^2
382 kn /= force->nktv2p;
383
384 return;
385
386 }
387
388