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