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 QString 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 QString 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 QString 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_ = NEWC(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 QString getFishCallEvents() const override {
125 return
126 "contact_activated"
127 ",slip_change";
128 }
129
130 // Return the specified contact model property.
131 QVariant 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 QVariant &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 QString &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 NEWC(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<QString,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_ = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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 QVariant 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 QVariant 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 QVariant();
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 QVariant &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_ = NEWC(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_ = NEWC(dpProps());
314 dpProps_->dp_sratio_ = val;
315 return true;
316 }
317 case kwDpMode: {
318 if (!v.canConvert<int>())
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_ = NEWC(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_ = NEWC(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<int>())
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*M_PI*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*M_PI*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*M_PI*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 QString &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 QStringList availableProperties = getProperties().simplified().replace(" ", "").split(",", Qt::SkipEmptyParts);
544 QRegExp rx(name, Qt::CaseInsensitive);
545 int idx = availableProperties.indexOf(rx) + 1;
546 bool ret = false;
547
548 if (idx <= 0)
549 return ret;
550
551 switch (idx) {
552 case kwShear: { //Shear
553 if (inheritanceField_ & shearMask)
554 ret = true; // return 'true' to ensure 'validate()' is called to affect the change in Shear Modulus through 'updateEndStiffCoef()'
555 break;
556 }
557 case kwPoiss: { //Poisson's
558 if (inheritanceField_ & poissMask)
559 ret =true; // return 'true' to ensure 'validate()' is called to affect the change in Poisson's ratio through 'updateEndStiffCoef()'
560 break;
561 }
562 case kwFric: { //fric
563 if (inheritanceField_ & fricMask)
564 updateEndFric(c); // friction does not influence any of the other parameters or time step size, so directly update the contact model friction
565 break;
566 }
567 case kwResFric: { //rr_fric
568 if (inheritanceField_ & resFricMask)
569 ret = true; // return 'true' to ensure 'validate()' is called to affect the change in rolling friction through 'updateEndResFric()' and 'fr_'
570 break;
571 }
572 case kwSurfAdh: { //surf_adh
573 if (inheritanceField_ & surfAdhMask)
574 ret = true; // return 'true' to ensure 'validate()' is called to affect the change in surface adhesion through 'updateEndSurfAdh()'
575 break;
576 }
577 //
578 }
579 return ret;
580 }
581
582 static const QString gstr("jkr_shear");
583 static const QString nustr("jkr_poiss");
584 bool ContactModelJKR::updateEndStiffCoef(const IContactMechanical *con) {
585 assert(con);
586 double g1 = shear_;
587 double g2 = shear_;
588 double nu1 = poiss_;
589 double nu2 = poiss_;
590 QVariant vg1 = con->getEnd1()->getProperty(gstr);
591 QVariant vg2 = con->getEnd2()->getProperty(gstr);
592 QVariant vnu1 = con->getEnd1()->getProperty(nustr);
593 QVariant vnu2 = con->getEnd2()->getProperty(nustr);
594 if (vg1.isValid() && vg2.isValid()) {
595 g1 = vg1.toDouble();
596 g2 = vg2.toDouble();
597 if (g1 <= 0.0 || g2 <= 0.0)
598 throw Exception("Negative or zero shear modulus not allowed in jkr contact model");
599 }
600 if (vnu1.isValid() && vnu2.isValid()) {
601 nu1 = vnu1.toDouble();
602 nu2 = vnu2.toDouble();
603 if (nu1 < 0.0 || nu1 > 0.5 || nu2 < 0.0 || nu2 > 0.5)
604 throw Exception("Poisson ratio should be in range [0,0.5] in jkr contact model");
605 }
606 if (g1*g2 == 0.0) return false;
607 double es = 1.0 / ((1.0 - nu1) / (2.0*g1) + (1.0 - nu2) / (2.0*g2));
608 double gs = 1.0 / ((2.0 - nu1) / g1 + (2.0 - nu2) / g2);
609 poiss_ = (4.0*gs - es) / (2.0*gs - es);
610 shear_ = 2.0*gs*(2 - poiss_);
611 if (shear_ <= 0.0)
612 throw Exception("Negative or zero shear modulus not allowed in jkr contact model");
613 if (poiss_ < 0.0 || poiss_ > 0.5)
614 throw Exception("Poisson ratio should be in range [0,0.5] in jkr contact model");
615 return true;
616 }
617
618 static const QString fricstr("fric");
619 bool ContactModelJKR::updateEndFric(const IContactMechanical *con) {
620 assert(con);
621 QVariant v1 = con->getEnd1()->getProperty(fricstr);
622 QVariant v2 = con->getEnd2()->getProperty(fricstr);
623 if (!v1.isValid() || !v2.isValid())
624 return false;
625 double fric1 = std::max(0.0, v1.toDouble());
626 double fric2 = std::max(0.0, v2.toDouble());
627 double val = fric_;
628 fric_ = std::min(fric1, fric2);
629 if (fric_ < 0.0)
630 throw Exception("Negative friction value not allowed in jkr contact model");
631 return ((fric_ != val));
632 }
633
634 static const QString rfricstr("rr_fric");
635 bool ContactModelJKR::updateEndResFric(const IContactMechanical *con) {
636 assert(con);
637 QVariant v1 = con->getEnd1()->getProperty(rfricstr);
638 QVariant v2 = con->getEnd2()->getProperty(rfricstr);
639 if (!v1.isValid() || !v2.isValid())
640 return false;
641 double rfric1 = std::max(0.0, v1.toDouble());
642 double rfric2 = std::max(0.0, v2.toDouble());
643 double val = res_fric_;
644 res_fric_ = std::min(rfric1, rfric2);
645 if (res_fric_ < 0.0)
646 throw Exception("Negative rolling friction value not allowed in jkr contact model");
647 return ((res_fric_ != val));
648 }
649
650 static const QString surfadhstr("surf_adh");
651 bool ContactModelJKR::updateEndSurfAdh(const IContactMechanical *con) {
652 assert(con);
653 QVariant v1 = con->getEnd1()->getProperty(surfadhstr);
654 QVariant v2 = con->getEnd2()->getProperty(surfadhstr);
655 if (!v1.isValid() || !v2.isValid())
656 return false;
657 double surfadh1 = std::max(0.0, v1.toDouble());
658 double surfadh2 = std::max(0.0, v2.toDouble());
659 if (surfadh1 <= 0.0 || surfadh2 <= 0.0)
660 throw Exception("Negative or zero surface adhesion not allowed in 'jkr' contact model");
661 double val = surf_adh_;
662 surf_adh_ = 0.5*(surfadh1 + surfadh2); // new surface adhesion enery
663 return ((surf_adh_ != val));
664 }
665
666 void ContactModelJKR::updateEffectiveStiffness(ContactModelMechanicalState *state) {
667
668 double overlap = rgap_ - state->gap_;
669 if (overlap <= -distance_active_) { // if contact was just created but not yet activated, distance_active_ == 0.
670 // The assumption below is based on the Hill contact model in PFC
671 // 0.01% of Diameter = 0.0001*D
672 // For monodisperse system, r_hertz_ = 0.5*R where R is the particle radius
673 // r_hertz_ = 0.5*(0.5*D) = 0.25*D where D is the particle diameter
674 // D = 4*r_hertz_
675 // 0.01%*D = 0.0001*D = 0.0001*4*r_hertz_ = 0.0004*r_hertz_
676 overlap = 0.0004*r_hertz_;
677 }
678 //
679 DVect2 force_radius = JKR_Analytical(overlap); // analytical patch radius and force for this (assumed) overlap
680 double patch_radius = force_radius.dof(1); // contact patch radius
681 // 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
682 double kn_tangent = kn_ * patch_radius; // tangent stiffness in normal direction
683 double ks_tangent = ks_ * patch_radius; // tangent stiffness in the shear direction - tis is not an assumption, but the correct solution for JKR
684 effectiveTranslationalStiffness_ = DVect2(kn_tangent, ks_tangent); // set using the tangent stiffness
685 if (res_fric_ > 0.0) { // rolling
686 kr_ = ks_tangent * rbar_square_; // based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
687 effectiveRotationalStiffness_ = DAVect(kr_); // effective rotational stiffness (bending only)
688#if DIM==3
689 effectiveRotationalStiffness_.rx() = 0.0;
690#endif
691 }
692 // correction if viscous damping active
693 if (dpProps_) {
694 DVect2 correct(1.0);
695 if (dpProps_->dp_nratio_)
696 correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
697 if (dpProps_->dp_sratio_)
698 correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
699 effectiveTranslationalStiffness_ /= (correct*correct);
700 }
701 }
702 //
703 // --------------------------------------------------------------------------------------------------------------------------------------------------------
704 //
705 DVect2 ContactModelJKR::JKR_Analytical(double overlap) { // analytical solution according to Parteli et al. 2014
706 DVect2 ret(0.0);
707 // reduces to Hertz model
708 if (surf_adh_ == 0.0) {
709 double patch_radius = std::sqrt(r_hertz_*overlap);
710 ret.rdof(0) = (2.0 / 3.0)*kn_*patch_radius*patch_radius*patch_radius / r_hertz_;
711 ret.rdof(1) = patch_radius;
712 return ret;
713 }
714 // otherwise jkr is active ...
715 double rad_over = r_hertz_ * overlap;
716 double c0 = rad_over * rad_over;
717 double c2 = -2.0*rad_over;
718 double P = -1.0 / 12.0*c2*c2 - c0;
719 double Q = -1.0 / 108.0*c2*c2*c2 + c0 * c2 / 3.0 - c1_*c1_ / 8.0;
720 double U = pow((-0.5*Q + sqrt(0.25*Q*Q + P*P*P / 27.0)), (1.0 / 3.0));
721 double s = 0.0;
722 if (P == 0.0) {
723 s = -5.0*c2 / 6.0 - pow(Q, (1.0 / 3.0));
724 }
725 else {
726 s = -5.0*c2 / 6.0 + U - P / (3.0*U);
727 }
728 double w = sqrt(c2 + 2.0*s);
729 double lambda = c1_ / (2.0*w);
730 double patch_radius = 0.5*(w + sqrt(w*w - 4.0*(c2 + s + lambda)));
731 // reduces to Hertz model due to very small values for surface adhesion
732 if (std::isnan(patch_radius)) {
733 patch_radius = std::sqrt(r_hertz_*overlap);
734 ret.rdof(0) = (2.0 / 3.0)*kn_*patch_radius*patch_radius*patch_radius / r_hertz_;
735 ret.rdof(1) = patch_radius;
736 return ret;
737 }
738 // force = 4*E*(a^3)/(3*R) - sqrt(16*pi*surf_adh*E*a^3);
739 // with E = kn/2 - where E is the effective contact Young's modulus - see updateStiffCoef(const IContactMechanical *con)
740 // the above reduces to the following for the JKR force ...
741 double patch_radius3 = patch_radius* patch_radius*patch_radius;
742 ret.rdof(0) = (2.0 / 3.0)*kn_*patch_radius3 / r_hertz_ - sqrt(8.0*M_PI*surf_adh_*kn_*patch_radius3); // set the total force
743 ret.rdof(1) = patch_radius; // set the patch contact radius
744 return ret;
745 }
746 //
747 // --------------------------------------------------------------------------------------------------------------------------------------------------------
748 //
749 bool ContactModelJKR::forceDisplacementLaw(ContactModelMechanicalState *state, const double ×tep) {
750 assert(state);
751 if (shear_ <= 0.0) {
752 throw Exception("'shear' must be specified using a value larger than zero in the 'jkr' model.");
753 }
754 if (surf_adh_ <= 0.0) {
755 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!");
756 }
757 //
758 // Current overlap
759 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.
760 // Relative translational increment
761 DVect trans = state->relativeTranslationalIncrement_;
762 //
763 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
764 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
765 }
766 //
767 // The contact was just activated from an inactive state
768 // Trigger the FISH callback if one is hooked up to the contact_activated event.
769 if (cmEvents_[fActivated] >= 0) {
770 // An FArray of QVariant is returned and these will be passed
771 // to the FISH function as an array of FISH symbols as the second
772 // argument to the FISH callback function.
773 auto c = state->getContact();
774 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
775 IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
776 fi->setCMFishCallArguments(c, arg, cmEvents_[fActivated]);
777
778 }
779 //
780 // Angular dispacement increment.
781 DAVect ang = state->relativeAngularIncrement_;
782 //
783 DVect jkr_F_old = jkr_F_; // store for energy computations only
784 //
785 // NORMAL FORCE ==============================================================================================================
786 //
787 DVect2 force_radius = JKR_Analytical(overlap); // analytical force and patch radius
788 jkr_F_.rdof(0) = force_radius.dof(0); // current contact force in the normal direction
789 double patch_radius = force_radius.dof(1); // current contact patch radius
790 //
791 double kn_tangent = kn_ * patch_radius;
792 if (trans.dof(0))
793 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
794 //
795 // SHEAR FORCE ==============================================================================================================
796 //
797 double ks_tangent = ks_* patch_radius; // tangent stiffness in the shear direction
798 DVect sforce(0.0);
799 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
800 // Loop over the shear components (note: the 0 component is the normal component) and calculate the shear force.
801 for (int i = 1; i < dim; ++i)
802 sforce.rdof(i) = jkr_F_.dof(i) - trans.dof(i) * ks_tangent; // shear force update
803 //
804 // The canFail flag corresponds to whether or not the contact can undergo non-linear
805 // force-displacement response. If the SOLVE ELASTIC command is given then the
806 // canFail state is set to FALSE. Otherwise it is always TRUE.
807 if (state->canFail_) {
808 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
809 double crit = fric_ * std::abs(jkr_F_.rdof(0) + 2.0*pull_off_f_);
810 // The is the magnitude of the shear force.
811 double sfmag = sforce.mag();
812 // Sliding occurs when the magnitude of the shear force is greater than the critical value.
813 if (sfmag > crit) {
814 // Lower the shear force to the critical value for sliding.
815 double rat = crit / sfmag;
816 sforce *= rat;
817 // Handle the slip_change event if one has been hooked up. Sliding has commenced.
818 if (!jkr_S_ && cmEvents_[fSlipChange] >= 0) {
819 auto c = state->getContact();
820 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
821 fish::Parameter((qint64)0) };
822 IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
823 fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
824 }
825 jkr_S_ = true;
826 }
827 else {
828 // Handle the slip_change event if one has been hooked up and
829 // the contact was previously sliding. Sliding has ceased.
830 if (jkr_S_) {
831 if (cmEvents_[fSlipChange] >= 0) {
832 auto c = state->getContact();
833 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
834 fish::Parameter((qint64)1) };
835 IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
836 fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
837 }
838 jkr_S_ = false;
839 }
840 }
841 }
842 //
843 // Set the shear components of the total force.
844 for (int i = 1; i < dim; ++i)
845 jkr_F_.rdof(i) = sforce.dof(i);
846 //
847 // ROLLING RESISTANCE ==============================================================================================================
848 //
849 DAVect res_M_old = res_M_; // used in energy calculation only
850 if (fr_ == 0.0) {
851 res_M_.fill(0.0);
852 kr_ = 0.0;
853 }
854 else {
855 DAVect angStiff(0.0);
856 DAVect MomentInc(0.0);
857 kr_ = ks_tangent * rbar_square_; // update rolling stiffness based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
858#if DIM==3
859 angStiff.rx() = 0.0;
860 angStiff.ry() = kr_;
861#endif
862 angStiff.rz() = kr_;
863 MomentInc = ang * angStiff * (-1.0);
864 res_M_ += MomentInc;
865 if (state->canFail_) {
866 // Account for bending strength
867 double rmax = std::abs(fr_*(jkr_F_.rdof(0) + 2.0*pull_off_f_)); // Using the same normal force as in the shear limit
868 double rmag = res_M_.mag();
869 if (rmag > rmax) {
870 double fac = rmax / rmag;
871 res_M_ *= fac;
872 res_S_ = true;
873 }
874 else {
875 res_S_ = false;
876 }
877 }
878 }
879 //
880 // STIFFNESS FOR TIME STEP UPDATES ============================================================================================
881 //
882 // set effective stiffness in the normal and shear directions for timestep calculation
883 effectiveTranslationalStiffness_ = DVect2(kn_tangent, ks_tangent);
884 // set the effective rotational stiffness (bending only)
885 effectiveRotationalStiffness_ = DAVect(kr_);
886#if DIM==3
887 effectiveRotationalStiffness_.rx() = 0.0;
888#endif
889 //
890 // DAMPING FORCE ==============================================================================================================
891 //
892 // Account for dashpot forces if the dashpot structure has been defined.
893 if (dpProps_) {
894 dpProps_->dp_F_.fill(0.0);
895 double vcn(0.0), vcs(0.0);
896 // Calculate the damping coefficients.
897 setDampCoefficients(state->inertialMass_, &kn_tangent, &ks_tangent, &vcn, &vcs); // this is based on the current tangent stiffness
898 // First damp the shear components
899 for (int i = 1; i < dim; ++i)
900 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
901 // Damp the normal component
902 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
903 //
904 if (jkr_S_ && dpProps_->dp_mode_ == 1) {
905 // Limit in shear if sliding.
906 double dfn = dpProps_->dp_F_.rx();
907 dpProps_->dp_F_.fill(0.0);
908 dpProps_->dp_F_.rx() = dfn;
909 }
910 // Correct effective translational stiffness
911 DVect2 correct(1.0);
912 if (dpProps_->dp_nratio_)
913 correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
914 if (dpProps_->dp_sratio_)
915 correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
916 effectiveTranslationalStiffness_ /= (correct*correct);
917 }
918
919 //Compute energies if energy tracking has been enabled.
920 if (state->trackEnergy_) {
921 assert(energies_);
922 // calculate the strain energy increment in the normal direction
923 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
924 if (ks_) {
925 DVect u_s_elastic = trans; // set the elastic displacement increment equal to the total displacement increment
926 u_s_elastic.rx() = 0.0; // set the normal component to zero: u_s_elatic = [0, trans_shear_1, trans_shear_2]
927 DVect shearF = jkr_F_; // set the shear force equal to the total jkr force (including the normal component)
928 shearF.rx() = 0.0; // set the normal component to zero: shearF = [0, shear_force_1, shear_force_2]
929 jkr_F_old.rx() = 0.0; // set normal component of the previous force equal to zero
930 DVect avg_F_s = (shearF + jkr_F_old)*0.5; // average shear force vector
931 if (jkr_S_) { // if sliding, calculate the slip energy and accumulate it
932 DVect u_s_total = u_s_elastic; // total shear displacement increment
933 u_s_elastic = (shearF - jkr_F_old) / ks_tangent; // elastic shear displacement increment
934 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)
935 }
936 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)
937 }
938 // Add the rolling resistance energy contributions. // done incrementally since the stiffness kr_ changes with the tangent shear stiffness (non-linearly)
939 if (kr_) {
940 DAVect t_s_elastic = ang; // set the elastic rotation increment equal to the total angle increment
941 DAVect avg_M = (res_M_ + res_M_old)*0.5; // average moment from this time step and the previous
942 if (res_S_) { // if sliding, calculate the slip energy and accumulate it
943 t_s_elastic = (res_M_ - res_M_old) / kr_; // elastic angle increment
944 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)
945 }
946 energies_->errstrain_ -= avg_M | t_s_elastic; // add the elastic component (if any - if previous force equals current force, the elastic incrment is zero)
947 }
948 // Calculate damping energy (accumulated) if the dashpots are active.
949 if (dpProps_) {
950 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
951 }
952 }
953
954 // This is just a sanity check to ensure, in debug mode, that the force isn't wonky.
955 assert(jkr_F_ == jkr_F_);
956 return true;
957 }
958
959 void ContactModelJKR::setForce(const DVect &, IContact *) {
960 //
961 }
962
963 DVect ContactModelJKR::getForce() const {
964 DVect ret(jkr_F_);
965 if (dpProps_)
966 ret += dpProps_->dp_F_;
967 return ret;
968 }
969
970 DAVect ContactModelJKR::getMomentOn1(const IContactMechanical *c) const {
971 DAVect ret(res_M_);
972 if (c) {
973 DVect force = getForce();
974 c->updateResultingTorqueOn1Local(force, &ret);
975 }
976 return ret;
977 }
978
979 DAVect ContactModelJKR::getMomentOn2(const IContactMechanical *c) const {
980 DAVect ret(res_M_);
981 if (c) {
982 DVect force = getForce();
983 c->updateResultingTorqueOn2Local(force, &ret);
984 }
985 return ret;
986 }
987
988 DAVect ContactModelJKR::getModelMomentOn1() const {
989 DAVect ret(res_M_);
990 return ret;
991 }
992
993 DAVect ContactModelJKR::getModelMomentOn2() const {
994 DAVect ret(res_M_);
995 return ret;
996 }
997
998 void ContactModelJKR::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
999 ret->clear();
1000 ret->push_back({"jkr_shear",ScalarInfo});
1001 ret->push_back({"jkr_poiss",ScalarInfo});
1002 ret->push_back({"fric",ScalarInfo});
1003 ret->push_back({"jkr_force",VectorInfo});
1004 ret->push_back({"jkr_slip",ScalarInfo});
1005 ret->push_back({"dp_nratio",ScalarInfo});
1006 ret->push_back({"dp_sratio",ScalarInfo});
1007 ret->push_back({"dp_mode",ScalarInfo});
1008 ret->push_back({"dp_force",VectorInfo});
1009 ret->push_back({"rr_fric",ScalarInfo});
1010 ret->push_back({"rr_moment",VectorInfo});
1011 ret->push_back({"rr_slip",ScalarInfo});
1012 ret->push_back({"ks_fac",ScalarInfo});
1013 ret->push_back({"surf_adh",ScalarInfo});
1014 ret->push_back({"active_mode",ScalarInfo});
1015 ret->push_back({"a0",ScalarInfo});
1016 ret->push_back({"pull_off_f",ScalarInfo});
1017 ret->push_back({"tear_off_d",ScalarInfo});
1018 }
1019
1020 void ContactModelJKR::objectPropValues(std::vector<double> *ret,const IContact *) const {
1021 ret->clear();
1022 ret->push_back(shear_);
1023 ret->push_back(poiss_);
1024 ret->push_back(fric_);
1025 ret->push_back(jkr_F_.mag());
1026 ret->push_back(jkr_S_);
1027 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
1028 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
1029 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
1030 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
1031 ret->push_back(res_fric_);
1032 ret->push_back(res_M_.mag());
1033 ret->push_back(res_S_);
1034 ret->push_back(ks_fac_);
1035 ret->push_back(surf_adh_);
1036 ret->push_back(active_mode_);
1037 ret->push_back(a0_);
1038 ret->push_back(pull_off_f_);
1039 ret->push_back(tear_off_d_);
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: Oct 31, 2024 |
