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