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_ = NEW 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 &timestep);
 82        virtual bool    endPropertyUpdated(const QString &name,const IContactMechanical *c);
 83        virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep);
 84        virtual DVect2  getEffectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_;}
 85        virtual DAVect  getEffectiveRotationalStiffness() const { return DAVect(0.0);}
 86        
 87        virtual ContactModelHertz *clone() const override { return NEW 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

Top

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 = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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        QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
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 &timestep) {
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_ = NEW 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

Top