Burger’s Contact Model Implementation

See this page for the documentation of this contact model.

contactmodelburger.h

  1#pragma once
  2// contactmodelburger.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef BURGER_LIB
  7#  define BURGER_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define BURGER_EXPORT
 10#else
 11#  define BURGER_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelBurger : public ContactModelMechanical {
 18    public:
 19        // Constructor: Set default values for contact model properties.
 20        BURGER_EXPORT ContactModelBurger();
 21        // Destructor, called when contact is deleted: free allocated memory, etc.
 22        BURGER_EXPORT virtual ~ContactModelBurger();
 23        // Contact model name (used as keyword for commands and FISH).
 24        virtual QString  getName() const { return "burger"; }
 25        // The index provides a quick way to determine the type of contact model.
 26        // Each type of contact model in PFC must have a unique index; this is assigned
 27        // by PFC when the contact model is loaded. This index should be set to -1
 28        virtual void     setIndex(int i) { index_=i;}
 29        virtual int      getIndex() const {return index_;}
 30        // Contact model version number (e.g., MyModel05_1). The version number can be
 31        // accessed during the save-restore operation (within the archive method,
 32        // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for 
 33        // future modifications to the contact model data structure.
 34        virtual uint32     getMinorVersion() const;
 35        // Copy the state information to a newly created contact model.
 36        // Provide access to state information, for use by copy method.
 37        virtual void     copy(const ContactModel *c) override;
 38        // Provide save-restore capability for the state information.
 39        virtual void     archive(ArchiveStream &); 
 40        // Enumerator for the properties.
 41        enum PropertyKeys { 
 42              kwKnK=1
 43            , kwCnK                            
 44            , kwKnM
 45            , kwCnM                            
 46            , kwKsK                            
 47            , kwCsK                            
 48            , kwKsM
 49            , kwCsM                            
 50            , kwMode
 51            , kwFric   
 52            , kwF
 53            , kwS
 54        };
 55        // Contact model property names in a comma separated list. The order corresponds with
 56        // the order of the PropertyKeys enumerator above. One can visualize any of these 
 57        // properties in PFC automatically. 
 58        virtual QString  getProperties() const { 
 59            return "bur_knk"
 60                   ",bur_cnk"
 61                   ",bur_knm"
 62                   ",bur_cnm"
 63                   ",bur_ksk"
 64                   ",bur_csk"
 65                   ",bur_ksm"
 66                   ",bur_csm"
 67                   ",bur_mode"
 68                   ",bur_fric"
 69                   ",bur_force"
 70                   ",bur_slip";
 71        }
 72        // Enumerator for contact model related FISH callback events. 
 73        enum FishCallEvents {
 74            fActivated=0
 75            ,fSlipChange
 76        };
 77        // Contact model FISH callback event names in a comma separated list. The order corresponds with
 78        // the order of the FishCallEvents enumerator above. 
 79        virtual QString  getFishCallEvents() const { 
 80            return 
 81                "contact_activated"
 82                ",slip_change"; 
 83        }
 84
 85        // Return the specified contact model property.
 86        virtual QVariant getProperty(uint32 i,const IContact *) const;
 87        // The return value denotes whether or not the property corresponds to the global
 88        // or local coordinate system (TRUE: global system, FALSE: local system). The
 89        // local system is the contact-plane system (nst) defined as follows.
 90        // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
 91        // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
 92        // and tc are unit vectors in directions of the nst axes.
 93        // This is used when rendering contact model properties that are vectors.
 94        virtual bool     getPropertyGlobal(uint32 i) const;
 95        // Set the specified contact model property, ensuring that it is of the correct type
 96        // and within the correct range --- if not, then throw an exception.
 97        // The return value denotes whether or not the update has affected the timestep
 98        // computation (by having modified the translational or rotational tangent stiffnesses).
 99        // If true is returned, then the timestep will be recomputed.
