Hysteretic Model Implementation

See this page for the documentation of this contact model.

contactmodelhysteretic.h

  1#pragma once
  2// contactmodelhysteretic.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef HYSTERETIC_LIB
  7#  define HYSTERETIC_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define HYSTERETIC_EXPORT
 10#else
 11#  define HYSTERETIC_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelHysteretic : public ContactModelMechanical {
 18    public:
 19        enum PropertyKeys { kwHzShear=1
 20                          , kwHzPoiss                            
 21                          , kwFric   
 22                          , kwHzF
 23                          , kwHzS
 24                          , kwHzSd
 25                          , kwHzAlpha
 26                          , kwDpMode
 27                          , kwDpEn
 28                          , kwDpEnMin
 29                          , kwDpF 
 30                         };
 31       
 32        HYSTERETIC_EXPORT ContactModelHysteretic();
 33        HYSTERETIC_EXPORT virtual ~ContactModelHysteretic();
 34        virtual void                copy(const ContactModel *c) override;
 35        virtual void                archive(ArchiveStream &); 
 36  
 37        virtual QString  getName() const { return "hysteretic"; }
 38        virtual void     setIndex(int i) { index_=i;}
 39        virtual int      getIndex() const {return index_;}
 40      
 41        virtual QString  getProperties() const { 
 42            return "hz_shear"
 43                   ",hz_poiss"
 44                   ",fric"
 45                   ",hz_force"
 46                   ",hz_slip"
 47                   ",hz_mode"
 48                   ",hz_alpha"
 49                   ",dp_mode"
 50                   ",dp_en"
 51                   ",dp_enmin"
 52                   ",dp_force"
 53            ;
 54        }
 55  
 56        enum EnergyKeys { kwEStrain=1,kwESlip,kwEDashpot};
 57        virtual QString  getEnergies() const { return "energy-strain,energy-slip,energy-dashpot";}
 58        virtual double   getEnergy(uint32 i) const;  // Base 1
 59        virtual bool     getEnergyAccumulate(uint32 i) const; // Base 1
 60        virtual void     setEnergy(uint32 i,const double &d); // Base 1
 61        virtual void     activateEnergy() { if (energies_) return; energies_ = NEWC(Energies());}
 62        virtual bool     getEnergyActivated() const {return (energies_ !=0);}
 63        
 64        enum FishCallEvents { fActivated=0, fSlipChange};
 65        virtual QString  getFishCallEvents() const { return "contact_activated,slip_change"; }
 66        virtual QVariant getProperty(uint32 i,const IContact *) const;
 67        virtual bool     getPropertyGlobal(uint32 i) const;
 68        virtual bool     setProperty(uint32 i,const QVariant &v,IContact *);
 69        virtual bool     getPropertyReadOnly(uint32 i) const;
 70        
 71        virtual bool     supportsInheritance(uint32 i) const; 
 72        virtual bool     getInheritance(uint32 i) const { assert(i<32); uint32 mask = to<uint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
 73        virtual void     setInheritance(uint32 i,bool b) { assert(i<32); uint32 mask = to<uint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
 74                
 75        virtual uint32     getMinorVersion() const;
 76        
 77        virtual bool    validate(ContactModelMechanicalState *state,const double &timestep);
 78        virtual bool    endPropertyUpdated(const QString &name,const IContactMechanical *c);
 79        virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep);
 80        virtual DVect2  getEffectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_;}
 81        virtual DAVect  getEffectiveRotationalStiffness() const { return DAVect(0.0);}
 82        
 83        virtual ContactModelHysteretic *clone() const override { return NEWC(ContactModelHysteretic()); }
 84        virtual double              getActivityDistance() const {return 0.0;}
 85        virtual bool                isOKToDelete() const { return !isBonded(); }
 86        virtual void                resetForcesAndMoments() { hz_F(DVect(0.0)); dp_F(DVect(0.0));  if (energies_) energies_->estrain_ = 0.0;}
 87        virtual void                setForce(const DVect &v,IContact *) { hz_F(v); }
 88        virtual void                setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }
 89        virtual double              getArea() const { return 0.0; }
 90
 91        virtual bool     checkActivity(const double &gap) { return gap <= 0.0; }
 92        
 93        virtual bool     isSliding() const { return hz_slip_; }
 94        virtual bool     isBonded() const { return false; }
 95        virtual void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes());
 96        virtual void     setNonForcePropsFrom(IContactModel *oldCM);
 97        
 98        const double & hz_shear() const {return hz_shear_;}
 99        void           hz_shear(const double &d) {hz_shear_=d;}
