JKR Model Implementation
See this page for the documentation of this contact model.
contactmodeljkr.h
1// contactmodeljkr.h
2#include "contactmodel/src/contactmodelmechanical.h"
3
4#include "module/interface/itablelist.h"
5#include "module/interface/itable.h"
6
7#ifdef JKR_LIB
8# define JKR_EXPORT EXPORT_TAG
9#elif defined(NO_MODEL_IMPORT)
10# define JKR_EXPORT
11#else
12# define JKR_EXPORT IMPORT_TAG
13#endif
14
15namespace cmodelsxd {
16 using namespace itasca;
17
18 class ContactModelJKR : public ContactModelMechanical {
19 public:
20 // Constructor: Set default values for contact model properties.
21 JKR_EXPORT ContactModelJKR();
22 JKR_EXPORT ContactModelJKR(const ContactModelJKR &) noexcept;
23 JKR_EXPORT const ContactModelJKR & operator=(const ContactModelJKR &);
24 JKR_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
25 // Destructor, called when contact is deleted: free allocated memory, etc.
26 JKR_EXPORT virtual ~ContactModelJKR();
27 // Contact model name (used as keyword for commands and FISH).
28 string getName() const override { return "jkr"; }
29 // The index provides a quick way to determine the type of contact model.
30 // Each type of contact model in PFC must have a unique index; this is assigned
31 // by PFC when the contact model is loaded. This index should be set to -1
32 void setIndex(int i) override { index_ = i; }
33 int getIndex() const override { return index_; }
34 // Contact model version number (e.g., MyModel05_1). The version number can be
35 // accessed during the save-restore operation (within the archive method,
36 // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for
37 // future modifications to the contact model data structure.
38 uint32 getMinorVersion() const override;
39 // Copy the state information to a newly created contact model.
40 // Provide access to state information, for use by copy method.
41 void copy(const ContactModel *c) override;
42 // Provide save-restore capability for the state information.
43 void archive(ArchiveStream &) override;
44 // Enumerator for the properties.
45 enum PropertyKeys {
46 kwShear = 1
47 , kwPoiss
48 , kwFric
49 , kwjkrF
50 , kwjkrS
51 , kwDpNRatio
52 , kwDpSRatio
53 , kwDpMode
54 , kwDpF
55 , kwResFric
56 , kwResMoment
57 , kwResS
58 , kwKsFac
59 , kwSurfAdh
60 , kwActiveMode
61 , kwA0
62 , kwPullOff
63 , kwTearOff
64 };
65 // Contact model property names in a comma separated list. The order corresponds with
66 // the order of the PropertyKeys enumerator above. One can visualize any of these
67 // properties in PFC automatically.
68 string getProperties() const override {
69 return "jkr_shear"
70 ",jkr_poiss"
71 ",fric"
72 ",jkr_force"
73 ",jkr_slip"
74 ",dp_nratio"
75 ",dp_sratio"
76 ",dp_mode"
77 ",dp_force"
78 ",rr_fric"
79 ",rr_moment"
80 ",rr_slip"
81 ",ks_fac"
82 ",surf_adh"
83 ",active_mode"
84 ",a0"
85 ",pull_off_f"
86 ",tear_off_d";
87 }
88 // Enumerator for the energies.
89 enum EnergyKeys {
90 kwEStrain = 1
91 , kwERRStrain
92 , kwESlip
93 , kwERRSlip
94 , kwEDashpot
95 };
96 // Contact model energy names in a comma separated list. The order corresponds with
97 // the order of the EnergyKeys enumerator above.
98 string getEnergies() const override {
99 return "energy-strain-jkr"
100 ",energy-rrstrain"
101 ",energy-slip"
102 ",energy-rrslip"
103 ",energy-dashpot";
104 }
105 // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
106 double getEnergy(uint32 i) const override;
107 // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1)
108 // returns whether or not the estrain energy is accumulated which is false).
109 bool getEnergyAccumulate(uint32 i) const override;
110 // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
111 void setEnergy(uint32 i, const double &d) override; // Base 1
112 // Activate the energy. This is only called if the energy tracking is enabled.
113 void activateEnergy() override { if (energies_) return; energies_ = NEW Energies(); }
114 // Returns whether or not the energy tracking has been enabled for this contact.
115 bool getEnergyActivated() const override { return (energies_ != 0); }
116
117 // Enumerator for contact model related FISH callback events.
118 enum FishCallEvents {
119 fActivated = 0
120 , fSlipChange
121 };
122 // Contact model FISH callback event names in a comma separated list. The order corresponds with
123 // the order of the FishCallEvents enumerator above.
124 string getFishCallEvents() const override {
125 return
126 "contact_activated"
127 ",slip_change";
128 }
129
130 // Return the specified contact model property.
131 base::Property getProperty(uint32 i, const IContact *) const override;
132 // The return value denotes whether or not the property corresponds to the global
133 // or local coordinate system (TRUE: global system, FALSE: local system). The
134 // local system is the contact-plane system (nst) defined as follows.
135 // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
136 // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
137 // and tc are unit vectors in directions of the nst axes.
138 // This is used when rendering contact model properties that are vectors.
139 bool getPropertyGlobal(uint32 i) const override;
140 // Set the specified contact model property, ensuring that it is of the correct type
141 // and within the correct range --- if not, then throw an exception.
142 // The return value denotes whether or not the update has affected the timestep
143 // computation (by having modified the translational or rotational tangent stiffnesses).
144 // If true is returned, then the timestep will be recomputed.
145 bool setProperty(uint32 i, const base::Property &v, IContact *c) override;
146 // The return value denotes whether or not the property is read-only
147 // (TRUE: read-only, FALSE: read-write).
148 bool getPropertyReadOnly(uint32 i) const override;
149
150 // The return value denotes whether or not the property is inheritable
151 // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
152 // the endPropertyUpdated method.
153 bool supportsInheritance(uint32 i) const override;
154 // Return whether or not inheritance is enabled for the specified property.
155 bool getInheritance(uint32 i) const override { assert(i < 32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
156 // Set the inheritance flag for the specified property.
157 void setInheritance(uint32 i, bool b) override { assert(i < 32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
158
159 // Prepare for entry into ForceDispLaw. The validate function is called when:
160 // (1) the contact is created, (2) a property of the contact that returns a true via
161 // the setProperty method has been modified and (3) when a set of cycles is executed
162 // via the {cycle N} command.
163 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
164 bool validate(ContactModelMechanicalState *state, const double ×tep) override;
165 // The endPropertyUpdated method is called whenever a surface property (with a name
166 // that matches an inheritable contact model property name) of one of the contacting
167 // pieces is modified. This allows the contact model to update its associated
168 // properties. The return value denotes whether or not the update has affected
169 // the time step computation (by having modified the translational or rotational
170 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
171 bool endPropertyUpdated(const string &name, const IContactMechanical *c) override;
172 // The forceDisplacementLaw function is called during each cycle. Given the relative
173 // motion of the two contacting pieces (via
174 // state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
175 // state->relativeAngularIncrement_ (Dtt, Dtbs, Dtbt)
176 // Ddn : relative normal-displacement increment, Ddn > 0 is opening
177 // Ddss : relative shear-displacement increment (s-axis component)
178 // Ddst : relative shear-displacement increment (t-axis component)
179 // Dtt : relative twist-rotation increment
180 // Dtbs : relative bend-rotation increment (s-axis component)
181 // Dtbt : relative bend-rotation increment (t-axis component)
182 // The relative displacement and rotation increments:
183 // Dd = Ddn*nc + Ddss*sc + Ddst*tc
184 // Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
185 // where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
186 // [see {Table 1: Contact State Variables} in PFC Model Components:
187 // Contacts and Contact Models: Contact Resolution]
188 // ) and the contact properties, this function must update the contact force and
189 // moment.
190 // The force_ is acting on piece 2, and is expressed in the local coordinate system
191 // (defined in getPropertyGlobal) such that the first component positive denotes
192 // compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
193 // in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to
194 // get the total moment.
195 // The return value indicates the contact activity status (TRUE: active, FALSE:
196 // inactive) during the next cycle.
197 // Additional information:
198 // * If state->activated() is true, then the contact has just become active (it was
199 // inactive during the previous time step).
200 // * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
201 // the forceDispLaw handle the case of {state->canFail_ == true}.
202 bool forceDisplacementLaw(ContactModelMechanicalState *state, const double ×tep) override;
203 // The getEffectiveXStiffness functions return the translational and rotational
204 // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
205 // the translational tangent shear stiffness is zero (but this stiffness reduction
206 // is typically ignored when computing a stable time step). If the contact model
207 // includes a dashpot, then the translational stiffnesses must be increased (see
208 // Potyondy (2009)).
209 // [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
210 // Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
211 // December 7, 2009.]
212 DVect2 getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
213 DAVect getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_; }
214
215 // Return a new instance of the contact model. This is used in the CMAT
216 // when a new contact is created.
217 ContactModelJKR *clone() const override { return NEW ContactModelJKR(); }
218 // The getActivityDistance function is called by the contact-resolution logic when
219 // the CMAT is modified. Return value is the activity distance used by the
220 // checkActivity function.
221 double getActivityDistance() const override { return (rgap_ + distance_active_); } // initially 'distance_active_' is zero so that the contat becomes active when the particles touch (zerp overlap)
222 // when the contact becomes active,it is set so that the contact only becomes inactive at this tear_off_ distance
223 // The isOKToDelete function is called by the contact-resolution logic when...
224 // Return value indicates whether or not the contact may be deleted.
225 // If TRUE, then the contact may be deleted when it is inactive.
226 // If FALSE, then the contact may not be deleted (under any condition).
227 bool isOKToDelete() const override { return !isBonded(); }
228 // Zero the forces and moments stored in the contact model. This function is called
229 // when the contact becomes inactive.
230 void resetForcesAndMoments() override {
231 jkr_F(DVect(0.0)); dp_F(DVect(0.0));
232 res_M(DAVect(0.0));
233 if (energies_) {
234 //energies_->estrain_ = 0.0;
235 energies_->errstrain_ = 0.0;
236 }
237 distance_active_ = 0.0; // reset
238 }
239 void setForce(const DVect &v, IContact *c) override;
240 void setArea(const double&) override { throw Exception("The setArea method cannot be used with the JKR contact model."); }
241 double getArea() const override { return 0.0; }
242 // The checkActivity function is called by the contact-resolution logic when...
243 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
244 // A contact with the arrlinear model is active if the surface gap is less than
245 // or equal to the attraction range (a_d0_).
246 bool checkActivity(const double &gap) override { return gap <= (rgap_ + distance_active_); } // initially 'distance_active_' is zero so that the contat becomes active when the particles touch (zerp overlap)
247 // when the contact becomes active,it is set so that the contact only becomes inactive at this tear_off_ distance
248 // Returns the sliding state (FALSE is returned if not implemented).
249 bool isSliding() const override { return jkr_S_; }
250 // Returns the bonding state (FALSE is returned if not implemented).
251 bool isBonded() const override { return false; }
252 // Both of these methods are called only for contacts with facets where the wall
253 // resolution scheme is set the full. In such cases one might wish to propagate
254 // contact state information (e.g., shear force) from one active contact to another.
255 // See the Faceted Wall section in the documentation.
256
257 // Return the total force that the contact model holds.
258 DVect getForce() const override;
259 // Return the total moment on 1 that the contact model holds
260 DAVect getMomentOn1(const IContactMechanical *) const override;
261 // Return the total moment on 1 that the contact model holds
262 DAVect getMomentOn2(const IContactMechanical *) const override;
263
264 DAVect getModelMomentOn1() const override;
265 DAVect getModelMomentOn2() const override;
266
267 // Used to efficiently get properties from the contact model for the object pane.
268 // List of properties for the object pane, comma separated.
269 // All properties will be cast to doubles for comparison. No other comparisons
270 // are supported. This may not be the same as the entire property list.
271 // Return property name and type for plotting.
272 void objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *) const override;
273 // All properties cast to doubles - this is what can be compared.
274 void objectPropValues(std::vector<double> *,const IContact *) const override;
275
276 // Methods to get and set properties.
277 const double & shear() const { return shear_; }
278 void shear(const double &d) { shear_ = d; }
279 const double & poiss() const { return poiss_; }
280 void poiss(const double &d) { poiss_ = d; }
281 const double & fric() const { return fric_; }
282 void fric(const double &d) { fric_ = d; }
283 const DVect & jkr_F() const { return jkr_F_; }
284 void jkr_F(const DVect &f) { jkr_F_ = f; }
285 bool jkr_S() const { return jkr_S_; }
286 void jkr_S(bool b) { jkr_S_ = b; }
287 const double & rgap() const { return rgap_; }
288 void rgap(const double &d) { rgap_ = d; }
289
290 const double & kn() const { return kn_; }
291 void kn(const double &d) { kn_ = d; }
292 const double & ks() const { return ks_; }
293 void ks(const double &d) { ks_ = d; }
294 const double & ks_fac() const { return ks_fac_; }
295 void ks_fac(const double &d) { ks_fac_ = d; }
296 const double & surf_adh() const { return surf_adh_; }
297 void surf_adh(const double &d) { surf_adh_ = d; }
298 const int & active_mode() const { return active_mode_; }
299 void active_mode(const int &d) { active_mode_ = d; }
300 const double & distance_active() const { return distance_active_; }
301 void distance_active(const double &d) { distance_active_ = d; }
302 const double & a0() const { return a0_; }
303 void a0(const double &d) { a0_ = d; }
304 const double & pull_off() const { return pull_off_f_; }
305 void pull_off(const double &d) { pull_off_f_ = d; }
306 const double & tear_off() const { return tear_off_d_; }
307 void tear_off(const double &d) { tear_off_d_ = d; }
308 const double & r_hertz() const { return r_hertz_; }
309 void r_hertz(const double &d) { r_hertz_ = d; }
310 const double & c1() const { return c1_; }
311 void c1(const double &d) { c1_ = d; }
312
313 bool hasDamping() const { return dpProps_ ? true : false; }
314 double dp_nratio() const { return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0); }
315 void dp_nratio(const double &d) { if (!hasDamping()) return; dpProps_->dp_nratio_ = d; }
316 double dp_sratio() const { return hasDamping() ? dpProps_->dp_sratio_ : 0.0; }
317 void dp_sratio(const double &d) { if (!hasDamping()) return; dpProps_->dp_sratio_ = d; }
318 int dp_mode() const { return hasDamping() ? dpProps_->dp_mode_ : -1; }
319 void dp_mode(int i) { if (!hasDamping()) return; dpProps_->dp_mode_ = i; }
320 DVect dp_F() const { return hasDamping() ? dpProps_->dp_F_ : DVect(0.0); }
321 void dp_F(const DVect &f) { if (!hasDamping()) return; dpProps_->dp_F_ = f; }
322
323 bool hasEnergies() const { return energies_ ? true : false; }
324 double estrain() const { return hasEnergies() ? energies_->estrain_ : 0.0; }
325 void estrain(const double &d) { if (!hasEnergies()) return; energies_->estrain_ = d; }
326 double errstrain() const { return hasEnergies() ? energies_->errstrain_ : 0.0; }
327 void errstrain(const double &d) { if (!hasEnergies()) return; energies_->errstrain_ = d; }
328 double eslip() const { return hasEnergies() ? energies_->eslip_ : 0.0; }
329 void eslip(const double &d) { if (!hasEnergies()) return; energies_->eslip_ = d; }
330 double errslip() const { return hasEnergies() ? energies_->errslip_ : 0.0; }
331 void errslip(const double &d) { if (!hasEnergies()) return; energies_->errslip_ = d; }
332 double edashpot() const { return hasEnergies() ? energies_->edashpot_ : 0.0; }
333 void edashpot(const double &d) { if (!hasEnergies()) return; energies_->edashpot_ = d; }
334
335 uint32 inheritanceField() const { return inheritanceField_; }
336 void inheritanceField(uint32 i) { inheritanceField_ = i; }
337
338 const DVect2 & effectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_; }
339 void effectiveTranslationalStiffness(const DVect2 &v) { effectiveTranslationalStiffness_ = v; }
340 const DAVect & effectiveRotationalStiffness() const { return effectiveRotationalStiffness_; }
341 void effectiveRotationalStiffness(const DAVect &v) { effectiveRotationalStiffness_ = v; }
342
343 // Rolling resistance methods
344 const double & res_fric() const { return res_fric_; }
345 void res_fric(const double &d) { res_fric_ = d; }
346 const DAVect & res_M() const { return res_M_; }
347 void res_M(const DAVect &f) { res_M_ = f; }
348 bool res_S() const { return res_S_; }
349 void res_S(bool b) { res_S_ = b; }
350 const double & kr() const { return kr_; }
351 void kr(const double &d) { kr_ = d; }
352 const double & fr() const { return fr_; }
353 void fr(const double &d) { fr_ = d; }
354 const double & rbar_square() const { return rbar_square_; }
355 void rbar_square(const double &d) { rbar_square_ = d; }
356 //
357
358 private:
359 // Index - used internally by PFC. Should be set to -1 in the cpp file.
360 static int index_;
361
362 // Structure to store the energies.
363 struct Energies {
364 Energies() : estrain_(0.0), errstrain_(0.0), eslip_(0.0), errslip_(0.0), edashpot_(0.0) {}
365 double estrain_; // elastic energy stored in linear group
366 double errstrain_; // elastic energy stored in rolling resistance group
367 double eslip_; // work dissipated by friction
368 double errslip_; // work dissipated by rolling resistance friction
369 double edashpot_; // work dissipated by dashpots
370 };
371
372 // Structure to store dashpot quantities.
373 struct dpProps {
374 dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
375 double dp_nratio_; // normal viscous critical damping ratio
376 double dp_sratio_; // shear viscous critical damping ratio
377 int dp_mode_; // for viscous mode (0-1) 0 = no limit or cutoff; 1 = shear cut-off when sliding
378 DVect dp_F_; // Force in the dashpots
379 };
380
381 bool updateEndStiffCoef(const IContactMechanical *con);
382 bool updateEndFric(const IContactMechanical *con);
383 bool updateEndResFric(const IContactMechanical *con);
384 bool updateEndSurfAdh(const IContactMechanical *con);
385 void updateStiffCoef(const IContactMechanical *con);
386 void updateEffectiveStiffness(ContactModelMechanicalState *state);
387 void setDampCoefficients(const double &mass, double *kn_tangent, double *ks_tangent, double *vcn, double *vcs);
388
389 bool SetJKRTable(ContactModelMechanicalState *state);
390 DVect2 InterpolateJKRTable_Method1(double overlap_ratio);
391 DVect2 InterpolateJKRTable_Method2(double overlap_ratio, double trans_n);
392 DVect2 JKR_Analytical(double overlap);
393
394 static const uint32 shearMask = 0x00000002; // Base 1!
395 static const uint32 poissMask = 0x00000004;
396 static const uint32 fricMask = 0x00000008;
397 static const uint32 resFricMask = 0x00004000;
398 static const uint32 surfAdhMask = 0x00004000;
399
400 // Contact model inheritance fields.
401 uint32 inheritanceField_ = shearMask | poissMask | fricMask | resFricMask | surfAdhMask;
402
403 // Effective translational stiffness.
404 DVect2 effectiveTranslationalStiffness_ = DVect2(0);
405 DAVect effectiveRotationalStiffness_ = DAVect(0); // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
406
407 // JKR model properties
408 double shear_ = 0.0; // Shear modulus
409 double poiss_ = 0.0; // Poisson's ratio
410 double fric_ = 0.0; // Coulomb friction coefficient
411 DVect jkr_F_ = DVect(0); // Force carried in the eepa model
412 bool jkr_S_ = false; // The current slip state
413 double rgap_ = 0.0; // Reference gap
414 dpProps * dpProps_ = nullptr; // The viscous properties
415 double kn_ = 0.0; // normal stiffness
416 double ks_ = 0.0; // shear stiffness
417 double ks_fac_ = 1.0; // shear stiffness scaling factor
418 double surf_adh_ = 0.0; // surface adhesion/cohesion energy (energy/area = J/m^2)
419 int active_mode_ = 1; // mode=0: no negtaive overlap is allowed - contact activated and deactivated when physical contact is made (Simplified SJKR-A model)
420 // mode=1: (default) negative overlap is allowed - contact activated on physical contact, but negative overlap (up to tear-off distance) allowed when unloading (full JKR model)
421 double distance_active_ = 0.0;// active distance - initially set to zero to establish contact, then after first time active it is set to the tear_off_ distance to allow for a seperation distance while under adhesion force
422 double a0_ = 0.0; // contact patch radius where the JKR force is equal to zero
423 double pull_off_f_ = 0.0; // the pull-off force which is the maximum tensile (adhesive) force
424 double tear_off_d_ = 0.0; // tear-off distance where the the adhesion snaps to zero under te
425 double r_hertz_ = 0.0; // effective contact radius: r_hertz_ = R1xR2/(R1+R2)
426 double c1_ = 0.0; // constant used in the analytical solution to the JKR contact patch radius
427
428 // rolling resistance properties
429 double res_fric_ = 0.0; // rolling friction coefficient
430 DAVect res_M_ = DAVect(0); // moment (bending only)
431 bool res_S_ = false; // The current rolling resistance slip state
432 double kr_ = 0.0; // bending rotational stiffness (read-only, calculated internaly)
433 double fr_ = 0.0; // rolling friction coefficient (rbar*res_fric_) (calculated internaly, not a property)
434 double rbar_square_ = 0.0; // used to update the rolling stiffness
435
436 Energies * energies_ = nullptr; // The energies
437 };
438}
439// EoF
contactmodeljkr.cpp
1// contactmodeljkr.cpp
2#include "contactmodeljkr.h"
3
4#include "module/interface/icontactmechanical.h"
5#include "module/interface/icontact.h"
6#include "module/interface/ipiece.h"
7#include "module/interface/ifishcalllist.h"
8
9#include "../version.txt"
10
11#include "utility/src/tptr.h"
12
13#include "kernel/interface/iprogram.h"
14#include "fish/src/parameter.h"
15
16#ifdef JKR_LIB
17int __stdcall DllMain(void *, unsigned, void *) {
18 return 1;
19}
20
21extern "C" EXPORT_TAG const char *getName() {
22#if DIM==3
23 return "contactmodelmechanical3dJKR";
24#else
25 return "contactmodelmechanical2dJKR";
26#endif
27}
28
29extern "C" EXPORT_TAG unsigned getMajorVersion() {
30 return MAJOR_VERSION;
31}
32
33extern "C" EXPORT_TAG unsigned getMinorVersion() {
34 return MINOR_VERSION;
35}
36
37extern "C" EXPORT_TAG void *createInstance() {
38 cmodelsxd::ContactModelJKR *m = NEW cmodelsxd::ContactModelJKR();
39 return (void *)m;
40}
41#endif
42
43namespace cmodelsxd {
44 using namespace itasca;
45
46 int ContactModelJKR::index_ = -1;
47 uint32 ContactModelJKR::getMinorVersion() const { return MINOR_VERSION; }
48
49 ContactModelJKR::ContactModelJKR() {
50 }
51
52 ContactModelJKR::ContactModelJKR(const ContactModelJKR &m) noexcept {
53 this->copy(&m);
54 }
55
56 const ContactModelJKR & ContactModelJKR::operator=(const ContactModelJKR &m) {
57 this->copy(&m);
58 return *this;
59 }
60
61 void ContactModelJKR::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
62 s->addToStorage<ContactModelJKR>(*this,ww);
63 }
64
65 ContactModelJKR::~ContactModelJKR() {
66 // Make sure to clean up after yourself!
67 if (dpProps_)
68 delete dpProps_;
69 if (energies_)
70 delete energies_;
71 }
72
73 void ContactModelJKR::archive(ArchiveStream &stream) {
74 // The stream allows one to archive the values of the contact model
75 // so that it can be saved and restored. The minor version can be
76 // used here to allow for incremental changes to the contact model too.
77 stream & shear_;
78 stream & poiss_;
79 stream & fric_;
80 stream & jkr_F_;
81 stream & jkr_S_;
82 stream & rgap_;
83 stream & res_fric_;
84 stream & res_M_;
85 stream & res_S_;
86 stream & kr_;
87 stream & fr_;
88 stream & rbar_square_;
89 stream & kn_;
90 stream & ks_;
91 stream & ks_fac_;
92 stream & surf_adh_;
93 stream & active_mode_;
94 stream & distance_active_;
95 stream & a0_;
96 stream & pull_off_f_;
97 stream & tear_off_d_;
98 stream & r_hertz_;
99 stream & c1_;
100
101 if (stream.getArchiveState() == ArchiveStream::Save) {
102 bool b = false;
103 if (dpProps_) {
104 b = true;
105 stream & b;
106 stream & dpProps_->dp_nratio_;
107 stream & dpProps_->dp_sratio_;
108 stream & dpProps_->dp_mode_;
109 stream & dpProps_->dp_F_;
110 }
111 else
112 stream & b;
113
114 b = false;
115 if (energies_) {
116 b = true;
117 stream & b;
118 stream & energies_->estrain_;
119 stream & energies_->errstrain_;
120 stream & energies_->eslip_;
121 stream & energies_->errslip_;
122 stream & energies_->edashpot_;
123 }
124 else
125 stream & b;
126 }
127 else {
128 bool b(false);
129 stream & b;
130 if (b) {
131 if (!dpProps_)
132 dpProps_ = NEW dpProps();
133 stream & dpProps_->dp_nratio_;
134 stream & dpProps_->dp_sratio_;
135 stream & dpProps_->dp_mode_;
136 stream & dpProps_->dp_F_;
137 }
138 stream & b;
139 if (b) {
140 if (!energies_)
141 energies_ = NEW Energies();
142 stream & energies_->estrain_;
143 stream & energies_->errstrain_;
144 stream & energies_->eslip_;
145 stream & energies_->errslip_;
146 stream & energies_->edashpot_;
147 }
148 }
149
150 stream & inheritanceField_;
151 stream & effectiveTranslationalStiffness_;
152 stream & effectiveRotationalStiffness_;
153 }
154
155 void ContactModelJKR::copy(const ContactModel *cm) {
156 // Copy all of the contact model properties. Used in the CMAT
157 // when a new contact is created.
158 ContactModelMechanical::copy(cm);
159 const ContactModelJKR *in = dynamic_cast<const ContactModelJKR*>(cm);
160 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
161 shear(in->shear());
162 poiss(in->poiss());
163 fric(in->fric());
164 jkr_F(in->jkr_F());
165 jkr_S(in->jkr_S());
166 rgap(in->rgap());
167 res_fric(in->res_fric());
168 res_M(in->res_M());
169 res_S(in->res_S());
170 kr(in->kr());
171 fr(in->fr());
172 rbar_square(in->rbar_square());
173 kn(in->kn());
174 ks(in->ks());
175 ks_fac(in->ks_fac());
176 surf_adh(in->surf_adh());
177 active_mode(in->active_mode());
178 distance_active(in->distance_active());
179 a0(in->a0());
180 pull_off(in->pull_off());
181 tear_off(in->tear_off());
182 r_hertz(in->r_hertz());
183 c1(in->c1());
184
185 if (in->hasDamping()) {
186 if (!dpProps_)
187 dpProps_ = NEW dpProps();
188 dp_nratio(in->dp_nratio());
189 dp_sratio(in->dp_sratio());
190 dp_mode(in->dp_mode());
191 dp_F(in->dp_F());
192 }
193 if (in->hasEnergies()) {
194 if (!energies_)
195 energies_ = NEW Energies();
196 estrain(in->estrain());
197 errstrain(in->errstrain());
198 eslip(in->eslip());
199 errslip(in->errslip());
200 edashpot(in->edashpot());
201 }
202 inheritanceField(in->inheritanceField());
203 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
204 effectiveRotationalStiffness(in->effectiveRotationalStiffness());
205 }
206
207 base::Property ContactModelJKR::getProperty(uint32 i, const IContact *) const {
208 // Return the property. The IContact pointer is provided so that
209 // more complicated properties, depending on contact characteristics,
210 // can be calcualted. The IContact pointer may be a nullptr!
211 base::Property var;
212 switch (i) {
213 case kwShear: return shear_;
214 case kwPoiss: return poiss_;
215 case kwFric: return fric_;
216 case kwjkrF: var.setValue(jkr_F_); return var;
217 case kwjkrS: return jkr_S_;
218 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
219 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
220 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
221 case kwDpF: {
222 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
223 return var;
224 }
225 case kwResFric: return res_fric_;
226 case kwResMoment: var.setValue(res_M_); return var;
227 case kwResS: return res_S_;
228
229 case kwKsFac: return ks_fac_;
230 case kwSurfAdh: return surf_adh_;
231 case kwActiveMode: return active_mode_;
232 case kwA0: return a0_;
233 case kwPullOff: return pull_off_f_;
234 case kwTearOff: return tear_off_d_;
235 }
236 assert(0);
237 return base::Property();
238 }
239
240 bool ContactModelJKR::getPropertyGlobal(uint32 i) const {
241 // Returns whether or not a property is held in the global axis system (TRUE)
242 // or the local system (FALSE). Used by the plotting logic.
243 switch (i) {
244 case kwjkrF:
245 case kwDpF:
246 case kwResMoment:
247 return false;
248 }
249 return true;
250 }
251
252 bool ContactModelJKR::setProperty(uint32 i, const base::Property &v, IContact *) {
253 // Set a contact model property. Return value indicates that the timestep
254 // should be recalculated.
255 dpProps dp;
256 switch (i) {
257 case kwShear: {
258 if (!v.canConvert<double>())
259 throw Exception("shear must be a double.");
260 double val(v.toDouble());
261 if (val <= 0.0)
262 throw Exception("zero or negative shear not allowed in the 'jkr' model.");
263 shear_ = val;
264 return true;
265 }
266 case kwPoiss: {
267 if (!v.canConvert<double>())
268 throw Exception("poiss must be a double.");
269 double val(v.toDouble());
270 if (val < 0.0 || val > 0.5)
271 throw Exception("poiss must be in the range [0, 0.5].");
272 poiss_ = val;
273 return true;
274 }
275 case kwFric: {
276 if (!v.canConvert<double>())
277 throw Exception("fric must be a double.");
278 double val(v.toDouble());
279 if (val < 0.0)
280 throw Exception("negative fric not allowed.");
281 fric_ = val;
282 return false;
283 }
284 case kwjkrF: {
285 if (!v.canConvert<DVect>())
286 throw Exception("jkr_force must be a vector.");
287 DVect val(v.value<DVect>());
288 jkr_F_ = val;
289 return false;
290 }
291 case kwDpNRatio: {
292 if (!v.canConvert<double>())
293 throw Exception("dp_nratio must be a double.");
294 double val(v.toDouble());
295 if (val < 0.0)
296 throw Exception("negative dp_nratio not allowed.");
297 if (val == 0.0 && !dpProps_)
298 return false;
299 if (!dpProps_)
300 dpProps_ = NEW dpProps();
301 dpProps_->dp_nratio_ = val;
302 return true;
303 }
304 case kwDpSRatio: {
305 if (!v.canConvert<double>())
306 throw Exception("dp_sratio must be a double.");
307 double val(v.toDouble());
308 if (val < 0.0)
309 throw Exception("negative dp_sratio not allowed.");
310 if (val == 0.0 && !dpProps_)
311 return false;
312 if (!dpProps_)
313 dpProps_ = NEW dpProps();
314 dpProps_->dp_sratio_ = val;
315 return true;
316 }
317 case kwDpMode: {
318 if (!v.canConvert<int64>())
319 throw Exception("the viscous mode dp_mode must be 0 or 1");
320 int val(v.toInt());
321 if (val == 0 && !dpProps_)
322 return false;
323 if (val < 0 || val > 1)
324 throw Exception("the viscous mode dp_mode must be 0 or 1.");
325 if (!dpProps_)
326 dpProps_ = NEW dpProps();
327 dpProps_->dp_mode_ = val;
328 return false;
329 }
330 case kwDpF: {
331 if (!v.canConvert<DVect>())
332 throw Exception("dp_force must be a vector.");
333 DVect val(v.value<DVect>());
334 if (val.fsum() == 0.0 && !dpProps_)
335 return false;
336 if (!dpProps_)
337 dpProps_ = NEW dpProps();
338 dpProps_->dp_F_ = val;
339 return false;
340 }
341 case kwResFric: {
342 if (!v.canConvert<double>())
343 throw Exception("res_fric must be a double.");
344 double val(v.toDouble());
345 if (val < 0.0)
346 throw Exception("negative res_fric not allowed.");
347 res_fric_ = val;
348 return true;
349 }
350 case kwResMoment: {
351 DAVect val(0.0);
352#ifdef TWOD
353 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
354 throw Exception("res_moment must be an angular vector.");
355 if (v.canConvert<DAVect>())
356 val = DAVect(v.value<DAVect>());
357 else
358 val = DAVect(v.toDouble());
359#else
360 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
361 throw Exception("res_moment must be an angular vector.");
362 if (v.canConvert<DAVect>())
363 val = DAVect(v.value<DAVect>());
364 else
365 val = DAVect(v.value<DVect>());
366#endif
367 res_M_ = val;
368 return false;
369 }
370 case kwKsFac: {
371 if (!v.canConvert<double>())
372 throw Exception("shear stiffness factor must be a double.");
373 double val(v.toDouble());
374 if (val <= 0.0)
375 throw Exception("negative or zero ks_fac not allowed.");
376 ks_fac_ = val;
377 return true;
378 }
379 case kwSurfAdh: {
380 if (!v.canConvert<double>())
381 throw Exception("surface adhesion must be a double.");
382 double val(v.toDouble());
383 if (val <= 0.0)
384 throw Exception("negative or zero surf_adh not allowed.");
385 surf_adh_ = val;
386 return true;
387 }
388 case kwActiveMode: {
389 if (!v.canConvert<int64>())
390 throw Exception("the active mode must be an integer: 0 or 1.");
391 int val(v.toInt());
392 if (val < 0 || val > 1)
393 throw Exception("active_mode must be 0 or 1.");
394 active_mode_ = val;
395 return true;
396 }
397 }//switch
398 return false;
399 }
400
401 bool ContactModelJKR::getPropertyReadOnly(uint32 i) const {
402 // Returns TRUE if a property is read only or FALSE otherwise.
403 switch (i) {
404 case kwDpF:
405 case kwjkrS:
406 case kwResS:
407 case kwA0:
408 case kwPullOff:
409 case kwTearOff:
410 return true;
411 default:
412 break;
413 }
414 return false;
415 }
416
417 bool ContactModelJKR::supportsInheritance(uint32 i) const {
418 // Returns TRUE if a property supports inheritance or FALSE otherwise.
419 switch (i) {
420 case kwShear:
421 case kwPoiss:
422 case kwFric:
423 case kwResFric:
424 case kwSurfAdh:
425 return true;
426 default:
427 break;
428 }
429 return false;
430 }
431
432 double ContactModelJKR::getEnergy(uint32 i) const {
433 // Return an energy value.
434 double ret(0.0);
435 if (!energies_)
436 return ret;
437 switch (i) {
438 case kwEStrain: return energies_->estrain_;
439 case kwERRStrain: return energies_->errstrain_;
440 case kwESlip: return energies_->eslip_;
441 case kwERRSlip: return energies_->errslip_;
442 case kwEDashpot: return energies_->edashpot_;
443 }
444 assert(0);
445 return ret;
446 }
447
448 bool ContactModelJKR::getEnergyAccumulate(uint32 i) const {
449 // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
450 switch (i) {
451 case kwEStrain: return true;
452 case kwERRStrain: return false;
453 case kwESlip: return true;
454 case kwERRSlip: return true;
455 case kwEDashpot: return true;
456 }
457 assert(0);
458 return false;
459 }
460
461 void ContactModelJKR::setEnergy(uint32 i, const double &d) {
462 // Set an energy value.
463 if (!energies_) return;
464 switch (i) {
465 case kwEStrain: energies_->estrain_ = d; return;
466 case kwERRStrain: energies_->errstrain_ = d; return;
467 case kwESlip: energies_->eslip_ = d; return;
468 case kwERRSlip: energies_->errslip_ = d; return;
469 case kwEDashpot: energies_->edashpot_ = d; return;
470 }
471 assert(0);
472 return;
473 }
474
475 bool ContactModelJKR::validate(ContactModelMechanicalState *state, const double &) {
476 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
477 // (1) the contact is created, (2) a property of the contact that returns a true via
478 // the setProperty method has been modified and (3) when a set of cycles is executed
479 // via the {cycle N} command.
480 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
481 assert(state);
482 const IContactMechanical *c = state->getMechanicalContact();
483 assert(c);
484 //
485 if (state->trackEnergy_)
486 activateEnergy();
487 //
488 if ((inheritanceField_ & shearMask) || (inheritanceField_ & poissMask))
489 updateEndStiffCoef(c);
490 if (inheritanceField_ & fricMask)
491 updateEndFric(c);
492 if (inheritanceField_ & resFricMask)
493 updateEndResFric(c);
494 if (inheritanceField_ & surfAdhMask)
495 updateEndSurfAdh(c);
496
497 updateStiffCoef(c); // calculate the stiffness values based on the material properties - specified directly or via inheritance
498 updateEffectiveStiffness(state); // effective stiffness for translation and rotation used in the time step estimation
499 //
500 return checkActivity(state->gap_);
501 }
502
503 void ContactModelJKR::updateStiffCoef(const IContactMechanical *con) {
504 //
505 double c12 = con->getEnd1Curvature().y();
506 double c22 = con->getEnd2Curvature().y();
507 r_hertz_ = c12 + c22;
508 if (r_hertz_ == 0.0)
509 throw Exception("jkr contact model undefined for 2 non-curved surfaces");
510 //
511 r_hertz_ = 1.0 / r_hertz_; // r_hertz = R1*R2/(R1 + R2)
512 double young_eff = shear_ / (1.0 - poiss_); // effective contact Young's modulus assuming the two contact pieces are identical in properties
513 double shear_eff = shear_ / (4.0 - 2.0*poiss_); // effective shear modulus for contact assuming the two contact pieces are identical in properties
514 kn_ = 2.0 * young_eff; // normal stiffness - needs to be multiplied by contact patch radius to obtain the tangent normal stiffness - (Marshall, 2009)
515 ks_ = ks_fac_ * 8.0*shear_eff; // shear stiffness - needs to be multiplied by contact patch radius to obtain the tangent shear stiffness - (Marshall, 2009)
516 //
517 a0_ = 9.0*dPi*surf_adh_*r_hertz_*r_hertz_ / young_eff; // contact patch radius where the JKR force is equal to zero
518 a0_ = pow(a0_, 1.0 / 3.0);
519 //
520 pull_off_f_ = 3.0*dPi*surf_adh_*r_hertz_; // pull-off force, the maximum tensile/adhesion force in JKR model (not the force at tear-off when the contact snaps)
521 tear_off_d_ = a0_* a0_ / (2.0*pow(6.0, (1.0 / 3.0))*r_hertz_); // tear-off distance where the the adhesion snaps to zero under tension
522 //
523 c1_ = -4.0*dPi*surf_adh_ * r_hertz_*r_hertz_ / young_eff; // constant used in the analytical solution for the patch radius
524 //
525 // rolling resistance
526 kr_ = 0.0;
527 fr_ = 0.0;
528 if (res_fric_ > 0.0) {
529 rbar_square_ = r_hertz_ * r_hertz_; // store this value with contact - used again when 'ks_tangent' is updated during loading-unloading to calculate the rolling stiffness
530 fr_ = res_fric_ * r_hertz_; // store this value with contact - used in force-displacement calculation
531 kr_ = ks_ * rbar_square_; // based on 'ks_' for now, this is updated and based on the shear tangent stiffness during loading-unloading
532 }
533 }
534
535 bool ContactModelJKR::endPropertyUpdated(const string &name, const IContactMechanical *c) {
536 // The endPropertyUpdated method is called whenever a surface property (with a name
537 // that matches an inheritable contact model property name) of one of the contacting
538 // pieces is modified. This allows the contact model to update its associated
539 // properties. The return value denotes whether or not the update has affected
540 // the time step computation (by having modified the translational or rotational
541 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
542 assert(c);
543 StringList availableProperties = split(replace(simplified(getProperties())," ", ""),",");
544 auto idx = findRegex(availableProperties,name);
545 if (idx==string::npos) return false;
546 idx += 1;
547 bool ret = false;
548
549 switch (idx) {
550 case kwShear: { //Shear
551 if (inheritanceField_ & shearMask)
552 ret = true; // return 'true' to ensure 'validate()' is called to affect the change in Shear Modulus through 'updateEndStiffCoef()'
553 break;
554 }
555 case kwPoiss: { //Poisson's
556 if (inheritanceField_ & poissMask)
557 ret =true; // return 'true' to ensure 'validate()' is called to affect the change in Poisson's ratio through 'updateEndStiffCoef()'
558 break;
559 }
560 case kwFric: { //fric
561 if (inheritanceField_ & fricMask)
562 updateEndFric(c); // friction does not influence any of the other parameters or time step size, so directly update the contact model friction
563 break;
564 }
565 case kwResFric: { //rr_fric
566 if (inheritanceField_ & resFricMask)
567 ret = true; // return 'true' to ensure 'validate()' is called to affect the change in rolling friction through 'updateEndResFric()' and 'fr_'
568 break;
569 }
570 case kwSurfAdh: { //surf_adh
571 if (inheritanceField_ & surfAdhMask)
572 ret = true; // return 'true' to ensure 'validate()' is called to affect the change in surface adhesion through 'updateEndSurfAdh()'
573 break;
574 }
575 //
576 }
577 return ret;
578 }
579
580 static const string gstr("jkr_shear");
581 static const string nustr("jkr_poiss");
582 bool ContactModelJKR::updateEndStiffCoef(const IContactMechanical *con) {
583 assert(con);
584 double g1 = shear_;
585 double g2 = shear_;
586 double nu1 = poiss_;
587 double nu2 = poiss_;
588 base::Property vg1 = con->getEnd1()->getProperty(gstr);
589 base::Property vg2 = con->getEnd2()->getProperty(gstr);
590 base::Property vnu1 = con->getEnd1()->getProperty(nustr);
591 base::Property vnu2 = con->getEnd2()->getProperty(nustr);
592 if (vg1.isValid() && vg2.isValid()) {
593 g1 = vg1.toDouble();
594 g2 = vg2.toDouble();
595 if (g1 <= 0.0 || g2 <= 0.0)
596 throw Exception("Negative or zero shear modulus not allowed in jkr contact model");
597 }
598 if (vnu1.isValid() && vnu2.isValid()) {
599 nu1 = vnu1.toDouble();
600 nu2 = vnu2.toDouble();
601 if (nu1 < 0.0 || nu1 > 0.5 || nu2 < 0.0 || nu2 > 0.5)
602 throw Exception("Poisson ratio should be in range [0,0.5] in jkr contact model");
603 }
604 if (g1*g2 == 0.0) return false;
605 double es = 1.0 / ((1.0 - nu1) / (2.0*g1) + (1.0 - nu2) / (2.0*g2));
606 double gs = 1.0 / ((2.0 - nu1) / g1 + (2.0 - nu2) / g2);
607 poiss_ = (4.0*gs - es) / (2.0*gs - es);
608 shear_ = 2.0*gs*(2 - poiss_);
609 if (shear_ <= 0.0)
610 throw Exception("Negative or zero shear modulus not allowed in jkr contact model");
611 if (poiss_ < 0.0 || poiss_ > 0.5)
612 throw Exception("Poisson ratio should be in range [0,0.5] in jkr contact model");
613 return true;
614 }
615
616 static const string fricstr("fric");
617 bool ContactModelJKR::updateEndFric(const IContactMechanical *con) {
618 assert(con);
619 base::Property v1 = con->getEnd1()->getProperty(fricstr);
620 base::Property v2 = con->getEnd2()->getProperty(fricstr);
621 if (!v1.isValid() || !v2.isValid())
622 return false;
623 double fric1 = std::max(0.0, v1.toDouble());
624 double fric2 = std::max(0.0, v2.toDouble());
625 double val = fric_;
626 fric_ = std::min(fric1, fric2);
627 if (fric_ < 0.0)
628 throw Exception("Negative friction value not allowed in jkr contact model");
629 return ((fric_ != val));
630 }
631
632 static const string rfricstr("rr_fric");
633 bool ContactModelJKR::updateEndResFric(const IContactMechanical *con) {
634 assert(con);
635 base::Property v1 = con->getEnd1()->getProperty(rfricstr);
636 base::Property v2 = con->getEnd2()->getProperty(rfricstr);
637 if (!v1.isValid() || !v2.isValid())
638 return false;
639 double rfric1 = std::max(0.0, v1.toDouble());
640 double rfric2 = std::max(0.0, v2.toDouble());
641 double val = res_fric_;
642 res_fric_ = std::min(rfric1, rfric2);
643 if (res_fric_ < 0.0)
644 throw Exception("Negative rolling friction value not allowed in jkr contact model");
645 return ((res_fric_ != val));
646 }
647
648 static const string surfadhstr("surf_adh");
649 bool ContactModelJKR::updateEndSurfAdh(const IContactMechanical *con) {
650 assert(con);
651 base::Property v1 = con->getEnd1()->getProperty(surfadhstr);
652 base::Property v2 = con->getEnd2()->getProperty(surfadhstr);
653 if (!v1.isValid() || !v2.isValid())
654 return false;
655 double surfadh1 = std::max(0.0, v1.toDouble());
656 double surfadh2 = std::max(0.0, v2.toDouble());
657 if (surfadh1 <= 0.0 || surfadh2 <= 0.0)
658 throw Exception("Negative or zero surface adhesion not allowed in 'jkr' contact model");
659 double val = surf_adh_;
660 surf_adh_ = 0.5*(surfadh1 + surfadh2); // new surface adhesion enery
661 return ((surf_adh_ != val));
662 }
663
664 void ContactModelJKR::updateEffectiveStiffness(ContactModelMechanicalState *state) {
665
666 double overlap = rgap_ - state->gap_;
667 if (overlap <= -distance_active_) { // if contact was just created but not yet activated, distance_active_ == 0.
668 // The assumption below is based on the Hill contact model in PFC
669 // 0.01% of Diameter = 0.0001*D
670 // For monodisperse system, r_hertz_ = 0.5*R where R is the particle radius
671 // r_hertz_ = 0.5*(0.5*D) = 0.25*D where D is the particle diameter
672 // D = 4*r_hertz_
673 // 0.01%*D = 0.0001*D = 0.0001*4*r_hertz_ = 0.0004*r_hertz_
674 overlap = 0.0004*r_hertz_;
675 }
676 //
677 DVect2 force_radius = JKR_Analytical(overlap); // analytical patch radius and force for this (assumed) overlap
678 double patch_radius = force_radius.dof(1); // contact patch radius
679 // Assume Hertz tangent stiffness in the normal direction - the tangent stiffness for JKR is different, but no closed-form exists - good enough of an assumption
680 double kn_tangent = kn_ * patch_radius; // tangent stiffness in normal direction
681 double ks_tangent = ks_ * patch_radius; // tangent stiffness in the shear direction - tis is not an assumption, but the correct solution for JKR
682 effectiveTranslationalStiffness_ = DVect2(kn_tangent, ks_tangent); // set using the tangent stiffness
683 if (res_fric_ > 0.0) { // rolling
684 kr_ = ks_tangent * rbar_square_; // based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
685 effectiveRotationalStiffness_ = DAVect(kr_); // effective rotational stiffness (bending only)
686#if DIM==3
687 effectiveRotationalStiffness_.rx() = 0.0;
688#endif
689 }
690 // correction if viscous damping active
691 if (dpProps_) {
692 DVect2 correct(1.0);
693 if (dpProps_->dp_nratio_)
694 correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
695 if (dpProps_->dp_sratio_)
696 correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
697 effectiveTranslationalStiffness_ /= (correct*correct);
698 }
699 }
700 //
701 // --------------------------------------------------------------------------------------------------------------------------------------------------------
702 //
703 DVect2 ContactModelJKR::JKR_Analytical(double overlap) { // analytical solution according to Parteli et al. 2014
704 DVect2 ret(0.0);
705 // reduces to Hertz model
706 if (surf_adh_ == 0.0) {
707 double patch_radius = std::sqrt(r_hertz_*overlap);
708 ret.rdof(0) = (2.0 / 3.0)*kn_*patch_radius*patch_radius*patch_radius / r_hertz_;
709 ret.rdof(1) = patch_radius;
710 return ret;
711 }
712 // otherwise jkr is active ...
713 double rad_over = r_hertz_ * overlap;
714 double c0 = rad_over * rad_over;
715 double c2 = -2.0*rad_over;
716 double P = -1.0 / 12.0*c2*c2 - c0;
717 double Q = -1.0 / 108.0*c2*c2*c2 + c0 * c2 / 3.0 - c1_*c1_ / 8.0;
718 double U = pow((-0.5*Q + sqrt(0.25*Q*Q + P*P*P / 27.0)), (1.0 / 3.0));
719 double s = 0.0;
720 if (P == 0.0) {
721 s = -5.0*c2 / 6.0 - pow(Q, (1.0 / 3.0));
722 }
723 else {
724 s = -5.0*c2 / 6.0 + U - P / (3.0*U);
725 }
726 double w = sqrt(c2 + 2.0*s);
727 double lambda = c1_ / (2.0*w);
728 double patch_radius = 0.5*(w + sqrt(w*w - 4.0*(c2 + s + lambda)));
729 // reduces to Hertz model due to very small values for surface adhesion
730 if (std::isnan(patch_radius)) {
731 patch_radius = std::sqrt(r_hertz_*overlap);
732 ret.rdof(0) = (2.0 / 3.0)*kn_*patch_radius*patch_radius*patch_radius / r_hertz_;
733 ret.rdof(1) = patch_radius;
734 return ret;
735 }
736 // force = 4*E*(a^3)/(3*R) - sqrt(16*pi*surf_adh*E*a^3);
737 // with E = kn/2 - where E is the effective contact Young's modulus - see updateStiffCoef(const IContactMechanical *con)
738 // the above reduces to the following for the JKR force ...
739 double patch_radius3 = patch_radius* patch_radius*patch_radius;
740 ret.rdof(0) = (2.0 / 3.0)*kn_*patch_radius3 / r_hertz_ - sqrt(8.0*dPi*surf_adh_*kn_*patch_radius3); // set the total force
741 ret.rdof(1) = patch_radius; // set the patch contact radius
742 return ret;
743 }
744 //
745 // --------------------------------------------------------------------------------------------------------------------------------------------------------
746 //
747 bool ContactModelJKR::forceDisplacementLaw(ContactModelMechanicalState *state, const double ×tep) {
748 assert(state);
749 if (shear_ <= 0.0) {
750 throw Exception("'shear' must be specified using a value larger than zero in the 'jkr' model.");
751 }
752 if (surf_adh_ <= 0.0) {
753 throw Exception("'surf_adh' must be specified using a value larger than zero in the 'jkr' model. For zero surface adhesion, use the Hertz model!");
754 }
755 //
756 // Current overlap
757 double overlap = rgap_ - state->gap_; // state->gap is negative in sign when the particles physically overlap. If rgap == 0, then the overlap is positive when the two particles physically overlap.
758 // Relative translational increment
759 DVect trans = state->relativeTranslationalIncrement_;
760 //
761 if (active_mode_ == 1 && distance_active_ == 0.0) { // only when in 'full jkr' mode = 1, for Simplified JKR-A (SJKR-A), negative overlap is not allowed
762 distance_active_ = tear_off_d_; // contact has become active - set the activity distance to rgap_ + distance_active = rgap_ + tear_off to allow for negative overlap under adhesion force
763 }
764 //
765 // The contact was just activated from an inactive state
766 // Trigger the FISH callback if one is hooked up to the contact_activated event.
767 if (cmEvents_[fActivated] >= 0) {
768 // An FArray of base::Property is returned and these will be passed
769 // to the FISH function as an array of FISH symbols as the second
770 // argument to the FISH callback function.
771 auto c = state->getContact();
772 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
773 IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
774 fi->setCMFishCallArguments(c, arg, cmEvents_[fActivated]);
775
776 }
777 //
778 // Angular dispacement increment.
779 DAVect ang = state->relativeAngularIncrement_;
780 //
781 DVect jkr_F_old = jkr_F_; // store for energy computations only
782 //
783 // NORMAL FORCE ==============================================================================================================
784 //
785 DVect2 force_radius = JKR_Analytical(overlap); // analytical force and patch radius
786 jkr_F_.rdof(0) = force_radius.dof(0); // current contact force in the normal direction
787 double patch_radius = force_radius.dof(1); // current contact patch radius
788 //
789 double kn_tangent = kn_ * patch_radius;
790 if (trans.dof(0))
791 kn_tangent = std::abs((jkr_F_old.dof(0) - jkr_F_.dof(0)) / trans.dof(0));// tangent stiffness in normal direction - used for timestep calculations and damping
792 //
793 // SHEAR FORCE ==============================================================================================================
794 //
795 double ks_tangent = ks_* patch_radius; // tangent stiffness in the shear direction
796 DVect sforce(0.0);
797 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
798 // Loop over the shear components (note: the 0 component is the normal component) and calculate the shear force.
799 for (int i = 1; i < dim; ++i)
800 sforce.rdof(i) = jkr_F_.dof(i) - trans.dof(i) * ks_tangent; // shear force update
801 //
802 // The canFail flag corresponds to whether or not the contact can undergo non-linear
803 // force-displacement response. If the SOLVE ELASTIC command is given then the
804 // canFail state is set to FALSE. Otherwise it is always TRUE.
805 if (state->canFail_) {
806 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
807 double crit = fric_ * std::abs(jkr_F_.rdof(0) + 2.0*pull_off_f_);
808 // The is the magnitude of the shear force.
809 double sfmag = sforce.mag();
810 // Sliding occurs when the magnitude of the shear force is greater than the critical value.
811 if (sfmag > crit) {
812 // Lower the shear force to the critical value for sliding.
813 double rat = crit / sfmag;
814 sforce *= rat;
815 // Handle the slip_change event if one has been hooked up. Sliding has commenced.
816 if (!jkr_S_ && cmEvents_[fSlipChange] >= 0) {
817 auto c = state->getContact();
818 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
819 fish::Parameter((int64)0) };
820 IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
821 fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
822 }
823 jkr_S_ = true;
824 }
825 else {
826 // Handle the slip_change event if one has been hooked up and
827 // the contact was previously sliding. Sliding has ceased.
828 if (jkr_S_) {
829 if (cmEvents_[fSlipChange] >= 0) {
830 auto c = state->getContact();
831 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
832 fish::Parameter((int64)1) };
833 IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
834 fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
835 }
836 jkr_S_ = false;
837 }
838 }
839 }
840 //
841 // Set the shear components of the total force.
842 for (int i = 1; i < dim; ++i)
843 jkr_F_.rdof(i) = sforce.dof(i);
844 //
845 // ROLLING RESISTANCE ==============================================================================================================
846 //
847 DAVect res_M_old = res_M_; // used in energy calculation only
848 if (fr_ == 0.0) {
849 res_M_.fill(0.0);
850 kr_ = 0.0;
851 }
852 else {
853 DAVect angStiff(0.0);
854 DAVect MomentInc(0.0);
855 kr_ = ks_tangent * rbar_square_; // update rolling stiffness based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
856#if DIM==3
857 angStiff.rx() = 0.0;
858 angStiff.ry() = kr_;
859#endif
860 angStiff.rz() = kr_;
861 MomentInc = ang * angStiff * (-1.0);
862 res_M_ += MomentInc;
863 if (state->canFail_) {
864 // Account for bending strength
865 double rmax = std::abs(fr_*(jkr_F_.rdof(0) + 2.0*pull_off_f_)); // Using the same normal force as in the shear limit
866 double rmag = res_M_.mag();
867 if (rmag > rmax) {
868 double fac = rmax / rmag;
869 res_M_ *= fac;
870 res_S_ = true;
871 }
872 else {
873 res_S_ = false;
874 }
875 }
876 }
877 //
878 // STIFFNESS FOR TIME STEP UPDATES ============================================================================================
879 //
880 // set effective stiffness in the normal and shear directions for timestep calculation
881 effectiveTranslationalStiffness_ = DVect2(kn_tangent, ks_tangent);
882 // set the effective rotational stiffness (bending only)
883 effectiveRotationalStiffness_ = DAVect(kr_);
884#if DIM==3
885 effectiveRotationalStiffness_.rx() = 0.0;
886#endif
887 //
888 // DAMPING FORCE ==============================================================================================================
889 //
890 // Account for dashpot forces if the dashpot structure has been defined.
891 if (dpProps_) {
892 dpProps_->dp_F_.fill(0.0);
893 double vcn(0.0), vcs(0.0);
894 // Calculate the damping coefficients.
895 setDampCoefficients(state->inertialMass_, &kn_tangent, &ks_tangent, &vcn, &vcs); // this is based on the current tangent stiffness
896 // First damp the shear components
897 for (int i = 1; i < dim; ++i)
898 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
899 // Damp the normal component
900 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
901 //
902 if (jkr_S_ && dpProps_->dp_mode_ == 1) {
903 // Limit in shear if sliding.
904 double dfn = dpProps_->dp_F_.rx();
905 dpProps_->dp_F_.fill(0.0);
906 dpProps_->dp_F_.rx() = dfn;
907 }
908 // Correct effective translational stiffness
909 DVect2 correct(1.0);
910 if (dpProps_->dp_nratio_)
911 correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
912 if (dpProps_->dp_sratio_)
913 correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
914 effectiveTranslationalStiffness_ /= (correct*correct);
915 }
916
917 //Compute energies if energy tracking has been enabled.
918 if (state->trackEnergy_) {
919 assert(energies_);
920 // calculate the strain energy increment in the normal direction
921 energies_->estrain_ -= 0.5*(jkr_F_old.dof(0) + jkr_F_.dof(0))*trans.dof(0);// under loading, trans.x < 0, hence the negative sign
922 if (ks_) {
923 DVect u_s_elastic = trans; // set the elastic displacement increment equal to the total displacement increment
924 u_s_elastic.rx() = 0.0; // set the normal component to zero: u_s_elatic = [0, trans_shear_1, trans_shear_2]
925 DVect shearF = jkr_F_; // set the shear force equal to the total jkr force (including the normal component)
926 shearF.rx() = 0.0; // set the normal component to zero: shearF = [0, shear_force_1, shear_force_2]
927 jkr_F_old.rx() = 0.0; // set normal component of the previous force equal to zero
928 DVect avg_F_s = (shearF + jkr_F_old)*0.5; // average shear force vector
929 if (jkr_S_) { // if sliding, calculate the slip energy and accumulate it
930 DVect u_s_total = u_s_elastic; // total shear displacement increment
931 u_s_elastic = (shearF - jkr_F_old) / ks_tangent; // elastic shear displacement increment
932 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s_total + u_s_elastic))); // where (u_s_total + u_s_elatic) is the slip displacment increment (due to the sign convention, the terms are added up)
933 }
934 energies_->estrain_ -= avg_F_s | u_s_elastic; // add the elastic component (if any - if previous force equals current force, the elastic incrment is zero)
935 }
936 // Add the rolling resistance energy contributions. // done incrementally since the stiffness kr_ changes with the tangent shear stiffness (non-linearly)
937 if (kr_) {
938 DAVect t_s_elastic = ang; // set the elastic rotation increment equal to the total angle increment
939 DAVect avg_M = (res_M_ + res_M_old)*0.5; // average moment from this time step and the previous
940 if (res_S_) { // if sliding, calculate the slip energy and accumulate it
941 t_s_elastic = (res_M_ - res_M_old) / kr_; // elastic angle increment
942 energies_->errslip_ -= std::min(0.0, (avg_M | (ang + t_s_elastic)));// where (ang + t_s_elatic) is the slip rotation increment (due to the sign convention, the terms are added up)
943 }
944 energies_->errstrain_ -= avg_M | t_s_elastic; // add the elastic component (if any - if previous force equals current force, the elastic incrment is zero)
945 }
946 // Calculate damping energy (accumulated) if the dashpots are active.
947 if (dpProps_) {
948 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
949 }
950 }
951
952 // This is just a sanity check to ensure, in debug mode, that the force isn't wonky.
953 assert(jkr_F_ == jkr_F_);
954 return true;
955 }
956
957 void ContactModelJKR::setForce(const DVect &, IContact *) {
958 //
959 }
960
961 DVect ContactModelJKR::getForce() const {
962 DVect ret(jkr_F_);
963 if (dpProps_)
964 ret += dpProps_->dp_F_;
965 return ret;
966 }
967
968 DAVect ContactModelJKR::getMomentOn1(const IContactMechanical *c) const {
969 DAVect ret(res_M_);
970 if (c) {
971 DVect force = getForce();
972 c->updateResultingTorqueOn1Local(force, &ret);
973 }
974 return ret;
975 }
976
977 DAVect ContactModelJKR::getMomentOn2(const IContactMechanical *c) const {
978 DAVect ret(res_M_);
979 if (c) {
980 DVect force = getForce();
981 c->updateResultingTorqueOn2Local(force, &ret);
982 }
983 return ret;
984 }
985
986 DAVect ContactModelJKR::getModelMomentOn1() const {
987 DAVect ret(res_M_);
988 return ret;
989 }
990
991 DAVect ContactModelJKR::getModelMomentOn2() const {
992 DAVect ret(res_M_);
993 return ret;
994 }
995
996 void ContactModelJKR::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
997 ret->clear();
998 ret->push_back({"jkr_shear",ScalarInfo});
999 ret->push_back({"jkr_poiss",ScalarInfo});
1000 ret->push_back({"fric",ScalarInfo});
1001 ret->push_back({"jkr_force",VectorInfo});
1002 ret->push_back({"jkr_slip",ScalarInfo});
1003 ret->push_back({"dp_nratio",ScalarInfo});
1004 ret->push_back({"dp_sratio",ScalarInfo});
1005 ret->push_back({"dp_mode",ScalarInfo});
1006 ret->push_back({"dp_force",VectorInfo});
1007 ret->push_back({"rr_fric",ScalarInfo});
1008 ret->push_back({"rr_moment",VectorInfo});
1009 ret->push_back({"rr_slip",ScalarInfo});
1010 ret->push_back({"ks_fac",ScalarInfo});
1011 ret->push_back({"surf_adh",ScalarInfo});
1012 ret->push_back({"active_mode",ScalarInfo});
1013 ret->push_back({"a0",ScalarInfo});
1014 ret->push_back({"pull_off_f",ScalarInfo});
1015 ret->push_back({"tear_off_d",ScalarInfo});
1016 }
1017
1018 void ContactModelJKR::objectPropValues(std::vector<double> *ret,const IContact *) const {
1019 FP_S;
1020 ret->clear();
1021 ret->push_back(shear_);
1022 ret->push_back(poiss_);
1023 ret->push_back(fric_);
1024 ret->push_back(jkr_F_.mag());
1025 ret->push_back(jkr_S_);
1026 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
1027 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
1028 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
1029 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
1030 ret->push_back(res_fric_);
1031 ret->push_back(res_M_.mag());
1032 ret->push_back(res_S_);
1033 ret->push_back(ks_fac_);
1034 ret->push_back(surf_adh_);
1035 ret->push_back(active_mode_);
1036 ret->push_back(a0_);
1037 ret->push_back(pull_off_f_);
1038 ret->push_back(tear_off_d_);
1039 FP_S;
1040 }
1041
1042 void ContactModelJKR::setDampCoefficients(const double &mass, double *kn_tangent, double *ks_tangent, double *vcn, double *vcs) {
1043 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(*kn_tangent));
1044 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(*ks_tangent));
1045 }
1046
1047} // namespace cmodelsxd
1048// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