100        virtual bool     setProperty(uint32 i,const QVariant &v,IContact *);
101        // The return value denotes whether or not the property is read-only
102        // (TRUE: read-only, FALSE: read-write).
103        virtual bool     getPropertyReadOnly(uint32 i) const;
104
105        // Prepare for entry into ForceDispLaw. The validate function is called when:
106        // (1) the contact is created, (2) a property of the contact that returns a true via
107        // the setProperty method has been modified and (3) when a set of cycles is executed
108        // via the {cycle N} command.
109        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
110        virtual bool    validate(ContactModelMechanicalState *state,const double &timestep);
111        // The forceDisplacementLaw function is called during each cycle. Given the relative
112        // motion of the two contacting pieces (via
113        //   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
114        //   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
115        //     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
116        //     Ddss : relative  shear-displacement increment (s-axis component)
117        //     Ddst : relative  shear-displacement increment (t-axis component)
118        //     Dtt  : relative twist-rotation increment
119        //     Dtbs : relative  bend-rotation increment (s-axis component)
120        //     Dtbt : relative  bend-rotation increment (t-axis component)
121        //       The relative displacement and rotation increments:
122        //         Dd = Ddn*nc + Ddss*sc + Ddst*tc
123        //         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
124        //       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
125        //       [see {Table 1: Contact State Variables} in PFC Model Components:
126        //       Contacts and Contact Models: Contact Resolution]
127        // ) and the contact properties, this function must update the contact force and
128        // moment.
129        //   The force_ is acting on piece 2, and is expressed in the local coordinate system
130        //   (defined in getPropertyGlobal) such that the first component positive denotes
131        //   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
132        //   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
133        //   get the total moment. 
134        // The return value indicates the contact activity status (TRUE: active, FALSE:
135        // inactive) during the next cycle.
136        // Additional information:
137        //   * If state->activated() is true, then the contact has just become active (it was
138        //     inactive during the previous time step).
139        //   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
140        //     the forceDispLaw handle the case of {state->canFail_ == true}.
141        virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep);
142        // The getEffectiveXStiffness functions return the translational and rotational
143        // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
144        // the translational tangent shear stiffness is zero (but this stiffness reduction
145        // is typically ignored when computing a stable time step). If the contact model
146        // includes a dashpot, then the translational stiffnesses must be increased (see
147        // Potyondy (2009)).
148        //   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
149        //   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
150        //   December 7, 2009.]
151        virtual DVect2  getEffectiveTranslationalStiffness() const { return DVect2(knk_,ksk_); }
152
153        // Return a new instance of the contact model. This is used in the CMAT
154        // when a new contact is created. 
155        virtual ContactModelBurger *clone() const override { return NEW ContactModelBurger(); }
156        // The getActivityDistance function is called by the contact-resolution logic when
157        // the CMAT is modified. Return value is the activity distance used by the
158        // checkActivity function.
159        virtual double              getActivityDistance() const {return 0.0;}
160        // The isOKToDelete function is called by the contact-resolution logic when...
161        // Return value indicates whether or not the contact may be deleted.
162        // If TRUE, then the contact may be deleted when it is inactive.
163        // If FALSE, then the contact may not be deleted (under any condition).
164        virtual bool                isOKToDelete() const { return !isBonded(); }
165        // Zero the forces and moments stored in the contact model. This function is called
166        // when the contact becomes inactive.
167        virtual void                resetForcesAndMoments() { force_.fill(0.0);}
168        virtual void                setForce(const DVect &v,IContact *) { force_ = v; }
169        virtual void                setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }
170        virtual double              getArea() const { return 0.0; }
171        // The checkActivity function is called by the contact-resolution logic when...
172        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
173        virtual bool     checkActivity(const double &gap) { return  gap <= 0.0; }
174
175        // Returns the sliding state (FALSE is returned if not implemented).
176        virtual bool     isSliding() const { return s_; }
177        // Returns the bonding state (FALSE is returned if not implemented).
178        virtual bool     isBonded() const { return false; }
179
180        virtual bool    endPropertyUpdated(const QString &,const IContactMechanical *) {return false;}
181       // Both of these methods are called only for contacts with facets where the wall 
182        // resolution scheme is set the full. In such cases one might wish to propagate 
183        // contact state information (e.g., shear force) from one active contact to another. 
184        // See the Faceted Wall section in the documentation. 
185        virtual void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes());
186        virtual void     setNonForcePropsFrom(IContactModel *oldCM);
187
188        /// Return the total force that the contact model holds.
189        virtual DVect    getForce(const IContactMechanical *) const;
190
191        /// Return the total moment on 1 that the contact model holds
192        virtual DAVect   getMomentOn1(const IContactMechanical *) const;
193
194        /// Return the total moment on 1 that the contact model holds
195        virtual DAVect   getMomentOn2(const IContactMechanical *) const;
196
197   
198        // Methods to get and set properties. 
199        const double & knk() const {return knk_;}
200        void           knk(const double &d) {knk_=d;}
201        const double & cnk() const {return cnk_;}
202        void           cnk(const double &d) {cnk_=d;}
203        const double & knm() const {return knm_;}
204        void           knm(const double &d) {knm_=d;}
205        const double & cnm() const {return cnm_;}
206        void           cnm(const double &d) {cnm_=d;}
207        const double & ksk() const {return ksk_;}
208        void           ksk(const double &d) {ksk_=d;}
209        const double & csk() const {return csk_;}
210        void           csk(const double &d) {csk_=d;}
211        const double & ksm() const {return ksm_;}
212        void           ksm(const double &d) {ksm_=d;}
213        const double & csm() const {return csm_;}
214        void           csm(const double &d) {csm_=d;}
215
216        const double & fric() const { return fric_;}
217        void           fric(const double &d) {fric_=d;}
218        uint32           bmode() const {return bmode_;}
219        void           bmode(uint32 b) {bmode_= b;}
220
221        const double & fn0() const { return fn0_;}
222        void           fn0(const double &d) {fn0_=d;}
223        const double & u_n0() const { return u_n0_;}
224        void           u_n0(const double &d) {u_n0_=d;}
225        const double & u_nk0() const { return u_nk0_;}
226        void           u_nk0(const double &d) {u_nk0_=d;}
227        const DVect &  u_sk() const { return u_sk_;}
228        void           u_sk(const DVect &v) {u_sk_=v;}
229
230        const double & conAn() const { return conAn_;}
231        void           conAn(const double &d) {conAn_=d;}
232        const double & conB_An() const { return conB_An_;}
233        void           conB_An(const double &d) {conB_An_=d;}
234        const double & conCn() const { return conCn_;}
235        void           conCn(const double &d) {conCn_=d;}
236        const double & conDn() const { return conDn_;}
237        void           conDn(const double &d) {conDn_=d;}
238
239        const double & conAs() const { return conAs_;}
240        void           conAs(const double &d) {conAs_=d;}
241        const double & conB_As() const { return conB_As_;}
242        void           conB_As(const double &d) {conB_As_=d;}
243        const double & conCs() const { return conCs_;}
244        void           conCs(const double &d) {conCs_=d;}
245        const double & conDs() const { return conDs_;}
246        void           conDs(const double &d) {conDs_=d;}
247
248        const double & tdel() const { return tdel_;}
249        void           tdel(const double &d) {tdel_=d;}
250        const DVect &  force() const { return force_;}
251        void           force(const DVect &v) {force_=v;}
252        bool           s() const {return s_;}
253        void           s(bool b) {s_= b;}
254
255    private:
256        // Index - used internally by PFC. Should be set to -1 in the cpp file. 
257        static int index_;
258
259        // Burger model properties
260        double  knk_;     // normal stiffness for Kelvin section  (bur_knk)
261        double  cnk_;     // normal viscosity for Kelvin section  (bur_cnk)
262        double  knm_;     // normal stiffness for Maxwell section (bur_knm)
263        double  cnm_;     // normal viscosity for Maxwell section (bur_cnm)
264        double  ksk_;     // shear stiffness for Kelvin section   (bur_ksk)
265        double  csk_;     // shear viscosity for Kelvin section   (bur_csk)
266        double  ksm_;     // shear stiffness for Maxwell section  (bur_ksm)
267        double  csm_;     // shear viscosity for Maxwell section  (bur_csm)
268        double  fric_;    // friction coefficient                 (bur_fric)
269        uint32    bmode_;   // FDLaw option, with or without tensile force, default false (with tensile)
270        double  fn0_;     // normal contact force 1 step before
271        double  u_n0_;    // normal total overlap 1 step before
272        double  u_nk0_;   // normal overlap of Kelvin part 1step before
273        DVect   u_sk_;    // shear relative displacement of Kelvin part 1step before
274        double  conAn_;   // constant A in eq.(), normal
275        double  conB_An_; // constant B/A in eq.(), normal
276        double  conCn_;   // constant C in eq.(), normal
277        double  conDn_;   // constant D in eq.(), normal
278        double  conAs_;   // constant A in eq.(), shear
279        double  conB_As_; // constant B/A in eq.(), shear
280        double  conCs_;   // constant C in eq.(), shear
281        double  conDs_;   // constant D in eq.(), shear
282        double  tdel_;    // current timestep
283        DVect   force_;   // current total force
284        bool    s_;       // current sliding state
285
286        // Constants
287        inline double A(const double k_k, const double c_k) {    return(1.0 + k_k*tdel_/(2.0*c_k)); }
288        inline double B(const double k_k, const double c_k) { return(1.0 - k_k*tdel_/(2.0*c_k)); }
289        inline double B_A(const double k_k, const double c_k) { return(B(k_k,c_k)/A(k_k,c_k)); }
290        inline double C(const double k_k, const double c_k, const double k_m, const double c_m) {
291          return(tdel_/(2.0*c_k*A(k_k,c_k)) + 1.0/k_m + tdel_/(2.0*c_m)); }
292        inline double D(const double k_k, const double c_k, const double k_m, const double c_m) {
293          return(tdel_/(2.0*c_k*A(k_k,c_k)) - 1.0/k_m + tdel_/(2.0*c_m)); }
294
295    };
296} // namespace cmodelsxd
297// EoF

