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