100        const double & hz_poiss() const {return hz_poiss_;}
101        void           hz_poiss(const double &d) {hz_poiss_=d;}
102        const double & fric() const {return fric_;}
103        void           fric(const double &d) {fric_=d;}
104        uint32           hz_mode() const {return hz_mode_;}
105        void           hz_mode(uint32 i) {hz_mode_=i;}
106        const DVect &  hz_F() const {return hz_F_;}
107        void           hz_F(const DVect &f) { hz_F_=f;}
108        bool           hz_S() const {return hz_slip_;}
109        void           hz_S(bool b) { hz_slip_=b;}
110        const double & hz_alpha() const {return hz_alpha_;}
111        void           hz_alpha(const double &d) {hz_alpha_=d;}
112        int            dp_mode() const {return dp_mode_;}
113        void           dp_mode(int i) {dp_mode_=i;}
114        const double & dp_en() const {return dp_en_;}
115        void           dp_en(const double &d) {dp_en_=d;}
116        const double & dp_enmin() const {return dp_enmin_;}
117        void           dp_enmin(const double &d) {dp_enmin_=d;}
118        const DVect &  dp_F() const {return dp_F_;}
119        void           dp_F(const DVect &f) { dp_F_=f;}
120        const double & hn() const {return hn_;}
121        void           hn(const double &d) {hn_=d;}
122        const double & hs() const {return hs_;}
123        void           hs(const double &d) {hs_=d;}
124        const double & vni() const {return vni_;}
125        void           vni(const double &d) {vni_=d;}
126        double         pfac() const {return pfac_;}
127        void           pfac(int i) {pfac_=i;}
128                
129        bool    hasEnergies() const {return energies_ ? true:false;}
130        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
131        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
132        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
133        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
134        double  edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
135        void    edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
136        
137        uint32 inheritanceField() const {return inheritanceField_;}
138        void inheritanceField(uint32 i) {inheritanceField_ = i;}
139        
140        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
141        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
142  
143        /// Return the total force that the contact model holds.
144        virtual DVect    getForce(const IContactMechanical *) const;
145
146        /// Return the total moment on 1 that the contact model holds
147        virtual DAVect   getMomentOn1(const IContactMechanical *) const;
148
149        /// Return the total moment on 1 that the contact model holds
150        virtual DAVect   getMomentOn2(const IContactMechanical *) const;
151
152    private:
153        static int index_;
154        
155        bool   updateStiffCoef(const IContactMechanical *con);
156        bool   updateEndStiffCoef(const IContactMechanical *con);
157        bool   updateEndFric(const IContactMechanical *con);
158        void   updateEffectiveStiffness(ContactModelMechanicalState *state);
159        // inheritance fields
160        uint32 inheritanceField_;
161        
162        // hertz model
163        double      hz_shear_;  // Shear modulus
164        double      hz_poiss_;  // Poisson ratio
165        double      fric_;       // Coulomb friction coefficient
166        DVect       hz_F_;      // Force carried in the hertz model
167        bool        hz_slip_;   // the current sliding state
168        uint32        hz_mode_;   // specifies down-scaling of the shear force when normal unloading occurs 
169        double      hz_alpha_;  // alpha exponent
170        int         dp_mode_;   // Damping scheme mode
171        double      dp_en_;     // normal restitution coefficient
172        double      dp_enmin_;  // minimal normal restitution coefficient in default mode
173        DVect       dp_F_;      // Damping Force
174                
175        // energies
176        struct Energies {
177            Energies() : estrain_(0.0), eslip_(0.0),edashpot_(0.0) {}
178            double estrain_;  // elastic energy stored in contact 
179            double eslip_;    // work dissipated by friction 
180            double edashpot_;    // work dissipated by dashpots
181        };
182        Energies *   energies_;    
183                
184        double hn_;                               // normal stiffness coefficient
185        double hs_;                               // shear stiffness coefficient
186        DVect2 effectiveTranslationalStiffness_;  // effective stiffness
187        double vni_;                              // normal impact velocity
188        double pfac_;                             // previous damping factor
189
190    };
191
192} // namespace cmodelsxd
193// EoF

Top

contactmodelhysteretic.cpp

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

Top