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