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 QString getName() const { return "hysteretic"; }
42 virtual void setIndex(int i) { index_=i;}
43 virtual int getIndex() const {return index_;}
44
45 virtual QString 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 QString 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_ = NEWC(Energies());}
66 virtual bool getEnergyActivated() const {return (energies_ !=0);}
67
68 enum FishCallEvents { fActivated=0, fSlipChange};
69 virtual QString getFishCallEvents() const { return "contact_activated,slip_change"; }
70 virtual QVariant getProperty(uint32 i,const IContact *) const;
71 virtual bool getPropertyGlobal(uint32 i) const;
72 virtual bool setProperty(uint32 i,const QVariant &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 QString &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 NEWC(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<QString,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 = NEWC(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_ = NEWC(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_ = NEWC(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 QVariant ContactModelHysteretic::getProperty(uint32 i,const IContact *) const {
153 // The IContact pointer may be a nullptr!
154 QVariant 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 QVariant();
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 QVariant &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<uint32>())
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<int>())
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 QString gstr("hz_shear");
348 static const QString 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 QVariant vg1 = con->getEnd1()->getProperty(gstr);
356 QVariant vg2 = con->getEnd2()->getProperty(gstr);
357 QVariant vnu1 = con->getEnd1()->getProperty(nustr);
358 QVariant 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 QString fricstr("fric");
384 bool ContactModelHysteretic::updateEndFric(const IContactMechanical *con) {
385 assert(con);
386 QVariant v1 = con->getEnd1()->getProperty(fricstr);
387 QVariant 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 QString &name,const IContactMechanical *c) {
398 assert(c);
399 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
400 QRegExp rx(name,Qt::CaseInsensitive);
401 int idx = availableProperties.indexOf(rx)+1;
402 bool ret=false;
403
404 if (idx<=0)
405 return ret;
406
407 switch(idx) {
408 case kwHzShear: {
409 if (inheritanceField_ & shearMask)
410 ret = updateEndStiffCoef(c);
411 break;
412 }
413 case kwHzPoiss: {
414 if (inheritanceField_ & poissMask)
415 ret = updateEndStiffCoef(c);
416 break;
417 }
418 case kwFric: {
419 if (inheritanceField_ & fricMask)
420 ret = updateEndFric(c);
421 break;
422 }
423 }
424 return ret;
425 }
426
427 void ContactModelHysteretic::updateEffectiveStiffness(ContactModelMechanicalState *state) {
428 effectiveTranslationalStiffness_ = DVect2(hn_,hs_);
429 if (state->gap_ >= 0.0) return;
430 double overlap = - state->gap_;
431 double kn = 1.5*hn_*sqrt(overlap);
432 double ks = hs_ * pow(hz_F_.x(),(1.0/3.0));
433 DVect2 ret(kn,ks);
434 effectiveTranslationalStiffness_ = ret;
435 }
436
437 bool ContactModelHysteretic::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
438 assert(state);
439 bool firstActive = false;
440 if (state->activated()) {
441 if (cmEvents_[fActivated] >= 0) {
442 auto c = state->getContact();
443 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
444 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
445 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
446 }
447 firstActive = true;
448 }
449
450 double overlap = - state->gap_;
451 DVect trans = state->relativeTranslationalIncrement_;
452#ifdef THREED
453 DVect norm(trans.x(),0.0,0.0);
454#else
455 DVect norm(trans.x(),0.0);
456#endif
457
458 //DAVect ang = state->relativeAngularIncrement_;
459 DVect hz_F_old = hz_F_;
460 hz_F_ = DVect(0.0);
461 dp_F_ = DVect(0.0);
462 double ks_old = hs_ * pow(hz_F_old.x(),(1.0/3.0));
463 DVect fs_old = hz_F_old;
464 fs_old.rx() = 0.0;
465 if (overlap > 0) {
466
467 double kn = 1.5 * hn_ * sqrt(overlap);
468
469 double vn = norm.x() / timestep;
470 if (firstActive) {
471 if (vn <= 0.0)
472 vni_= vn;
473 else
474 vni_ = 0.0;
475 fs_old = DVect(0.0);
476 ks_old = 0.0;
477 hz_F_old = DVect(0.0);
478 pfac_ = 0.0;
479 }
480 hz_F_.rx() = hn_*pow(overlap,hz_alpha_);
481
482 double ks = hs_ * pow(hz_F_.x(),(1.0/3.0));
483 effectiveTranslationalStiffness_ = DVect2(kn,ks);
484
485 //damping force
486 if (std::abs(vni_) <= limits<double>::epsilon()*1000.0)
487 dp_F_.rx() = 0.0;
488 else {
489 double ratio = vn/vni_;
490 double fac(0.0);
491 switch (dp_mode_) {
492 default:
493 {
494 double used = max(dp_en_,dp_enmin_);
495 fac = max(0.0,(1.0-used*used)/used)*ratio; //Gonthier et al.
496 }
497 break;
498 case 1:
499 fac = 0.75*(1.0-dp_en_*dp_en_)*ratio; // Lankarani et al (1989).
500 break;
501 case 2:
502 fac = 0.75*(1.0-dp_en_*dp_en_)*exp(2.0*(1-dp_en_))*ratio; // Zhiying et al.
503 break;
504 }
505
506 if (fac <= -1.0) { // sucking
507 if (pfac_ >= 0.0) { //switch in one timestep from pushing to sucking - instability
508 fac = 0.0;
509 pfac_ = 0.0;
510 } else
511 pfac_ = fac;
512 } else if (pfac_ != 0.0 && abs(fac) > 10.0*abs(pfac_)) { // growing too fast - instability
513 fac = 0.0;
514 pfac_ = 0.0;
515 } else
516 pfac_ = fac;
517
518 dp_F_.rx() = hz_F_.x()*fac;
519 }
520
521 DVect u_s = trans;
522 u_s.rx() = 0.0;
523 DVect vec = u_s * ks;
524 if (hz_mode_ && (hz_F_.x() < hz_F_old.x())) {
525 double rat = ks / ks_old;
526 fs_old *= rat;
527 }
528 DVect fs = fs_old - vec;
529
530 if (state->canFail_) {
531 // resolve sliding
532 double crit = hz_F_.x() * fric_;
533 double sfmag = fs.mag();
534 if (sfmag > crit) {
535 double rat = crit / sfmag;
536 fs *= rat;
537 if (!hz_slip_ && cmEvents_[fSlipChange] >= 0) {
538 auto c = state->getContact();
539 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
540 fish::Parameter() };
541 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
542 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
543 }
544 hz_slip_ = true;
545 } else {
546 if (hz_slip_) {
547 if (cmEvents_[fSlipChange] >= 0) {
548 auto c = state->getContact();
549 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
550 fish::Parameter((qint64)1) };
551 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
552 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
553 }
554 hz_slip_ = false;
555 }
556 }
557 }
558
559 fs.rx() = hz_F_.x();
560 hz_F_ = fs; // total force in hertz part
561
562 // 5) Compute energies
563 if (state->trackEnergy_) {
564 assert(energies_);
565 energies_->estrain_ = 0.0;
566 if (kn)
567 energies_->estrain_ = hz_alpha_*hz_F_.x()*hz_F_.x()/((hz_alpha_+1)*kn);
568 if (ks) {
569 DVect s = hz_F_;
570 s.rx() = 0.0;
571 double smag2 = s.mag2();
572 energies_->estrain_ += 0.5*smag2 / ks;
573 if (hz_slip_) {
574 hz_F_old.rx() = 0.0;
575 DVect avg_F_s = (s + hz_F_old)*0.5;
576 DVect u_s_el = (s - hz_F_old) / ks;
577 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
578 }
579 }
580 energies_->edashpot_ -= dp_F_ | trans;
581 }
582
583 } else {
584 hz_F_ = DVect(0.0);
585 dp_F_ = DVect(0.0);
586 pfac_ = 0.0;
587 }
588
589
590 return true;
591 }
592
593 void ContactModelHysteretic::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
594 // Only do something if the contact model is of the same type
595 if (old->getContactModel()->getName().compare("hysteretic",Qt::CaseInsensitive) == 0) {
596 ContactModelHysteretic *oldCm = (ContactModelHysteretic *)old;
597#ifdef THREED
598 // Need to rotate just the shear component from oldSystem to newSystem
599
600 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
601 DVect axis = oldSystem.e1() & newSystem.e1();
602 double c, ang, s;
603 DVect re2;
604 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
605 axis = axis.unit();
606 c = oldSystem.e1()|newSystem.e1();
607 if (c > 0)
608 c = std::min(c,1.0);
609 else
610 c = std::max(c,-1.0);
611 ang = acos(c);
612 s = sin(ang);
613 double t = 1. - c;
614 DMatrix<3,3> rm;
615 rm.get(0,0) = t*axis.x()*axis.x() + c;
616 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
617 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
618 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
619 rm.get(1,1) = t*axis.y()*axis.y() + c;
620 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
621 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
622 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
623 rm.get(2,2) = t*axis.z()*axis.z() + c;
624 re2 = rm*oldSystem.e2();
625 } else
626 re2 = oldSystem.e2();
627
628 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
629 axis = re2 & newSystem.e2();
630 DVect2 tpf;
631 DMatrix<2,2> m;
632 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
633 axis = axis.unit();
634 c = re2|newSystem.e2();
635 if (c > 0)
636 c = std::min(c,1.0);
637 else
638 c = std::max(c,-1.0);
639 ang = acos(c);
640 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
641 ang *= -1;
642 s = sin(ang);
643 m.get(0,0) = c;
644 m.get(1,0) = s;
645 m.get(0,1) = -m.get(1,0);
646 m.get(1,1) = m.get(0,0);
647 tpf = m*DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
648 } else {
649 m.get(0,0) = 1.;
650 m.get(0,1) = 0.;
651 m.get(1,0) = 0.;
652 m.get(1,1) = 1.;
653 tpf = DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
654 }
655 DVect pforce = DVect(0,tpf.x(),tpf.y());
656#else
657 oldSystem;
658 newSystem;
659 DVect pforce = DVect(0,oldCm->hz_F_.y());
660#endif
661 for (int i=1; i<dim; ++i)
662 hz_F_.rdof(i) += pforce.dof(i);
663 oldCm->hz_F_ = DVect(0.0);
664
665 if(oldCm->getEnergyActivated()) {
666 activateEnergy();
667 energies_->estrain_ = oldCm->energies_->estrain_;
668 energies_->eslip_ = oldCm->energies_->eslip_;
669 energies_->edashpot_ = oldCm->energies_->edashpot_;
670 oldCm->energies_->estrain_ = 0.0;
671 oldCm->energies_->eslip_ = 0.0;
672 oldCm->energies_->edashpot_ = 0.0;
673 }
674 }
675 }
676
677 void ContactModelHysteretic::setNonForcePropsFrom(IContactModel *old) {
678 // Only do something if the contact model is of the same type
679 if (old->getName().compare("hysteretic",Qt::CaseInsensitive) == 0 && !isBonded()) {
680 ContactModelHysteretic *oldCm = (ContactModelHysteretic *)old;
681 hn_ = oldCm->hn_;
682 hs_ = oldCm->hs_;
683 fric_ = oldCm->fric_;
684 vni_ = oldCm->vni_;
685 pfac_ = oldCm->pfac_;
686 }
687 }
688
689 DVect ContactModelHysteretic::getForce() const {
690 DVect ret(hz_F_);
691 ret += dp_F_;
692 return ret;
693 }
694
695 DAVect ContactModelHysteretic::getMomentOn1(const IContactMechanical *c) const {
696 DVect force = getForce();
697 DAVect ret(0.0);
698 c->updateResultingTorqueOn1Local(force,&ret);
699 return ret;
700 }
701
702 DAVect ContactModelHysteretic::getMomentOn2(const IContactMechanical *c) const {
703 DVect force = getForce();
704 DAVect ret(0.0);
705 c->updateResultingTorqueOn2Local(force,&ret);
706 return ret;
707 }
708
709 DAVect ContactModelHysteretic::getModelMomentOn1() const {
710 DAVect ret(0.0);
711 return ret;
712 }
713
714 DAVect ContactModelHysteretic::getModelMomentOn2() const {
715 DAVect ret(0.0);
716 return ret;
717 }
718
719 void ContactModelHysteretic::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
720 ret->clear();
721 ret->push_back({"hz_shear",ScalarInfo});
722 ret->push_back({"hz_poiss",ScalarInfo});
723 ret->push_back({"fric",ScalarInfo});
724 ret->push_back({"hz_force",VectorInfo});
725 ret->push_back({"hz_slip",ScalarInfo});
726 ret->push_back({"hz_mode",ScalarInfo});
727 ret->push_back({"hz_alpha",ScalarInfo});
728 ret->push_back({"dp_mode",ScalarInfo});
729 ret->push_back({"dp_en",ScalarInfo});
730 ret->push_back({"dp_enmin",ScalarInfo});
731 ret->push_back({"dp_force",VectorInfo});
732 }
733
734 void ContactModelHysteretic::objectPropValues(std::vector<double> *ret,const IContact *) const {
735 ret->clear();
736 ret->push_back(hz_shear_);
737 ret->push_back(hz_poiss_);
738 ret->push_back(fric_);
739 ret->push_back(hz_F_.mag());
740 ret->push_back(hz_slip_);
741 ret->push_back(hz_mode_);
742 ret->push_back(hz_alpha_);
743 ret->push_back(dp_mode_);
744 ret->push_back(dp_en_);
745 ret->push_back(dp_enmin_);
746 ret->push_back(dp_F_.mag());
747 }
748
749} // namespace cmodelsxd
750// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Dec 14, 2024 |
