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