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     Christoph Kloss (DCS Computing GmbH, Linz)
37     Christoph Kloss (JKU Linz)
38     Richard Berger (JKU Linz)
39 
40     Copyright 2012-     DCS Computing GmbH, Linz
41     Copyright 2009-2012 JKU Linz
42 ------------------------------------------------------------------------- */
43 
44 #ifdef NORMAL_MODEL
45 NORMAL_MODEL(HERTZ_STIFFNESS,hertz/stiffness,4)
46 #else
47 #ifndef NORMAL_MODEL_HERTZ_STIFFNESS_H_
48 #define NORMAL_MODEL_HERTZ_STIFFNESS_H_
49 #include "contact_models.h"
50 #include "normal_model_base.h"
51 #include "global_properties.h"
52 #include <cmath>
53 
54 namespace LIGGGHTS {
55 namespace ContactModels
56 {
57   template<>
58   class NormalModel<HERTZ_STIFFNESS> : public NormalModelBase
59   {
60   public:
61     NormalModel(LAMMPS * lmp, IContactHistorySetup * hsetup, class ContactModelBase * c) :
62       NormalModelBase(lmp, hsetup, c),
63       k_n(NULL),
64       k_t(NULL),
65       gamma_n(NULL),
66       gamma_t(NULL),
67       tangential_damping(false),
68       limitForce(false),
69       displayedSettings(false)
70     {
71 
72     }
73 
74     void registerSettings(Settings & settings)
75     {
76       settings.registerOnOff("tangential_damping", tangential_damping, true);
77       settings.registerOnOff("limitForce", limitForce);
78     }
79 
80     inline void postSettings(IContactHistorySetup * hsetup, ContactModelBase *cmb) {}
81 
82     void connectToProperties(PropertyRegistry & registry)
83     {
84       registry.registerProperty("k_n", &MODEL_PARAMS::createKn);
85       registry.registerProperty("k_t", &MODEL_PARAMS::createKt);
86       registry.registerProperty("gamma_n", &MODEL_PARAMS::createGamman);
87       registry.registerProperty("gamma_t", &MODEL_PARAMS::createGammat);
88 
89       registry.connect("k_n", k_n,"model hertz/stiffness");
90       registry.connect("k_t", k_t,"model hertz/stiffness");
91       registry.connect("gamma_n", gamma_n,"model hertz/stiffness");
92       registry.connect("gamma_t", gamma_t,"model hertz/stiffness");
93 
94       // error checks on coarsegraining
95       if(force->cg_active())
96         error->cg(FLERR,"model hertz/stiffness");
97     }
98 
99     // effective exponent for stress-strain relationship
100 
101     inline double stressStrainExponent()
102     {
103       return 1.5;
104     }
105 
106     inline void surfacesIntersect(SurfacesIntersectData & sidata, ForceData & i_forces, ForceData & j_forces)
107     {
108       const int itype = sidata.itype;
109       const int jtype = sidata.jtype;
110       const double meff = sidata.meff;
111       double reff = sidata.is_wall ? sidata.radi : (sidata.radi*sidata.radj/(sidata.radi+sidata.radj));
112 #ifdef SUPERQUADRIC_ACTIVE_FLAG
113       if(sidata.is_non_spherical && atom->superquadric_flag) {
114           reff = sidata.reff;
115       }
116 #endif
117 
118       const double polyhertz = sqrt(reff*sidata.deltan);
119       double kn = polyhertz*k_n[itype][jtype];
120       double kt = polyhertz*k_t[itype][jtype];
121       const double gamman = polyhertz*meff*gamma_n[itype][jtype];
122       const double gammat = tangential_damping ? polyhertz*meff*gamma_t[itype][jtype] : 0.0;
123 
124       if(!displayedSettings)
125       {
126         displayedSettings = true;
127 
128         /*
129         if(limitForce)
130             if(0 == comm->me) fprintf(screen," NormalModel<HERTZ_STIFFNESS>: will limit normal force.\n");
131         */
132       }
133       // convert Kn and Kt from pressure units to force/distance^2
134       kn /= force->nktv2p;
135       kt /= force->nktv2p;
136 
137       const double Fn_damping = -gamman*sidata.vn;
138       const double Fn_contact = kn*sidata.deltan;
139       double Fn               = Fn_damping + Fn_contact;
140 
141       //limit force to avoid the artefact of negative repulsion force
142       if(limitForce && (Fn<0.0) )
143       {
144           Fn = 0.0;
145       }
146 
147       sidata.Fn = Fn;
148       sidata.kn = kn;
149       sidata.kt = kt;
150       sidata.gamman = gamman;
151       sidata.gammat = gammat;
152 
153       #ifdef NONSPHERICAL_ACTIVE_FLAG
154           double torque_i[3] = {0.0, 0.0, 0.0}; //initialized here with zeros to avoid compiler warnings
155           double Fn_i[3] = { Fn * sidata.en[0], Fn * sidata.en[1], Fn * sidata.en[2]};
156           if(sidata.is_non_spherical) {
157             double xci[3];
158             vectorSubtract3D(sidata.contact_point, atom->x[sidata.i], xci);
159             vectorCross3D(xci, Fn_i, torque_i);
160           }
161       #endif
162       // apply normal force
163       if(sidata.is_wall) {
164         const double Fn_ = Fn * sidata.area_ratio;
165         i_forces.delta_F[0] += Fn_ * sidata.en[0];
166         i_forces.delta_F[1] += Fn_ * sidata.en[1];
167         i_forces.delta_F[2] += Fn_ * sidata.en[2];
168         #ifdef NONSPHERICAL_ACTIVE_FLAG
169                 if(sidata.is_non_spherical) {
170                   //for non-spherical particles normal force can produce torque!
171                   i_forces.delta_torque[0] += torque_i[0];
172                   i_forces.delta_torque[1] += torque_i[1];
173                   i_forces.delta_torque[2] += torque_i[2];
174                 }
175         #endif
176       } else {
177         i_forces.delta_F[0] += sidata.Fn * sidata.en[0];
178         i_forces.delta_F[1] += sidata.Fn * sidata.en[1];
179         i_forces.delta_F[2] += sidata.Fn * sidata.en[2];
180 
181         j_forces.delta_F[0] += -i_forces.delta_F[0];
182         j_forces.delta_F[1] += -i_forces.delta_F[1];
183         j_forces.delta_F[2] += -i_forces.delta_F[2];
184         #ifdef NONSPHERICAL_ACTIVE_FLAG
185                 if(sidata.is_non_spherical) {
186                   //for non-spherical particles normal force can produce torque!
187                   double xcj[3], torque_j[3];
188                   double Fn_j[3] = { -Fn_i[0], -Fn_i[1], -Fn_i[2]};
189                   vectorSubtract3D(sidata.contact_point, atom->x[sidata.j], xcj);
190                   vectorCross3D(xcj, Fn_j, torque_j);
191 
192                   i_forces.delta_torque[0] += torque_i[0];
193                   i_forces.delta_torque[1] += torque_i[1];
194                   i_forces.delta_torque[2] += torque_i[2];
195 
196                   j_forces.delta_torque[0] += torque_j[0];
197                   j_forces.delta_torque[1] += torque_j[1];
198                   j_forces.delta_torque[2] += torque_j[2];
199                 }
200         #endif
201       }
202     }
203 
204     void surfacesClose(SurfacesCloseData&, ForceData&, ForceData&){}
205     void beginPass(SurfacesIntersectData&, ForceData&, ForceData&){}
206     void endPass(SurfacesIntersectData&, ForceData&, ForceData&){}
207 
208   protected:
209     double ** k_n;
210     double ** k_t;
211     double ** gamma_n;
212     double ** gamma_t;
213 
214     bool tangential_damping;
215     bool limitForce;
216     bool displayedSettings;
217   };
218 }
219 }
220 #endif // NORMAL_MODEL_HERTZ_STIFFNESS_H_
221 #endif
222