Adhesive Rolling Resistance Linear Model
See this page for the documentation of this contact model.
contactmodelarrlinear.h
1#pragma once
2// contactmodelARRLinear.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef ARRLINEAR_LIB
7# define ARRLINEAR_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define ARRLINEAR_EXPORT
10#else
11# define ARRLINEAR_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelARRLinear : public ContactModelMechanical {
18 public:
19 // Constructor: Set default values for contact model properties.
20 ARRLINEAR_EXPORT ContactModelARRLinear();
21 ARRLINEAR_EXPORT ContactModelARRLinear(const ContactModelARRLinear &) noexcept;
22 ARRLINEAR_EXPORT const ContactModelARRLinear & operator=(const ContactModelARRLinear &);
23 ARRLINEAR_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
24 // Destructor, called when contact is deleted: free allocated memory, etc.
25 ARRLINEAR_EXPORT virtual ~ContactModelARRLinear();
26 // Contact model name (used as keyword for commands and FISH).
27 string getName() const override { return "arrlinear"; }
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 void setIndex(int i) override { index_=i;}
32 int getIndex() const override {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 uint32 getMinorVersion() const override;
38 // Copy the state information to a newly created contact model.
39 // Provide access to state information, for use by copy method.
40 void copy(const ContactModel *c) override;
41 // Provide save-restore capability for the state information.
42 void archive(ArchiveStream &) override;
43 // Enumerator for the properties.
44 enum PropertyKeys {
45 kwKn=1
46 , kwKs
47 , kwFric
48 , kwLinF
49 , kwLinS
50 , kwLinMode
51 , kwRGap
52 , kwEmod
53 , kwKRatio
54 , kwDpNRatio
55 , kwDpSRatio
56 , kwDpMode
57 , kwDpF
58 , kwResFric
59 , kwResMoment
60 , kwResS
61 , kwResKr
62 , kwAdhesiveF0
63 , kwAdhesiveD0
64 , kwAdhesiveF
65 , kwUserArea
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 string getProperties() const override {
71 return "kn"
72 ",ks"
73 ",fric"
74 ",lin_force"
75 ",lin_slip"
76 ",lin_mode"
77 ",rgap"
78 ",emod"
79 ",kratio"
80 ",dp_nratio"
81 ",dp_sratio"
82 ",dp_mode"
83 ",dp_force"
84 ",rr_fric"
85 ",rr_moment"
86 ",rr_slip"
87 ",rr_kr"
88 ",adh_f0"
89 ",adh_d0"
90 ",adh_force"
91 ",user_area";
92 }
93 // Enumerator for the energies.
94 enum EnergyKeys {
95 kwEStrain=1
96 , kwERRStrain
97 , kwESlip
98 , kwERRSlip
99 , kwEDashpot
100 , kwEAdhesive
101 };
102 // Contact model energy names in a comma separated list. The order corresponds with
103 // the order of the EnergyKeys enumerator above.
104 string getEnergies() const override {
105 return " energy-strain"
106 ",energy-rrstrain"
107 ",energy-slip"
108 ",energy-rrslip"
109 ",energy-dashpot"
110 ",energy-adhesive";
111 }
112 // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
113 double getEnergy(uint32 i) const override;
114 // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1)
115 // returns whether or not the estrain energy is accumulated which is false).
116 bool getEnergyAccumulate(uint32 i) const override;
117 // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
118 void setEnergy(uint32 i,const double &d) override; // Base 1
119 // Activate the energy. This is only called if the energy tracking is enabled.
120 void activateEnergy() override { if (energies_) return; energies_ = NEW Energies();}
121 // Returns whether or not the energy tracking has been enabled for this contact.
122 bool getEnergyActivated() const override {return (energies_ != 0);}
123
124 // Enumerator for contact model related FISH callback events.
125 enum FishCallEvents {
126 fActivated=0
127 ,fSlipChange
128 };
129 // Contact model FISH callback event names in a comma separated list. The order corresponds with
130 // the order of the FishCallEvents enumerator above.
131 string getFishCallEvents() const override {
132 return
133 "contact_activated"
134 ",slip_change";
135 }
136
137 // Return the specified contact model property.
138 base::Property getProperty(uint32 i,const IContact *) const;
139 // The return value denotes whether or not the property corresponds to the global
140 // or local coordinate system (TRUE: global system, FALSE: local system). The
141 // local system is the contact-plane system (nst) defined as follows.
142 // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
143 // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
144 // and tc are unit vectors in directions of the nst axes.
145 // This is used when rendering contact model properties that are vectors.
146 bool getPropertyGlobal(uint32 i) const override;
147 // Set the specified contact model property, ensuring that it is of the correct type
148 // and within the correct range --- if not, then throw an exception.
149 // The return value denotes whether or not the update has affected the timestep
150 // computation (by having modified the translational or rotational tangent stiffnesses).
151 // If true is returned, then the timestep will be recomputed.
152 bool setProperty(uint32 i,const base::Property &v,IContact *) override;
153 // The return value denotes whether or not the property is read-only
154 // (TRUE: read-only, FALSE: read-write).
155 bool getPropertyReadOnly(uint32 i) const override;
156
157 // The return value denotes whether or not the property is inheritable
158 // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
159 // the endPropertyUpdated method.
160 bool supportsInheritance(uint32 i) const override;
161 // Return whether or not inheritance is enabled for the specified property.
162 bool getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
163 // Set the inheritance flag for the specified property.
164 void setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
165
166 // Enumerator for contact model methods.
167 enum MethodKeys { kwDeformability=1, kwArea};
168 // Contact model methoid names in a comma separated list. The order corresponds with
169 // the order of the MethodKeys enumerator above.
170 string getMethods() const override { return "deformability,area";}
171 // Return a comma seprated list of the contact model method arguments (base 1).
172 string getMethodArguments(uint32 i) const override;
173 // Set contact model method arguments (base 1).
174 // The return value denotes whether or not the update has affected the timestep
175 // computation (by having modified the translational or rotational tangent stiffnesses).
176 // If true is returned, then the timestep will be recomputed.
177 bool setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con=0) override;
178
179 // Prepare for entry into ForceDispLaw. The validate function is called when:
180 // (1) the contact is created, (2) a property of the contact that returns a true via
181 // the setProperty method has been modified and (3) when a set of cycles is executed
182 // via the {cycle N} command.
183 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
184 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
185 // The endPropertyUpdated method is called whenever a surface property (with a name
186 // that matches an inheritable contact model property name) of one of the contacting
187 // pieces is modified. This allows the contact model to update its associated
188 // properties. The return value denotes whether or not the update has affected
189 // the time step computation (by having modified the translational or rotational
190 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
191 bool endPropertyUpdated(const string &name,const IContactMechanical *c) override;
192 // The forceDisplacementLaw function is called during each cycle. Given the relative
193 // motion of the two contacting pieces (via
194 // state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
195 // state->relativeAngularIncrement_ (Dtt, Dtbs, Dtbt)
196 // Ddn : relative normal-displacement increment, Ddn > 0 is opening
197 // Ddss : relative shear-displacement increment (s-axis component)
198 // Ddst : relative shear-displacement increment (t-axis component)
199 // Dtt : relative twist-rotation increment
200 // Dtbs : relative bend-rotation increment (s-axis component)
201 // Dtbt : relative bend-rotation increment (t-axis component)
202 // The relative displacement and rotation increments:
203 // Dd = Ddn*nc + Ddss*sc + Ddst*tc
204 // Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
205 // where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
206 // [see {Table 1: Contact State Variables} in PFC Model Components:
207 // Contacts and Contact Models: Contact Resolution]
208 // ) and the contact properties, this function must update the contact force and
209 // moment.
210 // The force_ is acting on piece 2, and is expressed in the local coordinate system
211 // (defined in getPropertyGlobal) such that the first component positive denotes
212 // compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
213 // in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to
214 // get the total moment.
215 // The return value indicates the contact activity status (TRUE: active, FALSE:
216 // inactive) during the next cycle.
217 // Additional information:
218 // * If state->activated() is true, then the contact has just become active (it was
219 // inactive during the previous time step).
220 // * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
221 // the forceDispLaw handle the case of {state->canFail_ == true}.
222 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
223 bool thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
224 // The getEffectiveXStiffness functions return the translational and rotational
225 // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
226 // the translational tangent shear stiffness is zero (but this stiffness reduction
227 // is typically ignored when computing a stable time step). If the contact model
228 // includes a dashpot, then the translational stiffnesses must be increased (see
229 // Potyondy (2009)).
230 // [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
231 // Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
232 // December 7, 2009.]
233 DVect2 getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
234 DAVect getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
235
236 // Return a new instance of the contact model. This is used in the CMAT
237 // when a new contact is created.
238 ContactModelARRLinear *clone() const override { return NEW ContactModelARRLinear(); }
239 // The getActivityDistance function is called by the contact-resolution logic when
240 // the CMAT is modified. Return value is the activity distance used by the
241 // checkActivity function.
242 double getActivityDistance() const override {return (rgap_ + a_d0_);}
243 // The isOKToDelete function is called by the contact-resolution logic when...
244 // Return value indicates whether or not the contact may be deleted.
245 // If TRUE, then the contact may be deleted when it is inactive.
246 // If FALSE, then the contact may not be deleted (under any condition).
247 bool isOKToDelete() const override { return !isBonded(); }
248 // Zero the forces and moments stored in the contact model. This function is called
249 // when the contact becomes inactive.
250 void resetForcesAndMoments() override { lin_F(DVect(0.0)); dp_F(DVect(0.0));
251 res_M(DAVect(0.0)); a_F(0.0);
252 if (energies_) { energies_->estrain_ = 0.0;
253 energies_->errstrain_ = 0.0; }
254 }
255 void setForce(const DVect &v,IContact *c) override;
256 void setArea(const double &d) override { userArea_ = d; }
257 double getArea() const override { return userArea_; }
258 // The checkActivity function is called by the contact-resolution logic when...
259 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
260 // A contact with the arrlinear model is active if the surface gap is less than
261 // or equal to the attraction range (a_d0_).
262 bool checkActivity(const double &gap) override { return gap <= (rgap_ + a_d0_); }
263
264 // Returns the sliding state (FALSE is returned if not implemented).
265 bool isSliding() const override { return lin_S_; }
266 // Returns the bonding state (FALSE is returned if not implemented).
267 bool isBonded() const override { return false; }
268
269 // Both of these methods are called only for contacts with facets where the wall
270 // resolution scheme is set the full. In such cases one might wish to propagate
271 // contact state information (e.g., shear force) from one active contact to another.
272 // See the Faceted Wall section in the documentation.
273 void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
274 void setNonForcePropsFrom(IContactModel *oldCM) override;
275
276 /// Return the total force that the contact model holds.
277 DVect getForce() const override;
278
279 /// Return the total moment on 1 that the contact model holds
280 DAVect getMomentOn1(const IContactMechanical *) const override;
281
282 /// Return the total moment on 1 that the contact model holds
283 DAVect getMomentOn2(const IContactMechanical *) const override;
284 /// Return moments without torque
285 DAVect getModelMomentOn1() const override;
286 DAVect getModelMomentOn2() const override;
287
288 // Used to efficiently get properties from the contact model for the object pane.
289 // List of properties for the object pane, comma separated.
290 // All properties will be cast to doubles for comparison. No other comparisons
291 // are supported. This may not be the same as the entire property list.
292 // Return property name and type for plotting.
293 void objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *) const override;
294 // All properties cast to doubles - this is what can be compared.
295 void objectPropValues(std::vector<double> *,const IContact *) const override;
296
297 // Methods to get and set properties.
298 const double & kn() const {return kn_;}
299 void kn(const double &d) {kn_=d;}
300 const double & ks() const {return ks_;}
301 void ks(const double &d) {ks_=d;}
302 const double & fric() const {return fric_;}
303 void fric(const double &d) {fric_=d;}
304 const DVect & lin_F() const {return lin_F_;}
305 void lin_F(const DVect &f) { lin_F_=f;}
306 bool lin_S() const {return lin_S_;}
307 void lin_S(bool b) { lin_S_=b;}
308 uint32 lin_mode() const {return lin_mode_;}
309 void lin_mode(uint32 i) { lin_mode_= i;}
310 const double & rgap() const {return rgap_;}
311 void rgap(const double &d) {rgap_=d;}
312 double emod(const IContact *c) const;
313 double kratio() const;
314
315 bool hasDamping() const {return dpProps_ ? true : false;}
316 double dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
317 void dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
318 double dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
319 void dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
320 int dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
321 void dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
322 DVect dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
323 void dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
324
325 bool hasEnergies() const {return energies_ ? true:false;}
326 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
327 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
328 double errstrain() const {return hasEnergies() ? energies_->errstrain_: 0.0;}
329 void errstrain(const double &d) { if(!hasEnergies()) return; energies_->errstrain_=d;}
330 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
331 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
332 double errslip() const {return hasEnergies() ? energies_->errslip_: 0.0;}
333 void errslip(const double &d) { if(!hasEnergies()) return; energies_->errslip_=d;}
334 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
335 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
336 double eadhesive() const {return hasEnergies() ? energies_->eadhesive_: 0.0;}
337 void eadhesive(const double &d) { if(!hasEnergies()) return; energies_->eadhesive_=d;}
338
339 uint32 inheritanceField() const {return inheritanceField_;}
340 void inheritanceField(uint32 i) {inheritanceField_ = i;}
341
342 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
343 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
344 const DAVect & effectiveRotationalStiffness() const {return effectiveRotationalStiffness_;}
345 void effectiveRotationalStiffness(const DAVect &v ) {effectiveRotationalStiffness_=v;}
346
347 // Rolling resistance methods
348 const double & res_fric() const {return res_fric_;}
349 void res_fric(const double &d) {res_fric_=d;}
350 const DAVect & res_M() const {return res_M_;}
351 void res_M(const DAVect &f) { res_M_=f;}
352 bool res_S() const {return res_S_;}
353 void res_S(bool b) { res_S_=b;}
354 const double & kr() const {return kr_;}
355 void kr(const double &d) {kr_=d;}
356 const double & fr() const {return fr_;}
357 void fr(const double &d) {fr_=d;}
358
359 // Adhesive methods
360 const double & a_f0() const {return a_f0_;}
361 void a_f0(const double &d) {a_f0_ = d;}
362 const double & a_d0() const {return a_d0_;}
363 void a_d0(const double &d) {a_d0_ = d;}
364 const double & a_F() const {return a_F_;}
365 void a_F(const double &d) {a_F_ = d;}
366
367 private:
368 // Index - used internally by PFC. Should be set to -1 in the cpp file.
369 static int index_;
370
371 // Structure to store the energies.
372 struct Energies {
373 Energies() : estrain_(0.0), errstrain_(0.0), eslip_(0.0), errslip_(0.0), edashpot_(0.0), eadhesive_(0.0) {}
374 double estrain_; // elastic energy stored in linear group
375 double errstrain_; // elastic energy stored in rolling resistance group
376 double eslip_; // work dissipated by friction
377 double errslip_; // work dissipated by rolling resistance friction
378 double edashpot_; // work dissipated by dashpots
379 double eadhesive_; // work done by attractive force on contacting pieces (positive or negative)
380 };
381
382 // Structure to store dashpot quantities.
383 struct dpProps {
384 dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
385 double dp_nratio_; // normal viscous critical damping ratio
386 double dp_sratio_; // shear viscous critical damping ratio
387 int dp_mode_; // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
388 DVect dp_F_; // Force in the dashpots
389 };
390
391 bool updateKn(const IContactMechanical *con);
392 bool updateKs(const IContactMechanical *con);
393 bool updateFric(const IContactMechanical *con);
394 bool updateResFric(const IContactMechanical *con);
395
396 void updateStiffness(ContactModelMechanicalState *state);
397
398 void setDampCoefficients(const double &mass,double *vcn,double *vcs);
399
400 // Contact model inheritance fields.
401 uint32 inheritanceField_;
402
403 // Effective translational stiffness.
404 DVect2 effectiveTranslationalStiffness_ = DVect2(0.0);
405 DAVect effectiveRotationalStiffness_ = DAVect(0.0); // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
406
407 // linear model properties
408 double kn_ = 0.0; // Normal stiffness
409 double ks_ = 0.0; // Shear stiffness
410 double fric_ = 0.0; // Coulomb friction coefficient
411 DVect lin_F_ = DVect(0.0); // Force carried in the linear model
412 bool lin_S_ = false; // The current slip state
413 uint32 lin_mode_ = 0; // Specifies absolute (0) or incremental (1) calculation mode
414 double rgap_ = 0; // Reference gap
415 dpProps * dpProps_ = nullptr; // The viscous properties
416
417 // rolling resistance properties
418 double res_fric_ = 0.0; // rolling friction coefficient
419 DAVect res_M_ = DAVect(0.0); // moment (bending only)
420 bool res_S_ = false; // The current rolling resistance slip state
421 double kr_ = 0.0; // bending rotational stiffness (read-only, calculated internaly)
422 double fr_ = 0.0; // rolling friction coefficient (rbar*res_fric_) (calculated internaly, not a property)
423
424 // Adhesive properties
425 double a_f0_ = 0.0; // maximum attractive force [force], "a_f0"
426 double a_d0_ = 0.0; // attraction range [length] , "a_d0"
427 double a_F_ = 0.0; // attractive force [force] , "a_force"
428 double userArea_ = 0.0; // Area as specified by the user
429 Energies * energies_ = nullptr; // The energies
430
431 };
432} // namespace cmodelsxd
433// EoF
contactmodelarrlinear.cpp
1// contactmodelARRLinear.cpp
2#include "contactmodelarrlinear.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 "shared/src/mathutil.h"
10
11#include "kernel/interface/iprogram.h"
12#include "module/interface/icontactthermal.h"
13#include "contactmodel/src/contactmodelthermal.h"
14#include "../version.txt"
15#include "fish/src/parameter.h"
16
17#ifdef ARRLINEAR_LIB
18#ifdef _WIN32
19 int __stdcall DllMain(void *,unsigned, void *) {
20 return 1;
21 }
22#endif
23
24 extern "C" EXPORT_TAG const char *getName() {
25#if DIM==3
26 return "contactmodelmechanical3dARRLinear";
27#else
28 return "contactmodelmechanical2dARRLinear";
29#endif
30 }
31
32 extern "C" EXPORT_TAG unsigned getMajorVersion() {
33 return MAJOR_VERSION;
34 }
35
36 extern "C" EXPORT_TAG unsigned getMinorVersion() {
37 return MINOR_VERSION; ;
38 }
39
40 extern "C" EXPORT_TAG void *createInstance() {
41 cmodelsxd::ContactModelARRLinear *m = NEW cmodelsxd::ContactModelARRLinear();
42 return (void *)m;
43 }
44#endif
45
46namespace cmodelsxd {
47 static const uint32 linKnMask = 0x00000002; // Base 1!
48 static const uint32 linKsMask = 0x00000004;
49 static const uint32 linFricMask = 0x00000008;
50 static const uint32 resFricMask = 0x00004000;
51
52 using namespace itasca;
53
54 int ContactModelARRLinear::index_ = -1;
55 uint32 ContactModelARRLinear::getMinorVersion() const { return MINOR_VERSION; }
56
57 ContactModelARRLinear::ContactModelARRLinear() : inheritanceField_(linKnMask|linKsMask|linFricMask|resFricMask)
58 { }
59
60 ContactModelARRLinear::ContactModelARRLinear(const ContactModelARRLinear &m) noexcept {
61 inheritanceField_ = linKnMask|linKsMask|linFricMask|resFricMask;
62 this->copy(&m);
63 }
64
65 const ContactModelARRLinear & ContactModelARRLinear::operator=(const ContactModelARRLinear &m) {
66 inheritanceField_ = linKnMask|linKsMask|linFricMask|resFricMask;
67 this->copy(&m);
68 return *this;
69 }
70
71 void ContactModelARRLinear::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
72 s->addToStorage<ContactModelARRLinear>(*this,ww);
73 }
74
75 ContactModelARRLinear::~ContactModelARRLinear() {
76 // Make sure to clean up after yourself!
77 if (dpProps_)
78 delete dpProps_;
79 if (energies_)
80 delete energies_;
81 }
82
83 void ContactModelARRLinear::archive(ArchiveStream &stream) {
84 // The stream allows one to archive the values of the contact model
85 // so that it can be saved and restored. The minor version can be
86 // used here to allow for incremental changes to the contact model too.
87 stream & kn_;
88 stream & ks_;
89 stream & fric_;
90 stream & lin_F_;
91 stream & lin_S_;
92 stream & lin_mode_;
93 stream & rgap_;
94 stream & res_fric_;
95 stream & res_M_;
96 stream & res_S_;
97 stream & kr_;
98 stream & fr_;
99 stream & a_f0_;
100 stream & a_d0_;
101 stream & a_F_;
102
103 if (stream.getArchiveState()==ArchiveStream::Save) {
104 bool b = false;
105 if (dpProps_) {
106 b = true;
107 stream & b;
108 stream & dpProps_->dp_nratio_;
109 stream & dpProps_->dp_sratio_;
110 stream & dpProps_->dp_mode_;
111 stream & dpProps_->dp_F_;
112 }
113 else
114 stream & b;
115
116 b = false;
117 if (energies_) {
118 b = true;
119 stream & b;
120 stream & energies_->estrain_;
121 stream & energies_->errstrain_;
122 stream & energies_->eslip_;
123 stream & energies_->errslip_;
124 stream & energies_->edashpot_;
125 stream & energies_->eadhesive_;
126 }
127 else
128 stream & b;
129 } else {
130 bool b(false);
131 stream & b;
132 if (b) {
133 if (!dpProps_)
134 dpProps_ = NEW dpProps();
135 stream & dpProps_->dp_nratio_;
136 stream & dpProps_->dp_sratio_;
137 stream & dpProps_->dp_mode_;
138 stream & dpProps_->dp_F_;
139 }
140 stream & b;
141 if (b) {
142 if (!energies_)
143 energies_ = NEW Energies();
144 stream & energies_->estrain_;
145 stream & energies_->errstrain_;
146 stream & energies_->eslip_;
147 stream & energies_->errslip_;
148 stream & energies_->edashpot_;
149 stream & energies_->eadhesive_;
150 }
151 }
152
153 stream & inheritanceField_;
154 stream & effectiveTranslationalStiffness_;
155 stream & effectiveRotationalStiffness_;
156
157 if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() > 1)
158 stream & userArea_;
159 }
160
161 void ContactModelARRLinear::copy(const ContactModel *cm) {
162 // Copy all of the contact model properties. Used in the CMAT
163 // when a new contact is created.
164 ContactModelMechanical::copy(cm);
165 const ContactModelARRLinear *in = dynamic_cast<const ContactModelARRLinear*>(cm);
166 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
167 kn(in->kn());
168 ks(in->ks());
169 fric(in->fric());
170 lin_F(in->lin_F());
171 lin_S(in->lin_S());
172 lin_mode(in->lin_mode());
173 rgap(in->rgap());
174 res_fric(in->res_fric());
175 res_M(in->res_M());
176 res_S(in->res_S());
177 kr(in->kr());
178 fr(in->fr());
179 a_f0(in->a_f0());
180 a_d0(in->a_d0());
181 a_F(in->a_F());
182
183 if (in->hasDamping()) {
184 if (!dpProps_)
185 dpProps_ = NEW dpProps();
186 dp_nratio(in->dp_nratio());
187 dp_sratio(in->dp_sratio());
188 dp_mode(in->dp_mode());
189 dp_F(in->dp_F());
190 }
191 if (in->hasEnergies()) {
192 if (!energies_)
193 energies_ = NEW Energies();
194 estrain(in->estrain());
195 errstrain(in->errstrain());
196 eslip(in->eslip());
197 errslip(in->errslip());
198 edashpot(in->edashpot());
199 eadhesive(in->eadhesive());
200 }
201 userArea_ = in->userArea_;
202 inheritanceField(in->inheritanceField());
203 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
204 effectiveRotationalStiffness(in->effectiveRotationalStiffness());
205 }
206
207
208 base::Property ContactModelARRLinear::getProperty(uint32 i,const IContact *con) const {
209 // Return the property. The IContact pointer is provided so that
210 // more complicated properties, depending on contact characteristics,
211 // can be calcualted. The IContact pointer may be a nullptr!
212 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
213 base::Property var;
214 switch (i) {
215 case kwKn: return kn_;
216 case kwKs: return ks_;
217 case kwFric: return fric_;
218 case kwLinF: var = lin_F_; return var;
219 case kwLinS: return lin_S_;
220 case kwLinMode: return lin_mode_;
221 case kwRGap: return rgap_;
222 case kwEmod: return emod(con);
223 case kwKRatio: return kratio();
224 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
225 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
226 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
227 case kwDpF: {
228 var = dpProps_ ? dpProps_->dp_F_ : DVect(0);
229 return var;
230 }
231 case kwResFric: return res_fric_;
232 case kwResMoment: var = res_M_; return var;
233 case kwResS: return res_S_;
234 case kwResKr: return kr_;
235 case kwAdhesiveF0: return a_f0_;
236 case kwAdhesiveD0: return a_d0_;
237 case kwAdhesiveF: return a_F_;
238 case kwUserArea : return userArea_;
239 }
240 assert(0);
241 return base::Property();
242 }
243
244 bool ContactModelARRLinear::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 kwLinF:
249 case kwDpF:
250 case kwResMoment:
251 return false;
252 }
253 return true;
254 }
255
256 bool ContactModelARRLinear::setProperty(uint32 i,const base::Property &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 kwKn: {
262 if (!v.canConvert<double>())
263 throw Exception("kn must be a double.");
264 double val(v.to<double>());
265 if (val<0.0)
266 throw Exception("Negative kn not allowed.");
267 kn_ = val;
268 return true;
269 }
270 case kwKs: {
271 if (!v.canConvert<double>())
272 throw Exception("ks must be a double.");
273 double val(v.to<double>());
274 if (val<0.0)
275 throw Exception("Negative ks not allowed.");
276 ks_ = 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.to<double>());
283 if (val<0.0)
284 throw Exception("Negative fric not allowed.");
285 fric_ = val;
286 return false;
287 }
288 case kwLinF: {
289 if (!v.canConvert<DVect>())
290 throw Exception("lin_force must be a vector.");
291 DVect val(v.to<DVect>());
292 lin_F_ = val;
293 return false;
294 }
295 case kwLinMode: {
296 if (!v.canConvert<int64>())
297 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
298 uint32 val(v.to<int64>());
299 if (val >1)
300 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
301 lin_mode_ = val;
302 return false;
303 }
304 case kwRGap: {
305 if (!v.canConvert<double>())
306 throw Exception("Reference gap must be a double.");
307 double val(v.to<double>());
308 rgap_ = val;
309 return false;
310 }
311 case kwDpNRatio: {
312 if (!v.canConvert<double>())
313 throw Exception("dp_nratio must be a double.");
314 double val(v.to<double>());
315 if (val<0.0)
316 throw Exception("Negative dp_nratio not allowed.");
317 if (val == 0.0 && !dpProps_)
318 return false;
319 if (!dpProps_)
320 dpProps_ = NEW dpProps();
321 dpProps_->dp_nratio_ = val;
322 return true;
323 }
324 case kwDpSRatio: {
325 if (!v.canConvert<double>())
326 throw Exception("dp_sratio must be a double.");
327 double val(v.to<double>());
328 if (val<0.0)
329 throw Exception("Negative dp_sratio not allowed.");
330 if (val == 0.0 && !dpProps_)
331 return false;
332 if (!dpProps_)
333 dpProps_ = NEW dpProps();
334 dpProps_->dp_sratio_ = val;
335 return true;
336 }
337 case kwDpMode: {
338 if (!v.canConvert<int64>())
339 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
340 int val(v.to<int64>());
341 if (val == 0 && !dpProps_)
342 return false;
343 if (val < 0 || val > 3)
344 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
345 if (!dpProps_)
346 dpProps_ = NEW dpProps();
347 dpProps_->dp_mode_ = val;
348 return false;
349 }
350 case kwDpF: {
351 if (!v.canConvert<DVect>())
352 throw Exception("dp_force must be a vector.");
353 DVect val(v.to<DVect>());
354 if (val.fsum() == 0.0 && !dpProps_)
355 return false;
356 if (!dpProps_)
357 dpProps_ = NEW dpProps();
358 dpProps_->dp_F_ = val;
359 return false;
360 }
361 case kwResFric: {
362 if (!v.canConvert<double>())
363 throw Exception("res_fric must be a double.");
364 double val(v.to<double>());
365 if (val<0.0)
366 throw Exception("Negative res_fric not allowed.");
367 res_fric_ = val;
368 return false;
369 }
370 case kwResMoment: {
371 DAVect val(0.0);
372#ifdef TWOD
373 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
374 throw Exception("res_moment must be an angular vector.");
375 if (v.canConvert<DAVect>())
376 val = DAVect(v.to<DAVect>());
377 else
378 val = DAVect(v.to<double>());
379#else
380 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
381 throw Exception("res_moment must be an angular vector.");
382 if (v.canConvert<DAVect>())
383 val = DAVect(v.value<DAVect>());
384 else
385 val = DAVect(v.value<DVect>());
386#endif
387 res_M_ = val;
388 return false;
389 }
390 case kwAdhesiveF0: {
391 if (!v.canConvert<double>())
392 throw Exception("a_f0 must be a double.");
393 double val(v.to<double>());
394 if (val<0.0)
395 throw Exception("Negative a_f0 not allowed.");
396 a_f0_ = val;
397 return true;
398 }
399 case kwAdhesiveD0: {
400 if (!v.canConvert<double>())
401 throw Exception("a_d0 must be a double.");
402 double val(v.to<double>());
403 if (val<0.0)
404 throw Exception("Negative a_d0 not allowed.");
405 a_d0_ = val;
406 return true;
407 }
408 case kwUserArea: {
409 if (!v.canConvert<double>())
410 throw Exception("user_area must be a double.");
411 double val(v.to<double>());
412 if (val < 0.0)
413 throw Exception("Negative user_area not allowed.");
414 userArea_ = val;
415 return true;
416 }
417 }
418 return false;
419 }
420
421 bool ContactModelARRLinear::getPropertyReadOnly(uint32 i) const {
422 // Returns TRUE if a property is read only or FALSE otherwise.
423 switch (i) {
424 case kwDpF:
425 case kwLinS:
426 case kwEmod:
427 case kwKRatio:
428 case kwResS:
429 case kwResKr:
430 case kwAdhesiveF:
431 return true;
432 default:
433 break;
434 }
435 return false;
436 }
437
438 bool ContactModelARRLinear::supportsInheritance(uint32 i) const {
439 // Returns TRUE if a property supports inheritance or FALSE otherwise.
440 switch (i) {
441 case kwKn:
442 case kwKs:
443 case kwFric:
444 case kwResFric:
445 return true;
446 default:
447 break;
448 }
449 return false;
450 }
451
452 string ContactModelARRLinear::getMethodArguments(uint32 i) const {
453 // Return a list of contact model method argument names.
454 switch (i) {
455 case kwDeformability:
456 return "emod,kratio";
457 case kwArea:
458 return {};
459 }
460 assert(0);
461 return {};
462 }
463
464 bool ContactModelARRLinear::setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con) {
465 FP_S;
466 // Apply the specified method.
467 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
468 switch (i) {
469 case kwDeformability: {
470 double emod;
471 double krat;
472 if (vl.at(0).isNull())
473 throw Exception("Argument emod must be specified with method deformability in contact model {0}.",getName());
474 emod = vl.at(0).to<double>();
475 if (emod<0.0)
476 throw Exception("Negative emod not allowed in contact model {0}.",getName());
477 if (vl.at(1).isNull())
478 throw Exception("Argument kratio must be specified with method deformability in contact model {0}.",getName());
479 krat = vl.at(1).to<double>();
480 if (krat<0.0)
481 throw Exception("Negative stiffness ratio not allowed in contact model {0}.",getName());
482 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
483 double rsum = calcRSum(c);
484 if (userArea_) {
485#ifdef THREED
486 rsq = std::sqrt(userArea_ / dPi);
487#else
488 rsq = userArea_ / 2.0;
489#endif
490 rsum = rsq + rsq;
491 rsq = safeDiv(1. , rsq);
492 }
493#ifdef TWOD
494 kn_ = safeDiv(2.0 * emod ,rsq * rsum);
495#else
496 kn_ = safeDiv(dPi * emod ,rsq * rsq * rsum);
497#endif
498 ks_ = safeDiv(kn_ , krat);
499 setInheritance(1,false);
500 setInheritance(2,false);
501 FP_S;
502 return true;
503 }
504 case kwArea: {
505 if (!userArea_) {
506 double rsq = calcRSQ(c);
507#ifdef THREED
508 userArea_ = rsq * rsq * dPi;
509#else
510 userArea_ = rsq * 2.0;
511#endif
512 }
513 FP_S;
514 return true;
515 }
516 }
517 FP_S;
518 return false;
519 }
520
521 double ContactModelARRLinear::getEnergy(uint32 i) const {
522 // Return an energy value.
523 double ret(0.0);
524 if (!energies_)
525 return ret;
526 switch (i) {
527 case kwEStrain: return energies_->estrain_;
528 case kwERRStrain: return energies_->errstrain_;
529 case kwESlip: return energies_->eslip_;
530 case kwERRSlip: return energies_->errslip_;
531 case kwEDashpot: return energies_->edashpot_;
532 case kwEAdhesive: return energies_->eadhesive_;
533 }
534 assert(0);
535 return ret;
536 }
537
538 bool ContactModelARRLinear::getEnergyAccumulate(uint32 i) const {
539 // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
540 switch (i) {
541 case kwEStrain: return false;
542 case kwERRStrain: return false;
543 case kwESlip: return true;
544 case kwERRSlip: return true;
545 case kwEDashpot: return true;
546 case kwEAdhesive: return true;
547 }
548 assert(0);
549 return false;
550 }
551
552 void ContactModelARRLinear::setEnergy(uint32 i,const double &d) {
553 // Set an energy value.
554 if (!energies_) return;
555 switch (i) {
556 case kwEStrain: energies_->estrain_ = d; return;
557 case kwERRStrain: energies_->errstrain_ = d; return;
558 case kwESlip: energies_->eslip_ = d; return;
559 case kwERRSlip: energies_->errslip_ = d; return;
560 case kwEDashpot: energies_->edashpot_= d; return;
561 case kwEAdhesive: energies_->eadhesive_= d; return;
562 }
563 assert(0);
564 return;
565 }
566
567 bool ContactModelARRLinear::validate(ContactModelMechanicalState *state,const double &) {
568 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
569 // (1) the contact is created, (2) a property of the contact that returns a true via
570 // the setProperty method has been modified and (3) when a set of cycles is executed
571 // via the {cycle N} command.
572 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
573 FP_S;
574 assert(state);
575 const IContactMechanical *c = state->getMechanicalContact();
576 assert(c);
577 FP_S;
578
579 if (state->trackEnergy_)
580 activateEnergy();
581 FP_S;
582
583 if (inheritanceField_ & linKnMask)
584 updateKn(c);
585 FP_S;
586 if (inheritanceField_ & linKsMask)
587 updateKs(c);
588 FP_S;
589 if (inheritanceField_ & linFricMask)
590 updateFric(c);
591 FP_S;
592 if (inheritanceField_ & resFricMask)
593 updateResFric(c);
594 FP_S;
595
596 updateStiffness(state);
597 FP_S;
598 auto ret = checkActivity(state->gap_);
599 FP_S;
600 return ret;
601 }
602
603 static const string knstr("kn");
604 bool ContactModelARRLinear::updateKn(const IContactMechanical *con) {
605 assert(con);
606 base::Property v1 = con->getEnd1()->getProperty(knstr);
607 base::Property v2 = con->getEnd2()->getProperty(knstr);
608 if (v1.isNull() || v2.isNull())
609 return false;
610 double kn1 = v1.to<double>();
611 double kn2 = v2.to<double>();
612 double val = kn_;
613 if (kn1 && kn2)
614 kn_ = kn1*kn2/(kn1+kn2);
615 else if (kn1)
616 kn_ = kn1;
617 else if (kn2)
618 kn_ = kn2;
619 return ( (kn_ != val) );
620 }
621
622 static const string ksstr("ks");
623 bool ContactModelARRLinear::updateKs(const IContactMechanical *con) {
624 assert(con);
625 base::Property v1 = con->getEnd1()->getProperty(ksstr);
626 base::Property v2 = con->getEnd2()->getProperty(ksstr);
627 if (v1.isNull() || v2.isNull())
628 return false;
629 double ks1 = v1.to<double>();
630 double ks2 = v2.to<double>();
631 double val = ks_;
632 if (ks1 && ks2)
633 ks_ = ks1*ks2/(ks1+ks2);
634 else if (ks1)
635 ks_ = ks1;
636 else if (ks2)
637 ks_ = ks2;
638 return ( (ks_ != val) );
639 }
640
641 static const string fricstr("fric");
642 bool ContactModelARRLinear::updateFric(const IContactMechanical *con) {
643 assert(con);
644 base::Property v1 = con->getEnd1()->getProperty(fricstr);
645 base::Property v2 = con->getEnd2()->getProperty(fricstr);
646 if (v1.isNull() || v2.isNull())
647 return false;
648 double fric1 = std::max(0.0,v1.to<double>());
649 double fric2 = std::max(0.0,v2.to<double>());
650 double val = fric_;
651 fric_ = std::min(fric1,fric2);
652 return ( (fric_ != val) );
653 }
654
655 static const string rfricstr("rr_fric");
656 bool ContactModelARRLinear::updateResFric(const IContactMechanical *con) {
657 assert(con);
658 base::Property v1 = con->getEnd1()->getProperty(rfricstr);
659 base::Property v2 = con->getEnd2()->getProperty(rfricstr);
660 if (v1.isNull() || v2.isNull())
661 return false;
662 double rfric1 = std::max(0.0,v1.to<double>());
663 double rfric2 = std::max(0.0,v2.to<double>());
664 double val = res_fric_;
665 res_fric_ = std::min(rfric1,rfric2);
666 return ( (res_fric_ != val) );
667 }
668
669 bool ContactModelARRLinear::endPropertyUpdated(const string &name,const IContactMechanical *c) {
670 // The endPropertyUpdated method is called whenever a surface property (with a name
671 // that matches an inheritable contact model property name) of one of the contacting
672 // pieces is modified. This allows the contact model to update its associated
673 // properties. The return value denotes whether or not the update has affected
674 // the time step computation (by having modified the translational or rotational
675 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
676 assert(c);
677 StringList availableProperties = split(replace(simplified(getProperties())," ",""),",");
678 auto idx = findRegex(availableProperties,name);
679 if (idx==string::npos) return false;
680 idx += 1;
681 bool ret = false;
682 switch(idx) {
683 case kwKn: { //kn
684 if (inheritanceField_ & linKnMask)
685 ret = updateKn(c);
686 break;
687 }
688 case kwKs: { //ks
689 if (inheritanceField_ & linKsMask)
690 ret =updateKs(c);
691 break;
692 }
693 case kwFric: { //fric
694 if (inheritanceField_ & linFricMask)
695 updateFric(c);
696 break;
697 }
698 case kwResFric: { //rr_fric
699 if (inheritanceField_ & resFricMask)
700 ret = updateResFric(c);
701 break;
702 }
703 }
704 return ret;
705 }
706
707 void ContactModelARRLinear::updateStiffness(ContactModelMechanicalState *state) {
708 // first compute rolling resistance stiffness
709 kr_ = 0.0;
710 if (res_fric_ > 0.0) {
711 double rbar = 0.0;
712 double r1 = 1.0/state->end1Curvature_.y();
713 rbar = r1;
714 double r2 = 0.0;
715 if (state->end2Curvature_.y()) {
716 r2 = 1.0 / state->end2Curvature_.y();
717 rbar = (r1*r2) / (r1+r2);
718 }
719 if (userArea_) {
720#ifdef THREED
721 r1 = std::sqrt(userArea_ / dPi);
722#else
723 r1 = userArea_ / 2.0;
724#endif
725 r2 = r1;
726 rbar = (r1*r2) / (r1+r2);
727 }
728 kr_ = ks_*rbar*rbar;
729 fr_ = res_fric_*rbar;
730 }
731 // Now calculate effective stiffness
732 DVect2 retT(kn_,ks_);
733 // correction if viscous damping active
734 if (dpProps_) {
735 DVect2 correct(1.0);
736 if (dpProps_->dp_nratio_)
737 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
738 if (dpProps_->dp_sratio_)
739 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
740 retT /= (correct*correct);
741 }
742 // Correction for adhesive group
743 if (a_d0_ != 0.0) { retT.rdof(0) += a_f0_ / a_d0_; }
744
745 effectiveTranslationalStiffness_ = retT;
746 // Effective rotational stiffness (bending only)
747 effectiveRotationalStiffness_ = DAVect(kr_);
748#if DIM==3
749 effectiveRotationalStiffness_.rx() = 0.0;
750#endif
751 }
752
753 bool ContactModelARRLinear::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
754 assert(state);
755
756 // Current overlap
757 double overlap = rgap_ - state->gap_;
758 // Relative translational increment
759 DVect trans = state->relativeTranslationalIncrement_;
760 // Correction factor to account for when the contact becomes newly active.
761 // We estimate the time of activity during the timestep when the contact has first
762 // become active and scale the forces accordingly.
763 double correction = 1.0;
764
765 // The contact was just activated from an inactive state
766 if (state->activated()) {
767 // Trigger the FISH callback if one is hooked up to the
768 // contact_activated event.
769 if (cmEvents_[fActivated] >= 0) {
770 auto c = state->getContact();
771 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
772 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
773 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
774 }
775 // Calculate the correction factor.
776 if (trans.x()) {
777 correction = -1.0*overlap / trans.x();
778 if (correction < 0)
779 correction = 1.0;
780 }
781 }
782
783 // Angular dispacement increment.
784 DAVect ang = state->relativeAngularIncrement_;
785 DVect lin_F_old = lin_F_;
786
787 if (lin_mode_ == 0)
788 lin_F_.rx() = overlap * kn_; // absolute mode for normal force calculation
789 else
790 lin_F_.rx() -= correction * trans.x() * kn_; // incremental mode for normal force calculation
791
792 // Normal force can only be positive.
793 lin_F_.rx() = std::max(0.0,lin_F_.x());
794
795 // Calculate the shear force.
796 DVect sforce(0.0);
797 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
798 // Loop over the shear components (note: the 0 component is the normal component)
799 // and calculate the shear force.
800 for (int i=1; i<dim; ++i)
801 sforce.rdof(i) = lin_F_.dof(i) - trans.dof(i) * ks_ * correction;
802
803 // The canFail flag corresponds to whether or not the contact can undergo non-linear
804 // force-displacement response. If the SOLVE ELASTIC command is given then the
805 // canFail state is set to FALSE. Otherwise it is always TRUE.
806 if (state->canFail_) {
807 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
808 double crit = lin_F_.x() * fric_;
809 // The is the magnitude of the shear force.
810 double sfmag = sforce.mag();
811 // Sliding occurs when the magnitude of the shear force is greater than the
812 // critical value.
813 if (sfmag > crit) {
814 // Lower the shear force to the critical value for sliding.
815 double rat = crit / sfmag;
816 sforce *= rat;
817 // Handle the slip_change event if one has been hooked up. Sliding has commenced.
818 if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
819 auto c = state->getContact();
820 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
821 fish::Parameter() };
822 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
823 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
824 }
825 lin_S_ = true;
826 } else {
827 // Handle the slip_change event if one has been hooked up and
828 // the contact was previously sliding. Sliding has ceased.
829 if (lin_S_) {
830 if (cmEvents_[fSlipChange] >= 0) {
831 auto c = state->getContact();
832 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
833 fish::Parameter((int64)1) };
834 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
835 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
836 }
837 lin_S_ = false;
838 }
839 }
840 }
841
842 // Set the shear components of the total force.
843 for (int i=1; i<dim; ++i)
844 lin_F_.rdof(i) = sforce.dof(i);
845
846 // Rolling resistance
847 DAVect res_M_old = res_M_;
848 if ((fr_ == 0.0) || (kr_==0.0)) {
849 res_M_.fill(0.0);
850 } else {
851 DAVect angStiff(0.0);
852 DAVect MomentInc(0.0);
853#if DIM==3
854 angStiff.rx() = 0.0;
855 angStiff.ry() = kr_;
856#endif
857 angStiff.rz() = kr_;
858 MomentInc = ang * angStiff * (-1.0);
859 res_M_ += MomentInc;
860 if (state->canFail_) {
861 // Account for bending strength
862 double rmax = std::abs(fr_*lin_F_.x());
863 double rmag = res_M_.mag();
864 if (rmag > rmax) {
865 double fac = rmax/rmag;
866 res_M_ *= fac;
867 res_S_ = true;
868 } else {
869 res_S_ = false;
870 }
871 }
872 }
873
874 // Account for dashpot forces if the dashpot structure has been defined.
875 if (dpProps_) {
876 dpProps_->dp_F_.fill(0.0);
877 double vcn(0.0), vcs(0.0);
878 // Calculate the damping coefficients.
879 setDampCoefficients(state->inertialMass_,&vcn,&vcs);
880 // First damp the shear components
881 for (int i=1; i<dim; ++i)
882 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
883 // Damp the normal component
884 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
885 // Need to change behavior based on the dp_mode.
886 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
887 // Limit in tension if not bonded.
888 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
889 dpProps_->dp_F_.rx() = - lin_F_.rx();
890 }
891 if (lin_S_ && dpProps_->dp_mode_ > 1) {
892 // Limit in shear if not sliding.
893 double dfn = dpProps_->dp_F_.rx();
894 dpProps_->dp_F_.fill(0.0);
895 dpProps_->dp_F_.rx() = dfn;
896 }
897 }
898
899 // Adhesive force
900 double a_F_old = a_F_;
901 double gs = state->gap_ - rgap_;
902 a_F_ = 0.0;
903 if ( gs <= 0.0 )
904 a_F_ = a_f0_;
905 else if ( gs < a_d0_ ) // if a_d0_ == 0.0, will not enter, so next line divide by a_d0_ is ok
906 a_F_ = a_f0_ * ( 1.0 - (gs/a_d0_) );
907
908 //Compute energies if energy tracking has been enabled.
909 if (state->trackEnergy_) {
910 assert(energies_);
911 energies_->estrain_ = 0.0;
912 if (kn_)
913 // Calculate the strain energy.
914 energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
915 if (ks_) {
916 DVect s = lin_F_;
917 s.rx() = 0.0;
918 double smag2 = s.mag2();
919 // Add the shear component of the strain energy.
920 energies_->estrain_ += 0.5*smag2 / ks_;
921
922 if (lin_S_) {
923 // If sliding calculate the slip energy and accumulate it.
924 lin_F_old.rx() = 0.0;
925 DVect avg_F_s = (s + lin_F_old)*0.5;
926 DVect u_s_el = (s - lin_F_old) / ks_;
927 DVect u_s(0.0);
928 for (int i=1; i<dim; ++i)
929 u_s.rdof(i) = trans.dof(i);
930 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
931 }
932 }
933 // Add the rolling resistance energy contributions.
934 energies_->errstrain_ = 0.0;
935 if (kr_) {
936 energies_->errstrain_ = 0.5*res_M_.mag2() / kr_;
937 if (res_S_) {
938 // If sliding calculate the slip energy and accumulate it.
939 DAVect avg_M = (res_M_ + res_M_old)*0.5;
940 DAVect t_s_el = (res_M_ - res_M_old) / kr_;
941 energies_->errslip_ -= std::min(0.0,(avg_M | (ang + t_s_el)));
942 }
943 }
944 // Add the adhesive energy contribution:
945 energies_->eadhesive_ -= 0.5*(a_F_old + a_F_) * trans.x();
946
947 if (dpProps_) {
948 // Calculate damping energy (accumulated) if the dashpots are active.
949 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
950 }
951 }
952
953 // This is just a sanity check to ensure, in debug mode, that the force isn't wonky.
954 assert(lin_F_ == lin_F_);
955 return true;
956 }
957
958 bool ContactModelARRLinear::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
959 // Account for thermal expansion in incremental mode
960 if (lin_mode_ == 0 || ts->gapInc_ == 0.0) return false;
961 DVect finc(0.0);
962 finc.rx() = kn_ * ts->gapInc_;
963 lin_F_ -= finc;
964 return true;
965 }
966
967 void ContactModelARRLinear::setForce(const DVect &v,IContact *c) {
968 lin_F(v);
969 if (v.x() > 0)
970 rgap_ = c->getGap() + v.x() / kn_;
971 }
972
973 void ContactModelARRLinear::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
974 // Only called for contacts with wall facets when the wall resolution scheme
975 // is set to full!
976 // Only do something if the contact model is of the same type
977 if (equal(old->getContactModel()->getName(),"arrlinear") && !isBonded()) {
978 ContactModelARRLinear *oldCm = (ContactModelARRLinear *)old;
979#ifdef THREED
980 // Need to rotate just the shear component from oldSystem to newSystem
981
982 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
983 DVect axis = oldSystem.e1() & newSystem.e1();
984 double c, ang, s;
985 DVect re2;
986 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
987 axis = axis.unit();
988 c = oldSystem.e1()|newSystem.e1();
989 if (c > 0)
990 c = std::min(c,1.0);
991 else
992 c = std::max(c,-1.0);
993 ang = acos(c);
994 s = sin(ang);
995 double t = 1. - c;
996 DMatrix<3,3> rm;
997 rm.get(0,0) = t*axis.x()*axis.x() + c;
998 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
999 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
1000 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
1001 rm.get(1,1) = t*axis.y()*axis.y() + c;
1002 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
1003 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
1004 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
1005 rm.get(2,2) = t*axis.z()*axis.z() + c;
1006 re2 = rm*oldSystem.e2();
1007 }
1008 else
1009 re2 = oldSystem.e2();
1010 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1011 axis = re2 & newSystem.e2();
1012 DVect2 tpf;
1013 DVect2 tpm;
1014 DMatrix<2,2> m;
1015 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1016 axis = axis.unit();
1017 c = re2|newSystem.e2();
1018 if (c > 0)
1019 c = std::min(c,1.0);
1020 else
1021 c = std::max(c,-1.0);
1022 ang = acos(c);
1023 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1024 ang *= -1;
1025 s = sin(ang);
1026 m.get(0,0) = c;
1027 m.get(1,0) = s;
1028 m.get(0,1) = -m.get(1,0);
1029 m.get(1,1) = m.get(0,0);
1030 tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1031 tpm = m*DVect2(oldCm->res_M_.y(),oldCm->res_M_.z());
1032 } else {
1033 m.get(0,0) = 1.;
1034 m.get(0,1) = 0.;
1035 m.get(1,0) = 0.;
1036 m.get(1,1) = 1.;
1037 tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1038 tpm = DVect2(oldCm->res_M_.y(),oldCm->res_M_.z());
1039 }
1040 DVect pforce = DVect(0,tpf.x(),tpf.y());
1041 //DVect pm = DVect(0,tpm.x(),tpm.y());
1042#else
1043 oldSystem;
1044 newSystem;
1045 DVect pforce = DVect(0,oldCm->lin_F_.y());
1046 //DVect pm = DVect(0,oldCm->res_M_.y());
1047#endif
1048 for (int i=1; i<dim; ++i)
1049 lin_F_.rdof(i) += pforce.dof(i);
1050 if (lin_mode_ && oldCm->lin_mode_)
1051 lin_F_.rx() = oldCm->lin_F_.x();
1052 oldCm->lin_F_ = DVect(0.0);
1053 oldCm->res_M_ = DAVect(0.0);
1054 if (dpProps_ && oldCm->dpProps_) {
1055#ifdef THREED
1056 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1057 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1058#else
1059 pforce = oldCm->dpProps_->dp_F_;
1060#endif
1061 dpProps_->dp_F_ += pforce;
1062 oldCm->dpProps_->dp_F_ = DVect(0.0);
1063 }
1064 if(oldCm->getEnergyActivated()) {
1065 activateEnergy();
1066 energies_->estrain_ = oldCm->energies_->estrain_;
1067 energies_->edashpot_ = oldCm->energies_->edashpot_;
1068 energies_->eslip_ = oldCm->energies_->eslip_;
1069 oldCm->energies_->estrain_ = 0.0;
1070 oldCm->energies_->edashpot_ = 0.0;
1071 oldCm->energies_->eslip_ = 0.0;
1072 }
1073 }
1074 assert(lin_F_ == lin_F_);
1075 }
1076
1077 void ContactModelARRLinear::setNonForcePropsFrom(IContactModel *old) {
1078 // Only called for contacts with wall facets when the wall resolution scheme
1079 // is set to full!
1080 // Only do something if the contact model is of the same type
1081 if (equal(old->getName(),"arrlinear") && !isBonded()) {
1082 ContactModelARRLinear *oldCm = (ContactModelARRLinear *)old;
1083 kn_ = oldCm->kn_;
1084 ks_ = oldCm->ks_;
1085 fric_ = oldCm->fric_;
1086 lin_mode_ = oldCm->lin_mode_;
1087 rgap_ = oldCm->rgap_;
1088 res_fric_ = oldCm->res_fric_;
1089 res_S_ = oldCm->res_S_;
1090 kr_ = oldCm->kr_;
1091 fr_ = oldCm->fr_;
1092 a_f0_ = oldCm->a_f0_;
1093 a_d0_ = oldCm->a_d0_;
1094 userArea_ = oldCm->userArea_;
1095
1096 if (oldCm->dpProps_) {
1097 if (!dpProps_)
1098 dpProps_ = NEW dpProps();
1099 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1100 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1101 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1102 }
1103 }
1104 }
1105
1106 DVect ContactModelARRLinear::getForce() const {
1107 DVect ret(lin_F_);
1108 if (dpProps_)
1109 ret += dpProps_->dp_F_;
1110 ret.rdof(0) -= a_F_;
1111 return ret;
1112 }
1113
1114 DAVect ContactModelARRLinear::getMomentOn1(const IContactMechanical *c) const {
1115 DVect force = getForce();
1116 DAVect ret(res_M_);
1117 c->updateResultingTorqueOn1Local(force,&ret);
1118 return ret;
1119 }
1120
1121 DAVect ContactModelARRLinear::getMomentOn2(const IContactMechanical *c) const {
1122 DVect force = getForce();
1123 DAVect ret(res_M_);
1124 c->updateResultingTorqueOn2Local(force,&ret);
1125 return ret;
1126 }
1127
1128 DAVect ContactModelARRLinear::getModelMomentOn1() const {
1129 DAVect ret(res_M_);
1130 return ret;
1131 }
1132
1133 DAVect ContactModelARRLinear::getModelMomentOn2() const {
1134 DAVect ret(res_M_);
1135 return ret;
1136 }
1137
1138 void ContactModelARRLinear::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
1139 ret->clear();
1140 ret->push_back({"kn",ScalarInfo});
1141 ret->push_back({"ks",ScalarInfo});
1142 ret->push_back({"fric",ScalarInfo});
1143 ret->push_back({"lin_force",VectorInfo});
1144 ret->push_back({"lin_slip",ScalarInfo});
1145 ret->push_back({"lin_mode",ScalarInfo});
1146 ret->push_back({"rgap",ScalarInfo});
1147 ret->push_back({"emod",ScalarInfo});
1148 ret->push_back({"kratio",ScalarInfo});
1149 ret->push_back({"dp_nratio",ScalarInfo});
1150 ret->push_back({"dp_sratio",ScalarInfo});
1151 ret->push_back({"dp_mode",ScalarInfo});
1152 ret->push_back({"dp_force",VectorInfo});
1153 ret->push_back({"rr_fric",ScalarInfo});
1154 ret->push_back({"rr_moment",VectorInfo});
1155 ret->push_back({"rr_slip",ScalarInfo});
1156 ret->push_back({"rr_kr",ScalarInfo});
1157 ret->push_back({"adh_f0",ScalarInfo});
1158 ret->push_back({"adh_d0",ScalarInfo});
1159 ret->push_back({"adh_force",VectorInfo});
1160 ret->push_back({"user_area",ScalarInfo});
1161 }
1162
1163 void ContactModelARRLinear::objectPropValues(std::vector<double> *ret,const IContact *c) const {
1164 FP_S;
1165 ret->clear();
1166 ret->push_back(kn());
1167 ret->push_back(ks());
1168 ret->push_back(fric());
1169 ret->push_back(safeMag(lin_F()));
1170 ret->push_back(lin_S());
1171 ret->push_back(lin_mode());
1172 ret->push_back(rgap());
1173 ret->push_back(emod(c));
1174 ret->push_back(kratio());
1175 ret->push_back(dp_nratio());
1176 ret->push_back(dp_sratio());
1177 ret->push_back(dp_mode());
1178 ret->push_back(safeMag(dp_F()));
1179 ret->push_back(res_fric());
1180 ret->push_back(safeMag(res_M()));
1181 ret->push_back(res_S());
1182 ret->push_back(kr());
1183 ret->push_back(a_f0());
1184 ret->push_back(a_d0());
1185 ret->push_back(a_F());
1186 ret->push_back(getArea());
1187 FP_S;
1188 }
1189
1190 double ContactModelARRLinear::emod(const IContact *con) const {
1191 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
1192 if (c ==nullptr) return 0.0;
1193 double rsq = calcRSQ(c);
1194 double rsum = calcRSum(c);
1195 if (userArea_) {
1196#ifdef THREED
1197 rsq = std::sqrt(userArea_ / dPi);
1198#else
1199 rsq = userArea_ / 2.0;
1200#endif
1201 rsum = rsq + rsq;
1202 rsq = safeDiv(1. , rsq);
1203 }
1204#ifdef TWOD
1205 return (kn_ * rsum * rsq / 2.0);
1206#else
1207 return (kn_ * rsum * rsq * rsq) / dPi;
1208#endif
1209 }
1210
1211 double ContactModelARRLinear::kratio() const {
1212 return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
1213 }
1214
1215 void ContactModelARRLinear::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
1216 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
1217 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
1218 }
1219
1220} // namespace cmodelsxd
1221// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