Top

contactmodelburger.cpp

  1// contactmodelburger.cpp
  2#include "contactmodelburger.h"
  3
  4#include "module/interface/icontactmechanical.h"
  5#include "module/interface/icontact.h"
  6#include "module/interface/ipiecemechanical.h"
  7#include "module/interface/ipiece.h"
  8#include "module/interface/ifishcalllist.h"
  9#include "fish/src/parameter.h"
 10
 11#include "../version.txt"
 12
 13#include "utility/src/tptr.h"
 14#include "shared/src/mathutil.h"
 15
 16#include "kernel/interface/iprogram.h"
 17#include "module/interface/icontactthermal.h"
 18#include "contactmodel/src/contactmodelthermal.h"
 19
 20#ifdef BURGER_LIB
 21#ifdef _WIN32
 22    int __stdcall DllMain(void *,unsigned, void *) {
 23        return 1;
 24    }
 25#endif
 26
 27    extern "C" EXPORT_TAG const char *getName() {
 28#if DIM==3
 29        return "contactmodelmechanical3dburger";
 30#else
 31        return "contactmodelmechanical2dburger";
 32#endif
 33    }
 34
 35    extern "C" EXPORT_TAG unsigned getMajorVersion() {
 36        return MAJOR_VERSION;
 37    }
 38
 39    extern "C" EXPORT_TAG unsigned getMinorVersion() {
 40        return MINOR_VERSION;
 41    }
 42
 43    extern "C" EXPORT_TAG void *createInstance() {
 44        cmodelsxd::ContactModelBurger *m = NEW cmodelsxd::ContactModelBurger();
 45        return (void *)m;
 46    }
 47#endif 
 48
 49namespace cmodelsxd {
 50    using namespace itasca;
 51
 52    int ContactModelBurger::index_ = -1;
 53    uint32 ContactModelBurger::getMinorVersion() const { return MINOR_VERSION; }
 54
 55    ContactModelBurger::ContactModelBurger() :knk_(0.0)      
 56                                             ,cnk_(0.0)      
 57                                             ,knm_(0.0)     
 58                                             ,cnm_(0.0)     
 59                                             ,ksk_(0.0)      
 60                                             ,csk_(0.0)   
 61                                             ,ksm_(0.0)  
 62                                             ,csm_(0.0)  
 63                                             ,fric_(0.0)  
 64                                             ,bmode_(0)  
 65                                             ,fn0_(0.0)  
 66                                             ,u_n0_(0.0)  
 67                                             ,u_nk0_(0.0)  
 68                                             ,u_sk_(0.0) 
 69                                             ,conAn_(0.0) 
 70                                             ,conB_An_(0.0)
 71                                             ,conCn_(0.0)   
 72                                             ,conDn_(0.0)   
 73                                             ,conAs_(0.0) 
 74                                             ,conB_As_(0.0)
 75                                             ,conCs_(0.0)  
 76                                             ,conDs_(0.0) 
 77                                             ,tdel_(0.0) 
 78                                             ,force_(0.0) 
 79                                             ,s_(false) 
 80    {  
 81
 82    }
 83
 84    ContactModelBurger::~ContactModelBurger() {
 85    }
 86
 87    void ContactModelBurger::archive(ArchiveStream &stream) {
 88        // The stream allows one to archive the values of the contact model
 89        // so that it can be saved and restored. The minor version can be
 90        // used here to allow for incremental changes to the contact model too. 
 91        stream & knk_;     
 92        stream & cnk_;     
 93        stream & knm_;     
 94        stream & cnm_;     
 95        stream & ksk_;     
 96        stream & csk_;     
 97        stream & ksm_;     
 98        stream & csm_;       
 99        stream & fric_;    
100        stream & bmode_;   
101        stream & fn0_;     
102        stream & u_n0_;    
103        stream & u_nk0_;   
104        stream & u_sk_;    
105        stream & conAn_;   
106        stream & conB_An_; 
107        stream & conCn_;   
108        stream & conDn_;   
109        stream & conAs_;   
110        stream & conB_As_; 
111        stream & conCs_;   
112        stream & conDs_;   
113        stream & tdel_;    
114        stream & force_;    
115        stream & s_;    
116    }
117
118    void ContactModelBurger::copy(const ContactModel *cm) {
119        // Copy all of the contact model properties. Used in the CMAT 
120        // when a new contact is created. 
121        ContactModelMechanical::copy(cm);
122        const ContactModelBurger *in = dynamic_cast<const ContactModelBurger*>(cm);
123        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
124        
125        knk(in->knk());     
126        cnk(in->cnk());     
127        knm(in->knm());     
128        cnm(in->cnm());     
129        ksk(in->ksk());     
130        csk(in->csk());     
131        ksm(in->ksm());     
132        csm(in->csm());     
133        fric(in->fric());    
134        bmode(in->bmode());   
135        fn0(in->fn0());     
136        u_n0(in->u_n0());    
137        u_nk0(in->u_nk0());   
138        u_sk(in->u_sk());    
139        conAn(in->conAn());   
140        conB_An(in->conB_An()); 
141        conCn(in->conCn());   
142        conDn(in->conDn());   
143        conAs(in->conAs());   
144        conB_As(in->conB_As()); 
145        conCs(in->conCs());   
146        conDs(in->conDs());   
147        tdel(in->tdel()); 
148        force(in->force());
149        s(in->s());
150    }
151
152
153    QVariant ContactModelBurger::getProperty(uint32 i,const IContact *) const {
154        // Return the property. The IContact pointer is provided so that 
155        // more complicated properties, depending on contact characteristics,
156        // can be calculated. Not used with the Burger model.
157        QVariant var;
158        switch (i) {
159        case kwKnK   : return knk_;
160        case kwCnK   : return cnk_;
161        case kwKnM   : return knm_;
162        case kwCnM   : return cnm_;
163        case kwKsK   : return ksk_;
164        case kwCsK   : return csk_;
165        case kwKsM   : return ksm_;
166        case kwCsM   : return csm_;
167        case kwMode  : return bmode_;
168        case kwFric  : return fric_;
169        case kwF     : var.setValue(force_); return var;
170        case kwS     : return s_;
171        }
172        assert(0);
173        return QVariant();
174    }
175
176    bool ContactModelBurger::getPropertyGlobal(uint32 i) const {
177        // Returns whether or not a property is held in the global axis system (TRUE)
178        // or the local system (FALSE). Used by the plotting logic.
179        switch (i) {
180        case kwF:   
181            return false;
182        }
183        return true;
184    }
185
186    bool ContactModelBurger::setProperty(uint32 i,const QVariant &v,IContact *) {
187        // Set a contact model property. Return value indicates that the timestep
188        // should be recalculated. 
189        switch (i) {
190        case kwKnK: {
191                if (!v.canConvert<double>())
192                    throw Exception("bur_knk must be a double.");
193                double val(v.toDouble());
194                if (val<0.0)
195                    throw Exception("Negative bur_knk not allowed.");
196                knk_ = val;
197                return true;
198            }
199        case kwCnK: {
200                if (!v.canConvert<double>())
201                    throw Exception("bur_cnk must be a double.");
202                double val(v.toDouble());
203                if (val<0.0)
204                    throw Exception("Negative bur_cnk not allowed.");
205                cnk_ = val;
206                return true;
207            }
208        case kwKnM: {
209                if (!v.canConvert<double>())
210                    throw Exception("bur_knm must be a double.");
211                double val(v.toDouble());
212                if (val<0.0)
213                    throw Exception("Negative bur_knm not allowed.");
214                knm_ = val;
215                return true;
216            }
217        case kwCnM: {
218                if (!v.canConvert<double>())
219                    throw Exception("bur_cnm must be a double.");
220                double val(v.toDouble());
221                if (val<0.0)
222                    throw Exception("Negative bur_cnm not allowed.");
223                cnm_ = val;
224                return true;
225            }
226        case kwKsK: {
227                if (!v.canConvert<double>())
228                    throw Exception("bur_ksk must be a double.");
229                double val(v.toDouble());
230                if (val<0.0)
231                    throw Exception("Negative bur_ksk not allowed.");
232                ksk_ = val;
233                return true;
234            }
235        case kwCsK: {
236                if (!v.canConvert<double>())
237                    throw Exception("bur_csk must be a double.");
238                double val(v.toDouble());
239                if (val<0.0)
240                    throw Exception("Negative bur_csk not allowed.");
241                csk_ = val;
242                return true;
243            }
244        case kwKsM: {
245                if (!v.canConvert<double>())
246                    throw Exception("bur_ksm must be a double.");
247                double val(v.toDouble());
248                if (val<0.0)
249                    throw Exception("Negative bur_ksm not allowed.");
250                ksm_ = val;
251                return true;
252            }
253        case kwCsM: {
254                if (!v.canConvert<double>())
255                    throw Exception("bur_csm must be a double.");
256                double val(v.toDouble());
257                if (val<0.0)
258                    throw Exception("Negative bur_csm not allowed.");
259                csm_ = val;
260                return true;
261            }
262        case kwMode: {
263                if (!v.canConvert<uint32>())
264                    throw Exception("bur_mode must be 0 (tensile) or 1 (no-tension).");
265                uint32 val(v.toUInt());
266                if (val >1)
267                    throw Exception("bur_mode must be 0 (tensile) or 1 (no-tension).");
268                bmode_ = val;
269                return false;
270            }
271        case kwFric: {
272                if (!v.canConvert<double>())
273                    throw Exception("bur_fric must be a double.");
274                double val(v.toDouble());
275                if (val<0.0)
276                    throw Exception("Negative bur_fric not allowed.");
277                fric_ = val;  
278                return false;
279            }
280        }
281        return false;
282    }
283
284    bool ContactModelBurger::getPropertyReadOnly(uint32 i) const {
285        // Returns TRUE if a property is read only or FALSE otherwise. 
286        switch (i) {
287        case kwF:
288        case kwS:
289            return true;
290        default:
291            break;
292        }
293        return false;
294    }
295
296    bool ContactModelBurger::validate(ContactModelMechanicalState *state,const double &) {
297        // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
298        // (1) the contact is created, (2) a property of the contact that returns a true via
299        // the setProperty method has been modified and (3) when a set of cycles is executed
300        // via the {cycle N} command.
301        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
302        return checkActivity(state->gap_);
303    }
304     
305    bool ContactModelBurger::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
306        assert(state);
307
308        // Current overlap
309        double overlap = - 1.0*state->gap_;
310        // Relative translational increment
311        DVect trans = state->relativeTranslationalIncrement_;
312        // Correction factor to account for when the contact becomes newly active.
313        // We estimate the time of activity during the timestep when the contact has first 
314        // become active and scale the forces accordingly.
315        double correction = 1.0;
316        // The contact was just activated from an inactive state
317        if (state->activated()) {
318            // Trigger the FISH callback if one is hooked up to the 
319            // contact_activated event.
320            if (cmEvents_[fActivated] >= 0) {
321                // An FArray of QVariant is returned and these will be passed
322                // to the FISH function as an array of FISH symbols as the second
323                // argument to the FISH callback function. 
324                auto c = state->getContact();
325                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
326                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
327                fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
328            }
329            // Calculate the correction factor.
330            if (trans.x()) {
331                correction = -1.0*overlap / trans.x();
332                if (correction < 0)
333                    correction = 1.0;
334            }
335        }
336
337        trans*=correction;
338
339        if (timestep!=tdel_) { // re-calculated constants.
340            tdel_ = timestep;
341            // need some protection for divided by zero (k_k c_k k_m c_m = zero)
342            conAn_ = A(knk_, cnk_);
343            double conBn = B(knk_, cnk_);
344            conB_An_ = conBn / conAn_;
345            conCn_ = C(knk_, cnk_, knm_, cnm_);
346            conDn_ = D(knk_, cnk_, knm_, cnm_);
347            conAs_ = A(ksk_, csk_);
348            double conBs = B(ksk_, csk_);
349            conB_As_ = conBs / conAs_;
350            conCs_ = C(ksk_, csk_, ksm_, csm_);
351            conDs_ = D(ksk_, csk_, ksm_, csm_);
352        }
353
354
355        // normal force 
356        force_.rx() = 1.0/conCn_*(overlap-u_n0_+(1.0-conB_An_)*u_nk0_-conDn_*fn0_);
357        if (bmode_ && force_.x()<0.0) force_.rx() = 0.0;
358        u_nk0_ = conB_An_*u_nk0_+timestep/(2.0*cnk_*conAn_)*(force_.x()+fn0_);
359        u_n0_  = overlap;
360        fn0_   = force_.x();
361
362        // Calculate the shear force.
363        DVect sforce(0.0);
364        DVect sforce_old = force_;
365        sforce_old.rx()=0.0;
366        DVect v1 = trans;
367        DVect v2 = u_sk_ * (1.0-conB_As_);
368        DVect v3 = sforce_old * conDs_;
369        sforce = (v1+v2+v3) / conCs_ * (-1.0);
370
371        double d1 = timestep / (2.0*csk_*conAs_);
372        sforce.rx() = 0.0;
373        v1 = sforce + sforce_old;
374        u_sk_ = u_sk_*conB_As_-v1*d1;
375
376        // The canFail flag corresponds to whether or not the contact can undergo non-linear
377        // force-displacement response. If the SOLVE ELASTIC command is given then the 
378        // canFail state is set to FALSE. Otherwise it is always TRUE. 
379        if (state->canFail_) {
380            // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
381            double crit = force_.x() * fric_;
382            crit = max(0.0,crit);
383            // This is the magnitude of the shear force.
384            double sfmag = sforce.mag();
385            // Sliding occurs when the magnitude of the shear force is greater than the 
386            // critical value.
387            if (sfmag > crit) {
388                // Lower the shear force to the critical value for sliding.
389                double rat = crit / sfmag;
390                sforce *= rat;
391                // Handle the slip_change event if one has been hooked up. Sliding has commenced.  
392                if (!s_ && cmEvents_[fSlipChange] >= 0) {
393                    auto c = state->getContact();
394                    std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
395                                                         fish::Parameter()                  };
396                    IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
397                    fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
398                }
399                s_ = true;
400            } else {
401                // Handle the slip_change event if one has been hooked up and
402                // the contact was previously sliding. Sliding has ceased.  
403                if (s_) {
404                    if (cmEvents_[fSlipChange] >= 0) {
405                        auto c = state->getContact();
406                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
407                                                             fish::Parameter((qint64)1)      };
408                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
409                        fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
410                    }
411                    s_ = false;
412                }
413            }
414        }
415        
416        // Set the shear components of the total force.
417        for (int i=1; i<dim; ++i)
418            force_.rdof(i) = sforce.dof(i);
419        // This is just a sanity check to ensure, in debug mode, that the force isn't wonky. 
420        assert(force_ == force_);
421        return true;
422    }
423
424    void ContactModelBurger::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
425        // Only called for contacts with wall facets when the wall resolution scheme
426        // is set to full!
427        // Only do something if the contact model is of the same type
428        if (old->getContactModel()->getName().compare("burger",Qt::CaseInsensitive) == 0 && !isBonded()) {
429            ContactModelBurger *oldCm = (ContactModelBurger *)old;
430#ifdef THREED
431            // Need to rotate just the shear component from oldSystem to newSystem
432
433            // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
434            DVect axis = oldSystem.e1() & newSystem.e1();
435            double c, ang, s;
436            DVect re2;
437            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
438                axis = axis.unit();
439                c = oldSystem.e1()|newSystem.e1();
440                if (c > 0)
441                    c = std::min(c,1.0);
442                else
443                    c = std::max(c,-1.0);
444                ang = acos(c);
445                s = sin(ang);
446                double t = 1. - c;
447                DMatrix<3,3> rm;
448                rm.get(0,0) = t*axis.x()*axis.x() + c;
449                rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
450                rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
451                rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
452                rm.get(1,1) = t*axis.y()*axis.y() + c;
453                rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
454                rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
455                rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
456                rm.get(2,2) = t*axis.z()*axis.z() + c;
457                re2 = rm*oldSystem.e2();
458            }
459            else
460                re2 = oldSystem.e2();
461            // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
462            axis = re2 & newSystem.e2();
463            DVect2 tpf;
464            DMatrix<2,2> m;
465            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
466                axis = axis.unit();
467                c = re2|newSystem.e2();
468                if (c > 0)
469                    c = std::min(c,1.0);
470                else
471                    c = std::max(c,-1.0);
472                ang = acos(c);
473                if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
474                    ang *= -1;
475                s = sin(ang);
476                m.get(0,0) = c;
477                m.get(1,0) = s;
478                m.get(0,1) = -m.get(1,0);
479                m.get(1,1) = m.get(0,0);
480                tpf = m*DVect2(oldCm->force_.y(),oldCm->force_.z());
481            } else {
482                m.get(0,0) = 1.;
483                m.get(0,1) = 0.;
484                m.get(1,0) = 0.;
485                m.get(1,1) = 1.;
486                tpf = DVect2(oldCm->force_.y(),oldCm->force_.z());
487            }
488            DVect pforce = DVect(0,tpf.x(),tpf.y());
489#else
490            oldSystem;
491            newSystem;
492            DVect pforce = DVect(0,oldCm->force_.y());
493#endif
494            for (int i=1; i<dim; ++i)
495                force_.rdof(i) += pforce.dof(i);
496            oldCm->force_ = DVect(0.0);
497        }
498        assert(force_ == force_);
499    }
500
501    void ContactModelBurger::setNonForcePropsFrom(IContactModel *old) {
502        // Only called for contacts with wall facets when the wall resolution scheme
503        // is set to full!
504        // Only do something if the contact model is of the same type
505        if (old->getName().compare("burger",Qt::CaseInsensitive)) {
506            ContactModelBurger *oldCm = (ContactModelBurger *)old;
507            knk_ = oldCm->knk_;
508            cnk_ = oldCm->cnk_;
509            knm_ = oldCm->knm_;
510            cnm_ = oldCm->cnm_;
511            ksk_ = oldCm->ksk_;
512            csk_ = oldCm->csk_;
513            ksm_ = oldCm->ksm_;
514            csm_ = oldCm->csm_;
515            fric_ = oldCm->fric_;
516            bmode_ = oldCm->bmode_;
517        }
518    }
519
520    DVect ContactModelBurger::getForce(const IContactMechanical *) const {
521        DVect ret(force_);
522        return ret;
523    }
524
525    DAVect ContactModelBurger::getMomentOn1(const IContactMechanical *c) const {
526        DVect force = getForce(c);
527        DAVect ret(0.0);
528        c->updateResultingTorqueOn1Local(force,&ret);
529        return ret;
530    }
531
532    DAVect ContactModelBurger::getMomentOn2(const IContactMechanical *c) const {
533        DVect force = getForce(c);
534        DAVect ret(0.0);
535        c->updateResultingTorqueOn2Local(force,&ret);
536        return ret;
537    }
538
539} // namespace cmodelsxd
540// EoF

Top