Hysteretic Model Implementation
See this page for the documentation of this contact model.
contactmodelhysteretic.h
1#pragma once
2// contactmodelhysteretic.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef HYSTERETIC_LIB
7# define HYSTERETIC_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define HYSTERETIC_EXPORT
10#else
11# define HYSTERETIC_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelHysteretic : public ContactModelMechanical {
18 public:
19 enum PropertyKeys { kwHzShear=1
20 , kwHzPoiss
21 , kwFric
22 , kwHzF
23 , kwHzS
24 , kwHzSd
25 , kwHzAlpha
26 , kwDpMode
27 , kwDpEn
28 , kwDpEnMin
29 , kwDpF
30 };
31
32 HYSTERETIC_EXPORT ContactModelHysteretic();
33 HYSTERETIC_EXPORT ContactModelHysteretic(const ContactModelHysteretic &) noexcept;
34 HYSTERETIC_EXPORT const ContactModelHysteretic & operator=(const ContactModelHysteretic &);
35 HYSTERETIC_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
36
37 HYSTERETIC_EXPORT virtual ~ContactModelHysteretic();
38 virtual void copy(const ContactModel *c) override;
39 virtual void archive(ArchiveStream &);
40
41 virtual string getName() const { return "hysteretic"; }
42 virtual void setIndex(int i) { index_=i;}
43 virtual int getIndex() const {return index_;}
44
45 virtual string getProperties() const {
46 return "hz_shear"
47 ",hz_poiss"
48 ",fric"
49 ",hz_force"
50 ",hz_slip"
51 ",hz_mode"
52 ",hz_alpha"
53 ",dp_mode"
54 ",dp_en"
55 ",dp_enmin"
56 ",dp_force"
57 ;
58 }
59
60 enum EnergyKeys { kwEStrain=1,kwESlip,kwEDashpot};
61 virtual string getEnergies() const { return "energy-strain,energy-slip,energy-dashpot";}
62 virtual double getEnergy(uint32 i) const; // Base 1
63 virtual bool getEnergyAccumulate(uint32 i) const; // Base 1
64 virtual void setEnergy(uint32 i,const double &d); // Base 1
65 virtual void activateEnergy() { if (energies_) return; energies_ = NEW Energies();}
66 virtual bool getEnergyActivated() const {return (energies_ !=0);}
67
68 enum FishCallEvents { fActivated=0, fSlipChange};
69 virtual string getFishCallEvents() const { return "contact_activated,slip_change"; }
70 virtual base::Property getProperty(uint32 i,const IContact *) const;
71 virtual bool getPropertyGlobal(uint32 i) const;
72 virtual bool setProperty(uint32 i,const base::Property &v,IContact *);
73 virtual bool getPropertyReadOnly(uint32 i) const;
74
75 virtual bool supportsInheritance(uint32 i) const;
76 virtual bool getInheritance(uint32 i) const { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
77 virtual void setInheritance(uint32 i,bool b) { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
78
79 virtual uint32 getMinorVersion() const;
80
81 virtual bool validate(ContactModelMechanicalState *state,const double ×tep);
82 virtual bool endPropertyUpdated(const string &name,const IContactMechanical *c);
83 virtual bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep);
84 virtual DVect2 getEffectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_;}
85 virtual DAVect getEffectiveRotationalStiffness() const { return DAVect(0.0);}
86
87 virtual ContactModelHysteretic *clone() const override { return NEW ContactModelHysteretic(); }
88 virtual double getActivityDistance() const {return 0.0;}
89 virtual bool isOKToDelete() const { return !isBonded(); }
90 virtual void resetForcesAndMoments() { hz_F(DVect(0.0)); dp_F(DVect(0.0)); if (energies_) energies_->estrain_ = 0.0;}
91 virtual void setForce(const DVect &v,IContact *) { hz_F(v); }
92 virtual void setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }
93 virtual double getArea() const { return 0.0; }
94
95 virtual bool checkActivity(const double &gap) { return gap <= 0.0; }
96
97 virtual bool isSliding() const { return hz_slip_; }
98 virtual bool isBonded() const { return false; }
99 virtual void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes());
100 virtual void setNonForcePropsFrom(IContactModel *oldCM);
101
102 const double & hz_shear() const {return hz_shear_;}
103 void hz_shear(const double &d) {hz_shear_=d;}
104 const double & hz_poiss() const {return hz_poiss_;}
105 void hz_poiss(const double &d) {hz_poiss_=d;}
106 const double & fric() const {return fric_;}
107 void fric(const double &d) {fric_=d;}
108 uint32 hz_mode() const {return hz_mode_;}
109 void hz_mode(uint32 i) {hz_mode_=i;}
110 const DVect & hz_F() const {return hz_F_;}
111 void hz_F(const DVect &f) { hz_F_=f;}
112 bool hz_S() const {return hz_slip_;}
113 void hz_S(bool b) { hz_slip_=b;}
114 const double & hz_alpha() const {return hz_alpha_;}
115 void hz_alpha(const double &d) {hz_alpha_=d;}
116 int dp_mode() const {return dp_mode_;}
117 void dp_mode(int i) {dp_mode_=i;}
118 const double & dp_en() const {return dp_en_;}
119 void dp_en(const double &d) {dp_en_=d;}
120 const double & dp_enmin() const {return dp_enmin_;}
121 void dp_enmin(const double &d) {dp_enmin_=d;}
122 const DVect & dp_F() const {return dp_F_;}
123 void dp_F(const DVect &f) { dp_F_=f;}
124 const double & hn() const {return hn_;}
125 void hn(const double &d) {hn_=d;}
126 const double & hs() const {return hs_;}
127 void hs(const double &d) {hs_=d;}
128 const double & vni() const {return vni_;}
129 void vni(const double &d) {vni_=d;}
130 double pfac() const {return pfac_;}
131 void pfac(int i) {pfac_=i;}
132
133 bool hasEnergies() const {return energies_ ? true:false;}
134 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
135 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
136 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
137 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
138 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
139 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
140
141 uint32 inheritanceField() const {return inheritanceField_;}
142 void inheritanceField(uint32 i) {inheritanceField_ = i;}
143
144 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
145 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
146
147 /// Return the total force that the contact model holds.
148 virtual DVect getForce() const;
149
150 /// Return the total moment on 1 that the contact model holds
151 virtual DAVect getMomentOn1(const IContactMechanical *) const;
152
153 /// Return the total moment on 1 that the contact model holds
154 virtual DAVect getMomentOn2(const IContactMechanical *) const;
155
156 virtual DAVect getModelMomentOn1() const;
157 virtual DAVect getModelMomentOn2() const;
158
159 // Used to efficiently get properties from the contact model for the object pane.
160 // List of properties for the object pane, comma separated.
161 // All properties will be cast to doubles for comparison. No other comparisons
162 // are supported. This may not be the same as the entire property list.
163 // Return property name and type for plotting.
164 void objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *) const override;
165 // All properties cast to doubles - this is what can be compared.
166 void objectPropValues(std::vector<double> *,const IContact *) const override;
167
168 private:
169 static int index_;
170
171 bool updateStiffCoef(const IContactMechanical *con);
172 bool updateEndStiffCoef(const IContactMechanical *con);
173 bool updateEndFric(const IContactMechanical *con);
174 void updateEffectiveStiffness(ContactModelMechanicalState *state);
175 // inheritance fields
176 uint32 inheritanceField_;
177
178 // hertz model
179 double hz_shear_ = 0.0; // Shear modulus
180 double hz_poiss_ = 0.0; // Poisson ratio
181 double fric_ = 0.0; // Coulomb friction coefficient
182 DVect hz_F_ = DVect(0.0); // Force carried in the hertz model
183 bool hz_slip_ = false; // the current sliding state
184 uint32 hz_mode_ = 0; // specifies down-scaling of the shear force when normal unloading occurs
185 double hz_alpha_ = 1.5; // alpha exponent
186 int dp_mode_ = 0; // Damping scheme mode
187 double dp_en_ = 1.0; // normal restitution coefficient
188 double dp_enmin_ = 0.0; // minimal normal restitution coefficient in default mode
189 DVect dp_F_ = DVect(0.0); // Damping Force
190
191 // energies
192 struct Energies {
193 Energies() : estrain_(0.0), eslip_(0.0),edashpot_(0.0) {}
194 double estrain_; // elastic energy stored in contact
195 double eslip_; // work dissipated by friction
196 double edashpot_; // work dissipated by dashpots
197 };
198 Energies * energies_ = nullptr;
199
200 double hn_ = 0.0; // normal stiffness coefficient
201 double hs_ = 0.0; // shear stiffness coefficient
202 DVect2 effectiveTranslationalStiffness_ = DVect2(0.0); // effective stiffness
203 double vni_ = 0.0; // normal impact velocity
204 double pfac_ = 0.0; // previous damping factor
205
206 };
207
208} // namespace cmodelsxd
209// EoF
contactmodelhysteretic.cpp
1// contactmodelhysteretic.cpp
2#include "contactmodelhysteretic.h"
3
4#include "../version.txt"
5#include "fish/src/parameter.h"
6#include "utility/src/tptr.h"
7#include "shared/src/mathutil.h"
8
9#include "module/interface/icontact.h"
10#include "module/interface/icontactmechanical.h"
11#include "module/interface/ifishcalllist.h"
12#include "module/interface/ipiece.h"
13#include "kernel/interface/iprogram.h"
14
15
16#ifdef HYSTERETIC_LIB
17#ifdef _WIN32
18 int __stdcall DllMain(void *,unsigned, void *) {
19 return 1;
20 }
21#endif
22
23 extern "C" EXPORT_TAG const char *getName() {
24#if DIM==3
25 return "contactmodelmechanical3dhysteretic";
26#else
27 return "contactmodelmechanical2dhysteretic";
28#endif
29 }
30
31 extern "C" EXPORT_TAG unsigned getMajorVersion() {
32 return MAJOR_VERSION;
33 }
34
35 extern "C" EXPORT_TAG unsigned getMinorVersion() {
36 return MINOR_VERSION;
37 }
38
39 extern "C" EXPORT_TAG void *createInstance() {
40 cmodelsxd::ContactModelHysteretic *m = NEW cmodelsxd::ContactModelHysteretic();
41 return (void *)m;
42 }
43#endif // HYSTERETIC_LIB
44
45namespace cmodelsxd {
46
47 static const uint32 shearMask = 0x00002;
48 static const uint32 poissMask = 0x00004;
49 static const uint32 fricMask = 0x00008;
50
51 using namespace itasca;
52
53 int ContactModelHysteretic::index_ = -1;
54 uint32 ContactModelHysteretic::getMinorVersion() const { return MINOR_VERSION; }
55
56 ContactModelHysteretic::ContactModelHysteretic() : inheritanceField_(shearMask|poissMask|fricMask) {
57 }
58
59 ContactModelHysteretic::ContactModelHysteretic(const ContactModelHysteretic &m) noexcept {
60 inheritanceField(shearMask|poissMask|fricMask);
61 this->copy(&m);
62 }
63
64 const ContactModelHysteretic & ContactModelHysteretic::operator=(const ContactModelHysteretic &m) {
65 inheritanceField(shearMask|poissMask|fricMask);
66 this->copy(&m);
67 return *this;
68 }
69
70 void ContactModelHysteretic::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
71 s->addToStorage<ContactModelHysteretic>(*this,ww);
72 }
73
74 ContactModelHysteretic::~ContactModelHysteretic() {
75 if (energies_)
76 delete energies_;
77 }
78
79 void ContactModelHysteretic::archive(ArchiveStream &stream) {
80 stream & hz_shear_;
81 stream & hz_poiss_;
82 stream & fric_;
83 stream & hz_F_;
84 stream & hz_slip_;
85 stream & hz_mode_;
86 stream & hz_alpha_;
87 stream & dp_mode_;
88 stream & dp_en_;
89 stream & dp_enmin_;
90 stream & dp_F_;
91 stream & hn_;
92 stream & hs_;
93 stream & vni_;
94 stream & pfac_;
95
96 if (stream.getArchiveState()==ArchiveStream::Save) {
97 bool b = false;
98 if (energies_) {
99 b = true;
100 stream & b;
101 stream & energies_->estrain_;
102 stream & energies_->eslip_;
103 stream & energies_->edashpot_;
104 } else
105 stream & b;
106 } else {
107 bool b(false);
108 stream & b;
109 if (b) {
110 if (!energies_)
111 energies_ = NEW Energies();
112 stream & energies_->estrain_;
113 stream & energies_->eslip_;
114 stream & energies_->edashpot_;
115 }
116 }
117
118 stream & inheritanceField_;
119 stream & effectiveTranslationalStiffness_;
120 }
121
122 void ContactModelHysteretic::copy(const ContactModel *cm) {
123 ContactModelMechanical::copy(cm);
124 const ContactModelHysteretic *in = dynamic_cast<const ContactModelHysteretic*>(cm);
125 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
126
127 hz_shear(in->hz_shear());
128 hz_poiss(in->hz_poiss());
129 fric(in->fric());
130 hz_F(in->hz_F());
131 hz_S(in->hz_S());
132 hz_mode(in->hz_mode());
133 hz_alpha(in->hz_alpha());
134 dp_mode(in->dp_mode());
135 dp_en(in->dp_en());
136 dp_enmin(in->dp_enmin());
137 hn(in->hn());
138 hs(in->hs());
139 vni(in->vni());
140 pfac(in->pfac());
141 if (in->hasEnergies()) {
142 if (!energies_)
143 energies_ = NEW Energies();
144 estrain(in->estrain());
145 eslip(in->eslip());
146 edashpot(in->edashpot());
147 }
148 inheritanceField(in->inheritanceField());
149 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
150 }
151
152 base::Property ContactModelHysteretic::getProperty(uint32 i,const IContact *) const {
153 // The IContact pointer may be a nullptr!
154 base::Property var;
155 switch (i) {
156 case kwHzShear: return hz_shear_;
157 case kwHzPoiss: return hz_poiss_;
158 case kwFric: return fric_;
159 case kwHzF: var.setValue(hz_F_); return var;
160 case kwHzS: return hz_slip_;
161 case kwHzSd: return hz_mode_;
162 case kwHzAlpha: return hz_alpha_;
163 case kwDpMode: return dp_mode_;
164 case kwDpEn: return dp_en_;
165 case kwDpEnMin: return dp_enmin_;
166 case kwDpF: var.setValue(dp_F_); return var;
167 }
168 assert(0);
169 return base::Property();
170 }
171
172 bool ContactModelHysteretic::getPropertyGlobal(uint32 i) const {
173 switch (i) {
174 case kwHzF: return false;
175 case kwDpF: return false;
176 }
177 return true;
178 }
179
180 bool ContactModelHysteretic::setProperty(uint32 i,const base::Property &v,IContact *) {
181 switch (i) {
182 case kwHzShear: {
183 if (!v.canConvert<double>())
184 throw Exception("hz_shear must be a double.");
185 double val(v.toDouble());
186 if (val<0.0)
187 throw Exception("Negative shear modulus (hz_shear) not allowed.");
188 hz_shear_ = val;
189 return true;
190 }
191 case kwHzPoiss: {
192 if (!v.canConvert<double>())
193 throw Exception("hz_poiss must be a double.");
194 double val(v.toDouble());
195 if (val<=-1.0 || val>0.5)
196 throw Exception("Poisson ratio (hz_poiss) must be in range (-1.0,0.5] ");
197 hz_poiss_ = val;
198 return true;
199 }
200 case kwFric: {
201 if (!v.canConvert<double>())
202 throw Exception("fric must be a double.");
203 double val(v.toDouble());
204 if (val<0.0)
205 throw Exception("Negative fric not allowed.");
206 fric_ = val;
207 return false;
208 }
209 case kwHzSd: {
210 if (!v.canConvert<int64>())
211 throw Exception("hz_mode must be 0 or 1.");
212 uint32 val(v.toUInt());
213 if (val >1)
214 throw Exception("hz_mode must be 0 or 1.");
215 hz_mode_ = val;
216 return false;
217 }
218 case kwHzAlpha: {
219 if (!v.canConvert<double>())
220 throw Exception("hz_alpha must be a double.");
221 double val(v.toDouble());
222 if (val<0.0)
223 throw Exception("Alpha exponent (hz_alpha) must be positive.");
224 hz_alpha_ = val;
225 return false;
226 }
227 case kwDpMode: {
228 if (!v.canConvert<int64>())
229 throw Exception("dp_mode must be an integer.");
230 int val(v.toInt());
231 dp_mode_ = val;
232 return false;
233 }
234 case kwDpEn: {
235 if (!v.canConvert<double>())
236 throw Exception("dp_en must be a double.");
237 double val(v.toDouble());
238 if (val<0.0 || val>1.0)
239 throw Exception("Restitution coefficient (dp_en) must be in range [0.0,1.0] ");
240 dp_en_ = val;
241 return false;
242 }
243 case kwDpEnMin: {
244 if (!v.canConvert<double>())
245 throw Exception("dp_enmin must be a double.");
246 double val(v.toDouble());
247 if (val<0.0 || val>1.0)
248 throw Exception("Minimal restitution coefficient (dp_enmin) must be in range [0.0,1.0] ");
249 dp_enmin_ = val;
250 return false;
251 }
252 }
253 return false;
254 }
255
256 bool ContactModelHysteretic::getPropertyReadOnly(uint32 i) const {
257 switch (i) {
258 case kwHzF:
259 case kwHzS:
260 case kwDpF:
261 return true;
262 default:
263 break;
264 }
265 return false;
266 }
267
268 bool ContactModelHysteretic::supportsInheritance(uint32 i) const {
269 switch (i) {
270 case kwHzShear:
271 case kwHzPoiss:
272 case kwFric:
273 return true;
274 default:
275 break;
276 }
277 return false;
278 }
279
280 double ContactModelHysteretic::getEnergy(uint32 i) const {
281 double ret(0.0);
282 if (!energies_)
283 return ret;
284 switch (i) {
285 case kwEStrain: return energies_->estrain_;
286 case kwESlip: return energies_->eslip_;
287 case kwEDashpot: return energies_->edashpot_;
288 }
289 assert(0);
290 return ret;
291 }
292
293 bool ContactModelHysteretic::getEnergyAccumulate(uint32 i) const {
294 switch (i) {
295 case kwEStrain: return false;
296 case kwESlip: return true;
297 case kwEDashpot: return true;
298 }
299 assert(0);
300 return false;
301 }
302
303 void ContactModelHysteretic::setEnergy(uint32 i,const double &d) {
304 if (!energies_) return;
305 switch (i) {
306 case kwEStrain: energies_->estrain_ = d; return;
307 case kwESlip: energies_->eslip_ = d; return;
308 case kwEDashpot: energies_->edashpot_ = d; return;
309 }
310 assert(0);
311 return;
312 }
313
314 bool ContactModelHysteretic::validate(ContactModelMechanicalState *state,const double &) {
315 assert(state);
316 const IContactMechanical *c = state->getMechanicalContact();
317 assert(c);
318
319 if (state->trackEnergy_)
320 activateEnergy();
321
322 updateStiffCoef(c);
323 if ((inheritanceField_ & shearMask) || (inheritanceField_ & poissMask))
324 updateEndStiffCoef(c);
325
326 if (inheritanceField_ & fricMask)
327 updateEndFric(c);
328
329 updateEffectiveStiffness(state);
330 return checkActivity(state->gap_);
331 }
332
333 bool ContactModelHysteretic::updateStiffCoef(const IContactMechanical *con) {
334 double hnold = hn_;
335 double hsold = hs_;
336 double c12 = con->getEnd1Curvature().y();
337 double c22 = con->getEnd2Curvature().y();
338 double reff = c12+c22;
339 if (reff == 0.0)
340 throw Exception("Hysteretic contact model undefined for 2 non-curved surfaces");
341 reff = 2.0 /reff;
342 hn_ = 2.0/3.0 * (hz_shear_/(1 -hz_poiss_)) * sqrt(2.0*reff);
343 hs_ = (2.0*pow(hz_shear_*hz_shear_*3.0*(1-hz_poiss_)*(reff),1.0/3.0)) / (2.0- hz_poiss_);
344 return ( (hn_ != hnold) || (hs_ != hsold) );
345 }
346
347 static const string gstr("hz_shear");
348 static const string nustr("hz_poiss");
349 bool ContactModelHysteretic::updateEndStiffCoef(const IContactMechanical *con) {
350 assert(con);
351 double g1 = hz_shear_;
352 double g2 = hz_shear_;
353 double nu1 = hz_poiss_;
354 double nu2 = hz_poiss_;
355 base::Property vg1 = con->getEnd1()->getProperty(gstr);
356 base::Property vg2 = con->getEnd2()->getProperty(gstr);
357 base::Property vnu1 = con->getEnd1()->getProperty(nustr);
358 base::Property vnu2 = con->getEnd2()->getProperty(nustr);
359 if (vg1.isValid() && vg2.isValid()) {
360 g1 = vg1.toDouble();
361 g2 = vg2.toDouble();
362 if (g1 < 0.0 || g2 < 0.0)
363 throw Exception("Negative shear modulus not allowed in Hysteretic contact model");
364 }
365 if (vnu1.isValid() && vnu2.isValid()) {
366 nu1 = vnu1.toDouble();
367 nu2 = vnu2.toDouble();
368 if (nu1 <= -1.0 || nu1 > 0.5 || nu2 <= -1.0 || nu2 > 0.5)
369 throw Exception("Poisson ratio should be in range (-1.0,0.5] in Hysteretic contact model");
370 }
371 if (g1*g2 == 0.0) return false;
372 double es = 1.0 / ((1.0-nu1) / (2.0*g1) + (1.0-nu2) / (2.0*g2));
373 double gs = 1.0 / ((2.0-nu1) / g1 + (2.0-nu2) /g2);
374 hz_poiss_ = (4.0*gs-es)/(2.0*gs-es);
375 hz_shear_ = 2.0*gs*(2-hz_poiss_);
376 if (hz_shear_ < 0.0)
377 throw Exception("Negative shear modulus not allowed in Hysteretic contact model");
378 if (hz_poiss_ <= -1.0 || hz_poiss_ > 0.5)
379 throw Exception("Poisson ratio should be in range (-1.0,0.5] in Hysteretic contact model");
380 return updateStiffCoef(con);
381 }
382
383 static const string fricstr("fric");
384 bool ContactModelHysteretic::updateEndFric(const IContactMechanical *con) {
385 assert(con);
386 base::Property v1 = con->getEnd1()->getProperty(fricstr);
387 base::Property v2 = con->getEnd2()->getProperty(fricstr);
388 if (!v1.isValid() || !v2.isValid())
389 return false;
390 double fric1 = std::max(0.0,v1.toDouble());
391 double fric2 = std::max(0.0,v2.toDouble());
392 double val = fric_;
393 fric_ = std::min(fric1,fric2);
394 return ( (fric_ != val) );
395 }
396
397 bool ContactModelHysteretic::endPropertyUpdated(const string &name,const IContactMechanical *c) {
398 assert(c);
399 StringList availableProperties = split(replace(simplified(getProperties())," ",""),",");
400 auto idx = findRegex(availableProperties,name);
401 if (idx==string::npos) return false;
402 idx += 1;
403 bool ret = false;
404
405 switch(idx) {
406 case kwHzShear: {
407 if (inheritanceField_ & shearMask)
408 ret = updateEndStiffCoef(c);
409 break;
410 }
411 case kwHzPoiss: {
412 if (inheritanceField_ & poissMask)
413 ret = updateEndStiffCoef(c);
414 break;
415 }
416 case kwFric: {
417 if (inheritanceField_ & fricMask)
418 ret = updateEndFric(c);
419 break;
420 }
421 }
422 return ret;
423 }
424
425 void ContactModelHysteretic::updateEffectiveStiffness(ContactModelMechanicalState *state) {
426 effectiveTranslationalStiffness_ = DVect2(hn_,hs_);
427 if (state->gap_ >= 0.0) return;
428 double overlap = - state->gap_;
429 double kn = 1.5*hn_*sqrt(overlap);
430 double ks = hs_ * pow(hz_F_.x(),(1.0/3.0));
431 DVect2 ret(kn,ks);
432 effectiveTranslationalStiffness_ = ret;
433 }
434
435 bool ContactModelHysteretic::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
436 assert(state);
437 bool firstActive = false;
438 if (state->activated()) {
439 if (cmEvents_[fActivated] >= 0) {
440 auto c = state->getContact();
441 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
442 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
443 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
444 }
445 firstActive = true;
446 }
447
448 double overlap = - state->gap_;
449 DVect trans = state->relativeTranslationalIncrement_;
450#ifdef THREED
451 DVect norm(trans.x(),0.0,0.0);
452#else
453 DVect norm(trans.x(),0.0);
454#endif
455
456 //DAVect ang = state->relativeAngularIncrement_;
457 DVect hz_F_old = hz_F_;
458 hz_F_ = DVect(0.0);
459 dp_F_ = DVect(0.0);
460 double ks_old = hs_ * pow(hz_F_old.x(),(1.0/3.0));
461 DVect fs_old = hz_F_old;
462 fs_old.rx() = 0.0;
463 if (overlap > 0) {
464
465 double kn = 1.5 * hn_ * sqrt(overlap);
466
467 double vn = norm.x() / timestep;
468 if (firstActive) {
469 if (vn <= 0.0)
470 vni_= vn;
471 else
472 vni_ = 0.0;
473 fs_old = DVect(0.0);
474 ks_old = 0.0;
475 hz_F_old = DVect(0.0);
476 pfac_ = 0.0;
477 }
478 hz_F_.rx() = hn_*pow(overlap,hz_alpha_);
479
480 double ks = hs_ * pow(hz_F_.x(),(1.0/3.0));
481 effectiveTranslationalStiffness_ = DVect2(kn,ks);
482
483 //damping force
484 if (std::abs(vni_) <= limits<double>::epsilon()*1000.0)
485 dp_F_.rx() = 0.0;
486 else {
487 double ratio = vn/vni_;
488 double fac(0.0);
489 switch (dp_mode_) {
490 default:
491 {
492 double used = max(dp_en_,dp_enmin_);
493 fac = max(0.0,(1.0-used*used)/used)*ratio; //Gonthier et al.
494 }
495 break;
496 case 1:
497 fac = 0.75*(1.0-dp_en_*dp_en_)*ratio; // Lankarani et al (1989).
498 break;
499 case 2:
500 fac = 0.75*(1.0-dp_en_*dp_en_)*exp(2.0*(1-dp_en_))*ratio; // Zhiying et al.
501 break;
502 }
503
504 if (fac <= -1.0) { // sucking
505 if (pfac_ >= 0.0) { //switch in one timestep from pushing to sucking - instability
506 fac = 0.0;
507 pfac_ = 0.0;
508 } else
509 pfac_ = fac;
510 } else if (pfac_ != 0.0 && abs(fac) > 10.0*abs(pfac_)) { // growing too fast - instability
511 fac = 0.0;
512 pfac_ = 0.0;
513 } else
514 pfac_ = fac;
515
516 dp_F_.rx() = hz_F_.x()*fac;
517 }
518
519 DVect u_s = trans;
520 u_s.rx() = 0.0;
521 DVect vec = u_s * ks;
522 if (hz_mode_ && (hz_F_.x() < hz_F_old.x())) {
523 double rat = ks / ks_old;
524 fs_old *= rat;
525 }
526 DVect fs = fs_old - vec;
527
528 if (state->canFail_) {
529 // resolve sliding
530 double crit = hz_F_.x() * fric_;
531 double sfmag = fs.mag();
532 if (sfmag > crit) {
533 double rat = crit / sfmag;
534 fs *= rat;
535 if (!hz_slip_ && cmEvents_[fSlipChange] >= 0) {
536 auto c = state->getContact();
537 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
538 fish::Parameter() };
539 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
540 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
541 }
542 hz_slip_ = true;
543 } else {
544 if (hz_slip_) {
545 if (cmEvents_[fSlipChange] >= 0) {
546 auto c = state->getContact();
547 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
548 fish::Parameter((int64)1) };
549 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
550 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
551 }
552 hz_slip_ = false;
553 }
554 }
555 }
556
557 fs.rx() = hz_F_.x();
558 hz_F_ = fs; // total force in hertz part
559
560 // 5) Compute energies
561 if (state->trackEnergy_) {
562 assert(energies_);
563 energies_->estrain_ = 0.0;
564 if (kn)
565 energies_->estrain_ = hz_alpha_*hz_F_.x()*hz_F_.x()/((hz_alpha_+1)*kn);
566 if (ks) {
567 DVect s = hz_F_;
568 s.rx() = 0.0;
569 double smag2 = s.mag2();
570 energies_->estrain_ += 0.5*smag2 / ks;
571 if (hz_slip_) {
572 hz_F_old.rx() = 0.0;
573 DVect avg_F_s = (s + hz_F_old)*0.5;
574 DVect u_s_el = (s - hz_F_old) / ks;
575 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
576 }
577 }
578 energies_->edashpot_ -= dp_F_ | trans;
579 }
580
581 } else {
582 hz_F_ = DVect(0.0);
583 dp_F_ = DVect(0.0);
584 pfac_ = 0.0;
585 }
586
587
588 return true;
589 }
590
591 void ContactModelHysteretic::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
592 // Only do something if the contact model is of the same type
593 if (equal(old->getContactModel()->getName(),"hysteretic")) {
594 ContactModelHysteretic *oldCm = (ContactModelHysteretic *)old;
595#ifdef THREED
596 // Need to rotate just the shear component from oldSystem to newSystem
597
598 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
599 DVect axis = oldSystem.e1() & newSystem.e1();
600 double c, ang, s;
601 DVect re2;
602 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
603 axis = axis.unit();
604 c = oldSystem.e1()|newSystem.e1();
605 if (c > 0)
606 c = std::min(c,1.0);
607 else
608 c = std::max(c,-1.0);
609 ang = acos(c);
610 s = sin(ang);
611 double t = 1. - c;
612 DMatrix<3,3> rm;
613 rm.get(0,0) = t*axis.x()*axis.x() + c;
614 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
615 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
616 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
617 rm.get(1,1) = t*axis.y()*axis.y() + c;
618 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
619 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
620 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
621 rm.get(2,2) = t*axis.z()*axis.z() + c;
622 re2 = rm*oldSystem.e2();
623 } else
624 re2 = oldSystem.e2();
625
626 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
627 axis = re2 & newSystem.e2();
628 DVect2 tpf;
629 DMatrix<2,2> m;
630 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
631 axis = axis.unit();
632 c = re2|newSystem.e2();
633 if (c > 0)
634 c = std::min(c,1.0);
635 else
636 c = std::max(c,-1.0);
637 ang = acos(c);
638 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
639 ang *= -1;
640 s = sin(ang);
641 m.get(0,0) = c;
642 m.get(1,0) = s;
643 m.get(0,1) = -m.get(1,0);
644 m.get(1,1) = m.get(0,0);
645 tpf = m*DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
646 } else {
647 m.get(0,0) = 1.;
648 m.get(0,1) = 0.;
649 m.get(1,0) = 0.;
650 m.get(1,1) = 1.;
651 tpf = DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
652 }
653 DVect pforce = DVect(0,tpf.x(),tpf.y());
654#else
655 oldSystem;
656 newSystem;
657 DVect pforce = DVect(0,oldCm->hz_F_.y());
658#endif
659 for (int i=1; i<dim; ++i)
660 hz_F_.rdof(i) += pforce.dof(i);
661 oldCm->hz_F_ = DVect(0.0);
662
663 if(oldCm->getEnergyActivated()) {
664 activateEnergy();
665 energies_->estrain_ = oldCm->energies_->estrain_;
666 energies_->eslip_ = oldCm->energies_->eslip_;
667 energies_->edashpot_ = oldCm->energies_->edashpot_;
668 oldCm->energies_->estrain_ = 0.0;
669 oldCm->energies_->eslip_ = 0.0;
670 oldCm->energies_->edashpot_ = 0.0;
671 }
672 }
673 }
674
675 void ContactModelHysteretic::setNonForcePropsFrom(IContactModel *old) {
676 // Only do something if the contact model is of the same type
677 if (equal(old->getName(),"hysteretic") && !isBonded()) {
678 ContactModelHysteretic *oldCm = (ContactModelHysteretic *)old;
679 hn_ = oldCm->hn_;
680 hs_ = oldCm->hs_;
681 fric_ = oldCm->fric_;
682 vni_ = oldCm->vni_;
683 pfac_ = oldCm->pfac_;
684 }
685 }
686
687 DVect ContactModelHysteretic::getForce() const {
688 DVect ret(hz_F_);
689 ret += dp_F_;
690 return ret;
691 }
692
693 DAVect ContactModelHysteretic::getMomentOn1(const IContactMechanical *c) const {
694 DVect force = getForce();
695 DAVect ret(0.0);
696 c->updateResultingTorqueOn1Local(force,&ret);
697 return ret;
698 }
699
700 DAVect ContactModelHysteretic::getMomentOn2(const IContactMechanical *c) const {
701 DVect force = getForce();
702 DAVect ret(0.0);
703 c->updateResultingTorqueOn2Local(force,&ret);
704 return ret;
705 }
706
707 DAVect ContactModelHysteretic::getModelMomentOn1() const {
708 DAVect ret(0.0);
709 return ret;
710 }
711
712 DAVect ContactModelHysteretic::getModelMomentOn2() const {
713 DAVect ret(0.0);
714 return ret;
715 }
716
717 void ContactModelHysteretic::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
718 ret->clear();
719 ret->push_back({"hz_shear",ScalarInfo});
720 ret->push_back({"hz_poiss",ScalarInfo});
721 ret->push_back({"fric",ScalarInfo});
722 ret->push_back({"hz_force",VectorInfo});
723 ret->push_back({"hz_slip",ScalarInfo});
724 ret->push_back({"hz_mode",ScalarInfo});
725 ret->push_back({"hz_alpha",ScalarInfo});
726 ret->push_back({"dp_mode",ScalarInfo});
727 ret->push_back({"dp_en",ScalarInfo});
728 ret->push_back({"dp_enmin",ScalarInfo});
729 ret->push_back({"dp_force",VectorInfo});
730 }
731
732 void ContactModelHysteretic::objectPropValues(std::vector<double> *ret,const IContact *) const {
733 FP_S;
734 ret->clear();
735 ret->push_back(hz_shear_);
736 ret->push_back(hz_poiss_);
737 ret->push_back(fric_);
738 ret->push_back(hz_F_.mag());
739 ret->push_back(hz_slip_);
740 ret->push_back(hz_mode_);
741 ret->push_back(hz_alpha_);
742 ret->push_back(dp_mode_);
743 ret->push_back(dp_en_);
744 ret->push_back(dp_enmin_);
745 ret->push_back(dp_F_.mag());
746 FP_S;
747 }
748
749} // namespace cmodelsxd
750// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
