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