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