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 BURGER_EXPORT ContactModelBurger(const ContactModelBurger &) noexcept;
22 BURGER_EXPORT const ContactModelBurger & operator=(const ContactModelBurger &);
23 BURGER_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
24 // Destructor, called when contact is deleted: free allocated memory, etc.
25 BURGER_EXPORT virtual ~ContactModelBurger();
26 // Contact model name (used as keyword for commands and FISH).
27 virtual QString getName() const { return "burger"; }
28 // The index provides a quick way to determine the type of contact model.
29 // Each type of contact model in PFC must have a unique index; this is assigned
30 // by PFC when the contact model is loaded. This index should be set to -1
31 virtual void setIndex(int i) { index_=i;}
32 virtual int getIndex() const {return index_;}
33 // Contact model version number (e.g., MyModel05_1). The version number can be
34 // accessed during the save-restore operation (within the archive method,
35 // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for
36 // future modifications to the contact model data structure.
37 virtual uint32 getMinorVersion() const;
38 // Copy the state information to a newly created contact model.
39 // Provide access to state information, for use by copy method.
40 virtual void copy(const ContactModel *c) override;
41 // Provide save-restore capability for the state information.
42 virtual void archive(ArchiveStream &);
43 // Enumerator for the properties.
44 enum PropertyKeys {
45 kwKnK=1
46 , kwCnK
47 , kwKnM
48 , kwCnM
49 , kwKsK
50 , kwCsK
51 , kwKsM
52 , kwCsM
53 , kwMode
54 , kwFric
55 , kwF
56 , kwS
57 };
58 // Contact model property names in a comma separated list. The order corresponds with
59 // the order of the PropertyKeys enumerator above. One can visualize any of these
60 // properties in PFC automatically.
61 virtual QString getProperties() const {
62 return "bur_knk"
63 ",bur_cnk"
64 ",bur_knm"
65 ",bur_cnm"
66 ",bur_ksk"
67 ",bur_csk"
68 ",bur_ksm"
69 ",bur_csm"
70 ",bur_mode"
71 ",bur_fric"
72 ",bur_force"
73 ",bur_slip";
74 }
75 // Enumerator for contact model related FISH callback events.
76 enum FishCallEvents {
77 fActivated=0
78 ,fSlipChange
79 };
80 // Contact model FISH callback event names in a comma separated list. The order corresponds with
81 // the order of the FishCallEvents enumerator above.
82 virtual QString getFishCallEvents() const {
83 return
84 "contact_activated"
85 ",slip_change";
86 }
87
88 // Return the specified contact model property.
89 virtual QVariant getProperty(uint32 i,const IContact *) const;
90 // The return value denotes whether or not the property corresponds to the global
91 // or local coordinate system (TRUE: global system, FALSE: local system). The
92 // local system is the contact-plane system (nst) defined as follows.
93 // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
94 // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
95 // and tc are unit vectors in directions of the nst axes.
96 // This is used when rendering contact model properties that are vectors.
97 virtual bool getPropertyGlobal(uint32 i) const;
98 // Set the specified contact model property, ensuring that it is of the correct type
99 // and within the correct range --- if not, then throw an exception.
100 // The return value denotes whether or not the update has affected the timestep
101 // computation (by having modified the translational or rotational tangent stiffnesses).
102 // If true is returned, then the timestep will be recomputed.
103 virtual bool setProperty(uint32 i,const QVariant &v,IContact *);
104 // The return value denotes whether or not the property is read-only
105 // (TRUE: read-only, FALSE: read-write).
106 virtual bool getPropertyReadOnly(uint32 i) const;
107
108 // Prepare for entry into ForceDispLaw. The validate function is called when:
109 // (1) the contact is created, (2) a property of the contact that returns a true via
110 // the setProperty method has been modified and (3) when a set of cycles is executed
111 // via the {cycle N} command.
112 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
113 virtual bool validate(ContactModelMechanicalState *state,const double ×tep);
114 // The forceDisplacementLaw function is called during each cycle. Given the relative
115 // motion of the two contacting pieces (via
116 // state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
117 // state->relativeAngularIncrement_ (Dtt, Dtbs, Dtbt)
118 // Ddn : relative normal-displacement increment, Ddn > 0 is opening
119 // Ddss : relative shear-displacement increment (s-axis component)
120 // Ddst : relative shear-displacement increment (t-axis component)
121 // Dtt : relative twist-rotation increment
122 // Dtbs : relative bend-rotation increment (s-axis component)
123 // Dtbt : relative bend-rotation increment (t-axis component)
124 // The relative displacement and rotation increments:
125 // Dd = Ddn*nc + Ddss*sc + Ddst*tc
126 // Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
127 // where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
128 // [see {Table 1: Contact State Variables} in PFC Model Components:
129 // Contacts and Contact Models: Contact Resolution]
130 // ) and the contact properties, this function must update the contact force and
131 // moment.
132 // The force_ is acting on piece 2, and is expressed in the local coordinate system
133 // (defined in getPropertyGlobal) such that the first component positive denotes
134 // compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
135 // in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to
136 // get the total moment.
137 // The return value indicates the contact activity status (TRUE: active, FALSE:
138 // inactive) during the next cycle.
139 // Additional information:
140 // * If state->activated() is true, then the contact has just become active (it was
141 // inactive during the previous time step).
142 // * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
143 // the forceDispLaw handle the case of {state->canFail_ == true}.
144 virtual bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep);
145 // The getEffectiveXStiffness functions return the translational and rotational
146 // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
147 // the translational tangent shear stiffness is zero (but this stiffness reduction
148 // is typically ignored when computing a stable time step). If the contact model
149 // includes a dashpot, then the translational stiffnesses must be increased (see
150 // Potyondy (2009)).
151 // [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
152 // Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
153 // December 7, 2009.]
154 virtual DVect2 getEffectiveTranslationalStiffness() const { return DVect2(knk_,ksk_); }
155
156 // Return a new instance of the contact model. This is used in the CMAT
157 // when a new contact is created.
158 virtual ContactModelBurger *clone() const override { return NEWC(ContactModelBurger()); }
159 // The getActivityDistance function is called by the contact-resolution logic when
160 // the CMAT is modified. Return value is the activity distance used by the
161 // checkActivity function.
162 virtual double getActivityDistance() const {return 0.0;}
163 // The isOKToDelete function is called by the contact-resolution logic when...
164 // Return value indicates whether or not the contact may be deleted.
165 // If TRUE, then the contact may be deleted when it is inactive.
166 // If FALSE, then the contact may not be deleted (under any condition).
167 virtual bool isOKToDelete() const { return !isBonded(); }
168 // Zero the forces and moments stored in the contact model. This function is called
169 // when the contact becomes inactive.
170 virtual void resetForcesAndMoments() { force_.fill(0.0);}
171 virtual void setForce(const DVect &v,IContact *) { force_ = v; }
172 virtual void setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }
173 virtual double getArea() const { return 0.0; }
174 // The checkActivity function is called by the contact-resolution logic when...
175 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
176 virtual bool checkActivity(const double &gap) { return gap <= 0.0; }
177
178 // Returns the sliding state (FALSE is returned if not implemented).
179 virtual bool isSliding() const { return s_; }
180 // Returns the bonding state (FALSE is returned if not implemented).
181 virtual bool isBonded() const { return false; }
182
183 virtual bool endPropertyUpdated(const QString &,const IContactMechanical *) {return false;}
184 // Both of these methods are called only for contacts with facets where the wall
185 // resolution scheme is set the full. In such cases one might wish to propagate
186 // contact state information (e.g., shear force) from one active contact to another.
187 // See the Faceted Wall section in the documentation.
188 virtual void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes());
189 virtual void setNonForcePropsFrom(IContactModel *oldCM);
190
191 /// Return the total force that the contact model holds.
192 virtual DVect getForce() const;
193
194 /// Return the total moment on 1 that the contact model holds
195 virtual DAVect getMomentOn1(const IContactMechanical *) const;
196
197 /// Return the total moment on 1 that the contact model holds
198 virtual DAVect getMomentOn2(const IContactMechanical *) const;
199
200 virtual DAVect getModelMomentOn1() const;
201 virtual DAVect getModelMomentOn2() const;
202
203 // Used to efficiently get properties from the contact model for the object pane.
204 // List of properties for the object pane, comma separated.
205 // All properties will be cast to doubles for comparison. No other comparisons
206 // are supported. This may not be the same as the entire property list.
207 // Return property name and type for plotting.
208 void objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *) const override;
209 // All properties cast to doubles - this is what can be compared.
210 void objectPropValues(std::vector<double> *,const IContact *) const override;
211
212 // Methods to get and set properties.
213 const double & knk() const {return knk_;}
214 void knk(const double &d) {knk_=d;}
215 const double & cnk() const {return cnk_;}
216 void cnk(const double &d) {cnk_=d;}
217 const double & knm() const {return knm_;}
218 void knm(const double &d) {knm_=d;}
219 const double & cnm() const {return cnm_;}
220 void cnm(const double &d) {cnm_=d;}
221 const double & ksk() const {return ksk_;}
222 void ksk(const double &d) {ksk_=d;}
223 const double & csk() const {return csk_;}
224 void csk(const double &d) {csk_=d;}
225 const double & ksm() const {return ksm_;}
226 void ksm(const double &d) {ksm_=d;}
227 const double & csm() const {return csm_;}
228 void csm(const double &d) {csm_=d;}
229
230 const double & fric() const { return fric_;}
231 void fric(const double &d) {fric_=d;}
232 uint32 bmode() const {return bmode_;}
233 void bmode(uint32 b) {bmode_= b;}
234
235 const double & fn0() const { return fn0_;}
236 void fn0(const double &d) {fn0_=d;}
237 const double & u_n0() const { return u_n0_;}
238 void u_n0(const double &d) {u_n0_=d;}
239 const double & u_nk0() const { return u_nk0_;}
240 void u_nk0(const double &d) {u_nk0_=d;}
241 const DVect & u_sk() const { return u_sk_;}
242 void u_sk(const DVect &v) {u_sk_=v;}
243
244 const double & conAn() const { return conAn_;}
245 void conAn(const double &d) {conAn_=d;}
246 const double & conB_An() const { return conB_An_;}
247 void conB_An(const double &d) {conB_An_=d;}
248 const double & conCn() const { return conCn_;}
249 void conCn(const double &d) {conCn_=d;}
250 const double & conDn() const { return conDn_;}
251 void conDn(const double &d) {conDn_=d;}
252
253 const double & conAs() const { return conAs_;}
254 void conAs(const double &d) {conAs_=d;}
255 const double & conB_As() const { return conB_As_;}
256 void conB_As(const double &d) {conB_As_=d;}
257 const double & conCs() const { return conCs_;}
258 void conCs(const double &d) {conCs_=d;}
259 const double & conDs() const { return conDs_;}
260 void conDs(const double &d) {conDs_=d;}
261
262 const double & tdel() const { return tdel_;}
263 void tdel(const double &d) {tdel_=d;}
264 const DVect & force() const { return force_;}
265 void force(const DVect &v) {force_=v;}
266 bool s() const {return s_;}
267 void s(bool b) {s_= b;}
268
269 private:
270 // Index - used internally by PFC. Should be set to -1 in the cpp file.
271 static int index_;
272
273 // Burger model properties
274 double knk_ = 0.0; // normal stiffness for Kelvin section (bur_knk)
275 double cnk_ = 0.0; // normal viscosity for Kelvin section (bur_cnk)
276 double knm_ = 0.0; // normal stiffness for Maxwell section (bur_knm)
277 double cnm_ = 0.0; // normal viscosity for Maxwell section (bur_cnm)
278 double ksk_ = 0.0; // shear stiffness for Kelvin section (bur_ksk)
279 double csk_ = 0.0; // shear viscosity for Kelvin section (bur_csk)
280 double ksm_ = 0.0; // shear stiffness for Maxwell section (bur_ksm)
281 double csm_ = 0.0; // shear viscosity for Maxwell section (bur_csm)
282 double fric_ = 0.0; // friction coefficient (bur_fric)
283 uint32 bmode_ = 0; // FDLaw option, with or without tensile force, default false (with tensile)
284 double fn0_ = 0.0; // normal contact force 1 step before
285 double u_n0_ = 0.0; // normal total overlap 1 step before
286 double u_nk0_ = 0.0; // normal overlap of Kelvin part 1step before
287 DVect u_sk_ = DVect(0.0); // shear relative displacement of Kelvin part 1step before
288 double conAn_ = 0.0; // constant A in eq.(), normal
289 double conB_An_ = 0.0; // constant B/A in eq.(), normal
290 double conCn_ = 0.0; // constant C in eq.(), normal
291 double conDn_ = 0.0; // constant D in eq.(), normal
292 double conAs_ = 0.0; // constant A in eq.(), shear
293 double conB_As_ = 0.0; // constant B/A in eq.(), shear
294 double conCs_ = 0.0; // constant C in eq.(), shear
295 double conDs_ = 0.0; // constant D in eq.(), shear
296 double tdel_ = 0.0; // current timestep
297 DVect force_ = DVect(0.0); // current total force
298 bool s_ = false; // current sliding state
299
300 // Constants
301 inline double A(const double k_k, const double c_k) { return(1.0 + k_k*tdel_/(2.0*c_k)); }
302 inline double B(const double k_k, const double c_k) { return(1.0 - k_k*tdel_/(2.0*c_k)); }
303 inline double B_A(const double k_k, const double c_k) { return(B(k_k,c_k)/A(k_k,c_k)); }
304 inline double C(const double k_k, const double c_k, const double k_m, const double c_m) {
305 return(tdel_/(2.0*c_k*A(k_k,c_k)) + 1.0/k_m + tdel_/(2.0*c_m)); }
306 inline double D(const double k_k, const double c_k, const double k_m, const double c_m) {
307 return(tdel_/(2.0*c_k*A(k_k,c_k)) - 1.0/k_m + tdel_/(2.0*c_m)); }
308
309 };
310} // namespace cmodelsxd
311// 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/ifishcalllist.h"
7#include "fish/src/parameter.h"
8
9#include "../version.txt"
10
11#include "shared/src/mathutil.h"
12
13#include "kernel/interface/iprogram.h"
14
15#ifdef BURGER_LIB
16#ifdef _WIN32
17 int __stdcall DllMain(void *,unsigned, void *) {
18 return 1;
19 }
20#endif
21
22 extern "C" EXPORT_TAG const char *getName() {
23#if DIM==3
24 return "contactmodelmechanical3dburger";
25#else
26 return "contactmodelmechanical2dburger";
27#endif
28 }
29
30 extern "C" EXPORT_TAG unsigned getMajorVersion() {
31 return MAJOR_VERSION;
32 }
33
34 extern "C" EXPORT_TAG unsigned getMinorVersion() {
35 return MINOR_VERSION;
36 }
37
38 extern "C" EXPORT_TAG void *createInstance() {
39 cmodelsxd::ContactModelBurger *m = NEWC(cmodelsxd::ContactModelBurger());
40 return (void *)m;
41 }
42#endif
43
44namespace cmodelsxd {
45 using namespace itasca;
46
47 int ContactModelBurger::index_ = -1;
48 uint32 ContactModelBurger::getMinorVersion() const { return MINOR_VERSION;}
49
50 ContactModelBurger::ContactModelBurger() {
51 }
52
53 ContactModelBurger::ContactModelBurger(const ContactModelBurger &m) noexcept {
54 this->copy(&m);
55 }
56
57 const ContactModelBurger & ContactModelBurger::operator=(const ContactModelBurger &m) {
58 this->copy(&m);
59 return *this;
60 }
61
62 void ContactModelBurger::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
63 s->addToStorage<ContactModelBurger>(*this,ww);
64 }
65
66 ContactModelBurger::~ContactModelBurger() {
67 }
68
69 void ContactModelBurger::archive(ArchiveStream &stream) {
70 // The stream allows one to archive the values of the contact model
71 // so that it can be saved and restored. The minor version can be
72 // used here to allow for incremental changes to the contact model too.
73 stream & knk_;
74 stream & cnk_;
75 stream & knm_;
76 stream & cnm_;
77 stream & ksk_;
78 stream & csk_;
79 stream & ksm_;
80 stream & csm_;
81 stream & fric_;
82 stream & bmode_;
83 stream & fn0_;
84 stream & u_n0_;
85 stream & u_nk0_;
86 stream & u_sk_;
87 stream & conAn_;
88 stream & conB_An_;
89 stream & conCn_;
90 stream & conDn_;
91 stream & conAs_;
92 stream & conB_As_;
93 stream & conCs_;
94 stream & conDs_;
95 stream & tdel_;
96 stream & force_;
97 stream & s_;
98 }
99
100 void ContactModelBurger::copy(const ContactModel *cm) {
101 // Copy all of the contact model properties. Used in the CMAT
102 // when a new contact is created.
103 ContactModelMechanical::copy(cm);
104 const ContactModelBurger *in = dynamic_cast<const ContactModelBurger*>(cm);
105 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
106
107 knk(in->knk());
108 cnk(in->cnk());
109 knm(in->knm());
110 cnm(in->cnm());
111 ksk(in->ksk());
112 csk(in->csk());
113 ksm(in->ksm());
114 csm(in->csm());
115 fric(in->fric());
116 bmode(in->bmode());
117 fn0(in->fn0());
118 u_n0(in->u_n0());
119 u_nk0(in->u_nk0());
120 u_sk(in->u_sk());
121 conAn(in->conAn());
122 conB_An(in->conB_An());
123 conCn(in->conCn());
124 conDn(in->conDn());
125 conAs(in->conAs());
126 conB_As(in->conB_As());
127 conCs(in->conCs());
128 conDs(in->conDs());
129 tdel(in->tdel());
130 force(in->force());
131 s(in->s());
132 }
133
134
135 QVariant ContactModelBurger::getProperty(uint32 i,const IContact *) const {
136 // Return the property. The IContact pointer is provided so that
137 // more complicated properties, depending on contact characteristics,
138 // can be calculated. Not used with the Burger model. The IContact pointer may be a nullptr!
139 QVariant var;
140 switch (i) {
141 case kwKnK : return knk_;
142 case kwCnK : return cnk_;
143 case kwKnM : return knm_;
144 case kwCnM : return cnm_;
145 case kwKsK : return ksk_;
146 case kwCsK : return csk_;
147 case kwKsM : return ksm_;
148 case kwCsM : return csm_;
149 case kwMode : return bmode_;
150 case kwFric : return fric_;
151 case kwF : var.setValue(force_); return var;
152 case kwS : return s_;
153 }
154 assert(0);
155 return QVariant();
156 }
157
158 bool ContactModelBurger::getPropertyGlobal(uint32 i) const {
159 // Returns whether or not a property is held in the global axis system (TRUE)
160 // or the local system (FALSE). Used by the plotting logic.
161 switch (i) {
162 case kwF:
163 return false;
164 }
165 return true;
166 }
167
168 bool ContactModelBurger::setProperty(uint32 i,const QVariant &v,IContact *) {
169 // Set a contact model property. Return value indicates that the timestep
170 // should be recalculated.
171 switch (i) {
172 case kwKnK: {
173 if (!v.canConvert<double>())
174 throw Exception("bur_knk must be a double.");
175 double val(v.toDouble());
176 if (val<0.0)
177 throw Exception("Negative bur_knk not allowed.");
178 knk_ = val;
179 return true;
180 }
181 case kwCnK: {
182 if (!v.canConvert<double>())
183 throw Exception("bur_cnk must be a double.");
184 double val(v.toDouble());
185 if (val<0.0)
186 throw Exception("Negative bur_cnk not allowed.");
187 cnk_ = val;
188 return true;
189 }
190 case kwKnM: {
191 if (!v.canConvert<double>())
192 throw Exception("bur_knm must be a double.");
193 double val(v.toDouble());
194 if (val<0.0)
195 throw Exception("Negative bur_knm not allowed.");
196 knm_ = val;
197 return true;
198 }
199 case kwCnM: {
200 if (!v.canConvert<double>())
201 throw Exception("bur_cnm must be a double.");
202 double val(v.toDouble());
203 if (val<0.0)
204 throw Exception("Negative bur_cnm not allowed.");
205 cnm_ = val;
206 return true;
207 }
208 case kwKsK: {
209 if (!v.canConvert<double>())
210 throw Exception("bur_ksk must be a double.");
211 double val(v.toDouble());
212 if (val<0.0)
213 throw Exception("Negative bur_ksk not allowed.");
214 ksk_ = val;
215 return true;
216 }
217 case kwCsK: {
218 if (!v.canConvert<double>())
219 throw Exception("bur_csk must be a double.");
220 double val(v.toDouble());
221 if (val<0.0)
222 throw Exception("Negative bur_csk not allowed.");
223 csk_ = val;
224 return true;
225 }
226 case kwKsM: {
227 if (!v.canConvert<double>())
228 throw Exception("bur_ksm must be a double.");
229 double val(v.toDouble());
230 if (val<0.0)
231 throw Exception("Negative bur_ksm not allowed.");
232 ksm_ = val;
233 return true;
234 }
235 case kwCsM: {
236 if (!v.canConvert<double>())
237 throw Exception("bur_csm must be a double.");
238 double val(v.toDouble());
239 if (val<0.0)
240 throw Exception("Negative bur_csm not allowed.");
241 csm_ = val;
242 return true;
243 }
244 case kwMode: {
245 if (!v.canConvert<uint32>())
246 throw Exception("bur_mode must be 0 (tensile) or 1 (no-tension).");
247 uint32 val(v.toUInt());
248 if (val >1)
249 throw Exception("bur_mode must be 0 (tensile) or 1 (no-tension).");
250 bmode_ = val;
251 return false;
252 }
253 case kwFric: {
254 if (!v.canConvert<double>())
255 throw Exception("bur_fric must be a double.");
256 double val(v.toDouble());
257 if (val<0.0)
258 throw Exception("Negative bur_fric not allowed.");
259 fric_ = val;
260 return false;
261 }
262 }
263 return false;
264 }
265
266 bool ContactModelBurger::getPropertyReadOnly(uint32 i) const {
267 // Returns TRUE if a property is read only or FALSE otherwise.
268 switch (i) {
269 case kwF:
270 case kwS:
271 return true;
272 default:
273 break;
274 }
275 return false;
276 }
277
278 bool ContactModelBurger::validate(ContactModelMechanicalState *state,const double &) {
279 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
280 // (1) the contact is created, (2) a property of the contact that returns a true via
281 // the setProperty method has been modified and (3) when a set of cycles is executed
282 // via the {cycle N} command.
283 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
284 return checkActivity(state->gap_);
285 }
286
287 bool ContactModelBurger::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
288 assert(state);
289
290 // Current overlap
291 double overlap = - 1.0*state->gap_;
292 // Relative translational increment
293 DVect trans = state->relativeTranslationalIncrement_;
294 // Correction factor to account for when the contact becomes newly active.
295 // We estimate the time of activity during the timestep when the contact has first
296 // become active and scale the forces accordingly.
297 double correction = 1.0;
298 // The contact was just activated from an inactive state
299 if (state->activated()) {
300 // Trigger the FISH callback if one is hooked up to the
301 // contact_activated event.
302 if (cmEvents_[fActivated] >= 0) {
303 // An FArray of QVariant is returned and these will be passed
304 // to the FISH function as an array of FISH symbols as the second
305 // argument to the FISH callback function.
306 auto c = state->getContact();
307 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
308 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
309 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
310 }
311 // Calculate the correction factor.
312 if (trans.x()) {
313 correction = -1.0*overlap / trans.x();
314 if (correction < 0)
315 correction = 1.0;
316 }
317 }
318
319 trans*=correction;
320
321 if (timestep!=tdel_) { // re-calculated constants.
322 tdel_ = timestep;
323 // need some protection for divided by zero (k_k c_k k_m c_m = zero)
324 conAn_ = A(knk_, cnk_);
325 double conBn = B(knk_, cnk_);
326 conB_An_ = conBn / conAn_;
327 conCn_ = C(knk_, cnk_, knm_, cnm_);
328 conDn_ = D(knk_, cnk_, knm_, cnm_);
329 conAs_ = A(ksk_, csk_);
330 double conBs = B(ksk_, csk_);
331 conB_As_ = conBs / conAs_;
332 conCs_ = C(ksk_, csk_, ksm_, csm_);
333 conDs_ = D(ksk_, csk_, ksm_, csm_);
334 }
335
336
337 // normal force
338 force_.rx() = 1.0/conCn_*(overlap-u_n0_+(1.0-conB_An_)*u_nk0_-conDn_*fn0_);
339 if (bmode_ && force_.x()<0.0) force_.rx() = 0.0;
340 u_nk0_ = conB_An_*u_nk0_+timestep/(2.0*cnk_*conAn_)*(force_.x()+fn0_);
341 u_n0_ = overlap;
342 fn0_ = force_.x();
343
344 // Calculate the shear force.
345 DVect sforce(0.0);
346 DVect sforce_old = force_;
347 sforce_old.rx()=0.0;
348 DVect v1 = trans;
349 DVect v2 = u_sk_ * (1.0-conB_As_);
350 DVect v3 = sforce_old * conDs_;
351 sforce = (v1+v2+v3) / conCs_ * (-1.0);
352
353 double d1 = timestep / (2.0*csk_*conAs_);
354 sforce.rx() = 0.0;
355 v1 = sforce + sforce_old;
356 u_sk_ = u_sk_*conB_As_-v1*d1;
357
358 // The canFail flag corresponds to whether or not the contact can undergo non-linear
359 // force-displacement response. If the SOLVE ELASTIC command is given then the
360 // canFail state is set to FALSE. Otherwise it is always TRUE.
361 if (state->canFail_) {
362 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
363 double crit = force_.x() * fric_;
364 crit = max(0.0,crit);
365 // This is the magnitude of the shear force.
366 double sfmag = sforce.mag();
367 // Sliding occurs when the magnitude of the shear force is greater than the
368 // critical value.
369 if (sfmag > crit) {
370 // Lower the shear force to the critical value for sliding.
371 double rat = crit / sfmag;
372 sforce *= rat;
373 // Handle the slip_change event if one has been hooked up. Sliding has commenced.
374 if (!s_ && cmEvents_[fSlipChange] >= 0) {
375 auto c = state->getContact();
376 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
377 fish::Parameter() };
378 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
379 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
380 }
381 s_ = true;
382 } else {
383 // Handle the slip_change event if one has been hooked up and
384 // the contact was previously sliding. Sliding has ceased.
385 if (s_) {
386 if (cmEvents_[fSlipChange] >= 0) {
387 auto c = state->getContact();
388 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
389 fish::Parameter((qint64)1) };
390 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
391 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
392 }
393 s_ = false;
394 }
395 }
396 }
397
398 // Set the shear components of the total force.
399 for (int i=1; i<dim; ++i)
400 force_.rdof(i) = sforce.dof(i);
401 // This is just a sanity check to ensure, in debug mode, that the force isn't wonky.
402 assert(force_ == force_);
403 return true;
404 }
405
406 void ContactModelBurger::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
407 // Only called for contacts with wall facets when the wall resolution scheme
408 // is set to full!
409 // Only do something if the contact model is of the same type
410 if (old->getContactModel()->getName().compare("burger",Qt::CaseInsensitive) == 0 && !isBonded()) {
411 ContactModelBurger *oldCm = (ContactModelBurger *)old;
412#ifdef THREED
413 // Need to rotate just the shear component from oldSystem to newSystem
414
415 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
416 DVect axis = oldSystem.e1() & newSystem.e1();
417 double c, ang, s;
418 DVect re2;
419 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
420 axis = axis.unit();
421 c = oldSystem.e1()|newSystem.e1();
422 if (c > 0)
423 c = std::min(c,1.0);
424 else
425 c = std::max(c,-1.0);
426 ang = acos(c);
427 s = sin(ang);
428 double t = 1. - c;
429 DMatrix<3,3> rm;
430 rm.get(0,0) = t*axis.x()*axis.x() + c;
431 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
432 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
433 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
434 rm.get(1,1) = t*axis.y()*axis.y() + c;
435 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
436 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
437 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
438 rm.get(2,2) = t*axis.z()*axis.z() + c;
439 re2 = rm*oldSystem.e2();
440 }
441 else
442 re2 = oldSystem.e2();
443 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
444 axis = re2 & newSystem.e2();
445 DVect2 tpf;
446 DMatrix<2,2> m;
447 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
448 axis = axis.unit();
449 c = re2|newSystem.e2();
450 if (c > 0)
451 c = std::min(c,1.0);
452 else
453 c = std::max(c,-1.0);
454 ang = acos(c);
455 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
456 ang *= -1;
457 s = sin(ang);
458 m.get(0,0) = c;
459 m.get(1,0) = s;
460 m.get(0,1) = -m.get(1,0);
461 m.get(1,1) = m.get(0,0);
462 tpf = m*DVect2(oldCm->force_.y(),oldCm->force_.z());
463 } else {
464 m.get(0,0) = 1.;
465 m.get(0,1) = 0.;
466 m.get(1,0) = 0.;
467 m.get(1,1) = 1.;
468 tpf = DVect2(oldCm->force_.y(),oldCm->force_.z());
469 }
470 DVect pforce = DVect(0,tpf.x(),tpf.y());
471#else
472 oldSystem;
473 newSystem;
474 DVect pforce = DVect(0,oldCm->force_.y());
475#endif
476 for (int i=1; i<dim; ++i)
477 force_.rdof(i) += pforce.dof(i);
478 oldCm->force_ = DVect(0.0);
479 }
480 assert(force_ == force_);
481 }
482
483 void ContactModelBurger::setNonForcePropsFrom(IContactModel *old) {
484 // Only called for contacts with wall facets when the wall resolution scheme
485 // is set to full!
486 // Only do something if the contact model is of the same type
487 if (old->getName().compare("burger",Qt::CaseInsensitive)) {
488 ContactModelBurger *oldCm = (ContactModelBurger *)old;
489 knk_ = oldCm->knk_;
490 cnk_ = oldCm->cnk_;
491 knm_ = oldCm->knm_;
492 cnm_ = oldCm->cnm_;
493 ksk_ = oldCm->ksk_;
494 csk_ = oldCm->csk_;
495 ksm_ = oldCm->ksm_;
496 csm_ = oldCm->csm_;
497 fric_ = oldCm->fric_;
498 bmode_ = oldCm->bmode_;
499 }
500 }
501
502 DVect ContactModelBurger::getForce() const {
503 DVect ret(force_);
504 return ret;
505 }
506
507 DAVect ContactModelBurger::getMomentOn1(const IContactMechanical *c) const {
508 DVect force = getForce();
509 DAVect ret(0.0);
510 c->updateResultingTorqueOn1Local(force,&ret);
511 return ret;
512 }
513
514 DAVect ContactModelBurger::getMomentOn2(const IContactMechanical *c) const {
515 DVect force = getForce();
516 DAVect ret(0.0);
517 c->updateResultingTorqueOn2Local(force,&ret);
518 return ret;
519 }
520
521 DAVect ContactModelBurger::getModelMomentOn1() const {
522 DAVect ret(0.0);
523 return ret;
524 }
525
526 DAVect ContactModelBurger::getModelMomentOn2() const {
527 DAVect ret(0.0);
528 return ret;
529 }
530
531 void ContactModelBurger::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
532 ret->clear();
533 ret->push_back({"bur_knk",ScalarInfo});
534 ret->push_back({"bur_cnk",ScalarInfo});
535 ret->push_back({"bur_knm",ScalarInfo});
536 ret->push_back({"bur_cnm",ScalarInfo});
537 ret->push_back({"bur_ksk",ScalarInfo});
538 ret->push_back({"bur_csk",ScalarInfo});
539 ret->push_back({"bur_ksm",ScalarInfo});
540 ret->push_back({"bur_csm",ScalarInfo});
541 ret->push_back({"bur_mode",ScalarInfo});
542 ret->push_back({"bur_fric",ScalarInfo});
543 ret->push_back({"bur_force",VectorInfo});
544 ret->push_back({"bur_slip",ScalarInfo});
545
546 }
547
548 void ContactModelBurger::objectPropValues(std::vector<double> *ret,const IContact *) const {
549 ret->clear();
550 ret->push_back(knk_);
551 ret->push_back(cnk_);
552 ret->push_back(knm_);
553 ret->push_back(cnm_);
554 ret->push_back(ksk_);
555 ret->push_back(csk_);
556 ret->push_back(ksm_);
557 ret->push_back(csm_);
558 ret->push_back(bmode_);
559 ret->push_back(fric_);
560 ret->push_back(force_.mag());
561 ret->push_back(s_);
562 }
563
564} // namespace cmodelsxd
565// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Oct 31, 2024 |
