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