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