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 ×tep);
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 ×tep);
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
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 ×tep) {
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
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Sep 26, 2024 |
