Soft-Bond Model Implementation
See this page for the documentation of this contact model.
contactmodelsoftbond.h
1#pragma once
2// contactmodelsoftbond.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef SOFTBOND_LIB
7# define SOFTBOND_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define SOFTBOND_EXPORT
10#else
11# define SOFTBOND_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelSoftBond : public ContactModelMechanical {
18 public:
19 // Constructor: Set default values for contact model properties.
20 SOFTBOND_EXPORT ContactModelSoftBond();
21 // Destructor, called when contact is deleted: free allocated memory, etc.
22 SOFTBOND_EXPORT virtual ~ContactModelSoftBond();
23 // Contact model name (used as keyword for commands and FISH).
24 QString getName() const override { return "softbond"; }
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 void setIndex(int i) override { index_=i;}
29 int getIndex() const override {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 uint32 getMinorVersion() const override;
35 // Copy the state information to a newly created contact model.
36 // Provide access to state information, for use by copy method.
37 void copy(const ContactModel *c) override;
38 // Provide save-restore capability for the state information.
39 void archive(ArchiveStream &) override;
40 // Enumerator for the properties.
41 enum PropertyKeys {
42 kwKn=1
43 , kwKs
44 , kwFric
45 , kwBMul
46 , kwTMul
47 , kwSBMode
48 , kwSBF
49 , kwSBM
50 , kwSBS
51 , kwSBBS
52 , kwSBTS
53 , kwSBRMul
54 , kwSBRadius
55 , kwEmod
56 , kwKRatio
57 , kwDpNRatio
58 , kwDpSRatio
59 , kwDpMode
60 , kwDpF
61 , kwSBState
62 , kwSBTStr
63 , kwSBSStr
64 , kwSBCoh
65 , kwSBFa
66 , kwSBMCF
67 , kwSBSig
68 , kwSBTau
69 , kwSBSoft
70 , kwSBCut
71 , kwSBArea
72 , kwUserArea
73 , kwRGap
74 };
75 // Contact model property names in a comma separated list. The order corresponds with
76 // the order of the PropertyKeys enumerator above. One can visualize any of these
77 // properties in PFC automatically.
78 QString getProperties() const override {
79 return "kn"
80 ",ks"
81 ",fric"
82 ",sb_bmul"
83 ",sb_tmul"
84 ",sb_mode"
85 ",sb_force"
86 ",sb_moment"
87 ",sb_slip"
88 ",sb_slipb"
89 ",sb_slipt"
90 ",sb_rmul"
91 ",sb_radius"
92 ",emod"
93 ",kratio"
94 ",dp_nratio"
95 ",dp_sratio"
96 ",dp_mode"
97 ",dp_force"
98 ",sb_state"
99 ",sb_ten"
100 ",sb_shear"
101 ",sb_coh"
102 ",sb_fa"
103 ",sb_mcf"
104 ",sb_sigma"
105 ",sb_tau"
106 ",sb_soft"
107 ",sb_cut"
108 ",sb_area"
109 ",user_area"
110 ",rgap"
111 ;
112 }
113 // Enumerator for the energies.
114 enum EnergyKeys {
115 kwEStrain=1
116 , kwESlip
117 , kwEDashpot
118 };
119 // Contact model energy names in a comma separated list. The order corresponds with
120 // the order of the EnergyKeys enumerator above.
121 QString getEnergies() const override {
122 return "energy-strain"
123 ",energy-slip"
124 ",energy-dashpot";
125 }
126 // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
127 double getEnergy(uint32 i) const override;
128 // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1)
129 // returns wther or not the estrain energy is accumulated which is false).
130 bool getEnergyAccumulate(uint32 i) const override;
131 // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
132 void setEnergy(uint32 i,const double &d) override; // Base 1
133 // Activate the energy. This is only called if the energy tracking is enabled.
134 void activateEnergy() override { if (energies_) return; energies_ = NEW Energies();}
135 // Returns whether or not the energy tracking has been enabled for this contact.
136 bool getEnergyActivated() const override {return (energies_ != 0);}
137
138 // Enumerator for contact model related FISH callback events.
139 enum FishCallEvents {
140 fActivated=0
141 ,fSlipChange
142 ,fBondBreak
143 };
144 // Contact model FISH callback event names in a comma separated list. The order corresponds with
145 // the order of the FishCallEvents enumerator above.
146 QString getFishCallEvents() const override {
147 return
148 "contact_activated"
149 ",slip_change"
150 ",bond_break";
151 }
152
153 // Return the specified contact model property.
154 QVariant getProperty(uint32 i,const IContact *) const override;
155 // The return value denotes whether or not the property corresponds to the global
156 // or local coordinate system (TRUE: global system, FALSE: local system). The
157 // local system is the contact-plane system (nst) defined as follows.
158 // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
159 // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
160 // and tc are unit vectors in directions of the nst axes.
161 // This is used when rendering contact model properties that are vectors.
162 bool getPropertyGlobal(uint32 i) const override;
163 // Set the specified contact model property, ensuring that it is of the correct type
164 // and within the correct range --- if not, then throw an exception.
165 // The return value denotes whether or not the update has affected the timestep
166 // computation (by having modified the translational or rotational tangent stiffnesses).
167 // If true is returned, then the timestep will be recomputed.
168 bool setProperty(uint32 i,const QVariant &v,IContact *) override;
169 // The return value denotes whether or not the property is read-only
170 // (TRUE: read-only, FALSE: read-write).
171 bool getPropertyReadOnly(uint32 i) const override;
172
173 // The return value denotes whether or not the property is inheritable
174 // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
175 // the endPropertyUpdated method.
176 bool supportsInheritance(uint32 i) const override;
177 // Return whether or not inheritance is enabled for the specified property.
178 bool getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
179 // Set the inheritance flag for the specified property.
180 void setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
181
182 // Enumerator for contact model methods.
183 enum MethodKeys { kwDeformability=1, kwBond, kwUnbond, kwArea};
184 // Contact model methoid names in a comma separated list. The order corresponds with
185 // the order of the MethodKeys enumerator above.
186 QString getMethods() const override { return "deformability,bond,unbond,area";}
187 // Return a comma seprated list of the contact model method arguments (base 1).
188 QString getMethodArguments(uint32 i) const override;
189 // Set contact model method arguments (base 1).
190 // The return value denotes whether or not the update has affected the timestep
191 // computation (by having modified the translational or rotational tangent stiffnesses).
192 // If true is returned, then the timestep will be recomputed.
193 bool setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override;
194
195 // Prepare for entry into ForceDispLaw. The validate function is called when:
196 // (1) the contact is created, (2) a property of the contact that returns a true via
197 // the setProperty method has been modified and (3) when a set of cycles is executed
198 // via the {cycle N} command.
199 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
200 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
201 // The endPropertyUpdated method is called whenever a surface property (with a name
202 // that matches an inheritable contact model property name) of one of the contacting
203 // pieces is modified. This allows the contact model to update its associated
204 // properties. The return value denotes whether or not the update has affected
205 // the time step computation (by having modified the translational or rotational
206 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
207 bool endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
208 // The forceDisplacementLaw function is called during each cycle. Given the relative
209 // motion of the two contacting pieces (via
210 // state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
211 // state->relativeAngularIncrement_ (Dtt, Dtbs, Dtbt)
212 // Ddn : relative normal-displacement increment, Ddn > 0 is opening
213 // Ddss : relative shear-displacement increment (s-axis component)
214 // Ddst : relative shear-displacement increment (t-axis component)
215 // Dtt : relative twist-rotation increment
216 // Dtbs : relative bend-rotation increment (s-axis component)
217 // Dtbt : relative bend-rotation increment (t-axis component)
218 // The relative displacement and rotation increments:
219 // Dd = Ddn*nc + Ddss*sc + Ddst*tc
220 // Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
221 // where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
222 // [see {Table 1: Contact State Variables} in PFC Model Components:
223 // Contacts and Contact Models: Contact Resolution]
224 // ) and the contact properties, this function must update the contact force and
225 // moment.
226 // The force_ is acting on piece 2, and is expressed in the local coordinate system
227 // (defined in getPropertyGlobal) such that the first component positive denotes
228 // compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
229 // in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to
230 // get the total moment.
231 // The return value indicates the contact activity status (TRUE: active, FALSE:
232 // inactive) during the next cycle.
233 // Additional information:
234 // * If state->activated() is true, then the contact has just become active (it was
235 // inactive during the previous time step).
236 // * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
237 // the forceDispLaw handle the case of {state->canFail_ == true}.
238 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
239 bool thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
240 // The getEffectiveXStiffness functions return the translational and rotational
241 // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
242 // the translational tangent shear stiffness is zero (but this stiffness reduction
243 // is typically ignored when computing a stable time step). If the contact model
244 // includes a dashpot, then the translational stiffnesses must be increased (see
245 // Potyondy (2009)).
246 // [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
247 // Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
248 // December 7, 2009.]
249 DVect2 getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
250 DAVect getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
251
252 // Return a new instance of the contact model. This is used in the CMAT
253 // when a new contact is created.
254 ContactModelSoftBond *clone() const override { return NEW ContactModelSoftBond(); }
255 // The getActivityDistance function is called by the contact-resolution logic when
256 // the CMAT is modified. Return value is the activity distance used by the
257 // checkActivity function.
258 double getActivityDistance() const override {return rgap_;}
259 // The isOKToDelete function is called by the contact-resolution logic when...
260 // Return value indicates whether or not the contact may be deleted.
261 // If TRUE, then the contact may be deleted when it is inactive.
262 // If FALSE, then the contact may not be deleted (under any condition).
263 bool isOKToDelete() const override { return !isBonded(); }
264 // Zero the forces and moments stored in the contact model. This function is called
265 // when the contact becomes inactive.
266 void resetForcesAndMoments() override {
267 sb_F(DVect(0.0));
268 dp_F(DVect(0.0));
269 sb_M(DAVect(0.0));
270 if (energies_) {
271 energies_->estrain_ = 0.0;
272 }
273 }
274 void setForce(const DVect &v,IContact *c) override;
275 void setArea(const double &d) override { userArea_ = d; }
276 double getArea() const override { return userArea_; }
277
278 // The checkActivity function is called by the contact-resolution logic when...
279 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
280 bool checkActivity(const double &gap) override { return gap <= rgap_ || isBonded();}
281
282 // Returns the sliding state (FALSE is returned if not implemented).
283 bool isSliding() const override { return sb_S_; }
284 // Returns the bonding state (FALSE is returned if not implemented).
285 bool isBonded() const override { return bProps_ ? (bProps_->sb_state_ >= 3) : false; }
286 void unbond() override { if (bProps_) bProps_->sb_state_ = 0; }
287
288 // Both of these methods are called only for contacts with facets where the wall
289 // resolution scheme is set the full. In such cases one might wish to propagate
290 // contact state information (e.g., shear force) from one active contact to another.
291 // See the Faceted Wall section in the documentation.
292 void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
293 void setNonForcePropsFrom(IContactModel *oldCM) override;
294 /// Return the total force that the contact model holds.
295 DVect getForce(const IContactMechanical *) const override;
296 /// Return the total moment on 1 that the contact model holds
297 DAVect getMomentOn1(const IContactMechanical *) const override;
298 /// Return the total moment on 1 that the contact model holds
299 DAVect getMomentOn2(const IContactMechanical *) const override;
300
301 // Methods to get and set properties.
302 const double & kn() const {return kn_;}
303 void kn(const double &d) {kn_=d;}
304 const double & ks() const {return ks_;}
305 void ks(const double &d) {ks_=d;}
306 const double & fric() const {return fric_;}
307 void fric(const double &d) {fric_=d;}
308 const double & sb_bmul() const { return sb_bmul_; }
309 void sb_bmul(const double &d) { sb_bmul_ = d; }
310 const double & sb_tmul() const { return sb_tmul_; }
311 void sb_tmul(const double &d) { sb_tmul_ = d; }
312 const DVect & sb_F() const {return sb_F_;}
313 void sb_F(const DVect &f) { sb_F_=f;}
314 const DAVect & sb_M() const { return sb_M_; }
315 void sb_M(const DAVect &f) { sb_M_ = f; }
316 bool sb_S() const {return sb_S_;}
317 void sb_S(bool b) { sb_S_=b;}
318 bool sb_BS() const { return sb_BS_; }
319 void sb_BS(bool b) { sb_BS_ = b; }
320 bool sb_TS() const { return sb_TS_; }
321 void sb_TS(bool b) { sb_TS_ = b; }
322 const double & sb_rmul() const { return sb_rmul_; }
323 void sb_rmul(const double &d) { sb_rmul_ = d; }
324 uint32 sb_mode() const {return sb_mode_;}
325 void sb_mode(uint32 i) { sb_mode_=i;}
326
327 bool hasBond() const { return bProps_ ? true : false; }
328 int sb_state() const { return (hasBond() ? bProps_->sb_state_ : 0); }
329 void sb_state(int i) { if (!hasBond()) return; bProps_->sb_state_ = i; }
330 double sb_Ten() const { return (hasBond() ? (bProps_->sb_ten_) : 0.0); }
331 void sb_Ten(const double &d) { if (!hasBond()) return; bProps_->sb_ten_ = d; }
332 double sb_Coh() const { return (hasBond() ? (bProps_->sb_coh_) : 0.0); }
333 void sb_Coh(const double &d) { if (!hasBond()) return; bProps_->sb_coh_ = d; }
334 double sb_FA() const { return (hasBond() ? (bProps_->sb_fa_) : 0.0); }
335 void sb_FA(const double &d) { if (!hasBond()) return; bProps_->sb_fa_ = d; }
336 double sb_MCF() const {return (hasBond() ? (bProps_->sb_mcf_) : 0.0);}
337 void sb_MCF(const double &d) { if(!hasBond()) return; bProps_->sb_mcf_=d;}
338 double sb_soft() const { return (hasBond() ? (bProps_->sb_soft_) : 0.0); }
339 void sb_soft(const double &d) { if (!hasBond()) return; bProps_->sb_soft_ = d; }
340 double sb_cut() const { return (hasBond() ? (bProps_->sb_cut_) : 0.0); }
341 void sb_cut(const double &d) { if (!hasBond()) return; bProps_->sb_cut_ = d; }
342 double sb_maxTen() const { return (hasBond() ? (bProps_->sb_maxTen_) : 0.0); }
343 void sb_maxTen(const double &d) { if (!hasBond()) return; bProps_->sb_maxTen_ = d; }
344 double sb_delu() const { return (hasBond() ? (bProps_->sb_delu_) : 0.0); }
345 void sb_delu(const double &d) { if (!hasBond()) return; bProps_->sb_delu_ = d; }
346 Quat sb_delo() const { return (hasBond() ? (bProps_->sb_delo_) : Quat::identity()); }
347 void sb_delo(const Quat &d) { if (!hasBond()) return; bProps_->sb_delo_ = d; }
348 double sb_maxu() const { return (hasBond() ? (bProps_->sb_maxu_) : 0.0); }
349 void sb_maxu(const double &d) { if (!hasBond()) return; bProps_->sb_maxu_ = d; }
350 double sb_critu() const { return (hasBond() ? (bProps_->sb_critu_) : 0.0); }
351 void sb_critu(const double &d) { if (!hasBond()) return; bProps_->sb_critu_ = d; }
352
353 bool hasDamping() const {return dpProps_ ? true : false;}
354 double dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
355 void dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
356 double dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
357 void dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
358 int dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
359 void dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
360 DVect dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
361 void dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
362
363 bool hasEnergies() const {return energies_ ? true:false;}
364 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
365 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
366 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
367 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
368 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
369 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
370
371 uint32 inheritanceField() const {return inheritanceField_;}
372 void inheritanceField(uint32 i) {inheritanceField_ = i;}
373
374 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
375 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
376 const DAVect & effectiveRotationalStiffness() const {return effectiveRotationalStiffness_;}
377 void effectiveRotationalStiffness(const DAVect &v ) {effectiveRotationalStiffness_=v;}
378
379 private:
380 // Index - used internally by PFC. Should be set to -1 in the cpp file.
381 static int index_;
382
383 bool FDLawBonded(ContactModelMechanicalState *state, const double ×tep);
384 bool FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep);
385
386 // Structure to compute stiffness
387 struct StiffData {
388 DVect2 trans_ = DVect2(0.0);
389 DAVect ang_ = DAVect(0.0);
390 double reff_ = 0.0;
391 };
392
393 // Structure to store the energies.
394 struct Energies {
395 double estrain_ = 0.0; // elastic energy
396 double eslip_ = 0.0; // work dissipated by friction
397 double edashpot_ = 0.0; // work dissipated by dashpots
398 };
399
400 // Structure to store dashpot quantities.
401 struct dpProps {
402 double dp_nratio_ = 0.0; // normal viscous critical damping ratio
403 double dp_sratio_ = 0.0; // shear viscous critical damping ratio
404 int dp_mode_ = 0; // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
405 DVect dp_F_ = DVect(0.0); // Force in the dashpots
406 };
407
408 // Structure to store bond-related quantities.
409 struct bProps {
410 int sb_state_ = 0; // bond mode - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B), 4 (B-Softening), 5 (B-Compression from Softening)
411 double sb_ten_ = 0.0; // normal strength
412 double sb_coh_ = 0.0; // cohesion
413 double sb_fa_ = 0.0; // friction angle
414 double sb_mcf_ = 1.0; // moment contribution factor
415 double sb_soft_ = 0.0; // softening factor
416 double sb_cut_ = 1.0; // critical bond length
417 double sb_maxTen_ = 0.0; // tensile strength one needs to reach for softening
418 double sb_delu_ = 0.0; // incremental elongation in softening
419 Quat sb_delo_ = Quat::identity(); // incremental orientation in softening
420 double sb_maxu_ = 0.0; // max elongation for softening
421 double sb_critu_ = 0.0; // critical elongation for softening
422 };
423
424
425 bool updateKn(const IContactMechanical *con);
426 bool updateKs(const IContactMechanical *con);
427 bool updateFric(const IContactMechanical *con);
428
429 StiffData computeStiffData(ContactModelMechanicalState *state) const;
430 DVect3 computeGeomData(const IContactMechanical *c) const;
431 DVect2 SMax(const IContactMechanical *con) const; // Maximum stress (tensile,shear) at bond periphery
432 double shearStrength(const double &pbArea) const; // Bond shear strength
433 double strainEnergy(double kn, double ks, double kb, double kt) const;
434
435 void updateStiffness(ContactModelMechanicalState *state);
436
437 // Contact model inheritance fields.
438 uint32 inheritanceField_;
439
440 // Effective translational stiffness.
441 DVect2 effectiveTranslationalStiffness_;
442 DAVect effectiveRotationalStiffness_; // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
443
444 // linear model properties
445 double kn_; // Normal stiffness
446 double ks_; // Shear stiffness
447 double fric_; // Coulomb friction coefficient
448 double sb_bmul_; // Bending friction multiplier
449 double sb_tmul_; // Twisting friction multiplier
450 uint32 sb_mode_; // specifies absolute (0) or incremental (1) behavior for the the normal force
451 DVect sb_F_; // Force carried in the model
452 DAVect sb_M_; // moment (bending + twisting in 3D)
453 bool sb_S_; // The current slip state
454 bool sb_BS_; // The bending slip state
455 bool sb_TS_; // The twisting slip state
456 double sb_rmul_; // Radius multiplier
457 double userArea_; // Area as specified by the user
458 double rgap_; // Reference gap
459
460 dpProps * dpProps_; // The viscous properties
461 bProps * bProps_; // The bond properties
462
463 Energies * energies_; // The energies
464
465 };
466} // namespace cmodelsxd
467// EoF
contactmodelsoftbond.cpp
1// contactmodelsoftbond.cpp
2#include "contactmodelsoftbond.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
10#include "utility/src/tptr.h"
11#include "shared/src/mathutil.h"
12
13#include "kernel/interface/iprogram.h"
14#include "module/interface/icontactthermal.h"
15#include "module/interface/icontactfluid.h"
16#include "contactmodel/src/contactmodelthermal.h"
17#include "contactmodel/src/contactmodelfluid.h"
18#include "../version.txt"
19#include "fish/src/parameter.h"
20
21#ifdef SOFTBOND_LIB
22#ifdef _WIN32
23 int __stdcall DllMain(void *,unsigned, void *) {
24 return 1;
25 }
26#endif
27 extern "C" EXPORT_TAG const char *getName() {
28#if DIM==3
29 return "contactmodelmechanical3dsoftbond";
30#else
31 return "contactmodelmechanical2dsoftbond";
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::ContactModelSoftBond *m = NEW cmodelsxd::ContactModelSoftBond();
45 return (void *)m;
46 }
47#endif
48
49namespace cmodelsxd {
50 static const uint32 KnMask = 0x00000002; // Base 1!
51 static const uint32 KsMask = 0x00000004;
52 static const uint32 FricMask = 0x00000008;
53
54 using namespace itasca;
55
56 int ContactModelSoftBond::index_ = -1;
57 uint32 ContactModelSoftBond::getMinorVersion() const { return MINOR_VERSION; }
58
59 ContactModelSoftBond::ContactModelSoftBond() : inheritanceField_(KnMask|KsMask|FricMask)
60 , effectiveTranslationalStiffness_(DVect2(0.0))
61 , effectiveRotationalStiffness_(DAVect(0.0))
62 , kn_(0.0)
63 , ks_(0.0)
64 , fric_(0.0)
65 , sb_bmul_(1.0)
66 , sb_tmul_(1.0)
67 , sb_mode_(0)
68 , sb_F_(DVect(0.0))
69 , sb_M_(DAVect(0.0))
70 , sb_S_(false)
71 , sb_BS_(false)
72 , sb_TS_(false)
73 , sb_rmul_(1.0)
74 , userArea_(0.0)
75 , rgap_(0.0)
76 , dpProps_(nullptr)
77 , bProps_(nullptr)
78 , energies_(nullptr) {
79 }
80
81 ContactModelSoftBond::~ContactModelSoftBond() {
82 // Make sure to clean up after yourself!
83 if (dpProps_)
84 delete dpProps_;
85 if (bProps_)
86 delete bProps_;
87 if (energies_)
88 delete energies_;
89 }
90
91 void ContactModelSoftBond::archive(ArchiveStream &stream) {
92 // The stream allows one to archive the values of the contact model
93 // so that it can be saved and restored. The minor version can be
94 // used here to allow for incremental changes to the contact model too.
95 stream & kn_;
96 stream & ks_;
97 stream & fric_;
98 stream & sb_mode_;
99 stream & sb_F_;
100 stream & sb_M_;
101 stream & sb_S_;
102 stream & sb_BS_;
103 stream & sb_TS_;
104 stream & sb_rmul_;
105
106 if (stream.getArchiveState()==ArchiveStream::Save) {
107 bool b = false;
108 if (dpProps_) {
109 b = true;
110 stream & b;
111 stream & dpProps_->dp_nratio_;
112 stream & dpProps_->dp_sratio_;
113 stream & dpProps_->dp_mode_;
114 stream & dpProps_->dp_F_;
115 }
116 else
117 stream & b;
118
119 b = false;
120 if (energies_) {
121 b = true;
122 stream & b;
123 stream & energies_->estrain_;
124 stream & energies_->eslip_;
125 stream & energies_->edashpot_;
126 }
127 else
128 stream & b;
129
130 b = false;
131 if (bProps_) {
132 b = true;
133 stream & b;
134 stream & bProps_->sb_state_;
135 stream & bProps_->sb_ten_;
136 stream & bProps_->sb_coh_;
137 stream & bProps_->sb_fa_;
138 stream & bProps_->sb_mcf_;
139 stream & bProps_->sb_soft_;
140 stream & bProps_->sb_cut_;
141 stream & bProps_->sb_maxTen_;
142 stream & bProps_->sb_delu_;
143 stream & bProps_->sb_delo_;
144 stream & bProps_->sb_maxu_;
145 stream & bProps_->sb_critu_;
146 }
147 else
148 stream & b;
149
150 } else {
151 bool b(false);
152 stream & b;
153 if (b) {
154 if (!dpProps_)
155 dpProps_ = NEW dpProps();
156 stream & dpProps_->dp_nratio_;
157 stream & dpProps_->dp_sratio_;
158 stream & dpProps_->dp_mode_;
159 stream & dpProps_->dp_F_;
160 }
161 stream & b;
162 if (b) {
163 if (!energies_)
164 energies_ = NEW Energies();
165 stream & energies_->estrain_;
166 stream & energies_->eslip_;
167 stream & energies_->edashpot_;
168 }
169 stream & b;
170 if (b) {
171 if (!bProps_)
172 bProps_ = NEW bProps();
173 stream & bProps_->sb_state_;
174 stream & bProps_->sb_ten_;
175 stream & bProps_->sb_coh_;
176 stream & bProps_->sb_fa_;
177 stream & bProps_->sb_mcf_;
178 stream & bProps_->sb_soft_;
179 stream & bProps_->sb_cut_;
180 stream & bProps_->sb_maxTen_;
181 stream & bProps_->sb_delu_;
182 stream & bProps_->sb_delo_;
183 stream & bProps_->sb_maxu_;
184 stream & bProps_->sb_critu_;
185 }
186
187 }
188
189 stream & inheritanceField_;
190 stream & effectiveTranslationalStiffness_;
191 stream & effectiveRotationalStiffness_;
192
193 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 1) {
194 stream & sb_bmul_;
195 stream & sb_tmul_;
196 }
197
198 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 2)
199 stream & userArea_;
200
201 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 3)
202 stream & rgap_;
203
204 }
205
206 void ContactModelSoftBond::copy(const ContactModel *cm) {
207 // Copy all of the contact model properties. Used in the CMAT
208 // when a new contact is created.
209 ContactModelMechanical::copy(cm);
210 const ContactModelSoftBond *in = dynamic_cast<const ContactModelSoftBond*>(cm);
211 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
212 kn(in->kn());
213 ks(in->ks());
214 fric(in->fric());
215 sb_bmul(in->sb_bmul());
216 sb_tmul(in->sb_tmul());
217 sb_mode(in->sb_mode());
218 sb_F(in->sb_F());
219 sb_S(in->sb_S());
220 sb_BS(in->sb_BS());
221 sb_TS(in->sb_TS());
222 sb_rmul(in->sb_rmul());
223 sb_M(in->sb_M());
224
225 if (in->hasDamping()) {
226 if (!dpProps_)
227 dpProps_ = NEW dpProps();
228 dp_nratio(in->dp_nratio());
229 dp_sratio(in->dp_sratio());
230 dp_mode(in->dp_mode());
231 dp_F(in->dp_F());
232 }
233 if (in->hasEnergies()) {
234 if (!energies_)
235 energies_ = NEW Energies();
236 estrain(in->estrain());
237 eslip(in->eslip());
238 edashpot(in->edashpot());
239 }
240 if (in->hasBond()) {
241 if (!bProps_)
242 bProps_ = NEW bProps();
243 sb_state(in->sb_state());
244 sb_Ten(in->sb_Ten());
245 sb_Coh(in->sb_Coh());
246 sb_FA(in->sb_FA());
247 sb_MCF(in->sb_MCF());
248 sb_soft(in->sb_soft());
249 sb_cut(in->sb_cut());
250 sb_maxTen(in->sb_maxTen());
251 sb_delu(in->sb_delu());
252 sb_delo(in->sb_delo());
253 sb_maxu(in->sb_maxu());
254 sb_critu(in->sb_critu());
255 }
256 userArea_ = in->userArea_;
257 rgap_ = in->rgap_;
258 inheritanceField(in->inheritanceField());
259 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
260 effectiveRotationalStiffness(in->effectiveRotationalStiffness());
261 }
262
263
264 QVariant ContactModelSoftBond::getProperty(uint32 i,const IContact *con) const {
265 // Return the property. The IContact pointer is provided so that
266 // more complicated properties, depending on contact characteristics,
267 // can be calcualted.
268 QVariant var;
269 switch (i) {
270 case kwKn: return kn_;
271 case kwKs: return ks_;
272 case kwFric: return fric_;
273 case kwBMul: return sb_bmul_;
274 case kwTMul: return sb_tmul_;
275 case kwSBMode: return sb_mode_;
276 case kwSBF: var.setValue(sb_F_); return var;
277 case kwSBM: var.setValue(sb_M_); return var;
278 case kwSBS: return sb_S_;
279 case kwSBBS: return sb_BS_;
280 case kwSBTS: return sb_TS_;
281 case kwSBRMul: return sb_rmul_;
282 case kwSBRadius: {
283 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
284 if (!c) return 0.0;
285 double Cmax1 = c->getEnd1Curvature().y();
286 double Cmax2 = c->getEnd2Curvature().y();
287 if (!userArea_)
288 return sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
289 else {
290#ifdef THREED
291 double rad = std::sqrt(userArea_ / dPi);
292#else
293 double rad = userArea_ / 2.0;
294#endif
295 return rad;
296 }
297
298 }
299 case kwEmod: {
300 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
301 if (!c) return 0.0;
302 double rsum(0.0);
303 if (c->getEnd1Curvature().y())
304 rsum += 1.0/c->getEnd1Curvature().y();
305 if (c->getEnd2Curvature().y())
306 rsum += 1.0/c->getEnd2Curvature().y();
307 if (userArea_)
308#ifdef THREED
309 rsum = 2.0 * std::sqrt(userArea_ / dPi);
310#else
311 rsum = userArea_;
312#endif
313 return kn_ * rsum;
314 }
315 case kwKRatio: return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
316 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
317 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
318 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
319 case kwDpF: {
320 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
321 return var;
322 }
323 case kwSBState: return bProps_ ? bProps_->sb_state_ : 0;
324 case kwSBTStr: return bProps_ ? bProps_->sb_ten_ : 0.0;
325 case kwSBSStr: {
326 if (!bProps_) return 0.0;
327 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
328 double area = computeGeomData(c).x();
329 return shearStrength(area);
330 }
331 case kwSBCoh: return bProps_ ? bProps_->sb_coh_ : 0;
332 case kwSBFa: return bProps_ ? bProps_->sb_fa_ : 0;
333 case kwSBMCF: return bProps_ ? bProps_->sb_mcf_ : 0;
334 case kwSBSig: {
335 if (!bProps_ || bProps_->sb_state_ < 3) return 0.0;
336 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
337 return SMax(c).x();
338 }
339 case kwSBTau: {
340 if (!bProps_ || bProps_->sb_state_ < 3) return 0.0;
341 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
342 return SMax(c).y();
343 }
344 case kwSBSoft:
345 if (!bProps_) return 0.0;
346 return bProps_->sb_soft_;
347 case kwSBCut:
348 if (!bProps_) return 0.0;
349 return bProps_->sb_cut_;
350 case kwSBArea: {
351 if (userArea_) return userArea_;
352 //if (!bProps_) return 0.0;
353 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
354 if (!c)
355 return 0.0;
356 return computeGeomData(c).x();
357 }
358 case kwUserArea:
359 return userArea_;
360 case kwRGap:
361 return rgap_;
362 }
363 assert(0);
364 return QVariant();
365 }
366
367 bool ContactModelSoftBond::getPropertyGlobal(uint32 i) const {
368 // Returns whether or not a property is held in the global axis system (TRUE)
369 // or the local system (FALSE). Used by the plotting logic.
370 switch (i) {
371 case kwSBF:
372 case kwSBM:
373 case kwDpF:
374 return false;
375 }
376 return true;
377 }
378
379 bool ContactModelSoftBond::setProperty(uint32 i,const QVariant &v,IContact *) {
380 // Set a contact model property. Return value indicates that the timestep
381 // should be recalculated.
382 dpProps dp;
383 switch (i) {
384 case kwKn: {
385 if (!v.canConvert<double>())
386 throw Exception("kn must be a double.");
387 double val(v.toDouble());
388 if (val<0.0)
389 throw Exception("Negative kn not allowed.");
390 kn_ = val;
391 return true;
392 }
393 case kwKs: {
394 if (!v.canConvert<double>())
395 throw Exception("ks must be a double.");
396 double val(v.toDouble());
397 if (val<0.0)
398 throw Exception("Negative ks not allowed.");
399 ks_ = val;
400 return true;
401 }
402 case kwFric: {
403 if (!v.canConvert<double>())
404 throw Exception("fric must be a double.");
405 double val(v.toDouble());
406 if (val<0.0)
407 throw Exception("Negative fric not allowed.");
408 fric_ = val;
409 return false;
410 }
411 case kwBMul: {
412 if (!v.canConvert<double>())
413 throw Exception("sb_bmul must be a double.");
414 double val(v.toDouble());
415 if (val<0.0)
416 throw Exception("Negative sb_bmul not allowed.");
417 sb_bmul_ = val;
418 return false;
419 }
420 case kwTMul: {
421 if (!v.canConvert<double>())
422 throw Exception("sb_tmul must be a double.");
423 double val(v.toDouble());
424 if (val<0.0)
425 throw Exception("Negative st_bmul not allowed.");
426 sb_tmul_ = val;
427 return false;
428 }
429 case kwSBMode: {
430 if (!v.canConvert<uint32>())
431 throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
432 double val(v.toUInt());
433 if (val>1)
434 throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
435 sb_mode_ = val;
436 return false;
437 }
438 case kwSBRMul: {
439 if (!v.canConvert<double>())
440 throw Exception("rmul must be a double.");
441 double val(v.toDouble());
442 if (val<0.0)
443 throw Exception("Negative rmul not allowed.");
444 sb_rmul_ = val;
445 return false;
446 }
447 case kwSBF: {
448 if (!v.canConvert<DVect>())
449 throw Exception("sb_force must be a vector.");
450 DVect val(v.value<DVect>());
451 sb_F_ = val;
452 return false;
453 }
454 case kwSBM: {
455 DAVect val(0.0);
456#ifdef TWOD
457 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
458 throw Exception("res_moment must be an angular vector.");
459 if (v.canConvert<DAVect>())
460 val = DAVect(v.value<DAVect>());
461 else
462 val = DAVect(v.toDouble());
463#else
464 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
465 throw Exception("res_moment must be an angular vector.");
466 if (v.canConvert<DAVect>())
467 val = DAVect(v.value<DAVect>());
468 else
469 val = DAVect(v.value<DVect>());
470#endif
471 sb_M_ = val;
472 return false;
473 }
474 case kwDpNRatio: {
475 if (!v.canConvert<double>())
476 throw Exception("dp_nratio must be a double.");
477 double val(v.toDouble());
478 if (val<0.0)
479 throw Exception("Negative dp_nratio not allowed.");
480 if (val == 0.0 && !dpProps_)
481 return false;
482 if (!dpProps_)
483 dpProps_ = NEW dpProps();
484 dpProps_->dp_nratio_ = val;
485 return true;
486 }
487 case kwDpSRatio: {
488 if (!v.canConvert<double>())
489 throw Exception("dp_sratio must be a double.");
490 double val(v.toDouble());
491 if (val<0.0)
492 throw Exception("Negative dp_sratio not allowed.");
493 if (val == 0.0 && !dpProps_)
494 return false;
495 if (!dpProps_)
496 dpProps_ = NEW dpProps();
497 dpProps_->dp_sratio_ = val;
498 return true;
499 }
500 case kwDpMode: {
501 if (!v.canConvert<int>())
502 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
503 int val(v.toInt());
504 if (val == 0 && !dpProps_)
505 return false;
506 if (val < 0 || val > 3)
507 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
508 if (!dpProps_)
509 dpProps_ = NEW dpProps();
510 dpProps_->dp_mode_ = val;
511 return false;
512 }
513 case kwDpF: {
514 if (!v.canConvert<DVect>())
515 throw Exception("dp_force must be a vector.");
516 DVect val(v.value<DVect>());
517 if (val.fsum() == 0.0 && !dpProps_)
518 return false;
519 if (!dpProps_)
520 dpProps_ = NEW dpProps();
521 dpProps_->dp_F_ = val;
522 return false;
523 }
524 case kwSBTStr: {
525 if (!v.canConvert<double>())
526 throw Exception("sb_ten must be a double.");
527 double val(v.toDouble());
528 if (val < 0.0)
529 throw Exception("Negative sb_ten not allowed.");
530 if (val == 0.0 && !bProps_)
531 return false;
532 if (!bProps_)
533 bProps_ = NEW bProps();
534 bProps_->sb_ten_ = val;
535 return false;
536 }
537 case kwSBCoh: {
538 if (!v.canConvert<double>())
539 throw Exception("sb_coh must be a double.");
540 double val(v.toDouble());
541 if (val<0.0)
542 throw Exception("Negative pb_coh not allowed.");
543 if (val == 0.0 && !bProps_)
544 return false;
545 if (!bProps_)
546 bProps_ = NEW bProps();
547 bProps_->sb_coh_ = val;
548 return false;
549 }
550 case kwSBFa: {
551 if (!v.canConvert<double>())
552 throw Exception("sb_fa must be a double.");
553 double val(v.toDouble());
554 if (val<0.0)
555 throw Exception("Negative sb_fa not allowed.");
556 if (val >= 90.0)
557 throw Exception("sb_fa must be lower than 90.0 degrees.");
558 if (val == 0.0 && !bProps_)
559 return false;
560 if (!bProps_)
561 bProps_ = NEW bProps();
562 bProps_->sb_fa_ = val;
563 return false;
564 }
565 case kwSBMCF: {
566 if (!v.canConvert<double>())
567 throw Exception("sb_mcf must be a double.");
568 double val(v.toDouble());
569 if (val<0.0)
570 throw Exception("Negative sb_mcf not allowed.");
571 if (val > 1.0)
572 throw Exception("sb_mcf must be lower or equal to 1.0.");
573 if (val == 1.0 && !bProps_)
574 return false;
575 if (!bProps_)
576 bProps_ = NEW bProps();
577 bProps_->sb_mcf_ = val;
578 return false;
579 }
580 case kwSBSoft: {
581 if (!v.canConvert<double>())
582 throw Exception("sb_soft must be a double.");
583 double val(v.toDouble());
584 if (val < 0.0)
585 throw Exception("Negative pb_soft not allowed.");
586 if (!bProps_)
587 bProps_ = NEW bProps();
588 bProps_->sb_soft_ = val;
589 return false;
590 }
591 case kwSBCut: {
592 if (!v.canConvert<double>())
593 throw Exception("sb_cut must be a double.");
594 double val(v.toDouble());
595 if (val < 0.0)
596 throw Exception("Negative sb_cut not allowed.");
597 if (!bProps_)
598 bProps_ = NEW bProps();
599 bProps_->sb_cut_ = val;
600 return false;
601 }
602 case kwSBArea:
603 case kwUserArea: {
604 if (!v.canConvert<double>())
605 throw Exception("area must be a double.");
606 double val(v.toDouble());
607 if (val < 0.0)
608 throw Exception("Negative area not allowed.");
609 userArea_ = val;
610 return true;
611 }
612 case kwRGap: {
613 if (!v.canConvert<double>())
614 throw Exception("Reference gap must be a double.");
615 double val(v.toDouble());
616 rgap_ = val;
617 return false;
618 }
619 }
620 return false;
621 }
622
623 bool ContactModelSoftBond::getPropertyReadOnly(uint32 i) const {
624 // Returns TRUE if a property is read only or FALSE otherwise.
625 switch (i) {
626 case kwDpF:
627 case kwSBS:
628 case kwSBBS:
629 case kwSBTS:
630 case kwEmod:
631 case kwKRatio:
632 case kwSBState:
633 case kwSBRadius:
634 case kwSBSStr:
635 case kwSBSig:
636 case kwSBTau:
637 return true;
638 default:
639 break;
640 }
641 return false;
642 }
643
644 bool ContactModelSoftBond::supportsInheritance(uint32 i) const {
645 // Returns TRUE if a property supports inheritance or FALSE otherwise.
646 switch (i) {
647 case kwKn:
648 case kwKs:
649 case kwFric:
650 return true;
651 default:
652 break;
653 }
654 return false;
655 }
656
657 QString ContactModelSoftBond::getMethodArguments(uint32 i) const {
658 // Return a list of contact model method argument names.
659 switch (i) {
660 case kwDeformability:
661 return "emod,kratio";
662 case kwBond:
663 return "gap,soft,cut";
664 case kwUnbond:
665 return "gap";
666 case kwArea:
667 return QString();
668 }
669 assert(0);
670 return QString();
671 }
672
673 bool ContactModelSoftBond::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
674 // Apply the specified method.
675 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
676 switch (i) {
677 case kwDeformability: {
678 double emod;
679 double krat;
680 if (vl.at(0).isNull())
681 throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
682 emod = vl.at(0).toDouble();
683 if (emod<0.0)
684 throw Exception("Negative emod not allowed in contact model %1.",getName());
685 if (vl.at(1).isNull())
686 throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
687 krat = vl.at(1).toDouble();
688 if (krat<0.0)
689 throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
690 double rsum(0.0);
691 if (c->getEnd1Curvature().y())
692 rsum += 1.0 / c->getEnd1Curvature().y();
693 if (c->getEnd2Curvature().y())
694 rsum += 1.0 / c->getEnd2Curvature().y();
695 if (userArea_)
696#ifdef THREED
697 rsum = 2.0 * std::sqrt(userArea_ / dPi);
698#else
699 rsum = userArea_;
700#endif
701 kn_ = emod / rsum;
702 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
703 setInheritance(1,false);
704 setInheritance(2,false);
705 return true;
706 }
707 case kwBond: {
708 if (bProps_ && bProps_->sb_state_ >= 3) return false;
709 double mingap = -1.0 * limits<double>::max();
710 double maxgap = 0;
711 if (vl.at(0).canConvert<double>())
712 maxgap = vl.at(0).toDouble();
713 else if (vl.at(0).canConvert<DVect2>()) {
714 DVect2 value = vl.at(0).value<DVect2>();
715 mingap = value.minComp();
716 maxgap = value.maxComp();
717 }
718 else if (!vl.at(0).isNull())
719 throw Exception("gap value %1 not recognized in method bond in contact model %2.", vl.at(1), getName());
720 double soft = bProps_ ? bProps_->sb_soft_ : 0.0;
721 if (!vl.at(1).isNull()) {
722 soft = vl.at(1).toDouble();
723 if (soft < 0.0)
724 throw Exception("Negative soft not allowed in contact model %1.", getName());
725 }
726 double cut = bProps_ ? bProps_->sb_cut_ : 1.0;
727 if (!vl.at(2).isNull()) {
728 if (vl.at(2).canConvert<double>())
729 cut = vl.at(2).toDouble();
730 if (cut < 0.0)
731 throw Exception("cut value %1 is negative, or not recognized in method bond in contact model %2.", vl.at(2), getName());
732 if (cut > 1.0)
733 throw Exception("cut value %1 must be in range [0,1] in method bond in contact model %2.", vl.at(2), getName());
734 }
735 double gap = c->getGap();
736 if (gap >= mingap && gap <= maxgap) {
737 if (!bProps_)
738 bProps_ = NEW bProps();
739 bProps_->sb_state_ = 3;
740 bProps_->sb_soft_ = soft;
741 // Update the critical distance
742 if (cut != -1)
743 bProps_->sb_cut_ = cut;
744 // seet to incremental normal force calculation
745 sb_mode_ = 1;
746 return true;
747 }
748 return false;
749 }
750 case kwUnbond: {
751 if (!bProps_ || bProps_->sb_state_ == 0) return false;
752 double mingap = -1.0 * limits<double>::max();
753 double maxgap = 0;
754 if (vl.at(0).canConvert<double>())
755 maxgap = vl.at(0).toDouble();
756 else if (vl.at(0).canConvert<DVect2>()) {
757 DVect2 value = vl.at(0).value<DVect2>();
758 mingap = value.minComp();
759 maxgap = value.maxComp();
760 }
761 else if (!vl.at(0).isNull())
762 throw Exception("gap value %1 not recognized in method unbond in contact model %2.", vl.at(0), getName());
763 double gap = c->getGap();
764 if (gap >= mingap && gap <= maxgap) {
765 bProps_->sb_state_ = 0;
766 return true;
767 }
768 return false;
769 }
770 case kwArea: {
771 if (!userArea_) {
772 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
773#ifdef THREED
774 userArea_ = rsq * rsq * dPi;
775#else
776 userArea_ = rsq * 2.0;
777#endif
778 }
779 return true;
780 }
781
782 }
783 return false;
784 }
785
786 double ContactModelSoftBond::getEnergy(uint32 i) const {
787 // Return an energy value.
788 double ret(0.0);
789 if (!energies_)
790 return ret;
791 switch (i) {
792 case kwEStrain: return energies_->estrain_;
793 case kwESlip: return energies_->eslip_;
794 case kwEDashpot: return energies_->edashpot_;
795 }
796 assert(0);
797 return ret;
798 }
799
800 bool ContactModelSoftBond::getEnergyAccumulate(uint32 i) const {
801 // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
802 switch (i) {
803 case kwEStrain: return false;
804 case kwESlip: return true;
805 case kwEDashpot: return true;
806 }
807 assert(0);
808 return false;
809 }
810
811 void ContactModelSoftBond::setEnergy(uint32 i,const double &d) {
812 // Set an energy value.
813 if (!energies_) return;
814 switch (i) {
815 case kwEStrain: energies_->estrain_ = d; return;
816 case kwESlip: energies_->eslip_ = d; return;
817 case kwEDashpot: energies_->edashpot_= d; return;
818 }
819 assert(0);
820 return;
821 }
822
823 bool ContactModelSoftBond::validate(ContactModelMechanicalState *state,const double &) {
824 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
825 // (1) the contact is created, (2) a property of the contact that returns a true via
826 // the setProperty method has been modified and (3) when a set of cycles is executed
827 // via the {cycle N} command.
828 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
829 assert(state);
830 const IContactMechanical *c = state->getMechanicalContact();
831 assert(c);
832
833 if (state->trackEnergy_)
834 activateEnergy();
835
836 if (inheritanceField_ & KnMask)
837 updateKn(c);
838 if (inheritanceField_ & KsMask)
839 updateKs(c);
840 if (inheritanceField_ & FricMask)
841 updateFric(c);
842
843 updateStiffness(state);
844 return checkActivity(state->gap_);
845 }
846
847 static const QString knstr("kn");
848 bool ContactModelSoftBond::updateKn(const IContactMechanical *con) {
849 assert(con);
850 QVariant v1 = con->getEnd1()->getProperty(knstr);
851 QVariant v2 = con->getEnd2()->getProperty(knstr);
852 if (!v1.isValid() || !v2.isValid())
853 return false;
854 double kn1 = v1.toDouble();
855 double kn2 = v2.toDouble();
856 double val = kn_;
857 if (kn1 && kn2)
858 kn_ = kn1*kn2/(kn1+kn2);
859 else if (kn1)
860 kn_ = kn1;
861 else if (kn2)
862 kn_ = kn2;
863 return ( (kn_ != val) );
864 }
865
866 static const QString ksstr("ks");
867 bool ContactModelSoftBond::updateKs(const IContactMechanical *con) {
868 assert(con);
869 QVariant v1 = con->getEnd1()->getProperty(ksstr);
870 QVariant v2 = con->getEnd2()->getProperty(ksstr);
871 if (!v1.isValid() || !v2.isValid())
872 return false;
873 double ks1 = v1.toDouble();
874 double ks2 = v2.toDouble();
875 double val = ks_;
876 if (ks1 && ks2)
877 ks_ = ks1*ks2/(ks1+ks2);
878 else if (ks1)
879 ks_ = ks1;
880 else if (ks2)
881 ks_ = ks2;
882 return ( (ks_ != val) );
883 }
884
885 static const QString fricstr("fric");
886 bool ContactModelSoftBond::updateFric(const IContactMechanical *con) {
887 assert(con);
888 QVariant v1 = con->getEnd1()->getProperty(fricstr);
889 QVariant v2 = con->getEnd2()->getProperty(fricstr);
890 if (!v1.isValid() || !v2.isValid())
891 return false;
892 double fric1 = std::max(0.0,v1.toDouble());
893 double fric2 = std::max(0.0,v2.toDouble());
894 double val = fric_;
895 fric_ = std::min(fric1,fric2);
896 return ( (fric_ != val) );
897 }
898
899 bool ContactModelSoftBond::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
900 // The endPropertyUpdated method is called whenever a surface property (with a name
901 // that matches an inheritable contact model property name) of one of the contacting
902 // pieces is modified. This allows the contact model to update its associated
903 // properties. The return value denotes whether or not the update has affected
904 // the time step computation (by having modified the translational or rotational
905 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
906 assert(c);
907 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
908 QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
909 int idx = availableProperties.indexOf(rx)+1;
910 bool ret=false;
911
912 if (idx<=0)
913 return ret;
914
915 switch(idx) {
916 case kwKn: { //kn
917 if (inheritanceField_ & KnMask)
918 ret = updateKn(c);
919 break;
920 }
921 case kwKs: { //ks
922 if (inheritanceField_ & KsMask)
923 ret =updateKs(c);
924 break;
925 }
926 case kwFric: { //fric
927 if (inheritanceField_ & FricMask)
928 updateFric(c);
929 break;
930 }
931 }
932 return ret;
933 }
934
935 ContactModelSoftBond::StiffData ContactModelSoftBond::computeStiffData(ContactModelMechanicalState *state) const {
936 // Update contact data
937 double Cmin1 = state->end1Curvature_.x();
938 double Cmax1 = state->end1Curvature_.y();
939 double Cmax2 = state->end2Curvature_.y();
940 double dthick = (Cmin1 == 0.0) ? 1.0 : 0.0;
941 double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
942 if (userArea_)
943#ifdef THREED
944 br = std::sqrt(userArea_ / dPi);
945#else
946 br = userArea_ / 2.0;
947#endif
948 double br2 = br * br;
949 double area = dthick <= 0.0 ? dPi * br2 : 2.0*br*dthick;
950 double bi = dthick <= 0.0 ? 0.25*area*br2 : 2.0*br*br2*dthick / 3.0;
951 StiffData ret;
952 ret.reff_ = br;
953 ret.trans_ = DVect2(kn_ * area , ks_ * area);
954 ret.ang_ = DAVect(kn_ * bi);
955#if DIM==3
956 ret.ang_.rx() = ks_ * 2.0*bi;
957#endif
958 return ret;
959 }
960
961 void ContactModelSoftBond::updateStiffness(ContactModelMechanicalState *state) {
962 // first compute stiffness data
963 StiffData stiff = computeStiffData(state);
964 // Now calculate effective stiffness
965 DVect2 retT = stiff.trans_;
966 // correction if viscous damping active
967 if (dpProps_) {
968 DVect2 correct(1.0);
969 if (dpProps_->dp_nratio_)
970 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
971 if (dpProps_->dp_sratio_)
972 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
973 retT /= (correct*correct);
974 }
975 effectiveTranslationalStiffness_ = retT;
976 // Effective rotational stiffness (bending and twisting)
977 effectiveRotationalStiffness_ = stiff.ang_;
978 }
979
980 bool ContactModelSoftBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
981 assert(state);
982
983 if (state->activated()) {
984 // The contact was just activated from an inactive state
985 // Trigger the FISH callback if one is hooked up to the
986 // contact_activated event.
987 if (cmEvents_[fActivated] >= 0) {
988 auto c = state->getContact();
989 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
990 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
991 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
992 }
993 }
994 updateStiffness(state);
995 if (isBonded()) return FDLawBonded(state, timestep);
996 else return FDLawUnBonded(state, timestep);
997
998 }
999
1000 bool ContactModelSoftBond::FDLawBonded(ContactModelMechanicalState *state, const double ×tep) {
1001 // Relative translational/rotational displacement increments
1002 DVect trans = state->relativeTranslationalIncrement_;
1003 DAVect ang = state->relativeAngularIncrement_;
1004
1005 // Store previous force and moment
1006 DVect sb_F_old = sb_F_;
1007 DAVect sb_M_old = sb_M_;
1008
1009 // Update stiffness data
1010 StiffData stiff = computeStiffData(state);
1011 DVect3 geom = computeGeomData(state->getMechanicalContact());
1012 double area = geom.x();
1013 double bi = geom.y();
1014 double br = geom.z();
1015 double kn = stiff.trans_.x();
1016 double ks = stiff.trans_.y();
1017 double kb = stiff.ang_.z();
1018#if DIM==3
1019 double kt = stiff.ang_.x();
1020#else
1021 double kt = 0.0;
1022#endif
1023
1024 double nsmax0 = -(sb_F_.x() / area) + bProps_->sb_mcf_* sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z()) * br / bi;
1025
1026 // incremental normal force calculation
1027 sb_F_.rx() -= trans.x() * kn;
1028
1029 // shear force calculation
1030 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1031 // Loop over the shear components (note: the 0 component is the normal component)
1032 // and calculate the shear force.
1033 for (int i = 1; i<dim; ++i)
1034 sb_F_.rdof(i) -= trans.dof(i) * ks;
1035
1036 // moment calculation
1037 sb_M_ -= ang * stiff.ang_;
1038 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1039
1040 // maximum tensile stress at bond periphery
1041 double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_* dbend * br / bi;
1042
1043 bool softened = false;
1044 // Mode check
1045 if (state->canFail_) {
1046 if (bProps_->sb_state_ == 3 || bProps_->sb_state_ == 5) {
1047 double compVal = bProps_->sb_state_ == 3 ? bProps_->sb_ten_ : bProps_->sb_maxTen_;
1048 if (nsmax >= compVal ) {
1049 // enter softening regime
1050 // current bond elongation when softening starts
1051 // This is the elongation at the bond periphery
1052 double ls = - sb_F_.x() / kn + bProps_->sb_mcf_*dbend* br / kb;
1053 bProps_->sb_maxTen_ = compVal;
1054 if (bProps_->sb_state_ == 3)
1055 bProps_->sb_critu_ = ls /**(1.0+bProps_->sb_soft_)*/;
1056 bProps_->sb_delu_ = 0.0;
1057 bProps_->sb_delo_ = Quat::identity();
1058 if (bProps_->sb_state_ == 5 && nsmax < bProps_->sb_maxTen_)
1059 softened = true;
1060 bProps_->sb_state_ = 4;
1061 }
1062 }
1063 }
1064
1065 if (bProps_->sb_state_ == 4 && !softened && !checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1066 double ls = bProps_->sb_critu_;
1067 double lc = ls * (1.0+bProps_->sb_soft_);
1068 DVect normal(0.0);
1069 normal.rx() = 1.0;
1070 DVect backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1071 double bend = acos(qBound(-1.0,normal|backNormal,1.0));
1072 double l0 = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1073 bProps_->sb_delu_ += trans.x();
1074 bProps_->sb_delo_.increment(ang);
1075 // Take the current contact normal and rotate it in the opposite direction of
1076 // the orientation - get the angle of bend from there
1077 backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1078 bend = acos(qBound(-1.0,normal|backNormal,1.0));
1079 double l = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1080 // target tensile stress
1081 double ns = bProps_->sb_ten_*(lc-l) / (bProps_->sb_soft_*ls);
1082 if (ns > 0) {
1083 if (nsmax >= ns) {
1084 double fac = ns / nsmax;
1085 sb_F_.rx() = fac*sb_F_.x();
1086#if DIM==3
1087 sb_M_.ry() = fac*sb_M_.y();
1088#endif
1089 sb_M_.rz() = fac*sb_M_.z();
1090 } else {
1091 bProps_->sb_state_ = 5;
1092 bProps_->sb_maxTen_ = nsmax0;
1093 bProps_->sb_maxu_ = (l0-ls);
1094 }
1095 } else {
1096 sb_F_.rx() = 0.0;
1097#if DIM==3
1098 sb_M_.ry() = 0.0;
1099#endif
1100 sb_M_.rz() = 0.0;
1101 }
1102 }
1103
1104 if (state->canFail_) {
1105 /* check for normal failure */
1106 bool failed = false;
1107 if (bProps_->sb_state_ == 4) {
1108 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1109 double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_*dbend * br / bi;
1110 if (nsmax <= bProps_->sb_ten_*bProps_->sb_cut_ || checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1111 // Failed in tension
1112 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1113 bProps_->sb_state_ = 1;
1114 sb_F_.fill(0.0);
1115 sb_M_.fill(0.0);
1116 failed = true;
1117 if (cmEvents_[fBondBreak] >= 0) {
1118 auto c = state->getContact();
1119 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1120 fish::Parameter((qint64)bProps_->sb_state_),
1121 fish::Parameter(nsmax),
1122 fish::Parameter(se)
1123 };
1124 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1125 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1126 }
1127 }
1128 }
1129
1130 if (!failed) {
1131 /* check for shear failure */
1132 double dtwist = sb_M_.x();
1133 DVect bfs(sb_F_);
1134 bfs.rx() = 0.0;
1135 double dbfs = bfs.mag();
1136 double ssmax = dbfs / area + bProps_->sb_mcf_*std::abs(dtwist) * 0.5* br / bi;
1137 double ss = shearStrength(area);
1138 if (ss < 0)
1139 ss = 0;
1140 if (ss <= ssmax) {
1141 // Failed in shear
1142 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1143 bProps_->sb_state_ = 2;
1144 if (cmEvents_[fBondBreak] >= 0) {
1145 auto c = state->getContact();
1146 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1147 fish::Parameter((qint64)bProps_->sb_state_),
1148 fish::Parameter(ss),
1149 fish::Parameter(se)
1150 };
1151 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1152 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1153 }
1154 // Resolve sliding.
1155 double crit = sb_F_.x() * fric_;
1156 if (crit < 0)
1157 crit = 0;
1158 DVect sforce = sb_F_; sforce.rx() = 0.0;
1159 // The is the magnitude of the shear force.
1160 double sfmag = sforce.mag();
1161 // Sliding occurs when the magnitude of the shear force is greater than the
1162 // critical value.
1163 if (sfmag > crit) {
1164 // Lower the shear force to the critical value for sliding.
1165 double rat = crit / sfmag;
1166 sforce *= rat;
1167 sforce.rx() = sb_F_.x();
1168 sb_F_ = sforce;
1169 sb_S_ = true;
1170 }
1171
1172 // Resolve bending
1173 crit = sb_bmul_*2.1*0.25*stiff.reff_*std::abs(sb_F_.x()); // Jiang 2015
1174 DAVect test = sb_M_;
1175#if DIM==3
1176 test.rx() = 0.0;
1177#endif
1178 double tmag = test.mag();
1179 if (tmag > crit) {
1180 // Lower the bending moment to the critical value for sliding.
1181 double rat = crit / tmag;
1182 test *= rat;
1183 sb_BS_ = true;
1184 }
1185 sb_M_.rz() = test.z();
1186#if DIM==3
1187 sb_M_.ry() = test.y();
1188 // Resolve twisting
1189 crit = sb_tmul_ * 0.65*fric_* stiff.reff_*std::abs(sb_F_.x()) ; // Jiang 2015
1190 tmag = std::abs(sb_M_.x());
1191 if (tmag > crit) {
1192 // Lower the shear force to the critical value for sliding.
1193 double rat = crit / tmag;
1194 tmag = sb_M_.x() * rat;
1195 sb_TS_ = true;
1196 } else
1197 tmag = sb_M_.x();
1198 sb_M_.rx() = tmag;
1199#endif
1200 }
1201 }
1202 }
1203
1204 // Account for dashpot forces if the dashpot structure has been defined.
1205 if (dpProps_) {
1206 dpProps_->dp_F_.fill(0.0);
1207 double vcn(0.0), vcs(0.0);
1208 // Calculate the damping coefficients.
1209 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1210 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1211 // First damp the shear components
1212 for (int i = 1; i<dim; ++i)
1213 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1214 // Damp the normal component
1215 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1216 // Need to change behavior based on the dp_mode.
1217 if (bProps_->sb_state_ < 3 && (dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1218 // Limit in tension if not bonded.
1219 if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1220 dpProps_->dp_F_.rx() = -sb_F_.rx();
1221 }
1222 if (sb_S_ && dpProps_->dp_mode_ > 1) {
1223 // Limit in shear if sliding.
1224 double dfn = dpProps_->dp_F_.rx();
1225 dpProps_->dp_F_.fill(0.0);
1226 dpProps_->dp_F_.rx() = dfn;
1227 }
1228 }
1229
1230 //Compute energies if energy tracking has been enabled.
1231 if (state->trackEnergy_) {
1232 assert(energies_);
1233 energies_->estrain_ = 0.0;
1234 if (kn)
1235 // Calculate the strain energy.
1236 energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1237 if (ks) {
1238 DVect s = sb_F_;
1239 s.rx() = 0.0;
1240 double smag2 = s.mag2();
1241 // Add the shear component of the strain energy.
1242 energies_->estrain_ += 0.5*smag2 / ks;
1243
1244 if (sb_S_) {
1245 // If sliding calculate the slip energy and accumulate it.
1246 sb_F_old.rx() = 0.0;
1247 DVect avg_F_s = (s + sb_F_old)*0.5;
1248 DVect u_s_el = (s - sb_F_old) / ks;
1249 DVect u_s(0.0);
1250 for (int i = 1; i < dim; ++i)
1251 u_s.rdof(i) = trans.dof(i);
1252 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1253 }
1254 }
1255 // Add the bending/twisting resistance energy contributions.
1256 if (kb) {
1257 DAVect tmp = sb_M_;
1258#ifdef THREED
1259 tmp.rx() = 0.0;
1260#endif
1261 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1262 if (sb_BS_) {
1263 // accumulate bending slip energy.
1264 DAVect tmp_old = sb_M_old;
1265#ifdef THREED
1266 tmp_old.rx() = 0.0;
1267#endif
1268 DAVect avg_M = (tmp + tmp_old)*0.5;
1269 DAVect t_s_el = (tmp - tmp_old) / kb;
1270 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1271 }
1272 }
1273#ifdef THREED
1274 if (kt) {
1275 double mt = std::abs(sb_M_.x());
1276 energies_->estrain_ += 0.5*mt*mt / kt;
1277 if (sb_TS_) {
1278 // accumulate twisting slip energy.
1279 DAVect tmp(0.0);
1280 DAVect tmp_old(0.0);
1281 tmp.rx() = sb_M_.x();
1282 tmp_old.rx() = sb_M_old.x();
1283 DAVect avg_M = (tmp + tmp_old)*0.5;
1284 DAVect t_s_el = (tmp - tmp_old) / kt;
1285 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1286 }
1287 }
1288#endif
1289
1290 if (dpProps_) {
1291 // Calculate damping energy (accumulated) if the dashpots are active.
1292 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1293 }
1294 }
1295
1296 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1297 assert(sb_F_ == sb_F_);
1298 assert(sb_M_ == sb_M_);
1299 return true;
1300 }
1301
1302 bool ContactModelSoftBond::FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep) {
1303
1304 // Relative translational/rotational displacement increments
1305 DVect trans = state->relativeTranslationalIncrement_;
1306 DAVect ang = state->relativeAngularIncrement_;
1307 double overlap = rgap_ - state->gap_;
1308 double correction = 1.0;
1309 if (state->activated() && sb_mode_ == 0 && trans.x()) {
1310 correction = -1.0*overlap / trans.x();
1311 if (correction < 0)
1312 correction = 1.0;
1313 }
1314
1315 // Store previous force and moment
1316 DVect sb_F_old = sb_F_;
1317 DAVect sb_M_old = sb_M_;
1318
1319 // Update stiffness data
1320 StiffData stiff = computeStiffData(state);
1321 double kn = stiff.trans_.x();
1322 double ks = stiff.trans_.y();
1323 double kb = stiff.ang_.z();
1324#if DIM==3
1325 double kt = stiff.ang_.x();
1326#endif
1327 // absolute/incremental normal force calculation
1328 if (sb_mode_==0)
1329 sb_F_.rx() = overlap * kn;
1330 else
1331 sb_F_.rx() -= trans.x() * kn;
1332 // Normal force can only be positive if unbonded
1333 sb_F_.rx() = std::max(0.0, sb_F_.x());
1334
1335 // Calculate the trial shear force.
1336 DVect sforce(0.0);
1337 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1338 // Loop over the shear components (note: the 0 component is the normal component)
1339 // and calculate the shear force.
1340 for (int i = 1; i<dim; ++i)
1341 sforce.rdof(i) = sb_F_.dof(i) - trans.dof(i) * ks;
1342
1343 // Calculate the trial moment.
1344 DAVect mom = sb_M_ - ang*stiff.ang_;
1345
1346 // If the SOLVE ELASTIC command is given then the
1347 // canFail state is set to FALSE. Otherwise it is always TRUE.
1348 if (state->canFail_) {
1349 bool changed = false;
1350 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1351 bool slip_changed = false;
1352 double crit = sb_F_.x() * fric_;
1353 // The is the magnitude of the shear force.
1354 double sfmag = sforce.mag();
1355 // Sliding occurs when the magnitude of the shear force is greater than the
1356 // critical value.
1357 if (sfmag > crit) {
1358 // Lower the shear force to the critical value for sliding.
1359 double rat = crit / sfmag;
1360 sforce *= rat;
1361 if (!sb_S_) {
1362 slip_changed = true;
1363 changed = true;
1364 }
1365 sb_S_ = true;
1366 }
1367 else {
1368 if (sb_S_) {
1369 slip_changed = true;
1370 changed = true;
1371 }
1372 sb_S_ = false;
1373 }
1374
1375 // Resolve bending
1376 bool bslip_changed = false;
1377 crit = sb_bmul_ * 2.1*0.25*sb_F_.x() * stiff.reff_; // Jiang 2015
1378 DAVect test = mom;
1379#if DIM==3
1380 test.rx() = 0.0;
1381#endif
1382 double tmag = test.mag();
1383 if (tmag > crit) {
1384 // Lower the bending moment to the critical value for sliding.
1385 double rat = crit / tmag;
1386 test *= rat;
1387 if (!sb_BS_) {
1388 bslip_changed = true;
1389 changed = true;
1390 }
1391 sb_BS_ = true;
1392 }
1393 else {
1394 if (sb_BS_) {
1395 bslip_changed = true;
1396 changed = true;
1397 }
1398 sb_BS_ = false;
1399 }
1400 mom.rz() = test.z();
1401#if DIM==3
1402 mom.ry() = test.y();
1403 // Resolve twisting
1404 bool tslip_changed = false;
1405 crit = sb_tmul_ * 0.65*fric_*sb_F_.x() * stiff.reff_; // Jiang 2015
1406 tmag = std::abs(mom.x());
1407 if (tmag > crit) {
1408 // Lower the twisting moment to the critical value for sliding.
1409 double rat = crit / tmag;
1410 mom.rx() *= rat;
1411 if (!sb_TS_) {
1412 tslip_changed = true;
1413 changed = true;
1414 }
1415 sb_TS_ = true;
1416 } else {
1417 if (sb_TS_) {
1418 tslip_changed = true;
1419 changed = true;
1420 }
1421 sb_TS_ = false;
1422 }
1423#endif
1424 if (changed && cmEvents_[fSlipChange] >= 0) {
1425 qint64 code = 0;
1426 if (slip_changed) {
1427 code = 1;
1428 if (bslip_changed) {
1429 code = 4;
1430#if DIM==3
1431 if (tslip_changed)
1432 code = 7;
1433#endif
1434 }
1435 }
1436 else if (bslip_changed) {
1437 code = 2;
1438#if DIM==3
1439 if (tslip_changed)
1440 code = 6;
1441#endif
1442 }
1443#if DIM==3
1444 else if (tslip_changed) {
1445 code = 3;
1446 if (slip_changed)
1447 code = 5;
1448 }
1449#endif
1450 auto c = state->getContact();
1451 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1452 fish::Parameter(code),
1453 fish::Parameter(sb_S_),
1454 fish::Parameter(sb_BS_)
1455#ifdef THREED
1456 ,fish::Parameter(sb_TS_)
1457#endif
1458 };
1459 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1460 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1461 }
1462 }
1463
1464 // Set the shear components of the total force.
1465 for (int i = 1; i<dim; ++i)
1466 sb_F_.rdof(i) = sforce.dof(i);
1467
1468 // Set the moment.
1469 sb_M_ = mom;
1470
1471 // Account for dashpot forces if the dashpot structure has been defined.
1472 if (dpProps_) {
1473 dpProps_->dp_F_.fill(0.0);
1474 double vcn(0.0), vcs(0.0);
1475 // Calculate the damping coefficients.
1476 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1477 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1478 // First damp the shear components
1479 for (int i = 1; i<dim; ++i)
1480 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1481 // Damp the normal component
1482 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1483 // Need to change behavior based on the dp_mode.
1484 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1485 // Limit in tension if not bonded.
1486 if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1487 dpProps_->dp_F_.rx() = -sb_F_.rx();
1488 }
1489 if (sb_S_ && dpProps_->dp_mode_ > 1) {
1490 // Limit in shear if not sliding.
1491 double dfn = dpProps_->dp_F_.rx();
1492 dpProps_->dp_F_.fill(0.0);
1493 dpProps_->dp_F_.rx() = dfn;
1494 }
1495 }
1496
1497 //Compute energies if energy tracking has been enabled.
1498 if (state->trackEnergy_) {
1499 assert(energies_);
1500 energies_->estrain_ = 0.0;
1501 if (kn_)
1502 // Calculate the strain energy.
1503 energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1504 if (ks_) {
1505 DVect s = sb_F_;
1506 s.rx() = 0.0;
1507 double smag2 = s.mag2();
1508 // Add the shear component of the strain energy.
1509 energies_->estrain_ += 0.5*smag2 / ks;
1510
1511 if (sb_S_) {
1512 // If sliding calculate the slip energy and accumulate it.
1513 sb_F_old.rx() = 0.0;
1514 DVect avg_F_s = (s + sb_F_old)*0.5;
1515 DVect u_s_el = (s - sb_F_old) / ks;
1516 DVect u_s(0.0);
1517 for (int i = 1; i < dim; ++i)
1518 u_s.rdof(i) = trans.dof(i);
1519 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1520 }
1521 }
1522 // Add the bending/twisting resistance energy contributions.
1523 if (kb) {
1524 DAVect tmp = sb_M_;
1525#ifdef THREED
1526 tmp.rx() = 0.0;
1527#endif
1528 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1529 if (sb_BS_) {
1530 // accumulate bending slip energy.
1531 DAVect tmp_old = sb_M_old;
1532#ifdef THREED
1533 tmp_old.rx() = 0.0;
1534#endif
1535 DAVect avg_M = (tmp + tmp_old)*0.5;
1536 DAVect t_s_el = (tmp - tmp_old) / kb;
1537 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1538 }
1539 }
1540#ifdef THREED
1541 if (kt) {
1542 double mt = std::abs(sb_M_.x());
1543 energies_->estrain_ += 0.5*mt*mt / kt;
1544 if (sb_TS_) {
1545 // accumulate twisting slip energy.
1546 DAVect tmp(0.0);
1547 DAVect tmp_old(0.0);
1548 tmp.rx() = sb_M_.x();
1549 tmp_old.rx() = sb_M_old.x();
1550 DAVect avg_M = (tmp + tmp_old)*0.5;
1551 DAVect t_s_el = (tmp - tmp_old) / kt;
1552 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1553 }
1554 }
1555#endif
1556
1557 if (dpProps_) {
1558 // Calculate damping energy (accumulated) if the dashpots are active.
1559 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1560 }
1561 }
1562
1563 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1564 assert(sb_F_ == sb_F_);
1565 assert(sb_M_ == sb_M_);
1566 return true;
1567 }
1568
1569 bool ContactModelSoftBond::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
1570 // Account for thermal expansion in incremental mode
1571 if (sb_mode_ == 0 || ts->gapInc_ == 0.0) return false;
1572 DVect finc(0.0);
1573 finc.rx() = kn_ * ts->gapInc_;
1574 sb_F_ -= finc;
1575 return true;
1576 }
1577
1578 void ContactModelSoftBond::setForce(const DVect &v,IContact *c) {
1579 sb_F(v);
1580 if (v.x() > 0)
1581 rgap_ = c->getGap() + v.x() / (kn_ * computeGeomData(convert_getcast<IContactMechanical>(c)).x());
1582 }
1583
1584 void ContactModelSoftBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
1585 // Only called for contacts with wall facets when the wall resolution scheme
1586 // is set to full!
1587 // Only do something if the contact model is of the same type
1588 if (old->getContactModel()->getName().compare("softbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1589 ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1590#ifdef THREED
1591 // Need to rotate just the shear component from oldSystem to newSystem
1592
1593 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
1594 DVect axis = oldSystem.e1() & newSystem.e1();
1595 double c, ang, s;
1596 DVect re2;
1597 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1598 axis = axis.unit();
1599 c = oldSystem.e1()|newSystem.e1();
1600 if (c > 0)
1601 c = std::min(c,1.0);
1602 else
1603 c = std::max(c,-1.0);
1604 ang = acos(c);
1605 s = sin(ang);
1606 double t = 1. - c;
1607 DMatrix<3,3> rm;
1608 rm.get(0,0) = t*axis.x()*axis.x() + c;
1609 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
1610 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
1611 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
1612 rm.get(1,1) = t*axis.y()*axis.y() + c;
1613 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
1614 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
1615 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
1616 rm.get(2,2) = t*axis.z()*axis.z() + c;
1617 re2 = rm*oldSystem.e2();
1618 }
1619 else
1620 re2 = oldSystem.e2();
1621 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1622 axis = re2 & newSystem.e2();
1623 DVect2 tpf;
1624 DVect2 tpm;
1625 DMatrix<2,2> m;
1626 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1627 axis = axis.unit();
1628 c = re2|newSystem.e2();
1629 if (c > 0)
1630 c = std::min(c,1.0);
1631 else
1632 c = std::max(c,-1.0);
1633 ang = acos(c);
1634 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1635 ang *= -1;
1636 s = sin(ang);
1637 m.get(0,0) = c;
1638 m.get(1,0) = s;
1639 m.get(0,1) = -m.get(1,0);
1640 m.get(1,1) = m.get(0,0);
1641 tpf = m*DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1642 tpm = m*DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1643 } else {
1644 m.get(0,0) = 1.;
1645 m.get(0,1) = 0.;
1646 m.get(1,0) = 0.;
1647 m.get(1,1) = 1.;
1648 tpf = DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1649 tpm = DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1650 }
1651 DVect pforce = DVect(0,tpf.x(),tpf.y());
1652 //DVect pm = DVect(0,tpm.x(),tpm.y());
1653#else
1654 oldSystem;
1655 newSystem;
1656 DVect pforce = DVect(0,oldCm->sb_F_.y());
1657 //DVect pm = DVect(0,oldCm->sb_M_.y());
1658#endif
1659 for (int i=1; i<dim; ++i)
1660 sb_F_.rdof(i) += pforce.dof(i);
1661 if (sb_mode_ && oldCm->sb_mode_)
1662 sb_F_.rx() = oldCm->sb_F_.x();
1663 oldCm->sb_F_ = DVect(0.0);
1664 oldCm->sb_M_ = DAVect(0.0);
1665 if (dpProps_ && oldCm->dpProps_) {
1666#ifdef THREED
1667 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1668 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1669#else
1670 pforce = oldCm->dpProps_->dp_F_;
1671#endif
1672 dpProps_->dp_F_ += pforce;
1673 oldCm->dpProps_->dp_F_ = DVect(0.0);
1674 }
1675 if(oldCm->getEnergyActivated()) {
1676 activateEnergy();
1677 energies_->estrain_ = oldCm->energies_->estrain_;
1678 energies_->edashpot_ = oldCm->energies_->edashpot_;
1679 energies_->eslip_ = oldCm->energies_->eslip_;
1680 oldCm->energies_->estrain_ = 0.0;
1681 oldCm->energies_->edashpot_ = 0.0;
1682 oldCm->energies_->eslip_ = 0.0;
1683 }
1684 rgap_ = oldCm->rgap_;
1685 }
1686 assert(sb_F_ == sb_F_);
1687 }
1688
1689 void ContactModelSoftBond::setNonForcePropsFrom(IContactModel *old) {
1690 // Only called for contacts with wall facets when the wall resolution scheme
1691 // is set to full!
1692 // Only do something if the contact model is of the same type
1693 if (old->getName().compare("softbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1694 ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1695 kn_ = oldCm->kn_;
1696 ks_ = oldCm->ks_;
1697 fric_ = oldCm->fric_;
1698 sb_bmul_ = oldCm->sb_bmul_;
1699 sb_tmul_ = oldCm->sb_tmul_;
1700 sb_mode_ = oldCm->sb_mode_;
1701 sb_rmul_ = oldCm->sb_rmul_;
1702 sb_S_ = oldCm->sb_S_;
1703 sb_BS_ = oldCm->sb_BS_;
1704 sb_TS_ = oldCm->sb_TS_;
1705 rgap_ = oldCm->rgap_;
1706 userArea_ = oldCm->userArea_;
1707
1708 if (oldCm->dpProps_) {
1709 if (!dpProps_)
1710 dpProps_ = NEW dpProps();
1711 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1712 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1713 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1714 }
1715 if (oldCm->bProps_) {
1716 if (!bProps_)
1717 bProps_ = NEW bProps();
1718 bProps_->sb_mcf_ = oldCm->bProps_->sb_mcf_;
1719 bProps_->sb_fa_ = oldCm->bProps_->sb_fa_;
1720 bProps_->sb_state_ = oldCm->bProps_->sb_state_;
1721 bProps_->sb_coh_ = oldCm->bProps_->sb_coh_;
1722 bProps_->sb_ten_ = oldCm->bProps_->sb_ten_;
1723 bProps_->sb_maxTen_ = oldCm->bProps_->sb_maxTen_;
1724 bProps_->sb_cut_ = oldCm->bProps_->sb_cut_;
1725 bProps_->sb_delu_ = oldCm->bProps_->sb_delu_;
1726 bProps_->sb_delo_ = oldCm->bProps_->sb_delo_;
1727 bProps_->sb_maxu_ = oldCm->bProps_->sb_maxu_;
1728 bProps_->sb_critu_ = oldCm->bProps_->sb_critu_;
1729 }
1730
1731 }
1732 }
1733
1734 DVect ContactModelSoftBond::getForce(const IContactMechanical *) const {
1735 DVect ret(sb_F_);
1736 if (dpProps_)
1737 ret += dpProps_->dp_F_;
1738 return ret;
1739 }
1740
1741 DAVect ContactModelSoftBond::getMomentOn1(const IContactMechanical *c) const {
1742 DVect force = getForce(c);
1743 DAVect ret(sb_M_);
1744 c->updateResultingTorqueOn1Local(force,&ret);
1745 return ret;
1746 }
1747
1748 DAVect ContactModelSoftBond::getMomentOn2(const IContactMechanical *c) const {
1749 DVect force = getForce(c);
1750 DAVect ret(sb_M_);
1751 c->updateResultingTorqueOn2Local(force,&ret);
1752 return ret;
1753 }
1754
1755 DVect3 ContactModelSoftBond::computeGeomData(const IContactMechanical *c) const {
1756 double Cmax1 = c->getEnd1Curvature().y();
1757 double Cmax2 = c->getEnd2Curvature().y();
1758 double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1759 if (userArea_)
1760#ifdef THREED
1761 br = std::sqrt(userArea_ / dPi);
1762#else
1763 br = userArea_ / 2.0;
1764#endif
1765 double br2 = br * br;
1766#ifdef TWOD
1767 double area = 2.0*br;
1768 double bi = 2.0*br*br2 / 3.0;
1769#else
1770 double area = dPi * br2;
1771 double bi = 0.25*area*br2;
1772#endif
1773 return DVect3(area, bi, br);
1774 }
1775
1776 DVect2 ContactModelSoftBond::SMax(const IContactMechanical *c) const {
1777 DVect3 data = computeGeomData(c);
1778 double area = data.x();
1779 double bi = data.y();
1780 double br = data.z();
1781 /* maximum stresses */
1782 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1783 double dtwist = sb_M_.x();
1784 DVect bfs(sb_F_);
1785 bfs.rx() = 0.0;
1786 double dbfs = bfs.mag();
1787 double nsmax = -(sb_F_.x() / area) + dbend * br / bi;
1788 double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
1789 return DVect2(nsmax, ssmax);
1790 }
1791
1792 double ContactModelSoftBond::shearStrength(const double &area) const {
1793 if (!bProps_) return 0.0;
1794 double sig = -1.0*sb_F_.x() / area;
1795 double nstr = bProps_->sb_state_ > 2 ? bProps_->sb_ten_ : 0.0;
1796 return sig <= nstr ? bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*sig
1797 : bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*nstr;
1798 }
1799
1800
1801 double ContactModelSoftBond::strainEnergy(double kn,double ks,double kb,double kt) const {
1802 double ret(0.0);
1803 if (kn)
1804 ret = 0.5 * sb_F_.x() * sb_F_.x() / kn;
1805 if (ks) {
1806 DVect tmp = sb_F_;
1807 tmp.rx() = 0.0;
1808 double smag2 = tmp.mag2();
1809 ret += 0.5 * smag2 / ks;
1810 }
1811
1812 if (kt)
1813 ret += 0.5 * sb_M_.x() * sb_M_.x() / kt;
1814 if (kb) {
1815 DAVect tmp = sb_M_;
1816#ifdef THREED
1817 tmp.rx() = 0.0;
1818 double smag2 = tmp.mag2();
1819#else
1820 double smag2 = tmp.z() * tmp.z();
1821#endif
1822 ret += 0.5 * smag2 / kb;
1823 }
1824 return ret;
1825 }
1826
1827} // namespace cmodelsxd
1828// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Sep 26, 2024 |
