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 QString 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 QString 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 QString 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_ = NEWC(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 QString getFishCallEvents() const override {
132 return
133 "contact_activated"
134 ",slip_change";
135 }
136
137 // Return the specified contact model property.
138 QVariant 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 QVariant &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 QString getMethods() const override { return "deformability,area";}
171 // Return a comma seprated list of the contact model method arguments (base 1).
172 QString 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 QVector<QVariant> &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 QString &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 NEWC(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<QString,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_ = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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 QVariant 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 QVariant var;
214 switch (i) {
215 case kwKn: return kn_;
216 case kwKs: return ks_;
217 case kwFric: return fric_;
218 case kwLinF: var.setValue(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 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
229 return var;
230 }
231 case kwResFric: return res_fric_;
232 case kwResMoment: var.setValue(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 QVariant();
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 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 kwKn: {
262 if (!v.canConvert<double>())
263 throw Exception("kn must be a double.");
264 double val(v.toDouble());
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.toDouble());
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.toDouble());
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.value<DVect>());
292 lin_F_ = val;
293 return false;
294 }
295 case kwLinMode: {
296 if (!v.canConvert<uint32>())
297 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
298 uint32 val(v.toUInt());
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.toDouble());
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.toDouble());
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_ = NEWC(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.toDouble());
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_ = NEWC(dpProps());
334 dpProps_->dp_sratio_ = val;
335 return true;
336 }
337 case kwDpMode: {
338 if (!v.canConvert<int>())
339 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
340 int val(v.toInt());
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_ = NEWC(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.value<DVect>());
354 if (val.fsum() == 0.0 && !dpProps_)
355 return false;
356 if (!dpProps_)
357 dpProps_ = NEWC(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.toDouble());
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.value<DAVect>());
377 else
378 val = DAVect(v.toDouble());
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.toDouble());
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.toDouble());
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.toDouble());
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 QString 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 QString();
459 }
460 assert(0);
461 return QString();
462 }
463
464 bool ContactModelARRLinear::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
465 // Apply the specified method.
466 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
467 switch (i) {
468 case kwDeformability: {
469 double emod;
470 double krat;
471 if (vl.at(0).isNull())
472 throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
473 emod = vl.at(0).toDouble();
474 if (emod<0.0)
475 throw Exception("Negative emod not allowed in contact model %1.",getName());
476 if (vl.at(1).isNull())
477 throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
478 krat = vl.at(1).toDouble();
479 if (krat<0.0)
480 throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
481 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
482 double rsum(0.0);
483 if (c->getEnd1Curvature().y())
484 rsum += 1.0/c->getEnd1Curvature().y();
485 if (c->getEnd2Curvature().y())
486 rsum += 1.0/c->getEnd2Curvature().y();
487 if (userArea_) {
488#ifdef THREED
489 rsq = std::sqrt(userArea_ / dPi);
490#else
491 rsq = userArea_ / 2.0;
492#endif
493 rsum = rsq + rsq;
494 rsq = 1. / rsq;
495 }
496#ifdef TWOD
497 kn_ = 2.0 * emod / (rsq * rsum);
498#else
499 kn_ = dPi * emod / (rsq * rsq * rsum);
500#endif
501 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
502 setInheritance(1,false);
503 setInheritance(2,false);
504 return true;
505 }
506 case kwArea: {
507 if (!userArea_) {
508 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
509#ifdef THREED
510 userArea_ = rsq * rsq * dPi;
511#else
512 userArea_ = rsq * 2.0;
513#endif
514 }
515 return true;
516 }
517 }
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 assert(state);
574 const IContactMechanical *c = state->getMechanicalContact();
575 assert(c);
576
577 if (state->trackEnergy_)
578 activateEnergy();
579
580 if (inheritanceField_ & linKnMask)
581 updateKn(c);
582 if (inheritanceField_ & linKsMask)
583 updateKs(c);
584 if (inheritanceField_ & linFricMask)
585 updateFric(c);
586 if (inheritanceField_ & resFricMask)
587 updateResFric(c);
588
589 updateStiffness(state);
590 return checkActivity(state->gap_);
591 }
592
593 static const QString knstr("kn");
594 bool ContactModelARRLinear::updateKn(const IContactMechanical *con) {
595 assert(con);
596 QVariant v1 = con->getEnd1()->getProperty(knstr);
597 QVariant v2 = con->getEnd2()->getProperty(knstr);
598 if (!v1.isValid() || !v2.isValid())
599 return false;
600 double kn1 = v1.toDouble();
601 double kn2 = v2.toDouble();
602 double val = kn_;
603 if (kn1 && kn2)
604 kn_ = kn1*kn2/(kn1+kn2);
605 else if (kn1)
606 kn_ = kn1;
607 else if (kn2)
608 kn_ = kn2;
609 return ( (kn_ != val) );
610 }
611
612 static const QString ksstr("ks");
613 bool ContactModelARRLinear::updateKs(const IContactMechanical *con) {
614 assert(con);
615 QVariant v1 = con->getEnd1()->getProperty(ksstr);
616 QVariant v2 = con->getEnd2()->getProperty(ksstr);
617 if (!v1.isValid() || !v2.isValid())
618 return false;
619 double ks1 = v1.toDouble();
620 double ks2 = v2.toDouble();
621 double val = ks_;
622 if (ks1 && ks2)
623 ks_ = ks1*ks2/(ks1+ks2);
624 else if (ks1)
625 ks_ = ks1;
626 else if (ks2)
627 ks_ = ks2;
628 return ( (ks_ != val) );
629 }
630
631 static const QString fricstr("fric");
632 bool ContactModelARRLinear::updateFric(const IContactMechanical *con) {
633 assert(con);
634 QVariant v1 = con->getEnd1()->getProperty(fricstr);
635 QVariant v2 = con->getEnd2()->getProperty(fricstr);
636 if (!v1.isValid() || !v2.isValid())
637 return false;
638 double fric1 = std::max(0.0,v1.toDouble());
639 double fric2 = std::max(0.0,v2.toDouble());
640 double val = fric_;
641 fric_ = std::min(fric1,fric2);
642 return ( (fric_ != val) );
643 }
644
645 static const QString rfricstr("rr_fric");
646 bool ContactModelARRLinear::updateResFric(const IContactMechanical *con) {
647 assert(con);
648 QVariant v1 = con->getEnd1()->getProperty(rfricstr);
649 QVariant v2 = con->getEnd2()->getProperty(rfricstr);
650 if (!v1.isValid() || !v2.isValid())
651 return false;
652 double rfric1 = std::max(0.0,v1.toDouble());
653 double rfric2 = std::max(0.0,v2.toDouble());
654 double val = res_fric_;
655 res_fric_ = std::min(rfric1,rfric2);
656 return ( (res_fric_ != val) );
657 }
658
659 bool ContactModelARRLinear::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
660 // The endPropertyUpdated method is called whenever a surface property (with a name
661 // that matches an inheritable contact model property name) of one of the contacting
662 // pieces is modified. This allows the contact model to update its associated
663 // properties. The return value denotes whether or not the update has affected
664 // the time step computation (by having modified the translational or rotational
665 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
666 assert(c);
667 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
668 QRegExp rx(name,Qt::CaseInsensitive);
669 int idx = availableProperties.indexOf(rx)+1;
670 bool ret=false;
671
672 if (idx<=0)
673 return ret;
674
675 switch(idx) {
676 case kwKn: { //kn
677 if (inheritanceField_ & linKnMask)
678 ret = updateKn(c);
679 break;
680 }
681 case kwKs: { //ks
682 if (inheritanceField_ & linKsMask)
683 ret =updateKs(c);
684 break;
685 }
686 case kwFric: { //fric
687 if (inheritanceField_ & linFricMask)
688 updateFric(c);
689 break;
690 }
691 case kwResFric: { //rr_fric
692 if (inheritanceField_ & resFricMask)
693 ret = updateResFric(c);
694 break;
695 }
696 }
697 return ret;
698 }
699
700 void ContactModelARRLinear::updateStiffness(ContactModelMechanicalState *state) {
701 // first compute rolling resistance stiffness
702 kr_ = 0.0;
703 if (res_fric_ > 0.0) {
704 double rbar = 0.0;
705 double r1 = 1.0/state->end1Curvature_.y();
706 rbar = r1;
707 double r2 = 0.0;
708 if (state->end2Curvature_.y()) {
709 r2 = 1.0 / state->end2Curvature_.y();
710 rbar = (r1*r2) / (r1+r2);
711 }
712 if (userArea_) {
713#ifdef THREED
714 r1 = std::sqrt(userArea_ / dPi);
715#else
716 r1 = userArea_ / 2.0;
717#endif
718 r2 = r1;
719 rbar = (r1*r2) / (r1+r2);
720 }
721 kr_ = ks_*rbar*rbar;
722 fr_ = res_fric_*rbar;
723 }
724 // Now calculate effective stiffness
725 DVect2 retT(kn_,ks_);
726 // correction if viscous damping active
727 if (dpProps_) {
728 DVect2 correct(1.0);
729 if (dpProps_->dp_nratio_)
730 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
731 if (dpProps_->dp_sratio_)
732 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
733 retT /= (correct*correct);
734 }
735 // Correction for adhesive group
736 if (a_d0_ != 0.0) { retT.rdof(0) += a_f0_ / a_d0_; }
737
738 effectiveTranslationalStiffness_ = retT;
739 // Effective rotational stiffness (bending only)
740 effectiveRotationalStiffness_ = DAVect(kr_);
741#if DIM==3
742 effectiveRotationalStiffness_.rx() = 0.0;
743#endif
744 }
745
746 bool ContactModelARRLinear::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
747 assert(state);
748
749 // Current overlap
750 double overlap = rgap_ - state->gap_;
751 // Relative translational increment
752 DVect trans = state->relativeTranslationalIncrement_;
753 // Correction factor to account for when the contact becomes newly active.
754 // We estimate the time of activity during the timestep when the contact has first
755 // become active and scale the forces accordingly.
756 double correction = 1.0;
757
758 // The contact was just activated from an inactive state
759 if (state->activated()) {
760 // Trigger the FISH callback if one is hooked up to the
761 // contact_activated event.
762 if (cmEvents_[fActivated] >= 0) {
763 auto c = state->getContact();
764 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
765 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
766 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
767 }
768 // Calculate the correction factor.
769 if (trans.x()) {
770 correction = -1.0*overlap / trans.x();
771 if (correction < 0)
772 correction = 1.0;
773 }
774 }
775
776 // Angular dispacement increment.
777 DAVect ang = state->relativeAngularIncrement_;
778 DVect lin_F_old = lin_F_;
779
780 if (lin_mode_ == 0)
781 lin_F_.rx() = overlap * kn_; // absolute mode for normal force calculation
782 else
783 lin_F_.rx() -= correction * trans.x() * kn_; // incremental mode for normal force calculation
784
785 // Normal force can only be positive.
786 lin_F_.rx() = std::max(0.0,lin_F_.x());
787
788 // Calculate the shear force.
789 DVect sforce(0.0);
790 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
791 // Loop over the shear components (note: the 0 component is the normal component)
792 // and calculate the shear force.
793 for (int i=1; i<dim; ++i)
794 sforce.rdof(i) = lin_F_.dof(i) - trans.dof(i) * ks_ * correction;
795
796 // The canFail flag corresponds to whether or not the contact can undergo non-linear
797 // force-displacement response. If the SOLVE ELASTIC command is given then the
798 // canFail state is set to FALSE. Otherwise it is always TRUE.
799 if (state->canFail_) {
800 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
801 double crit = lin_F_.x() * fric_;
802 // The is the magnitude of the shear force.
803 double sfmag = sforce.mag();
804 // Sliding occurs when the magnitude of the shear force is greater than the
805 // critical value.
806 if (sfmag > crit) {
807 // Lower the shear force to the critical value for sliding.
808 double rat = crit / sfmag;
809 sforce *= rat;
810 // Handle the slip_change event if one has been hooked up. Sliding has commenced.
811 if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
812 auto c = state->getContact();
813 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
814 fish::Parameter() };
815 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
816 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
817 }
818 lin_S_ = true;
819 } else {
820 // Handle the slip_change event if one has been hooked up and
821 // the contact was previously sliding. Sliding has ceased.
822 if (lin_S_) {
823 if (cmEvents_[fSlipChange] >= 0) {
824 auto c = state->getContact();
825 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
826 fish::Parameter((qint64)1) };
827 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
828 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
829 }
830 lin_S_ = false;
831 }
832 }
833 }
834
835 // Set the shear components of the total force.
836 for (int i=1; i<dim; ++i)
837 lin_F_.rdof(i) = sforce.dof(i);
838
839 // Rolling resistance
840 DAVect res_M_old = res_M_;
841 if ((fr_ == 0.0) || (kr_==0.0)) {
842 res_M_.fill(0.0);
843 } else {
844 DAVect angStiff(0.0);
845 DAVect MomentInc(0.0);
846#if DIM==3
847 angStiff.rx() = 0.0;
848 angStiff.ry() = kr_;
849#endif
850 angStiff.rz() = kr_;
851 MomentInc = ang * angStiff * (-1.0);
852 res_M_ += MomentInc;
853 if (state->canFail_) {
854 // Account for bending strength
855 double rmax = std::abs(fr_*lin_F_.x());
856 double rmag = res_M_.mag();
857 if (rmag > rmax) {
858 double fac = rmax/rmag;
859 res_M_ *= fac;
860 res_S_ = true;
861 } else {
862 res_S_ = false;
863 }
864 }
865 }
866
867 // Account for dashpot forces if the dashpot structure has been defined.
868 if (dpProps_) {
869 dpProps_->dp_F_.fill(0.0);
870 double vcn(0.0), vcs(0.0);
871 // Calculate the damping coefficients.
872 setDampCoefficients(state->inertialMass_,&vcn,&vcs);
873 // First damp the shear components
874 for (int i=1; i<dim; ++i)
875 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
876 // Damp the normal component
877 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
878 // Need to change behavior based on the dp_mode.
879 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
880 // Limit in tension if not bonded.
881 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
882 dpProps_->dp_F_.rx() = - lin_F_.rx();
883 }
884 if (lin_S_ && dpProps_->dp_mode_ > 1) {
885 // Limit in shear if not sliding.
886 double dfn = dpProps_->dp_F_.rx();
887 dpProps_->dp_F_.fill(0.0);
888 dpProps_->dp_F_.rx() = dfn;
889 }
890 }
891
892 // Adhesive force
893 double a_F_old = a_F_;
894 double gs = state->gap_ - rgap_;
895 a_F_ = 0.0;
896 if ( gs <= 0.0 )
897 a_F_ = a_f0_;
898 else if ( gs < a_d0_ ) // if a_d0_ == 0.0, will not enter, so next line divide by a_d0_ is ok
899 a_F_ = a_f0_ * ( 1.0 - (gs/a_d0_) );
900
901 //Compute energies if energy tracking has been enabled.
902 if (state->trackEnergy_) {
903 assert(energies_);
904 energies_->estrain_ = 0.0;
905 if (kn_)
906 // Calculate the strain energy.
907 energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
908 if (ks_) {
909 DVect s = lin_F_;
910 s.rx() = 0.0;
911 double smag2 = s.mag2();
912 // Add the shear component of the strain energy.
913 energies_->estrain_ += 0.5*smag2 / ks_;
914
915 if (lin_S_) {
916 // If sliding calculate the slip energy and accumulate it.
917 lin_F_old.rx() = 0.0;
918 DVect avg_F_s = (s + lin_F_old)*0.5;
919 DVect u_s_el = (s - lin_F_old) / ks_;
920 DVect u_s(0.0);
921 for (int i=1; i<dim; ++i)
922 u_s.rdof(i) = trans.dof(i);
923 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
924 }
925 }
926 // Add the rolling resistance energy contributions.
927 energies_->errstrain_ = 0.0;
928 if (kr_) {
929 energies_->errstrain_ = 0.5*res_M_.mag2() / kr_;
930 if (res_S_) {
931 // If sliding calculate the slip energy and accumulate it.
932 DAVect avg_M = (res_M_ + res_M_old)*0.5;
933 DAVect t_s_el = (res_M_ - res_M_old) / kr_;
934 energies_->errslip_ -= std::min(0.0,(avg_M | (ang + t_s_el)));
935 }
936 }
937 // Add the adhesive energy contribution:
938 energies_->eadhesive_ -= 0.5*(a_F_old + a_F_) * trans.x();
939
940 if (dpProps_) {
941 // Calculate damping energy (accumulated) if the dashpots are active.
942 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
943 }
944 }
945
946 // This is just a sanity check to ensure, in debug mode, that the force isn't wonky.
947 assert(lin_F_ == lin_F_);
948 return true;
949 }
950
951 bool ContactModelARRLinear::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
952 // Account for thermal expansion in incremental mode
953 if (lin_mode_ == 0 || ts->gapInc_ == 0.0) return false;
954 DVect finc(0.0);
955 finc.rx() = kn_ * ts->gapInc_;
956 lin_F_ -= finc;
957 return true;
958 }
959
960 void ContactModelARRLinear::setForce(const DVect &v,IContact *c) {
961 lin_F(v);
962 if (v.x() > 0)
963 rgap_ = c->getGap() + v.x() / kn_;
964 }
965
966 void ContactModelARRLinear::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
967 // Only called for contacts with wall facets when the wall resolution scheme
968 // is set to full!
969 // Only do something if the contact model is of the same type
970 if (old->getContactModel()->getName().compare("arrlinear",Qt::CaseInsensitive) == 0 && !isBonded()) {
971 ContactModelARRLinear *oldCm = (ContactModelARRLinear *)old;
972#ifdef THREED
973 // Need to rotate just the shear component from oldSystem to newSystem
974
975 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
976 DVect axis = oldSystem.e1() & newSystem.e1();
977 double c, ang, s;
978 DVect re2;
979 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
980 axis = axis.unit();
981 c = oldSystem.e1()|newSystem.e1();
982 if (c > 0)
983 c = std::min(c,1.0);
984 else
985 c = std::max(c,-1.0);
986 ang = acos(c);
987 s = sin(ang);
988 double t = 1. - c;
989 DMatrix<3,3> rm;
990 rm.get(0,0) = t*axis.x()*axis.x() + c;
991 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
992 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
993 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
994 rm.get(1,1) = t*axis.y()*axis.y() + c;
995 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
996 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
997 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
998 rm.get(2,2) = t*axis.z()*axis.z() + c;
999 re2 = rm*oldSystem.e2();
1000 }
1001 else
1002 re2 = oldSystem.e2();
1003 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1004 axis = re2 & newSystem.e2();
1005 DVect2 tpf;
1006 DVect2 tpm;
1007 DMatrix<2,2> m;
1008 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1009 axis = axis.unit();
1010 c = re2|newSystem.e2();
1011 if (c > 0)
1012 c = std::min(c,1.0);
1013 else
1014 c = std::max(c,-1.0);
1015 ang = acos(c);
1016 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1017 ang *= -1;
1018 s = sin(ang);
1019 m.get(0,0) = c;
1020 m.get(1,0) = s;
1021 m.get(0,1) = -m.get(1,0);
1022 m.get(1,1) = m.get(0,0);
1023 tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1024 tpm = m*DVect2(oldCm->res_M_.y(),oldCm->res_M_.z());
1025 } else {
1026 m.get(0,0) = 1.;
1027 m.get(0,1) = 0.;
1028 m.get(1,0) = 0.;
1029 m.get(1,1) = 1.;
1030 tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1031 tpm = DVect2(oldCm->res_M_.y(),oldCm->res_M_.z());
1032 }
1033 DVect pforce = DVect(0,tpf.x(),tpf.y());
1034 //DVect pm = DVect(0,tpm.x(),tpm.y());
1035#else
1036 oldSystem;
1037 newSystem;
1038 DVect pforce = DVect(0,oldCm->lin_F_.y());
1039 //DVect pm = DVect(0,oldCm->res_M_.y());
1040#endif
1041 for (int i=1; i<dim; ++i)
1042 lin_F_.rdof(i) += pforce.dof(i);
1043 if (lin_mode_ && oldCm->lin_mode_)
1044 lin_F_.rx() = oldCm->lin_F_.x();
1045 oldCm->lin_F_ = DVect(0.0);
1046 oldCm->res_M_ = DAVect(0.0);
1047 if (dpProps_ && oldCm->dpProps_) {
1048#ifdef THREED
1049 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1050 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1051#else
1052 pforce = oldCm->dpProps_->dp_F_;
1053#endif
1054 dpProps_->dp_F_ += pforce;
1055 oldCm->dpProps_->dp_F_ = DVect(0.0);
1056 }
1057 if(oldCm->getEnergyActivated()) {
1058 activateEnergy();
1059 energies_->estrain_ = oldCm->energies_->estrain_;
1060 energies_->edashpot_ = oldCm->energies_->edashpot_;
1061 energies_->eslip_ = oldCm->energies_->eslip_;
1062 oldCm->energies_->estrain_ = 0.0;
1063 oldCm->energies_->edashpot_ = 0.0;
1064 oldCm->energies_->eslip_ = 0.0;
1065 }
1066 }
1067 assert(lin_F_ == lin_F_);
1068 }
1069
1070 void ContactModelARRLinear::setNonForcePropsFrom(IContactModel *old) {
1071 // Only called for contacts with wall facets when the wall resolution scheme
1072 // is set to full!
1073 // Only do something if the contact model is of the same type
1074 if (old->getName().compare("arrlinear",Qt::CaseInsensitive) == 0 && !isBonded()) {
1075 ContactModelARRLinear *oldCm = (ContactModelARRLinear *)old;
1076 kn_ = oldCm->kn_;
1077 ks_ = oldCm->ks_;
1078 fric_ = oldCm->fric_;
1079 lin_mode_ = oldCm->lin_mode_;
1080 rgap_ = oldCm->rgap_;
1081 res_fric_ = oldCm->res_fric_;
1082 res_S_ = oldCm->res_S_;
1083 kr_ = oldCm->kr_;
1084 fr_ = oldCm->fr_;
1085 a_f0_ = oldCm->a_f0_;
1086 a_d0_ = oldCm->a_d0_;
1087 userArea_ = oldCm->userArea_;
1088
1089 if (oldCm->dpProps_) {
1090 if (!dpProps_)
1091 dpProps_ = NEWC(dpProps());
1092 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1093 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1094 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1095 }
1096 }
1097 }
1098
1099 DVect ContactModelARRLinear::getForce() const {
1100 DVect ret(lin_F_);
1101 if (dpProps_)
1102 ret += dpProps_->dp_F_;
1103 ret.rdof(0) -= a_F_;
1104 return ret;
1105 }
1106
1107 DAVect ContactModelARRLinear::getMomentOn1(const IContactMechanical *c) const {
1108 DVect force = getForce();
1109 DAVect ret(res_M_);
1110 c->updateResultingTorqueOn1Local(force,&ret);
1111 return ret;
1112 }
1113
1114 DAVect ContactModelARRLinear::getMomentOn2(const IContactMechanical *c) const {
1115 DVect force = getForce();
1116 DAVect ret(res_M_);
1117 c->updateResultingTorqueOn2Local(force,&ret);
1118 return ret;
1119 }
1120
1121 DAVect ContactModelARRLinear::getModelMomentOn1() const {
1122 DAVect ret(res_M_);
1123 return ret;
1124 }
1125
1126 DAVect ContactModelARRLinear::getModelMomentOn2() const {
1127 DAVect ret(res_M_);
1128 return ret;
1129 }
1130
1131 void ContactModelARRLinear::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
1132 ret->clear();
1133 ret->push_back({"kn",ScalarInfo});
1134 ret->push_back({"ks",ScalarInfo});
1135 ret->push_back({"fric",ScalarInfo});
1136 ret->push_back({"lin_force",VectorInfo});
1137 ret->push_back({"lin_slip",ScalarInfo});
1138 ret->push_back({"lin_mode",ScalarInfo});
1139 ret->push_back({"rgap",ScalarInfo});
1140 ret->push_back({"emod",ScalarInfo});
1141 ret->push_back({"kratio",ScalarInfo});
1142 ret->push_back({"dp_nratio",ScalarInfo});
1143 ret->push_back({"dp_sratio",ScalarInfo});
1144 ret->push_back({"dp_mode",ScalarInfo});
1145 ret->push_back({"dp_force",VectorInfo});
1146 ret->push_back({"rr_fric",ScalarInfo});
1147 ret->push_back({"rr_moment",VectorInfo});
1148 ret->push_back({"rr_slip",ScalarInfo});
1149 ret->push_back({"rr_kr",ScalarInfo});
1150 ret->push_back({"adh_f0",ScalarInfo});
1151 ret->push_back({"adh_d0",ScalarInfo});
1152 ret->push_back({"adh_force",VectorInfo});
1153 ret->push_back({"user_area",ScalarInfo});
1154 }
1155
1156 void ContactModelARRLinear::objectPropValues(std::vector<double> *ret,const IContact *c) const {
1157 ret->clear();
1158 ret->push_back(kn());
1159 ret->push_back(ks());
1160 ret->push_back(fric());
1161 ret->push_back(lin_F().mag());
1162 ret->push_back(lin_S());
1163 ret->push_back(lin_mode());
1164 ret->push_back(rgap());
1165 ret->push_back(emod(c));
1166 ret->push_back(kratio());
1167 ret->push_back(dp_nratio());
1168 ret->push_back(dp_sratio());
1169 ret->push_back(dp_mode());
1170 ret->push_back(dp_F().mag());
1171 ret->push_back(res_fric());
1172 ret->push_back(res_M().mag());
1173 ret->push_back(res_S());
1174 ret->push_back(kr());
1175 ret->push_back(a_f0());
1176 ret->push_back(a_d0());
1177 ret->push_back(a_F());
1178 ret->push_back(getArea());
1179 }
1180
1181 double ContactModelARRLinear::emod(const IContact *con) const {
1182 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
1183 if (c ==nullptr) return 0.0;
1184 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
1185 double rsum(0.0);
1186 if (c->getEnd1Curvature().y())
1187 rsum += 1.0/c->getEnd1Curvature().y();
1188 if (c->getEnd2Curvature().y())
1189 rsum += 1.0/c->getEnd2Curvature().y();
1190 if (userArea_) {
1191#ifdef THREED
1192 rsq = std::sqrt(userArea_ / dPi);
1193#else
1194 rsq = userArea_ / 2.0;
1195#endif
1196 rsum = rsq + rsq;
1197 rsq = 1. / rsq;
1198 }
1199#ifdef TWOD
1200 return (kn_ * rsum * rsq / 2.0);
1201#else
1202 return (kn_ * rsum * rsq * rsq) / dPi;
1203#endif
1204 }
1205
1206 double ContactModelARRLinear::kratio() const {
1207 return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
1208 }
1209
1210 void ContactModelARRLinear::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
1211 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
1212 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
1213 }
1214
1215} // namespace cmodelsxd
1216// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Dec 14, 2024 |
