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 QString getName() const { return "hertz"; }
43 virtual void setIndex(int i) { index_=i;}
44 virtual int getIndex() const {return index_;}
45
46 virtual QString 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 QString 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_ = NEWC(Energies());}
69 virtual bool getEnergyActivated() const {return (energies_ !=0);}
70
71 enum FishCallEvents { fActivated=0, fSlipChange};
72 virtual QString getFishCallEvents() const { return "contact_activated,slip_change"; }
73 virtual QVariant getProperty(uint32 i,const IContact *) const;
74 virtual bool getPropertyGlobal(uint32 i) const;
75 virtual bool setProperty(uint32 i,const QVariant &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 QString &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 NEWC(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<QString,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 = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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 QVariant ContactModelHertz::getProperty(uint32 i,const IContact *) const {
182 // The IContact pointer may be a nullptr!
183 QVariant 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 QVariant();
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 QVariant &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<uint32>())
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_ = NEWC(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_ = NEWC(dpProps());
292 dpProps_->dp_sratio_ = val;
293 return true;
294 }
295 case kwDpMode: {
296 if (!v.canConvert<int>())
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_ = NEWC(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_ = NEWC(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 QString gstr("hz_shear");
425 static const QString 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 QVariant vg1 = con->getEnd1()->getProperty(gstr);
433 QVariant vg2 = con->getEnd2()->getProperty(gstr);
434 QVariant vnu1 = con->getEnd1()->getProperty(nustr);
435 QVariant 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 QString fricstr("fric");
461 bool ContactModelHertz::updateEndFric(const IContactMechanical *con) {
462 assert(con);
463 QVariant v1 = con->getEnd1()->getProperty(fricstr);
464 QVariant 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 QString &name,const IContactMechanical *c) {
475 assert(c);
476 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
477 QRegExp rx(name,Qt::CaseInsensitive);
478 int idx = availableProperties.indexOf(rx)+1;
479 bool ret=false;
480
481 if (idx<=0)
482 return ret;
483
484 switch(idx) {
485 case kwHzShear: {
486 if (inheritanceField_ & shearMask)
487 ret = updateEndStiffCoef(c);
488 break;
489 }
490 case kwHzPoiss: {
491 if (inheritanceField_ & poissMask)
492 ret = updateEndStiffCoef(c);
493 break;
494 }
495 case kwFric: {
496 if (inheritanceField_ & fricMask)
497 ret = updateEndFric(c);
498 break;
499 }
500 }
501 return ret;
502 }
503
504 void ContactModelHertz::updateEffectiveStiffness(ContactModelMechanicalState *state) {
505 effectiveTranslationalStiffness_ = DVect2(hn_,hs_);
506 double overlap = rgap_ - state->gap_;
507 if (overlap <= 0.0) return;
508 double kn = hz_alpha_*hn_*pow(overlap,hz_alpha_-1.0);
509 double ks = hs_ * pow(hz_F_.x(),(hz_alpha_-1.0)/hz_alpha_);
510 DVect2 ret(kn,ks);
511 // correction if viscous damping active
512 if (dpProps_) {
513 DVect2 correct(1.0);
514 if (dpProps_->dp_nratio_)
515 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
516 if (dpProps_->dp_sratio_)
517 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
518 ret /= (correct*correct);
519 }
520 effectiveTranslationalStiffness_ = ret;
521 }
522
523 bool ContactModelHertz::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
524 assert(state);
525
526 if (state->activated()) {
527 if (cmEvents_[fActivated] >= 0) {
528 auto c = state->getContact();
529 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
530 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
531 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
532 }
533 }
534
535 double overlap = rgap_ - state->gap_;
536 DVect trans = state->relativeTranslationalIncrement_;
537#ifdef THREED
538 DVect norm(trans.x(),0.0,0.0);
539#else
540 DVect norm(trans.x(),0.0);
541#endif
542 //DAVect ang = state->relativeAngularIncrement_;
543 // normal force in Hertz part
544 double fn = hn_*pow(overlap,hz_alpha_);
545 // tangent normal stiffness
546 double kn = hz_alpha_ * hn_ * pow(overlap,hz_alpha_-1.0);
547 // initial tangent shear stiffness
548 double ks = hs_ * pow(fn,(hz_alpha_- 1.0)/hz_alpha_);
549
550 DVect fs_old = hz_F_;
551 fs_old.rx() = 0.0;
552
553 if (hz_mode_ && fn < hz_F_.x()) {
554 double ks_old = hs_ * pow(hz_F_.x(),(hz_alpha_- 1.0)/hz_alpha_);
555 double rat = ks/ks_old;
556 fs_old *= rat;
557 }
558
559 DVect u_s = trans;
560 u_s.rx() = 0.0;
561 DVect vec = u_s * ks;
562
563 DVect fs = fs_old - vec;
564
565 if (state->canFail_) {
566 // resolve sliding
567 double crit = fn * fric_;
568 double sfmag = fs.mag();
569 if (sfmag > crit) {
570 double rat = crit / sfmag;
571 fs *= rat;
572 if (!hz_slip_ && cmEvents_[fSlipChange] >= 0) {
573 auto c = state->getContact();
574 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
575 fish::Parameter((qint64)1) };
576 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
577 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
578 }
579 hz_slip_ = true;
580 } else {
581 if (hz_slip_) {
582 if (cmEvents_[fSlipChange] >= 0) {
583 auto c = state->getContact();
584 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
585 fish::Parameter((qint64)1) };
586 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
587 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
588 }
589 hz_slip_ = false;
590 }
591 }
592 }
593
594 hz_F_ = fs ; // total force in hertz part
595 hz_F_.rx() += fn;
596 effectiveTranslationalStiffness_ = DVect2(kn,ks);
597 // 3) Account for dashpot forces
598 if (dpProps_) {
599 dpProps_->dp_F_.fill(0.0);
600 double vcn(0.0), vcs(0.0);
601 setDampCoefficients(*state,&vcn,&vcs);
602 double fac = 1.0;
603 if (dpProps_->dp_alpha_ > 0.0) fac = pow(overlap,dpProps_->dp_alpha_);
604 // First damp all components
605 dpProps_->dp_F_ = u_s * (-1.0* vcs*fac) / timestep; // shear component
606 dpProps_->dp_F_ -= norm * vcn*fac / timestep; // normal component
607 // Need to change behavior based on the dp_mode
608 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
609 // limit the tensile if not bonded
610 if (dpProps_->dp_F_.x() + hz_F_.x() < 0)
611 dpProps_->dp_F_.rx() = - hz_F_.rx();
612 }
613 if (hz_slip_ && dpProps_->dp_mode_ > 1) {
614 // limit the shear if not sliding
615 double dfn = dpProps_->dp_F_.rx();
616 dpProps_->dp_F_.fill(0.0);
617 dpProps_->dp_F_.rx() = dfn;
618 }
619 // Correct effective translational stiffness
620 DVect2 correct(1.0);
621 if (dpProps_->dp_nratio_)
622 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
623 if (dpProps_->dp_sratio_)
624 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
625 effectiveTranslationalStiffness_ /= (correct*correct);
626 }
627
628 // 5) Compute energies
629 if (state->trackEnergy_) {
630 assert(energies_);
631 energies_->estrain_ = 0.0;
632 if (kn)
633 energies_->estrain_ = hz_alpha_*hz_F_.x()*hz_F_.x()/((hz_alpha_+1.0)*kn);
634 if (ks) {
635 double smag2 = fs.mag2();
636 energies_->estrain_ += 0.5*smag2 / ks;
637
638 if (hz_slip_) {
639 DVect avg_F_s = (fs + fs_old)*0.5;
640 DVect u_s_el = (fs - fs_old) / ks;
641 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
642 }
643 }
644 if (dpProps_) {
645 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
646 }
647 }
648
649 return true;
650 }
651
652 void ContactModelHertz::setForce(const DVect &v,IContact *c) {
653 hz_F(v);
654 if (v.x() > 0)
655 rgap_ = c->getGap() + pow(v.x()/hn_,1./hz_alpha_);
656 }
657
658 void ContactModelHertz::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
659 // Only do something if the contact model is of the same type
660 if (old->getContactModel()->getName().compare("hertz",Qt::CaseInsensitive) == 0) {
661 ContactModelHertz *oldCm = (ContactModelHertz *)old;
662#ifdef THREED
663 // Need to rotate just the shear component from oldSystem to newSystem
664
665 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
666 DVect axis = oldSystem.e1() & newSystem.e1();
667 double c, ang, s;
668 DVect re2;
669 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
670 axis = axis.unit();
671 c = oldSystem.e1()|newSystem.e1();
672 if (c > 0)
673 c = std::min(c,1.0);
674 else
675 c = std::max(c,-1.0);
676 ang = acos(c);
677 s = sin(ang);
678 double t = 1. - c;
679 DMatrix<3,3> rm;
680 rm.get(0,0) = t*axis.x()*axis.x() + c;
681 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
682 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
683 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
684 rm.get(1,1) = t*axis.y()*axis.y() + c;
685 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
686 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
687 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
688 rm.get(2,2) = t*axis.z()*axis.z() + c;
689 re2 = rm*oldSystem.e2();
690 } else
691 re2 = oldSystem.e2();
692
693 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
694 axis = re2 & newSystem.e2();
695 DVect2 tpf;
696 DMatrix<2,2> m;
697 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
698 axis = axis.unit();
699 c = re2|newSystem.e2();
700 if (c > 0)
701 c = std::min(c,1.0);
702 else
703 c = std::max(c,-1.0);
704 ang = acos(c);
705 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
706 ang *= -1;
707 s = sin(ang);
708 m.get(0,0) = c;
709 m.get(1,0) = s;
710 m.get(0,1) = -m.get(1,0);
711 m.get(1,1) = m.get(0,0);
712 tpf = m*DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
713 } else {
714 m.get(0,0) = 1.;
715 m.get(0,1) = 0.;
716 m.get(1,0) = 0.;
717 m.get(1,1) = 1.;
718 tpf = DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
719 }
720 DVect pforce = DVect(0,tpf.x(),tpf.y());
721#else
722 oldSystem;
723 newSystem;
724 DVect pforce = DVect(0,oldCm->hz_F_.y());
725#endif
726 for (int i=1; i<dim; ++i)
727 hz_F_.rdof(i) += pforce.dof(i);
728 oldCm->hz_F_ = DVect(0.0);
729 if (dpProps_ && oldCm->dpProps_) {
730#ifdef THREED
731 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
732 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
733#else
734 pforce = oldCm->dpProps_->dp_F_;
735#endif
736 dpProps_->dp_F_ += pforce;
737 oldCm->dpProps_->dp_F_ = DVect(0.0);
738 }
739
740 if(oldCm->getEnergyActivated()) {
741 activateEnergy();
742 energies_->estrain_ = oldCm->energies_->estrain_;
743 energies_->eslip_ = oldCm->energies_->eslip_;
744 energies_->edashpot_ = oldCm->energies_->edashpot_;
745 oldCm->energies_->estrain_ = 0.0;
746 oldCm->energies_->eslip_ = 0.0;
747 oldCm->energies_->edashpot_ = 0.0;
748 }
749 rgap_ = oldCm->rgap_;
750 }
751 }
752
753 void ContactModelHertz::setNonForcePropsFrom(IContactModel *old) {
754 // Only do something if the contact model is of the same type
755 if (old->getName().compare("hertz",Qt::CaseInsensitive) == 0 && !isBonded()) {
756 ContactModelHertz *oldCm = (ContactModelHertz *)old;
757 hn_ = oldCm->hn_;
758 hs_ = oldCm->hs_;
759 fric_ = oldCm->fric_;
760 rgap_ = oldCm->rgap_;
761
762 if (oldCm->dpProps_) {
763 if (!dpProps_)
764 dpProps_ = NEWC(dpProps());
765 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
766 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
767 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
768 }
769 }
770 }
771
772 DVect ContactModelHertz::getForce() const {
773 DVect ret(hz_F_);
774 if (dpProps_)
775 ret += dpProps_->dp_F_;
776 return ret;
777 }
778
779 DAVect ContactModelHertz::getMomentOn1(const IContactMechanical *c) const {
780 DVect force = getForce();
781 DAVect ret(0.0);
782 c->updateResultingTorqueOn1Local(force,&ret);
783 return ret;
784 }
785
786 DAVect ContactModelHertz::getMomentOn2(const IContactMechanical *c) const {
787 DVect force = getForce();
788 DAVect ret(0.0);
789 c->updateResultingTorqueOn2Local(force,&ret);
790 return ret;
791 }
792
793 DAVect ContactModelHertz::getModelMomentOn1() const {
794 DAVect ret(0.0);
795 return ret;
796 }
797
798 DAVect ContactModelHertz::getModelMomentOn2() const {
799 DAVect ret(0.0);
800 return ret;
801 }
802
803 void ContactModelHertz::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
804 ret->clear();
805 ret->push_back({"kwHzShear",ScalarInfo});
806 ret->push_back({"hz_poiss",ScalarInfo});
807 ret->push_back({"fric",ScalarInfo});
808 ret->push_back({"hz_alpha",ScalarInfo});
809 ret->push_back({"hz_slip",ScalarInfo});
810 ret->push_back({"hz_mode",ScalarInfo});
811 ret->push_back({"hz_force",VectorInfo});
812 ret->push_back({"dp_nratio",ScalarInfo});
813 ret->push_back({"dp_sratio",ScalarInfo});
814 ret->push_back({"dp_mode",ScalarInfo});
815 ret->push_back({"dp_force",VectorInfo});
816 ret->push_back({"dp_alpha",ScalarInfo});
817 ret->push_back({"rgap",ScalarInfo});
818 }
819
820 void ContactModelHertz::objectPropValues(std::vector<double> *ret,const IContact *) const {
821 ret->clear();
822 ret->push_back(hz_shear_);
823 ret->push_back(hz_poiss_);
824 ret->push_back(fric_);
825 ret->push_back(hz_alpha_);
826 ret->push_back(hz_slip_);
827 ret->push_back(hz_mode_);
828 ret->push_back(hz_F_.mag());
829 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0.0);
830 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0.0);
831 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
832 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
833 ret->push_back(dpProps_ ? dpProps_->dp_alpha_ : 0);
834 ret->push_back(rgap_);
835 }
836
837 void ContactModelHertz::setDampCoefficients(const ContactModelMechanicalState &state,double *vcn,double *vcs) {
838 double overlap = rgap_ - state.gap_;
839 double kn = hz_alpha_*hn_*pow(overlap,hz_alpha_-1.0);
840 double ks = hs_ * pow(hz_F_.x(),(hz_alpha_-1.0)/hz_alpha_);
841 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(state.inertialMass_*(kn));
842 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(state.inertialMass_*(ks));
843 }
844
845} // namespace cmodelsxd
846// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Oct 31, 2024 |
